ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/primitives/Torsion.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/primitives/Torsion.cpp (file contents):
Revision 1856 by tim, Sat Dec 4 19:29:54 2004 UTC vs.
Revision 1857 by tim, Mon Dec 6 20:15:02 2004 UTC

# Line 12 | Line 12 | void Torsion::calcForce() {
12      Vector3d pos3 = atom3_->getPos();
13      Vector3d pos4 = atom4_->getPos();
14  
15 <    Vector3d r12 = pos1 - pos2;
16 <    Vector3d r23 = pos2 - pos3;
17 <    Vector3d r34 = pos3 - pos4;
15 >    Vector3d r21 = pos1 - pos2;
16 >    Vector3d r32 = pos2 - pos3;
17 >    Vector3d r43 = pos3 - pos4;
18  
19      //  Calculate the cross products and distances
20 <    Vector3d A = cross(r12, r23);
20 >    Vector3d A = cross(r21, r32);
21      double rA = A.length();
22 <    Vector3d B = cross(r23, r34);
22 >    Vector3d B = cross(r32, r43);
23      double rB = B.length();
24 <    Vector3d C = cross(r23, A);
24 >    Vector3d C = cross(r32, A);
25      double rC = C.length();
26  
27      A.normalize();
# Line 49 | Line 49 | void Torsion::calcForce() {
49          Vector3d dcosdA = (cos_phi * A - B) /rA;
50          Vector3d dcosdB = (cos_phi * B - A) /rB;
51  
52 <        double dVdcosPhi = dVdPhi / sin_phi;
52 >        double dVdcosPhi = -dVdPhi / sin_phi;
53  
54 <        f1 = dVdcosPhi * cross(r23, dcosdA);
55 <        f2 = dVdcosPhi * ( cross(r34, dcosdB) - cross(r12, dcosdA));
56 <        f3 = dVdcosPhi * cross(r23, dcosdB);
54 >        f1 = dVdcosPhi * cross(r32, dcosdA);
55 >        f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
56 >        f3 = dVdcosPhi * cross(r32, dcosdB);
57  
58      } else {
59          //  This angle is closer to 0 or 180 than it is to
60          //  90, so use the cos version to avoid 1/sin terms
61  
62 <        double dVdsinPhi = -dVdPhi /cos_phi;
62 >        double dVdsinPhi = dVdPhi /cos_phi;
63          Vector3d dsindB = (sin_phi * B - C) /rB;
64          Vector3d dsindC = (sin_phi * C - B) /rC;
65  
66 <        f1.x() = dVdsinPhi*((r23.y()*r23.y() + r23.z()*r23.z())*dsindC.x() - r23.x()*r23.y()*dsindC.y() - r23.x()*r23.z()*dsindC.z());
66 >        f1.x() = dVdsinPhi*((r32.y()*r32.y() + r32.z()*r32.z())*dsindC.x() - r32.x()*r32.y()*dsindC.y() - r32.x()*r32.z()*dsindC.z());
67  
68 <        f1.y() = dVdsinPhi*((r23.z()*r23.z() + r23.x()*r23.x())*dsindC.y() - r23.y()*r23.z()*dsindC.z() - r23.y()*r23.x()*dsindC.x());
68 >        f1.y() = dVdsinPhi*((r32.z()*r32.z() + r32.x()*r32.x())*dsindC.y() - r32.y()*r32.z()*dsindC.z() - r32.y()*r32.x()*dsindC.x());
69  
70 <        f1.z() = dVdsinPhi*((r23.x()*r23.x() + r23.y()*r23.y())*dsindC.z() - r23.z()*r23.x()*dsindC.x() - r23.z()*r23.y()*dsindC.y());
70 >        f1.z() = dVdsinPhi*((r32.x()*r32.x() + r32.y()*r32.y())*dsindC.z() - r32.z()*r32.x()*dsindC.x() - r32.z()*r32.y()*dsindC.y());
71  
72 <        f2.x() = dVdsinPhi*(-(r23.y()*r12.y() + r23.z()*r12.z())*dsindC.x() + (2.0*r23.x()*r12.y() - r12.x()*r23.y())*dsindC.y()
73 <        + (2.0*r23.x()*r12.z() - r12.x()*r23.z())*dsindC.z() + dsindB.z()*r34.y() - dsindB.y()*r34.z());
72 >        f2.x() = dVdsinPhi*(-(r32.y()*r21.y() + r32.z()*r21.z())*dsindC.x() + (2.0*r32.x()*r21.y() - r21.x()*r32.y())*dsindC.y()
73 >        + (2.0*r32.x()*r21.z() - r21.x()*r32.z())*dsindC.z() + dsindB.z()*r43.y() - dsindB.y()*r43.z());
74  
75 <        f2.y() = dVdsinPhi*(-(r23.z()*r12.z() + r23.x()*r12.x())*dsindC.y() + (2.0*r23.y()*r12.z() - r12.y()*r23.z())*dsindC.z()
76 <        + (2.0*r23.y()*r12.x() - r12.y()*r23.x())*dsindC.x() + dsindB.x()*r34.z() - dsindB.z()*r34.x());
75 >        f2.y() = dVdsinPhi*(-(r32.z()*r21.z() + r32.x()*r21.x())*dsindC.y() + (2.0*r32.y()*r21.z() - r21.y()*r32.z())*dsindC.z()
76 >        + (2.0*r32.y()*r21.x() - r21.y()*r32.x())*dsindC.x() + dsindB.x()*r43.z() - dsindB.z()*r43.x());
77  
78 <        f2.z() = dVdsinPhi*(-(r23.x()*r12.x() + r23.y()*r12.y())*dsindC.z() + (2.0*r23.z()*r12.x() - r12.z()*r23.x())*dsindC.x()
79 <        +(2.0*r23.z()*r12.y() - r12.z()*r23.y())*dsindC.y() + dsindB.y()*r34.x() - dsindB.x()*r34.y());
78 >        f2.z() = dVdsinPhi*(-(r32.x()*r21.x() + r32.y()*r21.y())*dsindC.z() + (2.0*r32.z()*r21.x() - r21.z()*r32.x())*dsindC.x()
79 >        +(2.0*r32.z()*r21.y() - r21.z()*r32.y())*dsindC.y() + dsindB.y()*r43.x() - dsindB.x()*r43.y());
80  
81 <        f3 = dVdsinPhi * cross(dsindB, r23);
81 >        f3 = dVdsinPhi * cross(dsindB, r32);
82  
83      }
84  

Diff Legend

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