Article pubs.acs.org/JPCC
Anisotropic Thermal Transport in Organic−Inorganic Hybrid Crystal β‑ZnTe(en)0.5 Xin Qian, Xiaokun Gu, and Ronggui Yang* Department of Mechanical Engineering & Materials Science and Engineering Program, University of Colorado, Boulder, Colorado 80309-0427, United States ABSTRACT: By using the interatomic potential derived from ab initio simulations, equilibrium molecular dynamics simulations using the Green−Kubo relation were carried out to study the elastic constants and thermal conductivity of a layered organic−inorganic hybrid crystal β-ZnTe(en)0.5, whose inorganic ZnTe monolayers are connected by organic ligands ethylenediamine (en) with covalent bonds. As compared to inorganic ZnTe, the results of elastic constants showed that βZnTe(en)0.5 is much more flexible, especially in terms of the shearing stiffness. Low thermal conductivity values are found in β-ZnTe(en)0.5. At 300 K, the thermal conductivities are kST = 1.2 W/m·K in the stacking direction normal to the ZnTe monolayers, kAM = 0.8 W/m·K in the armchair direction, and kZZ = 1.8 W/m·K in the zigzag direction parallel to the ZnTe monolayers, respectively. The low thermal conductivity across the ZnTe monolayers is determined by the mismatch of phonon spectra between the organic ligands and ZnTe monolayers. The anisotropic thermal conductivity parallel to the ZnTe monolayers is dominated by the different interchain coupling strengths in the armchair and zigzag directions.
1. INTRODUCTION There has been growing interest in organic−inorganic hybrid crystals, which can combine the superior properties of both components, such as the good electronic properties from inorganic materials and superb flexibility from organic materials.1 The properties of organic−inorganic hybrid crystals could be easily tuned through using different organic ligands and changing the dimensionality of the lattice.2,3 As compared to organic−inorganic nanocomposites with a characteristic structure size in the order of 1−100 nm, the blending of organic and inorganic components in hybrid crystals happens at an atomic scale by forming strong chemical bonds.1 Because of the unique structural features, the electronic and optical properties of organic−inorganic hybrid crystals have been investigated extensively in the past two decades for their potential applications in electronics, photonics, and energy conversion devices.2−6 However, the thermal conductivity of organic−inorganic hybrid crystals is rarely studied, although understanding phonon transport in these novel materials is essential for functionality and reliability of electronic and optical devices using hybrid crystals. Thermal transport in hybrid crystals could exhibit quite different behavior from inorganic crystals and other hybrid composites. Previous research on thermal conductivity of organic− inorganic hybrid materials focused on either porous metal organic frameworks (MOFs)7,8 or hybrid nanocomposites with inorganic−organic interfaces, such as carbon nanotube− polymer composites,9,10 self-assembled materials,11−14 nanocrystal arrays,15−17 and superlattices.18−20 The dominant mechanism for thermal transport in these hybrid composites is usually phonon-boundary scattering, while recent exper© 2015 American Chemical Society
imental work on hybrid crystals suggests that the organic ligands can intrinsically affect the thermal transport properties.21 Although much progress has been made on the computational modeling on the thermal conductivity of nanostructured materials, the knowledge on the phonon transport in organic− inorganic hybrid crystals is still lacking, due to the complex crystal structure and the difficulty in describing interactions between organic and inorganic components. There are several methods that could potentially be used to predict thermal conductivity of hybrid crystals, including the solution of Boltzmann transport equation (BTE) and the molecular dynamics (MD) simulations. The BTE approach calculates the thermal conductivity using phonon properties including phonon dispersion relation and phonon−phonon scattering time. Although phonon dispersion can be calculated from the first-principles calculations, the computation time needed for the detailed phonon scattering rate22,23 increases dramatically with the number of atoms in a unit cell. On the other hand, the MD simulations suffer from the requirement of an appropriate potential field to accurately describe the interatomic forces. In this Article, we perform MD simulations to calculate the elastic constants, phonon properties, and thermal conductivity of layered organic−inorganic hybrid crystal β-ZnTe(en)0.5 by developing a potential field based on the first-principles calculations. The hybrid crystal β-ZnTe(en)0.5 with alternating ZnTe monolayers and ethylenediamine (en) ligands connected Received: September 30, 2015 Revised: November 22, 2015 Published: November 30, 2015 28300
DOI: 10.1021/acs.jpcc.5b09527 J. Phys. Chem. C 2015, 119, 28300−28308
Article
The Journal of Physical Chemistry C
Figure 1. Crystal structure of β-ZnTe(en)0.5. Two ZnTe monolayers are connected by ethylenediamine chains. Hydrogen atoms are omitted in the figure for clarity. Elastic constants at diagonal of stiffness matrix and the corresponding deformation of lattice are indicated by arrows. (a) Definition of unit cell is outlined in the figure. (b) Structure of ZnTe monolayer. (c) Lattice structure projected on AM−ST plane, with ZZ axis pointing out of the page. (d) Lattice structure projected on ST−ZZ plane, with AM axis pointing out of the page. Distances between carbon atoms of neighboring organic chains are indicated in (c) and (d).
cell, as shown in the rectangular box in Figure 1a, contains 32 atoms in total, including 4 Zn atoms, 4 Te atoms, 4 C atoms, 4 N atoms, and 16 H atoms. The ST axis is defined along the stacking direction normal to the ZnTe monolayers, and the AM and ZZ axes are defined along the armchair and the zigzag directions parallel to the ZnTe monolayer shown in Figure 1b, respectively. Figure 1c and d shows the side-views of the lattice in the ST−AM plane and in the ST−ZZ plane, respectively. This section presents a general simulation scheme we developed for the prediction of the thermal conductivity of organic−inorganic hybrid crystals. A potential field capable of accurately describing the lattice vibration is developed first using the first-principles-based DFT calculations and is used as inputs for the equilibrium molecular dynamics simulations. We describe the details for potential development in part A, followed by the calculation of elastic constants in part B, molecular dynamics simulation of thermal conductivity in part C, and the spectral analysis of phonon properties in part D. A. Potential Development. There are some models available for potential fields of both the inorganic ZnTe crystal and the organic (en) molecules. For example, Tersoff potential,24 Stillinger−Weber potential,25 Rhaman−Vashishta potential,26 and bond-order potential27 have been developed to describe the zinc-blende ZnTe crystal, and the organic component (en) molecules can be described by the COMPASS class II force field.28 However, these potential fields cannot be directly applied to organic−inorganic hybrid crystal β-ZnTe(en)0.5 for two reasons: (1) The structure of β-ZnTe(en)0.5 is very different from that of its inorganic and organic counterparts. For example, each Te atom in β-ZnTe(en)0.5 only forms three Zn−Te covalent bonds, but there are four Zn−Te bonds in a ZnTe crystal; each N atom forms an additional coordinative N−Zn bond in β-ZnTe(en)0.5, which
by covalent bonds is a member of the II−VI compound-based hybrid crystal family, which is finding promising applications in white light emitting diodes (LEDs).2 This work shows that the hybrid crystal β-ZnTe(en)0.5 is much more flexible than the inorganic crystal ZnTe, and the predicted thermal conductivity is one order-of-magnitude lower than that of inorganic ZnTe. We also found that thermal transport in β-ZnTe(en)0.5 is anisotropic while different mechanisms are responsible for phonon transport in different directions. Parallel to the ZnTe monolayers, phonon transport is determined by the difference of interchain coupling strengths in the armchair and zigzag directions, whereas the phonon transport in the direction across the ZnTe monolayers is dominated by the phonon spectrum mismatch between the ZnTe monolayers and (en) chains. This Article is organized as follows. In section 2, we describe the detailed atomistic simulation methods used for potential development, the calculation of elastic constants, the simulation for thermal conductivity, and the calculation of phonon properties. In section 3, we present the results and analysis, followed by a summary in section 4. The contribution of this work to the literature is two-fold: (1) developed a general modeling strategy that integrates first-principles-based density functional theory (DFT) calculations with molecular dynamics simulations for thermal conductivity prediction of complex organic−inorganic hybrid crystals; and (2) reported the low and anisotropic thermal conductivity values in organic− inorganic hybrid crystal ZnTe(en)0.5 and identified the origins of the anisotropicity.
2. SIMULATION METHODS Figure 1a shows the crystal structure of an organic−inorganic hybrid crystal β-ZnTe(en)0.5 whose inorganic ZnTe monolayers are connected by ethylenediamine (en) chains. The unit 28301
DOI: 10.1021/acs.jpcc.5b09527 J. Phys. Chem. C 2015, 119, 28300−28308
Article
The Journal of Physical Chemistry C Table 1. Interatomic Potentials of β-ZnTe(en)0.5 interactiona
potential form
Nearest Neighbor Zni−Tei Morse potential: E = D[1 − e−α(r−r0)]2 Zni−Tej Zn−N N−C C−Hc, N−Hnb COMPASS Class II Angular Interaction Tei−Zni−Tej cosine squared: E = K (cos θ − cos θ0)2 Tej−Zni−Tej Tei−Zni−N Tei−Znj−N Zni−Tei−Znj Znj−Tei−Znj Zn−N−C Zn−N−Hn N−C−C N−C−Hc C−N−Hn COMPASS Class II Hc−C−Hc Hn−N−Hn Second Nearest Neighbor Zni−Zni Morse potential: E = D[1 − e−α(r−r0)]2 Zni−Znj Zn−C Tei−Tei Tei−Tej
interactiona
parameters D = 49.4 kcal/mol, α Å−1, r0 = 2.668 Å D = 50.2 kcal/mol, α Å−1, r0 = 2.693 Å D = 57.3 kcal/mol, α Å−1, r0 = 2.138 Å D = 70.0 kcal/mol, α Å−1, r0 = 1.474 Å ref 28
potential form
parameters
Second Nearest Neighbor Te−N
= 1.184
= 1.673
LJ Couplingc N−N N−C
= 2.730
N−Hn
= 1.197
K = 34.0 kcal/mol, θ0 = 115.8° K = 34.0 kcal/mol, θ0 = 109.9° K = 28.8 kcal/mol, θ0 = 104.7° K = 34.5 kcal/mol, θ0 = 104.5° K = 28.8 kcal/mol, θ0 = 99.1° K = 34.0 kcal/mol, θ0 = 109.9° K = 33.4 kcal/mol, θ0 = 117.2° K = 6.21 kcal/mol, θ0 = 106.6°
r < rin:
N−Hc
E = 4ε[(σ/r)12 − (σ/ r)6] rin < r < rc:
C−C
F = ∑4n=1 Cn(r − rin)n−1
C−Hn C−Hc
rin = 6.0 Å
Hn−Hn
rout = 8.0 Å
Hn−Hc Hc−Hc
ref 28
Torsion Hn−N−C−Hc Hc−N−C−C N−C−C−Hc N−C−C−N
D = 16.75 kcal/mol, α = 0.645 Å−1, r0 = 4.407 Å D = 16.88 kcal/mol, α = 0.658 Å−1, r0 = 4.4079 Å D = 20.81 kcal/mol, α = 1.683 Å−1, r0 = 3.103 Å D = 15.73 kcal/mol, α = 0.631 Å−1, r0 = 4.407 Å D = 16.31 kcal/mol, α = 0.627 Å−1, r0 = 4.407 Å
D = 12.04 kcal/mol, α = 0.796 Å−1, r0 = 3.814 Å
E = K [1 + cos(nϕ)]
ε = 0.69 kcal/mol, σ = 3.412 Å ε = 0.085 kcal/mol, σ = 2.239 Å ε = 0.055 kcal/mol, σ = 1.892 Å ε = 0.055 kcal/mol, σ = 1.867 Å ε = 0.105 kcal/mol, σ = 4.262 Å ε = 0.0685 kcal/mol, σ = 1.845 Å ε = 0.068 kcal/mol, σ = 1.931 Å ε = 0.044 kcal/mol, σ = 1.461 Å ε = 0.044 kcal/mol, σ = 2.002 Å ε = 0.044 kcal/mol, σ = 1.585 Å K K K K
= = = =
0.1581 0.1581 0.1581 0.1581
kcal/mol, kcal/mol, kcal/mol, kcal/mol,
n n n n
= = = =
3 3 3 1
i,j = 1,2 and i ≠ j. bHn and Hc are hydrogen atoms bonded to N atoms and C atoms, respectively. N corresponds to “n3” type in ref 28. c The coefficients Cn of smoothing force are calculated by LAMMPS package. The LJ coupling for the bonded atoms is switched off. For the coupling between the nonadjacent atoms in the same (en) chain, the LJ interactions are switched on. For example, in the angular Hn−N−C, only the LJ coupling between the Hn and C atom is considered. a
θ0)2, where θ0 is the angle at the equilibrium position. Contributions to the ground-state energy from the interaction beyond third-nearest neighbors are integrated into the secondnearest neighbor interaction and the angular terms, because including them explicitly in the potential field would introduce too many atom types. For the organic (en) ligands, the Lennard-Jones (LJ) potential was used to describe the nonbonded interactions including the coupling between different (en) ligands (interchain coupling) and the nonbonded interactions inside the (en) chain (intrachain coupling). For the bonded interactions in (en) ligands, the LJ interaction is switched off, and we directly use the functional forms and parameters from the COMPASS force field28 for the bond stretching (C−C, C−H, and N−H) and angular interactions (N−C−C, H−C−H, H−C−N, and H−N−C). Parameters for bond-stretching terms of Zn−N and N−C and the three-body angular interactions of Zn−N−C and Te−Zn−N are also optimized to fit the interaction energy between ZnTe monolayers and organic ligands. We optimized the parameters for the N−C bonds rather than directly using parameters from the COMPASS28 force field because the state of N atom is changed by the additional Zn−N coordinative bond. Finally, we use a cosine form E = K(1 + cos nϕ) to describe the torsional energy in the (en) ligands, where K and n are the stiffness and angle of the dihedral, and n is an integer parameter. The dihedral torsional terms outside the (en) ligands are not
does not exist in (en) molecules. (2) The existing classical potential fields do not have the capability to describe the covalent bonding between the inorganic ZnTe monolayers and organic (en) chains. The first-principles-based DFT calculations have been used in this work to develop the potential field. This method has been used to construct potential fields for both inorganic and organic materials, even for those with complex unit cells.7,29,30 DFT calculations are performed for the ground-state energy change of the lattice with respect to a small displacement in different directions for individual atoms. These DFT calculations essentially sample an energy-atom displacement surface, which contains the information on the interatomic forces between the perturbed atom and its neighbors, that is then fitted using the selected analytical functional forms for the empirical potential field. Apparently, assigning an appropriate analytical functional form for various interactions in the crystal is the key to this fitting process. For the ZnTe monolayer, Morse potential (see Table 1) is chosen to describe the stretching of Zn−Te bonds and the pairwise interaction between the second-nearest neighbors, because the Morse potential intrinsically contains anharmonic terms necessary for phonon−phonon scattering with a very simple functional form that only needs a small number of empirical parameters.31 To describe the three-body angular interaction, we use a cosine-squared form: E = K(cos θ − cos 28302
DOI: 10.1021/acs.jpcc.5b09527 J. Phys. Chem. C 2015, 119, 28300−28308
Article
The Journal of Physical Chemistry C Table 2. Lattice Constants (Å) and Elastic Constants (GPa) β-ZnTe(en)0.5 parameters lattice constants
f
elastic constantsg
ZnTe(zinc-blende)
DFT ST AM ZZ C11 C22 C33 C44 C55 C66 C12 C13 C23
17.555 (16.665) 5.783 (5.594) 4.451 (4.299) 43.1 32.5 51.6 18.0 5.8 8.2 14.1 7.4 18.0
MD a
17.369 5.751a 4.397a
17.775 5.732 4.461 41.9 39.1 55.3 21.5 1.38 7.69 14.9 8.83 18.9
experiment b
17.198 5.679b 4.343b
DFT 6.195
61.5 61.5 61.5 29.0 29.0 29.0 33.6 33.6 33.6
experiment 6.171
64.7c 64.7c 64.7c 36.0c 36.0c 36.0c 40.1c 40.1c 40.1c
c
6.103d
71.1e 71.1e 71.1e 31.3e 31.3e 31.3e 40.7e 40.7e 40.7e
a
Reference 35. bReference 48. cReference 49. dReference 50. eReference 51. fLattice constants in brackets are predicted by DFT-D2 calculations. Subscripts of elastic constants, 1−3 denote ST, AM, and ZZ axes correspondingly. Corresponding deformations to C44, C55, and C66 are shown in Figure 1. g
included in the potential field to reduce the parameters that need to be fitted, because dihedrals usually make only a very small contribution (
