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 489 by mmeineke, Thu Apr 10 20:08:56 2003 UTC vs.
Revision 1097 by gezelter, Mon Apr 12 20:32:20 2004 UTC

# Line 1 | Line 1
1 < #include <cstdlib>
1 > #include <stdlib.h>
2  
3  
4   #include "Molecule.hpp"
# Line 12 | Line 12 | Molecule::Molecule( void ){
12    myBonds = NULL;
13    myBends = NULL;
14    myTorsions = NULL;
15 +  myRigidBodies = NULL;
16  
17   }
18  
# Line 40 | Line 41 | Molecule::~Molecule( void ){
41      delete[] myTorsions;
42    }
43  
44 <  if( myExcludes != NULL ){
45 <    for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
46 <    delete[] myExcludes;
44 >  if( myRigidBodies != NULL ){
45 >    for(i=0; i<nRigidBodies; i++) if(myRigidBodies[i] != NULL )
46 >      delete myRigidBodies[i];
47 >    delete[] myRigidBodies;
48    }
49 +  
50   }
51  
52  
# Line 54 | Line 57 | void Molecule::initialize( molInit &theInit ){
57    nBonds = theInit.nBonds;
58    nBends = theInit.nBends;
59    nTorsions = theInit.nTorsions;
60 <  nExcludes = theInit.nExcludes;
60 >  nRigidBodies = theInit.nRigidBodies;
61    nOriented = theInit.nOriented;
62  
63    myAtoms = theInit.myAtoms;
64    myBonds = theInit.myBonds;
65    myBends = theInit.myBends;
66    myTorsions = theInit.myTorsions;
67 <  myExcludes = theInit.myExcludes;
67 >  myRigidBodies = theInit.myRigidBodies;
68      
69   }
70  
# Line 69 | Line 72 | void Molecule::calcForces( void ){
72    
73    int i;
74  
75 +  for(i=0; i<nRigidBodies; i++) {
76 +    myRigidBodies[i]->updateAtoms();
77 +  }
78 +
79    for(i=0; i<nBonds; i++){
80      myBonds[i]->calc_forces();
81    }
# Line 80 | Line 87 | void Molecule::calcForces( void ){
87    for(i=0; i<nTorsions; i++){
88      myTorsions[i]->calc_forces();
89    }
90 +
91 +  // Rigid Body forces and torques are done after the fortran force loop
92 +
93   }
94  
95  
# Line 118 | Line 128 | void Molecule::printMe( void ){
128    for(i=0; i<nTorsions; i++){
129      myTorsions[i]->printMe();
130    }
131 +
132   }
133  
134   void Molecule::moveCOM(double delta[3]){
135 <  double x, y, z;
136 <  int i;
135 >  double aPos[3];
136 >  int i, j;
137  
138    for(i=0; i<nAtoms; i++) {
139      if(myAtoms[i] != NULL ) {
140 +      
141 +      myAtoms[i]->getPos( aPos );
142 +      
143 +      for (j=0; j< 3; j++)
144 +        aPos[j] += delta[j];
145  
146 <      x = myAtoms[i]->getX() + delta[0];
147 <      y = myAtoms[i]->getY() + delta[1];
148 <      z = myAtoms[i]->getZ() + delta[2];
146 >      myAtoms[i]->setPos( aPos );
147 >    }
148 >  }
149  
150 <      myAtoms[i]->setX(x);
151 <      myAtoms[i]->setY(y);
152 <      myAtoms[i]->setZ(z);
150 >  for(i=0; i<nRigidBodies; i++) {
151 >
152 >    if (myRigidBodies[i] != NULL) {
153 >      
154 >      myRigidBodies[i]->getPos( aPos );
155 >
156 >      for (j=0; j< 3; j++)
157 >        aPos[j] += delta[j];
158 >      
159 >      myRigidBodies[i]->setPos( aPos );
160      }
161 +  }  
162 + }
163 +
164 + void Molecule::atoms2rigidBodies( void ) {
165 +  int i;
166 +  for (i = 0; i < nRigidBodies; i++) {
167 +    if (myRigidBodies[i] != NULL) {
168 +      myRigidBodies[i]->calcForcesAndTorques();  
169 +    }
170    }
171   }
172  
173   void Molecule::getCOM( double COM[3] ) {
174  
175    double mass, mtot;
176 <  int i;
176 >  double aPos[3];
177 >  int i, j;
178  
179 <  COM[0] = 0.0;
180 <  COM[1] = 0.0;
181 <  COM[2] = 0.0;
179 >  for (j=0; j<3; j++)
180 >    COM[j] = 0.0;
181 >
182    mtot   = 0.0;
183  
184    for (i=0; i < nAtoms; i++) {
# Line 153 | Line 186 | void Molecule::getCOM( double COM[3] ) {
186  
187        mass = myAtoms[i]->getMass();
188        mtot   += mass;
189 <      COM[0] += myAtoms[i]->getX() * mass;
190 <      COM[1] += myAtoms[i]->getY() * mass;
158 <      COM[2] += myAtoms[i]->getZ() * mass;
189 >      
190 >      myAtoms[i]->getPos( aPos );
191  
192 +      for( j = 0; j < 3; j++)
193 +        COM[j] += aPos[j] * mass;
194 +
195      }
196    }
197  
198 <  COM[0] /= mtot;
199 <  COM[1] /= mtot;
165 <  COM[2] /= mtot;
166 <
198 >  for (j = 0; j < 3; j++)
199 >    COM[j] /= mtot;
200   }
201  
202   double Molecule::getCOMvel( double COMvel[3] ) {
203  
204    double mass, mtot;
205 <  int i;
205 >  double aVel[3];
206 >  int i, j;
207  
208 <  COMvel[0] = 0.0;
209 <  COMvel[1] = 0.0;
210 <  COMvel[2] = 0.0;
208 >
209 >  for (j=0; j<3; j++)
210 >    COMvel[j] = 0.0;
211 >
212    mtot   = 0.0;
213  
214    for (i=0; i < nAtoms; i++) {
# Line 181 | Line 216 | double Molecule::getCOMvel( double COMvel[3] ) {
216  
217        mass = myAtoms[i]->getMass();
218        mtot   += mass;
184 –      COMvel[0] += myAtoms[i]->get_vx() * mass;
185 –      COMvel[1] += myAtoms[i]->get_vy() * mass;
186 –      COMvel[2] += myAtoms[i]->get_vz() * mass;
219  
220 +      myAtoms[i]->getVel(aVel);
221 +
222 +      for (j=0; j<3; j++)
223 +        COMvel[j] += aVel[j]*mass;
224 +
225      }
226    }
227  
228 <  COMvel[0] /= mtot;
229 <  COMvel[1] /= mtot;
193 <  COMvel[2] /= mtot;
228 >  for (j=0; j<3; j++)
229 >    COMvel[j] /= mtot;
230  
231    return mtot;
232  
233   }
234  
235 < void Molecule::atomicRollCall(int* molMembership) {
236 <  int i, which;
237 <
238 <  for (i=0; i < nAtoms; i++) {
239 <    if (myAtoms[i] != NULL) {
240 <      which = myAtoms[i]->getIndex();
241 <      molMembership[which] = globalIndex;
242 <    }
235 > double Molecule::getTotalMass()
236 > {
237 >  int natoms;
238 >  Atom** atoms;
239 >  double totalMass;
240 >  
241 >  natoms = getNAtoms();
242 >  atoms = getMyAtoms();
243 >  totalMass = 0;
244 >  for(int i =0; i < natoms; i++){
245 >    totalMass += atoms[i]->getMass();
246    }
247 +
248 +  return totalMass;
249   }

Diff Legend

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