How to Improve The Accuracy of
Equilibrium Molecular Dynamics
For Computation of
Thermal Conductivity
CHEN Jie
Dept. of Physics & Centre for Computational
Science and Engineering
Department of Physics & CCSE
Outline
Part I: Equilibrium molecular dynamics:
how to improve accuracy
Part II: Non-equilibrium molecular dynamics:
effect of heat bath
Department of Physics & CCSE
Introduction
 Peierls-Boltzmann Phonon Gas Theory
καβ = ∑ Cλ vλα vλβτ λβ
λ

λ = ( q, j )
 Monte Carlo (MC) simulations
Metropolis probability p:
 Molecular Dynamics (MD) simulations
Atomistic simulation with empirical potential
Department of Physics & CCSE
Equilibrium MD Simulations
Heat current auto-correlation function (HCACF)
1
=
Green-Kubo
formula κ
3k BT 2V
∫
∞
0


dt J (0) • J (t )
Department of Physics & CCSE
Review of Different Approaches
1. Direct integration
Bulk Si at 1000 K
Time step: 0.55 fs
Total time: 1.1 ns
P. Schelling et al., Phys. Rev. B 65, 144306 (2002)
Department of Physics & CCSE
Thermal Conductivity of Si Crystal (Expt.)
Expt. value: ~32 W/m-K (at 1000 K)
C. Glassbrenner et al., Phys. Rev. 134, 1058 (1964)
Y. Magomedov et al., High Temperature 46, 422 (2008)
Department of Physics & CCSE
Review of Different Approaches
2. Exponential fitting of HCACF
Temporal decay of HCACF for
perfect β-SiC at 760 K. The dashed
line is the fit to exponential decay in
the range from 10 to 30 ps.
Fit HCACF according to Ae
− t /τ 0
in time domain [τ 1 ,τ 2 ]
J. Li et al., J. Nucl. Mater. 255, 139 (1998)
Department of Physics & CCSE
Review of Different Approaches
2. Exponential fitting of HCACF
Give an underestimated value of 38 W/m-K compared to direct integration
P. Schelling et al., Phys. Rev. B 65, 144306 (2002)
Department of Physics & CCSE
Review of Different Approaches
3. Double exponential fitting of HCACF
Expt. value: 2300 W/m-K (300K)
J. Che et al., J. Chem. Phys 113, 6880 (2000)
Department of Physics & CCSE
Thermal Conductivity Decomposition
Lennard-Jones argon
A. McGaughey et al., Int. J. Heat Mass Transfer 47, 1783 (2004)
Department of Physics & CCSE
Simulation Details
 Stillinger-Weber (SW) potential, bulk Si and Ge
 Canonical ensemble MD for thermal equilibrium
 Langevin heat bath at 1000 K
 Micro-canonical ensemble MD to record heat current

1
1

  
   
 Heat current J (=t ) ∑ v ε + 2 ∑ r ( F ⋅ v ) + 6 ∑ (r + r )( F ⋅ v )
i
i
i
ij i ≠ j
ij
ij
i
ijk i ≠ j
j ≠k
ij
ik
ijk
i
 Total simulation time 1.6 ns with time step 0.8 fs
 Average over 8 runs with different initial conditions
Department of Physics & CCSE
Determine Cut-off Time
Time dependence of HCACF
for 2 different realizations in
a 4*4*4 super cell (0-400ps).
Indicator
F (t ) =
σ (Cor (t ))
E (Cor (t ))
Accumulative conductivity
J. Chen et al., Phys. Lett. A
375, 2392 (2010)
Department of Physics & CCSE
Validity of Direct Integration
Thermal conductivity of the
crystalline Silicon (a) and
Germanium (Ge) at 1000 K
versus super cell size
(N*N*N unit cells) using
direct integration method.
J. Chen et al., Phys. Lett. A 375, 2392 (2010)
Department of Physics & CCSE
Correction to Double Exponential Fitting
Fitting formula
Cor (t ) / Cor (0) = A1e − t /τ1 + A2 e − t /τ 2 + Y0
Thermal conductivity
=
κ C (t )
Cor (0)
( A1τ 1 + A2τ 2 )
2
3k BT V
Cor (0)
( A1τ 1 + A2τ 2 + Y0τ c )
2
3k BT V
J. Chen et al., Phys. Lett. A 375, 2392 (2010)
=
κ F (t )
Department of Physics & CCSE
Validity of Double Exponential Fitting
J. Chen et al., Phys. Lett. A 375, 2392 (2010)
Department of Physics & CCSE
Conclusion for Part I
1.
Due to the finite number of ensemble average, HCACF is only
reliable up to a cut-off time and thus thermal conductivity can
only be calculated from the truncated HCACF.
2. We have proposed an efficient quantitative method to accurately
estimate the cut-off time. Using the cut-off time, direct
integration method can make a successful computation of
thermal conductivity of crystalline silicon and germanium.
3. Because of the finite cut-off time, a small nonzero correction
term has to be included to improve the accuracy of EMD
calculations based on the double-exponential-fitting of HCACF.
Department of Physics & CCSE
Non-equilibrium MD Simulations
Temperature profile of (100) SiNWs
Schematic setup
Fourier’s law
κ = − J /∇T
Department of Physics & CCSE
Canonical Ensemble (NVT)
 Stochastic heat bath
e.g. Langvin heat bath
Fluctuation-dissipation
ξ ~ N (0, 2mλ k BT )
 Deterministic heat bath
e.g. Nosé-Hoover heat bath
Dynamics of ζ
Department of Physics & CCSE
Simulation Details
 Stillinger-Weber (SW) potential
 3*3*10 SiNWs along (100) direction
 Heat bath temperature: T =310K, T =290K
 Total simulation time 40 ns with time step 0.8 fs
H
L
Department of Physics & CCSE
No. of Heat Bath Layers
(a)-(c) Temperature profile
(d) Heat flux
(e) Thermal conductivity
(f) Normalized correlation fuction of
velocity for atoms in heat bath
Berendsen: velocity scaling heat bath
J. Chen et al., J. Phys. Soc. Jpn. 79, 074604 (2010)
Department of Physics & CCSE
Effect of Heat Bath Parameter
(a)-(b) Temperature profile
(c) Heat flux
(d) Thermal conductivity
Fix NL=3
J. Chen et al., J. Phys. Soc. Jpn. 79, 074604 (2010)
Department of Physics & CCSE
Si-Ge Junction
(a)-(b) Heat flux
(c) Rectification Ratio
J+/J-: Si/Ge attached to high T
Temperature profile is correct in all cases
J. Chen et al., J. Phys. Soc. Jpn. 79, 074604 (2010)
Department of Physics & CCSE
Experiment Result
Larger heat flow in the direction of
decreasing mass density
C. W. Chang et al., Science 314, 1121 (2006)
Department of Physics & CCSE
Conclusion for Part II
1. Due to the existence of localized edge modes and its accumulation
effect because of the deterministic characteristic of Nosé-Hoover
heat bath, multiple heat bath layers are required in order to reduce
the temperature jump at the boundary, while one layer of Langevin
heat bath is good enough to reduce temperature jump due to its
stochastic excitation of all modes.
2. In order to obtain the correct temperature profile, intermediate
values of heat bath parameters are recommended for both NoséHoover and Langevin heat bath.
3. Deterministic heat bath can give rise to unphysical results.
Department of Physics & CCSE
Acknowledgement
 National University of Singapore:
Prof. Li Baowen
 Peking University:
Prof. Zhang Gang
 National University of Singapore:
Prof. Wang Jian-Sheng