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:
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 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
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 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
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 (121).
It is also possible to combine COSMO with other embedding schemes in InQuanto by setting do_cosmo=True
in the driver as above.