ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/primitives/Torsion.cpp
Revision: 3
Committed: Fri Sep 24 16:27:58 2004 UTC (21 years, 1 month ago) by tim
File size: 5729 byte(s)
Log Message: 
change the #include in source files

File Contents

# User Rev Content
1 tim 3 #include "primitives/SRI.hpp"
2     #include "primitives/Atom.hpp"
3 gezelter 2 #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     }