In [1]:
# opemm imports
import simtk.unit as u
import simtk.openmm as mm
import simtk.openmm.app as app
In [2]:
# ParmEd imports
from chemistry.charmm import CharmmPsfFile
from chemistry.charmm.parameters import CharmmParameterSet
In [4]:
import TorsionScanSet, TorsionFitModel
import glob
from pymc import MCMC
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-4-f6fa1eb33d3c> in <module>()
----> 1 import TorsionScanSet, TorsionFitModel
2 import glob
3 from pymc import MCMC
4 import matplotlib
5 import matplotlib.pyplot as plt
/Users/sternc1/src/ChayaSt/Torsions/torsions/TorsionFitModel.py in <module>()
5
6
----> 7 class TorsionFitModel(object):
8 """pymc model
9
/Users/sternc1/src/ChayaSt/Torsions/torsions/TorsionFitModel.py in TorsionFitModel()
64 lambda log_sigma=self.pymc_parameters['log_sigma']: np.exp(-2*log_sigma))
65
---> 66 @pymc.deterministic
67 def mm_potential(self, param):
68 self.update_param(param=param)
/Users/sternc1/anaconda/envs/py3k/lib/python3.4/site-packages/pymc/InstantiationDecorators.py in deterministic(__func__, **kwds)
248
249 if __func__:
--> 250 return instantiate_n(__func__)
251
252 return instantiate_n
/Users/sternc1/anaconda/envs/py3k/lib/python3.4/site-packages/pymc/InstantiationDecorators.py in instantiate_n(__func__)
240 def instantiate_n(__func__):
241 junk, parents = _extract(
--> 242 __func__, kwds, keys, 'Deterministic', probe=False)
243 return Deterministic(parents=parents, **kwds)
244
/Users/sternc1/anaconda/envs/py3k/lib/python3.4/site-packages/pymc/InstantiationDecorators.py in _extract(__func__, kwds, keys, classname, probe)
113 if i < arg_deficit - 1:
114 err_str += ','
--> 115 raise ValueError(err_str)
116
117 # Fill in parent dictionary
ValueError: Deterministic mm_potential: no parent provided for the following labels: self, param
In [4]:
param = CharmmParameterSet('../charmm_ff/top_all36_cgenff.rtf', '../charmm_ff/par_all36_cgenff.prm')
stream = '../structures/Pyrrol/pyrrol.str'
structure = '../structures/Pyrrol/pyrrol.psf'
#create pymc model
model = TorsionFitModel.TorsionFitModel(param, stream, ['pyrrol'])
mcmc = MCMC(model)
sample = mcmc.sample(iter=10000, burn=1000, thin=10)
[-----------------100%-----------------] 10000 of 10000 complete in 80.7 sec
In [5]:
plt.hist(mcmc.trace('CG251O_NG3C51_CG2R51_CG2R51_1_Phase')[:])
Out[5]:
(array([ 456., 0., 0., 0., 0., 0., 0., 0., 0., 444.]),
array([ 0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
<a list of 10 Patch objects>)
In [6]:
plt.hist(mcmc.trace('CG251O_NG3C51_CG2R51_CG2R51_1_K')[:])
Out[6]:
(array([ 95., 94., 79., 102., 82., 91., 104., 70., 86., 97.]),
array([ 0.02376539, 2.02082243, 4.01787947, 6.01493651,
8.01199356, 10.0090506 , 12.00610764, 14.00316469,
16.00022173, 17.99727877, 19.99433582]),
<a list of 10 Patch objects>)
In [7]:
plt.hist(mcmc.trace('CG251O_NG3C51_CG2R51_CG2R51_multiplicity_bitstring')[:])
Out[7]:
(array([ 110., 95., 70., 93., 94., 67., 111., 93., 77., 90.]),
array([ 0. , 6.3, 12.6, 18.9, 25.2, 31.5, 37.8, 44.1, 50.4,
56.7, 63. ]),
<a list of 10 Patch objects>)
In [8]:
plt.plot(mcmc.trace('CG251O_NG3C51_CG2R51_CG2R51_multiplicity_bitstring')[:])
Out[8]:
[<matplotlib.lines.Line2D at 0x113868e10>]
In [16]:
# Torsion Scan
logfiles = glob.glob('../structures/Pyrrol/torsion-scan/*.log')
scan = TorsionScanSet.read_scan_logfile(logfiles, structure)
# Remove all non optimized geometries
scan_opt = scan.extract_geom_opt()
loading ../structures/Pyrrol/torsion-scan/PRL.scan2.neg.log
loading ../structures/Pyrrol/torsion-scan/PRL.scan2.pos.log
loading ../structures/Pyrrol/torsion-scan/PRL.scan3.neg.log
loading ../structures/Pyrrol/torsion-scan/PRL.scan3.pos.log
loading ../structures/Pyrrol/torsion-scan/PRL.scan4.neg.log
loading ../structures/Pyrrol/torsion-scan/PRL.scan4.pos.log
loading ../structures/Pyrrol/torsion-scan/PRL.scan5.neg.log
loading ../structures/Pyrrol/torsion-scan/PRL.scan5.pos.log
In [17]:
# compute mm energy
scan_opt.compute_energy_from_positions(param)
scan_opt.to_dataframe()
Out[17]:
torsion
scan_direction
step_point_total
QM_energy KJ/mol
MM_energy KJ/mole
delta KJ/mole
0
[2, 5, 6, 7]
0
[2, 1, 7]
-964629.156046 kJ/mol
43.5042463908 kJ/mol
-964672.660293 kJ/mol
1
[2, 5, 6, 7]
0
[6, 2, 7]
-964627.250844 kJ/mol
46.1033652642 kJ/mol
-964673.35421 kJ/mol
2
[2, 5, 6, 7]
0
[5, 3, 7]
-964622.329681 kJ/mol
51.5999322751 kJ/mol
-964673.929614 kJ/mol
3
[2, 5, 6, 7]
0
[9, 4, 7]
-964616.653543 kJ/mol
56.5975354875 kJ/mol
-964673.251078 kJ/mol
4
[2, 5, 6, 7]
0
[6, 5, 7]
-964613.23588 kJ/mol
57.7700106875 kJ/mol
-964671.005891 kJ/mol
5
[2, 5, 6, 7]
0
[6, 6, 7]
-964614.826067 kJ/mol
55.0908226196 kJ/mol
-964669.916889 kJ/mol
6
[2, 5, 6, 7]
0
[6, 7, 7]
-964619.578995 kJ/mol
51.8808623717 kJ/mol
-964671.459858 kJ/mol
7
[2, 5, 6, 7]
1
[2, 1, 7]
-964629.156046 kJ/mol
43.5042463908 kJ/mol
-964672.660293 kJ/mol
8
[2, 5, 6, 7]
1
[7, 2, 7]
-964627.144126 kJ/mol
45.1683294131 kJ/mol
-964672.312455 kJ/mol
9
[2, 5, 6, 7]
1
[9, 3, 7]
-964622.217927 kJ/mol
50.9228599418 kJ/mol
-964673.140787 kJ/mol
10
[2, 5, 6, 7]
1
[6, 4, 7]
-964617.110361 kJ/mol
57.8374410231 kJ/mol
-964674.947802 kJ/mol
11
[2, 5, 6, 7]
1
[5, 5, 7]
-964615.038724 kJ/mol
62.9677195648 kJ/mol
-964678.006444 kJ/mol
12
[2, 5, 6, 7]
1
[6, 6, 7]
-964617.133104 kJ/mol
58.6760593568 kJ/mol
-964675.809163 kJ/mol
13
[2, 5, 6, 7]
1
[5, 7, 7]
-964622.19231 kJ/mol
52.8282586919 kJ/mol
-964675.020569 kJ/mol
14
[3, 2, 5, 6]
0
[2, 1, 7]
-964629.156033 kJ/mol
43.5053580595 kJ/mol
-964672.661391 kJ/mol
15
[3, 2, 5, 6]
0
[9, 2, 7]
-964628.536061 kJ/mol
44.5588922461 kJ/mol
-964673.094953 kJ/mol
16
[3, 2, 5, 6]
0
[4, 3, 7]
-964625.610193 kJ/mol
52.4528949674 kJ/mol
-964678.063088 kJ/mol
17
[3, 2, 5, 6]
0
[15, 4, 7]
-964622.622874 kJ/mol
73.5415501958 kJ/mol
-964696.164424 kJ/mol
18
[3, 2, 5, 6]
0
[5, 5, 7]
-964623.470248 kJ/mol
74.9560272558 kJ/mol
-964698.426276 kJ/mol
19
[3, 2, 5, 6]
0
[6, 6, 7]
-964623.912561 kJ/mol
73.7474625224 kJ/mol
-964697.660024 kJ/mol
20
[3, 2, 5, 6]
0
[5, 7, 7]
-964623.99701 kJ/mol
72.3960520045 kJ/mol
-964696.393062 kJ/mol
21
[3, 2, 5, 6]
1
[2, 1, 7]
-964629.156033 kJ/mol
43.5053580595 kJ/mol
-964672.661391 kJ/mol
22
[3, 2, 5, 6]
1
[9, 2, 7]
-964628.620725 kJ/mol
51.9533468434 kJ/mol
-964680.574072 kJ/mol
23
[3, 2, 5, 6]
1
[5, 3, 7]
-964626.630785 kJ/mol
65.7336098373 kJ/mol
-964692.364395 kJ/mol
24
[3, 2, 5, 6]
1
[6, 4, 7]
-964621.979177 kJ/mol
81.9840974038 kJ/mol
-964703.963275 kJ/mol
25
[3, 2, 5, 6]
1
[8, 5, 7]
-964613.183066 kJ/mol
99.3539395097 kJ/mol
-964712.537005 kJ/mol
26
[3, 2, 5, 6]
1
[9, 6, 7]
-964604.945784 kJ/mol
109.471779765 kJ/mol
-964714.417563 kJ/mol
27
[3, 2, 5, 6]
1
[8, 7, 7]
-964599.090568 kJ/mol
108.925492342 kJ/mol
-964708.01606 kJ/mol
28
[1, 2, 5, 6]
0
[2, 1, 7]
-964629.155941 kJ/mol
43.5044962267 kJ/mol
-964672.660437 kJ/mol
29
[1, 2, 5, 6]
0
[8, 2, 7]
-964628.014791 kJ/mol
50.213269356 kJ/mol
-964678.228061 kJ/mol
30
[1, 2, 5, 6]
0
[6, 3, 7]
-964625.594207 kJ/mol
57.7380604165 kJ/mol
-964683.332267 kJ/mol
31
[1, 2, 5, 6]
0
[5, 4, 7]
-964623.421102 kJ/mol
63.0697351262 kJ/mol
-964686.490837 kJ/mol
32
[1, 2, 5, 6]
0
[6, 5, 7]
-964622.358851 kJ/mol
67.0631558503 kJ/mol
-964689.422007 kJ/mol
33
[1, 2, 5, 6]
0
[7, 6, 7]
-964622.844125 kJ/mol
69.7917194045 kJ/mol
-964692.635844 kJ/mol
34
[1, 2, 5, 6]
0
[7, 7, 7]
-964623.811605 kJ/mol
72.0023810499 kJ/mol
-964695.813986 kJ/mol
35
[1, 2, 5, 6]
1
[2, 1, 7]
-964629.155941 kJ/mol
43.5044962267 kJ/mol
-964672.660437 kJ/mol
36
[1, 2, 5, 6]
1
[8, 2, 7]
-964627.462848 kJ/mol
49.4051977338 kJ/mol
-964676.868046 kJ/mol
37
[1, 2, 5, 6]
1
[6, 3, 7]
-964621.454364 kJ/mol
67.631821582 kJ/mol
-964689.086185 kJ/mol
38
[1, 2, 5, 6]
1
[7, 4, 7]
-964611.681792 kJ/mol
91.6266887133 kJ/mol
-964703.30848 kJ/mol
39
[1, 2, 5, 6]
1
[11, 5, 7]
-964602.246867 kJ/mol
95.3930775935 kJ/mol
-964697.639945 kJ/mol
40
[1, 2, 5, 6]
1
[9, 6, 7]
-964599.330244 kJ/mol
90.7388291669 kJ/mol
-964690.069073 kJ/mol
41
[1, 2, 5, 6]
1
[7, 7, 7]
-964599.581725 kJ/mol
90.5489223346 kJ/mol
-964690.130647 kJ/mol
42
[8, 1, 2, 5]
0
[2, 1, 7]
-964629.156062 kJ/mol
43.5039054453 kJ/mol
-964672.659967 kJ/mol
43
[8, 1, 2, 5]
0
[7, 2, 7]
-964628.120862 kJ/mol
49.1749208948 kJ/mol
-964677.295782 kJ/mol
44
[8, 1, 2, 5]
0
[6, 3, 7]
-964624.063945 kJ/mol
69.0303414596 kJ/mol
-964693.094286 kJ/mol
45
[8, 1, 2, 5]
0
[5, 4, 7]
-964616.670102 kJ/mol
99.0790731122 kJ/mol
-964715.749175 kJ/mol
46
[8, 1, 2, 5]
0
[4, 5, 7]
-964606.724381 kJ/mol
136.148721806 kJ/mol
-964742.873103 kJ/mol
47
[8, 1, 2, 5]
0
[5, 6, 7]
-964594.6486 kJ/mol
175.236597354 kJ/mol
-964769.885198 kJ/mol
48
[8, 1, 2, 5]
0
[5, 7, 7]
-964580.228057 kJ/mol
212.225905417 kJ/mol
-964792.453962 kJ/mol
49
[8, 1, 2, 5]
1
[2, 1, 7]
-964629.156062 kJ/mol
43.5039054453 kJ/mol
-964672.659967 kJ/mol
50
[8, 1, 2, 5]
1
[9, 2, 7]
-964628.227935 kJ/mol
56.3025307582 kJ/mol
-964684.530465 kJ/mol
51
[8, 1, 2, 5]
1
[6, 3, 7]
-964624.480874 kJ/mol
80.7821608896 kJ/mol
-964705.263035 kJ/mol
52
[8, 1, 2, 5]
1
[7, 4, 7]
-964617.060109 kJ/mol
113.278043483 kJ/mol
-964730.338153 kJ/mol
53
[8, 1, 2, 5]
1
[13, 5, 7]
-964607.998032 kJ/mol
149.214900082 kJ/mol
-964757.212932 kJ/mol
54
[8, 1, 2, 5]
1
[6, 6, 7]
-964597.645382 kJ/mol
188.025598544 kJ/mol
-964785.670981 kJ/mol
55
[8, 1, 2, 5]
1
[5, 7, 7]
-964585.572484 kJ/mol
224.818459925 kJ/mol
-964810.390944 kJ/mol
In [ ]:
Content source: hainm/Torsions
Similar notebooks: