..
: 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_wroking_with_orbitals:
Working with Orbitals in PyBEST
################################
The (molecular) orbitals in PyBEST contain information on the AO/MO coefficient
matrix, the orbital occupation numbers, and the orbital energies. Internally,
the orbitals
are an instance of the :py:class:`~pybest.linalg.dense.DenseOrbital` class.
Specifically, the AO/MO coefficients, occupation numbers, and energies are
attributes of the :py:class:`~pybest.linalg.dense.DenseOrbital` class.
In the following, some simple orbital manipulations are explained that might be
useful for users.
Creating an empty set of orbitals
=================================
Prior to any wave function calculation, the orbitals need to be created as an
instance of the :py:class:`~pybest.linalg.dense.DenseOrbital` class
.. code-block:: python
# Create an empty set of orbitals
orb_a = lf.create_orbital()
where ``lf`` is an instance of the :py:class:`~pybest.linalg.LinalgFactory`
class, where all linear algebra and tensor contraction operations are defined.
Printing orbital information
============================
The AO/MO coefficient matrix, the orbital occupation numbers, and the orbital
energies can be accessed as follows (assuming ``orb`` to be an instance of
:py:class:`~pybest.linalg.dense.DenseOrbital`)
.. code-block:: python
# Print AO/MO coefficients
print(orb.coeffs)
# Print occupations
print(orb.occupations)
# Print orbital energies
print(orb.energies)
The AO/MO coefficients can also be printed in a more readable format using the
:py:meth:`~pybest.orbital_utils.print_ao_mo_coeffs` function
.. code-block:: python
# Print all AO/MO coefficients in a more readable format
print_ao_mo_coeffs(gobasis, orb)
where ``gobasis`` is an instance of the ``GOBasis`` class. If only a subset of
orbitals is required, this can be achieved using optional arguments
.. code-block:: python
# Print AO/MO coefficients in a more readable format for orbitals
# no. 4 to no. 31 (including)
print_ao_mo_coeffs(gobasis, orb, begin=4, end=31)
where ``begin`` and ``end`` are the first and last index, respectively, of the
orbital to be printed. Note that the orbitals are indexed using Fortran
indexing, that is, beginning with 1. The order of the AOs is according to the list of atoms in the
geometry input. In the output, each column corresponds to one molecular orbital.
.. note::
In all electronic structure methods that optimize the orbitals, the current
orbitals are always overwritten with the most recent ones.
To save a copy of the current orbitals during the execution of the PyBEST code,
use the :py:meth:`~pybest.linalg.DenseOrbital.copy` method of
:py:class:`~pybest.linalg.dense.DenseOrbital`,
.. code-block:: python
# Create a true copy of the current orbitals
orb_output = orb.copy()
Writing orbitals to a checkpoint file
=====================================
In PyBEST, all electronic structure methods save their own copy of the molecular
orbitals and dump them in their method-dependent checkpoint files.
In case, you want to write the molecular orbitals to a separate file at any
stage of the program (represented in PyBEST's internal format), you can use
the :py:class:`~pybest.io.iodata.IOData` container (see also
:ref:`user_iodata` for more details on how to use the container),
.. code-block:: python
# Dump orbitals for restart
mol = IOData(olp=olp, orb_a=orb)
mol.to_file('checkpoint.h5')
where the overlap matrix is stored in ``olp`` and the orbitals in ``orb_a``.
The molecular orbitals (for restricted orbitals or alpha spin)
are always associated with ``orb_a``, while the
overlap matrix is associated with ``olp``. In case of unrestricted orbitals,
PyBEST uses the attribute name ``orb_b``. If the user chooses a different
naming convention, PyBEST might run into unexpected errors (for instance,
during restarts). See also :ref:`user_iodata` for more information on how
to use the IOData container and all default variable names.
.. note::
To use previously dumped orbitals for restart purposes, the overlap matrix
has to be written to file as well.
.. _user_orbitals_import:
Reading orbitals from a checkpoint file
=======================================
Most electronic structure methods in PyBEST feature a restart option, which
reads and projects some previously stored orbitals automatically. The
documentation below is mentioned for reasons of completeness to allow the user
to manually read in orbitals from disk.
.. _orbitals_restart_pybest_format:
PyBEST's internal format
-------------------------
Similar to the recipe of dumping molecular orbitals to file, they can be read
from disk using the :py:class:`~pybest.io.iodata.IOData` container functionality.
First, some previously stored PyBEST checkpoint file that contains the molecular
orbitals and overlap integrals has to be read in,
.. code-block:: python
# Read some PyBEST checkpoint file and store the corresponding data in 'old'
old = IOData.from_file('checkpoint.h5')
The molecular orbitals and overlap matrix can be accessed from ``old`` as
``old.orb_a`` and ``old.olp``, respectively.
If the molecular geometry of the restart orbitals differs from the current
molecular geometry, the molecular orbitals have to be projected against the
current molecular orbitals ``orb`` and current overlap matrix ``olp``,
.. code-block:: python
# Project orbitals read from file against current orbitals
project_orbitals(old.olp, olp, old.orb_a, orb)
The :py:meth:`~pybest.orbital_utils.project_orbitals` function overwrites
``orb`` and ``olp`` with the projected restart orbitals and overlap matrix,
respectively.
.. note::
If the restart orbitals are obtained for a different molecular geometry,
the :py:meth:`~pybest.orbital_utils.project_orbitals` function
has to be used, otherwise all electronic structure calculations will be
corrupted.
External file format
--------------------
Orbital information can also be imported from external programs. Supported
file formats are ``.mkl`` and ``.molden``. Similar to the procedure mentioned
above (:ref:`orbitals_restart_pybest_format`), the :py:class:`~pybest.io.iodata.IOData`
container is used to read in some external file. The file extension is handled
automatically.
.. code-block:: python
# Load molden file file
mol = IOData.from_file('water.molden')
# Print the total number of molecular orbitals
print(mol.orb_a.nfn)
Similarly, if the orbitals are not projected onto a new basis set, the
Hamiltonian has to be recomputed for the external basis set stored as
``mol.gobasis``.
The mol instance has the following attributes
:mol.gobasis: A :py:class:`~pybest.gbasis.pybasis.Basis` instance. It contains
all basis set information.
:mol.orb_a: A :py:class:`~pybest.linalg.dense.DenseOrbital` instance containing
the molecular orbitals for alpha spin
:mol.orb_b: A :py:class:`~pybest.linalg.dense.DenseOrbital` instance containing
the molecular orbitals for beta spin (if present)
:mol.coordinates: A np.array containing the xyz coordinates
:mol.atom: A list of strings with the element names
.. _orbitals_rotations:
Rotating orbitals
=================
The molecular orbitals can be perturbed by various rotations. The user has the
choice between (a series) of 2-by-2 rotations or a random rotation of all
molecular orbitals.
Givens rotation (2-by-2 rotations)
----------------------------------
Two molecular orbitals can be rotated using the
:py:meth:`~pybest.linalg.dense.DenseOrbital.rotate_2orbitals` method.
Depending on the angle of rotation, the selected orbitals are either mixed or
swapped. By default, the function rotates the HOMO and LUMO orbitals by 45 degrees.
.. code-block:: python
# Mix HOMO and LUMO orbitals
orb.rotate_2orbitals()
If a different pair of orbitals is to be rotated, the corresponding orbital
indices can be specified as optional arguments,
.. code-block:: python
# Rotate 7th and 12th orbital by 60 deg; note zero-based indexing
orb.rotate_2orbitals(60*deg, 6, 11)
Keep in mind that PyBEST uses radians as unit for angles, that is,
``60*deg == 60*np.pi/180``.
Since :py:meth:`~pybest.linalg.dense.DenseOrbital.rotate_2orbitals` is a class
method (defined in :py:class:`~pybest.linalg.dense.DenseOrbital`),
all quantities use Python-based indexing, starting with 0.
The :py:meth:`~pybest.linalg.dense.DenseOrbital.rotate_2orbitals` method can
be applied in sequence if several orbitals are to be rotated with each other,
.. code-block:: python
# Perform series of rotations: 7/12 by 60 degree, 8/14 by 20 degree, 9/26
# by 80 degree
orb.rotate_2orbitals(60*deg, 6, 11)
orb.rotate_2orbitals(20*deg, 7, 13)
orb.rotate_2orbitals(80*deg, 8, 25)
Note that an angle of 90 degrees is equivalent to swapping the corresponding
orbitals. This functionality can also be achieved using the
:py:meth:`~pybest.linalg.dense.DenseOrbital.swap_orbitals` method
(see :ref:`orbitals_swapping` below).
Random orbital rotations
------------------------
All molecular orbitals can be perturbed simultaneously using the
:py:meth:`~pybest.linalg.dense.DenseOrbital.rotate_random` method.
This function method performs a *random unitary rotation* to **all**
molecular orbitals.
.. code-block:: python
# perform random unitary rotation of orbitals orb
orb.rotate_random()
.. note::
All orbital rotation operations leave the occupation numbers and orbital
energies unchanged. Only the MO coefficients are perturbed. The AOs remain
always untouched.
.. _orbitals_swapping:
Swapping orbitals
=================
Two molecular orbitals can be swapped using the
:py:meth:`~pybest.linalg.dense.DenseOrbital.swap_orbitals` method.
The list of orbital swaps is passed as a 2-dimensional numpy array, where each
row corresponds to one transposition of orbitals. Since
:py:meth:`~pybest.linalg.dense.DenseOrbital.swap_orbitals` is again a class
method, the indexing starts with 0.
.. code-block:: python
# Swap some orbitals; note zero-based indexing
swaps = np.array([[0, 6], [5, 9]])
orb.swap_orbitals(swaps)
By default, the orbital occupations and orbital energies are swapped as well.
In some cases (like in SCF calculations), the swapping of orbital occupations
and/or orbital energies should be suppressed to provide a perturbed guess
for the molecular orbitals. These features can be switched on using optional
arguments,
.. code-block:: python
# Swap only AO/MO coeffients and leave occupations and energies unchanged;
# note zero-based indexing
swaps = np.array([[0, 6], [5, 9]])
orb.swap_orbitals(swaps, skip_occs=True, skip_energies=True)
.. _user_orbitals_localization:
Orbital localization
====================
The current version of PyBEST supports only Pipek-Mezey localization.
For more details on Pipek-Mezey localization please refer to :ref:`localization`.
Due to reasons of completeness, we will only briefly explain the
syntax on how to perform a localization of molecular orbitals.
In the Pipek-Mezey localization algorithm, the occupied and virtual block of the
(canoncial Hartree-Fock) orbitals are localized using Mulliken projectors.
In order to perform a localization, the Mulliken projectors need to be determined.
To do so, an instance of the :py:class:`~pybest.gbasis.pybasis.Basis` class
(here ``gobasis``) and of the :py:class:`~pybest.linalg.base.LinalgFactory`
class (here ``lf``) are required.
.. literalinclude:: ../src/pybest/data/examples/localization/water_pm.py
:lines: 46-53
As soon as the Mulliken projectors are calculated, they can be passed to the
:py:class:`~pybest.localization.localizaton.PipekMezey` class. The occupied and
virtual blocks can be localized individually by two separate function calls,
.. literalinclude:: ../src/pybest/data/examples/localization/water_pm.py
:lines: 54-61
.. note::
The orbitals have to passed as an ``IOData`` container, here ``hf_results``. The corresponding
attribute ``orb_a`` will be overwritten in the localization procedure.
Frozen core orbitals are also supported during the localization of the occupied
orbitals. To define a frozen core, the ``ncore`` argument has to be defined during
the initialization of the :py:class:`~pybest.localization.localizaton.PipekMezey`
class
.. literalinclude:: ../src/pybest/data/examples/localization/water_pm_fc.py
:lines: 50-53
The actual localization is invoked as in the case of no frozen core orbitals.
Example Python scripts
----------------------
The following, we list a complete example script for a restricted Hartree-Fock
calculation of the water molecule followed by orbital localization.
.. literalinclude:: ../src/pybest/data/examples/localization/water_pm.py
:caption: data/examples/localization/water_pm.py
:lines: 3-45
.. _user_orbitals_molden:
Generating molden files
=======================
PyBEST uses libint2 to dump molden files. As all I/O operations, the generation
of molden files is handled by the :py:class:`~pybest.io.iodata.IOData` container
class. The corresponding class instance has to contain an instance of the
orbitals, for which the molden file has to be generated, as well as an instance
of the :py:class:`~pybest.gbasis.pybasis.Basis` class. The latter contains
all necessary information to dump the molden file.
.. literalinclude:: ../src/pybest/data/examples/localization/water_pm.py
:lines: 46-
Specifying the file ending ``molden`` tells the
:py:meth:`~pybest.io.iodata.IOData.to_file` method to dump a molden file to
disk.
In the case of unrestricted orbitals, both alpha and beta orbitals have to be
passed to the :py:class:`~pybest.io.iodata.IOData` container.
We can also construct an :py:class:`~pybest.io.iodata.IOData` container
from scratch that contains all required information,
.. code-block:: python
molden = IOData(gobasis=gobasis,
orb_a=orb_a,
orb_b=orb_b
)
molden.to_file("water.molden")