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 1701 by tim, Wed Nov 3 16:08:43 2004 UTC vs.
Revision 1849 by gezelter, Sat Dec 4 19:24:16 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 <  Vector3d aR, bR, cR, dR;
43 <  Vector3d aF, bF, cF, dF;
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 < }
1 > <<<<<<< Torsion.cpp
2 > #include "primitives/Torsion.hpp"
3 >
4 > namespace oopse {
5 >
6 > Torsion::Torsion(Atom* atom1, Atom* atom2, Atom* atom3, Atom* atom4, TorsionType* tt)
7 >            : atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4) {
8 >
9 > }
10 >
11 > void Torsion::calcForce() {
12 >
13 >    Vector3d pos1 = atom1_->getPos();
14 >    Vector3d pos2 = atom2_->getPos();
15 >    Vector3d pos3 = atom3_->getPos();
16 >    Vector3d pos4 = atom4_->getPos();
17 >
18 >    Vector3d r12 = pos1 - pos2;
19 >    Vector3d r23 = pos2 - pos3;
20 >    Vector3d r34 = pos3 - pos4;
21 >
22 >    //  Calculate the cross products and distances
23 >    Vector3d A = cross(r12,r23);
24 >    double rA = A.length();
25 >    Vector3d B = cross(r23,r34);
26 >    double rB = B.length();
27 >    Vector3d C = cross(r23,A);
28 >    double rC = C.length();
29 >
30 >    //  Calculate the sin and cos
31 >    double cos_phi = (A*B)/(rA*rB);
32 >    double sin_phi = (C*B)/(rC*rB);
33 >
34 >    double phi= -atan2(sin_phi,cos_phi);
35 >
36 >    double firstDerivative;
37 >    
38 >    torsionType_->calcForce(cosPhi, sinPhi, Vtorsion, dVdCosPhi);
39 >
40 >
41 >    Vector3d f1,f2,f3;
42 >    
43 >    //  Normalize B
44 >    rB = 1.0/rB;
45 >    B *= rB;
46 >
47 >    //  Next, we want to calculate the forces.  In order
48 >    //  to do that, we first need to figure out whether the
49 >    //  sin or cos form will be more stable.  For this,
50 >    //  just look at the value of phi
51 >    if (fabs(sin_phi) > 0.1) {
52 >        //  use the sin version to avoid 1/cos terms
53 >
54 >        rA = 1.0/rA;
55 >        A *= rA;
56 >        Vector3d dcosdA = rA*(cos_phi*A-B);
57 >        Vector3d dcosdB = rB*(cos_phi*B-A);
58 >
59 >        K1 = K1/sin_phi;
60 >
61 >        //simple form
62 >        //f1 = K1 * cross(r23, dcosdA);
63 >        //f3 = K1 * cross(r23, dcosdB);
64 >        //f2 = K1 * ( cross(r34, dcosdB) - cross(r12, dcosdA));
65 >
66 >        f1.x = K1*(r23.y*dcosdA.z - r23.z*dcosdA.y);
67 >        f1.y = K1*(r23.z*dcosdA.x - r23.x*dcosdA.z);
68 >        f1.z = K1*(r23.x*dcosdA.y - r23.y*dcosdA.x);
69 >
70 >        f3.x = K1*(r23.z*dcosdB.y - r23.y*dcosdB.z);
71 >        f3.y = K1*(r23.x*dcosdB.z - r23.z*dcosdB.x);
72 >        f3.z = K1*(r23.y*dcosdB.x - r23.x*dcosdB.y);
73 >
74 >        f2.x = K1*(r12.z*dcosdA.y - r12.y*dcosdA.z + r34.y*dcosdB.z - r34.z*dcosdB.y);
75 >        f2.y = K1*(r12.x*dcosdA.z - r12.z*dcosdA.x + r34.z*dcosdB.x - r34.x*dcosdB.z);
76 >        f2.z = K1*(r12.y*dcosdA.x - r12.x*dcosdA.y + r34.x*dcosdB.y - r34.y*dcosdB.x);
77 >    } else {
78 >        //  This angle is closer to 0 or 180 than it is to
79 >        //  90, so use the cos version to avoid 1/sin terms
80 >
81 >        //  Normalize C
82 >        rC = 1.0/rC;
83 >        C *= rC;
84 >        Vector3d dsindC = rC*(sin_phi*C-B);
85 >        Vector3d dsindB = rB*(sin_phi*B-C);
86 >
87 >        K1 = -K1/cos_phi;
88 >
89 >        f1.x = K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x - r23.x*r23.y*dsindC.y - r23.x*r23.z*dsindC.z);
90 >        f1.y = K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y - r23.y*r23.z*dsindC.z - r23.y*r23.x*dsindC.x);
91 >        f1.z = K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z - r23.z*r23.x*dsindC.x - r23.z*r23.y*dsindC.y);
92 >
93 >        f3 = K1 *cross(dsindB,r23);
94 >
95 >        f2.x = K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x + (2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
96 >                 + (2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z + dsindB.z*r34.y - dsindB.y*r34.z);
97 >        f2.y = K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y + (2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
98 >                 + (2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x + dsindB.x*r34.z - dsindB.z*r34.x);
99 >        f2.z = K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z + (2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
100 >                  +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y + dsindB.y*r34.x - dsindB.x*r34.y);
101 >    }
102 >
103 >    atom1_->addFrc(f1);
104 >    atom2_->addFrc(f2 - f1);
105 >    atom3_->addFrc(f3 - f2);
106 >    atom4_->addFrc(-f3);
107 >    
108 > }
109 >
110 >
111 >      double K=0;       // energy
112 >      double K1=0;      // force
113 >
114 >  // get the dihedral information
115 >  int multiplicity = value->multiplicity;
116 >
117 >  //  Loop through the multiple parameter sets for this
118 >  //  bond.  We will only loop more than once if this
119 >  //  has multiple parameter sets from Charmm22
120 >  for (int mult_num=0; mult_num<multiplicity; mult_num++)
121 >  {
122 >    /* get angle information */
123 >    double k = value->values[mult_num].k * scale;
124 >    double delta = value->values[mult_num].delta;
125 >    int n = value->values[mult_num].n;
126 >
127 >    //  Calculate the energy
128 >    if (n)
129 >    {
130 >      //  Periodicity is greater than 0, so use cos form
131 >      K += k*(1+cos(n*phi + delta));
132 >      K1 += -n*k*sin(n*phi + delta);
133 >    }
134 >    else
135 >    {
136 >      //  Periodicity is 0, so just use the harmonic form
137 >      double diff = phi-delta;
138 >      if (diff < -PI)           diff += TWOPI;
139 >      else if (diff > PI)       diff -= TWOPI;
140 >
141 >      K += k*diff*diff;
142 >      K1 += 2.0*k*diff;
143 >    }
144 >  } /* for multiplicity */
145 >
146 >
147 > void Torsion::calc_forces(){
148 >  
149 >  /**********************************************************************
150 >   *
151 >   * initialize vectors
152 >   *
153 >   ***********************************************************************/
154 >  
155 >  vect r_ab; /* the vector whose origin is a and end is b */
156 >  vect r_cb; /* the vector whose origin is c and end is b */
157 >  vect r_cd; /* the vector whose origin is c and end is b */
158 >  vect r_cr1; /* the cross product of r_ab and r_cb */
159 >  vect r_cr2; /* the cross product of r_cb and r_cd */
160 >
161 >  double r_cr1_x2; /* the components of r_cr1 squared */
162 >  double r_cr1_y2;
163 >  double r_cr1_z2;
164 >  
165 >  double r_cr2_x2; /* the components of r_cr2 squared */
166 >  double r_cr2_y2;
167 >  double r_cr2_z2;
168 >
169 >  double r_cr1_sqr; /* the length of r_cr1 squared */
170 >  double r_cr2_sqr; /* the length of r_cr2 squared */
171 >  
172 >  double r_cr1_r_cr2; /* the length of r_cr1 * length of r_cr2 */
173 >  
174 >  Vector3d aR, bR, cR, dR;
175 >  Vector3d  aF, bF, cF, dF;
176 >
177 >  aR = c_p_a->getPos();
178 >  bR = c_p_b->getPos();
179 >  cR = c_p_c->getPos();
180 >  dR = c_p_d->getPos();
181 >
182 >  r_ab.x = bR[0] - aR[0];
183 >  r_ab.y = bR[1] - aR[1];
184 >  r_ab.z = bR[2] - aR[2];
185 >  r_ab.length  = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z));
186 >
187 >  r_cb.x = bR[0] - cR[0];
188 >  r_cb.y = bR[1] - cR[1];
189 >  r_cb.z = bR[2] - cR[2];
190 >  r_cb.length = sqrt((r_cb.x * r_cb.x + r_cb.y * r_cb.y + r_cb.z * r_cb.z));
191 >
192 >  r_cd.x = dR[0] - cR[0];
193 >  r_cd.y = dR[1] - cR[1];
194 >  r_cd.z = dR[2] - cR[2];
195 >  r_cd.length = sqrt((r_cd.x * r_cd.x + r_cd.y * r_cd.y + r_cd.z * r_cd.z));
196 >
197 >  r_cr1.x = r_ab.y * r_cb.z - r_cb.y * r_ab.z;
198 >  r_cr1.y = r_ab.z * r_cb.x - r_cb.z * r_ab.x;
199 >  r_cr1.z = r_ab.x * r_cb.y - r_cb.x * r_ab.y;
200 >  r_cr1_x2 = r_cr1.x * r_cr1.x;
201 >  r_cr1_y2 = r_cr1.y * r_cr1.y;
202 >  r_cr1_z2 = r_cr1.z * r_cr1.z;
203 >  r_cr1_sqr = r_cr1_x2 + r_cr1_y2 + r_cr1_z2;
204 >  r_cr1.length = sqrt(r_cr1_sqr);
205 >
206 >  r_cr2.x = r_cb.y * r_cd.z - r_cd.y * r_cb.z;
207 >  r_cr2.y = r_cb.z * r_cd.x - r_cd.z * r_cb.x;
208 >  r_cr2.z = r_cb.x * r_cd.y - r_cd.x * r_cb.y;
209 >  r_cr2_x2 = r_cr2.x * r_cr2.x;
210 >  r_cr2_y2 = r_cr2.y * r_cr2.y;
211 >  r_cr2_z2 = r_cr2.z * r_cr2.z;
212 >  r_cr2_sqr = r_cr2_x2 + r_cr2_y2 + r_cr2_z2;
213 >  r_cr2.length = sqrt(r_cr2_sqr);
214 >
215 >  r_cr1_r_cr2 = r_cr1.length * r_cr2.length;
216 >
217 >    //Vector3d pos1 = atom1_->getPos();
218 >    //Vector3d pos2 = atom2_->getPos();
219 >    //Vector3d pos3 = atom3_->getPos();
220 >    //Vector3d pos4 = atom4_->getPos();
221 >  
222 >    //Vector3d r12 = pos2 - pos1;
223 >    //Vector3d r32 = pos2 - pos3;
224 >    //Vector3d r34 = pos4 - pos3;
225 >
226 >    //A = cross(r12, r32);
227 >    //B = cross(r32, r34);
228 >    
229 >    //rA = A.length();  
230 >    //rB = B.length();
231 >  
232 >  /**********************************************************************
233 >   *
234 >   * dot product and angle calculations
235 >   *
236 >   ***********************************************************************/
237 >  
238 >  double cr1_dot_cr2; /* the dot product of the cr1 and cr2 vectors */
239 >  double cos_phi; /* the cosine of the torsion angle */
240 >
241 >  cr1_dot_cr2 = r_cr1.x * r_cr2.x + r_cr1.y * r_cr2.y + r_cr1.z * r_cr2.z;
242 >  
243 >  cos_phi = cr1_dot_cr2 / r_cr1_r_cr2;
244 >  
245 >   /* adjust for the granularity of the numbers for angles near 0 or pi */
246 >
247 >  if(cos_phi > 1.0) cos_phi = 1.0;
248 >  if(cos_phi < -1.0) cos_phi = -1.0;
249 >
250 >    //cos_phi = dot (A, B) / (rA * rB);
251 >    //if (cos_phi > 1.0) {
252 >    //  cos_phi = 1.0;
253 >    //}
254 >    //if (cos_phi < -1.0) {
255 >    //  cos_phi = -1.0;
256 >    //}
257 >
258 >
259 >
260 >  /********************************************************************
261 >   *
262 >   * This next section calculates derivatives needed for the force
263 >   * calculation
264 >   *
265 >   ********************************************************************/
266 >
267 >
268 >  /* the derivatives of cos phi with respect to the x, y,
269 >     and z components of vectors cr1 and cr2. */
270 >  double d_cos_dx_cr1;
271 >  double d_cos_dy_cr1;
272 >  double d_cos_dz_cr1;
273 >  double d_cos_dx_cr2;
274 >  double d_cos_dy_cr2;
275 >  double d_cos_dz_cr2;
276 >
277 >  d_cos_dx_cr1 = r_cr2.x / r_cr1_r_cr2 - (cos_phi * r_cr1.x) / r_cr1_sqr;
278 >  d_cos_dy_cr1 = r_cr2.y / r_cr1_r_cr2 - (cos_phi * r_cr1.y) / r_cr1_sqr;
279 >  d_cos_dz_cr1 = r_cr2.z / r_cr1_r_cr2 - (cos_phi * r_cr1.z) / r_cr1_sqr;
280 >
281 >  d_cos_dx_cr2 = r_cr1.x / r_cr1_r_cr2 - (cos_phi * r_cr2.x) / r_cr2_sqr;
282 >  d_cos_dy_cr2 = r_cr1.y / r_cr1_r_cr2 - (cos_phi * r_cr2.y) / r_cr2_sqr;
283 >  d_cos_dz_cr2 = r_cr1.z / r_cr1_r_cr2 - (cos_phi * r_cr2.z) / r_cr2_sqr;
284 >
285 >  //Vector3d dcosdA = B /(rA * rB) - cos_phi  /(rA * rA) * A;
286 >  //Vector3d dcosdA = 1.0 /rA * (B.normalize() - cos_phi * A.normalize());
287 >  //Vector3d dcosdB = 1.0 /rB * (A.normalize() - cos_phi * B.normalize());
288 >
289 >  /***********************************************************************
290 >   *
291 >   * Calculate the actual forces and place them in the atoms.
292 >   *
293 >   ***********************************************************************/
294 >
295 >  double force; /*the force scaling factor */
296 >
297 >  force = torsion_force(cos_phi);
298 >
299 >  aF[0] = force * (d_cos_dy_cr1 * r_cb.z - d_cos_dz_cr1 * r_cb.y);
300 >  aF[1] = force * (d_cos_dz_cr1 * r_cb.x - d_cos_dx_cr1 * r_cb.z);
301 >  aF[2] = force * (d_cos_dx_cr1 * r_cb.y - d_cos_dy_cr1 * r_cb.x);
302 >
303 >  bF[0] = force * (  d_cos_dy_cr1 * (r_ab.z - r_cb.z)
304 >           - d_cos_dy_cr2 *  r_cd.z  
305 >           + d_cos_dz_cr1 * (r_cb.y - r_ab.y)
306 >           + d_cos_dz_cr2 *  r_cd.y);
307 >  bF[1] = force * (  d_cos_dx_cr1 * (r_cb.z - r_ab.z)
308 >           + d_cos_dx_cr2 *  r_cd.z  
309 >           + d_cos_dz_cr1 * (r_ab.x - r_cb.x)
310 >           - d_cos_dz_cr2 *  r_cd.x);
311 >  bF[2] = force * (  d_cos_dx_cr1 * (r_ab.y - r_cb.y)
312 >           - d_cos_dx_cr2 *  r_cd.y  
313 >           + d_cos_dy_cr1 * (r_cb.x - r_ab.x)
314 >           + d_cos_dy_cr2 *  r_cd.x);
315 >
316 >  cF[0] = force * (- d_cos_dy_cr1 *  r_ab.z
317 >           - d_cos_dy_cr2 * (r_cb.z - r_cd.z)
318 >           + d_cos_dz_cr1 *  r_ab.y
319 >           - d_cos_dz_cr2 * (r_cd.y - r_cb.y));
320 >  cF[1] = force * (  d_cos_dx_cr1 *  r_ab.z
321 >           - d_cos_dx_cr2 * (r_cd.z - r_cb.z)
322 >           - d_cos_dz_cr1 *  r_ab.x
323 >           - d_cos_dz_cr2 * (r_cb.x - r_cd.x));
324 >  cF[2] = force * (- d_cos_dx_cr1 *  r_ab.y
325 >           - d_cos_dx_cr2 * (r_cb.y - r_cd.y)
326 >           + d_cos_dy_cr1 *  r_ab.x
327 >           - d_cos_dy_cr2 * (r_cd.x - r_cb.x));
328 >
329 >  dF[0] = force * (d_cos_dy_cr2 * r_cb.z - d_cos_dz_cr2 * r_cb.y);
330 >  dF[1] = force * (d_cos_dz_cr2 * r_cb.x - d_cos_dx_cr2 * r_cb.z);
331 >  dF[2] = force * (d_cos_dx_cr2 * r_cb.y - d_cos_dy_cr2 * r_cb.x);
332 >
333 >
334 >  c_p_a->addFrc(aF);
335 >  c_p_b->addFrc(bF);
336 >  c_p_c->addFrc(cF);
337 >  c_p_d->addFrc(dF);
338 >
339 >    //double firstDerivative;
340 >    //bondType_->calcForce(cos_phi, firstDerivative, potential_);
341 >    //f1 = force * cross (dcosdA, r32);
342 >    //f2 =
343 >    //f3 =
344 >    //f4 = force * cross(dcosdB, r32);
345 >    //atom1_->addFrc(f1);
346 >    //atom2_->addFrc(f2);
347 >    //atom3_->addFrc(f3);
348 >    //atom4_->addFrc(f4);
349 >
350 >  
351 > }
352 >
353 > }
354 > =======
355 > #include "primitives/Torsion.hpp"
356 >
357 > namespace oopse {
358 >
359 > Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
360 >                 TorsionType *tt) :
361 >    atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
362 >
363 > void Torsion::calcForce() {
364 >    Vector3d pos1 = atom1_->getPos();
365 >    Vector3d pos2 = atom2_->getPos();
366 >    Vector3d pos3 = atom3_->getPos();
367 >    Vector3d pos4 = atom4_->getPos();
368 >
369 >    Vector3d r12 = pos1 - pos2;
370 >    Vector3d r23 = pos2 - pos3;
371 >    Vector3d r34 = pos3 - pos4;
372 >
373 >    //  Calculate the cross products and distances
374 >    Vector3d A = cross(r12, r23);
375 >    double rA = A.length();
376 >    Vector3d B = cross(r23, r34);
377 >    double rB = B.length();
378 >    Vector3d C = cross(r23, A);
379 >    double rC = C.length();
380 >
381 >    A.normalize();
382 >    B.normalize();
383 >    C.normalize();
384 >    
385 >    //  Calculate the sin and cos
386 >    double cos_phi = dot(A, B) ;
387 >    double sin_phi = dot(C, B);
388 >
389 >    double dVdPhi;
390 >    torsionType_->calcForce(cos_phi, sin_phi, potential_, dVdPhi);
391 >
392 >    Vector3d f1;
393 >    Vector3d f2;
394 >    Vector3d f3;
395 >
396 >    //  Next, we want to calculate the forces.  In order
397 >    //  to do that, we first need to figure out whether the
398 >    //  sin or cos form will be more stable.  For this,
399 >    //  just look at the value of phi
400 >    if (fabs(sin_phi) > 0.1) {
401 >        //  use the sin version to avoid 1/cos terms
402 >
403 >        Vector3d dcosdA = (cos_phi * A - B) /rA;
404 >        Vector3d dcosdB = (cos_phi * B - A) /rB;
405 >
406 >        double dVdcosPhi = dVdPhi / sin_phi;
407 >
408 >        f1 = dVdcosPhi * cross(r23, dcosdA);
409 >        f2 = dVdcosPhi * ( cross(r34, dcosdB) - cross(r12, dcosdA));
410 >        f3 = dVdcosPhi * cross(r23, dcosdB);
411 >
412 >    } else {
413 >        //  This angle is closer to 0 or 180 than it is to
414 >        //  90, so use the cos version to avoid 1/sin terms
415 >
416 >        double dVdsinPhi = -dVdPhi /cos_phi;
417 >        Vector3d dsindB = (sin_phi * B - C) /rB;
418 >        Vector3d dsindC = (sin_phi * C - B) /rC;
419 >
420 >        f1.x() = dVdsinPhi*((r23.y()*r23.y() + r23.z()*r23.z())*dsindC.x() - r23.x()*r23.y()*dsindC.y() - r23.x()*r23.z()*dsindC.z());
421 >
422 >        f1.y() = dVdsinPhi*((r23.z()*r23.z() + r23.x()*r23.x())*dsindC.y() - r23.y()*r23.z()*dsindC.z() - r23.y()*r23.x()*dsindC.x());
423 >
424 >        f1.z() = dVdsinPhi*((r23.x()*r23.x() + r23.y()*r23.y())*dsindC.z() - r23.z()*r23.x()*dsindC.x() - r23.z()*r23.y()*dsindC.y());
425 >
426 >        f2.x() = dVdsinPhi*(-(r23.y()*r12.y() + r23.z()*r12.z())*dsindC.x() + (2.0*r23.x()*r12.y() - r12.x()*r23.y())*dsindC.y()
427 >        + (2.0*r23.x()*r12.z() - r12.x()*r23.z())*dsindC.z() + dsindB.z()*r34.y() - dsindB.y()*r34.z());
428 >
429 >        f2.y() = dVdsinPhi*(-(r23.z()*r12.z() + r23.x()*r12.x())*dsindC.y() + (2.0*r23.y()*r12.z() - r12.y()*r23.z())*dsindC.z()
430 >        + (2.0*r23.y()*r12.x() - r12.y()*r23.x())*dsindC.x() + dsindB.x()*r34.z() - dsindB.z()*r34.x());
431 >
432 >        f2.z() = dVdsinPhi*(-(r23.x()*r12.x() + r23.y()*r12.y())*dsindC.z() + (2.0*r23.z()*r12.x() - r12.z()*r23.x())*dsindC.x()
433 >        +(2.0*r23.z()*r12.y() - r12.z()*r23.y())*dsindC.y() + dsindB.y()*r34.x() - dsindB.x()*r34.y());
434 >
435 >        f3 = dVdsinPhi * cross(dsindB, r23);
436 >
437 >    }
438 >
439 >    atom1_->addFrc(f1);
440 >    atom2_->addFrc(f2 - f1);
441 >    atom3_->addFrc(f3 - f2);
442 >    atom4_->addFrc(-f3);
443 > }
444 >
445 > }
446 > >>>>>>> 1.2.2.4

Diff Legend

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