Revisiting the Broken Symmetry Phase of Solid Hydrogen: A Neural Network Variational Monte Carlo Study
Authors: Shengdu Chai, Chen Lin, Xinyang Dong, Yuqiang Li, Wanli Ouyang, Lei Wang, X. C. Xie
Organizations: Interdisciplinary Center for Theoretical Physics and Information Sciences (ICTPIS), Fudan University, Shanghai 200433, China · Shanghai Artificial Intelligence Laboratory, Shanghai 200232, China · Department of Engineering, University of Oxford, Oxford OX1 4BH, UK · Beijing National Laboratory for Condensed Matter Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China · Department of Information Engineering, The Chinese University of Hong Kong, Hong Kong SAR HKG, China · International Center for Quantum Materials, School of Physics, Peking University, Beijing 100871, China · Hefei National Laboratory, Hefei 230088, China
The crystal structure of high-pressure solid hydrogen remains a fundamental open problem. Although the research frontier has mostly shifted toward ultra-high pressure phases above 400 GPa, we show that even the broken symmetry phase observed around 130~GPa requires revisiting due to its intricate coupling of electronic and nuclear degrees of freedom. Here, we develop a first principle quantum Monte Carlo framework based on a deep neural network wave function that treats both electrons and nuclei quantum mechanically within the constant pressure ensemble. Our calculations reveal an unreported ground-state structure candidate for the broken symmetry phase with Cmcm space group symmetry, and we test its stability up to 96 atoms. The predicted structure quantitatively matches the experimental equation of state and gives the closest x-ray diffraction peak-position match among the tested candidates. Furthermore, our group-theoretical analysis provides a symmetry-counting compatibility check between the Cmcm structure and existing Raman and infrared spectroscopic data. Crucially, static density functional theory calculation reveals the Cmcm structure as a dynamically unstable saddle point on the Born-Oppenheimer potential energy surface, demonstrating that a full quantum many-body treatment of the problem is necessary. These results shed new light on the phase diagram of high-pressure hydrogen and call for further experimental verifications.
Figures & tables
Figure 1: The computational framework. The nuclear coordinates R and the electron coordinates r are fed into two neural networks that generate the nuclear wave function χ(R) and the electronic wave function φ(r,R) , respectively. Their product, Ψ(r,R)=χ(R)φ(r,R) , defines the total trial wave function. Together with the lattice parameters, they determine the enthalpy of the system G=E+PextΩ , where Pext is the external pressure and Ω the cell volume. Treating the enthalpy as the objective function, the lattice parameters L are optimized via simulated annealing, whereas the neural-network parameters θ are optimized by gradient descent.
Variational ansatz
E/N (Ha/atom)
σ2/N (Ha 2 /atom)
SJ-LDA (VMC) [ 49 ]
-0.5195(2)
SJ-LDA (DMC) [ 49 ]
-0.52415(5)
BF-PW (VMC) [ 49 ]
-0.52194(5)
0.025(1)
BF-PW (DMC) [ 49 ]
-0.52610(7)
NQS (VMC) [ 66 ]
-0.52854(9)
0.0088(1)
Present work
-0.529674(5)
0.01575(5)
Table 1: Ground-state energies for N=54 hydrogen atoms arranged into a BCC lattice with rs=1.31 treated in full quantum manner. The digits in the brackets represent the statistical uncertainty of the energy and corresponding variance. Details are provided in Sec. IV of the SM [ 1 ] .
Figure 2: The structure with the lowest enthalpy among those tested at 130 GPa for the broken symmetry phase of solid hydrogen. The figure shows the projection of H2 molecules onto the ab plane. Blue and red lines represent hydrogen molecules in the first and second layers along the c axis, respectively. The black dashed line outlines the primitive cell, while the black solid line indicates the conventional cell. 32 hydrogen atoms in the conventional cell fully occupy four 8g Wyckoff positions of the space group Cmcm .
Figure 3: Orientation order of the structure shown in Fig. 2 . The angular distribution g(θ,ϕ) illustrates that preferred orientations of H2 molecules are within the ab plane. The color bar denotes the normalized density.
Figure 4: The simulated XRD patterns for various candidate structures at 130 GPa, uniformly rescaled to the experimental volume, compared with the room-temperature observation of Ref. [ 51 ] ; the NNVMC calculation targets T=0 . For Cmcm , both the static profile and a coherent quantum-density-corrected profile obtained from ∣Ψ∣2 nuclear samples are shown. Details are provided in the SM [ 1 ] . To our knowledge, original XRD profiles (beyond reported peak positions) are unavailable in related literature, with the exception of Ref. [ 52 ] , which focuses on pressures above 200 GPa. This comparison highlights a gap that calls for further experimental validation.
Space group
IR vibron
Raman vibron
IR phonon
Raman phonon
Pca21 [ 23 ]
A1+B1+B2
A1+A2+B1+B2
2A1+2B1+2B2
2A1+3A2+2B1+2B2
P21/c [ 23 ]
Au+Bu
Ag+Bg
2Au+Bu
3Ag+3Bg
P63/m [ 23 ]
E1u
2Ag+E2g
Au+2E1u
2Ag+E1g+3E2g
Cmcm
3B2u+B3u
3Ag+B1g
3B2u+6B1u+3B3u
4B1g+7B3g+4Ag+5B2g
Experiments [ 77 , 45 , 43 ]
2
1
?
?
Table 2: Summary of theoretical predictions and experimental observations on IR and Raman modes. The comparison is based on symmetry-allowed mode counting and the observed features may form a subset of the allowed modes. Question marks (?) denote that the corresponding phonon modes are currently unreported or unresolvable. To our knowledge, the IR phonon is not reported in the literature, and the Raman phonon is reported as a single peak in [ 43 ] , while Hemley et al. [45] observed that this phonon mode becomes extremely weak.
Hyperparameter
Value
Hyperparameter
Value
Dimension of one-electron layer
256
Dimension of two-electron layer
32
Dimension of nuclear layer
32
Number of layers
3
Number of determinants
16
Optimizer
KFAC
Learning rate
3e-2
Learning rate decay
1
Learning rate delay
1e4
Damping
1e-3
Momentum of optimizer
0.0
Batch size
1024
Table S1: Hyperparameters for BCC Train
Hyperparameter
Value
Hyperparameter
Value
Dimension of one-electron layer
256
Dimension of two-electron layer
32
Dimension of nuclear layer
32
Number of layers
3
Number of determinants
8
Optimizer
KFAC
Learning rate
3e-3
Learning rate decay
1
Learning rate delay
1e4
Damping
1e-3
Momentum of optimizer
0.0
Batch size
1024
Table S2: Hyperparameters for Broken Symmetry Phase
Figure S1: Comparison of different annealing strategies. The high-noise strategy (single estimate, more steps) demonstrates faster convergence and lower variance compared to the low-noise strategy (averaged estimates), despite identical computational costs. The shadow in the figure represents the original data points, while the solid line is the smoothed curve via a moving average with a window size of 1000 steps.
Figure S2: Left : Training curves for the BCC dynamic lattice at rs=1.31 in the NVT ensemble. Right : Training curves for various initial structures at Pext=130 GPa in the NPT ensemble, highlighting the stability of the Cmcm structure. The legend denotes the initial structure. The shadow in the figure represents the original data points, while the solid line is the smoothed curve via a moving average with a window size of 1000 steps.
Figure S3: Training curves for various initial structures at Pext=130 GPa in the NPT ensemble with 2×2×2 twist averaged boundary conditions. The legend denotes the initial structure. Shaded regions represent the original data points, while solid lines show smoothed curves obtained via a moving average with a window size of 1000 steps.
Figure S4: The evolution of lattice parameters and volume per atoms as a function of training steps. The legend denotes the initial structure.
Figure S5: Comparison between two NPT training runs starting from a Cmcm structure with 64 atoms at 130 GPa: one with cell angles constrained to α=β=γ=90∘ (’Fixed angles’) and another with angles freely optimized along with other lattice parameters (’Free angles’). Shaded regions in enthalpy curves represent the original data points, while solid lines show smoothed curves obtained via a moving average with a window size of 1000 steps.
Initial Structure
Final Enthalpy (Ha/atom)
PBC
TABC
Cmcm
-0.483666(10)
-0.477739(6)
Pca21
-0.482817(9)
-0.475292(6)
P21c
-0.482795(9)
-0.475416(6)
P63/mmc with aligned molecule
-0.47791(1)
-0.45908(1)
P63m
-0.46987(2)
-0.45477(2)
Table S3: Final enthalpies based on inference with 10000 steps after training, including PBC and TABC.
Space group
Lattice parameters (Bohr, deg)
Wyckoff position
Fractional coordinates
x
y
z
Cmcm
a=8.48
8g
0.5992
0.3967
0.2500
b=9.60
8g
-0.0825
0.2402
0.2500
c=5.47
8g
0.7642
0.3967
0.2500
α=β=γ=90.0
8g
0.4175
0.0867
0.2500
Table S4: The structure identified at 130 GPa for the broken symmetry phase.
All selected peaks
Near experimental peaks
Candidate
Nall
Δsim→exp
Δexp→sim
Nnear
Δsim→exp
Δexp→sim
Cmcm
7
0.1248
0.0036
4
0.0039
0.0036
P21/c-8
9
0.1482
0.0116
6
0.0178
0.0116
P21/c-24
12
0.1447
0.0096
7
0.0286
0.0096
Pca21
8
0.1570
0.0078
5
0.0089
0.0078
P63/m
8
0.1603
0.0078
4
0.0220
0.0078
Table S5: RMS peak-position comparison after uniform volume rescaling. The left block includes all simulated peaks above the common intensity threshold; the right block is restricted to 1.35≤d≤1.65 Å. All RMS differences are in Å.
Figure S6: The structures used in the XRD patterns in Fig. 4 of the main text. The solid line represents the bonds between the hydrogen atoms, while the box represents the conventional cell for Cmcm and primitive cell for other structures.
Crystal structure prediction (CSP), which aims to predict the 3D atomic arrangement of a crystal from its composition, is central to materials discovery and mechanistic understanding. Crystal symmetry plays a crucial role in CSP, but given the composition in a unit cell, existing methods either struggle with the NP-hard combinatorial challenge of enforcing symmetry rigorously or rely on retrieving known templates, inherently limiting both physical fidelity and the discovery of genuinely new materials. To address this challenge, we introduce NextCrystal, a symmetry-driven generative framework that employs large language models to encode chemical semantics and directly generate fine-grained Wyckoff site patterns from atomic stoichiometry, eliminating reliance on database lookups. To overcome the combinatorial complexity of site assignments, we incorporate domain knowledge via an efficient, linear-complexity heuristic beam search, rigorously enforcing algebraic consistency between site multiplicities and atomic stoichiometry. By integrating this symmetry-consistent template into a diffusion backbone, the framework constrains the stochastic generative trajectory to a physically plausible geometric manifold. NextCrystal achieves state-of-the-art performance on stability, uniqueness, and novelty (SUN) benchmarks, as well as superior structural matching, establishing a rigorous paradigm for exploring previously unexplored crystallographic space without relying on prior structural templates. As a representative application, first-principles screening of HfO2 candidates generated by NextCrystal identifies a previously unreported dynamically stable Pnma phase, 0.056~eV/atom lower in energy than the conventional high-pressure Pnma phase.
Jinming Mu, Lixin He, Xudong Zhu +1
Laboratory of Quantum Information, University of Science and Technology of China, Hefei 230026, Anhui, China · Institute of Artificial Intelligence, Hefei Comprehensive National Science Center, Hefei 230088, Anhui, China
Simulating microstructure evolution requires quantum-mechanical accuracy and mesoscopic reach in length and time scales, a combination that no current method achieves. Classical phase-field models provide this reach, but their accuracy is limited by phenomenological free energies and mobilities. Here we develop a framework for learning ab initio phase-field models, where the mesoscopic equation is not postulated but derived from a Mori-Zwanzig projection of molecular dynamics onto species-density fields under explicit assumptions. The nonlocal free energy and mobility left unspecified by this equation are parametrized by neural networks and learned from short molecular dynamics trajectories generated with machine-learning interatomic potentials of ab initio accuracy. We demonstrate the framework on an iron-boron melt and on hydrogen-helium mixtures under planetary conditions. For iron-boron, the model shows that the melt at the FeB4 composition is spinodally unstable at ambient pressure but stabilized at 10 GPa, offering a thermodynamic rationale for why FeB4 has been synthesized only under high pressure. For hydrogen-helium, the model predicts the immiscibility boundary and captures droplet nucleation and growth in helium-rain simulations of a column corresponding to 2.2 million atoms, far beyond the scale of atomistic modeling at comparable accuracy. Trained across compositions and conditions, such models could provide a mesoscopic counterpart to ab initio molecular dynamics.
Mengyi Chen, Peichen Zhong, Zihan Zhang +1
Department of Mathematics, National University of Singapore, Singapore · Department of Materials Science and Engineering, National University of Singapore, Singapore · Institute for Functional Intelligent Materials, National University of Singapore, Singapore
We study the zero-temperature phase diagram of the 2D spin-imbalanced Fermi gas with short-ranged attractive interactions using the recently developed neural network variational Monte Carlo method with the AGPs FermiNet Ansatz. The Fulde-Ferrell-Larkin-Ovchinnikov phase is observed in the weakly interacting BCS limit and a polarised superfluid is seen in the strongly interacting BEC limit. When the interactions are strong, the minority-spin momentum density is reduced almost to zero in the momentum-space region occupied by the unpaired majority-spin electrons. When the interactions are very strong, phase separation occurs, with regions containing bosonic pairs and unpaired regions occupied by the remaining majority-spin particles. In addition, we observe translational symmetry breaking at intermediate interaction strengths, where the system forms an exotic crystal of Cooper pairs in a Fermi fluid of unpaired majority-spin particles. We provide a possible explanation for the formation of the crystalline phase, explain the origins of the k-space momentum-density hole when the pairs are tightly bound, and discuss how our approach opens new directions for future work.
Wan Tong Lou, Gino Cassella, Andres Perez Fadon +5
Department of Physics, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom · 2DeepMind, 6 Pancras Square, London N1C 4AG, United Kingdom · Department of Physics TQM, Technische Universität München, James-Franck-Straße 1, D-85748 Garching, Germany +1