ENH: Add an option for inverting the core vectors (#208)

CAT/attachment/ligand_attach.py

@@ -53,6 +53,7 @@
 from ..workflows import WorkFlow, HDF5_INDEX, MOL, OPT
 from ..settings_dataframe import SettingsDataFrame
 from ..data_handling import mol_to_file, WARN_MAP
+from ..utils import AllignmentTup, AllignmentEnum
 
 __all__ = ['init_qd_construction']
 
@@ -202,9 +203,13 @@ def _get_df(core_index: pd.MultiIndex,
     return SettingsDataFrame(data, index=index, columns=columns, settings=settings)
 
 
-def ligand_to_qd(core: Molecule, ligand: Molecule, path: str,
-                 allignment: str = 'sphere',
-                 idx_subset: Optional[Iterable[int]] = None) -> Molecule:
+def ligand_to_qd(
+    core: Molecule,
+    ligand: Molecule,
+    path: str,
+    allignment: AllignmentTup,
+    idx_subset: "None | Iterable[int]" = None,
+) -> Molecule:
     """Function that handles quantum dot (qd, *i.e.* core + all ligands) operations.
 
     Combine the core and ligands and assign properties to the quantom dot.
@@ -250,17 +255,16 @@ def get_name() -> str:
     ligand.properties.dummies.properties.anchor = True
 
     # Attach the rotated ligands to the core, returning the resulting strucutre (PLAMS Molecule).
-    if allignment == 'sphere':
-        vec2 = np.array(core.get_center_of_mass()) - sanitize_dim_2(core.properties.dummies)
+    anchor = sanitize_dim_2(core.properties.dummies)
+    if allignment.kind == AllignmentEnum.SPHERE:
+        vec2 = np.array(core.get_center_of_mass()) - anchor
         vec2 /= np.linalg.norm(vec2, axis=1)[..., None]
-    elif allignment == 'surface':
-        if isinstance(core.properties.dummies, np.ndarray):
-            anchor = core.properties.dummies
-        else:
-            anchor = core.as_array(core.properties.dummies)
+    elif allignment.kind == AllignmentEnum.SURFACE:
         vec2 = -get_surface_vec(np.array(core)[idx_subset_], anchor)
     else:
-        raise ValueError(repr(allignment))
+        raise ValueError(f"Unknown allignment kind: {allignment.kind}")
+    if allignment.invert:
+        vec2 *= -1
 
     # Rotate the ligands
     anchor_tup = ligand.properties.anchor_tup

CAT/base.py

@@ -35,6 +35,7 @@
 from .logger import logger
 from .mol_utils import to_symbol
 from .settings_dataframe import SettingsDataFrame
+from .utils import AllignmentEnum
 
 from .data_handling.mol_import import read_mol
 from .data_handling.update_qd_df import update_qd_df
@@ -236,7 +237,10 @@ def prep_core(core_df: SettingsDataFrame) -> SettingsDataFrame:
             raise MoleculeError(f"{repr(to_symbol(anchor))} was specified as core anchor atom, yet "
                                 f"no matching atoms were found in {core.properties.name} "
                                 f"(formula: {formula})")
-        elif len(dummies) < 4 and core_df.settings.optional.core.allignment == "surface":
+        elif (
+            len(dummies) < 4 and
+            core_df.settings.optional.core.allignment.kind == AllignmentEnum.SURFACE
+        ):
             raise NotImplementedError(
                 '`optional.core.allignment = "surface"` is not supported for cores with less '
                 f'than 4 anchor atoms ({core.get_formula()}); consider using '

CAT/data_handling/validation_schemas.py

@@ -89,6 +89,7 @@
 from .str_to_func import str_to_func
 from ..utils import get_template, validate_path, validate_core_atom, check_sys_var
 from ..mol_utils import to_atnum
+from ..utils import AllignmentEnum, AllignmentTup
 
 __all__ = ['mol_schema', 'core_schema', 'ligand_schema', 'qd_schema', 'database_schema',
            'mongodb_schema', 'bde_schema', 'ligand_opt_schema', 'qd_opt_schema', 'crs_schema',
@@ -287,6 +288,15 @@ def _get_crsjob() -> type:
         ),
 })
 
+allignment_mapping = {
+    'sphere': AllignmentTup(AllignmentEnum.SPHERE, False),
+    'surface': AllignmentTup(AllignmentEnum.SURFACE, False),
+    'sphere_invert': AllignmentTup(AllignmentEnum.SPHERE, True),
+    'sphere invert': AllignmentTup(AllignmentEnum.SPHERE, True),
+    'surface_invert': AllignmentTup(AllignmentEnum.SURFACE, True),
+    'surface invert': AllignmentTup(AllignmentEnum.SURFACE, True),
+}
+
 #: Schema for validating the ``['optional']['core']`` block.
 core_schema: Schema = Schema({
     'dirname':
@@ -312,9 +322,13 @@ def _get_crsjob() -> type:
     Optional_('subset', default=None):
         Or(None, dict, error="optional.core.subset epected 'None' or a dictionary"),
 
-    Optional_('allignment', default='surface'):
-        And(str, lambda n: n.lower() in {'sphere', 'surface'},
-            Use(str.lower), error="optional.core.allignment expected 'sphere' or 'surface'")
+    Optional_('allignment', default=AllignmentTup(AllignmentEnum.SURFACE, False)):
+        Or(
+            AllignmentTup,
+            And(str, Use(lambda n: allignment_mapping[n.lower()]),
+                error="optional.core.allignment expected "
+                      f"one of {list(allignment_mapping.keys())}"),
+        ),
 })
 
 #: Schema for validating the ``['optional']['core']['subset']`` block.

CAT/utils.py

@@ -64,7 +64,8 @@
 __all__ = [
     'JOB_MAP', 'check_sys_var', 'dict_concatenate', 'get_template',
     'cycle_accumulate', 'iter_repeat', 'get_nearest_neighbors',
-    'log_traceback_locals', 'KindEnum', 'AnchorTup',
+    'log_traceback_locals', 'KindEnum', 'AnchorTup', 'AllignmentEnum',
+    'AllignmentTup'
 ]
 
 JOB_MAP: Mapping[Type[Job], str] = MappingProxyType({
@@ -541,6 +542,13 @@ class KindEnum(enum.Enum):
     MEAN_TRANSLATE = 2
 
 
+class AllignmentEnum(enum.Enum):
+    """An enum with different core vector orientations (see :class:`AllignmentTup`)."""
+
+    SPHERE = 0
+    SURFACE = 1
+
+
 class AnchorTup(NamedTuple):
     """A named tuple with anchoring operation instructions."""
 
@@ -552,3 +560,10 @@ class AnchorTup(NamedTuple):
     kind: KindEnum = KindEnum.FIRST
     angle_offset: "None | float" = None
     dihedral: "None | float" = None
+
+
+class AllignmentTup(NamedTuple):
+    """A named tuple with core vector orientation options."""
+
+    kind: AllignmentEnum
+    invert: bool

docs/4_optional.rst

@@ -274,13 +274,17 @@ Core
 
         How the to-be attached ligands should be alligned with the core.
 
-        Has two allowed values:
+        Has four allowed values:
 
         * ``"surface"``: Define the core vectors as those orthogonal to the cores
           surface. Not this option requires at least four core anchor atoms.
           The surface is herein defined by a convex hull constructed from the core.
         * ``"sphere"``: Define the core vectors as those drawn from the core anchor
           atoms to the cores center.
+        * ``"surface invert"``/``"surface_invert"``: The same as ``"surface"``,
+          except the core vectors are inverted.
+        * ``"sphere invert"``/``"sphere_invert"``: The same as ``"sphere"``,
+          except the core vectors are inverted.
 
         Note that for a spherical core both approaches are equivalent.
 

tests/test_files/CAT_allignment.yaml

@@ -0,0 +1,20 @@
+path: .
+
+input_cores:
+    - hollow.xyz
+
+input_ligands:
+    - 'OCC'
+
+optional:
+    core:
+        dirname: core
+        anchor: Cl
+
+    ligand:
+        dirname: ligand
+        optimize: True
+        anchor: O[H]
+
+    qd:
+        dirname: QD

tests/test_files/core/hollow.xyz

@@ -0,0 +1,35 @@
+33
+ i =       47, E =     -4043.8485700335
+        Cl      -0.910491       7.214308       6.693878
+        Cl       1.590784       8.564860      -3.865445
+        Cl      -4.307037       0.009034      -8.420637
+        Cl       1.536902      -2.388272       9.101635
+        Cl      -5.993994      -6.274559      -4.734885
+        Cl       9.505902      -1.851327       1.294392
+        Cl      -7.679948       2.143190       4.355035
+        Cl       2.613744       9.025486       1.862328
+        Cl      -0.856856      -8.723727       3.735762
+        Cl       2.464013       3.714335       8.410423
+        Cl      -0.757143       9.412411      -0.456071
+        Cl      -2.100715      -3.454485      -8.689662
+        Cl      -4.361711      -8.298856       1.457025
+        Cl      -8.869807      -1.002493       1.688519
+        Cl       4.778033       0.002975       8.310429
+        Cl       8.360727      -4.888691      -1.281903
+        Cl       1.197377      -9.501254       0.268370
+        Cl       8.414037      -2.144278      -4.536689
+        Cl       4.915729       8.145936      -1.558279
+        Cl      -8.865258       1.562422      -1.352194
+        Cl      -2.305031      -8.941942      -2.291994
+        Cl      -0.806542       2.226697      -9.239422
+        Cl       1.379456      -1.202354      -9.470958
+        Cl       9.549569       0.864406      -1.983007
+        Cl      -0.815242       1.080078       9.393042
+        Cl      -7.676493       4.706107       1.310313
+		He       0.0            0.0            0.0
+		He       1.0            0.0            0.0
+		He      -1.0            0.0            0.0
+		He       0.0            1.0            0.0
+		He       0.0           -1.0            0.0
+		He       0.0            0.0            1.0
+		He       0.0            0.0           -1.0

tests/test_ligand_attach.py

@@ -1,12 +1,12 @@
 """Tests for :mod:`CAT.attachment.ligand_attach`."""
 
-import sys
 import shutil
 from pathlib import Path
 from typing import Generator, NamedTuple, TYPE_CHECKING
 
 import yaml
 import pytest
+import h5py
 import numpy as np
 from assertionlib import assertion
 from scm.plams import Settings, Molecule
@@ -77,6 +77,14 @@ def test_get_rotmat2() -> None:
 class DihedTup(NamedTuple):
     mol: Molecule
     ref: np.recarray
+    name: str
+
+
+class AllignmentTup(NamedTuple):
+    mol: Molecule
+    atoms_ref: np.recarray
+    bonds_ref: np.recarray
+    name: str
 
 
 class TestDihedral:
@@ -100,12 +108,12 @@ def run_cat(self, request: "_pytest.fixtures.SubRequest") -> Generator[DihedTup,
         qd = qd_df[MOL].iloc[0]
 
         ref = np.load(PATH / f"test_dihedral_{name}.npy").view(np.recarray)
-        yield DihedTup(qd, ref)
+        yield DihedTup(qd, ref, name)
 
         # Teardown
         files = [LIG_PATH, QD_PATH, DB_PATH]
-        for f in files:
-            shutil.rmtree(f, ignore_errors=True)
+        for file in files:
+            shutil.rmtree(file, ignore_errors=True)
 
     def test_atoms(self, output: DihedTup) -> None:
         dtype = [("symbols", "U2"), ("coords", "f8", 3)]
@@ -116,3 +124,55 @@ def test_atoms(self, output: DihedTup) -> None:
         assertion.eq(atoms.dtype, output.ref.dtype)
         np.testing.assert_array_equal(atoms.symbols, output.ref.symbols)
         np.testing.assert_allclose(atoms.coords, output.ref.coords)
+
+
+class TestAllignment:
+    PARAMS = ("sphere", "surface", "sphere_invert", "surface_invert")
+
+    @pytest.fixture(scope="class", name="output", params=PARAMS)
+    def run_cat(
+        self, request: "_pytest.fixtures.SubRequest"
+    ) -> Generator[AllignmentTup, None, None]:
+        # Setup
+        allignment: str = request.param
+        yaml_path = PATH / 'CAT_allignment.yaml'
+        with open(yaml_path, 'r') as f1:
+            arg = Settings(yaml.load(f1, Loader=yaml.FullLoader))
+
+        arg.path = PATH
+        arg.optional.core.allignment = allignment
+        qd_df, _, _ = prep(arg)
+        qd = qd_df[MOL].iloc[0]
+
+        with h5py.File(PATH / "test_allignment.hdf5", "r") as f2:
+            atoms_ref = f2[f"TestAllignment/{allignment}/atoms"][...].view(np.recarray)
+            bonds_ref = f2[f"TestAllignment/{allignment}/bonds"][...].view(np.recarray)
+        yield AllignmentTup(qd, atoms_ref, bonds_ref, allignment)
+
+        # Teardown
+        files = [LIG_PATH, QD_PATH, DB_PATH]
+        for file in files:
+            shutil.rmtree(file, ignore_errors=True)
+
+    def test_atoms(self, output: AllignmentTup) -> None:
+        dtype = [("symbols", "S2"), ("coords", "f8", 3)]
+        iterator = ((at.symbol, at.coords) for at in output.mol)
+        atoms = np.fromiter(iterator, dtype=dtype).view(np.recarray)
+
+        assertion.eq(atoms.dtype, output.atoms_ref.dtype)
+        np.testing.assert_array_equal(atoms.symbols, output.atoms_ref.symbols)
+        np.testing.assert_allclose(atoms.coords, output.atoms_ref.coords)
+
+    def test_bonds(self, output: AllignmentTup) -> None:
+        dtype = [("atom1", "i8"), ("atom2", "i8"), ("order", "f8")]
+        try:
+            output.mol.set_atoms_id()
+            iterator = ((b.atom1.id, b.atom2.id, b.order) for b in output.mol.bonds)
+            bonds = np.fromiter(iterator, dtype=dtype).view(np.recarray)
+        finally:
+            output.mol.unset_atoms_id()
+
+        assertion.eq(bonds.dtype, output.bonds_ref.dtype)
+        np.testing.assert_array_equal(bonds.atom1, output.bonds_ref.atom1)
+        np.testing.assert_array_equal(bonds.atom2, output.bonds_ref.atom2)
+        np.testing.assert_allclose(bonds.order, output.bonds_ref.order)

tests/test_schemas.py

@@ -9,7 +9,7 @@
 from scm.plams import AMSJob, ADFJob, Settings, CRSJob
 from assertionlib import assertion
 
-from CAT.utils import get_template
+from CAT.utils import get_template, AllignmentTup, AllignmentEnum
 from CAT.data_handling.str_to_func import str_to_func
 from CAT.data_handling.validation_schemas import (
     mol_schema, core_schema, ligand_schema, qd_schema, database_schema,
@@ -156,7 +156,7 @@ def test_core_schema() -> None:
         'dirname': '.',
         'anchor': None,
         'dummy': None,
-        'allignment': 'surface',
+        'allignment': AllignmentTup(AllignmentEnum.SURFACE, False),
         'subset': None
     }
 
@@ -175,10 +175,24 @@ def test_core_schema() -> None:
     assertion.assert_(core_schema.validate, core_dict, exception=SchemaError)
     core_dict['allignment'] = 'bob'  # Exception: incorrect value
     assertion.assert_(core_schema.validate, core_dict, exception=SchemaError)
+
     core_dict['allignment'] = 'SPHERE'
-    assertion.eq(core_schema.validate(core_dict)['allignment'], 'sphere')
+    assertion.eq(
+        core_schema.validate(core_dict)['allignment'],
+        AllignmentTup(AllignmentEnum.SPHERE, False),
+    )
+
     core_dict['allignment'] = 'surface'
-    assertion.eq(core_schema.validate(core_dict)['allignment'], 'surface')
+    assertion.eq(
+        core_schema.validate(core_dict)['allignment'],
+        AllignmentTup(AllignmentEnum.SURFACE, False),
+    )
+
+    core_dict['allignment'] = 'surface invert'
+    assertion.eq(
+        core_schema.validate(core_dict)['allignment'],
+        AllignmentTup(AllignmentEnum.SURFACE, True),
+    )
 
 
 def test_qd_schema() -> None:

tests/test_validate_input.py

@@ -14,6 +14,7 @@
 from nanoutils import delete_finally
 
 from CAT.data_handling.validate_input import validate_input
+from CAT.utils import AllignmentTup, AllignmentEnum
 
 try:
     from dataCAT import Database
@@ -46,7 +47,7 @@ def test_validate_input() -> None:
     ref = Settings()
     ref.core.dirname = join(PATH, 'core')
     ref.core.anchor = 35
-    ref.core.allignment = 'surface'
+    ref.core.allignment = AllignmentTup(AllignmentEnum.SURFACE, False)
     ref.core.subset = None
 
     ref.database.dirname = join(PATH, 'database')