Neural Canonical Transformations for Quantum Anharmonic Solids of Lithium

Introduction.– Accurate prediction of crystal structures has long been a central focus in materials science. At low temperatures, a deep understanding of the physical properties of crystals composed of light atoms typically requires proper treatment of nuclear quantum effects with anharmonicity Cazorla and Boronat (2017); Monserrat et al. (2013); Kapil et al. (2019); Flores-Livas et al. (2020); Monacelli et al. (2021a). These effects can play a crucial role in determining the crystal structure, as seen in solid hydrogen Azadi et al. (2014); Borinaga et al. (2016); Monacelli et al. (2021b, 2023), helium Pederiva and Chester (1998); Vitali et al. (2010); Monserrat et al. (2014), and some metal hydrides Errea et al. (2013, 2014). In this study, we explore one of the most notable examples: lithium, the lightest alkali metal, where the quantum effects of nuclei are pronounced in a wide range of pressures and temperatures Guillaume et al. (2011). Although lithium behaves as a nearly free-electron metal at low pressure and adopts simple, high-symmetry cubic structures, the free energy difference between its bcc (body-centered cubic) and fcc (face-centered cubic) phases is less than $1~{}\text{meV/atom}$ Hutcheon and Needs (2019); Jerabek et al. (2022); Riffe and Christensen (2024); Phuthi et al. (2024), making precise calculations challenging. Additionally, lithium exhibits several metastable structures that further complicate experimental measurements Ackland et al. (2017). As the pressure increases, lithium exhibits complex physical behaviors, such as anomalous melting curves Guillaume et al. (2011); Schaeffer et al. (2012); Frost et al. (2019), and intriguing phase transitions from metal to semiconductor and back Matsuoka and Shimizu (2009); Lv et al. (2011); Marqués et al. (2011). Moreover, some high-pressure phases consist of large unit cells with tens or even hundreds of atoms Marqués et al. (2011); Gorelli et al. (2012); Wang et al. (2023), posing substantial challenges for both theoretical and experimental studies.

Numerical approaches to studying quantum crystals at finite temperatures include the well-established path integral molecular dynamics Martyna et al. (1999) and path integral Monte Carlo Barker (1979). In recent years, inspired by the successful application of vibrational self-consistent field theory in molecular studies Bowman (1978, 1986); Gerber and Ratner (1979), efforts have been made to extend it to study crystals Monserrat et al. (2013, 2014); Hutcheon and Needs (2019); Kapil et al. (2019). However, it relies on a Taylor expansion of the Born-Oppenheimer energy surface (BOES), and the wave function is a simple Hartree product. The stochastic self-consistent harmonic approximation (SSCHA) Errea et al. (2013, 2014); Monacelli et al. (2021a) provides an alternative by accounting for both ionic quantum and anharmonic effects without assumptions on the specific function form of the BOES. Nevertheless, it still relies on the Gaussian variational density matrices to define the quantum probability distribution. Recent developments have extended SSCHA to non-Gaussian density matrices, yet the entropy is still restricted to Gaussian approximations Siciliano et al. (2024); Monacelli et al. (2024).

In this work, we utilize the recently developed neural canonical transformations (NCT) approach Xie et al. (2022, 2023); Zhang et al. (2024), which is an ab initio variational density matrix method based on deep generative models, to study quantum lattice dynamics of lithium. NCT constructs orthogonal variational wave functions to describe phonons through a normalizing flow model Dinh et al. (2014, 2016); Wang (2018); Papamakarios (2019); Papamakarios et al. (2021); Xie et al. (2022, 2023); Saleh et al. (2023); Zhang et al. (2024). Additionally, we create a probabilistic model to describe the classical energy occupation probabilities for these phonons, allowing for accurate entropy calculations. For electronic calculations, we employ the deep potential model Wang et al. (2018); Zeng et al. (2023); Wang et al. (2023), a machine learning BOES, offering significant computational efficiency improvements over density functional theory (DFT) calculations. NCT’s key advantage is its ability to integrate quantum and anharmonic effects of nuclei into the wave functions, which facilitates the determination of the phonon spectrum. Moreover, the independently optimized phonon energy occupation probabilities enable the computation of anharmonic entropy. The NCT codes for lithium are open-sourced and publicly available myc .

The vibrational Hamiltonian of quantum solids.– Due to the substantial mass difference between electrons and nuclei, typically spanning several orders of magnitude, the Born-Oppenheimer approximation can be applied to decouple their motions and treat them independently. Under this approximation, the vibrational Hamiltonian is expressed as $H_{\text{vib}}=-\sum_{i}\frac{1}{2M}\nabla_{i}^{2}+V_{\text{el}}(\bm{R}),$ where the mass of a lithium atom is taken as $M=6.941~{}\text{u}$. The term $V_{\text{el}}(\bm{R})$ is the BOES, derived from electronic structure calculations at nuclear positions $\bm{R}$. In this work, to ensure both accuracy and computational efficiency, we use the deep potential model to fit the BOES Wang et al. (2018); Zeng et al. (2023); Wang et al. (2023); sup , which is derived from DFT calculations using the Perdew-Burke-Ernzerhof (PBE) functional.

The dynamical matrix can be derived from the Hessian of the BOES at the equilibrium position $\bm{R_{0}}$ Martin (2020); Baroni et al. (2001); Souvatzis et al. (2008, 2009); Flores-Livas et al. (2020): $C_{(i\alpha),(j\beta)}=\frac{1}{M}\frac{\mathrm{\partial}^{2}V_{\text{el}}}{\mathrm{\partial}u_{i\alpha}\mathrm{\partial}u_{j\beta}},$ where $i,j$ index the nuclei, $\alpha,\beta$ represent the Cartesian components, and the displacement coordinates are defined as $\bm{u}=\bm{R}-\bm{R_{0}}$. Diagonalizing the matrix in a supercell containing $N$ atoms yields $D=3N-3$ non-zero eigenvalues, which correspond to the number of phonon modes. The eigenvalues are related to the squares of the phonon frequencies, $\omega_{k}^{2}$   ($k=1,2,\ldots,D$), and the associated eigenvectors define the unitary transformation from displacement coordinates $\bm{u}$ to phonon coordinates $\bm{q}$. Consequently, the Hamiltonian can be expressed in phonon coordinates:

where the term $V_{\text{anh}}$ represents the anharmonic contributions, as detailed in sup . In this representation, the separation of high and low-frequency modes greatly enhances the efficiency in the following calculations.

Neural canonical transformation for variational density matrix.– The solution for a many-body system in the canonical ensemble can be obtained by minimizing the Helmholtz free energy for a variational density matrix:

where $k_{B}$ is the Boltzmann constant and $T$ is the temperature. Assuming that the phonons occupy the states $|\Psi_{\bm{n}}\rangle$ with probability $p_{\bm{n}}$, the variational density matrix can be represented in terms of these quantum states as:

where $\bm{n}=(n_{1},n_{2},\ldots,n_{D})$ indexes the energy levels of the phonons. An unbiased estimate of the anharmonic free energy for the variational density matrix Eq. 3 can be written as nested thermal and quantum expectations:

where $\Psi_{\bm{n}}(\bm{q})=\langle\bm{q}|\Psi_{\bm{n}}\rangle$ is the wave function in phonon coordinates. The symbol ${\mathbb{E}}$ is the statistical expectation, which can be estimated through sampling Zhang et al. (2024); Becca and Sorella (2017). In this Letter, the variational parameters within the energy occupation probabilities and wave functions are denoted as $\bm{\mu}$ and $\bm{\theta}$, respectively, i.e., $p_{\bm{n}}=p_{\bm{n}}({\bm{\mu}})$, $\Psi_{\bm{n}}(\bm{q})=\Psi_{\bm{n}}(\bm{\theta},\bm{q})$. These parameters can be optimized via gradient descent Kingma and Ba (2014), with $F$ as the loss function. The gradients with respect to the probabilities $\nabla_{\bm{\mu}}F$ and the wave function $\nabla_{\bm{\theta}}F$ sup can be efficiently computed using automatic differentiation Bradbury et al. (2018).

In a supercell with $D$ vibrational modes, setting a cutoff of $K$ energy levels per phonon (i.e., $n_{k}=1,2,\ldots,K$) results in an exponentially large state space of $K^{D}$. For supercells containing hundreds of atoms, directly representing the energy occupation probabilities $p_{\bm{n}}$ becomes impractical in computations. In the study of lithium, we assume that the probability distributions take a product form Martyn and Swingle (2019): $p_{\bm{n}}=\prod_{k=1}^{D}p(n_{k})$, where $p(n_{k})$ represents the probability of the $k$-th phonon occupying state $n_{k}$, and they are governed by learnable parameters $\bm{\mu}$. We have checked that an even more powerful variational autoregressive network Wu et al. (2019); Liu et al. (2021); Xie et al. (2023) does not improve results, likely due to weak coupling between phonon modes. The entropy can be evaluated as the expectation of the probabilities:

We note the nonlinear SSCHA Siciliano et al. (2024); Monacelli et al. (2024) corresponds to even further simplification of $p_{\bm{n}}$, which assumes that the entropy is given by a set of independent harmonic oscillators.

To construct variational wave functions, we apply a unitary transformation to a set of orthogonal basis states Xie et al. (2022, 2023); Zhang et al. (2024): $|\Psi_{\bm{n}}\rangle=U_{\bm{\theta}}~{}|\Phi_{\bm{n}}\rangle$, where the basis states $|\Phi_{\bm{n}}\rangle$ are chosen as the wave functions of a $D$-dimensional harmonic oscillator with frequencies $\omega_{k}$. We implement the unitary transformation $U_{\bm{\theta}}$ using a normalizing flow Dinh et al. (2014, 2016); Wang (2018); Papamakarios (2019); Papamakarios et al. (2021); Xie et al. (2022, 2023); Saleh et al. (2023); Zhang et al. (2024), which establishes a learnable bijection between the phonon coordinates $\bm{q}$ and a set of quasi-phonon coordinates $\bm{\xi}$. The bijection is represented as a smooth, reversible function $\bm{\xi}=f_{\bm{\theta}}(\bm{q})$, where $f_{\bm{\theta}}$ consists of neural networks with learnable parameters $\bm{\theta}$, specifically, a real-valued non-volume preserving network Dinh et al. (2016). Accordingly, the orthogonal variational wave functions of all energy levels can be formulated as sup ; Zhang et al. (2024):

where $\Phi_{\bm{n}}(\bm{\xi})=\langle\bm{\xi}|\Phi_{\bm{n}}\rangle$ are basis states. Notably, in the study of lithium, the computation involves approximately ten million orthogonal states for each training. The Jacobian determinant in Eq. 6 captures phonon interactions and anharmonic effects, enabling a more flexible and accurate representation. This approach remains robust when imaginary phonons appear in strong anharmonicity systems (e.g., saddle points of BOES). In such cases, we can choose the corresponding basis states with real-valued frequencies, and the flow model will automatically optimize to find the most suitable wave functions. A detailed derivation of NCT can be found in Supplemental Material sup and our previous work Zhang et al. (2024).

We can extend NCT naturally to the isothermal-isobaric ensemble, where the goal is to minimize the Gibbs free energy at a target pressure $P^{*}$, defined as:

where $\Omega$ is the system volume. From the relation $dG=dF+\Omega\sum\sigma_{\alpha\beta}~{}d\varepsilon_{\alpha\beta}$, once the parameters $\bm{\mu}$ and $\bm{\theta}$ have converged under constant volume optimization (i.e., when $dF=0$), the gradient of the Gibbs free energy to the strain $\bm{\varepsilon}$ depends only on the stress tensor $\bm{\sigma}$. The stress tensor and pressure can be calculated using the virial theorem Nielsen and Martin (1985); Martin (2020); sup . Then, we can optimize the lattice constants $\bm{a}$ through the strain tensor $\varepsilon_{\alpha\beta}=\Omega(\sigma_{\alpha\beta}-P^{*}\delta_{\alpha\beta})$, which is similar to the structure relaxation in other methods Monserrat et al. (2013); Monacelli et al. (2021a).

Anharmonic and nuclear quantum effects on stability.– At ambient temperature and pressure, lithium adopts a simple bcc structure. As the temperature decreases, experiments have demonstrated that the true ground state of lithium is fcc Ackland et al. (2017). Some calculations have revealed that the free energies of these structures are extremely close Hutcheon and Needs (2019); Jerabek et al. (2022); Riffe and Christensen (2024); Phuthi et al. (2024), highlighting the necessity of fully accounting for quantum and anharmonic effects. To investigate the influence of anharmonicity, we first conducted NCT calculations for bcc and fcc using supercells with $250$ and $256$ atoms, respectively, at a fixed volume of $19.2~{}\AA^{3}/\text{atom}$ and temperature $300~{}\text{K}$. As a comparative study, we set $f_{\bm{\theta}}$ in Eq. 6 to an identity transformation, meaning the phonon wave functions are harmonic oscillators. In this case, only the phonon occupation probabilities $p_{\bm{n}}$ in Eq. 4 were optimized.

At a lower temperature of $100~{}\text{K}$, the free energies of fcc are lower than that of bcc, as expected. However, as the temperature increases to $300~{}\text{K}$, the impact of anharmonicity becomes evident, as shown in FIG. 1 (a). In fcc, the free energy difference between the two approaches remains small, about $0.11~{}\text{meV/atom}$. In contrast, the anharmonic effects are much stronger in bcc, and the free energy difference expands to $2.67~{}\text{meV/atom}$. It is also observed that when we only optimized $p_{\bm{n}}$, the free energy of fcc is lower than that of bcc. However, when the optimization of wave functions is included, i.e., when anharmonic effects are considered, the bcc structure becomes more stable. This phenomenon suggests that bcc exhibits stronger anharmonicity than fcc, underscoring the critical role of anharmonicity in determining the stability.

The findings are further supported by the radial distribution functions (RDF) of atoms, as shown in FIG. 1 (b). The RDF for fcc exhibits only slight influence from anharmonic effects. In contrast, the RDF for bcc exhibits a smoother curve when anharmonic effects are considered, indicating a reduction in atomic localization. This behavior suggests that anharmonicity softens the system, resulting in a lower zero-point energy (ZPE) than the harmonic approximation. Further insights are provided by the phonon density of states (DOS) depicted in FIG. 1 (c), where the ZPE is determined as half the sum of all phonon frequencies per atom, as detailed in sup . Although the bcc structure is more stable at high temperatures, the numerical results reveal that the ZPE of fcc remains lower than that of bcc. This phenomenon can be explained by the differences in coordination numbers: the coordination number of bcc is $8$, while that of fcc is $12$. Hence, atoms in bcc interact less strongly with their neighbors, leading to higher quantum fluctuations and greater anharmonicity.

To gain deeper insight into the influence of nuclear quantum effects, we performed calculations on the bcc and fcc structures under constant pressure. The Gibbs free energies of both structures are extremely close Hutcheon and Needs (2019); Jerabek et al. (2022); Riffe and Christensen (2024); Phuthi et al. (2024). An error of just $1~{}\text{meV}$ could lead to a shift of more than $100~{}\text{K}$ in the transition temperature Riffe and Christensen (2024). FIG. 2 (a) shows the Gibbs free energy difference between these structures at $1~{}\text{GPa}$, with the fcc structure used as the reference. The two curves intersect at $142~{}\text{K}$, indicating a phase transition at this temperature. We also calculated the transition temperature through classical MD simulations with thermodynamic integration on the same BOES, obtaining a value of $185~{}\text{K}$. The main difference between these methods is that NCT accounts for the quantum effects, while MD simulations do not. Similar results are also observed at $0$ and $2~{}\text{GPa}$, as detailed in Supplemental Material sup , where NCT consistently predicts lower transition temperatures ($84$ and $196~{}\text{K}$) compared to MD ($144$ and $218~{}\text{K}$).

We also calculated the ionic entropy of both structures, as shown in FIG. 2 (b). Notably, the anharmonic entropy obtained from NCT is derived directly from the probabilities of energy occupations in Eq. 5, beyond the harmonic oscillator assumption. Under the harmonic oscillator assumption, the entropy of fcc is higher than that of the bcc. However, when anharmonicity is considered, the relationship is reversed. The higher entropy of the bcc structure is a key factor contributing to its stability in finite temperatures Souvatzis et al. (2008, 2009); Söderlind et al. (2012); Hutcheon and Needs (2019); Phuthi et al. (2024). Furthermore, we quantified the free energy difference arising from the anharmonic effects of entropy as $-k_{B}T~{}(S_{\text{anh}}-S_{\text{har}})$ (inset of FIG. 2 (b)), estimating it to be on the order of several meV. This underscores the critical importance of accurately incorporating anharmonic effect in the calculations.

High-pressure structural stability of lithium.– Under high pressure, lithium exhibits more complex structures and larger unit cell sizes. We first applied the NCT method to calculate the cI16 (cubic I-centered, I-43d) structure at $100~{}\text{K}$, using a supercell of $432$ atoms under various pressures. The NCT method optimizes the equilibrium positions of atoms through coordinate transformations. Subsequently, we quenched the sampled structures to their ground state with the BOES and analyzed the fractional coordinates $(x,x,x)$ of Wyckoff position $16c$ as a function of atomic volume. As depicted in FIG. 3 (a), our results are consistent with the experiment reported in Ref. Hanfland et al. (2000), demonstrating the reliability of NCT in structure optimizations.

As illustrated in FIG. 3 (b), the phonon DOS and ZPE of the cI16 and oC88 (orthorhombic C-face centered, C2mb) structures are calculated at $70~{}\text{GPa}$. Under the harmonic approximation, the ZPE of both structures are found to be comparable, consistent with the results obtained using the finite displacement and density functional perturbation theory methods Wang et al. (2023). The anharmonic effects soften the phonon spectrum of cI16, reducing the ZPE by $4.92~{}\text{meV}$ and further enhancing its stability. In contrast, the ZPE of oC88 decreases by only $1.67~{}\text{meV}$ under anharmonic effects, indicating a smaller impact compared to cI16. This result suggests that when the anharmonic effect is considered, the stability of oC88 decreases, contrary to the expectations of previous studies Marqués et al. (2011); Lv et al. (2011); Wang et al. (2023).

Additionally, we calculated the Gibbs free energies at $100~{}\text{K}$, as shown in FIG. 3 (c). The free energy of oC88 remains consistently higher than that of cI16 across all pressure conditions, which contradicts previous experiments. It has been reported that the resistivity increases sharply by more than four orders of magnitude after the cI16 phase, ultimately transforming into a semiconductor state. Compression experiments in Ref. Guillaume et al. (2011) observed the cI16-oC88 transition, as evidenced by changes in crystal color and diffraction patterns. Raman spectra measurements in Ref. Gorelli et al. (2012) detected signals corresponding to the oC88 phase. Another experiment Frost et al. (2019) also observed a phase transition around $60~{}\text{GPa}$ through the changes in diffraction peaks.

After the oC88 structure was experimentally observed Guillaume et al. (2011), theoretical studies attempted to explain its stability. However, Ref. Marqués et al. (2011) concluded from their zero-temperature calculations that oC88 was only the second most stable phase with an enthalpy about $1~{}\text{meV}$ higher than cI16, attributing this to an insufficient consideration of ZPE and thermal effects. In another work Lv et al. (2011), the authors also failed to identify oC88 as a stable structure. In contrast, Ref. Gorelli et al. (2012) found that oC88 is stable when using the PBE functional with a harmonic ZPE at $200~{}\text{K}$. However, a recent study Wang et al. (2023) demonstrated that neither harmonic nor anharmonic approximations could reproduce the results of Ref. Gorelli et al. (2012) at various conditions, showing a free energy difference at least $1~{}\text{meV}$ with oC88 consistently higher than cI16. Surprisingly, as NCT captures nuclear quantum effects and anharmonic behaviors more accurately, the difference increases to $4~{}\text{meV}$. Therefore, we strongly suspect that the instability of oC88 arises from the limited accuracy of the DFT (PBE) calculations used in fitting the BOES Wang et al. (2023). It has been observed that DFT often over-stabilizes metallic states relative to non-metallic states Ma et al. (2009, 2007).

To validate our hypothesis, we employed the NCT-optimized structures at $70~{}\text{GPa}$ and conducted single-point electronic structure calculations using the high-accuracy Heyd-Scuseria-Ernzerhof (HSE06) functional Martin (2020). The HSE functional incorporates a hybrid exchange-correlation correction, enabling a clearer distinction between metallic and non-metallic states. Notably, HSE calculations are significantly more computationally demanding, requiring approximately two orders of magnitude more computational resources than PBE. Additional details of HSE calculations are provided in sup . Our results reveal that the relative energy of oC88 compared to cI16 decreases by $6.17~{}\text{meV}$ under HSE in comparison to PBE. This reduction is significantly larger than the contributions from ZPE, anharmonic, and finite temperature effects. The HSE correction, depicted as the thick line in FIG. 3 (c), predicts a cI16-oC88 phase transition at approximately $62~{}\text{GPa}$ and $100~{}\text{K}$. This finding is consistent with experimental observations, which report a narrow stability range for the oC88 phase, existing between $62$ and $70~{}\text{GPa}$, flanked by the cI16 and oC40 phases on either side, respectively Guillaume et al. (2011). The electronic DOS of oC88, depicted in FIG. 3 (d), shows that while the HSE correction lowers the potential energy, oC88 still behaves as a poor metal.

Conclusions.– In summary, we developed the NCT method Xie et al. (2022, 2023); Zhang et al. (2024) to study anharmonic quantum solids and applied it to lithium. It enables the calculation of excited-state wave functions of nuclear motions beyond the harmonic approximation, allowing for the extraction of anharmonic phonon spectra. Additionally, the independently optimized phonon occupation probabilities facilitate the computation of anharmonic entropy. The results demonstrate that anharmonic effects play a crucial role in the structural stability. Nuclear quantum effects also introduce significant corrections to the fcc-bcc phase transition temperature. The fractional coordinates of the cI16 structure have been determined and closely align with experimental findings. Moreover, we identified that the failure of previous numerical studies Marqués et al. (2011); Lv et al. (2011); Wang et al. (2023) to observe the stability of oC88 was due to the limitations of the PBE functional in accurately describing poor metallic states. To address this, we applied the HSE functional to refine the results and estimate the stability region of oC88. Looking ahead, both experimental and computational investigations suggest that the emergence of novel high-density lithium solid phases between the cI16 and liquid phases presents a promising avenue for future exploration Frost et al. (2019); Wang et al. (2023). In this regime, high-accuracy DFT calculations and precise treatment of nuclear quantum effects may be required. Overall, the NCT method holds great potential for advancing our understanding and addressing diverse challenges in quantum crystal studies.

Acknowledgements.– We are grateful for the useful discussions with Hao Xie, Jian Lv, Xinyang Dong, Zhendong Cao, Zihang Li, Ruisi Wang, and Peize Lin. This work is supported by the National Natural Science Foundation of China under Grants No. 92270107, No. T2225018, No. 12188101, No. T2121001, No. 12134012, and No. 12374067, and the Strategic Priority Research Program of the Chinese Academy of Sciences under Grants No. XDB0500000 and No. XDB30000000.

Supplemental Material: Neural Canonical Transformations for Quantum Anharmonic Solids of Lithium