KOKKOS: fix host/device coherence bugs and include-group binning - #50
Merged
Merged
Conversation
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<Kokkos::Serial>::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 65b2d13)
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 f8d10ee)
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 d50b393)
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 0edfa53)
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 a55c9cd)
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<Kokkos::Serial>::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 0dab7fa)
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 7303197)
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 7660543)
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<Kokkos::Serial>::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 a3a56dc)
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 52185ff)
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 cb16ae1)
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 26c0392)
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 fcbdc04)
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 3299c05)
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<Kokkos::Serial>::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 c98911b)
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
…hake 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.
Brings in the fixes for the orphaned base-class eatom/vatom arrays in the KOKKOS fixes and pair styles, the per-atom energy and virial of fix shake/kk, the energy accounting of fix addforce/kk and fix efield/kk, and the shared region contact list. Four conflicts, all in code both sides had touched: - region_block_kokkos.h, region_sphere_kokkos.h: this branch's openflag handling is kept, since running surface_exterior() first and surface_interior() only if it found nothing is the correct fix for a coordinate sitting exactly on a periodic boundary, where the other side's sum would have counted two contacts while both wrote index 0. It now takes the caller-owned contact list of the other side. The shared d_contact view is gone with that change, the coordinate remap member added here stays. - fix_addforce_kokkos.cpp, fix_efield_kokkos.cpp: same hunk, the comment explaining the host-to-device copy of atom-style variables from this branch, the condition from the other side, which also covers an atom-style energy or potential variable writing into the fourth column of the same array. The three files that merged without conflict were checked by hand: the sync fixes in fix shake/kk touch init(), grow_arrays() and the forward comm and do not overlap the per-atom energy and virial work, and k_remap.setup() in the two regions is untouched by the removal of their now empty destructors.
stanmoore1
force-pushed
the
kk_bugfixes
branch
from
August 31, 2026 15:22
180eda7 to
e12884f
Compare
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 lammps#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.
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#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.
stanmoore1
added a commit
that referenced
this pull request
Aug 31, 2026
The bugfix branch behind PR #50 and this one arrived at the same fixes independently, so the two would have collided once that PR reaches develop. Adopt its text verbatim for every hunk both carry -- the MinFireKokkos view refresh, and the include-group binning in NBinKokkos and NBinSSAKokkos, down to the member names and the local-variable structure -- so the merge is clean. Only two hunks stay: the guards in binatomsItem() and binIDAtomsItem() that stop a non-numeric coordinate or an out-of-range bin index before the atomic update writes past the end of the bin arrays.
stanmoore1
added a commit
that referenced
this pull request
Aug 31, 2026
The bugfix branch behind PR #50 and this one arrived at the same fixes independently, so the two would have collided once that PR reaches develop. Adopt its text verbatim for every hunk both carry -- the MinFireKokkos view refresh, and the include-group binning in NBinKokkos and NBinSSAKokkos, down to the member names and the local-variable structure -- so the merge is clean. Only two hunks stay: the guards in binatomsItem() and binIDAtomsItem() that stop a non-numeric coordinate or an out-of-range bin index before the atomic update writes past the end of the bin arrays.
MinFireKokkos takes k_mass.view_device() and divides by it in every one of its integration kernels, but never syncs it. The masses are written through the plain host array, so on a GPU backend, and in the mixed and single precision builds where the device copy is a separate allocation, that copy stays zero. dtfm = dtf/mass is then inf and the coordinates turn into NaN on the first iteration. In a double precision host build the device view aliases the host array, which is why this only shows up elsewhere. The per-type masses are not covered by the sync()/modified() masks, so styles that read them on the device have to sync them by hand. Do that in init(), the same way fix nve/kk and fix nh/kk do.
…tion dtv was declared KK_FLOAT but is the receive buffer of MPI_Allreduce(&dtvone, &dtv, 1, MPI_DOUBLE, MPI_MIN, world); so in the mixed and single precision builds MPI_Allreduce() wrote eight bytes into a four byte object. dtvone and the Kokkos::Min reducer feeding it are already doubles, so make dtv one as well. The integration kernels use dtv, so add a KK_FLOAT alias for them, following the dtf_kk / dtf_init_kk / l_dmax_kk pattern already used in this file, rather than promoting the kernel arithmetic to double. Together with the preceding commit this makes min_style fire usable again in mixed and single precision. Verified with lj/cut on an fcc lattice for all four integrators, with and without abcfire: double and mixed precision now both reproduce the non-accelerated reference. Before, mixed precision ran to completion but never moved the atoms, reporting the initial energy and force norm at every step.
Dump::write() temporarily replaces atom->x, atom->v, and atom->image by plain host copies before calling domain->x2lamda(), domain->pbc(), and domain->lamda2x() for "dump_modify pbc yes". The Kokkos versions of these three functions operate on the Kokkos views (atomKK->k_x etc.) rather than on atom->x, so they silently worked on the real atom data instead of on the copies. With a triclinic box this produced several wrong results: - the copy handed to the dump was never converted to lamda coordinates, so Domain::pbc() wrapped real space coordinates against the lamda bounds [0,1). Coordinates came out shifted by up to one box length with image flags off by one. - a second dump at the same timestep saw the stale host view left behind by the device/host sync sequence and wrote lamda coordinates instead of real space ones. - the image flags of the real atoms were modified as a side effect. Detect that atom->x no longer aliases the Kokkos managed array and fall back to the base class implementations on the host in that case. Nothing changes for all other callers, where atom->x is the Kokkos array. Fixes lammps#4923, fixes lammps#4940 (cherry picked from commit 3daee37)
(cherry picked from commit bf37766)
(cherry picked from commit 71d991a)
All six fix-timestep-wall_region tests defined the region as "block EDGE EDGE EDGE EDGE EDGE EDGE" (and a sphere of radius 10 for the sphere test), which puts every atom farther from the wall than the cutoff. The stored reference data was all zeros: zero forces, zero stress and a zero global scalar, so the tests passed no matter what fix wall/region computed. Size the regions so the atoms are within range of the walls. The reference data was regenerated on the current develop base rather than taken from the original commit, so that the skip_tests and input_coeffs entries added since are kept. The regenerated files carry the "generated" tag, which marks reference data that has not been reviewed and validated yet. (cherry picked from commit 9bd5a7c)
(cherry picked from commit 6cd3c96)
…fter The test as written passed both with and without the fix in DomainKokkos::pbc()/x2lamda()/lamda2x(), so it did not guard anything. Its single atom sat inside the box, where the remapping is a no-op, and it compared the dump against a coordinate transform computed in the test rather than against what the non-accelerated styles write. Run the same input twice in one process, once without and once with the KOKKOS package, and compare the two dump files atom by atom. Let the atoms drift out of the box with a short run first, which is the situation in which the two disagree. Without the fix every atom now comes out one full box length away from the reference. Also drop the std::set_terminate()/_Exit() workaround, which skipped MPI_Finalize() and would have swallowed any real abort, and the Kokkos target guessing in CMakeLists.txt: the test uses only the public LAMMPS interface, so linking against lammps is enough. Co-authored-by: Mitch Murphy <alphataubio@gmail.com>
Say in the comment that the caller owns the x, v and image buffers it has substituted, so syncing those masks would overwrite the copies with the Kokkos data they were made from. Without that note the missing masks read like an oversight and invite a "fix" that reintroduces the bug. Also make detached_atom_x() const, it only inspects.
stanmoore1
force-pushed
the
kk_bugfixes
branch
from
September 1, 2026 17:43
f557f41 to
cff7722
Compare
Conflicts and how they were resolved: src/KOKKOS/region_remap_kokkos.h, region_block_kokkos.h, region_sphere_kokkos.h develop grew its own device-side remap of the tested coordinate into the periodic box (commit "remap positions into the periodic box in region block/kk and sphere/kk"), which is the same fix as this branch carried, with different names: RegionRemapKokkos::capture() instead of setup(), the member called boxremap instead of k_remap. develop's version is kept, since the five new region styles it also added use it. This branch's data race fix on the contact list was re-applied on top of it: the list is passed in as a Contact* by the caller, MAXCONTACT gives its length, and the surface functions became const. The destructors went with the device-side contact view, and point_on_line_segment(), inside_face() and find_closest_point() are const because the now-const surface path calls them. src/KOKKOS/domain_kokkos.cpp both sides changed the group-bit overloads of x2lamda() and lamda2x(): develop added the missing MASK_MASK to the device sync, this branch added the fall back to the host implementations when atom->x has been detached. Both are kept. src/KOKKOS/nbin_kokkos.cpp both sides taught the binning to honor an include group. develop's version is kept whole, since the merge had already taken its includegroup_* members in nbin_kokkos.h; dropping the conflicting statement would have left them set but never read. The bin sizing from actual occupancy is untouched. src/KOKKOS/Install.sh develop's new region styles were added. The region_cone/cylinder/ellipsoid/plane/prism styles develop adds implement only match_kokkos(), not the surface and contact path, and fix wall/region/kk still accepts block and sphere only, so the contact list fix does not have to be extended to them. Verified after the merge: KOKKOS double, mixed and single all build and reproduce the non-accelerated reference for min_style fire with every integrator, with and without abcfire; dump_modify pbc yes on a triclinic box matches the plain styles; fix wall/region/kk gives the same energy on 1, 2 and 4 threads; fix nve/kk with lj/cut/kk reproduces the plain thermo output; the six fix wall/region tests and the dump_modify pbc test pass.
…ackage The "neigh" keyword was parsed in the KOKKOS pair styles, which stripped it from the argument list before forwarding the rest to the base class. An input written for pace/kk therefore stopped with "Unknown pair_style pace keyword: neigh" when run in a build without the KOKKOS package. Parse the keyword in the non-accelerated styles instead and store the request in a base class member the KOKKOS styles read, which is how pair style snap handles its accelerator-only chunksize and parallelthresh keywords. The keyword is validated everywhere, so a typo is still caught without the KOKKOS package, but it has no effect there. This removes the settings() overrides in both KOKKOS styles along with the argument copying they needed, which also leaked the copied array on the error paths since error->all() throws when exceptions are enabled.
bond_harmonic/kk had its read_restart() fixed to assign the class members instead of declaring locals with the same names. angle_charmm/kk, bond_class2/kk, dihedral_charmm/kk and dihedral_charmmfsw/kk were cloned from the same template and still declare the locals, so after reading a restart file the k_* members are left default constructed while coeff() would have filled them. Nothing reads those members again on the current code paths, and the d_* device views keep the data alive through their reference count, so no result changes: a restart round trip with angle charmm and dihedral charmm gives the same energies before and after this commit, with and without the KOKKOS package. This removes the inconsistent state rather than a visible failure, and brings the four styles in line with bond_harmonic/kk. Co-authored-by: Mitch Murphy <alphataubio@gmail.com>
The regression test for issues lammps#4923 and lammps#4940 was never executed. The workflows that run the full test suite build with the "most" preset, which does not include the KOKKOS package, so the test was not even compiled there, and unittest-kokkos.yml selects tests by name with a regular expression the test did not match. Add it to that expression. The test also could not have passed in four of the six configurations that workflow builds. It asked for two threads unconditionally, which the serial backend refuses to start; since the LAMMPS constructor ran inside a captured stdout block, the error message was swallowed and the test binary died without reporting anything. Ask for two threads only when a threading backend is available. Finally it compared a KOKKOS run against a run without the package, which only agree to the compared tolerance when both use double precision. The mixed and single precision builds store positions as float and drift apart during the preceding run. Compare two dumps of the same timestep inside a single run instead, one with "pbc yes" and one without: wrapping an atom into the box must not move it, so both have to write the same unwrapped coordinates. That is the property the two issues are about, it needs no second run to compare against, and the remapping is done on host copies in double precision, so it holds for every backend and precision. Verified that the test still fails without the fix in domain_kokkos.cpp for double/OpenMP, double/serial and single/OpenMP, with unwrapped coordinates off by several box lengths, and passes in all three with it.
Conflicts in pair_pace_kokkos.cpp and pair_pace.rst between the neigh scratch fallback on this branch and the CPU-backend kernels for pace/kk on develop: keep both. develop turned host_flag into a compile-time constant, so the constructor no longer assigns it and only initializes host_fallback and the two neigh_scratch counters.
- nbin/kk: the max-occupancy reduction ran in the default execution space while d_bincount lives in DeviceType, so the host binning of a GPU build read a HostSpace view from a device kernel as soon as a bin overflowed. Give it a RangePolicy<DeviceType>. - fix shake/kk: the angle statistics counted the central atom as well, where FixShake::bond_force() counts only the two outer atoms, so the printed angle counts disagreed with the plain style. - verlet/kk: force_clear() zeroed both host copies of the force arrays on every step, on a GPU two full host passes that nothing reads in a run whose force styles are all on the device. The host copies are only read without a preceding sync on the host/device overlap path, where F_MASK is excluded from the sync and the host styles accumulate into them, so zero them only when a force style runs on the host side. The scan of the force styles is factored into host_force_styles(), which run() now uses as well for the same decision. The fused path skips force_clear() entirely and is unchanged. - run_style respa: refuse only a pair, bonded or kspace style that runs on a device, instead of every run with the KOKKOS package. The KOKKOS styles of a build without a device backend share their memory with the plain arrays and worked with rRESPA before, as do the plain styles. - min sd/kk: nvcc refuses an extended host/device lambda in a member function with private access, which set_search_direction() had.
The batched CPU derivative kernel that pace/kk gained on develop gathers
its per-lane buffers from KK_FLOAT views and hands them to
pace_batched_derivative(), which takes KK_FLOAT pointers, but declared
the buffers as double. That only compiles when KK_FLOAT is double, so
the single and mixed precision builds failed with
error: cannot convert 'double*' to 'const KK_FLOAT*' {aka 'const float*'}
Declare the buffers as KK_FLOAT.
While there, fix the copies of the per-atom core repulsion factor and
extrapolation grade from the host mirrors into the base class arrays.
They used memcpy with sizeof(double)*chunk_size, but the mirrors are
KK_FLOAT, so in a single or mixed precision build that read past the end
of a float buffer and filled the double arrays with garbage. Copy
element by element, which converts.
- min quickmin/kk: sync the per-type masses to the device in init() - fix wall/harmonic/outside/kk: size the per-atom virial dual view in v_setup_peratom() with alloc = 0, after v_init() set vflag_atom - compute reaxff/atom/kk: TAG_MASK/EMPTY_MASK datamasks and a host sync of the tags before compute_local(), instead of the ALL_MASK defaults - compute temp/sphere/kk: sync and claim around the bias compute's remove_bias_all()/restore_bias_all(), as temp/deform/kk does - pair mliap/kk: launch the comm packing kernels in RangePolicy<DeviceType> so the host instantiation does not run them on the device - fix nvt/nph/npt/sphere/kk: remove leftover debug output - comm_style tiled (plain and KOKKOS): reverse_comm(Pair*) honored the size argument backwards; match CommBrick - fix shake/kk: pack_forward_comm() no longer claims a host write it does not make; drop the now unused i0 - pair uf3/kk: size the centroid virial with maxcvatom and free it in the destructor - fix propel/self/kk: v_init(vflag,0), the per-atom virial is a dual view - pair eam/kk: size the HIP team scratch with sizeof(KK_FLOAT) - fix addtorque/atom/kk, settorque/atom/kk: release the base class array before the dual view takes it over - fix qeq/reaxff/kk: struct/class tag mismatch on a friend declaration
- fix cmap/kk: tally the global virial and the per-atom energy and virial the way FixCMAP::post_force() does through Fix::ev_tally(): each crossterm split over its five atoms, only the owned ones counted. The per-atom arrays are dual views sized after ev_init(eflag,vflag,0). Before this the fix's virial was zero and compute pe/atom and stress/atom silently omitted the CMAP term. - pair hybrid/scaled under KOKKOS: a pair style that runs on the device but does not implement the KOKKOS packing communicated through a null KokkosBase pointer; both comm classes now take the host path for it - remap/kk: wait for all receives before copying the shared scratch buffer to the device (the grid3d fix), fence the pack kernel before MPI reads the device buffer, copy only the message's own stretch of the host send buffer, and MPI_STATUSES_IGNORE for the Waitall - pair dpd/kk: cache update->ntimestep on the host instead of dereferencing the Update object inside the team kernel - pair pod/kk: drop the commented-out debugging dumpers, whose fwrite of a KK_FLOAT mirror with sizeof(double) was wrong in single precision
- fix spring/kk, fix wall/flow/kk, compute ave/sphere/atom/kk: sync the per-type masses to the device before reading them, instead of relying on the integrator having done it - compute temp/profile/kk: sync the legacy host side of the masses, which is what the host loop reads; sync<LMPHostType>() refreshes the other side - atom map: map_one() claims its host write; map_clear() releases both sides before it overwrites the whole array - compute inertia/kk: the host sync before inertia_extended() now covers the radius, the ellipsoid indices and the bonus data it reads - fix nh/sphere/kk, compute temp/sphere/kk: sync the host copies before the radius/mask loops of init() and dof_compute() - fix rigid/small/kk and its Nose-Hoover variants: v_init(vflag,0), the per-atom virial is a dual view the fix zeroes itself - pair gran/hooke/kk, gran/hooke/history/kk, gran/hertz/history/kk: stop claiming and copying a per-atom energy that is never written - comm_style tiled/kk: MPI_STATUSES_IGNORE for Waitall, and take the send pointer after the pack that may resize the buffer - NBin gets a copymode guard, set by nbin/kk and nbin/ssa/kk while a copy of the object lives in a device functor; likewise MLIAP_SO3Kokkos - region ellipsoid/kk: rotate() locals no longer shadow the semi-axes - pair lj/gromacs/kk: drop the cut_inner dual views nothing reads
pair_style hybrid/scaled/kk kept its accumulators as plain host arrays, so every sub-style needed the forces (and torques) copied back and forth between host and device on each timestep. Do the work in Kokkos kernels instead: - fsum/tsum become dual views and the save/clear/scale-add/restore steps run as parallel_for kernels - init_style() picks the memory space for the sum: the device when every sub-style runs there, so the forces stay resident on the device, and the Kokkos host space otherwise - virial_fdotr is only synced to the device when it is actually needed Per-atom scale factors from an atom-style variable are now communicated through KokkosBase::pack/unpack_forward_comm_kokkos. Previously the forward communication of a Device pair style dispatched to the device path and dereferenced a null KokkosBase, which crashed in compute() and in single(). The variable itself is still evaluated on the host, so the local values are uploaded and the ghost values filled in by the exchange; single() and born_matrix() read them back from the dual view. Verified against pair_style lj/cut with sub-style scale factors summing to one (constant, equal-style, and spatially varying atom-style), with atom_style sphere for the torque path, with a mix of /kk and plain sub-styles, and against the CPU hybrid/scaled for compute group/group. The force-style unit tests pass in double/legacy-layout and in mixed-precision/default-layout builds.
Respa::init() called into the styles' KOKKOS execution spaces itself and refused any force style on the device. Move the decision into KokkosLMP::respa_check(), so that respa.cpp keeps a single guarded call and the dummy interface a no-op, and make the criterion what rRESPA actually needs: every participant on the host side of the atom data. That admits a build without a device backend as before, and now also a device build whose styles all use the /kk/host suffix with comm, sort and atom/map on the host; it refuses a style, fix or compute in the device execution space and device-side communication, sorting or atom map, and names what to change. Drop the versionchanged directive from the run_style documentation: the KOKKOS package never supported rRESPA, so nothing changed for users.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
A collection of KOKKOS bug fixes, most of them host/device data coherence
defects: an array written on one side and never claimed with
modified_*(),or read on the side that is not current. These are invisible on a host-only
Kokkos build, where the two sides share an allocation, and show up as wrong
forces or energies on a real device.
Also included are two fixes to the KOKKOS neighbor binning, which ignored
neigh_modify includeentirely, and a crash fix forfix qtpie/reaxff.The changes fall into these groups:
pair dpd/fdt/energy,min_style fire,fix deform,fix shake,fix nh,comm_tiled,grid3d,bond quartic,the
improperstyles,Atom::set_mass(), and thesetcommand.force_clear(): the host side of the forces was cleared only in thebranch without an include group, and only for the forces; the torques and
the SPIN forces were never cleared on the host in any branch.
neigh_modify include:NBinKokkosandNBinSSAKokkosbinned everyowned atom and every ghost, where the plain
NBinclasses bin only theatoms an include group's pairs are built from. Atoms outside the group
therefore appeared in the neighbor lists of atoms inside it.
~NeighBondKokkos()ran on the functor copy Kokkos handseach loop and released lists that belong to
Neighbor, so a bond, angle,dihedral or improper style without KOKKOS support crashed in setup.
DomainKokkos::remap_all()andimage_flip()were shadowing the baseclass methods rather than overriding them, so their kernels had never run.
fix qtpie/reaxffread the pair style's device neighbor list andsegfaulted; it now uses its own host-built list.
run_style respais refused with the KOKKOS package (see BackwardCompatibility).
Related Issue(s)
Several of the divergences reported in lammps#5037 are of the shapes fixed
here.
Author(s)
Stan Moore, Sandia National Laboratories
Licensing
By submitting this pull request, I agree, that my contribution will be
included in LAMMPS and redistributed under either the GNU General Public
License version 2 (GPL v2) or the GNU Lesser General Public License version
2.1 (LGPL v2.1).
Artificial Intelligence (AI) Tools Usage
AI tools were used for parts of this pull request and this section replaces
the default statement accordingly.
Most of the code changes here were generated by Claude Code, which was also
used to find the underlying defects: it ran a KOKKOS host/device
sync-debugging build over the
examples/tree, compared the thermodynamicoutput against a stock KOKKOS build, and root-caused each divergence. The
AI-generated changes are the ones touching
min_fire_kokkos,fix_deform_kokkosandfix_deform,fix_shake_kokkos,pair_dpd_fdt_energy_kokkos,comm_tiled_kokkos,domain_kokkosanddomain.h,modify_kokkos,atom_kokkosandatom.h,set.cpp,respa.cpp,fix_qtpie_reaxff.cpp,fix_numdiff*,compute_born_matrix,nbin_kokkos,nbin_ssa_kokkos, and theforce_clear()andinit_topology_kk()changes inverlet_kokkosandneigh_bond_kokkos.The remaining changes -- the
improperstyle syncing,bond_quartic_kokkos,grid3d_kokkos,fix_nh_kokkos::remap(), theNeighBondKokkosdestructor,and the first host-side clear in
force_clear()-- were written by a human.Every change, AI-generated or not, was reviewed and verified as described
below.
Backward Compatibility
One intentional break:
run_style respanow errors out when the KOKKOSpackage is active, where it previously ran and silently produced wrong
forces. rRESPA keeps per-level force copies and clears and sums them
through the plain LAMMPS arrays, without the host/device transfers that
run_style verletperforms, so on a GPU the results are wrong with noindication that anything is amiss. An input combining the two will now stop
with an error instead of producing bad numbers.
doc/src/run_style.rstisupdated with a
.. versionchanged:: TBDnote.No other input syntax changes.
Implementation Notes
How correctness was verified:
examples/were run on a stock KOKKOS buildand on a
KOKKOS_DEBUG_SYNCbuild and their thermodynamic output comparedcolumn by column. Before these fixes 47 inputs diverged; after them none
do.
force_clear()fixes were checked on a build whereOpenMP stands in for the device and Serial for the host, so the two sides
are genuinely distinct on a CPU. With the host clear removed the total
energy of a fourmol run climbs 371.81 -> 372.61 over four steps where the
non-KOKKOS run holds 371.68; with it in place the two agree to the last
digit. The same A/B confirms the include-group branch specifically.
atom_modify firstandneigh_modify include,E_pairwas -31960.18 against -31864.22 withoutKOKKOS, from the first step and roughly constant thereafter -- the shape of
a different set of pairs rather than of forces accumulating. The two runs
now agree to the last digit.
and the dpde/dpdrx shardlow examples are unchanged against their
pre-change output and match their non-KOKKOS runs.
make check-whitespace check-permissions check-fmtlib check-homepage check-errordocsall pass.Two notes for reviewers:
Domain::remap_all()andDomain::image_flip()becoming virtual activatesa
DomainKokkoskernel that had never run. That kernel omits thedeform_vremapvelocity correction, so a host fallback is used for thatcase; finishing the kernel instead may be preferable.
Atom::sync_host_arrays()/Atom::modified_host_arrays()are new virtualhooks that let a style declare a host write of the per-atom arrays without
depending on KOKKOS. They are currently undocumented in the developer
guide.
Testing has been on CPU only, with a single MPI rank; the fixes have not been
exercised on a GPU.
Post Submission Checklist
Generated by Claude Code