ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Molecule.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Molecule.cpp (file contents):
Revision 438 by chuckv, Mon Mar 31 21:50:59 2003 UTC vs.
Revision 1136 by tim, Tue Apr 27 16:26:44 2004 UTC

# Line 1 | Line 1
1 < #include <cstdlib>
1 > #include <stdlib.h>
2  
3  
4   #include "Molecule.hpp"
# Line 15 | Line 15 | Molecule::Molecule( void ){
15  
16   }
17  
18 –
19 –
18   Molecule::~Molecule( void ){
19    int i;
20    
# Line 40 | Line 38 | Molecule::~Molecule( void ){
38      delete[] myTorsions;
39    }
40  
43 –  if( myExcludes != NULL ){
44 –    for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
45 –    delete[] myExcludes;
46 –  }
41   }
42  
43  
44   void Molecule::initialize( molInit &theInit ){
45 <
45 >  double totMass;
46 >  
47    nAtoms = theInit.nAtoms;
48    nMembers = nAtoms;
49    nBonds = theInit.nBonds;
50    nBends = theInit.nBends;
51    nTorsions = theInit.nTorsions;
52 <  nExcludes = theInit.nExcludes;
52 >  nRigidBodies = theInit.nRigidBodies;
53    nOriented = theInit.nOriented;
54  
55    myAtoms = theInit.myAtoms;
56    myBonds = theInit.myBonds;
57    myBends = theInit.myBends;
58    myTorsions = theInit.myTorsions;
59 <  myExcludes = theInit.myExcludes;
60 <    
59 >  myRigidBodies = theInit.myRigidBodies;
60 >
61 >  myIntegrableObjects = theInit.myIntegrableObjects;
62 >
63 >  for (int i = 0; i < myRigidBodies.size(); i++)
64 >      myRigidBodies[i]->calcRefCoords();
65 >
66 >
67 >  //the mass ratio will never change during the simulation. Thus, we could
68 >  //just calculate it at the begining of the simulation
69 >  totMass = getTotalMass();
70 >  for(int i = 0; i < nAtoms; i ++)
71 >    myAtoms[i]->setMassRatio(myAtoms[i]->getMass()/totMass);  
72   }
73  
74   void Molecule::calcForces( void ){
75    
76    int i;
77 +  double com[3];
78  
79 +  for(i=0; i<myRigidBodies.size(); i++) {
80 +    myRigidBodies[i]->updateAtoms();
81 +  }
82 +
83 +  //calculate the center of mass of the molecule
84 +  getCOM(com);  
85 +  for(int i = 0; i < nAtoms; i ++)
86 +    myAtoms[i]->setRc(com);  
87 +  
88 +
89    for(i=0; i<nBonds; i++){
90      myBonds[i]->calc_forces();
91    }
# Line 80 | Line 97 | void Molecule::calcForces( void ){
97    for(i=0; i<nTorsions; i++){
98      myTorsions[i]->calc_forces();
99    }
100 +
101 +  // Rigid Body forces and torques are done after the fortran force loop
102 +
103   }
104  
105  
# Line 87 | Line 107 | double Molecule::getPotential( void ){
107    
108    int i;
109    double myPot = 0.0;
110 +
111 +  for(i=0; i<myRigidBodies.size(); i++) {
112 +    myRigidBodies[i]->updateAtoms();
113 +  }
114    
115    for(i=0; i<nBonds; i++){
116      myPot += myBonds[i]->get_potential();
# Line 118 | Line 142 | void Molecule::printMe( void ){
142    for(i=0; i<nTorsions; i++){
143      myTorsions[i]->printMe();
144    }
145 +
146   }
147 +
148 + void Molecule::moveCOM(double delta[3]){
149 +  double aPos[3];
150 +  int i, j;
151 +
152 +  for(i=0; i<myIntegrableObjects.size(); i++) {
153 +    if(myIntegrableObjects[i] != NULL ) {
154 +      
155 +      myIntegrableObjects[i]->getPos( aPos );
156 +      
157 +      for (j=0; j< 3; j++)
158 +        aPos[j] += delta[j];
159 +
160 +      myIntegrableObjects[i]->setPos( aPos );
161 +    }
162 +  }
163 +
164 +  for(i=0; i<myRigidBodies.size(); i++) {
165 +
166 +      myRigidBodies[i]->getPos( aPos );
167 +
168 +      for (j=0; j< 3; j++)
169 +        aPos[j] += delta[j];
170 +      
171 +      myRigidBodies[i]->setPos( aPos );
172 +    }
173 + }
174 +
175 + void Molecule::atoms2rigidBodies( void ) {
176 +  int i;
177 +  for (i = 0; i < myRigidBodies.size(); i++) {
178 +    myRigidBodies[i]->calcForcesAndTorques();  
179 +  }
180 + }
181 +
182 + void Molecule::getCOM( double COM[3] ) {
183 +
184 +  double mass, mtot;
185 +  double aPos[3];
186 +  int i, j;
187 +
188 +  for (j=0; j<3; j++)
189 +    COM[j] = 0.0;
190 +
191 +  mtot   = 0.0;
192 +
193 +  for (i=0; i < myIntegrableObjects.size(); i++) {
194 +    if (myIntegrableObjects[i] != NULL) {
195 +
196 +      mass = myIntegrableObjects[i]->getMass();
197 +      mtot   += mass;
198 +      
199 +      myIntegrableObjects[i]->getPos( aPos );
200 +
201 +      for( j = 0; j < 3; j++)
202 +        COM[j] += aPos[j] * mass;
203 +
204 +    }
205 +  }
206 +
207 +  for (j = 0; j < 3; j++)
208 +    COM[j] /= mtot;
209 + }
210 +
211 + double Molecule::getCOMvel( double COMvel[3] ) {
212 +
213 +  double mass, mtot;
214 +  double aVel[3];
215 +  int i, j;
216 +
217 +
218 +  for (j=0; j<3; j++)
219 +    COMvel[j] = 0.0;
220 +
221 +  mtot   = 0.0;
222 +
223 +  for (i=0; i < myIntegrableObjects.size(); i++) {
224 +    if (myIntegrableObjects[i] != NULL) {
225 +
226 +      mass = myIntegrableObjects[i]->getMass();
227 +      mtot   += mass;
228 +
229 +      myIntegrableObjects[i]->getVel(aVel);
230 +
231 +      for (j=0; j<3; j++)
232 +        COMvel[j] += aVel[j]*mass;
233 +
234 +    }
235 +  }
236 +
237 +  for (j=0; j<3; j++)
238 +    COMvel[j] /= mtot;
239 +
240 +  return mtot;
241 +
242 + }
243 +
244 + double Molecule::getTotalMass()
245 + {
246 +
247 +  double totalMass;
248 +  
249 +  totalMass = 0;
250 +  for(int i =0; i < myIntegrableObjects.size(); i++){
251 +    totalMass += myIntegrableObjects[i]->getMass();
252 +  }
253 +
254 +  return totalMass;
255 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines