1 |
mmeineke |
10 |
#include "SRI.hpp" |
2 |
|
|
#include "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 |
|
|
r_ab.x = c_p_b->getX() - c_p_a->getX(); |
43 |
|
|
r_ab.y = c_p_b->getY() - c_p_a->getY(); |
44 |
|
|
r_ab.z = c_p_b->getZ() - c_p_a->getZ(); |
45 |
|
|
r_ab.length = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z)); |
46 |
|
|
|
47 |
|
|
r_cb.x = c_p_b->getX() - c_p_c->getX(); |
48 |
|
|
r_cb.y = c_p_b->getY() - c_p_c->getY(); |
49 |
|
|
r_cb.z = c_p_b->getZ() - c_p_c->getZ(); |
50 |
|
|
r_cb.length = sqrt((r_cb.x * r_cb.x + r_cb.y * r_cb.y + r_cb.z * r_cb.z)); |
51 |
|
|
|
52 |
|
|
r_cd.x = c_p_d->getX() - c_p_c->getX(); |
53 |
|
|
r_cd.y = c_p_d->getY() - c_p_c->getY(); |
54 |
|
|
r_cd.z = c_p_d->getZ() - c_p_c->getZ(); |
55 |
|
|
r_cd.length = sqrt((r_cd.x * r_cd.x + r_cd.y * r_cd.y + r_cd.z * r_cd.z)); |
56 |
|
|
|
57 |
|
|
r_cr1.x = r_ab.y * r_cb.z - r_cb.y * r_ab.z; |
58 |
|
|
r_cr1.y = r_ab.z * r_cb.x - r_cb.z * r_ab.x; |
59 |
|
|
r_cr1.z = r_ab.x * r_cb.y - r_cb.x * r_ab.y; |
60 |
|
|
r_cr1_x2 = r_cr1.x * r_cr1.x; |
61 |
|
|
r_cr1_y2 = r_cr1.y * r_cr1.y; |
62 |
|
|
r_cr1_z2 = r_cr1.z * r_cr1.z; |
63 |
|
|
r_cr1_sqr = r_cr1_x2 + r_cr1_y2 + r_cr1_z2; |
64 |
|
|
r_cr1.length = sqrt(r_cr1_sqr); |
65 |
|
|
|
66 |
|
|
r_cr2.x = r_cb.y * r_cd.z - r_cd.y * r_cb.z; |
67 |
|
|
r_cr2.y = r_cb.z * r_cd.x - r_cd.z * r_cb.x; |
68 |
|
|
r_cr2.z = r_cb.x * r_cd.y - r_cd.x * r_cb.y; |
69 |
|
|
r_cr2_x2 = r_cr2.x * r_cr2.x; |
70 |
|
|
r_cr2_y2 = r_cr2.y * r_cr2.y; |
71 |
|
|
r_cr2_z2 = r_cr2.z * r_cr2.z; |
72 |
|
|
r_cr2_sqr = r_cr2_x2 + r_cr2_y2 + r_cr2_z2; |
73 |
|
|
r_cr2.length = sqrt(r_cr2_sqr); |
74 |
|
|
|
75 |
|
|
r_cr1_r_cr2 = r_cr1.length * r_cr2.length; |
76 |
|
|
|
77 |
|
|
/********************************************************************** |
78 |
|
|
* |
79 |
|
|
* dot product and angle calculations |
80 |
|
|
* |
81 |
|
|
***********************************************************************/ |
82 |
|
|
|
83 |
|
|
double cr1_dot_cr2; /* the dot product of the cr1 and cr2 vectors */ |
84 |
|
|
double cos_phi; /* the cosine of the torsion angle */ |
85 |
|
|
|
86 |
|
|
cr1_dot_cr2 = r_cr1.x * r_cr2.x + r_cr1.y * r_cr2.y + r_cr1.z * r_cr2.z; |
87 |
|
|
|
88 |
|
|
cos_phi = cr1_dot_cr2 / r_cr1_r_cr2; |
89 |
|
|
|
90 |
|
|
/* adjust for the granularity of the numbers for angles near 0 or pi */ |
91 |
|
|
|
92 |
|
|
if(cos_phi > 1.0) cos_phi = 1.0; |
93 |
|
|
if(cos_phi < -1.0) cos_phi = -1.0; |
94 |
|
|
|
95 |
|
|
|
96 |
|
|
/******************************************************************** |
97 |
|
|
* |
98 |
|
|
* This next section calculates derivatives needed for the force |
99 |
|
|
* calculation |
100 |
|
|
* |
101 |
|
|
********************************************************************/ |
102 |
|
|
|
103 |
|
|
|
104 |
|
|
/* the derivatives of cos phi with respect to the x, y, |
105 |
|
|
and z components of vectors cr1 and cr2. */ |
106 |
|
|
double d_cos_dx_cr1; |
107 |
|
|
double d_cos_dy_cr1; |
108 |
|
|
double d_cos_dz_cr1; |
109 |
|
|
double d_cos_dx_cr2; |
110 |
|
|
double d_cos_dy_cr2; |
111 |
|
|
double d_cos_dz_cr2; |
112 |
|
|
|
113 |
|
|
d_cos_dx_cr1 = r_cr2.x / r_cr1_r_cr2 - (cos_phi * r_cr1.x) / r_cr1_sqr; |
114 |
|
|
d_cos_dy_cr1 = r_cr2.y / r_cr1_r_cr2 - (cos_phi * r_cr1.y) / r_cr1_sqr; |
115 |
|
|
d_cos_dz_cr1 = r_cr2.z / r_cr1_r_cr2 - (cos_phi * r_cr1.z) / r_cr1_sqr; |
116 |
|
|
|
117 |
|
|
d_cos_dx_cr2 = r_cr1.x / r_cr1_r_cr2 - (cos_phi * r_cr2.x) / r_cr2_sqr; |
118 |
|
|
d_cos_dy_cr2 = r_cr1.y / r_cr1_r_cr2 - (cos_phi * r_cr2.y) / r_cr2_sqr; |
119 |
|
|
d_cos_dz_cr2 = r_cr1.z / r_cr1_r_cr2 - (cos_phi * r_cr2.z) / r_cr2_sqr; |
120 |
|
|
|
121 |
|
|
/*********************************************************************** |
122 |
|
|
* |
123 |
|
|
* Calculate the actual forces and place them in the atoms. |
124 |
|
|
* |
125 |
|
|
***********************************************************************/ |
126 |
|
|
|
127 |
|
|
double force; /*the force scaling factor */ |
128 |
|
|
|
129 |
|
|
force = torsion_force(cos_phi); |
130 |
|
|
|
131 |
|
|
c_p_a->addFx(force * (d_cos_dy_cr1 * r_cb.z - d_cos_dz_cr1 * r_cb.y)); |
132 |
|
|
c_p_a->addFy(force * (d_cos_dz_cr1 * r_cb.x - d_cos_dx_cr1 * r_cb.z)); |
133 |
|
|
c_p_a->addFz(force * (d_cos_dx_cr1 * r_cb.y - d_cos_dy_cr1 * r_cb.x)); |
134 |
|
|
|
135 |
|
|
c_p_b->addFx(force * ( d_cos_dy_cr1 * (r_ab.z - r_cb.z) |
136 |
|
|
- d_cos_dy_cr2 * r_cd.z |
137 |
|
|
+ d_cos_dz_cr1 * (r_cb.y - r_ab.y) |
138 |
|
|
+ d_cos_dz_cr2 * r_cd.y)); |
139 |
|
|
c_p_b->addFy(force * ( d_cos_dx_cr1 * (r_cb.z - r_ab.z) |
140 |
|
|
+ d_cos_dx_cr2 * r_cd.z |
141 |
|
|
+ d_cos_dz_cr1 * (r_ab.x - r_cb.x) |
142 |
|
|
- d_cos_dz_cr2 * r_cd.x)); |
143 |
|
|
c_p_b->addFz(force * ( d_cos_dx_cr1 * (r_ab.y - r_cb.y) |
144 |
|
|
- d_cos_dx_cr2 * r_cd.y |
145 |
|
|
+ d_cos_dy_cr1 * (r_cb.x - r_ab.x) |
146 |
|
|
+ d_cos_dy_cr2 * r_cd.x)); |
147 |
|
|
|
148 |
|
|
c_p_c->addFx(force * (- d_cos_dy_cr1 * r_ab.z |
149 |
|
|
- d_cos_dy_cr2 * (r_cb.z - r_cd.z) |
150 |
|
|
+ d_cos_dz_cr1 * r_ab.y |
151 |
|
|
- d_cos_dz_cr2 * (r_cd.y - r_cb.y))); |
152 |
|
|
c_p_c->addFy(force * ( d_cos_dx_cr1 * r_ab.z |
153 |
|
|
- d_cos_dx_cr2 * (r_cd.z - r_cb.z) |
154 |
|
|
- d_cos_dz_cr1 * r_ab.x |
155 |
|
|
- d_cos_dz_cr2 * (r_cb.x - r_cd.x))); |
156 |
|
|
c_p_c->addFz(force * (- d_cos_dx_cr1 * r_ab.y |
157 |
|
|
- d_cos_dx_cr2 * (r_cb.y - r_cd.y) |
158 |
|
|
+ d_cos_dy_cr1 * r_ab.x |
159 |
|
|
- d_cos_dy_cr2 * (r_cd.x - r_cb.x))); |
160 |
|
|
|
161 |
|
|
c_p_d->addFx(force * (d_cos_dy_cr2 * r_cb.z - d_cos_dz_cr2 * r_cb.y)); |
162 |
|
|
c_p_d->addFy(force * (d_cos_dz_cr2 * r_cb.x - d_cos_dx_cr2 * r_cb.z)); |
163 |
|
|
c_p_d->addFz(force * (d_cos_dx_cr2 * r_cb.y - d_cos_dy_cr2 * r_cb.x)); |
164 |
|
|
} |