In [1]:
import AlGDock.Nanopore
from AlGDock.Nanopore import *
# This will run a small molecule
# self = AlGDock.Nanopore.Simulation(\
# ligand_tarball='prmtopcrd/ligand.tar.gz', \
# ligand_database='ligand.db', \
# forcefield='prmtopcrd/gaff2.dat', \
# frcmodList=['ligand.frcmod'], \
# ligand_prmtop='ligand.prmtop', \
# starting_conf='prmtopcrd/anchor_and_grow_scored.mol2', \
# grid_LJr='grids/LJr.nc', \
# ef=1.0E8, \
# max_trials=10000, \
# report_interval=100)
# This will run a peptide
self = AlGDock.Nanopore.Simulation(\
ligand_database='prmtopcrd/YFF_14.db', \
forcefield='prmtopcrd/parm10.dat', \
frcmodList=['prmtopcrd/frcmod.ff14SB'], \
ligand_prmtop='prmtopcrd/YFF_14.prmtop', \
starting_conf='prmtopcrd/YFF_14.inpcrd', \
grid_LJr='grids/LJr.nc', \
ef=1.0E8, \
max_trials=10000, \
report_interval=100)
In [5]:
from AlGDock.Integrators.HamiltonianMonteCarlo.HamiltonianMonteCarlo \
import HamiltonianMonteCarloIntegrator
sampler = HamiltonianMonteCarloIntegrator(self.universe)
e_o = self.universe.energy()
T_LOW = 20.
T_START = 300.
delta_t = 1.5*MMTK.Units.fs
T_SERIES = T_LOW*(T_START/T_LOW)**(np.arange(30)/29.)
for T in T_SERIES:
attempts_left = 5
while attempts_left>0:
random_seed = int(T*10000) + attempts_left + \
int(self.universe.configuration().array[0][0]*10000) + \
int(time.time())
random_seed = random_seed%32767
(xs, energies, acc, ntrials, delta_t) = \
sampler(steps = 500, steps_per_trial = 50, T=T,\
delta_t=delta_t, random_seed=random_seed)
if (np.std(energies)>1E-3) and float(acc)/ntrials>0.4:
attempts_left = 0
else:
delta_t *= 0.9
attempts_left -= 1
self.universe.normalizePosition()
e_f = self.universe.energy()
In [2]:
print self.universe.energyTerms()
import AlGDock.IO
IO_dcd = AlGDock.IO.dcd(self.molecule,
ligand_atom_order = self.molecule.prmtop_atom_order)
IO_dcd.write(self.args['output_dcd'], [self.universe.configuration().array])
In [7]:
self.run()
self.view_trajectory()
In [ ]: