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 486 by mmeineke, Thu Apr 10 16:22:00 2003 UTC

# Line 38 | Line 38 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
38    myES = the_es;
39    isFirst = 1;
40  
41 –  std::cerr<< "calling symplectic constructor\n";
42 –
41    molecules = entry_plug->molecules;
42    nMols = entry_plug->n_mol;
43  
# Line 164 | Line 162 | void Symplectic::integrate( void ){
162    double dt2;                       // half the dt
163  
164    double vx, vy, vz;    // the velocities
165 < //  double vx2, vy2, vz2; // the square of the velocities
165 >  double vx2, vy2, vz2; // the square of the velocities
166    double rx, ry, rz;    // the postitions
167    
168    double ji[3];   // the body frame angular momentum
# Line 172 | Line 170 | void Symplectic::integrate( void ){
170    double Tb[3];   // torque in the body frame
171    double angle;   // the angle through which to rotate the rotation matrix
172    double A[3][3]; // the rotation matrix
173 +  double press[9];
174  
175    int time;
176  
# Line 192 | Line 191 | void Symplectic::integrate( void ){
191    StatWriter*  e_out;
192    DumpWriter*  dump_out;
193  
195 –  std::cerr << "about to call new thermo\n";
196 –
194    tStats   = new Thermo( entry_plug );
195    e_out    = new StatWriter( entry_plug );
199 –
200 –  std::cerr << "calling dumpWriter \n";
196    dump_out = new DumpWriter( entry_plug );
202 –  std::cerr << "called dumpWriter \n";
197  
198    Atom** atoms = entry_plug->atoms;
199    DirectionalAtom* dAtom;
# Line 243 | Line 237 | void Symplectic::integrate( void ){
237      for( tl=0; tl < n_loops; tl++ ){
238  
239        if (!strcasecmp( entry_plug->ensemble, "NVT"))
240 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
240 >        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
241        
242        for( j=0; j<nAtoms; j++ ){
243  
# Line 281 | Line 275 | void Symplectic::integrate( void ){
275        }
276  
277  
278 < //       for( i=0; i<nAtoms; i++ ){
279 < // //   if( atoms[i]->isDirectional() ){
278 >      for( i=0; i<nAtoms; i++ ){
279 >        if( atoms[i]->isDirectional() ){
280                    
281 < // //     dAtom = (DirectionalAtom *)atoms[i];
281 >          dAtom = (DirectionalAtom *)atoms[i];
282 >          
283 >          // get and convert the torque to body frame
284            
285 < // //     // get and convert the torque to body frame
285 >          Tb[0] = dAtom->getTx();
286 >          Tb[1] = dAtom->getTy();
287 >          Tb[2] = dAtom->getTz();
288            
289 < // //     Tb[0] = dAtom->getTx();
292 < // //     Tb[1] = dAtom->getTy();
293 < // //     Tb[2] = dAtom->getTz();
289 >          dAtom->lab2Body( Tb );
290            
291 < // //     dAtom->lab2Body( Tb );
296 <          
297 < // //     // get the angular momentum, and propagate a half step
291 >          // get the angular momentum, and propagate a half step
292            
293 < // //     ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
294 < // //     ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
295 < // //     ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
293 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
294 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
295 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
296            
297 < // //     // get the atom's rotation matrix
297 >          // get the atom's rotation matrix
298            
299 < // //     A[0][0] = dAtom->getAxx();
300 < // //     A[0][1] = dAtom->getAxy();
301 < // //     A[0][2] = dAtom->getAxz();
299 >          A[0][0] = dAtom->getAxx();
300 >          A[0][1] = dAtom->getAxy();
301 >          A[0][2] = dAtom->getAxz();
302            
303 < // //     A[1][0] = dAtom->getAyx();
304 < // //     A[1][1] = dAtom->getAyy();
305 < // //     A[1][2] = dAtom->getAyz();
303 >          A[1][0] = dAtom->getAyx();
304 >          A[1][1] = dAtom->getAyy();
305 >          A[1][2] = dAtom->getAyz();
306            
307 < // //     A[2][0] = dAtom->getAzx();
308 < // //     A[2][1] = dAtom->getAzy();
309 < // //     A[2][2] = dAtom->getAzz();
307 >          A[2][0] = dAtom->getAzx();
308 >          A[2][1] = dAtom->getAzy();
309 >          A[2][2] = dAtom->getAzz();
310            
311            
312 < // //     // use the angular velocities to propagate the rotation matrix a
313 < // //     // full time step
312 >          // use the angular velocities to propagate the rotation matrix a
313 >          // full time step
314            
315            
316 < // //     angle = dt2 * ji[0] / dAtom->getIxx();
317 < // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
316 >          angle = dt2 * ji[0] / dAtom->getIxx();
317 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
318            
319 < // //     angle = dt2 * ji[1] / dAtom->getIyy();
320 < // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
319 >          angle = dt2 * ji[1] / dAtom->getIyy();
320 >          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 = dt * ji[2] / dAtom->getIzz();
323 >          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-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 = dt2 * ji[0] / dAtom->getIxx();
329 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
330            
331            
332 < // //     dAtom->setA( A );
333 < // //     dAtom->setJx( ji[0] );
334 < // //     dAtom->setJy( ji[1] );
335 < // //     dAtom->setJz( ji[2] );
336 < // //   }
337 < //       }
332 >          dAtom->setA( A );
333 >          dAtom->setJx( ji[0] );
334 >          dAtom->setJy( ji[1] );
335 >          dAtom->setJz( ji[2] );
336 >        }
337 >      }
338        
339        // calculate the forces
340        
# Line 382 | Line 376 | void Symplectic::integrate( void ){
376          atoms[j]->set_vz(Vz[j]);
377        }
378        
379 < //       for( i=0; i< nAtoms; i++ ){
386 <
387 < //      if( atoms[i]->isDirectional() ){
379 >      for( i=0; i< nAtoms; i++ ){
380  
381 < //        dAtom = (DirectionalAtom *)atoms[i];
381 >        if( atoms[i]->isDirectional() ){
382 >
383 >          dAtom = (DirectionalAtom *)atoms[i];
384            
385 < //        // get and convert the torque to body frame
385 >          // get and convert the torque to body frame
386            
387 < //        Tb[0] = dAtom->getTx();
388 < //        Tb[1] = dAtom->getTy();
389 < //        Tb[2] = dAtom->getTz();
387 >          Tb[0] = dAtom->getTx();
388 >          Tb[1] = dAtom->getTy();
389 >          Tb[2] = dAtom->getTz();
390            
391 < //        dAtom->lab2Body( Tb );
391 >          dAtom->lab2Body( Tb );
392            
393 < //        // get the angular momentum, and complete the angular momentum
394 < //        // half step
393 >          // get the angular momentum, and complete the angular momentum
394 >          // half step
395            
396 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
397 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
398 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
396 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
397 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
398 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
399            
400 < //        dAtom->setJx( ji[0] );
401 < //        dAtom->setJy( ji[1] );
402 < //        dAtom->setJz( ji[2] );
403 < //      }
404 < //       }
400 >          dAtom->setJx( ji[0] );
401 >          dAtom->setJy( ji[1] );
402 >          dAtom->setJz( ji[2] );
403 >        }
404 >      }
405      
406  
407        if (!strcasecmp( entry_plug->ensemble, "NVT"))
408 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
408 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
409  
410 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
410 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
411 >        tStats->getPressureTensor(press);
412          myES->NoseHooverAndersonNPT( dt,
413                                       tStats->getKinetic(),
414 <                                     tStats->getPressure());
414 >                                     press);
415 >      }
416  
417        time = tl + 1;
418        
# Line 426 | Line 422 | void Symplectic::integrate( void ){
422        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
423        if( !((time+1) % status_n) ) {
424          calcPot = 1;
425 <        // bitwise masking in case we need it for NPT
430 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
425 >        calcStress = 1;
426        }
427        if( !(time % status_n) ){
428          e_out->writeStat( time * dt );
429          calcPot = 0;
430 <        // bitwise masking in case we need it for NPT
431 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
430 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
431 >        else calcStress = 0;
432        }
433      }
434    }
# Line 446 | Line 441 | void Symplectic::integrate( void ){
441        trans_kE = 0.0;
442  
443        if (!strcasecmp( entry_plug->ensemble, "NVT"))
444 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
444 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
445        
446        for( i=0; i<nAtoms; i++ ){
447          
# Line 473 | Line 468 | void Symplectic::integrate( void ){
468          atoms[i]->set_vy( vy );
469          atoms[i]->set_vz( vz );
470          
471 < //      if( atoms[i]->isDirectional() ){
471 >        if( atoms[i]->isDirectional() ){
472  
473 < //        dAtom = (DirectionalAtom *)atoms[i];
473 >          dAtom = (DirectionalAtom *)atoms[i];
474            
475 < //        // get and convert the torque to body frame
475 >          // get and convert the torque to body frame
476            
477 < //        Tb[0] = dAtom->getTx();
478 < //        Tb[1] = dAtom->getTy();
479 < //        Tb[2] = dAtom->getTz();
477 >          Tb[0] = dAtom->getTx();
478 >          Tb[1] = dAtom->getTy();
479 >          Tb[2] = dAtom->getTz();
480            
481 < //        dAtom->lab2Body( Tb );
481 >          dAtom->lab2Body( Tb );
482            
483 < //        // get the angular momentum, and propagate a half step
483 >          // get the angular momentum, and propagate a half step
484            
485 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
486 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
487 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
485 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
486 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
487 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
488            
489 < //        // get the atom's rotation matrix
489 >          // get the atom's rotation matrix
490            
491 < //        A[0][0] = dAtom->getAxx();
492 < //        A[0][1] = dAtom->getAxy();
493 < //        A[0][2] = dAtom->getAxz();
491 >          A[0][0] = dAtom->getAxx();
492 >          A[0][1] = dAtom->getAxy();
493 >          A[0][2] = dAtom->getAxz();
494            
495 < //        A[1][0] = dAtom->getAyx();
496 < //        A[1][1] = dAtom->getAyy();
497 < //        A[1][2] = dAtom->getAyz();
495 >          A[1][0] = dAtom->getAyx();
496 >          A[1][1] = dAtom->getAyy();
497 >          A[1][2] = dAtom->getAyz();
498            
499 < //        A[2][0] = dAtom->getAzx();
500 < //        A[2][1] = dAtom->getAzy();
501 < //        A[2][2] = dAtom->getAzz();
499 >          A[2][0] = dAtom->getAzx();
500 >          A[2][1] = dAtom->getAzy();
501 >          A[2][2] = dAtom->getAzz();
502            
503            
504 < //        // use the angular velocities to propagate the rotation matrix a
505 < //        // full time step
504 >          // use the angular velocities to propagate the rotation matrix a
505 >          // full time step
506            
507            
508 < //        angle = dt2 * ji[0] / dAtom->getIxx();
509 < //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
508 >          angle = dt2 * ji[0] / dAtom->getIxx();
509 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
510            
511 < //        angle = dt2 * ji[1] / dAtom->getIyy();
512 < //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
511 >          angle = dt2 * ji[1] / dAtom->getIyy();
512 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
513            
514 < //        angle = dt * ji[2] / dAtom->getIzz();
515 < //        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
514 >          angle = dt * ji[2] / dAtom->getIzz();
515 >          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
516            
517 < //        angle = dt2 * ji[1] / dAtom->getIyy();
518 < //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
517 >          angle = dt2 * ji[1] / dAtom->getIyy();
518 >          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
519            
520 < //        angle = dt2 * ji[0] / dAtom->getIxx();
521 < //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
520 >          angle = dt2 * ji[0] / dAtom->getIxx();
521 >          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
522            
523            
524 < //        dAtom->setA( A );
525 < //        dAtom->setJx( ji[0] );
526 < //        dAtom->setJy( ji[1] );
527 < //        dAtom->setJz( ji[2] );
528 < //      }
524 >          dAtom->setA( A );
525 >          dAtom->setJx( ji[0] );
526 >          dAtom->setJy( ji[1] );
527 >          dAtom->setJz( ji[2] );
528 >        }
529        }
530        
531        // calculate the forces
# Line 552 | Line 547 | void Symplectic::integrate( void ){
547          atoms[i]->set_vy( vy );
548          atoms[i]->set_vz( vz );
549          
550 < //      vx2 = vx * vx;
551 < //      vy2 = vy * vy;
552 < //      vz2 = vz * vz;
550 >        vx2 = vx * vx;
551 >        vy2 = vy * vy;
552 >        vz2 = vz * vz;
553          
554 < //      if( atoms[i]->isDirectional() ){
554 >        if( atoms[i]->isDirectional() ){
555  
556 < //        dAtom = (DirectionalAtom *)atoms[i];
556 >          dAtom = (DirectionalAtom *)atoms[i];
557            
558 < //        // get and convert the torque to body frame
558 >          // get and convert the torque to body frame
559            
560 < //        Tb[0] = dAtom->getTx();
561 < //        Tb[1] = dAtom->getTy();
562 < //        Tb[2] = dAtom->getTz();
560 >          Tb[0] = dAtom->getTx();
561 >          Tb[1] = dAtom->getTy();
562 >          Tb[2] = dAtom->getTz();
563            
564 < //        dAtom->lab2Body( Tb );
564 >          dAtom->lab2Body( Tb );
565            
566 < //        // get the angular momentum, and complete the angular momentum
567 < //        // half step
566 >          // get the angular momentum, and complete the angular momentum
567 >          // half step
568            
569 < //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
570 < //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
571 < //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
569 >          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
570 >          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
571 >          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
572            
573 < //        jx2 = ji[0] * ji[0];
574 < //        jy2 = ji[1] * ji[1];
575 < //        jz2 = ji[2] * ji[2];
573 >          jx2 = ji[0] * ji[0];
574 >          jy2 = ji[1] * ji[1];
575 >          jz2 = ji[2] * ji[2];
576            
577 < //        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
578 < //          + (jz2 / dAtom->getIzz());
577 >          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
578 >            + (jz2 / dAtom->getIzz());
579            
580 < //        dAtom->setJx( ji[0] );
581 < //        dAtom->setJy( ji[1] );
582 < //        dAtom->setJz( ji[2] );
583 < //      }
580 >          dAtom->setJx( ji[0] );
581 >          dAtom->setJy( ji[1] );
582 >          dAtom->setJz( ji[2] );
583 >        }
584 >
585        }
586      
587        if (!strcasecmp( entry_plug->ensemble, "NVT"))
588 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
588 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
589  
590 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
590 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
591 >        tStats->getPressureTensor(press);
592          myES->NoseHooverAndersonNPT( dt,
593                                       tStats->getKinetic(),
594 <                                     tStats->getPressure());
594 >                                     press);
595 >      }
596    
597        time = tl + 1;
598        
# Line 604 | Line 602 | void Symplectic::integrate( void ){
602        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
603        if( !((time+1) % status_n) ) {
604          calcPot = 1;
605 <        // bitwise masking in case we need it for NPT
608 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
605 >        calcStress = 1;
606        }
607        if( !(time % status_n) ){
608          e_out->writeStat( time * dt );
609          calcPot = 0;
610 <        // bitwise masking in case we need it for NPT
611 <        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
610 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
611 >        else calcStress = 0;
612        }      
613      }
614    }

Diff Legend

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