Added mdp option 'calc-lambda-neighbors'.

Added an option 'calc-lambda-neighbors' to limit the number of free
energy delta H evaluations and output operations: this value sets the
amount of neighboring lambda points to calculate and print the energy
for. If it's 1, it will only calculate directly adjacent lambda points'
energies (if it's -1, it will calculate all points). This tremendously
cuts down the number of (expensive) energy calculations done without
impacting the results for most users

Note that this commit only limits the number of outputs, not the number
of calculations. That is left for a later commit.

Change-Id: I10590dd1a2036b529946da86b414c2561d14bda1

diff --git a/include/types/inputrec.h b/include/types/inputrec.h
index 0a67146bdf8fd897fe52dcf5d204e50ea90449f6..359ed0e6f64512606c79b3b5a0f9650a6eeca5dc 100644 (file)
--- a/include/types/inputrec.h
+++ b/include/types/inputrec.h
@@ -142,6 +142,11 @@ typedef struct {
   gmx_bool bPrintEnergy; /* Whether to print the energy in the dhdl           */
   int  n_lambda;         /* The number of foreign lambda points               */
   double **all_lambda;   /* The array of all lambda values                    */
+  int lambda_neighbors;  /* The number of neighboring lambda states to
+                            calculate the energy for in up and down directions
+                            (-1 for all) */
+  int lambda_start_n;    /* The first lambda to calculate energies for */
+  int lambda_stop_n;     /* The last lambda +1 to calculate energies for */
   real sc_alpha;         /* free energy soft-core parameter                   */
   int  sc_power;         /* lambda power for soft-core interactions           */
   real sc_r_power;          /* r power for soft-core interactions                */

diff --git a/share/html/online/mdp_opt.html b/share/html/online/mdp_opt.html
index 804bc6d4547ce13f66a0246b46bbb179f9a9cfdf..86adb78f147103372c740579dc99d59e69d7a8c3 100644 (file)
--- a/share/html/online/mdp_opt.html
+++ b/share/html/online/mdp_opt.html
@@ -1529,6 +1529,15 @@ Only the particle masses are controlled with this component of the lambda vector
 be determined and written to dhdl.xvg every <b>nstdhdl</b> steps. Values must be between 0 and 1.
 Only the temperatures controlled with this component of the lambda vector.
 Note that these lambdas should not be used for replica exchange, only for simulated tempering.</dd>
+<dt><b>calc-lambda-neighbors (1)</b></dt>
+<dd>Controls the number of lambda values for which Delta H values will be
+calculated and written out, if <b>init-lambda-state</b> has been set. A
+positive value will limit the number of lambda points calculated to only the
+nth neighbors of <b>init-lambda-state</b>: for example, if
+<b>init-lambda-state</b> is 5 and this parameter has a value of 2, energies for
+lambda points 3-7 will be calculated and writen out. A value of -1 means all
+lambda points will be written out. For normal BAR such as with g_bar, a value
+of 1 is sufficient, while for MBAR -1 should be used.</dd>
 <dt><b>sc-alpha: (0)</b></dt>
 <dd>the soft-core alpha parameter, a value of 0 results in linear interpolation of the LJ and Coulomb interactions</dd>
 <dt><b>sc-r-power: (6)</b></dt>

diff --git a/src/gmxlib/tpxio.c b/src/gmxlib/tpxio.c
index c5a74c515b61e6512667f3a351b579bf805efae3..717c000919c006ed16b18f6f126501fd068c1a76 100644 (file)
--- a/src/gmxlib/tpxio.c
+++ b/src/gmxlib/tpxio.c
@@ -75,7 +75,7 @@
 static const char *tpx_tag = TPX_TAG_RELEASE;
 
 /* This number should be increased whenever the file format changes! */
-static const int tpx_version = 82;
+static const int tpx_version = 83;
 
 /* This number should only be increased when you edit the TOPOLOGY section
  * or the HEADER of the tpx format.
@@ -252,14 +252,18 @@ static void do_pullgrp(t_fileio *fio, t_pullgrp *pgrp, gmx_bool bRead,
   }
 }
 
-static void do_expandedvals(t_fileio *fio,t_expanded *expand,int n_lambda, gmx_bool bRead, int file_version)
+static void do_expandedvals(t_fileio *fio,t_expanded *expand,t_lambda *fepvals, gmx_bool bRead, int file_version)
 {
   /* i is used in the ndo_double macro*/
   int i;
   real fv;
   gmx_bool bDum=TRUE;
   real rdum;
+  int n_lambda=fepvals->n_lambda;
 
+  /* reset the lambda calculation window */
+  fepvals->lambda_start_n = 0;
+  fepvals->lambda_stop_n = n_lambda;
   if (file_version >= 79)
   {
       if (n_lambda>0)
@@ -327,6 +331,7 @@ static void do_fepvals(t_fileio *fio,t_lambda *fepvals,gmx_bool bRead, int file_
   real rdum;
 
   /* free energy values */
+
   if (file_version >= 79)
   {
       gmx_fio_do_int(fio,fepvals->init_fep_state);
@@ -501,6 +506,38 @@ static void do_fepvals(t_fileio *fio,t_lambda *fepvals,gmx_bool bRead, int file_
   {
       fepvals->bPrintEnergy = FALSE;
   }
+
+  /* handle lambda_neighbors */
+  if (file_version >= 83)
+  {
+      gmx_fio_do_int(fio,fepvals->lambda_neighbors);
+      if ( (fepvals->lambda_neighbors >= 0) && (fepvals->init_fep_state>=0) &&
+           (fepvals->init_lambda < 0) )
+      {
+          fepvals->lambda_start_n = (fepvals->init_fep_state -
+                                     fepvals->lambda_neighbors);
+          fepvals->lambda_stop_n = (fepvals->init_fep_state +
+                                    fepvals->lambda_neighbors + 1);
+          if (fepvals->lambda_start_n < 0)
+          {
+              fepvals->lambda_start_n = 0;;
+          }
+          if (fepvals->lambda_stop_n >= fepvals->n_lambda)
+          {
+              fepvals->lambda_stop_n = fepvals->n_lambda;
+          }
+      }
+      else
+      {
+          fepvals->lambda_start_n = 0;
+          fepvals->lambda_stop_n = fepvals->n_lambda;
+      }
+  }
+  else
+  {
+      fepvals->lambda_start_n = 0;
+      fepvals->lambda_stop_n = fepvals->n_lambda;
+  }
 }
 
 static void do_pull(t_fileio *fio, t_pull *pull,gmx_bool bRead, int file_version)
@@ -987,7 +1024,7 @@ static void do_inputrec(t_fileio *fio, t_inputrec *ir,gmx_bool bRead,
     }
     if (ir->bExpanded)
     {
-        do_expandedvals(fio,ir->expandedvals,ir->fepvals->n_lambda,bRead,file_version);
+        do_expandedvals(fio,ir->expandedvals,ir->fepvals,bRead,file_version);
     }
     if (file_version >= 57) {
       gmx_fio_do_int(fio,ir->eDisre); 

diff --git a/src/gmxlib/txtdump.c b/src/gmxlib/txtdump.c
index 39906a2fdd200ce147cb545b024553dee6b5af9d..77d3980cb0fd7e06f476056e3bec3e9f23293ae2 100644 (file)
--- a/src/gmxlib/txtdump.c
+++ b/src/gmxlib/txtdump.c
@@ -587,6 +587,7 @@ static void pr_fepvals(FILE *fp,int indent,t_lambda *fep, gmx_bool bMDPformat)
             fprintf(fp,"\n");
         }
     }
+    PI("calc-lambda-neighbors",fep->lambda_neighbors);
 
     PR("sc-alpha",fep->sc_alpha);
     PS("bScCoul",EBOOL(fep->bScCoul));

diff --git a/src/kernel/readir.c b/src/kernel/readir.c
index d47cef5dd2891456528dac161f69674be909579e..12b43489e746095b40c92b83085fad16a33a9a8f 100644 (file)
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@ -1914,6 +1914,7 @@ void get_ir(const char *mdparin,const char *mdparout,
   STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
   STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
   STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
+  ITYPE ("calc-lambda-neighbors",fep->lambda_neighbors, 1);
   STYPE ("init-lambda-weights",lambda_weights,NULL);
   EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
   RTYPE ("sc-alpha",fep->sc_alpha,0.0);

diff --git a/src/kernel/tpbcmp.c b/src/kernel/tpbcmp.c
index 48901e318063db8916e71f0a4a27e0c2192c3f4c..8d41b8a9cfad48ed923a7120edef712340a73422 100644 (file)
--- a/src/kernel/tpbcmp.c
+++ b/src/kernel/tpbcmp.c
@@ -606,6 +606,8 @@ static void cmp_fepvals(FILE *fp,t_lambda *fep1,t_lambda *fep2,real ftol, real a
           cmp_double(fp,"inputrec->fepvals->all_lambda",-1,fep1->all_lambda[i][j],fep2->all_lambda[i][j],ftol,abstol);
       }
   }
+  cmp_int(fp,"inputrec->fepvals->lambda_neighbors",1,fep1->lambda_neighbors,
+          fep2->lambda_neighbors);
   cmp_real(fp,"inputrec->fepvals->sc_alpha",-1,fep1->sc_alpha,fep2->sc_alpha,ftol,abstol);
   cmp_int(fp,"inputrec->fepvals->sc_power",-1,fep1->sc_power,fep2->sc_power);
   cmp_real(fp,"inputrec->fepvals->sc_r_power",-1,fep1->sc_r_power,fep2->sc_r_power,ftol,abstol);

diff --git a/src/mdlib/mdebin.c b/src/mdlib/mdebin.c
index d988e180e1db23603180b5a195fe978e19c71fe6..266cc19ff7a2f990e6e99d8b15dba8c31e2dfe39 100644 (file)
--- a/src/mdlib/mdebin.c
+++ b/src/mdlib/mdebin.c
@@ -745,7 +745,7 @@ extern FILE *open_dhdl(const char *filename,const t_inputrec *ir,
         nsets_dhdl = n_lambda_terms;
     }
     /* count the number of delta_g states */
-    nsets_de = fep->n_lambda;
+    nsets_de = fep->lambda_stop_n - fep->lambda_start_n;
 
     nsets = nsets_dhdl + nsets_de; /* dhdl + fep differences */
 
@@ -831,7 +831,7 @@ extern FILE *open_dhdl(const char *filename,const t_inputrec *ir,
         }
         nsetsbegin += nsets_dhdl;
 
-        for(i=0; i<fep->n_lambda; i++)
+        for(i=fep->lambda_start_n; i<fep->lambda_stop_n; i++)
         {
             print_lambda_vector(fep, i, FALSE, FALSE, lambda_vec_str);
             if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))
@@ -1135,7 +1135,7 @@ void upd_mdebin(t_mdebin *md,
                     }
                 }
             }
-            for(i=0;i<fep->n_lambda;i++)
+            for(i=fep->lambda_start_n;i<fep->lambda_stop_n;i++)
             {
                 fprintf(md->fp_dhdl," %#.8g",md->dE[i]);
             }
@@ -1170,7 +1170,7 @@ void upd_mdebin(t_mdebin *md,
                                     store_energy,
                                     pv,
                                     store_dhdl,
-                                    md->dE,
+                                    md->dE + fep->lambda_start_n,
                                     time);
         }
     }

diff --git a/src/mdlib/mdebin_bar.c b/src/mdlib/mdebin_bar.c
index bc62faff1f62affadba00f4d6d1d2979c65a00e8..2b1ccf5d4abdca4d69b7e7c9caf3d232d76375f8 100644 (file)
--- a/src/mdlib/mdebin_bar.c
+++ b/src/mdlib/mdebin_bar.c
@@ -446,7 +446,7 @@ void mde_delta_h_coll_init(t_mde_delta_h_coll *dhc, const t_inputrec *ir)
             }
         }
         /* add the lambdas */
-        dhc->nlambda = ir->fepvals->n_lambda;
+        dhc->nlambda = ir->fepvals->lambda_stop_n - ir->fepvals->lambda_start_n;
         dhc->ndh += dhc->nlambda;
         /* another compatibility check */
         if (dhc->start_lambda < 0)
@@ -523,7 +523,7 @@ void mde_delta_h_coll_init(t_mde_delta_h_coll *dhc, const t_inputrec *ir)
         /* add the lambdas */
         dhc->dh_du = dhc->dh + n;
         snew(lambda_vec, n_lambda_components);
-        for(i=0;i<ir->fepvals->n_lambda;i++)
+        for(i=ir->fepvals->lambda_start_n;i<ir->fepvals->lambda_stop_n;i++)
         {
             int k=0;