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 829 by gezelter, Tue Oct 28 16:03:37 2003 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 <
15 >  hasMassRatio = false;
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;
59 <    
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  
67   void Molecule::calcForces( void ){
68    
69    int i;
70 +  double com[3];
71  
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 80 | Line 99 | void Molecule::calcForces( void ){
99    for(i=0; i<nTorsions; i++){
100      myTorsions[i]->calc_forces();
101    }
102 +
103 +  // Rigid Body forces and torques are done after the fortran force loop
104 +
105   }
106  
107  
# Line 87 | 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 118 | Line 144 | void Molecule::printMe( void ){
144    for(i=0; i<nTorsions; i++){
145      myTorsions[i]->printMe();
146    }
147 +
148   }
149  
150   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<myRigidBodies.size(); i++) {
167 +
168 +      myRigidBodies[i]->getPos( aPos );
169 +
170 +      for (j=0; j< 3; j++)
171 +        aPos[j] += delta[j];
172 +      
173 +      myRigidBodies[i]->setPos( aPos );
174 +    }
175   }
176  
177 + void Molecule::atoms2rigidBodies( void ) {
178 +  int i;
179 +  for (i = 0; i < myRigidBodies.size(); i++) {
180 +    myRigidBodies[i]->calcForcesAndTorques();  
181 +  }
182 + }
183 +
184   void Molecule::getCOM( double COM[3] ) {
185  
186    double mass, mtot;
# Line 148 | 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 178 | 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 201 | Line 245 | double Molecule::getTotalMass()
245  
246   double Molecule::getTotalMass()
247   {
248 <  int natoms;
205 <  Atom** atoms;
248 >
249    double totalMass;
250    
208 –  natoms = getNAtoms();
209 –  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