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 658 by tim, Thu Jul 31 15:35:07 2003 UTC vs.
Revision 1136 by tim, Tue Apr 27 16:26:44 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  
44   void Molecule::initialize( molInit &theInit ){
45 <
45 >  double totMass;
46 >  
47    nAtoms = theInit.nAtoms;
48    nMembers = nAtoms;
49    nBonds = theInit.nBonds;
50    nBends = theInit.nBends;
51    nTorsions = theInit.nTorsions;
52 <  nExcludes = theInit.nExcludes;
52 >  nRigidBodies = theInit.nRigidBodies;
53    nOriented = theInit.nOriented;
54  
55    myAtoms = theInit.myAtoms;
56    myBonds = theInit.myBonds;
57    myBends = theInit.myBends;
58    myTorsions = theInit.myTorsions;
59 <  myExcludes = theInit.myExcludes;
60 <    
59 >  myRigidBodies = theInit.myRigidBodies;
60 >
61 >  myIntegrableObjects = theInit.myIntegrableObjects;
62 >
63 >  for (int i = 0; i < myRigidBodies.size(); i++)
64 >      myRigidBodies[i]->calcRefCoords();
65 >
66 >
67 >  //the mass ratio will never change during the simulation. Thus, we could
68 >  //just calculate it at the begining of the simulation
69 >  totMass = getTotalMass();
70 >  for(int i = 0; i < nAtoms; i ++)
71 >    myAtoms[i]->setMassRatio(myAtoms[i]->getMass()/totMass);  
72   }
73  
74   void Molecule::calcForces( void ){
75    
76    int i;
77 +  double com[3];
78  
79 +  for(i=0; i<myRigidBodies.size(); i++) {
80 +    myRigidBodies[i]->updateAtoms();
81 +  }
82 +
83 +  //calculate the center of mass of the molecule
84 +  getCOM(com);  
85 +  for(int i = 0; i < nAtoms; i ++)
86 +    myAtoms[i]->setRc(com);  
87 +  
88 +
89    for(i=0; i<nBonds; i++){
90      myBonds[i]->calc_forces();
91    }
# Line 80 | Line 97 | void Molecule::calcForces( void ){
97    for(i=0; i<nTorsions; i++){
98      myTorsions[i]->calc_forces();
99    }
100 +
101 +  // Rigid Body forces and torques are done after the fortran force loop
102 +
103   }
104  
105  
# Line 87 | Line 107 | double Molecule::getPotential( void ){
107    
108    int i;
109    double myPot = 0.0;
110 +
111 +  for(i=0; i<myRigidBodies.size(); i++) {
112 +    myRigidBodies[i]->updateAtoms();
113 +  }
114    
115    for(i=0; i<nBonds; i++){
116      myPot += myBonds[i]->get_potential();
# Line 118 | Line 142 | void Molecule::printMe( void ){
142    for(i=0; i<nTorsions; i++){
143      myTorsions[i]->printMe();
144    }
145 +
146   }
147  
148   void Molecule::moveCOM(double delta[3]){
149    double aPos[3];
150    int i, j;
151  
152 <  for(i=0; i<nAtoms; i++) {
153 <    if(myAtoms[i] != NULL ) {
152 >  for(i=0; i<myIntegrableObjects.size(); i++) {
153 >    if(myIntegrableObjects[i] != NULL ) {
154        
155 <      myAtoms[i]->getPos( aPos );
155 >      myIntegrableObjects[i]->getPos( aPos );
156        
157        for (j=0; j< 3; j++)
158          aPos[j] += delta[j];
159  
160 <      myAtoms[i]->setPos( aPos );
160 >      myIntegrableObjects[i]->setPos( aPos );
161      }
162    }
163 +
164 +  for(i=0; i<myRigidBodies.size(); i++) {
165 +
166 +      myRigidBodies[i]->getPos( aPos );
167 +
168 +      for (j=0; j< 3; j++)
169 +        aPos[j] += delta[j];
170 +      
171 +      myRigidBodies[i]->setPos( aPos );
172 +    }
173   }
174  
175 + void Molecule::atoms2rigidBodies( void ) {
176 +  int i;
177 +  for (i = 0; i < myRigidBodies.size(); i++) {
178 +    myRigidBodies[i]->calcForcesAndTorques();  
179 +  }
180 + }
181 +
182   void Molecule::getCOM( double COM[3] ) {
183  
184    double mass, mtot;
# Line 148 | Line 190 | void Molecule::getCOM( double COM[3] ) {
190  
191    mtot   = 0.0;
192  
193 <  for (i=0; i < nAtoms; i++) {
194 <    if (myAtoms[i] != NULL) {
193 >  for (i=0; i < myIntegrableObjects.size(); i++) {
194 >    if (myIntegrableObjects[i] != NULL) {
195  
196 <      mass = myAtoms[i]->getMass();
196 >      mass = myIntegrableObjects[i]->getMass();
197        mtot   += mass;
198        
199 <      myAtoms[i]->getPos( aPos );
199 >      myIntegrableObjects[i]->getPos( aPos );
200  
201        for( j = 0; j < 3; j++)
202          COM[j] += aPos[j] * mass;
# Line 178 | Line 220 | double Molecule::getCOMvel( double COMvel[3] ) {
220  
221    mtot   = 0.0;
222  
223 <  for (i=0; i < nAtoms; i++) {
224 <    if (myAtoms[i] != NULL) {
223 >  for (i=0; i < myIntegrableObjects.size(); i++) {
224 >    if (myIntegrableObjects[i] != NULL) {
225  
226 <      mass = myAtoms[i]->getMass();
226 >      mass = myIntegrableObjects[i]->getMass();
227        mtot   += mass;
228  
229 <      myAtoms[i]->getVel(aVel);
229 >      myIntegrableObjects[i]->getVel(aVel);
230  
231        for (j=0; j<3; j++)
232          COMvel[j] += aVel[j]*mass;
# Line 201 | Line 243 | double Molecule::getTotalMass()
243  
244   double Molecule::getTotalMass()
245   {
246 <  int natoms;
205 <  Atom** atoms;
246 >
247    double totalMass;
248    
208 –  natoms = getNAtoms();
209 –  atoms = getMyAtoms();
249    totalMass = 0;
250 <  for(int i =0; i < natoms; i++){
251 <    totalMass += atoms[i]->getMass();
250 >  for(int i =0; i < myIntegrableObjects.size(); i++){
251 >    totalMass += myIntegrableObjects[i]->getMass();
252    }
253  
254    return totalMass;

Diff Legend

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