167 |
|
void Verlet::integrate( void ){ |
168 |
|
|
169 |
|
int i, j; /* loop counters */ |
170 |
< |
int calcPot; |
170 |
> |
int calcPot, calcStress; |
171 |
|
|
172 |
|
double kE; |
173 |
|
|
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(); |
215 |
|
e_out->writeStat( 0.0 ); |
216 |
|
|
217 |
|
calcPot = 0; |
218 |
+ |
calcStress = 0; |
219 |
|
|
220 |
|
if( c_is_constrained ){ |
221 |
|
for(i = 0; i < n_loops; i++){ |
222 |
|
|
223 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
224 |
+ |
myES->NoseHooverNVT( dt , tStats->getKinetic() ); |
225 |
+ |
|
226 |
|
// fill R, V, and F arrays and RATTLE in fortran |
227 |
|
|
228 |
|
for( j=0; j<c_natoms; j++ ){ |
260 |
|
|
261 |
|
// calculate the forces |
262 |
|
|
263 |
< |
myFF->doForces(calcPot,0); |
263 |
> |
myFF->doForces(calcPot,calcStress); |
264 |
|
|
265 |
|
// finish the constrain move ( same as above. ) |
266 |
|
|
279 |
|
Fz[j] = c_atoms[j]->getFz(); |
280 |
|
} |
281 |
|
|
282 |
+ |
|
283 |
|
v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz, |
284 |
|
Fx, Fy, Fz, |
285 |
|
kE, c_n_constrained, c_constrained_dsqr, |
297 |
|
c_atoms[j]->set_vz(Vz[j]); |
298 |
|
} |
299 |
|
|
300 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
301 |
+ |
myES->NoseHooverNVT( dt , tStats->getKinetic() ); |
302 |
+ |
|
303 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) |
304 |
+ |
myES->NoseHooverAndersonNPT( dt, |
305 |
+ |
tStats->getKinetic(), |
306 |
+ |
tStats->getPressure()); |
307 |
+ |
|
308 |
|
time = i + 1; |
309 |
|
|
310 |
|
if( entry_plug->setTemp ){ |
311 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
312 |
|
} |
313 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
314 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
315 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
314 |
> |
if( !((time+1) % status_n) ) { |
315 |
> |
calcPot = 1; |
316 |
> |
calcStress = 1; |
317 |
> |
} |
318 |
> |
if( !(time % status_n) ){ |
319 |
> |
e_out->writeStat( time * dt ); |
320 |
> |
calcPot = 0; |
321 |
> |
calcStress = 0; |
322 |
> |
} |
323 |
|
} |
324 |
|
} |
325 |
|
else{ |
326 |
|
for(i = 0; i < n_loops; i++){ |
327 |
< |
|
327 |
> |
|
328 |
> |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
329 |
> |
myES->NoseHooverNVT( dt , tStats->getKinetic() ); |
330 |
> |
|
331 |
|
move_a( dt ); |
332 |
|
|
333 |
|
// calculate the forces |
334 |
|
|
335 |
< |
myFF->doForces(calcPot,0); |
335 |
> |
myFF->doForces(calcPot,calcStress); |
336 |
|
|
337 |
|
// complete the verlet move |
338 |
|
|
339 |
|
move_b( dt ); |
340 |
|
|
341 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
342 |
+ |
myES->NoseHooverNVT( dt , tStats->getKinetic() ); |
343 |
+ |
|
344 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) |
345 |
+ |
myES->NoseHooverAndersonNPT( dt, |
346 |
+ |
tStats->getKinetic(), |
347 |
+ |
tStats->getPressure()); |
348 |
+ |
|
349 |
|
time = i + 1; |
350 |
|
|
351 |
|
if( entry_plug->setTemp ){ |
352 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
353 |
|
} |
354 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
355 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
356 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
355 |
> |
if( !((time+1) % status_n) ) { |
356 |
> |
calcPot = 1; |
357 |
> |
calcStress = 1; |
358 |
> |
} |
359 |
> |
if( !(time % status_n) ){ |
360 |
> |
e_out->writeStat( time * dt ); |
361 |
> |
calcPot = 0; |
362 |
> |
calcStress = 0; |
363 |
> |
} |
364 |
|
} |
365 |
|
} |
366 |
|
|