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 1849 by gezelter, Sat Dec 4 19:24:16 2004 UTC vs.
Revision 1850 by tim, Sat Dec 4 19:29:54 2004 UTC

# Line 1 | Line 1
1 – <<<<<<< Torsion.cpp
1   #include "primitives/Torsion.hpp"
2  
3   namespace oopse {
4  
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 –
5   Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
6                   TorsionType *tt) :
7      atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
# Line 443 | Line 89 | void Torsion::calcForce() {
89   }
90  
91   }
446 – >>>>>>> 1.2.2.4

Diff Legend

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