Molecular dynamics study of thermal conductivities of cubic diamond,lonsdaleite,and nanotwinned diamond via machine-learned potential

Jia-Hao Xiong(熊佳豪), Zi-Jun Qi(戚梓俊), Kang Liang(梁康), Xiang Sun(孙祥),Zhan-Peng Sun(孙展鹏), Qi-Jun Wang(汪启军), Li-Wei Chen(陈黎玮),Gai Wu(吴改),‡, and Wei Shen(沈威),§

1The Institute of Technological Sciences,Wuhan University,Wuhan 430072,China

2School of Power and Mechanical Engineering,Wuhan University,Wuhan 430072,China

3Department of Mechanical Engineering,The University of Tokyo,Tokyo,113-8656,Japan

4Wuhan University Shenzhen Research Institute,Shenzhen 518057,China

5Hongyi Honor College,Wuhan University,Wuhan 430072,China

Keywords: diamond,neuroevolution potential,molecular dynamics,thermal conductivity,phonon transport

1.Intro ducti on

Diamond is a highly valuable material due to its optimal properties such as the wide bandgap,the high breakdown field strength, the remarkable mechanical properties, and the outstanding thermal conductivity in particular.[1-4]Diamond has some polytypes such as cubic diamond, lonsdaleite (hexagonal diamond), and nanotwinned diamond(NTD),which have different lattice structures.Among them, the lonsdaleite and the NTD both have important research values because of their special thermodynamic and mechanical properties.The direct synthesis of NTD typically involves with the use of a precursor of onion carbon nanoparticles at high-pressure and hightemperature,resulting in an average twin thickness of 5 nm.[5]The lonsdaleite can usually be regarded as a special type of NTD with a minimum twin thickness of 0.206 nm.[6]The twin boundaries(TBs)are interfaces that separate regions of a crystal with different crystallographic orientations.In the case of NTD, the presence of TBs can be attributed to the crystal’s formation, where the cubic and hexagonal diamond crystals grow in a partially misaligned fashion,resulting in distinct regions with different orientations.The formation of lonsdaleite through the emergence of crystal dislocations in the cubic diamond can be recognized as the existence of TBs to some extent.The TBs play a crucial role in optimizing mechanical properties in nanotwinned structures like lonsdaleite and NTD.[7]

Previous studies show that NTD has a Vickers hardness as high as 200 GPa, which has unprecedented hardness, and also has excellent thermal stability and an in-air oxidization temperature which exceeds 200 K higher than that of the cubic diamond.[5]Additionally, the mechanical and thermal properties exhibited by NTD are also basically evident in lonsdaleite.[6,8]The lonsdaleite is famous for its extraordinary hardness, exceeding the hardness of NTD, which is expected to exceed 200 GPa.[6]Based on the excellent properties of the NTD and lonsdaleite, such as the mechanical and thermodynamic properties, the NTD and lonsdaleite structures can exhibit promising potential applications in semiconductor field,such as the semiconductor devices in the high-temperature and high-pressure environment.However, the presence of TBs has a negative influence on the thermal conductivity of the NTD and lonsdaleite because of the rising possibility of phonon scattering,which is not conducive to their applications in semiconductor fields.Thus,the thermal conductivity evolutions of the NTD and lonsdaleite at different temperatures are worth further exploring for the application of the mentioned diamond structures in semiconductor fields.

Fig.1.Lattice structures of(a)cubic diamond(111),(b)lonsdaleite,and(c)NTD,studied in this work.

The density functional theory (DFT) and the molecular dynamics (MD) with empirical potentials and machinelearned (ML) potentials can both provide approaches to calculating the thermal conductivities of the diamond structures.[9,10]The DFT method can precisely calculate the thermal conductivity and phonon dispersion but it is only a time-consuming method for complex crystal structures.[11,12]The Tersoff potential (Tersoff-1989[13]) is a commonly-used empirical potential for the MD simulations of semiconductors such as diamond and silicon.Compared with the DFT calculations, the MD calculations based on the Tersoff potential are advantageous in terms of the simulation scale and speed.However, since it is based on empirical potential,there are limitations in accurately describing the energy and atomic forces,leading to inaccuracy in depicting the structural properties.[14,15]

In previous studies,[16,17]the thermal conductivities and phonon transport mechanisms of the three diamond structures have been explored based on other potentials.The thermal conductivities, phonon spectra, and other characteristics that describe the differences among their thermal properties have been well explored based on the molecular dynamics (MD) approach.[18-20]In order to better understand the thermal conductivities of three structures and save computing time, a machine-learned method and ab initio calculation method (from first principles) have been used to obtain a high-precision potential named neuroevolution potential(NEP).[21]The NEP method based on MD method outperforms the DFT approach in terms of the computational speed and simulation scale while providing greater accuracy than the empirical potentials.[22]In this paper, the NEP is trained in a neural network and adopted in calculating thermal properties in cubic diamond,lonsdaleite and NTD.The thermal conductivities of the mentioned three structures are obtained and the phonon dispersions are revealed based on the Tersoff-1989,DFT,and NEP.

2.Method

In this paper,the NEP of cubic diamond,lonsdaleite,and NTD are obtained by combining machine-learned potentials(MLPs)with the density functional theory(DFT).The datasets in the NEP training process are obtained through the DFT calculations.The thermal conductivities of the three diamond structures are calculated based on the equilibrium molecular dynamics(EMD).

2.1.Density functional theory

The DFT from ab initio calculations has been carried out by using the VASP method.[23]For the exchangecorrelation functional,the generalized gradient approximation with Perdew-Burke-Ernzerhof parametrization is applied.[24]The plane wave cutoff is set to be 520 eV.The optimized lattice parameters of the unit cell of diamond, lonsdaleite, and NTD are 2.53 ˚A×4.38 ˚A×6.19 ˚A, 2.51 ˚A×4.35 ˚A×4.18 ˚A,and 2.52 ˚A×4.37 ˚A×12.44 ˚A with 12, 8, and 24 atoms, respectively.During the cell optimization,the convergence criteria are 10-5eV for the energy difference in electronic selfconsistent calculation and 0.01 eV/˚A for the residual forces on atoms.Thek-point grids for the unit cell of diamond, lonsdaleite,and NTD are employed to be 11×7×5,11×7×7,and 11×7×2, respectively.The structure of the diamond, lonsdaleite,and NTD in the DFT simulation are 10×6×6,4×3×3,and 4×3×1 supercells,respectively.Thek-point grids for the supercells of diamond,lonsdaleite,and NTD are employed as the Gamma point.The three structures are first equilibrated for 0.2 ps in time steps of 0.5 fs, and then the snapshots are collected each step during production runs for 0.3 ps at 800 K.Each structure has 500 snapshots generated in the training process, and the total number of snapshots used for the training process is 1500.Each structure has 100 snapshots generated in the testing process,and the total number of snapshots used for the testing process is 300.After the calculations in DFT,we have had a training set and testing set,which are to be applied to the training of NEP.

2.2.Neuroevolution potential(NEP)

In graphics processing units molecular dynamics(GPUMD),the NEP3 has been adopted as the method of calculating the potential energy function.[21,25]Various molecular descriptors need setting and defining for the acquisition of the NEP.The NEP uses a simple feedforward neural network to represent the energy and force of atoms.The DFT method has been used to obtain precise atomic force and energy data.[24]All generated data in the DFT method, which are divided into the training set and the testing set, involve with nearly 3×104configurations and 9.5×105atoms.The training set has 8.43×105atoms and the testing set has 1.164×105atoms, which account for 88% and 12% of the total data,respectively.After obtaining the training set and the testing set, it is meant that the training process is completed in the GPUMD package.[26,27]

2.3.Equilibrium molecular dynamics(EMD)

The EMD is mainly based on the Green-Kubo method.The Green-Kubo method relies on the fact that the heat flow in a particle system in equilibrium fluctuates around zero and heat flux vectors and their dependencies are obtained by simulation.[28]The thermal conductivity is then predicted via the Green-Kubo relationship based on the time which takes for the thermal autocorrelation function to decay to zero.

In the Green-Kubo method,[29,30]the running thermal conductivity along thexdirection can be expressed as an integral of the heat current autocorrelation function(similar expressions are applied to the other directions):

wherekBis the Boltzmann’s constant,Vis the volume of the cubic diamond system or lonsdaleite systemor NTD system,Tis the absolute temperature,andtis the correlation time,Jx(0)andJx(t) are the total heat current of the system at two time points separated by an interval oft.The heat current at a given time depends on the location and velocities of the particles in the system.

The thermal conductivities in cubic diamond,lonsdaleite,and NTD are calculated by using the EMD method.The three directions of the three crystal lattices are established as [112] (X), [110] (Y), and [111] (Z), respectively, which correspond to three structures in Fig.1.For the EMD calculation, a 30.32 ˚A×30.63 ˚A×30.94 ˚A supercell consisting of 6408 atoms for the cubic diamond unit cell, a 30.32 ˚A×30.63 ˚A×28.88 ˚A supercell containing 4704 atoms for the lonsdaleite unit cell,and a 30.32 ˚A×30.63 ˚A×37.13 ˚A supercell composed of 5040 atoms for the NTD unit cell are utilized in the simulation.

Fig.2.Plots of the thermal conductivity versus atom number per cell for the cubic diamond (red line), lonsdaleite (orange line), and NTD(blue line).

It can be indicated from previous study that such a scale is large enough to neglect the finite size in the Green-Kubo method.[31]Additionally, we also calculate the thermal conductivities on different atomic scales at 300 K, such as the cubic diamond with 5080 atoms and 7300 atoms, the lonsdaleite with 3696 and 5824 atoms, and the NTD with 4034 and 6048 atoms, to exclude the size effect.It can be found that the differences in thermal conductivity among cubic diamond structure, lonsdaleite structure, and NTD structure for different cell sizes are within 5% as shown in Fig.2.Thus,it can be further conducted that the size effect on the thermal conductivity calculation for each of cubic diamond structure,lonsdaleite structure,and NTD structure in this work could be ignored.

The thermal conductivities of the three structures all have been calculated in thezdirection at 300 K, 400 K, 500 K,600 K, and 700 K, respectively.In the thermal conductivity calculations,the diamond structure should first reach the preset temperature in the NPT ensemble, and then the heat conduction is simulated in the NVT ensemble at a maximum correlation time oftmax=0.5 ps.The elastic constant or inverse compressibility parameter needed in the barostat is set to be 1050 GPa.In the NVT ensemble, we have a total of 8×106steps and output the thermodynamic information every 8×104steps.

3.Results and discussion

3.1.NEP model training and validating

After training the NEP,we calculate the root mean square error (RMSE) of the energy and atomic force by comparing the data in the testing set and the predictions from NEP.The comparison of the energy and atomic force between the NEP and DFT calculations are shown in Fig.3.In Fig.3(a), the RMSE of energy per atom is 0.24 meV/atom.In Fig.3(b),the RMSE of atomic force is 51.8 meV/˚A in the NEP model.The MLPs are recognized to have excellent performance when the RMSE of energy and atomic force falls in the range of several meV/atom and several hundred meV/˚A, respectively.[32,33]The small RMSE value obtained for energy and atomic force can illustrate the high accuracy of the NEP model,comparable to the accuracy of the DFT calculations.The linear approximation of the DFT data and NEP predictions also suggest that the potential has been satisfactorily trained.It has been shown that such a high accuracy can ensure high-precision predictions of the physical properties of materials.[34]

Fig.3.(a) NEP energy versus DFT energy and (b) NEP force versus DFT force, with comparison of between the training set and testing set based on the DFT calculated data.The RMSE is obtained by comparing the data between NEP and DFT calculations in training process,and yellow solid line represents ftiting curve for the training set.(c)Loss function of the training process.

Figure 3(c)shows the evolutions of the various loss functions with generation.With the increase of generation, the RMSE of energy and force are reduced and converged, although with some oscillations in the beginning.In contrast,the loss functions for L1 and L2 regularizations exhibit a trend of first increase and subsequent decrease, reflecting the influence of regularization on magnitude of weight and bias parameters in the neural network and highlighting its effectiveness.Without the regularization (settingλ1andλ2to zero or very small values),the weight and bias parameter in the neural network would increase,which can easily lead to overfitting.The final hyper-parameters obtained from the training process are listed in Table 1.

Table 1.Hyper-parameters in NEP model training process.

In Table 1,rRc(rAc)is the cutoff radius for the radial(angular)component of the descriptor,nRmax(nAmax)is Chebyshev polynomial expansion order for the radial (angular) components,lmaxis the Legendre polynomial expansion order for the angular components, neuron is the number of neurons in the hidden layer of the neural network,NRbas(NAbas)number of radial and angular basis functions, and generation is the maximum number of generations to be evolved.

3.2.Thermal conductivity of cubic diamond, lonsdaleite,and NTD

In the EMD method, the thermal conductivities (from 300 K to 700 K)are calculated as an integral of the heat current autocorrelation function.Owing to the well-known persistent noise in the autocorrelation function, 50 independent integrals often drift away in only 50 ps, and finally converge to obviously different values at 200 ps, which poses a great challenge to accurately deduce the thermal conductivity from such calculations.[35]As shown in Fig.4, the optimized thermal conductivity value could be achieved by carrying out 50 independent simulations and averaging the thermal conductivities from the simulation results.Based on the above method,the thermal conductivities of the cubic diamond, lonsdaleite,and NTD at 300 K are obtained.The differences in thermal conductivity among the three diamond structures are partially due to the existence of the TBs in NTD and lonsdaleite,which can increase the thermal resistance of material.It can be seen from Fig.5 that the thermal conductivities of the three diamond structures decrease in a similar trend with the increase of temperature.As the temperature increases from 300 K to 700 K,the phonon scattering will be enhanced,especially for the Umklapp(U)scattering,which will lead to the decrease in the phonon relaxation time, thus resulting in the decrease of thermal conductivity.[36]

Fig.4.Correlation-time-dependent thermal conductivities of(a)cubic diamond,(b)lonsdaleite,and(c)NTD at 300 K,obtained by EMD method,with blue dashed line denoting the result of each independent simulation, red dashed line representing the average curve of 50 independent simulations, and orange area referring to the 95%confidence interval band.

Fig.5.(a)Variations of thermal conductivities with temperature for cubic diamond,lonsdaleite,and NTD based on NEP,with green solid line representing experimental data of cubic diamond cited from Ref.[37],and dashed line denoting DFT results of cubic diamond cited from Ref.[38].(b)Thermal conductivities of lonsdaleite based on DFT(from Ref.[1]),NEP,and Tersoff-1989.

In the DFT method, the calculation of thermal conductivity frist involves computing the interatomic force constants,and then solving the Boltzmann transport equation(BTE).As shown in Fig.5, the full iterated solution in BTE (BTE Full)mainly considers the normal (N) scattering, so the calculated thermal conductivity is higher than the experimental value.The relaxation transport equation in BTE (BTE RTA) independently considers the Umklapp (U) scattering and normal(N)scattering,and the calculated thermal conductivity is lower than the experimental values.The thermal conductivity of the cubic diamond at 300 K based on the NEP and Tersoff-1989 are respectively 2507.3 W·m-1·K-1and 1508 W·m-1·K-1,which are about 18% and 51% lower than the DFT data(3086 W·m-1·K-1, BTE Full), respectively.The thermal conductivity of the lonsdaleite at 300 K based on the NEP and Tersoff-1989 are respectively 1557.2 W·m-1·K-1and 1178 W·m-1·K-1, which are 13% and 34% lower than the DFT result (1790 W·m-1·K-1).Anharmonic force constants have not been calculated by using DFT , because the calculation cost for NTD is too high.Therefore, it is difficult to obtain reliable DFT results for NTD in current researches.By using the NEP, the thermal conductivity of NTD at 300 K is obtained to be 985.6 W·m-1·K-1, which is 24% higher than that obtained by using Tersoff-1989 (794 W·m-1·K-1).The thermal conductivity of NTD at 300 K based on the NEP(985.6 W·m-1·K-1) is lower than that of the cubic diamond based on the NEP(2507.3 W·m-1·K-1).

Previous study[16]has shown that the thermal conductivity of NTD, evaluated via using an advanced potential with EMD,is lower than that of cubic diamond evaluated by using the same advanced potential with EMD, which is still higher than that of NTD evaluated by using an empirical potential with EMD at 300 K.Therefore, the thermal conductivity of NTD based on the NEP proposed in this paper is credible to some extent.Based on the NEP, the thermal conductivity of the cubic diamond and lonsdaleite derived from the EMD method are closer to the calculation results from DFT than the results based on the Tersoff potential.Thus,the MD simulations by utilizing the NEP can exhibit satisfactory accuracy and remarkable computational efficiency.The computational efficiency for the MD simulations based on the NEP can exceed that of the DFT method by nearly two or three orders of magnitude,while the thermal conductivity calculation based on the NEP will take more time than that based on the empirical Tersoff potential, usually differing by an order of magnitude.[22,39]

4.Phonon transport in cubic diamond, lonsdaleite,and NTD

When the temperature is higher than 300 K, the phonon scattering plays a more important role in thermal conductivity than the electron scattering.[40]Thus, the influence of the phonon dispersion on the thermal conductivity in the cubic diamond, lonsdaleite, and NTD are discussed in this subsection.The differences in thermal conductivity among the three diamond structures can be clarified by analyzing the phonon dispersions and comparing the previous studies.The NEP and Tersoff-1989 potential are used as potentials to calculate the phonon dispersions of the three diamond structures for the comparison with the DFT results.As shown in Fig.6,the phonon dispersions calculated by the NEP are more accurate and closer to the DFT results than those obtained by the Tersoff-1989 in three structures of diamond.The phonon dispersions predicted by the Tersoff-1989 (dashed line) deviate from the DFT calculations (dotted line) to a greater extent.The frequencies of the optical branches predicted by the Tersoff-1989 are overestimated by nearly 20% in comparison with the DFT results in the above-mentioned three diamond structures, thus affecting the phonon-phonon scattering processes, which will have some influence on thermal conductivity.[33]In Table 2, the transverse frequencies at high-symmetry points (Γ,K,X,L) in cubic diamond derived from the NEP are 39.32 THz, 29.50 THz, 22.96 THz,and 16.81 THz,which better match the DFT and experimental results than those from the Tersoff-1989.The phonon dispersion of cubic diamond, lonsdaleite and NTD based on NEP are in better agreement with the DFT calculation results,providing further evidence for the accuracy of the NEP than for the Tersoff potential.Thus, the NEP is capable of accurately characterizing the properties of the three diamond structures in comparison with the Tersoff-1989 potential, which in turn bolsters the validity of the derived thermal conductivities for cubic diamond,lonsdaleite and NTD.

Fig.6.(a) Phonon dispersion in cubic diamond along the high-symmetry direction based on DFT, Tersoff potential, and NEP, (b) phonon dispersion in lonsdaleite and(c)phonon dispersion in NTD.(d)Phonon participation ratio of cubic diamond,lonsdaleite,and NTD.

The comparison of the phonon dispersion among the three diamond structures in this subsection can reveal that the lattice vibration modes of the NTD and lonsdaleite structures are more complex than those of the cubic diamond.The optical phonon branches generated by Brillouin zone folding in lonsdaleite and NTD fall to the acoustic phonon range, resulting in no bandgap between the acoustic phonon branch and optical phonon branch in the lonsdaleite and NTD.This phenomenon highlights the limitations of traditional empirical potentials, which typically have a limited number of molecular descriptor symbols and may inaccurately describe phonon dispersions.[15]It can be found that the NEP is more accurate than the Tersoff-1989 in the phonon dispersion calculation for all the phonon branches, as validated by a high degree of consistency between the DFT result and NEP result, and this level of accuracy is proven to be essential in describing the energy and momentum conservation in phonon-phonon scattering processes for predicting the thermal conductivity.

The small mass of carbon atom and the strong covalent bond in the diamond structure can result in the atomic vibrations primarily near the minimum energy state,allowing a simple harmonic approximation of the phonon transport dominated by the normal (N) scattering, which leads to the low intergranular thermal resistance in diamond structures.[38]For any specific material, the thermal conductivity primarily depends on the phonon relaxation time and phonon velocity.It can be shown in previous researches that phonon velocity is regarded as the primary factor affecting thermal conductivity,[41]and the TBs will increase the possibility of Umklapp(U)scattering,which can cause the phonon energy and group velocity to decrease and lead the thermal conductivity to lessen.[42,43]It is shown in previous research[16]that the phonon velocity can decrease with the increase of temperature.It can be seen that the average phonon velocity in the lonsdaleite and NTD are lower than that in the cubic diamond due to the existence of TBs in NTD and lonsdaleite,eventually resulting in the lower thermal conductivity than the thermal conductivity of the cubic diamond.[16,17]Additionally,we calculated the phonon participation ratio for each of cubic diamond,lonsdaleite,and NTD.The participation ratio for lonsdaleite and NTD are lower than that of the cubic diamond in the low-frequency range and midfrequency range,respectively,as shown in Fig.6(d).Thus,the phonon in lonsdaleite structure and NTD structure become localized and are not effectively involved in the phonon transport process in a specific frequency range.This localization phenomenon leads to lower thermal conductivity in lonsdaleite and NTD than that in cubic diamond.

Table 2. Phonon frequencies (in units of THz) at high-symmetry points of Brillouin zone, calculated based on DFT, Tersoff-1989, and NEP along with experimental (exp.) data[44] for comparing with the results of cubic diamond.

5.Conclusions

In this work,it is demonstrated that the machine-learned neuroevolution potential (NEP) trained in GPUMD can describe the interatomic interactions of the cubic diamond,lonsdaleite,and nanotwinned diamond(NTD)with a much higher precision than the empirical potentials.By combining the NEP method with the EMD method, the phonon dispersions and thermal conductivities of the three diamond structures are well predicted,showing great agreement with the DFT results.Our approach can provide a more efficient alternative to the DFT calculations and a more accurate alternative to the widely-used Tersoff-1989 potential.

It can be revealed by analyzing the phonon spectra that the lower thermal conductivity of the lonsdaleite and NTD than that of the cubic diamond can be attributed to the lower phonon group velocity.Our findings can highlight the potential of machine learning methods to enhance the understanding of the thermodynamic properties of the three diamond structures and provide guidance for designing the diamond devices with tailored properties.

Acknowledgements

Project supported by the National Natural Science Foundation of China (Grant Nos.62004141 and 52202045), the Fundamental Research Funds for the Central Universities,China (Grant Nos.2042022kf1028 and 2042023kf0112), the Knowledge Innovation Program of Wuhan-Shuguang, China(Grant Nos.2023010201020243 and 2023010201020255),the Natural Science Foundation of Hubei Province,China(Grant No.2022CFB606),and the Guangdong Basic and Applied Basic Research Fund: Guangdong-Shenzhen Joint Fund, China(Grant No.2020B1515120005).