From 2bcac9efba857331a4821da4b7057241f6c1a0e0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 12 May 2025 12:13:05 -0400 Subject: [PATCH 01/18] fix spelling issue --- doc/utils/sphinx-config/false_positives.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index c81be91cbe..0c013abd56 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -397,6 +397,7 @@ Broglie brownian brownw Broyden +Bruenger Bruskin Brusselle Bryantsev From 75907ccf918333a106e084edd840e2ddceb8fb7d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 13 May 2025 02:32:27 -0400 Subject: [PATCH 02/18] add support to extract eflag/vflag_atom/global as global properties --- src/library.cpp | 48 +++++++++++++++++-- .../c-library/test_library_properties.cpp | 13 +++++ 2 files changed, 57 insertions(+), 4 deletions(-) diff --git a/src/library.cpp b/src/library.cpp index e6b14dafd3..a8acbade52 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -1622,6 +1622,11 @@ int lammps_extract_global_datatype(void * /*handle*/, const char *name) if (strcmp(name,"neigh_dihedrallist") == 0) return LAMMPS_INT_2D; if (strcmp(name,"neigh_improperlist") == 0) return LAMMPS_INT_2D; + if (strcmp(name,"eflag_global") == 0) return LAMMPS_BIGINT; + if (strcmp(name,"eflag_atom") == 0) return LAMMPS_BIGINT; + if (strcmp(name,"vflag_global") == 0) return LAMMPS_BIGINT; + if (strcmp(name,"vflag_atom") == 0) return LAMMPS_BIGINT; + if (strcmp(name,"map_style") == 0) return LAMMPS_INT; #if defined(LAMMPS_BIGBIG) if (strcmp(name,"map_tag_max") == 0) return LAMMPS_BIGINT; @@ -1707,6 +1712,7 @@ report the "native" data type. The following tables are provided: * :ref:`Simulation box settings ` * :ref:`System property settings ` * :ref:`Neighbor topology data ` +* :ref:`Energy and virial tally settings ` * :ref:`Git revision and version settings ` * :ref:`Unit settings ` @@ -1984,6 +1990,35 @@ Get length of lists with :ref:`lammps_extract_setting() atom->q_flag; - if (strcmp(name,"neigh_bondlist") == 0) return lmp->neighbor->bondlist; - if (strcmp(name,"neigh_anglelist") == 0) return lmp->neighbor->anglelist; - if (strcmp(name,"neigh_dihedrallist") == 0) return lmp->neighbor->dihedrallist; - if (strcmp(name,"neigh_improperlist") == 0) return lmp->neighbor->improperlist; + if (strcmp(name,"neigh_bondlist") == 0) return (void *) lmp->neighbor->bondlist; + if (strcmp(name,"neigh_anglelist") == 0) return (void *) lmp->neighbor->anglelist; + if (strcmp(name,"neigh_dihedrallist") == 0) return (void *) lmp->neighbor->dihedrallist; + if (strcmp(name,"neigh_improperlist") == 0) return (void *) lmp->neighbor->improperlist; + + if (strcmp(name,"eflag_global") == 0) return (void *) &lmp->update->eflag_global; + if (strcmp(name,"eflag_atom") == 0) return (void *) &lmp->update->eflag_atom; + if (strcmp(name,"vflag_global") == 0) return (void *) &lmp->update->vflag_global; + if (strcmp(name,"vflag_atom") == 0) return (void *) &lmp->update->vflag_atom; if (strcmp(name,"map_style") == 0) return (void *) &lmp->atom->map_style; if (strcmp(name,"map_tag_max") == 0) return (void *) &lmp->atom->map_tag_max; diff --git a/unittest/c-library/test_library_properties.cpp b/unittest/c-library/test_library_properties.cpp index e2149d98e6..033e9c041e 100644 --- a/unittest/c-library/test_library_properties.cpp +++ b/unittest/c-library/test_library_properties.cpp @@ -387,6 +387,19 @@ TEST_F(LibraryProperties, global) auto *b_ptr = (int64_t *)lammps_extract_global(lmp, "ntimestep"); EXPECT_EQ((*b_ptr), 2); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "eflag_global"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "eflag_global"); + EXPECT_EQ((*b_ptr), 2); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "eflag_atom"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "eflag_atom"); + EXPECT_EQ((*b_ptr), 0); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "vflag_global"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "vflag_global"); + EXPECT_EQ((*b_ptr), 2); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "vflag_atom"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "vflag_atom"); + EXPECT_EQ((*b_ptr), 0); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "dt"), LAMMPS_DOUBLE); auto *d_ptr = (double *)lammps_extract_global(lmp, "dt"); EXPECT_DOUBLE_EQ((*d_ptr), 0.1); From fdd91e597e8c508a5ce06b84ff75002d0ee63d25 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 13 May 2025 02:44:05 -0400 Subject: [PATCH 03/18] add a few more tests for extracted global properties --- .../c-library/test_library_properties.cpp | 29 +++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/unittest/c-library/test_library_properties.cpp b/unittest/c-library/test_library_properties.cpp index 033e9c041e..39b31e8217 100644 --- a/unittest/c-library/test_library_properties.cpp +++ b/unittest/c-library/test_library_properties.cpp @@ -383,10 +383,35 @@ TEST_F(LibraryProperties, global) char *c_ptr = (char *)lammps_extract_global(lmp, "units"); EXPECT_THAT(c_ptr, StrEq("real")); - EXPECT_EQ(lammps_extract_global_datatype(lmp, "ntimestep"), LAMMPS_INT64); - auto *b_ptr = (int64_t *)lammps_extract_global(lmp, "ntimestep"); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "ntimestep"), LAMMPS_BIGINT); + auto *b_ptr = (bigint *)lammps_extract_global(lmp, "ntimestep"); EXPECT_EQ((*b_ptr), 2); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "natoms"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "natoms"); + EXPECT_EQ((*b_ptr), 29); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "nbonds"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "nbonds"); + EXPECT_EQ((*b_ptr), 24); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "nangles"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "nangles"); + EXPECT_EQ((*b_ptr), 30); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "ndihedrals"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "ndihedrals"); + EXPECT_EQ((*b_ptr), 31); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "nimpropers"), LAMMPS_BIGINT); + b_ptr = (bigint *)lammps_extract_global(lmp, "nimpropers"); + EXPECT_EQ((*b_ptr), 2); + + EXPECT_EQ(lammps_extract_global_datatype(lmp, "neigh_bondlist"), LAMMPS_INT_2D); + EXPECT_NE(lammps_extract_global(lmp, "neigh_bondlist"), nullptr); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "neigh_anglelist"), LAMMPS_INT_2D); + EXPECT_NE(lammps_extract_global(lmp, "neigh_anglelist"), nullptr); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "neigh_dihedrallist"), LAMMPS_INT_2D); + EXPECT_NE(lammps_extract_global(lmp, "neigh_dihedrallist"), nullptr); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "neigh_improperlist"), LAMMPS_INT_2D); + EXPECT_NE(lammps_extract_global(lmp, "neigh_improperlist"), nullptr); + EXPECT_EQ(lammps_extract_global_datatype(lmp, "eflag_global"), LAMMPS_BIGINT); b_ptr = (bigint *)lammps_extract_global(lmp, "eflag_global"); EXPECT_EQ((*b_ptr), 2); From 83fa2cbc9313afb9b5d40258c6910c9c6264db01 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 13 May 2025 06:02:10 -0400 Subject: [PATCH 04/18] enable PotentialFileReader class to change line buffer size --- src/potential_file_reader.cpp | 10 ++++++++++ src/potential_file_reader.h | 1 + src/text_file_reader.cpp | 4 +++- 3 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/potential_file_reader.cpp b/src/potential_file_reader.cpp index 2a93a4a524..7bac388ba0 100644 --- a/src/potential_file_reader.cpp +++ b/src/potential_file_reader.cpp @@ -86,11 +86,21 @@ PotentialFileReader::~PotentialFileReader() /** Set comment (= text after '#') handling preference for the file to be read * * \param value Comment text is ignored if true, or not if false */ + void PotentialFileReader::ignore_comments(bool value) { reader->ignore_comments = value; } +/** Set line buffer size of the internal TextFileReader class instance. + * + * \param bufsize New size of the line buffer */ + +void PotentialFileReader::set_bufsize(int bufsize) +{ + reader->set_bufsize(bufsize); +} + /** Reset file to the beginning */ void PotentialFileReader::rewind() diff --git a/src/potential_file_reader.h b/src/potential_file_reader.h index c07b4b83f6..534457a2f8 100644 --- a/src/potential_file_reader.h +++ b/src/potential_file_reader.h @@ -41,6 +41,7 @@ class PotentialFileReader : protected Pointers { const int auto_convert = 0); ~PotentialFileReader() override; + void set_bufsize(int bufsize); void ignore_comments(bool value); void rewind(); diff --git a/src/text_file_reader.cpp b/src/text_file_reader.cpp index 715e82ab32..d598c274f0 100644 --- a/src/text_file_reader.cpp +++ b/src/text_file_reader.cpp @@ -93,7 +93,9 @@ TextFileReader::~TextFileReader() delete[] line; } -/** adjust line buffer size */ +/** adjust line buffer size + * + * \param newsize New size of the internal line buffer */ void TextFileReader::set_bufsize(int newsize) { From 9ba50df9d8e8c54fe045e0e641ce6a718ab0dbec Mon Sep 17 00:00:00 2001 From: Guillaume Fraux Date: Tue, 13 May 2025 16:09:13 +0200 Subject: [PATCH 05/18] Select one overload of log2 for the kokkos build --- src/KOKKOS/pair_kokkos.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h index 399142dfaf..9e04770ae3 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/src/KOKKOS/pair_kokkos.h @@ -961,7 +961,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, std::enable_if_t<(NEIGHFLAG&P lastcall = fpair->lmp->update->ntimestep; vectorsize = GetMaxNeighs(list); if (vectorsize == 0) vectorsize = 1; - vectorsize = MathSpecial::powint(2,(int(log2(vectorsize) + 0.5))); // round to nearest power of 2 + vectorsize = MathSpecial::powint(2,(int(log2(double(vectorsize)) + 0.5))); // round to nearest power of 2 #if defined(KOKKOS_ENABLE_HIP) int max_vectorsize = 64; From 01bde55e9a84d0508eed1697325fe21fd2c4114e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 13 May 2025 12:53:03 -0400 Subject: [PATCH 06/18] match all argument types for powint() --- src/KOKKOS/pair_kokkos.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h index 9e04770ae3..63e637c108 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/src/KOKKOS/pair_kokkos.h @@ -961,7 +961,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, std::enable_if_t<(NEIGHFLAG&P lastcall = fpair->lmp->update->ntimestep; vectorsize = GetMaxNeighs(list); if (vectorsize == 0) vectorsize = 1; - vectorsize = MathSpecial::powint(2,(int(log2(double(vectorsize)) + 0.5))); // round to nearest power of 2 + vectorsize = MathSpecial::powint(2.0,(int(log2(double(vectorsize)) + 0.5))); // round to nearest power of 2 #if defined(KOKKOS_ENABLE_HIP) int max_vectorsize = 64; From 32588f075e11e8fdfbffb4fa84bb8278512ce885 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 May 2025 02:53:12 -0400 Subject: [PATCH 07/18] remove assignments of variables to themselves --- src/KOKKOS/meam_dens_init_kokkos.h | 1 - src/KOKKOS/meam_force_kokkos.h | 1 - 2 files changed, 2 deletions(-) diff --git a/src/KOKKOS/meam_dens_init_kokkos.h b/src/KOKKOS/meam_dens_init_kokkos.h index dd63be96bd..23985c7082 100644 --- a/src/KOKKOS/meam_dens_init_kokkos.h +++ b/src/KOKKOS/meam_dens_init_kokkos.h @@ -236,7 +236,6 @@ MEAMKokkos::meam_dens_init(int inum_half, int ntype, typename AT::t_ this->d_neighbors_half = d_neighbors_half; this->d_neighbors_full = d_neighbors_full; this->d_offset = d_offset; - this->nlocal = nlocal; if (need_dup) { dup_rho0 = Kokkos::Experimental::create_scatter_view(d_rho0); diff --git a/src/KOKKOS/meam_force_kokkos.h b/src/KOKKOS/meam_force_kokkos.h index 1875e22dcf..703bee23d7 100644 --- a/src/KOKKOS/meam_force_kokkos.h +++ b/src/KOKKOS/meam_force_kokkos.h @@ -17,7 +17,6 @@ void MEAMKokkos::meam_force( { EV_FLOAT ev; - this->eflag_either = eflag_either; this->eflag_global = eflag_global; this->eflag_atom = eflag_atom; this->vflag_global = vflag_global; From 179d4f0148ce5d9ce3ff35f1bab513c3eaec96c5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 May 2025 02:53:51 -0400 Subject: [PATCH 08/18] work around C++ error --- src/LATBOLTZ/fix_lb_fluid.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LATBOLTZ/fix_lb_fluid.cpp b/src/LATBOLTZ/fix_lb_fluid.cpp index 286b56cab5..d5161bba42 100644 --- a/src/LATBOLTZ/fix_lb_fluid.cpp +++ b/src/LATBOLTZ/fix_lb_fluid.cpp @@ -2344,7 +2344,7 @@ void FixLbFluid::SetupBuffers() MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL, &dump_file_handle_raw); MPI_File_set_size(dump_file_handle_raw, 0); - MPI_File_set_view(dump_file_handle_raw, 0, MPI_DOUBLE, dump_file_mpitype, "native", + MPI_File_set_view(dump_file_handle_raw, 0, MPI_DOUBLE, dump_file_mpitype, (char *)"native", MPI_INFO_NULL); } } From 19cfd08eb87739c0532a74106a990b375cac358f Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 May 2025 08:15:01 -0400 Subject: [PATCH 09/18] correctly enable GPU package and make fix imd wait in background for bucky+cnt example --- examples/PACKAGES/imd/in.bucky-plus-cnt | 4 +-- examples/PACKAGES/imd/in.bucky-plus-cnt-gpu | 15 ++++---- examples/PACKAGES/imd/in.deca-ala_imd-gpu | 14 ++++---- examples/PACKAGES/imd/in.melt_imd-gpu | 38 +++++++++++---------- 4 files changed, 38 insertions(+), 33 deletions(-) diff --git a/examples/PACKAGES/imd/in.bucky-plus-cnt b/examples/PACKAGES/imd/in.bucky-plus-cnt index b3eeff3cc1..af511fe11f 100644 --- a/examples/PACKAGES/imd/in.bucky-plus-cnt +++ b/examples/PACKAGES/imd/in.bucky-plus-cnt @@ -46,8 +46,8 @@ fix integrate mobile nve fix thermostat mobile langevin 300.0 300.0 2000.0 234624 # IMD setup. -fix comm all imd 6789 unwrap on trate 10 -#fix comm all imd 6789 unwrap on trate 10 nowait on +#fix comm all imd 6789 unwrap on trate 10 +fix comm all imd 6789 unwrap on trate 10 nowait on # temperature is based on mobile atoms only compute mobtemp mobile temp diff --git a/examples/PACKAGES/imd/in.bucky-plus-cnt-gpu b/examples/PACKAGES/imd/in.bucky-plus-cnt-gpu index 5762ec68c8..f3e4b32cdc 100644 --- a/examples/PACKAGES/imd/in.bucky-plus-cnt-gpu +++ b/examples/PACKAGES/imd/in.bucky-plus-cnt-gpu @@ -1,16 +1,20 @@ # stick a buckyball into a nanotube + +# enable GPU package from within the input: +package gpu 0 pair/only on +suffix gpu + units real dimension 3 boundary f f f atom_style molecular -newton off processors * * 1 # read topology read_data data.bucky-plus-cnt -pair_style lj/cut/gpu 10.0 +pair_style lj/cut 10.0 bond_style harmonic angle_style charmm dihedral_style charmm @@ -29,9 +33,6 @@ neigh_modify delay 0 every 1 check yes timestep 2.0 -# required for GPU acceleration -fix gpu all gpu force 0 0 1.0 - # we only move some atoms. group mobile type 1 @@ -49,8 +50,8 @@ fix integrate mobile nve fix thermostat mobile langevin 300.0 300.0 2000.0 234624 # IMD setup. -fix comm all imd 6789 unwrap on trate 10 -#fix comm all imd 6789 unwrap on trate 10 nowait on +#fix comm all imd 6789 unwrap on trate 10 +fix comm all imd 6789 unwrap on trate 10 nowait on # temperature is based on mobile atoms only compute mobtemp mobile temp diff --git a/examples/PACKAGES/imd/in.deca-ala_imd-gpu b/examples/PACKAGES/imd/in.deca-ala_imd-gpu index 72c3f4aae9..9470f7c213 100644 --- a/examples/PACKAGES/imd/in.deca-ala_imd-gpu +++ b/examples/PACKAGES/imd/in.deca-ala_imd-gpu @@ -1,8 +1,12 @@ -# +# + +# enable GPU package from within the input: +package gpu 0 pair/only on +suffix gpu + units real neighbor 2.5 bin neigh_modify delay 1 every 1 -newton off atom_style full bond_style harmonic @@ -10,20 +14,18 @@ angle_style charmm dihedral_style charmm improper_style harmonic -pair_style lj/charmm/coul/long/gpu 8 10 +pair_style lj/charmm/coul/long 8 10 pair_modify mix arithmetic special_bonds charmm read_data data.deca-ala-solv -fix 0 all gpu force/neigh 0 0 1.0 - group peptide id <= 103 fix rigidh all shake 1e-6 100 1000 t 1 2 3 4 5 a 23 thermo 100 thermo_style multi timestep 2.0 -kspace_style pppm/gpu 1e-5 +kspace_style pppm 1e-5 fix ensemble all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 0.2 diff --git a/examples/PACKAGES/imd/in.melt_imd-gpu b/examples/PACKAGES/imd/in.melt_imd-gpu index 24904eb832..f1406befa6 100644 --- a/examples/PACKAGES/imd/in.melt_imd-gpu +++ b/examples/PACKAGES/imd/in.melt_imd-gpu @@ -1,30 +1,32 @@ -# 3d Lennard-Jones melt +# 3d Lennard-Jones melt with GPU package acceleration -units lj -atom_style atomic -newton off +# enable GPU package from within the input: +package gpu 0 +suffix gpu -lattice fcc 0.8442 -region box block 0 10 0 10 0 10 -create_box 1 box -create_atoms 1 box -mass 1 1.0 +units lj +atom_style atomic -velocity all create 3.0 87287 +lattice fcc 0.8442 +region box block 0 10 0 10 0 10 +create_box 1 box +create_atoms 1 box +mass 1 1.0 -pair_style lj/cut/gpu 2.5 -pair_coeff 1 1 1.0 1.0 2.5 +velocity all create 3.0 87287 -neighbor 0.3 bin -neigh_modify every 5 delay 10 check yes +pair_style lj/cut 2.5 +pair_coeff 1 1 1.0 1.0 2.5 + +neighbor 0.3 bin +neigh_modify every 5 delay 10 check yes thermo_style custom step pe ke spcpu -fix 0 all gpu force/neigh 0 0 1.0 -fix 1 all nve +fix 1 all nve # IMD setup. fix comm all imd 5678 unwrap off fscale 20.0 trate 20 nowait on -thermo 500 -run 5000000 +thermo 500 +run 5000000 From c3c01806490feaad10b719eb21b2a9ef6496d420 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 May 2025 08:15:35 -0400 Subject: [PATCH 10/18] correctly check for 32bit integer overflow --- src/ML-POD/compute_podd_atom.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ML-POD/compute_podd_atom.cpp b/src/ML-POD/compute_podd_atom.cpp index 4ab6e23393..c691b6b3f2 100644 --- a/src/ML-POD/compute_podd_atom.cpp +++ b/src/ML-POD/compute_podd_atom.cpp @@ -58,8 +58,8 @@ ComputePODDAtom::ComputePODDAtom(LAMMPS *lmp, int narg, char **arg) : pod = nullptr; elements = nullptr; - if (((((MAXBIGINT*3.0)*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0)) - error->all(FLERR, "Per-atom data too large"); + if ((((3.0*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0)) + error->all(FLERR, "Too many atoms ({}_ for compute {}", atom->natoms, style); size_peratom_cols = 3 * atom->natoms * podptr->Mdesc * podptr->nClusters; peratom_flag = 1; } @@ -110,8 +110,8 @@ void ComputePODDAtom::compute_peratom() if (atom->natoms > nmax) { memory->destroy(pod); nmax = atom->natoms; - if (((((MAXBIGINT*3.0)*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0)) - error->all(FLERR, "Per-atom data too large"); + if ((((3.0*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0)) + error->all(FLERR, "Too many atoms ({}) for compute {}", atom->natoms, style); int numdesc = 3 * atom->natoms * podptr->Mdesc * podptr->nClusters; memory->create(pod, nmax, numdesc,"podd/atom:pod"); array_atom = pod; From 5fd4d6bb8722d1c9a385acb529b7950af2d2d662 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 May 2025 23:15:26 -0400 Subject: [PATCH 11/18] remove directionCorrection and update docs for compute pair/local and compute bond/local --- doc/src/compute_bond_local.rst | 32 +++++++++++++++++++++----------- doc/src/compute_pair_local.rst | 27 +++++++++++++++++++-------- src/compute_bond_local.cpp | 9 +++------ src/compute_pair_local.cpp | 9 +++------ 4 files changed, 46 insertions(+), 31 deletions(-) diff --git a/doc/src/compute_bond_local.rst b/doc/src/compute_bond_local.rst index e070d507b1..9aa062f911 100644 --- a/doc/src/compute_bond_local.rst +++ b/doc/src/compute_bond_local.rst @@ -64,22 +64,32 @@ All these properties are computed for the pair of atoms in a bond, whether the two atoms represent a simple diatomic molecule, or are part of some larger molecule. -The value *dist* is the current length of the bond. -The values *dx*, *dy*, and *dz* are the xyz components of the -*distance* between the pair of atoms. This value is always the -distance from the atom of lower to the one with the higher id. +.. versionchanged:: TBD + + The sign of *dx*, *dy*, *dz* is no longer determined by the atom IDs + of the bonded atoms but by their order in the bond list to be + consistent with *fx*, *fy*, and *fz*. + +The value *dist* is the current length of the bond. The values *dx*, +*dy*, and *dz* are the :math:`(x,y,z)` components of the distance vector +:math:`\vec{x_i} - \vec{x_j}` between the atoms in the bond. The order +of the atoms is determined by the bond list and the respective atom-IDs +can be output with :doc:`compute property/local +`. The value *engpot* is the potential energy for the bond, based on the current separation of the pair of atoms in the bond. -The value *force* is the magnitude of the force acting between the -pair of atoms in the bond. +The value *force* is the magnitude of the force acting between the pair +of atoms in the bond, which is positive for a repulsive force and +negative for an attractive force. -The values *fx*, *fy*, and *fz* are the xyz components of -*force* between the pair of atoms in the bond. For bond styles that apply -non-central forces, such as :doc:`bond_style bpm/rotational -`, these values only include the :math:`(x,y,z)` -components of the normal force component. +The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of +the force vector on atom *i* obtained by projecting *force* on the +distance vector. For bond styles that apply non-central forces, such as +:doc:`bond_style bpm/rotational `, these values +only include the :math:`(x,y,z)` components of the normal force +component. The remaining properties are all computed for motion of the two atoms relative to the center of mass (COM) velocity of the two atoms in the diff --git a/doc/src/compute_pair_local.rst b/doc/src/compute_pair_local.rst index 31209f63f4..5e8ac73a22 100644 --- a/doc/src/compute_pair_local.rst +++ b/doc/src/compute_pair_local.rst @@ -56,19 +56,30 @@ force cutoff distance for that interaction, as defined by the :doc:`pair_style ` and :doc:`pair_coeff ` commands. -The value *dist* is the distance between the pair of atoms. -The values *dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the -*distance* between the pair of atoms. This value is always the -distance from the atom of higher to the one with the lower atom ID. +.. versionchanged:: TBD + + The sign of *dx*, *dy*, *dz* is no longer determined by the value of + their atom-IDs but by their order in the neighbor list to be + consistent with *fx*, *fy*, and *fz*. + +The value *dist* is the distance between the pair of atoms. The values +*dx*, *dy*, and *dz* are the :math:`(x,y,z)` components of the distance +vector :math:`\vec{x_i} - \vec{x_j}` between the pair of atoms. The +order of the atoms is determined by the neighbor list and the respective +atom-IDs can be output with :doc:`compute property/local +`. The value *eng* is the interaction energy for the pair of atoms. The value *force* is the force acting between the pair of atoms, which is positive for a repulsive force and negative for an attractive -force. The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of -*force* on atom I. For pair styles that apply non-central forces, -such as :doc:`granular pair styles `, these values only include -the :math:`(x,y,z)` components of the normal force component. +force. + +The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of +the force vector on atom *i* obtained by projecting *force* on the +distance vector. For pair styles that apply non-central forces, such as +:doc:`granular pair styles `, these values only include the +:math:`(x,y,z)` components of the normal force component. A pair style may define additional pairwise quantities which can be accessed as *p1* to *pN*, where :math:`N` is defined by the pair style. diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index e9632d254f..6354c67638 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -428,22 +428,19 @@ int ComputeBondLocal::compute_bonds(int flag) if (dstr) input->variable->internal_set(dvar, sqrt(rsq)); } - // to make sure dx, dy and dz are always from the lower to the higher id - double directionCorrection = tag[atom1] > tag[atom2] ? -1.0 : 1.0; - for (int n = 0; n < nvalues; n++) { switch (bstyle[n]) { case DIST: ptr[n] = sqrt(rsq); break; case DX: - ptr[n] = dx * directionCorrection; + ptr[n] = dx; break; case DY: - ptr[n] = dy * directionCorrection; + ptr[n] = dy; break; case DZ: - ptr[n] = dz * directionCorrection; + ptr[n] = dz; break; case ENGPOT: ptr[n] = engpot; diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index 57f15264f0..fa5d164844 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -277,22 +277,19 @@ int ComputePairLocal::compute_pairs(int flag) else ptr = alocal[m]; - // to make sure dx, dy and dz are always from the lower to the higher id - double directionCorrection = itag > jtag ? -1.0 : 1.0; - for (n = 0; n < nvalues; n++) { switch (pstyle[n]) { case DIST: ptr[n] = sqrt(rsq); break; case DX: - ptr[n] = delx * directionCorrection; + ptr[n] = delx; break; case DY: - ptr[n] = dely * directionCorrection; + ptr[n] = dely; break; case DZ: - ptr[n] = delz * directionCorrection; + ptr[n] = delz; break; case ENG: ptr[n] = eng; From 629ec2eabe21dd68aa4079a7081ec5a950ddc890 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 May 2025 23:16:31 -0400 Subject: [PATCH 12/18] update for consistency with docs --- src/read_data.cpp | 4 ++-- src/reset_atoms_mol.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/read_data.cpp b/src/read_data.cpp index 79d88148c5..cf7b224db2 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -174,13 +174,13 @@ void ReadData::command(int narg, char **arg) addflag = VALUE; bigint offset = utils::bnumeric(FLERR, arg[iarg + 1], false, lmp); if (offset > MAXTAGINT) - error->all(FLERR, "Read data add atomID offset {} is too big", offset); + error->all(FLERR, "Read data add IDoffset {} is too big", offset); id_offset = offset; if (atom->molecule_flag) { offset = utils::bnumeric(FLERR, arg[iarg + 2], false, lmp); if (offset > MAXTAGINT) - error->all(FLERR, "Read data add molID offset {} is too big", offset); + error->all(FLERR, "Read data add MOLoffset {} is too big", offset); mol_offset = offset; iarg++; } diff --git a/src/reset_atoms_mol.cpp b/src/reset_atoms_mol.cpp index 54d3bbcc76..363e2f08eb 100644 --- a/src/reset_atoms_mol.cpp +++ b/src/reset_atoms_mol.cpp @@ -211,7 +211,7 @@ void ResetAtomsMol::reset() // if offset < 0 (default), reset it // if group = all, offset = 0 - // else offset = largest molID of non-group atoms + // else offset = largest molecule ID of non-group atoms if (offset < 0) { if (groupbit != 1) { From 031fab210feb6d67bc8203af4c176a9ac91c1a5b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 May 2025 23:28:46 -0400 Subject: [PATCH 13/18] update unittest for change in compute pair/local --- unittest/python/python-numpy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unittest/python/python-numpy.py b/unittest/python/python-numpy.py index 4930527a61..be9109b9a3 100644 --- a/unittest/python/python-numpy.py +++ b/unittest/python/python-numpy.py @@ -153,7 +153,7 @@ class PythonNumpy(unittest.TestCase): self.assertEqual(values[0,0], 0.5) self.assertEqual(values[0,3], -0.5) self.assertEqual(values[1,0], 1.5) - self.assertEqual(values[1,3], 1.5) + self.assertEqual(values[1,3], -1.5) def testExtractAtom(self): self.lmp.command("units lj") From 63ee449dc1ec17f2c1403b9b24492e27b97dd8bd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 17 May 2025 21:15:48 -0400 Subject: [PATCH 14/18] reformulate description of force components --- doc/src/compute_bond_local.rst | 12 +++++++----- doc/src/compute_pair_local.rst | 7 +++++-- 2 files changed, 12 insertions(+), 7 deletions(-) diff --git a/doc/src/compute_bond_local.rst b/doc/src/compute_bond_local.rst index 9aa062f911..7ce6f9b15a 100644 --- a/doc/src/compute_bond_local.rst +++ b/doc/src/compute_bond_local.rst @@ -85,11 +85,13 @@ of atoms in the bond, which is positive for a repulsive force and negative for an attractive force. The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of -the force vector on atom *i* obtained by projecting *force* on the -distance vector. For bond styles that apply non-central forces, such as -:doc:`bond_style bpm/rotational `, these values -only include the :math:`(x,y,z)` components of the normal force -component. +the force on the first atom *i* in the bond due to the second atom *j*. +Mathematically, they are obtained by multiplying the value of *force* +from above with a unit vector created from the *dx*, *dy*, and *dz* +components of the distance vector also described above. For bond styles +that apply non-central forces, such as :doc:`bond_style bpm/rotational +`, these values only include the :math:`(x,y,z)` +components of the normal force component. The remaining properties are all computed for motion of the two atoms relative to the center of mass (COM) velocity of the two atoms in the diff --git a/doc/src/compute_pair_local.rst b/doc/src/compute_pair_local.rst index 5e8ac73a22..605bfc8e9e 100644 --- a/doc/src/compute_pair_local.rst +++ b/doc/src/compute_pair_local.rst @@ -76,8 +76,11 @@ is positive for a repulsive force and negative for an attractive force. The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of -the force vector on atom *i* obtained by projecting *force* on the -distance vector. For pair styles that apply non-central forces, such as +the force vector on the first atom *i* of a pair in the neighbor list +due to the second atom *j*. Mathematically, they are obtained by +multiplying the value of *force* from above with a unit vector created +from the *dx*, *dy*, and *dz* components of the distance vector also +described above. For pair styles that apply non-central forces, such as :doc:`granular pair styles `, these values only include the :math:`(x,y,z)` components of the normal force component. From e393b9803bebb86910af976769db1379cc801e88 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 17 May 2025 22:41:42 -0400 Subject: [PATCH 15/18] add information about the OPC model to the TIP4P howto --- doc/src/Howto_tip4p.rst | 129 ++++++++++++++++++-- doc/utils/sphinx-config/false_positives.txt | 3 + 2 files changed, 123 insertions(+), 9 deletions(-) diff --git a/doc/src/Howto_tip4p.rst b/doc/src/Howto_tip4p.rst index 47a1b9b578..5f38ea3c62 100644 --- a/doc/src/Howto_tip4p.rst +++ b/doc/src/Howto_tip4p.rst @@ -1,5 +1,5 @@ -TIP4P water model -================= +TIP4P and OPC water models +========================== The four-point TIP4P rigid water model extends the traditional :doc:`three-point TIP3P ` model by adding an additional @@ -9,9 +9,11 @@ the oxygen along the bisector of the HOH bond angle. A bond style of :doc:`harmonic ` and an angle style of :doc:`harmonic ` or :doc:`charmm ` should also be used. In case of rigid bonds also bond style :doc:`zero ` and angle -style :doc:`zero ` can be used. +style :doc:`zero ` can be used. Very similar to the TIP4P +model is the OPC water model. It can be realized the same way as TIP4P +but has different geometry and force field parameters. -There are two ways to implement TIP4P water in LAMMPS: +There are two ways to implement TIP4P-like water in LAMMPS: #. Use a specially written pair style that uses the :ref:`TIP3P geometry ` without the point M. The point M location is then @@ -21,7 +23,10 @@ There are two ways to implement TIP4P water in LAMMPS: computationally very efficient, but the charge distribution in space is only correct within the tip4p labeled styles. So all other computations using charges will "see" the negative charge incorrectly - on the oxygen atom. + located on the oxygen atom unless they are specially written for using + the TIP4P geometry internally as well, e.g. :doc:`compute dipole/tip4p + `, :doc:`fix efield/tip4p `, or + :doc:`kspace_style pppm/tip4p `. This can be done with the following pair styles for Coulomb with a cutoff: @@ -68,77 +73,90 @@ TIP4P/2005 model :ref:`(Abascal2) ` and a version of TIP4P parameters adjusted for use with a long-range Coulombic solver (e.g. Ewald or PPPM in LAMMPS). Note that for implicit TIP4P models the OM distance is specified in the :doc:`pair_style ` command, -not as part of the pair coefficients. +not as part of the pair coefficients. Also parameters for the OPC +model (:ref:`Izadi `) are provided. .. list-table:: :header-rows: 1 - :widths: 36 19 13 15 17 + :widths: 40 12 12 14 11 11 * - Parameter - TIP4P (original) - TIP4P/Ice - TIP4P/2005 - TIP4P (Ewald) + - OPC * - O mass (amu) - 15.9994 - 15.9994 - 15.9994 - 15.9994 + - 15.9994 * - H mass (amu) - 1.008 - 1.008 - 1.008 - 1.008 + - 1.008 * - O or M charge (:math:`e`) - -1.040 - -1.1794 - -1.1128 - -1.04844 + - -1.3582 * - H charge (:math:`e`) - 0.520 - 0.5897 - 0.5564 - 0.52422 + - 0.6791 * - LJ :math:`\epsilon` of OO (kcal/mole) - 0.1550 - 0.21084 - 0.1852 - 0.16275 + - 0.21280 * - LJ :math:`\sigma` of OO (:math:`\AA`) - 3.1536 - 3.1668 - 3.1589 - 3.16435 + - 3.1660 * - LJ :math:`\epsilon` of HH, MM, OH, OM, HM (kcal/mole) - 0.0 - 0.0 - 0.0 - 0.0 + - 0.0 * - LJ :math:`\sigma` of HH, MM, OH, OM, HM (:math:`\AA`) - 1.0 - 1.0 - 1.0 - 1.0 + - 1.0 * - :math:`r_0` of OH bond (:math:`\AA`) - 0.9572 - 0.9572 - 0.9572 - 0.9572 + - 0.8724 * - :math:`\theta_0` of HOH angle - 104.52\ :math:`^{\circ}` - 104.52\ :math:`^{\circ}` - 104.52\ :math:`^{\circ}` - 104.52\ :math:`^{\circ}` + - 103.60\ :math:`^{\circ}` * - OM distance (:math:`\AA`) - 0.15 - 0.1577 - 0.1546 - 0.1250 + - 0.1594 -Note that the when using the TIP4P pair style, the neighbor list cutoff +Note that the when using a TIP4P pair style, the neighbor list cutoff for Coulomb interactions is effectively extended by a distance 2 \* (OM distance), to account for the offset distance of the fictitious charges -on O atoms in water molecules. Thus it is typically best in an +on O atoms in water molecules. Thus, it is typically best in an efficiency sense to use a LJ cutoff >= Coulomb cutoff + 2\*(OM distance), to shrink the size of the neighbor list. This leads to slightly larger cost for the long-range calculation, so you can test the @@ -192,6 +210,94 @@ file changed): run 20000 write_data tip4p-implicit.data nocoeff +When constructing an OPC model, we cannot use the ``tip3p.mol`` file due +to the different geometry. Below is a molecule file providing the 3 +sites for an implicit OPC geometry with TIP4P styles. Note, that the +"Shake" and "Special" sections are missing here. Those will be +auto-generated since the molecule file is loaded *after* the simulation +box has been created. They are required only when the molecule file +is loaded *before*. + +.. _opc3p_molecule: +.. code-block:: + + # Water molecule. 3 point geometry for OPC model + + 3 atoms + 2 bonds + 1 angles + + Coords + + 1 0.00000 -0.06037 0.00000 + 2 0.68558 0.50250 0.00000 + 3 -0.68558 0.50250 0.00000 + + Types + + 1 1 # O + 2 2 # H + 3 2 # H + + Charges + + 1 -1.3582 + 2 0.6791 + 3 0.6791 + + Bonds + + 1 1 1 2 + 2 1 1 3 + + Angles + + 1 1 2 1 3 + +Below is a LAMMPS input file using the implicit method to implement +the OPC model using the molecule file from above and including the +PPPM long-range Coulomb solver. + +.. code-block:: LAMMPS + + units real + atom_style full + region box block -5 5 -5 5 -5 5 + create_box 2 box bond/types 1 angle/types 1 & + extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2 + + mass 1 15.9994 + mass 2 1.008 + + pair_style lj/cut/tip4p/long 1 2 1 1 0.1594 12.0 + pair_coeff 1 1 0.2128 3.166 + pair_coeff 2 2 0.0 1.0 + + bond_style zero + bond_coeff 1 0.8724 + + angle_style zero + angle_coeff 1 103.6 + + kspace_style pppm/tip4p 1.0e-5 + + molecule water opc3p.mol # this file has the OPC geometry but is without M + create_atoms 0 random 33 34564 NULL mol water 25367 overlap 1.33 + + fix rigid all shake 0.001 10 10000 b 1 a 1 + minimize 0.0 0.0 1000 10000 + + reset_timestep 0 + timestep 1.0 + velocity all create 300.0 5463576 + fix integrate all nvt temp 300 300 100.0 + + thermo_style custom step temp press etotal pe + + thermo 1000 + run 20000 + write_data opc-implicit.data nocoeff + Below is the code for a LAMMPS input file using the explicit method and a TIP4P molecule file. Because of using :doc:`fix rigid/small ` no bonds need to be defined and thus no extra storage needs @@ -279,3 +385,8 @@ Phys, 79, 926 (1983). **(Abascal2)** Abascal, J Chem Phys, 123, 234505 (2005) https://doi.org/10.1063/1.2121687 + +.. _Izadi: + +**(Izadi)** Izadi, Anandakrishnan, Onufriev, J. Phys. Chem. Lett., 5, 21, 3863 (2014) + https://doi.org/10.1021/jz501780a diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 0c013abd56..c87637a26e 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -103,6 +103,7 @@ Amit amsmath amu Amzallag +Anandakrishnan analytical Anders Andric @@ -1732,6 +1733,7 @@ Iyz iz izcm ized +Izadi Izrailev Izumi Izvekov @@ -2808,6 +2810,7 @@ oneMKL oneway onlysalt ons +Onufriev OO Oord opencl From 2968a62937ff06cd73e8061cc76be588163b6eb7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 17 May 2025 23:30:18 -0400 Subject: [PATCH 16/18] continue refactoring for bio force field and water moldel discussions --- doc/src/Howto.rst | 1 + doc/src/Howto_FFgeneral.rst | 55 +++++++++++++++++++++++++++++++++++++ doc/src/Howto_bioFF.rst | 38 ++++++++++--------------- doc/src/Howto_spc.rst | 4 +-- doc/src/Howto_tip4p.rst | 10 +++---- 5 files changed, 78 insertions(+), 30 deletions(-) create mode 100644 doc/src/Howto_FFgeneral.rst diff --git a/doc/src/Howto.rst b/doc/src/Howto.rst index ec90472c27..cdc4efd737 100644 --- a/doc/src/Howto.rst +++ b/doc/src/Howto.rst @@ -66,6 +66,7 @@ Force fields howto :name: force_howto :maxdepth: 1 + Howto_FFgeneral Howto_bioFF Howto_amoeba Howto_tip3p diff --git a/doc/src/Howto_FFgeneral.rst b/doc/src/Howto_FFgeneral.rst new file mode 100644 index 0000000000..1b96ae1119 --- /dev/null +++ b/doc/src/Howto_FFgeneral.rst @@ -0,0 +1,55 @@ +Some general force field considerations +======================================= + +A compact summary of the concepts, definitions, and properties of force +fields with explicit bonded interactions (like the ones discussed in +this HowTo) is given in :ref:`(Gissinger) `. + +A force field has 2 parts: the formulas that define its potential +functions and the coefficients used for a particular system. To assign +parameters it is first required to assign atom types. Those are not +only based on the elements, but also on the chemical environment due to +the atoms bound to them. This often follows the chemical concept of +*functional groups*. Example: a carbon atom bound with a single bond to +a single OH-group (alcohol) would be a different atom type than a carbon +atom bound to a methyl CH3 group (aliphatic carbon). The atom types +usually then determine the non-bonded Lennard-Jones parameters and the +parameters for bonds, angles, dihedrals, and impropers. On top of that, +partial charges have to be applied. Those are usually independent of +the atom types and are determined either for groups of atoms called +residues with some fitting procedure based on quantum mechanical +calculations, or based on some increment system that add or subtract +increments from the partial charge of an atom based on the types of +the neighboring atoms. + +Force fields differ in the strategies they employ to determine the +parameters and charge distribution in how generic or specific they are +which in turn has an impact on the accuracy (compare for example +CGenFF to CHARMM and GAFF to Amber). Because of the different +strategies, it is not a good idea to use a mix of parameters from +different force field *families* (like CHARMM, Amber, or GROMOS) +and that extends to the parameters for the solvent, especially +water. The publication describing the parameterization of a force +field will describe which water model to use. Changing the water +model usually leads to overall worse results (even if it may improve +on the water itself). + +In addition, one has to consider that *families* of force fields like +CHARMM, Amber, OPLS, or GROMOS have evolved over time and thus provide +different *revisions* of the force field parameters. These often +corresponds to changes in the functional form or the parameterization +strategies. This may also result in changes required for simulation +settings like the preferred cutoff or how Coulomb interactions are +computed (cutoff, smoothed/shifted cutoff, or long-range with Ewald +summation or equivalent). Unless explicitly stated in the publication +describing the force field, the Coulomb interaction cannot be chosen at +will but must match the revision of the force field. That said, +liberties may be taken during the initial equilibration of a system to +speed up the process, but not for production simulations. + +---------- + +.. _Typelabel2: + +**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024). + diff --git a/doc/src/Howto_bioFF.rst b/doc/src/Howto_bioFF.rst index 92dd45b9b6..cf8e4ab13e 100644 --- a/doc/src/Howto_bioFF.rst +++ b/doc/src/Howto_bioFF.rst @@ -1,22 +1,16 @@ CHARMM, AMBER, COMPASS, DREIDING, and OPLS force fields ======================================================= -A compact summary of the concepts, definitions, and properties of -force fields with explicit bonded interactions (like the ones discussed -in this HowTo) is given in :ref:`(Gissinger) `. - -A force field has 2 parts: the formulas that define it and the -coefficients used for a particular system. Here we only discuss -formulas implemented in LAMMPS that correspond to formulas commonly used -in the CHARMM, AMBER, COMPASS, and DREIDING force fields. Setting -coefficients is done either from special sections in an input data file -via the :doc:`read_data ` command or in the input script with -commands like :doc:`pair_coeff ` or :doc:`bond_coeff -` and so on. See the :doc:`Tools ` doc page for -additional tools that can use CHARMM, AMBER, or Materials Studio -generated files to assign force field coefficients and convert their -output into LAMMPS input. LAMMPS input scripts can also be generated by -`charmm-gui.org `_. +Here we only discuss formulas implemented in LAMMPS that correspond to +formulas commonly used in the CHARMM, AMBER, COMPASS, and DREIDING force +fields. Setting coefficients is done either from special sections in an +input data file via the :doc:`read_data ` command or in the +input script with commands like :doc:`pair_coeff ` or +:doc:`bond_coeff ` and so on. See the :doc:`Tools ` +doc page for additional tools that can use CHARMM, AMBER, or Materials +Studio generated files to assign force field coefficients and convert +their output into LAMMPS input. LAMMPS input scripts can also be +generated by `charmm-gui.org `_. CHARMM and AMBER ---------------- @@ -203,9 +197,11 @@ rather than individual force constants and geometric parameters that depend on the particular combinations of atoms involved in the bond, angle, or torsion terms. DREIDING has an :doc:`explicit hydrogen bond term ` to describe interactions involving a -hydrogen atom on very electronegative atoms (N, O, F). Unlike CHARMM -or AMBER, the DREIDING force field has not been parameterized for -considering solvents (like water). +hydrogen atom on very electronegative atoms (N, O, F). Unlike CHARMM or +AMBER, the DREIDING force field has not been parameterized for +considering solvents (like water) and has no rules for assigning +(partial) charges. That will seriously limit its accuracy when used for +simulating systems where those matter. See :ref:`(Mayo) ` for a description of the DREIDING force field @@ -272,10 +268,6 @@ compatible with a subset of OPLS interactions. ---------- -.. _Typelabel2: - -**(Gissinger)** J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz, A. Kohlmeyer, and E. B. Tadmor, J Phys Chem B, 128, 3282-3297 (2024). - .. _howto-MacKerell: **(MacKerell)** MacKerell, Bashford, Bellott, Dunbrack, Evanseck, Field, Fischer, Gao, Guo, Ha, et al (1998). J Phys Chem, 102, 3586 . https://doi.org/10.1021/jp973084f diff --git a/doc/src/Howto_spc.rst b/doc/src/Howto_spc.rst index 00bd8a1b10..f84d7797d2 100644 --- a/doc/src/Howto_spc.rst +++ b/doc/src/Howto_spc.rst @@ -1,5 +1,5 @@ -SPC water model -=============== +SPC and SPC/E water model +========================= The SPC water model specifies a 3-site rigid water molecule with charges and Lennard-Jones parameters assigned to each of the three atoms. diff --git a/doc/src/Howto_tip4p.rst b/doc/src/Howto_tip4p.rst index 5f38ea3c62..76c470d615 100644 --- a/doc/src/Howto_tip4p.rst +++ b/doc/src/Howto_tip4p.rst @@ -212,11 +212,11 @@ file changed): When constructing an OPC model, we cannot use the ``tip3p.mol`` file due to the different geometry. Below is a molecule file providing the 3 -sites for an implicit OPC geometry with TIP4P styles. Note, that the -"Shake" and "Special" sections are missing here. Those will be -auto-generated since the molecule file is loaded *after* the simulation -box has been created. They are required only when the molecule file -is loaded *before*. +sites of an implicit OPC geometry for use with TIP4P styles. Note, that +the "Shake" and "Special" sections are missing here. Those will be +auto-generated by LAMMPS when the molecule file is loaded *after* the +simulation box has been created. These sections are required only when +the molecule file is loaded *before*. .. _opc3p_molecule: .. code-block:: From a821654ef53a4650f61c56f512e5a046d570492a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 19 May 2025 11:30:51 -0400 Subject: [PATCH 17/18] correct error message --- src/ML-POD/compute_podd_atom.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ML-POD/compute_podd_atom.cpp b/src/ML-POD/compute_podd_atom.cpp index c691b6b3f2..364bfb54ac 100644 --- a/src/ML-POD/compute_podd_atom.cpp +++ b/src/ML-POD/compute_podd_atom.cpp @@ -59,7 +59,7 @@ ComputePODDAtom::ComputePODDAtom(LAMMPS *lmp, int narg, char **arg) : elements = nullptr; if ((((3.0*atom->natoms)*podptr->nClusters)*podptr->Mdesc) > (MAXSMALLINT*1.0)) - error->all(FLERR, "Too many atoms ({}_ for compute {}", atom->natoms, style); + error->all(FLERR, "Too many atoms ({}) for compute {}", atom->natoms, style); size_peratom_cols = 3 * atom->natoms * podptr->Mdesc * podptr->nClusters; peratom_flag = 1; } From b9b59bd23c93bedcfd003bafad91673d96960eb9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 19 May 2025 11:31:42 -0400 Subject: [PATCH 18/18] small clarifications and corrections. Sync with current state of affairs --- doc/src/Intro_authors.rst | 5 +++-- doc/src/Intro_portability.rst | 25 +++++++++++++++++-------- 2 files changed, 20 insertions(+), 10 deletions(-) diff --git a/doc/src/Intro_authors.rst b/doc/src/Intro_authors.rst index 38f1102595..730cd2e336 100644 --- a/doc/src/Intro_authors.rst +++ b/doc/src/Intro_authors.rst @@ -84,8 +84,9 @@ lammps.org". General questions about LAMMPS should be posted in the \normalsize -Past developers include Paul Crozier and Mark Stevens, both at SNL, -and Ray Shan, now at Materials Design. +Past core developers include Paul Crozier and Mark Stevens, both at SNL, +and Ray Shan while at SNL and later at Materials Design, now at Thermo +Fisher Scientific. ---------- diff --git a/doc/src/Intro_portability.rst b/doc/src/Intro_portability.rst index 63ae147b8c..564cdc47f4 100644 --- a/doc/src/Intro_portability.rst +++ b/doc/src/Intro_portability.rst @@ -28,8 +28,9 @@ Build systems LAMMPS can be compiled from source code using a (traditional) build system based on shell scripts, a few shell utilities (grep, sed, cat, tr) and the GNU make program. This requires running within a Bourne -shell (``/bin/sh``). Alternatively, a build system with different back -ends can be created using CMake. CMake must be at least version 3.16. +shell (``/bin/sh`` or ``/bin/bash``). Alternatively, a build system +with different back ends can be created using CMake. CMake must be +at least version 3.16. Operating systems ^^^^^^^^^^^^^^^^^ @@ -40,11 +41,18 @@ Also, compilation and correct execution on macOS and Windows (using Microsoft Visual C++) is checked automatically for the largest part of the source code. Some (optional) features are not compatible with all operating systems, either through limitations of the corresponding -LAMMPS source code or through incompatibilities of source code or -build system of required external libraries or packages. +LAMMPS source code or through incompatibilities or build system +limitations of required external libraries or packages. -Executables for Windows may be created natively using either Cygwin or -Visual Studio or with a Linux to Windows MinGW cross-compiler. +Executables for Windows may be created either natively using Cygwin, +MinGW, Intel, Clang, or Microsoft Visual C++ compilers, or with a Linux +to Windows MinGW cross-compiler. Native compilation is supported using +Microsoft Visual Studio or a terminal window (using the CMake build +system). + +Executables for macOS may be created either using Xcode or GNU compilers +installed with Homebrew. In the latter case, building of LAMMPS through +Homebrew instead of a manual compile is also possible. Additionally, FreeBSD and Solaris have been tested successfully to run LAMMPS and produce results consistent with those on Linux. @@ -61,8 +69,9 @@ CPU architectures ^^^^^^^^^^^^^^^^^ The primary CPU architecture for running LAMMPS is 64-bit x86, but also -32-bit x86, and 64-bit ARM and PowerPC (64-bit, Little Endian) are -regularly tested. +64-bit ARM and PowerPC (64-bit, Little Endian) are currently regularly +tested. Further architectures are tested by Linux distributions that +bundle LAMMPS. Portability compliance ^^^^^^^^^^^^^^^^^^^^^^