import math
from past.utils import old_div
from typing import List
import numpy as np
from schrodinger import structure
from schrodinger.application.desmond.constants import ORIGINAL_INDEX
from schrodinger.infra import fast3d
from schrodinger.infra import licensing
from schrodinger.infra import mm
from schrodinger.structutils import analyze
from schrodinger.structutils import measure
from schrodinger.structutils import minimize
from schrodinger.utils.env import EnvironmentContext
# ffld might not be installed
try:
from schrodinger.forcefield.packages import altcap
except ImportError:
altcap = None
CONF_ANGLE_PROP = 'r_m_conformer_angle'
MIN_ENERGY_PROP = 'r_m_min_energy'
MIN_OPT = {
'max_steps': 150,
'dielectric_constant': 1.0,
'energy_criterion': 5e-5,
'gradient_criterion': 0.1,
'cleanup': False,
'honor_pbc': False
}
[docs]class TorsionPotential:
"""A class to contain the atom numbers for each torsion"""
# Class variable that contains torsion parameters relevant to
# ligand torsions. This variable is created once when the first
# instance of Torsion class is initialized.
dih = None
[docs] def __init__(self,
a1,
a2,
a3,
a4,
cms,
lig_atoms,
results=None,
calc_tors=True):
(self.a1, self.a2, self.a3, self.a4) = (a1, a2, a3, a4)
self.cms = cms
self.lfrom = None
self.lto = None
self.col = None
self.rb_potential = None
self.nbins = 36
self.interval_list = list(
np.arange(-180.0, 180.0, old_div(360.0, self.nbins)))
# add last bead
self.interval_list.append(self.interval_list[-1] +
old_div(360.0, self.nbins))
self.results = results
self.lig_atoms = lig_atoms
# new_lig_atoms contains aid of atoms according to a CT the
# ligand is in
self.new_lig_atoms = None
tmp_list = cms.select_atom_comp(self._convert_index_to_asl(lig_atoms))
for i, alist in enumerate(tmp_list):
if alist:
self.new_lig_atoms = tmp_list[i]
self.atom_dict = {}
if len(self.lig_atoms) == len(self.new_lig_atoms):
for i, j in zip(self.lig_atoms, self.new_lig_atoms):
self.atom_dict[i] = j
if calc_tors:
self.potential = self._get_utt_potential(
[self.a1, self.a2, self.a3, self.a4])
else:
self.potential = (self.nbins + 1) * [0.0]
[docs] def set_color(self, col):
self.col = col
[docs] def get_color(self):
return self.col
[docs] def set_ligand_from(self, lfrom):
self.lfrom = lfrom
[docs] def set_ligand_to(self, lto):
self.lto = lto
[docs] def get_angles(self):
return self.interval_list
[docs] def set_rbpotential(self, rbp):
self.rb_potential = rbp
[docs] def get_potential(self):
if self.rb_potential:
return self.rb_potential
else:
return self.potential
def __str__(self):
ret = 'Torsion:'
ret += ' ' + str(self.a1) + '-' + str(self.a2) + '-' + str(
self.a3) + '-' + str(self.a4)
ret += '\n from:' + str(self.lfrom) + '; to:' + str(self.lto)
return ret
# functions needed to show torsion potentials on histograms
def _convert_index_to_asl(self, index_list):
asl = '(atom.num '
for i in index_list:
asl += '%s,' % i
asl = asl[0:len(asl) - 1] + ')'
return asl
def _get_utt_potential(self, atom_list):
'''
Given atom_list (i,j,k,l), to get a full potential around j-k in st,
we need to also look at all the atoms connected to j and k, and
reconstruct UTT
'''
st = self.cms
jatoms = []
katoms = []
# get atoms connected to j
for bond in st.atom[atom_list[1]].bond:
if bond.atom2.index != atom_list[2] and bond.order:
jatoms.append(bond.atom2.index)
# get atoms conntecte to k
for bond in st.atom[atom_list[2]].bond:
if bond.atom2.index != atom_list[1] and bond.order:
katoms.append(bond.atom2.index)
utt_tor = []
for j in jatoms:
for k in katoms:
utt_tor.append([j, atom_list[1], atom_list[2], k])
total_utt_potential = np.zeros(self.nbins + 1)
for utt in utt_tor:
try:
total_utt_potential += self._get_potential(st, utt)
except:
continue
return total_utt_potential
def _get_energy(self, tor, angle):
pot = 0.0
rad = angle / 180. * math.pi
pot += tor.c1 / 2. * (1 + math.cos(1 * rad))
pot += tor.c2 / 2. * (1 - math.cos(2 * rad))
pot += tor.c3 / 2. * (1 + math.cos(3 * rad))
pot += tor.c4 / 2. * (1 - math.cos(4 * rad))
return pot
def _get_potential(self, st, atom_list):
if not TorsionPotential.dih:
self.init_ff(st, atom_list)
dihed_p = None # dihedral parmeters
for e in TorsionPotential.dih:
ai = self.atom_dict[atom_list[0]]
aj = self.atom_dict[atom_list[1]]
ak = self.atom_dict[atom_list[2]]
al = self.atom_dict[atom_list[3]]
# match forward and backwards
if (e.ai == ai and e.aj == aj and e.ak == ak and e.al == al) or \
(e.ai == al and e.aj == ak and e.ak == aj and e.al == ai):
dihed_p = e
break
# check if dihedral parameters are either protected (INVALID) or
# all zeros
if not dihed_p or dihed_p.funct == 'INVALID' or \
len(set([dihed_p.c1, dihed_p.c2, dihed_p.c3, dihed_p.c4])) == 1:
return np.zeros(len(self.interval_list))
tor_potential = []
for a in self.interval_list:
tor_potential.append(self._get_energy(dihed_p, a))
return np.array(tor_potential)
[docs] def init_ff(self, st, atom_list):
"""
This function initializes torsion force field parameters
relevant to the ligand. Ligand is defined by a list of
atoms.
:param st: structure object for all molecules
:type st: `structure.Structure`
:param atom_list: list of ligand atom indices
:type atom_list: list
"""
selection = st.select_atom_comp(self._convert_index_to_asl(atom_list))
icomp = None
for i, element in enumerate(selection):
if element:
icomp = i
break
ct = st.comp_ct[icomp]
dih = ct.ffio.dihedral
TorsionPotential.dih = []
for e in dih:
if e.ai in self.new_lig_atoms and e.aj in self.new_lig_atoms and \
e.ak in self.new_lig_atoms and e.al in self.new_lig_atoms:
TorsionPotential.dih.append(e)
[docs]class BondRotator(object):
"""
perform RB scan
"""
[docs] def __init__(self,
st,
rb_atoms,
angle_diff=10,
ff=16,
debug=False,
canonicalize=True,
sampling=True):
self._fk = 30.0
self._st = st
self._rb_atoms = rb_atoms
self.dihedral_atoms = [st.atom[i] for i in rb_atoms]
self._angle_diff = angle_diff
self._ff = ff
self._sampling = sampling
self._min_opt = MIN_OPT.copy()
self._min_opt['ffld_version'] = self._ff
self._min_opt['struct'] = self._st
# Measure the current dihedral angle:
angle = measure.measure_dihedral_angle(*self.dihedral_atoms)
self.input_angle = angle
self.num_rotamers = 360 // angle_diff
self.angles = list(range(-180, 180 + int(angle_diff), angle_diff))
self.writer = None
if debug:
self.writer = structure.StructureWriter('rotator_%i_%i.mae' %
(rb_atoms[1], rb_atoms[2]))
with minimize.minimizer_context(**self._min_opt) as self.minimizer:
if sampling:
# Run with fast3d scan
self.search_and_scan(self._st)
else:
self._scan(self._st, canonicalize=canonicalize)
if self.writer:
self.writer.close()
[docs] def potential(self, conf=None, degree=True):
'''
get potential energy corresponding to 'conf'
:param conf: conformation in degrees to get corresponding potential
:type conf: float
:rtype float
'''
if not conf:
conf = self.input_angle
if not degree:
conf = conf * 180 / 3.14159
energy = self.getRotEnergy(offset=True)
for i, start in enumerate(self.angles[0:-1]):
end = self.angles[i + 1]
if start < conf and end > conf:
i_start = i
break
# Interpolate potential at 'conf' using simple line (y=mx+b) equation
e1 = energy[i_start]
e2 = energy[i_start + 1]
ang_diff = self.angles[i_start + 1] - self.angles[i_start]
x = (conf - self.angles[i_start]) / ang_diff
b = e1
m = e2 - e1
return (m * x) + b
def _generate_conformations(self,
max_confs: int = 10
) -> List[structure.Structure]:
"""
Run ConfGenX (fast3d) to generate `max_conf` desired cts.
"""
mode = 'default'
optimize = True
use_default_fragment_libraries = True
custom_fragment_libraries = []
# Set to use the FEP_GPGPU license
with EnvironmentContext('SCHRODINGER_CONFGENX_FEPGPGPU_MODE', "1"):
engine = fast3d.Engine(mode, optimize,
use_default_fragment_libraries,
custom_fragment_libraries)
engine.setDesiredNumberOfConformers(max_confs)
confs = engine.run(self._st)
ret_cts = []
for (energy, xyz, _) in confs:
ct = self._st.copy()
xyz = np.asarray(xyz).reshape(ct.atom_total, 3)
ct.setXYZ(xyz)
ret_cts.append(ct)
return ret_cts
[docs] def search_and_scan(self, ct):
"""
Get best profile by using fast3d to generate multiple initial
coordinates. `self.results` to contain the best profile.
"""
profiles = []
self._scan(ct)
profiles.append(self.results)
for ist, st in enumerate(self._generate_conformations()):
st.title = f'conf_{ist}'
self._scan(st, canonicalize=False)
profiles.append(self.results)
best_profile = []
for i in range(len(profiles[0])): # angles
es = [p[i] for p in profiles]
best_profile.append(min(es))
self.results = best_profile
def _scan(self, ct, canonicalize=True):
self.results = [None] * len(self.angles)
st = ct.copy()
if canonicalize:
# Set to use the FEP_GPGPU license
with EnvironmentContext('SCHRODINGER_CONFGENX_FEPGPGPU_MODE', "1"):
fast3d.SingleConformerEngine().run(st)
self.minimizer.updateCoordinates(st)
min_st = self.minimizer.getStructure()
angle2idx = {a: index for index, a in enumerate(self.angles)}
dihedral_atoms = [min_st.atom[i] for i in self._rb_atoms]
diff = np.abs(
np.array(self.angles) -
measure.measure_dihedral_angle(*dihedral_atoms))
idx = diff.argmin()
# Start with the angle closest to original angle of `min_st`.
angles = self.angles[idx:] + self.angles[:idx]
title = st.title
for rt in angles:
min_st.adjust(rt, *self._rb_atoms)
self.minimizer.addTorsionRestraint(*self._rb_atoms,
self._fk,
rt,
flat_bottom=1.0)
min_res = self.minimizer.minimize(
) # min_st's coordinates are updated.
self.results[angle2idx[rt]] = min_res.potential_energy
if self.writer:
min_st.title = f'{title}:{rt}deg'
min_st.property[CONF_ANGLE_PROP] = rt
min_st.property[MIN_ENERGY_PROP] = min_res.potential_energy
self.writer.append(min_st)
self.minimizer.deleteAllRestraints(rest_type=mm.MMFfldTorsionRest)
[docs] def getRotEnergy(self, offset=False):
if offset:
res = np.array(self.results)
e_min = res.min()
return np.round(res - e_min, decimals=3).tolist()
return self.results
[docs]def get_old2new(st, prop=ORIGINAL_INDEX):
"""
Get dictionary for atom-level property 'prop' as a key, and return
it's new atom index (aid)
"""
atom_map = {}
for a in st.atom:
if prop in a.property:
old_aid = a.property[prop]
atom_map[old_aid] = a.index
return atom_map
[docs]class IncompleteFragmentError(Exception):
pass
[docs]def get_rb_potential(cms_model, lig_aids, a1, a2, a3, a4):
# full_system atoms indices for the torsion
tor_sys_idx = [a1, a2, a3, a4]
lig_st = cms_model.extract(lig_aids)
return rb_potential_from_struct(lig_st, tor_sys_idx)
[docs]def rb_potential_from_struct(lig_st, tor_sys_idx):
"""
Creates bond rotator object to calculate torsional potential for a
given structure. The input structure should have 'ORIGINAL_INDEX'
property.
:param lig_st: ligand structure
:type lig_st: structure.Structure
:param tor_sys_idx: list of atom indices to define torsion
:type tor_sys_idx: list[int]
:return: bond rotator object
:rtype: BondRotator
"""
if not altcap or not licensing.licenseExists(licensing.FFLD_DEVELOPMENT):
return None
if not analyze.hydrogens_present(lig_st):
from schrodinger.structutils.build import add_hydrogens
lig_st = lig_st.copy()
add_hydrogens(lig_st)
ori2new_dict = get_old2new(lig_st)
# extracted ligand atoms indices for the torsion
tor_lig_idx = [ori2new_dict[tor_sys_idx[i]] for i in range(4)]
rb_key = (tor_lig_idx[1], tor_lig_idx[2])
# get fragment that has equivalent atom-types and preserves UTT type
with altcap.Capper() as capper:
frag_st = capper.get_fragment(lig_st, rot_bond_atoms=rb_key)
old2new_dict = get_old2new(frag_st, altcap.ORIG_ATOM_NUM_PROP)
# extracted fragment atoms indices for the torsion
tor_frag_idx = [old2new_dict[tor_lig_idx[i]] for i in range(4)]
energy_rotamer = BondRotator(frag_st, tor_frag_idx, sampling=True)
return energy_rotamer
[docs]def get_rb_torsions_potential(cms_model, lig_aids, a1, a2, a3, a4):
t = TorsionPotential(a1, a2, a3, a4, cms_model, lig_aids)
return np.round(t.potential, decimals=3).tolist()