File size: 3,057 Bytes
ce9690f |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 |
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]))
# Optimize conformations
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
|