Quantum Chemistry of Solids / 15-Electron Correlations in Molecules and Crystals

5.1 Electron Correlations in Molecules: Post-Hartree–Fock Methods

157

matrix elements between the ground state and any single-excited configuration vanish [8]. First, the second-order term is nonvanishing. This term is

		ˆ α,β α,β			ˆ
	Φ0	|∆H|Φi,j	−		|∆H|Φ0
							(5.30)
i,j α,β		E0		Eα,β
				i,j

where we have used the fact that only those configurations Φα,βi,j where exactly two electrons have been excited from the ground state Φ0 have nonvanishing matrix el-

ˆ

ements in the denominator. To prove this, one uses the fact that ∆H consists of two-electron operators and that there will be nonvanishing elements only when the two determinants di er by at most two orbitals. It can also be proved that there is no contribution from single-excited configurations. Furthermore,

N
i	(5.31)
E0 = i − E = EHF
=1
and
Ei,jα,β = EHF + α + β − i − j	(5.32)

are the eigenvalues of G . The most important aspect is now that the denominator in (5.30) equals

− α − β + i + j

(5.33)

due to (5.32). The two first energies are energies of orbitals that for the ground state are vacant, whereas the last two are those of occupied orbitals. This means that the denominator is small (and, hence, the correlation e ects are large) when there is a small energy di erence between occupied and unoccupied orbitals. This partly explains the results for the H + H system where the correlation is important for large interatomic distances, when the energies of the occupied and unoccupied orbitals approach each other. The result tells us also that, e.g., for compounds containing transition-metal atoms, where there are many (empty and occupied) d orbitals close to each other, we will expect that correlation e ects are important, or, alternatively, that the Hartree–Fock approximation (i.e. the single-Slater-determinant approximation) is not a very good one. Turning to the cyclic model of a crystalline solid we conclude that the smaller the bandgap the larger the correlation e ects on the ground-state energy.

The higher-order terms to the perturbation series can also be calculated. These will become increasingly complex, but already second-order MP2 perturbation theory provides a very important improvement over the pure Hartree–Fock approximation, and in very many cases results from this method are highly accurate. The higher-order calculations are given the labels MP3, MP4. We do not discuss here many numerical results for the ground state of molecular systems obtained in CI, MCSCF, CC or perturbation theory calculations, rather we refer the reader to monographs [5, 8, 102]. Meanwhile it is clear that direct extension of CI or CC methods to the molecular cluster models of periodic systems is a di cult task. The molecular cluster reasonably models the periodic system when it is chosen su ciently large to suppress the dangling bonds influence. Even for ionic crystals such a cluster can consist of up to 50–100 atoms. The one possibility of the calculation e orts decreasing is connected with the use of localized MOs in the Slater determinants. The local electron correlation methods are considered in the next subsection.

158 5 Electron Correlations in Molecules and Crystals

5.1.5 Local Electron-correlation Methods

In the wavefunction-based improvements upon mean-field HF method to account for electron correlation (CI, CC, MP) the demands on the computational resources (CPU time, memory or disk space) exhibit a strong dependence on the system size. Let the system consist of N atoms. If we assume a given basis set, the cost of MP2 calculations scales between N 3 and N 5 (asymptotic limit) [127]. Higher correlation methods CCSD and CCSD(T) scale as N 6 and N 7, respectively. The latter means that doubling the size of the molecular system increases the computational costs by a factor of 64 and 128! This high scaling behavior restricts the application range for these accurate and reliable methods to only small molecules, a traditional dilemma in quantum chemistry for a long time (especially, since density-functional theory became available).

Local electron-correlation methods are ab-initio wavefunction-based electronicstructure methods that exploit the short-range nature of dynamic correlation e ects and in this way allow linear scaling O(N ) in the electron-correlation calculations [128, 129, 131–135] to be attained. O(N ) methods are applied to the treatment of extended molecular systems at a very high level of accuracy and reliability as CPU time, memory and disk requirements scale linearly with increasing molecular size N .

Linear-scaling electron-correlation methods were developed by: 1) combining the Pulay–Saebo local correlation variant [136] with integral-direct techniques [131] and consequently exploiting the spatial locality of the electron-correlation e ect; 2) Laplace-transform techniques suggested in [137] and applied in [129].

In a first step, a linear-scaling local MP2 (LMP2) algorithm was implemented [129, 131, 132]. With this algorithm, MP2 calculations on spatially extended molecular systems including more than 200 atoms and 2000 basis functions can easily be performed on a current low-cost Pentium PC. The linear-scaling approach was also developed for high-level correlation methods, i.e. linear-scaling, local CCSD with (T) correction ( [133, 134, 142, 143]). In particular, for the CCSD(T) method, reductions in required CPU time were obtained by more than a factor of 1000 compared to the conventional (T) method.

The chemical systems well beyond 1000 basis functions, 100 atoms and some hundred correlated electrons are now within reach of the local electron-correlation methods. It is clear that these methods still are considerably more expensive than modern, linear-scaling DFT algorithms (see Chap. 7), but this now is just due to a higher prefactor, and no longer due to less favorable scaling behavior. Thus, the longstanding scaling problem of ab-initio electronic-structure theory clearly is defeated. On the other hand, local correlation methods combine all the well-known advantages of wavefunction-based methods over DFT, like reliability, possibilities for systematic improvements, etc.

At first, all these methods were developed for closed-shell systems only. Later research in this area was directed towards local methods for open-shell systems and excited states, local triples corrections beyond (T) (triples included in coupled cluster iterations), [138], local energy gradients for geometry optimizations of large molecules [139], combination of the local correlation method with explicitly correlated wavefunctions. It is evident from the discussion that these local O(N ) methods open the applications of coupled-cluster theory to entirely new classes of molecules, which were far out-of-scope for such an accurate treatment before. Possible applications lie, for example, in the determination of the thermochemistry of reactions involving

5.1 Electron Correlations in Molecules: Post-Hartree–Fock Methods

159

large closed-shell (or singlet biradiaclic) organic molecules including the computation of accurate barrier heights (which is an area, where DFT often fails badly). Furthermore, once local MP2 NMR shifts become available, NMR spectra of extended organic molecules could reliably be predicted, which might be helpful in the determination of the structure of such molecules.

Another potential application field are intermolecular complexes and clusters [140]. Here, the local correlation method has some conceptual advantages, insofar that it eliminates the notorious basis-set superposition error (BSSE) [141]. This allows for the BSSE-free calculation of interaction energies structures of large molecular clusters (i.e. involving many monomers) without the expensive counterpoise correction that is required in conjuction with canonical methods. There are many more application areas one could think of.

One of the most impressive applications of the local correlation methods is their extension to the periodic systems. Having in mind to discuss this extension on the MP2 level in the next section we restrict the discussion here to the consideration of local MP2 methods for molecular systems with closed shells, referring the reader to the papers discussing local CC methods [128, 133, 134, 142, 143].

In the local MP2 method [131] the occupied canonical MOs of the preceding HF calculation are localized using the Pipek–Mezey [37] localization procedure and keeping the localized orbitals orthonormal. Even though the localization is formally an O(N 3) step, the localization time is negligible for all the local MP2 calculations.

The use of local orbital bases opens the way for two distinct approximations. First, for the correlation of each electron pair ij an individual subset (domain) of the virtual orbitals that is independent of the molecular size can be selected. This reduces the scaling of the number of configuration-state functions (CSFs) and corresponding amplitudes from O(N 4) to O(N 2). Secondly, a hierarchical treatment of di erent electron pairs depending on the minimum distance R of the two correlated localized occupied MOs (LMOs) i and j can be devised. In implementation [131], a subset of atoms is assigned to each LMO according to the procedure of [144]. The pairs (ij) are then classified according to the minimum distance between any atoms in the two di erent subsets. Strong pairs (R ≤ 1 Bohr) are treated at the highest level, e.g., local CCSD. These involve pairs of local orbitals that share at least one atom and typically account for more than 90% of the correlation energy. Weak pairs (1 < R ≤ 8 Bohr) and distant pairs (8 < R ≤ 15 Bohr) are (optionally) treated by local MP2, while very distant pairs (R > 15 Bohr) can be entirely neglected (in [131], strong and weak pairs are both treated by LMP2). For the distant pairs, the required two-electron integrals can be obtained using multipole approximations [145] leading to substantial savings. The numbers of strong, weak, and distant pairs scale only linearly with molecular size N whereas the number of very distant pairs scales quadratically with N . The neglect of very distant pairs leads to an overall linear dependence of the number of CSFs and of the corresponding transformed integrals. The multipole treatment of distant pairs further reduces the number of integrals that must be constructed via the four-index integral transformation; this fact can be exploited to devise an integral transformation algorithm for which all computational resources (CPU, memory, and disk) scale only linearly with N [131].

In the basis {χµ} with overlap matrix Sµν = χµ|χν the localized occupied orbitals |φi are represented by a coe cient matrix L

160	5 Electron Correlations in Molecules and Crystals
			|χµ Lµi
		|φi =				(5.34)
		µ
					˜
The virtual space is spanned by a set of nonorthogonal functions {φr} that are ob-
tained directly from the AOs {χµ} by projecting out the occupied space, i.e.
	|φ˜r = 1 −	m
		i
		=1 |φi φi| |χµ =		µ	|χµ Pµν	(5.35)
where the projection matrix is defined as
		P = I − LL†S				(5.36)

Here indices i, j, . . . , are used for the occupied LMOs, r, s, . . . , for projected AOs, and µ, ν, ρ, σ for AOs. Quantities in the projected basis are capped by a tilde.

| |˜

To each LMO φi , a subset [i] (orbital domain) of the projected AOs φr is assigned. The orbital domains are selected in [131] as described by Boughton–Pulay [144] using a threshold of 0.02, but in addition discarding centers with a Mulliken gross charge below 0.01 (0.03 for H atoms). The electron pairs (ij) were classified as strong, weak or distant and very distant pairs (R > 15 Bohr) were neglected. The correlation space for the strong and weak electron pairs (ij) was spanned by pair domains [ij] = [i] [j], which are the union of the two related orbital domains. Distant pairs were treated by a multipole approximation, as described in [145]. This requires the use of asymmetric pair domains r [i], s [j] for this class of configurations, i.e. long-range ionic excitations are neglected. The error of this approximation is negligible if the cuto distance is chosen to be R > 8 Bohr. The linear dependencies are removed for each individual pair domain separately, as described in detail in [128].

In the local basis, the first-order wavefunction takes the form
	(1)		1			˜ij	rs	˜ij	˜ji
|Ψ		=		ij P rs [ij]		Trs	|Φij	with Trs	= Tsr	(5.37)
			2

where P represents the pair list and it is implicitly assumed that the pair domains [ij] are defined as described above. Note that the number of projected functions r, s [ij] for a given pair (ij) is independent of the molecular size. Therefore, the

˜ij

individual amplitude matrices Trs are very compact and their sizes are independent

˜ij

of the molecular size. The total number of amplitudes Trs depends linearly on the molecular size and it is assumed that they can be stored in high-speed memory.

Since the local orbital basis does not diagonalize the zeroth-order HF Hamiltonian,

˜ij

an iterative procedure is required to determine the amplitudes T . By minimizing the Hylleraas functional, one obtains the linear equations [130]

R˜ij = K˜ ij + F˜T˜ij S˜ + ST˜ ˜ij F˜ − k		S˜ FikT˜kj + Fkj T˜ik S˜	(5.38)
˜ij		must vanish for r, s [ij]. The quan-
For the desired solution, the residuals R	rs

˜ ˜

tities S and F are the overlap and Fock matrices in the projected basis, respectively,

5.1 Electron Correlations in Molecules: Post-Hartree–Fock Methods

161

and the exchange matrices		˜ ij		= (ri|sj) represent a small subset of the trans-
		K	rs
				˜ ij	,
formed two-electron integrals. For a given pair (ij), only the local blocks K
˜	˜				rs

Frs, and Trs for r, s [ij] are needed in the first three terms, while for the over-

lap matrices in the sum only the blocks connecting the domain [ij] with [ik] or [jk] are required. The sizes of all these matrix blocks are independent of the molecular size. Taking further into account that for a given pair (ij) the number of terms k in the summation becomes asymptotically independent of the molecular size if very distant pairs are neglected, it follows that the computational e ort scales linearly with molecular size. Further savings are possible by skipping terms with coupling matrix elements Fik or Fkj falling below a certain threshold, which can be progressively reduced as the iteration proceeds. In the calculations [131] the initial threshold was set to 10−2, and reduced in each iteration by a constant factor until the final value of 10−5 was reached. This threshold yields the energy to an accuracy better than 10−7 Hartree. This screening reduces the iteration time, but is not essential for achieving linear-scaling. There are two possible algorithms to compute the terms in the sum.

˜	˜kj	˜
In the first algorithm (I) the matrix multiplications S	× T	× S are performed inside

the summation, which seems wasteful at first glance but has the advantage that all matrix multiplications are over small local blocks only. Alternatively, in the second algorithm (II) the matrix multiplications are performed outside the summation, but then the accumulated matrix in the square brackets has a much larger dimension, which is the union of all pair domains having either the orbital i or j in common. If very distant pairs are neglected, these united-pair domains (UP2) also become independent of molecular size, and therefore algorithm II scales linearly with molecular size as well. A comparison of the measured CPU times for both algorithms for linear chains of polyglycine peptides as a function of the number of monomers [131] shows that both algorithms closely approach linear-scaling. It was found that Algorithm I has a significantly lower prefactor than Agorithm II, and has therefore been chosen as the default in the MOLPRO program [116]. Note that both algorithms would lead to O(N 3) scaling if the very distant pairs were not neglected. Typically, 5–7 iterations are needed to converge the energy to better than 10−7 Hartree, the number of iterations has been found to be essentially independent of the molecular size.

The potential bottleneck is the following: the overlap and Fock matrices are stored in memory as full M × M matrices, where M is the number of AOs. The subblocks needed for a given matrix multiplication are extracted on the fly from these matrices. Full square matrices are also needed for computing and storing the projected orbitals. Thus, the memory requirement formally scales with O(M 2) but since only very few full matrices are needed, the prefactor is so small that this causes no bottleneck in practice for up to 2000–3000 basis functions. In principle, this bottleneck could be avoided by storing for each pair only the local subblocks needed. This would lead to strict linear-scaling of the memory (or disk) requirements.

Once the amplitudes have converged, the second-order correlation energy is obtained by tracing the amplitudes with the exchange matrices, i.e.

	2T˜rsij − T˜srij K˜rsij
E2 = ij P rs [ij]		(5.39)

162 5 Electron Correlations in Molecules and Crystals

The overall cost to solve the linear equations (5.38) and to compute the second-order energy (5.39) depends linearly on the molecular size and does not constitute any computational or storage bottleneck, once the necessary quantities are available. As in the canonical case, the computational e ort of a local MP2 calculation is entirely dominated by the integral transformation required to construct the exchange matrices

˜ ij . The linear scaling integral transformation is considered in detail in [131, 132].

K

The optimization of atomic structures requires the availability of analytic energy gradients. In the LMP2 method, the bottleneck is the calculation of the 4-index electron-repulsion integrals (ERIs) and their transformation from the atomic orbitals into the local orbital basis. The computational e ort for this step scales with the fourth power of the basis-set size (for fixed molecular size), and therefore the cost quickly increases with basis-set quality. The same is true for the computation of LMP2 (or MP2) energy gradients, which requires a reverse transformation of the amplitudes from the molecular orbitals into the atomic orbital basis. It has been demonstrated for LMP2 [146] and LCCSD [147] that the computational e ort for ERI evaluation and transformation can be reduced by 1–2 orders of magnitude using density-fitting (DF) methods. In the DF approach the one-electron charge densities in the ERIs, which are binary products of orbitals, are approximated by linear expansions in an auxiliary basis set. This leads to a decomposition of the 4-index ERIs in terms of 2- and 3-index ERIs, and the O(M 4) AO dependence of the computational cost is reduced to O(M 3). linear-scaling of the computational cost with molecular size in the context of DF can be achieved by employing the concept of locality also for the fitting functions. The key here is the use of local fitting domains, as discussed in [146, 147]. Furthermore, local-density fitting can be used to speed up the calculation of the exchange part to the Fock matrix in Hartree–Fock calculations [146]. The fitting errors introduced by the DF approximation are systematic and substantially smaller than any other typical errors of the calculation (like basis-set truncation or the error due to the local approximation), provided that a suitable auxiliary basis set is employed. Energy gradients for canonical DF-MP2 on top of an ordinary Hartree–Fock (HF) reference function, i.e. without invoking the DF approximation already at the HF level, were first implemented in [148]. An e cient method to compute analytical energy derivatives for the local MP2 approach was presented in [139] where the advantages of the local and density-fitting approximations in the calculation of analytic energy gradients were combined. All 4-index integrals and their derivatives were replaced by products of 2-index and 3-index quantities, which can be manipulated very efficiently. Moreover, the expensive backtransformation of the e ective second-order density matrix into the basis was avoided. As a consequence, the method has a much lower prefactor than conventional MP2 gradient programs, and the e ort scales only cubically with the basis-set size for a given molecule.

Many applications of local MP2 and local CCSD methods to molecular systems can be found in [132–134, 139, 142, 149].

Another linear-scaling MP2 algorithm was proposed in [129] that is based on the atomic-orbitals Laplace-transform (LT) MP2 method [137]. In this method, the energy denominators (5.32) in (5.30) are eliminated by Laplace-transformation, which paves the way to express the MP2 energy directly in the basis. The price to pay is the additional Laplace integration, which is carried out by quadrature over a few (8–10) points. For each of the quadrature points an integral transformation has to

5.2 Incremental Scheme for Local Correlation in Periodic Systems

163

be performed, but the transformation matrices are much sparser than the canonical MO coe cients; therefore, more e cient prescreening is possible. In this respect, the algorithm [137] is related to the local MP2 method. It appears, however, that due to the multiple transformations the LT-MP2 method has a higher prefactor than the LMP2 method, which requires only a single transformation. On the other hand, the LT-MP2 method does not involve any truncation of the virtual space and therefore yields the exact canonical MP2 energy within a certain accuracy, which is primarily determined by the number of quadrature points and the prescreening thresholds used [131]. An LT linear-scaling coupled-cluster (LTCC) method was also developed and implemented for closed [143] and open [138] shells.

In Sections 5.3 and 5.4 we consider Laplace-transform MP2 (LTMP2) and local MP2 (LMP2) methods extensions to periodic systems, implemented in computer codes GAUSSIAN03 [107] and CRYSCOR [117], respectively. These two extensions are similar to their molecular counterparts in the preliminary calculation of delocalized canonical orbitals (Bloch functions for periodic systems) and further construction of localized crystalline orbitals (Wannier functions for periodic systems). However the correlation e ects for periodic systems may be taken into account in another way when the SCF calculations are made in Wannier-functions basis. Such an approach was named the incremental scheme [5, 49, 110] and is considered in the next section.

5.2 Incremental Scheme for Local Correlation in Periodic Systems

5.2.1 Weak and Strong Electron-correlation

In both molecules and solids the electron-correlation e ects may be strong or weak. Let us consider an example of a H2 molecule. For the two electrons involved, Heitler and London [150] suggested the following correlated ground-state wavefunction

ΨHL(r1σ, r2σ) =	1	[χa(r1)χb(r2) + χb(r1)χa(r2)] [α1β2 − α2β1]	(5.40)
	N

where the AOs χa(r) and χb(r) are centered on atoms a and b of the molecule, and the spin functions α and β denote spin-up and spin-down states, respectively. One notes that this wavefunction does not contain ionic configurations in which two electrons are in the same atomic orbital χa(r) or χb(r). Therefore, their mutual Coulomb repulsion is minimal at the expense of the kinetic energy. The Heitler–London ansatz should work well when the Coulomb repulsion completely dominates the kinetic-energy gain caused by spreading an electron over both atoms. In this case we speak of the strong correlation limit and we expect it to be applicable, e.g. when we stretch the H–H bond length to large values.

A very di erent ansatz applies when the opposite limit holds, i.e. when the kineticenergy gain due to spreading an electron over both atoms dominates the mutual electron repulsions. The molecular-orbital theory describes the two electrons within the independent electron approximation and the bonding molecular orbital of an H2 molecule is

	ϕ(r) =	1	[χa(r) + χb(r)]	(5.41)
		N

164 5 Electron Correlations in Molecules and Crystals

and it is occupied by the two electrons with opposite spins. Thus, the antisymmetric wavefunction is

ΨM O(r1, σ, r2, σ) =	1	[χa(r1)χa(r2) + χb(r1)χb(r2)
	N
+χa(r1)χb(r2) + χb(r1)χa(r2)] [α1β2 − α2β1]			(5.42)

Contrary to (5.40), the wavefunction (5.42) contains, with considerable weight, ionic configurations in which both electrons are at one H site. These configurations imply a large Coulomb repulsion between the two electrons. Therefore, the MO LCAO approximation will be a reasonable description of the true or exact wavefunction only when the energy gain due to hybridization of the atomic wavefunction is large as compared with the Coulomb repulsions. In reality, the ground-state wavefunction of H2 has a form between (5.40) and (5.42). There are ionic configurations present, but their weight is partially suppressed as compared with (5.42). The partial suppression of charge fluctuations, i.e. a reduction of the mean-square deviations of the electronic charges is a hallmark of electronic correlations. Charge fluctuations imply configurations in which the electronic charges deviate from the average charge distribution. When in a configuration of the H2 molecule both electrons are on one atomic site so that the other site is empty, the deviation from the average charge is as large as possible. In fact, we may use the degree of suppression of charge fluctuations as a measure of the strength of electron correlations. A complete suppression of charge fluctuations, as in the Heitler–London wavefunction, implies the strong correlation limit. From the above it is obvious that we can learn much about electron correlations and their influence by studying the wavefunction of a system.

The solids with weak and with strong correlations require di erent methods of treatment. For example, when the correlations are weak, a wavefunction of independent electrons, i.e. a self-consistent field (SCF) or Hartree–Fock (HF) wavefunction seems a good starting point for the implementation of correlation corrections. In a solid, the molecular orbitals are replaced by crystalline (Bloch) orbitals. The latter are solutions of the SCF(HF) equations (4.57).

As we have pointed out above, a possible measure of the strength of electronic correlations is the degree of suppression of charge fluctuations in the true or correlated ground-state wavefunction as compared with those in a corresponding SCF or HF wavefunction. In order to obtain a feeling of how strong those suppressions are, we can calculate them for di erent chemical bonds [110]. Let Σ denote the ratio

					Σ =	ΦSCF |δn2|ΦSCF − Ψ0|δn2				|Ψ0	(5.43)
						ΦSCF |δn2|ΦSCF − ΨP C |δn2				|ΨP C
where δn2	= n	2	− n¯	2	and n¯ is the		average electron number in an orbital being part
							sp	3	orbital of a C–C σ-bond. The wavefunctions
of a bond, e.g. n¯ = 1 in a carbon

|ΦSCF , |Ψ0 and |ΨP C denote the HF ground state, the correlated or true ground state and the ground state we would have if the energy gain due to delocalization, i.e. moving an electron from atom to atom, were arbitrarily small. In |ΨP C the charge fluctuations are reduced by the largest amount that is compatible with the charge distribution, i.e. with the one-particle density matrix. A value of Σ = 0 implies no suppression of charge fluctuations and hence the limit of uncorrelated electrons, while Σ = 1 corresponds to the limit of strong electron correlations. For the H2

5.2 Incremental Scheme for Local Correlation in Periodic Systems

165

molecule Σ = 0.15, while for a C–C σ-bond Σ 0.2, indicating that electrons in diamond are weakly correlated. The same holds true for silicon. For a C–C π-bond instead one finds Σ 0.5 and correlations, e.g. in polyacetylene are fairly strong. When we consider the Cu–O planes of the perovskite-like La2CuO4, we find that for a Cu site ΣCu = 0.8, while for an O site ΣO = 0.7. Therefore, all Cu configurations di erent from 3d10 and 3d9 are essentially suppressed in |Ψ0 , while they still have an appreciable weight in |ΨSCF . More examples can be found in [5]. Reductions of electronic charge fluctuations as quantified by Σ can be described within a minimal basis set, e.g. a set of s and p functions for the valence electrons of atoms of the second row. They are caused by interatomic correlations, and we refer to the latter when we speak of weakly or strongly correlated electrons. They should be distinguished from intraatomic correlations that influence electrons on a given atomic site. They require a larger basis set for their description. At least one set of polarization functions must be included in this set, e.g. a set of d functions in the case of atoms of the second row. For more details on these di erent types of correlations we refer to reference [5]. In order to determine the correlated ground-state wavefunction |Ψ0 for the solids with relatively weak correlations one may start from the HF wavefunction |ΦSCF .

The wavefunction-based (post-Hartree–Fock) methods have received little attention in the electron-correlation theory for crystals compared with those for molecular systems [110]. The progress in the electron-correlation study in solids is connected mainly with density-functional theory (DFT), determining directly various physical properties of a system without a knowledge of the many-electron wavefunction. This feature makes DFT theory simpler and more attractive in the calculations of the electronic structure of crystals. Unfortunately, for the density-functional-based methods there is no procedure for systematic improvement of the calculated results when a higher accuracy is desired. Although often ground-state properties of unexpected accuracy are obtained, this is not always the case, in particular when electron correlations in a crystal are strong. Various improvements that have been applied, like self-interaction corrections (SIC) or LDA + U [151–153], show that they remain essentially uncontrolled. On the other hand, today we can calculate electronic groundstate wavefunctions for small molecules with high precision and there are controlled approximation schemes available for various degrees of accuracy. There is no reason why similar calculations cannot be done for solids; in fact it has been demonstrated that they can [110]. It was already mentioned in Chap. 4, that accurate HF solutions for a variety of crystalline systems have been available for many years now [4, 154]. The clever exploitation of the full symmetry of the crystal and implementation of special algorithms for the exact or approximate estimate of infinite lattice sums, permit the solution of big systems like faujasite (144 atoms, 1728 AOs per unit cell) [155] to be obtained at quite low computational cost. However, the limitations of the HF approximation are well known, essentially due to its neglecting dynamical electroncorrelation.

The practical methods of post-HF calculations for solids were discussed during the workshop “Local correlation methods: from molecules to crystals” [156]. Computational strategies were considered and new developments suggested in this area of research (the texts of invited speaker talks are published on an Internet site [156]). In particular, it was stated that the development of post-HF methods for crystals is essentially connected with the progress in the localized Wannier-function (LWF)

166 5 Electron Correlations in Molecules and Crystals

generation. In the last decade, LWF have been very intensively used in practical calculations of solids, see [157] and references therein. LWF for correlation calculations can be derived from Bloch states by a localization method as long as one does not have to deal with partially filled bands as in metals. There is also the possibility to generate LWFs directly from HF calculations, as is done in the computer code WANNIER [158].

In order to calculate the correlated ground state |Ψ0 and its energy E0, we split the many-electron Hamiltonian H into

H = HSCF + Hres

(5.44)

where HSCF is the SCF part of the Hamiltonian and Hres is the residual interaction part. We may write

	E0	=	ΦSCF |H|Ψ0	(5.45)
			ΦSCF |Ψ0

provided ΦSCF |Ψ0 = 0. This will be the case as long as we deal with finite systems. It is worth mentioning that instead of dividing H as in (5.44) we could have split

it also into a one-electron part HKS corresponding to the Kohn-Sham equation in density-functional theory (see Chap. 7) and a remaining part. This has the advantage that the energy of the related ground state ΨKS , a Slater determinant constructed from solutions of the Kohn–Sham equation, is much closer to the true ground-state energy than is E0. Such a scheme may appear useful for the solids with the strong electron correlations and is used, for example, in the SIC and LDA+U approaches mentioned above. But, treating the remaining part Hres would not be interpretative as an correlation-energy calculation in the sense that is traditional for molecular systems.

Treating a solid implies dealing with essentially infinitely many-electrons. However, in practice, we can correlate only a relatively small number of them. Therefore, we first have to reduce the calculations to a small number of electrons. This can be done by means of the method of increments considered in the next section.

5.2.2 Method of Increments: Ground State

A method of increments [111, 159, 160] is a wavefunction-based ab-initio correlation method for solids. This method is closely related to the ideas of the local ansatz (LA), [5] where local operators acting on the SCF wavefunction are used to admix suitable oneand two-particle excitations to the mean-field HF ground state. The many-electron Hamiltonian is split according to (5.44) and the ground-state Hamiltonian HSCF and the corresponding wavefunction ΦSCF = Φ0 are assumed to be known.

A product of two operators A and B in the Liouville space is defined as follows:

(A|B) = Φ0|A†B|Φ0 c = (A†B c

(5.46)

The superscript c indicates that the cumulant of the expectation value is taken, which is given by

(A c = A	(5.47)
AB c = (AB − A B	(5.48)

.

.

.