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.9985212764983

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.9999923953867433