Reproducibility of real-time time-dependent density functional theory calculations of electronic stopping power in warm dense matter Open Access

Warm dense matter (WDM) is a challenging regime where plasma models falter and experiments remain extremely difficult. Recent efforts increasingly leverage first-principles approaches based on density functional theory (DFT) to benchmark, inform, and improve more computationally efficient models for WDM and beyond.^1–7 The credibility of the underlying first-principles data then underpins the predictive power of the more efficient models. The scarcity of experimental avenues for model validation heightens the importance of building confidence in first-principles results by ensuring reproducibility, scrutinizing methodological choices, and quantifying uncertainties.

DFT calculations boast a high standard of reproducibility for static or equilibrium properties of materials⁸ despite differences among numerical implementations such as basis sets, pseudopotentials, and eigensolver algorithms. Simulations of excited electron dynamics using real-time time-dependent DFT (TDDFT) involve additional choices including the time-stepping algorithm, perturbation properties, simulation length, and techniques for extracting observable information from time series data. While some of these choices may be viewed as convergence parameters, high computational costs and non-trivial interactions among different options can make errors difficult to characterize. Furthermore, compared to DFT, TDDFT has fewer practitioners and widely available full-featured implementations. Thus, a similar standard of reproducibility remains to be established.

Here, we consider TDDFT calculations of electronic stopping power, the rate at which a moving ion loses energy to electrons in a medium. In inertial confinement fusion, this quantity controls an important contribution to the hotspot energy balance as fusion products re-deposit their kinetic energy to heat the fuel.⁹ Therefore, accurate stopping powers are important for correctly predicting ignition requirements.¹⁰ For ambient conditions, decades of experimental research have informed broad databases¹¹ that enable validation of TDDFT results for cold materials.^12–18 However, the accuracy of stopping power predictions in high-energy density systems on the way to ignition remains difficult to ascertain.

In this work, we evaluate the reproducibility of TDDFT stopping power calculations and systematically examine the influence of different numerical and methodological choices in the warm dense regime. While previous studies have validated and/or cross-verified final stopping power results,^19–23 those comparisons focus on a single post-processed quantity and leave open the possibility of fortuitous error cancelations that obscure the true degree of reproducibility. Published comparisons of time-resolved TDDFT data^12,24 remain rare but reveal significant discrepancies in the core contribution to stopping power that were attributed to details of the pseudopotentials.¹² 

We consider a simpler case without the sensitivities of core electrons: alpha-particle stopping power in warm dense hydrogen at a density of 1 g/cm³ and temperature of 2 eV. We further focus on a single alpha-particle velocity of 5 a.u., a choice that lies above the Bragg peak and involves both localized and delocalized excitations. This test system enables rigorous assessment of pseudopotential effects by benchmarking against all-electron calculations with bare Coulomb potentials. We compare both time-resolved data and average stopping powers computed using four different TDDFT codes: a custom extension^25,26 of the Vienna Ab initio Simulation Package (VASP),^27–29 Stochastic and Hybrid Representation Electronic structure by Density functional theory (SHRED),^30–32 Qball,^33,34 and Inq.²⁴ 

We demonstrate good agreement across all TDDFT implementations, especially when matching methodological choices as much as possible. Further analysis of sensitivities to different choices reveals which differences among computational approaches influence predicted stopping powers most. This work clarifies and deepens understanding of best practices in TDDFT stopping power calculations and supports ongoing efforts to quantify uncertainties in large-scale simulations of high-energy-density systems.⁷ 

In addition to the usual approximations, challenges, and limitations of DFT calculations—a microscopic simulation cell, a finite basis representing $ϕ j$⁠, the form of $V XC$⁠, pseudopotentials entering $V ext$⁠, and rigorous formalisms for computing physical observables—real-time TDDFT simulations also choose an initial condition and perturbation form, discretize time, and access a finite temporal range. In this work, all TDDFT calculations used a plane wave basis, began from Mermin–Kohn–Sham equilibrium states,³⁶ and held ionic velocities fixed. A plane wave cutoff energy of 2000 eV sufficed to converge results computed with bare Coulomb potentials (see Sec. II A for details), though some SHRED calculations used an overly conservative cutoff. Other aspects of the simulations varied as described in Subsections II A–II D and summarized in Table I.

We obtain stopping powers from stopping forces $F ( x ) = − F ( x ) · v ̂$ against the projectile's direction of motion $v ̂$ that depend implicitly on time through the distance x(t) traveled by the alpha particle. Different workflows rely on either the total force $F tot$ or the difference in the electron–ion contribution relative to a fixed electron density $Δ F ei$ (see Table I). In either case, a time average of the stopping forces gives the average stopping power accessed by experiments and entering into larger-scale hydrodynamic simulations.

Calculations aimed at making direct comparisons of time-dependent data began from the same atomic configuration for the hydrogen host material and initial position and velocity for the alpha projectile. We selected two setups reflecting different approaches commonly used by our groups: (i) a cubic simulation cell and an off-axis alpha particle trajectory optimized according to the methods developed by Kononov et al.¹² [see Fig. 1(a)] and (ii) an elongated simulation cell and the alpha particle traveling parallel to the long axis, following the approach of earlier work by Nichols et al.⁵ [see Fig. 1(b)]. For both setups, the atomic configuration was generated from separate ab initio molecular dynamics (MD) simulations as described in the  Appendix. Our comparisons often focus on 512-atom hydrogen configurations, though we also consider other simulation cell sizes where indicated.

Although typical VASP^27–29 calculations use the projector augmented-wave (PAW) method,³⁷ the light elements involved in this work also allow all-electron calculations with no pseudopotentials. By eliminating computations involving nonlocal projectors while only modestly increasing the required cutoff energy, bare Coulomb potentials significantly reduce computational costs for this system. Thus, we used bare Coulomb potentials to maximize both accuracy and efficiency of VASP calculations, and we evaluate the influence of the pseudopotential approximation in Sec. III A.

The adiabatic local density approximation (LDA)^38,39 was used in VASP calculations throughout this work, and we compare VASP results for LDA and Perdew–Burke–Ernzerhof (PBE)⁴⁰ functionals in Sec. III A. About 1.3 electronic bands per hydrogen atom were needed to capture Fermi occupations of at least $10 − 5$ at a temperature of 2 eV. For the $v = 5 a . u .$ case scrutinized throughout this work, plane wave cutoffs of 2000 and 2500 eV produced average stopping powers differing by less than 1% and time-dependent forces differing by less than 2.3 eV/Å after the 4 Å transient regime.

The TDDFT extension^25,26 of VASP uses the Crank–Nicolson (CN) time-stepping algorithm and solves the associated update equations using the generalized minimal residual (GMRES) method⁴¹ with a GMRES relative tolerance of $10 − 9$⁠. For fast projectiles with $v > 3.5 a . u .$⁠, we set the time step inversely proportional to the projectile velocity such that the alpha particle traveled about 0.02 Å within each time step. Slower projectiles required a shorter time step of 0.25 as to adequately resolve the electron dynamics. These choices converge average stopping power within 1% (6%) and time-dependent forces within 1 eV/Å (5 eV/Å) relative to a time step half as long for an alpha particle with 5 a.u. (1.5 a.u.) of velocity.

In addition to the deterministic Kohn–Sham formulation used by the other codes, SHRED^30–32 also includes capabilities for orbital-free, stochastic, and mixed stochastic-deterministic TDDFT. In this work, we use both the deterministic and mixed stochastic-deterministic TDDFT implementations within SHRED. The latter approach generates an initial finite-temperature Mermin density matrix from a subset of occupied KS orbitals and the complementary (orthogonal to the deterministic) stochastic density matrix. The stochastic density matrix has a compressed low-rank form similar to the KS density matrix; thus, stochastic “orbitals” are propagated in time following the same procedure as the KS orbitals. The full procedure is described in Ref. 31. For the 1024 (512) atom configuration, we use 384 (192) deterministic orbitals and 64 stochastic “orbitals,” with the stochastic representation accounting for about 25% of the electrons. Comparisons to deterministic TDDFT with 640 KS orbitals at low (1.5 a.u.) and high (5 a.u.) velocities show a difference of less than 1%.

The SHRED calculations used the PBE exchange-correlation functional⁴⁰ and either bare Coulomb potentials or Hartwigsen–Goedeker–Hutter (HGH) pseudopotentials,^42,43 which have an analytic form that allows direct control of the cutoff radius. The short iterative Lanczos (SIL) algorithm is used to propagate the Kohn–Sham wavefunctions and stochastic vectors.⁴⁴ The propagation is fully explicit with a time step of 0.121 as or less.

Rather than computing the total force on the projectile as in the other TDDFT implementations, SHRED calculations typically aim to isolate force contributions that do not average to zero over time. TDDFT studies commonly assume that the force induced by host ions $F ii$ has a zero time average because of symmetry along periodic trajectories or because nuclear stopping power is small in the velocity regime of interest. Similarly, the force induced by the initial electron density matrix—if it were held fixed as the projectile moved through it—would also average to 0 due to the same symmetry. This force, denoted as $F ei ( 0 )$⁠, is expected to dominate instantaneous forces at high velocities, when the stopping power is small and the projectile only slightly perturbs the initial electron density. Therefore, the SHRED implementation typically computes electron–ion force differences $Δ F ei = F ei − F ei ( 0 )$⁠.

The Qball^33,34 calculations used Hamann–Schluter–Chiang–Vanderbilt (HSCV) norm-conserving pseudopotentials⁴⁵ with outermost cutoff radii of 0.8 a.u. (0.5 a.u.) for the alpha particle (hydrogen ions). The electronic states evolved over time according to the enforced time-reversal symmetry (ETRS) algorithm⁴⁶ with fourth-order Taylor expansions for the exponentials.³⁴ Rather than solving the nonlinear update equations to within a given tolerance, the ETRS implementation performs a two-step iteration that first estimates the electron density at the next time step, $n ( r , t + Δ t )$⁠, by evolving the KS states under H(t). The implicit portion of the propagator then uses the estimated $n ( r , t + Δ t )$ to construct $H ( t + Δ t )$⁠, with further self-consistent iteration not necessary for sufficiently small time steps.³⁴ 

Exchange and correlation was treated at the LDA^38,39 level, and other numerical parameters were identical to those listed in Sec. II A for VASP. Independent convergence tests verified that the cutoff energy and time step also suffice for the different time stepper and pseudopotentials used within Qball. For an alpha particle with 5 a.u. of velocity, increasing the cutoff energy to 2500 eV or halving the time step changes the average stopping power by less than 0.1% and alters time-dependent forces by less than 0.05 eV/Å.

The Inq²⁴ simulations follow very similar methodology as described for Qball in Sec. II C. However, we used the default pseudopotential set within Inq: optimized norm-conserving Vanderbilt (ONCV) pseudopotentials⁴⁷ with outermost cutoff radii of 1.25 a.u. (1.0 a.u.) for the alpha particle (hydrogen ions). Since these pseudopotentials are even softer than those used in the Qball calculations, we expect excellent convergence at the same cutoff energy. Also, the Inq ETRS implementation performs up to five self-consistent iterations to obtain $n ( r , t + Δ t )$ and $H ( t + Δ t )$⁠, further enhancing convergence with respect to time step.

We present the average stopping powers computed by different TDDFT codes and methodologies (see Table I) in Fig. 2. Overall, we find reasonable agreement, though some of the discrepancies exceed the stringent numerical convergence achieved within each set of calculations (see Sec. II). In particular, we see excellent agreement between VASP and SHRED results when both sets of calculations use bare Coulomb potentials and the same 512-atom simulation setup [blue squares and purple diamonds in Fig. 2(b)], but modest deviations arise when using pseudopotentials and/or changing the simulation cell shape and alpha particle trajectory. The spread among datasets in Fig. 2 represents overall uncertainties in the computed stopping powers.

Since stopping power is a highly integrated quantity, we also compare the instantaneous friction force experienced by the projectile along its path for the case of $v = 5 a . u .$ in Fig. 3, where all calculations used the same off-axis alpha particle trajectory through the same 512-atom hydrogen configuration within the cubic simulation setup [see Fig. 1(a)], i.e., the same setup used for Fig. 2(b). Overall, we find very good agreement in these stopping forces, particularly after an initial transient regime during the first few Å. However, small relative differences in the stopping forces accumulate into noticeable differences in average stopping power, which comes out to 28.8 eV/Å from VASP, 31.4 eV/Å from Qball, 30.9 eV/Å from SHRED, and 23.8 eV/Å from Inq when using the 512-atom cubic setup, identical post-processing methods to extract $S fit$ from $F tot$⁠, and other choices as listed in Table I.

In Subsections III A–III E, we examine the extent to which different methodological variations contribute to the modest discrepancies occurring in Figs. 2 and 3. Among the different choices listed in Table I, we show in Sec. III A that pseudopotentials can have a large and persistent effect on time-dependent stopping forces. Meanwhile, the XC functional (LDA or PBE) and the initial condition (including or excluding the alpha-particle projectile) primarily alter forces during the first few Å traversed and do not significantly affect average stopping powers for this system (see Secs. III A and III B). Sections III C and III D demonstrate close agreement between average stopping powers estimated from different possible force observables and post-processing techniques. Finally, Sec. III E reveals contrasting behavior of finite-size errors within the cubic and elongated simulation setups.

We begin by considering differences among the methods described in Sec. II that alter the underlying Hamiltonian [Eq. (3)]: the electron–ion potentials determining $V ext$ and the exchange-correlation (XC) functional $V XC$⁠. We expect the XC functional to have only a minor influence because LDA and PBE stopping powers differed by less than 0.5% in an earlier study of a metallic system.¹² On the other hand, it is well known that pseudopotentials can cause severe underestimation of high-velocity stopping powers by neglecting contributions from pseudized core electrons.¹⁵ Recent work also showed that softer pseudopotentials limit energy transfer during close collisions with host nuclei, reducing stopping power predictions by around 10% even for the same number of explicitly simulated electrons.¹² Although our light-ion test system contains no core electrons and all the calculations treated all of the electrons explicitly, pseudopotentials can still influence processes involving electrons localized around the nuclei.

In Fig. 4, we compare VASP predictions obtained with LDA and PBE XC functionals and using bare Coulomb or PAW potentials for both the projectile and host ions. We examine only the electron–ion contribution to the stopping forces because the ion–ion contribution dominating total forces remains identical across all cases. We find that LDA and PBE indeed produce very similar forces along the alpha particle's direction of motion, differing by at most 2.4 eV/Å after the alpha travels 1 Å within VASP. Although these deviations are significant relative to the average stopping power of 28.8 eV/Å for this case, their time-averages amount to less than 0.5% of the average stopping power. Therefore, we conclude that the choice of XC functional does not appreciably contribute to discrepancies among different codes in this work.

On the other hand, pseudopotentials can have a much stronger influence on stopping power: VASP calculations using PAW predict systematically weaker stopping forces than obtained with bare Coulomb potentials. These deviations persist throughout the simulation and lead to a 13% lower average stopping power with PAW than with bare Coulomb potentials. However, contrary to earlier findings for proton stopping in ambient aluminum,¹² we find less than 0.1% sensitivity to the hardness of the hydrogen PAW, suggesting that the projectile pseudopotential has a stronger influence in this case.

Similarly, different norm-conserving pseudopotentials within Qball or Inq produce significantly different stopping powers, either overestimating or underestimating by up to 18% relative to VASP with bare Coulomb potentials. Therefore, we attribute the larger deviations of Qball and Inq data from other codes in Figs. 2 and 3 to the pseudopotential approximation. Surprisingly, stopping powers predicted from SHRED appear less sensitive to pseudopotentials, with only a 2% difference between results computed using bare Coulomb potentials and HGH pseudopotentials with cutoff radii of 1 a.u. for both species.

In addition to the form of the Hamiltonian, the solution to the TDKS equations [Eqs. (1) and (2)] also depends on the initial KS orbitals $ϕ j ( r , t = 0 )$⁠. In particular, the density of states (DOS) corresponding to the initial orbitals relates to the distribution of electronic transitions that the projectile could excite. For this reason, the DOS enters some stopping power models based on the average-atom method,⁴ and close agreement for the initial DOS helps ensure agreement among stopping powers.

All the TDDFT simulations began from Mermin–KS equilibrium states,³⁶ and Fig. 5 demonstrates that the Mermin-DFT implementations predict very similar electronic properties for the static, unperturbed hydrogen host even under different pseudopotentials and XC functionals (see Table I). The largest deviations among the DOS in Fig. 5 amount to less than 1%. Furthermore, including the alpha particle within the Mermin-DFT calculation has only a minor influence on the DOS of the initial system, and we expect that this effect will further diminish with increasing system size.

Although we find excellent agreement for initial electronic properties, different treatments of the projectile within the initial condition determine initial screening of the alpha particle and thus influence stopping forces during the TDDFT simulations. The SHRED calculations exclude the projectile from the initial condition and insert it at t = 0 within the TDDFT calculation, thereby simulating an initially fully ionized alpha particle. On the other hand, VASP calculations typically include the projectile in the initial condition, resulting in an initially partially neutralized ion. Previous Qball studies have shown that the projectile's initial charge state can affect results for high-Z projectiles¹³ and pre-equilibrium behavior in 2D materials,⁴⁸ but light-ion stopping powers in bulk materials do not depend strongly on the initial condition¹⁵ because the projectile's charge state quickly equilibrates after traversing a few Ås.⁴⁹ 

Here, we revisit the sensitivity of stopping powers to the initial condition in the warm dense regime. Consistent with earlier results for ambient materials, we find that the VASP stopping forces do not depend strongly on whether the alpha particle was included in the initial condition (see Fig. 4). After a transient regime within the first 4 Å of the alpha particle's path—during which the ion dynamically captures or sheds electrons to reach an equilibrium charge state—the initial condition changes stopping forces by at most 3.8 eV/Å. On average, these deviations again amount to less than 0.5% of the average stopping power and thus do not meaningfully affect the code comparisons in this work.

In Fig. 6, we show that although the magnitude and time-dependent behavior of each force contribution exhibit stark differences, stopping powers inferred from $F tot , F ei$⁠, and $Δ F ei$ agree closely. Although the ion–ion contribution $F ii$ dominates the absolute scale of the total force [cf. Figs. 6(a) and 6(b)], its time average quickly decays toward 0 and amounts to 2.6% of the converged total stopping power [see Fig. 6(d)]. Similarly, including the force due to the unperturbed electron density $F ei ( 0 )$ changes the extracted stopping power by only 0.6% within VASP and 0.2% within SHRED.

Contributions from $F ii$ and $F ei ( 0 )$ do not depend on projectile velocity, but become more significant relative to lower stopping powers in the high-velocity regime. For the fastest alpha particle considered in this work, $v = 8 a . u .$⁠, including $F ii$ and $F ei ( 0 )$ changes the stopping power by 6% and 1.4%, respectively. Thus, excluding these terms becomes more important for precise calculations of electronic stopping power at high velocities.

Given a choice of observable, different analysis methods are possible for estimating an average stopping power [see Eqs. (5) and (6)]. Reference 12 previously showed that when analyzing the total force, $S fit$ behaves more smoothly and converges faster over time than $S avg$⁠. Here, we revisit the behavior of both methods when analyzing $F ei$ and $Δ F ei$⁠.

In Fig. 7(a), the stopping work computed from $F tot$ contains sharp features arising from strong occasional ion–ion force pulses [see Fig. 6(a)]. Consistent with Ref. 12, similar sharp features appear in $S avg$ because this estimate depends upon only $W ( x f )$ and $W ( x 0 )$⁠, whereas $S fit$ provides a smoother estimate because it is informed by a range of W(x) data [see Fig. 7(b)]. Excluding ion–ion forces and/or the force induced by the equilibrium density—i.e., analyzing $F ei$ or $Δ F ei$—dramatically improves the behavior of $S avg$⁠, but the qualitative behavior of $S fit$ is less sensitive to the choice of force and still maintains somewhat smoother evolution in all cases. All combinations of stopping force and averaging method agree within 4% by a projectile displacement of 80 Å, though these choices would contribute to uncertainties for shorter simulation lengths and may become more important for higher-Z host materials with core electrons contributing sharp features to $F ei$ and/or higher projectile velocities where $F ei ( 0 ) ≫ Δ F ei$⁠.

Finally, we compare data obtained from the two different simulation setups: cubic supercell with off-axis trajectory and elongated supercell with on-axis trajectory. The two approaches differ in terms of the environments sampled by the alpha particle and the behavior of finite-size errors, and these factors influence how well the computed stopping powers reflect practically meaningful values independently of the approximations and choices discussed above. To disentangle the two effects, we consider two types of elongated supercells: 256–1024-atom cells constructed by replicating a thermalized 256-atom configuration and a fully thermalized 1024-atom cell. The replicated cells provide information about finite-size errors under fixed environment sampling, while the thermalized cell represents a more converged production setup.

Naively, stopping powers in a hydrogen plasma might not be expected to depend strongly on environment sampling because the density of states resembles that of a Fermi gas (see Fig. 5). Sampling considerations become important in systems and regimes where core electrons (or electrons with strong spatial nonuniformities) significantly contribute to stopping powers. For example, earlier work^12,14,15 found that obtaining accurate stopping powers beyond the Bragg peak in crystalline aluminum requires that the projectile representatively sample close collisions with host nuclei and the associated 2s and 2p orbital excitations. In that system, no single on-axis trajectory achieves sufficiently representative sampling and only appropriately chosen off-axis trajectories produce stopping powers in agreement with empirical data.

Despite the nearly ideal DOS and lack of core electrons, the equilibrium electron density of the warm dense hydrogen host contains spatial variations spanning more than two orders of magnitude. This nonuniformity causes surprisingly large sensitivities to the alpha particle trajectory. Stopping powers computed using ten random on-axis trajectories through elongated supercells, i.e., using different initial positions for an alpha particle traveling in the long direction, deviated from their average by as much as 20% (see Fig. 8). While this trajectory-dependence appears weaker than in the aluminum case,^12,14,15 it remains a significant source of uncertainty in these TDDFT calculations.

The length of the simulation cell fundamentally limits the sampling fidelity achievable by a single on-axis trajectory, as the environment would become periodic after one pass through the cell. In contrast, an off-axis trajectory can continue sampling new environments and improving the sampling fidelity during each pass through the periodic cell while still avoiding its own wake.¹² Using an ensemble of several trajectories (either on-axis or off-axis) can further improve sampling fidelity.⁵⁰ 

The two simulation setups also differ in the behavior of finite-size errors, which TDDFT stopping power studies often attribute to the neglect of long-wavelength plasmonic excitations.⁵¹ The elongated simulation cell supports longer wavelength plasmons along the projectile's direction of motion, reducing this contribution to finite-size error. However, as the projectile crosses periodic boundary conditions, artificial interactions with earlier excitations cause significant deviations from the plasmonic model of finite-size errors.¹² These so-called ouroboros effects would be more severe for a long on-axis trajectory, where the projectile traverses the same path during each pass through the simulation cell.

Indeed, Fig. 8 shows only modest (up to 5%) stopping power sensitivity to supercell size for short on-axis trajectories that complete less than a single pass through the elongated cell (orange symbols), both for individual trajectories (empty symbols) and their ensemble average (filled symbols). We note that this observation only holds for supercells long enough that each trajectory begins to approach a steady-state stopping power within the single pass; we find much larger errors of up to 50% (23%) for individual trajectories (ensemble average) in a shorter, 256-atom supercell. In contrast, using a long on-axis trajectory (red right-filled symbols), which achieves steady-state stopping at the expense of making multiple passes over the same path through the elongated cell, incurs large and unsystematic finite-size errors of about 20% even for cell sizes that are well-converged over a single pass. We also find significant finite-size errors in the cubic setup with an optimized off-axis trajectory (green top-filled symbols), with a 15% difference between 512-atom and 1000-atom stopping powers.

Overall, we find good agreement within about 4% when appropriately controlling both environment sampling and finite-size errors using different strategies. Depending on the approach, both effects can significantly influence the accuracy of stopping power predictions. Long on-axis trajectories simultaneously suffer from strong finite-size effects and limited sampling fidelity. While using short on-axis trajectories through elongated supercells converges very quickly with respect to simulation cell size, achieving representative environment sampling with this method requires an ensemble of trajectories. Conversely, long off-axis trajectories can be optimized to achieve excellent sampling fidelity, but require larger supercells to converge finite-size errors.

In this study, both of the simulation setups that efficiently control environment sampling and finite-size errors—a single, optimized off-axis trajectory through a 1000-atom simulation cell or an ensemble of short, on-axis trajectories through a 512-atom simulation cell—incur comparable total computational costs. Future work may investigate the ability of a compromising approach relying on intermediate-length off-axis trajectories to better balance computational costs in reducing both types of errors.

Through detailed comparisons across four independent real-time TDDFT implementations, this work establishes reproducibility standards and quantifies uncertainties arising from different possible methodological choices in first-principles calculations of electronic stopping power. For the case of alpha particles in warm dense hydrogen, the choice of XC functional (adiabatic LDA or PBE) and initial condition (including or excluding the projectile in the Mermin–KS equilibrium states) only negligibly influence average stopping power because these options only weakly alter underlying time-dependent forces beyond a short transient regime. Similarly, for sufficiently long simulation times, we find close agreement among different post-processing techniques involving various possible force contributions.

The dominant sources of error in TDDFT stopping powers calculations typically arise from pseudopotentials, poor sampling of the environment by the projectile, and/or finite simulation cells. Even in the all-electron simulations considered here, pseudopotentials can systematically underestimate stopping powers by 13% compared to bare Coulomb potentials. The complex nature of environment sampling and finite-size errors has inspired different mitigation strategies, and direct comparison of their convergence behavior suggests that both effects remain significant even within this relatively simple test system. Infrared truncation of plasmonic excitations, limited sampling of different plasma environments, and artificial interaction with the projectile's wake would still limit accuracy even if exact simulation of the many-body electron dynamics⁵² became feasible, as obtaining experimentally relevant values requires expensive averages over multiple trajectories, large supercells containing thousands of particles, and/or complicated techniques for optimizing trajectories.

This work informs ongoing efforts to benchmark and improve more efficient stopping power models using first-principles calculations.^4,5,7,23,31 Our detailed sensitivity analysis not only probes the reliability of TDDFT results, but also enables careful interpretation of discrepancies within model comparisons. For example, using TDDFT simulations to resolve modest differences among stopping powers predicted by linear-response treatments^4,23,31 may require very strict convergence and rigorous control of methodological choices. Meanwhile, our work suggests that larger discrepancies between TDDFT data and certain average-atom models⁷ (around a factor of 2 in some cases) are relatively robust to details of the TDDFT calculations and therefore warrant increased scrutiny in future work.

We thank the organizers of the 2nd Charged-Particle Transport Coefficient Comparison Workshop—including Lucas J. Stanek, Brian M. Haines, Stephanie B. Hansen, Patrick F. Knapp, Michael S. Murillo, Liam G. Stanton, and Heather D. Whitley—for catalyzing this work. We are also grateful to Yifan Yao for technical advice, Joshua P. Townsend for sharing equilibrium H configurations, and Heath Hanshaw for pre-publication review.

A.K. and A.D.B. were partially supported by the U.S. Department of Energy Science Campaign 1 and Sandia National Laboratories' Laboratory Directed Research and Development (LDRD) Project No. 233196. A.J.W. was supported by the LDRD program of Los Alamos National Laboratory (LANL), under Project Nos. 20210233ER and 20230323ER, and U.S. Department of Energy Science Campaign 4. This research used computing resources provided by the LANL Institutional Computing and Advanced Scientific Computing programs. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). K.A.N. and S.X.H. were supported by the Department of Energy (National Nuclear Security Administration) University of Rochester “National Inertial Confinement Program” under Award No. DE-NA0004144.

Sandia National Laboratories is a multi-mission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC (NTESS), a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration (DOE/NNSA) under Contract No. DE-NA0003525. This written work is authored by an employee of NTESS. The employee, not NTESS, owns the right, title and interest in and to the written work and is responsible for its contents. Any subjective views or opinions that might be expressed in the written work do not necessarily represent the views of the U.S. Government. The publisher acknowledges that the U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this written work or allow others to do so, for U.S. Government purposes. The DOE will provide public access to results of federally sponsored research in accordance with the DOE Public Access Plan.

The authors have no conflicts to disclose.

Alina Kononov: Conceptualization (equal); Data curation (equal); Formal analysis (lead); Writing – original draft (lead); Writing – review & editing (equal). Alexander J. White: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Software (equal); Writing – original draft (supporting); Writing – review & editing (equal). Katarina A. Nichols: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Writing – review & editing (equal). S. X. Hu: Conceptualization (equal); Data curation (equal); Writing – original draft (supporting); Writing – review & editing (equal). Andrew D. Baczewski: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Software (equal); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In the main text, we consider cubic and elongated simulation cells containing a hydrogen configuration in thermal equilibrium (see Fig. 1). The cubic simulation cells were prepared with VASP^27–29 as described in Sec. 1 of the  Appendix, while the elongated simulation cells were prepared with SHRED^30–32 as described in Sec. 2 of the  Appendix.

The finite-temperature stopping power calculations used thermally equilibrated ionic configurations generated through separate ab initio MD simulations. The VASP MD simulations started from a disordered geometry where each hydrogen atom was randomly displaced from a periodic arrangement by up to 0.2 Å in each Cartesian direction. Then, the ionic system evolved under Hellmann-Feynman forces from the electronic system and the Nosé–Hoover thermostat⁵³ with a Nosé mass of 1 and time step of 0.2 fs. The equilibrated geometry after 1 ps of simulated time then served as the atomic configuration for TDDFT stopping power simulations.

The MD simulations treated the electronic system with many of the same choices and parameters listed in Sec. II A. However, a reduced cutoff energy of 1000 eV sufficed to converge the average pressure and the average free energy of the electronic system to within 3% (0.7%) of the 2000 eV cutoff used in TDDFT calculations. To further reduce computational costs, the MD simulations relaxed the electronic self-consistency criterion to energy convergence within 1 meV.

To generate the replicated configurations referenced in Sec. III E, we used ab initio MD with the deterministic DFT capability in SHRED. We started with random initial ion positions in a supercell containing 256 atoms. Their initial velocities are given by a Maxwell–Boltzmann distribution corresponding to the given temperature. In each MD step, SHRED started with a self-consistent-field (SCF) calculation to solve for the equilibrium electron density. Once a converged electron density is obtained, the electronic force is derived from the Hellman–Feynman theorem, which along with ionic forces drives the motion of classical ions in each MD step. These SCF-MD cycles are repeated for thousands of steps to sample different ionic configurations for a real warm-dense matter system.

Specifically for the warm dense hydrogen problem considered here, we used 480 bands, the PBE exchange-correlation functional,⁴⁰ and the Hartwigsen–Goedeker–Hutter (HGH) pseudopotential^42,43 for the electron–ion interaction in SHRED SCF calculations. The MD time step was 0.1 fs and the thermostat enforced isokinetic velocity scaling at every MD step. Overall, the SHRED MD run gave a converged pressure of $4.82 ± 0.066$ Mbar for the hydrogen system. A long trajectory from such an MD run can provide multiple uncorrelated snapshots for TDDFT stopping power calculations, though this work used only a single equilibrated MD configuration for each set of TDDFT simulations.

A similar procedure was used for the fully thermalized 1024-atom configuration, but the orbital-free DFT (OF-DFT) capability in SHRED was used for this case. The kinetic energy functional by Thomas, Fermi, and Perrot⁵⁴ and the nonlinear conjugate gradient algorithm outlined by Jiang and Yang⁵⁵ were used to solve for the electron density. The snapshot used for TDDFT stopping power calculations was taken after 1 ps (with an MD time step of about 0.25 fs) of OF-DFT-MD simulation time. The supercell geometry was rectilinear with an aspect ratio of 4:1:1 and 256 × 64 × 64 real-space grid points (678 eV cutoff energy). A pressure of $5.04 ± 0.0015$ Mbar was calculated after discarding the first 0.2 ps for equilibration.