{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Molecular dynamics simulations\n", "\n", "Classical molecular dynamics simulations involve numerically integrating Newton's equations of motion for $N$ particles interacting through a pair potential $V(\\mathbf{r}_1,\\mathbf{r}_2)$:\n", "$$\n", "m_i \\ddot{\\mathbf{r}}_i = - \\sum_{j \\neq i} \\nabla_i V(\\mathbf{r}_i,\\mathbf{r}_j)~.\n", "$$\n", "\n", "If the system is sufficiently large, one can use MD simulations to study its macroscopic (thermodynamic) properties such as the dependence of the pressure on temperature and density -- the equation of state.\n", "\n", "In many case the pair potential depends only on the distance $r_{ij} = |\\mathbf{r}_i - \\mathbf{r}_j|$.\n", "In this case the equations of motion read\n", "$$\n", "m_i \\ddot{\\mathbf{r}}_i = - \\sum_{j \\neq i} \\frac{d V(r_{ij})}{d r_{ij}} \\, \\frac{\\mathbf{r}_i - \\mathbf{r}_j}{r_{ij}}~.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lennard-Jones fluid\n", "\n", "Lennard-Jones fluid corresponds to a pair potential\n", "$$\n", "V_{LJ}(r) = 4 \\varepsilon \\left[ (\\sigma/r)^{12} - (\\sigma/r)^6 \\right].\n", "$$\n", "\n", "The Lennard-Jones system has a rich phase diagram containing various phases (gas, liquid, solid) and the associated phase transitions.\n", "\n", "\n", "\n", "To evaluate the properties of the LJ fluid, one would run molecular dynamics simulations for sufficiently long time to let the system equilibrate, and calculate its macroscopic properties like pressure as time average. The simulations are typically performed in a box with periodic boundary conditions to minimize finite-size effects." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Velocity Verlet method\n", "\n", "For long molecular dynamics simulations it makes sense to use stable, time-reversal numerical schemes that conserve energy. The leapfrog method is a good choice. When applied to Newton's equations of motion it is know as Verlet method.\n", "\n", "Suppose we are solving Newton's equations of motion\n", "$$\n", "\\frac{d^2 x}{dt^2} = f(x,t).\n", "$$\n", "These are rewritten as\n", "\\begin{align*}\n", "\\frac{dx}{dt} &= v, \\\\\n", "\\frac{dv}{dt} &= f(x,t).\n", "\\end{align*}\n", "\n", "If we apply the leapfrog scheme to this system of equations, it will look like\n", "\\begin{align*}\n", "x(t+h) & = x(t) + h v(t + h/2), \\\\\n", "v(t+3h/2) &= v(t+h/2) + h f[x(t+h), t+h)],\n", "\\end{align*}\n", "i.e. the coordinates are evaluated are full steps using velocity estimates at half-steps,\n", "and vice versa.\n", "\n", "This is the essence of the Verlet method.\n", "It's feature is that we do not need to keep track of coordinates at half-steps (only the velocities) and of velocities at full steps. In case we are interested in velocity values also at full step, these can be estimated as $v(t+h) = v(t+h/2) + (h/2) f[x(t+h),t+h]$.\n", "\n", "The full algorithm is as follows. Given initial values $x(t)$ and $v(t)$ one first computes\n", "$$\n", "v(t+h/2) = v(t) + \\frac{h}{2} f[x(t),t].\n", "$$\n", "Then each step of the Verlet algorithm corresponds to\n", "\\begin{align*}\n", "v(t+h/2) & = v(t) + \\frac{h}{2} f[x(t),t],\\\\\n", "x(t+h) & = x(t) + h v(t+h/2), \\\\\n", "v(t+h) & = v(t+h/2) + \\frac{h}{2} f[x(t+h),t+h)].\n", "\\end{align*}\n", "\n", "The algorithm is straightforwardly generalizable to more than one variable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us first write a code for integrating Newton's equations of motion for arbitrary central potential" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "# Simulation parameters (can be changed later)\n", "n_particles = 125 # Number of particles\n", "density = 0.8 # Density of the fluid\n", "box_length = (n_particles/density)**(1/3.) # Length of simulation box (determined by the number of particles and density)\n", "temperature0 = 1.0 # Initial (desired) temperature of the system\n", "time_step = 0.004 # Time step for the simulation\n", "num_steps = 1000 # Number of simulation steps\n", "\n", "\n", "# Initialization\n", "# Initial position at the nodes of the cube\n", "# Helps to avoid particle overlap\n", "def initial_positions():\n", " ret = np.zeros((n_particles,3))\n", " Nsingle = np.ceil(n_particles**(1/3.))\n", " dL = box_length / Nsingle\n", " for i in range(n_particles):\n", " ix = i % Nsingle\n", " iy = (np.trunc(i / Nsingle)) % Nsingle\n", " iz = np.trunc(i / (Nsingle * Nsingle))\n", " ret[i][0] = (ix + 0.5) * dL;\n", " ret[i][1] = (iy + 0.5) * dL;\n", " ret[i][2] = (iz + 0.5) * dL;\n", " return ret\n", "\n", "# Computes forces for a given vector of positions, interaction potential and its gradient\n", "# Return a tuple: positions, total potential energy, and the virial part of the pressure\n", "def compute_forces(forces, positions, potential, potential_gradient):\n", " # forces = np.zeros_like(positions)\n", " forces.fill(0.)\n", " potential_energy = 0.0\n", " virial = 0.0\n", " for i in range(n_particles):\n", " for j in range(i+1, n_particles):\n", " # Vector of relative distance\n", " r_ij = positions[i] - positions[j]\n", " # Periodic boundary conditions (minimum-image convention)\n", " r_ij = r_ij - box_length*np.round(r_ij/box_length)\n", " \n", " r_sq = np.sum(r_ij**2)\n", " f_ij = -potential_gradient(r_sq) * r_ij\n", " \n", " forces[i] += f_ij\n", " forces[j] -= f_ij\n", " potential_energy += potential(r_sq)\n", " virial += np.dot(f_ij, r_ij)\n", " \n", " virial = virial/(3.0*box_length**3)\n", " return forces, potential_energy, virial\n", "\n", "\n", "# Compute the kinetic temperature from velocities (kinetic energy) \n", "def compute_kinetic_temperature(velocities):\n", " return np.sum(velocities**2) / 3. / n_particles\n", "\n", "# Renormalize the velocities to match the desired kinetic temperature \n", "def renormalize_velocities(velocities, temperature = temperature0):\n", " Tkin = compute_kinetic_temperature(velocities)\n", " factor = np.sqrt(temperature / Tkin)\n", " velocities *= factor\n", "\n", "\n", "# Apply the velocity verlet time step\n", "# Returns the tuple of new positions, velocities, accelerations (forces), potential energy and pressure\n", "def velocity_verlet(positions, velocities, accelerations, \n", " time_step, potential, potential_gradient):\n", " # Update positions\n", " positions += velocities*time_step + 0.5*accelerations*time_step**2\n", " positions = positions - box_length*np.floor(positions/box_length)\n", " # Update velocities\n", " velocities_half = velocities + 0.5*accelerations*time_step\n", " # Compute new forces and potential energy\n", " accelerations, potential_energy, pressure = compute_forces(accelerations, positions, potential, potential_gradient)\n", " # Update velocities using new accelerations\n", " velocities = velocities_half + 0.5*accelerations*time_step\n", " # Add ideal gas contribution to the pressure\n", " kinetic_temperature = compute_kinetic_temperature(velocities)\n", " pressure += density * kinetic_temperature\n", " return positions, velocities, accelerations, potential_energy, pressure\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the Lennard-Jones potential properties" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Lennard-Jones potential as a function of squared distance\n", "def lj_potential(r_sq):\n", " r6 = r_sq**3\n", " r12 = r6**2\n", " return 4.0*(1./r12 - 1./r6)\n", "\n", "# The grandient term dV/dr / r in the rhs of Newton's equations\n", "# for the LJ potential\n", "def lj_potential_gradient(r_sq):\n", " r6 = r_sq**3\n", " r12 = r6**2\n", " return -24.0*(2./r12 - 1./r6) / r_sq" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Perform the simulation" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "