Adaptive time step in gpmdk - #374
Draft
mewall wants to merge 50 commits into
Draft
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
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 <noreply@anthropic.com>
…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 <noreply@anthropic.com>
Implements direct variable timestep coefficients for K=5 XLBO dissipation as an alternative to interpolation. Uses lookup tables with all 32 possible 5-step timestep history patterns (dt or dt/2). Key features: - Derives coefficients for each timestep history using moment constraints - Scales alpha by 3.0/d_K to maintain proper dissipation strength - Automatically enabled when AdaptiveTimeStep=T and ExtendedHistory=F - Successfully completes 1000 steps with excellent energy conservation Energy conservation results (1000 steps, TimeStep=0.4): - Total drift: -0.025 eV - Drift per step: -0.025 meV/step - Relative drift: 0.0018% - Std deviation: 0.017 eV Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Bug: When using variable timestep coefficients (use_variable_coeffs=T), the code loaded C0_use through C3_use from pattern-specific lookup tables, but C4 and C5 were missing from the lookup. This caused the integration formula to use the fixed constants C4=4.0 and C5=-1.0 instead of the correct pattern-specific values. This bug particularly affected patterns where positions n-3 and n-4 (multiplied by c_4 and c_5) contained timesteps that differed from the uniform case, leading to incorrect dynamics and instability. Fix: 1. Added C4_use and C5_use to coefficient lookup (lines 820-821) 2. Updated integration formula to use C4_use and C5_use (line 827) 3. Regenerated all K=5 coefficient arrays using Method 3 (min-norm rescaled to c_5=-1) which provides correct pattern-specific values for all six coefficients including C4 and C5 Note: This fix is necessary but not sufficient for full stability. The variable coefficient approach still requires pattern-specific alpha and kappa values (not yet implemented), and 56% of patterns remain fundamentally unstable. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
- Replace fixed Verlet coefficients (2, -1) with variable coefficients (1+r, -r) - Use proper scaling factor: 0.5*dt_n*(dt_n + dt_prev) for both kappa and alpha - Remove separate dt^2 scaling of kappa (now in unified scaling factor) - Fixes all three integration subroutines: nint, nint_kernel, nint_kernelTimesRes
Replace d_K-based alpha scaling with direct lookup of pattern-specific alpha values that have been verified stable for all 32 timestep patterns. Hybrid strategy: - Use α_scaled = 0.018 × 3.0/d_K when safe (23/32 patterns) - Use α_max/2 when scaling would exceed stability limits (9/32 patterns) All 32 patterns verified stable with max|λ| < 1.0 over γ ∈ [-1, 1] Safety margins range from 1.15× to 6.11× above actual usage Test results (1000 steps, TimeStep=0.4 fs): - Completed successfully: 1000/1000 steps - Energy slope: -2.14 × 10⁻⁵ eV/step (excellent) - Split-steps triggered: 0 Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Replace pattern-specific alpha values with conservative uniform scaling: alpha = 0.000796 × 3.0 / d_K This ensures alpha ≤ α_max/2 for all 32 patterns, providing a 20× average safety margin (vs 2.2× before). The more conservative values should prevent XLBO instability on large systems (2088+ atoms) while maintaining excellent energy conservation on small systems. Limiting pattern: 6 (d_K=0.015625, exactly at α_max/2 limit) Test results (696 atoms, 1000 steps, TimeStep=0.4 fs): - Completed successfully: 1000/1000 steps - Energy slope: -2.78 × 10⁻⁵ eV/step (excellent) - Previous slope: -2.14 × 10⁻⁵ eV/step - Trade-off: ~30% slower dissipation for much better stability Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
The stability analysis computed pattern-specific alpha values assuming alpha is NOT scaled by timestep factors. The code was incorrectly multiplying alpha by kappa_alpha_scale, which made alpha approximately constant across patterns instead of pattern-specific. This caused instability on large systems (2088 atoms) because patterns with small d_K were getting too little dissipation (alpha_effective was approximately constant instead of alpha ~ 1/d_K as intended). Fixed by removing kappa_alpha_scale from all 15 alpha term locations in: - prg_xlbo_nint (5 locations) - prg_xlbo_nint_kernel (5 locations) - prg_xlbo_nint_kernelTimesRes (5 locations) Kappa term correctly keeps kappa_alpha_scale. Alpha now varies properly with pattern as designed by the stability analysis. Tested on 696-atom water system: energy slope -4.94e-05 eV/step over 693 steps, resnorm stable at ~5e-06. Ready for large system testing. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Changes: 1. Revert alpha calculation to d_K scaling method: alpha_use = alpha * XLBO_K5_dK(31) / XLBO_K5_dK(hist_idx) Pattern-specific alpha table proved too aggressive for large systems. 2. Fix C4/C5 coefficients to use pattern-specific values: Use XLBO_K5_C4(hist_idx) and XLBO_K5_C5(hist_idx) instead of fixed C4/C5 values. This ensures all coefficients are consistent with the timestep pattern. 3. Rename variable for clarity: use_interpolation -> allow_adaptive_timestep Better reflects that this controls adaptive timestep behavior, not just interpolation method. Tested on 696-atom water system: energy slope -1.03e-04 eV/step, resnorm stable at ~1-3e-06. Ready for large system testing. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Critical bug: Alpha was only scaled by d_K after step 6 when allow_adaptive_timestep became true. During the warmup phase (steps 1-6), alpha used the base value without d_K scaling. This caused insufficient dissipation during early half-steps (d_K ~ 0.75) where alpha should be ~4× larger than for full steps (d_K ~ 3.0). The result was XLBO instability when taking the first full step transition. Fix: Apply d_K scaling to alpha immediately during warmup phase, before pattern-specific coefficients are available. This ensures proper dissipation from step 1 onwards. Applied to all three integration routines: - prg_xlbo_nint - prg_xlbo_nint_kernel - prg_xlbo_nint_kernelTimesRes Tested on 696-atom water system: resnorm stable at ~2-9e-06, energy slope -2.19e-04 eV/step. Ready for large system testing. Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
…tory(1) Critical bug: Verlet coefficients were using xl%dt_history(1) and xl%dt_history(2) as dt_n and dt_prev, which are the PREVIOUS two timesteps, not the CURRENT and PREVIOUS timesteps. This caused incorrect Verlet extrapolation, especially at pattern transitions like step 22 (first full step after half steps). Correct physics: - dt_n = dt (input parameter - CURRENT timestep we're taking) - dt_prev = xl%dt_history(1) (PREVIOUS timestep) - XLBO coefficients use xl%dt_history(1:5) for charge history pattern Changed in all three routines: dt_n = xl%dt_history(1) --> dt_n = dt dt_prev = xl%dt_history(2) --> dt_prev = xl%dt_history(1) This ensures: - Verlet P_n_coeff, P_n1_coeff computed with correct timestep ratio - kappa_alpha_scale = 0.5 * dt_n * (dt_n + dt_prev) uses correct values - XLBO dissipation coefficients still match charge history correctly Tested on 696-atom water: resnorm stable ~1-3e-06, energy slope -1.88e-04 eV/step, step 22 transition stable (resnorm 2.35e-06). Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>
Root cause: XLBO-integrated charges (n, n_0) were not synchronized across MPI ranks after integration, causing small numerical differences to accumulate in the charge history arrays. These differences amplified during timestep transitions, leading to instability on multi-rank runs. The fix adds MPI allreduce synchronization after every XLBO integration call: - prg_xlbo_nint - prg_xlbo_nint_kernel - prg_xlbo_nint_kernelTimesRes Both n and n_0 are synchronized since n_0=n is done inside the XLBO routine before returning. This ensures all ranks maintain identical charge history. Tested: Single rank always worked; multi-rank instability diagnosed via comparison showing problem disappeared on single rank.
Add rank ID to debug prints to diagnose MPI divergence: - Show which rank is making timestep splitting decisions - Show maxdisp calculation per rank - Confirm XLBO charge synchronization is happening This will help identify if different ranks are making different decisions or if the instability has another cause.
CRITICAL BUG: print_mdstep was only incremented on rank 1, causing it to stay at 0 on all other ranks. Since print_mdstep is used to decide when to stop forced timestep splitting (line 161), this caused: - Rank 1: stops forcing splits after print_mdstep > 4 (or 6 for K=10) - Other ranks: keep forcing splits forever (print_mdstep always 0) This caused different ranks to use different timesteps, leading to divergence and instability during the transition from forced to adaptive timesteps. The fix moves the print_mdstep increment outside the myRank==1 block so all ranks maintain the same counter value.
Replace dynamic alpha scaling (alpha * 3.0 / d_K) with pre-computed table values. The table is computed using the same formula but capped at alpha_max/2 for guaranteed stability. Benefits: - More conservative: 23 out of 32 patterns are capped below formula value - Guaranteed stability: all values <= alpha_max/2 - Same dissipation for uncapped patterns New alpha table computed as: alpha = min(0.018 * 3.0 / d_K, alpha_max/2) Add compute_alpha_table.py script to regenerate table if needed.
mewall
marked this pull request as draft
July 11, 2026 17:32
cnegre
approved these changes
Jul 11, 2026
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
3 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.
Add adaptive time step option in gpmdk
o Supports 5-timestep history for xlbo dissipative term