diff --git a/CAT/attachment/ligand_attach.py b/CAT/attachment/ligand_attach.py index d60c76f3..da587e6b 100644 --- a/CAT/attachment/ligand_attach.py +++ b/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 diff --git a/CAT/base.py b/CAT/base.py index c028b3d4..38600e36 100644 --- a/CAT/base.py +++ b/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 ' diff --git a/CAT/data_handling/validation_schemas.py b/CAT/data_handling/validation_schemas.py index 70e700f9..e9a75d22 100644 --- a/CAT/data_handling/validation_schemas.py +++ b/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. diff --git a/CAT/utils.py b/CAT/utils.py index 0c3e8a65..cb749c09 100644 --- a/CAT/utils.py +++ b/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 diff --git a/docs/4_optional.rst b/docs/4_optional.rst index 92fdb041..72c2b716 100644 --- a/docs/4_optional.rst +++ b/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. diff --git a/tests/test_files/CAT_allignment.yaml b/tests/test_files/CAT_allignment.yaml new file mode 100644 index 00000000..a93df542 --- /dev/null +++ b/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 diff --git a/tests/test_files/core/hollow.xyz b/tests/test_files/core/hollow.xyz new file mode 100644 index 00000000..5538d774 --- /dev/null +++ b/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 diff --git a/tests/test_files/test_allignment.hdf5 b/tests/test_files/test_allignment.hdf5 new file mode 100644 index 00000000..53b5a744 Binary files /dev/null and b/tests/test_files/test_allignment.hdf5 differ diff --git a/tests/test_ligand_attach.py b/tests/test_ligand_attach.py index b61e4335..548b9621 100644 --- a/tests/test_ligand_attach.py +++ b/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) diff --git a/tests/test_schemas.py b/tests/test_schemas.py index 11b7782e..bc97f76a 100644 --- a/tests/test_schemas.py +++ b/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: diff --git a/tests/test_validate_input.py b/tests/test_validate_input.py index ab25fa16..82eecb7e 100644 --- a/tests/test_validate_input.py +++ b/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')