Source code for iodata.formats.gaussianlog

# IODATA is an input and output module for quantum chemistry.
# Copyright (C) 2011-2019 The IODATA Development Team
#
# This file is part of IODATA.
#
# IODATA 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.
#
# IODATA 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/>
# --
"""Gaussian Log file format.

To write out the integrals in a Gaussian log file, which can be loaded with this module,
you need to use the following Gaussian command line:

.. code-block ::

    scf(conventional) iop(3/33=5) extralinks=l316 iop(3/27=999)

"""

import numpy as np
from numpy.typing import NDArray

from ..docstrings import document_load_one
from ..utils import LineIterator, set_four_index_element

__all__ = ()


PATTERNS = ["*.log"]


LOAD_ONE_NOTES = """\
Note that Gaussian writes four-center integrals in chemists' notation.
They will be reordered into physicists' notation when loading them with IOData.
"""




@document_load_one("Gaussian Log", [], ["one_ints", "two_ints"], notes=LOAD_ONE_NOTES)
def load_one(lit: LineIterator) -> dict:
    """Do not edit this docstring. It will be overwritten."""
    # First get the line with the number of orbital basis functions
    for line in lit:
        if line.startswith("    NBasis ="):
            nbasis = int(line[12:18])
            break

    # Then load the two- and four-index operators. This part is written such
    # that it does not make any assumptions about the order in which these
    # operators are printed.
    one_ints = {}
    two_ints = {}
    while True:
        line = next(lit)
        if line.startswith(" Normal termination of Gaussian"):
            break
        if line.startswith(" *** Overlap ***"):
            one_ints["olp"] = _load_twoindex_g09(lit, nbasis)
        elif line.startswith(" *** Kinetic Energy ***"):
            one_ints["kin_ao"] = _load_twoindex_g09(lit, nbasis)
        elif line.startswith(" ***** Potential Energy *****"):
            one_ints["na_ao"] = _load_twoindex_g09(lit, nbasis)
        elif line.startswith(" *** Dumping Two-Electron integrals ***"):
            two_ints["er_ao"] = _load_fourindex_g09(lit, nbasis)

    result = {}
    if one_ints:
        result["one_ints"] = one_ints
    if two_ints:
        result["two_ints"] = two_ints
    return result



def _load_twoindex_g09(lit: LineIterator, nbasis: int) -> NDArray[float]:
    """Load a two-index operator from a GAUSSIAN LOG file format.

    Parameters
    ----------
    lit
        The line iterator to read the data from.
    nbasis
        The number of atomic orbital basis functions.

    Returns
    -------
    Operator array with shape ``(nbasis, nbasis)``.

    """
    result = np.zeros((nbasis, nbasis))
    block_counter = 0
    while block_counter < nbasis:
        # skip the header line
        next(lit)
        # determine the number of rows in this part
        nrow = nbasis - block_counter
        for i in range(nrow):
            words = next(lit).split()[1:]
            for j, word in enumerate(words):
                value = float(word.replace("D", "E"))
                result[i + block_counter, j + block_counter] = value
                result[j + block_counter, i + block_counter] = value
        block_counter += 5
    return result


def _load_fourindex_g09(lit: LineIterator, nbasis: int) -> NDArray[float]:
    """Load a four-index operator from a GAUSSIAN LOG file.

    Parameters
    ----------
    lit
        The line iterator to read the data from.
    nbasis
        The number of atomic orbital basis functions.

    Returns
    -------
    Operator array with shape ``(nbasis, nbasis, nbasis, nbasis)``.

    """
    result = np.zeros((nbasis, nbasis, nbasis, nbasis))
    # Skip first six lines
    for _i in range(6):
        next(lit)
    # Start reading elements until a line is encountered that does not start
    # with ' I='
    for line in lit:
        if not line.startswith(" I="):
            break
        # print line[3:7], line[9:13], line[15:19], line[21:25], line[28:].replace('D', 'E')
        i0 = int(line[3:7]) - 1
        i1 = int(line[9:13]) - 1
        i2 = int(line[15:19]) - 1
        i3 = int(line[21:25]) - 1
        value = float(line[29:].replace("D", "E"))
        # Gaussian uses the chemists' notation for the 4-center indexes. IOdata
        # uses the physicists' notation.
        set_four_index_element(result, i0, i2, i1, i3, value)
    return result