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

.. _user_properties_expectation_multipole:

Expectation-value Multipole Moments
###################################

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

PyBEST provides a computationally cheaper alternative to calculate the multipole moments,
such as the dipole moment and quadrupole moment, by using the expectation-value approach.
This method exploits an alternate excitation operator :math:`\hat{S}` defined as

.. math::

   e^{\hat{S}} |\Phi_0\rangle
   =
   \frac{e^{\hat{T}^{\dagger}} e^{\hat{T}} |\Phi_0\rangle}
   {\langle \Phi_0 | e^{\hat{T}^{\dagger}} e^{\hat{T}} | \Phi_0 \rangle}

.. math::

   \langle \hat{O} \rangle_{\mathrm{XCC}}
   =
   \langle \Phi_0 |
   e^{S^\dagger}
   e^{-T}
   \hat{O}
   e^{T}
   e^{-S^\dagger}
   | \Phi_0 \rangle

Based on the order-by-order expansion of :math:`\hat{S}`, and using appropriate truncation scheme,
an approximated form of the 1-RDM can be obtained directly from the cluster amplitudes without the need for solving the :math:`\Lambda`-equations.
The different blocks of the 1-RDM, then, can be expressed as

.. math::

   \gamma_i^a = t_i^a + \sum_{jb} t_j^b \left(2 t_{ij}^{ab} - t_{ij}^{ba}\right)\qquad \text{[Occ--Virt]} \\

   \gamma_k^i = \delta_i^k - \sum_{jab}t_{ij}^{ab}\left(2 t_{kj}^{ab} - t_{kj}^{ba}\right) \qquad \text{[Occ--Occ]} \\

   \gamma_c^a =\sum_{ijb} t_{ij}^{ab} \left(2 t_{ij}^{cb} - t_{ij}^{bc}\right) \qquad \text{[Virt--Virt]}\\


The expectation-value property module can be called as

.. code-block:: python

    expectation_multipole = ExpectationMultipole(lf, occ_model)

where

:lf: A :py:class:`~pybest.linalg.base.LinalgFactory` instance (see :ref:`user_linalg_intro`).
:occ_model: An Aufbau occupation model of the :py:class:`~pybest.scf.occ.AufbauOccModel` class.

Thereafter, the available multipole moments can be calculated as

.. code-block:: python

    multipole = expectation_multipole(
                method_output,
                property_options={"moment": "dipole", "quadrupole"}
                )

The ``method_output`` is the output from the electronic structure calculation as :py:class:`~pybest.io.iodata.IOData` object.
The available options for the ``property_options`` dictionary are listed below.
The ``moment`` key is mandatory, while the other keys are optional and have default values.

:moment: (str, list, tuple) the type of multipole moment to calculate, "dipole" or "quadrupole" or both as a tuple or list of strings.
:molecular_orbitals: (boolean) if True, the 1-RDM is transformed back to the atomic orbital basis (default True).
                    This option is important for post-Hartree-Fock methods where the 1-RDM is stored in the molecular orbital basis.
                    For the Hartree-Fock method, the 1-RDM is already stored in the atomic orbital basis, so this option should be set to False.
:scale: (float) the scaling factor for the 1-RDM (default 2.0). All internal 1-RDMs in PyBEST are stored for one set of orbitals.
        Thus, even in the restricted case, only one spin-component is saved. The total 1-RDM (spin-free) is obtained by multiplying by the scaling factor.
:origin: (list, tuple) the coordinates of the reference point for the multipole moment calculation (in a.u.). If not provided, the center of nuclear charge (COC) is
        used as the reference point.

The current version of PyBEST supports expectation-value 1-RDMs for the following post-Hartree-Fock methods:

    - pCCD
    - LCC (pCCD-LCCD, pCCD-LCCSD)
    - fpCC (fpCCD, fpCCSD)
    - CCD
    - CCSD


Example calculation
-------------------

Here we show an example of calculating the dipole and quadrupole moments for the water molecule at the pCCD and various post-pCCD levels of theory.

.. literalinclude:: ../src/pybest/data/examples/properties/water_expectation_multipole.py
    :caption: data/examples/properties/water_expectation_multipole.py
    :lines: 3-

.. note::
    1. The same procedure can be followed for the optimized pCCD natural orbital basis.
    In that case, the following changes should be made in the above code:

    .. code-block:: python

        from pybest.geminals import ROOpCCD

        # Converge pCCD in the optimized natural orbital basis
        pccd_opt = ROOpCCD(lf, occ_model)
        pccd_output = pccd_opt(kin, ne, eri, hf_output)

    2. The origin can be changed by providing the coordinates (in a.u.) of the reference point in the ``property_options`` dictionary.
    For example, to use the center of mass (COM) as the reference point, the following code can be used:

    .. code-block:: python

        from pybest.utils import calculate_com

        # Calculate the center of mass (COM) coordinates
        com_coordinates = calculate_com(molecule)

        # Calculate the multipole moments with COM as the reference point
        multipole = response_multipole(
                    method_output,
                    property_options={"moment": "dipole", "quadrupole", "origin": com_coordinates}
                    )