src/primitives/Torsion.cpp

<<<<<<< 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) { }

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();

    A.normalize();
    B.normalize();
    C.normalize();
    
    //  Calculate the sin and cos
    double cos_phi = dot(A, B) ;
    double sin_phi = dot(C, B);

    double dVdPhi;
    torsionType_->calcForce(cos_phi, sin_phi, potential_, dVdPhi);

    Vector3d f1;
    Vector3d f2;
    Vector3d f3;

    //  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

        Vector3d dcosdA = (cos_phi * A - B) /rA;
        Vector3d dcosdB = (cos_phi * B - A) /rB;

        double dVdcosPhi = dVdPhi / sin_phi;

        f1 = dVdcosPhi * cross(r23, dcosdA);
        f2 = dVdcosPhi * ( cross(r34, dcosdB) - cross(r12, dcosdA));
        f3 = dVdcosPhi * cross(r23, dcosdB);

    } else {
        //  This angle is closer to 0 or 180 than it is to
        //  90, so use the cos version to avoid 1/sin terms

        double dVdsinPhi = -dVdPhi /cos_phi;
        Vector3d dsindB = (sin_phi * B - C) /rB;
        Vector3d dsindC = (sin_phi * C - B) /rC;

        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());

        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());

        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());

        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()
        + (2.0*r23.x()*r12.z() - r12.x()*r23.z())*dsindC.z() + dsindB.z()*r34.y() - dsindB.y()*r34.z());

        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()
        + (2.0*r23.y()*r12.x() - r12.y()*r23.x())*dsindC.x() + dsindB.x()*r34.z() - dsindB.z()*r34.x());

        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()
        +(2.0*r23.z()*r12.y() - r12.z()*r23.y())*dsindC.y() + dsindB.y()*r34.x() - dsindB.x()*r34.y());

        f3 = dVdsinPhi * cross(dsindB, r23);

    }

    atom1_->addFrc(f1);
    atom2_->addFrc(f2 - f1);
    atom3_->addFrc(f3 - f2);
    atom4_->addFrc(-f3);
}

}
>>>>>>> 1.2.2.4
Revision:	1849
Committed:	Sat Dec 4 19:24:16 2004 UTC (19 years, 9 months ago) by gezelter
File size:	14530 byte(s)
Log Message:	What?
#	User	Rev	Content
1	gezelter	1849	<<<<<<< Torsion.cpp
2	tim	1746	#include "primitives/Torsion.hpp"
3
4			namespace oopse {
5
6	gezelter	1849	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 = (AB)/(rArB);
32			double sin_phi = (CB)/(rCrB);
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_phiA-B);
57			Vector3d dcosdB = rB(cos_phiB-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.ydcosdA.z - r23.z*dcosdA.y);
67			f1.y = K1(r23.zdcosdA.x - r23.x*dcosdA.z);
68			f1.z = K1(r23.xdcosdA.y - r23.y*dcosdA.x);
69
70			f3.x = K1(r23.zdcosdB.y - r23.y*dcosdB.z);
71			f3.y = K1(r23.xdcosdB.z - r23.z*dcosdB.x);
72			f3.z = K1(r23.ydcosdB.x - r23.x*dcosdB.y);
73
74			f2.x = K1(r12.zdcosdA.y - r12.ydcosdA.z + r34.ydcosdB.z - r34.z*dcosdB.y);
75			f2.y = K1(r12.xdcosdA.z - r12.zdcosdA.x + r34.zdcosdB.x - r34.x*dcosdB.z);
76			f2.z = K1(r12.ydcosdA.x - r12.xdcosdA.y + r34.xdcosdB.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_phiC-B);
85			Vector3d dsindB = rB(sin_phiB-C);
86
87			K1 = -K1/cos_phi;
88
89			f1.x = K1((r23.yr23.y + r23.zr23.z)dsindC.x - r23.xr23.ydsindC.y - r23.xr23.zdsindC.z);
90			f1.y = K1((r23.zr23.z + r23.xr23.x)dsindC.y - r23.yr23.zdsindC.z - r23.yr23.xdsindC.x);
91			f1.z = K1((r23.xr23.x + r23.yr23.y)dsindC.z - r23.zr23.xdsindC.x - r23.zr23.ydsindC.y);
92
93			f3 = K1 *cross(dsindB,r23);
94
95			f2.x = K1(-(r23.yr12.y + r23.zr12.z)dsindC.x + (2.0r23.xr12.y - r12.xr23.y)dsindC.y
96			+ (2.0r23.xr12.z - r12.xr23.z)dsindC.z + dsindB.zr34.y - dsindB.yr34.z);
97			f2.y = K1(-(r23.zr12.z + r23.xr12.x)dsindC.y + (2.0r23.yr12.z - r12.yr23.z)dsindC.z
98			+ (2.0r23.yr12.x - r12.yr23.x)dsindC.x + dsindB.xr34.z - dsindB.zr34.x);
99			f2.z = K1(-(r23.xr12.x + r23.yr12.y)dsindC.z + (2.0r23.zr12.x - r12.zr23.x)dsindC.x
100			+(2.0r23.zr12.y - r12.zr23.y)dsindC.y + dsindB.yr34.x - dsindB.xr34.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(nphi + delta));
132			K1 += -nksin(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 += kdiffdiff;
142			K1 += 2.0kdiff;
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	tim	1746	Torsion::Torsion(Atom atom1, Atom atom2, Atom atom3, Atom atom4,
360			TorsionType *tt) :
361	tim	1847	atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
362	tim	1746
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	tim	1782	torsionType_->calcForce(cos_phi, sin_phi, potential_, dVdPhi);
391	tim	1746
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	tim	1782	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	tim	1746
422	tim	1782	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	tim	1746
424	tim	1782	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	tim	1746
426	tim	1782	f2.x() = dVdsinPhi(-(r23.y()r12.y() + r23.z()r12.z())dsindC.x() + (2.0r23.x()r12.y() - r12.x()r23.y())dsindC.y()
427			+ (2.0r23.x()r12.z() - r12.x()r23.z())dsindC.z() + dsindB.z()r34.y() - dsindB.y()r34.z());
428	tim	1746
429	tim	1782	f2.y() = dVdsinPhi(-(r23.z()r12.z() + r23.x()r12.x())dsindC.y() + (2.0r23.y()r12.z() - r12.y()r23.z())dsindC.z()
430			+ (2.0r23.y()r12.x() - r12.y()r23.x())dsindC.x() + dsindB.x()r34.z() - dsindB.z()r34.x());
431	tim	1746
432	tim	1782	f2.z() = dVdsinPhi(-(r23.x()r12.x() + r23.y()r12.y())dsindC.z() + (2.0r23.z()r12.x() - r12.z()r23.x())dsindC.x()
433			+(2.0r23.z()r12.y() - r12.z()r23.y())dsindC.y() + dsindB.y()r34.x() - dsindB.x()r34.y());
434	tim	1746
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	gezelter	1849	>>>>>>> 1.2.2.4