An investigation into the approximations used in wave packet molecular dynamics for the study of warm dense matter

Wave packet molecular dynamics (WPMD) has recently received a lot of attention as a computationally fast tool to study dynamical processes in warm dense matter beyond the Born-Oppenheimer approximation. These techniques, typically, employ many approximations to achieve computational efficiency while implementing semi-empirical scaling parameters to retain accuracy. We investigate three of the main approximations ubiquitous to WPMD: a restricted basis set, approximations to exchange, and the lack of correlation. We examine each of these approximations in atomic and molecular hydrogen in addition to a dense hydrogen plasma. We find that the biggest improvement to WPMD comes from combining a two Gaussian basis with a semi-empirical correction based on the valence-bond wave function. A single parameter scales this correction to match experimental pressures of dense hydrogen. Ultimately, we find that semi-empirical scaling parameters are necessary to correct for the main approximations in WPMD. However, reducing the scaling parameters for more ab-initio terms gives more accurate results and displays the underlying physics more readily.


Introduction
Warm Dense Matter (WDM) is a critically important physical regime that bridges the gap between condensed matter and classical plasma physics. The WDM state is found in several astrophysical environments (e.g., planetary interiors and white dwarfs) [1,2]. It also has practical applications for understanding controlled thermonuclear fusion, and material processing [3]. Typically described as a system of strongly coupled ions immersed in a degenerate electron sea, WDM may exist in either a compressed liquid or a highly excited solid state. In both states, the ions have a Coulomb energy comparable to the thermal energy, while the electrons, at temperatures below the Fermi temperature, exhibit strong quantum behavior [4]. Techniques that simulate WDM states must model the slow and long-time behavior of the strongly coupled ions while simultaneously capturing the electrons' quantum mechanical nature. These inherent complexities lead to the failure of perturbative techniques resulting in differences in predictions of important quantities. While different models generally agree on the thermodynamic and acoustic properties [5,6], important quantities such as transport coefficients can differ by up to an order-of-magnitude [7,8].
Atomistic models in which the ions, treated through classical molecular dynamics, are coupled with a quantum mechanical treatment of the electrons have had the most success. The ion trajectories from such simulations can provide transport properties, such as viscosity and thermal diffusivity [9], acoustic properties, such as the sound speed [10,11], and thermodynamic variables, including the equation of ion-ion correlation in dense aluminum plasmas at temperatures of a few electronvolts [5,39]. While there has been some effort to understand the effect of a pair-wise exchange in dense hydrogen [22][23][24][25]40,41], the basis set limitation has not previously been investigated. However, work involving a multiple-Gaussian basis has been used to investigate wave-packet spreading in electron-nuclear scattering [33] and ionization of a single hydrogen atom [42,43]; where improvements up to a five Gaussian basis were found. Figure  1 shows improvement to the ground state energy of the hydrogen atom with an increasing number of Gaussians in the electron basis; minimal improvement is observed beyond four Gaussians.  [44]. The insert shows the changes to the shape of the wave functions with an increasing number of basis functions. Beyond four Gaussians we find negligible improvement in the either the ground state energy or shape of the wavefunction.
We will focus our efforts on dense hydrogen, the prototypical test-bed for atomistic models, and investigate the accuracy of WPMD as the number of Gaussians in the basis is increased. For computational efficiency, we retain a pairwise Pauli potential; however, based on the work of H. Xiao [45], we extend the Pauli potential to include additional potential energy terms. Finally, a simple correction based on the Valence Bond (VB) wave function is introduced with a single scaling parameter to address the model's lack of correlation and pairwise exchange. With this correction and a two-gaussian basis, we are able to exactly match the low-temperature pressure curve of dense hydrogen. Ultimately, we hope to extend the parameter space where WPMD, and specifically eFF, is applicable. This paper is organized as follows. Section 2 details the theory and development of the WPMD equations for systems with a multiple Gaussian basis. Particular focus will be given to the development of the updated Pauli potential, including a comparison of the updated Pauli potential with past work. Section 3 details the improvements afforded to hydrogen-based systems. The results include geometry optimization, potential energy surfaces and dynamics of the hydrogen molecule, and hydrogen under high pressure. Where available, we benchmark our calculations to experimental results. Finally, section 4 details some of the effects of approximations in current WPMD techniques and suggests a path forward to provide large-scale simulations of dense plasmas in previously unexplored regions of phase space.

Introduction to Wave-packet Molecular Dynamics
In WPMD the many-electron trial wave function Ψ is parameterized by a set of variables (q(t)) [22][23][24][25]42]. Development of a fully anti-symmetric wave function requires a calculation of the Slater determinant of all single electron wave functions, a computationally expensive operation [24,42,43]. To decrease computational cost, many implementations of WPMD, including eFF, use the simpler Hartree product, constructed through the linear combination of the single-particle wave-functions, ψ k . The Hartree product wave function, Ψ H , is defined as: where X = { x 1 , ..., x k , ..., x N e }, and x k are used to indicate the space occupied by the k-th electron. Most WPMD techniques utilize a single isotropic Gaussian as a restricted wave function for each electron [22][23][24]32,34]. This wave function is typically paramaterized by a set of ten real physical variables q = { r, p, σ, p σ } i.e., The elegance of this definition is that r = r is the expectation of position, p = p is the expectation of momentum, and σ = r 2 − | r | 2 is the uncertainty in position with the corresponding conjugate momentum, p σ .
Equations-of-motion for the dynamical parameters are easily derived from variation of the time-dependent Schrodinger equation [31,33,46], where H is the total energy of the system and q is the set of all dynamic variables. The norm matrix, N, is defined as follows, where N ab is a matrix element of the norm matrix, q a represents a specific element of the set of time-dependent variational parameters, and q * a is the complex conjugate of q a . For the single Gaussian basis given in Equation 2 the equations-of-motion take on a simple Hamilton form [31].
The total Hamiltonian operator of the semi-classical many-electron system interacting with classical point-like ions is given byĤ where the operators take on their usual classical, and quantum mechanical, definitions. Here,T e is the electron kinetic energy operator,V ei is the electron-ion Coulombic potential energy operator,V ee is the electron-electron Coulombic potential energy operator,T i is the ion kinetic energy, andV ii is the ion-ion Coulombic potential energy. In these expressions, p I represents the classical momentum of the I-th ion, M I is the mass of the I-th ion, Z J is the number of protons in the J-th ion, and R J represents the position of the J-th ion. The final term in the Hamiltonian defined in Equation (5) is a harmonic energy term ubiquitous throughout WPMD. This term is typically used to constrain the electron's size, which may increase to the point that electron-ion interactions become negligible. For this potential we have used the form suggested by Zwicknagel et al. where H Harm = Ĥ Harm = ∑ N e k=1 9 8γ 4 0 σ 2 and γ 0 is set to to half the simulation box length [35,47]. It should be noted that this term represents a small contribution to the total energy as, in the WDM regime, most wavepackets are constrained enough from the ionic potentials present [5,36].

Extension to Multiple Gaussians
The primary aim of this work was to investigate the the improvements in describing dense plasmas with WPMD when the basis set is extended to include multiple Guassians. Following the framework of Morozov and Valuev [42,43] we extend the basis to include multiple Gaussians as follows, where ψ k represents the single particle wave function of the k-th electron, N g is the total number of Gaussian wave packets per electron, and n k = ∑ α,β ϕ * kα ϕ kβ d 3 x is the normalization factor. The term ϕ kα represents the Gaussian wave packet α in the k-th electron. For multiple Gaussians the simple relationship between the parameters used in Equation (2) and the expectation of the electron physical characteristics is lost. Thus, to simplify the analytic derivation of the energy terms in the Hamiltonian we use the following Gaussian representation, Here, the set of dynamical variables for each Gaussian are q kα = {a kα , b kα , d kα }. These five complex parameters provide a total of ten real dynamic variables for each wave packet. If necessary, these can be mapped directly to the ten physical parameters used in Equation (2). The parameter, c kα , is an not an independent variable and is used to ensure normalization of each Gaussian [43]. Within this framework of multiple Gaussians, the Hartree energy of the system may be expressed as the sum of the ion and electron kinetic energies, along with the electron-ion, electron-electron, and ion-ion potential energies: In each case, an analytical expression may be derived. For the semi-classical electron terms these are most easily expressed as the product of an overlap integral with a residual. i.e., where O kαlδ represents the overlap between the α Gaussian wave packet of the k-th electron with the δ Gaussian wave packet of the l-th electron. Expressions for the three residual terms T e kαkγ , V ei Ikαkγ , and V ee kαlβkγlδ and the overlap integral are easily derived in the Gaussian basis and are provided in Ref. [43].

Development of a Pairwise Pauli and Correlation Potential
The Hartree product defined in Equation (1) neglects exchange effects captured by the Slater determinant. Within this approximation, important effects necessary to describe a quantum mechanical system of interacting fermions, such as the Pauli exclusion principle, are neglected. Here we detail the development of a spatially anti-symmetrized pairwise exchange term added between electrons of like spin. This term is equal to the difference between the energy calculated with the Slater determinant and that calculated with a Hartree product. In line with the eFF method, we retain a pairwise Slater determinant to ensure a computational scaling comparable with classical techniques. We construct our pairwise Pauli potential from three terms, where, are the Pauli kinetic, Pauli electron-ion potential, and Pauli electron-electron potential energy terms, respectively. Using the usual definition for a two-particle Slater determinant, analytic expressions for the two-particle Pauli energy given in Equations (15)(16)(17) were derived. The majority of WPMD techniques that do not implement full exchange have opted to use a Pauli potential based solely on the kinetic energy component of the Pauli exchange [22,23,32,34]. In addition, while some authors retain the dependence on electron momentum [22], many models, including eFF, simplify the terms further by neglecting this dependence [25,34].  [25]. Due to the lack of ions in the system, our total Pauli term contains just one ee . This term acts to lowers the exchange energy when compared to both the Klakow potentials. The eFF potential, which lies above the Klakow model, performs quite poorly for this system. This is unsurprising as, with a simple form, it is attempting to correct for several approximations in the model.
In Figure 2b, we compare the energy of the anti-bonding of the H 2 molecule with that predicted by unrestricted Hartree Fock (UHF). For the anti-bonding case, there is no correlation, which makes it the most straightforward system to examine exchange calculations. It is clear that the kinetic energy term is dominant and primarily responsible for the Pauli exclusion principle; however, the other terms are not negligible and play an important role in stabilizing molecules and preventing Gaussian coalescence [45]. When all three Pauli energy terms are added together, the result matches the exact UHF result and establishes our Pauli potential's validity. Note that the eFF potential, denoted by the solid black line, does not match the exact UHF result.
One of the reasons for the success of the eFF technique over previous models is the inclusion of additional scaling parameters in the Pauli potential. These parameters, matched to a set of molecular test structures, account for the lack of full exchange, the limited basis set, and the neglect of correlation, where we define correlation as the difference between the exact energy and the Hartree-Fock (HF) energy. Despite our expanded basis, two approximations persist in our model; these are pairwise exchange and lack of correlation. The inclusion of correlation in such models is difficult and decades old problem in quantum mechanical many-body systems [48][49][50][51]. Motivated by eFF, we used a valence bond (VB) wave function to develop a simple pairwise correction term to account for these deficiencies. The approximation takes on a similar form as the Pauli potential, where As before, utilizing the usual definition of the two-particle VB wave function, we derived analytic expressions for the pairwise VB correction energy given in Equations (20)(21)(22). The VB correction was implemented between pairs of electrons with the same spin. Opposite spin electrons were assumed to be well separated spatially, due to the Pauli potential, and to not contribute significantly to this term [52]. The total energy of our system, H = H H + H P + H C may be written as follows, where δ ↑↑ equaled unity when the spins of the electrons were parallel and zero otherwise, while δ ↑↓ was set to zero when the spins of the electrons are parallel and unity when anti-parallel. The parameter, ρ, must be chosen a-priori for the system as interest. For a hydrogen molecule, a value of ρ = 1 approaches the exact potential energy surface. For calculating the pressure of dense hydrogen, we scaled this value to match experimental results. We note that the size of ρ as N g → ∞ gives an indication of the remaining approximations in the model.

Implementation of Periodic Boundary Conditions
The simulation of dense plasmas necessitates the use of a finite size simulation box utilizing periodic boundary conditions; this brings with it additional challenges. The first is the evaluation of long-range forces associated with electrostatic interactions. We utilize the same scheme as eFF where the electrostatic energies are multiplied by a seventh order spline that goes from 1 to 0 over a radial distance r cut , defined so that the first, second, and third derivatives at the endpoints are zero [32,53]. i.e., where x = x/r cut . The cutoff range, r cut , is chosen to be equal to half the simulation box length. This technique enables the use of the minimum-image convention, retaining one of the most desirable properties of WPMD, which is the ability to simulate large systems of particles. Finally, as demonstrated in Figure 2a, the newly defined exchange terms do not exhibit the same long-range characteristics as the electrostatic terms and thus are not multiplied by the spline.
When applying the minimum image convention to the electrons, the Gaussians that comprise that electron must be treated as a single particle; that is to say, they must be shifted together. Thus, to apply the minimum image convention, we use the expectation of each electron's position. Suppose care is not taken, and the Gaussians are individually shifted. In that case, the unphysical situation where a particle interacts with the same electron on both the left and the right can arise.
Finally, the Pauli and VB correction electron-ion energy terms are not genuine pairwise terms. Each term in the summation involves two electrons plus an ion, essentially making it a three-body potential. To ensure consistent calculation of this term, we periodically shift the ion location to position it within half of the box length of the mid-point between the two electrons. We note that this differs from consistently shifting the ion towards the electron only when the two electrons are spatially separated; in such situations, the Pauli and VB correction energy terms are negligible.

Results
Hydrogen plasmas are one of the simplest physical systems and, for this reason, have become the prototypical test-bed for atomistic models of dense plasmas. For example, the hydrogen atom has a well-defined analytical solution. Figure 1 demonstrates that we approach this analytic solution with an increasing Gaussian basis. The ground state energy can be obtained to within a percent, utilizing a four Gaussian basis. Beyond this number, we observe only minor changes in the ground state energy and the electron wavefunction shape. It should also be noted that the most significant improvement to both occur when increasing the basis from one to two Gaussians.
We now turn our attention to modeling the energy of the hydrogen molecule. In Figure 3 we plot the potential energy surface of H 2 and show improvement in representation as the Gaussian basis is increased. Figure 3a displays the potential energy curve calculated without including the VB correction; in this case, the binding energy approaches the unrestricted Hartree-Fock (UHF) result. In this case, negligible improvement is observed beyond six Gaussians, slightly more than needed to describe the hydrogen atom accurately. Figure 3b displays the potential energy curve calculated with an increasing number of Gaussians per electron, this time with the VB correction (ρ = 1). In this case, as the number of Gaussians is increased, the potential energy curve tends toward the generalized valence bond (GVB) result [52], validating the implementation of the multiple Gaussian basis. The eFF results, also shown in Figure 3b, lie between the one and two Gaussian lines, clearly indicating that the scaling factors in eFF somehow address the impediment from the limited basis and lack of correlation. The potential energy curves for the two Gaussian case in Figure 3a, and one Gaussian case in Figure  3b, both exhibit a kink around 1.6 Bohr. We attribute this to the restricted basis attempting to represent two phenomenologically different configurations to the left and right of the kink; that is, the expectation of the two electrons' position being located between the two ions when close together and centered near the ions when separated. This effect is less apparent when a larger basis is used.  Figure 3 represents the static use of the code. To validate the dynamic component, We calculated the fundamental vibrational frequency of H 2 through the dynamical equations-of-motion. In the simulations, the two hydrogen nuclei we separated by a small distance of 0.01 Bohr and released, the time step was 0.603 as, and the total run time was 73 fs. The energy was conserved to within 10 −12 Hartree. The V ii energy was plotted versus time, creating a sinusoidal graph where the peak to peak time was averaged over ten oscillations to obtain the fundamental vibrational frequency. As we increased the number of Gaussians used, the results tended towards a consistent result. Without the VB correction, this was found to be 4563 cm −1 , almost 10 % higher than the experimental value of 4161 cm −1 [54]. In the simulation in which correlation was accounted for, the frequency was found to be 4204 cm −1 , just 1% higher than the experimental value. As before, a significant improvement was found when the number of Gaussians was increased from one to two.
We have demonstrated the applicability of the method to the hydrogen trial systems presented thus far. We now turn our attention to modeling a periodic hydrogen plasma system. For this work, we utilized systems consisting of 1024 ions and an equal number of electrons. At this time, the code is implemented serially in MATLAB and is unable to be run for long times across all permutations of approximations and basis sets. Thus, for the majority of results presented in Figure 5, we compare the pressure of an energy-minimized system at 0 K, calculated according to the virial theorem [55], to experimental results and other models calculated at 300 K. However, at these densities, the difference in pressure between 300 K and 0 K is negligible. This was confirmed by comparing, in Figure 5a, the data from Klakow et al. (solid blue line) and our kinetic energy exchange (dashed blue line), both calculated with a single Gaussian basis; they are almost identical over the range of densities considered.
In Figure 5a, which shows only results calculated with a single-Gaussian basis, we find that results calculated with our pairwise exchange term diverge from those of Jakob et al. which were calculated using exact exchange. However, with the addition of the VB correction, in this case with the scaling parameter ρ = 0.33, we are able to better match experimental results obtained on diamond anvil cells [56]. In this case, the VB correction is accounting for the pairwise exchange, the limited basis, and lack of correlation. We note that our result follows more closely the experimental data than eFF. In Figure 5b the effect of an increased basis is shown for the same dense hydrogen system, both with and without the VB correction. Interestingly, simply with the pairwise exchange energy, an increased Gaussian basis does not significantly improve the results. This suggests the system's limitation can be traced to either the lack of correlation or the pairwise exchange. However, with the VB correction, we can scale ρ to lower the pressure. With a two Gaussian basis and a scaled VB correction (ρ = 0.52), we can exactly match the experimental results across the density regime tested. These results suggest that improved correlation and full exchange play a larger role than the expanded basis. With that said, if we investigate the crystal structure of the two minimized systems with the correction, we find the single-Gaussian basis has an orthorhombic structure (c.f. Figure 5c). In contrast, the two-Gaussian system has the correct hexagonal-close-packed structure (c.f. Figure 5d). For our best model, the two-Gaussian basis with ρ = 0.52, we were able to run a periodic dynamical simulation using a 0.24 as timestep for 130 as. During this time, we scaled the ion velocities to have a kinetic energy consistent with 300 K and let the electrons equilibrate with the ions. During this simulation, the pressure was observed to oscillate between 125 GPa and 160 GPa, denoted by the cross and error bars in Figure 5b. This simulation demonstrates the model's applicability but highlights the importance of parallelizing the code for modeling larger systems.

Discussion
Atomistic simulations such as WPMD that go beyond the BO approximation may be necessary to describe the dynamics of dense plasmas systems and, in particular, non-equilibrium matter. However, in many implementations of WPMD, multiple approximations are used to achieve computational efficiency. This efficiency has led to the widespread use of the eFF flavor of WPMD, which, despite the restrictive single-Gaussian basis, has achieved remarkable success. This can, in part, be attributed to the semi-empirical scaling parameters that simultaneously attempt to correct for deficiencies in the basis, approximations to exchange, and lack of correlation. However, the use of scaling parameters makes the method's applicability in untested regimes, such as for higher-Z or extremely dense plasmas, questionable.
To reduce the number of scaling parameters, we have implemented a version of WPMD with an extended multiple Gaussian basis to describe the behavior of periodic systems of dense hydrogen. Furthermore, we use the improved Pauli potential suggested by H. Xiao, which includes the electron-electron and electron-ion components of exchange. However, in order to achieve computational scaling comparable to classical methods, one of the distinct advantages of the eFF method, we retain the use of a pairwise Pauli potential. We find improvements in the description of hydrogen atoms, molecules, and plasma systems as the number of Gaussians are increased from one to four, in agreement with previous work on non-periodic systems of a few electrons [36,42,43]. Notably, increasing the basis from one to two Gaussians provides the greatest improvement.
With the improved basis and exchange, the most significant remaining error is attributed to the model's lack of correlation. We implemented a simple VB correction based on the VB wave function, which we demonstrate in the hydrogen molecule and plasma. However, to match experimental results, we must utilize a single-scaling parameter on this correction term. Future work will focus on improving the correlation part of this term, which could be modified to take into account local order [45], or exploiting the robust correlation functionals within DFT [49]. In addition, we are working on parallelization of the code, which will allow for the description of larger plasma systems for greater times.