| 167 |
|
void Verlet::integrate( void ){ |
| 168 |
|
|
| 169 |
|
int i, j; /* loop counters */ |
| 170 |
< |
int calcPot; |
| 170 |
> |
int calcPot, calcStress; |
| 171 |
|
|
| 172 |
|
double kE; |
| 173 |
|
|
| 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,0); |
| 208 |
> |
myFF->doForces(1,1); |
| 209 |
|
|
| 210 |
|
if( entry_plug->setTemp ){ |
| 211 |
|
tStats->velocitize(); |
| 217 |
|
|
| 218 |
|
calcPot = 0; |
| 219 |
|
|
| 220 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT")) { |
| 221 |
+ |
calcStress = 1; |
| 222 |
+ |
} else { |
| 223 |
+ |
calcStress = 0; |
| 224 |
+ |
} |
| 225 |
+ |
|
| 226 |
|
if( c_is_constrained ){ |
| 227 |
|
for(i = 0; i < n_loops; i++){ |
| 228 |
|
|
| 229 |
+ |
|
| 230 |
|
// fill R, V, and F arrays and RATTLE in fortran |
| 231 |
< |
|
| 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(); |
| 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,0); |
| 270 |
> |
myFF->doForces(calcPot,calcStress); |
| 271 |
|
|
| 272 |
|
// finish the constrain move ( same as above. ) |
| 273 |
|
|
| 286 |
|
Fz[j] = c_atoms[j]->getFz(); |
| 287 |
|
} |
| 288 |
|
|
| 289 |
+ |
|
| 290 |
|
v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz, |
| 291 |
|
Fx, Fy, Fz, |
| 292 |
|
kE, c_n_constrained, c_constrained_dsqr, |
| 304 |
|
c_atoms[j]->set_vz(Vz[j]); |
| 305 |
|
} |
| 306 |
|
|
| 307 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
| 308 |
+ |
myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() ); |
| 309 |
+ |
|
| 310 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) { |
| 311 |
+ |
tStats->getPressureTensor(press); |
| 312 |
+ |
myES->NoseHooverAndersonNPT( dt, |
| 313 |
+ |
tStats->getKinetic(), |
| 314 |
+ |
press); |
| 315 |
+ |
} |
| 316 |
+ |
|
| 317 |
|
time = i + 1; |
| 318 |
|
|
| 319 |
|
if( entry_plug->setTemp ){ |
| 320 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 321 |
|
} |
| 322 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 323 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
| 324 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 323 |
> |
|
| 324 |
> |
if( !((time+1) % status_n) ) { |
| 325 |
> |
calcPot = 1; |
| 326 |
> |
calcStress = 1; |
| 327 |
> |
} |
| 328 |
> |
if( !(time % status_n) ){ |
| 329 |
> |
e_out->writeStat( time * dt ); |
| 330 |
> |
calcPot = 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 |
< |
|
| 338 |
> |
|
| 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 |
|
|
| 347 |
< |
myFF->doForces(calcPot,0); |
| 347 |
> |
myFF->doForces(calcPot,calcStress); |
| 348 |
|
|
| 349 |
|
// complete the verlet move |
| 350 |
|
|
| 351 |
|
move_b( dt ); |
| 352 |
|
|
| 353 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
| 354 |
+ |
myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() ); |
| 355 |
+ |
|
| 356 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) { |
| 357 |
+ |
tStats->getPressureTensor(press); |
| 358 |
+ |
myES->NoseHooverAndersonNPT( dt, |
| 359 |
+ |
tStats->getKinetic(), |
| 360 |
+ |
press); |
| 361 |
+ |
} |
| 362 |
+ |
|
| 363 |
|
time = i + 1; |
| 364 |
|
|
| 365 |
|
if( entry_plug->setTemp ){ |
| 366 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 367 |
|
} |
| 368 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 369 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
| 370 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 369 |
> |
if( !((time+1) % status_n) ) { |
| 370 |
> |
calcPot = 1; |
| 371 |
> |
calcStress = 1; |
| 372 |
> |
} |
| 373 |
> |
if( !(time % status_n) ){ |
| 374 |
> |
e_out->writeStat( time * dt ); |
| 375 |
> |
calcPot = 0; |
| 376 |
> |
if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1; |
| 377 |
> |
else calcStress = 0; |
| 378 |
> |
} |
| 379 |
|
} |
| 380 |
|
} |
| 381 |
|
|
