Benchmarking vdW-DF first principle predictions against Coupled Electron-Ion Monte Carlo for high pressure liquid hydrogen
AArticle Type
Benchmarking vdW-DF first principle predictions againstCoupled Electron-Ion Monte Carlo for high pressure liquidhydrogen
Vitaly Gorelov, Carlo Pierleoni*, and David M. Ceperley Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay, 91191 Gif-sur-Yvette,France Department of Physical and Chemical Sciences, University of L’Aquila, Italy Department of Physics, University of Illinois Urbana-Champaign, Illinois, USA
Correspondence: *Carlo Pierleoni, Email: [email protected],[email protected]
Summary
We report first principle results for nuclear structure and optical responses of high pressure liquid hydrogen along twoisotherms in the region of molecular dissociation. We employ Density Functional Theory with the vdW-DFapproximation (vdW) and we benchmark the results against existing predictions from Coupling Electron-Ion MonteCarlo (CEIMC). At fixed density and temperature, we find that pressure from vdW is higher than pressure fromCEIMC by about 10 GPa in the molecular insulating phase and about 20 GPa in the dissociated metallic phase.Molecules are found to be overstabilized using vdW, with a slightly shorter bond length, and with a strongerresistance to compression. As a consequence, pressure dissociation along isotherms using vdW is more progressivethan computed with CEIMC. Below the critical point, the liquid-liquid phase transition is observed with boththeories in the same density region but the one predicted by vdW has a smaller density discontinuity, i.e. a smallerfirst order character. The optical conductivity computed using Kubo-Greenwood is rather similar for the two systemsand reflects the slightly more pronounced molecular character of vdW.
Keywords: high pressure hydrogen, liquid-liquid phase transition, liquid structure, quantum monte carlo methods
Condensed hydrogen at high pressure is interesting because of the interplay between molecular dissociation andmetallization in the liquid phase at high temperatures, and between other phase transitions in the solid phase at lowertemperatures [1] . First principle methods are the main theoretical tools for investigation of this system. The proximity a r X i v : . [ c ond - m a t . m t r l - s c i ] D ec Vitaly Gorelov et al of metallization and the dispersive character of the molecular interactions in the insulating phase have proven to bechallenging for the well-established Density Functional Theory (DFT) based methods [2–7] . Predictions from variousapproximations are in qualitative agreement but there are substantial quantitative differences which makes it difficultto make accurate predictions without benchmarking DFT results against more fundamental theories. At the same time,the simple electronic structure of hydrogen makes it the ideal candidate for developing new first principle methodsbased on more fundamental many-body theories such as those based on Quantum Monte Carlo techniques. Two similarmethods have appeared in the last decade and were applied to study compressed hydrogen: the Coupled Electron-IonMonte Carlo (CEIMC) method [2, 8–10] and the Quantum Monte Carlo Molecular Dynamics (QMCMD) [11, 12] . Theyare both based on solving the electronic problem by ground state Quantum Monte Carlo methods but they differ inthe way they sample the nuclear configuration space . The advantage of QMC over DFT is in its variational character:very different implementations of QMC, i.e. different forms of the many body wave function, can be directly comparedand ranked based on the value of their total energy . After several attempts, these two methods are now in agreementabout the occurrence and the location of the liquid-liquid phase transition (LLPT) in hydrogen [12] . This reinforcesthe idea that QMC can be used to assess the accuracy of DFT functionals for finite temperature predictions.Previous detailed investigations of the accuracy of several DFT approximations against ground state QMC forboth crystalline and liquid static configurations [5] found van-der-Waals (vdW-DF) [13] and Heyd-Scuseria-Ernzerhof(HSE) [14] approximations to be the best candidates among the ones considered, with energy biases within fractions ofmH/atoms and pressure biases within 10-20 GPa. A CEIMC investigation of the LLPT in hydrogen [10] confirmedthese expectations: first principle molecular dynamics (FPMD) with vdW-DF and HSE approximations predictLLPT transitions in closer agreement with the results from CEIMC than other approximations. Differences in theproton-proton pair correlation functions between CEIMC and FPMD are observed [2, 7, 15] . These differences couldprovide insight into the dissociation mechanism.In this paper we compare the thermodynamics and local nuclear structure from first principle simulations withvdW-DF and CEIMC, along two isotherms of hydrogen across the molecular dissociation region. These isotherms areabove and below the critical temperature of the LLPT. A comparison of the optical conductivity within the Kubo-Greenwood framework is made. We only consider vdW-DF and disregard HSE since the computational requirementsof HSE are much larger than of vdW-DF, while the expected accuracy is comparable [5] .The paper is organized as follows. First we briefly describe the methods employed in section 2. Then we report ourresults for the thermodynamics and for the local nuclear structure in section 3. In the following section 4 we comparethe optical conductivity and finally we discuss the results and draw our conclusions. CEIMC sample the nuclear configuration space by Metropolis Monte Carlo which requires computing electronic total energies, whileQMCMD implements a Langevin Dynamics based on electronic forces on the nuclei. The variance of the total energy is also an important quantity to compare since it decreases when the quality of the wave function isimproved. Therefore a better wave function has a lower energy and a smaller variance. italy Gorelov et al In this paper we used Metropolis Monte Carlo sampling for the nuclear degrees of freedom and, since the focus is onthe accuracy of the potential energy surface, we limit our discussion to classical protons. The electronic energies usedin the Metropolis sampling arise either from ground state QMC in the CEIMC [1, 8] , or from ground state DensityFunctional Theory with the vdW-DF approximation. In both methods the sampling is performed by a force-biasedMC algorithm with DFT-forces [16] . Given a nuclear configuration (cid:126)R of the N p protons (here vectors are in the N p dimensional space), a new configuration (cid:126)R (cid:48) is proposed by a drifted random displacement according to (cid:126)R (cid:48) = (cid:126)R + h (cid:126)F ( (cid:126)R ) + (cid:126)ξ (1)where (cid:126)F ( (cid:126)R ) is the total force acting on the protons, h is an adjustable parameter with suitable dimensions, and (cid:126)ξ agaussian distributed random vector with the properties (cid:104) (cid:126)ξ (cid:105) = 0 (cid:104) (cid:126)ξ (cid:126)ξ (cid:48) (cid:105) = 2 h δ (cid:126)ξ,(cid:126)ξ (cid:48) (2)where is the identity matrix and δ (cid:126)ξ,(cid:126)ξ (cid:48) indicates that the random noises at different steps are uncorrelated. Then thesampling probability is G ( (cid:126)R → (cid:126)R (cid:48) ; h ) ∝ e − [ (cid:126)R (cid:48)− (cid:126)R − h(cid:126)F ( (cid:126)R ) ] h (3)while the acceptance probability, satisfying detailed balance, is A ( (cid:126)R → (cid:126)R (cid:48) ) = min (cid:104) , q ( (cid:126)R → (cid:126)R (cid:48) ) (cid:105) q ( (cid:126)R → (cid:126)R (cid:48) ) = e − βE ( (cid:126)R (cid:48) ) G ( (cid:126)R (cid:48) → (cid:126)R ; h ) e − βE ( (cid:126)R ) G ( (cid:126)R → (cid:126)R (cid:48) ; h ) . (4)In the last equation E ( (cid:126)R ) is the total energy of nuclear configuration (cid:126)R . Within the DFT framework, this scheme canbe used straightforwardly, both forces and energies are from DFT-vdW-DF functional. This strategy allows one tohave a relatively high acceptance probability and also large moves even for system of hundreds of protons, hence anefficient sampling of the nuclear configurational space.When solving the electronic problem with QMC, energy and forces have statistical noise that must taken intoaccount to avoid biasing the sampling. In CEIMC the noise over electronic energy is accounted for with the Penalty Vitaly Gorelov et al method [8, 17] while forces are not computed. Our present CEIMC implementation uses a two-level Metropolis schemein which we first propose a new configuration according to the scheme described above based on DFT energy andforces, if the new configuration (cid:126)R (cid:48) is accepted we perform a second Metropolis acceptance test according to A ( (cid:126)R → (cid:126)R (cid:48) ) = min (cid:104) , q ( (cid:126)R → (cid:126)R (cid:48) ) (cid:105) q ( (cid:126)R → (cid:126)R (cid:48) ) = e − βE ( (cid:126)R ) e − βE ( (cid:126)R (cid:48) ) e − [ βδ ( (cid:126)R (cid:48) , (cid:126)R )+ β u ( (cid:126)R (cid:48) , (cid:126)R ) ] (5)where δ ( (cid:126)R (cid:48) , (cid:126)R ) indicates the estimate of the difference between the QMC energies of configurations (cid:126)R (cid:48) and (cid:126)R , and u ( (cid:126)R (cid:48) , (cid:126)R ) is a positive term, hence a penalty in the acceptance coming from the noise which, to first approximation,can be replaced by the estimate of the variance of the energy difference [8, 17] .We do not report new CEIMC calculations but use past calculations as a reference. Details about these calculationsare reported in refs. [10, 15, 18] . We do report new calculations using DFTMC employing the vdW-DF XC functionaland the plane augmented-wave (PAW) method [19] with an energy cutoff of 40 Rydberg. In order to have a directcomparison with the CEIMC results, the DFT calculations have been performed with exactly the same protocolemployed in CEIMC. Specifically we did simulations of periodic systems of N p = 54 protons, with a 4x4x4 Monkhorst-Pack grid of k-points. This choice ensures that both results are affected by the same finite size error at the single electronlevel. In CEIMC, many-body size effects are estimated and thermodynamic results corrected accordingly [10, 20]3 . Afurther check that residual finite size effects in the DFT calculations are small is provided by direct comparison withcalculations for a system of N p = 128 protons as detailed in the next section. The DFTMC calculations are performedwith our in-house CEIMC code (BOPIMC) which uses the PWSCF [21] package as a DFT solver.The electrical conductivity was computed using the Kubo-Greenwood (KG) [22, 23] formula on 40 randomly selectedconfigurations. The configurations are taken from DFTMC and CEIMC runs with N P = 54 protons for a given density( ρ ) and temperature ( T ). A post-processing tool in the QuantumEspresso suite [21] , KGEC [24] was employed. Forthese calculations we used an 8x8x8 Monkhorst-Pack grid of k-points to get a better convergence in the low- ω part ofthe conductivity and included 48 bands to converge in the high- ω part up to ω (cid:39) eV . We performed simulations along two isotherms T = 1500 K, K , one below and one above the critical temperatureof the LLPT. (The precise value of T c is yet to be determined accurately.) In fig. 1 we compare the equation of state(EOS) along the two isotherms from CEIMC and from DFTMC-vdW-DF. The data shown in fig 1 for CEIMC arefrom VMC and include finite size corrections and RQMC corrections as explained in the SM of ref. [10] (see also [20] ). Many-body size effects on the local nuclear structure are negligible. italy Gorelov et al
160 170 180 190 200 210 220 230 240 250 260 270 1.36 1.38 1.4 1.42 1.44 1.46 1.48 1.5 P ( G P a ) r s T=1500K
60 80 100 120 140 160 180 200 1.45 1.5 1.55 1.6 1.65 1.7 P ( G P a ) r s T=3000K P d ft - P q m c ( G P a ) r s T=1500K
10 12 14 16 18 20 22 24 1.45 1.5 1.55 1.6 1.65 1.7 P d ft - P q m c ( G P a ) r s T=3000K
Figure 1: EOS along the two isotherms investigated. The upper panels report the pressure versus r s from CEIMC(blue circles) and from DFTMC with N p = 54 (red squares). At T = 1500 K DFTMC for systems with N p = 128 arereported (black triangles). In the lower panels the difference between DFT and QMC pressures is shown. Symbols areas in the upper panels.The corresponding numerical values are given in the SM of ref. [10] . At the lower temperature, we observe a weaklyfirst order phase transition in both theories in the same range of densities. However the width of the density plateau,obtained as the shift between the two branches of the EOS from fitting the data, is about two times larger in CEIMCthan in DFT, signaling a stronger first order character of the transition in QMC than in DFT. Residual finite sizeeffects, that could arise from the limited size of our systems, in particular at the phase transition, are negligibly smallas shown in the upper left panel of the figure. The pressure from vdW-DF is always larger than from CEIMC. Moreoverthe accuracy of the vdW-DF predictions for the EOS is not uniform across molecular dissociation. This is shown in Vitaly Gorelov et al
012 r s =1.45 r s =1.44 r s =1.43012 0 1 2 3 4 5r s =1.42r (a ) 0 1 2 3 4 5r s =1.41r (a ) 0 1 2 3 4 5 6r s =1.40r (a ) Figure 2: Proton-proton radial distribution functions along the T = 1500 K isotherm. Comparison between CEIMC(blue lines) and DFTMC-vdW-DF with N p = 54 (red lines). At r s = 1 . , . , . , . results for systems of N p = 128 with DFTMC-vdW-DF are also reported (black dashed line).the two lower panels of the figure where we report the difference between DFT and QMC pressures. At T = 1500 K ,molecular dissociation occurs suddenly between r s = 1 . and r s = 1 . where the pressure difference jumps from ∼ − GPA to ∼ − GPa. At T = 3000 K molecular dissociation, as well as the change in pressure difference, aremore gradual with density. However, the observed pressure difference is in the same range of GP a ≤ ∆ P ≤ GP a .In fig. 2 we compare the proton-proton g ( r ) along the T = 1500 K isotherm. We observe a qualitative agreementbetween the two theories, but data from DFTMC exhibit a more pronounced molecular character in the entiredensity range. The molecular bond length in DFTMC theory is slightly shorter than in CEIMC, and the moleculardissociation is more progressive upon increasing density, in agreement with the smaller density discontinuity observedin the EOS. The same comparison along the T = 3000 K isotherm is reported in fig. 3. Here, in both theories, themolecular dissociation is a progressive process. As before, the molecular peak is slightly shorter and the molecularcharacter is slightly more pronounced (more correlation) in DFTMC than in CEIMC. The overall agreement is betterthan at lower temperature.In order to investigate in more detail the proton pairing in the system, we have used a cluster algorithm to identify,for any given nuclear configuration, the number of bound pairs. A similar analysis for CEIMC has been alreadyreported in ref. [15] . Specifically given a proton, the algorithm looks for the closest proton which has not been alreadypaired to another proton at shorter distance. The algorithm find all such pairs within a given cutoff distance.italy Gorelov et al
012 r s =1.70 r s =1.66 r s =1.65012 r s =1.64 r s =1.63 r s =1.62012 r s =1.61 r s =1.60 r s =1.59012 0 1 2 3 4 5r s =1.58r (a ) 0 1 2 3 4 5r s =1.55 r(a ) 0 1 2 3 4 5 6r s =1.52r (a ) Figure 3: Proton-proton radial distribution function along the T = 3000 K isotherm. Comparison between CEIMC(blue lines) and DFTMC-vdW-DF with N p = 54 (red lines).In the first analysis we have used a large cutoff distance, r c = 2 . a in order to compute the bond length distributionwithout artificially cutting off the large distance tail. Results of this analysis for both CEIMC and DFTMC trajectoriesalong the two isotherms are illustrated in the upper panels of figs. 4 and 5 where we report the bond distributionfunction. In each figure, the left panels correspond to the CEIMC results, while the right panels to DFTMC. At T = 1500 K we see a change in behavior for both theories at the LLPT, the discontinuity being more pronounced inCEIMC, but still clearly visible in DFTMC . In the molecular phase the observed distribution is rather independenton the specific density, the maximum of probability is at r = 1 . − . a in CEIMC and at r = 1 . − . a inDFTMC, and the shape of the distribution is rather similar in the two methods. A similar difference in the molecularbond length between vdW-DF and QMC was previously observed in the crystalline state (see fig. 4 of ref. [5] ). Atdensities beyond molecular dissociation the bond distribution exhibits a more pronounced density dependence in boththeories but remains slightly more localized in DFTMC than in CEIMC (the maximum of probability in DFTMC isat r ∼ . a while it is at r ∼ . a in CEIMC). At T = 3000 K a more progressive change of the bond distribution Convergence of results at the coexistence density is problematic since the system is observed to switch from one phase to the otherseveral times during the simulation.
Vitaly Gorelov et al a) P ( d ) d/a b) P ( d ) d/a c) P ( N pa i r s ) N pairs d) P ( N pa i r s ) N pairs Figure 4: T = 1500 K isotherm. Bond length distributions obtained with r c = 2 . (panels a) and b)) and the probabilityof the number of pairs at distance r ≤ . a (panels c) and d). CEIMC results panels a) and c); DFTMC results,panels b) and d).is observed over the investigated density range. The two theories have similar behavior but again the CEIMC resultsexhibit a stronger density dependence.In the second analysis we used a cutoff distance of r c = 1 . a (the first minimum of the g ( r ) ) and we computed thedistribution in the number of the bonded pairs found within this cutoff. Results for both CEIMC and DFTMC alongthe two isotherms are illustrated in the lower panels of figs. 4 and 5. Also this distribution exhibits a sudden changeof shape at the LLPT along the T = 1500 K isotherm, as seen in fig. 4. Again the two theories are in qualitativeagreement. More quantitatively, DFTMC shows a more localized distribution in the “dissociated” phase after the LLPTcompared to CEIMC, hence a stronger residual associated character. At higher temperature the change of the shapeof the distribution with density is gradual in both theories. However, as observed at lower temperature, the CEIMCdistributions become wider than the DFTMC distributions for increasing density signaling a less associated character.italy Gorelov et al a) P p ( d ) d/a b) P p ( d ) d/a c) P ( N pa i r s ) N pairs d) P ( N pa i r s ) N pairs Figure 5: T = 3000 K isotherm. Bond length distributions obtained with r c = 2 . (panels a) and b)) and the probabilityof the number of pairs at distance r ≤ . a (panels c) and d). CEIMC results panels a) and c), DFTMC results,panels b) and d). Optical properties are relevant since they are the main source of experimental information at these high pressure-high temperature states of liquid hydrogen. In particular recent experiments at NIF [25] provided support for theLLPT lines predicted by CEIMC and FPMD-vdW-DF for deuterium on the basis of the optical response of thesystem . A subsequent investigation, within the Kubo-Greenwood framework, of optical properties across moleculardissociation support experimental findings [26] , although other experiments suggest different scenarios [7] and thedebate is still ongoing. For the purposes of the present benchmark, we believe it is interesting to check the sensitivityof optical properties to the details of the local liquid structure. To this aim we have computed the optical conductivityat three distinct phase points, ( r s = 1 . , T = 1500 K ) , ( r s = 1 . , T = 1500 K ) and ( r s = 1 . , T = 3000 K ) for liquidstructure from both methods. At each state point and for each simulation method, we have considered 40 independentnuclear configurations generated during the simulation run, and for each nuclear configuration we have computed the The LLPT lines from CEIMC and from FPMD-vdW-DF are closer than the experimental resolution. et al , eV , ( c m ) a) vdW-DFCEIMC , eV D e n s i t y o f S t a t e s , ( a r b . un i t s ) b) vdW-DFCEIMC Figure 6: T = 1500 K, r s = 1 . . a) Real part of the average optical conductivity from CEIMC and DFTMC. b) Averageelectronic density of states from CEIMC and DFTMC. Statistical errors are represented by the thickness of the lines.optical conductivity with the Kubo-Greenwood single electron formalism, assuming Kohn-Sham eigenstates fromthe vdW-DF XC functional. Details are given in section 2. Statistical averages and errors are obtained over thesample of configurations. Figure 6a) compares the real part of the average optical conductivity from CEIMC andfrom DFTMC-vdW-DF at the first state point, while fig. 6b) reports the average electronic density of states fromboth theories. Results from the two methods are in good agreement showing, once more, the similarity betweenthe generated nuclear structures. However, we observe a slightly larger density of states at the Fermi level, whichcorrelates with a slightly higher conductivity at small energies, from DFTMC-vdW-DF than from CEIMC. Converselyaround the maximum of conductivity at ω (cid:39) eV , CEIMC data are slightly higher than DFTMC-vdW-DF data. Notethat this state point is at the edge of the metallic regimes: the energy gap is closed but the density of states at theFermi level is rather small, resulting in a very small DC conductivity.A similar analysis into the conducting regime is presented in figs. 7 where optical conductivities for the other twoinvestigated state points are reported. We observe similar comparisons at both state points: CEIMC conductivity islarger at small energies and the maximum is at slightly lower energy than the DFTMC-vdW-DF one. This differencecan be qualitatively justified as arising from the less “molecular” character of the CEIMC configurations (see figs.2 and 3), hence from a weaker electronic localization. Analysis of the conductivities in terms of the Drude-Smithmodel [27, 28] provides σ dft ( ω = 0) = 5450(50) , σ qmc ( ω = 0) = 7280(70)(Ω cm ) − at T = 1500 K r s = 1 . and σ dft ( ω = 0) = 4600(40) , σ qmc ( ω = 0) = 5900(40)(Ω cm ) − at T = 3000 K r s = 1 . for vdW-DF and CEIMCrespectively. In figure 7b) we also report a recent result for the DC conductivity from FPMD with vdW-DF XCapproximation [29] at a state point close to our present one . The comparison is shown to support our present analysis. We believe that the value is slightly larger than our own because thedensity of that state point is slightly larger. italy Gorelov et al , eV , ( c m ) a) vdW-DFCEIMC , eV , ( c m ) b) vdW-DFCEIMCKnudson et. al. Figure 7: a) Real part of the average optical conductivity at r s = 1 . and T = 1500 K from CEIMC and DFTMC. b)Real part of the average optical conductivity at r s = 1 . and T = 3000 K from CEIMC and DFTMC. Dashed linesrepresent a Drude-Smith fit to data. Statistical errors are represented by the thickness of the lines. In conclusion, we have presented a detailed comparison between CEIMC and DFTMC with vdW-DF approximationfor high pressure liquid hydrogen across the interesting region where pressure induced molecular dissociation occurs,either as a continuous process at high temperature, or through a first order phase transition at lower temperature. Wehave computed the EOS, the local nuclear structure of the fluid and the optical conductivity at several state points.We observe that pressure from vdW-DF is higher than from CEIMC by about GP a in the molecular phase andabout GP a in the dissociated phase, indicating a loss of accuracy of the approximation upon molecular dissociation.Inspection of the local nuclear structure, as monitored by the proton-proton g ( r ) , shows that the molecular characterof the fluid is more pronounced with DFTMC-vdW-DF and molecules have more resistance to compression and aremore difficult to dissociate. Despite these quantitive difference, our investigation demonstrate once more the goodquality of vdW-DF functional to study high-pressure hydrogen in this interesting region of its phase diagram. Acknowledgments
V. G. and C.P. were supported by the Agence Nationale de la Recherche (ANR) France, under the program“Accueil de Chercheurs de Haut Niveau 2015” project: HyLightExtreme. D.M.C. was supported by DOE Grant NADE-NA0002911 and by the Fondation NanoSciences (Grenoble). Computer time was provided by PRACE Project2016143296 and by an allocation of the Blue Waters sustained petascale computing project, supported by the NationalScience Foundation (Award OCI 07- 25070) and the State of Illinois, and by the HPC resources from GENCI-CINESunder the allocation 2018-A0030910282.2 REFERENCES
References [1] J. M. McMahon, M. A. Morales, C. Pierleoni, D. M. Ceperley,
Rev. Mod. Phys. , , 1607–1653.[2] M. A Morales, C. Pierleoni, E. Schwegler, D M Ceperley, Proc. Nat. Acad. Sc. , ( ), 12799–12803.[3] M. A. Morales, J. M. McMahon, C. Pierleoni, D. M. Ceperley, Phys. Rev. B , , 184107.[4] M. A. Morales, J. M. McMahon, C. Pierleoni, D. M. Ceperley, Phys. Rev. Letts. , ( ), 65702.[5] R. C. Clay, J. McMinis, J. M. McMahon, C. Pierleoni, D. M. Ceperley, M. A. Morales, Phys. Rev. B , ( ), 184106.[6] S. Azadi, W. M C Foulkes, T. D. Kuhne, New Journal of Physics , .[7] M. D. Knudson, M. P. Desjarlais, A. Becker, R. W. Lemke, K. R. Cochrane, M. E. Savage, D. E. Bliss, T. R. Mattsson, R. Redmer, Science , , 1455–1460.[8] C. Pierleoni, D. M. Ceperley , p. 641–683., Lecture notes in physics: Springer, Berlin.[9] N. M. Tubman, E. Liberatore, C. Pierleoni, M. Holzmann, D. M. Ceperley, Phys. Rev. Lett. , ( ), 1–7, available at 1408.6523.[10] C. Pierleoni, M. A Morales, G. Rillo, M. Holzmann, D. M Ceperley, Proc. Natl. Acad. Sci. , ( ), 4954–4957.[11] C. Attaccalite, S. Sorella, Phys. Rev. Letts. , ( March ), 1–4.[12] G. Mazzola, R. Helled, S. Sorella,
Physical Review Letters , , 025701.[13] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, B. I. Lundqvist, Phys. Rev. Letts. , ( ), 246401–1, available at 0402105.[14] J. Heyd, J. E Peralta, G. E Scuseria, R. L Martin, J. Chem. Phys. , ( ), 174101.[15] C. Pierleoni, M. Holzmann, D. M. Ceperley, Contrib. to Plasma Phys. , , 99–106.[16] M P Allen, D J Tildesley, Computer Simulation of liquids , Oxford University Press, New York .[17] D Ceperley, M Dewing,
J. Chem. Phys. , , 9812.[18] G. Rillo, M. A. Morales, D. M. Ceperley, C. Pierleoni, J. Chem Phys. , , 102314, available at 1708.07344.[19] P E Blochl, Phys. Rev. B , , 17953.[20] M. Holzmann, R. C. Clay, M. A. Morales, N. M. Tubman, D. M. Ceperley, C. Pierleoni, Phys. Rev. B , , 035126.[21] P. Giannozzi, O Andreussi, T Brumme, O Bunau, M Buongiorno Nardelli, M Calandra, R Car, C Cavazzoni, D Ceresoli, M Cococcioni,N Colonna, I Carnimeo, A Dal Corso, S de Gironcoli, P Delugas, R. A. J. DiStasio, A Ferretti, A Floris, G Fratesi, G Fugallo, RGebauer, U Gerstmann, F Giustino, T Gorni, J Jia, M Kawamura, H-Y Ko, A Kokalj, E Kukukbenli, M Lazzeri, M Marsili, NMarzari, F. Mauri, N L Nguyen, H-V Nguyen, A Oreto-de-la Roza, L Laulatto, S Ponce, D Rocca, R Sabatini, B Santra, M Schlipf,A. Seitsonen, A Smogunov, I Timrov, T. Thonhauser, P Umari, N Vast, X. Wu, S. Baroni, J. Phys. Cond Mat. , , 465901.[22] R. Kubo, Journal of the Physical Society of Japan , ( ), 570–586, available at 0211006.[23] D A Greenwood, Proc. Phys. Soc. , ( ), 585–596. EFERENCES 13 [24] L. Calderín, V. V. Karasiev, S. B. Trickey,
Comput. Phys. Commun. , , 118–142, available at 1707.08437.[25] P. M Celliers, M. Millot, S. Brygoo, R. S. McWilliams, D. E Fratanduono, J R. Rygg, A. F. Goncharov, P. Loubeyre, J. H Eggert,J. L. Peterson, N. B. Meezan, S. Le Pape, G. W Collins, R. Jeanloz, R. J. Hemley, Science , ( August ), 677–682.[26] G. Rillo, M. A. Morales, D. M. Ceperley, C. Pierleoni, submitted .[27] N. V. Smith,
Phys. Rev. B - Condens. Matter Mater. Phys. , ( ), 155106.[28] T. L. Cocker, D. Baillie, M. Buruma, L. V. Titova, R. D. Sydora, F. Marsiglio, F. A. Hegmann, Phys. Rev. B , ( ), 205439,available at 1705.10350.[29] M D Knudson, M P Desjarlais, M Preising, R Redmer, Physical Review B , , 174110., 174110.