Interactive Molecular Dynamics

Daniel V. Schroeder

Jaro Camphuijsen & Rahiel Kasim

Overview

  • Introduction
  • Goal
  • Molecular Dynamics
    • Ingredients
    • A Molecular Dynamics Program
  • Interactive MD

Introduction

“If, in some cataclysm, all of scientific knowledge were to be destroyed, and only one sentence passed on to the next generation of creatures, what statement would contain the most information in the fewest words? I believe it is the atomic hypothesis (or the atomic fact, or whatever you wish to call it) that all things are made of atoms—little particles that move around in perpetual motion, attracting each other when they are a little distance apart, but repelling upon being squeezed into one another.”

—Richard Feynman

Goal

Disadvantages of Conventional Education

  • No physics, only mathematics
  • Only ideal gasses, equilibria, thermodynamic limit

Why Interactive MD

  • Increase qualitative physical understanding
  • Perform experiments on inaccessible systems

Molecular Dynamics

  1. Position your atoms
  2. Evolve according to Newton
  3. Wait until it reaches equilibrium
  4. Start measuring $\left(E_\text{kin}, E_\text{pot}\right)$
  5. Compute quantities $\left(P, T, \dots\right)$
  6. Take the time average

Ingredients

Potential: Lennard-Jones

$$ u(r) = 4 \epsilon \left[ \left(\frac{\sigma}{r} \right)^{12} - \left(\frac{\sigma}{r}\right)^6\right]$$
$ u(r) = 4 \epsilon \left[ \left(\frac{\sigma}{r} \right)^{12} - \left(\frac{\sigma}{r}\right)^6\right]$

Force: Newton

  • $$ F = m a$$
  • $$ F = -\frac{du}{dr} $$

Integration: Velocity Verlet

\begin{align} r \left(t+\Delta t \right)&=r(t) + v(t) \Delta t + \frac{f(t)}{2m}\Delta t^2\\ \quad\\ v \left(t+\Delta t \right)&=v(t) + \frac{f(t+\Delta t) + f(t)}{2m}\Delta t \end{align}

Building a Molecular Dynamics program

A Simple Molecular Dynamics Program


def main():
    init()
    t = 0
    while (t < tmax):        # MD loop
        force(f, en)         # determine forces
        integrate(f, en)     # integrate equations of motion
        t = t + delt
        sample()             # sample averages

Force Calculation


def force(f, en):
    for i in range(npart - 1):
        for j in range(i + 1, npart):  # loop over all pairs
            xr = x[i] - x[j]
            xr = xr - box * round(xr / box)
            r2 = xr ** 2
            if (r2 < rc2):      # test cutoff
                r2i = 1 / r2
                r6i = r2i ** 3
                ff = 48 * r2i * r6i * (r6i - 0.5)  # L-J potential
                f[i] = f[i] + ff * xr
                f[j] = f[j] - ff * xr
                en = en + 4 * r6i * (r6i - 1) - ecut
    return

Questions...?

Follow the course Understanding Molecular Simulation in Period 3