# Fileset

[2502.02447v2.pdf](https://mdr.nims.go.jp/filesets/71a89a09-19a7-4d3b-b2b4-d9961e5e69db/download)

## Creator

Giacomo Tenti, Bastian Jäckl, [Kousuke Nakano](https://orcid.org/0000-0001-7756-4355), Matthias Rupp, Michele Casula

## Rights



## Other metadata

[Hydrogen liquid-liquid transition from first principles and machine learning](https://mdr.nims.go.jp/datasets/b3cac5fb-4d6a-47f2-8be4-9df0f77bf815)

## Fulltext

Hydrogen liquid-liquid transition from first principles and machine learningHydrogen liquid-liquid transition from first principles and machine learningGiacomo Tenti,1, ∗ Bastian Jäckl,2 Kousuke Nakano,3, 4 Matthias Rupp,2, 5 and Michele Casula6, †1International School for Advanced Studies (SISSA), Via Bonomea 265, 34136 Trieste, Italy2Department of Computer and Information Science,University of Konstanz, Box 188, 78457 Konstanz, Germany3Center for Basic Research on Materials, National Institute for MaterialsScience (NIMS), 1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, Japan4Center for Emergent Matter Science (CEMS), RIKEN,2-1 Hirosawa, Wako, Saitama, 351-0198, Japan.5Scientific Instrumentation and Process Technology unit,Luxembourg Institute of Science and Technology (LIST), Maison de l’Innovation,5 avenue des Hauts-Fourneaux, L-4362 Esch-sur-Alzette, Luxembourg6Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC),Sorbonne Université, CNRS UMR 7590, MNHN, 4 Place Jussieu, 75252 Paris, France(Dated: November 18, 2025)The molecular-to-atomic liquid-liquid transition (LLT) in high-pressure hydrogen is a fundamentaltopic touching domains from planetary science to materials modeling. Yet, the nature of theLLT is still under debate. To resolve it, numerical simulations must cover length and time scalesspanning several orders of magnitude. We overcome these size and time limitations by constructinga fast and accurate machine-learning interatomic potential (MLIP) built on the MACE neuralnetwork architecture. The MLIP is trained on Perdew-Burke-Ernzerhof (PBE) density functionalcalculations and uses a modified loss function correcting for an energy bias in the molecular phase.Classical and path-integral molecular dynamics driven by this MLIP show that the LLT is alwayssupercritical above the melting temperature. The position of the corresponding Widom line agreeswith previous ab initio PBE calculations, which in contrast predicted a first-order LLT. Accordingto our calculations, the crossover line becomes a first-order transition only inside the molecularcrystal region. These results call for a reconsideration of the LLT picture previously drawn.Introduction. Pristine hydrogen is one of the mostwidely studied materials, both theoretically andexperimentally [1]. Indeed, despite being made of thesimplest atomic constituent, it exhibits a surprisinglyrich phase diagram. The molecular-to-atomic transitionhappening upon liquid hydrogen compression is crucialfor planetary science, in particular for understandingthe interior of giant gas planets [2] and their magneticfields [3]. This liquid-liquid transition (LLT) has beenextensively investigated both experimentally [4–14] andvia numerical simulations [15–37].As is often the case for high-pressure hydrogen,experiments based on static [38, 39] and dynamic [40]compression give contrasting results, the dynamicexperiments predicting a larger transition pressure.Uncertainty is also present in the numerical simulations.Results obtained with ab initio molecular dynamics(AIMD) using density functional theory (DFT) showlarge variability with respect to the choice of theexchange-correlation functional. For instance, thetransition pressure can vary by 200 GPa when includinglong-range van der Waals corrections [41].The nature of the LLT has been debated, withmany first-principles simulations [15, 21, 22, 26, 27,30, 35] suggesting a first-order LLT below a critical∗ gtenti@sissa.it† michele.casula@impmc.upmc.frtemperature Tc and pressure pc, based on the observationof kinks in the equation of state (EOS). Given thelong correlations in both time and space expected nearthe transition, the outcome of first-principles moleculardynamics (MD) simulations using DFT or quantumMonte Carlo (QMC) methods has been questionedbecause of the short time scales and/or the small systemsizes considered. This could be reflected in the largevariability of the predicted Tc and pc, which range fromTc ∼ 4000K and pc ∼ 30GPa [27, 29], to Tc ∼ 1250Kand pc ∼ 150GPa for the most recent simulations [30].In the past two decades, large-scale simulations havebeen made possible by the introduction of machinelearning interatomic potentials [42] (MLIPs). Theseprovide results at much lower computational cost, oncetrained on datasets generated with ab initio methods.However, the derivation of an MLIP is a delicate step perse, possibly introducing a residual bias. A recent MLIPmodel, based on the NequIP [43] neural network (NN)and trained on DFT data, gave a Tc value smaller thanprevious estimates [44]. This model has a pressure biasof about 5 GPa in the EOS, once compared with thecorresponding ab initio predictions. Recent simulationsbased on another MLIP trained on the same DFTfunctional [34] suggested that the LLT is instead asmooth crossover, even though the accuracy of the modelhas been criticized [35, 45].In this work, we present results for the LLT obtainedwith MACE [46], a framework combining an equivariantmessage-passing NN with high-body-order messages,arXiv:2502.02447v2  [cond-mat.dis-nn]  16 Nov 2025mailto:gtenti@sissa.itmailto:michele.casula@impmc.upmc.frhttps://arxiv.org/abs/2502.02447v22which shows remarkable accuracy compared to othertypes of MLIPs. [47] In particular, we trained a MACEMLIP on the Perdew-Burke-Ernzerhof (PBE) functionaland used it to study the nature of the LLT as a function ofboth system size and simulation length. This also allowsus to directly compare our results with the large body ofliterature that used the same functional [1, 18, 19, 21–23, 25, 27, 29, 30, 33–35, 37, 44, 45, 47]. Our results reveala first-order transition between an atomic liquid and amolecular solid at low temperatures (T < 950K), due tothe melting line proximity. At higher temperatures, oursimulations indicate a crossover between molecular andatomic liquids in the thermodynamic limit.MACE model from first principles. To overcome thetime and size limitations of ab initio simulations, whilealso retaining a high accuracy across the LLT, weconstructed a MACE MLIP by applying a correctionto the loss function usually minimized to find optimalparameters of the model. During the optimization, theloss and its gradients are evaluated on a subset of trainingconfigurations (”batch”) of dimension Nb. For the µ-th configuration in the batch, having Nµ atoms, Epredµ ,Fpredµ,j , and σpredµ are the total energy, the force acting onthe j-th atom and the virial stress tensor, respectively,as predicted by the MLIP while Erefµ , Frefµ,j , and σrefµ arethe same quantities as calculated with DFT. A standardloss function L is the weighted sum of the mean squarederrors (MSEs) of the different observables:L = wE1NbNb∑µ=1[1Nµ(Epredµ − Erefµ)]2+ wF1NbNb∑µ=113NµNµ∑j=1∑α=x,y,z[F predµ,j,α − F refµ,j,α]2+ wσ1NbNb∑µ=119∑α=x,y,z∑β=x,y,z[σpredµ,α,β − σrefµ,α,β]2, (1)where wE , wF , and wσ are tunable weights. However,when employing the loss in Eq. (1) to train themodel [48], we noticed the appearance of molecular solid-like structures at high temperatures (T ∼ 1200K) duringthe dynamics. AIMD simulations at the same conditionsand system sizes show how these configurations arenot dynamically stable, and, therefore, are an artifactof the model. Even after the inclusion of thesestructures in the training set, the standard loss functionstill energetically favored these configurations (see theSupplemental Material (SM) for details [49]). Theatomic phase, on the other hand, seemed to be describedcorrectly.An analysis of the partial energy error distributionsfor this model, computed on subsets of the trainingset describing the selected phases, reveals how the lossfunction L produces a bias, i.e., a non-zero mean ofEpredµ − Erefµ . Indeed, being a function of the MSEs,L only guarantees that the total error distribution iscentered around zero.With this in mind, we modified the loss function byincluding a term that penalizes a global energy bias in theprediction error on the subset ofmolecular configurations,thus improving the description of this phase. Specifically,we used the modified loss L′ = L+∆L with∆L = wEλ1Nb∣∣∣∣∣∣∑µ∈mol1Nµ(Epredµ − Erefµ)∣∣∣∣∣∣, (2)where λ is a parameter controlling the relative weightof the penalty with respect to the energy term of thestandard loss function, wE is the same hyperparameteras in Eq. (1), and the sum only considers molecularconfigurations in a given training batch. To classifyconfigurations as either atomic or molecular, we used astatic criterion that solely depends on atomic positions:A configuration is classified as molecular if the first peakof its radial distribution function g(r) (estimated byfitting the g(r) with a Gaussian function for r ∈ [0, 1.3]Å)is larger than a threshold, here set to 1.8.We observed that the models trained using L′ havea much smaller energy error for the molecular solid-like structures, which do not appear anymore duringthe dynamics at high temperatures. As a consequenceof the modified loss, the energy error distribution onthe training set is slightly altered, with relative energiesbetween different types of structures in much betteragreement with respect to the ab initio PBE referencevalues, as discussed in the SM.To train the final MLIP, we compiled a datasetof 21 812 configurations. We extracted ∼ 17 000configurations of 128 hydrogen atoms each from theMD simulations of Refs. [50] and [51] (see the SM forthe density-temperature distribution of this set). Anadditional set of 3000 128-atoms configurations wasselected from a series of AIMD simulations at lowertemperatures (800 and 900K). We also added 500snapshots with 256 and 512 atoms extracted from MDruns with early iterations of the MACE model. Finally,we included solid structures (see Ref. [52]) and ∼ 100low-temperature 96-atoms configurations from Ref. [53].To ensure consistency across the dataset, we recomputedenergies, forces, and pressures at the PBE level, usingthe Quantum Espresso package [54–56]. A projectoraugmented-wave (PAW) pseudopotential [57] togetherwith a 60 Ry plane-waves cutoff was used. A sufficientlydense k-point grid for each system size was employed tocure finite-size effects. For instance, we used a 4× 4× 4grid for the configurations with 128 atoms.From this dataset, we randomly extracted 280configurations for testing. From the remainingconfigurations, 95% were used for training, and 5% forvalidation. This last group of configurations is usedduring the training to assess the performance of themodel. At the end of the optimization, the best model isselected as the one minimizing the loss on the validationset.Table I reports the accuracy of this MACE model,3FIG. 1. Results for MD simulations obtained with the MACE model for different temperatures: (a) T = 900 K, (b) T = 950 K,(c) T = 1000 K, (d) T = 1100 K. For each temperature, the left panel reports the EOS close to the LLT for the different systemsizes. For T = 900 K the red shaded area highlights the rs range for which hysteresis is observed at all values of N . The rightpanels show the average maxk S(k)/N along the trajectory (top right) and the average molecular fraction m (bottom right) asa function of the system size N for two different values of the Wigner-Seitz radius rs.RMSEE (meV/atom) RMSEF (meV/Å)Training 2.2 116Validation 2.1 118Test 2.0 116TABLE I. Root mean squared errors of the final MACE MLIPfor the energy per atom (RMSEE) and the forces (RMSEF )on the training, validation, and test sets.measured by the root mean squared error (RMSE) onthe energy per atom and the forces on the training,validation, and test sets. The test-set RMSE for the virialpressure is ∼ 1 GPa. Compared to the MLIP proposedin Ref. [34], our model has an error on energy 6−7 timessmaller and is twice as accurate on forces. The recentNequIP model of Ref. [44] shows a very similar RMSEon energy, but has a 40% larger RMSE on forces.LLT simulations: finite-size scaling. Thecomputational efficiency of the MACE model enabled usto study the behavior of the LLT as a function of bothtemperature and system size. We ran MD simulationsin the NVT ensemble using the LAMMPS code [58]interfaced with MACE. We considered systems withN = 128, 256, 512, and 2048 hydrogen atoms treatedas classical nuclei, and cubic supercells. We performedsimulations as long as 0.6 ns, about two orders ofmagnitude longer than what is usually achieved withAIMD. For the dynamics, we used a time step of0.2 fs and controlled the temperature via the stochasticvelocity rescaling scheme [59], with a characteristic timeof τ = 0.1 ps.To accurately study the nature of the transition,we calculated the EOS of the system in the vicinityof the LLT, using a dense grid of Wigner-Seitz radiirs, where 4π3(rsa0)3= VNwith V and a0 beingthe system volume and the Bohr radius, respectively.The results are presented in Fig. 1 for four differenttemperatures in the [900 − 1100] K range. From theEOS, we can identify three distinct regimes for the LLT:At the lowest temperature computed here, i. e., T =900 K (Fig. 1a), the model clearly predicts a first-ordertransition, signaled by the hysteresis present in the EOSfor all system sizes. At T = 950K and 1000K (Figs. 1band 1c, respectively), the EOS has a strong dependenceon the system size in the transition region: the smallersystems (i. e., N = 128 and 256) suggest again a first-order transition, while for larger N the pressure plateauand hysteresis are missing. Finally, at T ≥ 1100 K(Fig. 1d) the results indicate a smooth crossover, witha relatively small size dependence.To characterize the LLT, the sole observation ofthe EOS is not sufficient, since it is known that forthese transitions the density is not the primary orderparameter [60]. In particular, we computed the stablemolecular fraction [23, 30, 47] of the system, defined inthis case as the average number of hydrogen pairs whoseconstituent atoms stay within a distance of 1.05 Å for at4least a time τ ∼ 80 fs.The results corresponding to rs values slightlybefore/after the transition are also shown in Fig. 1. Inall cases, the molecular fraction rapidly increases fromvalues around 0.5 to values close to 1, signaling anatomic-to-molecular transition as rs increases.Besides the stable molecular fraction, we alsocomputed the structure factor to detect long-rangespatial correlations, i.e.,S(k) =〈1N∑i,jexp(ik · (Ri −Rj))〉,where ⟨· · · ⟩ indicates the average over the trajectory,R1, . . . ,RN are atom positions, and k = 2πL(n1, n2, n3),with L being the side of the cubic simulation box andn1, n2, n3 integer numbers. In particular, the maximumvalue of S(k) can be used as a proxy for the formationof crystalline structures, with maxk S(k) ∝ N in a solidsystem in the thermodynamic limit. On the contrary,in a liquid, the rotational invariance implies that thecontribution of each k to the S(k) will stay constantfor N → ∞. The value of the average maxk S(k)/Nas a function of the number of particles is reported inthe top right panels of Fig. 1. In the atomic phase,our simulations show that the system is a liquid at alltemperatures. In the molecular region, the behavior ofmaxk S(k)/N is instead strongly temperature dependent.Interestingly, the first-order transition observed in theEOS at T = 900 K is accompanied by a non-vanishingvalue of maxk S(k)/N in the thermodynamic limit,revealing a long-range spatial order of the molecularphase at this temperature (Fig. 1a). Moreover, at T =950 K and 1000 K, a large value of maxk S(k)/N ispresent for N = 128 and 256 on the molecular side, thenvanishing at larger N (Figs. 1(b-c)). This suggests thatthe first-order transition seen in smaller systems is due tofinite-size effects and it is between a molecular crystal andan atomic liquid. A genuine molecular liquid is insteadobserved for all system sizes at temperature T ≥ 1100 K(Fig. 1d).LLT simulations: finite-time effects. Our finite-sizescaling analysis highlights the importance of consideringsufficiently large systems to obtain converged results forthe LLT. This is further demonstrated in Fig. 2a, wherethe behavior of maxk S(k)/N is plotted as a function ofthe simulation time for two sizes, namely N = 256 and512, at T = 950 K and rs = 1.43. While the smallersystem is in a solid-like state for the majority of thetime, the system with N = 512 occasionally crystallizes,as indicated by the two jumps of maxk S(k)/N at thebeginning and the end of the MD run, but mostlyremains liquid. This behavior not only confirms the sizedependence already seen, but also shows that at least10 ps are necessary to melt the crystal for N = 512,suggesting that ∼ 100 ps are needed to obtain fullyconverged results near the transition.0 25 50 75 100 125 150 175 200Time (ps)0.10.20.3max k[S(k)]/N(a)T=950 K , rs=1.43N = 512N = 2561.48 1.49 1.50 1.51 1.52 1.53rs122126130134138142146Pressure (GPa)(b)T= 1400K, N = 128ab initioMACE 200 ps ±  (4 ps) ±2  (4 ps) ±3  (4 ps)1.425 1.430 1.435 1.440 1.445rs175177179181183185187Pressure (GPa)(c)T= 1000K, N = 512MACE 200 ps ±  (10 ps) ±2  (10 ps) ±3  (10 ps)FIG. 2. (a) Value of maxk S(k)/N as a function of thesimulation time for two systems, with N = 256 and N = 512atoms, respectively, at temperature T = 950 K and rs = 1.43.Colored lines represent running averages over a time windowof 8 ps. (b) Results for 200 ps long simulations, N = 128,T = 1400 K and corresponding confidence intervals, estimatedfrom a running average with time window τrun = 4 ps. TheAIMD result reported in Ref. [50] is also shown. (c) Same as(b) but with N = 512, T = 1000 K and τrun = 10 ps.To analyze the possibility that the kinks identified inthe EOS of previous AIMD results originate from theshort length of these simulations, we analyzed 0.2 ns longtrajectories obtained with the MACE model for values ofrs close to the transition at different temperatures aboveour Tc. The results are shown in Figs. 2b and 2c for thesystem ofN = 128 at T = 1400 K and the one ofN = 512at T = 1000 K, respectively. Using long trajectories, wecomputed the pressure running average using a variable-size time window τrun, corresponding to simulation timesachievable by AIMD, e. g., τrun ∼ 10 ps. The varianceof the running average gives an idea of the variabilityof the estimated equilibrium pressure when running MDsimulations of length equal to τrun. From Figs. 2b and2c, the estimated variance grows in the vicinity of the575 100 125 150 175 200Pressure (GPa)2468101214c v/KBT= 900KT= 950KT= 1000KT= 1100KT= 1200KT= 1400KT= 1600KT= 1800K256 512 1024 2048N102050max(cv)/k BFIG. 3. Specific heat per particle vs pressure along theisotherms, for N = 2048 obtained with the MACE model.The shaded areas indicate the uncertainty in the peakposition. In the inset, we report the size scaling of themaximum of cv/kB for four different temperatures (in log-log scale). The results confirm a first-order transition forT = 900 K, where cv scales with a behavior close tolinear, and a smooth crossover for higher temperatures in thethermodynamic limit.transition, and the presence of a pressure plateau in theEOS can be understood as an artifact due to the lack ofphase-space sampling.For the smallest system (N = 128, Fig. 2b), we alsocompared our results with the PBE AIMD ones reportedin Ref. [50] at T = 1400 K for the same size. Thepressure plateau here is within the 3σ uncertainty regionestimated from a running average of 4 ps, slightly longerthan the average length of the AIMD simulations. Theresults at lower temperature T = 1000 K for a largersystem of 512 atoms (Fig. 2c) show that even a 10 pslong dynamics can, in principle, produce artificial kinks ofthe size reported in the literature [21, 30]. Even though,in principle, this estimation depends on the algorithmicdetails, such as the choice of the thermostat, we do notexpect this to change our conclusions.Widom line. Both the first-order transition and thecrossover above the critical temperature Tc can be furthercharacterized by studying the behavior of the isochoricspecific heat (i. e., the heat capacity per particle), cv =1NkBT 2 (〈E2〉− ⟨E⟩2), as a function of temperature,density, and system size. The results for cv obtainedfrom simulations with N = 2048 hydrogen atoms attemperatures up to 1800 K are reported in Fig. 3. Foreach temperature, we estimated the location and valueof the maximum. For temperatures between 900 Kand 1100 K, we also studied the size scaling of the cvmaximum, as shown in the inset of Fig. 3. In a first-ordertransition, the cv peak presents a divergence scaling like∼ N in the thermodynamic limit [61]. As shown in thefigure, linear scaling is observed at T = 900 K, whilecv tends to a constant for higher temperatures. Thisis consistent with our EOS results locating the criticalpoint between T = 900 K and T = 950 K. The positionof the cv maximum at higher temperatures allows oneto determine the Widom line in the supercritical region.The results are reported in Fig. 4. Our Widom line showsremarkable agreement with previous simulations of theLLT obtained with the PBE functional. In particular,our results well reproduce the ones reported in Ref. [21]for temperatures up to 1400 K and those of Ref. [35] forhigher temperatures, even though they do not agree onthe LLT first-order character for T ≥ 900 K.100 125 150 175 200 225 250Pressure (GPa)800100012001400160018002000Temperature (K)Solid H2Liquid H2 Liquid HNiu et al. (Ref. 53), meltingLorenzen et al. (Ref. 21)Morales et al. (Ref. 22)Karasev et al. (Ref. 35)Cheng et al. (Ref. 34), max Cp Cheng et al. (Ref. 34), max  MACE PBE modelFIG. 4. Classical PBE-LLT location as computed withdifferent methods. AIMD results by Ref. [21] (light bluedashed line), Ref. [22] (green dashed line), Ref. [35] (bluedashed line). The results for the molecular-to-atomiccrossover obtained with an NN MLIP by Ref. [34] arereported with orange and violet markers, corresponding tothe maximum of the isobaric specific heat and density,respectively. The black markers and line indicate the recentlyproposed PBE melting line by Ref. [53], obtained with an NNMLIP and the two-phase method. The MACE model resultsare indicated with red markers. The filled point at T = 900 Kindicates the first-order character of the transition, while theempty points correspond to the location of the Widom linegiven by the cv maximum (see Fig. 3).Conclusions. Enabled by an accurate and efficientMLIP based on the MACE architecture, we show that theLLT between molecular and atomic hydrogen is alwayssupercritical above the melting temperature. For this, weperformed NVT MD simulations in a temperature rangebetween 900K and 1800K of systems with up to 2048hydrogen atoms and simulation times of up to 600 ps.We preserved the ab initio PBE quality of the MACEmodel across the transition by virtue of a modified lossfunction, yielding consistent energy differences betweenatomic and molecular configurations.For temperatures above 900K, the MACE modelpredicts a proper LLT in the thermodynamic limit (N ≥512), which turns out to be a smooth crossover. Up toT = 1000K, defective crystal structures are observed forsmall systems due to strong finite-size effects, giving riseto a first-order transition that disappears as N increases.At the lowest temperature (900K), the model predictsa first-order transition in the thermodynamic limit. The6structure factor analysis reveals the presence of long-range spatial order in the molecular phase, showing thatthis is a transition between a molecular solid-like system,frustrated by the cubic cell, and an atomic liquid. Thistransition point lies exactly on the PBE melting lineof Ref. [53] (see Fig. 4), calculated with an NN MLIPobtained using the DeePMD package [62], suggestingthat the LLT might become first-order inside the crystalregion for T ≤ 900 K.The physical picture provided by the MACE modelqualitatively agrees with Ref. [34] on the supercriticalnature of the LLT, even though their model couldnot accurately reproduce the crystallization regime, andpredicted a Widom line (estimated from the maximum ofthe isobaric specific heat cp) far from the AIMD results(see Fig. 4). The first-order character of the transitionobserved in other studies may be explained by the largefluctuations expected close to the Widom line, whichcould lead to the incorrect identification of density jumpsand pressure plateaux in the EOS for too short simulationtimes. Additional path integral MD simulations carriedout with the MACE model indicate that the conclusionsabout the crossover nature of the LLT are not changedby the inclusion of nuclear quantum effects, as reportedin the SM [49].These results urgently call for a reinvestigation ofthe LLT using ab initio methods beyond DFT, such asQMC. The ∆-learning scheme [50, 51] combined withthe present MACE model taken as baseline could allowone to access unprecedented system sizes and simulationlengths, by extending this work to QMC-based MLIPapplications.Data availability. A version of the MACE codeimplementing the modified loss function employed in thiswork can be found at Ref. [63]. The model, training set,and simulation results are available at Ref. [64].Acknowledgments. K.N. acknowledges financialsupport from MEXT Leading Initiative for ExcellentYoung Researchers (Grant No. JPMXS0320220025)and from by JST BOOST (Grant No. JPMJBY24F3).B.J. acknowledges support by the state of Baden-Württemberg through bwHPC. M. C. thanks theEuropean High Performance Computing JointUndertaking (JU) for the partial support throughthe ”EU-Japan Alliance in HPC” HANAMI project(Hpc AlliaNce for Applications and supercoMputingInnovation: the Europe - Japan collaboration).This work has received funding from the EuropeanCenter of Excellence in Exascale Computing TREX(Targeting Real chemical accuracy at the Exascale,grant 952165). The authors acknowledge discussionswith Prof. David Ceperley, Prof. Michele Ceriotti,Prof. Gábor Csányi, Prof. Carlo Pierleoni, Dr. ThomasBischoff, Dr. Marco Cherubini, and Dr. AbhishekRaghav. M.C. acknowledges GENCI for providingcomputational resources on the CEA-TGCC Irene andIDRIS Jean-Zay supercomputer clusters under projectnumber A0170906493 and the TGCC special session.[1] M. Bonitz, J. Vorberger, M. Bethkenhagen, M. P. Böhme,D. M. Ceperley, A. Filinov, T. Gawne, F. Graziani,G. Gregori, P. Hamann, S. B. Hansen, M. Holzmann,S. X. Hu, H. Kählert, V. V. Karasiev, U. Kleinschmidt,L. Kordts, C. Makait, B. Militzer, Z. A. Moldabekov,C. Pierleoni, M. Preising, K. Ramakrishna, R. Redmer,S. Schwalbe, P. Svensson, and T. Dornheim, Toward firstprinciples-based simulations of dense hydrogen, Physicsof Plasmas 31, 110501 (2024).[2] R. Helled, G. Mazzola, and R. Redmer, Understandingdense hydrogen at planetary conditions, Nature ReviewsPhysics 2, 562 (2020).[3] J. E. P. Connerney, M. Benn, J. B. Bjarno, T. Denver,J. Espley, J. L. Jorgensen, P. S. Jorgensen, P. Lawton,A. Malinnikova, J. M. Merayo, S. Murphy, J. Odom,R. Oliversen, R. Schnurr, D. Sheppard, and E. J. Smith,The Juno Magnetic Field Investigation, Space ScienceReviews 213, 39 (2017).[4] S. T. Weir, A. C. Mitchell, and W. J. Nellis, Metallizationof fluid molecular hydrogen at 140 GPa (1.4 Mbar),Physical Review Letters 76, 1860 (1996).[5] V. E. Fortov, R. I. Ilkaev, V. A. Arinin, V. V.Burtzev, V. A. Golubev, I. L. Iosilevskiy, V. V.Khrustalev, A. L. Mikhailov, M. A. Mochalov, V. Y.Ternovoi, and M. V. Zhernokletov, Phase transition in astrongly nonideal deuterium plasma generated by quasi-isentropical compression at megabar pressures, PhysicalReview Letters 99, 185001 (2007).[6] V. Dzyabura, M. Zaghoo, and I. F. Silvera, Evidence ofa liquid–liquid phase transition in hot dense hydrogen,Proceedings of the National Academy of Sciences 110,8040 (2013).[7] K. Ohta, K. Ichimaru, M. Einaga, S. Kawaguchi,K. Shimizu, T. Matsuoka, N. Hirao, and Y. Ohishi, Phaseboundary of hot dense fluid hydrogen, Scientific Reports5, 16560 (2015).[8] M. D. Knudson, M. P. Desjarlais, A. Becker, R. W.Lemke, K. R. Cochrane, M. E. Savage, D. E. Bliss,T. R. Mattsson, and R. Redmer, Direct observation ofan abrupt insulator-to-metal transition in dense liquiddeuterium, Science 348, 1455 (2015).[9] R. S. McWilliams, D. A. Dalton, M. F. Mahmood, andA. F. Goncharov, Optical properties of fluid hydrogenat the transition to a conducting state, Physical ReviewLetters 116, 255501 (2016).[10] M. Zaghoo, A. Salamat, and I. F. Silvera, Evidence of afirst-order phase transition to metallic hydrogen, PhysicalReview B 93, 155128 (2016).[11] M. Zaghoo and I. F. Silvera, Conductivity anddissociation in liquid metallic hydrogen and implicationsfor planetary interiors, Proceedings of the NationalAcademy of Sciences 114, 11873 (2017).[12] M. Zaghoo, R. J. Husband, and I. F. Silvera, Strikingisotope effect on the metallization phase lines of liquidhydrogen and deuterium, Physical Review B 98, 104102(2018).https://doi.org/10.1063/5.0219405https://doi.org/10.1063/5.0219405https://doi.org/10.1038/s42254-020-0223-3https://doi.org/10.1038/s42254-020-0223-3https://doi.org/10.1007/s11214-017-0334-zhttps://doi.org/10.1007/s11214-017-0334-zhttps://doi.org/10.1103/PhysRevLett.76.1860https://doi.org/10.1103/PhysRevLett.99.185001https://doi.org/10.1103/PhysRevLett.99.185001https://doi.org/10.1073/pnas.1300718110https://doi.org/10.1073/pnas.1300718110https://doi.org/10.1038/srep16560https://doi.org/10.1038/srep16560https://doi.org/10.1126/science.aaa7471https://doi.org/10.1103/PhysRevLett.116.255501https://doi.org/10.1103/PhysRevLett.116.255501https://doi.org/10.1103/PhysRevB.93.155128https://doi.org/10.1103/PhysRevB.93.155128https://doi.org/10.1073/pnas.1707918114https://doi.org/10.1073/pnas.1707918114https://doi.org/10.1103/PhysRevB.98.104102https://doi.org/10.1103/PhysRevB.98.1041027[13] 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. L. Pape, G. W. Collins, R. Jeanloz, andR. J. Hemley, Insulator-metal transition in dense fluiddeuterium, Science 361, 677 (2018).[14] S. Jiang, N. Holtgrewe, Z. M. Geballe, S. S. Lobanov,M. F. Mahmood, R. S. McWilliams, and A. F.Goncharov, A spectroscopic study of the insulator–metaltransition in liquid hydrogen and deuterium, AdvancedScience 7, 1901668 (2020).[15] S. Scandolo, Liquid–liquid phase transition in compressedhydrogen from first-principles simulations, Proceedingsof the National Academy of Sciences 100, 3051 (2003).[16] S. A. Bonev, B. Militzer, and G. Galli, Ab initiosimulations of dense liquid deuterium: Comparison withgas-gun shock-wave experiments, Physical Review B 69,014101 (2004).[17] K. T. Delaney, C. Pierleoni, and D. M. Ceperley,Quantum Monte Carlo simulation of the high-pressuremolecular-atomic crossover in fluid hydrogen, PhysicalReview Letters 97, 235702 (2006).[18] J. Vorberger, I. Tamblyn, B. Militzer, and S. A. Bonev,Hydrogen-helium mixtures in the interiors of giantplanets, Physical Review B 75, 024206 (2007).[19] B. Holst, R. Redmer, and M. P. Desjarlais,Thermophysical properties of warm dense hydrogenusing quantum molecular dynamics simulations,Physical Review B 77, 184201 (2008).[20] C. Attaccalite and S. Sorella, Stable liquid hydrogen athigh pressure by a novel ab initio molecular-dynamicscalculation, Physical Review Letters 100, 114501 (2008).[21] W. Lorenzen, B. Holst, and R. Redmer, First-orderliquid-liquid phase transition in dense hydrogen, PhysicalReview B 82, 195107 (2010).[22] M. A. Morales, C. Pierleoni, E. Schwegler, andD. M. Ceperley, Evidence for a first-order liquid-liquidtransition in high-pressure hydrogen from ab initiosimulations, Proceedings of the National Academy ofSciences 107, 12799 (2010).[23] I. Tamblyn and S. A. Bonev, Structure and phaseboundaries of compressed liquid hydrogen, PhysicalReview Letters 104, 065702 (2010).[24] G. Mazzola, S. Yunoki, and S. Sorella, Unexpectedly highpressure for molecular dissociation in liquid hydrogen byelectronic simulation, Nature Communications 5, 3487(2014).[25] J. Yang, C. L. Tian, F. S. Liu, L. C. Cai, H. K. Yuan,M. M. Zhong, and F. Xiao, A new evidence of first-orderphase transition for hydrogen at 3000 K, EurophysicsLetters 109, 36003 (2015).[26] C. Pierleoni, M. A. Morales, G. Rillo, M. Holzmann,and D. M. Ceperley, Liquid–liquid phase transitionin hydrogen by coupled electron–ion Monte Carlosimulations, Proceedings of the National Academy ofSciences 113, 4953 (2016).[27] G. E. Norman and I. M. Saitov, Critical point andmechanism of the fluid–fluid phase transition in warmdense hydrogen, Doklady Physics 62, 294 (2017).[28] G. Mazzola, R. Helled, and S. Sorella, Phase diagram ofhydrogen and a hydrogen-helium mixture at planetaryconditions by quantum Monte Carlo simulations,Physical Review Letters 120, 025701 (2018).[29] C. Tian, F. Liu, H. Yuan, H. Chen, and A. Kuan,First-order liquid-liquid phase transition in compressedhydrogen and critical point, The Journal of ChemicalPhysics 150, 204114 (2019).[30] H. Y. Geng, Q. Wu, M. Marqués, and G. J. Ackland,Thermodynamic anomalies and three distinct liquid-liquid transitions in warm dense liquid hydrogen,Physical Review B 100, 134109 (2019).[31] G. Rillo, M. A. Morales, D. M. Ceperley, and C. Pierleoni,Optical properties of high-pressure fluid hydrogen acrossmolecular dissociation, Proceedings of the NationalAcademy of Sciences 116, 9770 (2019).[32] J. Hinz, V. V. Karasiev, S. X. Hu, M. Zaghoo,D. Mej́ıa-Rodŕıguez, S. B. Trickey, and L. Caldeŕın,Fully consistent density functional theory determinationof the insulator-metal transition boundary in warm densehydrogen, Phys. Rev. Res. 2, 032065 (2020).[33] T. Bryk, C. Pierleoni, G. Ruocco, and A. P. Seitsonen,Characterization of molecular-atomic transformationin fluid hydrogen under pressure via long-wavelengthasymptote of charge density fluctuations, Journal ofMolecular Liquids 312, 113274 (2020).[34] B. Cheng, G. Mazzola, C. J. Pickard, and M. Ceriotti,Evidence for supercritical behaviour of high-pressureliquid hydrogen, Nature 585, 217 (2020).[35] V. V. Karasiev, J. Hinz, S. Hu, and S. Trickey, On theliquid–liquid phase transition of dense hydrogen, Nature600, E12 (2021).[36] A. Bergermann, L. Kleindienst, and R. Redmer,Nonmetal-to-metal transition in liquid hydrogenusing density functional theory and theHeyd–Scuseria–Ernzerhof exchange-correlationfunctional, The Journal of Chemical Physics 161,234303 (2024).[37] G. Gliaudelis, V. Lukyanchuk, N. Chtchelkatchev,I. Saitov, and N. Kondratyuk, Dynamical propertiesof hydrogen fluid at high pressures, The Journal ofChemical Physics 162, 024504 (2025).[38] A. F. Goncharov, I. I. Mazin, J. H. Eggert, R. J. Hemley,and H.-k. Mao, Invariant points and phase transitions indeuterium at megabar pressures, Physical Review Letters75, 2514 (1995).[39] W. A. Bassett, Diamond anvil cell, 50th birthday, HighPressure Research 29, 163 (2009).[40] W. J. Nellis, Dynamic compression of materials:metallization of fluid hydrogen at high pressures, Reportson Progress in Physics 69, 1479 (2006).[41] M. A. Morales, R. Clay, C. Pierleoni, and D. M. Ceperley,First principles methods: A perspective from quantumMonte Carlo, Entropy 16, 287 (2014).[42] J. Behler, Perspective: Machine learning potentials foratomistic simulations, The Journal of Chemical Physics145, 170901 (2016).[43] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P.Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, andB. Kozinsky, E(3)-equivariant graph neural networks fordata-efficient and accurate interatomic potentials, NatureCommunications 13, 2453 (2022).[44] M. Istas, S. Jensen, Y. Yang, M. Holzmann, C. Pierleoni,and D. M. Ceperley, Liquid-liquid phase transition ofhydrogen and its critical point: Analysis from ab initiosimulation and a machine-learned potential, Phys. Rev.E 111, 045307 (2025).https://doi.org/10.1126/science.aat0970https://doi.org/https://doi.org/10.1002/advs.201901668https://doi.org/https://doi.org/10.1002/advs.201901668https://doi.org/10.1073/pnas.0038012100https://doi.org/10.1073/pnas.0038012100https://doi.org/10.1103/PhysRevB.69.014101https://doi.org/10.1103/PhysRevB.69.014101https://doi.org/10.1103/PhysRevLett.97.235702https://doi.org/10.1103/PhysRevLett.97.235702https://doi.org/10.1103/PhysRevB.75.024206https://doi.org/10.1103/PhysRevB.77.184201https://doi.org/10.1103/PhysRevLett.100.114501https://doi.org/10.1103/PhysRevB.82.195107https://doi.org/10.1103/PhysRevB.82.195107https://doi.org/10.1073/pnas.1007309107https://doi.org/10.1073/pnas.1007309107https://doi.org/10.1103/PhysRevLett.104.065702https://doi.org/10.1103/PhysRevLett.104.065702https://doi.org/10.1038/ncomms4487https://doi.org/10.1038/ncomms4487https://doi.org/10.1209/0295-5075/109/36003https://doi.org/10.1209/0295-5075/109/36003https://doi.org/10.1073/pnas.1603853113https://doi.org/10.1073/pnas.1603853113https://doi.org/10.1134/S1028335817060088https://doi.org/10.1103/PhysRevLett.120.025701https://doi.org/10.1063/1.5096400https://doi.org/10.1063/1.5096400https://doi.org/10.1103/PhysRevB.100.134109https://doi.org/10.1073/pnas.1818897116https://doi.org/10.1073/pnas.1818897116https://doi.org/10.1103/PhysRevResearch.2.032065https://doi.org/https://doi.org/10.1016/j.molliq.2020.113274https://doi.org/https://doi.org/10.1016/j.molliq.2020.113274https://doi.org/https://doi.org/10.1038/s41586-020-2677-yhttps://doi.org/https://doi.org/10.1038/s41586-021-04078-xhttps://doi.org/https://doi.org/10.1038/s41586-021-04078-xhttps://doi.org/10.1063/5.0241111https://doi.org/10.1063/5.0241111https://doi.org/10.1063/5.0236394https://doi.org/10.1063/5.0236394https://doi.org/10.1103/PhysRevLett.75.2514https://doi.org/10.1103/PhysRevLett.75.2514https://doi.org/10.1080/08957950802597239https://doi.org/10.1080/08957950802597239https://doi.org/10.1088/0034-4885/69/5/r05https://doi.org/10.1088/0034-4885/69/5/r05https://doi.org/10.3390/e16010287https://doi.org/10.1063/1.4966192https://doi.org/10.1063/1.4966192https://doi.org/10.1038/s41467-022-29939-5https://doi.org/10.1038/s41467-022-29939-5https://doi.org/10.1103/PhysRevE.111.045307https://doi.org/10.1103/PhysRevE.111.0453078[45] B. Cheng, G. Mazzola, C. J. Pickard, and M. Ceriotti,Reply to: On the liquid–liquid phase transition of densehydrogen, Nature 600, E15 (2021).[46] I. Batatia, D. P. Kovacs, G. Simm, C. Ortner, andG. Csanyi, Mace: Higher order equivariant messagepassing neural networks for fast and accurate force fields,in Advances in Neural Information Processing Systems,Vol. 35, edited by S. Koyejo, S. Mohamed, A. Agarwal,D. Belgrave, K. Cho, and A. Oh (Curran Associates, Inc.,2022) pp. 11423–11436.[47] T. Bischoff, B. Jäckl, and M. Rupp, Hydrogenunder Pressure as a Benchmark for Machine-LearningInteratomic Potentials, arXiv , 2409.13390 (2024).[48] For building the model, we used 128 equivariantmessages, a correlation order of 3, and a cutoff radiusof rc = 3 Å.[49] See Supplemental Material at [URL will be insertedby publisher ] for further details about the training ofthe MACE model, a comparison between the errordistribution of two models trained using the standardand modified loss, an analysis of the defective solid foundin our simulations, a discussion of the path integralmolecular dynamics results, and a comparison betweenour results and those obtained from ab initio moleculardynamics. (see references [26, 35, 47, 50, 51, 53, 65–67]therein).[50] A. Tirelli, G. Tenti, K. Nakano, and S. Sorella, High-pressure hydrogen by machine learning and quantumMonte Carlo, Physical Review B 106, L041105 (2022).[51] G. Tenti, K. Nakano, A. Tirelli, S. Sorella, and M. Casula,Principal deuterium Hugoniot via quantum Monte Carloand ∆-learning, Physical Review B 110, L041107 (2024).[52] L. Monacelli, M. Casula, K. Nakano, S. Sorella, andF. Mauri, Quantum phase diagram of high-pressurehydrogen, Nature Physics 19, 845 (2023).[53] H. Niu, Y. Yang, S. Jensen, M. Holzmann, C. Pierleoni,and D. M. Ceperley, Stable solid molecular hydrogenabove 900 K from a machine-learned potential trainedwith diffusion quantum Monte Carlo, Physical ReviewLetters 130, 076102 (2023).[54] 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: amodular and open-source software project for quantumsimulations of materials, Journal of Physics: CondensedMatter 21, 395502 (2009).[55] 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. Küçükbenli,M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L.Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto,S. Poncé, 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, Advancedcapabilities for materials modelling with QUANTUMESPRESSO, Journal of Physics: Condensed Matter 29,465901 (2017).[56] P. Giannozzi, O. Baseggio, P. Bonfà, 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 towardthe exascale, The Journal of Chemical Physics 152,154105 (2020).[57] H.pbek-jpaw psl.1.0.0.UPF pseudopotentials availableat http://pseudopotentials.quantum-espresso.org/legacy_tables/ps-library/h.[58] A. P. Thompson, H. M. Aktulga, R. Berger, D. S.Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’tVeld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan,M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton,LAMMPS—a flexible simulation tool for particle-basedmaterials modeling at the atomic, meso, and continuumscales, Computer Physics Communications 271, 108171(2022).[59] G. Bussi, D. Donadio, and M. Parrinello, Canonicalsampling through velocity rescaling, The Journal ofChemical Physics 126, 014101 (2007).[60] H. Tanaka, Liquid–liquid transition and polyamorphism,The Journal of Chemical Physics 153, 130901 (2020).[61] K. Binder, Theory of first-order phase transitions,Reports on Progress in Physics 50, 783 (1987).[62] J. Zeng, D. Zhang, D. Lu, P. Mo, Z. Li, Y. Chen,M. Rynik, L. Huang, Z. Li, S. Shi, Y. Wang, H. Ye,P. Tuo, J. Yang, Y. Ding, Y. Li, D. Tisi, Q. Zeng, H. Bao,Y. Xia, J. Huang, K. Muraoka, Y. Wang, J. Chang,F. Yuan, S. L. Bore, C. Cai, Y. Lin, B. Wang, J. Xu, J.-X. Zhu, C. Luo, Y. Zhang, R. E. A. Goodall, W. Liang,A. K. Singh, S. Yao, J. Zhang, R. Wentzcovitch, J. Han,J. Liu, W. Jia, D. M. York, W. E, R. Car, L. Zhang, andH. Wang, DeePMD-kit v2: A software package for deeppotential models, The Journal of Chemical Physics 159,054801 (2023).[63] G. Tenti, https://github.com/giacomotenti/mace(2025).[64] G. Tenti, B. Jäckl, K. Nakano, M. Rupp, and M. Casula,Additional data for ”hydrogen liquid-liquid transitionfrom first principles and machine learning” (2025).[65] D. P. Kingma and J. Ba, Adam: A method forstochastic optimization, in 3rd International Conferenceon Learning Representations (ICLR), San Diego,California, USA, May 7–9 , edited by Y. Bengio andY. LeCun (2015).[66] K. Momma and F. Izumi, VESTA: a three-dimensionalvisualization system for electronic and structuralanalysis, Journal of Applied Crystallography 41, 653(2008).[67] F. Mouhat, S. Sorella, R. Vuilleumier, A. M. Saitta, andM. Casula, Fully quantum description of the Zundel ion:Combining variational quantum Monte Carlo with pathintegral Langevin dynamics, J. Chem. Theory Comput.13, 2400 (2017).https://doi.org/https://doi.org/10.1038/s41586-021-04079-whttps://proceedings.neurips.cc/paper_files/paper/2022/file/4a36c3c51af11ed9f34615b81edb5bbc-Paper-Conference.pdfhttps://doi.org/10.48550/arXiv.2409.13390https://doi.org/10.1103/PhysRevB.106.L041105https://doi.org/10.1103/PhysRevB.110.L041107https://doi.org/10.1038/s41567-023-01960-5https://doi.org/10.1103/PhysRevLett.130.076102https://doi.org/10.1103/PhysRevLett.130.076102https://doi.org/10.1088/0953-8984/21/39/395502https://doi.org/10.1088/0953-8984/21/39/395502https://doi.org/10.1088/1361-648x/aa8f79https://doi.org/10.1088/1361-648x/aa8f79https://doi.org/10.1063/5.0005082https://doi.org/10.1063/5.0005082http://pseudopotentials.quantum-espresso.org/legacy_tables/ps-library/hhttp://pseudopotentials.quantum-espresso.org/legacy_tables/ps-library/hhttps://doi.org/https://doi.org/10.1016/j.cpc.2021.108171https://doi.org/https://doi.org/10.1016/j.cpc.2021.108171https://doi.org/10.1063/1.2408420https://doi.org/10.1063/1.2408420https://doi.org/10.1063/5.0021045https://doi.org/10.1088/0034-4885/50/7/001https://doi.org/10.1063/5.0155600https://doi.org/10.1063/5.0155600https://github.com/giacomotenti/macehttps://doi.org/10.5281/zenodo.14790036https://doi.org/10.5281/zenodo.14790036http://arxiv.org/abs/1412.6980http://arxiv.org/abs/1412.6980http://arxiv.org/abs/1412.6980https://doi.org/10.1107/S0021889808012016https://doi.org/10.1107/S0021889808012016https://doi.org/10.1021/acs.jctc.7b00017https://doi.org/10.1021/acs.jctc.7b00017Supplemental Material for “Hydrogen liquid-liquid transition from first principles and machinelearning”Giacomo Tenti,1, ∗ Bastian Jäckl,2 Kousuke Nakano,3, 4 Matthias Rupp,2, 5 and Michele Casula6, †1International School for Advanced Studies (SISSA), Via Bonomea 265, 34136 Trieste, Italy2Department of Computer and Information Science, University of Konstanz, Box 188, 78457 Konstanz, Germany3Center for Basic Research on Materials, National Institute for Materials Science (NIMS), 1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, Japan4Center for Emergent Matter Science (CEMS), RIKEN, 2-1 Hirosawa, Wako, Saitama, 351-0198, Japan.5Scientific Instrumentation and Process Technology unit,Luxembourg Institute of Science and Technology (LIST), Maison de l’Innovation,5 avenue des Hauts-Fourneaux, L-4362 Esch-sur-Alzette, Luxembourg6Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC),Sorbonne Université, CNRS UMR 7590, MNHN, 4 Place Jussieu, 75252 Paris, France(Dated: November 18, 2025)I. MACE MODEL TRAINING DETAILSA. PBE dataset density-temperature distributionAs mentioned in the main text, the dataset used for both the model training and validation is a combination of configurationsextracted from multiple sources [1–3] and from other trajectories obtained with both ab initio molecular dynamics (AIMD) anddynamics using early iterations of the model. A majority of configurations (∼ 17000) were sampled from variation Monte Carlo-based trajectories of Ref. [1]. This set spans an rs range from 1.26 to 1.62, and temperatures from T = 1000 K to T = 1500 K.The dataset comprises both molecular and atomic configurations, as we can observe by computing the stable molecular fraction(as defined in the main text), with an H2 lifetime threshold of 20 ps. The distribution of the set in the rs − T space is reported inFig. 1.0 1FIG. 1: Distribution of configurations sampled from MD runs of Ref. [1] to generate part of the dataset of the MACE model.The dimension of the dots is proportional to the number of configurations extracted from each MD simulation, while the colorindicates the stable molecular fraction.∗ gtenti@sissa.it† michele.casula@impmc.upmc.fr2B. Traning parametersFor the final MACE model, we used 128 equivariant messages, a correlation order of 3, two message-passing layers, and acutoff radius of rc = 3 Å. The training was performed using the Adam optimizer [4] with a Nb = 16 batch size and initial learningrate of 0.01, and using the modified loss function L′ defined in the main text. The loss weights were taken equal to wE = 1,wF = 100, wσ = 100 for the first 320 epochs (defined as Nt/Nb optimization steps, with Nt being the total training set dimension)and wE = 200, wF = 10, wσ = 10 for the remaining 130 epochs. This is done to reduce the error on the energy after the forcesare converged. The relative weight of the energy penalty ∆L (see Eq. 2 of the main text) was set to λ = 50.C. Modified loss effect on the error distributionThe effect of using the modified loss function was studied by comparing the energy error distribution 1Nµ(Epredµ − Erefµ ) acrossthe training set for two different models. In particular, we considered the final MACE model employed for the simulations (usingthe modified loss), and a model trained with a standard loss function (Eq. 1 of the main text) and the same training parameters.Moreover, we split the training set into three partitions:1. Solid-like configurations extracted from molecular dynamics (MD) simulations with previous iterations of the model(trained using a standard loss);2. molecular configurations not belonging to the first partition (both solid and liquid);3. liquid atomic configurations.The classification between molecular and atomic configurations followed the static criterion used in the definition of ∆L,based on the analysis of the radial distribution function (RDF) first peak. The results for the error distributions are shown inFigs. 2 and 3 for the standard and modified loss, respectively. In Tab. I we also report, for each model and training set partition,the values for the root mean squared error (RMSE), the mean absolute error (MAE) and the average energy shift:∆E =1|S |∑µ∈S1Nµ(Epredµ − Erefµ),where S is one of the training set partitions and |S | is the number of configurations belonging to it.Atomic liquid Solid-like Remaining molecularRMSEE MAEE ∆E RMSEE MAEE ∆E RMSEE MAEE ∆EStandard 2.2 1.7 −0.3 7.2 6.8 −6.6 2.6 1.8 0.9Modified 2.3 1.8 −0.9 2.6 2.3 −2.2 1.9 1.4 −0.2TABLE I: Values of the energy RMSE, MAE and shift ∆E for the different training set partitions as obtained with two modelstrained with the standard and modified loss, respectively. All values are reported in meV/atom.We can immediately notice that the standard model dramatically underestimates the energy of the first partition, thus ex-plaining the appearance of these configurations during the dynamics even at relatively high temperatures (T ∼ 1200 K). Werecognized these configurations as unphysical since they do not appear in AIMD, contrary to the solid-like configurations de-scribed in the main text, which are at the origin of the first-order transition at T = 900 K. The use of the modified loss stronglyimproves the energy of this partition, decreasing the bias towards negative energies. The analysis of the error distribution of theother two partitions is also interesting. The RMSE, MAE and absolute value of the shift ∆E relative to the ”remaining molecular”partition are all reduced in the modified model with respect to the standard one. This is not surprising, since the extra term inthe loss function ∆L is directly proportional to ∆E evaluated on the molecular configurations. On the contrary, a larger absolutevalue of the bias is observed in the modified model for the ”atomic liquid” partition, even though the two models have almostthe same RMSE and MAE. Moreover, notice how the relative difference between the energy shift on the molecular and atomicpartitions is overall smaller in the model using the modified loss. We can expect this quantity to be crucial to accurately describethe PBE liquid-liquid phase transition (LLT), so that the modified model should necessarily be better for this task.Finally, we mention that we tested different forms for the energy penalty ∆L, also including an atomic term, e.g.,∆L = wEλ1Nb∣∣∣∣∣∣∣∣∑µ∈mol1Nµ(Epredµ − Erefµ)∣∣∣∣∣∣∣∣+∣∣∣∣∣∣∣∑µ∈atom1Nµ(Epredµ − Erefµ)∣∣∣∣∣∣∣.3The resulting model showed a similar accuracy, but was overall harder to optimize. In the end, we chose to use a correctiondepending only on the molecular energy shift.10 5 0 5 10Epred Eref (meV / atom) 0100200300400500600# countsStandard lossAtomic liquid configurationsRemaining molecular configurationsFrom solid-like trajectoriesFIG. 2: Histogram of the difference between the energy per atom predicted by a model trained using a ”standard” loss and thereference PBE value for different partitions of the training set.10 5 0 5 10Epred Eref (meV / atom) 0100200300400500600# countsModified lossAtomic liquid configurationsRemaining molecular configurationsFrom solid-like trajectoriesFIG. 3: Histogram of the difference between the energy per atom predicted by a model trained using a ”modified” loss(introduced in the main text) and the reference PBE value for different partitions of the training set.4II. SOLID-LIKE STRUCTURES ANALYSISTo further characterize the solid-like configurations at the lowest temperature T = 900 K, we computed the mean squareddisplacement (MSD) along the trajectory:MSD(t) =1NN∑i=1|Ri(t) − Ri(0)|2 . (1)The MSD computed for T = 900 K and T = 1000 K for different values of rs is shown in Fig. 4. Notice how the ”solid-like”system (corresponding to T = 900 and rs = 1.43) shows a non-zero diffusivity. In Ref. [5], this observation was used as anindication of the liquid phase. Here we obtained very similar values of the MSD (see the supplementary material of Ref. [5]),although our measure of the maxk S (k)/N clearly indicates the presence of long-range spatial correlations.FIG. 4: MSD as a function of time for (a) T = 900 K and rs = 1.40, rs = 1.43 (b) T = 1000 K and rs = 1.425, rs = 1.4425.The structure of these configurations can be appreciated in Fig. 5(a), where we show a snapshot taken from an MD simulationat T = 900 K and rs = 1.43 for a 2048 atoms system. The presence of ”planes” of hydrogen molecules in the system can benoticed. The analysis of the S (k) reveals multiple groups of non-collinear reciprocal vectors ±k with similar contributions (andscaling with N), confirming the solid nature of this phase. In Fig. 5(b) and 5(c), we also report structures at higher temperatureT = 1000 K in both the molecular and atomic phase. The radial distribution function and spherically averaged structure factorS (q) are also shown to highlight their difference. The atomic snapshots images are obtained using VESTA [6].To investigate the underlying crystal symmetry, we performed DFT structural optimization of the system using the QuantumEspresso package. In particular, we relaxed the atomic coordinates and cell under a constant pressure of p = 200 GPa. Analysisof the final structure reveals a Pbca-16 symmetry. We compared the final ground state DFT-PBE energy per atom of this systemwith the ones obtained considering a Fmmm-4 structure (proposed in Ref. [3] as the correct solid symmetry near the meltingline) and a C2c-24 structure, which should well represent the correct crystal at T = 0 K for this pressure. Also for these othersymmetries the structure was relaxed at p = 200 GPa. We computed the energies for each structure with both the modifiedMACE model and a model trained on a standard loss function (see Sec. I C). The modified model predicts a higher energyfor both the Fmmm-4 and Pbca-16 structures, while correctly reproducing their near degeneracy. In principle, this can lead tomelting via structural frustration. On the contrary, the model trained with the standard loss function has a smaller error for theFmmm-4 symmetry, but does not capture the energy degeneracy of the two structures. Finally, we note a higher accuracy of themodified model for the C2c-24 crystal.Energy per atom (eV/atom)PBE Standard Loss Modified LossC2c-24 −14.625 −14.611 −14.624Fmmm-4 −14.584 −14.582 −14.568Pbca-16 −14.585 −14.568 −14.569TABLE II: Energy per atom for different solid phase symmetries as computed with ab initio PBE, and the models trained with astandard loss function and a modified one, respectively.5(a)(b)(c)T = 900K,  = 1.43 rsT = 1000K,  = 1.45 rsT = 1000K,  = 1.425 rsFIG. 5: 2048 atoms structures extracted from classical MD simulations at different temperatures and rs values: (a) T = 900 Kand rs = 1.43 (molecular crystal); (b) T = 1000 K and rs = 1.45 (molecular liquid); (c) T = 1000 K and rs = 1.425 (atomicliquid). The dashed lines in panel (a) are a guide for the eye to highlight the intermolecular planes. For each structure, wereport the corresponding radial distribution function (RDF) (gray line), together with the approximant (red line) employed forthe classification of the molecular environments in the modified loss function (the threshold of 1.8 on the RDF peak is alsoshown). The spherically averaged structure factor S (q) is also depicted. The left images were obtained using the VESTAvisualization program [6]. Bonds are shown between atoms at a distance < 0.8 Å.61.44 1.45 1.46 1.47 1.48 1.49 1.50rs140150160Pressure (GPa)(a)0 2 4r (Å)0.00.51.01.52.0RDFm=0.492±0.011m=0.925±0.011(b)0 2 4Time (ps)0.00.10.20.3max k[S(k)]/N(c) rs=1.47rs=1.475T = 900 KFIG. 6: Results of PIOUD simulations on a 512 hydrogen atoms system at T = 900 K. (a) Equation of state close to the LLT.The vertical dashed lines highlight the rs values for which we show additional analysis. (b) Radial distribution function for thetwo rs values, i.e., rs = 1.47 (blue) and rs = 1.475 (orange). The value of the molecular fraction m (as defined in the main text)is also shown. (c) Value of maxk S (k)/N along the trajectories for the two rs values and its running average on a time windowof 375 fs. The different lines with the same color correspond to independent runs starting with different initial conditions.III. PIMD RESULTSDue to hydrogen light mass, nuclear quantum effects (NQEs) are relevant up to ∼ 3000 K and are known to shift the LLTtowards lower pressures [7]. Here we investigate whether the physical picture obtained from classical MD changes with theinclusion of NQEs. To do this, we performed constant temperature path integral Ornstein-Uhlenbeck dynamics (PIOUD) [8]using the MACE model as an energy and force driver. We ran simulations of a system of N = 512 hydrogen atoms, which issufficiently large to exclude finite-size effects as shown in the main text. For the dynamics, we used 12 replicas (beads), a timestep of 0.15 fs, and a friction term (used for controlling the temperature) equal to 0.06 fs−1. Given the larger computationaltime needed to perform the PIOUD, to accumulate more data we ran 5 independent runs for each density and temperature, eachstarting from different decorrelated initial configurations extracted from the classical trajectories. The equations of state obtainedat T = 900 K and T = 1000 K are shown in Figs. 6 and 7, respectively. We also calculated the radial distribution function, themolecular fraction, and the structure factor maximum maxk S (k), for values of the Wigner-Seitz radius rs before and after theLLT. Although with a different transition pressure, the PIOUD simulations confirm the presence of a first-order phase transitionbetween a molecular crystal and an atomic liquid at T = 900 K, followed by a liquid-liquid crossover at larger temperatures (i.e.,T = 1000 K).IV. COMPARISON OF THE MACE MODEL WITH AIMDTo validate our MLIP, we directly compared the results obtained with the MACE model and AIMD simulations. In Figs. (8-10), we report NVT results obtained at values of density and temperature roughly corresponding to p = 150 GPa taken fromRef. [5]. We performed simulations of length ∼ 2 ps with both AIMD and the MACE model on a system of 512 atoms, andcompared the structural properties of the system. The time-resolved value of maxk S (k) for each density and temperature isshown in Fig. 8 for AIMD and Fig. 9 for MACE. The two methods show good agreement and both predict solid-like structuresat T = 1050 K. Notice how this temperature matches the value of the melting line of Ref. [3] at p = 150 GPa. A comparison ofthe radial distribution function is reported in Fig. 10. At all temperatures, the MACE model shows a remarkable agreement withthe AIMD result, both in the molecular and atomic phases. Finally, in Fig. 11 we show the equations of state at different values71.45 1.46 1.47 1.48 1.49 1.50 1.51 1.52rs130140150Pressure (GPa)(a)0 2 4r (Å)0.00.51.01.52.0RDFm=0.269±0.007m=0.835±0.010(b)0 1 2Time (ps)0.00.10.20.3max k[S(k)]/N(c) rs=1.475rs=1.500T = 1000 KFIG. 7: Results of PIOUD simulations on a 512 hydrogen atoms system at T = 1000 K. (a) Equation of state close to the LLT.The vertical dashed lines highlight the rs values for which we show additional analysis. (b) Radial distribution function for thetwo rs values, i.e., rs = 1.475 (blue) and rs = 1.50 (orange). The value of the molecular fraction m (as defined in the main text)is also shown. (c) Value of maxk S (k)/N along the trajectories for the two rs values and its running average on a time windowof 375 fs. The different lines with the same color correspond to independent runs starting with different initial conditions.of temperature obtained with the final MACE MLIP and the AIMD ones obtained in Ref. [9]. Both results considered a systemof 128 hydrogen atoms. The two results are in excellent agreement, with a residual discrepancy which could be due to the shortlength of the AIMD simulations.These results further confirm the reliability of our MLIP./NFIG. 8: maxk S (k)/N resolved in time for different values of rs and temperatures, obtained from AIMD simulations for asystem of 512 atoms. The colored lines are the running averages over 40 fs.8/NFIG. 9: maxk S (k)/N resolved in time for different values of rs and temperatures, obtained from MD simulations using theMACE model for a system of 512 atoms. The colored lines are the running averages over 40 fs.90 1 2 3 4r (Å)0.51.01.52.0g(r)T = 1050K , rs = 1.46400 1 2 3 4r (Å)0.51.01.52.0g(r)T = 1150K , rs = 1.46400 1 2 3 4r (Å)0.51.01.52.0g(r)T = 1250K , rs = 1.46400 1 2 3 4r (Å)0.51.01.52.0g(r)T = 1350K , rs = 1.4640AIMDMACEFIG. 10: Comparison between AIMD and MACE for the radial distribution function g(r) at different values of rs andtemperatures.101.40 1.42 1.44 1.46 1.48 1.50 1.52 1.54rs120140160180200Pressure (GPa)T = 1000KT = 1100KT = 1200KT = 1300KT = 1400KFIG. 11: Comparison between the equations of state computed at different temperatures with AIMD (full symbols) and with thefinal MACE MLIP (empty symbols) for a system of 128 hydrogen atoms. The AIMD results are taken from [9].11[1] A. Tirelli, G. Tenti, K. Nakano, and S. Sorella, High-pressure hydrogen by machine learning and quantum Monte Carlo, Physical Review B106, L041105 (2022).[2] G. Tenti, K. Nakano, A. Tirelli, S. Sorella, and M. Casula, Principal deuterium Hugoniot via quantum Monte Carlo and ∆-learning, PhysicalReview B 110, L041107 (2024).[3] H. Niu, Y. Yang, S. Jensen, M. Holzmann, C. Pierleoni, and D. M. Ceperley, Stable solid molecular hydrogen above 900 K from amachine-learned potential trained with diffusion quantum Monte Carlo, Physical Review Letters 130, 076102 (2023).[4] D. P. Kingma and J. Ba, Adam: A method for stochastic optimization, in 3rd International Conference on Learning Representations (ICLR), San Diego, California,edited by Y. Bengio and Y. LeCun (2015).[5] V. V. Karasiev, J. Hinz, S. Hu, and S. Trickey, On the liquid–liquid phase transition of dense hydrogen, Nature 600, E12 (2021).[6] K. Momma and F. Izumi, VESTA: a three-dimensional visualization system for electronic and structural analysis, Journal of AppliedCrystallography 41, 653 (2008).[7] C. Pierleoni, M. A. Morales, G. Rillo, M. Holzmann, and D. M. Ceperley, Liquid–liquid phase transition in hydrogen by coupled elec-tron–ion Monte Carlo simulations, Proceedings of the National Academy of Sciences 113, 4953 (2016).[8] 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, J. Chem. Theory Comput. 13, 2400 (2017).[9] T. Bischoff, B. Jäckl, and M. Rupp, Hydrogen under Pressure as a Benchmark for Machine-Learning Interatomic Potentials, arXiv ,2409.13390 (2024). Hydrogen liquid-liquid transition from first principles and machine learning Abstract References