diff --git a/doc/src/Howto.rst b/doc/src/Howto.rst index ec90472c27..cdc4efd737 100644 --- a/doc/src/Howto.rst +++ b/doc/src/Howto.rst @@ -66,6 +66,7 @@ Force fields howto :name: force_howto :maxdepth: 1 + Howto_FFgeneral Howto_bioFF Howto_amoeba Howto_tip3p diff --git a/doc/src/Howto_FFgeneral.rst b/doc/src/Howto_FFgeneral.rst new file mode 100644 index 0000000000..1b96ae1119 --- /dev/null +++ b/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) `. + +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). + diff --git a/doc/src/Howto_bioFF.rst b/doc/src/Howto_bioFF.rst index 92dd45b9b6..cf8e4ab13e 100644 --- a/doc/src/Howto_bioFF.rst +++ b/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 -force fields with explicit bonded interactions (like the ones discussed -in this HowTo) is given in :ref:`(Gissinger) `. - -A force field has 2 parts: the formulas that define it and the -coefficients used for a particular system. Here we only discuss -formulas implemented in LAMMPS that correspond to formulas commonly used -in the CHARMM, AMBER, COMPASS, and DREIDING force fields. Setting -coefficients is done either from special sections in an input data file -via the :doc:`read_data ` command or in the input script with -commands like :doc:`pair_coeff ` or :doc:`bond_coeff -` and so on. See the :doc:`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 `_. +Here we only discuss formulas implemented in LAMMPS that correspond to +formulas commonly used in the CHARMM, AMBER, COMPASS, and DREIDING force +fields. Setting coefficients is done either from special sections in an +input data file via the :doc:`read_data ` command or in the +input script with commands like :doc:`pair_coeff ` or +:doc:`bond_coeff ` and so on. See the :doc:`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 `_. 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 ` to describe interactions involving a -hydrogen atom on very electronegative atoms (N, O, F). Unlike CHARMM -or AMBER, the DREIDING force field has not been parameterized for -considering solvents (like water). +hydrogen atom on very electronegative atoms (N, O, F). Unlike CHARMM or +AMBER, the DREIDING force field has not been parameterized for +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) ` 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 diff --git a/doc/src/Howto_spc.rst b/doc/src/Howto_spc.rst index 00bd8a1b10..f84d7797d2 100644 --- a/doc/src/Howto_spc.rst +++ b/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. diff --git a/doc/src/Howto_tip4p.rst b/doc/src/Howto_tip4p.rst index 47a1b9b578..76c470d615 100644 --- a/doc/src/Howto_tip4p.rst +++ b/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 ` 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 ` and an angle style of :doc:`harmonic ` or :doc:`charmm ` should also be used. In case of rigid bonds also bond style :doc:`zero ` and angle -style :doc:`zero ` can be used. +style :doc:`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 ` 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 + `, :doc:`fix efield/tip4p `, or + :doc:`kspace_style pppm/tip4p `. This can be done with the following pair styles for Coulomb with a cutoff: @@ -68,77 +73,90 @@ TIP4P/2005 model :ref:`(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 ` command, -not as part of the pair coefficients. +not as part of the pair coefficients. Also parameters for the OPC +model (:ref:`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 ` 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 diff --git a/doc/src/Intro_authors.rst b/doc/src/Intro_authors.rst index 38f1102595..730cd2e336 100644 --- a/doc/src/Intro_authors.rst +++ b/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, -and Ray Shan, now at Materials Design. +Past core developers include Paul Crozier and Mark Stevens, both at SNL, +and Ray Shan while at SNL and later at Materials Design, now at Thermo +Fisher Scientific. ---------- diff --git a/doc/src/Intro_portability.rst b/doc/src/Intro_portability.rst index 63ae147b8c..564cdc47f4 100644 --- a/doc/src/Intro_portability.rst +++ b/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 -ends can be created using CMake. CMake must be at least version 3.16. +shell (``/bin/sh`` or ``/bin/bash``). Alternatively, a build system +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 -build system of required external libraries or packages. +LAMMPS source code or through incompatibilities or build system +limitations of required external libraries or packages. -Executables for Windows may be created natively using either Cygwin or -Visual Studio or with a Linux to Windows MinGW cross-compiler. +Executables for Windows may be created either natively using Cygwin, +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 -regularly tested. +64-bit ARM and PowerPC (64-bit, Little Endian) are currently regularly +tested. Further architectures are tested by Linux distributions that +bundle LAMMPS. Portability compliance ^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/src/compute_bond_local.rst b/doc/src/compute_bond_local.rst index e070d507b1..7ce6f9b15a 100644 --- a/doc/src/compute_bond_local.rst +++ b/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. -The values *dx*, *dy*, and *dz* are the xyz components of the -*distance* between the pair of atoms. This value is always the -distance from the atom of lower to the one with the higher id. +.. versionchanged:: TBD + + The sign of *dx*, *dy*, *dz* is no longer determined by the atom IDs + 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 +`. 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 -pair of atoms in the bond. +The value *force* is the magnitude of the force acting between the pair +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 -*force* between the pair of atoms in the bond. For bond styles that apply -non-central forces, such as :doc:`bond_style bpm/rotational +The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of +the force on the first atom *i* in the bond 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 bond styles +that apply non-central forces, such as :doc:`bond_style bpm/rotational `, these values only include the :math:`(x,y,z)` components of the normal force component. diff --git a/doc/src/compute_pair_local.rst b/doc/src/compute_pair_local.rst index 31209f63f4..605bfc8e9e 100644 --- a/doc/src/compute_pair_local.rst +++ b/doc/src/compute_pair_local.rst @@ -56,19 +56,33 @@ force cutoff distance for that interaction, as defined by the :doc:`pair_style ` and :doc:`pair_coeff ` commands. -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* between the pair of atoms. This value is always the -distance from the atom of higher to the one with the lower atom ID. +.. versionchanged:: TBD + + The sign of *dx*, *dy*, *dz* is no longer determined by the value of + 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 +`. 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* on atom I. For pair styles that apply non-central forces, -such as :doc:`granular pair styles `, these values only include -the :math:`(x,y,z)` components of the normal force component. +force. + +The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of +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 `, 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. diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 40582263e2..af91745034 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/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 diff --git a/examples/PACKAGES/imd/in.bucky-plus-cnt b/examples/PACKAGES/imd/in.bucky-plus-cnt index b3eeff3cc1..af511fe11f 100644 --- a/examples/PACKAGES/imd/in.bucky-plus-cnt +++ b/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 nowait on +#fix comm all imd 6789 unwrap on trate 10 +fix comm all imd 6789 unwrap on trate 10 nowait on # temperature is based on mobile atoms only compute mobtemp mobile temp diff --git a/examples/PACKAGES/imd/in.bucky-plus-cnt-gpu b/examples/PACKAGES/imd/in.bucky-plus-cnt-gpu index 5762ec68c8..f3e4b32cdc 100644 --- a/examples/PACKAGES/imd/in.bucky-plus-cnt-gpu +++ b/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 nowait on +#fix comm all imd 6789 unwrap on trate 10 +fix comm all imd 6789 unwrap on trate 10 nowait on # temperature is based on mobile atoms only compute mobtemp mobile temp diff --git a/examples/PACKAGES/imd/in.deca-ala_imd-gpu b/examples/PACKAGES/imd/in.deca-ala_imd-gpu index 72c3f4aae9..9470f7c213 100644 --- a/examples/PACKAGES/imd/in.deca-ala_imd-gpu +++ b/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 diff --git a/examples/PACKAGES/imd/in.melt_imd-gpu b/examples/PACKAGES/imd/in.melt_imd-gpu index 24904eb832..f1406befa6 100644 --- a/examples/PACKAGES/imd/in.melt_imd-gpu +++ b/examples/PACKAGES/imd/in.melt_imd-gpu @@ -1,30 +1,32 @@ -# 3d Lennard-Jones melt +# 3d Lennard-Jones melt with GPU package acceleration -units lj -atom_style atomic -newton off +# enable GPU package from within the input: +package gpu 0 +suffix gpu -lattice fcc 0.8442 -region box block 0 10 0 10 0 10 -create_box 1 box -create_atoms 1 box -mass 1 1.0 +units lj +atom_style atomic -velocity all create 3.0 87287 +lattice fcc 0.8442 +region box block 0 10 0 10 0 10 +create_box 1 box +create_atoms 1 box +mass 1 1.0 -pair_style lj/cut/gpu 2.5 -pair_coeff 1 1 1.0 1.0 2.5 +velocity all create 3.0 87287 -neighbor 0.3 bin -neigh_modify every 5 delay 10 check yes +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +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 +fix 1 all nve # IMD setup. fix comm all imd 5678 unwrap off fscale 20.0 trate 20 nowait on -thermo 500 -run 5000000 +thermo 500 +run 5000000 diff --git a/src/KOKKOS/meam_dens_init_kokkos.h b/src/KOKKOS/meam_dens_init_kokkos.h index dd63be96bd..23985c7082 100644 --- a/src/KOKKOS/meam_dens_init_kokkos.h +++ b/src/KOKKOS/meam_dens_init_kokkos.h @@ -236,7 +236,6 @@ MEAMKokkos::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(d_rho0); diff --git a/src/KOKKOS/meam_force_kokkos.h b/src/KOKKOS/meam_force_kokkos.h index 1875e22dcf..703bee23d7 100644 --- a/src/KOKKOS/meam_force_kokkos.h +++ b/src/KOKKOS/meam_force_kokkos.h @@ -17,7 +17,6 @@ void MEAMKokkos::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; diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h index 399142dfaf..63e637c108 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/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; diff --git a/src/LATBOLTZ/fix_lb_fluid.cpp b/src/LATBOLTZ/fix_lb_fluid.cpp index 286b56cab5..d5161bba42 100644 --- a/src/LATBOLTZ/fix_lb_fluid.cpp +++ b/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); } } diff --git a/src/ML-POD/compute_podd_atom.cpp b/src/ML-POD/compute_podd_atom.cpp index 4ab6e23393..364bfb54ac 100644 --- a/src/ML-POD/compute_podd_atom.cpp +++ b/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)) - error->all(FLERR, "Per-atom data too large"); + if ((((3.0*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0)) + 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)) - error->all(FLERR, "Per-atom data too large"); + if ((((3.0*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0)) + 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; diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index e9632d254f..6354c67638 100644 --- a/src/compute_bond_local.cpp +++ b/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; diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index 57f15264f0..fa5d164844 100644 --- a/src/compute_pair_local.cpp +++ b/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; diff --git a/src/library.cpp b/src/library.cpp index e6b14dafd3..a8acbade52 100644 --- a/src/library.cpp +++ b/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 ` * :ref:`System property settings ` * :ref:`Neighbor topology data ` +* :ref:`Energy and virial tally settings ` * :ref:`Git revision and version settings ` * :ref:`Unit settings ` @@ -1984,6 +1990,35 @@ Get length of lists with :ref:`lammps_extract_setting() atom->q_flag; - if (strcmp(name,"neigh_bondlist") == 0) return lmp->neighbor->bondlist; - if (strcmp(name,"neigh_anglelist") == 0) return lmp->neighbor->anglelist; - if (strcmp(name,"neigh_dihedrallist") == 0) return lmp->neighbor->dihedrallist; - if (strcmp(name,"neigh_improperlist") == 0) return lmp->neighbor->improperlist; + if (strcmp(name,"neigh_bondlist") == 0) return (void *) lmp->neighbor->bondlist; + if (strcmp(name,"neigh_anglelist") == 0) return (void *) lmp->neighbor->anglelist; + if (strcmp(name,"neigh_dihedrallist") == 0) return (void *) lmp->neighbor->dihedrallist; + 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; diff --git a/src/potential_file_reader.cpp b/src/potential_file_reader.cpp index 2a93a4a524..7bac388ba0 100644 --- a/src/potential_file_reader.cpp +++ b/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() diff --git a/src/potential_file_reader.h b/src/potential_file_reader.h index c07b4b83f6..534457a2f8 100644 --- a/src/potential_file_reader.h +++ b/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(); diff --git a/src/read_data.cpp b/src/read_data.cpp index 79d88148c5..cf7b224db2 100644 --- a/src/read_data.cpp +++ b/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++; } diff --git a/src/reset_atoms_mol.cpp b/src/reset_atoms_mol.cpp index 54d3bbcc76..363e2f08eb 100644 --- a/src/reset_atoms_mol.cpp +++ b/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) { diff --git a/src/text_file_reader.cpp b/src/text_file_reader.cpp index 715e82ab32..d598c274f0 100644 --- a/src/text_file_reader.cpp +++ b/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) { diff --git a/unittest/c-library/test_library_properties.cpp b/unittest/c-library/test_library_properties.cpp index e2149d98e6..39b31e8217 100644 --- a/unittest/c-library/test_library_properties.cpp +++ b/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); - auto *b_ptr = (int64_t *)lammps_extract_global(lmp, "ntimestep"); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "ntimestep"), LAMMPS_BIGINT); + 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); diff --git a/unittest/python/python-numpy.py b/unittest/python/python-numpy.py index 4930527a61..be9109b9a3 100644 --- a/unittest/python/python-numpy.py +++ b/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")