# | Line 44 | Line 44 | namespace oopse { | |
---|---|---|
44 | namespace oopse { | |
45 | ||
46 | /**@todo still a lot left to improve*/ | |
47 | < | void Bend::calcForce() { |
47 | > | void Bend::calcForce(RealType& angle) { |
48 | Vector3d pos1 = atom1_->getPos(); | |
49 | Vector3d pos2 = atom2_->getPos(); | |
50 | Vector3d pos3 = atom3_->getPos(); | |
51 | ||
52 | Vector3d r21 = pos1 - pos2; | |
53 | < | double d21 = r21.length(); |
53 | > | RealType d21 = r21.length(); |
54 | ||
55 | < | double d21inv = 1.0 / d21; |
55 | > | RealType d21inv = 1.0 / d21; |
56 | ||
57 | Vector3d r23 = pos3 - pos2; | |
58 | < | double d23 = r23.length(); |
58 | > | RealType d23 = r23.length(); |
59 | ||
60 | < | double d23inv = 1.0 / d23; |
60 | > | RealType d23inv = 1.0 / d23; |
61 | ||
62 | < | double cosTheta = dot(r21, r23) / (d21 * d23); |
62 | > | RealType cosTheta = dot(r21, r23) / (d21 * d23); |
63 | ||
64 | //check roundoff | |
65 | if (cosTheta > 1.0) { | |
# | Line 68 | Line 68 | namespace oopse { | |
68 | cosTheta = -1.0; | |
69 | } | |
70 | ||
71 | < | double theta = acos(cosTheta); |
71 | > | RealType theta = acos(cosTheta); |
72 | ||
73 | < | double dVdTheta; |
73 | > | RealType dVdTheta; |
74 | ||
75 | bendType_->calcForce(theta, potential_, dVdTheta); | |
76 | + | //std::cout << atom1_->getType() << "\t" << atom2_->getType() << "\t" << atom3_->getType() << "\t"; |
77 | + | //std::cout << "theta = " << theta/M_PI * 180.0 <<", potential = " << potential_ << std::endl; |
78 | ||
79 | < | double sinTheta = sqrt(1.0 - cosTheta * cosTheta); |
79 | > | RealType sinTheta = sqrt(1.0 - cosTheta * cosTheta); |
80 | ||
81 | if (fabs(sinTheta) < 1.0E-6) { | |
82 | sinTheta = 1.0E-6; | |
83 | } | |
84 | ||
85 | < | double commonFactor1 = dVdTheta / sinTheta * d21inv; |
86 | < | double commonFactor2 = dVdTheta / sinTheta * d23inv; |
85 | > | RealType commonFactor1 = dVdTheta / sinTheta * d21inv; |
86 | > | RealType commonFactor2 = dVdTheta / sinTheta * d23inv; |
87 | ||
88 | Vector3d force1 = commonFactor1 * (r23 * d23inv - r21*d21inv*cosTheta); | |
89 | Vector3d force3 = commonFactor2 * (r21 * d21inv - r23*d23inv*cosTheta); | |
# | Line 93 | Line 95 | namespace oopse { | |
95 | atom1_->addFrc(force1); | |
96 | atom2_->addFrc(force2); | |
97 | atom3_->addFrc(force3); | |
98 | + | |
99 | + | angle = theta /M_PI * 180.0; |
100 | } | |
101 | ||
102 | } //end namespace oopse |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |