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 892 by chuckv, Mon Dec 22 21:27:04 2003 UTC vs.
Revision 1180 by chrisfen, Thu May 20 20:24:07 2004 UTC

# Line 31 | Line 31 | template<typename T> Integrator<T>::Integrator(SimInfo
31    }
32  
33    nAtoms = info->n_atoms;
34 <
34 >  integrableObjects = info->integrableObjects;
35 >
36    // check for constraints
37  
38    constrainedA = NULL;
# Line 44 | 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 68 | 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 114 | Line 119 | template<typename T> void Integrator<T>::checkConstrai
119      }
120    }
121  
122 +
123    if (nConstrained > 0){
124      isConstrained = 1;
125  
# Line 135 | 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 157 | 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;
# Line 176 | Line 182 | template<typename T> void Integrator<T>::integrate(voi
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();
# Line 207 | Line 222 | template<typename T> void Integrator<T>::integrate(voi
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      }
# Line 241 | Line 257 | template<typename T> void Integrator<T>::integrate(voi
257  
258      if (info->getTime() >= currStatus){
259        statOut->writeStat(info->getTime());
260 +      statOut->writeRaw(info->getTime());
261        calcPot = 0;
262        calcStress = 0;
263        currStatus += statusTime;
# Line 263 | 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  
267 –  // write the last frame
268 –  dumpOut->writeDump(info->getTime());
269 –
288    delete dumpOut;
289    delete statOut;
290   }
# Line 333 | 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 < nAtoms; i++){
361 <    atoms[i]->getVel(vel);
362 <    atoms[i]->getPos(pos);
363 <    atoms[i]->getFrc(frc);
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  
347 –    mass = atoms[i]->getMass();
348 –
367      for (j = 0; j < 3; j++){
368        // velocity half step
369        vel[j] += (dt2 * frc[j] / mass) * eConvert;
# Line 353 | 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()){
360 <      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 385 | Line 402 | template<typename T> void Integrator<T>::moveB(void){
402  
403   template<typename T> void Integrator<T>::moveB(void){
404    int i, j;
388 –  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()){
406 <      dAtom = (DirectionalAtom *) atoms[i];
421 >    if (integrableObjects[i]->isDirectional()){
422  
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 687 | 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);
714 >  sd->getA(A);
715 >  sd->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 );
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 <  // rotate about the y-axis
726 <  angle = dt2 * ji[1] / I[1][1];
707 <  this->rotate( 2, 0, angle, ji, A );
708 <
709 <  // rotate about the z-axis
710 <  angle = dt * ji[2] / I[2][2];
711 <  this->rotate( 0, 1, angle, ji, A);
725 >    angle = dt * ji[k] / I[k][k];
726 >    this->rotate( i, j, angle, ji, A);
727  
728 <  // rotate about the y-axis
729 <  angle = dt2 * ji[1] / I[1][1];
715 <  this->rotate( 2, 0, angle, ji, A );
728 >    angle = dt2 * ji[j] / I[j][j];
729 >    this->rotate( k, i, angle, ji, A );
730  
731 <  // rotate about the x-axis
732 <  angle = dt2 * ji[0] / I[0][0];
733 <  this->rotate( 1, 2, angle, ji, A );
734 <
735 <  dAtom->setA( A  );
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,

Diff Legend

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