7.2 Theoretical Background

Detailed description of methods implemented in the riper module is provided in Refs [9598] and references therein. Here, only a short summary of the underlying theory is provided.

7.2.1 Kohn-Sham DFT for Molecular and Periodic Systems

In periodic systems 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

 k        1   ∑   ikTL ∑   k
ψpσ(r) = √N----   e       CμpσμL(r)
           UC  L       μ

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 unrestricted Kohn-Sham equations

  k k    k  k k
F σCσ = S C σεσ,

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. Equation (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 = LeikTLF μνσ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,
FLμνσ = TμLν + JLμν + XLμνσ.

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

E = T + J + E  .

7.2.2 RI-CFMM Approach

The key component of riper is a combination of RI approximation and CFMM applied for the electronic Coulomb term [95,98]. In the RI scheme the total crystal electron density ρcryst is approximated by an auxiliary crystal electron density ˜ρ cryst

ρcryst ≈ ˜ρcryst =  ˜ρL,

composed of unit cell auxiliary densities ˜ρ L with

     ∑   T
˜ρL =    c αL,

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

    ∫∫         1    ∑       ′    ′  ∑             ∑
D =    δρ (r)|-r--r′ |  δρL(r)drdr =    (δρ | δρL) =  (ρ- ˜ρ | ρL - ˜ρL).
                     L               L             L

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

JLμν =   (μνL | ˜ρL′ - ρnL′),

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

with the real space density matrix elements obtained by integration

DLμνσ = 1--   DkμνσeikTLdk,
        Vk BZ

of the reciprocal space density matrix

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

over the BZ with volume V k.

Equations (7.9) and (7.10) as well as other expressions appearing in the RI scheme require calculation of infinite lattice sums of the form

   (ρ1 | ρ2L),

where the distribution ρ1 in the central cell interacts with an infinite number of distributions ρ2L, i.e., ρ2 translated by all possible direct lattice vectors L. In the RI-CFMM scheme [95] the sum in 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 based on their extents. Interactions between charges from well-separated boxes are calculated using a hierarchy of multipole expansions. 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. The remaining contribution to the Coulomb term is obtained from direct integration. This approach results in nearly linear scaling of the computational effort with the system size.

7.2.3 k Point Sampling Scheme

The integral in Eq. (7.11) is evaluated approximately using a set of sampling points k. riper uses a Γ-centered mesh of k points with weights wk, so that Eq. (7.11) can be written as

       ∑       T
DLμνσ ≈    wkeik LDkμνσ.

In 3D periodic systems each sampling point is defined by its components k1, k2 and k3 along the reciprocal lattice vectors b1, b2 and b3 as

k = k1b1 + k2b2 + k3b3.

For 2D periodic systems k3 = 0. In case of 1D periodicity k3 = 0 and k2 = 0. In riper the three components kj (j = 1,2,3) of k are given as

k  = - nj---1,- nj---1+ 1,..., nj---1- 1, nj---1.
 j     2nj     2nj           2nj      2nj

with nj (j = 1,2,3) as integer numbers. riper reduces the number of k points employed in actual calculation by a factor of two using time-inversion symmetry, i.e., the vectors k and -k are symmetry equivalent. The k point mesh can be specified providing the integer values nj within the data group $kpoints.

The number of k points required in a calculation critically depends on required accuracy. Generally, metallic systems require considerably more k points than insulators to reach the same precision. For metals, the number of k point also depends on parameters of the Gaussian smearing [99] used in riper. Please refer to Ref. [99] for more details.

7.2.4 Metals and Semiconductors: Gaussian Smearing

Achieving reasonable accuracy of DFT calculations for metals requires a higher number of k points than for semiconductors and insulators. The convergence with respect to the number of k points can be improved applying partial occupancies [99]. To achieve this, riper uses the Gaussian smearing method in which occupation numbers fnk are calculated as

   (       )     (      [       ])
     ϵnk --μ    1         ϵnk---μ-
fnk     σ    =  2  1- erf    σ      ,

where ϵnk are band (orbital) energies, μ is the Fermi energy and σ is the width of the smearing.

When smearing is applied the total energy E has to be replaced a generalized free energy F

F = E - ∑  σS(f  )
        nk     nk

in order to obtain variational functional. riper output file reports the values of F as “FREE ENERGY” and the term - nkσS(fnk) as “T*S”. In addition, the value of E for σ 0 is given as “ENERGY (sigma->0)”.

Gaussian smearing can be switched on for riper calculations by simply providing the value of σ > 0 within the $riper data group using the keyword sigma, e.g., sigma 0.01. The value of σ should be as large as possible, but small enough to yield negligible value of the “T*S” term. Please refer to Ref. [99] for a more detailed discussion.

The use of Gaussian smearing often requires much higher damping and orbital shifting. Please adjust the values for $scfdamp and $scforbitalshift if you encounter SCF convergence problems.

The optional keyword desnue can be used within the $riper data group to constrain the number of unpaired electrons. This can be used to force a certain multiplicity in case of an unrestricted calculation, e.g., desnue 0 for singlet and desnue 1 for dublet.

7.2.5 Low-Memory Iterative Density Fitting Method

For calculations on very large molecular systems a low-memory modification of the RI approximation has been implemented within the riper module [97]. In the RI approximation minimization of the Coulomb repulsion of the residual density, Eq. (7.8), yields a system of linear equations

Vc = γ,

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

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

In the LMIDF approach a conjugate gradient (CG) iterative method is used for solution of Eq. (7.19). In order to decrease the number of CG iterations a preconditioning is employed, i.e., Eq. (7.19) is transformed using a preconditioner P to an equivalent problem

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

with an improved condition number resulting in faster convergence of the CG method. The iterative CG solver in riper employs one of the following preconditioners that are formed from 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ββ   = (α | ˜ρ).
        β                 β

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.