# 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/>
# --
"""GAMESS punch file format."""
import numpy as np
from numpy.typing import NDArray
from ..docstrings import document_load_one
from ..utils import LineIterator, LoadError, angstrom
__all__ = ()
PATTERNS = ["*.dat"]
def _read_data(lit: LineIterator) -> tuple[str, str, list[str]]:
"""Extract ``title``, ``symmetry`` and ``symbols`` from the punch file."""
title = next(lit).strip()
symmetry = next(lit).split()[0]
# The dat file only contains symmetry-unique atoms, so we would be incapable of
# supporting non-C1 symmetry without significant additional coding.
if symmetry != "C1":
raise LoadError(f"Only C1 symmetry is supported, got {symmetry}.", lit)
symbols = []
line = True
while line != " $END \n":
line = next(lit)
if line[0] != " ":
symbols.append(line.split()[0])
return title, symmetry, symbols
def _read_coordinates(lit: LineIterator, result: dict[str]) -> tuple[NDArray[int], NDArray[float]]:
"""Extract ``numbers`` and ``coordinates`` from the punch file."""
for _ in range(2):
next(lit)
natom = len(result["symbols"])
# if the data are already read before, just overwrite them
numbers = result.get("atnums")
if numbers is None:
numbers = np.zeros(natom, int)
result["atnums"] = np.zeros(natom, int)
coordinates = result.get("atcoords")
if coordinates is None:
coordinates = np.zeros((natom, 3), float)
for i in range(natom):
words = next(lit).split()
numbers[i] = int(float(words[1]))
coordinates[i] = np.array([float(elem) for elem in words[2:5]]) * angstrom
return numbers, coordinates
def _read_energy(lit: LineIterator, result: dict[str]) -> tuple[float, NDArray[float]]:
"""Extract ``energy`` and ``gradient`` from the punch file."""
energy = float(next(lit).split()[1])
natom = len(result["symbols"])
# if the data are already read before, just overwrite them
gradient = result.get("gradient")
if gradient is None:
gradient = np.zeros((natom, 3), float)
for i in range(natom):
words = next(lit).split()
gradient[i] = words[2:5]
return energy, gradient
def _read_hessian(lit: LineIterator, result: dict[str]) -> NDArray[float]:
"""Extract ``hessian`` from the punch file."""
# check that $HESS is not already parsed
if "athessian" in result:
raise LoadError(
"Cannot parse $HESS twice. Make sure approximate hessian is not being parsed."
)
next(lit)
natom = len(result["symbols"])
hessian = np.zeros((3 * natom, 3 * natom), float)
tmp = hessian.ravel()
counter = 0
while True:
line = next(lit)
if line == " $END\n":
break
line = line[5:-1]
for j in range(len(line) // 15):
tmp[counter] = float(line[j * 15 : (j + 1) * 15])
counter += 1
return hessian
def _read_masses(lit: LineIterator, result: dict[str]) -> NDArray[float]:
"""Extract ``masses`` from the punch file."""
natom = len(result["symbols"])
masses = np.zeros(natom, float)
counter = 0
while counter < natom:
words = next(lit).split()
for word in words:
masses[counter] = float(word)
counter += 1
return masses
[docs]
@document_load_one(
"PUNCH",
["title", "energy", "grot", "atgradient", "athessian", "atmasses", "atnums", "atcoords"],
)
def load_one(lit: LineIterator) -> dict[str]:
"""Do not edit this docstring. It will be overwritten."""
result = {}
while True:
try:
line = next(lit)
except StopIteration:
break
if line == "$DATA\n":
result["title"], result["g_rot"], result["symbols"] = _read_data(lit)
elif line == " COORDINATES OF SYMMETRY UNIQUE ATOMS (ANGS)\n":
result["atnums"], result["atcoords"] = _read_coordinates(lit, result)
elif line == " $GRAD\n":
result["energy"], result["atgradient"] = _read_energy(lit, result)
elif line == "CAUTION, APPROXIMATE HESSIAN!\n":
# prevent the approximate Hessian from being parsed
while line != " $END\n":
line = next(lit)
elif line == " $HESS\n":
result["athessian"] = _read_hessian(lit, result)
elif line == "ATOMIC MASSES\n":
result["atmasses"] = _read_masses(lit, result)
result.pop("symbols")
return result