SciCodePile/SciCode-Domain-Code

1"keyword","repo_name","file_path","file_extension","file_size","line_count","content","language"
2"Conformation","ntampellini/prism_pruner","prism_pruner/algebra.py",".py","4951","164","""""""Algebra utilities.""""""3 4from typing import Sequence5 6import numpy as np7 8from prism_pruner.typing import Array1D_float, Array2D_float, Array3D_float9 10 11def normalize(vec: Array1D_float) -> Array1D_float:12    """"""Normalize a vector.""""""13    return vec / np.linalg.norm(vec)14 15 16def vec_angle(v1: Array1D_float, v2: Array1D_float) -> float:17    """"""Return the planar angle defined by two 3D vectors.""""""18    return float(19        np.degrees(20            np.arccos(21                np.clip(22                    np.dot(23                        v1 / np.linalg.norm(v1),24                        v2 / np.linalg.norm(v2),25                    ),26                    -1.0,27                    1.0,28                ),29            )30        )31    )32 33 34def dihedral(p: Array2D_float) -> float:35    """"""36    Find dihedral angle in degrees from 4 3D vecs.37 38    Praxeolitic formula: 1 sqrt, 1 cross product.39    """"""40    p0, p1, p2, p3 = p41 42    b0 = -1.0 * (p1 - p0)43    b1 = p2 - p144    b2 = p3 - p245 46    # normalize b1 so that it does not influence magnitude of vector47    # rejections that come next48    b1 /= np.linalg.norm(b1)49 50    # vector rejections51    # v = projection of b0 onto plane perpendicular to b152    #   = b0 minus component that aligns with b153    # w = projection of b2 onto plane perpendicular to b154    #   = b2 minus component that aligns with b155    v = b0 - np.dot(b0, b1) * b156    w = b2 - np.dot(b2, b1) * b157 58    # angle between v and w in a plane is the torsion angle59    # v and w may not be normalized but that's fine since tan is y/x60    x = np.dot(v, w)61    y = np.dot(np.cross(b1, v), w)62 63    return float(np.degrees(np.arctan2(y, x)))64 65 66def rot_mat_from_pointer(pointer: Array1D_float, angle: float) -> Array2D_float:67    """"""68    Get the rotation matrix from the rotation pivot using a quaternion.69 70    :param pointer: 3D vector representing the rotation pivot71    :param angle: rotation angle in degrees72    :return rotation_matrix: matrix that applied to a point, rotates it along the pointer73    """"""74    assert pointer.shape[0] == 375 76    angle_2 = np.radians(angle) / 277    sin = np.sin(angle_2)78    pointer = pointer / np.linalg.norm(pointer)79    return quaternion_to_rotation_matrix(80        [81            sin * pointer[0],82            sin * pointer[1],83            sin * pointer[2],84            np.cos(angle_2),85        ]86    )87 88 89def quaternion_to_rotation_matrix(quat: Array1D_float | Sequence[float]) -> Array2D_float:90    """"""91    Convert a quaternion into a full three-dimensional rotation matrix.92 93    This rotation matrix converts a point in the local reference frame to a94    point in the global reference frame.95 96    :param quat: 4-element array representing the quaternion (q0, q1, q2, q3)97    :return: 3x3 element array representing the full 3D rotation matrix98    """"""99    # Extract the values from Q (adjusting for scalar last in input)100    q1, q2, q3, q0 = quat101 102    # First row of the rotation matrix103    r00 = 2 * (q0 * q0 + q1 * q1) - 1104    r01 = 2 * (q1 * q2 - q0 * q3)105    r02 = 2 * (q1 * q3 + q0 * q2)106 107    # Second row of the rotation matrix108    r10 = 2 * (q1 * q2 + q0 * q3)109    r11 = 2 * (q0 * q0 + q2 * q2) - 1110    r12 = 2 * (q2 * q3 - q0 * q1)111 112    # Third row of the rotation matrix113    r20 = 2 * (q1 * q3 - q0 * q2)114    r21 = 2 * (q2 * q3 + q0 * q1)115    r22 = 2 * (q0 * q0 + q3 * q3) - 1116 117    # 3x3 rotation matrix118    return np.array([[r00, r01, r02], [r10, r11, r12], [r20, r21, r22]])119 120 121def get_inertia_moments(coords: Array3D_float, masses: Array1D_float) -> Array1D_float:122    """"""Compute the principal moments of inertia of a molecule.123 124    Returns a length-3 array [I_x, I_y, I_z], sorted ascending.125    """"""126    # Shift to center of mass127    com = np.sum(coords * masses[:, np.newaxis], axis=0) / np.sum(masses)128    coords = coords - com129 130    # Compute inertia tensor131    norms_sq = np.einsum(""ni,ni->n"", coords, coords)132    total = np.sum(masses * norms_sq)133    I_matrix = total * np.eye(3) - np.einsum(""n,ni,nj->ij"", masses, coords, coords)134 135    # Principal moments via symmetric eigendecomposition136    moments, _ = np.linalg.eigh(I_matrix)137 138    return np.sort(moments)139 140 141def diagonalize(a: Array2D_float) -> Array2D_float:142    """"""Build the diagonalized matrix.""""""143    eigenvalues_of_a, eigenvectors_of_a = np.linalg.eig(a)144    b = eigenvectors_of_a[:, np.abs(eigenvalues_of_a).argsort()]145    return np.dot(np.linalg.inv(b), np.dot(a, b))  # type: ignore[no-any-return]146 147 148def get_alignment_matrix(p: Array1D_float, q: Array1D_float) -> Array2D_float:149    """"""150    Build the rotation matrix that aligns vectors q to p (Kabsch algorithm).151 152    Assumes centered vector sets (i.e. their mean is the origin).153    """"""154    # calculate the covariance matrix155    cov_mat = p.T @ q156 157    # Compute the SVD158    v, _, w = np.linalg.svd(cov_mat)159 160    # Ensure proper rotation (det = 1, not -1)161    if np.linalg.det(v) * np.linalg.det(w) < 0.0:162        v[:, -1] *= -1163 164    return v @ w  # type: ignore[no-any-return]165","Python"
166"Conformation","ntampellini/prism_pruner","prism_pruner/__init__.py",".py","55","2","""""""PRISM - Pruning Interface for Similar Molecules.""""""167","Python"
168"Conformation","ntampellini/prism_pruner","prism_pruner/conformer_ensemble.py",".py","1860","58","""""""ConformerEnsemble class.""""""169 170import re171from dataclasses import dataclass, field172from pathlib import Path173from typing import Self174 175import numpy as np176 177from prism_pruner.typing import Array1D_float, Array1D_str, Array2D_float, Array3D_float178 179 180@dataclass181class ConformerEnsemble:182    """"""Class representing a conformer ensemble.""""""183 184    coords: Array3D_float185    atoms: Array1D_str186    energies: Array1D_float = field(default_factory=lambda: np.array([]))187 188    @classmethod189    def from_xyz(cls, file: Path | str, read_energies: bool = False) -> Self:190        """"""Generate ensemble from a multiple conformer xyz file.""""""191        coords = []192        atoms = []193        energies = []194        with Path(file).open() as f:195            for num in f:196                if read_energies:197                    energy = next(re.finditer(r""-*\d+\.\d+"", next(f))).group()198                    energies.append(float(energy))199                else:200                    _comment = next(f)201 202                conf_atoms = []203                conf_coords = []204                for _ in range(int(num)):205                    atom, *xyz = next(f).split()206                    conf_atoms.append(atom)207                    conf_coords.append([float(x) for x in xyz])208 209                atoms.append(conf_atoms)210                coords.append(conf_coords)211 212        return cls(coords=np.array(coords), atoms=np.array(atoms[0]), energies=np.array(energies))213 214    def to_xyz(self, file: Path | str) -> None:215        """"""Write ensemble to an xyz file.""""""216 217        def to_xyz(coords: Array2D_float) -> str:218            return f""{len(coords)}\n\n"" + ""\n"".join(219                f""{atom} {x:15.8f} {y:15.8f} {z:15.8f}""220                for atom, (x, y, z) in zip(self.atoms, coords, strict=True)221            )222 223        with Path(file).open(""w"") as f:224            f.write(""\n"".join(map(to_xyz, self.coords)))225","Python"
226"Conformation","ntampellini/prism_pruner","prism_pruner/torsion_module.py",".py","17259","507","""""""PRISM - Pruning Interface for Similar Molecules.""""""227 228from copy import deepcopy229from dataclasses import dataclass230from typing import Callable, Iterable, Sequence231 232import numpy as np233from networkx import (234    Graph,235    connected_components,236    has_path,237    is_isomorphic,238    minimum_spanning_tree,239    shortest_path,240    subgraph,241)242 243from prism_pruner.algebra import vec_angle244from prism_pruner.graph_manipulations import (245    get_phenyl_ids,246    get_sp_n,247    is_amide_n,248    is_ester_o,249)250from prism_pruner.rmsd import rmsd_and_max251from prism_pruner.typing import Array1D_bool, Array1D_str, Array2D_float, Array2D_int252from prism_pruner.utils import rotate_dihedral253 254 255@dataclass256class Torsion:257    """"""Torsion class.""""""258 259    i1: int260    i2: int261    i3: int262    i4: int263    mode: str | None = None264 265    @property266    def torsion(self) -> tuple[int, int, int, int]:267        """"""Return tuple of indices defining the torsion.""""""268        return (self.i1, self.i2, self.i3, self.i4)269 270 271def in_cycle(torsion: Torsion, graph: Graph) -> bool:272    """"""Return True if the torsion is part of a cycle.""""""273    graph.remove_edge(torsion.i2, torsion.i3)274    cyclical: bool = has_path(graph, torsion.i1, torsion.i4)275    graph.add_edge(torsion.i2, torsion.i3)276    return cyclical277 278 279def is_rotable(280    torsion: Torsion,281    graph: Graph,282    hydrogen_bonds: list[list[int]],283    keepdummy: bool = False,284) -> bool:285    """"""Return True if the Torsion object is rotatable.286 287    hydrogen bonds: iterable with pairs of sorted atomic indices.288    """"""289    if sorted((torsion.i2, torsion.i3)) in hydrogen_bonds:290        # self.n_fold = 6291        # # This has to be an intermolecular HB: rotate it292        # return True293        return False294 295    if _is_free(torsion.i2, graph) or (_is_free(torsion.i3, graph)):296        if keepdummy or (297            is_nondummy(torsion.i2, torsion.i3, graph)298            and (is_nondummy(torsion.i3, torsion.i2, graph))299        ):300            return True301 302    return False303 304 305def get_n_fold(torsion: Torsion, graph: Graph) -> int:306    """"""Return the n-fold of the rotation.""""""307    atoms = (graph.nodes[torsion.i2][""atoms""], graph.nodes[torsion.i3][""atoms""])308 309    if ""H"" in atoms:310        return 6  # H-N, H-O hydrogen bonds311 312    if is_amide_n(torsion.i2, graph, mode=2) or (is_amide_n(torsion.i3, graph, mode=2)):313        # tertiary amides rotations are 2-fold314        return 2315 316    if (""C"" in atoms) or (""N"" in atoms) or (""S"" in atoms):  # if C, N or S atoms317        sp_n_i2 = get_sp_n(torsion.i2, graph)318        sp_n_i3 = get_sp_n(torsion.i3, graph)319 320        if 3 == sp_n_i2 == sp_n_i3:321            return 3322 323        if 3 in (sp_n_i2, sp_n_i3):  # Csp3-X, Nsp3-X, Ssulfone-X324            if torsion.mode == ""csearch"":325                return 3326 327            elif torsion.mode == ""symmetry"":328                return sp_n_i3 or 2329 330        if 2 in (sp_n_i2, sp_n_i3):331            return 2332 333    return 4  # O-O, S-S, Ar-Ar, Ar-CO, and everything else334 335 336def is_linear(torsion: Torsion, coords: Array2D_float, max_dev_deg: float = 5.0) -> bool:337    """"""338    Return wether three or more of the four atoms involved in the torsion are in line.339 340    :type torsion: Torsion341    :type coords: Array2D_float342    :type max_dev_deg: float343    :rtype: bool344    """"""345    p1, p2, p3, p4 = coords[list(torsion.torsion)]346    v21 = p1 - p2347    v23 = p3 - p2348    a1 = vec_angle(v21, v23)349 350    if abs(180 - a1) < max_dev_deg:351        return True352 353    v34 = p4 - p3354    a2 = vec_angle(v34, -v23)355 356    if abs(180 - a2) < max_dev_deg:357        return True358 359    return False360 361 362def get_angles(torsion: Torsion, graph: Graph) -> tuple[int, ...]:363    """"""Return the angles associated with the torsion.""""""364    d = {365        1: (0,),  # in case some sp carbons make it to here366        2: (0, 180),367        3: (0, 120, 240),368        4: (0, 90, 180, 270),369        6: (0, 60, 120, 180, 240, 300),370    }371 372    n_fold = get_n_fold(torsion, graph)373 374    return d[n_fold]375 376 377def _is_free(index: int, graph: Graph) -> bool:378    """"""Return whether the torsion is free to rotate.379 380    Return True if the index specified381    satisfies all of the following:382    - Is not a sp2 carbonyl carbon atom383    - Is not the oxygen atom of an ester384    - Is not the nitrogen atom of a secondary amide (CONHR)385    """"""386    if all(387        (388            graph.nodes[index][""atoms""] == ""C"",389            2 == get_sp_n(index, graph),390            ""O"" in (graph.nodes[n][""atoms""] for n in graph.neighbors(index)),391        )392    ):393        return False394 395    if is_amide_n(index, graph, mode=1):396        return False397 398    if is_ester_o(index, graph):399        return False400 401    return True402 403 404def is_nondummy(i: int, root: int, graph: Graph) -> bool:405    """"""Return whether the torsion is not dummy.406 407    Checks that a molecular rotation along the dihedral408    angle (*, root, i, *) is non-dummy, that is the atom409    at index i, in the direction opposite to the one leading410    to root, has different substituents. i.e. methyl, CF3 and tBu411    rotations should return False.412    """"""413    if graph.nodes[i][""atoms""] not in (""C"", ""N""):414        return True415    # for now, we only discard rotations around carbon416    # and nitrogen atoms, like methyl/tert-butyl/triphenyl417    # and flat symmetrical rings like phenyl, N-pyrrolyl...418 419    G = deepcopy(graph)420    nb = list(G.neighbors(i))421    nb.remove(root)422 423    if len(nb) == 1:424        if len(list(G.neighbors(nb[0]))) == 2:425            return False426    # if node i has two bonds only (one with root and one with a)427    # and the other atom (a) has two bonds only (one with i)428    # the rotation is considered dummy: some other rotation429    # will account for its freedom (i.e. alkynes, hydrogen bonds)430 431    # check if it is a phenyl-like rotation432    if len(nb) == 2:433        # get the 6 indices of the aromatic atoms (i1-i6)434        phenyl_indices = get_phenyl_ids(i, G)435 436        # compare the two halves of the 6-membered ring (indices i2-i3 region with i5-i6 region)437        if phenyl_indices is not None:438            i1, i2, i3, i4, i5, i6 = phenyl_indices439            G.remove_edge(i3, i4)440            G.remove_edge(i4, i5)441            G.remove_edge(i1, i2)442            G.remove_edge(i1, i6)443 444            subgraphs = [445                subgraph(G, _set) for _set in connected_components(G) if i2 in _set or i6 in _set446            ]447 448            if len(subgraphs) == 2:449                return not is_isomorphic(450                    subgraphs[0],451                    subgraphs[1],452                    node_match=lambda n1, n2: n1[""atoms""] == n2[""atoms""],453                )454 455            # We should not end up here, but if we do, rotation should not be dummy456            return True457 458    # if not, compare immediate neighbors of i459    for n in nb:460        G.remove_edge(i, n)461 462    # make a set of each fragment around the chopped n-i bonds,463    # but only for fragments that are not root nor contain other random,464    # disconnected parts of the graph465    subgraphs_nodes = [466        _set for _set in connected_components(G) if root not in _set and any(n in _set for n in nb)467    ]468 469    if len(subgraphs_nodes) == 1:470        return True471        # if not, the torsion is likely to be rotable472        # (tetramethylguanidyl alanine C(β)-N bond)473 474    subgraphs = [subgraph(G, s) for s in subgraphs_nodes]475    for sub in subgraphs[1:]:476        if not is_isomorphic(477            subgraphs[0], sub, node_match=lambda n1, n2: n1[""atoms""] == n2[""atoms""]478        ):479            return True480    # Care should be taken because chiral centers are not taken into account: a rotation481    # involving an index where substituents only differ by stereochemistry, and where a482    # rotation is not an element of symmetry of the subsystem, the rotation is considered483    # dummy even if it would be more correct not to. For rotaionally corrected RMSD this484    # should only cause small inefficiencies and not lead to discarding any good conformer.485 486    return False487 488 489def get_hydrogen_bonds(490    coords: Array2D_float,491    atoms: Array1D_str,492    graph: Graph,493    d_min: float = 2.5,494    d_max: float = 3.3,495    max_angle: int = 45,496    elements: Sequence[Sequence[str]] | None = None,497    fragments: Sequence[Sequence[int]] | None = None,498) -> list[list[int]]:499    """"""Return a list of tuples with the indices of hydrogen bonding partners.500 501    An HB is a pair of atoms:502    - with one H and one X (N or O) atom503    - with an Y-X distance between d_min and d_max (i.e. N-O, Angstroms)504    - with an Y-H-X angle below max_angle (i.e. N-H-O, degrees)505 506    elements: iterable of two iterables with donor atomic symbols in the first507    element and acceptors in the second. default: ((""N"", ""O""), (""N"", ""O""))508 509    If fragments is specified (iterable of iterable of indices for each fragment)510    the function only returns inter-fragment hydrogen bonds.511    """"""512    hbs = []513    # initializing output list514 515    if elements is None:516        elements = ((""N"", ""O""), (""N"", ""O"", ""F""))517 518    het_idx_from = np.array([i for i, a in enumerate(atoms) if a in elements[0]], dtype=int)519    het_idx_to = np.array([i for i, a in enumerate(atoms) if a in elements[1]], dtype=int)520    # indices where N or O (or user-specified elements) atoms are present.521 522    for i1 in het_idx_from:523        for i2 in het_idx_to:524            # if inter-fragment HBs are requested, skip intra-HBs525            if fragments is not None:526                if any(((i1 in f and i2 in f) for f in fragments)):527                    continue528 529            # keep close pairs530            if d_min < np.linalg.norm(coords[i1] - coords[i2]) < d_max:531                # getting the indices of all H atoms attached to them532                Hs = [i for i in graph.neighbors(i1) if graph.nodes[i][""atoms""] == ""H""]533 534                # versor connectring the two Heteroatoms535                versor = coords[i2] - coords[i1]536                versor = versor / np.linalg.norm(versor)537 538                for iH in Hs:539                    # vectors connecting heteroatoms to H540                    v1 = coords[iH] - coords[i1]541                    v2 = coords[iH] - coords[i2]542 543                    # lengths of these vectors544                    d1 = np.linalg.norm(v1)545                    d2 = np.linalg.norm(v2)546 547                    # scalar projection in the heteroatom direction548                    l1 = v1 @ versor549                    l2 = v2 @ -versor550 551                    # largest planar angle between Het-H and Het-Het, in degrees (0 to 90°)552                    alfa = vec_angle(v1, versor) if l1 < l2 else vec_angle(v2, -versor)553 554                    # if the three atoms are not too far from being in line555                    if alfa < max_angle:556                        # adding the correct pair of atoms to results557                        if d1 < d2:558                            hbs.append(sorted((iH, i2)))559                        else:560                            hbs.append(sorted((iH, i1)))561 562                        break563 564    return hbs565 566 567def _get_rotation_mask(graph: Graph, torsion: Iterable[int]) -> Array1D_bool:568    """"""Return the rotation mask to be applied to coordinates before rotation.569 570    Get mask for the atoms that will rotate in a torsion:571    all the ones in the graph reachable from the last index572    of the torsion but not going through the central two573    atoms in the torsion quadruplet.574    """"""575    _, i2, i3, i4 = torsion576 577    graph.remove_edge(i2, i3)578    reachable_indices = shortest_path(graph, i4).keys()579    # get all indices reachable from i4 not going through i2-i3580 581    graph.add_edge(i2, i3)582    # restore modified graph583 584    mask = np.array([i in reachable_indices for i in graph.nodes], dtype=bool)585    # generate boolean mask586 587    # if np.count_nonzero(mask) > int(len(mask)/2):588    #     mask = ~mask589    # if we want to rotate more than half of the indices,590    # invert the selection so that we do less math591 592    mask[i3] = False593    # do not rotate i3: it would not move,594    # since it lies on the rotation axis595 596    return mask597 598 599def _get_quadruplets(graph: Graph) -> Array2D_int:600    """"""Return list of quadruplets that indicate potential torsions.""""""601    # Step 1: Find spanning tree602    spanning_tree = minimum_spanning_tree(graph)603 604    # Step 2: Add dihedrals for spanning tree605    dihedrals = []606 607    # For each edge in the spanning tree, we can potentially define a dihedral608    # We need edges that have at least 2 neighbors each to form a 4-point dihedral609    for edge in spanning_tree.edges():610        i, j = edge611 612        # Find neighbors of i and j in the original graph613        i_neighbors = [n for n in graph.neighbors(i) if n not in (i, j)]614        j_neighbors = [n for n in graph.neighbors(j) if n not in (i, j)]615 616        if len(i_neighbors) > 0 and len(j_neighbors) > 0:617            # Form dihedral: neighbor_of_i - i - j - neighbor_of_j618            k = i_neighbors[0]  # Choose first available neighbor619            m = j_neighbors[0]  # Choose first available neighbor620            dihedrals.append((k, i, j, m))621 622    return np.array(dihedrals)623 624 625def get_torsions(626    coords: Array2D_float,627    graph: Graph,628    hydrogen_bonds: list[list[int]],629    double_bonds: list[tuple[int, int]],630    keepdummy: bool = False,631    mode: str = ""csearch"",632) -> list[Torsion]:633    """"""Return list of Torsion objects.""""""634    torsions = []635    for path in _get_quadruplets(graph):636        _, i2, i3, _ = path637        bt = tuple(sorted((i2, i3)))638 639        if bt not in double_bonds:640            t = Torsion(*path)641            t.mode = mode642 643            # not including linear torsions (i.e. where three or more of the four atoms644            # are in line) will ignore all rotations involving alkynes and adjacent positions.645            # This will miss some potentially dummy rotations (i.e X-C#C-tBu) at the cost of646            # avoiding some more complex and potentially brittle way to account for these.647            # In any case, MOI-based pruning should account for such dummy rotations.648            if not is_linear(t, coords):649                # avoid torsions that are part of a cycle650                if not in_cycle(t, graph):651                    if is_rotable(t, graph, hydrogen_bonds, keepdummy=keepdummy):652                        torsions.append(t)653    # Create non-redundant torsion objects654    # Rejects (4,3,2,1) if (1,2,3,4) is present655    # Rejects torsions that do not represent a rotable bond656 657    return torsions658 659 660def rotationally_corrected_rmsd_and_max(661    ref: Array2D_float,662    coord: Array2D_float,663    atoms: Array1D_str,664    torsions: Array2D_int,665    graph: Graph,666    angles: Sequence[Sequence[int]],667    heavy_atoms_only: bool = True,668    debugfunction: Callable[..., object] | None = None,669    return_type: str = ""rmsd"",670) -> tuple[float, float] | Array2D_float:671    """"""Return RMSD and max deviation, corrected for degenerate torsions.672 673    Return a tuple with the RMSD between p and q674    and the maximum deviation of their positions.675    """"""676    assert return_type in (""rmsd"", ""coords"")677 678    torsion_corrections = [0 for _ in torsions]679 680    mask = (681        np.array([a != ""H"" for a in atoms]) if heavy_atoms_only else np.ones(len(atoms), dtype=bool)682    )683 684    # Now rotate every dummy torsion by the appropriate increment until we minimize local RMSD685    for i, torsion in enumerate(torsions):686        best_rmsd = 1e10687 688        # Look for the rotational angle set that minimizes the torsion RMSD and save it for later689        for angle in angles[i]:690            coord = rotate_dihedral(coord, torsion, angle, indices_to_be_moved=[torsion[3]])691 692            locally_corrected_rmsd, _ = rmsd_and_max(ref[torsion], coord[torsion])693 694            if locally_corrected_rmsd < best_rmsd:695                best_rmsd = locally_corrected_rmsd696                torsion_corrections[i] = angle697 698            # it is faster to undo the rotation rather than working with a copy of coords699            coord = rotate_dihedral(coord, torsion, -angle, indices_to_be_moved=[torsion[3]])700 701        # now rotate that angle to the desired orientation before going to the next angle702        if torsion_corrections[i] != 0:703            coord = rotate_dihedral(704                coord, torsion, torsion_corrections[i], mask=_get_rotation_mask(graph, torsion)705            )706 707        if debugfunction is not None:708            global_rmsd = rmsd_and_max(ref[mask], coord[mask])[0]709            debugfunction(710                f""    Torsion {i + 1} - {torsion}: best θ = {torsion_corrections[i]}°, ""711                + f""4-atom RMSD: {best_rmsd:.3f} Å, global RMSD: {global_rmsd:.3f} Å""712            )713 714    # we should have the optimal orientation on all torsions now:715    # calculate the RMSD716    rmsd, maxdev = rmsd_and_max(ref[mask], coord[mask])717 718    # since we could have segmented graphs, and therefore potentially only rotate719    # subsets of the graph where the torsion last two indices are,720    # we have to undo the final rotation too (would not be needed for connected graphs)721    for torsion, optimal_angle in zip(722        reversed(torsions), reversed(torsion_corrections), strict=False723    ):724        coord = rotate_dihedral(725            coord, torsion, -optimal_angle, mask=_get_rotation_mask(graph, torsion)726        )727 728    if return_type == ""rmsd"":729        return rmsd, maxdev730 731    return coord732","Python"
733"Conformation","ntampellini/prism_pruner","prism_pruner/utils.py",".py","4594","154","""""""PRISM - Pruning Interface for Similar Molecules.""""""734 735from typing import Any, Sequence736 737import numpy as np738from numpy.linalg import LinAlgError739from numpy.typing import ArrayLike740 741from prism_pruner.algebra import get_alignment_matrix, rot_mat_from_pointer742from prism_pruner.typing import Array1D_bool, Array1D_int, Array1D_str, Array2D_float, Array3D_float743 744EH_TO_EV = 27.211399745EH_TO_KCAL = 627.5096080305927746EV_TO_KCAL = 23.060541945329334747 748 749def align_structures(750    structures: Array3D_float, indices: Array1D_int | None = None751) -> Array3D_float:752    """"""Align structures.753 754    Aligns molecules of a structure array (shape is (n_structures, n_atoms, 3))755    to the first one, based on the indices. If not provided, all atoms are used756    to get the best alignment. Return is the aligned array.757    """"""758    reference = structures[0]759    targets = structures[1:]760    if isinstance(indices, (list, tuple)):761        indices = np.array(indices)762 763    indices = indices if indices is not None else np.array([i for i, _ in enumerate(structures[0])])764 765    reference -= np.mean(reference[indices], axis=0)766    for t, _ in enumerate(targets):767        targets[t] -= np.mean(targets[t, indices], axis=0)768 769    output = np.zeros(structures.shape)770    output[0] = reference771 772    for t, target in enumerate(targets):773        try:774            matrix = get_alignment_matrix(reference[indices], target[indices])775 776        except LinAlgError:777            # it is actually possible for the kabsch alg not to converge778            matrix = np.eye(3)779 780        # output[t+1] = np.array([matrix @ vector for vector in target])781        output[t + 1] = (matrix @ target.T).T782 783    return output784 785 786def time_to_string(total_time: float, verbose: bool = False, digits: int = 1) -> str:787    """"""Convert totaltime (float) to a timestring with hours, minutes and seconds.""""""788    timestring = """"789 790    names = (""days"", ""hours"", ""minutes"", ""seconds"") if verbose else (""d"", ""h"", ""m"", ""s"")791 792    if total_time > 24 * 3600:793        d = total_time // (24 * 3600)794        timestring += f""{int(d)} {names[0]} ""795        total_time %= 24 * 3600796 797    if total_time > 3600:798        h = total_time // 3600799        timestring += f""{int(h)} {names[1]} ""800        total_time %= 3600801 802    if total_time > 60:803        m = total_time // 60804        timestring += f""{int(m)} {names[2]} ""805        total_time %= 60806 807    timestring += f""{round(total_time, digits):{2 + digits}} {names[3]}""808 809    return timestring810 811 812double_bonds_thresholds_dict = {813    ""CC"": 1.4,814    ""CN"": 1.3,815}816 817 818def get_double_bonds_indices(coords: Array2D_float, atoms: Array1D_str) -> list[tuple[int, int]]:819    """"""Return a list containing 2-elements tuples of indices involved in any double bond.""""""820    mask = atoms != ""H""821    numbering = np.arange(len(coords))[mask]822    coords = coords[mask]823    atoms_masked = atoms[mask]824    output = []825 826    for i1, _ in enumerate(coords):827        for i2 in range(i1 + 1, len(coords)):828            dist = np.linalg.norm(coords[i1] - coords[i2])829            tag = """".join(sorted([atoms_masked[i1], atoms_masked[i2]]))830 831            threshold = double_bonds_thresholds_dict.get(tag)832            if threshold is not None and dist < threshold:833                output.append((numbering[i1], numbering[i2]))834 835    return output836 837 838def rotate_dihedral(839    coords: Array2D_float,840    dihedral: list[int] | tuple[int, ...],841    angle: float,842    mask: Array1D_bool | None = None,843    indices_to_be_moved: ArrayLike | None = None,844) -> Array2D_float:845    """"""Rotate a molecule around a given bond.846 847    Atoms that will move are the ones848    specified by mask or indices_to_be_moved.849    If both are None, only the first index of850    the dihedral iterable is moved.851 852    angle: angle, in degrees853    """"""854    i1, i2, i3, *_ = dihedral855 856    if indices_to_be_moved is not None:857        mask = np.isin(np.arange(len(coords)), indices_to_be_moved)858 859    if mask is None:860        mask = np.zeros(len(coords), dtype=bool)861        mask[i1] = True862 863    axis = coords[i2] - coords[i3]864    mat = rot_mat_from_pointer(axis, angle)865 866    center = coords[i3]867    coords[mask] = (coords[mask] - center) @ mat.T + center868 869    return coords870 871 872def flatten(array: Sequence[Any], typefunc: type = float) -> list[Any]:873    """"""Return the unraveled sequence, with items coerced into the typefunc type.""""""874    out = []875 876    def rec(_l: Any) -> None:877        """"""Recursive unraveling function.""""""878        for e in _l:879            if type(e) in [list, tuple, np.ndarray]:880                rec(e)881            else:882                out.append(typefunc(e))883 884    rec(array)885    return out886","Python"
887"Conformation","ntampellini/prism_pruner","prism_pruner/graph_manipulations.py",".py","6321","196","""""""Graph manipulation utilities for molecular structures.""""""888 889from functools import lru_cache890 891import numpy as np892from networkx import Graph, all_simple_paths, from_numpy_array, set_node_attributes893from periodictable import elements894from scipy.spatial.distance import cdist895 896from prism_pruner.algebra import dihedral897from prism_pruner.typing import Array1D_bool, Array1D_str, Array2D_float898 899 900@lru_cache()901def d_min_bond(a1: str, a2: str, factor: float = 1.2) -> float:902    """"""Return the bond distance between two atoms.""""""903    return factor * (elements.symbol(a1).covalent_radius + elements.symbol(a2).covalent_radius)  # type: ignore [no-any-return]904 905 906def graphize(907    atoms: Array1D_str,908    coords: Array2D_float,909    mask: Array1D_bool | None = None,910) -> Graph:911    """"""912    Return a NetworkX undirected graph of molecular connectivity.913 914    :param atoms: atomic symbols915    :param coords: atomic coordinates as 3D vectors916    :param mask: bool array, with False for atoms to be excluded in the bond evaluation917    :return: connectivity graph918    """"""919    mask = np.array([True for _ in atoms], dtype=bool) if mask is None else mask920    assert len(coords) == len(atoms)921    assert len(coords) == len(mask)922 923    matrix = np.zeros((len(coords), len(coords)))924    for i, mask_i in enumerate(mask):925        if not mask_i:926            continue927 928        for j, mask_j in enumerate(mask[i + 1 :], start=i + 1):929            if not mask_j:930                continue931 932            if np.linalg.norm(coords[i] - coords[j]) < d_min_bond(atoms[i], atoms[j]):933                matrix[i][j] = 1934 935    graph = from_numpy_array(matrix)936    set_node_attributes(graph, dict(enumerate(atoms)), ""atoms"")937 938    return graph939 940 941def get_sp_n(index: int, graph: Graph) -> int | None:942    """"""943    Get hybridization of selected atom.944 945    Return n, that is the apex of sp^n hybridization for CONPS atoms.946    This is just an assimilation to the carbon geometry in relation to sp^n:947    - sp¹ is linear948    - sp² is planar949    - sp³ is tetrahedral950    This is mainly used to understand if a torsion is to be rotated or not.951    """"""952    atom = graph.nodes[index][""atoms""]953 954    if atom not in {""C"", ""N"", ""O"", ""P"", ""S""}:955        return None956 957    # Relationship of number of neighbors to sp^n hybridization958    d: dict[str, dict[int, int | None]] = {959        ""C"": {2: 1, 3: 2, 4: 3},960        ""N"": {2: 2, 3: None, 4: 3},  # 3 could mean sp3 or sp2961        ""O"": {1: 2, 2: 3, 3: 3, 4: 3},962        ""P"": {2: 2, 3: 3, 4: 3},963        ""S"": {2: 2, 3: 3, 4: 3},964    }965    return d[atom].get(len(set(graph.neighbors(index))))966 967 968def is_amide_n(index: int, graph: Graph, mode: int = -1) -> bool:969    """"""970    Assess if the atom is an amide-like nitrogen.971 972    Note: carbamates and ureas are considered amides.973 974    mode:975    -1 - any amide976    0 - primary amide (CONH2)977    1 - secondary amide (CONHR)978    2 - tertiary amide (CONR2)979    """"""980    # Must be a nitrogen atom981    if graph.nodes[index][""atoms""] == ""N"":982        nb = set(graph.neighbors(index))983        nb_atoms = [graph.nodes[j][""atoms""] for j in nb]984 985        if mode != -1:986            # Primary amides need to have 1H, secondary amides none987            if nb_atoms.count(""H"") != (2, 1, 0)[mode]:988                return False989 990        for n in nb:991            # There must be at least one carbon atom next to N992            if graph.nodes[n][""atoms""] == ""C"":993                nb_nb = set(graph.neighbors(n))994                # Bonded to three atoms995                if len(nb_nb) == 3:996                    # and at least one of them has to be an oxygen997                    if ""O"" in {graph.nodes[i][""atoms""] for i in nb_nb}:998                        return True999    return False1000 1001 1002def is_ester_o(index: int, graph: Graph) -> bool:1003    """"""1004    Assess if the index is an ester-like oxygen.1005 1006    Note: carbamates and carbonates return True, carboxylic acids return False.1007    """"""1008    if graph.nodes[index][""atoms""] == ""O"":1009        if ""H"" in (nb := set(graph.neighbors(index))):1010            return False1011 1012        for n in nb:1013            if graph.nodes[n][""atoms""] == ""C"":1014                nb_nb = set(graph.neighbors(n))1015                if len(nb_nb) == 3:1016                    nb_nb_sym = [graph.nodes[i][""atoms""] for i in nb_nb]1017                    if nb_nb_sym.count(""O"") > 1:1018                        return True1019    return False1020 1021 1022def is_phenyl(coords: Array2D_float) -> bool:1023    """"""1024    Assess if the six atomic coords refer to a phenyl-like ring.1025 1026    Note: quinones evaluate to True1027 1028    :param coords: six coordinates of C/N atoms1029    :return: bool indicating if the six atoms look like part of a phenyl/naphtyl/pyridine1030             system, coordinates for the center of that ring1031    """"""1032    # if any atomic couple is more than 3 A away from each other, this is not a Ph1033    if np.max(cdist(coords, coords)) > 3:1034        return False1035 1036    threshold_delta: float = 1 - np.cos(10 * np.pi / 180)1037    flat_delta: float = 1 - np.abs(np.cos(dihedral(coords[[0, 1, 2, 3]]) * np.pi / 180))1038 1039    return flat_delta < threshold_delta1040 1041 1042def get_phenyl_ids(index: int, graph: Graph) -> list[int] | None:1043    """"""If index is part of a phenyl, return the six heavy atoms ids associated with the ring.""""""1044    for n in graph.neighbors(index):1045        for path in all_simple_paths(graph, source=index, target=n, cutoff=6):1046            if len(path) != 6 or any(graph.nodes[n][""atoms""] == ""H"" for n in path):1047                continue1048            if all(len(set(graph.neighbors(i))) == 3 for i in path):1049                return path  # type: ignore [no-any-return]1050 1051    return None1052 1053 1054def find_paths(1055    graph: Graph,1056    u: int,1057    n: int,1058    exclude_set: set[int] | None = None,1059) -> list[list[int]]:1060    """"""1061    Find paths in graph.1062 1063    Recursively find all paths of a NetworkX graph with length = n, starting from node u.1064 1065    :param graph: NetworkX graph1066    :param u: starting node1067    :param n: path length1068    :param exclude_set: set of nodes to exclude from the paths1069    :return: list of paths (each path is a list of node indices)1070    """"""1071    exclude_set = (exclude_set or set()) | {u}1072 1073    if n == 0:1074        return [[u]]1075 1076    return [1077        [u, *path]1078        for neighbor in graph.neighbors(u)1079        if neighbor not in exclude_set1080        for path in find_paths(graph, neighbor, n - 1, exclude_set)1081    ]1082","Python"
1083"Conformation","ntampellini/prism_pruner","prism_pruner/pruner.py",".py","24132","740","""""""PRISM - Pruning Interface for Similar Molecules.""""""1084 1085from dataclasses import dataclass, field1086from time import perf_counter1087from typing import Any, Callable, Sequence1088 1089import numpy as np1090from networkx import Graph, connected_components1091from periodictable import elements1092from scipy.spatial.distance import cdist1093 1094from prism_pruner.algebra import get_inertia_moments1095from prism_pruner.graph_manipulations import graphize1096from prism_pruner.rmsd import rmsd_and_max1097from prism_pruner.torsion_module import (1098    get_angles,1099    get_hydrogen_bonds,1100    get_torsions,1101    is_nondummy,1102    rotationally_corrected_rmsd_and_max,1103)1104from prism_pruner.typing import (1105    Array1D_bool,1106    Array1D_float,1107    Array1D_int,1108    Array1D_str,1109    Array2D_float,1110    Array2D_int,1111    Array3D_float,1112)1113from prism_pruner.utils import flatten, get_double_bonds_indices, time_to_string1114 1115 1116@dataclass1117class PrunerConfig:1118    """"""Configuration dataclass for Pruner.""""""1119 1120    structures: Array3D_float1121 1122    # Optional parameters that get initialized1123    energies: Array1D_float = field(default_factory=lambda: np.array([]))1124    max_dE: float = field(default=0.0)1125    debugfunction: Callable[[str], None] | None = field(default=None)1126 1127    # Computed fields1128    eval_calls: int = field(default=0, init=False)1129    cache_calls: int = field(default=0, init=False)1130    cache: set[tuple[int, int]] = field(default_factory=lambda: set(), init=False)1131 1132    def __post_init__(self) -> None:1133        """"""Validate inputs and initialize computed fields.""""""1134        # validate input types1135        assert type(self.structures) is np.ndarray1136 1137        self.mask = np.ones(shape=(self.structures.shape[0],), dtype=np.bool_)1138 1139        if len(self.energies) != 0:1140            assert self.max_dE > 0.0, (1141                ""If you provide energies, please also provide an appropriate energy window max_dE.""1142            )1143 1144        # Set defaults for optional parameters1145        if len(self.energies) == 0:1146            assert type(self.energies) is np.ndarray1147            self.energies = np.zeros(self.structures.shape[0], dtype=float)1148 1149        assert len(self.energies) == len(self.structures), (1150            ""Please make sure that the energies ""1151            + ""provided have the same len as the input structures.""1152        )1153 1154        if self.max_dE == 0.0:1155            self.max_dE = 1.01156 1157    def evaluate_sim(self, *args: Any, **kwargs: Any) -> bool:1158        """"""Stub method - override in subclasses as needed.""""""1159        raise NotImplementedError1160 1161 1162@dataclass1163class RMSDRotCorrPrunerConfig(PrunerConfig):1164    """"""Configuration dataclass for Pruner.""""""1165 1166    atoms: Array1D_str = field(kw_only=True)1167    max_rmsd: float = field(kw_only=True)1168    max_dev: float = field(kw_only=True)1169    angles: Sequence[Sequence[int]] = field(kw_only=True)1170    torsions: Array2D_int = field(kw_only=True)1171    graph: Graph = field(kw_only=True)1172    heavy_atoms_only: bool = True1173 1174    def __post_init__(self) -> None:1175        """"""Add type enforcing to the parent's __post_init__.""""""1176        super().__post_init__()1177 1178        # validate input types1179        assert type(self.atoms) is np.ndarray1180        assert type(self.graph) is Graph1181 1182    def evaluate_sim(self, i1: int, i2: int) -> bool:1183        """"""Return whether the structures are similar.""""""1184        rmsd, max_dev = rotationally_corrected_rmsd_and_max(1185            self.structures[i1],1186            self.structures[i2],1187            atoms=self.atoms,1188            torsions=self.torsions,1189            graph=self.graph,1190            angles=self.angles,1191            # debugfunction=self.debugfunction, # lots of printout1192            heavy_atoms_only=self.heavy_atoms_only,1193        )1194 1195        if rmsd > self.max_rmsd:1196            return False1197 1198        if max_dev > self.max_dev:1199            return False1200 

Showing the first 1,200 of 54741 lines. Download the file for the rest.