Quantum Chemistry of Solids / 14-Hartree-Fock LCAO Method for Periodic Systems

4.1 One-electron Approximation for Crystals

115

greater variational freedom allows the UHF method to produce wavefunctions that are energetically more stable, i.e. give lower electron energy. Another advantage of the UHF method is that it allows solutions with locally nonzero negative or positive spin density (i.e. ferromagnetic or antiferromagnetic systems). However, UHF solutions are formed by a mixture of spin states and are eigenfunctions only of the total spin-projection operator Sz .

In the UHF approach, the single-determinant wavefunction is computed using nα MOs ϕα(r) and nβ MOs ϕβ (r), corresponding to the electrons with α and β spin, respectively. The wavefunction of the system with n electrons (n = nα + nβ ) has the following form:

Ψ = (N !)−1/2 det

ϕ1α(r1)α(1) . . . ϕnαα (rnα )α(nα)

ϕβ

(r (

)β(n

α

+ 1) . . . ϕβ

(r

n

)β(n)

(4.40)

1

nα +1

nβ

The energy expression is given by the following formula:

E = Ψ |Hˆ |Ψ = i

ϕiα|Hˆ + 1/2(Jˆ − Kˆ α)|ϕiα + j

ϕjβ |Hˆ − 1/2(Jˆ − Kˆ β )|ϕiβ

(4.41)

ˆ

ˆ

α

ˆ

β

where the Coulomb J and the exchange K

, K

operators are defined as in the stan-

dard RHF equations, that is:

dr ϕi (r )|r − r |−1ϕi(r ) ϕ(r)

(4.42)

J(ˆ r)ϕ(r) =

i

Kˆ α(r)ϕ(r) = dr (ϕiα) (r )|r − r |−1ϕ(r ) ϕiα(r)

(4.43)

ˆ β (r)ϕ(r) K

=	i		dr (ϕiβ ) (r )|r − r |−1ϕ(r ) ϕiβ (r)	(4.44)
	i

In the LCAO approximation we write now MO for α- and β-spins with di erent

coe cients
ϕα(r) =	Cα		χµA(r); ϕβ (r) =		Cβ χµA(r)	(4.45)
i	iµ		i		iµ
	µA			µA
and introduce the total density Ptot and spin density matrices Pspin
Ptot = Pα + Pβ , Pspin = Pα − Pβ
Pλσα	= i	Ciαλ Ciασ, Pλσβ		= j	Cjβλ Cjβσ	(4.46)

The total density matrix is used to get the expectation value of one-electron operators that are spin independent

ˆ tot |ˆ |

O = Pµν χµ O χν (4.47)

µν

To get the equations for the optimal UHF orbitals we proceed as in the standard RHF case: we compute the first energy variation, then we introduce the orthonormality

116 4 Hartree–Fock LCAO Method for Periodic Systems

constraint of the basis set by means of the Lagrange multipliers method [68]. The Pople–Nesbet equations

FαCα = SCαEα	(4.48)
Fβ Cβ = SCβ Eβ	(4.49)

of the UHF LCAO method can be considered as the extension of the RHF LCAO method equations (4.33) to the case when MO occupied by α and β spins are varied independently. However, (4.48) and (4.49) can not be solved independently as the Fock matrix depends not only on the total density Ptot and has the form:

1					1
Fα = F −		K[Pspin], Fβ = F	+				K[Pspin]	(4.50)
	2					2
where the matrices F and G are defined as
F = H + G; G = J[Ptot] −			1		K[Ptot]			(4.51)
				2

The total electron energy expression now has the form

E = Sp PtotH +	1		1
		Sp PtotG −		Sp PspinK[Pspin]	(4.52)
	2		4

To solve the UHF method equations the iterative procedure is necessary as was already discussed for the RHF method. At first, an initial guess at the density matrices Pα and Pβ is made. Then, the Fock matrices Fα and Fβ are calculated from (4.46). The matrices of coe cients Cα and Cβ are determined by diagonalization of the Fock matrices for α- and β-spins and the density matrix for the next iteration is built. The whole process must be repeated to self-consistency; the convergence criteria consider either the density matrix or the total electron energy. In order to get spin-polarized solutions from an UHF calculation, it is necessary that the initial density matrices Pα and Pβ are chosen to be di erent and this di erence maintained during the whole SCF process [69]. The RHF closed-shell solution can be considered as a special case of the UHF method, where the number of α electrons coincide with the number of β electrons (nα = nβ ). In the RHF method Pα = Pβ for the whole SCF process. In the UHF method (when initial matrices Pα =Pβ ) we could have either spin-polarized solutions at the end of the SCF cycle or the spin polarization will vanish in the course of the SCF process. Reaching a closed-shell solution from a spin-polarized starting point is favored when the latter solution has lower energy.

The extension of the RHF, ROHF and UHF LCAO methods to periodic systems is considered in the next sections. We begin this consideration from the discussion of specific features of these methods when instead of a molecular system the cyclic model of a crystal is introduced.

4.1.4 Specific Features of the Hartree–Fock Method for a Cyclic Model of a Crystal

The extension of the HF method to a model of an infinite crystal leads to di culties connected with the one-electron Bloch functions (BF) behavior since they do not fall

4.1 One-electron Approximation for Crystals

117

o to zero at infinity and therefore can not be normalized to unity in a space [70]. To overcome this di culty, periodic boundary conditions (PBC) are introduced for the BF normalized in the main region of volume of the crystal containing a large number N of primitive unit cells of volume.

In the HF method for molecular systems a many-electron wavefunction in the form of a Slater determinant is written for the system containing a fixed number of electrons. PBC introduction for the infinite crystal allows the system of infinite number of electrons to be replaced by the system of a finite number of electrons (a cyclic cluster). But, in a cyclic cluster containing a large number N of primitive unit cells the number of electrons remains huge and from the very beginning it is necessary to take into account the translation symmetry that is an opportunity to generate a whole cyclic cluster by translations of a primitive unit cell. The concept of the cyclic cluster helps to understand why, on calculation of the Fock matrix for a crystal, it is necessary to make summations on a direct lattice and on the Brillouin zone. Thus, the dimension of the Fock matrix is determined by the number of the basis atomic functions referring not to an infinite crystal or its main region but to its small part – the primitive unit cell.

Let the infinite crystal be replaced with a cyclic system (cyclic cluster) from N = N1 × N2 ×N3 primitive cells, i.e. replaced with a parallelepiped with edges Niai(i = 1, 2, 3). Periodic boundary conditions mean that opposite sides of a parallelepiped are identified, i.e. translations on vectors Niai(i = 1, 2, 3) are equivalent to zero translation. The number of electrons in a cyclic system is equal N × n, where n is the number of electrons on a primitive unit cell. In the HF method the cyclic system under consideration can be considered as a huge molecular system with periodical boundary conditions imposed on molecular orbitals. The matrix elements of the Fock operator for the cyclic cluster in the LCAO basis have to include the summation on all atoms in the cyclic cluster. To calculate the one-electron density matrix elements the summation has to be performed on all the molecular orbitals occupied by N × n electrons. We shall notice that for essentially ionic systems the cyclic cluster under consideration should be placed in a Madelung field of the rest of an infinite crystal. This means that in the diagonal elements of the Fock matrix contributions should be added from the field of point charges of those atoms that are not included in the cyclic cluster. With such an approach interatomic interactions in a crystal are actually separated into short-range and long-range ones. The former are taken into account directly in the cyclic-cluster calculation, the latter are approximated by inclusion of the Madelung field acting on every atom of the cyclic cluster under consideration. The e cient techniques of the Madelung field calculation are developed in the solid-state theory [71] and we do not discuss them here. In [55] the symmetry group of the cyclic model of a crystal was examined and the connection of this group with the symmetry group of an infinite crystal was discussed. In particular, it was shown that the inner translations of a cyclic system form a subgroup of the translation group of an infinite crystal. This allows calculation of the electronic structure of the cyclic system using the atomic basis symmetrized over the operations of this translation subgroup. This means use of LCAO Bloch functions for those wavevectors that are connected with the group T(N ) of inner translations of cyclic system. By this the Fock matrix is split up into blocks whose dimension is determined by number of di erent one-dimensional representations of the translation group T(N ), i.e. those nonequivalent wavevectors

118 4 Hartree–Fock LCAO Method for Periodic Systems

in the Brillouin zone that correspond to these representations. These wavevectors are defined by the relation exp(−ikA) = 1, where A is the translation vector of the cyclic system as a whole. Comparing the HF method for crystals with that for molecules we conclude that instead of one the system of HF equations the separate equations for any of nonequivalent wavevectors is necessary to solve at each iteration of the self-consistent procedure. The preliminary atomic basis symmetrization in the Hartree–Fock calculations of molecular systems can also be carried out to split up the Fock matrix into blocks of smaller dimensions. However, the gain of such a splitting, as a rule, is insignificant especially for those point groups that have many dimensional irreducible representations. For a crystal, the direct account of translation symmetry is a unique way to calculate the electronic structure. The summation over a direct lattice in the Fock matrix appears as the sum over atoms included in the cyclic cluster. The use of the translation symmetry allows one to calculate only those matrix elements of the Fock matrix that are connected with the AOs of the reference primitive unit cell (the zero translation vector corresponds to the reference primitive cell). In the density-matrix calculation for a cyclic system the summation over wavevector appears. This summation means the use of an atomic basis symmetrized over the translation subgroup of a cyclic system symmetry group and replaces the summation over the occupied MO of a cyclic system by summation over occupied crystalline orbitals for di erent wavevectors. Their number is determined by the number of electrons in a primitive unit cell. The further splitting of the Fock matrix can be made using the connection between the irreducible representations of the symmetry group of a cyclic system and those of an infinite-crystal symmetry group. For each wavevector only the one ray of its star can be considered, and summation on a star is carried out by application of operations of the point group of a crystal to that crystalline orbital calculated for this ray. The full space group symmetry inclusion in the HF LCAO method is considered in [28, 29] and realized in the CRYSTAL code [23].

Naturally, there is a question: what size is it possible to be limited to in cyclic system modeling of each an actual crystal? This depends on the character of localization of electronic density in a crystal: in ionic crystals and semiconductors with a large energy gap (5–10 eV) the electronic density is well-enough localized on ions or interatomic bonds so the cyclic system of rather small sizes describes a short-range part of the interatomic interactions. This means that on calculation of the one-electron density for such crystals it is possible to be limited to summation on a rather small number of points of the Brillouin zone. In crystals with a small energy gap, or metals, the size of the needed cyclic system becomes essentially larger as the one-electron states are more delocalized (in metals these states are practically delocalized over the whole infinite crystal). Therefore, the summation over the Brillouin zone in calculations of crystals with a small energy gap requires many k-points to be included. We see that one-to-one correspondence between the cyclic system size in direct lattice and the finite set of wavevectors in the Brillouin zone allows an investigation of the convergence of the results of electronic-structure calculations on the size of the cyclic system su cient to describe short-range interatomic interactions in a crystal. Nevertheless, such an investigation is made by the primitive unit cell calculations.

In the next sections we consider in more detail the HF LCAO method for periodic systems.

4.1 One-electron Approximation for Crystals

119

4.1.5 Restricted Hartree–Fock LCAO Method for Crystals

The reformulation of Hartree–Fock LCAO molecular equations (4.33) to make them suitable for periodic systems was proposed in [74] about forty years ago. However, the first HF LCAO calculations were made mainly for the simple solids [72] or for one-periodic systems (polymers) [73, 74] in which the direct lattice summations are essentially simpler than in the twoor three-periodic systems. The periodic HF LCAO method for the systems with any periodicity (polymers, surfaces, bulk crystals) was implemented in the CRYSTAL code [23] after the detailed theoretical study of the exact-(nonlocal) exchange Hamiltonian for the periodic systems. The following main aspects were examined [75]: (i) the point and translational symmetry properties of the HF LCAO Hamiltonian and the transformation of all relevant monoelectronic and bielectronic terms under the symmetry operators, (ii) the use of general yet powerful criteria for the truncation of infinite sums of Coulomb and exchange integrals, (iii) the use of a procedure for performing integrals over wavevector k as needed in the reconstruction of the Fock matrix in the SCF calculation. The first successful applications of the CRYSTAL code to the calculations of the ground-state properties of the two-periodic [76] (hexagonal boron nitride) and three-periodic [77](cubic boron nitride and silicon) have opened a more than 25-year ab-initio HF LCAO study of the bulk electronic properties of crystals, both perfect and defective, the extension of the method to the surface modeling and magnetic systems. The density functional theory (DFT) LCAO and hybrid (HF-DFT) calculation schemes were implemented in the later versions of the CRYSTAL code [78]. The next step was made when the EMBED code [79] appeared with the implementation of the perturbed cluster model of point defects in crystals. The detailed discussion of many results obtained by HF and DFT LCAO methods we postpone until Chapters 9–11. The reader interested in the list of the actual systems calculations is referred to the Internet sites [80] and [81]. Recently, periodic boundary conditions were included in the molecular code Gaussian [82] and the first-principles HF and DFT calculations were made with the use of this code [83]. In Appendix C we give brief information about these two and some other HF and DFT LCAO computer codes for crystals.

In this section, we examine the main modifications in the Hartree–Fock–Roothaan equations (4.33), it being necessary to take into account the translation symmetry of periodic systems. The first most important di erence appears in the LCAO representation of the crystalline crbitals (CO) compared to the molecular orbitals (MO).

For the periodic system the AO symmetrization over the translation subgroup of

the space group is made giving the Bloch sums of AOs:
1		exp(ikRn)χµ(r − Rn)	(4.53)
χµk(r) = √N Rn
In (4.53), the index µ labels all AOs		in the reference primitive unit	cell (µ =

1, 2, . . . , M ) and Rn is the translation vector of the direct lattice (for the reference primitive cell Rn = 0). The summation in (4.53) is supposed to be made over the infinite direct lattice (in the model of the infinite crystal) or over the inner primitive translations R0n of the cyclic cluster (in the cyclic model of a crystal). In the latter case, the sum of the two inner translation vectors R0n + R0m = Rl may appear not to be the inner translation of the cyclic cluster. However, the subtraction of the translation vector A of the cyclic cluster as a whole (in the cyclic model the vector A is

120 4 Hartree–Fock LCAO Method for Periodic Systems

supposed to be equivalent to the zero translation) gives the inner translation again:

R0l = Rl − A.

The Bloch sum of AO (4.53) satisfies the Bloch theorem. To check this let us apply

						ˆ	(for the fixed translation Rm) to the sum (4.53)
the translation operator tRm
				ˆ
				tRm χµk(r) = χµk(r − Rm)
		1
= exp(−ikRm)			√			Rn exp(ik(Rn + Rm))χµ(r − Rm − Rn)
				N
= exp(	−	ikRm)		n			exp(kRn )χµ(r	−	−
									Rn ) = exp( ikRm)χµ(r) (4.54)
			R		=Rn +Rm

The replacing of the sum over translation vectors (Rn + Rm) by the sum over translation vectors Rn is evidently possible for the infinite crystal. In the cyclic cluster this replacement means the subtraction of the translation vector A.

In the LCAO approximation CO is the Bloch function as it is expanded over the

Bloch sums of AOs:
ϕik(r) =	Ciµ(k)χµk(r)	(4.55)

µ

In the MO LCAO approximation the index i in the expansion (4.32) numbers the MOs (their total number M is equal to the number of atomic orbitals used in the expansion (4.32)). In the case of closed shells the Ne electrons of the molecule occupy Ne/2 MOs and (M − Ne/2) MOs are empty. The total number of COs of the cyclic cluster equals M ×N (N is the number of the primitive unit cells in the cyclic cluster; M is the number of AO basis functions per primitive unit cell). For the cyclic cluster containing Ne = N ×n electrons (n is the number of electrons per primitive unit cell) Ne/2 crystalline orbitals are occupied by electrons and (M × N − Ne/2) orbitals are empty. The numbering of crystalline orbitals is made by two indices i and k: the oneelectron states of a crystal form the energy bands numbered by index i and joining the N states with the same i. The closed-shell case (the nonconducting crystals) means that all the energy bands are filled or empty (there are no partly filled energy bands). For N → ∞ the total number of COs also becomes infinite, the one-electron energy levels form the continious energy bands, but the total number of energy bands remains finite and is equal to M . For the closed-shell case n/2 electrons occupy the lowestenergy bands for each k-vector. As discussed in Sect. 3.2.5 MgO and Si crystals are examples of nonconducting systems with di erent values of the forbidden energy gap (Figs. 3.4–3.5). The latter is the crystalline analog of the HOMO (highest-occupied MO)–LUMO(lowest-unoccupied MO) one-electron energies di erence. In both cases, n = 8 (there are 8 valence electrons per primitive cell), 8 electrons occupy 4 valence energy bands (for MgO crystal the whole valence band is split into two subbands).

On the basis of Bloch functions, the Fock and overlap matrices F and S become

Fµν (k) =	exp(ikRn)Fµν (Rn); Sµν (k) = exp(ikRn)Sµν (Rn) (4.56)
Rn	Rn

where Fµν (Rn) is the matrix element of the Fock operator between the µth AO located in the reference (zero) cell and the ν-th AO located in the Rn cell. The matrix element Sµν (Rn) is the overlap integral of the same AOs. The row index can be limited to the

4.1 One-electron Approximation for Crystals

121

reference cell for translational symmetry. Matrices represented on the Bloch basis (or in k-space) take a block-diagonal form, as Bloch functions are bases for irreducible representations of the translation group T (T(N ) for the cyclic cluster of N primitive cells); each block has the dimensions of the AO basis in the primitive cell, M ,

F(k)C(k) = S(k)C(k)E(k)

(4.57)

In the HF LCAO method, (4.57) for the periodic systems replaces (4.33) written for the molecular systems. In principle, the above equation should be solved at each SCF procedure step for all the (infinite) k-points of the Brillouin zone. Usually, a finite set {kj }(j = 1, 2, . . . , L) of k-points is taken (this means the replacing the infinite crystal by the cyclic cluster of L primitive cells). The convergence of the results relative to the increase of the k-points set is examined in real calculations, for the convergent results the interpolation techniques are used for eigenvalues and eigenvectors as these are both continuous functions of k [84]. The convergence of the SCF calculation results is connected with the density matrix properties considered in Sect. 4.3

The overlap matrix elements in the AO basis are the lattice sums of the overlap integrals between AOs, numbered now by indices of AO in the zero cell and of the cell defined by the translation vector Rn:

Sµν (Rn) = χµ(r)χν (r − Rn)d3r

(4.58)

The Fock matrix elements (as in the case of molecules) are the sums of one-electron (the kinetic energy T and the nuclear attraction energy Z) and two-electron (Coulomb J and exchange energies) parts. The di erence of these matrices from the molecular analogs is the appearance of the sums over a direct lattice, containing one-electron (kinetic energy and nuclear attraction) and two-electron integrals:

Fµν (Rn) = Tµν (Rn) + Zµν (Rn) + Jµν (Rn) + Xµν (Rn)

			1
	Tµν (Rn) = −			χµ(r) r2 χν (r − Rn)d3r
			2
		NA			| − −		|
Zµν (Rn) = A=1 Rm						Z
			χµ(r)		A			χν (r − Rn)d3r
					r A Rm
		(A numbers atoms in the primitive cell)
	M
					(µ0νRn|λRmσ(Rm + Rm ))
Jµν (Rn) =		Pλσ(Rm )
	λσ Rm			Rm
1		M
Xµν (Rn) = −	2	Pλσ(Rm )				(µ0λRm|νRnσ(Rm + Rm ))
		λσ Rm		Rm

(4.59)

(4.60)

(4.61)

(4.62)

(4.63)

The two-electron integrals in (4.62) and (4.63) are defined as

(µ0νRn|λRmσ(Rm + Rm )) = χµ(r)χν (r − Rn)|r − r |−1

χλ(r − Rm)χσ(r − Rm − Rm )d3rd3r (4.64)

122 4 Hartree–Fock LCAO Method for Periodic Systems

The comparison of (4.36) and (4.64) shows the di erence in the two-electron integrals definitions for crystals and molecules.

The spinless one-electron density matrix (DM) elements are defined in the LCAO approximation as

Pλσ(Rm ) = 2		dkeikRm		Ciλ(k)Ciσ(k)θ( F − i(k))	(4.65)
			i

where F is the Fermi energy, the integration in (4.65) extends to the first Brillouin zone and corresponds to the model of the infinite crystal. In the cyclic-cluster model the integration is replaced by the summation over k-points numbering the irreps of the group of the inner cyclic-cluster translations. The DM elements depend on the eigenvectors Ci(k) of the F(k) matrix. The integration or summation over k in DM elements makes it impossible to solve (4.57) for di erent k independently.

The electron energy of the crystal (per primitive unit cell) as calculated within the HF LCAO approximation can be expressed in terms of the one-electron density matrix (DM) and includes the lattice sums

1		M
Ee =	2	Pµν (Rn) (Fµν (Rn) + Tµν (Rn) + Zµν (Rn))	(4.66)
		µν Rn

The direct lattice summations in the Fock matrix elements and the k dependence of the one-electron DM, energy levels and COs are the main di culties of the HF LCAO method for periodic systems, compared with molecules. A special strategy must be specified for the treatment of the infinite Coloumb and exchange series as well as for the substitution of the integral that appears in DM with a weighted sum extended to a finite set of k-points. The e cient solution of these problems has been implemented in the CRYSTAL code [23]. These problems are also valid for UHF and ROHF LCAO methods for periodic systems considered in the next subsection.

4.1.6 Unrestricted and Restricted Open-shell Hartree–Fock Methods for Crystals

The unrestricted and restricted open-shell Hartree–Fock Methods (UHF and ROHF) for crystals use a single-determinant wavefunction of type (4.40) introduced for molecules. The di erences appearing are common with those examined for the RHF LCAO method: use of Bloch functions for crystalline orbitals, the dependence of the Fock matrix elements on the lattice sums over the direct lattice and the Brillouinzone summation in the density matrix calculation. The use of one-determinant approaches is the only possibility of the first-principles wavefunction-based calculations for crystals as the many-determinant wavefunction approach (used for molecules) is practically unrealizable for the periodic systems. The UHF LCAO method allowed calculation of the bulk properties of di erent transition-metal compounds (oxides, perovskites) – the systems with open shells due to the transition-metal atom. We discuss the results of these calculations in Chap. 9. The point defects in crystals in many cases form the open-shell systems and also are interesting objects for UHF LCAO calculations (see Chap. 10).

4.1 One-electron Approximation for Crystals

123

In the UHF or ROHF cases two sets of matrix equations are solved self-consistently for α and β spins [69]:

Fα(k)Cα(k) = S(k)Cα(k)Eα(k)	(4.67)
Fβ (k)Cβ (k) = S(k)Cβ (k)Eβ (k)	(4.68)

where the Fα(k), Fβ (k) and S(k) matrices are obtained by Fourier transform from the corresponding direct space equivalent quantities as in (4.56).

The Fourier components of the density matrix for α and β spins are di erent as they are obtained by using the Cα(k) and Cβ (k) eigenvectors obtained from (4.67) and (4.68), respectively:

	Pλσα (Rm) =	dkeikRm		niαCiλα(k)Ciσα (k)	(4.69)
			i
	Pλσβ (Rm) =	dkeikRm		niβ Ciλβ (k)Ciσβ (k)	(4.70)
			i
where nα,β	are the occupations of CO by electrons with α and β spins.
i

The total density P(Rm) and spin density Pspin(Rm) Fourier components are

defined as
	P(Rm) = Pα(Rm) + Pβ (Rm)									(4.71)
Pspin(Rm) = Pα(Rm) − Pβ (Rm)										(4.72)
The Fα(Rm) and Fβ (Rm) matrices are defined as follows:
F α	(R	m	) = F	µν	(R	m	) + Xspin(R	m	)	(4.73)
µν							µν
Fµνβ (Rm) = Fµν (Rm) − Xµνspin(Rm)										(4.74)

The F(Rm) matrix is defined as in (4.59), where the total density matrix P(Rm) defined in (4.71) is used in the Coulomb and exchange terms; Xµνspin(Rm) is defined as Xµν (Rm) in (4.63), where, however, the spin density matrix Pspin(Rm) is used instead of the total density matrix P(Rm).

The ROHF LCAO method for crystals di ers from the UHF method in the equations defining the CO. Let n = nα + nβ electrons per primitive unit cell be considered (nα > nβ ). For the nβ CO the closed-shell RHF LCAO equations are solved, the CO of nα + nβ electrons with α spin are found from the equation analogous to (4.57). In the ROHF method the total density P(Rm) and spin density Pspin(Rm) matrices are defined as

P(Rm) = Pα,β (Rm) + Pα(Rm)	(4.75)
Pspin(Rm) = P(α−β)(Rm)	(4.76)

In (4.75) Pα,β is the DM calculated with the doubly occupied COs and in (4.76) P(α−β) is DM calculated for COs with α spin.

The UHF method many-electron wavefunction is the eigenfunction of the spin pro-

ˆ

jection operator Sz with zero eigenvalue, the ROHF method many-electron function

124 4 Hartree–Fock LCAO Method for Periodic Systems

describes the spin state with the maximally possible spin projection Sz = 12 (nα −nβ ). If we compare the UHF scheme with the RHF scheme for the closed-shell case, we see that the most relevant quantities (Fock and density matrices, eigenvalues and eigenvectors) are duplicated, so that it is possible to di erentiate α orbitals from β orbitals. As in the case of molecules it is important that at the starting guess Pα and Pβ density matrices are di erent. In the CRYSTAL code there exists an option to construct a DM formed by the superposition of atomic solutions, where it is possible to assign a specific polarization (or no polarization) to each atom of the system [69]. Once the SCF process is started, it is often useful to introduce some constraints on the populations of the α and β orbitals. This is performed in the CRYSTAL code by ad-hoc translations of the eigenvalues spectrum in such a way that the di erence between the number of α and β electrons is equal to the value introduced in the input.

It is well known from Hartree–Fock studies of molecular systems, that it is very common to have problems of SCF convergence when studying open-shell systems; similarly, convergence problems are not rare in the Hartree–Fock treatment of spinpolarized crystals. A well-known technique for the solution of convergence problems, in the case of open-shell molecules, is the so-called level-shifting method [85]; this approach has shown its e ectiveness in the periodic HF context also, especially in the case of crystals containing transition elements.

4.2 Special Points of Brillouin Zone

4.2.1 Supercells of Three-dimensional Bravais Lattices

Let ai(Γ1) (i=1,2,3) be the basic translation vectors of the initial direct lattice of type Γ1 and aj (Γ2) (j = 1, 2, 3) be the basic translation vectors of a new lattice of type Γ2 with the same point symmetry (symmetrical transformation) but composed of supercells. Then

aj (Γ2) = lji(Γ2Γ1) ai(Γ1) |det l| = L, j = 1, 2, 3 (4.77)

i

where lji(Γ2Γ1) are integer elements of the matrix l(Γ2Γ1) defining the transition from the lattice of type Γ1 to the lattice of type Γ2.

The vectors aj (Γ2) have well-defined orientation with respect to point-symmetry elements of the lattices that are the same for both lattices because of the symmetrical character of the transformation (4.77). Let us define the components of the vectors aj (Γ2) by the parameters sk assuring their correct orientation relative to the lattice symmetry elements and the correct relations between their lengths (if there are any). Then three vector relations (4.77) give nine linear nonhomogeneous equations to determine nine matrix elements lij (Γ2Γ1) as functions of the parameters sk. The requirements that these matrix elements must be integers define the possible values of the parameters sk giving the solution of the problem.

Let us demonstrate the procedure of finding the matrix of a symmetrical transformation (4.77) by the example of the rhombohedral crystal system where there is only one lattice type (R). The basic translation vectors of the initial lattice are the following: