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 768 by mmeineke, Wed Sep 17 14:22:15 2003 UTC vs.
Revision 781 by tim, Mon Sep 22 23:07:57 2003 UTC

# Line 177 | Line 177 | template<typename T> void Integrator<T>::integrate(voi
177    // initialize the forces before the first step
178  
179    calcForce(1, 1);
180 +
181 +  if (nConstrained){
182 +    preMove();
183 +    constrainA();
184 +    calcForce(1, 1);    
185 +    constrainB();
186 +  }
187    
188    if (info->setTemp){
189      thermalize();
# Line 290 | Line 297 | template<typename T> void Integrator<T>::moveA(void){
297    int i, j;
298    DirectionalAtom* dAtom;
299    double Tb[3], ji[3];
293 –  double A[3][3], I[3][3];
294 –  double angle;
300    double vel[3], pos[3], frc[3];
301    double mass;
302  
# Line 327 | Line 332 | template<typename T> void Integrator<T>::moveA(void){
332        for (j = 0; j < 3; j++)
333          ji[j] += (dt2 * Tb[j]) * eConvert;
334  
335 <      // use the angular velocities to propagate the rotation matrix a
331 <      // full time step
332 <
333 <      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);
335 >      this->rotationPropagation( dAtom, ji );
336  
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 –
337        dAtom->setJ(ji);
357 –      dAtom->setA(A);
338      }
339    }
340  
# Line 668 | Line 648 | template<typename T> void Integrator<T>::constrainB(vo
648    }
649   }
650  
651 + template<typename T> void Integrator<T>::rotationPropagation
652 + ( DirectionalAtom* dAtom, double ji[3] ){
653 +
654 +  double angle;
655 +  double A[3][3], I[3][3];
656 +
657 +  // use the angular velocities to propagate the rotation matrix a
658 +  // full time step
659 +
660 +  dAtom->getA(A);
661 +  dAtom->getI(I);
662 +  
663 +  // rotate about the x-axis      
664 +  angle = dt2 * ji[0] / I[0][0];
665 +  this->rotate( 1, 2, angle, ji, A );
666 +  
667 +  // rotate about the y-axis
668 +  angle = dt2 * ji[1] / I[1][1];
669 +  this->rotate( 2, 0, angle, ji, A );
670 +  
671 +  // rotate about the z-axis
672 +  angle = dt * ji[2] / I[2][2];
673 +  this->rotate( 0, 1, angle, ji, A);
674 +  
675 +  // rotate about the y-axis
676 +  angle = dt2 * ji[1] / I[1][1];
677 +  this->rotate( 2, 0, angle, ji, A );
678 +  
679 +  // rotate about the x-axis
680 +  angle = dt2 * ji[0] / I[0][0];
681 +  this->rotate( 1, 2, angle, ji, A );
682 +  
683 +  dAtom->setA( A  );    
684 + }
685 +
686   template<typename T> void Integrator<T>::rotate(int axes1, int axes2,
687                                                  double angle, double ji[3],
688                                                  double A[3][3]){

Diff Legend

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