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 763 by tim, Mon Sep 15 16:52:02 2003 UTC vs.
Revision 929 by tim, Tue Jan 13 15:46:49 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  
# Line 147 | Line 150 | template<typename T> void Integrator<T>::integrate(voi
150  
151  
152   template<typename T> void Integrator<T>::integrate(void){
150 –  int i, j;                         // loop counters
153  
154    double runTime = info->run_time;
155    double sampleTime = info->sampleTime;
# Line 160 | Line 162 | template<typename T> void Integrator<T>::integrate(voi
162    double currThermal;
163    double currStatus;
164    double currReset;
165 <  
165 >
166    int calcPot, calcStress;
165 –  int isError;
167  
168    tStats = new Thermo(info);
169    statOut = new StatWriter(info);
170    dumpOut = new DumpWriter(info);
171  
172    atoms = info->atoms;
172 –  DirectionalAtom* dAtom;
173  
174    dt = info->dt;
175    dt2 = 0.5 * dt;
176  
177 +  readyCheck();
178 +
179    // initialize the forces before the first step
180  
181    calcForce(1, 1);
182 +
183 +  if (nConstrained){
184 +    preMove();
185 +    constrainA();
186 +    calcForce(1, 1);
187 +    constrainB();
188 +  }
189    
190    if (info->setTemp){
191      thermalize();
# Line 192 | Line 201 | template<typename T> void Integrator<T>::integrate(voi
201    dumpOut->writeDump(info->getTime());
202    statOut->writeStat(info->getTime());
203  
195 –  readyCheck();
204  
205   #ifdef IS_MPI
206    strcpy(checkPointMsg, "The integrator is ready to go.");
# Line 205 | Line 213 | template<typename T> void Integrator<T>::integrate(voi
213        calcStress = 1;
214      }
215  
216 + #ifdef PROFILE
217 +    startProfile( pro1 );
218 + #endif
219 +    
220      integrateStep(calcPot, calcStress);
221  
222 + #ifdef PROFILE
223 +    endProfile( pro1 );
224 +
225 +    startProfile( pro2 );
226 + #endif // profile
227 +
228      info->incrTime(dt);
229  
230      if (info->setTemp){
# Line 222 | Line 240 | template<typename T> void Integrator<T>::integrate(voi
240      }
241  
242      if (info->getTime() >= currStatus){
243 <      statOut->writeStat(info->getTime());
244 <      calcPot = 0;
243 >      statOut->writeStat(info->getTime());
244 >      calcPot = 0;
245        calcStress = 0;
246        currStatus += statusTime;
247 <    }
247 >    }
248  
249      if (info->resetIntegrator){
250        if (info->getTime() >= currReset){
# Line 234 | Line 252 | template<typename T> void Integrator<T>::integrate(voi
252          currReset += resetTime;
253        }
254      }
255 +    
256 + #ifdef PROFILE
257 +    endProfile( pro2 );
258 + #endif //profile
259  
260   #ifdef IS_MPI
261      strcpy(checkPointMsg, "successfully took a time step.");
# Line 241 | Line 263 | template<typename T> void Integrator<T>::integrate(voi
263   #endif // is_mpi
264    }
265  
244 –  dumpOut->writeFinal(info->getTime());
245 –
266    delete dumpOut;
267    delete statOut;
268   }
# Line 250 | Line 270 | template<typename T> void Integrator<T>::integrateStep
270   template<typename T> void Integrator<T>::integrateStep(int calcPot,
271                                                         int calcStress){
272    // Position full step, and velocity half step
273 +
274 + #ifdef PROFILE
275 +  startProfile(pro3);
276 + #endif //profile
277 +
278    preMove();
279  
280 + #ifdef PROFILE
281 +  endProfile(pro3);
282 +
283 +  startProfile(pro4);
284 + #endif // profile
285 +
286    moveA();
287  
288 <  if (nConstrained){
289 <    constrainA();
290 <  }
288 > #ifdef PROFILE
289 >  endProfile(pro4);
290 >  
291 >  startProfile(pro5);
292 > #endif//profile
293  
294  
295   #ifdef IS_MPI
# Line 274 | Line 307 | template<typename T> void Integrator<T>::integrateStep
307    MPIcheckPoint();
308   #endif // is_mpi
309  
310 + #ifdef PROFILE
311 +  endProfile( pro5 );
312  
313 +  startProfile( pro6 );
314 + #endif //profile
315 +
316    // finish the velocity  half step
317  
318    moveB();
319  
320 <  if (nConstrained){
321 <    constrainB();
322 <  }
320 > #ifdef PROFILE
321 >  endProfile(pro6);
322 > #endif // profile
323  
324   #ifdef IS_MPI
325    strcpy(checkPointMsg, "Succesful moveB\n");
# Line 294 | Line 332 | template<typename T> void Integrator<T>::moveA(void){
332    int i, j;
333    DirectionalAtom* dAtom;
334    double Tb[3], ji[3];
297 –  double A[3][3], I[3][3];
298 –  double angle;
335    double vel[3], pos[3], frc[3];
336    double mass;
337  
# Line 331 | Line 367 | template<typename T> void Integrator<T>::moveA(void){
367        for (j = 0; j < 3; j++)
368          ji[j] += (dt2 * Tb[j]) * eConvert;
369  
370 <      // use the angular velocities to propagate the rotation matrix a
335 <      // full time step
370 >      this->rotationPropagation( dAtom, ji );
371  
372 <      dAtom->getA(A);
373 <      dAtom->getI(I);
372 >      dAtom->setJ(ji);
373 >    }
374 >  }
375  
376 <      // rotate about the x-axis      
377 <      angle = dt2 * ji[0] / I[0][0];
342 <      this->rotate(1, 2, angle, ji, A);
343 <
344 <      // rotate about the y-axis
345 <      angle = dt2 * ji[1] / I[1][1];
346 <      this->rotate(2, 0, angle, ji, A);
347 <
348 <      // rotate about the z-axis
349 <      angle = dt * ji[2] / I[2][2];
350 <      this->rotate(0, 1, angle, ji, A);
351 <
352 <      // rotate about the y-axis
353 <      angle = dt2 * ji[1] / I[1][1];
354 <      this->rotate(2, 0, angle, ji, A);
355 <
356 <      // rotate about the x-axis
357 <      angle = dt2 * ji[0] / I[0][0];
358 <      this->rotate(1, 2, angle, ji, A);
359 <
360 <      dAtom->setJ(ji);
361 <      dAtom->setA(A);
362 <    }
376 >  if (nConstrained){
377 >    constrainA();
378    }
379   }
380  
# Line 386 | Line 401 | template<typename T> void Integrator<T>::moveB(void){
401      if (atoms[i]->isDirectional()){
402        dAtom = (DirectionalAtom *) atoms[i];
403  
404 <      // get and convert the torque to body frame      
404 >      // get and convert the torque to body frame
405  
406        dAtom->getTrq(Tb);
407        dAtom->lab2Body(Tb);
# Line 402 | Line 417 | template<typename T> void Integrator<T>::moveB(void){
417        dAtom->setJ(ji);
418      }
419    }
420 +
421 +  if (nConstrained){
422 +    constrainB();
423 +  }
424   }
425  
426   template<typename T> void Integrator<T>::preMove(void){
# Line 420 | Line 439 | template<typename T> void Integrator<T>::constrainA(){
439   }
440  
441   template<typename T> void Integrator<T>::constrainA(){
442 <  int i, j, k;
442 >  int i, j;
443    int done;
444    double posA[3], posB[3];
445    double velA[3], velB[3];
# Line 560 | Line 579 | template<typename T> void Integrator<T>::constrainA(){
579      painCave.isFatal = 1;
580      simError();
581    }
582 +
583   }
584  
585   template<typename T> void Integrator<T>::constrainB(void){
586 <  int i, j, k;
586 >  int i, j;
587    int done;
588    double posA[3], posB[3];
589    double velA[3], velB[3];
# Line 572 | Line 592 | template<typename T> void Integrator<T>::constrainB(vo
592    int a, b, ax, ay, az, bx, by, bz;
593    double rma, rmb;
594    double dx, dy, dz;
595 <  double rabsq, pabsq, rvab;
576 <  double diffsq;
595 >  double rvab;
596    double gab;
597    int iteration;
598  
# Line 661 | Line 680 | template<typename T> void Integrator<T>::constrainB(vo
680      painCave.isFatal = 1;
681      simError();
682    }
683 + }
684 +
685 + template<typename T> void Integrator<T>::rotationPropagation
686 + ( DirectionalAtom* dAtom, double ji[3] ){
687 +
688 +  double angle;
689 +  double A[3][3], I[3][3];
690 +
691 +  // use the angular velocities to propagate the rotation matrix a
692 +  // full time step
693 +
694 +  dAtom->getA(A);
695 +  dAtom->getI(I);
696 +
697 +  // rotate about the x-axis
698 +  angle = dt2 * ji[0] / I[0][0];
699 +  this->rotate( 1, 2, angle, ji, A );
700 +
701 +  // rotate about the y-axis
702 +  angle = dt2 * ji[1] / I[1][1];
703 +  this->rotate( 2, 0, angle, ji, A );
704 +
705 +  // rotate about the z-axis
706 +  angle = dt * ji[2] / I[2][2];
707 +  this->rotate( 0, 1, angle, ji, A);
708 +
709 +  // rotate about the y-axis
710 +  angle = dt2 * ji[1] / I[1][1];
711 +  this->rotate( 2, 0, angle, ji, A );
712 +
713 +  // rotate about the x-axis
714 +  angle = dt2 * ji[0] / I[0][0];
715 +  this->rotate( 1, 2, angle, ji, A );
716 +
717 +  dAtom->setA( A  );
718   }
719  
720   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
# Line 728 | Line 782 | template<typename T> void Integrator<T>::rotate(int ax
782      }
783    }
784  
785 <  // rotate the Rotation matrix acording to:
785 >  // rotate the Rotation matrix acording to:
786    //            A[][] = A[][] * transpose(rot[][])
787  
788  
# Line 756 | Line 810 | template<typename T> double Integrator<T>::getConserve
810  
811   template<typename T> double Integrator<T>::getConservedQuantity(void){
812    return tStats->getTotalE();
813 < }
813 > }
814 > template<typename T> string Integrator<T>::getAdditionalParameters(void){
815 >  //By default, return a null string
816 >  //The reason we use string instead of char* is that if we use char*, we will
817 >  //return a pointer point to local variable which might cause problem
818 >  return string();
819 > }

Diff Legend

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