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

© Copyright 2020