# Fileset [2301.03570v2.pdf](https://mdr.nims.go.jp/filesets/0e29fb62-2770-4099-91b7-5e2215089e90/download) ## Creator Giacomo Tenti, [Kousuke Nakano](https://orcid.org/0000-0001-7756-4355), Andrea Tirelli, Sandro Sorella, Michele Casula ## Rights arXiv.org perpetual, non-exclusive license 1.0[In Copyright](http://rightsstatements.org/vocab/InC/1.0/) ## Other metadata [Principal deuterium Hugoniot via quantum Monte Carlo and

Δ

-learning](https://mdr.nims.go.jp/datasets/5e868547-9efe-4207-ba65-a111375f7de4) ## Fulltext Principal deuterium Hugoniot via Quantum Monte Carlo and ∆-learningGiacomo Tenti,1, ∗ Kousuke Nakano,2, † Andrea Tirelli,1 Sandro Sorella,1 and Michele Casula31International School for Advanced Studies (SISSA), Via Bonomea 265, 34136 Trieste, Italy2Center for Basic Research on Materials, National Institute for Materials Science (NIMS), Tsukuba, Ibaraki 305-0047, Japan3Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC),Sorbonne Université, CNRS UMR 7590, MNHN, 4 Place Jussieu, 75252 Paris, France(Dated: May 7, 2024)We present a study of the principal deuterium Hugoniot for pressures up to 150 GPa, usingMachine Learning potentials (MLPs) trained with Quantum Monte Carlo (QMC) energies, forcesand pressures. In particular, we adopted a recently proposed workflow based on the combinationof Gaussian kernel regression and ∆-learning. By fully taking advantage of this method, weexplicitly considered finite-temperature electrons in the dynamics, whose effects are highly relevantfor temperatures above 10 kK. The Hugoniot curve obtained by our MLPs shows a good agreementwith the most recent experiments, particularly in the region below 60 GPa. At larger pressures,our Hugoniot curve is slightly more compressible than the one yielded by experiments, whoseuncertainties generally increase, however, with pressure. Our work demonstrates that QMC canbe successfully combined with ∆-learning to deploy reliable MLPs for complex extended systemsacross different thermodynamic conditions, by keeping the QMC precision at the computational costof a mean-field calculation.I. INTRODUCTIONThe study of hydrogen under extreme conditions hasbeen a very active topic in condensed matter physics.Hydrogen is the most abundant element in the universeand the accurate knowledge of its phase diagram atpressures of the order of hundreds of GPa is extremelyimportant for a variety of applications, such as modelingthe interior of stars and giant gas planets [1–3], theinertial-confinement fusion [4], and the high-Tc hydrogen-based superconductors [5, 6]. Nevertheless, severalproperties of this system are still highly debated, evenat the qualitative level [7–10].One of the main reasons that hamper our fullunderstanding of high-pressure hydrogen is the difficultyof reproducing extreme pressures in a laboratory. Typicalshock-wave experiments [11] make use of acceleratedflyer plates to compress a material sample in a veryshort time, thus allowing the study of specimens athigh temperatures and pressures. In particular, the setof possible end-states that the system can reach fromsome given initial conditions, also named the principalHugoniot, must satisfy a set of equations, known asthe Rankine-Hugoniot (RH) relations [12], linking thethermodynamic properties of the final shocked state withthose of the starting one. During the years, the principaldeuterium Hugoniot has been measured for a wide rangeof pressures and with a great degree of accuracy [13–20],particularly in the lower pressure (≲ 60 GPa) region,where the relative error on the density is as small as 2%in recent experiments.In this context, numerical approaches, in particularAb Initio Molecular Dynamics (AIMD) simulations, are∗ gtenti@sissa.it† kousuke˙1123@icloud.comextremely valuable, since they are not constrained by anyexperimental setup and can thus give further insight intothis part of the phase diagram [21]. The Hugoniot regionis particularly important because of the availabilityof experimental data that can be used to benchmarkdifferent theoretical methods. Among them, DensityFunctional Theory (DFT) has been extensively applied tocompute the Hugoniot curve [22–28]. In this framework,the approximations behind the particular exchange-correlation functional often produce discrepancies acrossexisting DFT schemes, whose accuracy varies accordingto the thermodynamic conditions. Quantum MonteCarlo (QMC) simulations, which depend on morecontrollable approximations, have also been performed[29, 30]. Although in principle more accurate andsystematically improvable, these calculations have amuch larger computational cost than DFT, and theyare thus limited in system size and simulation length.Moreover, previous QMC calculations [29] seem to giveresults for the principal Hugoniot in disagreement withthe most recent experimental data, with the possibleorigin of this discrepancy being recently debated [31].To overcome the large computational cost of abinitio simulations, machine learning techniques, aimedat constructing accurate potential energy surfaces, havebecome increasingly popular. Within this approach, oneuses a dataset of configurations, i.e. the training set,to build a machine learning potential (MLP) that isable to reproduce energies and forces calculated withthe given target method [32]. Unlike DFT MLPs,the QMC ones are relatively less common, given thelarger computational cost and the consequent difficulty ofgenerating large datasets, usually necessary to constructaccurate MLPs.In this work, we have built an accurate MLP usingQMC energies, forces and pressures in the region ofthe deuterium Hugoniot, using the so-called ∆-learningarXiv:2301.03570v2 [cond-mat.str-el] 3 May 20242approach. The Hugoniot curve computed by our QMC-MLP shows an excellent agreement with the most recentexperiments at low density, while presenting a slightlylarger compressibility for temperatures above 10 kK,where the experimental points are also affected by largererror bars.II. COMPUTATIONAL DETAILSIn order to build an MLP with QMC references, weemployed a combination of Gaussian Kernel Regression(GKR), Smooth Overlap of Atomic Positions (SOAP)descriptors [33], and ∆-learning. The same approachhas been recently proposed in Ref. 34, where it wasapplied to the study of high-pressure hydrogen insimilar thermodynamic conditions. Following the ∆-learning approach, an MLP is trained on the differencebetween the target method and a usually much cheaperbaseline potential. Here, we trained several MLPs, usingVariational Monte Carlo (VMC) and Lattice RegularizedDiffusion Monte Carlo (LRDMC) [35, 36] datapointsas targets, and two different DFT baselines, withthe Perdew-Zunger Local Density Approximation (PZ-LDA) [37] and the Perdew-Burke-Ernzerhof (PBE) [38]functionals.To determine the principal Hugoniot, we made use ofthe RH jump equation:H(ρ, T ) = e(ρ, T )− e0 +12(ρ−1 − ρ−10 ) [p(ρ, T ) + p0] = 0,(1)where ρ, T , e(ρ, T ), p(ρ, T ) and ρ0, T0, e0, p0 are thedensity, temperature, energy per particle and pressure ofthe final and initial states, respectively. In particular, weran a first set of NV T simulations for a system of N =128 atoms at several temperatures in the [4 kK : 8 kK]range, and Wigner-Seitz radii between 1.80 Bohr and2.24 Bohr, corresponding to the range where the zero ofH(ρ, T ) was expected. These simulations were performedconsidering classical nuclei and ground-state electrons, asquantum corrections and thermal effects have been shownto be negligible for these temperatures [30]. At each step,the energy, forces and pressure were calculated using theQuantum Espresso package in its GPU acceleratedversion [39–41] with the chosen functional (PBE in mostcases), and then corrected with our MLP trained on thedifference between QMC and DFT data. The resultingdynamics has the same efficiency as a standard DFTAIMD simulation, which is roughly 100 times faster thanthe original QMC one. For the DFT simulations, weconsidered a 60 Ry plane-wave cutoff with a ProjectorAugmented Wave (PAW) pseudopotential [42] and a4 × 4 × 4 Monkhorst-Pack k-point grid, while for thedynamics we used a time step of 0.25 fs and a Langevinthermostat [43, 44] with damping γ = 13 ps−1. Foreach temperature, the Hugoniot (ρ∗, p∗) coordinates aredetermined by fitting the Hugoniot function H(ρ, T )and the pressure p(ρ, T ) with a spline function, and bynumerically finding ρ∗ and the corresponding p∗.The QMC calculations were performed using theTurboRVB package [45, 46]. When generating ourQMC-MLP model, we took particular care to the trainingset construction, based on energies, forces and pressuresall computed at the QMC (VMC or LRDMC) level. Ingeneral, QMC forces and pressure are not guaranteed tobe consistent with the relative potential energy surface,due to the practical difficulty of optimizing all theparameters within a given ansatz for the VMC wavefunction (WF). This is often called the self-consistencyerror [47–49]. Even if ML frameworks generally satisfythe consistency property by construction, the presenceof biased forces and pressure in the training can spoilthe accuracy of the model and produce, in principle,unintended results. Therefore, to avoid these issues,we mitigated the self-consistency error by directlyoptimizing not only the Jastrow factor but also thedeterminantal part of the VMC WF. The details ofour QMC simulations are reported in the SupplementalMaterial (SM) [50].Within our approach, we can fully take advantageof the ∆-learning method by estimating the effect ofthermalized electrons in our calculations. To do so, weran simulations at temperatures T = 10 kK, 15 kK,20 kK and 35 kK considering the effect of finite T inthe underlying PBE energies and forces. In this way, wecan include the effects of thermally excited electrons inour MLP without changing it, at least at the DFT levelof theory. Here, the effect of the explicit dependence ofthe DFT functional on T was not considered, since it hasbeen shown to be negligible for hydrogen systems and thetemperature range analyzed here [28, 51, 52]. We remarkthat this approach is also applicable when a DFT-MLPis used as a baseline in place of an ab initio calculation,where finite temperature effects can be estimated fromthe DFT density of states [53].III. RESULTS AND DISCUSSIONFig. 1a shows our results together with severalexperimental values for pressures below 150 GPa [16,19, 20]. We also report the principal Hugoniot obtainedby directly using the PBE baseline and the CoupledElectron Ion Monte Carlo (CEIMC) results of Ref. 30for comparison. For temperatures larger than T =10 kK the results refer to the MDs obtained using finitetemperature DFT as a baseline. Both the VMC andLRDMC models give a very similar Hugoniot line, wellreproducing the experimental points in the low density -low pressure region. With respect to the most accuratedata of Ref. 19, our estimate of the relative density ρ/ρ0at the compressibility peak is ∼ 3 − 4% larger, stillwithin the error bars. For larger pressures, we predict32.5 3.0 3.5 4.0 4.5 5.0/ 0020406080100120140160Pressure (GPa)PBE-FTVdW-DF1 (Ref.19)ExperimentsVMCLRDMCVMC (Ref. 30)RMC (Ref. 30)(a)1520253035Shock velocity Us (km/s) 10 15 20 25Particle velocity up (km/s)0.00.51.01.5Us (km/s)(b)FIG. 1: (1a) Principal Hugoniot in the density-pressure space. Red and yellow circles are the results obtained with ourMLPs trained on VMC and LRDMC datapoints, respectively, and a PBE baseline. Blue and pink triangles are the PBEresult calculated in this work and the VdW-DF1 result of Ref. 19 respectively. CEIMC results of Ref. 30 based on VariationalMonte Carlo (VMC) and Reptation Monte Carlo (RMC) are reported in green squares. Cyan diamonds are the experimentalresults of Refs. 16, 19, and 20. Dashed lines are guides for the eye. (1b) [top panel] Hugoniot in the up–Us space.Black-dashed line is the re-analyzed gas-gun fit reported in Ref. 19, i.e., the shock velocity extrapolated from measures ofmolecular deuterium at low pressure [13]. [bottom panel] Relative shock velocity with respect to the gas-gun fit. Only theexperimental points of Ref. 19 are reported.a Hugoniot mostly compatible within the experimentsbut systematically more compressible. However, in thisregime the correspondingly larger uncertainties in themeasures prevent a clear-cut assessment of our outcome.Our results also agree with the low-temperature CEIMCones reported in Ref. 30 within the statistical accuracy.Fig. 1b displays the same points in the up − Us space,where up is the particle velocity and Us is the shockvelocity, the two being calculated using the following RHrelations:up =√(p+ p0)(ρ−10 − ρ−1),Us = ρ−10√p+ p0ρ−10 − ρ−1.The difference ∆Us between these points and the linearfit on the gas-gun data re-analyzed in Ref. 19 is alsoshown (bottom panel of Fig. 1b). Notice that the dropin the slope of Us relative to up coincides with theonset of the molecular-atomic (MA) transition, while themagnitude of the ∆Us minimum relates to the positionof the relative compression peak. In particular, thePBE Hugoniot curve manifests a premature start of thedissociation, while it predicts correctly the magnitudeof the compressibility maximum. Our QMC resultscorrectly predict the position of the peak and startingslope, while showing some discrepancies for up ≳ 10 km/swith respect to the data of Ref. 19. In this regime, DFT,and in particular the result obtained with the VdW-DF1functional [54, 55], seems to be in better agreement withexperiments, thanks to a favorable error cancellation inthe Hugoniot [31]. We noticed how, in this study, thediscrepancy with the experiments is much milder thanthe value reported by previous QMC calculations atdensities and pressures close to the compressibility peak[29]. This can be due to the explicit optimization of theWF nodal surface provided by our WF ansatz, whichreduces the fixed node error mentioned in Ref. 31, theonly approximation left in any projective Monte Carlocalculation, such as LRDMC and RMC. The differencebetween the various methods is also apparent in theirequations of state, reported in the SM [50].The presence of an MA transition is investigatedin Fig. 2, where we report the radial distributionfunction, g(r), calculated on trajectories obtained withthe LRDMC model for several temperatures at densitiesclose to the Hugoniot curve. The inset of Fig. 2 displaysthe value of the molecular fraction m, defined as thepercentage of atoms that stay within a distance of 2 Bohr(roughly corresponding to the first g(r) minimum afterthe molecular peak) from another particle for longer thana characteristic time, here set to 6 fs. The results indicate4a distinct atomic character for T ≥ 10 kK and a clearmolecular peak at lower temperatures. The LRDMCmodel shows slower decay of the molecular fraction withtemperature than the PBE and VdW-DF1 ones, beingcompatible with the latter for temperatures above 10 kK.0 1 2 3 4 5 6 7 8 r (bohr)0.00.51.01.52.02.53.03.5 g(r)T = 4kKT = 6kKT = 7kKT = 8kKT = 10kKT = 15kKT = 20kKT = 35kK4 6 7 8 10 15Temperature (kK)0.20.40.60.8Molecular_fraction mPBEVdW-DF1DMC18 30 36 40 49 66Pressure (GPa)FIG. 2: g(r) for several temperatures and densities close tothe principal Hugoniot, obtained using the LRDMC model.The molecular fraction value, m, is reported in the inset foreach value of temperature up to 15kK. On the top axis thecorresponding pressure at the Hugoniot is also shown. Thevalues obtained with an ab initio DFT dynamics using thePBE and VdW-DF1 functionals are reported for comparisonat the same temperatures (notice that a pressure anddensity mismatch between methods can be present in thiscase due to different equations of state).IV. ERROR ANALYSISTo assess the quality of our principal Hugoniotdetermination, we analyzed the possible sources of errorsin relation to our machine learning scheme. There arethree main sources of errors: the uncertainties in thefit of H(ρ, T ), the prediction error of the MLP, and theuncertainties in the reference state energy estimate, i.e.e0 in Eq.(1). We verified that, in our case, the errorproduced by the fit is negligible compared to the othertwo sources, which we will discuss next.As mentioned before, we followed Ref. 34 to constructour MLPs and used a GKR model based on a modifiedversion of the SOAP kernel [33]. Our final dataset,including both training and test sets, comprises 561configurations selected through an iterative procedurewith 128 hydrogen atoms each, where we calculatedenergies, pressures and forces at the VMC and LRDMClevels. These configurations correspond to temperaturesfrom 4 kK up to 20 kK and Wigner-Seitz radii from1.80 Bohr to 2.12 Bohr. Finite size corrections (FSC)have also been estimated using the KZK functional [56].Details on the training set construction and the QMCcalculations, together with the performances of all MLPmodels can be found in the SM [50]. In particular wefound a final root mean square error (RMSE), calculatedon the test set, of the order of 7 meV/atom for theenergy, 250 meV/Å for the forces, and 0.75 GPa forthe pressures.At this point, it is worth to highlight some favourablefeatures of our machine learning approach, especially inapplications where it is coupled with computationallyexpensive methods such as QMC. They can be itemizedas follows:• transferability: the total energy of the system isexpressed as a sum of local terms [32], therefore ourmodels are capable of making accurate predictionson configurations whose size has never beenencountered in the training set. In particular, ourMLPs find their applicability to systems with anarbitrary number of atoms N .• efficiency and accuracy: within the ∆-learningframework, the machine learning task becomeseasier. Indeed, we obtained very accurate QMCpotentials, by training models on small datasetsand, thus, by reducing the amount of calculationsneeded. Moreover, since the computational costof the ML inference is negligible compared to thebaseline DFT calculation, we were able to performQMC-driven MD simulations at the cost of a DFTdynamics.• overfitting prevention: using a local sparsificationtechnique based on the farthest point sampling(see SM of Ref. 34), we discarded from eachconfiguration a possibly large fraction of thecorresponding N local environments, preventingoverfitting and allowing for an increased predictivepower of the model on unseen data. Sincethe computational cost of the predictions scaleswith the size of the training set, this proceduredrastically improves the efficiency of the finalmodel.We further validated the accuracy of our MLPsby comparing the Hugoniot curve obtained using twopotentials, independently trained with the same target,e.g. VMC, but with two different baselines. Taking intoaccount these results and the RMSE of the models, wecan estimate an uncertainty on the prediction of 0.06 onthe relative density ρ/ρ0 and 1 GPa on pressure.We now turn to the last source of error we identified,i.e. the one related to the calculation of e0 and p0.To estimate the reference state energy and pressure, wefollowed a procedure similar to Ref. 30. We performed a5path integral molecular dynamics (PIMD) simulation [57]on a system of N = 64 deuterium atoms at a temperatureT = 22 K and density ρ0 = 0.167 g/cm3 (correspondingto the initial conditions reported in Ref. 19), using DFT-PBE energy and forces. Details of this simulation arereported in the SM [50]. From the PIMD trajectory, weextracted 170 configurations and we calculated energiesand pressures with both DFT-PBE and QMC at VMCand LRDMC levels, adding the necessary finite sizecorrections. The reference sample was generated byextracting atomic positions from one of the 128 beadstaken at random, belonging to de-correlated snapshotsof the trajectory. Results for e0 for the various methodsare reported in Tab. I. The reference state pressure p0 isnot reported, since it is two orders of magnitude smallerthan the shocked pressure, and thus irrelevant for theHugoniot determination. Also in this case, we studiedthe effect of varying e0 within its confidence interval onthe Hugoniot density and pressure. In doing so, we alsotook into account the possible uncertainty on the energydifference e(ρ, T ) − e0 originating by the finite batchsize we used for estimating energy gradients in the WFoptimization. We estimated this uncertainty by runningoptimizations of increasing batch size on three different128-atom configurations. The results indicate an error≲ 0.5 mHa/atom on e(ρ, T )− e0. Taking everything intoaccount, varying the energy within standard deviationleads to shifts in the final principal Hugoniot whichstill fall in the error bars of our predictions estimatedpreviously.To summarize, we estimated the MLP prediction errorto be the most relevant source of uncertainty for theHugoniot, yielding, as discussed before, an absolute errorof 0.06 and 1 GPa on the relative density and pressure,respectively, reflected on the error bars reported in Fig. 1.V. CONCLUSIONSIn conclusion, using our recently proposed workflowfor the construction of MLPs, we have been able torun reliable VMC- and LRDMC-based MD simulationsand study the principal deuterium Hugoniot, in apressure range relevant for experiments. The accuracyepot (Ha/atom) e0 (Ha/atom)PBE -0.58217(2) -0.58055(2)VMC -0.58622(2) -0.58460(2)LRDMC -0.58660(2) -0.58498(2)VMC + FSC -0.58503(2) -0.58342(2)LRDMC + FSC -0.58542(2) -0.58380(2)TABLE I: Estimated potential (epot) and total (e0)energies per atom of the reference state at ρ0 = 0.167g/cm3 and T = 22 K for different methods, with andwithout finite size corrections.of the MLPs employed here has been extensively tested,supporting the validity of our calculations and estimatingtheir uncertainty. The resulting Hugoniot curveshows generally good agreement with the most recentexperimental measures, especially in the low temperaturemolecular regime. Exploiting the ∆-learning framework,we have also been able to treat FT electrons effects ina QMC-MLP, and we have thus managed to performaccurate simulations at higher temperatures. For thesetemperatures, the results systematically show a morecompressible Hugoniot curve than experiments, althoughthe experimental error bars are large in this regime. Theaforementioned discrepancy is milder than previouslyreported QMC calculations [29, 30], and falls withinthe measures uncertainty. In particular, this suggeststhat the use of optimized and more refined WFs hasa key role for obtaining good results in high-pressurehydrogen. We thus believe that our results will be usefulfor both future experimental research and numericalinvestigations of the Hugoniot. The efficiency of ourcomputational approach could be further improved, e.g.,by using cheaper baseline potentials than DFT. Longersimulations and larger systems will then be at reach.Other many-body methods, even more expensive thanQMC, can also be used as targets for this type ofMLPs, since the required size of the dataset is atleast one order of magnitude smaller compared to otherML approaches. Finally, our MLPs, and in particularthose trained on LRDMC data points, are promising forexploring the hydrogen phase diagram by keeping a highlevel of accuracy across a wide range of thermodynamicconditions.DATA AVAILABILITYThe machine learning code used in this work isavailable upon request. Additional information, such asdatasets, models, and detailed results of the simulationsare available at https://github.com/giacomotenti/QMC_hugoniot.ACKNOWLEDGMENTSThe computations in this work have mainly beenperformed using the Fugaku supercomputer providedby RIKEN through the HPCI System ResearchProject (Project ID: hp210038 and hp220060) andMarconi100 provided by CINECA through the ISCRAproject No. HP10BGJH1X and the SISSA three-year agreement 2022. A.T. acknowledges financialsupport from the MIUR Progetti di Ricerca diRilevante Interesse Nazionale (PRIN) Bando 2017 -grant 2017BZPKSZ. K.N. acknowledges financial supportfrom the JSPS Overseas Research Fellowships, fromGrant-in-Aid for Early Career Scientists (Grant No.6JP21K17752), from Grant-in-Aid for Scientific Research(Grant No. JP21K03400), and from MEXT LeadingInitiative for Excellent Young Researchers (Grant No.JPMXS0320220025). This work is also supported by theEuropean Centre of Excellence in Exascale ComputingTREX - Targeting Real Chemical Accuracy at theExascale. This project has received funding fromthe European Union’s Horizon 2020 - Research andInnovation program - under grant agreement no. 952165.The authors thank Dr. Guglielmo Mazzola and CesareCozza for the helpful suggestions and discussion. Wededicate this paper to the memory of Prof. SandroSorella (SISSA), who tragically passed away during thisproject, remembering him as one of the most influentialcontributors to the quantum Monte Carlo community,and in particular for deeply inspiring this work with thedevelopment of the ab initio QMC code, TurboRVB.[1] D. Saumon, G. Chabrier, and H. M. van Horn, Theastrophysical journal supplement series 99, 713 (1995).[2] J. J. Fortney and N. Nettelmann, Space Science Reviews152, 423 (2009).[3] Y. Miguel, T. Guillot, and L. Fayon, Astronomy &Astrophysics 596, A114 (2016).[4] S. Hu, V. Goncharov, T. Boehly, R. McCrory, S. Skupsky,L. A. Collins, J. D. Kress, and B. Militzer, Physics ofPlasmas 22, 056304 (2015).[5] A. P. Drozdov, M. I. Eremets, I. A. Troyan,V. Ksenofontov, and S. I. Shylin, Nature 525, 73 (2015).[6] M. Somayazulu, M. Ahart, A. K. Mishra, Z. M. Geballe,M. Baldini, Y. Meng, V. V. Struzhkin, and R. J. Hemley,Phys. Rev. Lett. 122, 027001 (2019).[7] J. M. McMahon, M. A. Morales, C. Pierleoni, and D. M.Ceperley, Rev. Mod. Phys. 84, 1607 (2012).[8] B. Cheng, G. Mazzola, C. J. Pickard, and M. Ceriotti,Nature 585, 217 (2020).[9] V. V. Karasiev, J. Hinz, S. X. Hu, and S. B. Trickey,Nature 600, E12 (2021).[10] B. Cheng, G. Mazzola, C. J. Pickard, and M. Ceriotti,Nature 600, E15 (2021).[11] W. J. Nellis, Reports on Progress in Physics 69, 1479(2006).[12] G. E. Duvall and R. A. Graham, Rev. Mod. Phys. 49,523 (1977).[13] W. J. Nellis, A. C. Mitchell, M. van Thiel,G. J. Devine, R. J. Trainor, and N. Brown,The Journal of Chemical Physics 79, 1480 (1983),https://doi.org/10.1063/1.445938.[14] M. D. Knudson, D. L. Hanson, J. E. Bailey, C. A. Hall,J. R. Asay, and W. W. Anderson, Phys. Rev. Lett. 87,225501 (2001).[15] G. V. Boriskov, A. I. Bykov, R. I. Il’kaev, V. D. Selemir,G. V. Simakov, R. F. Trunin, V. D. Urlin, A. N. Shuikin,and W. J. Nellis, Phys. Rev. B 71, 092104 (2005).[16] M. D. Knudson, D. L. Hanson, J. E. Bailey, C. A. Hall,J. R. Asay, and C. Deeney, Phys. Rev. B 69, 144209(2004).[17] D. G. Hicks, T. R. Boehly, P. M. Celliers, J. H. Eggert,S. J. Moon, D. D. Meyerhofer, and G. W. Collins, Phys.Rev. B 79, 014112 (2009).[18] P. Loubeyre, S. Brygoo, J. Eggert, P. M. Celliers, D. K.Spaulding, J. R. Rygg, T. R. Boehly, G. W. Collins, andR. Jeanloz, Phys. Rev. B 86, 144115 (2012).[19] M. D. Knudson and M. P. Desjarlais, Physical ReviewLetters 118, 1 (2017).[20] A. Fernandez-Pañella, M. Millot, D. E. Fratanduono,M. P. Desjarlais, S. Hamel, M. C. Marshall, D. J. Erskine,P. A. Sterne, S. Haan, T. R. Boehly, G. W. Collins, J. H.Eggert, and P. M. Celliers, Phys. Rev. Lett. 122, 255702(2019).[21] M. D. Knudson and M. P. Desjarlais, Journal of AppliedPhysics 129, 10.1063/5.0050878 (2021).[22] T. J. Lenosky, S. R. Bickham, J. D. Kress, and L. A.Collins, Phys. Rev. B 61, 1 (2000).[23] G. Galli, R. Q. Hood, A. U. Hazi, and F. m. c. Gygi,Phys. Rev. B 61, 909 (2000).[24] S. Bagnier, P. Blottiau, and J. Clérouin, Phys. Rev. E63, 015301 (2000).[25] S. A. Bonev, B. Militzer, and G. Galli, Phys. Rev. B 69,014101 (2004).[26] B. Holst, R. Redmer, and M. P. Desjarlais, Phys. Rev. B77, 184201 (2008).[27] L. Caillabet, S. Mazevet, and P. Loubeyre, Phys. Rev. B83, 094101 (2011).[28] V. V. Karasiev, S. X. Hu, M. Zaghoo, and T. R. Boehly,Physical Review B 99, 1 (2019).[29] N. M. Tubman, E. Liberatore, C. Pierleoni,M. Holzmann, and D. M. Ceperley, Physical ReviewLetters 115, 1 (2015).[30] M. Ruggeri, M. Holzmann, D. M. Ceperley, andC. Pierleoni, Physical Review B 102, 144108 (2020),2008.00269.[31] R. C. Clay, M. P. Desjarlais, and L. Shulenburger,Physical Review B 100, 75103 (2019).[32] J. Behler and M. Parrinello, Physical Review Letters 98,1 (2007).[33] S. De, A. P. Bartók, G. Csányi, and M. Ceriotti, Phys.Chem. Chem. Phys. 18, 13754 (2016).[34] A. Tirelli, G. Tenti, K. Nakano, and S. Sorella, Phys.Rev. B 106, L041105 (2022).[35] M. Casula, C. Filippi, and S. Sorella, Phys. Rev. Lett.95, 100201 (2005).[36] K. Nakano, R. Maezono, and S. Sorella, Phys. Rev. B101, 155106 (2020).[37] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048(1981).[38] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. 77, 3865 (1996).[39] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra,R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti,M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli,S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann,C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini,A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo,G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari,and R. M. Wentzcovitch, Journal of Physics: CondensedMatter 21, 395502 (2009).7[40] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau,M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni,D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo,A. D. Corso, S. de Gironcoli, P. Delugas, R. A.DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo,R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia,M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli,M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L.Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto,S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf,A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser,P. Umari, N. Vast, X. Wu, and S. Baroni, Journal ofPhysics: Condensed Matter 29, 465901 (2017).[41] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car,I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas,F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov,A. Urru, and S. Baroni, The Journal of Chemical Physics152, 154105 (2020), https://doi.org/10.1063/5.0005082.[42] H.pbek-jpaw psl.1.0.0.UPF pseudopotential availableat http://pseudopotentials.quantum-espresso.org/legacy_tables/ps-library/h.[43] A. Ricci and G. Ciccotti, Molecular Physics - MOL PHYS101, 1927 (2003).[44] C. Attaccalite and S. Sorella, Phys. Rev. Lett. 100,114501 (2008).[45] K. Nakano, C. Attaccalite, M. Barborini, L. Capriotti,M. Casula, E. Coccia, M. Dagrada, C. Genovese, Y. Luo,G. Mazzola, A. Zen, and S. Sorella, J. Chem. Phys. 152,204121 (2020).[46] K. Nakano, O. Kohulák, A. Raghav, M. Casula, andS. Sorella, J. Chem. Phys. 159, 224801 (2023).[47] J. Tiihonen, R. C. Clay III, and J. T. Krogel, J. Chem.Phys. 154, 204111 (2021).[48] K. Nakano, A. Raghav, and S. Sorella, The Journal ofChemical Physics 156, 034101 (2022).[49] K. Nakano, M. Casula, and G. Tenti, (2023),arXiv:2312.17608.[50] See Supplemental Material at [URL will be insertedby publisher ] for additional information about thecomputational details of QMC calculations, theMLP training and validation, the reference statecalculations, finite-size corrections, finite temperatureDFT simulations, and comparison with previous results[28, 30, 31, 34–37, 39–41, 44, 47–49, 56–69].[51] V. V. Karasiev, T. Sjostrom, J. Dufty, and S. B. Trickey,Phys. Rev. Lett. 112, 076403 (2014).[52] V. V. Karasiev, J. W. Dufty, and S. B. Trickey, Phys.Rev. Lett. 120, 076401 (2018).[53] C. Ben Mahmoud, F. Grasselli, and M. Ceriotti, Phys.Rev. B 106, L121116 (2022).[54] M. Dion, H. Rydberg, E. Schröder, D. C. Langreth, andB. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).[55] K. Berland, V. R. Cooper, K. Lee, E. Schröder,T. Thonhauser, P. Hyldgaard, and B. I. Lundqvist,Reports on Progress in Physics 78, 066501 (2015).[56] H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett.100, 126404 (2008).[57] F. Mouhat, S. Sorella, R. Vuilleumier, A. M. Saitta,and M. Casula, Journal of Chemical Theory andComputation 13, 2400 (2017).[58] M. Casula and S. Sorella, J. Chem. Phys. 119, 6500(2003).[59] S. Sorella, N. Devaux, M. Dagrada, G. Mazzola, andM. Casula, J. Chem. Phys. 143 (2015).[60] F. Becca and S. Sorella, Quantum Monte Carloapproaches for correlated systems (Cambridge UniversityPress, 2017).[61] K. Nakano, T. Morresi, M. Casula, R. Maezono, andS. Sorella, Phys. Rev. B 103, L121110 (2021).[62] C. J. Umrigar, Int. J. Quantum Chem 36, 217 (1989).[63] S. Sorella and L. Capriotti, J. Chem. Phys. 133, 234111(2010).[64] C. Filippi, R. Assaraf, and S. Moroni, J. Chem. Phys.144, 194105 (2016).[65] J. van Rhijn, C. Filippi, S. De Palo, and S. Moroni,Journal of chemical theory and computation 18, 118(2021).[66] S. Pathak and L. K. Wagner, AIP Adv. 10, 085213(2020).[67] P. Reynolds, R. Barnett, B. Hammond, R. Grimes, andW. Lester Jr, Int. J. Quantum Chem. 29, 589 (1986).[68] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785(1988).[69] N. D. Mermin, Phys. Rev. 137, A1441 (1965).Supplemental material: Principal deuterium Hugoniot via Quantum Monte Carlo and ∆-LearningGiacomo Tenti,1, ∗ Kousuke Nakano,2, † Andrea Tirelli,3 Sandro Sorella,1 and Michele Casula41International School for Advanced Studies (SISSA), Via Bonomea 265, 34136 Trieste, Italy2Center for Basic Research on Materials, National Institute for Materials Science (NIMS), Tsukuba, Ibaraki 305-0047, Japan3International School for Advanced Studies (SISSA)Via Bonomea 265, 34136 Trieste, Italy4Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC),Sorbonne Université, CNRS UMR 7590, MNHN, 4 Place Jussieu, 75252 Paris, France(Dated: May 7, 2024)I. COMPUTATIONAL DETAILS OF QMC CALCULATIONSThe Variational Monte Carlo (VMC) and lattice regularized diffusion Monte Carlo (LRDMC) [1] calculations in this studywere performed by the TurboRVB package [2]. The package employs a many-body WF ansatz Ψ which can be written asthe product of two terms, i.e., Ψ = ΦAS × exp J , where the term exp J and ΦAS are conventionally called Jastrow and anti-symmetric parts, respectively. The antisymmetric part is expressed as an Antisymmetrized Geminal Power (AGP) that reads:ΦAGP (r1, . . . , rN) = Â[Φ(r↑1, r↓1)Φ(r↑2, r↓2)· · ·Φ(r↑N/2, r↓N/2)],where Â is the antisymmetrization operator, andΦ(r↑, r↓)is calledthe paring (geminal) function [3]. The spatial part of the geminal function is expanded over the Gaussian-type atomic orbitals(GTOs): Φ(ri, r j)=∑a,b,l,mf{a,l},{b,m}ψa,l (ri)ψb,m(r j) where ψa,l and ψb,m are primitive Gaussian atomic orbitals, their indices land m indicate different orbitals centered on atoms a and b, i and j are coordinates of spin up and down electrons, respectively,and f{a,l},{b,m} are the variational parameters. The pairing function can be also written as Φ(ri, r j)=∑Mk=1 λkϕk(ri)ϕk(r j) withλk > 0, where ϕk(r) is a molecular orbital, i.e., ϕk(r) =∑Li=1 ci,kψi(r) and L is the total number of atomic orbitals. When theparing function is expanded over M molecular orbitals where M is equal to half of the total number of electrons (N/2), the AGPcoincides with the Slater-Determinant ansatz. In this work, we took M = N/2.The Jastrow term is composed of one-body, two-body and three/four-body factors (J = J1 + J2 + J3/4). The one-bodyand two-body factors are essentially used to fulfill the electron-ion and electron-electron cusp conditions, respectively, and thethree/four-body factor is employed to consider further electron-electron correlations (e.g., electron-nucleus-electron). The one-body Jastrow is decomposed into the so-called homogeneous and inhomogeneous parts, i.e., J1 = Jhom1 + Jinh1 . The homogeneousone-body Jastrow factor is J1hom (r1, . . . , rN) =∑i,I(−(2ZI)3/4u((2ZI)1/4 |ri − RI |))where ri are the electron positions, RI are theatomic positions with corresponding atomic number ZI , and u (r) is a short-range function containing a variational parameter b:u (r) = b2(1 − e−r/b). The inhomogeneous one-body Jastrow factor Jinh1 is represented as:Jinh1 (r1, . . . , rN) =∑Ni=1∑Natoma=1(∑lMa,lχa,l (ri)), where ri are the electron positions, Ra are the atomic positions with cor-responding atomic number Za, l runs over atomic orbitals χa,l (e.g., GTO) centered on the atom a, Natom is the total num-ber of atoms in a system, and {Ma,l} are variational parameters. The two-body Jastrow factor is defined as: J2 (r1, . . . rN) =exp(∑i< jv(∣∣∣ri − r j∣∣∣)), where v (r) = 12 r · (1 − F · r)−1 and F is a variational parameter. The three-body Jastrow factor is:J3/4 (r1, . . . rN) = exp(∑i< jΦJas(ri, r j)), and ΦJas(ri, r j)=∑l,m,a,bga,l,b,mχa,l (ri) χb,m(r j), where, as before, the indices l and mindicate different orbitals centered on corresponding atoms a and b. Here, the coefficients of the three/four-body Jastrow factorwere set to zero for a , b; this significantly decreases the number of variational parameters and rarely affects variational ener-gies. In this study, we used a basis set of [4s2p1d] and [2s2p1d] GTOs for the antisymmetric part and Jastrow part, respectively.The exponents in the gaussian orbitals were optimized beforehand for our particular system and thermodynamic conditions. Weverified that the optimal exponents do not vary significantly in the density and temperature range explored and thus we fixedtheir value for all the configurations in the dataset. For the antisymmetric part, we further reduced the number of variationalparameters thanks to contracted orbitals. In particular, the WF was first initialized with the Density Functional theory (DFT)built-in code of the TurboRVB package at the Γ point and then projected into a basis of 6 hybrid orbitals (geminal embedded∗ gtenti@sissa.it† kousuke 1123@icloud.comarXiv:2301.03570v2 [cond-mat.str-el] 3 May 20242orbitals or GEOs) [4]. The resulting geminal WF is then expressed as Φ(ri, r j)=∑a,b,µ,νf{a,µ},{b,ν}ψ̃a,µ (ri) ψ̃b,ν(r j) were the sumnow runs over the GEOs ψ̃a,µ (in general, for a , b, ψa,µ , ψb,µ ). The WF in the new basis was initialized by running a DFT-LDAcalculation at the Baldereschi point k = (0.25, 0.25, 0.25) (in crystal units). Both the Jastrow factor and the antisymmetric partof the WF were optimized using the Hessian matrix method. Finally, only the variational parameters f{a,µ},{b,ν} correspondingto atoms closer than a cutoff rc = 4.0 Bohr were optimized. The details of the constrained optimization technique are writtenin Ref. 5. For stable optimizations, we projected the AGP into N/2 molecular orbitals at each step [5], such that the resultingoptimized WF is equivalent to an optimal JSD WF. We found that the explicit optimization of the determinantal part of the WF isnot only important for getting more accurate energies and forces, but also mitigates the self-consistent error [6–8] present in theJSD ansatz. Total energies and forces are calculated at the VMC and the LRDMC levels with the optimized wavefunctions. TheLRDMC calculations were performed by the original single-grid scheme [1] with the discretization grid size a = 0.20 Bohr. Theresulting WF is a clear improvement with respect to JSD ansatz commonly employed in QMC, where only the Jastrow factorparameters are optimized, and the SD is defined with frozen DFT orbitals. This can be appreciated by comparing the referencestate average potential energy based on our optimized workflow and the one calculated using a standard JSD WF (see Tab.I). Toobtain a statistically meaningful value of VMC and LRDMC forces with finite variance [9], the so-called reweighting techniquesare needed because the Hellmann–Feynman (HF) and Pulay terms may diverge when the minimum electron–nucleus distancevanishes and when an electronic configuration is close to the nodal surface, respectively [7]. The infinite variance of the firstterm is cured by applying the so-called space-warp coordinate transformation (SWCT) algorithm [7, 10–12], whereas that ofthe second term can be alleviated by modifying the VMC sampling distribution using a modified trial wave function that differsfrom the original trial wave function only in the vicinity of the nodal surface [13], which we dub the Attaccalite and Sorella(AS) regularization. The AS regularization is not an optimal regularization for this purpose because it enforces a finite densityof walkers on the nodal surface [14]. Therefore, in this study, we employed the regularization technique recently proposed byPathak and Wagner [15] combined with mixed-averaged forces proposed by Reynolds [16].The final statistical noise on energies, forces and pressures was of the order of 7 × 10−6 Ha/atom, 1 mHa/Bohr and 0.05 GPa,for VMC, and 2 × 10−5 Ha/atom, 2 mHa/Bohr and 0.1 GPa, for LRDMC.II. MLP TRAINING AND VALIDATIONA. Dataset constructionTo construct our dataset, we performed a first set of ab initio DFT MD simulations with the PBE functional [17] on a system ofN = 128 atoms for temperatures in the range [4kK, 20kK] and densities in the range [1.80 Bohr, 2.12 Bohr], from which we ex-tracted an initial set of decorrelated snapshots. We then added other configurations according to an active learning scheme: witha model trained using this first dataset we ran MD simulations and iteratively selected new points where the MLP performanceswere expected to be poor. In particular we did this by monitoring, for each unseen configuration, the quantityχ =1NN∑i=1minµ∈training setK(Ri,Rµ) (1)where K(Ri,Rµ) is the normalized SOAP kernel between the i-th local environment of the configuration Ri and the µ-thlocal environment in the training set Rµ. The number χ defined in (1) gives a quantitative measure of ”how far” the unknownconfiguration is from what is already included in the training set. In particular, we stopped adding configurations when χ did notdrop under a fixed threshold of 0.80 during the dynamics. The final dataset comprised 561 configurations of 128 atoms in total.The final range of temperatures and Wigner-Seitz radii rs spanned by these configurations was [4 kK : 20 kK ] and [1.80 Bohr :2.12 Bohr ], respectively. The distribution of these configurations in the rs − T space is shown in Fig. 1.JSD optimized JSDeV MC (Ha / atom ) -0.58465(3) -0.58622(2)eLRDMC (Ha / atom) -0.58653(2) -0.58660(2)TABLE I: Average VMC and LRDMC energy for the reference state, calculated on a Jastrow-Slater ansatz with Jastrowoptimization and 4 × 4 × 4 k points integration and the optimized JSD WF obtained with our workflow at the Baldereschi point.31.80 1.84 1.88 1.92 1.96 2.00 2.04 2.08 2.12rs (Bohr)468101520Temperature (kK)19232503102624290200121417016013790140148177172141015222718350000028000002800000FIG. 1: Distribution in temperature and Wigner-Seitz radius of the configurations belonging to the final dataset used. The exactnumber of configurations is also shown.B. Training detailsFor the training procedure we followed the strategy outlined in [18, §I.B]. The energy is decomposed into local contributionsof the form:e(R, cµ) = e0 + e2body(R, {αi}) +Nenvs∑j=1β jK(R, R̄ j)where e2body is a pairwise term written as a linear combination of spline functions, K(R, R̄ j) is the SOAP kernel defined in [18]between the local environment to be predicted and one of the Nenvs environments belonging to training set, and (e0, {αi}, {β j}) ≡ cµare the variational parameter to optimize. The cost function C employed is the regularized weighted sum of the MSE on theobservables, i.e.C = C(cµ) = αMSE(E, Ê(cµ)) + βMSE(F, F̂(cµ)) + γMSE(P, P̂(cµ)) + λ||cµ||2,where E, F, P and Ê, F̂, P̂ are the vectors representing energy, forces and pressures obtained through QMC simulations andGKR, respectively. For the choice of the model hyperparameters, a cross-validation test led the following:• the cutoff radius used to compute local environments has been set to rc = 4.0 Bohr.• the term e2body was expressed as a linear combination of cubic splines defined on a grid of 10 equally spaced points,extending from rmin = 0.3 Bohr and rmax = 5.5 Bohr.4• the parameters α, β, γ and λ determining the cost function C(cµ) have been set to 20, 10, 107, 10−5 and and 20, 10, 2 ×107, 10−5 (in the appropriate inverse atomic units) for the VMC and LRDMC models respectively.• The number of environments selected among the full 128 atoms configurations was Nenvs = 6000.C. Models performanceThe performances of the models employed are measured through the root mean squared error (RMSE) on the observables onwhich the models were trained, avaluated on the test set. Such RMSEs are reported in Table II.RMSEE (mHa / atom) RMSE f (mHa / Bohr) RMSEp (GPa)VMC - PBE 0.26 4.9 0.76VMC - PBE (with FSC) 0.29 5.1 1.1VMC - LDA 0.35 5.5 0.80LRDMC - PBE (with FSC) 0.36 5.5 0.72TABLE II: Value of the RMSE on different observables as calculated on the test set for the final models used in the simulations.D. Effect of different baselinesIn order to further validate the accuracy of an MLP, a common strategy is to compare the results of the dynamics obtainedusing the trained model with those obtained with the target ab initio method directly, at least for some small system sizes. Inour case this is not an easy task, given the large computational time that would be needed for computing energies and forcesat each step with QMC. An alternative way to establish the performances of the models is to look at the variance of the resultsobtained with MLPs trained using different baselines. The Hugoniot function H(ρ,T ) and pressure at T = 8 kK, calculated usingmodels trained on the difference between VMC and two DFT functionals (i.e., PBE and LDA), are shown in Fig. 2. We observethat the variability of the Hugoniot relative density and pressure is compatible with the value estimated from the RMSE (i.e.,∆ρ/ρ0 = 0.06 and ∆p = 1 GPa).1.900 1.925 1.950 1.975rs (Bohr)0.0050.0000.0050.010H(r s,T) (Ha)1.900 1.925 1.950 1.975rs (Bohr)3436384042Pressure( GPa )T = 8000KVMC-PBEVMC-LDAFIG. 2: Results for the Hugoniot function H(rs,T ) and pressure for T = 8 kK with different MLPs trained on VMC data, usingdifferent baselines potentials.550 100 150 200 250# beads0.020.040.060.080.100.12Kinetic energy (Ha) virial PBEprimitive PBEvirial BLYPprimitive BLYPFIG. 3: Convergence of virial and primitive estimators for the quantum kinetic energy, as computed with Eqs. (2) (3), with thenumber of replicas used in the PIMD simulation.III. REFERENCE STATE CALCULATIONSAs explained in the main text, a crucial part in the numerical determination of the Hugoniot is to estimate the reference stateenergy per atom e0 and pressure p0. In particular, having a precise value of e0 within the target method is important to takeadvantage of possible error cancellation effects and remove biases related to finite basis sets. We considered a system of N = 64deuterium atoms at T0 = 22 K and ρ0 = 0.167 g/cm−3 and ran a path integral Ornstein-Uhlenbeck molecular dynamics [19](PIOUD) simulation to account for quantum effects, which are required because of the light deuterium mass and low temperature.Forces and energy were calculated with DFT through the Quantum Espresso package[20–22]. We checked the dependence ofthermodynamic quantities on the number of replicas (or beads) M and on the choice of the DFT functional by studying thequantum kinetic energy T for several values of M using the BLYP [23] and PBE functionals. In particular we considered twoestimators for T , namely the virial and primitive (or Barker) estimator, given respectively byTM,vir =N2β+12M3N∑i=1M∑j=1(x( j)i − x̄i)∂x( j)iV (2)TM,pri =3NM2β− mM2β2ℏ2M∑j=1(x( j)i − x( j−1)i)2(3)where M is the number of replicas used in the PIOUD simulation, x( j) =(x( j)1 , . . . , x( j)3N)are the coordinates of the systembelonging to the j-th bead, x̄i =1M∑Mj=1 x( j)i is the centroid position and β = kBT0. The results are shown in Fig. 3.We noticed that a very large number of replicas is necessary for having a sufficiently converged result, while the value obtainedwith the two functionals is extremely similar for all values of M. At the end we chose to use the PBE functional and M = 128replicas to have a reasonable trade-off between convergence and computational cost. For the DFT calculation we used a 60 Ryplane waves cutoff and a 2 × 2 × 2 Monkhorst-Pack k point mesh; for the dynamics we used a time step of 0.3 fm and let thesystem thermalize for 0.3 ps. We then extracted one configuration from a randomly chosen bead every 10 MD steps, for a totalof Nsample = 170 snapshots. Finally, the potential energy of these configurations was calculated using the appropriate method(PBE, VMC or LRDMC). We then estimated e0 for each method as6e0 =1N1Nsample∑sampleEpot (xi) + T PBE256,pri (4)using the value of the primitive estimator at M = 256 beads as the best guess for the converged value of the kinetic energy.The approximation for the potential energy was checked by running PBE simulations on this set and confirming that the ”true”mean value (as calculated by averaging over the beads and the trajectory) was consistent with our estimate obtained by averagingover the sample.IV. FINITE SIZE CORRECTIONSIn this section we investigate the effect of finite size corrections (as estimated using the KZK functional [24]) on our results.In Fig. 4 we show the Hugoniot function ( at T = 4 kK, T = 8 kK and T = 35 kK) given by two models trained on VMC andVMC with finite size corrections respectively, both with a PBE baseline potential. The difference between the two turns out tobe similar to the prediction error discussed in Sec. II D, for the system size used in the simulations (i.e., N = 128). At the endwe chose to apply finite size correction only for the model trained with LRDMC data.2.12 2.16 2.20rs (Bohr)0.0030.0020.0010.0000.0010.0020.003H(r s,T) (Ha)T = 4000KVMC (no FSC)VMC (with FSC)1.90 1.94 1.98rs (Bohr)0.00500.00250.00000.00250.00500.00750.0100H(r s,T) (Ha)T = 8000KVMC (no FSC)VMC (with FSC)1.90 1.94 1.98rs (Bohr)0.030.020.010.000.010.02H(r s,T) (Ha)T = 35000KVMC (no FSC)VMC (with FSC)FIG. 4: Hugoniot functions obtained using two MLPs trained on the difference between PBE and VMC with and without finitesize corrections respectively, for T = 4 kK, T = 8 kK and T = 35 kK.V. FINITE TEMPERATURE DFT SIMULATIONSUsing Mermin’s extension of the Hohenberg and Kohn theorems to non-zero temperature [25] we can treat finite temperatureelectrons in DFT by appropriately occupying the bands of the system according to the Fermi-Dirac distribution and minimizingthe Helmholtz free energy functional A = E − TS . In this work we performed finite temperature DFT (FT-DFT) simulationsto obtain the PBE Hugoniot and estimating the effect on the QMC Hugoniot. In both cases we used the zero temperaturePBE functional for the simulations. Even if this is not rigorous, recent FT-DFT results on the Hugoniot using a temperaturedependent GGA functional [26] have shown that for T ≲ 40 kK this approximation provides consistent results. In Fig. 5we show the convergence of the free energy and some force components with the number of bands calculated for the FT-PBEfunctional at two values of temperature. In the simulations we decided to use 120 bands for T = 10 kK, T = 15 kK and 150bands for T = 20 kK and T = 35 kK.VI. COMPARISON BETWEEN EQUATIONS OF STATEThe equations of state at T = 8 kK reported in Ref. [27] for both variational and reptation Monte Carlo are shown in Fig. 6,together with the VMC-MLP, LRDMC-MLP (with PBE baseline) and the ab initio PBE and vdW-DF1 ones. From this figurewe can observe a ∼ 4 GPa discrepancy between the pressure estimated with our MLPs and the one in Ref. [27], which is one7100 150 200nbands15.08015.07515.07015.065Free energy (eV / atom)100 150 200nbands0.050.000.050.10f i(ev / Å)T = 15000 K (a)100 150 200nbands16.116.015.915.8Free energy (eV / atom)100 150 200nbands0.40.20.0f i(ev / Å)T = 35000 K (b)FIG. 5: Convergence of the free energy A = E − TS and some force components vs the number of bands calculated withDFT-PBE for T = 15 kK (5a), and T = 35 kK (5b). For the force components the difference ∆ fi between the force obtainedwith N = nbands bands and N = 200 is shown. Bands are occupied using the Fermi-Dirac distribution at the appropriatetemperature.of the causes of the slightly different Hugoniot position. In fact, by shifting the data reported in [27] by this quantity, we cansee an almost perfect match of the Hugoniot positions (Fig. 7). In our case, the VMC and LRDMC pressures, which we usedfor training, were calculated using the adjoint algorithmic differentiation method to obtain directly the derivative of the totalenergy with respect to the cell parameters. As mentioned in the main text, the self consistency error can, in principle, spoilthe consistency between the pressure and the potential energy surface of the system. We mitigated this by directly optimizingthe WF nodal surface, as previously described. Nevertheless, we identified the pressure as a particularly sensitive quantity,difficult to reliably obtain, and we took particular care when including it in the training. In Ref. [27] a virial estimator wasused, which can in principle produce discrepancies of the order of magnitude observed here, as shown in Ref. [28]. We alsopoint out that both our models and the vdW-DF1 functional predict pressures slightly larger than those given by the PBE in theproximity of the Hugoniot position, in contrast with the EOS reported in [27], that predicts lower pressures. If it is true that8QMC pressures are consistently lower than PBE ones on a given configuration, here we found that the different dynamics of thesystem samples different points in the phase space, such that the resulting average pressure is, in some cases, even larger for ourMLPs. To demonstrate this, we extracted two samples of configurations from an ab initio PBE and a LRDMC-MLP trajectory,respectively, and then computed both PBE pressures and ∆ corrections. The results show how the increased PBE pressures forthe LRDMC sampled trajectory can compensate the negative contribution given by the MLP correction (Tab. III and Fig. 8).1.80 1.85 1.90 1.95 2.00rs (Bohr)30405060Pressure (GPa)PBEvdW-DF1LRDMCVMCVMC Ref. 27RMC Ref. 27FIG. 6: Average pressure vs rs for T = 8 kK. Results obtained with our MLPs are shown together with the ones reported in [27]and those calculated with the PBE and vdW-DF1 functionals.1.80 1.85 1.90 1.95 2.00rs (Bohr)0.030.020.010.000.01H(r s,T) (Ha)Ref. 27 (4 GPa shift)Ref. 27LRDMCFIG. 7: LRDMC (this work) and RMC (Ref.[27]) Hugoniot curves for T = 8 kK. For the latter, we also show the result ofapplying a costant 4 GPa positive pressure shift.9Virial pressure (GPa) PBE sample LRDMC samplePBE 14.6(4) 17.7(4)LRDMC 11.1(4) 14.5(5)∆ correction −3.57(4) −3.17(4)TABLE III: Average value of the virial pressure calculated by DFT-PBE and by the LRDMC MLP, on sampled trajectoriesobtained using the two methods, at thermodynamic conditions rs = 1.92 Bohr and T = 8 kK. The average ∆ correction is alsoreported.0 20 40# PBE samples5101520Pressure (GPa)0 20 40# LRDMC samples10.012.515.017.520.022.5Pressure (GPa)LRDMCPBEFIG. 8: Virial pressure calculated by DFT-PBE and by the LRDMC MLP, on sampled trajectories obtained using the twomethods, at thermodynamic conditions rs = 1.92 Bohr and T = 8 kK.[1] M. Casula, C. Filippi, and S. Sorella, Diffusion monte carlo method with lattice regularization, Phys. Rev. Lett. 95, 100201 (2005).[2] K. Nakano, C. Attaccalite, M. Barborini, L. Capriotti, M. Casula, E. Coccia, M. Dagrada, C. Genovese, Y. Luo, G. Mazzola, A. Zen,and S. Sorella, Turborvb: A many-body toolkit for ab initio electronic simulations by quantum monte carlo, J. Chem. Phys. 152, 204121(2020).[3] M. Casula and S. Sorella, Geminal wave functions with jastrow correlation: A first application to atoms, J. Chem. Phys. 119, 6500 (2003).[4] S. Sorella, N. Devaux, M. Dagrada, G. Mazzola, and M. Casula, Geminal embedding scheme for optimal atomic basis set construction incorrelated calculations, J. Chem. Phys. 143 (2015).[5] F. Becca and S. Sorella, Quantum Monte Carlo approaches for correlated systems (Cambridge University Press, 2017).[6] J. Tiihonen, R. C. Clay III, and J. T. Krogel, Toward quantum monte carlo forces on heavier ions: Scaling properties, J. Chem. Phys. 154,204111 (2021).[7] K. Nakano, A. Raghav, and S. Sorella, Space-warp coordinate transformation for efficient ionic force calculations in quantum montecarlo, The Journal of Chemical Physics 156, 034101 (2022).[8] K. Nakano, M. Casula, and G. Tenti, Unbiased and affordable atomic forces in ab initio variational monte carlo, (2023),arXiv:2312.17608.[9] K. Nakano, T. Morresi, M. Casula, R. Maezono, and S. Sorella, Atomic forces by quantum monte carlo: Application to phonon dispersioncalculations, Phys. Rev. B 103, L121110 (2021).[10] C. J. Umrigar, Two aspects of quantum monte carlo: determination of accurate wavefunctions and determination of potential energysurfaces of molecules, Int. J. Quantum Chem 36, 217 (1989).10[11] S. Sorella and L. Capriotti, Algorithmic differentiation and the calculation of forces by quantum monte carlo, J. Chem. Phys. 133, 234111(2010).[12] C. Filippi, R. Assaraf, and S. Moroni, Simple formalism for efficient derivatives and multi-determinant expansions in quantum montecarlo, J. Chem. Phys. 144, 194105 (2016).[13] C. Attaccalite and S. Sorella, Stable liquid hydrogen at high pressure by a novel ab initio molecular-dynamics calculation, Phys. Rev.Lett. 100, 114501 (2008).[14] J. van Rhijn, C. Filippi, S. De Palo, and S. Moroni, Energy derivatives in real-space diffusion monte carlo, Journal of chemical theory andcomputation 18, 118 (2021).[15] S. Pathak and L. K. Wagner, A light weight regularization for wave function parameter gradients in quantum monte carlo, AIP Adv. 10,085213 (2020).[16] P. Reynolds, R. Barnett, B. Hammond, R. Grimes, and W. Lester Jr, Quantum chemistry by quantum monte carlo: Beyond ground-stateenergy calculations, Int. J. Quantum Chem. 29, 589 (1986).[17] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B 23,5048 (1981).[18] A. Tirelli, G. Tenti, K. Nakano, and S. Sorella, High-pressure hydrogen by machine learning and quantum monte carlo, Phys. Rev. B 106,L041105 (2022).[19] F. Mouhat, S. Sorella, R. Vuilleumier, A. M. Saitta, and M. Casula, Fully quantum description of the zundel ion: Combining variationalquantum monte carlo with path integral langevin dynamics, Journal of Chemical Theory and Computation 13, 2400 (2017).[20] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso,S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov,P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations ofmaterials, Journal of Physics: Condensed Matter 21, 395502 (2009).[21] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni,N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo,R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili,N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf,A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materialsmodelling with quantum ESPRESSO, Journal of Physics: Condensed Matter 29, 465901 (2017).[22] P. Giannozzi, O. Baseggio, P. Bonfà, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino,A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, Quantum espresso toward the exascale, The Journal of Chemical Physics 152,154105 (2020), https://doi.org/10.1063/5.0005082.[23] C. Lee, W. Yang, and R. G. Parr, Development of the colle-salvetti correlation-energy formula into a functional of the electron density,Phys. Rev. B 37, 785 (1988).[24] H. Kwee, S. Zhang, and H. Krakauer, Finite-size correction in many-body electronic structure calculations, Phys. Rev. Lett. 100, 126404(2008).[25] N. D. Mermin, Thermal properties of the inhomogeneous electron gas, Phys. Rev. 137, A1441 (1965).[26] V. V. Karasiev, S. X. Hu, M. Zaghoo, and T. R. Boehly, Exchange-correlation thermal effects in shocked deuterium: Softening theprincipal Hugoniot and thermophysical properties, Physical Review B 99, 1 (2019).[27] M. Ruggeri, M. Holzmann, D. M. Ceperley, and C. Pierleoni, Quantum Monte Carlo determination of the principal Hugoniot of deuterium,Physical Review B 102, 144108 (2020), 2008.00269.[28] R. C. Clay, M. P. Desjarlais, and L. Shulenburger, Deuterium Hugoniot: Pitfalls of thermodynamic sampling beyond density functionaltheory, Physical Review B 100, 75103 (2019).