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

Geometry optimization with RHF and ROOpCCD wave functions
#########################################################

Before proceeding, please read the Quick Guide on geometry optimization
:ref:`here <user_geometry_intro>`.

The geometry optimization module uses geomeTRIC as the optimization backend.
At every optimization step, geomeTRIC proposes a new molecular geometry and
PyBEST evaluates the corresponding electronic energy and analytic nuclear
gradient.

.. _geometry_rhf:

RHF geometry optimization
=========================

RHF geometry optimizations are performed with
:py:class:`~pybest.geometry.geometry_optimizer.RHFGeometryOptimizer`.

The optimizer recomputes the RHF energy at every geometry and evaluates the
analytic RHF nuclear gradient. The one-particle density matrix required for the
gradient is generated from the converged RHF wave function.

A typical RHF geometry optimization has the following structure:

.. code-block:: python

    optimizer = RHFGeometryOptimizer(
        lf,
        occ_model,
        maxiter=100,
    )

    result = optimizer(
        hf_kwargs={
            "threshold": 1.0e-10,
        }
    )

The returned :py:class:`~pybest.io.iodata.IOData` container contains the optimized coordinates in the
``coordinates`` attribute.

.. _geometry_roopccd:

ROOpCCD geometry optimization
=============================

ROOpCCD geometry optimizations are performed with
:py:class:`~pybest.geometry.geometry_optimizer.ROOpCCDGeometryOptimizer`.

At each geometry, PyBEST first performs an RHF calculation to generate initial
orbitals. By default, the occupied and virtual orbitals are localized using the
Pipek--Mezey localization procedure before the ROOpCCD calculation is started.
The ROOpCCD energy and analytic gradient are then evaluated and passed to
geomeTRIC.

A typical ROOpCCD geometry optimization has the following structure:

.. code-block:: python

    optimizer = ROOpCCDGeometryOptimizer(
        lf,
        occ_model,
        maxiter=100,
    )

    result = optimizer(
        hf_kwargs={
            "threshold": 1.0e-10,
        },
        pccd_kwargs={
            "threshold": 1.0e-10,
            "localize": True,
        },
    )

To switch off Pipek--Mezey localization, set ``localize`` to ``False`` in
``pccd_kwargs``.

.. _geometry_keywords:

Summary of keyword arguments
============================

The geometry optimization classes accept the following geometry-related keyword
arguments during initialization:

    :coordsys: (str) Coordinate system passed to geomeTRIC. The default is
               ``"tric"``.

    :name: (str) Optional name used for output files and labels.

    :maxiter: (int) Maximum number of geometry optimization steps. The default
              is ``200``.

    :constraints: Optional geomeTRIC constraint specification. This can be
                  ``None`` or a path to a geomeTRIC constraints file.

    :hessian: (str) Hessian handling strategy passed to geomeTRIC. The default
              is ``"never"``.

    :transition: (boolean) If set to ``True``, geomeTRIC performs a
                 transition-state search. The default is ``False``.

    :check_interval: (int) Interval reserved for additional checks during the
                     optimization. The default is ``0``.

    :verbose: (boolean) If set to ``True``, more detailed optimization output is
              printed. The default is ``False``.

Additional keyword arguments passed when calling the optimizer are separated
into geometry-optimization options and wave-function options. Geometry options
are forwarded to geomeTRIC, while wave-function options are forwarded to the
energy calculation at each optimization step.

.. _geometry_wfn_keywords:

Wave-function keyword arguments
===============================

For RHF geometry optimization, SCF options can be passed through ``hf_kwargs``:

.. code-block:: python

    result = optimizer(
        hf_kwargs={
            "threshold": 1.0e-10,
            "maxiter": 128,
        }
    )

For ROOpCCD geometry optimization, options for the initial RHF calculation are
passed through ``hf_kwargs`` and options for the ROOpCCD calculation are passed
through ``pccd_kwargs``:

.. code-block:: python

    result = optimizer(
        hf_kwargs={
            "threshold": 1.0e-10,
        },
        pccd_kwargs={
            "threshold": 1.0e-10,
            "localize": True,
        },
    )

.. _geometry_constraints:

Constrained optimizations
=========================

Constraints can be provided through the ``constraints`` keyword. In typical
applications, this keyword points to a geomeTRIC constraints file.

.. code-block:: python

    optimizer = RHFGeometryOptimizer(
        lf,
        occ_model,
        constraints="constraints.txt",
    )

The exact syntax of the constraints file follows the
`geomeTRIC constraints format <https://github.com/leeping/geomeTRIC/blob/master/examples/constraints.txt>`_.

.. _geometry_transition_state:

Transition-state searches
=========================

Transition-state searches can be requested with the ``transition`` keyword:

.. code-block:: python

    optimizer = RHFGeometryOptimizer(
        lf,
        occ_model,
        transition=True,
        hessian="first+last",
    )

For transition-state optimizations, an appropriate initial structure and Hessian
strategy are important for convergence.

.. _geometry_output:

Output data
===========

The result of a geometry optimization is returned as a
:py:class:`~pybest.io.iodata.IOData` container. The optimized Cartesian
coordinates are available as

.. code-block:: python

    result.coordinates

During the optimization, PyBEST prints the energy and gradient norm for each
geometry step. After convergence, the final molecular geometry and internal
coordinates are printed.