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 483 by gezelter, Wed Apr 9 04:06:43 2003 UTC vs.
Revision 1109 by gezelter, Wed Apr 14 16:32:15 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  
# Line 54 | Line 48 | void Molecule::initialize( molInit &theInit ){
48    nBonds = theInit.nBonds;
49    nBends = theInit.nBends;
50    nTorsions = theInit.nTorsions;
51 <  nExcludes = theInit.nExcludes;
51 >  nRigidBodies = theInit.nRigidBodies;
52    nOriented = theInit.nOriented;
53  
54    myAtoms = theInit.myAtoms;
55    myBonds = theInit.myBonds;
56    myBends = theInit.myBends;
57    myTorsions = theInit.myTorsions;
58 <  myExcludes = theInit.myExcludes;
58 >  myRigidBodies = theInit.myRigidBodies;
59 >
60 >  myIntegrableObjects = theInit.myIntegrableObjects;
61      
62   }
63  
# Line 69 | Line 65 | void Molecule::calcForces( void ){
65    
66    int i;
67  
68 +  for(i=0; i<myRigidBodies.size(); i++) {
69 +    myRigidBodies[i]->updateAtoms();
70 +  }
71 +
72    for(i=0; i<nBonds; i++){
73      myBonds[i]->calc_forces();
74    }
# Line 80 | Line 80 | void Molecule::calcForces( void ){
80    for(i=0; i<nTorsions; i++){
81      myTorsions[i]->calc_forces();
82    }
83 +
84 +  // Rigid Body forces and torques are done after the fortran force loop
85 +
86   }
87  
88  
# Line 87 | Line 90 | double Molecule::getPotential( void ){
90    
91    int i;
92    double myPot = 0.0;
93 +
94 +  for(i=0; i<myRigidBodies.size(); i++) {
95 +    myRigidBodies[i]->updateAtoms();
96 +  }
97    
98    for(i=0; i<nBonds; i++){
99      myPot += myBonds[i]->get_potential();
# Line 118 | Line 125 | void Molecule::printMe( void ){
125    for(i=0; i<nTorsions; i++){
126      myTorsions[i]->printMe();
127    }
128 +
129   }
130  
131   void Molecule::moveCOM(double delta[3]){
132 <  double x, y, z;
133 <  int i;
132 >  double aPos[3];
133 >  int i, j;
134  
135    for(i=0; i<nAtoms; i++) {
136      if(myAtoms[i] != NULL ) {
137 +      
138 +      myAtoms[i]->getPos( aPos );
139 +      
140 +      for (j=0; j< 3; j++)
141 +        aPos[j] += delta[j];
142  
143 <      x = myAtoms[i]->getX() + delta[0];
144 <      y = myAtoms[i]->getY() + delta[1];
145 <      z = myAtoms[i]->getZ() + delta[2];
143 >      myAtoms[i]->setPos( aPos );
144 >    }
145 >  }
146  
147 <      myAtoms[i]->setX(x);
148 <      myAtoms[i]->setY(y);
149 <      myAtoms[i]->setZ(z);
147 >  for(i=0; i<myRigidBodies.size(); i++) {
148 >
149 >      myRigidBodies[i]->getPos( aPos );
150 >
151 >      for (j=0; j< 3; j++)
152 >        aPos[j] += delta[j];
153 >      
154 >      myRigidBodies[i]->setPos( aPos );
155      }
156 + }
157 +
158 + void Molecule::atoms2rigidBodies( void ) {
159 +  int i;
160 +  for (i = 0; i < myRigidBodies.size(); i++) {
161 +    myRigidBodies[i]->calcForcesAndTorques();  
162    }
163   }
164  
165   void Molecule::getCOM( double COM[3] ) {
166  
167    double mass, mtot;
168 <  int i;
168 >  double aPos[3];
169 >  int i, j;
170  
171 <  COM[0] = 0.0;
172 <  COM[1] = 0.0;
173 <  COM[2] = 0.0;
171 >  for (j=0; j<3; j++)
172 >    COM[j] = 0.0;
173 >
174    mtot   = 0.0;
175  
176    for (i=0; i < nAtoms; i++) {
# Line 153 | Line 178 | void Molecule::getCOM( double COM[3] ) {
178  
179        mass = myAtoms[i]->getMass();
180        mtot   += mass;
181 <      COM[0] += myAtoms[i]->getX() * mass;
182 <      COM[1] += myAtoms[i]->getY() * mass;
158 <      COM[2] += myAtoms[i]->getZ() * mass;
181 >      
182 >      myAtoms[i]->getPos( aPos );
183  
184 +      for( j = 0; j < 3; j++)
185 +        COM[j] += aPos[j] * mass;
186 +
187      }
188    }
189  
190 <  COM[0] /= mtot;
191 <  COM[1] /= mtot;
165 <  COM[2] /= mtot;
166 <
190 >  for (j = 0; j < 3; j++)
191 >    COM[j] /= mtot;
192   }
193  
194   double Molecule::getCOMvel( double COMvel[3] ) {
195  
196    double mass, mtot;
197 <  int i;
197 >  double aVel[3];
198 >  int i, j;
199  
200 <  COMvel[0] = 0.0;
201 <  COMvel[1] = 0.0;
202 <  COMvel[2] = 0.0;
200 >
201 >  for (j=0; j<3; j++)
202 >    COMvel[j] = 0.0;
203 >
204    mtot   = 0.0;
205  
206    for (i=0; i < nAtoms; i++) {
207      if (myAtoms[i] != NULL) {
208 <      
208 >
209        mass = myAtoms[i]->getMass();
210        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;
211  
212 +      myAtoms[i]->getVel(aVel);
213 +
214 +      for (j=0; j<3; j++)
215 +        COMvel[j] += aVel[j]*mass;
216 +
217      }
218    }
219  
220 <  COMvel[0] /= mtot;
221 <  COMvel[1] /= mtot;
193 <  COMvel[2] /= mtot;
220 >  for (j=0; j<3; j++)
221 >    COMvel[j] /= mtot;
222  
223    return mtot;
224  
225   }
226  
227 < void Molecule::atomicRollCall(int* molMembership) {
228 <  int i, which;
229 <
230 <  for (i=0; i < nAtoms; i++) {
231 <    if (myAtoms[i] != NULL) {
232 <      which = myAtoms[i]->getIndex();
233 <      molMembership[which] = myIndex;
234 <    }
227 > double Molecule::getTotalMass()
228 > {
229 >  int natoms;
230 >  Atom** atoms;
231 >  double totalMass;
232 >  
233 >  natoms = getNAtoms();
234 >  atoms = getMyAtoms();
235 >  totalMass = 0;
236 >  for(int i =0; i < natoms; i++){
237 >    totalMass += atoms[i]->getMass();
238    }
239 +
240 +  return totalMass;
241   }

Diff Legend

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