ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'lammps:develop' into bond_react_fixes_aug23

Browse Source
This commit is contained in:
Jacob Gissinger
2024-01-03 14:59:26 -05:00
committed by GitHub
parent 4d2aed8937 35214120ad
commit d9ea10ca24
13 changed files with 198 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

41
doc/src/Build_development.rst
View File

@ -255,16 +255,18 @@ A test run is then a a collection multiple individual test runs each
with many comparisons to reference results based on template input
files, individual command settings, relative error margins, and
reference data stored in a YAML format file with ``.yaml``
suffix. Currently the programs ``test_pair_style``, ``test_bond_style``, and
``test_angle_style`` are implemented. They will compare forces, energies and
(global) stress for all atoms after a ``run 0`` calculation and after a
few steps of MD with :doc:`fix nve <fix_nve>`, each in multiple variants
with different settings and also for multiple accelerated styles. If a
prerequisite style or package is missing, the individual tests are
skipped. All tests will be executed on a single MPI process, so using
the CMake option ``-D BUILD_MPI=off`` can significantly speed up testing,
since this will skip the MPI initialization for each test run.
Below is an example command and output:
suffix. Currently the programs ``test_pair_style``, ``test_bond_style``,
``test_angle_style``, ``test_dihedral_style``, and
``test_improper_style`` are implemented. They will compare forces,
energies and (global) stress for all atoms after a ``run 0`` calculation
and after a few steps of MD with :doc:`fix nve <fix_nve>`, each in
multiple variants with different settings and also for multiple
accelerated styles. If a prerequisite style or package is missing, the
individual tests are skipped. All force style tests will be executed on
a single MPI process, so using the CMake option ``-D BUILD_MPI=off`` can
significantly speed up testing, since this will skip the MPI
initialization for each test run. Below is an example command and
output:
.. code-block:: console
@ -416,15 +418,16 @@ When compiling LAMMPS with enabled tests, most test executables will
need to be linked against the LAMMPS library. Since this can be a very
large library with many C++ objects when many packages are enabled, link
times can become very long on machines that use the GNU BFD linker (e.g.
Linux systems). Alternatives like the ``lld`` linker of the LLVM project
or the ``gold`` linker available with GNU binutils can speed up this step
substantially. CMake will by default test if any of the two can be
enabled and use it when ``ENABLE_TESTING`` is active. It can also be
selected manually through the ``CMAKE_CUSTOM_LINKER`` CMake variable.
Allowed values are ``lld``, ``gold``, ``bfd``, or ``default``. The
``default`` option will use the system default linker otherwise, the
linker is chosen explicitly. This option is only available for the
GNU or Clang C++ compiler.
Linux systems). Alternatives like the ``mold`` linker, the ``lld``
linker of the LLVM project, or the ``gold`` linker available with GNU
binutils can speed up this step substantially (in this order). CMake
will by default test if any of the three can be enabled and use it when
``ENABLE_TESTING`` is active. It can also be selected manually through
the ``CMAKE_CUSTOM_LINKER`` CMake variable. Allowed values are
``mold``, ``lld``, ``gold``, ``bfd``, or ``default``. The ``default``
option will use the system default linker otherwise, the linker is
chosen explicitly. This option is only available for the GNU or Clang
C++ compilers.
Tests for other components and utility functions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

101
doc/src/Developer_unittest.rst
View File

@ -121,7 +121,7 @@ will be suppressed and only a summary printed, but adding
the '-V' option will then produce output from the tests
above like the following:
.. code-block::
.. code-block:: console
[...]
1: [ RUN ] Tokenizer.empty_string
@ -526,3 +526,102 @@ The ``unittest/tools`` folder contains tests for programs in the
shell, which are implemented as a python scripts using the ``unittest``
Python module and launching the tool commands through the ``subprocess``
Python module.
Troubleshooting failed unit tests
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The are by default no unit tests for newly added features (e.g. pair, fix,
or compute styles) unless your pull request also includes tests for the
added features. If you are modifying some features, you may see failures
for existing tests, if your modifications have some unexpected side effects
or your changes render the existing text invalid. If you are adding an
accelerated version of an existing style, then only tests for INTEL,
KOKKOS (with OpenMP only), OPENMP, and OPT will be run automatically.
Tests for the GPU package are time consuming and thus are only run
*after* a merge, or when a special label, ``gpu_unit_tests`` is added
to the pull request. After the test has started, it is often best to
remove the label since every PR activity will re-trigger the test (that
is a limitation of triggering a test with a label). Support for unit
tests with using KOKKOS with GPU acceleration is currently not supported.
When you see a failed build on GitHub, click on ``Details`` to be taken
to the corresponding LAMMPS Jenkins CI web page. Click on the "Exit"
symbol near the ``Logout`` button on the top right of that page to go to
the "classic view". In the classic view, there is a list of the
individual runs that make up this test run (they are shown but cannot be
inspected in the default view). You can click on any of those.
Clicking on ``Test Result`` will display the list of failed tests. Click
on the "Status" column to sort the tests based on their Failed or Passed
status. Then click on the failed test to expand its output.
For example, the following output snippet shows the failed unit test
.. code-block:: console
[ RUN ] PairStyle.gpu
/home/builder/workspace/dev/pull_requests/ubuntu_gpu/unit_tests/cmake_gpu_opencl_mixed_smallbig_clang_static/unittest/force-styles/test_main.cpp:63: Failure
Expected: (err) <= (epsilon)
Actual: 0.00018957912910606503 vs 0.0001
Google Test trace:
/home/builder/workspace/dev/pull_requests/ubuntu_gpu/unit_tests/cmake_gpu_opencl_mixed_smallbig_clang_static/unittest/force-styles/test_main.cpp:56: EXPECT_FORCES: init_forces (newton off)
/home/builder/workspace/dev/pull_requests/ubuntu_gpu/unit_tests/cmake_gpu_opencl_mixed_smallbig_clang_static/unittest/force-styles/test_main.cpp:64: Failure
Expected: (err) <= (epsilon)
Actual: 0.00022892713393549854 vs 0.0001
The failed assertions provide line numbers in the test source
(e.g. ``test_main.cpp:56``), from which one can understand what
specific assertion failed.
Note that the force style engine runs one of a small number of systems
in a rather off-equilibrium configuration with a few atoms for a few
steps, writes data and restart files, uses :doc:`the clear command
<clear>` to reset LAMMPS, and then runs from those files with different
settings (e.g. newton on/off) and integrators (e.g. verlet vs. respa).
Beyond potential issues/bugs in the source code, the mismatch between
the expected and actual values could be that force arrays are not
properly cleared between multiple run commands or that class members are
not correctly initialized or written to or read from a data or restart
file.
While the epsilon (relative precision) for a single, `IEEE 754 compliant
<https://en.wikipedia.org/wiki/IEEE_754>`_, double precision floating
point operation is at about 2.2e-16, the achievable precision for the
tests is lower due to most numbers being sums over intermediate results
and the non-associativity of floating point math leading to larger
errors. In some cases specific properties of the tested style. As a
rule of thumb, the test epsilon can often be in the range 5.0e-14 to
1.0e-13. But for "noisy" force kernels, e.g. those a larger amount of
arithmetic operations involving `exp()`, `log()` or `sin()` functions,
and also due to the effect of compiler optimization or differences
between compilers or platforms, epsilon may need to be further relaxed,
sometimes epsilon can be relaxed to 1.0e-12. If interpolation or lookup
tables are used, epsilon may need to be set to 1.0e-10 or even higher.
For tests of accelerated styles, the per-test epsilon is multiplied
by empirical factors that take into account the differences in the order
of floating point operations or that some or most intermediate operations
may be done using approximations or with single precision floating point
math.
To rerun the failed unit test individually, change to the ``build`` directory
and run the test with verbose output. For example,
.. code-block:: bash
env TEST_ARGS=-v ctest -R ^MolPairStyle:lj_cut_coul_long -V
``ctest`` with the ``-V`` flag also shows the exact command line
of the test. One can then use ``gdb --args`` to further debug and
catch exceptions with the test command, for example,
.. code-block:: bash
gdb --args /path/to/lammps/build/test_pair_style /path/to/lammps/unittest/force-styles/tests/mol-pair-lj_cut_coul_long.yaml
It is recommended to configure the build with ``-D
BUILD_SHARED_LIBS=on`` and use a custom linker to shorten the build time
during recompilation. Installing `ccache` in your development
environment helps speed up recompilation by caching previous
compilations and detecting when the same compilation is being done
again. Please see :doc:`Build_development` for further details.

20
doc/src/pair_ilp_tmd.rst
View File

@ -22,12 +22,12 @@ Examples
.. code-block:: LAMMPS
pair_style hybrid/overlay ilp/tmd 16.0 1
pair_coeff * * ilp/tmd MoS2.ILP Mo S S
pair_coeff * * ilp/tmd TMD.ILP Mo S S
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
pair_coeff * * ilp/tmd MoS2.ILP Mo S S Mo S S
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL W Se Se
pair_coeff * * ilp/tmd TMD.ILP Mo S S W Se Se
Description
"""""""""""
@ -36,7 +36,7 @@ Description
The *ilp/tmd* style computes the registry-dependent interlayer
potential (ILP) potential for transition metal dichalcogenides (TMD)
as described in :ref:`(Ouyang7) <Ouyang7>`.
as described in :ref:`(Ouyang7) <Ouyang7>` and :ref:`(Jiang) <Jiang>`.
.. math::
@ -69,7 +69,7 @@ calculating the normals.
each atom `i`, its six nearest neighboring atoms belonging to the same
sub-layer are chosen to define the normal vector `{\bf n}_i`.
The parameter file (e.g. MoS2.ILP), is intended for use with *metal*
The parameter file (e.g. TMD.ILP), is intended for use with *metal*
:doc:`units <units>`, with energies in meV. Two additional parameters,
*S*, and *rcut* are included in the parameter file. *S* is designed to
facilitate scaling of energies. *rcut* is designed to build the neighbor
@ -77,7 +77,7 @@ list for calculating the normals for each atom pair.
.. note::
The parameters presented in the parameter file (e.g. MoS2.ILP),
The parameters presented in the parameter file (e.g. TMD.ILP),
are fitted with taper function by setting the cutoff equal to 16.0
Angstrom. Using different cutoff or taper function should be careful.
These parameters provide a good description in both short- and long-range
@ -133,10 +133,10 @@ if LAMMPS was built with that package. See the :doc:`Build package
This pair style requires the newton setting to be *on* for pair
interactions.
The MoS2.ILP potential file provided with LAMMPS (see the potentials
The TMD.ILP potential file provided with LAMMPS (see the potentials
directory) are parameterized for *metal* units. You can use this
potential with any LAMMPS units, but you would need to create your own
custom MoS2.ILP potential file with coefficients listed in the appropriate
custom TMD.ILP potential file with coefficients listed in the appropriate
units, if your simulation does not use *metal* units.
Related commands
@ -164,3 +164,7 @@ tap_flag = 1
.. _Ouyang7:
**(Ouyang7)** W. Ouyang, et al., J. Chem. Theory Comput. 17, 7237 (2021).
.. _Jiang:
**(Jiang)** W. Jiang, et al., J. Phys. Chem. A, 127, 46, 9820–9830 (2023).

1
examples/PACKAGES/interlayer/ilp_tmds/MoS2.ILP
View File

@ -1 +0,0 @@
../../../../potentials/MoS2.ILP

1
examples/PACKAGES/interlayer/ilp_tmds/TMD.ILP Symbolic link
View File

@ -0,0 +1 @@
../../../../potentials/TMD.ILP

2
examples/PACKAGES/interlayer/ilp_tmds/in.mos2
View File

@ -12,7 +12,7 @@ mass 4 95.94
pair_style hybrid/overlay sw/mod sw/mod ilp/tmd 16.0
pair_coeff * * sw/mod 1 tmd.sw.mod Mo S S NULL NULL NULL
pair_coeff * * sw/mod 2 tmd.sw.mod NULL NULL NULL Mo S S
pair_coeff * * ilp/tmd MoS2.ILP Mo S S Mo S S
pair_coeff * * ilp/tmd TMD.ILP Mo S S Mo S S
# Calculate the pair potential
compute 0 all pair ilp/tmd

25
potentials/TMD.ILP Normal file
View File

@ -0,0 +1,25 @@
# DATE: 2021-12-02 UNITS: metal CONTRIBUTOR: Wengen Ouyang w.g.ouyang@gmail.com
# CITATION: W. Ouyang, et al., J. Chem. Theory Comput. 17, 7237 (2021).
# CITATION: W. Jiang, et al., J. Phys. Chem. A, 127, 46, 9820–9830 (2023).
# Interlayer Potential (ILP) for bilayer and bulk Group-VI Transition Metal Dichalcogenides.
# The parameters below are fitted against the HSE + MBD-NL DFT reference data.
#
# -------------------- Repulsion Potential -------------------++++++++++++++++ Vdw Potential ++++++++++++++++*********
# beta(A) alpha delta(A) epsilon(meV) C(meV) d sR reff(A) C6(meV*A^6) S rcut
Mo Mo 5.579450 9.377662 2.027222 144.151775 97.978570 89.437597 2.059031 5.122055 491850.316195 1.0 4.0
W W 5.530854 6.624992 1.983208 0.271792 140.174059 107.392585 1.356333 4.437591 691850.243962 1.0 4.0
S S 3.161402 8.093263 1.953140 4.586764 118.065466 58.809416 0.215367 4.299600 148811.243409 1.0 4.0
Se Se 3.938627 10.515924 2.415783 3.012583 22.400612 116.864517 0.151121 5.884241 112506.195626 1.0 4.0
Mo W 5.412298 8.647128 2.108665 51.177950 184.342860 201.281256 2.547743 2.492287 99996.913401 1.0 4.0
Mo S 3.627152 19.971375 7.585031 76.101931 3.317496 45.720328 0.947470 4.410425 150597.857716 1.0 4.0
Mo Se 6.196447 4.844134 14.362005 7.407221 0.058823 27.156223 0.976771 3.979186 786029.840651 1.0 4.0
W S 3.680136 11.163004 32.254117 110.019679 79.381335 138.340438 0.900750 8.875776 250600.809034 1.0 4.0
W Se 3.559392 20.638856 1.202717 20.478669 197.422484 10.005271 1.052738 3.815817 288321.561114 1.0 4.0
S Se 2.820092 7.491151 1.933323 141.532559 293.127817 90.470904 0.390492 4.170885 117688.987069 1.0 4.0
# Symmetric Atom Pair
W Mo 5.412298 8.647128 2.108665 51.177950 184.342860 201.281256 2.547743 2.492287 99996.913401 1.0 4.0
S Mo 3.627152 19.971375 7.585031 76.101931 3.317496 45.720328 0.947470 4.410425 150597.857716 1.0 4.0
Se Mo 6.196447 4.844134 14.362005 7.407221 0.058823 27.156223 0.976771 3.979186 786029.840651 1.0 4.0
S W 3.680136 11.163004 32.254117 110.019679 79.381335 138.340438 0.900750 8.875776 250600.809034 1.0 4.0
Se W 3.559392 20.638856 1.202717 20.478669 197.422484 10.005271 1.052738 3.815817 288321.561114 1.0 4.0
Se S 2.820092 7.491151 1.933323 141.532559 293.127817 90.470904 0.390492 4.170885 117688.987069 1.0 4.0

2
src/INTERLAYER/pair_ilp_tmd.cpp
View File

@ -210,7 +210,7 @@ void PairILPTMD::calc_FRep(int eflag, int /* vflag */)
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
ev_tally_xyz(k, i, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}

2
src/INTERLAYER/pair_kolmogorov_crespi_full.cpp
View File

@ -590,7 +590,7 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */)
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
ev_tally_xyz(k, i, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}

59
src/MOLECULE/fix_cmap.cpp
View File

@ -39,6 +39,7 @@
#include "memory.h"
#include "potential_file_reader.h"
#include "respa.h"
#include "text_file_reader.h"
#include "update.h"
#include <cmath>
@ -49,15 +50,14 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
#define MAXLINE 256
#define LISTDELTA 10000
#define LB_FACTOR 1.5
static constexpr int LISTDELTA = 10000;
static constexpr double LB_FACTOR = 1.5;
#define CMAPMAX 6 // max # of CMAP terms stored by one atom
#define CMAPDIM 24 // grid map dimension is 24 x 24
#define CMAPXMIN -360.0
#define CMAPXMIN2 -180.0
#define CMAPDX 15.0 // 360/CMAPDIM
static constexpr int CMAPMAX = 6; // max # of CMAP terms stored by one atom
static constexpr int CMAPDIM = 24; // grid map dimension is 24 x 24
static constexpr double CMAPXMIN = -360.0;
static constexpr double CMAPXMIN2 = -180.0;
static constexpr double CMAPDX = 15.0; // 360/CMAPDIM
/* ---------------------------------------------------------------------- */
@ -86,17 +86,15 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) :
wd_section = 1;
respa_level_support = 1;
ilevel_respa = 0;
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
eflag_caller = 1;
// allocate memory for CMAP data
memory->create(g_axis,CMAPDIM,"cmap:g_axis");
memory->create(cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:grid");
memory->create(d1cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d1grid");
memory->create(d2cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d2grid");
memory->create(d12cmapgrid,6,CMAPDIM,CMAPDIM,"cmap:d12grid");
memory->create(cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:grid");
memory->create(d1cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:d1grid");
memory->create(d2cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:d2grid");
memory->create(d12cmapgrid,CMAPMAX,CMAPDIM,CMAPDIM,"cmap:d12grid");
// read and setup CMAP data
@ -184,10 +182,6 @@ void FixCMAP::init()
for (i = 0; i < 6; i++)
set_map_derivatives(cmapgrid[i],d1cmapgrid[i],d2cmapgrid[i],d12cmapgrid[i]);
// define newton_bond here in case restart file was read (not data file)
newton_bond = force->newton_bond;
if (utils::strmatch(update->integrate_style,"^respa")) {
ilevel_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels-1;
if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
@ -238,6 +232,8 @@ void FixCMAP::min_setup(int vflag)
void FixCMAP::pre_neighbor()
{
int i,m,atom1,atom2,atom3,atom4,atom5;
const int me = comm->me;
const int nprocs = comm->nprocs;
// guesstimate initial length of local crossterm list
// if ncmap was not set (due to read_restart, no read_data),
@ -637,15 +633,22 @@ void FixCMAP::read_grid_map(char *cmapfile)
{
if (comm->me == 0) {
try {
memset(&cmapgrid[0][0][0], 0, 6*CMAPDIM*CMAPDIM*sizeof(double));
ncrosstermtypes = 0;
memset(&cmapgrid[0][0][0], 0, CMAPMAX*CMAPDIM*CMAPDIM*sizeof(double));
utils::logmesg(lmp, "Reading CMAP parameters from: {}\n", cmapfile);
PotentialFileReader reader(lmp, cmapfile, "cmap grid");
// there are six maps in this order.
// there may be up to six maps.
// the charmm36.cmap file has in this order.
// alanine, alanine-proline, proline, proline-proline, glycine, glycine-proline.
// read as one big blob of numbers while ignoring comments
reader.next_dvector(&cmapgrid[0][0][0],6*CMAPDIM*CMAPDIM);
// custom CMAP files created by charmm-gui may have fewer entries
// read one map at a time as a blob of numbers while ignoring comments
// and stop reading when whe have reached EOF.
for (ncrosstermtypes = 0; ncrosstermtypes < CMAPMAX; ++ncrosstermtypes)
reader.next_dvector(&cmapgrid[ncrosstermtypes][0][0],CMAPDIM*CMAPDIM);
} catch (EOFException &) {
utils::logmesg(lmp, " Read in CMAP data for {} crossterm types\n", ncrosstermtypes);
} catch (std::exception &e) {
error->one(FLERR,"Error reading CMAP potential file: {}", e.what());
}
@ -934,10 +937,6 @@ void FixCMAP::read_data_header(char *line)
} catch (std::exception &e) {
error->all(FLERR,"Invalid read data header line for fix cmap: {}", e.what());
}
// not set in constructor because this fix could be defined before newton command
newton_bond = force->newton_bond;
}
/* ----------------------------------------------------------------------
@ -957,10 +956,10 @@ void FixCMAP::read_data_section(char * /*keyword*/, int /*n*/, char *buf,
// loop over lines of CMAP crossterms
// tokenize the line into values
// add crossterm to one of my atoms, depending on newton_bond
// add crossterm to one of my atoms
for (const auto &line : lines) {
ValueTokenizer values(line);
ValueTokenizer values(utils::trim_comment(line));
try {
values.skip();
itype = values.next_int();

3
src/MOLECULE/fix_cmap.h
View File

@ -65,8 +65,7 @@ class FixCMAP : public Fix {
double memory_usage() override;
private:
int nprocs, me;
int newton_bond, eflag_caller;
int eflag_caller;
int ctype, ilevel_respa;
int ncrosstermtypes, crossterm_per_atom, maxcrossterm;
int ncrosstermlist;

2
src/potential_file_reader.cpp
View File

@ -144,6 +144,8 @@ void PotentialFileReader::next_dvector(double *list, int n)
{
try {
return reader->next_dvector(list, n);
} catch (EOFException &) {
throw EOFException("EOF reached");
} catch (FileReaderException &e) {
error->one(FLERR, e.what());
}

5
src/text_file_reader.cpp
View File

@ -189,8 +189,9 @@ void TextFileReader::next_dvector(double *list, int n)
char *ptr = next_line();
if (ptr == nullptr) {
// EOF
if (i < n) {
if (i == 0) { // EOF without any records
throw EOFException("EOF reached");
} else if (i < n) { // EOF with incomplete data
throw FileReaderException(
fmt::format("Incorrect format in {} file! {}/{} values", filetype, i, n));
}