Updated documentation and added MD results

doc/src/fix_langevin.txt

@@ -217,6 +217,10 @@ the particles.  As described below, this energy can then be printed
 out or added to the potential energy of the system to monitor energy
 conservation.
 
+NOTE: this accumulated energy does NOT include kinetic energy removed
+by the {zero} flag. LAMMPS will print a warning when both options are
+active.
+
 The keyword {zero} can be used to eliminate drift due to the
 thermostat. Because the random forces on different atoms are
 independent, they do not sum exactly to zero.  As a result, this fix
@@ -232,21 +236,34 @@ The keyword {gjf} can be used to run the "Gronbech-Jensen/Farago
 described in the papers cited below, the purpose of this method is to
 enable longer timesteps to be used (up to the numerical stability
 limit of the integrator), while still producing the correct Boltzmann
-distribution of atom positions.  It is implemented within LAMMPS, by
-changing how the random force is applied so that it is composed of
-the average of two random forces representing half-contributions from
-the previous and current time intervals.
+distribution of atom positions.
 
-In common with all methods based on Verlet integration, the
-discretized velocities generated by this method in conjunction with
-velocity-Verlet time integration are not exactly conjugate to the
-positions.  As a result the temperature (computed from the discretized
-velocities) will be systematically lower than the target temperature,
-by a small amount which grows with the timestep.  Nonetheless, the
-distribution of atom positions will still be consistent with the
+The current implementation provides the user with the option to output
+the velocity in one of two forms: {vfull} or {vhalf}, which replaces
+the outdated option {yes}. The {gjf} option {vfull} outputs the on-site
+velocity given in "Gronbech-Jensen/Farago"_#Gronbech-Jensen; this velocity
+is shown to be systematically lower than the target temperature by a small
+amount, which grows quadratically with the timestep.
+The {gjf} option {vhalf} outputs the 2GJ half-step velocity given in
+"Gronbech Jensen/Gronbech-Jensen"_#2Gronbech-Jensen; this velocity is shown
+to not have any linear statistical errors for any stable time step.
+An overview of statistically correct Boltzmann and Maxwell-Boltzmann
+sampling of true on-site and true half-step velocities is given in
+"Gronbech-Jensen_#1Gronbech-Jensen.
+Regardless of the choice of output velocity, the sampling of the configurational
+distribution of atom positions is the same, and linearly consistent with the
 target temperature.
 
-As an example of using the {gjf} keyword, for molecules containing C-H
+An example of a reason why to use the {gjf} keyword is the freedom to take a larger time step,
+up to the stability limit, while maintaining robust statistics. It is crucial to
+recall that while the equilibrium statistics is appropriately sampled, the correct dynamics
+of the trajectories may not be for large time steps, as is the case for all thermostats.
+All thermostats provide good statistics and dynamics for small time steps.
+The 2GJ half-step velocity {vhalf} samples the correct velocity distribution for the {gjf} trajectory.
+Results of simulations using the {gjf} option with both {vfull} and {vhalf} compared to
+other available thermostats are shown in the LAMMPS directory: examples/gjf.
+
+As an example of why to use the {gjf} keyword, for molecules containing C-H
 bonds, configurational properties generated with dt = 2.5 fs and tdamp
 = 100 fs are indistinguishable from dt = 0.5 fs.  Because the velocity
 distribution systematically decreases with increasing timestep, the
@@ -255,6 +272,7 @@ velocity distribution, such as the velocity auto-correlation function
 (VACF). In this example, the velocity distribution at dt = 2.5fs
 generates an average temperature of 220 K, instead of 300 K.
 
+
 :line
 
 Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are
@@ -312,7 +330,8 @@ This fix can ramp its target temperature over multiple runs, using the
 
 This fix is not invoked during "energy minimization"_minimize.html.
 
-[Restrictions:] none
+[Restrictions:] For {gjf} do not choose damp=dt/2. {gjf} is not compatible
+with run_style respa. 
 
 [Related commands:]
 
@@ -337,3 +356,11 @@ types, tally = no, zero = no, gjf = no.
 [(Gronbech-Jensen)] Gronbech-Jensen and Farago, Mol Phys, 111, 983
 (2013); Gronbech-Jensen, Hayre, and Farago, Comp Phys Comm,
 185, 524 (2014)
+
+:link(2Gronbech-Jensen)
+[(Gronbech-Jensen)] Gronbech Jensen and Gronbech-Jensen, Mol Phys, 117, 2511
+(2019)
+
+:link(1Gronbech-Jensen)
+[(Gronbech-Jensen)] Gronbech-Jensen, Mol Phys (2019);
+https://doi.org/10.1080/00268976.2019.1662506