Projection-based embedding

This tutorial aims to provide an introduction to projection-based embedding (Manby et al (2012)) for quantum chemistry calculations. Projection-based embedding, particularly within Density Functional Theory (DFT), is employed to address the computational challenges associated with large systems. This technique divides the system into a subsystem (e.g., a molecule) and an environment (e.g., crystal or solvent), treating them with different methods to optimize computational efficiency while maintaining accuracy for the region of interest. The fundamental concept involves projecting the wavefunction or electron density of the entire system onto the subsystem using projection operators. These operators define the interaction between the electronic structure of the subsystem and the surrounding environment.

The implementation involves partitioning the system, projecting the electronic information onto the subsystem using projection operators, and performing separate calculations for the subsystem (using DFT) and the environment (using a simpler method). The results are then combined to obtain the overall electronic structure of the system. The form of the projection operators depends on the specific details of the chosen embedding method. These operators encompass terms related to electron density, potential, or other relevant properties that describe the interaction between the subsystem and the environment.

Applications of projection-based embedding are prominent in materials science, especially in studying electronic properties in complex environments like surfaces or interfaces. The accuracy of this technique depends on the chosen methods for the subsystem and the approximations made in describing the environment.

Wavefunction-in-DFT Embedding

The projection-based approach readily allows for wavefunction-in-DFT (WF-in-DFT) embedding, in which subsystem \(\mathrm{A}\) is treated using a WF-level description and subsystem \(\mathrm{B}\) is described at the DFT level. Here, the WF-in-DFT energy is simply obtained by substituting the DFT energy of subsystem \(\mathrm{A}\) with the corresponding WF energy:

(111)\[\begin{split}\begin{aligned} E_{\mathrm{WF}-\mathrm{in}-\mathrm{DFT}}\left[\widetilde{\Psi}^\mathrm{A} ; \gamma^{\mathrm{A}}, \gamma^{\mathrm{B}}\right]= & E_{\mathrm{WF}}\left[\widetilde{\Psi}^{\mathrm{A}}\right]+\operatorname{tr}\left[\left(\widetilde{\gamma}^\mathrm{A} -\gamma^{\mathrm{A}}\right) \mathbf{v}_{\mathrm{emb}}\left[\gamma^{\Lambda}, \gamma^{\mathrm{B}}\right]\right]+ \\ & + E_{\mathrm{DFT}}\left[\gamma^{\Lambda}+\gamma^{\mathrm{B}}\right]-E_{\mathrm{DFT}}\left[\gamma^{\mathrm{A}}\right] +\mu \operatorname{tr}\left[\widetilde{\gamma}^{\Lambda} \mathbf{P}^{\mathrm{B}}\right]\end{aligned}\end{split}\]

where \(\widetilde{\Psi}^\mathrm{A}\) is the WF for subsystem \(\mathrm{A}\), \(\gamma^{\mathrm{A}}\) is the one-particle reduced density matrix corresponding to \(\widetilde{\Psi}^\mathrm{A}\), and \(E_{\mathrm{WF}}\left[\widetilde{\Psi}^{\mathrm{A}}\right]\) is the WF energy of subsystem \(\mathrm{A}\).

A projection-based WF-in-DFT embedding calculation proceeds as follows. A KS-DFT calculation is first performed over the full system. The resulting occupied MOs are localized and partitioned into two sets, corresponding to subsystems \(\mathrm{A}\) and \(\mathrm{B}\). These sets are used to construct \(\mathbf{h}^{\mathrm{A-in-B}}\),

(112)\[\mathbf{h}^{\mathrm{A-in-B}}\left[\gamma^{\mathrm{A}}, \gamma^{\mathrm{B}}\right]=\mathbf{h}+\mathbf{v}_{\mathrm{emb}}\left[\gamma^{\mathrm{A}}, \gamma^{\mathrm{B}}\right]+\mu \mathbf{P}^{\mathrm{B}}\]

which is an effective one-electron Hamiltonian containing the standard one-electron Hamiltonian, the embedding potential, and the projection operator. Finally, a correlated WF calculation is performed on subsystem \(\mathrm{A}\) wherein \(\mathbf{h}^{\mathrm{A-in-B}}\) replaces the standard one-electron Hamiltonian. The final WF-in-DFT energy is given by equation 1.

The WF calculation for subsystem \(\mathrm{A}\) consists of two steps: first, a set of reference orbitals is generated, and second, a correlated WF calculation is performed using those orbitals. The reference orbitals can be obtained either via Hartree–Fock (HF) or a multiconfigurational method. For the former case, the subsystem \(\mathrm{A}\) post-HF calculation begins with HF-in-DFT embedding. The HF-in-DFT Fock matrix, \(\mathbf{F}^\mathrm{A}\), is derived by inserting a Slater determinant for the subsystem \(\mathrm{A}\) WF into equation 1 and differentiating with respect to \({\widetilde{\gamma}_{\mathrm{HF}}^\mathrm{A}}\), giving

(113)\[\mathbf{F}^\mathrm{A}=\frac{\partial}{\partial \widetilde{\gamma}_{\mathrm{HF}}^\mathrm{A}} E_{\mathrm{HF}-\mathrm{in}-\mathrm{DFT}}\left[\tilde{\gamma}_{\mathrm{HF}}^\mathrm{A} ; \gamma^{\mathrm{A}}, \gamma^{\mathrm{B}}\right]=\mathbf{h}^{\mathrm{A}-\mathrm{in}-\mathrm{B}}\left[\gamma^{\mathrm{A}}, \gamma^{\mathrm{B}}\right]+\mathbf{g}\left[\tilde{\gamma}_{\mathrm{HF}}^{\mathrm{A}}\right ]\]

where \(\mathbf{g}\) includes all of the usual HF two-electron terms and \(\mathbf{h}^{\mathrm{A-in-B}}\) represents the effective one-electron Hamiltonian given by equation 2. Once the subsystem \(\mathrm{A}\) HF MOs are optimized in the presence of the DFT embedding potential, they are used for the correlated subsystem \(\mathrm{A}\) post-HF calculation. An analogous procedure holds for the case of multireference methods, wherein a multiconfigurational WF is substituted in place of the single Slater determinant in equation 1. In this way, projection-based WF-in-DFT embedding can be readily performed with any existing WF method (or quantum impurity solver) simply by modifying the one-electron Hamiltonian in the WF method to include the projection-based embedding terms. While projection-based embedding is exact for (same-functional) DFT-in-DFT embedding, projection-based WF-in-DFT embedding is necessarily approximate. For additional details, you can refer to Manby et al (2012).

WF-in-DFT in InQuanto

In InQuanto, WF-in-DFT is partially based on PsiEmbed, a computational chemistry software package tailored for quantum mechanical calculations, particularly focusing on projection-based embedding techniques. To demonstrate this embedding method, we employ the stretched ethane (C\(_2\)H\(_6\)) molecule as an example, as depicted below using the NGLView visualizer for InQuanto.

Here are the steps outlined:

• Define the system

• Define the embedded driver

• Run VQE to get the ground state energy using the embedded driver

import warnings
warnings.filterwarnings('ignore')

geometry = [
    ["C", [ 0.000000000, -0.664929641,  0.141563265]],
    ["C", [ 0.000000000,  0.664929641,  0.141563265]],
    ["H", [ 0.923341000, -1.237750972,  0.141563265]],
    ["H", [-0.923341000, -1.237750972,  0.141563265]],
    ["H", [ 0.923341000,  1.237750972,  0.141563265]],
    ["H", [-0.923341000,  1.237750972,  0.141563265]],
    ["H", [ 0.000000000,  0.000000000, -3.937392259]],
    ["H", [ 0.000000000,  0.000000000,  4.200000000]],
]

The NGLView visualizer for InQuanto can be used to display the system.

from inquanto.geometries import GeometryMolecular
from inquanto.extensions.nglview import VisualizerNGL

C2H6_geom = GeometryMolecular(geometry)
visualizer = VisualizerNGL(C2H6_geom)

#visualizer.visualize_molecule()

The initial step involves creating the InQuanto-PySCF projection-based embedding driver (inquanto.extensions.pyscf.ChemistryDriverPySCFEmbeddingROHF). This driver is essential for conducting molecular ROHF calculations and storing the outcomes in an InQuanto QubitOperator.

Parameters used:

geometry – Molecular geometry.

basis – Atomic basis set valid for Mole class.

multiplicity – Spin multiplicity of the total system.

frozen – Frozen orbital information.

functional – KS functional to use for the system calculation, or None if RHF is desired.

The FromActiveSpace function aids in determining the frozen orbital list based on the information provided about the active space. In this context, ncas represents the number of active orbitals, and nelecas represents the number of active electrons.

from inquanto.extensions.pyscf import ChemistryDriverPySCFEmbeddingROHF,FromActiveSpace

ncas, nelecas = 2, 2

driver = ChemistryDriverPySCFEmbeddingROHF(
    geometry=geometry,
    basis="3-21G*",
    multiplicity=3,
    frozen=FromActiveSpace(ncas, nelecas),
    functional="b3lyp5",
)

Additionally, it’s possible to compute the HF energy using the InQuanto driver.

driver.run_hf()

-79.09514941221715

The get_system function is responsible for computing the fermionic Hamiltonian operator, Fock space, and Hartree Fock state. The qubit_encode function carries out qubit encoding, utilizing the mapping class associated with the current integral operator. The default mapping approach is the Jordan-Wigner method.

chem_hamiltonian, space, state = driver.get_system()

qubit_hamiltonian = chem_hamiltonian.qubit_encode()

To construct our ansatz for the specified fermion space and fermion state, we have employed the Chemically Aware Unitary Coupled Cluster with singles and doubles excitations (UCCSD). The circuit is synthesized using Jordan-Wigner encoding.

from inquanto.ansatzes import FermionSpaceAnsatzChemicallyAwareUCCSD

ansatz = FermionSpaceAnsatzChemicallyAwareUCCSD(space, state)

Here, we carry out a straightforward VQE experiment to obtain the ground state energy of our system. For a more extensive guide on executing VQE calculations with InQuanto on quantum computers, we suggest referring to the VQE tutorial. To define where the computation is performed we set the backend to AerStateBackend().

from inquanto.express import run_vqe
from pytket.extensions.qiskit import AerStateBackend
from inquanto.minimizers import MinimizerRotosolve

backend = AerStateBackend()

vqe = run_vqe(
    ansatz,
    qubit_hamiltonian,
    backend=backend,
    with_gradient=False,
    minimizer=MinimizerRotosolve(),
)

# TIMER BLOCK-0 BEGINS AT 2024-09-23 15:34:32.188393
# TIMER BLOCK-0 ENDS - DURATION (s):  0.0667670 [0:00:00.066767]

Finally, we can print out the calculated energy value using VQE and compare the results with those obtained from CASCI and HF methods.

vqe.final_value

-79.15210375702853