7.2 Theoretical Background

The theory behind riper and implementation details are briefly presented as far as needed for applications and understanding of the main keywords. For more details see Refs [9799] and references therein.

In periodic systems the translational symmetry of solids leads to Bloch orbitals ψk and one-particle energies εk depending on the band index p, spin σ, and the wave vector k within the Brillouin zone (BZ), which is the unit cell of reciprocal space. The orbitals

           1  ∑      ∑
ψkpσ(r) = √----    eikL   CkμpσμL (r)
          NUC  L      μ
(7.1)

are expanded in GTO basis functions μ(r - Rμ - L) μL(r) centered at atomic positions Rμ in direct lattice cells L over all NUC unit cells. This results in the unrestricted Kohn-Sham equations,

FkσCkσ = SkCkσεkσ,
(7.2)

which may be solved separately for each k in the BZ. The same equations hold for the molecular case, where only L = k = 0 is a valid choice and NUC is one. Eq 7.2 contains the reciprocal space Kohn-Sham and the overlap matrices Fσk and Sk, respectively, obtained as Fourier transforms of real space matrices

Fμνσk = LeikTL FμνσL   S μνk = LeikTLS μνL. (7.3)
The elements FμνσL contain three contributions: elements TμνL of the kinetic energy matrix, elements JμνL of the Coulomb matrix, and elements XμνσL of the exchange-correlation matrix,
 L      L    L     L
Fμνσ = Tμν + Jμν + Xμνσ.
(7.4)

The total energy per unit cell E is calculated as the sum of the kinetic T, Coulomb J, and exchange-correlation contribution EXC,

E = T + J + E    .
             XC
(7.5)

Calculations of kinetic energy and corresponding matrix elements are similar as in the molecular case. For the exchange-correlation term an adaptive numerical integration scheme is used [98]. The core of this method is a hierarchical spatial grouping of basis functions based on their spatial extents using an octree.

For the Coulomb term the periodic Resolution of Identity approximation is applied [100]. In this approach the total crystal electron density ρcryst is approximated by auxiliary crystal electron density ˜ρ cryst

 cryst   cryst  ∑
ρ    ≈ ˜ρ    =    ˜ρL,
               L
(7.6)

composed of unit cell auxiliary densities ˜ρ L with

     ∑   T
˜ρL =    c αL,
      α
(7.7)

where α denotes the vector of auxiliary basis functions. The vector of expansion coefficients c is determined by minimizing the Coulomb repulsion D of the residual density δρ = ρ -ρ˜

    ∑             ∑
D =    (δρ | δρL) =  (ρ - ˜ρ | ρL - ˜ρL).
     L             L
(7.8)

The RI approximation allows to replace four-center electron repulsion integrals (ERIs) by two- and three-center ones. In this formalism, elements of Coulomb matrix JμνL are defined as

 L′  ∑
Jμν =   (μνL′ | ˜ρL - ρnL),
      L
(7.9)

where ρn denotes the unit cell nuclear charge distribution. The total Coulomb energy including the nuclear contribution is

J = ∑  DL ′JL′- 1 ∑  (˜ρ+ ρ |ρ˜ - ρ  ).
    μνL′ μν μν  2  L      n  L    nL
(7.10)

with the real space density matrix elements obtained by integration

          ∫
DLμνσ = 1--   DkμνσeikTLdk,
        Vk BZ
(7.11)

of the reciprocal space density matrix

       ∑
Dkμνσ =    fkpσ (Ckμpσ)*Ckνpσ
        p
(7.12)

over the BZ with volume V k. The reciprocal space integral in Eq. 7.11 is evaluated numerically on a finite grid of k-points.

Eqs. 7.9 and 7.10 as well as RI scheme require evaluating infinite lattice sums of the form

∑
   (ρ1 | ρ2L),
 L
(7.13)

where the distribution ρ1 in the central cell interacts with an infinite number of distributions ρ2L, i.e., ρ2 translated by all possible L. An important property of the RI-CFMM formulation in riper is that these sums are convergent. In the RI-CFMM scheme [86] the sum from Eq. 7.13 is partitioned into Crystal Far-Field (CFF) and Crystal Near-Field (CNF) parts. The CFF part contains summation over all direct space lattice vectors L for which the overlap between the distributions ρ1 and ρ2L is negligible. This part is very efficiently calculated using multipole expansions. The CNF contribution is evaluated using an octree based algorithm. In short, a cubic parent box enclosing all distribution centers of ρ1 and ρ2 is constructed that is large enough to yield a predefined number ntarg of distribution centers per lowest level box. The parent box is successively subdivided in half along all Cartesian axes yielding the octree. In the next step all charge distributions comprising ρ1 and ρ2 are sorted into boxes basing on their extents so that distributions from well-separated boxes do not overlap. Two boxes are considered well-separated if the distance between their centers is greater than sum of their lengths times 0.5×wsicl, where wsicl is a predefined parameter 2. Interactions between charges from well-separated boxes are calculated using the hierarchy of multipole expansions. The remaining contribution to the Coulomb term is obtained from direct integration. This approach results in linear scaling of the computational effort with the system size.

For the molecular case, a low-memory modification of the RI approximation is implemented within the riper module [99]. In the RI scheme minimization of the Coulomb repulsion of the residual density, Eq. 7.8, yields a system of linear equations

Vc = γ,
(7.14)

where V is the Coulomb metric matrix with elements V αβ = (α | β) representing Coulomb interaction between auxiliary basis functions and vector γ is defined as

     ∑
γα =    (α | μν)D μν.
      μν
(7.15)

In the new approach instead of attempting a direct solution of this inhomogeneous system of linear equations, the conjugate gradient (CG) iterative method is used. In order to decrease the number of CG iterations a preconditioning is employed, i.e., the original system of Eq. 7.14, is transformed using a preconditioner P to an equivalent problem

(P-1V )c = P-1γ
(7.16)

with an improved condition number resulting in faster convergence of the CG method. Iterative CG solver employs preconditioners that are formed from the blocks of the V matrix

corresponding to the strongest and most important interactions between the auxiliary basis functions such that P-1V I.

The costly matrix-vector products of the Vc type that need to be evaluated in each CG iteration are not calculated directly. Instead, the linear scaling CFMM implementation presented above is applied to carry out this multiplication since the elements of the Vc vector represent Coulomb interaction between auxiliary basis functions α and an auxiliary density ˜ρ

                     (         )
(Vc) =  ∑  (α | β)cβ = (α | ∑ cββ) = (α | ˜ρ).
    α   β                 β
(7.17)

Hence, in contrast to conventional RI neither the V matrix nor its Cholesky factors need to be stored and thus significant memory savings are achieved.