import os
from pathlib import Path
from typing import Dict, List, Optional, Tuple, Union
from networkx import DiGraph
import scm.reactions_discovery.utils
from scm.base import KFError, KFFile
from scm.inputs import ReactionsDiscovery
from scm.plams import AMSJob, Molecule, Results, Settings, SingleJob
from scm.plams.core.errors import FileError
from scm.reactions_discovery.combined_ct_results import CombinedMol, CombinedReaction
from scm.reactions_discovery.reaction_network import ReactionNetwork
from scm.utils.conversions import chemsys_to_plams_molecule
[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
scm.inputs model.
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[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.inputs.ReactionsDiscovery
scm.inputs model 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 a scm.inputs model. 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 = ReactionsDiscovery.from_input(text_input)
else:
self.settings.input = 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) -> ReactionsDiscovery:
"""scm.inputs model input"""
return self.settings.input
@input.setter
def input(self, input: ReactionsDiscovery):
self.settings.input = input
@property
def builder(self) -> ReactionsDiscovery.MolecularDynamicsBlock.BuildSystemBlock:
return self.input.MolecularDynamics.BuildSystem
@builder.setter
def builder(self, builder: ReactionsDiscovery.MolecularDynamicsBlock.BuildSystemBlock):
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)
def _serialize_single_molecule(self, name, mol):
return AMSJob._serialize_single_molecule(self, name, mol)
@staticmethod
def _atom_symbol(s):
return AMSJob._atom_symbol(s)
@staticmethod
def _atom_suffix(s):
return AMSJob._atom_suffix(s)
