| 160 |
|
double thermalTime = info->thermalTime; |
| 161 |
|
double resetTime = info->resetTime; |
| 162 |
|
|
| 163 |
< |
|
| 163 |
> |
double difference; |
| 164 |
|
double currSample; |
| 165 |
|
double currThermal; |
| 166 |
|
double currStatus; |
| 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); |
| 215 |
|
#endif // is_mpi |
| 216 |
|
|
| 217 |
|
while (info->getTime() < runTime && !stopIntegrator()){ |
| 218 |
< |
if ((info->getTime() + dt) >= currStatus){ |
| 218 |
> |
difference = info->getTime() + dt - currStatus; |
| 219 |
> |
if (difference > 0 || fabs(difference) < 1e-4 ){ |
| 220 |
|
calcPot = 1; |
| 221 |
|
calcStress = 1; |
| 222 |
|
} |
| 692 |
|
|
| 693 |
|
double angle; |
| 694 |
|
double A[3][3], I[3][3]; |
| 695 |
+ |
int i, j, k; |
| 696 |
|
|
| 697 |
|
// use the angular velocities to propagate the rotation matrix a |
| 698 |
|
// full time step |
| 700 |
|
sd->getA(A); |
| 701 |
|
sd->getI(I); |
| 702 |
|
|
| 703 |
< |
// rotate about the x-axis |
| 704 |
< |
angle = dt2 * ji[0] / I[0][0]; |
| 705 |
< |
this->rotate( 1, 2, angle, ji, A ); |
| 703 |
> |
if (sd->isLinear()) { |
| 704 |
> |
i = sd->linearAxis(); |
| 705 |
> |
j = (i+1)%3; |
| 706 |
> |
k = (i+2)%3; |
| 707 |
> |
|
| 708 |
> |
angle = dt2 * ji[j] / I[j][j]; |
| 709 |
> |
this->rotate( k, i, angle, ji, A ); |
| 710 |
|
|
| 711 |
< |
// rotate about the y-axis |
| 712 |
< |
angle = dt2 * ji[1] / I[1][1]; |
| 703 |
< |
this->rotate( 2, 0, angle, ji, A ); |
| 711 |
> |
angle = dt * ji[k] / I[k][k]; |
| 712 |
> |
this->rotate( i, j, angle, ji, A); |
| 713 |
|
|
| 714 |
< |
// rotate about the z-axis |
| 715 |
< |
angle = dt * ji[2] / I[2][2]; |
| 707 |
< |
this->rotate( 0, 1, angle, ji, A); |
| 714 |
> |
angle = dt2 * ji[j] / I[j][j]; |
| 715 |
> |
this->rotate( k, i, angle, ji, A ); |
| 716 |
|
|
| 717 |
< |
// rotate about the y-axis |
| 718 |
< |
angle = dt2 * ji[1] / I[1][1]; |
| 719 |
< |
this->rotate( 2, 0, angle, ji, A ); |
| 720 |
< |
|
| 721 |
< |
// rotate about the x-axis |
| 722 |
< |
angle = dt2 * ji[0] / I[0][0]; |
| 723 |
< |
this->rotate( 1, 2, angle, ji, A ); |
| 724 |
< |
|
| 717 |
> |
} else { |
| 718 |
> |
// rotate about the x-axis |
| 719 |
> |
angle = dt2 * ji[0] / I[0][0]; |
| 720 |
> |
this->rotate( 1, 2, angle, ji, A ); |
| 721 |
> |
|
| 722 |
> |
// rotate about the y-axis |
| 723 |
> |
angle = dt2 * ji[1] / I[1][1]; |
| 724 |
> |
this->rotate( 2, 0, angle, ji, A ); |
| 725 |
> |
|
| 726 |
> |
// rotate about the z-axis |
| 727 |
> |
angle = dt * ji[2] / I[2][2]; |
| 728 |
> |
this->rotate( 0, 1, angle, ji, A); |
| 729 |
> |
|
| 730 |
> |
// rotate about the y-axis |
| 731 |
> |
angle = dt2 * ji[1] / I[1][1]; |
| 732 |
> |
this->rotate( 2, 0, angle, ji, A ); |
| 733 |
> |
|
| 734 |
> |
// rotate about the x-axis |
| 735 |
> |
angle = dt2 * ji[0] / I[0][0]; |
| 736 |
> |
this->rotate( 1, 2, angle, ji, A ); |
| 737 |
> |
|
| 738 |
> |
} |
| 739 |
|
sd->setA( A ); |
| 740 |
|
} |
| 741 |
|
|
