NAMD
Public Member Functions | Protected Member Functions | List of all members
ADiheRestraint Class Referenceabstract

#include <FreeEnergyRestrain.h>

Inheritance diagram for ADiheRestraint:
ARestraint ABoundDiheRestraint AFixedDiheRestraint AForcingDiheRestraint

Public Member Functions

 ADiheRestraint ()
 
void PrintInfo ()
 
virtual Bool_t TwoTargets ()=0
 
virtual double GetDiheTarget1 ()=0
 
virtual double GetDiheTarget2 ()=0
 
- Public Member Functions inherited from ARestraint
 ARestraint ()
 
virtual ~ARestraint ()
 
int GetNumGroups ()
 
void SetKf (double Kf)
 
double GetKf ()
 
void SetLambdaKf (double LambdaKf)
 
void SetLambdaRef (double LambdaRef)
 
double GetLambdaKf ()
 
double GetLambdaRef ()
 
void SetGroup (AGroup &Group, int GroupIndex)
 
void SetGroups (AGroup &Group1)
 
void SetGroups (AGroup &Group1, AGroup &Group2)
 
void SetGroups (AGroup &Group1, AGroup &Group2, AGroup &Group3)
 
void SetGroups (AGroup &Group1, AGroup &Group2, AGroup &Group3, AGroup &Group4)
 
void UpdateCOMs (GlobalMasterFreeEnergy &CFE)
 
void DistributeForce (int WhichGroup, AVector Force, GlobalMasterFreeEnergy &CFE)
 
virtual AVector GetGradient (int WhichGroup)=0
 
virtual double GetEnergy ()=0
 
virtual void GetStr (char *Str)=0
 
virtual Bool_t IsForcing ()
 
virtual double Get_dU_dLambda ()
 
virtual void SetRefPos (AVector)
 
virtual void SetRefDist (double)
 
virtual void SetRefAngle (double)
 
virtual void SetBound (Bound_t)
 
virtual void SetLowerAngle (double)
 
virtual void SetUpperAngle (double)
 
virtual void SetIntervalAngle (double)
 
virtual void SetStartPos (AVector)
 
virtual void SetStopPos (AVector)
 
virtual void SetStartDist (double)
 
virtual void SetStopDist (double)
 
virtual void SetStartAngle (double)
 
virtual void SetStopAngle (double)
 

Protected Member Functions

double GetE (double RefDihe, double Const)
 
AVector GetGrad (int WhichGroup, double RefDihe, double Const)
 
AVector gradU (AVector &P1P2P3, AVector &P4P5P6, AVector &dP1, AVector &dP2, AVector &dP3, AVector &dP4, AVector &dP5, AVector &dP6)
 
- Protected Member Functions inherited from ARestraint
double GetAngle (AVector &A, AVector &B, AVector &C)
 
double GetDihe (AVector &A, AVector &B, AVector &C, AVector &D)
 
void EarlyExit (const char *Str, int AtomID)
 

Additional Inherited Members

- Protected Attributes inherited from ARestraint
double m_Kf
 
int m_NumGroups
 
AGroupm_pGroups
 
AVectorm_pCOMs
 
- Static Protected Attributes inherited from ARestraint
static double m_LambdaKf = 1.0
 
static double m_LambdaRef = 0.0
 

Detailed Description

Definition at line 146 of file FreeEnergyRestrain.h.

Constructor & Destructor Documentation

ADiheRestraint::ADiheRestraint ( )

Definition at line 537 of file FreeEnergyRestrain.C.

References ARestraint::m_NumGroups, ARestraint::m_pCOMs, and ARestraint::m_pGroups.

537  {
538 //----------------------------------------------------------------------------
539 // each ADiheRestraint restrains the dihedral angle between 4 groups of atoms
540 //----------------------------------------------------------------------------
541  m_NumGroups = 4;
543  m_pCOMs = new AVector[m_NumGroups];
544 }
AGroup * m_pGroups
AVector * m_pCOMs

Member Function Documentation

virtual double ADiheRestraint::GetDiheTarget1 ( )
pure virtual
virtual double ADiheRestraint::GetDiheTarget2 ( )
pure virtual
double ADiheRestraint::GetE ( double  RefDihe,
double  Const 
)
protected

Definition at line 587 of file FreeEnergyRestrain.C.

References ARestraint::GetDihe(), and ARestraint::m_pCOMs.

Referenced by AFixedDiheRestraint::GetEnergy(), ABoundDiheRestraint::GetEnergy(), and AForcingDiheRestraint::GetEnergy().

587  {
588 //---------------------------------------------------------------------------
589 // calculate and return the Energy for this angle restraint.
590 //
591 // E = (E0/2) * (1 - cos(Chi - ChiRef))
592 //
593 // where Chi is the dihedral angle between 4 centers-of-mass of restrained atoms,
594 // m_pCOMs[0] -- m_pCOMs[1] -- m_pCOMs[2] -- m_pCOMs[3].
595 // ChiRef is the reference angle.
596 //
597 // Note: COM's are calculated before this routine is called.
598 //---------------------------------------------------------------------------
599  double Angle;
600 
601  Angle = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
602  return( (Const/2.0) * (1.0 - cos(Angle-RefAngle)) );
603 }
struct angle Angle
double GetDihe(AVector &A, AVector &B, AVector &C, AVector &D)
AVector * m_pCOMs
AVector ADiheRestraint::GetGrad ( int  WhichGroup,
double  RefDihe,
double  Const 
)
protected

Definition at line 606 of file FreeEnergyRestrain.C.

References A, ASSERT, B, AVector::cross(), AVector::Dist(), AVector::dot(), ARestraint::GetDihe(), gradU(), kAlmostMinusOne, kAlmostOne, ARestraint::m_pCOMs, and AVector::Set().

Referenced by AFixedDiheRestraint::GetGradient(), ABoundDiheRestraint::GetGradient(), and AForcingDiheRestraint::GetGradient().

607  {
608 //---------------------------------------------------------------------------
609 // calculate and return the gradient for this dihedral angle restraint.
610 //
611 // E = (E0/2) * (1 - cos(Chi - ChiRef))
612 //
613 // return: grad(E)
614 //
615 // Notes: COM's are calculated before this routine is called.
616 //---------------------------------------------------------------------------
617  AVector A, B, C, D;
618 
619  ASSERT((WhichGroup==0)||(WhichGroup==1)||(WhichGroup==2)||(WhichGroup==3));
620 
621  if ((WhichGroup==0) || (WhichGroup==1)) {
622  A = m_pCOMs[0];
623  B = m_pCOMs[1];
624  C = m_pCOMs[2];
625  D = m_pCOMs[3];
626  }
627  // re-state the problem so the gradient is solved for either atoms 0 or 1
628  else {
629  A = m_pCOMs[3];
630  B = m_pCOMs[2];
631  C = m_pCOMs[1];
632  D = m_pCOMs[0];
633  if (WhichGroup==3) {WhichGroup=0;}
634  if (WhichGroup==2) {WhichGroup=1;}
635  }
636 
637  AVector CD(D - C);
638  AVector CB(B - C);
639  AVector BC(C - B);
640  AVector BA(A - B);
641  AVector AC(C - A);
642  AVector CDxCB, BCxBA;
643  AVector Vec;
644  double phi;
645  double top, bot, u;
646 
647  CDxCB = CD.cross(CB);
648  BCxBA = BC.cross(BA);
649 
650  top = CDxCB.dot(BCxBA);
651  bot = CDxCB.Dist() * BCxBA.Dist();
652 
653  u = top/bot;
654  // protect against acos(<-1.0), acos(>1.0), sqrt(<0), and divide by 0
655  if (u < kAlmostMinusOne) {u = kAlmostMinusOne;}
656  if (u > kAlmostOne) {u = kAlmostOne;}
657 
658  // get dihedral using atan
659  phi = GetDihe(A,B,C,D);
660 
661  AVector dP1, dP2, dP3;
662  AVector dP4, dP5, dP6;
663  ASSERT((WhichGroup==0) || (WhichGroup==1));
664  if (WhichGroup==0) {
665  dP1.Set( 0, 0, 0 );
666  dP2.Set( 0, 0, 0 );
667  dP3.Set( 0, 0, 0 );
668  dP4.Set( 0, -BC[2], BC[1]);
669  dP5.Set( BC[2], 0, -BC[0]);
670  dP6.Set(-BC[1], BC[0], 0 );
671  }
672  else {
673  dP1.Set( 0, -CD[2], CD[1]);
674  dP2.Set( CD[2], 0, -CD[0]);
675  dP3.Set(-CD[1], CD[0], 0 );
676  dP4.Set( 0, AC[2], -AC[1]);
677  dP5.Set(-AC[2], 0, AC[0]);
678  dP6.Set( AC[1], -AC[0], 0 );
679  }
680 
681  Vec = gradU(CDxCB, BCxBA, dP1, dP2, dP3, dP4, dP5, dP6);
682  Vec *= (Const/2.0) * sin(phi-RefAngle) * (-1.0/sqrt(1.0 - u*u));
683 
684  // flip gradient for negative angles
685  if (phi < 0) {
686  Vec *= -1.0;
687  }
688 
689  return(Vec);
690 }
Bool_t Set(double x, double y, double z)
const double kAlmostMinusOne
AVector cross(const AVector &Vector)
AVector gradU(AVector &P1P2P3, AVector &P4P5P6, AVector &dP1, AVector &dP2, AVector &dP3, AVector &dP4, AVector &dP5, AVector &dP6)
const BigReal A
double GetDihe(AVector &A, AVector &B, AVector &C, AVector &D)
double dot(const AVector &Vector)
#define ASSERT(E)
const BigReal B
double Dist()
const double kAlmostOne
AVector * m_pCOMs
AVector ADiheRestraint::gradU ( AVector P1P2P3,
AVector P4P5P6,
AVector dP1,
AVector dP2,
AVector dP3,
AVector dP4,
AVector dP5,
AVector dP6 
)
protected

Definition at line 693 of file FreeEnergyRestrain.C.

References AVector::Dist(), and AVector::dot().

Referenced by GetGrad().

695  {
696 //----------------------------------------------------------------
697 // calculate the gradient for ((P1P2P3.dot.P4P5P6)/(mag(P1P2P3)*mag(P4P5P6)))
698 // P1P2P3 = (P1)i + (P2)j + (P3)k
699 // P4P5P6 = (P4)i + (P5)j + (P6)k
700 // dP1 = (d(P1)/dx)i + (d(P1)/dy)j +(d(P1)/dz)k
701 // dP2 = (d(P2)/dx)i + (d(P2)/dy)j +(d(P2)/dz)k
702 // dP3 = (d(P3)/dx)i + (d(P3)/dy)j +(d(P3)/dz)k
703 // dP4 = (d(P4)/dx)i + (d(P4)/dy)j +(d(P4)/dz)k
704 // dP5 = (d(P5)/dx)i + (d(P5)/dy)j +(d(P5)/dz)k
705 // dP6 = (d(P6)/dx)i + (d(P6)/dy)j +(d(P6)/dz)k
706 //----------------------------------------------------------------
707  double Mag123, Mag456, Dot;
708  double Const1, Const2, Const3;
709  double P1, P2, P3, P4, P5, P6;
710  AVector RetVec;
711 
712  P1 = P1P2P3[0]; P2 = P1P2P3[1]; P3 = P1P2P3[2];
713  P4 = P4P5P6[0]; P5 = P4P5P6[1]; P6 = P4P5P6[2];
714 
715  Mag123 = P1P2P3.Dist();
716  Mag456 = P4P5P6.Dist();
717  Dot = P1P2P3.dot(P4P5P6);
718 
719  Const1 = 1.0 / (Mag123*Mag456);
720  Const2 = -Dot * (1.0 / (Mag123*Mag456*Mag456*Mag456));
721  Const3 = -Dot * (1.0 / (Mag456*Mag123*Mag123*Mag123));
722 
723  RetVec = (dP4*P1 + dP1*P4 + dP5*P2 + dP2*P5 + dP6*P3 + dP3*P6) * Const1 +
724  (dP4*P4 + dP5*P5 + dP6*P6) * Const2 +
725  (dP1*P1 + dP2*P2 + dP3*P3) * Const3;
726 
727  return(RetVec);
728 }
double dot(const AVector &Vector)
double Dist()
void ADiheRestraint::PrintInfo ( )
virtual

Implements ARestraint.

Definition at line 547 of file FreeEnergyRestrain.C.

References endi(), ARestraint::GetDihe(), GetDiheTarget1(), GetDiheTarget2(), iout, kPi, ARestraint::m_pCOMs, and TwoTargets().

547  {
548 //--------------------------------------------------------------------
549 // print the dihedral angle for this dihedral restraint
550 //--------------------------------------------------------------------
551  double Dihedral;
552  char Str1[20], Str2[20], Str3[20];
553 
554  Dihedral = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]) * (180/kPi);
555  sprintf(Str1, "%8.3f", Dihedral);
556  Dihedral = GetDiheTarget1() * (180/kPi);
557  sprintf(Str2, "%8.3f", Dihedral);
558  Dihedral = GetDiheTarget2() * (180/kPi);
559  sprintf(Str3, "%8.3f", Dihedral);
560 
561 #if defined(_VERBOSE_PMF)
562  iout << "Dihedral = ";
563  iout << Str1;
564  iout << " degrees";
565  iout << " Target = ";
566  iout << Str2;
567  iout << " degrees";
568  if (TwoTargets()) {
569  iout << " to ";
570  iout << Str3;
571  iout << " degrees";
572  }
573  iout << std::endl << endi;
574 #else
575  iout << Str1;
576  iout << " ";
577  iout << Str2;
578  if (TwoTargets()) {
579  iout << ", ";
580  iout << Str3;
581  }
582  iout << " | ";
583 #endif
584 }
struct dihedral Dihedral
double GetDihe(AVector &A, AVector &B, AVector &C, AVector &D)
#define iout
Definition: InfoStream.h:87
virtual double GetDiheTarget1()=0
virtual double GetDiheTarget2()=0
virtual Bool_t TwoTargets()=0
const double kPi
infostream & endi(infostream &s)
Definition: InfoStream.C:38
AVector * m_pCOMs
virtual Bool_t ADiheRestraint::TwoTargets ( )
pure virtual

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