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 596 by mmeineke, Fri Jul 11 22:34:48 2003 UTC vs.
Revision 597 by mmeineke, Mon Jul 14 21:28:54 2003 UTC

# Line 224 | Line 224 | void Integrator::integrate( void ){
224    pos  = Atom::getPosArray();
225    vel  = Atom::getVelArray();
226    frc  = Atom::getFrcArray();
227 –  trq  = Atom::getTrqArray();
228 –  Amat = Atom::getAmatArray();
227  
228    while( currTime < runTime ){
229  
# Line 234 | Line 232 | void Integrator::integrate( void ){
232        calcStress = 1;
233      }
234  
235 <    std::cerr << "calcPot = " << calcPot << "; calcStress = "
238 <              << calcStress << "\n";
235 >    std::cerr << currTime << "\n";
236  
237      integrateStep( calcPot, calcStress );
238        
# Line 282 | Line 279 | void Integrator::integrateStep( int calcPot, int calcS
279  
280    preMove();
281    moveA();
282 <  if( nConstrained ) constrainA();
282 >  //if( nConstrained ) constrainA();
283  
284    // calc forces
285  
# Line 304 | Line 301 | void Integrator::moveA( void ){
301    double Tb[3];
302    double ji[3];
303    double angle;
304 <  double A[3][3];
304 >  double A[3][3], At[3][3];
305  
306  
307    for( i=0; i<nAtoms; i++ ){
# Line 315 | Line 312 | void Integrator::moveA( void ){
312      for( j=atomIndex; j<(atomIndex+3); j++ )
313        vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
314  
318 –    std::cerr<< "MoveA vel[" << i << "] = "
319 –             << vel[atomIndex] << "\t"
320 –             << vel[atomIndex+1]<< "\t"
321 –             << vel[atomIndex+2]<< "\n";
315  
316      // position whole step    
317      for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j];
318      
319  
327 –    std::cerr<< "MoveA pos[" << i << "] = "
328 –             << pos[atomIndex] << "\t"
329 –             << pos[atomIndex+1]<< "\t"
330 –             << pos[atomIndex+2]<< "\n";
331 –
320      if( atoms[i]->isDirectional() ){
321  
322        dAtom = (DirectionalAtom *)atoms[i];
# Line 338 | Line 326 | void Integrator::moveA( void ){
326        Tb[0] = dAtom->getTx();
327        Tb[1] = dAtom->getTy();
328        Tb[2] = dAtom->getTz();
329 <      
329 >
330        dAtom->lab2Body( Tb );
331 <      
331 >
332        // get the angular momentum, and propagate a half step
333        
334        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert;
# Line 349 | Line 337 | void Integrator::moveA( void ){
337        
338        // use the angular velocities to propagate the rotation matrix a
339        // full time step
352 –      
353 –          // get the atom's rotation matrix
354 –          
355 –      A[0][0] = dAtom->getAxx();
356 –      A[0][1] = dAtom->getAxy();
357 –      A[0][2] = dAtom->getAxz();
340        
359 –      A[1][0] = dAtom->getAyx();
360 –      A[1][1] = dAtom->getAyy();
361 –      A[1][2] = dAtom->getAyz();
362 –      
363 –      A[2][0] = dAtom->getAzx();
364 –      A[2][1] = dAtom->getAzy();
365 –      A[2][2] = dAtom->getAzz();
366 –      
341        // rotate about the x-axis      
342        angle = dt2 * ji[0] / dAtom->getIxx();
343 <      this->rotate( 1, 2, angle, ji, A );
344 <      
343 >      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
344 >
345        // rotate about the y-axis
346        angle = dt2 * ji[1] / dAtom->getIyy();
347 <      this->rotate( 2, 0, angle, ji, A );
347 >      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
348        
349        // rotate about the z-axis
350        angle = dt * ji[2] / dAtom->getIzz();
351 <      this->rotate( 0, 1, angle, ji, A );
351 >      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
352        
353        // rotate about the y-axis
354        angle = dt2 * ji[1] / dAtom->getIyy();
355 <      this->rotate( 2, 0, angle, ji, A );
355 >      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
356        
357         // rotate about the x-axis
358        angle = dt2 * ji[0] / dAtom->getIxx();
359 <      this->rotate( 1, 2, angle, ji, A );
359 >      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
360        
361        dAtom->setJx( ji[0] );
362        dAtom->setJy( ji[1] );
363        dAtom->setJz( ji[2] );
364 +
365 +      std::cerr << "Amat[" << i << "]\n";
366 +      info->printMat9( &Amat[aMatIndex] );
367 +          
368 +      std::cerr << "ji[" << i << "]\t"
369 +                << ji[0] << "\t"
370 +                << ji[1] << "\t"
371 +                << ji[2] << "\n";
372 +          
373      }
374      
375    }
# Line 395 | Line 378 | void Integrator::moveB( void ){
378  
379   void Integrator::moveB( void ){
380    int i,j,k;
381 <  int atomIndex;
381 >  int atomIndex, aMatIndex;
382    DirectionalAtom* dAtom;
383    double Tb[3];
384    double ji[3];
385  
386    for( i=0; i<nAtoms; i++ ){
387      atomIndex = i * 3;
388 +    aMatIndex = i * 9;
389  
390      // velocity half step
391      for( j=atomIndex; j<(atomIndex+3); j++ )
392        vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert;
393  
394 <    std::cerr<< "MoveB vel[" << i << "] = "
411 <             << vel[atomIndex] << "\t"
412 <             << vel[atomIndex+1]<< "\t"
413 <             << vel[atomIndex+2]<< "\n";
414 <
415 <
394 >
395      if( atoms[i]->isDirectional() ){
396        
397        dAtom = (DirectionalAtom *)atoms[i];
# Line 423 | Line 402 | void Integrator::moveB( void ){
402        Tb[1] = dAtom->getTy();
403        Tb[2] = dAtom->getTz();
404        
405 +      std::cerr << "TrqB[" << i << "]\t"
406 +                << Tb[0] << "\t"
407 +                << Tb[1] << "\t"
408 +                << Tb[2] << "\n";
409 +
410        dAtom->lab2Body( Tb );
411        
412        // get the angular momentum, and complete the angular momentum
# Line 435 | Line 419 | void Integrator::moveB( void ){
419        dAtom->setJx( ji[0] );
420        dAtom->setJy( ji[1] );
421        dAtom->setJz( ji[2] );
422 +
423 +
424 +      std::cerr << "Amat[" << i << "]\n";
425 +      info->printMat9( &Amat[aMatIndex] );
426 +          
427 +      std::cerr << "ji[" << i << "]\t"
428 +                << ji[0] << "\t"
429 +                << ji[1] << "\t"
430 +                << ji[2] << "\n";
431      }
432    }
433  
# Line 690 | Line 683 | void Integrator::rotate( int axes1, int axes2, double
683  
684  
685   void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],
686 <                         double A[3][3] ){
686 >                         double A[9] ){
687  
688    int i,j,k;
689    double sinAngle;
# Line 702 | Line 695 | void Integrator::rotate( int axes1, int axes2, double
695    double tempA[3][3];
696    double tempJ[3];
697  
698 +
699    // initialize the tempA
700  
701    for(i=0; i<3; i++){
702      for(j=0; j<3; j++){
703 <      tempA[j][i] = A[i][j];
703 >      tempA[j][i] = A[3*i+j];
704      }
705    }
706  
# Line 763 | Line 757 | void Integrator::rotate( int axes1, int axes2, double
757  
758    for(i=0; i<3; i++){
759      for(j=0; j<3; j++){
760 <      A[j][i] = 0.0;
760 >      A[3*j+i] = 0.0;
761        for(k=0; k<3; k++){
762 <        A[j][i] += tempA[i][k] * rot[j][k];
762 >        A[3*j+i] += tempA[i][k] * rot[j][k];
763        }
764      }
765    }

Diff Legend

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