Merge pull request #1673 from charlessievers/lammps_gjf

Lammps fix langevin updates for gjf option

doc/src/Errors_messages.txt

@@ -4791,6 +4791,22 @@ Self-explanatory. :dd
 
 This fix option cannot be used with point particles. :dd
 
+{Fix langevin gjf and respa are not compatible} :dt
+
+Self-explanatory. :dd
+
+{Fix langevin gjf cannot have period equal to dt/2} :dt
+
+If the period is equal to dt/2 then division by zero will happen. :dd
+
+{Fix langevin gjf should come before fix nve} :dt
+
+Self-explanatory. :dd
+
+{Fix langevin gjf with tbias is not yet implemented with kokkos} :dt
+
+This option is not yet available. :dd
+
 {Fix langevin omega is not yet implemented with kokkos} :dt
 
 This option is not yet available. :dd

doc/src/Errors_warnings.txt

@@ -248,6 +248,10 @@ included one or more of the following: kspace, triclinic, a hybrid
 pair style, an eam pair style, or no "single" function for the pair
 style. :dd
 
+{Fix langevin gjf using random gaussians is not implemented with kokkos} :dt
+
+This will most likely cause errors in kinetic fluctuations.
+
 {Fix property/atom mol or charge w/out ghost communication} :dt
 
 A model typically needs these properties defined for ghost atoms. :dd

doc/src/fix_langevin.txt

@@ -24,9 +24,10 @@ keyword = {angmom} or {omega} or {scale} or {tally} or {zero} :l
   {angmom} value = {no} or factor
     {no} = do not thermostat rotational degrees of freedom via the angular momentum
     factor = do thermostat rotational degrees of freedom via the angular momentum and apply numeric scale factor as discussed below
-  {gjf} value = {no} or {yes}
+  {gjf} value = {no} or {vfull} or {vhalf}
     {no} = use standard formulation
-    {yes} = use Gronbech-Jensen/Farago formulation
+    {vfull} = use Gronbech-Jensen/Farago formulation
+    {vhalf} = use 2GJ formulation
   {omega} value = {no} or {yes}
     {no} = do not thermostat rotational degrees of freedom via the angular velocity
     {yes} = do thermostat rotational degrees of freedom via the angular velocity
@@ -217,6 +218,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,29 +237,24 @@ 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
-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
-method should not be used to generate properties that depend on the
-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 +312,10 @@ 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:]
 
@@ -335,5 +338,10 @@ types, tally = no, zero = no, gjf = no.
 
 :link(Gronbech-Jensen)
 [(Gronbech-Jensen)] Gronbech-Jensen and Farago, Mol Phys, 111, 983
-(2013); Gronbech-Jensen, Hayre, and Farago, Comp Phys Comm,
-185, 524 (2014)
+(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