Merge branch 'master' of github.com:alxvov/lammps into OSO

- solved conflict in doc/src/fix_nve_spin.txt

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_bond_react.txt

@@ -266,7 +266,7 @@ either 'none' or 'charges.' Further details are provided in the
 discussion of the 'update_edges' keyword. The fourth optional section
 begins with the keyword 'Constraints' and lists additional criteria
 that must be satisfied in order for the reaction to occur. Currently,
-there is one type of constraint available, as discussed below.
+there are two types of constraints available, as discussed below.
 
 A sample map file is given below:
 
@@ -300,14 +300,23 @@ Equivalences :pre
 :line
 
 Any number of additional constraints may be specified in the
-Constraints section of the map file. Currently there is one type of
-additional constraint, of type 'distance', whose syntax is as follows:
+Constraints section of the map file. The constraint of type 'distance'
+has syntax as follows:
 
 distance {ID1} {ID2} {rmin} {rmax} :pre
 
 where 'distance' is the required keyword, {ID1} and {ID2} are
 pre-reaction atom IDs, and these two atoms must be separated by a
-distance between {rmin} and {rmax} for the reaction to occur. This
+distance between {rmin} and {rmax} for the reaction to occur.
+
+The constraint of type 'angle' has the following syntax:
+
+angle {ID1} {ID2} {ID3} {amin} {amax} :pre
+
+where 'angle' is the required keyword, {ID1}, {ID2} and {ID3} are
+pre-reaction atom IDs, and these three atoms must form an angle
+between {amin} and {amax} for the reaction to occur (where {ID2} is
+the central atom). Angles must be specified in degrees. This
 constraint can be used to enforce a certain orientation between
 reacting molecules.
 

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

doc/src/fix_langevin_spin.txt

@@ -50,7 +50,7 @@ As an example:
 
 fix 1 all precession/spin zeeman 0.01 0.0 0.0 1.0
 fix 2 all langevin/spin 300.0 0.01 21
-fix 3 all nve/spin lattice yes :pre
+fix 3 all nve/spin lattice moving :pre
 
 is correct, but defining a force/spin command after the langevin/spin command
 would give an error message.

doc/src/fix_nve_spin.txt

@@ -15,28 +15,26 @@ fix ID group-ID nve/spin keyword values :pre
 ID, group-ID are documented in "fix"_fix.html command :ulb,l
 nve/spin = style name of this fix command :l
 keyword = {lattice} :l
-  {lattice} value = {no} or {yes} :pre
+  {lattice} value = {moving} or {frozen}
+    moving = integrate both spin and atomic degress of freedom
+    frozen = integrate spins on a fixed lattice :pre
 :ule
 
 [Examples:]
 
-fix 3 all nve/spin lattice yes
-fix 1 all nve/spin lattice no :pre
+fix 3 all nve/spin lattice moving
+fix 1 all nve/spin lattice frozen :pre
 
 [Description:]
 
 Perform a symplectic integration for the spin or spin-lattice system.
 
-The {lattice} keyword defines whether the spins are integrated on a 
-fixed or moving lattice. 
-
-If {lattice}=yes, the equations of motion of the atoms are integrated, 
-and a combined spin and lattice calculation is performed.
-This is the default option.
-
-If {lattice}=no, the equations of motion of the atoms are not
-integrated. The lattice degrees of freedom are frozen, and a
-spin dynamics only calculation is performed.
+The {lattice} keyword defines if the spins are integrated on a lattice
+of fixed atoms (lattice = frozen), or if atoms are moving
+(lattice = moving).
+The first case corresponds to a spin dynamics calculation, and
+the second to a spin-lattice calculation.
+By default a spin-lattice integration is performed (lattice = moving).
 
 The {nve/spin} fix applies a Suzuki-Trotter decomposition to
 the equations of motion of the spin lattice system, following the scheme:
@@ -69,8 +67,9 @@ instead of "array" is also valid.
 
 "atom_style spin"_atom_style.html, "fix nve"_fix_nve.html
 
-[Default:] By default (lattice = yes), a spin-lattice integration is 
-performed.
+[Default:] 
+
+The option default is lattice = moving.
 
 :line