Periodic Hartree−Fock Study of TiS2 - The Journal of Physical

J. Phys. Chem. 1996, 100, 15735-15747

15735

Periodic Hartree-Fock Study of TiS2 D. G. Clerc* and R. D. Poshusta Materials Science Program, Washington State UniVersity, Pullman, Washington 99164-4630

A. C. Hess EnVironmental Molecular Sciences Laboratory, MS K1-90, Pacific Northwest Laboratory, Richland, Washington 99352 ReceiVed: August 7, 1995; In Final Form: May 3, 1996X

The structural and electronic properties of titanium disulfide were investigated at the ab initio, all-electron, periodic Hartree-Fock level, using an extended basis set and a posteriori density functional correlation corrections to the total energy. Calculated lattice parameters, bulk modulus, linear moduli, cohesive energy, elastic constants, Mulliken populations, band structure, and density of states are reported. The pressureinduced semiconductor-to-semimetal phase transition in titanium disulfide is found to result from an indirect band overlap between a sulfur 3pz-based valence band at Γ and a titanium 3d-based conduction band at L. The pressure shift of this overlap in the metallic phase is predicted to be 30 meV/GPa. The phase transition is predicted to occur at a pressure of 84 GPa, which is far greater than the experimental value of 4.0 ( 0.5 GPa. This large error is due to the large overestimation of the zero-pressure optical band gap in the HartreeFock approximation. Correcting this energy gap by using empirical rigid-band shifts results in a predicted transition pressure between 1.5 and 4.0 GPa, which is consistent with experiment. Lattice parameters a, c, and z at the phase transition are predicted to be 5.42(3) Å, 3.284(5) Å and 0.2567(16), respectively.

1. Introduction TiS2 belongs to the family of layered transition metal dichalcogenides MCh2 (M ) transition metal, Ch ) S, Se, or Te), in which a hexagonally packed plane of M atoms is sandwiched between two similar planes of Ch atoms. The bonding within the Ch-M-Ch “layers” is of the fairly strong covalent type. Adjacent Ch-M-Ch layers are coupled to each other by weak van der Waals (VDW) forces, and the empty space between them is therefore called the “van der Waals gap”. The Ch-M-Ch layers are stacked in the c direction to give rise to the observed hexagonal structure analogous to that of graphite and many silicate clay minerals. The quasi two-dimensional nature of MCh2 compounds gives rise to a marked anisotropy in transport properties. The most important property, however, is that the weak interlayer bonding permits the intercalation of various guests between the layers. Depending on the guest species, unusual and dramatic changes in the physical properties of the host can occur, and this property has stimulated much research interest. There is also considerable technological interest in TiS2 as a host material for intercalation reactions with alkali metals.1 For example, the unusual affinity of TiS2 for lithium intercalation makes it one of the most promising solid cathode materials.2,3 TiS2 has also received attention as an oxygen catalyst for fuel cells,4 as a catalyst in the preparation of urethanes5 and linear polyesters,6 as a catalyst in the thermal transformation of dodecacarbonyltriiron (Fe3(CO)12),7 and as a solid lubricant.8 The electronic structure of TiS2 is fundamental to understanding its physical and chemical properties. For example, calculated energy bands and density of states can be used to interpret photoelectron spectra, to study the nonmetal-to-metal phase transition, and to predict the electronic properties of its intercalation compounds. These results can by analogy be used to interpret the properties of other group IV MCh2 compounds. X

Abstract published in AdVance ACS Abstracts, September 1, 1996.

S0022-3654(95)02264-7 CCC: $12.00

The electronic structure of TiS2 was disputed for many years because the observed metallic properties could be interpreted in terms of both semimetallic and extrinsic semiconductor models. More recent experimental evidence indicates that TiS2 is an extrinsic semiconductor.9 Because experimental results were long inconclusive in determining the electronic structure of TiS2, previous theoretical work was particularly relevant in elucidating its physical properties.10-14 This paper reports our theoretical investigation of TiS2. In this study we focus on the structural, elastic, and electronic properties of stoichiometric TiS2. The results are compared to both experimental data and to previous theoretical work. Future studies will focus on the intercalation of alkali metal atoms and other species into this layered material. 2. Method of Calculation Our calculations are performed using CRYSTAL92,15,16 which is a well-documented ab initio program previously applied to a range of ionic, covalent, and molecular crystals17-23 as well as to adsorbate/surface interactions.24 This method is an allelectron periodic Hartree-Fock (PHF) self-consistent-field (SCF) procedure that expands the ground state wave function as a linear combination of atomic orbitals. As is customary, the atomic orbitals are linear combinations of Gaussian functions collectively called the basis set. The titanium basis used in this study is derived from a 86411+3d set of Gaussian functions originally optimized in TiC.25 The sulfur basis is derived from a 8831 set originally optimized in the ground state sulfur atom.26 Both basis sets were then augmented with single uncontracted d-type polarization functions, which resulted in overall shell structures of 86411+(3-1)d for titanium and 8831+1d for sulfur. With the exception of core shells, the exponents and contraction coefficients of both basis sets were then reoptimized by minimizing the total PHF energy of TiS2 at its experimental geometry. The resultant exponents and coefficients are reported in Table 1. © 1996 American Chemical Society

15736 J. Phys. Chem., Vol. 100, No. 39, 1996

Clerc et al.

TABLE 1: Basis Set Gaussian Exponents (in au-2) and Contraction Coefficients for Ti and S in Zero-Pressure (Equilibrium) TiS2a atom orbital Ti

1s

2sp

3sp

4sp*° 5sp*° 1d°

S

2d*° 1s

2sp

3sp° 4sp*° 1d*°

exponent 228000 32450 6888.6 1802.4 543.2 187.44 73.19 30.45 553.4 132.18 43.61 17.02 7.260 2.376 28.3 11.24 4.656 1.865 0.7666 0.2960 6.201 1.567 0.4614NM-M 0.1345NM-M 11676.6 5866.1 4897.34 1428.5 476.364 164.255 58.9157 21.3337 703.25 175.182 63.9449 26.8889 11.6989 5.328 2.4613 0.8865 3.635 1.412 0.4797 0.1610NM-M 0.2980NM-M

s coefficient

p coefficient

0.000228 0.001929 0.011100 0.04999 0.17010 0.36916 0.4027 0.1445 -0.00546 -0.0704 -0.1177 0.2451 0.6708 0.286 0.0027 -0.1515 -0.7440 1.032 1.000 1.000

0.00853 0.06021 0.2133 0.3871 0.4021 0.239 -0.0271 -0.0767 0.1665 0.314 1.000 1.000

d coefficient

0.132 0.348 0.425 1.000 0.0007 0.001 0.0074 0.0283 0.0975 0.2721 0.4786 0.3011 0.0002 -0.0091 -0.0551 -0.1102 0.0632 0.4959 0.6041 0.1439 -0.0335 -0.0515 0.1985 1.000

0.0019 0.0133 0.0485 0.1448 0.311 0.4107 0.3519 0.092 0.0651 -0.2869 -1.5247 1.000 1.000

a

Asterisks denote atomic orbitals assumed to be unoccupied at the beginning of the calculations. Circles denote exponents and coefficients optimized in TiS2. Exponents reoptimized at the various volumes shown in the E(V) curve of Figure 8 are marked with “NM-M”.

Basis set completeness was estimated by computing changes in (1) the total PHF energy and (2) the computed c-axis, before and after adding the most diffuse titanium “d” shell to the basis set. In the first case, the basis set was optimized at the experimental geometry before and after adding the said function. The resultant decrease in energy, 6.60 millihartree, was small enough to indicate that the basis set in Table 1 is sufficiently complete. In the second case, the total correlation-corrected PHF energy as a function of c-axis length was computed. The minima in the two resultant E(c) curves differed by only 0.4 Å, or equivalently, by 0.70%. This small difference further indicated the non-necessity of additional basis set functions. In this study, the Fock matrix is diagonalized at 65 k points in the irreducible part of the first BZ (Brillouin zone). For metallic (high-pressure) TiS2, the Fermi energy (EF) is computed by interpolating between these eigenvalues at 417 k points. In all sections of this paper except the phase transition investigation, this interpolation is accomplished by fitting the bands to

four symmetrized plane waves. In the phase transition investigation, eight symmetrized plane waves are used. The present work couples the HF approach with density functional theory to produce TiS2 total energies, which are corrected for electron correlation. Such corrections are important in TiS2 since nonbonded orbital interactions play important roles in determining its physical properties. The correction is obtained by integrating three different correlation-only density functionals over all space using the converged HF electron density. These functionals are Wigner’s local correlation formula gradient-corrected following Levy’s method (WL),27 Perdew’s general gradient approximation (P),28 and the ColleSalvetti gradient-dependent correlation energy functional (CS).29,30 Results are reported as purely Hartree-Fock (HF) or HartreeFock with correlation correction (HF+CC), where CC may be WL, P, or CS. In this study, we compute TiS2 geometries using several different methods that differ by the particular degrees of freedom chosen to minimize the total HF+CC energy. In the case of TiS2, the degrees of freedom may include all symmetryindependent lattice parameters (a and c) and atomic coordinates (z). (Note: a, c, and z are defined in the next section.) The four methods we use are as follows. (1) The predicted zero-pressure (equilibrium) geometry is calculated by fitting a set of energies to quadratic. Each member of this set corresponds to a distinct geometry (i.e. distinct aand c-axes) having the minimum HF+CC total energy with respect to z. The zero-pressure geometrysdefined by a*, c*, and z*sthus includes the effects of internal relaxation of the atoms. Such optimizations are unconstrained and hence are referred to as “U-optimizations”. This calculation is fully described in ref 31. (2) The predicted geometries at various degrees of crystal dilation are calculated by holding the c-axis fixed while minimizing the HF+CC total energy with respect to the other degrees of freedom. These processes will be called “Coptimizations”. (Note: in the special case where c ) c*, the C-optimization produces essentially the same result as the U-optimization.) (3) The predicted geometry of TiS2 under pressure at volume V (V ) x3a2c/2) is calculated by minimizing the HF+CC total energy with respect to the a-axis and z while requiring c ) 2V/(x3a2). This process is referred to as a “V-optimization”. (4) We report energies obtained at various VDW gaps while requiring the TiS2 intralayer geometry to be independent of c for small compressions and arbitrary expansions. This stipulation is called the “rigid-layer (RL) approximation” and it will be fully described in a future section of this paper. 3. Crystal Structure X-ray crystallography shows that TiS2 is trigonal with the space group P3hm1.32 The hexagonal lattice parameters are a ) 3.4073(2) Å and c ) 5.6953(2) Å, with ratio c/a ) 1.6715. The ab plane coincides with the titanium plane, and there is one formula unit per unit cell (Z ) 1). The titanium and sulfur fractional atomic coordinates are Ti(0, 0, 0), S(1/3, 2/3, z), and S(2/3, 1/3, -z), where the internal coordinate z ) 0.2501(4). The projection of one TiS2 layer onto the ab plane is shown in Figure 1A. The thickness of the layer is 2cz, and the distance between adjacent layers is (1 - 2z)c. As seen in the figure, the titanium atoms lie at the center of a polyhedron whose vertices are defined by the six nearest neighbor sulfur atoms. The titanium site symmetry in this polyhedron would be octahedral (point group Oh) if the axial ratio were “ideal” (c/a

Periodic Hartree-Fock Study of TiS2

Figure 1. (A) Projection of one S-Ti-S layer of TiS2 onto the ab plane. Small spheres are coplanar Ti atoms and large + (-) spheres are sulfur atoms above (below) the Ti plane. Hexagonal lattice vectors for the primitive unit cell are a, b, and c, where |a| ) |b|, a ⊥ c, b ⊥ c, and ∠(a,b) ) 120°. Perpendicular vectors a′, b′, and c′ (c′ ) c) define a supercell for β′ ) ∠(a′,c′) ) 90°. Cartesian axes are x and y and x′ and y′ for the hexagonal lattice and the supercell, respectively. (B) TiS2 unit cell depicting a rigid-layer shear in the a′ direction corresponding to strain 13. Note that this unit cell belongs to a trigonal space group (P3hm1) for β′ ) 90° and to a monoclinic space group (C2/m) for arbitrary β′.

) x(8/3) ) 1.6330) and if the sulfur planes were equidistant (z ) 1/4). However, the titanium site symmetry actually corresponds to an elongated triangular antiprism (point group D3d) since the observed c/a ratio is larger than the ideal value. The lattice vectors of a monoclinic unit cell (Z ) 2, space group C2/m) are also shown in Figure 1A. Lattice vectors a′, b′, and c′ are shown for the special case where the monoclinic and hexagonal lattices are equivalent (i.e. a′ ) x3b′ and ∠(a′,c′) ) β′ ) 90°). In this case, the two coordinate systems are related by a′ ) a - b, b′ ) a + b, and c′ ) c. The fractional atomic positions in the nonprimitive system are Ti(0, 0, 0), Ti(1/2, 1/2, 0), S(1/3, 0, z), S(2/3, 0, -z), S(5/6, 1/2, z), and S(1/6, 1/2, -z). This unit cell is used to compute the elastic constant corresponding to shear parallel to the ab plane, which is depicted in Figure 1B. 4. Computational Results Our results are presented in three sections, which are entitled Structural Properties, Zero-Pressure Electronic Properties, and Nonmetal-to-Metal Phase Transition. 4.1. Structural Properties. In this section, we first define the interlayer binding energy, ILBE, and the interlayer potential energy, ILP(a,c,z). We then demonstrate the validity of the RLapproximation whereby ILP(a,c,z) at various c and ILBE are obtained. The remainder of this section will focus on the zeropressure lattice parameters and atomic coordinates, bulk modulus, linear moduli, elastic constants, and cohesive energy. 4.1.1. Interlayer Binding Energy. The interlayer binding energy is defined as ILBE ) ILP(a,c,z)L - ILP(a*,c*,z*), where ILP(a,c,z)L and ILP(a*,c*,z*) are the interlayer potential energies

J. Phys. Chem., Vol. 100, No. 39, 1996 15737

Figure 2. TiS2 interlayer potential as a function of c-axis length using the RL-approximation. Solid lines are fits of HF+CC data to modified Morse potentials: E(c) ) E0 + De[e-b(c-c0) - 1]2 + d(c - c0′), where E0, De, b, c0, d, and c0′ are constants. The inset diagram depicts the rigid-layer expansion of TiS2 layers, where S- and S+ denote the sulfur atomic planes in Figure 1A. The zero of energy is the sum of atomic energies (E(Ti) + 2E(S)). For clarity, the plots were set equal at c ) 10 Å.

of isolated layers and zero-pressure TiS2, respectively, per formula unit. In practice, ILP(a,c,z)L is assumed to be the energy that no longer changes with further increases in the c-axis; that is, ILP(a,c,z)L ≈ ILP(a,large c,z). Clearly one can compute ILP(a,c,z) by applying the said “Coptimization” procedure at successively larger values of c. In the case of TiS2, the resultant geometries satisfy the criteria of the RL-approximation, since 2cz differed from 2c*z* by less than 0.41% for 5.40 Å e c e 8.00 Å.31 As a result, considerable computational economies can be made without sacrificing overall accuracy. Using the RL-approximation we obtain ILP(c) by calculating the total energy at many c while requiring a ) a* and z ) c*z*/c. The HF, HF+CS, HF+P, and HF+WL ILP(c) are shown in Figure 2. The minima of these curves are the ILP(c*), and the well depths are the ILBE ) ILP(large c) - ILP(c*). The large differences between the HF and HF+CC ILBE dramatically illustrate the importance of (and uncertainties in) correlation corrections in TiS2. As seen in Figure 2, marked differences in the ILP(c) and ILBE exist depending upon the specific method used. Perhaps most importantly, the PHF calculations alone (without correlation corrections) exhibit no significant energy minimum and therefore predict that TiS2 is an unstable compound. In contrast, the HF+CS curve shows a minimum at c ) 6.3 Å. However, this value of c is 11% larger than the experimental (T ) 300 K) value of 5.6953 Å.32 The predicted c-axis improves using the HF+P method to c ) 6.0 Å (5.4% larger than experiment). The best results are obtained using the HF+WL method, which binds the crystal at c ) 5.616 Å (-1.4% different from experiment). The corresponding HF+WL interlayer binding energy (ILBE ) 3.8 kcal/mol) is a factor of ∼3 larger than

15738 J. Phys. Chem., Vol. 100, No. 39, 1996

Clerc et al.

TABLE 2: Resultsa HF

HF+P

HF+WL

experiment

3.574b N/A 0.2346b N/A 2.952b

3.417(0) 5.990(9) 0.2281(20) 2.732(24) 3.07(1)

3.321(1) 5.616(5) 0.2440(27) 2.740(30) 3.47(1)

3.3977(2)c 5.6629(2)c 0.2501(4) 2.832(4) 3.2854(4)c

∆E°coh,0K (kcal/mol) 282(5)b 398(5)

474(5)

340(1)

B (GPa) Ba (GPa) Bc (GPa)

34(3) 314(9) 44(5)

42(4) 311(21) 58(6)

41d 179,e 260f 58,e 61f

271(2) 6(1) 42(4) 0(2)

263(3) 9(3) 55(5) 15(5)

202g 14g 55(5)h 18(4)h

a (Å) c (Å) z 2cz (Å) F(g/cm3)

c11 + c12 (GPa) c13 (GPa) c33 (GPa) c44 (GPa) s11 + s12 (GPa-1) s13 (GPa-1) s33 (GPa-1)

251b

0.00371(3) 0.00385(4) 0.00513g 0.00054(9) -0.00064(19) -0.0013g 0.024(2) 0.019(1) 0.019g

a See ref 31 for general computational details. See Appendix for cohesive energy calculation. b Estimated by taking c ) 5.6953 Å. c 5 K value computed from room temperature values: a ) 3.4073(2) Å and c ) 5.6953(2) Å32 and 5-300 K thermal expansivity Ra ) 0.96 × 10-5 K-1 and Rc ) 1.94 × 10-5 K-1.34 d Estimated from ref 35, Figure 2. e From ref 36. f From ref 35. g Calculated from the experimental values Ba ) 260 GPa, Bc ) 61 GPa, and c33 ) 55 GPa. h From ref 37.

those obtained using the HF+CS and HF+P methods (ILBE ) 1.3 kcal/mol). However, in each case the small ILBE is indicative of weak binding. Evidently, correlation corrections to PHF theory are required to produce interlayer binding in TiS2. A similar result was previously noted by Harrison33 for the layered crystal MgCl2. However, in that case, empirical van der Waals pair potentials were used to stabilize the layered structure. 4.1.2. Lattice Parameters. In this section we discuss the calculated zero-pressure (equilibrium) geometry of TiS2. Predicted values of lattice parameters, internal coordinate, and density (F) are reported in Table 2. Since the HF method alone predicted ILBE ≈ 0 kcal/mol, only HF+CC geometries will be discussed. The data in Table 2 demonstrate that the predicted values of the a- and c-axes depend strongly on the specific correlation functional used. For example, the differences between the calculated and experimental values of the a-axis are 0.6% and 2.2% for HF+P and HF+WL, respectively. In contrast, the corresponding c-axis differences are 5.8% and -0.8%. In contrast to the above result, the accuracy of the predicted layer thickness is essentially independent of the functional used. This is illustrated by the small difference, 0.3%, between the HF+CC values of 2cz. 4.1.3. Elastic Properties. Calculated bulk modulus (B), linear moduli (Ba and Bc), stiffness constants (cij), and compliance constants (sij) as well as their experimental counterparts (where available) are listed in Table 2. As previously noted, the calculated elastic constants include the effects of internal relaxation, since U-optimizations are used in their computation. Details of these calculations and their associated uncertainties are given ref 31. In the following discussion, we use elastic properties to interpret the crystal bonding strength (stiffness) in various directions. For example, the stiffnesses in the b c direction, in directions parallel to the ab plane, and in linear combinations of these directions are c33, c11 + c12, and c13, respectively. In addition, the stiffness associated with shearing (13) of adjacent layers is c44. Many of the HF+WL elastic constants are in remarkably good agreement with the experimental values. For example,

the differences between the calculated and experimental B, Bc, c13, c33, and c44 are 2.4%, 3.3%, 36%, 0%, and 17%, respectively. (Note: the sizeable differences in c13 and c44 relative to experiment are mainly due to the small size of the numbers involved.) The largest discrepancy between the HF+WL results and experiment occurs in c11 + c12 and Ba. In these cases, calculated values of c11 + c12 and Ba are 30% and 20% (using ref 35) larger than experiment, respectively. A similar overestimation occurred using the HF+P approach. We find that relative to experiment, the HF+P elastic constants and moduli that strongly depend upon bonding in the c direction are underestimated, whereas those more dependent upon bonding in the ab plane are slightly overestimated. For example, c33, c44, B, and Bc are underestimated by 24%, 100%, 17%, and 27%, respectively, as compared to experiment. In contrast, c11 + c12 and Ba are overestimated by 34% and 21% (using ref 35), respectively. The magnitudes and trends among the TiS2 elastic constants are consistent with highly anisotropic bonding. For example, the value of c11 + c12 is greater than that of c33 by a factor of ∼5. This reflects the stronger intralayer bonding relative to the much weaker interlayer forces. The coupling between the layers is weak, as evidenced by the relatively small magnitude of c13. The small size of c44 reflects the relative ease with which the layers may be sheared or slid across one another. We find that the TiS2 layers remain essentially rigid during small shears in the a′ direction. This shear corresponds to strain 13 (see Figure 1B). Specifically, minimization of the HF+CC total energy with respect to the monoclinic unit cell internal coordinates at several a′, c′, and β′ results in the same intralayer geometry. (See ref 31 for details.) We note that the large difference (81 GPa) between the experimental values of Ba35,36 indicates that the HF+CC values of c11 + c12 and Ba may be more accurate than the said comparisons with experiment suggest. The large difference results mainly from the difficulty of growing single crystals large enough for accurate measurements. This has been a noted experimental limitation in the study of TiS2. Providing an accurate description of “nonbonded” interlayer forces is a stringent test for density functional correlation corrections used in conjunction with a HF Hamiltonian. On the basis of the predicted lattice parameters and elastic constants, we prefer the Wigner-Levy functional for the study of TiS2. However, more general statements concerning the overall applicability the WL functional are not warranted until additional systems have been tested. 4.1.4. CohesiVe Energy. The cohesive energy of TiS2 is defined as ∆Ecoh ) (ETi + 2ES) - E(a*,c*,z*), where ETi + 2ES is the total energy of the isolated atomic constituents. Calculated values of the cohesive energy at T ) 0 K (∆Ecoh,0K) are shown in Table 2. Details of the calculations appear in the Appendix. To provide crystal and atomic energies that differ by approximately the same amount from the exact HF limits, we follow the convention of using the same basis set shell structure for the crystal and atomic calculations.16 Specifically, only the diffuse Gaussian exponents of the basis set in Table 1 were reoptimized in the atomic calculation. The resultant total atomic energies were then correlation-corrected and used to compute ∆Ecoh,0K. See the Appendix for further details. As seen in Table 2, ∆Ecoh,0K is underestimated using the HF approach and is overestimated using both HF+CC methods. The HF value (282 kcal/mol) is 17% below the experimental value. In contrast, the values obtained using HF+P (398

Periodic Hartree-Fock Study of TiS2 kcal/mol) and HF+WL (474 kcal/mol) are greater by 41% and 39%, respectively. The underestimation of the cohesive energy using the HF approach was also observed in a recent pseudopotential CRYSTAL92 study of 17 III-V and IV-IV semiconductors.42 Errors ranging from -30% to -50% in the HF results, due to the lack of electron correlation, were found to be qualitatively similar to molecular results. However, using the P functional improved the agreement to within 4.2%. In another recent CRYSTAL92 study, the HF cohesive energies of Al2O3, SiO2, and the Li, Na, K, Be, Mg, Ca, and Mn oxides were found to be 20-50% smaller than the experimental values.43 The cohesive energies obtained after using the P and CS correlation corrections were improved, but still below experiment by 0-20%. In contrast to our results for TiS2, the HF+CC cohesive energies were not overestimated in either of the two cited studies. It is unlikely that this inconsistency can be attributed to the high anisotropy of TiS2, since the amount of additional binding in the c direction achieved using the correlation corrections (1-4 kcal/mol) is small relative to the total binding energy (280-480 kcal/mol). More data are required to find whether these results are characteristic of layered compounds. 4.2. Zero-Pressure Electronic Properties. In this section, the electronic structure of zero-pressure TiS2 is studied using computed Mulliken populations, band structures, total densities of states (DOS), and projected densities of states (PDOS). The energy bands and the DOS are compared with experimental spectra from angle-resolved photoelectron spectroscopy (ARPS), X-ray photoelectron spectroscopy (XPS), and ultraviolet photoelectron spectroscopy (UVPS). In the discussions to follow, the Ti(3d) and S(3p) orbital labels refer to the Table 1 notations Ti(1d+2d) and S(3p+4p), respectively. 4.2.1. Mulliken Population. Mulliken population analyses are used in this study as rough indicators of the partitioning of charge among the atoms in TiS2 and the charge rearrangements that take place when TiS2 is formed from the atoms. The Mulliken orbital populations of atomic titanium (3F), atomic sulfur (3P), and TiS2 were computed using Molpro MCSCF and CI, respectively.44-49 To accurately model the Ti atomic charge, a diffuse Ti(4s) orbital (exponent ) 0.05) was added to the basis set in Table 1. For further details, see ref 31. On the basis of our Mulliken populations, we find that 1.9 electrons are transferred from Ti(4s) to S(3p) when TiS2 is formed from the elements. Consequently, the net atomic charges in TiS2 are Ti+1.90S2-0.95. These oxidation states are similar to but slightly more ionic than the results of Zunger and Freeman (Z&F),10 which are Ti+1.19S2-0.60. Our results indicate that approximately two electrons occupy the Ti(3d) orbitals. 4.2.2. Electronic Band Structure. The first Brillouin zone for the hexagonal lattice of TiS2 is shown in Figure 3 after the notation of Bradley and Cracknell.50 The set of k points (BZ segments) used in computing the energy bands are also indicated along with their associated point groups in Scho¨nflies notation. In the discussions to follow, numeric band labels refer to the eigenvalue order rather than the symmetry. For example, band 27 is always the highest energy occupied eigenstate since there are 54 electrons per TiS2 formula unit. In this section we also state the orbital derivations of various energy bands as derived from the PDOS (discussed in the next section). Lastly, the zero of energy is at the highest energy occupied eigenstate and is called the “Fermi level” in this paper. The TiS2 energy bands are shown in Figure 4A. Evidently the bands are grouped into several distinct sets. The six highest

J. Phys. Chem., Vol. 100, No. 39, 1996 15739

Figure 3. TiS2 first Brillouin zone showing Cartesian axes (i, j, and k), direct lattice vectors (a, b, and c), reciprocal lattice vectors (g1, g2, and g3), and several point group elements: m ) mirror plane, C2 ) 2-fold rotation. The k point reciprocal lattice coordinates and their associated point groups (in Scho¨nflies notation) are also shown. The perimeter of the shaded planes corresponds to the segments given in the electronic band structure (Figure 4).

energy valence bands (VB1) form a group derived from S(3p) orbitals with a contribution from Ti(3d) orbitals occurring at ca. -5 eV. Below VB1 there are two valence bands (VB2), which have a predominant S(3s) character. Directly above VB1 we have a set of five conduction bands (CB1), which are derived from Ti(3d) orbitals. The minimum in CB1 occurs at L, and a local minimum 0.3 eV higher in energy occurs at M. The topology and relative dispersions of our energy bands are in good qualitative agreement with the density functional results of Z&F10 and the Green’s function (KKR) results of Myron and Freeman (M&F).12 For example, in each of the cited studies the locations of the TiS2 valence-to-conduction gaps are the indirect transition Γ f L (Eu f Ag), the indirect transition Γ f M (Eu f Ag), and the direct transition Γ f Γ (Eu f Eg), in order of increasing energy. However, our energy gaps and bandwidths are larger, which is a result typical of HartreeFock. In addition, there occur slight differences in the ordering of the bands. For example, our VB1 bands 26 and 27 at Γ (Eu) lie at EF, whereas band 25 (A2u) is ∼0.4 eV below EF. This order is the same as M&F but is reversed with respect to Z&F. This particular band order affects the TiS2 nonmetal-metal phase transition mechanism, as will be discussed in a later section. The large difference between the calculated optical band gap (OBG), which is 6.42 eV (Γ f L), and the experimental value of 0.18 ( 0.06 eV9 is typical of HF calculations of nonmetallic materials. This error results from the well-known fact that the differences between the eigenvalues of the one-electron Fock Hamiltonian are a poor approximation to single-particle excitation energies. Nevertheless, the toplogy of the occupied manifold and the first conduction bands in Hartree-Fock band structures is usually correct.16 The calculated bands are compared to an ARPS (He I 21.2 eV) spectrum51 of the Γ-M-A-L plane in Figure 5. The zero of energy in the theoretical and ARPS bands is the top of VB1 and the experimentally determined highest energy occupied

15740 J. Phys. Chem., Vol. 100, No. 39, 1996

Clerc et al.

Figure 4. (A) TiS2 Hartree-Fock electronic energy bands, (B) total density of states, and (C) Al KR (1486.6 eV) XPS spectrum (dashed line below EF),53 He II (40.8 eV) UVPS spectrum (solid line below EF),52 and He II (40.8 eV) IPES spectrum (solid line above EF).57 The labels VB2, VB1, and CB1 correspond to the groups of bands centered at -18, -5, and +10 eV, respectively. The VB1 DOS features are labeled following ref 52. The bands are labeled with the irreducible representation of the corresponding point group of k in Scho¨nflies notation (see Figure 3).

state, respectively. Because the component of the photoelectron wavevector perpendicular to the crystal surface (kz) is not conserved during transport through the surface,56 each ARPS band is an average of a given band over kz. To compare the ARPS data to our calculations, VB1 bands on parallel segments between Γ-M and A-L were calculated in order to match each band on Γ-M with the symmetryequivalent band on A-L. Specifically, VB1 bands were computed on the seven BZ segments parallel to Γ-M whose endpoints are the coordinates (0, 0, mπ/6c) - (2π/(x3a), 0, mπ/6c), where m ) 0, 1, ..., 6. The evolution of each band was followed as a function of m, and the symmetry-equivalent band pairs were then marked by connecting them with slanted lines. The bands at m ) 0 and m ) 6 were found to form envelopes for their respective bands at arbitrary m. Each ARPS band, therefore, corresponds (approximately) to the geometric average of a given band on Γ-M and A-L. For example, the set of band 22 positions (i.e. the lowest energy VB1 band) for arbitrary m lies between the positions on Γ-M (solid line) and A-L (dashed line). This set lies within the area denoted by the corresponding slanted lines. The center of gravity of this area is the predicted location of the corresponding ARPS band. As seen in Figure 5, the ARPS band directly below EF coincides well with the geometric center of the two highest energy VB1 bands (symmetry labels A′ and A′′). In contrast, the ARPS band centered at -3.8 eV differs in position by ∼1 eV relative to the nearest calculated bands. However, as deduced from comparing our VB1 DOS to photoelectron spectra (discussion following), this ARPS band actually corresponds to the flat calculated band (A′′) centered at -5 eV. Similarly,

the ARPS band at -5 eV corresponds to the parts of two calculated bands (both A′) which are centered at -7.1 eV. The ARPS data shown in Figure 5 indicates a conduction band near M and L. The investigators attributed this band to emission from the lowest energy conduction band (titanium d-based) by electrons derived from excess titanium, displaced titanium, or both.51 To demonstrate that the BZ location of this conduction band coincides with the calculated minimum in CB1, band 28 was rigidly shifted down on the Γ-M and A-L segments. The magnitude of this shift (6.24 eV) was such that the experimental zero-pressure optical band gap (Γ f L, 0.18 eV) was reproduced. As seen in Figure 5, the topology of the shifted CB1 band agrees well with the ARPS spectra, although its position is slightly high by ∼0.5 eV. 4.2.3. Density of States. The discussion of the DOS includes comparison to photoelectron spectra and it is divided into two parts. The first part presents the general features of the DOS and PDOS such as bandwidths, energy gaps, and orbital assignments. The second part is a detailed discussion of the VB1 DOS. In all DOS and PDOS, the zero of energy is the highest energy VB1 eigenstate. 4.2.3.1. General Features. The TiS2 DOS analysis is summarized in Table 3, where the theoretical VB2, VB1, and CB1 bandwidths, energy gaps, and VB1 peak separations are tabulated. The corresponding experimental values and the density functional results of Z&F10 are also included for comparison. The calculated DOS and the photoelectron spectra are shown in Figure 4B,C. Figure 4C shows XPS (Al KR: 1486.6 eV53), UVPS (He II: 40.8 eV52), and IPES (He II: 40.8 eV57) spectra

Periodic Hartree-Fock Study of TiS2

J. Phys. Chem., Vol. 100, No. 39, 1996 15741

Figure 5. TiS2 Hartree-Fock VB1 energy bands on the BZ segments Γ-M (solid lines) and A-L (dashed lines). He I (21.2 eV) ARPS data51 for the Γ-M-L-A plane is superimposed on the theoretical bands. Triangles represent data collected at oblique (20°) incidence. Open circles represent data having the wave vector in the first Brillouin zone. Filled circles represent data that were folded back from the second zone in the reduced zone scheme (after ref 51). Symmetry labels correspond to the irreducible representations of D3d (Γ and A), Cs (∑ and R), and C2h (M and L) in Scho¨nflies notation. Slanted lines connect a given band on Γ-M to the corresponding band on A-L. The lowest energy conduction band was shifted toward EF by 6.24 eV to set the OBG equal to 0.18 eV.9

TABLE 3: Theoretical DOS and Experimental Photoelectron Data theory (eV) peak a peak R peak β

-0.5 -3.6 -4.9

peak γ peak δ R-β β-γ γ-δ VB1 width VB2 width CB1 width VB1-VB2 gap CB1-VB1 gap

-7.1 -8.5 1.3 2.2 1.4 8.8 (5.5e) 3.2 (1.9e) 6.3 (3.6e) 8.7 (6.8e) 6.5 (0.3e)

assignment S(3px,y,z) S(3px,y,z) S(3px,y,z), Ti(3dxz,yz,x2-y2,xy) S(3px,y,z) S(3px,y,z)

experiment (eV) -0.5,a n.o.b -2.1,a -2.7b -3.9,a -3.5b -5.3,a -4.8b n.o.,a n.o.b 1.8,a 0.8b 1.4,a 1.3b 6-7,a 4-4.5,b 6-7c 2-3,a 1-2b 6-7,b 7-8d 0.18(6),f 0.3(2)g

UVPS (He II 40.8 eV).52 b XPS (Al KR 1486.6 eV);53 n.o. ) not observed. c Appearance potential.54 d X ray absorption.55 e Density functional calculation.10 f From pressure dependence of the Hall coefficient.9 g ARPS (He I 21.2 eV).51 a

of TiS2, where the features in VB1 are labeled (a, R, β, γ, and δ) following the notation of Shepard and Williams.52 The DOS is shown in Figure 4B, where the VB1 features are similarly labeled including subscripts (ac, Rc, βc, γc, and δc). Comparison of the energy bands (Figure 4A) to the DOS (Figure 4B) shows that the more prominent DOS features arise from several conspicuous regions of low band dispersion. VB1 peak ac derives from the nearly flat bands 26 and 27 (Eu) nearest EF on the BZ segment Γ-A and the flat portion of band 25 (A2u) near Γ on the Γ-M and K-Γ segments. The large DOS labeled βc is produced from the flat VB1 band located at -5

Figure 6. TiS2 projected densities of states. The dotted line corresponds to the Fermi energy. VB2, VB1, and CB1 are as defined in Figure 4. VB1 features are labeled after ref 52.

eV throughout most of the BZ. Similary, the large VB2 DOS peak derives from the nondisperse band centered at ca. -18 eV. There are two conspicuous differences between the photoelectron spectra (Figure 4C) and the DOS (Figure 4B), both of which are typical of HF theory without correlation corrections to the eigenvalues. The most apparent difference is the overestimation of the energy gaps. For example, the calculated VB1-VB2 gap is overestimated by 1.7-2.7 eV relative to the XPS spectrum. The second difference is the overestimation of the bandwidths. In particular, the bandwidths of VB1 and VB2 are overestimated by ∼4.5 and ∼1.7 eV, respectively, relative to the XPS results. When the effects of electron correlation are included, the bandwidths and energy gaps dramatically decrease. This is evident in Table 3, where our results are compared to the density functional results of Z&F.10 For example, the Z&F bandwidths of VB1 and VB2 (5.5 and 1.9 eV, respectively) are most consistent with the corresponding XPS results (4-4.5 and 1-2 eV). The orbital derivations of VB2, VB1, and CB1 can be inferred from the PDOS shown in Figure 6. Evidently VB1 derives mainly from S(3p) orbitals except for peak βc, where there is a Ti(3d) contribution. VB2 derives from S(3s) orbitals, and the lower energy states in CB1 are assigned to Ti(3d) orbitals. The five regions of high VB1 DOS derive from S(3px,y,z) orbitals, and peak βc contains an additional contribution from Ti(3dxz,yz,x2-y2,xy) orbitals. 4.2.3.2. VB1 Features. The VB1 DOS is compared to an XPS spectrum (Al KR, 1487 eV)53 in Figure 7A, where the DOS is convoluted with Gaussians having widths (hwhm) corresponding to typical instrument resolution. This convolution includes the effects of instrument resolution (0.6-0.8 eV58) as

15742 J. Phys. Chem., Vol. 100, No. 39, 1996

Figure 7. (A) TiS2 VB1 DOS convolution (hwhm ) 0.6 eV) scaled by calculated atomic X-ray cross sections (solid line) and XPS (Al KR, 1486.6 eV53) spectrum (triangles). (B) TiS2 VB1 DOS convolution (hwhm ) 0.1 eV) (solid line) and UVPS (He II, 40.8 eV52) spectrum (circles). In both figures, experimental peaks are labeled following ref 52. Calculated peaks are similarly labeled including subscript “c”.

well as differing X-ray absorption cross sections. Specifically, the total titanium and sulfur PDOS were scaled by the calculated atomic X-ray total cross sections (Ti, 169 700 b (barns); S, 45 290 b59). The resultant DOS was then convoluted with Gaussians (hwhm ) 0.6 eV). As seen in Figure 7A, the VB1 DOS features are very similar to those of the XPS spectrum. Both consist of a strong central peak, βc, and two smaller peaks, Rc and γc, where peak Rc occurs as a shoulder. However, the experimental peaks lie closer to EF than those calculated. The differences between peaks ac, Rc, βc, and γc and their experimental counterparts are 0.9, 1.4, and 2.3 eV, respectively (see Table 3). Here, the differences increase with distance from EF. There are marked differences between the VB1 DOSs shown in Figures 7A and 4B. Most conspicuously, there is a large accentuation of peak βc in Figure 7A due to the relatively high Ti(3d) content at that energy and the relatively large titanium X-ray cross section. In addition, peak broadening obscures peaks ac and δc in Figure 7A. In Figure 7B, the VB1 DOS is compared to a UVPS spectrum.52 The Gaussian convolution (hwhm ) 0.1 eV) corresponds to the reported instrument resolution ( 84 GPa . pexp) crystal geometry. The value pcalc ) 84 GPa is obtained by noting the pressure at which the OBG vanishes. This result is an order of magnitude greater than pexp ) 4.0 ( 0.5 GPa.9 However, the bulk of this large error is due to the considerable overestimation of OBG0,calc because PSM,calc and PSM,exp are similar. 4.3.3. Pressure Dependence of Structural Properties. The calculated c(p) data at low pressures (0 e p e 5.5 GPa) are compared to experiment in Figure 10A. The theoretical and experimental values agree quite well (within 0.8%) and both decrease linearly with pressure. A linear fit of the predicted data for 0 e p e 5.5 GPa gives c(p) ) 5.6653-0.06195p which yields c(pexp) ) 5.42(3) Å at the NM-M transition. The calculated a(p) data are similarly shown in Figure 10B. The difference between theoretical and experimental values (e2.5%) is reasonably small, although larger than for c(p). Linear fitting of the calculated data gives a(p) ) 3.31930.00873p and a(pexp) ) 3.284(5) Å. A similar fit of the experimental data results in aexp(p) ) 3.4095-0.0115p and the extrapolated aexp(pexp) ) 3.363(6) Å. An analogous plot of z versus pressure yields z(pexp) ) 0.2567(16).31

15744 J. Phys. Chem., Vol. 100, No. 39, 1996

Clerc et al.

band 25 [Γ-A]:

Figure 10. Calculated (open circles) and experimental (filled circles)35 pressure dependence of the TiS2 (A) c-axis and (B) a-axis. Solid lines are linear fits of the data, which yield c(pexp) ) 5.42(3) Å, a(pexp) ) 3.284(5) Å, and aexp(pexp) ) 3.363(6) Å. The vertical line is at pexp.

Changes in the intralayer dimensions are small at p ) pexp, and this further supports the results obtained using the RLapproximation (recall Figure 2).31 The layer thickness (2cz) has increased by only 0.6%. In addition, the a-axis has decreased by only 1.2%, and the Ti-S bond length (d(Ti-S) ) (c2z2 + a2/3)1/2) has decreased by only 0.5%. The largest pressure-induced structural change in TiS2 at pexp is the 9.0% decrease in the distance across the van der Waals gap (c-2cz).31 This agrees with the earlier finding that c is the “softest” direction in the crystal (c33 , c11 + c12). The planes of sulfur atoms that define the van der Waals gap are consequently moved closer together. As discussed in the next section, the proximity of these sulfur atoms plays an important role in the NM-M phase transition in TiS2. 4.3.4. Pressure Dependence of the Energy Bands. In this section, we identify the energy band that is most responsive to pressure and discuss the origin of this sensitivity. We then discuss the predicted behavior of this band during the increase in pressure p f pexp f pcalc. Following this, the VB1 band involved in the VB1-CB1 band overlap at pexp is predicted and compared to that inferred from the value of (mh/me)exp. This comparison will show the two results to be inconsistent. However, this inconsistency will be shown to result from a small error in the relative positions of the two highest energy zeropressure VB1 bands at Γ. In the following discussion, we refer to several VB1 crystal orbitals whose band structure at zero pressure is shown in Figure 4. These eigenvalues and the corresponding (unnormalized) crystal orbitals with k-dependence (eik‚R) omitted are

band 28 [M,L]: Ag(k), φAg(r,k) ) Ti(3dz2) + Ti(3dxz) + Ti(3dx2-y2) band 26-27[Γ-A]: Eu(k), φ1Eu(r,k) ) S1(3px) + S2(3px) and φ2Eu(r,k) ) S1(3py) + S2(3py)

A2u(k), φA2u(r,k) ) S1(3pz) + S2(3pz)

The bands and corresponding crystal orbitals are defined on the BZ segments and points given in square brackets. S1 and S2 are the sulfur atoms at (1/3, 2/3, z) and (2/3, 1/3, -z) in Figure 1, respectively. The forms of the crystal orbitals are found to be valid for 0 e p e 90 GPa by computing the eigenvectors at the volumes shown in Figure 8. Hence, the crystal orbital corresponding to Ag(k) (the lowest energy CB1 band) has the same composition in both the NM and M phases. The eigenstates most sensitive to pressure are A2u(Γ-A). The source of this unusual pressure sensitivity is the strong dependence of φA2u(r,Γ-A) on the overlap between sulfur atoms on opposite sides of the VDW gap. At k ) Γ, this overlap is antibonding, while at k ) A it is bonding. The dispersion of A2u(k) on Γ-A (A2u(Γ) - A2u(A)) increases with pressure because the said sulfur atoms are moved closer together. The pressure sensitivity of A2u(k) also contributes to the strong VB1 bandwidth pressure dependence (recall Figure 9). We observe that Eu(Γ) increases with pressure on Γ-A, albeit not nearly as rapidly as A2u(Γ). At sufficient pressure, A2u(Γ) surpasses Eu(Γ). At very high pressure A2u(Γ) . Eu(Γ). For example, at p ) 93 GPa, A2u(Γ) - Eu(Γ) ) 4.1 eV. The dispersion of Eu(k) on Γ-A is not strongly affected by the decrease in the van der Waals gap, since φEu(r,k) is localized within the layers. Even at very high pressure (90 GPa) the dispersion of Eu(k) is small (2.7 eV) relative to that of A2u(k) (9.9 eV). The highest energy occupied VB1 eigenstate at pexp is Eu(Γ), and it is calculated from the pressure, pcross, at which A2u(Γ) ) Eu(Γ). The value of pcross is found by fitting the energies of bands 25 and 26-27 at Γ to linear functions of pressure over the range 0 e p e 10.0 GPa. The resulting equations are A2u(Γ) ) 0.1869p - 0.4392 and Eu(Γ) ) 0.1171p + 0.0399, where the units of energy and pressure are electronvolts and gigapascals, respectively. Setting A2u(Γ) ) Eu(Γ) yields pcross ) 6.8 ( 1.8 GPa. Since pexp < pcross and since at zero pressure Eu(Γ) > A2u(Γ), then the highest energy occupied VB1 eigenstate at pexp is Eu(Γ). Alternately, the experimental result (mh/me)exp g 1 can be used to identify the highest energy VB1 band at pexp or equivalently, the VB1 band involved in the phase transition. Specifically, the effective mass ratio can be related to our calculations by comparing the curvature of bands 25-28 nearest EF on the BZ segments Γ-M-K-Γ (corresponding to the ab plane). This comparison shows that the band curvatures increase in the order A2u(Γ) < Ag(M) e Eu(Γ). Therefore, the corresponding effective masses increase in the order mhEu e meAg < mhA2u, where we assumed that the carriers in A2u(Γ) and Eu(Γ) are holes and that the carriers in Ag(M) are electrons. Therefore, if the highest energy occupied state at pexp is Eu(Γ), then (mh/me)calc e 1. However, if this state is A2u(Γ), then (mh/ me)calc > 1, which agrees with experiment. Therefore, the experimental effective mass ratio suggests that A2u(Γ) is the highest energy VB1 band at pexp, on the basis of our theoretical zero-pressure energy bands. The same result follows for 0 e p e 10 GPa since the said band curvatures are nearly invariant under low pressure. To summarize the previous arguments, the result pexp < pcross indicates that the VB1 band involved in the phase transition is Eu(k). In contrast, the experimental effective mass ratio and band curvatures are consistent with A2u(k). Comparison of our bands to those of Z&F10 indicates that the discrepancy pexp < pcross may be due to incorrect positions of A2u(Γ) and Eu(Γ) at p ) 0. Specifically, our band structure has A2u(Γ) < Eu(Γ), whereas the density functional results of

Periodic Hartree-Fock Study of TiS2

J. Phys. Chem., Vol. 100, No. 39, 1996 15745 The results of two empirical corrections, shown in Figure 11B, indicate that 1.5GPa e p′calc e 4.0 GPa. These two corrections correspond to different OBG0,exp obtained from Hall coefficient data (0.18 ( 0.06 eV9) and from ARPS (0.30 ( 0.10 eV64). As shown in Figure 11B, VB1(hi) and the two corrected CB1(lo) (labeled CB1(lo)′a and CB1(lo)′b) are plotted at low pressure. The correction that yields CB1(lo)′a is a rigid shift of CB1(lo) toward VB1(hi) equal to OBG0,calc - OBG0,exp ) 6.42 eV - 0.18 eV ) 6.24 eV. Similarly, the second correction, which results in CB1(lo)′b, corresponds to OBG0,exp ) 0.30 eV.64 The crossing point of VB1(hi) and the two corrected CB1(lo) thus occur within the interval 1.5 GPa e p′calc e 4.0 GPa, which is consistent with pexp ) 4.0 ( 0.5 GPa. 5. Conclusions

Figure 11. (A) TiS2 VB2 (diamonds), VB1 (circles), and CB1 (triangles) extrema as functions of pressure. The minima and maxima are labeled “lo” and “hi”, respectively. (B) Effects of empirical correction: CB1(lo) has been rigidly shifted toward VB1(max) by 6.24 eV (filled triangles: CB1(lo)′a) and 6.12 eV (open triangles: CB1(lo)′b) to set OBG0,calc equal to the experimental values of 0.18 ( 0.069 and 0.30 ( 0.10 eV,64 respectively.

Z&F10 have A2u(Γ) > Eu(Γ). If their ordering is correct, then the highest energy VB1 band is A2u(k) irrespective of pressure. 4.3.5. Empirical Rigid-Band Correlation Correction. In this section, we empirically correct pcalc by rigidly shifting CB1 toward VB1 at each pressure such that OBG0,calc ) OBG0,exp. This correction is justified by previously published Green’s function many-body correlation corrections to HF energy bands.67,68 In particular, the screened-exchange-plus-Coulombhole (SECH) method has been shown to affect the HF band structures of diamond67 and LiF68 much like a “scissor operator” which rigidly lowers conduction bands and raises valence bands. The effects of the SECH correction are recounted here for the case of diamond.67 For diamond, the lowest energy conduction band is shifted down by 3.61-4.29 eV, where the larger corrections occur for the eigenstates furthest from EF. Similarly, the highest energy valence band is shifted up by 2.45-3.68 eV. Thus the corrected bands tend to flatten out somewhat while essentially retaining their topology. The SECH correction to OBG0,calc is 6.5 eV, which yields OBG0,calc ) 5.6 eV, in excellent agreement with OBG0,exp ) 5.5-5.6 eV. Note that the correction required for TiS2, OBG0,calc - OBG0,exp ) 6.24 eV, is close to that derived for diamond. As was already discussed, the result pcalc . pexp is largely due to the result OBG0,calc . OBG0,exp. It will now be shown that empirically correcting OBG0,calc results in a corrected pcalc (p′calc), which is consistent with pexp. The effect of the empirical correction is illustrated in Figure 11A,B. The pressure dependences of VB2, VB1, and CB1 before the correction are shown in Figure 11A. The highest and lowest energy states of each group are labeled “hi” and “lo”, respectively. The most responsive states in the occupied manifold are the VB1 states near VB1(hi) due to the pressure sensitivity of A2u(k) for k near Γ, as previously discussed.

The object of the present study is to gain insight into the electronic properties of TiS2 by an ab initio theoretical investigation. The results of these calculations indicate the following. (1) The HF independent-electron approximation is not adequate for describing the van der Waals type bonding in TiS2. Post-SCF correlation corrections to the PHF total energies are required to bind the S-Ti-S layers with respect to displacements in the c direction. Two correlation-only density functional corrections were used to calculate structural properties: Perdew 91+ and Wigner-Levy. The results derived using the Perdew 91+ functional agree reasonably well with experiment. However, the results obtained using the Wigner-Levy functional agree most favorably. (2) The elastic constants indicate that the intralayer bonding is much stronger than the interlayer bonding and that the coupling of S-Ti-S sandwiches in the c direction is relatively weak. (3) On the basis of Mulliken populations, the formation of crystalline TiS2 from the atoms involves the transfer of ∼2e of titanium charge from the Ti(4s) orbital to the S(3p) orbitals. (4) The energy bands are in good agreement with previously published results aside from the overestimation of the energy gaps and bandwidths (which is typical of Hartree-Fock calculations). The location of the optical band gapsthe indirect transition Γ f Lsagrees with ARPS experiments and with previously published density functional results. (5) Projected densities of states indicate that the set of six highest energy valence bands are composed predominantly of S(3p) orbitals with a contribution from Ti(3d) orbitals occurring at peak βc. Except for the overestimation of the bandwidths and the VB1 peak separations, the features present in UVPS and XPS spectra are qualitatively well-reproduced by our calculations. In particular (a) the small peak a at the Fermi level is shown to be of S(3p) origin and not derived from the 3d states of excess titanium atoms, and (b) the strong peak β is found to contain the largest contribution from Ti(3d) orbitals, in agreement with photoelectron observations. (6) The pressure response of the optical band gap in both phases is predicted to be linear. Furthermore, the pressure response of the M phase (30 meV/GPa) is reasonably consistent with the experimental value (45 ( 15 meV/GPa), but is smaller than that of the NM phase at low pressure (86 meV/GPa). (7) Applying two different empirical rigid-band shifts of the conduction bands such that experimental zero-pressure optical band gaps are reproduced (0.18 ( 0.06 and 0.30 ( 0.10 eV) results in p′calc between 1.5 and 4.0 GPa, which is consistent with pexp. The energy bands involved in the band overlap at pexp are predicted to be A2u(Γ) and Ag(L). Lattice parameters

15746 J. Phys. Chem., Vol. 100, No. 39, 1996

Clerc et al.

a, c, and z are predicted to be 3.284(5) Å, 5.42(3) Å, and 0.2567(16), respectively, at pexp. Acknowledgment. The authors are grateful for helpful discussions with J. E. Jaffe. This work was supported by the Associated Western University, Inc. Northwest (AWU NW) Laboratory Graduate fellowship, and the research was carried out at Pacific Northwest National Laboratory. Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract No. DEAC06-76RLO#1830. Appendix The theoretical cohesive energy at 0 K, ∆E°coh,0K ) E(Ti) + 2E(S) - E(TiS2), was computed from the following data: E(Ti)1 (au) E(S)1 (au) E(Ti) + 2E(S) (au) E(TiS2) (au) ∆E°coh,0K (au) ∆E°coh,0K (kcal/mol)

HF -848.138 -397.489 -1643.115 -1643.5652 0.450 282

HF+P -849.074 -398.136 -1645.362 -1645.997 0.635 398

HF+WL -849.073 -398.147 -1645.376 -1646.132 0.756 474

1 Atomic HF and HF+CC energies are from CRYSTAL calculations using the basis set in Table 1 with the titanium 4sp, 5sp, 1d, 2d, and sulfur 3sp, 4sp, and 1d exponents optimized in the atom. Configurations: 3F (Ti) and 3P (S). 2 Taken as the HF energy occurring at the HF+WL geometry.

The experimental value of ∆E°coh,0K was computed from the following Born-Haber cycle:

(1)

TiS2,crys(0K) f TiS2,crys(298K)

E1 ) H°f,298[TiS2,crys] - H°f,0[TiS2,crys] ) 12.378 kJ/mol38,39 TiS2,crys(298K) f Ticrys(298K) + 2Srhom(298K)

(2)

E2 ) -∆H°f,298[TiS2,crys] ) 410.3 ( 2.4 kJ/mol38 (3)

Ticrys(298K) f Tigas(298K) E3 ) ∆H°f,298[Tigas] ) 469.21 ( 2.1 kJ/mol40

(4)

Tigas(298K) f Tigas(0K) E4 ) H°f,0[Tigas] - H°f,298[Tigas] ) -7.540 kJ/mol40

(5)

2Srhom(298K) f 2Sgas(298K) E5 ) -2∆H°f,298[Sgas] ) 2(278.805 kJ/mol)40

(6)

2Sgas(298K) f 2Sgas(0K) E6 ) 2[H°f,0[Sgas]] - H°f,298[Sgas]] ) 2(-6.657 kJ/mol)40

(Sum)

TiS2c(0K) f Tigas(0K) + 2Sgas(0K) ∆E°coh,0K ) 1428.64 kJ/mol ) 341.4 kcal/mol

The zero-point correction to the cohesive energy (Ezp) was found to be 0.5 kcal/mol by using the Debye temperature θD ) 234 K41 in Ezp ) 9RθD/8, where R ) 8.314 J‚mol-1‚K-1. References and Notes (1) Julien, C. Microionics 1991, 309. (2) Mao, Z.; White, R. E. J. Power Sources 1993, 43 (1-3), 181. (3) Shen, D. H.; Halpert, G. J. Power Sources 1993, 43 (1-3), 327. (4) Baresel, D.; Sarholz, W.; Scharner, P.; Schmitz, J. Ber. BunsenGes. Phys. Chem. 1974, 78 (6), 608. (5) Atlantic Richfield Co. Neth. Appl. 75 02,156 (Cl. C07C, B01J), 30 May 1975, U.S. Appl. 459998, 11 Apr 1974.

(6) Price, J. A.; Stewart, M. J. U.S. 3644291 (Cl. 260-75R; C08g), 22 Feb 1972, Appl. 777460, 20 Nov 1968. (7) Golub, A. S.; Protsenko, G. A.; Kozinkin, A. V.; Tatevosyan, M. M.; Shuvaev, A. T.; Novikov, Yu. N. IzV. Akad. Nauk SSSR, Ser. Khim. 1992, 5, 44. (8) Farr, J. P. G.; Franks, A. J.; Rowe, G. W. Lubr. Eng. 1980, 36, 586. (9) Klipstein, P. C.; Friend, R. H. J. Phys. C: Solid State Phys. 1984, 17, 2713. (10) Zunger, A.; Freeman, A. J. Phys. ReV. B 1977, 16 (2), 906. (11) Bullett, D. W. J. Phys. C: Solid State Phys. 1978, 11, 4501. (12) Myron, H. W.; Freeman, A. J. Phys. ReV. B 1974, 9 (2), 481. (13) Umrigar, C.; Ellis, D. E.; Wang, D.-s.; Krakauer, H.; Posternak, M. Phys. ReV. B 1982, 26 (9), 4935. (14) von Boehm, J.; Isoma¨ki, H. J. Phys. C: Solid State Phys. 1980, 13, 3181. (15) Dovesi, R.; Saunders, V. R.; Roetti, C. CRYSTAL92 An ab-initio Hartree-Fock LCAO program for periodic systems; 1992. (16) Pisani, C.; Dovesi, R.; Roetti, C. Hartree-Fock Ab Initio Treatment of Crystalline Systems Springer-Verlag: Berlin, 1988. (17) Causa`, M.; Dovesi, R.; Pisani, C.; Roetti, C. Phys. ReV. B 1986, 33, 1308. (18) Dovesi, R.; Ermondi, C.; Ferrero, E.; Pisani, C.; Roetti, C. Phys. ReV. B 1984, 29, 3591. (19) Causa`, M.; Dovesi, R.; Roetti, C.; Kotomin, E.; Saunders, V. R. Chem. Phys. Lett. 1987, 140, 120. (20) Dovesi, R.; Pisani, C.; Roetti, C. J. Chem. Phys. 1987, 86, 6967. (21) Catti, M.; Pavese, A.; Dovesi, R.; Roetti, C.; Causa`, M. Phys. ReV. B 1991, 44, 3509. (22) Hess, A. C.; Saunders, V. R. J. Phys. Chem. 1992, 96, 4367. (23) Dovesi, R.; Causa`, M.; Orlando, R.; Roetti, C.; Saunders, V. R. J. Chem. Phys. 1990, 92, 7402. (24) McCarthy, M. I.; Hess, A. C. J. Chem. Phys. 1992, 96, 6010. (25) Saunders, V. R. Private communication. (26) Harrison, N. M. Private communication. (27) Wilson, L. C.; Le´vy, M. Phys. ReV. B 1990, 41 (18), 12930. (28) Perdew, J. P.; Wang, Y. Phys. ReV. B 1991, 45, 13244. (29) Causa`, M.; Dovesi, R.; Pisani, C.; Colle, R.; Fortunelli, A. Phys. ReV. B 1987, 36 (2), 891. (30) Colle, R.; Salvetti, O. J. Chem. Phys. 1990, 93 (1), 534. (31) Clerc, D. G.; Ph.D. Thesis - Washington State University; May, 1995. (32) Chianelli, R. R.; Scanlon, J. C.; Thompson, A. H. Mat. Res. Bull. 1975, 10, 1379. (33) Harrison, N. M.; Saunders, V. R. J. Phys.: Condens. Matter 1992, 4, 3873. (34) Subba Rao, G. V., Shafer, M. W., Le´vy, F., Eds. Physics and Chemistry of Materials with Layered Structures, Vol. 6: Intercalated Layered Materials: Intercalation in Layered Transition Metal Dichalcogenides; Reidel Publishing Co.: Dordrecht, Holland, 1979; pp 123-137. (35) Negishi, H.; Ohara, S.; Inoue, M. Phys. Status Solidi B 1989, 151, 441. (36) Zhou, O.; Fischer, J. E.; Liang, K. S. Phys. ReV. B 1991, 44 (14), 7243. (37) Scha¨rli, M.; Le´vy, F. Phys. ReV. B 1986, 33 (6), 4317. (38) O’Hare, P. A. G.; Johnson, G. K. J. Chem. Thermodyn. 1986, 18, 189. (39) Beyer, R. P. J. Chem. Eng. Data 1983, 26, 347. (40) Weast, R. C., Ed.; Astle, Melvin J., Assoc. Ed. CRC Handbook of Chemistry and Physics, 63rd ed.; Chemical Rubber Publishing Co.: Boca Raton, FL, 1982. (41) Benda, J. A. Phys. ReV. B 1974, 10, 1409. (42) Causa`, M.; Dovesi, R.; Roetti, C. Phys. ReV. B 1991, 43 (14), 11937. (43) Dovesi, R.; Roetti, C.; Freyria-Fava, C.; Apra`, E.; Saunders, V. R.; Harrison, N. M. Phil. Trans. R. Soc. London A 1992, 341, 203. (44) MOLPRO is a package of ab initio programs written by H.-J. Werner and P. J. Knowles, with contributions from J. Almlof, R. D. Amos, M. J. O. Deegan, S. T. Elbert, C. Hampel, W. Meyer, K. A. Peterson, R. Pitzer, A. J. Stone, and P. R. Taylor. (45) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1985, 82, 5053. (46) Knowles, P. J.; Werner, H.-J. Chem. Phys. Lett. 1985, 115, 259. (47) Werner, H.-J.; Meyer, W. J. Chem. Phys. 1980, 73, 2342. (48) Werner, H.-J.; Meyer, W. J. Chem. Phys. 1981, 74, 5794. (49) Werner, H.-J. AdV. Chem. Phys. 1987, LXIX, 1. (50) Bradley, C. J.; Cracknell, A. P. The Mathematical Theory of Symmetry in Solids; Clarendon Press: Oxford, 1972; p 107. (51) Chen, C. H.; Fabian, W.; Brown, F. C.; Woo, K. C.; Davies, B.; DeLong, B.; Thompson, A. H. Phys. ReV. B 1980, 21 (2), 615. (52) Shephard, F. R.; Williams, P. M. J. Phys. C: Solid State Phys. 1974, 7, 4416. (53) Wertheim, G. K.; DiSalvo, F. J.; Buchanan, D. N. E. Solid State Commun. 1973, 13, 1225. (54) Webb, C.; Williams, P. M. Phys. ReV. B 1975, 11, 2082. (55) Fischer, D. W. Phys. ReV. B 1972, 8, 3576.

Periodic Hartree-Fock Study of TiS2 (56) Margaritondo, G.; Grasso, V., Ed. Physics and Chemistry of Materials with Low-Dimensional Structures: Series A: Layered Structures: Electronic Structure and Electronic Transitions in Layered Materials; D. Reidel Publishing Co.: Dordrecht, Holland, 1986; p 405. (57) Nametame, H.; Nohara, S.; Negishi, H.; Inoue, M.; Suga, S. Surf. Sci. 1991, 242, 294. (58) Fujimori, A.; Suga, S.; Negishi, H.; Inoue, M. Phys. ReV. B 1988, 38 (6), 3676. (59) Scofield, J. H. J. Electron Spectrosc. Relat. Phenom. 1976, 8, 129. (60) Williams, P. M., Lee, P. A., Eds. Physics and Chemistry of Materials with Layered Structures, Vol. 4: Optical and Electrical Properties: Photoemission Studies of Materials with Layered Structures; D. Reidel Publishing Co.: Dordrecht-Holland/Boston, 1976; pp 274, 292, 295.

J. Phys. Chem., Vol. 100, No. 39, 1996 15747 (61) Scha¨rli, M.; Brunner, J.; Vaterlaus, H. P.; Le´vy, F. J. Phys. C: Solid State Phys. 1983, 16, 1527. (62) Shephard, F. R.; Williams, P. M. J. Phys. C: Solid State Phys. 1974, 7, 4427. (63) Williams, P. M.; Shephard, F. R. J. Phys. C: Solid State Phys. 1973, 6, L36. (64) Barry, J. J.; Hughes, H. P.; Klipstein, P. C.; Friend, R. H. J. Phys. C: Solid State Phys. 1983, 16, 393. (65) Clerc, D. G. Unpublished. (66) Murnagahan, F. D. Am. J. Math. 1937, 59, 235. (67) Brener, N. E. Phys. ReV. B 1975, 11 (2), 929. (68) Brener, N. E. Phys. ReV. B 1975, 11 (4), 1600.

JP952264R