{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Metastability and Markov

\n",
"\n",
"(Sethna, \"Entropy, Order Parameters, and Complexity\", ex. 12.XXX)\n",
"\n",
"© 2019, James Sethna, all rights reserved."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The hints below will help you solve parts (e) and (f) of this question, where we will numerically evaluate the slowest decaying mode and barrier crossing time by computing eigenstates of the \"quantum\" Hamiltonian.\n",
"\n",
"In the Mathematica hints, we shall solve the differential equations directly. In the Python hints, we shall instead construct a Hamiltonian matrix by discretizing space into segments of length dx, and then finding the lowest eigenvalue of that Hamiltonian.\n",
"\n",
"Import packages"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%pylab inline\n",
"from scipy import *\n",
"from scipy.sparse import diags\n",
"from scipy.linalg import eig_banded\n",
"from scipy import special"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define functions to evaluate the potential $V(x)$ and effective quantum potential $V_\\rm{eff}(x)$ and compute a discretized Hamiltonian matrix $H$. We assume $\\eta=1$ and $k_B T = 1/2$. We will use this matrix to compute the slowest decaying mode for parts (e) and (f)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def V(x):\n",
" \"\"\"\n",
" Returns the cubic potential -x^3/3 + x\n",
" \"\"\"\n",
" return ...\n",
"\n",
"\n",
"def Veff(x, eta=1, kbT=1/2):\n",
" \"\"\"\n",
" Returns the effective quantum potential evaluated at x. \n",
" Calculate Veff by pluging the cubic potential into the \n",
" effective potential you found in part (g). You should find \n",
" (x^2 - 1)^2/(4 kbT eta) + x/eta (notice this is independent of eta)\n",
" \"\"\"\n",
" return ...\n",
"\n",
"\n",
"def hamiltonian(eta=1, kbT=1/2, x0=5, dx=0.001):\n",
" \"\"\"\n",
" Returns a discretized hamiltonian matrix for numerical eigendecomposition\n",
" +/-x0 - finite boundary conditions\n",
" dx - discrete grid size\n",
" \"\"\"\n",
" N = round(2*x0/dx) #total number of elements on grid (boundary at +/- x0)\n",
" x = arange(-x0,x0,dx) #make discrete grid\n",
" \n",
" ## Finite difference Kinetic Energy\n",
" #d^2/dx^2 is a matrix with -2 on the diagonal and 1 on the super/subdiagonal divided by dx^2\n",
" #In the first list put the elements: -2, 1, and 1. In the second list put the corresponding diagonal\n",
" #0 = diagonal, 1 = superdiagonal, -1 = subdiagonal\n",
" T = -kbT/eta*diags([..., ..., ...], [..., ..., ...], shape=(N, N))/(...)**2\n",
" \n",
" ## Potential Energy \n",
" #Effective quantum potential on discrete grid\n",
" Veff_x = Veff(x, eta, kbT)\n",
" #Matrix potential (Veff_x on diagonal)\n",
" Veff_mat = diags([...], [0])\n",
" \n",
" H = T + Veff_mat\n",
" \n",
" return H.toarray() \n",
"\n",
"def eigenSystem(eta=1, kbT=1/2, x0=5, dx=0.001, Nstates=1):\n",
" \"\"\"\n",
" Returns the Nstates lowest energy eigenvalues and associated eigenstates for \n",
" the quantum system corresponding to the Fokker-Plank equation with a cubic potential and noise g\n",
" \"\"\"\n",
" H = hamiltonian(eta, kbT, x0, dx)\n",
" \n",
" banded = array([diagonal(H), append(diagonal(H,1),0)])\n",
" eigs = eig_banded(banded, select=\"i\", select_range=[0,0], lower=True)\n",
" ##eigh_tridiagonal(diagonal(H), diagonal(H,1), select=\"i\", select_range=[0,0]) \n",
" # If you have scipy version 1.0 or newer, eigh_tridiagonal (add from scipy.linalg import eigh_tridiagonal)\n",
" # will give you slightly faster performance\n",
"\n",
" return eigs\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(e) * For the cubic potential (eqn 7),\n",
"numerically compute the eigenstate of \n",
"the transformed diffusion equation with smallest eigenvalue.\n",
"What does the eigenvalue predict for the lifetime?\n",
"How nearly does it agree with the classic calculation of Kramers:\n",
" \n",
"$\\lambda_0 \\approx \\sqrt{K {\\widetilde K}}/(2 \\pi \\eta) \\exp(-E/k_B T)$\n",
" *"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute eigenvalue and eigenstate\n",
"x0 = 5\n",
"dx=0.001\n",
"\n",
"val, vec = ...\n",
"val = val[0]\n",
"vec = transpose(vec)[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot the potential, notice where the well and barrier are located\n",
"x = arange(...)\n",
"\n",
"plot(x, ...)\n",
"xlabel(\"x\",size=20)\n",
"ylabel(\"V(x)\",size=20)\n",
"ylim([-10,10])\n",
"\n",
"show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compare eigenvalue approximation for escape time to analytical result\n",
"\n",
"tauNumerical = 1/(...)\n",
"tauKramers = 1/(sqrt(... * ...)/(2 * pi * ...) * exp(-(...)/(...)))\n",
"\n",
"# Percent difference between numerical and analytical result\n",
"diff = 100*abs(...)/tauNumerical\n",
"\n",
"print(tauNumerical)\n",
"print(tauKramers)\n",
"\n",
"print(str(round(diff,4))+\"% difference\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(e) * ... What does the eigenvalue predict for the lifetime? How nearly does it agree with $\\tau$ from the Kramers calculation? **\n",
"\n",
"(f) ** Using the corresponding eigenstate $\\rho_0$, plot the slowest decaying mode $\\rho_0(x) = \\sqrt{\\rho^*} \\sigma_0$, normalized to one, along with the Boltzmann distribution $\\rho^*(x)/Z$ and the Boltzmann probability distribution in the approximation that the well is quadratic. Explain how $\\rho_0$ differs from the quadratic approximation, and why this is physically sensible. Explain how $\\rho_0$ differs from $\\rho^*$, and how it resolves an important question about how to determine the metastable probability `in the well'.*"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute slowest decaying mode\n",
"\n",
"# Boltzmann distribution for cubic potential\n",
"boltz = array(exp(...))\n",
"\n",
"# Compute slowest decaying mode from eigenstate 'vec'\n",
"mode = (vec*sqrt(...))\n",
"\n",
"# Plot the unnormalized mode\n",
"plot(x,mode)\n",
"legend([r\"$\\rho_0$\"],prop={'size': 15})\n",
"xlabel(\"x\",size=20)\n",
"ylabel(\"Density\",size=20)\n",
"\n",
"show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Now normalize your slowest decaying mode\n",
"# We only know the slowest decaying mode on a discrete grid, so the norm is given by\n",
"# adding up the points\n",
"# Due to numerical errors your slowest decaying mode may blow up at one of the boundaries (the \n",
"# numerical eigenstate doesn't exactly cancel the blow up from the Boltzmann distribution)\n",
"# If this occurs, use cutoff to restrict the normalization calculation to [-xLim, xLim]\n",
"xLim = ...\n",
"# boltz, x, and mode are lists. 'cutoff' gives the range of list elements of interest\n",
"cutoff = round((x0-abs(xLim))/dx)+1\n",
"norm = sum(mode[cutoff:-cutoff])*dx\n",
"\n",
"# Adjust Z manually so that the Boltzmann distribution best matches the slowest decaying mode\n",
"# inside the well or approximate Z by normalizing over the inside of the potential well.\n",
"# If you choose the later option, you will want to restrict the normalization to the domain [xMin, xMax]\n",
"# Use your plot of the potential above to choose xMin and xMax\n",
"xMin = ...\n",
"xMax = ...\n",
"\n",
"# Pick range corresponding to inside of the well\n",
"cutoff1 = round((x0-abs(xMin))/dx)\n",
"cutoff2 = round((x0-abs(xMax))/dx)\n",
"\n",
"Z = sum(boltz[...:-...])*(...)\n",
"\n",
"\n",
"# Compute the Boltzmann distribution in the approximation that the well is quadratic\n",
"boltzQuadratic = exp(...)\n",
"ZQuadratic = sum(...)*(...)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compare Boltzmann distributions to mode\n",
"figure(1,[9,7])\n",
"plot(x[cutoff:-cutoff],V(x[cutoff:-cutoff]))\n",
"plot(x[cutoff:-cutoff],boltz[cutoff:-cutoff]/Z)\n",
"plot(x[cutoff:-cutoff],mode[cutoff:-cutoff]/norm)\n",
"plot(x[cutoff:-cutoff],boltzQuadratic[cutoff:-cutoff]/ZQuadratic,'--')\n",
"ylim([-1,1])\n",
"legend([r\"$V(x)$\",r\"$\\rho^*$\",r\"$\\rho_0$\",\"Quadratic Well\"],prop={'size': 15})\n",
"xlabel(\"x\",size=20)\n",
"ylabel(\"Density\",size=20)\n",
"\n",
"\n",
"show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}