ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/Torsion.cpp
Revision: 10
Committed: Tue Jul 9 18:40:59 2002 UTC (22 years ago) by mmeineke
Original Path: branches/mmeineke/mdtools/md_code/Torsion.cpp
File size: 5624 byte(s)
Log Message: 
everything you need to make libmdtools

File Contents

# User Rev Content
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     }