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

Comparing trunk/OOPSE-4/src/primitives/Torsion.cpp (file contents):
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 2391 by tim, Thu Oct 20 20:27:34 2005 UTC

# Line 80 | Line 80 | namespace oopse {
80      Vector3d f2;
81      Vector3d f3;
82  
83 <    //  Next, we want to calculate the forces.  In order
84 <    //  to do that, we first need to figure out whether the
85 <    //  sin or cos form will be more stable.  For this,
86 <    //  just look at the value of phi
87 <    //if (fabs(sin_phi) > 0.1) {
88 <    //  use the sin version to avoid 1/cos terms
83 >    if (fabs(sin_phi) > 0.5) {
84 >    //use the sin version to  prevent potential singularities
85  
86      Vector3d dcosdA = (cos_phi * A - B) /rA;
87      Vector3d dcosdB = (cos_phi * B - A) /rB;
# Line 96 | Line 92 | namespace oopse {
92      f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
93      f3 = dVdcosPhi * cross(dcosdB, r32);
94  
95 <    /** @todo fix below block, must be something wrong with the sign somewhere */
96 <    //} else {
101 <    //  This angle is closer to 0 or 180 than it is to
102 <    //  90, so use the cos version to avoid 1/sin terms
95 >    } else {
96 >    //use the cos version to  prevent potential singularities
97  
98 <    //double dVdsinPhi = dVdPhi /cos_phi;
99 <    //Vector3d dsindB = (sin_phi * B - C) /rB;
100 <    //Vector3d dsindC = (sin_phi * C - B) /rC;
98 >    double dVdsinPhi = dVdPhi /cos_phi;
99 >    Vector3d dsindB = (sin_phi * B - C) /rB;
100 >    Vector3d dsindC = (sin_phi * C - B) /rC;
101  
102 <    //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());
102 >    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());
103  
104 <    //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());
104 >    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());
105  
106 <    //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());
106 >    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());
107  
108 <    //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()
109 <    //+ (2.0*r32.x()*r21.z() - r21.x()*r32.z())*dsindC.z() + dsindB.z()*r43.y() - dsindB.y()*r43.z());
108 >    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()
109 >    + (2.0*r32.x()*r21.z() - r21.x()*r32.z())*dsindC.z() + dsindB.z()*r43.y() - dsindB.y()*r43.z());
110  
111 <    //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()
112 <    //+ (2.0*r32.y()*r21.x() - r21.y()*r32.x())*dsindC.x() + dsindB.x()*r43.z() - dsindB.z()*r43.x());
111 >    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()
112 >    + (2.0*r32.y()*r21.x() - r21.y()*r32.x())*dsindC.x() + dsindB.x()*r43.z() - dsindB.z()*r43.x());
113  
114 <    //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()
115 <    //+(2.0*r32.z()*r21.y() - r21.z()*r32.y())*dsindC.y() + dsindB.y()*r43.x() - dsindB.x()*r43.y());
114 >    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()
115 >    +(2.0*r32.z()*r21.y() - r21.z()*r32.y())*dsindC.y() + dsindB.y()*r43.x() - dsindB.x()*r43.y());
116  
117 <    //f3 = dVdsinPhi * cross(r32, dsindB);
117 >    f3 = dVdsinPhi * cross(dsindB, r32);
118 >    }
119  
125 –    //}
126 –
120      atom1_->addFrc(f1);
121      atom2_->addFrc(f2 - f1);
122      atom3_->addFrc(f3 - f2);

Diff Legend

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