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)


    

    




The net charge on the ligand is too low for the electric field to have an effect
  grid type LJa not found
  grid type ELE not found






    




Warning: multiple database entries for
         /Users/dminh/Installers/AlGDock-0.0.1/Nanopore_Example/prmtopcrd/yff_14.db,
         using first one
/Users/dminh/Installers/AlGDock-0.0.1/Nanopore_Example/prmtopcrd/yff_14.db
/Users/dminh/Installers/AlGDock-0.0.1/Nanopore_Example/prmtopcrd/yff_14.db

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()


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-c66980b33177> in <module>()
     18       delta_t *= 0.9
     19       attempts_left -= 1
---> 20 if normalize:
     21   self.universe.normalizePosition()
     22 e_f = self.universe.energy()

NameError: name 'normalize' is not defined

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])


    

    




[[ 1.11958016  1.35494716  0.51392345]
 [ 0.89105967  1.2255089   0.23902993]
 [ 1.37473363  1.26809565  0.78205746]
 [ 1.02385584  1.43040381  0.60557165]
 [ 1.48263862  1.19371162  0.70408454]
 [ 0.99896454  1.15112458  0.16105714]
 [ 1.19928848  1.00588673  0.2091959 ]
 [ 1.0122318   1.59398586  0.79700522]
 [ 1.68622158  1.0470314   0.72816574]
 [ 1.06991898  1.09162877  0.39550768]
 [ 1.22991835  1.52855898  0.7125535 ]
 [ 1.57179624  1.12155002  0.92855042]
 [ 1.28521359  0.94103894  0.29870361]
 [ 1.07381842  1.67688912  0.89152584]
 [ 1.77894478  0.97771488  0.80689084]
 [ 1.1558444   1.02678084  0.48501498]
 [ 1.29150465  1.61146239  0.80707425]
 [ 1.66451961  1.05223334  1.00727522]
 [ 1.09028143  1.51982079  0.70751899]
 [ 1.58264729  1.11894893  0.7889952 ]
 [ 1.09164108  1.08118166  0.25759753]
 [ 1.26349163  0.95148591  0.43661301]
 [ 1.21345469  1.68562736  0.8965604 ]
 [ 1.76809386  0.98031589  0.94644543]
 [ 0.96044729  1.32368201  0.3323693 ]
 [ 1.21139695  1.27096785  0.60157042]
 [ 1.4441211   1.36626873  0.87539697]
 [ 0.80386075  1.392799    0.17046762]
 [ 0.70811084  1.25944789  0.17339027]
 [ 0.83716307  1.27897796  0.05691228]
 [ 1.18010452  1.42744133  0.45949562]
 [ 1.31932384  1.1960517   0.84222735]
 [ 0.83564975  1.15346514  0.29919993]
 [ 0.95890209  1.49425796  0.5456997 ]
 [ 1.4360076   1.12003261  0.63867854]
 [ 0.9523334   1.07744556  0.09565124]
 [ 0.96318532  1.3591293   0.66142889]
 [ 1.53957544  1.26465748  0.64403553]
 [ 1.05590148  1.22207026  0.10100802]
 [ 1.21620073  0.99775273  0.10182361]
 [ 0.90351494  1.58718282  0.79308575]
 [ 1.69467013  1.04500605  0.61951253]
 [ 0.98610801  1.15025134  0.4331921 ]
 [ 1.29068582  1.47081623  0.64288221]
 [ 1.49115636  1.17754293  0.97591077]
 [ 1.36902488  0.88241653  0.26101942]
 [ 1.01305125  1.73463186  0.96119747]
 [ 1.85958496  0.92172204  0.75953082]
 [ 1.13893225  1.03491467  0.59238734]
 [ 1.40022149  1.61826562  0.81099404]
 [ 1.65607139  1.05425853  1.11592847]
 [ 1.41712941  0.83986984  0.47964859]
 [ 1.2614043   1.75017339  0.97015146]
 [ 1.84028547  0.9263481   1.00773857]
 [ 1.06377736  1.17239561  0.42148766]
 [ 1.27565395  1.43701002  0.70060474]
 [ 0.79996389  1.29405637  0.14959377]
 [ 1.04713816  1.27200045  0.41975327]
 [ 1.28363773  1.3366431   0.69262142]
 [ 1.3469619   0.88849098  0.5235633 ]
 [ 1.52528789  1.3264516   0.95864926]
 [ 0.93615526  1.44396833  0.32562386]
 [ 1.21693074  1.1490299   0.58725552]
 [ 1.41982887  1.48655502  0.86865171]]
{'LJr': 0.0, 'electrostatic/pair sum': -274.61491342746382, 'electric_field_z': -58.887352253229238, 'OBC': -674.3284442242807, 'cosine dihedral angle': 73.517879515929096, 'harmonic bond': 5.6658131091403723, 'electrostatic': -274.61491342746382, 'harmonic bond angle': 21.901990970253035, 'Lennard-Jones': 139737.9529031023}

In [7]:

    

self.run()
self.view_trajectory()


    

    




Simulating motion from 0.60000 to 5.40000
Trial      0, COM: [1.237132955847144, 1.243401069823231, 0.59694940253328]
Trial    100, COM: [1.36759375843505, 1.3158893624066776, 0.5442941825330433]
Trial    200, COM: [1.323765543017124, 1.4527882871233795, 0.6102412025874019]
Trial    300, COM: [1.340212081902348, 1.5783592586566035, 0.635318345615737]
Trial    400, COM: [1.303706707348885, 1.6175472092431704, 0.6880280465716864]
Trial    500, COM: [1.189586321137778, 1.5993760872530935, 0.7175061823042246]
Trial    600, COM: [1.135882553137019, 1.4314473813365838, 0.721785043088248]
Trial    700, COM: [1.1465917389862876, 1.4703042035283256, 0.7015148141291522]
Trial    800, COM: [1.077318332499994, 1.6693473795120095, 0.7945644701497677]
Trial    900, COM: [1.0859161945116171, 1.6262982953568021, 0.9188541590150485]
Trial   1000, COM: [1.0242707209326343, 1.6405355178820074, 0.918859279085983]
Trial   1100, COM: [1.0748728541845847, 1.7544639327444125, 0.9235990579844834]
Trial   1200, COM: [1.1532398605691891, 1.747892587785487, 1.114110744062306]
Trial   1300, COM: [1.1294350923602514, 1.7689024102138649, 1.0671233681039847]






    




---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-7-3814cdf80868> in <module>()
----> 1 self.run()
      2 self.view_trajectory()

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/site-packages/AlGDock/Nanopore.pyc in run(self)
    223     confs = []
    224     for n in range(self.args['max_trials']):
--> 225       sampler(T=300.0*MMTK.Units.K, steps=100)
    226       if (n%100==0):
    227         print 'Trial %6d, COM:'%n, self.universe.centerOfMass()

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/site-packages/AlGDock/Integrators/HamiltonianMonteCarlo/HamiltonianMonteCarlo.pyc in __call__(self, **options)
     94             (self.universe,
     95              self.universe.configuration().array,
---> 96              self.universe.velocities().array) + late_args)
     97 
     98           # Decide whether to accept the move

/Users/dminh/Applications/miniconda2/envs/algdock/lib/python2.7/site-packages/MMTK/Trajectory.pyc in run(self, function, args)
   1034                                       function, args, self.state_accessor)
   1035         else:
-> 1036             apply(function, args)
   1037 
   1038 #

KeyboardInterrupt: 

In [ ]: