diff --git a/.travis.yml b/.travis.yml index ca84d7d4..7828f45d 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 pylibefp -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..a1c6cea4 --- /dev/null +++ b/Tutorials/12_MD/12a_basics.ipynb @@ -0,0 +1,820 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "outputs": [], + "source": [ + "atoms_per_dim = 5\n", + "boxlength = 2.0*atoms_per_dim # atoms spaced by 2nm in each dimension\n", + "cutoff = boxlength/2.0 - 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)*0.02\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": { + "collapsed": true + }, + "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.0):\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.0, self.boxlength/2.0)\n", + " self.coordplot.set_ylim(-self.boxlength/2.0, self.boxlength/2.0)\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.0, self.boxlength/2.0, self.atoms_per_dim, endpoint=False)\n", + " spacing += self.boxlength/(2.0*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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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": { + "collapsed": true + }, + "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_MD/12b_ewald.ipynb b/Tutorials/12_MD/12b_ewald.ipynb new file mode 100644 index 00000000..5dab578f --- /dev/null +++ b/Tutorials/12_MD/12b_ewald.ipynb @@ -0,0 +1,796 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Ewald summation\n", + "\n", + "This notebook provides a detailed look at classical electrostatics, demonstrating problems related to convergence and how to overcome those issues. The techniques discussed are ubiquitous in periodic simulations, both classical and quantum." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "from mpl_toolkits.mplot3d import Axes3D\n", + "import matplotlib\n", + "import numpy as np\n", + "import functools\n", + "import scipy.special\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": [ + "## Classical electrostatics\n", + "\n", + "Molecular simulations of periodic systems and crystals are an extremely important area of study. To demonstrate the difficulties, and some popular methods to overcome them, we'll investigate a simple system of unit charges. Central to our discussion is the potential $\\phi$ at a point $\\mathbf{R}$ generated by a Gaussian charge density centered at $\\mathbf{R}_0$ with exponent $\\alpha$\n", + "\n", + "\\begin{equation}\n", + " q G(\\alpha, \\mathbf{R}_0, \\mathbf{R}) \\equiv q \\left(\\frac{\\alpha}{\\sqrt{\\pi}}\\right)^3\n", + " e^{-\\alpha^2(\\mathbf{R}-\\mathbf{R}_0)^2},\n", + "\\end{equation}\n", + "\n", + "which can be obtained as the solution of Poisson's equation:\n", + "\n", + "\\begin{equation}\n", + " \\nabla^2 \\phi(\\alpha, \\mathbf{R}) = -\\frac{q G(\\alpha, \\mathbf{0}, \\mathbf{R})}{\\epsilon_0}.\n", + "\\end{equation}\n", + "\n", + "The above involves the permittivity of free space $\\epsilon_0$ and we have assumed the Gaussian is placed at the origin. The Laplacian in spherical coordinates is\n", + "\n", + "\\begin{align}\n", + " \\nabla^2 & = \\frac{1}{r^2}\\frac{\\partial}{\\partial r}\\left(r^2\\frac{\\partial}{\\partial r}\\right)\n", + " +\\frac{1}{r^2 \\sin^2\\phi} \\frac{\\partial^2}{\\partial \\theta^2}\n", + " +\\frac{1}{r^2 \\sin \\theta}\\frac{\\partial}{\\partial \\phi}\\left(\\sin \\phi \\frac{\\partial}{\\partial \\phi}\\right),\n", + "\\end{align}\n", + "\n", + "but only the leading term will survive when applied to a spherically symmetric density. Recognizing that $\\frac{1}{r^2}\\frac{\\partial}{\\partial r}\\left(r^2\\frac{\\partial}{\\partial r}\\right)f = \\frac{1}{r}\\frac{\\partial^2}{\\partial r^2}(rf)$ Poisson's equation for a Gaussian potential is\n", + "\n", + "\\begin{align}\n", + " \\frac{1}{r}\\frac{\\partial^2}{\\partial r^2}\\left[r \\phi(\\alpha, \\mathbf{R})\\right] \n", + " &= -\\frac{q G(\\alpha, \\mathbf{0}, \\mathbf{R})}{\\epsilon_0} \\\\\n", + " \\frac{\\partial^2}{\\partial r^2}\\left[r \\phi(\\alpha, \\mathbf{R})\\right] \n", + " &= -\\frac{r q G(\\alpha, \\mathbf{0}, \\mathbf{R})}{\\epsilon_0} \\\\\n", + " \\frac{\\partial}{\\partial r}\\left[r \\phi(\\alpha, \\mathbf{R})\\right] \n", + " &= -\\int_r^\\infty \\frac{r q G(\\alpha, \\mathbf{0}, \\mathbf{R})}{\\epsilon_0}\\mathrm{d}r \n", + " = \\frac{q G(\\alpha, \\mathbf{0}, \\mathbf{R})}{2\\alpha^2\\epsilon_0} \\\\\n", + " r \\phi(\\alpha, \\mathbf{R}) \n", + " &= \\int_0^r \\frac{q G(\\alpha, \\mathbf{0}, \\mathbf{R})}{2\\alpha^2\\epsilon_0}\\mathrm{d}r \\\\\n", + " \\phi(\\alpha, \\mathbf{R}) \n", + " &= \\frac{q \\mathrm{Erf}(\\alpha r)}{4 \\pi \\epsilon_0 r},\\\\\n", + "\\end{align}\n", + "\n", + "which introduces the definition of the error function:\n", + "\\begin{equation}\n", + " \\mathrm{Erf}(z) \\equiv \\frac{2}{\\sqrt{\\pi}}\\int_0^z{-t^2}\\mathrm{d}t\n", + "\\end{equation}\n", + "\n", + "and the scalar distance $r=|\\mathbf{R}-\\mathbf{R}_0|$.\n", + "Recognizing that the infinite exponent limit of a Gaussian is simply a delta function, we get a very familiar result\n", + "\n", + "\\begin{equation}\n", + " \\lim_{\\alpha\\rightarrow\\infty}\\frac{q \\mathrm{Erf}(\\alpha r)}{4 \\pi \\epsilon_0r} = \\frac{q}{4 \\pi \\epsilon_0 r}\n", + "\\end{equation}\n", + "\n", + "*i.e.* the potential at a distance $r$ from a point charge $q$ is $\\phi(r) = \\frac{q}{4 \\pi \\epsilon_0 r}$. Placing a test charge $Q$ that same distance away then gets us to Coulomb's law:\n", + "\n", + "\\begin{equation}\n", + " U = Q\\phi(r) = \\frac{Qq}{4 \\pi \\epsilon_0 r}.\n", + "\\end{equation}\n", + "\n", + "If we want to embed a unit cell in an infinite number of replica (\"image\") unit cells in all directions, we just need to sum over all interactions within the primary unit cell, and over all primary-image cell interactions. The unit cell is described by a set of vectors, which for the orthorhombic case with lengths $a,b,c$ are simply \n", + "\n", + "$$\n", + " \\mathbf{a} = \\left(\\begin{eqnarray} &a&0&0\\\\ &0&b&0\\\\ &0&0&c\\\\ \\end{eqnarray}\\right).\n", + "$$\n", + "\n", + "With this notation, and labeling atoms $\\{i,j\\}$ possessing charges $\\{q_i,q_j\\}$, image cells can be compactly denoted by a vector $\\mathbf{n} = \\{n_1 \\mathbf{a}_1, n_2 \\mathbf{a}_2, n_3 \\mathbf{a}_3\\}$, where $n_1, n_2$ and $n_3$ are integers:\n", + "\n", + "\\begin{equation}\n", + " U = \\frac{1}{8 \\pi \\epsilon_0}\\sum_{\\mathbf{n}'}\\sum_{i,j}\\frac{q_i q_j}\n", + " {\\left|\\mathbf{r}_i - \\mathbf{r}_j + \\mathbf{n}\\right|} \n", + " = \\frac{1}{8 \\pi \\epsilon_0}\\sum_{\\mathbf{n}'}\\sum_{i,j}\\frac{q_i q_j}{r_{ij}}\n", + "\\end{equation}\n", + "\n", + "The prime on the summation of over image cells is a reminder that $\\mathbf{n}=\\mathbf{0}$ terms where $i=j$ should be omited and we have introduced the shorthand notation $r_{ij}=\\left|\\mathbf{r}_i - \\mathbf{r}_j + \\mathbf{n}\\right|$. Implementing electrostatics by brute force summation over many lattice vectors is not a good idea; to see why we'll take a brief diversion.\n", + "\n", + "## Conditional convergence\n", + "To get a feel for conditional convergence consider the classic example of the alternating harmonic series:\n", + "\n", + "\\begin{equation}\n", + " \\sum_{n=0}^\\infty \\frac{(-1)^{n+1}}{n} = \\left(1 - \\frac{1}{2} + \\frac{1}{3}\n", + " - \\frac{1}{4} + \\frac{1}{5} - \\frac{1}{6} + \\frac{1}{7} - \\frac{1}{8} + \\ldots\\right) = \\ln(2).\n", + "\\end{equation}\n", + "\n", + "Summing the terms in ascending order yields the result show. However, we can also reorder the terms in parentheses before summing\n", + "\n", + "\\begin{equation}\n", + " \\underbrace{1 - \\frac{1}{2}}_{\\frac{1}{2}} -\n", + " \\underbrace{\\frac{1}{4}}_{\\frac{1}{4} } +\n", + " \\underbrace{\\frac{1}{3} - \\frac{1}{6} }_{\\frac{1}{6}} -\n", + " \\underbrace{\\frac{1}{8}}_{\\frac{1}{8} } +\n", + " \\underbrace{\\frac{1}{5} - \\frac{1}{10} }_{\\frac{1}{10}} -\n", + " \\underbrace{\\frac{1}{12}}_{\\frac{1}{12} } +\n", + " \\underbrace{\\frac{1}{7} - \\frac{1}{14} }_{\\frac{1}{14}} -\n", + " \\underbrace{\\frac{1}{16}}_{\\frac{1}{16} } +\n", + " \\ldots,\n", + "\\end{equation}\n", + "\n", + "which yields exactly half of what we got above:\n", + "\n", + "\\begin{equation}\n", + " \\frac{1}{2}\\left(1 - \\frac{1}{2} + \\frac{1}{3} - \\frac{1}{4} + \\frac{1}{5} - \\frac{1}{6}\n", + " + \\frac{1}{7} - \\frac{1}{8} + \\ldots\\right) = \\frac{1}{2} \\ln(2).\n", + "\\end{equation}\n", + "\n", + "There exists some specific permutation of the summands that will generate any given result when the sum is conditionally convergent. A sum over some series $a_n$ is conditionally convergent if $\\sum_{n=1}^\\infty |a_n| = \\infty$. Returning to electrostatics, we introduce the integral analog of the potential expression above, avoiding the singularity by neglecting the infinitesimal point $\\mathbf{r}=\\mathbf{0}$:\n", + "$$\\int_{r>0} \\frac{1}{|\\mathbf{r}|} \\mathrm{d}^3\\mathbf{r} = 4\\pi \\int_{\\delta}^\\infty \\frac{r^2}{r} \\mathrm{d}r = \\infty$$\n", + "The factor of $r^2$ comes from the conversion to polar coordinates. We find that\n", + "\n", + "\\begin{equation}\n", + " \\int_{\\delta}^\\infty \\frac{1}{r^{n-2}} \\mathrm{d}r = \n", + " \\begin{cases}\n", + " \\frac{r_c^{3-n}}{n-3}, & \\text{if}\\ n>3 \\\\\n", + " \\infty, & \\text{otherwise}\n", + " \\end{cases}\n", + "\\end{equation}\n", + "\n", + "so, by analogy to the integrals, the Coulomb summation is conditionally convergent in three dimensional space. We don't arrive at an *absolutely* convergent series until we deal with an operator that decays as $\\frac{1}{r^4}$ in three dimensions. Not only is the Coulombic summation conditionally convergent, its convergence is very slow w.r.t. increasing cutoffs due to the slow decay of the operator.\n", + "\n", + "The purpose of this notebook is to demonstrate a clever trick to decompose the slow, conditionally convergent summation into two rapid, absolutely convergent terms. We use kJ/mol for energies, nm for length and atomic units for charge throughout." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Setup up a cubic simulation box\n", + "atoms_per_dim = 5\n", + "natoms = atoms_per_dim**3\n", + "epsilon0 = 5.72765E-4 # (e^2 mol) / (kJ nm)\n", + "coulomb = 1/(4*np.pi*epsilon0)\n", + "boxlength = 0.8 #nm\n", + "# Initial coordinates; regular lattice\n", + "spacing = np.linspace(-boxlength/2.0, boxlength/2.0, atoms_per_dim, endpoint=False)\n", + "spacing += boxlength/(2.0*atoms_per_dim)\n", + "# 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.0] = boxlength/2.0\n", + "coords[coords<-boxlength/2.0] = -boxlength/2.0\n", + "charges = np.zeros(len(coords))\n", + "if 1:\n", + " # Generate a large dipole\n", + " charges[:natoms//2] = 1.0\n", + " charges[(natoms+1)//2:] = -1.0\n", + "else:\n", + " # Generate a more homogeneous distribution\n", + " charges[1:natoms:2] = -1.0\n", + " charges[0:natoms-1:2] = 1.0\n", + "poscrd = coords[charges > 0]\n", + "negcrd = coords[charges < 0]\n", + "neucrd = coords[charges == 0]\n", + "print(\"Sum of charges\", np.sum(charges))\n", + "# Visualize the coordinates and charges in a rotatable 3D plot. Atoms are colored according to their charge.\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111, projection='3d')\n", + "ax.scatter(poscrd[:,0], poscrd[:,1], poscrd[:,2], 'bo')\n", + "ax.scatter(negcrd[:,0], negcrd[:,1], negcrd[:,2], 'ro')\n", + "ax.scatter(neucrd[:,0], neucrd[:,1], neucrd[:,2], 'ko')\n", + "ax.set_xlim(-boxlength/2.0, boxlength/2.0)\n", + "ax.set_ylim(-boxlength/2.0, boxlength/2.0)\n", + "ax.set_zlim(-boxlength/2.0, boxlength/2.0)\n", + "ax.set_xticks([])\n", + "ax.set_yticks([])\n", + "ax.set_zticks([])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the lattice energy using brute force summation, and evaluate over a range of lattice cutoffs.\n", + "\n", + "def coulomb_energy(qij, R):\n", + " \"\"\" Evaluate the conventional Coulomb energy, given products of charges and\n", + " distances between the centers bearing those charges. \"\"\"\n", + " return np.sum(qij/R)\n", + "\n", + "\n", + "def lattice_sum(charges, coords, nboxes, efun):\n", + " \"\"\" Compute the Coulomb energy of a system of charges, by considering interactions\n", + " between the primary box and all image boxes that are within nboxes of the\n", + " primary cell. The energy is evaluated by the provided function, efun. \"\"\"\n", + " natoms = coords.shape[0]\n", + " boxlist = range(-nboxes,nboxes+1)\n", + "\n", + " # Build some lists to describe interactions \n", + " # products of charges between two boxes\n", + " chargeij = np.prod(list(itertools.product(charges, charges)), axis=1)\n", + " dRij = np.repeat(coords, natoms, axis=0) - np.tile(coords, (natoms,1))\n", + " energy = 0.0\n", + " truncated_energy = 0.0\n", + " for nvecs in itertools.product(boxlist, boxlist, boxlist):\n", + " nvecs = np.array(nvecs)\n", + " shiftvec = nvecs*boxlength\n", + " absvecs = np.abs(nvecs)\n", + " if nvecs.any():\n", + " # This is an image box\n", + " shiftedRij = np.linalg.norm(dRij+shiftvec, axis=1)\n", + " eterm = 0.5*efun(qij=chargeij, R=shiftedRij)\n", + " else:\n", + " # Handle the 'home' box\n", + " eterm = 0.0\n", + " for i in range(natoms):\n", + " for j in range(i):\n", + " dR = np.linalg.norm(coords[i]-coords[j])\n", + " eterm += efun(qij=charges[i]*charges[j], R=dR)\n", + " eterm *= coulomb\n", + " if np.linalg.norm(absvecs) <= nboxes:\n", + " truncated_energy += eterm\n", + " energy += eterm \n", + " return energy,truncated_energy\n", + " \n", + "vals = []\n", + "for n in range(7):\n", + " vals.append(lattice_sum(charges, coords, n, efun=coulomb_energy))\n", + "vals = np.array(vals)\n", + "\n", + "# Plot the convergence of the energy w.r.t. cutoff using the two truncation schemes\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111)\n", + "xvals = range(1,len(vals))\n", + "ax.plot(xvals, vals[1:,0], 'bo', label ='Cubic')\n", + "ax.plot(xvals, vals[1:,1], 'ro', label = 'Spherical')\n", + "ax.set_xticks(xvals)\n", + "ax.set_xlabel('Number of image box layers')\n", + "ax.set_ylabel('Energy (kJ/mol)')\n", + "ax.legend(loc='lower right')\n", + "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]])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above code sets up a cell and performs the summation in two different ways 1) adding a complete layer of image boxes, which corresponds to a cubic macroscopic crystal and 2) adding image boxes in such a way that the macroscopic system is roughly spherical. The convergence is both slow and conditional; for this contrived example of a large net dipole the energies for the two summation schemes converge to slightly different results. When the second layer of vectors is added, we are introducing interaction up to three times the box length and the results are still not fully converged!\n", + "\n", + "To speed up the convergence, Ewald summation introduces Gaussian functions at the atomic positions, with the opposite sign to a given atom's charge, to screen the point charges, making them short range. To compensate, a Gaussian with the same sign as the charge is added in, but is handled separately. From the above arguments, this affects the potential by splitting the Coulomb operator into short- and long-range components:\n", + "\n", + "\\begin{align}\n", + " \\frac{1}{r} &= \\frac{1-\\mathrm{Erf}(\\alpha r)}{r} + \\frac{\\mathrm{Erf}(\\alpha r)}{r}\\\\\n", + " &= \\frac{\\mathrm{Erfc}(\\alpha r)}{r} + \\frac{\\mathrm{Erf}(\\alpha r)}{r}\n", + "\\end{align}\n", + "\n", + "In the plot below we show the density (left) where point charges and their screening Gaussians are shown in blue; the resulting $\\mathrm{Erfc}$ potential is shown in the right panel. Increasing the attenuation parameter $\\alpha$, which is inversely proportional to the width of the Gaussians, screens this short range potential more effectively, but also increases the overall contribution of the long-range potential (red) that results from the compensating Gaussians. The rapid decay of the short-range term means we can evaluate it in the standard real-space treatment using pairwise loops with short cutoffs. The red part is still long-range so will still decay slowly; however it is now free from the singularity that plagues the full (black) potential, and we can develop a treatment for it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define some point charges to demonstrate the density and potential\n", + "positions = [1, 3, 5]\n", + "qvals = [1, -1, 1]\n", + "\n", + "def plot_gaussians(gausplot, erfplot, alpha=3):\n", + " \"\"\" Generate a plot that shows point charges in a 1D coordinates system, whose\n", + " density is represented by blue delta functions, with blue screening Gaussians\n", + " introduced to attenuate these point charges. The Gaussian is also subtracted\n", + " out, which is shown in red, to make the method exact. The potential corresponding\n", + " to the regular Coulomb operator (black) is also shown, decomposed into its short\n", + " range (blue) component and its long range (red) term. \"\"\"\n", + " prefac = alpha/np.sqrt(np.pi) # 1D Gaussian normalization\n", + " def plot_gaussian(axis, xvals, origin, exponent, prefac, col):\n", + " axis.plot(xvals, prefac*np.exp(-np.square(exponent*(xvals-origin))), color=col)\n", + " npts = 150\n", + " xvals = np.linspace(0.001, 6, npts)\n", + " gausplot.cla()\n", + " erfplot.cla()\n", + " for p,c in zip(positions, qvals):\n", + " plot_gaussian(gausplot, xvals, p, alpha, c*prefac, 'r')\n", + " plot_gaussian(gausplot, xvals, p, alpha, -c*prefac, 'b')\n", + " gausplot.plot([p, p], [0, c], color='b')\n", + " gausplot.plot([0, 40], [0, 0], color='k', linestyle='-', linewidth=2)\n", + " erfvals = scipy.special.erf(alpha*xvals)\n", + " inv = 1/xvals\n", + " erfplot.set_title('Potential')\n", + " erfplot.plot(xvals, erfvals*inv, color='r', label='Erf(alpha R) / R')\n", + " erfplot.plot(xvals, (1-erfvals)*inv, color='b', label='Erfc(alpha R) / R')\n", + " erfplot.plot(xvals, inv, '--', color='k', label='1 / R')\n", + " erfplot.set_xlabel('r (nm)')\n", + " erfplot.set_ylabel('$\\\\phi(r)$')\n", + " erfplot.legend()\n", + " gausplot.set_title('Density')\n", + " gausplot.set_xlim(0, 6)\n", + " gausplot.set_ylim(-5,5)\n", + " gausplot.set_ylabel('$\\\\rho(r)$')\n", + " gausplot.set_xlabel('Position')\n", + " erfplot.set_xlim(0, 1.5)\n", + " erfplot.set_ylim(0, 10)\n", + " fig.canvas.draw()\n", + "\n", + "# 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)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Evaluate this cell, then play with the sliders to see the effect of\n", + "# varying alpha on the long and short range densities and potentials.\n", + "interact(plot_gaussians, gausplot=fixed(gausplot), erfplot=fixed(erfplot),\n", + " alpha=widgets.FloatSlider(min=0,max=9,step=.2,value=3),);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looking at the form of the long-range density in the left panel of the above plot, it is a periodic, smoothly oscillating function so any solution is likely to involve trigonometric functions somehow. Returning to the discussion at the top of the notebook, we remember that the potential comes from solving the Poisson equation:\n", + "\n", + "\\begin{equation}\n", + " \\nabla^2 \\phi(\\alpha, \\mathbf{R}) = -\\frac{q G(\\alpha, \\mathbf{0}, \\mathbf{R})}{\\epsilon_0}.\n", + "\\end{equation}\n", + "\n", + "Combining these considerations, it's easy to imagine that the trigonometric functions will be introduced via complex exponentials; because $\\frac{\\partial}{\\partial x}e^{-imx} = -im e^{-imx}$, introducing an exponential will effectively turn the derivative operator in the Poisson equation into a constant, making it easy to handle. As a reminder, the connection between complex exponentials and trigonometric functions is provided by Euler's formula: $e^{ix} = \\mathrm{cos}(x)+i\\mathrm{sin}(x)$. The transformation that will introduce the complex exponential is the Fourier transform.\n", + "\n", + "First, we normalize the coordinates into the range $[0,1]$, appropriate for integration, by introducing the reciprocal lattice vectors $\\mathbf{a}^*$, which are defined in terms of the unit cell vectors by the orthonormality relation $\\mathbf{a}_i^*\\cdot\\mathbf{a}_j=\\delta_{ij}$ and for the orthorhombic unit cell are\n", + "\n", + "$$ \\mathbf{a}^* = \\left(\\begin{eqnarray} &\\frac1a&0&0\\\\ &0&\\frac1b&0\\\\ &0&0&\\frac1c\\\\ \\end{eqnarray}\\right). $$\n", + "\n", + "The resulting reciprocal space coordinate system has orthogonal axes, even for non-orthorhombic crystals. Analogous to $\\mathbf{n}$ indexing summation over image unit cells, we use $\\mathbf{m}$ to index summation over *reciprocal* space lattice vectors, $\\mathbf{m} = \\{m_1 \\mathbf{a}_1^*, m_2 \\mathbf{a}_2^*, m_3 \\mathbf{a}_3^*\\}$; because we're evaluating the energy using these vectors, we call the long range energy the *reciprocal* energy, $U_\\mathbf{rec}$.\n", + "\n", + "The long-range density $\\rho^L(\\mathbf{r})$ and potential $\\phi^L(\\mathbf{r})$ have Fourier transforms defined by\n", + "\n", + "\\begin{align}\n", + " \\tilde\\rho^L(\\mathbf{m}) &= \\int_V \\rho^L(\\mathbf{r}) e^{-2\\pi i\\mathbf{m}\\cdot\\mathbf{r}} \\mathrm{d}^3\\mathbf{r}\\\\\n", + " \\tilde\\phi^L(\\mathbf{m}) &= \\int_V \\phi^L(\\mathbf{r}) e^{-2\\pi i\\mathbf{m}\\cdot\\mathbf{r}} \\mathrm{d}^3\\mathbf{r}\n", + "\\end{align}\n", + "\n", + "and the inverse transforms are defined by summation over reciprocal lattice vectors $\\mathbf{m}$\n", + "\n", + "\\begin{align}\n", + " \\rho^L(\\mathbf{r}) &= \\frac1V \\sum_{\\mathbf{m}}\\tilde\\rho^L(\\mathbf{m}) e^{2\\pi i\\mathbf{r}\\cdot\\mathbf{m}}\\\\\n", + " \\phi^L(\\mathbf{r}) &= \\frac1V \\sum_{\\mathbf{m}}\\tilde\\phi^L(\\mathbf{m}) e^{2\\pi i\\mathbf{r}\\cdot\\mathbf{m}}\n", + "\\end{align}\n", + "\n", + "Inserting $\\tilde\\rho^L(\\mathbf{m})$ into Poisson's equation yields the Fourier space expression\n", + "\n", + "\\begin{equation}\n", + " 4\\pi^2m^2 \\tilde\\phi^L(\\mathbf{m}) = \\frac{\\tilde\\rho^L(\\mathbf{m})}{\\epsilon_0}\n", + "\\end{equation}\n", + "\n", + "whose solution is trivial if we neglect the $\\mathbf{m}=\\mathbf{0}$ term, which is zero for neutral unit cells with $\\sum_i q_i = 0$. We define the scalar $m=|\\mathbf{m}|$. The reciprocal space long-range density is obtained from the Fourier transform:\n", + "\n", + "\\begin{align}\n", + " \\tilde\\rho(\\mathbf{m}) &= \\int_V \\sum_\\mathbf{n}\\sum_j q_j \n", + " G(\\alpha, \\mathbf{r}_j, \\mathbf{r+\\mathbf{n}}) \n", + " e^{-2\\pi i \\mathbf{r}\\cdot\\mathbf{m}}\\mathrm{d}\\mathbf{r}^3 \\\\\n", + " &=\\sum_\\mathbf{n}\\sum_j q_j \\int_\\mathrm{all\\ space} \n", + " G(\\alpha, \\mathbf{r}_j, \\mathbf{\\mathbf{w}}) \n", + " e^{-2\\pi i \\mathbf{(w-n)}\\cdot\\mathbf{m}}\\mathrm{d}\\mathbf{w}^3 \\\\\n", + " &=\\sum_\\mathbf{n}\\sum_j q_j \\int_\\mathrm{all\\ space}\n", + " G(\\alpha, \\mathbf{r}_j, \\mathbf{\\mathbf{w}})\n", + " e^{-2\\pi i \\mathbf{w}\\cdot\\mathbf{m}}\\mathrm{d}\\mathbf{w}^3 \\\\\n", + " &= \\sum_j q_j e^{-\\frac{m^2\\pi^2}{\\alpha^2}}\n", + " e^{-2\\pi i \\mathbf{r}_j\\cdot\\mathbf{m}}.\n", + "\\end{align}\n", + "\n", + "In developing this expression we first change variables $\\mathbf{w}=\\mathbf{r+n}$ before taking advantage of the fact that $e^{2\\pi i \\mathbf{n}\\cdot\\mathbf{m}} = 1$. The integral on the third line is the definition of the Fourier transform of a 3D Gaussian, which is evaluated for each dimension using the standard result $\\int e^{a^2x^2}e^{-2\\pi i m x}\\mathrm{d}x = \\frac{\\sqrt{\\pi}}{a}e^{-\\frac{\\pi^2m^2}{a^2}}$. With the Fourier space expression for the density in hand, we can now get the Fourier space potential from the simplified Poisson equation:\n", + "\n", + "\\begin{equation}\n", + " \\tilde\\phi(\\mathbf{m}) = \\frac{1}{4\\pi^2\\epsilon_0}\n", + " \\frac{e^{-\\frac{m^2\\pi^2}{\\alpha^2}}}{m^2} \\sum_j q_j e^{-2\\pi i \\mathbf{r}_j\\cdot\\mathbf{m}}\n", + "\\end{equation}\n", + "\n", + "This is inverse Fourier transformed to obtain the real space potential at some point $\\mathbf{r}$:\n", + "\n", + "\\begin{equation}\n", + " \\phi(\\mathbf{r}) = \\frac{1}{4\\pi^2\\epsilon_0 V}\n", + " \\sum_{\\mathbf{m}\\ne\\mathbf{0}}\\frac{e^{-\\frac{m^2\\pi^2}{\\alpha^2}}}{m^2}\n", + " \\sum_j q_j e^{-2\\pi i \\mathbf{r}_j\\cdot\\mathbf{m}}\n", + " e^{2\\pi i \\mathbf{r}\\cdot\\mathbf{m}}\n", + "\\end{equation}\n", + "\n", + "And finally, as above, we evaluate the potential at each atom, multiply by the charge of that atom, and divide by 2 to avoid double-counting:\n", + "\n", + "\\begin{align}\n", + " U_\\mathrm{rec} &= \\frac{1}{8\\pi^2\\epsilon_0 V}\n", + " \\sum_{\\mathbf{m}\\ne\\mathbf{0}}\\frac{e^{-\\frac{m^2\\pi^2}{\\alpha^2}}}{m^2}\n", + " \\sum_{i} q_i e^{2\\pi i \\mathbf{r}_i\\cdot\\mathbf{m}}\n", + " \\sum_{j} q_j e^{-2\\pi i \\mathbf{r}_j\\cdot\\mathbf{m}} \\\\\n", + " &= \\frac{1}{8\\pi^2\\epsilon_0 V}\\sum_{\\mathbf{m}\\ne\\mathbf{0}}\n", + " \\frac{e^{-\\frac{m^2\\pi^2}{\\alpha^2}}}{m^2}\n", + " S(\\mathbf{m})\n", + " S(-\\mathbf{m})\n", + "\\end{align}\n", + "\n", + "On the surface, this doesn't seem like a simplification. However, the presence of the factor $\\frac{e^{-\\frac{m^2\\pi^2}{\\alpha^2}}}{m^2}$ leads to a rapidly convergent summation involving a limited number of reciprocal lattice vectors $\\mathbf{m}$. The quantity $S(\\mathbf{m})$ is known as a structure factor. The only remaining issue to deal with is the fact that the real space term neglects the $i=j$ term in the primary unit cell summation, but the reciprocal space energy contains this term as a result of the reciprocal space potential containing contributions from all atoms. The spurious term that we need to remove is the interaction of the probe charge with its own contribution to the reciprocal space potential, which comes from its screening Gaussian:\n", + "\n", + "\\begin{split}\n", + " U_\\mathrm{self} &= -\\sum_i \\lim_{r\\rightarrow0} \\frac{q_i \\mathrm{Erf}(\\alpha r)q_i}{4 \\pi \\epsilon_0 r} \\\\\n", + " &= -\\frac{1}{4 \\pi \\epsilon_0}\\frac{\\alpha}{\\sqrt\\pi}\\sum_i q_i^2\n", + "\\end{split}\n", + "\n", + "In the code below, we investigate the effect of introducing the screening to the short-range *direct space* (as opposed to reciprocal space, used to evaluate the long-range term) energy:\n", + "\n", + "\\begin{equation}\n", + " U = \\frac{1}{2}\\sum_{\\mathbf{n}'}\\sum_{i,j}\n", + " \\frac{q_i q_j\\mathrm{Erfc}(\\alpha\\left|\\mathbf{r}_i - \\mathbf{r}_j + \\mathbf{n}\\right|)}\n", + " {\\left|\\mathbf{r}_i - \\mathbf{r}_j + \\mathbf{n}\\right|}.\n", + "\\end{equation}\n", + "\n", + "After choosing an attenuation parameter $\\alpha=4 \\mathrm{nm}^{-1}$ we repeat the lattice summation shown above:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def ewald_dir_energy(alpha, maxn):\n", + " \"\"\" An energy function, to be passed to the lattice sum code, that\n", + " evaluated the direct contribution to the Ewald energy. \"\"\"\n", + " # Define the direct space energy as the attenuated direct Ewald energy.\n", + " def ewald_dir_interaction(qij, R, alpha):\n", + " \"\"\" Evaluate the energy of pairs of charges, separated by R,\n", + " using the attenuated Coulomb operator. \"\"\"\n", + " return np.sum(qij*scipy.special.erfc(alpha*R)/R)\n", + " enerfunc = functools.partial(ewald_dir_interaction, alpha=alpha)\n", + " return lattice_sum(charges, coords, n, efun=enerfunc)\n", + " \n", + "vals = np.array([ewald_dir_energy(alpha=4, maxn=n) for n in range(4)])\n", + "\n", + "# Plot the convergence of the lattice sum, with respect to number of\n", + "# image cells included, using cubic and spherical layering schemes.\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111)\n", + "xvals = range(1,len(vals))\n", + "ax.plot(xvals, vals[1:,0], 'bo', label ='Cubic')\n", + "ax.plot(xvals, vals[1:,1], 'ro', label = 'Spherical')\n", + "ax.set_xticks(xvals)\n", + "ax.set_xlabel('Number of image box layers')\n", + "ax.set_ylabel('Energy (kJ/mol)')\n", + "ax.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now both summation schemes converge rapidly to the same results, showing that the direct space part of the lattice sum is well-behaved. Implementing the reciprocal and self terms, we find that the overall Ewald energy is also well-behaved:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def ewald_rec_energy(alpha, maxm):\n", + " \"\"\" Evaluate the long range (reciprocal space) Ewald energy using an\n", + " attenuation parameter of alpha, and using up to maxm wavevectors. \"\"\"\n", + " mdim = 2*maxm+1\n", + " twopii = np.pi*2j\n", + " pi2_a2 = np.pi**2/alpha**2 if alpha != 0.0 else 1e50\n", + " miter = [ 0 ]\n", + " for m in range(1,maxm+1):\n", + " miter.extend([m,-m])\n", + " miter = np.array(miter)\n", + " mvals = miter / boxlength\n", + " xvecs = np.exp(twopii*np.outer(coords[:,0], mvals))\n", + " yvecs = np.exp(twopii*np.outer(coords[:,1], mvals))\n", + " zvecs = np.exp(twopii*np.outer(coords[:,2], mvals))\n", + " Sm = np.einsum('n,nx,ny,nz->xyz', charges, xvecs, yvecs, zvecs)\n", + " m = np.array(list(itertools.product(mvals, mvals, mvals)))\n", + " m2 = np.einsum('nx,nx->n', m, m).reshape(mdim,mdim,mdim)\n", + " prefac = np.exp(-pi2_a2*m2)\n", + " m2[0,0,0] = 1e50 # avoid division by 0\n", + " prefac /= m2\n", + " V = boxlength**3\n", + " return coulomb*np.einsum('xyz,xyz', prefac, np.square(Sm.real) + np.square(Sm.imag)) / (2*np.pi*V)\n", + "\n", + "def ewald_self_energy(alpha):\n", + " \"\"\" Evaluate the Ewald self interaction energy. \"\"\"\n", + " return -coulomb*(alpha/np.sqrt(np.pi))*np.sum(np.square(charges))\n", + " \n", + "def ewald_energy(alpha, maxn, maxm):\n", + " \"\"\" Evaluate all components of the Ewald sum, to get the total energy. \"\"\"\n", + " Udir = ewald_dir_energy(alpha, maxn)[0]\n", + " Urec = ewald_rec_energy(alpha, maxm)\n", + " Uslf = ewald_self_energy(alpha)\n", + " U = Udir+Urec+Uslf\n", + " return Udir,Urec,Uslf,U\n", + "\n", + "# A list of alpha, maxn, maxm to loop over\n", + "params = np.array([ \n", + " [0, 4, 2],\n", + " [1, 3, 2],\n", + " [2, 2, 2],\n", + " [3, 2, 2],\n", + " [4, 1, 3],\n", + " [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", + "# convergence for each alpha value.\n", + "Uvals = np.array([ewald_energy(*p) for p in params])\n", + "\n", + "# Plot the total energy, and its components, as a function of the\n", + "# attenuation parameter alpha.\n", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111)\n", + "alphavals = params[:,0]\n", + "ax.plot(alphavals, Uvals[:,0], '--r', label = 'Udir')\n", + "ax.plot(alphavals, Uvals[:,1], '--b', label = 'Urec')\n", + "ax.plot(alphavals, Uvals[:,2], '--g', label = 'Uslf')\n", + "ax.plot(alphavals, Uvals[:,3], 'k', label = 'U')\n", + "ax.set_xticks(alphavals)\n", + "ax.set_xlabel('Alpha (1/nm)')\n", + "ax.set_ylabel('Energy (kJ/mol)')\n", + "ax.legend()\n", + "\n", + "plt.show()\n", + "\n", + "# Confirm that all high alpha values have the same total energy\n", + "assert np.allclose(Uvals[2,3], Uvals[3:,3])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For almost any nonzero screening parameter, the total energy is the same, although its composition in terms of reciprocal, direct and self terms changes; this useful fact is very handy when testing the correctness of an Ewald implementation. For small $\\alpha$ more $\\mathbf{n}$ vectors are needed, as the short-range operator is not well attenuated, while a large $\\alpha$ demands more $\\mathbf{m}$ vectors to evaluate the reciprocal space term; this is evident from the values chosen in the `params` array above. Because cutoffs can be used for the direct space term, it scales as $\\mathcal{O}(N)$. Being attenuated, the direct space term can be evaluated within the primary unit cell, using the minimum image convention.\n", + "\n", + "Constructing the stucture factor $S(\\mathbf{m})$ involves a summation over $N$ atoms as well as $\\mathbf{m}$; for a given $\\alpha$ value the latter scales with the volume of the box and is therefore proportional to $N$. Overall, the Ewald reciprocal energy therefore scales as $\\mathcal{O}(N^2)$, although it can be shown that optimally tuning $\\alpha$ to minimize computational cost yields a $\\mathcal{O}(N^{{3/2}})$ algorithm. A clever trick introduced by Darden, York and Pedersen in 1993 reduces this cost further, as we will now demonstrate.\n", + "\n", + "The key insight is that a discrete Fourier transform can be performed using a fast recursion algorithm if the atoms are arranged on a uniform mesh, hence the name particle mesh Ewald (PME). A method that only works for selected configurations would be useless, so PME uses spline interpolation to effectively coerce the atoms to predefined, regularly spaced locations. Consider some function $f$ we've tabulated at every integer. If we want either $f(12)$ or $f(13)$, we're in luck and can just look up the tabulated values. However, $f(12.5)$ would come from $\\frac{f(12)+f(13)}{2}$, while $f(12.1)$ would be best approximated with two points by $0.9f(12) + 0.1f(13)$. This simple order 2 example can be generalized to order $n$, with the coefficients of each pre-tabulated data point provided by an order $n$ *spline* function. The cardinal basis spline (B-spline) used in PME is defined for interpolation using the nearerst $n$ points by the recursion\n", + "\n", + "\\begin{equation}\n", + " M_n(u) = \\frac{u}{n-1}M_{n-1}(u) + \\frac{n-u}{n-1} M_{n-1}(u-1)\n", + "\\end{equation}\n", + "\n", + "where $u$ is the fractional distance from the nearest grid point to the left of the point of interest. Using this method to discretize the density to a regular lattice the fast Fourier transform (FFT) can then be used to perform the discrete Fourier transform step, and scales as $\\mathcal{O}(N log(N))$ instead of $\\mathcal{O}(N^2)$; the same is true of the inverse tranformation step. The PME algorithm is otherwise similar to conventional Ewald method, as demonstrated below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def PME_rec_energy(alpha, maxm, splineorder):\n", + " \"\"\" Evaluate the particle mesh Ewald long range (reciprocal space)\n", + " energy via discretization onto a regular cubic grid. \"\"\"\n", + " \n", + " def build_spline(order, u):\n", + " \"\"\" Construct a B spline of given order, from a seed value u \"\"\"\n", + " def one_pass_bspline(array, val, n):\n", + " div = 1.0 / float(n - 1)\n", + " array[n-1] = div * val * array[n-2]\n", + " for j in range(1, n-1): \n", + " array[n-j-1] = div * ((val+j) * array[n-j-2] + (n-j-val) * array[n-j-1])\n", + " array[0] *= div * (1 - val)\n", + " return array\n", + " \n", + " array = np.zeros(order)\n", + " array[0] = 1.0 - u\n", + " array[1] = u\n", + " for m_n in range(1, order-1):\n", + " array = one_pass_bspline(array, u, m_n+2)\n", + " return array\n", + " \n", + " def dftmod(splineorder, nfft):\n", + " \"\"\" Compute the norm of a B spline in Fourier space. \"\"\"\n", + " ftarr = np.zeros(nfft)\n", + " # modulus is invariant to position, so just make a spline with 0\n", + " ftarr[:splineorder] = build_spline(splineorder, 0.0)\n", + " ftarr = np.fft.fft(ftarr)\n", + " return np.square(ftarr.real) + np.square(ftarr.imag)\n", + "\n", + " mdim = 2*maxm+1\n", + " # A convenient array to allow iteration over the grid from any given\n", + " # starting point, correctly wrapping around the boundaries.\n", + " griditer = np.array([ [ ((p+n)%mdim) for p in range(splineorder) ] for n in range(mdim) ])\n", + " # Define a simple iterator array to bring m into the range -0.5 < m <= 0.5\n", + " shiftm = np.array([ k - mdim if k >= (mdim+1)/2 else k for k in range(mdim) ])\n", + "\n", + " splines = np.zeros((natoms, 3, splineorder))\n", + " gridstart = np.zeros((natoms, 3), dtype=int)\n", + " # Define coordinates in fractional coordinates\n", + " reduced_coords = coords/boxlength + 0.5\n", + " for atom,coord in enumerate(reduced_coords):\n", + " # Define gridpoint to the left in each dimension, small shift\n", + " # to handle atoms located on the rightmost grid point\n", + " gridstart[atom,:] = list(map(int, mdim*coord-1e-8))\n", + " wvals = mdim*coord - gridstart[atom]\n", + " for xyz in range(3):\n", + " splines[atom,xyz,:] = build_spline(splineorder, wvals[xyz])\n", + " \n", + " # Form the influence function, including spline normalization terms. Note\n", + " # that this is independent of particle positions and can therefore be formed\n", + " # at the start of the simulation and cached for reuse in each time step.\n", + " pi2_a2 = np.pi**2/alpha**2 if alpha != 0.0 else 1e50\n", + " mvals = shiftm / float(boxlength)\n", + " splinemod = 1/dftmod(splineorder, mdim)\n", + " splinemods = np.einsum('x,y,z->xyz', splinemod, splinemod, splinemod)\n", + " m = np.array(list(itertools.product(mvals, mvals, mvals)))\n", + " m2 = np.einsum('nx,nx->n', m, m).reshape(mdim,mdim,mdim)\n", + " gF = np.multiply(splinemods, np.exp(-pi2_a2*m2))\n", + " m2[0,0,0] = 1e50 # avoid division by 0\n", + " gF /= m2\n", + " \n", + " # Spread the charges onto the charge grid, using spline coefficients at O(N) cost\n", + " Qgrid = np.zeros((mdim, mdim, mdim))\n", + " for atom in range(natoms):\n", + " q = charges[atom]\n", + " for nx,xgrid in enumerate(griditer[gridstart[atom,0]]):\n", + " xspline = splines[atom,0,nx]\n", + " for ny,ygrid in enumerate(griditer[gridstart[atom,1]]):\n", + " yspline = splines[atom,1,ny]\n", + " for nz,zgrid in enumerate(griditer[gridstart[atom,2]]):\n", + " zspline = splines[atom,2,nz]\n", + " Qgrid[xgrid,ygrid,zgrid] += q*xspline*yspline*zspline\n", + " \n", + " # Forward FFT Q to form S(m), at O(Nlog(N)) cost\n", + " Sm = np.fft.fftn(Qgrid)\n", + " \n", + " # Convolution of gF with S(m), at O(N) cost\n", + " conv = np.multiply(gF, Sm)\n", + " \n", + " # Inverse FFT (S(m)*gF) real potential grid, at O(Nlog(N)) cost\n", + " phigrid = mdim**3*np.fft.ifftn(conv).real\n", + "\n", + " # Probe the potential grid using the splines, to get the energy at O(N) cost\n", + " e = 0.0\n", + " for atom in range(natoms):\n", + " q = charges[atom]\n", + " for nx,xgrid in enumerate(griditer[gridstart[atom,0]]):\n", + " xspline = splines[atom,0,nx]\n", + " for ny,ygrid in enumerate(griditer[gridstart[atom,1]]):\n", + " yspline = splines[atom,1,ny]\n", + " for nz,zgrid in enumerate(griditer[gridstart[atom,2]]):\n", + " zspline = splines[atom,2,nz]\n", + " e += q*xspline*yspline*zspline*phigrid[xgrid,ygrid,zgrid]\n", + " \n", + " return e*coulomb / (2*np.pi*boxlength**3)\n", + "\n", + " \n", + "def PME_energy(alpha, maxn, maxm, splineorder):\n", + " \"\"\" Compute the total PME energy by evaluating all components. \"\"\"\n", + " Udir = ewald_dir_energy(alpha, maxn)[0]\n", + " Urec = PME_rec_energy(alpha, maxm, splineorder)\n", + " Uslf = ewald_self_energy(alpha)\n", + " U = Udir+Urec+Uslf\n", + " 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=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", + "fig = plt.figure()\n", + "ax = fig.add_subplot(111)\n", + "alphavals = params[:,0]\n", + "ax.plot(alphavals, Uvals[:,0], '--r', label = 'Udir')\n", + "ax.plot(alphavals, Uvals[:,1], '--b', label = 'Urec')\n", + "ax.plot(alphavals, Uvals[:,2], '--g', label = 'Uslf')\n", + "ax.plot(alphavals, Uvals[:,3], 'k', label = 'U')\n", + "ax.set_xticks(alphavals)\n", + "ax.set_xlabel('Alpha (1/nm)')\n", + "ax.set_ylabel('Energy (kJ/mol)')\n", + "ax.legend()\n", + "\n", + "plt.show()\n", + "\n", + "# Confirm that all high alpha values have the same total energy\n", + "assert np.allclose(Uvals[2,3], Uvals[3:,3], atol=5)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "The PME energy performs similarly to regular Ewald summation, with significantly better asymptotic scaling. Forces may be evaluated easily by replacing the splines used to probe the potential grid with their derivatives, which are defined by a very simple recursion relation, as outlined by Darden and co-workers.\n", + "\n", + "In developing the absolutely convergent sums above, we seem to have magically removed the conditional convergence. The reason for this is somewhat complicated and beyond the scope of this discussion. The conditional convergence results from the mutual interaction of the image cells' dipole moments. By thinking of the periodic system not as an infinite system, but a large macroscopic crystal with a given shape (we considered cubes and spheres in the lattice sum code above) is is possible to derive a \"correction\" to the Ewald energy that describes the conditionally (and therefore crystal shape) dependent part of the lattice sum. Such corrections have been derived for various crystal shapes, and they all depend on the cell dipole. In dynamic simulations under periodic boundary condisions, a charge particle exiting the top of a unit cell will reappear at the bottom, possibly flipping the dipole of the cell. If the extra shape-dependent term were considered in this circumstance, the potential energy would be discontinuous, so these corrections are usually neglected.\n", + "\n", + "\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. U. Essmann, L. Perera, M. L. Berkowitz, T. A. Darden, H. Lee, and L. G. Pedersen, J. Chem. Phys. 103, 8577 (1995).\n", + "4. C. Sagui, L. G. Pedersen, and T. A. Darden, J. Chem. Phys. 120, 73 (2004)." + ] + } + ], + "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": { + "4b125575e3d641f392dce4b3b517ab38": { + "views": [ + { + "cell_index": 7 + } + ] + } + }, + "version": "1.2.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Tutorials/12_MD/README.md b/Tutorials/12_MD/README.md new file mode 100644 index 00000000..fa3f602a --- /dev/null +++ b/Tutorials/12_MD/README.md @@ -0,0 +1,12 @@ +Ewald Summation +=============== + +This module seeks to provide an overview of the theory and implementation of Ewald summation, which is ubiquitous in molecular simuations of periodic systems. The following topics are covered: + +- Derivation of classical electrostatics from the Poisson equation. + +- A discussion of conditional convergence, which plagues conventional summation methods. + +- Development of the Ewald summation method, with examples, to overcome conditional convergence. + +- Conversion of the Ewald sum code to use particle-mesh Ewald (PME), which is more efficient. diff --git a/tests/test_TU_12.py b/tests/test_TU_12.py new file mode 100644 index 00000000..3d014ad3 --- /dev/null +++ b/tests/test_TU_12.py @@ -0,0 +1,15 @@ +from addons import * +from utils import * + + +tdir = 'Tutorials/12_MD' + +@using_matplotlib +def test_12a(workspace): + exe_scriptified_ipynb(workspace, tdir, '12a_basics') + +@using_matplotlib +@using_scipy +def test_12b(workspace): + exe_scriptified_ipynb(workspace, tdir, '12b_ewald') +