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
|