DihedralElem Class Reference

#include <ComputeDihedrals.h>

Public Types
enum	{ size = 4 }
enum	{   dihedralEnergyIndex, dihedralEnergyIndex_f, dihedralEnergyIndex_ti_1, dihedralEnergyIndex_ti_2,   TENSOR =(virialIndex), reductionDataSize }
enum	{ reductionChecksumLabel = REDUCTION_DIHEDRAL_CHECKSUM }

Public Member Functions
int	hash () const
	DihedralElem ()
	DihedralElem (AtomID atom0, const TupleSignature *sig, const DihedralValue *v)
	DihedralElem (const Dihedral *a, const DihedralValue *v)
	DihedralElem (AtomID atom0, AtomID atom1, AtomID atom2, AtomID atom3)
	~DihedralElem ()
int	operator== (const DihedralElem &a) const
int	operator< (const DihedralElem &a) const

Static Public Member Functions
static void	computeForce (DihedralElem *, int, BigReal *, BigReal *)
static void	getMoleculePointers (Molecule *, int *, int32 ***, Dihedral **)
static void	getParameterPointers (Parameters *, const DihedralValue **)
static void	getTupleInfo (AtomSignature *sig, int *count, TupleSignature **t)
static void	submitReductionData (BigReal *, SubmitReduction *)

Public Attributes
AtomID	atomID [size]
int	localIndex [size]
TuplePatchElem *	p [size]
Real	scale
const DihedralValue *	value

Static Public Attributes
static int	pressureProfileSlabs = 0
static int	pressureProfileAtomTypes = 1
static BigReal	pressureProfileThickness = 0
static BigReal	pressureProfileMin = 0

Detailed Description

Definition at line 17 of file ComputeDihedrals.h.

Member Enumeration Documentation

anonymous enum

Enumerator
size

Definition at line 20 of file ComputeDihedrals.h.

{ size = 4 };

anonymous enum

Enumerator
dihedralEnergyIndex
dihedralEnergyIndex_f
dihedralEnergyIndex_ti_1
dihedralEnergyIndex_ti_2
TENSOR
reductionDataSize

Definition at line 47 of file ComputeDihedrals.h.

       { dihedralEnergyIndex, dihedralEnergyIndex_f, dihedralEnergyIndex_ti_1, 
         dihedralEnergyIndex_ti_2,TENSOR(virialIndex), reductionDataSize };

anonymous enum

Enumerator
reductionChecksumLabel

Definition at line 49 of file ComputeDihedrals.h.

{ reductionChecksumLabel = REDUCTION_DIHEDRAL_CHECKSUM };

Constructor & Destructor Documentation

DihedralElem::DihedralElem ( )

inline

Definition at line 12 of file ComputeDihedrals.inl.

{ ; }

DihedralElem::DihedralElem ( AtomID  atom0, const TupleSignature *  sig, const DihedralValue *  v  )

inline

Definition at line 14 of file ComputeDihedrals.inl.

References atomID, TupleSignature::offset, TupleSignature::tupleParamType, and value.

                                                                                                {
    atomID[0] = atom0;
    atomID[1] = atom0 + sig->offset[0];
    atomID[2] = atom0 + sig->offset[1];
    atomID[3] = atom0 + sig->offset[2];
    value = &v[sig->tupleParamType];
}

DihedralElem::DihedralElem ( const Dihedral *  a, const DihedralValue *  v  )

inline

Definition at line 22 of file ComputeDihedrals.inl.

References dihedral::atom1, dihedral::atom2, dihedral::atom3, dihedral::atom4, atomID, dihedral::dihedral_type, and value.

  {
    atomID[0] = a->atom1;
    atomID[1] = a->atom2;
    atomID[2] = a->atom3;
    atomID[3] = a->atom4;
    value = &v[a->dihedral_type];
  }

DihedralElem::DihedralElem ( AtomID  atom0, AtomID  atom1, AtomID  atom2, AtomID  atom3  )

inline

Definition at line 31 of file ComputeDihedrals.inl.

References atomID.

  {
    if (atom0 > atom3) {  // Swap end atoms so lowest is first!
      AtomID tmp = atom3; atom3 = atom0; atom0 = tmp; 
      tmp = atom1; atom1 = atom2; atom2 = tmp;
    }
    atomID[0] = atom0;
    atomID[1] = atom1;
    atomID[2] = atom2;
    atomID[3] = atom3;
  }

DihedralElem::~DihedralElem ( )

inline

Definition at line 56 of file ComputeDihedrals.h.

{};

Member Function Documentation

void DihedralElem::computeForce ( DihedralElem *  tuples, int  ntuple, BigReal *  reduction, BigReal *  pressureProfileData  )

static

Definition at line 40 of file ComputeDihedrals.C.

References A, TuplePatchElem::af, SimParameters::alchOn, atomID, B, cross(), DebugM, four_body_consts::delta, Lattice::delta(), dihedralEnergyIndex, dihedralEnergyIndex_f, dihedralEnergyIndex_ti_1, dihedralEnergyIndex_ti_2, TuplePatchElem::f, Patch::flags, Molecule::get_fep_bonded_type(), SimParameters::getBondLambda(), SimParameters::getCurrentLambda(), SimParameters::getCurrentLambda2(), four_body_consts::k, Patch::lattice, localIndex, Node::molecule, DihedralValue::multiplicity, four_body_consts::n, Node::Object(), p, TuplePatchElem::p, CompAtom::partition, PI, CompAtom::position, pp_clamp(), pp_reduction(), pressureProfileAtomTypes, pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, Vector::rlength(), scale, Node::simParameters, simParams, SimParameters::singleTopology, size, Flags::step, TWOPI, value, DihedralValue::values, TuplePatchElem::x, Vector::x, Vector::y, and Vector::z.

{
 const Lattice & lattice = tuples[0].p[0]->p->lattice;

 //fepb BKR
 SimParameters *const simParams = Node::Object()->simParameters;
 const int step = tuples[0].p[0]->p->flags.step;
 const BigReal alchLambda = simParams->getCurrentLambda(step);
 const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
 const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
 const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
 const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
 const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
 Molecule *const mol = Node::Object()->molecule;
 //fepe

 for ( int ituple=0; ituple<ntuple; ++ituple ) {
  const DihedralElem &tup = tuples[ituple];
  enum { size = 4 };
  const AtomID (&atomID)[size](tup.atomID);
  const int    (&localIndex)[size](tup.localIndex);
  TuplePatchElem * const(&p)[size](tup.p);
  const Real (&scale)(tup.scale);
  const DihedralValue * const(&value)(tup.value);

  DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
               << localIndex[1] << " " << localIndex[2] << std::endl);

  //  Calculate the vectors between atoms
  const Position & pos0 = p[0]->x[localIndex[0]].position;
  const Position & pos1 = p[1]->x[localIndex[1]].position;
  const Vector r12 = lattice.delta(pos0,pos1);
  const Position & pos2 = p[2]->x[localIndex[2]].position;
  const Vector r23 = lattice.delta(pos1,pos2);
  const Position & pos3 = p[3]->x[localIndex[3]].position;
  const Vector r34 = lattice.delta(pos2,pos3);

  //  Calculate the cross products and distances
  Vector A = cross(r12,r23);
  register  BigReal rAinv = A.rlength();
  Vector B = cross(r23,r34);
  register  BigReal rBinv = B.rlength();
  Vector C = cross(r23,A);
  register  BigReal rCinv = C.rlength();

  //  Calculate the sin and cos
  BigReal cos_phi = (A*B)*(rAinv*rBinv);
  BigReal sin_phi = (C*B)*(rCinv*rBinv);

  BigReal phi= -atan2(sin_phi,cos_phi);

  BigReal K=0;          // energy
  BigReal 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 */
    Real k = value->values[mult_num].k * scale;
    Real 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
      BigReal 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 */

  //fepb - BKR scaling of alchemical bonded terms
  //       NB: TI derivative is the _unscaled_ energy.
  if ( simParams->alchOn && !simParams->singleTopology) {
    switch ( mol->get_fep_bonded_type(atomID, 4) ) {
    case 1:
      reduction[dihedralEnergyIndex_ti_1] += K;
      reduction[dihedralEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*K;
      K *= bond_lambda_1;
      K1 *= bond_lambda_1;
      break;
    case 2:
      reduction[dihedralEnergyIndex_ti_2] += K;
      reduction[dihedralEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*K;
      K *= bond_lambda_2;
      K1 *= bond_lambda_2;
      break;
    }
  }
  //fepe

  Force f1,f2,f3;

  //  Normalize B
  //rB = 1.0/rB;
  B *= rBinv;

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

      //  Normalize A
      A *= rAinv;
      Vector dcosdA;
      Vector dcosdB;

      dcosdA.x = rAinv*(cos_phi*A.x-B.x);
      dcosdA.y = rAinv*(cos_phi*A.y-B.y);
      dcosdA.z = rAinv*(cos_phi*A.z-B.z);
            
      dcosdB.x = rBinv*(cos_phi*B.x-A.x);
      dcosdB.y = rBinv*(cos_phi*B.y-A.y);
      dcosdB.z = rBinv*(cos_phi*B.z-A.z);

      K1 = K1/sin_phi;

      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 *= rCinv;
      
      Vector dsindC;
      Vector dsindB;

      dsindC.x = rCinv*(sin_phi*C.x-B.x);
      dsindC.y = rCinv*(sin_phi*C.y-B.y);
      dsindC.z = rCinv*(sin_phi*C.z-B.z);

      dsindB.x = rBinv*(sin_phi*B.x-C.x);
      dsindB.y = rBinv*(sin_phi*B.y-C.y);
      dsindB.z = rBinv*(sin_phi*B.z-C.z);

      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 = cross(K1,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);
    }

  /* store the forces */
  //  p[0]->f[localIndex[0]] += f1;
  //  p[1]->f[localIndex[1]] += f2 - f1;
  //  p[2]->f[localIndex[2]] += f3 - f2;
  //  p[3]->f[localIndex[3]] += -f3;

  p[0]->f[localIndex[0]].x += f1.x;
  p[0]->f[localIndex[0]].y += f1.y;
  p[0]->f[localIndex[0]].z += f1.z;

  p[1]->f[localIndex[1]].x += f2.x - f1.x;
  p[1]->f[localIndex[1]].y += f2.y - f1.y;
  p[1]->f[localIndex[1]].z += f2.z - f1.z;

  p[2]->f[localIndex[2]].x += f3.x - f2.x;
  p[2]->f[localIndex[2]].y += f3.y - f2.y;
  p[2]->f[localIndex[2]].z += f3.z - f2.z;

  p[3]->f[localIndex[3]].x += -f3.x;
  p[3]->f[localIndex[3]].y += -f3.y;
  p[3]->f[localIndex[3]].z += -f3.z;  

    /* store the force for dihedral-only accelMD */
  if ( p[0]->af ) {
    p[0]->af[localIndex[0]].x += f1.x;
    p[0]->af[localIndex[0]].y += f1.y;
    p[0]->af[localIndex[0]].z += f1.z;

    p[1]->af[localIndex[1]].x += f2.x - f1.x;
    p[1]->af[localIndex[1]].y += f2.y - f1.y;
    p[1]->af[localIndex[1]].z += f2.z - f1.z;

    p[2]->af[localIndex[2]].x += f3.x - f2.x;
    p[2]->af[localIndex[2]].y += f3.y - f2.y;
    p[2]->af[localIndex[2]].z += f3.z - f2.z;

    p[3]->af[localIndex[3]].x += -f3.x;
    p[3]->af[localIndex[3]].y += -f3.y;
    p[3]->af[localIndex[3]].z += -f3.z;
  }

  DebugM(3, "::computeForce() -- ending with delta energy " << K << std::endl);
  reduction[dihedralEnergyIndex] += K;
  reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
  reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
  reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
  reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
  reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
  reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
  reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
  reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
  reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );

  if (pressureProfileData) {
    BigReal z1 = p[0]->x[localIndex[0]].position.z;
    BigReal z2 = p[1]->x[localIndex[1]].position.z;
    BigReal z3 = p[2]->x[localIndex[2]].position.z;
    BigReal z4 = p[3]->x[localIndex[3]].position.z;
    int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
    int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
    int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
    int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
    pp_clamp(n1, pressureProfileSlabs);
    pp_clamp(n2, pressureProfileSlabs);
    pp_clamp(n3, pressureProfileSlabs);
    pp_clamp(n4, pressureProfileSlabs);
    int p1 = p[0]->x[localIndex[0]].partition;
    int p2 = p[1]->x[localIndex[1]].partition;
    int p3 = p[2]->x[localIndex[2]].partition;
    int p4 = p[3]->x[localIndex[3]].partition;
    int pn = pressureProfileAtomTypes;
    pp_reduction(pressureProfileSlabs, n1, n2,
                p1, p2, pn,
                f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
                pressureProfileData);
    pp_reduction(pressureProfileSlabs, n2, n3,
                p2, p3, pn,
                f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
                pressureProfileData);
    pp_reduction(pressureProfileSlabs, n3, n4,
                p3, p4, pn,
                f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
                pressureProfileData);
  }

 }
}

void DihedralElem::getMoleculePointers ( Molecule *  mol, int *  count, int32 ***  byatom, Dihedral **  structarray  )

static

Definition at line 25 of file ComputeDihedrals.C.

References NAMD_die(), and Molecule::numDihedrals.

{
#ifdef MEM_OPT_VERSION
  NAMD_die("Should not be called in DihedralElem::getMoleculePointers in memory optimized version!");
#else
  *count = mol->numDihedrals;
  *byatom = mol->dihedralsByAtom;
  *structarray = mol->dihedrals;
#endif
}

void DihedralElem::getParameterPointers ( Parameters *  p, const DihedralValue **  v  )

static

Definition at line 36 of file ComputeDihedrals.C.

References Parameters::dihedral_array.

                                                                              {
  *v = p->dihedral_array;
}

static void DihedralElem::getTupleInfo ( AtomSignature *  sig, int *  count, TupleSignature **  t  )

inlinestatic

Definition at line 29 of file ComputeDihedrals.h.

References AtomSignature::dihedralCnt, and AtomSignature::dihedralSigs.

                                                                                 {
        *count = sig->dihedralCnt;
        *t = sig->dihedralSigs;
    }

int DihedralElem::hash ( void  ) const

inline

Definition at line 43 of file ComputeDihedrals.h.

References atomID.

                   { 
    return 0x7FFFFFFF &((atomID[0]<<24) + (atomID[1]<<16) + (atomID[2]<<8) + atomID[3]);
  }

int DihedralElem::operator< ( const DihedralElem &  a ) const

inline

Definition at line 50 of file ComputeDihedrals.inl.

References atomID.

  {
    return  (atomID[0] < a.atomID[0] ||
            (atomID[0] == a.atomID[0] &&
            (atomID[1] < a.atomID[1] ||
            (atomID[1] == a.atomID[1] &&
            (atomID[2] < a.atomID[2] ||
            (atomID[2] == a.atomID[2] &&
             atomID[3] < a.atomID[3] 
             ))))));
  }

int DihedralElem::operator== ( const DihedralElem &  a ) const

inline

Definition at line 44 of file ComputeDihedrals.inl.

References atomID.

  {
    return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
        a.atomID[2] == atomID[2] && a.atomID[3] == atomID[3]);
  }

void DihedralElem::submitReductionData ( BigReal *  data, SubmitReduction *  reduction  )

static

Definition at line 328 of file ComputeDihedrals.C.

References ADD_TENSOR, dihedralEnergyIndex, dihedralEnergyIndex_f, dihedralEnergyIndex_ti_1, dihedralEnergyIndex_ti_2, SubmitReduction::item(), REDUCTION_BONDED_ENERGY_F, REDUCTION_BONDED_ENERGY_TI_1, REDUCTION_BONDED_ENERGY_TI_2, and REDUCTION_DIHEDRAL_ENERGY.

{
  reduction->item(REDUCTION_DIHEDRAL_ENERGY) += data[dihedralEnergyIndex];
  reduction->item(REDUCTION_BONDED_ENERGY_F) += data[dihedralEnergyIndex_f];
  reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[dihedralEnergyIndex_ti_1];
  reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[dihedralEnergyIndex_ti_2];
  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
  ADD_TENSOR(reduction,REDUCTION_VIRIAL_AMD_DIHE,data,virialIndex);
}

Member Data Documentation

AtomID DihedralElem::atomID[size]

Definition at line 21 of file ComputeDihedrals.h.

Referenced by computeForce(), DihedralElem(), hash(), operator<(), and operator==().

int DihedralElem::localIndex[size]

Definition at line 22 of file ComputeDihedrals.h.

Referenced by computeForce().

TuplePatchElem* DihedralElem::p[size]

Definition at line 23 of file ComputeDihedrals.h.

Referenced by computeForce().

int DihedralElem::pressureProfileAtomTypes = 1

static

Definition at line 36 of file ComputeDihedrals.h.

Referenced by computeForce().

BigReal DihedralElem::pressureProfileMin = 0

static

Definition at line 38 of file ComputeDihedrals.h.

Referenced by computeForce().

int DihedralElem::pressureProfileSlabs = 0

static

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 35 of file ComputeDihedrals.h.

Referenced by computeForce().

BigReal DihedralElem::pressureProfileThickness = 0

static

Definition at line 37 of file ComputeDihedrals.h.

Referenced by computeForce().

Real DihedralElem::scale

Definition at line 24 of file ComputeDihedrals.h.

Referenced by computeForce().

const DihedralValue* DihedralElem::value

Definition at line 41 of file ComputeDihedrals.h.

Referenced by computeForce(), and DihedralElem().

---

The documentation for this class was generated from the following files:

• ComputeDihedrals.h
• ComputeDihedrals.C
• ComputeDihedrals.inl