In [1]:
import AlGDock.BindingPMF_plots
from AlGDock.BindingPMF_plots import *
import os, shutil, glob
phases = ['NAMD_Gas', 'NAMD_OBC', 'OpenMM_OBC2']
# phases = ['sander_Gas', \
# 'sander_HCT', 'sander_OBC1', 'sander_OBC2', 'sander_GBn', 'sander_GBn2', \
# 'sander_PBSA', 'sander_ALPB_HCT', 'sander_ALPB_OBC1', 'sander_ALPB_OBC2', \
# 'sander_ALPB_GBn', 'gbnsr6_Still', 'gbnsr6_CHA']
self = AlGDock.BindingPMF_plots.BPMF_plots(\
dir_dock='dock', dir_cool='cool',\
ligand_database='prmtopcrd/ligand.db', \
forcefield='prmtopcrd/gaff2.dat', \
ligand_prmtop='prmtopcrd/ligand.prmtop', \
ligand_inpcrd='prmtopcrd/ligand.trans.inpcrd', \
ligand_mol2='prmtopcrd/ligand.mol2', \
ligand_rb='prmtopcrd/ligand.rb', \
receptor_database='prmtopcrd/receptor.db', \
receptor_prmtop='prmtopcrd/receptor.prmtop', \
receptor_inpcrd='prmtopcrd/receptor.trans.inpcrd', \
receptor_fixed_atoms='prmtopcrd/receptor.pdb', \
complex_prmtop='prmtopcrd/complex.prmtop', \
complex_inpcrd='prmtopcrd/complex.trans.inpcrd', \
complex_fixed_atoms='prmtopcrd/complex.pdb', \
score = 'prmtopcrd/xtal_plus_dock6_scored.mol2', \
temperature_scaling = 'Quadratic', \
pose = -1, \
rmsd=True, \
dir_grid='grids', \
protocol='Adaptive', cool_therm_speed=25.0, dock_therm_speed=0.25, \
T_HIGH=450.0, T_SIMMIN=300.0, T_TARGET=300.0, \
sampler='HMC', \
MCMC_moves=1, \
sampling_importance_resampling = True, \
solvation = 'Full', \
seeds_per_state=10, steps_per_seed=200, darts_per_seed=0, \
sweeps_per_cycle=50, snaps_per_cycle=25, attempts_per_sweep=100, \
steps_per_sweep=50, darts_per_sweep=0, \
cool_repX_cycles=3, dock_repX_cycles=4, \
site='Sphere', site_center=[1.7416, 1.7416, 1.7416], \
site_max_R=1.0, \
site_density=10., \
phases=phases, \
cores=-1, \
random_seed=-1, \
max_time=240, \
keep_intermediate=True)
In [2]:
(confs, Es) = self._get_confs_to_rescore(site=False, \
minimize=False, sort=False)
if self.params['dock']['pose']<len(confs):
starting_pose = np.copy(confs[self.params['dock']['pose']])
self.confs['dock']['starting_poses'] = [np.copy(starting_pose)]
else:
self._clear('dock')
self._store_infinite_f_RL()
raise Exception('Pose index greater than number of poses')
Xo = np.copy(self.universe.configuration().array)
self.universe.setConfiguration(Configuration(self.universe, starting_pose))
import AlGDock.RigidBodies
rb = AlGDock.RigidBodies.identifier(self.universe, self.molecule)
(TorsionRestraintSpecs, ExternalRestraintSpecs) = rb.poseInp()
self.universe.setConfiguration(Configuration(self.universe, Xo))
natoms = self.universe.numberOfAtoms()
In [3]:
TorsionRestraintSpecs
Out[3]:
In [4]:
self._set_universe_evaluator(
{'k_angular_ext':self.params['dock']['k_pose'], \
'k_spatial_ext':self.params['dock']['k_pose'], \
'k_angular_int':self.params['dock']['k_pose'],
'T':300.})
self.universe.energyTerms()
Out[4]:
In [5]:
BAT = rb.BAT(Xo, extended=True)
print BAT
In [6]:
%matplotlib inline
import matplotlib.pyplot as plt
In [7]:
# Check internal torsions
for n in range(len(rb._softTorsionInd)):
softTorsionInd = rb._softTorsionInd[n]
BAT_ind = np.array(range(6+5,natoms*3,3))[softTorsionInd]
Es = []
for torsion_offset in np.linspace(-2*np.pi,2*np.pi,100):
BAT_n = [BAT[ind] if ind!=BAT_ind else BAT[ind] + torsion_offset \
for ind in range(len(BAT))]
XYZ = rb.Cartesian(BAT_n)
self.universe.setConfiguration(Configuration(self.universe, XYZ))
Es.append(self.universe.energy())
# self.universe.energyAndGradients()
plt.plot(np.linspace(-2*np.pi,2*np.pi,100), Es)
In [8]:
# Check external translations
for n in range(3):
# print '*'*20
# print n
# print '*'*20
Es = []
BATs = []
for spatial_offset in np.linspace(-2.0,2.0,100):
BAT_n = [BAT[ind] if ind!=n else BAT[ind] + spatial_offset \
for ind in range(len(BAT))]
XYZ = rb.Cartesian(BAT_n)
self.universe.setConfiguration(Configuration(self.universe, XYZ))
Es.append(self.universe.energy())
BATs.append(rb.BAT(XYZ, extended=True)[:6])
# self.universe.energyAndGradients()
# terms = self.universe.energyTerms()
# print terms # ['pose external distance'], terms['pose external angle']
plt.plot(np.linspace(-2.0,2.0,100), Es)
# print '\n'.join(['\t'.join(['%3.2f'%b for b in BAT_b]) for BAT_b in BATs])
In [11]:
# Check external angles
for n in range(3,6):
Es = []
BATs = []
for angle_offset in np.linspace(-2*np.pi,2*np.pi,100):
BAT_n = [BAT[ind] if ind!=n else BAT[ind] + angle_offset \
for ind in range(len(BAT))]
XYZ = rb.Cartesian(BAT_n)
self.universe.setConfiguration(Configuration(self.universe, XYZ))
Es.append(self.universe.energy())
BATs.append(BAT_n[3:6] + list(rb.BAT(XYZ, extended=True)[3:6]))
# print self.universe.energyTerms()
# self.universe.energyAndGradients()
plt.plot(np.linspace(-2*np.pi,2*np.pi,100), Es)
# print '\n'.join(['\t'.join(['%3.2f'%b for b in BAT_b]) for BAT_b in BATs])
In [10]:
# The discontinuities in the potential as a function of theta are due to gimbal lock.
In [ ]: