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

Comparing branches/new_design/OOPSE-3.0/src/primitives/Bend.cpp (file contents):
Revision 1782 by tim, Wed Nov 24 20:55:03 2004 UTC vs.
Revision 1849 by gezelter, Sat Dec 4 19:24:16 2004 UTC

# Line 1 | Line 1
1 + <<<<<<< Bend.cpp
2   #include "primitives/Bend.hpp"
3 + namespace oopse {
4  
5 + /**@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   namespace oopse {
71  
72   /**@todo still a lot left to improve*/
# Line 56 | Line 123 | void Bend::calcForce() {
123   }
124  
125   } //end namespace oopse
126 + >>>>>>> 1.2.2.4

Diff Legend

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