7 |
|
#include "SimInfo.hpp" |
8 |
|
#include "Thermo.hpp" |
9 |
|
#include "ReadWrite.hpp" |
10 |
+ |
#include "ExtendedSystem.hpp" |
11 |
|
|
12 |
|
extern "C"{ |
13 |
|
|
30 |
|
} |
31 |
|
|
32 |
|
|
33 |
< |
Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){ |
33 |
> |
Verlet::Verlet( SimInfo &info, ForceFields* the_ff, ExtendedSystem* the_es ){ |
34 |
|
|
35 |
|
// get what information we need from the SimInfo object |
36 |
|
|
37 |
|
entry_plug = &info; |
38 |
|
myFF = the_ff; |
39 |
< |
|
39 |
> |
myES = the_es; |
40 |
|
|
41 |
|
c_natoms = info.n_atoms; |
42 |
|
c_atoms = info.atoms; |
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(); |
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 / 2.0 , tStats->getKinetic() ); |
230 |
+ |
|
231 |
|
// fill R, V, and F arrays and RATTLE in fortran |
232 |
< |
|
232 |
> |
|
233 |
|
for( j=0; j<c_natoms; j++ ){ |
234 |
< |
|
234 |
> |
|
235 |
|
Rx[j] = c_atoms[j]->getX(); |
236 |
|
Ry[j] = c_atoms[j]->getY(); |
237 |
|
Rz[j] = c_atoms[j]->getZ(); |
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 |
|
|
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, |
302 |
|
c_atoms[j]->set_vz(Vz[j]); |
303 |
|
} |
304 |
|
|
305 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
306 |
+ |
myES->NoseHooverNVT( dt / 2.0, 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 / 2.0, 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 / 2.0 , 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 |
|
|