"""Module for GAMESS US IO."""
from __future__ import annotations
import os
import re
from copy import deepcopy
from io import TextIOBase
from subprocess import TimeoutExpired, call
import numpy as np
from ase import Atoms
from ase.calculators.singlepoint import SinglePointCalculator
from ase.units import Bohr, Debye, Hartree
from ase.utils import reader, workdir, writer
def _format_value(val: bool | str) -> str:
if isinstance(val, bool):
return '.t.' if val else '.f.'
return str(val).upper()
def _write_block(name: str, args: dict[str, bool | str]) -> str:
out = [f' ${name.upper()}']
for key, val in args.items():
out.append(f' {key.upper()}={_format_value(val)}')
out.append(' $END')
return '\n'.join(out)
def _write_geom(atoms: Atoms, basis_spec: dict[str | int, str] | None) -> str:
out = [' $DATA', atoms.get_chemical_formula(), 'C1']
for i, atom in enumerate(atoms):
out.append(
'{:<3} {:>3} {:20.13e} {:20.13e} {:20.13e}'.format(
atom.symbol, atom.number, *atom.position
)
)
if basis_spec is not None:
basis = basis_spec.get(i)
if basis is None:
basis = basis_spec.get(atom.symbol)
if basis is None:
msg = f'Could not find an appropriate basis set for atom {i}'
raise ValueError(msg)
out += [basis, '']
out.append(' $END')
return '\n'.join(out)
def _write_ecp(atoms: Atoms, ecp: dict[int | str, str]) -> str:
out = [' $ECP']
for i, symbol in enumerate(atoms.symbols):
if i in ecp:
out.append(ecp[i])
elif symbol in ecp:
out.append(ecp[symbol])
else:
msg = f'Could not find an appropriate ECP for atom {i}'
raise ValueError(msg)
out.append(' $END')
return '\n'.join(out)
_xc = dict(LDA='SVWN')
[docs]
@writer
def write_gamess_us_in(fd: TextIOBase, atoms: Atoms, properties=None, **params):
params = deepcopy(params)
if properties is None:
properties = ['energy']
# set RUNTYP from properties iff value not provided by the user
contrl = params.pop('contrl', {})
if 'runtyp' not in contrl:
if 'forces' in properties:
contrl['runtyp'] = 'gradient'
else:
contrl['runtyp'] = 'energy'
# Support xc keyword for functional specification
xc = params.pop('xc', None)
if xc is not None and 'dfttyp' not in contrl:
contrl['dfttyp'] = _xc.get(xc.upper(), xc.upper())
# Automatically determine multiplicity from magnetic moment
magmom_tot = int(round(atoms.get_initial_magnetic_moments().sum()))
if 'mult' not in contrl:
contrl['mult'] = abs(magmom_tot) + 1
# Since we're automatically determining multiplicity, we also
# need to automatically switch to UHF when the multiplicity
# is not 1
if 'scftyp' not in contrl:
contrl['scftyp'] = 'rhf' if contrl['mult'] == 1 else 'uhf'
# effective core potentials
ecp = params.pop('ecp', None)
if ecp is not None and 'pp' not in contrl:
contrl['pp'] = 'READ'
# If no basis set is provided, use 3-21G by default.
basis_spec = None
if 'basis' not in params:
params['basis'] = dict(gbasis='N21', ngauss=3)
else:
keys = set(params['basis'])
# Check if the user is specifying a literal per-atom basis.
# We assume they are passing a per-atom basis if the keys of the
# basis dict are atom symbols, or if they are atom indices, or
# a mixture of both.
if keys.intersection(set(atoms.symbols)) or any(
map(lambda x: isinstance(x, int), keys)
):
basis_spec = params.pop('basis')
out = [_write_block('contrl', contrl)]
out += [_write_block(*item) for item in params.items()]
out.append(_write_geom(atoms, basis_spec))
if ecp is not None:
out.append(_write_ecp(atoms, ecp))
fd.write('\n\n'.join(out))
_geom_re = re.compile(r'^\s*ATOM\s+ATOMIC\s+COORDINATES')
_atom_re = re.compile(r'^\s*(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s*\n')
_energy_re = re.compile(r'^\s*FINAL [\S\s]+ ENERGY IS\s+(\S+) AFTER')
_grad_re = re.compile(r'^\s*GRADIENT OF THE ENERGY\s*')
_charges_re = re.compile(r'^\s*TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS\s*')
_dipole_re = re.compile(r'^\s+DX\s+DY\s+DZ\s+\/D\/\s+\(DEBYE\)')
[docs]
@reader
def read_gamess_us_out(fd):
atoms = None
energy = None
forces = None
charges = None
dipole = None
for line in fd:
# Geometry
if _geom_re.match(line):
fd.readline()
symbols = []
pos = []
while True:
atom = _atom_re.match(fd.readline())
if atom is None:
break
symbol, _, x, y, z = atom.groups()
symbols.append(symbol.capitalize())
pos.append(list(map(float, [x, y, z])))
atoms = Atoms(symbols, np.array(pos) * Bohr)
continue
# Energy
ematch = _energy_re.match(line)
if ematch is not None:
energy = float(ematch.group(1)) * Hartree
# MPn energy. Supplants energy parsed above.
elif line.strip().startswith('TOTAL ENERGY'):
energy = float(line.strip().split()[-1]) * Hartree
# Higher-level energy (e.g. coupled cluster)
# Supplants energies parsed above.
elif line.strip().startswith('THE FOLLOWING METHOD AND ENERGY'):
energy = float(fd.readline().strip().split()[-1]) * Hartree
elif _charges_re.match(line):
fd.readline()
charges = np.array(
[float(fd.readline().split()[3]) for _ in symbols],
)
# Gradients
elif _grad_re.match(line):
for _ in range(3):
fd.readline()
grad = []
while True:
atom = _atom_re.match(fd.readline())
if atom is None:
break
grad.append(list(map(float, atom.groups()[2:])))
forces = -np.array(grad) * Hartree / Bohr
elif _dipole_re.match(line):
dipole = np.array(list(map(float, fd.readline().split()[:3])))
dipole *= Debye
atoms.calc = SinglePointCalculator(
atoms,
energy=energy,
forces=forces,
charges=charges,
dipole=dipole,
)
return atoms
[docs]
@reader
def read_gamess_us_punch(fd):
atoms = None
energy = None
forces = None
charges = None
dipole = None
for line in fd:
if line.strip() == '$DATA':
symbols = []
pos = []
while line.strip() != '$END':
line = fd.readline()
atom = _atom_re.match(line)
if atom is None:
# The basis set specification is interlaced with the
# molecular geometry. We don't care about the basis
# set, so ignore lines that don't match the pattern.
continue
symbols.append(atom.group(1).capitalize())
pos.append(list(map(float, atom.group(3, 4, 5))))
atoms = Atoms(symbols, np.array(pos))
elif line.startswith('E('):
energy = float(line.split()[1][:-1]) * Hartree
elif line.strip().startswith('POPULATION ANALYSIS'):
# Mulliken charges
charges = np.array(
[float(fd.readline().split()[2]) for _ in symbols],
)
elif line.strip().startswith('DIPOLE'):
dipole = np.array(list(map(float, line.split()[1:]))) * Debye
elif line.strip() == '$GRAD':
# The gradient block also contains the energy, which we prefer
# over the energy obtained above because it is more likely to
# be consistent with the gradients. It probably doesn't actually
# make a difference though.
energy = float(fd.readline().split()[1]) * Hartree
grad = []
while line.strip() != '$END':
line = fd.readline()
atom = _atom_re.match(line)
if atom is None:
continue
grad.append(list(map(float, atom.group(3, 4, 5))))
forces = -np.array(grad) * Hartree / Bohr
atoms.calc = SinglePointCalculator(
atoms,
energy=energy,
forces=forces,
charges=charges,
dipole=dipole,
)
return atoms
def clean_userscr(userscr: str, prefix: str) -> None:
for fname in os.listdir(userscr):
tokens = fname.split('.')
if tokens[0] == prefix and tokens[-1] != 'bak':
fold = os.path.join(userscr, fname)
os.rename(fold, fold + '.bak')
def get_userscr(prefix: str, command: str) -> str | None:
prefix_test = prefix + '_test'
command = command.replace('PREFIX', prefix_test)
with workdir(prefix_test, mkdir=True):
try:
call(command, shell=True, timeout=2)
except TimeoutExpired:
pass
try:
with open(f'{prefix_test}.log', encoding='utf-8') as fd:
for line in fd:
if line.startswith('GAMESS supplementary output files'):
return ' '.join(line.split(' ')[8:]).strip()
except FileNotFoundError:
return None
return None
