..
: 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/>
:
: --
Molecular Hamiltonians
######################
A molecular Hamiltonian is defined in three steps:
1. Set up a molecular geometry (see :ref:`setup-molgeometry`)
2. Define a Gaussian basis set (see :ref:`user_molecularham_basis`)
3. Calculate the matrix representation of all contributions separately that
are defined in the molecular Hamiltonian (see :ref:`user_molecularham_matrix_elements`)
Before a HF or post-HF calculation can be performed, the overlap matrix
of the atomic basis set as well as the molecular orbitals need to be specified
(see :ref:`user_molecularham_overlap`)
.. _setup-molgeometry:
Specifying the molecular geometry
=================================
The basis set reader in PyBEST requires the molecular coordinates to be defined
in an xyz file. If an xyz file is present, only the (relative) path to the file
(in most cases just the file name) needs to be specified. If the molecular
geometry is not present as an xyz file, PyBEST will internally convert it to
a file in the xyz format and store it in a temporary directory.
Note that PyBEST can load/dump the molecular geometry from/to different file
formats. The file format is determined automatically using the extension or
prefix of the filename. For more information on supported file formats, refer to
:ref:`user_iodata`.
In addition to providing the filename, the molecular geometry can be specified
as described below.
.. _read-molgeometry:
Reading the molecular geometry from file
----------------------------------------
The molecular geometry can be read from file using the method
:py:meth:`~pybest.io.iodata.IOData.from_file` of the
:py:class:`~pybest.io.iodata.IOData` class,
.. code-block:: python
mol = IOData.from_file(filename)
For instance, the molecular coordinates can be provided in the conventional
xyz file format,
.. code-block:: python
mol = IOData.from_file('molecule.xyz')
If just the filename is given, PyBEST assumes that the file is present in the
working directory. If the coordinate file is lying elsewhere, the absolute or
relative path have to be provided, for instance,
.. code-block:: python
mol = IOData.from_file('./coords/molecule.xyz')
Constructing a molecular geometry from scratch
----------------------------------------------
A molecular geometry can also be generated using Python code. The
following example constructs the atomic coordinates for a water molecule
and writes the corresponding geometry information to an xyz file.
Note that the coordinates have to be defined in bohr (atomic units).
The xyz file is represented in Angstrom.
.. literalinclude:: ../src/pybest/data/examples/hamiltonian/water.py
:lines: 3-
:caption: data/examples/hamiltonian/water.py
.. _user_molecularham_basis:
Specifying the atomic basis set
===============================
PyBEST uses the Libint2 basis set reader. Thus, the supported basis set format
is the **g94** format. Note that Libint2 needs to be compiled to support up to
a maximum angular momentum quantum number of seven, that is, I functions.
PyBEST is distributed with its own basis set library. An atomic basis set can
be read in from the PyBEST basis set library using the function call
:py:func:`~pybest.gbasis.gobasis.get_gobasis`.
.. _bs-library:
Supported atomic basis sets in the PyBEST basis set library
------------------------------------------------------------
All supported basis sets are located in ``src/pybest/data/basis`` and include
1. Minimal basis sets
* STO-3G, STO-6G
2. Pople basis sets
* 3-21g
* 6-31g, 6-31g*, 6-31g**, 6-31++g**
3. Correlation consistent basis sets
* cc-pVDZ, cc-pVTZ, cc-pVQZ, cc-pV5Z
* cc-pcVDZ, cc-pcVTZ, cc-pcVQZ
* aug-cc-pVDZ, aug-cc-pVTZ, aug-cc-pVQZ, aug-cc-pV5Z
4. ANO basis sets
* ANO-RCC (large)
* ANO-RCC-VDZ, ANO-RCC-VDZP, ANO-RCC-VTZ, ANO-RCC-VTZP, ANO-RCC-VQZP
The same basis set for different atoms
--------------------------------------
By default, all atoms in the molecule are assigned the same atomic basis set.
To define an atomic basis set, the molecular geometry has to be provided to
the :py:func:`~pybest.gbasis.gobasis.get_gobasis` (see :ref:`read-molgeometry`).
The atomic basis set can then be constructed by only one function call,
.. code-block:: python
gobasis = get_gobasis(basisname, coordinate_filename)
where ``basisname`` is the name of the basis set (a string) and ``coordinate_filename``
is the file containing the coordinates. If a ``cc-pVDZ`` basis is used and the
coordinates are stored in the ``mol.xyz`` file, the basis set can be defined
as follows
.. code-block:: python
gobasis = get_gobasis('cc-pvdz', 'mol.xyz')
Loading user-defined basis sets from file
-----------------------------------------
If an atomic basis set is not present in the PyBEST basis set library, it can
be read from file. Note that the basis set file has to be available in the
g94 file format.
To read in a user-defined basis set file, the ``basisname`` argument of the
:py:func:`~pybest.gbasis.gobasis.get_gobasis` function has to be substituted
by the path to the basis set file.
.. code-block:: python
gobasis = get_gobasis(path_to_basis, coordinate_filename)
``path_to_basis`` can be either the absolute or relative path. If the user-specified
basis set is stored in the ``my_own_dz.g94`` file, the above example becomes
.. code-block:: python
gobasis = get_gobasis('./my_own_dz.g94', 'mol.xyz')
.. note::
The file extension ``.g94`` has to be explicitely provided.
Only with the file extension, PyBEST will read a user-defined basis set
instead of searching for the basis set in its basis set library.
.. _user_molecularham_geom_and_basis:
Loading geometry and basis set information from another file or external program
--------------------------------------------------------------------------------
The molecular geometry and basis set information can also be read in from
another file generated by external software like (Molpro, ORCA, PSI4, etc).
Supported file formats are ``.mkl`` and ``.molden`` (in addition to
PyBEST's internal ``.h5`` format).
As for molecular geometries, the :py:meth:`~pybest.io.iodata.IOData.from_file`
method has to be used to load the file.
.. code-block:: python
# Load the geometry, orbital basis, and wave function from a molden file:
mol = IOData.from_file('water.molden')
# Print the number of basis functions
print(mol.gobasis.nbasis)
The corresponding information on the basis set or molecular orbitals is stored
in the corresponding attributes. For instance, the basis set object is available
in the ``gobasis`` attribute as the return value.
See also :ref:`user_iodata` for more details on how to handle I/O operations.
.. _setup-dummy:
Specifying dummy or ghost atoms and active fragments
====================================================
In PyBEST, dummy or ghost atoms and active molecular fragments are defined
using the :py:func:`~pybest.gbasis.gobasis.get_gobasis` function, that is,
during the construction of the atomic basis set.
Dummy or ghost atoms
--------------------
All dummy or ghost atoms are passed as a list to the :py:func:`~pybest.gbasis.gobasis.get_gobasis`
function that contains the indices of all dummy/ghost atoms as its elements.
The elements are sorted according to their order in the xyz coordinate file.
.. code-block:: python
# set first (0) and second (1) element in mol.xyz as dummy/ghost atom
gobasis = get_gobasis('cc-pvdz', 'mol.xyz', dummy=[0,1])
.. note::
In PyBEST, all indexing is performed using the Python convention, that is,
beginning with 0.
Active fragments
----------------
PyBEST allows us to divide one xyz coordinate file in an active and inactive
fragment. The inactive fragment will be neglected during the construction
of the basis set and molecular Hamiltonian.
.. code-block:: python
# define molecule containing atoms 3, 4, and 5 as active fragment contained
# in the mol.xyz file
gobasis = get_gobasis('cc-pvdz', 'mol.xyz', active_fragment=[3,4,5])
Combining ghost atoms and active fragments
------------------------------------------
The ``dummy`` option and the ``active_fragment`` option can be combined so that additional ghost or
dummy atoms can be defined in the active fragment.
.. code-block:: python
# define molecule containing atoms 3, 4, and 5 as active fragment contained
# in the mol.xyz file, set atom index 3 to be a ghost atom
gobasis = get_gobasis('cc-pvdz', 'mol.xyz', active_fragment=[3,4,5], dummy=[3])
.. note::
The dummy index has to be included in the ``active_fragment`` list. If this is not the case,
an exception will be raised.
.. _user_molecularham_matrix_elements:
Computing the matrix representation of the Hamiltonian
======================================================
As soon as the molecular geometry and atomic basis sets are defined, the matrix
representation of the molecular Hamiltonian can be calculated.
In the current version of PyBEST, each term has to be evaluated separately.
The molecular Hamiltonian contains only one- and two-electron integrals as well
as a constant term due to the nuclei-nuclei repulsion.
Various one- and two-electron integrals are implemented in PyBEST (see also
:ref:`hamiltonian_list_of_ints` for a list of supported integrals).
Most integrals are calculated using the LibInt package [valeev2019]_ using the
modern C++ API, except the Cholesky-decomposed two-electron repulsion integrals,
which are evaluated using the LibChol Library, the pVp integrals, which are
determined by interfacing specific features from the LibInt package [valeev2019]_,
and the DKH2 Hamiltonian, which is calculated using Python code.
All one-electron integrals as well as the overlap matrix of the atomic basis
set are stored as two-index objects (that is, in its full or dense matrix
representation), which are implemented in the ``DenseLinalgFactory`` class.
The two-electron integrals are either stored as a dense four-index object
(``DenseLinalgFactory`` class, memory intensive) or Cholesky decomposed
(``CholeskyLinalgFactory`` class, more memory efficient).
See :ref:`user_linalg` for more information on the ``LinalgFactory`` and its features.
The evaluation of all one- and two-electron integrals is performed by
various PyBEST function calls.
.. _user_schrodinger_hamiltonian:
The Schrödinger Hamiltonian
---------------------------
The electronic Schrödigner Hamiltonian,
.. math::
:label: schrodinger
{\hat{H}}_{\mathrm{el}} = \sum_i^{N_{\mathrm{el}}} \hat{T}_i +
\sum_i^{N_{\mathrm{el}}} \sum_I^{N_{\mathrm{nuc}}} \hat{V}_{iI} +
\frac{1}{2} \sum_{i,j}^{N_{\mathrm{el}}} \hat{V}_{i,j} +
\frac{1}{2} \sum_{I,J}^{N_{\mathrm{nuc}}} \hat{V}_{I,J}
which contains the kinetic energy of the electrons, the electron-nucleus
attraction, the electron-electron repulsion, and the nucleus-nucleus repulsion
terms.
In PyBEST, each term has to be evaluated explicitly.
Before computing all one- and two-electron integrals, the Linalg factory
has to be specified.
Currently, two Linalg factories are supported: ``DenseLinalgFactory`` and
``CholeskyLinalgFactory``, which differ only in the representation of the
two-electron repulsion integrals. The one-electron integrals are represented
in their dense representation using the ``DenseLinalgFactory`` class,
irrespective of the choice between ``DenseLinalgFactory`` or
``CholeskyLinalgFactory``.
For both cases, the number of atomic basis functions has to be passed as an
argument. It is stored in the ``nbasis`` attribute of the ``PyBasis`` instance
(here ``gobasis``).
.. code-block:: python
# use DenseLinalgFactory
lf = DenseLinalgFactory(gobasis.nbasis)
# use CholeskyLinalgFactory
lf = CholeskyLinalgFactory(gobasis.nbasis)
.. note::
The ``DenseLinalgFactory`` and ``CholeskyLinalgFactory`` differ only in their
representation of the two-electron integrals. All smaller dimensional tensors are
represented using the DenseLinalgFactory.
The kinetic energy can be calculated using the
:py:meth:`~pybest.gbasis.gobasis.compute_kinetic` function, which takes an
instance of ``PyBasis`` as argument,
.. code-block:: python
# evaluate kinetic energy of electrons and store in kin
kin = compute_kinetic(gobasis)
Similarly, the electron repulsion integrals are determined
.. code-block:: python
# evaluate electron repulsion integrals and store in eri
eri = compute_eri(gobasis)
# or evaluate Cholesky-decomposed two-electron integals, the default threshold is
# set to 1e-3
eri = compute_cholesky_eri(gobasis, threshold=1e-8)
The potential energy associated with the nuclei is determined in a similar
fashion,
.. code-block:: python
# Nuclear-electron attraction term stored in na
na = compute_nuclear(gobasis)
# Nuclear-nuclear repulsion term stored as dictionary in nuc
nuc = compute_nuclear_repulsion(gobasis)
A complete example can be found here:
.. literalinclude:: ../src/pybest/data/examples/hamiltonian/schrodinger.py
:lines: 3-38
:caption: data/examples/hamiltonian/schrodinger.py, lines 3-39
.. _user_dkh_hamiltonian:
The DKH Hamiltonian
-------------------
The one-electron terms in eq. :eq:`schrodinger` can be substituted by the
DKH Hamiltonian to account for (scalar) relativistic effects. PyBEST currently
supports DKH up to second order (DKH2).
The DKH2 Hamiltonian can be calculated using the
:py:meth:`~pybest.gbasis.dkhn.compute_dkh` function, which takes an
instance of ``PyBasis`` as argument,
.. code-block:: python
# evaluate the DKH2 Hamiltonian and store in dkh
dkh = compute_dkh(gobasis)
A complete example can be found here:
.. literalinclude:: ../src/pybest/data/examples/hamiltonian/dkh.py
:lines: 3-32
:caption: data/examples/hamiltonian/dkh.py, lines 3-32
.. _user_molecularham_overlap:
Computing the overlap matrix and defining the orbitals
======================================================
To perform a HF or post-HF calculation, the overlap matrix of the atomic
orbital basis set has to be calculated and a set of (molecular) orbitals has
to be defined.
The overlap matrix can be calculated using the :py:meth:`~pybest.gbasis.gobasis.compute_overlap`
method, which requires an instance of ``PyBasis`` as argument. The return value is
a dense two-index object of ``DenseLinalgFactory`` (as all one-electron integrals).
.. code-block:: python
# overlap matrix of the atomic basis set
olp = compute_overlap(gobasis)
The (molecular) orbitals are initialized using the :py:meth:`~pybest.linalg.dense.create_orbital`
method of the ``DenseLinalgFactory`` class,
.. code-block:: python
# create molecular orbitals
orb = lf.create_orbital()
The orbitals are an instance of the ``DenseOrbital`` class. The instance
``orb`` contains the AO/MO coefficients (``orb.coeffs``), the occupation numbers
(``orb.occupations``), and orbital energies (``orb.energies``),
which are stored as two- (coefficients) or one-dimensional (occupations, energies)
numpy arrays. See also :ref:`user_wroking_with_orbitals` for more details on
orbitals and how to work with orbitals in PyBEST.
.. _user_molecularham_dipole:
Computing the multipole moment integrals
========================================
PyBEST supports the calculation of the electric dipole and quadrupole moment
integrals.
The electric dipole integrals can be calculated using the
:py:meth:`~pybest.gbasis.gobasis.compute_dipole` method, which requires an
instance of ``PyBasis`` as argument. The return value is a list containing the
overlap matrix and each component {x, y, z} of the (Cartesian) electric dipole operator,
each being a dense two-index object of ``DenseLinalgFactory`` (as all one-electron integrals).
.. code-block:: python
# electric dipole moment integrals of the atomic basis set
olp, mux, muy, muz = compute_dipole(gobasis)
# or store as tuple
dipole = compute_dipole(gobasis)
By default, the multipole origin is set to ``{0,0,0}``. This can be adjusted by
setting the ``x``, ``y``, and ``z`` arguments to some new origin.
If you wish to take the center of mass (COM) of the molecule as the new origin,
you first determine the COM and then pass it explicitly to the
:py:meth:`~pybest.gbasis.gobasis.compute_dipole` method,
.. code-block:: python
# determine COM
x, y, z = get_com(gobasis)
# electric dipole moment integrals of the atomic basis set wrt COM
olp, mux, muy, muz = compute_dipole(gobasis, x=x, y=y, z=z)
# or store as tuple
dipole = compute_dipole(gobasis, x=x, y=y, z=z)
The electric quadrupole moment integrals can be determined in a similar way
using the :py:meth:`~pybest.gbasis.gobasis.compute_quadrupole` method.
In addition to the overlap matrix and each component of the (Cartesian) electric
dipole operator, all 6 components, {xx, xy, xz, yy, yz, zz}, of the Cartesian
electric quadrupole operator are returned (thus, 10 objects in total),
.. code-block:: python
# electric quadrupole moment integrals of the atomic basis set
# ------------------------------------------------------------
olp, mux, muy, muz, muxx, muxy, muxz, muyy, muyz, muzz = compute_quadrupole(gobasis)
# or store as tuple
quadrupole = compute_quadrupole(gobasis)
Similar to the electric dipole moment integrals, the quadrupole moment integrals
are determined with respect to the origin of the Cartesian coordinate system,
``{0,0,0}``. This can be adjusted by setting the ``x``, ``y``, and ``z`` arguments
to some new origin,
.. code-block:: python
# determine COM
x, y, z = get_com(gobasis)
# electric quadrupole moment integrals of the atomic basis set wrt COM
# --------------------------------------------------------------------
olp, mux, muy, muz, muxx, muxy, muxz, muyy, muyz, muzz = compute_quadrupole(gobasis, x=x, y=y, z=z)
# or store as tuple
quadrupole = compute_quadrupole(gobasis, x=x, y=y, z=z)
PyBEST provides additional wrappers that can be used to determine the electric
dipole moment (see :ref:`dipole` for more details).
Example Python script
=====================
The molecular Hamiltonian for the water molecule
------------------------------------------------
This example shows how to set up the molecular Hamiltonian, orbitals, and
overlap matrix for the water molecule.
.. literalinclude:: ../src/pybest/data/examples/hamiltonian/schrodinger.py
:lines: 3-
:caption: data/examples/hamiltonian/schrodinger.py, lines 3-ff