参数修改0PhysRevB.81.205441

Optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene L. Lindsay Department of Physics, Boston College, Chestnut Hill, Massachusetts 02467, USA and Department of Physics, Computer Science, and Engineering, Christopher Newport University, Newport News, Virginia 23606, USA D. A. Broido Department of Physics, Boston College, Chestnut Hill, Massachusetts 02467, USA �Received 3 March 2010; revised manuscript received 10 May 2010; published 27 May 2010� We have examined the commonly used Tersoff and Brenner empirical interatomic potentials in the context of the phonon dispersions in graphene. We have found a parameter set for each empirical potential that provides improved fits to some structural data and to the in-plane phonon-dispersion data for graphite. These optimized parameter sets yield values of the acoustic-phonon velocities that are in better agreement with measured data. They also provide lattice thermal conductivity values in single-walled carbon nanotubes and graphene that are considerably improved compared to those obtained from the original parameter sets. DOI: 10.1103/PhysRevB.81.205441 PACS number�s�: 63.22.Rc, 65.80.�g, 63.20.�e, 65.40.�b I. INTRODUCTION It is well known that carbon-based materials such as dia- mond, graphite, graphene, and carbon nanotubes possess the highest known thermal conductivities.1–9 This has motivated intense theoretical effort to understand the thermal transport properties of these systems.10–30 Such investigations com- monly employ either molecular-dynamics simulations or Boltzmann transport equation �BTE� approaches. In both cases accurate representation of the interactions between at- oms is required. Harmonic interatomic force constants �IFCs� are required to obtain phonon frequencies and eigen- vectors, and acoustic-phonon velocities. Real-space anhar- monic IFCs are also needed to describe the phonon-phonon scattering that is responsible for the intrinsic thermal resis- tance. Rigorous first-principles and force-constant approaches have been employed in graphene and single-walled carbon nanotubes �SWCNTs� to calculate the harmonic IFCs in these systems, and highly accurate phonon dispersions have been obtained over the entire Brillouin zone.31–33 Anhar- monic IFCs have also been obtained recently in these sys- tems using density-functional perturbation theory �DFPT� and used successfully in phonon lifetime calculations.34,35 However, DFPT calculations of the real-space anharmonic IFCs are considerably more difficult so no commonly avail- able ab initio DFPT package, such as QUANTUM ESPRESSO,36 VASP,37 and ABINIT �Ref. 38� generates them. In contrast, empirical interatomic potentials �EIPs� provide these anhar- monic interactions in a conveniently extractable form and so are frequently used in thermal transport calculations in SWCNTs and more recently in graphene.23,24,30 The most commonly used EIPs are those developed by Tersoff39,40 and Brenner.41,42 The Tersoff and Brenner EIPs are short range and so can- not accurately fit the graphene dispersion over the entire Brillouin zone. However, the thermal conductivity depends more sensitively on the accuracy of acoustic phonon fre- quencies near the zone center where the longitudinal- and transverse-acoustic �LA and TA� velocities and the quadratic curvature of the out-of-plane acoustic �ZA� branch are deter- mined. On the other hand, only weak thermal excitation of the optic phonons and acoustic phonons near the Brillouin- zone boundary occurs around room temperature because of the large graphene Debye temperature ��2000 K�. There- fore, accurate fitting to the corresponding phonon frequen- cies is less important. The original parameter sets of the Tersoff and Brenner EIPs do not accurately reproduce the phonon dispersions of graphene, as has been noted previously.43 In particular, they do not accurately obtain the velocities of the three acoustic branches near the center of the Brillouin zone, thus, also misrepresenting these properties in SWCNTs. We present here optimized parameter sets for the Tersoff and Brenner EIPs, which better represent the lattice dynamical properties of graphene. We test these optimized parameter sets on the lattice thermal conductivities of SWCNTs and graphene, and we find them to yield values that are in much better agree- ment with available data. In Sec. II we briefly describe the Tersoff and Brenner EIPs and the approach we took in opti- mizing the empirical parameters. Section III presents the op- timized parameter sets and compares the graphene phonon dispersions obtained from the original and optimized param- eter sets. It also examines the phonon thermal conductivities of SWCNTs and graphene obtained from old and new pa- rameter sets for both EIPs. Section IV presents a summary and conclusions. II. THEORY The convenience of the Tersoff and Brenner EIPs comes from their rather simple, analytical forms and the short range of atomic interactions. For carbon-based systems, the Tersoff model has nine adjustable parameters �listed in Table I� that were originally fit to cohesive energies of various carbon systems, the lattice constant of diamond, and the bulk modu- lus of diamond. The Brenner EIP is based directly on the Tersoff EIP but has additional terms and parameters which PHYSICAL REVIEW B 81, 205441 �2010� 1098-0121/2010/81�20�/205441�6� ©2010 The American Physical Society205441-1 sodom 高亮 sodom 高亮 sodom 高亮 sodom 高亮 sodom 高亮 allow it to better describe various chemical reactions in hy- drocarbons and include nonlocal effects �parameters listed in Table II�. A. Tersoff model The analytical form for the pair potential, Vij, of the Ter- soff model is given by the following functions with corre- sponding parameters listed in Table I:39,40 Vij = f ijC�aijf ijR − bijf ijA� , �1a� f ijR = Ae−�1rij , �1b� f ijA = Be−�2rij , �1c� where rij is the distance between atoms i and j, f ijA and f ijR are competing attractive and repulsive pairwise terms, and f ijC is a cutoff term which ensures only nearest-neighbor interac- tions. The term, aij, is a range-limiting term on the repulsive potential that is typically set equal to 1. We do so here. The bond angle term, bij, depends on the local coordination of atoms around atom i and the angle between atoms i, j, and k, bij = �1 + �n�ij n �−1/2n, �2a� �ij = � k�i,j f ikCgijke�3 3�rij − rik�3, �2b� gijk = 1 + c2 d2 − c2 d2 + �h − cos��ijk��2 , �2c� where �ijk is the angle between atoms i, j, and k. This bond angle term allows the Tersoff model to describe the strong covalent bonding that occurs in carbon systems, which can- not be represented by purely central potentials. This angle- dependent term also allows for description of carbon systems that bond in different geometries, such as tetrahedrally bonded diamond and the 120° tribonded graphene. B. Brenner model The Brenner potential for solid-state carbon structures is given by the following functions with corresponding param- eters listed in Table II:41,42 Vij = f ijC�f ijR − b¯ijf ijA� , �3a� f ijR = �1 + Q rij �Ae−�rij , �3b� f ijA = � n=1 3 Bne−�nrij , �3c� where many of the terms are similar to the Tersoff model described above and the bond angle term, b¯ij, is given by b¯ij = 1 2 �bij �−� + bji �−�� + ij RC + bij DH , �4a� bij �−� = �1 + � k�i,j f ikCgijk�−1/2, �4b� gijk = � i=0 5 �i cos i��ijk� . �4c� Here, bij �−� depends on the local coordination of atoms around atom i and the angle between atoms i, j, and k, �ijk. This pi-bond function is symmetric for graphene, graphite, and diamond, bij �−� =bji �−� . The coefficients, �i, in the bond- bending spline function, gijk, were fit to experimental data for graphite and diamond and are also listed in Table II. The term, ij RC , accounts for various radical energetics, such as vacancies, which are not considered here; thus, this term is taken to be zero. The term, bij DH , is a dihedral bending func- tion that depends on the local conjugation and is zero for diamond but important for describing graphene and SWCNTs. This dihedral function involves third-nearest- neighbor atoms and is given by41,42,44 bij DH = T0 2 �k,l�i,j f ik C f jlC�1 − cos2� ijkl�� , �5� where T0 is a parameter, f ijC is the cutoff function, and ijkl is the dihedral angle of four atoms identified by the indices, i, j, k, and l, and is given by TABLE I. Original parameters for the Tersoff EIP for carbon- based systems given in Ref. 40. A=1393.6 eV B=346.74 eV �1=3.4879 Å−1 �2=2.2119 Å−1 �3=0.0000 Å−1 n=0.72751 c=38049.0 �=1.5724�10−7 d=4.3484 h=−0.57058 R=1.95 Å D=0.15 Å TABLE II. Original parameters for the Brenner EIP for solid- state carbon-based structures given in Ref. 42. Also listed are the coefficients for the fifth-order polynomial spline, gijk, described in Ref. 42. A=10953.544162170 eV B1=12388.79197798 eV B2=17.56740646509 eV B3=30.71493208065 eV �=4.7465390606595 Å−1 �1=4.7204523127 Å−1 �2=1.4332132499 Å−1 �3=1.3826912506 Å−1 Q=0.3134602960833 Å R=2.0 Å D=1.7 Å T0=−0.00809675 �0=0.7073 �1=5.6774 �2=24.0970 �3=57.5918 �4=71.8829 �5=36.2789 L. LINDSAY AND D. A. BROIDO PHYSICAL REVIEW B 81, 205441 �2010� 205441-2 cos� ijkl� = �� jik · �� ijl, �6a� �� jik = r� ji � r�ik r� ji r�ik sin��ijk� , �6b� where �� jik and �� ijl are unit vectors normal to the triangles formed by the atoms given by the subscripts, r�ij is the vector from atom i to atom j, and �ijk is the angle between atoms i, j, and k. In flat graphene, the dihedral angle, ijkl, is either 0 or � and the dihedral term is subsequently zero.44 Bending of the graphene layer leads to a contribution from this term. Some of the main differences when compared to the Ter- soff EIP are: the Brenner EIP includes two additional expo- nential terms with corresponding adjustable parameters in the attractive pairwise term, it includes a screened Coulomb term in the repulsive pairwise term, it uses a fifth-order poly- nomial spline between bond orders for diamond and graph- ite, and it includes a dihedral bending term for bond energies which plays a role in SWCNTs and graphene. C. Parameter optimization We have implemented a 2 minimization procedure for each of these EIPs. A numerical algorithm was developed to minimize 2 given by45,46 2 = � i ��i − �exp�2 �exp 2 �i, �7� where �exp are experimental parameters used in the fitting process, �i are the corresponding values obtained using each potential, and �i are weighting factors that determine the relative importance of �i in the fitting procedure. In minimizing 2, the greatest significance was given to the phonon frequencies, ��, and the zone-center acoustic ve- locities, v¯�=d�� /dq� , of graphene in the high-symmetry di- rections. Here, �= �q� , j� designates a phonon with wave vec- tor, q� , in branch, j. The phonon frequencies are determined by diagonalization of the dynamical matrix for a given q� in the two-dimensional graphene Brillouin zone. The dynamical matrix is D�� ����q�� = 1 m�m�� � �� ��� 0�,����eiq� ·R � ��, �8� where �� designates the �th atom in the �th unit cell, m� is the mass of the �th atom, R� � is the lattice vector for the �th unit cell, and � and � are Cartesian components. In Eq. �8�, ��� 0�,���� are second-order IFCs which are determined by each EIP. We attached the most significance to the phonon frequen- cies and zone-center acoustic velocities because of the im- portant roles that they play in thermal transport calculations. The parameters calculated by each EIP for graphene were compared to the corresponding experimental parameters for in-plane graphite. First-principles calculations of the phonon dispersions in graphene have found excellent agreement with measured in-plane dispersion of graphite.31,32 This is consis- tent with the known weak coupling between the graphene layers in graphite. The cohesive energies and lattice con- stants of graphite and diamond were also considered in the minimization procedure but were given lesser weight. Using this minimization procedure for the Tersoff EIP, we found that simply modifying the h parameter, which helps to adjust the strength of the bond-angle term, provided vast improvement to the optical branches of the phonon disper- sion, while also improving the fit to the TA zone-center ve- locity. Adjustment of the B parameter associated with the strength of the attractive term was required to retain a decent fit to the experimental lattice constants for both graphite and diamond. Simple adjustment of the h parameter has previ- ously been found to significantly improve the fits to the mea- sured linear expansion coefficient47 and thermal conductivity of bulk silicon.48 A slightly different approach was taken in optimizing the Brenner model to better fit the phonon spectrum. The Bren- ner potential includes a dihedral bonding term, not included in the Tersoff model, which plays a role in graphene systems. The dihedral term, Eq. �5�, has a single adjustable parameter, T0, which was determined by fitting the lattice constant of a hypothetical three-dimensional, hexagonal system whose di- hedral bond angles are � /2.42,49 Changing this parameter alone simply alters the out-of-plane acoustic and optic �ZA and ZO� modes in graphite but leaves the other phonon modes unaltered and has no effect on diamond since T0 is zero for the tetrahedral configuration. We have chosen to adjust this parameter to better fit the phonon frequencies for the ZA branch in graphite. The six coefficients, �i, in the fifth-order polynomial spline, gijk, �Eq. �4c�� used to represent the bond-bending term in the Brenner EIP introduce additional flexibility not available in the Tersoff potential. These coefficients were originally determined by fixing the values of the gijk and their first and second derivatives at 109° and 120° �corre- sponding to diamond and graphite bond angles� to match various experimental data.42 The values for gijk were deter- mined by fitting the cohesive bond energies of graphite and diamond while the second derivatives were chosen to fit the elastic constants, c11, for diamond and in-plane graphite. The first derivatives were simply chosen to suppress oscillations of the spline function. The coefficients, �i, are determined from these values for gijk. These coefficients fix the spline and its derivatives at the bond angles for diamond and graph- ite, and the function is interpolated between these angles. Thus, they can be adjusted to separately fit experimental data for graphite while leaving the representation for diamond unaltered. The bond angles for SWCNTs, though close to graphite, fall between these two structures and thus depend on both. We note that the bond angles for large diameter SWCNTs approach that of graphite. We restricted the parameter optimization of the Brenner potential to T0 and �i to avoid significantly altering the pre- vious fits to the extensive structural experimental data sets.41,42 The values for gijk at 109° and 120° were altered slightly to better fit the given experimental lattice constants for diamond and graphite. We then adjusted the values of the second derivatives of gijk in these materials to better fit the zone-center acoustic velocities and corresponding phonon frequencies. Upon plotting the fifth-order polynomial spline OPTIMIZED TERSOFF AND BRENNER EMPIRICAL… PHYSICAL REVIEW B 81, 205441 �2010� 205441-3 sodom 高亮 sodom 高亮 sodom 高亮 sodom 高亮 sodom 高亮 sodom 高亮 sodom 高亮 sodom 高亮 given by the old coefficients versus the optimized coeffi- cients no visual difference can be seen. However, these small changes introduce noticeable changes in the corresponding phonon dispersions in graphene. III. RESULTS AND DISCUSSION The optimized parameter sets for the Tersoff and Brenner EIPs are listed in Table III. These parameters provide im- proved fits to experimental phonon-acoustic velocities and frequencies without significantly altering fits to other struc- tural data. The calculated phonon dispersion for graphene as given by the Tersoff �Brenner� EIP is shown in Fig. 1 �Fig. 2� along with the corresponding measured in-plane phonon dis- persion for graphite.50,51 In each figure the black �red� lines correspond to the optimized �original� parameter sets. The original set of parameters for the Tersoff EIP gives higher values for the TA branch velocities in the high- symmetry directions compared to the measured data while giving values for the quadratic ZA frequencies that fall be- low experiment around the M point. The most obvious fail- ure of the Tersoff EIP in describing the phonon dispersion of graphene comes in the highest optical branches, as seen in Fig. 1. The measured in-plane upper optic modes at the � point are degenerate with a value of 300 THz, while the Tersoff potential gives a value of 470 THz, a discrepancy of nearly 40%. The Tersoff model using the optimized param- eter set more accurately describes these upper optic phonon branches while providing a decent fit to the acoustic veloci- ties and phonon frequencies. However, using these param- eters provides a poorer fit to the experimental dispersion for the out-of-plane ZO branch. The inability to simultaneously fit the acoustic branches and all of the optic branches is a consequence of the Tersoff potential’s short range with only second-nearest-neighbor interactions represented. Tewary and Yang43 obtained a better fit to the phonon dispersion in graphene using a longer-range EIP that extended to fourth- nearest neighbors. Los et al.52 have also developed a some- what more complicated long-range EIP based on the Brenner EIP, though, phonon dispersion results have not yet been published. The original set of Brenner parameters, while providing a much better description of the optic branches, does not ac- curately represent the zone-center velocities for all of the acoustic modes. With the original Brenner EIP parameters, the velocities of the TA branch in the high-symmetry direc- tions are too low by 30%, those of the LA branch 12% too small, and the dispersion of the ZA branch undershoots the data by as much as 20%. The optimized parameter set im- proves the fit to the zone-center acoustic velocities and most of the experimental ZA, TA, and LA dispersion data. This optimized set provides a somewhat worse fit to the optic dispersion. The inability to simultaneously fit both acoustic and optic branches is again most likely a consequence of the short range of this EIP. Table IV lists measured lattice constants, cohesive ener- gies, and acoustic velocities in the �→M direction for graphite and the �→X direction for diamond compared with those obtained from the Tersoff and Brenner EIPs using both the original and optimized parameter sets. For both the Ter- soff and Brenner EIPs, the optimized parameter sets provide TABLE III. Optimized parameters and coefficients for the Ter- soff and Brenner EIPs. All parameters not listed are unaltered from the original sets. Tersoff h=−0.930 B=430.0 eV Brenner T0=−0.0165 �0=0.0000 �1=−3.1822 �2=−19.9928 �3=−51.4108 �4=−61.9925 �5=−29.0523 0 100 200 300 400 500 � (T H z ) � �� � FIG. 1. �Color online� Phonon dispersion for graphene along high-symmetry directions obtained using the Tersoff EIP. Thick black lines correspond to the optimized parameter set �this work�; thin red lines correspond to the original parameter set. Squares �tri- angles� are in-plane experimental data points for graphite from Ref. 43 �Ref. 44�. 0 50 100 150 200 250 300 350 � (T H z ) � � � � FIG. 2. �Color online� Phonon dispersion for graphene as given by the Brenner EIP for high-symmetry directions. Thick black lines correspond to the optimized parameter set �this work�; thin red lines correspond to the original parameter set. Squares �triangles� are in-plane experimental data