ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/primitives/Bend.cpp
Revision: 1849
Committed: Sat Dec 4 19:24:16 2004 UTC (20 years, 7 months ago) by gezelter
File size: 2975 byte(s)
Log Message: 
What?

File Contents

# User Rev Content
1 gezelter 1849 <<<<<<< Bend.cpp
2 tim 1746 #include "primitives/Bend.hpp"
3 gezelter 1849 namespace oopse {
4 tim 1746
5 gezelter 1849 /**@todo still a lot left to improve*/
6     void Bend::calcForce() {
7    
8     Vector3d pos1 = atom1_->getPos();
9     Vector3d pos2 = atom2_->getPos();
10     Vector3d pos3 = atom3_->getPos();
11    
12     Vector3d r12 = pos1 - pos2;
13     double d12 = r12.length();
14     double d12inv = 1. 0 / d12;
15    
16     Vector3d r32 = pos3 - pos2;
17     double d32 = r32.length();
18     double d32inv = 1. 0 / d32;
19    
20     double cosTheta = (r12*r32)/(d12*d32);
21    
22     //check roundoff
23     if (cosTheta > 1.0) {
24     cosTheta = 1.0;
25     } else if (cosTheta < -1.0) {
26     cosTheta = -1.0;
27     }
28    
29     double theta = acos(cosTheta);
30    
31     double firstDerivative;
32    
33     bendType_->calcForce(theta, firstDerivative, potential_);
34    
35     double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
36    
37     if(fabs(sinTheta) < 1.0E-12 ) {
38     sinTheta = 1.0E-12;
39     }
40    
41     double commonFactor1 = - firstDerivative / sinTheta * d12inv;
42     double commonFactor2 = - firstDerivative / sinTheta * d32inv;
43    
44     Vector3d force1 = commonFactor1*(r12*(d12inv*cosTheta) - r32*d32inv);
45     Vector3d force3 = commonFactor2*(r32*(d32inv*cosTheta) - r12*d12inv);
46    
47     //total force in current bend is zero
48     Vector3d force2 = force1 + force3;
49     force2 *= -1.0;
50    
51     atom1_->addFrc( force1);
52     atom2_->addFrc(force2);
53     atom3_->addFrc(force3);
54    
55     }
56    
57     double k = value->k * scale;
58     double theta0 = value->theta0;
59    
60     double diff = theta - theta0;
61    
62     double energy = k *diff*diff;
63    
64     double firstDerivative = 2.0 * k * diff;
65    
66     }//end namespace oopse
67     =======
68     #include "primitives/Bend.hpp"
69    
70 tim 1746 namespace oopse {
71    
72     /**@todo still a lot left to improve*/
73     void Bend::calcForce() {
74     Vector3d pos1 = atom1_->getPos();
75     Vector3d pos2 = atom2_->getPos();
76     Vector3d pos3 = atom3_->getPos();
77    
78     Vector3d r12 = pos1 - pos2;
79     double d12 = r12.length();
80    
81 tim 1782 double d12inv = 1.0 / d12;
82 tim 1746
83     Vector3d r32 = pos3 - pos2;
84     double d32 = r32.length();
85    
86 tim 1782 double d32inv = 1.0 / d32;
87 tim 1746
88 tim 1782 double cosTheta = dot(r12, r32) / (d12 * d32);
89 tim 1746
90     //check roundoff
91     if (cosTheta > 1.0) {
92     cosTheta = 1.0;
93     } else if (cosTheta < -1.0) {
94     cosTheta = -1.0;
95     }
96    
97     double theta = acos(cosTheta);
98    
99     double firstDerivative;
100    
101     bendType_->calcForce(theta, firstDerivative, potential_);
102    
103     double sinTheta = sqrt(1.0 - cosTheta * cosTheta);
104    
105 tim 1782 if (fabs(sinTheta) < 1.0E-12) {
106     sinTheta = 1.0E-12;
107 tim 1746 }
108    
109     double commonFactor1 = -firstDerivative / sinTheta * d12inv;
110     double commonFactor2 = -firstDerivative / sinTheta * d32inv;
111    
112     Vector3d force1 = commonFactor1*(r12*(d12inv*cosTheta) - r32*d32inv);
113    
114     Vector3d force3 = commonFactor2*(r32*(d32inv*cosTheta) - r12*d12inv);
115    
116     //total force in current bend is zero
117     Vector3d force2 = force1 + force3;
118     force2 *= -1.0;
119    
120     atom1_->addFrc(force1);
121     atom2_->addFrc(force2);
122     atom3_->addFrc(force3);
123     }
124    
125     } //end namespace oopse
126 gezelter 1849 >>>>>>> 1.2.2.4