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 1057 by tim, Tue Feb 17 19:23:44 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 65 | Line 68 | template<typename T> void Integrator<T>::checkConstrai
68  
69    SRI** theArray;
70    for (int i = 0; i < nMols; i++){
71 <    theArray = (SRI * *) molecules[i].getMyBonds();
71 >
72 >          theArray = (SRI * *) molecules[i].getMyBonds();
73      for (int j = 0; j < molecules[i].getNBonds(); j++){
74        constrained = theArray[j]->is_constrained();
75  
# Line 111 | Line 115 | template<typename T> void Integrator<T>::checkConstrai
115      }
116    }
117  
118 +
119    if (nConstrained > 0){
120      isConstrained = 1;
121  
# Line 132 | Line 137 | template<typename T> void Integrator<T>::checkConstrai
137      }
138  
139  
140 <    // save oldAtoms to check for lode balanceing later on.
140 >    // save oldAtoms to check for lode balancing later on.
141  
142      oldAtoms = nAtoms;
143  
# Line 147 | Line 152 | template<typename T> void Integrator<T>::integrate(voi
152  
153  
154   template<typename T> void Integrator<T>::integrate(void){
150 –  int i, j;                         // loop counters
155  
156    double runTime = info->run_time;
157    double sampleTime = info->sampleTime;
158    double statusTime = info->statusTime;
159    double thermalTime = info->thermalTime;
160 +  double resetTime = info->resetTime;
161  
162 +
163    double currSample;
164    double currThermal;
165    double currStatus;
166 +  double currReset;
167  
168    int calcPot, calcStress;
162 –  int isError;
169  
170    tStats = new Thermo(info);
171    statOut = new StatWriter(info);
172    dumpOut = new DumpWriter(info);
173  
174    atoms = info->atoms;
169 –  DirectionalAtom* dAtom;
175  
176    dt = info->dt;
177    dt2 = 0.5 * dt;
178  
179 +  readyCheck();
180 +
181    // initialize the forces before the first step
182  
183    calcForce(1, 1);
177 –  // myFF->doForces(1,1);
184  
185 +  //temp test
186 +  tStats->getPotential();
187 +  
188 +  if (nConstrained){
189 +    preMove();
190 +    constrainA();
191 +    calcForce(1, 1);
192 +    constrainB();
193 +  }
194 +  
195    if (info->setTemp){
196      thermalize();
197    }
198  
183 –  calcPot = 0;
184 –  calcStress = 0;
185 –  currSample = sampleTime;
186 –  currThermal = thermalTime;
187 –  currStatus = statusTime;
188 –  
199    calcPot     = 0;
200    calcStress  = 0;
201    currSample  = sampleTime + info->getTime();
202    currThermal = thermalTime+ info->getTime();
203    currStatus  = statusTime + info->getTime();
204 +  currReset   = resetTime  + info->getTime();
205  
206    dumpOut->writeDump(info->getTime());
207    statOut->writeStat(info->getTime());
208  
198 –  readyCheck();
209  
210   #ifdef IS_MPI
211    strcpy(checkPointMsg, "The integrator is ready to go.");
# Line 208 | Line 218 | template<typename T> void Integrator<T>::integrate(voi
218        calcStress = 1;
219      }
220  
221 + #ifdef PROFILE
222 +    startProfile( pro1 );
223 + #endif
224 +    
225      integrateStep(calcPot, calcStress);
226  
227 + #ifdef PROFILE
228 +    endProfile( pro1 );
229 +
230 +    startProfile( pro2 );
231 + #endif // profile
232 +
233      info->incrTime(dt);
234  
235      if (info->setTemp){
# Line 225 | Line 245 | template<typename T> void Integrator<T>::integrate(voi
245      }
246  
247      if (info->getTime() >= currStatus){
248 <      statOut->writeStat(info->getTime());
249 <      calcPot = 0;
248 >      statOut->writeStat(info->getTime());
249 >      calcPot = 0;
250        calcStress = 0;
251        currStatus += statusTime;
252 <    }
252 >    }
253  
254 +    if (info->resetIntegrator){
255 +      if (info->getTime() >= currReset){
256 +        this->resetIntegrator();
257 +        currReset += resetTime;
258 +      }
259 +    }
260 +    
261 + #ifdef PROFILE
262 +    endProfile( pro2 );
263 + #endif //profile
264 +
265   #ifdef IS_MPI
266      strcpy(checkPointMsg, "successfully took a time step.");
267      MPIcheckPoint();
268   #endif // is_mpi
269    }
270  
240 –  dumpOut->writeFinal(info->getTime());
241 –
271    delete dumpOut;
272    delete statOut;
273   }
# Line 246 | Line 275 | template<typename T> void Integrator<T>::integrateStep
275   template<typename T> void Integrator<T>::integrateStep(int calcPot,
276                                                         int calcStress){
277    // Position full step, and velocity half step
278 +
279 + #ifdef PROFILE
280 +  startProfile(pro3);
281 + #endif //profile
282 +
283    preMove();
284  
285 + #ifdef PROFILE
286 +  endProfile(pro3);
287 +
288 +  startProfile(pro4);
289 + #endif // profile
290 +
291    moveA();
292  
293 <  if (nConstrained){
294 <    constrainA();
295 <  }
293 > #ifdef PROFILE
294 >  endProfile(pro4);
295 >  
296 >  startProfile(pro5);
297 > #endif//profile
298  
299  
300   #ifdef IS_MPI
# Line 270 | Line 312 | template<typename T> void Integrator<T>::integrateStep
312    MPIcheckPoint();
313   #endif // is_mpi
314  
315 + #ifdef PROFILE
316 +  endProfile( pro5 );
317  
318 +  startProfile( pro6 );
319 + #endif //profile
320 +
321    // finish the velocity  half step
322  
323    moveB();
324  
325 <  if (nConstrained){
326 <    constrainB();
327 <  }
325 > #ifdef PROFILE
326 >  endProfile(pro6);
327 > #endif // profile
328  
329   #ifdef IS_MPI
330    strcpy(checkPointMsg, "Succesful moveB\n");
# Line 290 | Line 337 | template<typename T> void Integrator<T>::moveA(void){
337    int i, j;
338    DirectionalAtom* dAtom;
339    double Tb[3], ji[3];
293 –  double A[3][3], I[3][3];
294 –  double angle;
340    double vel[3], pos[3], frc[3];
341    double mass;
342  
# Line 327 | Line 372 | template<typename T> void Integrator<T>::moveA(void){
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
331 <      // full time step
375 >      this->rotationPropagation( dAtom, ji );
376  
377 <      dAtom->getA(A);
378 <      dAtom->getI(I);
379 <
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);
377 >      dAtom->setJ(ji);
378 >    }
379 >  }
380  
381 <      // rotate about the x-axis
382 <      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);
359 <    }
381 >  if (nConstrained){
382 >    constrainA();
383    }
384   }
385  
# Line 383 | Line 406 | template<typename T> void Integrator<T>::moveB(void){
406      if (atoms[i]->isDirectional()){
407        dAtom = (DirectionalAtom *) atoms[i];
408  
409 <      // get and convert the torque to body frame      
409 >      // get and convert the torque to body frame
410  
411        dAtom->getTrq(Tb);
412        dAtom->lab2Body(Tb);
# Line 399 | Line 422 | template<typename T> void Integrator<T>::moveB(void){
422        dAtom->setJ(ji);
423      }
424    }
425 +
426 +  if (nConstrained){
427 +    constrainB();
428 +  }
429   }
430  
431   template<typename T> void Integrator<T>::preMove(void){
# Line 417 | Line 444 | template<typename T> void Integrator<T>::constrainA(){
444   }
445  
446   template<typename T> void Integrator<T>::constrainA(){
447 <  int i, j, k;
447 >  int i, j;
448    int done;
449    double posA[3], posB[3];
450    double velA[3], velB[3];
# Line 557 | Line 584 | template<typename T> void Integrator<T>::constrainA(){
584      painCave.isFatal = 1;
585      simError();
586    }
587 +
588   }
589  
590   template<typename T> void Integrator<T>::constrainB(void){
591 <  int i, j, k;
591 >  int i, j;
592    int done;
593    double posA[3], posB[3];
594    double velA[3], velB[3];
# Line 569 | Line 597 | template<typename T> void Integrator<T>::constrainB(vo
597    int a, b, ax, ay, az, bx, by, bz;
598    double rma, rmb;
599    double dx, dy, dz;
600 <  double rabsq, pabsq, rvab;
573 <  double diffsq;
600 >  double rvab;
601    double gab;
602    int iteration;
603  
# Line 660 | Line 687 | template<typename T> void Integrator<T>::rotate(int ax
687    }
688   }
689  
690 + template<typename T> void Integrator<T>::rotationPropagation
691 + ( DirectionalAtom* dAtom, double ji[3] ){
692 +
693 +  double angle;
694 +  double A[3][3], I[3][3];
695 +
696 +  // use the angular velocities to propagate the rotation matrix a
697 +  // full time step
698 +
699 +  dAtom->getA(A);
700 +  dAtom->getI(I);
701 +
702 +  // rotate about the x-axis
703 +  angle = dt2 * ji[0] / I[0][0];
704 +  this->rotate( 1, 2, angle, ji, A );
705 +
706 +  // rotate about the y-axis
707 +  angle = dt2 * ji[1] / I[1][1];
708 +  this->rotate( 2, 0, angle, ji, A );
709 +
710 +  // rotate about the z-axis
711 +  angle = dt * ji[2] / I[2][2];
712 +  this->rotate( 0, 1, angle, ji, A);
713 +
714 +  // rotate about the y-axis
715 +  angle = dt2 * ji[1] / I[1][1];
716 +  this->rotate( 2, 0, angle, ji, A );
717 +
718 +  // rotate about the x-axis
719 +  angle = dt2 * ji[0] / I[0][0];
720 +  this->rotate( 1, 2, angle, ji, A );
721 +
722 +  dAtom->setA( A  );
723 + }
724 +
725   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
726                                                  double angle, double ji[3],
727                                                  double A[3][3]){
# Line 725 | Line 787 | template<typename T> void Integrator<T>::rotate(int ax
787      }
788    }
789  
790 <  // rotate the Rotation matrix acording to:
790 >  // rotate the Rotation matrix acording to:
791    //            A[][] = A[][] * transpose(rot[][])
792  
793  
# Line 750 | Line 812 | template<typename T> void Integrator<T>::thermalize(){
812   template<typename T> void Integrator<T>::thermalize(){
813    tStats->velocitize();
814   }
815 +
816 + template<typename T> double Integrator<T>::getConservedQuantity(void){
817 +  return tStats->getTotalE();
818 + }
819 + template<typename T> string Integrator<T>::getAdditionalParameters(void){
820 +  //By default, return a null string
821 +  //The reason we use string instead of char* is that if we use char*, we will
822 +  //return a pointer point to local variable which might cause problem
823 +  return string();
824 + }

Diff Legend

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