Diagonal Coulomb Hamiltonians¶

This page explains the diagonal Coulomb Hamiltonian, a restricted but important class of fermionic Hamiltonians.

Definition¶

A diagonal Coulomb Hamiltonian has the form

\[H = \sum_{\sigma, pq} h_{pq} a^\dagger_{\sigma, p} a_{\sigma, q} + \frac12 \sum_{\sigma \tau, pq} J^{\sigma \tau}_{pq} n_{\sigma, p} n_{\tau, q} + \text{constant}\]

where \(n_{\sigma, p} = a^\dagger_{\sigma, p} a_{\sigma, p}\) is the number operator on orbital \(p\) with spin \(\sigma\).

The Hamiltonian is specified by the following data (\(N\) is the number of spatial orbitals):

The one-body tensor \(h_{pq}\), which is an \(N \times N\) Hermitian matrix.
The diagonal Coulomb matrices \(J^{\sigma \tau}\), given as a pair of \(N \times N\) real symmetric matrices specifying the independent coefficients for alpha-alpha and alpha-beta interactions. We require that \(J^{\alpha\alpha} = J^{\beta\beta}\) and \(J^{\alpha\beta} = J^{\beta\alpha}\), so only two matrices are needed.
The constant, which is a real number.

Compared to the general molecular Hamiltonian, the two-body interaction is restricted to density-density form (\(n_{\sigma, p} n_{\tau, q}\) terms). This restriction means the two-body part is specified by \(O(N^2)\) parameters rather than \(O(N^4)\), which enables more efficient simulation.

Data representation¶

A diagonal Coulomb Hamiltonian is represented in ffsim using the DiagonalCoulombHamiltonian class. You can initialize it by passing the one-body tensor, the diagonal Coulomb matrices, and an optional constant.

[1]:

import numpy as np

import ffsim

# Use 4 spatial orbitals, as an example.
norb = 4

rng = np.random.default_rng(12345)
one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
diag_coulomb_mats = np.array(
    [
        ffsim.random.random_real_symmetric_matrix(norb, seed=rng),
        ffsim.random.random_real_symmetric_matrix(norb, seed=rng),
    ]
)
constant = rng.standard_normal()

hamiltonian = ffsim.DiagonalCoulombHamiltonian(
    one_body_tensor=one_body_tensor,
    diag_coulomb_mats=diag_coulomb_mats,
    constant=constant,
)

print(f"Number of orbitals: {hamiltonian.norb}")
print(f"One-body tensor shape: {hamiltonian.one_body_tensor.shape}")
print(f"Diagonal Coulomb matrices shape: {hamiltonian.diag_coulomb_mats.shape}")
print(f"Constant: {hamiltonian.constant}")

Number of orbitals: 4
One-body tensor shape: (4, 4)
Diagonal Coulomb matrices shape: (2, 4, 4)
Constant: -1.6404178369858733

The diag_coulomb_mats attribute has shape (2, N, N). The first matrix (diag_coulomb_mats[0]) contains the alpha-alpha (equivalently, beta-beta) interactions, and the second matrix (diag_coulomb_mats[1]) contains the alpha-beta (equivalently, beta-alpha) interactions.

Example: Fermi-Hubbard model¶

The Fermi-Hubbard model is a widely studied example of a Hamiltonian with diagonal Coulomb form, since its two-body interaction is an onsite density-density interaction. The two-dimensional Fermi-Hubbard Hamiltonian on a square lattice is given by

\[H = -t \sum_{\sigma, \langle pq \rangle} (a^\dagger_{\sigma, p} a_{\sigma, q} + a^\dagger_{\sigma, q} a_{\sigma, p}) + U \sum_p n_{\alpha, p} n_{\beta, p} - \mu \sum_p (n_{\alpha, p} + n_{\beta, p})\]

where \(t\) is the tunneling amplitude, \(U\) is the onsite interaction strength, and \(\mu\) is the chemical potential. The index \(\langle pq \rangle\) runs over pairs of neighboring orbitals on the lattice.

The two-dimensional Fermi-Hubbard Hamiltonian can be constructed as a FermionOperator using the fermi_hubbard_2d function, and then converted to a DiagonalCoulombHamiltonian using the from_fermion_operator method.

[2]:

norb_x = 2
norb_y = 2
norb = norb_x * norb_y
nelec = (2, 2)

# Build the Fermi-Hubbard Hamiltonian as a FermionOperator
hubbard_op = ffsim.fermi_hubbard_2d(
    norb_x=norb_x,
    norb_y=norb_y,
    tunneling=1.0,
    interaction=4.0,
    chemical_potential=2.0,
)

# Convert to a DiagonalCoulombHamiltonian
hubbard_ham = ffsim.DiagonalCoulombHamiltonian.from_fermion_operator(hubbard_op)

print(f"One-body tensor:\n{hubbard_ham.one_body_tensor.real}")
print(f"\nDiagonal Coulomb matrix (alpha-alpha):\n{hubbard_ham.diag_coulomb_mats[0]}")
print(f"\nDiagonal Coulomb matrix (alpha-beta):\n{hubbard_ham.diag_coulomb_mats[1]}")
print(f"\nConstant: {hubbard_ham.constant}")

One-body tensor:
[[-2. -1. -1.  0.]
 [-1. -2.  0. -1.]
 [-1.  0. -2. -1.]
 [ 0. -1. -1. -2.]]

Diagonal Coulomb matrix (alpha-alpha):
[[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]

Diagonal Coulomb matrix (alpha-beta):
[[4. 0. 0. 0.]
 [0. 4. 0. 0.]
 [0. 0. 4. 0.]
 [0. 0. 0. 4.]]

Constant: 0

We can verify that the DiagonalCoulombHamiltonian representation yields the same ground state energy as the original FermionOperator.

[3]:

import scipy.sparse.linalg

# Compute the ground state energy using the DiagonalCoulombHamiltonian
linop = ffsim.linear_operator(hubbard_ham, norb=norb, nelec=nelec)
eigs, _ = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
energy_diag_coulomb = eigs[0]

# Compute the ground state energy using the FermionOperator
linop_op = ffsim.linear_operator(hubbard_op, norb=norb, nelec=nelec)
eigs_op, _ = scipy.sparse.linalg.eigsh(linop_op, k=1, which="SA")
energy_op = eigs_op[0]

print(f"Energy (DiagonalCoulombHamiltonian): {energy_diag_coulomb}")
print(f"Energy (FermionOperator):            {energy_op}")

np.testing.assert_allclose(energy_diag_coulomb, energy_op)

Energy (DiagonalCoulombHamiltonian): -10.10274848346205
Energy (FermionOperator):            -10.10274848346206

Time evolution via Trotter-Suzuki formulas¶

The diagonal Coulomb Hamiltonian can be decomposed into two terms for the purpose of Trotter-Suzuki simulation:

\[H = H_0 + H_1 + \text{constant},\]

where

\[H_0 = \sum_{\sigma, pq} h_{pq} a^\dagger_{\sigma, p} a_{\sigma, q}\]

is a quadratic Hamiltonian and

\[H_1 = \frac12 \sum_{\sigma \tau, pq} J^{\sigma \tau}_{pq} n_{\sigma, p} n_{\tau, q}\]

is a diagonal Coulomb operator. Note,

\(H_0\) can be simulated exactly using orbital rotations.
\(H_1\) can be simulated exactly using apply_diag_coulomb_evolution.

This split-operator approach is implemented in ffsim by the function simulate_trotter_diag_coulomb_split_op. As with other Trotter functions in ffsim, the order parameter controls the order of the Trotter-Suzuki formula: order=0 is the first-order asymmetric formula (the default), order=1 is the first-order symmetric (second-order) formula, and so on.

In the following example, we compare the Trotter-approximated time evolution against the exact time evolution (computed using expm_multiply).

[4]:

# Initial state
vec = ffsim.hartree_fock_state(norb, nelec)

# Evolution time
time = 1.0

# Compute exact time evolution
linop = ffsim.linear_operator(hubbard_ham, norb=norb, nelec=nelec)
trace = ffsim.trace(hubbard_ham, norb=norb, nelec=nelec)
exact = scipy.sparse.linalg.expm_multiply(
    -1j * time * linop, vec, traceA=-1j * time * trace
)

# Compute Trotter-approximated time evolution with increasing number of steps
for n_steps in [1, 2, 5, 10]:
    result = ffsim.simulate_trotter_diag_coulomb_split_op(
        vec,
        hubbard_ham,
        time,
        norb=norb,
        nelec=nelec,
        n_steps=n_steps,
        order=1,
    )
    fidelity = abs(np.vdot(result, exact))
    print(f"n_steps = {n_steps:2d}, fidelity = {fidelity:.8f}")

n_steps =  1, fidelity = 0.45702529
n_steps =  2, fidelity = 0.95880093
n_steps =  5, fidelity = 0.99915103
n_steps = 10, fidelity = 0.99994861

As expected, the fidelity increases with the number of Trotter steps.