ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Integrator.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Integrator.cpp (file contents):
Revision 829 by gezelter, Tue Oct 28 16:03:37 2003 UTC vs.
Revision 1180 by chrisfen, Thu May 20 20:24:07 2004 UTC

# Line 7 | Line 7
7   #include <unistd.h>
8   #endif //is_mpi
9  
10 + #ifdef PROFILE
11 + #include "mdProfile.hpp"
12 + #endif // profile
13 +
14   #include "Integrator.hpp"
15   #include "simError.h"
16  
# Line 25 | Line 29 | template<typename T> Integrator<T>::Integrator(SimInfo
29    if (info->the_integrator != NULL){
30      delete info->the_integrator;
31    }
28 –  
29 –  nAtoms = info->n_atoms;
32  
33 +  nAtoms = info->n_atoms;
34 +  integrableObjects = info->integrableObjects;
35 +
36    // check for constraints
37  
38    constrainedA = NULL;
# Line 40 | Line 45 | template<typename T> Integrator<T>::Integrator(SimInfo
45    nConstrained = 0;
46  
47    checkConstraints();
48 +
49 +  for (i=0; i<nMols; i++)
50 +    zAngle[i] = 0.0;
51   }
52  
53   template<typename T> Integrator<T>::~Integrator(){
# Line 64 | Line 72 | template<typename T> void Integrator<T>::checkConstrai
72  
73    SRI** theArray;
74    for (int i = 0; i < nMols; i++){
75 <    theArray = (SRI * *) molecules[i].getMyBonds();
75 >
76 >          theArray = (SRI * *) molecules[i].getMyBonds();
77      for (int j = 0; j < molecules[i].getNBonds(); j++){
78        constrained = theArray[j]->is_constrained();
79  
# Line 110 | Line 119 | template<typename T> void Integrator<T>::checkConstrai
119      }
120    }
121  
122 +
123    if (nConstrained > 0){
124      isConstrained = 1;
125  
# Line 131 | Line 141 | template<typename T> void Integrator<T>::checkConstrai
141      }
142  
143  
144 <    // save oldAtoms to check for lode balanceing later on.
144 >    // save oldAtoms to check for lode balancing later on.
145  
146      oldAtoms = nAtoms;
147  
# Line 153 | Line 163 | template<typename T> void Integrator<T>::integrate(voi
163    double thermalTime = info->thermalTime;
164    double resetTime = info->resetTime;
165  
166 <
166 >  double difference;
167    double currSample;
168    double currThermal;
169    double currStatus;
170    double currReset;
171 <  
171 >
172    int calcPot, calcStress;
173  
174    tStats = new Thermo(info);
# Line 171 | Line 181 | template<typename T> void Integrator<T>::integrate(voi
181    dt2 = 0.5 * dt;
182  
183    readyCheck();
184 +
185 +  // remove center of mass drift velocity (in case we passed in a configuration
186 +  // that was drifting
187 +  tStats->removeCOMdrift();
188 +
189 +  // initialize the retraints if necessary
190 +  if (info->useThermInt) {
191 +    myFF->initRestraints();
192 +  }
193  
194    // initialize the forces before the first step
195  
196    calcForce(1, 1);
197 <
197 >  
198    if (nConstrained){
199      preMove();
200      constrainA();
201 <    calcForce(1, 1);    
201 >    calcForce(1, 1);
202      constrainB();
203    }
204    
# Line 198 | Line 217 | template<typename T> void Integrator<T>::integrate(voi
217    statOut->writeStat(info->getTime());
218  
219  
201 –
220   #ifdef IS_MPI
221    strcpy(checkPointMsg, "The integrator is ready to go.");
222    MPIcheckPoint();
223   #endif // is_mpi
224  
225 <  while (info->getTime() < runTime){
226 <    if ((info->getTime() + dt) >= currStatus){
225 >  while (info->getTime() < runTime && !stopIntegrator()){
226 >    difference = info->getTime() + dt - currStatus;
227 >    if (difference > 0 || fabs(difference) < 1e-4 ){
228        calcPot = 1;
229        calcStress = 1;
230      }
231  
232 + #ifdef PROFILE
233 +    startProfile( pro1 );
234 + #endif
235 +    
236      integrateStep(calcPot, calcStress);
237  
238 + #ifdef PROFILE
239 +    endProfile( pro1 );
240 +
241 +    startProfile( pro2 );
242 + #endif // profile
243 +
244      info->incrTime(dt);
245  
246      if (info->setTemp){
# Line 227 | Line 256 | template<typename T> void Integrator<T>::integrate(voi
256      }
257  
258      if (info->getTime() >= currStatus){
259 <      statOut->writeStat(info->getTime());
260 <      calcPot = 0;
259 >      statOut->writeStat(info->getTime());
260 >      statOut->writeRaw(info->getTime());
261 >      calcPot = 0;
262        calcStress = 0;
263        currStatus += statusTime;
264 <    }
264 >    }
265  
266      if (info->resetIntegrator){
267        if (info->getTime() >= currReset){
# Line 239 | Line 269 | template<typename T> void Integrator<T>::integrate(voi
269          currReset += resetTime;
270        }
271      }
272 +    
273 + #ifdef PROFILE
274 +    endProfile( pro2 );
275 + #endif //profile
276  
277   #ifdef IS_MPI
278      strcpy(checkPointMsg, "successfully took a time step.");
# Line 246 | Line 280 | template<typename T> void Integrator<T>::integrate(voi
280   #endif // is_mpi
281    }
282  
283 +  // dump out a file containing the omega values for the final configuration
284 +  if (info->useThermInt)
285 +    myFF->dumpzAngle();
286 +  
287  
250 –  // write the last frame
251 –  dumpOut->writeDump(info->getTime());
252 –
288    delete dumpOut;
289    delete statOut;
290   }
# Line 257 | Line 292 | template<typename T> void Integrator<T>::integrateStep
292   template<typename T> void Integrator<T>::integrateStep(int calcPot,
293                                                         int calcStress){
294    // Position full step, and velocity half step
295 +
296 + #ifdef PROFILE
297 +  startProfile(pro3);
298 + #endif //profile
299 +
300    preMove();
301  
302 <  moveA();
302 > #ifdef PROFILE
303 >  endProfile(pro3);
304  
305 +  startProfile(pro4);
306 + #endif // profile
307  
308 +  moveA();
309  
310 + #ifdef PROFILE
311 +  endProfile(pro4);
312 +  
313 +  startProfile(pro5);
314 + #endif//profile
315  
316 +
317   #ifdef IS_MPI
318    strcpy(checkPointMsg, "Succesful moveA\n");
319    MPIcheckPoint();
# Line 279 | Line 329 | template<typename T> void Integrator<T>::integrateStep
329    MPIcheckPoint();
330   #endif // is_mpi
331  
332 + #ifdef PROFILE
333 +  endProfile( pro5 );
334  
335 +  startProfile( pro6 );
336 + #endif //profile
337 +
338    // finish the velocity  half step
339  
340    moveB();
341  
342 + #ifdef PROFILE
343 +  endProfile(pro6);
344 + #endif // profile
345  
288 –
346   #ifdef IS_MPI
347    strcpy(checkPointMsg, "Succesful moveB\n");
348    MPIcheckPoint();
# Line 294 | Line 351 | template<typename T> void Integrator<T>::moveA(void){
351  
352  
353   template<typename T> void Integrator<T>::moveA(void){
354 <  int i, j;
354 >  size_t i, j;
355    DirectionalAtom* dAtom;
356    double Tb[3], ji[3];
357    double vel[3], pos[3], frc[3];
358    double mass;
359 +
360 +  for (i = 0; i < integrableObjects.size() ; i++){
361 +    integrableObjects[i]->getVel(vel);
362 +    integrableObjects[i]->getPos(pos);
363 +    integrableObjects[i]->getFrc(frc);
364 +    
365 +    mass = integrableObjects[i]->getMass();
366  
303 –  for (i = 0; i < nAtoms; i++){
304 –    atoms[i]->getVel(vel);
305 –    atoms[i]->getPos(pos);
306 –    atoms[i]->getFrc(frc);
307 –
308 –    mass = atoms[i]->getMass();
309 –
367      for (j = 0; j < 3; j++){
368        // velocity half step
369        vel[j] += (dt2 * frc[j] / mass) * eConvert;
# Line 314 | Line 371 | template<typename T> void Integrator<T>::moveA(void){
371        pos[j] += dt * vel[j];
372      }
373  
374 <    atoms[i]->setVel(vel);
375 <    atoms[i]->setPos(pos);
374 >    integrableObjects[i]->setVel(vel);
375 >    integrableObjects[i]->setPos(pos);
376  
377 <    if (atoms[i]->isDirectional()){
321 <      dAtom = (DirectionalAtom *) atoms[i];
377 >    if (integrableObjects[i]->isDirectional()){
378  
379        // get and convert the torque to body frame
380  
381 <      dAtom->getTrq(Tb);
382 <      dAtom->lab2Body(Tb);
381 >      integrableObjects[i]->getTrq(Tb);
382 >      integrableObjects[i]->lab2Body(Tb);
383  
384        // get the angular momentum, and propagate a half step
385  
386 <      dAtom->getJ(ji);
386 >      integrableObjects[i]->getJ(ji);
387  
388        for (j = 0; j < 3; j++)
389          ji[j] += (dt2 * Tb[j]) * eConvert;
390  
391 <      this->rotationPropagation( dAtom, ji );
391 >      this->rotationPropagation( integrableObjects[i], ji );
392  
393 <      dAtom->setJ(ji);
393 >      integrableObjects[i]->setJ(ji);
394      }
395    }
396  
# Line 346 | Line 402 | template<typename T> void Integrator<T>::moveB(void){
402  
403   template<typename T> void Integrator<T>::moveB(void){
404    int i, j;
349 –  DirectionalAtom* dAtom;
405    double Tb[3], ji[3];
406    double vel[3], frc[3];
407    double mass;
408  
409 <  for (i = 0; i < nAtoms; i++){
410 <    atoms[i]->getVel(vel);
411 <    atoms[i]->getFrc(frc);
409 >  for (i = 0; i < integrableObjects.size(); i++){
410 >    integrableObjects[i]->getVel(vel);
411 >    integrableObjects[i]->getFrc(frc);
412  
413 <    mass = atoms[i]->getMass();
413 >    mass = integrableObjects[i]->getMass();
414  
415      // velocity half step
416      for (j = 0; j < 3; j++)
417        vel[j] += (dt2 * frc[j] / mass) * eConvert;
418  
419 <    atoms[i]->setVel(vel);
419 >    integrableObjects[i]->setVel(vel);
420  
421 <    if (atoms[i]->isDirectional()){
367 <      dAtom = (DirectionalAtom *) atoms[i];
421 >    if (integrableObjects[i]->isDirectional()){
422  
423 <      // get and convert the torque to body frame      
423 >      // get and convert the torque to body frame
424  
425 <      dAtom->getTrq(Tb);
426 <      dAtom->lab2Body(Tb);
425 >      integrableObjects[i]->getTrq(Tb);
426 >      integrableObjects[i]->lab2Body(Tb);
427  
428        // get the angular momentum, and propagate a half step
429  
430 <      dAtom->getJ(ji);
430 >      integrableObjects[i]->getJ(ji);
431  
432        for (j = 0; j < 3; j++)
433          ji[j] += (dt2 * Tb[j]) * eConvert;
434  
435  
436 <      dAtom->setJ(ji);
436 >      integrableObjects[i]->setJ(ji);
437      }
438    }
439  
# Line 648 | Line 702 | template<typename T> void Integrator<T>::rotationPropa
702   }
703  
704   template<typename T> void Integrator<T>::rotationPropagation
705 < ( DirectionalAtom* dAtom, double ji[3] ){
705 > ( StuntDouble* sd, double ji[3] ){
706  
707    double angle;
708    double A[3][3], I[3][3];
709 +  int i, j, k;
710  
711    // use the angular velocities to propagate the rotation matrix a
712    // full time step
713  
714 <  dAtom->getA(A);
715 <  dAtom->getI(I);
716 <  
717 <  // rotate about the x-axis      
718 <  angle = dt2 * ji[0] / I[0][0];
719 <  this->rotate( 1, 2, angle, ji, A );
720 <  
721 <  // rotate about the y-axis
722 <  angle = dt2 * ji[1] / I[1][1];
723 <  this->rotate( 2, 0, angle, ji, A );
724 <  
725 <  // rotate about the z-axis
726 <  angle = dt * ji[2] / I[2][2];
727 <  this->rotate( 0, 1, angle, ji, A);
728 <  
729 <  // rotate about the y-axis
730 <  angle = dt2 * ji[1] / I[1][1];
731 <  this->rotate( 2, 0, angle, ji, A );
732 <  
733 <  // rotate about the x-axis
734 <  angle = dt2 * ji[0] / I[0][0];
735 <  this->rotate( 1, 2, angle, ji, A );
736 <  
737 <  dAtom->setA( A  );    
714 >  sd->getA(A);
715 >  sd->getI(I);
716 >
717 >  if (sd->isLinear()) {
718 >    i = sd->linearAxis();
719 >    j = (i+1)%3;
720 >    k = (i+2)%3;
721 >    
722 >    angle = dt2 * ji[j] / I[j][j];
723 >    this->rotate( k, i, angle, ji, A );
724 >
725 >    angle = dt * ji[k] / I[k][k];
726 >    this->rotate( i, j, angle, ji, A);
727 >
728 >    angle = dt2 * ji[j] / I[j][j];
729 >    this->rotate( k, i, angle, ji, A );
730 >
731 >  } else {
732 >    // rotate about the x-axis
733 >    angle = dt2 * ji[0] / I[0][0];
734 >    this->rotate( 1, 2, angle, ji, A );
735 >    
736 >    // rotate about the y-axis
737 >    angle = dt2 * ji[1] / I[1][1];
738 >    this->rotate( 2, 0, angle, ji, A );
739 >    
740 >    // rotate about the z-axis
741 >    angle = dt * ji[2] / I[2][2];
742 >    this->rotate( 0, 1, angle, ji, A);
743 >    
744 >    // rotate about the y-axis
745 >    angle = dt2 * ji[1] / I[1][1];
746 >    this->rotate( 2, 0, angle, ji, A );
747 >    
748 >    // rotate about the x-axis
749 >    angle = dt2 * ji[0] / I[0][0];
750 >    this->rotate( 1, 2, angle, ji, A );
751 >    
752 >  }
753 >  sd->setA( A  );
754   }
755  
756   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
# Line 747 | Line 818 | template<typename T> void Integrator<T>::rotate(int ax
818      }
819    }
820  
821 <  // rotate the Rotation matrix acording to:
821 >  // rotate the Rotation matrix acording to:
822    //            A[][] = A[][] * transpose(rot[][])
823  
824  
# Line 776 | Line 847 | template<typename T> double Integrator<T>::getConserve
847   template<typename T> double Integrator<T>::getConservedQuantity(void){
848    return tStats->getTotalE();
849   }
850 + template<typename T> string Integrator<T>::getAdditionalParameters(void){
851 +  //By default, return a null string
852 +  //The reason we use string instead of char* is that if we use char*, we will
853 +  //return a pointer point to local variable which might cause problem
854 +  return string();
855 + }

Diff Legend

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