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 H2 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:

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 Li2 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.7956772445189  EFR=-114.78345837272859  Q=-0.04983881865476647 

# FRAGMENT 1 - None: <H>=-301.9065410150268  EFR=-152.88306683468306  Q=-0.005690013606173494 
# FRAGMENT 2 - None: <H>=-301.910153581972  EFR=-150.83827026648896  Q=-0.0014970562586942293 

# FRAGMENT 3 - None: <H>=-301.944486539135  EFR=-153.441357294514  Q=-0.004749782126747704 
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0001
# FRAGMENT 0 - None: <H>=-301.79657624354684  EFR=-114.78406015411913  Q=-0.04973186128714602 

# FRAGMENT 1 - None: <H>=-301.90783538929065  EFR=-152.8847864254808  Q=-0.005423346447429367 
# FRAGMENT 2 - None: <H>=-301.91155394010946  EFR=-150.83961285075648  Q=-0.0012572548814002715 

# FRAGMENT 3 - None: <H>=-301.94588710086833  EFR=-153.4426208516083  Q=-0.004532549349359982 
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.007436949983905783
# FRAGMENT 0 - None: <H>=-301.86256397736656  EFR=-114.82804988189295  Q=-0.04191500418040128 

# FRAGMENT 1 - None: <H>=-302.0028761208732  EFR=-153.0108811892352  Q=0.014138160576115055 
# FRAGMENT 2 - None: <H>=-302.01436028232143  EFR=-150.93805407032585  Q=0.016332731510379617 

# FRAGMENT 3 - None: <H>=-302.04870529436187  EFR=-153.5353016245028  Q=0.01140796450071413 
# CHEMICAL POTENTIAL 0.00744130421000467
# FINAL PARAMETERS: []
# FINAL CHEMICAL POTENTIAL: 0.00744130421000467
# FINAL ENERGY: -302.2586992530255
(-302.2586992530255, 0.00744130421000467, 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-07-19 14:20:08.853072

# TIMER BLOCK-0 ENDS - DURATION (s):  0.5225244 [0:00:00.522524]
# TIMER BLOCK-1 BEGINS AT 2024-07-19 14:20:09.431718

# TIMER BLOCK-1 ENDS - DURATION (s):  0.3224666 [0:00:00.322467]
# TIMER BLOCK-2 BEGINS AT 2024-07-19 14:20:09.836926

# TIMER BLOCK-2 ENDS - DURATION (s):  0.5195846 [0:00:00.519585]
# FINAL MONOMER ENERGIES: [-1.8244970693394003, -1.8373292776161896, -1.824497051962278]

# TIMER BLOCK-3 BEGINS AT 2024-07-19 14:20:11.026779

# TIMER BLOCK-3 ENDS - DURATION (s): 12.4277569 [0:00:12.427757]

# TIMER BLOCK-4 BEGINS AT 2024-07-19 14:20:23.982529

# TIMER BLOCK-4 ENDS - DURATION (s): 10.3037380 [0:00:10.303738]

# TIMER BLOCK-5 BEGINS AT 2024-07-19 14:20:34.955788

# TIMER BLOCK-5 ENDS - DURATION (s): 12.4024963 [0:00:12.402496]
# FINAL DIMER ENERGIES: [-4.7454566892151115, -4.189591290931863, -4.745456662426452]
# FINAL DIMER CORRECTION ENERGIES: [-1.083630342259522, -0.5405971696301846, -1.0836303328479846]
# TOTAL ENERGY: -3.351264804432345

-3.351264804432345

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

(116)\[E_{\text{tot}} = E_{\text{QM}} + E_{\text{QM/MM}} + E_{\text{MM}}\]

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

(117)\[h_{i,j} \rightarrow h_{i,j}^{\text{QM/MM}} = \int{d\mathbf{x}\phi_{i}^{*}(\mathbf{x})\bigg(\frac{-\nabla^{2}}{2}} - \sum_{A}\frac{Z_{A}}{|\mathbf{R}_{A} - \mathbf{r}|} - \sum_M \frac{q_M}{|\mathbf{R}_{M} - \mathbf{r}|} \bigg)\phi_{j}(\mathbf{x})\]

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

(118)\[V_{\text{nuc}} \rightarrow V_{\text{nuc}}^{\text{QM/MM}} = \sum_{A, B}\frac{Z_{A}Z_{B}}{|\mathbf{R}_{A} - \mathbf{R}_{B}|} + \sum_{A, M} \frac{Z_{A} q_M}{|\mathbf{R}_{A} - \mathbf{R}_{M}|}\]

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 H2 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

(119)\[E_{\text{SCF}} = \langle \Psi|\hat{H} + \hat{V}_{\text{ddCOSMO}} | \Psi \rangle\]

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

(120)\[E_{\text{ddCOSMO}} = \frac{1}{2}f(\epsilon_s) \int_{\Omega} \rho(\mathbf{r})W(\mathbf{r}) \,d\mathbf{r}\]

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:

(121)\[E_{\text{tot}} = \langle \Psi|\hat{H} + \hat{V}_{\text{ddCOSMO}} | \Psi \rangle - \text{Tr}(D \hat{V}_{\text{ddCOSMO}}) + E_{\text{ddCOSMO}}\]

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

It is also possible to combine COSMO with other embedding schemes in InQuanto by setting do_cosmo=True in the driver as above.