| import copy, time |
| import numpy as np |
| from collections import defaultdict |
| from rdkit import Chem, RDLogger |
| from rdkit.Chem import AllChem, rdMolTransforms |
| from rdkit import Geometry |
| import networkx as nx |
| from scipy.optimize import differential_evolution |
|
|
| RDLogger.DisableLog('rdApp.*') |
|
|
| """ |
| Conformer matching routines from Torsional Diffusion |
| """ |
|
|
| def GetDihedral(conf, atom_idx): |
| return rdMolTransforms.GetDihedralRad(conf, atom_idx[0], atom_idx[1], atom_idx[2], atom_idx[3]) |
|
|
|
|
| def SetDihedral(conf, atom_idx, new_vale): |
| rdMolTransforms.SetDihedralRad(conf, atom_idx[0], atom_idx[1], atom_idx[2], atom_idx[3], new_vale) |
|
|
|
|
| def apply_changes(mol, values, rotatable_bonds, conf_id): |
| opt_mol = copy.copy(mol) |
| [SetDihedral(opt_mol.GetConformer(conf_id), rotatable_bonds[r], values[r]) for r in range(len(rotatable_bonds))] |
| return opt_mol |
|
|
|
|
| def optimize_rotatable_bonds(mol, true_mol, rotatable_bonds, probe_id=-1, ref_id=-1, seed=0, popsize=15, maxiter=500, |
| mutation=(0.5, 1), recombination=0.8): |
| opt = OptimizeConformer(mol, true_mol, rotatable_bonds, seed=seed, probe_id=probe_id, ref_id=ref_id) |
| max_bound = [np.pi] * len(opt.rotatable_bonds) |
| min_bound = [-np.pi] * len(opt.rotatable_bonds) |
| bounds = (min_bound, max_bound) |
| bounds = list(zip(bounds[0], bounds[1])) |
|
|
| |
| result = differential_evolution(opt.score_conformation, bounds, |
| maxiter=maxiter, popsize=popsize, |
| mutation=mutation, recombination=recombination, disp=False, seed=seed) |
| opt_mol = apply_changes(opt.mol, result['x'], opt.rotatable_bonds, conf_id=probe_id) |
|
|
| return opt_mol |
|
|
|
|
| class OptimizeConformer: |
| def __init__(self, mol, true_mol, rotatable_bonds, probe_id=-1, ref_id=-1, seed=None): |
| super(OptimizeConformer, self).__init__() |
| if seed: |
| np.random.seed(seed) |
| self.rotatable_bonds = rotatable_bonds |
| self.mol = mol |
| self.true_mol = true_mol |
| self.probe_id = probe_id |
| self.ref_id = ref_id |
|
|
| def score_conformation(self, values): |
| for i, r in enumerate(self.rotatable_bonds): |
| SetDihedral(self.mol.GetConformer(self.probe_id), r, values[i]) |
| return AllChem.AlignMol(self.mol, self.true_mol, self.probe_id, self.ref_id) |
|
|
|
|
| def get_torsion_angles(mol): |
| torsions_list = [] |
| G = nx.Graph() |
| for i, atom in enumerate(mol.GetAtoms()): |
| G.add_node(i) |
| nodes = set(G.nodes()) |
| for bond in mol.GetBonds(): |
| start, end = bond.GetBeginAtomIdx(), bond.GetEndAtomIdx() |
| G.add_edge(start, end) |
| for e in G.edges(): |
| G2 = copy.deepcopy(G) |
| G2.remove_edge(*e) |
| if nx.is_connected(G2): continue |
| l = list(sorted(nx.connected_components(G2), key=len)[0]) |
| if len(l) < 2: continue |
| n0 = list(G2.neighbors(e[0])) |
| n1 = list(G2.neighbors(e[1])) |
| torsions_list.append( |
| (n0[0], e[0], e[1], n1[0]) |
| ) |
| return torsions_list |
|
|
|
|
