"""
Classes for parsing Jaguar output files and accessing output properties
programmatically.
Copyright Schrodinger, LLC. All rights reserved.
"""
import contextlib
import copy
import gzip
import os.path
import re
import textwrap
from math import ceil
from past.utils import old_div
import numpy as np
import schrodinger.infra.mm as mm
import schrodinger.structure as structure
import schrodinger.utils.fileutils as futils
from schrodinger.application.jaguar import constants
from schrodinger.application.jaguar import jaguar_diff
from schrodinger.application.jaguar.results import JaguarOptions
from schrodinger.application.jaguar.results import JaguarResults
from schrodinger.application.jaguar.results import _Attribute
from schrodinger.application.jaguar.textparser import JaguarParseError
from schrodinger.application.jaguar.textparser import TextParser
from schrodinger.utils import units
restart_re = re.compile(r"^(.+)\.(\d\d)$")
#-----------------------------------------------------------------------------
[docs]def join_structure_files(name, files):
"""
Pull the CTs out of each .mae file in the 'files' list and put
them into a new .mae file called 'name', which is expected to
already contain the .mae extension. If a file called 'name'
already exists, it will be overwritten.
"""
all_cts = structure.MaestroWriter(name)
for f in files:
if futils.get_structure_file_format(f) == futils.MAESTRO:
if os.path.isfile(f):
for st in structure.StructureReader(f):
all_cts.append(st)
else:
print("Could not find structure file: %s" % f)
all_cts.close()
return
[docs]def restart_name(name):
"""
Generate a restart name from a job name. If there is no number suffix
for the job this will be jobname.01; otherwise the suffix will be
incremented.
"""
m = restart_re.search(name)
if m:
basename = m.group(1)
count = int(m.group(2))
count += 1
return "%s.%02d" % (basename, count)
else:
return "%s.01" % (name)
[docs]def get_section(fname, sect_start, sect_end=None):
"""
Generator of lines for a section of a (jaguar) file.
The section is defined by a start and optionally by an end.
Note that if sect_end is not specified the remainder of the file
after a sect_start will be returned by the generator.
Note that, if multiple sections defined by sect_start and sect_end
are present in the file, all will appear in the order they appear
:type fname: string
:param fname: Name of file to parse for section
:type sect_start: string
:param sect_start: Indicator of where section begins, i.e. the
first line of desired section
:type sect_end: string (or None)
:param sect_end: If specified, indicator of the end of a section, i.e.
the last line in the desired section
"""
in_sec = False
with open(fname, "r") as fin:
for line in fin:
if line.strip() == sect_start:
in_sec = True
elif sect_end is not None and line.strip() == sect_end:
in_sec = False
else:
if in_sec:
yield line
[docs]def get_overlap_dim(fname):
"""
Get size of overlap matrix from output file.
Unfortunately, this does not necessarily match the nbasis of
JaguarOutput (possibly the overlap matrix is always in caretesians,
i.e. numd=6, but since we aren't 100% sure, we will just check how many
rows are in the printed overlap matrix.
:type fname: string
:param fname: name of output file
:rtype: int
:return: dimension of overlap matrix
"""
olap_sec = get_section(fname, 'X*S*X')
# first line is a blank space
next(olap_sec)
data = next(olap_sec).split()
# Iterate through rows until we encounter a blank line
# The row number on the previous line is the dimension
while data:
nbas = int(data[0])
data = next(olap_sec).split()
return nbas
[docs]def parse_overlap_mat(fname):
"""
Parse overlap matrix from output file
:type fname: string
:param fname: name of output file
:return: nbas x nbas overlap matrix as numpy array
"""
if not fname.endswith(".out"):
raise ValueError("%s is not an output file" % fname)
#Get number of AOs
nbas = get_overlap_dim(fname)
olap = np.zeros((nbas, nbas))
# top triangle printed 5 cols at a time
ncycles = ceil(nbas / 5)
data = []
olap_sec = get_section(fname, 'X*S*X')
for icycle in range(ncycles):
# first line in each data block is a blank space
next(olap_sec)
# second line in each data block has the column indices
data = next(olap_sec).split()
j = int(data[0]) - 1
# each following line has a row index and entries in the matrix
for row in range(nbas):
data = next(olap_sec).split()
i = int(data[0])
for idx, value in enumerate(map(float, data[1:])):
olap[i - 1, j + idx] = value
return olap
#-----------------------------------------------------------------------------
# JaguarOutput attributes, defined as simple objects to allow some sanity
# checking and programmatic manipulation of a big list.
_tmp_attributes = [
_Attribute("job_id",
description="Job ID",
datatype="string",
comparison=None),
_Attribute("opts",
description="Job Calculation Options",
datatype="JaguarOptions",
init=JaguarOptions,
comparison=None),
_Attribute("host",
description="Job Host",
datatype="string",
comparison=None),
_Attribute("mae_out",
description="Maestro output file",
datatype="string",
comparison=None),
_Attribute("mae_in",
description="Maestro input file",
datatype="string",
comparison=None),
_Attribute(
"status",
description=
"Job status - set to 0, 1, or 2 corresponding to UNKNOWN, OK, or SPLAT respectively",
pretty_name="Job Status",
datatype="int",
init=lambda: JaguarOutput.UNKNOWN),
_Attribute("fatal_error",
description="Error message in the event the job failed",
pretty_name="Fatal Error Msg",
datatype="string"),
_Attribute("fatal_errorno",
description="Error number in the event the job failed",
pretty_name='Fatal Error Number',
datatype="string"),
_Attribute("glibc",
description="Reported glibc version",
datatype="string",
comparison=None),
_Attribute("lastexe",
description="Last Jaguar Executable Used",
datatype="string",
comparison=None),
_Attribute(
"start_time",
description="Starting time of the calculation - NOT CURRENTLY SUPPORTED",
pretty_name="Start Time",
datatype="datetime.datetime",
comparison=None),
_Attribute(
"end_time",
description="Ending time of the calculation - NOT CURRENTLY SUPPORTED",
pretty_name="End Time",
datatype="datetime.datetime",
comparison=None),
_Attribute("symmetrized",
description="Whether the geometry has been symmetrized or not",
pretty_name='Is Symmetrized',
datatype="bool",
comparison=None),
_Attribute("geopt_stuck",
description="Whether the geopt or tsopt got stuck",
pretty_name='Geopt Is Stuck',
datatype="bool",
comparison=None),
#CURRENTLY SKIPPING BELOW ATTRS ALAN
_Attribute(
"sm_geopt_step",
description=
"results for individual SM geopt steps (see note below) - NOT CURRENTLY SUPPORTED",
datatype="list of lists of JaguarResults",
init=list,
comparison=None),
_Attribute(
"sm_n_steps",
description=
"results for individual string method steps (see note below) - NOT CURRENTLY SUPPORTED",
datatype="list of lists of JaguarResults",
comparison=None),
_Attribute(
"geopt_step",
description="results for each step in a geopt - NOT CURRENTLY SUPPORTED",
datatype="list of JaguarResults",
init=list,
comparison=None),
_Attribute(
"gas_phase_geopt_step",
description=
"results for the gas phase geometry optimization in a solvated geopt - NOT CURRENTLY SUPPORTED",
datatype="list of JaguarResults",
init=list,
comparison=None),
_Attribute(
"solution_phase_geopt_step",
description=
"results for the solution phase geometry optimization in a solvated geopt - NOT CURRENTLY SUPPORTED",
datatype="list of JaguarResults",
init=list,
comparison=None),
_Attribute(
"scan_geopt_step",
description=
"results for individual scan geopt steps (see note below) - NOT CURRENTLY SUPPORTED",
datatype="list of lists of JaguarResults",
init=list,
comparison=None),
_Attribute(
"irc_geopt_step",
description=
"results for individual IRC geopt steps (see note below) - NOT CURRENTLY SUPPORTED",
datatype="list of lists of JaguarResults",
init=list,
comparison=None),
_Attribute(
"input_geometry",
description=
"the input geometry, unmodified by symmetrization - NOT CURRENTLY SUPPORTED",
datatype="schrodinger.structure.Structure",
comparison=None),
_Attribute(
"input_geometry2",
description=
"the input geometry from zmat2, unmodified by symmetrization - NOT CURRENTLY SUPPORTED",
datatype="schrodinger.structure.Structure",
comparison=None),
_Attribute(
"input_geometry3",
description=
"the input geometry from zmat3, unmodified by symmetrization - NOT CURRENTLY SUPPORTED",
datatype="schrodinger.structure.Structure",
comparison=None),
_Attribute("convergence_category",
description=
"Results for the convergence category - NOT CURRENTLY SUPPORTED",
datatype="list of JaguarResults",
init=list,
comparison=None)
]
_jo_attributes = [
_Attribute("method",
description="Calculation Type",
datatype="string",
init=lambda: constants.ctHF),
_Attribute("functional", description="DFT Functional", datatype="string"),
_Attribute("basis", description="Basis Set", datatype="string"),
_Attribute("nbasis",
description="Number of Basis Functions",
pretty_name='Num Basis Fxns',
datatype="int"),
_Attribute("nelectron",
description="Number of Electrons",
pretty_name='Num Electrons',
datatype="int"),
_Attribute("qm_atoms",
description="Number of QM Atoms",
pretty_name='Num QM Atoms',
datatype="int"),
_Attribute("stoichiometry",
description="Stoichiometry of input geometry",
pretty_name='Input Stoichiometry',
datatype="string"),
_Attribute("point_group",
description="Molecular point group of the input molecule",
pretty_name='Mol Point Group',
datatype="string"),
_Attribute("point_group_used",
description="Point group used in the calculation",
pretty_name='Used Point Group',
datatype="string"),
_Attribute("charge",
description="Molecular charge of Input Structure",
pretty_name='Mol Charge',
datatype="int",
init=lambda: 0),
_Attribute("multiplicity",
description="Spin Multiplicity of Input Structure",
pretty_name='Multiplicity',
datatype="int",
init=lambda: 1),
_Attribute("mol_weight",
description="Molecular weight of input geometry",
pretty_name='Mol Weight',
datatype="float",
units=constants.unAtomicMassUnit,
comparison=jaguar_diff._FLOAT_PRECISION,
precision="mol_weight_precision"),
_Attribute("coords_opt",
description="Number of optimization coordinates",
pretty_name='Num Opt Coords',
datatype="int"),
_Attribute("coords_ind",
description="Number of independent coordinates",
pretty_name='Num Independent Coords',
datatype="int"),
_Attribute("coords_nred",
description="Number of non-redundant coordinates",
pretty_name='Num Non-Redundant Coords',
datatype="int"),
_Attribute("coords_frozen",
description="Number of frozen coordinates",
pretty_name='Num Frozen Coords',
datatype="int"),
_Attribute("coords_harmonic",
description="number of harmonic constraints",
pretty_name='Num Harmonic Constraints',
datatype="int"),
_Attribute(
"ts_component_descriptions",
description="Descriptions of the transition state vector components",
pretty_name='TS Component Descriptions',
datatype="string"),
_Attribute("_sm_n_points",
description="number of string method points",
datatype="int"),
#BELOW ATTRS NOT CURRENTLY HANDLED ALAN
_Attribute("scan_coords",
description="scan coordinates - NOT CURRENTLY SUPPORTED",
datatype="list of Scan objects",
init=list),
_Attribute(
"path_structures",
description="path structures for IRC/RSM jobs - NOT CURRENTLY SUPPORTED",
datatype="list of Structure instances",
comparison=jaguar_diff._PATH_PRECISION,
precision="rmsd_precision",
init=list),
]
_jo_attributes.sort()
## Create a list of attributes for the JaguarOutput docstring.
#
_tmp_description_list = ["A class for accessing Jaguar output properties.\n"]
_tmp_description_list.extend(
textwrap.wrap(
"""Stores a JaguarResults object for each geometry (multiple geometry
optimization steps in the geopt_step list, multiple scan steps in the
scan_step list). Each JaguarResults object holds all properties specific
to a given geometry. (See the JaguarResults documentation for more
detail.) Each JaguarResults object also holds a list of
JaguarAtomicResults objects for atomic specific properties."""))
_tmp_description_list.append("\nGeneral Attributes\n")
_tmp_description_list.extend([item.summary() for item in _tmp_attributes])
_tmp_description_list.extend(
textwrap.wrap(
"""The scan_geopt_step and irc_geopt_step attributes are both lists of
lists of JaguarResults. Each element of the list represents the geopt
steps taken within an individual scan or IRC step. In the case where
geometry optimization is not done with the scan or IRC calculation, each
element is a list with a single JaguarResults element."""))
_tmp_description_list.append("")
_tmp_description_list.extend(
textwrap.wrap("""The properties scan_step
and irc_step are probably more useful. They provide a list of the final
geometries for each scan/IRC step."""))
_tmp_description_list.append("\nProperty Attributes\n")
_tmp_description_list.extend([item.summary() for item in _jo_attributes])
_jo_docstring = "\n".join(_tmp_description_list)
# Add the attributes I wanted stuck at the head of the docstring into the
# _jo_attributes list for later manipulation.
_jo_attributes.extend(_tmp_attributes)
del (_tmp_description_list)
del (_tmp_attributes)
#-----------------------------------------------------------------------------
[docs]class JaguarOutput(object):
# Define the docstring above, somewhat automatically from the list of
# attributes.
__doc__ = _jo_docstring
UNKNOWN, OK, SPLAT = list(range(3))
_attributes = _jo_attributes
# Attributes that will be created in the __init__ method.
#_attributes = _jo_attributes
# This precision specification is a reasonable default but can be
# overridden if you wish to compare to a different value.
mol_weight_precision = 1.0e-2
rmsd_precision = 1.0e-4
# _text_parser_class is set after the TextParser class definition to
# keep from shifting the order of the whole file around.
_text_parser_class = TextParser
[docs] def __init__(self, output=None, partial_ok=False):
"""
Initialize from an output filename or output name.
:param str output: The name of the Jaguar output file to read. Can be
the .out or .out.gz name or a file with the same basename as the
desired .out/.out.gz file.
:param bool partial_ok: Do not raise an exception if parsing fails
partway through the file. Instead, process what has been read and
store the parsing error in the parsing_error attribute. The
resulting JaguarOutput object may be in an incomplete state.
Exceptions:
:raise IOError: Raised if output file cannot be found.
:raise JaguarParseError: Raised if the output file can't be parsed.
If this is raised, the state of the resulting object is not
guaranteed to be useful.
:raise StopIteration: Raised in some cases if the file ends unexpectedly
If this is raised, the state of the resulting object is not
guaranteed to be useful.
"""
#Make each JaguarOutput instance have its own self._attributes
self._attributes = copy.deepcopy(JaguarOutput._attributes)
# See the _tmp_attributes and _jo_attributes lists above (or the
# generated docstring) for descriptions of all the attributes that
# are present in a JaguarOutput instance.
for attr in self._attributes:
self.__dict__[attr.name] = attr.init()
self._scan_step = None
self._irc_step = None
self._path_step = []
self._sm_geopt_step = []
self._sm_n_points = 0
self.parsing_error = None
# self._results is the current set of results to set properties etc.
# in parsing.
self._results = JaguarResults()
# self.last_results is the last complete set of results
self.last_results = self._results
self._mmjag_initialized = False
try:
mm.mmjag_initialize(mm.error_handler)
self._mmjag_initialized = True
except:
pass
# Allow construction of an empty object that doesn't get loaded from
# a file.
if not output:
return
if output.endswith(".out"):
self.name = os.path.basename(output)[:-4]
self.filename = futils.extended_windows_path(output)
elif output.endswith(".out.gz"):
self.name = os.path.basename(output)[:-7]
self.filename = futils.extended_windows_path(output)
else:
self.name = os.path.basename(output)
for suffix in [".out", ".out.gz"]:
filename = futils.extended_windows_path(output + suffix)
if os.path.exists(filename):
self.filename = filename
break
else:
raise IOError("No output file could be found for '%s'" % output)
# Load up the object.
if self.filename.endswith(".gz"):
file_iter = gzip.open(self.filename, 'rt')
else:
file_iter = open(self.filename, 'r')
with contextlib.closing(file_iter):
tp = self._text_parser_class(self, file_iter)
try:
tp.parse()
except (JaguarParseError, StopIteration) as err:
if isinstance(err, StopIteration):
msg = 'Unexpected end of .out file'
else:
msg = str(err)
if partial_ok:
# Caller requested to process the partial output rather than
# fail on an error
self.parsing_error = msg
# Try to finish processing what output was parsed
try:
tp._finalize(tp.jag_out)
except JaguarParseError as err2:
self.parsing_error += '\n' + str(err2)
tp.jag_out = None
else:
# Partial results are not OK, raise the error that occurred
raise
# Construct derived attributes from parsed values
self._results.derived_attrs.buildDerivedAttrs(
self._results,
self) #this line seems not very kosher, or at least inelegant
def __del__(self):
if self._mmjag_initialized:
mm.mmjag_terminate()
self._mmjag_initialized = False
def __getattr__(self, attr):
"""
Fall back to the last_results object for attribute access.
If an attribute isn't present, try getting it from the last_results
attribute.
"""
return getattr(self.last_results, attr)
# Properties
def _getRestart(self):
# See restart property doc.
return restart_name(self.name)
restart = property(_getRestart,
doc="""
Return the restart name for this output object.
""")
[docs] def getDuration(self):
# See duration property doc.
if self.start_time and self.end_time:
return self.end_time - self.start_time
else:
return None
duration = property(getDuration,
doc="""
Return the duration of the job as a datetime.timedelta object.
""")
[docs] def getScanStep(self):
# See scan_step property doc.
if self._scan_step is None:
self._scan_step = [step[-1] for step in self.scan_geopt_step]
return self._scan_step
scan_step = property(getScanStep,
doc="""
Return a list of final scan geometries for each scan step.
""")
[docs] def getIrcStep(self):
# See irc_step property doc.
if self._irc_step is None:
self._irc_step = [step[-1] for step in self.irc_geopt_step]
return self._irc_step
irc_step = property(getIrcStep,
doc="""
Return a list of final IRC geometries for each IRC step.
""")
def __eq__(self, other):
if self.diff(other, short_circuit=True):
return False
return True
def __ne__(self, other):
return not self.__eq__(other)
[docs] def diff(self, other, factor=1.0, short_circuit=False):
"""
Return a list of all differing attributes.
Each item is a tuple of (property name, self value, other value).
Note that the property names are not necessarily usable in getattr;
some may be properties of atoms, such as "atom[1].forces".
Parameters
other (JaguarOutput)
The JaguarOutput instance to compare against.
factor (float)
A constant factor to multiply all float comparison tolerances
by.
short_circuit (boolean)
If true, will return immediately upon finding a difference. The
values in the tuple will both be None in this case.
"""
different = jaguar_diff._diff(self,
other,
short_circuit=short_circuit,
factor=factor)
if short_circuit and different:
return [(item, None, None) for item in different]
different.update(
self.last_results.diff(other.last_results, factor=factor))
if short_circuit and different:
return [(item, None, None) for item in different]
difflist = []
for prop in different:
# Handle list (of list, of list ...) attributes
if isinstance(prop, tuple):
if isinstance(prop[0], tuple):
if isinstance(prop[0][0], tuple):
# e.g. prop = ((('atom', 2), 'forces'), 2)
((attr1, idx1), attr2), idx2 = prop
name = '%s[%d].%s[%d]' % (attr1, idx1, attr2, idx2)
t = (name, getattr(getattr(self, attr1)[idx1],
attr2)[idx2],
getattr(getattr(other, attr1)[idx1], attr2)[idx2])
difflist.append(t)
else:
# e.g. prop = (('atom', 2), 'charge_esp')
(attr1, idx1), attr2 = prop
name = '%s[%d].%s' % (attr1, idx1, attr2)
t = (name, getattr(getattr(self, attr1)[idx1], attr2),
getattr(getattr(other, attr1)[idx1], attr2))
difflist.append(t)
else:
# e.g. prop = ('excitation_energies', 2)
attr, idx = prop
name = '%s[%d]' % prop
t = (name, getattr(self, attr)[idx], getattr(other,
attr)[idx])
difflist.append(t)
else:
t = (prop, getattr(self, prop), getattr(other, prop, None))
difflist.append(t)
difflist.sort()
return difflist
def _defaultStructureProperties(self):
"""
Return a list of name, value tuples that are the default properties
that should be added to all written mae structures for a given
output file.
"""
props = []
mae_in = None
if self.mae_in:
if os.path.exists(self.mae_in):
mae_in = self.mae_in
else:
dir_, file_ = os.path.split(self.filename)
mae_in = os.path.join(dir_, self.mae_in)
if not os.path.exists(mae_in):
mae_in = None
if mae_in:
st = structure.Structure.read(mae_in)
for key in ['s_m_title', 's_m_entry_name', 's_m_entry_id']:
value = st.property.get(key)
if value is not None:
props.append((key, value))
props.append(('i_m_Molecular_charge', self.charge))
props.append(('i_m_Spin_multiplicity', self.multiplicity))
if self.opts.solvation:
props.append(('s_j_QM_Method', self.method + "/SOLV"))
else:
props.append(('s_j_QM_Method', self.method))
props.append(('s_j_QM_Basis', self.basis))
return props
@property
def path_structures(self):
"""
List of structures along path for IRC or RSM jobs, empty list otherwise
"""
if self._sm_geopt_step:
npts = self._sm_n_points
return self.getStructures()[-npts:]
elif self.irc_geopt_step:
return self.getStructures()
else:
return list()
[docs] def getStructures(self):
"""
Get Structure objects for the geometries in the output file.
If this job is a geometry optimization, it will contain geometries
for all steps. If it's a scan, it will contain the geometries for
each scan point (but only the end geometries if it's a relaxed
scan).
Return a list of Structure objects.
"""
structures = []
props = self._defaultStructureProperties()
# will retrieve all geometries
if self._sm_geopt_step:
for result in self._sm_geopt_step:
st = result.getStructure(props)
structures.append(st)
elif self._path_step:
for result in self._path_step:
structures.append(result.getStructure(props))
elif self.scan_geopt_step:
# Just create a slot for the property we will update in each
# step of the iterator below.
props.append(None)
ix = 0
for step in self.scan_geopt_step:
ix += 1
props[-1] = ('i_j_ScanPt_Num', ix)
st = step[-1].getStructure(props)
structures.append(st)
elif self.geopt_step:
# Just create a slot for the property we will update in each
# step of the iterator below.
props.append(None)
ix = 0
for step in self.geopt_step:
ix += 1
props[-1] = ('i_j_Geom_Iter_Num', ix)
st = step.getStructure(props)
structures.append(st)
elif self.irc_geopt_step:
for step in self.irc_geopt_step:
st = step[-1].getStructure(props)
structures.append(st)
else:
st = self.last_results.getStructure(props)
structures = [st]
return structures
[docs] def getStructure(self):
"""
Return a structure object for the last geometry in the file.
"""
structures = self.getStructures()
return structures[-1]
[docs] def write(self,
filename=None,
mimic_backend=False,
add_title=False,
add_entry=False):
"""
Write a maestro file for the structure in the output file.
Note that this method overwrites any file with the same pathname.
If this job is a geometry optimization, it will contain geometries
for all steps. If it's a scan, it will contain the geometries for
each scan point (but only the end geometries if it's a relaxed
scan).
filename (str)
The filename to write to; if not specified, defaults to the
restart name with the '.mae' suffix.
mimic_backend (bool)
If false, all geometry optimization structures will be written.
If true, the geometry optimization structures will be written as
in regular jobs; by default, only the last geometry will be
used, but if ip472 is greater than 1, all geometries will be
included.
add_title (bool)
If true, then an empty title will be replaced with the output
file's jobname.
add_entry (bool)
If the entry name is empty or starts with 'Scratch' it will be
replaced with the output file's jobname.
"""
mae_format = "maestro"
if not filename:
filename = os.path.join(os.path.dirname(self.filename),
self.restart + ".mae")
if os.path.exists(filename):
os.remove(filename)
structures = self.getStructures()
if mimic_backend and self.opts.ip[472] <= 1:
st = structures[-1]
if add_title:
if not st.title or st.title.startswith("Scratch"):
st.title = self.name
if add_entry:
entry_name = st.property['s_m_entry_name']
if not entry_name or entry_name.startswith("Scratch"):
st.property['s_m_entry_name'] = self.name
st.write(filename, format=mae_format)
else:
for st in structures:
st.append(filename, format=mae_format)
return
[docs] def writeGrd(self, filename):
"""
Write a .grd file for 1D or 2D visualization of scans in maestro to
file 'filename'.
If the job is not a scan job, this will raise a RuntimeError.
"""
#FIXME: this function is too complex.
if len(self.scan_coords) < 1 or len(self.scan_coords) > 2:
raise RuntimeError("Can only write a .grd file for a scan job "
"with one or two scan coordinates.")
grd = open(filename, "w")
initial_values = [0.0, 0.0, 0.0]
steps = [0, 0, 0]
increments = [0.0, 0.0, 0.0]
for index, scan_coord in enumerate(self.scan_coords):
if scan_coord.at_values:
raise RuntimeError(
"Cannot write a .grd file for a job with an 'at values' scan."
)
increments[index] = scan_coord.step_size
if scan_coord.step_size > 0:
initial_values[index] = scan_coord.initial
else:
initial_values[index] = scan_coord.final
steps[index] = scan_coord.steps
# Make two header lines
grd.write("BMIN 0 0 0 0\n")
grd.write("BMIN 0 0 0 0\n")
struct = self.scan_geopt_step[0][0].getStructure()
# Add "origins" line
fmt0 = "%5d%12.6f%12.6f%12.6f\n"
grd.write(fmt0 % (len(struct.atom), initial_values[0],
initial_values[1], initial_values[2]))
# Add "increments" lines
grd.write(fmt0 % (steps[0], abs(increments[0]), 0, 0))
grd.write(fmt0 % (steps[1], 0, abs(increments[1]), 0))
grd.write(fmt0 % (steps[2], 0, 0, abs(increments[2])))
# Add molecule geometry
fmt1 = "%5d%12.6f%12.6f%12.6f%12.6f\n"
for at in struct.atom:
grd.write(fmt1 %
(at.atom_type, at.partial_charge, at.x, at.y, at.z))
# Print scan energy values
# Deal with scan points where the calculation was skipped due to
# close contact between atoms
energies = []
valid_energies = []
for scanpt in self.scan_step:
energy = scanpt.energy
energies.append(energy)
if energy is not None:
valid_energies.append(energy)
max_energy = max(valid_energies)
min_energy = min(valid_energies)
missing_energy = max_energy + old_div(100.0, constants.kcal_per_hartree)
if len(self.scan_coords) == 1:
if len(valid_energies) != len(self.scan_step):
for ix, energy in enumerate(energies):
if energy is None:
energies[ix] = missing_energy
# Deal with reversed coordinates for the 1-d case
if increments[0] < 0:
energies.reverse()
else:
# Deal with the zig-zag scanning algorithm used in the back end.
x = self.scan_coords[0].var_name
y = self.scan_coords[1].var_name
x_init = initial_values[0]
y_init = initial_values[1]
y_steps = self.scan_coords[1].steps
x_step_size = abs(self.scan_coords[0].step_size)
y_step_size = abs(self.scan_coords[1].step_size)
energies = [missing_energy] * len(self.scan_step)
for scanpt in self.scan_step:
if scanpt.scan_value:
x_index = int(
old_div(round(scanpt.scan_value[x] - x_init),
x_step_size))
y_index = int(
old_div(round(scanpt.scan_value[y] - y_init),
y_step_size))
energies[y_index + x_index * y_steps] = scanpt.energy
for e in energies:
grd.write("%12.6f\n" %
((e - min_energy) * units.KJOULES_PER_MOL_PER_HARTREE,))
grd.close()