FreeEnergyLambda.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 "InfoStream.h"
00011 #include "FreeEnergyEnums.h"
00012 #include "FreeEnergyAssert.h"
00013 #include "FreeEnergyLambda.h"
00014 #include "Vector.h"
00015 #include "FreeEnergyVector.h"
00016 
00017 // initialize static member variables
00018 int  ALambdaControl::m_CurrStep = 0;
00019 
00020 ALambdaControl::ALambdaControl() {
00021 //------------------------------------------------------------------------
00022 // initialize member variables
00023 //------------------------------------------------------------------------
00024   // initialize these for the first block.
00025   m_StartStep = 0;
00026   m_LambdaKf  = 1.0;         // applied to Kf
00027   m_LambdaRef = 0.0;         // applied to reference (pos, dist, angle, dihe)
00028   m_Sum_dU_dLambda = 0.0;    // for free energy measurement
00029   m_Num_dU_dLambda = 0;      // "
00030   m_MCTI_Integration = 0.0;  // "
00031 
00032   // set the rest to illegal values.
00033   m_Task =           kUnknownTask;
00034   m_NumSteps =      -1;
00035   m_NumEquilSteps = -1;
00036   m_NumAccumSteps = -1;
00037   m_NumPrintSteps = -1;
00038   m_NumRepeats =    -1;
00039   m_StopStep =      -1;
00040 }
00041 
00042 
00043 void ALambdaControl::Init(ALambdaControl& PriorBlock) {
00044 //------------------------------------------------------------------------
00045 // initialize this object using the settings for the prior object.
00046 //------------------------------------------------------------------------
00047   m_Task =          PriorBlock.m_Task;
00048   m_NumSteps =      PriorBlock.m_NumSteps;
00049   m_NumEquilSteps = PriorBlock.m_NumEquilSteps;
00050   m_NumAccumSteps = PriorBlock.m_NumAccumSteps;
00051   m_NumPrintSteps = PriorBlock.m_NumPrintSteps;
00052   m_NumRepeats =    PriorBlock.m_NumRepeats;
00053   switch (PriorBlock.m_Task) {
00054 
00055     case kUp:       m_LambdaKf=1.0;  m_LambdaRef=1.0; break;
00056     case kStepUp:   m_LambdaKf=1.0;  m_LambdaRef=1.0; break;
00057     case kDown:     m_LambdaKf=1.0;  m_LambdaRef=0.0; break;
00058     case kStepDown: m_LambdaKf=1.0;  m_LambdaRef=0.0; break;
00059     case kStop:     m_LambdaKf=1.0;  m_LambdaRef=PriorBlock.m_LambdaRef; break;
00060 
00061     case kGrow:     m_LambdaKf=1.0;  m_LambdaRef=PriorBlock.m_LambdaRef; break;
00062     case kStepGrow: m_LambdaKf=1.0;  m_LambdaRef=PriorBlock.m_LambdaRef; break;
00063     case kFade:     m_LambdaKf=0.0;  m_LambdaRef=PriorBlock.m_LambdaRef; break;
00064     case kStepFade: m_LambdaKf=0.0;  m_LambdaRef=PriorBlock.m_LambdaRef; break;
00065     case kNoGrow:   m_LambdaKf =  PriorBlock.m_LambdaKf;
00066                     m_LambdaRef = PriorBlock.m_LambdaRef;  break;
00067     default:        ASSERT(kFalse); break;  //should never get here
00068 
00069   }
00070   m_StartStep = PriorBlock.GetLastStep();
00071 }
00072 
00073 
00074 ALambdaControl& ALambdaControl::operator= (ALambdaControl& PmfBlock) {
00075 //------------------------------------------------------------------------
00076 // copy everything from PmfBlock to this block.
00077 //------------------------------------------------------------------------
00078   m_NumSteps =      PmfBlock.m_NumSteps;
00079   m_NumEquilSteps = PmfBlock.m_NumEquilSteps;
00080   m_NumAccumSteps = PmfBlock.m_NumAccumSteps;
00081   m_NumRepeats =    PmfBlock.m_NumRepeats;
00082   m_NumPrintSteps = PmfBlock.m_NumPrintSteps;
00083   m_StartStep =     PmfBlock.m_StartStep;
00084   m_StopStep =      PmfBlock.m_StopStep;
00085   m_LambdaKf =      PmfBlock.m_LambdaKf;
00086   m_LambdaRef =     PmfBlock.m_LambdaRef;
00087   m_Task =          PmfBlock.m_Task;
00088   return(*this);
00089 }
00090 
00091 
00092 void ALambdaControl::Accumulate(double dU_dLambda) {
00093 //------------------------------------------------------------------------
00094 // if lambda is constant, sum dU/dLambda
00095 // if lambda is changing, sum dU/dLambda * dLambda
00096 //------------------------------------------------------------------------
00097 
00098   m_Num_dU_dLambda++;
00099   switch (m_Task) {
00100     // lambda is constant
00101     case kStop:      case kNoGrow:    case kStepUp:
00102     case kStepDown:  case kStepGrow:  case kStepFade:
00103       m_Sum_dU_dLambda += dU_dLambda;
00104       break;
00105     // lambda is increasing
00106     case kUp:
00107     case kGrow:
00108       m_Sum_dU_dLambda += dU_dLambda / (double)m_NumSteps;
00109       break;
00110     // lambda is decreasing
00111     case kDown:
00112     case kFade:
00113       m_Sum_dU_dLambda -= dU_dLambda / (double)m_NumSteps;
00114       break;
00115     // should never get here
00116     default:
00117       ASSERT(kFalse);
00118       break;
00119   }
00120 }
00121 
00122 
00123 void ALambdaControl::Integrate_MCTI() {
00124 //------------------------------------------------------------------------
00125 // integrate MCTI:  <dU/dLambda> * dLambda
00126 //------------------------------------------------------------------------
00127   ASSERT(m_Task==kStepUp   || m_Task==kStepDown || 
00128          m_Task==kStepGrow || m_Task==kStepFade);
00129 
00130   switch (m_Task) {
00131     // lambda is increasing
00132     case kStepUp:
00133     case kStepGrow:
00134       m_MCTI_Integration += GetAccumulation() / (double)m_NumRepeats;
00135       break;
00136     // lambda is decreasing
00137     case kStepDown:
00138     case kStepFade:
00139       m_MCTI_Integration -= GetAccumulation() / (double)m_NumRepeats;
00140       break;
00141     // should never get here
00142     default:
00143       ASSERT(kFalse);
00144       break;
00145   }
00146 }
00147 
00148 
00149 double ALambdaControl::GetIntegration() {
00150 //------------------------------------------------------------------------
00151 // get MCTI integral.  integral(dU/dLambda> * dLambda
00152 //------------------------------------------------------------------------
00153   ASSERT(m_Task==kStepUp   || m_Task==kStepDown || 
00154          m_Task==kStepGrow || m_Task==kStepFade);
00155 
00156   return(m_MCTI_Integration);
00157 }
00158 
00159 
00160 double ALambdaControl::GetAccumulation() {
00161 //------------------------------------------------------------------------
00162 // if lambda is constant, return (sum/N)
00163 // if lambda is changing, return (sum)
00164 //------------------------------------------------------------------------
00165   switch (m_Task) {
00166     // lambda is constant
00167     case kStop:    case kNoGrow:
00168     case kStepUp:  case kStepDown:  case kStepGrow:  case kStepFade:
00169       ASSERT(m_Num_dU_dLambda != 0);
00170       return(m_Sum_dU_dLambda / (double)m_Num_dU_dLambda);
00171     // lambda is changing
00172     case kUp:  case kDown:  case kGrow:  case kFade:
00173       return(m_Sum_dU_dLambda);
00174     // should never get here
00175     default:
00176       ASSERT(kFalse);
00177       return(0);
00178   }
00179 }
00180 
00181 
00182 int ALambdaControl::GetNumAccumStepsSoFar() {
00183 //------------------------------------------------------------------------
00184 // return the total number of steps dU/dLambda has been accumulated
00185 // this is only called during MCTI
00186 //------------------------------------------------------------------------
00187   ASSERT(m_Task==kStepUp   || m_Task==kStepDown || 
00188          m_Task==kStepGrow || m_Task==kStepFade);
00189 
00190   int Count, Maybe;
00191 
00192   Count = (m_CurrStep-m_StartStep) / (m_NumAccumSteps+m_NumEquilSteps);
00193   Count *= m_NumAccumSteps;
00194   Maybe = (m_CurrStep-m_StartStep) % (m_NumAccumSteps+m_NumEquilSteps);
00195   Maybe -= m_NumEquilSteps;
00196   if (Maybe > 0) {
00197     Count += Maybe;
00198   }
00199   return(Count);
00200 }
00201 
00202 
00203 int ALambdaControl::GetNumSteps() {
00204 //------------------------------------------------------------------------
00205 // get the number of steps needed for this pmf or mcti block
00206 //------------------------------------------------------------------------
00207   // make sure m_StopStep is calculated
00208   GetLastStep();
00209   
00210   return( (m_StopStep - m_StartStep) + 1 );
00211 }
00212 
00213 
00214 Bool_t ALambdaControl::IsLastStep() {
00215 //------------------------------------------------------------------------
00216 // return true if we're on the last step of this pmf or mcti block
00217 //------------------------------------------------------------------------
00218         if (m_CurrStep == GetLastStep()) {
00219     return(kTrue);
00220   }
00221   else {
00222     return(kFalse);
00223   }
00224 }
00225 
00226 
00227 int ALambdaControl::GetLastStep() {
00228 //------------------------------------------------------------------------
00229 // get the last step of this task
00230 //------------------------------------------------------------------------
00231   // if it's already calculated, just return it
00232   if (m_StopStep > 0) {
00233     return(m_StopStep);
00234   }
00235   // otherwise calculate it
00236   switch (m_Task) {
00237     case kStepUp:
00238     case kStepDown:
00239     case kStepGrow:
00240     case kStepFade:
00241       m_StopStep = m_StartStep +
00242                   (m_NumAccumSteps+m_NumEquilSteps) * m_NumRepeats;
00243       break;
00244     default:
00245       m_StopStep = m_StartStep + m_NumSteps;
00246       break;
00247   }
00248   // and return it
00249   return(m_StopStep);
00250 }
00251 
00252 
00253 void ALambdaControl::GetPaddedTaskStr(char* Str) {
00254 //------------------------------------------------------------------------
00255 // get a string that describes this task
00256 //------------------------------------------------------------------------
00257   switch (m_Task) {
00258     case kUp:        strcpy(Str, "      Up");      break;
00259     case kDown:      strcpy(Str, "    Down");      break;
00260     case kStop:      strcpy(Str, "    Stop");      break;
00261     case kGrow:      strcpy(Str, "    Grow");      break;
00262     case kFade:      strcpy(Str, "    Fade");      break;
00263     case kNoGrow:    strcpy(Str, "  NoGrow");      break;
00264     case kStepUp:    strcpy(Str, "  StepUp");      break;
00265     case kStepDown:  strcpy(Str, "StepDown");      break;
00266     case kStepGrow:  strcpy(Str, "StepGrow");      break;
00267     case kStepFade:  strcpy(Str, "StepFade");      break;
00268     default:         strcpy(Str, "Bug Alert!!!");  break;
00269   }
00270 }
00271 
00272 
00273 void ALambdaControl::GetTaskStr(char* Str) {
00274 //------------------------------------------------------------------------
00275 // get a string that describes this task
00276 //------------------------------------------------------------------------
00277   switch (m_Task) {
00278     case kUp:        strcpy(Str, "Up");            break;
00279     case kDown:      strcpy(Str, "Down");          break;
00280     case kStop:      strcpy(Str, "Stop");          break;
00281     case kGrow:      strcpy(Str, "Grow");          break;
00282     case kFade:      strcpy(Str, "Fade");          break;
00283     case kNoGrow:    strcpy(Str, "NoGrow");        break;
00284     case kStepUp:    strcpy(Str, "StepUp");        break;
00285     case kStepDown:  strcpy(Str, "StepDown");      break;
00286     case kStepGrow:  strcpy(Str, "StepGrow");      break;
00287     case kStepFade:  strcpy(Str, "StepFade");      break;
00288     default:         strcpy(Str, "Bug Alert!!!");  break;
00289   }
00290 }
00291 
00292 
00293 Bool_t ALambdaControl::IsTimeToClearAccumulator() {
00294 //------------------------------------------------------------------------
00295 // ASSUMING that this object is currently active, decide if it's time
00296 // to start accumulating dU/dLambda from zero.
00297 // (clear the accumulator on the 1st step)
00298 //------------------------------------------------------------------------
00299   Bool_t  RetVal=kFalse;
00300 
00301   switch (m_Task) {
00302     // for pmf blocks, clear accumulator on first step
00303     case kUp:      case kDown:    case kStop:
00304     case kGrow:    case kFade:    case kNoGrow:
00305       if (m_CurrStep == m_StartStep) {
00306         RetVal = kTrue;
00307       }
00308       break;
00309     // for mcti blocks, clear accumulator after each equilibration
00310     case kStepUp:  case kStepDown:  case kStepGrow:  case kStepFade:
00311       if ( (m_CurrStep-(m_StartStep+m_NumEquilSteps)) % 
00312              (m_NumEquilSteps+m_NumAccumSteps) == 1) {
00313         RetVal = kTrue;
00314       }
00315       break;
00316     // should never get here
00317     default:
00318       ASSERT(kFalse);
00319       break;
00320   }
00321   return(RetVal);
00322 }
00323 
00324 
00325 Bool_t ALambdaControl::IsFirstStep() {
00326 //------------------------------------------------------------------------
00327 // ASSUMING that this object is currently active, decide if it's the
00328 // first step of the control object.
00329 //------------------------------------------------------------------------
00330   Bool_t  RetVal=kFalse;
00331 
00332   if (m_CurrStep == (m_StartStep+1)) {
00333     RetVal = kTrue;
00334   }
00335   return(RetVal);
00336 }
00337 
00338 
00339 Bool_t ALambdaControl::IsTimeToPrint() {
00340 //------------------------------------------------------------------------
00341 // ASSUMING that this object is currently active, decide if it's time
00342 // to print out restraint information.
00343 //------------------------------------------------------------------------
00344   Bool_t  RetVal=kFalse;
00345 
00346   // if printing is required
00347   if (m_NumPrintSteps > 0) {
00348     // if number-of-steps from StartStep is an even multiple of NumPrintSteps
00349     // or if it's the last step of this pmf or mcti block,
00350     // then it's time to print
00351     if ( IsLastStep() || (((m_CurrStep-m_StartStep)%m_NumPrintSteps)==0) ) {
00352       RetVal = kTrue;
00353     }
00354   }
00355   return(RetVal);
00356 }
00357 
00358 
00359 Bool_t ALambdaControl::IsEndOf_MCTI() {
00360 //------------------------------------------------------------------------
00361 // ASSUMING that this object is currently active, decide if this is
00362 // the last time step of an mcti block.
00363 //------------------------------------------------------------------------
00364   Bool_t  RetVal=kFalse;
00365 
00366   // if an MCTI block is currently active
00367   if ( (m_Task==kStepUp)   || (m_Task==kStepDown) ||
00368        (m_Task==kStepGrow) || (m_Task==kStepFade) ) {
00369     // if this is the last step
00370     if (m_CurrStep == GetLastStep()) {
00371       RetVal = kTrue;
00372     }
00373   }
00374   return(RetVal);
00375 }
00376 
00377 
00378 Bool_t ALambdaControl::IsEndOf_MCTI_Step() {
00379 //------------------------------------------------------------------------
00380 // ASSUMING that this object is currently active, decide if this is
00381 // the last time step of an mcti step.
00382 //------------------------------------------------------------------------
00383   Bool_t  RetVal=kFalse;
00384 
00385   // if an MCTI block is currently active
00386   if ( (m_Task==kStepUp)   || (m_Task==kStepDown) ||
00387        (m_Task==kStepGrow) || (m_Task==kStepFade) ) {
00388     // if this is the last step of accumulation
00389     if (((m_CurrStep-m_StartStep)%(m_NumEquilSteps+m_NumAccumSteps))==0) {
00390       // then this is the last time step of this mcti step
00391       RetVal = kTrue;
00392     }
00393   }
00394   return(RetVal);
00395 }
00396 
00397 
00398 Bool_t ALambdaControl::IsTimeToPrint_dU_dLambda() {
00399 //------------------------------------------------------------------------
00400 // ASSUMING that this object is currently active, decide if it's time
00401 // to print out dU/dLambda information.
00402 //------------------------------------------------------------------------
00403   Bool_t  RetVal=kFalse;
00404 
00405   // if printing is required
00406   if (m_NumPrintSteps > 0) {
00407     // if number-of-steps from StartStep is an even multiple of NumPrintSteps
00408     // or if it's the last step of this pmf or mcti block,
00409     // then it might be time to print du/dLambda
00410     if ( IsLastStep() || (((m_CurrStep-m_StartStep)%m_NumPrintSteps)==0) ) {
00411       // for mcti blocks
00412       if ((m_Task==kStepUp)   || (m_Task==kStepDown) ||
00413           (m_Task==kStepGrow) || (m_Task==kStepFade)) {
00414         // it's only time to print if we're no longer equilibrating
00415         if ( ((m_CurrStep-m_StartStep-1) % (m_NumEquilSteps+m_NumAccumSteps))
00416              > m_NumEquilSteps) {
00417           RetVal = kTrue;
00418         }
00419       }
00420       // for other blocks (up, down, grow, fade) it's time
00421       else {
00422         RetVal = kTrue;
00423       }
00424     }
00425   }
00426   return(RetVal);
00427 }
00428 
00429 
00430 void ALambdaControl::PrintLambdaHeader(double dT) {
00431 //----------------------------------------------------------------------------
00432 // ASSUMING that this object is currently active, write out the header for
00433 // a new lambda control object
00434 //----------------------------------------------------------------------------
00435   // double Time;
00436   
00437   // calculate current time in femto-seconds
00438   // Time = (double)m_CurrStep * dT;
00439 
00440 iout << "FreeEnergy: ";
00441 #if !defined(_VERBOSE_PMF)
00442   iout << "nstep  time(ps)  ";
00443   iout << "    task  lambdaKf  lambdaRef     delta-G  #steps  n*{value  target |}" << std::endl;
00444   iout << "FreeEnergy: -----  --------  ";
00445   iout << "--------  --------  ---------  ----------  ------  -------------------" << std::endl;
00446   iout << endi;
00447 #endif
00448 }
00449 
00450 
00451 void ALambdaControl::PrintHeader(double dT) {
00452 //----------------------------------------------------------------------------
00453 // ASSUMING that this object is currently active, write out the current time
00454 //----------------------------------------------------------------------------
00455   double  Time;
00456   char    Str[100], Str2[100];
00457 
00458   // calculate current time in femto-seconds
00459   Time = (double)m_CurrStep * dT;
00460 
00461 #if defined(_VERBOSE_PMF)
00462   iout << "FreeEnergy: " << std::endl << endi;
00463   iout << "FreeEnergy: ";
00464   iout << "Time Step = "  << m_CurrStep            <<    ",  ";
00465   iout << "Time = ";
00466   // write out time in ps
00467   iout << Time/1000.0     << " ps,  ";
00468   iout << "Lambda_Kf = "  << m_LambdaKf            <<    ",  ";
00469   iout << "Lambda_Ref = " << m_LambdaRef           <<     "  ";
00470   GetTaskStr(Str);
00471   iout << "(" << Str << ")";
00472   iout << std::endl << endi;
00473   iout << "FreeEnergy: ";
00474   iout << "------------------------------------------------";
00475   iout << "-------------------";
00476   iout << std::endl << endi;
00477 #else
00478   sprintf(Str, "%5d", m_CurrStep);
00479   // write out time in ps
00480   sprintf(Str2, "%8.3f", Time/1000.0);
00481   iout << "FreeEnergy: ";
00482   iout << Str << "  " << Str2 << "  ";
00483   GetPaddedTaskStr(Str);
00484   iout << Str << "  ";
00485   sprintf(Str, "%8.5f", m_LambdaKf);
00486   iout << Str << "  ";
00487   sprintf(Str, "%9.5f", m_LambdaRef);
00488   iout << Str << "  ";
00489 #endif
00490 }
00491 
00492 
00493 Bool_t ALambdaControl::IsActive() {
00494 //------------------------------------------------------------------------
00495 // determine if this object is currently active
00496 //------------------------------------------------------------------------
00497   if ( (m_CurrStep>=m_StartStep) && (m_CurrStep<=GetLastStep()) ) {
00498     return(kTrue);
00499   }
00500   else {
00501     return(kFalse);
00502   }
00503 }
00504 
00505 
00506 double ALambdaControl::GetLambdaKf() {
00507 //------------------------------------------------------------------------
00508 // calculate LambdaKf for grow, fade, stepgrow, and stepfade.
00509 // LambdaKf=1.0, for up, down, stepup, stepdown, and stop.
00510 // for nogrow, LambdaKf is Lambda from the config file.
00511 //------------------------------------------------------------------------
00512   int     N;
00513   double  CurrStep  = m_CurrStep;
00514   double  StartStep = m_StartStep;
00515   double  NumSteps  = m_NumSteps;
00516   double  NumEquilSteps = m_NumEquilSteps;
00517   double  NumAccumSteps = m_NumAccumSteps;
00518   double  NumRepeats = m_NumRepeats;
00519 
00520   if (IsActive()) {
00521     switch (m_Task) {
00522       case kGrow:
00523         m_LambdaKf = (CurrStep-StartStep)/NumSteps;
00524         break;
00525       case kFade:
00526         m_LambdaKf = 1.0-(CurrStep-StartStep)/NumSteps;
00527         break;
00528       case kStepGrow:
00529         N = (int) ( (CurrStep-StartStep-1) / (NumEquilSteps+NumAccumSteps) );
00530         m_LambdaKf = (N+1)/NumRepeats;
00531         break;
00532       case kStepFade:
00533         N = (int) ( (CurrStep-StartStep-1) / (NumEquilSteps+NumAccumSteps) );
00534         m_LambdaKf = 1.0 - (N+1)/NumRepeats;
00535         break;
00536       case kNoGrow:
00537         break;              // return prior setting of m_LambdaKf
00538       default:
00539         m_LambdaKf=1.0;
00540     }
00541   }
00542   else {
00543     m_LambdaKf=1.0;
00544   }
00545   return(m_LambdaKf);
00546 }
00547 
00548 
00549 double ALambdaControl::GetLambdaRef() {
00550 //------------------------------------------------------------------------
00551 // calculate LambdaRef for up, down, stepup, and stepdown.
00552 // for stop, LambdaRef is Lambda from the config file.
00553 // for grow, fade, stepgrow, stepfade, and nogrow,
00554 //   LambdaRef is LambdaT from the config file.
00555 //------------------------------------------------------------------------
00556   int     N;
00557   double  CurrStep  = m_CurrStep;
00558   double  StartStep = m_StartStep;
00559   double  NumSteps  = m_NumSteps;
00560   double  NumEquilSteps = m_NumEquilSteps;
00561   double  NumAccumSteps = m_NumAccumSteps;
00562   double  NumRepeats = m_NumRepeats;
00563 
00564   if (IsActive()) {
00565     switch (m_Task) {
00566       case kUp:
00567         m_LambdaRef = (CurrStep-StartStep)/NumSteps;
00568         break;
00569       case kDown:
00570         m_LambdaRef = 1.0-(CurrStep-StartStep)/NumSteps;
00571         break;
00572       case kStepUp:
00573         N = (int) ( (CurrStep-StartStep-1) / (NumEquilSteps+NumAccumSteps) );
00574         m_LambdaRef = (N+1)/NumRepeats;
00575         break;
00576       case kStepDown:
00577         N = (int) ( (CurrStep-StartStep-1) / (NumEquilSteps+NumAccumSteps) );
00578         m_LambdaRef = 1.0 - (N+1)/NumRepeats;
00579       default: 
00580         break;             // return prior setting of m_LambdaRef
00581     }
00582   }
00583   else {
00584     m_LambdaRef=0.0;
00585   }
00586   return(m_LambdaRef);
00587 }
00588 

Generated on Fri Sep 22 01:17:12 2017 for NAMD by  doxygen 1.4.7