..
    : PyBEST: Pythonic Black-box Electronic Structure Tool
    : Copyright (C) 2016-- The PyBEST Development Team
    :
    : This file is part of PyBEST.
    :
    : PyBEST is free software; you can redistribute it and/or
    : modify it under the terms of the GNU General Public License
    : as published by the Free Software Foundation; either version 3
    : of the License, or (at your option) any later version.
    :
    : PyBEST is distributed in the hope that it will be useful,
    : but WITHOUT ANY WARRANTY; without even the implied warranty of
    : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    : GNU General Public License for more details.
    :
    : You should have received a copy of the GNU General Public License
    : along with this program; if not, see <http://www.gnu.org/licenses/>
    : --

Orbital entanglement analysis
#############################

.. _orbital_entanglement:

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 :math:`s(1)_i` (or one-orbital entropy).
It is calculated from the eigenvalues :math:`\omega_{\alpha,i}` of the
one-orbital reduced density matrix (1-RDM)

.. math::

    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 :math:`s(2)_{i,j}` quantifies the interaction
of an orbital pair :math:`i,j` and all other orbitals. It is calculated from the
eigenvalues of the two-orbital RDM :math:`\omega_{\alpha, i, j}` (with 16
possible states for spatial orbitals)

.. math::

    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 :math:`(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,

.. math::

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

where we impose that :math:`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]_.


Supported features
==================

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

    1. :ref:`Calculation of single-orbital entropy and orbital-pair mutual
       information for a seniority-zero wavefunction <orbital_entanglementseniorityzero>`
    2. :ref:`Calculation of single-orbital entropy and orbital-pair mutual information for pCCD-LCC methods <orbital_entanglementpCCD-LCC>`

Supported wave function models are

    1. :ref:`pCCD <user_pccd>` (seniority-zero wavefunction)
    2. :ref:`pCCD-LCCD <pccdlccd>`
    3. :ref:`pCCD-LCCSD <pccdlccsd>`


.. _orbital_entanglementseniorityzero:

Seniority zero wavefunctions
============================

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

.. _qg_oe_pccd:

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.

pCCD
^^^^

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

    :lf: A :py:class:`~pybest.linalg.base.LinalgFactory` instance (see :ref:`user_linalg_intro`).

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

.. literalinclude:: ../src/pybest/data/examples/orbital_entanglement/water.py
    :lines: 69-73



.. _orbital_entanglementpccd-lcc:

pCCD-LCC wavefunctions
======================

If you use this module, please cite [nowak2021]_.

.. _qg_oe_pccdlcc:

Quick Guide
-----------
To evaluate the single-orbital entropy and orbital-pair mutual information for a given LCC correction,
you have to perform a pCCD calculation followed by some pCCD-LCC calculation with the following keyword argument ``lambda_equations=True`` (see :math:`\Lambda`-:ref:`equations <lcc_lambda_equations>`) to solve for the :math:`\Lambda`-amplitudes that are exploited to determine the reduced density matrices. After the pCCD-LCC calculations converged, the :py:class:`~pybest.io.iodata.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


.. code-block:: bash

    entanglement = OrbitalEntanglementRpCCDLCC(lf, pccdlccd_output)
    entanglement()

for pCCD-LCCD and

.. code-block:: bash

    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 :py:class:`~pybest.io.iodata.IOData` container is passed to the function call (either ``pccdlccd_output`` or ``pccdlccsd_output``).

.. _output_results:

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 :ref:`correlation_diagrams` for more details).



.. _correlation_diagrams:

Correlation diagrams
====================

.. note::

    To use the visualization scripts shipped with PyBEST, you need to
    install the ``matplotlib`` package. Make sure that you have at least version 3.7.

    .. code-block:: bash

        # On Linux
        python3 -m pip install matplotlib~=3.7 --user

        # On MacOS
        pip install matplotlib~=3.7

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

.. code-block:: bash

    pybest-entanglement.py [--threshold 0.001 --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,

.. code-block:: bash

    pybest-entanglement.py -h

The most important arguments for controlling the functionality of the script are

    :threshold:  (float) lower cutoff value of the mutual information
    :iname:    (str) filename or path containing the mutual information data
    :sname:    (str) filename or path containing the single-orbital entropy data
    :indices:    (int [int int ...]) orbital indices to be considered in the mutual information plot
    :order:    (int [int int ...]) defines the new order of orbitals in the mutual information plot. Note that the indices start from 1, irrespective of the actual values for the orbital `indices` in the corresponding input files.
    :zoom:    (float) scaling factor used to scale MO pictures that are added to the mutual information plot

The script automatically arranges the pictures of the molecular orbitals around the mutual information plot.
It assumes that the orbital pictures are stored as ``mo_[index].png``, where ``[index]`` is the orbital index of the orbital in question (as labeled in the input files).
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``.
For additional examples on how to use the ``pybest-entanglement.py`` script, see the Examples below.


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

Example Python scripts
======================

Orbital entanglement analysis of a 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.

.. literalinclude:: ../src/pybest/data/examples/orbital_entanglement/water.py
    :caption: data/examples/orbital_entanglement/water.py
    :lines: 3-

After the calculations finish, you can plot the orbital correlation diagram (orbital-pair mutual information),
for instance for the pCCD-LCCSD wave function, executing the following command

.. code-block:: bash

    pybest-entanglement.py --threshold 0.001 --iname pybest-results/i12-pccd-lcc.dat --sname pybest-results/s1-pccd-lcc.dat

We can also print only a specific number of orbitals that feature the largest values of the mutual information
using the ``--largest-i [int]`` option. For instance, the snippet below restricts the final plot to at most 15 orbitals,

.. code-block:: bash

    pybest-entanglement.py --threshold 0.001 --largest-i 15 --iname pybest-results/i12-pccd-lcc.dat --sname pybest-results/s1-pccd-lcc.dat