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

Real-time time-dependent density functional theory (TDDFT) is widely considered to be the most accurate available method for calculating electronic stopping powers from first principles, but there have been relatively few assessments of the consistency of its predictions across different implementations. This problem is particularly acute in the warm dense regime, where computational costs are high and experimental validation is rare and resource intensive. We report a comprehensive cross-verification of stopping power calculations in conditions relevant to inertial confinement fusion conducted using four different TDDFT implementations. We find excellent agreement among both the post-processed stopping powers and relevant time-resolved quantities for alpha particles in warm dense hydrogen. We also analyze sensitivities to a wide range of methodological details, including the exchange-correlation model, pseudopotentials, initial conditions, observable from which the stopping power is extracted, averaging procedures, projectile trajectory, and finite-size effects. We show that among these details, pseudopotentials, trajectory-dependence, and finite-size effects have the strongest influence, and we discuss different strategies for controlling the latter two considerations.


I. INTRODUCTION
Warm dense matter (WDM) is a challenging regime where plasma models falter and experiments remain extremely difficult.2][3][4][5][6][7] The credibility of the underlying firstprinciples 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 8 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 timestepping 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.
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 hot spot energy balance as fusion products redeposit their kinetic energy to heat the fuel. 9Therefore, accurate stopping powers are important for correctly predicting ignition requirements. 10 this work, we evaluate the reproducibility of TDDFT stopping power calculations and systematically examine the influence of different numerical and methodological choices.While previous studies have validated and/or cross-verified final stopping power results, [11][12][13][14][15] those comparisons focus on a single post-processed quantity and leave open the possibility of fortuitous error cancellations that obscure the true degree of reproducibility.Published comparisons of time-resolved TDDFT data 16,17 remain rare but reveal significant discrepancies in the core contribution to stopping power that were attributed to details of the pseudopotentials. 17 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 3 and temperature of 2 eV.We further focus on a single alpha-particle velocity of 5 at.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 allelectron calculations with bare Coulomb potentials.We compare both time-resolved data and average stopping powers computed using four different TDDFT codes: a custom extension 18,19 of the Vienna Ab initio Simulation Package (VASP), [20][21][22] Stochastic and Hybrid Represen-tation Electronic structure by Density functional theory (SHRED), [23][24][25] Qball, 26,27 and Inq. 16e 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 highenergy-density systems. 7

II. METHODS
The central numerical task of a real-time TDDFT calculation 28 consists of solving the time-dependent Kohn-Sham (TDKS) equations, subject to some initial condition ϕ j (r, t = 0) and closure where ϕ j are Kohn-Sham (KS) electronic orbitals with temperature-dependent occupations f j (T e ) and n(r, t) is the electron density.The Kohn-Sham Hamiltonian H is a functional of the electron density and a sum of kinetic energy, external potential, Hartree, and exchangecorrelation (XC) terms: where V ext includes electron-ion interactions, V Har captures the classical Coulombic electron-electron interaction, and V XC contains quantum-mechanical corrections to the electron-electron interaction.Within a TDDFT stopping power calculation, electrons respond to a moving particle that introduces explicit time-dependence within V ext .
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 with a cutoff energy of at least 2000 eV, began from Mermin-Kohn-Sham equilibrium states, 29 and held ionic velocities fixed.Other aspects of the simulations varied as described in the following subsections and summarized in Table I.
We obtain stopping powers from stopping forces F (x) = −F(x) • v against the projectiles'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.
Again, different choices exist for estimating the average stopping power from time-dependent forces.One method first computes the stopping work W (i.e., energy deposited by the projectile) as a function of x and then fits W (x) to a linear model, where S fit and W 0 are model parameters and [x 0 , x f ] is the range of projectile displacements included in the analysis.Another approach averages the stopping force directly, effectively computing 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 offaxis alpha particle trajectory optimized according to the methods developed by Kononov et al. 17 (see Fig. 1a) and (ii) an elongated simulation cell and the alpha particle travelling parallel to the long axis, following the approach of earlier work by Nichols et al. 5 (see Fig. 1b).For both setups, the atomic configuration was generated from separate ab initio molecular dynamics (MD) simulations as described in Appendix A. Our comparisons often focus on 512-atom hydrogen configurations, though we also consider other simulation cell sizes where indicated.

A. VASP calculations
Although typical VASP [20][21][22] calculations use the projector augmented-wave (PAW) method, 30 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.  Sec.III, we consider several variations upon the choices listed here as described in the text.For the initial condition, α ∈ IC (α ̸ ∈ IC) indicates that the alpha particle projectile was included in (excluded from) from the Mermin-DFT calculation for the equilibrium electronic orbitals.In some cases, the time step and total simulation time varied inversely with the projectile velocity v. Ftot denotes the total force on the projectile and ∆Fei denotes the electron-ion contribution excluding the force on the projectile under a fixed initial electron density.S fit and Savg estimate average stopping power through a linear fit and direct averaging, respectively (see Eqs. ( 4) -( 6)).The adiabatic local density approximation (LDA) 31,32 was used in VASP calculations throughout this work, and we compare VASP results for LDA and Perdew-Burke-Ernzerhof (PBE) 33 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 at.u. case scrutinized throughout this work, plane-wave cutoffs of 2000 eV 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 18,19 of VASP uses the Crank-Nicolson (CN) time-stepping algorithm and solves the associated update equations using the Generalized Minimal Residual (GMRES) method 34 with a GMRES relative tolerance of 10 −9 .For fast projectiles with v > 3.5 at.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 at.u. (1.5 at.u.) of velocity.

B. SHRED calculations
In addition to the deterministic Kohn-Sham formulation used by the other codes, SHRED [23][24][25] 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 finitetemperature 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. 24.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 DFT with 640 KS orbitals at low (1.5 at.u.) and high (5 at.u.) velocities show a difference of less than 1% percent.
The SHRED calculations used the PBE exchangecorrelation functional 33 and either bare Coulomb potentials or Hartwigsen-Goedeker-Hutter (HGH) pseudopotentials, 35,36 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. 37he 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 itwould also average to 0 due to the same symmetry.This force, denoted as F (0) ei , 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 (0) ei .

C. Qball calculations
The Qball 26,27 calculations used Hamann-Schluter-Chiang-Vanderbilt (HSCV) norm-conserving pseudopotentials 38 with outermost cutoff radii of 0.8 at.u. (0.5 at.u.) for the alpha particle (hydrogen ions).The electronic states evolved over time according to the enforced time-reversal symmetry (ETRS) algorithm 39 with fourth-order Taylor expansions for the exponentials. 27Rather 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. 27xchange and correlation was treated at the LDA 31,32 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 at.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/ Å.

D. Inq calculations
The Inq 16 simulations follow very similar methodology as described for Qball in Sec.II C.However, we used the default pseudopotential set within Inq: optimized normconserving Vanderbilt (ONCV) pseudopotentials 40 with outermost cutoff radii of 1.25 at.u. (1.0 at.u.) for the alpha particle (hydrogen ions).Since these pseudopotentials are even softer than those used in the Qball calcula- Results are labelled by the TDDFT code (with corresponding methods summarized in Table I) and the number of hydrogen atoms in the simulation.mSHRED and dSHRED denote mixed stochastic-deterministic and deterministic (KS) TDDFT variants within SHRED.In (a), the VASP (SHRED) calculations used cubic (elongated) simulations cells with offaxis (on-axis) trajectories.In (b), all simulations used the same cubic / off-axis setup.
tions, we expect excellent convergence at the same cutoff energy.Also, the Inq ETRS implementation performs up to 5 self-consistent iterations to obtain n(r, t + ∆t) and H(t + ∆t), further enhancing convergence with respect to time step.

III. RESULTS AND DISCUSSION
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. 2b), but modest deviations arise when using pseudopotentials and/or changing the simulation cell shape and alpha particle trajectory.
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 at.u. in Fig. 3, where all calculations used the same Total stopping forces on an alpha particle moving through warm dense hydrogen (1 g/cm 3 , 2 eV) with 5 atomic units (at.u.) of velocity as computed by each TDDFT code for the 512-atom cubic test system.(a) shows the force against the direction of projectile motion and (b) shows differences from the VASP data (solid blue).off-axis alpha particle trajectory through the same 512atom hydrogen configuration within the cubic simulation setup (see Fig. 1a), i.e., the same setup used for Fig. 2b.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 the following subsections, 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.

A. Influence of model approximations
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. 17n 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. 41Recent 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. 17Although 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 entire simulation length 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, 17 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 Electron-ion stopping forces on an alpha particle moving through warm dense hydrogen (1 g/cm 3 , 2 eV) with 5 at.u. of velocity for the cubic test system.We compare results from a series of VASP calculations using bare Coulomb or PAW potentials, LDA or PBE, and the alpha particle included in the initial condition (α ∈ IC) or excluded from it (α / ∈ IC).(a) shows the force against the direction of projectile motion and (b) shows differences from the case with bare Coulomb potentials, LDA, and the alpha particle included in the intial condition (solid blue).
and HGH pseudopotentials with cutoff radii of 1 at.u. for both species.

B. Influence of initial condition
In addition to the form of the Hamiltonian, the solution to the TDKS equations (Eqs.( 1) -( 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 into some stopping power models based on the average-atom method, 4 and close agreement for the initial DOS helps ensure agreement among stopping powers.
All the TDDFT simulations began from Mermin-KS equilibrium states, 29 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 DsOS 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.(a) Densities of states (DsOS) for hydrogen at 1 g/cm 3 and 2 eV computed using each code with the pseudopotentials and XC functionals listed in Table I.To obtain smooth curves, these Mermin-DFT calculations used a denser 2 3 Γ-centered reciprocal-space quadrature and the DsOS were obtained by applying a Gaussian broadening of 2 eV to the Mermin-KS eigenvalues.The "α ∈ IC" curve (dotted red) was computed with VASP and included the alpha particle projectile.The DFT data departs modestly from the ideal Fermi gas (dotted black) at the corresponding electron density.Gray shading indicates the occupied DOS.In (b), differences are taken relative to the VASP DOS (solid blue).
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 42 and pre-equilibrium behavior in 2D materials, 43 but light-ion stopping powers in bulk materials do not depend strongly on the initial condition 41 because the projectile's charge state quickly equilibrates after traversing a few Ås. 44ere, 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.

C. Influence of observable choice
Even with the same underlying electron dynamics, it is possible to extract an average stopping power from different observables.For example, Ref. 41 previously demonstrated the equivalence of analyzing the total energy under fixed ion velocities and analyzing the total force on the projectile; these observables are related through the work done by maintaining a constant projectile velocity.Here, we consider the role of different contributions to the total force on the projectile, where F ei is the electron-ion force and F ii is the ion-ion force.We also consider where ei represents the force that a fixed initial electronic system would induce, a contribution often excluded from SHRED calculations as described in Sec.II B.
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 ionion contribution F ii dominates the absolute energy scale of the total force (cf.Figs.6a and b), its time average quickly decays toward 0 and amounts to 2.6% of the converged total stopping power (see Fig. 6d).Similarly, including the force due to the unperturbed electron density F (0) ei changes the extracted stopping power by only 0.6% within VASP and 0.2% within SHRED.
Contributions from F ii and F (0) ei 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 at.u., including F ii and F (0) ei 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.

D. Influence of averaging procedure
Given a choice of observable, different analysis methods are possible for estimating an average stopping power (see Eqs. ( 5) and ( 6)).Ref. 17 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 .Different contributions to the stopping forces on an alpha particle moving through warm dense hydrogen (1 g/cm 3 , 2 eV) with 5 at.u. of velocity as computed by VASP and SHRED for the 512-atom cubic test system.Panels (a) -(c) compare the total force Ftot, ion-ion force Fii, electron-ion force Fei, and the difference between electron-ion forces under time-evolving and static electron densities ∆Fei = Fei − F In Fig. 7a, the stopping work computed from F tot contains sharp features arising from strong occasional ionion force pulses (see Fig. 6a).Consistent with Ref. 17, 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's informed by a range of W (x) data (see Fig. 7b).Excluding ionion 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  4)) as computed from different forces within VASP on an alpha particle moving through warm dense hydrogen (1 g/cm 3 , 2 eV) with 5 at.u. of velocity.(b) Average stopping power estimated as S fit (solid, see Eq. ( 5)) and Savg (dotted, see Eq. ( 6)).
contributing sharp features to F ei and/or higher projectile velocities where F (0) ei ≫ ∆F ei .

E. Influence of simulation setup
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 17,41,45 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, 17,41,45 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. 17Using an ensemble of several trajectories (either on-axis or off-axis) can further improve sampling fidelity. 46sing the quantitative metric proposed by Ref. 17 to assess trajectory quality, the ten random on-axis trajectories deviate from ideal sampling by D H between 0.1 and 0.3 when considered individually.When taken in aggregate, the ensemble of on-axis trajectories achieves a much lower D H value of 0.05 -0.06, comparable to the sampling fidelity of the optimized off-axis trajectories in the cubic setup which achieve D H values around 0.04.Interestingly, the distribution of D H values in the fully thermalized 1024-atom supercell (turquoise squares in Fig. 8) appears similar to the replicated supercells (orange squares) despite the less periodic environment sampled in the former case.For the largest supercells considered in this work, stopping powers predicted from the optimized offaxis trajectory and an ensemble average over on-axis trajectories agree within 4%, confirming the validity of both approaches.
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. 47The 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. 17These so-called "ouroboros effects" would be more severe for a long onaxis 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 onaxis 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 512atom 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 approaches.Depending on the approach, both effects can significantly influence the accuracy of stopping power predictions.Long onaxis trajectories suffer from both 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.Future work may investigate the ability of a compromising approach relying on intermediatelength off-axis trajectories to better balance computational costs in reducing both types of errors.

IV. CONCLUSIONS
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 48 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.
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 exchangecorrelation functional, 33 and the Hartwigsen-Goedeker-Hutter (HGH) pseudopotential 35,36 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 50 and the nonlinear conjugate gradient algorithm outlined by Jiang and Yang 51 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

FIG. 1 .
FIG. 1.(a) Cubic and (b) elongated simulation setups considered in this work.Both contain one alpha particle projectile (red) moving through a 512-atom hydrogen (white) configuration with trajectories indicated by red arrows.The hydrogen geometries were obtained from MD simulations as described in Appendix A. In (a), the alpha particle's initial position and direction of motion were optimized according to the methods proposed by Ref. 17.In (b), the alpha particle begins at an arbitrary point on the left edge of the supercell and travels in the long direction.

FIG. 2 .
FIG.2.Average stopping powers in warm dense hydrogen (1 g/cm 3 , 2 eV) as a function of the alpha particle velocity.Results are labelled by the TDDFT code (with corresponding methods summarized in TableI) and the number of hydrogen atoms in the simulation.mSHRED and dSHRED denote mixed stochastic-deterministic and deterministic (KS) TDDFT variants within SHRED.In (a), the VASP (SHRED) calculations used cubic (elongated) simulations cells with offaxis (on-axis) trajectories.In (b), all simulations used the same cubic / off-axis setup.
FIG. 3.Total stopping forces on an alpha particle moving through warm dense hydrogen (1 g/cm 3 , 2 eV) with 5 atomic units (at.u.) of velocity as computed by each TDDFT code for the 512-atom cubic test system.(a) shows the force against the direction of projectile motion and (b) shows differences from the VASP data (solid blue).
FIG.4.Electron-ion stopping forces on an alpha particle moving through warm dense hydrogen (1 g/cm 3 , 2 eV) with 5 at.u. of velocity for the cubic test system.We compare results from a series of VASP calculations using bare Coulomb or PAW potentials, LDA or PBE, and the alpha particle included in the initial condition (α ∈ IC) or excluded from it (α / ∈ IC).(a) shows the force against the direction of projectile motion and (b) shows differences from the case with bare Coulomb potentials, LDA, and the alpha particle included in the intial condition (solid blue).
FIG. 5.(a) Densities of states (DsOS) for hydrogen at 1 g/cm 3 and 2 eV computed using each code with the pseudopotentials and XC functionals listed in TableI.To obtain smooth curves, these Mermin-DFT calculations used a denser 2 3 Γ-centered reciprocal-space quadrature and the DsOS were obtained by applying a Gaussian broadening of 2 eV to the Mermin-KS eigenvalues.The "α ∈ IC" curve (dotted red) was computed with VASP and included the alpha particle projectile.The DFT data departs modestly from the ideal Fermi gas (dotted black) at the corresponding electron density.Gray shading indicates the occupied DOS.In (b), differences are taken relative to the VASP DOS (solid blue).
FIG. 6.Different contributions to the stopping forces on an alpha particle moving through warm dense hydrogen (1 g/cm 3 , 2 eV) with 5 at.u. of velocity as computed by VASP and SHRED for the 512-atom cubic test system.Panels (a) -(c) compare the total force Ftot, ion-ion force Fii, electron-ion force Fei, and the difference between electron-ion forces under time-evolving and static electron densities ∆Fei = Fei − F ) shows the average stopping power generated by each force contribution with colors corresponding to panels (a) -(c).

TABLE I .
Differences among TDDFT simulations performed with each code to produce the data in Figs. 2. Throughout DH quantifies deviations from ideal environment sampling, and DO is the minimum distance between the projectile and periodic images of its previously traversed path.Symbol shapes indicate supercell size, while colors and fill styles correspond to different schemes for obtaining an average stopping power.Empty symbols indicate random on-axis trajectories with a single pass through an elongated simulation cell such that the projectile stops 2 Å before reaching its starting point.Data points with a solid orange outline used simulation cells comprised of a 256-atom cell replicated in the direction of projectile motion, while dashed turquoise outlines used a fully equilibrated 1024-atom cell.Filled orange and turquoise points indicate aggregate results for corresponding trajectory ensembles.Semi-filled points use a single, 80 Å-long trajectory in the on-axis, elongated setup (red, right-filled) or the optimized off-axis, cubic setup (green, top-filled).