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

5.2 Incremental Scheme for Local Correlation in Periodic Systems

167

By using (5.46) the exact ground-state energy E is written in the following way:

E = (H|Ω) = EHF + (HI |Ω)

(5.49)

where Ω plays the role of the wave operator that describes the transformation from the HF ground-state to the exact ground-state. For a solid with well-defined bonds, one can express the HF ground-state Φ0 in terms of localized orbitals and label those orbitals by a bond index i. For ionic systems this index numbers the ions whose states form the occupied bands (for example, oxygen ions in MgO crystal).

The following operators are defined: Ai, where i should be considered as a compact index that includes the bond i as well as the oneand two-particle excitations of bond i, and Aij , which describes the two-particle excitations where one excitation is out of bond i while the other is out of bond j. Within the restricted operator subspace

spanned by Ai and Aij the operator Ω can be written in the form
|	Ω) =	1 +			niAi +	nij Aij		(5.50)
				i	ij
					(i	=j)
The parameters ni and nij in (5.50)				are determined from the set of equations
0 = (Ak|H) +							nij (Ak|HAij )
		i	ni(Ak|HAi) +			ij
						(i =j)
			ni(Akl|HAi) +			(i	nij (Akl|HAij )	(5.51)
0 = (Akl|H) +		i				ij
						=j)

where (5.50) was used.

The method of increments provides a scheme in which the set of equations (5.51) and hence the correlation energy is evaluated in a hierarchical order [159].

(a) First, all electrons are kept frozen except for those, e.g., in bond i. The operators Ai describe the corresponding excitations of these two electrons and (5.51) reduces to

0 = (Ai|H) + ni(1)(Ai|HAi)			(5.52)
Within this approximation, the n(1)	are independent of each other and the correlation
i
energy becomes	Ecorr(1) = i
with		i	(5.53)
i = ni(1)(HI |Ai)			(5.54)

(b) In the next step the electrons in two bonds, e.g., i and j are correlated. The corresponding n(2) parameters are determined from the coupled equations

0= (Ai|H) + n(2)i (Ai|HAi) + n(2)j (Ai|HAj ) + n(2)ij (Ai|HAij )

0= (Aj |H) + n(2)i (Aj |HAi) + n(2)j (Aj |HAj ) + n(2)ij (Aj |HAij )

0 = (Aij |H) + n(2)i (Aij |HAi) + n(2)j (Aij |HAj ) + n(2)ij (Aij |HAij ) (5.55)

168 5 Electron Correlations in Molecules and Crystals

Again, the increments δni = n(2)i − n(1)i and δnj = n(2)j − n(1)j are treated as independent of each other in this approximation, and we have

		1	(i
Ecorr(2) =	i +		∆ ij	(5.56)
	i	2	ij
			=j)
where	= ij − ( i + j )			(5.57)
∆ ij
and
ij = (HI |ni(2)Ai + nj(2)Aj + nij(2)Aij )				(5.58)

(c) Analogously, one calculates the three-bond energy increment, which is defined

as
∆ ijk = ijk − (∆ i + ∆ j + ∆ k) + (∆ ij + ∆ jk + ∆ ik)	(5.59)

The correlation energy ijk is that obtained when all electrons are kept frozen except those in bonds i, j, and k. Again, the increments ∆ ijk are treated as being independent of each other.

The total correlation energy within this approximation is the sum of all increments,

		1			1			(i
Ecorrsolid =	i +	2		∆ ij +		6		∆ ijk + · · ·	(5.60)
i			ij					ijk
			(i =j)					=j =k)

It is obvious that by calculating higher and higher increments the exact correlation energy within the coupled electron-pair approximation at level zero [161] is obtained.

The method of increments is useful only if the incremental expansion is well convergent, i.e. if increments up to, say, triples are su cient, and if increments become rapidly small with increasing distance between localized orbitals. These conditions were shown to be well met in the case of di erent solids [110, 162]. Ideally, the increments should be local entities not sensitive to the surroundings.

The theory described above has been applied to a great variety of materials, thus demonstrating the feasibility of calculations of that kind. They include the elemental semiconductors [159], III-V [163] and II-VI compounds [164], ionic crystals like MgO [165], CaO [166], NiO [167], alkali halides [168], TiO2 (with a sizeable amount of covalency) [169], rare-gas crystals [170, 171], solid mercury [172, 173] and the rareearth compound GdN [174] with the 4f electrons kept within the core. The method of increments allows the CCSD local correlation scheme to be extended from molecules to solids. In most cases the program package CRYSTAL [23] was used for the SCF part including a localization procedure for determining the Wannier functions.

An alternative approach to CRYSTAL is an embedded-cluster approach called WANNIER, where the localization procedure is part of the SCF calculations. Thereby, the solid is modeled by a central cluster embedded in a field created by the remaining part of the infinite solid [158]. Test calculations at the SCF level for LiH agree very well with those obtained from CRYSTAL. For example, the di erence in the SCF ground-state energy is only 0.05 eV per Li when a Li Wannier function extending over four and over three nearest-neighbors is used. Embedded-cluster calculations seem to be a very promising scheme for the future in order to perform high-quality

5.2 Incremental Scheme for Local Correlation in Periodic Systems

169

electronic-structure calculations not only for solids, but also for very large molecules. In the CeO2 coupled-cluster calculations [175] several embeddings were tested: from purely point charges up to pseudopotential – surrounding of oxygens to imitate the Ce ions. Embedded large nearly neutral clusters (Ce4O7) as well as several individual clusters, each of which consists of only correlated parts of the solid, were used in the calculations. Analysis of the obtained results allowed the conclusion that all tested embedding schemes should result in the same value of correlation energy. At the same time a large neutral cluster seems to be the best choice for the embedding.

A third way of obtaining the SCF ground-state wavefunction is by means of an ordinary cluster calculation. In this case, a fragment of the crystal lattice is used. Its dangling bonds are saturated with hydrogen atoms when dealing with covalent solids and it is surrounded by a large number of point charges when an ionic solid is treated. The wavefunction of the innermost bond or ion can be used to a good approximation for every bond or ion of the solid. An example is X35H36 (X = C, Si, Ge, Sn), where the localized orbital for the central bond agrees with that obtained from CRYSTAL, after a localization procedure has been applied [110].

5.2.3 Method of Increments: Valence-band Structure

and Bandgap

The valence-band energies for solids can be regarded as energy di erences between the N electron ground-state and a state with (N −1) electrons, where one delocalized Bloch electron has been removed. The bandgap can be regarded as the energy di erence between the N -electron ground state and the state with (N + 1) electrons, where one delocalized Bloch electron has been added. Fixing the number of electrons N in the ground-state means that PBC are introduced for the main region of a crystal, i.e. that the finite size system is considered.

At the SCF level, Koopmans’ theorem holds, and the task of determining the band energies and energy gap can be reduced to finding the one-electron energies of the respective Bloch states. It is well known that SCF (HF) calculations tend to overestimate the energy gap, giving too low an energy for the top of the upper valence band and too high an energy of the conduction-band bottom. DFT, on the other hand, yields gaps that are far too small when using the common functionals such as LDA or GGA. To overcome this problem, several DFT-based or hybrid HF-DFT approaches were developed (see Chap. 7). As was already mentioned, DFT-based approaches do not provide clear concepts of systematically improving accuracy and the numerical results essentially depend on the density-functional chosen.

The incremental scheme based on the wavefunction HF method was extended to the calculation of valence-band energies when the electron-correlation is taken into account. In [176, 177] an e ective Hamiltonian for the (N − 1)-electron system was set up in terms of local matrix elements derived from multireference configurationinteraction (MRCI) calculations for finite clusters. This allowed correlation corrections to a HF band structure to be expressed and reliable results obtained for the valenceband structure of covalent semiconductors. A related method based on an e ective Hamiltonian in localized Wannier-type orbitals has also been proposed and applied to polymers [178, 179]. Later, the incremental scheme was used to estimate the relative energies of valence-band states and also yield absolute positions of such states [180].

170 5 Electron Correlations in Molecules and Crystals

This is important as a prerequisite for the determination of bandgaps. Unfortunately, a further extension of the formalism to also cover conduction-band states (which would complete the information on bandgaps) is not possible along the same lines.

This is because an e ective Hamiltonian for (N + 1)-particle states cannot be extracted from finite-cluster calculations (which do not lead to stable states for the extra electron). Instead, a scheme for the calculation of localized states in infinite periodic surroundings is necessary [112, 181]. However, the information on the correlationinduced shift of the absolute position of the upper valence-band edge, already allows for an estimate of correlation e ects on bandgaps. In a reasonable approximation, one can assume that correlation e ects for the process of removing one electron, and the inverse process of adding an extra electron, are almost symmetric. This approximation seems plausible in view of the fact that the dominant e ect of correlations caused by an added electron (hole) is the generation of a long-ranged polarization cloud. The latter moves with the added particle and together with it forms a quasiparticle. A simple estimate of the corresponding gap correction is obtained by calculating the classical polarization-energy gain in a continuum approximation. This is given in [5]:

2				·		−	2 0	lc
δE =	1		d3rP		E =		0 −	1 e2	(5.61)

where P is the macroscopic polarization, E is the electric field generated by the extra particle and 0 is the dielectric constant of the medium. The cuto lc is approximately equal to the length at which the dielectric function, (r), reaches its asymptotic value,0, as a function of the distance, r, from the extra particle. The contribution, δE is approximately the same for an added electron and an added hole because lc should be nearly equal in both cases. Therefore, this assumption can be considered as physically justified, and indeed final estimates for the correlation-induced reduction of bandgaps compare reasonably well with experiment [180].

The incremental scheme for the valence band rests on the observation that a HFhole can be either described in reciprocal or in direct space. In reciprocal space, the (N -1)-electron states are introduced

|kνσ = akνσ|Φscf

(5.62)

where akνσ annihilates an electron with wave vector k and spin σ in band ν of the HFdeterminant |Φscf of the neutral system in the ground state (the spin index will be suppressed as the spinless HF Hamiltonian is supposed to be considered). According to Koopmans’ Theorem, its HF one-electron energies can be written in the form

εkscfν = E0scf − kν|H|kν

(5.63)

where E0scf is the HF-ground-state energy.

The states |kν can be approximated in terms of basis functions adapted to trans-

lational symmetry, with creation – annihilation operators bk† n, bkn:
akν =		(5.64)
	dνnbkn

n

The coe cients, dνn are found by diagonalizing the matrix

5.2 Incremental Scheme for Local Correlation in Periodic Systems	171
Hnnscf (k) = kn|E0scf − H|kn	(5.65)

where |kn = bkn|Φscf are assumed to be orthogonal, for simplicity.

To base considerations on states with localized extra electrons or holes, the Wannier-type transformation is introduced, yielding localized states within each unit cell. The number of these states equals the number of bands included in the Wannier transformation. In the case of elemental semiconductors such as diamond or silicon these are one-electron states localized at the middle of the nearest atoms bonds so that n=1,2,3,4 (see Chap. 3 for the discussion of the connection between Bloch and Wannier states). A local description arises through the decomposition of the bkn into a sum of local operators marked with a cell index R (translation vector of the direct

lattice):
		N		−
			R
bkn =	√1		exp(		ikR)cRn	(5.66)

The operator cRn annihilates an electron in a local basis function centered in cell R. The sum in (5.66) is precisely the projector onto the irreducible representation k of the translation group. In this basis the matrix (5.65) becomes

Hnnscf (k) = exp(−ikR) 0n|E0scf − H|Rn	(5.67)

R

which can be obtained by the direct-space quantities 0n|H|Rn . Without loss of generality, the functions annihilated by the cRn can be chosen as localized orbitals spanning the occupied HF space of the neutral system, i.e. as localized two-center bond orbitals for covalent semiconductors. In this case, 0n|H|Rn represents the hopping matrix element between a HF-hole located in the central cell in bond n to a bond n in cell R. The direct lattice formalism allows for inclusion of local correlation e ects. Such e ects are obtained when the HF states |Rn are replaced by their correlated counterparts |Rn }, and E0scf is replaced by the true ground-state energy,

E0.

An incremental scheme application to valence bands requires computing the correlated nondiagonal matrix elements 0n|H|Rn in terms of a rapidly converging series of approximations. Figure (5.1) shows the correlated valence-band structure for Si with and without correlation, taken from [177]. The energy at the top of the valence band in both cases has been set to zero. The local matrix elements were extracted from calculations of a set of small molecules XnHm and used to set up an incremental expansion of the bulk band structure.

For computing the absolute shift of the valence bands one needs to consider also the diagonal matrix element

{0n|H|0n }

(5.68)

Such a calculation is discussed in detail in [180]. First, a HF-hole is considered in a specific localized bond orbital Bi of a finite cluster simulating the crystalline semi-

conductor, with B		i representing |0i .
scf
Let i	denote the HF-value of a hole in this particular bond Bi:
		iscf = Bi|E0scf − H|Bi	(5.69)

172 5 Electron Correlations in Molecules and Crystals

Energy, eV

5.0
0.0
- 5.0
-10.0
-15.0
-20.0 L		X

Fig. 5.1. The SCF (dashed lines) and correlated (solid lines) band structure for silicon, [177].

As a first step correlation calculations are performed both for the neutral cluster and for the cluster with a hole in bond Bi by restricting electron excitations to excitations out of this bond only. The energy di erence, intrai between the correlated energies of the neutral and the charged cluster is then an excitation energy improved with respect to the HF-value (5.69) by the inclusion of correlation e ects. The di erence,

∆iin = iintra − iscf

(5.70)

constitutes the “intrabond” increment of the correlation corrections to the diagonal element. In the next step, correlations are taken into account by releasing electrons out of two bonds, namely the one that contains the hole, Bi, and one additional bond Bj . Again, a reference calculation is needed for the neutral cluster, where electrons are excited out of these two bonds. Let us denote the corresponding energy with Ei0;j . For the case with a hole, a multireference configuration-interaction (MRCI) calculation is perfomed where the two HF-determinants that contain the HF-hole in bond Bi and bond Bj , respectively, are the reference functions forming the model space. The calculation results in a diagonal (2 × 2)-eigenenergy matrix with eigenvalues λ1 and λ2 and corresponding eigenvectors |Ψ1 and |Ψ2 . Due to the diagonalization in the course of the correlation treatment the two eigenvectors are delocalized over the region spanned by the two bonds considered. To get the matrix element of the Hamiltonian between the local states, as in (5.69), the diagonal representation of H given by the eigenenergies is rotated back into a local matrix representation. To do so, the eigenvectors are projected onto the model space to give the reduced eigenvector

matrix, D :		Bj |Ψ1	Bj |Ψ2
	D =	Bi|Ψ1	Bi|Ψ2		(5.71)

5.2 Incremental Scheme for Local Correlation in Periodic Systems

173

Due to this projection, the two eigenvectors contained in (5.71) are not strictly orthonormal, but can be made so by means of L¨owdin symmetric orthonormalization

( → ) ). The matrix

D D

M = D	Ei0;j − λ1	0	0		D†	(5.72)
)	0	Ei;j	−	λ2	)

gives the matrix representation of H in local states. The diagonal element of this matrix corresponding to the bond Bi is the correlated value of one-electron energy (5.69) taking into account the e ect of correlation of one additional bond Bj . The procedure can be extended to include further increments. Double-bond increments contain the e ect of considering two additional bonds Bj , Bk, while the triple-bond increments include three further bonds and Bj , Bk, Bl.

The mainly local character of the correlation e ects manifests itself in a rapid decrease of the magnitude of the increments, both with regard to the number of bonds and distances between them. This allows the numerical e ort for the calculation to be reduced to just a few increments, which can be evaluated for a finite (embedded) cluster. In [180] it is demonstrated that for diamond and silicon crystals the largest correlation correction comes from the single excitations of the increment called SB1 that takes into account 66% and 64% of the overall correlation e ect, respectively. Thus the polarization of the first coordination sphere, stabilizing the ionized system, gives by far the most important separate correlation e ect.

The bonds Bi used in the formulation of the incremental scheme above, can be obtained by a HF-calculation for such a finite cluster with an a-posteriori Boys localization [38], which directly yields the bonds as localized HF-orbitals. All increments were derived in [180] from the cluster calculations. C and Si clusters dangling bonds were suppressed by the introduction of hydrogen atoms. The electron correlations dramatically reduce the diagonal elements from their HF eigenvalues –0.68, –0.51 hartree to –0.45, –0.34 Hartree, for C and Si, respectively. When these values of the diagonal matrix element are fed into a band-structure calculation the absolute shift of the valence bands due to electronic correlations is obtained. The above results lead to a shift of the upper valence band at the Γ -point of 3.90 eV for C and 2.63 eV for Si. To compare these results with experiment, we note that in diamond the HF-gap at the Γ -point has been calculated as 13.8 eV [182], while the experimental one is 7.3 eV. If one attributes the reduction of 6.5 eV in equal parts to the shift of the conduction and the valence band, each shift should amount to 3.3 eV, so we overestimate the shift by 0.6 eV. In the case of Si the (direct) HF-gap at the Γ -point of 8.4 eV [182], is too large by 5.0 eV, as compared to the corresponding experimental value of 3.4 eV, which would attribute a share of 2.5 eV to the valence band shift. The result obtained in [180] of 2.6 eV nearly coincides with this shift. The deviation of these estimates for diamond from the experimental value can be explained by the crude assumption of a symmetric contribution, both for the valence and for the conduction band, to the closing of the direct HF-gap.

Di erent approaches to the problem of excitations in correlated electron systems are considered in [183]. They are based on a quasiparticle description when electron correlations are weak and on a Green’s function or projection operator approach when they are strong. In both cases, intersite correlation contributions to the energy bands require special attention.

174 5 Electron Correlations in Molecules and Crystals

It is shown that for both small and large correlation-energy contributions the method of increments reduces the amount of computations decisively.

In [162] a critical review of the incremental scheme of obtaining the energetic properties of extended systems from wavefunction-based ab-initio calculations of small (embedded) building blocks is given, which has been applied to a variety of van der Waals-bound (solid Xe), ionic (MgO), and covalent (C and Si) solids in the past few years. Its accuracy is assessed by means of model calculations for finite systems, and the prospects for applying it to delocalized systems are given. The individual determination of the increments of the many-body expansion leads, on the other hand, to the possibility of determining these increments in small finite model systems where a suitable embedding makes the orbitals of a given building block (or pairs, triples, etc. of building blocks) similar to those of the extended one. This e ciently reduces the computational e ort and obviates the necessity of restricting excitations to the virtual space (by introducing a domain structure) as in the local-correlation methods. This also means that the many-body expansion approach can be used with every quantum-chemical program package.

The correlated band structure of the one-periodic systems (polymers) was e - ciently studied using the AO Laplace-transformed MP2 theory discussed in next section.

5.3 Atomic Orbital Laplace-transformed MP2 Theory for Periodic Systems

5.3.1 Laplace MP2 for Periodic Systems: Unit-cell Correlation Energy

In [184] explicit expressions were presented for electron-correlation at the MP2 level and implemented for the total energy per unit cell and for the band structure of extended systems. Using MP2 for EN , EN −1, EN +1 (the total energies of one determinant states for N , (N − 1)- and (N + 1)-electron systems) a formula was presented for a direct evaluation of the bandgap rather than obtaining it as a di erence of two large numbers. The formulation was given in the conventional crystalline orbitals delocalized over the whole extended system. Later [185] this theory for periodic systems was reformulated. The unit cell MP2 energy was evaluated by Laplace transforming the energy denominator and heavily screening negligible contributions in the AO formulation. This is the extension to periodic systems of the AO Laplace-transformed formulation applied earlier for molecules in [186, 187]. This Laplace MP2 real-space approach avoids the high computational cost of multidimensional wavevector k integrations that is critical in order to obtain accurate estimates. It is e cient to treat one-, two and three-periodic crystals (1D, 2D, and 3D systems, respectively). The benchmark calculations [185, 188] were carried out for polymers (1D) and hexagonal BN plane (2D) as well as stacked 1D polymers. The HF portion of calculations was carried out using a PBC option of the molecular computer code Gaussian03 [107].

Here, we briefly present the main points of Laplace MP2 theory for periodic systems.

Equation (5.30) for the correlation energy in the MP2 approximation for molecules can be written in the form

5.3 Atomic Orbital Laplace-transformed MP2 Theory

for Periodic Systems

175

Occ,V ir (ia|jb)[2(ia|jb) − (ib|ja)]

E2 = − (5.73)

ij,ab

a + b − i − j

where i and a are HF energies of the canonical occupied and virtual MOs, respectively. The two-electron integrals in the canonical MO basis, (ia|jb), are obtained by transformation of the Coulomb two-electron integrals in the atom-centered Gaussiantype orbital basis,

(ia|jb) =	1		a(1)b(2)dr1dr2
	i(1)j(2)
		r12
=	(µν|λσ)CµiCνaCλj Cσb			(5.74)

µνλσ

For periodic systems the MP2 correlation energy correction is given by an expression very similar to the molecular case [184]; the only di erence is the use of the crystalline Bloch (canonical) orbitals instead of MOs:

	1
E2 = −Re			dk1dk2dk3dk4 (i(k1)a(k3)|j(k2)b(k4))
		Vk4
×	2(i(k1)a(k3) j(k2)b(k4)) − (i(k1)b(k4)|j(k2)b(k3))			(5.75)
			a(k| 3) + b(k4) − i(k1) − j (k2)

The generalization to a spin unrestricted formalism is done in a similar fashion as in the molecular case.

The two-electron integrals in the energy formula are now complex and given by

(i(k1)a(k3)|j(k2)b(k4)) =
	(µpνq|λrσs)CµiCνaCλj Cσb
	pqrs µνλσ
× exp [i(−k1p + k3q − k2r + k4s) · g]		(5.76)

Here, the subscripts p, q, r, s indicate which cell the particular AO belongs to. Provided that the energy denominator is a strictly positive quantity, the Laplace-

transform can be used to write

1
a(k3) + b(k4) − i(k1) − j (k2)
= exp [(− a(k3) − b(k4))t] exp [( i(k1) + j (k2))t] dt	(5.77)

The central feature of the transformed energy denominator is that the indices and attached quantities that were coupled by the denominator are now decoupled, allowing a reordering of the summations. In particular, the four di erent k-integrations can be carried out independently of each other, yielding the following unit-cell MP2 correlation energy formula (where repeated greek indices are contracted):

E2 =		νpσs	(t) [2(µ0νp|λrσs) − (µ0σs|λrνp)] dt
	prs Tµ0λr
			=	νpσs	(t)vνµp0σλsr dt	(5.78)
				Tµ0λr

prs

176 5 Electron Correlations in Molecules and Crystals

where the amplitudes T of the Coulomb integrals transformed by the Fouriertransformed weighted density and complimentary matrices X and Y , respectively:

νpσs

Tµ0λr (t) = Xµ0γt Yνpδu (γtδu|κv w)Xλrκv Yσs w

= (

µ

0

p|λr

s)

(5.79)

where

ν

σ

Xµpγs =

1

C(k)µj exp(− j (k)t)C(k)γj

Vk

j

× exp [ik · (p − s)] dk,

1

Yµpγs =

C(k)µa exp( a(k)t)C(k)γa

Vk

a

× exp [ik · (p − s)] dk

(5.80)

In matrix elements (5.80) the summation runs over the occupied states (j) and virtual states (a), respectively. C(k)µj and C(k)µa are the expansion coe cients on the atomic basis for states j and a at point k. Vk is the volume of the first Brillouin zone, where integration over k is performed.

By comparing the expression (5.78) for the MP2 correlation energy to the canonical reciprocal-space expression (5.75), it becomes clear that the multidimensional quasimomentum integration has been reduced to a series of independent Fourier transforms (5.80). For all practical purposes, the computational cost of calculating the MP2 correlation energy in the AO basis is independent of the number of k-points used in the Brillouin-zone integration. Very dense k-point meshes can be used at little extra cost and thus, convergence issues found in earlier [184] reciprocal-space MP2 implementations are avoided. As was shown in [189] the interactions between atoms in the central cell (cell 0) and atoms in the neighboring cells can be quantified by considering the atoms the basis functions in (5.78) are centered on. In the case of periodic MP2 theory, it is particularly interesting to consider a cell-pair partitioning,

				E2 =	e20r	(5.81)
					r
0r	=		νpσs	(t) [2(µ0	νp|λrσs) − (µ0σs|λrνp)] dt	(5.82)
e2		ps Tµ0λr

Whereas for e02r the contributions from cells p and s decay exponentially with intercell distance between 0 and r, the correlation contribution from e02r itself decays quickly, with respect to |o − r| distance.

It has been shown for molecules [189] that the expected 1/R6 decay of distant correlation contributions is closely modeled by the atom-pair partitioning. The same type of decay is also observed in the periodic case and e02r decays as 1/|o − r|6. Of course, the extent of the lattice summations over prs in (5.78) greatly a ects both the computational cost and the accuracy of the calculation. By taking advantage of the spatial decay properties mentioned above, the energy contraction can be carried out only for cells such that |o − r| ≤ rmax, while the p and s lattice summation truncations must be controlled [185].