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 733 by tim, Wed Aug 27 19:23:29 2003 UTC vs.
Revision 1108 by tim, Wed Apr 14 15:37:41 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 +  integrableObjects = info->integrableObjects;
35  
36    // check for constraints
37  
# Line 65 | Line 69 | template<typename T> void Integrator<T>::checkConstrai
69  
70    SRI** theArray;
71    for (int i = 0; i < nMols; i++){
72 <    theArray = (SRI * *) molecules[i].getMyBonds();
72 >
73 >          theArray = (SRI * *) molecules[i].getMyBonds();
74      for (int j = 0; j < molecules[i].getNBonds(); j++){
75        constrained = theArray[j]->is_constrained();
76  
# Line 111 | Line 116 | template<typename T> void Integrator<T>::checkConstrai
116      }
117    }
118  
119 +
120    if (nConstrained > 0){
121      isConstrained = 1;
122  
# Line 132 | Line 138 | template<typename T> void Integrator<T>::checkConstrai
138      }
139  
140  
141 <    // save oldAtoms to check for lode balanceing later on.
141 >    // save oldAtoms to check for lode balancing later on.
142  
143      oldAtoms = nAtoms;
144  
# Line 147 | Line 153 | template<typename T> void Integrator<T>::integrate(voi
153  
154  
155   template<typename T> void Integrator<T>::integrate(void){
150 –  int i, j;                         // loop counters
156  
157    double runTime = info->run_time;
158    double sampleTime = info->sampleTime;
159    double statusTime = info->statusTime;
160    double thermalTime = info->thermalTime;
161 +  double resetTime = info->resetTime;
162  
163 +
164    double currSample;
165    double currThermal;
166    double currStatus;
167 +  double currReset;
168  
169    int calcPot, calcStress;
162 –  int isError;
170  
171    tStats = new Thermo(info);
172    statOut = new StatWriter(info);
173    dumpOut = new DumpWriter(info);
174  
175    atoms = info->atoms;
169 –  DirectionalAtom* dAtom;
176  
177    dt = info->dt;
178    dt2 = 0.5 * dt;
179  
180 +  readyCheck();
181 +
182    // initialize the forces before the first step
183  
184    calcForce(1, 1);
185 +
186 +  //temp test
187 +  tStats->getPotential();
188    
189 +  if (nConstrained){
190 +    preMove();
191 +    constrainA();
192 +    calcForce(1, 1);
193 +    constrainB();
194 +  }
195 +  
196    if (info->setTemp){
197      thermalize();
198    }
199  
182 –  calcPot = 0;
183 –  calcStress = 0;
184 –  currSample = sampleTime;
185 –  currThermal = thermalTime;
186 –  currStatus = statusTime;
187 –  
200    calcPot     = 0;
201    calcStress  = 0;
202    currSample  = sampleTime + info->getTime();
203    currThermal = thermalTime+ info->getTime();
204    currStatus  = statusTime + info->getTime();
205 +  currReset   = resetTime  + info->getTime();
206  
207    dumpOut->writeDump(info->getTime());
208    statOut->writeStat(info->getTime());
209  
197 –  readyCheck();
210  
211   #ifdef IS_MPI
212    strcpy(checkPointMsg, "The integrator is ready to go.");
213    MPIcheckPoint();
214   #endif // is_mpi
215  
216 <  while (info->getTime() < runTime){
216 >  while (info->getTime() < runTime && stopIntegrator()){
217      if ((info->getTime() + dt) >= currStatus){
218        calcPot = 1;
219        calcStress = 1;
220      }
221  
222 + #ifdef PROFILE
223 +    startProfile( pro1 );
224 + #endif
225 +    
226      integrateStep(calcPot, calcStress);
227  
228 + #ifdef PROFILE
229 +    endProfile( pro1 );
230 +
231 +    startProfile( pro2 );
232 + #endif // profile
233 +
234      info->incrTime(dt);
235  
236      if (info->setTemp){
# Line 224 | Line 246 | template<typename T> void Integrator<T>::integrate(voi
246      }
247  
248      if (info->getTime() >= currStatus){
249 <      statOut->writeStat(info->getTime());
250 <      calcPot = 0;
249 >      statOut->writeStat(info->getTime());
250 >      calcPot = 0;
251        calcStress = 0;
252        currStatus += statusTime;
253 <    }
253 >    }
254 >
255 >    if (info->resetIntegrator){
256 >      if (info->getTime() >= currReset){
257 >        this->resetIntegrator();
258 >        currReset += resetTime;
259 >      }
260 >    }
261 >    
262 > #ifdef PROFILE
263 >    endProfile( pro2 );
264 > #endif //profile
265  
266   #ifdef IS_MPI
267      strcpy(checkPointMsg, "successfully took a time step.");
# Line 236 | Line 269 | template<typename T> void Integrator<T>::integrate(voi
269   #endif // is_mpi
270    }
271  
239 –  dumpOut->writeFinal(info->getTime());
240 –
272    delete dumpOut;
273    delete statOut;
274   }
# Line 245 | Line 276 | template<typename T> void Integrator<T>::integrateStep
276   template<typename T> void Integrator<T>::integrateStep(int calcPot,
277                                                         int calcStress){
278    // Position full step, and velocity half step
279 +
280 + #ifdef PROFILE
281 +  startProfile(pro3);
282 + #endif //profile
283 +
284    preMove();
285  
286 + #ifdef PROFILE
287 +  endProfile(pro3);
288 +
289 +  startProfile(pro4);
290 + #endif // profile
291 +
292    moveA();
293  
294 <  if (nConstrained){
295 <    constrainA();
296 <  }
294 > #ifdef PROFILE
295 >  endProfile(pro4);
296 >  
297 >  startProfile(pro5);
298 > #endif//profile
299  
300  
301   #ifdef IS_MPI
# Line 269 | Line 313 | template<typename T> void Integrator<T>::integrateStep
313    MPIcheckPoint();
314   #endif // is_mpi
315  
316 + #ifdef PROFILE
317 +  endProfile( pro5 );
318  
319 +  startProfile( pro6 );
320 + #endif //profile
321 +
322    // finish the velocity  half step
323  
324    moveB();
325  
326 <  if (nConstrained){
327 <    constrainB();
328 <  }
326 > #ifdef PROFILE
327 >  endProfile(pro6);
328 > #endif // profile
329  
330   #ifdef IS_MPI
331    strcpy(checkPointMsg, "Succesful moveB\n");
# Line 286 | Line 335 | template<typename T> void Integrator<T>::moveA(void){
335  
336  
337   template<typename T> void Integrator<T>::moveA(void){
338 <  int i, j;
338 >  size_t i, j;
339    DirectionalAtom* dAtom;
340    double Tb[3], ji[3];
292 –  double A[3][3], I[3][3];
293 –  double angle;
341    double vel[3], pos[3], frc[3];
342    double mass;
343 +
344 +  for (i = 0; i < integrableObjects.size() ; i++){
345 +    integrableObjects[i]->getVel(vel);
346 +    integrableObjects[i]->getPos(pos);
347 +    integrableObjects[i]->getFrc(frc);
348 +    
349 +    mass = integrableObjects[i]->getMass();
350  
297 –  for (i = 0; i < nAtoms; i++){
298 –    atoms[i]->getVel(vel);
299 –    atoms[i]->getPos(pos);
300 –    atoms[i]->getFrc(frc);
301 –
302 –    mass = atoms[i]->getMass();
303 –
351      for (j = 0; j < 3; j++){
352        // velocity half step
353        vel[j] += (dt2 * frc[j] / mass) * eConvert;
# Line 308 | Line 355 | template<typename T> void Integrator<T>::moveA(void){
355        pos[j] += dt * vel[j];
356      }
357  
358 <    atoms[i]->setVel(vel);
359 <    atoms[i]->setPos(pos);
358 >    integrableObjects[i]->setVel(vel);
359 >    integrableObjects[i]->setPos(pos);
360  
361 <    if (atoms[i]->isDirectional()){
315 <      dAtom = (DirectionalAtom *) atoms[i];
361 >    if (integrableObjects[i]->isDirectional()){
362  
363        // get and convert the torque to body frame
364  
365 <      dAtom->getTrq(Tb);
366 <      dAtom->lab2Body(Tb);
365 >      integrableObjects[i]->getTrq(Tb);
366 >      integrableObjects[i]->lab2Body(Tb);
367  
368        // get the angular momentum, and propagate a half step
369  
370 <      dAtom->getJ(ji);
370 >      integrableObjects[i]->getJ(ji);
371  
372        for (j = 0; j < 3; j++)
373          ji[j] += (dt2 * Tb[j]) * eConvert;
374  
375 <      // use the angular velocities to propagate the rotation matrix a
330 <      // full time step
375 >      this->rotationPropagation( integrableObjects[i], ji );
376  
377 <      dAtom->getA(A);
333 <      dAtom->getI(I);
334 <
335 <      // rotate about the x-axis      
336 <      angle = dt2 * ji[0] / I[0][0];
337 <      this->rotate(1, 2, angle, ji, A);
338 <
339 <      // rotate about the y-axis
340 <      angle = dt2 * ji[1] / I[1][1];
341 <      this->rotate(2, 0, angle, ji, A);
342 <
343 <      // rotate about the z-axis
344 <      angle = dt * ji[2] / I[2][2];
345 <      this->rotate(0, 1, angle, ji, A);
346 <
347 <      // rotate about the y-axis
348 <      angle = dt2 * ji[1] / I[1][1];
349 <      this->rotate(2, 0, angle, ji, A);
350 <
351 <      // rotate about the x-axis
352 <      angle = dt2 * ji[0] / I[0][0];
353 <      this->rotate(1, 2, angle, ji, A);
354 <
355 <
356 <      dAtom->setJ(ji);
357 <      dAtom->setA(A);
377 >      integrableObjects[i]->setJ(ji);
378      }
379    }
380 +
381 +  if (nConstrained){
382 +    constrainA();
383 +  }
384   }
385  
386  
387   template<typename T> void Integrator<T>::moveB(void){
388    int i, j;
365 –  DirectionalAtom* dAtom;
389    double Tb[3], ji[3];
390    double vel[3], frc[3];
391    double mass;
392  
393 <  for (i = 0; i < nAtoms; i++){
394 <    atoms[i]->getVel(vel);
395 <    atoms[i]->getFrc(frc);
393 >  for (i = 0; i < integrableObjects.size(); i++){
394 >    integrableObjects[i]->getVel(vel);
395 >    integrableObjects[i]->getFrc(frc);
396  
397 <    mass = atoms[i]->getMass();
397 >    mass = integrableObjects[i]->getMass();
398  
399      // velocity half step
400      for (j = 0; j < 3; j++)
401        vel[j] += (dt2 * frc[j] / mass) * eConvert;
402  
403 <    atoms[i]->setVel(vel);
403 >    integrableObjects[i]->setVel(vel);
404  
405 <    if (atoms[i]->isDirectional()){
383 <      dAtom = (DirectionalAtom *) atoms[i];
405 >    if (integrableObjects[i]->isDirectional()){
406  
407 <      // get and convert the torque to body frame      
407 >      // get and convert the torque to body frame
408  
409 <      dAtom->getTrq(Tb);
410 <      dAtom->lab2Body(Tb);
409 >      integrableObjects[i]->getTrq(Tb);
410 >      integrableObjects[i]->lab2Body(Tb);
411  
412        // get the angular momentum, and propagate a half step
413  
414 <      dAtom->getJ(ji);
414 >      integrableObjects[i]->getJ(ji);
415  
416        for (j = 0; j < 3; j++)
417          ji[j] += (dt2 * Tb[j]) * eConvert;
418  
419  
420 <      dAtom->setJ(ji);
420 >      integrableObjects[i]->setJ(ji);
421      }
422    }
423 +
424 +  if (nConstrained){
425 +    constrainB();
426 +  }
427   }
428  
429   template<typename T> void Integrator<T>::preMove(void){
# Line 416 | Line 442 | template<typename T> void Integrator<T>::constrainA(){
442   }
443  
444   template<typename T> void Integrator<T>::constrainA(){
445 <  int i, j, k;
445 >  int i, j;
446    int done;
447    double posA[3], posB[3];
448    double velA[3], velB[3];
# Line 556 | Line 582 | template<typename T> void Integrator<T>::constrainA(){
582      painCave.isFatal = 1;
583      simError();
584    }
585 +
586   }
587  
588   template<typename T> void Integrator<T>::constrainB(void){
589 <  int i, j, k;
589 >  int i, j;
590    int done;
591    double posA[3], posB[3];
592    double velA[3], velB[3];
# Line 568 | Line 595 | template<typename T> void Integrator<T>::constrainB(vo
595    int a, b, ax, ay, az, bx, by, bz;
596    double rma, rmb;
597    double dx, dy, dz;
598 <  double rabsq, pabsq, rvab;
572 <  double diffsq;
598 >  double rvab;
599    double gab;
600    int iteration;
601  
# Line 659 | Line 685 | template<typename T> void Integrator<T>::rotate(int ax
685    }
686   }
687  
688 + template<typename T> void Integrator<T>::rotationPropagation
689 + ( StuntDouble* sd, double ji[3] ){
690 +
691 +  double angle;
692 +  double A[3][3], I[3][3];
693 +
694 +  // use the angular velocities to propagate the rotation matrix a
695 +  // full time step
696 +
697 +  sd->getA(A);
698 +  sd->getI(I);
699 +
700 +  // rotate about the x-axis
701 +  angle = dt2 * ji[0] / I[0][0];
702 +  this->rotate( 1, 2, angle, ji, A );
703 +
704 +  // rotate about the y-axis
705 +  angle = dt2 * ji[1] / I[1][1];
706 +  this->rotate( 2, 0, angle, ji, A );
707 +
708 +  // rotate about the z-axis
709 +  angle = dt * ji[2] / I[2][2];
710 +  this->rotate( 0, 1, angle, ji, A);
711 +
712 +  // rotate about the y-axis
713 +  angle = dt2 * ji[1] / I[1][1];
714 +  this->rotate( 2, 0, angle, ji, A );
715 +
716 +  // rotate about the x-axis
717 +  angle = dt2 * ji[0] / I[0][0];
718 +  this->rotate( 1, 2, angle, ji, A );
719 +
720 +  sd->setA( A  );
721 + }
722 +
723   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
724                                                  double angle, double ji[3],
725                                                  double A[3][3]){
# Line 724 | Line 785 | template<typename T> void Integrator<T>::rotate(int ax
785      }
786    }
787  
788 <  // rotate the Rotation matrix acording to:
788 >  // rotate the Rotation matrix acording to:
789    //            A[][] = A[][] * transpose(rot[][])
790  
791  
# Line 749 | Line 810 | template<typename T> void Integrator<T>::thermalize(){
810   template<typename T> void Integrator<T>::thermalize(){
811    tStats->velocitize();
812   }
813 +
814 + template<typename T> double Integrator<T>::getConservedQuantity(void){
815 +  return tStats->getTotalE();
816 + }
817 + template<typename T> string Integrator<T>::getAdditionalParameters(void){
818 +  //By default, return a null string
819 +  //The reason we use string instead of char* is that if we use char*, we will
820 +  //return a pointer point to local variable which might cause problem
821 +  return string();
822 + }

Diff Legend

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