TorsionForce.h

#ifndef SIMP_TORSION_FORCE_H
#define SIMP_TORSION_FORCE_H

/*
* Simpatico - Simulation Package for Polymeric and Molecular Liquids
*
* Copyright 2010 - 2017, The Regents of the University of Minnesota
* Distributed under the terms of the GNU General Public License.
*/

#include "Torsion.h"

#include <cmath>

namespace Simp
{

   using namespace Util;

   struct TorsionForce : public Torsion
   {
   
   public:

      Vector d1;

      Vector d2;
   
      Vector d3;
   
      bool computeDerivatives(const Vector& b1, const Vector& b2, const Vector& b3);
 
   };

   // Inline method definitions

   /* 
   * Calculate cosPhi.
   */ 
   inline bool 
   TorsionForce::computeDerivatives(
                    const Vector& b1, const Vector& b2, const Vector& b3)
   {
      Vector u1, u2, t1, t2;
      double r1, r2;

      u1.cross(b1, b2);
      r1 = u1.square();
      if (r1 < 1.0E-10) {
         return 1; // Error code
      }
      r1 = sqrt(r1);
      u1 /= r1;

      u2.cross(b2, b3);
      r2 = u2.square();
      if (r2 < 1.0E-10) {
         return 1; // Error code
      }
      r2 = sqrt(r2);
      u2 /= r2;

      cosPhi = u1.dot(u2);

      t1.multiply(u1, -cosPhi);
      t1 += u2;
      t1 *= 1.0 / r1;

      t2.multiply(u2, -cosPhi);
      t2 += u1;
      t2 *= 1.0 / r2;

      d1.cross(b2, t1);
      d3.cross(t2, b2);

      d2.cross(t1, b1);
      t1.cross(b3, t2);
      d2 += t1;

      return 0; // Normal completion
   }

} 
#endif