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 466 by gezelter, Mon Apr 7 14:30:36 2003 UTC vs.
Revision 497 by chuckv, Mon Apr 14 21:16:37 2003 UTC

# Line 162 | 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 170 | 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 184 | Line 185 | void Symplectic::integrate( void ){
185    int status_n = (int)( statusTime / dt );
186    int vel_n    = (int)( thermalTime / dt );
187  
188 <  int calcPot;
188 >  int calcPot, calcStress;
189  
190 <   Thermo *tStats = new Thermo( entry_plug );
190 >  Thermo *tStats;
191 >  StatWriter*  e_out;
192 >  DumpWriter*  dump_out;
193  
194 <  StatWriter*  e_out    = new StatWriter( entry_plug );
195 <  DumpWriter*  dump_out = new DumpWriter( entry_plug );
194 >  tStats   = new Thermo( entry_plug );
195 >  e_out    = new StatWriter( entry_plug );
196 >  dump_out = new DumpWriter( entry_plug );
197  
198    Atom** atoms = entry_plug->atoms;
199    DirectionalAtom* dAtom;
# Line 197 | Line 201 | void Symplectic::integrate( void ){
201  
202    // initialize the forces the before the first step
203  
204 <  myFF->doForces(1,0);
204 >  myFF->doForces(1,1);
205    
206    if( entry_plug->setTemp ){
207      
# Line 209 | Line 213 | void Symplectic::integrate( void ){
213    
214    calcPot = 0;
215  
216 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
217 +    calcStress = 1;
218 +  } else {
219 +    calcStress = 0;
220 +  }
221 +
222    if( n_constrained ){
223  
224      double *Rx = new double[nAtoms];
# Line 225 | Line 235 | void Symplectic::integrate( void ){
235      
236  
237      for( tl=0; tl < n_loops; tl++ ){
238 +
239        
240        for( j=0; j<nAtoms; j++ ){
241  
# Line 323 | Line 334 | void Symplectic::integrate( void ){
334          }
335        }
336        
337 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
338 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
339 +
340        // calculate the forces
341        
342 <      myFF->doForces(calcPot, 0);
342 >      myFF->doForces(calcPot, calcStress);
343        
344        // move b
345  
# Line 390 | Line 404 | void Symplectic::integrate( void ){
404          }
405        }
406      
407 +
408 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
409 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
410 +
411 +      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
412 +        tStats->getPressureTensor(press);
413 +        myES->NoseHooverAndersonNPT( dt,
414 +                                     tStats->getKinetic(),
415 +                                     press);
416 +      }
417 +
418        time = tl + 1;
419        
420        if( entry_plug->setTemp ){
421          if( !(time % vel_n) ) tStats->velocitize();
422        }
423        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
424 <      if( !((time+1) % status_n) ) calcPot = 1;
425 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
424 >      if( !((time+1) % status_n) ) {
425 >        calcPot = 1;
426 >        calcStress = 1;
427 >      }
428 >      if( !(time % status_n) ){
429 >        e_out->writeStat( time * dt );
430 >        calcPot = 0;
431 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
432 >        else calcStress = 0;
433 >      }
434      }
435    }
436    else{
# Line 407 | Line 440 | void Symplectic::integrate( void ){
440        kE = 0.0;
441        rot_kE= 0.0;
442        trans_kE = 0.0;
443 <      
443 >
444        for( i=0; i<nAtoms; i++ ){
445          
446          // velocity half step
# Line 492 | Line 525 | void Symplectic::integrate( void ){
525            dAtom->setJz( ji[2] );
526          }
527        }
528 +
529 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
530 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
531        
532 +      
533        // calculate the forces
534        
535 <      myFF->doForces(calcPot,0);
535 >      myFF->doForces(calcPot,calcStress);
536        
537        for( i=0; i< nAtoms; i++ ){
538          
# Line 512 | Line 549 | void Symplectic::integrate( void ){
549          atoms[i]->set_vy( vy );
550          atoms[i]->set_vz( vz );
551          
552 < //      vx2 = vx * vx;
553 < //      vy2 = vy * vy;
554 < //      vz2 = vz * vz;
552 >        vx2 = vx * vx;
553 >        vy2 = vy * vy;
554 >        vz2 = vz * vz;
555          
556          if( atoms[i]->isDirectional() ){
557  
# Line 546 | Line 583 | void Symplectic::integrate( void ){
583            dAtom->setJy( ji[1] );
584            dAtom->setJz( ji[2] );
585          }
586 +
587        }
588 <      
588 >    
589 >      if (!strcasecmp( entry_plug->ensemble, "NVT"))
590 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
591 >
592 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
593 >        tStats->getPressureTensor(press);
594 >        myES->NoseHooverAndersonNPT( dt,
595 >                                     tStats->getKinetic(),
596 >                                     press);
597 >      }
598 >  
599        time = tl + 1;
600        
601        if( entry_plug->setTemp ){
602          if( !(time % vel_n) ) tStats->velocitize();
603        }
604        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
605 <      if( !((time+1) % status_n) ) calcPot = 1;
606 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
605 >      if( !((time+1) % status_n) ) {
606 >        calcPot = 1;
607 >        calcStress = 1;
608 >      }
609 >      if( !(time % status_n) ){
610 >        e_out->writeStat( time * dt );
611 >        calcPot = 0;
612 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
613 >        else calcStress = 0;
614 >      }      
615      }
616    }
617  

Diff Legend

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