MemoryKernelFitting.jl

FDT-Safe Prony-Series Fits for Memory-Dependent Electronic Friction

---

MemoryKernelFitting.jl fits underdamped Prony-series memory kernels to configuration-resolved spectral data. The resulting extended-variable stochastic differential equation reproduces memory-dependent, configuration-dependent electronic friction while preserving the fluctuation–dissipation structure required to sample the correct Gibbs measure.

The input is a real-valued spectral density $\Lambda(x,\omega)$ on a configuration grid $x$ and frequency window $\omega \in [0,\omega_{\max})$. The output is a compact set of real and complex-conjugate modes from which the quasi-Markovian blocks $\widetilde{\boldsymbol\Gamma}(x)$ and $\widetilde{\boldsymbol\Sigma}(x)$ can be constructed.

$$\begin{aligned} \Lambda_\theta(x,\omega) &= a_0(x) + \sum_{k=1}^{K_R} a_k(x)\frac{\lambda_k}{\lambda_k^2+\omega^2} \\\ &\quad + \sum_{j=1}^{K_C} A_j(x)L(\omega;\alpha_j,\beta_j), \end{aligned}$$ $$\begin{aligned} K(x,\tau) &= 2a_0(x)\delta(\tau) + \sum_k a_k(x)e^{-\lambda_k\tau} \\\ &\quad + \sum_j A_j(x)e^{-\alpha_j\tau}\cos(\beta_j\tau). \end{aligned}$$

Features

Component	Description
Spectral model	Zero-frequency white part, overdamped Lorentzians, and mirrored underdamped Lorentzian pairs through spectral_density
Time-domain kernel	Matching exponential and damped-oscillatory Prony series through memory_kernel
Basis functions	lorentzian and mirrored_lorentzian, implemented as an explicit one-sided cosine-transform pair
Resonant modes	Real $2\times2$ rotation blocks for complex-conjugate poles $\alpha\pm i\beta$, allowing peaks at nonzero frequency
FDT-safe parameterisation	Non-negative amplitudes represented structurally as squares, with stable positive decay rates
Pole ladders	Log-spaced and uniform fixed ladders with the identifiability box of the frequency window, through log_ladder and uniform_ladder
Configuration dependence	GridAmplitudes for the convex stage, plus polynomial and warped-Chebyshev bases for the smooth stage, behind one AmplitudeModel interface
QGLE readout	Extended-variable blocks and the closed-form block exponential through qgle_blocks and exp_gamma22
Spectral data	Atomic-unit loading of the configuration-resolved spectra, graded-grid quadrature weights, and logarithmic frequency subsampling
Validation	Tests for symmetry, resonance onset, Markovian limits, dimension checks, numerical kernel–spectrum consistency, and fluctuation–dissipation structure
Fitting pipeline	Weighted objective and fixed-ladder per-configuration NNLS implemented; the coupled roughness solve and optional diagnostic VARPRO remain planned

Installation

Prerequisites: Julia ≥ 1.10 (download). The package is in early development and is not yet registered in the Julia General registry.

Setup:

julia> using Pkg
julia> Pkg.develop(url="https://github.com/Louhokseson/MemoryKernelFitting.jl")
julia> Pkg.test("MemoryKernelFitting")

For local development from a clone:

julia --project=. -e 'using Pkg; Pkg.instantiate(); Pkg.test()'

Quick Start

using MemoryKernelFitting

a₀ = 0.0                                      # white / instantaneous part
a, λ = [0.75], [10.0]                         # overdamped modes
A, α, β = [1.5, 0.4], [0.5, 1.2], [3.0, 7.0] # underdamped pairs

Λω = spectral_density(3.0, a₀, a, λ, A, α, β)
Kτ = memory_kernel(0.2, a, λ, A, α, β)

# Markovian anchor: the zero-frequency friction used by standard MDEF
Λ₀ = spectral_density(0.0, a₀, a, λ, A, α, β)

The white contribution $2a_0\delta(\tau)$ is intentionally excluded from memory_kernel; it enters the dynamics as instantaneous friction rather than through the regular memory convolution.

Development Status

The model, data, objective, and fixed-ladder stage-1 fit are implemented and tested: spectral and time-domain evaluation, the two basis functions, pole ladders, the three configuration-dependent amplitude representations, QGLE readout, atomic-unit spectral data loading, the weighted objective, and independent per-configuration NNLS.

The remaining fitting work is the $\mu$-coupled roughness solve and the stage-2 coefficient fit that compresses the grid amplitudes into a smooth representation. Outer variable projection is an optional pole diagnostic rather than part of the fixed-ladder production path.

Design decisions with lasting consequences are recorded in docs/design/, and progress is logged in docs/DEVLOG.md.

Repository Layout

Path	Contents
src/	Package implementation and public model-evaluation API
test/	Unit and numerical cosine-transform consistency tests
docs/PronySeries/	QGLE derivation, optimisation problem, and theory handover notes
assets/	Static and animated package logos plus the animation generator
data/	Local spectral and simulation data; large outputs are not tracked
plots/	Generated visualisations and fitting diagnostics
notebooks/	Exploratory analysis notebooks

The Fitting Model

The mirrored underdamped basis is

$$\begin{aligned} L(\omega;\alpha,\beta) &= \frac{\alpha}{2} \left[ \frac{1}{\alpha^2+(\omega-\beta)^2} + \frac{1}{\alpha^2+(\omega+\beta)^2} \right], \\\ a_0,\,a_k,\,A_j &\ge 0. \end{aligned}$$

The underdamped pairs are mandatory, not a refinement. The target spectra are not monotone in $\omega$: they carry a resonance at $\omega_{\mathrm{pk}}\approx1.02$–$1.15,|h(x)|$ that sweeps from zero up to $9.4,\mathrm{eV}$ as the electronic level crosses the Fermi energy. A sum of zero-centred Lorentzians is strictly decreasing and cannot produce such a peak; complex-conjugate pole pairs can.

The mirror in $L$ follows from realness rather than modelling preference. Because $\boldsymbol\Gamma_{2,2}$ is real, its complex eigenvalues $\alpha\pm i\beta$ occur in conjugate pairs and generate Lorentzians at $\omega=\pm\beta$ with half-width $\alpha$. The pair is a genuine resonance when $\beta>\alpha/\sqrt{3}$; below that threshold it merges continuously into a monotone, zero-centred response.

FDT-Safe Construction

Set $m=K_R+2K_C$ and retain real rotation blocks:

$$\begin{aligned} \boldsymbol\Gamma_{2,2} &= \mathrm{blockdiag} \left( \lambda_1,\ldots,\lambda_{K_R}, R(\alpha_1,\beta_1),\ldots \right), \\\ R(\alpha,\beta) &= \begin{pmatrix} \alpha & \beta\\\ -\beta & \alpha \end{pmatrix}. \end{aligned}$$ $$\begin{aligned} \widetilde{\boldsymbol\Gamma}_{1,1}(x) &=a_0(x), \\\ \boldsymbol g(x) &= \left( \sqrt{a_1},\ldots,\sqrt{A_1},0,\ldots \right)^T, \\\ \widetilde{\boldsymbol\Gamma}_{2,1} &=\boldsymbol g, \qquad \widetilde{\boldsymbol\Gamma}_{1,2}=-\boldsymbol g^T. \end{aligned}$$

Two structural consequences make the fit well posed:

• Amplitudes are squares. The fluctuation–dissipation theorem reduces the coupling to one fitted function $\boldsymbol g(x)$ and enforces non-negative $a_k$ and $A_j$. A negative fitted amplitude is therefore a hard signal that no Gibbs-consistent QGLE within this ansatz reproduces the data.
• $\boldsymbol Q=\boldsymbol I$ survives. Each real rotation block gives $\boldsymbol g^T e^{-\tau R}\boldsymbol g =|\boldsymbol g|^2e^{-\alpha\tau}\cos(\beta\tau)$ and $R+R^T=2\alpha\boldsymbol I$, so $\widetilde{\boldsymbol\Sigma}2=\sqrt{2\alpha},\boldsymbol I_2$. The stability requirement is $\boldsymbol\Gamma{2,2}\boldsymbol Q +\boldsymbol Q\boldsymbol\Gamma_{2,2}^T\succ0$, not symmetry of the drift matrix.

Fitting Strategy

The model is linear in the amplitudes and nonlinear only in the poles, so variable projection separates the optimisation into an inner convex non-negative least-squares problem and an outer low-dimensional pole search. Relative residual weights are needed because $\Lambda$ spans approximately $10^{-12}$ to $10^4$; configuration-space quadrature weights correct a grid where 83% of the columns cover only 6.7% of the range; and roughness regularisation prevents near-degenerate poles from exchanging amplitude arbitrarily.

Rates are shared across configurations while amplitudes vary with $x$. This is the regime in which the two-time path-ordered kernel collapses exactly to a Prony series in the lag. Allowing poles to track $x$ instead is an adiabatic approximation that fails most strongly near the crossing, where $\tau_{\mathrm{mem}}/\tau_{\mathrm{conf}}\approx34$ even at 300 K. A tiling argument gives $K_C=30$ pairs ($m=60$), but the measured fixed-ladder fit plateaus near $K_C\approx10$ ($m=20$); the final ladder should therefore be chosen from measured accuracy and integration cost rather than the tiling estimate alone. Once the ladder is fixed, the outer search disappears and the amplitude fit is a convex NNLS problem with a global optimum.

Documentation

The detailed theory and optimisation specification live in docs/PronySeries/:

File	Contents
Prony series (corrected).md	Canonical QGLE derivation, two-time kernel, and FDT pairing
Prony_fit_optimisation_problem.md	Full fitting objective, calibration, and underdamped-series construction
HANDOVER.md	Index and development handover

Citation

This package builds on the quasi-Markovian GLE and underdamped memory-kernel constructions described in:

• Sachs, PhD thesis (University of Edinburgh, 2017), Eq. (3.7) and §4.7.2.
• Leimkuhler and Sachs (2019), Ergodic Properties of Quasi-Markovian GLEs with Configuration Dependent Noise.
• Stella, Lorenz, and Kantorovich, Physical Review B 89, 134303 (2014).
• Lei, Baker, and Li, PNAS 113, 14183 (2016).

If you use MemoryKernelFitting.jl in published work, cite the relevant theory above and this repository:

@software{MemoryKernelFitting_jl,
  author = {Lu, Xuexun},
  title  = {{MemoryKernelFitting.jl}: FDT-safe Prony-series fits for memory-dependent electronic friction},
  year   = {2026},
  url    = {https://github.com/Louhokseson/MemoryKernelFitting.jl}
}

License

This project is licensed under the MIT License. See the full text for details.

---

Xuexun Lu (Hokseon)
PhD Candidate, The Maurer Computational Surface Science Group
The University of Warwick, UK

Copyright © 2026 Xuexun Lu. Licensed under the MIT License.