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 1683, Thu Oct 28 22:34:02 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>
1 > <<<<<<< Torsion.cpp
2 > #include "primitives/Torsion.hpp"
3  
4 < void Torsion::set_atoms( Atom &a, Atom &b, Atom &c, Atom &d){
5 <  c_p_a = &a;
6 <  c_p_b = &b;
7 <  c_p_c = &c;
8 <  c_p_d = &d;
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    /**********************************************************************
# Line 39 | Line 171 | void Torsion::calc_forces(){
171    
172    double r_cr1_r_cr2; /* the length of r_cr1 * length of r_cr2 */
173    
174 <  double aR[3], bR[3], cR[3], dR[3];
175 <  double aF[3], bF[3], cF[3], dF[3];
174 >  Vector3d aR, bR, cR, dR;
175 >  Vector3d  aF, bF, cF, dF;
176  
177 <  c_p_a->getPos( aR );
178 <  c_p_b->getPos( bR );
179 <  c_p_c->getPos( cR );
180 <  c_p_d->getPos( dR );
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];
# Line 81 | Line 213 | void Torsion::calc_forces(){
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
# Line 100 | Line 247 | void Torsion::calc_forces(){
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
# Line 126 | Line 282 | void Torsion::calc_forces(){
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.
# Line 141 | Line 301 | void Torsion::calc_forces(){
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);
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);
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);
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));
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));
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));
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);
# Line 175 | Line 335 | void Torsion::calc_forces(){
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