# InQuanto-PySCF¶

The `inquanto-pyscf`

extension provides an interface to
PySCF, a classical computational chemistry package.
This extension allows a user to input chemistry information to be processed by using classical computational methods, and construct the data to be mapped
onto quantum computers.

## Basic usage¶

The main feature of this extension is to construct the fermionic second-quantized Hamiltonian in the molecular orbital basis,
along with the fermionic Fock space and state by running a classical computational mean-field (Hartree-Fock) calculation. An
`inquanto-pyscf`

driver may be used to set up the calculation, and the
`get_system()`

method generates all of these ingredients, shown
below for an H_{2} molecule at 0.75 Å separation, with the STO-3G atomic basis set:

```
"""Minimal basis H2."""
from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHF
# Initialize the chemistry driver.
driver = ChemistryDriverPySCFMolecularRHF(
geometry="H 0 0 0; H 0 0 0.75",
basis="sto3g",
)
# Get the data to be mapped onto quantum computers.
hamiltonian, fock_space, fock_state = driver.get_system()
print("Fock state:")
fock_space.print_state(fock_state)
print("Hamiltonian:")
print(hamiltonian.df())
```

```
Fock state:
0 0a : 1
1 0b : 1
2 1a : 0
3 1b : 0
Hamiltonian:
Coefficients Terms
0 0.705570
1 -1.247285 F0^ F0
2 0.336424 F1^ F0^ F0 F1
3 0.336424 F1^ F0^ F0 F1
4 0.090886 F1^ F0^ F2 F3
5 0.090886 F1^ F0^ F2 F3
6 -1.247285 F1^ F1
7 -0.090886 F2^ F0^ F0 F2
8 -0.090886 F2^ F0^ F0 F2
9 0.330989 F2^ F0^ F0 F2
10 0.330989 F2^ F0^ F0 F2
11 -0.090886 F2^ F1^ F0 F3
12 -0.090886 F2^ F1^ F0 F3
13 0.330989 F2^ F1^ F1 F2
14 0.330989 F2^ F1^ F1 F2
15 -0.481273 F2^ F2
16 0.330989 F3^ F0^ F0 F3
17 0.330989 F3^ F0^ F0 F3
18 -0.090886 F3^ F0^ F1 F2
19 -0.090886 F3^ F0^ F1 F2
20 -0.090886 F3^ F1^ F1 F3
21 -0.090886 F3^ F1^ F1 F3
22 0.330989 F3^ F1^ F1 F3
23 0.330989 F3^ F1^ F1 F3
24 0.090886 F3^ F2^ F0 F1
25 0.090886 F3^ F2^ F0 F1
26 0.347908 F3^ F2^ F2 F3
27 0.347908 F3^ F2^ F2 F3
28 -0.481273 F3^ F3
```

With the Hamiltonian, Fock space, and Hartree-Fock state, one may construct the quantum problem as usual (see here to get started).

Above, the `hamiltonian`

is an integral operator of type `ChemistryRestrictedIntegralOperator`

.
For molecular systems, PySCF drivers may produce more specialized integral operators.
For example, with the `symmetry`

option we generate a
`ChemistryRestrictedIntegralOperatorCompact`

object, which stores symmetry-compacted
representations of the two-body integral tensor for reducing memory requirements:

```
compact_hamiltonian, space, state = driver.get_system(symmetry=8) # Maximum symmetry reduction for restricted spins
```

Or, with the `get_system_ao()`

method, we get a
`PySCFChemistryRestrictedIntegralOperator`

, which wraps a PySCF `SCF`

object:

```
pyscf_hamiltonian, space, state = driver.get_system_ao()
```

The latter case is similarly useful for reducing classical memory requirements, and is also an important component in Fragment Molecular Orbital (FMO) calculations (see the embedding sections below).

## Generating a driver from a PySCF object¶

Experienced PySCF users may want to start their calculations with a PySCF mean-field object directly, before
preparing the quantum calculation. The `from_mf()`

method
initializes an `inquanto-pyscf`

driver in this way, from which we can prepare the
hamiltonian, Fock space, and Fock state as usual:

```
"""Generating a driver from a pyscf mean-field object."""
from pyscf import gto, scf
# PySCF calculation.
mol = gto.Mole(
atom='H 0 0 0; H 0 0 0.75',
basis='sto3g',
verbose=0,
).build()
mf = scf.RHF(mol)
mf.kernel()
# Initialize InQuanto driver
driver = ChemistryDriverPySCFMolecularRHF.from_mf(mf)
hamiltonian, fock_space, fock_state = driver.get_system()
```

Warning

Although this interface is useful in some cases, it is recommended to use the standard driver interface where possible. This advanced feature may cause a runtime error, as it is difficult to guarantee that the wrapper is fully consistent with all PySCF options.

## Periodic systems¶

Extended solid-state simulations at any level of dimensionality reduction can be performed by taking advantage of
PySCF’s periodic boundary conditions (PBC) capabilities. The `inquanto-pyscf`

extension
has two types of drivers for periodic calculations: Gamma point drivers, for \(k=0\)
calculations, and momentum space drivers, where a \(k\)-point grid is specified.

Below, we show a Gamma point example with the `ChemistryDriverPySCFGammaROHF`

driver
for a Pd (111) surface with a 2x2x1 slab and 20 Å of vacuum:

```
from pyscf.pbc import gto as pgto
from inquanto.extensions.pyscf import ChemistryDriverPySCFGammaRHF
cell_pd221=[
[5.50129075763134, 0.0, 0.0],
[2.75064537881567, 4.764257549713281, 0.0],
[0.0, 0.0, 20.0]
]
geometry_pd221=[
['Pd', [ 0., 0., 10.]],
['Pd', [ 2.75064538, 0. , 10. ]],
['Pd', [ 1.37532269, 2.38212877, 10. ]],
['Pd', [ 4.12596807, 2.38212877, 10. ]]
]
driver_pd221 = ChemistryDriverPySCFGammaRHF(
basis="lanl2dz",
ecp="lanl2dz",
geometry=geometry_pd221,
charge=0,
cell=cell_pd221,
dimension=2, # SLAB 2D PBC system
df="GDF",
output=None,
verbose=1,
exp_to_discard = 0.1
)
```

where we have also used Gaussian density fitting (GDF) to reduce the computational cost of operations with the
two-body electronic integrals (see here). From this point, one
may use the `get_system()`

method as usual to
build the fermionic hamiltonian and construct the quantum problem.

Calculations with a non-trivial \(k\)-point grid proceed similarly. Below, we consider an H2 chain (1D system) with PBC and a [4,1,1] \(k\)-point grid:

```
from inquanto.extensions.pyscf import ChemistryDriverPySCFMomentumRHF
cell_h2 = pgto.Cell()
cell_h2.atom=[['H', [0., 5., 5.]], ['H', [0.75, 5. , 5. ]]]
cell_h2.a=[(1.875,0.,0.), (0.,10.0,0.), (0.,0.,10.0)]
driver_h2 = ChemistryDriverPySCFMomentumRHF(
basis="sto-3g",
geometry=cell_h2.atom,
charge=0,
cell=cell_h2.a,
nks=[4, 1, 1],
dimension=1, # 1D linear chain PBC system
precision=1e-15,
df="GDF",
output=None,
verbose=1,
exp_to_discard = 0.1
)
print("Hartree-Fock energy: ", driver_h2.run_hf())
```

```
Hartree-Fock energy: 14080944.281806948
```

In this case, `get_system()`

generates a special
Fock space class: `FermionSpaceBrillouin`

, which manages the \(k\)-point quantum number
in addition to orbital and spin indices. Otherwise, construction of the quantum problem proceeds as normal.

## Classical post-HF calculations¶

Classical post-HF energies are useful to compare against quantum computational results for quick benchmarking. Each chemistry driver has an interface to the post-HF calculators provided by PySCF:

```
driver = ChemistryDriverPySCFMolecularRHF(
geometry="H 0 0 0; H 0 0 0.75",
basis="631g",
)
driver.get_system()
# Classical Post-HF energies for benchmarking.
print('HF\t', driver.mf_energy)
print('MP2\t', driver.run_mp2())
print('CCSD\t', driver.run_ccsd())
print('CASCI\t', driver.run_casci())
```

```
HF -1.1265450345356904
MP2 -1.1440347834365332
CCSD -1.1516885473648508
CASCI -1.151688547516609
```

## InQuanto-PySCF driver classes¶

A variety of drivers are available for different systems and applications, many of which are covered in more detail in the following sections. A full list of all drivers is given below:

For molecular systems:

`inquanto.extensions.pyscf.ChemistryDriverPySCFMolecularRHFQMMMCOSMO`

`inquanto.extensions.pyscf.ChemistryDriverPySCFMolecularROHFQMMMCOSMO`

`inquanto.extensions.pyscf.ChemistryDriverPySCFMolecularUHFQMMMCOSMO`

`inquanto.extensions.pyscf.ChemistryDriverPySCFEmbeddingROHF_UHF`

For periodic systems:

## Active space specification and AVAS¶

It is standard practice to select an orbital active space to reduce the resource requirements of quantum chemistry
calculations. The active space may be specified manually using the `frozen`

argument in an
`inquanto-pyscf`

driver. Below we consider the example of an H:sub:`2`O
molecule with the 6-31G basis set. The full Fock space is given by:

```
"""H2O with 6-31G basis set."""
from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHF
h2o_geom = """
O 0.00000 0.00000 0.11779;
H 0.00000 0.75545 -0.47116;
H 0.00000 -0.75545 -0.47116;"""
# Initialize the chemistry driver.
driver = ChemistryDriverPySCFMolecularRHF(
geometry=h2o_geom,
basis="631g",
)
hamiltonian, fock_space, fock_state = driver.get_system()
print("Fock state:")
fock_space.print_state(fock_state)
```

```
Fock state:
0 0a : 1
1 0b : 1
2 1a : 1
3 1b : 1
4 2a : 1
5 2b : 1
6 3a : 1
7 3b : 1
8 4a : 1
9 4b : 1
10 5a : 0
11 5b : 0
12 6a : 0
13 6b : 0
14 7a : 0
15 7b : 0
16 8a : 0
17 8b : 0
18 9a : 0
19 9b : 0
20 10a : 0
21 10b : 0
22 11a : 0
23 11b : 0
24 12a : 0
25 12b : 0
```

We may freeze molecular orbitals by passing a list of frozen orbital indices:

```
"""H2O with 6-31G basis set with 4-orbital 4-electron active space."""
driver = ChemistryDriverPySCFMolecularRHF(
geometry=h2o_geom,
basis="631g",
frozen=[0, 1, 2, 7, 8, 9, 10, 11, 12],
)
hamiltonian, fock_space, fock_state = driver.get_system()
print("Fock state (reduced by setting the active space):")
fock_space.print_state(fock_state)
```

```
Fock state (reduced by setting the active space):
0 0a : 1
1 0b : 1
2 1a : 1
3 1b : 1
4 2a : 0
5 2b : 0
6 3a : 0
7 3b : 0
```

or, equivalently, by using the `inquanto.extensions.pyscf.FromActiveSpace`

callable class, where we specify the number
of active spatial orbitals and electrons:

```
from inquanto.extensions.pyscf import FromActiveSpace
driver = ChemistryDriverPySCFMolecularRHF(
geometry=h2o_geom,
basis="631g",
frozen=FromActiveSpace(ncas=4, nelecas=4),
)
hamiltonian, fock_space, fock_state = driver.get_system()
print("Fock state (reduced by setting the active space):")
fock_space.print_state(fock_state)
```

```
Fock state (reduced by setting the active space):
0 0a : 1
1 0b : 1
2 1a : 1
3 1b : 1
4 2a : 0
5 2b : 0
6 3a : 0
7 3b : 0
```

InQuanto also supports the use of Atomic Valence Active Space (AVAS), a technique for selecting an active space for a localized domain of a molecule. AVAS constructs the active space by projecting the target molecular orbitals onto a set of atomic orbitals.

To use AVAS, instantiate the `inquanto.extensions.pyscf.AVAS`

class with the atomic orbital labels, and build
the PySCF driver as follows:

```
from inquanto.extensions.pyscf import AVAS
avas= AVAS(aolabels=['Li 2s', 'H 1s'], threshold=0.8, threshold_vir=0.8, verbose=5)
driver = ChemistryDriverPySCFMolecularRHF(
geometry='Li 0 0 0; H 0 0 1.75',
basis='631g',
transf=avas,
frozen=avas.frozenf
)
hamiltonian, fock_space, fock_state = driver.get_system()
print("Fock state (reduced by setting the active space with AVAS):")
fock_space.print_state(fock_state)
```

```
```

```
******** AVAS flags ********
```

```
aolabels = ['Li 2s', 'H 1s']
```

```
aolabels_vir = ['Li 2s', 'H 1s']
```

```
frozen = []
```

```
minao = minao
```

```
threshold = 0.8
```

```
threshold (virtual) = 0.8
```

```
with_iao = False
```

```
force_halves_active = False
```

```
canonicalize = True
```

```
** AVAS **
```

```
Total number of electrons: 4.0
```

```
Total number of HF MOs is equal to 11
```

```
Number of occupied HF MOs is equal to 2
```

```
reference AO indices for minao ['Li 2s', 'H 1s']: [1 2]
```

```
(full basis) reference AO indices for 631g ['Li 2s', 'H 1s']: [1 9]
```

```
Threshold 0.8, threshold_vir 0.8
```

```
Active from occupied = 1 , eig [0.96867354]
```

```
Inactive from occupied = 1
```

```
Active from unoccupied = 1 , eig [0.98544091]
```

```
Inactive from unoccupied = 8
```

```
Number of active orbitals 2
```

```
# of alpha electrons 1
```

```
# of beta electrons 1
```

```
Fock state (reduced by setting the active space with AVAS):
0 0a : 1
1 0b : 1
2 1a : 0
3 1b : 0
```

where the `transf`

argument specifies the AVAS orbital transformation, and the `AVAS.frozenf`

attribute
returns the list of frozen orbitals.

## Energy correction with NEVPT2 and AC0¶

InQuanto introduces an innovative approach merging quantum and classical techniques for strongly contracted N-electron Valence State 2nd-order Perturbation Theory (SC-NEVPT2), replacing CASCI’s static correlation with quantum computer-simulated n-particle Reduced Density Matrices (n-RDMs) measurements [61]. These measurements are then employed in a classical SC-NEVPT2 calculation to approximate the remaining dynamic electron correlation effects, marking a pioneering hybrid quantum-classical calculation for multi-reference perturbation theory, utilizing a quantum computer for the quantum component.

CASSCF calculations are initially performed using the PySCF package to obtain Pre-Density Matrices (PDMs) [62] and RDMs. The `inquanto-pyscf`

driver
sets up the calculation, and `from_mf()`

initializes the driver from a PySCF mean-field object. The example
involves an Li_{2} molecule with 6-31g atomic basis set. HF orbital coefficients are replaced with those from CASSCF to enhance the orbital
description and system representation. The `get_nevpt2_correction`

function is employed to calculate the strongly contracted NEVPT2 correction to
the energy using the provided density matrices. The `make_dm1234`

function in PySCF calculates the PDMs, which do not directly correspond to
the 1-, 2-, 3-, and 4-particle density matrices.

```
from inquanto.extensions.pyscf import FromActiveSpace,ChemistryDriverPySCFMolecularRHF
from pyscf import gto, scf, fci, mcscf
geometry = [["Li", [0.0, 0.0, 0.0]], ["Li", [2.63, 0.0, 0.0]]]
ncas, nelecas = 4, 2
mf = gto.M(atom=geometry, basis="6-31g").apply(scf.RHF).run()
mc = mcscf.CASSCF(mf, ncas, nelecas )
mc.run()
mf.mo_coeff=mc.mo_coeff
(pdm1, pdm2, pdm3, pdm4) = fci.rdm.make_dm1234(
"FCI4pdm_kern_sf", mc.ci, mc.ci, ncas, nelecas
)
driver_casscf = ChemistryDriverPySCFMolecularRHF.from_mf(
mf, frozen=FromActiveSpace(ncas, nelecas)(mf)
)
nevpt2_energy = driver_casscf.get_nevpt2_correction(rdms=(pdm1, pdm2, pdm3, pdm4))
```

```
converged SCF energy = -14.8652654987828
```

```
CASSCF energy = -14.8887816132809
```

```
CASCI E = -14.8887816132809 E(CI) = -0.619242395586962 S^2 = 0.0000000
```

```
Sr (-1)', E = 0.00000000000062
```

```
Si (+1)', E = -0.00001592231055
```

```
Sijrs (0) , E = -0.00022778405930
```

```
Sijr (+1) , E = -0.00007024183829
```

```
Srsi (-1) , E = -0.00023752541083
```

```
Srs (-2) , E = -0.00127421569351
```

```
Sij (+2) , E = -0.00002668012276
```

```
Sir (0)' , E = -0.00013325630380
```

```
Nevpt2 Energy = -0.001985625738416
```

The first-order expansion of the Adiabatic Connection for Complete Active Space wave functions (AC0-CAS) method [63], incorporating 1- and 2-particle reduced density matrices, utilizes the AC method and CAS approach within quantum chemistry.
By employing the same driver, one can calculate the AC0 energy correction, utilizing the `get_ac0_correction`

function where
we specify the actual 1- and 2-particle density matrices converted from PDMs using the `reorder_dm12`

function.

```
rdm1, rdm2 = fci.rdm.reorder_dm12(pdm1, pdm2)
ac0_energy = driver_casscf.get_ac0_correction((rdm1, rdm2))
```

Note

To perform the AC0 calculation, it is necessary to install pyscf-AC0. See here: here.

An alternative approach would be to perform VQE calculations to obtain the PDMs and RDMs. A detailed tutorial including both approaches can be found here.

## Hamiltonians for Embedding methods¶

PySCF drivers provide an interface to various orbital localization schemes and a wide range of high-level chemistry methods. These can be utilized in a fragment solver for Density Matrix Embedding Theory (DMET) and Fragment Molecular Orbital (FMO) embedding methods.

Prior to an embedding algorithm one needs to generate a Hamiltonian
with a localized and orthonormal basis or atomic orbitals, in which spatial fragments can
be specified. To this end, PySCF drivers can generate the Löwdin representation of a system
with the `get_lowdin_system`

method:

```
driver = ChemistryDriverPySCFMolecularRHF(basis="sto-3g", geometry=g, charge=0)
hamiltonian_operator_lowdin, space, rdm1_hf_lowdin = driver.get_lowdin_system()
```

which performs a full Hartree-Fock calculation, and returns the 1-RDM in the localized basis. This representation is required for DMET calculations.

Additionally, for FMO calculations, one may use the `get_system_ao`

method to generate a `PySCFChemistryRestrictedIntegralOperator`

which,
internally, has access to the atomic orbitals for generating fragment Hamiltonians:

```
driver = ChemistryDriverPySCFMolecularRHF(basis="sto-3g", geometry=g, charge=0)
hamiltonian_operator_pyscf, space, rdm1_hf_pyscf = driver.get_system_ao(run_hf=False)
```

where the `run_hf=False`

option prevents a full Hartree-Fock calculation over the entire system.

The general API of both DMET and FMO has two types of classes that are necessary to perform a calculation:
one that is responsible for the total system with the embedding method itself,
and another that is responsible to compute the high-level solution of a particular fragment. For example,
as it is discussed in the Density Matrix Embedding Theory section, there is an
`ImpurityDMETROHF`

class driving the embedding method, and there are
fragment solver classes such as `ImpurityDMETROHFFragmentED`

. Similarly
for FMO we have the `inquanto.extensions.pyscf.fmo.FMO`

and
`inquanto.extensions.pyscf.fmo.FMOFragmentPySCFCCSD`

, for example.

## DMET with PySCF fragment solvers¶

Here we outline the use of PySCF calculators as fragment solvers for both Impurity DMET and full DMET with InQuanto. The reader is referred to the main DMET section for an introduction.

### Impurity DMET¶

InQuanto’s core packages provide a simple exact diagonalization solver and a Hartree-Fock fragment solver for impurity
DMET. `inquanto-pyscf`

provides additional fragment solvers:

Provided we have the Hamiltonian and the density matrix in a localized basis we can run the
single fragment Impurity DMET method. Below we consider a Hamiltonian computed with the Löwdin transformation, with
6 spatial orbitals. The `ImpurityDMETROHF`

class is initialized as follows:

```
from inquanto.embeddings import ImpurityDMETROHF
dmet = ImpurityDMETROHF(
hamiltonian_operator_lowdin, rdm1_hf_lowdin
)
```

The class name of `ImpurityDMETROHF`

denotes that this method assumes the Hamiltonian
operator and the 1-RDM are spin restricted. After this initialization, a fragment and a corresponding fragment
solver are specified in the following way:

```
from numpy import array
from inquanto.extensions.pyscf import (
ImpurityDMETROHFFragmentPySCFFCI,
ImpurityDMETROHFFragmentPySCFMP2,
ImpurityDMETROHFFragmentPySCFCCSD,
ImpurityDMETROHFFragmentPySCFROHF,
ImpurityDMETROHFFragmentPySCFActive
)
mask = array([True, True, False, False, False, False])
fragment = ImpurityDMETROHFFragmentPySCFFCI(dmet, mask)
```

`ImpurityDMETROHFFragmentPySCFFCI`

defines the high-level solver for the
fragment, which in this case is an FCI method. The fragment solver is initialized with the `dmet`

object and the
fragment mask `mask`

. The `mask`

is a boolean array with the length of the number of localized spatial orbitals.
The `True`

values in the mask select the indices of the localized spatial orbitals that belongs to the fragment.

Finally the embedding method can be run as usual:

```
result = dmet.run(fragment)
```

### One-shot DMET and full DMET¶

In practice, we maybe want to deal with larger molecules, and with multiple fragments. To handle multiple fragments
we need to use the one-shot DMET or full DMET methods, both covered by the class `inquanto.embeddings.dmet.DMETRHF`

.
Similarly to impurity DMET, `inquanto-pyscf`

includes a range of PySCF-driven
fragment solvers:

Here we consider the example of phenol:

```
from inquanto.geometries import GeometryMolecular
from inquanto.embeddings.dmet import DMETRHF
from inquanto.extensions.pyscf import (
get_fragment_orbital_masks,
get_fragment_orbitals,
ChemistryDriverPySCFMolecularRHF,
DMETRHFFragmentPySCFCCSD,
DMETRHFFragmentPySCFMP2,
)
# Phenol (C6H5OH)
# Setup the system - Phenol geometry
geometry = [['C', [-0.921240800, 0.001254500, 0.000000000]],
['C', [-0.223482600, 1.216975000, 0.000000000]],
['C', [ 1.176941000, 1.209145000, 0.000000000]],
['C', [ 1.882124000, 0.000763689, 0.000000000]],
['C', [ 1.171469000, -1.208183000, 0.000000000]],
['C', [-0.225726600, -1.216305000, 0.000000000]],
['O', [-2.284492000, -0.060545780, 0.000000000]],
['H', [-0.771286100, 2.161194000, 0.000000000]],
['H', [ 1.715459000, 2.156595000, 0.000000000]],
['H', [ 2.970767000, -0.000448048, 0.000000000]],
['H', [ 1.709985000, -2.155694000, 0.000000000]],
['H', [-0.792751600, -2.145930000, 0.000000000]],
['H', [-2.630400000, 0.901564000, 0.000000000]]]
g = GeometryMolecular(geometry)
driver = ChemistryDriverPySCFMolecularRHF(basis="sto-3g", geometry=g.xyz, charge=0)
hamiltonian_operator, space, rdm1 = driver.get_lowdin_system()
maskOH, maskCCH, maskCHCH1, maskCHCH2 = get_fragment_orbital_masks(
driver, [6, 12], [0, 1, 7], [2, 8, 3, 9], [4, 10, 5, 11]
)
dmet = DMETRHF(hamiltonian_operator, rdm1)
frOH = DMETRHFFragmentPySCFCCSD(dmet, maskOH)
frCCH = DMETRHFFragmentPySCFCCSD(dmet, maskCCH)
frCHCH1 = DMETRHFFragmentPySCFMP2(dmet, maskCHCH1)
frCHCH2 = DMETRHFFragmentPySCFCCSD(dmet, maskCHCH2)
fragments = [frOH, frCCH, frCHCH1, frCHCH2]
result = dmet.run(fragments)
print(result)
```

```
# STARTING CHEMICAL POTENTIAL 0.0
# STARTING CORR POTENTIAL PARAMETERS []
# FULL SCF ITERATION 0
```

```
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0
# FRAGMENT 0 - None: <H>=-301.7956772445182 EFR=-114.78345837272866 Q=-0.04983881865471318
```

```
# FRAGMENT 1 - None: <H>=-301.90654101502685 EFR=-152.8830668346832 Q=-0.005690013606114874
# FRAGMENT 2 - None: <H>=-301.9101535819737 EFR=-150.83827026648876 Q=-0.0014970562587865999
```

```
# FRAGMENT 3 - None: <H>=-301.9444865391344 EFR=-153.44135729451506 Q=-0.004749782126630464
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0001
# FRAGMENT 0 - None: <H>=-301.7965762435462 EFR=-114.78406015411932 Q=-0.04973186128709628
```

```
# FRAGMENT 1 - None: <H>=-301.9078353892906 EFR=-152.88478642548102 Q=-0.005423346447376076
# FRAGMENT 2 - None: <H>=-301.9115539401111 EFR=-150.83961285075645 Q=-0.0012572548814944184
```

```
# FRAGMENT 3 - None: <H>=-301.9458871008678 EFR=-153.4426208516093 Q=-0.004532549349235637
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.007436949983921135
# FRAGMENT 0 - None: <H>=-301.86256397736605 EFR=-114.82804988189318 Q=-0.04191500418033023
```

```
# FRAGMENT 1 - None: <H>=-302.0028761208733 EFR=-153.01088118923553 Q=0.014138160576203873
# FRAGMENT 2 - None: <H>=-302.01436028232354 EFR=-150.93805407032627 Q=0.016332731510335208
```

```
# FRAGMENT 3 - None: <H>=-302.0487052943609 EFR=-153.5353016245042 Q=0.011407964500872225
# CHEMICAL POTENTIAL 0.0074413042099870695
# FINAL PARAMETERS: []
# FINAL CHEMICAL POTENTIAL: 0.0074413042099870695
# FINAL ENERGY: -302.2586992530279
(-302.2586992530279, 0.0074413042099870695, array([], dtype=float64))
```

The utility function `get_fragment_orbital_masks()`

helps to create
the orbital masks from atom indices. In this particular examples, the 4
fragments are specified by lists of atom indices defined in the `geometry`

.
Consequently, `get_fragment_orbital_masks()`

returns 4
masks. Note however, that one could in principle define their own masks
which selects spatial orbitals within atoms or across multiple intra
atomic orbitals. Howeverm for `DMETRHF`

it is important the the set of
fragments completely cover the entire molecule.

## DMET with a custom solver¶

Running a VQE calculation on a DMET fragment of a large system requires the definition of a bespoke fragment solver.
There are many type of ansatze and strategies to perform VQE; instead of providing a corresponding range of many
black-box VQE fragment solvers, `inquanto-pyscf`

provides an easy way for a user to
define a fragment solver class.

To create a DMET fragment, we subclass the `inquanto.extensions.pyscf.DMETRHFFragmentPySCFActive`

class, which
requires us to implement the `solve_active()`

method. This
method assumes that the Hamiltonian, Fragment Energy operator and other arguments are only for the active
space specified when the fragment solver is constructed. The mandatory return is the expectation value of the
Hamiltonian and the Fragment Energy operator with the ground state (`energy`

and `fragment_energy`

) and the 1-RDM
of the ground state. Below, we outline the structure of these calculations, but full examples can be found on the examples
page.

```
class MyFragment(DMETRHFFragmentPySCFActive):
def solve_active( self,
hamiltonian_operator, fragment_energy_operator, fermion_space, fermion_state,
):
# ... your VQE solution ...
return energy, fragment_energy, vqe_rdm1
```

The active space can be specified with the CAS notation via the `FromActiveSpace`

class:

```
fr = MyFragment(dmet, mask, frozen=FromActiveSpace(2, 2))
```

or one can also provide the explicit list of indices of the frozen orbitals of the fragment embedding system. Note that the fragment solver in the case of DMET is solving the embedding system that is the fragment orbitals and the bath orbitals together.

For Impurity DMET the implementation is even simpler, only the energy is required as an output:

```
from inquanto.extensions.pyscf import ImpurityDMETROHFFragmentPySCFActive
from inquanto.operators import ChemistryRestrictedIntegralOperator
from inquanto.spaces import FermionSpace
from inquanto.states import FermionState
from inquanto.mappings import QubitMappingJordanWigner
from inquanto.ansatzes import FermionSpaceAnsatzUCCSD
from inquanto.computables import ExpectationValue
from inquanto.minimizers import MinimizerScipy
from inquanto.algorithms import AlgorithmVQE
from inquanto.protocols import SparseStatevectorProtocol
from pytket.extensions.qiskit import AerStateBackend
class MyFragment(ImpurityDMETROHFFragmentPySCFActive):
def solve_active(
self,
hamiltonian_operator: ChemistryRestrictedIntegralOperator,
fermion_space: FermionSpace,
fermion_state: FermionState,
):
jw = QubitMappingJordanWigner()
ansatz = FermionSpaceAnsatzUCCSD(fermion_space, fermion_state, jw)
h_op = jw.operator_map(hamiltonian_operator)
objective_expression = ExpectationValue(ansatz, h_op)
vqe = AlgorithmVQE(
objective_expression=objective_expression,
minimizer=MinimizerScipy(method="COBYLA", disp=True),
initial_parameters=ansatz.state_symbols.construct_random(0, 0.0, 0.01),
)
vqe.build(backend=AerStateBackend(), protocol_objective=SparseStatevectorProtocol())
vqe.run()
energy = vqe.final_value
return energy
```

## Other PySCF Hamiltonians for DMET¶

The API of DMET based embedding methods does not restrict the origin
of the Hamiltonian. As long as it is in a localized orthonormal basis
and its type is `ChemistryRestrictedIntegralOperator`

the Hamiltonian can be
constructed for a periodic system, or it can be a model Hamiltonian from the
Hubbard driver, or it can be generated with COSMO.

For example, if one would like to perform an embedding with a COSMO
calculation, then the Hamiltonian and the 1-RDM should be computed with
the `ChemistryDriverPySCFMolecularRHFQMMMCOSMO`

driver:

```
from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHFQMMMCOSMO
driver = ChemistryDriverPySCFMolecularRHFQMMMCOSMO(
basis="sto-3g", geometry=geometry, charge=0, do_cosmo=True
)
hamiltonian, fermion_space, dm = driver.get_lowdin_system()
# ... dmet and fragment definitions ...
energy, chemical_potential, parameters = dmet.run(fragments)
energy_water_frozen = energy + driver.cosmo_correction
```

Note

On the last line above, the COSMO energy is calculated by correcting the ground state energy of the Hamiltonian (which in this case is computed approximately by DMET) with the COSMO correction. This type of calculation assumes non self-consistent COSMO correction.

The COSMO method will be discussed in more detail below.

## FMO with a custom solver¶

In addition to DMET embedding methods, InQuanto supports FMO based embedding calculations with the important difference that the Hamiltonian needs access to the atomic orbital basis:

```
hamiltonian_operator_ao, _, _ = driver.get_system_ao(run_hf=False)
```

To run an FMO simulation, one needs to instantiate an `inquanto.extensions.pyscf.fmo.FMO`

class and, similarly to
DMET, define fragment solvers before finally calling the `run()`

method. In the
example below we create a custom FMO class and perform VQE on the fragment Hamiltonians.

```
from inquanto.extensions.pyscf.fmo import FMOFragmentPySCFActive, FMO
from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHF
from pytket.extensions.qiskit import AerStateBackend
class MyFMOFragmentVQE(FMOFragmentPySCFActive):
def solve_final_active(
self,
hamiltonian_operator,
fermion_space,
fermion_state,
) -> float:
from inquanto.ansatzes import FermionSpaceAnsatzUCCSD
qubit_operator = hamiltonian_operator.qubit_encode()
ansatz = FermionSpaceAnsatzUCCSD(fermion_space, fermion_state)
from inquanto.express import run_vqe
vqe = run_vqe(ansatz, hamiltonian_operator, AerStateBackend())
return vqe.final_value
driver = ChemistryDriverPySCFMolecularRHF(
geometry=[
["H", [0, 0, 0]],
["H", [0, 0, 0.8]],
["H", [0, 0, 2.0]],
["H", [0, 0, 2.8]],
["H", [0, 0, 4.0]],
["H", [0, 0, 4.8]],
],
basis="sto3g",
verbose=0,
)
full_integral_operator, _, _ = driver.get_system_ao(run_hf=False)
fmo = FMO(full_integral_operator)
fragments = [
MyFMOFragmentVQE(
fmo, [True, True, False, False, False, False], 2, "H2-1"
),
MyFMOFragmentVQE(
fmo, [False, False, True, True, False, False], 2, "H2-2"
),
MyFMOFragmentVQE(
fmo, [False, False, False, False, True, True], 2, "H2-3"
),
]
fmo.run(fragments)
```

```
# FIRST MONOMER ENERGIES: [-3.5180388948642447, -2.88601395248236, -1.686345871669464]
```

```
# CONVERGENCE REACHED!
# CONVERGED MONOMER ENERGIES: [-1.8009644337834718, -1.8134981363881528, -1.8009642542761402]
# TIMER BLOCK-0 BEGINS AT 2024-10-30 13:37:24.533702
```

```
# TIMER BLOCK-0 ENDS - DURATION (s): 0.5228625 [0:00:00.522863]
# TIMER BLOCK-1 BEGINS AT 2024-10-30 13:37:25.113326
```

```
# TIMER BLOCK-1 ENDS - DURATION (s): 0.3234876 [0:00:00.323488]
# TIMER BLOCK-2 BEGINS AT 2024-10-30 13:37:25.512686
```

```
# TIMER BLOCK-2 ENDS - DURATION (s): 0.5232339 [0:00:00.523234]
# FINAL MONOMER ENERGIES: [-1.8244970693394007, -1.8373292776161896, -1.824497051962279]
```

```
# TIMER BLOCK-3 BEGINS AT 2024-10-30 13:37:26.729120
```

```
# TIMER BLOCK-3 ENDS - DURATION (s): 12.8545473 [0:00:12.854547]
```

```
# TIMER BLOCK-4 BEGINS AT 2024-10-30 13:37:40.134625
```

```
# TIMER BLOCK-4 ENDS - DURATION (s): 10.5295081 [0:00:10.529508]
```

```
# TIMER BLOCK-5 BEGINS AT 2024-10-30 13:37:51.354352
```

```
# TIMER BLOCK-5 ENDS - DURATION (s): 12.9686133 [0:00:12.968613]
# FINAL DIMER ENERGIES: [-4.745456689215107, -4.18959129093186, -4.745456662426452]
# FINAL DIMER CORRECTION ENERGIES: [-1.0836303422595166, -0.5405971696301801, -1.0836303328479837]
# TOTAL ENERGY: -3.3512648044323345
```

```
-3.3512648044323345
```

## QM/MM¶

In QM/MM methods, the energy is partitioned into quantum mechanical (QM) and molecular mechanical (MM) contributions [64] corresponding to different parts of the system. Thus the QM/MM approach combines the accuracy of QM with the speed of MM. Typically, the QM subsystem is a small region where the interesting chemistry takes place, while the MM part corresponds to some larger environment, represented by simple point charges, surrounding the QM region. In general, there are two ways of expressing the total energy in QM/MM, “subtractive” and “additive” [65]. In InQuanto the additive approach is adopted, which means the total energy takes the form

where \(E_{\text{QM}}\) represents the contributions internal to the QM subsystem (calculated using high-level methods which include a two-body interaction or approximation thereof), \(E_{\text{MM}}\) comes from interactions between particles in the MM environment (a low cost classical term like the Coulomb interaction between point charges), and \(E_{\text{QM/MM}}\) represents the interactions between the QM and MM regions (also a classical interaction). The latter involves a modification of the one-body component of the electronic Hamiltonian. The modified one-body part of the Hamiltonian becomes

where \(\mathbf{x} = \{\mathbf{r}, s\}\) with \(\mathbf{r}\) (\(s\)) the electronic position (spin), and the modified classical Coulomb interaction for nuclei becomes

where \(M\) labels the MM point charges, with charge \(q_M\) and position \(\mathbf{R}_{M}\).

In `inquanto-pyscf`

the specialized driver class: `ChemistryDriverPySCFMolecularRHFQMMMCOSMO`

,
uses the QM/MM scheme of PySCF to modify the one-body integrals and nuclear Coulomb
interaction. With this driver, QM/MM embedding is performed as follows; we consider
an H_{2} molecule surrounded by ten random point charges:

```
from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHFQMMMCOSMO
# define the MM region
mm_charges = 0.1 # this will set all MM particles to have charge = +0.1 in units of electron charge.
# lists and numpy arrays are also accepted to specify charges individually.
mm_geometry = [ # pregenerated random geometry using nunmpy.random.seed(1)
[-0.49786797, 1.32194696, -2.99931375],
[-1.18600456, -2.11946466, -2.44596843],
[-1.88243873, -0.92663564, -0.61939515],
[0.2329004, -0.48483291, 1.111317],
[-1.7732865, 2.26870462, -2.83567444],
[1.02280506, -0.49617119, 0.35213897],
[-2.15767837, -1.81139107, 1.80446741],
[2.80956945, -1.11945493, 1.15393569],
[2.25833491, 2.36763998, -2.48973473],
[-2.7656713, -1.98101748, 2.26885502],
]
# e_mm_coulomb should be 0.068638356203663
# Execute the chemistry drivers.
driver_qmmm = ChemistryDriverPySCFMolecularRHFQMMMCOSMO(
zmatrix="""
H
H 1 0.7122
""",
basis="sto3g",
mm_charges=mm_charges,
mm_geometry=mm_geometry,
do_mm_coulomb=True,
do_qmmm=True,
)
```

The Hamiltonian returned by `get_system()`

will now
contain the modified one-body terms due to the QM/MM interaction.
To calculate \(E_{\text{tot}}\) as above including the Coulomb interaction between MM point
charges, the optional boolean keyword `do_mm_coulomb`

must be set to `True`

, which stores the MM Coulomb interaction
energy \(E_{\text{MM}}\) as a driver property: `e_mm_coulomb`

.
However, this does not add \(E_{\text{MM}}\) to the total energy, which must be done manually to recover
\(E_{\text{QM}} + E_{\text{QM/MM}} + E_{\text{MM}}\).

QM/MM can also be used with ROHF and UHF methods with the corresponding drivers
`ChemistryDriverPySCFMolecularROHFQMMMCOSMO`

and
`ChemistryDriverPySCFMolecularUHFQMMMCOSMO`

. In addition, QM/MM can be combined with other embedding schemes.

## COSMO¶

The COSMO (COnductor-like Screening MOdel) scheme [66] is a method for modelling the interaction of a molecule with a solvent. The solvent is approximated as a continuous medium (implicit solvation), and the interaction between the molecule and solvent corresponds to a modification of the one-body interaction. In InQuanto, the PySCF implementation of ddCOSMO (domain-decomposition COSMO) is utilized. In this scheme, the mean-field SCF problem is solved with a modified Hamiltonian

Where \(\hat{V}_{\text{ddCOSMO}}\) is a functional of the electronic density. This allows the orbitals to respond to the implicit solvation during the SCF procedure. However, \(E_{\text{SCF}}\) would not correspond to the total energy according to the ddCOSMO model. For this, a contribution \(E_{\text{ddCOSMO}}\) calculated from a classical dielectric problem is needed, in which the molecule is contained in a cavity defined by overlapping spheres centred on its atoms, and outside the cavity there is a continuous medium with dielectric constant fixed to that of the solvent

where \(f(\epsilon_s)\) is an empirical function of the solvent dielectric constant \(\epsilon_s\), \(\rho(\mathbf{r})\) is the electronic density, \(W(\mathbf{r})\) is a surface potential obtained from spherical harmonics, and the integral is over the surface defined by the cavity. For more details see [67].

Since \(E_{\text{ddCOSMO}}\) is treated as a classical constant term in the expectation value of the energy, the
quantity \(\langle \Psi|\hat{H}| \Psi \rangle + E_{\text{ddCOSMO}}\) would not contain the response of orbitals to
the solvent environment. Hence, the following expression for the total energy is adopted in PySCF and used in
the `inquanto-pyscf`

extension:

where \(D\) is the one-body reduced density matrix. Hence, the SCF problem is solved in the presence of the solvent potential which affects the orbitals, but this contribution is subtracted out and the classical COSMO energy \(E_{\text{ddCOSMO}}\) is added back in, so that the energy corresponds to the ddCOSMO model while the orbitals are also consistent with the background solvent potential.

The following example shows the instantiation of the `ChemistryDriverPySCFMolecularRHFQMMMCOSMO`

PySCF driver class for COSMO solvation:

```
from inquanto.extensions.pyscf import ChemistryDriverPySCFMolecularRHFQMMMCOSMO
# instantiation of driver with COSMO solvent
driver_cosmo = ChemistryDriverPySCFMolecularRHFQMMMCOSMO(
zmatrix="""
H
H 1 0.7122
""",
basis="sto3g",
do_cosmo = True
)
h, fsp, fst = driver_cosmo.get_system()
mf_solvent_e_tot = driver_cosmo.mf_energy
```

The Hamiltonian `h`

returned by `get_system()`

contains the modified one-body term, and construction of the quantum problem can procede as usual. The energy
`mf_solvent_e_tot`

corresponds to the mean field energy with the COSMO corrections as specified in (113).

It is also possible to combine COSMO with other embedding schemes in InQuanto by setting `do_cosmo=True`

in the driver as above.