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 726 by tim, Tue Aug 26 20:37:30 2003 UTC vs.
Revision 1180 by chrisfen, Thu May 20 20:24:07 2004 UTC

# Line 1 | Line 1
1   #include <iostream>
2 < #include <cstdlib>
3 < #include <cmath>
2 > #include <stdlib.h>
3 > #include <math.h>
4  
5   #ifdef IS_MPI
6   #include "mpiSimulation.hpp"
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 –  info->the_integrator = this;
32  
33    nAtoms = info->n_atoms;
34 <
34 >  integrableObjects = info->integrableObjects;
35 >
36    // check for constraints
37  
38    constrainedA = NULL;
# Line 41 | 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 65 | 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 111 | Line 119 | template<typename T> void Integrator<T>::checkConstrai
119      }
120    }
121  
122 +
123    if (nConstrained > 0){
124      isConstrained = 1;
125  
# Line 132 | 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 147 | Line 156 | template<typename T> void Integrator<T>::integrate(voi
156  
157  
158   template<typename T> void Integrator<T>::integrate(void){
150 –  int i, j;                         // loop counters
159  
160    double runTime = info->run_time;
161    double sampleTime = info->sampleTime;
162    double statusTime = info->statusTime;
163    double thermalTime = info->thermalTime;
164 +  double resetTime = info->resetTime;
165  
166 +  double difference;
167    double currSample;
168    double currThermal;
169    double currStatus;
170 +  double currReset;
171  
172    int calcPot, calcStress;
162 –  int isError;
173  
174    tStats = new Thermo(info);
175    statOut = new StatWriter(info);
176    dumpOut = new DumpWriter(info);
177  
178    atoms = info->atoms;
169 –  DirectionalAtom* dAtom;
179  
180    dt = info->dt;
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 <  // myFF->doForces(1,1);
198 <
197 >  
198 >  if (nConstrained){
199 >    preMove();
200 >    constrainA();
201 >    calcForce(1, 1);
202 >    constrainB();
203 >  }
204 >  
205    if (info->setTemp){
206      thermalize();
207    }
208  
183 –  calcPot = 0;
184 –  calcStress = 0;
185 –  currSample = sampleTime;
186 –  currThermal = thermalTime;
187 –  currStatus = statusTime;
188 –  
209    calcPot     = 0;
210    calcStress  = 0;
211    currSample  = sampleTime + info->getTime();
212    currThermal = thermalTime+ info->getTime();
213    currStatus  = statusTime + info->getTime();
214 +  currReset   = resetTime  + info->getTime();
215  
216    dumpOut->writeDump(info->getTime());
217    statOut->writeStat(info->getTime());
218  
198 –  readyCheck();
219  
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 225 | 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){
268 +        this->resetIntegrator();
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.");
279      MPIcheckPoint();
280   #endif // is_mpi
281    }
282  
283 <  dumpOut->writeFinal(info->getTime());
283 >  // dump out a file containing the omega values for the final configuration
284 >  if (info->useThermInt)
285 >    myFF->dumpzAngle();
286 >  
287  
288    delete dumpOut;
289    delete statOut;
# Line 246 | 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 + #ifdef PROFILE
303 +  endProfile(pro3);
304 +
305 +  startProfile(pro4);
306 + #endif // profile
307 +
308    moveA();
309  
310 <  if (nConstrained){
311 <    constrainA();
312 <  }
310 > #ifdef PROFILE
311 >  endProfile(pro4);
312 >  
313 >  startProfile(pro5);
314 > #endif//profile
315  
316  
317   #ifdef IS_MPI
# Line 270 | 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 <  if (nConstrained){
343 <    constrainB();
344 <  }
342 > #ifdef PROFILE
343 >  endProfile(pro6);
344 > #endif // profile
345  
346   #ifdef IS_MPI
347    strcpy(checkPointMsg, "Succesful moveB\n");
# Line 287 | 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];
293 –  double A[3][3], I[3][3];
294 –  double angle;
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  
298 –  for (i = 0; i < nAtoms; i++){
299 –    atoms[i]->getVel(vel);
300 –    atoms[i]->getPos(pos);
301 –    atoms[i]->getFrc(frc);
302 –
303 –    mass = atoms[i]->getMass();
304 –
367      for (j = 0; j < 3; j++){
368        // velocity half step
369        vel[j] += (dt2 * frc[j] / mass) * eConvert;
# Line 309 | 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()){
316 <      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 <      // use the angular velocities to propagate the rotation matrix a
331 <      // full time step
391 >      this->rotationPropagation( integrableObjects[i], ji );
392  
393 <      dAtom->getA(A);
334 <      dAtom->getI(I);
335 <
336 <      // rotate about the x-axis      
337 <      angle = dt2 * ji[0] / I[0][0];
338 <      this->rotate(1, 2, angle, ji, A);
339 <
340 <      // rotate about the y-axis
341 <      angle = dt2 * ji[1] / I[1][1];
342 <      this->rotate(2, 0, angle, ji, A);
343 <
344 <      // rotate about the z-axis
345 <      angle = dt * ji[2] / I[2][2];
346 <      this->rotate(0, 1, angle, ji, A);
347 <
348 <      // rotate about the y-axis
349 <      angle = dt2 * ji[1] / I[1][1];
350 <      this->rotate(2, 0, angle, ji, A);
351 <
352 <      // rotate about the x-axis
353 <      angle = dt2 * ji[0] / I[0][0];
354 <      this->rotate(1, 2, angle, ji, A);
355 <
356 <
357 <      dAtom->setJ(ji);
358 <      dAtom->setA(A);
393 >      integrableObjects[i]->setJ(ji);
394      }
395    }
396 +
397 +  if (nConstrained){
398 +    constrainA();
399 +  }
400   }
401  
402  
403   template<typename T> void Integrator<T>::moveB(void){
404    int i, j;
366 –  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()){
384 <      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 +
440 +  if (nConstrained){
441 +    constrainB();
442 +  }
443   }
444  
445   template<typename T> void Integrator<T>::preMove(void){
# Line 417 | Line 458 | template<typename T> void Integrator<T>::constrainA(){
458   }
459  
460   template<typename T> void Integrator<T>::constrainA(){
461 <  int i, j, k;
461 >  int i, j;
462    int done;
463    double posA[3], posB[3];
464    double velA[3], velB[3];
# Line 557 | Line 598 | template<typename T> void Integrator<T>::constrainA(){
598      painCave.isFatal = 1;
599      simError();
600    }
601 +
602   }
603  
604   template<typename T> void Integrator<T>::constrainB(void){
605 <  int i, j, k;
605 >  int i, j;
606    int done;
607    double posA[3], posB[3];
608    double velA[3], velB[3];
# Line 569 | Line 611 | template<typename T> void Integrator<T>::constrainB(vo
611    int a, b, ax, ay, az, bx, by, bz;
612    double rma, rmb;
613    double dx, dy, dz;
614 <  double rabsq, pabsq, rvab;
573 <  double diffsq;
614 >  double rvab;
615    double gab;
616    int iteration;
617  
# Line 660 | Line 701 | template<typename T> void Integrator<T>::constrainB(vo
701    }
702   }
703  
704 + template<typename T> void Integrator<T>::rotationPropagation
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 +  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,
757                                                  double angle, double ji[3],
758                                                  double A[3][3]){
# Line 725 | 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 750 | Line 843 | template<typename T> void Integrator<T>::thermalize(){
843   template<typename T> void Integrator<T>::thermalize(){
844    tStats->velocitize();
845   }
846 +
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