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.79567724451863 EFR=-114.7834583727295 Q=-0.049838818654846406
# FRAGMENT 1 - None: <H>=-301.9065410150269 EFR=-152.88306683468858 Q=-0.0056900136046422745
# FRAGMENT 2 - None: <H>=-301.91015358197444 EFR=-150.8382702664842 Q=-0.0014970562597600434
# FRAGMENT 3 - None: <H>=-301.9444865391351 EFR=-153.44135729451307 Q=-0.0047497821270710006
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0001
# FRAGMENT 0 - None: <H>=-301.7965762435465 EFR=-114.78406015412008 Q=-0.04973186128722773
# FRAGMENT 1 - None: <H>=-301.9078353892907 EFR=-152.88478642548625 Q=-0.005423346445912358
# FRAGMENT 2 - None: <H>=-301.911553940112 EFR=-150.839612850752 Q=-0.0012572548824536511
# FRAGMENT 3 - None: <H>=-301.9458871008683 EFR=-153.44262085160713 Q=-0.004532549349701043
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0074369499840891366
# FRAGMENT 0 - None: <H>=-301.86256397736804 EFR=-114.82804988189508 Q=-0.04191500418028582
# FRAGMENT 1 - None: <H>=-302.00287612087556 EFR=-153.01088118924383 Q=0.014138160578141878
# FRAGMENT 2 - None: <H>=-302.0143602823266 EFR=-150.93805407032397 Q=0.016332731509780984
# FRAGMENT 3 - None: <H>=-302.04870529436437 EFR=-153.53530162450397 Q=0.011407964500783407
# NEWTON ITERATION - CHEMICAL POTENTIAL 0.0074413042099937335
# FRAGMENT 0 - None: <H>=-301.8626031555281 EFR=-114.82807589346457 Q=-0.041910382966385384
# FRAGMENT 1 - None: <H>=-302.00293256711984 EFR=-153.01095598149362 Q=0.014149767620423859
# FRAGMENT 2 - None: <H>=-302.0144213310982 EFR=-150.9381124557574 Q=0.016343168386436346
# FRAGMENT 3 - None: <H>=-302.0487663483698 EFR=-153.53535661374167 Q=0.011417426261168373
# CHEMICAL POTENTIAL 0.007441306704681578
# FINAL PARAMETERS: []
# FINAL CHEMICAL POTENTIAL: 0.007441306704681578
# FINAL ENERGY: -302.2589134315259
(-302.2589134315259, np.float64(0.007441306704681578), 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. However 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: [np.float64(-3.518038894864243), np.float64(-2.8860139524823594), np.float64(-1.6863458716694648)]
# CONVERGENCE REACHED!
# CONVERGED MONOMER ENERGIES: [np.float64(-1.8009644337834712), np.float64(-1.8134981363881546), np.float64(-1.8009642542761397)]
-1.754715627117691
0.0
# TIMER BLOCK-0 BEGINS AT 2026-11-04 18:06:17.235474
# TIMER BLOCK-0 ENDS - DURATION (s): 0.5999092 [0:00:00.599909]
-1.785296660017343
0.0
# TIMER BLOCK-1 BEGINS AT 2026-11-04 18:06:18.437019
# TIMER BLOCK-1 ENDS - DURATION (s): 0.2589431 [0:00:00.258943]
-1.7547156112929558
0.0
# TIMER BLOCK-2 BEGINS AT 2026-11-04 18:06:19.816201
# TIMER BLOCK-2 ENDS - DURATION (s): 0.5787629 [0:00:00.578763]
# FINAL MONOMER ENERGIES: [np.float64(-1.8244970693393996), np.float64(-1.837329277616191), np.float64(-1.8244970519622776)]
-5.10457505374208
0.0
# TIMER BLOCK-3 BEGINS AT 2026-11-04 18:06:22.011214
# TIMER BLOCK-3 ENDS - DURATION (s): 6.4059473 [0:00:06.405947]
-4.455053374646447
0.0
# TIMER BLOCK-4 BEGINS AT 2026-11-04 18:06:29.417171
# TIMER BLOCK-4 ENDS - DURATION (s): 4.9807997 [0:00:04.980800]
-5.1045750116606
0.0
# TIMER BLOCK-5 BEGINS AT 2026-11-04 18:06:35.266001
# TIMER BLOCK-5 ENDS - DURATION (s): 5.8850437 [0:00:05.885044]
# FINAL DIMER ENERGIES: [np.float64(-4.745456689215109), np.float64(-4.189591290931863), np.float64(-4.745456662426452)]
# FINAL DIMER CORRECTION ENERGIES: [np.float64(-1.0836303422595184), np.float64(-0.5405971696301854), np.float64(-1.0836303328479837)]
# TOTAL ENERGY: -3.3512648044323416
np.float64(-3.3512648044323416)