Quantum Chemistry of Solids / 22-LCAO Calculations of Perfect-crystal Properties

9.2 Electron Properties of Crystals in LCAO Methods

387

applied, [373], which uses in the exchange part the mixture of the Fock (20%) and Becke’s (80%) exchange, whereas in the correlation part the Perdew–Wang (PWGGA) nonlocal correlation functional is employed. As was demonstrated in [606], the B3PW functional gives the best description of the atomic and electronic structures, as well as elastic properties of several ABO3 perovskite materials. In some cases, for a comparison the B3LYP hybrid functional was also used (the exchange part is the same as in the B3PW functional and the correlation is the Lee–Yang–Parr nonlocal functional). For comparison, the DFT PW calculations with the PW91 GGA nonlocal functional were made in [639]. As was noted earlier DFT PW calculations with hybrid functionals are not implemented in modern computer codes.

The relative smallness of the total energy di erences between di erent magnetic phases requires a high numerical accuracy in the direct-lattice and Brillouin-zone summations. Therefore, in [639] the cuto threshold parameters of the CRYSTAL code for Coulomb and exchange integrals evaluation have been set to 7, 7, 7, 7, and 14, respectively. The integration BZ has been carried out on the Monkhorst–Pack grid of shrinking factor 8 (its increase up to 16 gave only a small change in the total energy per unit cell).

The self-consistent procedure was considered as converged when the total energy in the two successive steps di ers by less than 10−6 a.u. In order to check how the results depend on the basis-set (BS) choice, B3LYP and B3PW LCAO calculations were performed for the cubic LaMnO3 with one formula unit per primitive cell with the experimental lattice constant of 3.95 ˚A. Five di erent BS (BS1-BS5) were used, including the BS optimized both in all electron (AE) and pseudopotential calculations (largeand small-core Hay–Wadt pseudopotentials were taken for La and Mn atoms). The details of BS optimization and choice are given in [639]. The total energies for the cubic structure (with the five atom primitive unit cell) obtained in nonspin-polarized and spin-polarized B3PW calculations were compared with di erent magnetic ordering of four d electrons on the Mn3+ ion (the d-levels are split in the cubic field into t2g and eg levels, the latter are higher in energy than the former): total spin projection Sz = 2 (four electrons occupy t2g and eg levels), Sz = 1 (three α electrons and one β electron occupy the t2g level), and Sz = 0 (two α electrons and two β electrons occupy the t2g level). The calculated self-consistent Mn atom magnetic moment is close to 0, 2 and 4 for the three cases under consideration. Since the number of electrons per cell depends on the basis used (AE or without core electrons), the absolute values of these energies are quite di erent. The basis optimization for the same BS choice results in a lower energy for the same spin projection. It was shown that the order of relative energies for di erent magnetic configurations is the same for all BS chosen and even the absolute values of energy di erences are close. Independently of the BS choice the lowest total energy per cell corresponds to the maximal spin projection Sz = 2, i.e. for the Mn3+ ion in a crystal the Hund rule holds and the lowest energy corresponds to the maximal spin. This result demonstrates that the magnetic behavior of insulators can usually be discussed in terms of the physics that applies to isolated atoms or ions that have a magnetic moment. The Mulliken atomic charges remain practically the same for di erent magnetic orderings, provided the BS is fixed. At the same time, their absolute values show the BS dependence for the same spin projection value. The relative energies of di erent magnetic orderings reveal the same sign in PW DFT and LCAO calculations (the lowest energy corresponds to Sz = 2). The Bader atomic

388 9 LCAO Calculations of Perfect-crystal Properties

charges, calculated in the PW DFT method, compared with the Mulliken charges, are much smaller and also show a weak dependence on the magnetic ordering.

To study the magnetic ordering in LaMnO3, the so-called broken symmetry approach was adopted [638]. It allows one to deduce the magnitude of the magnetic coupling data making spin-polarized calculations for di erent magnetic orderings of transition-metal atoms. LaMnO3 is stabilized at moderate temperatures in the orthorhombic structure comprising four formula units, space group D216h (in P bnm and P nma settings the largest orthorhombic lattice translation vector is directed along the z or y axis, respectively; in what follows, the P bnm setting is chosen). The real

structure can be viewed as a distorted cubic perovskite structure with a quadrupled

√ √

tetragonal unit cell (ap 2, ap 2, 2ap) where ap is the lattice parameter of the cubic perovskite structure. This tetragonal cell is slightly distorted in the real orthorhombic structure. The calculations on the tetragonal model (both in LCAO and PW basis), i.e. four formula units without the structural distortion have shown that LaMnO3 is metallic in all magnetic states and the ground state is FM. This contradicts the experimental data, for both the energy gap (LaMnO3 is believed to be a spin-controlled Mott–Hubbard insulator with the lowest-energy d–d transitions around 2 eV [640]) and magnetic AAF ordering in the ground state below 140 K. The results obtained in [639] for cubic LaMnO3 with the primitive cell of five atoms explain this fact: the tetragonal structure remains in fact cubic, with the tetragonal supercell; since in the primitive unit cell the FM configuration corresponds to the metallic ground state, the same is true for the undistorted tetragonal structure. However, when the orthorhombic atomic distortions are taken into account, the AAF structure turns out to be the ground state, in agreement with experiment.

This is seen from results given in Table 9.23, [639], where the relative energies of di erent magnetic phases are presented (the energy for the FM phase is taken as zero energy). The calculations were made for orthorhombic structure with the structural parameters from a neutron-di raction study [641]. Magnetic moments µ on the Mn atom and magnetic and coupling constants Jab and Jc are also presented there.

The Ising model Hamiltonian

H = −Jab SziSzj − Jc	SzkSzl	(9.79)
ij	kl

was used, where Jab and Jc are exchange integrals (magnetic coupling constants) between nearest neighbors in the basal plane (xy) and between nearest neighbors along the c-axis, respectively, Szi stands for the z-component of the total spin on the magnetic center i, and (ij) and (kl) indicate summation over intraplane and interplane nearest magnetic centers, respectively. Due to two di erent possible choices of the Ising Hamiltonian presentation, we stress that in (9.79) that the form of (9.78) is used that contains double summation over each pair of centers and gives positive values for Jab and negative for Jc. The latter must be taken into account in a comparison with the experimental data obtained by fitting to the Ising Hamiltonian. In particular, experimental Jab and Jc values [641] have to be multiplied by a factor of 2 for the Hamiltonian (9.79).

A set of equations relating the energy di erences for the FM, AAF, GAF and CAF configurations with the magnetic coupling constants sought for was used:

E(FM) − E(AAF) = E(CAF) − E(GAF) = −32Jc

(9.80)

9.2 Electron Properties of Crystals in LCAO Methods

389

Table 9.23. The energy (in meV per unit cell) of the di erent magnetic phases for orthorhombic LaMnO3 for the experimental structure [641]. The energy of the FM configuration is taken as zero energy. Magnetic moments µ on the Mn atom in µB , magnetic coupling constants Jab and Jc in meV. Experimental data: µ =3.87 for AAF [632], Jc=–1.2, Jab=1.6 [640]

Method/basis set	AAF		GAF		CAF		Jc	Jab
UHF(BS1)a	–4.8	—	55.6	—	—	—	–0.15	0.94
B3LYP(BS1)b	–33.0	3.80	112.0	3.72	120.0	3.73	–0.64	2.07
UHF(BS2)c	–8.0	—	48.0	—	56.0	—	–0.25	0.88
B3LYP(BS2)b	–33.0	3.78	84.0	3.71	101.0	3.71	–0.78	1.70
UHF(BS3)d	–5.2	3.96	51.2	—	55.9	—	–0.15	0.88
Fock-50(BS3)d	–12.2	3.89	89.22	—	93.64	—	–0.26	1.52
B3LYP(BS3)d	–32.2	3.80	114.0	—	121.5	—	–0.62	2.09
B3LYP(BS4)b	–32.0	3.81	103.0	3.73	117.0	3.76	–0.72	1.97
B3LYP(BS5)b	–30.0	3.82	106.0	3.74	117.0	3.77	–0.64	1.98
B3PW(BS5)b	–19.0	3.86	153.0	3.75	152.0	3.79	–0.28	2.50
GGAPW91(BS5)b	–40.0	3.71	248.0	3.58	234.0	3.64	–0.83	4.08
(GGA-PW)b	–59.0	3.55	189.0	3.35	224.0	3.46	–1.47	3.69
(LAPW)e	–72.0	—	96.0	—	136.0	—	–1.75	2.38
(LMTO)f	–62.0	3.46	243.0	3.21	—	—	–1.94	4.76
FLMTO(GGA)g	–98.8	—	142.8	—	167.2	—	–1.92	3.19
(LDA+U)h	–34.0	—	170.0	—	—	—	–1.06	2.66

a Reference [635], b Reference [639], c Reference [634], d Reference [638], e Reference [642], f Reference [643], g Reference [644], h Reference [645].

E(FM) − E(CAF) = E(AAF) − E(GAF) = −64Jab

(9.81)

Again, to avoid misunderstanding, (9.80) and (9.81) are written for a quadruple cell, corresponding to four formula units. For calculating the coupling constants both (9.80)and (9.81) were used and an averaging was performed. Unfortunately, E(CAF) is not always presented in the published results (see Table 9.23). The calculated magnetic coupling constants are compared with the experimental data in the last two columns of Table 9.23. This table demonstrates how the calculated magnetic coupling constants depend on the Hamiltonian choice (LDA, UHF, hybrid) and the BS choice (LCAO, PW, LAPW).

These results can be briefly summarized as follows. Independently from the BS and Hamiltonian used, all calculations mentioned in Table 9.23 correctly reproduce the sign of the experimental exchange integrals and their relative values (|Jab| > |Jc|). In all UHF and hybrid (B3LYP, B3PW) calculations the exchange integrals agree better with the experimental data than in DFT calculations. We explain this by the incorporation of the Fock exchange into UHF and the hybrid methods. For example, if we fix the AO basis as BS3 and analyze a series of the UHF (pure Fock exchange), Fock-50 (50% of Fock exchange) and B3LYP (20% of Fock exchange included), the coupling constants Jab and Jc in this series are getting closer to the experimental values. The e ect of the correlation part is smaller (compare B3LYP and B3PW results for BS5). The lack of Fock exchange in LCAO GGA and GGA-PW calculations leads to overestimated values of both magnetic coupling constants. This overestimate is

390 9 LCAO Calculations of Perfect-crystal Properties

well observed also in previous LAPW, FLMTO, and LDA+U calculations. Therefore, in agreement with the conclusions, made by other researchers it is demonstrated that when calculating the experimentally observable magnetic coupling constants, the nonlocal exchange plays an important role. Such hybrid DFT or UHF calculations are practically possible only for the LCAO BS. These results confirm again the conclusion [632] that the CRYSTAL code is a valuable tool for the study of magnetic properties for open-shell transition-metal compounds. The optical gap for the AAF orthorhombic phase was also calculated in [639]. The B3PW LCAO gives 2.9 eV and 4 eV for the Mn d-d and O2p-Mnd transitions, in good agreement with the experiment [640], whereas the DFT PW gap of 0.6 eV is an underestimate typical for the DFT. These results, in line with those in [638], demonstrate that it is a combination of the nonlocal exchange with correlation e ects realized in the hybrid functionals that considerably narrows the gap between calculated and experimental magnetic coupling constants. This questions the generally accepted idea that DFT is better suited than UHF and related methods for the study of transition-metal oxides, in particular, LaMnO3. The surface calculations of this crystal made in [639] we discuss in Chap. 11.

We noted that the hexagonal manganites RMnO3 (R = Sc, Y, Ho–Lu) are an interesting group of compounds because of their unusual combination of electric and magnetic properties. At low temperatures, they show coexistence of ferroelectric and magnetic orderings. ScMnO3 plays a prominent role in the series of hexagonal manganites. It has the smallest distance between magnetic Mn3+ centers along the c-direction and has the highest Neel temperature. Below the temperature 1220 K its structure belongs to space group P 63cm [646]. The Mn3+ ions are surrounded by five oxygen ions that form a distorted trigonal bipyramid, see Fig. 2.18. In this environment the 5D state of Mn3+ with 3d4 occupancy splits into three states. The magnetic moment of the Mn3+ ions in ScMnO3 has been measured as 4.11 µB . Its Neel temperature is 130 K [647], much higher than that of YMnO3 (71 K), [648]. The assignment of the strong optical absorption band at about 1.6 eV for hexagonal manganites is a matter of debate [649]. From local spin-density approximation (LSDA) and LSDA+U calculations of YMnO3 [650] this band was interpreted as an interatomic charge-transfer transition from O 2p to Mn 3d orbitals. In another theoretical study of YMnO3 based on LSDA and LSDA+U [651], it was found that the states near the top of the valence band have predominantly Mn 3dx2−y2 character, while the lowest-unoccupied subband has Mn 3dz2 character. The d-band centers of the occupied and empty bands are separated by 2.0 eV. The latter theoretical result is in agreement with a careful analysis of experimental spectra [649] that revealed that the band at 1.6 eV is due to an intra-atomic 3d − 3d transition in the Mn3+ ions. The O(2p)→ Mn(3d) charge-transfer transition in LaMnO3 and related compounds occurs at 3.1 eV or higher. Later experiments on ScMnO3 [652] were interpreted di erently, partly based on previous LSDA calculations. In particular, the strong band in ScMnO3 (at 1.57 eV), ErMnO3 (1.59 eV), and YMnO3 (1.61 eV) was considered now as a charge-transfer band. In this context we have to point out that LSDA values for bandgaps are usually much too small. In [638] a charge-transfer transition in LaMnO3 was found at 2.3 eV with B3LYP and at 4.0 eV with Fock35 (35 per cents of the HF exchange in the hybrid functional) functionals. Considering that one-particle energies obtained by band-structure calculations do not describe local d–d transitions, the agreement between LSDA charge-transfer and the experimentally determined band may be fortuitous.

9.2 Electron Properties of Crystals in LCAO Methods

391

The LCAO calculations [653] with the computer code CRYSTAL03, [23] is the first attempt to establish the electronic and magnetic structure of the low-temperature phase of ScMnO3 from first principles. The solid phase was modeled with threedimensional periodic supercells of stoichiometry Sc6Mn6O18. Three di erent Hamiltonians were used: UHF, hybrid HF-DFT (B3LYP) and DFT (BLYP). In a previous study of magnetic coupling in CuGeO3 [632] the same Hamiltonians had been compared and it was found that the hybrid HF-DFT B3LYP Hamiltonian gives the most accurate values for the magnetic coupling constant J. Optimized basis functions were used for Mn and O as obtained in [639]. The valence basis of Sc was optimized for Sc2O3 bulk crystal, Sc atom inner (1s 2sp) electrons were described with the Hay– Wadt small-core e ective core potential. The Sc6Mn6O18 cell contains two hexagonal atomic planes. Three magnetic phases were considered: ferromagnetic (FM) when the four unpaired electrons of each of six Mn ions have parallel spins; two antiferromagnetic (AFc and AFa). In the AFc phase the spins of the Mn ions in di erent atomic planes along the c-direction are antiparallel, while all spins are parallel in the same plane (between two planes, the Mn–Mn distance is about 5.8 ˚A, while it is only 3.4 ˚A in the plane). In the AFa state each plane contains two types of Mn with antiparallel spin in a ratio of 1:2. The local spin arrangements are shown in Fig. 9.11.

DE

Fig. 9.11. Two local spin configurations in the AFa state in a hexagonal plane; in case (a) the central spin-down Mn ion is surrounded by 4 spin-down and 2 spin-up Mn ions, in (b) the central spin-up Mn ion is surrounded by 6 spin-down Mn ions.

This is a simplification of the real situation in ScMnO3 where the Mn spins are not collinear [647]: the projections of neighboring Mn spins in the same plane are in fact rotated by 120◦.

If the FM state and the AFc state are compared, the two neighbors along the c-direction in the planes above and below each magnetic center have di erent spin. It was found that the corresponding interplane coupling constant Jc ≤ 0.2 meV for all three Hamiltonians under consideration, i.e.the magnetic coupling between Mn ions of di erent planes is negligible. The results of the calculation of intraplane coupling constant are given in Table 9.24. The magnetic coupling constants of ScMnO3 have

392 9 LCAO Calculations of Perfect-crystal Properties

Table 9.24. Energy di erence (eV) between the AF and the FM states of ScMnO3, estimated magnetic coupling constant (meV), Mulliken charges q (a.u.), magnetic moment of Mn ions (µB ) obtained with three di erent Hamiltonians [653]

Hamiltonian		Ja	Qsc	QMn	QO	e
UHFa	–0.163	4.1	+2.3	+2.4	–1.5	4.0
B3LYPa	–0.555	13.9	+1.9	+2.0	–1.37	3.7
B3LYPb	–0.449	11.2	+1.9	+2.0	–1.3	3.8
BLYPa	–0.925	23.1	+1.87	+1.9	–1.2	3.5
Exper.		5.5–7.0c				4.11d

a with fixed structural parameters from experiment [654] b with fully optimized in [653] structures

c estimated from a comparison with YMnO3 (the Neel temperature 70 K and Ja =3.8 meV) d Reference [648]

not yet been measured at the time of calculations [653]. The simplest assumption is a linear relation between J and the Neel temperature TN [655]. It is then possible to estimate Ja of ScMnO3 from the coupling constant of the isostructural compound YMnO3 (3.0–3.8 meV, [655]). With TN = 130 K for ScMnO3 and TN =71 K for YMnO3 [648] one obtains Ja =5.5–7.0 meV. In Table 9.24 the calculated coupling constants are compared to this crude estimate. The UHF method gives Ja =4.1 meV, which is slightly smaller than the estimated value. When electron correlation is taken into account with B3LYP, a much larger value (Ja =13.9 meV) is obtained. But also the description of electron exchange has a substantial e ect on the magnetic structure. The BLYP method, which has no contributions from exact HF exchange, gives the largest coupling constant, Ja =23.1 meV. The same ordering of magnetic coupling constants obtained with these methods has been obtained in studies of CuGeO3 [632], where the B3LYP result favorably agreed with the measured value. Since for ScMnO3 the deviation from the estimated crude value is substantial with B3LYP, it is likely that the linear extrapolation scheme is inadequate. But before any conclusion was drawn the e ect of changes in the local structure on the calculated coupling constant was investigated. A full optimization of hexagonal lattice parameters and fractional coordinates of ScMnO3 was performed with B3LYP for both the FM and the AFa states, [653]. The energy di erence between these two states is reduced by 20% (compared with that obtained for the experimental structure) and Ja accordingly decreases to 11.2 meV. Since this value is still considerably larger than the estimated one the above statement about the inaccuracy of the linear interpolation scheme holds. In Table 9.24 there are also compared the calculated self-consistent magnetic moments µef f with the experimental value 4.11 µB [648]. The obtained values range from 3.5 µB (BLYP) to 4.0 µB (UHF) and are thus in agreement with experiment. The small change from 3.7 µB to 3.8 µB due to optimization at the B3LYP level indicates that the local atomic structure has a small e ect on this property. Atomic charges were calculated based on the Mulliken population analysis. Table 9.24 shows that the ionicity is reduced in the series from UHF to BLYP. The oxygen and manganese charges range from –1.2 to –1.5 and +1.9 to +2.4, respectively, indicating substantial covalent contributions to the bonding. These values do not essentially di er from those (–1.6 to –1.8 and +1.5 to +1.8), obtained in the LaMnO3 calculations [639]. The

9.3 Total Energy and Related Observables in LCAO Methods for Solids

393

PDOS calculated for ScMnO3 with B3LYP for the fully optimized structure of the AFa state shows the following. The highest-occupied states between –0.5 eV and 0 eV (the Fermi level is taken as the zero energy), and the lowest-unoccupied bands between 3.3 eV and 4.5 eV consist of a mixture of O(2p) and Mn(3d) orbitals. Sc orbitals only contribute to higher unoccupied states above 5 eV. In agreement with experimental estimates [649], the interband transition at 3.3 eV with B3LYP was found. On the basis of the calculated PDOS the excitation is neither a pure d–d nor a charge-transfer transition, rather it represents an intermediate situation.

As is seen from Table 9.24 the structure-parameter optimization influences the results of the electronic-structure calculations. Such an optimization also allows the calculation of the such important properties of solids as the cohesive energy, relative stability of phases with di erent structure, bulk moduli, elastic constants. LCAO calculations of these total energy related properties are considered in the next section.

9.3 Total Energy and Related Observables in LCAO Methods for Solids

9.3.1 Equilibrium Structure and Cohesive Energy

The theoretical determination of the equilibrium structure of molecules and periodic systems and their cohesive energy is of primary importance in quantum chemistry of molecules and solids.

In the molecular systems, the parameters that must be optimized are the nuclear coordinates. Such an optimization requires the numerical or analytic calculation of forces – the total energy derivatives with respect to nuclear coordinates. The methodology for calculating analytic first derivatives at the HF (SCF) LCAO level of theory was suggested by Pulay [656]. Extensions of the HF derivatives to density-functional theory methods were straightforward. Derivatives of correlated (post-HF) energies followed a number of years after the SCF derivatives [657]. Analytical gradients have become a standard tool for the optimization of molecular structures and are implemented in modern molecular electronic-structure codes. Calculations with the structure optimization are now possible on molecules with hundreds of atoms, and energy derivatives provide an enormous amount of information about the potential-energy surface at very little additional cost. Analytic derivatives with respect to geometrical parameters are used not only for the structure optimization but also for exploring potential-energy surfaces and the study of transition states and following reaction paths.

Such molecular properties as vibrational frequencies, IR and Raman intensities, NMR shielding constants, etc., can be formulated in terms of second and higher derivatives with respect to geometry and applied fields. Such calculations are now practical and routine using analytic derivatives at the SCF level and a few correlated methods, [657].

In periodic systems, the cell dimensions are a set of optimized structure parameters additional to nuclear coordinates in the primitive unit cell. Nowadays, the majority of solid-state codes compute the total energy with KS PW methods. In KS PW calculations the analytical gradients of total energy for nuclear coordinates and cellparameter optimization are implemented in computer codes and widely used in the

394 9 LCAO Calculations of Perfect-crystal Properties

electronic-structure calculations of crystals [377]. The more complicated task is the analytical gradients theory in LCAO methods for periodic systems. The total energy E of the crystal (per primitive cell) is the sum of electron energy Ee and nuclear repulsion energy EN : E = Ee + EN . The electron energy Ee as calculated within the HF LCAO approximation can be expressed in terms of the one-electron density matrix and includes the direct lattice sums, see (7.46). In the KS LCAO method only the electron-density is involved in the Hamiltonian, see (7.47).

The analytical gradients in HF and MP2 LCAO calculations were applied for a long time only for one-periodic (1D) systems [658, 659]. The e cient method of analytic energy-gradient calculation in the KS LCAO method for periodic systems with one-, twoand three-dimensional periodicity was presented only recently in [82] and implemented in the code GAUSSIAN03 [107]. A real space approach is used, where all the summations are performed in direct space. Energy derivatives with respect to geometrical parameters (atomic forces and stress tensor) are computed analytically to high accuracy using techniques based on the fast multipole method (see Chap. 7). This allows for the accurate lattice parameters and nuclear-coordinate optimization and evaluation of vibrational frequencies using finite di erences of analytic forces. The extension of this technique to the hybrid HF-DFT LCAO Hamiltonians allowed optimization of the structure of di erent solids in good agreement with experiment (as was demonstrated in Chap. 7).

Simultaneously, progress has been achieved to implement analytical HF gradients in the LCAO code CRYSTAL [23, 660]. The formulas for analytic total-energy gradients with respect to nuclear coordinates in periodic systems were presented in [661, 662] and illustrated in test calculations. It was demonstrated that an e - cient nuclear-coordinates optimization of large systems with any periodicity can be performed. Later, the theory of analytical HF gradients with respect to the cell parameter for three-periodic systems was presented [663], so that full structure optimization with the help of analytical gradients is now possible. It is important to note that the CRYSTAL code is based on the Ewald method in three dimensions, so that computing analytical gradients with respect to the cell parameter requires various additional derivatives that were not yet available with the implementation of nuclear gradients, and this has been documented in great detail in [663]. The oneand two-dimensional cases are again di erent because di erent potentials are used in a real-space approach. The article [664] complements [663] by analytical HF gradients with respect to cell parameters for systems periodic in one and two dimensions. The implementations mentioned include the cases of spin-restricted and unrestricted polarization.

We intentionally do not give here the mathematical details of the approaches allowing implementation of analytical gradients in the calculations of periodic systems (for these details readers are referred to the cited publications and references therein). We note that such an implementation is essentially more complicated than in the case of molecular systems and requires high accuracy in the total-energy calculation. One can find the detailed analysis of the accuracy in gradients calculations on numerical examples in [660,661,663]. The comparison of the numerical and analytical derivatives values can also be found. In what follows we turn to the results obtained for the equilibrium structure and cohesion energy in the crystalline metal oxides. The LCAO calculations discussed were made with the CRYSTAL code and use of HF, KS and hybrid Hamiltonians.

9.3 Total Energy and Related Observables in LCAO Methods for Solids

395

The cohesive (binding) energy Ecoh is defined as the energy required to dissociate a solid into noninteracting atomic species A with energies EA

Ecoh = EA − Esolid	(9.82)
A

In this definition, Ecoh is positive for any thermodynamically stable crystal. The summation is made over the atoms included in the primitive unit cell, Esolid is the total energy per primitive cell. Calculated values of the cohesive energy are compared with experimental results that can be obtained by measuring the latent heat of sublimation at various low temperatures, and extrapolating to zero Kelvin.

The binding-energy definition depends on what the constituent parts are considered to be: reference can be made to the ions in ionic crystals or to the molecules for molecular crystals. The expression “lattice energy” is also in use either as a synonym of cohesive energy or to denote the energy di erence relative to the free ions or molecules [568]. In the calculations of Ecoh we are interested not in the total energy of a system as such, but rather in energy di erences that might be as small as a few kcal/mol. It is with respect to this scale of energy that the overall accuracy of a calculation must be verified. The molecular codes now allow an accuracy of 1 kcal/mol for thermochemical data, which is still far from being attained in solid-state chemistry, although attention to the quantitative aspects of the calculation is increasing rapidly, [568]. The di erent computational aspects of the ground-state total-energy evaluation for crystals are considered in [665]. In particular, some sources of errors in the total-energy calculations are analyzed: 1) the choice of Hamiltonian (HF, DFT, HF+ a-posteriori correlation correction, description of the core electrons by one or other pseudopotential); 2) the use of a finite variational basis; 3)the numerical approximations, introduced in the solution of LCAO equations (truncations of infinite lattice sums, reciprocal-space integration). In the evaluation of the cohesive energy, partially di erent basis sets for the valence or semicore states need to be used for the atoms or ions and the bulk crystal. The bulk basis set can be improved by the addition of polarization functions, which do not contribute to the energy of atoms and ions in the gas phase because of AO orthogonality, but they can be important in the expansion of the bulk wavefunction. The atomic function tails become unnecessary in the description of the crystalline orbitals. Moreover, the use of di use AOs is normally to be avoided (see Chap. 8). Thus, separate optimizations of basis set for the bulk and the isolated atoms are necessary, i.e. not equal, but variationally-equivalent basis sets are to be used. The handling of these problems in the cohesive-energy calculation is illustrated in [568] and the general influence of the Hamiltonian choice on the calculated cohesive energy is discussed. An important aspect of the comparison between the calculated and experimental cohesive energy is related to the di erence between the definition of cohesive energy and the crystal-formation energy that is reported in thermodynamic tables, the main point probably being that quantum-mechanical calculations refer to the static limit (T = 0 K and frozen nuclei), whereas experiments refer to some finite temperature. In fact, the comparison is never straightforward, and the original experimental datum is linked to the calculated cohesive-energy values through a chain of thermodynamic transformations. As an example, the Born–Haber cycle for the formation of a crystal has to be considered (see below).

396 9 LCAO Calculations of Perfect-crystal Properties

A general feature of the HF LCAO method is the underestimation of the cohesive energy (the error varies between –20% and –45% for the series of compounds considered in [568]). DFT calculations of Ecoh allow one to recover part of the contributions that are disregarded with HF, at about the same computational cost. LDA even tends to overestimate the cohesive energy, whereas GGA and B3LYP results are closer to the experimental measurements. At any rate, the performance of none of the Hamiltonians used is fully satisfactory, and the correct answer is always somewhere in between the two extremes represented by HF and LDA, but LDA results are generally improved when gradient corrections are included. Properly correlated wavefunctions perform better than HF and DFT methods: the binding energies for the incremental scheme and local MP2 methods (see Chap. 5) agree better with the experimental data.

Table 9.25. The equilibrium structure and cohesive energy for rutile TiO2 ( [666]). Lattice parameters are given in ˚A, cohesive energy in eV

Observable	HF (AE)*	HF(PP)*	HF+correlation	Experiment
a	4.559	4.555	4.548	4.592
c	3.027	3.024	2.993	2.958
u	0.3048	0.3061	0.3046	0.3048
Ecoh	13.95	15.99	19.35	20.00

*Reference [596]

As an example, we give in Table 9.25 the results obtained for equilibrium structure and cohesive energy in rutile TiO2 by di erent LCAO methods [666]. For the correlation calculations in [666] an incremental scheme was used, one-, twoand threeatom clusters in a large slab of Madelung point charges were included in the CCSD calculations. For the Ti ions nearest to the cluster atoms Ti4+ pseudopotentials were used instead of the bare point charges, in order to simulate the Pauli repulsion on the O2− electrons of the inner cluster. The cohesive energy was calculated including zeropoint vibrations (see Sect. 9.3.3). As is seen from Table 9.25 the electron correlation for TiO2 crystal does not influence the equilibrium structure parameters in a major way but is important for the cohesive energy.

In Table 9.26, a summary of the results from a selection of experimental and recent ab-initio studies of the structural parameters of rutile is presented, see [597]. Early LCAO-HF and PW-LDA studies of rutile yielded lattice parameters to within 2% of experiment. More recent studies that have taken advantage of improvements in the theoretical techniques and available computing power to perform calculations with improved treatments of exchange and correlation (DFT calculations based on the GGA) and higher numerical accuracy improved basis sets have yielded results consistent with those from the earlier work.

Generally speaking, all the ab-initio studies of rutile have yielded structural parameters to within a few per cent of experiment. In particular, LCAO methods give results with the same accuracy as PW methods.