Energy drift in split_step branch due to XLBO multistep predictor assuming uniform timesteps

[mewall]

Description

The split_step branch implements adaptive timestepping for MD simulations, but exhibits systematic energy drift when timestep splitting occurs. Energy
conservation is correct when no splits happen, confirming the issue is specific to the splitting mechanism.

Root Cause

The Extended Lagrangian Born-Oppenheimer (XLBO) integration uses a 6th-order multistep predictor with fixed coefficients that assume uniform timestep
spacing:

File: src/prg_xlbo_mod.F90, lines 18-23, 149-150

real(dp), parameter :: C0 = -6.0_dp
real(dp), parameter :: C1 = 14.0_dp
real(dp), parameter :: C2 = -8.0_dp
real(dp), parameter :: C3 = -3.0_dp
real(dp), parameter :: C4 = -4.0_dp
real(dp), parameter :: C5 = -1.0_dp

n = 2.0_dp*n_0 - n_1 + xl%cc*kappa*(charges-n) &
     + alpha*(C0*n_0+C1*n_1+C2*n_2+C3*n_3+C4*n_4+C5*n_5)

These coefficients are derived assuming the history values n_0, n_1, n_2, n_3, n_4, n_5 are at uniformly-spaced time intervals: t, t-dt, t-2dt, t-3dt, t-4dt, 
t-5dt.

When splitting occurs, the actual spacing becomes non-uniform:
Step 10: dt/2 (first half)  → n_0 at t=0
Step 11: dt/2 (second half) → n_1 at t=-dt/2
Step 12: dt (full step)     → n_2 at t=-dt
                              n_3 at t=-2dt
                              n_4 at t=-3dt
                              n_5 at t=-4dt
                              
The fixed coefficients are incorrect for non-uniform spacing, causing the electronic degrees of freedom to drift from the correct trajectory. This error
propagates into the forces and causes atomic energy drift.

Affected Code

The issue affects all XLBO integration routines:
- prg_xlbo_nint() - line 118
- prg_xlbo_nint_kernel() - line 158
- prg_xlbo_nint_kernelTimesRes() - line 210

Called from examples/gpmdk/src/gpmdcov_mdloop.F90:
- Lines 311, 319, 375, 381, 389

Proposed Solutions

Option 1: Reset History After Split Pairs (Quick Fix)

Pros: Simple 5-line change, confirms diagnosis
Cons: Temporarily reduces accuracy during splits

After completing a split pair in gpmdcov_mdloop.F90:
if (.not. first_substep_taken .and. half_timestep_flag) then
    ! Reset XLBO history to prevent drift from non-uniform spacing
    n_1 = n_0
    n_2 = n_0
    n_3 = n_0
    n_4 = n_0
    n_5 = n_0
endif

Option 2: Dynamic Coefficient Computation (Recommended)

Pros: Maintains full accuracy, mathematically correct
Cons: More complex implementation

1. Add timestep history tracking to xlbo_type:
type, public :: xlbo_type
    ! ... existing fields ...
    real(dp) :: dt_history(5)  ! Last 5 timesteps
    integer :: nsteps_taken
end type

2. Compute coefficients dynamically based on actual time intervals using backward finite differences on non-uniform grids
3. Update prg_xlbo_nint() signature to accept current timestep:
subroutine prg_xlbo_nint(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,xl,dt)
    real(dp), intent(in) :: dt
    
4. Compute adaptive coefficients each iteration based on xl%dt_history

Option 3: Fallback to Lower-Order Method

Pros: Moderate complexity, maintains stability
Cons: Reduced accuracy during splits

Detect non-uniform spacing and fall back to 2nd-order predictor:
if (history_is_uniform) then
    ! Use full 6th order formula
else
    ! Use simpler 2nd order: n = 2.0*n_0 - n_1 + xl%cc*kappa*(charges-n)
endif

Testing Plan

1. Verify Option 1 eliminates drift (confirms diagnosis)
2. Implement Option 2 for production use
3. Compare energy conservation with/without splits
4. Validate against systems with frequent velocity-triggered splits

References

- Split step implementation: remotes/origin/split_step
- Main MD loop: examples/gpmdk/src/gpmdcov_mdloop.F90
- XLBO module: src/prg_xlbo_mod.F90

Related Commits

Recent work on split_step branch:
- 47ba959 - Fix timestep for .pdb trajectory output
- be1897d - Make adaptive time step compatible with minimization
- 86b4c9e - Added lines for substeps (si) in gpmdcov_mdloop