31 |
|
} |
32 |
|
|
33 |
|
nAtoms = info->n_atoms; |
34 |
+ |
integrableObjects = info->integrableObjects; |
35 |
|
|
36 |
|
// check for constraints |
37 |
|
|
179 |
|
|
180 |
|
readyCheck(); |
181 |
|
|
182 |
+ |
// remove center of mass drift velocity (in case we passed in a configuration |
183 |
+ |
// that was drifting |
184 |
+ |
tStats->removeCOMdrift(); |
185 |
+ |
|
186 |
|
// initialize the forces before the first step |
187 |
|
|
188 |
|
calcForce(1, 1); |
184 |
– |
|
185 |
– |
//temp test |
186 |
– |
tStats->getPotential(); |
189 |
|
|
190 |
|
if (nConstrained){ |
191 |
|
preMove(); |
214 |
|
MPIcheckPoint(); |
215 |
|
#endif // is_mpi |
216 |
|
|
217 |
< |
while (info->getTime() < runTime){ |
217 |
> |
while (info->getTime() < runTime && !stopIntegrator()){ |
218 |
|
if ((info->getTime() + dt) >= currStatus){ |
219 |
|
calcPot = 1; |
220 |
|
calcStress = 1; |
336 |
|
|
337 |
|
|
338 |
|
template<typename T> void Integrator<T>::moveA(void){ |
339 |
< |
int i, j; |
339 |
> |
size_t i, j; |
340 |
|
DirectionalAtom* dAtom; |
341 |
|
double Tb[3], ji[3]; |
342 |
|
double vel[3], pos[3], frc[3]; |
343 |
|
double mass; |
344 |
< |
|
345 |
< |
for (i = 0; i < nAtoms; i++){ |
346 |
< |
atoms[i]->getVel(vel); |
347 |
< |
atoms[i]->getPos(pos); |
348 |
< |
atoms[i]->getFrc(frc); |
349 |
< |
|
350 |
< |
mass = atoms[i]->getMass(); |
344 |
> |
|
345 |
> |
for (i = 0; i < integrableObjects.size() ; i++){ |
346 |
> |
integrableObjects[i]->getVel(vel); |
347 |
> |
integrableObjects[i]->getPos(pos); |
348 |
> |
integrableObjects[i]->getFrc(frc); |
349 |
> |
|
350 |
> |
mass = integrableObjects[i]->getMass(); |
351 |
|
|
352 |
|
for (j = 0; j < 3; j++){ |
353 |
|
// velocity half step |
356 |
|
pos[j] += dt * vel[j]; |
357 |
|
} |
358 |
|
|
359 |
< |
atoms[i]->setVel(vel); |
360 |
< |
atoms[i]->setPos(pos); |
359 |
> |
integrableObjects[i]->setVel(vel); |
360 |
> |
integrableObjects[i]->setPos(pos); |
361 |
|
|
362 |
< |
if (atoms[i]->isDirectional()){ |
361 |
< |
dAtom = (DirectionalAtom *) atoms[i]; |
362 |
> |
if (integrableObjects[i]->isDirectional()){ |
363 |
|
|
364 |
|
// get and convert the torque to body frame |
365 |
|
|
366 |
< |
dAtom->getTrq(Tb); |
367 |
< |
dAtom->lab2Body(Tb); |
366 |
> |
integrableObjects[i]->getTrq(Tb); |
367 |
> |
integrableObjects[i]->lab2Body(Tb); |
368 |
|
|
369 |
|
// get the angular momentum, and propagate a half step |
370 |
|
|
371 |
< |
dAtom->getJ(ji); |
371 |
> |
integrableObjects[i]->getJ(ji); |
372 |
|
|
373 |
|
for (j = 0; j < 3; j++) |
374 |
|
ji[j] += (dt2 * Tb[j]) * eConvert; |
375 |
|
|
376 |
< |
this->rotationPropagation( dAtom, ji ); |
376 |
> |
this->rotationPropagation( integrableObjects[i], ji ); |
377 |
|
|
378 |
< |
dAtom->setJ(ji); |
378 |
> |
integrableObjects[i]->setJ(ji); |
379 |
|
} |
380 |
|
} |
381 |
|
|
387 |
|
|
388 |
|
template<typename T> void Integrator<T>::moveB(void){ |
389 |
|
int i, j; |
389 |
– |
DirectionalAtom* dAtom; |
390 |
|
double Tb[3], ji[3]; |
391 |
|
double vel[3], frc[3]; |
392 |
|
double mass; |
393 |
|
|
394 |
< |
for (i = 0; i < nAtoms; i++){ |
395 |
< |
atoms[i]->getVel(vel); |
396 |
< |
atoms[i]->getFrc(frc); |
394 |
> |
for (i = 0; i < integrableObjects.size(); i++){ |
395 |
> |
integrableObjects[i]->getVel(vel); |
396 |
> |
integrableObjects[i]->getFrc(frc); |
397 |
|
|
398 |
< |
mass = atoms[i]->getMass(); |
398 |
> |
mass = integrableObjects[i]->getMass(); |
399 |
|
|
400 |
|
// velocity half step |
401 |
|
for (j = 0; j < 3; j++) |
402 |
|
vel[j] += (dt2 * frc[j] / mass) * eConvert; |
403 |
|
|
404 |
< |
atoms[i]->setVel(vel); |
404 |
> |
integrableObjects[i]->setVel(vel); |
405 |
|
|
406 |
< |
if (atoms[i]->isDirectional()){ |
407 |
< |
dAtom = (DirectionalAtom *) atoms[i]; |
406 |
> |
if (integrableObjects[i]->isDirectional()){ |
407 |
|
|
408 |
|
// get and convert the torque to body frame |
409 |
|
|
410 |
< |
dAtom->getTrq(Tb); |
411 |
< |
dAtom->lab2Body(Tb); |
410 |
> |
integrableObjects[i]->getTrq(Tb); |
411 |
> |
integrableObjects[i]->lab2Body(Tb); |
412 |
|
|
413 |
|
// get the angular momentum, and propagate a half step |
414 |
|
|
415 |
< |
dAtom->getJ(ji); |
415 |
> |
integrableObjects[i]->getJ(ji); |
416 |
|
|
417 |
|
for (j = 0; j < 3; j++) |
418 |
|
ji[j] += (dt2 * Tb[j]) * eConvert; |
419 |
|
|
420 |
|
|
421 |
< |
dAtom->setJ(ji); |
421 |
> |
integrableObjects[i]->setJ(ji); |
422 |
|
} |
423 |
|
} |
424 |
|
|
687 |
|
} |
688 |
|
|
689 |
|
template<typename T> void Integrator<T>::rotationPropagation |
690 |
< |
( DirectionalAtom* dAtom, double ji[3] ){ |
690 |
> |
( StuntDouble* sd, double ji[3] ){ |
691 |
|
|
692 |
|
double angle; |
693 |
|
double A[3][3], I[3][3]; |
694 |
+ |
int i, j, k; |
695 |
|
|
696 |
|
// use the angular velocities to propagate the rotation matrix a |
697 |
|
// full time step |
698 |
|
|
699 |
< |
dAtom->getA(A); |
700 |
< |
dAtom->getI(I); |
699 |
> |
sd->getA(A); |
700 |
> |
sd->getI(I); |
701 |
|
|
702 |
< |
// rotate about the x-axis |
703 |
< |
angle = dt2 * ji[0] / I[0][0]; |
704 |
< |
this->rotate( 1, 2, angle, ji, A ); |
702 |
> |
if (sd->isLinear()) { |
703 |
> |
i = sd->linearAxis(); |
704 |
> |
j = (i+1)%3; |
705 |
> |
k = (i+2)%3; |
706 |
> |
|
707 |
> |
angle = dt2 * ji[j] / I[j][j]; |
708 |
> |
this->rotate( k, i, angle, ji, A ); |
709 |
|
|
710 |
< |
// rotate about the y-axis |
711 |
< |
angle = dt2 * ji[1] / I[1][1]; |
708 |
< |
this->rotate( 2, 0, angle, ji, A ); |
710 |
> |
angle = dt * ji[k] / I[k][k]; |
711 |
> |
this->rotate( i, j, angle, ji, A); |
712 |
|
|
713 |
< |
// rotate about the z-axis |
714 |
< |
angle = dt * ji[2] / I[2][2]; |
712 |
< |
this->rotate( 0, 1, angle, ji, A); |
713 |
> |
angle = dt2 * ji[j] / I[j][j]; |
714 |
> |
this->rotate( k, i, angle, ji, A ); |
715 |
|
|
716 |
< |
// rotate about the y-axis |
717 |
< |
angle = dt2 * ji[1] / I[1][1]; |
718 |
< |
this->rotate( 2, 0, angle, ji, A ); |
719 |
< |
|
720 |
< |
// rotate about the x-axis |
721 |
< |
angle = dt2 * ji[0] / I[0][0]; |
722 |
< |
this->rotate( 1, 2, angle, ji, A ); |
723 |
< |
|
724 |
< |
dAtom->setA( A ); |
716 |
> |
} else { |
717 |
> |
// rotate about the x-axis |
718 |
> |
angle = dt2 * ji[0] / I[0][0]; |
719 |
> |
this->rotate( 1, 2, angle, ji, A ); |
720 |
> |
|
721 |
> |
// rotate about the y-axis |
722 |
> |
angle = dt2 * ji[1] / I[1][1]; |
723 |
> |
this->rotate( 2, 0, angle, ji, A ); |
724 |
> |
|
725 |
> |
// rotate about the z-axis |
726 |
> |
angle = dt * ji[2] / I[2][2]; |
727 |
> |
this->rotate( 0, 1, angle, ji, A); |
728 |
> |
|
729 |
> |
// rotate about the y-axis |
730 |
> |
angle = dt2 * ji[1] / I[1][1]; |
731 |
> |
this->rotate( 2, 0, angle, ji, A ); |
732 |
> |
|
733 |
> |
// rotate about the x-axis |
734 |
> |
angle = dt2 * ji[0] / I[0][0]; |
735 |
> |
this->rotate( 1, 2, angle, ji, A ); |
736 |
> |
|
737 |
> |
} |
738 |
> |
sd->setA( A ); |
739 |
|
} |
740 |
|
|
741 |
|
template<typename T> void Integrator<T>::rotate(int axes1, int axes2, |