{ "cells": [ { "cell_type": "markdown", "id": "4881c197", "metadata": {}, "source": [ "# Chemistry to qubit Hamiltonians\n", "\n", "This tutorial describes how to load a molecular Hamiltonian from [$\\texttt{pyscf}$](https://pyscf.org/) and convert it to a qubit Hamiltonian object of the $\\texttt{qpe-toolbox}$'s `Hamiltonian` class.\n", "The `Hamiltonian` is designed to be compatible with $\\texttt{quimb}$.\n", "\n", "Here we show how to use $\\texttt{quimb}$ to find the ground state energy with a classical algorithm: the Density-Matrix Renormalization Group (DMRG). In the other tutorials we encode the Hamiltonian spectrum into unitary operators using Trotterization or Block Encoding, and compute the energy with a quantum algorithm: the Quantum Phase Estimation algorithm, using $\\texttt{quimb}$'s quantum circuit simulation mode.\n", "\n", "This notebook assumes the reader is already familiar with basic concepts of DMRG and Matrix Product States (MPS). For a review, we refer to e.g. [Schollwöck, The density-matrix renormalization group in the age of matrix product states](https://arxiv.org/abs/1008.3477)" ] }, { "cell_type": "code", "execution_count": 1, "id": "9a082168", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:46:46.635435Z", "iopub.status.busy": "2026-04-03T03:46:46.635265Z", "iopub.status.idle": "2026-04-03T03:46:51.536924Z", "shell.execute_reply": "2026-04-03T03:46:51.536169Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from IPython.display import display\n", "from pyscf import gto\n", "from quimb.tensor import DMRG2\n", "\n", "from qpe_toolbox.hamiltonian import chemistry_hamiltonian" ] }, { "cell_type": "markdown", "id": "86942b9b", "metadata": {}, "source": [ "The `chemistry_hamiltonian` function is built on the `Hamiltonian` class. It takes a molecule described with a [`Mole`](https://pyscf.org/user/gto.html) object from $\\texttt{pyscf}$, performs the Hartree-Fock calculation and converts the molecular Hamiltonian into a `Hamiltonian` object describing the qubit Hamiltonian." ] }, { "cell_type": "markdown", "id": "183787ce", "metadata": {}, "source": [ "## $H_2$\n", "\n", "Let us start with the $H_2$ molecule in the minimal atomic orbital basis set STO-3G. We give the geometry and basis set as keyword arguments to the function `pyscf.gto.M` to initialize the molecular structure object (spin and charge can also be specified, see below $O_2$).\n", "Then we call the `chemistry_hamiltonian`.\n", "\n", "* The `hf_mode` parameter describes how spin is assigned to the spatial molecular orbitals in the Hartree Fock calculations. Restricted Hartree-Fock (`\"rhf\"`) uses the same molecular orbital twice, one for each state in the spin doublet, while [Unrestricted Hartree-Fock](https://en.wikipedia.org/wiki/Unrestricted_Hartree%E2%80%93Fock) (`\"uhf\"`) uses different molecular orbitals for the different eigenstates of $S_z$. Unrestricted Hartree-Fock is mostly used for open-shell molecules.\n", "* Optionally, the FCI or CCSD energies can be computed." ] }, { "cell_type": "code", "execution_count": 2, "id": "d24f65ad", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:46:51.539032Z", "iopub.status.busy": "2026-04-03T03:46:51.538711Z", "iopub.status.idle": "2026-04-03T03:46:52.359239Z", "shell.execute_reply": "2026-04-03T03:46:52.358229Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged SCF energy = -1.116998996754\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "nOrb : 2\n", "nElec : 2\n", "E_HF : -1.1169989968\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E_CI : -1.1373060358\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E(CCSD) = -1.137306193391969 E_corr = -0.02030719663796497\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E_CCSD: -1.1373061934\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "n_qubits = 4\n" ] } ], "source": [ "mol = gto.M(\n", " atom=[(\"H\", (0.0, 0.0, 0.0)), (\"H\", (0.0, 0.0, 0.735))],\n", " basis=\"sto-3g\",\n", ")\n", "\n", "h2_ham_sto3g = chemistry_hamiltonian(mol, \"rhf\", do_fci=True, do_ccsd=True)\n", "print(f\"n_qubits = {h2_ham_sto3g.n_qubits}\")" ] }, { "cell_type": "markdown", "id": "ca9eea10", "metadata": {}, "source": [ "With this choice of atomic orbital basis set, the molecule is described with $2$ molecular orbitals (a single orbital per atom), giving $4$ qubits when spin is taken into account.\n", "\n", "The FCI energy is given by exact diagonalization of the Hamiltonian. It is the best ground state energy estimation within this choice of basis set.\n", "In chemistry, an energy estimation must reach the chemical accuracy standard, defined to be $1.6 \\text{mHa} \\approx 500 \\text{K}$.\n", "The HF energy is far from the FCI energy by $20 \\text{mHa}$, while the CCSD energy is close within $10^{-6} \\text{Ha}$.\n", "\n", "However, the basis set itself only gives an approximate description of the true solutions to the Schrödinger equation. Indeed, note that with this small basis set, the FCI energy is $-1.13$ Ha, while reported values in the literature ([here](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.83.2541)) for $H_2$ are around $-1.16$ Ha. Even with exact diagonalization, the accuracy is at best $30 \\text{mHa}$ in this basis.\n", "Hence, reaching chemical accuracy will always mean using much larger basis sets.\n", "For a more precise discussion on the importance of atomic orbitals basis sets for the quality of energy estimation, see [E.V. Elfving et al., arXiv:2009.12472](https://arxiv.org/pdf/2009.12472).\n", "\n", "In the following, we consider the cc-pvdz basis set with $5$ orbitals per Hydrogen atom." ] }, { "cell_type": "code", "execution_count": 3, "id": "51e4f5fa", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:46:52.360799Z", "iopub.status.busy": "2026-04-03T03:46:52.360638Z", "iopub.status.idle": "2026-04-03T03:46:53.210415Z", "shell.execute_reply": "2026-04-03T03:46:53.209635Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged SCF energy = -1.12862276992659\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "nOrb : 10\n", "nElec : 2\n", "E_HF : -1.1286227699\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E_CI : -1.1632095576\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E(CCSD) = -1.163209563562649 E_corr = -0.03458679363605879\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E_CCSD: -1.1632095636\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "n_qubits = 20\n" ] } ], "source": [ "mol = gto.M(\n", " atom=[(\"H\", (0.0, 0.0, 0.0)), (\"H\", (0.0, 0.0, 0.735))],\n", " basis=\"cc-pvdz\",\n", ")\n", "\n", "h2_ham = chemistry_hamiltonian(mol, \"rhf\", do_fci=True, do_ccsd=True)\n", "print(f\"n_qubits = {h2_ham.n_qubits}\")" ] }, { "cell_type": "markdown", "id": "65e8aaf8", "metadata": {}, "source": [ "The FCI energy is now $-1.1632 \\text{Ha}$, within chemical accuracy compared to the reported value of $-1.1640 \\text{Ha}$.\n", "\n", "Let us now build the MPO representation of the qubit Hamiltonian $H$ using the `to_mpo` method from the `Hamiltonian` class, and perform a DMRG calculation of the energy using $\\texttt{quimb}$'s `DMRG2` solver." ] }, { "cell_type": "code", "execution_count": 4, "id": "ae8fee30", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:46:53.211865Z", "iopub.status.busy": "2026-04-03T03:46:53.211706Z", "iopub.status.idle": "2026-04-03T03:46:54.730368Z", "shell.execute_reply": "2026-04-03T03:46:54.729614Z" } }, "outputs": [ { "data": { "text/html": [ "
MatrixProductOperator(tensors=20, indices=59, L=20, max_bond=166)
Tensor(shape=(4, 2, 2), inds=[_f713bfAAAAB, k0, b0], tags={I0}),backend=numpy, dtype=complex128, data=array([[[ 1.+0.j, 0.+0.j],\n", " [ 0.+0.j, -1.+0.j]],\n", "\n", " [[ 1.+0.j, 0.+0.j],\n", " [ 0.+0.j, 1.+0.j]],\n", "\n", " [[ 0.+0.j, 1.+0.j],\n", " [ 1.+0.j, 0.+0.j]],\n", "\n", " [[ 0.+0.j, 0.-1.j],\n", " [ 0.+1.j, 0.+0.j]]])
Tensor(shape=(4, 17, 2, 2), inds=[_f713bfAAAAB, _f713bfAAAAC, k1, b1], tags={I1}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(17, 33, 2, 2), inds=[_f713bfAAAAC, _f713bfAAAAD, k2, b2], tags={I2}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(33, 47, 2, 2), inds=[_f713bfAAAAD, _f713bfAAAAE, k3, b3], tags={I3}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(47, 67, 2, 2), inds=[_f713bfAAAAE, _f713bfAAAAF, k4, b4], tags={I4}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(67, 83, 2, 2), inds=[_f713bfAAAAF, _f713bfAAAAG, k5, b5], tags={I5}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(83, 89, 2, 2), inds=[_f713bfAAAAG, _f713bfAAAAH, k6, b6], tags={I6}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(89, 91, 2, 2), inds=[_f713bfAAAAH, _f713bfAAAAI, k7, b7], tags={I7}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(91, 106, 2, 2), inds=[_f713bfAAAAI, _f713bfAAAAJ, k8, b8], tags={I8}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(106, 104, 2, 2), inds=[_f713bfAAAAJ, _f713bfAAAAK, k9, b9], tags={I9}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(104, 127, 2, 2), inds=[_f713bfAAAAK, _f713bfAAAAL, k10, b10], tags={I10}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(127, 132, 2, 2), inds=[_f713bfAAAAL, _f713bfAAAAM, k11, b11], tags={I11}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(132, 141, 2, 2), inds=[_f713bfAAAAM, _f713bfAAAAN, k12, b12], tags={I12}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(141, 138, 2, 2), inds=[_f713bfAAAAN, _f713bfAAAAO, k13, b13], tags={I13}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(138, 152, 2, 2), inds=[_f713bfAAAAO, _f713bfAAAAP, k14, b14], tags={I14}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(152, 166, 2, 2), inds=[_f713bfAAAAP, _f713bfAAAAQ, k15, b15], tags={I15}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(166, 88, 2, 2), inds=[_f713bfAAAAQ, _f713bfAAAAR, k16, b16], tags={I16}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(88, 60, 2, 2), inds=[_f713bfAAAAR, _f713bfAAAAS, k17, b17], tags={I17}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(60, 5, 2, 2), inds=[_f713bfAAAAS, _f713bfAAAAT, k18, b18], tags={I18}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(5, 2, 2), inds=[_f713bfAAAAT, k19, b19], tags={I19}),backend=numpy, dtype=complex128, data=array([[[ 1. +0.j, 0. +0.j],\n", " [ 0. +0.j, 1. +0.j]],\n", "\n", " [[ 1. +0.j, 0. +0.j],\n", " [ 0. +0.j, -1. +0.j]],\n", "\n", " [[ 0. +0.j, 1. +0.j],\n", " [ 1. +0.j, 0. +0.j]],\n", "\n", " [[ 0. +0.j, 0. -1.j],\n", " [ 0. +1.j, 0. +0.j]],\n", "\n", " [[ 0.14679212+0.j, 0. +0.j],\n", " [ 0. +0.j, -0.14679212+0.j]]])
" ], "text/plain": [ "MatrixProductOperator(tensors=20, indices=59, L=20, max_bond=166)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "h2_mpo = h2_ham.to_mpo()\n", "display(h2_mpo)" ] }, { "cell_type": "markdown", "id": "0b114bb4", "metadata": {}, "source": [ "Here we use directly $\\texttt{quimb}$'s solver to introduce the different settings of a DMRG calculation. Note that the `do_dmrg` function of the `hamiltonian` module performs a basic DMRG calculation with default settings.\n", "\n", "In DMRG, the wavefunction is represented as an MPS. The energy is obtained as `(MPS.H) @ MPO @ MPS` where `@` indicates a contraction and `MPS.H` is the hermitian conjugate. The DMRG algorithm optimizes successively the different tensors of the MPS so as to minimize the energy. Once all the tensors have been optimized in a forward pass, the pass is carried backward: this forth-and-back optimization constitutes a full DMRG sweep. The maximum number of sweeps can be set with `max_sweeps`. Compression is performed along the algorithm, controlled by setting a maximal bond dimension and a cutoff. The `bond_dims` and `cutoffs` arguments take either a single value or a sequence of values (if one wants to assign different parameters to different sweeps).\n", "\n", "For an introduction on DMRG, see [this presentation](https://tensornetwork.org/mps/algorithms/dmrg/) and references therein." ] }, { "cell_type": "code", "execution_count": 5, "id": "cd14812d", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:46:54.732201Z", "iopub.status.busy": "2026-04-03T03:46:54.732024Z", "iopub.status.idle": "2026-04-03T03:47:17.407741Z", "shell.execute_reply": "2026-04-03T03:47:17.407024Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1, R, max_bond=(8/64), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-28.56424105232588-1.1142155991253344e-14j) ... not converged.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "2, R, max_bond=(16/64), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-29.124950365625956-1.109442399060967e-14j) ... not converged.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "3, R, max_bond=(32/64), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-29.125458255393504-3.953434801751143e-15j) ... not converged.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "4, R, max_bond=(58/64), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-29.125494836608233+6.4557734158476876e-15j) ... converged!\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1min 15s, sys: 762 ms, total: 1min 16s\n", "Wall time: 22.7 s\n" ] } ], "source": [ "%%time\n", "\n", "dmrg = DMRG2(h2_mpo)\n", "dmrg.solve(max_sweeps=16, bond_dims=64, verbosity=1, cutoffs=1e-12);" ] }, { "cell_type": "code", "execution_count": 6, "id": "0b2c9c7c", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:17.412288Z", "iopub.status.busy": "2026-04-03T03:47:17.412094Z", "iopub.status.idle": "2026-04-03T03:47:17.419319Z", "shell.execute_reply": "2026-04-03T03:47:17.418895Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6.084470541622977e-06\n" ] } ], "source": [ "e_dmrg = dmrg.energy + h2_ham.e_const\n", "print(abs(h2_ham.e_fci - e_dmrg))" ] }, { "cell_type": "markdown", "id": "27403fe6", "metadata": {}, "source": [ "## $O_2$\n", "\n", "Let us now consider a bigger molecule, dioxygen. We take the minimal STO-3G basis for simplicity.\n", "Since $O_2$ is an open-shell molecule, we must run Hartree-Fock in unrestricted mode. We must also specify the spin of the molecule since it is non trivial." ] }, { "cell_type": "code", "execution_count": 7, "id": "42709476", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:17.422681Z", "iopub.status.busy": "2026-04-03T03:47:17.422502Z", "iopub.status.idle": "2026-04-03T03:47:18.425438Z", "shell.execute_reply": "2026-04-03T03:47:18.424677Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged SCF energy = -147.63345273345 = 2.0034094 2S+1 = 3.0022721\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "nOrb : 10\n", "nElec : 16\n", "E_HF : -147.6334527334\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E_CI : -147.7415968576\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E(UCCSD) = -147.7396029424888 E_corr = -0.1061502090391701\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E_CCSD: -147.7396029425\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "n_qubits = 20\n" ] } ], "source": [ "basis = \"STO-3G\"\n", "\n", "mol = gto.M(\n", " atom=[(\"O\", (0.0, 0.0, 0.0)), (\"O\", (0.0, 0.0, 1.2))],\n", " basis=basis,\n", " spin=2,\n", ")\n", "\n", "o2_ham = chemistry_hamiltonian(mol, \"uhf\", do_fci=True, do_ccsd=True)\n", "print(f\"n_qubits = {o2_ham.n_qubits}\")" ] }, { "cell_type": "markdown", "id": "d4dfef66", "metadata": {}, "source": [ "Let us build the Hamiltonian MPO representation and compress it to reduce its bond dimension (hence reducing the cost of running DMRG)." ] }, { "cell_type": "code", "execution_count": 8, "id": "d5fe824d", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:18.427272Z", "iopub.status.busy": "2026-04-03T03:47:18.427067Z", "iopub.status.idle": "2026-04-03T03:47:19.975448Z", "shell.execute_reply": "2026-04-03T03:47:19.974768Z" } }, "outputs": [ { "data": { "text/html": [ "
MatrixProductOperator(tensors=20, indices=59, L=20, max_bond=156)
Tensor(shape=(4, 2, 2), inds=[_f713bfAAACz, k0, b0], tags={I0}),backend=numpy, dtype=complex128, data=array([[[ 1.+0.j, 0.+0.j],\n", " [ 0.+0.j, -1.+0.j]],\n", "\n", " [[ 1.+0.j, 0.+0.j],\n", " [ 0.+0.j, 1.+0.j]],\n", "\n", " [[ 0.+0.j, 1.+0.j],\n", " [ 1.+0.j, 0.+0.j]],\n", "\n", " [[ 0.+0.j, 0.-1.j],\n", " [ 0.+1.j, 0.+0.j]]])
Tensor(shape=(4, 17, 2, 2), inds=[_f713bfAAACz, _f713bfAAADA, k1, b1], tags={I1}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(17, 33, 2, 2), inds=[_f713bfAAADA, _f713bfAAADB, k2, b2], tags={I2}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(33, 47, 2, 2), inds=[_f713bfAAADB, _f713bfAAADC, k3, b3], tags={I3}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(47, 67, 2, 2), inds=[_f713bfAAADC, _f713bfAAADD, k4, b4], tags={I4}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(67, 83, 2, 2), inds=[_f713bfAAADD, _f713bfAAADE, k5, b5], tags={I5}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(83, 89, 2, 2), inds=[_f713bfAAADE, _f713bfAAADF, k6, b6], tags={I6}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(89, 92, 2, 2), inds=[_f713bfAAADF, _f713bfAAADG, k7, b7], tags={I7}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(92, 106, 2, 2), inds=[_f713bfAAADG, _f713bfAAADH, k8, b8], tags={I8}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(106, 104, 2, 2), inds=[_f713bfAAADH, _f713bfAAADI, k9, b9], tags={I9}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(104, 127, 2, 2), inds=[_f713bfAAADI, _f713bfAAADJ, k10, b10], tags={I10}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(127, 132, 2, 2), inds=[_f713bfAAADJ, _f713bfAAADK, k11, b11], tags={I11}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(132, 143, 2, 2), inds=[_f713bfAAADK, _f713bfAAADL, k12, b12], tags={I12}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(143, 138, 2, 2), inds=[_f713bfAAADL, _f713bfAAADM, k13, b13], tags={I13}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(138, 144, 2, 2), inds=[_f713bfAAADM, _f713bfAAADN, k14, b14], tags={I14}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(144, 156, 2, 2), inds=[_f713bfAAADN, _f713bfAAADO, k15, b15], tags={I15}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(156, 88, 2, 2), inds=[_f713bfAAADO, _f713bfAAADP, k16, b16], tags={I16}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(88, 66, 2, 2), inds=[_f713bfAAADP, _f713bfAAADQ, k17, b17], tags={I17}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(66, 5, 2, 2), inds=[_f713bfAAADQ, _f713bfAAADR, k18, b18], tags={I18}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(5, 2, 2), inds=[_f713bfAAADR, k19, b19], tags={I19}),backend=numpy, dtype=complex128, data=array([[[ 1. +0.j, 0. +0.j],\n", " [ 0. +0.j, 1. +0.j]],\n", "\n", " [[ 1. +0.j, 0. +0.j],\n", " [ 0. +0.j, -1. +0.j]],\n", "\n", " [[ 0. +0.j, 1. +0.j],\n", " [ 1. +0.j, 0. +0.j]],\n", "\n", " [[ 0. +0.j, 0. -1.j],\n", " [ 0. +1.j, 0. +0.j]],\n", "\n", " [[ 0.15462394+0.j, 0. +0.j],\n", " [ 0. +0.j, -0.15462394+0.j]]])
" ], "text/plain": [ "MatrixProductOperator(tensors=20, indices=59, L=20, max_bond=156)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "o2_mpo_original = o2_ham.to_mpo()\n", "display(o2_mpo_original)" ] }, { "cell_type": "markdown", "id": "2399ee63", "metadata": {}, "source": [ "To make sure that we don't lose physical information when compressing, we measure the MPO norm" ] }, { "cell_type": "code", "execution_count": 9, "id": "0f6495c9", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:19.977367Z", "iopub.status.busy": "2026-04-03T03:47:19.977145Z", "iopub.status.idle": "2026-04-03T03:47:20.008503Z", "shell.execute_reply": "2026-04-03T03:47:20.007453Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "28188.084034552536\n" ] } ], "source": [ "norm = o2_mpo_original.norm()\n", "print(norm)" ] }, { "cell_type": "markdown", "id": "48561d32", "metadata": {}, "source": [ "Let us compress using a cutoff (alternatively, one can set a maximal bond dimension with the option `max_bond`)." ] }, { "cell_type": "code", "execution_count": 10, "id": "95f2e3f5", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:20.010613Z", "iopub.status.busy": "2026-04-03T03:47:20.010444Z", "iopub.status.idle": "2026-04-03T03:47:21.130161Z", "shell.execute_reply": "2026-04-03T03:47:21.129684Z" } }, "outputs": [ { "data": { "text/html": [ "
MatrixProductOperator(tensors=20, indices=59, L=20, max_bond=51)
Tensor(shape=(4, 2, 2), inds=[_f713bfAAACz, k0, b0], tags={I0}),backend=numpy, dtype=complex128, data=array([[[-1.85428891e+04+6.06143346e-15j,\n", " 2.98447151e-15-1.32158765e-14j],\n", " [ 8.26574809e-15-1.32525034e-15j,\n", " -1.64513560e+04+6.06143346e-15j]],\n", "\n", " [[ 8.86628638e+03+1.95246622e-16j,\n", " 6.97938055e-16-1.44001441e-15j],\n", " [-5.30118220e-16-1.45867397e-16j,\n", " -9.99349634e+03+1.95246622e-16j]],\n", "\n", " [[ 9.32501954e-17+5.39520448e-17j,\n", " 3.76486550e+02-6.34550978e+02j],\n", " [-4.98205251e+02-6.45486251e+01j,\n", " 2.27831474e-16+5.39520448e-17j]],\n", "\n", " [[ 2.00919160e-16-3.11856220e-16j,\n", " -3.47618871e+02-3.62679083e+02j],\n", " [ 2.87295045e+02-6.79601813e+02j,\n", " 2.00919160e-16-3.11856220e-16j]]])
Tensor(shape=(4, 7, 2, 2), inds=[_f713bfAAACz, _f713bfAAADA, k1, b1], tags={I1}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(7, 14, 2, 2), inds=[_f713bfAAADA, _f713bfAAADB, k2, b2], tags={I2}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(14, 21, 2, 2), inds=[_f713bfAAADB, _f713bfAAADC, k3, b3], tags={I3}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(21, 32, 2, 2), inds=[_f713bfAAADC, _f713bfAAADD, k4, b4], tags={I4}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(32, 51, 2, 2), inds=[_f713bfAAADD, _f713bfAAADE, k5, b5], tags={I5}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(51, 46, 2, 2), inds=[_f713bfAAADE, _f713bfAAADF, k6, b6], tags={I6}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(46, 39, 2, 2), inds=[_f713bfAAADF, _f713bfAAADG, k7, b7], tags={I7}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(39, 38, 2, 2), inds=[_f713bfAAADG, _f713bfAAADH, k8, b8], tags={I8}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(38, 20, 2, 2), inds=[_f713bfAAADH, _f713bfAAADI, k9, b9], tags={I9}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(20, 37, 2, 2), inds=[_f713bfAAADI, _f713bfAAADJ, k10, b10], tags={I10}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(37, 35, 2, 2), inds=[_f713bfAAADJ, _f713bfAAADK, k11, b11], tags={I11}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(35, 35, 2, 2), inds=[_f713bfAAADK, _f713bfAAADL, k12, b12], tags={I12}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(35, 37, 2, 2), inds=[_f713bfAAADL, _f713bfAAADM, k13, b13], tags={I13}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(37, 24, 2, 2), inds=[_f713bfAAADM, _f713bfAAADN, k14, b14], tags={I14}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(24, 25, 2, 2), inds=[_f713bfAAADN, _f713bfAAADO, k15, b15], tags={I15}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(25, 23, 2, 2), inds=[_f713bfAAADO, _f713bfAAADP, k16, b16], tags={I16}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(23, 12, 2, 2), inds=[_f713bfAAADP, _f713bfAAADQ, k17, b17], tags={I17}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(12, 4, 2, 2), inds=[_f713bfAAADQ, _f713bfAAADR, k18, b18], tags={I18}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(4, 2, 2), inds=[_f713bfAAADR, k19, b19], tags={I19}),backend=numpy, dtype=complex128, data=array([[[-7.22793422e-01+0.00000000e+00j,\n", " 0.00000000e+00+1.63944940e-40j],\n", " [ 1.02532712e-39-1.66797690e-40j,\n", " -6.91064157e-01+2.90681316e-21j]],\n", "\n", " [[-6.91064157e-01+0.00000000e+00j,\n", " 0.00000000e+00-1.71472247e-40j],\n", " [-1.07240361e-39+1.74455978e-40j,\n", " 7.22793422e-01-3.04027551e-21j]],\n", "\n", " [[-0.00000000e+00+0.00000000e+00j,\n", " 9.88172485e-01+1.53346469e-01j],\n", " [-7.25932805e-06-1.70388667e-05j,\n", " 0.00000000e+00+0.00000000e+00j]],\n", "\n", " [[ 0.00000000e+00+0.00000000e+00j,\n", " -7.79038627e-26+1.85208213e-05j],\n", " [-8.48998360e-01+5.28395480e-01j,\n", " -1.38901188e-39+0.00000000e+00j]]])
" ], "text/plain": [ "MatrixProductOperator(tensors=20, indices=59, L=20, max_bond=51)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "o2_mpo = o2_mpo_original.copy()\n", "o2_mpo.compress(cutoff=1e-6)\n", "display(o2_mpo)" ] }, { "cell_type": "markdown", "id": "abb60f3f", "metadata": {}, "source": [ "We then make sure that the compression doesn't affect the norm and measure the overlap with the original uncompressed MPO.\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "f9447b49", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:21.137570Z", "iopub.status.busy": "2026-04-03T03:47:21.137392Z", "iopub.status.idle": "2026-04-03T03:47:21.218678Z", "shell.execute_reply": "2026-04-03T03:47:21.218015Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.9999951650744955\n", "(0.9999951650744959+3.800224397730763e-19j)\n" ] } ], "source": [ "print(o2_mpo.norm() / norm)\n", "print(np.sqrt(o2_mpo.H @ o2_mpo_original) / norm)" ] }, { "cell_type": "markdown", "id": "a764846a", "metadata": {}, "source": [ "We now run DMRG" ] }, { "cell_type": "code", "execution_count": 12, "id": "2f7f55f5", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:21.220805Z", "iopub.status.busy": "2026-04-03T03:47:21.220134Z", "iopub.status.idle": "2026-04-03T03:47:28.789831Z", "shell.execute_reply": "2026-04-03T03:47:28.789317Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1, R, max_bond=(8/200), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-58.283859790056184-6.8863318308243215e-06j) ... not converged.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "2, R, max_bond=(16/200), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-58.29091362411727-3.867180399019787e-06j) ... not converged.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "3, R, max_bond=(32/200), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-58.29103581768155-3.3722779117228985e-06j) ... not converged.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "4, R, max_bond=(60/200), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-58.29103805557315-3.376817343414551e-06j) ... converged!\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 30.1 s, sys: 7.77 ms, total: 30.1 s\n", "Wall time: 7.56 s\n" ] } ], "source": [ "%%time\n", "dmrg = DMRG2(o2_mpo)\n", "dmrg.solve(max_sweeps=16, bond_dims=[200], verbosity=1, cutoffs=1e-12);" ] }, { "cell_type": "markdown", "id": "76ef1bdf", "metadata": {}, "source": [ "Below we print the results and computation time for running DMRG in the larger 6-31G basis set.\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "ed821ccc", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:28.794851Z", "iopub.status.busy": "2026-04-03T03:47:28.794671Z", "iopub.status.idle": "2026-04-03T03:47:28.800016Z", "shell.execute_reply": "2026-04-03T03:47:28.798988Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6-31G DMRG energy with max_bond 128 = -149.78252800735535 in 2 hours w MPO bond dim = 660\n", "6-31G DMRG energy with max bond 200 = -149.78507049242708 in 1 hour w MPO bond dim = 190\n", "STO-3G DMRG energy = -147.73948105719143\n", "STO-3G UCCSD energy = -147.7396029424888\n" ] } ], "source": [ "print(\n", " \"6-31G DMRG energy with max_bond 128 = -149.78252800735535 in 2 hours w MPO bond dim = 660\"\n", ")\n", "print(\n", " \"6-31G DMRG energy with max bond 200 = -149.78507049242708 in 1 hour w MPO bond dim = 190\"\n", ")\n", "# Current results\n", "print(f\"{basis} DMRG energy = {np.real(dmrg.energy) + o2_ham.e_const}\")\n", "print(f\"{basis} UCCSD energy = {o2_ham.e_ccsd}\")" ] }, { "cell_type": "markdown", "id": "180a9c53", "metadata": {}, "source": [ "## $H_2O$\n", "\n", "We are of course not restricted to molecules with a single atomic element. A lot of molecules' descriptions including geometry, spin, charge, symmetry arguments compatible with $\\texttt{pyscf}$ can be found on the internet. We are only limited by computational power. In the following we take the water molecule in the minimal STO-3G basis." ] }, { "cell_type": "code", "execution_count": 14, "id": "42738860", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:28.801666Z", "iopub.status.busy": "2026-04-03T03:47:28.801513Z", "iopub.status.idle": "2026-04-03T03:47:29.282589Z", "shell.execute_reply": "2026-04-03T03:47:29.281904Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "converged SCF energy = -74.963146775618\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "nOrb : 7\n", "nElec : 10\n", "E_HF : -74.9631467756\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E_CI : -75.0127761764\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E(CCSD) = -75.0126602641099 E_corr = -0.04951348849192778\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "E_CCSD: -75.0126602641\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "n_qubits = 14\n" ] } ], "source": [ "basis = \"STO-3G\"\n", "mol = gto.M(\n", " atom=[\n", " (\"O\", (0.0, 0.0, 0.117790)),\n", " (\"H\", (0.0, 0.755453, -0.471161)),\n", " (\"H\", (0.0, -0.755453, -0.471161)),\n", " ],\n", " basis=basis,\n", ")\n", "\n", "h2o_ham = chemistry_hamiltonian(mol, \"rhf\", do_fci=True, do_ccsd=True)\n", "print(f\"n_qubits = {h2o_ham.n_qubits}\")" ] }, { "cell_type": "code", "execution_count": 15, "id": "3d778d5c", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:29.284400Z", "iopub.status.busy": "2026-04-03T03:47:29.284245Z", "iopub.status.idle": "2026-04-03T03:47:29.531955Z", "shell.execute_reply": "2026-04-03T03:47:29.530986Z" } }, "outputs": [ { "data": { "text/html": [ "
MatrixProductOperator(tensors=14, indices=41, L=14, max_bond=83)
Tensor(shape=(4, 2, 2), inds=[_f713bfAAAHC, k0, b0], tags={I0}),backend=numpy, dtype=complex128, data=array([[[ 1.+0.j, 0.+0.j],\n", " [ 0.+0.j, -1.+0.j]],\n", "\n", " [[ 1.+0.j, 0.+0.j],\n", " [ 0.+0.j, 1.+0.j]],\n", "\n", " [[ 0.+0.j, 1.+0.j],\n", " [ 1.+0.j, 0.+0.j]],\n", "\n", " [[ 0.+0.j, 0.-1.j],\n", " [ 0.+1.j, 0.+0.j]]])
Tensor(shape=(4, 17, 2, 2), inds=[_f713bfAAAHC, _f713bfAAAHD, k1, b1], tags={I1}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(17, 33, 2, 2), inds=[_f713bfAAAHD, _f713bfAAAHE, k2, b2], tags={I2}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(33, 38, 2, 2), inds=[_f713bfAAAHE, _f713bfAAAHF, k3, b3], tags={I3}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(38, 46, 2, 2), inds=[_f713bfAAAHF, _f713bfAAAHG, k4, b4], tags={I4}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(46, 54, 2, 2), inds=[_f713bfAAAHG, _f713bfAAAHH, k5, b5], tags={I5}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(54, 52, 2, 2), inds=[_f713bfAAAHH, _f713bfAAAHI, k6, b6], tags={I6}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(52, 69, 2, 2), inds=[_f713bfAAAHI, _f713bfAAAHJ, k7, b7], tags={I7}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(69, 74, 2, 2), inds=[_f713bfAAAHJ, _f713bfAAAHK, k8, b8], tags={I8}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(74, 83, 2, 2), inds=[_f713bfAAAHK, _f713bfAAAHL, k9, b9], tags={I9}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(83, 68, 2, 2), inds=[_f713bfAAAHL, _f713bfAAAHM, k10, b10], tags={I10}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(68, 48, 2, 2), inds=[_f713bfAAAHM, _f713bfAAAHN, k11, b11], tags={I11}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(48, 5, 2, 2), inds=[_f713bfAAAHN, _f713bfAAAHO, k12, b12], tags={I12}),backend=numpy, dtype=complex128, data=...
Tensor(shape=(5, 2, 2), inds=[_f713bfAAAHO, k13, b13], tags={I13}),backend=numpy, dtype=complex128, data=array([[[ 1. +0.j, 0. +0.j],\n", " [ 0. +0.j, 1. +0.j]],\n", "\n", " [[ 1. +0.j, 0. +0.j],\n", " [ 0. +0.j, -1. +0.j]],\n", "\n", " [[ 0. +0.j, 1. +0.j],\n", " [ 1. +0.j, 0. +0.j]],\n", "\n", " [[ 0. +0.j, 0. -1.j],\n", " [ 0. +1.j, 0. +0.j]],\n", "\n", " [[ 0.112834+0.j, 0. +0.j],\n", " [ 0. +0.j, -0.112834+0.j]]])
" ], "text/plain": [ "MatrixProductOperator(tensors=14, indices=41, L=14, max_bond=83)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "h2o_mpo = h2o_ham.to_mpo()\n", "display(h2o_mpo)" ] }, { "cell_type": "code", "execution_count": 16, "id": "4cf96f6d", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:29.533624Z", "iopub.status.busy": "2026-04-03T03:47:29.533300Z", "iopub.status.idle": "2026-04-03T03:47:32.306029Z", "shell.execute_reply": "2026-04-03T03:47:32.305519Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1, R, max_bond=(8/64), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-28.58687041704924-1.539220140234221e-14j) ... not converged.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "2, R, max_bond=(16/64), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-28.58984744190927-5.35682609381638e-15j) ... not converged.\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "3, R, max_bond=(32/64), cutoff:1e-12\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Energy: (-28.589859362032215+1.6112111644872584e-14j) ... converged!\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 11 s, sys: 976 μs, total: 11 s\n", "Wall time: 2.77 s\n" ] } ], "source": [ "%%time\n", "dmrg = DMRG2(h2o_mpo)\n", "dmrg.solve(max_sweeps=16, bond_dims=64, verbosity=1, cutoffs=1e-12);" ] }, { "cell_type": "code", "execution_count": 17, "id": "b24cc925", "metadata": { "execution": { "iopub.execute_input": "2026-04-03T03:47:32.308836Z", "iopub.status.busy": "2026-04-03T03:47:32.307989Z", "iopub.status.idle": "2026-04-03T03:47:32.313435Z", "shell.execute_reply": "2026-04-03T03:47:32.312559Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "===== H2O -- STO-3G =====\n", "CCSD energy = -75.0126602641\n", "FCI energy = -75.0127761764\n", "DMRG energy = -75.0127759531\n" ] } ], "source": [ "print(f\"===== H2O -- {basis} =====\")\n", "print(f\"CCSD energy = {h2o_ham.e_ccsd:.10f}\")\n", "print(f\"FCI energy = {h2o_ham.e_fci:.10f}\")\n", "print(f\"DMRG energy = {np.real(dmrg.energy + h2o_ham.e_const):.10f}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "b29764c9", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.12" } }, "nbformat": 4, "nbformat_minor": 5 }