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