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
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.
[1]:
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):
[2]:
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:
[3]:
# 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:
[4]:
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:
[5]:
import scipy.sparse.linalg
eigs, vecs = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
energy = eigs[0]
energy
[5]:
np.float64(-99.55717072551563)
Time evolution by the Hamiltonian can be computed using expm_multiply
:
[6]:
time = 1.0
evolved_vec = scipy.sparse.linalg.expm_multiply(-1j * time * linop, vec)
/tmp/ipykernel_4237/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):
[7]:
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
)