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 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC vs.
Revision 1140 by tim, Wed Apr 28 22:34:02 2004 UTC

# Line 12 | Line 12 | Molecule::Molecule( void ){
12    myBonds = NULL;
13    myBends = NULL;
14    myTorsions = NULL;
15 <  myRigidBodies = NULL;
16 <
15 >  hasMassRatio = false;
16   }
17  
19 –
20 –
18   Molecule::~Molecule( void ){
19    int i;
20    
# Line 41 | Line 38 | Molecule::~Molecule( void ){
38      delete[] myTorsions;
39    }
40  
44 –  if( myRigidBodies != NULL ){
45 –    for(i=0; i<nRigidBodies; i++) if(myRigidBodies[i] != NULL )
46 –      delete myRigidBodies[i];
47 –    delete[] myRigidBodies;
48 –  }
49 –  
41   }
42  
43  
# Line 65 | Line 56 | void Molecule::initialize( molInit &theInit ){
56    myBends = theInit.myBends;
57    myTorsions = theInit.myTorsions;
58    myRigidBodies = theInit.myRigidBodies;
59 <    
59 >
60 >  myIntegrableObjects = theInit.myIntegrableObjects;
61 >
62 >  for (int i = 0; i < myRigidBodies.size(); i++)
63 >      myRigidBodies[i]->calcRefCoords();
64 >
65   }
66  
67   void Molecule::calcForces( void ){
68    
69    int i;
70 +  double com[3];
71  
72 <  for(i=0; i<nRigidBodies; i++) {
72 >  for(i=0; i<myRigidBodies.size(); i++) {
73      myRigidBodies[i]->updateAtoms();
74    }
75  
76 +  //the mass ratio will never change during the simulation. Thus, we could
77 +  //just calculate it at the begining of the simulation
78 +  if (!hasMassRatio){
79 +    double totMass = getTotalMass();
80 +    for(int i = 0; i < nAtoms; i ++)
81 +      myAtoms[i]->setMassRatio(myAtoms[i]->getMass()/totMass);  
82 +    hasMassRatio =  true;
83 +  }
84 +
85 +  //calculate the center of mass of the molecule
86 +  getCOM(com);  
87 +  for(int i = 0; i < nAtoms; i ++)
88 +    myAtoms[i]->setRc(com);  
89 +  
90 +
91    for(i=0; i<nBonds; i++){
92      myBonds[i]->calc_forces();
93    }
# Line 97 | Line 109 | double Molecule::getPotential( void ){
109    
110    int i;
111    double myPot = 0.0;
112 +
113 +  for(i=0; i<myRigidBodies.size(); i++) {
114 +    myRigidBodies[i]->updateAtoms();
115 +  }
116    
117    for(i=0; i<nBonds; i++){
118      myPot += myBonds[i]->get_potential();
# Line 135 | Line 151 | void Molecule::moveCOM(double delta[3]){
151    double aPos[3];
152    int i, j;
153  
154 <  for(i=0; i<nAtoms; i++) {
155 <    if(myAtoms[i] != NULL ) {
154 >  for(i=0; i<myIntegrableObjects.size(); i++) {
155 >    if(myIntegrableObjects[i] != NULL ) {
156        
157 <      myAtoms[i]->getPos( aPos );
157 >      myIntegrableObjects[i]->getPos( aPos );
158        
159        for (j=0; j< 3; j++)
160          aPos[j] += delta[j];
161  
162 <      myAtoms[i]->setPos( aPos );
162 >      myIntegrableObjects[i]->setPos( aPos );
163      }
164    }
165  
166 <  for(i=0; i<nRigidBodies; i++) {
166 >  for(i=0; i<myRigidBodies.size(); i++) {
167  
152 –    if (myRigidBodies[i] != NULL) {
153 –      
168        myRigidBodies[i]->getPos( aPos );
169  
170        for (j=0; j< 3; j++)
# Line 158 | Line 172 | void Molecule::moveCOM(double delta[3]){
172        
173        myRigidBodies[i]->setPos( aPos );
174      }
161 –  }  
175   }
176  
177   void Molecule::atoms2rigidBodies( void ) {
178    int i;
179 <  for (i = 0; i < nRigidBodies; i++) {
180 <    if (myRigidBodies[i] != NULL) {
168 <      myRigidBodies[i]->calcForcesAndTorques();  
169 <    }
179 >  for (i = 0; i < myRigidBodies.size(); i++) {
180 >    myRigidBodies[i]->calcForcesAndTorques();  
181    }
182   }
183  
# Line 181 | Line 192 | void Molecule::getCOM( double COM[3] ) {
192  
193    mtot   = 0.0;
194  
195 <  for (i=0; i < nAtoms; i++) {
196 <    if (myAtoms[i] != NULL) {
195 >  for (i=0; i < myIntegrableObjects.size(); i++) {
196 >    if (myIntegrableObjects[i] != NULL) {
197  
198 <      mass = myAtoms[i]->getMass();
198 >      mass = myIntegrableObjects[i]->getMass();
199        mtot   += mass;
200        
201 <      myAtoms[i]->getPos( aPos );
201 >      myIntegrableObjects[i]->getPos( aPos );
202  
203        for( j = 0; j < 3; j++)
204          COM[j] += aPos[j] * mass;
# Line 211 | Line 222 | double Molecule::getCOMvel( double COMvel[3] ) {
222  
223    mtot   = 0.0;
224  
225 <  for (i=0; i < nAtoms; i++) {
226 <    if (myAtoms[i] != NULL) {
225 >  for (i=0; i < myIntegrableObjects.size(); i++) {
226 >    if (myIntegrableObjects[i] != NULL) {
227  
228 <      mass = myAtoms[i]->getMass();
228 >      mass = myIntegrableObjects[i]->getMass();
229        mtot   += mass;
230  
231 <      myAtoms[i]->getVel(aVel);
231 >      myIntegrableObjects[i]->getVel(aVel);
232  
233        for (j=0; j<3; j++)
234          COMvel[j] += aVel[j]*mass;
# Line 234 | Line 245 | double Molecule::getTotalMass()
245  
246   double Molecule::getTotalMass()
247   {
248 <  int natoms;
238 <  Atom** atoms;
248 >
249    double totalMass;
250    
241 –  natoms = getNAtoms();
242 –  atoms = getMyAtoms();
251    totalMass = 0;
252 <  for(int i =0; i < natoms; i++){
253 <    totalMass += atoms[i]->getMass();
252 >  for(int i =0; i < myIntegrableObjects.size(); i++){
253 >    totalMass += myIntegrableObjects[i]->getMass();
254    }
255  
256    return totalMass;

Diff Legend

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