"""
This module contains classes for launching different MOPAC versions.
"""
# Contributor: Mark A. Watson
import abc
import os
import sys
import time
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 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.infra import mm
from schrodinger.infra.mmcheck import MmException
from schrodinger.utils import fileutils
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')
[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), outfile
"""
print(f"Starting semi-empirical calculation for input file {mopfile}.")
sys.stdout.flush()
outfile = mopfile + '.out' # Assume this is true
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)
# Store the results
self._results = MopacResults71()
self._results.set_return_code(rc)
self._results.set_count(self._last_numcal)
self._results.set_output_file(outfile)
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)
status = [self._results.statusOk]
if len(structures) == 1 and self._results.statusOk:
# Populate and store a Structure object with the output data.
self._results.set_final_structure(structures[0], mopfile)
elif len(structures) > 1:
if not self._energy_only:
msg = f'{self.__class__} must be initialised with'
msg += ' energy_only=True if using multiple structures.'
raise RuntimeError(msg)
# The f2py interface cannot handle more than one result, so
# we simply parse the .out file to populate the structure data.
t1 = time.time()
status = self._results.populate_from_outfile(structures)
print(f"{time.time() - t1:.2f} secs to parse energies.")
# Raise error message if problem was found
if rc == 1:
print("MOPAC7.1 failure with input (error code 1).")
if self._results is not None:
print(self._results.get_error_text())
elif rc != 0:
print(f"MOPAC7.1 failure with error code {rc}.")
if self._results is not None:
print(self._results.get_error_text())
return status, outfile
[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
self._tmpdir = None
def __del__(self):
if isinstance(self._tmpdir, str) and os.path.exists(self._tmpdir):
fileutils.force_rmtree(self._tmpdir, ignore_errors=True)
@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 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')
def _execute_mopac_binary(self, inputfile):
"""
Execute currently supported MOPAC binary as a subprocess.
:type inputfile: str
:param inputfile: name of .mop input file (without suffix)
"""
if not isinstance(self._tmpdir, str):
msg = f"\nCould not acquire a {MOPAC_MAIN} license.\n"
msg += f"If you bought a {MOPAC_MAIN} license, please contact "
msg += "Schrodinger at help@schrodinger.com\n"
try:
self._tmpdir = mm.mmjag_mopac_startup()
except MmException:
print(msg)
raise MopacLicenseError(msg)
if isinstance(self._tmpdir, str) and os.path.exists(self._tmpdir):
pass
else:
print(msg)
raise MopacLicenseError(msg)
os.environ["MOPAC_LICENSE"] = self._tmpdir
t0 = time.time()
rc = subprocess.call([self.MOPAC_EXEC, inputfile + '.mop'])
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, structures=()):
"""
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 (without the suffix)
:type structures: list
:param structures: list of Structure objects
:return: list(job status), outfile
"""
print(f"Starting semi-empirical calculation for input file {mopfile}.")
sys.stdout.flush()
# Run calculation with MOPAC external binary
rc = self._execute_mopac_binary(mopfile)
t0 = time.time()
status = []
full_outfile = mopfile + '.out'
if self._energy_only:
# Parse .out file only
if os.path.exists(full_outfile):
with open(full_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 ct, properties in zip_longest(
structures, all_results, fillvalue={}):
i += 1
if properties:
for k, v in properties.items():
if k[0] == 'mae':
ct.property[k[1]] = v
status.append(True)
else:
status.append(False)
if ct:
ct.property['s_j_jname'] = f'{mopfile}_{i}'
else:
# Parse .out and .aux files and populate structure properties
outfiles = split_outfile(full_outfile)
auxfiles = split_outfile(mopfile + '.aux')
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(f'The results for structure {idx} ({ct.title}) '
'was not found in Mopac .aux file. The '
'descriptors might be incomplete.')
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:
print(f'The results for structure {idx} ({ct.title}) was'
' not found in Mopac .out file. Skipping descriptors'
' 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)
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, mopfile + '.out'
#=============================================================================
# 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