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 1745 by tim, Tue Nov 16 20:36:18 2004 UTC vs.
Revision 1746 by tim, Wed Nov 17 06:37:56 2004 UTC

# Line 1 | Line 1
1 < #include "primitives/Bend.hpp"
2 < namespace oopse {
3 <
4 < /**@todo still a lot left to improve*/
5 < void Bend::calcForce() {
6 <
7 <    Vector3d pos1 = atom1_->getPos();
8 <    Vector3d pos2 = atom2_->getPos();
9 <    Vector3d pos3 = atom3_->getPos();
10 <
11 <    Vector3d r12 = pos1 - pos2;
12 <    double d12 = r12.length();
13 <    double d12inv = 1. 0 / d12;
14 <    
15 <    Vector3d r32 = pos3 - pos2;
16 <    double d32 = r32.length();
17 <    double d32inv = 1. 0 / d32;
18 <    
19 <    double cosTheta = (r12*r32)/(d12*d32);
20 <
21 <    //check roundoff    
22 <    if (cosTheta > 1.0) {
23 <        cosTheta = 1.0;
24 <    } else if (cosTheta < -1.0) {
25 <        cosTheta = -1.0;    
26 <    }
27 <
28 <    double theta = acos(cosTheta);
29 <
30 <    double firstDerivative;
31 <    
32 <    bendType_->calcForce(theta, firstDerivative, potential_);
33 <
34 <    double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
35 <
36 <    if(fabs(sinTheta) < 1.0E-12 ) {
37 <        sinTheta = 1.0E-12;
38 <    }
39 <    
40 <    double  commonFactor1 = - firstDerivative / sinTheta * d12inv;
41 <    double commonFactor2 = - firstDerivative / sinTheta * d32inv;
42 <
43 <    Vector3d force1 = commonFactor1*(r12*(d12inv*cosTheta) - r32*d32inv);
44 <    Vector3d force3 = commonFactor2*(r32*(d32inv*cosTheta) - r12*d12inv);
45 <
46 <    //total force in current bend is zero
47 <    Vector3d force2 = force1 +  force3;
48 <    force2 *= -1.0;
49 <
50 <    atom1_->addFrc( force1);
51 <    atom2_->addFrc(force2);
52 <    atom3_->addFrc(force3);
53 <
54 < }
55 <
56 <    double k = value->k * scale;
57 <    double theta0 = value->theta0;
58 <
59 <    double diff = theta - theta0;
60 <
61 <    double energy = k *diff*diff;
62 <
63 <    double firstDerivative = 2.0 * k * diff;
64 <    
65 < }//end namespace oopse
1 > #include "primitives/Bend.hpp"
2 >
3 > namespace oopse {
4 >
5 > /**@todo still a lot left to improve*/
6 > void Bend::calcForce() {
7 >    Vector3d pos1 = atom1_->getPos();
8 >    Vector3d pos2 = atom2_->getPos();
9 >    Vector3d pos3 = atom3_->getPos();
10 >
11 >    Vector3d r12 = pos1 - pos2;
12 >    double d12 = r12.length();
13 >
14 >    double d12inv = 1. 0 / d12;
15 >
16 >    Vector3d r32 = pos3 - pos2;
17 >    double d32 = r32.length();
18 >
19 >    double d32inv = 1. 0 / d32;
20 >
21 >    double cosTheta = (r12 * r32) / (d12 * d32);
22 >
23 >    //check roundoff    
24 >    if (cosTheta > 1.0) {
25 >        cosTheta = 1.0;
26 >    } else if (cosTheta < -1.0) {
27 >        cosTheta = -1.0;
28 >    }
29 >
30 >    double theta = acos(cosTheta);
31 >
32 >    double firstDerivative;
33 >
34 >    bendType_->calcForce(theta, firstDerivative, potential_);
35 >
36 >    double sinTheta = sqrt(1.0 - cosTheta * cosTheta);
37 >
38 >    if (fabs(sinTheta) < 1.0E - 12) {
39 >        sinTheta = 1.0E - 12;
40 >    }
41 >
42 >    double commonFactor1 = -firstDerivative / sinTheta * d12inv;
43 >    double commonFactor2 = -firstDerivative / sinTheta * d32inv;
44 >
45 >    Vector3d force1 = commonFactor1*(r12*(d12inv*cosTheta) - r32*d32inv);
46 >
47 >    Vector3d force3 = commonFactor2*(r32*(d32inv*cosTheta) - r12*d12inv);
48 >
49 >    //total force in current bend is zero
50 >    Vector3d force2 = force1 + force3;
51 >    force2 *= -1.0;
52 >
53 >    atom1_->addFrc(force1);
54 >    atom2_->addFrc(force2);
55 >    atom3_->addFrc(force3);
56 > }
57 >
58 >    double k = value->k * scale;
59 >    double theta0 = value->theta0;
60 >
61 >    double diff = theta - theta0;
62 >
63 >    double energy = k * diff * diff;
64 >
65 >    double firstDerivative = 2.0 * k * diff;
66 >
67 > } //end namespace oopse

Diff Legend

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