Add alchIDWSFreq fallback to outputEnergies 70/4770/7
authorradakb <brian.radak@gmail.com>
Fri, 2 Nov 2018 16:35:52 +0000 (12:35 -0400)
committerJim Phillips <jim@ks.uiuc.edu>
Wed, 7 Nov 2018 20:01:48 +0000 (14:01 -0600)
If alchOutFreq = 0, fallback to alchIDWSFreq = outputEnergies. This
enables IDWS output to stdout without sending any output to
alchOutFile.

Also, because the value of alchLambda2 is step-dependent (not event
driven), the use of two IDWS output streams may be inherently
incompatible if the output frequencies are different (outputEnergies
must be an odd multiple of alchOutFreq). Users can still use whatever
settings they would like, but there is a warning issued when
alchOutFreq > 0 and alchOutFreq != outputEnergies.

Minor change:

alchIDWSfreq --> alchIDWSFreq, to match camel casing convention
throughout the rest of the code

Change-Id: I5d6611b6eb9374cd1cc97e1e87169e4beba83322

src/Controller.C
src/SimParameters.C
src/SimParameters.h

index 7908a19..8ef5c41 100644 (file)
@@ -3630,7 +3630,7 @@ void Controller::outputFepEnergy(int step) {
        electEnergy - electEnergySlow - ljEnergy;
   BigReal RT = BOLTZMANN * simParams->alchTemp;
 
-  if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. || (step / simParams->alchIDWSfreq) % 2 == 1 )){
+  if (alchEnsembleAvg && (simParams->alchLambdaIDWS < 0. || (step / simParams->alchIDWSFreq) % 2 == 1 )){
     // with IDWS, only accumulate stats on those timesteps where target lambda is "forward"
     FepNo++;
     exp_dE_ByRT += exp(-dE/RT);
@@ -3870,7 +3870,7 @@ void Controller::writeFepEnergyData(int step, ofstream_namd &file) {
 
   if(stepInRun){
     if(!FepWhamOn){
-      if ( simParams->alchLambdaIDWS >= 0. && (step / simParams->alchIDWSfreq) % 2 == 0 ) {
+      if ( simParams->alchLambdaIDWS >= 0. && (step / simParams->alchIDWSFreq) % 2 == 0 ) {
         // IDWS is active and we are on a "backward" timestep
         fepFile << FEPTITLE_BACK(step);
       } else {
index f19db39..06b15f2 100644 (file)
@@ -362,15 +362,7 @@ void SimParameters::scriptSet(const char *param, const char *value) {
 
   if ( ! strncasecmp(param,"alchLambdaIDWS",MAX_SCRIPT_PARAM_SIZE) ) {
     alchLambdaIDWS = atof(value);
-    if ( alchLambdaIDWS > 1. ) {
-      NAMD_die("alchLambdaIDWS should be either in the range [0.0, 1.0], or negative (disabled).\n");
-    }
-    // Switch lambda2 every other cycle of fullElectFrequency steps
-    // or every other nonbondedFrequency steps if undefined
-    // or every alchOutFreq steps if larger (no need to switch faster that we output)
-    alchIDWSfreq = fullElectFrequency > 0 ? fullElectFrequency : nonbondedFrequency;
-    if ( alchOutFreq > alchIDWSfreq )
-      alchIDWSfreq = alchOutFreq;
+    setupIDWS();
     ComputeNonbondedUtil::select();
     return;
   }
@@ -3582,6 +3574,8 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
        if (alchLambda2 < 0.0 || alchLambda2 > 1.0)
          NAMD_die("alchLambda2 values should be in the range [0.0, 1.0]\n");
 
+       setupIDWS(); // setup IDWS if it was activated.
+       
        if (!opts.defined("alchoutfile")) {
          strcpy(alchOutFile, outputFilename);
          strcat(alchOutFile, ".fep");
@@ -3649,17 +3643,6 @@ void SimParameters::check_config(ParseOptions &opts, ConfigList *config, char *&
          strcat(alchOutFile, ".ti"); 
        } 
      }
-     if (alchLambdaIDWS >= 0.) {
-       if ( alchLambdaIDWS > 1. ) {
-         NAMD_die("alchLambdaIDWS should be either in the range [0.0, 1.0], or negative (disabled).\n");
-      }
-       // Switch lambda2 every other cycle of fullElectFrequency steps
-       // or every other nonbondedFrequency steps if undefined
-       // or every alchOutFreq steps if larger (no need to switch faster that we output)
-       alchIDWSfreq = fullElectFrequency > 0 ? fullElectFrequency : nonbondedFrequency;
-       if ( alchOutFreq > alchIDWSfreq )
-         alchIDWSfreq = alchOutFreq;
-     }
    }
         
 //fepe
@@ -7172,14 +7155,40 @@ void SimParameters::receive_SimParameters(MIStream *msg)
 //fepb IDWS
 BigReal SimParameters::getCurrentLambda2(const int step) {
   if ( alchLambdaIDWS >= 0. ) {
-    const BigReal lambda2 = ( (step / alchIDWSfreq) % 2 == 1 ) ? alchLambda2 : alchLambdaIDWS;
+    const BigReal lambda2 = ( (step / alchIDWSFreq) % 2 == 1 ) ? alchLambda2 : alchLambdaIDWS;
     return lambda2;
   } else {
     return alchLambda2;
   }
 }
-//fepe IDWS
 
+/* Return true if IDWS is active, else return false. */
+int SimParameters::setupIDWS() {
+  if (alchLambdaIDWS < 0.) return 0;
+  if (alchLambdaIDWS > 1.) {
+    NAMD_die("alchLambdaIDWS should be either in the range [0.0, 1.0], or negative (disabled).\n");
+  }
+ /* 
+  * The internal parameter alchIDWSFreq determines the number of steps of MD
+  * before each switch of the value of alchLambda2. At most this occurs every
+  * time the energy is evaluated and thus the default is the greater of
+  * fullElectFrequency and nonbondedFrequency. However, this choice fails to
+  * report alternating values if output is printed less often than every step
+  * (which is almost certainly true). Thus the frequency is reset to match
+  * alchOutFreq or, if that is zero, outputEnergies. Note that, if 
+  * alchOutFreq > 0 but != outputEnergies, then the data going to stdout
+  * are likely not useful since the comparison value is difficult to infer.
+  */
+  alchIDWSFreq = fullElectFrequency > 0 ? fullElectFrequency : nonbondedFrequency;
+  if ( !alchOutFreq && outputEnergies > alchIDWSFreq ) alchIDWSFreq = outputEnergies;
+  if ( alchOutFreq > alchIDWSFreq ) alchIDWSFreq = alchOutFreq;
+  if ( alchOutFreq && alchOutFreq != outputEnergies) {
+    iout << iWARN << "alchOutFreq and outputEnergies do not match. IDWS ouput"
+         << " to stdout may not be useful!\n" << endi;
+  }
+  return 1;
+}
+//fepe IDWS
 
 //fepb BKR
 BigReal SimParameters::getCurrentLambda(const int step) {
index 4954da7..dcaaca8 100644 (file)
@@ -386,11 +386,12 @@ public:
   BigReal alchLambda;       //  lambda for dynamics
   BigReal alchLambda2;      //  lambda for comparison
   BigReal alchLambdaIDWS;   //  alternate lambda for interleaved double-wide sampling
-  int alchIDWSfreq;         //  freq with which lambda2 changes to lambdaIDWS
+  int alchIDWSFreq;         //  freq with which lambda2 changes to lambdaIDWS
   int alchLambdaFreq;       //  freq. (in steps) with which lambda changes 
                             //  from alchLambda to alchLambda2
   BigReal getCurrentLambda(const int); // getter for changing lambda
   BigReal getCurrentLambda2(const int); // getter for alternating lambda2 in IDWS
+  int setupIDWS();          //  activates IDWS and sets alchIDWSFreq
   BigReal getLambdaDelta(void); // getter for lambda increment
   BigReal alchRepLambda;    //  lambda for WCA repulsive interaction
   BigReal alchDispLambda;   //  lambda for WCA dispersion interaction