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 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 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 | Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3
43  
44   namespace oopse {
45  
46 < Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
47 <                 TorsionType *tt) :
46 >  Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
47 >                   TorsionType *tt) :
48      atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
49  
50 < void Torsion::calcForce() {
50 >  void Torsion::calcForce() {
51      Vector3d pos1 = atom1_->getPos();
52      Vector3d pos2 = atom2_->getPos();
53      Vector3d pos3 = atom3_->getPos();
# Line 85 | Line 85 | void Torsion::calcForce() {
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
88 >    //  use the sin version to avoid 1/cos terms
89  
90 <        Vector3d dcosdA = (cos_phi * A - B) /rA;
91 <        Vector3d dcosdB = (cos_phi * B - A) /rB;
90 >    Vector3d dcosdA = (cos_phi * A - B) /rA;
91 >    Vector3d dcosdB = (cos_phi * B - A) /rB;
92  
93 <        double dVdcosPhi = -dVdPhi / sin_phi;
93 >    double dVdcosPhi = -dVdPhi / sin_phi;
94  
95 <        f1 = dVdcosPhi * cross(r32, dcosdA);
96 <        f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
97 <        f3 = dVdcosPhi * cross(dcosdB, r32);
95 >    f1 = dVdcosPhi * cross(r32, dcosdA);
96 >    f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
97 >    f3 = dVdcosPhi * cross(dcosdB, r32);
98  
99      /** @todo fix below block, must be something wrong with the sign somewhere */
100      //} 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
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
103  
104 <        //double dVdsinPhi = dVdPhi /cos_phi;
105 <        //Vector3d dsindB = (sin_phi * B - C) /rB;
106 <        //Vector3d dsindC = (sin_phi * C - B) /rC;
104 >    //double dVdsinPhi = dVdPhi /cos_phi;
105 >    //Vector3d dsindB = (sin_phi * B - C) /rB;
106 >    //Vector3d dsindC = (sin_phi * C - B) /rC;
107  
108 <        //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());
108 >    //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());
109  
110 <        //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());
110 >    //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());
111  
112 <        //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());
112 >    //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());
113  
114 <        //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()
115 <        //+ (2.0*r32.x()*r21.z() - r21.x()*r32.z())*dsindC.z() + dsindB.z()*r43.y() - dsindB.y()*r43.z());
114 >    //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()
115 >    //+ (2.0*r32.x()*r21.z() - r21.x()*r32.z())*dsindC.z() + dsindB.z()*r43.y() - dsindB.y()*r43.z());
116  
117 <        //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()
118 <        //+ (2.0*r32.y()*r21.x() - r21.y()*r32.x())*dsindC.x() + dsindB.x()*r43.z() - dsindB.z()*r43.x());
117 >    //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()
118 >    //+ (2.0*r32.y()*r21.x() - r21.y()*r32.x())*dsindC.x() + dsindB.x()*r43.z() - dsindB.z()*r43.x());
119  
120 <        //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()
121 <        //+(2.0*r32.z()*r21.y() - r21.z()*r32.y())*dsindC.y() + dsindB.y()*r43.x() - dsindB.x()*r43.y());
120 >    //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()
121 >    //+(2.0*r32.z()*r21.y() - r21.z()*r32.y())*dsindC.y() + dsindB.y()*r43.x() - dsindB.x()*r43.y());
122  
123 <        //f3 = dVdsinPhi * cross(r32, dsindB);
123 >    //f3 = dVdsinPhi * cross(r32, dsindB);
124  
125      //}
126  
# Line 128 | Line 128 | void Torsion::calcForce() {
128      atom2_->addFrc(f2 - f1);
129      atom3_->addFrc(f3 - f2);
130      atom4_->addFrc(-f3);
131 < }
131 >  }
132  
133   }

Diff Legend

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