Tutorial: Molkit

Tutorial: Molkit#

Introduction#

This tutorial will be focus in tools available in the new ModelHamiltonian module to manipulate molecular hamiltonians.

1. Instantiating a Molecular Hamiltonian#

import numpy as np
from moha import HamHeisenberg
from moha.molkit.hamiltonians import MolHam
from moha.molkit.utils.tools import to_geminal

# load a hamiltonian from the fcidump file
hami = MolHam()
hami.from_fcidump('example.fcidump')
print(hami)

<moha.molkit.hamiltonians.MolHam object at 0x7f7a91a24130>

/Users/valeriichuiko/Desktop/projects/ModelHamiltonian/moha/molkit/hamiltonians.py:176: UserWarning: Duplicate entries in the FCIDUMP file are ignored
  data = load_fcidump(open(path))

2. Converting the Hamiltonian to the Spin-Orbital Basis#

You can easily perform this conversion by using the spinize_H method, which returns the one-body and two-body integrals in the spin-orbital basis.

Conversion Rules#

  • One-body term:

    • The spatial one-body integral is mapped to the spin-orbital basis as:

      • \(h_{pq}^{\alpha\alpha} = h_{pq}^{\beta\beta} = h_{pq}^{\text{spatial}}\)

      • \(h_{pq}^{\alpha\beta} = h_{pq}^{\beta\alpha} = 0\)

  • Two-body term:

    • The spatial two-body integral is mapped as:

      • \(V_{pqrs}^{\alpha\alpha\alpha\alpha} = V_{pqrs}^{\alpha\beta\alpha\beta} = V_{pqrs}^{\beta\alpha\beta\alpha} = V_{pqrs}^{\text{spatial}}\)

one_body_spin, two_body_spin = hami.spinize_H()
print("Shape of one-body integrals (spin basis):\n", one_body_spin.shape)
print("Shape of two-body integrals (spin basis):\n", two_body_spin.shape)

Shape of one-body integrals (spin basis):
 (8, 8)
Shape of two-body integrals (spin basis):
 (8, 8, 8, 8)

3. Antisymmetrizing the Two-Electron Integrals#

You can apply proper antisymmetrization to the two-electron integrals using the antisymmetrize method.

This ensures the integrals obey the required permutation symmetry for fermionsm following the rule: \(V_{pqrs} = -V_{pqsr} = -V_{qprs} = V_{qpsr}\)

hami.antisymmetrize()
hami.two_body_spin[0, 1, 0, 1]
two_body_aaaa = hami.get_spin_blocks()['aaaa']
two_body_abab = hami.get_spin_blocks()['abab']

print("Elements 0, 1, 0, 1 and 1, 0, 0, 1 of aaaa part:\n", two_body_aaaa[0, 1, 0, 1], two_body_aaaa[1, 0, 0, 1])
print("--"*10)
print("Elements 0, 1, 0, 1 and 1, 0, 0, 1 of abab part:\n", two_body_abab[0, 1, 0, 1], two_body_abab[1, 0, 0, 1])
print("--"*10)
print("Elements 0, 0, 1, 1 for aaaa and abab parts:\n", two_body_aaaa[0, 0, 1, 1], two_body_abab[0, 0, 1, 1])

Elements 0, 1, 0, 1 and 1, 0, 0, 1 of aaaa part:
 0.1544687847893533 -0.1544687847893533
--------------------
Elements 0, 1, 0, 1 and 1, 0, 0, 1 of abab part:
 0.3187342023639616 0.1642654175746083
--------------------
Elements 0, 0, 1, 1 for aaaa and abab parts:
 0.0 0.1642654175746083

4. Converting the Two-Body Term to the Geminal Basis#

The to_geminal function in the utils module converts the two-body term from the spin-orbital basis to the geminal basis.
You can specify the type of two-body term: 'h2' (Hamiltonian) ; we also added conversion to the geminal basis for reduced density matrices: 'rdm2'.

The Hamiltonian in the geminal basis is obtained by the following formula: \(V_{A B} =\frac{1}{2}\left(V_{p q r s}-V_{q p r s}-V_{p q r s}+V_{qprs}\right)\)

geminal = to_geminal(two_body_spin, type='h2')

5. Constructing the Reduced 4-Index Hamiltonian Tensor#

The to_reduced method returns the reduced 4-index Hamiltonian tensor \(k_{pqrs}\) in the spin-orbital basis, given the number of electrons \(N\).

Formula#

\[ k_{pqrs} = \frac{1}{2(N-1)} \left( h_{pq} \delta_{rs} + h_{rs} \delta_{pq} \right) + \frac{1}{2} V_{pqrs} \]

where:

  • \(h_{pq}\): one-electron integrals (spin-orbital basis)

  • \(V_{pqrs}\): two-electron integrals (spin-orbital basis)

  • \(\delta_{pq}\): Kronecker delta

reduced = hami.to_reduced(n_elec=4)
print(reduced.shape)

(8, 8, 8, 8)

6. Extracting Spin Blocks from the Two-Body Spin-Orbital Tensor#

After converting the Hamiltonian to the spin-orbital basis, you can extract the main spin blocks from the two-body tensor using the get_spin_blocks method, by accesing the dictionaries keys: ‘aaaa’, ‘bbbb’, ‘abab’

spin_blocks = hami.get_spin_blocks()
print(spin_blocks.keys())

dict_keys(['aaaa', 'bbbb', 'abab'])

7. Computing Energy of a system#

In this section, we will compute the ground-state energy of a toy molecular system using the Full Configuration Interaction (FCI) method using the pyscf library.

  1. Perform an FCI calculation to obtain the wavefunction and reduced density matrices (RDMs).

  2. Compute the energy using both spin-orbital and geminal representations to verify consistency.

from pyscf import fci
from moha.molkit.utils.tools import to_geminal

# One-electron integrals (spatial orbital basis)
h1 = hami.one_body
h2 = hami.two_body
core_energy = hami.zero_body

N_ELEC, norb  = 4, hami.n_spatial

# Perform Full CI calculation
e_fci, ci_vec = fci.direct_spin1.kernel(h1, np.einsum('pqrs -> prqs', h2), norb, N_ELEC, ecore=core_energy, tol=1e-12)
rdm1s, rdm2s = fci.direct_spin1.make_rdm12s(ci_vec, norb, N_ELEC)

# construct 2-rdm
rdm2s_new = np.zeros((8, 8, 8, 8))
for i in range(4):
    for j in range(4):
        for k in range(4):
            for l in range(4):
                rdm2s_new[i, j, k, l] = rdm2s[0][i, k, j, l] # aaaa
                rdm2s_new[i+4, j+4, k+4, l+4] = rdm2s[2][i, k, j, l] #bbbb
                rdm2s_new[i, j+4, k, l+4] = rdm2s[1][i, k, j, l] #abab
                rdm2s_new[i+4, j, k+4, l] = rdm2s[1][i, k, j, l] #baba



# Convert to reduced form
H = hami.to_reduced(n_elec=N_ELEC)  # H = 0.5 * (h \outer h)/(N - 1) + 0.5 * V
# Compute the energy in the MO basis
E_spin = np.sum(H * rdm2s_new)

# Convert to geminal form
H_gem = to_geminal(H, type='h2')
rdm_gem = to_geminal(rdm2s_new, type='rdm2')

# Energy from geminal formulation
E_gem = np.einsum('pq,pq', H_gem, rdm_gem)

# Print the energies from all formulations
print("E_spin-orbital :", E_spin)
print("E_geminal      :", E_gem)
print("FCI energy     :", e_fci - core_energy)
assert np.allclose(E_spin, E_gem), "Energies do not match!"
assert np.allclose(E_spin, e_fci - core_energy), "Energies do not match!"

E_spin-orbital : -3.0622781536522004
E_geminal      : -3.0622781536522012
FCI energy     : -3.062278153652201
Contents