Merge branch 'grem-feature' of https://github.com/dstelter92/lammps into grem-feature

Resolved merge conflicts and adapted logic to most recent changes in feature branch

Closes #46

doc/src/Eqs/fix_grem.tex

@@ -0,0 +1,9 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+
+$$
+  T_{eff} = \lambda + \eta (H - H_0)
+$$
+
+\end{document}

doc/src/fix_grem.txt

@@ -0,0 +1,108 @@
+"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
+
+:link(lws,http://lammps.sandia.gov)
+:link(ld,Manual.html)
+:link(lc,Section_commands.html#comm)
+
+:line
+
+fix grem command :h3
+
+[Syntax:]
+
+fix ID group-ID grem lambda eta H0 thermostat-ID
+
+ID, group-ID are documented in "fix"_fix.html command
+grem = style name of this fix command
+lambda = intercept parameter of linear effective temperature function
+eta = slope parameter of linear effective temperature function
+H0 = shift parameter of linear effective temperature function
+thermostat-ID = ID of thermostat used in simulation
+
+[Examples:]
+
+fix             fxgREM all grem 400 -0.01 -30000 fxnpt
+thermo_modify   press fxgREM_press
+
+fix             fxgREM all grem 502 -0.15 -80000 fxnvt
+
+[Description:]
+
+This fix implements the molecular dynamics version
+of the generalized replica 
+exchange method (gREM) originally developed by "(Kim)"_#Kim,
+which uses non-Boltzmann ensembles to sample over first
+order phase transitions.
+
+The is done by defining replicas with an enthalpy dependent
+effective temperature
+
+:c,image(Eqs/fix_grem.jpg)
+
+with {eta} negative and steep enough to only intersect 
+the characteristic microcanonical temperature (Ts) of the
+system once, ensuring a unimodal enthalpy distribution in 
+that replica. {Lambda} is the intercept and effects the 
+generalized ensemble similar to how temperature effects 
+a Boltzmann ensemble. {H0} is a reference enthalpy, and
+is typically set as the lowest desired sampled enthalpy.
+Further explanation can be found in our recent papers 
+"(Malolepsza)"_#Malolepsza.
+
+This fix requires a thermostat, with ID passed to fix_grem
+by {thermostat-ID}. Two distinct temperatures exist in this
+generalized ensemble, the effective temperature defined above,
+and a kinetic temperature that controls the velocity 
+distribution of particles as usual. Either constant volume
+or constant pressure algorithms can be used.
+
+The fix enforces a generalized ensemble in a single replica
+only. Typically, different replicas only differ by {lambda}
+for simplicity, but this is not necessary. Multi-replica 
+runs need to be run outside of LAMMPS. An example of this 
+can be found in examples/USER/misc/grem
+
+In general, defining the generalized ensembles is unique for 
+every system. When starting a many-replica simulation without 
+any knowledge of the underlying microcanonical temperature, 
+there are several tricks we have utilized to optimize the process. 
+Choosing a less-steep {eta} yields broader distributions, 
+requiring fewer replicas to map the microcanonical temperature. 
+While this likely struggles from the same sampling problems 
+gREM was built to avoid, it provides quick insight to Ts. 
+Initially using an evenly-spaced {lambda} distribution identifies
+regions where small changes in enthalpy lead to large temperature
+changes. Replicas are easily added where needed. 
+
+:line
+
+[Restart, fix_modify, output, run start/stop, minimize info:]
+
+No information about this fix is written to "binary restart
+files"_restart.html.
+
+The "thermo_modify"_thermo_modify.html {press} option is supported
+by this fix to add the rescaled kinetic pressure as part of 
+"thermodynamic output"_thermo_style.html.
+
+[Restrictions:]
+
+This fix is part of the USER-MISC package. It is only enabled if 
+LAMMPS was built with that package. See the "Making 
+LAMMPS"_Section_start.html#start_3 section for more info.
+
+[Related commands:]
+
+"fix_nh"_fix_nh.html
+"thermo_modify"_thermo_modify.html
+
+[Default:] none
+
+:line
+
+:link(Kim)
+[(Kim)] Kim, Keyes, Straub, J. Chem. Phys., 132, 224107 (2010).
+
+:link(Malolepsza)
+[(Malolepsza)] Malolepsza, Secor, Keyes, J. Phys. Chem. B. 119 (42), 
+13379-13384 (2015).

src/USER-MISC/fix_grem.cpp

@@ -65,8 +65,8 @@ FixGrem::FixGrem(LAMMPS *lmp, int narg, char **arg) :
   h0 = force->numeric(FLERR,arg[5]);
 
   int n = strlen(arg[6])+1;
-  id_npt = new char[n];
-  strcpy(id_npt,arg[6]);
+  id_nh = new char[n];
+  strcpy(id_nh,arg[6]);
 
   // create a new compute temp style
   // id = fix-ID + temp
@@ -133,15 +133,24 @@ FixGrem::FixGrem(LAMMPS *lmp, int narg, char **arg) :
   modify->add_compute(3,newarg);
   delete [] newarg;
 
-  int ifix = modify->find_fix(id_npt);
+  int ifix = modify->find_fix(id_nh);
   if (ifix < 0)
-    error->all(FLERR,"Fix id for npt fix does not exist");
-  Fix *npt = modify->fix[ifix];
+    error->all(FLERR,"Fix id for nvt or npt fix does not exist");
+  Fix *nh = modify->fix[ifix];
 
-  char *modargs[2];
-  modargs[0] = (char *) "press";
-  modargs[1] = id_press;
-  npt->modify_param(2,modargs);
+  pressflag = 0;
+  int *p_flag = (int *)nh->extract("p_flag",ifix);
+  if ((p_flag == NULL) || (ifix != 1) || (p_flag[0] == 0)
+      || (p_flag[1] == 0) || (p_flag[2] == 0)) {
+    pressflag = 0;
+  } else if ((p_flag[0] == 1) && (p_flag[1] == 1)
+             && (p_flag[2] == 1) && (ifix == 1)) {
+    pressflag = 1;
+    char *modargs[2];
+    modargs[0] = (char *) "press";
+    modargs[1] = id_press;
+    nh->modify_param(2,modargs);
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -158,7 +167,7 @@ FixGrem::~FixGrem()
   delete [] id_press;
   delete [] id_ke;
   delete [] id_pe;
-  delete [] id_npt;
+  delete [] id_nh;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -195,36 +204,39 @@ void FixGrem::init()
     error->all(FLERR,"PE compute ID for fix grem does not exist");
   pe = modify->compute[icompute];
 
-  int ifix = modify->find_fix(id_npt);
+  int ifix = modify->find_fix(id_nh);
   if (ifix < 0)
-    error->all(FLERR,"Fix id for npt fix does not exist");
-  Fix *npt = modify->fix[ifix];
+    error->all(FLERR,"Fix id for nvt or npt fix does not exist");
+  Fix *nh = modify->fix[ifix];
 
-  double *t_start = (double *)npt->extract("t_start",ifix);
-  double *t_stop = (double *)npt->extract("t_stop",ifix);
+  double *t_start = (double *)nh->extract("t_start",ifix);
+  double *t_stop = (double *)nh->extract("t_stop",ifix);
   if ((t_start != NULL) && (t_stop != NULL) && (ifix == 0)) {
     tbath = *t_start;
     if (*t_start != *t_stop)
       error->all(FLERR,"Thermostat temperature ramp not allowed");
   } else
-    error->all(FLERR,"Problem extracting target temperature from fix npt");
+    error->all(FLERR,"Problem extracting target temperature from fix nvt or npt");
 
-  int *p_flag = (int *)npt->extract("p_flag",ifix);
-  double *p_start = (double *) npt->extract("p_start",ifix);
-  double *p_stop = (double *) npt->extract("p_stop",ifix);
-  if ((p_flag != NULL) && (p_start != NULL) && (p_stop != NULL)
-      && (ifix == 1)) {
-    ifix = 0;
-    pressref = p_start[0];
-    if ((p_start[0] != p_stop[0]) || (p_flag[0] != 1)) ++ ifix;
-    if ((p_start[1] != p_stop[1]) || (p_flag[0] != 1)) ++ ifix;
-    if ((p_start[2] != p_stop[2]) || (p_flag[0] != 1)) ++ ifix;
-    if ((p_start[0] != p_start[1]) || (p_start[1] != p_start[2])) ++ifix;
-    if ((p_flag[3] != 0) || (p_flag[4] != 0) || (p_flag[5] != 0)) ++ifix;
-    if (ifix > 0)
-      error->all(FLERR,"Unsupported pressure settings in fix npt");
-  } else
-    error->all(FLERR,"Problem extracting target pressure from fix npt");
+  pressref = 0.0;
+  if (pressflag) {
+    int *p_flag = (int *)nh->extract("p_flag",ifix);
+    double *p_start = (double *) nh->extract("p_start",ifix);
+    double *p_stop = (double *) nh->extract("p_stop",ifix);
+    if ((p_flag != NULL) && (p_start != NULL) && (p_stop != NULL)
+        && (ifix == 1)) {
+      ifix = 0;
+      pressref = p_start[0];
+      if ((p_start[0] != p_stop[0]) || (p_flag[0] != 1)) ++ ifix;
+      if ((p_start[1] != p_stop[1]) || (p_flag[0] != 1)) ++ ifix;
+      if ((p_start[2] != p_stop[2]) || (p_flag[0] != 1)) ++ ifix;
+      if ((p_start[0] != p_start[1]) || (p_start[1] != p_start[2])) ++ifix;
+      if ((p_flag[3] != 0) || (p_flag[4] != 0) || (p_flag[5] != 0)) ++ifix;
+      if (ifix > 0)
+        error->all(FLERR,"Unsupported pressure settings in fix npt");
+    } else
+      error->all(FLERR,"Problem extracting target pressure from fix npt");
+  }
 }
 
 /* ---------------------------------------------------------------------- */

src/USER-MISC/fix_grem.h

@@ -40,9 +40,9 @@ class FixGrem : public Fix {
   double lambda,eta,h0,tbath,pressref;
 
  protected:
-  char *id_temp,*id_press,*id_ke,*id_pe,*id_npt;
+  char *id_temp,*id_press,*id_ke,*id_pe,*id_nh;
   class Compute *temperature,*pressure,*ke,*pe;
-  int pflag,tflag,keflag,peflag;
+  int pressflag;
 };
 
 }