In [7]:
import numpy as np
import h5py
import scipy
import scipy.stats
from scipy.spatial.distance import pdist, squareform
import os
import sys
import errno
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
we_dir = os.path.abspath('we_custom')
bf_dir = os.path.abspath('bruteforce')
if bf_dir not in sys.path:
sys.path.append(bf_dir)
from calculate_statistics import calc_stats as bf_calc_stats
In [8]:
import seaborn as sns
sns.set(context="notebook", style="ticks", font="Arial")
mpl.rcParams.update({"xtick.major.width": 0.5, "ytick.major.width": 0.5, "savefig.dpi": 150,
"xtick.labelsize": 10, "ytick.labelsize": 10, "axes.labelsize": 11})
#mpl.rcParams.update({"xtick.direction": "in",
# "ytick.direction": "in",})
double_width = 7 # inches
single_width = 3.25 # inches
sc = sns.color_palette()
In [9]:
def mkdir_p(path):
try:
os.makedirs(path)
except OSError as exc: # Python >2.5
if exc.errno == errno.EEXIST and os.path.isdir(path):
pass
else:
raise
def save_figure(fig, figname, suppl=False, bbox_extra_artists=None):
fig_dir = 'supplemental_figures' if suppl else 'figures'
mkdir_p(fig_dir)
fig.savefig("%s/%s.pdf" % (fig_dir, figname), dpi=300, bbox_extra_artists=bbox_extra_artists, bbox_inches='tight')
fig.savefig("%s/%s.eps" % (fig_dir, figname), dpi=300, bbox_extra_artists=bbox_extra_artists, bbox_inches='tight')
fig.savefig("%s/%s.tiff" % (fig_dir, figname), dpi=300, bbox_extra_artists=bbox_extra_artists, bbox_inches='tight')
In [10]:
bf_data = np.load(os.path.join(bf_dir, 'data', 'md-solvent-langevin-distance.npy'))
n_a2b, n_b2a, time_in_a, time_in_b, tt_a2b, tt_b2a = bf_calc_stats(bf_data)
total_time = time_in_a + time_in_b
frac_a = (1.0 * time_in_a) / total_time
frac_b = (1.0 * time_in_b) / total_time
flux_a2b = (1.0 * n_a2b) / total_time
flux_b2a = (1.0 * n_b2a) / total_time
mfpt_a2b = frac_a / flux_a2b
mfpt_b2a = frac_b / flux_b2a
rate_a2b = flux_a2b / frac_a
rate_b2a = flux_b2a / frac_b
# Get confidence itervals for the rates
def mean_confidence_interval(data, confidence=0.95):
a = 1.0*np.array(data)
n = len(a)
m, se = np.mean(a), scipy.stats.sem(a)
h = se * scipy.stats.t.ppf((1+confidence)/2., n-1)
return m, m-h, m+h
def calc_ci(data, ngrps, confidence=0.95):
ds = np.split(data, N)
rates_a2b_splits = []
rates_b2a_splits = []
for d in ds:
n_a2b, n_b2a, time_in_a, time_in_b, tt_a2b, tt_b2a = bf_calc_stats(d)
total_time = time_in_a + time_in_b
flux_a2b = (1.0 * n_a2b) / total_time
flux_b2a = (1.0 * n_b2a) / total_time
rate_a2b = flux_a2b / frac_a
rate_b2a = flux_b2a / frac_b
rates_a2b_splits.append(rate_a2b)
rates_b2a_splits.append(rate_b2a)
return mean_confidence_interval(rates_a2b_splits, confidence), mean_confidence_interval(rates_b2a_splits, confidence)
N = 5
r_a2b_ci, r_b2a_ci = calc_ci(bf_data, 10)
print rate_a2b
print r_a2b_ci
print rate_b2a
print r_b2a_ci
In [11]:
def plot_we_direct_rate(ax, h5file, istate, jstate, color, sf=1.0):
with h5py.File(h5file, 'r') as f:
data = f['rate_evolution'][:,istate, jstate]
iters = data['iter_stop'] - 1
ax.plot(iters, data['expected']/sf, color=color)
ax.fill_between(iters, data['ci_lbound']/sf, data['ci_ubound']/sf,
color=color, facecolor=color, alpha=0.5)
def plot_bf_reference(ax, value, max_time, sf=1.0):
t = np.linspace(0, max_time, 20)
ax.plot(t, (value[0]/sf)*np.ones_like(t), color='black', lw=1)
ax.plot(t, (value[1]/sf)*np.ones_like(t), '--', color='black', lw=1)
ax.plot(t, (value[2]/sf)*np.ones_like(t), '--', color='black', lw=1)
In [12]:
def rate_evolution_plot():
fig = plt.figure(figsize=(single_width, 5))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212, sharex=ax1)
plt.setp( ax1.get_xticklabels(), visible=False)
kinavg_file = os.path.join(we_dir, 'kinavg.h5')
sf = 1.0E-3
plot_we_direct_rate(ax1, kinavg_file, 0, 1, sc[0], sf)
plot_we_direct_rate(ax2, kinavg_file, 1, 0, sc[0], sf)
plot_bf_reference(ax1, r_a2b_ci, 4000, sf)
plot_bf_reference(ax2, r_b2a_ci, 4000, sf )
ax1.set_ylim(0, 0.0025/sf)
ax2.set_ylim(0, 0.0025/sf)
ax2.set_xticks([0, 1000, 2000, 3000, 4000])
ax1.set_ylabel(r'$k_{A{\to}B}{\times}10^3\, \left(\tau^{-1}\right)$')
ax2.set_ylabel(r'$k_{B{\to}A}{\times}10^3\, \left(\tau^{-1}\right)$')
ax2.set_xlabel('Iteration')
plt.tight_layout()
save_figure(fig, 'wca_rate_convergence')
rate_evolution_plot()
In [13]:
def downsample_points(x, y, cutoff):
'''Downsample points so that adjacent points are at least cutoff apart'''
ii = [0,]
indx = 0
dist = squareform(pdist(np.column_stack([x,y])))
while indx <= x.shape[0]:
for k in xrange(indx + 1, dist.shape[1]):
if dist[indx, k] > cutoff:
ii.append(k)
indx = k
break
else:
indx = x.shape[0] + 1
return ii
def prob_dist_plot():
scale = 'log'
fig = plt.figure(figsize=(single_width, 3))
ax = fig.add_subplot(111)
r0 = 2.**(1./6.) * 3.4
# WE data
start_iter = 1000
with h5py.File(os.path.join(we_dir, 'pdist.h5'), 'r') as f:
h_we = f['histograms'][start_iter:].sum(axis=0)
bin_edges_we = f['binbounds_0'][:]
bin_midpts_we = f['midpoints_0'][:]
# Brute force data
h_bf, bin_edges_bf = np.histogram(bf_data, bins=bin_edges_we, density=False)
bin_midpts_bf = bin_edges_bf[0:-1] + 0.5*(bin_edges_bf[1] - bin_edges_bf[0])
ax.plot(bin_midpts_bf/r0, (1.0*h_bf)/h_bf.sum(), lw=1)
x = bin_midpts_we/r0
y = h_we/h_we.sum()
ax.plot(x, y, 'o', mfc='none', mec=sc[0], mew=1, ms=3)
ax.set_yscale(scale)
if scale == 'log':
ax.set_ylim(1E-6, 0.1)
else:
ax.set_ylim(0, 0.1)
ax.set_xlabel(r'$r/r_0$')
ax.set_ylabel(r'$p\left(r/r_0\right)$')
plt.tight_layout()
save_figure(fig, 'wca_prob_dist')
prob_dist_plot()
In [15]:
def prob_dist_plot2():
scale = 'linear'
fig = plt.figure(figsize=(single_width, 3))
ax = fig.add_subplot(111)
r0 = 2.**(1./6.) * 3.4
# WE data
start_iter = 1000
stop_iter = 4000
with h5py.File(os.path.join(we_dir, 'pdist.h5'), 'r') as f:
h_we = f['histograms'][start_iter:stop_iter].sum(axis=0)
bin_edges_we = f['binbounds_0'][:]
bin_midpts_we = f['midpoints_0'][:]
# Brute force data
h_bf, bin_edges_bf = np.histogram(bf_data, bins=bin_edges_we, density=False)
bin_midpts_bf = bin_edges_bf[0:-1] + 0.5*(bin_edges_bf[1] - bin_edges_bf[0])
ax.plot(bin_midpts_bf/r0, (1.0*h_bf)/h_bf.sum(), lw=1)
x = bin_midpts_we/r0
y = h_we/h_we.sum()
ax.plot(x, y, 'o', mfc='none', mec=sc[0], mew=1, ms=3, lw=1)
ax.set_yscale(scale)
if scale == 'log':
ax.set_ylim(1E-6, .1)
else:
ax.set_ylim(0, 0.1)
ax.set_xlabel(r'$r/r_0$')
ax.set_ylabel(r'$p\left(r/r_0\right)$')
plt.tight_layout()
save_figure(fig, 'wca_prob_dist2')
prob_dist_plot2()
In [16]:
from matplotlib.patches import Ellipse
def plot_potential():
plot_bins = False
fig = plt.figure(figsize=(single_width, 2))
ax = fig.add_subplot(111)
h = 5.0
sigma = 3.4
r0 = 2**(1./6) * sigma
s = 0.5 * r0
r = np.linspace(0,10,200)
V = h * (1 - (r - r0 - s)**2/s**2)**2
ax.plot(r/r0, V)
# Bins
if plot_bins:
bins = np.arange(3.3, 8.1, 0.1) / r0
y = np.linspace(0,10,10)
for b in bins:
ax.plot(b*np.ones_like(y), y, 'gray', lw=0.5)
# Draw cartoon
# dimer
rad = 100
y = 0.7
x1 = 0.25
x2 = x1 + 0.08
x3 = 0.60
x4 = x3 + 2*(x2 - x1)
pts = np.array([[x1, y],
[x2, y],
[x3, y],
[x4, y]])
ax.scatter(pts[:,0], pts[:,1], s=rad, edgecolors='black', facecolors='white', lw=1, transform=ax.transAxes)
ax.plot([x1, x2], [y, y], '-', lw=2, c='black', transform=ax.transAxes, zorder=0)
ax.plot([x3, x4], [y, y], '-', lw=2, c='black', transform=ax.transAxes, zorder=0)
#solvent
x0a = x1 + 0.5*(x2 - x1)
x0b = x3 + 0.5*(x4 - x3)
pts = np.array([[x0a + 0.01, y + 0.17],
[x0a - 0.06, y + 0.18],
[x0a + 0.09, y + 0.11],
[x0a + 0.03, y - 0.16],
[x0a - 0.04, y - 0.167],
[x0a - 0.1, y - 0.09],
[x0a - 0.115, y + 0.063],
[x0a + 0.1, y - 0.09]])
ax.scatter(pts[:,0], pts[:,1], s=rad, edgecolors='black', facecolors='gray', lw=1, transform=ax.transAxes, alpha=0.8)
pts = np.array([[x0b + 0.03, y + 0.11],
[x0b - 0.045, y + 0.16],
[x0b + 0.11, y + 0.138],
[x0b + 0.02, y - 0.12],
[x0b - 0.05, y - 0.163],
[x0b - 0.12, y - 0.125],
[x0b - 0.15, y + 0.02],
[x0b - 0.115, y + 0.15],
[x0b + 0.09, y - 0.15],
[x0b + 0.145, y - 0.03],])
#[x0b + 0.1, y - 0.1]])
ax.scatter(pts[:,0], pts[:,1], s=rad, edgecolors='black', facecolors='gray', lw=1, transform=ax.transAxes, alpha=0.8)
ax.set_xticks([1,2])
ax.set_ylim(0,12)
ax.set_xlim(0.5, 2.5)
ax.set_xlabel(r'$r/r_0$')
ax.set_ylabel(r'$U_{bond}\left(r\right) / k_BT$')
sns.despine()
fig.subplots_adjust(hspace=0.01, bottom=0.05)
plt.tight_layout()
save_figure(fig, 'potential')
plot_potential()
In [17]:
from matplotlib.patches import Ellipse
def plot_potential():
plot_bins = False
fig = plt.figure(figsize=(single_width, 2))
ax = fig.add_subplot(111)
h = 5.0
sigma = 3.4
r0 = 2**(1./6) * sigma
s = 0.5 * r0
r = np.linspace(0,10,200)
V = h * (1 - (r - r0 - s)**2/s**2)**2
ax.plot(r/r0, V)
# Bins
if plot_bins:
bins = np.arange(3.3, 8.1, 0.1) / r0
y = np.linspace(0,10,10)
for b in bins:
ax.plot(b*np.ones_like(y), y, 'gray', lw=0.5)
# Draw cartoon
# dimer
rad = 100
y = 0.55
x1 = 0.22
x2 = x1 + 0.08
x3 = 0.64
x4 = x3 + 2*(x2 - x1)
pts = np.array([[x1, y],
[x2, y],
[x3, y],
[x4, y]])
ax.scatter(pts[:,0], pts[:,1], s=rad, edgecolors='black', facecolors='white', lw=1, transform=ax.transAxes)
ax.plot([x1, x2], [y, y], '-', lw=2, c='black', transform=ax.transAxes, zorder=0)
ax.plot([x3, x4], [y, y], '-', lw=2, c='black', transform=ax.transAxes, zorder=0)
#solvent
x0a = x1 + 0.5*(x2 - x1)
x0b = x3 + 0.5*(x4 - x3)
pts = np.array([[x0a + 0.01, y + 0.17],
[x0a - 0.06, y + 0.18],
[x0a + 0.09, y + 0.11],
[x0a + 0.03, y - 0.16],
[x0a - 0.04, y - 0.167],
[x0a - 0.1, y - 0.09],
[x0a - 0.115, y + 0.063],
[x0a + 0.1, y - 0.09]])
#ax.scatter(pts[:,0], pts[:,1], s=rad, edgecolors='black', facecolors='gray', lw=1, transform=ax.transAxes, alpha=0.8)
pts = np.array([[x0b + 0.03, y + 0.11],
[x0b - 0.045, y + 0.16],
[x0b + 0.11, y + 0.138],
[x0b + 0.02, y - 0.12],
[x0b - 0.05, y - 0.163],
[x0b - 0.12, y - 0.125],
[x0b - 0.15, y + 0.02],
[x0b - 0.115, y + 0.15],
[x0b + 0.09, y - 0.15],
[x0b + 0.145, y - 0.03],])
#[x0b + 0.1, y - 0.1]])
#ax.scatter(pts[:,0], pts[:,1], s=rad, edgecolors='black', facecolors='gray', lw=1, transform=ax.transAxes, alpha=0.8)
ax.text(x0a - 0.1, y+0.12, 'Compact', transform=ax.transAxes)
ax.text(x0b - 0.12, y+0.12, 'Extended', transform=ax.transAxes)
ax.set_xticks([1,2])
ax.set_ylim(0,12)
ax.set_xlim(0.5, 2.5)
ax.set_xlabel(r'$r/r_0$')
ax.set_ylabel(r'$U_{bond}\left(r\right) / k_BT$')
sns.despine()
fig.subplots_adjust(hspace=0.01, bottom=0.05)
plt.tight_layout()
save_figure(fig, 'potential_nosolvent')
plot_potential()