2. Chemistry to qubit Hamiltonians¶

This tutorial describes how to load a molecular Hamiltonian from \(\texttt{pyscf}\) and convert it to a qubit Hamiltonian object of the \(\texttt{qpe-toolbox}\)’s Hamiltonian class. The Hamiltonian is designed to be compatible with \(\texttt{quimb}\).

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.

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

import numpy as np
from IPython.display import display
from pyscf import gto
from quimb.tensor import DMRG2

from qpe_toolbox.hamiltonian import chemistry_hamiltonian

The chemistry_hamiltonian function is built on the Hamiltonian class. It takes a molecule described with a Mole object from \(\texttt{pyscf}\), performs the Hartree-Fock calculation and converts the molecular Hamiltonian into a Hamiltonian object describing the qubit Hamiltonian.

2.1. \(H_2\)¶

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\)). Then we call the chemistry_hamiltonian.

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 ("uhf") uses different molecular orbitals for the different eigenstates of \(S_z\). Unrestricted Hartree-Fock is mostly used for open-shell molecules.
Optionally, the FCI or CCSD energies can be computed.

mol = gto.M(
    atom=[("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.735))],
    basis="sto-3g",
)

h2_ham_sto3g = chemistry_hamiltonian(mol, "rhf", do_fci=True, do_ccsd=True)
print(f"n_qubits = {h2_ham_sto3g.n_qubits}")

converged SCF energy = -1.116998996754

nOrb   : 2
nElec  : 2
E_HF   : -1.1169989968

E_CI   : -1.1373060358

E(CCSD) = -1.137306193391969  E_corr = -0.02030719663796497

E_CCSD: -1.1373061934

n_qubits = 4

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.

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. In chemistry, an energy estimation must reach the chemical accuracy standard, defined to be \(1.6 \text{mHa} \approx 500 \text{K}\). The HF energy is far from the FCI energy by \(20 \text{mHa}\), while the CCSD energy is close within \(10^{-6} \text{Ha}\).

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) for \(H_2\) are around \(-1.16\) Ha. Even with exact diagonalization, the accuracy is at best \(30 \text{mHa}\) in this basis. Hence, reaching chemical accuracy will always mean using much larger basis sets. 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.

In the following, we consider the cc-pvdz basis set with \(5\) orbitals per Hydrogen atom.

mol = gto.M(
    atom=[("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.735))],
    basis="cc-pvdz",
)

h2_ham = chemistry_hamiltonian(mol, "rhf", do_fci=True, do_ccsd=True)
print(f"n_qubits = {h2_ham.n_qubits}")

converged SCF energy = -1.12862276992659

nOrb   : 10
nElec  : 2
E_HF   : -1.1286227699

E_CI   : -1.1632095576

E(CCSD) = -1.163209563562649  E_corr = -0.03458679363605879

E_CCSD: -1.1632095636

n_qubits = 20

The FCI energy is now \(-1.1632 \text{Ha}\), within chemical accuracy compared to the reported value of \(-1.1640 \text{Ha}\).

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.

h2_mpo = h2_ham.to_mpo()
display(h2_mpo)

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], [ 0.+0.j, -1.+0.j]], [[ 1.+0.j, 0.+0.j], [ 0.+0.j, 1.+0.j]], [[ 0.+0.j, 1.+0.j], [ 1.+0.j, 0.+0.j]], [[ 0.+0.j, 0.-1.j], [ 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], [ 0. +0.j, 1. +0.j]], [[ 1. +0.j, 0. +0.j], [ 0. +0.j, -1. +0.j]], [[ 0. +0.j, 1. +0.j], [ 1. +0.j, 0. +0.j]], [[ 0. +0.j, 0. -1.j], [ 0. +1.j, 0. +0.j]], [[ 0.14679212+0.j, 0. +0.j], [ 0. +0.j, -0.14679212+0.j]]])

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.

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

For an introduction on DMRG, see this presentation and references therein.

%%time

dmrg = DMRG2(h2_mpo)
dmrg.solve(max_sweeps=16, bond_dims=64, verbosity=1, cutoffs=1e-12);

1, R, max_bond=(8/64), cutoff:1e-12

Energy: (-28.56424105232588-1.1142155991253344e-14j) ... not converged.

2, R, max_bond=(16/64), cutoff:1e-12

Energy: (-29.124950365625956-1.109442399060967e-14j) ... not converged.

3, R, max_bond=(32/64), cutoff:1e-12

Energy: (-29.125458255393504-3.953434801751143e-15j) ... not converged.

4, R, max_bond=(58/64), cutoff:1e-12

Energy: (-29.125494836608233+6.4557734158476876e-15j) ... converged!

CPU times: user 1min 15s, sys: 762 ms, total: 1min 16s
Wall time: 22.7 s

e_dmrg = dmrg.energy + h2_ham.e_const
print(abs(h2_ham.e_fci - e_dmrg))

6.084470541622977e-06

2.2. \(O_2\)¶

Let us now consider a bigger molecule, dioxygen. We take the minimal STO-3G basis for simplicity. 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.

basis = "STO-3G"

mol = gto.M(
    atom=[("O", (0.0, 0.0, 0.0)), ("O", (0.0, 0.0, 1.2))],
    basis=basis,
    spin=2,
)

o2_ham = chemistry_hamiltonian(mol, "uhf", do_fci=True, do_ccsd=True)
print(f"n_qubits = {o2_ham.n_qubits}")

converged SCF energy = -147.63345273345  <S^2> = 2.0034094  2S+1 = 3.0022721

nOrb   : 10
nElec  : 16
E_HF   : -147.6334527334

E_CI   : -147.7415968576

E(UCCSD) = -147.7396029424888  E_corr = -0.1061502090391701

E_CCSD: -147.7396029425

n_qubits = 20

Let us build the Hamiltonian MPO representation and compress it to reduce its bond dimension (hence reducing the cost of running DMRG).

o2_mpo_original = o2_ham.to_mpo()
display(o2_mpo_original)

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], [ 0.+0.j, -1.+0.j]], [[ 1.+0.j, 0.+0.j], [ 0.+0.j, 1.+0.j]], [[ 0.+0.j, 1.+0.j], [ 1.+0.j, 0.+0.j]], [[ 0.+0.j, 0.-1.j], [ 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], [ 0. +0.j, 1. +0.j]], [[ 1. +0.j, 0. +0.j], [ 0. +0.j, -1. +0.j]], [[ 0. +0.j, 1. +0.j], [ 1. +0.j, 0. +0.j]], [[ 0. +0.j, 0. -1.j], [ 0. +1.j, 0. +0.j]], [[ 0.15462394+0.j, 0. +0.j], [ 0. +0.j, -0.15462394+0.j]]])

To make sure that we don’t lose physical information when compressing, we measure the MPO norm

norm = o2_mpo_original.norm()
print(norm)

28188.084034552536

Let us compress using a cutoff (alternatively, one can set a maximal bond dimension with the option max_bond).

o2_mpo = o2_mpo_original.copy()
o2_mpo.compress(cutoff=1e-6)
display(o2_mpo)

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, 2.98447151e-15-1.32158765e-14j], [ 8.26574809e-15-1.32525034e-15j, -1.64513560e+04+6.06143346e-15j]], [[ 8.86628638e+03+1.95246622e-16j, 6.97938055e-16-1.44001441e-15j], [-5.30118220e-16-1.45867397e-16j, -9.99349634e+03+1.95246622e-16j]], [[ 9.32501954e-17+5.39520448e-17j, 3.76486550e+02-6.34550978e+02j], [-4.98205251e+02-6.45486251e+01j, 2.27831474e-16+5.39520448e-17j]], [[ 2.00919160e-16-3.11856220e-16j, -3.47618871e+02-3.62679083e+02j], [ 2.87295045e+02-6.79601813e+02j, 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, 0.00000000e+00+1.63944940e-40j], [ 1.02532712e-39-1.66797690e-40j, -6.91064157e-01+2.90681316e-21j]], [[-6.91064157e-01+0.00000000e+00j, 0.00000000e+00-1.71472247e-40j], [-1.07240361e-39+1.74455978e-40j, 7.22793422e-01-3.04027551e-21j]], [[-0.00000000e+00+0.00000000e+00j, 9.88172485e-01+1.53346469e-01j], [-7.25932805e-06-1.70388667e-05j, 0.00000000e+00+0.00000000e+00j]], [[ 0.00000000e+00+0.00000000e+00j, -7.79038627e-26+1.85208213e-05j], [-8.48998360e-01+5.28395480e-01j, -1.38901188e-39+0.00000000e+00j]]])

We then make sure that the compression doesn’t affect the norm and measure the overlap with the original uncompressed MPO.

print(o2_mpo.norm() / norm)
print(np.sqrt(o2_mpo.H @ o2_mpo_original) / norm)

0.9999951650744955
(0.9999951650744959+3.800224397730763e-19j)

We now run DMRG

%%time
dmrg = DMRG2(o2_mpo)
dmrg.solve(max_sweeps=16, bond_dims=[200], verbosity=1, cutoffs=1e-12);

1, R, max_bond=(8/200), cutoff:1e-12

Energy: (-58.283859790056184-6.8863318308243215e-06j) ... not converged.

2, R, max_bond=(16/200), cutoff:1e-12

Energy: (-58.29091362411727-3.867180399019787e-06j) ... not converged.

3, R, max_bond=(32/200), cutoff:1e-12

Energy: (-58.29103581768155-3.3722779117228985e-06j) ... not converged.

4, R, max_bond=(60/200), cutoff:1e-12

Energy: (-58.29103805557315-3.376817343414551e-06j) ... converged!

CPU times: user 30.1 s, sys: 7.77 ms, total: 30.1 s
Wall time: 7.56 s

Below we print the results and computation time for running DMRG in the larger 6-31G basis set.

print(
    "6-31G DMRG energy with max_bond 128 = -149.78252800735535 in 2 hours w MPO bond dim = 660"
)
print(
    "6-31G DMRG energy with max bond 200 = -149.78507049242708 in 1 hour w MPO bond dim = 190"
)
# Current results
print(f"{basis} DMRG energy = {np.real(dmrg.energy) + o2_ham.e_const}")
print(f"{basis} UCCSD energy = {o2_ham.e_ccsd}")

6-31G DMRG energy with max_bond 128 = -149.78252800735535 in 2 hours w MPO bond dim = 660
6-31G DMRG energy with max bond 200 = -149.78507049242708 in 1 hour w MPO bond dim = 190
STO-3G DMRG energy = -147.73948105719143
STO-3G UCCSD energy = -147.7396029424888

2.3. \(H_2O\)¶

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.

basis = "STO-3G"
mol = gto.M(
    atom=[
        ("O", (0.0, 0.0, 0.117790)),
        ("H", (0.0, 0.755453, -0.471161)),
        ("H", (0.0, -0.755453, -0.471161)),
    ],
    basis=basis,
)

h2o_ham = chemistry_hamiltonian(mol, "rhf", do_fci=True, do_ccsd=True)
print(f"n_qubits = {h2o_ham.n_qubits}")

converged SCF energy = -74.963146775618

nOrb   : 7
nElec  : 10
E_HF   : -74.9631467756

E_CI   : -75.0127761764

E(CCSD) = -75.0126602641099  E_corr = -0.04951348849192778

E_CCSD: -75.0126602641

n_qubits = 14

h2o_mpo = h2o_ham.to_mpo()
display(h2o_mpo)

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], [ 0.+0.j, -1.+0.j]], [[ 1.+0.j, 0.+0.j], [ 0.+0.j, 1.+0.j]], [[ 0.+0.j, 1.+0.j], [ 1.+0.j, 0.+0.j]], [[ 0.+0.j, 0.-1.j], [ 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], [ 0. +0.j, 1. +0.j]], [[ 1. +0.j, 0. +0.j], [ 0. +0.j, -1. +0.j]], [[ 0. +0.j, 1. +0.j], [ 1. +0.j, 0. +0.j]], [[ 0. +0.j, 0. -1.j], [ 0. +1.j, 0. +0.j]], [[ 0.112834+0.j, 0. +0.j], [ 0. +0.j, -0.112834+0.j]]])

%%time
dmrg = DMRG2(h2o_mpo)
dmrg.solve(max_sweeps=16, bond_dims=64, verbosity=1, cutoffs=1e-12);

1, R, max_bond=(8/64), cutoff:1e-12

Energy: (-28.58687041704924-1.539220140234221e-14j) ... not converged.

2, R, max_bond=(16/64), cutoff:1e-12

Energy: (-28.58984744190927-5.35682609381638e-15j) ... not converged.

3, R, max_bond=(32/64), cutoff:1e-12

Energy: (-28.589859362032215+1.6112111644872584e-14j) ... converged!

CPU times: user 11 s, sys: 976 μs, total: 11 s
Wall time: 2.77 s

print(f"===== H2O -- {basis} =====")
print(f"CCSD energy = {h2o_ham.e_ccsd:.10f}")
print(f"FCI energy = {h2o_ham.e_fci:.10f}")
print(f"DMRG energy = {np.real(dmrg.energy + h2o_ham.e_const):.10f}")

===== H2O -- STO-3G =====
CCSD energy = -75.0126602641
FCI energy = -75.0127761764
DMRG energy = -75.0127759531