a b/unimol/utils/geom.py
1
import random
2
import numpy as np
3
import rdkit
4
import rdkit.Chem as Chem
5
import rdkit.Chem.AllChem as AllChem
6
import rdkit.Chem.Descriptors as Descriptors
7
import tqdm
8
9
def change_torsion(pointset, index, angle=None):
10
    # pointset: (N, 3)
11
    #return the new pointset
12
    # angle: rad
13
    # index: (2)
14
    if angle is None:
15
        angle = random.uniform(-np.pi, np.pi)
16
    pointset = np.array(pointset)
17
    #get rotate axis
18
    axis = pointset[index[0][-1]] - pointset[index[1][0]]
19
    axis = axis / np.linalg.norm(axis)
20
    #get rotate matrix
21
    cos = np.cos(angle)
22
    sin = np.sin(angle)
23
    R = np.array([[cos + axis[0] ** 2 * (1 - cos), axis[0] * axis[1] * (1 - cos) - axis[2] * sin, axis[0] * axis[2] * (1 - cos) + axis[1] * sin],
24
                    [axis[1] * axis[0] * (1 - cos) + axis[2] * sin, cos + axis[1] ** 2 * (1 - cos), axis[1] * axis[2] * (1 - cos) - axis[0] * sin],
25
                    [axis[2] * axis[0] * (1 - cos) - axis[1] * sin, axis[2] * axis[1] * (1 - cos) + axis[0] * sin, cos + axis[2] ** 2 * (1 - cos)]])
26
    #rotate
27
    pointset[index[1]] = np.dot(pointset[index[1]] - pointset[index[1][0]], R) + pointset[index[1][0]]
28
    return pointset, angle
29
30
def gen_conformation(mol, num_conf=20, num_worker=8):
31
    try:
32
        mol = Chem.AddHs(mol)
33
        AllChem.EmbedMultipleConfs(mol, numConfs=num_conf, numThreads=num_worker, pruneRmsThresh=0.1, maxAttempts=5, useRandomCoords=False)
34
        AllChem.MMFFOptimizeMoleculeConfs(mol, numThreads=num_worker)
35
        mol = Chem.RemoveHs(mol)
36
        return mol
37
    except:
38
        return None
39
40
def RotatableBond(mol):
41
    """Get the rotatable bond index of a molecule.
42
    Args:
43
        mol: rdkit.Chem.rdchem.Mol
44
    Returns:
45
        rotatable_bond: list of tuple
46
    """
47
    rotatable_bond = []
48
    for bond in mol.GetBonds():
49
        if bond.IsInRing():
50
            continue
51
        if bond.GetBondType() == Chem.rdchem.BondType.SINGLE:
52
            rotatable_bond.append((bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()))
53
    return rotatable_bond