15.2. Orbital entanglement analysis

15.2.1. Orbital entanglement and orbital correlation

In quantum chemistry, the interaction between orbitals is commonly used to understand chemical processes. For example, molecular orbital diagrams, frontier orbital theory, and ligand field theory all use orbitals to understand and justify chemical phenomena. The interaction between orbitals can be quantified using concepts of quantum information theory.

Specifically, the interaction between one orbital and all the other orbitals can be measured by the single-orbital entropy \(s(1)_i\) (or one-orbital entropy). It is calculated from the eigenvalues \(\omega_{\alpha,i}\) of the one-orbital reduced density matrix (1-RDM)

\[s(1)_i = -\sum_{\alpha=1}^4 \omega_{\alpha,i}\ln \omega_{\alpha,i}.\]

The one-orbital RDM is obtained from the N-particle RDM by tracing out all other orbital degrees of freedom except those of orbital i. This leads to an RDM whose dimension is equal to the dimension of the one-orbital Fock space. For spatial orbitals, the one-orbital Fock space has a dimension of 4 (one for unoccupied, one for doubly occupied and two for singly occupied).

Similarly, the two-orbital entropy \(s(2)_{i,j}\) quantifies the interaction of an orbital pair \(i,j\) and all other orbitals. It is calculated from the eigenvalues of the two-orbital RDM \(\omega_{\alpha, i, j}\) (with 16 possible states for spatial orbitals)

\[s(2)_{i,j} =-\sum_{\alpha=1}^{16} \omega_{\alpha, i, j} \ln \omega_{\alpha, i, j}.\]

The total amount of correlation between any pair of orbitals \((i,j)\) can be evaluated from the orbital-pair mutual information. The orbital-pair mutual information is calculated using the single- and two-orbital entropy and thus requires the one- and two-orbital RDMs,

\[I_{i|j} = s(2)_{i,j} - s(1)_{i} - s(1)_{j},\]

where we impose that \(i\neq j\), that is, we exclude self-interactions.

Note that a correlated wave function is required to have non-zero orbital entanglement and correlation. In the case of an uncorrelated wave function (for instance, a single Slater determinant) the (orbital) entanglement entropy is zero.

For more information on orbital entanglement and correlation, see refs. [boguslawski2015a] and [boguslawski2016b].

15.2.2. Supported features

Unless mentioned otherwise, the orbital entanglement module only supports restricted orbitals, DenseLinalgFactory and CholeskyLinalgFactory. The current version of PyBEST offers

Supported wave function models are

  1. pCCD (seniority-zero wavefunction)

  2. pCCD-LCCD

  3. pCCD-LCCSD

15.2.3. Seniority zero wavefunctions

If you use this module, please cite [boguslawski2015a], [boguslawski2016b], and [boguslawski2017b].

15.2.3.1. Quick Guide

To evaluate the single-orbital entropy and orbital-pair mutual information for a given (seniority-zero) wave function model, the corresponding one and two-particle RDMs need to be calculated first. Then the one- and two-particle RDMs are used to evaluate the one- and two-orbital RDMs whose eigenvalues are needed to calculate the single-orbital and two-orbital entropy.

In PyBEST, the corresponding RDMs (if available) are by default determined in each module and hence no additional steps are necessary.

15.2.3.1.1. pCCD

The single-orbital entropy and the orbital-pair mutual information for a pCCD wave function can only be determined if the ROOpCCD module is used as the (response) 1- and 2-RDM are only determined if the orbital-optimization is switched on. For the RpCCD module, the \(\Lambda\)-equations are not solved and hence no response density matrices can be calculated (see also The pCCD module) We thus assume that you have performed a ROOpCCD calculation, whose results are stored in the IOData container pccd_output. Furthermore, we will assume the following names for all PyBEST objects

lf:

A LinalgFactory instance (see Preliminaries).

The code snippet below shows how the single-orbital entropy and the orbital-pair mutual information for a pCCD wave function are calculated,

###############################################################################
entanglement = OrbitalEntanglementRpCCD(lf, pccd_output)

15.2.4. pCCD-LCC wavefunctions

If you use this module, please cite [nowak2020].

15.2.4.1. Quick Guide

To evaluate the single-orbital entropy and orbital-pair mutual information for a given LCC correction, you have to be perform a pCCD calculation followed by some pCCD-LCC calculation with the following keyword argument l=True (see \(\Lambda\)-equations) to solve for the \(\Lambda\)-amplitudes that are exploited to determine the reduced density matrices. After the pCCD-LCC calculations converged, the IOData container contains all necessary components of the 1-, 2-, 3-, 4-RDMs of the LCC correction and in addition the 1-, 2- RDMs of the pCCD method. To perform an orbital entanglement analysis for pCCD-LCC, execute

entanglement = OrbitalEntanglementRpCCDLCC(lf, pccdlccd_output)
entanglement()

for pCCD-LCCD and

entanglement = OrbitalEntanglementRpCCDLCC(lf, pccdlccsd_output)
entanglement()

for pCCD-LCCSD, respectively.

Note

Both the LCCD and LCCSD correction use the same orbital entanglement module. They differ in which IOData container is passed to the function call (either pccdlccd_output or pccdlccsd_output).

15.2.5. Output data generated by the Orbital Entanglement module

The Orbital Entanglement module will generate .dat files, which contain all unique values for the single-orbital entropy and the orbital-pair mutual information. For pCCD, all results are stored in s1-pccd.dat, while the latter is dumped to the i12-pccd.dat file. All pCCD-LCC results are saved in the s1-pccd-lcc.dat and i12-pccd-lcc.dat files. PyBEST supplies a script that will allow you to generate (separate) diagrams for both the single-orbital entropy and the orbital-pair mutual information (see Correlation diagrams for more details).

15.2.6. Correlation diagrams

To generate the single-orbital entropy diagram and the orbital-pair mutual information plot, execute the pybest-entanglement.py script shipped together with PyBEST,

pybest-entanglement.py threshold [--iname pybest-results/i12-pccd.dat --sname pybest-results/s1-pccd.dat]

where threshold determines the lower cutoff value of the mutual information and must be given in orders of magnitude (1, 0.1, 0.01, 0.001, etc.). Orbital correlations that are smaller than cutoff will not be displayed in the mutual information diagram.

This script offers additional optional arguments. To obtain more information on the pybest-entanglement.py script, you can display the help messages,

pybest-entanglement.py -h

By default, the single-orbital entropy is written to pybest-results/s1.pdf, while the orbital-pair mutual information plot is saved under pybest-results/i12.pdf.

Note

By default, the pybest-entanglement.py script reads the s1.dat and i12.dat data files. To use the data from a specific method, you have to use the –sname and –iname options when executing the script. See also the -h option for additional information on all possible options.

15.2.7. Example Python scripts

15.2.7.1. Orbital entanglement analysis of an pCCD wave function

This is a basic example of how to perform an orbital entanglement analysis in PyBEST. This script performs an orbital-optimized pCCD calculation, followed by an orbital entanglement analysis of the pCCD wave function for the water molecule using the cc-pVDZ basis set.

Listing 15.3 data/examples/orbital_entanglement/water.py
from pybest import context
from pybest.gbasis import (
    get_gobasis,
    compute_dipole,
    compute_eri,
    compute_kinetic,
    compute_nuclear,
    compute_nuclear_repulsion,
    compute_overlap,
)
from pybest.geminals import ROOpCCD
from pybest.cc import RpCCDLCCD, RpCCDLCCSD
from pybest.linalg import DenseLinalgFactory
from pybest.scf import AufbauOccModel
from pybest.wrappers import RHF
from pybest.orbital_entanglement import (
    OrbitalEntanglementRpCCD,
    OrbitalEntanglementRpCCDLCC,
)

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# Use the XYZ file from PyBEST's test data directory.
fn_xyz = context.get_fn("test/water.xyz")
obasis = get_gobasis("cc-pvdz", fn_xyz)
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb_a = lf.create_orbital(obasis.nbasis)
olp = compute_overlap(obasis)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = compute_kinetic(obasis)
ne = compute_nuclear(obasis)
er = compute_eri(obasis)
external = compute_nuclear_repulsion(obasis)

###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
hf = RHF(lf, occ_model)
hf_output = hf(kin, ne, er, external, olp, orb_a)

###############################################################################
## Do OO-pCCD optimization ####################################################
###############################################################################
pccd = ROOpCCD(lf, occ_model)
pccd_output = pccd(kin, ne, er, hf_output)

###############################################################################
## pCCD-LCCD calculation ######################################################
###############################################################################
lccd = RpCCDLCCD(lf, occ_model)
lccd_output = lccd(kin, ne, er, pccd_output, l=True)

###############################################################################
## pCCD-LCCSD calculation ######################################################
###############################################################################
lccsd = RpCCDLCCSD(lf, occ_model)
lccsd_output = lccsd(kin, ne, er, pccd_output, l=True)

###############################################################################
## Do orbital entanglement analysis for pCCD ###################################
###############################################################################
entanglement = OrbitalEntanglementRpCCD(lf, pccd_output)
entanglement()

###############################################################################
## Do orbital entanglement analysis for pCCD-LCCD ##############################
###############################################################################
entanglement = OrbitalEntanglementRpCCDLCC(lf, lccd_output)
entanglement()

###############################################################################
## Do orbital entanglement analysis for pCCD-LCCSD #############################
###############################################################################
entanglement = OrbitalEntanglementRpCCDLCC(lf, lccsd_output)
entanglement()