ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/Torsion.cpp
Revision: 11
Committed: Tue Jul 9 18:40:59 2002 UTC (21 years, 11 months ago) by mmeineke
File size: 5624 byte(s)
Log Message: 
This commit was generated by cvs2svn to compensate for changes in r10, which
included commits to RCS files with non-trunk default branches.

File Contents

# Content
1 #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 }