Reconstruction of release isentropes based on first-principle simulations

In this work we compare different methods of reconstruction of isentropic expansion curves based on quantum molecular dynamics simulation data. The VASP code is used. We analyze accuracy and computational complexity of three methods: the Zel’dovich approach, a method of re-shock Hugoniots and a direct calculation of entropy using the 2PT model. In the first method, an ordinary differential equation for temperature is solved numerically. The second method is based on the feature of the second-order contact of a Hugoniot and isentrope at the same initial point. The density of vibrational states which is reconstructed from a velocity autocorrelation function is used to calculate entropy by the third method. For aluminum the three methods turned out to give very close results for a release isentrope in agreement with experimental data. For molybdenum an experimental release isentrope is well reproduced by the second method.


Introduction
One of the most important problems in modern physics is the study of thermodynamic properties of materials under extreme conditions. This is necessary both for the construction of wide-range equations of state (EOSs) of materials at high temperatures and pressures, and for solving of a lot of applied and fundamental tasks in power engineering, space research, military industry, astrophysics, and solid state physics. EOSs of matter are required to complete the set of hydrodynamic equations, and, in particular, in numerical simulations of many processes and phenomena of high-energy density physics (inertial confinement fusion, hypervelocity impact, and many others). Dynamic shock-wave experiments are the main source of information about the EOS at high pressure and temperature. They also provide a basis for the verification of numerous theoretical models of condensed matter and plasma.
As experiments on shock compression and isentropic expansion are complicated and very expensive, a method of estimation of measured parameters is required. However, self-similar analysis of wave processes is based upon a generally unknown EOS. On the other hand, theoretical description of shock-wave experiments in condensed matter remains a serious problem since the appearance of first publications more than 60 years ago [1][2][3].
The development of computational methods and supercomputers [4][5][6] made it possible to obtain reliable data on thermodynamic properties of substances as a result of ab initio calculations. During the last two decades the quantum molecular dynamics (QMD) method 2 1234567890 ''"" became very popular. It is based upon the density functional theory (DFT), but also takes into account the movement of ions. Currently more than 1000 atoms can be used in calculations [7][8][9][10], this allows one to study unordered systems and even phase transitions [7,10]. This makes QMD a promising method for the description of thermophysical properties of matter in a wide range of parameters. Indeed, in our previous works we demonstrated high accuracy of QMD computations for Al, Cu, and LiD [11][12][13]. In this work we apply QMD to reproduce release isentropes of Al and Mo using three different methods.
2. Methods of reconstruction of isentropic expansion curves 2.1. Zel'dovich's approach (ZA) To calculate release isentropes, we use the Zel'dovich's approach [14] and integrate an ordinary differential equation for temperature T : Here, P -the pressure; V -the specific volume; E-the specific internal energy; S-the specific entropy. Formally speaking the ZA is an exact method of isentrope reconstruction: the accuracy is determined only by the right-hand side of equation (1). In the present work, we use a global reconstruction of the derivative (∂E/∂P ) V . For this purpose we compute the interpolation functions E(T, V ) and P (T, V ) from the VASP calculations on the grid of isotherms and isochores using a bicubic spline. In this case one can use an integration scheme of any order with any step on V for the solution of the ordinary differential equation (1); also the accuracy of calculation of the partial derivative (∂E/∂P ) V may be improved by the refinement of the grid of isotherms and isochors.
To determine the particle velocity (u p ) along the isentrope we use Riemann invariants, according to which [15] where ρ-the density.

Method of re-shock Hugoniots (RH)
This method is based on the feature of the second-order contact of a Hugoniot and isentrope at the same initial point [15]. We reconstruct a release isentrope by a series of re-shock Hugoniots. Each successive point on the calculated curve is related to the previous one by the Hugoniot equation (3), where the previous state parameters are denoted by index 0. To obtain a point on the re-shock Hugoniot, we make a series of calculations along an isotherm and then we use the bisection method to get parameters for which the function H(P, V, E) equals zero with a desired accuracy, As a result of our calculations we get pressure and internal energy of matter on the Hugoniot (3) at a given density. Shock (u s ) and particle (u p ) velocity can be easily restored using the laws of conservation of mass and momentum: where u p0 is the particle velocity in the previous state [16].

Method of direct isentrope reconstruction (DR)
The 2PT model is based upon the velocity autocorrelation function (VACF) and its frequency Fourier transform known as a vibrational density of states (VDOS). The VDOS function F (ν) can be obtained by the Fourier transform of the VACF Z(t): with Here ν is frequency, v(t) is the velocity of a particle at moment t, a dot sign means the scalar product, angular brackets denote the averaging over the trajectories of all particles. As can be seen from (6) the definition of Z(t) contains a normalization factor such that Z(0) ≡ 1. The integration of the VDOS function over positive frequencies gives the total number of degrees of freedom of the system. Thus F (ν) as defined here has the following normalization Following Lin et al [17] the total spectrum is decomposed into the solid-like F s (ν) and gas-like where f g is the gas-like fraction of the system. The total ionic entropy is decomposed into two terms as well and each of them in turn can be written as integrals over the components with appropriate weighting functions W s and W g , where k B is Boltzmann's constant, and weighting function for the solid-like component W s is given by a quantum-corrected one-phase thermodynamic model The gas-like component can be described as the system of hard spheres (HS) [17] or using a memory function (MF) representation [18].
Total entropy includes also the electronic contribution to the entropy where f i is the Fermi occupation of the ith electronic state obtained from the FT-DFT electronic minimization at at each ionic time step.

Simulation parameters
For our simulations we use the QMD approach also known as ab initio molecular dynamics (AIMD), which is based upon the Born-Oppenheimer approximation in which the electrons are immediately adjusted to the current spatial arrangement of the ions. For QMD simulation we use the plane-wave pseudopotential code VASP [20][21][22][23] and a projector augmented-wave (PAW) [24] potential. The generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) [25,26] corrections for the exchange-correlation functional is used. The Fermi smearing for electron occupancies is applied.
The density of matter ρ is defined by the lattice parameter and remains constant during the simulation process. The dynamics of Al and Mo atoms was simulated within no less than 1 ps with 2 fs time step. We use the Nosé-Hoover thermostat to maintain a given temperature. Then, by the Hellmann-Feynman theorem, the forces acting on the ions were computed and from the Newton's equations new coordinates and velocities of the ions were obtained. One point (Γ-point) in the Brillouin zone was used to determine the band structure. Thermodynamic parameters of the system were found by averaging of the corresponding values at the equilibrium stage of the simulation, the estimated statistical error was less than 1%. The kinetic contribution of ions to pressure was added to the QMD pressure obtained by VASP.
In this work, the QMD simulation was performed for 108 atoms of Al and 128 atoms of Mo, initially arranged in the ideal fcc and bcc lattices, respectively, with periodic boundary conditions. For Al simulations were performed with a 3-electron pseudopotential in the range of temperatures 3-10 kK and densitites 1.5-5 g/cm 3 . For Mo simulations were carried out with a 6-electron pseudopotential in the range of temperatures 3-5 kK and densitites 9.6-15 g/cm 3 .

Aluminum
In figure 1 presented are the principle Hugoniot of Al together with one isentrope which corresponds to experiment [19]. The QMD Hugoniot of Al was obtained in our previous work [11] from the numerical solution of equation (3). The initial point of the isentrope corresponded to a shocked state with an experimental mass velocity u pH = 6.44 km/s; other parameters obtained by QMD were as follows: density ρ H = 5 g/cm 3 , temperature T H ≈ 10 kK, and pressure P H ≈ 2500 kbar. The release isentrope was restored using the three methods described in section 2. As the ZA is an exact method among the two others, we consider the ZA isentrope as a standard. In the RH method three different step sizes on density were used: 2 g/cm 3 , 1 g/cm3, and 0.5 g/cm 3 . A smaller step should provide a more accurate result. It can be seen from figure 1(b) that the variation of pressure between the points obtained with the three steps on density is noticeable but small (see the signs at ρ ≈ 3 g/cm 3 ); the point corresponding to ∆ρ = 2 g/cm 3 shows the biggest deviation to lower pressures from the ZA isentrope. However, in the pressure-particle velocity plane with the linear scale on pressure (figure 2) the RH points are situated on both sides of the isentrope, see figure 2(b). The DR method with the HS approximation gives slightly higher pressure on the isentrope than with the MF one, but both curves are very close to the ZA isentrope. The experimental point on the adiabatic release of Al into aerogel [19] is quite close to the calculated isentropes; the difference is remarkable only in large-scale figure 2(b).

Molybdenum
We calculated a release isentrope of molybdenum from the shocked initial state with an experimental particle velocity u pH = 2.73 km/s [27]. Other parameters at the initial point were determined by QMD: ρ H = 15 g/cm 3 , T H = 5.3 kK, and P H ≈ 2500 kbar. The QMD shock Hugoniot and release isentrope are shown in figure 3. The isentrope was restored by only one RH method; good agreement with experiment [27] was obtained. We also show the release isentrope with the same initial point computed from a multiphase EOS by Fortov and Lomonosov [28]. It can be easily seen that the experimental points and theoretical curves are very close to 6 1234567890 ''""  Figure 3. Release isentrope of Mo in pressure-particle velocity plane. Initial point corresponds to experiment [27]. (a) Blue line-multiphase EOS release isentrope [28]; black line with black circles-QMD Hugoniot; blue open squares-release isentrope calculated by the RH method; red triangles-experimental data [27]. (b) Linear pressure scale. Denotations are the same as in figure 3(a).

Conclusions
• All the three methods of calculation of release isentropes for Al showed good agreement with experiment and with each other. • The ZA is the only exact method among the other two; however, it is computationally effective only if many isentropes are required. • The RH method is the fastest; its use is justified if only one or several isentropes are needed.
• The DR method is computationally expensive but gets direct information about entropy and free energy. • The RH method gives satisfactory results even for a small number of intermediate points along a release isentrope. • The DR method with the MF approximation for the gas part gives a better match with experiment in comparison with the HS approximation. • QMD calculations provide for good accuracy of thermodynamic properties of matter and can be used for calibration of semiempirical EOSs in case of lack of experimental data.