FreeEnergyRestrain.C

Go to the documentation of this file.
00001 
00007 // written by David Hurwitz, March to May 1998.
00008 
00009 #include <string.h>
00010 //#include <iomanip.h>
00011 #include "common.h"
00012 #include "InfoStream.h"
00013 #include "FreeEnergyAssert.h"
00014 #include "FreeEnergyEnums.h"
00015 #include "Vector.h"
00016 #include "FreeEnergyVector.h"
00017 #include "FreeEnergyGroup.h"
00018 #include "FreeEnergyRestrain.h"
00019 #include "FreeEnergyRMgr.h"
00020 #include "FreeEnergyLambda.h"
00021 #include "FreeEnergyLambdMgr.h"
00022 
00023 #include "NamdTypes.h"
00024 #include "GlobalMaster.h"
00025 #include "GlobalMasterFreeEnergy.h"
00026 
00027 // initialize static member variables
00028 // (lambda is the same for all forcing restraints)
00029 double  ARestraint::m_LambdaKf  = 1.0;
00030 double  ARestraint::m_LambdaRef = 0.0;
00031 
00032 #if defined(_DEBUG)
00033 void ARestraint::SetCOM(int Index, AVector& Pos) {
00034 //-----------------------------------------------------------------
00035 // only for testing!!!!!!
00036 //-----------------------------------------------------------------
00037   ASSERT( (Index>=0) && (Index<m_NumGroups) );
00038   m_pCOMs[Index] = Pos;
00039 }
00040 #endif
00041 
00042 ARestraint::ARestraint() {
00043 //-----------------------------------------------------------------
00044 // constructor for base class
00045 //-----------------------------------------------------------------
00046   m_pGroups = NULL;
00047   m_pCOMs = NULL;
00048   m_NumGroups = 0;
00049 }
00050 
00051 
00052 ARestraint::~ARestraint() {
00053 //-----------------------------------------------------------------
00054 // free space that may have been allocated for Groups and COM's
00055 //-----------------------------------------------------------------
00056   if (m_pGroups != NULL) {
00057     ASSERT(m_pCOMs != NULL);
00058     delete []m_pGroups;
00059     delete []m_pCOMs;
00060   }
00061 }
00062 
00063 
00064 void ARestraint::EarlyExit(const char* Str, int AtomID) {
00065 //-----------------------------------------------------------------
00066 // unrecoverable error
00067 //-----------------------------------------------------------------
00068   iout << "FreeEnergy: " << std::endl << endi;
00069   iout << "FreeEnergy: " << Str;
00070   iout << " for AtomID: " << AtomID;
00071   iout << std::endl << endi;
00072   NAMD_die("FreeEnergy: Fatal Error with Fixed or Forcing Restraints");
00073 }
00074 
00075 
00076 void ARestraint::SetGroup(AGroup& Group, int GroupIndex) {
00077 //-----------------------------------------------------------------
00078 // set one group of atoms
00079 //-----------------------------------------------------------------
00080   ASSERT( (GroupIndex>=0) && (GroupIndex<m_NumGroups) );
00081   m_pGroups[GroupIndex] = Group;
00082 }
00083 
00084 
00085 void ARestraint::SetGroups(AGroup& Group1) {
00086 //-----------------------------------------------------------------
00087 // set one group of atoms
00088 //-----------------------------------------------------------------
00089   ASSERT(m_NumGroups >= 1);
00090   m_pGroups[0] = Group1;
00091 }
00092 
00093 
00094 void ARestraint::SetGroups(AGroup& Group1, AGroup& Group2) {
00095 //-----------------------------------------------------------------
00096 // set two groups of atoms
00097 //-----------------------------------------------------------------
00098   ASSERT(m_NumGroups >= 2);
00099   m_pGroups[0] = Group1;
00100   m_pGroups[1] = Group2;
00101 }
00102 
00103 
00104 void ARestraint::SetGroups(AGroup& Group1, AGroup& Group2, AGroup& Group3) {
00105 //-----------------------------------------------------------------
00106 // set three groups of atoms
00107 //-----------------------------------------------------------------
00108   ASSERT(m_NumGroups >= 3);
00109   m_pGroups[0] = Group1;
00110   m_pGroups[1] = Group2;
00111   m_pGroups[2] = Group3;
00112 }
00113 
00114 
00115 void ARestraint::SetGroups(AGroup& Group1, AGroup& Group2, AGroup& Group3, AGroup& Group4) {
00116 //-----------------------------------------------------------------
00117 // set four groups of atoms
00118 //-----------------------------------------------------------------
00119   ASSERT(m_NumGroups >= 4);
00120   m_pGroups[0] = Group1;
00121   m_pGroups[1] = Group2;
00122   m_pGroups[2] = Group3;
00123   m_pGroups[3] = Group4;
00124 }
00125 
00126   
00127 double ARestraint::GetAngle(AVector& A, AVector& B, AVector& C) {
00128 //-----------------------------------------------------------------
00129 // determine the angle formed by the points A-B-C
00130 //-----------------------------------------------------------------
00131   double u;
00132   double a = B.Dist(C);
00133   double b = A.Dist(C);
00134   double c = A.Dist(B);
00135 
00136   u = (a*a + c*c - b*b) / (2.0*a*c);
00137   // protect against acos(<-1.0) and acos(>1.0)
00138   if (u < -1.0) {u = -1.0;}
00139   if (u >  1.0) {u =  1.0;}
00140   return(acos(u));
00141 }
00142 
00143 
00144 double ARestraint::GetDihe(AVector& A, AVector& B, AVector& C, AVector& D) {
00145 //-----------------------------------------------------------------
00146 // determine the dihedral angle formed by the points A-B-C-D
00147 //-----------------------------------------------------------------
00148   AVector CD(D - C);
00149   AVector CB(B - C);
00150   AVector BC(C - B);
00151   AVector BA(A - B);
00152   AVector CDxCB, BCxBA;
00153   double  top, bot, cos_u, sin_u, Angle;
00154   AVector topVec;
00155 
00156   CDxCB = CD.cross(CB);
00157   BCxBA = BC.cross(BA);
00158 
00159   top = CDxCB.dot(BCxBA);
00160   bot = CDxCB.Dist() * BCxBA.Dist();
00161   cos_u = top/bot;
00162 
00163   // protect against acos(<-1.0) and acos(>1.0)
00164   if (cos_u < -1.0) {cos_u = -1.0;}
00165   if (cos_u >  1.0) {cos_u =  1.0;}
00166 
00167   topVec = CDxCB.cross(BCxBA);
00168   sin_u = (topVec/bot).dot(CB/CB.Dist());
00169 
00170   // protect against asin(<-1.0) and asin(>1.0)
00171   if (sin_u < -1.0) {sin_u = -1.0;}
00172   if (sin_u >  1.0) {sin_u =  1.0;}
00173 
00174   Angle = atan2(sin_u, cos_u);
00175   return(Angle);
00176 }
00177 
00178 
00179 void ARestraint::DistributeForce(int WhichGroup, AVector Force,
00180                                  GlobalMasterFreeEnergy& CFE) {
00181 //----------------------------------------------------------------------
00182 // Distribute Force among the group of atoms specified by WhichGroup
00183 //
00184 // note:  m_pGroups points to an array of Groups
00185 //        m_pGroups[WhichGroup] references one of the Groups
00186 //        m_pGroups[WhichGroup][i] returns an AtomID from the Group
00187 //        (operator[] is defined to return an item from the Group)
00188 //----------------------------------------------------------------------
00189   int     i, AtomID, NumAtoms, RetVal;
00190   double  Mass, TotalMass=0;
00191   AVector SmallForce;
00192   Vector  NAMD_Vector;
00193 
00194   ASSERT( (WhichGroup>=0) && (WhichGroup<m_NumGroups) );
00195 
00196   // calculate the total mass for the group
00197   NumAtoms = m_pGroups[WhichGroup].GetNumInGroup();
00198   for (i=0; i<NumAtoms; i++) {
00199     AtomID = m_pGroups[WhichGroup][i];
00200     Mass = CFE.getMass(AtomID);
00201     if (Mass < 0) {EarlyExit("Negative Mass", AtomID);};
00202     TotalMass += Mass;
00203   }
00204 
00205   // distribute Force according to mass of each atom in the group
00206   for (i=0; i<NumAtoms; i++) {
00207     AtomID = m_pGroups[WhichGroup][i];
00208     Mass = CFE.getMass(AtomID);
00209     if (Mass < 0) {EarlyExit("Negative Mass", AtomID);}
00210     SmallForce = Force * (Mass/TotalMass);
00211     // cast SmallForce to a NAMD-type vector (addForce uses Vector class)
00212     SetEqual(NAMD_Vector, SmallForce);
00213     RetVal = CFE.addForce(AtomID, NAMD_Vector);
00214     if (RetVal < 0) {EarlyExit("Can't add Force", AtomID);}
00215   }
00216 }
00217 
00218 
00219 void ARestraint::UpdateCOMs(GlobalMasterFreeEnergy& CFE) {
00220 //-----------------------------------------------------------------
00221 // calculate the center-of-mass of each group of atoms
00222 //
00223 // note:  m_pGroups points to an array of Groups
00224 //        m_pGroups[i] references one of the Groups
00225 //        m_pGroups[i][j] returns an AtomID from the Group
00226 //        (operator[] is defined to return an item from the Group)
00227 //-----------------------------------------------------------------
00228   int      i, j, AtomID, RetVal;
00229   Vector   NAMD_Vector;
00230   AVector  COM, Pos;
00231   double   Mass, TotalMass;
00232 
00233   ASSERT(m_NumGroups > 0);
00234   // for each group of atoms
00235   for (i=0; i<m_NumGroups; i++) {
00236     TotalMass = 0;
00237     COM.Set(0,0,0);
00238     // for each atom in the group
00239     for (j=0; j<m_pGroups[i].GetNumInGroup(); j++) {
00240       AtomID = m_pGroups[i][j];
00241       // get its position, weight position with atom's mass
00242       RetVal = CFE.getPosition(AtomID, NAMD_Vector);
00243       if (RetVal < 0) {EarlyExit("Can't get Position", AtomID);}
00244       // cast NAMD_Vector to AVector (getPosition uses Vector class)
00245       SetEqual(Pos, NAMD_Vector);
00246       Mass = CFE.getMass(AtomID);
00247       if (Mass < 0) {EarlyExit("Negative Mass", AtomID);}
00248       TotalMass += Mass;
00249       COM += Pos * Mass;
00250     }
00251     m_pCOMs[i] = COM / TotalMass;
00252   }
00253 }
00254 
00255 
00256 APosRestraint::APosRestraint() {
00257 //-----------------------------------------------------------------
00258 // each APosRestraint restrains 1 group of atoms to a location
00259 //-----------------------------------------------------------------
00260   m_NumGroups = 1;
00261   m_pGroups = new AGroup[m_NumGroups];
00262   m_pCOMs = new AVector[m_NumGroups];
00263 }
00264 
00265 
00266 void APosRestraint::PrintInfo() {
00267 //--------------------------------------------------------------------
00268 // print the position for this position restraint
00269 //--------------------------------------------------------------------
00270   double  Distance;
00271   char    Str[20];
00272 
00273   Distance = GetDistance();
00274   sprintf(Str, "%7.3f", Distance);
00275 
00276 #if defined(_VERBOSE_PMF)
00277   iout << "Position = ";
00278   m_pCOMs[0].Out();
00279   iout << "  Target = ";
00280   GetPosTarget().Out();
00281   iout << "  Distance = ";
00282   iout << Str;
00283   iout << std::endl << endi;
00284 #else
00285   m_pCOMs[0].Out();
00286   iout << "  ";
00287   GetPosTarget().Out();
00288   iout << "  ";
00289   iout << Str;
00290   iout << " | ";
00291 #endif
00292 }
00293 
00294 
00295 double APosRestraint::GetE(AVector RefPos, double LambdaKf) {
00296 //--------------------------------------------------------------------
00297 // calculate and return the Energy for this position restraint.
00298 //
00299 //     E = (Kf/2) * (|ri - rref|)**2
00300 //
00301 // where |ri - rref| is the distance between a) the center-of-mass
00302 // of the restrained atoms and b) the reference position.
00303 //
00304 // Note:  COM is calculated before this routine is called.
00305 //--------------------------------------------------------------------
00306   return( ((m_Kf*LambdaKf)/2.0) * m_pCOMs[0].DistSqr(RefPos) );
00307 }
00308 
00309 
00310 AVector APosRestraint::GetGrad(int /* WhichGroup */, 
00311                                AVector RefPos, double LambdaKf) {
00312 //-------------------------------------------------------------------------
00313 // calculate and return the gradient for this position restraint.
00314 //
00315 //     E = (Kf/2) * (|ri - rref|)**2
00316 //
00317 // return:  grad(E)
00318 //
00319 // Notes: COM is calculated before this routine is called.
00320 //        m_pCOMs points to an array of vectors
00321 //        m_pCOMS[0] references the only COM for a position restraint
00322 //        m_pCOMS[0][0] returns the x value from the vector (x,y,z)
00323 //        (operator[] is defined to return an item from the vector)
00324 //-------------------------------------------------------------------------
00325   // WhichGroup = 0;  // don't care -- there's only 1 atom restrained
00326   AVector Vec(m_pCOMs[0][0] - RefPos[0],
00327               m_pCOMs[0][1] - RefPos[1],
00328               m_pCOMs[0][2] - RefPos[2]);
00329   return(Vec*m_Kf*LambdaKf);
00330 }
00331 
00332 
00333 ADistRestraint::ADistRestraint() {
00334 //-----------------------------------------------------------------------
00335 // each ADistRestraint restrains the distance between 2 groups of atoms
00336 //-----------------------------------------------------------------------
00337   m_NumGroups = 2;
00338   m_pGroups = new AGroup[m_NumGroups];
00339   m_pCOMs = new AVector[m_NumGroups];
00340 }
00341 
00342 
00343 void ADistRestraint::PrintInfo() {
00344 //--------------------------------------------------------------------
00345 // print the distance for this distance restraint
00346 //--------------------------------------------------------------------
00347   double  Distance;
00348   char    Str1[20], Str2[20];
00349 
00350   Distance = m_pCOMs[0].Dist(m_pCOMs[1]);
00351   sprintf(Str1, "%7.3f", Distance);
00352   Distance = GetDistTarget();
00353   sprintf(Str2, "%7.3f", Distance);
00354 
00355 #if defined(_VERBOSE_PMF)
00356   iout << "Distance = ";
00357   iout << Str1;
00358   iout << "  Target = ";
00359   iout << Str2;
00360   iout << std::endl << endi;
00361 #else
00362   iout << Str1;
00363   iout << "  ";
00364   iout << Str2;
00365   iout << " | ";
00366 #endif
00367 }
00368 
00369 
00370 double ADistRestraint::GetE(double RefDist, double LambdaKf) {
00371 //---------------------------------------------------------------------------
00372 // calculate and return the Energy for this distance restraint.
00373 //
00374 //     E = (Kf/2) * (di-dref)**2
00375 //
00376 // where di is the distance between 2 centers-of-mass of restrained atoms,
00377 // and dref is the reference distance.
00378 //
00379 // Note:  COM's are calculated before this routine is called.
00380 //---------------------------------------------------------------------------
00381   double Dist, Diff;
00382 
00383   Dist = m_pCOMs[0].Dist(m_pCOMs[1]);
00384   Diff = Dist - RefDist;
00385   return( ((m_Kf*LambdaKf)/2.0) * (Diff*Diff) );
00386 }
00387 
00388 
00389 AVector ADistRestraint::GetGrad(int WhichGroup,
00390                                 double RefDist, double LambdaKf) {
00391 //---------------------------------------------------------------------------
00392 // calculate and return the gradient for this distance restraint.
00393 //
00394 //     E = (Kf/2) * (di-dref)**2
00395 //
00396 // return:  grad(E)
00397 //
00398 // Notes: COM is calculated before this routine is called.
00399 //        m_pCOMS[0 & 1] reference the COM's of each group of atoms
00400 //---------------------------------------------------------------------------
00401   double  Dist;
00402   AVector Vec, A, B;
00403 
00404   ASSERT( (WhichGroup==0) || (WhichGroup==1) );
00405   if (WhichGroup == 0) {
00406     A = m_pCOMs[0];
00407     B = m_pCOMs[1];
00408   }
00409   else {
00410     A = m_pCOMs[1];
00411     B = m_pCOMs[0];
00412   }
00413 
00414   Dist = A.Dist(B);
00415   Vec.Set(A[0]-B[0], A[1]-B[1], A[2]-B[2]);
00416   Vec *= m_Kf * LambdaKf * (Dist - RefDist) / Dist;
00417   return(Vec);
00418 }
00419 
00420 
00421 AnAngleRestraint::AnAngleRestraint() {
00422 //-----------------------------------------------------------------------
00423 // each AnAngleRestraint restrains the angle between 3 groups of atoms
00424 //-----------------------------------------------------------------------
00425   m_NumGroups = 3;
00426   m_pGroups = new AGroup[m_NumGroups];
00427   m_pCOMs = new AVector[m_NumGroups];
00428 }
00429 
00430 
00431 void AnAngleRestraint::PrintInfo() {
00432 //--------------------------------------------------------------------
00433 // print the angle for this angle restraint
00434 //--------------------------------------------------------------------
00435   double  Angle;
00436   char    Str1[20], Str2[20];
00437 
00438   Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]) * (180/kPi);
00439   sprintf(Str1, "%8.3f", Angle);
00440   Angle = GetAngleTarget() * (180/kPi);
00441   sprintf(Str2, "%8.3f", Angle);
00442 
00443 #if defined(_VERBOSE_PMF)
00444   iout << "Angle = ";
00445   iout << Str1;
00446   iout << " degrees";
00447   iout << "  Target = ";
00448   iout << Str2;
00449   iout << " degrees";
00450   iout << std::endl << endi;
00451 #else
00452   iout << Str1;
00453   iout << "  ";
00454   iout << Str2;
00455   iout << " | ";
00456 #endif
00457 }
00458 
00459 
00460 double AnAngleRestraint::GetE(double RefAngle, double LambdaKf) {
00461 //---------------------------------------------------------------------------
00462 // calculate and return the Energy for this angle restraint.
00463 //
00464 //     E = (Kf/2) * (Theta-ThetaRef)**2
00465 //
00466 // where Theta is the angle between 3 centers-of-mass of restrained atoms,
00467 // m_pCOMs[0] -- m_pCOMs[1] -- m_pCOMs[2].
00468 // ThetaRef is the reference angle.
00469 //
00470 // Note:  COM's are calculated before this routine is called.
00471 //---------------------------------------------------------------------------
00472   double Angle, Diff;
00473 
00474   Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]);
00475   Diff = Angle - RefAngle;
00476   return( ((m_Kf*LambdaKf)/2.0) * (Diff*Diff) );
00477 }
00478 
00479 
00480 AVector AnAngleRestraint::GetGrad(int WhichGroup,
00481                                   double RefAngle, double LambdaKf) {
00482 //---------------------------------------------------------------------------
00483 // calculate and return the gradient for this angle restraint.
00484 //
00485 //     E = (Kf/2) * (Theta-ThetaRef)**2
00486 //
00487 // return:  grad(E)
00488 //
00489 // Notes: COM's are calculated before this routine is called.
00490 //---------------------------------------------------------------------------
00491   AVector A, B, C;
00492   double  Angle;
00493   double  a, b, c;
00494   double  u;
00495   double  Const1, Const2, Const3, Const4, Const5;
00496   AVector Vec1, Vec2;
00497 
00498   ASSERT( (WhichGroup==0) || (WhichGroup==1) || (WhichGroup==2) );
00499   A = m_pCOMs[0];
00500   B = m_pCOMs[1];
00501   C = m_pCOMs[2];
00502 
00503   a = B.Dist(C);
00504   b = A.Dist(C);
00505   c = A.Dist(B);
00506 
00507   u = (a*a + c*c - b*b) / (2.0*a*c);
00508 
00509   // protect against acos(<-1.0), acos(>1.0), sqrt(<0), and divide by 0
00510   if (u < -1.0)       {u = -1.0;}
00511   if (u > kAlmostOne) {u = kAlmostOne;}
00512   Angle = acos(u);
00513 
00514   Const1 = -1.0 / sqrt(1.0 - u*u);
00515   Const2 = m_Kf * LambdaKf * (Angle-RefAngle) * Const1;
00516 
00517   Const3 = -a/(2.0*c*c) + 1.0/(a+a) + (b*b)/(2.0*a*c*c);
00518   Const4 = -c/(2.0*a*a) + 1.0/(c+c) + (b*b)/(2.0*c*a*a);
00519   Const5 = -b/(a*c);
00520 
00521   if (WhichGroup == 0) {
00522     Vec1 = (A-C) * (Const5/b);
00523     Vec2 = (A-B) * (Const3/c);
00524   }
00525   else if (WhichGroup == 1) {
00526     Vec1 = (B-A) * (Const3/c);
00527     Vec2 = (B-C) * (Const4/a);
00528   }
00529   else {
00530     Vec1 = (C-A) * (Const5/b);
00531     Vec2 = (C-B) * (Const4/a);
00532   }
00533   return( (Vec1+Vec2)*Const2);
00534 }
00535 
00536 
00537 ADiheRestraint::ADiheRestraint() {
00538 //----------------------------------------------------------------------------
00539 // each ADiheRestraint restrains the dihedral angle between 4 groups of atoms
00540 //----------------------------------------------------------------------------
00541   m_NumGroups = 4;
00542   m_pGroups = new AGroup[m_NumGroups];
00543   m_pCOMs = new AVector[m_NumGroups];
00544 }
00545 
00546 
00547 void ADiheRestraint::PrintInfo() {
00548 //--------------------------------------------------------------------
00549 // print the dihedral angle for this dihedral restraint
00550 //--------------------------------------------------------------------
00551   double  Dihedral;
00552   char    Str1[20], Str2[20], Str3[20];
00553 
00554   Dihedral = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]) * (180/kPi);
00555   sprintf(Str1, "%8.3f", Dihedral);
00556   Dihedral = GetDiheTarget1() * (180/kPi);
00557   sprintf(Str2, "%8.3f", Dihedral);
00558   Dihedral = GetDiheTarget2() * (180/kPi);
00559   sprintf(Str3, "%8.3f", Dihedral);
00560 
00561 #if defined(_VERBOSE_PMF)
00562   iout << "Dihedral = ";
00563   iout << Str1;
00564   iout << " degrees";
00565   iout << "  Target = ";
00566   iout << Str2;
00567   iout << " degrees";
00568   if (TwoTargets()) {
00569     iout << " to ";
00570     iout << Str3;
00571     iout << " degrees";
00572   }
00573   iout << std::endl << endi;
00574 #else
00575   iout << Str1;
00576   iout << "  ";
00577   iout << Str2;
00578   if (TwoTargets()) {
00579     iout << ", ";
00580     iout << Str3;
00581   }
00582   iout << " | ";
00583 #endif
00584 }
00585 
00586 
00587 double ADiheRestraint::GetE(double RefAngle, double Const) {
00588 //---------------------------------------------------------------------------
00589 // calculate and return the Energy for this angle restraint.
00590 //
00591 //    E = (E0/2) * (1 - cos(Chi - ChiRef))
00592 //
00593 // where Chi is the dihedral angle between 4 centers-of-mass of restrained atoms,
00594 // m_pCOMs[0] -- m_pCOMs[1] -- m_pCOMs[2] -- m_pCOMs[3].
00595 // ChiRef is the reference angle.
00596 //
00597 // Note:  COM's are calculated before this routine is called.
00598 //---------------------------------------------------------------------------
00599   double  Angle;
00600 
00601   Angle = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
00602   return( (Const/2.0) * (1.0 - cos(Angle-RefAngle)) );
00603 }
00604 
00605 
00606 AVector ADiheRestraint::GetGrad(int WhichGroup,
00607                                 double RefAngle, double Const) {
00608 //---------------------------------------------------------------------------
00609 // calculate and return the gradient for this dihedral angle restraint.
00610 //
00611 //    E = (E0/2) * (1 - cos(Chi - ChiRef))
00612 //
00613 // return:  grad(E)
00614 //
00615 // Notes: COM's are calculated before this routine is called.
00616 //---------------------------------------------------------------------------
00617   AVector A, B, C, D;
00618 
00619   ASSERT((WhichGroup==0)||(WhichGroup==1)||(WhichGroup==2)||(WhichGroup==3));
00620 
00621   if ((WhichGroup==0) || (WhichGroup==1)) {
00622     A = m_pCOMs[0];
00623     B = m_pCOMs[1];
00624     C = m_pCOMs[2];
00625     D = m_pCOMs[3];
00626   }
00627   // re-state the problem so the gradient is solved for either atoms 0 or 1
00628   else {
00629     A = m_pCOMs[3];
00630     B = m_pCOMs[2];
00631     C = m_pCOMs[1];
00632     D = m_pCOMs[0];
00633     if (WhichGroup==3) {WhichGroup=0;}
00634     if (WhichGroup==2) {WhichGroup=1;}
00635   }
00636 
00637   AVector CD(D - C);
00638   AVector CB(B - C);
00639   AVector BC(C - B);
00640   AVector BA(A - B);
00641   AVector AC(C - A);
00642   AVector CDxCB, BCxBA;
00643   AVector Vec;
00644   double  phi;
00645   double  top, bot, u;
00646 
00647   CDxCB = CD.cross(CB);
00648   BCxBA = BC.cross(BA);
00649 
00650   top = CDxCB.dot(BCxBA);
00651   bot = CDxCB.Dist() * BCxBA.Dist();
00652 
00653   u = top/bot;
00654   // protect against acos(<-1.0), acos(>1.0), sqrt(<0), and divide by 0
00655   if (u < kAlmostMinusOne) {u = kAlmostMinusOne;}
00656   if (u > kAlmostOne)      {u = kAlmostOne;}
00657 
00658   // get dihedral using atan
00659   phi = GetDihe(A,B,C,D);
00660 
00661   AVector dP1, dP2, dP3;
00662   AVector dP4, dP5, dP6;
00663   ASSERT((WhichGroup==0) || (WhichGroup==1));
00664   if (WhichGroup==0) {
00665     dP1.Set( 0,      0,      0    );
00666     dP2.Set( 0,      0,      0    );
00667     dP3.Set( 0,      0,      0    );
00668     dP4.Set( 0,     -BC[2],  BC[1]);
00669     dP5.Set( BC[2],  0,     -BC[0]);
00670     dP6.Set(-BC[1],  BC[0],  0    );
00671   }
00672   else {
00673     dP1.Set( 0,     -CD[2],  CD[1]);
00674     dP2.Set( CD[2],  0,     -CD[0]);
00675     dP3.Set(-CD[1],  CD[0],  0    );
00676     dP4.Set( 0,      AC[2], -AC[1]);
00677     dP5.Set(-AC[2],  0,      AC[0]);
00678     dP6.Set( AC[1], -AC[0],  0    );
00679   }
00680 
00681   Vec = gradU(CDxCB, BCxBA, dP1, dP2, dP3, dP4, dP5, dP6);
00682   Vec *= (Const/2.0) * sin(phi-RefAngle) * (-1.0/sqrt(1.0 - u*u));
00683 
00684   // flip gradient for negative angles
00685   if (phi < 0) {
00686     Vec *= -1.0;
00687   }
00688   
00689   return(Vec);
00690 }
00691 
00692 
00693 AVector ADiheRestraint::gradU(AVector& P1P2P3, AVector& P4P5P6,
00694                               AVector& dP1,    AVector& dP2,    AVector& dP3,
00695                               AVector& dP4,    AVector& dP5,    AVector& dP6) {
00696 //----------------------------------------------------------------
00697 // calculate the gradient for ((P1P2P3.dot.P4P5P6)/(mag(P1P2P3)*mag(P4P5P6)))
00698 // P1P2P3 = (P1)i + (P2)j + (P3)k
00699 // P4P5P6 = (P4)i + (P5)j + (P6)k
00700 // dP1 = (d(P1)/dx)i + (d(P1)/dy)j +(d(P1)/dz)k
00701 // dP2 = (d(P2)/dx)i + (d(P2)/dy)j +(d(P2)/dz)k
00702 // dP3 = (d(P3)/dx)i + (d(P3)/dy)j +(d(P3)/dz)k
00703 // dP4 = (d(P4)/dx)i + (d(P4)/dy)j +(d(P4)/dz)k
00704 // dP5 = (d(P5)/dx)i + (d(P5)/dy)j +(d(P5)/dz)k
00705 // dP6 = (d(P6)/dx)i + (d(P6)/dy)j +(d(P6)/dz)k
00706 //----------------------------------------------------------------
00707   double  Mag123, Mag456, Dot;
00708   double  Const1, Const2, Const3;
00709   double  P1, P2, P3, P4, P5, P6;
00710   AVector RetVec;
00711 
00712   P1 = P1P2P3[0];  P2 = P1P2P3[1];  P3 = P1P2P3[2];
00713   P4 = P4P5P6[0];  P5 = P4P5P6[1];  P6 = P4P5P6[2];
00714 
00715   Mag123 = P1P2P3.Dist();
00716   Mag456 = P4P5P6.Dist();
00717   Dot    = P1P2P3.dot(P4P5P6);
00718 
00719   Const1 =         1.0 / (Mag123*Mag456);
00720   Const2 = -Dot * (1.0 / (Mag123*Mag456*Mag456*Mag456));
00721   Const3 = -Dot * (1.0 / (Mag456*Mag123*Mag123*Mag123));
00722 
00723   RetVec = (dP4*P1 + dP1*P4 + dP5*P2 + dP2*P5 + dP6*P3 + dP3*P6) * Const1 +
00724            (dP4*P4 + dP5*P5 + dP6*P6)                            * Const2 +
00725            (dP1*P1 + dP2*P2 + dP3*P3)                            * Const3;
00726 
00727   return(RetVec);
00728 }
00729 
00730 
00731 double AFixedPosRestraint::GetEnergy() {
00732 //--------------------------------------------------------------------
00733 // return the Energy for this fixed position restraint.
00734 //--------------------------------------------------------------------
00735   return(GetE(m_RefPos));
00736 }
00737 
00738 
00739 AVector AFixedPosRestraint::GetGradient(int WhichGroup) {
00740 //-------------------------------------------------------------------------
00741 // return the Gradient for this fixed position restraint.
00742 //-------------------------------------------------------------------------
00743   return(GetGrad(WhichGroup, m_RefPos));
00744 }
00745 
00746 
00747 double ABoundPosRestraint::GetEnergy() {
00748 //--------------------------------------------------------------------
00749 // calculate and return the Energy for this bound position restraint.
00750 //
00751 // Note:  This is an exception because the form for the E term is
00752 //        different from the other postion restraints.
00753 //--------------------------------------------------------------------
00754   double  E, Dist, Diff;
00755 
00756   E = 0.0;
00757   Dist = m_pCOMs[0].Dist(m_RefPos);
00758   if (((m_Bound==kUpper) && (Dist>m_RefDist)) ||
00759       ((m_Bound==kLower) && (Dist<m_RefDist))) {
00760     Diff = Dist - m_RefDist;
00761     E = (m_Kf/2.0) * (Diff*Diff);
00762   }
00763   return(E);
00764 }
00765 
00766 
00767 AVector ABoundPosRestraint::GetGradient(int /* WhichGroup */) {
00768 //---------------------------------------------------------------------------
00769 // calculate and return the gradient for this bound position restraint.
00770 //
00771 // Note:  This is an exception because the form for the E term is
00772 //        different from the other postion restraints.
00773 //---------------------------------------------------------------------------
00774   double  Dist;
00775   AVector Vec;   // Vec is initialized to (0,0,0)
00776 
00777   // WhichGroup = 0;  // don't care -- there's only 1 atom restrained
00778   Dist = m_pCOMs[0].Dist(m_RefPos);
00779   if (((m_Bound==kUpper) && (Dist>m_RefDist)) ||
00780       ((m_Bound==kLower) && (Dist<m_RefDist))) {
00781     Vec.Set(m_pCOMs[0][0] - m_RefPos[0],
00782             m_pCOMs[0][1] - m_RefPos[1],
00783             m_pCOMs[0][2] - m_RefPos[2]);
00784     Vec *= m_Kf * (Dist - m_RefDist) / Dist;
00785   }
00786   return(Vec);
00787 }
00788 
00789 
00790 double AForcingPosRestraint::GetEnergy() {
00791 //--------------------------------------------------------------------
00792 // return the Energy for this forcing position restraint.
00793 //
00794 // rref = lambda*r1 + (1-lambda)*r0.
00795 // where r0 is the starting position and r1 is the final position
00796 //--------------------------------------------------------------------
00797   AVector  RefPos;
00798 
00799   RefPos = m_StopPos*m_LambdaRef + m_StartPos*(1.0-m_LambdaRef);
00800   return(GetE(RefPos, m_LambdaKf));
00801 }
00802 
00803 
00804 AVector AForcingPosRestraint::GetGradient(int WhichGroup) {
00805 //---------------------------------------------------------------------------
00806 // return the gradient for this forcing position restraint.
00807 //
00808 // rref = lambda*r1 + (1-lambda)*r0.
00809 // where r0 is the starting position and r1 is the final position
00810 //---------------------------------------------------------------------------
00811   AVector  RefPos;
00812 
00813   RefPos = m_StopPos*m_LambdaRef + m_StartPos*(1.0-m_LambdaRef);
00814   return(GetGrad(WhichGroup, RefPos, m_LambdaKf));
00815 }
00816 
00817 
00818 double AForcingPosRestraint::Get_dU_dLambda() {
00819 //---------------------------------------------------------------------------
00820 // return dU/dLambda for this forcing position restraint
00821 //---------------------------------------------------------------------------
00822   AVector  RefPos;
00823   double   T1, T2, T3;
00824 
00825   RefPos = m_StopPos*m_LambdaRef + m_StartPos*(1.0-m_LambdaRef);
00826   T1 = (m_pCOMs[0][0] - RefPos[0]) * (m_StartPos[0] - m_StopPos[0]);
00827   T2 = (m_pCOMs[0][1] - RefPos[1]) * (m_StartPos[1] - m_StopPos[1]);
00828   T3 = (m_pCOMs[0][2] - RefPos[2]) * (m_StartPos[2] - m_StopPos[2]);
00829   return( m_Kf * m_LambdaKf * (T1+T2+T3) );
00830 }
00831 
00832 
00833 double AFixedDistRestraint::GetEnergy() {
00834 //---------------------------------------------------------------------------
00835 // return the Energy for this fixed distance restraint.
00836 //---------------------------------------------------------------------------
00837   return(GetE(m_RefDist));
00838 }
00839 
00840 
00841 AVector AFixedDistRestraint::GetGradient(int WhichGroup) {
00842 //---------------------------------------------------------------------------
00843 // return the gradient for this fixed distance restraint.
00844 //---------------------------------------------------------------------------
00845   return(GetGrad(WhichGroup, m_RefDist));
00846 }
00847 
00848 
00849 double ABoundDistRestraint::GetEnergy() {
00850 //---------------------------------------------------------------------------
00851 // return the Energy for this bound distance restraint.
00852 //---------------------------------------------------------------------------
00853   double Dist, E;
00854 
00855   E = 0.0;
00856   Dist = m_pCOMs[0].Dist(m_pCOMs[1]);
00857   if (((m_Bound==kUpper) && (Dist>m_RefDist)) ||
00858       ((m_Bound==kLower) && (Dist<m_RefDist))) {
00859     E = GetE(m_RefDist);
00860   }
00861   return(E);
00862 }
00863 
00864 
00865 AVector ABoundDistRestraint::GetGradient(int WhichGroup) {
00866 //---------------------------------------------------------------------------
00867 // return the gradient for this bound distance restraint.
00868 //---------------------------------------------------------------------------
00869   double  Dist;
00870   AVector Vec;
00871 
00872   Dist = m_pCOMs[0].Dist(m_pCOMs[1]);
00873   if (((m_Bound==kUpper) && (Dist>m_RefDist)) ||
00874       ((m_Bound==kLower) && (Dist<m_RefDist))) {
00875     Vec = GetGrad(WhichGroup, m_RefDist);
00876   }
00877   return(Vec);
00878 }
00879 
00880 
00881 double AForcingDistRestraint::GetEnergy() {
00882 //---------------------------------------------------------------------------
00883 // return the Energy for this forcing distance restraint.
00884 //---------------------------------------------------------------------------
00885   double  RefDist;
00886 
00887   RefDist = m_StopDist*m_LambdaRef + m_StartDist*(1.0-m_LambdaRef);
00888   return(GetE(RefDist, m_LambdaKf));
00889 }
00890 
00891 
00892 AVector AForcingDistRestraint::GetGradient(int WhichGroup) {
00893 //---------------------------------------------------------------------------
00894 // return the gradient for this forcing distance restraint.
00895 //---------------------------------------------------------------------------
00896   double  RefDist;
00897   
00898   RefDist = m_StopDist*m_LambdaRef + m_StartDist*(1.0-m_LambdaRef);
00899   return(GetGrad(WhichGroup, RefDist, m_LambdaKf));
00900 }
00901 
00902 
00903 double AForcingDistRestraint::Get_dU_dLambda() {
00904 //---------------------------------------------------------------------------
00905 // return dU/dLambda for this forcing distance restraint
00906 //---------------------------------------------------------------------------
00907   double  Dist;
00908   double  RefDist;
00909   
00910   Dist = m_pCOMs[0].Dist(m_pCOMs[1]);
00911   RefDist = m_StopDist*m_LambdaRef + m_StartDist*(1.0-m_LambdaRef);
00912   return( m_Kf * m_LambdaKf * (Dist-RefDist)*(m_StartDist-m_StopDist) );
00913 }
00914 
00915 
00916 double AFixedAngleRestraint::GetEnergy() {
00917 //---------------------------------------------------------------------------
00918 // return the Energy for this fixed angle restraint.
00919 //---------------------------------------------------------------------------
00920   return(GetE(m_RefAngle));
00921 }
00922 
00923 
00924 AVector AFixedAngleRestraint::GetGradient(int WhichGroup) {
00925 //---------------------------------------------------------------------------
00926 // return the gradient for this fixed angle restraint.
00927 //---------------------------------------------------------------------------
00928   return(GetGrad(WhichGroup, m_RefAngle));
00929 }
00930 
00931 
00932 double ABoundAngleRestraint::GetEnergy() {
00933 //---------------------------------------------------------------------------
00934 // return the Energy for this bound angle restraint.
00935 //---------------------------------------------------------------------------
00936   double  E, Angle;
00937 
00938   E = 0.0;
00939   Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]);
00940   if (((m_Bound==kUpper) && (Angle>m_RefAngle)) ||
00941       ((m_Bound==kLower) && (Angle<m_RefAngle))) {
00942     E = GetE(m_RefAngle);
00943   }
00944   return(E);
00945 }
00946 
00947 
00948 AVector ABoundAngleRestraint::GetGradient(int WhichGroup) {
00949 //---------------------------------------------------------------------------
00950 // return the gradient for this bound angle restraint
00951 //---------------------------------------------------------------------------
00952   double  Angle;
00953   AVector Vec;
00954 
00955   Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]);
00956   if (((m_Bound==kUpper) && (Angle>m_RefAngle)) ||
00957       ((m_Bound==kLower) && (Angle<m_RefAngle))) {
00958     Vec = GetGrad(WhichGroup, m_RefAngle);
00959   }
00960   return(Vec);
00961 }
00962 
00963 
00964 double AForcingAngleRestraint::GetEnergy() {
00965 //---------------------------------------------------------------------------
00966 // return the Energy for this forcing angle restraint.
00967 //---------------------------------------------------------------------------
00968   double  RefAngle;
00969 
00970   RefAngle = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
00971   return(GetE(RefAngle, m_LambdaKf));
00972 }
00973 
00974 
00975 AVector AForcingAngleRestraint::GetGradient(int WhichGroup) {
00976 //---------------------------------------------------------------------------
00977 // return the gradient for this forcing angle restraint.
00978 //---------------------------------------------------------------------------
00979   double  RefAngle;
00980 
00981   RefAngle = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
00982   return(GetGrad(WhichGroup, RefAngle, m_LambdaKf));
00983 }
00984 
00985 
00986 double AForcingAngleRestraint::Get_dU_dLambda() {
00987 //---------------------------------------------------------------------------
00988 // return dU/dLambda for this forcing angle restraint
00989 //---------------------------------------------------------------------------
00990   double  Angle;
00991   double  RefAngle;
00992 
00993   Angle = GetAngle(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2]);
00994   RefAngle = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
00995   return( m_Kf * m_LambdaKf * (Angle-RefAngle)*(m_StartAngle-m_StopAngle) );
00996 }
00997 
00998 
00999 double AFixedDiheRestraint::GetEnergy() {
01000 //---------------------------------------------------------------------------
01001 // return the Energy for this fixed dihedral angle restraint.
01002 //---------------------------------------------------------------------------
01003   return(GetE(m_RefAngle, m_Kf));
01004 }
01005 
01006 
01007 AVector AFixedDiheRestraint::GetGradient(int WhichGroup) {
01008 //---------------------------------------------------------------------------
01009 // return the gradient for this fixed dihedral angle restraint.
01010 //---------------------------------------------------------------------------
01011   return(GetGrad(WhichGroup, m_RefAngle, m_Kf));
01012 }
01013 
01014 
01015 double ABoundDiheRestraint::GetEnergy() {
01016 //---------------------------------------------------------------------------
01017 // return the Energy for this bound dihedral angle restraint.
01018 //---------------------------------------------------------------------------
01019   double  E, Dihe, Const;
01020 
01021   Const = m_Kf / (1.0 - cos(m_IntervalAngle));
01022   Dihe = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
01023   // dihedral angle is between LowerAngle and UpperAngle
01024   if ( (Dihe>m_LowerAngle) && (Dihe<m_UpperAngle) ) {
01025     E = 0.0;
01026   }
01027   // dihedral angle is between LowerAngle and LowerAngle-IntervalAngle
01028   else if ( (Dihe<m_LowerAngle) && (Dihe>(m_LowerAngle-m_IntervalAngle)) ) {
01029     E = GetE(m_LowerAngle, Const);
01030   }
01031   // dihedral angle is between UpperAngle and UpperAngle+IntervalAngle
01032   else if ( (Dihe>m_UpperAngle) && (Dihe<(m_UpperAngle+m_IntervalAngle)) ) {
01033     E = GetE(m_UpperAngle, Const);
01034   }
01035   // dihedral angle is more than UpperAngle or less than LowerAngle
01036   else {
01037     E = Const;
01038   }
01039   return(E);
01040 }
01041 
01042 
01043 AVector ABoundDiheRestraint::GetGradient(int WhichGroup) {
01044 //---------------------------------------------------------------------------
01045 // return the gradient for this bound dihedral angle restraint.
01046 //---------------------------------------------------------------------------
01047   AVector Vec;
01048   double  Dihe, Const;
01049 
01050   Const = m_Kf / (1.0 - cos(m_IntervalAngle));
01051   Dihe = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
01052   // dihedral angle is between LowerAngle and LowerAngle-IntervalAngle
01053   if ( (Dihe<m_LowerAngle) && (Dihe>(m_LowerAngle-m_IntervalAngle)) ) {
01054     Vec = GetGrad(WhichGroup, m_LowerAngle, Const);
01055   }
01056   // dihedral angle is between UpperAngle and UpperAngle+IntervalAngle
01057   else if ( (Dihe>m_UpperAngle) && (Dihe<(m_UpperAngle+m_IntervalAngle)) ) {
01058     Vec = GetGrad(WhichGroup, m_UpperAngle, Const);
01059   }
01060   return(Vec);
01061 }
01062 
01063 
01064 double AForcingDiheRestraint::GetEnergy() {
01065 //---------------------------------------------------------------------------
01066 // return the Energy for this forcing dihedral angle restraint.
01067 //---------------------------------------------------------------------------
01068   double  RefDihe;
01069 
01070   RefDihe = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
01071   return(GetE(RefDihe, m_Kf*m_LambdaKf));
01072 }
01073 
01074 
01075 AVector AForcingDiheRestraint::GetGradient(int WhichGroup) {
01076 //---------------------------------------------------------------------------
01077 // return the gradient for this forcing dihedral angle restraint.
01078 //---------------------------------------------------------------------------
01079   double  RefDihe;
01080   
01081   RefDihe = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
01082   return(GetGrad(WhichGroup, RefDihe, m_Kf*m_LambdaKf));
01083 }
01084 
01085 
01086 double AForcingDiheRestraint::Get_dU_dLambda() {
01087 //---------------------------------------------------------------------------
01088 // return dU/dLambda for this forcing dihedral angle restraint
01089 //---------------------------------------------------------------------------
01090   double  Dihe;
01091   double  RefDihe;
01092   
01093   Dihe = GetDihe(m_pCOMs[0], m_pCOMs[1], m_pCOMs[2], m_pCOMs[3]);
01094   RefDihe = m_StopAngle*m_LambdaRef + m_StartAngle*(1.0-m_LambdaRef);
01095   return((m_Kf/2)*m_LambdaKf * sin(Dihe-RefDihe) * (m_StartAngle-m_StopAngle));
01096 }
01097 

Generated on Sat Sep 22 01:17:13 2018 for NAMD by  doxygen 1.4.7