ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/primitives/Torsion.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/primitives/Torsion.cpp (file contents):
Revision 1691, Thu Oct 28 22:34:02 2004 UTC vs.
Revision 1692 by tim, Mon Nov 1 20:15:58 2004 UTC

# Line 1 | Line 1
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 <  double aR[3], bR[3], cR[3], dR[3];
43 <  double aF[3], bF[3], cF[3], dF[3];
44 <
45 <  c_p_a->getPos( aR );
46 <  c_p_b->getPos( bR );
47 <  c_p_c->getPos( cR );
48 <  c_p_d->getPos( dR );
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 < }
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 >  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 > }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines