Split step clean - #375
Open
mewall wants to merge 71 commits into
Open
Conversation
Ran for t2, t3, t34, t4.
Added new logic to top of mdstep loop.
Added input file for usinf kernel and LangeDynamics.
…r than minimization steps.
Scale kappa with (timestep)^2 as required by κ = Δt²ω². Result: Stable through 100+ split-step cycles. Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Changed MD loop to continue until print_mdstep reaches the requested number of steps, rather than using a fixed internal step count. This ensures that MDSteps= in input.in controls the number of output steps at the user timestep, regardless of how many split-steps occur. Before: MDSteps=50 with frequent splits would give < 50 output steps After: MDSteps=50 always gives exactly 50 output steps Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Add interpolation-based approach for handling non-uniform timesteps. When timesteps vary, interpolate historical charges n_0...n_5 from non-uniform grid to uniform grid, then apply fixed coefficients C0-C5. This is mathematically equivalent to recomputing adaptive coefficients but simpler to implement. Key changes: - Add dt_history(5) and nsteps_taken to xlbo_type for tracking timesteps - Implement prg_xlbo_interpolate_charges() using 5th-order Lagrange polynomial - Modify prg_xlbo_nint, prg_xlbo_nint_kernel, prg_xlbo_nint_kernelTimesRes to use interpolation when nsteps >= 6 and dt provided - Automatic activation after 6 MD steps, graceful fallback for early steps Test results: - Uniform timesteps (0.25 fs): Identical to original (0.0 eV difference) - Variable timesteps (0.5 fs, 60 split-steps): Same accuracy as recomputing coefficients (0.0 eV difference), confirming mathematical equivalence Builds on kappa scaling fix from commit 7f62eec. Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Replace static dt_base variable with timestep ratio parameter to fix critical bug in adaptive timestepping. The dt parameter now represents the ratio of current timestep to base timestep (1.0 or 0.5), not an absolute timestep value. Changes: - Remove static dt_base from all three XLBO integration subroutines (prg_xlbo_nint, prg_xlbo_nint_kernel, prg_xlbo_nint_kernelTimesRes) - Simplify kappa scaling: kappa_use = kappa * dt^2 (dt is now ratio) - Update dt_history to store ratios instead of absolute timesteps - Modify mdloop to pass lt%timestep/user_timestep ratio to XLBO calls Impact: - Fixes incorrect kappa scaling when adaptive timestepping triggers - Enables proper energy conservation with variable timesteps - Tested with water example: energy drift < 0.005% over 50 steps Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Add a logical input parameter AdaptiveTimeStep to the GPMD parser to
allow users to enable/disable adaptive timestep splitting behavior.
When disabled, simulations run with fixed timestep throughout.
Changes:
- Add adaptive_timestep field to gpmd_type in gpmdcov_parser.F90
- Add AdaptiveTimeStep= keyword to logical parser keys (default: false)
- Gate timestep splitting logic in mdloop with adaptive_timestep flag
Usage in input.in:
GPMD{
AdaptiveTimeStep= T ! Enable adaptive timestep (default: F)
...
}
Benefits:
- Allows direct comparison of fixed vs adaptive timestep methods
- Provides control for benchmarking and validation studies
- Maintains backward compatibility (default is off)
Tested:
- AdaptiveTimeStep=F: No splitting occurs (fixed timestep)
- AdaptiveTimeStep=T: Splitting occurs when max displacement > 0.02 Å
Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Critical bug fix: kappa was being incorrectly scaled by dt^2 even when using Lagrange interpolation. When interpolation is active, charges are already interpolated to a uniform grid, so kappa should NOT be scaled. Changes: - Move kappa scaling logic after interpolation decision - Only scale kappa when NOT using interpolation - Add check: if (present(dt) .and. .not. use_interpolation) Impact: - Fixed timestep: K=10 works correctly (slight improvement vs K=5) - Adaptive timestep: K=10 now stable and completes 1930 steps * Before fix: catastrophic failure at step 315 * After fix: 27.6% lower energy drift than K=5 (-75.7 vs -104.5 meV/ps) Testing (1000 step water simulation, dt=0.5 fs base): - K=5 fixed: drift = +33.95 meV, RMS = 26.1 meV - K=10 fixed: drift = +33.58 meV, RMS = 25.9 meV - K=5 adaptive: drift slope = -0.052 meV/step - K=10 adaptive: drift slope = -0.038 meV/step (28% better) Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Bug: prg_xlbo_nint_kernelTimesRes call at line 381 was missing the dt parameter (timestep ratio). This caused incorrect history tracking when adaptive timestepping was used. Fixed: - Added lt%timestep/user_timestep parameter after xl - Changed n_6=n_6 keyword syntax to positional n_6 syntax for consistency This ensures dt_history is properly maintained for both K=5 and K=10 XLBO with variable timesteps. Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Only allow timestep splitting after K=5 has 6 steps or K=10 has 11 steps of history. This ensures Lagrange interpolation is always available when variable timesteps are used, avoiding fallback to kappa scaling which caused instability. Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
Kappa should always scale with dt^2 when variable timesteps are used, regardless of whether interpolation is active. The kappa term cc*kappa*(charges-n) provides extended Lagrangian coupling and must scale with the actual timestep being taken in the current MD step. Previous logic incorrectly disabled kappa scaling when interpolation was active, but interpolation only affects the dissipation history term, not the coupling term. Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
…dt/2 uniform grid Replaced high-order Lagrange polynomial interpolation with cubic spline interpolation to reduce oscillations when interpolating charge history from non-uniform to uniform time grids during adaptive timestepping. Key changes: - Added cubic_spline_coeffs() and cubic_spline_eval() helper functions - Modified prg_xlbo_interpolate_charges() (K=5) to use cubic splines - Modified prg_xlbo_interpolate_charges_K10() (K=10) to use cubic splines - Fixed binary search in cubic_spline_eval() to handle descending arrays - Use fixed dt/2 = 0.5 uniform target grid for consistent interpolation - Force initial timestep splitting for first 4 (K=5) or 6 (K=10) print_mdsteps to build history at dt/2 spacing before normal adaptive behavior Cubic splines provide C² continuity and avoid the oscillations that can occur with 10th-degree Lagrange polynomials for K=10. The fixed dt/2 uniform grid ensures consistent interpolation target regardless of whether current step is full or half timestep. Testing shows stable residuals (~2×10⁻⁶) but increased energy drift (~10⁻⁴ eV/step) compared to uniform timesteps. System remains stable for 1000+ steps with K=10 and adaptive timestepping enabled. Co-Authored-By: Claude Sonnet 4.5 <[email protected]>
De-duplicate coefficient setup into xlbo_compute_coeffs/xlbo_push_dt_history, remove dead vars and stale interpolation comments, drop debug prints in mdloop. Numerics unchanged. Co-Authored-By: Claude Opus 4.8 <[email protected]>
The neighbor-list builders sized their lists from a hardcoded
MAX_DENSITY=0.25 atoms/Ang^3 macro, too low for dense solid-state
systems. Replace the macro with a module variable in
gpmdcov_neighbor_mod set once from a new GPMD{} MaxDensity= input
keyword (default 0.25, so existing runs are unchanged). The value is
propagated in gpmdcov_init right after parsing; no per-routine argument
threading is needed since gpmdcov_init already uses the neighbor module
and there is no use-dependency between the parser and neighbor modules.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Both prg_replicate and prg_replicate_system used the 1-based Fortran loop
counters (i,j,k in 1..n) directly as cell multipliers, so replica cells were
placed at offsets {L,2L,...,nL} instead of {0,L,...,(n-1)L}. This shifted the
entire supercell by one full lattice vector in every direction: cell (1,1,1)
landed at +L rather than the origin, and the last layer stuck out of the
[0,nL box. The error was masked whenever TranslateAndFoldToBox (default TRUE)
folded atoms back under PBC, so it surfaced only when folding was off or a
downstream step assumed atoms in [0,Lbox).
Use (i-1),(j-1),(k-1) as the multipliers. Also fix prg_replicate_system's
coordinate formula to sum lattice_vector(v,:) over vectors (matching
prg_replicate and the supercell scaling below), so non-orthogonal cells
replicate correctly too.
Verified with a 2x2x2 replicate of an 8 Ang cell: the first replica now sits
exactly at the input coordinates and the 8 cells tile (0,0,0)..(8,8,8) within
the 16 Ang supercell.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
EOF
)
Halve alpha for the 16 half-ending patterns (even indices) in XLBO_K5_alpha so alpha*d_K per unit physical time matches full steps. Previously two half-steps applied ~2x the reference dissipation over the same interval as one full step, an asymmetry that drives slow energy drift on adaptive runs. Coefficient tables untouched (moment constraints preserved); halved values stay within stability caps (min margin 2.4x). Co-Authored-By: Claude Opus 4.8 <[email protected]>
Replace the ad-hoc hybrid XLBO_K5_alpha table with a single rule applied
uniformly to all 32 patterns:
alpha_idx = min(0.054 * dt_current / d_K(idx), 0.15)
so the table supports any split history rather than being tuned to a
particular split mix. 28/32 patterns are rate-matched; idx 1,3,5,6,7 hit
the 0.15 ceiling.
Ceiling chosen empirically (scripts/Niklasson_JCP_2009_table_I/empirical_stability.py):
stability is conditional on SCF convergence quality. For a stable run
|q[n]-n| < 1e-3 => gamma >= ~0.999, where this table is stable (peak amplitude
~1.3 over 20000 steps) with margin as gamma degrades. Periodic-orbit eigenvalue
caps are pathological here, so validation uses quasi-random histories.
The previous hybrid table is retained, commented out, for easy reversion.
Note: aggregate split-boundary friction is barely changed (idx 7 is
stability-locked below its rate-match target), so this is a robustness/
code-quality change and is unlikely to reduce energy drift on its own.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Add examples/gpmdk/docs/adaptive_time_step.md consolidating the adaptive
time-step documentation: what the split-step does and how to enable it
(GPMD{ AdaptiveTimeStep=T }); how the K=5 integrator handles the mixed
full/half sequence via per-pattern coefficient tables with constant kappa
and a symmetric ka_scale (which keeps the Verlet backbone reversible);
the history-independent general alpha rule and the physical rate-match
rationale behind it; energy-conservation methodology and results for
300-atom water and the large Mac1 system; the fixed-schedule control
experiment identifying state-dependent step selection as the drift
mechanism; NVT-vs-NVE usage guidance; a "Rejected approaches" section
recording dead ends (charge interpolation, per-pattern kappa, kappa-scaling
removal with C0-C5 at 11 points, periodic-orbit eigenvalue caps); and the
planned reversible backtracking scheme.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Restore examples/gpmdk/run/water/input.in to the master baseline (drop local test toggles: partition counts, coords file, TimeStep/MDSteps) and remove scratch my_waterInput*.in files. How to enable the adaptive time step is documented in examples/gpmdk/docs/adaptive_time_step.md. Co-Authored-By: Claude Opus 4.8 <[email protected]>
Include the scripts used to compute the values that appear in src/prg_xlbo_mod.F90, so the coefficients and alpha ceiling are reproducible: - K=5 variable-timestep coefficients: derive_k5_variable_timesteps.py and generate_k5_fortran_lookup.py produce the 32-pattern XLBO_K5_C0..C5 / XLBO_K5_dK lookup tables. - K=10 coefficients: Niklasson_JCP_2009_table_I/ derive_xlbo_table_I_normalizations.py reproduces Table I (K=3-9) and extends to K=10 (verified: outputs C0_K10..C10_K10 = -858, 2652, -3094, 1496, 272, -952, 731, -322, 88, -14, 1; d_K=286). - Stability / alpha determination: empirical_stability.py (quasi-random history amplitude growth), stability_faithful.py (physical-band recurrence), reversibility_test.py (constant vs per-pattern kappa), find_max_alpha_per_pattern.py (per-pattern stability caps). Plus the READMEs and generated coefficient tables. Only the authoritative scripts are included; exploratory/dead-end variants are not committed. Co-Authored-By: Claude Opus 4.8 <[email protected]>
compute_alpha_table.py used the earlier alpha formula (0.018*3.0/d_K capped at alpha_max/2), which was replaced by the history-independent rule min(0.054*dt/d_K, 0.15) now active in prg_xlbo_mod.F90. Nothing references it: the adaptive-time-step doc documents the final rule and cites empirical_stability.py / find_max_alpha_per_pattern.py for the ceiling. Its only non-derivable content was a stale alpha_max snapshot, whose live source (find_max_alpha_per_pattern.py) is already committed. Co-Authored-By: Claude Opus 4.8 <[email protected]>
Group the three top-level K5 derivation files (derive_k5_variable_timesteps.py, generate_k5_fortran_lookup.py, and their k5_variable_timestep_coefficients.txt output) with the rest of the coefficient/stability scripts under the Niklasson subdirectory. Sibling import and output paths are unaffected. Co-Authored-By: Claude Opus 4.8 <[email protected]>
The AdaptiveTimeStep split criterion computed the projected max displacement as maxval(user_timestep*sy%velocity) -- the maximum over signed velocity components. That understates the displacement for atoms whose fastest motion is in a negative direction (-x/-y/-z), so genuinely fast atoms could fail to trigger a split. Use user_timestep*maxval(abs(sy%velocity)) so the test is on the velocity magnitude, consistent with the MAXVAL(ABS(...)) already used for the "Maximum Velocity" diagnostic just below. Co-Authored-By: Claude Opus 4.8 <[email protected]>
mewall
marked this pull request as draft
August 25, 2026 21:35
Reconstructs the performance/optimization history from ~2024-04 (pre-hackathon CPU vectorization campaign) through HEAD, organized by theme: CPU vectorization, GPU offload (OpenACC), neighbor-list rewrite + GPU residency, graph-partition / subgraph-update MPI reduction, memory/single-precision, core numerics, and the adaptive-timestep/XLBO scientific thread. Includes directive-level technique extracted from commit diffs and an honest list of benchmarks to regenerate (the history proves the changes but rarely their magnitudes). Co-Authored-By: Claude Opus 4.8 <[email protected]>
Add a "Cross-cutting optimization principles" section highlighting three themes the reviewer flagged as under-surfaced, each traced to concrete commits: (P1) minimize allocation cost via buffer reuse/resize in place (9033070, a20cb00, 7e8e21b, 599e723, 285a541); (P2) prioritize whole-kernel GPU offload to keep data resident and cut host<->device traffic (60f0ad3, baec46c, 9989a03); (P3) tune OpenACC pragmas -- loop level, collapse(2), gang/vector sizing, drop redundant worker/private clauses (c692490, c052c49, e890f53, a05ee53). Cross-reference the principles from the per-subsystem catalog and the framing section. Co-Authored-By: Claude Opus 4.8 <[email protected]>
Extend the cross-cutting principles section from three to seven: - P4 minimize MPI communication (low-comm subgraph graph update 46a1396, 6847b7e, d828f83; multiple parts per rank 2fd9a86; multi-rank consistency prereqs 67a98ce, b044b16, 32be3f8) - P5 more efficient algorithms (rankN update d44b400/0b49975, mu bisection 2e8e70d/2294304, new subgraph method be23670/6802493/f8474d0, partitioners 9391d45) - P6 increase accuracy to widen the stable time-step window (factor-of-2 DM fix 3b49c1a, SCF/kernel fix 91d6a57, FL convergence 6e96b54) - P7 adaptive time step to overcome longer-uniform-step instabilities (f0df335, 1576e3d, c318aa2, 296f24a) Wire P6->P7 (accuracy enables longer steps) and update the framing spine to reference all seven. Co-Authored-By: Claude Opus 4.8 <[email protected]>
Extract per-NVTX-region timings from the SEDACS baseline-vs-offload profile pairs (.nsys-rep -> exported .sqlite, read natively with sqlite3 -- no nsys binary needed) and add a "Measured performance" section to the notes: prg_get_charges 818.7 -> 65.6 ms (12.5x, 9989a03) prg_get_pulayforce 631.9 -> 135.4 ms (4.7x, 2851bb4/4a94d6a) get_skforce 415.8 -> 95.4 ms (4.4x, a6c45c5) DM build/diag DM_min 1578 -> 891 ms (1.77x, 00408fe/b25f46d) Also record the honest Ewald scale crossover (loss on small per-rank water-2088, win on larger TrpCage) as concrete evidence for principle P2. Add examples/gpmdk/tools/nvtx_regions.sh to regenerate any row. Softens the "magnitudes missing" caveat now that real numbers exist. Co-Authored-By: Claude Opus 4.8 <[email protected]>
Add examples/gpmdk/docs/profile_log.csv (the 58-run curated Nsight log with per-step + per-NVTX-tag means, STDEVs, baseline/after comparison sets, and the qmd-progress commit per run) and rewrite the "Measured performance" section around it. This log is authoritative and also resolves the previously-missing profiled commit SHAs (504b7d5, f41917c, 77afcce, d1d45d0, dc22aa2, ...). Highlights now cited from the log (NVTX-tag ms per MD step, before -> after): - P2 offload: charges 81.9->6.6 (12.4x), pulay 123->27 (4.5x), skforce 83->19 (4.4x), nonortho 108->45 (2.4x), DM_min 277->178, buildzdiag 273->202, Ewald TrpCage 517->414 - Honest scale-dependence: Ewald loses on small water-2088 (9.3->15.7) but wins on TrpCage; initial get_hsmat offload was slower (111->2040) before tuning - P3 pragma tuning cumulative on get_hsmat/InitParts: 1845->197 ms (9.3x) - P5 algorithms: mu-search NR removal 219->0.83 ms (algorithm removal), get_dH_or_dS inplace 279->98 (2.9x) Keep nvtx_regions.sh as an independent cross-check that reproduces the ratios. Co-Authored-By: Claude Opus 4.8 <[email protected]>
Instrument the MD coord/velocity reduction (halfVerlet -> updatecoords ->
6x prg_sumRealReduceN) in gpmdcov_MDloop with per-rank work/wait/reduce
timers to diagnose intermittent MPI_Allreduce slowness.
A diagnostic prg_barrierParallel separates the two failure modes: a
straggler shows large "work" on one rank and large "wait" on the rest,
while a genuinely slow collective shows uniformly large "reduce". Timings
are gathered to rank 0 via allGatherRealParallel and printed as a per-rank
vector plus min/max/mean and imbalance spreads.
Gated on a new opt-in input flag GPMD{ RankTiming=T } (default off).
Co-Authored-By: Claude Opus 4.8 <[email protected]>
… dumps
Extends GPMD{ RankTiming=T } to distinguish the two candidate causes of the
intermittent every-50-step slowness on the 512-rank run:
- Upstream per-phase per-rank timers (Part / InitParts / DM_min SCF /
EnergAndForces) gathered via allGatherRealParallel and reported by
gpmdcov_report_phase_timing. Reveals which phase a straggler burns time in,
since the coord/vel-reduction "wait" already showed the cost is arrival skew,
not the collective.
- Rank-0 file-I/O timers around the trajectory write and restart dump, printed
as "[RankTiming] fileio ..." to test the ~100s filesystem-stall hypothesis.
- Suppress restart dumps during annealing: print_mdstep is pinned at 0 until
mdstep > minimization_steps, so mod(0,DumpEach)==0 rewrote restart.dmp every
anneal step. Now gated on mdstep > minimization_steps.
- Fix straggler identification in gpmdcov_report_rank_timing: the last-arriving
rank has the MINIMUM wait (everyone waited on it), not the max work (work~0).
Co-Authored-By: Claude Opus 4.8 <[email protected]>
… energy best
Adds three opt-in integrator variants in the gpmdk MD loop, each gated
behind a GPMD{} flag and off by default, to test whether multirate
schemes can beat plain single-rate XLBO on adaptive-timestep energy drift:
- Respa=T / RespaInnerSteps=n : reversible r-RESPA, electronic step on a
uniform outer dt, nuclei substepped with the cheap rho-free pair force.
- UniformElectronicDt=T : full-force nuclear substeps with the XLBO
electronic propagation frozen to a uniform grid on split halves.
- RespaShadowCoul=T : RESPA with the Coulomb force refreshed each
substep from the frozen XLBO shadow charges (Ewald, no DM solve).
Best energy conservation overall: plain single-rate XLBO at a FIXED timestep
(AdaptiveTimeStep=F, or Respa=T with RespaInnerSteps=1, which is
bit-identical to single-rate).
Best among the adaptive / variable-timestep schemes: the existing plain
AdaptiveTimeStep split-step scheme. It stays near the energy noise floor
for small systems (300-atom water) and drifts appreciably only for large
systems (the per-split-per-atom injection scales with N). All three
multirate variants added here drift even for small systems, so none of
them is preferred; they are retained only as validated references and are
documented as such in the code.
Also adds verbose>=2 force-split diagnostics (FORCECOMP per-component
force max/rms, SPLITDIAG per-step displacement and split flag) used to
diagnose the drift. No src/ library changes.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
mewall
marked this pull request as ready for review
September 8, 2026 21:08
mewall
force-pushed
the
split_step_clean
branch
from
September 8, 2026 21:39
1e34d80 to
0efacc7
Compare
Adds a "Multirate follow-on schemes" section covering the three opt-in
integrator variants now in the code (Respa/RespaInnerSteps,
UniformElectronicDt, RespaShadowCoul), each with its GPMD{} invocation,
mutual-exclusivity rules, and the measured energy result. Records that
none improved on single-rate: fixed-dt single-rate conserves energy best
overall, and plain AdaptiveTimeStep is best among variable-timestep
schemes. Results are stated with their measurement windows; mechanism
explanations are marked as interpretation tied to the FORCECOMP force
magnitudes, not asserted as proven.
Also: correct the adaptive-trigger snippet to match the source
(user_timestep * maxval(abs(sy%velocity))), note the multirate route
under the alpha "structural levers" discussion, and mark the reversible-
backtracking section as a design note not implemented on this branch.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Co-Authored-By: Claude Opus 4.8 <[email protected]>
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.
Adaptive time step studies.
The default one with AdaptiveTimeStep=True appears to exhibit the best stability.
In this default mode, the time step is split whenever the maximum distance an atom moves along any direction is too high.
No noticeable drift for 300 atom systems, but drift is seen for ~30,000 atom systems. It might be a reasonable choice for NVE using small systems, or NVT (Langevin) simulations of larger systems.
Also includes documentation of performance improvements since the Venado hackathon.