{ "cells": [ { "attachments": {}, "cell_type": "markdown", "id": "d29537f8-fd49-411a-a4bd-78d07c2f4044", "metadata": {}, "source": [ "Tackling larger systems with fragmentation\n", "===========================================\n", "\n", "In the [basic VQE](InQ_tut_vqe_1.ipynb) and [extended VQE](InQ_tut_vqe_2.ipynb) tutorials, we covered how to run a simple VQE calculation using InQuanto and some optimizations that can be performed. Here, we look at using Density Matrix Embedding Theory (DMET) to examine a larger system by fragmenting it. As an example, we consider HCOOH (formic acid). Without considering symmetry or active space reductions, this system would require 34 qubits to simulate using an STO-3G basis. This requires more resources than are available on current quantum computers, and will be extremely expensive to simulate on a classical device.\n", "\n", "DMET is a method of studying large molecules by partitioning the system into fragments containing a smaller number of atoms. Each fragment is treated independently in a bath corresponding to the molecular environment. Crucially, DMET allows for different fragments to be treated using different electronic structure methods. For example, we could imagine using VQE on the quantum computer to treat one particular fragment of interest. We focus here on a simplified implementation of DMET -- the so-called one-shot DMET. More examples are in the examples/embeddings folder. For discussion of the theory underpinning DMET, see [Knizia & Chan (2012)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.109.186404). As DMET relies heavily on performing classical electronic structure calculations in addition to any quantum computations, we need to import a driver and the fragment solvers from the `inquanto-pyscf` extension." ] }, { "cell_type": "code", "execution_count": null, "id": "dfc0b93e-14ab-4cbf-b411-a2f2b16814c2", "metadata": {}, "outputs": [], "source": [ "from inquanto.geometries import GeometryMolecular\n", "from inquanto.embeddings import DMETRHF\n", "from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHF\n", "from inquanto.extensions.pyscf import DMETRHFFragmentPySCFCCSD, get_fragment_orbital_masks, get_fragment_orbitals\n", "from pytket.extensions.qiskit import AerStateBackend" ] }, { "attachments": {}, "cell_type": "markdown", "id": "be2e7a06-a884-4e37-9b21-d349055fca0c", "metadata": {}, "source": [ "In order to use DMET to study our system, we must choose a scheme to split the molecule into fragments. In general, this is a task which requires chemical intuition and an awareness of the resource implications of the size of each fragment. As our goal here is to perform a simulation that runs quickly (and not to obtain highly accurate results), we choose a very fine fragmentation scheme with several small fragments.\n", "\n", "![](_images/tutorial_3_formic.png)\n", "\n", "The above figure shows the fragmentation scheme graphically. We see that other than the hydroxyl, each atom is in its own fragment. This ensures that the maximum number of qubits required to simulate an individual fragment would be 12 (the hydroxyl fragment), a reasonable number of qubits to simulate on a large classical computer." ] }, { "cell_type": "code", "execution_count": null, "id": "ba6b44d0-2ae6-4eb1-9396-283362f870a3", "metadata": {}, "outputs": [], "source": [ "# ################# #\n", "# MOLECULE & DRIVER #\n", "# ################# #\n", "\n", "xyz = [['C', [0.000, 0.442, 0.000]],\n", " ['O', [-1.046, -0.467, 0.000]],\n", " ['O', [1.171, 0.120, 0.000]],\n", " ['H', [-0.389, 1.475, 0.000]],\n", " ['H', [-0.609, -1.355, 0.000]]]\n", "\n", "geometry = GeometryMolecular(xyz)\n", "\n", "basis = 'sto-3g'\n", "charge = 0\n", "driver = ChemistryDriverPySCFMolecularRHF(basis='sto-3g',\n", " geometry=geometry,\n", " charge=0,\n", " verbose=0)\n", "\n", "hamiltonian_operator, space, rdm1 = driver.get_lowdin_system()\n", "\n", "dmet = DMETRHF(hamiltonian_operator, rdm1)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "0a40b1ef", "metadata": {}, "source": [ "As before, we first initialise the driver. The initialization of the driver is much the same as in standard VQE with regards to the molecule geometry, basis and charge specification. However, because of the spatial fragmentation, DMET requires the localisation of the molecular orbitals. Therefore, when we compute the Hamiltonian operator, instead of `get_system()` we call the `get_lowdin_system()` method. This will perform an RHF simulation of the molecule, and will return the mean-field 1e-RDM as `rdm1`, the `space` and the `hamiltonian_operator` in the Löwdin basis, that is computed by Löwdin symmetric orthogonalization of the atomic orbitals.\n", "\n", "Finally, we initialize the DMET method with the `hamiltonian_operator` and the `rdm1`.\n", "\n", "At this stage, we have not specified any particular fragmentation scheme. Fragments in InQuanto are associated with the level of electronic structure theory that will be used to simulate them. As a first test, we try specifying that each fragment should be treated with classical CCSD." ] }, { "attachments": {}, "cell_type": "markdown", "id": "c34f707e", "metadata": {}, "source": [ "In order to specify a fragment, we need to determine the corresponding Löwdin orbitals. The easiest way to specify the fragments is by atoms. Based on the indices in the `xyz` table, we can make four sets of atom indices, and using the `get_fragment_orbitals()` utility function we can tabulate the orbitals and the four orbital fragment masks that select the Löwdin orbitals corresponding to the fragments:" ] }, { "cell_type": "code", "execution_count": null, "id": "2b4a747d", "metadata": {}, "outputs": [ { "data": { "application/vnd.microsoft.datawrangler.viewer.v0+json": { "columns": [ { "name": "index", "rawType": "int64", "type": "integer" }, { "name": "0", "rawType": "object", "type": "string" }, { "name": "1", "rawType": "bool", "type": "boolean" }, { "name": "2", "rawType": "bool", "type": "boolean" }, { "name": "3", "rawType": "bool", "type": "boolean" }, { "name": "4", "rawType": "bool", "type": "boolean" } ], "ref": "e01d9970-b0fb-4ebe-894a-d84591324186", "rows": [ [ "0", "0 C 1s ", "True", "False", "False", "False" ], [ "1", "0 C 2s ", "True", "False", "False", "False" ], [ "2", "0 C 2px ", "True", "False", "False", "False" ], [ "3", "0 C 2py ", "True", "False", "False", "False" ], [ "4", "0 C 2pz ", "True", "False", "False", "False" ], [ "5", "1 O 1s ", "False", "False", "True", "False" ], [ "6", "1 O 2s ", "False", "False", "True", "False" ], [ "7", "1 O 2px ", "False", "False", "True", "False" ], [ "8", "1 O 2py ", "False", "False", "True", "False" ], [ "9", "1 O 2pz ", "False", "False", "True", "False" ], [ "10", "2 O 1s ", "False", "True", "False", "False" ], [ "11", "2 O 2s ", "False", "True", "False", "False" ], [ "12", "2 O 2px ", "False", "True", "False", "False" ], [ "13", "2 O 2py ", "False", "True", "False", "False" ], [ "14", "2 O 2pz ", "False", "True", "False", "False" ], [ "15", "3 H 1s ", "False", "False", "False", "True" ], [ "16", "4 H 1s ", "False", "False", "True", "False" ] ], "shape": { "columns": 5, "rows": 17 } }, "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01234
00 C 1sTrueFalseFalseFalse
10 C 2sTrueFalseFalseFalse
20 C 2pxTrueFalseFalseFalse
30 C 2pyTrueFalseFalseFalse
40 C 2pzTrueFalseFalseFalse
51 O 1sFalseFalseTrueFalse
61 O 2sFalseFalseTrueFalse
71 O 2pxFalseFalseTrueFalse
81 O 2pyFalseFalseTrueFalse
91 O 2pzFalseFalseTrueFalse
102 O 1sFalseTrueFalseFalse
112 O 2sFalseTrueFalseFalse
122 O 2pxFalseTrueFalseFalse
132 O 2pyFalseTrueFalseFalse
142 O 2pzFalseTrueFalseFalse
153 H 1sFalseFalseFalseTrue
164 H 1sFalseFalseTrueFalse
\n", "
" ], "text/plain": [ " 0 1 2 3 4\n", "0 0 C 1s True False False False\n", "1 0 C 2s True False False False\n", "2 0 C 2px True False False False\n", "3 0 C 2py True False False False\n", "4 0 C 2pz True False False False\n", "5 1 O 1s False False True False\n", "6 1 O 2s False False True False\n", "7 1 O 2px False False True False\n", "8 1 O 2py False False True False\n", "9 1 O 2pz False False True False\n", "10 2 O 1s False True False False\n", "11 2 O 2s False True False False\n", "12 2 O 2px False True False False\n", "13 2 O 2py False True False False\n", "14 2 O 2pz False True False False\n", "15 3 H 1s False False False True\n", "16 4 H 1s False False True False" ] }, "execution_count": null, "metadata": {}, "output_type": "execute_result" } ], "source": [ "get_fragment_orbitals(driver, [0], [2], [1, 4], [3])" ] }, { "cell_type": "markdown", "id": "2bcf7d26", "metadata": {}, "source": [ "With another utility function `get_fragment_orbital_masks(...)` we can obtain the orbital fragment masks as arrays. Once we have the orbital masks we can complete the DMET simulation. " ] }, { "cell_type": "code", "execution_count": null, "id": "a5cb88b2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# STARTING CHEMICAL POTENTIAL 0.0\n", "# STARTING CORR POTENTIAL PARAMETERS []\n", "\n", "# FULL SCF ITERATION 0\n", "# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0\n", "# FRAGMENT 0 - C: =-186.3133843261541 EFR=-58.838018023683354 Q=0.020105631085782782 \n", "# FRAGMENT 1 - O: =-186.31338532077373 EFR=-94.07923540836087 Q=-0.08915722472854348 \n", "# FRAGMENT 2 - OH: =-186.29269508290224 EFR=-98.21603860323472 Q=-0.07076155306003606 \n", "# FRAGMENT 3 - H: =-186.23216575139492 EFR=-3.8054069201358613 Q=0.0029867861860670697 \n", "# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0001\n", "# FRAGMENT 0 - C: =-186.31396619851375 EFR=-58.83898988857627 Q=0.020314471024494907 \n", "# FRAGMENT 1 - O: =-186.3141965670411 EFR=-94.07970951859826 Q=-0.08903855658758708 \n", "# FRAGMENT 2 - OH: =-186.29359289137597 EFR=-98.21647309737251 Q=-0.07064482022363094 \n", "# FRAGMENT 3 - H: =-186.23226272837834 EFR=-3.805683950277591 Q=0.0030663357332485486 \n", "# NEWTON ITERATION - CHEMICAL POTENTIAL 0.02612234664724773\n", "# FRAGMENT 0 - C: =-186.4660997004068 EFR=-59.091294391329406 Q=0.07462903596371717 \n", "# FRAGMENT 1 - O: =-186.52571212209148 EFR=-94.20327799732226 Q=-0.05804063711589613 \n", "# FRAGMENT 2 - OH: =-186.52760819716138 EFR=-98.32803927525933 Q=-0.04067365994052352 \n", "# FRAGMENT 3 - H: =-186.25776880475843 EFR=-3.8776501906496676 Q=0.0237672802581691 \n", "# NEWTON ITERATION - CHEMICAL POTENTIAL 0.026183195699930065\n", "# FRAGMENT 0 - C: =-186.46645711276136 EFR=-59.091882980913724 Q=0.0747559732777825 \n", "# FRAGMENT 1 - O: =-186.5262076771368 EFR=-94.20356742638096 Q=-0.05796786917957952 \n", "# FRAGMENT 2 - OH: =-186.52815629411197 EFR=-98.32829670367698 Q=-0.0406045087346385 \n", "# FRAGMENT 3 - H: =-186.2578290778784 EFR=-3.8778181679791044 Q=0.02381568221018371 \n", "# NEWTON ITERATION - CHEMICAL POTENTIAL 0.026183334258757535\n", "# FRAGMENT 0 - C: =-186.46645792663077 EFR=-59.09188432117837 Q=0.07475626232523247 \n", "# FRAGMENT 1 - O: =-186.52620880556586 EFR=-94.20356808543914 Q=-0.057967703478901456 \n", "# FRAGMENT 2 - OH: =-186.52815754218213 EFR=-98.32829728986279 Q=-0.04060435127064288 \n", "# FRAGMENT 3 - H: =-186.2578292151292 EFR=-3.877818550467012 Q=0.023815792422891557 \n", "# CHEMICAL POTENTIAL 0.02618333425902995\n", "# FINAL PARAMETERS: []\n", "# FINAL CHEMICAL POTENTIAL: 0.02618333425902995\n", "# FINAL ENERGY: -186.53116435311793\n" ] } ], "source": [ "maskC, maskO, maskOH, maskH = get_fragment_orbital_masks(driver, [0], [2], [1, 4], [3])\n", "\n", "fragments = [DMETRHFFragmentPySCFCCSD(dmet, maskC, name=\"C\")]\n", "fragments += [DMETRHFFragmentPySCFCCSD(dmet, maskO, name=\"O\")]\n", "fragments += [DMETRHFFragmentPySCFCCSD(dmet, maskOH, name=\"OH\")]\n", "fragments += [DMETRHFFragmentPySCFCCSD(dmet, maskH, name=\"H\")]\n", "\n", "result = dmet.run(fragments)" ] }, { "attachments": {}, "cell_type": "markdown", "id": "9e4306ad", "metadata": {}, "source": [ "In the first block, we have specified our fragments as a list. We use the `inquanto.extensions.pyscf.DMETRHFFragmentPySCFCCSD` class as we want to look at each fragment using classical coupled cluster. Each fragment takes the `dmet` as an argument in addition to an arbitrary string giving the fragment name. It also requires the masks for the orbitals within the fragment to be specified. These are given as an array of booleans, marking the index of orbitals with `True` if it is in the fragment and `False` if it is outside it.\n", "\n", "The `dmet.run()` method is then invoked passing the fragments as a list. During the execution we can observe details about the calculation. The FINAL ENERGY line in the end gives us the final ground state energy of the system calculated by DMET." ] }, { "cell_type": "code", "execution_count": null, "id": "ff55637a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "REFERENCE MP2 ENERGY: -186.3768719260956\n", "REFERENCE CCSD ENERGY: -186.40300888507952\n" ] } ], "source": [ "print(\"REFERENCE MP2 ENERGY: \", driver.run_mp2())\n", "print(\"REFERENCE CCSD ENERGY: \", driver.run_ccsd())" ] }, { "attachments": {}, "cell_type": "markdown", "id": "ecf371c9", "metadata": {}, "source": [ "As a reference we computed MP2 and CCSD energies. We can see that in this instance DMET obtains about 0.13 Ha more correlation energy than non-DMET classical CCSD. Although defeating the point of fragmenting the system to reduce resource requirements, benchmarking a DMET calculation with the same level of theory for each fragment against a non-DMET calculation is a good way to estimate the error incurred with the fragmentation scheme. Note that DMET is non-variational and thus can yield energies lower than the exact (i.e. FCI-level) energy.\n", "\n", "DMET allows the use of different levels of theory for each fragment. By using the [examples](https://docs.quantinuum.com/inquanto/tutorials/examples_overview.html), we encourage the reader to modify this notebook such that some of the fragments (for example the lone hydrogen atom) are using the VQE method with a state vector simulator. " ] } ], "metadata": {}, "nbformat": 4, "nbformat_minor": 5 }