diff --git a/.travis.yml b/.travis.yml index f1d58963..1ea50de2 100644 --- a/.travis.yml +++ b/.travis.yml @@ -30,7 +30,7 @@ install: - conda update -q conda - conda info -a - if [[ "$PYTHON_VER" == "2.7" ]]; then - conda create -q -n p4env python=$PYTHON_VER ci-psi4 psi4 numpy=1.12 matplotlib jupyter -c psi4; + conda create -q -n p4env python=$PYTHON_VER ci-psi4 psi4 numpy=1.12 matplotlib jupyter scipy -c psi4; else conda create -q -n p4env python=$PYTHON_VER ci-psi4 psi4 numpy=1.13 matplotlib jupyter scipy -c psi4/label/dev -c psi4; fi diff --git a/Tutorials/12_MD/12a_basics.ipynb b/Tutorials/12_MD/12a_basics.ipynb new file mode 100644 index 00000000..dc7f6c93 --- /dev/null +++ b/Tutorials/12_MD/12a_basics.ipynb @@ -0,0 +1,792 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "\n", + "import matplotlib\n", + "import numpy as np\n", + "import itertools\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.animation as animation\n", + "from ipywidgets import interact, interactive, fixed, FloatSlider, widgets" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Overview\n", + "======\n", + "\n", + "This notebook introduces a number of concepts from molecular mechanics: Lennard-Jones potentials, Verlet integration algorithms and periodic boundary conditions. The internal consistency of SI units make life easier when implementing equations, but we'll use kJ/mol, nm, a.m.u. and ps units for energy, distance, mass and time, respectively. In these units, forces on particle $i$ at the resulting from $F_i(\\mathbf{r})=m_ia_i(\\mathbf{r})$ and $F_i(\\mathbf{r})=-\\nabla_i U(\\mathbf{r})$, where $U(\\mathbf{r})$ is the potential energy at configuration $\\mathbf{r}$, are consistent from the perspective of units.\n", + "\n", + "The Lennard-Jones Potential\n", + "-------------------------\n", + "Given the well depth $\\epsilon$ (units of kJ/mol) and a parameter $\\sigma$ that is related to the radius of an atom, the Lennard-Jones (LJ) potential for a pair of atoms separated by $r$ is given by\n", + "\n", + "\\begin{equation}\n", + "U_{LJ}(r) = 4\\epsilon\\left[\\left(\\frac{\\sigma}{r}\\right)^{12}-\\left(\\frac{\\sigma}{r}\\right)^6\\right].\n", + "\\end{equation}\n", + "\n", + "The following plot shows the LJ energy, as a function of distance." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def showpotential(ax, sigma=0.25, epsilon=0.05):\n", + " \"\"\" Plot the Lennard-Jones potential for a pair,\n", + " given the sigma and epsilon values. \"\"\"\n", + " npts = 200\n", + " sigmapair = 2*sigma # combine the sigmas on each atom\n", + " # Go up to 15A, using npts points\n", + " maxr = 1.5\n", + " rvals = np.linspace(0.015, maxr, npts)\n", + " sig_r = sigmapair / rvals\n", + " sig_r2 = np.square(sig_r)\n", + " sig_r6 = sig_r2*sig_r2*sig_r2\n", + " sig_r12 = sig_r6*sig_r6\n", + " energies = 4*epsilon*(sig_r12 - sig_r6)\n", + " ax.cla()\n", + " ax.plot(rvals,energies,lw=3)\n", + " ax.set_xlabel('R (nm)')\n", + " ax.set_ylabel('Energy (kJ/mol)')\n", + " ax.set_xlim(0, maxr)\n", + " ax.set_ylim(-0.4,0.2)\n", + " fig.canvas.draw()\n", + "\n", + "# Generate the plots; evaluate the next box to generate\n", + "# sliders that allow sigma and epsilon to be varied\n", + "fig,ax = plt.subplots(1,1,figsize=(4,4))\n", + "showpotential(ax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Evaluate this cell, then play with the sliders to see the effect of the well\n", + "# depth and radius on the shape of the potential curve for a dimer.\n", + "interact(showpotential,ax=fixed(ax),sigma=widgets.FloatSlider(min=0,max=0.5,step=.01,value=0.25),\n", + " epsilon=widgets.FloatSlider(min=0,max=.3,step=.01,value=.05),);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Periodic Boundary Conditions\n", + "--------------------------------\n", + "\n", + "Condensed phase simulations typically seek behavior of a given substance in the bulk. However, the presence of container walls can cause significantly different behavior from the bulk, so it's desirable to avoid confining substances with hard walls. One solution is to use periodic boundary conditions (PBC) whereby a unit cell is defined. A particle leaving a face of the unit cell will then simply reappear in the opposing wall with the same velocity and momentum. This is demonstrated below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rval = 0.5\n", + "dr = 0.02\n", + "boxlength = 1.5\n", + "\n", + "fig, ax = plt.subplots(1,1,figsize=(4,4))\n", + "line, = ax.plot(rval, 0, 'ko')\n", + "\n", + "def move_particle(i):\n", + " global rval\n", + " rval += dr\n", + " # Impose PBC\n", + " if(rval > boxlength):\n", + " rval -= boxlength\n", + " if(rval < 0):\n", + " rval += boxlength\n", + " line.set_xdata(rval)\n", + " return line,\n", + "\n", + "def init():\n", + " line.set_xdata(rval)\n", + " return line,\n", + "\n", + "ani = animation.FuncAnimation(fig, move_particle, 500, init_func=init,interval=25, blit=True, repeat=False)\n", + "ax.set_xlim(0, boxlength)\n", + "ax.get_xaxis().set_ticks([])\n", + "ax.get_yaxis().set_ticks([])\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Minimum Image Convention\n", + "-----------------------------\n", + "\n", + "Under PBC, we can think of the system as an infinite array of tessellated, identical unit cells; we are interested in studying the energy and properties of the primary unit cell in the presence of its cloned (\"image\") unit cells. Given the infinite nature of the system, we're obviously not going to compute an infinite number of interactions, so we use a cutoff to truncate the calculation to a range of pairs whose distance is shorter than some threshold.\n", + "\n", + "For each atom in the primary unit cell, we have to interact with all other atoms and images within the cutoff. To make this easy to implement, we settle on the *minimum image convention*: each primary atom $i$ will interact with only the closest version of atom $j$, whether the closest version is a primary or an image. This makes bookkeeping very simple, as we just have to shift interatomic distances to account for the PBC. However, this imposes the restriction that the cutoff must be less than half of the box length, or we will include multiple copies of $j$ for a given $i$, which would require extra bookkeeping to correctly implement.\n", + "\n", + "Below we setup a simple 2D box with particles on a regular grid, perturbed by a very small amount to break symmetry. We then implement periodic boundary conditions, and a function to compute the LJ energy and its gradient. The gradient expression is tested by numerical differentiation, which is always a good sanity check to perform." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "atoms_per_dim = 5\n", + "boxlength = 2*atoms_per_dim # atoms spaced by 2nm in each dimension\n", + "cutoff = boxlength/2 - 1e-8 # Longest cutoff possible\n", + "epsilon = 0.25 # kJ/mol\n", + "sigma = 0.2 # nm\n", + "\n", + "# intermediates just for convenience\n", + "cutoff2 = cutoff*cutoff\n", + "foureps = 4.0*epsilon\n", + "natoms = atoms_per_dim**2\n", + "pairs = np.array(list(itertools.combinations(range(natoms),2)))\n", + "iindices = pairs[:,0]\n", + "jindices = pairs[:,1]\n", + "sigmaij = 2.0*sigma # combine the sigma on each atom\n", + "sig2 = sigmaij*sigmaij\n", + "boxinv = 1.0/boxlength\n", + "\n", + "# Initial coordinates; regular lattice\n", + "spacing = np.linspace(-boxlength/2, boxlength/2, atoms_per_dim, endpoint=False)\n", + "spacing += boxlength/(2*atoms_per_dim)\n", + "# Center the array\n", + "coords = np.array([np.repeat(spacing, atoms_per_dim), np.tile(spacing, atoms_per_dim)]).T\n", + "# Move the atoms very slightly away from their regular grid positions\n", + "np.random.seed(0)\n", + "coords += (np.random.rand(natoms,2)-0.5)/50\n", + "\n", + "\n", + "def apply_periodicity(dR):\n", + " \"\"\" Apply minimum image convention. \"\"\"\n", + " dR -= boxlength*np.floor(dR*boxinv+0.5)\n", + " return dR\n", + "\n", + "def compute_energy_and_gradient_cutoff(coords):\n", + " \"\"\" Use a simple cutoff method to compute the LJ energy and gradient. \"\"\"\n", + " forces = np.zeros((natoms, 2))\n", + " # Apply minimum image convention\n", + " dR = apply_periodicity(coords[jindices] - coords[iindices])\n", + " r2 = np.einsum('ab,ab->a',dR, dR)\n", + " r2inv = 1.0/r2\n", + " # A cheesy way of applying the cutoffs\n", + " r2inv[r2>cutoff2] = 0.0\n", + " sig_r2 = sig2*r2inv\n", + " sig_r6 = sig_r2*sig_r2*sig_r2\n", + " sig_r12 = sig_r6*sig_r6\n", + " repulsive = foureps*sig_r12\n", + " attractive = foureps*sig_r6\n", + " pairenergies = repulsive - attractive\n", + " pairforces = np.einsum('ab,a->ab',dR,6*r2inv*(repulsive+pairenergies))\n", + " energy = pairenergies.sum()\n", + "\n", + " for n,(iind,jind) in enumerate(pairs):\n", + " forces[iind] += pairforces[n]\n", + " forces[jind] -= pairforces[n]\n", + " return energy, forces\n", + "#\n", + "# Test the gradient by finite differences\n", + "#\n", + "\n", + "# Displace the atoms in positive and negative directions for each\n", + "# degree of freedom and use the approximation\n", + "#\n", + "# E(r+delta) - E(r-delta)\n", + "# grad = -----------------------\n", + "# 2 delta\n", + "#\n", + "# to verify that the force expression has been coded up correctly.\n", + "fdgrad = np.zeros((natoms,2))\n", + "delta = 1e-6 # The step size can be varied\n", + "refE, refF = compute_energy_and_gradient_cutoff(coords)\n", + "for atom in range(natoms):\n", + " # Plus x\n", + " coords[atom,0] += delta\n", + " Epx, Fpx = compute_energy_and_gradient_cutoff(coords)\n", + " coords[atom,0] -= delta\n", + " # Minus x\n", + " coords[atom,0] -= delta\n", + " Emx, Fmx = compute_energy_and_gradient_cutoff(coords)\n", + " coords[atom,0] += delta\n", + " # Plus y\n", + " coords[atom,1] += delta\n", + " Epy, Fpy = compute_energy_and_gradient_cutoff(coords)\n", + " coords[atom,1] -= delta\n", + " # Minus y\n", + " coords[atom,1] -= delta\n", + " Emy, Fmy = compute_energy_and_gradient_cutoff(coords)\n", + " coords[atom,1] += delta\n", + " \n", + " fdgrad[atom] = [(Epx-Emx)/(2*delta), (Epy-Emy)/(2*delta)]\n", + "\n", + "print(\"Analytic:\")\n", + "print(refF)\n", + "print(\"Finite Difference:\")\n", + "print(fdgrad)\n", + "print(\"Difference:\")\n", + "print(refF-fdgrad)\n", + "print(\"Ratio:\")\n", + "print(refF/fdgrad)\n", + "\n", + "# Verify the energies and gradients are correct\n", + "assert np.allclose(-0.00370179210166, refE)\n", + "assert np.allclose(fdgrad, refF)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Running a trajectory\n", + "-------------------\n", + "\n", + "The aim of molecular dynamics is to propagate a system in time, according to Newton's equations of motion, in increments of $\\Delta$t. To do this, we use the beautifully simple Verlet integrator. Start by recognizing that the velocities are the time derivatives of the positions, $\\mathbf{v}(t) = \\mathbf{r}'(t)$, the accelerations are the second time derivatives, $\\mathbf{a}(t) = \\mathbf{r}''(t)$, and the jerk is the corresponding third derivative, $\\mathbf{b}(t) = \\mathbf{r}'''(t)$. Now we can expand the positions as a power series in time, in forwards and backwards directions:\n", + "\n", + "\\begin{equation}\n", + " \\begin{split}\n", + " \\mathbf{r}(t+\\Delta t) &= \\mathbf{r}(t) + \\mathbf{v}(t)\\Delta t + \\frac{1}{2}\\mathbf{a}(t)\\Delta t^2 + \\frac{1}{6}\\mathbf{b}(t)\\Delta t^3 + \\ldots \\\\\n", + " \\\\\n", + " \\mathbf{r}(t-\\Delta t) &= \\mathbf{r}(t) - \\mathbf{v}(t)\\Delta t + \\frac{1}{2}\\mathbf{a}(t)\\Delta t^2 - \\frac{1}{6}\\mathbf{b}(t)\\Delta t^3 + \\ldots \\\\\n", + " \\end{split}\n", + "\\end{equation}\n", + "\n", + "Adding these expressions and rearranging gives us the Verlet integration algorithm:\n", + "\n", + "\\begin{equation}\n", + " \\mathbf{r}(t+\\Delta t) = 2\\mathbf{r}(t) - \\mathbf{r}(t-\\Delta t) + \\mathbf{a}(t)\\Delta t^2 + \\mathcal{O}(\\Delta t^4).\n", + "\\end{equation}\n", + "\n", + "There are a few noteworthy features of this expression. First, we note that the velocities are not needed (they can be computed by finite differences of the current, previous and next positions) although variants of the algorithm do exist that use velocities. We need the previous step's position information, which means that we just use the simple forwards power series expansion above to get started. Note that the next and previous positions appear in a symmetrical way, which means that this algorithm can be run in reverse to recover the starting positions.\n", + "\n", + "We want the potential energy function, $U(\\mathbf{r})$, to drive the dynamics. To make that connection, we remember the result $F_i(\\mathbf{r})=m_i\\mathbf{a}_i(t)$ from classical mechanics, and rearrange \n", + "\n", + "\\begin{equation}\n", + " \\mathbf{a}_i(t) = \\frac{F_i(\\mathbf{r})}{m_i} = -\\frac{1}{m_i}\\frac{\\partial U(\\mathbf{r})}{\\partial \\mathbf{r}_i},\n", + "\\end{equation}\n", + "\n", + "remembering that the $i$ subscript labels each atom. Integrating the equations this way should maintain total energy, because no energy or particles are exchanged with the surroundings; in statistical mechanics this is known as a microcanonical ensemble, or NVE to reflect the that that number of particles, volume and energy are all constant. In the microcanonical ensemble, although total energy is conserved, kinetic and potential may interconvert through time.\n", + "\n", + "To test different simulation conditions easily, we encapsulate all of the simulation machinery we've developed so far into a class, which can take simulation setup as input to allow us to experiment with different conditions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Rgas = 0.0083144621 # Molar gas constant in kJ/(mol K)\n", + "\n", + "class MDSimulation(object):\n", + " \"\"\" A class to encapsulate coordinate and potential energy information about\n", + " a Lennard-Jones simulation in a cubic unit cell. Plotting functions are\n", + " provided to keep track of the simulation in real time, using Matplotlib. \"\"\"\n", + " def __init__(self, sigma, epsilon, mass, temp, dt, atoms_per_dim, boxlength, cutoff, nsteps):\n", + " \"\"\" Set up the simulation conditions and LJ potential information. \"\"\"\n", + " self.sigma = sigma\n", + " self.epsilon = epsilon\n", + " self.mass = mass\n", + " self.dt = dt\n", + " self.T = temp\n", + " self.RT = Rgas*self.T\n", + " self.atoms_per_dim = atoms_per_dim\n", + " self.boxlength = boxlength\n", + " self.cutoff = cutoff\n", + " self.nsteps = nsteps\n", + " \n", + " # intermediates just for convenience\n", + " self.cutoff2 = self.cutoff**2\n", + " self.foureps = 4.0*self.epsilon\n", + " self.natoms = self.atoms_per_dim**2\n", + " self.pairs = np.array(list(itertools.combinations(range(self.natoms),2)))\n", + " self.iindices = self.pairs[:,0]\n", + " self.jindices = self.pairs[:,1]\n", + " self.sigmaij = 2.0*self.sigma # combine the sigma on each atom\n", + " self.sig2 = self.sigmaij**2\n", + " self.boxinv = 1.0/self.boxlength\n", + " self.minv = 1.0/self.mass\n", + " \n", + " if(cutoff >= boxlength/2):\n", + " raise ValidationError(\"The cutoff must be less than boxlength/2\")\n", + "\n", + " np.random.seed(0)\n", + " self.init_coords()\n", + " self.last_coords = np.copy(self.coords)\n", + " # Generate velocities\n", + " velocities = np.sqrt(self.RT/self.mass)*np.random.normal(size=(self.natoms,2))\n", + " KE = 0.5*self.mass*np.sum(np.square(velocities))\n", + " # Take a step x1 = x0 + v dt + 1/2 a dt^2\n", + " PE, gradient = self.compute_energy_and_gradient(self.coords)\n", + " \n", + " self.timestep = 0.0\n", + " self.timeseries = [ self.timestep ]\n", + " self.PEs = [ PE ]\n", + " self.KEs = [ KE ]\n", + " self.TEs = [ KE + PE ]\n", + " a = -self.minv*gradient\n", + " self.coords += velocities*self.dt + 0.5*self.dt*self.dt*a\n", + " \n", + " # Plot objects\n", + " self.fig, (self.enerplot, self.coordplot) = plt.subplots(1,2,figsize=(8.2,4))\n", + " self.coorddata, = self.coordplot.plot(self.coords[:,0], self.coords[:,1], 'ko')\n", + " self.kline, = self.enerplot.plot(0, 0, 'bo', ms=0.5, label='kinetic')\n", + " self.vline, = self.enerplot.plot(0, 0, 'ro', ms=0.5, label='potential')\n", + " self.eline, = self.enerplot.plot(0, 0, 'ko', ms=0.5, label='total')\n", + " self.coordplot.set_xlim(-self.boxlength/2, self.boxlength/2)\n", + " self.coordplot.set_ylim(-self.boxlength/2, self.boxlength/2)\n", + " self.coordplot.set_aspect('equal')\n", + " self.coordplot.get_xaxis().set_ticks([])\n", + " self.coordplot.get_yaxis().set_ticks([])\n", + " self.coordplot.set_ylabel('simulation box')\n", + " self.enerplot.set_xlim(0, nsteps*self.dt)\n", + " self.enerplot.set_xlabel(\"Time (ps)\")\n", + " self.enerplot.set_ylabel(\"Energy (kJ/mol)\")\n", + " self.enerplot.legend()\n", + " self.animator = None\n", + "\n", + "\n", + " def compute_energy_and_gradient(self, coords):\n", + " \"\"\" Method to compute the LJ energy and gradient, which can be overloaded. \"\"\"\n", + " return self.compute_energy_and_gradient_cutoff(coords)\n", + " \n", + " def propagate(self):\n", + " \"\"\" Update the particle positions, using Verlet integration. \"\"\"\n", + " # Take a step r(t+dt) = 2 r(t) - r(t-dt) + dt^2 a(t)\n", + " PE, gradient = self.compute_energy_and_gradient(self.coords)\n", + " a = -self.minv*gradient\n", + " newcoords = 2*self.coords - self.last_coords + self.dt*self.dt*a\n", + " v = (newcoords - self.last_coords)/(2*self.dt)\n", + " KE = 0.5*self.mass*np.sum(np.square(v))\n", + " # Update coordinates\n", + " self.last_coords = self.coords.copy()\n", + " self.coords = newcoords.copy()\n", + " displaycoords = self.apply_periodicity(self.coords.copy())\n", + " self.coorddata.set_data(displaycoords[:,0], displaycoords[:,1]) # update the data\n", + " self.timestep += self.dt\n", + " self.timeseries.append(self.timestep)\n", + " self.PEs.append(PE)\n", + " self.KEs.append(KE)\n", + " self.TEs.append(PE + KE)\n", + " self.vline.set_data(self.timeseries, self.PEs)\n", + " self.kline.set_data(self.timeseries, self.KEs)\n", + " self.eline.set_data(self.timeseries, self.TEs)\n", + " en = (PE, KE, PE+KE)\n", + " minlim = np.min((1.2*np.min(en), self.enerplot.get_ylim()[0]))\n", + " maxlim = np.max((1.2*np.max(en), self.enerplot.get_ylim()[1]))\n", + " self.enerplot.set_ylim(minlim, maxlim)\n", + " return self.coorddata,\n", + " \n", + " def init_coords(self):\n", + " \"\"\" Set up the initial simulation coordinates on a regular grid, perturbed slightly. \"\"\"\n", + " # Initial coordinates; regular lattice\n", + " spacing = np.linspace(-self.boxlength/2, self.boxlength/2, self.atoms_per_dim, endpoint=False)\n", + " spacing += self.boxlength/(2*self.atoms_per_dim)\n", + " # Center the array\n", + " self.coords = np.array([np.repeat(spacing, self.atoms_per_dim), np.tile(spacing, self.atoms_per_dim)]).T\n", + " # Move the atoms very slightly away from their regular grid positions\n", + " self.coords += 0.2*(np.random.rand(self.natoms,2)-0.5)\n", + "\n", + " def testgrad(self, show=0):\n", + " \"\"\" Use finite differences to verify that the energy and gradient expressions are consistent. \"\"\"\n", + " fdgrad = np.zeros((self.natoms,2))\n", + " delta = 1e-6 # The step size can be varied\n", + " refE, refF = self.compute_energy_and_gradient(self.coords)\n", + " for atom in range(self.natoms):\n", + " # Plus x\n", + " self.coords[atom,0] += delta\n", + " Epx, Fpx = self.compute_energy_and_gradient(self.coords)\n", + " self.coords[atom,0] -= delta\n", + " # Minus x\n", + " self.coords[atom,0] -= delta\n", + " Emx, Fmx = self.compute_energy_and_gradient(self.coords)\n", + " self.coords[atom,0] += delta\n", + " # Plus y\n", + " self.coords[atom,1] += delta\n", + " Epy, Fpy = self.compute_energy_and_gradient(self.coords)\n", + " self.coords[atom,1] -= delta\n", + " # Minus y\n", + " self.coords[atom,1] -= delta\n", + " Emy, Fmy = self.compute_energy_and_gradient(self.coords)\n", + " self.coords[atom,1] += delta\n", + " #\n", + " fdgrad[atom] = [(Epx-Emx)/(2*delta), (Epy-Emy)/(2*delta)]\n", + " if show:\n", + " print(\"Analytic:\")\n", + " print(refF)\n", + " print(\"Finite Difference:\")\n", + " print(fdgrad)\n", + " print(\"Difference:\")\n", + " print(refF-fdgrad) \n", + " print(\"Ratio:\")\n", + " print(refF/fdgrad) \n", + " # Confirm the gradients are correct\n", + " assert np.allclose(refF, fdgrad)\n", + "\n", + " def apply_periodicity(self, dR):\n", + " \"\"\" Apply minimum image convention to an internuclear vector. \"\"\"\n", + " dR -= self.boxlength*np.floor(dR*self.boxinv+0.5)\n", + " return dR\n", + "\n", + " def compute_energy_and_gradient_cutoff(self, coords):\n", + " \"\"\" Compute the LJ energy and gradient, using a simple truncation scheme. \"\"\"\n", + " forces = np.zeros((self.natoms, 2))\n", + " # Apply minimum image convention\n", + " dR = self.apply_periodicity(coords[self.jindices] - coords[self.iindices])\n", + " r2 = np.einsum('ab,ab->a',dR, dR)\n", + " r2inv = 1.0/r2\n", + " # A cheesy way of applying the cutoffs\n", + " r2inv[r2>self.cutoff2] = 0.0\n", + " sig_r2 = self.sig2*r2inv\n", + " sig_r6 = sig_r2*sig_r2*sig_r2\n", + " sig_r12 = sig_r6*sig_r6\n", + " repulsive = self.foureps*sig_r12\n", + " attractive = self.foureps*sig_r6\n", + " pairenergies = repulsive - attractive\n", + " pairforces = np.einsum('ab,a->ab',dR,6*r2inv*(repulsive+pairenergies))\n", + " energy = pairenergies.sum()\n", + " for n,(iind,jind) in enumerate(self.pairs):\n", + " forces[iind] += pairforces[n]\n", + " forces[jind] -= pairforces[n]\n", + " return energy, forces\n", + " \n", + " def run_simulation(self):\n", + " \"\"\" A hook to propagate the trajectory and animate in real time using\n", + " Matplotlib's Funcanimation routine. \"\"\"\n", + " def take_step(i, self):\n", + " # A dummy function because FuncAnimation wants to pass an integer in\n", + " self.propagate()\n", + " # N.B. The animator must be persistent, which is why we cache it as a class variable here\n", + " self.animator = animation.FuncAnimation(self.fig, take_step, self.nsteps, fargs=(self,),\n", + " interval=25, repeat=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run a short simulation with a short cutoff\n", + "\n", + "This contrived example shows what can happen if insufficient cutoffs are used in a simulation. First, we run with a very short 0.5 nm cutoff." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim = MDSimulation(sigma=0.27, epsilon=0.3, mass=2, temp=300, dt=0.005, atoms_per_dim=5, boxlength=3,\n", + " cutoff=0.5, nsteps=400)\n", + "# Verify that the first energy is correct\n", + "assert np.allclose(13.7148290021, sim.compute_energy_and_gradient(sim.coords)[0])\n", + "# Verify that the first gradient is correct\n", + "sim.testgrad(show=0)\n", + "sim.run_simulation()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run a short simulation with a long cutoff\n", + "\n", + "Running with a longer cutoff restores energy conservation. Although kinetic and potential energy are interchanged, the total energy stays roughly constant with a 1.2 nm cutoff. However, the computational cost is increased, as more pairs must be considered at each time step." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim = MDSimulation(sigma=0.27, epsilon=0.3, mass=2, temp=300, dt=0.005, atoms_per_dim=5, boxlength=3, cutoff=1.2, nsteps=500)\n", + "# Verify that the first energy is correct\n", + "assert np.allclose(2.52645225873, sim.compute_energy_and_gradient(sim.coords)[0])\n", + "# Verify that the first gradient is correct\n", + "sim.testgrad(show=0)\n", + "sim.run_simulation()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Dealing with cutoff errors\n", + "--------------------------\n", + "The origin of the energy deviations is a step in the potential, at the cutoff. As atoms cross this discontinuity, energy is gained or lost, so the microcanonical ensemble is not properly preserved. The longer cutoffs remedied the situation, because the magnitude of the step is diminished. To visualize this, we plot the truncated LJ potential below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def showpotential(ax, sigma=0.25, epsilon=0.11, cutoff=0.8):\n", + " npts = 200\n", + " sigmapair = 2*sigma # combine the sigmas on each atom\n", + " # Go up to 1.5nm, using npts points\n", + " maxr = 1.5\n", + " rvals = np.linspace(0.3, maxr, npts)\n", + " sig_r = sigmapair / rvals\n", + " sig_r[rvals>cutoff] = 0\n", + " sig_r2 = np.square(sig_r)\n", + " sig_r6 = sig_r2*sig_r2*sig_r2\n", + " sig_r12 = sig_r6*sig_r6\n", + " energies = 4*epsilon*(sig_r12 - sig_r6)\n", + " ax.cla()\n", + " ax.plot(rvals,energies,lw=3)\n", + " ax.set_xlabel('R (nm)')\n", + " ax.set_ylabel('Energy (kJ/mol)')\n", + " ax.set_xlim(0, maxr)\n", + " ax.set_ylim(-0.4,0.2)\n", + " fig.canvas.draw()\n", + "\n", + "# Plot the initial LJ potential. Running the next box will\n", + "# provide slider to vary the different parameters.\n", + "fig,ax = plt.subplots(1,1,figsize=(3,3))\n", + "showpotential(ax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Evaluate this cell, then play with the sliders to see the effect of the well\n", + "# depth, radius and cutoff on the shape of the potential curve for a dimer.\n", + "interact(showpotential,ax=fixed(ax),sigma=widgets.FloatSlider(min=0,max=0.5,step=.01,value=0.25),\n", + " epsilon=widgets.FloatSlider(min=0,max=.3,step=.01,value=.11),\n", + " cutoff=widgets.FloatSlider(min=0,max=1.5,step=.05,value=0.8));" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Switching techniques\n", + "---------------------\n", + "\n", + "To avoid the 'brute force' remedy of using longer cutoffs, we can instead introduce a switching function to smooth the potential around the cutoff. To do this, we define a window immediately inside the cutoff and smoothly attenuate the potential from its value at the start of the window, to zero at the end of the window, which occurs at the cutoff. We now plot the switched LJ potential." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def showpotential(ax, sigma=0.25, epsilon=0.11, cutoff=0.8, window=0.2):\n", + " \"\"\" Plot a switched LJ potential, provided the parameters and the switching window. \"\"\"\n", + " npts = 200\n", + " sigmapair = 2*sigma # combine the sigmas on each atom\n", + " # Go up to 15A, using npts points\n", + " maxr = 1.5\n", + " rvals = np.linspace(0.3, maxr, npts)\n", + " ron2 = (cutoff-window)**2\n", + " roff2 = cutoff**2\n", + " r2vals = np.square(rvals)\n", + " switch = np.where(r2vals>ron2, (roff2 + 2*r2vals - 3*ron2)*(roff2-r2vals)**2/(roff2-ron2)**3, 1)\n", + " # Introduce cutoff\n", + " switch[rvals>cutoff] = 0\n", + " # Introduce switch\n", + " sig_r = (sigmapair/rvals)\n", + " sig_r2 = np.square(sig_r)\n", + " sig_r6 = sig_r2*sig_r2*sig_r2\n", + " sig_r12 = sig_r6*sig_r6\n", + " energies = 4*epsilon*switch*(sig_r12 - sig_r6)\n", + " ax.cla()\n", + " ax.plot(rvals,energies,lw=3)\n", + " ax.set_xlabel('R (nm)')\n", + " ax.set_ylabel('Energy (kJ/mol)')\n", + " ax.set_xlim(0, maxr)\n", + " ax.set_ylim(-0.4,0.2)\n", + " fig.canvas.draw()\n", + "\n", + "# Plot the switched LJ potential; run the next box to get\n", + "# sliders to vary the parameters.\n", + "fig,ax = plt.subplots(1,1,figsize=(3,3))\n", + "showpotential(ax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Evaluate this cell, then play with the sliders to see the effect of the well\n", + "# depth, radius, window and cutoff on the shape of the potential curve for a dimer.\n", + "interact(showpotential,ax=fixed(ax),sigma=widgets.FloatSlider(min=0,max=0.5,step=.01,value=0.25),\n", + " epsilon=widgets.FloatSlider(min=0,max=.3,step=.01,value=.11),\n", + " cutoff=widgets.FloatSlider(min=0,max=1.5,step=.05,value=0.8),\n", + " window=widgets.FloatSlider(min=0.001,max=0.4,step=.05,value=0.2));" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Updating the simulation class\n", + "------------------------------\n", + "\n", + "We can modify the simulation class by deriving a new type from it and overriding the `compute_energy_and_gradient` method, replacing it with the switched LJ function; all other functions will be the same as those in the parent `MDSimulation` class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class MDSimulationSwitch(MDSimulation):\n", + " \"\"\" Overloaded version of MDSimulation to use a switched LJ potential. \"\"\"\n", + " \n", + " def __init__(self, sigma, epsilon, mass, temp, dt, atoms_per_dim, boxlength, cutoff, window, nsteps):\n", + " \"\"\" Instantiate the underlying MDSimulation object, with additional switching window information. \"\"\"\n", + " self.window = window\n", + " super(MDSimulationSwitch, self).__init__(sigma, epsilon, mass, temp, dt, atoms_per_dim, boxlength, cutoff, nsteps)\n", + "\n", + " \n", + " def compute_energy_and_gradient(self, coords):\n", + " \"\"\" Compute the LJ energy and gradient, with optional switching to smooth the potential. \"\"\"\n", + " if self.window == 0:\n", + " return self.compute_energy_and_gradient_cutoff(coords)\n", + " else:\n", + " return self.compute_energy_and_gradient_switch(coords)\n", + " \n", + " def compute_energy_and_gradient_switch(self, coords):\n", + " \"\"\" Compute the switched LJ potential energy and gradient. \"\"\"\n", + " ron2 = (self.cutoff-self.window)**2\n", + " roff2 = self.cutoff2\n", + " denom = 1.0 / (roff2 - ron2)**3\n", + " forces = np.zeros((self.natoms, 2))\n", + " # Apply minimum image convention\n", + " dR = self.apply_periodicity(coords[self.jindices] - coords[self.iindices])\n", + " r2 = np.einsum('ab,ab->a',dR, dR)\n", + " r2inv = 1.0/r2\n", + " # Switching functions\n", + " switch = np.where(r2>ron2, (roff2 + 2.0*r2 - 3.0*ron2)*(roff2 - r2)**2 * denom, 1.0)\n", + " dswitch = np.where(r2>ron2, -12.0*(r2 - roff2)*(r2 - ron2) * denom, 0.0)\n", + " # Introduce cutoff\n", + " switch[r2>self.cutoff2] = 0.0\n", + " dswitch[r2>self.cutoff2] = 0.0\n", + " sig_r2 = self.sig2*r2inv\n", + " sig_r6 = sig_r2*sig_r2*sig_r2\n", + " sig_r12 = sig_r6*sig_r6\n", + " rep = self.foureps*sig_r12\n", + " att = self.foureps*sig_r6\n", + " pairenergies = switch*(rep - att)\n", + " pairforces = np.einsum('ab,a->ab', dR, 6.0*r2inv*switch*(2.0*rep-att) + dswitch*(rep-att))\n", + " energy = pairenergies.sum()\n", + " for n,(iind,jind) in enumerate(self.pairs):\n", + " forces[iind] += pairforces[n]\n", + " forces[jind] -= pairforces[n]\n", + " return energy, forces" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run the same short cutoff simulation, but with a switching function\n", + "Hint: to run the switch-free simulation again, set the window to zero." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim = MDSimulationSwitch(sigma=0.27, epsilon=0.3, mass=2, temp=300, dt=0.005, atoms_per_dim=5,\n", + " boxlength=3, cutoff=0.5, nsteps=500, window=0.1)\n", + "# Verify that the first energy is correct\n", + "assert np.allclose(5.76688248814, sim.compute_energy_and_gradient(sim.coords)[0])\n", + "# Verify that the first gradient is correct\n", + "sim.testgrad(show=0)\n", + "sim.run_simulation()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "By introducing switching, energy is conserved even with a 0.5 nm cutoff.\n", + "\n", + "References\n", + "=====\n", + "\n", + "1. M. P. Allen and D. J. Tildesley, *Computer Simulation of Liquids*, Oxford University Press, New York 1987.\n", + "2. D. Frenkel and B. Smit, *Understanding Molecular Simulation: From Algorithms to Applications*, Academic Press, San Diego 1996.\n", + "3. P. J. Steinbach and B. R. Brooks, *J. Comp. Chem.*, 15 667 (1994)." + ] + } + ], + "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.6.2" + }, + "widgets": { + "state": { + "2bc6e6bc0a2b46238a57008b002c11ea": { + "views": [ + { + "cell_index": 18 + } + ] + }, + "7025377924dc4ba2a440de1981fc6d3c": { + "views": [ + { + "cell_index": 21 + } + ] + }, + "a3295e3afd834bfc829e694165a948bd": { + "views": [ + { + "cell_index": 3 + } + ] + } + }, + "version": "1.2.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Tutorials/12_Ewald_Summation/12a_ewald.ipynb b/Tutorials/12_MD/12b_ewald.ipynb similarity index 97% rename from Tutorials/12_Ewald_Summation/12a_ewald.ipynb rename to Tutorials/12_MD/12b_ewald.ipynb index 32a4c018..b87655de 100644 --- a/Tutorials/12_Ewald_Summation/12a_ewald.ipynb +++ b/Tutorials/12_MD/12b_ewald.ipynb @@ -12,7 +12,9 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "collapsed": true + }, "outputs": [], "source": [ "%matplotlib notebook\n", @@ -168,6 +170,7 @@ "# Center the array\n", "coords = np.array(list(itertools.product(spacing, spacing, spacing)))\n", "# Move the atoms very slightly away from their regular grid positions\n", + "np.random.seed(0)\n", "coords += 0.2*(np.random.rand(natoms,3)-0.5)\n", "coords[coords>boxlength/2] = boxlength/2\n", "coords[coords<-boxlength/2] = -boxlength/2\n", @@ -262,7 +265,17 @@ "ax.set_xlabel('Number of image box layers')\n", "ax.set_ylabel('Energy (kJ/mol)')\n", "ax.legend(loc='lower right')\n", - "plt.show()" + "plt.show()\n", + "\n", + "# Check the reference values\n", + "assert np.allclose(vals,\n", + "[[ 361515.2359571, 361515.2359571 ],\n", + " [ 528282.46725449, 557057.25818972],\n", + " [ 534335.79047581, 536496.90616012],\n", + " [ 535962.70789396, 536005.76078745],\n", + " [ 536633.65606731, 537475.71261986],\n", + " [ 536973.60561882, 537518.58787992],\n", + " [ 537169.21808044, 537527.68088663]])" ] }, { @@ -331,8 +344,7 @@ "# Draw the plots with a default alpha value. Run the box below\n", "# to get a slider to vary alpha and update the plots in real time.\n", "fig, (gausplot, erfplot) = plt.subplots(1,2,figsize=(8,4))\n", - "plot_gaussians(gausplot, erfplot)\n", - "\n" + "plot_gaussians(gausplot, erfplot)" ] }, { @@ -540,6 +552,15 @@ " [5, 1, 3],\n", " [6, 1, 4],\n", "])\n", + "params = np.array([ \n", + " [0, 4, 2],\n", + " [1, 3, 3],\n", + " [2, 2, 4],\n", + " [3, 2, 4],\n", + " [4, 1, 5],\n", + " [5, 1, 5],\n", + " [6, 1, 6],\n", + "])\n", "\n", "# Generate Ewald energies, varying the attentuation paramter alpha,\n", "# as well as the real and recpiprocal space cutoffs needed to reach\n", @@ -560,7 +581,10 @@ "ax.set_ylabel('Energy (kJ/mol)')\n", "ax.legend()\n", "\n", - "plt.show()" + "plt.show()\n", + "\n", + "# Confirm that all high alpha values have the same total energy\n", + "assert np.allclose(Uvals[1,3], Uvals[2:,3])" ] }, { @@ -692,7 +716,7 @@ " return Udir,Urec,Uslf,U\n", "\n", "# Repeat the same scan performed above for Ewald, using PME with 4th order spline.\n", - "efunc = functools.partial(PME_energy, splineorder=4)\n", + "efunc = functools.partial(PME_energy, splineorder=5)\n", "Uvals = np.array([efunc(*p) for p in params])\n", "\n", "# Plot the PME energy and its components as a function of attenuation parameter.\n", @@ -708,7 +732,10 @@ "ax.set_ylabel('Energy (kJ/mol)')\n", "ax.legend()\n", "\n", - "plt.show()" + "plt.show()\n", + "\n", + "# Confirm that all high alpha values have the same total energy\n", + "assert np.allclose(Uvals[1,3], Uvals[2:,3], atol=5)" ] }, { diff --git a/Tutorials/12_Ewald_Summation/README.md b/Tutorials/12_MD/README.md similarity index 100% rename from Tutorials/12_Ewald_Summation/README.md rename to Tutorials/12_MD/README.md diff --git a/tests/test_TU_12.py b/tests/test_TU_12.py index 9e0180b3..3d014ad3 100644 --- a/tests/test_TU_12.py +++ b/tests/test_TU_12.py @@ -2,10 +2,14 @@ from utils import * -tdir = 'Tutorials/12_Ewald_Summation' - +tdir = 'Tutorials/12_MD' @using_matplotlib def test_12a(workspace): - exe_scriptified_ipynb(workspace, tdir, '12a_ewald') + exe_scriptified_ipynb(workspace, tdir, '12a_basics') + +@using_matplotlib +@using_scipy +def test_12b(workspace): + exe_scriptified_ipynb(workspace, tdir, '12b_ewald')