| 43 |
|
|
| 44 |
|
namespace oopse { |
| 45 |
|
|
| 46 |
< |
void DLM::doRotate(StuntDouble* sd, Vector3d& ji, double dt) { |
| 47 |
< |
double dt2 = 0.5 * dt; |
| 48 |
< |
double angle; |
| 46 |
> |
void DLM::doRotate(StuntDouble* sd, Vector3d& ji, RealType dt) { |
| 47 |
> |
RealType dt2 = 0.5 * dt; |
| 48 |
> |
RealType angle; |
| 49 |
|
|
| 50 |
|
RotMat3x3d A = sd->getA(); |
| 51 |
|
Mat3x3d I = sd->getI(); |
| 94 |
|
} |
| 95 |
|
|
| 96 |
|
|
| 97 |
< |
void DLM::rotateStep(int axes1, int axes2, double angle, Vector3d& ji, RotMat3x3d& A) { |
| 97 |
> |
void DLM::rotateStep(int axes1, int axes2, RealType angle, Vector3d& ji, RotMat3x3d& A) { |
| 98 |
|
|
| 99 |
< |
double sinAngle; |
| 100 |
< |
double cosAngle; |
| 101 |
< |
double angleSqr; |
| 102 |
< |
double angleSqrOver4; |
| 103 |
< |
double top, bottom; |
| 99 |
> |
RealType sinAngle; |
| 100 |
> |
RealType cosAngle; |
| 101 |
> |
RealType angleSqr; |
| 102 |
> |
RealType angleSqrOver4; |
| 103 |
> |
RealType top, bottom; |
| 104 |
|
|
| 105 |
|
RotMat3x3d tempA(A); // initialize the tempA |
| 106 |
|
Vector3d tempJ(0.0); |
| 109 |
|
|
| 110 |
|
// use a small angle aproximation for sin and cosine |
| 111 |
|
|
| 112 |
< |
//angleSqr = angle * angle; |
| 113 |
< |
//angleSqrOver4 = angleSqr / 4.0; |
| 114 |
< |
//top = 1.0 - angleSqrOver4; |
| 115 |
< |
//bottom = 1.0 + angleSqrOver4; |
| 112 |
> |
angleSqr = angle * angle; |
| 113 |
> |
angleSqrOver4 = angleSqr / 4.0; |
| 114 |
> |
top = 1.0 - angleSqrOver4; |
| 115 |
> |
bottom = 1.0 + angleSqrOver4; |
| 116 |
|
|
| 117 |
< |
//cosAngle = top / bottom; |
| 118 |
< |
//sinAngle = angle / bottom; |
| 119 |
< |
cosAngle = cos(angle); |
| 120 |
< |
sinAngle = sin(angle); |
| 117 |
> |
cosAngle = top / bottom; |
| 118 |
> |
sinAngle = angle / bottom; |
| 119 |
> |
|
| 120 |
> |
// or don't use the small angle approximation: |
| 121 |
> |
//cosAngle = cos(angle); |
| 122 |
> |
//sinAngle = sin(angle); |
| 123 |
|
rot(axes1, axes1) = cosAngle; |
| 124 |
|
rot(axes2, axes2) = cosAngle; |
| 125 |
|
|
| 129 |
|
// rotate the momentum acoording to: ji[] = rot[][] * ji[] |
| 130 |
|
ji = rot * ji; |
| 131 |
|
|
| 132 |
< |
// rotate the Rotation matrix acording to: |
| 133 |
< |
// A[][] = A[][] * transpose(rot[][]) |
| 134 |
< |
// transpose(A[][]) = transpose(A[][]) * transpose(rot[][]) |
| 132 |
> |
// This code comes from converting an algorithm detailed in |
| 133 |
> |
// J. Chem. Phys. 107 (15), pp. 5840-5851 by Dullweber, |
| 134 |
> |
// Leimkuhler and McLachlan (DLM) for use in our code. |
| 135 |
> |
// In Appendix A, the DLM paper has the change to the rotation |
| 136 |
> |
// matrix as: Q = Q * rot.transpose(), but our rotation matrix |
| 137 |
> |
// A is actually equivalent to Q.transpose(). This fact can be |
| 138 |
> |
// seen on page 5849 of the DLM paper where a lab frame |
| 139 |
> |
// dipole \mu_i(t) is expressed in terms of a body-fixed |
| 140 |
> |
// reference orientation, \bar{\mu_i} and the rotation matrix, Q: |
| 141 |
> |
// \mu_i(t) = Q * \bar{\mu_i} |
| 142 |
> |
// Our code computes lab frame vectors from body-fixed reference |
| 143 |
> |
// vectors using: |
| 144 |
> |
// v_{lab} = A.transpose() * v_{body} |
| 145 |
> |
// (See StuntDouble.hpp for confirmation of this fact). |
| 146 |
> |
// |
| 147 |
> |
// So, using the identity: |
| 148 |
> |
// (A * B).transpose() = B.transpose() * A.transpose(), we |
| 149 |
> |
// get the equivalent of Q = Q * rot.transpose() for our code to be: |
| 150 |
|
|
| 151 |
< |
A = rot * A; //? A = A* rot.transpose(); |
| 151 |
> |
A = rot * A; |
| 152 |
|
|
| 153 |
|
} |
| 154 |
|
|
