[ViewVC] Diff of: group/branches/new_design/OOPSE-3.0/src/primitives/Torsion.cpp

Comparing branches/new_design/OOPSE-3.0/src/primitives/Torsion.cpp (file contents):
Revision 1847 by tim, Sat Dec 4 05:24:07 2004 UTC vs.
Revision 1849 by gezelter, Sat Dec 4 19:24:16 2004 UTC

-+
+<<<<<<< Torsion.cpp
+#include "primitives/Torsion.hpp"
+namespace oopse {
-+
+Torsion::Torsion(Atom* atom1, Atom* atom2, Atom* atom3, Atom* atom4, TorsionType* tt)
-+
+            : atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4) {
-+
-+
+}
-+
-+
+void Torsion::calcForce() {
-+
-+
+    Vector3d pos1 = atom1_->getPos();
-+
+    Vector3d pos2 = atom2_->getPos();
-+
+    Vector3d pos3 = atom3_->getPos();
-+
+    Vector3d pos4 = atom4_->getPos();
-+
-+
+    Vector3d r12 = pos1 - pos2;
-+
+    Vector3d r23 = pos2 - pos3;
-+
+    Vector3d r34 = pos3 - pos4;
-+
-+
+    //  Calculate the cross products and distances
-+
+    Vector3d A = cross(r12,r23);
-+
+    double rA = A.length();
-+
+    Vector3d B = cross(r23,r34);
-+
+    double rB = B.length();
-+
+    Vector3d C = cross(r23,A);
-+
+    double rC = C.length();
-+
-+
+    //  Calculate the sin and cos
-+
+    double cos_phi = (A*B)/(rA*rB);
-+
+    double sin_phi = (C*B)/(rC*rB);
-+
-+
+    double phi= -atan2(sin_phi,cos_phi);
-+
-+
+    double firstDerivative;
-+
-+
+    torsionType_->calcForce(cosPhi, sinPhi, Vtorsion, dVdCosPhi);
-+
-+
-+
+    Vector3d f1,f2,f3;
-+
-+
+    //  Normalize B
-+
+    rB = 1.0/rB;
-+
+    B *= rB;
-+
-+
+    //  Next, we want to calculate the forces.  In order
-+
+    //  to do that, we first need to figure out whether the
-+
+    //  sin or cos form will be more stable.  For this,
-+
+    //  just look at the value of phi
-+
+    if (fabs(sin_phi) > 0.1) {
-+
+        //  use the sin version to avoid 1/cos terms
-+
-+
+        rA = 1.0/rA;
-+
+        A *= rA;
-+
+        Vector3d dcosdA = rA*(cos_phi*A-B);
-+
+        Vector3d dcosdB = rB*(cos_phi*B-A);
-+
-+
+        K1 = K1/sin_phi;
-+
-+
+        //simple form
-+
+        //f1 = K1 * cross(r23, dcosdA);
-+
+        //f3 = K1 * cross(r23, dcosdB);
-+
+        //f2 = K1 * ( cross(r34, dcosdB) - cross(r12, dcosdA));
-+
-+
+        f1.x = K1*(r23.y*dcosdA.z - r23.z*dcosdA.y);
-+
+        f1.y = K1*(r23.z*dcosdA.x - r23.x*dcosdA.z);
-+
+        f1.z = K1*(r23.x*dcosdA.y - r23.y*dcosdA.x);
-+
-+
+        f3.x = K1*(r23.z*dcosdB.y - r23.y*dcosdB.z);
-+
+        f3.y = K1*(r23.x*dcosdB.z - r23.z*dcosdB.x);
-+
+        f3.z = K1*(r23.y*dcosdB.x - r23.x*dcosdB.y);
-+
-+
+        f2.x = K1*(r12.z*dcosdA.y - r12.y*dcosdA.z + r34.y*dcosdB.z - r34.z*dcosdB.y);
-+
+        f2.y = K1*(r12.x*dcosdA.z - r12.z*dcosdA.x + r34.z*dcosdB.x - r34.x*dcosdB.z);
-+
+        f2.z = K1*(r12.y*dcosdA.x - r12.x*dcosdA.y + r34.x*dcosdB.y - r34.y*dcosdB.x);
-+
+    } else {
-+
+        //  This angle is closer to 0 or 180 than it is to
-+
+        //  90, so use the cos version to avoid 1/sin terms
-+
-+
+        //  Normalize C
-+
+        rC = 1.0/rC;
-+
+        C *= rC;
-+
+        Vector3d dsindC = rC*(sin_phi*C-B);
-+
+        Vector3d dsindB = rB*(sin_phi*B-C);
-+
-+
+        K1 = -K1/cos_phi;
-+
-+
+        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);
-+
+        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);
-+
+        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);
-+
-+
+        f3 = K1 *cross(dsindB,r23);
-+
-+
+        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
-+
+                 + (2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z + dsindB.z*r34.y - dsindB.y*r34.z);
-+
+        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
-+
+                 + (2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x + dsindB.x*r34.z - dsindB.z*r34.x);
-+
+        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
-+
+                  +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y + dsindB.y*r34.x - dsindB.x*r34.y);
-+
+    }
-+
-+
+    atom1_->addFrc(f1);
-+
+    atom2_->addFrc(f2 - f1);
-+
+    atom3_->addFrc(f3 - f2);
-+
+    atom4_->addFrc(-f3);
-+
-+
+}
-+
-+
-+
+      double K=0;       // energy
-+
+      double K1=0;      // force
-+
-+
+  // get the dihedral information
-+
+  int multiplicity = value->multiplicity;
-+
-+
+  //  Loop through the multiple parameter sets for this
-+
+  //  bond.  We will only loop more than once if this
-+
+  //  has multiple parameter sets from Charmm22
-+
+  for (int mult_num=0; mult_num<multiplicity; mult_num++)
-+
+  {
-+
+    /* get angle information */
-+
+    double k = value->values[mult_num].k * scale;
-+
+    double delta = value->values[mult_num].delta;
-+
+    int n = value->values[mult_num].n;
-+
-+
+    //  Calculate the energy
-+
+    if (n)
-+
+    {
-+
+      //  Periodicity is greater than 0, so use cos form
-+
+      K += k*(1+cos(n*phi + delta));
-+
+      K1 += -n*k*sin(n*phi + delta);
-+
+    }
-+
+    else
-+
+    {
-+
+      //  Periodicity is 0, so just use the harmonic form
-+
+      double diff = phi-delta;
-+
+      if (diff < -PI)           diff += TWOPI;
-+
+      else if (diff > PI)       diff -= TWOPI;
-+
-+
+      K += k*diff*diff;
-+
+      K1 += 2.0*k*diff;
-+
+    }
-+
+  } /* for multiplicity */
-+
-+
-+
+void Torsion::calc_forces(){
-+
-+
+  /**********************************************************************
-+
+   *
-+
+   * initialize vectors
-+
+   *
-+
+   ***********************************************************************/
-+
-+
+  vect r_ab; /* the vector whose origin is a and end is b */
-+
+  vect r_cb; /* the vector whose origin is c and end is b */
-+
+  vect r_cd; /* the vector whose origin is c and end is b */
-+
+  vect r_cr1; /* the cross product of r_ab and r_cb */
-+
+  vect r_cr2; /* the cross product of r_cb and r_cd */
-+
-+
+  double r_cr1_x2; /* the components of r_cr1 squared */
-+
+  double r_cr1_y2;
-+
+  double r_cr1_z2;
-+
-+
+  double r_cr2_x2; /* the components of r_cr2 squared */
-+
+  double r_cr2_y2;
-+
+  double r_cr2_z2;
-+
-+
+  double r_cr1_sqr; /* the length of r_cr1 squared */
-+
+  double r_cr2_sqr; /* the length of r_cr2 squared */
-+
-+
+  double r_cr1_r_cr2; /* the length of r_cr1 * length of r_cr2 */
-+
-+
+  Vector3d aR, bR, cR, dR;
-+
+  Vector3d  aF, bF, cF, dF;
-+
-+
+  aR = c_p_a->getPos();
-+
+  bR = c_p_b->getPos();
-+
+  cR = c_p_c->getPos();
-+
+  dR = c_p_d->getPos();
-+
-+
+  r_ab.x = bR[0] - aR[0];
-+
+  r_ab.y = bR[1] - aR[1];
-+
+  r_ab.z = bR[2] - aR[2];
-+
+  r_ab.length  = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z));
-+
-+
+  r_cb.x = bR[0] - cR[0];
-+
+  r_cb.y = bR[1] - cR[1];
-+
+  r_cb.z = bR[2] - cR[2];
-+
+  r_cb.length = sqrt((r_cb.x * r_cb.x + r_cb.y * r_cb.y + r_cb.z * r_cb.z));
-+
-+
+  r_cd.x = dR[0] - cR[0];
-+
+  r_cd.y = dR[1] - cR[1];
-+
+  r_cd.z = dR[2] - cR[2];
-+
+  r_cd.length = sqrt((r_cd.x * r_cd.x + r_cd.y * r_cd.y + r_cd.z * r_cd.z));
-+
-+
+  r_cr1.x = r_ab.y * r_cb.z - r_cb.y * r_ab.z;
-+
+  r_cr1.y = r_ab.z * r_cb.x - r_cb.z * r_ab.x;
-+
+  r_cr1.z = r_ab.x * r_cb.y - r_cb.x * r_ab.y;
-+
+  r_cr1_x2 = r_cr1.x * r_cr1.x;
-+
+  r_cr1_y2 = r_cr1.y * r_cr1.y;
-+
+  r_cr1_z2 = r_cr1.z * r_cr1.z;
-+
+  r_cr1_sqr = r_cr1_x2 + r_cr1_y2 + r_cr1_z2;
-+
+  r_cr1.length = sqrt(r_cr1_sqr);
-+
-+
+  r_cr2.x = r_cb.y * r_cd.z - r_cd.y * r_cb.z;
-+
+  r_cr2.y = r_cb.z * r_cd.x - r_cd.z * r_cb.x;
-+
+  r_cr2.z = r_cb.x * r_cd.y - r_cd.x * r_cb.y;
-+
+  r_cr2_x2 = r_cr2.x * r_cr2.x;
-+
+  r_cr2_y2 = r_cr2.y * r_cr2.y;
-+
+  r_cr2_z2 = r_cr2.z * r_cr2.z;
-+
+  r_cr2_sqr = r_cr2_x2 + r_cr2_y2 + r_cr2_z2;
-+
+  r_cr2.length = sqrt(r_cr2_sqr);
-+
-+
+  r_cr1_r_cr2 = r_cr1.length * r_cr2.length;
-+
-+
+    //Vector3d pos1 = atom1_->getPos();
-+
+    //Vector3d pos2 = atom2_->getPos();
-+
+    //Vector3d pos3 = atom3_->getPos();
-+
+    //Vector3d pos4 = atom4_->getPos();
-+
-+
+    //Vector3d r12 = pos2 - pos1;
-+
+    //Vector3d r32 = pos2 - pos3;
-+
+    //Vector3d r34 = pos4 - pos3;
-+
-+
+    //A = cross(r12, r32);
-+
+    //B = cross(r32, r34);
-+
-+
+    //rA = A.length();
-+
+    //rB = B.length();
-+
-+
+  /**********************************************************************
-+
+   *
-+
+   * dot product and angle calculations
-+
+   *
-+
+   ***********************************************************************/
-+
-+
+  double cr1_dot_cr2; /* the dot product of the cr1 and cr2 vectors */
-+
+  double cos_phi; /* the cosine of the torsion angle */
-+
-+
+  cr1_dot_cr2 = r_cr1.x * r_cr2.x + r_cr1.y * r_cr2.y + r_cr1.z * r_cr2.z;
-+
-+
+  cos_phi = cr1_dot_cr2 / r_cr1_r_cr2;
-+
-+
+   /* adjust for the granularity of the numbers for angles near 0 or pi */
-+
-+
+  if(cos_phi > 1.0) cos_phi = 1.0;
-+
+  if(cos_phi < -1.0) cos_phi = -1.0;
-+
-+
+    //cos_phi = dot (A, B) / (rA * rB);
-+
+    //if (cos_phi > 1.0) {
-+
+    //  cos_phi = 1.0;
-+
+    //}
-+
+    //if (cos_phi < -1.0) {
-+
+    //  cos_phi = -1.0;
-+
+    //}
-+
-+
-+
-+
+  /********************************************************************
-+
+   *
-+
+   * This next section calculates derivatives needed for the force
-+
+   * calculation
-+
+   *
-+
+   ********************************************************************/
-+
-+
-+
+  /* the derivatives of cos phi with respect to the x, y,
-+
+     and z components of vectors cr1 and cr2. */
-+
+  double d_cos_dx_cr1;
-+
+  double d_cos_dy_cr1;
-+
+  double d_cos_dz_cr1;
-+
+  double d_cos_dx_cr2;
-+
+  double d_cos_dy_cr2;
-+
+  double d_cos_dz_cr2;
-+
-+
+  d_cos_dx_cr1 = r_cr2.x / r_cr1_r_cr2 - (cos_phi * r_cr1.x) / r_cr1_sqr;
-+
+  d_cos_dy_cr1 = r_cr2.y / r_cr1_r_cr2 - (cos_phi * r_cr1.y) / r_cr1_sqr;
-+
+  d_cos_dz_cr1 = r_cr2.z / r_cr1_r_cr2 - (cos_phi * r_cr1.z) / r_cr1_sqr;
-+
-+
+  d_cos_dx_cr2 = r_cr1.x / r_cr1_r_cr2 - (cos_phi * r_cr2.x) / r_cr2_sqr;
-+
+  d_cos_dy_cr2 = r_cr1.y / r_cr1_r_cr2 - (cos_phi * r_cr2.y) / r_cr2_sqr;
-+
+  d_cos_dz_cr2 = r_cr1.z / r_cr1_r_cr2 - (cos_phi * r_cr2.z) / r_cr2_sqr;
-+
-+
+  //Vector3d dcosdA = B /(rA * rB) - cos_phi  /(rA * rA) * A;
-+
+  //Vector3d dcosdA = 1.0 /rA * (B.normalize() - cos_phi * A.normalize());
-+
+  //Vector3d dcosdB = 1.0 /rB * (A.normalize() - cos_phi * B.normalize());
-+
-+
+  /***********************************************************************
-+
+   *
-+
+   * Calculate the actual forces and place them in the atoms.
-+
+   *
-+
+   ***********************************************************************/
-+
-+
+  double force; /*the force scaling factor */
-+
-+
+  force = torsion_force(cos_phi);
-+
-+
+  aF[0] = force * (d_cos_dy_cr1 * r_cb.z - d_cos_dz_cr1 * r_cb.y);
-+
+  aF[1] = force * (d_cos_dz_cr1 * r_cb.x - d_cos_dx_cr1 * r_cb.z);
-+
+  aF[2] = force * (d_cos_dx_cr1 * r_cb.y - d_cos_dy_cr1 * r_cb.x);
-+
-+
+  bF[0] = force * (  d_cos_dy_cr1 * (r_ab.z - r_cb.z)
-+
+           - d_cos_dy_cr2 *  r_cd.z
-+
+           + d_cos_dz_cr1 * (r_cb.y - r_ab.y)
-+
+           + d_cos_dz_cr2 *  r_cd.y);
-+
+  bF[1] = force * (  d_cos_dx_cr1 * (r_cb.z - r_ab.z)
-+
+           + d_cos_dx_cr2 *  r_cd.z
-+
+           + d_cos_dz_cr1 * (r_ab.x - r_cb.x)
-+
+           - d_cos_dz_cr2 *  r_cd.x);
-+
+  bF[2] = force * (  d_cos_dx_cr1 * (r_ab.y - r_cb.y)
-+
+           - d_cos_dx_cr2 *  r_cd.y
-+
+           + d_cos_dy_cr1 * (r_cb.x - r_ab.x)
-+
+           + d_cos_dy_cr2 *  r_cd.x);
-+
-+
+  cF[0] = force * (- d_cos_dy_cr1 *  r_ab.z
-+
+           - d_cos_dy_cr2 * (r_cb.z - r_cd.z)
-+
+           + d_cos_dz_cr1 *  r_ab.y
-+
+           - d_cos_dz_cr2 * (r_cd.y - r_cb.y));
-+
+  cF[1] = force * (  d_cos_dx_cr1 *  r_ab.z
-+
+           - d_cos_dx_cr2 * (r_cd.z - r_cb.z)
-+
+           - d_cos_dz_cr1 *  r_ab.x
-+
+           - d_cos_dz_cr2 * (r_cb.x - r_cd.x));
-+
+  cF[2] = force * (- d_cos_dx_cr1 *  r_ab.y
-+
+           - d_cos_dx_cr2 * (r_cb.y - r_cd.y)
-+
+           + d_cos_dy_cr1 *  r_ab.x
-+
+           - d_cos_dy_cr2 * (r_cd.x - r_cb.x));
-+
-+
+  dF[0] = force * (d_cos_dy_cr2 * r_cb.z - d_cos_dz_cr2 * r_cb.y);
-+
+  dF[1] = force * (d_cos_dz_cr2 * r_cb.x - d_cos_dx_cr2 * r_cb.z);
-+
+  dF[2] = force * (d_cos_dx_cr2 * r_cb.y - d_cos_dy_cr2 * r_cb.x);
-+
-+
-+
+  c_p_a->addFrc(aF);
-+
+  c_p_b->addFrc(bF);
-+
+  c_p_c->addFrc(cF);
-+
+  c_p_d->addFrc(dF);
-+
-+
+    //double firstDerivative;
-+
+    //bondType_->calcForce(cos_phi, firstDerivative, potential_);
-+
+    //f1 = force * cross (dcosdA, r32);
-+
+    //f2 =
-+
+    //f3 =
-+
+    //f4 = force * cross(dcosdB, r32);
-+
+    //atom1_->addFrc(f1);
-+
+    //atom2_->addFrc(f2);
-+
+    //atom3_->addFrc(f3);
-+
+    //atom4_->addFrc(f4);
-+
-+
-+
+}
-+
-+
+}
-+
+=======
-+
+#include "primitives/Torsion.hpp"
-+
-+
+namespace oopse {
-+
+Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
+                 TorsionType *tt) :
+    atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
-+
+>>>>>>> 1.2.2.4

Diff Legend

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

Comparing branches/new_design/OOPSE-3.0/src/primitives/Torsion.cpp (file contents): Revision 1847 by tim, Sat Dec 4 05:24:07 2004 UTC vs. Revision 1849 by gezelter, Sat Dec 4 19:24:16 2004 UTC

Diff Legend

Comparing branches/new_design/OOPSE-3.0/src/primitives/Torsion.cpp (file contents):
Revision 1847 by tim, Sat Dec 4 05:24:07 2004 UTC vs.
Revision 1849 by gezelter, Sat Dec 4 19:24:16 2004 UTC