Merge remote-tracking branch 'github/develop' into cmake-adjustments

doc/src/Howto.rst

@@ -66,6 +66,7 @@ Force fields howto
    :name: force_howto
    :maxdepth: 1
 
+   Howto_FFgeneral
    Howto_bioFF
    Howto_amoeba
    Howto_tip3p

doc/src/Howto_FFgeneral.rst

@@ -0,0 +1,55 @@
+Some general force field considerations
+=======================================
+
+A compact summary of the concepts, definitions, and properties of force
+fields with explicit bonded interactions (like the ones discussed in
+this HowTo) is given in :ref:`(Gissinger) <Typelabel2>`.
+
+A force field has 2 parts: the formulas that define its potential
+functions and the coefficients used for a particular system.  To assign
+parameters it is first required to assign atom types.  Those are not
+only based on the elements, but also on the chemical environment due to
+the atoms bound to them.  This often follows the chemical concept of
+*functional groups*.  Example: a carbon atom bound with a single bond to
+a single OH-group (alcohol) would be a different atom type than a carbon
+atom bound to a methyl CH3 group (aliphatic carbon).  The atom types
+usually then determine the non-bonded Lennard-Jones parameters and the
+parameters for bonds, angles, dihedrals, and impropers.  On top of that,
+partial charges have to be applied.  Those are usually independent of
+the atom types and are determined either for groups of atoms called
+residues with some fitting procedure based on quantum mechanical
+calculations, or based on some increment system that add or subtract
+increments from the partial charge of an atom based on the types of
+the neighboring atoms.
+
+Force fields differ in the strategies they employ to determine the
+parameters and charge distribution in how generic or specific they are
+which in turn has an impact on the accuracy (compare for example
+CGenFF to CHARMM and GAFF to Amber).  Because of the different
+strategies, it is not a good idea to use a mix of parameters from
+different force field *families* (like CHARMM, Amber, or GROMOS)
+and that extends to the parameters for the solvent, especially
+water.  The publication describing the parameterization of a force
+field will describe which water model to use.  Changing the water
+model usually leads to overall worse results (even if it may improve
+on the water itself).
+
+In addition, one has to consider that *families* of force fields like
+CHARMM, Amber, OPLS, or GROMOS have evolved over time and thus provide
+different *revisions* of the force field parameters.  These often
+corresponds to changes in the functional form or the parameterization
+strategies.  This may also result in changes required for simulation
+settings like the preferred cutoff or how Coulomb interactions are
+computed (cutoff, smoothed/shifted cutoff, or long-range with Ewald
+summation or equivalent).  Unless explicitly stated in the publication
+describing the force field, the Coulomb interaction cannot be chosen at
+will but must match the revision of the force field.  That said,
+liberties may be taken during the initial equilibration of a system to
+speed up the process, but not for production simulations.
+
+----------
+
+.. _Typelabel2:
+
+**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).
+

doc/src/Howto_bioFF.rst

@@ -1,22 +1,16 @@
 CHARMM, AMBER, COMPASS, DREIDING, and OPLS force fields
 =======================================================
 
-A compact summary of the concepts, definitions, and properties of
+Here we only discuss formulas implemented in LAMMPS that correspond to
-force fields with explicit bonded interactions (like the ones discussed
+formulas commonly used in the CHARMM, AMBER, COMPASS, and DREIDING force
-in this HowTo) is given in :ref:`(Gissinger) <Typelabel2>`.
+fields.  Setting coefficients is done either from special sections in an
-
+input data file via the :doc:`read_data <read_data>` command or in the
-A force field has 2 parts: the formulas that define it and the
+input script with commands like :doc:`pair_coeff <pair_coeff>` or
-coefficients used for a particular system.  Here we only discuss
+:doc:`bond_coeff <bond_coeff>` and so on.  See the :doc:`Tools <Tools>`
-formulas implemented in LAMMPS that correspond to formulas commonly used
+doc page for additional tools that can use CHARMM, AMBER, or Materials
-in the CHARMM, AMBER, COMPASS, and DREIDING force fields.  Setting
+Studio generated files to assign force field coefficients and convert
-coefficients is done either from special sections in an input data file
+their output into LAMMPS input. LAMMPS input scripts can also be
-via the :doc:`read_data <read_data>` command or in the input script with
+generated by `charmm-gui.org <https://charmm-gui.org/>`_.
-commands like :doc:`pair_coeff <pair_coeff>` or :doc:`bond_coeff
-<bond_coeff>` and so on.  See the :doc:`Tools <Tools>` doc page for
-additional tools that can use CHARMM, AMBER, or Materials Studio
-generated files to assign force field coefficients and convert their
-output into LAMMPS input. LAMMPS input scripts can also be generated by
-`charmm-gui.org <https://charmm-gui.org/>`_.
 
 CHARMM and AMBER
 ----------------
@@ -203,9 +197,11 @@ rather than individual force constants and geometric parameters that
 depend on the particular combinations of atoms involved in the bond,
 angle, or torsion terms.  DREIDING has an :doc:`explicit hydrogen bond
 term <pair_hbond_dreiding>` to describe interactions involving a
-hydrogen atom on very electronegative atoms (N, O, F).  Unlike CHARMM
+hydrogen atom on very electronegative atoms (N, O, F).  Unlike CHARMM or
-or AMBER, the DREIDING force field has not been parameterized for
+AMBER, the DREIDING force field has not been parameterized for
-considering solvents (like water).
+considering solvents (like water) and has no rules for assigning
+(partial) charges.  That will seriously limit its accuracy when used for
+simulating systems where those matter.
 
 See :ref:`(Mayo) <howto-Mayo>` for a description of the DREIDING force field
 
@@ -272,10 +268,6 @@ compatible with a subset of OPLS interactions.
 
 ----------
 
-.. _Typelabel2:
-
-**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024).
-
 .. _howto-MacKerell:
 
 **(MacKerell)** MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, et al (1998).  J Phys Chem, 102, 3586 . https://doi.org/10.1021/jp973084f

doc/src/Howto_spc.rst

@@ -1,5 +1,5 @@
-SPC water model
+SPC and SPC/E water model
-===============
+=========================
 
 The SPC water model specifies a 3-site rigid water molecule with
 charges and Lennard-Jones parameters assigned to each of the three atoms.

doc/src/Howto_tip4p.rst

@@ -1,5 +1,5 @@
-TIP4P water model
+TIP4P and OPC water models
-=================
+==========================
 
 The four-point TIP4P rigid water model extends the traditional
 :doc:`three-point TIP3P <Howto_tip3p>` model by adding an additional
@@ -9,9 +9,11 @@ the oxygen along the bisector of the HOH bond angle.  A bond style of
 :doc:`harmonic <bond_harmonic>` and an angle style of :doc:`harmonic
 <angle_harmonic>` or :doc:`charmm <angle_charmm>` should also be used.
 In case of rigid bonds also bond style :doc:`zero <bond_zero>` and angle
-style :doc:`zero <angle_zero>` can be used.
+style :doc:`zero <angle_zero>` can be used.  Very similar to the TIP4P
+model is the OPC water model.  It can be realized the same way as TIP4P
+but has different geometry and force field parameters.
 
-There are two ways to implement TIP4P water in LAMMPS:
+There are two ways to implement TIP4P-like water in LAMMPS:
 
 #. Use a specially written pair style that uses the :ref:`TIP3P geometry
    <tip3p_molecule>` without the point M. The point M location is then
@@ -21,7 +23,10 @@ There are two ways to implement TIP4P water in LAMMPS:
    computationally very efficient, but the charge distribution in space
    is only correct within the tip4p labeled styles.  So all other
    computations using charges will "see" the negative charge incorrectly
-   on the oxygen atom.
+   located on the oxygen atom unless they are specially written for using
+   the TIP4P geometry internally as well, e.g. :doc:`compute dipole/tip4p
+   <compute_dipole>`, :doc:`fix efield/tip4p <fix_efield>`, or
+   :doc:`kspace_style pppm/tip4p <kspace_style>`.
 
    This can be done with the following pair styles for Coulomb with a cutoff:
 
@@ -68,77 +73,90 @@ TIP4P/2005 model :ref:`(Abascal2) <Abascal2>` and a version of TIP4P
 parameters adjusted for use with a long-range Coulombic solver
 (e.g. Ewald or PPPM in LAMMPS).  Note that for implicit TIP4P models the
 OM distance is specified in the :doc:`pair_style <pair_style>` command,
-not as part of the pair coefficients.
+not as part of the pair coefficients. Also parameters for the OPC
+model (:ref:`Izadi <Izadi>`) are provided.
 
 .. list-table::
       :header-rows: 1
-      :widths: 36 19 13 15 17
+      :widths: 40 12 12 14 11 11
 
       * - Parameter
         - TIP4P (original)
         - TIP4P/Ice
         - TIP4P/2005
         - TIP4P (Ewald)
+        - OPC
       * - O mass (amu)
         - 15.9994
         - 15.9994
         - 15.9994
         - 15.9994
+        - 15.9994
       * - H mass (amu)
         - 1.008
         - 1.008
         - 1.008
         - 1.008
+        - 1.008
       * - O or M charge (:math:`e`)
         - -1.040
         - -1.1794
         - -1.1128
         - -1.04844
+        - -1.3582
       * - H charge (:math:`e`)
         - 0.520
         - 0.5897
         - 0.5564
         - 0.52422
+        - 0.6791
       * - LJ :math:`\epsilon` of OO (kcal/mole)
         - 0.1550
         - 0.21084
         - 0.1852
         - 0.16275
+        - 0.21280
       * - LJ :math:`\sigma` of OO (:math:`\AA`)
         - 3.1536
         - 3.1668
         - 3.1589
         - 3.16435
+        - 3.1660
       * - LJ :math:`\epsilon` of HH, MM, OH, OM, HM (kcal/mole)
         - 0.0
         - 0.0
         - 0.0
         - 0.0
+        - 0.0
       * - LJ :math:`\sigma` of HH, MM, OH, OM, HM (:math:`\AA`)
         - 1.0
         - 1.0
         - 1.0
         - 1.0
+        - 1.0
       * - :math:`r_0` of OH bond (:math:`\AA`)
         - 0.9572
         - 0.9572
         - 0.9572
         - 0.9572
+        - 0.8724
       * - :math:`\theta_0` of HOH angle
         - 104.52\ :math:`^{\circ}`
         - 104.52\ :math:`^{\circ}`
         - 104.52\ :math:`^{\circ}`
         - 104.52\ :math:`^{\circ}`
+        - 103.60\ :math:`^{\circ}`
       * - OM distance (:math:`\AA`)
         - 0.15
         - 0.1577
         - 0.1546
         - 0.1250
+        - 0.1594
 
-Note that the when using the TIP4P pair style, the neighbor list cutoff
+Note that the when using a TIP4P pair style, the neighbor list cutoff
 for Coulomb interactions is effectively extended by a distance 2 \* (OM
 distance), to account for the offset distance of the fictitious charges
-on O atoms in water molecules.  Thus it is typically best in an
+on O atoms in water molecules.  Thus, it is typically best in an
 efficiency sense to use a LJ cutoff >= Coulomb cutoff + 2\*(OM
 distance), to shrink the size of the neighbor list.  This leads to
 slightly larger cost for the long-range calculation, so you can test the
@@ -192,6 +210,94 @@ file changed):
     run 20000
     write_data tip4p-implicit.data nocoeff
 
+When constructing an OPC model, we cannot use the ``tip3p.mol`` file due
+to the different geometry.  Below is a molecule file providing the 3
+sites of an implicit OPC geometry for use with TIP4P styles.  Note, that
+the "Shake" and "Special" sections are missing here.  Those will be
+auto-generated by LAMMPS when the molecule file is loaded *after* the
+simulation box has been created.  These sections are required only when
+the molecule file is loaded *before*.
+
+.. _opc3p_molecule:
+.. code-block::
+
+   # Water molecule. 3 point geometry for OPC model
+
+   3 atoms
+   2 bonds
+   1 angles
+
+   Coords
+
+   1    0.00000  -0.06037   0.00000
+   2    0.68558   0.50250   0.00000
+   3   -0.68558   0.50250   0.00000
+
+   Types
+
+   1        1   # O
+   2        2   # H
+   3        2   # H
+
+   Charges
+
+   1       -1.3582
+   2        0.6791
+   3        0.6791
+
+   Bonds
+
+   1   1      1      2
+   2   1      1      3
+
+   Angles
+
+   1   1      2      1      3
+
+Below is a LAMMPS input file using the implicit method to implement
+the OPC model using the molecule file from above and including the
+PPPM long-range Coulomb solver.
+
+.. code-block:: LAMMPS
+
+    units real
+    atom_style full
+    region box block -5 5 -5 5 -5 5
+    create_box 2 box bond/types 1 angle/types 1 &
+                extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
+
+    mass 1 15.9994
+    mass 2 1.008
+
+    pair_style lj/cut/tip4p/long 1 2 1 1 0.1594 12.0
+    pair_coeff 1 1 0.2128 3.166
+    pair_coeff 2 2 0.0    1.0
+
+    bond_style zero
+    bond_coeff 1 0.8724
+
+    angle_style zero
+    angle_coeff 1 103.6
+
+    kspace_style pppm/tip4p 1.0e-5
+
+    molecule water opc3p.mol  # this file has the OPC geometry but is without M
+    create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33
+
+    fix rigid all shake 0.001 10 10000 b 1 a 1
+    minimize 0.0 0.0 1000 10000
+
+    reset_timestep 0
+    timestep 1.0
+    velocity all create 300.0 5463576
+    fix integrate all nvt temp 300 300 100.0
+
+    thermo_style custom step temp press etotal pe
+
+    thermo 1000
+    run 20000
+    write_data opc-implicit.data nocoeff
+
 Below is the code for a LAMMPS input file using the explicit method and
 a TIP4P molecule file.  Because of using :doc:`fix rigid/small
 <fix_rigid>` no bonds need to be defined and thus no extra storage needs
@@ -279,3 +385,8 @@ Phys, 79, 926 (1983).
 
 **(Abascal2)** Abascal, J Chem Phys, 123, 234505 (2005)
    https://doi.org/10.1063/1.2121687
+
+.. _Izadi:
+
+**(Izadi)** Izadi, Anandakrishnan, Onufriev, J. Phys. Chem. Lett., 5, 21, 3863 (2014)
+   https://doi.org/10.1021/jz501780a

doc/src/Intro_authors.rst

@@ -84,8 +84,9 @@ lammps.org".  General questions about LAMMPS should be posted in the
 
    \normalsize
 
-Past developers include Paul Crozier and Mark Stevens, both at SNL,
+Past core developers include Paul Crozier and Mark Stevens, both at SNL,
-and Ray Shan, now at Materials Design.
+and Ray Shan while at SNL and later at Materials Design, now at Thermo
+Fisher Scientific.
 
 ----------
 

doc/src/Intro_portability.rst

@@ -28,8 +28,9 @@ Build systems
 LAMMPS can be compiled from source code using a (traditional) build
 system based on shell scripts, a few shell utilities (grep, sed, cat,
 tr) and the GNU make program. This requires running within a Bourne
-shell (``/bin/sh``).  Alternatively, a build system with different back
+shell (``/bin/sh`` or ``/bin/bash``).  Alternatively, a build system
-ends can be created using CMake.  CMake must be at least version 3.16.
+with different back ends can be created using CMake.  CMake must be
+at least version 3.16.
 
 Operating systems
 ^^^^^^^^^^^^^^^^^
@@ -40,11 +41,18 @@ Also, compilation and correct execution on macOS and Windows (using
 Microsoft Visual C++) is checked automatically for the largest part of
 the source code.  Some (optional) features are not compatible with all
 operating systems, either through limitations of the corresponding
-LAMMPS source code or through incompatibilities of source code or
+LAMMPS source code or through incompatibilities or build system
-build system of required external libraries or packages.
+limitations of required external libraries or packages.
 
-Executables for Windows may be created natively using either Cygwin or
+Executables for Windows may be created either natively using Cygwin,
-Visual Studio or with a Linux to Windows MinGW cross-compiler.
+MinGW, Intel, Clang, or Microsoft Visual C++ compilers, or with a Linux
+to Windows MinGW cross-compiler.  Native compilation is supported using
+Microsoft Visual Studio or a terminal window (using the CMake build
+system).
+
+Executables for macOS may be created either using Xcode or GNU compilers
+installed with Homebrew.  In the latter case, building of LAMMPS through
+Homebrew instead of a manual compile is also possible.
 
 Additionally, FreeBSD and Solaris have been tested successfully to
 run LAMMPS and produce results consistent with those on Linux.
@@ -61,8 +69,9 @@ CPU architectures
 ^^^^^^^^^^^^^^^^^
 
 The primary CPU architecture for running LAMMPS is 64-bit x86, but also
-32-bit x86, and 64-bit ARM and PowerPC (64-bit, Little Endian) are
+64-bit ARM and PowerPC (64-bit, Little Endian) are currently regularly
-regularly tested.
+tested.  Further architectures are tested by Linux distributions that
+bundle LAMMPS.
 
 Portability compliance
 ^^^^^^^^^^^^^^^^^^^^^^

doc/src/compute_bond_local.rst

@@ -64,20 +64,32 @@ All these properties are computed for the pair of atoms in a bond,
 whether the two atoms represent a simple diatomic molecule, or are part
 of some larger molecule.
 
-The value *dist* is the current length of the bond.
+.. versionchanged:: TBD
-The values *dx*, *dy*, and *dz* are the xyz components of the
+
-*distance* between the pair of atoms. This value is always the
+   The sign of *dx*, *dy*, *dz* is no longer determined by the atom IDs
-distance from the atom of lower to the one with the higher id.
+   of the bonded atoms but by their order in the bond list to be
+   consistent with *fx*, *fy*, and *fz*.
+
+The value *dist* is the current length of the bond.  The values *dx*,
+*dy*, and *dz* are the :math:`(x,y,z)` components of the distance vector
+:math:`\vec{x_i} - \vec{x_j}` between the atoms in the bond.  The order
+of the atoms is determined by the bond list and the respective atom-IDs
+can be output with :doc:`compute property/local
+<compute_property_local>`.
 
 The value *engpot* is the potential energy for the bond,
 based on the current separation of the pair of atoms in the bond.
 
-The value *force* is the magnitude of the force acting between the
+The value *force* is the magnitude of the force acting between the pair
-pair of atoms in the bond.
+of atoms in the bond, which is positive for a repulsive force and
+negative for an attractive force.
 
-The values *fx*, *fy*, and *fz* are the xyz components of
+The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of
-*force* between the pair of atoms in the bond. For bond styles that apply
+the force on the first atom *i* in the bond due to the second atom *j*.
-non-central forces, such as :doc:`bond_style bpm/rotational
+Mathematically, they are obtained by multiplying the value of *force*
+from above with a unit vector created from the *dx*, *dy*, and *dz*
+components of the distance vector also described above.  For bond styles
+that apply non-central forces, such as :doc:`bond_style bpm/rotational
 <bond_bpm_rotational>`, these values only include the :math:`(x,y,z)`
 components of the normal force component.
 

doc/src/compute_pair_local.rst

@@ -56,19 +56,33 @@ force cutoff distance for that interaction, as defined by the
 :doc:`pair_style <pair_style>` and :doc:`pair_coeff <pair_coeff>`
 commands.
 
-The value *dist* is the distance between the pair of atoms.
+.. versionchanged:: TBD
-The values *dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the
+
-*distance* between the pair of atoms. This value is always the
+   The sign of *dx*, *dy*, *dz* is no longer determined by the value of
-distance from the atom of higher to the one with the lower atom ID.
+   their atom-IDs but by their order in the neighbor list to be
+   consistent with *fx*, *fy*, and *fz*.
+
+The value *dist* is the distance between the pair of atoms.  The values
+*dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the distance
+vector :math:`\vec{x_i} - \vec{x_j}` between the pair of atoms.  The
+order of the atoms is determined by the neighbor list and the respective
+atom-IDs can be output with :doc:`compute property/local
+<compute_property_local>`.
 
 The value *eng* is the interaction energy for the pair of atoms.
 
 The value *force* is the force acting between the pair of atoms, which
 is positive for a repulsive force and negative for an attractive
-force.  The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of
+force.
-*force* on atom I. For pair styles that apply non-central forces,
+
-such as :doc:`granular pair styles <pair_gran>`, these values only include
+The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of
-the :math:`(x,y,z)` components of the normal force component.
+the force vector on the first atom *i* of a pair in the neighbor list
+due to the second atom *j*.  Mathematically, they are obtained by
+multiplying the value of *force* from above with a unit vector created
+from the *dx*, *dy*, and *dz* components of the distance vector also
+described above.  For pair styles that apply non-central forces, such as
+:doc:`granular pair styles <pair_gran>`, these values only include the
+:math:`(x,y,z)` components of the normal force component.
 
 A pair style may define additional pairwise quantities which can be
 accessed as *p1* to *pN*, where :math:`N` is defined by the pair style.

doc/utils/sphinx-config/false_positives.txt

@@ -103,6 +103,7 @@ Amit
 amsmath
 amu
 Amzallag
+Anandakrishnan
 analytical
 Anders
 Andric
@@ -397,6 +398,7 @@ Broglie
 brownian
 brownw
 Broyden
+Bruenger
 Bruskin
 Brusselle
 Bryantsev
@@ -1731,6 +1733,7 @@ Iyz
 iz
 izcm
 ized
+Izadi
 Izrailev
 Izumi
 Izvekov
@@ -2808,6 +2811,7 @@ oneMKL
 oneway
 onlysalt
 ons
+Onufriev
 OO
 Oord
 opencl

examples/PACKAGES/imd/in.bucky-plus-cnt

@@ -46,8 +46,8 @@ fix     integrate   mobile nve
 fix     thermostat  mobile langevin 300.0 300.0 2000.0 234624
 
 # IMD setup.
-fix  comm       all imd 6789 unwrap on trate 10
+#fix  comm       all imd 6789 unwrap on trate 10
-#fix  comm       all imd 6789 unwrap on trate 10 nowait on
+fix  comm       all imd 6789 unwrap on trate 10 nowait on
 
 # temperature is based on mobile atoms only
 compute mobtemp mobile temp

examples/PACKAGES/imd/in.bucky-plus-cnt-gpu

@@ -1,16 +1,20 @@
 # stick a buckyball into a nanotube
+
+# enable GPU package from within the input:
+package gpu 0 pair/only on
+suffix gpu
+
 units           real
 dimension       3
 boundary       f f f
 atom_style      molecular
-newton          off
 
 processors * * 1
 
 # read topology 
 read_data       data.bucky-plus-cnt
 
-pair_style  lj/cut/gpu  10.0
+pair_style  lj/cut  10.0
 bond_style  harmonic
 angle_style charmm
 dihedral_style charmm
@@ -29,9 +33,6 @@ neigh_modify    delay 0 every 1 check yes
 
 timestep        2.0
 
-# required for GPU acceleration
-fix   gpu  all      gpu  force 0 0 1.0
-
 # we only move some atoms.
 group mobile type 1
 
@@ -49,8 +50,8 @@ fix     integrate   mobile nve
 fix     thermostat  mobile langevin 300.0 300.0 2000.0 234624
 
 # IMD setup.
-fix  comm       all imd 6789 unwrap on trate 10
+#fix  comm       all imd 6789 unwrap on trate 10
-#fix  comm       all imd 6789 unwrap on trate 10 nowait on
+fix  comm       all imd 6789 unwrap on trate 10 nowait on
 
 # temperature is based on mobile atoms only
 compute mobtemp mobile temp

examples/PACKAGES/imd/in.deca-ala_imd-gpu

@@ -1,8 +1,12 @@
 #
+
+# enable GPU package from within the input:
+package gpu 0 pair/only on
+suffix gpu
+
 units           real
 neighbor        2.5 bin
 neigh_modify    delay 1 every 1 
-newton          off
 
 atom_style      full
 bond_style      harmonic
@@ -10,20 +14,18 @@ angle_style     charmm
 dihedral_style  charmm
 improper_style  harmonic
 
-pair_style      lj/charmm/coul/long/gpu 8 10
+pair_style      lj/charmm/coul/long 8 10
 pair_modify     mix arithmetic
 special_bonds   charmm
 read_data       data.deca-ala-solv
 
-fix             0 all gpu force/neigh 0 0 1.0    
-
 group peptide   id <= 103
 fix             rigidh all shake 1e-6 100 1000 t 1 2 3 4 5 a 23
 
 thermo          100
 thermo_style    multi
 timestep        2.0
-kspace_style    pppm/gpu 1e-5
+kspace_style    pppm 1e-5
 
 fix             ensemble all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 0.2
 

examples/PACKAGES/imd/in.melt_imd-gpu

@@ -1,8 +1,11 @@
-# 3d Lennard-Jones melt
+# 3d Lennard-Jones melt with GPU package acceleration
+
+# enable GPU package from within the input:
+package gpu 0
+suffix gpu
 
 units           lj
 atom_style      atomic
-newton      off
 
 lattice         fcc 0.8442
 region          box block 0 10 0 10 0 10
@@ -12,7 +15,7 @@ mass		1 1.0
 
 velocity        all create 3.0 87287
 
-pair_style	lj/cut/gpu 2.5
+pair_style      lj/cut 2.5
 pair_coeff      1 1 1.0 1.0 2.5
 
 neighbor        0.3 bin
@@ -20,7 +23,6 @@ neigh_modify	every 5 delay 10 check yes
 
 thermo_style    custom step pe ke spcpu
 
-fix     0 all gpu force/neigh 0 0 1.0
 fix             1 all nve
 
 # IMD setup.

src/KOKKOS/meam_dens_init_kokkos.h

@@ -236,7 +236,6 @@ MEAMKokkos<DeviceType>::meam_dens_init(int inum_half, int ntype, typename AT::t_
   this->d_neighbors_half = d_neighbors_half;
   this->d_neighbors_full = d_neighbors_full;
   this->d_offset = d_offset;
-  this->nlocal = nlocal;
 
   if (need_dup) {
     dup_rho0 = Kokkos::Experimental::create_scatter_view<Kokkos::Experimental::ScatterSum, Kokkos::Experimental::ScatterDuplicated>(d_rho0);

src/KOKKOS/meam_force_kokkos.h

@@ -17,7 +17,6 @@ void MEAMKokkos<DeviceType>::meam_force(
 {
   EV_FLOAT ev;
 
-  this->eflag_either = eflag_either;
   this->eflag_global = eflag_global;
   this->eflag_atom = eflag_atom;
   this->vflag_global = vflag_global;

src/KOKKOS/pair_kokkos.h

@@ -961,7 +961,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, std::enable_if_t<(NEIGHFLAG&P
       lastcall = fpair->lmp->update->ntimestep;
       vectorsize = GetMaxNeighs(list);
       if (vectorsize == 0) vectorsize = 1;
-      vectorsize = MathSpecial::powint(2,(int(log2(vectorsize) + 0.5))); // round to nearest power of 2
+      vectorsize = MathSpecial::powint(2.0,(int(log2(double(vectorsize)) + 0.5))); // round to nearest power of 2
 
   #if defined(KOKKOS_ENABLE_HIP)
       int max_vectorsize = 64;

src/LATBOLTZ/fix_lb_fluid.cpp

@@ -2344,7 +2344,7 @@ void FixLbFluid::SetupBuffers()
                   MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &dump_file_handle_raw);
 
     MPI_File_set_size(dump_file_handle_raw, 0);
-    MPI_File_set_view(dump_file_handle_raw, 0, MPI_DOUBLE, dump_file_mpitype, "native",
+    MPI_File_set_view(dump_file_handle_raw, 0, MPI_DOUBLE, dump_file_mpitype, (char *)"native",
                       MPI_INFO_NULL);
   }
 }

src/ML-POD/compute_podd_atom.cpp

@@ -58,8 +58,8 @@ ComputePODDAtom::ComputePODDAtom(LAMMPS *lmp, int narg, char **arg) :
   pod = nullptr;
   elements = nullptr;
 
-  if (((((MAXBIGINT*3.0)*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0))
+  if ((((3.0*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0))
-    error->all(FLERR, "Per-atom data too large");
+      error->all(FLERR, "Too many atoms ({}) for compute {}", atom->natoms, style);
   size_peratom_cols = 3 * atom->natoms * podptr->Mdesc * podptr->nClusters;
   peratom_flag = 1;
 }
@@ -110,8 +110,8 @@ void ComputePODDAtom::compute_peratom()
   if (atom->natoms > nmax) {
     memory->destroy(pod);
     nmax = atom->natoms;
-    if (((((MAXBIGINT*3.0)*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0))
+    if ((((3.0*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0))
-      error->all(FLERR, "Per-atom data too large");
+      error->all(FLERR, "Too many atoms ({}) for compute {}", atom->natoms, style);
     int numdesc = 3 * atom->natoms * podptr->Mdesc * podptr->nClusters;
     memory->create(pod, nmax, numdesc,"podd/atom:pod");
     array_atom = pod;

src/compute_bond_local.cpp

@@ -428,22 +428,19 @@ int ComputeBondLocal::compute_bonds(int flag)
           if (dstr) input->variable->internal_set(dvar, sqrt(rsq));
         }
 
-        // to make sure dx, dy and dz are always from the lower to the higher id
-        double directionCorrection = tag[atom1] > tag[atom2] ? -1.0 : 1.0;
-
         for (int n = 0; n < nvalues; n++) {
           switch (bstyle[n]) {
             case DIST:
               ptr[n] = sqrt(rsq);
               break;
             case DX:
-              ptr[n] = dx * directionCorrection;
+              ptr[n] = dx;
               break;
             case DY:
-              ptr[n] = dy * directionCorrection;
+              ptr[n] = dy;
               break;
             case DZ:
-              ptr[n] = dz * directionCorrection;
+              ptr[n] = dz;
               break;
             case ENGPOT:
               ptr[n] = engpot;

src/compute_pair_local.cpp

@@ -277,22 +277,19 @@ int ComputePairLocal::compute_pairs(int flag)
         else
           ptr = alocal[m];
 
-        // to make sure dx, dy and dz are always from the lower to the higher id
-        double directionCorrection = itag > jtag ? -1.0 : 1.0;
-
         for (n = 0; n < nvalues; n++) {
           switch (pstyle[n]) {
             case DIST:
               ptr[n] = sqrt(rsq);
               break;
             case DX:
-              ptr[n] = delx * directionCorrection;
+              ptr[n] = delx;
               break;
             case DY:
-              ptr[n] = dely * directionCorrection;
+              ptr[n] = dely;
               break;
             case DZ:
-              ptr[n] = delz * directionCorrection;
+              ptr[n] = delz;
               break;
             case ENG:
               ptr[n] = eng;

src/library.cpp

@@ -1622,6 +1622,11 @@ int lammps_extract_global_datatype(void * /*handle*/, const char *name)
   if (strcmp(name,"neigh_dihedrallist") == 0) return LAMMPS_INT_2D;
   if (strcmp(name,"neigh_improperlist") == 0) return LAMMPS_INT_2D;
 
+  if (strcmp(name,"eflag_global") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"eflag_atom") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"vflag_global") == 0) return LAMMPS_BIGINT;
+  if (strcmp(name,"vflag_atom") == 0) return LAMMPS_BIGINT;
+
   if (strcmp(name,"map_style") == 0) return LAMMPS_INT;
 #if defined(LAMMPS_BIGBIG)
   if (strcmp(name,"map_tag_max") == 0) return LAMMPS_BIGINT;
@@ -1707,6 +1712,7 @@ report the "native" data type.  The following tables are provided:
 * :ref:`Simulation box settings <extract_box_settings>`
 * :ref:`System property settings <extract_system_settings>`
 * :ref:`Neighbor topology data <extract_neighbor_lists>`
+* :ref:`Energy and virial tally settings <extract_tally_settings>`
 * :ref:`Git revision and version settings <extract_git_settings>`
 * :ref:`Unit settings <extract_unit_settings>`
 
@@ -1984,6 +1990,35 @@ Get length of lists with :ref:`lammps_extract_setting() <extract_neighbor_settin
      - nimproperlist
      - list of impropers (atom1, atom2, atom3, atom4, type)
 
+.. _extract_tally_settings:
+
+**Energy and virial tally settings**
+
+.. list-table::
+   :header-rows: 1
+   :widths: 20 12 16 52
+
+   * - Name
+     - Type
+     - Length
+     - Description
+   * - eflag_global
+     - bigint
+     - 1
+     - timestep global energy is tallied on
+   * - eflag_atom
+     - bigint
+     - 1
+     - timestep per-atom energy is tallied on
+   * - vflag_global
+     - bigint
+     - 1
+     - timestep global virial is tallied on
+   * - vflag_atom
+     - bigint
+     - 1
+     - timestep per-atom virial is tallied on
+
 .. _extract_git_settings:
 
 **Git revision and version settings**
@@ -2194,10 +2229,15 @@ void *lammps_extract_global(void *handle, const char *name)
 
   if (strcmp(name,"q_flag") == 0) return (void *) &lmp->atom->q_flag;
 
-  if (strcmp(name,"neigh_bondlist") == 0) return lmp->neighbor->bondlist;
+  if (strcmp(name,"neigh_bondlist") == 0) return (void *) lmp->neighbor->bondlist;
-  if (strcmp(name,"neigh_anglelist") == 0) return lmp->neighbor->anglelist;
+  if (strcmp(name,"neigh_anglelist") == 0) return (void *) lmp->neighbor->anglelist;
-  if (strcmp(name,"neigh_dihedrallist") == 0) return lmp->neighbor->dihedrallist;
+  if (strcmp(name,"neigh_dihedrallist") == 0) return (void *) lmp->neighbor->dihedrallist;
-  if (strcmp(name,"neigh_improperlist") == 0) return lmp->neighbor->improperlist;
+  if (strcmp(name,"neigh_improperlist") == 0) return (void *) lmp->neighbor->improperlist;
+
+  if (strcmp(name,"eflag_global") == 0) return (void *) &lmp->update->eflag_global;
+  if (strcmp(name,"eflag_atom") == 0) return (void *) &lmp->update->eflag_atom;
+  if (strcmp(name,"vflag_global") == 0) return (void *) &lmp->update->vflag_global;
+  if (strcmp(name,"vflag_atom") == 0) return (void *) &lmp->update->vflag_atom;
 
   if (strcmp(name,"map_style") == 0) return (void *) &lmp->atom->map_style;
   if (strcmp(name,"map_tag_max") == 0) return (void *) &lmp->atom->map_tag_max;

src/potential_file_reader.cpp

@@ -86,11 +86,21 @@ PotentialFileReader::~PotentialFileReader()
 /** Set comment (= text after '#') handling preference for the file to be read
  *
  * \param   value   Comment text is ignored if true, or not if false */
+
 void PotentialFileReader::ignore_comments(bool value)
 {
   reader->ignore_comments = value;
 }
 
+/** Set line buffer size of the internal TextFileReader class instance.
+ *
+ * \param   bufsize   New size of the line buffer */
+
+void PotentialFileReader::set_bufsize(int bufsize)
+{
+  reader->set_bufsize(bufsize);
+}
+
 /** Reset file to the beginning */
 
 void PotentialFileReader::rewind()

src/potential_file_reader.h

@@ -41,6 +41,7 @@ class PotentialFileReader : protected Pointers {
                       const int auto_convert = 0);
   ~PotentialFileReader() override;
 
+  void set_bufsize(int bufsize);
   void ignore_comments(bool value);
 
   void rewind();

src/read_data.cpp

@@ -174,13 +174,13 @@ void ReadData::command(int narg, char **arg)
         addflag = VALUE;
         bigint offset = utils::bnumeric(FLERR, arg[iarg + 1], false, lmp);
         if (offset > MAXTAGINT)
-          error->all(FLERR, "Read data add atomID offset {} is too big", offset);
+          error->all(FLERR, "Read data add IDoffset {} is too big", offset);
         id_offset = offset;
 
         if (atom->molecule_flag) {
           offset = utils::bnumeric(FLERR, arg[iarg + 2], false, lmp);
           if (offset > MAXTAGINT)
-            error->all(FLERR, "Read data add molID offset {} is too big", offset);
+            error->all(FLERR, "Read data add MOLoffset {} is too big", offset);
           mol_offset = offset;
           iarg++;
         }

src/reset_atoms_mol.cpp

@@ -211,7 +211,7 @@ void ResetAtomsMol::reset()
 
     // if offset < 0 (default), reset it
     // if group = all, offset = 0
-    // else offset = largest molID of non-group atoms
+    // else offset = largest molecule ID of non-group atoms
 
     if (offset < 0) {
       if (groupbit != 1) {

src/text_file_reader.cpp

@@ -93,7 +93,9 @@ TextFileReader::~TextFileReader()
   delete[] line;
 }
 
-/** adjust line buffer size */
+/** adjust line buffer size
+ *
+ * \param  newsize  New size of the internal line buffer */
 
 void TextFileReader::set_bufsize(int newsize)
 {

unittest/c-library/test_library_properties.cpp

@@ -383,10 +383,48 @@ TEST_F(LibraryProperties, global)
     char *c_ptr = (char *)lammps_extract_global(lmp, "units");
     EXPECT_THAT(c_ptr, StrEq("real"));
 
-    EXPECT_EQ(lammps_extract_global_datatype(lmp, "ntimestep"), LAMMPS_INT64);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "ntimestep"), LAMMPS_BIGINT);
-    auto *b_ptr = (int64_t *)lammps_extract_global(lmp, "ntimestep");
+    auto *b_ptr = (bigint *)lammps_extract_global(lmp, "ntimestep");
     EXPECT_EQ((*b_ptr), 2);
 
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "natoms"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "natoms");
+    EXPECT_EQ((*b_ptr), 29);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "nbonds"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "nbonds");
+    EXPECT_EQ((*b_ptr), 24);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "nangles"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "nangles");
+    EXPECT_EQ((*b_ptr), 30);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "ndihedrals"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "ndihedrals");
+    EXPECT_EQ((*b_ptr), 31);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "nimpropers"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "nimpropers");
+    EXPECT_EQ((*b_ptr), 2);
+
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "neigh_bondlist"), LAMMPS_INT_2D);
+    EXPECT_NE(lammps_extract_global(lmp, "neigh_bondlist"), nullptr);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "neigh_anglelist"), LAMMPS_INT_2D);
+    EXPECT_NE(lammps_extract_global(lmp, "neigh_anglelist"), nullptr);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "neigh_dihedrallist"), LAMMPS_INT_2D);
+    EXPECT_NE(lammps_extract_global(lmp, "neigh_dihedrallist"), nullptr);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "neigh_improperlist"), LAMMPS_INT_2D);
+    EXPECT_NE(lammps_extract_global(lmp, "neigh_improperlist"), nullptr);
+
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "eflag_global"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "eflag_global");
+    EXPECT_EQ((*b_ptr), 2);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "eflag_atom"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "eflag_atom");
+    EXPECT_EQ((*b_ptr), 0);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "vflag_global"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "vflag_global");
+    EXPECT_EQ((*b_ptr), 2);
+    EXPECT_EQ(lammps_extract_global_datatype(lmp, "vflag_atom"), LAMMPS_BIGINT);
+    b_ptr = (bigint *)lammps_extract_global(lmp, "vflag_atom");
+    EXPECT_EQ((*b_ptr), 0);
+
     EXPECT_EQ(lammps_extract_global_datatype(lmp, "dt"), LAMMPS_DOUBLE);
     auto *d_ptr = (double *)lammps_extract_global(lmp, "dt");
     EXPECT_DOUBLE_EQ((*d_ptr), 0.1);

unittest/python/python-numpy.py

@@ -153,7 +153,7 @@ class PythonNumpy(unittest.TestCase):
         self.assertEqual(values[0,0], 0.5)
         self.assertEqual(values[0,3], -0.5)
         self.assertEqual(values[1,0], 1.5)
-        self.assertEqual(values[1,3], 1.5)
+        self.assertEqual(values[1,3], -1.5)
 
     def testExtractAtom(self):
         self.lmp.command("units lj")