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 466 by gezelter, Mon Apr 7 14:30:36 2003 UTC vs.
Revision 475 by gezelter, Tue Apr 8 12:44:18 2003 UTC

# Line 167 | Line 167 | void Verlet::integrate( void ){
167   void Verlet::integrate( void ){
168    
169    int i, j; /* loop counters */
170 <  int calcPot;
170 >  int calcPot, calcStress;
171  
172    double kE;
173  
# Line 204 | Line 204 | void Verlet::integrate( void ){
204    // the first time integrate is called, the forces need to be initialized
205  
206    
207 <  myFF->doForces(1,0);
207 >  myFF->doForces(1,1);
208    
209    if( entry_plug->setTemp ){
210      tStats->velocitize();
# Line 215 | Line 215 | void Verlet::integrate( void ){
215    e_out->writeStat( 0.0 );
216  
217    calcPot = 0;
218 +
219 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
220 +    calcStress = 1;
221 +  } else {
222 +    calcStress = 0;
223 +  }
224  
225    if( c_is_constrained ){
226      for(i = 0; i < n_loops; i++){
227        
228 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
229 +        myES->NoseHooverNVT( dt , tStats->getKinetic() );
230 +
231        // fill R, V, and F arrays and RATTLE in fortran
232  
233        for( j=0; j<c_natoms; j++ ){
# Line 256 | Line 265 | void Verlet::integrate( void ){
265  
266        // calculate the forces
267        
268 <      myFF->doForces(calcPot,0);
268 >      myFF->doForces(calcPot,calcStress);
269        
270        // finish the constrain move ( same as above. )
271  
# Line 275 | Line 284 | void Verlet::integrate( void ){
284          Fz[j] = c_atoms[j]->getFz();
285        }
286          
287 +
288        v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz,
289                        Fx, Fy, Fz,
290                        kE, c_n_constrained, c_constrained_dsqr,
# Line 292 | Line 302 | void Verlet::integrate( void ){
302          c_atoms[j]->set_vz(Vz[j]);
303        }
304        
305 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
306 +        myES->NoseHooverNVT( dt , tStats->getKinetic() );
307 +
308 +      if (!strcasecmp( entry_plug->ensemble, "NPT") )
309 +        myES->NoseHooverAndersonNPT( dt,
310 +                                     tStats->getKinetic(),
311 +                                     tStats->getPressure());
312 +
313        time = i + 1;
314        
315        if( entry_plug->setTemp ){
316          if( !(time % vel_n) ) tStats->velocitize();
317        }
318        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
319 <      if( !((time+1) % status_n) ) calcPot = 1;
320 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
319 >      if( !((time+1) % status_n) ) {
320 >        calcPot = 1;
321 >        // bitwise masking in case we need it for NPT
322 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
323 >      }
324 >      if( !(time % status_n) ){
325 >        e_out->writeStat( time * dt );
326 >        calcPot = 0;
327 >        // bitwise masking in case we need it for NPT
328 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
329 >      }
330      }
331    }
332    else{
333      for(i = 0; i < n_loops; i++){
334 <      
334 >
335 >      if (!strcasecmp( entry_plug->ensemble, "NVT"))
336 >        myES->NoseHooverNVT( dt , tStats->getKinetic() );
337 >    
338        move_a( dt );
339        
340        // calculate the forces
341        
342 <      myFF->doForces(calcPot,0);
342 >      myFF->doForces(calcPot,calcStress);
343              
344        // complete the verlet move
345        
346        move_b( dt );
347  
348 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
349 +        myES->NoseHooverNVT( dt , tStats->getKinetic() );
350 +
351 +      if (!strcasecmp( entry_plug->ensemble, "NPT") )
352 +        myES->NoseHooverAndersonNPT( dt,
353 +                                     tStats->getKinetic(),
354 +                                     tStats->getPressure());
355 +
356        time = i + 1;
357        
358        if( entry_plug->setTemp ){
359          if( !(time % vel_n) ) tStats->velocitize();
360        }
361        if( !(time % sample_n) )  dump_out->writeDump( time * dt );
362 <      if( !((time+1) % status_n) ) calcPot = 1;
363 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
362 >      if( !((time+1) % status_n) ) {
363 >        calcPot = 1;
364 >        // bitwise masking in case we need it for NPT
365 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
366 >      }
367 >      if( !(time % status_n) ){
368 >        e_out->writeStat( time * dt );
369 >        calcPot = 0;
370 >        // bitwise masking in case we need it for NPT
371 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
372 >      }
373      }
374    }
375    

Diff Legend

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