"""
This module contains classes for launching different MOPAC versions.
"""
# Contributor: Mark A. Watson
import abc
import os
import re
import sys
import time
import traceback
from itertools import zip_longest
import schrodinger.infra.mopac as mopaclib
from schrodinger.application.mopac.results71 import MOPAC71
from schrodinger.application.mopac.results71 import MopacResults71
from schrodinger.application.mopac.results_main import MOPAC_MAIN
from schrodinger.application.mopac.results_main import MOPAC_HEAT_OF_FORMATION_PROP
from schrodinger.application.mopac.results_main import MopacMainTextParser
from schrodinger.application.mopac.results_main import MopacResultsMain
from schrodinger.application.mopac.results_main import split_outfile
from schrodinger.application.mopac.structure_launchers import StructureLauncher
from schrodinger.job import jobcontrol
from schrodinger.utils import subprocess
[docs]class MopacLicenseError(Exception):
pass
[docs]class MopacLauncher(object, metaclass=abc.ABCMeta):
"""
This abstract base class (ABC) is designed to guide developers writing
code to support future MOPAC releases.
It shouldn't be instantiated.
The intention is that a new release will require a new MopacLauncher
subclass and this ABC documents the required interface to be automatically
compliant with the legacy code.
There are currently two subclasses:
(1) MopacLauncher71 - launches a customized open source Mopac 7.1
(2) MopacLauncherMain - launches the currently supported version of Mopac
"""
#
# Required attributes
#
@abc.abstractproperty
def default_method(self):
raise RuntimeError
@abc.abstractproperty
def valid_methods(self):
raise RuntimeError
@abc.abstractproperty
def method_synonyms(self):
raise RuntimeError
@abc.abstractproperty
def extra_keywords(self):
raise RuntimeError
@abc.abstractproperty
def results(self):
raise RuntimeError
#
# Required methods
#
[docs] @abc.abstractmethod
def run(self, inputfile, structures=()):
"""
Run a MOPAC calculation in the local directory.
The optional Structure object is updated with the output data.
:type inputfile: str
:param inputfile: name of MOPAC .mop input file
:type structures: list
:param structures: list of Structure objects
:return: list(job status), outfile
"""
raise RuntimeError
[docs] @abc.abstractmethod
def write_mop_file(self,
cts,
mopfile,
method=None,
geopt=True,
keywords='',
plotMO=None,
gridres=None,
gridext=None):
"""
Write a new .mop MOPAC input file based on a list of Structure objects
and input keywords and settings to be applied to all structures.
:type cts: list of schrodinger.structure.Structure
:param cts: structures to use in writing the file.
:type mopfile: str
:param mopfile: name of .mop file to write.
:type method: str
:param method: The semi-empirical method to use for the calculation.
:type geopt: bool
:param geopt: If True, find the minimum energy geometry.
:type keywords: str
:param keywords: Space-separated keywords to use in MOPAC input file.
:type plotMO: int
:param plotMO: Plot <n> MOs around the HOMO/LUMO gap.
:type gridres: float
:param gridres: Grid resolution for plots.
:type gridext: float
:param gridext: Grid size beyond the nuclei.
"""
raise NotImplementedError
[docs]class MopacLauncher71(MopacLauncher):
"""
This is the API for executing the MOPAC7.1 backend compiled from source,
where we link to a shared library.
"""
# Official MOPAC7.1 method keywords.
_default_method = 'MNDO'
_valid_methods = ['AM1', 'MNDO', 'MNDOD', 'PM3', 'RM1']
# This dictionary maps alternative user-input method names
# to official MOPAC7.1 keywords.
# (These synonyms should all be upper case.)
_method_synonyms = {'MNDO-D': 'MNDOD', 'MNDO/D': 'MNDOD'}
[docs] def __init__(self, energy_only=False):
"""
:type energy_only: bool
:param energy_only: parse only energies from the MOPAC jobs.
"""
self._energy_only = energy_only
self._last_numcal = 0
self._results = None
self._extra_keywords = []
@property
def default_method(self):
return self._default_method
@property
def valid_methods(self):
return self._valid_methods
@property
def method_synonyms(self):
return self._method_synonyms
@property
def extra_keywords(self):
return self._extra_keywords
@property
def results(self):
return self._results
[docs] def write_mop_file(self,
cts,
mopfile,
method=_default_method,
geopt=True,
keywords='',
plotMO=None,
gridres=None,
gridext=None):
"""
Write a new .mop MOPAC input file based on a list of Structure objects
and input keywords and settings to be applied to all structures.
:type cts: list of schrodinger.structure.Structure
:param cts: structures to use in writing the file.
:type mopfile: str
:param mopfile: name of .mop file to write.
:type method: str
:param method: The semi-empirical method to use for the calculation.
:type geopt: bool
:param geopt: If True, find the minimum energy geometry.
:type keywords: str
:param keywords: Space-separated keywords to use in MOPAC input file.
:type plotMO: int
:param plotMO: Plot <n> MOs around the HOMO/LUMO gap.
:type gridres: float
:param gridres: Grid resolution for plots.
:type gridext: float
:param gridext: Grid size beyond the nuclei.
"""
# For MOPAC7.1, plotting is requested through MOPAC keywords
if gridres is not None:
keywords += ' GRIDRES=' + str(gridres)
if gridext is not None:
keywords += ' GRIDEXT=' + str(gridext)
if plotMO is not None:
keywords += ' PLOTMO=' + str(plotMO)
sl = StructureLauncher(self, method, minimize=geopt, keywords=keywords)
with open(mopfile, 'w') as fh:
for ct in cts:
buff = sl.get_mopfile_text(ct)
fh.write(buff.getvalue())
# A blank line between input datasets is required
fh.write('\n')
# Register mopfile with jobcontrol
jobbe = jobcontrol.get_backend()
if jobbe:
jobbe.addOutputFile(mopfile)
[docs] def run(self, mopfile, structures=()):
"""
Run a MOPAC calculation in the local directory with a .mop input file.
Simple support for multiple structures in the input file is provided.
It is important that the number of jobs in the input file matches the
number of structures in the "structures" argument, as structure
properties will be updated assuming a 1:1 mapping. Currently only
total energies are updated when running with multiple structures.
:type mopfile: str
:param mopfile: name of .mop input file (without the suffix)
:type structures: list
:param structures: list of Structure objects
:return: list(job status)
"""
print(f"Starting semi-empirical calculation for input file {mopfile}.")
sys.stdout.flush()
# Register MOPAC7.1 output file with jobcontrol
outfile = mopfile + '.out' # Assume this is true
jobbe = jobcontrol.get_backend()
if jobbe:
jobbe.addOutputFile(outfile)
try:
# Update the stored value of numcal.
self._last_numcal = int(mopaclib.molkst_c.numcal)
# Run a MOPAC7.1 calculation
t0 = time.time()
rc = mopaclib.mopac_main(mopfile)
print(f"{time.time() - t0:.2f} secs to run MOPAC7.1 binary.")
finally:
# Close all files up to 31 except 5, and 6 (stdin, stdout).
# (File 0 is connected to jobname.log, not stderr.)
py_close = mopaclib.pysupport.py_close
for i in range(0, 5):
py_close(i)
for i in range(7, 32):
py_close(i)
# Raise error message if problem was found with MOPAC backend
if rc == 1:
print(
f"\nMOPAC7.1 job {mopfile}.mop failed with user-input (error code 1)."
)
elif rc != 0:
print(f"\nMOPAC7.1 job {mopfile}.mop failed with error code {rc}.")
# Read and store the results
self._results = MopacResults71(mopfile)
self._results.set_return_code(rc)
self._results.set_count(self._last_numcal)
self._results.set_output_file(outfile)
if rc != 0:
print(
f"MOPAC7.1 job {mopfile}.mop failed with the following error message:"
)
print('---')
print(self._results.get_error_text())
print('---')
print(f"Please check {outfile} for additional information.")
if not structures:
status = self._populate_results_from_memory()
elif len(structures) == 1:
status = self._populate_results_from_memory(structures[0])
else:
status = self._populate_results_from_outfile(structures)
return status
def _populate_results_from_outfile(self, structures):
"""
Parse MOPAC7.1 output file to gather results for multiple
structures.
The f2py interface cannot handle more than one result, so
instead of reading data directly from memory we instead parse
the .out file to populate the structures.
:type structures: list(Structure)
:param structures: structures to populate with output data
:return: success status of jobs
"""
assert len(structures) > 1
if not self._energy_only:
msg = f'{self.__class__} must be initialised with'
msg += ' energy_only=True if using multiple structures'
msg += ' because it can currently only parse energies.'
raise RuntimeError(msg)
t1 = time.time()
print(f'\nParsing {self._results._output_file} output file:')
try:
status = self._results.populate_from_outfile(structures)
except Exception as e:
# Bare Exception because we don't know how the parser will fail
print('\nERROR: failed to parse the MOPAC7.1 output'
f' in {self._results._output_file}.')
print('Some results may be missing in the final output.')
print('The parsing error occurred here:')
traceback.print_exc(file=sys.stdout)
# Undefined status, so assume all results are potentially bad
status = [False] * len(structures)
print(f"{time.time() - t1:.2f} secs to parse energies.")
return status
def _populate_results_from_memory(self, st=None):
"""
Read MOPAC7.1 results from memory and optionally populate
a Structure instance with associated properties.
This only works for single-structure MOPAC7.1 jobs due to
problems with memory allocation/re-use etc.
:type st: Structure instance
:param st: structure to optionally populate with output data
:return: success status of job
"""
try:
self._results.populate_from_memory()
except Exception as e:
# Bare Exception because we don't know how the interface will fail
print('\nERROR: failed to interface with the MOPAC7.1 output.')
print('Some results may be missing in the final output.')
print('The interface error occurred here:')
traceback.print_exc(file=sys.stdout)
for m in self.valid_methods:
if self._results.check_key(m):
self._results.set_method(m)
break
else:
self._results.set_method(self.default_method)
# Populate and store a Structure object with the output data.
if st is not None and self._results.statusOk:
self._results.set_final_structure(st)
return [self._results.statusOk]
[docs]class MopacLauncherMain(MopacLauncher):
"""
This is the API for executing the backend of the currently supported
version of MOPAC as an external binary.
"""
MOPAC_EXEC = 'MOPAC2016.exe'
# Official MOPAC_MAIN method keywords
_default_method = 'RM1'
_valid_methods = [
'AM1', 'MNDO', 'MNDOD', 'PM3', 'PM6', 'PM6-D3', 'PM6-DH+', 'PM6-DH2',
'PM6-DH2X', 'PM6-D3H4', 'PM6-D3H4X', 'PM7', 'PM7-TS', 'RM1'
]
# This dictionary maps alternative user-input method names
# to official MOPAC keywords.
# (These synonyms should all be upper case.)
_method_synonyms = {'MNDO-D': 'MNDOD', 'MNDO/D': 'MNDOD'}
[docs] def __init__(self, energy_only=False):
"""
:type energy_only: bool
:param energy_only: parse only energies from the MOPAC jobs.
"""
self._energy_only = energy_only
self._results = None
self._extra_keywords = ['THREADS=1']
if not self._energy_only:
self._extra_keywords += ['AUX(PRECISION=10)']
self._plotMO = None
self._gridres = 0.0
self._gridext = 0.0
@property
def default_method(self):
return self._default_method
@property
def valid_methods(self):
return self._valid_methods
@property
def method_synonyms(self):
return self._method_synonyms
@property
def extra_keywords(self):
return self._extra_keywords
@property
def results(self):
return self._results
[docs] def write_mop_file(self,
cts,
mopfile,
method=_default_method,
geopt=True,
keywords='',
plotMO=None,
gridres=None,
gridext=None):
"""
Write a new .mop MOPAC input file based on a list of Structure objects
and input keywords and settings to be applied to all structures.
:type cts: list of schrodinger.structure.Structure
:param cts: Structures to use in writing the file.
:type mopfile: str
:param mopfile: name of .mop file to write (eg - infile.mop)
:type method: str
:param method: The semi-empirical method to use for the calculation.
:type geopt: bool
:param geopt: If True, find the minimum energy geometry.
:type keywords: str
:param keywords: Space-separated keywords to use in MOPAC input file.
:type plotMO: int
:param plotMO: Plot <n> MOs around the HOMO/LUMO gap.
:type gridres: float
:param gridres: Grid resolution for plots.
:type gridext: float
:param gridext: Grid size beyond the nuclei.
"""
# For MOPAC_MAIN, plotting is done post-calculation via this class
if gridres is not None:
self._gridres = gridres
if gridext is not None:
self._gridext = gridext
if plotMO is not None:
self._plotMO = plotMO
# Need to print all the eigenvectors etc in
# the .aux file even for large molecules.
self._extra_keywords.append('LARGE')
sl = StructureLauncher(self, method, minimize=geopt, keywords=keywords)
with open(mopfile, 'w') as fh:
for ct in cts:
buff = sl.get_mopfile_text(ct)
fh.write(buff.getvalue())
# A blank line between input datasets is required
fh.write('\n')
# Register mopfile with jobcontrol
jobbe = jobcontrol.get_backend()
if jobbe:
jobbe.addOutputFile(mopfile)
def _execute_mopac_binary(self, inputfile):
"""
Execute currently supported MOPAC binary as a subprocess.
:type inputfile: str
:param inputfile: name of .mop input file (eg - infile.mop)
"""
t0 = time.time()
rc = subprocess.call([self.MOPAC_EXEC, inputfile])
print(f"{time.time() - t0:.2f} secs to run {MOPAC_MAIN} binary.")
if rc != 0:
print(f"{MOPAC_MAIN} exited with nonzero status ({rc})")
return rc
[docs] def run(self, mopfile, structs=()):
"""
Run a MOPAC calculation in the local directory with a .mop input file.
If optional list of Structure instances is given, then parse the output
files for data with which to populate the structure properties.
Support for multiple structures in the input file is provided.
It is important that the number of jobs in the input file matches the
number of structures in the "structures" argument, as structure
properties will be updated assuming a 1:1 mapping.
:type mopfile: str
:param mopfile: name of .mop input file (eg - infile.mop)
:type structs: list
:param structs: list of Structure objects
:return: list(job status)
"""
print(f"Starting {MOPAC_MAIN} calculation for input file {mopfile}.")
sys.stdout.flush()
# Register MOPAC_MAIN binary output file with jobcontrol
mopfile_base = mopfile.split('.mop')[0]
outfile = mopfile_base + '.out' # Assume this is true
auxfile = mopfile_base + '.aux' # Assume this is true
jobbe = jobcontrol.get_backend()
if jobbe:
jobbe.addOutputFile(outfile)
jobbe.addOutputFile(auxfile)
# Run calculation with MOPAC_MAIN external binary
rc = self._execute_mopac_binary(mopfile)
# Parse MOPAC_MAIN output files
status = [True]
t0 = time.time()
if structs:
if self._energy_only:
status = self._populate_energies_only(outfile, structs)
else:
status = self._populate_all_results(outfile, auxfile, structs)
if self._plotMO:
print(
f"{time.time() - t0:.2f} secs to parse results and plot MOs."
)
else:
print(f"{time.time() - t0:.2f} secs to parse results.")
return status
def _populate_energies_only(self, outfile, structures):
"""
Parse .out file only for energies and no other properties
and populate structure properties.
:type outfile: str
:param outfile: name of MOPAC output file
:type structures: list
:param structures: list of Structure objects
:return: list(job status)
"""
status = []
if os.path.exists(outfile):
with open(outfile, 'r') as fh:
parser = MopacMainTextParser(fh)
all_results = parser.parse_multiple_jobs(minimal=True)
# Store results as .mae properties.
# If lists have different lengths, we assume results has been
# truncated due to e.g. job crash, so we only update CT's as
# far as we can.
i = 0
for st, properties in zip_longest(structures,
all_results,
fillvalue={}):
i += 1
stat = False
if properties:
for k, v in properties.items():
if k[0] == 'mae':
st.property[k[1]] = v
# Only mark as successful if the energy was found
if MOPAC_HEAT_OF_FORMATION_PROP in st.property:
stat = True
status.append(stat)
if st:
mopfile = re.sub('.out$', '', outfile)
st.property['s_j_jname'] = f'{mopfile}_{i}'
else:
status = [False] * len(structures)
print('\nNo successful MOPAC output file {outfile} was found'
' so no results will be reported.')
return status
def _populate_all_results(self, full_outfile, full_auxfile, structures):
"""
# Parse .out and .aux files and populate structure properties
:type full_outfile: str
:param full_outfile: name of MOPAC output file
:type full_auxfile: str
:param full_auxfile: name of MOPAC auxiliary output file
:type structures: list
:param structures: list of Structure objects
:return: list(job status)
"""
outfiles = split_outfile(full_outfile)
auxfiles = split_outfile(full_auxfile)
status = []
for idx, (ct, outfile, auxfile) in enumerate(
zip_longest(structures, outfiles, auxfiles), 1):
# Grep data from the MOPAC output files and store it
if outfile is not None:
if auxfile is None:
print('\nNo MOPAC auxiliary output was found for'
f' structure {idx} ({ct.title}).\nThis usually'
' indicates a problem with the calculation and'
' it is likely some properties will be missing'
' in the output .mae file.\nRecommend inspecting'
f' the MOPAC {full_outfile} file for additional'
' error messages.')
self._results = MopacResultsMain(outfile, auxfile, ct)
status.append(self._results.statusOk)
if ct is not None:
ct.property['s_j_jname'] = self._results.name
else:
status.append(False)
print('\nNo successful MOPAC output was found for'
f' structure {idx} ({ct.title}).\nNo results will be'
' reported for this structure.')
if self._plotMO:
# Build .vis files from MOPAC output files.
try:
self._results.write_vis_files(self._plotMO, self._gridres,
self._gridext)
except Exception:
msg = "\nERROR: failed to plot requested MOs.\n"
msg += "Note: plotting of UHF wavefunctions is currently"
msg += f" unsupported with -{MOPAC_MAIN}"
print(msg)
return status
#=============================================================================
# When new versions of MOPAC are supported, add calls to return
# their MopacLauncher subclasses here
[docs]def get_mopac_launcher(version, energy_only=False):
"""
Return a MopacLauncher object for the requested MOPAC version.
The optional argument <energy_only> may be used to tell the launcher
to only parse energies from the MOPAC job when populating any structure
objects that might be provided.
:type version: str
:param version: MOPAC version number
:type energy_only: bool
:param energy_only: parse only energies from the MOPAC jobs.
"""
if version == MOPAC_MAIN:
mopac_launcher = MopacLauncherMain(energy_only)
elif version == MOPAC71:
mopac_launcher = MopacLauncher71(energy_only)
else:
msg = f'MOPAC version {version} is not supported'
raise ValueError(msg)
if not isinstance(mopac_launcher, MopacLauncher):
raise TypeError('mopac_launcher must be a MopacLauncher object')
return mopac_launcher