diff --git a/doc/src/Build_basics.rst b/doc/src/Build_basics.rst index 761c743e2a..be4b312578 100644 --- a/doc/src/Build_basics.rst +++ b/doc/src/Build_basics.rst @@ -203,7 +203,7 @@ LAMMPS. check if the detected or selected compiler is compatible with the C++ support requirements of LAMMPS and stop with an error, if this is not the case. A C++11 compatible compiler is currently - required, but a transition to require C++17 is in progess and + required, but a transition to require C++17 is in progress and planned to be completed in Summer 2025. Currently, setting ``-DLAMMPS_CXX11=yes`` is required when configuring with CMake while using a C++11 compatible compiler that does not support C++17, @@ -329,7 +329,7 @@ LAMMPS. either as a binary package or through compiling from source. While a C++11 compatible compiler is currently sufficient to compile - LAMMPS, a transition to require C++17 is in progess and planned to + LAMMPS, a transition to require C++17 is in progress and planned to be completed in Summer 2025. Currently, setting ``-DLAMMPS_CXX11`` in the ``LMP_INC =`` line in the machine makefile is required when using a C++11 compatible compiler that does not support C++17. diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 75a7bb83cb..048a54ed37 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -115,7 +115,9 @@ OPT. * :doc:`gw/zbl ` * :doc:`harmonic/cut (o) ` * :doc:`hbond/dreiding/lj (o) ` + * :doc:`hbond/dreiding/lj/angleoffset (o) ` * :doc:`hbond/dreiding/morse (o) ` + * :doc:`hbond/dreiding/morse/angleoffset (o) ` * :doc:`hdnnp ` * :doc:`hippo (g) ` * :doc:`ilp/graphene/hbn (t) ` diff --git a/doc/src/Developer_code_design.rst b/doc/src/Developer_code_design.rst index 97cf0fc10e..9213efa18f 100644 --- a/doc/src/Developer_code_design.rst +++ b/doc/src/Developer_code_design.rst @@ -203,6 +203,7 @@ processed in the expected order before types are removed from dynamic dispatch. .. admonition:: Important Notes + :class: note In order to be able to detect incompatibilities at compile time and to avoid unexpected behavior, it is crucial that all member functions diff --git a/doc/src/Fortran.rst b/doc/src/Fortran.rst index bc641a237f..53dc5b8897 100644 --- a/doc/src/Fortran.rst +++ b/doc/src/Fortran.rst @@ -323,6 +323,12 @@ of the contents of the :f:mod:`LIBLAMMPS` Fortran interface to LAMMPS. :ftype set_internal_variable: subroutine :f eval: :f:func:`eval` :ftype eval: function + :f clearstep_compute: :f:subr:`clearstep_compute` + :ftype clearstep_compute: subroutine + :f addstep_compute: :f:subr:`addstep_compute` + :ftype addstep_compute: subroutine + :f addstep_compute_all: :f:subr:`addstep_compute_all` + :ftype addstep_compute_all: subroutine :f gather_atoms: :f:subr:`gather_atoms` :ftype gather_atoms: subroutine :f gather_atoms_concat: :f:subr:`gather_atoms_concat` @@ -956,6 +962,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type :f:func:`extract_atom` between runs. .. admonition:: Array index order + :class: tip Two-dimensional arrays returned from :f:func:`extract_atom` will be **transposed** from equivalent arrays in C, and they will be indexed @@ -1068,6 +1075,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type you based on data from the :cpp:class:`Compute` class. .. admonition:: Array index order + :class: tip Two-dimensional arrays returned from :f:func:`extract_compute` will be **transposed** from equivalent arrays in C, and they will be indexed @@ -1326,6 +1334,7 @@ Procedures Bound to the :f:type:`lammps` Derived Type :rtype data: polymorphic .. admonition:: Array index order + :class: tip Two-dimensional global, per-atom, or local array data from :f:func:`extract_fix` will be **transposed** from equivalent arrays in @@ -1450,11 +1459,62 @@ Procedures Bound to the :f:type:`lammps` Derived Type an internal-style variable, an error is generated. :p character(len=*) name: name of the variable - :p read(c_double) val: new value to assign to the variable + :p real(c_double) val: new value to assign to the variable :to: :cpp:func:`lammps_set_internal_variable` -------- +.. f:function:: eval(expr) + + This function is a wrapper around :cpp:func:`lammps_eval` that takes a + LAMMPS equal style variable string, evaluates it and returns the resulting + scalar value as a floating-point number. + + .. versionadded:: TBD + + :p character(len=\*) expr: string to be evaluated + :to: :cpp:func:`lammps_eval` + :r value [real(c_double)]: result of the evaluated string + +-------- + +.. f:subroutine:: clearstep_compute() + + Clear whether a compute has been invoked + + .. versionadded:: TBD + + :to: :cpp:func:`lammps_clearstep_compute` + +-------- + +.. f:subroutine:: addstep_compute(nextstep) + + Add timestep to list of future compute invocations + if the compute has been invoked on the current timestep + + .. versionadded:: TBD + + overloaded for 32-bit and 64-bit integer arguments + + :p integer(kind=8 or kind=4) nextstep: next timestep + :to: :cpp:func:`lammps_addstep_compute` + +-------- + +.. f:subroutine:: addstep_compute_all(nextstep) + + Add timestep to list of future compute invocations + + .. versionadded:: TBD + + overloaded for 32-bit and 64-bit integer arguments + + :p integer(kind=8 or kind=4) nextstep: next timestep + :to: :cpp:func:`lammps_addstep_compute_all` + +-------- + .. f:subroutine:: gather_atoms(name, count, data) This function calls :cpp:func:`lammps_gather_atoms` to gather the named diff --git a/doc/src/Howto_barostat.rst b/doc/src/Howto_barostat.rst index 0c25e2c53c..89ff4d183f 100644 --- a/doc/src/Howto_barostat.rst +++ b/doc/src/Howto_barostat.rst @@ -10,20 +10,21 @@ and/or pressure (P) is specified by the user, and the thermostat or barostat attempts to equilibrate the system to the requested T and/or P. -Barostatting in LAMMPS is performed by :doc:`fixes `. Two +Barostatting in LAMMPS is performed by :doc:`fixes `. Three barostatting methods are currently available: Nose-Hoover (npt and -nph) and Berendsen: +nph), Berendsen, and various linear controllers in deform/pressure: * :doc:`fix npt ` * :doc:`fix npt/sphere ` * :doc:`fix npt/asphere ` * :doc:`fix nph ` * :doc:`fix press/berendsen ` +* :doc:`fix deform/pressure ` The :doc:`fix npt ` commands include a Nose-Hoover thermostat and barostat. :doc:`Fix nph ` is just a Nose/Hoover barostat; -it does no thermostatting. Both :doc:`fix nph ` and :doc:`fix press/berendsen ` can be used in conjunction -with any of the thermostatting fixes. +it does no thermostatting. The fixes :doc:`nph `, :doc:`press/berendsen `, and :doc:`deform/pressure ` +can be used in conjunction with any of the thermostatting fixes. As with the :doc:`thermostats `, :doc:`fix npt ` and :doc:`fix nph ` only use translational motion of the @@ -44,9 +45,9 @@ a temperature or pressure compute to a barostatting fix. .. note:: As with the thermostats, the Nose/Hoover methods (:doc:`fix npt ` and :doc:`fix nph `) perform time integration. - :doc:`Fix press/berendsen ` does NOT, so it should - be used with one of the constant NVE fixes or with one of the NVT - fixes. + :doc:`Fix press/berendsen ` and :doc:`fix deform/pressure ` + do NOT, so they should be used with one of the constant NVE fixes or with + one of the NVT fixes. Thermodynamic output, which can be setup via the :doc:`thermo_style ` command, often includes pressure diff --git a/doc/src/Install_git.rst b/doc/src/Install_git.rst index 5108009a73..d01bc6a4c5 100644 --- a/doc/src/Install_git.rst +++ b/doc/src/Install_git.rst @@ -52,6 +52,7 @@ your machine and "release" is one of the 3 branches listed above. between them at any time using "git checkout ".) .. admonition:: Saving time and disk space when using ``git clone`` + :class: note The complete git history of the LAMMPS project is quite large because it contains the entire commit history of the project since fall 2006, diff --git a/doc/src/Library_objects.rst b/doc/src/Library_objects.rst index 53edfce9e6..86b6d253f0 100644 --- a/doc/src/Library_objects.rst +++ b/doc/src/Library_objects.rst @@ -13,6 +13,9 @@ fixes, or variables in LAMMPS using the following functions: - :cpp:func:`lammps_set_internal_variable` - :cpp:func:`lammps_variable_info` - :cpp:func:`lammps_eval` +- :cpp:func:`lammps_clearstep_compute` +- :cpp:func:`lammps_addstep_compute_all` +- :cpp:func:`lammps_addstep_compute` ----------------------- @@ -61,6 +64,21 @@ fixes, or variables in LAMMPS using the following functions: ----------------------- +.. doxygenfunction:: lammps_clearstep_compute + :project: progguide + +----------------------- + +.. doxygenfunction:: lammps_addstep_compute_all + :project: progguide + +----------------------- + +.. doxygenfunction:: lammps_addstep_compute + :project: progguide + +----------------------- + .. doxygenenum:: _LMP_DATATYPE_CONST .. doxygenenum:: _LMP_STYLE_CONST diff --git a/doc/src/fix_python_invoke.rst b/doc/src/fix_python_invoke.rst index ad55882270..4f33f5483b 100644 --- a/doc/src/fix_python_invoke.rst +++ b/doc/src/fix_python_invoke.rst @@ -66,6 +66,15 @@ gives access to the LAMMPS state from Python. from these callbacks, trying to execute input script commands will in the best case not work or in the worst case result in undefined behavior. +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options +are relevant to this fix. No global or per-atom quantities are stored +by this fix for access by various :doc:`output commands `. +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. This fix is not invoked during :doc:`energy minimization `. + Restrictions """""""""""" diff --git a/doc/src/fix_reaxff_species.rst b/doc/src/fix_reaxff_species.rst index 107695f0f6..068d1de995 100644 --- a/doc/src/fix_reaxff_species.rst +++ b/doc/src/fix_reaxff_species.rst @@ -200,8 +200,8 @@ The 2 values in the global vector are as follows: The per-atom vector stores the molecule ID for each atom as identified by the fix. If an atom is not in a molecule, its ID will be 0. For atoms in the same molecule, the molecule ID for all of them -will be the same and will be equal to the smallest atom ID of -any atom in the molecule. +will be the same, and molecule IDs will range from 1 to the number +of molecules. No parameter of this fix can be used with the *start/stop* keywords of the :doc:`run ` command. diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst index 0d864cce88..81a411ffc8 100644 --- a/doc/src/fix_wall_gran.rst +++ b/doc/src/fix_wall_gran.rst @@ -248,11 +248,11 @@ listed in the following table. | 8 | Radius :math:`r` of atom | distance units | +-------+----------------------------------------------------+----------------+ -If a granular submodel calculates additional contact information (e.g. the -heat submodels calculate the amount of heat exchanged), these quantities +If a granular sub-model calculates additional contact information (e.g. the +heat sub-models calculate the amount of heat exchanged), these quantities are appended to the end of this array. First, any extra values from the -normal submodel are appended followed by the damping, tangential, rolling, -twisting, then heat models. See the descriptions of granular submodels in +normal sub-model are appended followed by the damping, tangential, rolling, +twisting, then heat models. See the descriptions of granular sub-models in the :doc:`pair granular ` page for information on any extra quantities. diff --git a/doc/src/fix_wall_gran_region.rst b/doc/src/fix_wall_gran_region.rst index 0e5e98f1a8..4ad3b9d6c5 100644 --- a/doc/src/fix_wall_gran_region.rst +++ b/doc/src/fix_wall_gran_region.rst @@ -269,11 +269,11 @@ listed in the following table. | 8 | Radius :math:`r` of atom | distance units | +-------+----------------------------------------------------+----------------+ -If a granular submodel calculates additional contact information (e.g. the -heat submodels calculate the amount of heat exchanged), these quantities +If a granular sub-model calculates additional contact information (e.g. the +heat sub-models calculate the amount of heat exchanged), these quantities are appended to the end of this array. First, any extra values from the -normal submodel are appended followed by the damping, tangential, rolling, -twisting, then heat models. See the descriptions of granular submodels in +normal sub-model are appended followed by the damping, tangential, rolling, +twisting, then heat models. See the descriptions of granular sub-models in the :doc:`pair granular ` page for information on any extra quantities. diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 070e845c7f..f54079daf1 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -258,7 +258,7 @@ in damping model. The definition of multiple *mdr* models in the *pair_style* is currently not supported. Similarly, the *mdr* model cannot be combined with a different normal -model in the *pair_style*. Physically this means that only one homogenous +model in the *pair_style*. Physically this means that only one homogeneous collection of particles governed by a single *mdr* model is allowed. The *mdr* model currently only supports *fix wall/gran/region*, not @@ -303,6 +303,7 @@ radius in the *mdr* model, the keyword/arg pair *cutoff radius* must be specifie simulation involving 200 particles named *in.tableting.200*. The second is a triaxial compaction simulation involving 12 particles named *in.triaxial.compaction.12*. + ---------- In addition, the normal force is augmented by a damping term of the @@ -941,16 +942,16 @@ particle I. The next entry (8) is the magnitude of the rolling torque. The next entry (9) is the magnitude of the twisting torque acting about the vector connecting the two particle centers. The next 3 (10-12) are the components of the vector connecting -the centers of the two particles (x_I - x_J). If a granular submodel -calculates additional contact information (e.g. the heat submodels +the centers of the two particles (x_I - x_J). If a granular sub-model +calculates additional contact information (e.g. the heat sub-models calculate the amount of heat exchanged), these quantities are appended -to the end of this list. First, any extra values from the normal submodel +to the end of this list. First, any extra values from the normal sub-model are appended followed by the damping, tangential, rolling, twisting, then -heat models. See the descriptions of specific granular submodels above +heat models. See the descriptions of specific granular sub-models above for information on any extra quantities. If two or more models are defined by pair coefficients, the size of the array is set by the maximum number of extra quantities in a model but the order of quantities -is determined by each model's specific set of submodels. Any unused +is determined by each model's specific set of sub-models. Any unused quantities are zeroed. These extra quantities can be accessed by the :doc:`compute pair/local ` command, as *p1*, *p2*, ..., @@ -1046,7 +1047,7 @@ a bulk elastic response. Journal of the Mechanics and Physics of Solids, **(Zunker et al, 2025)** Zunker, W., Dunatunga, S., Thakur, S., Tang, P., & Kamrin, K. (2025). Experimentally validated DEM for large deformation powder compaction: mechanically-derived contact model and -screening of non-physical contacts. engrXiv. +screening of non-physical contacts. .. _Luding2008: diff --git a/doc/src/pair_hbond_dreiding.rst b/doc/src/pair_hbond_dreiding.rst index ce19ff9e38..e059fdf2ba 100644 --- a/doc/src/pair_hbond_dreiding.rst +++ b/doc/src/pair_hbond_dreiding.rst @@ -1,30 +1,46 @@ .. index:: pair_style hbond/dreiding/lj .. index:: pair_style hbond/dreiding/lj/omp +.. index:: pair_style hbond/dreiding/lj/angleoffset +.. index:: pair_style hbond/dreiding/lj/angleoffset/omp .. index:: pair_style hbond/dreiding/morse .. index:: pair_style hbond/dreiding/morse/omp +.. index:: pair_style hbond/dreiding/morse/angleoffset +.. index:: pair_style hbond/dreiding/morse/angleoffset/omp pair_style hbond/dreiding/lj command ==================================== Accelerator Variants: *hbond/dreiding/lj/omp* +pair_style hbond/dreiding/lj/angleoffset command +================================================ + +Accelerator Variants: *hbond/dreiding/lj/angleoffset/omp* + pair_style hbond/dreiding/morse command ======================================= Accelerator Variants: *hbond/dreiding/morse/omp* +pair_style hbond/dreiding/morse/angleoffset command +=================================================== + +Accelerator Variants: *hbond/dreiding/morse/angleoffset/omp* + + Syntax """""" .. code-block:: LAMMPS - pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff + pair_style style N inner_distance_cutoff outer_distance_cutoff angle_cutoff equilibrium_angle -* style = *hbond/dreiding/lj* or *hbond/dreiding/morse* +* style = *hbond/dreiding/lj* or *hbond/dreiding/morse* or *hbond/dreiding/lj/angleoffset* or *hbond/dreiding/morse/angleoffset* * N = power of cosine of angle theta (integer) * inner_distance_cutoff = global inner cutoff for Donor-Acceptor interactions (distance units) * outer_distance_cutoff = global cutoff for Donor-Acceptor interactions (distance units) * angle_cutoff = global angle cutoff for Acceptor-Hydrogen-Donor interactions (degrees) +* (with style angleoffset) equilibrium_angle = global equilibrium angle for Acceptor-Hydrogen-Donor interactions (degrees) Examples """""""" @@ -40,6 +56,9 @@ Examples labelmap atom 1 C 2 O 3 H pair_coeff C O hbond/dreiding/morse H i 3.88 1.7241379 2.9 2 9.0 11.0 90.0 + pair_style hybrid/overlay lj/cut 10.0 hbond/dreiding/lj 4 9.0 11.0 90 170.0 + pair_coeff 1 2 hbond/dreiding/lj 3 i 9.5 2.75 4 9.0 11.0 90.0 + Description """"""""""" @@ -74,42 +93,53 @@ hydrogen (H) and the donor atoms: .. image:: JPG/dreiding_hbond.jpg :align: center -These 3-body interactions can be defined for pairs of acceptor and -donor atoms, based on atom types. For each donor/acceptor atom pair, -the third atom in the interaction is a hydrogen permanently bonded to -the donor atom, e.g. in a bond list read in from a data file via the +These 3-body interactions can be defined for pairs of acceptor and donor +atoms, based on atom types. For each donor/acceptor atom pair, the +third atom in the interaction is a hydrogen permanently bonded to the +donor atom, e.g. in a bond list read in from a data file via the :doc:`read_data ` command. The atom types of possible hydrogen atoms for each donor/acceptor type pair are specified by the :doc:`pair_coeff ` command (see below). Style *hbond/dreiding/lj* is the original DREIDING potential of -:ref:`(Mayo) `. It uses a LJ 12/10 functional for the Donor-Acceptor -interactions. To match the results in the original paper, use n = 4. +:ref:`(Mayo) `. It uses a LJ 12/10 functional for the +Donor-Acceptor interactions. To match the results in the original paper, +use n = 4. Style *hbond/dreiding/morse* is an improved version using a Morse potential for the Donor-Acceptor interactions. :ref:`(Liu) ` showed that the Morse form gives improved results for Dendrimer simulations, when n = 2. +.. versionadded:: TBD + +The style variants *hbond/dreiding/lj/angleoffset* and +*hbond/dreiding/lj/angleoffset* take the equilibrium angle of the AHD as +input, allowing it to reach 180 degrees. This variant option was added +to account for cases (especially in some coarse-grained models) in which +the equilibrium state of the bonds may equal the minimum energy state. + See the :doc:`Howto bioFF ` page for more information on the DREIDING force field. .. note:: - Because the Dreiding hydrogen bond potential is only one portion - of an overall force field which typically includes other pairwise - interactions, it is common to use it as a sub-style in a :doc:`pair_style hybrid/overlay ` command, where another pair style - provides the repulsive core interaction between pairs of atoms, e.g. a - 1/r\^12 Lennard-Jones repulsion. + Because the Dreiding hydrogen bond potential is only one portion of + an overall force field which typically includes other pairwise + interactions, it is common to use it as a sub-style in a + :doc:`pair_style hybrid/overlay ` command, where another + pair style provides the repulsive core interaction between pairs of + atoms, e.g. a 1/r\^12 Lennard-Jones repulsion. .. note:: - When using the hbond/dreiding pair styles with :doc:`pair_style hybrid/overlay `, you should explicitly define pair + When using the hbond/dreiding pair styles with :doc:`pair_style + hybrid/overlay `, you should explicitly define pair interactions between the donor atom and acceptor atoms, (as well as between these atoms and ALL other atoms in your system). Whenever - :doc:`pair_style hybrid/overlay ` is used, ordinary mixing - rules are not applied to atoms like the donor and acceptor atoms - because they are typically referenced in multiple pair styles. + :doc:`pair_style hybrid/overlay ` is used, ordinary + mixing rules are not applied to atoms like the donor and acceptor + atoms because they are typically referenced in multiple pair styles. Neglecting to do this can cause difficult-to-detect physics problems. .. note:: @@ -119,6 +149,13 @@ on the DREIDING force field. special_bonds command (e.g. "special_bonds lj 0.0 0.0 1.0") to turn these interactions on. +.. note:: + + For the *angleoffset* variants, the referenced angle offset is the + supplementary angle of the equilibrium angle parameter. It means if + the equilibrium angle is 166.6 degrees, the calculated angle offset + is 13.4 degrees. + ---------- The following coefficients must be defined for pairs of eligible @@ -169,7 +206,10 @@ follows: * distance cutoff :math:`r_{out}` (distance units) * angle cutoff (degrees) -A single hydrogen atom type K can be specified, or a wild-card asterisk +For both the *hbond/dreiding/lj/angleoffset* and *hbond/dreiding/morse/angleoffset* styles an additional parameter is added: +* equilibrium angle (degrees) + +For all styles, a single hydrogen atom type K can be specified, or a wild-card asterisk can be used in place of or in conjunction with the K arguments to select multiple types as hydrogen atoms. This takes the form "\*" or "\*n" or "n\*" or "m\*n". See the :doc:`pair_coeff ` @@ -245,8 +285,7 @@ heading) the following commands could be included in an input script: Restrictions """""""""""" -This pair style can only be used if LAMMPS was built with the -MOLECULE package. See the :doc:`Build package ` doc page +The base pair styles can only be used if LAMMPS was built with the MOLECULE package. The *angleoffset* variant also requires the EXTRA-MOLECULE package. See the :doc:`Build package ` doc page for more info. Related commands diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst index e325de0aa6..7f8b36bb83 100644 --- a/doc/src/pair_mliap.rst +++ b/doc/src/pair_mliap.rst @@ -145,6 +145,7 @@ per line. The detail of *nn* module implementation can be found at :ref:`(Yanxon) `. .. admonition:: Notes on mliappy models + :class: note When the *model* keyword is *mliappy*, if the filename ends in '.pt', or '.pth', it will be loaded using pytorch; otherwise, it will be diff --git a/doc/src/pair_reaxff.rst b/doc/src/pair_reaxff.rst index 495572dc0e..45532bc2a6 100644 --- a/doc/src/pair_reaxff.rst +++ b/doc/src/pair_reaxff.rst @@ -158,11 +158,36 @@ drops to zero. Optional keywords *safezone*, *mincap*, and *minhbonds* are used for allocating reaxff arrays. Increasing these values can avoid memory problems, such as segmentation faults and bondchk failed errors, that -could occur under certain conditions. These keywords are not used by +could occur under certain conditions. These keywords are **not** used by the Kokkos version, which instead uses a more robust memory allocation scheme that checks if the sizes of the arrays have been exceeded and automatically allocates more memory. +.. admonition:: Memory management problems with ReaxFF + :class: tip + + The LAMMPS implementation of ReaxFF is adapted from a standalone MD + program written in C called `PuReMD + `_. It inherits from this code + a heuristic memory management that is different from what the rest of + LAMMPS uses. It assumes that a system is dense and already well + equilibrated, so that there are no large changes in how many and what + types of neighbors atoms have. However, not all systems are like + that, and thus there can be errors or segmentation faults if the + system changes too much. If you run into problems, here are three + options to avoid them: + + - Use the KOKKOS version of ReaxFF (KOKKOS is not only for GPUs, + but can also be compiled for serial or OpenMP execution) which + uses a different memory management approach. + - Break down a run command during which memory related errors happen + into multiple smaller segments so that the memory management + heuristics are re-initialized for each segment before they become + invalid. + - Increase the values for *safezone*, *mincap*, and *minhbonds* as + needed. This can lead to significant increase of memory consumption + through. + The keyword *tabulate* controls the size of interpolation table for Lennard-Jones and Coulomb interactions. Tabulation may also be set in the control file (see below). If tabulation is set in both the input script and the diff --git a/doc/src/pair_style.rst b/doc/src/pair_style.rst index 93e02422d2..bdf06d6b66 100644 --- a/doc/src/pair_style.rst +++ b/doc/src/pair_style.rst @@ -207,7 +207,9 @@ accelerated styles exist. * :doc:`gw/zbl ` - Gao-Weber potential with a repulsive ZBL core * :doc:`harmonic/cut ` - repulsive-only harmonic potential * :doc:`hbond/dreiding/lj ` - DREIDING hydrogen bonding LJ potential +* :doc:`hbond/dreiding/lj/angleoffset ` - DREIDING hydrogen bonding LJ potential with offset for hbond angle * :doc:`hbond/dreiding/morse ` - DREIDING hydrogen bonding Morse potential +* :doc:`hbond/dreiding/morse/angleoffset ` - DREIDING hydrogen bonding Morse potential with offset for hbond angle * :doc:`hdnnp ` - High-dimensional neural network potential * :doc:`hippo ` - * :doc:`ilp/graphene/hbn ` - registry-dependent interlayer potential (ILP) diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 6746d591c6..f9fb3a0f8a 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -108,6 +108,7 @@ Andrienko Andzelm Ang anglegrad +angleoffset angletangrad angmom angmomx @@ -1583,6 +1584,7 @@ Impropers imulator includelink incompressible +incompressibility incrementing indenter indenters @@ -1762,6 +1764,7 @@ Kadiri Kai Kalia Kamberaj +Kamrin Kantorovich Kapfer Kapil @@ -2532,6 +2535,7 @@ Nevery Nevins newfile Newns +newstep newtype nextsort Neyts @@ -3731,6 +3735,7 @@ tgnpt tgnvt th Thakkar +Thakur Thaokar thb thei @@ -3769,6 +3774,7 @@ Tigran Tij Tildesley Timan +timeflag timeI timespan timestamp @@ -3831,6 +3837,7 @@ Tref Tretyakov tri triangleflag +triaxial Tribello triclinic Triclinic diff --git a/examples/COUPLE/plugin/liblammpsplugin.c b/examples/COUPLE/plugin/liblammpsplugin.c index 99e38df32b..619b8828fc 100644 --- a/examples/COUPLE/plugin/liblammpsplugin.c +++ b/examples/COUPLE/plugin/liblammpsplugin.c @@ -118,6 +118,9 @@ liblammpsplugin_t *liblammpsplugin_load(const char *lib) ADDSYM(set_internal_variable); ADDSYM(variable_info); ADDSYM(eval); + ADDSYM(clearstep_compute); + ADDSYM(addstep_compute); + ADDSYM(addstep_compute_all); ADDSYM(gather_atoms); ADDSYM(gather_atoms_concat); diff --git a/examples/COUPLE/plugin/liblammpsplugin.h b/examples/COUPLE/plugin/liblammpsplugin.h index dd1c9628f5..c765b0adc3 100644 --- a/examples/COUPLE/plugin/liblammpsplugin.h +++ b/examples/COUPLE/plugin/liblammpsplugin.h @@ -164,6 +164,9 @@ struct _liblammpsplugin { int (*set_internal_variable)(void *, const char *, double); int (*variable_info)(void *, int, char *, int); double (*eval)(void *, const char *); + void (*clearstep_compute)(void *); + void (*addstep_compute)(void *, void *); + void (*addstep_compute_all)(void *, void *); void (*gather_atoms)(void *, const char *, int, int, void *); void (*gather_atoms_concat)(void *, const char *, int, int, void *); diff --git a/examples/PACKAGES/stressprofile/in.flat b/examples/PACKAGES/stressprofile/in.flat index 8b484a423f..a0aea07d83 100644 --- a/examples/PACKAGES/stressprofile/in.flat +++ b/examples/PACKAGES/stressprofile/in.flat @@ -32,7 +32,7 @@ fix 1 all nvt temp 0.7 0.7 0.2 #dump_modify 3 pad 3 fix 2 all recenter NULL NULL 15 units lattice -compute p1 all stress/cartesian z 0.5 +compute p1 all stress/cartesian z 0.5 NULL 0 fix 3 all ave/time 100 1 100 c_p1[*] file flat.out mode vector thermo 50 diff --git a/examples/multi/in.granular b/examples/multi/in.granular index 468d0dcbf1..6abbd8da3b 100644 --- a/examples/multi/in.granular +++ b/examples/multi/in.granular @@ -1,4 +1,4 @@ -# Big colloid particles and small LJ particles +# Binary granular system units lj atom_style sphere diff --git a/fortran/lammps.f90 b/fortran/lammps.f90 index 552b3dfad3..2cfd4422b0 100644 --- a/fortran/lammps.f90 +++ b/fortran/lammps.f90 @@ -127,6 +127,16 @@ MODULE LIBLAMMPS PROCEDURE :: set_string_variable => lmp_set_string_variable PROCEDURE :: set_internal_variable => lmp_set_internal_variable PROCEDURE :: eval => lmp_eval + + PROCEDURE :: clearstep_compute => lmp_clearstep_compute + PROCEDURE, PRIVATE :: lmp_addstep_compute_smallint + PROCEDURE, PRIVATE :: lmp_addstep_compute_bigint + GENERIC :: addstep_compute => lmp_addstep_compute_smallint, lmp_addstep_compute_bigint + PROCEDURE, PRIVATE :: lmp_addstep_compute_all_smallint + PROCEDURE, PRIVATE :: lmp_addstep_compute_all_bigint + GENERIC :: addstep_compute_all => lmp_addstep_compute_all_smallint, & + lmp_addstep_compute_all_bigint + PROCEDURE, PRIVATE :: lmp_gather_atoms_int PROCEDURE, PRIVATE :: lmp_gather_atoms_double GENERIC :: gather_atoms => lmp_gather_atoms_int, & @@ -626,6 +636,24 @@ MODULE LIBLAMMPS REAL(c_double) :: lammps_eval END FUNCTION lammps_eval + SUBROUTINE lammps_clearstep_compute(handle) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle + END SUBROUTINE lammps_clearstep_compute + + SUBROUTINE lammps_addstep_compute(handle, step) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, step + END SUBROUTINE lammps_addstep_compute + + SUBROUTINE lammps_addstep_compute_all(handle, step) BIND(C) + IMPORT :: c_ptr + IMPLICIT NONE + TYPE(c_ptr), VALUE :: handle, step + END SUBROUTINE lammps_addstep_compute_all + SUBROUTINE lammps_gather_atoms(handle, name, TYPE, count, DATA) BIND(C) IMPORT :: c_int, c_ptr IMPLICIT NONE @@ -1846,6 +1874,80 @@ CONTAINS CALL lammps_free(Cexpr) END FUNCTION lmp_eval + ! equivalent subroutine to lammps_clearstep_compute + SUBROUTINE lmp_clearstep_compute(self) + CLASS(lammps), INTENT(IN) :: self + CALL lammps_clearstep_compute(self%handle) + END SUBROUTINE lmp_clearstep_compute + + ! equivalent subroutine to lammps_addstep_compute + SUBROUTINE lmp_addstep_compute_bigint(self, nextstep) + CLASS(lammps), INTENT(IN) :: self + INTEGER(kind=8), INTENT(IN) :: nextstep + INTEGER(c_int), TARGET :: smallstep + INTEGER(c_int64_t), TARGET :: bigstep + TYPE(c_ptr) :: ptrstep + IF (SIZE_BIGINT == 4_c_int) THEN + smallstep = INT(nextstep,kind=c_int) + ptrstep = C_LOC(smallstep) + ELSE + bigstep = nextstep + ptrstep = C_LOC(bigstep) + END IF + CALL lammps_addstep_compute(self%handle, ptrstep) + END SUBROUTINE lmp_addstep_compute_bigint + + ! equivalent subroutine to lammps_addstep_compute + SUBROUTINE lmp_addstep_compute_smallint(self, nextstep) + CLASS(lammps), INTENT(IN) :: self + INTEGER(kind=4), INTENT(IN) :: nextstep + INTEGER(c_int), TARGET :: smallstep + INTEGER(c_int64_t), TARGET :: bigstep + TYPE(c_ptr) :: ptrstep + IF (SIZE_BIGINT == 4_c_int) THEN + smallstep = nextstep + ptrstep = C_LOC(smallstep) + ELSE + bigstep = nextstep + ptrstep = C_LOC(bigstep) + END IF + CALL lammps_addstep_compute(self%handle, ptrstep) + END SUBROUTINE lmp_addstep_compute_smallint + + ! equivalent subroutine to lammps_addstep_compute_all + SUBROUTINE lmp_addstep_compute_all_bigint(self, nextstep) + CLASS(lammps), INTENT(IN) :: self + INTEGER(kind=8), INTENT(IN) :: nextstep + INTEGER(c_int), TARGET :: smallstep + INTEGER(c_int64_t), TARGET :: bigstep + TYPE(c_ptr) :: ptrstep + IF (SIZE_BIGINT == 4_c_int) THEN + smallstep = INT(nextstep,kind=c_int) + ptrstep = C_LOC(smallstep) + ELSE + bigstep = nextstep + ptrstep = C_LOC(bigstep) + END IF + CALL lammps_addstep_compute_all(self%handle, ptrstep) + END SUBROUTINE lmp_addstep_compute_all_bigint + + ! equivalent subroutine to lammps_addstep_compute_all + SUBROUTINE lmp_addstep_compute_all_smallint(self, nextstep) + CLASS(lammps), INTENT(IN) :: self + INTEGER(kind=4), INTENT(IN) :: nextstep + INTEGER(c_int), TARGET :: smallstep + INTEGER(c_int64_t), TARGET :: bigstep + TYPE(c_ptr) :: ptrstep + IF (SIZE_BIGINT == 4_c_int) THEN + smallstep = nextstep + ptrstep = C_LOC(smallstep) + ELSE + bigstep = nextstep + ptrstep = C_LOC(bigstep) + END IF + CALL lammps_addstep_compute_all(self%handle, ptrstep) + END SUBROUTINE lmp_addstep_compute_all_smallint + ! equivalent function to lammps_gather_atoms (for integers) SUBROUTINE lmp_gather_atoms_int(self, name, count, data) CLASS(lammps), INTENT(IN) :: self diff --git a/python/lammps/core.py b/python/lammps/core.py index 89a3dd2bf0..089517bab7 100644 --- a/python/lammps/core.py +++ b/python/lammps/core.py @@ -422,6 +422,10 @@ class lammps(object): self.lib.lammps_extract_variable_datatype.argtypes = [c_void_p, c_char_p] self.lib.lammps_extract_variable_datatype.restype = c_int + self.lib.lammps_clearstep_compute.argtype = [c_void_p] + self.lib.lammps_addstep_compute.argtype = [c_void_p, c_void_p] + self.lib.lammps_addstep_compute_all.argtype = [c_void_p, c_void_p] + self.lib.lammps_eval.argtypes = [c_void_p, c_char_p] self.lib.lammps_eval.restype = c_double @@ -1594,6 +1598,26 @@ class lammps(object): # ------------------------------------------------------------------------- + def clearstep_compute(self, nextstep): + with ExceptionCheck(self): + return self.lib.lammps_clearstep_compute(self.lmp) + + # ------------------------------------------------------------------------- + + def addstep_compute(self, nextstep): + with ExceptionCheck(self): + nextstep = self.c_bigint(nextstep) + return self.lib.lammps_addstep_compute(self.lmp, POINTER(nextstep)) + + # ------------------------------------------------------------------------- + + def addstep_compute_all(self, nextstep): + with ExceptionCheck(self): + nextstep = self.c_bigint(nextstep) + return self.lib.lammps_addstep_compute_all(self.lmp, POINTER(nextstep)) + + # ------------------------------------------------------------------------- + def flush_buffers(self): """Flush output buffers diff --git a/src/.gitignore b/src/.gitignore index 9431d55913..3630d940de 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -1281,6 +1281,10 @@ /pair_hbond_dreiding_lj.h /pair_hbond_dreiding_morse.cpp /pair_hbond_dreiding_morse.h +/pair_hbond_dreiding_lj_angleoffset.cpp +/pair_hbond_dreiding_lj_angleoffset.h +/pair_hbond_dreiding_morse_angleoffset.cpp +/pair_hbond_dreiding_morse_angleoffset.h /pair_hdnnp.cpp /pair_hdnnp.h /pair_ilp_graphene_hbn.cpp diff --git a/src/EXTRA-FIX/fix_electron_stopping_fit.cpp b/src/EXTRA-FIX/fix_electron_stopping_fit.cpp index 2657ddc85f..c9cca2679d 100644 --- a/src/EXTRA-FIX/fix_electron_stopping_fit.cpp +++ b/src/EXTRA-FIX/fix_electron_stopping_fit.cpp @@ -68,7 +68,8 @@ FixElectronStoppingFit::FixElectronStoppingFit(LAMMPS *lmp, int narg, char **arg error->all(FLERR,"Incorrect number of fix electron/stopping/fit arguments"); } - scalar_flag = 1; + scalar_flag = 1; // intensive total energy loss since start of run + extscalar = 0; global_freq = 1; energy_coh_in = new double[atom->ntypes+1]; diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp new file mode 100644 index 0000000000..02ac009773 --- /dev/null +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.cpp @@ -0,0 +1,131 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Tod A Pascal (Caltech), Don Xu/EiPi Fun +------------------------------------------------------------------------- */ + +#include "pair_hbond_dreiding_lj_angleoffset.h" + +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "math_special.h" +#include "memory.h" +#include "molecule.h" +#include "neigh_list.h" +#include "neighbor.h" + +#include +#include + +static constexpr double SMALL = 0.001; +static constexpr int CHUNK = 8; + +using namespace LAMMPS_NS; +using namespace MathConst; +// using namespace MathSpecial; + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingLJAngleoffset::PairHbondDreidingLJAngleoffset(LAMMPS *lmp) + : PairHbondDreidingLJ(lmp) { + + angle_offset_flag = 1; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHbondDreidingLJAngleoffset::coeff(int narg, char **arg) +{ + if (narg < 6 || narg > 11) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi,klo,khi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM); + + int donor_flag; + if (strcmp(arg[3],"i") == 0) donor_flag = 0; + else if (strcmp(arg[3],"j") == 0) donor_flag = 1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + + double epsilon_one = utils::numeric(FLERR, arg[4], false, lmp); + double sigma_one = utils::numeric(FLERR, arg[5], false, lmp); + + int ap_one = ap_global; + if (narg > 6) ap_one = utils::inumeric(FLERR, arg[6], false, lmp); + double cut_inner_one = cut_inner_global; + double cut_outer_one = cut_outer_global; + if (narg > 8) { + cut_inner_one = utils::numeric(FLERR, arg[7], false, lmp); + cut_outer_one = utils::numeric(FLERR, arg[8], false, lmp); + } + if (cut_inner_one>cut_outer_one) + error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); + double cut_angle_one = cut_angle_global; + if (narg > 9) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0; + double angle_offset_one = angle_offset_global; + if (narg == 11) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[10], false, lmp)) * MY_PI/180.0; + if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) + error->all(FLERR,"Illegal angle offset"); + + // grow params array if necessary + + if (nparams == maxparam) { + maxparam += CHUNK; + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), + "pair:params"); + + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind + + memset(params + nparams, 0, CHUNK*sizeof(Param)); + } + + params[nparams].epsilon = epsilon_one; + params[nparams].sigma = sigma_one; + params[nparams].ap = ap_one; + params[nparams].cut_inner = cut_inner_one; + params[nparams].cut_outer = cut_outer_one; + params[nparams].cut_innersq = cut_inner_one*cut_inner_one; + params[nparams].cut_outersq = cut_outer_one*cut_outer_one; + params[nparams].cut_angle = cut_angle_one; + params[nparams].angle_offset = angle_offset_one; + params[nparams].denom_vdw = + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq); + + // flag type2param with either i,j = D,A or j,i = D,A + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + for (int k = klo; k <= khi; k++) { + if (donor_flag == 0) type2param[i][j][k] = nparams; + else type2param[j][i][k] = nparams; + count++; + } + nparams++; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h new file mode 100644 index 0000000000..5ae1fcba28 --- /dev/null +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_lj_angleoffset.h @@ -0,0 +1,38 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(hbond/dreiding/lj/angleoffset,PairHbondDreidingLJAngleoffset); +// clang-format on +#else + +#ifndef LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_H +#define LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_H + +#include "pair_hbond_dreiding_lj.h" + +namespace LAMMPS_NS { + +class PairHbondDreidingLJAngleoffset : public PairHbondDreidingLJ { + + public: + PairHbondDreidingLJAngleoffset(class LAMMPS *); + void coeff(int, char **) override; + +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp new file mode 100644 index 0000000000..21f26ea8d1 --- /dev/null +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.cpp @@ -0,0 +1,129 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Tod A Pascal (Caltech), Don Xu/EiPi Fun +------------------------------------------------------------------------- */ + +#include "pair_hbond_dreiding_morse_angleoffset.h" + +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "math_special.h" +#include "memory.h" +#include "molecule.h" +#include "neigh_list.h" +#include "neighbor.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace MathConst; + +static constexpr int CHUNK = 8; + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingMorseAngleoffset::PairHbondDreidingMorseAngleoffset(LAMMPS *lmp) : + PairHbondDreidingMorse(lmp) { + + angle_offset_flag = 1; +} + +/* ---------------------------------------------------------------------- + * set coeffs for one or more type pairs + * ---------------------------------------------------------------------- */ + +void PairHbondDreidingMorseAngleoffset::coeff(int narg, char **arg) +{ + if (narg < 7 || narg > 12) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi,klo,khi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, lmp, Atom::ATOM); + + int donor_flag; + if (strcmp(arg[3],"i") == 0) donor_flag = 0; + else if (strcmp(arg[3],"j") == 0) donor_flag = 1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + + double d0_one = utils::numeric(FLERR, arg[4], false, lmp); + double alpha_one = utils::numeric(FLERR, arg[5], false, lmp); + double r0_one = utils::numeric(FLERR, arg[6], false, lmp); + + int ap_one = ap_global; + if (narg > 7) ap_one = utils::inumeric(FLERR, arg[7], false, lmp); + double cut_inner_one = cut_inner_global; + double cut_outer_one = cut_outer_global; + if (narg > 9) { + cut_inner_one = utils::numeric(FLERR, arg[8], false, lmp); + cut_outer_one = utils::numeric(FLERR, arg[9], false, lmp); + } + if (cut_inner_one>cut_outer_one) + error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); + double cut_angle_one = cut_angle_global; + if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, lmp) * MY_PI/180.0; + double angle_offset_one = angle_offset_global; + if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[11], false, lmp)) * MY_PI/180.0; + if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) + error->all(FLERR,"Illegal angle offset {}", angle_offset_one); + + // grow params array if necessary + + if (nparams == maxparam) { + maxparam += CHUNK; + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),"pair:params"); + + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind + + memset(params + nparams, 0, CHUNK*sizeof(Param)); + } + + params[nparams].d0 = d0_one; + params[nparams].alpha = alpha_one; + params[nparams].r0 = r0_one; + params[nparams].ap = ap_one; + params[nparams].cut_inner = cut_inner_one; + params[nparams].cut_outer = cut_outer_one; + params[nparams].cut_innersq = cut_inner_one*cut_inner_one; + params[nparams].cut_outersq = cut_outer_one*cut_outer_one; + params[nparams].cut_angle = cut_angle_one; + params[nparams].angle_offset = angle_offset_one; + params[nparams].denom_vdw = (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq); + + // flag type2param with either i,j = D,A or j,i = D,A + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + for (int k = klo; k <= khi; k++) { + if (donor_flag == 0) type2param[i][j][k] = nparams; + else type2param[j][i][k] = nparams; + count++; + } + nparams++; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} diff --git a/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h new file mode 100644 index 0000000000..20d2a8698d --- /dev/null +++ b/src/EXTRA-MOLECULE/pair_hbond_dreiding_morse_angleoffset.h @@ -0,0 +1,38 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(hbond/dreiding/morse/angleoffset,PairHbondDreidingMorseAngleoffset); +// clang-format on +#else + +#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H +#define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_H + +#include "pair_hbond_dreiding_morse.h" + +namespace LAMMPS_NS { + +class PairHbondDreidingMorseAngleoffset : public PairHbondDreidingMorse { + + public: + PairHbondDreidingMorseAngleoffset(class LAMMPS *); + void coeff(int, char **) override; + +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/MISC/pair_srp.cpp b/src/MISC/pair_srp.cpp index 8decae5aca..1fe19bc30c 100644 --- a/src/MISC/pair_srp.cpp +++ b/src/MISC/pair_srp.cpp @@ -397,10 +397,11 @@ void PairSRP::coeff(int narg, char **arg) error->all(FLERR,"PairSRP: Incorrect args for pair coeff"); if (!allocated) allocate(); - if (btype_str.size() > 0) + if (btype_str.size() > 0) { btype = utils::expand_type_int(FLERR, btype_str, Atom::BOND, lmp); - if ((btype > atom->nbondtypes) || (btype <= 0)) - error->all(FLERR,"Invalid bond type {} for pair style srp", btype); + if ((btype > atom->nbondtypes) || (btype <= 0)) + error->all(FLERR,"Invalid bond type {} for pair style srp", btype); + } if (bptype_str.size() > 0) bptype = utils::expand_type_int(FLERR, bptype_str, Atom::ATOM, lmp); diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index 4536cc8e05..fd0e61edd2 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -13,7 +13,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing author: Tod A Pascal (Caltech) + Contributing authors: Tod A Pascal (Caltech), Don Xu/EiPi Fun ------------------------------------------------------------------------- */ #include "pair_hbond_dreiding_lj.h" @@ -55,6 +55,9 @@ PairHbondDreidingLJ::PairHbondDreidingLJ(LAMMPS *lmp) : Pair(lmp) nextra = 2; pvector = new double[2]; + + angle_offset_flag = 0; + angle_offset_global = 0.0; } /* ---------------------------------------------------------------------- */ @@ -82,7 +85,7 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) tagint tagprev; double delx,dely,delz,rsq,rsq1,rsq2,r1,r2; double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2,d; + double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; double fi[3],fj[3],delr1[3],delr2[3]; double r2inv,r10inv; double switch1,switch2; @@ -178,6 +181,13 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -186,15 +196,12 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) r2inv = 1.0/rsq; r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; - force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * - powint(c,pm.ap); - force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * - powint(c,pm.ap-1)*s; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * powint(c,pm.ap-1)*s; + force_switch = 0.0; eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); - force_switch=0.0; - if (rsq > pm.cut_innersq) { switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / @@ -300,12 +307,19 @@ void PairHbondDreidingLJ::allocate() void PairHbondDreidingLJ::settings(int narg, char **arg) { - if (narg != 4) error->all(FLERR,"Illegal pair_style command"); + + // narg = 4 for standard form, narg = 5 or 6 if angleoffset LJ or Morse variants respectively (from EXTRA-MOLECULE) + if (narg != 4 && narg != 5) error->all(FLERR,"Illegal pair_style command"); ap_global = utils::inumeric(FLERR,arg[0],false,lmp); cut_inner_global = utils::numeric(FLERR,arg[1],false,lmp); cut_outer_global = utils::numeric(FLERR,arg[2],false,lmp); cut_angle_global = utils::numeric(FLERR,arg[3],false,lmp) * MY_PI/180.0; + + // update when using angleoffset variant + if (angle_offset_flag) { + angle_offset_global = (180.0 - utils::numeric(FLERR, arg[4], false, lmp)) * MY_PI/180.0; + } } /* ---------------------------------------------------------------------- @@ -314,8 +328,14 @@ void PairHbondDreidingLJ::settings(int narg, char **arg) void PairHbondDreidingLJ::coeff(int narg, char **arg) { - if (narg < 6 || narg > 10) + // account for angleoffset variant in EXTRA-MOLECULE + int maxarg = 10; + if (angle_offset_flag == 1) maxarg = 11; + + // check settings + if (narg < 6 || narg > maxarg) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); int ilo,ihi,jlo,jhi,klo,khi; @@ -343,16 +363,15 @@ void PairHbondDreidingLJ::coeff(int narg, char **arg) error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); double cut_angle_one = cut_angle_global; if (narg == 10) cut_angle_one = utils::numeric(FLERR, arg[9], false, lmp) * MY_PI/180.0; + // grow params array if necessary if (nparams == maxparam) { maxparam += CHUNK; - params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), - "pair:params"); + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), "pair:params"); // make certain all addional allocated storage is initialized // to avoid false positives when checking with valgrind - memset(params + nparams, 0, CHUNK*sizeof(Param)); } @@ -540,6 +559,13 @@ double PairHbondDreidingLJ::single(int i, int j, int itype, int jtype, if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.h b/src/MOLECULE/pair_hbond_dreiding_lj.h index 91a906c5cb..131f3946b3 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.h +++ b/src/MOLECULE/pair_hbond_dreiding_lj.h @@ -25,6 +25,7 @@ PairStyle(hbond/dreiding/lj,PairHbondDreidingLJ); namespace LAMMPS_NS { class PairHbondDreidingLJ : public Pair { + public: PairHbondDreidingLJ(class LAMMPS *); ~PairHbondDreidingLJ() override; @@ -38,6 +39,8 @@ class PairHbondDreidingLJ : public Pair { protected: double cut_inner_global, cut_outer_global, cut_angle_global; int ap_global; + int angle_offset_flag; // 1 if angle offset variant used + double angle_offset_global; // updated if angle offset variant used struct Param { double epsilon, sigma; @@ -45,7 +48,7 @@ class PairHbondDreidingLJ : public Pair { double d0, alpha, r0; double morse1; double denom_vdw; - double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset; + double cut_inner, cut_outer, cut_innersq, cut_outersq, cut_angle, offset, angle_offset; int ap; }; diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index d976b66460..8506765984 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -148,6 +148,13 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -158,6 +165,7 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) dr = r - pm.r0; dexp = exp(-pm.alpha * dr); eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; force_switch = 0.0; @@ -196,12 +204,12 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) vz1 = a11*delr1[2] + a12*delr2[2]; vz2 = a22*delr2[2] + a12*delr1[2]; - fi[0] = vx1 + (b+d)*delx; - fi[1] = vy1 + (b+d)*dely; - fi[2] = vz1 + (b+d)*delz; - fj[0] = vx2 - (b+d)*delx; - fj[1] = vy2 - (b+d)*dely; - fj[2] = vz2 - (b+d)*delz; + fi[0] = vx1 + b*delx + d*delx; + fi[1] = vy1 + b*dely + d*dely; + fi[2] = vz1 + b*delz + d*delz; + fj[0] = vx2 - b*delx - d*delx; + fj[1] = vy2 - b*dely - d*dely; + fj[2] = vz2 - b*delz - d*delz; f[i][0] += fi[0]; f[i][1] += fi[1]; @@ -238,7 +246,9 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) void PairHbondDreidingMorse::coeff(int narg, char **arg) { - if (narg < 7 || narg > 11) + int maxarg = 12; + if (angle_offset_flag == 1) maxarg = 12; + if (narg < 7 || narg > maxarg) error->all(FLERR,"Incorrect args for pair coefficients"); if (!allocated) allocate(); @@ -443,6 +453,13 @@ double PairHbondDreidingMorse::single(int i, int j, int itype, int jtype, if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac < pm.cut_angle || ac > (2.0*MY_PI - pm.cut_angle)) return 0.0; s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp new file mode 100644 index 0000000000..11c09ed549 --- /dev/null +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.cpp @@ -0,0 +1,127 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun +------------------------------------------------------------------------- */ + +#include "pair_hbond_dreiding_lj_angleoffset_omp.h" + +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "domain.h" +#include "force.h" +#include "math_const.h" +#include "math_special.h" +#include "memory.h" +#include "molecule.h" +#include "neigh_list.h" +#include "suffix.h" + +#include + +#include "omp_compat.h" +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecial; + +static constexpr int CHUNK = 8; + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingLJAngleoffsetOMP::PairHbondDreidingLJAngleoffsetOMP(LAMMPS *lmp) : + PairHbondDreidingLJOMP(lmp) { + angle_offset_flag = 1; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairHbondDreidingLJAngleoffsetOMP::coeff(int narg, char **arg) +{ + auto mylmp = PairHbondDreidingLJ::lmp; + if (narg < 6 || narg > 11) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi,klo,khi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, mylmp, Atom::ATOM); + + int donor_flag; + if (strcmp(arg[3],"i") == 0) donor_flag = 0; + else if (strcmp(arg[3],"j") == 0) donor_flag = 1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + + double epsilon_one = utils::numeric(FLERR, arg[4], false, mylmp); + double sigma_one = utils::numeric(FLERR, arg[5], false, mylmp); + + int ap_one = ap_global; + if (narg > 6) ap_one = utils::inumeric(FLERR, arg[6], false, mylmp); + double cut_inner_one = cut_inner_global; + double cut_outer_one = cut_outer_global; + if (narg > 8) { + cut_inner_one = utils::numeric(FLERR, arg[7], false, mylmp); + cut_outer_one = utils::numeric(FLERR, arg[8], false, mylmp); + } + if (cut_inner_one>cut_outer_one) + error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); + double cut_angle_one = cut_angle_global; + if (narg > 9) cut_angle_one = utils::numeric(FLERR, arg[9], false, mylmp) * MY_PI/180.0; + double angle_offset_one = angle_offset_global; + if (narg == 11) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[10], false, mylmp)) * MY_PI/180.0; + if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) + error->all(FLERR,"Illegal angle offset"); + + // grow params array if necessary + + if (nparams == maxparam) { + maxparam += CHUNK; + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param), + "pair:params"); + + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind + + memset(params + nparams, 0, CHUNK*sizeof(Param)); + } + + params[nparams].epsilon = epsilon_one; + params[nparams].sigma = sigma_one; + params[nparams].ap = ap_one; + params[nparams].cut_inner = cut_inner_one; + params[nparams].cut_outer = cut_outer_one; + params[nparams].cut_innersq = cut_inner_one*cut_inner_one; + params[nparams].cut_outersq = cut_outer_one*cut_outer_one; + params[nparams].cut_angle = cut_angle_one; + params[nparams].angle_offset = angle_offset_one; + params[nparams].denom_vdw = + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq); + + // flag type2param with either i,j = D,A or j,i = D,A + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + for (int k = klo; k <= khi; k++) { + if (donor_flag == 0) type2param[i][j][k] = nparams; + else type2param[j][i][k] = nparams; + count++; + } + nparams++; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} diff --git a/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h new file mode 100644 index 0000000000..03d3392e4d --- /dev/null +++ b/src/OPENMP/pair_hbond_dreiding_lj_angleoffset_omp.h @@ -0,0 +1,41 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(hbond/dreiding/lj/angleoffset/omp,PairHbondDreidingLJAngleoffsetOMP); +// clang-format on +#else + +#ifndef LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_OMP_H +#define LMP_PAIR_HBOND_DREIDING_LJ_ANGLEOFFSET_OMP_H + +#include "pair_hbond_dreiding_lj_omp.h" + +namespace LAMMPS_NS { + +class PairHbondDreidingLJAngleoffsetOMP : public PairHbondDreidingLJOMP { + + public: + PairHbondDreidingLJAngleoffsetOMP(class LAMMPS *); + void coeff(int, char **) override; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp b/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp index b0f6dcfb5b..60ec8938fe 100644 --- a/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_lj_omp.cpp @@ -120,15 +120,16 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,iatom,imol; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2; - double factor_hb,force_angle,force_kernel,evdwl,eng_lj; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; + double factor_hb,force_angle,force_kernel,evdwl,eng_lj,ehbond,force_switch; + double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; double fi[3],fj[3],delr1[3],delr2[3]; double r2inv,r10inv; double switch1,switch2; int *ilist,*jlist,*numneigh,**firstneigh; const tagint *klist; - evdwl = 0.0; + evdwl = ehbond = 0.0; + int hbcount = 0; const auto * _noalias const x = (dbl3_t *) atom->x[0]; auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; @@ -151,9 +152,6 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) // jj = loop over acceptors // kk = loop over hydrogens bonded to donor - int hbcount = 0; - double hbeng = 0.0; - for (ii = iifrom; ii < iito; ++ii) { i = ilist[ii]; itype = type[i]; @@ -222,6 +220,13 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -230,30 +235,34 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) r2inv = 1.0/rsq; r10inv = r2inv*r2inv*r2inv*r2inv*r2inv; - force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * - powint(c,pm.ap); - force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * - powint(c,pm.ap-1)*s; + force_kernel = r10inv*(pm.lj1*r2inv - pm.lj2)*r2inv * powint(c,pm.ap); + force_angle = pm.ap * r10inv*(pm.lj3*r2inv - pm.lj4) * powint(c,pm.ap-1)*s; + force_switch = 0.0; eng_lj = r10inv*(pm.lj3*r2inv - pm.lj4); + if (rsq > pm.cut_innersq) { switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * (pm.cut_outersq + 2.0*rsq - 3.0*pm.cut_innersq) / pm.denom_vdw; switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * (rsq-pm.cut_innersq) / pm.denom_vdw; - force_kernel = force_kernel*switch1 + eng_lj*switch2/rsq; - force_angle *= switch1; - eng_lj *= switch1; + + force_kernel *= switch1; + force_angle *= switch1; + force_switch = eng_lj*switch2/rsq; + eng_lj *= switch1; } if (EFLAG) { evdwl = eng_lj * powint(c,pm.ap); evdwl *= factor_hb; + ehbond += evdwl; } a = factor_hb*force_angle/s; b = factor_hb*force_kernel; + d = factor_hb*force_switch; a11 = a*c / rsq1; a12 = -a / (r1*r2); @@ -266,12 +275,12 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) vz1 = a11*delr1[2] + a12*delr2[2]; vz2 = a22*delr2[2] + a12*delr1[2]; - fi[0] = vx1 + b*delx; - fi[1] = vy1 + b*dely; - fi[2] = vz1 + b*delz; - fj[0] = vx2 - b*delx; - fj[1] = vy2 - b*dely; - fj[2] = vz2 - b*delz; + fi[0] = vx1 + b*delx + d*delx; + fi[1] = vy1 + b*dely + d*dely; + fi[2] = vz1 + b*delz + d*delz; + fj[0] = vx2 - b*delx - d*delx; + fj[1] = vy2 - b*dely - d*dely; + fj[2] = vz2 - b*delz - d*delz; fxtmp += fi[0]; fytmp += fi[1]; @@ -290,7 +299,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr); if (EFLAG) { hbcount++; - hbeng += evdwl; + ehbond += evdwl; } } } @@ -302,7 +311,7 @@ void PairHbondDreidingLJOMP::eval(int iifrom, int iito, ThrData * const thr) } const int tid = thr->get_tid(); hbcount_thr[tid] = static_cast(hbcount); - hbeng_thr[tid] = hbeng; + hbeng_thr[tid] = ehbond; } /* ---------------------------------------------------------------------- */ diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp new file mode 100644 index 0000000000..e7c75f29e4 --- /dev/null +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.cpp @@ -0,0 +1,127 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + This software is distributed under the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun +------------------------------------------------------------------------- */ + +#include "pair_hbond_dreiding_morse_angleoffset_omp.h" + +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "domain.h" +#include "force.h" +#include "math_const.h" +#include "math_special.h" +#include "memory.h" +#include "molecule.h" +#include "neigh_list.h" +#include "suffix.h" + +#include + +#include "omp_compat.h" +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecial; + +static constexpr int CHUNK = 8; + +/* ---------------------------------------------------------------------- */ + +PairHbondDreidingMorseAngleoffsetOMP::PairHbondDreidingMorseAngleoffsetOMP(LAMMPS *lmp) : + PairHbondDreidingMorseOMP(lmp) { + angle_offset_flag = 1; +} + +/* ---------------------------------------------------------------------- + * set coeffs for one or more type pairs + * ---------------------------------------------------------------------- */ + +void PairHbondDreidingMorseAngleoffsetOMP::coeff(int narg, char **arg) +{ + auto mylmp = PairHbondDreidingMorse::lmp; + if (narg < 7 || narg > 12) + error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi,klo,khi; + utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); + utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error); + utils::bounds_typelabel(FLERR, arg[2], 1, atom->ntypes, klo, khi, mylmp, Atom::ATOM); + + int donor_flag; + if (strcmp(arg[3],"i") == 0) donor_flag = 0; + else if (strcmp(arg[3],"j") == 0) donor_flag = 1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + + double d0_one = utils::numeric(FLERR, arg[4], false, mylmp); + double alpha_one = utils::numeric(FLERR, arg[5], false, mylmp); + double r0_one = utils::numeric(FLERR, arg[6], false, mylmp); + + int ap_one = ap_global; + if (narg > 7) ap_one = utils::inumeric(FLERR, arg[7], false, mylmp); + double cut_inner_one = cut_inner_global; + double cut_outer_one = cut_outer_global; + if (narg > 9) { + cut_inner_one = utils::numeric(FLERR, arg[8], false, mylmp); + cut_outer_one = utils::numeric(FLERR, arg[9], false, mylmp); + } + if (cut_inner_one>cut_outer_one) + error->all(FLERR,"Pair inner cutoff >= Pair outer cutoff"); + double cut_angle_one = cut_angle_global; + if (narg > 10) cut_angle_one = utils::numeric(FLERR, arg[10], false, mylmp) * MY_PI/180.0; + double angle_offset_one = angle_offset_global; + if (narg == 12) angle_offset_one = (180.0 - utils::numeric(FLERR, arg[11], false, mylmp)) * MY_PI/180.0; + if (angle_offset_one < 0.0 || angle_offset_one > 90.0 * MY_PI/180.0) + error->all(FLERR,"Illegal angle offset {}", angle_offset_one); + + // grow params array if necessary + + if (nparams == maxparam) { + maxparam += CHUNK; + params = (Param *) memory->srealloc(params, maxparam*sizeof(Param),"pair:params"); + + // make certain all addional allocated storage is initialized + // to avoid false positives when checking with valgrind + + memset(params + nparams, 0, CHUNK*sizeof(Param)); + } + + params[nparams].d0 = d0_one; + params[nparams].alpha = alpha_one; + params[nparams].r0 = r0_one; + params[nparams].ap = ap_one; + params[nparams].cut_inner = cut_inner_one; + params[nparams].cut_outer = cut_outer_one; + params[nparams].cut_innersq = cut_inner_one*cut_inner_one; + params[nparams].cut_outersq = cut_outer_one*cut_outer_one; + params[nparams].cut_angle = cut_angle_one; + params[nparams].angle_offset = angle_offset_one; + params[nparams].denom_vdw = (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq) * + (params[nparams].cut_outersq-params[nparams].cut_innersq); + + // flag type2param with either i,j = D,A or j,i = D,A + + int count = 0; + for (int i = ilo; i <= ihi; i++) + for (int j = MAX(jlo,i); j <= jhi; j++) + for (int k = klo; k <= khi; k++) { + if (donor_flag == 0) type2param[i][j][k] = nparams; + else type2param[j][i][k] = nparams; + count++; + } + nparams++; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} diff --git a/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h new file mode 100644 index 0000000000..40a797ac33 --- /dev/null +++ b/src/OPENMP/pair_hbond_dreiding_morse_angleoffset_omp.h @@ -0,0 +1,42 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (Temple U), Don Xu/EiPi Fun +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(hbond/dreiding/morse/angleoffset/omp,PairHbondDreidingMorseAngleoffsetOMP); +// clang-format on +#else + +#ifndef LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_OMP_H +#define LMP_PAIR_HBOND_DREIDING_MORSE_ANGLEOFFSET_OMP_H + +#include "pair_hbond_dreiding_morse_omp.h" + +namespace LAMMPS_NS { + +class PairHbondDreidingMorseAngleoffsetOMP : + public PairHbondDreidingMorseOMP { + + public: + PairHbondDreidingMorseAngleoffsetOMP(class LAMMPS *); + void coeff(int, char **) override; +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp b/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp index 0e43e2a037..7772ea69fa 100644 --- a/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp +++ b/src/OPENMP/pair_hbond_dreiding_morse_omp.cpp @@ -120,14 +120,14 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) int i,j,k,m,ii,jj,kk,jnum,knum,itype,jtype,ktype,imol,iatom; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq,rsq1,rsq2,r1,r2; - double factor_hb,force_angle,force_kernel,evdwl; - double c,s,a,b,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; + double factor_hb,force_angle,force_kernel,force_switch,evdwl,ehbond; + double c,s,a,b,d,ac,a11,a12,a22,vx1,vx2,vy1,vy2,vz1,vz2; double fi[3],fj[3],delr1[3],delr2[3]; double r,dr,dexp,eng_morse,switch1,switch2; int *ilist,*jlist,*numneigh,**firstneigh; const tagint *klist; - evdwl = 0.0; + evdwl = ehbond = 0.0; const auto * _noalias const x = (dbl3_t *) atom->x[0]; auto * _noalias const f = (dbl3_t *) thr->get_f()[0]; @@ -151,7 +151,6 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) // kk = loop over hydrogens bonded to donor int hbcount = 0; - double hbeng = 0.0; for (ii = iifrom; ii < iito; ++ii) { @@ -222,6 +221,13 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) if (c < -1.0) c = -1.0; ac = acos(c); + if (angle_offset_flag){ + ac = ac + pm.angle_offset; + c = cos(ac); + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + } + if (ac > pm.cut_angle && ac < (2.0*MY_PI - pm.cut_angle)) { s = sqrt(1.0 - c*c); if (s < SMALL) s = SMALL; @@ -232,8 +238,10 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) dr = r - pm.r0; dexp = exp(-pm.alpha * dr); eng_morse = pm.d0 * (dexp*dexp - 2.0*dexp); + force_kernel = pm.morse1*(dexp*dexp - dexp)/r * powint(c,pm.ap); force_angle = pm.ap * eng_morse * powint(c,pm.ap-1)*s; + force_switch = 0.0; if (rsq > pm.cut_innersq) { switch1 = (pm.cut_outersq-rsq) * (pm.cut_outersq-rsq) * @@ -241,18 +249,22 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) pm.denom_vdw; switch2 = 12.0*rsq * (pm.cut_outersq-rsq) * (rsq-pm.cut_innersq) / pm.denom_vdw; - force_kernel = force_kernel*switch1 + eng_morse*switch2/rsq; + + force_kernel *= switch1; force_angle *= switch1; + force_switch = eng_morse*switch2/rsq; eng_morse *= switch1; } if (EFLAG) { evdwl = eng_morse * powint(c,pm.ap); evdwl *= factor_hb; + ehbond += evdwl; } a = factor_hb*force_angle/s; b = factor_hb*force_kernel; + d = factor_hb*force_switch; a11 = a*c / rsq1; a12 = -a / (r1*r2); @@ -265,12 +277,12 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) vz1 = a11*delr1[2] + a12*delr2[2]; vz2 = a22*delr2[2] + a12*delr1[2]; - fi[0] = vx1 + b*delx; - fi[1] = vy1 + b*dely; - fi[2] = vz1 + b*delz; - fj[0] = vx2 - b*delx; - fj[1] = vy2 - b*dely; - fj[2] = vz2 - b*delz; + fi[0] = vx1 + b*delx + d*delx; + fi[1] = vy1 + b*dely + d*dely; + fi[2] = vz1 + b*delz + d*delz; + fj[0] = vx2 - b*delx - d*delx; + fj[1] = vy2 - b*dely - d*dely; + fj[2] = vz2 - b*delz - d*delz; fxtmp += fi[0]; fytmp += fi[1]; @@ -289,7 +301,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) if (EVFLAG) ev_tally3_thr(this,k,i,j,evdwl,0.0,fi,fj,delr1,delr2,thr); if (EFLAG) { hbcount++; - hbeng += evdwl; + ehbond += evdwl; } } } @@ -301,7 +313,7 @@ void PairHbondDreidingMorseOMP::eval(int iifrom, int iito, ThrData * const thr) } const int tid = thr->get_tid(); hbcount_thr[tid] = static_cast(hbcount); - hbeng_thr[tid] = hbeng; + hbeng_thr[tid] = ehbond; } /* ---------------------------------------------------------------------- */ diff --git a/src/PYTHON/fix_python_invoke.cpp b/src/PYTHON/fix_python_invoke.cpp index 7fd3ad88f7..6f27831438 100644 --- a/src/PYTHON/fix_python_invoke.cpp +++ b/src/PYTHON/fix_python_invoke.cpp @@ -22,6 +22,7 @@ #include "lmppython.h" #include "python_compat.h" #include "python_utils.h" +#include "modify.h" #include "update.h" #include @@ -70,6 +71,12 @@ FixPythonInvoke::FixPythonInvoke(LAMMPS *lmp, int narg, char **arg) : } lmpPtr = PY_VOID_POINTER(lmp); + + // nvalid = next step on which end_of_step or post_force does something + // add nextvalid() to all computes that store invocation times + // since we don't know a priori which are invoked by python code + nvalid = nextvalid(); + modify->addstep_compute_all(nvalid); } /* ---------------------------------------------------------------------- */ @@ -89,8 +96,23 @@ int FixPythonInvoke::setmask() /* ---------------------------------------------------------------------- */ +void FixPythonInvoke::init() +{ + // need to reset nvalid if nvalid < ntimestep b/c minimize was performed + + if (nvalid < update->ntimestep) { + nvalid = nextvalid(); + modify->addstep_compute_all(nvalid); + } +} + +/* ---------------------------------------------------------------------- */ + void FixPythonInvoke::end_of_step() { + // python code may invoke computes so wrap with clear/add + modify->clearstep_compute(); + PyUtils::GIL lock; PyObject * result = PyObject_CallFunction((PyObject*)pFunc, (char *)"O", (PyObject*)lmpPtr); @@ -101,6 +123,9 @@ void FixPythonInvoke::end_of_step() } Py_CLEAR(result); + + nvalid = nextvalid(); + modify->addstep_compute(nvalid); } /* ---------------------------------------------------------------------- */ @@ -116,6 +141,9 @@ void FixPythonInvoke::post_force(int vflag) { if (update->ntimestep % nevery != 0) return; + // python code may invoke computes so wrap with clear/add + modify->clearstep_compute(); + PyUtils::GIL lock; char fmt[] = "Oi"; @@ -127,4 +155,14 @@ void FixPythonInvoke::post_force(int vflag) } Py_CLEAR(result); + + nvalid = nextvalid(); + modify->addstep_compute(nvalid); +} + +/* ---------------------------------------------------------------------- */ + +bigint FixPythonInvoke::nextvalid() +{ + return (update->ntimestep/nevery + 1)*nevery; } diff --git a/src/PYTHON/fix_python_invoke.h b/src/PYTHON/fix_python_invoke.h index 09382e5780..fdb931b69c 100644 --- a/src/PYTHON/fix_python_invoke.h +++ b/src/PYTHON/fix_python_invoke.h @@ -31,6 +31,7 @@ class FixPythonInvoke : public Fix { ~FixPythonInvoke() override; int setmask() override; void setup(int) override; + void init() override; void end_of_step() override; void post_force(int) override; @@ -38,6 +39,8 @@ class FixPythonInvoke : public Fix { void *lmpPtr; void *pFunc; int selected_callback; + bigint nextvalid(); + bigint nvalid; }; } // namespace LAMMPS_NS diff --git a/src/REAXFF/fix_reaxff_species.cpp b/src/REAXFF/fix_reaxff_species.cpp index 1916597a69..09c3884c34 100644 --- a/src/REAXFF/fix_reaxff_species.cpp +++ b/src/REAXFF/fix_reaxff_species.cpp @@ -26,6 +26,7 @@ #include "domain.h" #include "error.h" #include "fix_ave_atom.h" +#include "fix_property_atom.h" #include "force.h" #include "group.h" #include "input.h" @@ -141,13 +142,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : } x0 = nullptr; - clusterID = nullptr; - - int ntmp = atom->nmax; - memory->create(x0, ntmp, "reaxff/species:x0"); - memory->create(clusterID, ntmp, "reaxff/species:clusterID"); - memset(clusterID, 0, sizeof(double) * ntmp); - vector_atom = clusterID; nmax = 0; setupflag = 0; @@ -304,7 +298,6 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) : FixReaxFFSpecies::~FixReaxFFSpecies() { memory->destroy(BOCut); - memory->destroy(clusterID); memory->destroy(x0); memory->destroy(nd); @@ -330,6 +323,7 @@ FixReaxFFSpecies::~FixReaxFFSpecies() try { modify->delete_compute(fmt::format("SPECATOM_{}", id)); modify->delete_fix(fmt::format("SPECBOND_{}", id)); + modify->delete_fix(fmt::format("clusterID_{}", id)); } catch (std::exception &) { } } @@ -372,9 +366,6 @@ void FixReaxFFSpecies::init() reaxff->fixspecies_flag = 1; - // reset next output timestep if not yet set or timestep has been reset - if (nvalid != update->ntimestep) nvalid = update->ntimestep + nfreq; - if (!setupflag) { // create a compute to store properties modify->add_compute(fmt::format("SPECATOM_{} all SPEC/ATOM q x y z vx vy vz abo01 abo02 " @@ -387,11 +378,25 @@ void FixReaxFFSpecies::init() auto fixcmd = fmt::format("SPECBOND_{} all ave/atom {} {} {}", id, nevery, nrepeat, nfreq); for (int i = 1; i < 32; ++i) fixcmd += fmt::format(" c_SPECATOM_{}[{}]", id, i); f_SPECBOND = dynamic_cast(modify->add_fix(fixcmd)); + + // create a fix to point to fix_property_atom for storing clusterID + fixcmd = fmt::format("clusterID_{} all property/atom d_clusterID ghost yes", id); + f_clusterID = dynamic_cast(modify->add_fix(fixcmd)); + + // per-atom property for clusterID + int flag,cols; + int index1 = atom->find_custom("clusterID",flag,cols); + clusterID = atom->dvector[index1]; + vector_atom = clusterID; + + int ntmp = atom->nmax; + memory->create(x0, ntmp, "reaxff/species:x0"); + setupflag = 1; } // check for valid variable name for delete Nlimit keyword - if (delete_Nsteps > 0) { + if (delete_Nsteps > 0 && delete_Nlimit_varid > -1) { delete_Nlimit_varid = input->variable->find(delete_Nlimit_varname.c_str()); if (delete_Nlimit_varid < 0) error->all(FLERR, "Fix reaxff/species: Variable name {} does not exist", @@ -424,10 +429,16 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/) { int Nmole, Nspec; + // per-atom property for clusterID + int flag,cols; + int index1 = atom->find_custom("clusterID",flag,cols); + clusterID = atom->dvector[index1]; + vector_atom = clusterID; + // point to fix_ave_atom f_SPECBOND->end_of_step(); - if (ntimestep != nvalid) { + if (ntimestep != nvalid && nvalid != -1) { // push back delete_Tcount on every step if (delete_Nsteps > 0) for (int i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1]; @@ -439,11 +450,7 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/) if (atom->nmax > nmax) { nmax = atom->nmax; memory->destroy(x0); - memory->destroy(clusterID); memory->create(x0, nmax, "reaxff/species:x0"); - memory->create(clusterID, nmax, "reaxff/species:clusterID"); - memset(clusterID, 0, sizeof(double) * nmax); - vector_atom = clusterID; } for (int i = 0; i < nmax; i++) { x0[i].x = x0[i].y = x0[i].z = 0.0; } @@ -464,9 +471,14 @@ void FixReaxFFSpecies::Output_ReaxFF_Bonds(bigint ntimestep, FILE * /*fp*/) if (comm->me == 0) fflush(pos); } - if (delflag) DeleteSpecies(Nmole, Nspec); + if (delflag && nvalid != -1) { + DeleteSpecies(Nmole, Nspec); - nvalid += nfreq; + // reset molecule ID to index from 1 + SortMolecule(Nmole); + } + + nvalid = ntimestep + nfreq; } /* ---------------------------------------------------------------------- */ @@ -826,7 +838,8 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) int count, count_tmp, m, n, k; int *Nameall; int *mask = atom->mask; - double avq, avq_tmp, avx[3], avx_tmp, box[3], halfbox[3]; + double *rmass = atom->rmass; + double totq, totq_tmp, com[3], com_tmp, thism, totm, box[3], halfbox[3]; double **spec_atom = f_SPECBOND->array_atom; if (multipos) OpenPos(); @@ -844,7 +857,7 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) update->ntimestep, Nmole, Nspec, domain->boxlo[0], domain->boxhi[0], domain->boxlo[1], domain->boxhi[1], domain->boxlo[2], domain->boxhi[2]); - fprintf(pos, "ID\tAtom_Count\tType\tAve_q\t\tCoM_x\t\tCoM_y\t\tCoM_z\n"); + fprintf(pos, "ID\tAtom_Count\tType\tTot_q\t\tCoM_x\t\tCoM_y\t\tCoM_z\n"); } Nameall = nullptr; @@ -853,8 +866,9 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) for (m = 1; m <= Nmole; m++) { count = 0; - avq = 0.0; - for (n = 0; n < 3; n++) avx[n] = 0.0; + totq = 0.0; + totm = 0.0; + for (n = 0; n < 3; n++) com[n] = 0.0; for (n = 0; n < nutypes; n++) Name[n] = 0; for (i = 0; i < nlocal; i++) { @@ -864,30 +878,37 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) itype = ele2uele[atom->type[i] - 1]; Name[itype]++; count++; - avq += spec_atom[i][0]; + totq += spec_atom[i][0]; if ((x0[i].x - spec_atom[i][1]) > halfbox[0]) spec_atom[i][1] += box[0]; if ((spec_atom[i][1] - x0[i].x) > halfbox[0]) spec_atom[i][1] -= box[0]; if ((x0[i].y - spec_atom[i][2]) > halfbox[1]) spec_atom[i][2] += box[1]; if ((spec_atom[i][2] - x0[i].y) > halfbox[1]) spec_atom[i][2] -= box[1]; if ((x0[i].z - spec_atom[i][3]) > halfbox[2]) spec_atom[i][3] += box[2]; if ((spec_atom[i][3] - x0[i].z) > halfbox[2]) spec_atom[i][3] -= box[2]; - for (n = 0; n < 3; n++) avx[n] += spec_atom[i][n + 1]; + if (rmass) thism = rmass[i]; + else thism = atom->mass[atom->type[i]]; + for (n = 0; n < 3; n++) com[n] += spec_atom[i][n+1]*thism; + totm += thism; } } - avq_tmp = 0.0; - MPI_Allreduce(&avq, &avq_tmp, 1, MPI_DOUBLE, MPI_SUM, world); - avq = avq_tmp; + totq_tmp = 0.0; + MPI_Allreduce(&totq, &totq_tmp, 1, MPI_DOUBLE, MPI_SUM, world); + totq = totq_tmp; for (n = 0; n < 3; n++) { - avx_tmp = 0.0; - MPI_Reduce(&avx[n], &avx_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world); - avx[n] = avx_tmp; + com_tmp = 0.0; + MPI_Reduce(&com[n], &com_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world); + com[n] = com_tmp; } MPI_Reduce(&count, &count_tmp, 1, MPI_INT, MPI_SUM, 0, world); count = count_tmp; + com_tmp = 0.0; + MPI_Reduce(&totm, &com_tmp, 1, MPI_DOUBLE, MPI_SUM, 0, world); + totm = com_tmp; + MPI_Reduce(Name, Nameall, nutypes, MPI_INT, MPI_SUM, 0, world); for (n = 0; n < nutypes; n++) Name[n] = Nameall[n]; @@ -900,16 +921,15 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) } } if (count > 0) { - avq /= count; for (k = 0; k < 3; k++) { - avx[k] /= count; - if (avx[k] >= domain->boxhi[k]) avx[k] -= box[k]; - if (avx[k] < domain->boxlo[k]) avx[k] += box[k]; + com[k] /= totm; + if (com[k] >= domain->boxhi[k]) com[k] -= box[k]; + if (com[k] < domain->boxlo[k]) com[k] += box[k]; - avx[k] -= domain->boxlo[k]; - avx[k] /= box[k]; + com[k] -= domain->boxlo[k]; + com[k] /= box[k]; } - fprintf(pos, "\t%.8f \t%.8f \t%.8f \t%.8f", avq, avx[0], avx[1], avx[2]); + fprintf(pos, "\t%.8f \t%.8f \t%.8f \t%.8f", totq, com[0], com[1], com[2]); } fprintf(pos, "\n"); } @@ -922,21 +942,29 @@ void FixReaxFFSpecies::WritePos(int Nmole, int Nspec) void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) { - int ndeletions; + int i, ndeletions; int headroom = -1; if (delete_Nsteps > 0) { - if (delete_Tcount[delete_Nsteps - 1] == -1) return; + if (delete_Tcount[delete_Nsteps - 1] == -1) { + for (i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1]; + return; + } ndeletions = delete_Tcount[0] - delete_Tcount[delete_Nsteps - 1]; if (delete_Nlimit_varid > -1) delete_Nlimit = input->variable->compute_equal(delete_Nlimit_varid); headroom = MAX(0, delete_Nlimit - ndeletions); - if (headroom == 0) return; + if (headroom == 0) { + for (i = delete_Nsteps - 1; i > 0; i--) delete_Tcount[i] = delete_Tcount[i - 1]; + return; + } } - int i, j, m, n, itype, cid; + int j, m, n, itype, cid; int ndel, ndelone, count, count_tmp; int *Nameall; int *mask = atom->mask; + double *mass = atom->mass; + double *rmass = atom->rmass; double localmass, totalmass; std::string species_str; @@ -989,7 +1017,8 @@ void FixReaxFFSpecies::DeleteSpecies(int Nmole, int Nspec) Name[itype]++; count++; marklist[nmarklist++] = i; - localmass += atom->mass[atom->type[i]]; + if (rmass) localmass += rmass[i]; + else localmass += atom->mass[atom->type[i]]; } } @@ -1177,7 +1206,7 @@ double FixReaxFFSpecies::memory_usage() { double bytes; - bytes = 4 * nmax * sizeof(double); // clusterID + x0 + bytes = 3 * nmax * sizeof(double); // x0 return bytes; } diff --git a/src/REAXFF/fix_reaxff_species.h b/src/REAXFF/fix_reaxff_species.h index b9afc5466a..d378065a82 100644 --- a/src/REAXFF/fix_reaxff_species.h +++ b/src/REAXFF/fix_reaxff_species.h @@ -88,6 +88,7 @@ class FixReaxFFSpecies : public Fix { class NeighList *list; class FixAveAtom *f_SPECBOND; + class FixPropertyAtom *f_clusterID; class PairReaxFF *reaxff; }; } // namespace LAMMPS_NS diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index d972f19ffa..121473f7e3 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -273,11 +273,13 @@ double FixRHEOPressure::calc_rho(double p, int i) error->one(FLERR, "Rho calculation from pressure not yet supported for cubic pressure equation"); } else if (pressure_style[type] == TAITWATER) { - rho = pow(7.0 * p + csq[type] * rho0[type], SEVENTH); + double tmp = 7.0 * p + csq[type] * rho0[type]; + rho = pow(MAX(0.0, tmp), SEVENTH); rho *= pow(rho0[type], 6.0 * SEVENTH); rho *= pow(csq[type], -SEVENTH); } else if (pressure_style[type] == TAITGENERAL) { - rho = pow(tpower[type] * p + csq[type] * rho0[type], 1.0 / tpower[type]); + double tmp = tpower[type] * p + csq[type] * rho0[type]; + rho = pow(MAX(0.0, tmp), 1.0 / tpower[type]); rho *= pow(rho0[type], 1.0 - 1.0 / tpower[type]); rho *= pow(csq[type], -1.0 / tpower[type]); } else if (pressure_style[type] == IDEAL) { diff --git a/src/library.cpp b/src/library.cpp index 2cd9879e76..07ed9184bc 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -2925,6 +2925,8 @@ int lammps_variable_info(void *handle, int idx, char *buffer, int buf_size) { return 0; } +/* ---------------------------------------------------------------------- */ + /** Evaluate an immediate variable expression * \verbatim embed:rst @@ -2958,6 +2960,86 @@ double lammps_eval(void *handle, const char *expr) return result; } +/* ---------------------------------------------------------------------- */ + +/** Clear whether a compute has been invoked. + * +\verbatim embed:rst + +.. versionadded:: TBD + + This function clears the invoked flag of all computes. + Called everywhere that computes are used, before computes are invoked. + The invoked flag is used to avoid re-invoking same compute multiple times + and to flag computes that store invocation times as having been invoked + +*See also* + :cpp:func:`lammps_addstep_compute_all` + :cpp:func:`lammps_addstep_compute` + +\endverbatim + + * \param handle pointer to a previously created LAMMPS instance cast to ``void *``. + */ +void lammps_clearstep_compute(void *handle) { + auto lmp = (LAMMPS *) handle; + lmp->modify->clearstep_compute(); +} + +/* ---------------------------------------------------------------------- */ + +/** Add next timestep to all computes + * +\verbatim embed:rst + +.. versionadded:: TBD + + loop over all computes + schedule next invocation for those that store invocation times + called when not sure what computes will be needed on newstep + do not loop only over n_timeflag, since may not be set yet + +*See also* + :cpp:func:`lammps_clearstep_compute` + :cpp:func:`lammps_addstep_compute` + +\endverbatim + + * \param handle pointer to a previously created LAMMPS instance cast to ``void *``. + * \param newstep pointer to bigint of next timestep the compute will be invoked + */ +void lammps_addstep_compute_all(void *handle, void *newstep) { + auto lmp = (LAMMPS *) handle; + auto ns = (bigint *) newstep; + if (lmp && lmp->modify && ns) lmp->modify->addstep_compute_all(*ns); +} +/* ---------------------------------------------------------------------- */ + +/** Add next timestep to compute if it has been invoked in the current timestep + * +\verbatim embed:rst + +.. versionadded:: TBD + + loop over computes that store invocation times + if its invoked flag set on this timestep, schedule next invocation + called everywhere that computes are used, after computes are invoked + +*See also* + :cpp:func:`lammps_addstep_compute_all` + :cpp:func:`lammps_clearstep_compute` + +\endverbatim + + * \param handle pointer to a previously created LAMMPS instance cast to ``void *``. + * \param newstep next timestep the compute will be invoked + */ +void lammps_addstep_compute(void *handle, void *newstep) { + auto lmp = (LAMMPS *) handle; + auto ns = (bigint *) newstep; + if (lmp && lmp->modify && ns) lmp->modify->addstep_compute(*ns); +} + // ---------------------------------------------------------------------- // Library functions for scatter/gather operations of data // ---------------------------------------------------------------------- diff --git a/src/library.h b/src/library.h index f64ae9e7c2..99b251ee85 100644 --- a/src/library.h +++ b/src/library.h @@ -191,6 +191,10 @@ int lammps_set_internal_variable(void *handle, const char *name, double value); int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); double lammps_eval(void *handle, const char *expr); +void lammps_clearstep_compute(void *handle); +void lammps_addstep_compute_all(void *handle, void * nextstep); +void lammps_addstep_compute(void *handle, void * nextstep); + /* ---------------------------------------------------------------------- * Library functions for scatter/gather operations of data * ---------------------------------------------------------------------- */ diff --git a/tools/swig/lammps.i b/tools/swig/lammps.i index bc9f46f407..476e6ad17d 100644 --- a/tools/swig/lammps.i +++ b/tools/swig/lammps.i @@ -143,6 +143,9 @@ extern int lammps_set_string_variable(void *, const char *, const char *); extern int lammps_set_internal_variable(void *, const char *, double); extern int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); extern double lammps_eval(void *handle, const char *expr); +extern void lammps_clearstep_compute(void *handle); +extern void lammps_addstep_compute(void *handle, void *nstep); +extern void lammps_addstep_compute_all(void *handle, void *nstep); extern void lammps_gather_atoms(void *, char *, int, int, void *); extern void lammps_gather_atoms_concat(void *, char *, int, int, void *); @@ -336,6 +339,9 @@ extern int lammps_set_string_variable(void *, const char *, const char *); extern int lammps_set_internal_variable(void *, const char *, double); extern int lammps_variable_info(void *handle, int idx, char *buf, int bufsize); extern double lammps_eval(void *handle, const char *expr); +extern void lammps_clearstep_compute(void *handle); +extern void lammps_addstep_compute(void *handle, void *nstep); +extern void lammps_addstep_compute_all(void *handle, void *nstep); extern void lammps_gather_atoms(void *, char *, int, int, void *); extern void lammps_gather_atoms_concat(void *, char *, int, int, void *); diff --git a/unittest/c-library/test_library_properties.cpp b/unittest/c-library/test_library_properties.cpp index 737015ccdc..1feccf5cb0 100644 --- a/unittest/c-library/test_library_properties.cpp +++ b/unittest/c-library/test_library_properties.cpp @@ -3,9 +3,12 @@ #include "library.h" #include "atom.h" +#include "compute.h" #include "lammps.h" #include "lmptype.h" +#include "modify.h" #include "platform.h" + #include #include @@ -668,6 +671,77 @@ TEST_F(LibraryProperties, neighlist) } }; +static constexpr char lj_setup[] = "lattice fcc 0.8442\n" + "region box block 0 10 0 10 0 10\n" + "create_box 1 box\n" + "create_atoms 1 box\n" + "mass 1 1.0\n" + "pair_style lj/cut 2.5\n" + "pair_coeff 1 1 1.0 1.0\n" + "fix 1 all nve\n"; + +TEST_F(LibraryProperties, step_compute) +{ + ::testing::internal::CaptureStdout(); + lammps_commands_string(lmp, lj_setup); + lammps_command(lmp, "compute pr all pressure thermo_temp"); + lammps_command(lmp, "fix av all ave/time 2 1 2 c_pr mode scalar"); + lammps_command(lmp, "run 2 post no"); + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + if (lammps_has_error(lmp)) { + char buf[2048]; + lammps_get_last_error_message(lmp, buf, 2048); + FAIL() << buf << "\n"; + } + auto lammps = (LAMMPS_NS::LAMMPS *)lmp; + auto icomp = lammps->modify->get_compute_by_id("pr"); + EXPECT_EQ(icomp->ntime, 2); + EXPECT_EQ(icomp->tlist[0], 4); + EXPECT_EQ(icomp->tlist[1], 2); + EXPECT_EQ(icomp->invoked_flag, 0); + EXPECT_EQ(icomp->invoked_scalar, 2); + EXPECT_EQ(icomp->invoked_vector, -1); + lammps_clearstep_compute(lmp); + EXPECT_EQ(icomp->invoked_flag, 0); + EXPECT_EQ(icomp->invoked_scalar, 2); + EXPECT_EQ(icomp->invoked_vector, -1); + bigint nextstep = 6; + lammps_addstep_compute(lmp, (void *)&nextstep); + EXPECT_EQ(icomp->ntime, 3); + EXPECT_EQ(icomp->tlist[0], 6); + EXPECT_EQ(icomp->tlist[1], 4); + EXPECT_EQ(icomp->tlist[2], 2); + EXPECT_EQ(icomp->invoked_flag, 0); + EXPECT_EQ(icomp->invoked_scalar, 2); + EXPECT_EQ(icomp->invoked_vector, -1); + lammps_command(lmp, "run 4 post no"); + EXPECT_EQ(icomp->ntime, 2); + EXPECT_EQ(icomp->tlist[0], 8); + EXPECT_EQ(icomp->tlist[1], 6); + EXPECT_EQ(icomp->invoked_flag, 0); + EXPECT_EQ(icomp->invoked_scalar, 6); + EXPECT_EQ(icomp->invoked_vector, -1); + lammps_command(lmp, "run 2 post no"); + EXPECT_EQ(icomp->ntime, 2); + EXPECT_EQ(icomp->tlist[0], 10); + EXPECT_EQ(icomp->tlist[1], 8); + EXPECT_EQ(icomp->invoked_flag, 0); + EXPECT_EQ(icomp->invoked_scalar, 8); + EXPECT_EQ(icomp->invoked_vector, -1); + nextstep = 9; + lammps_addstep_compute(lmp, (void *)&nextstep); + lammps_command(lmp, "run 1 post no"); + EXPECT_EQ(icomp->ntime, 2); + EXPECT_EQ(icomp->tlist[0], 10); + EXPECT_EQ(icomp->tlist[1], 9); + EXPECT_EQ(icomp->invoked_flag, 0); + EXPECT_EQ(icomp->invoked_scalar, -1); + EXPECT_EQ(icomp->invoked_vector, -1); + icomp->compute_scalar(); + EXPECT_EQ(icomp->invoked_scalar, 9); +} + TEST_F(LibraryProperties, has_error) { EXPECT_EQ(lammps_has_error(lmp), 0); diff --git a/unittest/force-styles/CMakeLists.txt b/unittest/force-styles/CMakeLists.txt index 043b628fd3..fb6cd1a63b 100644 --- a/unittest/force-styles/CMakeLists.txt +++ b/unittest/force-styles/CMakeLists.txt @@ -246,7 +246,7 @@ if(MLIAP_ENABLE_PYTHON AND (NOT WIN32)) add_executable(test_mliappy_unified test_mliappy_unified.cpp) target_link_libraries(test_mliappy_unified PRIVATE lammps GTest::GMockMain) add_test(NAME TestMliapPyUnified COMMAND test_mliappy_unified) - set_tests_properties(${TNAME} PROPERTIES ENVIRONMENT "${FORCE_TEST_ENVIRONMENT}") + set_tests_properties(TestMliapPyUnified PROPERTIES ENVIRONMENT "${FORCE_TEST_ENVIRONMENT}") endif() add_executable(test_pair_list test_pair_list.cpp) diff --git a/unittest/force-styles/test_angle_style.cpp b/unittest/force-styles/test_angle_style.cpp index ebd8170b7c..e24e5401df 100644 --- a/unittest/force-styles/test_angle_style.cpp +++ b/unittest/force-styles/test_angle_style.cpp @@ -53,15 +53,16 @@ using ::testing::StartsWith; using namespace LAMMPS_NS; -void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg) +void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg) { platform::unlink(cfg.basename + ".restart"); platform::unlink(cfg.basename + ".data"); platform::unlink(cfg.basename + "-coeffs.in"); delete lmp; + lmp = nullptr; } -LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true) +LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton) { LAMMPS *lmp; @@ -92,21 +93,7 @@ LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton // utility lambdas to improve readability auto command = [&](const std::string &line) { - try { - lmp->input->one(line); - } catch (LAMMPSAbortException &ae) { - fprintf(stderr, "LAMMPS Error: %s\n", ae.what()); - exit(2); - } catch (LAMMPSException &e) { - fprintf(stderr, "LAMMPS Error: %s\n", e.what()); - exit(3); - } catch (fmt::format_error &fe) { - fprintf(stderr, "fmt::format_error: %s\n", fe.what()); - exit(4); - } catch (std::exception &e) { - fprintf(stderr, "General exception: %s\n", e.what()); - exit(5); - } + lmp->input->one(line); }; auto parse_input_script = [&](const std::string &filename) { lmp->input->file(filename.c_str()); @@ -230,7 +217,12 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) // initialize system geometry LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite"}; - LAMMPS *lmp = init_lammps(args, config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, config, true); + } catch (std::exception &e) { + FAIL() << e.what(); + } if (!lmp) { std::cerr << "One or more prerequisite styles are not available " "in this LAMMPS configuration:\n"; @@ -321,8 +313,14 @@ TEST(AngleStyle, plain) LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -371,7 +369,12 @@ TEST(AngleStyle, plain) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -439,8 +442,14 @@ TEST(AngleStyle, omp) "-pk", "omp", "4", "-sf", "omp"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -493,7 +502,12 @@ TEST(AngleStyle, omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -540,14 +554,21 @@ TEST(AngleStyle, kokkos_omp) // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only. if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") || Info::has_accelerator_feature("KOKKOS", "api", "hip") || - Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP(); + Info::has_accelerator_feature("KOKKOS", "api", "sycl")) + GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled"; LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite", "-k", "on", "t", "4", "-sf", "kk"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -597,7 +618,12 @@ TEST(AngleStyle, kokkos_omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -664,8 +690,14 @@ TEST(AngleStyle, numdiff) LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -717,7 +749,13 @@ TEST(AngleStyle, single) // create a LAMMPS instance with standard settings to detect the number of atom types if (!verbose) ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); if (!lmp) { @@ -865,7 +903,13 @@ TEST(AngleStyle, extract) LAMMPS::argv args = {"AngleStyle", "-log", "none", "-echo", "screen", "-nocite"}; if (!verbose) ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); if (!lmp) { diff --git a/unittest/force-styles/test_bond_style.cpp b/unittest/force-styles/test_bond_style.cpp index 4a3a985ffc..660435cf49 100644 --- a/unittest/force-styles/test_bond_style.cpp +++ b/unittest/force-styles/test_bond_style.cpp @@ -53,15 +53,16 @@ using ::testing::StartsWith; using namespace LAMMPS_NS; -void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg) +void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg) { platform::unlink(cfg.basename + ".restart"); platform::unlink(cfg.basename + ".data"); platform::unlink(cfg.basename + "-coeffs.in"); delete lmp; + lmp = nullptr; } -LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true) +LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton) { LAMMPS *lmp; @@ -92,21 +93,7 @@ LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton // utility lambdas to improve readability auto command = [&](const std::string &line) { - try { - lmp->input->one(line); - } catch (LAMMPSAbortException &ae) { - fprintf(stderr, "LAMMPS Error: %s\n", ae.what()); - exit(2); - } catch (LAMMPSException &e) { - fprintf(stderr, "LAMMPS Error: %s\n", e.what()); - exit(3); - } catch (fmt::format_error &fe) { - fprintf(stderr, "fmt::format_error: %s\n", fe.what()); - exit(4); - } catch (std::exception &e) { - fprintf(stderr, "General exception: %s\n", e.what()); - exit(5); - } + lmp->input->one(line); }; auto parse_input_script = [&](const std::string &filename) { lmp->input->file(filename.c_str()); @@ -230,7 +217,13 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) // initialize system geometry LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite"}; - LAMMPS *lmp = init_lammps(args, config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, config, true); + } catch (std::exception &e) { + FAIL() << e.what(); + } + if (!lmp) { std::cerr << "One or more prerequisite styles are not available " "in this LAMMPS configuration:\n"; @@ -321,8 +314,14 @@ TEST(BondStyle, plain) LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -371,7 +370,12 @@ TEST(BondStyle, plain) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -441,8 +445,14 @@ TEST(BondStyle, omp) "-pk", "omp", "4", "-sf", "omp"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -495,7 +505,12 @@ TEST(BondStyle, omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -542,14 +557,21 @@ TEST(BondStyle, kokkos_omp) // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only. if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") || Info::has_accelerator_feature("KOKKOS", "api", "hip") || - Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP(); + Info::has_accelerator_feature("KOKKOS", "api", "sycl")) + GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled"; LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite", "-k", "on", "t", "4", "-sf", "kk"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -603,7 +625,12 @@ TEST(BondStyle, kokkos_omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -652,8 +679,14 @@ TEST(BondStyle, numdiff) LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -705,7 +738,13 @@ TEST(BondStyle, single) // create a LAMMPS instance with standard settings to detect the number of atom types if (!verbose) ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); if (!lmp) { @@ -959,7 +998,13 @@ TEST(BondStyle, extract) LAMMPS::argv args = {"BondStyle", "-log", "none", "-echo", "screen", "-nocite"}; if (!verbose) ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); if (!lmp) { diff --git a/unittest/force-styles/test_dihedral_style.cpp b/unittest/force-styles/test_dihedral_style.cpp index b538c45f42..5d0bd86d2b 100644 --- a/unittest/force-styles/test_dihedral_style.cpp +++ b/unittest/force-styles/test_dihedral_style.cpp @@ -53,15 +53,16 @@ using ::testing::StartsWith; using namespace LAMMPS_NS; -void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg) +void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg) { platform::unlink(cfg.basename + ".restart"); platform::unlink(cfg.basename + ".data"); platform::unlink(cfg.basename + "-coeffs.in"); delete lmp; + lmp = nullptr; } -LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true) +LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton) { auto *lmp = new LAMMPS(args, MPI_COMM_WORLD); @@ -237,7 +238,13 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) // initialize system geometry LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite"}; - LAMMPS *lmp = init_lammps(args, config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, config, true); + } catch (std::exception &e) { + FAIL() << e.what(); + } + if (!lmp) { std::cerr << "One or more prerequisite styles are not available " "in this LAMMPS configuration:\n"; @@ -322,8 +329,14 @@ TEST(DihedralStyle, plain) LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -372,7 +385,12 @@ TEST(DihedralStyle, plain) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -442,8 +460,14 @@ TEST(DihedralStyle, omp) "-pk", "omp", "4", "-sf", "omp"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -497,7 +521,12 @@ TEST(DihedralStyle, omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -544,15 +573,22 @@ TEST(DihedralStyle, kokkos_omp) // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only. if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") || Info::has_accelerator_feature("KOKKOS", "api", "hip") || - Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP(); + Info::has_accelerator_feature("KOKKOS", "api", "sycl")) + GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled"; LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite", "-k", "on", "t", "4", "-sf", "kk"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -606,7 +642,12 @@ TEST(DihedralStyle, kokkos_omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -655,8 +696,14 @@ TEST(DihedralStyle, numdiff) LAMMPS::argv args = {"DihedralStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; diff --git a/unittest/force-styles/test_fix_timestep.cpp b/unittest/force-styles/test_fix_timestep.cpp index 957226d22b..b952a3c045 100644 --- a/unittest/force-styles/test_fix_timestep.cpp +++ b/unittest/force-styles/test_fix_timestep.cpp @@ -56,13 +56,14 @@ using ::testing::StartsWith; using namespace LAMMPS_NS; -void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg) +void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg) { platform::unlink(cfg.basename + ".restart"); delete lmp; + lmp = nullptr; } -LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool use_respa = false) +LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool use_respa) { LAMMPS *lmp; @@ -178,7 +179,12 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) { // initialize system geometry LAMMPS::argv args = {"FixIntegrate", "-log", "none", "-echo", "screen", "-nocite"}; - LAMMPS *lmp = init_lammps(args, config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, config, false); + } catch (std::exception &e) { + FAIL() << e.what(); + } if (!lmp) { std::cerr << "One or more prerequisite styles are not available " "in this LAMMPS configuration:\n"; @@ -272,7 +278,14 @@ TEST(FixTimestep, plain) LAMMPS::argv args = {"FixTimestep", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -430,14 +443,20 @@ TEST(FixTimestep, plain) // fix nve/limit cannot work with r-RESPA ifix = lmp->modify->get_fix_by_id("test"); if (ifix && !utils::strmatch(ifix->style, "^rigid") && - !utils::strmatch(ifix->style, "^nve/limit") && - !utils::strmatch(ifix->style, "^recenter")) { + !utils::strmatch(ifix->style, "^nve/limit") && !utils::strmatch(ifix->style, "^recenter")) { if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); if (!verbose) ::testing::internal::GetCapturedStdout(); ::testing::internal::CaptureStdout(); - lmp = init_lammps(args, test_config, true); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -572,15 +591,26 @@ TEST(FixTimestep, omp) "-pk", "omp", "4", "-sf", "omp"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; if (!lmp) { - std::cerr << "One or more prerequisite styles are not available " + std::cerr << "One or more prerequisite styles with /omp suffix are not available " "in this LAMMPS configuration:\n"; for (auto &prerequisite : test_config.prerequisites) { - std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n"; + if (prerequisite.first == "atom") { + std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n"; + } else { + std::cerr << prerequisite.first << "_style " << prerequisite.second << "/omp\n"; + } } GTEST_SKIP(); } @@ -731,7 +761,13 @@ TEST(FixTimestep, omp) if (!verbose) ::testing::internal::GetCapturedStdout(); ::testing::internal::CaptureStdout(); - lmp = init_lammps(args, test_config, true); + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -861,13 +897,21 @@ TEST(FixTimestep, kokkos_omp) // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only. if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") || Info::has_accelerator_feature("KOKKOS", "api", "hip") || - Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP(); - + Info::has_accelerator_feature("KOKKOS", "api", "sycl")) { + GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled"; + } LAMMPS::argv args = {"FixTimestep", "-log", "none", "-echo", "screen", "-nocite", "-k", "on", "t", "4", "-sf", "kk"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -875,7 +919,7 @@ TEST(FixTimestep, kokkos_omp) std::cerr << "One or more prerequisite styles with /kk suffix\n" "are not available in this LAMMPS configuration:\n"; for (auto &prerequisite : test_config.prerequisites) { - std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n"; + std::cerr << prerequisite.first << "_style " << prerequisite.second << "/kk\n"; } GTEST_SKIP(); } diff --git a/unittest/force-styles/test_improper_style.cpp b/unittest/force-styles/test_improper_style.cpp index cddb27a2f7..e272033e85 100644 --- a/unittest/force-styles/test_improper_style.cpp +++ b/unittest/force-styles/test_improper_style.cpp @@ -53,15 +53,16 @@ using ::testing::StartsWith; using namespace LAMMPS_NS; -void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg) +void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg) { platform::unlink(cfg.basename + ".restart"); platform::unlink(cfg.basename + ".data"); platform::unlink(cfg.basename + "-coeffs.in"); delete lmp; + lmp = nullptr; } -LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true) +LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton) { LAMMPS *lmp; @@ -92,21 +93,7 @@ LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton // utility lambdas to improve readability auto command = [&](const std::string &line) { - try { - lmp->input->one(line); - } catch (LAMMPSAbortException &ae) { - fprintf(stderr, "LAMMPS Error: %s\n", ae.what()); - exit(2); - } catch (LAMMPSException &e) { - fprintf(stderr, "LAMMPS Error: %s\n", e.what()); - exit(3); - } catch (fmt::format_error &fe) { - fprintf(stderr, "fmt::format_error: %s\n", fe.what()); - exit(4); - } catch (std::exception &e) { - fprintf(stderr, "General exception: %s\n", e.what()); - exit(5); - } + lmp->input->one(line); }; auto parse_input_script = [&](const std::string &filename) { lmp->input->file(filename.c_str()); @@ -230,7 +217,12 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) // initialize system geometry LAMMPS::argv args = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite"}; - LAMMPS *lmp = init_lammps(args, config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, config, true); + } catch (std::exception &e) { + FAIL() << e.what(); + } if (!lmp) { std::cerr << "One or more prerequisite styles are not available " "in this LAMMPS configuration:\n"; @@ -315,8 +307,14 @@ TEST(ImproperStyle, plain) LAMMPS::argv args = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -365,7 +363,12 @@ TEST(ImproperStyle, plain) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -435,8 +438,14 @@ TEST(ImproperStyle, omp) "-pk", "omp", "4", "-sf", "omp"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -490,7 +499,12 @@ TEST(ImproperStyle, omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -536,15 +550,22 @@ TEST(ImproperStyle, kokkos_omp) // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only. if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") || Info::has_accelerator_feature("KOKKOS", "api", "hip") || - Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP(); + Info::has_accelerator_feature("KOKKOS", "api", "sycl")) + GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled"; LAMMPS::argv args = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite", "-k", "on", "t", "4", "-sf", "kk"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -597,7 +618,12 @@ TEST(ImproperStyle, kokkos_omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton bond is forced to be on @@ -644,8 +670,14 @@ TEST(ImproperStyle, numdiff) LAMMPS::argv args = {"ImproperStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; diff --git a/unittest/force-styles/test_pair_style.cpp b/unittest/force-styles/test_pair_style.cpp index bb1346fe70..6dabeab7a4 100644 --- a/unittest/force-styles/test_pair_style.cpp +++ b/unittest/force-styles/test_pair_style.cpp @@ -54,15 +54,16 @@ using ::testing::StartsWith; using namespace LAMMPS_NS; -void cleanup_lammps(LAMMPS *lmp, const TestConfig &cfg) +void cleanup_lammps(LAMMPS *&lmp, const TestConfig &cfg) { platform::unlink(cfg.basename + ".restart"); platform::unlink(cfg.basename + ".data"); platform::unlink(cfg.basename + "-coeffs.in"); delete lmp; + lmp = nullptr; } -LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton = true) +LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton) { LAMMPS *lmp; @@ -93,21 +94,7 @@ LAMMPS *init_lammps(LAMMPS::argv &args, const TestConfig &cfg, const bool newton // utility lambdas to improve readability auto command = [&](const std::string &line) { - try { - lmp->input->one(line); - } catch (LAMMPSAbortException &ae) { - fprintf(stderr, "LAMMPS Error: %s\n", ae.what()); - exit(2); - } catch (LAMMPSException &e) { - fprintf(stderr, "LAMMPS Error: %s\n", e.what()); - exit(3); - } catch (fmt::format_error &fe) { - fprintf(stderr, "fmt::format_error: %s\n", fe.what()); - exit(4); - } catch (std::exception &e) { - fprintf(stderr, "General exception: %s\n", e.what()); - exit(5); - } + lmp->input->one(line); }; auto parse_input_script = [&](const std::string &filename) { @@ -242,8 +229,12 @@ void generate_yaml_file(const char *outfile, const TestConfig &config) { // initialize system geometry LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"}; - - LAMMPS *lmp = init_lammps(args, config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, config, true); + } catch (std::exception &e) { + FAIL() << e.what(); + } if (!lmp) { std::cerr << "One or more prerequisite styles are not available " "in this LAMMPS configuration:\n"; @@ -340,8 +331,14 @@ TEST(PairStyle, plain) LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -399,7 +396,12 @@ TEST(PairStyle, plain) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // skip over these tests if newton pair is forced to be on @@ -480,9 +482,14 @@ TEST(PairStyle, plain) if (pair->respa_enable) { if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); - lmp->input->one("run_style respa 2 1 inner 1 4.8 5.5 outer 2"); - run_lammps(lmp); + try { + lmp = init_lammps(args, test_config, false); + lmp->input->one("run_style respa 2 1 inner 1 4.8 5.5 outer 2"); + run_lammps(lmp); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); // need to relax error by a large amount with tabulation, since @@ -519,8 +526,14 @@ TEST(PairStyle, omp) if (utils::strmatch(test_config.pair_style, "^dpd")) args[8] = "1"; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -578,7 +591,12 @@ TEST(PairStyle, omp) if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); pair = lmp->force->pair; @@ -636,7 +654,9 @@ TEST(PairStyle, kokkos_omp) // if KOKKOS has GPU support enabled, it *must* be used. We cannot test OpenMP only. if (Info::has_accelerator_feature("KOKKOS", "api", "cuda") || Info::has_accelerator_feature("KOKKOS", "api", "hip") || - Info::has_accelerator_feature("KOKKOS", "api", "sycl")) GTEST_SKIP(); + Info::has_accelerator_feature("KOKKOS", "api", "sycl")) { + GTEST_SKIP() << "Cannot test KOKKOS/OpenMP with GPU support enabled"; + } LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-k", "on", "t", "4", "-sf", "kk"}; @@ -655,8 +675,14 @@ TEST(PairStyle, kokkos_omp) args[9] = "1"; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -713,7 +739,12 @@ TEST(PairStyle, kokkos_omp) if (lmp->force->newton_pair == 0) { if (!verbose) ::testing::internal::CaptureStdout(); cleanup_lammps(lmp, test_config); - lmp = init_lammps(args, test_config, false); + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); pair = lmp->force->pair; @@ -788,8 +819,14 @@ TEST(PairStyle, gpu) } ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, false); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, false); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -868,8 +905,14 @@ TEST(PairStyle, intel) if (utils::strmatch(test_config.pair_style, "^dpd")) args[12] = "1"; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -948,8 +991,14 @@ TEST(PairStyle, opt) LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite", "-sf", "opt"}; ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config); - + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + std::string output = ::testing::internal::GetCapturedStdout(); + if (verbose) std::cout << output; + FAIL() << e.what(); + } std::string output = ::testing::internal::GetCapturedStdout(); if (verbose) std::cout << output; @@ -1032,7 +1081,13 @@ TEST(PairStyle, single) // create a LAMMPS instance with standard settings to detect the number of atom types if (!verbose) ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); if (!lmp) { @@ -1276,7 +1331,13 @@ TEST(PairStyle, extract) LAMMPS::argv args = {"PairStyle", "-log", "none", "-echo", "screen", "-nocite"}; if (!verbose) ::testing::internal::CaptureStdout(); - LAMMPS *lmp = init_lammps(args, test_config, true); + LAMMPS *lmp = nullptr; + try { + lmp = init_lammps(args, test_config, true); + } catch (std::exception &e) { + if (!verbose) ::testing::internal::GetCapturedStdout(); + FAIL() << e.what(); + } if (!verbose) ::testing::internal::GetCapturedStdout(); if (!lmp) { diff --git a/unittest/force-styles/tests/fix-timestep-efield_lepton.yaml b/unittest/force-styles/tests/fix-timestep-efield_lepton.yaml index c6f80af274..7928c81c45 100644 --- a/unittest/force-styles/tests/fix-timestep-efield_lepton.yaml +++ b/unittest/force-styles/tests/fix-timestep-efield_lepton.yaml @@ -6,7 +6,7 @@ epsilon: 2e-13 skip_tests: prerequisites: ! | atom full - fix efield + fix efield/lepton pre_commands: ! "" post_commands: ! | region half block 0 EDGE EDGE EDGE EDGE EDGE diff --git a/unittest/force-styles/tests/fix-timestep-efield_lepton_region.yaml b/unittest/force-styles/tests/fix-timestep-efield_lepton_region.yaml index b2275f047b..f0a6b121a4 100644 --- a/unittest/force-styles/tests/fix-timestep-efield_lepton_region.yaml +++ b/unittest/force-styles/tests/fix-timestep-efield_lepton_region.yaml @@ -6,7 +6,7 @@ epsilon: 2e-13 skip_tests: prerequisites: ! | atom full - fix efield + fix efield/lepton pre_commands: ! "" post_commands: ! | region half block 0 EDGE EDGE EDGE EDGE EDGE