Development of a new quantum trajectory molecular dynamics framework

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’.


Introduction
The establishment of high power lasers facilities during the last decades has been instrumental in the achievements towards inertial confinement fusion (ICF) [1][2][3], but also for the creation of highdensity 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][10][11][12][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 ionion interactions [27][28][29], classical electrons with effective quantum statistical potentials (QSP) [30][31][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][39][40], time-dependent DFT [41][42][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][47][48][49][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 formulationsapplicable 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.

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, whereĤ 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 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 Fermions are described by states antisymmetric under exchange and Slater-determinants have been considered in [58][59][60]. However, this approach scales unfavourably with particle number N. Instead, here we employ a product state |Q = |q 1 ⊗ |q 2 ⊗ · · · ⊗ |q N , (2.4) of single-particle orbitals, |q i , 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 |q i 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 where ξ i = x − r i . The wave packet is parameterized by 18 degrees of freedom, the position r i , momentum p i 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]. where Σ iαβ = (Σ i ) αβ and Π iαβ = (Π i ) αβ are the components of the symmetric matrices. The prefactor τ αβ 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 R I , momentum P I , charge Z I e and mass M I , and quantum electrons with positionx i and momentump i operators as well as charge −e and mass m. The system is described by the Hamiltonian the state average of which is required for the time evolution. The average kinetic energy 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] 1 r = erfc(gr) r + erf(gr) r , (2.9) 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, erfc(gr) r p c p e −α p r 2 , (2.10) where the coefficients c p 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 r i + Ln for all n ∈ Z 3 and L being the length of a cubic simulation cell. The interaction energy is 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 E k and the self-energy E s . In appendix B, we evaluate E k in reciprocal space to be and k = 2π n/L. Equation (2.12) converges rapidly due to the exponential factor. The self-energy term, E s , is independent of r i 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][76][77][78][79][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 with a background interaction from all other particles in the system 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 written here in terms of single-particle orbitals |q i . For equal spin particles only the spatially antisymmetric state is allowed and the Pauli potential V P ij is the difference between Ĥ 2 for the triplet and the product state |q i ⊗ |q j . In the case of opposite spin particles, the spatial state depends on the spin structure, but along the lines  [81] a correlation potential is introduced based on the singlet state, ρV C ij , multiplied by the parameter ρ. Therefore, the potentials are which can be shown to scale as | q i |q j | 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 V bg 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 |q j → |q i 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][83][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 where σ 2 i,α is the αth eigenvalue of Σ i and θ (x) is the Heaviside step function. The parameters l w and A w set the width of a free particle σ free (l w , A w ) l w by balancing the shape-kinetic energy. This potential is rotationally invariant and acts only on wave packets with a width larger than l w . Furthermore, the confinement reduces to the commonly used potential based on a harmonic potential centred at the particle position in the limit of l w → 0. In this specific limit, the potential has also been used to address the heat capacity in the classical limit [87].

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.
Furthermore, figure 1 a also demonstrates the scaling of computational cost with particle number for a test system where r s = (4π/3na 3 B ) −1/3 = 2 where n is the proton number density and a B 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.

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 (H 2 ), where the ability for the wave packet to stretch is believed to be crucial for molecular binding [53].  [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 ground state energy of H 2 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.8a B , 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.8a B , compared with our model at δ ≈ 2.8a B .
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, H 2 , 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][91][92]. The resulting trajectory data for impact parameters in the range of 0.6-3.0a B 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.

(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 θ = k B T/E F = 1 and a density r s = 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 A w = 3 E h /a 2 B and l w = 1a B . 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. 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], where · is a thermal average and J is the total charge current. The current is dominated by the electron contribution, J e , 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, 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, σ elo 0 = (34 300 ± 1 900) Ω −1 cm −1 and λ elo = (1.48 ± 0.06) fs −1 , while for the isotropic model σ iso 0 = (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].

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.

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 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 r max = ∞ was chosen. On the other hand, for the long-range correction, V(r) = erf(gr)/r, a finite r max is needed and has been chosen to r max = 15a B . The latter decomposition is only applied to terms in the Pauli interaction, see §2c, which are exponentially suppressed by the overlap | q i |q j | 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 r max does not constitute a significant limitation. The decomposition in terms of the coefficients c p 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/∂c p = 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 c p = c p ({α p }). For the shortrange interaction the search was limited to α p > 0, where in the long-range case r −2 max /2 < α p < g 2 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.

Appendix B. Ewald summation
The main contribution E k 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 , 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 Σ i k) × exp(−0.5k Σ j k) 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 where the integration could be performed if the kernel were Gaussian. The integrand is exponentially small if x 1 ≈ x 2 and a Gaussian decomposition of the long-range interaction kernel is used, using erf(gr)/r = pc p exp(−α p r 2 ). The self-energy contribution is yielding the correct classical limit, Σ i → 0, as pc p = 2g/ √ π, if the decomposition retain the correct value at r = 0.
The background terms, V bg , 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 where μ P/C (x) = ± i<j e 1 ∓ | q i ||q j | 2 2 q j |x x|q j q j |q j − | x|q i | 2 + | x|q j | 2 | q i ||q j | 2 (C 2) and the sum is restricted to electrons with relevant spins. A self-energy term V P/C s is introduced when the background terms do not have terms where k = i, j, 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, 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 ν = p i , and the internal ones, 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.