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 446 by gezelter, Thu Apr 3 20:19:50 2003 UTC vs.
Revision 1113 by tim, Thu Apr 15 16:18:26 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 >  for (int i = 0; i < myRigidBodies.size(); i++)
63 >      myRigidBodies[i]->calcRefCoords();
64      
65   }
66  
# Line 69 | Line 68 | void Molecule::calcForces( void ){
68    
69    int i;
70  
71 +  for(i=0; i<myRigidBodies.size(); i++) {
72 +    myRigidBodies[i]->updateAtoms();
73 +  }
74 +
75    for(i=0; i<nBonds; i++){
76      myBonds[i]->calc_forces();
77    }
# Line 80 | Line 83 | void Molecule::calcForces( void ){
83    for(i=0; i<nTorsions; i++){
84      myTorsions[i]->calc_forces();
85    }
86 +
87 +  // Rigid Body forces and torques are done after the fortran force loop
88 +
89   }
90  
91  
# Line 87 | Line 93 | double Molecule::getPotential( void ){
93    
94    int i;
95    double myPot = 0.0;
96 +
97 +  for(i=0; i<myRigidBodies.size(); i++) {
98 +    myRigidBodies[i]->updateAtoms();
99 +  }
100    
101    for(i=0; i<nBonds; i++){
102      myPot += myBonds[i]->get_potential();
# 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){
135 <  double x, y, z;
134 > void Molecule::moveCOM(double delta[3]){
135 >  double aPos[3];
136 >  int i, j;
137 >
138 >  for(i=0; i<myIntegrableObjects.size(); i++) {
139 >    if(myIntegrableObjects[i] != NULL ) {
140 >      
141 >      myIntegrableObjects[i]->getPos( aPos );
142 >      
143 >      for (j=0; j< 3; j++)
144 >        aPos[j] += delta[j];
145 >
146 >      myIntegrableObjects[i]->setPos( aPos );
147 >    }
148 >  }
149 >
150 >  for(i=0; i<myRigidBodies.size(); i++) {
151 >
152 >      myRigidBodies[i]->getPos( aPos );
153 >
154 >      for (j=0; j< 3; j++)
155 >        aPos[j] += delta[j];
156 >      
157 >      myRigidBodies[i]->setPos( aPos );
158 >    }
159 > }
160 >
161 > void Molecule::atoms2rigidBodies( void ) {
162    int i;
163 +  for (i = 0; i < myRigidBodies.size(); i++) {
164 +    myRigidBodies[i]->calcForcesAndTorques();  
165 +  }
166 + }
167  
168 <  for(i=0; i<nAtoms; i++) {
128 <    if(myAtoms[i] != NULL ) {
168 > void Molecule::getCOM( double COM[3] ) {
169  
170 <      x = myAtoms[i]->getX() + delta[0];
171 <      y = myAtoms[i]->getY() + delta[1];
172 <      z = myAtoms[i]->getZ() + delta[2];
170 >  double mass, mtot;
171 >  double aPos[3];
172 >  int i, j;
173  
174 <      myAtoms[i]->setX(x);
175 <      myAtoms[i]->setY(y);
176 <      myAtoms[i]->setZ(z);
174 >  for (j=0; j<3; j++)
175 >    COM[j] = 0.0;
176 >
177 >  mtot   = 0.0;
178 >
179 >  for (i=0; i < myIntegrableObjects.size(); i++) {
180 >    if (myIntegrableObjects[i] != NULL) {
181 >
182 >      mass = myIntegrableObjects[i]->getMass();
183 >      mtot   += mass;
184 >      
185 >      myIntegrableObjects[i]->getPos( aPos );
186 >
187 >      for( j = 0; j < 3; j++)
188 >        COM[j] += aPos[j] * mass;
189 >
190      }
191    }
192 +
193 +  for (j = 0; j < 3; j++)
194 +    COM[j] /= mtot;
195   }
196  
197 < double* Molecule::getCOM() {
197 > double Molecule::getCOMvel( double COMvel[3] ) {
198  
199    double mass, mtot;
200 <  int i;
200 >  double aVel[3];
201 >  int i, j;
202  
203 <  COM[0] = 0.0;
204 <  COM[1] = 0.0;
205 <  COM[2] = 0.0;
203 >
204 >  for (j=0; j<3; j++)
205 >    COMvel[j] = 0.0;
206 >
207    mtot   = 0.0;
208  
209 <  for (i=0; i < nAtoms; i++) {
210 <    if (myAtoms[i] != NULL) {
209 >  for (i=0; i < myIntegrableObjects.size(); i++) {
210 >    if (myIntegrableObjects[i] != NULL) {
211  
212 <      mass = myAtoms[i]->getMass();
212 >      mass = myIntegrableObjects[i]->getMass();
213        mtot   += mass;
156 –      COM[0] += myAtoms[i]->getX() * mass;
157 –      COM[1] += myAtoms[i]->getY() * mass;
158 –      COM[2] += myAtoms[i]->getZ() * mass;
214  
215 +      myIntegrableObjects[i]->getVel(aVel);
216 +
217 +      for (j=0; j<3; j++)
218 +        COMvel[j] += aVel[j]*mass;
219 +
220      }
221    }
222  
223 <  COM[0] /= mtot;
224 <  COM[1] /= mtot;
165 <  COM[2] /= mtot;
223 >  for (j=0; j<3; j++)
224 >    COMvel[j] /= mtot;
225  
226 <  return COM;
226 >  return mtot;
227 >
228   }
229 +
230 + double Molecule::getTotalMass()
231 + {
232 +
233 +  double totalMass;
234 +  
235 +  totalMass = 0;
236 +  for(int i =0; i < myIntegrableObjects.size(); i++){
237 +    totalMass += myIntegrableObjects[i]->getMass();
238 +  }
239 +
240 +  return totalMass;
241 + }

Diff Legend

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