ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/integrators/DLM.cpp
(Generate patch)

Comparing trunk/src/integrators/DLM.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1339 by gezelter, Thu Apr 23 18:31:05 2009 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 43 | Line 43 | namespace oopse {
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();
# Line 53 | Line 53 | void DLM::doRotate(StuntDouble* sd, Vector3d& ji, doub
53      // use the angular velocities to propagate the rotation matrix a full time step
54      if (sd->isLinear()) {
55  
56 <        int i = sd->linearAxis();
57 <        int j = (i+1)%3;
58 <        int k = (i+2)%3;
56 >      int i = sd->linearAxis();
57 >      int j = (i+1)%3;
58 >      int k = (i+2)%3;
59  
60 <        angle = dt2 * ji[j] / I(j, j);
61 <        rotateStep( k, i, angle, ji, A );
60 >      angle = dt2 * ji[j] / I(j, j);
61 >      rotateStep( k, i, angle, ji, A );
62  
63 <        angle = dt * ji[k] / I(k, k);
64 <        rotateStep( i, j, angle, ji, A);
63 >      angle = dt * ji[k] / I(k, k);
64 >      rotateStep( i, j, angle, ji, A);
65  
66 <        angle = dt2 * ji[j] / I(j, j);
67 <        rotateStep( k, i, angle, ji, A );
66 >      angle = dt2 * ji[j] / I(j, j);
67 >      rotateStep( k, i, angle, ji, A );
68  
69      } else {
70 <        // rotate about the x-axis
71 <        angle = dt2 * ji[0] / I(0, 0);
72 <        rotateStep( 1, 2, angle, ji, A );
70 >      // rotate about the x-axis
71 >      angle = dt2 * ji[0] / I(0, 0);
72 >      rotateStep( 1, 2, angle, ji, A );
73  
74 <        // rotate about the y-axis
75 <        angle = dt2 * ji[1] / I(1, 1);
76 <        rotateStep( 2, 0, angle, ji, A );
74 >      // rotate about the y-axis
75 >      angle = dt2 * ji[1] / I(1, 1);
76 >      rotateStep( 2, 0, angle, ji, A );
77  
78 <        // rotate about the z-axis
79 <        angle = dt * ji[2] / I(2, 2);
80 <        sd->addZangle(angle);
81 <        rotateStep( 0, 1, angle, ji, A);
78 >      // rotate about the z-axis
79 >      angle = dt * ji[2] / I(2, 2);
80 >      sd->addZangle(angle);
81 >      rotateStep( 0, 1, angle, ji, A);
82  
83 <        // rotate about the y-axis
84 <        angle = dt2 * ji[1] / I(1, 1);
85 <        rotateStep( 2, 0, angle, ji, A );
83 >      // rotate about the y-axis
84 >      angle = dt2 * ji[1] / I(1, 1);
85 >      rotateStep( 2, 0, angle, ji, A );
86  
87 <        // rotate about the x-axis
88 <        angle = dt2 * ji[0] / I(0, 0);
89 <        rotateStep( 1, 2, angle, ji, A );
87 >      // rotate about the x-axis
88 >      angle = dt2 * ji[0] / I(0, 0);
89 >      rotateStep( 1, 2, angle, ji, A );
90  
91      }
92  
93      sd->setA( A  );
94 < }
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);
# Line 117 | Line 117 | void DLM::rotateStep(int axes1, int axes2, double angl
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  
# Line 126 | Line 129 | void DLM::rotateStep(int axes1, int axes2, double angl
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 < }
153 >  }
154  
155  
156   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines