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 427 by mmeineke, Thu Mar 27 20:48:37 2003 UTC vs.
Revision 1140 by tim, Wed Apr 28 22:34:02 2004 UTC

# Line 1 | Line 1
1 < #include <cstdlib>
1 > #include <stdlib.h>
2  
3  
4   #include "Molecule.hpp"
5 + #include "simError.h"
6  
7  
8  
# Line 11 | Line 12 | Molecule::Molecule( void ){
12    myBonds = NULL;
13    myBends = NULL;
14    myTorsions = NULL;
15 <
15 >  hasMassRatio = false;
16   }
17  
17 –
18 –
18   Molecule::~Molecule( void ){
19    int i;
20    
# Line 39 | Line 38 | Molecule::~Molecule( void ){
38      delete[] myTorsions;
39    }
40  
42 –  if( myExcludes != NULL ){
43 –    for(i=0; i<nExcludes; i++) if(myExcludes[i] != NULL ) delete myExcludes[i];
44 –    delete[] myExcludes;
45 –  }
41   }
42  
43  
# Line 53 | 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.myExcludses;
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 79 | 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  
108 < void Molecule::getPotential( void ){
108 > 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 101 | Line 128 | void Molecule::getPotential( void ){
128  
129    return myPot;
130   }
131 +
132 + void Molecule::printMe( void ){
133 +  
134 +  int i;
135 +
136 +  for(i=0; i<nBonds; i++){
137 +    myBonds[i]->printMe();
138 +  }
139 +
140 +  for(i=0; i<nBends; i++){
141 +    myBends[i]->printMe();
142 +  }
143 +
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<myIntegrableObjects.size(); i++) {
155 +    if(myIntegrableObjects[i] != NULL ) {
156 +      
157 +      myIntegrableObjects[i]->getPos( aPos );
158 +      
159 +      for (j=0; j< 3; j++)
160 +        aPos[j] += delta[j];
161 +
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;
187 +  double aPos[3];
188 +  int i, j;
189 +
190 +  for (j=0; j<3; j++)
191 +    COM[j] = 0.0;
192 +
193 +  mtot   = 0.0;
194 +
195 +  for (i=0; i < myIntegrableObjects.size(); i++) {
196 +    if (myIntegrableObjects[i] != NULL) {
197 +
198 +      mass = myIntegrableObjects[i]->getMass();
199 +      mtot   += mass;
200 +      
201 +      myIntegrableObjects[i]->getPos( aPos );
202 +
203 +      for( j = 0; j < 3; j++)
204 +        COM[j] += aPos[j] * mass;
205 +
206 +    }
207 +  }
208 +
209 +  for (j = 0; j < 3; j++)
210 +    COM[j] /= mtot;
211 + }
212 +
213 + double Molecule::getCOMvel( double COMvel[3] ) {
214 +
215 +  double mass, mtot;
216 +  double aVel[3];
217 +  int i, j;
218 +
219 +
220 +  for (j=0; j<3; j++)
221 +    COMvel[j] = 0.0;
222 +
223 +  mtot   = 0.0;
224 +
225 +  for (i=0; i < myIntegrableObjects.size(); i++) {
226 +    if (myIntegrableObjects[i] != NULL) {
227 +
228 +      mass = myIntegrableObjects[i]->getMass();
229 +      mtot   += mass;
230 +
231 +      myIntegrableObjects[i]->getVel(aVel);
232 +
233 +      for (j=0; j<3; j++)
234 +        COMvel[j] += aVel[j]*mass;
235 +
236 +    }
237 +  }
238 +
239 +  for (j=0; j<3; j++)
240 +    COMvel[j] /= mtot;
241 +
242 +  return mtot;
243 +
244 + }
245 +
246 + double Molecule::getTotalMass()
247 + {
248 +
249 +  double totalMass;
250 +  
251 +  totalMass = 0;
252 +  for(int i =0; i < myIntegrableObjects.size(); i++){
253 +    totalMass += myIntegrableObjects[i]->getMass();
254 +  }
255 +
256 +  return totalMass;
257 + }

Diff Legend

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