11 |
|
Vector3d r12 = pos1 - pos2; |
12 |
|
double d12 = r12.length(); |
13 |
|
|
14 |
< |
double d12inv = 1. 0 / d12; |
14 |
> |
double d12inv = 1.0 / d12; |
15 |
|
|
16 |
|
Vector3d r32 = pos3 - pos2; |
17 |
|
double d32 = r32.length(); |
18 |
|
|
19 |
< |
double d32inv = 1. 0 / d32; |
19 |
> |
double d32inv = 1.0 / d32; |
20 |
|
|
21 |
< |
double cosTheta = (r12 * r32) / (d12 * d32); |
21 |
> |
double cosTheta = dot(r12, r32) / (d12 * d32); |
22 |
|
|
23 |
|
//check roundoff |
24 |
|
if (cosTheta > 1.0) { |
35 |
|
|
36 |
|
double sinTheta = sqrt(1.0 - cosTheta * cosTheta); |
37 |
|
|
38 |
< |
if (fabs(sinTheta) < 1.0E - 12) { |
39 |
< |
sinTheta = 1.0E - 12; |
38 |
> |
if (fabs(sinTheta) < 1.0E-12) { |
39 |
> |
sinTheta = 1.0E-12; |
40 |
|
} |
41 |
|
|
42 |
|
double commonFactor1 = -firstDerivative / sinTheta * d12inv; |
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 |
– |
|
58 |
|
} //end namespace oopse |