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 483 by gezelter, Wed Apr 9 04:06:43 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 219 | 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 235 | Line 235 | void Symplectic::integrate( void ){
235      
236  
237      for( tl=0; tl < n_loops; tl++ ){
238 +
239 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
240 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
241        
242        for( j=0; j<nAtoms; j++ ){
243  
# Line 400 | Line 403 | void Symplectic::integrate( void ){
403          }
404        }
405      
406 +
407 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
408 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
409 +
410 +      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
411 +        tStats->getPressureTensor(press);
412 +        myES->NoseHooverAndersonNPT( dt,
413 +                                     tStats->getKinetic(),
414 +                                     press);
415 +      }
416 +
417        time = tl + 1;
418        
419        if( entry_plug->setTemp ){
# Line 408 | 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 <        calcStress = 1;
425 >        // bitwise masking in case we need it for NPT
426 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
427        }
428        if( !(time % status_n) ){
429          e_out->writeStat( time * dt );
430          calcPot = 0;
431 <        calcStress = 0;
431 >        // bitwise masking in case we need it for NPT
432 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
433        }
434      }
435    }
# Line 424 | Line 440 | void Symplectic::integrate( void ){
440        kE = 0.0;
441        rot_kE= 0.0;
442        trans_kE = 0.0;
443 +
444 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
445 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
446        
447        for( i=0; i<nAtoms; i++ ){
448          
# Line 529 | Line 548 | void Symplectic::integrate( void ){
548          atoms[i]->set_vy( vy );
549          atoms[i]->set_vz( vz );
550          
551 < //      vx2 = vx * vx;
552 < //      vy2 = vy * vy;
553 < //      vz2 = vz * vz;
551 >        vx2 = vx * vx;
552 >        vy2 = vy * vy;
553 >        vz2 = vz * vz;
554          
555          if( atoms[i]->isDirectional() ){
556  
# Line 562 | Line 581 | void Symplectic::integrate( void ){
581            dAtom->setJx( ji[0] );
582            dAtom->setJy( ji[1] );
583            dAtom->setJz( ji[2] );
584 <        }
584 >        }
585        }
586 <      
586 >    
587 >      if (!strcasecmp( entry_plug->ensemble, "NVT"))
588 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
589 >
590 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
591 >        tStats->getPressureTensor(press);
592 >        myES->NoseHooverAndersonNPT( dt,
593 >                                     tStats->getKinetic(),
594 >                                     press);
595 >      }
596 >  
597        time = tl + 1;
598        
599        if( entry_plug->setTemp ){
# Line 573 | 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 <        calcStress = 1;
605 >        // bitwise masking in case we need it for NPT
606 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
607        }
608        if( !(time % status_n) ){
609          e_out->writeStat( time * dt );
610          calcPot = 0;
611 <        calcStress = 0;
611 >        // bitwise masking in case we need it for NPT
612 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
613        }      
614      }
615    }

Diff Legend

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