| 31 |
|
} |
| 32 |
|
|
| 33 |
|
nAtoms = info->n_atoms; |
| 34 |
+ |
integrableObjects = info->integrableObjects; |
| 35 |
|
|
| 36 |
|
// check for constraints |
| 37 |
|
|
| 69 |
|
|
| 70 |
|
SRI** theArray; |
| 71 |
|
for (int i = 0; i < nMols; i++){ |
| 72 |
< |
theArray = (SRI * *) molecules[i].getMyBonds(); |
| 72 |
> |
|
| 73 |
> |
theArray = (SRI * *) molecules[i].getMyBonds(); |
| 74 |
|
for (int j = 0; j < molecules[i].getNBonds(); j++){ |
| 75 |
|
constrained = theArray[j]->is_constrained(); |
| 76 |
|
|
| 116 |
|
} |
| 117 |
|
} |
| 118 |
|
|
| 119 |
+ |
|
| 120 |
|
if (nConstrained > 0){ |
| 121 |
|
isConstrained = 1; |
| 122 |
|
|
| 138 |
|
} |
| 139 |
|
|
| 140 |
|
|
| 141 |
< |
// save oldAtoms to check for lode balanceing later on. |
| 141 |
> |
// save oldAtoms to check for lode balancing later on. |
| 142 |
|
|
| 143 |
|
oldAtoms = nAtoms; |
| 144 |
|
|
| 183 |
|
|
| 184 |
|
calcForce(1, 1); |
| 185 |
|
|
| 186 |
+ |
//temp test |
| 187 |
+ |
tStats->getPotential(); |
| 188 |
+ |
|
| 189 |
|
if (nConstrained){ |
| 190 |
|
preMove(); |
| 191 |
|
constrainA(); |
| 335 |
|
|
| 336 |
|
|
| 337 |
|
template<typename T> void Integrator<T>::moveA(void){ |
| 338 |
< |
int i, j; |
| 338 |
> |
size_t i, j; |
| 339 |
|
DirectionalAtom* dAtom; |
| 340 |
|
double Tb[3], ji[3]; |
| 341 |
|
double vel[3], pos[3], frc[3]; |
| 342 |
|
double mass; |
| 343 |
< |
|
| 344 |
< |
for (i = 0; i < nAtoms; i++){ |
| 345 |
< |
atoms[i]->getVel(vel); |
| 346 |
< |
atoms[i]->getPos(pos); |
| 347 |
< |
atoms[i]->getFrc(frc); |
| 348 |
< |
|
| 349 |
< |
mass = atoms[i]->getMass(); |
| 343 |
> |
|
| 344 |
> |
for (i = 0; i < integrableObjects.size() ; i++){ |
| 345 |
> |
integrableObjects[i]->getVel(vel); |
| 346 |
> |
integrableObjects[i]->getPos(pos); |
| 347 |
> |
integrableObjects[i]->getFrc(frc); |
| 348 |
> |
|
| 349 |
> |
mass = integrableObjects[i]->getMass(); |
| 350 |
|
|
| 351 |
|
for (j = 0; j < 3; j++){ |
| 352 |
|
// velocity half step |
| 355 |
|
pos[j] += dt * vel[j]; |
| 356 |
|
} |
| 357 |
|
|
| 358 |
< |
atoms[i]->setVel(vel); |
| 359 |
< |
atoms[i]->setPos(pos); |
| 358 |
> |
integrableObjects[i]->setVel(vel); |
| 359 |
> |
integrableObjects[i]->setPos(pos); |
| 360 |
|
|
| 361 |
< |
if (atoms[i]->isDirectional()){ |
| 356 |
< |
dAtom = (DirectionalAtom *) atoms[i]; |
| 361 |
> |
if (integrableObjects[i]->isDirectional()){ |
| 362 |
|
|
| 363 |
|
// get and convert the torque to body frame |
| 364 |
|
|
| 365 |
< |
dAtom->getTrq(Tb); |
| 366 |
< |
dAtom->lab2Body(Tb); |
| 365 |
> |
integrableObjects[i]->getTrq(Tb); |
| 366 |
> |
integrableObjects[i]->lab2Body(Tb); |
| 367 |
|
|
| 368 |
|
// get the angular momentum, and propagate a half step |
| 369 |
|
|
| 370 |
< |
dAtom->getJ(ji); |
| 370 |
> |
integrableObjects[i]->getJ(ji); |
| 371 |
|
|
| 372 |
|
for (j = 0; j < 3; j++) |
| 373 |
|
ji[j] += (dt2 * Tb[j]) * eConvert; |
| 374 |
|
|
| 375 |
< |
this->rotationPropagation( dAtom, ji ); |
| 375 |
> |
this->rotationPropagation( integrableObjects[i], ji ); |
| 376 |
|
|
| 377 |
< |
dAtom->setJ(ji); |
| 377 |
> |
integrableObjects[i]->setJ(ji); |
| 378 |
|
} |
| 379 |
|
} |
| 380 |
|
|
| 386 |
|
|
| 387 |
|
template<typename T> void Integrator<T>::moveB(void){ |
| 388 |
|
int i, j; |
| 384 |
– |
DirectionalAtom* dAtom; |
| 389 |
|
double Tb[3], ji[3]; |
| 390 |
|
double vel[3], frc[3]; |
| 391 |
|
double mass; |
| 392 |
|
|
| 393 |
< |
for (i = 0; i < nAtoms; i++){ |
| 394 |
< |
atoms[i]->getVel(vel); |
| 395 |
< |
atoms[i]->getFrc(frc); |
| 393 |
> |
for (i = 0; i < integrableObjects.size(); i++){ |
| 394 |
> |
integrableObjects[i]->getVel(vel); |
| 395 |
> |
integrableObjects[i]->getFrc(frc); |
| 396 |
|
|
| 397 |
< |
mass = atoms[i]->getMass(); |
| 397 |
> |
mass = integrableObjects[i]->getMass(); |
| 398 |
|
|
| 399 |
|
// velocity half step |
| 400 |
|
for (j = 0; j < 3; j++) |
| 401 |
|
vel[j] += (dt2 * frc[j] / mass) * eConvert; |
| 402 |
|
|
| 403 |
< |
atoms[i]->setVel(vel); |
| 403 |
> |
integrableObjects[i]->setVel(vel); |
| 404 |
|
|
| 405 |
< |
if (atoms[i]->isDirectional()){ |
| 402 |
< |
dAtom = (DirectionalAtom *) atoms[i]; |
| 405 |
> |
if (integrableObjects[i]->isDirectional()){ |
| 406 |
|
|
| 407 |
|
// get and convert the torque to body frame |
| 408 |
|
|
| 409 |
< |
dAtom->getTrq(Tb); |
| 410 |
< |
dAtom->lab2Body(Tb); |
| 409 |
> |
integrableObjects[i]->getTrq(Tb); |
| 410 |
> |
integrableObjects[i]->lab2Body(Tb); |
| 411 |
|
|
| 412 |
|
// get the angular momentum, and propagate a half step |
| 413 |
|
|
| 414 |
< |
dAtom->getJ(ji); |
| 414 |
> |
integrableObjects[i]->getJ(ji); |
| 415 |
|
|
| 416 |
|
for (j = 0; j < 3; j++) |
| 417 |
|
ji[j] += (dt2 * Tb[j]) * eConvert; |
| 418 |
|
|
| 419 |
|
|
| 420 |
< |
dAtom->setJ(ji); |
| 420 |
> |
integrableObjects[i]->setJ(ji); |
| 421 |
|
} |
| 422 |
|
} |
| 423 |
|
|
| 686 |
|
} |
| 687 |
|
|
| 688 |
|
template<typename T> void Integrator<T>::rotationPropagation |
| 689 |
< |
( DirectionalAtom* dAtom, double ji[3] ){ |
| 689 |
> |
( StuntDouble* sd, double ji[3] ){ |
| 690 |
|
|
| 691 |
|
double angle; |
| 692 |
|
double A[3][3], I[3][3]; |
| 694 |
|
// use the angular velocities to propagate the rotation matrix a |
| 695 |
|
// full time step |
| 696 |
|
|
| 697 |
< |
dAtom->getA(A); |
| 698 |
< |
dAtom->getI(I); |
| 697 |
> |
sd->getA(A); |
| 698 |
> |
sd->getI(I); |
| 699 |
|
|
| 700 |
|
// rotate about the x-axis |
| 701 |
|
angle = dt2 * ji[0] / I[0][0]; |
| 717 |
|
angle = dt2 * ji[0] / I[0][0]; |
| 718 |
|
this->rotate( 1, 2, angle, ji, A ); |
| 719 |
|
|
| 720 |
< |
dAtom->setA( A ); |
| 720 |
> |
sd->setA( A ); |
| 721 |
|
} |
| 722 |
|
|
| 723 |
|
template<typename T> void Integrator<T>::rotate(int axes1, int axes2, |
