Article

Evidence for supercritical behaviour of

high-pressure liquid hydrogen

https://doi.org/10.1038/s41586-020-2677-y Bingqing Cheng1,2,3 ✉, Guglielmo Mazzola4, Chris J. Pickard5,6 & Michele Ceriotti7,8

Received: 28 June 2019

Accepted: 10 July 2020 Hydrogen, the simplest and most abundant element in the Universe, develops a
Published online: 9 September 2020 remarkably complex behaviour upon compression1. Since Wigner predicted the
dissociation and metallization of solid hydrogen at megabar pressures almost a
Check for updates
century ago2, several efforts have been made to explain the many unusual properties
of dense hydrogen, including a rich and poorly understood solid polymorphism1,3–5,
an anomalous melting line6 and the possible transition to a superconducting state7.
Experiments at such extreme conditions are challenging and often lead to
hard-to-interpret and controversial observations, whereas theoretical investigations
are constrained by the huge computational cost of sufficiently accurate quantum
mechanical calculations. Here we present a theoretical study of the phase diagram of
dense hydrogen that uses machine learning to ‘learn’ potential-energy surfaces and
interatomic forces from reference calculations and then predict them at low
computational cost, overcoming length- and timescale limitations. We reproduce
both the re-entrant melting behaviour and the polymorphism of the solid phase.
Simulations using our machine-learning-based potentials provide evidence for a
continuous molecular-to-atomic transition in the liquid, with no first-order transition
observed above the melting line. This suggests a smooth transition between
insulating and metallic layers in giant gas planets, and reconciles existing
discrepancies between experiments as a manifestation of supercritical behaviour.

Liquid hydrogen constitutes the interior of giant planets and brown approximations16–18,21,22. Whereas early simulations gave contradic-
dwarf stars, and it is commonly assumed to undergo a first-order tory results16,19,21, the most recent calculations identify small density
phase transition between an insulating molecular fluid and a con- discontinuities below 1,500 K (refs. 17,18,20,22), which have been
ducting metallic fluid1. Understanding the nature of this liquid–liquid interpreted as the signatures of a first-order LLT.
transition (LLT) is crucial for accurately modelling the structure Even for DFT simulations, which offer a balance between computa-
and evolution of giant planets, including Jupiter, Saturn and many tional cost and efficiency, the systems studied are limited to sizes of
exoplanets8. Standard planetary models assume a sharp LLT that few hundreds of atoms and timescales of a few picoseconds16–18,21,22.
is accompanied by a discontinuity in density, and therefore give a Given the subtlety of phase-transition phenomena, it is important to
clear-cut transition between an inner metallic mantle and an outer overcome the size and timescale limitations, as well as to elucidate
insulating mantle9. the effect of the details of the electronic-structure methods on the
Probing the nature of the LLT in the laboratory faces the challenges location of the LLT, the melting line and the stabilities of the different
of creating controllable high-pressure and -temperature environ- solid phases17,18,23.
ments and of confining hydrogen specimens while making measure- To address these issues, we constructed three sets of
ments. Consequently, experimental studies have not yet reached a machine-learning potentials (MLPs), using the Behler–Parrinello
consensus on whether the LLT is a first-order or a smooth transition10. artificial neural network architecture24. The three MLPs are based
Furthermore, there are considerable discrepancies of up to 100 GPa on different electronic-structure references: DFT with the Perdew–
(see Fig. 1a) between experiments at the transition pressure of the Burke–Ernzerhof (PBE) exchange-correlation functional, DFT with the
LLT10–15 in the phase diagram. Becke88–Lee–Yang–Parr (BLYP) functional, and variational QMC20.
Given the experimental difficulties, computer simulations have More details and benchmarks are provided in Supplementary Informa-
played a fundamental role in characterizing the phase diagram of tion. The results from the three methods are in qualitative agreement.
hydrogen6,16–18, using a quantum mechanical treatment of electrons In what follows we report the results generated by the MLP based on
to describe atomic interactions. Different levels of electronic-structure PBE, and describe the other results in Supplementary Information. For
theories have been employed, ranging from the accurate quantum both solid and liquid structures of small sizes, the MLP shows excellent
Monte Carlo (QMC) methods17,19,20, to density functional theory (DFT) agreement with the underlying ab initio method. Moreover, the low
1
Department of Chemistry, University of Cambridge, Cambridge, UK. 2TCM Group, Cavendish Laboratory, University of Cambridge, Cambridge, UK. 3Trinity College, Cambridge, UK.
4
IBM Quantum, IBM Research – Zurich, Rüschlikon, Switzerland. 5Department of Materials Science and Metallurgy, University of Cambridge, Cambridge, UK. 6Advanced Institute for Materials
Research, Tohoku University, Sendai, Japan. 7Laboratory of Computational Science and Modeling, Institute of Materials, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.
8
National Centre for Computational Design and Discovery of Novel Materials (MARVEL), École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland. ✉e-mail: bc509@cam.ac.uk

Nature | Vol 585 | 10 September 2020 | 217

Article
Umax C Pmax Tm T0.5 Ts U(g cm–3) CP (J K−1 mol−1)
a b
1.0 1.24 40
3,000 Celliers et al.10 1.22
20
Knudson et al.11 1.20
350 GPa
0
Ohta et al.12 0.8
1.16 40
Zaghoo et al.14,15 1.14
20

Average order parameter

McWilliams et al.13 1.12
1.10 300 GPa 0
Temperature (K)

2,000 0.6
1.08 40
1.06
1.04 20
1.02 250 GPa 0
0.4
0.98 40
0.96
1,000 20
0.94
200 GPa
0.2 0
0.86
40
0.84 20
0 0.82 150 GPa 0
100 200 300 400 0.71
c P (GPa) 40
0.70 20
0.69 100 GPa 0
1,000 2,000 3,000
Temperature (K)
P21/c-8 C2-24 C2/c-24 Cmce-12 Cmce-4

Fig. 1 | Thermodynamic properties of high-pressure hydrogen predicted by polyamorphic solution model, respectively. The intersection between the two
the MLP based on PBE DFT. The results from MLPs based on BLYP DFT and QMC green curves, marked by a green star, is the predicted location of the critical
are shown in Supplementary Information. a, The colour scheme indicates the point of the LLT. The experimental results are taken from refs. 10–15. b, The purple
molecular fraction defined by the order parameter. The black curve is the curves show the density isobar, and the orange curves show the molar heat
estimated solid–liquid coexistence line, with the upper and the lower bound of capacity at different pressures. The shaded regions indicate the conditions under
hysteresis indicated by the error bars. The density (ρ) and molar heat capacity which solid phases are stable, corresponding to the solid–liquid coexistence line
(CP) maxima at different pressures are indicated by purple and orange dots, shown in a. Error bars indicate statistical uncertainties. c, At each given pressure
respectively. The dashed and dotted green curves are the coexistence and (black lines), the crystalline structure, the space group and the size of the
phase-separation lines of atomic and molecular fluids predicted by the primitive cell of the solid hydrogen phase with the lowest enthalpy are shown.

cost of the MLP allowed us to investigate hydrogen phase transitions

using large system sizes and long simulation times. If performed using Liquid–liquid transition
DFT, the total computational cost of this study would have required We performed MD simulations across a broad range of temperatures
several hundred millions of CPU years, exceeding the capacity of the and pressures using 1,728-atom simulation cells and 0.4 ns simula-
world’s fastest supercomputers. tion time for each run. We employed an order parameter defined as
the fraction of atoms that have one neighbour within a smooth cutoff
function that is equal to 1 up to 0.8 Å and decays to 0 at 1.1 Å. As evident
Solid–liquid transition from Fig. 1a, the molecular fraction varies smoothly across the liquid
Solid hydrogen exhibits complex polymorphism, and only a few of phase diagram, with the transition region becoming narrower at low
its crystal structures and phase boundaries have been characterized temperature (T) and high pressure (P). Other observables, including
conclusively1,3–5,25,26. Encouragingly, the MLP correctly captures the the density (ρ), the molar heat capacity CP (Fig. 1c), the pair correla-
ground-state crystals (Fig. 1c) during random searches27 over a wide tion function and the electronic density of states (see Supplementary
pressure range between 100 GPa and 400 GPa (see details in Sup- Information), also show the absence of discontinuities. Both ρ and CP
plementary Information). Furthermore, it reproduces the subtle exhibit anomalous behaviours, namely, smooth peaks that become
enthalpy differences of a few millielectronvolts between the competing sharper at higher pressures (Fig. 1b). The loci of these maxima, as
polymorphs28,29, consistently with the underlying DFT reference. The well as the atomic–molecular transition region, converge towards the
accuracy of these predictions demonstrates the immense promise of melting line at above 350 GPa.
MLPs for crystal structure discovery.
To estimate the melting curve without prior knowledge of the
solid-phase diagram, we applied a simple hysteresis method: Polyamorphic solution model
1,728-atom hydrogen systems were first cooled from the liquid phase The MD simulation results (Fig. 1a) suggest the lack of a sharp LLT. To
until solidification, and subsequently reheated until melting, for a total provide a more quantitative analysis, we determined the parameters
of 0.8 ns of molecular dynamics (MD) simulation time. At each pressure of a polyamorphic solution model22,31 that describes a mixture of two
we performed eight simulations. Owing to the presence of nuclea- interconvertible liquid states. At each thermodynamic state point, the
tion barriers, the freezing and melting temperatures are affected by regular-solution molar free energy g(x) as a function of the molecular
hysteresis. Therefore, the true melting point Tm of the system lies fraction x is
between these two temperatures and can be estimated using their
mean value. The shape of the estimated melting line (the black curve in g (x) = x Δg + kBTx ln x + kBT (1 − x)ln(1 − x) + ωx (1 − x ). (1)
Fig. 1a), which has a peak at around 125 GPa and 850 K and then declines
at higher pressures, is in excellent agreement with recent experimental The term Δg = gm − ga is the chemical potential difference between the
measurements30 and previous PBE results6. atomic and the molecular phases, ω is an enthalpic term that accounts

218 | Nature | Vol 585 | 10 September 2020

 a b
0.5
3

0
100 GPa
2

Z/kB (×103 K)
0.5

1
0
150 GPa

0
g (K mol–1)

0.5

0
200 GPa –1

c
0.5 100
3 150
0 200
250 GPa
250
2

Δg/kB (×103 K)
300
0.5
350

0 1
300 GPa

0.5 Fit 0

0 Simulation
350 GPa –1
0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1,000 2,000 3,000
Molecular fraction, x Temperature (K)

Fig. 2 | Polyamorphic solution model fits of the high-pressure hydrogen increases. b, The dots represent values of Δg = gm − ga fitted to the solution
system. a, The dots show the computed Gibbs free energy profiles g(x) as model, and the lines are linear fits to Δg. c, The dots are the individual values of
functions of the molecular fraction order parameter. The results are from one ω obtained from fitting g(x) to the solution model at different pressures and
of the eight sets of metadynamics simulations. The smooth curves show the temperatures, and the curves are fits to those values. The dotted green line
individual fits to the polyamorphic solution model. The series correspond to corresponds to ω = 2T, which corresponds to the phase-separation line, and the
results obtained at T = 600, 800, 900, 1,000, 1,200, 1,500, 1,700, 2,000, 2,500 dashed line to Δg = 0, that is, the coexistence line. The error bars were
and 3,000 K, plotted in shades of red, from dark to bright. As expected, the estimated from the error of the mean for the eight sets of simulations.
minimum of g(x) shifts to lower molecular fractions as the temperature

for the non-ideality of mixing and kB is the Boltzmann constant. To lines, respectively. The two lines cross at the critical point (marked
obtain a free-energy profile that can be compared to equation (1) from by a green star) of the fluid–fluid phase transition, which is located at
simulations, we performed a set of metadynamics32 simulations in (Pc, Tc) ≈ (350 ± 1 GPa, 416 ± 2 K), coinciding approximately with the melt-
which we enhanced the spontaneous fluctuations of the order param- ing line. At T > Tc the system exhibits supercritical behaviour, without
eter x (see Supplementary Information). phase separation and with anomalies in the thermodynamic properties
As shown in Fig. 2a, each g(x) obtained from simulations has a single of the mixture following different Widom lines that emanate from the
minimum and matches perfectly the form predicted by the solution critical point. At T < Tc, the system shows a first-order phase transition,
model. This indicates perfect mixing of the two liquids and absence and the T0.5 coexistence line becomes the phase boundary. Because for
of an LLT throughout the range of temperatures and pressures that this system Tc approximately coincides with Tm, no sharp LLT can be
we explored. In addition, we used the simple empirical expressions observed. Instead, the anomalous behaviours induced by this hidden
Δg = a0 + a1P + a2T + a3PT and ω = b0 + b1P + b2/T + b3P2 to describe the critical point can be observed throughout the liquid phase diagram,
model parameters Δg and ω. As shown in Fig. 2b, c, these expressions much like the case of water33.
agree well with the Δg and ω obtained by independent fits at each state Our observation of the supercritical hydrogen fluid above the
point, at temperature and pressure conditions above the melting line. melting line contradicts several recent DFT and QMC simulations17,18,20,
The small discrepancy at low temperatures and low pressures is caused which showed a sharp LLT, suggested by small discontinuities in
by solidification, which the solution model does not consider. We then density up to around 1,000–1,500 K. The probable origin of this
used the analytic expressions to estimate the x = 0.5 coexistence line T0.5 discrepancy, which we traced to finite-size effects on the solid–liq-
(that is, the temperature at which atomic and molecular fluids become uid transition, is discussed in detail in Supplementary Information.
equally stable, determined by Δg(P, T) = 0)) and the phase-separation We performed explicit DFT MD simulations with a system size of
line Ts (that is, the temperature below which the two fluids start demix- 128 atoms. We reproduced the pressure–density relations of the
ing, determined as Ts = ω(P, Ts)/2kB). We note that being at T < Ts is not previous studies17,18, albeit using a coarser density grid along each
sufficient to observe the demixing behaviour; it is also necessary that isotherm. However, we observed the formation of solid phases
the atomic and molecular phases are equally stable, because the two at temperatures up to T = 1,250 K in constant-pressure simulations
phases are interconvertible. It is also worth noting that, although there and the appearance of defective solids up to T = 1,000 K in constant-
are different ways to define the molecular fraction x using an order volume simulations. Such solidification-like events caused disconti-
parameter, the curves Ts and T0.5 are rather insensitive to such defini- nuities in the density and radial-distribution functions, similar to a
tions. These two curves are plotted in Fig. 1a as dashed and dotted green sharp LLT.

Nature | Vol 585 | 10 September 2020 | 219

Article
The results above are based on the reference PBE potential energy 1. McMahon, J. M., Morales, M. A., Pierleoni, C. & Ceperley, D. M. The properties of
hydrogen and helium under extreme conditions. Rev. Mod. Phys. 84, 1607–1653 (2012).
surface, but our conclusion concerning supercriticality is insensitive to 2. Wigner, E. & Huntington, H. B. On the possibility of a metallic modification of hydrogen.
the electronic-structure method used to compute the atomic interac- J. Chem. Phys. 3, 764–770 (1935).
tions. The analogous results based on the MLPs trained on BLYP DFT 3. Howie, R. T., Guillaume, C. L., Scheler, T., Goncharov, A. F. & Gregoryanz, E. Mixed
molecular and atomic phase of dense hydrogen. Phys. Rev. Lett. 108, 125501 (2012).
and variational QMC datasets are shown in Supplementary Informa- 4. Zha, C., Liu, Z., Ahart, M., Boehler, R. & Hemley, R. J. High-pressure measurements of
tion. The different methods used to obtain the electronic structure hydrogen phase IV using synchrotron infrared spectroscopy. Phys. Rev. Lett. 110, 217402
lead to different locations of the coexistence and phase-separation (2013).
5. Dalladay-Simpson, P., Howie, R. T. & Gregoryanz, E. Evidence for a new phase of dense
lines, but the qualitative picture remains the same. In all cases, the hydrogen above 325 gigapascals. Nature 529, 63–67 (2016).
critical point is below the melting line. Our conclusion is reinforced 6. Bonev, S. A., Schwegler, E., Ogitsu, T. & Galli, G. A quantum fluid of metallic hydrogen
by the analysis based on the polyamorphic solution model. Although suggested by first-principles calculations. Nature 431, 669–672 (2004).
7. Ashcroft, N. W. Metallic hydrogen: a high-temperature superconductor? Phys. Rev. Lett.
a change in the electronic structure method, the inclusion of nuclear 21, 1748–1749 (1968).
quantum effects, and the residual difference between the MLP and 8. Guillot, T. The interiors of giant planets: models and outstanding questions. Annu. Rev.
the reference method can influence the details of the phase diagram, Earth Planet. Sci. 33, 493–530 (2005).
9. Hubbard, W. B. & Militzer, B. A preliminary Jupiter model. Astrophys. J. 820, 80 (2016).
the solution model shows that in the region that is usually proposed 10. Celliers, P. M. et al. Insulator-metal transition in dense fluid deuterium. Science 361,
for the critical point, the system is far from critical. At P < 200 GPa, 677–682 (2018).
T0.5 is above Ts, which means that the atomic fluid is unstable at the 11. Knudson, M. D. et al. Direct observation of an abrupt insulator-to-metal transition in dense
liquid deuterium. Science 348, 1455–1460 (2015).
conditions at which phase separation can happen. At about P > 200 GPa, 12. Ohta, K. et al. Phase boundary of hot dense fluid hydrogen. Sci. Rep. 5, 16560 (2015).
ω is negative at T ≳ 800 K, suggesting that mixing is enthalpically—and 13. McWilliams, R. S., Dalton, D. A., Mahmood, M. F. & Goncharov, A. F. Optical properties of
not only entropically—favourable. fluid hydrogen at the transition to a conducting state. Phys. Rev. Lett. 116, 255501 (2016).
14. Zaghoo, M., Salamat, A. & Silvera, I. F. Evidence of a first-order phase transition to metallic
The predicted supercritical behaviour of fluid hydrogen can hydrogen. Phys. Rev. B 93, 155128 (2016).
explain the discrepancies between different experiments. If the LLT 15. Zaghoo, M. & Silvera, I. F. Conductivity and dissociation in liquid metallic hydrogen and
is indeed first-order, all observables should undergo an abrupt change implications for planetary interiors. Proc. Natl Acad. Sci. USA 114, 11873–11877 (2017).
16. Scandolo, S. Liquid–liquid phase transition in compressed hydrogen from first-principles
when crossing the coexistence line. Instead, the supercriticality of simulations. Proc. Natl Acad. Sci. USA 100, 3051–3053 (2003).
fluid hydrogen means that the boundary of the LLT is blurred and its 17. Morales, M. A., Pierleoni, C., Schwegler, E. & Ceperley, D. M. Evidence for a first-order
liquid–liquid transition in high-pressure hydrogen from ab initio simulations. Proc. Natl
location depends on the specific criterion used to define it. In other
Acad. Sci. USA 107, 12799–12803 (2010).
words, different observables may exhibit anomalous behaviours that 18. Lorenzen, W., Holst, B. & Redmer, R. First-order liquid–liquid phase transition in dense
follow different Widom lines, as we observed for the density and heat hydrogen. Phys. Rev. B 82, 195107 (2010).
19. Delaney, K. T., Pierleoni, C. & Ceperley, D. M. Quantum Monte Carlo simulation of the
capacity. Indeed, the LLT boundaries measured by different teams at
high-pressure molecular-atomic crossover in fluid hydrogen. Phys. Rev. Lett. 97, 235702
1,000 < T < 2,000 K all qualitatively extrapolate towards the proposed (2006).
critical point10–15 (see Fig. 1a). The observation of a sharper transition 20. Mazzola, G., Helled, R. & Sorella, S. Phase diagram of hydrogen and a hydrogen–helium
mixture at planetary conditions by quantum Monte Carlo simulations. Phys. Rev. Lett. 120,
in the low-temperature compression experiments of Knudson et al.11,
025701 (2018).
in comparison to those performed by Celliers et al.10, is also consistent 21. Vorberger, J., Tamblyn, I., Militzer, B. & Bonev, S. A. Hydrogen–helium mixtures in the
with supercritical behaviour. interiors of giant planets. Phys. Rev. B 75, (2007).
22. Geng, H. Y., Wu, Q., Marqués, M. & Ackland, G. J. Thermodynamic anomalies and three
The polyamorphic solution model, which we validated in our
distinct liquid–liquid transitions in warm dense liquid hydrogen. Phys. Rev. B 100, 134109
simulations by combining an MLP trained on electronic-structure (2019).
calculations and thorough statistical sampling, models the stability 23. Clay, R. C. III et al. Benchmarking exchange-correlation functionals for hydrogen at high
pressures using quantum Monte Carlo. Phys. Rev. B 89, 184106 (2014).
and miscibility of atomic and molecular hydrogen and quantitatively
24. Behler, J. & Parrinello, M. Generalized neural network representation of high-dimensional
describes the molecular-to-atomic transition in dense liquid hydrogen. potential energy surfaces. Phys. Rev. Lett. 98, 146401 (2007).
This model thus provides a thermodynamic picture of the LLT, and 25. Magdău, I. B., Marqués, M., Borgulya, B. & Ackland, G. J. Simple thermodynamic model
for the hydrogen phase diagram. Phys. Rev. B 95, 094107 (2017).
can be directly employed to interpret experiments and astrophysi- 26. Loubeyre, P., Occelli, F. & Dumas, P. Synchrotron infrared spectroscopic evidence of the
cal observations. Our general approach can be used to quantitatively probable transition to metal hydrogen. Nature 577, 631–635 (2020).
assess the properties of mixtures of hydrogen and heavier elements, 27. Pickard, C. J. & Needs, R. J. Ab initio random structure searching. J. Phys. Condens. Matter
23, 053201 (2011).
as well as to address the long-standing questions concerning Jupiter’s 28. Pickard, C. J. & Needs, R. J. Structure of phase III of solid hydrogen. Nat. Phys. 3, 473
core solubility and the anomalous-luminosity problem of Saturn8. (2007).
The accuracy of the MLP for the solid phases demonstrates enormous 29. Monserrat, B., Needs, R. J., Gregoryanz, E. & Pickard, C. J. Hexagonal structure of phase III
of solid hydrogen. Phys. Rev. B 94, 134101 (2016).
promise to answer the many open questions concerning the solid-phase 30. Zha, C., Liu, H., Tse John, S. & Hemley, R. J. Melting and high P–T transitions of hydrogen
diagram of dense hydrogen26. up to 300 GPa. Phys. Rev. Lett. 119, 075302 (2017).
31. Anisimov, M. A. et al. Thermodynamics of fluid polyamorphism. Phys. Rev. X 8, 011004
(2018).
32. Laio, A. & Parrinello, M. Escaping free energy minima. Proc. Natl Acad. Sci. USA 99,
Online content 12562–12566 (2002).
Any methods, additional references, Nature Research reporting sum- 33. Soper, A. K. & Ricci, M. A. Structures of high-density and low-density water. Phys. Rev.
Lett. 84, 2881 (2000).
maries, source data, extended data, supplementary information,
acknowledgements, peer review information; details of author con- Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
tributions and competing interests; and statements of data and code published maps and institutional affiliations.

availability are available at https://doi.org/10.1038/s41586-020-2677-y. © The Author(s), under exclusive licence to Springer Nature Limited 2020

220 | Nature | Vol 585 | 10 September 2020

 s957. G.M. acknowledges financial support from the Swiss National Science Foundation
Data availability through grant number 200021-179312. C.J.P. is supported by the Royal Society through a
Royal Society Wolfson Research Merit award and the EPSRC through grant EP/P022596/
The data supporting the findings of this study are available within the 1. M.C. acknowledges funding from the Swiss National Science Foundation (project 200021-
paper, and all input files that are necessary to reproduce the reported 182057).

results are included in Supplementary Information. All data gener-

Author contributions B.C., G.M. and M.C. conceptualized the research; B.C., C.J.P. and M.C.
ated for the study are available upon request from the corresponding performed the research and analysed the data; B.C., G.M., C.J.P. and M.C. wrote the paper.
author, and the MLP for hydrogen constructed here are available at
https://github.com/BingqingCheng/MLP-highP-H. Competing interests The authors declare no competing interests.

Additional information
Acknowledgements We are thankful to G. Ackland, H. Geng. and R. Redmer, who shared their Supplementary information is available for this paper at https://doi.org/10.1038/s41586-020-
AIMD trajectories for us to benchmark the MLP. We thank S. Sorella for providing the VMC 2677-y.
training dataset. We acknowledge D. Frenkel, B. Monserrat, M. Casula, A. M. Saitta, R. Helled, G. Correspondence and requests for materials should be addressed to B.C.
Carleo and S. Sorella for discussions. B.C. acknowledges funding from the Swiss National Peer review information Nature thanks Graeme Ackland and the other, anonymous,
Science Foundation (project P2ELP2-184408), resources provided by the Cambridge Tier-2 reviewer(s) for their contribution to the peer review of this work.
system funded by EPSRC Tier-2 capital grant EP/P020259/1 and by CSCS under project ID Reprints and permissions information is available at http://www.nature.com/reprints.