Tutorial: Richardson-Gaudin model

Tutorial: Richardson-Gaudin model#

This tutorial demonstrates how to use the Richardson-Gaudin model with the ModelHamiltonian package. Specifically, it shows how to set up the connectivity matrix and perform the FCI calculation using the pyscf library.

The primary focus of this tutorial, however, it to demostrate the flexibility behind the Model Hamiltonin package. To do so, we will show how to use this library to build a Richardson Model that’s given by:

\[ \hat{H}_{B C S}=\frac{1}{2} \sum_{i=1}^N \varepsilon_i\left(a_{i \uparrow}^{\dagger} a_{1 \uparrow}+a_{i \downarrow}^{\dagger} a_{i \downarrow}\right)-\frac{g}{2} \sum_{i, j=1}^N a_{i \uparrow}^{\dagger} a_{i \downarrow}^{\dagger} a_{j \downarrow} a_{j \uparrow} \]

This definition is given in a paper and can be used to validate numerical results that were kindly provided by Paul Johnson.

1. Match the definition#

In our case, we define the Richardson model as follows:

\[ \hat{H}_{R G}=\sum_p\left(\mu_p^Z-J_{p p}^{\mathrm{eq}}\right) S_p^Z+\sum_{p q} J_{p q}^{\mathrm{eq}} S_p^{+} S_q^{-} \]

To match the above definition we need to move around some of the terms. Skipping the details, we can match the definition by setting the following parameters:

\[\begin{split} \begin{align*} \mu_p &= \varepsilon_p + J_{p p}^{\mathrm{eq}} \\ J_{p q}^{\mathrm{eq}} &= -\frac{g}{2} \\ \end{align*} \end{split}\]

2. Provided data#

For the values of \(\varepsilon_i\) = [0, 1, 2, 3, 4, 5] and ground state energy of 6 electron system, the following data was provided by the author:

g (unitless)

Energy (unitless)

3.536

-14.20687553

460.5454545454545

-2755.77652711

120.8888888888889

-717.84780933

8.182

-41.80555486

6.869

-33.96821191

480.5252525252525

-2875.655157

5.556

-26.14992638

1.011

-0.14612126

6.263

-30.35668599

4.849

-21.95332479

Calculation#

import numpy as np
from moha.hamiltonians import HamRG
from pyscf import fci
from tabulate import tabulate

# Define the connectivity matrix assuming a fully connected system
connectivity = np.ones((6, 6))
g_values = [
    3.536,
    460.5454545454545,
    120.88888888888889,
    8.182,
    6.869000000000001,
    480.52525252525254,
    5.556000000000001,
    1.011,
    6.263000000000001,
    4.849000000000001
]

energy_values = [
    -14.20687553,
    -2755.77652711,
    -717.84780933,
    -41.80555486,
    -33.96821191,
    -2875.655157,
    -26.14992638,
    -0.14612126,
    -30.35668599,
    -21.95332479
]
energy_computed = []
eps = np.arange(6)
# Create the Richardson-Gaudin Hamiltonian
for g in g_values:
    # convert parameters to match the definition
    J = -g / 2
    mu = eps + J   

    # Create the Hamiltonian
    ham = HamRG(J_eq = J*connectivity,
                mu = mu)
    e0 = ham.generate_zero_body_integral()

    h1 = ham.generate_one_body_integral(dense=True, basis='spatial basis')
    h2 = ham.generate_two_body_integral(dense=True, basis='spatial basis', sym=4)

    # convert to chemists notation
    h2_ch = np.einsum('ijkl->ikjl', h2)


    n_sites = 6
    occ = (n_sites//2, n_sites//2)

    # Solve the FCI problem
    e, ci = fci.direct_spin0.kernel(1*h1, 2*h2_ch, nelec=occ, norb=n_sites, nroots=6, ecore=0)
    energy_computed.append(e[0])

# Print the results in a table format

table_data = []
for g, energy, computed in zip(g_values, energy_values, energy_computed):
    error = abs(computed - energy)
    table_data.append([f"{g:.3f}", f"{energy:.6f}", f"{computed:.6f}", f"{error:.6f}"])

headers = ["g", "Energy", "Computed", "Error"]
print(tabulate(table_data, headers=headers, tablefmt="github"))

|       g |       Energy |     Computed |   Error |
|---------|--------------|--------------|---------|
|   3.536 |   -14.2069   |   -14.2069   |       0 |
| 460.545 | -2755.78     | -2755.78     |       0 |
| 120.889 |  -717.848    |  -717.848    |       0 |
|   8.182 |   -41.8056   |   -41.8056   |       0 |
|   6.869 |   -33.9682   |   -33.9682   |       0 |
| 480.525 | -2875.66     | -2875.66     |       0 |
|   5.556 |   -26.1499   |   -26.1499   |       0 |
|   1.011 |    -0.146121 |    -0.146121 |       0 |
|   6.263 |   -30.3567   |   -30.3567   |       0 |
|   4.849 |   -21.9533   |   -21.9533   |       0 |
Contents