from typing import List, Optional, Union, Dict, Tuple
import os
from pathlib import Path
from networkx import DiGraph
from scm.plams import Results, AMSJob, SingleJob, Settings, Molecule
from scm.plams.core.errors import FileError
from scm.input_classes import drivers
from scm.libbase import KFFile, KFError
import scm.reactions_discovery.utils
from scm.reactions_discovery.reaction_network import ReactionNetwork
from scm.reactions_discovery.combined_ct_results import CombinedMol, CombinedReaction
[docs]class ReactionsDiscoveryResults(Results):
"""
Results class for ReactionsDiscoveryJob
"""
[docs] def get_errormsg(self) -> str:
"""Returns the error message of this calculation if any were raised.
:return: String containing the error message.
:rtype: str
"""
return self.job.get_errormsg()
[docs] def get_md_jobs(self) -> List[AMSJob]:
"""Returns the AMSJobs used during the calculation.
:return: List of AMSJobs used during the calculation.
:rtype: List[AMSJob]
"""
return self.job.get_md_jobs()
[docs] def rkfpath(self) -> str:
"""Returns path to reactions_discovery.rkf
:return: Path to reactions_discovery.rkf
:rtype: str
"""
return str((Path(self.job.path) / "reactions_discovery.rkf").resolve())
[docs] def get_network_rd(self) -> Tuple[DiGraph, Dict[str, List[CombinedMol]], Dict[int, CombinedReaction], List[str]]:
"""Returns the reaction network represented by Reactions Discovery CombinedMol and CombinedReaction classes.
:raises KFError: If the KF file associated with this result does not contain the required information.
:return: Graph representing the reaction network, a dictionairy of categories and lists of CombinedMol, a
dictionairy of categories and CombinedReaction and a list of categories.
:rtype: Tuple[DiGraph, Dict[str, List[CombinedMol]], Dict[int, CombinedReaction], List[str]]
"""
reaction_network = ReactionNetwork()
reaction_network.read_extraction_results(self.job.path)
reaction_network.molecules_to_graph()
molecules = reaction_network.molecules
graph = reaction_network.graph
assert isinstance(graph, DiGraph)
with KFFile(self.rkfpath()) as kf:
try:
categories = [category for category in kf.get_skeleton().get("Categories", []) if category[:3] != "Num"]
assert isinstance(molecules, list)
rd_molecules: Dict[str, List[CombinedMol]] = {}
reactions: Dict[int, CombinedReaction] = {}
for category in categories:
mol_idxs = kf.read_ints_np("Categories", category)
rd_molecules[category] = [molecules[idx - 1] for idx in mol_idxs]
for molecule in molecules:
for property in ["smiles", "name", "net_flux", "string_description", "molecular_formula", "id"]:
molecule.plams_mol.properties[property] = getattr(molecule, property)
for reaction in molecule.reactions_as_product + molecule.reactions_as_reactant:
reactions[reaction.hash] = reaction
except KFError:
raise KFError(f"Could not read reaction network from {self.rkfpath()}")
return graph, rd_molecules, reactions, categories
[docs] def get_network(self) -> Tuple[DiGraph, Dict[str, List[Molecule]], List[str]]:
"""Returns the reaction network represented by a DiGraph and a dictionairy of lists of PLAMS molecules.
Each key in the dictionary is a category.
:return: graph of the reaction network, dictionary of categories and lists of Molecules, and a list of
categories.
:rtype: Tuple[DiGraph, Dict[str, List[Molecule]], List[str]]
"""
graph, molecules, reactions, categories = self.get_network_rd()
molecules = {category: [molecule.plams_mol for molecule in molecules[category]] for category in molecules}
return graph, molecules, categories
[docs] def get_num_md_simulations(self) -> int:
"""Returns the number of MD simulations used during the Molecular Dynamics stage.
:raises KFError: If the KF file associated with this result does not contain the right information.
:return: The number of MD simulations used during the Molecular Dynamics stage.
:rtype: int
"""
with KFFile(self.rkfpath()) as kf:
try:
njobs = kf.read_int("MolecularDynamicsResults", "NumSimulations")
except KFError:
raise KFError(f'"MolecularDynamicsResults%NumSimulations" not found on {self.rkfpath()}') from None
return njobs
def recreate_molecule(self) -> Dict[str, Molecule]:
"""Obtain the input molecule(s) used to create this result from the rkf file associated with this result.
:return: The input molecule(s).
:rtype: Dict[str, Molecule]
"""
molecule = ReactionsDiscoveryJob.from_rkf(self.rkfpath()).molecule
if isinstance(molecule, Molecule):
molecule = {"": molecule}
return molecule
def recreate_settings(self) -> Settings:
"""Obtain the original input used to create this result from the rkf file associated with this result.
:return: The original input Settings.
:rtype: Settings
"""
return ReactionsDiscoveryJob.from_rkf(self.rkfpath()).settings
[docs]class ReactionsDiscoveryJob(SingleJob):
"""
PLAMS Job class for running Reactions Discovery.
This class inherits from the PLAMS SingleJob class. For usage, see the SingleJob documentation.
If you supply a Settings object to the constructor, it will be converted to a
PISA (Python Input System for AMS) object.
Attributes:
* ``input``: an alias for self.settings.input
* ``builder``: an alias for self.settings.input.MolecularDynamics.BuildSystem
"""
_result_type = ReactionsDiscoveryResults
results: ReactionsDiscoveryResults
_command = "reactions_discovery"
_json_definitions = "reactions_discovery"
_subblock_end = "End"
[docs] def __init__(
self,
name: str = "reactions_discovery_job",
driver: Optional[drivers.ReactionsDiscovery] = None,
settings: Optional[Settings] = None,
molecule: Optional[Union[Molecule, Dict[str, Molecule]]] = None,
**kwargs,
):
"""
Initialize the ReactionsDiscoveryJob.
name : str
The name of the job
driver : scm.input_classes.drivers.ReactionsDiscovery
PISA object describing the input to the ReactionsDiscovery program
settings: scm.plams.Settings
All settings for the job. Input settings in the PLAMS settings format under ``settings.input`` are
automatically converted to the PISA format. You can specify ``settings.runscript.nproc`` to
set the total number of cores to run on.
molecule: scm.plams.Molecule or Dict[str, scm.plams.Molecule]
Two possibilities:
* ``molecule`` is of type Molecule - it should then be the *complete* system as a PLAMS Molecule .
Cannot be combined with the ``driver.input.MolecularDynamics.BuildSystem`` or
``settings.input.ams.MolecularDynamics.BuildSystem``. It will be written to the
main System block in the input.
* ``molecule`` is a dictionary with string keys and Molecule values - the keys should then be given in
the ``driver.input.MolecularDynamics.BuildSystem.Molecule[i].SystemID`` input option. The molecules will
then be used to build the system before the MD.
"""
super().__init__(name=name, settings=settings, molecule=molecule, **kwargs)
if driver is not None:
self.settings.input = driver
elif self.settings.input:
text_input = AMSJob(settings=self.settings).get_input()
self.settings.input = drivers.ReactionsDiscovery.from_text(text_input)
else:
self.settings.input = drivers.ReactionsDiscovery()
[docs] @classmethod
def from_rkf(cls, path: str) -> "ReactionsDiscoveryJob":
"""Initialize a job from a reactions_discovery.rkf file.
:param path: Path to a reactions_discovery.rkf file
:type path: str
:return: A new ReactionsDiscoveryJob instance based on the information found in path.
:rtype: ReactionsDiscoveryJob
"""
with KFFile(path) as kf:
text_input = kf.read_string("General", "user input")
return cls.from_input(text_input)
[docs] def get_errormsg(self) -> str:
"""Returns the contents of the jobname.err file if it exists. If the file does not exist an
empty string is returned.
:return: The error message
:rtype: str
"""
try:
with open(self.results["$JN.err"], "r") as err:
errlines = err.read()
return errlines
except FileNotFoundError:
return ""
[docs] def get_runscript(self) -> str:
"""
Generates the runscript. Use ``self.settings.runscript.preamble_lines = ['line1', 'line2']``
or similarly for ``self.settings.runscript.postamble_lines`` to set custom settings.
``self.settings.runscript.nproc`` controls the total number of cores to run on.
"""
filename = self._filename("inp")
ret = ""
for line in self.settings.runscript.get("preamble_lines", ""):
ret += f"{line}\n"
# need to use `pwd` here and not "." since "." doesn't expand
ret += f'AMS_JOBNAME="{self.name}" AMS_RESULTSDIR="`pwd`" "$AMSBIN/{self._command}" '
nproc = self.settings.runscript.get("nproc", None)
if nproc:
ret += f"-n {nproc} "
ret += f'< "{filename}"\n'
for line in self.settings.runscript.get("postamble_lines", ""):
ret += f"{line}\n"
return ret
[docs] def check(self) -> bool:
"""Returns True if "NORMAL TERMINATION" is given in the General section of reactions_discovery.rkf,
AND all molecular dynamics jobs also have finished successfully.
"""
with KFFile(self.results.rkfpath()) as kf:
termination = kf.read_string("General", "termination status")
i_am_ok = "NORMAL TERMINATION" in termination
if not i_am_ok:
return False
md_ok = all(job.check() for job in self.get_md_jobs())
return md_ok
[docs] def ok(self) -> bool:
"""Synonym for check()"""
return self.check()
[docs] def get_md_jobs(self) -> List[AMSJob]:
"""Returns: List of AMSJob"""
# this is very ugly because the user may create new rkf files inside the
# directory and then those are counted....
results = scm.reactions_discovery.utils.load_results(self.path)
return [x.job for x in results]
@property
def input(self) -> drivers.ReactionsDiscovery:
"""PISA format input"""
return self.settings.input
@input.setter
def input(self, input: drivers.ReactionsDiscovery):
self.settings.input = input
@property
def builder(self) -> drivers.ReactionsDiscovery._MolecularDynamics._BuildSystem:
return self.input.MolecularDynamics.BuildSystem
@builder.setter
def builder(self, builder: drivers.ReactionsDiscovery._MolecularDynamics._BuildSystem):
self.input.MolecularDynamics.BuildSystem = builder
[docs] @classmethod
def load_external(cls, path: Union[str, Path], finalize: bool = False) -> "ReactionsDiscoveryJob":
"""Load a previous ReactionsDiscovery job from disk.
:param path: A reactions discovery results folder.
:type path: Union[str, Path]
:param finalize: See SingleJob, defaults to False
:type finalize: bool, optional
:raises FileError: When the path does not exist.
:return: An initialized ReactionsDiscoveryJob
:rtype: ReactionsDiscoveryJob
"""
path = Path(path)
if not os.path.isdir(path):
if os.path.exists(path):
path = os.path.dirname(os.path.abspath(path))
elif os.path.isdir(path / ".results"):
path = path / ".results"
elif os.path.isdir(path / "results"):
path = path / "results"
else:
raise FileError("Path {} does not exist, cannot load from it.".format(path))
job = super(ReactionsDiscoveryJob, cls).load_external(path, finalize=finalize)
if job.name.endswith(".results") and len(job.name) > 8:
job.name = job.name[:-8]
return job
def _serialize_input(self, s):
return AMSJob._serialize_input(self, s)
def _serialize_molecule(self):
return AMSJob._serialize_molecule(self)
@staticmethod
def _atom_symbol(s):
return AMSJob._atom_symbol(s)
@staticmethod
def _atom_suffix(s):
return AMSJob._atom_suffix(s)