Review fixes for upstream PR #5056 (dispersion/d3 OPENMP and KOKKOS variants) - #52
Closed
stanmoore1 wants to merge 17 commits into
Closed
Review fixes for upstream PR #5056 (dispersion/d3 OPENMP and KOKKOS variants)#52stanmoore1 wants to merge 17 commits into
stanmoore1 wants to merge 17 commits into
Conversation
- remove trailing whitespace in MathSpecial::powauto()
- register pair_dispersion_d3_kokkos.{cpp,h} in src/KOKKOS/Install.sh so the
traditional make build installs the style (the base style lives in the
optional EXTRA-PAIR package, hence the dependency argument)
- add the accelerator code letters to the dispersion/d3 entry in
doc/src/Commands_pair.rst
- add the missing ".. include:: accel_styles.rst" block to
doc/src/pair_dispersion_d3.rst, which doc/utils/check-styles.py requires
once a page lists accelerated style index entries
The KOKKOS version was missing the "n == 0" shortcut, so powauto(0, 0) returned 0 on the device and 1 on the host, and it called round() and pow() unqualified, which resolves to the global double overloads in device code and promotes for a KK_FLOAT argument. Use Kokkos::round() and Kokkos::pow() as required for device kernels, and take the arguments by const reference like the neighboring powint(). Add the doxygen comment block that every other function in these headers has, the matching entry in the programmer guide, and unit tests covering the integer, fractional and zero cases against powint() and std::pow().
The /omp variant tallied the virial incorrectly. compute() passed vflag_either to reduce_thr(), which bit-tests its argument against the raw VIRIAL_* flags, while ev_setup() reduces vflag_either to 0/1. So "1 & (VIRIAL_ATOM|VIRIAL_CENTROID)" was always 0 and the per-atom virial was never reduced across threads, and "1 & (VIRIAL_PAIR|VIRIAL_FDOTR)" was 1, so the per-thread fdotr virial was added on top of the trailing virial_fdotr_compute() call. With "compute stress/atom" active the global pressure came out 3.4x too large and the per-atom stress lost everything computed by threads 1..N-1. With only the fdotr virial requested the two errors cancelled, which is why the force styles test did not catch it. Pass vflag and drop the trailing virial_fdotr_compute(), as in every other OPENMP pair style. Energy, pressure tensor and summed per-atom stress are now bit identical to the plain style for 1, 2, 4 and 8 threads. The thread private cn and dc6 accumulators were std::vector<double> of atom->nmax, allocated and zeroed per thread on every timestep and reduced with an "omp atomic" loop over all of nmax, which also read the elements past nall that no memset had initialized. Switch to the scheme pair_adp/omp uses: grow cn and dc6 to nthreads*nmax, let each thread accumulate into its own nall long slice, and reduce with data_reduce_thr(). PairDispersionD3::get_dC6() returned a pointer to a function local "static double c6_res[3]", so it was not usable from a threaded or a device loop. Both accelerated variants worked around that by carrying a verbatim ~60 line copy of it. Give the base version an output parameter instead and delete both copies, so a future fix to the interpolation cannot miss a variant. The KOKKOS variant also carried a copy of coeff() with the device transfer appended, which pulled the 2.9 MB d3_parameters.h - and with it a second copy of the 32385 x 5 reference C6 table - into a second translation unit. Call the base version and keep only the device transfer. The constants K1, K3, AUTOANG and AUTOEV existed in three places, and in pair_dispersion_d3_kokkos.h they sat at global namespace scope in a style header, which is what akohlmey objected to for the pow() helpers. Define them once in namespace LAMMPS_NS::DispersionD3 in the base header, together with AUTOANG6, which all three variants recomputed inside the innermost loop of get_dC6(). Declare allocate() and calc_coordination_number() virtual in the base header and mark the overrides, as required for a KOKKOS style; they were hidden, not overridden, and only happened to work because coeff() and compute() are overridden as well. The "unknown damping code" check moves to PairDispersionD3::init_style(), so neither the threaded loop nor the device kernel has to report an error from a context that cannot. Further fixes found while reviewing: - coeff() reallocated c6ab and re-read the reference table whenever the largest CN grid index exceeded 4. The largest encoded index in d3_parameters.h is 4, so the branch was unreachable - and had it ever fired, the first read_c6ab() call would already have written out of bounds on the fixed 5x5 buffer. Drop the branch, keep the max_mxci value the KOKKOS view sizing needs, and bounds check the grid indices in read_c6ab(), where they are produced. - coeff() leaked the malloc()ed element list when an unknown element name triggered error->all(); use std::vector. - k_cutsq was a plain dual view of KK_FLOAT while cutsq is double **, which only compiles when KK_FLOAT is double. Use DAT::ttransform_kkfloat_2d, as in pair_lj_cut_kokkos. - The kernel init()/join() reduction hooks zeroed and summed only evdwl and v[6] of EV_FLOAT, leaving ecoul indeterminate. Drop them and use the default constructor and operator+= of s_EV_FLOAT. - The "den > 1.0E-99" fallback threshold in the device version of get_dC6() underflows to zero in a single precision build; use the smallest normalized KK_FLOAT instead. - reverse_comm_device was never set, so the device pack/unpack methods were dead code and every newton step round-tripped cn and dc6 through the host. - exp() was called unqualified in three device kernels. - Drop the unused f, d_rcov, nall and cn_thr members shipped with the kernel functors, the unused sbmask() of the pair class, the dead grid dimension constants, the unused includes, the unused thr_evflag member, and the constructor assignments the base class already makes. - Collapse the six way copy of the coordination number kernel launch into a template, and the byte identical bj/bjm damping cases into one.
…d3/kk The kernels forced ScatterNonDuplicated whenever DeviceType was LMPDeviceType. In a GPU build that matches NeedDup_v anyway, but in a CPU only build LMPDeviceType and LMPHostType are the same type, so the "half" neighbor style with more than one thread - which maps to HALFTHREAD - fell back to atomics instead of the duplicated scatter views every other pair style uses there. Use NeedDup_v<NEIGHFLAG,DeviceType> directly, as PairComputeFunctor does, and drop the matching special case from the need_dup computation in compute() so the two stay consistent (they must agree, or ScatterViewHelper hands the functor the wrong half of the dup_/ndup_ pair). Also drop two variables left unused by the previous commit, so both changed OPENMP and EXTRA-PAIR files compile clean under -Wall -Wextra. Checked against the plain style with "compute stress/atom" active: energy, pressure tensor and summed per-atom stress are bit identical for neigh half and neigh full, newton on and off, at 1 and 4 threads.
stanmoore1
force-pushed
the
claude/pr-5056-review-conflicts-hvlvgo
branch
2 times, most recently
from
August 25, 2026 18:53
f6cdb46 to
edfbfb3
Compare
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
Upstream lammps/lammps#5056 ("Add OPENMP and KOKKOS package variants of the pair_style dispersion/d3", by Marc L. Descoteaux et al.) has not been updated since June. This branch is a strict fast-forward from that pull request's head (
b454b455): the 13 original commits unchanged, with 4 review-fix commits appended and nothing else.git merge --ff-onlyfrom lammps#5056's branch onto this one works.Against lammps#5056's head the delta is 13 files, +296/-506 — just the fixes. Against
developit is the full feature, 14 files, +2567/-80.This is a staging branch in this fork, meant as a reviewable delta to hand back to the author of lammps#5056, not as a competing submission to
lammps/lammps.Merging
developafterwards. One four-line conflict, insrc/EXTRA-PAIR/pair_dispersion_d3.h:Resolution:
virtual void allocate();, withmemory_usage()keeping the new home in the public section thatdevelopgave it. The conflict is unavoidable —developdeletes thememory_usage()declaration from the protected block while this branch edits the line directly beneath it, so any three-way merge sees one overlapping hunk. Everything else auto-merges.Related Issue(s)
Companion to lammps#5056. Does not close it.
Author(s)
The
dispersion/d3/ompanddispersion/d3/kkstyles are by Marc L. Descoteaux, Yizhong R. Hu, Ulrik Unneberg and William C. Witt (Harvard University), per the contributing-author blocks in the new files. The four review-fix commits are by 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, in two separate places.
The original lammps#5056 commits carry the authors' own disclosure: the KOKKOS pair style was heavily revised from a human-written version using the Codex VSCode extension (Codex-5.2, Codex-5.3 and GPT 5.4 models), then re-examined and revised by a human developer.
The four review-fix commits appended here were produced with Claude Code, driven and reviewed by a human developer. That covers the changes to
src/EXTRA-PAIR/pair_dispersion_d3.{cpp,h},src/OPENMP/pair_dispersion_d3_omp.{cpp,h},src/KOKKOS/pair_dispersion_d3_kokkos.{cpp,h},src/KOKKOS/math_special_kokkos.h,src/math_special.h,src/KOKKOS/Install.sh, the threedoc/src/*.rstfiles andunittest/utils/test_math_special.cpp. Every behavioral claim below was verified by building and running the code, not by inspection alone.Backward Compatibility
No input scripts break. The
zerodamping keyword remains accepted as an alias fororiginal.MathSpecial::powauto()andMathSpecialKokkos::powauto()are new; the KOKKOS one changes a signature added earlier in this same branch and never released.The only user-visible change to an existing style is that
pair_style dispersion/d3/ompnow reports the correct virial (see below); results that were being computed wrongly will change.Implementation Notes
Two real bugs in
dispersion/d3/omp.compute()passedvflag_eithertoreduce_thr(), which bit-tests its argument against the rawVIRIAL_*flags whilePair::ev_setup()reducesvflag_eitherto 0/1. So1 & (VIRIAL_ATOM|VIRIAL_CENTROID)was always 0 and the per-atom virial was never reduced across threads, and1 & (VIRIAL_PAIR|VIRIAL_FDOTR)was 1, so the per-thread fdotr virial was added on top of a trailingvirial_fdotr_compute()call. With only the fdotr virial requested the two errors cancel, which is why the force-styles test did not catch it. Withcompute stress/atomactive, on a 64-atom Si cell:/ompas in lammps#5056/ompafter this branchThe fix passes
vflagand drops the trailingvirial_fdotr_compute(), as in every other OPENMP pair style (pair_reaxff_ompis the sole exception).Thread-private accumulators.
cnanddc6werestd::vector<double>ofatom->nmax, allocated and zeroed per thread on every timestep, then reduced with anomp atomicloop over all ofnmax— which also read the elements pastnallthat no memset had initialized. Replaced with the schemepair_adp/ompuses:cnanddc6are grown tonthreads*nmax, each thread accumulates into its ownnall-long slice, anddata_reduce_thr()sums the slices.A KOKKOS build break in single/mixed precision.
k_cutsqwas declaredtypename AT::tdual_kkfloat_2dwhilecutsqisdouble **, somemoryKK->create_kokkos()only binds whenKK_FLOATisdouble. NowDAT::ttransform_kkfloat_2d, as inpair_lj_cut_kokkos.Other KOKKOS fixes. Unqualified
exp()in three device kernels; theden > 1.0E-99fallback threshold underflowing to zero in a float build;init()/join()leavingEV_FLOAT::ecoulindeterminate;reverse_comm_devicenever set, so the device pack/unpack methods were dead code and every newton step round-trippedcnanddc6through the host; and scatter-view duplication forced off for HALFTHREAD in a CPU-only build, whereLMPDeviceTypeandLMPHostTypeare the same type.Removing the copy-paste.
PairDispersionD3::get_dC6()returned a pointer to a function-localstatic double c6_res[3], so it was not usable from a threaded or a device loop — which is why both accelerated variants carried a verbatim ~60-line copy of it. It takes an output parameter now and both copies are gone. The KOKKOS variant also carried a copy ofcoeff(), which pulled the 2.9 MBd3_parameters.h(and with it a second copy of the 32385 x 5 reference C6 table) into a second translation unit; it calls the base version now.K1,K3,AUTOANGandAUTOEVexisted in three places, and inpair_dispersion_d3_kokkos.hthey sat at global namespace scope in a style header; they are defined once innamespace LAMMPS_NS::DispersionD3in the base header, together withAUTOANG6, which all three variants were recomputing inside the innermost loop ofget_dC6().allocate()andcalc_coordination_number()are nowvirtualin the base header, as required for a KOKKOS style.An unreachable branch.
coeff()reallocatedc6aband re-read the reference table whenever the largest CN grid index exceeded 4. Decoding the table shows the largest encoded value ind3_parameters.his 482, i.e. grid index 4, so the branch never fires — and had it ever fired, the firstread_c6ab()call would already have written out of bounds on the fixed 5x5 buffer. The branch is dropped,max_mxci(which the KOKKOS view sizing needs) is kept, and the grid indices are bounds-checked inread_c6ab()where they are produced.Packaging and docs. Trailing whitespace in
math_special.h(failsmake check-whitespace); missingsrc/KOKKOS/Install.shentries;dispersion/d3missing its(ko)accelerator letters inCommands_pair.rst; missing.. include:: accel_styles.rstinpair_dispersion_d3.rst, whichdoc/utils/check-styles.pyrequires once a page lists accelerated index entries.Two things deliberately not done. The unit-test YAML was not regenerated: the maximum relative error against the stored reference is 2e-12 against
epsilon: 7.5e-08, and thezero->originaledit is a numerical no-op (both map todampingCode 1), so the stored numbers remain valid and refreshingdate_generatedwould be less accurate, not more. Andev_setup_thr()was not added to the second parallel region: it memsetseatom_pair/vatom_pairand would erase the first phase's tallies. The existing structure is correct and is now documented as such.The package-wide
-Wimplicit-float-conversion/-Wdouble-promotionsweep that landed indevelopafter lammps#5056 branched (commit708e1dba18) has not been applied to the new KOKKOS style. It still has roughly 70 narrowing sites (1.0f/3.0/0.5literals mixed withKK_FLOAT,AUTOANG/AUTOEV/K1doubles folded into kernel expressions, host-sidespecial_lj/virial/comm-buffer assignments). Only the subset that is a build break or a behavior change rather than a warning was fixed here. This is the main outstanding item, and it only becomes relevant oncedevelopis merged in.Verification, on this exact branch (EXTRA-PAIR + OPENMP + KOKKOS, Serial and OpenMP backends, CMake): build clean, no warnings on the changed files.
PairStyle.plain,.omp,.kokkos_omp,.extractand.extract_omppass;MathSpecial35/35 including the newpowautocases. Energy, pressure tensor and summed per-atom stress are bit-identical to the plain style for/ompat 4 threads and for/kkwithneigh half(newton on) andneigh full(newton off) at 4 threads.make check-whitespaceandcheck-permissionspass.The same change set was additionally verified on top of current
developbefore this branch was rebased onto lammps#5056 — the two differ only in hunk offsets, the changed lines are byte-identical. There it also passed 353/353 ctest, the/ompsweep at 1, 2, 4 and 8 threads, the/kksweep with newton on and off, the traditional make build, and the doccheck-styles.py, char, role and package checks.Sphinx and doxygen are not installed in the environment this was prepared in, so
make htmlandmake spellingcould not be run. That needs checking before this goes upstream.Post Submission Checklist
Further Information, Files, and Links
pow()helpers) and thestd::vectorscratch buffers, which were switched away fromnew[]/delete[]as asked but were stillnmax-sized and reallocated every timestep.