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

Comparing branches/new_design/OOPSE-2.0/src/primitives/Bend.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 8 | Line 8 | void Bend::calcForce() {
8      Vector3d pos2 = atom2_->getPos();
9      Vector3d pos3 = atom3_->getPos();
10  
11 <    Vector3d r12 = pos1 - pos2;
12 <    double d12 = r12.length();
11 >    Vector3d r21 = pos1 - pos2;
12 >    double d21 = r21.length();
13  
14 <    double d12inv = 1.0 / d12;
14 >    double d21inv = 1.0 / d21;
15  
16 <    Vector3d r32 = pos3 - pos2;
17 <    double d32 = r32.length();
16 >    Vector3d r23 = pos3 - pos2;
17 >    double d23 = r23.length();
18  
19 <    double d32inv = 1.0 / d32;
19 >    double d23inv = 1.0 / d23;
20  
21 <    double cosTheta = dot(r12, r32) / (d12 * d32);
21 >    double cosTheta = dot(r21, r23) / (d21 * d23);
22  
23      //check roundoff    
24      if (cosTheta > 1.0) {
# Line 29 | Line 29 | void Bend::calcForce() {
29  
30      double theta = acos(cosTheta);
31  
32 <    double firstDerivative;
32 >    double dVdTheta;
33  
34 <    bendType_->calcForce(theta, firstDerivative, potential_);
34 >    bendType_->calcForce(theta, potential_, dVdTheta);
35  
36      double sinTheta = sqrt(1.0 - cosTheta * cosTheta);
37  
38 <    if (fabs(sinTheta) < 1.0E-12) {
39 <        sinTheta = 1.0E-12;
38 >    if (fabs(sinTheta) < 1.0E-6) {
39 >        sinTheta = 1.0E-6;
40      }
41  
42 <    double commonFactor1 = -firstDerivative / sinTheta * d12inv;
43 <    double commonFactor2 = -firstDerivative / sinTheta * d32inv;
42 >    double commonFactor1 = dVdTheta / sinTheta * d21inv;
43 >    double commonFactor2 = dVdTheta / sinTheta * d23inv;
44  
45 <    Vector3d force1 = commonFactor1*(r12*(d12inv*cosTheta) - r32*d32inv);
45 >    Vector3d force1 = commonFactor1 * (r23 * d23inv - r21*d21inv*cosTheta);
46 >    Vector3d force3 = commonFactor2 * (r21 * d21inv - r23*d23inv*cosTheta);
47  
47 –    Vector3d force3 = commonFactor2*(r32*(d32inv*cosTheta) - r12*d12inv);
48 –
48      //total force in current bend is zero
49      Vector3d force2 = force1 + force3;
50      force2 *= -1.0;

Diff Legend

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