In [1]:
%matplotlib inline
import pymc
import seaborn
import numpy as np
import glob
import matplotlib as mpl
import matplotlib.pyplot as plt
from parmed.charmm import CharmmParameterSet
import torsionfit.TorsionFitModel as TorsionFitModel
import torsionfit.TorsionScanSet as TorsionScanSet
import torsionfit.utils as utils
import simtk.openmm as omm
import pymbar
In [2]:
param = CharmmParameterSet('../charmm_ff/top_all36_cgenff.rtf', '../charmm_ff/par_all36_cgenff.prm')
stream = 'pyrrol.str'
structure = 'pyrrol.psf'
scan = ['torsion-scan/PRL.scan2.neg.log',
'torsion-scan/PRL.scan4.pos.log',
'torsion-scan/PRL.scan2.pos.log',
'torsion-scan/PRL.scan4.neg.log',
'torsion-scan/PRL.scan3.neg.log',
'torsion-scan/PRL.scan5.pos.log',
'torsion-scan/PRL.scan3.pos.log',
'torsion-scan/PRL.scan5.neg.log']
#scan = glob.glob('torsion-scan/*.log')
pyrrol = TorsionScanSet.read_scan_logfile(scan, structure)
# extract only qm optimized structures
pyrrol_opt = pyrrol.extract_geom_opt()
In [3]:
# create pymc model
model = TorsionFitModel.TorsionFitModel(param, stream, pyrrol_opt)
mcmc = pymc.MCMC(model.pymc_parameters)
In [4]:
# import database
db = pymc.database.sqlite.load('phase.test.db.sqlite')
db2 = pymc.database.sqlite.load('pyrrol.database_3')
In [5]:
seaborn.tsplot(db.trace('mm_energy')[100:], err_style='unit_traces'), plt.plot(pyrrol_opt.qm_energy, color='r')
Out[5]:
In [6]:
sigma_plot = plt.plot(np.exp(db.trace('log_sigma')[:]), label = 'sigma')
plt.legend()
Out[6]:
In [7]:
error = (db.trace('mm_energy')[:,:] - np.tile(pyrrol_opt.qm_energy, [10000,1]))
rms_error = np.sqrt((error**2).mean(1))
In [8]:
plt.plot(rms_error[100:],'.', label='RMSE')
Out[8]:
In [9]:
deviance = plt.plot(db.trace('deviance')[100:], label = 'deviance')
In [42]:
multiplicities = (1,2,3,4,6)
torsion_parameters = {}
multiplicity_traces = {}
torsions = model.parameters_to_optimize
for name in torsions:
torsion_name = name[0] + '_' + name[1] + '_' + name[2] + '_' + name[3]
torsion_parameters[torsion_name] = []
multiplicity_bitstring = torsion_name + '_multiplicity_bitstring'
torsion_parameters[torsion_name].append(multiplicity_bitstring)
for m in multiplicities:
k = torsion_name + '_' + str(m) + '_K'
torsion_parameters[torsion_name].append(k)
phase = torsion_name + '_' + str(m) + '_Phase'
torsion_parameters[torsion_name].append(phase)
multiplicity_traces[torsion_name + '_' + str(m)] = []
for i in db._traces[multiplicity_bitstring][:]:
if 2**(m-1) & int(i):
multiplicity_traces[torsion_name + '_' + str(m)].append(1)
else:
multiplicity_traces[torsion_name + '_' + str(m)].append(0)
In [95]:
for parameters in torsion_parameters:
for param in torsion_parameters[parameters]:
if param[-1] == 'K':
statistics[param] = pymbar.timeseries.detectEquilibration(db._traces[param][:])
In [157]:
import pandas as pd
In [186]:
stats = pd.DataFrame(data=statistics, index=['t', 'g: ' r'$1+2\tau$', 'Neff_max'])
stats.to_csv('equil_detection_output.csv')
In [156]:
for i in multiplicities:
print 'CG251O_NG3C51_CG2R51_CG2R51_' + str(i) + '_K'
print(statistics['CG251O_NG3C51_CG2R51_CG2R51_' + str(i) + '_K'])
print 'CG251O_NG3C51_CG2R51_CG2R51_' + str(i) + '_Phase'
print(statistics['CG251O_NG3C51_CG2R51_CG2R51_' + str(i) + '_Phase'])
In [131]:
k = []
phase = []
multiplicity = []
for i in multiplicities:
k.append([db._traces['CG251O_NG3C51_CG2R51_CG2R51_' + str(i) + '_K'][:]])
phase.append([db._traces['CG251O_NG3C51_CG2R51_CG2R51_' + str(i) + '_Phase'][:]])
multiplicity.append([multiplicity_traces['CG251O_NG3C51_CG2R51_CG2R51_' + str(i)]])
In [146]:
axes_k = plt.subplot(3, 1, 1)
plt.plot(k[0][0])
plt.title('CG251O_NG3C51_CG2R51_CG2R51_1')
axes_k.axvline((2070), color='red', lw=3)
plt.legend('effective_samples = 758')
plt.subplot(3, 1, 2)
plt.plot(phase[0][0])
plt.legend('Phase')
plt.subplot(3, 1, 3)
plt.plot(multiplicity[0][0])
plt.legend('multiplicity')
plt.show()
In [148]:
axes_2 = plt.subplot(3, 1, 1)
plt.plot(k[1][0])
plt.title('CG251O_NG3C51_CG2R51_CG2R51_2')
axes_2.axvline((428), color='red', lw=3)
plt.legend('K', 'effective_samples = 758')
plt.subplot(3, 1, 2)
plt.plot(phase[1][0])
plt.legend('Phase')
plt.subplot(3, 1, 3)
plt.plot(multiplicity[1][0])
plt.legend('multiplicity')
plt.show()
In [150]:
axes_3 = plt.subplot(3, 1, 1)
plt.plot(k[2][0])
plt.title('CG251O_NG3C51_CG2R51_CG2R51_3')
axes_3.axvline((3299), color='red', lw=3)
plt.legend('K', 'effective_samples = 758')
plt.subplot(3, 1, 2)
plt.plot(phase[2][0])
plt.legend('Phase')
plt.subplot(3, 1, 3)
plt.plot(multiplicity[2][0])
plt.legend('multiplicity')
plt.show()
In [151]:
axes_4 = plt.subplot(3, 1, 1)
plt.plot(k[3][0])
plt.title('CG251O_NG3C51_CG2R51_CG2R51_4')
axes_4.axvline((7677), color='red', lw=3)
plt.legend('K')
plt.subplot(3, 1, 2)
plt.plot(phase[3][0])
plt.legend('Phase')
plt.subplot(3, 1, 3)
plt.plot(multiplicity[3][0])
plt.legend('multiplicity')
plt.show()
In [153]:
axes_5 = plt.subplot(3, 1, 1)
plt.plot(k[4][0])
plt.title('CG251O_NG3C51_CG2R51_CG2R51_6')
axes_5.axvline((7373), color='red', lw=3)
plt.legend('K')
plt.subplot(3, 1, 2)
plt.plot(phase[4][0])
plt.legend('Phase')
plt.subplot(3, 1, 3)
plt.plot(multiplicity[4][0])
plt.legend('multiplicity')
Out[153]:
In [ ]: