In [1]:

import matplotlib.pyplot as plt
import matplotlib.cm as cm
%matplotlib inline
import h5py
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 [2]:

import numpy as np
import itertools
from scipy.stats import norm
import time
import cossio
import kinetics



### Molecular potential of mean force



In [3]:

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)

#ax[2].set_ylabel('$\partial^2 G(x)/\partial x^2$', fontsize=20)

plt.tight_layout()







In [4]:

x = np.linspace(-10,10,100)
q = np.linspace(-12,12,100)
kl_all = np.linspace(0,2,20)
G2d_all = []

for kl in kl_all:
G2d = np.ones((100, 100), float)*[cossio.Gx(y, barrier=5.) for y in x]
for i, j in itertools.product(range(100), range(100)):
G2d[i,j] += cossio.V(q[i], x[j], kl)
G2d = np.array(G2d)
G2d -= np.min(G2d)
G2d_all.append(G2d)




In [5]:

i = 0
for G2d in G2d_all:
fig, ax = plt.subplots(figsize=(6,5))
cs = ax.contourf(x, q,  G2d, cmap=cm.rainbow, levels=np.arange(0,10,1), alpha=0.9)
cbar = plt.colorbar(cs)
ax.text(-10,9.5,'$\kappa_l=%2.1f$'%kl_all[i], fontsize=20)
ax.set_xlim(-12,12)
ax.set_ylim(-12,12)
ax.set_xlabel('$x$', fontsize=20)
ax.set_ylabel('$q$', fontsize=20)
i+=1







In [6]:

barrier = []
i = 0
for G2d in G2d_all:
expGM = np.trapz(np.exp(-G2d), q, axis=0)
GM = -np.log(expGM)
expGA = np.trapz(np.exp(-G2d), x, axis=1)
GA = -np.log(expGA)
fig, ax = plt.subplots()
ax.plot(x, GM - np.min(GM), label='$G_M$', lw=4)
ax.plot(q, GA  - np.min(GA), label='$G_A$', lw=4)
barrier.append(np.max((GA  - np.min(GA))[45:55]))
#ax.plot(bin_centers, [cossio.Gx(y) for y in bin_centers], '--', c='red', lw=3)
ax.set_xlim(-10,10)
ax.set_ylim(-1,7)
ax.text(-8,6,'$\kappa_l=%2.1f$'%kl_all[i], fontsize=20)
ax.set_xlabel('Extension', fontsize=20)
ax.set_ylabel('Free Energy', fontsize=20)
ax.legend(loc=1)
i+=1
fig.tight_layout()