Anomalous Thermal Transport across the Superionic Transition in Ice

Rong Qiu(邱荣)$^{1,2}$, Qiyu Zeng(曾启昱)$^{1,2}$, Han Wang(王涵)$^{3}$, Dongdong Kang(康冬冬)$^{1,2}$,\\ Xiaoxiang Yu(余晓翔)$^{1,2\hyperlink{s}{*}}$, and Jiayu Dai(戴佳钰)$^{1,2\hyperlink{s}{*}}$

doi:10.1088/0256-307X/40/11/116301

⇦Chinese Physics Letters, 2023, Vol. 40, No. 11, Article code 116301Express Letter⇨ Anomalous Thermal Transport across the Superionic Transition in Ice Rong Qiu (邱荣)1,2, Qiyu Zeng (曾启昱)1,2, Han Wang (王涵)3, Dongdong Kang (康冬冬)1,2, Xiaoxiang Yu (余晓翔)1,2*, and Jiayu Dai (戴佳钰)1,2* Affiliations 1Department of Physics, National University of Defense Technology, Changsha 410073, China 2Hunan Key Laboratory of Extreme Matter and Applications, National University of Defense Technology, Changsha 410073, China 3Laboratory of Computational Physics, Institute of Applied Physics and Computational Mathematics, Beijing 100088, China Received 19 September 2023; accepted manuscript online 27 October 2023; published online 6 November 2023 ^*Corresponding authors. Email: xxyu@nudt.edu.cn; jydai@nudt.edu.cn Citation Text: Qiu R, Zeng Q Y, Wang H et al. 2023 Chin. Phys. Lett. 40 116301 Abstract Superionic ices with highly mobile protons within stable oxygen sub-lattices occupy an important proportion of the phase diagram of ice and widely exist in the interior of icy giants and throughout the Universe. Understanding the thermal transport in superionic ice is vital for the thermal evolution of icy planets. However, it is highly challenging due to the extreme thermodynamic conditions and dynamical nature of protons, beyond the capability of the traditional lattice dynamics and empirical potential molecular dynamics approaches. By utilizing the deep potential molecular dynamics approach, we investigate the thermal conductivity of ice-VII and superionic ice-VII$''$ along the isobar of $P = 30$ GPa. A non-monotonic trend of thermal conductivity with elevated temperature is observed. Through heat flux decomposition and trajectory-based spectra analysis, we show that the thermally activated proton diffusion in ice-VII and superionic ice-VII$''$ contribute significantly to heat convection, while the broadening in vibrational energy peaks and significant softening of transverse acoustic branches lead to a reduction in heat conduction. The competition between proton diffusion and phonon scattering results in anomalous thermal transport across the superionic transition in ice. This work unravels the important role of proton diffusion in the thermal transport of high-pressure ice. Our approach provides new insights into modeling the thermal transport and atomistic dynamics in superionic materials.

DOI:10.1088/0256-307X/40/11/116301 © 2023 Chinese Physics Society Article Text As one of the most abundant substances in the Earth and Universe, ice is of vital importance from a scientific perspective and attracts wide research interests. Especially, under the interior conditions of icy moons, where pressure ranges from 2 GPa to hundreds of GPa and temperature ranges from 300 K to 4000 K, high-pressure ice phases (VII/VII$''$/X) is expected to widely exist.[1] These phases present the same body-centered cubic (BCC) oxygen sub-lattice but differ in the dynamics of hydrogen atoms (protons).[2] For the molecular crystal ice-VII, the orientation of hydrogen bonding is disordered and continually changing as in hexagonal ice, obeying the ‘ice rule’.[3] When temperature grows above thousands of Kelvin, ice-VII transforms into the superionic phase VII$''$.[4] It has been suggested that suitable conditions for superionic ice lie deep inside the watery giants Uranus and Neptune and may be common throughout the Universe.[1] VII$''$ is characterized by highly mobile hydrogen ions (protons), behaving like a liquid and moving within BCC oxygen sub-lattices. The difference in the behavior of protons can result in anomalies in thermodynamic and transport properties of ice. The occurrence and geodynamic behaviors of these high-pressure ice polymorphs (MPa–GPa range) have important effects on the thermal evolution of icy planets. On this issue, thermal conductivity serves a key role for in-depth understanding. However, despite the enormous phase transition regime and proton transfer dynamics explored in previous efforts,[1,2,5-8] it seems likely that we still know little about the thermal transport properties of dense ice across superionic transition. Existing experimental efforts have been pursued to measure the thermal conductivity of ice-VII up to 20 GPa,[9,10] but still far from the condition of superionic regime due to the limitation of experimental techniques under extreme conditions. From a theoretical point of view, the dynamical nature of protons prevents the most commonly used tool, the lattice dynamics approach, from tackling these issue. Another way to obtain thermal conductivity is molecular dynamics simulation.[11-14] However, the ab initio method requires an expensive computational cost to reach the long-time trajectories with large simulation size required for estimation of correlation function.[15] Moreover, the diverse local environments that characterize the different relevant phases of water make classical force fields unfit for an accurate simulation of their properties. Until now, the microscopic mechanism determining the thermal conductivity of superionic ice at high pressure remains unclear. Recent advances in machine-learning potential surface allow a full quantum-mechanical, ab initio treatment of the interatomic interactions efficiently.[16-22] The deep potential water model is reported to predict a phase diagram close to experiments,[23] and its following applications have demonstrated its success in estimating the transport properties of water.[24,25] Therefore, in this work, we adopt the deep-potential (DP) water model and conduct a series of deep potential molecular dynamics (DPMD) simulations to obtain the thermal conductivity of VII and VII$''$, as well as diffusion coefficient, spectral energy density, and dynamic structure factor, to understand the heat transport and to unravel the impact of mobile protons across superionic phase transition. Computational Details. The DP model was trained with the DeePMD-kit package[26,27] using diverse ice crystal and liquid phase covering from ambient condition to extreme thermodynamic state ($P=50$ GPa, $T=2400$ K). The training data were obtained from density functional theory calculations using the strongly constrained and appropriated normed (SCAN) exchange-correlation functional. More details can be found in the Supplemental Information (SI) and in Ref. [23]. With DPMD simulation, the lattice thermal conductivity $\kappa$ is obtained from the integration of the heat current autocorrelation function (HCACF), known as the Green–Kubo formula.[28,29] Since the ice-VII and ice-VII$''$ are isotropic, the $\kappa$ results are the averaged values in the $x$, $y$, and $z$ directions. We performed a series of DPMD simulations with the LAMMPS package.[30] A large supercell containing 1296 atoms is used to overcome the size effect (see Fig. S2 in the SI). The timestep was set to 0.5 fs and the Nosé–Hoover thermostat[31,32] was employed in the isothermal ensemble. For each ($P,T$) case, the equilibrium lattice constant is firstly obtained under $NPT$ ensemble with 1296 atoms. Then the simulation box is fixed to maintain the pressure of 30 GPa in the following simulation. After a thermalization stage of 20 ps under $NVT$ ensemble, the ensemble is switched into the $NVE$ ensemble to calculate the HCACF during the next 220 ps with the correlation time set to 22 ps. Each ($P,T$) case repeats 20 times, with independent initial velocity distribution to provide a representative sample for the relevant statistical analysis. In general, each ($P,T$) case requires a total MD trajectory of 4.8 ns.

Fig. 1. (a) Atomic trajectories of ice-VII and VII$''$ at different temperatures during a 20-ps long run. The oxygen and hydrogen atoms are orange and blue respectively. The cyan color is used to highlight the selected hydrogen atoms that undergo transitions from bonded states between two adjacent oxygen atoms to superionic states around different oxygen atoms. (b) Temperature-dependent thermal conductivity $\kappa$ of ice-VII and VII$''$ along the isobar of $P = 30$ GPa. The gray vertical dashed line denotes the VII–VII$''$ phase boundary obtained from previous work.[23]

Non-Monotonic Behavior of Thermal Conductivity at Elevated Temperature. The isobar of $P = 30$ GPa is chosen to investigate the thermal and proton transport of ice-VII and VII$''$. As the temperature increases from 800 K to 1600 K, our DPMD simulations reproduce the superionic phase transition of ice reported in previous works.[6,8,23] The atomic trajectories at temperatures near the phase boundary are shown in Fig. 1(a). More atomic trajectories can be seen in Fig. S3 in the SI. We can easily identify the different behaviors of protons in the oxygen sub-lattice. At low temperatures, the hydrogen atoms are bonded and only vibrate within the O–H$\cdots$O bonds. At 1200 K close to the phase boundary, the hydrogen atoms begin to initiate hopping between different O–H$\cdots$O bonds but remain bonded with oxygen atoms. The proton can migrate from an O–H$\cdots$O bond to another, leading to a fast change in the orientation of water molecules. At higher temperatures, the protons diffuse freely out of the bcc oxygen sub-lattice. Namely, the system transits into a superionic phase. Correspondingly, $\kappa$ exhibits a non-monotonic trend, as shown in Fig. 1(b). Firstly, $\kappa$ decreases from $8.65\pm 0.76$ W$\cdot$m$^{-1}\cdot$K$^{-1}$ at 800 K to a minimum value of $7.24\pm 0.77$ W$\cdot$m$^{-1}\cdot$K$^{-1}$ at 1000 K. Then a significant increase in $\kappa$ to a three-fold value of $26.22\pm 2.69$ W$\cdot$m$^{-1}\cdot$K$^{-1}$ at $T = 1200$ K is observed. As the ice transits from VII phase into VII$''$ phase, $\kappa$ gradually increases and finally reaches a plateau value of $38.99\pm 3.76$ W$\cdot$m$^{-1}\cdot$K$^{-1}$ at 1600 K. Dominant Role of Proton Diffusion in Heat Convection. The anomalous non-monotonic trend of $\kappa$ is attributed to the significant change in proton transport across the superionic transition. We extract the contributions of heat convection and conduction to $\kappa$ ($\kappa_{\rm conv}$ and $\kappa_{\rm cond}$) by decomposing the heat current into a heat convection term and a heat conduction term, respectively.[33] As depicted in Fig. 2(a), $\kappa_{\rm conv}$ shows a monotonic increasing trend with increasing temperature. At low temperatures, $\kappa_{\rm conv}$ is close to zero. At temperatures near the phase boundary, $\kappa_{\rm conv}$ increases sharply. At higher temperatures, $\kappa_{\rm conv}$ gradually converges to a plateau, similar to $\kappa$. Here $\kappa_{\rm conv}$ dominates $\kappa$ after the superionic transition, because of the fast diffusion of protons. The behavior of $\kappa_{\rm conv}$ dominates the dramatic increase and saturation of $\kappa$, highlighting the importance of proton transfer dynamics in understanding the heat transport in ice-VII and superionic VII$''$. We plot the mean square displacements (MSD) of oxygen and hydrogen atoms in Fig. S5 in the SI. The diffusion coefficient can be estimated from the slope of MSD from DPMD trajectories. Oxygen atoms have a flat curve of MSD and thus a diffusion coefficient close to zero. In comparison, the MSD curves of proton show diffusion characteristics and different slopes at different temperatures. As shown in Fig. 2(b), the diffusion coefficient of proton $D_{\rm H}$ increases by two orders of magnitude, as the temperature increases from 1000 K to 1600 K. In addition, the diffusion behavior of proton at elevated temperatures can be well described as an Arrhenius-like process, which is given as $D(T) = A e^{-E_{\rm a}/k_{\rm B} T}$, where $T$ is temperature, $A$ is a prefactor, $E_{\rm a}$ is the activation energy of the hopping mechanism leading to particle diffusion, and $k_{\rm B}$ is Boltzmann constant. By fitting the $D_{\rm H}$ with the Arrhenius model, the $E_{\rm a}$ values of ice-VII and VII$''$ are 1.77 eV and 0.36 eV, respectively. The much smaller $E_{\rm a}$ of ice-VII$''$ indicates the weaker confinement and stronger diffusion of protons in superionic phase. The temperature where the diffusion behavior of proton changes is consistent with the phase boundary of superionic transition.

Fig. 2. Temperature dependence of (a) thermal conductivity contributed by heat convection $\kappa_{\rm conv}$ and (b) diffusion coefficient along the isobar of $P = 30$ GPa. Blue and orange markers denote the results of ice-VII and superionic VII$''$, respectively. The gray vertical dashed line denotes the VII–VII$''$ phase boundary obtained from the previous work.[23]

Except for the mass diffusion, we further investigate the thermal diffusion based on the dynamic structure factor $S(k,\omega)$ (see more details in the SI). The central Rayleigh peak encodes the thermal diffusion process as the wavenumber is small enough to reach the hydrodynamic regime.[12,19] At the hydrodynamic limit, the shape of this central peak can be described by a Lorentzian function with peak width related to the thermal diffusivity $D_{\rm T}$, which gives $S(k,\omega) \propto 2D_{\rm T}k^2/[\omega^2+(D_{\rm T}k^2)^2]$. We calculate the $S(k,\omega)$ of oxygen and hydrogen atoms. By fitting the DPMD results with the hydrodynamic expression, we present the normalized Rayleigh–Brillouin triplets for both H-contributed and O-contributed $S(k,\omega)$ at a small wavenumber of $k = 0.07$ Å$^{-1}$. As presented in Fig. 3, with increasing temperature, the broadening of central Rayleigh peaks for protons is more significant than that for oxygen atoms, indicating a more significant increase in thermal diffusion of protons compared with that of oxygen atoms. The sharp increase in $\kappa$ across the superionic transition stems from the dramatically enhanced diffusion of protons and thus sharply increased thermal diffusion of protons.

Fig. 3. Normalized dynamic structure factor $S(k,\omega)$ of hydrogen and oxygen atoms at $k = 0.07$ Å$^{-1}$ along the isobar of $P = 30$ GPa.

Transverse Mode Softening across Superionic Transition. We also extract the contributions of heat conduction to $\kappa$ ($\kappa_{\rm cond}$). As presented in Fig. 4(a), $\kappa_{\rm cond}$ shows a monotonic decrease trend for ice-VII. The $1/T$ dependence of $\kappa_{\rm cond}$ is a typical characteristic that reveals the dominant role of three-phonon scatterings.[34] At low temperatures, $\kappa_{\rm cond}$ is much larger than $\kappa_{\rm conv}$, leading to a decreasing trend. Therefore, the anomalous non-monotonic trend of $\kappa$ originates from the competing mechanism of heat conduction and convection. We also note a discontinuity of $\kappa_{\rm cond}$ near the phase boundary. On the one hand, it can be attributed to the density decrease accompanied by the VII–VII$''$ phase transition (see Fig. S8 in the SI). On the other hand, a softening of the transverse acoustic modes is observed. Here we calculate the spectral energy density $C(k,\omega)$ from DPMD trajectories (see details in the SI). The $C({\boldsymbol k},\omega)$ provides information of collective vibrational modes (group velocity, lifetime) inside the complex ice polymorph. As shown in Fig. 4(b), we present the $C(k,\omega)$ of low-frequency regime ($\nu \le 20$ THz), which contains the longitudinal and transverse acoustic branches that make dominant contributions to $\kappa_{\rm cond}$. As the wavenumber approaches the center of the first Brillouin zone, the dispersion relationship exhibits a linear behavior, and the sound velocity can be extracted. For ice-VII, the sound velocities of both longitudinal and transverse acoustic branches decrease slightly with increasing temperature. For ice-VII$''$, the sound velocity of longitudinal acoustic branches maintains a slight decreasing trend while the sound velocity of transverse acoustic branches shows a sudden larger decrease across the superionic transition. A significant softening of the transverse acoustic modes across the superionic transition is observed with decreased group velocity from $0.96$ km/s to $0.74$ km/s. Moreover, as temperature increases, the vibrational energy peaks of both longitudinal and transverse acoustic branches exhibit broadening, corresponding to the reduction in phonon lifetimes.

Fig. 4. Temperature dependence of (a) thermal conductivity contributed by heat conduction $\kappa_{\rm cond}$ and (b) normalized spectra energy density $C(k,\omega)$, where the direction of wavevector ${\boldsymbol k} = (n_x \Delta k_x,0,0)$ is set to $x$ direction with the wavenumber resolution of $\Delta k_x = 2 \pi /L_x \sim 0.07$ Å$^{-1}$. The red and blue dotted lines denote the linear dispersion relationship for longitudinal and transverse acoustic branches, respectively.

Overall, by summarizing our results, we can understand the anomalous behavior of $\kappa$ in ice-VII across the superionic transition. At moderate temperatures, the lattice vibration dominates the $\kappa$ and exhibits a typical $1/T$ dependence due to the three-phonon scattering process. At $T = 1000$ K, the onset of hydrogen diffusion between O–H$\cdots$O pairs leads to an exponential increase by two orders of magnitude in the diffusion coefficient of protons. These protons hop within the oxygen-formed BCC sub-lattice and carry heat, creating a non-negligible contribution via heat convection. As temperature increases above the superionic transition threshold ($T = 1250$ K), the first-order phase transition occurs. The significantly enhanced mobility of hydrogen, combined with softening in transverse acoustic branches, leads to a saturated value above $T = 1300$ K. Under such conditions, the diffusing protons can experience a stronger thermal diffusion process as compared with oxygen atoms. In conclusion, by utilizing the newly developed DP-SCAN water model, we investigate the microscopic mechanism behind the anomalous thermal and proton transport of high-pressure ice across the superionic transition. We explain the anomalous increasing trend of $\kappa$ with elevated temperature and illustrate the important role of proton diffusion in superionic ice. To overcome the limitation of the traditional lattice dynamics approach which requires high-order force constants to correct the quasi-harmonic approximation, here we extract all the mass transport, heat transport, and collective dynamics from long-time large-scale molecular dynamics trajectories, without any assumptions or approximations made. These approaches, combined with ab initio accurate deep neural network potential energy surface model, can be applied to various complex materials, including proton-disorder ice polymorphs, proton-diffusion superionic crystal, and amorphous materials. Acknowledgments. This work was supported by the National Natural Science Foundation of China (Grant Nos. 11874424, 12122103, and 12304307), and the Science and Technology Innovation Program of Hunan Province (Grant No. 2021RC4026). References Structure and properties of two superionic ice phasesQuantum simulation of thermally-driven phase transition and oxygen K-edge x-ray absorption of high-pressure iceA Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl IonsSuperionic and Metallic States of Water and Ammonia at Giant Planet ConditionsMelting of ice under pressureSuperionic-Superionic Phase Transitions in Body-Centered Cubic

H_{2} O

IceProton dynamics and the phase diagram of dense water iceMelting Curve and Isostructural Solid Transition in Superionic IceThermal conductivity of crystalline and amorphous ices and its implications on amorphization and glassy waterThermal conductivity of compressed

H_{2}

O to 22 GPa: A test of the Leibfried-Schlömann equationNanoscale Topological Morphology Transition and Controllable Thermal Conductivity of Wrinkled Hexagonal Boron Nitride: Implications for Thermal Manipulation and ManagementHow Does van der Waals Confinement Enhance Phonon Transport?Anomalous Impact of Surface Wettability on Leidenfrost Effect at NanoscaleEffects of Localized Interface Phonons on Heat Conductivity in Ingredient Heterogeneous SolidsHeat and charge transport in H2O at ice-giant conditions from ab initio molecular dynamics simulationsDeep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum MechanicsAccuracy of Machine Learning Potential for Predictions of Multiple-Target Physical PropertiesAtomistic Mechanism of Phase Transition in Shock Compressed Gold Revealed by Deep PotentialAb initio validation on the connection between atomistic and hydrodynamic description to unravel the ion dynamics of warm dense matterLattice Thermal Conductivity of MgSiO₃ Perovskite and Post-Perovskite under Lower Mantle Conditions Calculated by Deep Potential Molecular DynamicsThree-step Formation of Diamonds in Shock-compressed Hydrocarbons: Decomposition, Species Separation, and NucleationPhase Diagram of a Deep Potential Water ModelThermal Conductivity of Water at Extreme ConditionsViscosity in water from first-principles and deep-neural-network simulationsDeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamicsDeePMD-kit v2: A software package for deep potential modelsMarkoff Random Processes and the Statistical Mechanics of Time-Dependent PhenomenaStatistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction ProblemsFast Parallel Algorithms for Short-Range Molecular DynamicsA unified formulation of the constant temperature molecular dynamics methodsCanonical dynamics: Equilibrium phase-space distributionsHybrid Thermal Transport Characteristics of Doped Organic Semiconductor Poly(3,4-ethylenedioxythiophene):TosylateThermal Conductivity of Complex Dielectric Crystals

[1]	Prakapenka V B, Holtgrewe N, Lobanov S S, and Goncharov A F 2021 Nat. Phys. 17 1233
[2]	Kang D D, Dai J Y, Sun H Y, Hou Y, and Yuan J M 2013 Sci. Rep. 3 3272
[3]	Bernal J D and Fowler R H 1933 J. Chem. Phys. 1 515
[4]	Cavazzoni C, Chiarotti G L, Scandolo S, Tosatti E, Bernasconi M, and Parrinello M 1999 Science 283 44
[5]	Schwegler E, Sharma M, Gygi F, and Galli G 2008 Proc. Natl. Acad. Sci. USA 105 14779
[6]	Hernandez J A and Caracas R 2016 Phys. Rev. Lett. 117 135503
[7]	Hernandez J A and Caracas R 2018 J. Chem. Phys. 148 214501
[8]	Queyroux J A, Hernandez J A, Weck G, Ninet S, Plisson T, Klotz S, Garbarino G, Guignot N, Mezouar M, Hanfland M, Itié J P, and Datchi F 2020 Phys. Rev. Lett. 125 195501
[9]	Andersson O and Inaba A 2005 Phys. Chem. Chem. Phys. 7 1441
[10]	Chen B, Hsieh W P, Cahill D G, Trinkle D R, and Li J 2011 Phys. Rev. B 83 132301
[11]	Qiu R, Yu X X, Wang D, Zhang S, Kang D D, and Dai J Y 2021 ACS Appl. Nano Mater. 4 10665
[12]	Yu X, Ma D, Deng C, Wan X, An M, Meng H, Li X, Huang X, and Yang N 2021 Chin. Phys. Lett. 38 014401
[13]	Wang Y, Yu X X, Wan X, Yang N, and Deng C C 2021 Chin. Phys. Lett. 38 094401
[14]	Wu M, Shi R, Qi R, Li Y, Feng T, Liu B, Yan J, Li X, Liu Z, Wang T, Wei T, Liu Z, Du J, Chen J, and Gao P 2023 Chin. Phys. Lett. 40 036801
[15]	Grasselli F, Stixrude L, and Baroni S 2020 Nat. Commun. 11 3605
[16]	Zhang L F, Han J Q, Wang H, Car R, and E W N 2018 Phys. Rev. Lett. 120 143001
[17]	Ouyang Y L, Zhang Z W, Yu C Q, He J, Yan G, and Chen J 2020 Chin. Phys. Lett. 37 126301
[18]	Chen B, Zeng Q, Wang H, Kang D, Dai J, and Ng A 2021 arXiv:2006.13136 [cond-mat.mtrl-sci]
[19]	Zeng Q Y, Yu X X, Yao Y P, Gao T Y, Chen B, Zhang S, Kang D D, Wang H, and Dai J Y 2021 Phys. Rev. Res. 3 033116
[20]	Yang F H, Zeng Q Y, Chen B, Kang D D, Zhang S, Wu J H, Yu X X, and Dai J Y 2022 Chin. Phys. Lett. 39 116301
[21]	Chen B, Zeng Q, Yu X, Chen J, Zhang S, Kang D, and Dai J 2022 arXiv:2208.01830 [astro-ph.EP]
[22]	Han J, Zeng Q, Chen K, Yu X, and Dai J 2023 Nanomaterials 13 1576
[23]	Zhang L F, Wang H, Car R, and E W N 2021 Phys. Rev. Lett. 126 236001
[24]	Zhang C Z, Puligheddu M, Zhang L F, Car R, and Galli G 2023 J. Phys. Chem. B 127 7011
[25]	Malosso C, Zhang L F, Car R, Baroni S, and Tisi D 2022 npj Comput. Mater. 8 139
[26]	Wang H, Zhang L F, Han J Q, and E W N 2018 Comput. Phys. Commun. 228 178
[27]	Zeng J Z, Zhang D, Lu D H, Mo P H, Li Z Y, Chen Y X, Rynik M, Huang L, Li Z, Shi S, Wang Y, Ye H, Tuo P, Yang J, Ding Y, Li Y, Tisi D, Zeng Q, Bao H, Xia Y, Huang J, Muraoka K, Wang Y, Chang J, Yuan F, Bore S L, Cai C, Lin Y, Wang B, Xu J, Zhu J X, Luo C, Zhang Y, Goodall R E A, Liang W, Singh A K, Yao S, Zhang J, Wentzcovitch R, Han J, Liu J, Jia W, York D M, E W N, Car R, Zhang L, and Wang H 2023 J. Chem. Phys. 159 054801
[28]	Green M S 1952 J. Chem. Phys. 20 1281
[29]	Kubo R 1957 J. Phys. Soc. Jpn. 12 570
[30]	Plimpton S 1995 J. Comput. Phys. 117 1
[31]	Nosé S 1984 J. Chem. Phys. 81 511
[32]	Hoover W G 1985 Phys. Rev. A 31 1695
[33]	Yu X X, Li R Y, Shiga T, Feng L, An M, Zhang L F, Shiomi J, and Yang N 2019 J. Phys. Chem. C 123 26735
[34]	Roufosse M and Klemens P G 1973 Phys. Rev. B 7 5379