03 Hartree-Fock Closed-Shell SCF Procedure¶

The purpose of this project is to provide a deeper understanding of (restricted) Hartree-Fock theory by demonstrating a simple implementation of the self-consistent-field (SCF) method. The theoretical background can be found in Ch. 3 of the text by Szabo and Ostlund or in the nice set of on-line notes written by David Sherrill. A concise set of instructions for this project may be found here.

Original authors (Crawford, et al.) thank Dr. Yukio Yamaguchi of the University of Georgia for the original version of this project.

This project could be a challenging one. Although it is somehow easy to achieve each steps, the number of steps is quite large. You may well possibly have done all tasks, but confuse on how and why iterative SCF could be written in this way. You are encouraged to review this project once you've done all steps.

The test case used in the following discussion is for a water molecule with a bond-length of $1.1 {\buildrel_{\circ} \over{\mathsf{A}}}$ and a bond angle of $104.0 {}^\circ$ with an STO-3G basis set. The input to the project consists of the nuclear repulsion energy ({download}input/h2o/STO-3G/enuc.dat) and pre-computed sets of one- and two-electron integrals: overlap integrals ({download}input/h2o/STO-3G/s.dat), kinetic-energy integrals ({download}input/h2o/STO-3G/t.dat), nuclear-attraction integrals ({download}input/h2o/STO-3G/v.dat), electron-electron repulsion integrals ({download}input/h2o/STO-3G/eri.dat).

{admonition}
:class: warning

From this project, it could be multiple ways to achieve a task. What *multiple* means here is that one could
- either read integrals from input file, or from PySCF's integral interface;
- either do (high-dimension) tensor contraction by extremely powerful `numpy.einsum` ([NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.einsum.html)), or by simply for-looping, or hybrid for-looping and 2-dimension matrix manipulation;
- format outputs by his own style (with most important values unchanged to the reference solution output);
- etc.

Multiple approaches could arrive to the almost same solution, and every approach is appreciated. The reader may develop his/hers own coding style. However, as author of the documentation, there's still some partiality for using PySCF's integral interface and using `numpy.einsum`.

In this project, using module pyscf.scf is prohibited. Using pyscf.gto is permitted. The structure of this project is slightly changed compared to the original project 03.

{admonition}
:class: warning, dropdown

It is important to make notation convention clear at the very beginning, or otherwise one may could have taken the wrong comprehension all the time reading article or documentation even without recognizing that. Here we make some essential, but not all, notations that we use in the documents thereafter. Meaning of some notations to some readers maybe not clear or familiar. These notations would be discussed in later sections or documents.

- $\mu, \nu, \kappa, \lambda$: Atomic orbitals. These notations are different to Szabo and Ostlund's $\mu, \nu, \sigma, \lambda$. Decision of changing $\sigma$ to $\kappa$ based on the convention that $\sigma$ is more abused to be notation of orbital spin.
- $i, j, k, l, ...$: Occupied orbitals. These notations are more accepted nowadays, which are different to Szabo and Ostlund's $a, b, c, ...$.
- $a, b, c, d, ...$: Virtual orbitals (un-occupied orbitals). The same reason to above.
- $p, q, r, s, m$: All accepted orbitals.
- $n_\mathrm{occ}, n_\mathrm{vir}, n_\mathrm{AO}, n_\mathrm{MO}$: Number of occupied, virtual, atomic and molecular orbitals, respectively.
- $A, B, C, M$: Atom subscripts. These notations are different to the previous projects' $i, j, k, ...$
- $\boldsymbol{r}, \boldsymbol{R}_A$: Coordinate of electron and nucleus of atom $A$, respectively.
- $Z_A$: Atomic charge of atom $A$.

In the documentation, we will not learn frozen orbitals. So division of occupied and virtual is sufficient. 

Molecule Object Initialization¶

In this project, we may use the updated Molecule initialization. Since these are technical details, we toggle the code and explanation below.

{toggle}

The meaning of attributes of `Molecule` are:

- `atom_charges`: Array to store atomic charges
- `atom_coords`: Matrix to store coordinates of atoms
- `natm`: Number of atoms
- `mol`: `pyscf.gto.Mole` instance, which describes molecular and basis information; *only used in PySCF approaches*
- `path_dict`: Dictionary of file path which containing integral, energy, ...; *only used in reading-input approaches*
- `nao`: Number of atomic orbitals $n_\mathrm{AO}$
- `charge`: Molecular charge ($\mathsf{NH_4^+}$ contains +1 charge and $\mathsf{NO_3^-}$ contains -1 charge); in most cases we just use neutral molecules
- `nocc`: Number of occupied orbitals $n_\mathrm{occ}$

Readers may be familier with the first three attributes which we have extensively used in [Project 01](../Project_01/Project_01.ipynb) and [Project 02](../Project_02/Project_02.ipynb). One may call `Molecule.construct_from_dat_file` to initialize these three attributes. For demonstration, we use STO-3G water molecule which is instantiated as variable `mole`.

{admonition}
:class: dropdown

In PySCF, molecule and basis information is stored in `pyscf.gto.Mole` ([PySCF API](http://pyscf.org/pyscf/gto.html)) instance. For example, the water molecule and basis STO-3G (the same to de demonstrated molecule) could be instantiated by

```python
mol = gto.Mole()
mol.unit = "Bohr"
mol.atom = """
8   0.000000000000  -0.143225816552   0.000000000000
1   1.638036840407   1.136548822547  -0.000000000000
1  -1.638036840407   1.136548822547  -0.000000000000
"""
mol.basis = "STO-3G"
mol.build()
```

Though the molecule can be manually instantiated, we still need a method function in `Molecule`, using existing `mole.atom_charges` and `mole.atom_coords` and external basis `"STO-3G"`, to build the `pyscf.gto.Mole` instance and store in `mole.mol` attribute.

Note that we just discuss closed-shell situation, so we set `mol.spin = 0`. The convention of molecular spin in different softwares might be quite deviated. In PySCF, `mol.spin` means number of alpha electrons minus beta electrons.

Class `pyscf.gto.Mole` can be extremely useful and convenient for this documentation, as almost all electron integrals can be called by method function `pyscf.gto.Mole.intor` in only one short line of code.

Implementation

Reader should fill all NotImplementedError in the following code:

Solution

{admonition}
:class: dropdown

`Molecule.charge` and `Molecule.path_dict` should be expicitly defined by caller. If you decided to use PySCF generating integrals, and not using pre-generated integrals in directory `input/h2o/STO-3G`, you are prefered to set this value to `None` or `NotImplemented`. Otherwise, you may explicitly state paths manually as the following toggled code:

Step 1: Nuclear Repulsion Energy¶

Nuclear Repulsion Energy can be determined as (Szabo, eq 3.185)

$$ E_\mathrm{nuc} = \sum_{B > A} \sum_{A} \frac{Z_A Z_B}{R_{AB}} $$

{admonition}
:class: dropdown

PySCF provides a method function that directly calculate nuclear repulsion energy `pyscf.gto.Mole.energy_nuc`.

{admonition}
:class: dropdown

File {download}`input/h2o/STO-3G/enuc.dat` contains nuclear repulsion energy.

{admonition}
:class: dropdown

Actually programming nuclear repulsion energy is not a difficult job. One can simply write two for loops to achieve the task. However, there could be a bonus challenge, only using matrix/tensor operation without for loop, to calculate the energy. The point is

$$
E_\mathrm{nuc} = \frac{1}{2} \sum_{AB} \frac{Z_A Z_B}{R_{AB}}
$$

where we (*unphysically*) define $R_{AB} = +\infty$ if $A = B$.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

Step 2: One-Electron Integrals¶

Read various atomic orbital (AO) electron integrals.

Overlap s

$$ S_{\mu \nu} = \int \phi_\mu (\boldsymbol{r}) \phi_\nu (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} $$

Kinetic Integral t

$$ T_{\mu \nu} = \int \phi_\mu (\boldsymbol{r}) \left( - \frac{1}{2} \nabla_{\boldsymbol{r}}^2 \right) \phi_\nu (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} $$

Nuclear-Attraction Integral v

$$ V_{\mu \nu} = \int \phi_\mu (\boldsymbol{r}) \left( - \sum_A \frac{Z_A}{|\boldsymbol{r} - \boldsymbol{R}_A|} \right) \phi_\nu (\boldsymbol{r}) \, \mathrm{d} \boldsymbol{r} $$

Then form the core Hamiltonian:

$$ H_{\mu \nu}^\mathrm{core} = T_{\mu \nu} + V_{\mu \nu} $$

{admonition}
:class: dropdown

It could be important to know how large the matrix is, before we handle electron integrals.

To get the dimension of matrix from fomular, for example overlap $S_{\mu \nu}$, we should take a look at subscripts $\mu \nu$, where it indicates dimension is related to $(\mu, \nu)$. Since $\mu, \nu$ are atomic orbitals, dimension of overlap should be $(n_\mathrm{AO}, n_\mathrm{AO})$.

However, $n_\mathrm{AO}$ can not be apparently got (since it requires summation of basis; however not everyone is required to recite how many basis functions are for Oxygen or Hydrogen DZP, or whether Uraine could be described by STO-3G).

For PySCF approach, one can use `pyscf.gto.Mole.nao` attribute to obtain $n_\mathrm{AO}$. For read-input approach, one could read the last line of file, for example in {download}`input/h2o/STO-3G/t.dat`:
```
    7     7    0.760031883566609
```
Then we should know there are 7 atomic orbitals, i.e. $n_\mathrm{AO} = 7$.

One may implement a short method function `obtain_nao` to store number of atomic orbitals.

{admonition}
:class: dropdown

Here we need to read three integrals. Since these three integrals are stored in exactly the same way, so we can make some kind of *abstraction*, i.e. a function which could not only read overlap, but also read kinetic or nuclear attraction integrals.

Recall that all these 1-electron integrals are symmetric matrices (or *permutationally unique*), so in .dat files, only lower-triangle part is stored. However, we prefer to store full matrices in computer memory for coding convenience. Also note that the AO indices on the integrals in the files start with 1 rather than 0.

One may implement a function `read_integral_1e` to read 1-electron integrals.

{admonition}
:class: dropdown

In PySCF, integrals can be generated by elegent function `pyscf.gto.Mole.intor`. For example, overlap integral can be generated in following code.

```python
>>> sol_mole.mol.intor("int1e_ovlp")
array([[ 1.       ,  0.2367039,  0.       , -0.       ,  0.       ,  0.0384056,  0.0384056],
       [ 0.2367039,  1.       ,  0.       ,  0.       ,  0.       ,  0.3861388,  0.3861388],
       [ 0.       ,  0.       ,  1.       ,  0.       ,  0.       ,  0.2684382, -0.2684382],
       [-0.       ,  0.       ,  0.       ,  1.       ,  0.       ,  0.2097269,  0.2097269],
       [ 0.       ,  0.       ,  0.       ,  0.       ,  1.       ,  0.       ,  0.       ],
       [ 0.0384056,  0.3861388,  0.2684382,  0.2097269,  0.       ,  1.       ,  0.1817599],
       [ 0.0384056,  0.3861388, -0.2684382,  0.2097269,  0.       ,  0.1817599,  1.       ]])
```

Full list of supported integrals are listd in [PySCF API](http://pyscf.org/pyscf/gto.html#pyscf.gto.moleintor.getints). What we need in this step is
- `int1e_ovlp`: overlap
- `int1e_kin`: kinetic integral
- `int1e_nuc`: nuclear-attraction integral

Implementation¶

You may either want to store the integrals by Molecule attributes (in computer memory), or generate integrals on-the-fly (i.e. generate integrals when we aquire them, otherwise they are not stored in computer memory). The inplementation of reference solution uses the latter way.

Reader should fill all NotImplementedError in the following code:

Reader should choose one approach and remove NotImplementedError in the following code:

Solution¶

Step 3: Two-Electron Integrals¶

Read the two-electron repulsion integrals (electron repulsion integral, ERI) from the {download}input/h2o/STO-3G/eri.dat file. The integrals in this file are provided in Mulliken notation over real AO basis functions:

$$ (\mu \nu | \kappa \lambda) = \int \phi_\mu (\boldsymbol{r}_1) \phi_\nu (\boldsymbol{r}_1) \frac{1}{|\boldsymbol{r}_1 - \boldsymbol{r}_2|} \phi_\kappa (\boldsymbol{r}_2) \phi_\lambda (\boldsymbol{r}_2) \, \mathrm{d} \boldsymbol{r}_1 \, \mathrm{d} \boldsymbol{r}_2 $$

Hence, the integrals obey the eight-fold permutational symmetry relationships:

$$ \begin{align} (\mu \nu | \kappa \lambda) &= (\mu \nu | \lambda \kappa) = (\nu \mu | \kappa \lambda) = (\nu \mu | \lambda \kappa) \\ = (\kappa \lambda | \mu \nu) &= (\lambda \kappa | \mu \nu) = (\kappa \lambda | \nu \mu) = (\lambda \kappa | \nu \mu) \end{align} $$

and only the permutationally unique integrals are provided in the file, with the restriction that, for each integral, the following relationships hold:

$$ \mu \geqslant \nu , \quad \kappa \geqslant \lambda \quad \mu \nu \geqslant \kappa \lambda $$

where

$$ \mu \nu := \mu (\mu + 1) / 2 + \nu \, \quad \kappa \lambda = \kappa (\kappa + 1) / 2 + \lambda $$

Although two-electron integrals may be stored efficiently in a one-dimensional array and the above relationship, we still use full dimension (i.e. dimension $(\mu, \nu, \kappa, \lambda)$ or $(n_\mathrm{AO}, n_\mathrm{AO}, n_\mathrm{AO}, n_\mathrm{AO})$) to store the two-electron integral tensor, for coding convenience.

If you prefer smaller memory buffer but somehow greater cost of computational efficiency cost and code complexity, original project 03 provides some insights on this.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Reader should choose one approach and remove NotImplementedError in the following code:

Solution¶

Here we only print one sub-matrix in the electron repulsion integral:

Step 4: Build the Orthogonalization Matrix¶

Diagonalize the overlap matrix:

$$ \mathbf{S} \mathbf{L} = \mathbf{L} \mathbf{\Lambda} $$

where $\mathbf{L}$ is the matrix of eigenvectors (columns) and $\mathbf{\Lambda}$ is the diagonal matrix of corresponding eigenvalues.

Build the symmetric orthogonalization matrix using (Szabo and Ostlund, eq 3.167):

$$ \mathbf{S}^{-1/2} = \mathbf{L} \mathbf{\Lambda}^{-1/2} \mathbf{L}^\dagger $$

where dagger $\dagger$ denotes the matrix transpose.

{admonition}
:class: dropdown

Actually this process is exactly the same to matrix power (with fractional exponent). This could be done simply by `scipy.linalg.fractional_matrix_power` ([SciPy API](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.fractional_matrix_power.html))

For proper definition of matrix power, see [Matrix function: Diagonalizable matrices](https://en.wikipedia.org/wiki/Matrix_function#Diagonalizable_matrices). Another useful definition could be [Cauchy integral](https://en.wikipedia.org/wiki/Matrix_function#Cauchy_integral) from complex analysis; actually this definition is applied to -1/2 power of matrix in direct random phase approximation (dRPA) correlation energy calculation. Well, in this project, diagonalizable matrix scheme should be sufficient.

{admonition}
:class: dropdown

For diagonal matrix, -1/2 power of matrix is equal to -1/2 power to its diagonal values. For example,

$$
\begin{pmatrix} 2 & 0 \\ 0 & 3 \end{pmatrix}^{-1/2} = \begin{pmatrix} 1/\sqrt{2} & 0 \\ 0 & 1/\sqrt{3} \end{pmatrix}
$$

So, not all matrices could be implied with -1/2 power. Only those matrices which eigenvalues are non-zero could be implied with -1/2 power. If negative eigenvalues exists, the matrix implied with -1/2 power could be complex.

{admonition}
:class: dropdown

Our cases in this project is somehow well-behaved. So, the overlap matrices $\mathbf{S}$ would have positive eigenvalues, and orthogonalization matrix is simply symmetric with dimension $(n_\mathrm{AO}, n_\mathrm{AO})$. Thus number of molecular orbital $n_\mathrm{MO}$ is equal to $n_\mathrm{AO}$.

However, when basis become larger, overlap matrices $\mathbf{S}$ could be ill-behaved (linear dependent, or some eigenvalues very close to zero), thus $\mathbf{S}^{-1/2}$ can be numerical instable. Solution to this behavior is discussed in Szabo and Ostlund roughly, around eq 3.172. In this situation, $n_\mathrm{MO} < n_\mathrm{AO}$.

It is important to know that, though in this documentation, $n_\mathrm{MO} = n_\mathrm{AO}$ in most situations, they are essentially two different concepts and does not need to be equal to each other. Molecular orbitals is linear combinition of atomic orbitals (LCAO). Computationally or say in practice, these two values are related by orthogonalization matrix.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

Step 5: Build Fock Matrix¶

Fock matrix of restricted RHF could be calculated as (Szabo and Ostlund, eq 3.154)

$$ F_{\mu \nu} = H_{\mu \nu}^\mathrm{core} + \sum_{\kappa \lambda} D_{\kappa \lambda} \big[ (\mu \nu | \kappa \lambda) - \frac{1}{2} (\mu \kappa | \nu \lambda) \big] $$

where the double-summation runs over all the atomic orbitals.

Fock matrix could be regarded as function of density matrix $D_{\kappa \lambda}$. For reproducible and arbitrariness consideration, we may use an arbitary density matrix dm_ramdom in the toggled code (notice that density matrix should be symmetry):

{admonition}
:class: dropdown

There are multiple ways to achieve Coulomb integral

$$
\sum_{\kappa \lambda} D_{\kappa \lambda} (\mu \nu | \kappa \lambda)
$$

One simple way could be for-loop:
```python
eri = mole.integral_eri()
coulomb = np.zeros((mole.nao, mole.nao))
for mu in range(mole.nao):
    for nu in range(mole.nao):
        coulomb[mu, nu] = (eri[mu, nu] * dm_random).sum()
```

Another way is `numpy.einsum` ([NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.einsum.html), [Einstein Summation Convention](https://en.wikipedia.org/wiki/Einstein_notation)):
```python
eri = mole.integral_eri()
coulomb = np.einsum("uvkl, kl -> uv", eri, dm_random)
```

A third way is utilizing numpy boardcasting ([NumPy Docs](https://numpy.org/doc/stable/user/theory.broadcasting.html)):
```python
eri = mole.integral_eri()
coulomb = (eri * dm_random).sum(axis=(-1, -2))
```

Result of the above variable `coulomb` should be
```python
array([[23.9330989,  4.234324 ,  0.0829256,  0.202269 ,  0.3299715,  0.7005698,  0.6955185],
       [ 4.234324 ,  8.9051099,  0.1324784,  0.4010859,  0.9716293,  2.6311948,  2.6251234],
       [ 0.0829256,  0.1324784,  9.0270751, -0.1065442,  0.0349767,  1.596199 , -1.5653815],
       [ 0.202269 ,  0.4010859, -0.1065442,  9.0874751, -0.1440637,  1.2958253,  1.3189783],
       [ 0.3299715,  0.9716293,  0.0349767, -0.1440637,  8.5657791,  0.2010907,  0.1926869],
       [ 0.7005698,  2.6311948,  1.596199 ,  1.2958253,  0.2010907,  3.7626551,  0.8928712],
       [ 0.6955185,  2.6251234, -1.5653815,  1.3189783,  0.1926869,  0.8928712,  3.7152703]])
```

{admonition}
:class: dropdown

There are also multiple ways to achieve Exchange integral

$$
\sum_{\kappa \lambda} D_{\kappa \lambda} (\mu \kappa | \nu \lambda)
$$

One simple way could be for-loop:
```python
eri = mole.integral_eri()
exchange = np.zeros((mole.nao, mole.nao))
for mu in range(mole.nao):
    for nu in range(mole.nao):
        exchange[mu, nu] = (eri[mu, :, nu] * dm_random).sum()
```

Another way is `numpy.einsum` ([NumPy API](https://numpy.org/doc/stable/reference/generated/numpy.einsum.html), [Einstein Summation Convention](https://en.wikipedia.org/wiki/Einstein_notation)):
```python
eri = mole.integral_eri()
exchange = np.einsum("ukvl, kl -> uv", eri, dm_random)
```

A third way is utilizing numpy boardcasting ([NumPy Docs](https://numpy.org/doc/stable/user/theory.broadcasting.html)):
```python
eri = mole.integral_eri()
exchange = (eri * dm_random[:, None, :]).sum(axis=(-1, -3))
```

Result of the above variable `exchange` should be
```python
array([[17.1408607,  2.9642433,  1.6954651,  3.4316066,  4.5482841,  1.0611853,  0.5458695],
       [ 2.9642433,  2.8747645,  0.352244 ,  1.8035899,  2.7594346,  1.3116692,  1.3224965],
       [ 1.6954651,  0.352244 ,  4.1142193, -1.0576397,  0.2781629,  0.6138329, -1.6697351],
       [ 3.4316066,  1.8035899, -1.0576397,  4.0318516, -1.39224  ,  0.5225643,  0.8203703],
       [ 4.5482841,  2.7594346,  0.2781629, -1.39224  , -0.8581923, -0.220426 , -0.137229 ],
       [ 1.0611853,  1.3116692,  0.6138329,  0.5225643, -0.220426 , -0.2555552,  0.0139439],
       [ 0.5458695,  1.3224965, -1.6697351,  0.8203703, -0.137229 ,  0.0139439,  0.006697 ]])
```

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

Step 6: Orbital Coefficients¶

Orbital coefficients $C_{\mu p}$ can be obtained from Fock matrix diagonalization (Roothaan equation, Szabo and Ostlund, eq 3.138-139):

$$ \mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C} \boldsymbol{\varepsilon} $$

or in element of matrix form:

$$ \sum_\nu F_{\mu \nu} C_{\nu i} = \varepsilon_i \sum_v S_{\mu \nu} C_{\nu i} $$

{admonition}
:class: dropdown

The conventional approach is described in Szabo and Ostlund, eq 3.174 - 3.178. First, we need to obtain transformed Fock matrix $\mathbf{F}'$ :

$$
\mathbf{F}' = \mathbf{S}^{-1/2} \mathbf{F} \mathbf{S}^{-1/2}
$$

Then obtain transformed orbital coefficient $\mathbf{C}'$ from the following transformed Roothaan equation:

$$
\mathbf{F}' \mathbf{C}' = \mathbf{C}' \boldsymbol{\varepsilon}
$$

Then we may get the orbital coefficient $\mathbf{C}$ by:

$$
\mathbf{C} = \mathbf{S}^{-1/2} \mathbf{C}'
$$

{admonition}
:class: dropdown

Actually, SciPy can directly solve Roothaan equation, i.e. $\mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C} \boldsymbol{\varepsilon}$. See `scipy.linalg.eigh` ([SciPy API](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigh.html)).

Although this function is very similar to `numpy.linalg.eigh`, `scipy.linalg.eigh` can be more powerful.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

We may use the Fock matrix fock generated in Step 5 to test whether our function works correctly.

Step 7: Build Density Matrix¶

Build the (restricted) density matrix using the occupied molecular orbital coefficients (Szabo and Ostlund, eq 3.145):

$$ D_{\mu \nu} = 2 \sum_i C_{\mu i} C_{\nu i} $$

{admonition}
:class: dropdown

<img src="https://raw.githubusercontent.com/ajz34/PyCrawfordProgProj/master/source/Project_03/assets/Duv.png" style="max-width:300px" />

Notice that not all orbital coefficient are taken into account. Only occupied molecular orbital coefficients is included.

In the figure above, matrix in the middle is the usual representation of orbital coefficients $C_{\mu p}$, where row denotes atomic orbital subscript $\mu$, and column denotes molecular orbital subscript $\nu$. For water molecule, there are 10 electrons, so number of occupied orbitals is 10/2 = 5. So only the first 5 columns (area painted in blue) is considered as $C_{\mu i}$.

We don't have a function determining number of occupied orbitals so far. So the reader is encouraged to implement a simple function to do this job.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

We may use the Fock matrix coeff generated in Step 6 to test whether our function works correctly.

Step 8: Update Density Matrix¶

Review procedure from Step 5 to Step 7:

• Step 5: $D_{\mu \nu} \rightarrow F_{\mu \nu}$
• Step 6: $F_{\mu \nu} \rightarrow C_{\mu p}$
• Step 7: $C_{\mu p} \rightarrow D_{\mu \nu}$

Run these processes in series, and the density updated. If we just simply loop these three steps infinitely, then we except the density converged to stable, i.e. self-consistent.

Write a function that inputs a guess density matrix, output an updated density matrix.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

We use the random-generated density matrix dm_random in Step 5 to test whether our function works correctly.

Step 9: Compute the Total Energy¶

Compute the total energy as (Szabo and Ostlund, eq 3.184 - 3.185)

$$ E_\mathrm{tot} = \frac{1}{2} \sum_{\mu \nu} D_{\mu \nu} (H_{\mu \nu}^\mathrm{core} + F_{\mu \nu}) + E_\mathrm{nuc} $$

where $F_{\mu \nu}$ is obtained from the the same density matrix $D_{\mu \nu}$ in the fomular above. Indeed, $E_\mathrm{tot}$ could be regarded as a function of density matrix $D_{\mu \nu}$. $E_\mathrm{nuc}$ has been already obtained in Step 1.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

We use the random-generated density matrix dm_random in Step 5 to test whether our function works correctly.

{warning}

You may find the converged energy is about -75 Hartree (a.u.). Maybe to you, it is very suspicious that some density could yield much lower energy (-127) than the variational minimum (-75). This is because the density `dm_random` does not satisfy boundary condition $\mathrm{tr} (\mathbf{DS}) = n_\mathrm{elec}$.

Step 10: SCF Loop and Convergence¶

Given initial (empty) density $D_{\mu \nu}^0 = 0$ and initial energy value $E_\mathrm{tot}^0 = 0$. Use the function described in Step 9 to update density matrix until convergence.

What defines convergence is energy and density criterion:

$$ \begin{align} \Delta_E &= |E_\mathrm{tot}^t - E_\mathrm{tot}^{t-1}| < \delta_1 \\ \mathrm{rms}_\mathbf{D} &= \left[ \sum_{\mu \nu} (D_{\mu \nu}^t - D_{\mu \nu}^{t-1})^2 \right]^{-1/2} < \delta_2 \end{align} $$

where superscript $t$ means self-consistent field loop iteration. $\mathrm{rms}_\mathbf{D}$ is actually Frobenius norm of the difference in consecutive density matrix.

In our implementation, we set criterion

$$ \delta_1 = 10^{-10} \, \mathsf{a.u.}, \quad \delta_2 = 10^{-8} \, \mathsf{a.u.} $$

and maximum iteration number set to 64.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

Addition 1: Dipole Moment¶

As discussed in detail in Ch. 3 of the text by Szabo and Ostlund, the calculation of one-electron properties requires density matrix and the relevant property integrals. The electronic contribution to the electric-dipole moment may be computed using (Szabo and Ostlund, eq 3.190 - 3.192)

$$ \boldsymbol{\mu} = - \sum_{\mu \nu} P_{\mu \nu} (\mu | \boldsymbol{r} | \nu) + \sum_A Z_A \boldsymbol{R}_A $$

{admonition}
:class: dropdown

Integral $(\mu | \boldsymbol{r} | \nu)$ is stored in $\boldsymbol{r}$'s three components $(x, y, z)$. $x$ component is stored in file {download}`input/h2o/STO-3G/mux.dat`. For PySCF users, simply feed `"int1e_r"` to `pyscf.gto.Mole.intor` and integral $(\mu | \boldsymbol{r} | \nu)$ is returned in dimension $(3, n_\mathrm{AO}, n_\mathrm{AO})$.

Implementation¶

Reader should fill all NotImplementedError in the following code:

Solution¶

Addition 2: Mulliken Population Analysis¶

A Mulliken population analysis (also described in Szabo and Ostlund, Ch. 3) requires the overlap integrals and the electron density, in addition to information about the number of basis functions centered on each atom. The charge on atom A may be computed as (Szabo and Ostlund, eq 3.196)

$$ q_A = Z_A - \sum_{\mu \in A} \sum_\nu P_{\mu \nu} S_{\nu \mu} $$

where the summation is limited to only those basis functions centered on atom $A$.

{admonition}
:class: dropdown

Due to limited file input, this task could only be accomplished by the help from PySCF. Function `pyscf.gto.Mole.aoslice_by_atom` could help group atomic orbitals by its atom center.

Implementation¶

This task is a little difficult, due to $\mu \in A$ is hard to figure out by program. So this task is not required to be accomplished by one's own. However, hard-coding is possible for a specified system, for example water/STO-3G.

Reader may try to fill NotImplementedError in the following code:

Solution¶

Test Cases¶

Water/STO-3G (Directory)

Water/DZ (Directory)

Water/DZP (Directory) For this molecule, read-input approach is prefered. For explanation, see comment in code.

Methane/STO-3G (Directory)

References¶

• Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory Dover Publication Inc., 1996.

  ISBN-13: 978-0486691862