Added a short recipe for calculating and rescaling ligand charges

nanoCAT/ff/ff_assignment.py

@@ -87,7 +87,9 @@ def start_ff_assignment(mol_list: Iterable[Molecule], jobs: Tuple[Type[Job], ...
         run_match_job(mol, s, job)
 
 
-def run_match_job(mol: Molecule, s: Settings, job_type: Type[Job] = MatchJob) -> None:
+def run_match_job(mol: Molecule, s: Settings,
+                  job_type: Type[Job] = MatchJob,
+                  action: str = 'warn') -> None:
     """Assign atom types and charges to **mol** based on the results of MATCH_.
 
     Performs an inplace update of :attr:`Atom.properties` ``["symbol"]``,
@@ -106,6 +108,10 @@ def run_match_job(mol: Molecule, s: Settings, job_type: Type[Job] = MatchJob) ->
     job_type : :class:`type` [|plams.Job|]
         The type of Job.
 
+    action : :class:`str`
+        The to-be undertaken action when the Job crashes.
+        Accepted values are ``"raise"``, ``"warn"`` and ``"ignore"``.
+
     See Also
     --------
     :class:`.MatchJob`
@@ -125,10 +131,13 @@ def run_match_job(mol: Molecule, s: Settings, job_type: Type[Job] = MatchJob) ->
         logger.info(f'{job.__class__.__name__}: {mol.properties.name} parameter assignment '
                     f'({job.name}) is successful')
     except Exception as ex:
-        logger.info(f'{job.__class__.__name__}: {mol.properties.name} parameter assignment '
-                    f'({job.name}) has failed')
-        logger.debug(f'{ex.__class__.__name__}: {ex}', exc_info=True)
-        return None
+        if action == 'raise':
+            raise ex
+        elif 'action' == 'warm':
+            logger.info(f'{job.__class__.__name__}: {mol.properties.name} parameter assignment '
+                        f'({job.name}) has failed')
+            logger.debug(f'{ex.__class__.__name__}: {ex}', exc_info=True)
+            return None
 
     # Update properties with new symbols, charges and the consntructed parameter (.prm) file
     mol.properties.prm = prm = results['$JN.prm']

nanoCAT/recipes/__init__.py

@@ -9,14 +9,19 @@
 .. code:: python
 
     >>> from nanoCAT.recipes import bulk_workflow
+    >>> from nanoCAT.recipes import get_lig_charge
     >>> from nanoCAT.recipes import replace_surface
     >>> from nanoCAT.recipes import dissociate_surface, row_accumulator
     ...
 
 """
 
 from .bulk import bulk_workflow
+from .charges import get_lig_charge
 from .mark_surface import replace_surface
 from .surface_dissociation import dissociate_surface, row_accumulator
 
-__all__ = ['bulk_workflow', 'replace_surface', 'dissociate_surface', 'row_accumulator']
+__all__ = [
+    'bulk_workflow', 'replace_surface', 'dissociate_surface',
+    'row_accumulator', 'get_lig_charge'
+]

nanoCAT/recipes/charges.py

@@ -0,0 +1,119 @@
+"""
+nanoCAT.recipes.charges
+=======================
+
+A short recipe for calculating and rescaling ligand charges.
+
+Index
+-----
+.. currentmodule:: nanoCAT.recipes.charges
+.. autosummary::
+    get_lig_charge
+
+API
+---
+.. autofunction:: get_lig_charge
+
+"""
+
+from os import PathLike
+from typing import Optional, Union
+
+import pandas as pd
+from scm.plams import Molecule, Atom, Settings, MoleculeError, JobError, init, finish
+
+from nanoCAT.ff import run_match_job
+
+__all__ = ['get_lig_charge']
+
+
+def get_lig_charge(ligand: Molecule,
+                   ligand_anchor: int,
+                   core_anchor_charge: float,
+                   settings: Optional[Settings] = None,
+                   path: Union[None, str, PathLike] = None,
+                   folder: Union[None, str, PathLike] = None) -> pd.Series:
+    """Calculate and rescale the **ligand** charges using MATCH_.
+
+    The atomic charge of **ligand_anchor**, as deterimined by MATCH, is rescaled such that
+    the molecular charge of **ligand** will be equal to **core_anchor_charge**.
+
+    .. _MATCH: http://brooks.chem.lsa.umich.edu/index.php?page=match&subdir=articles/resources/software
+
+    Examples
+    --------
+    .. code:: python
+
+        >>> import pandas as pd
+        >>> from scm.plams import Molecule
+
+        >>> from CAT.recipes import get_lig_charge
+
+        >>> ligand = Molecule(...)
+        >>> ligand_anchor = int(...)  # 1-based index
+        >>> core_anchor_charge = float(...)
+
+        >>> charge_series: pd.Series = get_lig_charge(
+        ...     ligand, ligand_anchor, core_anchor_charge
+        ... )
+
+        >>> charge_series.sum() == core_anchor_charge
+        True
+
+
+    Parameters
+    ----------
+    ligand : :class:`~scm.plams.core.mol.molecule.Molecule`
+        The input ligand.
+
+    ligand_anchor : :class:`int`
+        The (1-based) atomic index of the ligand anchor atom.
+        The charge of this atom will be rescaled.
+
+    core_anchor_charge : :class:`float`
+        The atomic charge of the core anchor atoms.
+
+    settings : :class:`~scm.plams.core.settings.Settings`, optional
+        The input settings for :class:`~nanoCAT.ff.match_job.MatchJob`.
+        Will default to the ``"top_all36_cgenff_new"`` forcefield if not specified.
+
+    path : :class:`str` or :class:`~os.PathLike`, optional
+        The path to the PLAMS workdir as passed to :func:`~scm.plams.core.functions.init`.
+        Will default to the current working directory if ``None``.
+
+    folder : :class:`str` or :class:`~os.PathLike`, optional
+        The name of the to-be created to the PLAMS working directory
+        as passed to :func:`~scm.plams.core.functions.init`.
+        Will default to ``"plams_workdir"`` if ``None``.
+
+    Returns
+    -------
+    :class:`pd.Series` [:class:`str`, :class:`float`]
+        A Series with the atom types of **ligand** as keys and atomic charges as values.
+
+    See Also
+    --------
+    :class:`MatchJob`
+        A :class:`~scm.plams.core.basejob.Job` subclass for interfacing with MATCH_:
+        Multipurpose Atom-Typer for CHARMM.
+
+    """  # noqa
+    if settings is None:
+        settings = Settings()
+        settings.input.forcefield = 'top_all36_cgenff_new'
+
+    # Run the MATCH Job
+    init(path, folder)
+    run_match_job(ligand, settings, action='raise')
+    finish()
+
+    # Extract the charges and new atom types
+    charge = [at.properties.charge_float for at in ligand]
+    symbol = [at.properties.symbol for at in ligand]
+
+    # Correct the charge of ligand[i]
+    i = ligand_anchor
+    charge[i - 1] -= sum(charge) - core_anchor_charge
+    ligand[i].properties.charge_float = charge[i - 1]
+
+    return pd.Series(charge, index=symbol, name='charge')