In [1]:

    

import AlGDock.BindingPMF_plots
from AlGDock.BindingPMF_plots import *
import os, shutil, glob

phases = ['NAMD_Gas', 'NAMD_OBC']

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_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, 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)


    

    




###########
# AlGDock #
###########
Molecular docking with adaptively scaled alchemical interaction grids

in /Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/site-packages/AlGDock/BindingPMF.py
last modified Fri Jul  7 04:39:28 2017
    
using 4/4 available cores
using random number seed of -1

*** Directories ***
  start: /Users/dminh/Installers/AlGDock-0.0.1/Example
  cool: /Users/dminh/Installers/AlGDock-0.0.1/Example/cool
  dock: /Users/dminh/Installers/AlGDock-0.0.1/Example/dock

*** Files ***
previously stored in cool directory:
  prmtop:
    L: prmtopcrd/ligand.prmtop
  ligand_database: prmtopcrd/ligand.db
  forcefield: ../Data/gaff2.dat
  inpcrd:
    L: prmtopcrd/ligand.trans.inpcrd

previously stored in dock directory:
  namd: ../../../Applications/NAMD_2.10/namd2
  grids:
    LJr: grids/LJr.nc
    LJa: grids/LJa.nc
    ELE: grids/pbsa.nc
    desolv: grids/desolv.nc
  prmtop:
    L: prmtopcrd/ligand.prmtop
    R: prmtopcrd/receptor.prmtop
    RL: prmtopcrd/complex.prmtop
  ligand_database: prmtopcrd/ligand.db
  fixed_atoms:
    R: prmtopcrd/receptor.pdb
    RL: prmtopcrd/complex.pdb
  mol2:
    L: prmtopcrd/ligand.mol2
  score: prmtopcrd/xtal_plus_dock6_scored.mol2
  forcefield: ../Data/gaff2.dat
  dir_cool: cool
  inpcrd:
    L: prmtopcrd/ligand.trans.inpcrd
    R: prmtopcrd/receptor.trans.inpcrd
    RL: prmtopcrd/complex.trans.inpcrd

from arguments and defaults:
  failed to find file in: prmtopcrd/ligand.rb, /Users/dminh/Installers/AlGDock-0.0.1/Example/dock/prmtopcrd/ligand.rb
  ligand_database: prmtopcrd/ligand.db
  forcefield: ../Data/gaff2.dat
  prmtop:
    L: prmtopcrd/ligand.prmtop
    R: prmtopcrd/receptor.prmtop
    RL: prmtopcrd/complex.prmtop
  inpcrd:
    L: prmtopcrd/ligand.trans.inpcrd
    R: prmtopcrd/receptor.trans.inpcrd
    RL: prmtopcrd/complex.trans.inpcrd
  mol2:
    L: prmtopcrd/ligand.mol2
  fixed_atoms:
    R: prmtopcrd/receptor.pdb
    RL: prmtopcrd/complex.pdb
  grids:
    LJr: grids/LJr.nc
    LJa: grids/LJa.nc
    ELE: grids/pbsa.nc
    desolv: grids/desolv.nc
  score: prmtopcrd/xtal_plus_dock6_scored.mol2
  dir_cool: cool

to be used:
  ligand_database: prmtopcrd/ligand.db
  forcefield: ../Data/gaff2.dat
  frcmodList: ['/Users/dminh/Installers/AlGDock-0.0.1/Example/prmtopcrd/ligand.frcmod']
  prmtop:
    L: prmtopcrd/ligand.prmtop
    R: prmtopcrd/receptor.prmtop
    RL: prmtopcrd/complex.prmtop
  inpcrd:
    L: prmtopcrd/ligand.trans.inpcrd
    R: prmtopcrd/receptor.trans.inpcrd
    RL: prmtopcrd/complex.trans.inpcrd
  mol2:
    L: prmtopcrd/ligand.mol2
  fixed_atoms:
    R: prmtopcrd/receptor.pdb
    RL: prmtopcrd/complex.pdb
  grids:
    LJr: grids/LJr.nc
    LJa: grids/LJa.nc
    ELE: grids/pbsa.nc
    desolv: grids/desolv.nc
  score: prmtopcrd/xtal_plus_dock6_scored.mol2
  dir_cool: cool
  namd: ../../../Applications/NAMD_2.10/namd2

>>> Setting up the simulation
  considering a residue named "iga" as the ligand

*** Simulation parameters and constants ***

for cool:
  protocol: Adaptive
  therm_speed: 25.0
  T_HIGH: 450.0
  T_SIMMIN: 300.0
  T_TARGET: 300.0
  H_mass: 4.0
  fraction_CD: 0.5
  CD_steps_per_trial: 5
  delta_t: 3.0
  sampler: HMC
  steps_per_seed: 200
  seeds_per_state: 10
  darts_per_seed: 0
  repX_cycles: 3
  min_repX_acc: 0.4
  sweeps_per_cycle: 50
  attempts_per_sweep: 100
  steps_per_sweep: 50
  darts_per_sweep: 0
  snaps_per_independent: 3.0
  phases: ['NAMD_Gas', 'NAMD_OBC']
  sampling_importance_resampling: True
  solvation: Full
  keep_intermediate: True
  GMC_attempts: 0
  GMC_tors_threshold: 0.0
  delta_t_CD: 4.0

for dock:
  protocol: Adaptive
  therm_speed: 0.25
  T_HIGH: 450.0
  T_SIMMIN: 300.0
  T_TARGET: 300.0
  H_mass: 4.0
  fraction_CD: 0.5
  CD_steps_per_trial: 5
  delta_t: 3.0
  sampler: HMC
  steps_per_seed: 200
  seeds_per_state: 10
  darts_per_seed: 0
  repX_cycles: 4
  min_repX_acc: 0.4
  sweeps_per_cycle: 50
  attempts_per_sweep: 100
  steps_per_sweep: 50
  darts_per_sweep: 0
  snaps_per_independent: 20.0
  phases: ['NAMD_Gas', 'NAMD_OBC']
  sampling_importance_resampling: True
  solvation: Full
  keep_intermediate: True
  GMC_attempts: 0
  GMC_tors_threshold: 0.0
  temperature_scaling: Quadratic
  site: Sphere
  site_center: array([ 1.7416,  1.7416,  1.7416])
  site_max_R: 1.0
  site_density: 10.0
  pose: -1
  k_pose: 200.0
  MCMC_moves: 1
  rmsd: True
  receptor_NAMD_Gas: array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])
  receptor_NAMD_OBC: array([[     0.      ,      0.      ,      0.      ,      0.      ,
        -56747.052264,      0.      ,      0.      , -56747.052264]])
  delta_t_CD: 4.0

In [4]:

    

minimize=False

    # Determine the name of the file
    prefix = 'xtal' if self._FNs['score']=='default' else \
      os.path.basename(self._FNs['score']).split('.')[0]
    if minimize:
      prefix = 'min_' + prefix
    energyFN = os.path.join(self.dir['dock'],prefix+'.pkl.gz')

    # Set the force field to fully interacting
    lambda_o = self._lambda(1.0, 'dock')
    self._set_universe_evaluator(lambda_o)

    # Load the configurations
    if os.path.isfile(energyFN):
      (confs, Es) = self._load_pkl_gz(energyFN)
    else:
      (confs, Es) = self._get_confs_to_rescore(site=False, \
        minimize=minimize, sort=False)


    

    




  LJr grid loaded from LJr.nc in 0.40 s
  LJa grid loaded from LJa.nc in 0.16 s
  ELE grid loaded from pbsa.nc in 0.03 s
  keeping 693 configurations out of
  0 from xtal, 693 from dock6, 0 from initial docking, and 0 duplicated

In [5]:

    

# Calculate MM energies
    if not 'MM' in Es.keys():
      Es = self._energyTerms(confs, Es)

    # Calculate RMSD
    if (self.params['dock']['rmsd'] is not False):
      Es['rmsd'] = np.array([np.sqrt(((confs[c][self.molecule.heavy_atoms,:] - \
        self.confs['rmsd'])**2).sum()/self.molecule.nhatoms) \
          for c in range(len(confs))])


    

    




  sLJr grid loaded from LJr.nc in 0.29 s
  sELE grid loaded from pbsa.nc in 0.02 s

In [13]:

    

%matplotlib inline
import matplotlib.pyplot as plt


    

In [10]:

    

E_test = self._energyTerms(confs, {})


    

In [12]:

    

E_test.keys()


    

    
Out[12]:





['MM', 'sLJr', 'OBC', 'ELE', 'LJa', 'sELE', 'site', 'LJr']

In [14]:

    

Es.keys()


    

    
Out[14]:





['Grid_vdw',
 'Cluster Size',
 'HA_RMSDh',
 'LJr',
 'rmsd',
 'Int_energy',
 'MM',
 'Grid Score',
 'HA_RMSDm',
 'Grid_es',
 'OBC',
 'ELE',
 'LJa',
 'sLJr',
 'site',
 'HA_RMSDs',
 'sELE',
 'total']

In [40]:

    

scores = {}
for solvation in ['Desolvated','Fractional','Full']:
  self.params['dock']['solvation'] = solvation
  lambda_o = self._lambda(1.0, 'dock')
  self._set_universe_evaluator(lambda_o)
  Eterms = self._energyTerms(confs, Es)
  scores[solvation] = self._u_kln([Eterms], [lambda_o])


    

In [41]:

    

scores['Full']-scores['Fractional']


    

    
Out[41]:





array([ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.])

In [ ]: