Projection-based embedding with energy corrections

This can be viewed as a continuation of the projection-based embedding tutorial. In this tutorial, we will once again utilize projection-based embedding for quantum chemistry calculations within the context of Density Functional Theory (DFT). However, the key distinction here is the incorporation of NEVPT2 and AC0 energy correction methods to enhance the accuracy of our ground-state energy calculation. Incorporating these methods into our workflow requires us to redefine the final InQuanto driver used for this example because our previous approach did not yield the correct energy results. To demonstrate this, we have selected ethanol (C\(_2\)H\(_6\)O) as our test case.

To gain a better undersanding of the application of NEVPT2 and AC0 to enhance the precision of quantum computations beyond mean-field methods like Hartree-Fock or density functional theory (DFT), please consult the NEVPT2+AC0 tutorial. For further insights into Projection-based embedding, please refer to Manby et al (2012).

Here are the steps outlined:

• Define the system.

• Define the embedded RHF driver using AVAS.

• Generate a driver object wrapping the active space Hamiltonian.

• Run VQE to get the ground state energy and final parameters.

• Compute the NEVPT2 correction to the energy.

• Compute the AC0 correction to the energy.

The initial step involves creating the InQuanto-PySCF projection-based embedding driver (inquanto.extensions.pyscf.ChemistryDriverPySCFEmbeddingRHF). This driver is essential for conducting and storing the outcome of molecular RHF calculations.

Parameters used:

geometry – Molecular geometry.

basis – Atomic basis set valid for Mole class.

frozen – Frozen orbital information.

transf – Orbital transformer.

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

In order to reduce hardware resource requirements in multi-configuration and multireference electronic structure calculations, the Atomic Valence Active Space (AVAS) approximation in the PySCF extension of InQuanto was used. AVAS is an automated technique for creating active orbital spaces, particularly useful for describing electronic configurations arising from specific atomic valence orbitals, like metal d orbitals in coordination complexes. The variable aolabels contains the atomic orbital (AO) labels for the AO active space. The variables threshold and threshold_vir specify the truncation thresholds for the AO-projector. These thresholds determine which AOs are retained within the active space for occupied and virtual orbitals, respectively, based on their values. More detailed information can be found in Sayfutyarova et al (2017) and in the Fe4N2 system preparation tutorial.

import warnings
warnings.filterwarnings('ignore')

geometry = [
    ["O", [-1.1867, -0.2472, 0.0000]],
    ["H", [-1.9237, 0.3850, 0.0000]],
    ["H", [-0.0227, 1.1812, 0.8852]],
    ["C", [0.0000, 0.5526, 0.0000]],
    ["H", [-0.0227, 1.1812, -0.8852]],
    ["C", [1.1879, -0.3829, 0.0000]],
    ["H", [2.0985, 0.2306, 0.0000]],
    ["H", [1.1184, -1.0093, 0.8869]],
    ["H", [1.1184, -1.0093, -0.8869]],
]

basis="3-21G",
functional="b3lyp5"

from inquanto.extensions.pyscf import ChemistryDriverPySCFEmbeddingRHF, AVAS, FromActiveSpace

avas = AVAS(
    aolabels=["0 O 2p", "0 O 3p", "0 O 2s", "0 O 3s", "0 O 3d", "1 H"],
    threshold=0.8,
    threshold_vir=0.5,
)


driver = ChemistryDriverPySCFEmbeddingRHF(
    geometry=geometry,
    basis=basis,
    functional=functional,
    transf=avas,
    frozen=avas.frozenf,
)

To incorporate dynamic energy correlation methods, it is necessary to create a new driver object derived from the effective WF-in-DFT Hamiltonian. Hence, a PySCF driver object that encapsulates the Hamiltonian of the current active space is created using the InQuanto get_subsystem_driver function. This function accepts two arguments: the frozen orbital information (frozen), and the orbital transformer (transf).

FromActiveSpace 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. We could also perform postprocessing of orbitals using the Complete Active Space Self-Consistent Field (CASSCF) method to construct molecular integrals by enabling the casscf_transform parameter.

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 is Jordan-Wigner (used throughout this tutorial).

from inquanto.extensions.pyscf._transf import CASSCF

ncas, nelecas = 4, 4
casscf_transform = False

if casscf_transform:
    driver_embedded = driver.get_subsystem_driver(frozen=FromActiveSpace(ncas, nelecas), transf=CASSCF(ncas, nelecas))
else:
    driver_embedded = driver.get_subsystem_driver(frozen=FromActiveSpace(ncas, nelecas))

chem_hamiltonian, space, state = driver_embedded.get_system()
qubit_hamiltonian = chem_hamiltonian.qubit_encode()

When we examine the imaginary-coefficient terms of our Hamiltonian by using the hermitian_factorisation property, we observe that they are negligible for this active space. Consequently, we have constructed our qubit Hamiltonian by utilizing the hermitian part which exclusively includes terms with real coefficients.

hermitian_part, antihermitian_part = qubit_hamiltonian.hermitian_factorisation()
print(antihermitian_part)
qubit_hamiltonian_hermitian = qubit_hamiltonian.hermitian_part()

(-1.1641573924477122e-10, Y4 Z5 X6), (1.1641573924477122e-10, X4 Z5 Y6), (-1.1641573946161166e-10, Y5 Z6 X7), (1.1641573946161166e-10, X5 Z6 Y7)

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

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_hermitian,
    backend=backend,
    with_gradient=False,
    minimizer=MinimizerRotosolve(),
)

# TIMER BLOCK-0 BEGINS AT 2024-09-23 15:34:29.635203
# TIMER BLOCK-0 ENDS - DURATION (s): 143.6662325 [0:02:23.666232]

The InQuanto symmetry_operators_z2 function is employed to retrieve a list of symmetry operators applicable to our system. These symmetry operators are associated with the point group, spin parity, and particle number parity \(\mathbb{Z}_2\) symmetries that uphold a specific symmetry sector.

from inquanto.spaces import QubitSpace

symmetry_operators = QubitSpace(space.n_spin_orb).symmetry_operators_z2(
    qubit_hamiltonian_hermitian
)

In InQuanto, we can use a quantum computer to measure n-particle Reduced Density Matrices (n-RDMs) and employ them in a classical SC-NEVPT2 calculation, replacing the static correlation effects typically handled by CASCI with quantum simulations. The remaining dynamic electron correlation effects are approximated.

The class inquanto.computables.composite.PDM1234RealComputable computes the 1st, 2nd, 3rd, and 4th Pre-Density Matrices (PDMs) for a specific state, as determined by the provided ansatz and parameters.

Parameters used:

space – Fermion occupation space spanned by this RDM.

ansatz – Ansatz state with respect to which expectation values are computed.

encoding – Fermion to qubit mapping.

symmetry_operators – Z2 symmetries of the Hamiltonian.

cas_elec – Number of active electrons.

cas_orbs – Number of active orbitals.

The evaluation of PDMs is carried out using inquanto.protocols.SparseStatevectorProtocol, which is designed for sparse statevector calculations while utilizing caching. The get_evaluator method is employed to generate and provide a function (evaluator) that receives a specific quantum computable and computes it based on its type.

The get_nevpt2_correction function is used to compute the strongly contracted NEVPT2 correction to the energy using the provided density matrices.

For more extensive information about this method, please refer to Krompiec & Muñoz Ramo (2022).

Note that the next cell can take about 15 mins or more to run.

from inquanto.computables.composite import PDM1234RealComputable
from inquanto.mappings import QubitMappingJordanWigner
from inquanto.protocols import SparseStatevectorProtocol

pdm_computables = PDM1234RealComputable(
    space=space,
    ansatz=ansatz,
    encoding=QubitMappingJordanWigner(),
    symmetry_operators=symmetry_operators,
    cas_elec=ncas,
    cas_orbs=nelecas,
)

pdms = pdm_computables.evaluate(
    SparseStatevectorProtocol(backend).get_evaluator(vqe.final_parameters)
)

nevpt2_energy = driver_embedded.get_nevpt2_correction(pdms)

As an alternative approach, we can calculate the AC0 correction to the energy using the provided density matrices. Here, one- and two-particle reduced density matrices are employed to account for electron correlation effects. The inquanto.computables.composite.SpinlessNBodyRDMArrayRealComputable class was employed to calculate a general n-body RDM.

Parameters used:

n – n-body RDM.

fermion_space – Fermion space where the operators are defined.

ansatz – Ansatz with respect to which the expectation values are computed.

encoding – Qubit encoding from fermion space to qubit space.

symmetry_operators – List of Z2 symmetries of the Hamiltonian.

taperer – Optional taperer object.

The inquanto.protocols.SparseStatevectorProtocol class and its get_evaluator method are utilized in the same way as previously.

The get_ac0_correction function is used to calculate the AC0 correction to the energy from the provided density matrices.

To obtain both practical and theoretical insights, you can refer to the NEVPT2+AC0 tutorial, which provides an in-depth exploration of NEVPT2 and AC0 corrections within InQuanto.

from inquanto.computables.composite import SpinlessNBodyRDMArrayRealComputable
from inquanto.computables import ComputableList

rdm_computable = ComputableList(
    [
        SpinlessNBodyRDMArrayRealComputable(
            n=n,
            fermion_space=space,
            ansatz=ansatz,
            encoding=QubitMappingJordanWigner(),
            symmetry_operators=symmetry_operators,
            taperer=None,
        )
        for n in (1, 2)
    ]
)

rdm_computed = rdm_computable.evaluate(
    SparseStatevectorProtocol(backend).get_evaluator(vqe.final_parameters)
)


dm1 = rdm_computed[0]
dm2 = rdm_computed[1]

ac0_energy = driver_embedded.get_ac0_correction((dm1, dm2))

As the final step, we can conduct a comparison between our results and those obtained from classical calculations. Upon applying energy correction methods to the VQE result, it becomes evident that there is a very slight difference between the quantum and classical ground state energies.

print("FCI-in-DFT reference: {:.6f} Ha".format(-153.76095170040134))

print("VQE Energy: {:.6f} Ha".format(vqe.final_value))

print("NEVPT2 correction: {:.6f} Ha".format(nevpt2_energy))
print("Total VQE+NEVPT2 energy: {:.6f} Ha".format(vqe.final_value+nevpt2_energy))

print("AC0 correction using VQE PDMs: {:.6f} Ha".format(ac0_energy))
print("Total VQE+AC0 energy: {:.6f} Ha".format(vqe.final_value+ac0_energy))

FCI-in-DFT reference: -153.760952 Ha
VQE Energy: -153.754823 Ha
NEVPT2 correction: -0.004565 Ha
Total VQE+NEVPT2 energy: -153.759388 Ha
AC0 correction using VQE PDMs: -0.004790 Ha
Total VQE+AC0 energy: -153.759613 Ha