Thermal conductivity calculation in anisotropic crystals by molecular dynamics : application to alpha-Fe2O3

TL;DR

This paper develops a molecular dynamics methodology to accurately calculate the thermal conductivity of anisotropic hematite crystals, emphasizing size effects, directional dependence, and temperature behavior in simulations.

Contribution

It introduces an adapted potential and a comprehensive approach for calculating anisotropic thermal conductivity in hematite using non-equilibrium molecular dynamics.

Findings

Both longitudinal and transverse sizes affect thermal conductivity results.

Thermal conductivity varies with crystallographic direction, indicating anisotropy.

Non-linear temperature dependence of thermal conductivity is observed in simulations.

Abstract

In this work, we aim to study the thermal properties of materials using classical molecular dynamics simulations and specialized numerical methods. We focus primarily on the thermal conductivity \kappa using non-equilibrium molecular dynamics to study the response of a crystalline solid, namely hematite (alpha-Fe2O3), to an imposed heat flux as is the case in real life applications. We present a methodology for the calculation of \kappa as well as an adapted potential for hematite. Taking into account the size of the simulation box, we show that not only the longitudinal size (in the direction of the heat flux) but also the transverse size plays a role in the determination of \kappa and should be converged properly in order to have reliable results. Moreover we propose a comparison of thermal conductivity calculations in two different crystallographic directions to highlight the spatial anisotropy and we investigate the non-linear temperature behavior typically observed in NEMD methods.

Tables

Table 1. Table 1: Crystal structure validation

Cell parameters.	Experimental (Å)	Calculated.111calculated values from Pedone et al. (2006) (Å)	Calculated222this work (Å)	Difference (%)
a	5.039	4.95	5.066	0.54
b	5.039	4.95	5.066	0.54
c	13.77	13.42	13.74	0.23
c/a	2.73	2.71	2.71	0.77

Equations

$${\overline{\mathbf{v}}}_{\mathbf{i}}={\mathbf{v}}_{\mathbf{G}}+\alpha ({\mathbf{v}}_{\mathbf{i}}-{\mathbf{v}}_{\mathbf{G}})$$

$$\alpha =\sqrt{1+\frac{\Delta ϵ}{E_{P_{+}}}}.$$

$$E_{P_{+}}=\frac{1}{2}\sum _i{m}_i{{\mathbf{v}}_{\mathbf{i}}}^2-\frac{1}{2}\sum _i{m}_i{{\mathbf{v}}_{\mathbf{G}}}^2.$$

$$\vec{J}=-\kappa \vec{\nabla }T$$

$$U(r)=D_{ij}([1-e^{-a_{ij}(r-r_0)}{]}^2-1)+\frac{C_{ij}}{r^{12}}+\frac{z_i{z}_j{e}^2}{r}$$

$$\frac{C_{xx}}{C_{zz}}=\frac{362}{326}\sim 1.1$$

Full text

thermal conductivity, heat transport, NEMD, molecular dynamics, anisotropy, hematite

I Introduction

From ab initio calculations of material properties to finite-element models of macroscopic systems, numerical studies of the thermal transport are motivated by the wide range of technological applications. Examples of such applications are found in nanoelectronics Sinha and Goodson (2005), aerospace Lee et al. (2012), automotive Yang (2005) and building sector Jelle (2011). Therefore, being able to model the heat transfer and predict thermal properties presents a definite advantage for designing or improving materials and industrial processes.

NEMD methods have been applied to thermal conductivity calculations since the beginning of the 1980s Evans (1982) and a number of different algorithms have been developed Schelling, Phillpot, and Keblinski (2002). Two main approaches can be identified: imposing the temperature gradient or imposing the heat flux. The latter approach has the advantage of a faster convergence Heino and Ristolainen (2003) and two of its algorithms are now commonly applied Schelling, Phillpot, and Keblinski (2002): the particle velocity exchange presented in  Müller-Plathe (1997) and the heat source and heat sink method described in  Jund and Jullien (1999),  Ikeshoji and Hafskjold (1994) and  Aubry et al. (2004). Here we make use of the source and sink method which was first developed for amorphous materials and simple Lennard-Jones (LJ) liquids and has been generalized to simple crystals, e.g. FCC Cu, Ag, Au, etc. in  Bracht et al. (2014); Ge and Chen (2013); Zhang, Hu, and Tang (2013); Zhou, Anglin, and Strachan (2007) and, in some cases, to more complex structures such as tetragonal ZrO2Schelling and Phillpot (2001) or zeolitic imidazolate frameworks in  Zhang and Jiang (2013).

The work reported here is part of a wider effort which aims to predict the material-lubricant interactions involved in the operation of a car engine by determining numerically the thermal conductivity of the different parts. As such, we focused in a first step on developing a methodology that we then applied to iron oxide (hematite) as a model material for the surface of the solid parts of the engine and we considered a temperature range between 300 and 500 K. In addition to its industrial and technological importance Dzade, Roldan, and de Leeuw (2014), in recent years hematite has been subject to a renewed interest due to discoveries concerning the geological structure and mineral properties of Mars Christensen et al. (2000); Morris et al. (2006); Fraeman et al. (2013). Bulk hematite ($\alpha\textrm{-Fe}_{2}\textrm{O}_{3}$), of space group $R\overline{3}c(167)$ , has a crystal structure that can be indexed as hexagonal with a unit cell consisting of 6 formula units in which the oxygen ions lie approximately in a hexagonal close-packed framework while the iron ions are positioned symetrically in two-thirds of the octahedral interstices Morrish (1994); Nolze and Winkelmann (2014). We chose to model this material using an empirical interatomic potential that we adapted and tested.

The main objective of the study was to put together a set of methods and tools for the calculation of the thermal conductivity of a crystal based on a previous work on amorphous materials by non equilibrium molecular dynamics (NEMD) Jund and Jullien (1999). We considered the influence of sample size, temperature and crystallographic orientation and we present a step by step methodology.

The paper is divided into three sections. In the first section we describe the NEMD algorithm and methodology, the empirical potential used to model the material and the numerical methods applied to circumvent sample size effects. In the second part we present and discuss our results on the thermal conductivity of hematite single crystals, its temperature dependence, the spatial anisotropy and the nonlinear behavior observed in specific regions of the simulation box. Finally, we draw the major conclusions in the last section.

II Numerical methods

II.1 Molecular dynamics and NEMD

Classical molecular dynamics (MD) is a very popular method providing atomic-scale information on material properties and processes. Based on the integration of the Newton equations of motion at the atomic level, it allows to use relatively large simulation boxes compared to first-principles calculations. In non equilibrium molecular dynamics (NEMD) one or several internal variables are constrained to keep the system out of the equilibrium state. In our case we use the standard velocity-Verlet algorithm for time integration with a time step of 0.6 fs and we apply the method described in  Jund and Jullien (1999) for the determination of the thermal conductivity $\kappa$ . The first step is choosing the direction of the heat flux. Here we will start by considering the z direction, parallel to the [001] crystallographic direction in the hexagonal lattice of hematite. Then we define two plates positioned at 1/4 and 3/4 of the simulation box length and orthogonal to the z axis. Periodic boundary conditions (PBC) are applied in every direction. At each time step, one of these plates receives a fixed amount of energy $\Delta\epsilon$ while the same amount is subtracted from the other. This is done by constantly rescaling the velocities of the atoms within the plates. This effectively results in imposing a heat flux across the box and parallel to the z direction. Care is taken to prevent a drift of the center of mass caused by the velocity rescaling. For a particle $i$ in the hot plate $P_{+}$, the rescaling can be expressed as

[TABLE]

where $\mathbf{\overline{v}_{i}}$ is the updated velocity, $\mathbf{v_{G}}$ the velocity of the center of mass of the ensemble of particles in $P_{+}$ and

[TABLE]

Of course, in the case of the cold plate $P_{-}$, $\Delta\epsilon$ should be subtracted. The non-translational kinetic energy $E_{P_{+}}$ is given by

[TABLE]

To keep track of the temperature profile along the z axis, the box is divided into a set of slices in which the local temperature is computed, as summarized in figure 1. The width $a$ of those slices, and therefore the number of atoms in each of them, affects the temperature calculation. A larger number of atoms provides better statistics and so allows to reduce the time over which the values need to be averaged. But to obtain a detailed temperature profile the number of slices has to be large enough and therefore their width is limited. In order to eliminate this constraint we introduce several sets of overlapping bins, as described in figure 2. This allows to increase the temperature profile resolution without decreasing the number of atoms per slice. We find that 4 sets of 12 slices is a good compromise between precision and computational effort.

The initial state of the simulated system is obtained by replicating a smaller cell previously optimized by energy minimization. After a first run of relaxation in a (NPT) ensemble at the desired temperature, the heat flux is switched on. We have checked that small variations of the amount of energy added to or subtracted from the plates don’t significantly affect the value of $\kappa$. Nevertheless, large values of $\Delta\epsilon$ will cause a large temperature difference between the source and the sink. A value of 1.5 % of $k_{B}T$ seems reasonable for a simulation box containing 60000 atoms. Once the heat flux is switched on, it is necessary to let the dynamics run long enough to produce a steady state. Typically, this stationarization is of the order of 1 ns as reported in  Schelling, Phillpot, and Keblinski (2002); Maiti, Mahan, and Pantelides (1997) for Stillinger-Weber silicon. We observe an average time of 0.8 ns before stable temperatures are observed in the source and sink plates. Then we start collecting and averaging the temperature values in every bin. Depending on the size of the system, we find that the averaging time necessary to remove the numerical noise and obtain a smooth temperature profile, on which a linear function can be fitted, varies from 2 to 6 ns. Except for the first relaxation run, all the calculations are performed in a (NVE) microcanonical ensemble. Finally, the temperature gradient is calculated from the time-averaged temperature profile and the Fourier’s law (4) is used to obtain $\kappa$. Recently, an adapted version of the heat exchange algorithm was proposed by Wirnsberger et al. Wirnsberger, Frenkel, and Dellago (2015) to correct a large drift of the total energy of the system observed in long-term simulations. However, in the present work we observe a moderate difference between the initial and final values of the total energy (of the order of $10^{-4}$ %). This is due to the length of our simulations which is less than 10 ns and to the relatively small value of the time step chosen to ensure energy conservation.

[TABLE]

The application of such a method takes its toll on numerical performance. Indeed the permanent temperature rescaling in the heat source and heat sink and the local temperature computation in each slice are time consuming. And with the constraints of stationarization and averaging, the time cost is a primary issue. Some of the simulations performed as part of this work required more than 100000 CPU hours individually. All molecular dynamics simulations were conducted using a modified version of the LAMMPS package Plimpton (1995) (based on version may 2015).

II.2 Interaction potential

To perform molecular dynamics on an ionic crystal such as hematite one needs to apply an interatomic potential defined to model the structure of the material as well as its physical properties. As mentioned in section I, we consider a hexagonal unit cell with 6 formula units (30 atoms) and the following lattice constants, respectively a, b and c: $5.039$, $5.039$, $13.77$ ÅSadykov et al. (1996). We use a modified version of the potential $U(r)$ published in 2006 by Pedone et al. Pedone et al. (2006) to model our material. Three terms contribute to the equation of $U$: a short-range Morse potential, a short-range $r^{-12}$ repulsion and a long range Coulomb interaction. A cutoff distance of 7 Å was applied to the Morse and short-range repulsion terms and the Coulomb interaction was implemented by way of an Ewald summation.

[TABLE]

The parameters of the potential were initially fitted on the experimental lattice constants, atomic positions and elastic constants using free energy minimization Catlow (1997) and other empirical fitting methods implemented in the GULP package Gale (1997, 1996); Bush et al. (1994). These methods allow to fit the potential on a structure relaxed at a finite temperature taking into acount quantities such as mechanical and dielectric properties. Nevertheless the experimental lattice constants used for this fit are for some reason different from the values obtained in most of the studies Sadykov et al. (1996); Sorescu et al. (2004); Sawada (1996); Barinov et al. (1995); Morrish (1994); Maslen et al. (1994). A fine-tuning of the Morse equilibrium parameter ($r_{0}$) for the Fe-O and the O-O interactions allowed us to accurately fit the “correct” experimental values. Indeed differences between lattice constants and interatomic distances experimentally observed and obtained with the modified potential are well under 1%, as presented in table 1. The new values of $r_{0}$ are 2.41810 Å and 3.65455 Å respectively for the Fe-O and O-O interactions.

In addition to the structural properties, the bulk modulus was investigated in two different ways in order to validate further the modified potential function. First the elastic constants were calculated by deforming the hexagonal box in the 6 degrees of freedom and observing the change in the stress tensor at zero temperature. The Voigt-Reuss approximation applied to the calculated elastic tensor gives a value of 217 GPa for the bulk modulus. Additionnally, we performed single-point energy calculations at different volumes around the equilibrium volume of the cell and fitted the curve E=f(V) with the Vinet equation of state Vinet et al. (1987). We obtained this way a bulk modulus of 219 GPa. These calculated results fall within the measured values of 203 Liebermann and Schreiber (1968) and 230 GPa Catti, Valerio, and Dovesi (1995) available in the literature. The small difference between the two calculated values (less than 1%) may be used as an indication of the precision of this kind of calculations.

II.3 Finite size effects

With a band gap of 2.1 eV Morrish (1994), the heat transport in hematite is primarily described by the lattice thermal conductivity which, in turn, is determined by the phonon transport. When performing NEMD on a finite simulation box, the phonon mean free path (MFP), which is the average distance a phonon travels before being scattered, is of particular importance. Phonons with different MFPs contribute differently to the thermal transport and in some materials the MFP of the contributing phonons can take values larger than several hundred nanometers Sellan et al. (2010); Hermet and Jund (2016). The NEMD method, as described in the previous sections, requires that two plates be defined at one quarter and three quarters of the box length to serve as heat source and sink. It would therefore be necessary to have a simulation box twice as long as the length of the largest contributing MFP. These scales are very difficult to handle in molecular dynamics simulations since the simulations would require enormous computing resources and very long execution times. We studied the effect of finite box dimensions in the direction perpendicular to the heat flux and we applied the method proposed by Schelling et al. in  Schelling, Phillpot, and Keblinski (2002) to the finite size effects parallel to the heat flux. For the latter, several simulations with increasing box sizes $L$ are performed and their results are combined with the use of the so-called linear extrapolation procedure. As the size of the simulation box grows, more phonon mean free paths are taken into account which increases the calculated value of the thermal conductivity. As described in  Sellan et al. (2010), within a first order approximation a linear dependence can be expected between $1/\kappa$ and $1/L$. As a consequence, $\kappa_{\infty}$ can be evaluated by plotting $1/\kappa$ against $1/L$ and extrapolating the curve to $1/L_{\infty}=0$ with a linear function. To obtain a reasonable precision of the curve to be fitted, a large number of individual simulations has to be performed. Nevertheless, the computational effort needed is considerably lower than for a real-size simulation with $L$ larger than the largest phonon mean free path. This method has also the advantage of providing a thermal conductivity value that relies on a significant number of different simulations with different initial states, thus reducing the statistical errors that could be attributed to an individual simulation.

Sellan et al. Sellan et al. (2010) mentioned a number of limitations of this method, due primarily to the first order approximation. For example, when considering a large collection of different box sizes, they reported good results for the application of the extrapolation method to Lennard-Jones argon, while an underestimation was observed in the case of Stillinger-Weber silicon. In the latter case, the calculated values and the extrapolated curve would diverge for the largest sizes. To address this issue, Hu et al. Hu, Evans, and Keblinski (2011) suggested a relation between the divergence and the aspect ratio of the simulation box (length/width). Studying Lennard-Jones solid argon, Lennard-Jones WSe2 and graphite, their first conclusion is that a divergence can be observed at very high aspect ratios (200 to 300) for an elastically isotropic material such as LJ argon. Nevertheless these authors add that such high limits don’t have practical consequences since a converged value of $\kappa$ can be computed well below those limits. Their second conclusion is that a clear divergent behavior is observed for low aspect ratios ($\sim$ 30) for elastically anisotropic materials such as LJ WSe2 and even earlier for graphite. With the smallest lateral size considered for the initial tests, we reached a maximum aspect ratio of 133. With the lateral size used in most of our simulations the largest aspect ratios range from 83 to 116. Since no divergence can be seen in the results presented hereafter, we can conclude that our material is probably elastically isotropic. Moreover, Hu et al. proposed a criterion to predict if the thermal conductivity of a material can be modeled by NEMD methods or not. This criterion is based on the ratio of the elastic constants along two crystallographic directions ($\frac{C_{xx}}{C_{zz}}$) which should be reasonably low (< 5). We computed the 6x6 matrix of the elastic constants as described in section II.2 and obtained the following ratio:

[TABLE]

which shows that our material is perfectly suited for NEMD methods and is not prone to the divergence issues described in Hu, Evans, and Keblinski (2011).

III Results and discussion

III.1 Size-dependent simulations

As explained in the previous section, the calculation of $\kappa$ requires two steps. First a number of size-dependent simulations, then a size-independent extrapolation. Figure 3 shows the typical result for an individual simulation for a given box size $L$. The time-averaged temperature profile is presented along with the position and width of the slabs where heat is added or subtracted (heat source and heat sink). In the central part of the simulation box, the temperature exhibits a linear profile which can be fitted to calculate the thermal gradient. Non-linear areas are noticeable close to the heat source and sink.

In figure 4 the thermal conductivity at 300 K is shown as a function of the box size in the directions both orthogonal and parallel to the heat flux. One can observe a quick convergence of $\kappa$ as a function of the lateral size, and we fixed this size to 30 Å for the rest of the study. As for the evolution of $\kappa$ as a function of $L$, the size parallel to the heat flux (horizontal axis in figure 4), it requires the application of the extrapolation method discussed in section II.3.

A plot of $1/\kappa$ against $1/L$ is presented in figure 5 showing the calculated values and the extrapolated curve with a logscale x axis. The corresponding value of $\kappa$ is calculated to be 16 W.m*-1*.K*-1* for a pure single crystal of hematite at 300 K in the [001] crystallographic direction. Several thermal conductivity measurements of polycrystalline hematite at room temperature are reported in the litterature. The values presented in  Diment and Pratt (1988); Clark (1966); Horai (1971) range between 11 and 13 W.m*-1*.K*-1* while Akiyama et al. report a much larger value of 17 W.m*-1*.K*-1* in  Akiyama et al. (1992). However, experimental values of $\kappa$ for single-crystal hematite samples are less common. In  Clauser and Huenges (1995), references are made to a 1974 work Dreyer (1974) where a value of 12.1 W.m*-1*.K*-1* was reported for a single crystal in the [001] direction and 14.7 W.m*-1*.K*-1* (22 % more) for a second direction orthogonal to [001].

III.2 Temperature dependence

In order to determine the temperature dependence of $\kappa$ in the range 300 to 500 K, it is in principle necessary to apply the previously described procedure for several temperatures. But taking into account the computational cost of such a brute force method, we propose an optimized approach where a full study is conducted at 300 K and 500 K but only a partial analysis is done in between. Applying the full extrapolation procedure at the limits of the range provides a validation of its applicability. Assuming that the thermal conductivity follows the same behavior at the intermediate temperatures, we perform only a limited number of simulations at the lowest and highest box sizes for those temperatures providing thus the start and end points for the extrapolation. The results are coherent as can be observed in figure 6 and validate the approach. We estimate that proceeding in this way required 30 to 40 % less computational time than a full study.

From the infinite box size extrapolations, the temperature dependence of the thermal conductivity can be assessed. The results are presented in figure 7 along with a comparison with measurements made on polycristalline hematite in  Akiyama et al. (1992). As expected, $\kappa$ decreases with the temperature in the investigated range. In this temperature range it is reasonable to think that Umklapp processes are dominant and thus $\kappa$ should decrease like 1/T: this behavior is not obvious from the calculated values of figure 7. However, we observed that the individual, size-dependent, conductivity values such as the ones shown in figure 4 may vary by up to 7.5% when differences are introduced in the initial state of the simulations (e.g. random initial velocity distributions). And even though the extrapolated, size-independent, values of $\kappa$ rely on several different simulations, the lack of precision makes it difficult to assess a precise mathematical function from the curve of figure 7.

III.3 Spatial Anisotropy

The results presented in the previous sections were obtained for a heat flux in the [001] crystallographic direction. To investigate the spatial anisotropy of $\kappa$, we applied the same methods to the [100] direction at 300 K. Figure 8 shows the evolution of $\kappa$ as a function of system size together with the extrapolated curve. A spatial anisotropy is observed for hematite since we obtained a value of 20 W.m*-1*.K*-1* for $\kappa$ with the heat flux orientated in the [100] crystallographic direction. This is 25 % more than in the [001] direction which is consistent with the experimental values mentioned earlier in section III.1.

III.4 Non-linearity

As can be observed in figure 3, the temperature profile exhibits a nonlinear behavior near the source and the sink. This behavior has been partially explained by phonon scattering in previous studies Schelling, Phillpot, and Keblinski (2002); Sellan et al. (2010). In order to calculate the temperature gradient and apply Fourier’s law, it is thus necessary to consider the profile far enough from those non-linear areas. Our experience led us to choose an “exclusion” distance equal to the width $a$ of the temperature bins (heat source and sink included). Indeed, we found that this specific choice was appropriate in most of the simulations performed in this work, independently of the actual value of $a$. Moreover, investigating the evolution of the non-linearity as a function of time, we find that in the case of long enough simulations this non-linearity decreases significantly even after the system is considered to have reached the steady state. Figure 9 shows the evolution of the time-averaged temperature profile in a 140 nm long simulation box made of a 6x6x105 supercell (113000 atoms) during simulations lasting up to 7.8 ns.

The area between the temperature curve and the linear regression (colored area in figure 9) was calculated to quantify the non-linearity. Figure 10 shows the evolution of this quantity, starting after the stationarization. The value peaks around 3 ns and decreases afterwards, until a minimum value is reached. We have observed this evolution with time for different system sizes. Thus it appears that the non-linearity of the temperature gradient close to the source and the sink can partly be characterized as a slow transient phenomenon.

III.5 Conclusion

A detailed methodology has been presented for the determination of the thermal conductivity of crystals and applied to pure single-crystalline hematite. Calculated values for different temperatures are in reasonable agreement with available experimental data on polycrystals with values ranging from 16 to 10 W.m*-1*.K*-1* between 300 and 500 K. Moreover, an investigation of the spatial anisotropy has been undertaken and shows, at 300 K, a 25% increase of the thermal conductivity in the [100] direction with respect to the [001] direction in agreement with measurements on single crystals. Finally, some specific elements of the calculation procedure, such as the width of the temperature bins or the nature of the nonlinear behavior, have been analyzed highlighting new aspects in the application of the NEMD scheme for the determination of the thermal conductivity of crystalline solids.

Acknowledgements.

This work was granted access to the HPC resources of CINES under the allocation 2016-c2016097598 made by GENCI. We gratefully acknowledge Total S.A. and Total M.S. for their financial support and we thank Sophie Loehlé and Vincent Lacour for fruitful discussions.