| 162 |
|
double dt2; // half the dt |
| 163 |
|
|
| 164 |
|
double vx, vy, vz; // the velocities |
| 165 |
< |
// double vx2, vy2, vz2; // the square of the velocities |
| 165 |
> |
double vx2, vy2, vz2; // the square of the velocities |
| 166 |
|
double rx, ry, rz; // the postitions |
| 167 |
|
|
| 168 |
|
double ji[3]; // the body frame angular momentum |
| 170 |
|
double Tb[3]; // torque in the body frame |
| 171 |
|
double angle; // the angle through which to rotate the rotation matrix |
| 172 |
|
double A[3][3]; // the rotation matrix |
| 173 |
+ |
double press[9]; |
| 174 |
|
|
| 175 |
|
int time; |
| 176 |
|
|
| 185 |
|
int status_n = (int)( statusTime / dt ); |
| 186 |
|
int vel_n = (int)( thermalTime / dt ); |
| 187 |
|
|
| 188 |
< |
int calcPot; |
| 188 |
> |
int calcPot, calcStress; |
| 189 |
|
|
| 190 |
< |
Thermo *tStats = new Thermo( entry_plug ); |
| 190 |
> |
Thermo *tStats; |
| 191 |
> |
StatWriter* e_out; |
| 192 |
> |
DumpWriter* dump_out; |
| 193 |
|
|
| 194 |
< |
StatWriter* e_out = new StatWriter( entry_plug ); |
| 195 |
< |
DumpWriter* dump_out = new DumpWriter( entry_plug ); |
| 194 |
> |
tStats = new Thermo( entry_plug ); |
| 195 |
> |
e_out = new StatWriter( entry_plug ); |
| 196 |
> |
dump_out = new DumpWriter( entry_plug ); |
| 197 |
|
|
| 198 |
|
Atom** atoms = entry_plug->atoms; |
| 199 |
|
DirectionalAtom* dAtom; |
| 201 |
|
|
| 202 |
|
// initialize the forces the before the first step |
| 203 |
|
|
| 204 |
< |
myFF->doForces(1,0); |
| 204 |
> |
myFF->doForces(1,1); |
| 205 |
|
|
| 206 |
|
if( entry_plug->setTemp ){ |
| 207 |
|
|
| 213 |
|
|
| 214 |
|
calcPot = 0; |
| 215 |
|
|
| 216 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT")) { |
| 217 |
+ |
calcStress = 1; |
| 218 |
+ |
} else { |
| 219 |
+ |
calcStress = 0; |
| 220 |
+ |
} |
| 221 |
+ |
|
| 222 |
|
if( n_constrained ){ |
| 223 |
|
|
| 224 |
|
double *Rx = new double[nAtoms]; |
| 235 |
|
|
| 236 |
|
|
| 237 |
|
for( tl=0; tl < n_loops; tl++ ){ |
| 238 |
+ |
|
| 239 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
| 240 |
+ |
myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() ); |
| 241 |
|
|
| 242 |
|
for( j=0; j<nAtoms; j++ ){ |
| 243 |
|
|
| 338 |
|
|
| 339 |
|
// calculate the forces |
| 340 |
|
|
| 341 |
< |
myFF->doForces(calcPot, 0); |
| 341 |
> |
myFF->doForces(calcPot, calcStress); |
| 342 |
|
|
| 343 |
|
// move b |
| 344 |
|
|
| 403 |
|
} |
| 404 |
|
} |
| 405 |
|
|
| 406 |
+ |
|
| 407 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
| 408 |
+ |
myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() ); |
| 409 |
+ |
|
| 410 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) { |
| 411 |
+ |
tStats->getPressureTensor(press); |
| 412 |
+ |
myES->NoseHooverAndersonNPT( dt, |
| 413 |
+ |
tStats->getKinetic(), |
| 414 |
+ |
press); |
| 415 |
+ |
} |
| 416 |
+ |
|
| 417 |
|
time = tl + 1; |
| 418 |
|
|
| 419 |
|
if( entry_plug->setTemp ){ |
| 420 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 421 |
|
} |
| 422 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 423 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
| 424 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 423 |
> |
if( !((time+1) % status_n) ) { |
| 424 |
> |
calcPot = 1; |
| 425 |
> |
calcStress = 1; |
| 426 |
> |
} |
| 427 |
> |
if( !(time % status_n) ){ |
| 428 |
> |
e_out->writeStat( time * dt ); |
| 429 |
> |
calcPot = 0; |
| 430 |
> |
if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1; |
| 431 |
> |
else calcStress = 0; |
| 432 |
> |
} |
| 433 |
|
} |
| 434 |
|
} |
| 435 |
|
else{ |
| 439 |
|
kE = 0.0; |
| 440 |
|
rot_kE= 0.0; |
| 441 |
|
trans_kE = 0.0; |
| 442 |
+ |
|
| 443 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
| 444 |
+ |
myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() ); |
| 445 |
|
|
| 446 |
|
for( i=0; i<nAtoms; i++ ){ |
| 447 |
|
|
| 530 |
|
|
| 531 |
|
// calculate the forces |
| 532 |
|
|
| 533 |
< |
myFF->doForces(calcPot,0); |
| 533 |
> |
myFF->doForces(calcPot,calcStress); |
| 534 |
|
|
| 535 |
|
for( i=0; i< nAtoms; i++ ){ |
| 536 |
|
|
| 547 |
|
atoms[i]->set_vy( vy ); |
| 548 |
|
atoms[i]->set_vz( vz ); |
| 549 |
|
|
| 550 |
< |
// vx2 = vx * vx; |
| 551 |
< |
// vy2 = vy * vy; |
| 552 |
< |
// vz2 = vz * vz; |
| 550 |
> |
vx2 = vx * vx; |
| 551 |
> |
vy2 = vy * vy; |
| 552 |
> |
vz2 = vz * vz; |
| 553 |
|
|
| 554 |
|
if( atoms[i]->isDirectional() ){ |
| 555 |
|
|
| 581 |
|
dAtom->setJy( ji[1] ); |
| 582 |
|
dAtom->setJz( ji[2] ); |
| 583 |
|
} |
| 584 |
+ |
|
| 585 |
|
} |
| 586 |
< |
|
| 586 |
> |
|
| 587 |
> |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
| 588 |
> |
myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() ); |
| 589 |
> |
|
| 590 |
> |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) { |
| 591 |
> |
tStats->getPressureTensor(press); |
| 592 |
> |
myES->NoseHooverAndersonNPT( dt, |
| 593 |
> |
tStats->getKinetic(), |
| 594 |
> |
press); |
| 595 |
> |
} |
| 596 |
> |
|
| 597 |
|
time = tl + 1; |
| 598 |
|
|
| 599 |
|
if( entry_plug->setTemp ){ |
| 600 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 601 |
|
} |
| 602 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 603 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
| 604 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 603 |
> |
if( !((time+1) % status_n) ) { |
| 604 |
> |
calcPot = 1; |
| 605 |
> |
calcStress = 1; |
| 606 |
> |
} |
| 607 |
> |
if( !(time % status_n) ){ |
| 608 |
> |
e_out->writeStat( time * dt ); |
| 609 |
> |
calcPot = 0; |
| 610 |
> |
if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1; |
| 611 |
> |
else calcStress = 0; |
| 612 |
> |
} |
| 613 |
|
} |
| 614 |
|
} |
| 615 |
|
|
