# Direct Inversion of the Iterative Subspace

When solving systems of linear (or nonlinear) equations, iterative methods are often employed.  Unfortunately, such methods often suffer from convergence issues such as numerical instability, slow convergence, and significant computational expense when applied to difficult problems.  In these cases, convergence accelleration methods may be applied to both speed up, stabilize and/or reduce the cost for the convergence patterns of these methods, so that solving such problems become computationally tractable.  One such method is known as the direct inversion of the iterative subspace (DIIS) method, which is commonly applied to address convergence issues within self consistent field computations in Hartree-Fock theory (and other iterative electronic structure methods).  In this tutorial, we'll introduce the theory of DIIS for a general iterative procedure, before integrating DIIS into our previous implementation of RHF.

## I. Theory

DIIS is a widely applicable convergence acceleration method, which is applicable to numerous problems in linear algebra and the computational sciences, as well as quantum chemistry in particular.  Therefore, we will introduce the theory of this method in the general sense, before seeking to apply it to SCF.  

Suppose that for a given problem, there exist a set of trial vectors $\{\mid{\bf p}_i\,\rangle\}$ which have been generated iteratively, converging toward the true solution, $\mid{\bf p}^f\,\rangle$.  Then the true solution can be approximately constructed as a linear combination of the trial vectors,
$$\mid{\bf p}\,\rangle = \sum_ic_i\mid{\bf p}_i\,\rangle,$$
where we require that the residual vector 
$$\mid{\bf r}\,\rangle = \sum_ic_i\mid{\bf r}_i\,\rangle\,;\;\;\; \mid{\bf r}_i\,\rangle 
=\, \mid{\bf p}_{i+1}\,\rangle - \mid{\bf p}_i\,\rangle$$
is a least-squares approximate to the zero vector, according to the constraint
$$\sum_i c_i = 1.$$
This constraint on the expansion coefficients can be seen by noting that each trial function ${\bf p}_i$ may be represented as an error vector applied to the true solution, $\mid{\bf p}^f\,\rangle + \mid{\bf e}_i\,\rangle$.  Then
\begin{align}
\mid{\bf p}\,\rangle &= \sum_ic_i\mid{\bf p}_i\,\rangle\\
&= \sum_i c_i(\mid{\bf p}^f\,\rangle + \mid{\bf e}_i\,\rangle)\\
&= \mid{\bf p}^f\,\rangle\sum_i c_i + \sum_i c_i\mid{\bf e}_i\,\rangle
\end{align}
Convergence results in a minimization of the error (causing the second term to vanish); for the DIIS solution vector $\mid{\bf p}\,\rangle$ and the true solution vector $\mid{\bf p}^f\,\rangle$ to be equal, it must be that $\sum_i c_i = 1$.  We satisfy our condition for the residual vector by minimizing its norm,
$$\langle\,{\bf r}\mid{\bf r}\,\rangle = \sum_{ij} c_i^* c_j \langle\,{\bf r}_i\mid{\bf r}_j\,\rangle,$$
using Lagrange's method of undetermined coefficients subject to the constraint on $\{c_i\}$:
$${\cal L} = {\bf c}^{\dagger}{\bf Bc} - \lambda\left(1 - \sum_i c_i\right)$$
where $B_{ij} = \langle {\bf r}_i\mid {\bf r}_j\rangle$ is the matrix of residual vector overlaps.  Minimization of the Lagrangian with respect to the coefficient $c_k$ yields (for real values)
\begin{align}
\frac{\partial{\cal L}}{\partial c_k} = 0 &= \sum_j c_jB_{jk} + \sum_i c_iB_{ik} - \lambda\\
&= 2\sum_ic_iB_{ik} - \lambda
\end{align}
which has matrix representation
\begin{equation}
\begin{pmatrix}
  B_{11} & B_{12} & \cdots & B_{1m} & -1 \\
  B_{21} & B_{22} & \cdots & B_{2m} & -1 \\
  \vdots  & \vdots  & \ddots & \vdots  & \vdots \\
  B_{n1} & B_{n2} & \cdots & B_{nm} & -1 \\
  -1 & -1 & \cdots & -1 & 0
\end{pmatrix}
\begin{pmatrix}
c_1\\
c_2\\
\vdots \\
c_n\\
\lambda
\end{pmatrix}
=
\begin{pmatrix}
0\\
0\\
\vdots\\
0\\
-1
\end{pmatrix},
\end{equation}

which we will refer to as the Pulay equation, named after the inventor of DIIS.  It is worth noting at this point that our trial vectors, residual vectors, and solution vector may in fact be tensors of arbitrary rank; it is for this reason that we have used the generic notation of Dirac in the above discussion to denote the inner product between such objects.

## II. Algorithms for DIIS
The general DIIS procedure, as described above, has the following structure during each iteration:
#### Algorithm 1: Generic DIIS procedure
1. Compute new trial vector, $\mid{\bf p}_{i+1}\,\rangle$, append to list of trial vectors
2. Compute new residual vector, $\mid{\bf r}_{i+1}\,\rangle$, append to list of trial vectors
3. Check convergence criteria
    - If RMSD of $\mid{\bf r}_{i+1}\,\rangle$ sufficiently small, and
    - If change in DIIS solution vector $\mid{\bf p}\,\rangle$ sufficiently small, break
4. Build **B** matrix from previous residual vectors
5. Solve Pulay equation for coefficients $\{c_i\}$
6. Compute DIIS solution vector $\mid{\bf p}\,\rangle$

For SCF iteration, the most common choice of trial vector is the Fock matrix **F**; this choice has the advantage over other potential choices (e.g., the density matrix **D**) of **F** not being idempotent, so that it may benefit from extrapolation.  The residual vector is commonly chosen to be the orbital gradient in the AO basis,
$$g_{\mu\nu} = ({\bf FDS} - {\bf SDF})_{\mu\nu},$$
however the better choice (which we will make in our implementation!) is to orthogonormalize the basis of the gradient with the inverse overlap metric ${\bf A} = {\bf S}^{-1/2}$:
$$r_{\mu\nu} = ({\bf A}^{\rm T}({\bf FDS} - {\bf SDF}){\bf A})_{\mu\nu}.$$
Therefore, the SCF-specific DIIS procedure (integrated into the SCF iteration algorithm) will be:
#### Algorithm 2: DIIS within an SCF Iteration
1. Compute **F**, append to list of previous trial vectors
2. Compute AO orbital gradient **r**, append to list of previous residual vectors
3. Compute RHF energy
3. Check convergence criteria
    - If RMSD of **r** sufficiently small, and
    - If change in SCF energy sufficiently small, break
4. Build **B** matrix from previous AO gradient vectors
5. Solve Pulay equation for coefficients $\{c_i\}$
6. Compute DIIS solution vector **F_DIIS** from $\{c_i\}$ and previous trial vectors
7. Compute new orbital guess with **F_DIIS**

## III. Implementation

In order to implement DIIS, we're going to integrate it into an existing RHF program.  Since we just-so-happened to write such a program in the last tutorial, let's re-use the part of the code before the SCF integration which won't change when we include DIIS:


```python
# ==> Basic Setup <==
# Import statements
import psi4
import numpy as np

# Memory specification
psi4.set_memory(int(5e8))
numpy_memory = 2

# Set output file
psi4.core.set_output_file('output.dat', False)

# Define Physicist's water -- don't forget C1 symmetry!
mol = psi4.geometry("""
O
H 1 1.1
H 1 1.1 2 104
symmetry c1
""")

# Set computation options
psi4.set_options({'basis': 'cc-pvdz',
                  'scf_type': 'pk',
                  'e_convergence': 1e-8})

# Maximum SCF iterations
MAXITER = 40
# Energy convergence criterion
E_conv = 1.0e-6
D_conv = 1.0e-3
```


```python
# ==> Static 1e- & 2e- Properties <==
# Class instantiation
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option('basis'))
mints = psi4.core.MintsHelper(wfn.basisset())

# Overlap matrix
S = np.asarray(mints.ao_overlap())

# Number of basis Functions & doubly occupied orbitals
nbf = S.shape[0]
ndocc = wfn.nalpha()

print('Number of occupied orbitals: %d' % ndocc)
print('Number of basis functions: %d' % nbf)

# Memory check for ERI tensor
I_size = (nbf**4) * 8.e-9
print('\nSize of the ERI tensor will be %4.2f GB.' % I_size)
if I_size > numpy_memory:
    psi4.core.clean()
    raise Exception("Estimated memory utilization (%4.2f GB) exceeds allotted memory \
                     limit of %4.2f GB." % (I_size, numpy_memory))

# Build ERI Tensor
I = np.asarray(mints.ao_eri())

# Build core Hamiltonian
T = np.asarray(mints.ao_kinetic())
V = np.asarray(mints.ao_potential())
H = T + V
```

    Number of occupied orbitals: 5
    Number of basis functions: 24
    
    Size of the ERI tensor will be 0.00 GB.



```python
# ==> CORE Guess <==
# AO Orthogonalization Matrix
A = mints.ao_overlap()
A.power(-0.5, 1.e-16)
A = np.asarray(A)

# Transformed Fock matrix
F_p = A.dot(H).dot(A)

# Diagonalize F_p for eigenvalues & eigenvectors with NumPy
e, C_p = np.linalg.eigh(F_p)

# Transform C_p back into AO basis
C = A.dot(C_p)

# Grab occupied orbitals
C_occ = C[:, :ndocc]

# Build density matrix from occupied orbitals
D = np.einsum('pi,qi->pq', C_occ, C_occ, optimize=True)

# Nuclear Repulsion Energy
E_nuc = mol.nuclear_repulsion_energy()
```

Now let's put DIIS into action.  Before our iterations begin, we'll need to create empty lists to hold our previous residual vectors (AO orbital gradients) and trial vectors (previous Fock matrices), along with setting starting values for our SCF energy and previous energy:


```python
# ==> Pre-Iteration Setup <==
# SCF & Previous Energy
SCF_E = 0.0
E_old = 0.0
```

Now we're ready to write our SCF iterations according to Algorithm 2.  Here are some hints which may help you along the way:

#### Starting DIIS
Since DIIS builds the approximate solution vector $\mid{\bf p}\,\rangle$ as a linear combination of the previous trial vectors $\{\mid{\bf p}_i\,\rangle\}$, there's no need to perform DIIS on the first SCF iteration, since there's only one trial vector for DIIS to use!

#### Building **B**
1. The **B** matrix in the Lagrange equation is really $\tilde{\bf B} = \begin{pmatrix} {\bf B} & -1\\ -1 & 0\end{pmatrix}$.
2. Since **B** is the matrix of residual overlaps, it will be a square matrix of dimension equal to the number of residual vectors.  If **B** is an $N\times N$ matrix, how big is $\tilde{\bf B}$?
3. Since our residuals are real, **B** will be a symmetric matrix.
4. To build $\tilde{\bf B}$, make an empty array of the appropriate dimension, then use array indexing to set the values of the elements.

#### Solving the Pulay equation
1. Use built-in NumPy functionality to make your life easier.
2. The solution vector for the Pulay equation is $\tilde{\bf c} = \begin{pmatrix} {\bf c}\\ \lambda\end{pmatrix}$, where $\lambda$ is the Lagrange multiplier, and the right hand side is $\begin{pmatrix} {\bf 0}\\ -1\end{pmatrix}$.  


```python
# Start from fresh orbitals
F_p =  A.dot(H).dot(A)
e, C_p = np.linalg.eigh(F_p)
C = A.dot(C_p)
C_occ = C[:, :ndocc]
D = np.einsum('pi,qi->pq', C_occ, C_occ, optimize=True)

# Trial & Residual Vector Lists
F_list = []
DIIS_RESID = []

# ==> SCF Iterations w/ DIIS <==
print('==> Starting SCF Iterations <==\n')

# Begin Iterations
for scf_iter in range(1, MAXITER + 1):
    # Build Fock matrix
    J = np.einsum('pqrs,rs->pq', I, D, optimize=True)
    K = np.einsum('prqs,rs->pq', I, D, optimize=True)
    F = H + 2*J - K
    
    # Build DIIS Residual
    diis_r = A.dot(F.dot(D).dot(S) - S.dot(D).dot(F)).dot(A)
    
    # Append trial & residual vectors to lists
    F_list.append(F)
    DIIS_RESID.append(diis_r)
    
    # Compute RHF energy
    SCF_E = np.einsum('pq,pq->', (H + F), D, optimize=True) + E_nuc
    dE = SCF_E - E_old
    dRMS = np.mean(diis_r**2)**0.5
    print('SCF Iteration %3d: Energy = %4.16f dE = % 1.5E dRMS = %1.5E' % (scf_iter, SCF_E, dE, dRMS))
    
    # SCF Converged?
    if (abs(dE) < E_conv) and (dRMS < D_conv):
        break
    E_old = SCF_E
    
    if scf_iter >= 2:
        # Build B matrix
        B_dim = len(F_list) + 1
        B = np.empty((B_dim, B_dim))
        B[-1, :] = -1
        B[:, -1] = -1
        B[-1, -1] = 0
        for i in range(len(F_list)):
            for j in range(len(F_list)):
                B[i, j] = np.einsum('ij,ij->', DIIS_RESID[i], DIIS_RESID[j], optimize=True)

        # Build RHS of Pulay equation 
        rhs = np.zeros((B_dim))
        rhs[-1] = -1
        
        # Solve Pulay equation for c_i's with NumPy
        coeff = np.linalg.solve(B, rhs)
        
        # Build DIIS Fock matrix
        F = np.zeros_like(F)
        for x in range(coeff.shape[0] - 1):
            F += coeff[x] * F_list[x]
    
    # Compute new orbital guess with DIIS Fock matrix
    F_p =  A.dot(F).dot(A)
    e, C_p = np.linalg.eigh(F_p)
    C = A.dot(C_p)
    C_occ = C[:, :ndocc]
    D = np.einsum('pi,qi->pq', C_occ, C_occ, optimize=True)
    
    # MAXITER exceeded?
    if (scf_iter == MAXITER):
        psi4.core.clean()
        raise Exception("Maximum number of SCF iterations exceeded.")

# Post iterations
print('\nSCF converged.')
print('Final RHF Energy: %.8f [Eh]' % SCF_E)
```

    ==> Starting SCF Iterations <==
    
    SCF Iteration   1: Energy = -68.9800327333871337 dE = -6.89800E+01 dRMS = 1.16551E-01
    SCF Iteration   2: Energy = -69.6472544393141675 dE = -6.67222E-01 dRMS = 1.07430E-01
    SCF Iteration   3: Energy = -75.7919291462249021 dE = -6.14467E+00 dRMS = 2.89274E-02
    SCF Iteration   4: Energy = -75.9721892296711019 dE = -1.80260E-01 dRMS = 7.56446E-03
    SCF Iteration   5: Energy = -75.9893690602363563 dE = -1.71798E-02 dRMS = 8.74982E-04
    SCF Iteration   6: Energy = -75.9897163367029691 dE = -3.47276E-04 dRMS = 5.35606E-04
    SCF Iteration   7: Energy = -75.9897932415930200 dE = -7.69049E-05 dRMS = 6.21200E-05
    SCF Iteration   8: Energy = -75.9897956274068633 dE = -2.38581E-06 dRMS = 2.57879E-05
    SCF Iteration   9: Energy = -75.9897957845313670 dE = -1.57125E-07 dRMS = 1.72817E-06
    
    SCF converged.
    Final RHF Energy: -75.98979578 [Eh]


Congratulations! You've written your very own Restricted Hartree-Fock program with DIIS convergence accelleration!  Finally, let's check your final RHF energy against <span style='font-variant: small-caps'> Psi4</span>:


```python
# Compare to Psi4
SCF_E_psi = psi4.energy('SCF')
psi4.compare_values(SCF_E_psi, SCF_E, 6, 'SCF Energy')
```

    	SCF Energy........................................................PASSED





    True



## References
1. P. Pulay. *Chem. Phys. Lett.* **73**, 393-398 (1980)
2. C. David Sherrill. *"Some comments on accellerating convergence of iterative sequences using direct inversion of the iterative subspace (DIIS)".* Available at: vergil.chemistry.gatech.edu/notes/diis/diis.pdf. (1998)

03bcb8fffae1af959427f5b5fa9d42500276db0c
ipynb
Jupyter Notebook
Tutorials/03_Hartree-Fock/3b_rhf-diis.ipynb
andyj10224/psi4numpy
cbef6ddcb32ccfbf773befea6dc4aaae2b428776
BSD-3-Clause
2017-03-01T08:04:48.000Z
2022-03-23T08:52:04.000Z
Tutorials/03_Hartree-Fock/3b_rhf-diis.ipynb
andyj10224/psi4numpy
cbef6ddcb32ccfbf773befea6dc4aaae2b428776
BSD-3-Clause
2017-03-03T13:20:20.000Z
2022-03-05T18:20:27.000Z
Tutorials/03_Hartree-Fock/3b_rhf-diis.ipynb
andyj10224/psi4numpy
cbef6ddcb32ccfbf773befea6dc4aaae2b428776
BSD-3-Clause
2017-02-17T19:44:47.000Z
2022-03-22T05:52:43.000Z
Qwen/Qwen-72B
1. YES
2. YES