Hamiltonians

This page explains how ffsim represents Hamiltonians.

Data representation

ffsim includes several classes for representing Hamiltonians in different forms. To take a concrete example, we’ll use the molecular Hamiltonian of the form

\[H = \sum_{\sigma, pq} h_{pq} a^\dagger_{\sigma, p} a_{\sigma, q} + \frac12 \sum_{\sigma \tau, pqrs} h_{pqrs} a^\dagger_{\sigma, p} a^\dagger_{\tau, r} a_{\tau, s} a_{\sigma, q} + \text{constant}.\]

To represent this Hamiltonian, the information that needs to be stored consists of (\(N\) is the number of spatial orbitals):

• The one-body tensor \(h_{pq}\), which is an \(N \times N\) matrix (stored as a NumPy array).

• The two-body tensor \(h_{pqrs}\), which is an \(N \times N \times N \times N\) tensor (stored as a NumPy array).

• The constant, which is a real number (stored as a float).

In order to have some objects to work with, we sample some random instances of these data in the following code cell.

import numpy as np

import ffsim

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

rng = np.random.default_rng(1234)
one_body_tensor = ffsim.random.random_real_symmetric_matrix(norb, seed=rng)
# Pass dtype=float to obtain a real-valued two-body tensor.
# Complex tensors are not fully supported yet.
two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng, dtype=float)
constant = rng.standard_normal()

The molecular Hamiltonian is represented in ffsim using the MolecularHamiltonian class. You initialize the class by passing the three pieces of information (the constant is optional and defaults to zero if not specified):

mol_hamiltonian = ffsim.MolecularHamiltonian(
    one_body_tensor=one_body_tensor,
    two_body_tensor=two_body_tensor,
    constant=constant,
)

The arguments passed to the initialization are now available as attributes of the same name, i.e., mol_hamiltonian.one_body_tensor accesses the one-body tensor.

Operator action

The basic operation that a Hamiltonian represention should support is applying its action, as a linear operator, to a vector. This basic operation can be used to implement more complex ones, such as operator exponentiation and eigenvalue computation. ffsim uses Scipy’s LinearOperator class to support these operations for Hamiltonians. To obtain a LinearOperator, call the linear_operator function:

# Use 2 alpha electrons and 2 beta electrons, as an example.
nelec = (2, 2)

linop = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec)

The linear_operator function requires the number of orbitals and the number of alpha and beta electrons to be passed because this information is needed to fully specify the vector space on which the linear operator acts.

LinearOperators support matrix multiplication with a vector:

vec = ffsim.hartree_fock_state(norb, nelec)

result = linop @ vec

They also work with many of Scipy’s sparse linear algebra routines. For example, you can compute the ground state energy using the eigsh eigenvalue routine:

import scipy.sparse.linalg

eigs, vecs = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
energy = eigs[0]

energy

np.float64(-99.5571707255156)

Time evolution by the Hamiltonian can be computed using expm_multiply:

time = 1.0
evolved_vec = scipy.sparse.linalg.expm_multiply(-1j * time * linop, vec)

/tmp/ipykernel_7680/2190712273.py:2: UserWarning: Trace of LinearOperator not available, it will be estimated. Provide `traceA` to ensure performance.
  evolved_vec = scipy.sparse.linalg.expm_multiply(-1j * time * linop, vec)

When passing a LinearOperator to expm_multiply, Scipy issues a warning if an estimate of the trace is not provided via the traceA argument. You can avoid this warning by passing an estimate of the trace. For Hamiltonians with real-valued tensors, the trace function can compute the exact trace (complex-valued tensors are not supported yet):

trace = ffsim.trace(mol_hamiltonian, norb=norb, nelec=nelec)
evolved_vec_2 = scipy.sparse.linalg.expm_multiply(
    -1j * time * linop, vec, traceA=-1j * time * trace
)

Hamiltonian classes

ffsim includes several classes for representing Hamiltonians of different forms.

MolecularHamiltonian

The MolecularHamiltonian class stores the Hamiltonian in the standard second-quantized form with a one-body tensor, two-body tensor, and constant:

\[H = \sum_{\sigma, pq} h_{pq} a^\dagger_{\sigma, p} a_{\sigma, q} + \frac12 \sum_{\sigma \tau, pqrs} h_{pqrs} a^\dagger_{\sigma, p} a^\dagger_{\tau, r} a_{\tau, s} a_{\sigma, q} + \text{constant}.\]

The data stored by this class consists of:

• one_body_tensor: The \(N \times N\) matrix \(h_{pq}\).

• two_body_tensor: The \(N \times N \times N \times N\) tensor \(h_{pqrs}\).

• constant: A scalar constant energy offset.

DoubleFactorizedHamiltonian

The DoubleFactorizedHamiltonian class stores the Hamiltonian in the double-factorized representation, which decomposes the two-body tensor into a sum of diagonal Coulomb operators acting in rotated orbital bases:

\[H = \sum_{\sigma, pq} \kappa_{pq} a^\dagger_{\sigma, p} a_{\sigma, q} + \frac12 \sum_t \sum_{\sigma\tau, ij} J^{(t)}_{ij} n^{(t)}_{\sigma, i} n^{(t)}_{\tau, j} + \text{constant}'\]

where \(n^{(t)}_{\sigma, i} = \sum_{pq} U^{(t)}_{pi} a^\dagger_{\sigma, p} a_{\sigma, q} U^{(t)*}_{qi}\) is the number operator in the rotated basis defined by the unitary matrix \(U^{(t)}\).

The data stored by this class consists of:

• one_body_tensor: The \(N \times N\) modified one-body tensor \(\kappa_{pq}\).

• diag_coulomb_mats: An array of shape \((L, N, N)\) containing the \(L\) diagonal Coulomb matrices \(J^{(t)}\).

• orbital_rotations: An array of shape \((L, N, N)\) containing the \(L\) unitary matrices \(U^{(t)}\).

• constant: A scalar constant energy offset.

• z_representation: Whether the Hamiltonian uses the “Z” representation (an alternative form that can yield simpler quantum circuits).

See Double-factorized representation of the molecular Hamiltonian for more information.

DiagonalCoulombHamiltonian

The DiagonalCoulombHamiltonian class stores a Hamiltonian whose two-body interaction is diagonal in the computational basis:

\[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_{\sigma, p}^\dagger a_{\sigma, p}\) is the number operator.

The data stored by this class consists of:

• one_body_tensor: The \(N \times N\) matrix \(h_{pq}\).

• diag_coulomb_mats: An array of shape \((2, N, N)\) containing the diagonal Coulomb matrices \(J^{\alpha\alpha}\) and \(J^{\alpha\beta}\). The constraint \(J^{\alpha\alpha} = J^{\beta\beta}\) and \(J^{\alpha\beta} = J^{\beta\alpha}\) is assumed, so only two matrices are needed.

• constant: A scalar constant energy offset.

See Diagonal Coulomb Hamiltonians for more information.