1. Introduction

Perovskite materials, as prominent thermoelectric materials, excel in solar cells, microchips and thermoelectric applications due to their excellent photoelectric and tunable thermal properties^[1-7]. Antiperovskites, the electrically inverted counterparts of perovskites, remain relatively underexplored in materials science, especially regarding their thermal properties. Recently, a new family of antiperovskite materials, double antiperovskites, has shown significant potential in energy-related applications. For example, X₆SOA₂ (X = Na, K; A = Cl, Br and I) exhibits stable zT (dimensionless figure of merit) over a wide temperature range, indicating outstanding thermoelectric performance suitable for energy storage and harvesting applications^[8,9]. Li₆OS(BH₄)₂ demonstrates stable cyclic performance and long lifespan, making it a promising novel solid-state electrolyte^[10]. Simulations reveal that BiPCa₃^[11] and Mg₃BiN^[12] possess high zT values at elevated temperatures with good thermal stability, rendering them promising candidates for thermoelectric devices. Experimental measurements show that Ca₃AO (A = Pb, Sn) also exhibits great potential for thermoelectric applications due to its large thermoelectric power and low electrical resistivity^[13]. IIt has been demonstrated that both perovskite materials including oxide perovskite SrTiO₃^[14], double perovskites Cs₂NaYbCl₆^[15] and Cs₂AgBiBr₆^[16], whose strong lattice anharmonicity mainly originates from octahedral rotations and tilting, loosely bonded cations, or strong four-phonon scattering, and antiperovskite materials such as Li₃OCl^[17,18], Rb₃X(Se & Te)I^[19], and X(Ba & Sr)₃BiN^[20], where strong anharmonicity mainly arises from weak bonding, exhibit pronounced anharmonicity due to their unique crystal structures. Given their similar lattice structures, double antiperovskite materials are expected to exhibit strong anharmonicity as well. While extensive research has been devoted to thermoelectric materials for optimizing device performance^[21,22], fundamental studies on the thermal properties of double antiperovskites remain scarce.

Many theoretical and experimental studies have been conducted to understand the thermal transport properties of strongly anharmonic crystals. Zhao et al.^[23] investigated thermal transport in A₃BO (A = K, Rb; B = Br, Au) using first-principles calculations combined with the Boltzmann transport equation (BTE), revealing strong quartic phonon scattering and a resultant weak temperature dependence of lattice thermal conductivity. Wang et al.^[24] employed first-principles-based machine-learning potentials (MLP), the Wigner transport equation (WTE), and experimental measurements to study the thermal properties of halide double perovskites Cs₂AgB(II)X₅ (B(II): Pd or Pt; X: Cl or Br), demonstrating that unique atomic vibrations and atom rattling contribute to their strong anharmonicity. Liu et al.^[25] examined thermal transport in Cu₆Te₃S through density functional theory (DFT) calculations and identified glass-like thermal conductivity caused by anharmonic vibrations of Cu atoms and overall weak bonding. In theoretical computational studies of strongly anharmonic materials, lattice dynamics calculations combined with the BTE or WTE based on first-principles are commonly used to obtain mode-resolved phonon transport properties. Although replacing first-principles calculations with MLP partially enhances the capacity to handle complex crystals, challenges remain due to the lower symmetry of strongly anharmonic materials, the breakdown of perturbation theory, and the cumbersome process of obtaining temperature-dependent high-order interatomic force constants (IFCs) and solving phonon scattering rates beyond three-phonon processes. Molecular dynamics (MD) simulations using MLP derived from first-principles calculations can capture full lattice anharmonicity while balancing computational accuracy and efficiency. These approaches have been applied extensively to thermal property calculations of materials such as metal-organic frameworks^[26-28], water^[29-33], and graphene^[34-39].

Since the introduction of neural network potentials in 2007^[40], MLPs have achieved great success in materials science. Leveraging data from first-principles calculations, methods including Gaussian approximation potential^[41], moment tensor potential^[42,43], neuroevolution potential (NEP)^[44,45], recursive embedded-atom neural network^[46,47] and various others employ diverse atomic environment descriptors and machine learning models to achieve both high accuracy and fast computational speed. Recently, the integration of high-performance graphics processing units (GPUs) has further accelerated MD simulations. For example, the NEP within the GPUMD package achieves a computational speed of approximately 1 × 10⁷ atom steps per second on a single V100 GPU, outperforming other MLPs^[26,48] and providing a powerful tool for investigating thermal transport properties in complex crystals^[49].

In this work, we develop accurate machine-learning NEPs for the recently synthesized stable double antiperovskites Li₆NII₂ and Li₆NBrBr₂^[50] using an active learning (AL) strategy^[49] combined with DFT calculations. By incorporating thermal expansion and phonon anharmonic renormalization effects, we find that both harmonic and temperature-dependent phonon dispersion relations exhibit a significant portion of optical phonons with imaginary frequencies, indicating the failure of perturbation theory in describing phonon transport properties of these double antiperovskites. The calculated acoustic phonon branches of both materials at 300 K across the entire Brillouin zone remain below 2.5 THz, demonstrating low phonon energies and ultralow acoustic phonon group velocities. Large-scale MD simulations with long runtimes confirm the thermodynamic stability of Li₆NII₂ and Li₆NBrBr₂ and reveal confined Li-ion movements among lattice points within a single conventional unit cell. Thermal transport in these materials is further studied using both homogeneous nonequilibrium molecular dynamics (HNEMD) and EMD methods. Ultralow lattice thermal conductivities with glass-like temperature-dependent behavior are observed for both materials. Detailed phonon spectral analysis confirms that these thermal transport characteristics arise from the confined Li-ion motions, suppressed contributions, and weak temperature dependence of acoustic phonons.

2. Computational Models and Details

2.1 Construction of neuroevolution potentials with active learning

Previous studies have demonstrated the excellent performance of the NEP approach in investigating thermal transport in complex perovskite materials. In this work, we construct accurate NEPs for both Li₆NII₂ and Li₆NBrBr₂ based on first-principles calculations using an AL strategy^[49]. DFT calculations were performed with the projector augmented wave (PAW) method^[51] and the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation functional^[52], implemented in the Vienna Ab initio Simulation Package (VASP)^[53]. Initially, ab initio molecular dynamics trajectories of a 2 × 2 × 2 primitive cells for both materials were obtained at 300 K using VASP^[54] under an isothermal-isobaric (NPT) ensemble with a timestep of 1 fs and 10,000 steps. Only the Г point was used for k-point integration, and the energy cutoff was set to 500 eV (convergence tests for the cutoff energy are presented in Figure S1). Configurations for the initial training dataset were sampled every 25 steps, resulting in a total of 400 configurations. Additional configurations were obtained through independent MD simulations at various temperatures using the AL strategy^[55,56] with a timestep of 1 fs for both Li₆NII₂ and Li₆NBrBr₂. Furthermore, configurations subjected to strains ranging from -3% to 3% were also included following the same AL procedure implemented in the PYNEP framework^[49]. The selection distance threshold for new configurations was set at 0.025. This select-and-train process was iterated until no new structures were selected after 2 ns of MD simulation. To further enhance the stability of the NEPs, manually generated uniaxial strain, biaxial strain, and rattled configurations were added to the training set. Ultimately, 356 and 320 configurations were selected to train the NEPs for Li₆NII₂ and Li₆NBrBr₂, respectively (the distribution of sampled structures obtained via PYNEP is shown in Figure S2). Accurate virials, energies, and forces for these configurations were computed via precise DFT calculations using a 7 × 7 × 7 k-point mesh and a total energy convergence criterion of 10^-8 eV. The radial and angular cutoffs and other hyperparameters used for NEP training are summarized in Table S1.

Our NEPs for both Li₆NII₂ and Li₆NBrBr₂ exhibit excellent agreement with DFT results in predicting energies, atomic forces, and virials for both training and testing datasets, as illustrated in Figure 1, Figure S3 and Figure S4, respectively. For the Li₆NII₂ training dataset, the root mean square errors (RMSE) for energies, atomic forces, and virials are 0.39 meV, 14.41 meV/Å and 2.64 meV per atom, respectively. The convergence trends of each term in the loss function as well as phonon dispersion validation are presented in Figure S5 and Figure S6.

Display Full Size

Figure 1. For Li₆NII₂, (a) Energy; (b) Force; (c) Virial calculated by NEP as compared to the relevant results calculated from DFT in the training dataset. The solid lines in (a), (b) and (c) are the identity function applied to guide the eyes. NEP: neuroevolution potential; RMSE: root mean square errors.

2.2 Crystal structures and phonon dispersion relations

Determining temperature-dependent lattice constants is fundamental for accurately describing thermal transport in strongly anharmonic materials. In this study, we extracted the temperature-dependent lattice constants of both Li₆NII₂ and Li₆NBrBr₂ using our NEPs combined with NPT (MD simulations. The calculated lattice constants forLi₆NII₂ and Li₆NBrBr₂ are 9.714 Å and 9.131 Å, respectively, which are close to the experimental values of 9.536 Å and 8.925 Å at room temperature^[50]. The primitive and conventional unit cells of both materials are presented in Figure 2a and Figure S7, respectively, while the temperature-dependent lattice constants are shown in Figure 2b and Figure S8.

Display Full Size

Figure 2. (a) Primitive cell of Li₆NII₂ or Li₆NBrBr₂; (b) The temperature-dependent lattice constants of Li₆NII₂ and Li₆NBrBr₂.

Figure 3a and Figure 3b depict the harmonic phonon dispersion relations of Li₆NII₂ and Li₆NBrBr₂, calculated using the finite displacement method and DFT. A 5 × 5 × 5 k-point mesh, an energy convergence criterion of 10^-8 eV, and 2 × 2 × 2 supercells were employed. Large portions of both acoustic and optical phonons with imaginary frequencies appear in the harmonic phonon dispersions. Anharmonic renormalized IFCs at finite temperatures were also obtained using the temperature-dependent effective potential method^[57] and trajectories from isothermal (NVT) ensemble MD simulations with NEPs for both materials (details are provided in Figure S9). The resulting anharmonic phonon dispersion relations are shown in Figure 3c and Figure 3d, revealing clear temperature dependencies. Although acoustic phonon branches become stable at finite temperatures, optical phonon branches exhibit large imaginary frequencies (approximately -6 THz) across the entire Brillouin zone, indicating the failure of perturbation theory to accurately describe phonon transport in both Li₆NII₂ and Li₆NBrBr₂.

Display Full Size

Figure 3. Harmonic phonon dispersion of (a) Li₆NII₂ and (b) Li₆NBrBr₂; The temperature-dependent phonon dispersion of (c) Li₆NII₂ from 200 to 500 K and (d) Li₆NBrBr₂ from 100 to 400 K.

To clarify the origin of optical phonons with imaginary frequencies, we further analyzed the vibrational patterns of atoms in both materials. The atomic vibrational modes corresponding to the optical phonons with imaginary frequencies are shown in Figure 4a, Figure 4b, and Figure S10. Notably, only Li atoms participate in these imaginary phonon modes for both compounds. Additionally, atomic trajectories from NVT MD simulations are presented in Figure 4c and Figure S11. Unlike the vibrations of N, Br, and I atoms near their equilibrium positions, Li atoms exhibit stochastic movements among adjacent Li lattice sites confined within a single octahedron. These stochastic Li-ion motions invalidate perturbation theory, suggesting that MD simulations are more appropriate for calculating the thermal conductivities of such materials. Accordingly, we evaluate thermal conductivity using EMD and HNEMD methods, and analyze frequency-dependent phonon contributions via the spectral heat current (SHC) method^[58]. Additional verifications of thermodynamic stability and related discussions for both materials are provided in Figure S12, Figure S13, Figure S14 and Figure S15.

Display Full Size

Figure 4. Phonon vibrational modes of Li₆NII₂ and Li₆NBrBr₂ at the (a) Γ and (b) X points (red arrows indicate the atomic vibrational directions); (c) The projection on the xz plane of two adjacent Li-ion stochastic movements in the conventional cell at 300 K. Blue: Li-ion hopping trajectories; Yellow: Li-ion atomic vibration trajectories; green: Li-ion equilibrium position).

3. Results and Discussion

3.1 Room-temperature thermal properties using molecular dynamics

Figure 5a and Figure 5b present the EMD results at 300 K, obtained using orthogonal supercells sized 5.8 nm × 5.8 nm × 5.8 nm and 5.4 nm × 5.4 nm × 5.4 nm for Li₆NII₂ and Li₆NBrBr₂, respectively (convergence tests with respect to supercell size are shown in Figure S16 and Figure S17). Figure 5c and Figure 5d show that for both double antiperovskites, a supercell with a characteristic lattice length of approximately 5 nm is required to converge the thermal conductivity based on EMD at 300 K. Forty-two independent simulations were conducted for each material. Each simulation began with a 300 ps equilibration stage, followed by a 10 ns production stage performed separately under the NVT and microcanonical (NVE) ensembles. The mean thermal conductivity values along the x, y, and z directions converge over the correlation time. The final reported thermal conductivities are averaged over these three directions throughout the entire correlation period. The predicted thermal conductivities of Li₆NII₂ and Li₆NBrBr₂ at 300 K are 0.367 W/mK and 0.536 W/mK, respectively. Results and convergence tests for EMD simulations at other temperatures are presented in Figures S18, Figure S19, Figure S20 and Figure S21.

Display Full Size

Figure 5. Running lattice thermal conductivity as a function of correlation time using the EMD method for (a) Li₆NII₂ and (b) Li₆NBrBr₂ at 300 K along the x, y and z directions, as well as their averages; Lattice thermal conductivity κ as a function of the characteristic length for (c) Li₆NII₂ and (d) Li₆NBrBr₂ at 300 K from the EMD simulations. EMD: equilibrium molecular dynamics.

Figure 6a and Figure 6b show the HNEMD results at 300 K, obtained using orthogonal supercells of the same sizes as above for 5.8 nm × 5.8 nm × 5.8 nm and 5.4 nm × 5.4 nm × 5.4 nm for Li₆NII₂ and Li₆NBrBr₂, respectively (convergence tests with respect to supercell size are included in the Supplementary Materials).Figure 6c and Figure 6d indicate that a supercell with a characteristic length of about 5 nm is also required to converge the thermal conductivity from HNEMD at 300 K. Given the isotropic nature of both materials, lattice thermal conductivity was computed along the y-direction using HNEMD. The simulations involved a 300 ps equilibration stage followed by a 6 ns production stage under the NVT ensemble at 300 K. The driving force magnitude was optimized around 2.0×10^-4 Å^-1 for both materials to maintain the system within the linear response regime and ensure a sufficiently large signal-to-noise ratio. Thermal conductivities were averaged over the last 2 ns of the production stage, based on eight independent simulations, resulting in values of 0.35 W/mK and 0.531 W/mK for Li₆NII₂ and Li₆NBrBr₂, respectively. These ultralow thermal conductivities are comparable to those of metal oxides (Bi₃Ir₃O₁₁: 1 W/mK)^[59], nitride perovskite (LaWN₃: 2.98 W/mK)^[60], double perovskites (Cs₂AgBiBr₆: 0.33 W/mK^[16], Cs₂PTi₆: 0.15 W/mK^[61], and Cs₂SnI₆: 0.29 W/mK^[62]) and antiperovskites (K₃BrO: 0.73 W/mK^[23], K₃ITe: 0.18 W/mK^[63], and K₃ISe: 0.6 W/mK^[63]) at room temperature. This comparison highlights the strong lattice anharmonicity of both nitrohalide double antiperovskites. Further details of HNEMD simulations are provided in Figure S22, Figure S23, Figure S24 and Figure S25.

Display Full Size

Figure 6. Running lattice thermal conductivity as a function of simulation time using the HNEMD method for (a) Li₆NII₂ and (b) Li₆NBrBr₂ at 300 K along the y direction, as well as their average; Lattice thermal conductivity κ as a function of the characteristic length for (c) Li₆NII₂ and (d) Li₆NBrBr₂ at 300 K from the HNEMD simulations; Spectral lattice thermal conductivity as a function of phonon frequency for (e) Li₆NII₂ and (f) Li₆NBrBr₂ at 300 K. HNEMD: homogeneous nonequilibrium molecular dynamics.

Additionally, we calculated the spectral lattice thermal conductivity of both Li₆NII₂ and Li₆NBrBr₂ using the HNEMD-based SHC method, as shown in Figure 6e and Figure 6f. The dominant phonon contributions to thermal conductivity arise from phonons with frequencies below 5 THz and around 11 THz.

3.2 Temperature-dependent thermal transport properties

Figure 7a and Figure 7b show the temperature-dependent lattice thermal conductivities of Li₆NII₂ and Li₆NBrBr₂, obtained using both EMD and HNEMD methods. The results from the two methods exhibit good agreement. Compared with other antiperovskites such as Li₃OCl and Li₃OBr^[18], as well as Li₃HS, Li₃HSe and Li₃HTe^[64], both Li₆NII₂ and Li₆NBrBr₂ exhibit significantly lower lattice thermal conductivities and reduced sensitivity to temperature over the studied range. For Li₆NII₂, the thermal conductivity remains below 0.4 W/mK from 200 K to 500 K, with a total variation of only 16.44% across this temperature range. Li₆NBrBr₂ displays slightly higher values, decreasing from 0.64 W/mK (0.59 W/mK) to 0.51 W/mK (0.54 W/mK) as temperature increases from 100 K to 400 K, based on EMD (HNEMD) results. The temperature dependence of thermal conductivity for Li₆NII₂ and Li₆NBrBr₂ follows a T^-0.13 and T^-0.11 trend. These temperature-insensitive and ultralow thermal conductivities indicate the glass-like thermal transport behavior of both materials. Figure 7c and Figure 7d display the spectral thermal conductivities of both compounds at different temperatures. Phonons with frequencies below 5 THz are the main contributors to thermal transport, which aligns with the observed flattening of phonon branches above 3 THz and the near-zero phonon group velocities, as shown in the temperature-dependent phonon dispersions and group velocity plots in Figure S26 and Figure S27. An additional substantial contribution comes from phonons near 11 THz. As revealed in Figure 3, Figure S10 and Figure S28, the vibrations of N, I, Br, and some Li atoms correspond to well-defined phonon modes, while the stochastic motion of a subset of Li atoms gives rise to imaginary phonons. Because no distinct phonon branches exist near 11 THz, the thermal conductivity in this region can be attributed to the stochastic motion of Li atoms and their interactions with surrounding atoms. The defining lattice feature of both Li₆NII₂ and Li₆NBrBr₂ is the stochastic hopping of Li atoms among neighboring lattice sites, confined within a single octahedron. This behavior violates the assumptions of perturbation theory, which requires atoms to vibrate around their equilibrium positions, allowing potential energy to be expanded in terms of small displacements. The hopping distances of Li atoms can reach 4.689 Å in Li₆NII₂ and 4.408 Å in Li₆NBrBr₂, which renders perturbation theory inapplicable. In other materials such as (BA)₂PbI₄^[65], Cs₂SnI₆^[66,67], Cs₂PbI₂Cl₂^[68], LaWN₃^[51] and antiperovskites like Li₃OCl^[17] and Na₃FS^[69], strong lattice anharmonicity generally originates from octahedral rotations, tilting, weakly bonded cations (rattler atoms), soft phonon branches, and significant four-phonon scattering. In contrast, the anharmonicity in Li₆NII₂ and Li₆NBrBr₂ is dominated by suppressed acoustic phonons and the presence of imaginary phonons. Mode-projected phonon eigenvector analysis indicates that these imaginary branches span the entire Brillouin zone and result from the stochastic motion of Li atoms. This invalidates the use of conventional transport formalisms such as the phonon Boltzmann transport equation and the Wigner transport equation^[70] in describing phonon transport in these materials. Instead, MD-based simulations and spectral analysis show that phonons below 5 THz and around 11 THz make the largest contributions to thermal conductivity. The significant contribution of phonons above 3 THz, especially those near 11 THz, is mainly associated with the stochastic motion of Li atoms.

Display Full Size

Figure 7. The lattice thermal conductivity of (a) Li₆NII₂ from 200 to 500 K and (b) Li₆NBrBr₂ from 100 to 400 K; Corresponding spectral thermal conductivity of (c) Li₆NII₂ and (d) Li₆NBrBr₂. EMD: equilibrium molecular dynamics; HNEMD: homogeneous nonequilibrium molecular dynamics.

Given their ultralow lattice thermal conductivities across a wide temperature range and their composition from earth-abundant elements, Li₆NII₂ and Li₆NBrBr₂ are promising candidates for low-cost thermoelectric applications, particularly at room temperature. The stochastic motion of Li ions confined within a single conventional unit cell may also be tuned to enhance carrier mobility in thermoelectric devices, thereby improving overall energy conversion efficiency.

4. Conclusions

In summary, we investigated phonon transport in nitrohalide double antiperovskites Li₆NII₂ and Li₆NBrBr₂ using first-principles machine-learning potentials. The confined stochastic motion of Li ions within a single conventional unit cell leads to the emergence of imaginary phonon modes throughout the entire Brillouin zone and causes the breakdown of perturbation theory. Both equilibrium and nonequilibrium molecular dynamics simulations confirm that Li₆NII₂ and Li₆NBrBr₂ exhibit ultralow lattice thermal conductivities with glass-like temperature dependence. Spectral thermal conductivity analysis reveals that phonons with frequencies below 5 THz and around 11 THz dominate the overall thermal transport. Notably, the significant contribution near 11 THz arises from the stochastic, confined movement of Li ions. This study uncovers unconventional atomic dynamics and a distinct phonon transport mechanism in double antiperovskites, offering valuable insights for the design of low-thermal-conductivity materials for thermoelectric applications.

Supplementary materials

The supplementary material for this article is available at: Supplementary materials.

Acknowledgements

Qi Wang acknowledges the support from the start-up funding of the Suzhou Institute for Advanced Research at the University of Science and Technology of China.

Authors contribution

Li Y: Data curation, formal analysis, investigation, methodology and writing-original draft.

Zheng W: Data analysis and interpretation.

Pu JH: Resources and supervision.

Wang Q: Conceptualization, formal analysis, methodology, resources, software, supervision, writing-review and editing.

Jiang JH: Supervision, writing-review and editing.

Conflict of interest

Jian-Hua Jiang is an Editorial Board member of Thermo-X. The other authors declared that there are no conflicts of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The data that support the findings of this study are available within the article and its supplementary materials. Other questions regarding this study can be directed to the corresponding authors.

Funding

This study was supported by the Natural Science Foundation of Shandong Province (Grant No. ZR2023QE223)

Copyright

© The Author(s) 2025.