Next Article in Journal
Pulsed Discharge Plasma in High-Pressure Environment for Water Pollutant Degradation and Nanoparticle Synthesis
Next Article in Special Issue
Parametric Study of Proton Acceleration from Laser-Thin Foil Interaction
Previous Article in Journal
Multi-Device Piezoelectric Direct Discharge for Large Area Plasma Treatment
Previous Article in Special Issue
Method for Measuring the Laser Field and the Opacity of Spectral Lines in Plasmas
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Investigation into the Approximations Used in Wave Packet Molecular Dynamics for the Study of Warm Dense Matter

Department of Physics, University of Nevada, Reno, Reno, NV 89557, USA
*
Author to whom correspondence should be addressed.
Submission received: 31 March 2021 / Revised: 18 May 2021 / Accepted: 20 May 2021 / Published: 27 May 2021
(This article belongs to the Special Issue Laser–Plasma Interactions and Applications)

Abstract

:
Wave packet molecular dynamics (WPMD) has recently received a lot of attention as a computationally fast tool with which 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 investigated three of the main approximations ubiquitous to WPMD: a restricted basis set, approximations to exchange, and the lack of correlation. We examined each of these approximations in regard to atomic and molecular hydrogen in addition to a dense hydrogen plasma. We found 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 found 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.

1. 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 state [10,12]. The most prevalent of these techniques is density functional theory molecular dynamics (DFT-MD), in which the electrons are treated within the framework of either orbital-free [10] or Kohn–Sham density functional theory [13]. DFT-MD employs the Born–Oppenheimer (BO) approximation, in which the electrons are considered to respond instantaneously to the ion dynamics, usually justified by the disparate time-scales of the electron and ion motion. Although applicable for equilibrium properties, such as the equation of state [14], the BO treatment of the ions may not be suitable for calculating dynamic properties such as sound-speed and transport coefficients [6,11]. Furthermore, by its very nature, the BO approximation prohibits direct energy transfer between electrons and ions, and is therefore problematic for modeling non-equilibrium matter [15]. Ultimately, in the WDM regime, it is still unknown how the BO approximation impacts atomistic calculations.
Recently, several theoretical techniques have been developed that go beyond the BO approximation. The simplest technique couples the ions to a Langevin thermostat which models the electron–ion collisions through an additional stochastic Gaussian force added to the equations of motion [11,16,17,18,19]. While efficient, this phenomenological approach uses a single collision frequency that must be determined a priori. Other techniques that go beyond the BO approximation include using a thermally averaged, linearized Bohm potential, which has successfully described acoustic oscillations in warm dense aluminum [20,21], and the most widespread, wave packet molecular dynamics (WPMD) [5,6,22,23,24,25,26,27,28,29]. Time-dependent DFT, often used successfully in quantum chemistry [30], remains too computationally intensive to study the large systems of interest. In this work, we will focus on the applicability and approximations used within WPMD.
WPMD is a time-dependent quantum mechanical technique that simultaneously simulates, (1) the propagation of the ions as classical point particles, and (2) the electrons as quantum-mechanical entities. In WPMD, each electron is represented as a quantum wave packet, a spatially-localized complex function often implemented on a Gaussian basis [31]. A wave packet uniquely defines the state of a single electron, with the total many-body wave function being constructed from either a Hartree product or a Slater determinant [32,33], a choice that is driven by the importance of balancing exchange effects with computational cost. Equations of motion for the dynamical parameters are easily derived from variation of the time-dependent Schrodinger equation, where, for a single-Gaussian basis, they take on a simple Hamilton form [25,31]. The direct inclusion of electrons, and thus the effects of electron–ion interactions, means that WPMD intrinsically goes beyond the BO approximation. It is capable of computing electron–ion energy exchange in non-equilibrium systems, the effects of electron–ion collisions, and more generally calculating observables in quantum many-body systems [29,33].
Many flavors of WPMD exist that utilize varying degrees of approximation. The three most common approximations are a restricted basis, consisting of a single Gaussian per electron; a pairwise exchange interaction, often identified as a Pauli potential; and the exclusion of correlation. Furthermore, the Pauli potential itself is often assumed to depend only on the kinetic energy component of exchange in addition to the dependence on electron momentum being ignored [22,25,34]. With these simplifications, WPMD can obtain the same computational efficiency as in many classical methods [35,36]. These efficiencies have led to the widespread use of a semi-empirical WPMD method known as the electron force field (eFF). In eFF, several experimentally derived scaling parameters are used, with remarkable success, to correct for deficiencies in basis, approximations to exchange, and lacking correlations. To date, eFF has been used to investigate material properties in extreme environments [27,28,32,37,38], temperature relaxation rates in warm dense hydrogen [5], sound-speed in warm dense aluminum [5], and diffusion in warm dense hydrogen [6].
Despite the success of eFF, exactly how the scaling parameters address these approximations is not well understood. The predictive capability is limited, and its use in new regimes should always be corroborated with other models or experimental data [33]. For example, eFF was recently shown to underestimate the ion–ion correlation in dense aluminum plasma at temperatures of a few electronvolts [5,39]. While there has been some effort made to understand the effect of a pairwise 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 improvements in the ground state energy of the hydrogen atom with an increasing number of Gaussians in the electron basis; minimal improvements are observed beyond four Gaussians.
We 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 to 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 is 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 with 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.

2. Materials and Methods

2.1. 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 the calculation of the Slater determinant of all single electron wave functions, a computationally expensive operation [24,42,43]. To decrease the 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:
Ψ H ( X , t ) = k = 1 N e ψ k ( x k , t ) ,
where X = { x 1 , , x k , , x N e } , and x k are used to indicate the space occupied by the k-th electron. Using the Hartree product for the many-body electron wave function omits the correct anti-symmetry needed to ensure fermionic particles, such as electrons, do not violate the Pauli exclusion principle. To make practical use of the Hartree product, the use of Pauli potentials helps approximate the missing anti-symmetry. We will discuss the use of Pauli potentials in more detail in Section 2.3.
Another approximation implemented, for computational speed, in most WPMD techniques, is the utilization of a single isotropic Gaussian as a restricted wave function for each electron [22,23,24,32,34]. Gaussian bases have been used since the 1950s to calculate electronic structure [46]; however, to get the most accurate wave function requires the use of an increased basis, as was demonstrated in Figure 1. We will discuss the changes to the WPMD theory due to a larger basis set in more detail in Section 2.2.
The Gaussian wave function, used in WPMD, is typically paramaterized by a set of ten real physical variables q = { r , p , σ , p σ } , i.e.,
ψ k ( x k ) = 3 2 π σ 2 3 / 4 exp 3 4 σ 2 i p σ 2 σ r x k 2 + i p · r x k .
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,47],
N q ˙ = H q
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.
N q k q l = q k * q l ln Ψ H q * Ψ H ( q ) ,
where N q k q l is a matrix element of the norm matrix, q k represents a specific element of the set of time-dependent variational parameters, and q k * is the complex conjugate of q k . 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
H ^ = ( T ^ i + V ^ i i ) + ( T ^ e + V ^ e i + V ^ e e ) + H ^ H a r m ,
where the operators take on their usual classical,
T ^ i = I p I 2 2 M I , V ^ i i = 1 2 I , J Z I Z J | R I R J | ,
and quantum mechanical,
T ^ e = 1 2 k k 2 , V ^ e i = I k Z I | x k R I | , V ^ e e = 1 2 k , l 1 | x k x l | ,
definitions. T ^ e is the electron kinetic energy operator, V ^ e i is the electron–ion Coulombic potential energy operator, V ^ e e is the electron-electron Coulombic potential energy operator, T ^ i is the ion kinetic energy, and V ^ i i 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 (see [35,36]). This term is typically used to constrain the electron’s size, which may increase to the point that electron–ion interactions become negligible. As such, the harmonic energy term is not an energy term born out of physical arguments, but it comes from a computational necessity. For this potential we have used the form suggested by Zwicknagel et al. where H Harm = H ^ H a r m = k = 1 N e 9 8 γ 0 4 σ 2 and γ 0 is set to to half the simulation box length [48]. It should be noted that this term represents a small contribution to the total energy, as—in the WDM regime—most wave packets are constrained enough from the ionic potentials present [5,36].

2.2. Extension to Multiple Gaussians

The primary aim of this work was to investigate 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:
ψ k ( x k , t ) = n k 1 / 2 α = 1 N g φ k α ( x k , t ) ,
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, n k = α , γ O k α k γ = α , γ φ k α * φ k γ d 3 x is the normalization factor, and O k α k γ is the overlap of Gaussians. 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:
φ k α ( x k , t ) = d k α ( t ) e a k α ( t ) ( x k · x k ) + b k α ( t ) · x k + c k α ( t ) .
Here, the set of dynamical variables for each Gaussian is 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 not an independent variable and is used to ensure normalization of each Gaussian [43].
The norm matrix given in Equation (4) becomes a block diagonal in the Hartree Product ansatz, where the block for the k-th electron is given by
N q k α q k γ = 2 Im ( ψ k ) q k α | ( ψ k ) q k γ ( ψ k ) q k α | ψ k ψ k | ( ψ k ) q k γ ,
which within the framework of multiple Gaussians, can be expressed in terms of the derivatives of each individual Gaussian [43]:
N q k α q k γ = 2 Im n k 1 ( φ k α ( x ) ) q k α | ( φ k γ ( x ) ) q k γ L k , q k α * L k , q k γ
where L k , q k γ = ( n k ) 1 σ φ k σ | ( φ k γ ) q k γ , and ( φ k γ ) q k γ is the derivative with respect to one of the q k γ dynamic parameters.
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:
H H = T i + T e + V i i + V e i + V e e .
In each case, an analytical expression may be derived. For the semi-classical electron terms, they are most easily expressed as the products of overlapping an integral with a residual, i.e.,
T e = Ψ H | T ^ e | Ψ H = k α , γ O k α k γ T k α k γ e ,
V e i = Ψ H | V ^ e i | Ψ H = I k α , γ O k α k γ V I k α k γ e i ,
V e e = Ψ H | V ^ e e | Ψ H = k , l α , β , γ , δ O k α k γ O l β l δ V k α l β k γ l δ e e ,
where O k α l δ represents the overlap between the α Gaussian wave packet of the k-th electron and the δ Gaussian wave packet of the l-th electron. Expressions for the three residual terms T k α k γ e , V I k α k γ e i , V k α l β k γ l δ e e , and the overlap integral are easily derived in the Gaussian basis and are provided in reference [43].

2.3. 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:
H P = T e P + V e i P + V e e P
where
T e P = Ψ S | T ^ e | Ψ S Ψ H | T ^ e | Ψ H ,
V e i P = Ψ S | V ^ e i | Ψ S Ψ H | V ^ e i | Ψ H ,
V e e P = Ψ S | V ^ e e | Ψ S Ψ H | V ^ e e | Ψ H ,
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,
Ψ S ( x 1 , x 2 ) = 1 2 [ ψ 1 ( x 1 ) ψ 2 ( x 2 ) ψ 1 ( x 2 ) ψ 2 ( x 1 ) ] ,
analytic expressions for the two-particle Pauli energy given in Equations (17)–(19) were derived. The analytical expressions are as follows:
T e P = k , l α , β , γ , δ ( O k l O l k ) O k α k γ O l β l δ T k α k γ e + T l β l δ e O k α l δ O l β k γ T k α l δ e + T l β k γ e ( n k n l ) ( 1 ( O k l O l k ) ) ,
V e i P = I k , l α , β , γ , δ ( O k l O l k ) O k α k γ O l β l δ V I k α k γ e i O k α l δ O l β k γ V I k α l δ e i ( n k n l ) ( 1 O k l O l k ) ,
V e e P = k , l α , β , γ , δ ( O k l O l k ) O k α k γ O l β l δ V k α l β k γ l δ e e O k α l δ O l β k γ V k α l β l δ k γ e e ( n k n l ) ( 1 O k l O l k ) ,
where O k l = ( n k n l ) 1 / 2 α , δ φ k α | φ l δ , O l k = O k l * , and the residual terms, T k α l δ e , V I k α l δ e i , and V k α l β l δ k γ e e , are provided in reference [43].
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]. Figure 2a compares our pairwise Pauli potential with other published results. The results by Klakow et al. contain only the kinetic energy’s contribution to exchange, and for this system, agree with our T e P term. The eFF model also makes use of the kinetic energy’s Pauli potential but differs from Klakow et al.’s potential due to the incorporation of empirical scaling parameters [25]. Due to the lack of ions in the system, our total Pauli term contains just one additional term, V e e ( P ) . This term acts to lower 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 important roles 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 with 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 a decades-old problem in quantum mechanical many-body systems [49,50,51,52]. 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
H C = T e C + V e i C + V e e C ,
where
T e C = Ψ V B | T ^ e | Ψ V B Ψ H | T ^ e | Ψ H ,
V e i C = Ψ V B | V ^ e i | Ψ V B Ψ H | V ^ e i | Ψ H ,
V e e C = Ψ V B | V ^ e e | Ψ V B Ψ H | V ^ e e | Ψ H .
As before, utilizing the usual definition of the two-particle VB wave function,
Ψ V B ( x 1 , x 2 ) = 1 2 [ ψ 1 ( x 1 ) ψ 2 ( x 2 ) + ψ 1 ( x 2 ) ψ 2 ( x 1 ) ] ,
we derived analytic expressions for the pairwise VB correction energy given in Equations (25)–(27). The analytical expressions are as follows:
T e C = k , l α , β , γ , δ O k α l δ O l β k γ T k α l δ e + T l β k γ e ( O k l O l k ) O k α k γ O l β l δ T k α k γ e + T l β l δ e ( n k n l ) ( 1 + ( O k l O l k ) ) ,
V e i C = I k , l α , β , γ , δ O k α l δ O l β k γ V I k α l δ e i ( O k l O l k ) O k α k γ O l β l δ V I k α k γ e i ( n k n l ) ( 1 + O k l O l k ) ,
V e e C = k , l α , β , γ , δ O k α l δ O l β k γ V k α l β l δ k γ e e ( O k l O l k ) O k α k γ O l β l δ V k α l β k γ l δ e e ( n k n l ) ( 1 + O k l O l k ) .
The VB correction was implemented between pairs of electrons with the opposite spin. The same spin electrons were assumed to be well separated spatially, due to the Pauli potential, and to not contribute significantly to this term [53].
The total energy of our system, H = H H + H P + H C may be written as follows:
H = ( T i + V i i ) + ( T e + V e i + V e e ) + ( T e P + V e i P + V e e P ) δ + ρ ( T e C + V e i C + V e e C ) δ ,
where δ equals unity when the spins of the electrons are parallel and zero otherwise; δ is 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 of 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.

2.4. The Implementation of Periodic Boundary Conditions

The simulation of dense plasmas necessitates the use of a finite-in-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 c u t , defined so that the first, second, and third derivatives at the endpoints are zero [32,54]. i.e.,
f cutoff = 20 x 7 70 x 6 + 84 x 5 35 x 4 + 1 ,
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. There are no restrictions on the positions of the Gaussians within one electron; however, when applying 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.

3. Results

Dense hydrogen is one of the simplest physical systems, and for this reason, it has become the prototypical test-bed for atomistic models of dense plasmas. At this time, our code is implemented serially in MATLAB and is unable to be run for long periods of time for all permutations of approximations, basis sets, and parameter space. Thus, we limit our investigation to cold dense hydrogen, which is often used to benchmark WPMD approximations before being extended to higher temperature conditions [23,24,25].
Figure 1 demonstrates that WPMD approaches the analytic solution of the hydrogen atom with an increasing Gaussian basis. The ground state energy can be obtained to within one percent, utilizing a four-Gaussian basis. Beyond this number, we observed only minor changes in the ground state energy and the electron wave function shape. It should also be noted that the most significant improvement occurs 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 an 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, only negligible improvements are observable 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 [53], 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’ positions 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 results shown in Figure 4, the two hydrogen nuclei were 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 i i 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 [55]. 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. 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 [56], to experimental results and other models calculated at 300 K. However, at these densities, the difference in pressure between 300 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 can see that the results calculated with our pairwise exchange term diverge from those of Jakob et al. which were calculated using an exact exchange. However, with the addition of the VB correction, in this case with the scaling parameter ρ = 0.33 , we were able to better match experimental results obtained on diamond anvil cells [57]. In this case, the VB correction is accounting for the pairwise exchange, the limited basis, and the 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 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.

4. 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. Inclusion of the electron–electron exchange in the Pauli potential, and correlation approximation, allow for spatial correlations among dynamically moving and interacting electrons. 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 achieved improvements in the descriptions of hydrogen atoms, molecules, and plasma systems as the number of Gaussians was 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 provided the greatest improvement.
With the improved basis and exchange, the most significant remaining error was attributed to the model’s lack of correlation. We implemented a simple VB correction based on the VB wave function, which we demonstrated 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 [50]. In addition, this work has been restricted to investigating low-temperature dense systems due to the serial implementation in MATLAB. To overcome computational limitations and explore higher temperature regions of phase space, we plan to port the code over to C++ with CPU scalability in mind. However, the success of simpler WPMD models in the WD [5,14,32] leads us to believe our improvements will extend the region of validity; however, that is yet to be seen.

Author Contributions

W.A.A. and T.G.W. contributed equally to this manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This material is partially based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under grant number DE-SC0019268.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

Special thanks to Nuvraj Bilkhu, Jacob Molina, and Cameron Allen.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Guillot, T. A Comparison of the Interiors of Jupiter and Saturn. Planet. Space Sci. 1999, 47, 1183–1200. [Google Scholar] [CrossRef] [Green Version]
  2. Paquette, C.; Pelletier, C.; Fontaine, G.; Michaud, G. Diffusion Coefficients for Stellar Plasmas. ApJS 1986, 61, 177–195. [Google Scholar] [CrossRef]
  3. Atzeni, S.; Meyer-ter-Vehn, J. Thermonuclear Fusion and Confinement. In The Physics of Inertial Fusion: Beam Plasma Interaction, Hydrodynamics, Hot Dense Matter, 1st ed.; Clarendon Press: Oxford, UK, 2004; Volume 125, pp. 31–36. [Google Scholar]
  4. Ichimaru, S. Strongly Coupled Plasmas: High-Density Classical Plasmas and Degenerate Electron Liquids. Rev. Mod. Phys. 1982, 54, 1017–1059. [Google Scholar] [CrossRef]
  5. Davis, R.A.; Angermeier, W.A.; Hermsmeier, R.K.T.; White, T.G. Ion Modes in Dense Ionized Plasmas through Nonadiabatic Molecular Dynamics. Phys. Rev. Res. 2020, 2, 043139. [Google Scholar] [CrossRef]
  6. Yao, Y.; Zeng, Q.; Chen, K.; Kang, D.; Hou, Y.; Ma, Q.; Dai, J. Reduced Ionic Diffusion by the Dynamic Electron–Ion Collisions in Warm Dense Hydrogen. Phys. Plasmas 2021, 28, 012704. [Google Scholar] [CrossRef]
  7. Lorenzen, W.; Becker, A.; Redmer, R. Progress in Warm Dense Matter and Planetary Physics. In Frontiers and Challenges in Warm Dense Matter; Graziani, F., Desjarlais, M.P., Redmer, R., Trickey, S.B., Eds.; Lecture Notes in Computational Science and Engineering; Springer International Publishing: Cham, Switzerland, 2014; pp. 203–234. [Google Scholar]
  8. Grabowski, P.E.; Hansen, S.B.; Murillo, M.S.; Stanton, L.G.; Graziani, F.R.; Zylstra, A.B.; Baalrud, S.D.; Arnault, P.; Baczewski, A.D.; Benedict, L.X.; et al. Review of the First Charged-Particle Transport Coefficient Comparison Workshop. High Energy Density Phys. 2020, 37, 100905. [Google Scholar] [CrossRef]
  9. Salin, G.; Caillol, J.M. Equilibrium Molecular Dynamics Simulations of the Transport Coefficients of the Yukawa One Component Plasma. Phys. Plasmas 2003, 10, 1220–1230. [Google Scholar] [CrossRef]
  10. White, T.G.; Richardson, S.; Crowley, B.J.B.; Pattison, L.K.; Harris, J.W.O.; Gregori, G. Orbital-Free Density-Functional Theory Simulations of the Dynamic Structure Factor of Warm Dense Aluminum. Phys. Rev. Lett. 2013, 111, 175002. [Google Scholar] [CrossRef] [Green Version]
  11. Mabey, P.; Richardson, S.; White, T.G.; Fletcher, L.B.; Glenzer, S.H.; Hartley, N.J.; Vorberger, J.; Gericke, D.O.; Gregori, G. A Strong Diffusive Ion Mode in Dense Ionized Matter Predicted by Langevin Dynamics. Nat. Commun. 2017, 8, 14125. [Google Scholar] [CrossRef] [Green Version]
  12. Graziani, F.R.; Batista, V.S.; Benedict, L.X.; Castor, J.I.; Chen, H.; Chen, S.N.; Fichtl, C.A.; Glosli, J.N.; Grabowski, P.E.; Graf, A.T.; et al. Large-Scale Molecular Dynamics Simulations of Dense Plasmas: The Cimarron Project. High Energy Density Phys. 2012, 8, 105–131. [Google Scholar] [CrossRef]
  13. Rüter, H.R.; Redmer, R. Ab Initio Simulations for the Ion-Ion Structure Factor of Warm Dense Aluminum. Phys. Rev. Lett. 2014, 112, 145007. [Google Scholar] [CrossRef]
  14. Kang, D.; Hou, Y.; Zeng, Q.; Dai, J. Unified First-Principles Equations of State of Deuterium-Tritium Mixtures in the Global Inertial Confinement Fusion Region. Matter Radiat. Extrem. 2020, 5, 055401. [Google Scholar] [CrossRef]
  15. Clérouin, J.; Robert, G.; Arnault, P.; Ticknor, C.; Kress, J.D.; Collins, L.A. Evidence for Out-of-Equilibrium States in Warm Dense Matter Probed by x-Ray Thomson Scattering. Phys. Rev. E 2015, 91, 011101. [Google Scholar] [CrossRef] [Green Version]
  16. Dai, J.; Yuan, J. Large-Scale Efficient Langevin Dynamics, and Why It Works. EPL 2009, 88, 20001. [Google Scholar] [CrossRef]
  17. Tamm, A.; Caro, M.; Caro, A.; Samolyuk, G.; Klintenberg, M.; Correa, A.A. Langevin Dynamics with Spatial Correlations as a Model for Electron-Phonon Coupling. Phys. Rev. Lett. 2018, 120, 185501. [Google Scholar] [CrossRef] [Green Version]
  18. Duffy, D.M.; Rutherford, A.M. Including the Effects of Electronic Stopping and Electron–Ion Interactions in Radiation Damage Simulations. J. Phys. Condens. Matter 2007, 19, 016207. [Google Scholar] [CrossRef]
  19. Rutherford, A.M.; Duffy, D.M. The Effect of Electron–Ion Interactions on Radiation Damage Simulations. J. Phys. Condens. Matter 2007, 19, 496201. [Google Scholar] [CrossRef]
  20. Larder, B.; Gericke, D.O.; Richardson, S.; Mabey, P.; White, T.G.; Gregori, G. Fast Nonadiabatic Dynamics of Many-Body Quantum Systems. Sci. Adv. 2019, 5, eaaw1634. [Google Scholar] [CrossRef] [Green Version]
  21. Moldabekov, Z.A.; Dornheim, T.; Gregori, G.; Graziani, F.; Bonitz, M.; Cangi, A. The Quantum Bohm Potential for Many-Fermion Systems. arXiv 2021, arXiv:2103.08523. [Google Scholar]
  22. Klakow, D.; Toepffer, C. Hydrogen under Extreme Conditions. Phys. Lett. A 1994, 192, 55–59. [Google Scholar] [CrossRef]
  23. Knaup, M.; Zwicknagel, G.; Reinhard, P.G.; Toepffer, C. Wave Packet Molecular Dynamics Simulations of Hydrogen under Extreme Conditions. Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers Detect. Assoc. Equip. 2001, 464, 267–270. [Google Scholar] [CrossRef]
  24. Jakob, B.; Reinhard, P.G.; Toepffer, C.; Zwicknagel, G. Wave Packet Simulation of Dense Hydrogen. Phys. Rev. E 2007, 76, 036406. [Google Scholar] [CrossRef]
  25. Su, J.T.; Goddard, W.A. Excited Electron Dynamics Modeling of Warm Dense Matter. Phys. Rev. Lett. 2007, 99, 185003. [Google Scholar] [CrossRef] [Green Version]
  26. Lavrinenko, Y.S.; Morozov, I.V.; Valuev, I.A. Thermodynamic Properties of the Nonideal Hydrogen Plasmas: Comparison of Different Simulation Techniques. J. Phys. Conf. Ser. 2018, 946, 012097. [Google Scholar] [CrossRef]
  27. Jaramillo-Botero, A.; Su, J.; Qi, A.; Goddard, W.A. Large-Scale, Long-Term Nonadiabatic Electron Molecular Dynamics for Describing Material Properties and Phenomena in Extreme Environments. J. Comput. Chem. 2011, 32, 497–512. [Google Scholar] [CrossRef]
  28. Kim, H.; Su, J.T.; Goddard, W.A. High-Temperature High-Pressure Phases of Lithium from Electron Force Field (eFF) Quantum Electron Dynamics Simulations. Proc. Natl. Acad. Sci. USA 2011, 108, 15101–15105. [Google Scholar] [CrossRef] [Green Version]
  29. Ma, Q.; Dai, J.; Kang, D.; Murillo, M.S.; Hou, Y.; Zhao, Z.; Yuan, J. Extremely Low Electron-Ion Temperature Relaxation Rates in Warm Dense Hydrogen: Interplay between Quantum Electrons and Coupled Ions. Phys. Rev. Lett. 2019, 122, 015001. [Google Scholar] [CrossRef] [Green Version]
  30. Li, X.; Tully, J.C.; Schlegel, H.B.; Frisch, M.J. Ab Initio Ehrenfest Dynamics. J. Chem. Phys. 2005, 123, 084106. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  31. Feldmeier, H.; Schnack, J. Molecular Dynamics for Fermions. Rev. Mod. Phys. 2000, 72, 655–688. [Google Scholar] [CrossRef] [Green Version]
  32. Su, J.T.; Goddard, W.A. The Dynamics of Highly Excited Electronic Systems: Applications of the Electron Force Field. J. Chem. Phys. 2009, 131, 244501. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Grabowski, P.E. A Review of Wave Packet Molecular Dynamics. In Frontiers and Challenges in Warm Dense Matter; Graziani, F., Desjarlais, M.P., Redmer, R., Trickey, S.B., Eds.; Lecture Notes in Computational Science and Engineering; Springer International Publishing: Cham, Switzerland, 2014; pp. 265–282. [Google Scholar]
  34. Boal, D.H.; Glosli, J.N. Quasiparticle Model for Nuclear Dynamics Studies: Ground-State Properties. Phys. Rev. C 1988, 38, 1870–1878. [Google Scholar] [CrossRef] [PubMed]
  35. Morozov, I.V.; Valuev, I.A. Localization Constraints in Gaussian Wave Packet Molecular Dynamics of Nonideal Plasmas. J. Phys. A Math. Theory 2009, 42, 214044. [Google Scholar] [CrossRef]
  36. Grabowski, P.E.; Markmann, A.; Morozov, I.V.; Valuev, I.A.; Fichtl, C.A.; Richards, D.F.; Batista, V.S.; Graziani, F.R.; Murillo, M.S. Wave Packet Spreading and Localization in Electron-Nuclear Scattering. Phys. Rev. E 2013, 87, 063104. [Google Scholar] [CrossRef] [Green Version]
  37. Xiao, H.; Jaramillo-Botero, A.; Theofanis, P.L.; Goddard, W.A. Non-Adiabatic Dynamics Modeling Framework for Materials in Extreme Conditions. Mech. Mater. 2015, 90, 243–252. [Google Scholar] [CrossRef] [Green Version]
  38. Lan, M.; Yang, Z.H.; Wang, X. Displacement Damage in Silicon Studied by the Electronic Force Field Method in the keV Regime. Comput. Mater. Sci. 2020, 179, 109697. [Google Scholar] [CrossRef]
  39. Fletcher, L.B.; Lee, H.J.; Döppner, T.; Galtier, E.; Nagler, B.; Heimann, P.; Fortmann, C.; Lepape, S.; Ma, T.; Millot, M.; et al. Ultrabright X-Ray Laser Scattering for Dynamic Warm Dense Matter Physics. Nat. Photonics 2015, 9, 274–279. [Google Scholar] [CrossRef]
  40. Knaup, M.; Reinhard, P.G.; Toepffer, C.; Zwicknagel, G. Wave Packet Molecular Dynamics Simulations of Hydrogen at Mbar Pressures. Comput. Phys. Commun. 2002, 147, 202–204. [Google Scholar] [CrossRef]
  41. Knaup, M.; Reinhard, P.G.; Toepffer, C.; Zwicknagel, G. Wave Packet Molecular Dynamics Simulations of Warm Dense Hydrogen. J. Phys. A Math. Gen. 2003, 36, 6165–6171. [Google Scholar] [CrossRef]
  42. Morozov, I.V.; Valuev, I.A. Improvement of Wave Packet Molecular Dynamics Using Packet Splitting. Contrib. Plasma Phys. 2012, 52, 140–144. [Google Scholar] [CrossRef]
  43. Valuev, I.A.; Morozov, I.V. Extension of the Wave Packet Molecular Dynamics Method towards the Accurate Quantum Simulations of Electron Dynamics. J. Phys. Conf. Ser. 2015, 653, 012153. [Google Scholar] [CrossRef] [Green Version]
  44. Cohen-Tannoudji, C.; Diu, B.; Laloë, F. Quantum Mechanics, Volume 1: Basic Concepts, Tools, and Applications, 1st ed.; Wiley: New York, NY, USA, 1977; Volume 1. [Google Scholar]
  45. Xiao, H. First Principles Based Multiparadigm Modeling of Electronic Structures and Dynamics. Ph.D. Thesis, California Institute of Technology, Pasadena, CA, USA, 2015. [Google Scholar] [CrossRef]
  46. Boys, S.F.; Egerton, A.C. Electronic Wave Functions—I. A General Method of Calculation for the Stationary States of Any Molecular System. Proc. R. Soc. London. Ser. A Math. Phys. Sci. 1950, 200, 542–554. [Google Scholar] [CrossRef]
  47. Littlejohn, R.G. The Semiclassical Evolution of Wave Packets. Phys. Rep. 1986, 138, 193–291. [Google Scholar] [CrossRef] [Green Version]
  48. Zwicknagel, G.; Pschiwul, T. WPMD Simulations of a Two-Component Plasma. J. Phys. A Math. Gen. 2006, 39, 4359–4364. [Google Scholar] [CrossRef]
  49. Löwdin, P.O. Quantum Theory of Many-Particle Systems. III. Extension of the Hartree-Fock Scheme to Include Degenerate Systems and Correlation Effects. Phys. Rev. 1955, 97, 1509–1520. [Google Scholar] [CrossRef]
  50. Lavrinenko, Y.S.; Morozov, I.V.; Valuev, I.A. Wave Packet Molecular Dynamics–Density Functional Theory Method for Non-Ideal Plasma and Warm Dense Matter Simulations. Contrib. Plasma Phys. 2019, 59, e201800179. [Google Scholar] [CrossRef]
  51. Chachiyo, T.; Chachiyo, H. Understanding Electron Correlation Energy through Density Functional Theory. Comput. Theor. Chem. 2020, 1172, 112669. [Google Scholar] [CrossRef] [Green Version]
  52. Lavrinenko, Y.S.; Morozov, I.V.; Valuev, I.A. High Performance Wave Packet Molecular Dynamics with Density Functional Exchange-Correlation Term for Non-Ideal Plasma Simulations. J. Phys. Conf. Ser. 2021, 1787, 012043. [Google Scholar] [CrossRef]
  53. Su, J.T.l. An Electron Force Field for Simulating Large Scale Excited Electron Dynamics. Ph.D. Thesis, California Institute of Technology, Pasadena, CA, USA, 2007. [Google Scholar] [CrossRef]
  54. Izvekov, S.; Swanson, J.M.J.; Voth, G.A. Coarse-Graining in Interaction Space: A Systematic Approach for Replacing Long-Range Electrostatics with Short-Range Potentials. J. Phys. Chem. B 2008, 112, 4711–4724. [Google Scholar] [CrossRef]
  55. Dickenson, G.D.; Niu, M.L.; Salumbides, E.J.; Komasa, J.; Eikema, K.S.E.; Pachucki, K.; Ubachs, W. Fundamental Vibration of Molecular Hydrogen. Phys. Rev. Lett. 2013, 110, 193601. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Landau, L.D.; Lifshitz, E.M. Statistical Physics: Volume 5, 3rd ed.; Butterworth-Heinemann: Oxford, UK, 2013. [Google Scholar]
  57. Loubeyre, P.; LeToullec, R.; Hausermann, D.; Hanfland, M.; Hemley, R.J.; Mao, H.K.; Finger, L.W. X-Ray Diffraction and Equation of State of Hydrogen at Megabar Pressures. Nature 1996, 383, 702–704. [Google Scholar] [CrossRef]
Figure 1. The ground state energy of the hydrogen atom. The crosses were the calculated from minimization of the energy equations as the number of Gaussian basis functions were increased. The dotted line shows the exact ground state energy of hydrogen [44]. The insert shows the changes to the shape of the wave functions with an increasing number of basis functions. Beyond four Gaussians, improvements in the ground state energy and shape of the wave function are negligible.
Figure 1. The ground state energy of the hydrogen atom. The crosses were the calculated from minimization of the energy equations as the number of Gaussian basis functions were increased. The dotted line shows the exact ground state energy of hydrogen [44]. The insert shows the changes to the shape of the wave functions with an increasing number of basis functions. Beyond four Gaussians, improvements in the ground state energy and shape of the wave function are negligible.
Plasma 04 00020 g001
Figure 2. Validation of the three-term Pauli potential. (a) The Pauli potential consisting of all three energy terms is compared to the potential used in eFF and that developed by Klakow et al. which uses only the kinetic energy term [32]. Results are given for the exchange energy between two electrons of fixed-width σ = 3 / 4 . Unlike the electrostatic energy, shown by the dashed line, the Pauli potentials do not exhibit long-ranged behavior. (b) Verification of the Pauli potential following the method suggested by Xiao et al. [45]. The different energy contributions upon anti-symmetrization for the triplet H 2 system. In this system, the unrestricted Hartree–Fock (UHF) result is exact. Each electron is represented by a single Gaussian of radius 1.7 Bohr centered on the atom. Note, the eFF potential does not match the exact UHF result.
Figure 2. Validation of the three-term Pauli potential. (a) The Pauli potential consisting of all three energy terms is compared to the potential used in eFF and that developed by Klakow et al. which uses only the kinetic energy term [32]. Results are given for the exchange energy between two electrons of fixed-width σ = 3 / 4 . Unlike the electrostatic energy, shown by the dashed line, the Pauli potentials do not exhibit long-ranged behavior. (b) Verification of the Pauli potential following the method suggested by Xiao et al. [45]. The different energy contributions upon anti-symmetrization for the triplet H 2 system. In this system, the unrestricted Hartree–Fock (UHF) result is exact. Each electron is represented by a single Gaussian of radius 1.7 Bohr centered on the atom. Note, the eFF potential does not match the exact UHF result.
Plasma 04 00020 g002
Figure 3. Binding energy of the H 2 molecule. (a) VB correction is excluded. Circles represent the unrestricted Hartree–Fock (UHF) solution and crosses represent the exact solution. (b) VB correction ( ρ = 1 ) is included. Circles represent the generalized valence bond (GVB) solution and crosses denote the exact solution. The thick solid line represents the solution obtained using the eFF method.
Figure 3. Binding energy of the H 2 molecule. (a) VB correction is excluded. Circles represent the unrestricted Hartree–Fock (UHF) solution and crosses represent the exact solution. (b) VB correction ( ρ = 1 ) is included. Circles represent the generalized valence bond (GVB) solution and crosses denote the exact solution. The thick solid line represents the solution obtained using the eFF method.
Plasma 04 00020 g003
Figure 4. Vibrational frequency of the H 2 molecule with an increasing Gaussian basis set. Blue stars are the frequencies with no VB correction. Red crosses are the frequencies calculated with the VB correction ( ρ = 1 ). The black dotted line represents the experimental value of the fundamental frequency of H 2 . In each case the fundamental frequency was calculated by displacing two H ions by 0.01 Bohr and observing the resulting oscillations.
Figure 4. Vibrational frequency of the H 2 molecule with an increasing Gaussian basis set. Blue stars are the frequencies with no VB correction. Red crosses are the frequencies calculated with the VB correction ( ρ = 1 ). The black dotted line represents the experimental value of the fundamental frequency of H 2 . In each case the fundamental frequency was calculated by displacing two H ions by 0.01 Bohr and observing the resulting oscillations.
Plasma 04 00020 g004
Figure 5. The pressure of dense hydrogen. (a) A comparison of our results to those in the literature for a single-Gaussian basis. Our results were calculated at 0 K, whereas the published results, both experimental and theoretical, were calculated at 300 K. The small offset between the Klakow et al. results and our results using purely the kinetic energy component of exchange demonstrate the applicability of this approximation. (b) Extension to multiple Gaussians. The red lines, calculated with pairwise exchange only, show that increasing the number of Gaussians has only a small effect on the pressure curve. The VB correction is needed to best match experimental results. (c) Orthorhombic crystal structure of hydrogen at a Wigner–Seitz radius of 1.6 Bohr for a single-Gaussian basis and scaled VB correction ( ρ = 0.33 ). (d) Hexagonal-close-packed crystal structure of hydrogen at a Wigner–Seitz radius of 1.6 with a two-Gaussian basis and scaled VB correction ( ρ = 0.52 ).
Figure 5. The pressure of dense hydrogen. (a) A comparison of our results to those in the literature for a single-Gaussian basis. Our results were calculated at 0 K, whereas the published results, both experimental and theoretical, were calculated at 300 K. The small offset between the Klakow et al. results and our results using purely the kinetic energy component of exchange demonstrate the applicability of this approximation. (b) Extension to multiple Gaussians. The red lines, calculated with pairwise exchange only, show that increasing the number of Gaussians has only a small effect on the pressure curve. The VB correction is needed to best match experimental results. (c) Orthorhombic crystal structure of hydrogen at a Wigner–Seitz radius of 1.6 Bohr for a single-Gaussian basis and scaled VB correction ( ρ = 0.33 ). (d) Hexagonal-close-packed crystal structure of hydrogen at a Wigner–Seitz radius of 1.6 with a two-Gaussian basis and scaled VB correction ( ρ = 0.52 ).
Plasma 04 00020 g005
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Angermeier, W.A.; White, T.G. An Investigation into the Approximations Used in Wave Packet Molecular Dynamics for the Study of Warm Dense Matter. Plasma 2021, 4, 294-308. https://0-doi-org.brum.beds.ac.uk/10.3390/plasma4020020

AMA Style

Angermeier WA, White TG. An Investigation into the Approximations Used in Wave Packet Molecular Dynamics for the Study of Warm Dense Matter. Plasma. 2021; 4(2):294-308. https://0-doi-org.brum.beds.ac.uk/10.3390/plasma4020020

Chicago/Turabian Style

Angermeier, William A., and Thomas G. White. 2021. "An Investigation into the Approximations Used in Wave Packet Molecular Dynamics for the Study of Warm Dense Matter" Plasma 4, no. 2: 294-308. https://0-doi-org.brum.beds.ac.uk/10.3390/plasma4020020

Article Metrics

Back to TopTop