Compute expectation values

Computing the expectation value of an operator on a state vector involves the following steps:

1. Convert your operator to a SciPy LinearOperator using the ffsim.linear_operator protocol.

2. Use matrix multiplication and np.vdot to compute the expectation value.

If you already have a SciPy LinearOperator, you can skip Step 1. Otherwise, you can get a LinearOperator from any object that implements ffsim’s SupportsLinearOperator protocol. Such objects include the Hamiltonian classes as well as the FermionOperator class.

The following code cell demonstrates these steps for a molecular Hamiltonian with randomly generated one- and two-body tensors.

import numpy as np

import ffsim

# Initialize pseudorandom number generator
rng = np.random.default_rng(12345)

# Set numbers of orbitals and electrons to use for example
norb = 8
nelec = (5, 5)

# Generate a random state vector
vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)

# Generate random one- and two-body tensors
one_body_tensor = ffsim.random.random_hermitian(norb, seed=rng)
two_body_tensor = ffsim.random.random_two_body_tensor(norb, seed=rng)

# Form MolecularHamiltonian object
mol_ham = ffsim.MolecularHamiltonian(
    one_body_tensor=one_body_tensor, two_body_tensor=two_body_tensor
)

# Convert the Hamiltonian to a SciPy LinearOperator
linop = ffsim.linear_operator(mol_ham, norb=norb, nelec=nelec)

# The expectation value of a Hamiltonian is real, so we can discard the imaginary part
expectation_value = np.vdot(vec, linop @ vec).real

print(expectation_value)

-998.3155443224805

For a molecular system, the expectation value of the Hamiltonian on the Hartree-Fock state equals the Hartree-Fock energy.

import pyscf
import pyscf.data.elements

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]], basis="sto-6g", symmetry="Dooh")
scf = pyscf.scf.RHF(mol).run()

# Use the MolecularData class to obtain the MolecularHamiltonian for the molecule
mol_data = ffsim.MolecularData.from_scf(scf)
norb, nelec = mol_data.norb, mol_data.nelec
mol_hamiltonian = mol_data.hamiltonian

# Compute the expectation value on the Hartree-Fock state
linop = ffsim.linear_operator(mol_hamiltonian, norb=norb, nelec=nelec)
vec = ffsim.hartree_fock_state(norb, nelec)
expectation_value = np.vdot(vec, linop @ vec).real

# It should match the Hartree-Fock energy
print(expectation_value)

converged SCF energy = -108.464957764796
-108.46495776479557

As mentioned above, you can obtain LinearOperators from FermionOperators as well as other Hamiltonian classes. The following code example constructs a Fermi-Hubbard model Hamiltonian as a FermionOperator, then converts it to an instance of the DiagonalCoulombHamiltonian class. Both classes can be converted to LinearOperators, yielding the same expectation value.

# Set Fermi-Hubbard model parameters
norb_x = 3
norb_y = 2
tunneling = 1.0
interaction = 4.0
norb = norb_x * norb_y
nelec = (norb // 2, norb // 2)

# Create the Hamiltonian as a FermionOperator
op = ffsim.fermi_hubbard_2d(
    norb_x=norb_x, norb_y=norb_y, tunneling=tunneling, interaction=interaction
)

# Generate a random state vector
vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=rng)

# Compute the expectation value
linop = ffsim.linear_operator(op, norb=norb, nelec=nelec)
print(np.vdot(vec, linop @ vec).real)

# Convert to a DiagonalCoulombHamiltonian and compute the expectation value
dc_ham = ffsim.DiagonalCoulombHamiltonian.from_fermion_operator(op)
linop = ffsim.linear_operator(dc_ham, norb=norb, nelec=nelec)
print(np.vdot(vec, linop @ vec).real)

5.926393653639861
5.926393653639866