Merge branch 'master' into fixes-for-stable

doc/src/Eqs/fix_gcmc1.tex

@@ -0,0 +1,9 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+
+\begin{eqnarray*}
+\mu &=&\mu^{id} + \mu^{ex}
+\end{eqnarray*}
+
+\end{document}

doc/src/Eqs/fix_gcmc2.tex

@@ -0,0 +1,10 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+
+\begin{eqnarray*}
+\mu^{id} &=& k T \ln{\rho \Lambda^3} \\
+&=& k T \ln{\frac{\phi P \Lambda^3}{k T}} 
+\end{eqnarray*}
+
+\end{document}

doc/src/Eqs/fix_gcmc3.tex

@@ -0,0 +1,9 @@
+\documentclass[12pt]{article}
+
+\begin{document}
+
+\begin{eqnarray*}
+\Lambda &=& \sqrt{ \frac{h^2}{2 \pi m k T}}
+\end{eqnarray*}
+
+\end{document}

doc/src/dihedral_charmm.txt

@@ -40,8 +40,11 @@ field.
 
 NOTE: The newer {charmmfsh} style was released in March 2017.  We
 recommend it be used instead of the older {charmm} style when running
-a simulation with the CHARMM force field.  See the discussion below
-and more details on the "pair_style charmm"_pair_charmm.html doc page.
+a simulation with the CHARMM force field and Coulomb cutoffs, via the
+"pair_style lj/charmmfsw/coul/charmmfsh"_pair_charmm.html command.
+Otherwise the older {charmm} style is fine to use.  See the discussion
+below and more details on the "pair_style charmm"_pair_charmm.html doc
+page.
 
 The following coefficients must be defined for each dihedral type via the
 "dihedral_coeff"_dihedral_coeff.html command as in the example above, or in
@@ -82,13 +85,18 @@ special_bonds 1-4 scaling factor to 0.0 (which is the
 default). Otherwise 1-4 non-bonded interactions in dihedrals will be
 computed twice.
 
-For simulations using the CHARMM force field, the difference between
-the {charmm} and {charmmfsh} styles is in the computation of the 1-4
-non-bond interactions, if the distance between the two atoms is within
-the switching distance of the pairwise potential defined by the
-corresponding CHARMM pair style, i.e. between the inner and outer
-cutoffs specified for the pair style.  See the discussion on the
-"CHARMM pair_style"_pair_charmm.html doc page for details.
+For simulations using the CHARMM force field with a Coulomb cutoff,
+the difference between the {charmm} and {charmmfsh} styles is in the
+computation of the 1-4 non-bond interactions, though only if the
+distance between the two atoms is within the switching region of the
+pairwise potential defined by the corresponding CHARMM pair style,
+i.e. between the inner and outer cutoffs specified for the pair style.
+The {charmmfsh} style should only be used when using the "pair_style
+lj/charmmfsw/coul/charmmfsh"_pair_charmm.html to make the Coulombic
+pairwise calculations consistent.  Use the {charmm} style with
+long-range Coulombics or the older "pair_style
+lj/charmm/coul/charmm"_pair_charmm.html command.  See the discussion
+on the "CHARMM pair_style"_pair_charmm.html doc page for details.
 
 Note that for AMBER force fields, which use pair styles with "lj/cut",
 the special_bonds 1-4 scaling factor should be set to the AMBER

doc/src/fix_gcmc.txt

@@ -122,11 +122,11 @@ If used with "fix nvt"_fix_nh.html, the temperature of the imaginary
 reservoir, T, should be set to be equivalent to the target temperature
 used in fix nvt. Otherwise, the imaginary reservoir
 will not be in thermal equilibrium with the simulation cell. Also,
-it is important that the temperature used by fix nvt be dynamic,
+it is important that the temperature used by fix nvt be dynamic/dof,
 which can be achieved as follows:
 
 compute mdtemp mdatoms temp
-compute_modify mdtemp dynamic yes
+compute_modify mdtemp dynamic/dof yes
 fix mdnvt mdatoms nvt temp 300.0 300.0 10.0
 fix_modify mdnvt temp mdtemp :pre
 
@@ -204,10 +204,43 @@ atoms/molecules are assigned to two groups: the default group "all"
 and the group specified in the fix gcmc command (which can also be
 "all").
 
-The gas reservoir pressure can be specified using the {pressure}
+The chemical potential is a user-specified input parameter defined
+as:
+
+:c,image(Eqs/fix_gcmc1.jpg)
+
+The second term mu_ex is the excess chemical potential due to
+energetic interactions and is formally zero for the fictitious gas
+reservoir but is non-zero for interacting systems. So, while the chemical
+potential of the reservoir and the simulation cell are equal, mu_ex is not,
+and as a result, the densities of the two are generally quite different.
+The first term mu_id is the ideal gas contribution to the chemical potential.
+mu_id can be related to the density or pressure of the fictitious gas reservoir by:
+
+:c,image(Eqs/fix_gcmc2.jpg)
+
+where k is Boltzman's constant,
+T is the user-specified temperature, rho is the number density,
+P is the pressure, and phi is the fugacity coefficient.
+The constant Lambda is required for dimensional consistency.
+For all unit styles except {lj} it is defined as the thermal
+de Broglie wavelength
+
+:c,image(Eqs/fix_gcmc3.jpg)
+
+where h is Planck's constant, and m is the mass of the exchanged atom or molecule.
+For unit style {lj}, Lambda is simply set to the unity. Note that prior to March 2017
+Lambda for unit style {lj} was calculated using the above formula with h set to
+the rather specific value of 0.18292026. Chemical potential
+under the old definition can be converted to an equivalent value under the new
+definition by subtracting 3kTln(Lambda_old).
+
+As an alternative to specifying mu directly, the ideal gas reservoir
+can be defined by its pressure P using the {pressure}
 keyword, in which case the user-specified chemical potential is
-ignored. For non-ideal gas reservoirs, the user may also specify the
-fugacity coefficient using the {fugacity_coeff} keyword.
+ignored. The user may also specify the
+fugacity coefficient phi using the {fugacity_coeff} keyword, which
+defaults to unity.
 
 The {full_energy} option means that fix GCMC will compute the total
 potential energy of the entire simulated system. The total system
@@ -224,7 +257,8 @@ potential energy calculations, including the following:
   many-body pair styles
   hybrid pair styles
   eam pair styles
-  triclinic systems
+  tail corrections
+
   need to include potential energy contributions from other fixes :ul
 
 In these cases, LAMMPS will automatically apply the {full_energy}
@@ -276,6 +310,13 @@ therefore, you will want to use the
 current number of atoms is used as a normalizing factor each time
 temperature is computed.  Here is the necessary command:
 
+NOTE: If the density of the cell is initially very small or zero,
+and increases to a much larger density after a period of equilibration,
+then certain quantities that are only calculated once at the start
+(kspace parameters, tail corrections) may no longer be accurate.
+The solution is to start a new simulation after the equilibrium
+density has been reached.
+
 With some pair_styles, such as "Buckingham"_pair_buck.html,
 "Born-Mayer-Huggins"_pair_born.html and "ReaxFF"_pair_reax_c.html,
 two atoms placed close to each other may have an arbitrary large, 
@@ -366,7 +407,7 @@ referenced by the user for each subsequent fix gcmc command.
 [Default:]
 
 The option defaults are mol = no, maxangle = 10, overlap_cutoff = 0.0,
-and full_energy = no,
+fugacity_coeff = 1, and full_energy = no, 
 except for the situations where full_energy is required, as
 listed above.
 

doc/src/pair_charmm.txt

@@ -84,9 +84,9 @@ CHARMM force field.
 The styles with {charmm} (not {charmmfsw} or {charmmfsh}) in their
 name are the older, original LAMMPS implementations.  They compute the
 LJ and Coulombic interactions with an energy switching function (esw,
-a cubic polynomial, shown in the formula below), which ramps the
-energy smoothly to zero between the inner and outer cutoff.  This can
-cause irregularities in pair-wise forces (due to the discontinuous 2nd
+shown in the formula below as S(r)), which ramps the energy smoothly
+to zero between the inner and outer cutoff.  This can cause
+irregularities in pair-wise forces (due to the discontinuous 2nd
 derivative of energy at the boundaries of the switching region), which
 in some cases can result in detectable artifacts in an MD simulation.