==== Front
Philos Trans A Math Phys Eng Sci
Philos Trans A Math Phys Eng Sci
RSTA
roypta
Philosophical transactions. Series A, Mathematical, physical, and engineering sciences
1364-503X
1471-2962
The Royal Society

10.1098/rsta.2022.0325
rsta20220325
100915215989Articles
Research Articles
Development of a new quantum trajectory molecular dynamics framework
Development of a new quantum trajectory molecular dynamics framework
http://orcid.org/0000-0002-7013-208X
Svensson Pontus Conceptualization Data curation Formal analysis Investigation Methodology Project administration Resources Software Validation Visualization Writing – original draft Writing – review & editing pontus.svensson@physics.ox.ac.uk
1
http://orcid.org/0000-0002-0104-1424
Campbell Thomas Conceptualization Investigation Writing – review & editing 1
http://orcid.org/0000-0002-7204-0509
Graziani Frank Conceptualization Investigation Writing – review & editing 2
http://orcid.org/0000-0002-9725-9208
Moldabekov Zhandos Data curation Methodology Writing – review & editing 3 4
Lyu Ningyi Data curation 5
http://orcid.org/0000-0002-3262-1237
Batista Victor S. 5 6
Richardson Scott Supervision 7
http://orcid.org/0000-0003-1016-0975
Vinko Sam M. Conceptualization Funding acquisition Investigation Resources Supervision Writing – review & editing 1 8
http://orcid.org/0000-0002-4153-0628
Gregori Gianluca Conceptualization Funding acquisition Investigation Resources Supervision Writing – review & editing 1
1 Department of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, UK
2 Lawrence Livermore National Laboratory, Livermore, CA 94550, USA
3 Center of Advanced Systems Understanding (CASUS), D-02826 Görlitz, Germany
4 Helmholtz-Zentrum Dresden-Rossendorf (HZDR), D-01328 Dresden, Germany
5 Department of Chemistry, Yale University, New Haven, CT 06520, USA
6 Yale Quantum Institute, Yale University, New Haven, CT 06511, USA
7 AWE, Aldermaston, Reading, Berkshire RG7 4PR, UK
8 Central Laser Facility, STFC Rutherford Appleton Laboratory, Didcot OX11 0QX, UK
One contribution of 11 to a theme issue ‘Dynamic and transient processes in warm dense matter’.

24 7 2023 July 24, 2023
24 7 2023 July 24, 2023
381 2253 2022032514 11 2022 November 14, 2022
19 1 2023 January 19, 2023
© 2023 The Authors.
2023
https://creativecommons.org/licenses/by/4.0/ Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

An extension to the wave packet description of quantum plasmas is presented, where the wave packet can be elongated in arbitrary directions. A generalized Ewald summation is constructed for the wave packet models accounting for long-range Coulomb interactions and fermionic effects are approximated by purpose-built Pauli potentials, self-consistent with the wave packets used. We demonstrate its numerical implementation with good parallel support and close to linear scaling in particle number, used for comparisons with the more common wave packet employing isotropic states. Ground state and thermal properties are compared between the models with differences occurring primarily in the electronic subsystem. Especially, the electrical conductivity of dense hydrogen is investigated where a 15% increase in DC conductivity can be seen in our wave packet model compared with other models.

This article is part of the theme issue ‘Dynamic and transient processes in warm dense matter’.

warm dense matter
, wave packet molecular dynamics
, non-adiabatic electron dynamics
Atomic Weapons Establishment http://dx.doi.org/10.13039/501100000766 Oxford Physics Endowment for Graduates Center for Advanced Systems Understanding cover-dateJuly 24, 2023
==== Body
pmc1. Introduction

The establishment of high power lasers facilities during the last decades has been instrumental in the achievements towards inertial confinement fusion (ICF) [1–3], but also for the creation of high-density and high-temperature conditions [4] otherwise only found in astrophysical objects [5,6]. Furthermore, X-ray lasers are now able to reach complementary high-pressure regions in phase space [7,8]. One of the exotic states now accessible is warm dense matter (WDM), which exists in gas giants [9–13], brown [14] and white dwarf stars [15,16], the crust of neutron stars [17,18] and during the compression of an ICF capsule [19]. WDM is a strongly coupled quantum plasma, with ions moving in a partially degenerate electron fluid with kinetic energy comparable with the ion–ion interaction energy [20]. Consequently, WDM inherits properties from both condensed matter systems and classical plasmas, a challenging combination to model. Various computational techniques are commonly used to describe these systems, yielding similar thermodynamic [21] and acoustic properties [22,23], although dynamic properties differ by orders of magnitude [24]. These uncertainties limit our understanding of, for example, the Jovian interior [25], or the modelling of ICF implosions [26].

The three main complications in modelling WDM are electron degeneracy, strong ion correlations and the separation in time scales between the electron and ion dynamics. A full solution would require a quantum mechanical treatment of the electrons, resolving electron dynamics while considering phenomena on the ion time scale. Consequently, explicit models of ionic motion span a wide range of theories, including classical systems with effective ion–ion interactions [27–29], classical electrons with effective quantum statistical potentials (QSP) [30–32], Bohmian mechanics [33], density functional theory molecular dynamics (DFT-MD) using both orbital-free [34,35] and Kohn–Sham [12,36,37] DFT-variants, phenomenological quantum hydrodynamics based on DFT-functionals [38–40], time-dependent DFT [41–43] and quantum Monte-Carlo and path integral Monte-Carlo [13,44,45] approaches. Coarse-grained models with effective interactions are fundamentally based on reconstructing some equilibrium property, the choice of which is arbitrary and limited to a specific thermodynamic condition, whereas experimental realizations are commonly non-stationary [7,46–50]. Furthermore, models rooted in the Born–Oppenheimer approximation—where the electrons are treated adiabatically—e.g. DFT-MD, cannot capture a dynamic electron response, believed to be important for the description of dynamic properties such as some transport coefficients [51], stopping power [39] and energy transfer between the electronic and ionic subsystems. However, time-dependent approaches are computationally costly, and are typically limited in terms of particle numbers and time scales of studied phenomena.

Wave packet molecular dynamics (WPMD) [52,53] is a family of models in which the electron dynamics are computed explicitly, while simulating hundreds to thousands of particles over ionic time scales. This is made possible by restricting the wave function of each electron to a parameterized functional form. We present an extension to existing wave packet formulations—applicable to the WDM regime—in which the wave packets can be elongated in arbitrary directions. The model accounts for the long-range behaviour of electrostatic interactions and of fermionic properties by effective Pauli interactions, while implemented within the scalable molecular dynamics framework LAMMPS [54] to treat systems with thousands of particles.

In the following section, the theoretical model is described, after which §3 outlines the numerical details and performance of the implementation. The model is compared with other computational techniques in §4, where we apply it to ground state and dynamic properties of a dense hydrogen plasma. We compute some structural and transport properties, which are compared with an isotropic wave packet model. We conclude with a summary of our results.

2. Theoretical description

Originally proposed in the 1970s as an approximate solution to Schrödinger’s equation [55,56], wave packet models can systematically be derived from variations of the action, 2.1 S=∫ dt ⟨Q|iℏddt−H^|Q⟩,

where H^ is the system Hamiltonian and the state, |Q⟩=|Q(Qμ)⟩, is restricted to some manifold, M, defined by the adopted wave packets and parameterized by its parameters Qμ. The resulting time evolution reproduces the true quantum dynamics to the best of its ability being restricted to the manifold, M, during short time scales of length δt. Concretely, it can be shown ⟨Δ(t,δt)|Δ(t,δt)⟩ is minimized to O(δt3), where |Δ(t,δt)⟩=|Ψ(t+δt)⟩−|Q(t+δt)⟩ and |Ψ(t+δt)⟩ is the true solution to Schrödinger’s equation starting from |Ψ(t)⟩=|Q(t)⟩ [52]. The long-time evolution is constrained by appropriate conservation laws, most notably energy conservation [57].

In general, the equations of motion are quasi-Hamiltonian 2.2 iℏ∑νCμνdQνdt=∂H∂Qμ∗,

where 2.3 H=⟨Q|H^|Q⟩⟨Q|Q⟩≡⟨H^⟩andCμν=∂2∂Qμ∗∂Qνln⁡(⟨Q(Qμ∗)|Q(Qν)⟩).

Fermions are described by states antisymmetric under exchange and Slater-determinants have been considered in [58–60]. However, this approach scales unfavourably with particle number N. Instead, here we employ a product state 2.4 |Q⟩=|q1⟩⊗|q2⟩⊗⋯⊗|qN⟩,

of single-particle orbitals, |qi⟩, and ⊗ is the tensor product. Exchange effects are approximated by Pauli potentials of the type first introduced by Klakow et al. [61,62]. This structure simplifies Cμν, which becomes block diagonal, and the orbitals |qi⟩ only couple through the energy H.

(a) Wave packets

The choice of wave packet shape is central to the model, dictating the states that can be described [63]. Most commonly, isotropic Gaussians are used—primarily motivated by computational ease—yet other variants exist, see Grabowski [53] and references therein. To account for local gradients, an anisotropic wave packets form is introduced 2.5 ⟨x|qi⟩=((2π)3det(Σi))−1/4×exp⁡[−ξi⊺(14Σi−1−iℏΠi)ξi+iℏpi⊺ξi],

where ξi=x−ri. The wave packet is parameterized by 18 degrees of freedom, the position ri, momentum pi and two symmetrical 3×3 matrices, Σi, describing the elongation and orientation and Πi the associated momentum to Σi. A similar type of wave packet has been treated previously, describing molecular binding in water molecules [64], and is generally believed to improve the description of molecular states [53].

The equations of motion for the functional form (2.5) have a classical-looking structure, for the ‘classical’ degrees of freedom [64] 2.6a dridt=∂H∂pianddpidt=−∂H∂ri,

and the ‘internal’ dynamics of the wave packet follow as 2.6b ddtΣiαβ=ταβ∂H∂ΠiαβandddtΠiαβ=−ταβ∂H∂Σiαβ,

where Σiαβ=(Σi)αβ and Πiαβ=(Πi)αβ are the components of the symmetric matrices. The pre-factor ταβ is unity if α=β and one half otherwise, which accounts for the symmetric structure of Σi and Πi where Σiαβ (Πiαβ) and Σiβα (Πiβα) are treated as symbolically the same. Specifically, we consider a charged system of classical ions, with position RI, momentum PI, charge ZIe and mass MI, and quantum electrons with position x^i and momentum p^i operators as well as charge −e and mass m. The system is described by the Hamiltonian 2.7 H^=∑IPI22MI+∑I<JZIZJe2|RI−RJ|+∑ip^i22m+∑i<je2|x^i−x^j|−∑i∑IZI e2|x^i−RI|,

the state average of which is required for the time evolution. The average kinetic energy 2.8 ⟨p^i22m⟩=pi22m+2mTr{ΠiΣiΠi}+ℏ28mTr{Σi−1}

includes both a classical contribution and a part internal to the wave packet. The last term in equation (2.8) is the so-called shape-kinetic energy [65], which keeps Σi positive definite and the wave packet well defined during the time evolution. The interaction terms have not been evaluated explicitly and the following section is dedicated to the treatment of these terms.

(b) Generalized Ewald summation

Within molecular dynamics, it is desirable to truncate pair-interactions at some distance such that the computation formally scales as O(N) [66]. However, in our case the electrostatic interaction is long-range [67] and it is beneficial to perform the split [68,69] 2.9 1r=erfc(gr)r+erf(gr)r,

chosen so that the first term can be truncated at a distance of order g−1, while the second term is regular as r→0 and efficiently evaluated in Fourier space. The Ewald parameter g is chosen to optimize performance. Below we present a self-consistent treatment of both terms, as the long-range part has only been mentioned once for isotropic wave packets [70] and is commonly neglected.

(i) Short-range forces

In the case of a Gaussian interaction kernel, the required state average can promptly be evaluated [71,72]. Therefore, we construct a Gaussian decomposition of the interaction kernel, 2.10 erfc(gr)r≃∑pcp e−αpr2,

where the coefficients cp and αp are fitting parameters. A robust numerical scheme to perform the decomposition is described in appendix A, where typically only 5–15 modes are required. By approximating the potential form, the notion of energy conservation is retained. We note that this is not the case for methods based on Taylor expansions [64] or on direct numerical evaluations [73], due to either truncation errors or numerical noise.

(ii) Long-range forces

To limit surface effects in a finite size simulation, the simulation box is periodically repeated. Periodic images are included in accordance with the standard treatment of Ewald summation with particles positioned at ri+Ln for all n∈Z3 and L being the length of a cubic simulation cell. The interaction energy is 2.11 ∑i<j⟨V(xij)⟩→12∑i,j′∑n∈Z3⟨V(xij-Ln)⟩=12∑i,j∑n∈Z3⟨V(xij-Ln)⟩⏟Ek−12∑i=j⟨V(xij)⟩⏟Es,

where V is the long-range part of the Coulomb interaction in equation (2.9). The special case i=j is excluded when n=0 (denoted by the primed sum) resulting in two distinct terms, the main contribution Ek and the self-energy Es. In appendix B, we evaluate Ek in reciprocal space to be 2.12 Ek=12L3∑k≠04πk2 e−k2/(4g2)|ρ~uc(k)|2,

where ρ~uc(k) is the charge density 2.13 ρ~uc(k)=e∑iZi e−ik⋅ri e−12k⊺Σik,

and k=2πn/L. Equation (2.12) converges rapidly due to the exponential factor. The self-energy term, Es, is independent of ri and does not influence the classical degrees of freedom directly. Although for distributed particles, it does depend on Σi and influences the internal dynamics. Based on the decomposition of the long-range interaction kernel (appendix A) the self-energy is evaluated in appendix B.

(c) Pauli interactions

The difference in the kinetic energy between a pair-wise antisymmetrized state and the product state has often been used as a Pauli potential to correct for the fermionic structure of electrons [61,62,70]. The electron force field (eFF) model introduced fitting parameters in the Pauli interaction to achieve stable bounds for elements with Z≤6 [74] and has been widely used with minor modifications [22,23,75–80]. Angermeier & White [81] considered exchange contributions to the interaction terms while including a correlation potential with a free parameter. This treatment of the Pauli interaction is extended here to anisotropic Gaussian states.

We construct the potential by considering two electrons, i and j, with the Hamiltonian 2.14a H^2=p^i2+p^j22m+e2|x^i−x^j|+Vbg(x^i,t)+Vbg(x^j,t),

with a background interaction from all other particles in the system 2.14b Vbg(x,t)=∑k≠i,j e2∫ d3xk|⟨xk|qk⟩|2|x−xk|−∑IZI e2|x−RI|,

where k (I) runs over all other electrons (ions) in the system. The two-electron system can be characterized based on its spin structure either as a singlet or a triplet state, requiring either a symmetric or antisymmetric spatial state 2.15 Singlet: |qi⟩⊗|qj⟩+|qj⟩⊗|qi⟩Triplet: |qi⟩⊗|qj⟩−|qj⟩⊗|qi⟩,

written here in terms of single-particle orbitals |qi⟩.

For equal spin particles only the spatially antisymmetric state is allowed and the Pauli potential VijP is the difference between ⟨H^2⟩ for the triplet and the product state |qi⟩⊗|qj⟩. In the case of opposite spin particles, the spatial state depends on the spin structure, but along the lines of Angermeier & White [81] a correlation potential is introduced based on the singlet state, ρVijC, multiplied by the parameter ρ. Therefore, the potentials are 2.16 VijP/C=∓⟨qiqj|H^2|qjqi⟩−⟨qiqj|H^2|qiqj⟩|⟨qi|qj⟩|21∓|⟨qi|qj⟩|2,

which can be shown to scale as |⟨qi|qj⟩|2 for large particle separations. Due to Gaussians being localized states, this gives a short-range interaction between particles i and j. The background term Vbg introduces a long-range dependence in terms of the third particle and is accounted for by a type of Ewald summation, see appendix C. The state averages of the interaction terms in equation (2.16) are evaluated based on the Gaussian mode decomposition.

The correlation potential, the spin interaction between opposite spin electrons, is constructed based on the same premise as the Pauli potential in a pair-wise approximation, however, with an additional parameter ρ which needs to be chosen a priori. In the case of a ground state helium or molecular hydrogen, the electronic structure is well described by a single state and ρ=1 is an appropriate choice. For a free electron gas, this would overestimate the correlation effects as the appropriate two-particle state is not simply the singlet state, and therefore ρ<1 is more suitable.

Lastly, the above-presented schemes—although widely used—are in general limited to the weakly and moderately degenerate systems, cause considering the example of Pauli blocking. The Pauli potential in equation (2.16) appears to be divergent as the orbital overlap tends to unity; however, the numerator vanishes as well when |qj⟩→|qi⟩ resulting in only a finite energy barrier. The remaining part of the Pauli exclusion should be accounted for by the left-hand side of equation (2.2) by a complete antisymmetrization scheme.

(d) Confining potentials

It has been well documented that at sufficiently high temperatures wave packets tend to expand indefinitely [22,58,70,82–84] and the wave packet may extend over multiple ions without the ability to localize on multiple sites [85]. If a wave packet is spread too large, it effectively ceases to interact with other particles as the charge density effectively becomes flat. Multiple approaches to counter this expansion have been proposed, see [79,83,86]. Currently, we employ an additional potential energy term of the form 2.17 VΣ=12Aw∑i=1N∑α=13(σi,α−lw)2θ (σi,α2−lw2),

where σi,α2 is the αth eigenvalue of Σi and θ (x) is the Heaviside step function. The parameters lw and Aw set the width of a free particle σfree(lw,Aw)≳lw by balancing the shape-kinetic energy. This potential is rotationally invariant and acts only on wave packets with a width larger than lw. Furthermore, the confinement reduces to the commonly used potential based on a harmonic potential centred at the particle position in the limit of lw→0. In this specific limit, the potential has also been used to address the heat capacity in the classical limit [87].

3. Numerical realization

The standard velocity-Verlet integrator almost exclusively used for MD simulations is not appropriate for our model because of the momentum-dependent Pauli potentials [52] resulting in a non-separable Hamiltonian. This prevents a straightforward generalization of the velocity-Verlet algorithm which is based on the ability to separate the Hamiltonian into terms where the dynamics following from each term in isolation can be solved exactly [88]. Explicit Runge–Kutta methods of orders 2 and 4 are employed instead. Furthermore, this momentum dependence of the potential affects the interpretation of temperature in the system, further described in appendix D.

The time integrator, the generalized Ewald summation and the Pauli interaction, are all natively implemented in LAMMPS [54], a MD framework written in C++ which uses MPI to distribute the computation [89]. Figure 1 b shows the computational time for a varying degree of parallelization for a fixed system of 2000 protons and an equal number of electrons. In particular, a good scaling of the pair interaction and the Ewald summation is established as the computation is distributed. The synchronization time, the time different processes need to wait for each other due to an unbalanced load caused by statistical fluctuations in the number of particles in the region assigned to each processor, limits the efficiency of the parallelization when only a few particles are assigned to each process. In the future, dynamic load balancing could potentially address this issue; however, it should be noted that the point of the plateau moves further out as the size of the system is increased. Figure 1. Scaling of computational cost for a quasi-neutral system with N protons, density corresponding to rs=2 and temperature equal to the Fermi temperature. (a) Scaling with particle number for different levels of parallelization, in the range between N1.1 and N1.3. The Ewald parameter g was optimized to within ±0.5aB−1 for each case. (b) MPI parallelization of the computation for a system of N=2000 ions with a fixed Ewald parameter. Showing the full computation (circles), pair interaction (squares), Ewald summation (diamonds), communication between MPI processes (upwards triangles) and the synchronization time different processors need to wait for each other due to an unbalanced load (downwards triangles), which is compared with the optimal scaling based on the single thread performance (dashed).

Furthermore, figure 1 a also demonstrates the scaling of computational cost with particle number for a test system where rs=(4π/3naB3)−1/3=2 where n is the proton number density and aB is the Bohr radius. Close to linear scaling in particle number is demonstrated, showing the feasibility of employing this modelling technique for large systems of particles. The exact exponent varies in the range 1.1–1.3, depending on the degree of parallelization, caused by different limiting factors in the computation. The synchronization time is most likely one of these factors for the high parallelization case.

4. Test systems

(a) Ground state properties

Ground state properties of the wave packet models can be obtained by the introduction of a generalized friction term into the equations of motion (2.6). Some care is needed to guarantee continuous energy loss due to the momentum dependence of the Pauli potential. This is further described in appendix E.

The ground state of isolated atoms is spherically symmetric and does not use the additional degrees of freedom of the elongated wave packets. One of the simplest physical systems which naturally breaks this symmetry are diatomic molecules and in particular diatomic hydrogen (H2), where the ability for the wave packet to stretch is believed to be crucial for molecular binding [53]. The ground state energy of H2 within the wave packet model for both elongated and isotropic Gaussians is shown in figure 2 for a varying nuclear separation δ. The elongated wave packets demonstrate an improvement over the isotropic ones when δ<2.8aB, above which the electron density is localized on each ion and close to spherically symmetric. Furthermore, the isotropic model transitions to an electron density localized on each nucleus at a significantly shorter nuclear separation, δ≈1.8aB, compared with our model at δ≈2.8aB. Figure 2. Ground state energy of H2-molecule, E, as a function of nuclear separation δ, computed for the full wave packet model (elongated), the isotropic one (isotropic) and the virtually exact result from Angermeier & White [81] (exact). The ground state is modelled as a triplet state, with ρ=1 in the correlation potential (2.16). The equilibrium displacement is marked for the two wave packet models (vertical lines) and the energy of two isolated hydrogen atoms for the wave packets as well as the exact result are shown with horizontal markers.

The energy difference between our model and the exact result is close to constant for nuclear separation larger than the equilibrium position. In this respect, the presented wave packet model is as good a descriptor of the molecular ground state, H2, as the atomic one, 2×H, suggesting that for further improvements one needs a more versatile description of the isotropic atomistic limit first. The present model has the ability to treat each direction separately and what is limiting the agreement is the restriction on functional form.

(b) Binary collisions

The dynamical properties of the wave packet model have been tested for electron–ion scattering against a full numerical solution of the Schrödinger equation realized by the SOFT code [90–92]. The resulting trajectory data for impact parameters in the range of 0.6–3.0aB is shown in figure 3, alongside the time evolution of the extent of the wave function, which is compared with the result from both isotropic and elongated wave packets. The centre of mass trajectories agree well between all three different sets of simulations over long-length scales. In the full numerical solution, the electron density can split and partially bind to the ion core, a qualitative feature the wave packets cannot reproduce due to their limited functional form. However, this is of minimal importance for the trajectories discussed here. At these energies, only a minor fraction of the electron wave packet gets bound so that the centre of mass agrees well with the mode position in the full numerical simulation. The two wave packet models differ in the internal degrees of freedom, where the isotropic wave packet does not have the flexibility of the complete model. The isotropic wave packet model cannot reproduce the internal dynamics of the SOFT computation as well as the elongated wave packet model for the larger set of impact parameters. In the case of smaller impact parameters, the full numerical solution has wave packets with non-Gaussian structure during the close approach between the electron and the ion, impacting the subsequent evolution of the wave function. One of the more important characteristics for determining the dynamic properties of the system is the energy transfer in particle collisions, which is shown for the electron–proton scattering in figure 4. The energy transfer is calculated by treating the proton dynamically in contrast to the case shown in figure 3. Differences between the wave packet models can in particular be seen for intermediate impact parameters, between the classical behaviour for large impact parameters and the symmetrical configuration when b=0. Furthermore, the difference is the most pronounced for smaller wave packets where the asymmetry on the scale of the wave packet is more pronounced. The result suggests there might be an appreciable difference in the dynamical properties of the two models. Figure 3. Electron scattering of a fixed proton with different impact parameters b computed using elongated wave packets, isotropic wave packets and SOFT. The trajectory of electron mass-centre (a) and the width of the electron wave function along y-direction, ⟨y2⟩−⟨y⟩2, for two different impact parameters (b). The wave functions were initialized as isotropic Gaussians with a width 1aB and velocity such that the total classical energy of each trajectory were (3/2) kBT and kBT=10 eV.

Figure 4. Transfer of (classical) kinetic energy ΔE/E from the electron to the ion in a single electron–ion scattering event as a function of impact parameter b. Initial velocity set as in figure 3 with a temperature (a) kBT=5 eV and (b) kBT=10 eV. The energy transferred differs between the isotropic and elongated models, especially for intermediate impact parameters and smaller wave packets. The result is shown for different free particle widths w set by the confining potential in equation (2.17) by varying lw with Aw=3 Eh/aB2.

(c) Transport properties

The previously shown demonstrations of the model were limited to the dynamics of a few particles; however, one of the strengths of the wave packet models is in the treatment of large collections of particles. We demonstrate this by studying a hydrogen plasma with degeneracy parameter θ=kBT/EF=1 and a density rs=2 under WDM conditions. The system is modelled by a thousand protons and an equal number of electron wave packets describing a spin un-polarized electron fluid (with an equal number of spin-up and spin-down electrons). For this initial test we set ρ=1 in equation (2.16) and the width of the wave packets are regularized by the width confinement where Aw=3 Eh/aB2 and lw=1aB. The initial random configurations of particles were allowed to thermalize under the influence of periodic velocity re-scaling (see appendix D) for 75 fs, after which data were collected for 225 fs, a procedure repeated five times for each wave packet model.

Figure 5a shows the static pair-correlation functions [27] for the two wave packet models under consideration, computed classically without accounting for antisymmetrization of the electrons, and such even equal spin electrons have a contribution at no separation where otherwise the antisymmetrization completely set this limit. The static structures primarily differ in electronic structure. The isotropic wave packets are seen to interact more strongly both with ions and between themselves. A part of this stronger interaction is during a collision between a wave packet and a proton. The isotropic wave function expands slower, as the appropriate expansion is averaged over all directions, and the smaller wave packets have a more localized charge and a stronger interaction. Figure 5. (a) Radial distribution function for proton–proton (p–p), electron–electron (e–e) and proton–electron (p–e) in both wave packet models described, elongated (solid) and isotropic (dashed). Qualitatively similar structures are seen between the two models, where for isotropic wave packets the electrons have a stronger interaction between themselves as well as the ions. The typical errors—estimated from the distinct runs—are shown in the shaded area only notable at small r. (b) Electron current auto correlation function for the wave packet model (solid), described well by an exponential decay shown by the corresponding fit (dashed). The system considered in both (a) and (b) is the one described in §4c.

One of the properties of interest for WDM systems is the electrical conductivity σ. The conductivity relates to the microscopic dynamics in an atomistic simulation via [93], 4.1 σ(ω)=13kB T V∫0∞⟨⟨J(t)⋅J(0)⟩⟩ eiωt dt,

where ⟨⟨⋅⟩⟩ is a thermal average and J is the total charge current. The current is dominated by the electron contribution, Je, 4.2 J≈Je=−em∑i=1Ne⟨p^i⟩,

where we sum over all electrons. The current–current correlation function has roughly an exponentially decaying behaviour and therefore in the vein of Mithen et al. [94] an exponential form is fitted to reduce the influence of noise. Both the numerical data and the fit are shown in figure 5 b with a high level of agreement. An exponential current–current correlation function results in a Drude-like conductivity, 4.3 σ(ω)=σ01+(ω/λ)2+i σ0ω/λ1+(ω/λ)2,

where σ0 and λ are related to the amplitude and time constant of the decaying correlation function. Physically, σ0 is the DC conductivity while λ corresponds to the inverse of the mean free scattering time. The wave packet models differ in their prediction of both constants. For elongated wave packets, σ0elo=(34 300±1 900) Ω−1 cm−1 and λelo=(1.48±0.06) fs−1, while for the isotropic model σ0iso=(29 800±1 200) Ω−1 cm−1 and λiso=(1.73±0.06) fs−1. This corresponds roughly to a 15% increase in DC conductivity and a 15% decrease in the scattering frequency as we extend the wave packet formulation. For a Drude-like conductivity, σ0 and λ are strongly correlated, which can be understood by the constraints set by the fluctuation–dissipation theorem [95]. Within this wave packet formulation the electron–electron scattering is explicitly included in the dynamic formulation, not the case in the commonly used Kubo–Greenwood formulation of conductivity for DFT-MD, preventing the DFT-based technique to achieve the appropriate high-temperature limit [96].

5. Conclusion

Non-equilibrium modelling of quantum many-body systems is a formidable task that further suffers from the large proton-to-electron mass ratio, resulting in vastly different time scales for the evolution of electron and ion dynamics. Any dynamical treatment must, therefore, resolve the electron motion, while extending over comparatively long-time scales to investigate ion dynamics. Wave packet models address this with an ansatz for the electronic wave function allowing the time evolution of a large number of particles to be performed over long-time scales, while retaining quantum mechanical properties dynamically within the model.

We have extended the functional form of the wave packets used for modelling WDM by allowing the wave packet to be elongated with arbitrary rotation. This in turn allows for a dynamic response to gradients across the wave packet that should better represent quantum dynamics. As a consequence of the non-isotropic states, explicit evaluation of the interaction Hamiltonian has not been possible, and a generalized Ewald summation has been used to appropriately evaluate both the short- and long-range effect of the Coulomb interaction. Furthermore, a decomposition of the short-range interaction kernel into Gaussian modes has been constructed and used for an explicit evaluation scheme.

Crucially, WDM systems are partially degenerate and the exchange interaction contributes significantly to their evolution. A scalable approximation for this in terms of Pauli potentials has been derived as they fundamentally depend on the wave packets used. The interactions are implemented along with the complete dynamical description in LAMMPS, with good parallel support and close to a linear scaling with particle number.

The elongated wave packet model is seen to improve the description of ground state hydrogen molecules compared to isotropic wave packets where spherical symmetry is naturally broken. The added degrees of freedom also further improve the dynamical description when compared to full quantum mechanical treatments, which we demonstrate for electron–proton scattering. Furthermore, in such collisions different energy transfers are observed between the models, illustrating an important variation in the predicted behaviour. Finally, the model is used to investigate a partially degenerate hydrogen plasma. Differences between wave packet models are seen both in the electronic structure and in the dynamics. A property of fundamental interest for this type of system is the electrical conductivity, which we extracted from our wave packet models. The electrical conductivity is dominated by electron motion and the DC conductivity was seen to increase by approximately 15% when extending the functional form of the wave packet.

Acknowledgements

The authors would like to thank S. Azadi and T. Gawne for valuable discussions.

Data accessibility

The data for all figures and a statically compiled version of the code are provided via the link: https://ora.ox.ac.uk/objects/uuid:5717ba1f-9bd4-4fd9-9a3e-04ad23ff4259 (DOI: https://dx.doi.org/10.5287/ora-eojjdyvw9).

Authors' contributions

P.S.: conceptualization, data curation, formal analysis, investigation, methodology, project administration, resources, software, validation, visualization, writing—original draft, writing—review and editing; T.C.: conceptualization, investigation, writing—review and editing; F.G.: conceptualization, investigation, writing—review and editing; Z.M.: data curation, methodology, writing—review and editing; N.L.: data curation; V.S.B.: data curation; S.R.: supervision; S.M.V.: conceptualization, funding acquisition, investigation, resources, supervision, writing—review and editing; G.G.: conceptualization, funding acquisition, investigation, resources, supervision, writing—review and editing. All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

This work was financially supported by AWE UK via Oxford Centre for High Energy Density Science (OxCHEDS) scholarship jointly founded by Oxford Physics Endowment for Graduates (OxPEG). Z.M. gratefully acknowledges funding by the Center for Advanced Systems Understanding (CASUS) which is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon state government out of the State budget approved by the Saxon State Parliament. V.S.B. acknowledges support from National Science Foundation Grant CHE-1900160 for the calculations performed by N.L.

Appendix A. Gaussian decomposition

The decomposition of the Coulomb interaction kernels described in §2b is solely introduced to efficiently evaluate state averages and only needed to be computed once. The state averages of the interaction kernel V can be reduced to a volume integration weighted by some Gaussian profile. Therefore, the error in Gaussian approximation was quantified by the volume average of the square errors A 1 L=∫0rmax dr r2|V(r)−∑pcp e−αpr2|2,

which is finite without a small r cut-off. In the case of the short-range contribution, V(r)=erfc(gr)/r, the integrals are convergent and rmax=∞ was chosen. On the other hand, for the long-range correction, V(r)=erf(gr)/r, a finite rmax is needed and has been chosen to rmax=15aB. The latter decomposition is only applied to terms in the Pauli interaction, see §2c, which are exponentially suppressed by the overlap |⟨qi|qj⟩|2 for large distances and the self-energy for the Ewald summation which is also exponentially suppressed for large distances, see appendix B. Therefore, a finite rmax does not constitute a significant limitation.

The decomposition in terms of the coefficients cp and αp is based on the minimization of L. For a fixed set of αp, the minimization is solved by inversion of the linear system of equations ∂L/∂cp=0. However, instead of using a particular set of αp [72], a general minimization algorithm was used to find the optimal αp’s, preconditioned with the optimal cp=cp({αp}). For the short-range interaction the search was limited to αp>0, where in the long-range case rmax−2/2<αp<g2 was needed for a stable decomposition. Figure 6 demonstrates the final decomposition resulting from this procedure described above, typical for the ones used in the main text. The singular nature of the short-range interaction is approximated by a series of modes with successively larger amplitude and shorter range. Figure 6. The Gaussian mode decomposition of the short-range interaction (a) and the long range interaction (b) with a Ewald parameter g=0.5.

Appendix B. Ewald summation

The main contribution Ek to the Ewald summation (2.11) is the interaction of the charge distribution in the simulation box ρuc and the total surrounding charge distribution ρtot, B 1 Ek=∫ d3x1 d3x2erf(g|x1−x2|)2|x1−x2|ρuc(x1) ρtot(x2),

seen by exchanging the particle summation and state average in equation (2.11). The charge distributions are based on the functional form of the wave packet and the image charges, namely B 2 ρuc(x)=∑ieZi|⟨x|qi⟩|2andρtot(x)=∑n∈Z3ρuc(x−Ln),

the Fourier transform of which, ρ~uc(k), has been evaluated in equation (2.13). With the use of Plancherel’s theorem and the Fourier convolution theorem, one can arrive at the expression (2.12), which converges rapidly in K-space due to the exponential dependence on k. The k=0 term vanishes due to charge neutrality ρ~uc(0)=0.

In any given computation only a finite number of k-vectors may be used and the associated error incurred by truncating equation (2.12) was described for the classical point particle case by Kolafa & Perram [69]. It can be shown that the resulting force from the Ewald summation for distributed particles is suppressed by factors exp⁡(−0.5k⊺Σik)×exp⁡(−0.5k⊺Σjk) for each k-vector and the force converges more rapidly—as the charge distribution now is more smooth—such that the classical estimate still can be used as an upper bound. The calculations in the manuscript were performed such that this error estimate was one thousand of the force two ions experience an average distance apart.

The self-energy term can be formulated similarly B 3 Es=−∑i(eZi)2 ∫ d3x1|⟨x1|qi⟩|2∫ d3x2 erf(g|x1−x2|)2|x1−x2||⟨x2|qi⟩|2,

where the integration could be performed if the kernel were Gaussian. The integrand is exponentially small if x1≉x2 and a Gaussian decomposition of the long-range interaction kernel is used, using erf(gr)/r=∑pc~pexp⁡(−α~pr2). The self-energy contribution is B 4 Es=−12∑i(eZi)2∑pc~pdet(I+4α~pΣi),

yielding the correct classical limit, Σi→0, as ∑pc~p=2g/π, if the decomposition retain the correct value at r=0.

Appendix C. Log-range Pauli interaction

The background terms, Vbg, in the Pauli potential (2.14b) have a long-range behaviour and by splitting the interaction in accordance with the Ewald summation discussed in §2b it may be truncated and any residual long-range term is then accounted for. Once again the long-range interaction can be formulated as the interaction between two distributions; one is the particle overlap and other is the charge density from the background particles C 1 VlongP/C=∫ d3x1 d3x2erf(g|x1−x2|)2|x1−x2|μP/C(x1)ρtot(x2)+VsP/C,

where C 2 μP/C(x) =±∑i<je1∓|⟨qi||qj⟩|2[2ℜ{⟨qj|x⟩⟨x|qj⟩⟨qj|qj⟩} −(|⟨x|qi⟩|2+|⟨x|qj⟩|2)|⟨qi||qj⟩|2]

and the sum is restricted to electrons with relevant spins. A self-energy term VsP/C is introduced when the background terms do not have terms where k=i,j, C 3 VsP/C =±∑i<j∫ d3x1 d3x2erf(g|x1−x2|)2|x1−x2|e21∓|⟨qi|qj⟩|2(|⟨x2|qi⟩|2+|⟨x2|qj⟩|2) ×[2ℜ{⟨qj|x1⟩⟨x1|qj⟩⟨qj|qj⟩}−(|⟨x1|qi⟩|2+|⟨x1|qj⟩|2)|⟨qi|qj⟩|2].

While the above expression has been evaluated in K-space and included in the force computation, it has only a limited impact on the overall dynamics, and it could be omitted if computational speed is needed.

Appendix D. Temperature measurements

The system temperature, T, is defined through the equipartition theorem, D 1 ⟨ ⟨Qμ∂H∂Qν⟩ ⟩=δμνkBT,

where δμν is the Kronecker delta and ⟨⟨⋅⟩⟩ represents the thermal average in the microcanonical ensemble commonly referred to as the NVE ensemble, due to a constant particle number, N, volume, V and energy E. This corresponds to a time average in molecular dynamics in accordance with the ergodic hypothesis [66]. We may define a temperature based on the classical degrees of freedom Qμ=Qν=pi, D 2 ⟨ ⟨pi2m+pi⊺∂V∂pi⟩ ⟩=3kBT,

and the internal ones, D 3 ⟨ ⟨4mTr{ΠiΣiΠi}+∑αβΠiαβ τβα ∂V∂Πiβα⟩ ⟩=6kBT,

where V is the state average of the interaction terms of the Hamiltonian. The above two definitions are combined, averaged over the particle index i and evaluated instantaneously to monitor the temperature evolution of the system. Instantaneously, equations (D 2) and (D 3) may differ, nonetheless converge to the same value by time averaging. Previously, Ma et al. [79] has incorporated the width kinetic energy in the temperature measurement for isotropic wave packets. When we employ isotropic wave packets we use equation (D 3) accounting for the reduced degrees of freedom.

Appendix E. Generalized friction

To achieve a continuous energy loss from friction, the equations of motion are modified accordingly E 1 dpidt=−∂H∂ri+ΓipandddtΠiαβ=−ταβ∂H∂Σiαβ+(ΓiΠ)αβ.

The resulting energy loss E 2 ddtH=−∑i(∂H∂pi⊺ Γip+∑αβταβ∂H∂Πiαβ(ΓiΠ)βα),

can be guaranteed to result in a monotonic energy decrease with the choice E 3 Γip=mγp ∂H∂pi(ΓiΠ)αβ=m8γΠ∑γ((Σi−1)αγτγβ∂H∂Πiγβ+ταγ∂H∂Πiαγ(Σi−1)γβ),

which reduces the classical case in the absence of a Pauli potential and ΣiΠi=ΠiΣi.
==== Refs
References

1. Nuckolls J, Wood L, Thiessen A, Zimmerman G. 1972 Laser compression of matter to super-high densities: thermonuclear (CTR) applications. Nature 239 , 139-142. (10.1038/239139a0)
2. Betti R, Hurricane O. 2016 Inertial-confinement fusion with lasers. Nat. Phys. 12 , 435-448. (10.1038/nphys3736)
3. Abu-Shawareb H et al. 2022 Lawson criterion for ignition exceeded in an inertial fusion experiment. Phys. Rev. Lett. 129 , 075001. (10.1103/PhysRevLett.129.075001)36018710
4. Smith RF et al. 2018 Equation of state of iron under core conditions of large rocky exoplanets. Nat. Astron. 2 , 452-458. (10.1038/s41550-018-0437-9)
5. Falk K. 2018 Experimental methods for warm dense matter research. High Power Laser Sci. Eng. 6 , e59. (10.1017/hpl.2018.53)
6. Frank A et al. 2019 Exoplanets and high energy density plasma science. Bull. AAS 51 , 36.
7. Levy A et al. 2015 The creation of large-volume, gradient-free warm dense matter with an X-ray free-electron laser. Phys. Plasmas 22 , 030703. (10.1063/1.4916103)
8. Fletcher L et al. 2015 Ultrabright X-ray laser scattering for dynamic warm dense matter physics. Nat. Photonics 9 , 274-279. (10.1038/nphoton.2015.41)
9. Schlanges M, Bonitz M, Tschttschjan A. 1995 Plasma phase transition in fluid hydrogen-helium mixtures. Contrib. Plasma Phys. 35 , 109-125. (10.1002/ctpp.2150350203)
10. Bezkrovniy V, Filinov V, Kremp D, Bonitz M, Schlanges M, Kraeft W, Levashov P, Fortov V. 2004 Monte Carlo results for the hydrogen Hugoniot. Phys. Rev. E 70 , 057401. (10.1103/PhysRevE.70.057401)
11. Vorberger J, Tamblyn I, Militzer B, Bonev S. 2007 Hydrogen-helium mixtures in the interiors of giant planets. Phys. Rev. B 75 , 024206. (10.1103/PhysRevB.75.024206)
12. Lorenzen W, Becker A, Redmer R. 2014 Progress in warm dense matter and planetary physics. In Frontiers and challenges in warm dense matter (eds Graziani F, Desjarlais MP, Redmer R, Trickey SB), pp. 203–234. Switzerland: Springer International Publishing.
13. Militzer B, González-Cataldo F, Zhang S, Driver KP, Soubiran F. 2021 First-principles equation of state database for warm dense matter computation. Phys. Rev. E 103 , 013203. (10.1103/PhysRevE.103.013203)33601631
14. Saumon D, Hubbard WB, Chabrier G, Van Horn HM. 1992 The role of the molecular-metallic transition of hydrogen in the evolution of Jupiter, Saturn, and brown dwarfs. Astrophys. J. 391 , 827-831. (10.1086/171391)
15. Chabrier G. 1993 Quantum effects in dense Coulumbic matter-Application to the cooling of white dwarfs. Astrophys. J. 414 , 695-700. (10.1086/173115)
16. Chabrier G, Brassard P, Fontaine G, Saumon D. 2000 Cooling sequences and color-magnitude diagrams for cool white dwarfs with hydrogen atmospheres. Astrophys. J. 543 , 216. (10.1086/317092)
17. Haensel P, Potekhin AY, Yakovlev DG. 2007 Neutron stars 1: equation of state and structure. New York, NY: Springer New York.
18. Daligault J, Gupta S. 2009 Electron–ion scattering in dense multi-component plasmas: application to the outer crust of an accreting neutron star. Astrophys. J. 703 , 994. (10.1088/0004-637X/703/1/994)
19. Hu S et al. 2018 A review on ab initio studies of static, transport, and optical properties of polystyrene under extreme conditions for inertial confinement fusion applications. Phys. Plasmas 25 , 056306. (10.1063/1.5017970)
20. Bonitz M, Dornheim T, Moldabekov ZA, Zhang S, Hamann P, Kählert H, Filinov A, Ramakrishna K, Vorberger J. 2020 Ab initio simulation of warm dense matter. Phys. Plasmas 27 , 042710. (10.1063/1.5143225)
21. Gaffney J et al. 2018 A review of equation-of-state models for inertial confinement fusion materials. High Energy Density Phys. 28 , 7-24. (10.1016/j.hedp.2018.08.001)
22. Davis RA, Angermeier W, Hermsmeier R, White TG. 2020 Ion modes in dense ionized plasmas through nonadiabatic molecular dynamics. Phys. Rev. Res. 2 , 043139. (10.1103/PhysRevResearch.2.043139)
23. Yao Y, Zeng Q, Chen K, Kang D, Hou Y, Ma Q, Dai J. 2021 Reduced ionic diffusion by the dynamic electron–ion collisions in warm dense hydrogen. Phys. Plasmas 28 , 012704. (10.1063/5.0028925)
24. Grabowski P et al. 2020 Review of the first charged-particle transport coefficient comparison workshop. High Energy Density Phys. 37 , 100905. (10.1016/j.hedp.2020.100905)
25. Wahl SM et al. 2017 Comparing Jupiter interior structure models to Juno gravity measurements and the role of a dilute core. Geophys. Res. Lett. 44 , 4649-4659. (10.1002/2017GL073160)
26. Hu S, Goncharov V, Boehly T, McCrory R, Skupsky S, Collins LA, Kress JD, Militzer B. 2015 Impact of first-principles properties of deuterium–tritium on inertial confinement fusion target designs. Phys. Plasmas 22 , 056304. (10.1063/1.4917477)
27. Ichimaru S. 1982 Strongly coupled plasmas: high-density classical plasmas and degenerate electron liquids. Rev. Mod. Phys. 54 , 1017. (10.1103/RevModPhys.54.1017)
28. Mithen JP, Daligault J, Gregori G. 2011 Extent of validity of the hydrodynamic description of ions in dense plasmas. Phys. Rev. E 83 , 015401. (10.1103/PhysRevE.83.015401)
29. Kählert H. 2020 Thermodynamic and transport coefficients from the dynamic structure factor of Yukawa liquids. Phys. Rev. Res. 2 , 033287. (10.1103/PhysRevResearch.2.033287)
30. Minoo H, Gombert M, Deutsch C. 1981 Temperature-dependent Coulomb interactions in hydrogenic systems. Phys. Rev. A 23 , 924. (10.1103/PhysRevA.23.924)
31. Hansen J, McDonald I. 1983 Thermal relaxation in a strongly coupled two-temperature plasma. Phys. Lett. A 97 , 42-44. (10.1016/0375-9601(83)90097-X)
32. Glosli J et al. 2008 Molecular dynamics simulations of temperature equilibration in dense hydrogen. Phys. Rev. E 78 , 025401. (10.1103/PhysRevE.78.025401)
33. Larder B, Gericke D, Richardson S, Mabey P, White T, Gregori G. 2019 Fast nonadiabatic dynamics of many-body quantum systems. Sci. Adv. 5 , eaaw1634. (10.1126/sciadv.aaw1634)31803829
34. Lambert F, Clérouin J, Mazevet S, Gilles D. 2007 Properties of hot dense plasmas by orbital-free molecular dynamics. Contrib. Plasma Phys. 47 , 272-280. (10.1002/ctpp.200710037)
35. White T, Richardson S, Crowley B, Pattison L, Harris J, Gregori G. 2013 Orbital-free density-functional theory simulations of the dynamic structure factor of warm dense aluminum. Phys. Rev. Lett. 111 , 175002. (10.1103/PhysRevLett.111.175002)24206498
36. Collins L, Kwon I, Kress J, Troullier N, Lynch D. 1995 Quantum molecular dynamics simulations of hot, dense hydrogen. Phys. Rev. E 52 , 6202. (10.1103/PhysRevE.52.6202)
37. Rüter HR, Redmer R. 2014 Ab initio simulations for the ion-ion structure factor of warm dense aluminum. Phys. Rev. Lett. 112 , 145007. (10.1103/PhysRevLett.112.145007)24765982
38. Moldabekov ZA, Bonitz M, Ramazanov T. 2018 Theoretical foundations of quantum hydrodynamics for plasmas. Phys. Plasmas 25 , 031903. (10.1063/1.5003910)
39. Michta D. 2020 Quantum hydrodynamics: theory and computation with applications to charged particle stopping in warm dense matter. PhD thesis, Princeton University, Princeton.
40. Moldabekov Z, Dornheim T, Gregori G, Graziani F, Bonitz M, Cangi A. 2022 Towards a quantum fluid theory of correlated many-fermion systems from first principles. SciPost Phys. 12 , 062. (10.21468/SciPostPhys.12.2.062)
41. Baczewski AD, Shulenburger L, Desjarlais M, Hansen S, Magyar R. 2016 X-ray Thomson scattering in warm dense matter without the Chihara decomposition. Phys. Rev. Lett. 116 , 115004. (10.1103/PhysRevLett.116.115004)27035307
42. Magyar RJ, Shulenburger L, Baczewski AD. 2016 Stopping of deuterium in warm dense deuterium from Ehrenfest time-dependent density functional theory. Contrib. Plasma Phys. 56 , 459-466. (10.1002/ctpp.201500143)
43. Frydrych S et al. 2020 Demonstration of X-ray Thomson scattering as diagnostics for miscibility in warm dense matter. Nat. Commun. 11 , 1-7. (10.1038/s41467-020-16426-y)31911652
44. Militzer B, Ceperley DM. 2000 Path integral Monte Carlo calculation of the deuterium Hugoniot. Phys. Rev. Lett. 85 , 1890. (10.1103/PhysRevLett.85.1890)10970640
45. Hu SX, Militzer B, Goncharov VN, Skupsky S. 2011 First-principles equation-of-state table of deuterium for inertial confinement fusion applications. Phys. Rev. B 84 , 224109. (10.1103/PhysRevB.84.224109)
46. Chapman D, Gericke D. 2011 Analysis of Thomson scattering from nonequilibrium plasmas. Phys. Rev. Lett. 107 , 165004. (10.1103/PhysRevLett.107.165004)22107396
47. Zastrau U et al. 2014 Resolving ultrafast heating of dense cryogenic hydrogen. Phys. Rev. Lett. 112 , 105002. (10.1103/PhysRevLett.112.105002)24679300
48. White T et al. 2014 Electron-ion equilibration in ultrafast heated graphite. Phys. Rev. Lett. 112 , 145005. (10.1103/PhysRevLett.112.145005)24765980
49. Plagemann KU, Rüter HR, Bornath T, Shihab M, Desjarlais MP, Fortmann C, Glenzer SH, Redmer R. 2015 Ab initio calculation of the ion feature in X-ray Thomson scattering. Phys. Rev. E 92 , 013103. (10.1103/PhysRevE.92.013103)
50. Clérouin J, Robert G, Arnault P, Ticknor C, Kress JD, Collins LA. 2015 Evidence for out-of-equilibrium states in warm dense matter probed by X-ray Thomson scattering. Phys. Rev. E 91 , 011101. (10.1103/PhysRevE.91.011101)
51. Mabey P, Richardson S, White T, Fletcher L, Glenzer S, Hartley N, Vorberger J, Gericke D, Gregori G. 2017 A strong diffusive ion mode in dense ionized matter predicted by Langevin dynamics. Nat. Commun. 8 , 1-6. (10.1038/ncomms14125)28232747
52. Feldmeier H, Schnack J. 2000 Molecular dynamics for fermions. Rev. Mod. Phys. 72 , 655. (10.1103/RevModPhys.72.655)
53. Grabowski PE. 2014 A review of wave packet molecular dynamics. In Frontiers and challenges in warm dense matter (eds F Graziani, MP Desjarlais, R Redmer, SB Trickey), pp. 265–282. Cham, Switzerland. Springer International Publishing.
54. Plimpton S. 1995 Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117 , 1-19. (10.1006/jcph.1995.1039)
55. Heller EJ. 1975 Time-dependent approach to semiclassical dynamics. J. Chem. Phys. 62 ,1544-1555. (10.1063/1.430620)
56. Corbin N, Singer K. 1982 Semiclassical molecular dynamics of wave packets. Mol. Phys. 46 , 671-677. (10.1080/00268978200101511)
57. Broeckhove J, Lathouwers L, Van Leuven P. 1989 Time-dependent variational principles and conservation laws in wavepacket dynamics. J. Phys. A: Math. Gen. 22 , 4395. (10.1088/0305-4470/22/20/016)
58. Knaup M, Reinhard P, Toepffer C, Zwicknagel G. 2003 Wave packet molecular dynamics simulations of warm dense hydrogen. J. Phys. A: Math. Gen. 36 , 6165. (10.1088/0305-4470/36/22/344)
59. Jakob B, Reinhard PG, Toepffer C, Zwicknagel G. 2007 Wave packet simulation of dense hydrogen. Phys. Rev. E 76 , 036406. (10.1103/PhysRevE.76.036406)
60. Jakob B, Reinhard P, Toepffer C, Zwicknagel G. 2009 Wave packet simulations for the insulator–metal transition in dense hydrogen. J. Phys. A: Math. Theor. 42 , 214055. (10.1088/1751-8113/42/21/214055)
61. Klakow D, Toepffer C, Reinhard PG. 1994 Hydrogen under extreme conditions. Phys. Lett. A 192 , 55-59. (10.1016/0375-9601(94)91015-4)
62. Klakow D, Toepffer C, Reinhard PG. 1994 Semiclassical molecular dynamics for strongly coupled Coulomb systems. J. Chem. Phys. 101 , 10 766-10 774. (10.1063/1.467889)
63. Murillo M, Timmermans E. 2003 Self-similar wavepackets for describing warm dense matter. Contrib. Plasma Phys. 43 , 333-337. (10.1002/ctpp.200310040)
64. Ono J, Ando K. 2012 Semiquantal molecular dynamics simulations of hydrogen-bond dynamics in liquid water using multi-dimensional Gaussian wave packets. J. Chem. Phys. 137 , 174503. (10.1063/1.4762840)23145735
65. Wyatt RE. 2005 Quantum dynamics with trajectories: introduction to quantum hydrodynamics, vol. 28 . New York, NY: Springer Science & Business Media.
66. Rapaport DC. 2004 The art of molecular dynamics simulation. Cambridge, UK: Cambridge University Press.
67. de Leeuw SW, Perram JW, Smith ER. 1980 Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric constants. Proc. R. Soc. Lond. A 373 , 27-56. (10.1098/rspa.1980.0135)
68. Deserno M, Holm C. 1998 How to mesh up Ewald sums. I. A theoretical and numerical comparison of various particle mesh routines. J. Chem. Phys. 109 , 7678-7693. (10.1063/1.477414)
69. Kolafa J, Perram JW. 1992 Cutoff errors in the Ewald summation formulae for point charge systems. Mol. Simul. 9 , 351-368. (10.1080/08927029208049126)
70. Knaup M, Reinhard PG, Toepffer C. 1999 Wave packet molecular dynamics simulations of hydrogen near the transition to a metallic fluid. Contrib. Plasma Phys. 39 , 57-60. (10.1002/ctpp.2150390114)
71. Frantsuzov PA, Mandelshtam VA. 2004 Quantum statistical mechanics with Gaussians: equilibrium properties of van der Waals clusters. J. Chem. Phys. 121 , 9247-9256. (10.1063/1.1804495)15538845
72. Shigeta Y, Inui T, Baba T, Okuno K, Kuwabara H, Kishi R, Nakano M. 2013 Quantal cumulant mechanics and dynamics for multidimensional quantum many-body clusters. Int. J. Quantum Chem. 113 , 348-355. (10.1002/qua.24052)
73. Buch V. 2002 Exploration of multidimensional variational Gaussian wave packets as a simulation tool. J. Chem. Phys. 117 , 4738-4750. (10.1063/1.1497968)
74. Su JT, Goddard III WA. 2007 Excited electron dynamics modeling of warm dense matter. Phys. Rev. Lett. 99 , 185003. (10.1103/PhysRevLett.99.185003)17995416
75. Su JT, Goddard III WA. 2009 The dynamics of highly excited electronic systems: applications of the electron force field. J. Chem. Phys. 131 , 244501. (10.1063/1.3272671)20059073
76. Jaramillo-Botero A, Su J, Qi A, Goddard III WA. 2011 Large-scale, long-term nonadiabatic electron molecular dynamics for describing material properties and phenomena in extreme environments. J. Comput. Chem. 32 , 497-512. (10.1002/jcc.21637)20812325
77. Kim H, Su JT, Goddard WA. 2011 High-temperature high-pressure phases of lithium from electron force field (eFF) quantum electron dynamics simulations. Proc. Natl Acad. Sci. USA 108 , 15 101-15 105. (10.1073/pnas.1110322108)
78. Xiao H, Jaramillo-Botero A, Theofanis PL, Goddard III WA. 2015 Non-adiabatic dynamics modeling framework for materials in extreme conditions. Mech. Mater. 90 , 243-252. (10.1016/j.mechmat.2015.02.008)
79. Ma Q, Dai J, Kang D, Murillo M, Hou Y, Zhao Z, Yuan J. 2019 Extremely low electron-ion temperature relaxation rates in warm dense hydrogen: interplay between quantum electrons and coupled ions. Phys. Rev. Lett. 122 , 015001. (10.1103/PhysRevLett.122.015001)31012692
80. Liu Y, Liu X, Zhang S, Liu H, Mo C, Fu Z, Dai J. 2021 Molecular dynamics investigation of the stopping power of warm dense hydrogen for electrons. Phys. Rev. E 103 , 063215. (10.1103/PhysRevE.103.063215)34271766
81. Angermeier WA, White TG. 2021 An investigation into the approximations used in wave packet molecular dynamics for the study of warm dense matter. Plasma 4 , 294-308. (10.3390/plasma4020020)
82. Knaup M, Reinhard PG, Toepffer C. 2001 Wave packet molecular dynamics simulations of deuterium in the region of laser shock-wave experiments. Contrib. Plasma Phys. 41 , 159-162. (10.1002/1521-3986(200103)41:2/3<159::AID-CTPP159>3.0.CO;2-P)
83. Ebeling W, Filinov A, Bonitz M, Filinov V, Pohl T. 2006 The method of effective potentials in the quantum-statistical theory of plasmas. J. Phys. A: Math. Gen. 39 , 4309. (10.1088/0305-4470/39/17/S01)
84. Morozov I, Valuev I. 2009 Localization constraints in Gaussian wave packet molecular dynamics of nonideal plasmas. J. Phys. A: Math. Theor. 42 , 214044. (10.1088/1751-8113/42/21/214044)
85. Grabowski PE, Markmann A, Morozov IV, Valuev IA, Fichtl CA, Richards DF, Batista VS, Graziani FR, Murillo MS. 2013 Wave packet spreading and localization in electron-nuclear scattering. Phys. Rev. E 87 , 063104. (10.1103/PhysRevE.87.063104)
86. Lavrinenko YS, Morozov I, Valuev I. 2016 Reflecting boundary conditions for classical and quantum molecular dynamics simulations of nonideal plasmas. Contrib. Plasma Phys. 56 ,448-458. (10.1002/ctpp.201500139)
87. Ebeling W, Militzer B. 1997 Quantum molecular dynamics of partially ionized plasmas. Phys. Lett. A 226 , 298-304. (10.1016/S0375-9601(96)00948-6)
88. Leimkuhler B, Reich S. 2004 Simulating hamiltonian dynamics, vol. 14 . Cambridge, UK: Cambridge university press.
89. Gabriel E et al. 2004 Open MPI: goals, concept, and design of a next generation MPI implementation. In Proc. 11th European PVM/MPI Users’ Group Meeting, Budapest, Hungary, pp. 97–104. Berlin, Germany: Springer.
90. Feit M, Fleck Jr J, Steiger A. 1982 Solution of the Schrödinger equation by a spectral method. J. Comput. Phys. 47 , 412-433. (10.1016/0021-9991(82)90091-2)
91. Markmann A, Graziani F, Batista VS. 2012 Kepler predictor–corrector algorithm: scattering dynamics with one-over-R singular potentials. J. Chem. Theory Comput. 8 , 24-35. (10.1021/ct200452h)26592868
92. Greene SM, Batista VS. 2017 Tensor-train split-operator fourier transform (TT-SOFT) method: multidimensional nonadiabatic quantum dynamics. J. Chem. Theory Comput. 13 , 4034-4042. (10.1021/acs.jctc.7b00608)28763215
93. Vashishta P, Kalia RK, Nakano A, Rino JP. 2008 Interaction potentials for alumina and molecular dynamics simulations of amorphous and liquid alumina. J. Appl. Phys. 103 , 083504. (10.1063/1.2901171)
94. Mithen J, Daligault J, Gregori G. 2012 Molecular dynamics simulations for the shear viscosity of the one-component plasma. Contrib. Plasma Phys. 52 , 58-61. (10.1002/ctpp.201100050)
95. Kubo R. 1966 The fluctuation-dissipation theorem. Rep. Prog. Phys. 29 , 255. (10.1088/0034-4885/29/1/306)
96. French M, Röpke G, Schörner M, Bethkenhagen M, Desjarlais MP, Redmer R. 2022 Electronic transport coefficients from density functional theory across the plasma plane. Phys. Rev. E 105 , 065204. (10.1103/PhysRevE.105.065204)35854489