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

Response Multipole Moments
##########################

PyBEST provides pathways to calculate various single-electron properties,
such as the dipole moment and quadrupole moment, by solving the response
:math:`\Lambda`-equations at pCCD and various post-pCCD levels of theory.
If you use this module, please cite [pccd-dipole-moments-jctc-2024]_.
For response properties with the pCCD-LCC methods, please also cite [nowak2021]_.

The response property module can be called as

.. code-block:: python

    response_multipole = ResponseMultipole(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 = response_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,
which contains the 1-RDM.
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 response 1-RDMs for the following post-Hartree-Fock methods:

    - MP2 (relaxed and unrelaxed)
    - pCCD
    - LCC (LCCD, LCCSD, pCCD-LCCD, pCCD-LCCSD)

.. note::

    For the LCC methods, the response 1-RDMs are only available when the :math:`\Lambda`-equations are solved by setting the keyword argument
    ``lambda_equations=True`` in the LCC function call.

    .. code-block:: python

        lcc=LCC(lf, occ_model)
        lcc_result = lcc(kin, ne, eri, reference_output, lambda_equations=True)

    where

    :kin:   the kinetic energy integrals.
    :ne:    the nucleus-electron attraction integrals.
    :eri:   the two-electron repulsion integrals.

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

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

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


Octupole moments can also be calculated by setting the ``moment`` key in the ``property_options`` dictionary to "octupole".
For example, the following code calculates the octupole moment at the HF and pCCD level of theory.

.. literalinclude:: ../src/pybest/data/examples/properties/water_response_octupole_pccd.py
    :caption: data/examples/properties/water_response_multipole_octupole.py
    :lines: 50-71


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

        # Calculate the center of mass (COM) coordinates
        com_coordinates = get_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}
                    )