<jupyter_start><jupyter_text>[View in Colaboratory](https://colab.research.google.com/github/ktaletsk/CPF/blob/master/1D_Example/CPF_1D_toy.ipynb)  Coarse-graining of Brownian dynamics of toy modelPredicting the first passage time using the constant probability flux (CPF) method. This is Jupyter notebook written in Python language. You can run each cell by pressing CTRL+SHIFT on your keyboard. Reproduce results pulished in paper and feel free to make changes and use this code in your research  Define variables and functions and generate simulation code<jupyter_code>from sympy import *
from sympy.utilities.codegen import codegen
import numpy as np
from IPython.display import Markdown
import matplotlib.pyplot as plt
# import scipy.integrate as integrate
from IPython.display import HTML, Math
def latex_print(expr):
  display(HTML("<script src='https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/" "latest.js?config=default'></script>"))
  return Math(latex(expr))<jupyter_output><empty_output><jupyter_text>Consider a detailed dynamics, say MD, with the trajectory of some interesting species can reasonably be mapped to dynamics along a single reaction coordinate $x$.It might then be desirable to coarse grain these simulations to a nonlinear, stochastic differential equation.It is well established that such dynamics are equivalently described by a Fokker-Planck (or Smoluchowski) equation\begin{equation}\frac{\partial p}{\partial t} =    \frac{\partial}{\partial x}\frac{1}{\zeta} \left[ \left(\frac{dF^\text{CG}}{d x}\right) p    + k_\text{B}T \frac{\partial p}{\partial x}    \right],\end{equation}where $p\equiv p(x;t|x_0;t_0)$ is the probability density of finding the species in the neighborhood of $x$ at time $t$, given that it had coordinate $x_0$ at some earlier time $t_0$.The equivalent Ito-form stochastic differential equation is\begin{equation}dx = \left( \frac{d \zeta^{-1}}{d x} - \zeta^{-1} \frac{d F}{dx} \right) dt + \sqrt{2 \zeta^{-1}} d W_t\end{equation}We chose free energy of our toy system in a form:\begin{equation}\frac{F}{k_B T} = H_0 (1-x_0)^2 + \sum_{i=1}^{n-1} H_i (x_i - \mu_i)^2\end{equation}where $\mu_i=\frac{\sin A(x_0-C_i)}{B_i}$<jupyter_code>#Number of dimensions
n=1
A = []
B = []
C = []
H = [4]
#Define variables
x = symbols('x0:%d'%n)

def floor(x,A,B,C):
    return sin(A*(x-C))/B

def wall(x,y,A,B,C):
    return (y - floor(x,A,B,C))**2

def free_energy(n, x):
    result = H[0]*(1-x[0])**2
    for i in range(1, n):
        result += H[i] * wall(x[0], x[i], A[i-1], B[i-1], C[i-1])
    return result

def free_energy_CG(x):
    return H[0]*(1-x[0])**2

def zeta_inv(n, x):
    return 0.5*(1+x[0]**2)**2 #friction

def tshift(n,x,i):
    return diff(zeta_inv(n,x), x[i])-zeta_inv(n, x)*diff(free_energy(n,x), x[i])

def wshift(n,x):
    return sqrt(2.0*zeta_inv(n, x))<jupyter_output><empty_output><jupyter_text>**Free energy of our one-dimensional toy model**<jupyter_code>latex_print(free_energy(n, x))
xx = MatrixSymbol('x', 1, n)
dx = MatrixSymbol('dx', 1, n)
dWt = MatrixSymbol('dWt', 1, n)
dt = symbols('dt')
subs_rule = list(zip(x,xx))
array = [(tshift(n,x,i) * dt + wshift(n,x) * dWt[i]).subs(subs_rule) for i in range(n)]
m = Matrix(1, n, array)
equation = Eq(dx, m)
latex_print(equation)<jupyter_output><empty_output><jupyter_text>**Calculate equilibrium distribution** \begin{equation}p(x_i | x_0) = \frac{p(x_i, x_0)}{p(x_0)}\end{equation}<jupyter_code>p = symbols('p0:%d'%(n-1))
pp = MatrixSymbol('p', 1, n)
xx = MatrixSymbol('x', 1, n)
xa = symbols('xa')

subs_rule2 = list(zip(p,pp))
subs_rule2.insert(0,(x[0],xa))

array2 = [(floor(x[0],A[i],B[i],C[i]) + p[i]/sqrt(2*H[i+1])).subs(subs_rule2)  for i in range(0,n-1)]
array2.insert(0,xa)
m2 = Matrix(1, n, array2)
equation2 = Eq(xx, m2)
latex_print(equation2)<jupyter_output><empty_output><jupyter_text>Generate code for timestep and initialization<jupyter_code>#Put source files 'FPT-gen.cpp' and 'CPF-gen.cpp' in the same folder as this notebook (for Google Colab)
!git clone https://github.com/ktaletsk/CPF.git; cp CPF/1D_Example/CPF-gen.cpp .; cp CPF/1D_Example/FPT-gen.cpp .; rm -rf CPF
!ls
#Compile Ito equation source files for given free energy and friction
codegen(("shift",equation), "C", "shift", header=False, empty=False, to_files=True)
codegen(("init",equation2), "C", "init", header=False, empty=False, to_files=True)
#Complie CPF code
! g++ FPT-gen.cpp shift.c init.c -o FPT-gen
! g++ CPF-gen.cpp shift.c init.c -o CPF-gen
#Run CPF simulation
!./CPF-gen<jupyter_output><empty_output><jupyter_text>Processing CPF results<jupyter_code>with open('cpf.dat') as f:
    lines = f.readlines()
    x0 = np.array([float(line.split()[0]) for line in lines])
    x0.sort()
x0.shape[0]
P = np.arange(0.0, 1.0, 1.0/x0.shape[0])
fig = plt.figure(figsize=(12, 9))
ax1 = fig.add_subplot(111)
ax1.set_title("Steady state probability density distribution")
ax1.set_xlabel(r'$x$')
ax1.set_ylabel(r'$\mathbb{P}^{ss}$')
ax1.plot(x0,P, c='r', label=r'$f_d(t)$')
plt.axvline(x=0, color='k', linestyle='--')
plt.axvline(x=1,  color='k', linestyle='--')
plt.show()<jupyter_output><empty_output><jupyter_text>Using [Kernel Density Estimation](http://scikit-learn.org/stable/modules/density.html) to find PDF in steady state<jupyter_code>#Rule of thumb bandwith for kernel density estimstor
h = 1.06*np.std(x0)*np.size(x0)**(-0.2)
h
def myKernelDensity(x, h, xi):
    return np.sum(np.exp(-(x-xi)**2/(2*h**2)))/len(xi)
def myKernelDensityDerivative(x, h, xi):
    return np.sum(-(x-xi)/h**2*np.exp(-(x-xi)**2/(2*h**2)))/len(xi)
X_plot = np.linspace(0.0, 1.0, 1000)
mydens = [myKernelDensity(i, h, x0) for i in X_plot]
mydensder = [myKernelDensityDerivative(i, h, x0) for i in X_plot]
fig2 = plt.figure(figsize=(12, 9))
ax2 = fig2.add_subplot(111)
ax2.set_title("Steady state probability density distribution")
ax2.set_xlabel(r'$x$')
ax2.set_ylabel(r'$p^{ss}$')
ax2.plot(X_plot, np.array(mydens)/np.sqrt(2*np.pi*h**2), c='r')
plt.axvline(x=0, color='k', linestyle='--')
plt.axvline(x=1,  color='k', linestyle='--')
plt.show()
pCG = lambdify(x[0], exp(-free_energy_CG(x))/simplify(integrate(exp(-free_energy_CG(x)),(x[0],-oo,oo))), "sympy")
PCG = lambdify(x[0], simplify(integrate(exp(-free_energy_CG(x)),(x[0],-oo,x[0])))/simplify(integrate(exp(-free_energy_CG(x)),(x[0],-oo,oo))), "sympy")
with open('flux.dat') as f:
    lines = f.readlines()
J0 = float(lines[-1].split()[1])/float(lines[-1].split()[0])
xr=1
arr1 = np.array([(N(PCG(xr))-N(PCG(xT)))/N(pCG(xT))/J0 for xT in X_plot], dtype="float")
N(PCG(2))<jupyter_output><empty_output><jupyter_text>Friction  1. Friction from KDE of $p^{ss}$<jupyter_code>dFdx = lambdify(x[0],diff(free_energy_CG(x)))
actual_friction = lambdify(x[0],zeta_inv(n, x))
dFdX_arr = [dFdx(i) for i in X_plot]
friction = -(np.multiply(mydens,dFdX_arr)+mydensder)/J0/np.sqrt(2*np.pi*h**2)
act_f = [1/actual_friction(i) for i in X_plot]
fig5 = plt.figure(figsize=(12, 9))

ax5 = fig5.add_subplot(111)

ax5.set_title("Friction")
ax5.set_xlabel(r'$x$')
ax5.set_ylabel(r'$\zeta(x)/k_B T$')

ax5.plot(X_plot, friction, c='k',label='KDE estimate')
ax5.plot(X_plot, act_f, c='g', label='Actual friction')
plt.axvline(x=h, color='r', linestyle='--')
plt.axvline(x=xr-h, color='r', linestyle='--')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>2. Fitting friction with a sum of Chebyshev polynomials \begin{equation}\mathbb{P}^\text{ss}(x) =\int_{x_s}^{x} \mathbb{P}^\text{eq}(x') \frac{\tilde{\zeta}(x')}{p^\text{eq}(x')} dx' + \mathbb{P}^\text{eq}(x) \int_{x}^{x_a} \frac{\tilde{\zeta}(x')}{p^\text{eq}(x')} dx'\end{equation}where \begin{equation}\tilde{\zeta}(x) := \frac{J_0 \zeta(x)}{k_\text{B}T}\end{equation}We fit friction with the sum of Chebychev's polynomials, so that\begin{equation}\tilde{\zeta}(x) = \sum_i c_i T_i(x)\end{equation}\begin{equation}\mathbb{P}^\text{ss}(x) = \sum_i c_i \left( \int_{x_s}^{x} \mathbb{P}^\text{eq}(x') \frac{T_i(x')}{p^\text{eq}(x')} dx' + \mathbb{P}^\text{eq}(x) \int_{x}^{x_a} \frac{T_i(x')}{p^\text{eq}(x')} dx' \right) = \sum_i c_i I_i(x)\end{equation}We define the set of reference points $\{ x_j \}$ and solve matrix equation\begin{equation}c_i I_i(x_j) = \mathbb{P}^{ss}(x_j)\end{equation}to find the set of weights $\{c_i\}$<jupyter_code>import scipy.integrate
def integral1(n, lower, upper):
    return scipy.integrate.quad(lambda x: PCG(x)*chebyshevt_poly(n,x)/pCG(x), lower, upper)[0]

def integral2(n, lower, upper):
    return scipy.integrate.quad(lambda x: chebyshevt_poly(n,x)/pCG(x), lower, upper)[0]

def Pss_cheb(x, n, L0, La):
    return integral1(n,L0,x)+N(PCG(x))*integral2(n,x,La)
  
# def PSS_func(x):
  
from scipy.interpolate import interp1d
PSS_func = interp1d(x0, P)
PSS_func(0.5)
nmodes=7 #number of Chebyshev polynomials (starting from 0-th)
npoints=500 #reference points in domain
III = [[Pss_cheb(float(j)/npoints, i, 0, 1) for i in range(nmodes)] for j in range(npoints)]
rhs = [PSS_func(float(j)/npoints) for j in range(npoints)]
fit = np.linalg.lstsq(III, rhs, rcond=-1)[0]
fit
fig = plt.figure(figsize=(12, 9))

ax1 = fig.add_subplot(111)

ax1.set_title("Steady state probability density distribution (Chebychev fit)")
# ax1.set_xlabel(r'$t/\tau_c$')
# ax1.set_ylabel(r'$f_d(t)$')

ax1.plot(x0,P, c='r', label=r'$f_d(t)$')
ax1.plot(np.arange(float(npoints))/npoints, [sum([fit[i]*Pss_cheb(float(j)/npoints, i, 0, 1) for i in range(nmodes)]) for j in range(npoints)], c='k', label='Chebychev approximation')


#plt.xlim(xmin=-1.1, xmax=-0.9)
#plt.ylim(ymin=0, ymax=0.05)
# leg = ax1.legend()
#ax1.set_xscale('log')
# ax1.set_yscale('log')

plt.show()
cheb_friction = [sum([fit[i]*chebyshevt_poly(i,j)/J0 for i in range(nmodes)]) for j in X_plot]
fig5 = plt.figure(figsize=(12, 9))

ax5 = fig5.add_subplot(111)

ax5.set_title("Friction")
ax5.set_xlabel(r'$x_0$')
ax5.set_ylabel(r'$\zeta(x_0)/k_B T$')

ax5.plot(X_plot, friction, c='k', label='KDE friction')
ax5.plot(X_plot, cheb_friction, c='r', label='Chebychev friction')
ax5.plot(X_plot, act_f, c='g', label='Actual friction')
leg = ax5.legend()

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19672293942403177, 'fraction_numerical': 0.025124131082423038, 'mean_word_length': 5.931176875430144, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 39, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '453305', 'license': 'MIT', 'n_tokens_mistral': 4341, 'n_tokens_neox': 4203, 'n_words': 855, 'path': '1D_Example/CPF_1D_toy.ipynb', 'repo_name': 'ktaletsk/CPF', 'size': '10070'}