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

Comparing trunk/OOPSE/libmdtools/Verlet.cpp (file contents):
Revision 471 by gezelter, Mon Apr 7 20:51:59 2003 UTC vs.
Revision 497 by chuckv, Mon Apr 14 21:16:37 2003 UTC

# Line 185 | Line 185 | void Verlet::integrate( void ){
185    
186    int time;
187  
188 +  double press[9];
189 +
190    double dt = entry_plug->dt;
191    double runTime = entry_plug->run_time;
192    double sampleTime = entry_plug->sampleTime;
# Line 202 | Line 204 | void Verlet::integrate( void ){
204    DumpWriter*  dump_out = new DumpWriter( entry_plug );
205  
206    // the first time integrate is called, the forces need to be initialized
205 –
207    
208    myFF->doForces(1,1);
209    
# Line 215 | Line 216 | void Verlet::integrate( void ){
216    e_out->writeStat( 0.0 );
217  
218    calcPot = 0;
218 –  calcStress = 0;
219  
220 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
221 +    calcStress = 1;
222 +  } else {
223 +    calcStress = 0;
224 +  }
225 +
226    if( c_is_constrained ){
227      for(i = 0; i < n_loops; i++){
228        
229 <      if (!strcasecmp( entry_plug->ensemble, "NVT"))
230 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
231 <
226 <      // fill R, V, and F arrays and RATTLE in fortran
227 <
229 >      
230 >      // fill R, V, and F arrays and RATTLE in fortran
231 >      
232        for( j=0; j<c_natoms; j++ ){
233 <
233 >        
234          Rx[j] = c_atoms[j]->getX();
235          Ry[j] = c_atoms[j]->getY();
236          Rz[j] = c_atoms[j]->getZ();
# Line 258 | Line 262 | void Verlet::integrate( void ){
262          c_atoms[j]->set_vz(Vz[j]);
263        }
264  
265 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
266 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
267 +  
268        // calculate the forces
269        
270        myFF->doForces(calcPot,calcStress);
# Line 298 | Line 305 | void Verlet::integrate( void ){
305        }
306        
307        if (!strcasecmp( entry_plug->ensemble, "NVT"))
308 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
309 <
310 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
308 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
309 >      
310 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {        
311 >        tStats->getPressureTensor(press);
312          myES->NoseHooverAndersonNPT( dt,
313                                       tStats->getKinetic(),
314 <                                     tStats->getPressure());
314 >                                     press);
315 >      }
316  
317        time = i + 1;
318        
# Line 311 | Line 320 | void Verlet::integrate( void ){
320          if( !(time % vel_n) ) tStats->velocitize();
321        }
322        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
323 +
324        if( !((time+1) % status_n) ) {
325          calcPot = 1;
326 <        calcStress = 1;
326 >        calcStress = 1;
327        }
328        if( !(time % status_n) ){
329          e_out->writeStat( time * dt );
330          calcPot = 0;
331 <        calcStress = 0;
331 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
332 >        else calcStress = 0;
333        }
334      }
335    }
336    else{
337      for(i = 0; i < n_loops; i++){
338  
328 –      if (!strcasecmp( entry_plug->ensemble, "NVT"))
329 –        myES->NoseHooverNVT( dt , tStats->getKinetic() );
339      
340        move_a( dt );
341 +
342 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
343 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
344        
345        // calculate the forces
346        
# Line 339 | Line 351 | void Verlet::integrate( void ){
351        move_b( dt );
352  
353        if (!strcasecmp( entry_plug->ensemble, "NVT"))
354 <        myES->NoseHooverNVT( dt , tStats->getKinetic() );
354 >        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
355  
356 <      if (!strcasecmp( entry_plug->ensemble, "NPT") )
356 >      if (!strcasecmp( entry_plug->ensemble, "NPT") ) {
357 >        tStats->getPressureTensor(press);
358          myES->NoseHooverAndersonNPT( dt,
359                                       tStats->getKinetic(),
360 <                                     tStats->getPressure());
360 >                                     press);
361 >      }
362  
363        time = i + 1;
364        
# Line 354 | Line 368 | void Verlet::integrate( void ){
368        if( !(time % sample_n) )  dump_out->writeDump( time * dt );
369        if( !((time+1) % status_n) ) {
370          calcPot = 1;
371 <        calcStress = 1;
371 >        calcStress = 1;
372        }
373        if( !(time % status_n) ){
374          e_out->writeStat( time * dt );
375          calcPot = 0;
376 <        calcStress = 0;
376 >        if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1;
377 >        else calcStress = 0;
378        }
379      }
380    }

Diff Legend

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