From 0df7a5d44fc0bfadd894a1573325af0c2cd41526 Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Fri, 19 Nov 2021 14:08:22 +0100 Subject: [PATCH 1/5] ENH: Add the `dihedral` option --- CAT/data_handling/anchor_parsing.py | 10 +++++++++- CAT/utils.py | 1 + 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/CAT/data_handling/anchor_parsing.py b/CAT/data_handling/anchor_parsing.py index 9eec4ba2..9f67bb94 100644 --- a/CAT/data_handling/anchor_parsing.py +++ b/CAT/data_handling/anchor_parsing.py @@ -80,7 +80,7 @@ def _parse_kind(typ: "None | str | KindEnum") -> KindEnum: def _parse_angle_offset( offset: "None | SupportsFloat | SupportsIndex | bytes | str" ) -> "None | float": - """Parse the ``angle_offset`` option; convert the offset to radians.""" + """Parse the ``angle_offset`` and ``dihedral`` options; convert the offset to radians.""" if offset is None: return None elif not isinstance(offset, str): @@ -103,6 +103,7 @@ def _parse_angle_offset( Optional("remove", default=None): Use(_parse_remove), Optional("kind", default=KindEnum.FIRST): Use(_parse_kind), Optional("angle_offset", default=None): Use(_parse_angle_offset), + Optional("dihedral", default=None): Use(_parse_angle_offset) }) @@ -162,6 +163,13 @@ def parse_anchors( raise ValueError("`group_idx` must contain at least 3 atoms when " "`angle_offset` is specified") + # Check that at least 2 atoms are available for `dihedral` + # (so the third dihedral-defining vector can be defined) + dihedral = kwargs["dihedral"] + if dihedral is not None and len(group_idx) < 2: + raise ValueError("`group_idx` must contain at least 3 atoms when " + "`dihedral` is specified") + mol = _smiles_to_rdmol(group) ret.append(AnchorTup(**kwargs, mol=mol)) return tuple(ret) diff --git a/CAT/utils.py b/CAT/utils.py index 765e8005..0c3e8a65 100644 --- a/CAT/utils.py +++ b/CAT/utils.py @@ -551,3 +551,4 @@ class AnchorTup(NamedTuple): remove: "None | Tuple[int, ...]" = None kind: KindEnum = KindEnum.FIRST angle_offset: "None | float" = None + dihedral: "None | float" = None From f97fb744ce3e11fbd0cda890e9e82cdad54fe63e Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Fri, 19 Nov 2021 15:14:29 +0100 Subject: [PATCH 2/5] ENH: Implement the `dihedral` option --- CAT/attachment/ligand_attach.py | 60 +++++++++++++++++++++++++++++---- 1 file changed, 53 insertions(+), 7 deletions(-) diff --git a/CAT/attachment/ligand_attach.py b/CAT/attachment/ligand_attach.py index e4247729..d60c76f3 100644 --- a/CAT/attachment/ligand_attach.py +++ b/CAT/attachment/ligand_attach.py @@ -49,6 +49,7 @@ from .perp_surface import get_surface_vec from ..mol_utils import get_index, round_coords # noqa: F401 +from ..utils import AnchorTup from ..workflows import WorkFlow, HDF5_INDEX, MOL, OPT from ..settings_dataframe import SettingsDataFrame from ..data_handling import mol_to_file, WARN_MAP @@ -261,7 +262,19 @@ def get_name() -> str: else: raise ValueError(repr(allignment)) - lig_array = rot_mol(ligand, vec1, vec2, atoms_other=core.properties.dummies, core=core, idx=idx) + # Rotate the ligands + anchor_tup = ligand.properties.anchor_tup + if anchor_tup.dihedral is None: + lig_array = rot_mol( + ligand, vec1, vec2, atoms_other=core.properties.dummies, core=core, idx=idx + ) + else: + lig_array = rot_mol( + ligand, vec1, vec2, atoms_other=core.properties.dummies, core=core, + idx=idx, anchor_tup=anchor_tup, + ) + + # Combine the ligands and core qd = core.copy() array_to_qd(ligand, lig_array, mol_out=qd) qd.round_coords() @@ -281,6 +294,18 @@ def get_name() -> str: return qd +def _get_dihedrals(mat1: np.ndarray, mat2: np.ndarray, vec3: np.ndarray) -> np.ndarray: + """Get the dihedral angle between three vectors in radian.""" + v1v2 = np.cross(-mat1, mat2) + v2v3 = np.cross(vec3, mat2) + v1v2_v2v3 = np.cross(v1v2, v2v3) + v2_norm_v2 = mat2 / np.linalg.norm(mat2, axis=1)[..., None] + return np.arctan2( + np.einsum("ij,ij->i", v1v2_v2v3, v2_norm_v2), + np.einsum("ij,ij->i", v1v2, v2v3), + ) + + def _get_rotmat1(vec1: np.ndarray, vec2: np.ndarray) -> np.ndarray: """Calculate the rotation matrix for rotating **vec1** to **vec2**. @@ -345,7 +370,11 @@ def _get_rotmat1(vec1: np.ndarray, vec2: np.ndarray) -> np.ndarray: return ret -def _get_rotmat2(vec: np.ndarray, step: float = (1/16)) -> np.ndarray: +def _get_rotmat2( + vec: np.ndarray, + step: float = 1/16, + angle_vec: "None | np.ndarray" = None, +) -> np.ndarray: r"""Calculate the rotation matrix for rotating m vectors along their axes, each vector yielding :math:`k = (2 / step)` possible rotations. Paramaters @@ -367,10 +396,13 @@ def _get_rotmat2(vec: np.ndarray, step: float = (1/16)) -> np.ndarray: [v3, zero, -v1], [-v2, v1, zero]]).T - step_range = np.pi * np.arange(0.0, 2.0, step) - a1 = np.sin(step_range)[:, None, None, None] - a2 = (1 - np.cos(step_range))[:, None, None, None] - + if angle_vec is None: + step_range = np.pi * np.arange(0.0, 2.0, step) + a1 = np.sin(step_range)[:, None, None, None] + a2 = 1 - np.cos(step_range)[:, None, None, None] + else: + a1 = np.sin(angle_vec)[:, None, None] + a2 = 1 - np.cos(angle_vec)[:, None, None] return np.identity(3) + a1 * w + a2 * w@w @@ -383,7 +415,8 @@ def rot_mol(xyz_array: np.ndarray, bond_length: Optional[int] = None, step: float = 1/16, dist_to_self: bool = True, - ret_min_dist: bool = False) -> np.ndarray: + ret_min_dist: bool = False, + anchor_tup: "None | AnchorTup" = None) -> np.ndarray: r"""Rotate **xyz_array**. Paramaters @@ -441,6 +474,19 @@ def rot_mol(xyz_array: np.ndarray, rotmat1 = _get_rotmat1(vec1, vec2) xyz = xyz@rotmat1 + # Code-path for fixed-angle dihedrals + if anchor_tup is not None: + i, j, *_ = anchor_tup.anchor_idx + vec3 = xyz[:, i] - xyz[:, j] + dihedral_vec = _get_dihedrals(vec3, vec2, vec1) + dihedral_vec -= anchor_tup.dihedral + rotmat2 = _get_rotmat2(vec2, angle_vec=dihedral_vec) + xyz = np.matmul(xyz, rotmat2, out=xyz) + + at_other = sanitize_dim_2(atoms_other) + xyz += (at_other - xyz[:, idx])[:, None] + return xyz + # Create all k possible rotations of all m ligands rotmat2 = _get_rotmat2(vec2, step=step) xyz = np.swapaxes(xyz@rotmat2, 0, 1) From e2ea3a0d1ca3ce72730cb0ff997451abd3fac404 Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Fri, 19 Nov 2021 15:20:56 +0100 Subject: [PATCH 3/5] MAINT: Simplify the `group_idx`/`remove` out-of-bounds check --- CAT/data_handling/anchor_parsing.py | 34 +++++++++++------------------ 1 file changed, 13 insertions(+), 21 deletions(-) diff --git a/CAT/data_handling/anchor_parsing.py b/CAT/data_handling/anchor_parsing.py index 9f67bb94..f07ccce8 100644 --- a/CAT/data_handling/anchor_parsing.py +++ b/CAT/data_handling/anchor_parsing.py @@ -5,7 +5,7 @@ from typing import Union, Tuple, Collection, Iterable, SupportsFloat from rdkit.Chem import Mol -from scm.plams import Units, PT +from scm.plams import Units from schema import Schema, Use, Optional from typing_extensions import TypedDict, SupportsIndex @@ -33,15 +33,6 @@ class _AnchorDict(TypedDict): angle_offset: "None | float" -def _pt_sort_func(symbol: str) -> Tuple[int, str]: - return len(symbol), symbol - - -_SYMBOL = "|".join(sorted(PT.symtonum, key=_pt_sort_func)) -_SYMBOL_PATTERN = re.compile(f"({_SYMBOL})") -_UNIT_PATTERN = re.compile(r"([\.\_0-9]+)(\s+)?(\w+)?") - - def _parse_group_idx(item: "SupportsIndex | Iterable[SupportsIndex]") -> Tuple[int, ...]: """Parse the ``group_idx`` option.""" try: @@ -77,6 +68,9 @@ def _parse_kind(typ: "None | str | KindEnum") -> KindEnum: return KindEnum[typ.upper()] +_UNIT_PATTERN = re.compile(r"([\.\_0-9]+)(\s+)?(\w+)?") + + def _parse_angle_offset( offset: "None | SupportsFloat | SupportsIndex | bytes | str" ) -> "None | float": @@ -146,16 +140,6 @@ def parse_anchors( if remove is not None and not set(group_idx).isdisjoint(remove): raise ValueError("`group_idx` and `remove` must be disjoint") - # Check that the indices in `group_idx` and `remove` are not out of bounds - group = kwargs["group"] - symbols_list = _SYMBOL_PATTERN.findall(group) - if len(symbols_list) <= max(group_idx): - raise IndexError(f"`group_idx` index {max(group_idx)} is out of bounds " - f"for a `group` with {len(symbols_list)} atoms") - elif remove is not None and len(symbols_list) <= max(remove): - raise IndexError(f"`remove` index {max(remove)} is out of bounds " - f"for a `group` with {len(symbols_list)} atoms") - # Check that at least 3 atoms are available for `angle_offset` # (so a plane can be defined) angle_offset = kwargs["angle_offset"] @@ -170,6 +154,14 @@ def parse_anchors( raise ValueError("`group_idx` must contain at least 3 atoms when " "`dihedral` is specified") - mol = _smiles_to_rdmol(group) + # Check that the indices in `group_idx` and `remove` are not out of bounds + mol = _smiles_to_rdmol(kwargs["group"]) + atom_count = len(mol.GetAtoms()) + if atom_count <= max(group_idx): + raise IndexError(f"`group_idx` index {max(group_idx)} is out of bounds " + f"for a `group` with {atom_count} atoms") + elif remove is not None and atom_count <= max(remove): + raise IndexError(f"`remove` index {max(remove)} is out of bounds " + f"for a `group` with {atom_count} atoms") ret.append(AnchorTup(**kwargs, mol=mol)) return tuple(ret) From 1c3204bd0dee7b9754165ab020e78d761f4fe8a1 Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Fri, 19 Nov 2021 17:08:49 +0100 Subject: [PATCH 4/5] DOC: Document the new `dihedral` option --- docs/4_optional.rst | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/docs/4_optional.rst b/docs/4_optional.rst index 111c9369..92fdb041 100644 --- a/docs/4_optional.rst +++ b/docs/4_optional.rst @@ -630,6 +630,7 @@ Ligand * :attr:`anchor.remove` * :attr:`anchor.kind` * :attr:`anchor.angle_offset` + * :attr:`anchor.dihedral` .. note:: @@ -716,6 +717,23 @@ Ligand but if so desired one can explicitly pass the unit: ``angle_offset: "0.25 rad"``. + .. attribute:: optional.ligand.anchor.dihedral + + :Parameter: * **Type** - :data:`None`, :class:`float` or :class:`str` + * **Default value** – :data:`None` + + Manually specify the ligands vector dihedral angle, rather than optimizing it w.r.t. the inter-ligand distance. + + The dihedral angle is defined by three vectors: + + * The first two in dices in :attr:`anchor.group_idx `. + * The core vector(s). + * The Cartesian X-axis as defined by the core. + + By default the angle unit is assumed to be in degrees, + but if so desired one can explicitly pass the unit: ``dihedral: "0.5 rad"``. + + .. attribute:: optional.ligand.split :Parameter: * **Type** - :class:`bool` From 7c17563531079a45f53cfd661f9f8ac2328a48c2 Mon Sep 17 00:00:00 2001 From: Bas van Beek Date: Fri, 19 Nov 2021 17:40:21 +0100 Subject: [PATCH 5/5] TST: Add tests for the new `dihedral` option --- tests/test_files/CAT_dihedral.yaml | 24 ++++++ tests/test_files/test_dihedral_dihed_180.npy | Bin 0 -> 9888 bytes tests/test_files/test_dihedral_dihed_45.npy | Bin 0 -> 9888 bytes .../test_files/test_dihedral_dihed_45_deg.npy | Bin 0 -> 9888 bytes tests/test_ligand_attach.py | 71 ++++++++++++++++-- 5 files changed, 90 insertions(+), 5 deletions(-) create mode 100644 tests/test_files/CAT_dihedral.yaml create mode 100644 tests/test_files/test_dihedral_dihed_180.npy create mode 100644 tests/test_files/test_dihedral_dihed_45.npy create mode 100644 tests/test_files/test_dihedral_dihed_45_deg.npy diff --git a/tests/test_files/CAT_dihedral.yaml b/tests/test_files/CAT_dihedral.yaml new file mode 100644 index 00000000..a8378063 --- /dev/null +++ b/tests/test_files/CAT_dihedral.yaml @@ -0,0 +1,24 @@ +path: . + +input_cores: + - Cd68Se55.xyz + +input_ligands: + - 'CC(O)=O' + +optional: + core: + dirname: core + anchor: Cl + + ligand: + dirname: ligand + optimize: True + anchor: + group: "[H]OC=O" + group_idx: [1, 2, 3] + remove: 0 + dihedral: 45 + + qd: + dirname: QD diff --git a/tests/test_files/test_dihedral_dihed_180.npy b/tests/test_files/test_dihedral_dihed_180.npy new file mode 100644 index 0000000000000000000000000000000000000000..db315cf6180aaef6477497867b747e9003c46eab GIT binary patch literal 9888 zcmZWv36zx8m2SXbXc4;`#D;F@s_NR;|NpD1yQ&)FRR`Hd1w=wrU{Km|L7FWZL}Xl; zj1pp$QGz0h8Y1J8gcuTuhbTs)#EiQRNf=jj99INIjBw_?``vo?WbAV|9{cP6ulw%u z-S6I4Yfd}!)XCrJ*KldWvc%lp#pf+b)Iw(`6N{H#a&F(e#ffYvG3Lx5k;;aWiSzpU z7R|-4&aW6GS$493kG(IPdTvvcdR{wy_1XuswDSEazZw06;x8I>#?k@}$!)p)VdhPs{$32P-yz1tYkIFcF z@dY=2c*TkW&Du6$aR2+FG`PR{gms?7t=UJslBSkx9^JX|i41*U%v<@s;^y>`zbuUt zoA<9i{pH}57dHf8d4OIKszms#Y@5v(F zy3YJm6=!=s^~hRvO+yBL9#?elH{Mz`@1JuN{`n%~DUZI|Yy(!LGfeDz&A(?4*n?#d z|Jez;bEJv*DM^R_=^69$gE;Jp7pRjRq_~ahldIy$n_P-)cvgK-aL@Dm&-MvzW5jsu>hsn&tdwu+ij zbulJHpPTaDakIwH($fjw?Zmkb#KnJsv!C;&p#yk#9`Wg~IP%*k!^*iWqjJRgw7bnc z&-uPS4_;UoNS(Jk+@{~|Xh10&$a@pSd|d5xIG<-s8oz(Ih&(<`{l|IL8h~Ffj!1qf zzj6KlF^2kIU6MRh{ziRE{dTsKZ6A$mHU3<;wFq&1UnvoX`M;Iln0K$uQ=4r>{X)d` z__XrdLPHC8dujSdn_jUv`qbQ#JKslFoKT*%Ek>L+2eaU@O^Tx*MUbDb`^4oqM{(vY zuS^^#Q`g|bhVNL z7d8L-YypK+#CiB=33aen{eQ1vpX=s$^4=KiIK9h#KG)6OEHTesng##vQ2YtQsu$dt zq)#l3P4LN(oD-pEF9tJkwH^KA<@GtxM`U4?DWlKYYH& zO6V_-2a+d>b6$-Ngtxkg?QK^7@%a{Z3ICMt^QgTkkS#N%0G<|7g+mbV)=wnA}e`P%r7u`5127cKda_{GQou8o*HV}d6 zkey+Nvwz0Lp$k?9qQeviFMFt)Vh8GfzUmj|x68xCJUqAP^0ThLxZhotAL=E1xd?ReL=XJb&vW*H&*SKCV>?kF7plI)=Y~ zb3+k4ZTPPqAD$~RqMt8`5cj+HwLUmL4$^X6V(z)jBYTRT*^0Bi z_^?E$aAWhNa7ey<9Mb=a7a^bf_E?_VxK*PNY^KAv%(&%6@xz`Om+GKZ=ESU2r1 zVeWY-Mz)Ra@w9$;{<=ZlyQ@sxubx$$pEvvk9hN{Jd9CC=pYtN=f!DU>C4bd_eE!uL z)Wx;{_9YbOdFk9P%q8y>F=u|H_0RF0*+s1L&D_TR{hiha=7J*n_A@2ujL#M4xnp9N z*q6op(608gF5O=wo?8z4l7EUrfBC5Id&-zkmO(#Q%>Z6WfqzF=r2f(0EY9)Wn1MdK zpd0#bc^92-LG#0N<*W?&XKWYd%U_lqj(ICb+-F|tMjbxYMdR^!1J7$;C7{nP>86A& zq+7ck&Smwt1o51;q=L^MrO$JJH+hZuvI-tw+s*yOegxN>nMa}DYly=V=xd9!?)bn% z|9+|hKfN7u_H*78Wd6CMTXaC&;XGI61L8UIw-q<1)X&HBe5N3}qbB^7aOZp8j+^O% z=<1rpsng+n?)@Hk>)jgqYgX;&`Z*~gyq_cFYqL8)=f(D@@Lvw`d2yJtpZn69C_dlJ zo$UX!Mmik2Hi7y0xg5@UTZ^*~dyz2mXF~xz?kf&H3u8YrAWO@w!=F|h_C>&p>vF;~ zjn00a3ru~oewv!Xp5oS;!@0lh_K5rRnJMW*RcAl+p(p&F!W`6IaX9xEqlY*@r%N4t ztLo0je0o$tcr8IZ4?a+HIM-)m0sS$Y0ss6=an?s~d6?g(q=@H%1wGDw*7pZSK>z-r zgLofw3-r6qPway_2t3^Ze%Kmu_VfHw3`>6r%X^{@$9^L$`nR1}Z)d~Ke%7n<`JNo1|FyCAIu|KHE+>w+$81$*d`TfhQvUf_u z&Z(ZWAM-;6`=<-jm>*A79Q&>s`ofR01U~SbeZ1E@vWB{SI)OfWeeP>G_D*r+`2u~O z=e=1y;E#V#3XiG%$meESVHc3!^M!vff4K9_JkHH;tmi%|y7OV~Z9`u#`rx;jigVv6 zwSu=^3NS}MTypmF{1IzL{4Vt1$2W^#!?}KLlK!(Ia5(e-mZQ-Z+5)mRH-xl5aDEN- zWxlInfAA~JzxGc#FLujZd0JI;jK&B0-$#8OTS1>P^K-p_p#F-;>*-a=L-jxUwuks{ zuAzR{`_6vMryg|Asv7#^`hvr`FFfiauNL(X>*b%S{m^&9DfhI-t*U ze_WPETs}Jz`o-*n?S9<;aDBe*A#aBl!NYG=9nN{N#Rso_7(mC~qx{JJHTz%AvnK-V z*)HyO=ZF3XVE-q9==jRt`=317FPd{%L1F5#-XGaWgEJ2ts_P)-U*2E;(3kos%AQpD zANnYRzP7I>dot|5?7LBKS()3a=)a$5o&CIjd>|wBTM_=tIh^Oqg$2?7RoJmF>+I+L z+dnOP;U4JAJMs=k|4fT6uSp-wIGpv#x-4|@t{&lcwV&sUC*$Csm*S$&n_T=jKdLD_ ze=v^v-!atL&-HmIDSMX$>^x+_xL(Kj_a`>_9rp>N!ulsaf~IPaf+6hr@M zj$uFdLw%n0`s4(iZ+|w$p(B&uboN7EMP$Bhk^QgMCyyPQqww1`$I2c<`GfoCWiq!l zwxK@E98>2P+$)IxrVT|OyIuJO`ZI*S`$~(<$I2hjpCQbZ=MBXk;9BJ$&W~S2Wbbtp z^!KmTe(Zw>$-PWg?x7UNy;ckOy_CZo5%-<_?6-G@i|)_jJzMmA=npN>jnM+;{C#R4 z@4IIVLVx}^1%CTN&tJ!}XU)odtmo(Pv~d9LrOcek^Ku0HV_RRizbVRddGNq<$}g-x z+DhyA9~Yl`v<=BF}vDVoRm7AwEt-Y-D?&d8!p_b7khJ|d90 zsv!MQaqjz(5_Cx>51#vp@(b`Pbg$7xZ0`!}kL~&4zY5|yKZ^a$l^GX5zK+`8<@O64 zuO*5jKPv5M_c|Mzf^13-(4M%I#&B} zuORmmhX%?0r|KK%U(w$$G)evE-T9zb3z$!9!_h}pD~|gxANzrGhM*5NX?|h;5dNz+ zLPyV0ea8CkrY!h#(qQoOGRi3j7 z|GjLkb^P%)$;+hhklGLb73Kb<3Vl9Uo4ya8sv`b>(){E3@Z$iyxVDP@@)*ru?k`hH=wC}~a*vPu z1M5=WZ|@l-bDxL!Ou+q@#o4}Si|~^#^R@PO%umCK=gbWr^!;S5Kd!r;7Mv> zxZkVucLDbduj~2hIPU=-^q>!irksCxt{gr<_C`MDKTrF6eSgb&IQAIP*Lr?@uOKOV zr52gLQ!akIA24$X?jL0i>r(HIQNSaJlSh!M29vx`}sa=Ujcc$H6!0+818VckA?O9T}I~cAibrj_#jnH{=c!u@wM_~*zHc>H?RpV;>whrYa^g!%Ce)t8u$+T?p3MbT%9 zbN&`v(Vu$D;NkmJfAalDs#)ff0D1PR_U}5*d203!e2?653N*{5o8)#dS7vlpJfAYN`c+lwz;a^ZQs^A^U_E%7m@hw((Z#fzWc z+q+;Eel(}cSiF2XkvP}<#+=>-ix$l6nPJSc?epVj&RaO!{Cwfune%6x|9g4Am^QCi z+VcOt1|DBkRW-}}yZ(DW-1p)dk4hW=uy{*zidGLlV%ngO+o@vD?QQFS)vF#|@YJH! z!*{gPz}ME!AA3XE*?)ELPm|9Z+(zARowDkVN7JU-;cz9UQam_0)h&d%rl#HVI^@WUb5 z<9)SS?fm{yo@{vUJ1arkzwpkv>FS(&ewqzTZ9pfioVWkGlYIKzgyTSe^^fB_B?JGS z$!Ss6p!WD;c4v+{2__dLJyB+0VIY+ZW!RM)cjJEc$DN=j@1Q1buN!D}83+V7S)z zh341?|I|p&PULGrk)z0S08NxXe#OL|%QSSYf=S&t&3;A>$H`Hin=XmvIXt&AY zp|MfvQ^gnV`{hye|MEbdS9}}B2~X4aZ2(_tLLHSAUo0%0QqxR_ORL*P|l_9RD zy^Zd9uJ?Pgz=gX)>GLLMx8=9DEu@4^VU6c7Ov{E8ecn(UhnY{x79`T?=rPNXU4+F-XiM% zv^es7uJ#xE|7{xhadU^%x5l5(HKGXI=#3E9TQTeWTfzn6k&VJhM|T3pKGXQ~eE%X8 z{BNVY{d!sJuh%wEYl65Aw-(R`2Q~f=8}zw;rZ4YphaabOxX4pX zG^jfFW{*C#Fm_Tvwq1%_)Xwd>yn`%^Z+tkw?_IC);rjeRff#3o6lD&m9r)HEeVv6L z9U31#-_r%mm+hg{iQ2iY#)X1g9mM|DY5e$n7j+2!C4l428o!B_kZdZTZWBK6;~uqh zUoQ=b)3YW|obNZ(&hsVeW3E2ZA>Y-``70Js|Fct42SpcuuIJVeaXsu~4(?OCrA?_7 zIhn5oobRQgn@^ZqGS?qzkvvA6oq6MoHewuID>!&X)Y+Mr&}hLoAM@c(wX^@XwV|$$ zh@yUflXTBBpP1K-^X*H>In~a5^It~Ow}!SYI<*aR>`0xj%x7Z48z;2`U$(Tk_jAAY zq-caqL?}FDSF5vgd`^sk7c38jhp8R7?4xh;&FKFg~{9gF! z8wu%;(eC|>Z^s3QkKvh^m(9-3eK{dR_gH%LL4lZW#x=go54Poj`)9P_9J|`x^Ne4c zJ>lm?)aN*jAM^4D^JG|E(^<1dksZgY&QN~Ey!k*Bb!+gE`Tv>9@9cP(KQ;ir+|@#i z$Agr=0Us*5ArG83_*cmf>xz`{=L;jm^X?<<56;i8JK+qWn> zOykGAX>S3#XJb3rF}ly!{$c%flf3u!BJsR>LGAp!!7uQzIOfP(1^4+}7oHDXdonBa ztMTLWuS%gWo($n%TKL8n-w`C?r;Jq7$3*8zR`MA6xyk37Wh?hes+oiqWDSMl8bWs;U! z0&rA``1fBaIy=9AbOv#~u?sw|MxW=o{9HnKX$kmndneBq`w`rhjiT>fED_`7v*3@` z4*emzv#&hAIOhE4y2*+Dxvx`rfZ93#+p@^-ttHX_p7Wpa=L4UhV@slMBF@fxhj{$Z~o?|(T-uNoN=TvzApyf28wfPYgnz|ZBi z&W?OU2y-?=%$Mg6bpA7*RR^NKl9+FE%FfPwY_t#j-{2ANQx=z;|KNwd;Cli(sHw}@ zp^vlRuLonmp?6EpKgPYk_`tQ})55>Y&R)?UtgC*LK>s$Vo%6fjhhNhY*dz6JJO4R9 z*E9jQ9&3{M3jS{Mi@nY$(PgdZ?_p8rKkJu#tIU^H@aUhlIXlnS2SyU};671_v(_Ku zWv&Tx<_a=*P*vvvdCcwS@UT*MmKiasIKcm{*p%_K=^& zfwQxJyTpEX3VpE6cXs90paQIhR)7<_O*KOqJDYd@damR zJ-elbtglrqA?ACl&IjgiMmJ*bR0dBtL-WVDX7H7D*}@X^$%w$65Bxts-5t|~I2iq0 znIGV1KI-$G65ex2^TYV@Od$NSj6LhJtn>dRyRI9U2Y*VIfxjy?e+*|k0_;yN>IOew zr~WhlzBYxvTU(1bFV*~k_axDuxq7!>tN73Lxg;g{G6ME3p*tVX7sJ08C%T5>^LOg` zSpTL|sLwYBgP-51`Qdfu#ti1#ra{1|shWS@w=GLb{~n3>9IyT}Kin<+ctIZXbz{ll z5995p!t3r2v47sJ_{jBV_P?C}ogw1;<4$)z<`+#q=JWFbbdu@k%6#Ye-5FqhU0cBX zU7>d7w|hd&ukJkN-%r)fd(B&7e!Nn0 zbEc9%{(zw+(Cr|1mnIl3$(=(_4vg;=^OWe}AX= zf&F(B`@_@_*c<-|KZnI+ZVZOb*`xTv{l7UPy1g2E&drK1%txkGgSUR1K_0$PJMJ6m zB!4O3+=uEvkHz;Jg!iYh$JwBE+z-@(cTGs3eugT(V1HE&T z%wdDKmHx$jQa_1*M)-;52lqk+(brkw`HC;B3(v^|->=GmpFE=Yg87<0&vUVKJVFXQ9JLqrWJ*kX62qh?Hu1&0r-{8f}gHaeBu54 z>=Neg+$i>pZ)5*%>y!0dbC>K#qw+oJ>W};9^A5LPU>&$9ecotn;7i`}Q34X}N*;zGRi!!5`D29|i-TZ&m-GPb!~Zjrmrm^@sYG`O-5$ z_Cs2~e7?Iye{Jj|c~*Xg`=%0f>C2w%(N+JkznLzHUjp@Yx$rwfw6k27Pg|FCzKx$=-l%opZiKX3B~ zeO3T(pILG^Ucpc3i#%}g4)Jra)*tV`uFgw;cjErxeXU=<&)poNZ__2gDa{}DZw2gq zjh*Mwn7H$w_5DwaqFYMP!Efq(WdBbuVDEZq+3l$+{`33KtH$1CMIieD?LVBa4t;aJ z+%s?0^Wi>S`1id&?s=!E9rL|b?vJwI=M8#3==N&xM{~c%dhtE=ANX4X+$)EO<09=( z@S}dnTROl!M5yzf_3{7@y6UDH@Z&4BKXLyWm%Z3fxySHa|M9-$){3qf4xZYZaCYXM zm$n1ndWMS5(ENa3#?dzeYSE_;>-jN%lfpl%<(@t1^2__m=VNle*I(|pasOlMlk>MG z1%0&Fm-DH<$NnJ++`KgqKZZE}@%IWT_%$H_&l)<+*|ER$P#23r;g@yJKg2IZ*yH77 zPh0Kmyia^FBYP#`uRD%%cD_HjFA09u72xkAM)q_3g!jou`oPb>Cc%@}DZhb!%Hs1U zryY(~_zCN#Lph0WO5)q#&d2+etjx`$v)~6eDE|Wg%;E0?9`{85sQ=6_jP3*vO#$bQ z7~sxVxi7-pOvyb+A9p^UpSPEVpJs(WXnt8gOeqQP%MpHG>yz(4+PmcW9O6);_`>}; zsf79XS`Ob0->UErzOOg9fWLD{<9z1tYbti$f4)){zkS?qU#xcEWjFHs>#WoV?%!;E z@_g!V`qj?IXXT!7gU)}xcb}RPU7C~o72PlKKE&*wu_u+kkC~(VhWXL@gy@E>;G^z` zct2(MKl7wL0qXoZwKH!`=48L)qmP!V9eXpmH@;SM^wYXuV*hG~%DsC5@!PHQ9e>}@ zfP05W^Qf<1DZhk1YQ$b=dtUUX+L_N3>hb)Ip}QZi653N*{5o8)#dS7vlpJfAYN`c+lwz;a^ZQs^A^U_E%7m@hw((Z#fzWc z+q+;Eel(}cSiF2XkvP}<#+=>-ix$l6nPJSc?epVj&RaO!{Cwfune%6x|9g4Am^QCi z+VcOt1|DBkRW-}}yZ(DW-1p)dk4hW=uy{*zidGLlV%ngO+o@vD?QQFS)vF#|@YJH! z!*{gPz}ME!AA3XE*?)ELPm|9Z+(zARowDkVN7JU-;cz9UQam_0)h&d%rl#HVI^@WUb5 z<9)SS?fm{yo@{vUJ1arkzwpkv>FS(&ewqzTZ9pfioVWkGlYIKzgyTSe^^fB_B?JGS z$!Ss6p!WD;c4v+{2__dLJyB+0VIY+ZW!RM)cjJEc$DN=j@1Q1buN!D}83+V7S)z zh341?|I|p&PULGrk)z0S08NxXe#OL|%QSSYf=S&t&3;A>$H`Hin=XmvIXt&AY zp|MfvQ^gnV`{hye|MEbdS9}}B2~X4aZ2(_tLLHSAUo0%0QqxR_ORL*P|l_9RD zy^Zd9uJ?Pgz=gX)>GLLMx8=9DEu@4^VU6c7Ov{E8ecn(UhnY{x79`T?=rPNXU4+F-XiM% zv^es7uJ#xE|7{xhadU^%x5l5(HKGXI=#3E9TQTeWTfzn6k&VJhM|T3pKGXQ~eE%X8 z{BNVY{d!sJuh%wEYl65Aw-(R`2Q~f=8}zw;rZ4YphaabOxX4pX zG^jfFW{*C#Fm_Tvwq1%_)Xwd>yn`%^Z+tkw?_IC);rjeRff#3o6lD&m9r)HEeVv6L z9U31#-_r%mm+hg{iQ2iY#)X1g9mM|DY5e$n7j+2!C4l428o!B_kZdZTZWBK6;~uqh zUoQ=b)3YW|obNZ(&hsVeW3E2ZA>Y-``70Js|Fct42SpcuuIJVeaXsu~4(?OCrA?_7 zIhn5oobRQgn@^ZqGS?qzkvvA6oq6MoHewuID>!&X)Y+Mr&}hLoAM@c(wX^@XwV|$$ zh@yUflXTBBpP1K-^X*H>In~a5^It~Ow}!SYI<*aR>`0xj%x7Z48z;2`U$(Tk_jAAY zq-caqL?}FDSF5vgd`^sk7c38jhp8R7?4xh;&FKFg~{9gF! z8wu%;(eC|>Z^s3QkKvh^m(9-3eK{dR_gH%LL4lZW#x=go54Poj`)9P_9J|`x^Ne4c zJ>lm?)aN*jAM^4D^JG|E(^<1dksZgY&QN~Ey!k*Bb!+gE`Tv>9@9cP(KQ;ir+|@#i z$Agr=0Us*5ArG83_*cmf>xz`{=L;jm^X?<<56;i8JK+qWn> zOykGAX>S3#XJb3rF}ly!{$c%flf3u!BJsR>LGAp!!7uQzIOfP(1^4+}7oHDXdonBa ztMTLWuS%gWo($n%TKL8n-w`C?r;Jq7$3*8zR`MA6xyk37Wh?hes+oiqWDSMl8bWs;U! z0&rA``1fBaIy=9AbOv#~u?sw|MxW=o{9HnKX$kmndneBq`w`rhjiT>fED_`7v*3@` z4*emzv#&hAIOhE4y2*+Dxvx`rfZ93#+p@^-ttHX_p7Wpa=L4UhV@slMBF@fxhj{$Z~o?|(T-uNoN=TvzApyf28wfPYgnz|ZBi z&W?OU2y-?=%$Mg6bpA7*RR^NKl9+FE%FfPwY_t#j-{2ANQx=z;|KNwd;Cli(sHw}@ zp^vlRuLonmp?6EpKgPYk_`tQ})55>Y&R)?UtgC*LK>s$Vo%6fjhhNhY*dz6JJO4R9 z*E9jQ9&3{M3jS{Mi@nY$(PgdZ?_p8rKkJu#tIU^H@aUhlIXlnS2SyU};671_v(_Ku zWv&Tx<_a=*P*vvvdCcwS@UT*MmKiasIKcm{*p%_K=^& zfwQxJyTpEX3VpE6cXs90paQIhR)7<_O*KOqJDYd@damR zJ-elbtglrqA?ACl&IjgiMmJ*bR0dBtL-WVDX7H7D*}@X^$%w$65Bxts-5t|~I2iq0 znIGV1KI-$G65ex2^TYV@Od$NSj6LhJtn>dRyRI9U2Y*VIfxjy?e+*|k0_;yN>IOew zr~WhlzBYxvTU(1bFV*~k_axDuxq7!>tN73Lxg;g{G6ME3p*tVX7sJ08C%T5>^LOg` zSpTL|sLwYBgP-51`Qdfu#ti1#ra{1|shWS@w=GLb{~n3>9IyT}Kin<+ctIZXbz{ll z5995p!t3r2v47sJ_{jBV_P?C}ogw1;<4$)z<`+#q=JWFbbdu@k%6#Ye-5FqhU0cBX zU7>d7w|hd&ukJkN-%r)fd(B&7e!Nn0 zbEc9%{(zw+(Cr|1mnIl3$(=(_4vg;=^OWe}AX= zf&F(B`@_@_*c<-|KZnI+ZVZOb*`xTv{l7UPy1g2E&drK1%txkGgSUR1K_0$PJMJ6m zB!4O3+=uEvkHz;Jg!iYh$JwBE+z-@(cTGs3eugT(V1HE&T z%wdDKmHx$jQa_1*M)-;52lqk+(brkw`HC;B3(v^|->=GmpFE=Yg87<0&vUVKJVFXQ9JLqrWJ*kX62qh?Hu1&0r-{8f}gHaeBu54 z>=Neg+$i>pZ)5*%>y!0dbC>K#qw+oJ>W};9^A5LPU>&$9ecotn;7i`}Q34X}N*;zGRi!!5`D29|i-TZ&m-GPb!~Zjrmrm^@sYG`O-5$ z_Cs2~e7?Iye{Jj|c~*Xg`=%0f>C2w%(N+JkznLzHUjp@Yx$rwfw6k27Pg|FCzKx$=-l%opZiKX3B~ zeO3T(pILG^Ucpc3i#%}g4)Jra)*tV`uFgw;cjErxeXU=<&)poNZ__2gDa{}DZw2gq zjh*Mwn7H$w_5DwaqFYMP!Efq(WdBbuVDEZq+3l$+{`33KtH$1CMIieD?LVBa4t;aJ z+%s?0^Wi>S`1id&?s=!E9rL|b?vJwI=M8#3==N&xM{~c%dhtE=ANX4X+$)EO<09=( z@S}dnTROl!M5yzf_3{7@y6UDH@Z&4BKXLyWm%Z3fxySHa|M9-$){3qf4xZYZaCYXM zm$n1ndWMS5(ENa3#?dzeYSE_;>-jN%lfpl%<(@t1^2__m=VNle*I(|pasOlMlk>MG z1%0&Fm-DH<$NnJ++`KgqKZZE}@%IWT_%$H_&l)<+*|ER$P#23r;g@yJKg2IZ*yH77 zPh0Kmyia^FBYP#`uRD%%cD_HjFA09u72xkAM)q_3g!jou`oPb>Cc%@}DZhb!%Hs1U zryY(~_zCN#Lph0WO5)q#&d2+etjx`$v)~6eDE|Wg%;E0?9`{85sQ=6_jP3*vO#$bQ z7~sxVxi7-pOvyb+A9p^UpSPEVpJs(WXnt8gOeqQP%MpHG>yz(4+PmcW9O6);_`>}; zsf79XS`Ob0->UErzOOg9fWLD{<9z1tYbti$f4)){zkS?qU#xcEWjFHs>#WoV?%!;E z@_g!V`qj?IXXT!7gU)}xcb}RPU7C~o72PlKKE&*wu_u+kkC~(VhWXL@gy@E>;G^z` zct2(MKl7wL0qXoZwKH!`=48L)qmP!V9eXpmH@;SM^wYXuV*hG~%DsC5@!PHQ9e>}@ zfP05W^Qf<1DZhk1YQ$b=dtUUX+L_N3>hb)Ip}QZi None: @@ -50,11 +64,58 @@ def test_get_rotmat2() -> None: vec1 = np.array([1, 0, 0], dtype=float) vec2 = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=float) - ref1 = np.load(join(PATH, 'rotmat2_1.npy')) - ref2 = np.load(join(PATH, 'rotmat2_2.npy')) + ref1 = np.load(PATH / 'rotmat2_1.npy') + ref2 = np.load(PATH / 'rotmat2_2.npy') rotmat1 = _get_rotmat2(vec1) np.testing.assert_allclose(rotmat1, ref1) rotmat2 = _get_rotmat2(vec2) np.testing.assert_allclose(rotmat2, ref2) + + +class DihedTup(NamedTuple): + mol: Molecule + ref: np.recarray + + +class TestDihedral: + PARAMS = { + "dihed_45": 45, + "dihed_45_deg": "45 deg", + "dihed_180": 180.0, + } + + @pytest.fixture(scope="class", autouse=True, name="output", params=PARAMS.items(), ids=PARAMS) + def run_cat(self, request: "_pytest.fixtures.SubRequest") -> Generator[DihedTup, None, None]: + # Setup + name, dihed = request.param # type: str, str | float + yaml_path = PATH / 'CAT_dihedral.yaml' + with open(yaml_path, 'r') as f: + arg = Settings(yaml.load(f, Loader=yaml.FullLoader)) + + arg.path = PATH + arg.optional.ligand.anchor.dihedral = dihed + qd_df, _, _ = prep(arg) + qd = qd_df[MOL].iloc[0] + + ref = np.load(PATH / f"test_dihedral_{name}.npy").view(np.recarray) + yield DihedTup(qd, ref) + + # Teardown + files = [LIG_PATH, QD_PATH, DB_PATH] + for f in files: + shutil.rmtree(f, ignore_errors=True) + + def test_atoms(self, output: DihedTup) -> None: + dtype = [("symbols", "U2"), ("coords", "f8", 3)] + atoms = np.fromiter( + [(at.symbol, at.coords) for at in output.mol], dtype=dtype + ).view(np.recarray) + + assertion.eq(atoms.dtype, output.ref.dtype) + np.testing.assert_array_equal(atoms.symbols, output.ref.symbols) + + if sys.version_info >= (3, 9): + pytest.xfail("Geometries must be updated for RDKit >2019.09.2") + np.testing.assert_allclose(atoms.coords, output.ref.coords)