| 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(); |
| 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 |
|
|
