| 31 |
|
} |
| 32 |
|
|
| 33 |
|
nAtoms = info->n_atoms; |
| 34 |
< |
|
| 34 |
> |
integrableObjects = info->integrableObjects; |
| 35 |
> |
|
| 36 |
|
// check for constraints |
| 37 |
|
|
| 38 |
|
constrainedA = NULL; |
| 45 |
|
nConstrained = 0; |
| 46 |
|
|
| 47 |
|
checkConstraints(); |
| 48 |
+ |
|
| 49 |
+ |
for (i=0; i<nMols; i++) |
| 50 |
+ |
zAngle[i] = 0.0; |
| 51 |
|
} |
| 52 |
|
|
| 53 |
|
template<typename T> Integrator<T>::~Integrator(){ |
| 72 |
|
|
| 73 |
|
SRI** theArray; |
| 74 |
|
for (int i = 0; i < nMols; i++){ |
| 75 |
< |
theArray = (SRI * *) molecules[i].getMyBonds(); |
| 75 |
> |
|
| 76 |
> |
theArray = (SRI * *) molecules[i].getMyBonds(); |
| 77 |
|
for (int j = 0; j < molecules[i].getNBonds(); j++){ |
| 78 |
|
constrained = theArray[j]->is_constrained(); |
| 79 |
|
|
| 119 |
|
} |
| 120 |
|
} |
| 121 |
|
|
| 122 |
+ |
|
| 123 |
|
if (nConstrained > 0){ |
| 124 |
|
isConstrained = 1; |
| 125 |
|
|
| 163 |
|
double thermalTime = info->thermalTime; |
| 164 |
|
double resetTime = info->resetTime; |
| 165 |
|
|
| 166 |
< |
|
| 166 |
> |
double difference; |
| 167 |
|
double currSample; |
| 168 |
|
double currThermal; |
| 169 |
|
double currStatus; |
| 181 |
|
dt2 = 0.5 * dt; |
| 182 |
|
|
| 183 |
|
readyCheck(); |
| 184 |
+ |
|
| 185 |
+ |
// remove center of mass drift velocity (in case we passed in a configuration |
| 186 |
+ |
// that was drifting |
| 187 |
+ |
tStats->removeCOMdrift(); |
| 188 |
+ |
|
| 189 |
+ |
// initialize the retraints if necessary |
| 190 |
+ |
if (info->useThermInt) { |
| 191 |
+ |
myFF->initRestraints(); |
| 192 |
+ |
} |
| 193 |
|
|
| 194 |
|
// initialize the forces before the first step |
| 195 |
|
|
| 196 |
|
calcForce(1, 1); |
| 182 |
– |
|
| 183 |
– |
//temp test |
| 184 |
– |
tStats->getPotential(); |
| 197 |
|
|
| 198 |
|
if (nConstrained){ |
| 199 |
|
preMove(); |
| 222 |
|
MPIcheckPoint(); |
| 223 |
|
#endif // is_mpi |
| 224 |
|
|
| 225 |
< |
while (info->getTime() < runTime){ |
| 226 |
< |
if ((info->getTime() + dt) >= currStatus){ |
| 225 |
> |
while (info->getTime() < runTime && !stopIntegrator()){ |
| 226 |
> |
difference = info->getTime() + dt - currStatus; |
| 227 |
> |
if (difference > 0 || fabs(difference) < 1e-4 ){ |
| 228 |
|
calcPot = 1; |
| 229 |
|
calcStress = 1; |
| 230 |
|
} |
| 257 |
|
|
| 258 |
|
if (info->getTime() >= currStatus){ |
| 259 |
|
statOut->writeStat(info->getTime()); |
| 260 |
+ |
statOut->writeRaw(info->getTime()); |
| 261 |
|
calcPot = 0; |
| 262 |
|
calcStress = 0; |
| 263 |
|
currStatus += statusTime; |
| 280 |
|
#endif // is_mpi |
| 281 |
|
} |
| 282 |
|
|
| 283 |
+ |
// dump out a file containing the omega values for the final configuration |
| 284 |
+ |
if (info->useThermInt) |
| 285 |
+ |
myFF->dumpzAngle(); |
| 286 |
+ |
|
| 287 |
+ |
|
| 288 |
|
delete dumpOut; |
| 289 |
|
delete statOut; |
| 290 |
|
} |
| 351 |
|
|
| 352 |
|
|
| 353 |
|
template<typename T> void Integrator<T>::moveA(void){ |
| 354 |
< |
int i, j; |
| 354 |
> |
size_t i, j; |
| 355 |
|
DirectionalAtom* dAtom; |
| 356 |
|
double Tb[3], ji[3]; |
| 357 |
|
double vel[3], pos[3], frc[3]; |
| 358 |
|
double mass; |
| 359 |
< |
|
| 360 |
< |
for (i = 0; i < nAtoms; i++){ |
| 361 |
< |
atoms[i]->getVel(vel); |
| 362 |
< |
atoms[i]->getPos(pos); |
| 363 |
< |
atoms[i]->getFrc(frc); |
| 364 |
< |
|
| 365 |
< |
mass = atoms[i]->getMass(); |
| 359 |
> |
|
| 360 |
> |
for (i = 0; i < integrableObjects.size() ; i++){ |
| 361 |
> |
integrableObjects[i]->getVel(vel); |
| 362 |
> |
integrableObjects[i]->getPos(pos); |
| 363 |
> |
integrableObjects[i]->getFrc(frc); |
| 364 |
> |
|
| 365 |
> |
mass = integrableObjects[i]->getMass(); |
| 366 |
|
|
| 367 |
|
for (j = 0; j < 3; j++){ |
| 368 |
|
// velocity half step |
| 371 |
|
pos[j] += dt * vel[j]; |
| 372 |
|
} |
| 373 |
|
|
| 374 |
< |
atoms[i]->setVel(vel); |
| 375 |
< |
atoms[i]->setPos(pos); |
| 374 |
> |
integrableObjects[i]->setVel(vel); |
| 375 |
> |
integrableObjects[i]->setPos(pos); |
| 376 |
|
|
| 377 |
< |
if (atoms[i]->isDirectional()){ |
| 359 |
< |
dAtom = (DirectionalAtom *) atoms[i]; |
| 377 |
> |
if (integrableObjects[i]->isDirectional()){ |
| 378 |
|
|
| 379 |
|
// get and convert the torque to body frame |
| 380 |
|
|
| 381 |
< |
dAtom->getTrq(Tb); |
| 382 |
< |
dAtom->lab2Body(Tb); |
| 381 |
> |
integrableObjects[i]->getTrq(Tb); |
| 382 |
> |
integrableObjects[i]->lab2Body(Tb); |
| 383 |
|
|
| 384 |
|
// get the angular momentum, and propagate a half step |
| 385 |
|
|
| 386 |
< |
dAtom->getJ(ji); |
| 386 |
> |
integrableObjects[i]->getJ(ji); |
| 387 |
|
|
| 388 |
|
for (j = 0; j < 3; j++) |
| 389 |
|
ji[j] += (dt2 * Tb[j]) * eConvert; |
| 390 |
|
|
| 391 |
< |
this->rotationPropagation( dAtom, ji ); |
| 391 |
> |
this->rotationPropagation( integrableObjects[i], ji ); |
| 392 |
|
|
| 393 |
< |
dAtom->setJ(ji); |
| 393 |
> |
integrableObjects[i]->setJ(ji); |
| 394 |
|
} |
| 395 |
|
} |
| 396 |
|
|
| 402 |
|
|
| 403 |
|
template<typename T> void Integrator<T>::moveB(void){ |
| 404 |
|
int i, j; |
| 387 |
– |
DirectionalAtom* dAtom; |
| 405 |
|
double Tb[3], ji[3]; |
| 406 |
|
double vel[3], frc[3]; |
| 407 |
|
double mass; |
| 408 |
|
|
| 409 |
< |
for (i = 0; i < nAtoms; i++){ |
| 410 |
< |
atoms[i]->getVel(vel); |
| 411 |
< |
atoms[i]->getFrc(frc); |
| 409 |
> |
for (i = 0; i < integrableObjects.size(); i++){ |
| 410 |
> |
integrableObjects[i]->getVel(vel); |
| 411 |
> |
integrableObjects[i]->getFrc(frc); |
| 412 |
|
|
| 413 |
< |
mass = atoms[i]->getMass(); |
| 413 |
> |
mass = integrableObjects[i]->getMass(); |
| 414 |
|
|
| 415 |
|
// velocity half step |
| 416 |
|
for (j = 0; j < 3; j++) |
| 417 |
|
vel[j] += (dt2 * frc[j] / mass) * eConvert; |
| 418 |
|
|
| 419 |
< |
atoms[i]->setVel(vel); |
| 419 |
> |
integrableObjects[i]->setVel(vel); |
| 420 |
|
|
| 421 |
< |
if (atoms[i]->isDirectional()){ |
| 405 |
< |
dAtom = (DirectionalAtom *) atoms[i]; |
| 421 |
> |
if (integrableObjects[i]->isDirectional()){ |
| 422 |
|
|
| 423 |
|
// get and convert the torque to body frame |
| 424 |
|
|
| 425 |
< |
dAtom->getTrq(Tb); |
| 426 |
< |
dAtom->lab2Body(Tb); |
| 425 |
> |
integrableObjects[i]->getTrq(Tb); |
| 426 |
> |
integrableObjects[i]->lab2Body(Tb); |
| 427 |
|
|
| 428 |
|
// get the angular momentum, and propagate a half step |
| 429 |
|
|
| 430 |
< |
dAtom->getJ(ji); |
| 430 |
> |
integrableObjects[i]->getJ(ji); |
| 431 |
|
|
| 432 |
|
for (j = 0; j < 3; j++) |
| 433 |
|
ji[j] += (dt2 * Tb[j]) * eConvert; |
| 434 |
|
|
| 435 |
|
|
| 436 |
< |
dAtom->setJ(ji); |
| 436 |
> |
integrableObjects[i]->setJ(ji); |
| 437 |
|
} |
| 438 |
|
} |
| 439 |
|
|
| 702 |
|
} |
| 703 |
|
|
| 704 |
|
template<typename T> void Integrator<T>::rotationPropagation |
| 705 |
< |
( DirectionalAtom* dAtom, double ji[3] ){ |
| 705 |
> |
( StuntDouble* sd, double ji[3] ){ |
| 706 |
|
|
| 707 |
|
double angle; |
| 708 |
|
double A[3][3], I[3][3]; |
| 709 |
+ |
int i, j, k; |
| 710 |
|
|
| 711 |
|
// use the angular velocities to propagate the rotation matrix a |
| 712 |
|
// full time step |
| 713 |
|
|
| 714 |
< |
dAtom->getA(A); |
| 715 |
< |
dAtom->getI(I); |
| 714 |
> |
sd->getA(A); |
| 715 |
> |
sd->getI(I); |
| 716 |
|
|
| 717 |
< |
// rotate about the x-axis |
| 718 |
< |
angle = dt2 * ji[0] / I[0][0]; |
| 719 |
< |
this->rotate( 1, 2, angle, ji, A ); |
| 717 |
> |
if (sd->isLinear()) { |
| 718 |
> |
i = sd->linearAxis(); |
| 719 |
> |
j = (i+1)%3; |
| 720 |
> |
k = (i+2)%3; |
| 721 |
> |
|
| 722 |
> |
angle = dt2 * ji[j] / I[j][j]; |
| 723 |
> |
this->rotate( k, i, angle, ji, A ); |
| 724 |
|
|
| 725 |
< |
// rotate about the y-axis |
| 726 |
< |
angle = dt2 * ji[1] / I[1][1]; |
| 706 |
< |
this->rotate( 2, 0, angle, ji, A ); |
| 707 |
< |
|
| 708 |
< |
// rotate about the z-axis |
| 709 |
< |
angle = dt * ji[2] / I[2][2]; |
| 710 |
< |
this->rotate( 0, 1, angle, ji, A); |
| 725 |
> |
angle = dt * ji[k] / I[k][k]; |
| 726 |
> |
this->rotate( i, j, angle, ji, A); |
| 727 |
|
|
| 728 |
< |
// rotate about the y-axis |
| 729 |
< |
angle = dt2 * ji[1] / I[1][1]; |
| 714 |
< |
this->rotate( 2, 0, angle, ji, A ); |
| 728 |
> |
angle = dt2 * ji[j] / I[j][j]; |
| 729 |
> |
this->rotate( k, i, angle, ji, A ); |
| 730 |
|
|
| 731 |
< |
// rotate about the x-axis |
| 732 |
< |
angle = dt2 * ji[0] / I[0][0]; |
| 733 |
< |
this->rotate( 1, 2, angle, ji, A ); |
| 734 |
< |
|
| 735 |
< |
dAtom->setA( A ); |
| 731 |
> |
} else { |
| 732 |
> |
// rotate about the x-axis |
| 733 |
> |
angle = dt2 * ji[0] / I[0][0]; |
| 734 |
> |
this->rotate( 1, 2, angle, ji, A ); |
| 735 |
> |
|
| 736 |
> |
// rotate about the y-axis |
| 737 |
> |
angle = dt2 * ji[1] / I[1][1]; |
| 738 |
> |
this->rotate( 2, 0, angle, ji, A ); |
| 739 |
> |
|
| 740 |
> |
// rotate about the z-axis |
| 741 |
> |
angle = dt * ji[2] / I[2][2]; |
| 742 |
> |
this->rotate( 0, 1, angle, ji, A); |
| 743 |
> |
|
| 744 |
> |
// rotate about the y-axis |
| 745 |
> |
angle = dt2 * ji[1] / I[1][1]; |
| 746 |
> |
this->rotate( 2, 0, angle, ji, A ); |
| 747 |
> |
|
| 748 |
> |
// rotate about the x-axis |
| 749 |
> |
angle = dt2 * ji[0] / I[0][0]; |
| 750 |
> |
this->rotate( 1, 2, angle, ji, A ); |
| 751 |
> |
|
| 752 |
> |
} |
| 753 |
> |
sd->setA( A ); |
| 754 |
|
} |
| 755 |
|
|
| 756 |
|
template<typename T> void Integrator<T>::rotate(int axes1, int axes2, |
