| 1 |
tim |
1692 |
#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 |
|
|
double aR[3], bR[3], cR[3], dR[3];
|
| 43 |
|
|
double aF[3], bF[3], cF[3], dF[3];
|
| 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 |
|
|
}
|
