Hamiltonians

This page explains how ffsim represents Hamiltonians.

Data representation

ffsim includes several classes for representing Hamiltonians in different forms. In this page we’ll focus on 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 concrete 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 via SciPy LinearOperators

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.55717072551563)

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_4257/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
)