Tackling larger systems with fragmentation

In the basic VQE and extended VQE 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.

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

[1]:
from inquanto.geometries import GeometryMolecular
from inquanto.embeddings import DMETRHF
from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHF
from inquanto.extensions.pyscf import DMETRHFFragmentPySCFCCSD, get_fragment_orbital_masks, get_fragment_orbitals
from pytket.extensions.qiskit import AerStateBackend

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.

image0

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.

[2]:
# ################# #
# MOLECULE & DRIVER #
# ################# #

xyz = [['C', [0.000, 0.442, 0.000]],
            ['O', [-1.046, -0.467, 0.000]],
            ['O', [1.171, 0.120, 0.000]],
            ['H', [-0.389, 1.475, 0.000]],
            ['H', [-0.609, -1.355, 0.000]]]

geometry = GeometryMolecular(xyz)

basis = 'sto-3g'
charge = 0
driver = ChemistryDriverPySCFMolecularRHF(basis='sto-3g',
                                        geometry=geometry,
                                              charge=0,
                                              verbose=0)

hamiltonian_operator, space, rdm1 = driver.get_lowdin_system()

dmet = DMETRHF(hamiltonian_operator, rdm1)

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.

Finally, we initialize the DMET method with the hamiltonian_operator and the rdm1.

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.

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:

[3]:
get_fragment_orbitals(driver, [0], [2], [1, 4], [3])
[3]:
0 1 2 3 4
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

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.

[4]:
maskC, maskO, maskOH, maskH = get_fragment_orbital_masks(driver, [0], [2], [1, 4], [3])

fragments = [DMETRHFFragmentPySCFCCSD(dmet, maskC, name="C")]
fragments += [DMETRHFFragmentPySCFCCSD(dmet, maskO, name="O")]
fragments += [DMETRHFFragmentPySCFCCSD(dmet, maskOH, name="OH")]
fragments += [DMETRHFFragmentPySCFCCSD(dmet, maskH, name="H")]

result = dmet.run(fragments)
# STARTING CHEMICAL POTENTIAL 0.0
# STARTING CORR POTENTIAL PARAMETERS []

# FULL SCF ITERATION 0
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0
# FRAGMENT 0 - C: <H>=-186.31338432615382  EFR=-58.83801802368308  Q=0.020105631085693076
# FRAGMENT 1 - O: <H>=-186.31338532077524  EFR=-94.07923540833242  Q=-0.08915722473148513
# FRAGMENT 2 - OH: <H>=-186.29269508290201  EFR=-98.21603860323464  Q=-0.07076155306002363
# FRAGMENT 3 - H: <H>=-186.2321657513947  EFR=-3.805406920135852  Q=0.0029867861860859435
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0001
# FRAGMENT 0 - C: <H>=-186.31396619851347  EFR=-58.83898988857602  Q=0.0203144710244052
# FRAGMENT 1 - O: <H>=-186.31419656704261  EFR=-94.0797095185704  Q=-0.08903855659046123
# FRAGMENT 2 - OH: <H>=-186.29359289137574  EFR=-98.21647309737244  Q=-0.07064482022361496
# FRAGMENT 3 - H: <H>=-186.23226272837812  EFR=-3.805683950277546  Q=0.0030663357332569863
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.02612234664479768
# FRAGMENT 0 - C: <H>=-186.4660997003921  EFR=-59.09129439130545  Q=0.07462903595851689
# FRAGMENT 1 - O: <H>=-186.5257121220714  EFR=-94.20327799731163  Q=-0.058040637118688565
# FRAGMENT 2 - OH: <H>=-186.5276081971391  EFR=-98.32803927524893  Q=-0.04067365994329997
# FRAGMENT 3 - H: <H>=-186.25776880475584  EFR=-3.8776501906428535  Q=0.023767280256224987
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.02618319569991146
# FRAGMENT 0 - C: <H>=-186.46645711276096  EFR=-59.09188298091328  Q=0.0747559732776546
# FRAGMENT 1 - O: <H>=-186.5262076771365  EFR=-94.20356742638153  Q=-0.05796786917950669
# FRAGMENT 2 - OH: <H>=-186.5281562941116  EFR=-98.32829670367677  Q=-0.04060450873464205
# FRAGMENT 3 - H: <H>=-186.25782907787814  EFR=-3.8778181679788246  Q=0.02381568221012531
# CHEMICAL POTENTIAL 0.02618333425876141
# FINAL PARAMETERS: []
# FINAL CHEMICAL POTENTIAL: 0.02618333425876141
# FINAL ENERGY: -186.531161385121

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.

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.

[5]:
print("REFERENCE MP2  ENERGY: ", driver.run_mp2())
print("REFERENCE CCSD ENERGY: ", driver.run_ccsd())
REFERENCE MP2  ENERGY:  -186.3768719260956
REFERENCE CCSD ENERGY:  -186.4030088850797

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.

DMET allows the use of different levels of theory for each fragment. By using the examples, 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.