diff --git a/build.sh b/build.sh
index 2bea0dda..1e7217e3 100755
--- a/build.sh
+++ b/build.sh
@@ -1,6 +1,6 @@
 #!/bin/bash
 
-TOP_DIR=$(readlink --canonicalize $(dirname $0))
+TOP_DIR=$(readlink -f $(dirname $0))
 
 : ${BUILD_DIR:=${TOP_DIR}/build}
 : ${INSTALL_DIR:=${TOP_DIR}/install}
diff --git a/build_scaling.sh b/build_scaling.sh
index e2e1a9cf..070c8a2f 100755
--- a/build_scaling.sh
+++ b/build_scaling.sh
@@ -15,11 +15,18 @@ else
     export CXX=${CXX:=g++}
 fi
 
+#export FFLAGS="-I$CONDA_PREFIX/lib"
+#export FCFLAGS="-I$CONDA_PREFIX/lib"
+#export LD_LIBRARY_FLAGS=$LD_LIBRARY_FLAGS:"$CONDA_PREFIX/lib"
+
+export PROGRESS_MPI=no
+export BML_DIR=${BML_DIR:=/Users/mewall/packages/gpmd/bml/install}
 export PROGRESS_OPENMP=${PROGRESS_OPENMP:=yes}
-export PROGRESS_GRAPHLIB=${PROGRESS_GRAPHLIB:=yes}
+export PROGRESS_GRAPHLIB=${PROGRESS_GRAPHLIB:=no}
 export PROGRESS_TESTING=${PROGRESS_TESTING:=yes}
 export PROGRESS_EXAMPLES=${PROGRESS_EXAMPLES:=yes}
 export CMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE:=Release}
+export CMAKE_PREFIX_PATH=${CMAKE_PREFIX_PATH}:${BML_DIR}
 export VERBOSE_MAKEFILE=${VERBOSE_MAKEFILE:=yes}
 export COMMAND=${1:-compile}
 
diff --git a/examples/gpmdk/compare_branches.sh b/examples/gpmdk/compare_branches.sh
new file mode 100755
index 00000000..edc6c437
--- /dev/null
+++ b/examples/gpmdk/compare_branches.sh
@@ -0,0 +1,311 @@
+#!/bin/bash
+#
+# Compare MD simulation behavior between two branches
+#
+# Usage: ./compare_branches.sh <branch1> <branch2> <mdsteps> <timestep> [run_dir]
+#
+# Example: ./compare_branches.sh split_step xlbo_adapt 50 0.4 run/water
+#
+
+set -e
+
+# Check arguments
+if [ "$#" -lt 4 ] || [ "$#" -gt 5 ]; then
+    echo "Usage: $0 <branch1> <branch2> <mdsteps> <timestep> [run_dir]"
+    echo "Example: $0 split_step xlbo_adapt 50 0.4 run/water"
+    echo ""
+    echo "Arguments:"
+    echo "  branch1, branch2: Git branch names to compare"
+    echo "  mdsteps:          Number of MD steps to run"
+    echo "  timestep:         Timestep in femtoseconds"
+    echo "  run_dir:          Directory containing input.in (default: run/water)"
+    exit 1
+fi
+
+BRANCH1="$1"
+BRANCH2="$2"
+MDSTEPS="$3"
+TIMESTEP="$4"
+RUN_SUBDIR="${5:-run/water}"
+
+# Script is in examples/gpmdk/, so repo root is ../..
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+REPO_ROOT="$(cd "${SCRIPT_DIR}/../.." && pwd)"
+BUILD_DIR="${REPO_ROOT}/build"
+RUN_DIR="${SCRIPT_DIR}/${RUN_SUBDIR}"
+INPUT_FILE="${RUN_DIR}/input.in"
+
+# Check that run directory exists
+if [ ! -d "${RUN_DIR}" ]; then
+    echo "ERROR: Run directory not found: ${RUN_DIR}"
+    exit 1
+fi
+
+if [ ! -f "${INPUT_FILE}" ]; then
+    echo "ERROR: input.in not found in: ${RUN_DIR}"
+    exit 1
+fi
+
+# Save original branch
+ORIGINAL_BRANCH=$(git rev-parse --abbrev-ref HEAD)
+
+echo "=========================================="
+echo "   Branch Comparison Tool"
+echo "=========================================="
+echo "Branch 1:   ${BRANCH1}"
+echo "Branch 2:   ${BRANCH2}"
+echo "MD Steps:   ${MDSTEPS}"
+echo "TimeStep:   ${TIMESTEP} fs"
+echo "Run Dir:    ${RUN_SUBDIR}"
+echo ""
+
+# Function to run simulation on a branch
+run_branch() {
+    local BRANCH="$1"
+    local OUTPUT_FILE="$2"
+
+    echo "=========================================="
+    echo "Running ${BRANCH}"
+    echo "=========================================="
+
+    # Checkout branch
+    echo "Checking out ${BRANCH}..."
+    cd "${REPO_ROOT}"
+    git checkout "${BRANCH}" 2>&1 | grep -v "^M\s" || true
+
+    # Build
+    echo "Building ${BRANCH}..."
+    cd "${BUILD_DIR}"
+    make -j4 install > /dev/null 2>&1 || {
+        echo "ERROR: Build failed on ${BRANCH}"
+        cd "${REPO_ROOT}"
+        git checkout "${ORIGINAL_BRANCH}" 2>&1 | grep -v "^M\s" || true
+        exit 1
+    }
+
+    # Modify input.in
+    echo "Setting TimeStep=${TIMESTEP} and MDSteps=${MDSTEPS} in input.in..."
+    cd "${RUN_DIR}"
+
+    # Backup original input.in if not already backed up
+    if [ ! -f input.in.backup ]; then
+        cp input.in input.in.backup
+    fi
+
+    # Restore from backup and modify
+    cp input.in.backup input.in
+    sed -i.tmp "s/TimeStep=.*/TimeStep= ${TIMESTEP}/" input.in
+    sed -i.tmp "s/MDSteps=.*/MDSteps= ${MDSTEPS}/" input.in
+    rm -f input.in.tmp
+
+    # Run simulation
+    echo "Running simulation on ${BRANCH}..."
+    OMP_NUM_THREADS=4 "${BUILD_DIR}/gpmdk" input.in > "${OUTPUT_FILE}" 2>&1 || {
+        echo "ERROR: Simulation failed on ${BRANCH}"
+        cd "${REPO_ROOT}"
+        git checkout "${ORIGINAL_BRANCH}" 2>&1 | grep -v "^M\s" || true
+        exit 1
+    }
+
+    echo "${BRANCH} complete!"
+    echo ""
+}
+
+# Run both branches
+cd "${REPO_ROOT}"
+OUTPUT1="${RUN_DIR}/out_${BRANCH1}_comparison"
+OUTPUT2="${RUN_DIR}/out_${BRANCH2}_comparison"
+
+run_branch "${BRANCH1}" "${OUTPUT1}"
+run_branch "${BRANCH2}" "${OUTPUT2}"
+
+# Return to original branch
+echo "Returning to ${ORIGINAL_BRANCH}..."
+cd "${REPO_ROOT}"
+git checkout "${ORIGINAL_BRANCH}" 2>&1 | grep -v "^M\s" || true
+
+# Restore original input.in
+if [ -f "${INPUT_FILE}.backup" ]; then
+    cp "${INPUT_FILE}.backup" "${INPUT_FILE}"
+fi
+
+# Analyze results
+echo "=========================================="
+echo "   Analyzing Results"
+echo "=========================================="
+echo ""
+
+cd "${RUN_DIR}"
+
+# Export variables for Python
+export BRANCH1="${BRANCH1}"
+export BRANCH2="${BRANCH2}"
+export TIMESTEP="${TIMESTEP}"
+export MDSTEPS="${MDSTEPS}"
+
+python << PYTHON_ANALYSIS
+import numpy as np
+import sys
+import os
+
+def analyze_with_splits(filename, label):
+    """Analyze energy data and split-steps from output file"""
+    energies = []
+    split_steps = []
+
+    if not os.path.exists(filename):
+        print(f"ERROR: Output file not found: {filename}")
+        sys.exit(1)
+
+    with open(filename, 'r') as f:
+        for line in f:
+            if line.startswith("Mdstep, Energy"):
+                parts = line.split()
+                step = int(parts[5])
+                energy = float(parts[6])
+                energies.append((step, energy))
+            elif "Splitting mdstep" in line:
+                split_step = int(line.split()[-1])
+                split_steps.append(split_step)
+
+    if len(energies) == 0:
+        print(f"ERROR: No energy data found in {filename}")
+        sys.exit(1)
+
+    energy_arr = np.array([e[1] for e in energies])
+
+    e_mean = np.mean(energy_arr)
+    e_std = np.std(energy_arr)
+    e_min = np.min(energy_arr)
+    e_max = np.max(energy_arr)
+    e_drift = e_max - e_min
+    e_drift_pct = (e_drift / abs(e_mean)) * 100
+
+    timesteps = np.arange(len(energy_arr))
+    coeffs = np.polyfit(timesteps, energy_arr, 1)
+    slope = coeffs[0]
+
+    return {
+        'label': label,
+        'steps': len(energy_arr),
+        'mean': e_mean,
+        'std': e_std,
+        'drift': e_drift,
+        'drift_pct': e_drift_pct,
+        'slope': slope,
+        'energies': energy_arr,
+        'split_count': len(split_steps),
+        'split_range': (min(split_steps), max(split_steps)) if split_steps else (None, None)
+    }
+
+# Get branch names from environment
+branch1 = os.environ.get('BRANCH1', 'branch1')
+branch2 = os.environ.get('BRANCH2', 'branch2')
+timestep = os.environ.get('TIMESTEP', '?')
+mdsteps = os.environ.get('MDSTEPS', '?')
+
+# Analyze both outputs
+data1 = analyze_with_splits(f'out_{branch1}_comparison', branch1)
+data2 = analyze_with_splits(f'out_{branch2}_comparison', branch2)
+
+# Print comparison table
+print("="*75)
+print(f"   COMPARISON at TimeStep={timestep} fs, MDSteps={mdsteps}")
+print("="*75)
+print("")
+print(f"{'Metric':<35} {branch1:>18} {branch2:>18}")
+print("-"*75)
+print(f"{'Total MD Steps':<35} {data1['steps']:>18} {data2['steps']:>18}")
+print(f"{'Split-steps triggered':<35} {data1['split_count']:>18} {data2['split_count']:>18}")
+if data1['split_count'] > 0:
+    range1 = f"{data1['split_range'][0]}-{data1['split_range'][1]}"
+    range2 = f"{data2['split_range'][0]}-{data2['split_range'][1]}"
+    print(f"{'Split-step range':<35} {range1:>18} {range2:>18}")
+print("")
+print(f"{'Mean Energy (eV)':<35} {data1['mean']:>18.6f} {data2['mean']:>18.6f}")
+print(f"{'Std Dev (eV)':<35} {data1['std']:>18.6f} {data2['std']:>18.6f}")
+print(f"{'Total Drift (eV)':<35} {data1['drift']:>18.6f} {data2['drift']:>18.6f}")
+print(f"{'Drift (% of E)':<35} {data1['drift_pct']:>18.4f} {data2['drift_pct']:>18.4f}")
+print(f"{'Linear Drift (eV/step)':<35} {data1['slope']:>18.6e} {data2['slope']:>18.6e}")
+print("")
+
+# Calculate differences
+max_diff = np.max(np.abs(data1['energies'] - data2['energies']))
+mean_diff = np.mean(np.abs(data1['energies'] - data2['energies']))
+print(f"{'Maximum energy difference':<35} {max_diff:>18.6e} eV")
+print(f"{'Mean absolute difference':<35} {mean_diff:>18.6e} eV")
+print("")
+
+# Conclusion
+print("="*75)
+print("                            CONCLUSION")
+print("="*75)
+print("")
+
+if data1['split_count'] > 0:
+    print(f"Split-steps triggered: {data1['split_count']} times", end="")
+    if data1['split_range'][0]:
+        print(f" (steps {data1['split_range'][0]}-{data1['split_range'][1]})")
+    else:
+        print()
+    print("")
+
+if max_diff < 1e-10:
+    print("✓ Both methods produce IDENTICAL results")
+elif max_diff < 1e-6:
+    print("✓ Both methods produce essentially identical results")
+    print(f"  (max difference {max_diff:.2e} eV - likely numerical noise)")
+elif max_diff < 0.001:
+    print(f"~ Methods show small differences (max {max_diff:.5f} eV)")
+    print(f"  Mean difference: {mean_diff:.2e} eV")
+else:
+    print(f"⚠️  Methods show measurable differences:")
+    print(f"  Max difference:  {max_diff:.5f} eV")
+    print(f"  Mean difference: {mean_diff:.5f} eV")
+
+print("")
+print(f"Both branches maintain {'excellent' if max(data1['drift_pct'], data2['drift_pct']) < 0.01 else 'good'} energy conservation")
+print(f"({branch1}: {data1['drift_pct']:.4f}%, {branch2}: {data2['drift_pct']:.4f}%)")
+print("")
+
+# Save results to file
+output_file = f"comparison_{branch1}_vs_{branch2}_ts{timestep}_md{mdsteps}.txt"
+with open(output_file, 'w') as f:
+    f.write("="*75 + "\n")
+    f.write(f"   COMPARISON: {branch1} vs {branch2}\n")
+    f.write("="*75 + "\n")
+    f.write(f"TimeStep: {timestep} fs\n")
+    f.write(f"MDSteps:  {mdsteps}\n")
+    f.write(f"Date:     {os.popen('date').read().strip()}\n")
+    f.write("\n")
+    f.write(f"{'Metric':<35} {branch1:>18} {branch2:>18}\n")
+    f.write("-"*75 + "\n")
+    f.write(f"{'Total MD Steps':<35} {data1['steps']:>18} {data2['steps']:>18}\n")
+    f.write(f"{'Split-steps triggered':<35} {data1['split_count']:>18} {data2['split_count']:>18}\n")
+    if data1['split_count'] > 0:
+        range1 = f"{data1['split_range'][0]}-{data1['split_range'][1]}"
+        range2 = f"{data2['split_range'][0]}-{data2['split_range'][1]}"
+        f.write(f"{'Split-step range':<35} {range1:>18} {range2:>18}\n")
+    f.write("\n")
+    f.write(f"{'Mean Energy (eV)':<35} {data1['mean']:>18.6f} {data2['mean']:>18.6f}\n")
+    f.write(f"{'Std Dev (eV)':<35} {data1['std']:>18.6f} {data2['std']:>18.6f}\n")
+    f.write(f"{'Total Drift (eV)':<35} {data1['drift']:>18.6f} {data2['drift']:>18.6f}\n")
+    f.write(f"{'Drift (% of E)':<35} {data1['drift_pct']:>18.4f} {data2['drift_pct']:>18.4f}\n")
+    f.write(f"{'Linear Drift (eV/step)':<35} {data1['slope']:>18.6e} {data2['slope']:>18.6e}\n")
+    f.write("\n")
+    f.write(f"{'Maximum energy difference':<35} {max_diff:>18.6e} eV\n")
+    f.write(f"{'Mean absolute difference':<35} {mean_diff:>18.6e} eV\n")
+
+print(f"Results saved to: {output_file}")
+print("")
+PYTHON_ANALYSIS
+
+echo "=========================================="
+echo "   Comparison Complete"
+echo "=========================================="
+echo ""
+echo "Output files:"
+echo "  ${OUTPUT1}"
+echo "  ${OUTPUT2}"
+echo "  comparison_${BRANCH1}_vs_${BRANCH2}_ts${TIMESTEP}_md${MDSTEPS}.txt"
+echo ""
diff --git a/examples/gpmdk/compare_branches_README.md b/examples/gpmdk/compare_branches_README.md
new file mode 100644
index 00000000..3fb2d0a0
--- /dev/null
+++ b/examples/gpmdk/compare_branches_README.md
@@ -0,0 +1,142 @@
+# Branch Comparison Script
+
+## Overview
+
+The `compare_branches.sh` script automatically compares MD simulation behavior between two branches, building each branch, running simulations with specified parameters, and generating a detailed comparison report.
+
+## Location
+
+This script is located in `examples/gpmdk/` and works with any run directory under `examples/gpmdk/`.
+
+## Usage
+
+```bash
+./compare_branches.sh <branch1> <branch2> <mdsteps> <timestep> [run_dir]
+```
+
+**Arguments:**
+- `branch1`: First branch name (e.g., split_step)
+- `branch2`: Second branch name (e.g., xlbo_adapt)
+- `mdsteps`: Number of MD steps to run
+- `timestep`: Timestep in femtoseconds
+- `run_dir`: Directory containing input.in, relative to examples/gpmdk/ (default: run/water)
+
+**Examples:**
+```bash
+# Use default run/water directory
+./compare_branches.sh split_step xlbo_adapt 50 0.6
+
+# Specify a different run directory
+./compare_branches.sh split_step xlbo_adapt 100 0.35 run/ammonia
+
+# Run from examples/gpmdk/ directory
+cd examples/gpmdk
+./compare_branches.sh split_step xlbo_adapt 50 0.4
+```
+
+## What the Script Does
+
+1. **Checks out and builds each branch**
+   - Automatically switches between branches
+   - Builds each branch using `make -j4 install`
+   - Returns to original branch when complete
+
+2. **Modifies input parameters**
+   - Backs up original `input.in`
+   - Sets requested `TimeStep` and `MDSteps`
+   - Restores original after completion
+
+3. **Runs simulations**
+   - Executes GPMD with `OMP_NUM_THREADS=4`
+   - Captures full output for analysis
+
+4. **Analyzes and compares results**
+   - Extracts energy data from both runs
+   - Counts split-step occurrences
+   - Calculates energy conservation metrics
+   - Computes differences between branches
+
+5. **Generates reports**
+   - Prints comparison table to console
+   - Saves detailed results to file
+   - Preserves output files for inspection
+
+## Output Files
+
+The script generates files in the specified run directory:
+- `out_<branch1>_comparison` - Full simulation output for branch 1
+- `out_<branch2>_comparison` - Full simulation output for branch 2
+- `comparison_<branch1>_vs_<branch2>_ts<timestep>_md<mdsteps>.txt` - Detailed comparison report
+- `input.in.backup` - Backup of original input.in (created if needed)
+
+## Comparison Metrics
+
+The script reports:
+- **Total MD Steps**: Number of output steps produced
+- **Split-steps triggered**: How many times timestep was split
+- **Split-step range**: Which steps had splits
+- **Mean Energy**: Average total energy
+- **Std Dev**: Energy fluctuation magnitude
+- **Total Drift**: Maximum energy change
+- **Drift (% of E)**: Relative energy drift
+- **Linear Drift**: Systematic energy drift rate
+- **Maximum/Mean energy difference**: How much branches differ
+
+## Interpretation
+
+**Energy Conservation Quality:**
+- < 0.01% drift = Excellent
+- 0.01-0.1% drift = Good
+- \> 0.1% drift = Poor (investigate)
+
+**Branch Differences:**
+- < 1e-10 eV: Identical (within machine precision)
+- < 1e-6 eV: Essentially identical (numerical noise)
+- < 0.001 eV: Small differences
+- \> 0.001 eV: Measurable differences
+
+## Example Output
+
+```
+===========================================================================
+   COMPARISON at TimeStep=0.6 fs, MDSteps=50
+===========================================================================
+
+Metric                                      split_step         xlbo_adapt
+---------------------------------------------------------------------------
+Total MD Steps                                      50                 50
+Split-steps triggered                               96                 96
+Split-step range                                  3-98               3-98
+
+Mean Energy (eV)                          -1360.058794       -1360.060119
+Std Dev (eV)                                  0.009284           0.008485
+Total Drift (eV)                              0.045340           0.043950
+Drift (% of E)                                  0.0033             0.0032
+Linear Drift (eV/step)                    1.604643e-04       1.332754e-04
+
+Maximum energy difference                 3.390000e-03 eV
+Mean absolute difference                  1.535800e-03 eV
+```
+
+## Requirements
+
+- **Python**: Must have numpy installed
+- **Git**: Repository must be a git repo
+- **Build system**: CMake/Make setup must work
+- **Both branches**: Must exist and be buildable
+
+## Notes
+
+- Script automatically backs up and restores `input.in`
+- Returns to original branch on completion
+- Handles build failures gracefully
+- Safe to run multiple times
+- Uses conda/system Python (not python3)
+
+## Troubleshooting
+
+**"Build failed"**: Check that both branches compile successfully manually first
+
+**"No energy data found"**: Simulation may have crashed - check output files directly
+
+**"operands could not be broadcast"**: Branches produced different numbers of output steps - ensure both have MDSteps fix applied
diff --git a/examples/gpmdk/run/water/input.in b/examples/gpmdk/run/water/input.in
index 75c825a4..71c27b13 100644
--- a/examples/gpmdk/run/water/input.in
+++ b/examples/gpmdk/run/water/input.in
@@ -29,17 +29,17 @@ Latte{
   MaxSCFIter= 500
   CoulAcc= 1.0d-5
   TimeRatio= 10.0
-  TimeStep= 0.2
+  TimeStep= 0.6
   #TimeStep= 0.00
-  MDSteps= 20
+  MDSteps= 50
   #ParamPath= "../sulfurTBparam"
   ParamPath= "../../tests/latteTBparams"
   #ParamPath= "../latteTBparams_orig"
   #CoordsFile= coords.ltt
   #CoordsFile= coords_300New.dat
   #CoordsFile= coords_300_sort.dat
-  #CoordsFile= coords_300.dat
-  CoordsFile= coords_2088.dat
+  CoordsFile= coords_300.dat
+  #CoordsFile= coords_2088.dat
   #CoordsFile= "./polyaniline.pdb"
   #CoordsFile= graphite2048.pdb
   #CoordsFile= carbon_2197.pdb
@@ -99,16 +99,16 @@ GSP2{
   BMLType= Ellpack
   #GraphElement=  Orbital
   GraphElement=         Atom
-  #PartitionType=        Block
+  PartitionType=        Block
   #NodesPerPart=         333
   #NodesPerPart=         18
   #NodesPerPart=         27
   #NodesPerPart=         512
   #NodesPerPart=         17
   #NodesPerPart=         48
-  #NodesPerPart=         150
+  NodesPerPart=         300
   #NodesPerPart=         1331
-  PartitionType=       Sedacs
+  #PartitionType=       Sedacs
   #PartitionType=        METIS+SA
   #PartitionType=        METIS+KL
   #PartitionRefinement=  None
@@ -121,10 +121,10 @@ GSP2{
   #PartitionCount=      256
   #PartitionCount=      512
   #PartitionCount=      16
-  PartitionCount=      8
-  PartitionCountX=     2
-  PartitionCountY=     2
-  PartitionCountZ=     2
+  PartitionCount=      1
+  PartitionCountX=     1
+  PartitionCountY=     1
+  PartitionCountZ=     1
   #PartitionCount=     8
   #PartitionCount=      1024
   #PartitionCount=      32
@@ -161,6 +161,9 @@ XLBO{
 
 
 KERNEL{
+   XLBOLevel1= T
+   ScaledDelta= T
+   ScaledDeltaConstant= 0.2
    KernelType= ByParts
    #KernelType= Full
    #KernelType= ByBlocks
@@ -174,6 +177,7 @@ KERNEL{
 }
 
 GPMD{
+   AdaptiveTimeStep= T
    DoVelocityRescale= F
    #VRFactor= 1.0
    WriteTrajectory= T
diff --git a/examples/gpmdk/run/water/my_waterInput.in b/examples/gpmdk/run/water/my_waterInput.in
new file mode 100644
index 00000000..c89034bc
--- /dev/null
+++ b/examples/gpmdk/run/water/my_waterInput.in
@@ -0,0 +1,130 @@
+INPUT FILE FOR THE GPMD PROGRAM
+===============================
+
+#LATTE parameters
+Latte{
+  JobName= GPMD
+  #BMLType= Ellpack
+  BMLType= Dense
+  #Method= GSP2
+  #Method= SP2
+  #Method= Diag
+  #Method= DiagEf
+  Method= DiagEfFull
+  MDim= -1
+  #Threshold= 1.0d-5
+  Threshold= 0.0
+  Verbose= 2  #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all)
+  #Verbose= 10  #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all)
+  #SCF variables#
+  #StopAt= "gpmdcov_Energ"
+  #StopAt= "gpmdcov_DM_Min"
+  #StopAt= "gpmdcov_FirstCharges"
+  MPulay= 10
+  #ZMat= ZSP
+  ZMat= Diag
+  PulayCoeff= 0.1
+  #MixCoeff= 0.6  #VALID FOR WAT
+  MixCoeff= 0.2
+  SCFTol= 1.0d-8
+  MaxSCFIter= 500
+  CoulAcc= 1.0d-5
+  TimeRatio= 10.0
+  #TimeStep= 0.2
+  TimeStep= 0.2
+  #TimeStep= 0.00
+  MDSteps= 2000
+  ParamPath= "../../tests/latteTBparams"
+  CoordsFile= coords_300.dat
+  #CoordsFile= coords_2088.dat
+  NlistEach= 10
+  MuCalcType= FromParts
+  EFermi= -0.0
+  #kBT= 0.04308695
+  #kBT= 0.025
+  kBT= 0.2
+  Entropy= T
+  DoKernel= F
+}
+
+#SP2 Solver
+SP2{
+  MinSP2Iter=           10
+  MaxSP2Iter=           200
+  SP2Tol=   1.0d-5
+  SP2Conv= Rel
+}
+
+#Graph-based SP2 parameters
+GSP2{
+
+  BMLType= Ellpack
+  GraphElement=         Atom
+  #PartitionType=        Box
+  PartitionType=        Block
+  #PartitionType=       Sedacs
+  NLGraphCut= 4.5
+  CovGraphFact= 4.5
+  NodesPerPart= 300
+  PartitionCount=      1
+  PartitionCountX=     1
+  PartitionCountY=     1
+  PartitionCountZ=     1
+  GraphThreshold=      0.00001
+  ErrLimit=             1.0e-12
+  PartEach= 1000
+  SmallSubgraphs= T
+  Alpha= 10
+  Mdim= -1
+}
+
+
+#Sparse propagation of the inverse overlap
+ZSP{
+  Verbose= 1
+  NFirst= 8
+  NRefI= 3
+  NRefF= 1
+  Int= .true.
+  NumthreshI= 1.0d-8
+  NumthreshF= 1.0d-5
+}
+
+#Extended Lagrangian parameters
+XLBO{
+  JobName= XLBO
+  Verbose= 1
+  Mprg_init= 2
+  MaxSCFIter= 0
+  MaxSCFInitIter= 50
+  NumThresh= 0.0
+}
+
+
+KERNEL{
+   XLBOLevel1= T
+   ScaledDelta= T
+   ScaledDeltaConstant= 0.2
+   KernelType= ByParts
+   #KernelType= Full
+   #KernelType= ByBlocks
+   BuildAlways= F
+   RankNUpdate= 2
+   KernelMixing= T
+   InitialMixingWith= DIIS
+   UpdateEach= 1
+   UpdateAfterBuild= T
+   Verbose= 1
+}
+
+GPMD{
+   DoVelocityRescale= F
+   #VRFactor= 1.0
+   WriteTrajectory= F
+   WriteCoordsEach= 10
+   LangevinMethod= Siva
+   LangevinDynamics= F
+   LangevinGamma= 0.01
+   InitialTemperature= 300.0
+   SymmetrizeGraph= T
+}
diff --git a/examples/gpmdk/run/water/my_waterInput_mod.in b/examples/gpmdk/run/water/my_waterInput_mod.in
new file mode 100644
index 00000000..2ea8ad00
--- /dev/null
+++ b/examples/gpmdk/run/water/my_waterInput_mod.in
@@ -0,0 +1,130 @@
+INPUT FILE FOR THE GPMD PROGRAM
+===============================
+
+#LATTE parameters
+Latte{
+  JobName= GPMD
+  #BMLType= Ellpack
+  BMLType= Dense
+  #Method= GSP2
+  #Method= SP2
+  #Method= Diag
+  #Method= DiagEf
+  Method= DiagEfFull
+  MDim= -1
+  #Threshold= 1.0d-5
+  Threshold= 0.0
+  Verbose= 2  #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all)
+  #Verbose= 10  #Verbosity levels: Basic info(0), 1(Basic routines info), 2(Print Physics data), 3(Print Relevant Matrices), 5(Print auxiliary matrices), 10(Print all)
+  #SCF variables#
+  #StopAt= "gpmdcov_Energ"
+  #StopAt= "gpmdcov_DM_Min"
+  #StopAt= "gpmdcov_FirstCharges"
+  MPulay= 10
+  #ZMat= ZSP
+  ZMat= Diag
+  PulayCoeff= 0.1
+  #MixCoeff= 0.6  #VALID FOR WAT
+  MixCoeff= 0.2
+  SCFTol= 1.0d-8
+  MaxSCFIter= 500
+  CoulAcc= 1.0d-5
+  TimeRatio= 10.0
+  #TimeStep= 0.2
+  TimeStep= 0.4
+  #TimeStep= 0.00
+  MDSteps= 2000
+  ParamPath= "../../tests/latteTBparams"
+  CoordsFile= coords_300.dat
+  #CoordsFile= coords_2088.dat
+  NlistEach= 10
+  MuCalcType= FromParts
+  EFermi= -0.0
+  #kBT= 0.04308695
+  #kBT= 0.025
+  kBT= 0.2
+  Entropy= T
+  DoKernel= T
+}
+
+#SP2 Solver
+SP2{
+  MinSP2Iter=           10
+  MaxSP2Iter=           200
+  SP2Tol=   1.0d-5
+  SP2Conv= Rel
+}
+
+#Graph-based SP2 parameters
+GSP2{
+
+  BMLType= Ellpack
+  GraphElement=         Atom
+  #PartitionType=        Box
+  PartitionType=        Block
+  #PartitionType=       Sedacs
+  NLGraphCut= 4.5
+  CovGraphFact= 4.5
+  NodesPerPart= 300
+  PartitionCount=      1
+  PartitionCountX=     1
+  PartitionCountY=     1
+  PartitionCountZ=     1
+  GraphThreshold=      0.00001
+  ErrLimit=             1.0e-12
+  PartEach= 1000
+  SmallSubgraphs= T
+  Alpha= 10
+  Mdim= -1
+}
+
+
+#Sparse propagation of the inverse overlap
+ZSP{
+  Verbose= 1
+  NFirst= 8
+  NRefI= 3
+  NRefF= 1
+  Int= .true.
+  NumthreshI= 1.0d-8
+  NumthreshF= 1.0d-5
+}
+
+#Extended Lagrangian parameters
+XLBO{
+  JobName= XLBO
+  Verbose= 1
+  Mprg_init= 2
+  MaxSCFIter= 0
+  MaxSCFInitIter= 50
+  NumThresh= 0.0
+}
+
+
+KERNEL{
+   XLBOLevel1= T
+   ScaledDelta= T
+   ScaledDeltaConstant= 0.2
+   KernelType= ByParts
+   #KernelType= Full
+   #KernelType= ByBlocks
+   BuildAlways= F
+   RankNUpdate= 2
+   KernelMixing= T
+   InitialMixingWith= DIIS
+   UpdateEach= 1
+   UpdateAfterBuild= T
+   Verbose= 1
+}
+
+GPMD{
+   DoVelocityRescale= F
+   #VRFactor= 1.0
+   WriteTrajectory= F
+   WriteCoordsEach= 10
+   LangevinMethod= Siva
+   LangevinDynamics= T
+   LangevinGamma= 0.01
+   InitialTemperature= 300.0
+   SymmetrizeGraph= T
+}
diff --git a/examples/gpmdk/src/gpmdcov_init.F90 b/examples/gpmdk/src/gpmdcov_init.F90
index a39f56d6..0dc1d0a3 100644
--- a/examples/gpmdk/src/gpmdcov_init.F90
+++ b/examples/gpmdk/src/gpmdcov_init.F90
@@ -80,7 +80,7 @@ subroutine gpmdcov_Init(lib_on)
 
     !> Parsing specific variales for the gpmd code
     call gpmdcov_parse(trim(adjustl(inputfile)),gpmdt)
-    
+
     !> Parsing specific variales for controlling electronic structure output
     call gpmdcov_estructout_parse(trim(adjustl(inputfile)),estrout)
     
diff --git a/examples/gpmdk/src/gpmdcov_mdloop.F90 b/examples/gpmdk/src/gpmdcov_mdloop.F90
index 05142d3e..7c52c23b 100644
--- a/examples/gpmdk/src/gpmdcov_mdloop.F90
+++ b/examples/gpmdk/src/gpmdcov_mdloop.F90
@@ -37,7 +37,11 @@ subroutine gpmdcov_MDloop()
     real(dp) :: pressure_tensor(3,3)
     real(dp), allocatable :: saved_velocities(:,:)
     real(dp), allocatable :: saved_forces(:,:)
-    integer :: total_steps
+
+    real(dp) :: user_timestep,this_maxdisp,user_half_timestep
+    real(dp), parameter :: maxdist = 0.02
+    logical  :: first_substep_taken,half_timestep_flag
+    integer :: total_steps, print_mdstep
     integer :: cuda_error
     logical                           ::  newnl ! Indicates new neighbor list
     type(neighlist_type)              ::  nl2
@@ -75,7 +79,7 @@ end function cudaProfilerStop
     endif
     
     call gpmdcov_msI("gpmdcov_MDloop","In gpmdcov_MDloop ...",lt%verbose,myRank)
-    savets = lt%timestep
+    !savets = lt%timestep
     !do mdstep = -1,lt%mdsteps
     if(gpmdt%minimization_steps.ne.0)then
        saved_velocities = sy%velocity
@@ -87,12 +91,29 @@ end function cudaProfilerStop
     call gpmdcov_get_vol(sy%lattice_vector,sy%volr)
 
     total_steps = lt%mdsteps + gpmdt%minimization_steps
-    
-    if(gpmdt%freeze) then 
+
+    if(gpmdt%freeze) then
        call freeze(gpmdt%freezef,freeze_list,sy%velocity)
     endif
-    
-    do mdstep = 1,total_steps
+
+    ! user_timestep is a full timestep
+    ! user_half_timestep is a half timestep
+    ! first_substep_taken indicates that the first of 2 half timesteps was taken
+    ! half_timestep_flag indicates that 2 half timesteps were used
+    !   an output message is printed after the mdsteps line
+    ! print_mdstep is the mdstep used for output
+    !
+    user_timestep = lt%timestep
+    user_half_timestep = lt%timestep/2.0
+    first_substep_taken = .false.
+    half_timestep_flag = .false.
+    print_mdstep = 0
+    Time = 0.0
+
+    ! Loop continues until we've completed the requested number of user timesteps
+    mdstep = 0
+    do while (print_mdstep < total_steps)
+      mdstep = mdstep + 1
       !    if(mdstep < 0)then
       !            savets = lt%timestep
       !            lt%timestep = 0
@@ -125,6 +146,33 @@ end function cudaProfilerStop
         write(*,*)""
       endif
 
+      this_maxdisp = maxval(user_timestep*sy%velocity)
+      write(*,*)"Rank ", myRank, " for mdstep ", mdstep, "this_maxdisp = ", this_maxdisp
+
+      ! For dt/2 grid approach: force timestep splitting during initial history building
+      ! K=5: split first 4 print_mdsteps (gives 8 mdsteps at dt/2, >= 6 needed)
+      ! Then allow normal adaptive timestepping
+      if (gpmdt%adaptive_timestep .and. &
+          (print_mdstep <= 4 .or. &
+           (first_substep_taken .or.(this_maxdisp > maxdist)) .and. mdstep.gt.gpmdt%minimization_steps)) then
+        ! Only print when starting a new split (not when taking second half)
+        if (.not. first_substep_taken) then
+          write(*,*)"Rank ", myRank, " Splitting print_mdstep ", print_mdstep
+        endif
+        lt%timestep = user_half_timestep
+        half_timestep_flag = .true.
+
+        if (first_substep_taken) then
+          first_substep_taken = .false.
+        else
+          first_substep_taken = .true.
+        endif
+
+      else
+       lt%timestep = user_timestep
+       half_timestep_flag = .false.
+      endif
+
       maxv_atom_axis = MAXLOC(ABS(sy%velocity))
       call gpmdcov_msI("gpmdcov_MDloop","Maximum Velocity "//to_string(MAXVAL(ABS(sy%velocity)))//" &
         &for (atom,axis) = ("//to_string(maxv_atom_axis(2))//","//to_string(maxv_atom_axis(1))//")",lt%verbose,myRank)
@@ -147,8 +195,6 @@ end function cudaProfilerStop
       endif
       !! Total Energy in eV
       Energy = EKIN + EPOT;
-      !! Time in fs
-      Time = mdstep*lt%timestep;
 
       !! Statistical pressure
       do i = 1,3
@@ -162,6 +208,7 @@ end function cudaProfilerStop
       
       if(myRank == 1)then
         write(*,*)"Time [fs] = ",Time
+        write(*,*)"Time Step [fs] = ",lt%timestep 
         write(*,*)"Energy Kinetic [eV] = ",EKIN
         write(*,*)"Energy Potential [eV] = ",EPOT
         write(*,*)"Energy Total [eV] = ",Energy
@@ -187,6 +234,7 @@ end function cudaProfilerStop
                write(*,*)i,sy%velocity(1,i),sy%velocity(2,i),sy%velocity(3,i)
             enddo
          endif
+
          !> Update positions
          call gpmdcov_msMem("gpmdcov_mdloop", "Before updatecoords",lt%verbose,myRank)
          if(myRank == 1 .and. lt%verbose >= 1) call prg_timer_start(dyn_timer,"Update positions")
@@ -260,7 +308,17 @@ end function cudaProfilerStop
                 n = sy%net_charge
                 call gpmdcov_applyKernel(sy%net_charge,n,syprtk,KK0Res)
                 call prg_xlbo_nint_kernelTimesRes(sy%net_charge,n,n_0,&
-                     &n_1,n_2,n_3,n_4,n_5,mdstep,KK0Res,xl)
+                     &n_1,n_2,n_3,n_4,n_5,mdstep,KK0Res,xl,lt%timestep/user_timestep)
+
+                ! Synchronize XLBO charges across MPI ranks
+#ifdef DO_MPI
+                if (numRanks .gt. 1) then
+                  call prg_sumRealReduceN(n, sy%nats)
+                  call prg_sumRealReduceN(n_0, sy%nats)
+                  n = n / real(numRanks, dp)
+                  n_0 = n_0 / real(numRanks, dp)
+                endif
+#endif
               endif
               if(mdstep > 1 .and. kernel%rankNUpdate > 0 .and. &
                    & mod(mdstep,kernel%updateEach) == 0)then
@@ -268,7 +326,17 @@ end function cudaProfilerStop
 
                 !call gpmdcov_applyKernel(sy%net_charge,n,syprtk,KK0Res)
                 call prg_xlbo_nint_kernelTimesRes(sy%net_charge,n,n_0,&
-                     &n_1,n_2,n_3,n_4,n_5,mdstep,KK0Res,xl)
+                     &n_1,n_2,n_3,n_4,n_5,mdstep,KK0Res,xl,lt%timestep/user_timestep)
+
+                ! Synchronize XLBO charges across MPI ranks
+#ifdef DO_MPI
+                if (numRanks .gt. 1) then
+                  call prg_sumRealReduceN(n, sy%nats)
+                  call prg_sumRealReduceN(n_0, sy%nats)
+                  n = n / real(numRanks, dp)
+                  n_0 = n_0 / real(numRanks, dp)
+                endif
+#endif
                 !Use n > H >  to get q_min
                 ! call gpmdcov_DM_Min_Eig(1,sy%net_charge,.false.)
                 !Compute KK0Res
@@ -280,9 +348,9 @@ end function cudaProfilerStop
               deallocate(kernelTimesRes)
             else
               STOP "XLBOLevel1 not implemented for other than kernelType= ByParts"
-            endif
+            endif ! if by parts
 
-          else
+          else ! if XLBO level 1
             if(kernel%kernelType == "ByParts")then
               allocate(kernelTimesRes(sy%nats))
               if(mdstep.le.1)then
@@ -324,20 +392,53 @@ end function cudaProfilerStop
               endif
               call gpmdcov_msMem("gpmdcov_mdloop", "Before prg_xlbo_nint_kernelTimesRes",lt%verbose,myRank)
               call prg_xlbo_nint_kernelTimesRes(sy%net_charge,n,n_0,&
-                   &n_1,n_2,n_3,n_4,n_5,mdstep,KK0Res,xl)
+                   &n_1,n_2,n_3,n_4,n_5,mdstep,KK0Res,xl,lt%timestep/user_timestep)
               call gpmdcov_msMem("gpmdcov_mdloop", "After prg_xlbo_nint_kernelTimesRes",lt%verbose,myRank)
+
+              ! Synchronize XLBO charges across MPI ranks
+#ifdef DO_MPI
+              if (numRanks .gt. 1) then
+                call prg_sumRealReduceN(n, sy%nats)
+                call prg_sumRealReduceN(n_0, sy%nats)
+                n = n / real(numRanks, dp)
+                n_0 = n_0 / real(numRanks, dp)
+              endif
+#endif
               deallocate(kernelTimesRes)
-            else
+            else ! if byparts
               call gpmdcov_msMem("gpmdcov_mdloop", "Before prg_xlbo_nint_kernel",lt%verbose,myRank)
-              call prg_xlbo_nint_kernel(sy%net_charge,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,Ker,xl)
+              call prg_xlbo_nint_kernel(sy%net_charge,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,Ker,xl,lt%timestep/user_timestep)
               call gpmdcov_msMem("gpmdcov_mdloop", "After prg_xlbo_nint_kernel",lt%verbose,myRank)
-            endif
-          endif
-        else
+
+              ! Synchronize XLBO charges across MPI ranks
+#ifdef DO_MPI
+              if (numRanks .gt. 1) then
+                if (myRank == 1) write(*,*) "DEBUG: Rank 1 syncing XLBO charges at mdstep ", mdstep
+                call prg_sumRealReduceN(n, sy%nats)
+                call prg_sumRealReduceN(n_0, sy%nats)
+                n = n / real(numRanks, dp)
+                n_0 = n_0 / real(numRanks, dp)
+              endif
+#endif
+            endif ! byparts
+          endif ! if XLBO level 1
+        else ! if kernel
           call gpmdcov_msMem("gpmdcov_mdloop", "Before prg_xlbo_nint",lt%verbose,myRank)
-          
+
           if(gpmdt%xlboon)then
-                call prg_xlbo_nint(sy%net_charge,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,xl)
+                call prg_xlbo_nint(sy%net_charge,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,xl,lt%timestep/user_timestep)
+
+                ! Synchronize XLBO charges across MPI ranks to prevent divergence
+                ! Both n and n_0 need sync since n_0=n is done inside prg_xlbo_nint
+#ifdef DO_MPI
+                if (numRanks .gt. 1) then
+                  if (myRank == 1) write(*,*) "DEBUG: Rank 1 syncing XLBO charges at mdstep ", mdstep
+                  call prg_sumRealReduceN(n, sy%nats)
+                  call prg_sumRealReduceN(n_0, sy%nats)
+                  n = n / real(numRanks, dp)
+                  n_0 = n_0 / real(numRanks, dp)
+                endif
+#endif
           else
                 n = sy%net_charge
           endif
@@ -359,7 +460,7 @@ end function cudaProfilerStop
 
       !> Update neighbor list (Actialized every nlisteach times steps)
       mls_md1 = mls()
-      if(mod(mdstep,lt%nlisteach) == 0 .or. mdstep == 0 .or. mdstep == 1)then
+      if((mod(mdstep,lt%nlisteach) == 0 .or. mdstep == 0 .or. mdstep == 1))then
            call gpmdcov_msMemGPU("mdloop","Before NeighborList",lt%verbose,myRank)
         call gpmdcov_msMem("gpmdcov_mdloop", "Before build_nlist_int",lt%verbose,myRank)
         !call gpmdcov_destroy_nlist(nl,lt%verbose)
@@ -369,13 +470,8 @@ end function cudaProfilerStop
 #ifdef USE_NVTX
            call gpmdStartRange("build_nlist_sparse_sedacs",3)
 #endif
-           !call gpmdcov_destroy_nlist(nl2,lt%verbose)
-#ifdef USE_OFFLOAD
            call gpmdcov_build_nlist_sedacs(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks)
-#else
-           call gpmdcov_build_nlist_sedacs(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks)
-           !call gpmdcov_build_nlist_sparse_v2(sy%coordinate,sy%lattice_vector,coulcut,nl,lt%verbose,myRank,numRanks)
-#endif
+
            ! if(any(nl2%nrnnstruct.ne.nl%nrnnstruct))then
            !    write(*,*)"DEBUG: nrnnstruct not equal"
            !    do k = 1,size(nl%nrnnstruct)
@@ -433,12 +529,14 @@ end function cudaProfilerStop
 #endif
 
            call gpmdcov_Part(2)
+
 #ifdef USE_NVTX
            call gpmdEndRange
 #endif
       call gpmdcov_msMem("gpmdcov_mdloop", "After gpmdcov_Part",lt%verbose,myRank)
       call gpmdcov_msI("gpmdcov_MDloop","Time for gpmdcov_Part &
            &"//to_string(mls() - mls_i)//" ms",lt%verbose,myRank)
+
       !> Reprg_initialize parts.
       mls_i = mls()
       call gpmdcov_msMem("gpmdcov_mdloop", "Before gpmdcov_InitParts",lt%verbose,myRank)
@@ -460,7 +558,8 @@ end function cudaProfilerStop
       mls_md1 = mls()
       resnorm = 0.0_dp
 
-      if((mdstep >= 2) .and. (.not. (kernel%xlbolevel1.and.lt%doKernel))) resnorm =  norm2(sy%net_charge - n)/sqrt(dble(sy%nats))
+      !if((mdstep >= 2) .and. (.not. (kernel%xlbolevel1.and.lt%doKernel))) resnorm =  norm2(sy%net_charge - n)/sqrt(dble(sy%nats))
+      if((mdstep >= 2) .and. (.not. kernel%xlbolevel1)) resnorm =  norm2(sy%net_charge - n)/sqrt(dble(sy%nats))
 
       Nr_SCF_It = xl%maxscfiter;
       !> Use SCF the first MD steps
@@ -513,7 +612,7 @@ end function cudaProfilerStop
 #ifdef USE_NVTX
            call gpmdEndRange
 #endif
-      if(kernel%xlbolevel1.and.lt%doKernel)then
+       if(kernel%xlbolevel1.and.lt%doKernel)then
         allocate(n1(sy%nats))
         if(mdstep > 1)then
           !sy%net_charge = n
@@ -554,6 +653,14 @@ end function cudaProfilerStop
 
       mls_md1 = mls()
       call gpmdcov_msI("gpmdcov_MDloop","ResNorm = "//to_string(resnorm),lt%verbose,myRank)
+
+      ! Update print_mdstep counter on all ranks (used for forced splitting decision)
+      if(mdstep.gt.gpmdt%minimization_steps)then
+         if (.not.first_substep_taken)then
+           print_mdstep = print_mdstep + 1
+         endif
+      endif
+
       if(myRank == 1)then
          if(mdstep.le.gpmdt%minimization_steps)then
             if(.not.gpmdt%anneal_graph)then
@@ -564,10 +671,15 @@ end function cudaProfilerStop
                     &mdstep," ", Energy," ", egap_glob," ", resnorm," ", Temp
             endif
          else
-            write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", &
-                 &mdstep-gpmdt%minimization_steps," ", Energy," ", egap_glob," ", resnorm," ", Temp
+            ! Write output (rank 1 only)
+            if (.not.first_substep_taken)then
+              write(*,'(A35,I15,A1,F18.5,A1,ES12.5,A1,ES12.5,A1,ES12.5)')"Mdstep, Energy, Egap, Resnorm, Temp", &
+                 &print_mdstep," ",  Energy," ", egap_glob," ", resnorm," ", Temp
+              if (half_timestep_flag)then
+                write(*,*) "WARNING: Two half timesteps were performed for step ", print_mdstep
+              endif
+            endif
          endif
-        !write(*,*)"Step, Energy, EGap, Resnorm", mdstep, Energy, egap_glob, resnorm
       endif
 #ifdef USE_NVTX
       call gpmdStartRange("EnergAndForces",7)
@@ -699,16 +811,16 @@ end function cudaProfilerStop
 #ifdef USE_NVTX
       call gpmdStartRange("Write trajectory",3)
 #endif
-      if(gpmdt%writetraj .and. myRank == 1 .and. mdstep.ge.gpmdt%minimization_steps)then
+      if(gpmdt%writetraj .and. myRank == 1 .and. mdstep.ge.gpmdt%minimization_steps .and. first_substep_taken .eqv. .false.)then
         if((gpmdt%traj_format .eq. "XYZ").and. &
-          (mod(mdstep-gpmdt%minimization_steps,gpmdt%writetreach).eq.0.or. &
-           (mdstep-gpmdt%minimization_steps).eq.1))then
-           call prg_write_trajectory(sy,mdstep-gpmdt%minimization_steps,gpmdt%writetreach,&
-                &lt%timestep,adjustl(trim(lt%jobname))//"_trajectory","xyz")
+          (mod(print_mdstep,gpmdt%writetreach).eq.0.or. &
+           (print_mdstep).eq.1))then
+           call prg_write_trajectory(sy,print_mdstep,gpmdt%writetreach,&
+                &user_timestep,adjustl(trim(lt%jobname))//"_trajectory","xyz")
            call prg_write_system(sy,adjustl(trim(lt%jobname))//"_latest","pdb")
         else
-           call prg_write_trajectory(sy,mdstep-gpmdt%minimization_steps,gpmdt%writetreach,&
-             &lt%timestep,adjustl(trim(lt%jobname))//"_trajectory","pdb")
+           call prg_write_trajectory(sy,print_mdstep,gpmdt%writetreach,&
+             &user_timestep,adjustl(trim(lt%jobname))//"_trajectory","pdb")
         endif
      endif
 #ifdef USE_NVTX
@@ -725,7 +837,7 @@ end function cudaProfilerStop
       
       ! Save MD state each 120 steps
       if(gpmdt%dumpeach .gt. 0)then
-         if(mod(mdstep-gpmdt%minimization_steps,gpmdt%dumpeach) == 0)call gpmdcov_dump()
+         if(mod(print_mdstep,gpmdt%dumpeach) == 0)call gpmdcov_dump()
       endif
       
       if(mdstep.eq.gpmdt%minimization_steps)then
@@ -737,6 +849,10 @@ end function cudaProfilerStop
          endif
       endif
             
+      !! Time in fs
+      !Time = mdstep*lt%timestep;
+      Time = Time + lt%timestep
+      
     enddo
     ! End of MD loop.
 
diff --git a/examples/gpmdk/src/gpmdcov_parser.F90 b/examples/gpmdk/src/gpmdcov_parser.F90
index fb35cc37..d629d917 100644
--- a/examples/gpmdk/src/gpmdcov_parser.F90
+++ b/examples/gpmdk/src/gpmdcov_parser.F90
@@ -180,7 +180,10 @@ module gpmdcov_parser_mod
 
     !> Rescale velocities from restart file to match initial temperature
     logical :: rescale_restart_vel
-    
+
+    !> Use adaptive timestep splitting when max displacement exceeds threshold
+    logical :: adaptive_timestep
+
   end type gpmd_type
 
   !> electrontic structure output type
@@ -243,7 +246,7 @@ subroutine gpmdcov_parse(filename,gpmdt)
     implicit none 
     character(len=*), intent(in) :: filename
     type(gpmd_type), intent(inout) :: gpmdt
-    integer, parameter :: nkey_char = 6, nkey_int = 15, nkey_re = 7, nkey_log = 22
+    integer, parameter :: nkey_char = 6, nkey_int = 15, nkey_re = 7, nkey_log = 23
     integer :: i
     real(dp) :: realtmp
     character(20) :: dummyc
@@ -275,11 +278,11 @@ subroutine gpmdcov_parse(filename,gpmdt)
          &'ComputeCurrents=', 'TranslateAndFoldToBox=', 'UseVectSKBlock=', 'ApplyVoltage=','XLBO=',&
          'CoarseQMD=',&
          &'UseDispersion=','UseFreeze=','SymmetrizeGraph=','AnnealGraph=',&
-         &'UseCustomSeed=','UseRandomSeed=','RescaleRestartVelocities=']
+         &'UseCustomSeed=','UseRandomSeed=','RescaleRestartVelocities=','AdaptiveTimeStep=']
     logical :: valvector_log(nkey_log) = (/&
          &.false.,.false.,.false.,.false.,.false.,.false.,.false.,.false.,.false., &
          &.false.,.True.,.false.,.false.,.true.,.false.,.false.,.false.,.false.,.false.,&
-         &.false.,.false.,.false./)
+         &.false.,.false.,.false.,.false./)
 
     !Start and stop characters
     character(len=50), parameter :: startstop(2) = [character(len=50) :: &
@@ -412,7 +415,8 @@ subroutine gpmdcov_parse(filename,gpmdt)
     gpmdt%usecustomseed = valvector_log(20)
     gpmdt%userandomseed = valvector_log(21)
     gpmdt%rescale_restart_vel = valvector_log(22)
-    
+    gpmdt%adaptive_timestep = valvector_log(23)
+
     if(gpmdt%applyv)then 
         gpmdt%voltagef = valvector_char(5)
     endif
diff --git a/scripts/build_mac.sh b/scripts/build_mac.sh
new file mode 100644
index 00000000..26eec553
--- /dev/null
+++ b/scripts/build_mac.sh
@@ -0,0 +1 @@
+BML_DIR=/Users/mewall/packages/gpmd/bml/install PROGRESS_EXAMPLES=yes bash build.sh install
diff --git a/scripts/compute_alpha_table.py b/scripts/compute_alpha_table.py
new file mode 100644
index 00000000..c244cf3c
--- /dev/null
+++ b/scripts/compute_alpha_table.py
@@ -0,0 +1,55 @@
+#!/usr/bin/env python
+"""Compute alpha values using formula: alpha = 0.018 * 3.0 / d_K, capped at alpha_max/2"""
+
+import numpy as np
+
+# d_K values from XLBO_K5_dK table
+d_K = [
+    0.75000000, 0.28947368, 0.38888889, 0.41666667,
+    0.32926829, 0.42765957, 0.35294118, 0.44444444,
+    6.11111111, 4.39090909, 2.25000000, 3.46153846,
+    11.90476190, 6.13636364, 4.70000000, 2.82692308,
+    7.82608696, 3.00000000, 4.89130435, 2.53846154,
+    8.24137931, 3.72972973, 5.77777778, 3.08571429,
+    11.74193548, 4.56756757, 7.00000000, 3.32954545,
+    10.66666667, 4.34782609, 6.25000000, 3.00000000
+]
+
+# Alpha_max values from stability analysis
+alpha_max = [
+    0.144000, 0.372414, 0.276923, 0.258992,
+    0.327273, 0.252323, 0.305512, 0.242969,
+    0.017660, 0.024610, 0.048000, 0.031206,
+    0.009067, 0.011402, 0.022979, 0.038162,
+    0.015428, 0.036000, 0.022087, 0.042521,
+    0.013091, 0.028954, 0.018685, 0.034951,
+    0.009192, 0.023633, 0.015428, 0.032432,
+    0.010122, 0.024941, 0.017280, 0.036000
+]
+
+base_alpha = 0.018
+d_K_base = 3.0  # Pattern 31
+
+print("Pattern | d_K      | formula_alpha | alpha_max | alpha_max/2 | Use (min)  | Capped?")
+print("--------|----------|---------------|-----------|-------------|------------|--------")
+
+alpha_values = []
+for i in range(32):
+    formula = base_alpha * d_K_base / d_K[i]
+    half_max = alpha_max[i] / 2.0
+    use_alpha = min(formula, half_max)
+    capped = "YES" if formula > half_max else "NO"
+    alpha_values.append(use_alpha)
+    print(f"{i:7d} | {d_K[i]:8.5f} | {formula:13.6f} | {alpha_max[i]:9.6f} | {half_max:11.6f} | {use_alpha:10.6f} | {capped:7s}")
+
+print("\n" + "="*80)
+print("Fortran array for prg_xlbo_mod.F90:")
+print("="*80)
+print()
+print("  real(dp), parameter :: XLBO_K5_alpha(0:31) = [ &")
+for i in range(0, 32, 4):
+    values = ", ".join([f"{alpha_values[j]:14.6f}_dp" for j in range(i, min(i+4, 32))])
+    if i < 28:
+        print(f"    {values}, &")
+    else:
+        print(f"    {values}]")
diff --git a/src/prg_xlbo_mod.F90 b/src/prg_xlbo_mod.F90
index 98c8416f..ddd2086e 100644
--- a/src/prg_xlbo_mod.F90
+++ b/src/prg_xlbo_mod.F90
@@ -15,17 +15,117 @@ module prg_xlbo_mod
 
   integer, parameter :: dp = kind(1.0d0)
 
-  !> Coefficients for modified Verlet integration
+  !> Coefficients for K=5 (6-point history) XLBO dissipation - uniform timestep
   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;
+  real(dp), parameter :: C5 = -1.0_dp
+  real(dp), parameter :: kappa = 1.82_dp
+  real(dp), parameter :: alpha = 0.018_dp
+
+  !> K=5 Variable Timestep Coefficient Lookup Tables
+  !! Fixed normalization: c_4 = 4.0, c_5 = -1.0 for all patterns
+  !! Indexed by 5-bit pattern: bit 0 = most recent dt, bit 4 = oldest dt
+  !! Bit value: 0 = half step (dt/2), 1 = full step (dt)
+  real(dp), parameter :: XLBO_K5_C0(0:31) = [ &
+    -6.00000000_dp, -0.85714286_dp, -1.92857143_dp, -0.32539683_dp, &
+    -1.33333333_dp,  0.08571429_dp,  0.87500000_dp,  0.40625000_dp, &
+    22.00000000_dp,  4.92857143_dp, 16.00000000_dp,  4.23809524_dp, &
+    44.33333333_dp,  9.57142857_dp, 28.37500000_dp,  7.47656250_dp, &
+   -62.00000000_dp,-10.50000000_dp,-29.42857143_dp, -6.58730159_dp, &
+   -63.00000000_dp,-11.34285714_dp,-30.75000000_dp, -7.11718750_dp, &
+   -98.00000000_dp,-14.71428571_dp,-39.00000000_dp, -7.83333333_dp, &
+   -75.66666667_dp,-11.80000000_dp,-30.00000000_dp, -6.00000000_dp]
+
+  real(dp), parameter :: XLBO_K5_C1(0:31) = [ &
+    14.00000000_dp,  3.00000000_dp,  3.00000000_dp,  0.52380952_dp, &
+     1.33333333_dp, -1.94285714_dp, -2.12500000_dp, -1.64062500_dp, &
+   -64.00000000_dp,-31.00000000_dp,-31.00000000_dp,-14.28571429_dp, &
+  -110.33333333_dp,-46.28571429_dp,-50.62500000_dp,-21.72265625_dp, &
+   160.00000000_dp, 53.00000000_dp, 53.00000000_dp, 19.23809524_dp, &
+   147.00000000_dp, 47.77142857_dp, 52.25000000_dp, 18.63671875_dp, &
+   245.00000000_dp, 68.00000000_dp, 68.00000000_dp, 21.00000000_dp, &
+   170.66666667_dp, 44.80000000_dp, 49.00000000_dp, 14.00000000_dp]
+
+  real(dp), parameter :: XLBO_K5_C2(0:31) = [ &
+    -8.00000000_dp, -0.57142857_dp,  0.71428571_dp,  2.05555556_dp, &
+     1.66666667_dp,  3.70000000_dp,  3.06250000_dp,  3.06250000_dp, &
+    67.00000000_dp, 47.28571429_dp, 34.00000000_dp, 26.66666667_dp, &
+    79.33333333_dp, 48.00000000_dp, 32.81250000_dp, 23.51562500_dp, &
+  -133.00000000_dp,-63.00000000_dp,-40.28571429_dp,-23.77777778_dp, &
+   -94.00000000_dp,-42.80000000_dp,-27.12500000_dp,-15.42187500_dp, &
+  -192.00000000_dp,-73.14285714_dp,-44.00000000_dp,-19.83333333_dp, &
+  -102.66666667_dp,-35.70000000_dp,-21.00000000_dp, -8.00000000_dp]
+
+  real(dp), parameter :: XLBO_K5_C3(0:31) = [ &
+    -3.00000000_dp, -4.57142857_dp, -4.78571429_dp, -5.25396825_dp, &
+    -4.66666667_dp, -4.84285714_dp, -4.81250000_dp, -4.82812500_dp, &
+   -28.00000000_dp,-24.21428571_dp,-22.00000000_dp,-19.61904762_dp, &
+   -16.33333333_dp,-14.28571429_dp,-13.56250000_dp,-12.26953125_dp, &
+    32.00000000_dp, 17.50000000_dp, 13.71428571_dp,  8.12698413_dp, &
+     7.00000000_dp,  3.37142857_dp,  2.62500000_dp,  0.90234375_dp, &
+    42.00000000_dp, 16.85714286_dp, 12.00000000_dp,  3.66666667_dp, &
+     4.66666667_dp, -0.30000000_dp, -1.00000000_dp, -3.00000000_dp]
+
+  real(dp), parameter :: XLBO_K5_C4(0:31) = [ &
+     4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp, &
+     4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp, &
+     4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp, &
+     4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp,  4.0_dp]
+
+  real(dp), parameter :: XLBO_K5_C5(0:31) = [ &
+    -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, &
+    -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, &
+    -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, &
+    -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp, -1.0_dp]
+
+  real(dp), parameter :: XLBO_K5_dK(0:31) = [ &
+     0.75000000_dp,  0.28571429_dp,  0.39285714_dp,  0.17063492_dp, &
+     0.33333333_dp,  0.06785714_dp,  0.01562500_dp,  0.07812500_dp, &
+     2.00000000_dp,  1.14285714_dp,  2.25000000_dp,  1.38095238_dp, &
+     5.08333333_dp,  2.71428571_dp,  4.70312500_dp,  2.83203125_dp, &
+     7.00000000_dp,  3.00000000_dp,  4.89285714_dp,  2.53968254_dp, &
+     8.25000000_dp,  3.72857143_dp,  5.78125000_dp,  3.08984375_dp, &
+    11.75000000_dp,  4.57142857_dp,  7.00000000_dp,  3.33333333_dp, &
+    10.66666667_dp,  4.32500000_dp,  6.25000000_dp,  3.00000000_dp]
+
+  !> Pattern-specific alpha values for K=5 XLBO dissipation
+  !! Conservative scaling: alpha = 0.000796 × 3.0 / d_K
+  !! All values ≤ α_max/2, ensuring stability for all 32 patterns
+  !! Limiting pattern: 6 (d_K=0.015625, exactly at α_max/2)
+  !! Average safety margin: 20× above usage (suitable for large systems)
+  !! Indexed by 5-bit pattern: bit 0 = most recent dt, bit 4 = oldest dt
+  !> Pattern-specific alpha values for K=5 variable timesteps
+  !! Computed as: alpha = 0.018 * 3.0 / d_K, capped at alpha_max/2
+  !! This ensures stability while maximizing dissipation for each pattern
+  real(dp), parameter :: XLBO_K5_alpha(0:31) = [ &
+        0.072000_dp,     0.186207_dp,     0.138461_dp,     0.129496_dp, &
+        0.163636_dp,     0.126162_dp,     0.152756_dp,     0.121484_dp, &
+        0.008830_dp,     0.012298_dp,     0.024000_dp,     0.015600_dp, &
+        0.004534_dp,     0.005701_dp,     0.011489_dp,     0.019081_dp, &
+        0.006900_dp,     0.018000_dp,     0.011040_dp,     0.021261_dp, &
+        0.006546_dp,     0.014477_dp,     0.009343_dp,     0.017476_dp, &
+        0.004596_dp,     0.011817_dp,     0.007714_dp,     0.016216_dp, &
+        0.005061_dp,     0.012420_dp,     0.008640_dp,     0.018000_dp]
+
+  !> Coefficients for K=10 (11-point history) XLBO dissipation
+  !> From Niklasson et al. JCP 2009 Table I extended
+  real(dp), parameter :: C0_K10 = -858.0_dp
+  real(dp), parameter :: C1_K10 = 2652.0_dp
+  real(dp), parameter :: C2_K10 = -3094.0_dp
+  real(dp), parameter :: C3_K10 = 1496.0_dp
+  real(dp), parameter :: C4_K10 = 272.0_dp
+  real(dp), parameter :: C5_K10 = -952.0_dp
+  real(dp), parameter :: C6_K10 = 731.0_dp
+  real(dp), parameter :: C7_K10 = -322.0_dp
+  real(dp), parameter :: C8_K10 = 88.0_dp
+  real(dp), parameter :: C9_K10 = -14.0_dp
+  real(dp), parameter :: C10_K10 = 1.0_dp
+  real(dp), parameter :: kappa_K10 = 1.88_dp
+  real(dp), parameter :: alpha_K10 = 0.036e-3_dp
 
-  !> Coefficients for modified Verlet integration
-  real(dp), parameter :: kappa = 1.82_dp;
-  real(dp), parameter :: alpha = 0.018_dp;
   real(dp), parameter :: cc = 0.9_dp;        ! Scaled prg_delta kernel
 
   !> General xlbo solver type
@@ -50,6 +150,11 @@ module prg_xlbo_mod
     !> Scaled prg_delta Kernel
     real(dp) :: cc
 
+    !> Timestep history for adaptive time step
+    !> Size 10 to support K=10 (11-point history); only first 5 used for K=5
+    real(dp) :: dt_history(10)
+    integer :: nsteps_taken
+
   end type xlbo_type
 
   public :: prg_parse_xlbo, prg_xlbo_nint, prg_xlbo_nint_kernel, prg_xlbo_fcoulupdate
@@ -102,27 +207,54 @@ subroutine prg_parse_xlbo(xlbo,filename)
     xlbo%threshold = valvector_re(1)
     xlbo%cc = valvector_re(2)
 
-    !Logicals
-
     !Integers
     xlbo%verbose = valvector_int(1)
     xlbo%minit = valvector_int(2)
     xlbo%maxscfiter = valvector_int(3)
     xlbo%maxscfinititer = valvector_int(4)
 
+    !Initialize timestep history
+    xlbo%dt_history = 0.0_dp
+    xlbo%nsteps_taken = 0
+
   end subroutine prg_parse_xlbo
 
 
+
+  !> Compute K=5 variable timestep lookup index from dt_history
+  !! \brief Converts dt_history into a 5-bit integer for coefficient lookup
+  !! \param dt_history Timestep history (most recent first, 5 elements)
+  !! \return Bit pattern: 0 = half step (dt/2), 1 = full step (dt)
+  function get_K5_history_index(dt_history) result(index)
+    implicit none
+    real(dp), intent(in) :: dt_history(5)
+    integer :: index
+    integer :: k
+
+    index = 0
+    do k = 1, 5
+      ! If timestep is close to 1.0 (full step), set bit k-1
+      if (abs(dt_history(k) - 1.0_dp) < 0.1_dp) then
+        index = ibset(index, k-1)
+      endif
+    end do
+  end function get_K5_history_index
+
   !> This routine integrates the dynamical variable "n"
   !! \param charges
-  subroutine prg_xlbo_nint(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,xl)
+  subroutine prg_xlbo_nint(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,xl,dt)
     implicit none
     real(dp), allocatable, intent(inout) :: n(:), n_0(:), n_1(:), n_2(:), n_3(:), n_4(:), n_5(:)
     real(dp), allocatable, intent(in) :: charges(:)
-    type(xlbo_type), intent(in) :: xl
-
+    type(xlbo_type), intent(inout) :: xl
     integer, intent(in) :: mdstep
+    real(dp), intent(in), optional :: dt
     integer :: nats
+    real(dp) :: kappa_use, alpha_use
+    logical :: allow_adaptive_timestep, use_K10
+    integer :: hist_idx
+    real(dp) :: C0_use, C1_use, C2_use, C3_use, C4_use, C5_use, d_K_use
+    real(dp) :: dt_n, dt_prev, r, P_n_coeff, P_n1_coeff, kappa_alpha_scale
 
     nats = size(charges,dim=1)
 
@@ -134,8 +266,8 @@ subroutine prg_xlbo_nint(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,xl)
       allocate(n_3(nats))
       allocate(n_4(nats))
       allocate(n_5(nats))
-    endif
-
+   endif
+   
     if(mdstep.le.1)then
       n = charges;
       n_0 = charges;
@@ -144,26 +276,114 @@ subroutine prg_xlbo_nint(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,xl)
       n_3 = charges;
       n_4 = charges;
       n_5 = charges;
+      xl%dt_history = 0.0_dp
+      xl%nsteps_taken = 0
     endif
 
-    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);
-    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n;
+    allow_adaptive_timestep = present(dt) .and. xl%nsteps_taken >= 6
+
+    kappa_use = kappa
+    alpha_use = alpha
+
+    ! Use pattern-specific alpha for early steps (during warmup before full history)
+    if (present(dt) .and. .not. allow_adaptive_timestep) then
+      ! Look up pattern-specific alpha based on current dt_history
+      hist_idx = get_K5_history_index(xl%dt_history(1:5))
+      alpha_use = XLBO_K5_alpha(hist_idx)
+    endif
+
+    ! Compute variable timestep Verlet coefficients
+    if (present(dt)) then
+      ! Use current timestep (dt) and previous timestep (xl%dt_history(1))
+      dt_n = dt                    ! Current timestep (input parameter)
+      dt_prev = xl%dt_history(1)   ! Previous timestep
+
+      ! Check if we have valid history (dt_prev > 0)
+      if (dt_prev > 1.0e-12_dp) then
+        ! Compute Verlet coefficients for variable timesteps
+        r = dt_n / dt_prev  ! Timestep ratio
+        P_n_coeff = 1.0_dp + r
+        P_n1_coeff = r
+
+        ! Compute kappa and alpha scaling factor
+        ! This comes from: 0.5 * dt_n * (dt_n + dt_{n-1})
+        kappa_alpha_scale = 0.5_dp * dt_n * (dt_n + dt_prev)
+      else
+        ! First step or dt_prev not set yet - use uniform coefficients
+        P_n_coeff = 2.0_dp
+        P_n1_coeff = 1.0_dp
+        kappa_alpha_scale = dt_n * dt_n  ! dt^2 for first step
+      endif
+    else
+      ! Uniform timestep: standard Verlet coefficients
+      P_n_coeff = 2.0_dp
+      P_n1_coeff = 1.0_dp
+      kappa_alpha_scale = 1.0_dp
+    endif
+
+    if (allow_adaptive_timestep) then
+          ! New method: Use variable timestep coefficients directly (no interpolation)
+
+          ! Get coefficient index from timestep history
+          hist_idx = get_K5_history_index(xl%dt_history(1:5))
+
+          ! Lookup coefficients for this specific history pattern
+          C0_use = XLBO_K5_C0(hist_idx)
+          C1_use = XLBO_K5_C1(hist_idx)
+          C2_use = XLBO_K5_C2(hist_idx)
+          C3_use = XLBO_K5_C3(hist_idx)
+          C4_use = XLBO_K5_C4(hist_idx)
+          C5_use = XLBO_K5_C5(hist_idx)
+          d_K_use = XLBO_K5_dK(hist_idx)
+
+          ! Use pattern-specific alpha value (capped at alpha_max/2 for stability)
+          alpha_use = XLBO_K5_alpha(hist_idx)
+
+          ! Integration using raw charges with variable coefficients and variable timestep Verlet
+          n = P_n_coeff*n_0 - P_n1_coeff*n_1 + xl%cc*kappa_alpha_scale*kappa_use*(charges-n) &
+               + alpha_use*(C0_use*n_0+C1_use*n_1+C2_use*n_2+C3_use*n_3+C4_use*n_4+C5_use*n_5)
+
+    else
+      ! Integration using raw charges with variable timestep Verlet
+        n = P_n_coeff*n_0 - P_n1_coeff*n_1 + xl%cc*kappa_alpha_scale*kappa_use*(charges-n) &
+             + alpha_use*(C0*n_0+C1*n_1+C2*n_2+C3*n_3+C4*n_4+C5*n_5)
+    endif
+
+    ! Shift history arrays
+    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n
+
+    ! Update timestep history if dt provided (dt is ratio: 1.0 for full step, 0.5 for half step)
+    ! History stores ratios that indicate spacing relative to user timestep
+    ! Interpolation always maps to fixed dt/2 uniform grid for stability
+    if (present(dt)) then
+      xl%dt_history(5) = xl%dt_history(4)
+      xl%dt_history(4) = xl%dt_history(3)
+      xl%dt_history(3) = xl%dt_history(2)
+      xl%dt_history(2) = xl%dt_history(1)
+      xl%dt_history(1) = dt  ! Store ratio (1.0 = full user timestep, 0.5 = half user timestep)
+      xl%nsteps_taken = xl%nsteps_taken + 1
+    endif
 
   end subroutine prg_xlbo_nint
 
 
   !> This routine integrates the dynamical variable "n"
   !! \param charges
-  subroutine prg_xlbo_nint_kernel(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,kernel,xl)
+  subroutine prg_xlbo_nint_kernel(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,kernel,xl,dt)
     implicit none
     real(dp), allocatable, intent(inout) :: n(:), n_0(:), n_1(:), n_2(:), n_3(:), n_4(:), n_5(:)
     real(dp), allocatable, intent(in) :: charges(:)
     real(dp), allocatable, intent(in) :: kernel(:,:)
-    type(xlbo_type), intent(in) :: xl
-
+    type(xlbo_type), intent(inout) :: xl
     integer, intent(in) :: mdstep
+    real(dp), intent(in), optional :: dt
     integer :: nats
+    real(dp) :: kappa_use, alpha_use
+    real(dp), allocatable :: KK0n(:)
+    logical :: allow_adaptive_timestep
+    integer :: hist_idx
+    real(dp) :: C0_use, C1_use, C2_use, C3_use, C4_use, C5_use, d_K_use
+    real(dp) :: dt_n, dt_prev, r, P_n_coeff, P_n1_coeff, kappa_alpha_scale
 
     nats = size(charges,dim=1)
 
@@ -185,6 +405,49 @@ subroutine prg_xlbo_nint_kernel(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,kernel,
       n_3 = charges;
       n_4 = charges;
       n_5 = charges;
+      xl%dt_history = 0.0_dp
+      xl%nsteps_taken = 0
+    endif
+
+    allow_adaptive_timestep = present(dt) .and. xl%nsteps_taken >= 6
+
+    kappa_use = kappa
+    alpha_use = alpha
+
+    ! Use pattern-specific alpha for early steps (during warmup before full history)
+    if (present(dt) .and. .not. allow_adaptive_timestep) then
+      ! Look up pattern-specific alpha based on current dt_history
+      hist_idx = get_K5_history_index(xl%dt_history(1:5))
+      alpha_use = XLBO_K5_alpha(hist_idx)
+    endif
+
+    ! Compute variable timestep Verlet coefficients
+    if (present(dt)) then
+      ! Use current timestep (dt) and previous timestep (xl%dt_history(1))
+      dt_n = dt                    ! Current timestep (input parameter)
+      dt_prev = xl%dt_history(1)   ! Previous timestep
+
+      ! Check if we have valid history (dt_prev > 0)
+      if (dt_prev > 1.0e-12_dp) then
+        ! Compute Verlet coefficients for variable timesteps
+        r = dt_n / dt_prev  ! Timestep ratio
+        P_n_coeff = 1.0_dp + r
+        P_n1_coeff = r
+
+        ! Compute kappa and alpha scaling factor
+        ! This comes from: 0.5 * dt_n * (dt_n + dt_{n-1})
+        kappa_alpha_scale = 0.5_dp * dt_n * (dt_n + dt_prev)
+      else
+        ! First step or dt_prev not set yet - use uniform coefficients
+        P_n_coeff = 2.0_dp
+        P_n1_coeff = 1.0_dp
+        kappa_alpha_scale = dt_n * dt_n  ! dt^2 for first step
+      endif
+    else
+      ! Uniform timestep: standard Verlet coefficients
+      P_n_coeff = 2.0_dp
+      P_n1_coeff = 1.0_dp
+      kappa_alpha_scale = 1.0_dp
     endif
 
     ! From developper's code
@@ -193,12 +456,48 @@ subroutine prg_xlbo_nint_kernel(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,kernel,
     !   alpha*(C0*n_0+C1*n_1+C2*n_2+C3*n_3+C4*n_4+C5*n_5+C6*n_6)
     !   n_6 = n_5; n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n
 
-    !call bml_print_matrix("ker",kernel,1,10,1,10)
-    !write(*,*)matmul(kernel,(charges-n))
-    !n = 2.0_dp*n_0 - n_1 + xl%cc*kappa*(charges-n) &
-    n = 2.0_dp*n_0 - n_1 - 1.0_dp*kappa*matmul(kernel,(charges-n)) &
-         + alpha*(C0*n_0+C1*n_1+C2*n_2+C3*n_3+C4*n_4+C5*n_5);
-    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n;
+    if (allow_adaptive_timestep) then
+          ! New method: Use variable timestep coefficients directly (no interpolation)
+
+          ! Get coefficient index from timestep history
+          hist_idx = get_K5_history_index(xl%dt_history(1:5))
+
+          ! Lookup coefficients for this specific history pattern
+          C0_use = XLBO_K5_C0(hist_idx)
+          C1_use = XLBO_K5_C1(hist_idx)
+          C2_use = XLBO_K5_C2(hist_idx)
+          C3_use = XLBO_K5_C3(hist_idx)
+          C4_use = XLBO_K5_C4(hist_idx)
+          C5_use = XLBO_K5_C5(hist_idx)
+          d_K_use = XLBO_K5_dK(hist_idx)
+
+          ! Use pattern-specific alpha value (capped at alpha_max/2 for stability)
+          alpha_use = XLBO_K5_alpha(hist_idx)
+
+          ! Integration using raw charges with variable coefficients and variable timestep Verlet
+          n = P_n_coeff*n_0 - P_n1_coeff*n_1 - kappa_alpha_scale*kappa_use*matmul(kernel,(charges-n)) &
+               + alpha_use*(C0_use*n_0+C1_use*n_1+C2_use*n_2+C3_use*n_3+C4_use*n_4+C5_use*n_5)
+
+    else
+      ! Integration using raw charges with variable timestep Verlet
+        n = P_n_coeff*n_0 - P_n1_coeff*n_1 - kappa_alpha_scale*kappa_use*matmul(kernel,(charges-n)) &
+             + alpha_use*(C0*n_0+C1*n_1+C2*n_2+C3*n_3+C4*n_4+C5*n_5)
+    endif
+
+    ! Shift history arrays
+    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n
+
+    ! Update timestep history if dt provided (dt is ratio: 1.0 for full step, 0.5 for half step)
+    ! History stores ratios that indicate spacing relative to user timestep
+    ! Interpolation always maps to fixed dt/2 uniform grid for stability
+    if (present(dt)) then
+      xl%dt_history(5) = xl%dt_history(4)
+      xl%dt_history(4) = xl%dt_history(3)
+      xl%dt_history(3) = xl%dt_history(2)
+      xl%dt_history(2) = xl%dt_history(1)
+      xl%dt_history(1) = dt  ! Store ratio (1.0 = full user timestep, 0.5 = half user timestep)
+      xl%nsteps_taken = xl%nsteps_taken + 1
+    endif
 
   end subroutine prg_xlbo_nint_kernel
 
@@ -207,15 +506,20 @@ end subroutine prg_xlbo_nint_kernel
   !! \brief In this case we are passing a premultiplied ressidue x kernel
   !! tis is done to avoid rank-specific multiplication within this routine.
   !! \param charges
-  subroutine prg_xlbo_nint_kernelTimesRes(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,kernelTimesRes,xl)
+  subroutine prg_xlbo_nint_kernelTimesRes(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,kernelTimesRes,xl,dt)
     implicit none
     real(dp), allocatable, intent(inout) :: n(:), n_0(:), n_1(:), n_2(:), n_3(:), n_4(:), n_5(:)
     real(dp), allocatable, intent(in) :: charges(:)
     real(dp), allocatable, intent(in) :: kernelTimesRes(:)
-    type(xlbo_type), intent(in) :: xl
-
+    type(xlbo_type), intent(inout) :: xl
     integer, intent(in) :: mdstep
+    real(dp), intent(in), optional :: dt
     integer :: nats
+    real(dp) :: kappa_use, alpha_use
+    logical :: allow_adaptive_timestep
+    integer :: hist_idx
+    real(dp) :: C0_use, C1_use, C2_use, C3_use, C4_use, C5_use, d_K_use
+    real(dp) :: dt_n, dt_prev, r, P_n_coeff, P_n1_coeff, kappa_alpha_scale
 
     nats = size(charges,dim=1)
 
@@ -237,11 +541,94 @@ subroutine prg_xlbo_nint_kernelTimesRes(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep
       n_3 = charges;
       n_4 = charges;
       n_5 = charges;
+      xl%dt_history = 0.0_dp
+      xl%nsteps_taken = 0
+    endif
+
+    ! Determine if we should allow adaptive time step
+    allow_adaptive_timestep = present(dt) .and. xl%nsteps_taken >= 6
+
+    kappa_use = kappa
+    alpha_use = alpha
+
+    ! Use pattern-specific alpha for early steps (during warmup before full history)
+    if (present(dt) .and. .not. allow_adaptive_timestep) then
+      ! Look up pattern-specific alpha based on current dt_history
+      hist_idx = get_K5_history_index(xl%dt_history(1:5))
+      alpha_use = XLBO_K5_alpha(hist_idx)
     endif
 
-    n = 2.0_dp*n_0 - n_1 - 1.0_dp*kappa*kernelTimesRes &
-         & + alpha*(C0*n_0+C1*n_1+C2*n_2+C3*n_3+C4*n_4+C5*n_5);
-    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n;
+    ! Compute variable timestep Verlet coefficients
+    if (present(dt)) then
+      ! Use current timestep (dt) and previous timestep (xl%dt_history(1))
+      dt_n = dt                    ! Current timestep (input parameter)
+      dt_prev = xl%dt_history(1)   ! Previous timestep
+
+      ! Check if we have valid history (dt_prev > 0)
+      if (dt_prev > 1.0e-12_dp) then
+        ! Compute Verlet coefficients for variable timesteps
+        r = dt_n / dt_prev  ! Timestep ratio
+        P_n_coeff = 1.0_dp + r
+        P_n1_coeff = r
+
+        ! Compute kappa and alpha scaling factor
+        ! This comes from: 0.5 * dt_n * (dt_n + dt_{n-1})
+        kappa_alpha_scale = 0.5_dp * dt_n * (dt_n + dt_prev)
+      else
+        ! First step or dt_prev not set yet - use uniform coefficients
+        P_n_coeff = 2.0_dp
+        P_n1_coeff = 1.0_dp
+        kappa_alpha_scale = dt_n * dt_n  ! dt^2 for first step
+      endif
+    else
+      ! Uniform timestep: standard Verlet coefficients
+      P_n_coeff = 2.0_dp
+      P_n1_coeff = 1.0_dp
+      kappa_alpha_scale = 1.0_dp
+    endif
+
+    if (allow_adaptive_timestep) then
+          ! New method: Use variable timestep coefficients directly (no interpolation)
+
+          ! Get coefficient index from timestep history
+          hist_idx = get_K5_history_index(xl%dt_history(1:5))
+
+          ! Lookup coefficients for this specific history pattern
+          C0_use = XLBO_K5_C0(hist_idx)
+          C1_use = XLBO_K5_C1(hist_idx)
+          C2_use = XLBO_K5_C2(hist_idx)
+          C3_use = XLBO_K5_C3(hist_idx)
+          C4_use = XLBO_K5_C4(hist_idx)
+          C5_use = XLBO_K5_C5(hist_idx)
+          d_K_use = XLBO_K5_dK(hist_idx)
+
+          ! Use pattern-specific alpha value (capped at alpha_max/2 for stability)
+          alpha_use = XLBO_K5_alpha(hist_idx)
+
+          ! Integration using raw charges with variable coefficients and variable timestep Verlet
+          n = P_n_coeff*n_0 - P_n1_coeff*n_1 - kappa_alpha_scale*kappa_use*kernelTimesRes &
+               & + alpha_use*(C0_use*n_0+C1_use*n_1+C2_use*n_2+C3_use*n_3+C4_use*n_4+C5_use*n_5)
+
+    else
+      ! Integration using raw charges with variable timestep Verlet
+        n = P_n_coeff*n_0 - P_n1_coeff*n_1 - kappa_alpha_scale*kappa_use*kernelTimesRes &
+             & + alpha_use*(C0*n_0+C1*n_1+C2*n_2+C3*n_3+C4*n_4+C5*n_5)
+    endif
+
+    ! Shift history arrays
+    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n
+
+    ! Update timestep history if dt provided (dt is ratio: 1.0 for full step, 0.5 for half step)
+    ! History stores ratios that indicate spacing relative to user timestep
+    ! Interpolation always maps to fixed dt/2 uniform grid for stability
+    if (present(dt)) then
+      xl%dt_history(5) = xl%dt_history(4)
+      xl%dt_history(4) = xl%dt_history(3)
+      xl%dt_history(3) = xl%dt_history(2)
+      xl%dt_history(2) = xl%dt_history(1)
+      xl%dt_history(1) = dt  ! Store ratio (1.0 = full user timestep, 0.5 = half user timestep)
+      xl%nsteps_taken = xl%nsteps_taken + 1
+    endif
 
   end subroutine prg_xlbo_nint_kernelTimesRes