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

9.3 Total Energy and Related Observables in LCAO Methods for Solids

397

Table 9.26. The theoretical and experimental structural parameters for rutile (in ˚A), numbers in parentheses indicate the percent deviation from low-temperature neutron-di raction experiments. The references to the theoretical and experimental data, given in this table, can be found in [597]. SC,LC mean small-core and large-core pseudopotentials, AE –all electron calculations

Method		a		c		u
PW-LDA SC	4.555	(–0.70)	2.922	(–1.08)	0.304	(–0.33)
PW-LDA SC	4.528	(–1.29)	2.918	(–1.21)	0.303	(–0.66)
PW-LDA SC	4.567	(–0.44)	2.932	(–0.74)	0.305	(0.00)
PW-LDA LC	4.536	(–1.11)	2.915	(–1.32)	0.304	(–0.33)
PW-LDA LC	4.653	(1.44)	2.966	(0.41)	0.305	(0.00)
PW-LDA LC	4.603	(0.35)	2.976	(0.74)	0.304	(–0.33)
PW-LDA LC	4.638	(1.11)	2.923	(–1.05)	0.305	(0.00)
PW-GGA LC	4.624	(0.81)	2.992	(1.29 )	0.305	(0.00)
PW-GGA LC	4.690	(2.25)	2.990	(1.22)	0.306	(0.33)
OLCAO-LD AE	4.622	(0.77)	2.983	(0.99)	0.304	(–0.33)
LCAO-HF LC	4.555	(–0.70)	2.998	(1.49)	0.306	(0.33)
LCAO-HF AE	4.560	(–0.59)	3.022	(2.30)	0.305	(0.00)
LCAO-HF AE	4.529	(–1.26)	3.088	(4.54)	0.305	(0.00)
LCAO-HF AE	4.548	(–0.85)	2.993	(1.32)	0.305	(0.00)
LCAO-LDA AE	4.529	(–1.26)	2.942	(–0.41)	0.304	(–0.33)
X ray 298 K	4.594	(0.15)	2.958	(0.14)	0.305	(0.00)
Neutron 295 K	4.593	(0.13)	2.959	(0.17)	0.305	(0.00)
Neutron 15 K	4.587		2.954		0.305

Table 9.27 shows the results of binding-energy HF LCAO calculations [598] in corundum-like Me2O3 oxides (M=Ti,V,Cr,Fe,Co,Ni). The a-posteriori DFT correlation correction of total energy was added in the UHF LCAO AE calculations with numerical optimization of corundum structure parameters.

Table 9.27. Calculated and experimental binding energies (eV) of corundum-like oxides [598]

Oxides	UHF	UHF+corr	Experiment
Al2O3	–23.64	–30.08	–32.14
Ti2O3	–21.41	–28.05	–33.45
V2O3	–19.28	–26.23	–31.21
Cr2O3	–18.10	–25.08	–28.05
Fe2O3	–17.46	–22.00	–25.17
Co2O3	–19.80	–24.90	—
Ni2O3	–17.41	–22.74	—

The atomic energies of Al, Ti, V, Cr, Co, Ni and oxygen were derived by using the basis sets adopted for solid oxides, but supplemented by two more di use

Ecoh(expt) =

398 9 LCAO Calculations of Perfect-crystal Properties

shells, the exponents of which were optimized. The binding energy is computed as the di erence between the total crystal energy per formula unit (for the most stable antiferromagnetic configuration) and the energies of the constituent isolated atoms. For comparison, experimental binding energies are also given. These were obtained by applying a suitable Born–Haber thermochemical cycle of type

Ha,i0 + Hi298 − Hi0 + ∆f H298 − H298 − H0 − Evib0 (9.83)

i

Here, the sum is extended to all chemical elements in the formula unit. The formation enthalpy ∆f H298 of Me2O3 compounds, sublimation (Me) and dissociation (O2) enthalpies Ha,i0 and heating enthalpies Hi298 − Hi0 for all chemical species involved were taken from [667]. The zero-point vibrational energies Evib0 of all oxides were estimated by the Debye model, using an isotropic approximation for the mean acoustic-wave velocity derived from the bulk modulus in order to obtain the Debye temperature.

Table 9.27 shows that the general trend of experimental data (from Al2O3 to Fe2O3) is simulated correctly by theoretical binding energies, with deviations of the order of –30% for UHF energies, which are reduced to –12% by including the DFTbased correction for electron correlation. A secondary feature of the experimental data is not reproduced by the calculated results: the peculiar stabilization energy of Ti2O3 with respect to other oxides. Also, V2O3 seems to be a ected by a larger error than the average. This is consistent with the di culties found in achieving SCF convergence to stable insulating states for these two oxides: oxides of early-row transition metals, with their very di use d orbitals and tendency to metal states, are harder to simulate by HF methods than those of end-row metals. As we have seen in Sect. 9.1 the UHF calculations of Ti2O3 reproduce the electron-correlation e ects on the local properties of chemical bonding. It appears that a-posteriori correction of UHF total energy is not enough to reproduce the correlation e ects on the binding energy, i.e. self-consistent hybrid HF-DFT calculations are required.

In the next section we consider the LCAO calculations of other total-energy-related observables.

9.3.2 Bulk Modulus, Elastic Constants and Phase Stability of Solids: LCAO ab-initio Calculations

We considered above the calculation of observables requiring the first-order derivatives of the total energy. For the equilibrium structure of crystal (equilibrium nuclear coordinates and unit-cell parameters) the total-energy derivatives with respect to nuclear coordinates r and lattice parameters ai equal zero at constant temperature T :

	∂E	T	= 0,	∂E	T	= 0.
	∂ri			∂ai

An important variable is pressure as under pressure interatomic distances in crystals show larger variations than those induced by temperature. At constant temper-

ature T pressure P is related to the rate of energy change with the unit-cell volume

V by relation P = − ∂V∂E T , including the first-order derivative of total energy with respect to the cell volume. Such observables as the bulk modulus, the elastic and force constants depend on the second-order derivatives of the total energy.

The bulk modulus B measures the response of a crystal to isotropic lattice expansion or compression and can be related to the second-order derivative of the total

9.3 Total Energy and Related Observables in LCAO Methods for Solids						399
energy with respect to the volume, V , evaluated at the equilibrium volume V0:
	∂P			∂2E
B = −V		T	= −V		V0	(9.84)
	∂V			∂V 2

The anisotropic response of a crystal to a mechanical force can be described by the elastic constants, Cij , which are defined as the second derivatives of the total energy with respect to the components εi and εj of the strain tensor, ε:

Cij =	∂2E	eq	(9.85)
	∂εi∂εj

	The strain-tensor components εij =	1	∂ui	+ (	∂uj	are first-order derivatives of the
					∂xi
		2	∂xj

displacements u. The elastic constants provide a full description of the mechanical properties of crystalline materials. The bulk modulus is related to the elastic tensor. In the case of a cubic system, where only three independent components of the elastic tensor di er from zero, B can be obtained from C11 and C12 as 13 (C11 + 2C12). Table 9.28 shows the results of the HF LCAO optimization of the Cu2O crystal structure with the all-electron basis set [622]. These results show the feature noted above of the HF method – an overestimation of the lattice size, mainly due to the neglect of electronic correlation. The a-posteriori inclusion of the electronic correlation considerably improves the results of the structure optimization.

			˚			˚
Table 9.28. Calculated lattice constant a (A), interatomic distances (A), bulk modulus
(GPa) and elastic constants (GPa) for Cu2O [622]
		HF	HF+LYP	Expt.
	a	4.435	4.277	4.2696
	Cu–O	1.920	1.852	1.849
	Cu–Cu	3.136	3.025	3.012
	B	100	93	112
	C11	118	103	116	– 126
	C12	90	88	105	– 108
	C44	59	65	11	– 13

The isotropic variation of the volume of the cubic unit cell was used also to evaluate the bulk modulus. The calculated bulk modulus is in good agreement with the available experimental data. The three independent components of the elastic tensor were derived numerically from the changes in energy obtained by applying adequate deformations to the unit cell. Good agreement has been obtained for the C11 and C12 elastic constants, while there is definite disagreement for C44. This disagreement is certainly to be attributed to the lack of precision of these methods to estimate small values of the elastic constants. Good-quality calculations of the elastic constants would require a larger accuracy, in contrast to those properties like the cell parameters.

Various possible structures of a crystal and the associated relative energies can be determined experimentally as a function of pressure. This structure evolution as

400 9 LCAO Calculations of Perfect-crystal Properties

a function of pressure can involve one (polymorphism) or several (solid-state reaction driven by pressure) types of systems. For example, CaO crystal phase transition CaO(B1) CaO(B2) includes two polymorphic phases B1 (fcc lattice, Ca atom is sixfold coordinated) and B2 (sc lattice, Ca atom is eightfold coordinated). The decomposition reaction of MgAl2O4 spinel into its oxide components MgAl2O4 MgO+Al2O3 includes three types of system. The knowledge of pressure enables investigation of phase stability and transitions, [568]. In fact, enthalpy is immediately obtained from the total energy by H = E + P V at T =0 K, where any transformation of a pure substance tends to be isotropic, phase stability can be related to the enthalpy and a phase transition occurs at those points in the phase diagram where two phases have equal enthalpy. From the computational point of view, it is possible to explore a range of crystalline volumes by isometric lattice deformations and obtain the corresponding values of pressure and, consequently, of enthalpy. It is intended that nuclei are allowed to relax to their equilibrium geometry after every lattice deformation.

The theoretical study of phase transitions is made as follows. For each crystal phase, the total energy E is computed at a number of unit-cell volumes V ; at each volume, the structure parameters (nuclear coordinates and lattice parameters) that minimize E are determined. An analytical representation of E vs. V is obtained by using a polynomial expression or the Murnaghan equation of state (or any other fitting function). The Murnaghan function, by far the most universally adopted, is as follows:

"

−

B −1

−

#

1

V

1

E(V ) = E0 + B0V0

B (B 1)

V

+

B V0

−

B

1

(9.86)

The four fitting parameters are V0 (equilibrium volume), B0 (zero-pressure bulk modulus) and B (pressure derivative of the bulk modulus B at P = 0), and E0 (equilibrium energy). From the P (V ) = −∂V∂E relationship, we get

		B	P + 1	−	1
					B
V (P ) = V0				(9.87)
		B0

Inserting (9.87) in (9.86), one obtains the analytic E vs. P dependence; by adding

the P V term, the enthalpy as a function of pressure is obtained:
	−	"		B			1−	1	− 1#
	B							B
H(P ) = E + P V = E0 +	0V0				P + 1					(9.88)
	B 1			B0

At T =0 K, the transition pressure, Pt, corresponds to the point where all the systems have the same enthalpy:

∆H(Pt) = 0

(9.89)

This equation is solved numerically yielding the transition or the decomposition pressure. Knowing this pressure, the equation of state followed by the system during the process can be deduced.

A rough estimate of the transition pressure, Pt , can be obtained just from the knowledge of the equilibrium values (E0 and V0) for each phase. At T =0 K, the enthalpy as a function of pressure can be evaluated as follows:

9.3 Total Energy and Related Observables in LCAO Methods for Solids

401

H(P ) = H0 + P dH = E0 + P V dP		(9.90)
0	0
Using this relation, (9.89) becomes
Pt
∆H(Pt) = ∆E0 + 0	∆V dP = 0	(9.91)

If the pressure dependence of ∆V is negligible, one obtains

Pt Pt = −∆E0/∆V0

(9.92)

For fitting of one or multiphase transitions MULFAS (multiphase transition analysis program) by Llunell can be used. The information about this program can be found at the site www.crystal.unito.it in tutorials (phase transitions) directory. An input file for this program contains the volume and total energy data of one phase obtained by using an external quantum-mechanical or shell-model program.

The first-principles all-electron LCAO calculations of the crystal structures, bulk moduli, and relative stabilities of seven TiO2 polymorphs (anatase, rutile, columbite, baddeleyite, cotunnite, pyrite, and fluorite structures) have been carried out in [597]. From the optimal crystal structures obtained with the Hartree–Fock theory at various pressures, the bulk modulus and phase-transition pressures of various high-pressure polymorphs have been derived at the athermal limit. In most cases, the calculated unit-cell data agree to within 2% of the corresponding experimental determination. In Table 9.29 the calculated bulk moduli of the various phases of TiO2 are given and compared with the existing experimental data. The calculated bulk moduli are within 10% of the most reliable experimental results. It was shown in [597] that the com-

Table 9.29. The bulk moduli (GPa) of various phases of TiO2 compound [597]

	Rutile	Anatase	Columbite	Baddeleyite	Pyrite			Fluorite			Cotunnite
Method	P 42/mnm	I4/amd	P bcn	P 21/c	P a3			F m3m			P nma
	Z=2	Z=2	Z=4	Z=4	Z=4			Z=1			Z=4
LCAO-HF	239±10	202±10	264±10	300±10	318±10			331±10			380±10
LCAO-LDA	241±10	195±10	—	—	273±10			308±10			–
PW-LDA	243	194	247	249	—			282			—
experiment	211±7	178±1	253±4	304±6	—			—			431±10

puted anatase–columbite, rutile–columbite, columbite–baddeleyite, and baddeleyite– cotunnite phase transitions appear in the same order as observed in experiments, and the transition pressures agree semiquantitatively with those measured. The pyrite and fluorite structures are predicted to be less stable than other polymorphs at pressures below 70 GPa, in agreement with experiments.

The detailed LCAO calculations for bulk properties and the electronic structure of the cubic phase of SrTiO3 (STO), BaTiO3 (BTO), and PbTiO3 (PTO) perovskite

402 9 LCAO Calculations of Perfect-crystal Properties

crystals with optimization of basis set (BS) are presented in [606]. The results given in Table 9.30 are obtained using it ab-initio Hartree–Fock (HF) and density-functional theory (DFT) with Hay–Wadt pseudopotentials. A number of di erent exchangecorrelation functionals including hybrid (B3PW and B3LYP) exchange techniques are used. Results, obtained for seven types of Hamiltonians, are compared in Table 9.30 with available experimental data. On average, the disagreement between the lattice constants computed using hybrid HF-DFT functionals and experimental values for all three perovskites is less than 0.5%.

The bulk moduli were calculated in two ways – as the total energy second derivatives (9.84), B1 and using the elastic constants, B2. The results for both ways of bulk-moduli evaluation di er by no more than 10–15%. Especially good agreement with the experimental data has been achieved for hybrid functionals. With the polarization orbitals added to the BS of oxygen atom, the calculated optical bandgaps are 3.57, 3.42 and 2.87 eV for STO, BTO and PTO respectively, in very good agreement with experimental data.

Table 9.30. The optimized lattice constants a (˚A), bulk moduli B (GPa) and elastic constants Cij (in 1011 dyn/cm2) for perovskites SrTiO3 (STO), BaTiO3(BTO), PbTiO3(TO), [606]

Method		LDA	PWGGA	PBE	BLYP	P3PW	B3LYP	HF	Exper.
					STO
a		3.86	3.95	3.94	3.98	3.90	3.94	3.92	3.89
C11		42.10	31.29	31.93	29.07	31.60	32.83	41.7	31.72
C12		12.21	9.80	9.75	9.39	9.27	10.57	7.11	10.25
C44		13.32	11.34	11.30	11.09	12.01	12.46	10.5	12.35
B1		222	170	171	159	167	180	186	174
B2		214	167	169	164	177	177	219	179
					BTO
a		3.96	4.03	4.03	4.08	4.01	4.04	4.01	4.00
C11		35.81	30.11	31.04	28.22	31.12	29.75	30.1	20.60
C12		11.52	10.35	10.72	10.78	11.87	11.57	13.46	14.00
C44		14.98	13.22	13.98	12.24	14.85	14.54	17.34	12.60
B1		196	169	175	166	183	176	190	162
B2		204	175	180	154	188	172	194	195
					PTO
a		3.93	3.96	3.96	4.02	3.93	3.96	3.94	3.97
C11		45.03	32.47	34.25	23.03	43.04	34.42	39.83	22.9
C12		26.14	15.81	15.52	9.93	24.95	18.08	16.90	10.1
C44		11.28	10.69	10.96	8.25	10.93	10.35	17.20	10.0
B1		324	213	217	143	310	235	245	144
B2		321	246	252	140	279	242	299

We note that the results given in Table 9.30, were obtained with the BS optimization (the details of this optimization can be found in [606]). In our opinion, it is this optimization that allowed good agreement of calculated and experimental data to be obtained.

9.3 Total Energy and Related Observables in LCAO Methods for Solids

403

Strontium zirconate SrZrO3 is of interest because of possible applications in fuel cells, hydrogen gas sensors and steam electrolysis. The powder neutron-di raction data suggest the following sequence of phase transitions in SrZrO3: from the orthorhombic P bnm to the orthorhombic Cmcm at 970 K, to the tetragonal I4/mcm at 1100 K and lastly, to the cubic

P m3m

at 1400 K [668]. The first ab-initio DFT PBE calculations for all four phases of SrZrO3 using both PW and LCAO basis sets are presented in [614]. The structural parameter optimization was performed using DFT PBE-PW calculations and good agreement with experimental data was obtained. LCAO PBE calculations were made both for the experimental and optimized structures. Table 9.31 shows the calculated relative energies (per formula unit) of four SrZrO3 polymorphs.

Table 9.31. Calculated relative energies (eV per formula unit) of SrZrO3 polymorphs [614]

			Structure parameters
			Exper.	Optimized
phase			LCAO	LCAO	PW
			0.000	0.000	0.000
P m3m
I4/mcm			–0.040	0.002	0.004
Cmcm			–0.177	–0.252	–0.235
P bnm			–0.241	–0.303	–0.274

It is seen from Table 9.31 that both LCAO and PW results agree with the experimental sequence of the SrZrO3 phases: the most stable turn out to be two orthorhombic phases, whereas the tetragonal and cubic high temperature phases are close in energy.

The sequence of phase transitions in ABO3-type oxides with the temperature increase could be described in terms of the nearest BO6 octahedral tilts/rotations [668]. In the SrZrO3 case, the orthorhombic structure with three rotations is changed by two rotations for Cmcm and finally by one rotaton for the tetragonal structure. The combination of three rotations produces the lowest-energy structure. The electrostatic energy is lowest for the lowest symmetry because of the inherent instability of the ideal corner-shared octahedral network. To compensate the loss in Madelung energy, the repulsive energy should also be lowered, in order to achieve the equilibrium. Thermal expansion or entropy is the driving force for the octahedral tilts.

9.3.3 Lattice Dynamics and LCAO Calculations

of Vibrational Frequencies

The second-order derivatives of the energy with respect to the nuclear coordinates are involved in lattice dynamics, in particular, in the calculation of vibrational (phonon) spectra. We shall begin from the molecular case, [669]. The decoupling of the nuclear from the electronic motion is made in the adiabatic approximation (see Chap. 4).

404 9 LCAO Calculations of Perfect-crystal Properties

Let ui represents a displacement of the ith cartesian coordinate from its equilibrium

value (i = 1, 2, . . . , 3N ), where N is the number of nuclei in a molecule, and q =

√ i

Miui are the generalized coordinate (Mi is the mass of the atom associated with

.

the ith coordinate) and its derivative with respect to time qi = pi. In the harmonic approximation the classical vibrational Hamiltonian of a polyatomic molecule becomes

	1	.	2	+ Hij Uj + V0 =	1
H = T + V =	2	Miui			2	( p|p + q|W |q ) + V0	(9.93)

iij

Here V0 is the electron energy for the equilibrium atomic coordinates, and Hij are the Hessian matrix elements

Hij =	1		∂2V		0		(9.94)
	2		∂ui∂uj
evaluated at equilibrium. The relation Wij =				√		Hij	defines the elements of the
						Mi Mj

weighted Hessian.

The eigenvalues κj of the Hermitian matrix W are the generalized force constants. The Hamiltonian (9.93) then can be factorized into 3N one-dimensional harmonic

	Hamiltonians		1
	H = hν =			Pν2 + ων2Qν2	(9.95)
			2
	ν	ν

Thus each of the 3N −6 vibrational modes can be interpreted as a collective oscillatory movement with frequency ων = √κν /2π and the problem of calculating vibrational spectra reduces to the diagonalization of matrix W to find the set of eigenvalues κj .

For periodic systems, the translation invariance of the potential energy and Hessian matrix should be used. The generalized coordinates obey the Bloch theorem and are

written in the form
qi(k) = N	exp(−ikg)qig	(9.96)

g

The vibrational problem is block-factorized into a set of problems (one for each k point in BZ) of dimension 3N − 6 where N is the number of atoms in the primitive unit cell.

The k-block of the k-factorized W matrix takes the form

	Hij0g
	4
Wij (k) = exp(ikg)		MiMj	(9.97)
g

where Hij0g is the second derivative of potential energy at equilibrium with respect to atom i in the reference cell 0 and atom j in cell g. The number of equations (9.97) to be solved equals the number of k-points in the BZ, i.e. is infinite for the infinite crystal. In practice, the calculations of phonon frequencies are made for a finite number of k-points and the interpolation is used to obtain so-called phonon branches ω1(k), . . . , ωi(k), . . . , ω3N (k) (like one-electron energies are obtained in SCF calculations for a finite set of k-points and then interpolated to form the electronenergy bands). The relationship between phonon frequencies ω and wavevector k determines the phonon dispersion.

9.3 Total Energy and Related Observables in LCAO Methods for Solids

405

Comparing the vibrational branches and electronic bands calculations we note that in the former case the equations for di erent k values are solved independently while in the latter case the self-consistent calculation is necessary due to the BZ summation in the HF or KS Hamiltonian (see Chapters 4 and 7). Once the phonon dispersion in a crystal is known, thermodynamic functions can be calculated on the basis of statistical mechanics equations. As an example, the Helmholtz free energy, F , can be

obtained as:	ik		2 ωik + kB T ln 1 − exp		−kBiT		(9.98)
F =
			1			ω k

where the sum is extended to all lattice vibrations, ωik, and kB is the Boltzmann’s constant. Another way of computing thermodynamic functions is based on the use of the phonon density of states. The evolution of the crystal structure as a function of temperature and pressure can also be simulated by minimizing G = F + pV . The procedure requires a sequence of geometry optimizations, and lattice-vibration calculations [568].

Lattice vibrations can be measured experimentally by means of classical vibration spectroscopic techniques (infrared and Raman) or neutron inelastic scattering. However, only the latter technique allows one to measure the full spectrum in a range of k vectors, whereas with infrared and Raman spectroscopy, only lattice vibrations at Γ (k = 0) are usually detected (the second-order spectra, corresponding to nonzero wavevector k =0 are demanding). The calculations of the vibrational frequencies only at Γ point require the solution of only one equation

	Hij0g		(9.99)
	4
det |W (0)| = 0, Wij (0) =		MiMj
g

To analyze the symmetry of phonon states the method of induced reps of space groups can be used [13]. The procedure of analysis of the phonon symmetry is the following. First, for the space group G of a given crystal the simple induced reps are generated (see Chap. 3). Second, using the simple induced reps together with the compatibility relations and arranging atoms in the primitive cell over the Wycko positions one can determine the symmetry of the phonons. Only those of the induced reps that are induced by the irreps of the site-symmetry groups according to which the components of the vectors of the local atomic displacements transform are used. The total dimension n of the induced rep (called the mechanical representation) equals 3N (N is the number of atoms in the primitive cell).

As an example, Table 9.32 shows the phonon symmetry in rutile TiO2 crystal. The results given are easily obtained from the simple induced reps of the space group D414h (see Table 4.5 in [13]) and the atomic arrangement in rutile TiO2 (see Chap. 2).

There is a one-to-one correspondence between irreps of crystal point group D4h and irreps of the space group at Γ point: A1g,u − 1±, A2g,u − 3±, B1g,u − 2±, B2g,u − 4±, Eg,u − 5±. The atomic displacements of six atoms in the primitive cell generate the 18-dimensional reducible represntation, which contains three acoustic modes and 15 optical modes (A1g + A2g + A2u + B1g + 2B1u + B2g + Eg + 3Eu). Three acoustic modes have zero frequency at the Γ point and are associated with the translation of the entire crystal along any direction in space. These branches are called acoustic modes as the corresponding vibrations behave as acoustic waves. The translation of

406	9 LCAO Calculations of Perfect-crystal Properties
	Table 9.32. Phonon symmetry in rutile TiO2 crystal with space group D414h
		q-basis			k-basis
		Gq	Rep. of Gq	Γ	M	Z	A	X	R
		Ti(2a)	b1u(z)	2−,3−	2−,3−	1	1	1	1−
	(000)		b2u(x − y)	5−	5−	3	3	2	1−
		D2h	b3u(x + y)	5−	5−	4	4	2	1−
		O(4f)	a1(x + y)	1+,4+,5−	1+,4+,5−	1,4	1,4	2,2	1+,1−
		(xx0)	b2(x − y)	2+,3+,5−	2+,3+,5−	2,3	2,3	2,2	1+1−
		C2v	b1(z)	2−,3−,5+	2−,3−,5+	1,3	1,3	1,1	1+1−

the entire crystal along the z-axis corresponds to irrep 3−(A2u) at Γ point, translation in the plane xy – 5−(Eu).

Both acoustic modes are polar and split into transverse A2u (TO) and longitudinal Eu (LO) with di erent frequencies due to macroscopic electric field. All other branches show finite nonzero frequencies at Γ and are known as optical modes, because they correspond to unit-cell dipole moment oscillations that can interact with an electromagnetic radiation.

In the model of a finite crystal the rotation of the entire crystal around the z- and xy-axis corresponds to the irreps 3+ and 5+, respectively, including the displacements of only oxygen atoms, see Table 9.32. As seen from Table 9.32, the phonons of even symmetry (1+, 2+, 3+, 4+, 5+) are connected only with the oxygen-atom displacements. Due to the di erent atomic masses of oxygen and titanium atoms, the corresponding lines can appear in di erent parts of the vibrational spectra, making its interpretation easier. The knowledge of phonon symmetry is useful in the analysis of infrared and Raman spectra of solids as the symmetries of active phonons in these spectra are governed by selection rules following from the symmetry restrictions imposed on the transitions matrix elements [13]. Using the symmetry of the phonons found one can establish which vibrational modes are active in the firstand second-order infrared and Raman spectra.

The calculation of phonon frequencies of the crystalline structure is one of the fundamental subjects when considering the phase stability, phase transitions, and thermodynamics of crystalline materials. The approaches of ab-initio calculations fall into two classes: the linear response method [670] and the direct method, see [671] and references therein.

In the first approach, the dynamical matrix is expressed in terms of the inverse dielectric matrix describing the response of the valence electron-density to a periodic lattice perturbation. For a number of systems the linear-response approach is di cult, since the dielectric matrix must be calculated in terms of the electronic eigenfunctions and eigenvalues of the perfect crystal.

There are two variants of the direct method. In the frozen-phonon approach the phonon energy is calculated as a function of the displacement amplitude in terms of the di erence in the energies of the distorted and ideal lattices. This approach is restricted to phonons whose wavelength is compatible with the periodic boundary conditions

applied to the supercell used in the calculations. Another approach of the direct

method uses the forces ∂E calculated in the total-energy calculations, derives

∂ui eq