ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Symplectic.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Symplectic.cpp (file contents):
Revision 475 by gezelter, Tue Apr 8 12:44:18 2003 UTC vs.
Revision 482 by chuckv, Tue Apr 8 22:38:43 2003 UTC

# Line 164 | Line 164 | void Symplectic::integrate( void ){
164    double dt2;                       // half the dt
165  
166    double vx, vy, vz;    // the velocities
167 < //  double vx2, vy2, vz2; // the square of the velocities
167 >  double vx2, vy2, vz2; // the square of the velocities
168    double rx, ry, rz;    // the postitions
169    
170    double ji[3];   // the body frame angular momentum
# Line 243 | Line 243 | void Symplectic::integrate( void ){
243      for( tl=0; tl < n_loops; tl++ ){
244  
245        if (!strcasecmp( entry_plug->ensemble, "NVT"))
246 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
246 >        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
247        
248        for( j=0; j<nAtoms; j++ ){
249  
# Line 281 | Line 281 | void Symplectic::integrate( void ){
281        }
282  
283  
284 < //       for( i=0; i<nAtoms; i++ ){
285 < // //   if( atoms[i]->isDirectional() ){
284 >      for( i=0; i<nAtoms; i++ ){
285 >        if( atoms[i]->isDirectional() ){
286                    
287 < // //     dAtom = (DirectionalAtom *)atoms[i];
287 >          dAtom = (DirectionalAtom *)atoms[i];
288            
289 < // //     // get and convert the torque to body frame
290 <          
291 < // //     Tb[0] = dAtom->getTx();
292 < // //     Tb[1] = dAtom->getTy();
293 < // //     Tb[2] = dAtom->getTz();
289 >          // get and convert the torque to body frame
290            
291 < // //     dAtom->lab2Body( Tb );
291 >          Tb[0] = dAtom->getTx();
292 >          Tb[1] = dAtom->getTy();
293 >          Tb[2] = dAtom->getTz();
294            
295 < // //     // get the angular momentum, and propagate a half step
295 >          dAtom->lab2Body( Tb );
296            
297 < // //     ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
300 < // //     ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
301 < // //     ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
297 >          // get the angular momentum, and propagate a half step
298            
299 < // //     // get the atom's rotation matrix
299 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
300 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
301 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
302            
303 < // //     A[0][0] = dAtom->getAxx();
306 < // //     A[0][1] = dAtom->getAxy();
307 < // //     A[0][2] = dAtom->getAxz();
303 >          // get the atom's rotation matrix
304            
305 < // //     A[1][0] = dAtom->getAyx();
306 < // //     A[1][1] = dAtom->getAyy();
307 < // //     A[1][2] = dAtom->getAyz();
305 >          A[0][0] = dAtom->getAxx();
306 >          A[0][1] = dAtom->getAxy();
307 >          A[0][2] = dAtom->getAxz();
308            
309 < // //     A[2][0] = dAtom->getAzx();
310 < // //     A[2][1] = dAtom->getAzy();
311 < // //     A[2][2] = dAtom->getAzz();
309 >          A[1][0] = dAtom->getAyx();
310 >          A[1][1] = dAtom->getAyy();
311 >          A[1][2] = dAtom->getAyz();
312            
313 <          
314 < // //     // use the angular velocities to propagate the rotation matrix a
315 < // //     // full time step
313 >          A[2][0] = dAtom->getAzx();
314 >          A[2][1] = dAtom->getAzy();
315 >          A[2][2] = dAtom->getAzz();
316            
317            
318 < // //     angle = dt2 * ji[0] / dAtom->getIxx();
319 < // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
318 >          // use the angular velocities to propagate the rotation matrix a
319 >          // full time step
320            
325 – // //     angle = dt2 * ji[1] / dAtom->getIyy();
326 – // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
321            
322 < // //     angle = dt * ji[2] / dAtom->getIzz();
323 < // //     this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
322 >          angle = dt2 * ji[0] / dAtom->getIxx();
323 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
324            
325 < // //     angle = dt2 * ji[1] / dAtom->getIyy();
326 < // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
325 >          angle = dt2 * ji[1] / dAtom->getIyy();
326 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
327            
328 < // //     angle = dt2 * ji[0] / dAtom->getIxx();
329 < // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
328 >          angle = dt * ji[2] / dAtom->getIzz();
329 >          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
330            
331 +          angle = dt2 * ji[1] / dAtom->getIyy();
332 +          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
333            
334 < // //     dAtom->setA( A );
335 < // //     dAtom->setJx( ji[0] );
336 < // //     dAtom->setJy( ji[1] );
337 < // //     dAtom->setJz( ji[2] );
338 < // //   }
339 < //       }
334 >          angle = dt2 * ji[0] / dAtom->getIxx();
335 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
336 >          
337 >          
338 >          dAtom->setA( A );
339 >          dAtom->setJx( ji[0] );
340 >          dAtom->setJy( ji[1] );
341 >          dAtom->setJz( ji[2] );
342 >        }
343 >      }
344        
345        // calculate the forces
346        
# Line 382 | Line 382 | void Symplectic::integrate( void ){
382          atoms[j]->set_vz(Vz[j]);
383        }
384        
385 < //       for( i=0; i< nAtoms; i++ ){
385 >      for( i=0; i< nAtoms; i++ ){
386 >
387 >        if( atoms[i]->isDirectional() ){
388  
389 < //      if( atoms[i]->isDirectional() ){
388 <
389 < //        dAtom = (DirectionalAtom *)atoms[i];
389 >          dAtom = (DirectionalAtom *)atoms[i];
390            
391 < //        // get and convert the torque to body frame
391 >          // get and convert the torque to body frame
392            
393 < //        Tb[0] = dAtom->getTx();
394 < //        Tb[1] = dAtom->getTy();
395 < //        Tb[2] = dAtom->getTz();
393 >          Tb[0] = dAtom->getTx();
394 >          Tb[1] = dAtom->getTy();
395 >          Tb[2] = dAtom->getTz();
396            
397 < //        dAtom->lab2Body( Tb );
397 >          dAtom->lab2Body( Tb );
398            
399 < //        // get the angular momentum, and complete the angular momentum
400 < //        // half step
399 >          // get the angular momentum, and complete the angular momentum
400 >          // half step
401            
402 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
403 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
404 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
402 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
403 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
404 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
405            
406 < //        dAtom->setJx( ji[0] );
407 < //        dAtom->setJy( ji[1] );
408 < //        dAtom->setJz( ji[2] );
409 < //      }
410 < //       }
406 >          dAtom->setJx( ji[0] );
407 >          dAtom->setJy( ji[1] );
408 >          dAtom->setJz( ji[2] );
409 >        }
410 >      }
411      
412  
413        if (!strcasecmp( entry_plug->ensemble, "NVT"))
414 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
414 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
415  
416        if (!strcasecmp( entry_plug->ensemble, "NPT") )
417          myES->NoseHooverAndersonNPT( dt,
# Line 446 | Line 446 | void Symplectic::integrate( void ){
446        trans_kE = 0.0;
447  
448        if (!strcasecmp( entry_plug->ensemble, "NVT"))
449 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
449 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
450        
451        for( i=0; i<nAtoms; i++ ){
452          
# Line 473 | Line 473 | void Symplectic::integrate( void ){
473          atoms[i]->set_vy( vy );
474          atoms[i]->set_vz( vz );
475          
476 < //      if( atoms[i]->isDirectional() ){
476 >        if( atoms[i]->isDirectional() ){
477  
478 < //        dAtom = (DirectionalAtom *)atoms[i];
478 >          dAtom = (DirectionalAtom *)atoms[i];
479            
480 < //        // get and convert the torque to body frame
480 >          // get and convert the torque to body frame
481            
482 < //        Tb[0] = dAtom->getTx();
483 < //        Tb[1] = dAtom->getTy();
484 < //        Tb[2] = dAtom->getTz();
482 >          Tb[0] = dAtom->getTx();
483 >          Tb[1] = dAtom->getTy();
484 >          Tb[2] = dAtom->getTz();
485            
486 < //        dAtom->lab2Body( Tb );
486 >          dAtom->lab2Body( Tb );
487            
488 < //        // get the angular momentum, and propagate a half step
488 >          // get the angular momentum, and propagate a half step
489            
490 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
491 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
492 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
493 <          
494 < //        // get the atom's rotation matrix
490 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
491 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
492 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
493            
494 < //        A[0][0] = dAtom->getAxx();
497 < //        A[0][1] = dAtom->getAxy();
498 < //        A[0][2] = dAtom->getAxz();
494 >          // get the atom's rotation matrix
495            
496 < //        A[1][0] = dAtom->getAyx();
497 < //        A[1][1] = dAtom->getAyy();
498 < //        A[1][2] = dAtom->getAyz();
496 >          A[0][0] = dAtom->getAxx();
497 >          A[0][1] = dAtom->getAxy();
498 >          A[0][2] = dAtom->getAxz();
499            
500 < //        A[2][0] = dAtom->getAzx();
501 < //        A[2][1] = dAtom->getAzy();
502 < //        A[2][2] = dAtom->getAzz();
500 >          A[1][0] = dAtom->getAyx();
501 >          A[1][1] = dAtom->getAyy();
502 >          A[1][2] = dAtom->getAyz();
503            
504 +          A[2][0] = dAtom->getAzx();
505 +          A[2][1] = dAtom->getAzy();
506 +          A[2][2] = dAtom->getAzz();
507            
509 – //        // use the angular velocities to propagate the rotation matrix a
510 – //        // full time step
508            
509 +          // use the angular velocities to propagate the rotation matrix a
510 +          // full time step
511            
513 – //        angle = dt2 * ji[0] / dAtom->getIxx();
514 – //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
512            
513 < //        angle = dt2 * ji[1] / dAtom->getIyy();
514 < //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
513 >          angle = dt2 * ji[0] / dAtom->getIxx();
514 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
515            
516 < //        angle = dt * ji[2] / dAtom->getIzz();
517 < //        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
516 >          angle = dt2 * ji[1] / dAtom->getIyy();
517 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
518            
519 < //        angle = dt2 * ji[1] / dAtom->getIyy();
520 < //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
519 >          angle = dt * ji[2] / dAtom->getIzz();
520 >          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
521            
522 < //        angle = dt2 * ji[0] / dAtom->getIxx();
523 < //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
522 >          angle = dt2 * ji[1] / dAtom->getIyy();
523 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
524            
525 +          angle = dt2 * ji[0] / dAtom->getIxx();
526 +          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
527            
528 < //        dAtom->setA( A );
529 < //        dAtom->setJx( ji[0] );
530 < //        dAtom->setJy( ji[1] );
531 < //        dAtom->setJz( ji[2] );
532 < //      }
528 >          
529 >          dAtom->setA( A );
530 >          dAtom->setJx( ji[0] );
531 >          dAtom->setJy( ji[1] );
532 >          dAtom->setJz( ji[2] );
533 >        }
534        }
535        
536        // calculate the forces
# Line 552 | Line 552 | void Symplectic::integrate( void ){
552          atoms[i]->set_vy( vy );
553          atoms[i]->set_vz( vz );
554          
555 < //      vx2 = vx * vx;
556 < //      vy2 = vy * vy;
557 < //      vz2 = vz * vz;
555 >        vx2 = vx * vx;
556 >        vy2 = vy * vy;
557 >        vz2 = vz * vz;
558          
559 < //      if( atoms[i]->isDirectional() ){
559 >        if( atoms[i]->isDirectional() ){
560  
561 < //        dAtom = (DirectionalAtom *)atoms[i];
561 >          dAtom = (DirectionalAtom *)atoms[i];
562            
563 < //        // get and convert the torque to body frame
563 >          // get and convert the torque to body frame
564            
565 < //        Tb[0] = dAtom->getTx();
566 < //        Tb[1] = dAtom->getTy();
567 < //        Tb[2] = dAtom->getTz();
565 >          Tb[0] = dAtom->getTx();
566 >          Tb[1] = dAtom->getTy();
567 >          Tb[2] = dAtom->getTz();
568            
569 < //        dAtom->lab2Body( Tb );
569 >          dAtom->lab2Body( Tb );
570            
571 < //        // get the angular momentum, and complete the angular momentum
572 < //        // half step
571 >          // get the angular momentum, and complete the angular momentum
572 >          // half step
573            
574 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
575 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
576 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
574 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
575 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
576 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
577            
578 < //        jx2 = ji[0] * ji[0];
579 < //        jy2 = ji[1] * ji[1];
580 < //        jz2 = ji[2] * ji[2];
578 >          jx2 = ji[0] * ji[0];
579 >          jy2 = ji[1] * ji[1];
580 >          jz2 = ji[2] * ji[2];
581            
582 < //        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
583 < //          + (jz2 / dAtom->getIzz());
582 >          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
583 >            + (jz2 / dAtom->getIzz());
584            
585 < //        dAtom->setJx( ji[0] );
586 < //        dAtom->setJy( ji[1] );
587 < //        dAtom->setJz( ji[2] );
588 < //      }
585 >          dAtom->setJx( ji[0] );
586 >          dAtom->setJy( ji[1] );
587 >          dAtom->setJz( ji[2] );
588 >        }
589        }
590      
591        if (!strcasecmp( entry_plug->ensemble, "NVT"))
592 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
592 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
593  
594        if (!strcasecmp( entry_plug->ensemble, "NPT") )
595          myES->NoseHooverAndersonNPT( dt,

Diff Legend

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