TiC lattice dynamics from ab initio calculations

0.16 MB·English

by  Paweł T. Jochym

EPJ manuscript No.
(will be inserted by the editor)
TiC Lattice Dynamics from ab initio Calculations
P.T. Jochym, K. Parlinski, and M. Sternik
Instituteof Nuclear Physics, ul. Radzikowskiego 152, 31–342 Cracow, Poland
3
1 the dateof receipt and acceptance should beinserted later
0
2
Abstract. Ab initio calculations and a direct method have been applied to derive the phonon dispersion

n curvesandphonondensityofstatesfortheTiCcrystal.Theresultsarecomparedandfoundtobeinagood
a agreementwiththeexperimentalneutronscatteringdata.Theforceconstantshavebeendeterminedfrom
J the Hellmann-Feynman forces induced by atomic displacements in the 2×2×2 supercell. The calculated

5 phonondensityofstatessuggests thatvibrationsofTiatomsform acoustic branches,whereas themotion
2 of C atoms is confined to optic branches. The elastic constants have been found using the deformation
method and compared with theresults obtained from acoustic phonon slopes.
]
i
c PACS. 63.20.–e Phonons and vibrations in crystal lattices – 71.15.N Total energy calculations
s
-
l
r
t Thetransition-metalcarbidecompounds,ofwhichTiC values are generally in very good agreement with experi-
m
is a representative, are of big scientific and technologi- ment.
. cal interest because of their striking mechanical proper-
t
a ties, extreme hardness combined additionally with metal- In this paper we intend to extend first-principle cal-
m lic electrical and thermalconductivity. They have usually culations to describe the phonon dispersion curves and
- a rock-salt crystal structure – like some ionic crystals – phonon density of TiC. The method which was is based
d whilehavingpropertiestypicalformaterialswithcovalent on the total energy calculation and Hellmann–Feynman
n
bonding. It suggests that the bonding have some mixed (HF) forces. Dispersion relations are calculated by direct
o
nature, which makes these materials very interesting as method [9,10,11,12,13], in which the force constants of
c
[ an object of study. thedynamicalmatrixarecalculatedfromHFforces.Elas-
tic constants and bulk modulus have been estimated by
Over a large range of fractional carbon content from
1
straight evaluation of energy derivatives with respect to
v TiC to TiC0.5 [1] the TiC is stable in the NaCl structure,
7 with space group Fm¯3m. In practice, it is usually non- deformation.
7 stoichiometric and contains carbon vacancy defects. The
The energies and forces of the TiC crystal were cal-
0 phonon dispersioncurves of nearly stoichiometric TiC0.95 culated by the method of total energy minimization, us-
6 have been measured along main symmetry directions by
1. Pintschoviusetal[2]usingtheinelasticneutronscattering inganorm–conservingpseudopotentials asanapproxima-
tion of the atomic core – valence electron interaction [14,
0 technique. In the same paper, the influence of the carbon
15]. For an excellent review of ab initio total–energy cal-
3 content on the phonon dispersion curves has been addi-
culations see Ref.[16] and further references given there.
1 tionallycheckedbymeasuringasamplewithlowercarbon
v: concentration TiC0.89. The results shown that the optic wFoitrhthpeerliaotdtiiccebdoyunnadmariycsccoanldcuitliaotniosnasntdh6e42a×t2o×m2sshuapsebreceenll
i modes could differ by 3.5%, while the acoustic modes are
X used. For optimizing the structure and for direct calcu-
less sensitive to carbon content. The elastic constants de-
lation of energy deformations 1 1 1 supercell has been
r rived from the measured acoustic dispersion curves agree × ×
a utilized.The ab initio totalenergycalculationsweredone
verywellwiththoseobtainedbyultrasonicmeasurements
with the CASTEP package [17] and standard pseudopo-
for TiC0.91 crystal [3]. tentialsprovidedwithinthispackage.Thesepseudopoten-
VeryinterestingpropertiesandsimplestructureofTiC
tialswereconstructedwithinLDAapproximation.TheTi
causedthateffortstounderstanditspropertiesusingfirst-
pseudopotential treats 3p electrons as belonging to the
principle calculations have been undertaken severaltimes
core.We haveemployednon–localvariantofthose poten-
[4,5,6,7].Theelectronicstructure,bulkmodulus,andelas-
tialsparametrizedinthereciprocalspacewith900eVand
tic constants have been found [8] by means of accurate
600 eV cut–off energy for 1 1 1 and 2 2 2 supercell,
first-principles total-energy calculations using the full po- × × × ×
respectively. Tests which had been made with the cut-off
tential linear muffin-tin orbital method. The calculated
energy of 600 eV applied to the 1 1 1 supercell, showed
× ×
thatinthiscasethecohesiveenergyisonly0.2eVhigher,
Correspondence to: [email protected] andthelatticeconstantislongerby0.0006˚A(0.01%),com-
2 P.T. Jochym et al.: TiC Lattice Dynamicsfrom ab initio Calculations
paring to the 900 eV cut-off. Therefore,we have used 600 Cummulant force constants appear as a result of peri-
eV cut-off energy for a larger supercell calculation. odic boundary conditions imposed on the supercell. To
Ageneralized gradient approximation (GGA)hasbeen calculate HF forces an atomic configuration with a sin-
used for the exchange energy term of the valence states gledisplacedatommustbeminimizedwithrespecttothe
of the Hamiltonian [18]. Although this procedure is not electronic part only. Each of such runs provides 3n HF
quite consistent, there are indications [19] that LDA gen- forces, where n = 8 or 64 is the number of atoms in the
erated pseudopotentials used with GGA exchange term 1 1 1 or 2 2 2 supercell, respectively. Two runs with
× × × ×
incrystalsaresufficiently accuratefor energiesandstruc- displacementsalongz direction,oneforTiandsecondfor
tures. The integration over the Brillouin zone has been C were performed.To minimize the systematic errorstwo
performed with weighted summation over wave vectors other runs with negative displacements were carried on.
generated by Monkhorst-Pack scheme [20] using 0.1 ˚A−1 Thedisplacementamplitudesweresetto1%ofthelattice
and 0.06 ˚A−1 grid sizes which leads to 4 and 32 wave constant. This amplitude has been selected after number
vectors for 1 1 1 supercell and to 1 and 4 wave vectors of careful tests in which displacements ware varied from
for2 2 2su×per×cell,respectively.Nosignificantdifference 0.1% to 1.6%. It has been proven that the anharmonic
due t×o d×ifferent number of wave vectors, 1 or 4, has been contributionsarenegligableupto 1%of displacement.At
foundincalculatedphononfrequencies.Theresultsofdis- smallerdisplacements,closeto0.1%theHFforcesbecame
persion curve presented here are generated from the 0.06 toosmallandtheuncertaintyofforceconstantevaluation
˚A−1 data set with 4 wave vectors. is not acceptable.
SinceTiCisaweaklymetalliccompounditcouldhave
been appropriate to use a smearing of electron levels in
105
the Brillouin zone integration. This procedure is compu-
tationally more expensive. Thus we have performed a se-
riesofteststocomparelatticeconstants,bulkmoduliand 104
phononsinΓ andXpointsobtainedusinggaussiansmear-
ingintherangefrom4to0.1eVanddenserk-spacegrids 103
(4,14,32,63k-points)foroneunitcellconfiguration.The
comparison of results shown in the Table 1 indicate that, nt
for sufficiently dense k-space grid used, these quantities sta 102
n
are not very sensitive (specially the lattice constant and co
bulk modulus) to the choice of method of Brillouin zone orce 101
integration. Thus we have decided to use much simpler F
and faster non-metallic approach in the final 2×2×2 su- 100
percell calculation. Note also that results obtained in the
metallic regime lead to larger divergency of phonon fre-
quencies from experimental data. 10−1
The minimization of the total energy with respect to
the lattice constant was performed within the CASTEP 10−2
0.5 1.0 1.5
minimizermodulewith1 1 1supercellandcross-checked
Interaction distance [lattice units]
× ×
with 2 2 2supercell.The equilibriumlattice constantis
a = 4.×344×8˚A which could be compared with the exper- Fig. 1. Absolute eigenvaluesof theforce constant matrices of
imental value a = 4.3269˚A[21]. Such a good agreament theTiC crystal as a function of interatomic distance.
may be coincidental, since inclusion of 3p electrons in the
valence band as well as gradientcorrectionsin core states
The data of 4 3n HF forces F (n,ν) and four dis-
tends to change the equilibrium lattice constant [22,23, × i
placements u (m,µ) form an overdeterminedset of equa-
24]. j
tions, Eq.(2), for the force constants. This system was
The HF forces are defined as
solvedby the singularvalue decompositionalgorithm[11,
25], which automatically provides the least-squares solu-
F (n,ν)= ∂E /∂R (n,ν) (1)
i tot i tion. For 2 2 2 supercell the 768 components of the HF
−
× ×
forces lead to 33 non-zero independent parameters of the
wheren,ν aretheindicesoftheunitcellandatomwithin force constants. Knowledge of these force constants al-
the unit cell, respectively, and i is the Cartesian compo- lows one to test the range of interaction potential. For
nent. At extremum all HF forces vanish. Non-zeroth HF that, similarly to procedure used in [13], we have calcu-
forces arise, when a single atom (m,µ) is displaced by lated the absolute values of the eigenvalues of the 3 3
uj(m,µ) from its equilibrium position. The arising forces Φij(n,ν;m,µ) matrices and plotted them against a d×is-
are related with the cummulant force constants Φij by a tance between atoms (n,ν) and (m,µ). The result for
relation [11,12] 2 2 2supercellisdepictedinFig.1.Itshowsthatthein-
× ×
teractiondropsdownatleastbytwoordersofmagnitudes
Fi(n,ν)= X Φij(n,ν;m,µ)uj(m,µ) (2) withinthesizeofthesupercell.Thescatterofpointsistoo
− largetoassignanyrulewhichcouldgovernthedecayingof
m,µ,j
P.T. Jochym et al.: TiC Lattice Dynamicsfrom ab initio Calculations 3
Quantity Exper. Non-metallic Metallic
# k 4 4 14 32 63
a0 4.327 4.345 4.35 4.342 4.345 4.342
B 2.4 2.42 2.41 3.527 2.418 2.427
ωΓ 16.3 16.5 19.0 11.8 15.6 15.6
ω1 8.5 8.5 10.2 5.1 9.2 8.5
X
ω2 10.8 10.0 11.1 10.3 10.7 10.7
X
ω3 16.5 17.0 18.2 12.8 15.7 15.3
X
ω4 18.5 17.8 19.3 16.7 18.1 17.7
X
Table 1. Comparison of results from metallic and non-metallic approach
theforceconstantsparameterswithdistance.FromFig.1 Thedispersioncurvescalculatedonthebasisof2 2 2
× ×
itisclearthatthe1 1 1supercellistoosmalltobeused supercellareshowninFig.2.Theycorrespondtotemper-
× ×
for calculation of dispersion curves since the considered ature T=0. At the same Figure 2 we compare the cal-
force constants terminate at radius of a√3/2=0.866a. culated dispersion curves with the experimental phonon
The knowledge of the cummulant force constants al- frequencies measured by the inelastic neutron scattering
lows one to define an approximate dynamical matrix, in [2].Differentsymbolsofexperimentalpointscorrespondto
which the summation over all atoms is confined to those differentphononpolarizationsestablishedexperimentally.
atoms which reside within the volume of the supercell. The experimental points are about ∆ω =0.05ω lower
The approximatedynamical matrix becomes equal to the then the calculated values. One of the reasons for such
conventionalone atdiscrete wavevectorsk givenbe the situation could be that experimental values are measured
L
equation exp(2πikL L) = 1 [11], where L = (L1,L2,L3) at room temperature and anharmonic effects might dem-
arethelatticevectors· ofthesupercell.Atthewavevectors inish the phonon frequencies. Another reason may be re-
k the phonon frequencies ω2(k ), calculated by diago- lated with a non-stoichiometric concentration of carbon
L L
nalization of the approximate dynamical matrix, are the in measured samples, which diminish in average phonon
same as those calculated from the exact dynamical ma- frequencies as well. The efect that experimental phonon
trix. For 2 2 2 supercell the exact wave vectors are at points are at lower frequencies then the calculated ones
Γ, X, L, W×, ×and two other wave vector points, namely, appears as well in other ab-initio calculations [9,11].
the midpoint betweenΓ andX along 1,0,0 and 1,1,1 The theory of lattice dynamics provides summation
directions. The advantage of the abohve desicribedhdirecti ruleswhichfollowfromtranslationalandrotationalinvari-
method is that it does not impose any limit to the range ancesofthe crystal.ForTiCtheseinvariancesleadtotwo
ofinteraction.When the supercellsize issmaller thanthe equations, which can be used to set the values of on-site
range of interaction, the direct method interpolates the forceconstantparameters.Iftheseconditionsaresatisfied
dispersion curves between the exact points. the acoustic phonon branches at k=0 would point to ω =
0. Our dispersion curves, without imposed translational-
rotational invariances point at k= 0 to ω = 0.31THz.
Γ ∆ X K Σ Γ Λ L Q W Z X Since the invariance conditions should be satisfied within
25 the ab initio and direct methods, we have not corrected
ourdispersioncurvesdisplayedinFig.2 accordingto this
effect.
20 Bysamplingthedynamicalmatrixformanywavevec-
tors, one can calculate the phonon density function g(ω),
z] and the partial phonon density functions gx,Ti(ω) and
H15
T g (ω). The g(ω) describes the number of phonon fre-
y [ x,C
c quencies in an interval around ω, while g and g
n x,Ti x,C
e
u specify the number of phonon frequencies in the interval
eq10
Fr aroundωbutmodifiedbythepopulationofthexdisplace-
ment component of either Ti or C atom. The density of
statesfunctionsareshowninFig.3.Theyareconvention-
5
ally normalized to g(ω)dω = 1 and g (ω)dω = 1/6,
R R x,α
where α= Ti or C. The curves show that motions within
0 acoustic dispersion curves are almost entirely due to Ti
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
atoms.ThesumofTitaniumdensityfunctions(g (ω)+
Wave vector [1/A] x,Ti
g (ω)+g (ω)) fits to the totaldensity function g(ω)
y,Ti z,Ti
Fig. 2. Phonon dispersion relations of TiC crystal calculated below 14 THz. Part of g(ω) above the gap is mainly due
from 2×2×2 supercell (full line). The same dispersion curves to C atoms. Thus, TiC form a rather special crystal, in
divided by a factor 1.05 (dashed line). Experimental points which the Ti atoms, as heavier, form a frame for elas-
taken from Ref. [2]. tic motion, and the C atoms vibrate within the optical
4 P.T. Jochym et al.: TiC Lattice Dynamicsfrom ab initio Calculations
0.25 interactionsexceeding the size of the supercell,except for
Total the special points Γ, X, L, W mentioned earlier.
Carbon
Titanium
0.20
Result B c11 c12 c44
Exp.TiC0.91 [3] (2.422) 5.145 1.060 1.788
sity 0.15 Exp.TiC [27] 2.4 5.00 1.13 1.75
en Exp.TiC0.95 [2] (2.5) 5.4 1.1 1.8
n d Exp.TiC0.89 [2] (2.5) 5.2 1.1 1.8
o
n Calculations [8] 2.2 4.7 0.97 1.67
o
Ph 0.10 Ours,acoustic modes (3.57) 5.7 2.5 1.51
Ours,deform. energy 2.42 5.39 0.94 1.92
Ours,stress–strain (2.46) 5.21 1.09 1.93
0.05 Table 2. Experimental and calculated bulk modulus B and
elasticconstantscij oftheTiCcrystal,inunitsof1011Nm−2 =
1MBar.ValuesinparenthesisarecalculatedfromB= 31(c11+
0.00 2c12).
0 5 10 15 20 25
Frequency [THz]
Fig. 3. Phonon densityspectrumg(ω)oftheTiC crystaland
Toruleoutthatthediscrepancyforc12doesnotfollow
partialphonondensityspectragi,α(ω),whereα=CorTiand fromab-initiocalculations,wehavederivedthebulkmod-
i = x, y or z are the displacement directions. Experimental
ulus and the elastic constants by two conventional meth-
points taken from Ref. [2].
ods [8]. The bulk modulus was calculated in the 1 1 1
× ×
supercellfromthe totalenergyasa functionofthe lattice
constant.Thisdatahavebeenfittedtoathirdorderpoly-
modes. We add, however,that the magnitude ofthe force nomial. Hence, the derivative ∂E, and the bulk modulus
constants between Ti - C and between Ti - Ti and C - C ∂V
B were calculated from the relation: B = V ∂E, where
at the same distances remains of the same order. Thus, − ∂V
V is the volume of TiC crystallographic unit cell. The B
TiC does not consist of two weakly bounded subsystems.
values aregiveninTable 2 andthey agreequite wellwith
InFig.3theexperimentallydeterminedphonondensityof
experimental data.
states[2,26]isalsoshown.Theagreementwithinacoustic
region is quite qood. In optic part the experimental reso- A calculationof c11, c12, c44 elastic constantsinvolved
asimilarprocedure,inwhichthe crystalwasdeformedby
lution of 2.5THz and the nonsoichiometry of carbon in
≈ lengthening and sheer deformations. Deformations from
the sample lead to broadening of the optic band.
0.5%to3%inlengthandfrom1to5degreesinanglehave
Theelasticconstantscanbeevaluatedfromthe slopes
been used. The results of small and large deformations
of acoustic branches of phonon dispersion curves. Here,
are consistent. The deformed lattices have space groups
we haveusedthe invarianceconditionsto set the acoustic Fm¯3m,I4/mmmandI/mmm,forbulkexpansion,elonga-
modes to ω = 0 at k= 0. For the fcc structure, the set of
tion along z, and sheer modes. In all these deformed lat-
seven equations for three unknown variables c has to be
ij tices all atoms remain in the high-symmetry sites which
solved:ρvα2,L =c11;ρvα2,T =c44;ρvβ2,L =c11+4c44+2c12; guaranties an extremum of potential energy at those po-
ρvβ2,T =c11+c44−c12; ρvγ2,L =c11+2c44+c12; ρvγ2,Th = sitions. Thus, in the case of cubic TiC, one is entitled
c11 − c12; ρvγ2,T3 = c44, where the velocities vi are the not to carry on the relaxation of internal degrees of free-
slopes of appropriate acoustic branches. Here, α,β and domduringsupercelldeformations.Thecalculatedelastic
γ indexes describe three main directions in the recipro- constants are compared to experimental data in Table 2.
cal space 0,0,1 , 1,1,1 and 0,1,1 , respectively. The TheyagreewithcalculatedvaluesgiveninRef.[8]andare
h i h i h3 i
densityofTiCρ=4849.6kg/m wasfoundusingthe cal- in good agreement with experiment.
culated lattice constant. The solution of this set of equa- We have also checked the above results using stresses
tions,bymeansoftheleast-squaresfit,resultsinc11 =5.7, calculated by CASTEP. The elastic constants were found
c12 = 2.5 and c44 = 1.51 MBar. Comparing these values as first derivatives of the stress–strain relationship. For
withexperimentaldata(Table2)onenoticesconsiderable that,wehaveusedthesamerunsofelectroniccalculations
discrepancies, especially for c12. Closer look to dispersion asfortheelasticconstantderivationgivenaboveusingde-
curves,Fig.2,assuresusthatthetransverseacousticmode formation dependence of energy. The results are included
in 1,1,0 direction and the longitudinalin 1,1,1 direc- in Table 2 and they show slightly better agreement with
h i h i
tion do not follow exactly the experimental points. We the experimental data.
have checked, that such discrepancy of slope is sufficient Insummary,wehaveshownthatthe ab initio calcula-
to produce the observeddeviations of the calculated elas- tions of HF forces, together with the direct method lead
tic constants, specially of c12. The discrepancy in slopes to satisfactory description of phonons in the TiC crys-
is so small that it could be assigned to the missing force tal. As usual the frequency of the calculated dispersion
constants describing mainly the long range Coulomb in- curves are slightly too high. The acoustic modes are not
teraction. Indeed, the direct method did not handle well sufficiently accurate, thus the elastic constants calculated
P.T. Jochym et al.: TiC Lattice Dynamicsfrom ab initio Calculations 5
from them differ from those found through the supercell 25. W.H.Press,B.P.Flamery,S.A.Tenkolsky,andW.T.Vet-
deformations. The acoustic part of density of states fits terling,NumericalRecipes(CambridgeUniv.Press,1988),
well to neutron experimental data. p.301
26. F. Gompf and W. Reichardt,1977, unpublished.
27. M. M. Choy at. al., Numerical Data and Functional Re-
The computational advices of J. Oleszkiewicz and the use of lationships in Science and Technology (Springer, Berlin,
facilities of the Computer and Network Center of the Jagel- 1979), vol.11.
lonian University, Cracow, where the calculations have been
done, are kindly acknowledged. This work was partially sup-
ported by the Polish State Committee of Scientific Research
(KBN),grant No 2 PO3B 145 10.
References
1. L.E. Toth, Transition Metal Carbides and Nitrides (Aca-
demicPress, New York 1971).
2. L. Pintschovius, W. Reichardt and B. Scheerer, J. Phys.
C: Solid State Phys.11, 1557 (1978).
3. R. Chang and L. J. Graham, J. Appl. Phys. 37, 3778
(1966).
4. A. Neckel, R. Eibler, P.Weinberger and K. Schwartz, J.
PhysC 9, 579 (1976).
5. V.P.Zhukov,V.A.Gubanov,O.Jepsen,N.E.Christensen
andO.K.Andersen,J.Phys.Chem.Solids49,841(1988).
6. D.L.PriceandB.R.Cooper,Phys.Rev.B39,4945(1989).
7. J.H¨aglund,A.F.Guillerment,G.GrimvallandM.K¨orling,
Phys.Rev.B 48, 11685 (1993).
8. R. Ahuja, O. Eriksson, J. M. Wills, B. Johansson, Phys.
Rev.B 53, 3073 (1996).
9. W.Frank, C.Els¨asser, and M,Fa¨hnle, Phys.Rev.Lett. 74,
1791 (1995).
10. G.Kresse,J.Furthmu¨ller,andJ.Hafner,Europhys.Lett.32,
729 (1995).
11. K. Parlinski, Z. Q. Li and Y. Kawazoe, Phys. Rev. Lett.
78, 4063 (1997).
12. G. J. Ackland, M. C. Warren and S. J. Clark, J. Phys.:
Condens. Matter 9, 7861–7872 (1997).
13. M.H.F. Sluiter,M. Weinert, Y.Kawazoe, Europhys.Lett.
43, 183 (1998).
14. J.S.Lin,A.Qteish,M.C. PayneandV.Heine,Phys.Rev.
B 47, 4174 (1993).
15. J. Goniakowski, J.M. Holender, L.N. Kantorovich, M.J.
Gillan and J.A. White, Phys. Rev.B 53, 953 (1996).
16. M.C. Payene, M.P. Teter, D.C. Allan, T.A. Arias, J.D.
Joannopoulos, Rev. of Mod. Phys. 64, 1045 (1992)
17. CASTEPisapartofCerius2systemdevelopedbyMolec-
ular Simulations.
18. J.P.Perdew,J.A.Chevary,S.H.Vosko,K.A.Jackson,M.R.
Pederson, D.J. Singh, C. Fiolhais, Phys. Rev. B 46, 6671
(1992).
19. J. Goniakowski, J. M. Holender, L.N.Kantorovich,M. J.
Gillan, J. A. White, Phys.Rev.B 53, 957 (1996)
20. H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13, 5188
(1976).
21. I. N´aray-Szab´o, Inorganic Crystal Chemistry (Akad´emiai
kiad´o, Budapest, 1969), p.301.
22. A. Garc´ıa, C. Els¨asser, J. Zhu, S.G. Louie, M.L. Cohen,
Phys.Rev.B 46, 9829 (1992)
23. Y.M.JuanandE.Kaxiras,Phys.Rev.B48,14994(1993)
24. K. Seifert, J. Hafner, J. Furthmu¨ller, G. Kresse, J. Phys.
Condens. Matter 7, 3683 (1995)