{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Basis-Set-Free VQEs with the Tequila - Madness interface"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This tutorial covers the basic usage of the `tequila`-`madness` interface that allows to construct basis-set-free qubit Hamiltonians described in [arxiv:2008.02819](https://arxiv.org/abs/2008.02819/) using the MP2-PNO surrogate described in [doi:10.1063/1.5141880](https://doi.org/10.1063/1.5141880).\n",
"\n",
"## Content\n",
"\n",
"* [Installation](#installation) \n",
"* [A First Example: The Hydrogen Molecule](#first_example)\n",
" * [Compute More Orbitals](#first_example)\n",
" * [VQE: Separable Pair Approximation](#spa)\n",
" * [VQE: Other Methods](#vqe)\n",
"* [A Serious Example: The BeH2 Molecule](#second_example)\n",
"* [Data I/O](#data_io)\n",
"* [For Madness Experts](#madness_experts)\n",
"* [More Information](#more)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Installation \n",
"\n",
"\n",
"In order to use madness as chemistry backend in tequila we need a slightly altered version. \n",
"We can conveniently install it over the anaconda cloud with:\n",
"```bash\n",
"conda install madtequila -c kottmann\n",
"```\n",
"Unfortunately this is currently only supported for Linux-64. Mac users can however resort to manual compilation: follow the instructions from this [fork](https://github.com/kottmanj/madness). \n",
"Madness is quite sensitive about changes in it's dependencies (e.g. MKL or MPICH), so we recommend to install this in a fresh environment without many other packages. See the last cell of this notebook for an example."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### A First Example\n",
"\n",
"You should see output like\n",
"```bash\n",
"Starting madness calculation with executable: /PATH/TO/bin/pno_integrals\n",
"output redirected to he_pno_integrals.out logfile\n",
"```\n",
"and the time should be a few seconds (~10 seconds on an intel i5 processor)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals\n",
"output redirected to h2_pno_integrals.out logfile\n",
"finished after 9.780884265899658s\n",
"took 9.79s\n",
"created a qubit Hamiltonian with 4 qubits\n",
"basis_type : custom\n",
"basis_name : custom\n",
"orthogonal : True\n",
"functions : 2\n",
"reference : [0]\n",
"Current Orbitals\n",
"{idx_total:0, idx:0, occ:2.0, pair:(0, 0)}\n",
"coefficients: [1. 0.]\n",
"{idx_total:1, idx:1, occ:0.0425791, pair:(0, 0)}\n",
"coefficients: [0. 1.]\n"
]
}
],
"source": [
"import tequila as tq\n",
"import time\n",
"\n",
"geometry= \"H 0.0 0.0 0.0\\nH 0.0 0.0 1.5\"\n",
"start=time.time()\n",
"h2_2_4 = tq.Molecule(geometry=geometry)\n",
"print(\"took {:4.2f}s\".format(time.time()-start))\n",
"H = h2_2_4.make_hamiltonian()\n",
"print(\"created a qubit Hamiltonian with {} qubits\".format(H.n_qubits))\n",
"h2_2_4.print_basis_info()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the previous cell we created an H2 Molecule in a system adapted orbital basis. Note that we did not specify how many orbitals shall be computed. The default for this is: N-Orbitals = N-electrons. The H2 molecule has 2 electrons, so 2 orbitals (leading to 4 spin-orbitals/qubits) were computed. \n",
"Following the notation of [arxiv:2008.02819](https://arxiv.org/abs/2008.02819/) we denote this with H2/MRA-PNO(2,4).\n",
"\n",
"The printed information shows us, that we have one Hartree-Fock orbital with occupation number 2.0 and one pair-natural orbital (PNO) with occupation number 0.04. These occupation numbers are from the surrogate model (MP2-PNO) which is based on perturbation theory. From the 0.04 occupation number in the PNO we see that the perturbation in the surrogate is not large meaning that it is justified. We can therefore expect, that these two orbitals are close to optimal."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### VQE: Separable Pair Approximation\n",
"\n",
"Through the 'mol.make_hamiltonian()' function we get the qubit Hamiltonian. \n",
"We can create expectation values for tequila and optimize them in the same manner as with any other Hamiltonian. \n",
"The chemistry module offers some convenience in creating pre-defined circuits. One example is the so-called separable pair approximation (SPA) described in [arxiv:2105.03836](https://arxiv.org/abs/2105.03836). "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SPA Energy of H2/(2,4) is -1.05358\n"
]
}
],
"source": [
"H = h2_2_4.make_hamiltonian()\n",
"U = h2_2_4.make_ansatz(name=\"SPA\")\n",
"E = tq.ExpectationValue(H=H, U=U)\n",
"result = tq.minimize(E, silent=True)\n",
"print(\"SPA Energy of H2/({},{}) is {:+2.5f}\".format(mol.n_electrons, H.n_qubits, result.energy))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### VQE: Other Methods\n",
"\n",
"Other circuits can be initialized in the same manner and the expectatiovalues can be modified and combined like all other tequila expectation values (see the [main tutorials](https://github.com/tequilahub/tequila-tutorials) or [here](https://kottmanj.github.io/tequila-in-a-nutshell/#/)). \n",
"Here are some examples (all perform well since H2 on 4 qubits is quite easy)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"energy = -1.0535782052050917\n",
"variance = 1.1920928955078125e-07\n",
"energy = -1.0535782052050917\n",
"variance = 1.1920928955078125e-07\n"
]
}
],
"source": [
"U1 = mol.make_ansatz(name=\"UpCCGSD\")\n",
"U2 = tq.gates.X(0) + tq.gates.Ry(angle=\"a\",target=2) + tq.gates.CNOT(2,0) + tq.gates.CNOT(0,1) + tq.gates.CNOT(2,3)\n",
"for U in [U1,U2]:\n",
" E = tq.ExpectationValue(H=H, U=U)\n",
" result = tq.minimize(E, silent=True)\n",
" V = E**2 - tq.ExpectationValue(H=H**2, U=U)\n",
" # make it executable\n",
" V = tq.compile(V)\n",
" # evaluate with optimized variables\n",
" variance = V(result.variables)\n",
" print(\"energy = \", result.energy)\n",
" print(\"variance = \", variance)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compute More Orbitals\n",
"\n",
"We can of course also compute more orbitals with the keyword `n_pno` that defines the total number of PNOs that shall be computed. Let's compute 2 and 3 PNOs (so in total 3 and 4 orbitals leading to 6 and 8 qubits) and compare the VQE energies with an separable pair ansatz (SPA)."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals\n",
"output redirected to h2_pno_integrals.out logfile\n",
"finished after 21.300454139709473s\n",
"Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals\n",
"output redirected to h2_pno_integrals.out logfile\n",
"finished after 21.94425892829895s\n"
]
}
],
"source": [
"h2_2_6 = tq.Molecule(geometry=geometry, n_pno=2)\n",
"h2_2_8 = tq.Molecule(geometry=geometry, n_pno=3)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 4 qubits, energy = -1.05358, fci = -1.05358\n",
" 6 qubits, energy = -1.05683, fci = -1.05700\n",
" 8 qubits, energy = -1.05969, fci = -1.05988\n"
]
}
],
"source": [
"for mol in [h2_2_4, h2_2_6, h2_2_8]:\n",
" H = mol.make_hamiltonian()\n",
" U = mol.make_ansatz(name=\"SPA\")\n",
" E = tq.ExpectationValue(H=H, U=U)\n",
" result = tq.minimize(E, silent=True)\n",
" fci = mol.compute_energy(\"fci\") # needs pyscf installed\n",
" print(\"{:3} qubits, energy = {:2.5f}, fci = {:2.5f}\".format(H.n_qubits, result.energy, fci))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the SPA wavefunction produces identical energie than full-CI (FCI) which is equivalent to the exact groundstate in the given orbital basis. This is not further surprising as it was shown that SPA is exact for two-electron systems in [arxiv:2207.12421](https://arxiv.org/abs/2207.12421) provided that the orbitals that form the basis are in an optimal linear combination. In this case it looks like the PNOs are already optimal.\n",
"\n",
"Lets do the same computation with standard basis-sets of similar size. As tequila can not exploit PNO-structure to automatically create the SPA ansatz we need to provide the information which orbitals are assigned to which edge of the molecular graph (see [arxiv:2207.12421](https://arxiv.org/abs/2207.12421)) - for H2 this is trivial as there is just one edge/bond."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"converged SCF energy = -0.910873554594386\n",
" 4 qubits, energy = -0.99815, fci = -0.99815\n",
"converged SCF energy = -0.997497294328357\n",
" 8 qubits, energy = -1.04200, fci = -1.05435\n"
]
}
],
"source": [
"# we need pyscf installed for this cell to execute\n",
"for basis_set in [\"sto-3g\", \"6-31G\"]:\n",
" mol = tq.Molecule(geometry=geometry, basis_set=basis_set)\n",
" H = mol.make_hamiltonian()\n",
" U = mol.make_ansatz(name=\"SPA\", edges=[[i for i in range(mol.n_orbitals)],])\n",
" E = tq.ExpectationValue(H=H, U=U)\n",
" result = tq.minimize(E, silent=True)\n",
" fci = mol.compute_energy(\"fci\")\n",
" print(\"{:3} qubits, energy = {:2.5f}, fci = {:2.5f}\".format(H.n_qubits, result.energy, fci))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When comparing the SPA and FCI energies we see, that the orbitals of 6-31G (canonical Hartree-Fock) are not yet optimal (this was also observed in [arxiv:2207.12421](https://arxiv.org/abs/2207.12421)). We can optimize them with `tq.chemistry.optimize_orbitals` (see paper).\n",
"\n",
"Comparing the FCI energies of the MRA-PNO molecules above and the standard basis sets used here we see that the MRA-PNOs lead to lower energies for the same qubit sizes. Due to the variational principle, this means that we created a better basis with the system adapted approach."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Serious Example: BeH2\n",
"\n",
"Let's do a similar calculation for the BeH2 molecule. The default is again to compute one spatial orbital for each electron (or in other words: two spatial orbitals for each electron pair). In [arxiv:2207.12421](https://arxiv.org/abs/2207.12421) this was called a minimally correlated basis. \n",
"\n",
"Let's compute BeH2 close to it's equilibrium geometry and run an VQE with the SPA ansatz"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals\n",
"output redirected to beh2_pno_integrals.out logfile\n",
"finished after 45.55256175994873s\n",
" 4 electrons\n",
" 4 orbitals\n",
" 8 qubits\n"
]
}
],
"source": [
"beh2_geom = \"Be 0.0 0.0 0.0\\nH 0.0 0.0 {R}\\nH 0.0 0.0 -{R}\"\n",
"beh2_4_8_15 = tq.Molecule(geometry=beh2_geom.format(R=1.5))\n",
"H1 = beh2_4_8_15.make_hamiltonian()\n",
"\n",
"print(\"{:2} electrons\".format(beh2_4_8_15.n_electrons))\n",
"print(\"{:2} orbitals\".format(beh2_4_8_15.n_orbitals))\n",
"print(\"{:2} qubits\".format(H1.n_qubits))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A BeH2 molecule with 4 electrons and 8 qubits (spin-orbitals) was computed by default. But shouldn't it be 6 electrons? This is correct, the so called core-electrons (the electron pair 'sitting' close to the Be atom) is automatically frozen (the so-called 'frozen-core' approximation, standard in most quantum chemistry applications). We can deactivate this with tq.Molecule(...,frozen_core=False), but for most cases it won't lead to much (see [arxiv:2207.12421](https://arxiv.org/abs/2207.12421) for a detailed analysis using LiH). \n",
"\n",
"So let's compute the SPA-VQE. We will see that it is almost the same as the exact (FCI) energy (see https://arxiv.org/abs/2105.03836)."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"BeH2 with R=1.5\n",
"SPA error = -0.00170\n"
]
}
],
"source": [
"U = beh2_4_8_15.make_ansatz(name=\"SPA\")\n",
"H = beh2_4_8_15.make_hamiltonian()\n",
"E = tq.ExpectationValue(H=H, U=U)\n",
"result = tq.minimize(E, silent=True)\n",
"fci = beh2_4_8_15.compute_energy(\"fci\")\n",
"\n",
"print(\"BeH2 with R=1.5\")\n",
"print(\"SPA error = {:2.5f}\".format(fci-result.energy))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's do this again with a larger bond distance"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals\n",
"output redirected to beh2_pno_integrals.out logfile\n",
"finished after 82.53703570365906s\n",
"BeH2 with R=1.5\n",
"SPA error = -0.19778\n"
]
}
],
"source": [
"beh2_4_8_45 = tq.Molecule(geometry=beh2_geom.format(R=4.5))\n",
"H = beh2_4_8_45.make_hamiltonian()\n",
"U = beh2_4_8_45.make_ansatz(name=\"SPA\")\n",
"E = tq.ExpectationValue(H=H, U=U)\n",
"result = tq.minimize(E, silent=True)\n",
"fci = beh2_4_8_45.compute_energy(\"fci\") # needs pyscf\n",
"\n",
"print(\"BeH2 with R=1.5\")\n",
"print(\"SPA error = {:2.5f}\".format(fci-result.energy))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The SPA error is much larger (also shown in [arxiv:2105.03836](https://arxiv.org/abs/2105.03836)) but we can now try to re-optimize the orbitals to be optimal for an SPA ansatz in the given orbital basis. Note that the orbital basis is not changed but that we just form new linear combinations out of the existing MRA-PNO orbitals."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"opt-SPA error = -0.00029\n"
]
}
],
"source": [
"# needs pyscf\n",
"opt = tq.chemistry.optimize_orbitals(circuit=U, molecule=beh2_4_8_45, initial_guess=\"random\", silent=True)\n",
"print(\"opt-SPA error = {:2.5f}\".format(fci-opt.energy))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data I/O\n",
"\n",
"The datafiles created by the madness backend can be used in order to read in the orbital-data (in form of the corresponding one- and two-body integrals). This is useful to avoid recomputation of the orbitals which is costly in madness or to read in pre-computed orbitals when madness is not installed on the system. \n",
"\n",
"The data input/output is controlled over two keywords:\n",
"- `name`: the name of the molecule. If not set this is automatically deduced from the geometry\n",
"- `datadir`: the directory where the data should be stored or loaded from (needs tq version 1.8.1 or larger)\n",
"- `n_pno`: if not set a minimally correlated set of orbitals is computed, if set to `n_pno=\"read\"` data is read in\n",
"\n",
"The three files generated by madness are\n",
"- 'name_htensor.npy': The one-body integrals (kinetic energy + potential energy from the nuclear-potential)\n",
"- 'name_gtensor.npy': The two-body integrals (electronic repulsion integrals - or \"eri\") \n",
"in \"Mulliken\" (or \"Chemist\") notation (see the tequila paper [arxiv:2011.03057](https://arxiv.org/abs/2011.03057)). \n",
"- 'name_pnoinfo.txt': Information about the orbitals (which ones are HF which ones PNOs and which pairs are the PNOs assigned to)\n",
"- 'name_pno_integrals.out': File is not needed by tequila but contains the output of madness\n",
"\n",
"The gitub repo [github.com/kottmanj/moldata](https://github.com/kottmanj/moldata) contains for example a small collection of precomputed orbitals.\n",
"\n",
"Here is a small example"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals\n",
"output redirected to my_molecule_pno_integrals.out logfile\n",
"finished after 13.88163685798645s\n",
"contents of directory: my_data_dir\n",
"['my_molecule_pno_integrals.out', 'my_molecule_gtensor.npy', 'my_molecule_pnoinfo.txt', 'my_molecule_htensor.npy']\n",
"new molecule with 2 orbitals\n"
]
}
],
"source": [
"import tequila as tq\n",
"mol1 = tq.Molecule(name=\"my_molecule\", geometry=\"he 0.0 0.0 0.0\", datadir=\"my_data_dir\")\n",
"# see in the printout where the two integral tensors came from\n",
"# print(mol)\n",
"print(\"contents of directory: my_data_dir\")\n",
"import os\n",
"print(os.listdir(\"my_data_dir\"))\n",
"# now read it in again\n",
"mol2 = tq.Molecule(name=\"my_molecule\", geometry=\"he 0.0 0.0 0.0\", datadir=\"my_data_dir\", n_pno=\"read\")\n",
"print(\"new molecule with \", mol2.n_orbitals, \" orbitals\")\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## For Madness Experts\n",
"\n",
"If you know how the mandess input structure works you can modify the keywords like this (HF section `dft` and pno section). \n",
"Take for example a look into the 'he_pno_integrals.out' file to see all possible keywords. \n",
"Consider however, that not all keywords will have an effect on the specialized pno_integrals routine."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Starting madness calculation with executable: /home/jsk/anaconda3/envs/tq-1.8.1/bin/pno_integrals\n",
"output redirected to he_pno_integrals.out logfile\n",
"finished after 30.72634220123291s\n"
]
}
],
"source": [
"mol = tq.Molecule(geometry=\"he 0.0 0.0 0.0\", dft={\"k\":9, \"L\":25.0}, pno={\"maxrank\":4}, frozen_core=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## More Information\n",
"\n",
"One of `tequila`s primary aims is to simplify usage of many specialized algorithms and programs. In the following you will find a list of articles that describe some of the technology that is applied behind the scenes when you are using tequila with madness as chemistry backend:\n",
"\n",
"#### Technology automatically applied in the background within this tutorial\n",
"- [arxiv:2008.02819](https://arxiv.org/abs/2008.02819): initial article about basis-set-free VQEs\n",
"- [doi:10.1063/1.5141880](https://doi.org/10.1063/1.5141880): the MP2-PNO surrogate model in MRA representation\n",
"- [arxiv:1507.01888](https://arxiv.org/abs/1507.01888): `madness` overview\n",
"- [arxiv:2011.03057](https://arxiv.org/abs/2011.03057): `tequila` overview\n",
"- [arxiv:2011.05938](https://arxiv.org/abs/2011.05938): affordable gradients for UCC operations automatically available in tequila\n",
"- [arxiv:2105.03836](https://arxiv.org/abs/2105.03836): Separable Pair Ansatz with automatic orbital construction though MRA-PNOs from madness\n",
"- [arxiv:2207.12421](https://arxiv.org/abs/2207.12421): Graph-Based circuit design\n",
"\n",
"\n",
"#### Other Backends used in this tutorial\n",
"- [`OpenFermion`](https://github.com/quantumlib/OpenFermion): Fermion to qubit mappings\n",
"- [`qulacs`](https://github.com/qulacs/qulacs): Qulacs quantum circuit simulator\n",
"- [scipy](https://github.com/scipy/scipy): Optimizers\n",
"- [jax](https://github.com/google/jax)): Gradients of python functions\n",
"\n",
"#### Last stable run of this notebook\n",
"\n",
"Used anaconda3 with conda v4.12.0\n",
"```bash\n",
"conda create -n test_env python=3.8\n",
"conda activate test_env\n",
"\n",
"conda install madtequila -c kottmann\n",
"python -m pip install --upgrade pip\n",
"python -m pip install \"tequila-basic==1.8.1\"\n",
"python -m pip install \"qulacs==0.3.1\"\n",
"python -m pip install \"pyscf==2.0.1\"\n",
"```\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "tq-1.8.1",
"language": "python",
"name": "tq-1.8.1"
},
"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.8.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}