Chemistry simulations with Majorana propagation

This notebook demonstrates how to simulate a quantum circuit (expressed in Majoranas) using the monoprop package.

Background

The Majorana propagation can efficiently simulate fermionic systems without the need for mapping to qubits, allowing larger system sizes and higher accuracy through controlled truncation of the operator algebra.

Two variants are explored below:

Schrödinger picture – the state is propagated forward.
Heisenberg picture – the Hamiltonian is propagated backward (time-reversed).

1. Imports and Environment

The cell below imports all required packages.

from __future__ import annotations

import numpy as np
from pyscf.tools import fcidump
from pyscf import ao2mo

from monoprop import Circuit, ExpGate, MajoranaPropagator
from monoprop.fermi import MajoranaOperator
from monoprop.integral_conversion import integrals_to_fermion

2. The chemistry setup

We will use the chromium dimer as a toy model, a famous example of a strongly correlated molecule. Here we will use pyscf to obtain the integrals, though in principle any tool can be used. We choose an active space of 12 electrons in 12 orbitals, representing all 3d and 4s orbitals of the chromium atoms. The integrals used were generated with the following code

from pyscf import gto, mcscf, scf

# Define the molecule (Cr2 at 1.65 Å)
mol = gto.M(
    atom="""
    Cr 0.0 0.0 0.0
    Cr 0.0 0.0 1.65
    """,
    basis="def2-svp",
    spin=0,  # singlet
)

# Hartree-Fock
mf = scf.RHF(mol)
mf.kernel()

# Choose active orbitals from canonical HF orbitals (indices are 0-based)
cas_list = [18, 19, 20, 21, 22, 23, 24, 25, 26, 31, 32, 33]

# CASCI to use as reference expectation value (impractical for active spaces > 16 orbitals)
mc = mcscf.CASCI(mf, 12, 12)
mo_ord = mc.sort_mo(cas_list, base=0)
mc.mo_coeff = mo_ord

ref_energy, *_ = mc.kernel()
fcidump.from_mcscf(mc, "Cr2_CASCI.FCIDUMP")

Note: since pyscf may generate different integrals depending on the environment on available hardware, for reproducibility reasons we are going to use precomputed FCIDUMP file for the rest of the notebook.

# Load the precomputed molecule
fcidump_path = "Cr2_CASCI.FCIDUMP"
loaded = fcidump.read(fcidump_path)

ref_energy = -2086.12157623417

# Retrieve integrals
h1e = np.asarray(loaded["H1"])
g2e = np.asarray(loaded["H2"])
ecore = loaded.get("ECORE", 0.0)

norb = loaded.get("NORB", h1e.shape[0])
g2e = ao2mo.restore(1, g2e, norb)

# Convert to unrestricted format
h1 = np.array([h1e, h1e])
g2 = np.array([g2e, g2e, 0.5 * (g2e + g2e.transpose(2, 3, 0, 1))])

# Transform to Fermionic Operator
hamiltonian = integrals_to_fermion(ecore, h1, g2)

n_modes = 2 * norb  # number of spin-orbitals
# occupation in the Hartree-Fock state (the first 6 alpha and 6 beta orbitals)
hf = [0, 1, 2, 3, 4, 5, 12, 13, 14, 15, 16, 17]

Parsing Cr2_CASCI.FCIDUMP

3. Monoprop simulations

Helper Functions

plot_energies – plots $|E_k - E_\text{FCI}|$ on a log scale against circuit gate.

def plot_energies(energies: list[float], ref_energy: float | None = None) -> None:
    import matplotlib.pyplot as plt

    if ref_energy is None or ref_energy == 0:
        ref_energy = min(energies)
    diff = np.abs(np.array(energies) - ref_energy)
    plt.plot(diff, marker="o")
    plt.yscale("log")
    plt.xlabel("Gate")
    plt.ylabel("Energy - Ground State Energy")
    plt.title("Energy convergence in Monoprop")
    plt.grid()
    plt.show()

The Majorana circuit

Here we will use a 10-gate circuit, i.e. a circuit containing 10 Majorana operators $\hat{M}_j$ (that are related to excitation operators) to add correlation. The expectation value or any other properties can then be evaluated as

$E = \langle \Psi_\mathrm{HF} | \hat{U}^\dagger \hat{H} \hat{U} | \Psi_\mathrm{HF} \rangle$

with $\hat{U} = \prod_{j=1}^{10} e^{-i\theta_j \hat{M}_j}$

In monoprop, this can be done by either applying the Majorana operators on the state (Schrödinger picture) or on the operator (Heisenberg picture). The choice of Majorana operator $M_j$ and angle $\theta_j$ can be done iteratively through an ADAPT procedure. Here we just used a pre-optimized circuit.

gates = [
    (7, 8, 38, 40),
    (5, 22, 28, 46),
    (1, 2, 18, 20),
    (25, 26, 42, 44),
    (3, 20, 26, 44),
    (1, 18, 24, 42),
    (11, 12, 34, 36),
    (5, 22, 26, 44),
    (3, 20, 28, 46),
    (1, 18, 28, 46),
]
angles = [
    -1.57079633,
    -0.20831642,
    -0.08161787,
    -0.16063185,
    -0.11953465,
    -0.1301597,
    -0.15917059,
    -0.12542137,
    -0.12696826,
    0.12298328,
]

# Pure generator gates (no angle); the angles are supplied separately to the
# Circuit constructor, gate i pairing with angle i.
circuit = [ExpGate(MajoranaOperator({tuple(gate): 1.0}, n_modes)) for gate in gates]

Schrödinger Picture

In the Schrödinger picture the ansatz acts on the initial state. The simulator propagates the state forward through the ansatz layers.

Key parameters for this run:

Parameter	Value	Meaning
`cutoff`	8	Maximum Majorana-string length kept in the Hamiltonian
`schrodinger_cutoff`	`cutoff + 2`	Maximum string length in the propagated state — set slightly larger than `cutoff` to capture leading corrections

The initial expectation value (before any gates added) gives the Hartree-Fock expectation value within the truncation.

# First, we set the cutoffs for the simulation. The `cutoff` parameter determines
# the maximum length of the Majorana strings that will be included in the simulation, while the `schrodinger_cutoff`
# parameter determines the maximum length of the Majorana strings that will be included in the
# initial state.
cutoff = 8
schrodinger_cutoff = cutoff + 2

quantum_circuit = Circuit(gates=circuit, parameters=angles, initial_state=hf)

simulator = MajoranaPropagator(
    hamiltonian,
    quantum_circuit.initial_state,
    cutoff=cutoff,
    schrodinger_cutoff=schrodinger_cutoff,
)

initial_expval = simulator.expval()
print(f"Initial energy: {initial_expval:.4f} (diff: {initial_expval - ref_energy:.4f})")

Initial energy: -2085.5984 (diff: 0.5232)

We can now evolve the state using the gates one by one to see the successive improvement in the expectation value.

expvals_schrodinger = [initial_expval]
for i in range(10):
    # Apply one Majorana gate at a time to the state (immediate contraction). Each single-gate
    # Circuit pairs gate i with its own angle value.
    gate_circuit = Circuit(
        gates=[quantum_circuit.gates[i]],
        parameters=[quantum_circuit.parameters[i]],
    )
    simulator.propagate(gate_circuit)
    simulator_expval = simulator.expval()
    print(
        f"Energy after gate {i + 1:1d}: {simulator_expval:.4f} (diff: {simulator_expval - ref_energy:.4f})"
    )
    expvals_schrodinger.append(simulator_expval)
plot_energies(expvals_schrodinger, ref_energy)

Energy after gate 1: -2085.8354 (diff: 0.2862)
Energy after gate 2: -2085.8800 (diff: 0.2416)
Energy after gate 3: -2085.8937 (diff: 0.2279)
Energy after gate 4: -2085.9139 (diff: 0.2076)
Energy after gate 5: -2085.9266 (diff: 0.1950)
Energy after gate 6: -2085.9406 (diff: 0.1810)
Energy after gate 7: -2085.9464 (diff: 0.1752)
Energy after gate 8: -2085.9546 (diff: 0.1670)
Energy after gate 9: -2085.9633 (diff: 0.1582)
Energy after gate 10: -2085.9730 (diff: 0.1485)

Heisenberg Picture

In the Heisenberg picture, the roles of state and operator are swapped: instead of propagating the state forward, we propagate the Hamiltonian backward through the ansatz:

$\hat{H}(\boldsymbol{\theta}) = e^{i\theta_1 \hat{M}_1 } \cdots e^{i\theta_k \hat{M}_k}\, \hat{H}\, e^{-i\theta_k \hat{M}_k} \cdots e^{-i\theta_1 \hat{M}_1}$

The expectation value is then evaluated as $E = \langle \Psi_\text{HF} | \hat{H}(\boldsymbol{\theta}) | \Psi_\text{HF} \rangle$ .

Advantages over the Schrödinger picture:

The operator algebra stays within the truncated Majorana space without a growing state vector.
A larger cutoff can typically be used, capturing more correlations.
Memory scaling is generally more favourable for deeper ansätze.

# Heisenberg picture: schrodinger_cutoff=None tells the simulator to propagate
# the Hamiltonian operator rather than the state. A larger cutoff is viable here
# because the operator space does not grow with ansatz depth.
cutoff = 10
schrodinger_cutoff = None

simulator = MajoranaPropagator(
    hamiltonian, quantum_circuit.initial_state, cutoff=cutoff
)

initial_expval = simulator.expectation_value()
print(f"Initial energy: {initial_expval:.4f} (diff: {initial_expval - ref_energy:.4f})")

energies_heisenberg = []
for i in range(11):
    simulator = MajoranaPropagator(
        hamiltonian, quantum_circuit.initial_state, cutoff=cutoff
    )
    # Propagate the length-i prefix of the circuit (gates and their angle values).
    prefix = Circuit(
        gates=quantum_circuit.gates[:i],
        parameters=quantum_circuit.parameters[:i],
    )
    simulator.propagate(prefix)
    simulator_expval = simulator.expectation_value()
    print(
        f"Energy after gate {i:1d}: {simulator_expval:.4f} (diff: {simulator_expval - ref_energy:.4f})"
    )
    energies_heisenberg.append(simulator_expval)
plot_energies(energies_heisenberg, ref_energy)

Initial energy: -2085.5984 (diff: 0.5232)


Energy after gate 0: -2085.5984 (diff: 0.5232)


Energy after gate 1: -2085.8354 (diff: 0.2862)


Energy after gate 2: -2085.8800 (diff: 0.2416)


Energy after gate 3: -2085.8937 (diff: 0.2279)


Energy after gate 4: -2085.9139 (diff: 0.2076)


Energy after gate 5: -2085.9266 (diff: 0.1950)


Energy after gate 6: -2085.9406 (diff: 0.1810)


Energy after gate 7: -2085.9464 (diff: 0.1752)


Energy after gate 8: -2085.9546 (diff: 0.1670)


Energy after gate 9: -2085.9634 (diff: 0.1582)


Energy after gate 10: -2085.9729 (diff: 0.1486)

Comparison: Schrödinger vs Heisenberg

Both runs used the same circuits, but with different cutoffs and propagation strategies. The plot below overlays their convergence curves on a single log-scale axis so the trade-off is immediately visible.

Note that the two runs are not apples-to-apples: the Heisenberg picture uses a larger cutoff (10 vs 8), which generally yields a better starting expectation value and faster convergence per iteration at the cost of more compute per expectation value evaluation.

import matplotlib.pyplot as plt

diff_s = np.abs(np.array(expvals_schrodinger) - ref_energy)
diff_h = np.abs(np.array(energies_heisenberg) - ref_energy)

fig, ax = plt.subplots()
ax.plot(diff_s, marker="o", label=f"Schrödinger (cutoff={8})")
ax.plot(diff_h, marker="s", label=f"Heisenberg (cutoff={10})")
ax.set_yscale("log")
ax.set_xlabel("Gates")
ax.set_ylabel(r"$|E - E_\mathrm{FCI}|$")
ax.set_title("Energy convergence in monoprop")
ax.legend()
ax.grid(True, which="both", ls="--", alpha=0.5)
plt.tight_layout()
plt.show()

4. Summary and Next Steps

This notebook demonstrated Majorana simulations on a 12-electron/12-orbital molecular active space.

Things to try:

Vary cutoff and schrodinger_cutoff to study the accuracy/cost trade-off.
Try different properties ( $S^2$ , natural orbitals)

Chemistry simulations with Majorana propagation

On this page