5.4. Orbital entanglement analysis¶

5.4.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 phenomenon. The interaction between orbitals can be quantified through ideas in entanglement.

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 (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} = \frac{1}{2} \big(s(2)_{i,j} - s(1)_{i} - s(1)_{j} \big) \big(1 - \delta_{ij}\big),\]

where \(\delta_{ij}\) is the Kronecker delta.

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

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

5.4.2. Supported features¶

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

1. Calculation of single-orbital entropy and orbital-pair mutual information for a seniority-zero wavefunction

Supported wavefunction models are

AP1roG (seniority-zero wavefunction)

5.4.3. Orbital entanglement and correlation for a seniority zero wavefunction¶

If you use this module, please cite [boguslawski2015a]

5.4.3.1. How to set-up a calculation¶

To evaluate the single-orbital entropy and orbital-pair mutual information for a given wavefunction model, the corresponding one and two-particle RDMs need to be calculated first. Then one- and two-particle RDMs are used to evaluate the one- and two-orbital RDM’s from which the eigenvalues are needed to calculate the single-orbital and two-orbital entropy.

Given the one- and two-particle RDMs, an instance of the orbital_entanglement can be created. A function call of this instance, __call__(), calculates the entanglement and correlation measures.

5.4.3.2. How to generate correlation diagrams¶

To generate the single-orbital entropy diagram and mutual information plot, run in the terminal:

horton-entanglement.py threshold init_index final_index

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. The arguments init_index and final_index are optional arguments. If provided, the single-orbital entropy will be plotted for orbital indices in the interval [init_index, final_index].

5.4.4. Example Python scripts¶

5.4.4.1. Orbital entanglement analysis of an AP1roG wavefunction¶

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

data/examples/orbital_entanglement/water.py¶

from horton import *
from horton.test.common import numpy_seed

###############################################################################
## Set up molecule, define basis set ##########################################
###############################################################################
# Use the XYZ file from HORTON's test data directory.
fn_xyz = context.get_fn('test/water.xyz')
mol = IOData.from_file(fn_xyz)
obasis = get_gobasis(mol.coordinates, mol.numbers, 'cc-pvdz')
###############################################################################
## Define Occupation model, expansion coefficients and overlap ################
###############################################################################
lf = DenseLinalgFactory(obasis.nbasis)
occ_model = AufbauOccModel(5)
orb = lf.create_expansion(obasis.nbasis)
olp = obasis.compute_overlap(lf)
###############################################################################
## Construct Hamiltonian ######################################################
###############################################################################
kin = obasis.compute_kinetic(lf)
na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf)
er = obasis.compute_electron_repulsion(lf)
external = {'nn': compute_nucnuc(mol.coordinates, mol.pseudo_numbers)}
terms = [
    RTwoIndexTerm(kin, 'kin'),
    RDirectTerm(er, 'hartree'),
    RExchangeTerm(er, 'x_hf'),
    RTwoIndexTerm(na, 'ne'),
]
ham = REffHam(terms, external)
###############################################################################
## Perform initial guess ######################################################
###############################################################################
guess_core_hamiltonian(olp, kin, na, orb)
###############################################################################
## Do a Hartree-Fock calculation ##############################################
###############################################################################
scf_solver = PlainSCFSolver(1e-6)
scf_solver(ham, lf, olp, occ_model, orb)
###############################################################################
## Combine one-electron integrals to single Hamiltonian #######################
###############################################################################
one = kin.copy()
one.iadd(na)

###############################################################################
## Do OO-AP1roG optimization ##################################################
###############################################################################
ap1rog = RAp1rog(lf, occ_model)
with numpy_seed():  # reproducible 'random' numbers to make sure it always works
    energy, g, l = ap1rog(one, er, external['nn'], orb, olp, True)

###############################################################################
## Calculate response density matrices ########################################
###############################################################################
one_dm = lf.create_one_index()
one_dm.assign(orb.occupations)
twoppqq = lf.create_two_index()
twopqpq = lf.create_two_index()
ap1rog.compute_2dm(twoppqq, one_dm, g, l, 'ppqq')
ap1rog.compute_2dm(twoppqq, one_dm, g, l, 'pqpq')

###############################################################################
## Do orbital entanglement analysis ###########################################
###############################################################################
entanglement = OrbitalEntanglementAp1rog(lf, one_dm, [twoppqq,twopqpq])
entanglement()