--- jupytext: encoding: '# -*- coding: utf-8 -*-' formats: md:myst,ipynb text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- ```{code-cell} :tags: [hide-input] import mmf_setup; mmf_setup.nbinit() import os, sys from pathlib import Path FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/' os.makedirs(FIG_DIR, exist_ok=True) import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL) %matplotlib inline import numpy as np, matplotlib.pyplot as plt try: from myst_nb import glue except: glue = None from matplotlib.animation import FuncAnimation from phys_581.contexts import FPS ``` (sec:MolecularDynamics)= # Molecular Dynamics Another way to simulate material properties is with [molecular dynamics][] (MD). This simply uses Newton's laws to track the positions $\vect{r}_{i}$ and velocities $\vect{v}_{i} = \dot{\vect{r}}_{i}$ of many particles: \begin{gather*} m_i\ddot{\vect{r}}_i = \vect{F}_{i}(\{\vect{r},\vect{v}\}). \end{gather*} What makes MD simulations difficult is calculating the force $\vect{F}_{i}$, which may depend on many or all of the other particles. For example, if the particles interact through a potential $V(r)$, then we have formally \begin{gather*} \vect{F}_{i}= -\sum_{j\neq i}\vect{r}_{ij}\frac{V'(\norm{\vect{r}_{ij}})} {2\norm{\vect{r}_{ij}}}, \qquad \vect{r}_{ij} = \vect{r}_{i} - \vect{r}_{j}. \end{gather*} :::{admonition} Details: Getting the right sign. :class: dropdown This comes from \begin{gather*} \vect{F}_{i} = -\vect{\nabla}_{i}V(\{\vect{r}_{i}\}) = -\vect{\nabla}_{i}\sum_{i\neq j}\frac{1}{2}V(\norm{\vect{r}_{ij}})\\ = -\sum_{i\neq j}\frac{1}{2}V'(\norm{\vect{r}_{ij}})\vect{\nabla}_{i}\norm{\vect{r}_{ij}} = -\sum_{i\neq j}\frac{1}{2}V'(\norm{\vect{r}_{ij}}) \frac{\vect{r}_{ij}}{\norm{\vect{r}_{ij}}}. \end{gather*} If $V(r)$ is repulsive, $V'(r)$ will be negative, meaning that $\vect{F}_{i}$ will point in the direction of $\vect{r}_{ij}$, i.e., from $\vect{r}_{i}$ (tip) to $\vect{r}_{j}$ (tail) -- away from the particle as expected for a repulsive potential. ::: Naïvely, this requires summing over all the particles -- $O(N)$ operations -- for every particle, leading to a computational complexity of $O(N^2)$. Part of the art of MD is using clever approximations to reduce this complexity. See {cite}`Griebel:2007` for details. ## Lennard-Jones Potential A common example is to use a [Lennard-Jones potential][]: \begin{gather*} V(r) = \epsilon f\left(\frac{r}{\sigma}\right), \qquad f(q) = 4\frac{1-q^{6}}{q^{12}}, \qquad f'(q) = 24\frac{r^{6} - 2}{q^{13}} \end{gather*} ```{code-cell} :tags: [hide-input, margin] r = np.linspace(0, 3, 500)[1:] f = 4/r**6*(1/r**6 - 1) fig, ax = plt.subplots(figsize=(3,2), layout="constrained") ax.plot(r, f) ax.set(xlabel="$r/\sigma$", ylabel="$V(r)/\epsilon$", title="Lennard-Jones Potential", ylim=(-1.1, 1)); ``` This has a minimum at $r_\min = 2^{1/6}\sigma$ where $V'(r_\min) = 0$, which can be used to create a stable crystalline solid, thereby allowing us to simulate the elastic properties of materials in extreme conditions (i.e., where they explode). ## Time Evolution Care must be taken when evolving an MD simulation in time, otherwise, numerical errors will rapidly invalidate your simulation. Here we use the [Velocity-Størmer-Verlet][] integration scheme \begin{align*} \vect{r}(t+\delta t) &= \vect{r}(t) + \vect{v}(t) \delta t + \tfrac{1}{2}\vect{a}(t)\delta t^2,\\ \vect{a}(t+\delta t) &= \vect{F}\bigl(\{\vect{r}(t+\delta t)\}\bigr)/m,\\ \vect{v}(t+\delta t) &= \vect{v}(t) + \frac{\vect{a}(t) + \vect{a}(t+\delta t)}{2}\delta t. \end{align*} This requires storing the previous acceleration $\vect{a}(t)$ with each step. Here we present a simple simulation of a projectile impacting a slab (see Fig. 3.12 of `Griebel:2007`). This is a simple code but runs slowly as it suffers from the $O(N^2)$ complexity (hence the small resolution shown here). The philosophy here is to first code the simplest thing you can to make sure we understand the essential pieces of the algorithm, and that there are no roadblocks. Once you have this, it can form a test-bed for optimizing and exploring more complicated algorithms. ```{code-cell} :tags: [hide-input] _EPS = np.finfo(float).eps class MolecularDynamics: L1 = 250 L2 = 40 epsilon = 5 sigma = 1 m = 1 N1xy = (40, 40) N2xy = (160, 40) v1 = 0 - 10j v2 = 0 + 0j r_cut_sigma = 2.5 dt = 0.00005 beta = 0.0 sep = 0+5j # Initial separation (units of dx) # For colouring vmax = 5 cmap = plt.colormaps['viridis'] def __init__(self, **kw): for key in kw: if not hasattr(self, key): raise ValueError(f"Unknown {key=}") setattr(self, key, kw[key]) self.init() def init(self): self.dx = 2**(1/6) * self.sigma # Equilibrium spacing self.z0_1 = (-self.N1xy[0]/2 + (self.N1xy[1]-1)*1j + self.sep) * self.dx self.z0_2 = (-self.N2xy[0]/2 + 0j) * self.dx def F_r(self, r2): """Return F/r where F is the Lennard-Jones force for separation r. I.e., F_ij = (F/r)*r where r = r_j - r_i. Arguments ========= r2 : array-like Square of the separation distance. """ # Prevent divide-by-zero errors r2_ = r2 + _EPS # s = (sigma/r)**6 s = np.where(r2 > 0, (self.sigma**2 / r2_)**2, 0) s *= s**2 f = 24 * self.epsilon * s / r2_ * (1 - 2*s) return f def get_initial_state(self): (X1, Y1), (X2, Y2) = [ np.meshgrid(*[np.arange(N)* self.dx + x0 for (N, x0) in zip(Nxy, (z0.real, z0.imag))], indexing='ij') for Nxy, z0 in [(self.N1xy, self.z0_1), (self.N2xy, self.z0_2)] ] Z1 = X1 + 1j*Y1 Z2 = X2 + 1j*Y2 V1 = 0*Z1 + self.v1 V2 = 0*Z2 + self.v2 Z, V = [np.concatenate(list(map(np.ravel, A))) for A in ((Z1, Z2), (V1, V2))] return (0, Z, V, 0*Z) def plot(self, state, ax=None): if ax is None: fig, ax = plt.subplots(layout="constrained") else: fig = ax.figure t, z, v, a = state ax.scatter(z.real, z.imag, c=self.cmap(abs(v)/self.vmax)) ax.set(title=f"{t=:.4f}", aspect=1) return ax def get_F(self, z): """Return the total force on particles at z.""" r = z[None, :] - z[:, None] return (self.F_r(r2=abs(r)**2) * r).sum(axis=1) def step(self, state): """Return the acceleration on each particle. This uses the Velocity-Størmer-Verlet method. """ dt = self.dt t0, z0, v0, a0 = state z1 = z0 + dt * v0 + dt**2/2 * a0 a1 = self.get_F(z1) / self.m v1 = v0 + dt / 2 * (a0 + a1) return (t0 + dt, z1, v1, a1) md = MolecularDynamics() md = MolecularDynamics(N1xy=(4*3, 4*3), N2xy=(16*3, 4*3), sep=0.4+4j) ``` ```{code-cell} :tags: [hide-input] %%time md = MolecularDynamics( v1=-5j, beta=0.0, N1xy=(4*2, 4*2), N2xy=(16*2, 4*2), sep=-10.3+2j, dt=0.00005*40) s = md.get_initial_state() fig, ax = plt.subplots(layout="constrained") for frame in FPS(100, fig=fig, embed=True): for n in range(20): s = md.step(s) ax.cla() ax = md.plot(s, ax=ax) ax.set(xlim=(-12*2, 12*2), ylim=(-5, 10*2)) ``` [Molecular dynamics]: [Velocity-Størmer-Verlet]: [Lennard-Jones potential]: [Van der Waals interaction]: