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

Generated on Mon Nov 20 01:17:12 2017 for NAMD by  doxygen 1.4.7