In [2]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
import time
import itertools
import h5py
import numpy as np
import multiprocessing as mp
from scipy.stats import norm
from scipy.stats import expon
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
sns.set(style="ticks", color_codes=True, font_scale=1.5)
sns.set_style({"xtick.direction": "in", "ytick.direction": "in"})
In [3]:
import cossio
import kinetics
In [4]:
def pmf2d(xq, qk):
H, x_edges, y_edges = np.histogram2d(data[:,1],data[:,2], \
bins=[np.linspace(-12,12,100), np.linspace(-12,12,100)])
fig, ax = plt.subplots(figsize=(6,5))
pmf = -np.log(H.transpose())
pmf -= np.min(pmf)
#pmf = np.max(pmf) - pmf
cs = ax.contourf(pmf, extent=[x_edges.min(), x_edges.max(), \
y_edges.min(), y_edges.max()], \
cmap=cm.rainbow, levels=np.arange(0, 8,0.5))
#pmf = -np.log(counts.transpose())
#cs = ax.contourf(counts.transpose()) #, \
#extent=[xbins.min(), xbins.max(), ybins.min(), ybins.max()], \
#cmap=cm.rainbow, levels=np.arange(0,10 ,1))
cbar = plt.colorbar(cs)
ax.set_xlim(-12,12)
ax.set_ylim(-12,12)
ax.set_xlabel('$x$', fontsize=20)
ax.set_ylabel('$q$', fontsize=20)
plt.tight_layout()
In [5]:
def calc_rates(y):
lifeA, lifeB = kinetics.calc_life([y])
meanA = 1./np.exp(np.mean(np.log([x for x in lifeA if x>0])))
meanB = 1./np.exp(np.mean(np.log([x for x in lifeB if x>0])))
errorA = meanA/np.sqrt(len(lifeA))
errorB = meanA/np.sqrt(len(lifeB))
return np.mean([meanA, meanB]), np.mean([errorA, errorB])
In [6]:
def cossio_runner(inp):
np.random.seed()
kl = inp[0]
sc = inp[1]
numsteps = inp[2]
Dq = sc*Dx
x, q = [5., 5.]
tt, xk, qk = cossio.run_brownian(x0=x, barrier=5., kl=kl, \
Dx=Dx, Dq=Dq, numsteps=numsteps, \
fwrite=int(1./dt))
data = np.column_stack((tt,xk,qk))
h5file = "data/cossio_kl%g_Dx%g_Dq%g.h5"%(kl, Dx, Dq)
with h5py.File(h5file, "w") as hf:
hf.create_dataset("data", data=data)
return h5file
In [7]:
# Globals
dt = 5e-4
Dx = 1. # Diffusion coefficient for molecular coordinate
In [154]:
x = np.linspace(-10,10,1000)
fig, ax = plt.subplots(2,1, figsize=(7,5), sharex=True)
Gx = [cossio.Gx(y, barrier=5.) for y in x]
dGqxdx = [cossio.dGqxdx(0, y, barrier=5.) for y in x]
ax[0].plot(x, Gx)
ax[0].set_ylabel('$G(x)$', fontsize=20)
ax[0].set_ylim(-1.1*5,0.5*5)
ax[1].plot(x, dGqxdx)
ax[1].set_xlabel('$x$', fontsize=20)
ax[1].set_ylabel('$\partial G(x)/\partial x$', fontsize=20)
ax[1].hlines(0, -10, 10, linestyle='dashed', linewidth=0.5)
plt.tight_layout(h_pad=0)
In [117]:
x = np.linspace(-8,8,1000)
fig, ax = plt.subplots(2,1, figsize=(7,5), sharex=True)
Gx = np.array([cossio.Gx(y, barrier=5.) for y in x])
ax[0].plot(x, Gx)
ax[0].plot(x[450:550], Gx[450:550])
ax[0].set_ylabel('$G(x)$', fontsize=20)
ax[0].set_ylim(-1.1*5,1)
Px = np.exp(-np.array(Gx))
ax[1].plot(x, Px)
ax[1].plot(x[450:550], Px[450:550])
ax[1].set_xlabel('$x$', fontsize=20)
ax[1].set_ylabel('$P(x)$', fontsize=20)
ax[1].hlines(0, -10, 10, linestyle='dashed', linewidth=0.5)
plt.tight_layout(h_pad=0)
In [58]:
print "Barrier height: %g kBT"%-np.log(np.trapz(Px[430:570],x[430:570])/np.trapz(Px, x))
print "Barrier height: %g kBT"%-np.log(np.trapz(Px[430:570],x[430:570])/np.trapz(Px, x))
In [293]:
fig, ax = plt.subplots(3,1, figsize=(7,7), sharex=True)
x = np.linspace(-8,8,1000)
Gx = np.array([cossio.Gx(y, barrier=5.) for y in x])
ax[0].plot(x, Gx)
ax[0].set_ylim(-1.1*5,1)
Px = np.exp(-np.array(Gx))
ax12 = ax[0].twinx()
ax12.plot(x, Px, 'g')
i2 = np.argmin(np.abs(x-2.5))
i1 = np.argmin(np.abs(x+2.5))
print "x1 : ", i1, x[i1]
print "x2 :", i2, x[i2]
pfold = np.array([np.trapz(np.exp(Gx[i:i2]),x[i:i2])/np.trapz(np.exp(Gx[i1:i2]), x[i1:i2]) \
if x[i]>-2.5 else 1. for i in range(len(x))])
#ax.set_ylim(-1,1)
Px = np.exp(-np.array(Gx))
ax[1].plot(x, pfold, label=r"$\phi_{A}(x)$")
ax[1].plot(x, 1.-pfold, label=r"$\phi_{B}(x)$")
ax[1].set_ylim(-0.1,1.1)
ax22 = ax[1].twinx()
ax22.plot(x, pfold*(1.-pfold), 'r', label=r"$\phi_{A}(x)\phi_{B}(x)$")
ax22.set_ylim(-0.025,0.3)
ax[2].plot(x,Px*pfold*(1-pfold)/np.trapz(Px,x))
ax[2].set_xlabel('$x$', fontsize=20)
for i in range(3):
ax[i].axvline(-2.5, -6,1, c='gray', ls='dashed', lw=0.5)
ax[i].axvline(2.5, -6,1, c='gray', ls='dashed', lw=0.5)
ax[1].legend(loc=1)
#ax22.legend(loc=4)
ax[0].set_ylabel('$G(x)$', fontsize=20, color='b')
ax12.set_ylabel('$P(x)$', fontsize=20, color='g')
ax12.set_yticklabels([])
ax[1].set_ylabel(r'$\phi_{A/B}(x)$', fontsize=20)
ax22.set_ylabel(r'$\phi_{A}(x)\phi_{B}(x)$', fontsize=20, color='r')
ax[2].set_ylabel(r'$P(x)\phi_{A}(x)\phi_{B}(x)$', fontsize=20)
ax[2].set_yticklabels([])
plt.tight_layout(h_pad=0)
In [272]:
pTP = 2.*np.trapz(Px[i1:i2]*pfold[i1:i2]*(1-pfold[i1:i2]), x[i1:i2])/np.trapz(Px,x)
print "pTP(Hummer) = %f"%pTP
tTP = np.trapz(Px[i1:i2]*pfold[i1:i2]*(1-pfold[i1:i2]), x[i1:i2])/np.trapz(Px,x)
In [277]:
h5file = "data/cossio_kl0_Dx1_Dq0.h5"
print h5file
f = h5py.File(h5file, 'r')
data = np.array(f['data'])
f.close()
fig, ax = plt.subplots(figsize=(12,2), sharex=True,sharey=True)
ax.plot(data[:,0],data[:,1],'.', markersize=1)
ax.set_ylim(-10,10)
ax.set_xlim(0,25e4)
ax.set_ylabel('x')
ax.set_xlabel('time')
plt.tight_layout(h_pad=0)
fig, ax = plt.subplots(figsize=(6,4))
hist, bin_edges = np.histogram(data[:,1], bins=np.linspace(-10,10,50), normed=True)
bin_centers = [0.5*(bin_edges[i]+bin_edges[i+1]) \
for i in range(len(bin_edges)-1)]
ax.plot(bin_centers, -np.log(hist))
ax.set_xlim(-7,7)
ax.set_ylim(0,10)
ax.set_xlabel('x')
ax.set_ylabel('F ($k_BT$)')
tau_f, tau_u, data_f, data_u, tp_f, tp_u, recrossings = \
kinetics.lifetimes(data, f_bound=-2.5, u_bound=2.5)
tpt_u = []
tpt_f = []
for td in tp_u:
try:
tpt_u.append(td[:,0][-1] - td[:,0][0])
except IndexError:
pass
for td in tp_f:
try:
tpt_f.append(td[:,0][-1] - td[:,0][0])
except IndexError:
pass
fig, ax = plt.subplots(2,1,figsize=(6,4))
histogram, bin_edges, patches = \
ax[0].hist(np.log(tau_f + tau_u), \
bins=np.linspace(0,12,20), \
normed=True, cumulative=False, alpha=0.75)
param = expon.fit([tau_f + tau_u], floc=0)
tau_fold = param[1]
tau_error = param[1]/np.sqrt(len(tau_f + tau_u))
lx = np.linspace(0,12,200)
pdf = 1/param[1]*np.exp(lx)*np.exp(-1/param[1]*(np.exp(lx)))
ax[0].plot(lx, pdf, lw=3, color='b')
histogram, bin_edges, patches = \
ax[1].hist([x for x in (tpt_f + tpt_u) ], \
bins=np.linspace(0,14,15), \
normed=True, cumulative=False, alpha=0.75)
tpt = np.mean([x for x in (tpt_f + tpt_u) if x >= 0])
tpt_error = tpt/np.sqrt(len([x/2. for x in (tpt_f + tpt_u) if x >= 0]))
ax[0].set_yticklabels([])
ax[1].set_yticklabels([])
ax[0].set_xlabel(r'ln($\tau$)')
ax[1].set_xlabel(r'$\tau_\mathrm{TP}$')
ax[0].set_ylabel(r'P(ln$\tau$)')
ax[1].set_ylabel(r'P($\tau_\mathrm{TP}$)')
plt.tight_layout(h_pad=0.5, w_pad=0)
In [274]:
print "Folding time: %4.1f +/- %4.1f"%(tau_fold, tau_error)
print "Transition path time: %4.1f +/- %4.1f"%(tpt, tpt_error)
print "Barrier Estimate : %4.1f +/-%4.1f kBT"%(-np.log(tpt/tau_fold), -np.log(tpt/(tau_fold+tau_error)) + np.log(tpt/tau_fold) )
In [275]:
pTP = (np.sum(tpt_f) + np.sum(tpt_u))/float(len(data))
print "pTP (direct) = %4.2e"%(pTP)
In [323]:
from scipy.optimize import curve_fit
def parabola(x, a, b, c):
return a*(x-b)**2 + c
fig, ax = plt.subplots()
x = np.linspace(-10,10,1000)
Gx = [cossio.Gx(y, barrier=5.) for y in x]
ax.plot(x, Gx)
ax.set_ylabel('$G(x)$', fontsize=20)
ax.set_ylim(-1.1*5,0.5*5)
ax.set_xlabel('x')
popt_f, pcov = curve_fit(parabola, x[0:400], Gx[0:400])
print " kappa_f: %g kBT/nm**2"%popt_f[0]
popt_b, pcov = curve_fit(parabola, x[400:600], Gx[400:600])
print " kappa_b: %g kBT/nm**2"%popt_b[0]
popt_u, pcov = curve_fit(parabola, x[600:], Gx[600:])
#print popt_u
print " kappa_u: %g kBT/nm**2"%popt_u[0]
ax.plot(np.linspace(3,5,50),parabola(np.linspace(3,5,50), popt_u[0], popt_u[1], popt_u[2]),\
'hotpink', lw=3)
ax.plot(np.linspace(-5,-3,50),parabola(np.linspace(-5,-3,50), popt_f[0], popt_f[1], popt_f[2]), \
'lightgreen', lw=3)
ax.plot(np.linspace(-1,1,50),parabola(np.linspace(-1,1,50), popt_b[0], popt_b[1], popt_b[2]), \
'cyan', lw=3)
In [342]:
from scipy.special import erf
from sympy import *
fig, ax = plt.subplots()
histogram, bin_edges, patches = \
ax.hist([x for x in (tpt_f + tpt_u) ], \
bins=np.linspace(0,14,15), \
normed=True, cumulative=False, alpha=0.75)
tpt = np.mean([x for x in (tpt_f + tpt_u) ])
tpt_error = tpt/np.sqrt(len([x/2. for x in (tpt_f + tpt_u) if x >= 0]))
t=np.linspace(0.1,14,100)
ax.plot(t,2*np.sqrt(5)/(1-erf(5))*np.exp(-5*mpmath.coth(popt_f[0]*t)), lw=4, ls='dashed')
print mpmath.coth(popt_f[0]*1)
ax.plot(t,2*5*1/0.625*np.exp(-1/0.625*t), lw=4, ls='dashed')
ax.set_ylim(0,1.5*np.max(histogram))
ax.set_xlim(0,8)
In [ ]: