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; |
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 |
|
|
226 |
|
if( c_is_constrained ){ |
227 |
|
for(i = 0; i < n_loops; i++){ |
228 |
|
|
228 |
– |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
229 |
– |
myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() ); |
229 |
|
|
230 |
|
// fill R, V, and F arrays and RATTLE in fortran |
231 |
|
|
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); |
307 |
|
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
308 |
|
myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() ); |
309 |
|
|
310 |
< |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) |
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 |
|
|
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 |
< |
// bitwise masking in case we need it for NPT |
322 |
< |
calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1; |
326 |
> |
calcStress = 1; |
327 |
|
} |
328 |
|
if( !(time % status_n) ){ |
329 |
|
e_out->writeStat( time * dt ); |
330 |
|
calcPot = 0; |
331 |
< |
// bitwise masking in case we need it for NPT |
332 |
< |
calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 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 |
|
|
335 |
– |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
336 |
– |
myES->NoseHooverNVT( dt / 2.0, 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 |
|
|
353 |
|
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
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 |
|
|
368 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
369 |
|
if( !((time+1) % status_n) ) { |
370 |
|
calcPot = 1; |
371 |
< |
// bitwise masking in case we need it for NPT |
365 |
< |
calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1; |
371 |
> |
calcStress = 1; |
372 |
|
} |
373 |
|
if( !(time % status_n) ){ |
374 |
|
e_out->writeStat( time * dt ); |
375 |
|
calcPot = 0; |
376 |
< |
// bitwise masking in case we need it for NPT |
377 |
< |
calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0; |
376 |
> |
if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1; |
377 |
> |
else calcStress = 0; |
378 |
|
} |
379 |
|
} |
380 |
|
} |