How to simulate Trotterized time evolution for the molecular Hamiltonian

In this guide, we show how to use ffsim to simulate Trotterized time evolution for the molecular Hamiltonian in the double-factorized representation.

First, let’s build a molecule to use as an example.

import pyscf

import ffsim

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

# Define active space
n_frozen = pyscf.data.elements.chemcore(mol)
active_space = range(n_frozen, mol.nao_nr())

# Get molecular data and Hamiltonian
scf = pyscf.scf.RHF(mol).run()
mol_data = ffsim.MolecularData.from_scf(scf, active_space=active_space)
norb, nelec = mol_data.norb, mol_data.nelec
mol_hamiltonian = mol_data.hamiltonian

print(f"norb = {norb}")
print(f"nelec = {nelec}")

converged SCF energy = -108.464957764796
norb = 8
nelec = (5, 5)

Next, let’s get the Hamiltonian in double-factorized form.

df_hamiltonian = ffsim.DoubleFactorizedHamiltonian.from_molecular_hamiltonian(
    mol_hamiltonian
)

print(
    f"Number of terms in double factorization: {len(df_hamiltonian.diag_coulomb_mats)}"
)

Number of terms in double factorization: 35

Now, let’s apply time evolution to an initial state, which we’ll take as the Hartree-Fock state here.

# Construct the initial state.
initial_state = ffsim.hartree_fock_state(norb, nelec)

# Set the evolution time.
time = 1.0

# Perform the Trotterized time evolution
final_state = ffsim.simulate_trotter_double_factorized(
    initial_state,
    df_hamiltonian,
    time,
    norb=norb,
    nelec=nelec,
    n_steps=5,
    order=0,
)

To check the error from Trotterization, we can compute the result of exact time evolution and compare it with the final state we got.

import numpy as np
import scipy.sparse.linalg

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

# Compute the exact result of time evolution
exact_state = scipy.sparse.linalg.expm_multiply(
    -1j * time * linop,
    initial_state,
    traceA=-1j * time * ffsim.trace(mol_hamiltonian, norb=norb, nelec=nelec),
)

# Compute fidelity between results from Trotterized evolution and exact evolution
fidelity = abs(np.vdot(final_state, exact_state))

print(f"Fidelity of Trotter-evolved state with exact state: {fidelity}")

Fidelity of Trotter-evolved state with exact state: 0.9985210982781135

To reduce the Trotter error, we can increase the number of steps or use a higher-order formula.

final_state = ffsim.simulate_trotter_double_factorized(
    initial_state,
    df_hamiltonian,
    time,
    norb=norb,
    nelec=nelec,
    n_steps=5,
    order=1,
)

fidelity = abs(np.vdot(final_state, exact_state))

print(f"Fidelity of Trotter-evolved state with exact state: {fidelity}")

Fidelity of Trotter-evolved state with exact state: 0.9999923953888441