diff --git a/doc/src/Eqs/fix_grem.jpg b/doc/src/Eqs/fix_grem.jpg new file mode 100644 index 0000000000..666b4344a6 Binary files /dev/null and b/doc/src/Eqs/fix_grem.jpg differ diff --git a/doc/src/Eqs/fix_grem.tex b/doc/src/Eqs/fix_grem.tex new file mode 100644 index 0000000000..c00ea696a3 --- /dev/null +++ b/doc/src/Eqs/fix_grem.tex @@ -0,0 +1,9 @@ +\documentclass[12pt]{article} + +\begin{document} + +$$ + T_{eff} = \lambda + \eta (H - H_0) +$$ + +\end{document} \ No newline at end of file diff --git a/doc/src/fix_grem.txt b/doc/src/fix_grem.txt new file mode 100644 index 0000000000..1f62b00104 --- /dev/null +++ b/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). diff --git a/examples/USER/misc/grem/lj-single/in.gREM b/examples/USER/misc/grem/lj-single/in.gREM-npt similarity index 100% rename from examples/USER/misc/grem/lj-single/in.gREM rename to examples/USER/misc/grem/lj-single/in.gREM-npt diff --git a/src/USER-MISC/fix_grem.cpp b/src/USER-MISC/fix_grem.cpp index 339b41c616..84e2fb5929 100644 --- a/src/USER-MISC/fix_grem.cpp +++ b/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"); + } } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-MISC/fix_grem.h b/src/USER-MISC/fix_grem.h index 21a9cb65a1..c080d07efe 100644 --- a/src/USER-MISC/fix_grem.h +++ b/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; }; }