..
: 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_scf_preliminaries:
Preliminaries
=============
The SCF module in PyBEST contains various solvers to solve the Schrödinger
equation self-constistently. The same solvers are used in *restricted* and
*unrestricted* Hartree-Fock calculations.
To perform an SCF calculation, PyBEST requires
- In case of molecular calculation, the molecular geometry, basis set, and molecular
Hamiltonian need to be present.
- For model Hamiltonians, the corresponding Hamiltonian matrix elements have to be
available.
See also :ref:`user_hamiltonian` for more details.
.. note::
Before setting up an SCF calculation, make sure that you have followed the
instructions in :ref:`user_hamiltonian`.
In PyBEST, all electronic structure methods feature a similar structure.
The individual steps of an SCF calculation are
1. create an instance of the :py:class:`~pybest.wrapper.hf.RHF` or
:py:class:`~pybest.wrapper.hf.UHF` class
2. initiate the SCF optimization using a function call
After the SCF is converged, the final data can be processed to
1. dump orbitals (for instance into molden files)
2. pass output to post-HF methods
Throughout this documentation, all electronic structure methods will label
all objects according to the PyBEST name convention. See also
:ref:`user_iodata` for a list of used acronyms in PyBEST.
Note that all variables listed below have to be initialized and/or evaluated
beforehand (see also :ref:`user_hamiltonian`). Below, we assume that the following
quantities have been defined prior entering the SCF module
:lf: the LinalgFactory used to handle all tensor contractions
:occ_model: the occupation model used
:kin: the kinetic energy integrals
:ne: the nucleus-electron attraction integrals
:eri: the two-electron repulsion integrals
:external: some external potential, for instance the nuclear repulsion term
:olp: the overlap integrals of the AO basis
:orb_a: the molecular orbitals for alpha spin (or for restricted orbitals)
:orb_b: the molecular orbitals for beta spin (if present)
.. _user_scf_hf:
The Hartree-Fock wrapper
========================
PyBEST provides a wrapper that allows the user to easily define and perform
Hartree-Fock calculations. It automatically generates some initial guess orbitals,
chooses some DIIS flavor, and dumps and returns all Hartree-Fock output data
required for restarts and post-processing.
.. _user_scf_hf_rhf:
Restricted Hartree-Fock calculations
------------------------------------
The code snippet below shows how to perform a restricted Hartree-Fock calculation
in PyBEST,
.. literalinclude:: ../src/pybest/data/examples/hf/rhf_water_dense.py
:lines: 51-54
.. note::
All arguments (all integrals and orbitals) can be passed in any order. The
:py:class:`~pybest.wrapper.hf.RHF` class determines their type automatically.
By default, PyBEST performs a core Hamiltonian guess, where the one-electron
Hamiltonian is diagonalized. Specifically, the Hamiltonian contains only all
one-electron terms that have been passed to the function call ``hf()``,
while all two-electron interactions are neglected.
Due to the lack of the electron-electron interaction terms, the eigenfunctions
of the one-electron Hamiltonian can be obtained by simply diagonalizing the
core Hamiltonian. These orbitals can then be used as guess orbitals for any
SCF procedure.
.. note::
The Hartree-Fock method overwrites the initial guess orbitals ``orb_a``
with the corresponding canonical HF orbitals. It also updates the
occupation numbers and the orbital energies.
The :py:class:`~pybest.wrapper.hf.RHF` class returns an instance of the
:py:class:`~pybest.io.iodata.IOData` container class, which stores all important
Hartree-Fock results as attributes. This container can be passed as argument
to other (post-Hartree-Fock or post-processing) methods to provide for instance
(guess) orbitals, reference energies, etc.
Specifically, the ``hf_`` output container possesses the following attributes
by default
:hf\_.converged: if converged, returns True (boolean)
:hf\_.olp: the overlap integrals
:hf\_.dm_1: the Hartree-Fock 1-particle reduced density matrix (RDM) in the AO
basis
:hf\_.orb_a: the canonical Hartree-Fock orbitals (stored as a true/deep copy)
:hf\_.e_tot: the total Hartree-Fock energy
:hf\_.e_ref: the energy of the Hartree-Fock reference determinant (equivalent
to e_tot)
:hf\_.e_kin: the kinetic energy
:hf\_.e_ne: the nuclear attraction energy
:hf\_.e_hartree: the Coulomb energy
:hf\_.e_x\_hf: the Hartree-Fock exchange energy
:hf\_.e_core: the energy associated with some external potential
By default, the results of a Hartree-Fock calculation are stored to disk
using PyBEST's internal format and can be used for restarts.
The default filename is ``checkpoint_scf.h5``.
.. note::
All checkpoint files are stored in PyBEST's results directory called
``pybest-results``. This directory can be changed globally to any directory specified
by the user using the ``filemanager``,
.. code-block:: python
from pybest.filemanager import filemanager
# "some_other_directory" will be created related to the PyBEST run script
filemanager.result_dir = "some_other_directory"
.. _user_scf_hf_uhf:
Unrestricted Hartree-Fock calculations
--------------------------------------
To perform unrestricted Hartree-Fock calculation in PyBEST, the above code
snippet has to be modified only slightly,
.. literalinclude:: ../src/pybest/data/examples/hf/uhf_methyl_dense.py
:lines: 48-50
Thus, an instance of the :py:class:`~pybest.wrapper.hf.UHF` class has to be
created and the beta orbitals have to be passed as an additional argument.
.. note::
PyBEST assumes that the first orbitals passed as arguments correspond to
alpha electrons, while the second set of orbitals is associated with
beta electrons.
The corresponding :py:class:`~pybest.io.iodata.IOData` container (return value)
also contains the beta orbitals and the 1-particle RDM of the alpha and beta electrons
as attribute,
:hf\_.dm_1_a: the Hartree-Fock 1-particle RDM in the AO basis for alpha electrons
:hf\_.dm_1_b: the Hartree-Fock 1-particle RDM in the AO basis for beta electrons
:hf\_.orb_b: the canonical Hartree-Fock orbitals (stored as a true/deep copy)
for beta electrons
.. _user_scf_convergence:
Steering the SCF optimization
=============================
The convergence behaviour of the SCF algorithm can be steered using keyword
arguments. The procedure is similar for both **unrestricted** and **restricted**
Hartree-Fock calculations.
.. _user_scf_convergence_threshold:
Setting the convergence threshold
---------------------------------
By default the convergence threshold is set to 1e-8. To modify this value,
adjust the :py:class:`~pybest.wrapper.hf.RHF` or :py:class:`~pybest.wrapper.hf.UHF`
function call by adding the ``threshold`` keyword argument as follows,
.. code:: python
# set a tighter convergence threshold in RHF
hf_ = hf(kin, na, er, external, olp, orb_a, threshold=1e-10)
The maximum number of SCF iterations steps (default 128) can be adjusted using
the ``maxiter`` keyword,
.. code:: python
# set a tighter convergence threshold and increase number of iterations in RHF
hf_ = hf(kin, na, er, external, olp, orb_a, threshold=1e-10, maxiter=200)
.. _user_scf_convergence_diis:
Choosing the DIIS solver
------------------------
PyBEST supports four different SCF algorithms:
* A simple SCF solver without any convergence acceleration: ``plain`` or ``None``
* The traditional SCF solver employing Pulay mixing: ``cdiis`` (default)
* An SCF solver using energy-based DIIS: ``ediis``
* An SCF solver combing Pulay mixing and energy-based DIIS: ``ediis2``
.. * An optimal damping SCF solver: ``damping``
If ``plain`` or ``None`` is selected, an ordinary SCF solver is chosen that
does not feature any convergence acceleration techniques. Instead, in each
iteration, the Fock matrix is built and then diagonalized.
By default, the ``cdiis`` options is selected, which corresponds to the (traditional)
commutator-based direct inversion of the iterative subspace (CDIIS) algorithm,
which is also known as Pulay mixing [pulay1980]_. This SCF flavour is typically
very efficient, albeit sensitive to the initial guess.
``ediis`` invokes the energy direct inversion of the iterative subspace (EDIIS) method [kudin2002]_.
This DIIS flavour works well for the first few iterations, becomes, however,
numerically unstable close to convergence. It is the preferred SCF algorithm
in cases the initial guess is poor.
Finally, the ``ediis2`` method combines the CDIIS and EDIIS algorithms [kudin2002]_.
This hybrid method aims at joining the benefits of both approaches.
The DIIS flavor can be changed using the keyword ``diis`` in the function call
(similar for both *restricted* and *unrestricted* orbitals),
.. code:: python
# select plain solver
hf_ = hf(kin, na, er, external, olp, orb_a, diis='plain') # or diis=None
# select CDIIS (default)
hf_ = hf(kin, na, er, external, olp, orb_a, diis='cdiis')
# select EDIIS
hf_ = hf(kin, na, er, external, olp, orb_a, diis='ediis')
# select EDIIS2
hf_ = hf(kin, na, er, external, olp, orb_a, diis='ediis2')
The behavior of the DIIS algorithm can be adjusted using additional keywords.
For instance the convergence threshold and number of iterations can be controlled
as described above (:ref:`user_scf_convergence_threshold`). In addition,
the number of DIIS vectors can be adjusted (default value is 6).
To do so, use the keyword ``nvector``,
.. code:: python
# select CDIIS (default) with at most 8 vectors
hf_ = hf(kin, na, er, external, olp, orb_a, nvector=8)
.. _user_scf_to_file:
Dumping SCF results
===================
At each SCF iteration step, the current results (orbitals, densities, etc.) are
dumped to disk and stored for later restarts, visualization, or post-processing.
If converged, those results are updated one last time and stored in the
checkpoint file ``checkpoint_scf.h5`` of the ``pybest-results`` directory (by default).
After SCF convergence, you can also dump all Hartree-Fock results to user-specified
file and file format. Currently, PyBEST supports the following file formats
1. PyBEST's internal format based on HDF5
2. the Molden format: :ref:`user_orbitals_molden`
The internal format has the advantage that as much information as provided by
the user can be stored in full precision.
For internal storage, PyBEST uses the :py:class:`~pybest.io.iodata.IOData`
container (see :ref:`user_iodata` for more details).
For visualization purposes, the Molden format is preferred. Furthermore,
the Molden format can be read by other quantum
chemistry codes. Unfortunately, Molden only supports basis functions up to
including G functions.
Since the return value of the :py:class:`~pybest.wrapper.hf.RHF` and
:py:class:`~pybest.wrapper.hf.UHF` function call is an
:py:class:`~pybest.io.iodata.IOData` container, dumping the corresponding
(RHF or UHF) results to different file formats is particularly easy.
To dump results to disk the :py:meth:`~pybest.io.iodata.IOData.to_file` method
of the :py:class:`~pybest.io.iodata.IOData` container has to be used, where
the file ending specifies the file format.
PyBEST's internal format
------------------------
To dump into PyBEST's internal file format, use the extension ``.h5``
.. code-block:: python
# Dump data to file using PyBEST's internal format
hf_.to_file('water-scf.h5')
The Molden format
-----------------
To dump Molden files, we have to update the :py:class:`~pybest.io.iodata.IOData`
container to include also the ``gobasis`` information (by default, it is
omitted from the container)
.. literalinclude:: ../src/pybest/data/examples/hf/rhf_water_dense.py
:lines: 56-60
See also :ref:`user_iodata_dumping_filling` on how to modify an
:py:class:`~pybest.io.iodata.IOData` container on the fly.
.. _user_scf_restart:
Restarting SCF calculations
===========================
Restarting from a checkpoint file
---------------------------------
One straightforward way to restart a Hartree-Fock calculation from a PyBEST
checkpoint file is to use the ``restart`` keyword, which specifies the file name
of the previously dumped SCF solution.
PyBEST will automatically read in the new orbitals and perform a projection.
.. literalinclude:: ../src/pybest/data/examples/hf/rhf_water_dense.py
:lines: 62-66
.. _user_scf_restart_perturbed:
Restarting from perturbed orbitals
----------------------------------
In case you want to perturb the orbitals prior restart, like swapping their order
(:ref:`orbitals_swapping`) or 2x2 rotations (:ref:`orbitals_rotations`), you
have to read in the corresponding checkpoint file form disk manually, perform
the orbital manipulation you desire, and then pass the modified orbitals to the
:py:class:`~pybest.wrapper.hf.RHF` or :py:class:`~pybest.wrapper.hf.UHF`
function call,
.. literalinclude:: ../src/pybest/data/examples/hf/rhf_water_dense.py
:lines: 68-
.. note::
In order to guarantee that the :py:class:`~pybest.wrapper.hf.RHF` or
:py:class:`~pybest.wrapper.hf.UHF` solvers swap the orbitals in question,
``skip_occs`` **and** ``skip_energies`` have to be set to ``True``.
Otherwise, orbital swapping has **no** effect at all.
.. _user_scf_examples:
Example Python scripts
======================
The following, we list a basic example script for a restricted Hartree-Fock
calculation for the water molecule.
.. literalinclude:: ../src/pybest/data/examples/hf/rhf_water_dense.py
:caption: data/examples/hf/rhf_water_dense.py
:lines: 8-