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 1849 by gezelter, Sat Dec 4 19:24:16 2004 UTC vs.
Revision 1850 by tim, Sat Dec 4 19:29:54 2004 UTC

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

Diff Legend

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