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 469 by mmeineke, Mon Apr 7 20:06:31 2003 UTC vs.
Revision 472 by mmeineke, Mon Apr 7 21:16:35 2003 UTC

# Line 272 | Line 272 | void Symplectic::integrate( void ){
272        }
273  
274  
275 <      for( i=0; i<nAtoms; i++ ){
276 <        if( atoms[i]->isDirectional() ){
275 > //       for( i=0; i<nAtoms; i++ ){
276 > // //   if( atoms[i]->isDirectional() ){
277                    
278 <          dAtom = (DirectionalAtom *)atoms[i];
279 <          
280 <          // get and convert the torque to body frame
278 > // //     dAtom = (DirectionalAtom *)atoms[i];
279            
280 <          Tb[0] = dAtom->getTx();
283 <          Tb[1] = dAtom->getTy();
284 <          Tb[2] = dAtom->getTz();
280 > // //     // get and convert the torque to body frame
281            
282 <          dAtom->lab2Body( Tb );
283 <          
284 <          // get the angular momentum, and propagate a half step
282 > // //     Tb[0] = dAtom->getTx();
283 > // //     Tb[1] = dAtom->getTy();
284 > // //     Tb[2] = dAtom->getTz();
285            
286 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
291 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
292 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
286 > // //     dAtom->lab2Body( Tb );
287            
288 <          // get the atom's rotation matrix
288 > // //     // get the angular momentum, and propagate a half step
289 >          
290 > // //     ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
291 > // //     ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
292 > // //     ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
293 >          
294 > // //     // get the atom's rotation matrix
295            
296 <          A[0][0] = dAtom->getAxx();
297 <          A[0][1] = dAtom->getAxy();
298 <          A[0][2] = dAtom->getAxz();
296 > // //     A[0][0] = dAtom->getAxx();
297 > // //     A[0][1] = dAtom->getAxy();
298 > // //     A[0][2] = dAtom->getAxz();
299            
300 <          A[1][0] = dAtom->getAyx();
301 <          A[1][1] = dAtom->getAyy();
302 <          A[1][2] = dAtom->getAyz();
300 > // //     A[1][0] = dAtom->getAyx();
301 > // //     A[1][1] = dAtom->getAyy();
302 > // //     A[1][2] = dAtom->getAyz();
303            
304 <          A[2][0] = dAtom->getAzx();
305 <          A[2][1] = dAtom->getAzy();
306 <          A[2][2] = dAtom->getAzz();
307 <          
308 <          
309 <          // use the angular velocities to propagate the rotation matrix a
310 <          // full time step
304 > // //     A[2][0] = dAtom->getAzx();
305 > // //     A[2][1] = dAtom->getAzy();
306 > // //     A[2][2] = dAtom->getAzz();
307            
308            
309 <          angle = dt2 * ji[0] / dAtom->getIxx();
310 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
309 > // //     // use the angular velocities to propagate the rotation matrix a
310 > // //     // full time step
311            
316 –          angle = dt2 * ji[1] / dAtom->getIyy();
317 –          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
312            
313 <          angle = dt * ji[2] / dAtom->getIzz();
314 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
313 > // //     angle = dt2 * ji[0] / dAtom->getIxx();
314 > // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
315            
316 <          angle = dt2 * ji[1] / dAtom->getIyy();
317 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
316 > // //     angle = dt2 * ji[1] / dAtom->getIyy();
317 > // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
318            
319 <          angle = dt2 * ji[0] / dAtom->getIxx();
320 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
319 > // //     angle = dt * ji[2] / dAtom->getIzz();
320 > // //     this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
321            
322 + // //     angle = dt2 * ji[1] / dAtom->getIyy();
323 + // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
324            
325 <          dAtom->setA( A );
326 <          dAtom->setJx( ji[0] );
327 <          dAtom->setJy( ji[1] );
328 <          dAtom->setJz( ji[2] );
329 <        }
330 <      }
325 > // //     angle = dt2 * ji[0] / dAtom->getIxx();
326 > // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
327 >          
328 >          
329 > // //     dAtom->setA( A );
330 > // //     dAtom->setJx( ji[0] );
331 > // //     dAtom->setJy( ji[1] );
332 > // //     dAtom->setJz( ji[2] );
333 > // //   }
334 > //       }
335        
336        // calculate the forces
337        
# Line 373 | Line 373 | void Symplectic::integrate( void ){
373          atoms[j]->set_vz(Vz[j]);
374        }
375        
376 <      for( i=0; i< nAtoms; i++ ){
376 > //       for( i=0; i< nAtoms; i++ ){
377  
378 <        if( atoms[i]->isDirectional() ){
378 > //      if( atoms[i]->isDirectional() ){
379  
380 <          dAtom = (DirectionalAtom *)atoms[i];
380 > //        dAtom = (DirectionalAtom *)atoms[i];
381            
382 <          // get and convert the torque to body frame
382 > //        // get and convert the torque to body frame
383            
384 <          Tb[0] = dAtom->getTx();
385 <          Tb[1] = dAtom->getTy();
386 <          Tb[2] = dAtom->getTz();
384 > //        Tb[0] = dAtom->getTx();
385 > //        Tb[1] = dAtom->getTy();
386 > //        Tb[2] = dAtom->getTz();
387            
388 <          dAtom->lab2Body( Tb );
388 > //        dAtom->lab2Body( Tb );
389            
390 <          // get the angular momentum, and complete the angular momentum
391 <          // half step
390 > //        // get the angular momentum, and complete the angular momentum
391 > //        // half step
392            
393 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
394 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
395 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
393 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
394 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
395 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
396            
397 <          dAtom->setJx( ji[0] );
398 <          dAtom->setJy( ji[1] );
399 <          dAtom->setJz( ji[2] );
400 <        }
401 <      }
397 > //        dAtom->setJx( ji[0] );
398 > //        dAtom->setJy( ji[1] );
399 > //        dAtom->setJz( ji[2] );
400 > //      }
401 > //       }
402      
403        time = tl + 1;
404        
# Line 450 | Line 450 | void Symplectic::integrate( void ){
450          atoms[i]->set_vy( vy );
451          atoms[i]->set_vz( vz );
452          
453 <        if( atoms[i]->isDirectional() ){
453 > //      if( atoms[i]->isDirectional() ){
454  
455 <          dAtom = (DirectionalAtom *)atoms[i];
455 > //        dAtom = (DirectionalAtom *)atoms[i];
456            
457 <          // get and convert the torque to body frame
457 > //        // get and convert the torque to body frame
458            
459 <          Tb[0] = dAtom->getTx();
460 <          Tb[1] = dAtom->getTy();
461 <          Tb[2] = dAtom->getTz();
459 > //        Tb[0] = dAtom->getTx();
460 > //        Tb[1] = dAtom->getTy();
461 > //        Tb[2] = dAtom->getTz();
462            
463 <          dAtom->lab2Body( Tb );
463 > //        dAtom->lab2Body( Tb );
464            
465 <          // get the angular momentum, and propagate a half step
465 > //        // get the angular momentum, and propagate a half step
466            
467 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
468 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
469 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
470 <          
471 <          // get the atom's rotation matrix
467 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
468 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
469 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
470            
471 <          A[0][0] = dAtom->getAxx();
474 <          A[0][1] = dAtom->getAxy();
475 <          A[0][2] = dAtom->getAxz();
471 > //        // get the atom's rotation matrix
472            
473 <          A[1][0] = dAtom->getAyx();
474 <          A[1][1] = dAtom->getAyy();
475 <          A[1][2] = dAtom->getAyz();
473 > //        A[0][0] = dAtom->getAxx();
474 > //        A[0][1] = dAtom->getAxy();
475 > //        A[0][2] = dAtom->getAxz();
476            
477 <          A[2][0] = dAtom->getAzx();
478 <          A[2][1] = dAtom->getAzy();
479 <          A[2][2] = dAtom->getAzz();
477 > //        A[1][0] = dAtom->getAyx();
478 > //        A[1][1] = dAtom->getAyy();
479 > //        A[1][2] = dAtom->getAyz();
480            
481 + //        A[2][0] = dAtom->getAzx();
482 + //        A[2][1] = dAtom->getAzy();
483 + //        A[2][2] = dAtom->getAzz();
484            
486 –          // use the angular velocities to propagate the rotation matrix a
487 –          // full time step
485            
486 + //        // use the angular velocities to propagate the rotation matrix a
487 + //        // full time step
488            
490 –          angle = dt2 * ji[0] / dAtom->getIxx();
491 –          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
489            
490 <          angle = dt2 * ji[1] / dAtom->getIyy();
491 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
490 > //        angle = dt2 * ji[0] / dAtom->getIxx();
491 > //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
492            
493 <          angle = dt * ji[2] / dAtom->getIzz();
494 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
493 > //        angle = dt2 * ji[1] / dAtom->getIyy();
494 > //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
495            
496 <          angle = dt2 * ji[1] / dAtom->getIyy();
497 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
496 > //        angle = dt * ji[2] / dAtom->getIzz();
497 > //        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
498            
499 <          angle = dt2 * ji[0] / dAtom->getIxx();
500 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
499 > //        angle = dt2 * ji[1] / dAtom->getIyy();
500 > //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
501            
502 + //        angle = dt2 * ji[0] / dAtom->getIxx();
503 + //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
504            
505 <          dAtom->setA( A );
506 <          dAtom->setJx( ji[0] );
507 <          dAtom->setJy( ji[1] );
508 <          dAtom->setJz( ji[2] );
509 <        }
505 >          
506 > //        dAtom->setA( A );
507 > //        dAtom->setJx( ji[0] );
508 > //        dAtom->setJy( ji[1] );
509 > //        dAtom->setJz( ji[2] );
510 > //      }
511        }
512        
513        // calculate the forces
# Line 533 | Line 533 | void Symplectic::integrate( void ){
533   //      vy2 = vy * vy;
534   //      vz2 = vz * vz;
535          
536 <        if( atoms[i]->isDirectional() ){
536 > //      if( atoms[i]->isDirectional() ){
537  
538 <          dAtom = (DirectionalAtom *)atoms[i];
538 > //        dAtom = (DirectionalAtom *)atoms[i];
539            
540 <          // get and convert the torque to body frame
540 > //        // get and convert the torque to body frame
541            
542 <          Tb[0] = dAtom->getTx();
543 <          Tb[1] = dAtom->getTy();
544 <          Tb[2] = dAtom->getTz();
542 > //        Tb[0] = dAtom->getTx();
543 > //        Tb[1] = dAtom->getTy();
544 > //        Tb[2] = dAtom->getTz();
545            
546 <          dAtom->lab2Body( Tb );
546 > //        dAtom->lab2Body( Tb );
547            
548 <          // get the angular momentum, and complete the angular momentum
549 <          // half step
548 > //        // get the angular momentum, and complete the angular momentum
549 > //        // half step
550            
551 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
552 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
553 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
551 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
552 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
553 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
554            
555 <          jx2 = ji[0] * ji[0];
556 <          jy2 = ji[1] * ji[1];
557 <          jz2 = ji[2] * ji[2];
555 > //        jx2 = ji[0] * ji[0];
556 > //        jy2 = ji[1] * ji[1];
557 > //        jz2 = ji[2] * ji[2];
558            
559 <          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
560 <            + (jz2 / dAtom->getIzz());
559 > //        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
560 > //          + (jz2 / dAtom->getIzz());
561            
562 <          dAtom->setJx( ji[0] );
563 <          dAtom->setJy( ji[1] );
564 <          dAtom->setJz( ji[2] );
565 <        }
562 > //        dAtom->setJx( ji[0] );
563 > //        dAtom->setJy( ji[1] );
564 > //        dAtom->setJz( ji[2] );
565 > //      }
566        }
567        
568        time = tl + 1;

Diff Legend

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