# 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

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

### Why Interactive MD

• Increase qualitative physical understanding
• Perform experiments on inaccessible systems

# Molecular Dynamics

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