From 0e57af7222daff075001bbd7ab8245436b3bb983 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 01:49:59 +0000 Subject: [PATCH 01/43] KOKKOS: claim the host write the dpd/fdt/energy reverse comm makes pair_style dpd/fdt/energy/kk finishes its compute by copying duCond and duMech to the host and running the reverse communication, which adds the ghost contributions into the plain host arrays. That write was never declared, so the two copies were left apart with counters that said they agreed: the sync in fix dpd/energy/kk then had nothing to copy and the half step integrated the pre-communication values instead. Found with the sync-debugging build, which reports [watch] pair:duCond: the host side was written without a claim and this sync_device has nothing to copy -- the device keeps stale data [stale] pair:duCond: device side read while host side is newer, from FixDPDenergyKokkos::take_half_step() and shows up as diverging thermo output from the tenth step on in examples/PACKAGES/dpd-react (dpde-vv, dpde-shardlow, dpdh-shardlow). (cherry picked from commit 65b2d13311b6a98f1ea6c13c74c6d5aa7eb26f2d) --- src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp b/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp index 2ba9e1f79e8..d4f194b20cf 100644 --- a/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp +++ b/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp @@ -328,6 +328,14 @@ void PairDPDfdtEnergyKokkos::compute(int eflag_in, int vflag_in) k_duMech.template modify(); k_duMech.sync_host(); comm->reverse_comm(this); + + // the reverse communication adds the ghost contributions through the plain + // host arrays, so claim that write: without it the device copy keeps the + // pre-communication values and the sync in fix dpd/energy/kk has nothing + // left to copy + + k_duCond.modify_host(); + k_duMech.modify_host(); } if (eflag_global) eng_vdwl += ev.evdwl; From e47236c9ac67885cd93964ef964f3ea8117800a0 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 02:00:52 +0000 Subject: [PATCH 02/43] KOKKOS: re-take the atom data and views around energy_force in min fire MinFireKokkos::run_iterate() synced the per-atom arrays to the device once, took the device views and nlocal, and then ran the whole iteration loop on them. energy_force() inside that loop runs the communication, the neighbor build and the fixes, so it can leave the newest data on the host and can grow or reorder the arrays. The integration kernels then worked from a device copy the host had overtaken, and the modified() after them found the host side claimed as well; nlocal could also be out of date after a migration. The sync-debugging build stops the run at that point: LAMMPS::DualView::modify_device ERROR: concurrent modification of host and device views in DualView "atom:x" LAMMPS_NS::MinFireKokkos::run_iterate<0, false>(int) reproduced by examples/PACKAGES/pafi, where fix pafi has no KOKKOS support and does its work on the host. (cherry picked from commit f8d10ee8b01c10e3649e9d5ee6e35ae5ed27e0f5) --- src/KOKKOS/min_fire_kokkos.cpp | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/KOKKOS/min_fire_kokkos.cpp b/src/KOKKOS/min_fire_kokkos.cpp index 4ab31ed9202..1c0e4a43012 100644 --- a/src/KOKKOS/min_fire_kokkos.cpp +++ b/src/KOKKOS/min_fire_kokkos.cpp @@ -120,9 +120,28 @@ int MinFireKokkos::run_iterate(int maxiter) { auto l_type = atomKK->k_type.view_device(); int nlocal = atom->nlocal; + // energy_force() runs the communication, the neighbor build and the fixes. + // Any of those can leave the newest per-atom data on the host -- a fix + // without KOKKOS support writes there -- and can grow or reorder the arrays, + // so the data and the views both have to be taken again after every call. + // Without this the kernels below integrate a device copy that the host has + // since overtaken, and the modified() that follows finds both sides claimed. + + auto refresh = [&]() { + atomKK->sync(Device, X_MASK | V_MASK | F_MASK | RMASS_MASK | TYPE_MASK); + l_x = atomKK->k_x.view_device(); + l_v = atomKK->k_v.view_device(); + l_f = atomKK->k_f.view_device(); + l_rmass = atomKK->k_rmass.view_device(); + l_mass = atomKK->k_mass.view_device(); + l_type = atomKK->k_type.view_device(); + nlocal = atom->nlocal; + }; + if constexpr (INTEGRATOR == LEAPFROG) { energy_force(0); neval++; + refresh(); double dtf = -0.5 * dt * force->ftm2v; Kokkos::parallel_for("min_fire/leapfrog_init", atom->nlocal, LAMMPS_LAMBDA(const int i) { KK_FLOAT dtfm = dtf / (l_rmass.data() ? l_rmass(i) : l_mass(l_type(i))); @@ -135,6 +154,8 @@ int MinFireKokkos::run_iterate(int maxiter) { for (int iter = 0; iter < maxiter; iter++) { if (timer->check_timeout(niter)) return TIMEOUT; + refresh(); + bigint ntimestep = ++update->ntimestep; niter++; @@ -215,6 +236,7 @@ int MinFireKokkos::run_iterate(int maxiter) { if (!ABCFLAG && flagv0) { energy_force(0); neval++; + refresh(); double dtf_init = dt * force->ftm2v; Kokkos::parallel_for("min_fire/v_init", nlocal, LAMMPS_LAMBDA(const int i) { KK_FLOAT dtfm = dtf_init / (l_rmass.data() ? l_rmass(i) : l_mass(l_type(i))); @@ -315,9 +337,9 @@ int MinFireKokkos::run_iterate(int maxiter) { eprevious = ecurrent; ecurrent = energy_force(0); neval++; + refresh(); if constexpr (INTEGRATOR == VERLET) { - atomKK->sync(Device, V_MASK | F_MASK); Kokkos::parallel_for("min_fire/verlet_v_final", nlocal, LAMMPS_LAMBDA(const int i) { KK_FLOAT dtfm_half = dtf_half / (l_rmass.data() ? l_rmass(i) : l_mass(l_type(i))); l_v(i,0) += dtfm_half * l_f(i,0); From 0a6056c421b42b2eebb81937d37929da896aefe1 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 02:10:36 +0000 Subject: [PATCH 03/43] KOKKOS: release the host side of the SPIN forces in force_clear force_clear() overwrites atom->fm and atom->fm_long on the device and claims them there, the same way it does for f and torque, but only f and torque had their host side released first. A host side left claimed by the setup was therefore still claimed when the device claim came, giving two claimed sides with nothing to say which one is current. Every input under examples/SPIN stops in the sync-debugging build with LAMMPS::DualView::modify_device ERROR: concurrent modification of host and device views in DualView "atom:fm" LAMMPS_NS::VerletKokkos::force_clear() (cherry picked from commit d50b393313d6afacc8af4d97d88254f5768ca255) --- src/KOKKOS/verlet_kokkos.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index fbabd6c2bb0..a1fe8e38a3f 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -577,6 +577,16 @@ void VerletKokkos::force_clear() atomKK->k_f.clear_sync_state(); // ignore host forces/torques since device views atomKK->k_torque.clear_sync_state(); // will be cleared below + // the SPIN forces below are overwritten in the same way, so their host side + // has to be released here as well -- without this a host side left claimed + // from the setup is still claimed when the device side is claimed below, and + // the two disagree with nothing to say which one is current + + if (extraflag) { + atomKK->k_fm.clear_sync_state(); + atomKK->k_fm_long.clear_sync_state(); + } + // clear force on all particles // if either newton flag is set, also include ghosts // when using threads always clear all forces. From bb93230e515442db90edb722c0081e23db51269a Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 03:20:21 +0000 Subject: [PATCH 04/43] KOKKOS: carry the comm buffer claim fixes into the tiled twin comm_kokkos.cpp drops the claim that resizing a send or scratch buffer leaves behind, because those buffers are filled through raw pointers on whichever side packs them and the claim otherwise stands until something claims the other side and the two collide. comm_tiled_kokkos.cpp has the same two functions and did not get the same treatment. examples/balance/in.balance.neigh.rcb, which is what selects the tiled communication, stops in the sync-debugging build with LAMMPS::DualView::modify_host ERROR: concurrent modification of host and device views LAMMPS_NS::CommTiledKokkos::grow_send_kokkos(int, int, ExecutionSpace) (cherry picked from commit 0edfa539df67a5557b00c5bb8e7ad4681c4571db) --- src/KOKKOS/comm_tiled_kokkos.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/KOKKOS/comm_tiled_kokkos.cpp b/src/KOKKOS/comm_tiled_kokkos.cpp index abfc1884bfe..62ecde2b7fa 100644 --- a/src/KOKKOS/comm_tiled_kokkos.cpp +++ b/src/KOKKOS/comm_tiled_kokkos.cpp @@ -668,6 +668,12 @@ void CommTiledKokkos::grow_buf_pair(int n) max_buf_pair = n * BUFFACTOR; k_buf_send_pair.resize(max_buf_pair); k_buf_recv_pair.resize(max_buf_pair); + + // resizing claims a side; these are scratch buffers that are filled + // before they are read, so drop the claim rather than leave it for the + // next modify_host() to collide with + k_buf_send_pair.clear_sync_state(); + k_buf_recv_pair.clear_sync_state(); } /* ---------------------------------------------------------------------- @@ -855,6 +861,12 @@ void CommTiledKokkos::grow_send_kokkos(int n, int flag, ExecutionSpace space) atomKK->avecKK->size_border + atomKK->avecKK->size_velocity); else k_buf_send.resize(maxsend_border,atomKK->avecKK->size_border); + + // the claim above only steers the resize to the side whose contents have + // to survive; after it this is a scratch buffer again, filled through raw + // pointers on whichever side does the packing, so drop the claim rather + // than leave it standing forever + k_buf_send.clear_sync_state(); } else { if (ghost_velocity) MemoryKokkos::realloc_kokkos(k_buf_send,"comm:k_buf_send",maxsend_border, From b65be47205b8eb019fc35763fceae7385b323ce0 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 03:58:50 +0000 Subject: [PATCH 05/43] KOKKOS: claim the dpd/fdt/energy forces after the kernels write them pair_style dpd/fdt/energy/kk claimed the force array on its execution space at the top of compute(), before the kernels that fill it. Between the two it copies the energy changes to the host and runs their reverse communication, and that copies the forces to the host as well and takes the claim with it. The forces the kernels then wrote were left unclaimed, and because this style declares an empty datamask there is no second claim from the integrator to cover them: the reverse communication of the forces afterwards found nothing to copy and summed the ghost contributions into the previous step's values. [watch] atom:f: the device side was written without a claim and this sync_host has nothing to copy -- the host keeps stale data LAMMPS_NS::CommKokkos::reverse_comm() LAMMPS_NS::VerletKokkos::setup(int) Claim them after the kernels instead, as pair_style dpd/kk does. With this and the reverse-communication claim of the energy changes, the three examples/PACKAGES/dpd-react inputs reproduce the reference build exactly; before, they parted from it at the tenth step. (cherry picked from commit a55c9cd163f7e37b8841229f81fa302682b9c0a3) --- src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp b/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp index d4f194b20cf..9fc8e3943c6 100644 --- a/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp +++ b/src/KOKKOS/pair_dpd_fdt_energy_kokkos.cpp @@ -164,8 +164,6 @@ void PairDPDfdtEnergyKokkos::compute(int eflag_in, int vflag_in) k_cutsq.template sync(); k_params.template sync(); atomKK->sync(execution_space,X_MASK | F_MASK | TYPE_MASK | ENERGY_MASK | VIRIAL_MASK); - if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK); - else atomKK->modified(execution_space,F_MASK); special_lj[0] = force->special_lj[0]; special_lj[1] = force->special_lj[1]; @@ -338,6 +336,15 @@ void PairDPDfdtEnergyKokkos::compute(int eflag_in, int vflag_in) k_duMech.modify_host(); } + // claim the forces here rather than before the kernels above: this style + // declares an empty datamask, so this is the only claim they get, and the + // reverse communication of the energy changes syncs the host in between, + // which would take a claim made up there with it and leave the forces the + // kernels wrote unclaimed + + if (evflag) atomKK->modified(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK); + else atomKK->modified(execution_space,F_MASK); + if (eflag_global) eng_vdwl += ev.evdwl; if (vflag_global) { virial[0] += ev.v[0]; From 8e870a5cee84ba1f6fbe7d0efc291a1917e17c4b Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 05:26:14 +0000 Subject: [PATCH 06/43] KOKKOS: claim the per-type masses when they are set Atom::set_mass() writes the per-type masses through the plain host array, which leaves the device copy holding whatever it was allocated with and nothing in the state to say so. Two styles worked around it by claiming and copying the masses in their own init() -- fix nve/kk and fix nh/kk -- and the other fourteen KOKKOS styles that read the masses on the device relied on one of those two being in the input. examples/SPIN has neither: fix nve/spin has no KOKKOS version, so nobody claimed the write and compute temp/kk read masses that were never copied. The sync-debugging build reports [watch] atom::mass: the host side was written without a claim and this sync_device has nothing to copy -- the device keeps stale data [stale] atom::mass: device side read while host side is newer, from ComputeTempKokkos::compute_scalar() and all twelve examples/SPIN inputs come out different from the reference build; running them with the two copies of the masses forced into one reproduces it exactly, which is what identifies the array. Claim the write where it happens instead, by making the four set_mass() overloads virtual and overriding them in AtomKokkos. The two init() workarounds are left as they are; they are harmless now. (cherry picked from commit 0dab7fae4197e2d6c75f7ec46e04fcd0866d7d08) --- src/KOKKOS/atom_kokkos.cpp | 30 ++++++++++++++++++++++++++++++ src/KOKKOS/atom_kokkos.h | 10 ++++++++++ src/atom.h | 8 ++++---- 3 files changed, 44 insertions(+), 4 deletions(-) diff --git a/src/KOKKOS/atom_kokkos.cpp b/src/KOKKOS/atom_kokkos.cpp index 993d4aef431..67b55d5418b 100644 --- a/src/KOKKOS/atom_kokkos.cpp +++ b/src/KOKKOS/atom_kokkos.cpp @@ -275,6 +275,36 @@ void AtomKokkos::sync_pinned(const ExecutionSpace space, uint64_t mask, int asyn avecKK->sync_pinned(space, mask, async_flag); for (int n = 0; n < nprop_atom; n++) fix_prop_atom[n]->sync_pinned(space, mask, async_flag); } +/* ---------------------------------------------------------------------- + the four ways of setting the per-type masses all write the plain host + array, so claim that write for the device copy +------------------------------------------------------------------------- */ + +void AtomKokkos::set_mass(const char *file, int line, const char *str, + int type_offset, int labelflag, int *ilabel) +{ + Atom::set_mass(file, line, str, type_offset, labelflag, ilabel); + k_mass.modify_host(); +} + +void AtomKokkos::set_mass(const char *file, int line, int itype, double value) +{ + Atom::set_mass(file, line, itype, value); + k_mass.modify_host(); +} + +void AtomKokkos::set_mass(const char *file, int line, int narg, char **arg) +{ + Atom::set_mass(file, line, narg, arg); + k_mass.modify_host(); +} + +void AtomKokkos::set_mass(double *values) +{ + Atom::set_mass(values); + k_mass.modify_host(); +} + /* ---------------------------------------------------------------------- */ void AtomKokkos::allocate_type_arrays() diff --git a/src/KOKKOS/atom_kokkos.h b/src/KOKKOS/atom_kokkos.h index 981b8dee91f..50336e1101b 100644 --- a/src/KOKKOS/atom_kokkos.h +++ b/src/KOKKOS/atom_kokkos.h @@ -189,6 +189,16 @@ class AtomKokkos : public Atom { void init() override; void update_property_atom(); void allocate_type_arrays() override; + + // the per-type masses are written through the plain host array, which leaves + // the device copy behind with nothing to say so. Claim the write here, at + // the one place all four spellings of the mass command go through, rather + // than in each of the styles that read the masses on the device. + + void set_mass(const char *, int, const char *, int, int, int *) override; + void set_mass(const char *, int, int, double) override; + void set_mass(const char *, int, int, char **) override; + void set_mass(double *) override; void *extract(const char *) override; void sync(const ExecutionSpace space, uint64_t mask); void modified(const ExecutionSpace space, uint64_t mask); diff --git a/src/atom.h b/src/atom.h index 29418989cb5..a6a7dcac43b 100644 --- a/src/atom.h +++ b/src/atom.h @@ -354,10 +354,10 @@ class Atom : protected Pointers { void data_fix_compute_variable(int, int); virtual void allocate_type_arrays(); - void set_mass(const char *, int, const char *, int, int, int *); - void set_mass(const char *, int, int, double); - void set_mass(const char *, int, int, char **); - void set_mass(double *); + virtual void set_mass(const char *, int, const char *, int, int, int *); + virtual void set_mass(const char *, int, int, double); + virtual void set_mass(const char *, int, int, char **); + virtual void set_mass(double *); void check_mass(const char *, int); int radius_consistency(int, double &); From d76cde0bdcd3cb6d313c846a96cea8f6dad9189a Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 05:47:02 +0000 Subject: [PATCH 07/43] KOKKOS: bracket only the host step of fix deform's pre_exchange FixDeformKokkos::pre_exchange() wrapped the whole base class call in a host sync and a host claim. Everything that call does to the per-atom data goes through DomainKokkos -- image_flip(), remap_all(), x2lamda() and lamda2x() -- which runs on the device and declares itself, so the wrapper claimed the host side of arrays the device work had just claimed, and the run stopped: LAMMPS::DualView::modify_host ERROR: concurrent modification of host and device views in DualView "atom:x" LAMMPS_NS::ModifyKokkos::pre_exchange() reproduced by examples/VISCOSITY/in.nemd.2d at the first box flip. The host sync at the front was also too early to be of use to the one step that needs it, the atom migration, which runs after two device kernels have moved the coordinates on. Give the migration its own virtual hook in FixDeform and bracket that instead, leaving the device work to declare itself as it already does. (cherry picked from commit 7303197d18a5d0e6786b0099fe5a589082fcdf5a) --- src/KOKKOS/fix_deform_kokkos.cpp | 17 ++++++++++++++++- src/KOKKOS/fix_deform_kokkos.h | 1 + src/fix_deform.cpp | 9 ++++++++- src/fix_deform.h | 4 ++++ 4 files changed, 29 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/fix_deform_kokkos.cpp b/src/KOKKOS/fix_deform_kokkos.cpp index 5494748f3a3..b0c5c271dd8 100644 --- a/src/KOKKOS/fix_deform_kokkos.cpp +++ b/src/KOKKOS/fix_deform_kokkos.cpp @@ -46,8 +46,23 @@ FixDeformKokkos::FixDeformKokkos(LAMMPS *lmp, int narg, char **arg) : FixDeform( void FixDeformKokkos::pre_exchange() { - atomKK->sync(Host,ALL_MASK); + // everything the base class does here goes through DomainKokkos, which runs + // on the device and declares itself, except the atom migration below. + // Bracketing the whole call for the host claimed a side the device work had + // already claimed, and left the migration reading host data that the device + // work had since overtaken. + FixDeform::pre_exchange(); +} + +/* ---------------------------------------------------------------------- + the migration is the one host only step in pre_exchange() +------------------------------------------------------------------------- */ + +void FixDeformKokkos::migrate_atoms() +{ + atomKK->sync(Host,ALL_MASK); + FixDeform::migrate_atoms(); atomKK->modified(Host,ALL_MASK); } diff --git a/src/KOKKOS/fix_deform_kokkos.h b/src/KOKKOS/fix_deform_kokkos.h index a306a291ef4..0c30f85ffe2 100644 --- a/src/KOKKOS/fix_deform_kokkos.h +++ b/src/KOKKOS/fix_deform_kokkos.h @@ -32,6 +32,7 @@ class FixDeformKokkos : public FixDeform { FixDeformKokkos(class LAMMPS *, int, char **); void pre_exchange() override; + void migrate_atoms() override; void update_box() override; }; diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index c95c57b9a44..259d283d4e8 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -778,6 +778,13 @@ void FixDeform::init() image flags to new values, making eqs in doc of Domain:image_flip incorrect ------------------------------------------------------------------------- */ +void FixDeform::migrate_atoms() +{ + irregular->migrate_atoms(); +} + +/* ---------------------------------------------------------------------- */ + void FixDeform::pre_exchange() { if (flip == 0) return; @@ -821,7 +828,7 @@ void FixDeform::pre_exchange() domain->remap_all(); domain->x2lamda(atom->nlocal); - irregular->migrate_atoms(); + migrate_atoms(); domain->lamda2x(atom->nlocal); flip = 0; diff --git a/src/fix_deform.h b/src/fix_deform.h index f36f4d85424..d79972ca75c 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -37,6 +37,10 @@ class FixDeform : public Fix { int setmask() override; void init() override; void pre_exchange() override; + + // the atom migration is a separate step so that an accelerator version can + // move the per-atom data to the side this runs on and back + virtual void migrate_atoms(); void end_of_step() override; void write_restart(FILE *) override; void restart(char *buf) override; From 2460560dde0103beaa46f7a91f24f307a332d46b Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 05:50:24 +0000 Subject: [PATCH 08/43] KOKKOS: claim the coordinates and velocities min fire's kernels write Four of the kernels in MinFireKokkos::run_iterate() wrote the per-atom coordinates or velocities on the device and never declared it. The one that matters most is the inertia reset, which steps the coordinates back by half a step and is followed immediately by energy_force(): the forward communication there syncs the host, finds nothing to copy and builds the ghost positions from coordinates the device had already moved on from. [watch] atom:x: the device side was written without a claim and this sync_host has nothing to copy -- the host keeps stale data LAMMPS_NS::CommKokkos::forward_comm(int) LAMMPS_NS::MinKokkos::energy_force(int) LAMMPS_NS::MinFireKokkos::run_iterate<0, false>(int) With this and the refresh around energy_force, examples/fire (in.fire, in.fire_mod, in.meam.fire) and examples/PACKAGES/pafi reproduce the reference build exactly. (cherry picked from commit 76605431e70072ae17ba0cbeeafd4185660dc402) --- src/KOKKOS/min_fire_kokkos.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/KOKKOS/min_fire_kokkos.cpp b/src/KOKKOS/min_fire_kokkos.cpp index 1c0e4a43012..8ddc685a5aa 100644 --- a/src/KOKKOS/min_fire_kokkos.cpp +++ b/src/KOKKOS/min_fire_kokkos.cpp @@ -149,6 +149,7 @@ int MinFireKokkos::run_iterate(int maxiter) { l_v(i,1) = dtfm * l_f(i,1); l_v(i,2) = dtfm * l_f(i,2); }); + atomKK->modified(Device, V_MASK); } for (int iter = 0; iter < maxiter; iter++) { @@ -230,6 +231,7 @@ int MinFireKokkos::run_iterate(int maxiter) { } l_v(i,0) = l_v(i,1) = l_v(i,2) = 0.0; }); + atomKK->modified(Device, X_MASK | V_MASK); flagv0 = 1; } @@ -244,6 +246,7 @@ int MinFireKokkos::run_iterate(int maxiter) { l_v(i,1) = dtfm * l_f(i,1); l_v(i,2) = dtfm * l_f(i,2); }); + atomKK->modified(Device, V_MASK); } // cannot use "if constexpr" below because CUDA device lambdas @@ -268,6 +271,7 @@ int MinFireKokkos::run_iterate(int maxiter) { Kokkos::parallel_for("min_fire/final_v_zero", nlocal, LAMMPS_LAMBDA(const int i) { l_v(i,0) = l_v(i,1) = l_v(i,2) = 0.0; }); + atomKK->modified(Device, V_MASK); } KK_FLOAT dtf_final = dtv * force->ftm2v; From eaca3f045ff76c98ca9b9ea512f6c813506bade0 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 06:30:53 +0000 Subject: [PATCH 09/43] KOKKOS: grow the shake clusters from the side the plain pointers use FixShakeKokkos::grow_arrays() synced the cluster tables to the device before growing them. grow_kokkos() grows the host side and hands the plain shake_flag, shake_atom and shake_type pointers back to the base class, which reads and writes the clusters through them in copy_arrays(), the exchange packers and the destructor, so it is the host side that has to be current there. Syncing the other way left the grown host arrays holding what they had before the last device update. [stale] shake:shake_flag: host side read while device side is newer, from FixShakeKokkos::grow_arrays(int) examples/rdf-adf/in.spce and in.spce.hbond stopped at the first step with 'Out of range atoms - cannot compute PPPM' and a pressure of -7e+300; both now reproduce the reference build exactly, and examples/peptide and examples/micelle, which also use fix shake, are unchanged. (cherry picked from commit a3a56dc2899ad4df2c9ac6f67a8f9d906c80dcc1) --- src/KOKKOS/fix_shake_kokkos.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index 13ee3341cc2..4c2a4e19866 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -1617,9 +1617,13 @@ void FixShakeKokkos::stats() template void FixShakeKokkos::grow_arrays(int nmax) { - k_shake_flag.sync_device(); - k_shake_atom.sync_device(); - k_shake_type.sync_device(); + // grow_kokkos() grows the host side and hands the plain pointers below it + // back to the base class, which reads and writes the clusters through them, + // so the host side is the one that has to be current here + + k_shake_flag.sync_host(); + k_shake_atom.sync_host(); + k_shake_type.sync_host(); memoryKK->grow_kokkos(k_shake_flag,shake_flag,nmax,"shake:shake_flag"); memoryKK->grow_kokkos(k_shake_atom,shake_atom,nmax,4,"shake:shake_atom"); From 6289387d5659a3fdc60e323708e506e8cfd2d7f6 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 14:19:39 +0000 Subject: [PATCH 10/43] Refuse run_style respa with the KOKKOS package, and claim nve_x's write rRESPA keeps its own copy of the forces of each level and clears and sums them through the plain LAMMPS arrays, and it calls the force computations directly, without the transfers between the host and the device that run_style verlet does for KOKKOS. There is no rRESPA version in the package to supply them. Where the two sides have separate memory the forces are simply wrong, with nothing to say so: [watch] atom:f: the device side was written without a claim and this sync_host has nothing to copy -- the host keeps stale data LAMMPS_NS::CommKokkos::reverse_comm() LAMMPS_NS::Respa::setup(int) examples/relres/in.22DMH.respa ran to completion with the temperature at 2930 K instead of 292 K by the fiftieth step. Refuse the combination. Along the way, fix nvt/kk claimed the coordinates before the kernel that writes them rather than after, so any copy running in between took the claim and left the new coordinates unclaimed -- the same shape as the dpd/fdt/energy claim. Its two siblings, nve_v() and nh_v_press(), claim after their kernels already. (cherry picked from commit 52185ff254a265390b6eed1d83240cee38a8d703) --- doc/src/run_style.rst | 9 +++++++++ src/KOKKOS/fix_nh_kokkos.cpp | 7 ++++++- src/respa.cpp | 10 ++++++++++ 3 files changed, 25 insertions(+), 1 deletion(-) diff --git a/doc/src/run_style.rst b/doc/src/run_style.rst index af43a9b244d..f95e8928901 100644 --- a/doc/src/run_style.rst +++ b/doc/src/run_style.rst @@ -375,6 +375,15 @@ REPLICA package. Correspondingly the *respa/omp* style is available only if the OPENMP package was included. See the :doc:`Build package ` page for more info. +.. versionchanged:: TBD + +The *respa* style cannot be used together with the KOKKOS package. It +keeps a copy of the forces of each level and clears and sums those +through the plain LAMMPS arrays, and it calls the force computations +without the transfers between the host and the device that the *verlet* +style performs for KOKKOS. On a GPU this gives wrong forces with no +indication that anything is amiss, so the combination is refused. + Run style *verlet/split* is not compatible with kspace styles from the INTEL package and it is not compatible with any tip4p, dipole, or spin kspace styles. diff --git a/src/KOKKOS/fix_nh_kokkos.cpp b/src/KOKKOS/fix_nh_kokkos.cpp index 1f6fee0014f..b7228570aea 100644 --- a/src/KOKKOS/fix_nh_kokkos.cpp +++ b/src/KOKKOS/fix_nh_kokkos.cpp @@ -613,7 +613,6 @@ template void FixNHKokkos::nve_x() { atomKK->sync(execution_space,X_MASK | V_MASK | MASK_MASK); - atomKK->modified(execution_space,X_MASK); x = atomKK->k_x.view(); v = atomKK->k_v.view(); @@ -626,6 +625,12 @@ void FixNHKokkos::nve_x() copymode = 1; Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); copymode = 0; + + // claim the coordinates after the kernel has written them, not before: a + // claim made up front is taken by any copy that runs in between, and under + // rRESPA one does, which leaves the new coordinates on the device unclaimed + + atomKK->modified(execution_space,X_MASK); } template diff --git a/src/respa.cpp b/src/respa.cpp index 0b08799d9b6..1e920f5891b 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -17,6 +17,7 @@ #include "respa.h" +#include "accelerator_kokkos.h" #include "angle.h" #include "atom.h" #include "atom_vec.h" @@ -288,6 +289,15 @@ void Respa::init() { Integrate::init(); + // rRESPA keeps its own copy of the forces of each level and clears and sums + // them through the plain LAMMPS arrays, and it calls the force computations + // without the transfers between the host and the device that run_style + // verlet does for the KOKKOS package. On a GPU that silently gives the + // wrong forces, so refuse the combination rather than run it. + + if (lmp->kokkos) + error->all(FLERR, "Run style respa is not supported by the KOKKOS package"); + // warn if no fixes if (modify->nfix == 0 && comm->me == 0) From 7c8efe6dcf132810b8f7bed578832b5847deb393 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 15:03:00 +0000 Subject: [PATCH 11/43] Let a style declare a host write of the per-atom arrays, and use it fix numdiff, fix numdiff/virial and compute born/matrix numdiff all take a finite difference by displacing the atoms through the plain coordinate array and then calling the force computations directly. With the KOKKOS package those computations work from their own copies, and the copy that brings the host side up to date runs after the displacement was written, so it overwrites it: every difference was taken between two evaluations of the undisplaced configuration. The same happens to the forces, which these styles clear, save and restore through the plain array. [watch] atom:x: the host side was written, never claimed, and is now lost; the write is between modify_device and sync_host, which discards it element 0 of 49152 changed from 0 to 0.0001 examples/numdiff/in.numdiff reported a relative force error of 2.5e-03 where the reference build gives 5.6e-09, and the restored coordinates were off by the displacement, so the trajectory drifted and the total energy stopped being conserved. Add Atom::sync_host_arrays() and Atom::modified_host_arrays(), which do nothing without the package and which AtomKokkos maps onto its existing host transfers, and bracket the six places where these three styles read or write the arrays directly. This avoids a KOKKOS version of each of them, and the same two calls are available to any other style that reaches for the plain arrays. The input now reproduces the reference build exactly. (cherry picked from commit cb16ae12875c7268e73b3891c009993668e060d3) --- src/EXTRA-COMPUTE/compute_born_matrix.cpp | 27 +++++++++++++++++++++++ src/EXTRA-FIX/fix_numdiff.cpp | 24 ++++++++++++++++++++ src/EXTRA-FIX/fix_numdiff_virial.cpp | 19 ++++++++++++++++ src/KOKKOS/atom_kokkos.h | 3 +++ src/atom.h | 10 +++++++++ 5 files changed, 83 insertions(+) diff --git a/src/EXTRA-COMPUTE/compute_born_matrix.cpp b/src/EXTRA-COMPUTE/compute_born_matrix.cpp index 00e17923b57..2787fbe8000 100644 --- a/src/EXTRA-COMPUTE/compute_born_matrix.cpp +++ b/src/EXTRA-COMPUTE/compute_born_matrix.cpp @@ -19,6 +19,7 @@ #include "angle.h" #include "atom.h" +#include "atom_masks.h" #include "atom_vec.h" #include "bond.h" #include "comm.h" @@ -480,6 +481,8 @@ void ComputeBornMatrix::compute_numdiff() // store copy of current forces for owned and ghost atoms + atom->sync_host_arrays(X_MASK | F_MASK); + double **x = atom->x; double **f = atom->f; @@ -533,8 +536,12 @@ void ComputeBornMatrix::compute_numdiff() // restore original forces for owned and ghost atoms + atom->sync_host_arrays(F_MASK); + for (int i = 0; i < nall; i++) for (int k = 0; k < 3; k++) f[i][k] = temp_f[i][k]; + + atom->modified_host_arrays(F_MASK); } /* ---------------------------------------------------------------------- @@ -543,6 +550,12 @@ void ComputeBornMatrix::compute_numdiff() void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) { + // the strain goes into the plain coordinate array, and the energy evaluation + // that follows works from the KOKKOS copies, so bring the host side up to + // date first and hand the write over afterwards + + atom->sync_host_arrays(X_MASK); + double **x = atom->x; // NOTE: virial_addon() expressions predicated on @@ -563,6 +576,8 @@ void ComputeBornMatrix::displace_atoms(int nall, int idir, double magnitude) x[i][k] = temp_x[i][k] + 0.5 * numdelta * magnitude * (temp_x[i][l] - fixedpoint[l]); x[i][l] = temp_x[i][l] + 0.5 * numdelta * magnitude * (temp_x[i][k] - fixedpoint[k]); } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -576,6 +591,9 @@ void ComputeBornMatrix::restore_atoms(int nall, int idir) int k = dirlist[idir][0]; int l = dirlist[idir][1]; + + atom->sync_host_arrays(X_MASK); + double **x = atom->x; if (l == k) for (int i = 0; i < nall; i++) x[i][k] = temp_x[i][k]; @@ -584,6 +602,8 @@ void ComputeBornMatrix::restore_atoms(int nall, int idir) x[i][l] = temp_x[i][l]; x[i][k] = temp_x[i][k]; } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -673,9 +693,16 @@ void ComputeBornMatrix::virial_addon() void ComputeBornMatrix::force_clear(int nall) { + // the forces are cleared through the plain array and accumulated again by + // the force computations, which work from the KOKKOS copies + + atom->sync_host_arrays(F_MASK); + double **forces = atom->f; size_t nbytes = 3 * sizeof(double) * nall; if (nbytes) memset(&forces[0][0], 0, nbytes); + + atom->modified_host_arrays(F_MASK); } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-FIX/fix_numdiff.cpp b/src/EXTRA-FIX/fix_numdiff.cpp index fbebe159e61..acd7e222e24 100644 --- a/src/EXTRA-FIX/fix_numdiff.cpp +++ b/src/EXTRA-FIX/fix_numdiff.cpp @@ -19,6 +19,7 @@ #include "angle.h" #include "atom.h" +#include "atom_masks.h" #include "bond.h" #include "compute.h" #include "dihedral.h" @@ -191,6 +192,8 @@ void FixNumDiff::calculate_forces() // store copy of current forces for owned and ghost atoms + atom->sync_host_arrays(X_MASK | F_MASK); + double **x = atom->x; double **f = atom->f; int nlocal = atom->nlocal; @@ -245,8 +248,12 @@ void FixNumDiff::calculate_forces() // restore original forces for owned and ghost atoms + atom->sync_host_arrays(F_MASK); + for (i = 0; i < nall; i++) for (j = 0; j < 3; j++) { f[i][j] = temp_f[i][j]; } + + atom->modified_host_arrays(F_MASK); } /* ---------------------------------------------------------------------- @@ -257,6 +264,12 @@ void FixNumDiff::displace_atoms(int ilocal, int idim, int magnitude) { if (ilocal < 0) return; + // the displacement goes into the plain coordinate array, and the energy + // evaluation that follows works from the KOKKOS copies, so bring the host + // side up to date first and hand the write over afterwards + + atom->sync_host_arrays(X_MASK); + double **x = atom->x; int *sametag = atom->sametag; int j = ilocal; @@ -266,6 +279,8 @@ void FixNumDiff::displace_atoms(int ilocal, int idim, int magnitude) j = sametag[j]; x[j][idim] += delta * magnitude; } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -276,6 +291,8 @@ void FixNumDiff::restore_atoms(int ilocal, int idim) { if (ilocal < 0) return; + atom->sync_host_arrays(X_MASK); + double **x = atom->x; int *sametag = atom->sametag; int j = ilocal; @@ -285,6 +302,8 @@ void FixNumDiff::restore_atoms(int ilocal, int idim) j = sametag[j]; x[j][idim] = temp_x[j][idim]; } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -294,7 +313,12 @@ void FixNumDiff::restore_atoms(int ilocal, int idim) double FixNumDiff::update_energy() { + // the forces are cleared through the plain array and accumulated again by + // the force computations, which work from the KOKKOS copies + + atom->sync_host_arrays(F_MASK); force_clear(atom->f); + atom->modified_host_arrays(F_MASK); // flag that we only need to compute the global energy int eflag = ENERGY_GLOBAL | ENERGY_ONLY; diff --git a/src/EXTRA-FIX/fix_numdiff_virial.cpp b/src/EXTRA-FIX/fix_numdiff_virial.cpp index beff4b616d2..e6b8ba14ef4 100644 --- a/src/EXTRA-FIX/fix_numdiff_virial.cpp +++ b/src/EXTRA-FIX/fix_numdiff_virial.cpp @@ -19,6 +19,7 @@ #include "angle.h" #include "atom.h" +#include "atom_masks.h" #include "bond.h" #include "compute.h" #include "dihedral.h" @@ -201,6 +202,8 @@ void FixNumDiffVirial::calculate_virial() // store copy of current forces for owned and ghost atoms + atom->sync_host_arrays(X_MASK | F_MASK); + double **x = atom->x; double **f = atom->f; int nall = atom->nlocal + atom->nghost; @@ -237,8 +240,12 @@ void FixNumDiffVirial::calculate_virial() // restore original forces for owned and ghost atoms + atom->sync_host_arrays(F_MASK); + for (int i = 0; i < nall; i++) for (int k = 0; k < 3; k++) f[i][k] = temp_f[i][k]; + + atom->modified_host_arrays(F_MASK); } /* ---------------------------------------------------------------------- @@ -247,11 +254,19 @@ void FixNumDiffVirial::calculate_virial() void FixNumDiffVirial::displace_atoms(int nall, int idir, double magnitude) { + // the strain goes into the plain coordinate array, and the energy evaluation + // that follows works from the KOKKOS copies, so bring the host side up to + // date first and hand the write over afterwards + + atom->sync_host_arrays(X_MASK); + double **x = atom->x; int k = dirlist[idir][0]; int l = dirlist[idir][1]; for (int i = 0; i < nall; i++) x[i][k] = temp_x[i][k] + delta * magnitude * (temp_x[i][l] - fixedpoint[l]); + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- @@ -260,9 +275,13 @@ void FixNumDiffVirial::displace_atoms(int nall, int idir, double magnitude) void FixNumDiffVirial::restore_atoms(int nall, int idir) { + atom->sync_host_arrays(X_MASK); + double **x = atom->x; int k = dirlist[idir][0]; for (int i = 0; i < nall; i++) { x[i][k] = temp_x[i][k]; } + + atom->modified_host_arrays(X_MASK); } /* ---------------------------------------------------------------------- diff --git a/src/KOKKOS/atom_kokkos.h b/src/KOKKOS/atom_kokkos.h index 50336e1101b..ddd0a508bec 100644 --- a/src/KOKKOS/atom_kokkos.h +++ b/src/KOKKOS/atom_kokkos.h @@ -190,6 +190,9 @@ class AtomKokkos : public Atom { void update_property_atom(); void allocate_type_arrays() override; + void sync_host_arrays(uint64_t mask) override { sync(Host, mask); } + void modified_host_arrays(uint64_t mask) override { modified(Host, mask); } + // the per-type masses are written through the plain host array, which leaves // the device copy behind with nothing to say so. Claim the write here, at // the one place all four spellings of the mass command go through, rather diff --git a/src/atom.h b/src/atom.h index a6a7dcac43b..f4da148e86d 100644 --- a/src/atom.h +++ b/src/atom.h @@ -354,6 +354,16 @@ class Atom : protected Pointers { void data_fix_compute_variable(int, int); virtual void allocate_type_arrays(); + + // A style that reads or writes the plain per-atom arrays on the host says so + // with these, so that the KOKKOS package can bring that side up to date + // first and can carry the write over to the device afterwards. They do + // nothing without the package. Use them where a style touches the arrays + // directly and then calls something -- a force computation, the + // communication -- that works from the KOKKOS copies. + + virtual void sync_host_arrays(uint64_t) {} + virtual void modified_host_arrays(uint64_t) {} virtual void set_mass(const char *, int, const char *, int, int, int *); virtual void set_mass(const char *, int, int, double); virtual void set_mass(const char *, int, int, char **); From e37de9bb98179562f07b136a428d9370b01da5cf Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 16:05:02 +0000 Subject: [PATCH 12/43] KOKKOS: let DomainKokkos::remap_all and image_flip actually be reached Domain::remap_all() and Domain::image_flip() were not virtual, so the versions in DomainKokkos hid them instead of overriding them. Everything reaches these through a Domain pointer -- fix deform, fix nh, fix bocs, fix npt/cauchy -- so the host implementations ran and read and wrote the plain per-atom arrays while the current copy was on the device, and the KOKKOS versions were never called at all. AddressSanitizer names the read at the faulting line: use-after-poison LAMMPS_NS::Domain::x2lamda(double*, double*) domain.cpp:2433 LAMMPS_NS::Domain::remap_all() domain.cpp:1713 LAMMPS_NS::FixDeform::pre_exchange() fix_deform.cpp:828 Make both virtual. The image flip kernel is a faithful copy of the host loop, but remap_all's is not: it leaves out the velocity correction that fix deform's 'remap v' applies to a wrapped atom. Hand that case back to the host implementation, with the arrays brought over and the write declared, rather than silently drop the correction. examples/VISCOSITY/in.nemd.2d now reproduces the reference build exactly, and the reference build itself is unchanged, which is what says the host path kept its meaning. ModifyKokkos::setup() also ran the computes with no transfer around them, unlike every other loop in that file; compute chunk/atom read the coordinates on the host there. Give them the same treatment as the fixes. (cherry picked from commit 26c03921ada17604ef0b6d247fc8a7545bdf428d) --- src/KOKKOS/domain_kokkos.cpp | 11 +++++++++++ src/KOKKOS/domain_kokkos.h | 4 ++-- src/KOKKOS/modify_kokkos.cpp | 13 ++++++++++++- src/domain.h | 4 ++-- 4 files changed, 27 insertions(+), 5 deletions(-) diff --git a/src/KOKKOS/domain_kokkos.cpp b/src/KOKKOS/domain_kokkos.cpp index c256fca2167..d7efce81243 100644 --- a/src/KOKKOS/domain_kokkos.cpp +++ b/src/KOKKOS/domain_kokkos.cpp @@ -416,6 +416,17 @@ void DomainKokkos::pbc() void DomainKokkos::remap_all() { + // the kernel below does not carry the velocity correction that fix deform's + // "remap v" applies when an atom is wrapped, so leave that case to the host + // implementation rather than drop the correction + + if (deform_vremap) { + atomKK->sync(Host,X_MASK | V_MASK | IMAGE_MASK | MASK_MASK); + Domain::remap_all(); + atomKK->modified(Host,X_MASK | V_MASK | IMAGE_MASK); + return; + } + atomKK->sync(Device,X_MASK | IMAGE_MASK); x = atomKK->k_x.view_device(); diff --git a/src/KOKKOS/domain_kokkos.h b/src/KOKKOS/domain_kokkos.h index ad98acbcab0..16d6a442b0e 100644 --- a/src/KOKKOS/domain_kokkos.h +++ b/src/KOKKOS/domain_kokkos.h @@ -34,8 +34,8 @@ class DomainKokkos : public Domain { ~DomainKokkos() override = default; void reset_box() override; void pbc() override; - void remap_all(); - void image_flip(int, int, int); + void remap_all() override; + void image_flip(int, int, int) override; void x2lamda(int) override; void x2lamda(int,int) override; void lamda2x(int) override; diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp index 26ee88ff513..06cc778ce14 100644 --- a/src/KOKKOS/modify_kokkos.cpp +++ b/src/KOKKOS/modify_kokkos.cpp @@ -53,7 +53,18 @@ void ModifyKokkos::setup(int vflag) } } - for (int i = 0; i < ncompute; i++) compute[i]->setup(); + // the computes get the same treatment as the fixes above: several of them + // read the per-atom arrays in setup(), through the plain pointers when they + // have no KOKKOS version, and were reaching them without a transfer + + for (int i = 0; i < ncompute; i++) { + atomKK->sync(compute[i]->execution_space,compute[i]->datamask_read); + int prev_auto_sync = lmp->kokkos->auto_sync; + if (!compute[i]->kokkosable) lmp->kokkos->auto_sync = 1; + compute[i]->setup(); + lmp->kokkos->auto_sync = prev_auto_sync; + atomKK->modified(compute[i]->execution_space,compute[i]->datamask_modify); + } if (update->whichflag == 1) for (int i = 0; i < nfix; i++) { diff --git a/src/domain.h b/src/domain.h index bd6affa123a..89e1a2c9d49 100644 --- a/src/domain.h +++ b/src/domain.h @@ -142,13 +142,13 @@ class Domain : protected Pointers { void closest_image(const double *const, const double *const, double *const); void remap(double *, imageint &); void remap(double *); - void remap_all(); + virtual void remap_all(); void remap_near(double *, double *); void unmap_inv(double *x, imageint); void unmap(double *, imageint); void unmap(const double *, imageint, double *); void unmap(const double *, const double *, imageint, int, double *, double *); - void image_flip(int, int, int); + virtual void image_flip(int, int, int); int ownatom(int, double *, imageint *, int); void define_general_triclinic(double *, double *, double *, double *); From a87bcc04b36071881fae7eb05eda6d79b0fa080b Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 22 Aug 2026 16:10:16 +0000 Subject: [PATCH 13/43] Declare the set command's writes to the per-atom arrays The set command reads and writes the per-atom arrays through the plain pointers. With the KOKKOS package the current copy can be on the device, so the values it read could be stale and the values it wrote were dropped by the next transfer, with nothing to say so. examples/MC-LOOP/in.mc moves one atom per Monte Carlo step with 'set atom $i x' and takes the energy of the result, and one move in two thousand came out at the energy of the unmoved configuration. Bracket the action loop rather than each of the thirty keywords: the command already walks the arrays it was given, and the two calls do nothing without the package. (cherry picked from commit fcbdc04714bb749f41b46ee8f8dbd9852aa9051d) --- src/set.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/set.cpp b/src/set.cpp index 5fae968fb8d..4ca43ef5cac 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -16,6 +16,7 @@ #include "arg_info.h" #include "atom.h" +#include "atom_masks.h" #include "atom_vec.h" #include "atom_vec_body.h" #include "atom_vec_ellipsoid.h" @@ -534,6 +535,12 @@ void Set::selection(int n) void Set::invoke_actions() { + // every action below reads and writes the per-atom arrays through the plain + // pointers, so bring the host side up to date first and hand the writes over + // afterwards; without the KOKKOS package these do nothing + + atom->sync_host_arrays(ALL_MASK); + // reallocate per-atom variable storage if needed if (varflag && atom->nlocal > maxvariable) { @@ -593,6 +600,8 @@ void Set::invoke_actions() action->count_select = count_select; action->count_action = count_action; } + + atom->modified_host_arrays(ALL_MASK); } /* ---------------------------------------------------------------------- */ From 1ef451554d40eb66295443a566aefd294787e1cd Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sun, 23 Aug 2026 00:26:59 +0000 Subject: [PATCH 14/43] Stop fix qtpie/reaxff from reading the pair style's device neighbor list fix qtpie/reaxff borrows pair reaxff's neighbor list when there is one, to avoid building a second one, and reads it through the plain ilist, numneigh and firstneigh. The KOKKOS version of the pair style builds its list on the device, where those are not filled in, so the fix indexed on uninitialized values: AddressSanitizer: SEGV on unknown address LAMMPS_NS::FixQtpieReaxFF::init_storage() fix_qtpie_reaxff.cpp:665 LAMMPS_NS::FixQtpieReaxFF::setup_pre_force(int) LAMMPS_NS::ModifyKokkos::setup_pre_force(int) examples/reaxff/water/in.water.qtpie and in.water.qtpie.field both crashed this way. This fix has no KOKKOS version -- fix qeq/reaxff does, which is why the same code path is never reached there. Fall back to the list this fix requests for itself, which is the list it already uses when there is no reaxff pair style at all, whenever the pair style is a KOKKOS one. Both inputs now run and reproduce the non-KOKKOS results exactly. (cherry picked from commit 3299c05f55305ee212c3d6003e416d85028162e5) --- src/REAXFF/fix_qtpie_reaxff.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/REAXFF/fix_qtpie_reaxff.cpp b/src/REAXFF/fix_qtpie_reaxff.cpp index caacd4a297e..17b4fc5e2dd 100644 --- a/src/REAXFF/fix_qtpie_reaxff.cpp +++ b/src/REAXFF/fix_qtpie_reaxff.cpp @@ -155,6 +155,15 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) : // register with Atom class reaxff = dynamic_cast(force->pair_match("^reaxff",0)); + + // this fix borrows the pair style's neighbor list when it can, to avoid + // building a second one. The KOKKOS version of pair reaxff builds its list + // on the device, where the plain ilist, numneigh and firstneigh below are + // not filled in, so fall back to the list this fix requests for itself -- + // the same one it uses when there is no reaxff pair style at all. + + if (reaxff && reaxff->kokkosable) reaxff = nullptr; + reaxflag = 0; nlevels_respa = 1; From 2a4e6b88e2ae486546e78069f06d313715ce0067 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sun, 23 Aug 2026 00:58:22 +0000 Subject: [PATCH 15/43] KOKKOS: name fix ilves among the fixes that turn bonds off Neighbor::init_topology() chooses between the all and the partial variant of the topology list builders. The partial one skips a bond whose type has been made negative; the all one copies the type into the list without looking at its sign. A fix that turns bonds off does so after this choice is taken, so it cannot be found by the scan of the types that follows and has to be named instead. The list in the host version names shake, rattle and ilves; the KOKKOS version was never given ilves. fix ilves therefore got the all variant, its negative bond types reached the list, and the bond style indexed its coefficient arrays with them: AddressSanitizer: heap-buffer-overflow BondHarmonicKokkos::operator()( TagBondHarmonicCompute<1, 1>, int const&, s_EV_FLOAT&) const bond_harmonic_kokkos.cpp:170 -- d_r0[type] examples/PACKAGES/ilves/in.rhodo-ilves reported E_bond = 1.6e9 where the run without the package gives 2537.99, every other term agreeing, and later stopped with 'Non-numeric box dimensions'. It now reproduces the non-KOKKOS energies, and AddressSanitizer is quiet. (cherry picked from commit c98911b21094bdd54db2919ede1c8b3c5ebdedce) --- src/KOKKOS/neigh_bond_kokkos.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/KOKKOS/neigh_bond_kokkos.cpp b/src/KOKKOS/neigh_bond_kokkos.cpp index 4d1e5988b6e..f7cf9a6f5f6 100644 --- a/src/KOKKOS/neigh_bond_kokkos.cpp +++ b/src/KOKKOS/neigh_bond_kokkos.cpp @@ -124,6 +124,12 @@ void NeighBondKokkos::init_topology_kk() { int i,m; int bond_off = 0; int angle_off = 0; + // keep this list the same as the one in Neighbor::init_topology(): a fix + // that turns bonds off by making their type negative has to be named here, + // because it does so after this decision is taken and the scan below cannot + // see it yet. Without the name the all variant is chosen, and that one + // copies the type into the list without looking at its sign. + for (const auto &ifix : modify->get_fix_list()) if (utils::strmatch(ifix->style,"^shake") || utils::strmatch(ifix->style,"^rattle") || utils::strmatch(ifix->style,"^ilves")) From cd9ddd2dd88e2fdd5483eccf430f23df9e055008 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 15:25:32 +0000 Subject: [PATCH 16/43] KOKKOS: leave the atom syncing for improper/harmonic to the integrator The style synced and claimed the per-atom arrays itself at the top of compute(). When host executing and device executing styles overlap, the integrator deliberately keeps the force mask out of its own sync and modify calls: the host force buffer is zeroed, the host styles accumulate into it alone, and it is merged into the device buffer afterwards. A style that syncs for itself does not know that, and pulls the force array back from the device in the middle of that accumulation, so the pair force lands in the host buffer as well and the merge adds it twice. In rhodopsin with the bonded styles on the host, that dies a few steps later with atoms out of range in PPPM. The bond, angle and dihedral styles already leave this to the integrator; the impropers were missed. --- src/KOKKOS/improper_harmonic_kokkos.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/KOKKOS/improper_harmonic_kokkos.cpp b/src/KOKKOS/improper_harmonic_kokkos.cpp index 6fc8ebdf530..cd95e2e0d15 100644 --- a/src/KOKKOS/improper_harmonic_kokkos.cpp +++ b/src/KOKKOS/improper_harmonic_kokkos.cpp @@ -88,11 +88,8 @@ void ImproperHarmonicKokkos::compute(int eflag_in, int vflag_in) } else Kokkos::deep_copy(d_vatom,0.0); } - atomKK->sync(execution_space,datamask_read); k_k.template sync(); k_chi.template sync(); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); From 6f9671dd3bf13feab247e41f9eec63bbee2f1fa3 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 15:25:33 +0000 Subject: [PATCH 17/43] KOKKOS: leave the atom syncing for the improper styles to the integrator Same as improper/harmonic: each of these synced and claimed the per-atom arrays itself at the top of compute(), which puts the force array back in play in the middle of the overlap path, where the host force buffer is zeroed, filled by the host styles alone and merged into the device buffer afterwards. The pair force then lands in the host buffer as well and the merge counts it twice. The integrator's own calls are the masked ones, which reduce to the plain masks when nothing is excluded, so this loses nothing on the ordinary path. For the pair style it was issuing both the plain and the masked form, and the plain one re-armed exactly what the masked one leaves out; keep the masked one alone. bond/quartic keeps its calls: it copies to the legacy host arrays in the middle of compute() for the bond breaking that follows, which is a different pattern and is being looked at separately (lammps#5037). With rhodopsin and the bonded styles on the host, the run no longer dies with atoms out of range in PPPM, and its trajectory now matches the same run with every style on the device. --- src/KOKKOS/improper_class2_kokkos.cpp | 3 --- src/KOKKOS/improper_cossq_kokkos.cpp | 3 --- src/KOKKOS/improper_cvff_kokkos.cpp | 3 --- src/KOKKOS/improper_distance_kokkos.cpp | 3 --- src/KOKKOS/improper_distharm_kokkos.cpp | 3 --- src/KOKKOS/improper_fourier_kokkos.cpp | 3 --- src/KOKKOS/improper_inversion_harmonic_kokkos.cpp | 3 --- src/KOKKOS/improper_ring_kokkos.cpp | 3 --- src/KOKKOS/improper_sqdistharm_kokkos.cpp | 3 --- src/KOKKOS/improper_umbrella_kokkos.cpp | 3 --- src/KOKKOS/verlet_kokkos.cpp | 6 ++++-- 11 files changed, 4 insertions(+), 32 deletions(-) diff --git a/src/KOKKOS/improper_class2_kokkos.cpp b/src/KOKKOS/improper_class2_kokkos.cpp index b3139d67c48..346356a813d 100644 --- a/src/KOKKOS/improper_class2_kokkos.cpp +++ b/src/KOKKOS/improper_class2_kokkos.cpp @@ -82,7 +82,6 @@ void ImproperClass2Kokkos::compute(int eflag_in, int vflag_in) d_vatom = k_vatom.template view(); } - atomKK->sync(execution_space,datamask_read); k_k0.template sync(); k_chi0.template sync(); k_aa_k1.template sync(); @@ -95,8 +94,6 @@ void ImproperClass2Kokkos::compute(int eflag_in, int vflag_in) k_setflag_i.template sync(); k_setflag_aa.template sync(); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_cossq_kokkos.cpp b/src/KOKKOS/improper_cossq_kokkos.cpp index 17e0f155d9d..0c0f06dfac4 100644 --- a/src/KOKKOS/improper_cossq_kokkos.cpp +++ b/src/KOKKOS/improper_cossq_kokkos.cpp @@ -88,9 +88,6 @@ void ImproperCossqKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_cvff_kokkos.cpp b/src/KOKKOS/improper_cvff_kokkos.cpp index d76208c3705..de84807a94d 100644 --- a/src/KOKKOS/improper_cvff_kokkos.cpp +++ b/src/KOKKOS/improper_cvff_kokkos.cpp @@ -88,12 +88,9 @@ void ImproperCvffKokkos::compute(int eflag_in, int vflag_in) } else Kokkos::deep_copy(d_vatom,0.0); } - atomKK->sync(execution_space,datamask_read); k_k.template sync(); k_sign.template sync(); k_multiplicity.template sync(); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_distance_kokkos.cpp b/src/KOKKOS/improper_distance_kokkos.cpp index f229ee4504c..e088ee256b0 100644 --- a/src/KOKKOS/improper_distance_kokkos.cpp +++ b/src/KOKKOS/improper_distance_kokkos.cpp @@ -81,9 +81,6 @@ void ImproperDistanceKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_distharm_kokkos.cpp b/src/KOKKOS/improper_distharm_kokkos.cpp index 8f56fbaccee..2c444b6de34 100644 --- a/src/KOKKOS/improper_distharm_kokkos.cpp +++ b/src/KOKKOS/improper_distharm_kokkos.cpp @@ -81,9 +81,6 @@ void ImproperDistHarmKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_fourier_kokkos.cpp b/src/KOKKOS/improper_fourier_kokkos.cpp index bae981e220e..32e96d228f1 100644 --- a/src/KOKKOS/improper_fourier_kokkos.cpp +++ b/src/KOKKOS/improper_fourier_kokkos.cpp @@ -87,9 +87,6 @@ void ImproperFourierKokkos::compute(int eflag_in, int vflag_in) k_C1.template sync(); k_C2.template sync(); k_all.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_inversion_harmonic_kokkos.cpp b/src/KOKKOS/improper_inversion_harmonic_kokkos.cpp index d0ed0674985..4dd1453bdf3 100644 --- a/src/KOKKOS/improper_inversion_harmonic_kokkos.cpp +++ b/src/KOKKOS/improper_inversion_harmonic_kokkos.cpp @@ -81,9 +81,6 @@ void ImproperInversionHarmonicKokkos::compute(int eflag_in, int vfla k_kw.template sync(); k_w0.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_ring_kokkos.cpp b/src/KOKKOS/improper_ring_kokkos.cpp index 59b26b1278b..baa01987752 100644 --- a/src/KOKKOS/improper_ring_kokkos.cpp +++ b/src/KOKKOS/improper_ring_kokkos.cpp @@ -83,9 +83,6 @@ void ImproperRingKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_sqdistharm_kokkos.cpp b/src/KOKKOS/improper_sqdistharm_kokkos.cpp index 022bc578719..f994ef4b500 100644 --- a/src/KOKKOS/improper_sqdistharm_kokkos.cpp +++ b/src/KOKKOS/improper_sqdistharm_kokkos.cpp @@ -81,9 +81,6 @@ void ImproperSQDistHarmKokkos::compute(int eflag_in, int vflag_in) k_k.template sync(); k_chi.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/improper_umbrella_kokkos.cpp b/src/KOKKOS/improper_umbrella_kokkos.cpp index 403b6de1657..86e8663bfe6 100644 --- a/src/KOKKOS/improper_umbrella_kokkos.cpp +++ b/src/KOKKOS/improper_umbrella_kokkos.cpp @@ -85,9 +85,6 @@ void ImproperUmbrellaKokkos::compute(int eflag_in, int vflag_in) k_kw.template sync(); k_w0.template sync(); k_C.template sync(); - atomKK->sync(execution_space,datamask_read); - if (eflag || vflag) atomKK->modified(execution_space,datamask_modify); - else atomKK->modified(execution_space,F_MASK); x = atomKK->k_x.view(); f = atomKK->k_f.view(); diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index fbabd6c2bb0..8b3365ee2bf 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -442,11 +442,13 @@ void VerletKokkos::run(int n) if (pair_compute_flag) { int prev_auto_sync = lmp->kokkos->auto_sync; if (!force->pair->kokkosable) lmp->kokkos->auto_sync = 1; - atomKK->sync(force->pair->execution_space,force->pair->datamask_read); + // the masked form only: the mask is the plain one when nothing is + // excluded, and syncing or claiming the full one first would put the + // force array back in play exactly where the overlap path is trying to + // keep it out atomKK->sync(force->pair->execution_space,~(~force->pair->datamask_read|datamask_exclude)); force->pair->compute(eflag,vflag); lmp->kokkos->auto_sync = prev_auto_sync; - atomKK->modified(force->pair->execution_space,force->pair->datamask_modify); atomKK->modified(force->pair->execution_space,~(~force->pair->datamask_modify|datamask_exclude)); timer->stamp(Timer::PAIR); } From c14e90b906d229104a77ec2d4ec4923742227389 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 15:25:33 +0000 Subject: [PATCH 18/43] KOKKOS: sync the forces into bond/quartic before its host pass bond/quartic runs its kernel on the execution space, then copies the positions, forces and bond topology to the legacy host arrays, breaks bonds and corrects the 1-4 pair interaction there, and copies the forces back. That round trip starts from whatever the force array holds on this side, and compute() never made sure it was current: on any path that reaches it without the integrator's sync, the host pass reads one side and writes the other. Claiming the device side afterwards then collides with the still pending host modification, which is the concurrent modification abort reported in lammps#5037; suppressing that abort alone leaves the wrong forces behind. With the fourmol bond/quartic case, four steps under -sf kk came out with a force error of rms 8e3 against the same run without KOKKOS; with the sync they agree exactly. Also turn off the host/device overlap while this style is in use, as bond/hybrid does for a host executing sub-style: the overlap keeps a separate host force buffer that is zeroed, filled by the host executing styles alone and merged into the device buffer afterwards, and this style's round trip would land in the middle of it. --- src/KOKKOS/bond_quartic_kokkos.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/KOKKOS/bond_quartic_kokkos.cpp b/src/KOKKOS/bond_quartic_kokkos.cpp index 19ee0b26097..eae2407d6cf 100644 --- a/src/KOKKOS/bond_quartic_kokkos.cpp +++ b/src/KOKKOS/bond_quartic_kokkos.cpp @@ -19,6 +19,7 @@ #include "bond_quartic_kokkos.h" #include "atom_kokkos.h" +#include "kokkos.h" #include "atom_masks.h" #include "comm.h" #include "force.h" @@ -44,6 +45,15 @@ BondQuarticKokkos::BondQuarticKokkos(LAMMPS *lmp) : BondQuartic(lmp) execution_space = ExecutionSpaceFromDevice::space; datamask_read = X_MASK | F_MASK | ENERGY_MASK | VIRIAL_MASK; datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; + + // this style breaks bonds and corrects the 1-4 pair interaction on the host, + // in the middle of its own compute(), which means copying the forces to the + // host and back again. The overlap path keeps a separate host force buffer, + // zeroed and filled by the host executing styles alone and merged into the + // device buffer afterwards, and a copy in either direction from here lands + // in the middle of that. Turn the overlap off while this style is in use, + // as bond/hybrid does for a host executing sub-style. + lmp->kokkos->allow_overlap = 0; } /* ---------------------------------------------------------------------- */ @@ -90,6 +100,13 @@ void BondQuarticKokkos::compute(int eflag_in, int vflag_in) k_rc.template sync(); k_u0.template sync(); + // the kernel adds to the forces and the host pass below then copies them + // away and back, so they have to be current on this side before it runs + // whatever the caller did; on any path that reaches compute() without the + // integrator's own sync, claiming the device side after the kernel would + // otherwise collide with a pending host modification + atomKK->sync(execution_space,X_MASK|F_MASK); + x = atomKK->k_x.template view(); f = atomKK->k_f.template view(); neighborKK->k_bondlist.template sync(); From 36aa1bebffc8a2d9ead87c3859102933e72c9cfa Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 15:25:33 +0000 Subject: [PATCH 19/43] KOKKOS: finish the tiled grid receives before copying the buffer The tiled forward and reverse comm post every receive into one host buffer, each into its own stretch of it, and then walk the completions with MPI_Waitany. Inside that loop they copy the whole buffer to the device, which is started while the receives that have not completed are still writing into it: the copy reads those stretches early and puts them on the device, and the later iterations copy them again. The brick path does not do this -- it waits for its one receive before copying. Wait for all of them, copy once, then unpack. Found by reading the two paths against each other while looking at the tiled pppm failures in lammps#5037; on a CPU the copy is synchronous and the window is narrow, so this needs checking on a device before the tiled skips come off. --- src/KOKKOS/grid3d_kokkos.cpp | 35 +++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/src/KOKKOS/grid3d_kokkos.cpp b/src/KOKKOS/grid3d_kokkos.cpp index b0b20020944..e200a1e8b80 100644 --- a/src/KOKKOS/grid3d_kokkos.cpp +++ b/src/KOKKOS/grid3d_kokkos.cpp @@ -711,7 +711,7 @@ void Grid3dKokkos:: forward_comm_kspace_tiled(KSpace *kspace, int which, int nper, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { - int i,m,offset; + int m,offset; KokkosBaseFFT* kspaceKKBase = dynamic_cast(kspace); FFT_SCALAR* buf1; @@ -755,14 +755,19 @@ forward_comm_kspace_tiled(KSpace *kspace, int which, int nper, // unpack all received data - for (i = 0; i < nrecv; i++) { - MPI_Waitany(nrecv,requests,&m,MPI_STATUS_IGNORE); + // every receive writes into its own stretch of the one host buffer, and the + // copy to the device takes the whole buffer, so it cannot be started while + // any of them is still in flight: it would read the stretches that have not + // arrived and copy them over the device side. Collect them all first, copy + // once, then unpack. The brick path already waits before it copies. + MPI_Waitall(nrecv,requests,MPI_STATUSES_IGNORE); - if (!lmp->kokkos->gpu_aware_flag) { - k_buf2.modify_host(); - k_buf2.sync(); - } + if (!lmp->kokkos->gpu_aware_flag) { + k_buf2.modify_host(); + k_buf2.sync(); + } + for (m = 0; m < nrecv; m++) { offset = nper * recv[m].offset; kspaceKKBase->unpack_forward_grid_kokkos(which,k_buf2,offset, recv[m].nunpack,k_recv_unpacklist,m); @@ -852,7 +857,7 @@ void Grid3dKokkos:: reverse_comm_kspace_tiled(KSpace *kspace, int which, int nper, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf1, FFT_DAT::tdual_FFT_SCALAR_1d &k_buf2, MPI_Datatype datatype) { - int i,m,offset; + int m,offset; KokkosBaseFFT* kspaceKKBase = dynamic_cast(kspace); @@ -896,14 +901,16 @@ reverse_comm_kspace_tiled(KSpace *kspace, int which, int nper, } // unpack all received data - for (i = 0; i < nsend; i++) { - MPI_Waitany(nsend,requests,&m,MPI_STATUS_IGNORE); + // as in the forward direction: one buffer, one copy, and only once every + // receive into it has landed + MPI_Waitall(nsend,requests,MPI_STATUSES_IGNORE); - if (!lmp->kokkos->gpu_aware_flag) { - k_buf2.modify_host(); - k_buf2.sync(); - } + if (!lmp->kokkos->gpu_aware_flag) { + k_buf2.modify_host(); + k_buf2.sync(); + } + for (m = 0; m < nsend; m++) { offset = nper * send[m].offset; kspaceKKBase->unpack_reverse_grid_kokkos(which,k_buf2,offset, send[m].npack,k_send_packlist,m); From 12d2991990b1796938b4b23afaecd728aab7040c Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 15:25:34 +0000 Subject: [PATCH 20/43] KOKKOS: read the host mask in the host loops of fix nh remap Remapping with a dilate group walks the atoms on the host, and tested each one against "mask", which is the view for the execution space. On a device that is device memory read from host code, which is the inaccessible memory space error for atom:mask reported against fix npt/kk in lammps#5037. Take the host view instead, and sync it once for the loop rather than syncing the positions again for every atom. On a CPU the device view is host readable, so this reads the right values either way and the forces are unchanged (2.9e-15 relative on an npt run with a dilate group, i.e. summation order); the fault it removes only shows on a real device. --- src/KOKKOS/fix_nh_kokkos.cpp | 32 ++++++++++++++++++++------------ 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/src/KOKKOS/fix_nh_kokkos.cpp b/src/KOKKOS/fix_nh_kokkos.cpp index 1f6fee0014f..4c6fd70d2bc 100644 --- a/src/KOKKOS/fix_nh_kokkos.cpp +++ b/src/KOKKOS/fix_nh_kokkos.cpp @@ -325,13 +325,17 @@ void FixNHKokkos::remap() if (allremap) domainKK->x2lamda(nlocal); else { - for ( int i = 0; i < nlocal; i++) - if (mask[i] & dilate_group_bit) { - auto h_x = atomKK->k_x.view_host(); - atomKK->sync(Host,X_MASK); + // this loop runs on the host, so it needs the host side of both arrays: + // "mask" is the view for the execution space, and reading it here is an + // access to device memory from host code. Sync once for the whole loop + // rather than once per atom as well. + atomKK->sync(Host,X_MASK|MASK_MASK); + auto h_x = atomKK->k_x.view_host(); + auto h_mask = atomKK->k_mask.view_host(); + for (int i = 0; i < nlocal; i++) + if (h_mask[i] & dilate_group_bit) domainKK->x2lamda(&h_x(i,0), &h_x(i,0)); - atomKK->modified(Host,X_MASK); - } + atomKK->modified(Host,X_MASK); } if (rfix.size() > 0) @@ -476,13 +480,17 @@ void FixNHKokkos::remap() if (allremap) domainKK->lamda2x(nlocal); else { - for ( int i = 0; i < nlocal; i++) - if (mask[i] & dilate_group_bit) { - auto h_x = atomKK->k_x.view_host(); - atomKK->sync(Host,X_MASK); + // this loop runs on the host, so it needs the host side of both arrays: + // "mask" is the view for the execution space, and reading it here is an + // access to device memory from host code. Sync once for the whole loop + // rather than once per atom as well. + atomKK->sync(Host,X_MASK|MASK_MASK); + auto h_x = atomKK->k_x.view_host(); + auto h_mask = atomKK->k_mask.view_host(); + for (int i = 0; i < nlocal; i++) + if (h_mask[i] & dilate_group_bit) domainKK->lamda2x(&h_x(i,0), &h_x(i,0)); - atomKK->modified(Host,X_MASK); - } + atomKK->modified(Host,X_MASK); } // for (auto &ifix : rfix) ifix->deform(1); From aab4d7a3076210a9728e23fa1cc524cd6c6a4259 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 15:25:34 +0000 Subject: [PATCH 21/43] KOKKOS: clear the host forces too, not just the device ones force_clear() zeroed the device force view and left the host side alone. Anything that adds its forces there instead -- a style running on the host, or a style without KOKKOS support adding into the plain LAMMPS array -- therefore started each step from the forces it added last step, and when the two sides were brought together those were counted again. The total energy climbs step by step: with the fourmol improper case and pair_style zero, which is what the force-style tests use, four steps came out with a relative force error of 4 against the same run without KOKKOS, and the total energy drifted from 95.617 to 96.108 where the reference holds it exactly. With the host side cleared as well the two agree to the last digit and the energy is flat. A run whose styles are all on the device never reads that buffer, which is why this only shows when something runs on the host: the ewald and pppm divergences in lammps#5037 are of that shape, since kspace_style ewald has no KOKKOS variant at all and runs as a plain style inside the KOKKOS run. --- src/KOKKOS/verlet_kokkos.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index 8b3365ee2bf..c7312f25049 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -588,6 +588,16 @@ void VerletKokkos::force_clear() if (force->newton) nall += atomKK->nghost; Kokkos::parallel_for(nall, Zero(atomKK->k_f.view_device())); + + // clear the host side as well. A style that runs on the host adds its + // forces into that buffer, and a style without KOKKOS support adds into + // the plain LAMMPS array behind it; neither is cleared by the loop above, + // so whatever they added last step is still there and gets counted again + // when the two sides are brought together. With everything on the device + // this costs one clear of an array that is about to be overwritten anyway. + Kokkos::deep_copy(LMPHostType(),atomKK->k_f.view_hostkk(),0.0); + Kokkos::deep_copy(LMPHostType(),atomKK->k_f.view_host(),0.0); + atomKK->modified(Device,F_MASK); if (torqueflag) { From f912fb14bdb9090b1ce7465128e6d3cf64d72242 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 15:25:34 +0000 Subject: [PATCH 22/43] KOKKOS: keep the topology lists alive through the neighbor loops Kokkos passes the loop body a copy of NeighBondKokkos and destroys that copy when the loop finishes, so ~NeighBondKokkos() ran on every bond, angle, dihedral and improper build. It released the topology lists that belong to Neighbor, leaving neighbor->bondlist and its three siblings as null pointers from the first neighbor loop onward. A bond, angle, dihedral or improper style without KOKKOS support reads those pointers directly and crashed in setup. Skip the release when the destructor runs on such a copy, the way the other KOKKOS styles already do. Also move the transfer of the topology lists to the host out of the device branch, so the lists reach a style without KOKKOS support on a host run as well. --- src/KOKKOS/neigh_bond_kokkos.cpp | 30 ++++++++++++++++++++++++++++++ src/KOKKOS/neigh_bond_kokkos.h | 1 + src/KOKKOS/neighbor_kokkos.cpp | 23 ++++++++++++----------- 3 files changed, 43 insertions(+), 11 deletions(-) diff --git a/src/KOKKOS/neigh_bond_kokkos.cpp b/src/KOKKOS/neigh_bond_kokkos.cpp index 4d1e5988b6e..aac473fbbaf 100644 --- a/src/KOKKOS/neigh_bond_kokkos.cpp +++ b/src/KOKKOS/neigh_bond_kokkos.cpp @@ -66,6 +66,8 @@ NeighBondKokkos::NeighBondKokkos(LAMMPS *lmp) : Pointers(lmp) maxangle = 0; maxdihedral = 0; maximproper = 0; + + copymode = 0; } /* ---------------------------------------------------------------------- */ @@ -73,6 +75,12 @@ NeighBondKokkos::NeighBondKokkos(LAMMPS *lmp) : Pointers(lmp) template NeighBondKokkos::~NeighBondKokkos() { + // Kokkos hands the loop bodies below a copy of this class as the functor and + // destroys that copy when the loop is done. The topology lists belong to + // Neighbor, not to the copy, so only the original may release them. + + if (copymode) return; + memoryKK->destroy_kokkos(k_bondlist,neighbor->bondlist); memoryKK->destroy_kokkos(k_anglelist,neighbor->anglelist); memoryKK->destroy_kokkos(k_dihedrallist,neighbor->dihedrallist); @@ -256,7 +264,9 @@ void NeighBondKokkos::bond_all() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -337,7 +347,9 @@ void NeighBondKokkos::bond_partial() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -399,7 +411,9 @@ void NeighBondKokkos::bond_check() atomKK->sync(execution_space, X_MASK); k_bondlist.sync(); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,neighbor->nbondlist),*this,flag); + copymode = 0; int flag_all; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); @@ -444,7 +458,9 @@ void NeighBondKokkos::angle_all() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -531,7 +547,9 @@ void NeighBondKokkos::angle_partial() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -601,7 +619,9 @@ void NeighBondKokkos::angle_check() atomKK->sync(execution_space, X_MASK); k_anglelist.sync(); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,neighbor->nanglelist),*this,flag); + copymode = 0; int flag_all; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); @@ -658,7 +678,9 @@ void NeighBondKokkos::dihedral_all() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -750,7 +772,9 @@ void NeighBondKokkos::dihedral_partial() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -825,7 +849,9 @@ void NeighBondKokkos::dihedral_check(int nlist, typename AT::t_int_2 atomKK->sync(execution_space, X_MASK); k_dihedrallist.sync(); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlist),*this,flag); + copymode = 0; int flag_all; MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world); @@ -899,7 +925,9 @@ void NeighBondKokkos::improper_all() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); @@ -991,7 +1019,9 @@ void NeighBondKokkos::improper_partial() Kokkos::deep_copy(d_scalars,0); + copymode = 1; Kokkos::parallel_reduce(Kokkos::RangePolicy(0,nlocal),*this,nmissing); + copymode = 0; Kokkos::deep_copy(h_scalars,d_scalars); diff --git a/src/KOKKOS/neigh_bond_kokkos.h b/src/KOKKOS/neigh_bond_kokkos.h index 38208b99f64..daea6b92db4 100644 --- a/src/KOKKOS/neigh_bond_kokkos.h +++ b/src/KOKKOS/neigh_bond_kokkos.h @@ -90,6 +90,7 @@ class NeighBondKokkos : protected Pointers { int maxbond,maxangle,maxdihedral,maximproper; // size of bond lists int me,nprocs; + int copymode; // 1 while a Kokkos loop holds a copy private: int map_style; diff --git a/src/KOKKOS/neighbor_kokkos.cpp b/src/KOKKOS/neighbor_kokkos.cpp index 41c5386bb8c..e0a81b63c27 100644 --- a/src/KOKKOS/neighbor_kokkos.cpp +++ b/src/KOKKOS/neighbor_kokkos.cpp @@ -375,17 +375,6 @@ void NeighborKokkos::build_topology() { k_dihedrallist = neighbond_device.k_dihedrallist; k_improperlist = neighbond_device.k_improperlist; - // Transfer topology neighbor lists to Host for non-Kokkos styles - - if (force->bond && force->bond->execution_space == Host) - k_bondlist.sync_host(); - if (force->angle && force->angle->execution_space == Host) - k_anglelist.sync_host(); - if (force->dihedral && force->dihedral->execution_space == Host) - k_dihedrallist.sync_host(); - if (force->improper && force->improper->execution_space == Host) - k_improperlist.sync_host(); - } else { neighbond_host.build_topology_kk(); @@ -394,4 +383,16 @@ void NeighborKokkos::build_topology() { k_dihedrallist = neighbond_host.k_dihedrallist; k_improperlist = neighbond_host.k_improperlist; } + + // transfer topology neighbor lists to the host for non-Kokkos styles, + // which read them through the plain pointers in Neighbor + + if (force->bond && force->bond->execution_space == Host) + k_bondlist.sync_host(); + if (force->angle && force->angle->execution_space == Host) + k_anglelist.sync_host(); + if (force->dihedral && force->dihedral->execution_space == Host) + k_dihedrallist.sync_host(); + if (force->improper && force->improper->execution_space == Host) + k_improperlist.sync_host(); } From 0405a66c23fff19153a8ebd067ac7787a836ea35 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 16:16:06 +0000 Subject: [PATCH 23/43] KOKKOS: drop the rRESPA case from the nve_x comment The same change that moved this claim also refuses run_style respa with the KOKKOS package, so the case the comment named as the one that hits it can no longer happen. Claiming after the kernel is right either way: a claim taken up front is taken by whatever copy runs in between. --- src/KOKKOS/fix_nh_kokkos.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/fix_nh_kokkos.cpp b/src/KOKKOS/fix_nh_kokkos.cpp index 30803272579..96cfd7e14d7 100644 --- a/src/KOKKOS/fix_nh_kokkos.cpp +++ b/src/KOKKOS/fix_nh_kokkos.cpp @@ -635,8 +635,8 @@ void FixNHKokkos::nve_x() copymode = 0; // claim the coordinates after the kernel has written them, not before: a - // claim made up front is taken by any copy that runs in between, and under - // rRESPA one does, which leaves the new coordinates on the device unclaimed + // claim made up front can be taken by a copy that runs in between, which + // leaves the coordinates this kernel writes unclaimed atomKK->modified(execution_space,X_MASK); } From 50245ccbabeaef66347719635f3412d4a1224185 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 17:36:15 +0000 Subject: [PATCH 24/43] KOKKOS: clear the host forces in the branches that were missed The host side of the forces was cleared only in the branch that runs without an include group, and there only for the forces themselves. A run with neigh_modify include takes the other branch, where the loops zero the device copy and leave the host side untouched, so the stale host forces that branch was meant to remove survive after all. The torques and the SPIN forces are overwritten the same way and were never cleared on the host in any branch. Clear the host side next to every device loop that goes with it, for all four arrays. Each clear takes the same range as its device loop rather than the whole view, so the atoms outside the include group are left alone exactly as Verlet::force_clear() leaves them; that also makes the range in the first branch agree with the loop above it, which reaches nall and not the end of the array. Both host views have to be cleared: the Kokkos host view that a style running on the host adds into, and the plain LAMMPS array behind it that a style without KOKKOS support adds into. On a host-only Kokkos build the two sides share an allocation, so this is redundant there and the results are unchanged; the runs it matters for are the ones with a real device. --- src/KOKKOS/verlet_kokkos.cpp | 48 ++++++++++++++++++++++++++++-------- 1 file changed, 38 insertions(+), 10 deletions(-) diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index f8adf3c604c..736f866d664 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/src/KOKKOS/verlet_kokkos.cpp @@ -62,6 +62,32 @@ struct Zero { } }; +/* ---------------------------------------------------------------------- + zero count entries of a per-atom array on both host sides, from first + + the loops that zero the device copy leave the host side alone. A style + that runs on the host adds its forces into the Kokkos host view, and a + style without KOKKOS support adds into the plain LAMMPS array behind it, + so whatever either of them added last step is still there and is counted + again when the two sides are brought together. Each call takes the same + range as the device loop it goes with, so that the atoms the plain code + leaves alone -- those outside an include group -- are left alone here too. +------------------------------------------------------------------------- */ + +template +static void zero_host_view(const View &v, int first, int count) +{ + Kokkos::parallel_for(Kokkos::RangePolicy(first,first+count),Zero(v)); +} + +template +static void zero_host(const DualView &k, int first, int count) +{ + if (count <= 0) return; + zero_host_view(k.view_hostkk(),first,count); + zero_host_view(k.view_host(),first,count); +} + /* ---------------------------------------------------------------------- */ VerletKokkos::VerletKokkos(LAMMPS *lmp, int narg, char **arg) : @@ -598,20 +624,12 @@ void VerletKokkos::force_clear() if (force->newton) nall += atomKK->nghost; Kokkos::parallel_for(nall, Zero(atomKK->k_f.view_device())); - - // clear the host side as well. A style that runs on the host adds its - // forces into that buffer, and a style without KOKKOS support adds into - // the plain LAMMPS array behind it; neither is cleared by the loop above, - // so whatever they added last step is still there and gets counted again - // when the two sides are brought together. With everything on the device - // this costs one clear of an array that is about to be overwritten anyway. - Kokkos::deep_copy(LMPHostType(),atomKK->k_f.view_hostkk(),0.0); - Kokkos::deep_copy(LMPHostType(),atomKK->k_f.view_host(),0.0); - + zero_host(atomKK->k_f,0,nall); atomKK->modified(Device,F_MASK); if (torqueflag) { Kokkos::parallel_for(nall, Zero(atomKK->k_torque.view_device())); + zero_host(atomKK->k_torque,0,nall); atomKK->modified(Device,TORQUE_MASK); } @@ -619,8 +637,10 @@ void VerletKokkos::force_clear() if (extraflag) { Kokkos::parallel_for(nall, Zero(atomKK->k_fm.view_device())); + zero_host(atomKK->k_fm,0,nall); atomKK->modified(Device,FM_MASK); Kokkos::parallel_for(nall, Zero(atomKK->k_fm_long.view_device())); + zero_host(atomKK->k_fm_long,0,nall); atomKK->modified(Device,FML_MASK); } @@ -630,10 +650,12 @@ void VerletKokkos::force_clear() } else { Kokkos::parallel_for(atomKK->nfirst, Zero(atomKK->k_f.view_device())); + zero_host(atomKK->k_f,0,atomKK->nfirst); atomKK->modified(Device,F_MASK); if (torqueflag) { Kokkos::parallel_for(atomKK->nfirst, Zero(atomKK->k_torque.view_device())); + zero_host(atomKK->k_torque,0,atomKK->nfirst); atomKK->modified(Device,TORQUE_MASK); } @@ -641,18 +663,22 @@ void VerletKokkos::force_clear() if (extraflag) { Kokkos::parallel_for(atomKK->nfirst, Zero(atomKK->k_fm.view_device())); + zero_host(atomKK->k_fm,0,atomKK->nfirst); atomKK->modified(Device,FM_MASK); Kokkos::parallel_for(atomKK->nfirst, Zero(atomKK->k_fm_long.view_device())); + zero_host(atomKK->k_fm_long,0,atomKK->nfirst); atomKK->modified(Device,FML_MASK); } if (force->newton) { auto range = Kokkos::RangePolicy(atomKK->nlocal, atomKK->nlocal + atomKK->nghost); Kokkos::parallel_for(range, Zero(atomKK->k_f.view_device())); + zero_host(atomKK->k_f,atomKK->nlocal,atomKK->nghost); atomKK->modified(Device,F_MASK); if (torqueflag) { Kokkos::parallel_for(range, Zero(atomKK->k_torque.view_device())); + zero_host(atomKK->k_torque,atomKK->nlocal,atomKK->nghost); atomKK->modified(Device,TORQUE_MASK); } @@ -660,8 +686,10 @@ void VerletKokkos::force_clear() if (extraflag) { Kokkos::parallel_for(range, Zero(atomKK->k_fm.view_device())); + zero_host(atomKK->k_fm,atomKK->nlocal,atomKK->nghost); atomKK->modified(Device,FM_MASK); Kokkos::parallel_for(range, Zero(atomKK->k_fm_long.view_device())); + zero_host(atomKK->k_fm_long,atomKK->nlocal,atomKK->nghost); atomKK->modified(Device,FML_MASK); } } From e197052edd81a8eb32000fb1ca4ca39017d88d43 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 18:00:35 +0000 Subject: [PATCH 25/43] KOKKOS: honor the include group when binning atoms NBinKokkos::bin_atoms() put every owned atom and every ghost into the bins. NBinStandard::bin_atoms() bins only the atoms an include group's pairs are built from: the owned atoms of the group, which sorting has put first, and the ghosts that are in the group. The atoms outside the group therefore turned up in the neighbor lists of the atoms inside it, and a run with neigh_modify include computed pairs that the same run without KOKKOS does not. On the peptide example with atom_modify first and neigh_modify include the difference is there from the first step, with E_pair at -31960.18 against -31864.22 without KOKKOS, and it stays about that size for the whole run rather than growing, which is the shape of a different set of pairs rather than of forces accumulating. With the group honored the two runs agree to the last digit. atom2bin is still set for every atom. Nothing reads it for an atom outside the group, and leaving it alone would leave a stale bin behind instead. --- src/KOKKOS/nbin_kokkos.cpp | 25 +++++++++++++++++++++++++ src/KOKKOS/nbin_kokkos.h | 5 +++++ 2 files changed, 30 insertions(+) diff --git a/src/KOKKOS/nbin_kokkos.cpp b/src/KOKKOS/nbin_kokkos.cpp index 45b25bf2ce8..ec78a4c55bc 100644 --- a/src/KOKKOS/nbin_kokkos.cpp +++ b/src/KOKKOS/nbin_kokkos.cpp @@ -17,6 +17,7 @@ #include "atom_kokkos.h" #include "atom_masks.h" #include "comm.h" +#include "group.h" #include "kokkos.h" #include "memory_kokkos.h" #include "update.h" @@ -85,6 +86,12 @@ void NBinKokkos::bin_atoms() { last_bin = update->ntimestep; + // an include group restricts which atoms go into the bins, see below + + includegroup_bitmask = includegroup ? group->bitmask[includegroup] : 0; + includegroup_nfirst = atom->nfirst; + includegroup_nlocal = atom->nlocal; + k_bins.template sync(); k_bincount.template sync(); k_atom2bin.template sync(); @@ -102,6 +109,11 @@ void NBinKokkos::bin_atoms() atomKK->sync(ExecutionSpaceFromDevice::space,X_MASK); x = atomKK->k_x.view(); + if (includegroup_bitmask) { + atomKK->sync(ExecutionSpaceFromDevice::space,MASK_MASK); + mask = atomKK->k_mask.view(); + } + bboxlo_[0] = bboxlo[0]; bboxlo_[1] = bboxlo[1]; bboxlo_[2] = bboxlo[2]; bboxhi_[0] = bboxhi[0]; bboxhi_[1] = bboxhi[1]; bboxhi_[2] = bboxhi[2]; @@ -134,6 +146,19 @@ void NBinKokkos::binatomsItem(const int &i) const const int ibin = coord2bin(x(i, 0), x(i, 1), x(i, 2)); atom2bin(i) = ibin; + + // with an include group only the atoms the group's pairs are built from + // belong in the bins: the owned atoms of the group, which sorting has put + // first, and the ghosts that are in the group. Binning the rest would put + // them in the neighbor lists of the group's atoms, which is what the plain + // NBinStandard::bin_atoms() leaves out. + + if (includegroup_bitmask) { + if (i < includegroup_nlocal) { + if (i >= includegroup_nfirst) return; + } else if (!(mask(i) & includegroup_bitmask)) return; + } + const int ac = Kokkos::atomic_fetch_add(&bincount[ibin], (int)1); if (ac < (int)bins.extent(1)) { bins(ibin, ac) = i; diff --git a/src/KOKKOS/nbin_kokkos.h b/src/KOKKOS/nbin_kokkos.h index e99ba90bf77..6c63ef68499 100644 --- a/src/KOKKOS/nbin_kokkos.h +++ b/src/KOKKOS/nbin_kokkos.h @@ -55,6 +55,11 @@ class NBinKokkos : public NBinStandard { typename AT::t_int_scalar d_resize; HAT::t_int_scalar h_resize; typename AT::t_kkfloat_1d_3_lr_randomread x; + typename AT::t_int_1d mask; + + int includegroup_bitmask; // bit of the include group, 0 if there is none + int includegroup_nfirst; // # of owned atoms in the include group + int includegroup_nlocal; // # of owned atoms // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION From 4e08c04d852efedc89f68ce965ce7e22bdef0430 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 24 Aug 2026 18:02:46 +0000 Subject: [PATCH 26/43] KOKKOS: honor the include group when binning atoms for SSA The same omission as in NBinKokkos, in the twin that the SSA neighbor build uses: every owned atom and every ghost went into the bins, where NBinSSA::bin_atoms() takes only the owned atoms of the include group, which are the first nfirst, and the ghosts that are in the group. The ghosts still start at the number of owned atoms rather than at nfirst. Without an include group the group bit is zero and the ghost range is the number of owned atoms, so nothing changes there; the dpde/shardlow and dpdrx/shardlow examples come out the same as without KOKKOS either way. --- src/KOKKOS/nbin_ssa_kokkos.cpp | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/KOKKOS/nbin_ssa_kokkos.cpp b/src/KOKKOS/nbin_ssa_kokkos.cpp index 596cd55077d..5cceb628548 100644 --- a/src/KOKKOS/nbin_ssa_kokkos.cpp +++ b/src/KOKKOS/nbin_ssa_kokkos.cpp @@ -21,6 +21,7 @@ #include "atom_kokkos.h" #include "domain.h" +#include "group.h" #include "update.h" #include "atom_masks.h" #include "kokkos_type.h" @@ -105,6 +106,22 @@ void NBinSSAKokkos::bin_atoms() int nghost = atom->nghost; int nall = nlocal + nghost; + // an include group leaves out the atoms its pairs are not built from: the + // owned atoms of the group are the first nfirst, and of the ghosts only + // those in the group are binned. This is what NBinSSA::bin_atoms() does, + // and without it the atoms outside the group turn up in the lists of the + // atoms inside it. The ghosts still start at the number of owned atoms. + + const int nowned = nlocal; + const int group_bitmask = includegroup ? group->bitmask[includegroup] : 0; + if (includegroup) nlocal = atom->nfirst; + + typename AT::t_int_1d mask_; + if (includegroup) { + atomKK->sync(ExecutionSpaceFromDevice::space,MASK_MASK); + mask_ = atomKK->k_mask.view(); + } + atomKK->sync(ExecutionSpaceFromDevice::space,X_MASK); x = atomKK->k_x.view(); @@ -141,7 +158,7 @@ void NBinSSAKokkos::bin_atoms() k_gbincount.sync(); ghosts_per_gbin = 0; NPairSSAKokkosBinIDGhostsFunctor f(*this); - Kokkos::parallel_reduce(Kokkos::RangePolicy(nlocal,nall), f, ghosts_per_gbin); + Kokkos::parallel_reduce(Kokkos::RangePolicy(nowned,nall), f, ghosts_per_gbin); } // actually bin the ghost atoms @@ -158,8 +175,9 @@ void NBinSSAKokkos::bin_atoms() auto gbincount_ = gbincount; auto gbins_ = gbins; - Kokkos::parallel_for(Kokkos::RangePolicy(nlocal,nall), + Kokkos::parallel_for(Kokkos::RangePolicy(nowned,nall), LAMMPS_LAMBDA (const int i) { + if (group_bitmask && !(mask_(i) & group_bitmask)) return; const int iAIR = binID_(i); if (iAIR > 0) { // include only ghost atoms in an AIR const int ac = Kokkos::atomic_fetch_add(&gbincount_[iAIR], (int)1); From 5405f6f35324c7c4f0308cf8394c7a1f82daed4a Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 16:46:31 +0000 Subject: [PATCH 27/43] KOKKOS: sync and claim xshake around the device forward comm in fix shake The host pack/unpack pair for the SHAKE unconstrained coordinates brings xshake to the host before reading it and claims the host side after writing it. The device pair did neither: the pack read the device copy without bringing it up to date, and the unpack wrote the ghost coordinates there without claiming them. Today the sequence in post_force() hides this, because unconstrained_update() claims the device side immediately before the communication, so the device copy happens to be the current one and the sync that follows the communication has nothing to do. Nothing in either routine relies on that, though, and if the host side is the current one on entry the pack sends whatever the device copy last held and the sync after the communication then discards the ghosts the unpack just wrote. Results are unchanged on a host-only build, where the two sides share an allocation. --- src/KOKKOS/fix_shake_kokkos.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index 4c2a4e19866..dd8bc374cd9 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -1992,6 +1992,13 @@ int FixShakeKokkos::pack_forward_comm_kokkos(int n, DAT::tdual_int_1 DAT::tdual_double_1d &k_buf, int pbc_flag, int* pbc) { + // the host variant below syncs before it reads and claims after it writes; + // this one has to do the same, or the buffer is packed from whichever side + // the device copy last held and the ghosts it unpacks are discarded by the + // next sync + + k_xshake.sync(); + d_sendlist = k_sendlist.view(); d_buf = k_buf.view(); @@ -2053,6 +2060,8 @@ void FixShakeKokkos::unpack_forward_comm_kokkos(int n, int first_in, first = first_in; d_buf = buf.view(); Kokkos::parallel_for(Kokkos::RangePolicy(0,n),*this); + + k_xshake.modify(); } template From 9ba00b501fc0faf83b6a6c9bddbdd6cf4bfef973 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 19:18:41 +0000 Subject: [PATCH 28/43] KOKKOS: sync atom data in ModifyKokkos::energy_couple/energy_global ModifyKokkos::energy_couple() and ModifyKokkos::energy_global() marked the atom data listed in a fix's datamask_modify as modified after calling compute_scalar(), but they were the only ModifyKokkos wrappers that did not sync that fix's datamask_read first. When the data was still dirty on the host, marking it modified on the device triggered: Kokkos::DualView::modify_device ERROR: Concurrent modification of host and device views in DualView "atom:f" This is reached with "fix langevin" plus "econserve" in thermo_style and a "minimize" command: MinKokkos::setup() turns auto_sync off, then ModifyKokkos::setup() marks all atom data as modified on the host for the internal MINIMIZE/kk fix (it uses the default Host execution space and ALL_MASK datamasks), and the following thermo output calls energy_couple() which marks atom->f as modified on the device without syncing. Add the missing sync to both routines, matching all other ModifyKokkos wrappers. --- src/KOKKOS/modify_kokkos.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp index 06cc778ce14..4ff0fc1b48a 100644 --- a/src/KOKKOS/modify_kokkos.cpp +++ b/src/KOKKOS/modify_kokkos.cpp @@ -462,6 +462,8 @@ double ModifyKokkos::energy_couple() { double energy = 0.0; for (int i = 0; i < n_energy_couple; i++) { + atomKK->sync(fix[list_energy_couple[i]]->execution_space, + fix[list_energy_couple[i]]->datamask_read); int prev_auto_sync = lmp->kokkos->auto_sync; if (!fix[list_energy_couple[i]]->kokkosable) lmp->kokkos->auto_sync = 1; energy += fix[list_energy_couple[i]]->compute_scalar(); @@ -482,6 +484,8 @@ double ModifyKokkos::energy_global() { double energy = 0.0; for (int i = 0; i < n_energy_global; i++) { + atomKK->sync(fix[list_energy_global[i]]->execution_space, + fix[list_energy_global[i]]->datamask_read); int prev_auto_sync = lmp->kokkos->auto_sync; if (!fix[list_energy_global[i]]->kokkosable) lmp->kokkos->auto_sync = 1; energy += fix[list_energy_global[i]]->compute_scalar(); From 46ae8890a82be89489a1ad77f28a7511b3258f42 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 19:18:41 +0000 Subject: [PATCH 29/43] KOKKOS: set EMPTY_MASK datamasks for the internal MINIMIZE/kk fix FixMinimizeKokkos left datamask_read and datamask_modify at the Fix defaults of ALL_MASK, with the default Host execution space. The fix has setmask() == 0 and overrides none of the callbacks Modify invokes, and it syncs the atom data it does touch itself (reset_coords() syncs X_MASK; grow_arrays(), copy_arrays() and the exchange routines only touch its own vectors). So the masks describe atom data the fix never accesses through Modify. The consequence is that ModifyKokkos::setup(), which loops over all fixes, syncs all atom data to the host and then marks all of it as modified on the host for this fix. MINIMIZE/kk is added by Min::init() and is therefore always the last fix in the list, so nothing in that loop syncs the data back, and MinKokkos::setup() has turned auto_sync off. Every atom array is then left dirty on the host for the rest of minimization setup, which is how the concurrent host/device modification of atom->f reported in issue #5080 came about. Declare both masks EMPTY_MASK, as all other Kokkos fixes that manage their own syncing do. --- src/KOKKOS/fix_minimize_kokkos.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/KOKKOS/fix_minimize_kokkos.cpp b/src/KOKKOS/fix_minimize_kokkos.cpp index b1d423c89d7..741e46e7fc9 100644 --- a/src/KOKKOS/fix_minimize_kokkos.cpp +++ b/src/KOKKOS/fix_minimize_kokkos.cpp @@ -31,6 +31,12 @@ FixMinimizeKokkos::FixMinimizeKokkos(LAMMPS *lmp, int narg, char **arg) : { kokkosable = 1; atomKK = (AtomKokkos *) atom; + + // this fix only stores per-atom data of its own; it syncs the atom data it + // touches itself, so Modify must not sync or invalidate any of it + + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; } /* ---------------------------------------------------------------------- */ From a979cd5e1cf9a790f9d9a58a3e3d29ca7b3694ce Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 20:18:15 +0000 Subject: [PATCH 30/43] KOKKOS: initialize cached per-atom pointers in atom_vec_*_kokkos constructors The atom_vec_*_kokkos classes cache the per-atom arrays of the corresponding atom->... fields and refresh them in grow_pointers(). Their constructors left those members uninitialized, so the pointers are indeterminate if the object is destroyed or re-grown before grow_pointers() first runs. Coverity flags each constructor with UNINIT (uninitialized pointer field). Initialize the members to nullptr in the constructor initializer lists, in declaration order. AtomVecHybridKokkos already did this for stylesKK. Fixes #5055 --- src/KOKKOS/atom_vec_angle_kokkos.cpp | 5 +++-- src/KOKKOS/atom_vec_bond_kokkos.cpp | 5 +++-- src/KOKKOS/atom_vec_dpd_kokkos.cpp | 4 ++-- src/KOKKOS/atom_vec_full_kokkos.cpp | 8 ++++++-- src/KOKKOS/atom_vec_molecular_kokkos.cpp | 8 ++++++-- src/KOKKOS/atom_vec_sphere_kokkos.cpp | 4 ++-- 6 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/KOKKOS/atom_vec_angle_kokkos.cpp b/src/KOKKOS/atom_vec_angle_kokkos.cpp index 2b6ed3a3d1e..c0e9045d123 100644 --- a/src/KOKKOS/atom_vec_angle_kokkos.cpp +++ b/src/KOKKOS/atom_vec_angle_kokkos.cpp @@ -27,8 +27,9 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecAngleKokkos::AtomVecAngleKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecAngle(lmp) +AtomVecAngleKokkos::AtomVecAngleKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecAngle(lmp), molecule(nullptr), special(nullptr), + bond_atom(nullptr), angle_atom1(nullptr), angle_atom2(nullptr), angle_atom3(nullptr) { } diff --git a/src/KOKKOS/atom_vec_bond_kokkos.cpp b/src/KOKKOS/atom_vec_bond_kokkos.cpp index 5e8a924c9e9..cdf6e46a1f9 100644 --- a/src/KOKKOS/atom_vec_bond_kokkos.cpp +++ b/src/KOKKOS/atom_vec_bond_kokkos.cpp @@ -26,8 +26,9 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecBondKokkos::AtomVecBondKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecBond(lmp) +AtomVecBondKokkos::AtomVecBondKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecBond(lmp), molecule(nullptr), special(nullptr), + bond_atom(nullptr) { } diff --git a/src/KOKKOS/atom_vec_dpd_kokkos.cpp b/src/KOKKOS/atom_vec_dpd_kokkos.cpp index be7a5d4f697..6b6995356b7 100644 --- a/src/KOKKOS/atom_vec_dpd_kokkos.cpp +++ b/src/KOKKOS/atom_vec_dpd_kokkos.cpp @@ -28,8 +28,8 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecDPDKokkos::AtomVecDPDKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecDPD(lmp) +AtomVecDPDKokkos::AtomVecDPDKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecDPD(lmp), duChem(nullptr) { } diff --git a/src/KOKKOS/atom_vec_full_kokkos.cpp b/src/KOKKOS/atom_vec_full_kokkos.cpp index e1e5de88f0d..6caaba7d007 100644 --- a/src/KOKKOS/atom_vec_full_kokkos.cpp +++ b/src/KOKKOS/atom_vec_full_kokkos.cpp @@ -26,8 +26,12 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecFullKokkos::AtomVecFullKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecFull(lmp) +AtomVecFullKokkos::AtomVecFullKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecFull(lmp), q(nullptr), molecule(nullptr), + special(nullptr), bond_atom(nullptr), angle_atom1(nullptr), angle_atom2(nullptr), + angle_atom3(nullptr), dihedral_atom1(nullptr), dihedral_atom2(nullptr), + dihedral_atom3(nullptr), dihedral_atom4(nullptr), improper_atom1(nullptr), + improper_atom2(nullptr), improper_atom3(nullptr), improper_atom4(nullptr) { } diff --git a/src/KOKKOS/atom_vec_molecular_kokkos.cpp b/src/KOKKOS/atom_vec_molecular_kokkos.cpp index ae72ea7b164..5c74cb917f6 100644 --- a/src/KOKKOS/atom_vec_molecular_kokkos.cpp +++ b/src/KOKKOS/atom_vec_molecular_kokkos.cpp @@ -27,8 +27,12 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -AtomVecMolecularKokkos::AtomVecMolecularKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecMolecular(lmp) +AtomVecMolecularKokkos::AtomVecMolecularKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecMolecular(lmp), molecule(nullptr), special(nullptr), + bond_atom(nullptr), angle_atom1(nullptr), angle_atom2(nullptr), angle_atom3(nullptr), + dihedral_atom1(nullptr), dihedral_atom2(nullptr), dihedral_atom3(nullptr), + dihedral_atom4(nullptr), improper_atom1(nullptr), improper_atom2(nullptr), + improper_atom3(nullptr), improper_atom4(nullptr) { } diff --git a/src/KOKKOS/atom_vec_sphere_kokkos.cpp b/src/KOKKOS/atom_vec_sphere_kokkos.cpp index 83a8ff37ff4..a40a246daec 100644 --- a/src/KOKKOS/atom_vec_sphere_kokkos.cpp +++ b/src/KOKKOS/atom_vec_sphere_kokkos.cpp @@ -30,8 +30,8 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -AtomVecSphereKokkos::AtomVecSphereKokkos(LAMMPS *lmp) : AtomVec(lmp), -AtomVecKokkos(lmp), AtomVecSphere(lmp) +AtomVecSphereKokkos::AtomVecSphereKokkos(LAMMPS *lmp) : + AtomVec(lmp), AtomVecKokkos(lmp), AtomVecSphere(lmp), torque(nullptr) { } From 811895e7ad3ed572866e38628b3f066df58bceee Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 20:18:24 +0000 Subject: [PATCH 31/43] KOKKOS: fix per-atom energy/virial memory leak in tip4p/kk pair styles PairTIP4PKokkos::prepare() called ev_init() with the default alloc = 1, so Pair::ev_setup() allocated the plain base-class eatom and vatom arrays. prepare() then replaced those pointers with Kokkos dual views via destroy_kokkos()/create_kokkos(); since destroy_kokkos() only clears the pointer, the plain allocations were orphaned and never freed. On a later reallocation Pair::ev_setup() would in addition call memory->destroy() on a Kokkos-owned pointer. Pass alloc = 0, as all other KOKKOS pair styles do, so the per-atom arrays are managed only through k_eatom/k_vatom. This affects tip4p/cut/kk, tip4p/long/kk, lj/cut/tip4p/cut/kk and lj/cut/tip4p/long/kk. Verified with valgrind on all four styles with per-atom energy and virial output enabled; forces, energies and stresses are unchanged. The corresponding leak in pace/kk was already fixed separately. Fixes #5062 --- src/KOKKOS/pair_tip4p_kokkos.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/KOKKOS/pair_tip4p_kokkos.h b/src/KOKKOS/pair_tip4p_kokkos.h index dac73aa7f68..bca63c65e25 100644 --- a/src/KOKKOS/pair_tip4p_kokkos.h +++ b/src/KOKKOS/pair_tip4p_kokkos.h @@ -389,7 +389,9 @@ class PairTIP4PKokkos : public PairCPUBase { { this->eflag = eflag_in; this->vflag = vflag_in; - this->ev_init(this->eflag,this->vflag); + // alloc = 0: the per-atom energy/virial arrays are allocated below through + // the Kokkos dual views, so Pair::ev_setup() must not allocate plain arrays + this->ev_init(this->eflag,this->vflag,0); this->atomKK->sync(this->execution_space,this->datamask_read); From 00d4c8c59a3bb9fac5f47b1b3672c842c79a7fef Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 20:52:24 +0000 Subject: [PATCH 32/43] KOKKOS: set EMPTY_MASK datamasks for fix property/atom FixPropertyAtomKokkos left datamask_read and datamask_modify at the Fix defaults of ALL_MASK, with the default Host execution space, the same way the internal MINIMIZE/kk fix did. The fix has setmask() == 0 and overrides none of the callbacks Modify invokes, and the per-atom arrays it owns are synced through AtomKokkos::sync() and AtomKokkos::modified(), which call into its own sync() and modified() rather than consulting the datamasks. So the masks describe atom data the fix never accesses through Modify. ModifyKokkos::setup() loops over all fixes and therefore marks every atom array as modified on the host for this fix. When "fix property/atom" is the last fix defined, nothing later in that loop syncs the data back, and the atom data stays dirty on the host through thermo output. Reading econserve then marks atom->f as modified on the device and aborts with a concurrent host/device modification, the same failure as issue #5080 but reachable from a plain run, with no minimize involved. Declare both masks EMPTY_MASK, as all other Kokkos fixes that manage their own syncing do. --- src/KOKKOS/fix_property_atom_kokkos.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/KOKKOS/fix_property_atom_kokkos.cpp b/src/KOKKOS/fix_property_atom_kokkos.cpp index ecb60495b14..fc003b71711 100644 --- a/src/KOKKOS/fix_property_atom_kokkos.cpp +++ b/src/KOKKOS/fix_property_atom_kokkos.cpp @@ -31,6 +31,12 @@ FixPropertyAtomKokkos::FixPropertyAtomKokkos(LAMMPS *lmp, int narg, char **arg) atomKK = (AtomKokkos *) atom; kokkosable = 1; + // this fix syncs the atom data it owns itself, through AtomKokkos::sync() + // and AtomKokkos::modified(); Modify must not sync or invalidate any of it + + datamask_read = EMPTY_MASK; + datamask_modify = EMPTY_MASK; + dvector_flag = 0; ivector_flag = 0; iarray_flag = 0; From 95cf1400b0d3da679f2c0658db03da90c4f1d6b5 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 21:28:42 +0000 Subject: [PATCH 33/43] KOKKOS: fix per-atom virial memory leak and out-of-bounds access in fixes The KOKKOS fixes that keep the per-atom virial in a Kokkos dual view first let the base class allocate the plain vatom array and then replaced the pointer with destroy_kokkos()/create_kokkos(). Since destroy_kokkos() only clears the pointer, the plain allocation was orphaned and never freed. This is the same defect as in the tip4p/kk and pace/kk pair styles, where Pair::ev_setup() already takes an alloc argument to suppress the allocation. Give Fix::ev_setup() and Fix::v_setup() (and the ev_init()/v_init() wrappers) the same alloc argument, defaulting to 1 so that all other fixes are unaffected, and pass alloc = 0 from fix shake/kk, fix efield/kk and fix wall/region/kk. For the fix wall styles the per-atom virial was additionally allocated before FixWall::post_force() called v_init(), i.e. before vflag_atom was known. On the first invocation with per-atom virial output the dual view was therefore still empty while the kernels tallied into it, which segfaulted: fix wf all wall/lj93 zlo -9.0 1.0 1.0 6.0 units box compute st all stress/atom NULL fix Add a v_setup_peratom() hook to FixWall that the accelerated styles override, so the dual view is (re)allocated after v_init() has set the flags. fix wall/harmonic/kk, wall/lj93/kk, wall/lj126/kk, wall/lj1043/kk and wall/morse/kk now reproduce the per-atom virial of their CPU counterparts exactly. Also drop a duplicated v_init() call in fix wall/region/kk and correct two copy-and-paste memory labels in fix shake/kk. Verified with valgrind on all eight styles: no bytes definitely lost, no invalid reads or writes. Per-atom and global virial, energy and thermo output are identical to the plain CPU styles, and the KOKKOS results of examples/{melt,micelle,indent,min,crack} are unchanged. --- src/KOKKOS/fix_efield_kokkos.cpp | 7 ++-- src/KOKKOS/fix_shake_kokkos.cpp | 26 ++++++++++++--- src/KOKKOS/fix_wall_harmonic_kokkos.cpp | 15 ++++++++- src/KOKKOS/fix_wall_harmonic_kokkos.h | 1 + src/KOKKOS/fix_wall_lj1043_kokkos.cpp | 15 ++++++++- src/KOKKOS/fix_wall_lj1043_kokkos.h | 1 + src/KOKKOS/fix_wall_lj126_kokkos.cpp | 15 ++++++++- src/KOKKOS/fix_wall_lj126_kokkos.h | 1 + src/KOKKOS/fix_wall_lj93_kokkos.cpp | 14 ++++++-- src/KOKKOS/fix_wall_lj93_kokkos.h | 1 + src/KOKKOS/fix_wall_morse_kokkos.cpp | 15 ++++++++- src/KOKKOS/fix_wall_morse_kokkos.h | 1 + src/KOKKOS/fix_wall_region_kokkos.cpp | 11 +++---- src/fix.cpp | 44 +++++++++++++++---------- src/fix.h | 18 ++++++---- src/fix_wall.cpp | 2 +- src/fix_wall.h | 8 +++++ 17 files changed, 153 insertions(+), 42 deletions(-) diff --git a/src/KOKKOS/fix_efield_kokkos.cpp b/src/KOKKOS/fix_efield_kokkos.cpp index 609045cb3f9..1fcf88a4769 100644 --- a/src/KOKKOS/fix_efield_kokkos.cpp +++ b/src/KOKKOS/fix_efield_kokkos.cpp @@ -93,9 +93,12 @@ void FixEfieldKokkos::post_force(int vflag) // virial setup - v_init(vflag); + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) - // reallocate per-atom arrays if necessary + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index 4fff360a394..0c38c2be666 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -345,7 +345,22 @@ template void FixShakeKokkos::min_post_force(int vflag) { int eflag = eflag_pre_reverse; - ev_init(eflag, vflag); + + // the per-atom virial is accumulated into a dual view, so the plain + // base-class arrays must not be allocated here (alloc = 0) + + ev_init(eflag, vflag, 0); + + // reallocate the per-atom virial dual view if necessary. minimization + // only tallies the global virial, but the freshly created view keeps the + // per-atom virial zeroed for computes reading it + + if (vflag_atom) { + memoryKK->destroy_kokkos(k_vatom,vatom); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"shake:vatom"); + d_vatom = k_vatom.template view(); + } + ebond = 0.0; atomKK->sync(execution_space, X_MASK | F_MASK); @@ -568,13 +583,16 @@ void FixShakeKokkos::post_force(int vflag) // virial setup - v_init(vflag); + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); - // reallocate per-atom arrays if necessary + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); - memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"improper:vatom"); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"shake:vatom"); d_vatom = k_vatom.template view(); } diff --git a/src/KOKKOS/fix_wall_harmonic_kokkos.cpp b/src/KOKKOS/fix_wall_harmonic_kokkos.cpp index 49b9dc299f5..fb3c1a90b29 100644 --- a/src/KOKKOS/fix_wall_harmonic_kokkos.cpp +++ b/src/KOKKOS/fix_wall_harmonic_kokkos.cpp @@ -55,14 +55,27 @@ void FixWallHarmonicKokkos::precompute(int /*m_in*/) /* ---------------------------------------------------------------------- */ template -void FixWallHarmonicKokkos::post_force(int vflag) +void FixWallHarmonicKokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary + if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom, vatom); memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "wall_harmonic:vatom"); d_vatom = k_vatom.template view(); } +} +/* ---------------------------------------------------------------------- */ + +template +void FixWallHarmonicKokkos::post_force(int vflag) +{ FixWallHarmonic::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_harmonic_kokkos.h b/src/KOKKOS/fix_wall_harmonic_kokkos.h index 582a4b664a4..97d0f4955ff 100644 --- a/src/KOKKOS/fix_wall_harmonic_kokkos.h +++ b/src/KOKKOS/fix_wall_harmonic_kokkos.h @@ -40,6 +40,7 @@ class FixWallHarmonicKokkos : public FixWallHarmonic { ~FixWallHarmonicKokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_lj1043_kokkos.cpp b/src/KOKKOS/fix_wall_lj1043_kokkos.cpp index ba2b5c34e4e..d18baf32028 100644 --- a/src/KOKKOS/fix_wall_lj1043_kokkos.cpp +++ b/src/KOKKOS/fix_wall_lj1043_kokkos.cpp @@ -120,14 +120,27 @@ void FixWallLJ1043Kokkos::precompute(int m_in) /* ---------------------------------------------------------------------- */ template -void FixWallLJ1043Kokkos::post_force(int vflag) +void FixWallLJ1043Kokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary + if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom, vatom); memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "wall_lj1043:vatom"); d_vatom = k_vatom.template view(); } +} +/* ---------------------------------------------------------------------- */ + +template +void FixWallLJ1043Kokkos::post_force(int vflag) +{ FixWallLJ1043::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_lj1043_kokkos.h b/src/KOKKOS/fix_wall_lj1043_kokkos.h index 31ed5e97a49..d254e86b31e 100644 --- a/src/KOKKOS/fix_wall_lj1043_kokkos.h +++ b/src/KOKKOS/fix_wall_lj1043_kokkos.h @@ -40,6 +40,7 @@ class FixWallLJ1043Kokkos : public FixWallLJ1043 { ~FixWallLJ1043Kokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_lj126_kokkos.cpp b/src/KOKKOS/fix_wall_lj126_kokkos.cpp index 1ced7256191..a2704bc4010 100644 --- a/src/KOKKOS/fix_wall_lj126_kokkos.cpp +++ b/src/KOKKOS/fix_wall_lj126_kokkos.cpp @@ -98,14 +98,27 @@ void FixWallLJ126Kokkos::precompute(int m_in) /* ---------------------------------------------------------------------- */ template -void FixWallLJ126Kokkos::post_force(int vflag) +void FixWallLJ126Kokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary + if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom, vatom); memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "wall_lj126:vatom"); d_vatom = k_vatom.template view(); } +} +/* ---------------------------------------------------------------------- */ + +template +void FixWallLJ126Kokkos::post_force(int vflag) +{ FixWallLJ126::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_lj126_kokkos.h b/src/KOKKOS/fix_wall_lj126_kokkos.h index e86dca95065..6120790419d 100644 --- a/src/KOKKOS/fix_wall_lj126_kokkos.h +++ b/src/KOKKOS/fix_wall_lj126_kokkos.h @@ -40,6 +40,7 @@ class FixWallLJ126Kokkos : public FixWallLJ126 { ~FixWallLJ126Kokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_lj93_kokkos.cpp b/src/KOKKOS/fix_wall_lj93_kokkos.cpp index d63ebfc4caf..f7c4c6f8539 100644 --- a/src/KOKKOS/fix_wall_lj93_kokkos.cpp +++ b/src/KOKKOS/fix_wall_lj93_kokkos.cpp @@ -105,17 +105,27 @@ void FixWallLJ93Kokkos::precompute(int m) /* ---------------------------------------------------------------------- */ template -void FixWallLJ93Kokkos::post_force(int vflag) +void FixWallLJ93Kokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); - // reallocate per-atom arrays if necessary + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"wall_lj93:vatom"); d_vatom = k_vatom.template view(); } +} + +/* ---------------------------------------------------------------------- */ +template +void FixWallLJ93Kokkos::post_force(int vflag) +{ FixWallLJ93::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_lj93_kokkos.h b/src/KOKKOS/fix_wall_lj93_kokkos.h index 51bdd7dbd2c..334348490a5 100644 --- a/src/KOKKOS/fix_wall_lj93_kokkos.h +++ b/src/KOKKOS/fix_wall_lj93_kokkos.h @@ -40,6 +40,7 @@ class FixWallLJ93Kokkos : public FixWallLJ93 { ~FixWallLJ93Kokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_morse_kokkos.cpp b/src/KOKKOS/fix_wall_morse_kokkos.cpp index 1f547f381f0..396422eff4a 100644 --- a/src/KOKKOS/fix_wall_morse_kokkos.cpp +++ b/src/KOKKOS/fix_wall_morse_kokkos.cpp @@ -100,14 +100,27 @@ void FixWallMorseKokkos::precompute(int m_in) /* ---------------------------------------------------------------------- */ template -void FixWallMorseKokkos::post_force(int vflag) +void FixWallMorseKokkos::v_setup_peratom(int vflag) { + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary + if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom, vatom); memoryKK->create_kokkos(k_vatom, vatom, maxvatom, "wall_morse:vatom"); d_vatom = k_vatom.template view(); } +} +/* ---------------------------------------------------------------------- */ + +template +void FixWallMorseKokkos::post_force(int vflag) +{ FixWallMorse::post_force(vflag); if (vflag_atom) { diff --git a/src/KOKKOS/fix_wall_morse_kokkos.h b/src/KOKKOS/fix_wall_morse_kokkos.h index 3495549e35d..f1cf9197578 100644 --- a/src/KOKKOS/fix_wall_morse_kokkos.h +++ b/src/KOKKOS/fix_wall_morse_kokkos.h @@ -40,6 +40,7 @@ class FixWallMorseKokkos : public FixWallMorse { ~FixWallMorseKokkos() override; void precompute(int) override; void post_force(int) override; + void v_setup_peratom(int) override; void wall_particle(int, int, double) override; int m; diff --git a/src/KOKKOS/fix_wall_region_kokkos.cpp b/src/KOKKOS/fix_wall_region_kokkos.cpp index cc8d876d043..5c062b4ffc8 100644 --- a/src/KOKKOS/fix_wall_region_kokkos.cpp +++ b/src/KOKKOS/fix_wall_region_kokkos.cpp @@ -78,9 +78,12 @@ void FixWallRegionKokkos::post_force(int vflag) // virial setup - v_init(vflag); + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) - // reallocate per-atom arrays if necessary + v_init(vflag,0); + + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); @@ -96,10 +99,6 @@ void FixWallRegionKokkos::post_force(int vflag) region->prematch(); - // virial setup - - v_init(vflag); - // region->match() ensures particle is in region or on surface, else error // if returned contact dist r = 0, is on surface, also an error // in COLLOID case, r <= radius is an error diff --git a/src/fix.cpp b/src/fix.cpp index 884ffa854be..e3cc98a6cca 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -203,7 +203,7 @@ void Fix::set_molecule(int, tagint, int, double *, double *, double *) energy is *much* faster. ------------------------------------------------------------------------- */ -void Fix::ev_setup(int eflag, int vflag) +void Fix::ev_setup(int eflag, int vflag, int alloc) { int i,n; @@ -234,18 +234,24 @@ void Fix::ev_setup(int eflag, int vflag) if (eflag_atom && atom->nlocal > maxeatom) { maxeatom = atom->nmax; - memory->destroy(eatom); - memory->create(eatom,maxeatom,"fix:eatom"); + if (alloc) { + memory->destroy(eatom); + memory->create(eatom,maxeatom,"fix:eatom"); + } } if (vflag_atom && atom->nlocal > maxvatom) { maxvatom = atom->nmax; - memory->destroy(vatom); - memory->create(vatom,maxvatom,6,"fix:vatom"); + if (alloc) { + memory->destroy(vatom); + memory->create(vatom,maxvatom,6,"fix:vatom"); + } } if (cvflag_atom && atom->nlocal > maxcvatom) { maxcvatom = atom->nmax; - memory->destroy(cvatom); - memory->create(cvatom,maxcvatom,9,"fix:cvatom"); + if (alloc) { + memory->destroy(cvatom); + memory->create(cvatom,maxcvatom,9,"fix:cvatom"); + } } // zero accumulators @@ -253,11 +259,11 @@ void Fix::ev_setup(int eflag, int vflag) // fixes tally it individually via fix_modify energy yes and compute_scalar() if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (eflag_atom) { + if (eflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) eatom[i] = 0.0; } - if (vflag_atom) { + if (vflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; @@ -268,7 +274,7 @@ void Fix::ev_setup(int eflag, int vflag) vatom[i][5] = 0.0; } } - if (cvflag_atom) { + if (cvflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) { cvatom[i][0] = 0.0; @@ -291,7 +297,7 @@ void Fix::ev_setup(int eflag, int vflag) if thermo_virial is not set, virial tallying is disabled ------------------------------------------------------------------------- */ -void Fix::v_setup(int vflag) +void Fix::v_setup(int vflag, int alloc) { int i,n; @@ -310,19 +316,23 @@ void Fix::v_setup(int vflag) if (vflag_atom && atom->nlocal > maxvatom) { maxvatom = atom->nmax; - memory->destroy(vatom); - memory->create(vatom,maxvatom,6,"fix:vatom"); + if (alloc) { + memory->destroy(vatom); + memory->create(vatom,maxvatom,6,"fix:vatom"); + } } if (cvflag_atom && atom->nlocal > maxcvatom) { maxcvatom = atom->nmax; - memory->destroy(cvatom); - memory->create(cvatom,maxcvatom,9,"fix:cvatom"); + if (alloc) { + memory->destroy(cvatom); + memory->create(cvatom,maxcvatom,9,"fix:cvatom"); + } } // zero accumulators if (vflag_global) for (i = 0; i < 6; i++) virial[i] = 0.0; - if (vflag_atom) { + if (vflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) { vatom[i][0] = 0.0; @@ -333,7 +343,7 @@ void Fix::v_setup(int vflag) vatom[i][5] = 0.0; } } - if (cvflag_atom) { + if (cvflag_atom && alloc) { n = atom->nlocal; for (i = 0; i < n; i++) { cvatom[i][0] = 0.0; diff --git a/src/fix.h b/src/fix.h index 5f259a666c3..184b5409d21 100644 --- a/src/fix.h +++ b/src/fix.h @@ -282,25 +282,31 @@ class Fix : protected Pointers { int dynamic; // recount atoms for temperature computes - void ev_init(int eflag, int vflag) + // alloc = 0 tells ev_setup()/v_setup() to only update the flags and the + // maxeatom/maxvatom sizes, but not to allocate or zero the plain per-atom + // arrays. Styles that manage eatom/vatom themselves (e.g. the KOKKOS + // variants, which store them in dual views) must pass alloc = 0, since + // otherwise the plain arrays allocated here are orphaned by the style. + + void ev_init(int eflag, int vflag, int alloc = 1) { if ((eflag && thermo_energy) || (vflag && thermo_virial)) - ev_setup(eflag, vflag); + ev_setup(eflag, vflag, alloc); else evflag = eflag_either = eflag_global = eflag_atom = eflag_only = vflag_either = vflag_global = vflag_atom = cvflag_atom = 0; } - void ev_setup(int, int); + void ev_setup(int, int, int alloc = 1); void ev_tally(int, int *, double, double, double *); - void v_init(int vflag) + void v_init(int vflag, int alloc = 1) { if (vflag && thermo_virial) - v_setup(vflag); + v_setup(vflag, alloc); else evflag = vflag_either = vflag_global = vflag_atom = cvflag_atom = 0; } - void v_setup(int); + void v_setup(int, int alloc = 1); void v_tally(int, int *, double, double *); void v_tally(int, int *, double, double *, int, int, int[][2], double *, double[][3]); void v_tally(int, int *, double, double *, double[][3], double[][3], double[]); diff --git a/src/fix_wall.cpp b/src/fix_wall.cpp index e104714db88..47b0e791435 100644 --- a/src/fix_wall.cpp +++ b/src/fix_wall.cpp @@ -465,7 +465,7 @@ void FixWall::post_force(int vflag) { // virial setup - v_init(vflag); + v_setup_peratom(vflag); // energy intialize. // eflag is used to track whether wall energies have been communicated. diff --git a/src/fix_wall.h b/src/fix_wall.h index f863ccf182b..dd2b3aa84a9 100644 --- a/src/fix_wall.h +++ b/src/fix_wall.h @@ -47,6 +47,14 @@ class FixWall : public Fix { virtual void precompute(int) = 0; virtual void wall_particle(int, int, double) = 0; + + // set up the per-atom virial storage at the beginning of post_force(). + // accelerator styles that keep the per-atom virial in their own arrays + // override this and pass alloc = 0 to v_init(), so that the plain + // base-class vatom array is not allocated behind their back. + + virtual void v_setup_peratom(int vflag) { v_init(vflag); } + static void update_image_plane(int, int, double, double **, class Domain *); protected: From ac6f7634bd47c4120e29562f05a34f089f5ca06a Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 21:28:53 +0000 Subject: [PATCH 34/43] KOKKOS: fix missing virial, every keyword and sforce sizing in fix addforce/kk Three defects in fix addforce/kk, found while auditing the per-atom virial handling of the KOKKOS fixes: - The kernels never called v_tally(), which was therefore dead code. Neither the global nor the per-atom virial of the added force was tallied, so the pressure was wrong and compute stress/atom returned zero for the fix contribution. - The every keyword was ignored: post_force() was missing the ntimestep % nevery test, so the force was added on every step. - The sforce array was only grown for varflag == ATOM. With an atom-style energy variable and equal-style force components (varflag == EQUAL, estyle == ATOM) it kept its initial length of one atom while compute_atom() wrote nlocal entries into it, and the values it did receive were never synced to the device. Tally the virial as FixAddForce::post_force() does, honor nevery, and grow and sync sforce for estyle == ATOM as well. Also correct a copy-and-paste memory label. Constant force, equal-style and atom-style force variables, atom-style energy variables and every N now all reproduce the pressure, energy and per-atom stress of fix addforce exactly, and valgrind reports no invalid accesses for the atom-style energy variable case, which previously overflowed the sforce allocation. --- src/KOKKOS/fix_addforce_kokkos.cpp | 42 ++++++++++++++++++++++++++---- 1 file changed, 37 insertions(+), 5 deletions(-) diff --git a/src/KOKKOS/fix_addforce_kokkos.cpp b/src/KOKKOS/fix_addforce_kokkos.cpp index 8944eae61a7..902b6e5cfcb 100644 --- a/src/KOKKOS/fix_addforce_kokkos.cpp +++ b/src/KOKKOS/fix_addforce_kokkos.cpp @@ -74,6 +74,8 @@ void FixAddForceKokkos::init() template void FixAddForceKokkos::post_force(int vflag) { + if (update->ntimestep % nevery) return; + atomKK->sync(execution_space, X_MASK | F_MASK | IMAGE_MASK | MASK_MASK); x = atomKK->k_x.view(); @@ -84,13 +86,16 @@ void FixAddForceKokkos::post_force(int vflag) // virial setup - v_init(vflag); + // the per-atom virial is accumulated into a dual view, so the plain + // base-class vatom array must not be allocated here (alloc = 0) + + v_init(vflag,0); - // reallocate per-atom arrays if necessary + // reallocate the per-atom virial dual view if necessary if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); - memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"efield:vatom"); + memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"addforce:vatom"); d_vatom = k_vatom.template view(); } @@ -109,7 +114,7 @@ void FixAddForceKokkos::post_force(int vflag) // reallocate sforce array if necessary - if (varflag == ATOM && atom->nmax > maxatom) { + if (((varflag == ATOM) || (estyle == ATOM)) && atom->nmax > maxatom) { maxatom = atom->nmax; memoryKK->destroy_kokkos(k_sforce,sforce); memoryKK->create_kokkos(k_sforce,sforce,maxatom,4,"addforce:sforce"); @@ -149,7 +154,9 @@ void FixAddForceKokkos::post_force(int vflag) modify->addstep_compute(update->ntimestep + 1); - if (varflag == ATOM) { // this can be removed when variable class is ported to Kokkos + // this can be removed when the variable class is ported to Kokkos + + if ((varflag == ATOM) || (estyle == ATOM)) { k_sforce.modify_host(); k_sforce.sync(); } @@ -204,6 +211,17 @@ void FixAddForceKokkos::operator()(TagFixAddForceConstant, const int if (xstyle) f(i,0) += static_cast(xvalue_kk); if (ystyle) f(i,1) += static_cast(yvalue_kk); if (zstyle) f(i,2) += static_cast(zvalue_kk); + + if (evflag) { + KK_FLOAT v[6]; + v[0] = static_cast(xvalue * unwrapKK[0]); + v[1] = static_cast(yvalue * unwrapKK[1]); + v[2] = static_cast(zvalue * unwrapKK[2]); + v[3] = static_cast(xvalue * unwrapKK[1]); + v[4] = static_cast(xvalue * unwrapKK[2]); + v[5] = static_cast(yvalue * unwrapKK[2]); + v_tally(result,i,v); + } } } @@ -242,6 +260,20 @@ void FixAddForceKokkos::operator()(TagFixAddForceNonConstant, const else if (ystyle) f(i,1) += static_cast(yvalue_kk); if (zstyle == ATOM) f(i,2) += static_cast(d_sforce(i,2)); else if (zstyle) f(i,2) += static_cast(zvalue_kk); + + if (evflag) { + const double xv = (xstyle == ATOM) ? static_cast(d_sforce(i,0)) : xvalue; + const double yv = (ystyle == ATOM) ? static_cast(d_sforce(i,1)) : yvalue; + const double zv = (zstyle == ATOM) ? static_cast(d_sforce(i,2)) : zvalue; + KK_FLOAT v[6]; + v[0] = xstyle ? static_cast(xv * unwrapKK[0]) : static_cast(0.0); + v[1] = ystyle ? static_cast(yv * unwrapKK[1]) : static_cast(0.0); + v[2] = zstyle ? static_cast(zv * unwrapKK[2]) : static_cast(0.0); + v[3] = xstyle ? static_cast(xv * unwrapKK[1]) : static_cast(0.0); + v[4] = xstyle ? static_cast(xv * unwrapKK[2]) : static_cast(0.0); + v[5] = ystyle ? static_cast(yv * unwrapKK[2]) : static_cast(0.0); + v_tally(result,i,v); + } } } From 988b96350ac47553ebc378d3bc4c7f3c0a28a8e9 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 25 Aug 2026 15:31:10 -0600 Subject: [PATCH 35/43] Add missing data transfer in fix shake/kk --- src/KOKKOS/fix_shake_kokkos.cpp | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index 734f5f69869..faa0235b8b3 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -152,6 +152,15 @@ FixShakeKokkos::~FixShakeKokkos() template void FixShakeKokkos::init() { + // FixShake::init() reads shake_flag/shake_atom/shake_type through the + // plain host pointers to recompute angle_distance[]; comm->exchange() + // updates the device copies only, so the host side must be brought + // current before FixShake::init() can see this run's cluster data + + k_shake_flag.sync_host(); + k_shake_atom.sync_host(); + k_shake_type.sync_host(); + FixShake::init(); if (utils::strmatch(update->integrate_style,"^respa")) From 257b8364f37b427e9adbc2147f2d5513b01fafb7 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Wed, 26 Aug 2026 13:00:36 +0000 Subject: [PATCH 36/43] KOKKOS: compute the per-atom energy and virial in fix shake/kk fix shake/kk reported zero per-atom energy and virial where fix shake reports real values, and the per-atom virial was additionally lost whenever more than one thread was used. - min_post_force() only accumulated the global virial. Tally the per-atom energy and virial of the substituted bond restraints as FixShake::bond_force() does, through a new device ev_tally() that mirrors Fix::ev_tally(). The energy and virial of a restraint are shared out over the atoms of the pair that this processor owns, so the image indices the kernel uses for the positions are first mapped back onto their owners; a cluster reaching across a periodic boundary is then credited to both atoms, as on the CPU. - post_force() and min_post_force() never contributed the duplicated per-atom virial scatter view back into the dual view, so with more than one OpenMP thread compute stress/atom saw only zeros. - min_post_force() dispatched the kernel on vflag rather than evflag, which skipped the reduction, and with it the restraint energy, whenever only per-atom energy was requested. Verified against the non-accelerated style with 1, 2 and 4 OpenMP threads, for both the MD and the minimization path: global pressure, global energy, and the per-atom energy and stress of every atom agree to round-off (largest relative deviation 3e-11), and valgrind reports no leaks or invalid accesses. The Kokkos kernels do not compute the centroid virial of the constraint forces, but the style inherited centroidstressflag = CENTROID_AVAIL from FixShake and so silently contributed zeros to compute centroid/stress/atom. Set CENTROID_NOTAVAIL, as the KOKKOS angle and dihedral styles do, so that combination is rejected with a clear error instead, and document the restriction. --- doc/src/compute_stress_atom.rst | 2 + doc/src/fix_shake.rst | 8 ++ src/KOKKOS/fix_shake_kokkos.cpp | 179 +++++++++++++++++++++++++------- src/KOKKOS/fix_shake_kokkos.h | 15 ++- 4 files changed, 161 insertions(+), 43 deletions(-) diff --git a/doc/src/compute_stress_atom.rst b/doc/src/compute_stress_atom.rst index 6c4e0b690c7..15e5fe4539a 100644 --- a/doc/src/compute_stress_atom.rst +++ b/doc/src/compute_stress_atom.rst @@ -276,6 +276,8 @@ i.e. the kspace_style command in LAMMPS. It also does not implement the following fixes which add rigid-body constraints: :doc:`fix rigid/* ` and the OpenMP accelerated version of :doc:`fix rigid/small `, while all other :doc:`fix rigid/*/small ` are implemented. +The KOKKOS version of :doc:`fix shake ` is not implemented +either, while the non-accelerated version is. LAMMPS will generate an error if one of these options is included in your model. Extension of centroid stress calculations to these force diff --git a/doc/src/fix_shake.rst b/doc/src/fix_shake.rst index 268135d945a..0ef293a8ee8 100644 --- a/doc/src/fix_shake.rst +++ b/doc/src/fix_shake.rst @@ -304,6 +304,14 @@ can make minimization very inefficient and also cause stability problems with some minimization algorithms. Sometimes those can be avoided by reducing the :doc:`timestep `. +.. versionchanged:: TBD + +The *shake/kk* style does not compute the centroid virial of the +constraint forces, so it cannot be used with :doc:`compute +centroid/stress/atom `. Use the non-accelerated +*shake* style for that. Previously this combination was accepted but +silently contributed zero. + Related commands """""""""""""""" diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index 0c38c2be666..c3174f13c3f 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -49,6 +49,12 @@ FixShakeKokkos::FixShakeKokkos(LAMMPS *lmp, int narg, char **arg) : if (store_flag) error->all(FLERR, "Option 'store yes' is not (yet) supported by fix {}/kk", style); + // the centroid virial of the SHAKE constraint forces is not (yet) computed + // by the Kokkos kernels, so do not claim it is available. compute + // centroid/stress/atom then errors out instead of silently reporting zeros. + + centroidstressflag = CENTROID_NOTAVAIL; + datamask_read = EMPTY_MASK; datamask_modify = EMPTY_MASK; @@ -141,6 +147,7 @@ FixShakeKokkos::~FixShakeKokkos() memoryKK->destroy_kokkos(k_list,list); memoryKK->destroy_kokkos(k_closest_list,closest_list); + memoryKK->destroy_kokkos(k_eatom,eatom); memoryKK->destroy_kokkos(k_vatom,vatom); } @@ -351,10 +358,13 @@ void FixShakeKokkos::min_post_force(int vflag) ev_init(eflag, vflag, 0); - // reallocate the per-atom virial dual view if necessary. minimization - // only tallies the global virial, but the freshly created view keeps the - // per-atom virial zeroed for computes reading it + // reallocate the per-atom energy and virial dual views if necessary + if (eflag_atom) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"shake:eatom"); + d_eatom = k_eatom.template view(); + } if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"shake:vatom"); @@ -362,15 +372,29 @@ void FixShakeKokkos::min_post_force(int vflag) } ebond = 0.0; + nlocal = atomKK->nlocal; atomKK->sync(execution_space, X_MASK | F_MASK); k_shake_flag.sync(); + k_shake_atom.sync(); k_shake_type.sync(); k_list.sync(); k_closest_list.sync(); k_bond_distance.sync(); k_angle_distance.sync(); + // the restraint energy and virial are credited to the owner of each atom, + // which requires mapping image atoms back through the atom map + + map_style = atom->map_style; + if (map_style == Atom::MAP_ARRAY) { + k_map_array = atomKK->k_map_array; + k_map_array.template sync(); + } else if (map_style == Atom::MAP_HASH) { + k_map_hash = atomKK->k_map_hash; + k_map_hash.template sync(); + } + // Assign class member views for Tagged Kernel access this->d_x = atomKK->k_x.view(); this->d_f = atomKK->k_f.view(); @@ -384,10 +408,15 @@ void FixShakeKokkos::min_post_force(int vflag) if (neighflag != HALF) need_dup = std::is_same_v, Kokkos::Experimental::ScatterDuplicated>; - if (need_dup) - dup_f = Kokkos::Experimental::create_scatter_view(d_f); - else - ndup_f = Kokkos::Experimental::create_scatter_view(d_f); + if (need_dup) { + dup_f = Kokkos::Experimental::create_scatter_view(d_f); + dup_eatom = Kokkos::Experimental::create_scatter_view(d_eatom); + dup_vatom = Kokkos::Experimental::create_scatter_view(d_vatom); + } else { + ndup_f = Kokkos::Experimental::create_scatter_view(d_f); + ndup_eatom = Kokkos::Experimental::create_scatter_view(d_eatom); + ndup_vatom = Kokkos::Experimental::create_scatter_view(d_vatom); + } copymode = 1; @@ -411,20 +440,27 @@ void FixShakeKokkos::min_post_force(int vflag) EV_FLOAT ev; if (neighflag == HALF) { - if(vflag) + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy>(0, nlist), *this, ev); else - Kokkos::parallel_for(Kokkos::RangePolicy>(0, nlist), *this); + Kokkos::parallel_reduce(Kokkos::RangePolicy>(0, nlist), *this, ev); } else { - if(vflag) + if (evflag) Kokkos::parallel_reduce(Kokkos::RangePolicy>(0, nlist), *this, ev); else - Kokkos::parallel_for(Kokkos::RangePolicy>(0, nlist), *this); + Kokkos::parallel_reduce(Kokkos::RangePolicy>(0, nlist), *this, ev); } copymode = 0; - if (need_dup) Kokkos::Experimental::contribute(d_f, dup_f); + // reduction over duplicated memory + + if (need_dup) { + Kokkos::Experimental::contribute(d_f, dup_f); + if (eflag_atom) Kokkos::Experimental::contribute(d_eatom, dup_eatom); + if (vflag_atom) Kokkos::Experimental::contribute(d_vatom, dup_vatom); + } + comm->reverse_comm(this); this->ebond = static_cast(ev.evdwl); @@ -437,22 +473,37 @@ void FixShakeKokkos::min_post_force(int vflag) virial[5] += static_cast(ev.v[5]); } + if (eflag_atom) { + k_eatom.template modify(); + k_eatom.sync_host(); + } + if (vflag_atom) { + k_vatom.template modify(); + k_vatom.sync_host(); + } + atomKK->modified(execution_space, F_MASK); + // free duplicated memory + + if (need_dup) { + dup_f = {}; + dup_eatom = {}; + dup_vatom = {}; + } + if (update->ntimestep == next_output) { atomKK->modified(execution_space, X_MASK); stats(); } - - if (need_dup) dup_f = {}; } /* ---------------------------------------------------------------------- */ template -template +template KOKKOS_INLINE_FUNCTION -void FixShakeKokkos::operator()(TagFixShakeMinPostForce, const int &i, EV_FLOAT &ev) const +void FixShakeKokkos::operator()(TagFixShakeMinPostForce, const int &i, EV_FLOAT &ev) const { auto v_f = ScatterViewHelper, decltype(dup_f), decltype(ndup_f)>::get(dup_f, ndup_f); auto a_f = v_f.template access>(); @@ -461,7 +512,12 @@ void FixShakeKokkos::operator()(TagFixShakeMinPostForce::operator()(TagFixShakeMinPostForce(delx * fbond); a_f(idx1, 1) -= static_cast(dely * fbond); a_f(idx1, 2) -= static_cast(delz * fbond); - ev.evdwl += static_cast(eb); - if (VFLAG) { - ev.v[0] += static_cast(static_cast(0.5) * delx * delx * fbond); - ev.v[1] += static_cast(static_cast(0.5) * dely * dely * fbond); - ev.v[2] += static_cast(static_cast(0.5) * delz * delz * fbond); - ev.v[3] += static_cast(static_cast(0.5) * delx * dely * fbond); - ev.v[4] += static_cast(static_cast(0.5) * delx * delz * fbond); - ev.v[5] += static_cast(static_cast(0.5) * dely * delz * fbond); + + // energy and virial are shared out over the owned atoms of the pair, as + // FixShake::bond_force() does. the closest image of an atom owned by + // this processor is mapped back onto the owner first, so that a cluster + // reaching across a periodic boundary is still credited to both atoms. + + int atomlist[2]; + int count = 0; + const int own0 = AtomKokkos::map_kokkos(d_shake_atom(m,slot0),map_style, + k_map_array,k_map_hash); + const int own1 = AtomKokkos::map_kokkos(d_shake_atom(m,slot1),map_style, + k_map_array,k_map_hash); + if ((own0 >= 0) && (own0 < nlocal)) atomlist[count++] = own0; + if ((own1 >= 0) && (own1 < nlocal)) atomlist[count++] = own1; + + const KK_FLOAT total = static_cast(2.0); + ev.evdwl += static_cast((static_cast(count)/total) * eb); + + if (EVFLAG) { + KK_FLOAT v[6]; + v[0] = static_cast(0.5) * delx * delx * fbond; + v[1] = static_cast(0.5) * dely * dely * fbond; + v[2] = static_cast(0.5) * delz * delz * fbond; + v[3] = static_cast(0.5) * delx * dely * fbond; + v[4] = static_cast(0.5) * delx * delz * fbond; + v[5] = static_cast(0.5) * dely * delz * fbond; + ev_tally(ev,count,atomlist,total,eb,v); } if (output_every && !is_angle) { Kokkos::atomic_add(&d_b_stats(type_idx, 0), 1.0); @@ -499,20 +574,20 @@ void FixShakeKokkos::operator()(TagFixShakeMinPostForce(2.0)*r1*r2)) * static_cast(180.0)/static_cast(MY_PI); int mt = d_shake_type(m, 2); @@ -528,11 +603,11 @@ void FixShakeKokkos::operator()(TagFixShakeMinPostForce -template +template KOKKOS_INLINE_FUNCTION -void FixShakeKokkos::operator()(TagFixShakeMinPostForce, const int &i) const { +void FixShakeKokkos::operator()(TagFixShakeMinPostForce, const int &i) const { EV_FLOAT ev; - this->template operator()(TagFixShakeMinPostForce(), i, ev); + this->template operator()(TagFixShakeMinPostForce(), i, ev); } /* ---------------------------------------------------------------------- @@ -657,8 +732,10 @@ void FixShakeKokkos::post_force(int vflag) // reduction over duplicated memory - if (need_dup) + if (need_dup) { Kokkos::Experimental::contribute(d_f,dup_f); + if (vflag_atom) Kokkos::Experimental::contribute(d_vatom,dup_vatom); + } atomKK->modified(execution_space,F_MASK); @@ -2201,6 +2278,32 @@ void FixShakeKokkos::correct_coordinates(int vflag) { atomKK->modified(Host,X_MASK|V_MASK|F_MASK); } +/* ---------------------------------------------------------------------- + tally energy and virial into global and per-atom accumulators + n = # of local owned atoms involved, with local indices in list + eng = total energy for the interaction involving total atoms + increment per-atom energy of each atom in list by 1/total fraction + mirrors Fix::ev_tally() +------------------------------------------------------------------------- */ + +template +template +// NOLINTNEXTLINE +KOKKOS_INLINE_FUNCTION +void FixShakeKokkos::ev_tally(EV_FLOAT &ev, int n, int *atomlist, KK_FLOAT total, + KK_FLOAT eng, KK_FLOAT *v) const +{ + if (eflag_atom) { + auto v_eatom = ScatterViewHelper,decltype(dup_eatom),decltype(ndup_eatom)>::get(dup_eatom,ndup_eatom); + auto a_eatom = v_eatom.template access>(); + const KK_FLOAT fraction = eng/total; + for (int i = 0; i < n; i++) + a_eatom(atomlist[i]) += static_cast(fraction); + } + + v_tally(ev,n,atomlist,total,v); +} + /* ---------------------------------------------------------------------- tally virial into global and per-atom accumulators n = # of local owned atoms involved, with local indices in list diff --git a/src/KOKKOS/fix_shake_kokkos.h b/src/KOKKOS/fix_shake_kokkos.h index 13a4dd6f896..afea28d331f 100644 --- a/src/KOKKOS/fix_shake_kokkos.h +++ b/src/KOKKOS/fix_shake_kokkos.h @@ -32,7 +32,7 @@ namespace LAMMPS_NS { struct TagFixShakePreNeighbor{}; -template +template struct TagFixShakeMinPostForce{}; template @@ -86,13 +86,13 @@ class FixShakeKokkos : public FixShake, public KokkosBase { KOKKOS_INLINE_FUNCTION void operator()(TagFixShakePreNeighbor, const int&) const; - template + template KOKKOS_INLINE_FUNCTION - void operator()(TagFixShakeMinPostForce, const int&, EV_FLOAT&) const; + void operator()(TagFixShakeMinPostForce, const int&, EV_FLOAT&) const; - template + template KOKKOS_INLINE_FUNCTION - void operator()(TagFixShakeMinPostForce, const int&) const; + void operator()(TagFixShakeMinPostForce, const int&) const; template // NOLINTNEXTLINE @@ -234,6 +234,11 @@ class FixShakeKokkos : public FixShake, public KokkosBase { KOKKOS_INLINE_FUNCTION void v_tally(EV_FLOAT&, int, int *, KK_FLOAT, KK_FLOAT *) const; + template +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void ev_tally(EV_FLOAT&, int, int *, KK_FLOAT, KK_FLOAT, KK_FLOAT *) const; + int first,nsend; typename AT::t_int_1d d_sendlist; From 67c3ddbfb820cdfeb3cb9ccd9d75729a134ebe4b Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Wed, 26 Aug 2026 13:18:13 +0000 Subject: [PATCH 37/43] KOKKOS: fix data race on the region contact list RegBlockKokkos and RegSphereKokkos kept the list of contacts between a particle and the region surface in a single view owned by the region. surface_kokkos() is called from inside a parallel region, once per particle, so every thread wrote into the same list and read back whichever contacts happened to survive. fix wall/region/kk therefore computed wrong forces as soon as more than one thread was used, and the result was not reproducible from run to run. With four OpenMP threads, 637 of the 1827 per-atom force and stress components of a small test differed from the non-accelerated style, by up to 170%. Pass the contact list in from the caller instead, so that it can live on the stack of the thread that uses it, and add a MAXCONTACT constant to each region class for its length. surface_kokkos(), surface_interior_kokkos(), surface_exterior_kokkos() and the helpers they call no longer modify the region and are now const, so that a shared buffer cannot be reintroduced unnoticed. The device-side contact views and the destructors that freed them are gone. Verified against the non-accelerated style with 1, 2 and 4 OpenMP threads: forces and per-atom stresses of every atom are now bit for bit identical for a block and a sphere region, inside and outside, with open faces, with a moving region, and for the harmonic wall style, and valgrind reports no leaks or invalid accesses. The same race affects the GPU backends, which could not be tested here. --- src/KOKKOS/fix_wall_region_kokkos.cpp | 19 +++-- src/KOKKOS/region_block_kokkos.cpp | 11 --- src/KOKKOS/region_block_kokkos.h | 117 +++++++++++++------------- src/KOKKOS/region_sphere_kokkos.cpp | 11 --- src/KOKKOS/region_sphere_kokkos.h | 73 ++++++++-------- 5 files changed, 112 insertions(+), 119 deletions(-) diff --git a/src/KOKKOS/fix_wall_region_kokkos.cpp b/src/KOKKOS/fix_wall_region_kokkos.cpp index 5c062b4ffc8..167ac394660 100644 --- a/src/KOKKOS/fix_wall_region_kokkos.cpp +++ b/src/KOKKOS/fix_wall_region_kokkos.cpp @@ -28,6 +28,8 @@ #include "region_block_kokkos.h" #include "region_sphere_kokkos.h" +#include + using namespace LAMMPS_NS; using namespace MathSpecialKokkos; @@ -163,14 +165,21 @@ void FixWallRegionKokkos::wall_particle(T regionKK, const int i, val else tooclose = 0.0; - int n = regionKK->surface_kokkos(static_cast(d_x(i,0)), static_cast(d_x(i,1)), static_cast(d_x(i,2)), cutoff); + // the contact list lives on the stack, so that concurrently running + // threads do not overwrite each other's contacts + + Region::Contact contact[std::remove_pointer_t::MAXCONTACT]; + + int n = regionKK->surface_kokkos(static_cast(d_x(i,0)), + static_cast(d_x(i,1)), + static_cast(d_x(i,2)), cutoff, contact); for ( int m = 0; m < n; m++) { - KK_FLOAT r = static_cast(regionKK->d_contact[m].r); - KK_FLOAT delx = static_cast(regionKK->d_contact[m].delx); - KK_FLOAT dely = static_cast(regionKK->d_contact[m].dely); - KK_FLOAT delz = static_cast(regionKK->d_contact[m].delz); + KK_FLOAT r = static_cast(contact[m].r); + KK_FLOAT delx = static_cast(contact[m].delx); + KK_FLOAT dely = static_cast(contact[m].dely); + KK_FLOAT delz = static_cast(contact[m].delz); if (r <= tooclose) Kokkos::abort("Particle outside surface of region used in fix wall/region"); diff --git a/src/KOKKOS/region_block_kokkos.cpp b/src/KOKKOS/region_block_kokkos.cpp index 3802c22f968..336bbe99a5e 100644 --- a/src/KOKKOS/region_block_kokkos.cpp +++ b/src/KOKKOS/region_block_kokkos.cpp @@ -16,7 +16,6 @@ #include "atom_kokkos.h" #include "atom_masks.h" -#include "memory_kokkos.h" using namespace LAMMPS_NS; using namespace MathSpecialKokkos; @@ -28,16 +27,6 @@ RegBlockKokkos::RegBlockKokkos(LAMMPS *lmp, int narg, char **arg) : RegBlock(lmp, narg, arg) { atomKK = (AtomKokkos*) atom; - memoryKK->create_kokkos(d_contact,6,"region_block:d_contact"); -} - -/* ---------------------------------------------------------------------- */ - -template -RegBlockKokkos::~RegBlockKokkos() -{ - if (copymode) return; - memoryKK->destroy_kokkos(d_contact); } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/region_block_kokkos.h b/src/KOKKOS/region_block_kokkos.h index 24285e2448f..191fd386f68 100644 --- a/src/KOKKOS/region_block_kokkos.h +++ b/src/KOKKOS/region_block_kokkos.h @@ -43,8 +43,12 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { typedef DeviceType device_type; typedef ArrayTypes AT; + // maximum number of contacts a particle can have with this region, i.e. the + // length of the contact list the caller has to provide. matches RegBlock::cmax + + static constexpr int MAXCONTACT = 6; + RegBlockKokkos(class LAMMPS *, int, char **); - ~RegBlockKokkos() override; void match_all_kokkos(int, DAT::tdual_int_1d) override; @@ -63,7 +67,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_kokkos(double x, double y, double z, double cutoff) + int surface_kokkos(double x, double y, double z, double cutoff, Contact *contact) const { int ncontact; double xs, ys, zs; @@ -78,32 +82,31 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { if (!openflag) { if (interior) - ncontact = surface_interior_kokkos(xnear, cutoff); + ncontact = surface_interior_kokkos(xnear, cutoff, contact); else - ncontact = surface_exterior_kokkos(xnear, cutoff); + ncontact = surface_exterior_kokkos(xnear, cutoff, contact); } else { // one of surface_int/ext() will return 0 // so no need to worry about offset of contact indices - ncontact = surface_exterior_kokkos(xnear, cutoff) + surface_interior_kokkos(xnear, cutoff); + ncontact = surface_exterior_kokkos(xnear, cutoff, contact) + + surface_interior_kokkos(xnear, cutoff, contact); } if (rotateflag && ncontact) { for (int i = 0; i < ncontact; i++) { - xs = xnear[0] - d_contact[i].delx; - ys = xnear[1] - d_contact[i].dely; - zs = xnear[2] - d_contact[i].delz; + xs = xnear[0] - contact[i].delx; + ys = xnear[1] - contact[i].dely; + zs = xnear[2] - contact[i].delz; forward_transform(xs, ys, zs); - d_contact[i].delx = xorig[0] - xs; - d_contact[i].dely = xorig[1] - ys; - d_contact[i].delz = xorig[2] - zs; + contact[i].delx = xorig[0] - xs; + contact[i].dely = xorig[1] - ys; + contact[i].delz = xorig[2] - zs; } } return ncontact; } - Kokkos::View d_contact; - private: int groupbit; typename AT::t_int_1d d_match; @@ -112,7 +115,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_interior_kokkos(double *x, double cutoff) + int surface_interior_kokkos(double *x, double cutoff, Contact *contact) const { double delta; @@ -126,58 +129,58 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { delta = x[0] - xlo; if (delta < cutoff && !open_faces[0]) { - d_contact[n].r = delta; - d_contact[n].delx = delta; - d_contact[n].dely = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 0; + contact[n].r = delta; + contact[n].delx = delta; + contact[n].dely = contact[n].delz = 0.0; + contact[n].radius = 0; + contact[n].iwall = 0; n++; } delta = xhi - x[0]; if (delta < cutoff && !open_faces[1]) { - d_contact[n].r = delta; - d_contact[n].delx = -delta; - d_contact[n].dely = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 1; + contact[n].r = delta; + contact[n].delx = -delta; + contact[n].dely = contact[n].delz = 0.0; + contact[n].radius = 0; + contact[n].iwall = 1; n++; } delta = x[1] - ylo; if (delta < cutoff && !open_faces[2]) { - d_contact[n].r = delta; - d_contact[n].dely = delta; - d_contact[n].delx = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 2; + contact[n].r = delta; + contact[n].dely = delta; + contact[n].delx = contact[n].delz = 0.0; + contact[n].radius = 0; + contact[n].iwall = 2; n++; } delta = yhi - x[1]; if (delta < cutoff && !open_faces[3]) { - d_contact[n].r = delta; - d_contact[n].dely = -delta; - d_contact[n].delx = d_contact[n].delz = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 3; + contact[n].r = delta; + contact[n].dely = -delta; + contact[n].delx = contact[n].delz = 0.0; + contact[n].radius = 0; + contact[n].iwall = 3; n++; } delta = x[2] - zlo; if (delta < cutoff && !open_faces[4]) { - d_contact[n].r = delta; - d_contact[n].delz = delta; - d_contact[n].delx = d_contact[n].dely = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 4; + contact[n].r = delta; + contact[n].delz = delta; + contact[n].delx = contact[n].dely = 0.0; + contact[n].radius = 0; + contact[n].iwall = 4; n++; } delta = zhi - x[2]; if (delta < cutoff && !open_faces[5]) { - d_contact[n].r = delta; - d_contact[n].delz = -delta; - d_contact[n].delx = d_contact[n].dely = 0.0; - d_contact[n].radius = 0; - d_contact[n].iwall = 5; + contact[n].r = delta; + contact[n].delz = -delta; + contact[n].delx = contact[n].dely = 0.0; + contact[n].radius = 0; + contact[n].iwall = 5; n++; } @@ -186,7 +189,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_exterior_kokkos(double *x, double cutoff) + int surface_exterior_kokkos(double *x, double cutoff, Contact *contact) const { double xp, yp, zp; double xc, yc, zc, dist, mindist; @@ -241,24 +244,24 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { if (mindist == MAXDOUBLEINT) return 0; } - add_contact(0, x, xp, yp, zp); - d_contact[0].iwall = 0; - if (d_contact[0].r < cutoff) return 1; + add_contact(0, x, xp, yp, zp, contact); + contact[0].iwall = 0; + if (contact[0].r < cutoff) return 1; return 0; } // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - void add_contact(int n, double *x, double xp, double yp, double zp) + void add_contact(int n, double *x, double xp, double yp, double zp, Contact *contact) const { double delx = x[0] - xp; double dely = x[1] - yp; double delz = x[2] - zp; - d_contact[n].r = sqrt(delx * delx + dely * dely + delz * delz); - d_contact[n].radius = 0; - d_contact[n].delx = delx; - d_contact[n].dely = dely; - d_contact[n].delz = delz; + contact[n].r = sqrt(delx * delx + dely * dely + delz * delz); + contact[n].radius = 0; + contact[n].delx = delx; + contact[n].dely = dely; + contact[n].delz = delz; } // NOLINTNEXTLINE @@ -325,7 +328,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - void point_on_line_segment(double *a, double *b, double *c, double *d) + void point_on_line_segment(const double *a, const double *b, const double *c, double *d) const { double ba[3], ca[3]; @@ -349,7 +352,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - double inside_face(double *xproj, int iface) + double inside_face(double *xproj, int iface) const { if (iface < 2) { if (xproj[1] > 0 && (xproj[1] < yhi - ylo) && xproj[2] > 0 && (xproj[2] < zhi - zlo)) return 1; @@ -365,7 +368,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - double find_closest_point(int i, double *x, double &xc, double &yc, double &zc) + double find_closest_point(int i, double *x, double &xc, double &yc, double &zc) const { double dot, d2, d2min; double xr[3], xproj[3], p[3]; diff --git a/src/KOKKOS/region_sphere_kokkos.cpp b/src/KOKKOS/region_sphere_kokkos.cpp index b279891e1ba..62f3c463563 100644 --- a/src/KOKKOS/region_sphere_kokkos.cpp +++ b/src/KOKKOS/region_sphere_kokkos.cpp @@ -20,7 +20,6 @@ #include "atom_kokkos.h" #include "atom_masks.h" -#include "memory_kokkos.h" using namespace LAMMPS_NS; @@ -31,16 +30,6 @@ RegSphereKokkos::RegSphereKokkos(LAMMPS *lmp, int narg, char **arg) : RegSphere(lmp, narg, arg) { atomKK = (AtomKokkos*) atom; - memoryKK->create_kokkos(d_contact,1,"region_sphere:d_contact"); -} - -/* ---------------------------------------------------------------------- */ - -template -RegSphereKokkos::~RegSphereKokkos() -{ - if (copymode) return; - memoryKK->destroy_kokkos(d_contact); } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/region_sphere_kokkos.h b/src/KOKKOS/region_sphere_kokkos.h index bde0a9f6193..e05281dea08 100644 --- a/src/KOKKOS/region_sphere_kokkos.h +++ b/src/KOKKOS/region_sphere_kokkos.h @@ -40,8 +40,12 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { typedef DeviceType device_type; typedef ArrayTypes AT; + // maximum number of contacts a particle can have with this region, i.e. the + // length of the contact list the caller has to provide. matches RegSphere::cmax + + static constexpr int MAXCONTACT = 1; + RegSphereKokkos(class LAMMPS *, int, char **); - ~RegSphereKokkos() override; void match_all_kokkos(int, DAT::tdual_int_1d) override; @@ -60,7 +64,7 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_kokkos(double x, double y, double z, double cutoff) + int surface_kokkos(double x, double y, double z, double cutoff, Contact *contact) const { int ncontact; double xs, ys, zs; @@ -73,32 +77,31 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { xnear[0] = x; xnear[1] = y; xnear[2] = z; if (!openflag) { - if (interior) ncontact = surface_interior_kokkos(xnear, cutoff); + if (interior) ncontact = surface_interior_kokkos(xnear, cutoff, contact); else - ncontact = surface_exterior_kokkos(xnear, cutoff); + ncontact = surface_exterior_kokkos(xnear, cutoff, contact); } else { // one of surface_int/ext() will return 0 // so no need to worry about offset of contact indices - ncontact = surface_exterior_kokkos(xnear, cutoff) + surface_interior_kokkos(xnear, cutoff); + ncontact = surface_exterior_kokkos(xnear, cutoff, contact) + + surface_interior_kokkos(xnear, cutoff, contact); } if (rotateflag && ncontact) { for (int i = 0; i < ncontact; i++) { - xs = xnear[0] - d_contact[i].delx; - ys = xnear[1] - d_contact[i].dely; - zs = xnear[2] - d_contact[i].delz; + xs = xnear[0] - contact[i].delx; + ys = xnear[1] - contact[i].dely; + zs = xnear[2] - contact[i].delz; forward_transform(xs, ys, zs); - d_contact[i].delx = xorig[0] - xs; - d_contact[i].dely = xorig[1] - ys; - d_contact[i].delz = xorig[2] - zs; + contact[i].delx = xorig[0] - xs; + contact[i].dely = xorig[1] - ys; + contact[i].delz = xorig[2] - zs; } } return ncontact; } - Kokkos::View d_contact; - private: int groupbit; typename AT::t_int_1d d_match; @@ -107,7 +110,7 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_interior_kokkos(double *x, double cutoff) + int surface_interior_kokkos(double *x, double cutoff, Contact *contact) const { double delx = x[0] - xc; double dely = x[1] - yc; @@ -117,13 +120,13 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { double delta = radius - r; if (delta < cutoff) { - d_contact[0].r = delta; - d_contact[0].delx = delx * (1.0 - radius / r); - d_contact[0].dely = dely * (1.0 - radius / r); - d_contact[0].delz = delz * (1.0 - radius / r); - d_contact[0].radius = -radius; - d_contact[0].iwall = 0; - d_contact[0].varflag = 1; + contact[0].r = delta; + contact[0].delx = delx * (1.0 - radius / r); + contact[0].dely = dely * (1.0 - radius / r); + contact[0].delz = delz * (1.0 - radius / r); + contact[0].radius = -radius; + contact[0].iwall = 0; + contact[0].varflag = 1; return 1; } return 0; @@ -131,7 +134,7 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - int surface_exterior_kokkos(double *x, double cutoff) + int surface_exterior_kokkos(double *x, double cutoff, Contact *contact) const { double delx = x[0] - xc; double dely = x[1] - yc; @@ -141,13 +144,13 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { double delta = r - radius; if (delta < cutoff) { - d_contact[0].r = delta; - d_contact[0].delx = delx * (1.0 - radius / r); - d_contact[0].dely = dely * (1.0 - radius / r); - d_contact[0].delz = delz * (1.0 - radius / r); - d_contact[0].radius = radius; - d_contact[0].iwall = 0; - d_contact[0].varflag = 1; + contact[0].r = delta; + contact[0].delx = delx * (1.0 - radius / r); + contact[0].dely = dely * (1.0 - radius / r); + contact[0].delz = delz * (1.0 - radius / r); + contact[0].radius = radius; + contact[0].iwall = 0; + contact[0].varflag = 1; return 1; } return 0; @@ -155,16 +158,16 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION - void add_contact(int n, double *x, double xp, double yp, double zp) + void add_contact(int n, double *x, double xp, double yp, double zp, Contact *contact) const { double delx = x[0] - xp; double dely = x[1] - yp; double delz = x[2] - zp; - d_contact[n].r = sqrt(delx * delx + dely * dely + delz * delz); - d_contact[n].radius = 0; - d_contact[n].delx = delx; - d_contact[n].dely = dely; - d_contact[n].delz = delz; + contact[n].r = sqrt(delx * delx + dely * dely + delz * delz); + contact[n].radius = 0; + contact[n].delx = delx; + contact[n].dely = dely; + contact[n].delz = delz; } // NOLINTNEXTLINE From 26fbb5759d1b635d5737968f9afe5906751babb7 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 29 Aug 2026 17:53:39 +0000 Subject: [PATCH 38/43] KOKKOS: match Region::match()/surface() coordinate remapping on device Region::match() and Region::surface() map the tested coordinate back into the periodic box before handing it to the region subclass, because not every region handles a shape that reaches past a periodic edge. The KOKKOS regions run the same tests inside a kernel and did not, so any atom that had drifted outside the box since the last reneighboring was tested at its unwrapped position and could be matched differently than on the CPU. With "fix setforce ... region" on a region whose test is sensitive to the wrap, an LJ melt of 2048 atoms saw the force on 62 atoms change in the first timestep and the thermo output diverge by 1e-2 relative after 20 steps. Add region_remap_kokkos.h with a device-side copy of Domain::remap() and the few box quantities it needs, refreshed from prematch() and from match_all_kokkos(). Also bring surface_kokkos() in line with the CPU version for open regions, where the interior test now runs only if the exterior test found no contact. Fix a second, unrelated crash found by the same runs: AtomKokkos::map_delete() frees sametag but, unlike Atom::map_delete(), left max_same at its old value, so the reallocation guard in map_set_device()/map_set_host() did not fire and the next atom map build dereferenced the freed view. This segfaulted on any input that grows the maximum atom ID after the first atom map build, e.g. two create_atoms commands, with "package kokkos atom/map device" and a map array. --- src/KOKKOS/Install.sh | 1 + src/KOKKOS/atom_map_kokkos.cpp | 1 + src/KOKKOS/region_block_kokkos.cpp | 1 + src/KOKKOS/region_block_kokkos.h | 24 +++++- src/KOKKOS/region_remap_kokkos.h | 117 ++++++++++++++++++++++++++++ src/KOKKOS/region_sphere_kokkos.cpp | 1 + src/KOKKOS/region_sphere_kokkos.h | 29 +++++-- 7 files changed, 166 insertions(+), 8 deletions(-) create mode 100644 src/KOKKOS/region_remap_kokkos.h diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 0db125589e7..81515a29e06 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -577,6 +577,7 @@ action rand_pool_wrap_kokkos.cpp action rand_pool_wrap_kokkos.h action region_block_kokkos.cpp action region_block_kokkos.h +action region_remap_kokkos.h action region_sphere_kokkos.cpp action region_sphere_kokkos.h action remap_kokkos.cpp remap.cpp diff --git a/src/KOKKOS/atom_map_kokkos.cpp b/src/KOKKOS/atom_map_kokkos.cpp index 67cf8dfd218..fc636f4f6ca 100644 --- a/src/KOKKOS/atom_map_kokkos.cpp +++ b/src/KOKKOS/atom_map_kokkos.cpp @@ -426,6 +426,7 @@ void AtomKokkos::map_delete() { memoryKK->destroy_kokkos(k_sametag, sametag); sametag = nullptr; + max_same = 0; if (map_style == MAP_ARRAY) { memoryKK->destroy_kokkos(k_map_array, map_array); diff --git a/src/KOKKOS/region_block_kokkos.cpp b/src/KOKKOS/region_block_kokkos.cpp index 3802c22f968..b6ce93041c4 100644 --- a/src/KOKKOS/region_block_kokkos.cpp +++ b/src/KOKKOS/region_block_kokkos.cpp @@ -47,6 +47,7 @@ void RegBlockKokkos::match_all_kokkos(int groupbit_in, DAT::tdual_in { groupbit = groupbit_in; d_match = k_match_in.template view(); + k_remap.setup(domain); auto execution_space = ExecutionSpaceFromDevice::space; atomKK->sync(execution_space, X_MASK | MASK_MASK); d_x = atomKK->k_x.view(); diff --git a/src/KOKKOS/region_block_kokkos.h b/src/KOKKOS/region_block_kokkos.h index 24285e2448f..1fc6ee5a153 100644 --- a/src/KOKKOS/region_block_kokkos.h +++ b/src/KOKKOS/region_block_kokkos.h @@ -28,6 +28,7 @@ RegionStyle(block/kk/host,RegBlockKokkos); #include "kokkos_base.h" #include "kokkos_type.h" #include "math_special_kokkos.h" +#include "region_remap_kokkos.h" namespace LAMMPS_NS { @@ -48,6 +49,12 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { void match_all_kokkos(int, DAT::tdual_int_1d) override; + void prematch() override + { + RegBlock::prematch(); + k_remap.setup(domain); + } + // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION void operator()(TagRegBlockMatchAll, const int&) const; @@ -56,6 +63,9 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { KOKKOS_INLINE_FUNCTION int match_kokkos(double x, double y, double z) const { + // Region::match() maps the coordinate back into the box if periodic, since + // not all subclasses/methods treat a region extending beyond a periodic edge + k_remap.remap(x,y,z); if (dynamic) inverse_transform(x,y,z); if (openflag) return 1; return !(k_inside(x,y,z) ^ interior); @@ -69,6 +79,10 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { double xs, ys, zs; double xnear[3], xorig[3]; + // Region::surface() maps the coordinate back into the box if periodic, since + // not all subclasses/methods treat a region extending beyond a periodic edge + k_remap.remap(x, y, z); + if (dynamic) { xorig[0] = x; xorig[1] = y; xorig[2] = z; inverse_transform(x, y, z); @@ -82,9 +96,12 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { else ncontact = surface_exterior_kokkos(xnear, cutoff); } else { - // one of surface_int/ext() will return 0 - // so no need to worry about offset of contact indices - ncontact = surface_exterior_kokkos(xnear, cutoff) + surface_interior_kokkos(xnear, cutoff); + // most of the time, one of surface_int/ext() will return 0 + // however, when exactly on top of a periodic boundary + // both could return 1, so run exterior then interior + ncontact = surface_exterior_kokkos(xnear, cutoff); + if (ncontact == 0) + ncontact = surface_interior_kokkos(xnear, cutoff); } if (rotateflag && ncontact) { @@ -103,6 +120,7 @@ class RegBlockKokkos : public RegBlock, public KokkosBase { } Kokkos::View d_contact; + RegionRemapKokkos k_remap; private: int groupbit; diff --git a/src/KOKKOS/region_remap_kokkos.h b/src/KOKKOS/region_remap_kokkos.h new file mode 100644 index 00000000000..113295f5141 --- /dev/null +++ b/src/KOKKOS/region_remap_kokkos.h @@ -0,0 +1,117 @@ +/* -*- 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. +------------------------------------------------------------------------- */ + +#ifndef LMP_REGION_REMAP_KOKKOS_H +#define LMP_REGION_REMAP_KOKKOS_H + +#include "domain.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +/* ---------------------------------------------------------------------- + device-side copy of Domain::remap(double *) + + Region::match() and Region::surface() map a coordinate back into the + periodic box before testing it, because not every region subclass handles + a region that reaches past a periodic edge. The KOKKOS regions run the + same test inside a kernel, where the Domain instance is not reachable, so + the few box quantities remap() needs are copied into this struct on the + host and captured by value with the region object. +------------------------------------------------------------------------- */ + +struct RegionRemapKokkos { + int triclinic, xperiodic, yperiodic, zperiodic; + double lo[3], hi[3], period[3]; + double boxlo[3], h[6], h_inv[6]; + + RegionRemapKokkos() : triclinic(0), xperiodic(0), yperiodic(0), zperiodic(0) + { + for (int i = 0; i < 3; i++) lo[i] = hi[i] = period[i] = boxlo[i] = 0.0; + for (int i = 0; i < 6; i++) h[i] = h_inv[i] = 0.0; + } + + // refresh from the Domain instance; call whenever the box may have changed + + void setup(Domain *domain) + { + triclinic = domain->triclinic; + xperiodic = domain->xperiodic; + yperiodic = domain->yperiodic; + zperiodic = domain->zperiodic; + for (int i = 0; i < 3; i++) { + if (triclinic) { + lo[i] = domain->boxlo_lamda[i]; + hi[i] = domain->boxhi_lamda[i]; + period[i] = domain->prd_lamda[i]; + } else { + lo[i] = domain->boxlo[i]; + hi[i] = domain->boxhi[i]; + period[i] = domain->prd[i]; + } + boxlo[i] = domain->boxlo[i]; + } + for (int i = 0; i < 6; i++) { + h[i] = domain->h[i]; + h_inv[i] = domain->h_inv[i]; + } + } + +// NOLINTNEXTLINE + KOKKOS_INLINE_FUNCTION + void remap(double &x, double &y, double &z) const + { + double coord[3]; + + if (triclinic == 0) { + coord[0] = x; coord[1] = y; coord[2] = z; + } else { + const double d0 = x - boxlo[0]; + const double d1 = y - boxlo[1]; + const double d2 = z - boxlo[2]; + coord[0] = h_inv[0]*d0 + h_inv[5]*d1 + h_inv[4]*d2; + coord[1] = h_inv[1]*d1 + h_inv[3]*d2; + coord[2] = h_inv[2]*d2; + } + + if (xperiodic) { + while (coord[0] < lo[0]) coord[0] += period[0]; + while (coord[0] >= hi[0]) coord[0] -= period[0]; + if (coord[0] < lo[0]) coord[0] = lo[0]; + } + + if (yperiodic) { + while (coord[1] < lo[1]) coord[1] += period[1]; + while (coord[1] >= hi[1]) coord[1] -= period[1]; + if (coord[1] < lo[1]) coord[1] = lo[1]; + } + + if (zperiodic) { + while (coord[2] < lo[2]) coord[2] += period[2]; + while (coord[2] >= hi[2]) coord[2] -= period[2]; + if (coord[2] < lo[2]) coord[2] = lo[2]; + } + + if (triclinic == 0) { + x = coord[0]; y = coord[1]; z = coord[2]; + } else { + x = h[0]*coord[0] + h[5]*coord[1] + h[4]*coord[2] + boxlo[0]; + y = h[1]*coord[1] + h[3]*coord[2] + boxlo[1]; + z = h[2]*coord[2] + boxlo[2]; + } + } +}; + +} // namespace LAMMPS_NS + +#endif diff --git a/src/KOKKOS/region_sphere_kokkos.cpp b/src/KOKKOS/region_sphere_kokkos.cpp index b279891e1ba..f77e7e58c66 100644 --- a/src/KOKKOS/region_sphere_kokkos.cpp +++ b/src/KOKKOS/region_sphere_kokkos.cpp @@ -50,6 +50,7 @@ void RegSphereKokkos::match_all_kokkos(int groupbit_in, DAT::tdual_i { groupbit = groupbit_in; d_match = k_match_in.template view(); + k_remap.setup(domain); auto execution_space = ExecutionSpaceFromDevice::space; atomKK->sync(execution_space, X_MASK | MASK_MASK); d_x = atomKK->k_x.view(); diff --git a/src/KOKKOS/region_sphere_kokkos.h b/src/KOKKOS/region_sphere_kokkos.h index bde0a9f6193..395403a577c 100644 --- a/src/KOKKOS/region_sphere_kokkos.h +++ b/src/KOKKOS/region_sphere_kokkos.h @@ -27,6 +27,7 @@ RegionStyle(sphere/kk/host,RegSphereKokkos); #include "kokkos_base.h" #include "kokkos_type.h" +#include "region_remap_kokkos.h" namespace LAMMPS_NS { @@ -45,6 +46,12 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { void match_all_kokkos(int, DAT::tdual_int_1d) override; + void prematch() override + { + RegSphere::prematch(); + k_remap.setup(domain); + } + // NOLINTNEXTLINE KOKKOS_INLINE_FUNCTION void operator()(TagRegSphereMatchAll, const int&) const; @@ -53,6 +60,9 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { KOKKOS_INLINE_FUNCTION int match_kokkos(double x, double y, double z) const { + // Region::match() maps the coordinate back into the box if periodic, since + // not all subclasses/methods treat a region extending beyond a periodic edge + k_remap.remap(x,y,z); if (dynamic) inverse_transform(x,y,z); if (openflag) return 1; return !(k_inside(x,y,z) ^ interior); @@ -66,9 +76,14 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { double xs, ys, zs; double xnear[3], xorig[3]; - xorig[0] = x; xorig[1] = y; xorig[2] = z; - if (dynamic) + // Region::surface() maps the coordinate back into the box if periodic, since + // not all subclasses/methods treat a region extending beyond a periodic edge + k_remap.remap(x, y, z); + + if (dynamic) { + xorig[0] = x; xorig[1] = y; xorig[2] = z; inverse_transform(x, y, z); + } xnear[0] = x; xnear[1] = y; xnear[2] = z; @@ -77,9 +92,12 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { else ncontact = surface_exterior_kokkos(xnear, cutoff); } else { - // one of surface_int/ext() will return 0 - // so no need to worry about offset of contact indices - ncontact = surface_exterior_kokkos(xnear, cutoff) + surface_interior_kokkos(xnear, cutoff); + // most of the time, one of surface_int/ext() will return 0 + // however, when exactly on top of a periodic boundary + // both could return 1, so run exterior then interior + ncontact = surface_exterior_kokkos(xnear, cutoff); + if (ncontact == 0) + ncontact = surface_interior_kokkos(xnear, cutoff); } if (rotateflag && ncontact) { @@ -98,6 +116,7 @@ class RegSphereKokkos : public RegSphere, public KokkosBase { } Kokkos::View d_contact; + RegionRemapKokkos k_remap; private: int groupbit; From a06c0e2c3da957b12eef0aff3fbfb38c1e72771a Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 29 Aug 2026 17:53:50 +0000 Subject: [PATCH 39/43] KOKKOS: claim the per-atom buck energy and retire stale topology claims Two more host/device coherence defects, both found by building with -D KOKKOS_DEBUG_SYNC=on and running the unit tests, where the host side gets its own allocation and a missing sync() or modify() fails on a CPU instead of only corrupting results on a GPU. pair buck/coul/long/kk computed eatom and vatom on the device but never claimed or copied them back, unlike every other pair style on the same template, so a per-atom energy or virial read on the host was whatever the buffer held before. MolPairStyle:buck_coul_long failed the compute-sum check by 369 with a tolerance of 2e-12. NeighBondKokkos rebuilds the bond, angle, dihedral and improper lists from the atom topology on the device and then claims the device side. bond quartic edits neighbor->bondlist in place on the host and claims it, so the next rebuild collided with that claim and BondStyle:quartic aborted. Retire the claim before the rebuild, which overwrites every entry it uses anyway. --- src/KOKKOS/neigh_bond_kokkos.cpp | 40 +++++++++++++++++++++++ src/KOKKOS/pair_buck_coul_long_kokkos.cpp | 10 ++++++ 2 files changed, 50 insertions(+) diff --git a/src/KOKKOS/neigh_bond_kokkos.cpp b/src/KOKKOS/neigh_bond_kokkos.cpp index 80ae98fc47c..f083afdcb56 100644 --- a/src/KOKKOS/neigh_bond_kokkos.cpp +++ b/src/KOKKOS/neigh_bond_kokkos.cpp @@ -255,6 +255,11 @@ template void NeighBondKokkos::bond_all() { atomKK->sync(execution_space, BOND_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_bondlist.clear_sync_state(); v_bondlist = k_bondlist.view(); num_bond = atomKK->k_num_bond.view(); bond_atom = atomKK->k_bond_atom.view(); @@ -338,6 +343,11 @@ template void NeighBondKokkos::bond_partial() { atomKK->sync(execution_space, BOND_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_bondlist.clear_sync_state(); v_bondlist = k_bondlist.view(); num_bond = atomKK->k_num_bond.view(); bond_atom = atomKK->k_bond_atom.view(); @@ -447,6 +457,11 @@ template void NeighBondKokkos::angle_all() { atomKK->sync(execution_space, ANGLE_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_anglelist.clear_sync_state(); v_anglelist = k_anglelist.view(); num_angle = atomKK->k_num_angle.view(); angle_atom1 = atomKK->k_angle_atom1.view(); @@ -536,6 +551,11 @@ template void NeighBondKokkos::angle_partial() { atomKK->sync(execution_space, ANGLE_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_anglelist.clear_sync_state(); v_anglelist = k_anglelist.view(); num_angle = atomKK->k_num_angle.view(); angle_atom1 = atomKK->k_angle_atom1.view(); @@ -666,6 +686,11 @@ template void NeighBondKokkos::dihedral_all() { atomKK->sync(execution_space, DIHEDRAL_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_dihedrallist.clear_sync_state(); v_dihedrallist = k_dihedrallist.view(); num_dihedral = atomKK->k_num_dihedral.view(); dihedral_atom1 = atomKK->k_dihedral_atom1.view(); @@ -760,6 +785,11 @@ template void NeighBondKokkos::dihedral_partial() { atomKK->sync(execution_space, DIHEDRAL_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_dihedrallist.clear_sync_state(); v_dihedrallist = k_dihedrallist.view(); num_dihedral = atomKK->k_num_dihedral.view(); dihedral_atom1 = atomKK->k_dihedral_atom1.view(); @@ -913,6 +943,11 @@ template void NeighBondKokkos::improper_all() { atomKK->sync(execution_space, IMPROPER_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_improperlist.clear_sync_state(); v_improperlist = k_improperlist.view(); num_improper = atomKK->k_num_improper.view(); improper_atom1 = atomKK->k_improper_atom1.view(); @@ -1007,6 +1042,11 @@ template void NeighBondKokkos::improper_partial() { atomKK->sync(execution_space, IMPROPER_MASK); + // the loop below rebuilds the whole list from the atom topology, so retire any + // outstanding claim first: a host style that edits the list in place (bond + // quartic breaks bonds there) leaves one behind, and the claim taken at the + // end of this function would then collide with it + k_improperlist.clear_sync_state(); v_improperlist = k_improperlist.view(); num_improper = atomKK->k_num_improper.view(); improper_atom1 = atomKK->k_improper_atom1.view(); diff --git a/src/KOKKOS/pair_buck_coul_long_kokkos.cpp b/src/KOKKOS/pair_buck_coul_long_kokkos.cpp index 4ee2bb0526b..927335c9eea 100644 --- a/src/KOKKOS/pair_buck_coul_long_kokkos.cpp +++ b/src/KOKKOS/pair_buck_coul_long_kokkos.cpp @@ -142,6 +142,16 @@ void PairBuckCoulLongKokkos::compute(int eflag_in, int vflag_in) virial[5] += static_cast(ev.v[5]); } + if (eflag_atom) { + k_eatom.template modify(); + k_eatom.sync_host(); + } + + if (vflag_atom) { + k_vatom.template modify(); + k_vatom.sync_host(); + } + if (vflag_fdotr) pair_virial_fdotr_compute(this); copymode = 0; From e9c3eda55e08a4687e42f03197dd5c3011796952 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Sat, 29 Aug 2026 19:37:33 +0000 Subject: [PATCH 40/43] KOKKOS: correct energy accounting in fix addforce/kk, efield/kk and shake/kk Four defects found while reviewing the preceding commits. - fix addforce/kk enumerated the EQUAL and ATOM cases separately when accumulating the energy of a variable-controlled force, and so dropped every component given as a plain number. With 'fix ID group addforce 1.0 v_fy 0.5' the reported energy was 0.15 instead of -30.54. FixAddForce::post_force() tests 'if (xstyle)', which covers the constant case as well; use the same form, and share the per-component values with the virial, which already had it right but was therefore inconsistent with the energy next to it. - fix efield/kk sized and synchronized its efield array only for varflag == ATOM, but an atom-style energy or potential variable also writes into efield[i][3]. The device copy was then never refreshed from the host, and the array was never grown for that case. Include estyle and pstyle in both conditions, and in the sizing condition of FixEfield::post_force(), which has the same omission (cf. the corresponding condition in FixAddForce::post_force()). - fix shake/kk reset the per-atom virial in post_force() but not the per-atom energy. That was harmless while nothing filled the latter, but now that min_post_force() does, a run following a minimization kept reporting the per-atom energies of the last minimization step. Follow FixShake::post_force() and use ev_init(), reallocating the per-atom energy dual view alongside the virial one. - the SHAKE bond statistics of min_post_force() counted two atoms per restraint whether or not this processor owns them, unlike FixShake::min_post_force(), which over-reports the bond count and distorts the average for a cluster reaching across a boundary. Verified against the non-accelerated style with 1, 2 and 4 OpenMP threads. The first and third defect are reproduced by the inputs above and are gone; the second only shows up when host and device memory are distinct, and the fourth only when a cluster spans a boundary, so neither could be exercised on this host-only single-rank build. --- src/KOKKOS/fix_addforce_kokkos.cpp | 19 ++++++++++--------- src/KOKKOS/fix_efield_kokkos.cpp | 10 +++++++--- src/KOKKOS/fix_shake_kokkos.cpp | 30 ++++++++++++++++++++---------- src/fix.h | 11 +++++++---- src/fix_efield.cpp | 2 +- 5 files changed, 45 insertions(+), 27 deletions(-) diff --git a/src/KOKKOS/fix_addforce_kokkos.cpp b/src/KOKKOS/fix_addforce_kokkos.cpp index 902b6e5cfcb..09b4008562f 100644 --- a/src/KOKKOS/fix_addforce_kokkos.cpp +++ b/src/KOKKOS/fix_addforce_kokkos.cpp @@ -241,15 +241,19 @@ void FixAddForceKokkos::operator()(TagFixAddForceNonConstant, const x_i[2] = static_cast(x(i,2)); auto unwrapKK = DomainKokkos::unmap(prd,h,triclinic,x_i,image(i)); + // an atom-style variable supplies the value per atom, any other style + // (constant or equal-style variable) the same value for all of them + + const double xv = (xstyle == ATOM) ? static_cast(d_sforce(i,0)) : xvalue; + const double yv = (ystyle == ATOM) ? static_cast(d_sforce(i,1)) : yvalue; + const double zv = (zstyle == ATOM) ? static_cast(d_sforce(i,2)) : zvalue; + if (estyle == ATOM) { result[0] += static_cast(d_sforce(i,3)); } else { - if (xstyle == EQUAL) result[0] -= xvalue * unwrapKK[0]; - if (ystyle == EQUAL) result[0] -= yvalue * unwrapKK[1]; - if (zstyle == EQUAL) result[0] -= zvalue * unwrapKK[2]; - if (xstyle == ATOM) result[0] -= static_cast(d_sforce(i,0)) * unwrapKK[0]; - if (ystyle == ATOM) result[0] -= static_cast(d_sforce(i,1)) * unwrapKK[1]; - if (zstyle == ATOM) result[0] -= static_cast(d_sforce(i,2)) * unwrapKK[2]; + if (xstyle) result[0] -= xv * unwrapKK[0]; + if (ystyle) result[0] -= yv * unwrapKK[1]; + if (zstyle) result[0] -= zv * unwrapKK[2]; } result[1] += static_cast(f(i,0)); result[2] += static_cast(f(i,1)); @@ -262,9 +266,6 @@ void FixAddForceKokkos::operator()(TagFixAddForceNonConstant, const else if (zstyle) f(i,2) += static_cast(zvalue_kk); if (evflag) { - const double xv = (xstyle == ATOM) ? static_cast(d_sforce(i,0)) : xvalue; - const double yv = (ystyle == ATOM) ? static_cast(d_sforce(i,1)) : yvalue; - const double zv = (zstyle == ATOM) ? static_cast(d_sforce(i,2)) : zvalue; KK_FLOAT v[6]; v[0] = xstyle ? static_cast(xv * unwrapKK[0]) : static_cast(0.0); v[1] = ystyle ? static_cast(yv * unwrapKK[1]) : static_cast(0.0); diff --git a/src/KOKKOS/fix_efield_kokkos.cpp b/src/KOKKOS/fix_efield_kokkos.cpp index 1fcf88a4769..75a572480cd 100644 --- a/src/KOKKOS/fix_efield_kokkos.cpp +++ b/src/KOKKOS/fix_efield_kokkos.cpp @@ -119,9 +119,11 @@ void FixEfieldKokkos::post_force(int vflag) d_match = k_match.template view(); } - // reallocate sforce array if necessary + // reallocate efield array if necessary + // an atom-style energy or potential variable writes into efield[i][3], + // so the array is needed for those as well, not only for varflag == ATOM - if (varflag == ATOM && atom->nmax > maxatom) { + if (((varflag == ATOM) || (estyle == ATOM) || (pstyle == ATOM)) && atom->nmax > maxatom) { maxatom = atom->nmax; memoryKK->destroy_kokkos(k_efield,efield); memoryKK->create_kokkos(k_efield,efield,maxatom,4,"efield:efield"); @@ -157,7 +159,9 @@ void FixEfieldKokkos::post_force(int vflag) FixEfield::update_efield_variables(); - if (varflag == ATOM) { // this can be removed when variable class is ported to Kokkos + // this can be removed when the variable class is ported to Kokkos + + if ((varflag == ATOM) || (estyle == ATOM) || (pstyle == ATOM)) { k_efield.modify_host(); k_efield.sync(); } diff --git a/src/KOKKOS/fix_shake_kokkos.cpp b/src/KOKKOS/fix_shake_kokkos.cpp index c3174f13c3f..838e172770a 100644 --- a/src/KOKKOS/fix_shake_kokkos.cpp +++ b/src/KOKKOS/fix_shake_kokkos.cpp @@ -563,10 +563,12 @@ void FixShakeKokkos::operator()(TagFixShakeMinPostForce(ev,count,atomlist,total,eb,v); } if (output_every && !is_angle) { - Kokkos::atomic_add(&d_b_stats(type_idx, 0), 1.0); - Kokkos::atomic_add(&d_b_stats(type_idx, 1), (double)r); - Kokkos::atomic_add(&d_b_stats(type_idx, 0), 1.0); - Kokkos::atomic_add(&d_b_stats(type_idx, 1), (double)r); + // only atoms owned by this processor are counted, as on the CPU + const double nown = (double) ((idx0 < nlocal) + (idx1 < nlocal)); + if (nown > 0.0) { + Kokkos::atomic_add(&d_b_stats(type_idx, 0), nown); + Kokkos::atomic_add(&d_b_stats(type_idx, 1), nown * (double)r); + } Kokkos::atomic_max(&d_b_stats(type_idx, 2), (double)r); Kokkos::atomic_min(&d_b_stats(type_idx, 3), (double)r); } @@ -656,15 +658,23 @@ void FixShakeKokkos::post_force(int vflag) comm->forward_comm(this); k_xshake.sync(); - // virial setup - - // the per-atom virial is accumulated into a dual view, so the plain - // base-class vatom array must not be allocated here (alloc = 0) + // energy and virial setup, as FixShake::post_force() does. the per-atom + // energy and virial are accumulated into dual views, so the plain + // base-class arrays must not be allocated here (alloc = 0) - v_init(vflag,0); + int eflag = eflag_pre_reverse; + ev_init(eflag,vflag,0); - // reallocate the per-atom virial dual view if necessary + // reallocate the per-atom energy and virial dual views if necessary. the + // constraint forces contribute no per-atom energy during dynamics, but the + // freshly created view keeps it zeroed rather than holding on to whatever + // a preceding minimization left behind + if (eflag_atom) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"shake:eatom"); + d_eatom = k_eatom.template view(); + } if (vflag_atom) { memoryKK->destroy_kokkos(k_vatom,vatom); memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"shake:vatom"); diff --git a/src/fix.h b/src/fix.h index 184b5409d21..ab4d6897b96 100644 --- a/src/fix.h +++ b/src/fix.h @@ -283,10 +283,13 @@ class Fix : protected Pointers { int dynamic; // recount atoms for temperature computes // alloc = 0 tells ev_setup()/v_setup() to only update the flags and the - // maxeatom/maxvatom sizes, but not to allocate or zero the plain per-atom - // arrays. Styles that manage eatom/vatom themselves (e.g. the KOKKOS - // variants, which store them in dual views) must pass alloc = 0, since - // otherwise the plain arrays allocated here are orphaned by the style. + // maxeatom/maxvatom/maxcvatom sizes, but not to allocate or zero any of the + // plain per-atom arrays. Styles that manage eatom/vatom themselves (e.g. + // the KOKKOS variants, which store them in dual views) must pass alloc = 0, + // since otherwise the plain arrays allocated here are orphaned by the style. + // Note that this covers cvatom as well: a style that passes alloc = 0 has to + // provide the centroid virial itself or set centroidstressflag to + // CENTROID_NOTAVAIL, since cvatom is then left unallocated. void ev_init(int eflag, int vflag, int alloc = 1) { diff --git a/src/fix_efield.cpp b/src/fix_efield.cpp index 24a5724cd69..655bf46949b 100644 --- a/src/fix_efield.cpp +++ b/src/fix_efield.cpp @@ -305,7 +305,7 @@ void FixEfield::post_force(int vflag) // reallocate efield array if necessary - if ((varflag == ATOM) && (atom->nmax > maxatom)) { + if (((varflag == ATOM) || (estyle == ATOM) || (pstyle == ATOM)) && (atom->nmax > maxatom)) { maxatom = atom->nmax; memory->destroy(efield); memory->create(efield, maxatom, 4, "efield:efield"); From c34cb6e939a5612132d74ce82f6fe449f75fd4cb Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 31 Aug 2026 15:02:24 +0000 Subject: [PATCH 41/43] STUBS: add MPI_STATUSES_IGNORE Since 36aa1bebff the tiled grid communication of grid3d_kokkos.cpp waits on its receives with MPI_Waitall() and MPI_STATUSES_IGNORE, which the MPI stubs library does not define, so a serial build of the KOKKOS package with KSPACE fails to compile. Define it alongside MPI_STATUS_IGNORE rather than working around it in the caller. --- src/STUBS/mpi.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/STUBS/mpi.h b/src/STUBS/mpi.h index 03b8d1fc263..d9a92dc778f 100644 --- a/src/STUBS/mpi.h +++ b/src/STUBS/mpi.h @@ -58,6 +58,7 @@ #define MPI_ANY_SOURCE (-1) #define MPI_STATUS_IGNORE NULL +#define MPI_STATUSES_IGNORE NULL #define MPI_Comm int #define MPI_Request int From d636174b2a660a4915605fa9f2c613c1999acdbc Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Mon, 31 Aug 2026 15:30:50 +0000 Subject: [PATCH 42/43] KOKKOS: size neighbor bins from actual occupancy to avoid setup blowup NBinKokkos::bin_atoms() grew atoms_per_bin by a fixed increment of 16 and re-binned every atom once per increment whenever a bin overflowed. For skewed atom distributions - such as the large-cutoff bins the KOKKOS package uses on GPUs, especially with pair_style hybrid - the busiest bin can hold thousands of atoms, so this loop performs O(max bin occupancy) reallocations and full re-bins. That makes the first neighbor build during run setup take minutes and grow GPU memory monotonically (see issue #4988). Because the atomic increment in binatomsItem() runs for every atom even after a bin is full, bincount holds the true occupancy of every bin after a pass. Use a parallel reduction over bincount to size atoms_per_bin from the actual maximum (plus a small margin) in a single step, so the loop converges in at most two passes regardless of how skewed the distribution is. Results are unchanged; only the setup cost is fixed. --- src/KOKKOS/nbin_kokkos.cpp | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/src/KOKKOS/nbin_kokkos.cpp b/src/KOKKOS/nbin_kokkos.cpp index 6d214bb51ba..3652b0bf338 100644 --- a/src/KOKKOS/nbin_kokkos.cpp +++ b/src/KOKKOS/nbin_kokkos.cpp @@ -122,7 +122,32 @@ void NBinKokkos::bin_atoms() Kokkos::deep_copy(h_resize, d_resize); if (h_resize()) { - atoms_per_bin += 16; + // A bin overflowed its capacity. bincount now holds the true + // occupancy of every bin, because the atomic increment in + // binatomsItem() runs for every binned atom whether or not the bin + // was already full. Size atoms_per_bin from the actual maximum + // occupancy in a single step instead of growing by a fixed increment + // and re-binning all atoms once per increment. The latter costs + // O(max bin occupancy) reallocations and re-bins and dominates + // neighbor setup for skewed distributions, e.g. the large cutoff bins + // used by the KOKKOS package on GPUs. + + auto d_bincount = k_bincount.view(); + int max_bincount = 0; + Kokkos::parallel_reduce(mbins, + LAMMPS_LAMBDA(const int i, int &max_val) { + max_val = MAX(max_val,d_bincount[i]); + },Kokkos::Max(max_bincount)); + + // grow to the true maximum occupancy plus ~10% headroom (at least 16) + // so small density fluctuations on later steps do not immediately + // force another regrow. Reaching this branch means a bin overflowed, + // so max_bincount > atoms_per_bin and the new capacity strictly + // exceeds both the old one and the true occupancy: the next pass + // cannot overflow, bounding the loop at one more re-bin. + + atoms_per_bin = max_bincount + MAX(16,max_bincount/10); + k_bins = DAT::tdual_int_2d("Neighbor::bins", mbins, atoms_per_bin); bins = k_bins.view(); c_bins = bins; From fc2cf8c8a3b435342a171e3382e65abf4e944079 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Wed, 26 Aug 2026 04:14:28 +0000 Subject: [PATCH 43/43] KOKKOS PACE: fall back to global scratch when shared memory is exceeded The ComputeNeigh kernel in pair pace/kk and pace/extrapolation/kk caches the short neighbor list in level-0 (on-chip shared) team scratch memory, sized team_size*maxneigh*sizeof(int). On GPUs shared memory is a very limited resource, so runs with many neighbors and/or many atomic species could abort with "Requested too much scratch memory on level 0" (CUDA) or "could not find a valid team size" (HIP). See lammps/lammps#5063. Query the maximum available level-0 scratch from Kokkos via TeamPolicy::scratch_size_max(0) and transparently fall back to level-1 (global memory) scratch when the request does not fit, printing a warning the first time. The limit is queried rather than hard-coded (e.g. 48 KiB), so larger shared-memory limits such as the opt-in >48 KiB shared memory in newer Kokkos are picked up automatically. Add a "neigh" pair_style keyword (auto|shared|global) so the user can override the automatic choice, and document it. Checked on a Kokkos Serial build of the fcc-Cu deck: all three keyword values run and give forces identical to the non-accelerated evaluator, an unknown value is rejected by the KOKKOS styles, and the non-accelerated styles stop with "Unknown pair_style pace keyword: neigh" as the docs now state. The scratch-level selection itself only takes effect on a GPU backend and needs a CUDA or HIP build to exercise. --- doc/src/pair_pace.rst | 30 ++++++- src/KOKKOS/pair_pace_extrapolation_kokkos.cpp | 88 ++++++++++++++++++- src/KOKKOS/pair_pace_extrapolation_kokkos.h | 14 +++ src/KOKKOS/pair_pace_kokkos.cpp | 88 ++++++++++++++++++- src/KOKKOS/pair_pace_kokkos.h | 14 +++ 5 files changed, 229 insertions(+), 5 deletions(-) diff --git a/doc/src/pair_pace.rst b/doc/src/pair_pace.rst index 04060afbb2d..90aa8c39461 100644 --- a/doc/src/pair_pace.rst +++ b/doc/src/pair_pace.rst @@ -24,10 +24,14 @@ Syntax .. parsed-literal:: - keyword = *product* or *recursive* or *chunksize* + keyword = *product* or *recursive* or *chunksize* or *neigh* *product* = use product algorithm for basis functions *recursive* = use recursive algorithm for basis functions *chunksize* value = number of atoms in each pass + *neigh* value = *auto* or *shared* or *global* + *auto* = automatically select team scratch memory level for the neighbor list build (default) + *shared* = force on-chip (level 0) shared memory scratch + *global* = force global (level 1) memory scratch .. code-block:: LAMMPS @@ -90,6 +94,30 @@ For example if there are 8192 atoms in the simulation and the *chunksize* is set to 4096, the ACE calculation will be broken up into two passes (running on a single GPU). +.. versionadded:: TBD + +The keyword *neigh* is only recognized by the KOKKOS versions of the pair +styles (*pace/kk* and *pace/extrapolation/kk*); the non-accelerated styles +stop with an "unknown keyword" error if it is given. It only has an effect +on GPU backends. It controls which level of Kokkos team scratch memory is +used to build the short neighbor list. Level 0 is fast on-chip shared +memory, but it is a limited resource that can be exceeded when atoms +have many neighbors and/or when there are many atomic species, which +would otherwise abort the run with an error such as "Requested too much +scratch memory on level 0". Level 1 is (much larger) global memory, +which avoids the limit at the cost of slower access. + +With the default value *auto*, the pair style queries the amount of +shared memory available on the device (rather than assuming a fixed +value such as 48 KiB, so that larger limits available in newer versions +of the Kokkos library are used automatically) and transparently falls +back to level 1 when the request does not fit into level 0, printing a +warning the first time this happens. The value *shared* forces the use +of level 0 (on-chip) scratch memory, and *global* forces the use of +level 1 (global) scratch memory; the latter can be used to silence the +fallback warning or to force global memory when the automatic heuristic +is too conservative. + Extrapolation grade """"""""""""""""""" diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp index 21d5d786e46..9f707c25d54 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp @@ -28,6 +28,7 @@ #include "memory_kokkos.h" #include "neighbor_kokkos.h" #include "neigh_request.h" +#include "utils.h" #include "ace-evaluator/ace_version.h" #include "ace-evaluator/ace_radial.h" @@ -64,6 +65,10 @@ PairPACEExtrapolationKokkos::PairPACEExtrapolationKokkos(LAMMPS *lmp datamask_modify = EMPTY_MASK; host_flag = (execution_space == HostKK); + + neigh_scratch_request = NEIGH_SCRATCH_AUTO; + neigh_scratch_level = 0; + neigh_scratch_warned = 0; } /* ---------------------------------------------------------------------- @@ -501,6 +506,76 @@ double PairPACEExtrapolationKokkos::init_one(int i, int j) return cutone; } +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +template +void PairPACEExtrapolationKokkos::settings(int narg, char **arg) +{ + // intercept the KOKKOS-only "neigh" keyword, which selects the team scratch + // memory level used to build the short neighbor list, then forward the + // remaining keywords to the CPU base class + + auto base_arg = new char*[narg]; + int base_narg = 0; + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg], "neigh") == 0) { + if (iarg+2 > narg) + utils::missing_cmd_args(FLERR, "pair_style pace/extrapolation neigh", error); + if (strcmp(arg[iarg+1], "auto") == 0) + neigh_scratch_request = NEIGH_SCRATCH_AUTO; + else if (strcmp(arg[iarg+1], "shared") == 0) + neigh_scratch_request = NEIGH_SCRATCH_SHARED; + else if (strcmp(arg[iarg+1], "global") == 0) + neigh_scratch_request = NEIGH_SCRATCH_GLOBAL; + else + error->all(FLERR, "Unknown pair_style pace/extrapolation neigh keyword: {}", arg[iarg+1]); + iarg += 2; + } else { + base_arg[base_narg++] = arg[iarg]; + iarg++; + } + } + + PairPACEExtrapolation::settings(base_narg, base_arg); + + delete[] base_arg; +} + +/* ---------------------------------------------------------------------- + select the team scratch memory level for the ComputeNeigh short neighbor + list build; falls back from level 0 (fast on-chip shared memory) to level 1 + (global memory) when the request does not fit into the available shared + memory, unless the user forced a level +------------------------------------------------------------------------- */ + +template +int PairPACEExtrapolationKokkos::neigh_scratch_level_select(int scratch_size, int max_level0) +{ + // honor an explicit user request + if (neigh_scratch_request == NEIGH_SCRATCH_SHARED) return 0; + if (neigh_scratch_request == NEIGH_SCRATCH_GLOBAL) return 1; + + // automatic: use fast level-0 (shared) scratch when it fits, otherwise fall + // back to level-1 (global) scratch. max_level0 is queried from Kokkos rather + // than hard-coded, so larger shared-memory limits (e.g. the opt-in >48 KiB + // shared memory available in newer Kokkos) are used automatically. + if (scratch_size <= max_level0) return 0; + + if (!neigh_scratch_warned && comm->me == 0) { + error->warning(FLERR, + "Pair pace/extrapolation/kk short neighbor list needs {} bytes of team " + "scratch memory but only {} bytes of on-chip (level-0) shared memory are " + "available; falling back to slower global (level-1) memory. Reduce the " + "neighbor count or use the pair_style 'neigh global' keyword to silence " + "this warning.", scratch_size, max_level0); + neigh_scratch_warned = 1; + } + return 1; +} + /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ @@ -676,7 +751,16 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in check_team_size_for(chunk_size,team_size,vector_length); int scratch_size = scratch_size_helper(team_size * maxneigh); typename Kokkos::TeamPolicy policy_neigh(chunk_size,team_size,vector_length); - policy_neigh = policy_neigh.set_scratch_size(0, Kokkos::PerTeam(scratch_size)); + + // The ComputeNeigh kernel caches the short neighbor list in team scratch + // memory. On GPUs level-0 scratch is fast on-chip shared memory but is a + // scarce resource: with many neighbors and/or atom types the request can + // exceed what the device provides and abort the run (see + // https://github.com/lammps/lammps/issues/5063). Query the level-0 limit + // from Kokkos (never hard-coded) and transparently fall back to level-1 + // (global memory) scratch when the request does not fit. + neigh_scratch_level = neigh_scratch_level_select(scratch_size, policy_neigh.scratch_size_max(0)); + policy_neigh = policy_neigh.set_scratch_size(neigh_scratch_level, Kokkos::PerTeam(scratch_size)); Kokkos::parallel_for("ComputeNeigh",policy_neigh,*this); } @@ -835,7 +919,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig // If it is, inside is assigned to 1, otherwise -1 const int team_rank = team.team_rank(); const int scratch_shift = team_rank * maxneigh; // offset into pointer for entire team - int* inside = (int*)team.team_shmem().get_shmem(team.team_size() * maxneigh * sizeof(int), 0) + scratch_shift; + int* inside = (int*)team.team_shmem().get_shmem(team.team_size() * maxneigh * sizeof(int), neigh_scratch_level) + scratch_shift; // loop over list of all neighbors within force cutoff // distsq[] = distance sq to each diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.h b/src/KOKKOS/pair_pace_extrapolation_kokkos.h index 4aec48a7c3f..bc76c7fc0bb 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.h +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.h @@ -56,6 +56,7 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { ~PairPACEExtrapolationKokkos() override; void compute(int, int) override; + void settings(int, char **) override; void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; @@ -110,6 +111,19 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, total_num_functions_max, idx_sph_max; int host_flag; + // team scratch memory level used by the ComputeNeigh short neighbor list build: + // NEIGH_SCRATCH_AUTO - automatically use level 0 (fast on-chip shared + // memory) when it fits, else fall back to level 1 + // (global memory) + // NEIGH_SCRATCH_SHARED - always use level 0 + // NEIGH_SCRATCH_GLOBAL - always use level 1 + enum { NEIGH_SCRATCH_AUTO = 0, NEIGH_SCRATCH_SHARED, NEIGH_SCRATCH_GLOBAL }; + int neigh_scratch_request; // user preference (pair_style "neigh" keyword) + int neigh_scratch_level; // level actually used by ComputeNeigh (0 or 1) + int neigh_scratch_warned; // whether the auto-fallback warning was printed + + int neigh_scratch_level_select(int scratch_size, int max_level0); + int eflag, vflag; int neighflag, max_ndensity; diff --git a/src/KOKKOS/pair_pace_kokkos.cpp b/src/KOKKOS/pair_pace_kokkos.cpp index f954fd4b69c..e72af91da85 100644 --- a/src/KOKKOS/pair_pace_kokkos.cpp +++ b/src/KOKKOS/pair_pace_kokkos.cpp @@ -28,6 +28,7 @@ #include "memory_kokkos.h" #include "neighbor_kokkos.h" #include "neigh_request.h" +#include "utils.h" #include "ace-evaluator/ace_version.h" #include "ace-evaluator/ace_radial.h" @@ -65,6 +66,10 @@ PairPACEKokkos::PairPACEKokkos(LAMMPS *lmp) : PairPACE(lmp) datamask_modify = EMPTY_MASK; host_flag = (execution_space == HostKK); + + neigh_scratch_request = NEIGH_SCRATCH_AUTO; + neigh_scratch_level = 0; + neigh_scratch_warned = 0; } /* ---------------------------------------------------------------------- @@ -473,6 +478,76 @@ double PairPACEKokkos::init_one(int i, int j) return cutone; } +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +template +void PairPACEKokkos::settings(int narg, char **arg) +{ + // intercept the KOKKOS-only "neigh" keyword, which selects the team scratch + // memory level used to build the short neighbor list, then forward the + // remaining keywords to the CPU base class + + auto base_arg = new char*[narg]; + int base_narg = 0; + int iarg = 0; + while (iarg < narg) { + if (strcmp(arg[iarg], "neigh") == 0) { + if (iarg+2 > narg) + utils::missing_cmd_args(FLERR, "pair_style pace neigh", error); + if (strcmp(arg[iarg+1], "auto") == 0) + neigh_scratch_request = NEIGH_SCRATCH_AUTO; + else if (strcmp(arg[iarg+1], "shared") == 0) + neigh_scratch_request = NEIGH_SCRATCH_SHARED; + else if (strcmp(arg[iarg+1], "global") == 0) + neigh_scratch_request = NEIGH_SCRATCH_GLOBAL; + else + error->all(FLERR, "Unknown pair_style pace neigh keyword: {}", arg[iarg+1]); + iarg += 2; + } else { + base_arg[base_narg++] = arg[iarg]; + iarg++; + } + } + + PairPACE::settings(base_narg, base_arg); + + delete[] base_arg; +} + +/* ---------------------------------------------------------------------- + select the team scratch memory level for the ComputeNeigh short neighbor + list build; falls back from level 0 (fast on-chip shared memory) to level 1 + (global memory) when the request does not fit into the available shared + memory, unless the user forced a level +------------------------------------------------------------------------- */ + +template +int PairPACEKokkos::neigh_scratch_level_select(int scratch_size, int max_level0) +{ + // honor an explicit user request + if (neigh_scratch_request == NEIGH_SCRATCH_SHARED) return 0; + if (neigh_scratch_request == NEIGH_SCRATCH_GLOBAL) return 1; + + // automatic: use fast level-0 (shared) scratch when it fits, otherwise fall + // back to level-1 (global) scratch. max_level0 is queried from Kokkos rather + // than hard-coded, so larger shared-memory limits (e.g. the opt-in >48 KiB + // shared memory available in newer Kokkos) are used automatically. + if (scratch_size <= max_level0) return 0; + + if (!neigh_scratch_warned && comm->me == 0) { + error->warning(FLERR, + "Pair pace/kk short neighbor list needs {} bytes of team scratch memory " + "but only {} bytes of on-chip (level-0) shared memory are available; " + "falling back to slower global (level-1) memory. Reduce the neighbor " + "count or use the pair_style 'neigh global' keyword to silence this " + "warning.", scratch_size, max_level0); + neigh_scratch_warned = 1; + } + return 1; +} + /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ @@ -636,7 +711,16 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) check_team_size_for(chunk_size,team_size,vector_length); int scratch_size = scratch_size_helper(team_size * maxneigh); typename Kokkos::TeamPolicy policy_neigh(chunk_size,team_size,vector_length); - policy_neigh = policy_neigh.set_scratch_size(0, Kokkos::PerTeam(scratch_size)); + + // The ComputeNeigh kernel caches the short neighbor list in team scratch + // memory. On GPUs level-0 scratch is fast on-chip shared memory but is a + // scarce resource: with many neighbors and/or atom types the request can + // exceed what the device provides and abort the run (see + // https://github.com/lammps/lammps/issues/5063). Query the level-0 limit + // from Kokkos (never hard-coded) and transparently fall back to level-1 + // (global memory) scratch when the request does not fit. + neigh_scratch_level = neigh_scratch_level_select(scratch_size, policy_neigh.scratch_size_max(0)); + policy_neigh = policy_neigh.set_scratch_size(neigh_scratch_level, Kokkos::PerTeam(scratch_size)); Kokkos::parallel_for("ComputeNeigh",policy_neigh,*this); } @@ -781,7 +865,7 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen // If it is, inside is assigned to 1, otherwise -1 const int team_rank = team.team_rank(); const int scratch_shift = team_rank * maxneigh; // offset into pointer for entire team - int* inside = (int*)team.team_shmem().get_shmem(team.team_size() * maxneigh * sizeof(int), 0) + scratch_shift; + int* inside = (int*)team.team_shmem().get_shmem(team.team_size() * maxneigh * sizeof(int), neigh_scratch_level) + scratch_shift; // loop over list of all neighbors within force cutoff // distsq[] = distance sq to each diff --git a/src/KOKKOS/pair_pace_kokkos.h b/src/KOKKOS/pair_pace_kokkos.h index cd3cd6b037f..f1203b8386c 100644 --- a/src/KOKKOS/pair_pace_kokkos.h +++ b/src/KOKKOS/pair_pace_kokkos.h @@ -55,6 +55,7 @@ class PairPACEKokkos : public PairPACE { ~PairPACEKokkos() override; void compute(int, int) override; + void settings(int, char **) override; void coeff(int, char **) override; void init_style() override; double init_one(int, int) override; @@ -105,6 +106,19 @@ class PairPACEKokkos : public PairPACE { int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, idx_sph_max; int host_flag; + // team scratch memory level used by the ComputeNeigh short neighbor list build: + // NEIGH_SCRATCH_AUTO - automatically use level 0 (fast on-chip shared + // memory) when it fits, else fall back to level 1 + // (global memory) + // NEIGH_SCRATCH_SHARED - always use level 0 + // NEIGH_SCRATCH_GLOBAL - always use level 1 + enum { NEIGH_SCRATCH_AUTO = 0, NEIGH_SCRATCH_SHARED, NEIGH_SCRATCH_GLOBAL }; + int neigh_scratch_request; // user preference (pair_style "neigh" keyword) + int neigh_scratch_level; // level actually used by ComputeNeigh (0 or 1) + int neigh_scratch_warned; // whether the auto-fallback warning was printed + + int neigh_scratch_level_select(int scratch_size, int max_level0); + int eflag, vflag; int neighflag, max_ndensity;