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._cleanStProperties()
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()
def _cleanStProperties(self):
'''
Clean Fast3D (f3d) properties whose presence can cause f3d call to work
improperly in this class. For the purposes we're using f3d here, no
specialized f3d properties are needed. (see DESMOND-11503).
'''
for prop in self._st.getPropertyNames():
if structure.PropertyName(prop).family == 'f3d':
del self._st.property[prop]
[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()