1 |
#include "primitives/SRI.hpp"
|
2 |
#include "primitives/Atom.hpp"
|
3 |
#include <math.h>
|
4 |
#include <iostream>
|
5 |
#include <stdlib.h>
|
6 |
|
7 |
void Torsion::set_atoms( Atom &a, Atom &b, Atom &c, Atom &d){
|
8 |
c_p_a = &a;
|
9 |
c_p_b = &b;
|
10 |
c_p_c = &c;
|
11 |
c_p_d = &d;
|
12 |
}
|
13 |
|
14 |
|
15 |
void Torsion::calc_forces(){
|
16 |
|
17 |
/**********************************************************************
|
18 |
*
|
19 |
* initialize vectors
|
20 |
*
|
21 |
***********************************************************************/
|
22 |
|
23 |
vect r_ab; /* the vector whose origin is a and end is b */
|
24 |
vect r_cb; /* the vector whose origin is c and end is b */
|
25 |
vect r_cd; /* the vector whose origin is c and end is b */
|
26 |
vect r_cr1; /* the cross product of r_ab and r_cb */
|
27 |
vect r_cr2; /* the cross product of r_cb and r_cd */
|
28 |
|
29 |
double r_cr1_x2; /* the components of r_cr1 squared */
|
30 |
double r_cr1_y2;
|
31 |
double r_cr1_z2;
|
32 |
|
33 |
double r_cr2_x2; /* the components of r_cr2 squared */
|
34 |
double r_cr2_y2;
|
35 |
double r_cr2_z2;
|
36 |
|
37 |
double r_cr1_sqr; /* the length of r_cr1 squared */
|
38 |
double r_cr2_sqr; /* the length of r_cr2 squared */
|
39 |
|
40 |
double r_cr1_r_cr2; /* the length of r_cr1 * length of r_cr2 */
|
41 |
|
42 |
Vector3d aR, bR, cR, dR;
|
43 |
Vector3d aF, bF, cF, dF;
|
44 |
|
45 |
aR = c_p_a->getPos();
|
46 |
bR = c_p_b->getPos();
|
47 |
cR = c_p_c->getPos();
|
48 |
dR = c_p_d->getPos();
|
49 |
|
50 |
r_ab.x = bR[0] - aR[0];
|
51 |
r_ab.y = bR[1] - aR[1];
|
52 |
r_ab.z = bR[2] - aR[2];
|
53 |
r_ab.length = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z));
|
54 |
|
55 |
r_cb.x = bR[0] - cR[0];
|
56 |
r_cb.y = bR[1] - cR[1];
|
57 |
r_cb.z = bR[2] - cR[2];
|
58 |
r_cb.length = sqrt((r_cb.x * r_cb.x + r_cb.y * r_cb.y + r_cb.z * r_cb.z));
|
59 |
|
60 |
r_cd.x = dR[0] - cR[0];
|
61 |
r_cd.y = dR[1] - cR[1];
|
62 |
r_cd.z = dR[2] - cR[2];
|
63 |
r_cd.length = sqrt((r_cd.x * r_cd.x + r_cd.y * r_cd.y + r_cd.z * r_cd.z));
|
64 |
|
65 |
r_cr1.x = r_ab.y * r_cb.z - r_cb.y * r_ab.z;
|
66 |
r_cr1.y = r_ab.z * r_cb.x - r_cb.z * r_ab.x;
|
67 |
r_cr1.z = r_ab.x * r_cb.y - r_cb.x * r_ab.y;
|
68 |
r_cr1_x2 = r_cr1.x * r_cr1.x;
|
69 |
r_cr1_y2 = r_cr1.y * r_cr1.y;
|
70 |
r_cr1_z2 = r_cr1.z * r_cr1.z;
|
71 |
r_cr1_sqr = r_cr1_x2 + r_cr1_y2 + r_cr1_z2;
|
72 |
r_cr1.length = sqrt(r_cr1_sqr);
|
73 |
|
74 |
r_cr2.x = r_cb.y * r_cd.z - r_cd.y * r_cb.z;
|
75 |
r_cr2.y = r_cb.z * r_cd.x - r_cd.z * r_cb.x;
|
76 |
r_cr2.z = r_cb.x * r_cd.y - r_cd.x * r_cb.y;
|
77 |
r_cr2_x2 = r_cr2.x * r_cr2.x;
|
78 |
r_cr2_y2 = r_cr2.y * r_cr2.y;
|
79 |
r_cr2_z2 = r_cr2.z * r_cr2.z;
|
80 |
r_cr2_sqr = r_cr2_x2 + r_cr2_y2 + r_cr2_z2;
|
81 |
r_cr2.length = sqrt(r_cr2_sqr);
|
82 |
|
83 |
r_cr1_r_cr2 = r_cr1.length * r_cr2.length;
|
84 |
|
85 |
/**********************************************************************
|
86 |
*
|
87 |
* dot product and angle calculations
|
88 |
*
|
89 |
***********************************************************************/
|
90 |
|
91 |
double cr1_dot_cr2; /* the dot product of the cr1 and cr2 vectors */
|
92 |
double cos_phi; /* the cosine of the torsion angle */
|
93 |
|
94 |
cr1_dot_cr2 = r_cr1.x * r_cr2.x + r_cr1.y * r_cr2.y + r_cr1.z * r_cr2.z;
|
95 |
|
96 |
cos_phi = cr1_dot_cr2 / r_cr1_r_cr2;
|
97 |
|
98 |
/* adjust for the granularity of the numbers for angles near 0 or pi */
|
99 |
|
100 |
if(cos_phi > 1.0) cos_phi = 1.0;
|
101 |
if(cos_phi < -1.0) cos_phi = -1.0;
|
102 |
|
103 |
|
104 |
/********************************************************************
|
105 |
*
|
106 |
* This next section calculates derivatives needed for the force
|
107 |
* calculation
|
108 |
*
|
109 |
********************************************************************/
|
110 |
|
111 |
|
112 |
/* the derivatives of cos phi with respect to the x, y,
|
113 |
and z components of vectors cr1 and cr2. */
|
114 |
double d_cos_dx_cr1;
|
115 |
double d_cos_dy_cr1;
|
116 |
double d_cos_dz_cr1;
|
117 |
double d_cos_dx_cr2;
|
118 |
double d_cos_dy_cr2;
|
119 |
double d_cos_dz_cr2;
|
120 |
|
121 |
d_cos_dx_cr1 = r_cr2.x / r_cr1_r_cr2 - (cos_phi * r_cr1.x) / r_cr1_sqr;
|
122 |
d_cos_dy_cr1 = r_cr2.y / r_cr1_r_cr2 - (cos_phi * r_cr1.y) / r_cr1_sqr;
|
123 |
d_cos_dz_cr1 = r_cr2.z / r_cr1_r_cr2 - (cos_phi * r_cr1.z) / r_cr1_sqr;
|
124 |
|
125 |
d_cos_dx_cr2 = r_cr1.x / r_cr1_r_cr2 - (cos_phi * r_cr2.x) / r_cr2_sqr;
|
126 |
d_cos_dy_cr2 = r_cr1.y / r_cr1_r_cr2 - (cos_phi * r_cr2.y) / r_cr2_sqr;
|
127 |
d_cos_dz_cr2 = r_cr1.z / r_cr1_r_cr2 - (cos_phi * r_cr2.z) / r_cr2_sqr;
|
128 |
|
129 |
/***********************************************************************
|
130 |
*
|
131 |
* Calculate the actual forces and place them in the atoms.
|
132 |
*
|
133 |
***********************************************************************/
|
134 |
|
135 |
double force; /*the force scaling factor */
|
136 |
|
137 |
force = torsion_force(cos_phi);
|
138 |
|
139 |
aF[0] = force * (d_cos_dy_cr1 * r_cb.z - d_cos_dz_cr1 * r_cb.y);
|
140 |
aF[1] = force * (d_cos_dz_cr1 * r_cb.x - d_cos_dx_cr1 * r_cb.z);
|
141 |
aF[2] = force * (d_cos_dx_cr1 * r_cb.y - d_cos_dy_cr1 * r_cb.x);
|
142 |
|
143 |
bF[0] = force * ( d_cos_dy_cr1 * (r_ab.z - r_cb.z)
|
144 |
- d_cos_dy_cr2 * r_cd.z
|
145 |
+ d_cos_dz_cr1 * (r_cb.y - r_ab.y)
|
146 |
+ d_cos_dz_cr2 * r_cd.y);
|
147 |
bF[1] = force * ( d_cos_dx_cr1 * (r_cb.z - r_ab.z)
|
148 |
+ d_cos_dx_cr2 * r_cd.z
|
149 |
+ d_cos_dz_cr1 * (r_ab.x - r_cb.x)
|
150 |
- d_cos_dz_cr2 * r_cd.x);
|
151 |
bF[2] = force * ( d_cos_dx_cr1 * (r_ab.y - r_cb.y)
|
152 |
- d_cos_dx_cr2 * r_cd.y
|
153 |
+ d_cos_dy_cr1 * (r_cb.x - r_ab.x)
|
154 |
+ d_cos_dy_cr2 * r_cd.x);
|
155 |
|
156 |
cF[0] = force * (- d_cos_dy_cr1 * r_ab.z
|
157 |
- d_cos_dy_cr2 * (r_cb.z - r_cd.z)
|
158 |
+ d_cos_dz_cr1 * r_ab.y
|
159 |
- d_cos_dz_cr2 * (r_cd.y - r_cb.y));
|
160 |
cF[1] = force * ( d_cos_dx_cr1 * r_ab.z
|
161 |
- d_cos_dx_cr2 * (r_cd.z - r_cb.z)
|
162 |
- d_cos_dz_cr1 * r_ab.x
|
163 |
- d_cos_dz_cr2 * (r_cb.x - r_cd.x));
|
164 |
cF[2] = force * (- d_cos_dx_cr1 * r_ab.y
|
165 |
- d_cos_dx_cr2 * (r_cb.y - r_cd.y)
|
166 |
+ d_cos_dy_cr1 * r_ab.x
|
167 |
- d_cos_dy_cr2 * (r_cd.x - r_cb.x));
|
168 |
|
169 |
dF[0] = force * (d_cos_dy_cr2 * r_cb.z - d_cos_dz_cr2 * r_cb.y);
|
170 |
dF[1] = force * (d_cos_dz_cr2 * r_cb.x - d_cos_dx_cr2 * r_cb.z);
|
171 |
dF[2] = force * (d_cos_dx_cr2 * r_cb.y - d_cos_dy_cr2 * r_cb.x);
|
172 |
|
173 |
|
174 |
c_p_a->addFrc(aF);
|
175 |
c_p_b->addFrc(bF);
|
176 |
c_p_c->addFrc(cF);
|
177 |
c_p_d->addFrc(dF);
|
178 |
}
|