**1. Introduction**

16 Will-be-set-by-IN-TECH

72 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Radoul, M.; Sundararajan, M.; Potapov, A.; Riplinger, C.; Neese, F.; Goldfarb, D. (2010).

Rezá ˘ c, J.; Fanfrlík, J.; Salahub, D.; Hobza, P. (2009) ˘ *Journal of Chemical Theory and Computation*

Rosas-Garcia, V. M.; Lopez-Cortina, S. T.; Chavarria-Resendez, A. *Journal of Molecular*

Rosas-Garcia, V. M.; Saenz-Tavera, I.; Cantu-Morales, D. accepted manuscript,

Saiz, E.; Guzmán, J.; Iglesias, M. T.; Riande, E. (1996). Dynamics and Polarity of Substituted

Saunders, M.; Houk, K. N.; Wu, Y. D.; Still, W. C.; Lipton, M.; Chang, G.; Guida, W. C. (1990).

Schenker, S.; Schneider, C.; Tsogoeva, S. B.; Clark, T. (2011) *Journal of Chemical Theory and*

Shiroishi, H.; Oda, T.; Hamada, I.; Fujima, N. (2005). Structure and magnetism of anion

Stewart, J. J. P. (1989). Optimization of parameters for semiempirical methods I. Method. *Journal of Computational Chemistry* Vol. 10, No. 2, (March 1989) 209–220 Stewart, J. J. P. (2007). Optimization of parameters for semiempirical methods V: Modification

Stewart, J. J. P. 2008. MOPAC2009. Colorado Springs, CO, USA: Stewart Computational

Takayanagi, T.; Yoshikawa, T.; Kakizaki, A.; Shiga, M.; Tachikawa, M. (2008). Molecular

Tosco, P.; Lolli, M. L. (2008). Hydroxy-1,2,5-oxadiazolyl moiety as bioisoster of the

Yang, Z.; Xiong, S.-J. (2008). Structures and electronic properties of small FeB*<sup>n</sup>* (n = 1-10) clusters. *Journal of Chemical Physics* Vol. 128, No. 18, (May 2008) 184310-8

*Computation* Vol. 7, No. 11, (October 2011) 3586–3595

*Modeling* Vol. 13, (September 9) 1173–1213

Chemistry. URL: http://OpenMOPAC.net

No. 1-3, (November 2008) 29–36

16, (December 2005) 1701–1718

1,3-Dioxacyclohexanes. *Journal of Physical Chemistry* Vol. 100, No. 47, (January 1996)

Conformations of cycloheptadecane. A comparison of methods for conformational searching. *Journal of the American Chemical Society* Vol. 112, No. 4, (February 1990)

iron oxide clusters, Fe*n*O*m*-(n=3,4). *Polyhedron* Vol. 24, No. 16-17, (November 2005)

of NDDO approximations and application to 70 elements. *Journal of Molecular*

dynamics simulations of small glycine-(H2O)*<sup>n</sup>* (n = 2-7) clusters on semiempirical PM6 potential energy surfaces. *Journal of Molecular Structure: THEOCHEM* Vol. 869,

carboxy function. A computational study on *γ*-aminobutyric acid (GABA) related compounds. *Journal of Molecular Modeling* Vol. 14, No. 4, February 2008, 279–291. Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. (2005).

GROMACS: Fast, flexible, and free. *Journal of Computational Chemistry*, Vol. 26, No.

*Structure: THEOCHEM* Vol. 958, No. 1-3, (October 2010) 133–136

Ponder, J. (2011). *TINKER* URL: http://dasher.wustl.edu/tinker/

No. 26, (May 2010) 7276–7289

Vol. 5, No. 7, (July 2009) 1749–1760

DOI:10.1007/s10876-011-0420-4

18345–18350

1419–1427

2472–2476

dynamics and free energy calculations to simulate the structural and energetic properties of molecules. *Comp. Phys. Commun.* Vol. 91, No. 1-3, (September 1995) 1–41

Revisiting the nitrosyl complex of myoglobin by high-field pulse EPR spectroscopy and quantum mechanical calculations. *Physical Chemistry Chemical Physics*, Vol. 12,

> This chapter presents an overview of the Molecular Dynamics (MD) simulation technique to predict thermal transport properties of nanostructured materials. This covers systems having characteristic lengths of the order of a few nanometers like carbon nanotubes, nanowires and also superlattices, i.e. composite materials made of submicronic thickness of solid layers. The common features of these systems is the small ratio between their characteristic system size and the phonon mean free path, which leads to ballistic heat transport and deviations from the classical Fourier law. Also when the density of interfaces gets large, the energy transport properties of the materials can not longer be described solely by the thermal conductivities of the constituents of the material, but depend also on the thermal boundary resistance which measures the transmission of phonons across an interface. In this context, molecular dynamics was proven to be a very useful technique to study heat transport in nanostructured materials. The main reasons are; the length scale probed by the method is in the nanometer range, and it does not make any assumption on the phonons dynamics except their classical nature.

> In this contribution, we present two MD methods, the equilibrium and the non-equilibrium method, which are now commonly used to determine both the thermal conductivity and the thermal boundary resistance of nanostructured materials. We focus on superlattices and discuss how the structural features of the interfaces like height, shape, inter-diffusion phenomena and the layer thickness affect the thermal conductivity of the superlattice. We show how these complex phenomena can be predicted by simple models of Lennard-Jones crystals with a mass ratio corresponding to the acoustic impedance ratio of Si/Ge and GaAs/AlAs superlattices.

#### **2. Molecular dynamics**

The development of molecular simulations began in the early fifties after the considerable development of computer facilities in the United States during World War II. A few years after the first Monte-Carlo simulation, Alder and Wainwright, first introduced in the late 50's the Molecular Dynamics (MD) method (Alder & Wainwright, 1957, 1959). The aim of the first simulations both Monte-Carlo or MD was to probe the different phases of model hard spheres. The necessity to model liquids motivated the development of realistic potentials. Rahman was the first to model Argon using Lennard-Jones potential , which is still considered as a standard potential for MD (Rahman, 1964). This opened the way to consider a broad range of condensed matter systems, ranging from liquid water first modelled by Rahman and Stillinger to silicon (Rahman & Stillinger, 1974). The efforts towards realistic modelling have permitted to applicate MD to characterize the collective excitations in solids (Hensen 1976). The last step came with the implementation of MD thermostats opening the way to probe the phonon dynamics in rare-gas solids (Ladd 1986). More recently, there has been an increasing number of MD studies on nano-scale heat transport motivated by fine measurements of energy transport in nano-materials (Volz 1999). We refer the reader to the classical textbooks on MD (Allain & Tildesley 1987 and Frenkel & Smit 1996) for a thorough introduction.

MD is a simulation method based on the numerical integration of Newton's equation of motion:

$$
\rho m \frac{d^2 \vec{r}\_i}{dt^2} = -\frac{\partial V}{\partial \vec{r}\_i} \tag{1}
$$

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 75

systems or reach orders of magnitude longer times . The system represented in a MD simulation has microscopic dimensions, and usually one is not so much interested in simulating a nano system with free surfaces, but rather a part of a bulk system. To embed the small system simulated in a bulk-like system, one can use *periodic boundary conditions*. These conditions assume that the system is repeated periodically in all space directions and thus any particle close to the boundaries of the simulation box interacts with an image particle, and when a particle crosses one of the face of the central simulation box, one of its images will enter the central box through the opposite face. One of the consequences of the use of periodic boundary conditions is the cut-off of long wavelengths fluctuations, i.e. those having a wavelength larger than the simulation box length. In the context of heat transport

The strength of MD lies on the versatility and the flexibility in the choie of the atomic potential (eq.1). Generally speaking, the total potential can be decomposed in the sum

12 3 ... *ii i j i <sup>j</sup> <sup>k</sup> i i j i <sup>j</sup> <sup>k</sup> V r = V r + V r ,r + V r ,r ,r +*

where the first term represents the one body contribution to the potential and is often associated to an external field like e.g. an electric field. The second and third terms represent respectively two-body and three-body potentials of the inter-atomic interaction. The choice of the inter-atomic potential depends on the problem at hand (type of atomic structure and on the properties to be studied), and on the necessity to model a real system or an experiment. One of the simplest potential used in condensed matter physics is the Lennard-Jones (LJ) potential, where the total potential decomposes in a sum of pair potentials:

> <sup>12</sup> <sup>6</sup> <sup>2</sup> *Vr= r r ij* 4/ / *ij ij*

This is the LJ pair potential, which depends only on the distance between neighbouring

atoms of the simulation box interact. From a practical point of view however, it would be computationally costly to estimate <sup>2</sup> *N* / 2 inter-atomic forces at each time step, with *N* 10000 . To bring back the computational cost of the force calculation to *O N* , the inter-

represents typically a neighbourhood of 50 atoms per particle. It should be mentioned at this point, that even if the potential at the radius of truncation is small compared to the well

thermodynamics of the system modelled. As an example, let us estimate the effect of truncating the inter-atomic potential respectively on the internal energy per particle and on

 is the atomic diameter. Although very simple, the LJ potential is often thought to be a good potential to model rare gas and in particular Argon, using the following parameters:

atomic pair potential is often truncated at a cut-off radius, usually 2,5 *cr =*

 

*= m* 3,40 10 . Due to the algebraic decay of the LJ potential, all the

, truncating the potential may have a non negligible effect on the

(5)

<sup>2</sup> *ij i j V= V r* (6)

identifies with the depth well of the potential, and

, which

(7)

simulations this will imply that phonons with long wavelengths are not present.

with

particles *ij ij r =r* . The interaction energy

 1,67 10 *J* and <sup>10</sup> 

the pressure. The respective errors are given by:

21

depth: 0,016 *V rc*

where: *m* , *ir* and *V* denote respectively the mass of the particle i, its position vector and the inter-atomic potential. The integration scheme commonly used is the Verlet algorithm (Verlet 1967a, 1967b) which predicts the positions at time *t + dt* given the positions at the earlier times *t* and *t dt* :

$$
\vec{r}\_i(t+dt) = 2\vec{r}\_i(t) - \vec{r}\_i(t-dt) + \frac{d^2\vec{r}\_i}{dt^2}dt^2 + O\left(dt^4\right) \tag{2}
$$

where the acceleration is calculated using the inter-atomic forces. The velocities can be computed using

$$
\vec{v}\_i(t) = \frac{\vec{r}\_i(t+dt) - \vec{r}\_i(t-dt)}{2dt} + O\left(dt^2\right) \tag{3}
$$

and the error is larger than the error on the positions, but as soon as the velocities are only used to compute the instantaneous kinetic energy the consequences are minor. Of course there are other integration schemes, but the Verlet algorithm has the advantage to be simple, easy to implement, stable and time-reversible. The typical time step used in the integration algorithms is a fraction of the characteristic atomic period

$$\tau = \sqrt{\frac{m\sigma^2}{k\_B T}} \approx 1ps \tag{4}$$

where *m* is the mass of the particles and denotes a typical atomic diameter.

Given the power of modern computers, it takes a few hours on a mono-processor machine to follow the trajectory of a set of 10000 particles over a time comparable to 1 ns. Using parallel MD codes or GPUs (graphics processing units) opens the way to model larger systems or reach orders of magnitude longer times . The system represented in a MD simulation has microscopic dimensions, and usually one is not so much interested in simulating a nano system with free surfaces, but rather a part of a bulk system. To embed the small system simulated in a bulk-like system, one can use *periodic boundary conditions*. These conditions assume that the system is repeated periodically in all space directions and thus any particle close to the boundaries of the simulation box interacts with an image particle, and when a particle crosses one of the face of the central simulation box, one of its images will enter the central box through the opposite face. One of the consequences of the use of periodic boundary conditions is the cut-off of long wavelengths fluctuations, i.e. those having a wavelength larger than the simulation box length. In the context of heat transport simulations this will imply that phonons with long wavelengths are not present.

The strength of MD lies on the versatility and the flexibility in the choie of the atomic potential (eq.1). Generally speaking, the total potential can be decomposed in the sum

$$V\left(\vec{r\_i}\right) = \sum\_{i} V\_1\left(\vec{r\_i}\right) + \sum\_{i$$

where the first term represents the one body contribution to the potential and is often associated to an external field like e.g. an electric field. The second and third terms represent respectively two-body and three-body potentials of the inter-atomic interaction. The choice of the inter-atomic potential depends on the problem at hand (type of atomic structure and on the properties to be studied), and on the necessity to model a real system or an experiment. One of the simplest potential used in condensed matter physics is the Lennard-Jones (LJ) potential, where the total potential decomposes in a sum of pair potentials:

$$V = \sum\_{i$$

with

74 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

the first simulations both Monte-Carlo or MD was to probe the different phases of model hard spheres. The necessity to model liquids motivated the development of realistic potentials. Rahman was the first to model Argon using Lennard-Jones potential , which is still considered as a standard potential for MD (Rahman, 1964). This opened the way to consider a broad range of condensed matter systems, ranging from liquid water first modelled by Rahman and Stillinger to silicon (Rahman & Stillinger, 1974). The efforts towards realistic modelling have permitted to applicate MD to characterize the collective excitations in solids (Hensen 1976). The last step came with the implementation of MD thermostats opening the way to probe the phonon dynamics in rare-gas solids (Ladd 1986). More recently, there has been an increasing number of MD studies on nano-scale heat transport motivated by fine measurements of energy transport in nano-materials (Volz 1999). We refer the reader to the classical textbooks on MD (Allain & Tildesley 1987 and

MD is a simulation method based on the numerical integration of Newton's equation of

*d r V m = dt r* 

the inter-atomic potential. The integration scheme commonly used is the Verlet algorithm (Verlet 1967a, 1967b) which predicts the positions at time *t + dt* given the positions at the

where the acceleration is calculated using the inter-atomic forces. The velocities can be

 <sup>2</sup> 2 *i i*

and the error is larger than the error on the positions, but as soon as the velocities are only used to compute the instantaneous kinetic energy the consequences are minor. Of course there are other integration schemes, but the Verlet algorithm has the advantage to be simple, easy to implement, stable and time-reversible. The typical time step used in the integration

> 2 1

*B m = ps k T* 

Given the power of modern computers, it takes a few hours on a mono-processor machine to follow the trajectory of a set of 10000 particles over a time comparable to 1 ns. Using parallel MD codes or GPUs (graphics processing units) opens the way to model larger

*r t + dt r t dt v t= + O dt dt*

<sup>2</sup> <sup>2</sup> *<sup>i</sup>*

*d r r t + dt = r t r t dt + dt + O dt dt*

*i ii*

*i*

algorithms is a fraction of the characteristic atomic period

where *m* is the mass of the particles and

*i*

and *V* denote respectively the mass of the particle i, its position vector and

2

2 4

(3)

(4)

denotes a typical atomic diameter.

(2)

(1)

2 2 *i*

Frenkel & Smit 1996) for a thorough introduction.

motion:

where: *m* , *ir*

computed using

earlier times *t* and *t dt* :

$$V\_2\left(\vec{r}\_{\vec{\eta}}\right) = 4\varepsilon \left( \left(\sigma \;/\ r\_{\vec{\eta}}\right)^{12} - \left(\sigma \;/\ r\_{\vec{\eta}}\right)^6 \right) \tag{7}$$

This is the LJ pair potential, which depends only on the distance between neighbouring particles *ij ij r =r* . The interaction energy identifies with the depth well of the potential, and is the atomic diameter. Although very simple, the LJ potential is often thought to be a good potential to model rare gas and in particular Argon, using the following parameters: 21 1,67 10 *J* and <sup>10</sup> *= m* 3,40 10 . Due to the algebraic decay of the LJ potential, all the atoms of the simulation box interact. From a practical point of view however, it would be computationally costly to estimate <sup>2</sup> *N* / 2 inter-atomic forces at each time step, with *N* 10000 . To bring back the computational cost of the force calculation to *O N* , the interatomic pair potential is often truncated at a cut-off radius, usually 2,5 *cr =* , which represents typically a neighbourhood of 50 atoms per particle. It should be mentioned at this point, that even if the potential at the radius of truncation is small compared to the well depth: 0,016 *V rc* , truncating the potential may have a non negligible effect on the thermodynamics of the system modelled. As an example, let us estimate the effect of truncating the inter-atomic potential respectively on the internal energy per particle and on the pressure. The respective errors are given by:

$$
\Delta \mathcal{U} = \rho \int\_{r\_c}^{\ast \infty} 4\pi r^2 V\_2(r) dr \approx 1,07\rho\sigma^3 \varepsilon \tag{8}
$$

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 77

As we have already briefly mentioned, the output of MD simulations can be compared to experiments if we know how to relate the microscopic state of the material under study to the macroscopic observables typically measured in an experiment. Generally speaking, one can distinguish three classes of information that can be extracted from MD. On the one hand, MD allows to compute global quantities which correspond usually to the thermodynamics of the system modelled. For instance, the internal energy can be

1 <sup>2</sup>

2

*i*

*B*

where the brackets denote an ensemble average. If sufficient time has been left to the system to reach equilibrium, the ensemble average may be performed by averaging over a sufficient long time thanks to ergodicity. On the other hand, in non-equilibrium situations, different initial conditions (usually the positions and the velocities of the atoms) may be used to generate independent trajectories in phase space. The latter definition of the kinetic energy

*mv Tr=*

*iVr*

where the sum runs over particles i in a small volume V(r) centred around r and which contains n(r) particles. Another important quantity in heat transport simulations is the

( ) 3

*i i <sup>i</sup>*

<sup>1</sup>

For simulations in crystals at low temperatures, the previous energy flux displays large oscillations due to optical phonons which carry a negligible amount of heat, and sometimes

<sup>1</sup> <sup>0</sup>

where the superscript "0" denotes the equilibrium positions of the atoms. From the knowledge of the local temperature and the energy flux, it is possible to measure thermal conductivity and we would come back to the measurement of transport coefficients in MD

*nr k*

is intimately related to the definition of the local temperature :

where *Ei* is the total energy of particle i, yielding for pair potentials:

it can be more suited to work with the equivalent definition of the flux:

in the next section, with a particular emphasis on heat transport.

*<sup>i</sup> <sup>i</sup> U = V r KE* (11)

<sup>2</sup> *<sup>i</sup> <sup>i</sup> KE = mv* (12)

(13)

*<sup>d</sup> J= Er dt* (14)

<sup>2</sup> *i i ij <sup>i</sup> <sup>j</sup> ij i i <sup>j</sup> J= Ev + F v +v r* (15)

<sup>2</sup> *ij <sup>i</sup> <sup>j</sup> ij i j J= F v +v r* (16)

expressed as:

and the total kinetic energy:

energy flux, defined as

and

$$
\Delta P = \frac{-\rho^2}{6} \int\_{r\_c}^{+\upsilon} r \frac{\partial V\_2}{\partial r} 4\pi r^2 dr \approx -1.07 \,\rho^2 \sigma^3 \varepsilon \tag{9}
$$

where we have used the previous value of the cut-off radius and we have assumed that beyond *rc* , the medium is structureless, i.e. the pair correlation function may be considered close to unity. The previous errors are found to be typically 10% the values of the internal energy and pressure for a condensed phase. This is not worrying since the calculation of thermodynamic quantities can be corrected using the previous estimations, but this has to be kept in mind when calculating a phase diagram for instance.

Although the LJ potential is quite simple, some situations require to use more sophisticated semi-empirical potentials whose parameters are chosen to reproduce either microscopic or macroscopic properties of a model system. These semi-empirical potentials are usually no longer pair potentials like the LJ potential but many-body, where the many-body terms describes how does the potential energy of an atom depends on its coordination.

As an illustration, let us mention the embedded atom model (EAM) used in the context of metals (Daw & Baskes 1984), where the potential

$$\sum\_{i$$

consists of a pair term *r* accounting for repulsion at short distances and a many body term *u r* which accounts for the cohesion of the metallic bond and which depends on the local density measured using the positions of neighbouring particles *i i <sup>j</sup> <sup>j</sup> r = wr* . The choice of the functional *u* and the weight function *w r* depends on the microscopic and macroscopic properties to be reproduced, usually the lattice constant, the elastic constants and the sublimation energy of the metal (Daw & Baskes 1984).

Another limitation of LJ potentials is the description of solids crystalizing in non-compact structures. Indeed, Lennard Jones atoms form close-packed structures at low temperatures (usually fcc at low pressure), and thus LJ potentials are also not adapted to model materials which crystallize in opened structures such as diamond. Semi-conductors like silicon or germanium are usually modelled using many-body semi-empirical potentials, in which the many body terms account for the local preferential bond ordering of the semi-conductor. Among the most famous potentials for silicon, let us quote the Tersoff potential (Tersoff 1986, 1988a, 1988b) and the Stillinger-Weber potential (Stillinger & Weber 1985). The parameters of the former are chosen so as to reproduce the elastic properties of bulk silicon, while the latter describes satisfactorily the structural properties of liquid silicon. Of course the list of potentials is non exhaustive and there exists plethora of empirical potentials to model as disparate systems as charged systems, liquid crystals, polymers, surfactants, granular media,... and mixtures of them!

As we have already briefly mentioned, the output of MD simulations can be compared to experiments if we know how to relate the microscopic state of the material under study to the macroscopic observables typically measured in an experiment. Generally speaking, one can distinguish three classes of information that can be extracted from MD. On the one hand, MD allows to compute global quantities which correspond usually to the thermodynamics of the system modelled. For instance, the internal energy can be expressed as:

$$\text{CLI} = \left\langle \sum\_{i} V\left(r\_{i}\right) \right\rangle + \text{KE} \tag{11}$$

and the total kinetic energy:

76 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

where we have used the previous value of the cut-off radius and we have assumed that beyond *rc* , the medium is structureless, i.e. the pair correlation function may be considered close to unity. The previous errors are found to be typically 10% the values of the internal energy and pressure for a condensed phase. This is not worrying since the calculation of thermodynamic quantities can be corrected using the previous estimations, but this has to

Although the LJ potential is quite simple, some situations require to use more sophisticated semi-empirical potentials whose parameters are chosen to reproduce either microscopic or macroscopic properties of a model system. These semi-empirical potentials are usually no longer pair potentials like the LJ potential but many-body, where the many-body terms

As an illustration, let us mention the embedded atom model (EAM) used in the context of

*ij <sup>i</sup> i j <sup>i</sup>*

the local density measured using the positions of neighbouring particles *i i <sup>j</sup> <sup>j</sup>*

microscopic and macroscopic properties to be reproduced, usually the lattice constant, the

Another limitation of LJ potentials is the description of solids crystalizing in non-compact structures. Indeed, Lennard Jones atoms form close-packed structures at low temperatures (usually fcc at low pressure), and thus LJ potentials are also not adapted to model materials which crystallize in opened structures such as diamond. Semi-conductors like silicon or germanium are usually modelled using many-body semi-empirical potentials, in which the many body terms account for the local preferential bond ordering of the semi-conductor. Among the most famous potentials for silicon, let us quote the Tersoff potential (Tersoff 1986, 1988a, 1988b) and the Stillinger-Weber potential (Stillinger & Weber 1985). The parameters of the former are chosen so as to reproduce the elastic properties of bulk silicon, while the latter describes satisfactorily the structural properties of liquid silicon. Of course the list of potentials is non exhaustive and there exists plethora of empirical potentials to model as disparate systems as charged systems, liquid crystals, polymers, surfactants,

 

which accounts for the cohesion of the metallic bond and which depends on

*r+ u r* (10)

and the weight function *w r* depends on the

*r = wr* .

*r* accounting for repulsion at short distances and a many body

describes how does the potential energy of an atom depends on its coordination.

elastic constants and the sublimation energy of the metal (Daw & Baskes 1984).

*cr U = r V r dr*

*<sup>V</sup> P = r r dr r*

2

be kept in mind when calculating a phase diagram for instance.

metals (Daw & Baskes 1984), where the potential

consists of a pair term

The choice of the functional *u*

granular media,... and mixtures of them!

term *u r* 

<sup>6</sup> *cr*

and

 2 3 <sup>2</sup> 4 1,07

<sup>2</sup> <sup>2</sup> 2 3 4 1.07

 

(8)

(9)

$$KE = \left\langle \sum\_{i} \frac{1}{2} m \vec{v}\_i^2 \right\rangle \tag{12}$$

where the brackets denote an ensemble average. If sufficient time has been left to the system to reach equilibrium, the ensemble average may be performed by averaging over a sufficient long time thanks to ergodicity. On the other hand, in non-equilibrium situations, different initial conditions (usually the positions and the velocities of the atoms) may be used to generate independent trajectories in phase space. The latter definition of the kinetic energy is intimately related to the definition of the local temperature :

$$T\left(\vec{r}\right) = \left\langle \sum\_{i \in V(\vec{r})} \frac{m\vec{v}\_i^2}{\Im m(\vec{r})k\_B} \right\rangle \tag{13}$$

where the sum runs over particles i in a small volume V(r) centred around r and which contains n(r) particles. Another important quantity in heat transport simulations is the energy flux, defined as

$$\vec{J} = \frac{d}{dt} \left(\sum\_{i} E\_i \vec{\tau}\_i\right) \tag{14}$$

where *Ei* is the total energy of particle i, yielding for pair potentials:

$$\vec{J} = \sum\_{i} E\_{i}\vec{v}\_{i} + \frac{1}{2} \sum\_{i$$

For simulations in crystals at low temperatures, the previous energy flux displays large oscillations due to optical phonons which carry a negligible amount of heat, and sometimes it can be more suited to work with the equivalent definition of the flux:

$$\vec{J} = \frac{1}{2} \sum\_{i$$

where the superscript "0" denotes the equilibrium positions of the atoms. From the knowledge of the local temperature and the energy flux, it is possible to measure thermal conductivity and we would come back to the measurement of transport coefficients in MD in the next section, with a particular emphasis on heat transport.

The other type of information that one can extract from MD is structural. Since the output of a MD simulation is the set of positions *ij r* , it is possible to compute the pair distribution function

$$\rho \,\rho \mathbf{g}\left(\vec{r}\right) = \left\langle \frac{1}{N} \sum\_{i \neq j} \mathcal{S}\left(\vec{r} - \vec{r}\_{ij}\right) \right\rangle \tag{17}$$

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 79

 , , 0 *k k S tS* 

*k*

the grand-canonical ensemble for instance (Frenkel & Smit, 1996).

It is a straightforward exercize to show that

/ and group velocities *v k= <sup>g</sup> <sup>k</sup>*

Secondly, from the time decay of *S tS k, k,*

relation

mode *k,*

 

(Ladd 1986):

*v k= k*

 

 

 2 2 , *S t kT w k k*

 , , , 2 ,

 

*S t S dt*

*= k* for the two polarizations. From the dispersion relation, the phonon phase

may be predicted.

0 0

Hence, from the previous autocorrelation function, one can determine the dispersion

> 

*k k*

 

*S t* 

The wavelength dependance of the phonon lifetime may then be compared to the theoretical predictions, given either by Callaway 1959 or Holland 1963, see e.g. (McGaughey 2004).

As we have already mentioned, the conversion of the microscopic information given by a MD simulation to macroscopic observables requires to evaluate averages over phase space. The phase space is a multidimensional space generated by the positions and the momenta of a classical system, and has dimension of 6N for a system of N particles. A molecular dynamics simulation generates a sequence of points in phase space as a function of time, whose points belong to the same thermodynamic ensemble. The first MD simulations were performed in the micro-canonical ensemble or NVE ensemble, which corresponds to a fixed number of atoms N, a fixed volume V, and a fixed energy E. Although these simulations are very easy to implement, they do not represent a realistic system as the set of molecules studied is completely isolated. The need to model real situations led to the development of thermostats for MD, which allow to simulate the dynamics of a system in the canonical or

Although MD can be used to calculate the thermodynamics, structural and vibrational properties of systems at equilibrium, it is more designed to study the non-equilibrium situations where the system under study is driven by an external force; the latter can be of thermodynamic nature. In the language of out of equilibrium statistical physics, MD can help in determining the relation between the forces and the fluxes. Statistical physics predicts that there is proportionality between the forces and the fluxes if the system under study is driven weakly out of equilibrium. The coefficients of proportionality are called transport coefficients and measure the susceptibility of the system to respond to a given thermodynamic force. For instance, the thermal conductivity quantifies the amount of energy flowing in a material submitted to a temperature gradient. MD is a very powerful tool to calculate the transport coefficients of a model system. To this end, two routes can be traditionally followed: either the model system is submitted to an external force or the thermally induced fluctuations of an internal variable are probed at equilibrium. Statistical physics states indeed that the typical time decay of the spontaneous thermal fluctuations is

*k*

ű ű (23)

0 / *<sup>B</sup>* ű ű (24)

0 , one can compute the lifetime of the

ű ű (25)

where is the mean number density of the system. Physically, g(r) represents the probability to find two particles separated by a distance r, and <sup>2</sup> 4*r g r dr* is the average number of particles located at a distance between r and r+dr from a given particle. Note that sometimes, it may be more useful to compute the structure factor:

$$S(\vec{q}) = \left\langle \frac{1}{N} \sum\_{i,j} \exp\left(i\vec{q} \cdot \left(\vec{r}\_i - \vec{r}\_j\right)\right) \right\rangle \tag{18}$$

which is simply related to the pair distribution function through a Fourier transform:

$$S(\vec{q}) = 1 + \rho \int \mathcal{g}(\vec{r}) \exp(-i\vec{q} \cdot \vec{r}) d\vec{r} \tag{19}$$

Finally, MD may be used to compute the vibrational properties of a model system. In particular, the vibrational density of states (DOS)

$$\log\left(\phi\right) = \frac{1}{V} \sum\_{k} \delta\left(\phi - \phi\_{k}\right) \tag{20}$$

where the sum runs over the eigenmodes of the system can be estimated using the Fourier transform of the velocity autocorrelation function:

$$\lg\left(\rho\right) = \bigwedge\_{0}^{\circ} \langle \vec{v}\big(t\big) \cdot \vec{v}\big(0\big) \rangle \exp\big(-i\alpha t\big) dt \tag{21}$$

where here *v t* is the velocity of an atom. One has to keep in mind that the previous equation is only approximate, as it assumes small deviations from harmonicity, and thus is only valid strictly speaking at low temperatures. Note also that the resulting DOS mixes information about the different polarizations of the atomic vibrations, e.g. transversal or longitudinal, and to determine the DOS relative to a particular polarization, a more refined analysis of the atomic displacement should be undertaken. Apart from the DOS, MD is a common technique to characterise the propagation of phonons in a crystal.

The vibrations of a crystal containing N atoms may be decomposed in normal modes:

$$S\_{\vec{k},\nu}\left(t\right) = N^{-1/2} \sum\_{i} m\_i^{1/2} \exp\left(-i\vec{k}\cdot\vec{r}\_i^0\right) \vec{e}\left(\vec{k},\nu\right) \cdot \vec{u}\_i \tag{22}$$

where *u t <sup>i</sup>* denotes the displacement of atom i in the crystal vibrating around the equilibrium position <sup>0</sup> *ir* , is an index designating the polarization of the vibration, e.g. longitudinal or transversal and *e k,* is the corresponding polarization vector. From the normal mode amplitude, one can compute the dispersion relation and the phonon lifetime as follows. First consider the autocorrelation of the atomic displacement:

$$
\langle \mathfrak{dS}\_{\vec{k},\nu}(t) \overline{\mathfrak{S}}\_{\vec{k},\nu}(0) \mathfrak{d} \rangle \tag{23}
$$

It is a straightforward exercize to show that

78 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

The other type of information that one can extract from MD is structural. Since the output of

 <sup>1</sup> *ij i j gr = r r <sup>N</sup>*

number of particles located at a distance between r and r+dr from a given particle. Note that

 

is the mean number density of the system. Physically, g(r) represents the

sometimes, it may be more useful to compute the structure factor:

,

particular, the vibrational density of states (DOS)

transform of the velocity autocorrelation function:

*ir* , 

longitudinal or transversal and *e k,*

equilibrium position <sup>0</sup>

probability to find two particles separated by a distance r, and <sup>2</sup> 4

1

which is simply related to the pair distribution function through a Fourier transform:

*S q = + g r iq r dr* 1 exp 

Finally, MD may be used to compute the vibrational properties of a model system. In

<sup>1</sup>

*k*

where the sum runs over the eigenmodes of the system can be estimated using the Fourier

*g v t v exp i t dt*

where here *v t* is the velocity of an atom. One has to keep in mind that the previous equation is only approximate, as it assumes small deviations from harmonicity, and thus is only valid strictly speaking at low temperatures. Note also that the resulting DOS mixes information about the different polarizations of the atomic vibrations, e.g. transversal or longitudinal, and to determine the DOS relative to a particular polarization, a more refined analysis of the atomic displacement should be undertaken. Apart from the DOS, MD is a

 

*<sup>g</sup> <sup>V</sup>* 

0

common technique to characterise the propagation of phonons in a crystal.

as follows. First consider the autocorrelation of the atomic displacement:

The vibrations of a crystal containing N atoms may be decomposed in normal modes:

*i*

 1/2 1/2 <sup>0</sup> , , *<sup>k</sup> <sup>i</sup> i i*

where *u t <sup>i</sup>* denotes the displacement of atom i in the crystal vibrating around the

normal mode amplitude, one can compute the dispersion relation and the phonon lifetime

*S t N m exp ik r e k u*

*k*

 

0

, it is possible to compute the pair distribution

(17)

exp *i j i j Sq = iq r r <sup>N</sup>* (18)

(19)

(21)

(22)

is the corresponding polarization vector. From the

is an index designating the polarization of the vibration, e.g.

*r g r dr* 

(20)

is the average

a MD simulation is the set of positions *ij r*

function

where

$$
\langle \text{d}\mathbb{S}\_{\vec{k},\nu} \left( t = 0 \right) \vec{\text{d}} \rangle = k\_B T \;/\; w\_{\vec{k}}^2 \tag{24}
$$

Hence, from the previous autocorrelation function, one can determine the dispersion relation *= k* for the two polarizations. From the dispersion relation, the phonon phase *v k= k* / and group velocities *v k= <sup>g</sup> <sup>k</sup>* may be predicted.

Secondly, from the time decay of *S tS k, k,* 0 , one can compute the lifetime of the mode *k,* (Ladd 1986):

$$\tau\_{\vec{k}\_{,\nu}} = \frac{\int \langle S\_{\vec{k},\nu} \left( t \right) \overline{S}\_{\vec{k},\nu} \left( 0 \right) dt}{\langle \text{US}\_{\vec{k},\nu} \left( t = 0 \right) \hat{\mathbf{f}} \rangle} \tag{25}$$

The wavelength dependance of the phonon lifetime may then be compared to the theoretical predictions, given either by Callaway 1959 or Holland 1963, see e.g. (McGaughey 2004).

As we have already mentioned, the conversion of the microscopic information given by a MD simulation to macroscopic observables requires to evaluate averages over phase space. The phase space is a multidimensional space generated by the positions and the momenta of a classical system, and has dimension of 6N for a system of N particles. A molecular dynamics simulation generates a sequence of points in phase space as a function of time, whose points belong to the same thermodynamic ensemble. The first MD simulations were performed in the micro-canonical ensemble or NVE ensemble, which corresponds to a fixed number of atoms N, a fixed volume V, and a fixed energy E. Although these simulations are very easy to implement, they do not represent a realistic system as the set of molecules studied is completely isolated. The need to model real situations led to the development of thermostats for MD, which allow to simulate the dynamics of a system in the canonical or the grand-canonical ensemble for instance (Frenkel & Smit, 1996).

Although MD can be used to calculate the thermodynamics, structural and vibrational properties of systems at equilibrium, it is more designed to study the non-equilibrium situations where the system under study is driven by an external force; the latter can be of thermodynamic nature. In the language of out of equilibrium statistical physics, MD can help in determining the relation between the forces and the fluxes. Statistical physics predicts that there is proportionality between the forces and the fluxes if the system under study is driven weakly out of equilibrium. The coefficients of proportionality are called transport coefficients and measure the susceptibility of the system to respond to a given thermodynamic force. For instance, the thermal conductivity quantifies the amount of energy flowing in a material submitted to a temperature gradient. MD is a very powerful tool to calculate the transport coefficients of a model system. To this end, two routes can be traditionally followed: either the model system is submitted to an external force or the thermally induced fluctuations of an internal variable are probed at equilibrium. Statistical physics states indeed that the typical time decay of the spontaneous thermal fluctuations is

Molecular Dynamics Simulations and Thermal Transport at the Nano-Scale 81

semiconductors (Yang, 2004). There are three principal techniques used to evaluate the thermal conductivity with molecular dynamic simulations (Allen & Tildesley 1987, Chantrenne, 2007); a. the equilibrium approach based on the Green-Kubo formulae (Frenkel & Smit, 1996), b. the non-equilibrium MD, also called the direct method, which uses a heat source and a heat sink, or the so-called direct method, based on the creation of a temperature gradient, and c. the homogeneous non-equilibrium MD, where a heat flux is induced (Evans, 1982, 1986). The comparison between the two methods has been undertaken by many researchers, who have concluded that the two methods give consistent results (Shelling et al, 2002; Mahajan et al, 2007; Termentzidis et

Kotake and Wakuri proposed the direct method (Kotake & Wakuri, 1994), which is similar to the hot-plate experiment setup. A temperature gradient is imposed across the structure under study by allowing thermal power exchange between the heat source and sink and measure the resulting heat flux (Chantrenne & Barrat, 2004a, 2004b, Termentzidis et al, 2009). The thermal conductivity is then obtained as the ratio of the heat flux and the temperature gradient. An alternative, but equivalent way consists in inducing a heat flux and to measure the resulting temperature gradient (Muller-Plathe, 1997). In both cases the system is first allowed to reach a steady state, after which prolong simulations are conducted allowing to obtain correct statistical measurements (Stevens et al, 2007). The NEMD method is often the method of choice for studies of nanomaterials (Poetzsch & Botter, 1994) while for bulk thermal conductivity, particularly of high-conductivity materials, equilibrium method is typically preferred due to less severe size effects. The size effects in the determination of the thermal conductivity using NEMD can be understood using the following simple line of thought (Schelling & al 2002). Imagine a phonon travelling in the crystal between the two heat reservoirs. This phonon can experience at least two kinds of scattering events: either it can be scattered by other phonons travelling in the material, or it can be scattered by the reservoirs which are seen by the considered phonon as a different material with an almost infinite thermal conductivity. The mean free path associated with this former mechanism can be roughly approximated by / 2 *reservoir = L* , where L is the distance between the reservoirs, because on average a phonon will travel ballistically a distance L/2 before being scattered by the heat source/sink. Now invoking Mathiessen rule to predict the effect of reservoir scattering on the overall thermal conductivity, the effective mean free path is

where the last term *L* describes the phonon-phonon interaction in a bulklike medium. The length dependence of the thermal conductivity can now be estimated

volume, *v* is the mean group velocity. Of course, this analysis is simplified because we have assumed that the effect of phonon-phonon scattering can be described in terms of a

1/ 2/ 1/ *L = L+ L* (26)

*L = cv L* /3 where *c* is the heat capacity per unit

**3.1 Non-equilibrium molecular dynamics method (NEMD)** 

al, 2011b).

given by:

using kinetic theory of gas:

proportional to the appropriate transport coefficient. This will be illustrated in the context of the heat transfer in the next sextion. Before focusing on energy transfer applications of MD, let us mention that MD can be used also beyond the linear response domain of out of equilibrium systems. This includes situations of large external forcing such as e.g. polymer melts flowing at large shear rates (Vladkov 2006b). But this includes also systems for which intrinsically the linear relationship between forces and fluxes is violated. As an example consider a nano-structured system at low temperatures for which the phonon mean free path is comparable with the characteristic system size. Heat transfer in such systems is no longer described by Fourier law, but is rather described by an effective conductivity which depends on the strength of the thermal flux flowing across the system. Molecular dynamics has been shown to predict the effective conductivity under conditions where the heat carriers travel ballistically in the system and being scattered only by the boundaries of the nanostructure.
