{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Molecular Dynamics" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ ":::{figure-md} markdown-fig \n", "\n", "\"dp\"\n", "\n", "Simulation of Cu atoms \n", "\n", "::: " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ ":::{figure-md} markdown-fig\n", "\n", "![MD vs MC](./figs/mcmd.png){width=30%}\n", "\n", "**MD vs MC:** Both sample microstates. The former follows the natural motion (dynamics); the latter samples from the Boltzmann distribution using rules designed to improve efficiency.\n", ":::\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Timescales and Lengthscales\n", "\n", "- Classical Molecular Dynamics can access a hiearrchy of time-scales from pico seconds to microseconds. \n", "- It is also possible to go beyond the time scale of brute force MD byb emplying clever enhanced sampling techniques." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\n", ":::{figure-md} markdown-fig \n", "\n", "\"dp\"\n", "\n", "Different time-scales underlying different leng-scales/motions in molecules\n", "::: \n", "\n", "\n", ":::{figure-md} markdown-fig \n", "\n", "\"dp\"\n", "\n", "Simulation of water box\n", "\n", "::: \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Is MD just Newton's laws applied on big systems?\n", "\n", "**Not quite: Noble prize in Chemistry 2013**\n", "\n", "- Classical molecular dynamics (MD) is a powerful computational technique for studying complex molecular systems.\n", "- Applications span wide range including proteins, polymers, inorganic and organic materials. \n", "- Alos molecular dynamics simulation is being used in a complimentary way to the analysis of experimental data coming from NMR, IR, UV spectroscopies and elastic-scattering techniques, such as small angle scattering or diffraction.\n", "\n", "- [2013 Noble Lectures by M Karplus, A Warshell, M Levitt](https://www.youtube.com/watch?v=NuaeD9xYBtY)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrating equations of motion Numerically\n", "\n", "\n", "- The **Euler method** is the simplest numerical integrator for ordinary differential equations (ODEs).\n", "- More accurate integrators that include higher-order terms are known as **Runge-Kutta (RK)** methods — e.g., RK2, RK4, RK6.\n", "\n", "- Given an ODE in standard form:\n", "\n", "$$\n", "\\frac{d \\mathbf{y}}{dt} = \\mathbf{f}(t, \\mathbf{y}),\n", "$$\n", "\n", "- we approximate the derivative using finite differences:\n", "\n", "$$\n", "\\frac{d \\mathbf{y}(t)}{dt} \\approx \\frac{\\mathbf{y}_{n+1} - \\mathbf{y}_n}{\\Delta t},\n", "$$\n", "\n", "- leading to the **Euler update rule**:\n", "\n", "$$\n", "\\mathbf{y}_{n+1} \\approx \\mathbf{y}_n + \\Delta t \\cdot \\mathbf{f}(t_n, \\mathbf{y}_n).\n", "$$\n", "\n", "\n", "#### Example: Harmonic Oscillator\n", "\n", "The harmonic oscillator with parameters $m=1, k=1$ is defined by the state vector $\\mathbf{y} = \\begin{bmatrix} x \\\\ v \\end{bmatrix}$ with components for $\\dot{y}$ given by:\n", "\n", "$$\n", "\\dot{x} = v, \\quad \\dot{v} = -x.\n", "$$\n", "\n", "$$\n", "\\frac{d \\mathbf{y}}{dt} = \\mathbf{f}(\\mathbf{y}) = \\begin{bmatrix} v \\\\ -x \\end{bmatrix}.\n", "$$\n", "\n", "Using the Euler method, the system evolves as:\n", "\n", "$$\n", "\\mathbf{y}_{n+1} = \\mathbf{y}_n + \\Delta t \\cdot \\mathbf{f}(\\mathbf{y}_n).\n", "$$\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "y = np.array([1.0, 1.0]) # Initial [x, v]\n", "pos, vel = [], []\n", "t = 0\n", "dt = 0.1 \n", "\n", "for i in range(1000):\n", " \n", " dydt = np.array([y[1], -y[0]]) # [v, -x]\n", " y += dt * dydt # Euler update\n", " \n", " t += dt\n", " \n", " pos.append(y[0])\n", " vel.append(y[1])\n", "\n", "# Convert to arrays\n", "pos, vel = np.array(pos), np.array(vel)\n", "\n", "# Plot results\n", "fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(8,3))\n", "\n", "# Phase space\n", "ax[0].plot(pos, vel)\n", "ax[0].set_title(\"Phase space (x vs v)\")\n", "\n", "# Time series\n", "ax[1].plot(pos, label=\"x(t)\")\n", "ax[1].plot(vel, label=\"v(t)\")\n", "ax[1].legend()\n", "ax[1].set_title(\"Time evolution\")\n", "\n", "# Energy\n", "ax[2].plot(0.5*pos**2 + 0.5*vel**2)\n", "ax[2].set_title(\"Energy (should be constant)\")\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Verlet algortihm" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "Taylor expansion of position $\\vec{r}(t)$ after timestep $\\Delta t$ we obtain forward and backward Euler schems\n", "\n", "$$r_{t+\\Delta t} = r_t +v_t\\Delta t +\\frac{1}{2}a_t \\Delta t^2 + O(\\Delta t^3)$$\n", "\n", "$$r_{t-\\Delta t} = r_t -v_t \\Delta t +\\frac{1}{2}a_t \\Delta t^2 + O(\\Delta t^3)$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- In 1967 Loup Verlet introduced a new algorithm into molecular dynamics simulations which preserves energy is accurate and efficient.\n", "\n", "- Summing the two taylor expansion above we get a updating scheme which is an order of mangnitude more accurate\n", "\n", "$$r_{t+\\Delta t} = 2r_t - r_{t-\\Delta t} +a_t \\Delta t^2+O(\\Delta t^4)$$\n", "\n", "$$v_t = \\frac{r_{t+\\Delta t}-r_{t-\\Delta t}}{2\\Delta t} +O(\\Delta t^2) $$\n", "\n", "- Velocity is not needed to update the positions. But we still need them to set the temperature. \n", "\n", "- Terms of order $O(\\Delta t^3)$ cancel in position giving position an accuracy of order $O(\\Delta t^4)$\n", "\n", "- To update the position we need positions in the past at two different time points! This is is not very efficient." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Velocity Verlet updating scheme\n", "\n", "- A better updating scheme has been proposed known as **Velocity-Verlet (VV)** which stores **positions, velocities and accelerations at the same time.** Time stepping backward expansion $r(t-\\Delta t + \\Delta t)$ and summing with the forward Tayloer expansions we get Velocity Verlet updating scheme:\n", "\n", "$$v_{t+\\Delta t} = v_t + \\frac{1}{2}(a_t+a_{t+\\Delta t})\\Delta t +O(\\Delta t^3)$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- Substituting forces $a=\\frac{F}{m}$ instead of acelration we get \n", "\n", "$$r_{t+\\Delta t} = r_t + v_t\\Delta t + \\frac{F_t}{2m}\\Delta t^2$$\n", "\n", "\n", "$$v_{t+\\Delta t} = v_t + \\frac{F_t+F_{t+\\Delta t}}{2m}\\Delta t$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ ":::{admonition} **Velocity Verlet Algorithm** \n", ":class: tip\n", "\n", "**1. Evaluate the initial force from the current position:**\n", "\n", "$$\n", "F_t = -\\frac{\\partial U(r)}{\\partial r} \\Bigg|_{r(t)}\n", "$$\n", "\n", "**2. Update the position:**\n", "\n", "$$\n", "r_{t+\\Delta t} = r_t + v_t \\Delta t + \\frac{F_t}{2m} \\Delta t^2\n", "$$\n", "\n", "**3. Partially update the velocity:**\n", "\n", "$$\n", "v_{t+\\Delta t/2} = v_t + \\frac{F_t}{2m} \\Delta t\n", "$$\n", "\n", "**4. Evaluate the force at the new position:**\n", "\n", "$$\n", "F_{t+\\Delta t} = -\\frac{\\partial U(r)}{\\partial r} \\Bigg|_{r(t+\\Delta t)}\n", "$$\n", "\n", "**5. Complete the velocity update:**\n", "\n", "$$\n", "v_{t+\\Delta t} = v_{t+\\Delta t/2} + \\frac{F_{t+\\Delta t}}{2m} \\Delta t\n", "$$\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Molecular Dynamics of Classical Harmonic Oscillator (NVE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Velocity Verlet integration of harmonic oscillator" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "import numpy as np\n", "\n", "def run_md_nve_1d(x, v, dt, t_max, en_force):\n", " \"\"\"\n", " Minimalistic 1D Molecular Dynamics simulation (NVE ensemble) \n", " using Velocity Verlet integration.\n", " Simulates a particle moving in a 1D potential without thermal noise or friction\n", " (i.e., energy-conserving dynamics).\n", " \n", " Parameters\n", " ----------\n", " x : float\n", " Initial position.\n", " v : float\n", " Initial velocity.\n", " dt : float\n", " Time step for integration.\n", " t_max : float\n", " Total simulation time.\n", " en_force : callable\n", " Function that takes position `x` as input and returns a tuple (potential energy, force).\n", " \n", " Returns\n", " -------\n", " pos : ndarray\n", " Array of particle positions over time.\n", " vel : ndarray\n", " Array of particle velocities over time.\n", " KE : ndarray\n", " Array of kinetic energies over time.\n", " PE : ndarray\n", " Array of potential energies over time.\n", " \n", " Example\n", " -------\n", " >>> def harmonic_force(x):\n", " >>> k = 1.0\n", " >>> return 0.5 * k * x**2, -k * x\n", " >>> pos, vel, KE, PE = run_md_nve_1d(1.0, 0.0, 0.01, 10.0, harmonic_force)\n", " \"\"\"\n", " \n", " times, pos, vel, KE, PE = [], [], [], [], []\n", " \n", " # Initialize force and potential energy\n", " pe, F = en_force(x)\n", " \n", " t = 0.0\n", " for step in range(int(t_max / dt)):\n", " \n", " # Velocity Verlet Integration\n", " \n", " # Half-step velocity update\n", " v += 0.5 * F * dt\n", " \n", " # Full-step position update\n", " x += v * dt\n", " \n", " # Update force at new position\n", " pe, F = en_force(x)\n", " \n", " # Half-step velocity update\n", " v += 0.5 * F * dt\n", " \n", " # Save results\n", " times.append(t)\n", " pos.append(x)\n", " vel.append(v)\n", " KE.append(0.5 * v * v)\n", " PE.append(pe)\n", " \n", " # Advance time\n", " t += dt\n", " \n", " return np.array(pos), np.array(vel), np.array(KE), np.array(PE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Run NVE simulation of harmonic oscillator" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "#----parameters of simulation----\n", "k = 3 \n", "x0 = 1 \n", "v0 = 0\n", "dt = 0.01 * 2*np.pi/np.sqrt(k) #A good timestep determined by using oscillator frequency\n", "t_max = 1000\n", "\n", "### Define Potential Energy function\n", "def ho_en_force(x, k=k):\n", " '''Force field of harmonic oscillator:\n", " returns potential energy and force'''\n", " \n", " return k*x**2, -k*x\n", "\n", "### Run simulation\n", "pos, vel, KE, PE = run_md_nve_1d(x0, v0, dt, t_max, ho_en_force)\n", "\n", "# Plot results\n", "fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(8,3))\n", "\n", "# Phase space\n", "ax[0].plot(pos, vel)\n", "ax[0].set_title(\"Phase space (x vs v)\")\n", "\n", "# Time series\n", "ax[1].plot(pos, label=\"x(t)\")\n", "ax[1].plot(vel, label=\"v(t)\")\n", "ax[1].legend()\n", "ax[1].set_title(\"Time evolution\")\n", "\n", "# Energy\n", "ax[2].plot(0.5*pos**2 + 0.5*vel**2)\n", "ax[2].set_title(\"Energy (should be constant)\")\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot Distribution in phase-space" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "### Visualize\n", "fig, ax =plt.subplots(ncols=2)\n", "\n", "ax[0].hist(pos);\n", "ax[1].hist(vel, color='orange');\n", "ax[0].set_ylabel('P(x)')\n", "ax[1].set_ylabel('P(v)')\n", "\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Ensemble averages\n", "\n", ":::{admonition} **Ergodic hypothesis**\n", ":class: important\n", "\n", "$$\\langle A \\rangle = \\frac{1}{t} \\int^{ \\tau_{eq} +t}_{ \\tau_{eq} } A(t)$$\n", "\n", "- $\\tau_{eq}$ time it teks for system to \"settle into\" equilibrium and $t$ time of simulation. \n", "- Time averages are equal to ensemble averages when $t\\lim \\infty$. But in real life we need to determine how long to sample. \n", "- Ergodic assumption is: during the coruse of MD we visit all the microstates (or small but representative subset) that go into the ensemble average!\n", ":::" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ ":::{admonition} **Energy conservation in simulations**\n", ":class: tip\n", "\n", "- Kinetic, potential and total energies are flcutuation quantities in the Molecualr dynamics simulations.\n", "- Energy in (NVE) or its average $\\langle E \\rangle$ (in NVT, NPT, etc) must remain constant!\n", "\n", "$$KE_t = \\sum_i \\frac{p_t^2}{2m}$$\n", "\n", "$$PE_t = \\sum_{ij}u_{ij}(r_t)$$\n", "\n", "$$\\langle E \\rangle = \\langle KE \\rangle +\\langle PE \\rangle \\approx const $$\n", "\n", ":::\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ ":::{admonition} **Temperature control in simulations**\n", ":class: tip\n", "\n", "According to equipariting result of equilibrium statistical mechanics in the NVT ensmeble\n", "\n", "$$\\Big\\langle \\sum_i \\frac{p_t^2}{2m} \\Big\\rangle =\\frac{3 }{2}N k_B T$$\n", "\n", "\n", "$$T_t = \\frac{2}{3 N k_B} \\big\\langle KE \\big\\rangle_t $$\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ ":::{admonition} **Pressure control in simulations**\n", ":class: tip\n", "\n", "$$P = \\frac{1}{3V} \\sum^N_i \\Big [\\frac{p^2_i}{m} +\\vec{F}_i(r) \\cdot \\vec{r}_i(t) \\Big ]$$\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Langevin equation\n", "\n", "- A particle of mass $m$ moves under the force derived from a potential energy $U(x)$. The motion is purely deterministic.\n", "\n", "$$\n", "m \\ddot{x} = -\\nabla_x U\n", "$$\n", "\n", "- **Challenge:** What should we do when we have only a one or few particles, and cannot explicitly simulate the vast surrounding environment in order to assign temperature?\n", "- **Solution:** We model the surrounding medium (e.g., a solvent) as an **implicit thermal bath** that interacts with the particle.\n", " - The particle exchanges energy with the bath, maintaining **thermal equilibrium** at a fixed temperature $T$.\n", " - This motivates **Langevin dynamics**, where the effects of the solvent are captured by **friction** (dissipation) and **random thermal kicks** (fluctuations), without simulating solvent molecules explicitly.\n", "\n", ":::{admonition} **Langevin Equation**\n", ":class: important \n", "\n", "$$\n", "m \\ddot{x} = -\\nabla_x U - \\gamma \\dot{x} + \\eta(t)\n", "$$\n", "\n", ":::\n", "\n", ":::{admonition} **Overdamped Limit of Langevin Dynamics ($m \\ddot{x} = 0$)**\n", ":class: important \n", "\n", "$$\n", "\\gamma \\dot{x} = \\nabla_x U + \\eta(t)\n", "$$\n", "\n", ":::\n", "\n", "- The friction $\\gamma$ and thermal noise $\\eta(t)$ are clearly connected becasue the faster the particle movies (more noise) the more it also dissipates energy. \n", "- the connection is known by the name of Fluctuation-Dissipation Theorem:\n", "\n", ":::{admonition} **Fluctuation-Dissipation Theorem (FDT)**\n", ":class: important \n", "\n", "$$\n", "\\langle \\eta(t) \\eta(t') \\rangle = 2 \\gamma k_B T \\, \\delta(t - t')\n", "$$\n", "\n", ":::\n", "\n", "- The environment \"forgets\" what happened almost immediately after a collision very short memory hence noise terms are uncorrelated (independent). Just like what we had in brownian motion. \n", "- The FDT ensures that the strength of random thermal kicks is precisely tuned to the amount of viscous damping, so that the system reaches and maintains thermal equilibrium at temperature $T$\n", "- FDT connects diffusion (random spreading) and viscosity (resistance to motion), both fundamentally controlled by temperature.\n", "\n", ":::{admonition} **Derivation of FDT**\n", ":class: tip, dropdown\n", "\n", "\n", "We consider the **underdamped Langevin equation** for a free particle ($ U(x) = 0 $):\n", "\n", "$$\n", "m \\dot{v}(t) = -\\gamma v(t) + \\eta(t)\n", "$$\n", "\n", "\n", "\n", "Here:\n", "- $ v(t) $ is the velocity,\n", "- $ \\gamma $ is the friction coefficient,\n", "- $ \\eta(t) $ is Gaussian white noise with:\n", "\n", "$$\n", "\\langle \\eta(t) \\eta(t') \\rangle = C \\delta(t - t')\n", "$$\n", "\n", "We want to compute the **steady-state** value of $ \\langle v^2(t) \\rangle $ and relate it to $ C $, $ \\gamma $, and $ T $.\n", "\n", "\n", " **Step 1: Solve the Langevin Equation**\n", "\n", "We write the equation in standard ODE form:\n", "\n", "$$\n", "\\dot{v}(t) + \\frac{\\gamma}{m} v(t) = \\frac{1}{m} \\eta(t)\n", "$$\n", "\n", "This is a **linear inhomogeneous ODE**, solvable by integrating factor:\n", "\n", "$$\n", "v(t) = v(0) e^{-\\gamma t / m} + \\frac{1}{m} \\int_0^t e^{-\\gamma (t - s) / m} \\eta(s) \\, ds\n", "$$\n", "\n", "\n", "\n", "**Step 2: Compute $ \\langle v^2(t) \\rangle $**\n", "\n", "We take the expectation:\n", "\n", "$$\n", "\\langle v^2(t) \\rangle\n", "= \\left\\langle \\left[ v(0) e^{-\\gamma t / m} + \\frac{1}{m} \\int_0^t e^{-\\gamma (t - s)/m} \\eta(s) \\, ds \\right]^2 \\right\\rangle\n", "$$\n", "\n", "This expands to:\n", "\n", "$$\n", "\\langle v^2(t) \\rangle\n", "= v^2(0) e^{-2\\gamma t / m}\n", "+ \\frac{2 v(0) e^{-\\gamma t / m}}{m} \\left\\langle \\int_0^t e^{-\\gamma (t - s)/m} \\eta(s) \\, ds \\right\\rangle\n", "+ \\frac{1}{m^2} \\left\\langle \\left( \\int_0^t e^{-\\gamma (t - s)/m} \\eta(s) \\, ds \\right)^2 \\right\\rangle\n", "$$\n", "\n", "\n", "- The second term is zero since $ \\langle \\eta(s) \\rangle = 0 $,\n", "- For the third term, use the noise correlation:\n", "\n", "$$\n", "\\left\\langle \\left( \\int_0^t e^{-\\gamma (t - s)/m} \\eta(s) \\, ds \\right)^2 \\right\\rangle\n", "= \\int_0^t \\int_0^t e^{-\\gamma (2t - s - s')/m} \\langle \\eta(s) \\eta(s') \\rangle \\, ds \\, ds'\n", "= C \\int_0^t e^{-2\\gamma (t - s)/m} \\, ds\n", "$$\n", "\n", "Change variable: $ u = t - s \\Rightarrow du = -ds $\n", "\n", "$$\n", "= C \\int_0^t e^{-2\\gamma (t - s)/m} \\, ds\n", "= C \\int_0^t e^{-2\\gamma u/m} \\, du\n", "= \\frac{C m}{2\\gamma} \\left( 1 - e^{-2\\gamma t/m} \\right)\n", "$$\n", "\n", "So the full result is:\n", "\n", "$$\n", "\\langle v^2(t) \\rangle = v^2(0) e^{-2\\gamma t/m} + \\frac{C}{2 m \\gamma} \\left( 1 - e^{-2\\gamma t/m} \\right)\n", "$$\n", "\n", "\n", "\n", "**Long-Time Limit**\n", "\n", "As $ t \\to \\infty $, the exponential terms vanish:\n", "\n", "$$\n", "\\boxed{ \\langle v^2 \\rangle = \\lim_{t \\to \\infty} \\langle v^2(t) \\rangle = \\frac{C}{2 m \\gamma} }\n", "$$\n", "\n", "\n", "**Apply Equipartition**\n", "\n", "At thermal equilibrium, equipartition gives:\n", "\n", "$$\n", "\\left\\langle \\frac{1}{2} m v^2 \\right\\rangle = \\frac{1}{2} k_B T\n", "\\quad \\Rightarrow \\quad\n", "\\langle v^2 \\rangle = \\frac{k_B T}{m}\n", "$$\n", "\n", "Matching with the derived expression:\n", "\n", "$$\n", "\\frac{C}{2 m \\gamma} = \\frac{k_B T}{m}\n", "\\quad \\Rightarrow \\quad\n", "C = 2 \\gamma k_B T\n", "$$\n", "\n", "$$\n", "\\boxed{ \\langle \\eta(t) \\eta(t') \\rangle = 2 \\gamma k_B T \\, \\delta(t-t') }\n", "$$\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Molecualr Dynamics of Harmonic oscillator (NVT)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Our goal is to Numerically solve the **Langevin equation**:\n", "\n", "$$\n", "m \\ddot{x} = -\\nabla V(x) - \\gamma \\dot{x} + \\sqrt{2\\gamma k_BT} \\, \\eta(t),\n", "$$\n", "\n", "- The target is to sample the **canonical distribution**:\n", "\n", "$$\n", "P(x, v) \\propto \\exp\\left( -\\beta \\left( \\frac{1}{2} m v^2 + V(x) \\right) \\right).\n", "$$\n", "\n", "- There are many ways of simulating langevin. Simplest is to use Euler method just like we did for diffusion.\n", "- Here we will use **BAOAB** scheme which splits the Langevin evolution operator into three pieces:\n", " - **B**: Free drift in position ($x$ evolution under $v$).\n", " - **A**: Velocity change under conservative force $F = -\\nabla V(x)$.\n", " - **O**: Stochastic process (friction + random kicks in $v$).\n", "\n", "- Thus, **BAOAB** updates **position and velocity** via a **B → A → O → A → B** sequence at each timestep $\\Delta t$.\n", "\n", "| Feature | Explanation |\n", "|:--------|:------------|\n", "| Stability | Handles both large and small friction $\\gamma$ robustly |\n", "| Accuracy | Samples correct Boltzmann distribution up to $\\mathcal{O}(\\Delta t^2)$ errors |\n", "| Efficiency | Large stable timesteps allowed compared to naive Euler schemes |\n", "| Universality | Reduces to velocity Verlet (Hamiltonian MD) if $\\gamma=0$ and to overdamped Langevin if $\\gamma \\to \\infty$ |\n", "\n", ":::{admonition} **BAOAB scheme**\n", ":class: tip, dropdown\n", "\n", "1. **B: Drift half-step in position**\n", "\n", " $$\n", " x \\to x + \\frac{\\Delta t}{2} v\n", " $$\n", "\n", "2. **A: Kick velocity half-step using force**\n", "\n", " $$\n", " v \\to v + \\frac{\\Delta t}{2m} F(x)\n", " $$\n", "\n", "3. **O: Apply stochastic friction and noise to velocity**\n", "\n", " $$\n", " v \\to e^{-\\gamma \\Delta t} v + \\sqrt{k_B T (1 - e^{-2\\gamma \\Delta t})/m} \\, R,\n", " $$\n", " where $R$ is a standard normal random number.\n", "\n", "4. **A: Kick velocity half-step again**\n", "\n", " $$\n", " v \\to v + \\frac{\\Delta t}{2m} F(x)\n", " $$\n", "\n", "5. **B: Drift half-step in position**\n", "\n", " $$\n", " x \\to x + \\frac{\\Delta t}{2} v\n", " $$\n", "\n", ":::" ] }, { "cell_type": "code", "execution_count": 197, "metadata": { "slideshow": { "slide_type": "slide" }, "tags": [ "hide-input" ] }, "outputs": [], "source": [ "def langevin_md_1d(x, v, dt, kBT, gamma, t_max, en_force):\n", " '''Langevin dynamics applied to 1D potentials\n", " Using integration scheme known as BA-O-AB.\n", " INPUT: Any 1D function with its parameters\n", " '''\n", " \n", " times, pos, vel, KE, PE = [], [], [], [], []\n", " \n", " t = 0 \n", " for step in range(int(t_max/dt)):\n", " \n", " #B-step\n", " pe, F = en_force(x)\n", " v += F*dt/2\n", " \n", " #A-step\n", " x += v*dt/2\n", "\n", " #O-step\n", " v = v*np.exp(-gamma*dt) + np.sqrt(1-np.exp(-2*gamma*dt)) * np.sqrt(kBT) * np.random.normal()\n", " \n", " #A-step\n", " x += v*dt/2\n", " \n", " #B-step\n", " pe, F = en_force(x)\n", " v += F*dt/2\n", " \n", " ### Save output \n", " times.append(t), pos.append(x), vel.append(v), KE.append(0.5*v*v), PE.append(pe) \n", " \n", " return np.array(times), np.array(pos), np.array(vel), np.array(KE), np.array(PE)\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Run lagevin dynamics of 1d harmonic oscillator" ] }, { "cell_type": "code", "execution_count": 205, "metadata": { "slideshow": { "slide_type": "slide" }, "tags": [ "hide-input" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAEYCAYAAADMEEeQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABO/ElEQVR4nO3deXjU9bn//+edmawkEFaBBAir7AQIyCKu1YJaleNysFZrq8dD9/Yc+5XTU7vor6d67KKn1VKrVqtFtK4ouOGCIEsJCIYdZE1YEnZC9uT+/TEzIYSEJCQz71nux3XNNdsnM68h4ZM771VUFWOMMcYYEzviXAcwxhhjjDGhZQWgMcYYY0yMsQLQGGOMMSbGWAFojDHGGBNjrAA0xhhjjIkxXtcBWqpLly6alZXlOoYxJsKsWrXqoKp2dZ3jXNm5zxhzLho790VcAZiVlUVubq7rGMaYCCMiu1xnaA079xljzkVj5z7rAjbGGGOMiTFWABpjjDHGxBgrAI0xxhhjYowVgMYYY4wxMSZoBaCIPC0ihSKyronjxolItYjcGKwsxhhjjDHmlGC2AD4DTD3bASLiAR4C3g1iDmOMMcYYU0fQCkBV/QQ43MRh3wNeAQqDlcMYY0JJRKaKyGYR2SYisxp4voOIvCkia0VkvYh8w0VOY0xsczYGUEQygOnA7GYce7eI5IpIblFRUfDDGWPMOfD3ajwGTAOGAreIyNB6h30H2KCqo4BLgN+KSEJIgxpjYp7LSSCPAPeqanVTB6rqE6qao6o5XbtG7EL+xpjoNx7YpqrbVbUCmAtcV+8YBdJERIBUfD0lVaGNaYyJdS4LwBxgrojsBG4EHheR6x3mMVFq06ZN/PjHP2bmzJkcPHjQdRwT3TKAPXXu5/sfq+uPwBBgL5AH/EBVaxp6Mev9MPVVVFTw8ssvc+utt/L222+7jmMimLOt4FS1b+C2iDwDvKWqr7vKY6JLaWkpL7/8Mn/5y19YvHhx7eNvvvkmzz//PJdeeqnDdCaKSQOPab37XwbWAJcB/YH3RWSxqh4/4wtVnwCeAMjJyan/OiaGbN26lSeffJK//vWvBP4YeO2111i+fDkjR450nM5EomAuA/MCsAw4X0TyReROEZkpIjOD9Z7GADzyyCP07NmT22+/ncWLF9OuXTvuuusuJk+ezN69e7n88sv56U9/SmVlpeuoJvrkA73q3M/E19JX1zeAV9VnG7ADGByifCbCHDp0iKlTpzJo0CD+93//l6KiIoYPH87ll19OaWkpN9xwA0ePHnUd00SgoLUAquotLTj2jmDlMLHlgw8+4Ec/+hEAOTk53H333TywIZ33E1PQSV+hg87lxPIX+dWvfsWHH37InDlzyMrKchvaRJOVwEAR6QsUADOAr9Y7ZjdwObBYRM4Dzge2hzSliQxzhPv/Bu++CymJ8K8T4O5L4YIB6yirXMekbbBm2zZuv/12Xn/9deLibG8H03z202KixokTJ7jzzjsB+MUvfsHKlSv5t3/7N+ISUwCQOA/pU27lww8/JDMzk2XLljF27Fjy8/NdxjZRRFWrgO/iW9t0I/CSqq6v1/vxADBJRPKAD/BNhrPBqeYM+47AEx/6bi/9OTx9N0wYCCKQnACv/BA6tvMNbfn1r3/tNKuJPFYAmqjx4x//mF27djFmzBh+8pOfNHrc198uhn95mMTeIzh8+DCDp90RupAm6qnqAlUdpKr9VfVX/sdmq+ps/+29qnqlqo5Q1eGq+rzbxCZc/WY+lFXC9BwY1efM5/t1g79/G0SE++67j/feey/0IU3EsgLQRIX33nuPP//5z8THx/Pss88SHx9/1uM9yWl0nvp9iPNyct2HfPbZZyFKaowxTSssLORPH/hu3ze98eOmZft6PFSVW265hZ07d4YinokCVgCaiHfs2DHuuusuAH75y18yfPjwZn1dfMcepI25GlDuueceVG2SpTEmPPz2t7+ltAK+MgZGZ5392J/+9KdcddVVHD58mO9///shyWcinxWAJuLdc8897Nmzh3HjxvHjH/+4RV/bYdIM4pJS+fDDD1mwYEGQEhpjTPMdPHiQxx57DICfnaX1LyBuroenr1mAJw4WzH+TotkNrUZkzOmsADQR7Z133uHJJ58kISGBZ555Bq+3ZRPbPclpdJg0A/CNIayqsg0ZjDFu/f73v+fkyZNclQ05/Zr3Ned1gCtHQHUN/GNFUOOZKGEFoIlY1dXVfOtb3wLggQceYOjQ+luuNk/amKvp378/Gzdu5Mknn2zLiMYY0yKHDx/mD3/4AwD3Xd+yr711su/675+2bSYTnawANBGrx7/ez86dO/Gm9+Cxg4PJmjX/nF5HPPE8+OCDAPz85z/n+PEzNmQwxpiQePTRRzlx4gRXXnklEwa27GuvG+tbL3DpVtixY0dwApqoYQWgiVjFa3z7YKZmfxmJ87TqtW644QYmTZpEYWEhDz30UFvEM8aYFjl69CiPPvooAD/72c9a/PWpSXD9WN/tOXPmtGU0E4WsADQRac+ePZRuXwVxXlKHf6nVr9f3vxawvd+/APDr//0Nvb5rS7MZY0LrhRde4NixY1x66aVMnjz5nF7jq5N813//+99tZQNzVlYAmoj05JNPgtaQMmginnbpbfKaiRmDSR4wHq2q4OSGj9rkNY0xplnmCAv+8m0Abhv4Ecw5t5m8V46ALmmwceNG1qxZ04YBTbSxAtBEnKqqqtrJGmmjp7XpawdaE4vXWwFojAmd0gr4YL3v9tSR5/468V64+QLfbesGNmfTsjUzjAkDb731Fnv37sXbKZPEXiNOe+5cJ4IEJPcfR1xSKpWFO8jLy2PEiBFNf5ExxrTSoo2+InBMFvTo2LrXunUyPL7Q16X84IMP4vG0boy0iU7WAmgizp///GcA0kZ9GZG2XfBUvPGkDJ4CwHPPPdemr22MMY1ZsMZ3fVV2619r4kDo27cvBQUFfPLJJ61/QROVrAA0EWXnzp28++67JCYm0m7E5UF5j3bDLgN8g6irq6uD8h7GGBOgqsxf47vdFgWgCHz1q18FfOcxYxpiBaCJKH/5y19QVW688UY8ye2D8h6JGYPxpvdg7969fPjhh0F5D2OMCdi6dSvbC6FTKozv3zaveeuttwLw8ssvU1ZW1jYvaqKKFYAmYlRWVvLUU08BMHPmzKC9j4jQbtilgHUDG2OCb/5839jlqSPB00a/lYcMGUJ2djbHjh2zfc5Ng6wANBHjjTfe4MCBAwwdOvSc18hqrkAB+Morr1BcXBzU9zLGxLZAgdYW3b91BVoB586d27YvbKKCFYAmYgQmf/z7v/97m0/+qC++Yw8mT55MSUkJr732WlDfyxgTu4qLi1m0aBEirVv+5QxzhOmeHwPw4dv/oOb54J4zTeSxAtBEhKKiIj744AMSEhK47bbbQvKegfexbmBjTLB88MEHVFZWMmEAdE5r29fu1w0yOsKhYthQ0LavbSKfFYAmIgy5439QVeIyhjP6oaWtXu+vOW6++WYSEhJYuHAhBQV29jTGtL3a7t9Rbf/aInDxEN/tTza1/eubyGYFoIkIpdv+CUDKgAtC9p6jH1qKNysHVWX41+4LSdFpjIkdqhq08X8BgQJw0cbgvL6JXFYAmrBXVlZG6Y7VACQPGBfS92433Lcm4EnbGs4Y08by8vLIz8+ne/fuZPcJzntcNNh3vWiTr+A0JsAKQBP2Pv74Y7SyjPhu/fC27xbS907uN5a45PZUFu2komhnSN/bRCYRmSoim0Vkm4jMauD5H4vIGv9lnYhUi0gnF1mNW4HWv2nTphEXpN/G5/eA8zrAgWOwZcuW4LyJiUhWAJqwN2/ePABSBowP+XuLJ55k//uWfpEb8vc3kUVEPMBjwDRgKHCLiAyte4yqPqyq2aqaDfwXsEhVD4c8rHGutvv3qquC9h4idVoBFy0K2vuYyBO0AlBEnhaRQhFZ18jzt4rI5/7LUhEJwhBYE+lUlTfffBOA5BCO/6sruV8OAKXbrQA0TRoPbFPV7apaAcwFrjvL8bcAL4QkmQkrR44cYenSpXg8Hq644oqgvtfFVgCaBgSzBfAZYOpZnt8BXKyqI4EHgCeCmMVEqDVr1pCfn48ntRMJ3dtoj6QWSs7KBomjvGAjx48fd5LBRIwMYE+d+/n+x84gIin4zpGvNPZiInK3iOSKSG5RUVGbBjVuLfxFJ6qrq7lwUDUd5qcH9b1qJ4IsWmTjAE2toBWAqvoJ0Gi3hqouVdUj/rvLgcxgZTGRq7b1r/94RNyMWIhLSiUxYzDUVLNw4UInGUzEaGi13cZ+434F+PRs3b+q+oSq5qhqTteuXdskoAkPn/qH41029OzHtYWhGdA5FQoKCtixY0fw39BEhHAZA3gn8HZjT9pfwbHrVPdv6Mf/1RXoBn777UZ/TI0BX4tfrzr3M4G9jRw7A+v+jVnLtvquJw4M/nvFxdk4QHMm5wWgiFyKrwC8t7Fj7K/g2LR3715yc3NJTk4mqY/bIaLJ/cYCvgLQulDMWawEBopIXxFJwFfkzat/kIh0AC4G3ghxPhMGysrK+Gynb4LG+BCNbLEC0NTntAAUkZHAk8B1qnrIZRYTft566y0ArrjiCuLiE51mie/WD0+7jhQUFLBuXYPzmoxBVauA7wLvAhuBl1R1vYjMFJGZdQ6dDrynqidd5DRurV69mspqX9dsh5TQvGfdcYDGgMMCUER6A68Ct6mqLU5kzhBY/uXaa691nAREhKS+p1oBjWmMqi5Q1UGq2l9Vf+V/bLaqzq5zzDOqOsNdSuPS8uXLAZg4IHTvObI3dOjQgZ07d7J79+7QvbEJW8FcBuYFYBlwvojki8id9f4K/hnQGXjcvyCqrbFhap08eZIPPvgAgKuvvtpxGp+63cDGGHOuli1bBsCEEIz/C/DEwZQpUwBrBTQ+wZwFfIuq9lDVeFXNVNWn6v4VrKp3qWrHwIKoqpoTrCwm8ixcuJCysjIuuOACunfv7joOAEl9RxMXF8eSJUtsORhjzDlz0QIIcPHFFwPwySefhPaNTVhyPgnEmIYEZv9+5StfcZzkFE9SKhMnTqSqqqq2ddIYY1oiPz+f/Px8OqTA4J6hfe9AAWgtgAasADRhSFVrJ4CEw/i/uqZNmwZYN7Ax5twEWv8u6E/Q9v9tzOjRo0lLS2Pr1q3s27cvtG9uwo4VgCbsrF+/ngMHDtCzZ0+GDx/uOs5p6haAthyMMaalAgXghBB3/wJ4vV4mT54MWCugsQLQhKGPPvoIgMsuuwyRhjZWcGf6iwXEtUsnPz+fjLseJ2vWfNeRjDERJDABJBQLQDfkoosuAqwANFYAmjB07x/mAPD2oU5kzZofVkWWSBzJ/uVgSrevcpzGGBNJKioqWLXKd964wEELIHOEi8p+AsCn82fDnPD6A9uElhWAJqxUV1dTvse30HJi75GO0zQssByMFYDGmJZYs2YN5eXlDB48mI7t3GQYneVbEmZDAZSUu8lgwoMVgCasrF27lpqyYjwdziM+PTyWf6kvKWs0SBzl+RuoKS9xHccYEyFqx/9NmOAsQ0oiDMuE6hpYa+tBxzQrAE1YCYz/S+o9wnGSxnmS00jsMQhqqijbk+c6jjEmQtSO/5s40WmOnL6+69ztTmMYx6wANGHlww8/BCCpzyjHSc4usY+ve7p8t+0LbIxpnnBoAQTI6ee7XrXDaQzjmBWAJmxUVlbWrlAfzi2AAEm9fPnK9lgBaIxp2v79+9m5cyepqakMGzbMaZaxgRZAKwBjmhWAJmysWrWK4uJivJ0y8KZ1cR3nrBIzhkCch4oDX9i2cMaYJgVa/8aPH4/H43GaZWQv8HpgYwEUFxc7zWLcsQLQhI1IGP8XEJeQREL3AaA1fPrpp67jGGPCXLiM/wNISoARvaBGfTOTTWyyAtCEjdrxf2G6/Et9gW7gjz/+2G0QY0x4myMsn/e/AEwo+1VYrL9XOxEkN9dtEOOMFYAmLJSXl9e2pEVCCyCcymkr6htjzqaqGlb6Z9xOcLQDSH2BiSBWAMYuKwBNWFixYgWlpaUMGzYMT7uOruM0S2LGEJA4cnNzbRyNMaZRn++G0goYcB50SXOdxsdaAI0VgCYs1N3/N1LEJaaQ0L0/1dXVLF261HUcY0yYWvGF73qCi+3fGjG8FyR4YfPmzTaRLUZZAWjCQmD836WXXuo4ScvYOEBjTFNW+5dbCXS7hoMEL4zq7bu9evVqt2GME1YAGudKSkpYvnw5IsLFF1/sOk6LJNo4QGNMEz7b5bsek+U0xhlqF4ReZfuaxyIrAI1zS5cupaKiguzsbDp16uQ6ToskZQ4lLi6OlStXUlJi+wIbY05XWVlJ3h7f7UCLW7gYa+MAY5oVgMa5SBz/FxCX2I7s7GwqKytr1/kysU1EporIZhHZJiKzGjnmEhFZIyLrRcSaj6PYhg0bqKjyTQBpn+I6zelsIkhsswLQOBep4/8CLrnkEsDGARoQEQ/wGDANGArcIiJD6x2TDjwOXKuqw4CbQp3ThM5nn30GwOgstzkaMjQDkpKS2LZtG0eOHHEdx4SYFYDGqZMnT7Jy5Uri4uKYMmWK6zjnJDBu0cYBGmA8sE1Vt6tqBTAXuK7eMV8FXlXV3QCqWhjijCaEAhMswm38H0C8F7KzswGbCBKLrAA0TvW761Gqq6vxdu3LyP9ZTNas+a4jtdiUKVMQkdq1DE1MywD21Lmf73+srkFARxH5WERWicjtjb2YiNwtIrkikltUVBSEuCbYwrkFECAnJwewbuBYZAWgcao8fwMAiZlDmzgyfI1+aCnerllUVFTQ95uPRGQRa9pMQ3t8ab37XmAscDXwZeA+ERnU0Iup6hOqmqOqOV27dm3bpCboampqavfaHd3HbZbGWAEYu6wANE6dKgCHOU7SOoH1AMv25DlOYhzLB3rVuZ8J7G3gmHdU9aSqHgQ+AUaFKJ8JoS+++ILi4mIyOkK3Dq7TNMwKwNhlBaBxpqqqivK9mwD/tmoRLKn3cADK9qxznMQ4thIYKCJ9RSQBmAHMq3fMG8AUEfGKSApwAbAxxDlNCATG1YVr9y/A4MGDSUlJYefOnRw6dMh1HBNCQSsAReRpESkUkQZ/I4rP//mXSvhcRMYEK4sJT3l5eWhFKd707njTOruO0yqBFsyKvZvRqgrHaYwrqloFfBd4F19R95KqrheRmSIy03/MRuAd4HPgn8CTqmp/OUShcB//B+DxeBgzxvfr1xaEji3eIL72M8Afgb818vw0YKD/cgHwJ/+1iRFLliwBIr/1D8CT0oH4Ln2oPLiL8n1bXMcxDqnqAmBBvcdm17v/MPBwKHOZ0AsUgOE4A7jWHGFsKiwBcp/6MlceBL5af9iqiUZBawFU1U+Aw2c55Drgb+qzHEgXkR7BymPCT20BGOHj/wISe/k+R3mB9eYZE+tU9VQXcJhOAAkILAi9crvbHCa0XI4BbM5yCYAthRCNVLVOARi5M4DrCrRkBia2GGNiV0FBAQcPHqRjx4707uI6zdkFtoT7bKfTGCbEXBaAzVkuwfegLYUQdXbt2sXevXuJS0ojvnOm6zhtIslfyJYXbKKmpsZxGmOMS7Xj/0aPRhr6bRdGBvWA5ATYdRCOnHSdxoSKywKwOcslmCh1qvVvCCLRMRnd074bntRO1JSdYPPmza7jGGMcqt0BZEz4z2/0xMEI/2/jz3e7zWJCx+Vv3nnA7f7ZwBOAY6q6z2EeE0KnJoBER/cvgIjUdgN/+umnjtMYY1yq2wIYCUb19l2v2eU2hwmdYC4D8wKwDDhfRPJF5M66SyHgmyW3HdgG/AX4drCymPATKJCiZfxfQKCgtQLQmNgWaQVgtn+iihWAsSNoy8Co6i1NPK/Ad4L1/iZ8HTlyhHXr1pGYmEhi94Gu47SpxExrATQm1h06dIjdu3eTkpLCoEGD4DPXiZpmBWDsiY7BVyaiLF26FIBx48Yh3njHadpWQrd+SHwiW7duxWasGxObAq1/o0aNwuPxOE7TPIExgBsKoKLCFrOPBVYAmpALjP+bPHmy4yRtTzxeEnoMAk4VusaY2BJp3b8Aackw4DyoqIJNmza5jmNCwApAE3KB7tELL7zQcZLgsHGAxsS2SJoBXFdtN/CaNU5zmNCwAtCEVHl5Of/85z8BmDRpkuM0wZGUMRiwAtCYmDRH+GzRXABG778L5oT5IoB1WAEYW6wANCG1atUqysvLGTZsGJ06dXIdJygS/EvB5ObmUlZW5jiNMSaUistgy37wemBYhK1xH1gKZu3atW6DmJCwAtCEVGD8X7R2/wJ4klIZNmwYFRUVrFq1ynUcY0wIfb4bVGF4JiRG2By3ui2AvoU6TDSzAtCEVKBbNBongNQV+HzWDWxMbAnspzs6y2WKc5PRCTqnwuHDh8nPz3cdxwSZFYAmZFQ16ieABAQKQJsJbExsCayjl93bbY5zIQKjbBxgzLAC0ITMli1bOHToED169CArK8t1nKAKTHBZunSpdaUYE0PW+vfSDRRSkSbbxgHGDCsATcgEWsMmTZqESOTMjDsX/fv3p1u3bhQVFbF161bXcYwxIVBVVUXeHt/tkRHYAgg2EziWWAFoQqZuARjtRMTGARoTY7Zu3UpZJfTpAh3buU5zbqwLOHZYAWhCJpYKwKxZ8/n4aEcAvv/IXLJmzXecyBgTbIFu01ER2voHMLgnJCQk8MUXX3D8+HHXcUwQWQFoQuLIkSNs2LCBxMTEiNoeqTUCO4KU5W9wnMQYEwq1BWCEjv8DSPDCsGHDAMjLy3OcxgSTFYAmJJYvXw5ATk4OiYmJjtOERkL3/uCJp+pwPtWl9pe0MdEu0G0ayS2AANnZ2YB1A0c7KwBNSMRS92+AeOJJ7DEQgPKCjY7TmFARkakisllEtonIrAaev0REjonIGv/lZy5ymrYXDV3AAKNGjQKsAIx2XtcBTGxYtmwZAM9si+elGBoPl5gxlPL8DZQXbHIdxYSAiHiAx4ArgHxgpYjMU9X64wAWq+o1IQ9ogqawsJB9+/aRmgT9urlO0zrWAhgbrAXQBF1VVRUrVqwAILHnEMdpQivRvy+wtQDGjPHANlXdrqoVwFzgOseZTAgEWv9G9oK4CP/NGmgBXLduHVVVVY7TmGCJ8B9TEwnWrVtHcXEx3vTueFI7uo4TUokZgwGo2LeVyspKx2lMCGQAe+rcz/c/Vt9EEVkrIm+LyLDGXkxE7haRXBHJLSoqauuspg1FwwSQgPT0dPr06UNZWRlbtmxxHccEiRWAJugC4/8CrWGxxJPSAW/HnmhVua2sHxsaWuG8/lYwq4E+qjoK+APwemMvpqpPqGqOquZ07dq17VKaNhct4/8CrBs4+lkBaIIulgtAOPW5A+MgTVTLB3rVuZ8J7K17gKoeV9Vi/+0FQLyIdAldRBMMgQIwOwpaAJkjZHveAGDNnFthTnTv3BSrrAA0QXeqABzsOIkbgQIw8O9gotpKYKCI9BWRBGAGMK/uASLSXfx7IYrIeHzn4UMhT2raTHl5ORs3bkREGJ7pOk3bqN0SbpfbHCZ4bBawCar9+/ezY8cOUlNTie8SDX8at1yg8LUCMPqpapWIfBd4F/AAT6vqehGZ6X9+NnAj8C0RqQJKgRmqWr+b2ESQDRs2UFVVxaBBg2iXFB1j5gJd2Wt2g2rDYxtMZLMC0ARVoNtzwoQJbI3zOE7jRnyX3khCCrt37yY/P5/MzChpIjAN8nfrLqj32Ow6t/8I/DHUuUzw1I7/GzUKiI4CMKsrtE+GouOw/yj0cB3ItDnrAjZBFWj1mjhxouMk7ojEkdjzfMDGARoTjQITJQITJ6KByKlWwLW73WYxwWEFoAmqWNwBpCE2DtCY6HV6C2D0sHGA0S2oBWAztkTqICJv+tfDWi8i3whmHhNa5eXl5ObmAr4u4FhmBaAx0UlVo74AtBbA6BS0ArDOlkjTgKHALSIytN5h3wE2+NfDugT4rX/mnIkCq1evpqKigmHDhpGenu46jlOJPc9HRFi9ejWlpaWu4xhj2kh+fj5HjhyhU6dOZGQ0tOZ35BplLYBRLZgtgM3ZEkmBNP+SCKnAYcD2nYkS1v17SlxiCiNGjKCqqqq2VdQYE/nqjv/zr+4TNYZlgCcOtuyDkpIS13FMGwtmAdicLZH+CAzBt1BqHvADVa2p/0K2HVJksgLwdIF/B+sGNiZ6RGv3L0BSAgzpCTXq29LTRJdgFoDN2RLpy8AaoCeQDfxRRNqf8UW2HVLEUVUrAOsJ/DvYTGBjokc0F4BQpxvYtoSLOsEsAJvcEgn4BvCq+mwDdgCxuV1ElNm5cyf79++nc+fODBw40HWcsDBrSTkAb77/MX3ufYusWfMdJzLGtFagMIrWAjDbCsCoFcwCsMktkYDdwOUAInIecD6wPYiZTIh8+umngK/VK9rGxZwrb3p34lLSqSk5RtXRfa7jGGNaqfhp4YsvthHvgaGfj47KPXNrZwL7WzpN9AhaAaiqVUBgS6SNwEuBLZEC2yIBDwCTRCQP+AC4V1UPBiuTCZ1AATh58mTHScKHiNRuC1desNFxGmNMa+Xt8W2TNiQDEqJ0X63axaDXrqWm5owh+iaCBfVHthlbIu0FrgxmBuOGFYANS8wYQunW5ZQXbCR1+OWu4xhjWmGtf3mUQJEUjbq2h54dYe+Rk2zfvp0BAwa4jmTaiO0EYtrc0aNHWbduHfHx8eTk5LiOE1YCC0KXF2xynMQY01qf+QvAQDdptAoUuDYOMLo0qwVQRDLxjeGbgm/GbimwDpgPvN3Q0i0mdi1fvhxVRbr2Y/AvPnAdJ6wkdh8AcV4qi3ZRU37SdRxzFnbeM00JLJA8OsoLwOw+8PZaXwF44403uo5j2kiTLYAi8lfgaaACeAi4Bfg2sBCYCiwRkYuCGdJElkD3b1JG/Y1fjHgTSOjeH1BrBQxjdt4zTamqquJz/xZpo2KgAASbCBJtmtMC+FtVbWgFyHXAq/4ZvlE8AsK0VKAATMwc4jhJeErMGELF3s1WAIY3O++Zs9q8eTNlldCnC3RKdZ0muKwLODo12QIYOAmKSLf6z4nI+apa4V/DzxiqqqpYsWIFcGq8mzldUqavZbS8YL3jJKYxdc5700UksYHn7bwX4wLF0OgspzFCYkB3SElJIT8/n0OHDrmOY9pISyaBLBaRmwN3ROQ/gdfaPpKJZGvXrqWkpARvxx542nV0HScsJfq7xsv3bqaystJxGtOEa4EtIvKciFwtIlG62Idpqc8++wyI/vF/4NsPeOTIkYB1A0eTlhSAlwC3icg/ROQTYBAwPiipTMSq7f618X+N8rRLx9uxJ1pZbifTMKeq3wAGAP8Avgp8ISJPuk1lwkGgAIz2GcABgZ1OrBs4ejS7AFTVfcA7wEQgC/ibqhYHKZeJUKcKQOv+PZtEfzfwkiVLHCcxTVHVSuBtYC6wCrjObSLjmqrGVBcwQHZ2NmAtgNGk2QWgiLwPXAAMB64Cfi8ivwlWMBN5VNUKwGZKsgIwIojIVBF5BtgG3Ag8CfRwGso4t2fPHg4fPkznVMjs5DpNaAQKQGsBjB4t6QJ+TFVvV9Wj/gHSk4BjQcplItDu3bspKCggPT2d+C69XMcJa4Eu8k8//RRVdZzGnMUdwOvAIFX9uqou8G9z2Sh/0bhZRLaJyKyzHDdORKpFxBZWizB1u39jZavzESNGICJs2LCB8vJy13FMG2jOOoACoKqv131cVatU9YG6x5jYFmj9mzRpEiK2yczZeDtlEJfcnv3797N9+3bXcUw9dc57M1T1dVU94zdeQ+c9EfEAjwHTgKHALSJyxoBY/3EP4dsr3USY2gkgWW5zhFK7du0YOHAgVVVVbNxoe5lHg+b8lv5IRL4nIqeteSUiCSJymYg8C3w9OPFMJLH9f5tPRGwcYHg71/PeeGCbqm5X1Qp84wYbGjP4PeAVoLCtg5vgqx3/FyMTQACYI4zquAWANX8aDXOs3SfSNacAnApUAy+IyF4R2SAi24Gt+FbH/72qPhPEjCZCWAHYMnW7gU3YqXve29eC814GsKfO/Xz/Y7VEJAOYDswORnATfLHYAgh1dgTZ7TaHaRtNrmmlqmXA48DjIhIPdAFKVfVokLOZCHL8+HHy8vLwer2MGzcO3v7IdaSwZxNBwlcrznsNNYvUH+T5CHCvqlY3NXpGRO4G7gbo3ds2HgkHhw4dYvfu3SQnJzOoR6nrOCEVKAA/2+U2h2kbTRaAIpIEzMS3FtbnwNNNDYI2sWfFihXU1NSQk5NDSkqK6zgRIaF7f5KSkti4cSOHDh2ic+fOriMZv1ac9/KBujOgMoG99Y7JAeb6i78uwFUiUlV/nDWAqj4BPAGQk5Njs4XCQGAZlJEjR+KJW+E4TWgFCsA1u0C14b92TORoThfws/hOWHn4ln/5bVATmYhUdwKIaR7xxKNdBwAw5N8fJWvWfLJmzXecyvid63lvJTBQRPr69wueAcyre4Cq9lXVLFXNAl4Gvt1Q8WfCU2337+jRjpOEXs+O0D0djpXAFwdcpzGt1ZwCcKiqfk1V/4xvHawpQc5kIpCN/zs3gYkg5fkbHCcx9ZzTec/fSvhdfLN7NwIvqep6EZkpIjODF9eESu0SMP518WLN2Czf9eqdLlOYttCcfS1rNytV1Spb8cXUV1VVxfLlywErAFsqKWMIx7ECMAyd83lPVRcAC+o91uCED1W94xzzGUdqZwCPHu1bHjzGjOkL89fAqh1ws+swplWaUwCOEpHj/tsCJPvvC6Cq2j5o6UxEWLt2LcXFxfTt25cePWyThJZIyBgCCOX7t6JVFYg3wXUk42PnPXOG0tJSNm3ahMfjYcSIETFZAAZaAFftcBrDtIHmzAL2hCKIiVyLFi0C4OKLL3acJPJ4klKJ79qHyqKdlO/fSlLmMNeRDHbeMw3Ly8ujurqaYcOGkZyc7DqOE2P7+q5X7/Rt/2m9gpHLtmswrfbJJ58AVgCeq8C+ydYNbEx4C3T/xur4P4CMTtCtPRw5CTt37nQdx7SCFYCmVWpqali8eDFgBeC5sokgxkSGWJ4BHCACY7J8t1etWuU0i2md5owBNKZR69ev5/Dhw3jSunDJ7PWIWBHTUoEFocsLNqJa4ziNMaZBc4TP3vfdHH34Hphzj9s8Do3tC+987isAb7zxRtdxzDmyFkDTKoHxf0m9httYkHPkad8NT2pnasqKqTy4p+kvMMaEXHUNfO7/75kdS3sAN6B2HODq1W6DmFaxAtC0SmD8X2Kv4Y6TRC4RqdMNvN5xGmNMQ7bsg9IK6N0ZOqW6TuNWoABctWoVqrZBTaQKagEoIlNFZLOIbBORWY0cc4mIrBGR9SKyKJh5TNtS1dNaAM25S+o9AoCy3XmOkxhjGvLZTt/16CyXKcJDr87QOfXUvsgmMgWtABQRD/AYMA0YCtwiIkPrHZOOb8P1a1V1GHBTsPKYtrdlyxYKCwuJa5eOt1OG6zgRrbYA3JNnf1EbE4YCO1+MjvHuX/BNBKnbCmgiUzBbAMcD21R1u6pWAHOB6+od81XgVVXdDaCqhUHMY9pYbetfpo3/ay1vp0zi2qVTc/Iomzdvdh3HGFPPyu2+63H93eYIFzYOMPIFswDMAOqOaM/3P1bXIKCjiHwsIqtE5PYg5jFt7NT4P1u8uLVEhKRevlbAjz/+2G0YY8xpqqura3e+GNfPbZZwYS2AkS+YBWBDTUL1+7a8wFjgauDLwH0iMuiMFxK5W0RyRSS3qKio7ZOaFjtt/J+/+9K0TuDf8aOPPnKcxBhT18aNGzlZDlldoattAgicvhagDVuJTMEsAPOBXnXuZwJ7GzjmHVU9qaoHgU+AUfVfSFWfUNUcVc3p2rVr0AKb5tu5cyf5+fl06tSJ+C69XceJCoEC8OOPP7YTqjFh5J///CdgrX91ZXWFjh07UlRURH5+vus45hwEswBcCQwUkb4ikgDMAObVO+YNYIqIeEUkBbgA2BjETKaNBFr/pkyZgoitJtQWvJ0y8bTrSGFhIZs2bXIdxxjjt3LlSsAKwLpEYOzYsYCNA4xUQfvNrapVwHeBd/EVdS+p6noRmSkiM/3HbATeAT4H/gk8qarrgpXJtB3b/7ftiQiJvW0coDHhJlAAjrcJIKcJFIA2DjAyBbXpRlUXqOogVe2vqr/yPzZbVWfXOeZhVR2qqsNV9ZFg5jFtJ9ACeNFFFzlOEl2SrAA0JqyUlZWxdu3a0/bANT5jxowBrACMVNZ3Z1osPz+f7du3k5aWRnZ2tus4UaXuTGAbB2iMe2vXrqWqqoqhGZCW7DpNeKnbAmjnq8hjBaBpsUD374UXXojH43GcJrp4O2XQvXt3GwdoTJiwCSCN69evH+np6Rw4cIB9+/a5jmNayApA02I2/i94RIRLLrkEsG5gY8KBTQBpnIhYN3AEswLQtJiN/wsuKwCNCR82AeQs5ghjkj8EYNUz18Ic2xEqknhdBzCRpdf3nid/0yYkPpEZrx5A3pjvOlLUuT/Xd/3yW+/R5963EBF2Pni121DGxKBjx46xadMmEhISGNm7wnWcsFS7I8gOtzlMy1kLoGmR0p1rAEjMHIZ44t2GiVLeThl4UjtRU3KUykN7mv4CE1ZEZKqIbBaRbSIyq4HnrxORz0VkjX+Howtd5DRNC3Rrjho1igRrLmlQoABcuR1sHkhksQLQtEiZvwBMzsp2miOaiQiJ/tnA5bvzHKcxLSEiHuAxYBowFLhFRIbWO+wDYJSqZgPfBJ4MaUjTbLXdv+PHO04SvgacB51T4cAx2Gk7tUYUKwBNs6kqZTs/AyApa7TjNNEtsB5gmRWAkWY8sE1Vt6tqBTAXuK7uAaparKfWzGjHmXukmzBROwN43DjHScKXCEwc6Lu9bKvbLKZlrAA0zbZx40aqiw8T1y6d+K5ZruNEtdoCcE+era8VWTKAuv32+f7HTiMi00VkEzAfXytgg0Tkbn83cW5RkTWvhFrtDGArAM+qtgDc5jaHaRkrAE2zvf/++wAk98lGxGZ7BZO3Y0//OMBjNg4wsjT0H+OMCl5VX1PVwcD1wAONvZiqPqGqOaqa07Vr17ZLaZq0f/9+9uzZQ1paGueff77rOGFt4gDf9dItbnOYlrEC0DRboAC07t/gqzsOMDDu0kSEfKBXnfuZwN7GDlbVT4D+ItIl2MFMywRa/8aOHWsL3jdhXH/wxMHa3XDy5EnXcUwzWQFomqWioqJ2XbqkrFFuw8SI5L6+Qrtsx2rHSUwLrAQGikhfEUkAZgDz6h4gIgPE34QuImOABOBQyJOas7IJIM2XmgQje0N1DeTm5rqOY5rJCkDTLMuXL+fkyZPEd+6NN80aK0Ihqa9vhf2yPXmUlZU5TmOaQ1WrgO8C7wIbgZdUdb2IzBSRmf7DbgDWicgafDOG/1VtoGfYsQkgLRPoBl62bJnbIKbZrAA0zXKq+9da/0LFm9qJ+K5ZaGU5S5YscR3HNJOqLlDVQaraX1V/5X9stqrO9t9+SFWHqWq2qk5UVfvmhhn9u7Dy03cBGFdwk+1w0QyTBvmuly5d6jaIaTYrAE2z2Pg/N5L7jQXg3XffdZzEmNixvRAOF0O39tDbOjyapXYm8LJltnJBhLAC0DTp6NGjrFy5Eq/XS1Kv4a7jxJRAN7AVgMaEzsrtvutx/Xzr3Jmm9e3qK5gPHjzIF1984TqOaQYrAE2TPvroI2pqapg4cSJxiSmu48SUpIyhSHwSeXl5FBQUuI5jTExY7l/Pbnx/tzkiiYh1A0caKwBNkwLdv1dccYXjJLFHvPG1i0K/9957jtMYExs+2eS7vtCW/2sRmwgSWawANE2yAtCtQDfwO++84ziJMdHv2LFjrNkFXg9MGOA6TWSpOw7QhD8rAM1Z7dy5k23bttGhQwdycnJcx4lJgYkg77//PtXV1Y7TGBPdli5diqpv/F9Kous0kSWnH3i9XvLy8jhx4oTrOKYJVgCaswq0/l122WV4vV7HaWJTfMeeeNO7c+TIETLveISsWfPJmjXfdSxjotInn3wCwEWDHQeJQMkJMHr0aGpqamrXUTThywpA06isWfP54W//BsCikz2s6HAoqa+vFbDUdgUxJqisAGydiRMnAtYNHAmsADSN0ppqynatBSApK9ttmBiX7B8HWLpjleMkxkSvkpISVq5ciQhMHuQ6TWQKFIA2Ezj8WQFoGlWxbws1ZSfwdjgPb3oP13FiWlLvERDnpWLfVqpLbWyNMcGwYsUKKisrye4DHWzFq3MyadIkwLd9aE1NjeM05mysADSNKtm2AoDkAeMRWw3VqbjEFBIzh4DWULZzjes4xkSlQPfvFFv+5Zz1WtyHnh3hyJEjbPmdx7bRC2NWAJpGlW71F4ADJzhOYqBuN7CNAzQmGBYvXgzY+L/WEKmzHuBWt1nM2QW1ABSRqSKyWUS2icissxw3TkSqReTGYOYxzbd161YqD+0hLrEdSZnDXMcxQLJ/IkjZjlW216YxbayioqJ23Jq1ALZO7Y4gW9zmMGcXtAJQRDzAY8A0YChwi4gMbeS4hwDb7DSMvPHGGwAk9c9BPLb8SziI75ZFXLt0qosPU3lwl+s4xkSV1atXU1payuDBg+nWwXWayBZYEHqptQCGtWC2AI4HtqnqdlWtAOYC1zVw3PeAV4DCIGYxLRQoAFMGWPdvuBCJO9UNvM3W2DKmLdUu/3LRRY6TRL6xfX2LaG8ogP1HXacxjQlmAZgB7KlzP9//WC0RyQCmA7PP9kIicreI5IpIblFRUZsHNacrKirydYXEeWt3oTDhIWWgb4mFki22xpYxbckKwLaT4IWL/N3oH6x3m8U0LpgFYENTf+oPXHoEuFdVz7q/lao+oao5qprTtWvXtspnGjF//nxqampI6j2CuERbCyGcJPUdjcQnUrF/K7t373Ydx5ioUF1dzZIlSwCYMmWK4zTR4YoRvuuF69zmMI0LZgGYD/Sqcz8T2FvvmBxgrojsBG4EHheR64OYyTRDbfevzf4NO3HxSbWTQV577TXHaYyJDnl5eRw7dow+ffrQu3dv13GiwpeG+67fz8MmrYWpYBaAK4GBItJXRBKAGcC8ugeoal9VzVLVLOBl4Nuq+noQM5kmlJaW8t577wGQPOACx2lMQ1LOnwzAK6+84jiJMdGhdvkX6/5tM8MzoVt7KDgCmzdvdh3HNCBoBaCqVgHfxTe7dyPwkqquF5GZIjIzWO9rWmfhwoWUlJQwduxYvO27uI5jGpDcfxx4vCxZsoQDBw64jmNMZJsjfPL89wG4KPE5W7i4jcTFweX+FcQWLlzoNoxpUFDXAVTVBao6SFX7q+qv/I/NVtUzJn2o6h2q+nIw85imzZvna6S97rqGJmybcBCXmEJy1mhUlddff911HGMimip8ssl32xaAbluBbmArAMOT7QRiatXU1PDmm28CcO211zpOY84mZZBvNvCrr77qOImpr6kF8EXkVhH53H9ZKiKjXOQ0Plv2QeFxOK8DDOzuOk10CRSAH330EVVVVW7DmDNYAWhqrVixggMHDtCnTx9GjhzpOo45i+QBF+DxePjwww85cuSI6zjGr5kL4O8ALlbVkcADwBOhTWnqCrT+TTnft42ZaTu9u8CgHnD8+HFWrlzpOo6pxwpAA0DWrPlM/dHvADjcZSR9/2uB40TmbDwpHbj44oupqqqqbbU1YaHJBfBVdamqBqr25fhWSDCOvJfnu75kiNsc0epLNg4wbFkBaGqVblsBQLIt/xIR1nh8K63OfOBxsmbNJ2vWfMeJDM1YAL+eO4G3G3vSFsEPrsrKytoC8Kpsp1Gilo0DDF9WABoAKg/uofLQHiSxHUmZw1zHMc0QKNRLd6ympqLUcRrj15wF8H0HilyKrwC8t7EXs0Xwg2vJkiUcL4UhPaFvN9dpotMlQyEuLo5ly5ZRXFzsOo6pwwpAA0Dxho8ASBk0CfF4HacxzeFN60xixhCorqR0+yrXcYxPcxbAR0RGAk8C16nqoRBlM/W89dZbAFwz2nGQKNaxHeTk5FBZWVm73Z4JD1YAGmpqaji53lcApg6/1HEa0xKB2cAlmz91nMT4NbkAvoj0Bl4FblPVLQ4yGr/5833DJq62AjCorrjiCsC6gcONFYCGxYsXU328CE/7riT2Gu46jmmB5EGTACjdnotWVThOY5q5AP7PgM74tr5cIyK5juLGtG3btrF582Y6pMCkga7TRLcvfelLgBWA4cb6+gzPPfccAO2GXoKI/U0QSeLTu5NwXn8qDnxB6Y7PgOmuI8U8VV0ALKj32Ow6t+8C7gp1LnO6QOvf1JEQb78Jg2rirktJTvDtubz/caF7OvBV2x/YNfttH+NKS0v5xz/+AUDqsMscpzHnImXwFABOrvvAcRJjIkdt92+22xyxIDH+1C4rH6x3m8WcYgVgjHvzzTc5fvw4Cd0HEt+lV9NfYMJOu2GXgsRRsm0FtlSIMU07ceIEixYtQkSYlu06TWyoXQ5mndsc5hQrAGNcbffvMJv8Eam8aZ1J7jcWaqprv5/GmMYtXLiQiooKJkyYQJc012liwxX+AvC9PKipcZvF+FgBGMMKCwt555138Hg8tBtykes4phVSR/pm2T311FOo2tgaY86mtvv36qsdJ4kdI3tDr86w9wgs3+Y6jQErAGPa3LlzqaqqYurUqXjapbuOY1ohuf844lI6sGHDBttz05izqKmpYcEC3xyda665xnGa2CECN4333f7HCrdZjI8VgDEs0F142223OU5iWks88bWTeJ5++mnHaYwJX5999hn79u0jMzOTkSNHuo4TU266wHf9jxW+Qty4ZQVgjNq0aRO5ubm0b9+ea6+91nUc0wbajfCttfXCCy9QUlLiOI0x4SnQ/XvVVVch0tDOfSZYLhjg6wYuOALLly93HSfmWQEYowKtfzfeeCPJycmO05i2kNC1Dwk9zuf48eP0vvk+smbNJ2vWfNexjAkfc4T5f/s5ANekPQFzrAAMpdO6gf3Ljxl3rACMQTU1NTz//POAdf9Gm8BkkOK89x0nMSb8HDgG//zCty7dZUNdp4lNgW7gl19+2bqBHbMCMAb1mPH/sXv3bjztu/L1t09YK1EUaTfkIsSbSPnuPCqP7HMdx5iw8vYa3/WlQ6BdktMoMSvQDZyfn2/dwI5ZARiDjq98HYC00Vfb1m9RJi4xhZTBkwEozrN9N42p60V/vfGVMW5zxDIRuNG6gcOC/faPMevWraNs52dIfCKp2VNdxzFBkDrySgBO5i1Ea6odpzEmPBQUFPBeHsR74OYJrtPEtputGzgsWAEYYx555BEA2g3/Ep6kVLdhTFAkZg7D27EH1cWHKN2x2nUcY8LC888/T43CtWOw3T8cu2AA9OrVi/z8fFassEUBXbECMIYUFhb6J38I7XNs6ZdoJSKkjvoyAMf/+arjNMa4p6o888wzANxhmx45J+JbgQLgpZdecpwmdlkBGEP+9Kc/UV5eTvKAccR3ynAdxwRRWvY0JCGF8t15NtDaxLwVK1awadMmzusAX7a1n8PCTTfdBFg3sEtWAMaIsrIyHn/8cQDaj7vebRgTdHGJ7Ugb49vn9Ne//rXjNMa4FWj9+9pkiPe6zWJ8LvhiUu1s4BUPeFzHiUlBLQBFZKqIbBaRbSIyq4HnbxWRz/2XpSIyKph5YtmcOXMoLCwkOzubxF4jXMcxIdA+5zrEm8C8efNYt26d6zjGOFFaWsrcuXMB6/4NJ3FxdWYD2zBAJ4JWAIqIB3gMmAYMBW4RkfpLb+4ALlbVkcADwBPByhPLVJXf//73APzoRz+y7Y9ihKddeu3C0A899JDjNMa48cYbb3Ds2DFycnIY3st1GlNXYFHoF5dDVVWV2zAxKJgtgOOBbaq6XVUrgLnAdXUPUNWlqnrEf3c5kBnEPDFr4cKFrFu3jh49ejBjxgzXcUwItR//L3g8Hl544QV27NjhOo4xIVc7+eOOO5zmMGe6oD8M6gF7j8Brr73mOk7MCWYBmAHsqXM/3/9YY+4E3m7oCRG5W0RyRSS3qKioDSPGhmtn/hcApQO+xKCf2RZhscTb4TySBl9EdXU1I2/4nu36YmJKfn4+7733HgkJCfbHbxiKi4Mf+BYsqF2izIROMAvAhvoZtcEDRS7FVwDe29DzqvqEquaoak7Xrl3bMGL0W7NmDWXbVyFeW/g5VrW/wLfcwsm896k+eaSJo42JEnOE5/6zF6rKtdkVdH63i+tEpgG3T4EOKbB06VJWrlzpOk5MCWYBmA/UHXGRCeytf5CIjASeBK5T1UNBzBNzVJV77rkHgNTsqXhSOjhOZFxI6NqH5IET0KoKjue+4TpO1GvG5LfBIrJMRMpF5B4XGWOBKjyz2HfbJn+Er9Qk+LdLfbcfffRRt2FiTDALwJXAQBHpKyIJwAxgXt0DRKQ38Cpwm6puCWKWmPT222/zwQcfEJfYjg6TrPsjlnWY4Ftz68TqBRw7dsxxmujVzMlvh4HvA78JcbyYsnwbbNkH3dNt7b9w990rIC4ujpdeeom9e89oJzJBErQCUFWrgO8C7wIbgZdUdb2IzBSRmf7DfgZ0Bh4XkTUikhusPLGmqqqqtvWvw6QZeJJt76NYltjzfBJ7j0QrSnj44Yddx4lmzZn8VqiqK4FKFwFjxeP+4c5fmwxeW2YurPXpCtOnT6eyspI//elPruPEjKCuA6iqC1R1kKr2V9Vf+R+braqz/bfvUtWOqprtv+QEM08seeqpp9i4cSP9+vUjbcw1ruOYMJA+5TYAHn74YbZt2+Y4TdRq6eS3s7IJcOdm8+bNzFnqK/y+/SXXaUxz/OAHPwBg9uzZlJWVOU4TG2wnkCh0/PhxfvaznwHw4IMPIt54x4lMOEjKHEK74V+ioqKC73//+6g2OCfLtE6zJ781h02AOzcPPPAANQrfuAj6dnOdxjTHhRdeyJgxYzh48CBz5sxxHScmWAEYZbJmzafX1LspLCwksedg7slNdh3JhJGOl9xBhw4dePvtt5k3b17TX2BaqlmT30zwbNq0iRdeeAGvB35yXdPHm/AgIvzwhz8EfJNB7A/U4LMCMMpUHS/ixErfgpodL7vTdv0wp/G0S+eBBx4AfF0uJSUljhNFnSYnv5ngeuCBB6ipqeGbF0OWNZpGjjnCzTW3c14H+Pzzz/n4p1aeBJv9C0eZo4ufQ6sqSDn/QhIzhriOY8LQb/b0Ib5bX3bt2kXPK+60xaHbUHMmv4lIdxHJB/4D+KmI5ItIe3epo0eg9S8+Pp6fXOs6jWmpxPhTYzZ//47bLLHACsAosnjxYk6u+wjivKRfcofrOCZMSZyHTld8C4BjK16h8oj1ULalZkx+26+qmaraXlXT/bePu00dHe6//35UlW9+85v0sda/iDTzcl8h+OZqWLJkies4Uc0KwChx5MgRbr31VkBpf8ENxKd3dx3JhLGkzKG0G345VFdyeOGfbbyNiXgbN25k7ty5vta/n/zEdRxzjrp1gHv9C1d85zvfoaqqym2gKGYFYBRQVe6++2727NlDQo/zSZ98i+tIJgJ0vOQOJLEdZdtX8fLLL7uOY8y5myPcf+dQVJU7L6qk95I+rhOZVph1rW/85ueff87s2bNdx4laVgBGgaeeeoqXX36ZtLQ0ulz7Y8TjdR3JRABPu450vOh2AO688062bLHNeExk2pAPLy6HBK/N/I0GyQnw+6/5bv/0pz+lsLDQbaAoZQVghNu0aVPtApqPP/64df2aFkkdfRUp50/mxIkTTJ8+neLiYteRjGmRmpoafvCcb+/fuy6BXp1dJzJt4bqxMHXqVI4dO8asWWdsqW3agBWAEazPPa8z8uKrKSkpod2wS/npuo6uI5kIIyJ0nvYD4jv3YsOGDZw37ir63PuWzQw2EeMPf/gDC9dBlzS4b7rrNKatiPjWA4yPj+evf/0ry5cvdx0p6lgBGMGOfPIslYXb8aZ3r53VaUxLxSWm0HX6T5CEZEo2LeZEri1bZyLDunXruPfeewH4y13QPd1tHtO2Bg0aVLun/Xe+8x2qq6sdJ4ouVgBGqGeffZYTK18HiaPLV35MXGKK60gmgsV37kXnq34IwJGPnqJszzq3gYxpQnl5OV/72tcoLy/nrkvgettJPvrMEf57wK/J7ASrV6/mL/9m49vbkhWAEejFF1/km9/8JuDb7SOx5/mOE5lo0O78ybQf/y+gNRS98SAFBQWuIxnTqPvuu4+1a9fSv39/fn+b6zQmWNolwe/8E0Lunetb7se0DSsAI8y8efP42te+Rk1NDR0uvJX2OTblzbSd9Iu/TmLvkdScPMoll1zCrl27XEcy5gwfffQRv/nNb/B4PDz//POkJrlOZILpxvFwwzg4XgrXXHMNRUVFriNFBSsAI8h7773HTTfdRFVVFffeey8dJs1wHclEGYnz0PX6WSSc159t27bRf8Q4Mu5+gqxZ821iiAkLR/8ifP2my1BVfnpdNRO2T3QdyQSZCPztW5DTD7Zv3871119PWVmZ61gRzwrACLFo0SKuv/56Kioq+N73vsevf/1rRMR1LBOFPMntOW/Gr0jMGEL1iSL2z7mXiqKdrmMZQ0lJCTf9H+w5BBf0h/+2DpCYkZII8/4DMjMzWbp0Kd/85jdtB6NWsgIwArz66qtcc801lJaWctddd/HII49Y8WeCKi4plW43309SH1938IE5/0X5vq2uY5kYdvLkSa655hoWroNu7eH5b0O8zQmIKT06wvz580lNTeWFF17gl7/8petIEc0KwDBWUVHBf/zHf3DDDTdQXFxMu+GX8V7Hr9DvJ29bd5wJuriEZLrd+AuSB4ynpuwEB+b+hPnz7efOhN6JEyeYNm0aH330ET3SYdF9MMDWvI9JI0eO5MUXXyQuLo5f/vKXPP/8864jRSz7+ylM7dmzh3/9139l2bJleL1e0i76Bmk511rLnwkp8SbQ9fqfcPCt31KyaTHXXHMNqaOm0vGyO4lLSAZg54NXO05potmxY8eYNm0ay5YtIzMzkw9/lM9AK/5i1xzhKuDR2+B7z8LXb7+NgoIC/t//+3/2+7GFrAUwDJ138/1knT+cZcuW4UnrQpcZv6b9uOvsh9s4IR4vXb5yD+mXfAM8XorXvsO+v36PsvwNrqOZKHfkL8IVY9JZtmwZvTvDov+04s/4fPdKuP9GqFGYNWsWN910EydOnHAdK6JYARhGtm7dys0330zhP35OTelxkvqOpccdj5KYMcR1NBPjJM5DhwtuoMfXHyG+W1+qju7nwJxZHFn0jM3GM0Hx5ptvMuJeWLkd+nb1dfv26+Y6lQkn902Hef8J7du355VXXmH8+PFs2rTJdayIYQVgGNi/fz/f/va3GTJkCP/4xz8QbwLpF91Ot5t+jielg+t4xtRK6JpFj9t/R/uJNwNwfPnL9O3bl9/97necPHnScToTDQoLC5kxYwbXXnstBUdgfH9f8ZfV1XUyE46+MgZWrlzJsGHD2LRpE+PHj+fFF1+0GcLNIJH2j5STk6O5ubmuY7SJ/Px8ht18DydyX0cry0HiSB3xJTpc+FW8aV1cxzPmrMoLNnL4/dlUHPgCgLjk9rQfdz1pY65m9+9vdpzuTCKySlUjdsOwaDr3NURVee655/jRj37E4cOHSUlJ4Vf/UsL3vgwea6owTSgugzufgJdW+O6PHz+e+++/nyuvvDLmh081du6zAjDESkpKeO2113j22WdZuHBh7V8pyQMn0PGirxPfpZfjhMY0n6pSuj2XY5/OpWLfZgAksR3thkwh5fwLSeo9AonzAO4ni1gBGJ6Ki4uZM2cOf/qff2eNf+OZK4bDn++Evtbla1pAFWZ/AL94BQqP+x6bNGkS999/P5dddlnMFoJWADq0f/9+Fi1axJ2/epKSzZ+iFaW+JzxeUgZMIC3nWpIyh7oNaUwrqCplu9ZybNmLlO/Oq308LqUDKQMnkjL4QnY/ew+JiYnOMloBGF4+//xzZs+ezfPPP187eL9re3j4Frh9im/3B2POxckyeHwhPPQmHCr2PTauH9w083+ZPn06AwYMcBswxKwADJGKigq2bNnCpT/5G2V71lG2O4+qw/mnHZPQ43xSR1xOyuApeJLTHCU1JjgqCrdzctOnlGxeQtXhgtrHExISGD16NBMnTmTixIlMmDCBXr16heyvcisA3crPz2fRokV8/PHHLHrrSbbuP/XclPNh5uVww3hIjHeX0USXE6Xwx/fg4flwpM4Q5WHDhjF9+nSmTp3KqFGjSE1NdRcyBJwUgCIyFXgU8ABPquqD9Z4X//NXASXAHaq6+myv6fokWFNTQ2FhIQUFBRQUFLB3717u/dtHVB3Kp+LgbqqO7AWtOe1rJD6RxIyhJPUeQcrAidbNa2KCqlJZtJOTm5ZQunU5lQd3A6efb9LS0hgwYAADBw6svWRkZNC9e3fOO+88OnfuTFxc2wwAC1UBGIzzHrg/9zVHTU0NBw4cYOvWradd1qxZw/bt2087tkMK3HYh/PtlMNxOiSaIisvgnbXwWi689RkcLz31nIgwcOBAsrOzyc7OZujQoWRmZpKRkUG3bt3a7PzjUsgLQBHxAFuAK4B8YCVwi6puqHPMVcD38J0ILwAeVdULzva6LTkJfvHFF2zdupWqqiqqq6tPu66srKSioqL2Ul5eTklJyWmXkydPcvToUY4cOVJ7ffjI0TMKvHqfHG/H7sR37kViz8Ek9R5BQveBiMfW3Daxrab8JOV7t1C+dxPlBZuo2LeFmrKzr9vl8Xjo2rUrHTt2pEOHDrRv3772euzYsXzrW99q9vuHogAM1nkPWnbuW7NmDeXl5YCvEA+c51WVmpqaM66rq6tPuzR0jiwrK+PkyZOnXY4dO8bBgwcpKiqiqKiIQ4cOUVPT8PmxfbKvpe/iIb7LmCzwepr1cYxpMxVV8PEGXzG4dCtsKICq6oaP9Xq99OzZk/POO4/09HQ6dOhQe0lLSyMpKem0S2JiIl6v97SLx+MhLi7utIuINHoBTusVaaiHZODAgaSnpzf7Mzd27gtmVTIe2Kaq2/0B5gLXAXVXj70O+Jv6zk7LRSRdRHqo6r62CPD888/zi1/8oi1e6jRxye3xpHbCk9YZb2pnPGmdie+USXyXXng7ZhAX726ckzHhKi6xHcl9R5Pcd3TtY9Wlx6k6vJfKI3upOrKXyiP7qD55mJqTR6k+eYTqsmL279/P/v37z3i9gwcPtqgADBHn5z2Aq6++mr1797bVy7VI51QYcB4M7H7qMrgnjOxts3mNewleuHKk7wJQXgkbC2DNLvhsF2w7AAWHIf8wHCquYvfu3ezevdtt6Hrmz5/PVVdd1erXCWYBmAHsqXM/H99fu00dkwGcdiIUkbuBu/13i0Vkc9tGbZYuwEGAmtLj1JQep7Jop4MYQVP7+aJUtH8+iP7PeNrne+ONN1o6frBPmyc6U5ud98DJua/VP0OHin2XFV+0UaLWi8b/F/aZIkebf66rr27xigoNnvuCWQA2dGau39/cnGNQ1SeAJ9oi1LkSkdxIHkDeFPt8kS/aP2OEfL42O+9B6M99EfJv3CL2mSJDNH4mCO/PFcwG+Xyg7tDeTKB+n0RzjjHGmEhh5z1jTEQIZgG4EhgoIn1FJAGYAcyrd8w84HbxmQAca8txMMYYE2J23jPGRISgdQGrapWIfBd4F99yCE+r6noRmel/fjawAN9MuG34lkP4RrDytAGnXdAhYJ8v8kX7Zwz7zxcF572w/zc+B/aZIkM0fiYI488VcQtBG2OMMcaY1rFJ+cYYY4wxMcYKQGOMMcaYGGMFYCNE5CYRWS8iNSLS6BRuEZkqIptFZJuIzAplxtYQkU4i8r6IbPVfd2zkuJ0ikicia0QkvPehounvh3/g/f/5n/9cRMa4yHmumvH5LhGRY/7v1xoR+ZmLnOdKRJ4WkUIRWdfI8xH9/YsEIvKwiGzy//u+JiLprjO1VnPP55EgUn/nNKap//ORSER6ichHIrLR/3P3A9eZGmIFYOPWAf8CfNLYAf5tnx4DpgFDgVtEZGho4rXaLOADVR0IfOC/35hLVTU7XNcyCmjm92MaMNB/uRv4U0hDtkILft4W+79f2ap6f0hDtt4zwNSzPB+x378I8j4wXFVH4tvW7r8c52kLTZ7PI0GE/85pzDOc/f98JKoC/lNVhwATgO+E4/fJCsBGqOpGVW1q1f3abZ9UtQIIbPsUCa4DnvXffha43l2UNtOc70ftNlyquhxIF5EeoQ56jiL5561ZVPUT4PBZDonk719EUNX3VLXKf3c5vnUKI1ozz+eRIOrOAc34Px9xVHWfqq723z4BbMS3209YsQKwdRrb0ikSnBdYe8x/3a2R4xR4T0RW+belCmfN+X5E8vesudknishaEXlbRIaFJlrIRPL3LxJ9E3jbdQhTy37+I4yIZAGjgRWOo5whmFvBhT0RWQh0b+Cp/1bVN5rzEg08Fjbr6pzt87XgZSar6l4R6Qa8LyKb/H+xhaM23YYrDDUn+2qgj6oWi8hVwOv4ukujRSR//8JGc859IvLf+Lqy/h7KbOeqDc7nkcB+/iOIiKQCrwA/VNXjrvPUF9MFoKp+qZUvEdZbOp3t84nIARHpoar7/F1ohY28xl7/daGIvIavCyJcC8Bo34aryex1TzKqukBEHheRLqoaLZusR/L3L2w0de4Tka8D1wCXa4QsFtsG5/NIYD//EUJE4vEVf39X1Vdd52mIdQG3TnO2fQpX84Cv+29/HTjjL2QRaSciaYHbwJX4BlOHq2jfhqvJzyci3UVE/LfH4/s/fijkSYMnkr9/EUFEpgL3AteqaonrPOY0kfw7J2b4z8FPARtV9Xeu8zTGCsBGiMh0EckHJgLzReRd/+M9RWQB+LZ9AgLbPm0EXlLV9a4yt9CDwBUishW4wn//tM8HnAcsEZG1wD+B+ar6jpO0zdDY90NEZop/Ky5823Btx7cN11+AbzsJew6a+fluBNb5v2f/B8yIlBYcABF5AVgGnC8i+SJyZ7R8/yLIH4E0fEM+1ojIbNeBWqux83mkifDfOQ1q6P+860xtYDJwG3CZnFqS6yrXoeqzreCMMcYYY2KMtQAaY4wxxsQYKwCNMcYYY2KMFYDGGGOMMTHGCkBjjDHGmBhjBaAxxhhjTIyxAtAYY4wxJsZYAWiMMcYYE2P+f8NBZDL0gb+9AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# Initial conditions\n", "x0 = 0.0\n", "v0 = 0.5\n", "\n", "# Input parameters\n", "kBT = 0.25\n", "gamma = 0.01\n", "dt = 0.01\n", "t_max = 10000\n", "freq = 10\n", "k = (2 * np.pi * freq)**2 # Spring constant (harmonic oscillator)\n", "\n", "### Define Potential: Energy and Force\n", "def ho_en_force(x, k=k):\n", " energy = 0.5 * k * x**2\n", " force = -k * x\n", " return energy, force\n", "\n", "### Langevin dynamics (assuming you have this function correctly defined)\n", "# Should return arrays: times, pos, vel, KE, PE\n", "times, pos, vel, KE, PE = langevin_md_1d(x0, v0, dt, kBT, gamma, t_max, ho_en_force)\n", "\n", "### Plotting\n", "fig, ax = plt.subplots(nrows=1, ncols=4, figsize=(10, 4))\n", "\n", "bins = 50\n", "\n", "# Theoretical distributions\n", "def gaussian_x(x, k, kBT):\n", " return np.exp(-k*x**2/(2*kBT)) / np.sqrt(2*np.pi*kBT/k)\n", "\n", "def gaussian_v(v, kBT):\n", " return np.exp(-v**2/(2*kBT)) / np.sqrt(2*np.pi*kBT)\n", "\n", "# Plot histograms\n", "ax[0].hist(pos, bins=bins, density=True, alpha=0.6, color='skyblue')\n", "ax[1].hist(vel, bins=bins, density=True, alpha=0.6, color='salmon')\n", "\n", "# Plot theoretical curves\n", "x_grid = np.linspace(min(pos), max(pos), 300)\n", "v_grid = np.linspace(min(vel), max(vel), 300)\n", "\n", "ax[0].plot(x_grid, gaussian_x(x_grid, k, kBT), 'k-', lw=2, label='Theory')\n", "ax[0].set_xlabel('Position x')\n", "ax[0].set_ylabel('P(x)')\n", "ax[0].legend()\n", "\n", "\n", "ax[1].plot(v_grid, gaussian_v(v_grid, kBT), 'k-', lw=2, label='Theory')\n", "ax[1].set_xlabel('Velocity v')\n", "ax[1].set_ylabel('P(v)')\n", "ax[1].legend()\n", "\n", "ax[2].plot(pos[-1000:], vel[-1000:], alpha=0.3)\n", "ax[2].set_xlabel('Position x')\n", "ax[2].set_ylabel('Velocity v')\n", "ax[2].set_title('Phase Space')\n", "\n", "E_tot=KE+PE\n", "ax[3].plot(E_tot)\n", "ax[3].set_xlabel('Time')\n", "ax[3].set_ylabel('Energy')\n", "ax[3].set_title('Total Energy')\n", "\n", "fig.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Double well potential" ] }, { "cell_type": "code", "execution_count": 206, "metadata": { "tags": [ "hide-input" ] }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "5c8646d733e44926a421c0911e97fd39", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(FloatSlider(value=1.0, description='k', max=1.0, min=0.1), FloatSlider(value=3.0, descri…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def double_well(x, k=1, a=3):\n", " \n", " energy = 0.25*k*((x-a)**2) * ((x+a)**2)\n", " force = -k*x*(x-a)*(x+a)\n", " \n", " return energy, force\n", " \n", "x = np.linspace(-6,6,1000)\n", "energy, force = double_well(x) \n", "\n", "plt.plot(x, energy, '-o',lw=3)\n", "plt.plot(x, force, '-', lw=3, alpha=0.5)\n", "\n", "plt.ylim(-20,40)\n", "plt.grid(True)\n", "plt.legend(['$U(x)$', '$F=-\\partial_x U(x)$'], fontsize=15)" ] }, { "cell_type": "code", "execution_count": 207, "metadata": { "tags": [ "hide-input" ] }, "outputs": [], "source": [ "# Potential \n", "def double_well(x, k=1, a=3):\n", " \n", " energy = 0.25*k*((x-a)**2) * ((x+a)**2)\n", " force = -k*x*(x-a)*(x+a)\n", " \n", " return energy, force\n", "\n", "# Ininital conditions\n", "x = 0.1\n", "v = 0.5\n", "\n", "# Input parameters of simulation\n", "kBT = 5 # vary this\n", "gamma = 0.1 # vary this\n", "dt = 0.05\n", "t_max = 10000\n", "freq = 10\n", "\n", "#### Run the simulation\n", "times, pos, vel, KE, PE = langevin_md_1d(x, v, dt, kBT, gamma, t_max, double_well)\n", "\n", "#### Plotting \n", "fig, ax =plt.subplots(nrows=1, ncols=2, figsize=(13,5))\n", "\n", "x = np.linspace(min(pos), max(pos), 50)\n", "\n", "ax[0].plot(pos)\n", "ax[1].hist(pos, bins=50, density=True, alpha=0.5);\n", "\n", "v = np.linspace(min(vel), max(vel),50)\n", "\n", "ax[0].set_xlabel('t')\n", "ax[0].set_ylabel('x(t)')\n", "\n", "ax[1].set_xlabel('Computed P(x)')\n", "ax[1].set_ylabel('x')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Additional resoruces for learning MD\n", "\n", "- [Goran Wahnstrom's MD and MC lectures](http://fy.chalmers.se/~tfsgw/CompPhys/)\n", "- [Online course materials: **principles of modern molecular simulation methods**](https://sites.engineering.ucsb.edu/~shell/che210d/assignments.html)\n", "- [\"Understanding Molecular Simulation: From Algorithms to Applications\" a book by Daan Frankel](https://www.elsevier.com/books/understanding-molecular-simulation/frenkel/978-0-12-267351-1)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Problems\n", "\n", "#### 1D potential\n", "\n", "- In molecular simulations, we often want to sample from the **Boltzmann distribution** $P(x) \\propto e^{-V(x)/k_BT}$. \n", "- Langevin dynamics enables us to do this by modeling the effects of both **deterministic forces** (from the potential) and **random thermal fluctuations**.\n", "- Write a Python function using NumPy to simulate a 1D system under Langevin dynamics for the following potential:\n", "\n", "$$\n", "V(x) = Ax^2 + Bx^3 + Cx^4\n", "$$\n", "\n", "You should:\n", "\n", "1. Implement Langevin dynamics to evolve position and velocity over time.\n", "2. Allow inputs for:\n", " - Coefficients A, B, C\n", " - Initial position and velocity\n", " - Time step `dt`, number of steps\n", " - Temperature `kBT` and friction coefficient `gamma`\n", "3. Store all positions and plot a **normalized histogram** (probability density) after the simulation.\n", "4. Try the following setups:\n", " - **Double well**: A = -1, B = 0, C = 1\n", " - **Asymmetric well**: A = -1, B = -1, C = 1\n", " - Vary initial positions and explain what changes\n", "\n", "\n", "5. Compare your histograms with $e^{-V(x)/k_BT}$. Do they match?\n", "6. What happens if `gamma` is very small? What about very large?" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }