MATERIALS TRANSACTIONS
Online ISSN : 1347-5320
Print ISSN : 1345-9678
ISSN-L : 1345-9678
Materials Chemistry
High-Precision Prediction of Thermal Conductivity of Metals by Molecular Dynamics Simulation in Combination with Machine Learning Approach
Qi KongYasushi Shibuta
Author information
JOURNAL FREE ACCESS FULL-TEXT HTML

2023 Volume 64 Issue 6 Pages 1241-1249

Details
Abstract

Molecular dynamics (MD) simulation is a powerful tool to estimate materials properties from atomistic viewpoint. However, the scope of application of MD simulations is limited to problems where the Newton’s equation of motion for atoms is dominant. Therefore, it is inherently insufficient to estimate thermal conductivity of metallic materials, which consists of phononic and electronic components. In this study, machine learning (ML)-based regression model is employed to predict thermal conductivity of metals with high accuracy using deficient results from MD simulations. A regression analysis with the least absolute shrinkage and selection operator (Lasso) including electrical conductivity as predictor variables successfully predict the thermal conductivity of metals with negative temperature dependence, which indicates a significant contribution of electrons to thermal conduction in metals. It should be stressed that our prediction is better than the conventional estimation from the Wiedemann–Franz law. This study shows us new possibilities of new ML approach for improving the accuracy of physical properties obtained from MD simulations.

1. Introduction

Molecular dynamics (MD) is a powerful computer simulation method to study properties of materials. Due to recent improvements in high-performance computing environment1) combined with various methods for accelerating simulations such as metadynamics,2) hyperdynamics3) and coarse-grained modeling,4) efficiency of MD simulations has been increased significantly and spatio-temporal scale that MD simulation can handles has expanded greatly, making it possible for atomistic simulations close to realistic environment. For example, multi-graphics processing unit (GPU) parallel computation on high-performance supercomputer enabled billion-atom MD simulations to capture the nature of homogeneous nucleation during solidification of pure metal.5)

However, the scope of application of MD simulations is limited to problems where the classical mechanics of atomistic scale is dominant68) because MD numerically solves classical Newton’s equations of motion to describe behavior trajectory of atoms. For example, MD simulation has been widely used to study thermal conductivity, which is an important property for designing and processing of materials. Although thermal conductivity of various materials has been successfully estimated by MD simulations,914) there are yet few successful reports on metallic materials. It is because that MD not accounting for the behaviors of electrons, is only able to calculate lattice thermal conductivity attributing to phonons,9) although electrons contribute greatly to thermal conduction of metals.15) In order to consider the behavior of electrons, it is essential to perform first-principles calculations,16) which requires much higher computation cost than MD simulations.

Recently, machine learning (ML) approaches have become popular in the field of materials science.1721) In general, the basic idea of ML applied in materials science is to predict properties of materials quantitatively from existing databases8,19) while accepting the decrease of physical insights in most cases. Because MD simulation is an efficient way to generate preferable number of physically meaningful data of desirable properties, MD simulation combined with ML approach has become a powerful method to predict properties of materials.1921) Another application of ML approach combined with MD simulation is to construct interatomic potentials called ML interatomic potentials.2224) This method has contributed to improve the accuracy of MD simulation from energetic viewpoint. However, it is not yet straightforward to tackle properties such as thermal conductivity as far as solving the equation of motion classically even though the accurate interatomic potential is employed as described above. Although there exists discrepancy between results of MD simulations and the true materials properties, it should be noted that there must be some unknown or hidden physical relationship between them. Therefore, it is expected that ML approach is useful to reveal such hidden physical relationship.19,21) In this study, we proposed new ML approach which combines the result of MD simulation and existing data to predict thermal conductivity of metallic materials with high accuracy. We firstly estimate thermal conductivity of various metals systematically by the conventional MD simulation and compare simulation results with reference experimental data. Then, we apply a ML-based regression model to predict thermal conductivity from MD results and relevant predictor variables. Our results indicate a feasible strategy to find a hidden physical relationship between MD simulation and experimental data and to improve the accuracy of physical properties obtained from MD simulation.

2. Molecular Dynamics Simulation

2.1 Simulation methodologies

Equilibrium molecular dynamics (EMD) simulation is one of the most commonly used atomistic simulation methods to study transport properties of materials due to its relatively high computational efficiency.25) In this study, thermal conductivities of various representative metals with face-centered cubic (fcc) crystal structure (Al, Ni, Cu, Pd, Ag, Pt, Au) and body-centered cubic (bcc) structure (V, Cr, Fe, W) are calculated through EMD simulation at various temperatures. The embedded atom method (EAM) potential26,27) and the Lennard-Jones (LJ) potential with parameters shown in Table 128) are employed as interatomic potentials of those fcc metals. EAM potentials fitted by Mendelev et al.29) and Han et al.,30) and the modified embedded atom method (MEAM) potential fitted by Lee et al.31) are employed as interatomic potentials of those bcc metals. Only MEAM potential is employed for the Cr case. All simulations are performed using Large-scale Atomic/Molecular Massive Parallel Simulator (LAMMPS).32)

Table 1 Parameters of the Lennard-Jones potential.28)

In the EMD simulation, thermal conductivity tensor κuv of materials can be obtained based on Green-Kubo formula33,34) as follows   

\begin{equation} \kappa_{uv} = \frac{V}{k_{B}T^{2}} \int_{0}^{\infty}C_{uv}(t)dt \end{equation} (1)
where V is the volume of simulation system, kB is the Boltzmann constant, T is the absolute temperature of the system and Cuv(t) is the heat flux autocorrelation function (HFACF), which is defined as   
\begin{equation} C_{uv}(t) = \langle J_{u}(0)J_{v}(t)\rangle \end{equation} (2)
where Ji refers to the i-component of heat flux. For an isotropic 3D material, the thermal conductivity κ can be considered as the average of the three components of thermal conductivity tensor as follows   
\begin{equation} \kappa = \frac{(\kappa_{xx} + \kappa_{yy} + \kappa_{zz})}{3}. \end{equation} (3)

In this study, each simulation system consists of 10 × 10 × 10 unit cells (4000 atoms of fcc metals or 2000 atoms of bcc metals). The system is simulated firstly in an NVT-constant ensemble for 20 ps to obtain an equilibrium structure at the specified temperature and then in an NVE-constant ensemble for calculation of the heat flux and the HFACF. The velocity-Verlet method is used to integrate the classical equation of motion with a time step of 10.0 fs for fcc metals and 1.0 fs for bcc metals. The Nose-Hoover thermostat35,36) is employed to control temperature. It is impracticable to calculate the HFACF of the infinite time, but we found that the average of running thermal conductivity calculated from the HFACF according to Green-Kubo formula could become convergent when the integration time is long enough. Figure 1 shows a calculation example of the normalized HFACF (NHFACF) and the corresponding running thermal conductivity of Cu at 300 K under EAM potentials, where the autocorrelation time is from 0 to 20 ps. It is found that the average of running thermal conductivity converges to 12.6 W/(mK).

Fig. 1

(a) Normalized heat flux autocorrelation function (NHFACF) of Cu calculated by equilibrium molecular dynamics simulation at 300 K with the EAM potential, where the autocorrelation time is from 0 to 20 ps. Gray and blue lines represent results of ten simulations and their averages. (b) Running thermal conductivity calculated using NHFACF above based on Green-Kubo formula. Gray and red lines represent the results of ten simulations and their averages.

2.2 Thermal conductivity obtained from molecular dynamics simulation

The above method is employed to calculate the thermal conductivities between 100 K and the melting point at 100 K intervals for fcc metals, and between 300 K and the melting point of bcc metals at 100 K intervals, respectively. Calculations are performed ten times for each condition to gather statistics. Figures 2 and 3 show the thermal conductivity of fcc and bcc metals at various temperatures calculated by MD simulations with different interatomic potentials, respectively. It can be observed that thermal conductivities calculated by MD simulation tend to decrease with increasing temperature for all cases. Such negative temperature dependence is the characteristic of lattice thermal conductivity, which can be simply expressed as eq. (4) derived from kinetic theory of phonon37)   

\begin{equation} \kappa = \frac{1}{3}C_{V}vl \end{equation} (4)
where CV is the specific heat at constant volume, v is the sound velocity and l is the mean free path of phonons. CV and v are approximately constant above Debye temperature while l is proportional to 1/T,38) and therefore lattice thermal conductivity decreases with increasing temperature.

Fig. 2

Thermal conductivity of fcc metals at various temperatures calculated by molecular dynamics simulation using embedded atom method (EAM) and Lennard-Jones (LJ) potentials. Corresponding experimental values (Exp)15) are also plotted.

Fig. 3

Thermal conductivity of bcc metals at various temperatures calculated by molecular dynamics simulation using embedded atom method (EAM) and modified EAM (MEAM) potentials. Corresponding experimental values (Exp)15) are also plotted.

Furthermore, the calculation results are compared with the experimental data.15) As expected, simulation results for all metals deviate significantly from the experimental values and obvious gaps between them can be observed except for few cases of low temperatures. In metallic materials, principal carriers of heat are electrons in addition to phonons and therefore the total thermal conductivity of metals should consist of two parts: phononic and electronic components. However, only phononic component can be obtained from the MD simulation because MD makes use of the classical Newtonian equations of motion to describe the trajectory of atoms, thus information of electronic component is not considered explicitly in the simulation. Note that there is not a problem of the accuracy of interatomic potentials but an inherent problem of scope of application of the MD simulation as described above. Therefore, it is not surprising to find the discrepancy between results of MD simulation and experimental values. In the next section, an ML approach is proposed for constructing the unknown relationship between them and for predicting accurate thermal conductivity.

3. Machine Learning Approach

3.1 Least absolute shrinkage and selection operator (Lasso)

A regression analysis, which is one of the most widely used ML algorithms, is usually employed to estimate the relationship between a set of variables (Y1, Y2, …, Yi) and another set of variables (X1, X2, …, Xi) and predict the value Y based on the established relationship when the value of X is provided. Here, Xi is called independent variables, also known as predictor variables in regression analyses, and Yi is called dependent or response variables. Among many regression analysis methods, the least absolute shrinkage and selection operator (Lasso) is a useful one which can enhance the prediction accuracy and the interpretability of the statistical model.39) The characteristic of lasso is the introduction of L1 regularization added to the general residual sum of squares when calculating the loss function, aiming to overcome overfitting of linear regression. When estimating a linear model using Lasso, one should minimize the following loss function   

\begin{equation} \sum\nolimits_{i = 1}^{n}\left(y_{i} - \sum\nolimits_{j}x_{ij}\beta_{j}\right)^{2} {}+ \lambda \sum\nolimits_{j = 1}^{p}|\beta_{j}| \end{equation} (5)
where β is the linear parameter of linear function xijβj, and λ is the tuning parameter which determines the effect of L1 regularization.

3.2 Regression analyses for the prediction of thermal conductivity

Two Lasso regression models with different sets of predictor variables are employed for thermal conductivity prediction of metals (hereafter, referred to as ML models). Figure 4 shows thermal conductivities of fcc metals Al, Cu, Ag, Au at various temperatures predicted by two ML models respectively versus corresponding experimental values in scatter plots. For both ML models (a) and (b), thermal conductivity (and its square) obtained from MD simulations with the EAM potential, and temperatures are employed as the predictor variables. In addition, electrical conductivities at room temperature40) are also included in the predictor variables of model (b). The dependent variables of model (a) and (b) can be expressed as following linear regression equations, respectively:   

\begin{equation} (\text{a})\ y = w_{0} + w_{1}\kappa_{\textit{MD}} + w_{2}\kappa_{\textit{MD}}{}^{2} + w_{3}T \end{equation} (6)
  
\begin{equation} (\text{b})\ y = w_{0} + w_{1}\kappa_{\textit{MD}} + w_{2}\kappa_{\textit{MD}}{}^{2} + w_{3}T + w_{4}\sigma \end{equation} (7)
where κMD, T and σ refers to thermal conductivity obtained from molecular dynamics simulation, temperature and electrical conductivity, respectively.

Fig. 4

Comparison of predicted and experimental values for thermal conductivities of fcc metals, Al, Cu, Ag and Au at various temperatures predicted by the Lasso regression analysis with predictor variables of (a) (κMD, κMD2, T) and (b) (κMD, κMD2, T, σ), where κMD, T and σ refers to thermal conductivity obtained from molecular dynamics simulation, temperature and electrical conductivity, respectively.

From Fig. 4, good agreement between experimental values and thermal conductivities predicted by the ML model including electrical conductivities is observed while thermal conductivities predicted by ML model without electrical conductivities deviate from experimental values obviously. As mentioned above, MD simulation is inherently lacking the information of electronic thermal conductivity. Thus, a good prediction of thermal conductivity is obtained by the ML model including electrical conductivities because it compensates the shortage of information in MD simulation successfully.

The weights of each component in predictor variables of the model represent their contributions to the prediction result. Table 2 shows weights of each component in predictor variables for the model including electrical conductivities. The greatest weight of electrical conductivity is observed, which indicates again the significant contribution of electrons to the thermal conductivity of metallic materials. This conclusion is in accordance with the general principle of thermal conductivity of metals15) and thus can be regarded as the physical interpretation of this ML model. In general, thermal conductivity κ can be approximately estimated from electrical conductivity and temperature according to Wiedemann–Franz law41) as following:   

\begin{equation} \kappa = L\sigma T \end{equation} (8)
where L is the Lorenz number as L = 2.44 × 10−8 WΩK−2. Thermal conductivity predicted by machine learning is compared with Wiedemann–Franz law in Fig. 5. It is found that machine learning results is more consistent with experimental results than Wiedemann-Franz law when temperature is relatively low. The weight of electrical conductivity as a predictor variable in the ML models is actually much higher than thermal conductivity of MD, but at the same time thermal conductivity of MD also becomes larger at lower temperature. Therefore, we infer that MD results are still important for more accurate prediction of the thermal conductivity of metallic materials although electrical conductivity is included.

Table 2 Weights of standardized components in predictor variables.
Fig. 5

Comparison of thermal conductivity derived from experiment data (red), machine learning (blue) and Wiedemann–Franz law (green), respectively.

In the same manner, the thermal conductivity of other metals is predicted by the ML model including the electrical conductivity. In order to ensure the amount of data for training the ML model, two or more metals are grouped together at training and predicting procedures. We expect that such method can expand to more practical prediction for thermal conductivity of alloys as further application. Tables 3 and 4 show the contents of each group in this work and the electrical conductivities used for thermal conductivity prediction, respectively. Figures 6 and 7 show thermal conductivities of different groups of fcc and bcc metals at various temperatures predicted by the ML model and the corresponding experimental values. It can be seen from Fig. 6 that the ML model successfully predicts the thermal conductivities of Al, Cu, Ag, Au when making use of results of MD simulations with both potentials. However, for other fcc metal groups containing Ni, Pd or Pt, values of thermal conductivity predicted by the ML model often deviate from experimental values for both the cases of EAM and LJ potentials.

Table 3 Content of metal groups for regression.
Table 4 Electrical conductivity of metals used in predictor variables.40)
Fig. 6

Comparison of predicted and experimental values for thermal conductivities of various groups of fcc metals at various temperatures predicted by the machine learning model with electrical conductivity.

Fig. 7

Comparison of predicted and experimental values for thermal conductivities of various groups of bcc metals at various temperatures predicted by the machine learning model with electrical conductivity.

Similarly, among bcc metal groups in Fig. 7, thermal conductivities of Fe and W predicted by the ML model agree well with corresponding experimental values both for EAM and MEAM potentials. On the other hand, when the related data of V is added to the group of Fe and W, the ML model fails to predict the thermal conductivities of them. Furthermore, for the case of Fe, W and Cr, relatively good predictions of thermal conductivity can be obtained when using two metals among them whereas the group of all the three metals leads to unsuccessful prediction.

From the regression analysis, thermal conductivity of Al, Cu, Ag, Au (fcc metals), Fe, W and Cr (bcc metals) can be predicted well by the ML model while thermal conductivity of Ni, Pt, Pd (fcc metals) and V (bcc metal) cannot be predicted exactly in the same manner. It is found that experimental values of the thermal conductivity of Al, Cu, Ag, Au, Fe and W, which are predicted well the by ML model, decrease monotonically with increasing temperature (i.e., the negative temperature dependence) as shown in Figs. 2 and 3. On the contrary, experimental values of the thermal conductivity of Ni, Pd, Pt and V tend to increase with increasing temperature (i.e., the positive temperature dependence) at the high temperature range and thermal conductivity of these metals are not predicted well by the ML model. Therefore, it is expected that the temperature dependence of thermal conductivity is a key factor for the precise prediction.

3.3 Performance of machine learning model

In order to evaluate the performance of ML model quantitatively, the coefficient of determination (R2) is calculated as follows:   

\begin{equation} R^{2} = 1 - \frac{\displaystyle\sum\nolimits_{i}(y_{i} - f_{i})^{2}}{\displaystyle\sum\nolimits_{i}(y_{i} - \bar{y})^{2}} \end{equation} (9)
where fi is thermal conductivity predicted by the ML model, yi is the corresponding experimental values, and $\bar{y}$ is the average of the experimental values. According to the equation, R2 value ranges from −∞ to 1 (mostly in the range of 0 to 1) and it represents that the closer R2 value is to 1, the smaller the relative residual error. The calculated R2 values of all cases of fcc metals are summarized in Fig. 8. From Fig. 6(a-1) and (a-2), it is observed that in Al–Cu–Ag–Au group, thermal conductivities of each metal predicted by the ML model agree well with their experimental values, and the R2 values of them are all close to 1, which indicates high accuracy of the prediction. In Fig. 6(b-1) and (b-2), thermal conductivity of Ni predicted by the ML model decreases with increasing temperature while the experimental value increases when increasing temperature at the high temperature range. The R2 value of Ni is very low and R2 of other metals are also lower than that of Fig. 6(a-1) and (a-2). In general, a ML model makes prediction based on the learned features from input data. In Al–Cu–Ag–Au group, the temperature dependence of thermal conductivity of these four metals is similar to each other and this common feature is correctly learned by the ML model. However, in Al–Cu–Ag–Au–Ni group, negative temperature dependence obtained from the regression analysis is found in all elements because negative dependence by the influence of Al, Cu, Ag and Au is dominant. This negative trend is learned dominantly and cause a prediction of negative temperature dependence for Ni as well by the ML model although it is different from the trend in experimental values. Similarly, for Pd in Fig. 6(c-1) and (c-2), and Pt in Figure (d-1) and (d-2), due to the different temperature dependency of thermal conductivity compared with Al, Cu, Ag and Au, it is difficult for the ML model employed here to learn the correct feature of thermal conductivity of all metals at the same time. As a results, it caused the undesirable prediction with negative values of R2 for Pd and Pt. In general, it is not straightforward to find physical background in results of regression by machine learning. We can infer that the performance for Ni is better than those of Pt and Pd because the positive temperature dependence in experimental values of Ni appears only at high temperature region whereas such positive temperature dependence is found in almost all temperature region for Pt and Pd as shown in Fig. 2. However, it is also possible that this is due to the characteristics or defect of the linear regression model.

Fig. 8

Summary of the coefficient of determination (R2) for the regression analyses of fcc metals in Fig. 6.

The calculated R2 values of all cases of bcc metals are summarized in Fig. 9. In all combinations of two bcc metals shown in Fig. 7(a), (b) and (c), thermal conductivities of each metal predicted by the ML model agree well with their experimental values, and the R2 of them are all close to 1. However, the performance of the ML model decreases when predicting thermal conductivity of metals in the group of various temperature dependences as shown in Fig. 7(d-1) and (d-2). On the other hand, the Fe–W–Cr group in Fig. 7(e) makes a successful prediction for only W and failed for Fe and Cr. It is considered that the different decreasing rate of thermal conductivity with respect to temperature makes it difficult for the ML model to learn suitable features for prediction for all the three metals but possible for two of them, as shown in Fig. 7(a-1), (a-2), (b) and (c) because only the negative temperature dependence is needed to be learned.

Fig. 9

Summary of the coefficient of determination (R2) for the regression analyses of bcc metals in Fig. 7.

4. Conclusion

In this study, MD simulations with commonly used interatomic potentials are performed to calculate the thermal conductivity of representative metals at various temperatures. All results from MD simulations are found to deviate significantly from the experimental values. It is because that the thermal conductivity of metals consists of components contributed by phonons and electrons respectively while the MD simulation is inherently insufficient to simulate the contribution of electrons in the thermal conductivity. To make up for the deficiency of the MD simulation, the ML model is employed to predict accurate values of thermal conductivity using deficient MD results. A regression model, Lasso with predictor variables including results of MD simulations, temperature and electrical conductivity successfully predicts the thermal conductivity of metals with negative temperature dependence. Our prediction is better than the estimation from the Wiedemann–Franz law, which is the advantage of this study compared to the conventional estimation of the effect of electrical conductivity on the thermal conductivity of metallic materials. Although the current ML model has room for improvement to predict thermal conductivity of metals with various temperature dependence, it shows us new possibilities of new ML approach for improving the accuracy of physical properties obtained from MD simulation.

Acknowledgments

This work was supported by Grant-in-Aid for Scientific Research (B) (No. 22H01754) from Japan Society for the Promotion of Science (JSPS).

REFERENCES
 
© 2023 The Japan Institute of Metals and Materials
feedback
Top