Ctrl+K

Molecular Dynamics : simulation of a crystal formation in 2D

%matplotlib widget
from IPython.display import Image
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import animation, rc
import sympy as sp
import pandas as pd
from PMD import PMD, distances

sp.init_printing(use_latex="mathjax")
rc("animation", html="html5")

Molecular Dynamics : simulation of a crystal formation in 2D#

Atomic simulation using the Morse potential

https://en.wikipedia.org/wiki/Morse_potential

This notebook uses the (Point Mass Dynamics) PMD class to simulate gravitational interaction between massive objects.

Required files

Before using this notebook, download the module PMD.py and put it in your working directory.

Image(
    url="https://upload.wikimedia.org/wikipedia/commons/7/7a/Morse-potential.png",
    width=500,
    height=500,
)

PMD module#

PMD module available here: PMD.py

Potential#

De, a, re, r = sp.symbols("D_e a r_e r")
V = De * (1 - sp.exp(-a * (r - re))) ** 2
V
\[\displaystyle D_{e} \left(1 - e^{- a \left(r - r_{e}\right)}\right)^{2}\]

Force#

F = -V.diff(r)
F
\[\displaystyle - 2 D_{e} a \left(1 - e^{- a \left(r - r_{e}\right)}\right) e^{- a \left(r - r_{e}\right)}\]

Plotting#

Hide code cell source 
values = {De: 1.0, a: 2.0, re: 1}
Vf = sp.lambdify(r, V.subs(values), "numpy")
Ff = sp.lambdify(r, F.subs(values), "numpy")
vr = np.linspace(0.3 * values[re], 5 * values[re], 100)
fig = plt.figure()
ax = fig.add_subplot(2, 1, 1)
plt.ylim(0, 5)
plt.plot(vr, Vf(vr))
plt.grid()
plt.ylabel("Potential Value, $V$")
ax = fig.add_subplot(2, 1, 2)
plt.ylim(-2.5, 5)
plt.plot(vr, Ff(vr))
plt.grid()
plt.xlabel("Interatomic Distance, $r$")
plt.ylabel("Force Value, $F$")
plt.show()

class MD(PMD):
    """
    Molecular dynamics w. Morse potential.
    """

    def __init__(self, De=1.0, a=1.0, re=1.0, mu=0.0, **kwargs):
        self.De = De
        self.a = a
        self.re = re
        self.mu = mu
        super().__init__(**kwargs)

    def derivative(self, X, t, cutoff_radius=1.0e-2):
        De, a, re, mu = self.De, self.a, self.re, self.mu
        n = len(m)
        P = X[: 2 * n].reshape(n, 2)
        V = X[2 * n :].reshape(n, 2)
        M = m * m[:, np.newaxis]
        D, R, U = distances(P)
        np.fill_diagonal(R, np.inf)
        if cutoff_radius > 0.0:
            R = np.where(R > cutoff_radius, R, cutoff_radius)
        # F =((G * M * R**-2)[:,:,np.newaxis] * U).sum(axis = 0)
        F = (
            (2.0 * De * a * (1.0 - np.exp(-a * (R - re))) * np.exp(-a * (R - re)))[
                :, :, np.newaxis
            ]
            * U
        ).sum(axis=0) - mu * (V**2).sum(axis=0) ** 1.5 * V
        A = (F.T / m).T
        X2 = X.copy()
        X2[: 2 * n] = V.flatten()
        X2[2 * n :] = A.flatten()
        return X2

    def potential(self, P):
        De, a, re, mu = self.De, self.a, self.re, self.mu
        D, R, U = distances(P)
        return (
            np.where(R != 0.0, De * a * (1 - np.exp(-a * (R - re))) ** 2, 0.0).sum()
            / 2.0
        )

np.random.seed(666)
re = 0.9
nmx = 6
nmy = 6
nm = nmx * nmy
x0 = np.arange(nmx) * re * 1.2
y0 = np.arange(nmy) * re * 1.2
X0, Y0 = np.meshgrid(x0, y0)
P0 = np.array([X0.flatten(), Y0.flatten()]).T
P0 *= 1.0 + np.random.rand(*P0.shape) * 0.1
P0 -= P0.mean(axis=0)
V0 = np.zeros_like(P0)
pcolors = "r"
tcolors = "b"
m = np.ones(nm) * 1.0e0
s = MD(m=m, P=P0, V=V0, mu=0.05, re=re, a=10, De=0.5, nk=1000)
dt = 1.0e-1
nt = 100


fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect("equal")
ax.axis("off")
margin = 1.0
plt.xlim(P0[:, 0].min() - margin, P0[:, 0].max() + margin)
plt.ylim(P0[:, 1].min() - margin, P0[:, 1].max() + margin)
# plt.grid()
# ax.axis("off")
points = []

msize = 10.0 * (s.m / s.m.max()) ** (1.0 / 6.0)
for i in range(nm):
    plc = len(pcolors)
    pc = pcolors[i % plc]
    tlc = len(tcolors)
    tc = tcolors[i % tlc]
    (trail,) = ax.plot([], [], "-" + tc)
    (point,) = ax.plot([], [], "o" + pc, markersize=msize[i])
    points.append(point)
    points.append(trail)


def init():
    for i in range(2 * nm):
        points[i].set_data([], [])
    return points


def animate(i):
    s.solve(dt, nt)  # , rtol = 1.e-8, atol = 1.e-8)
    x, y = s.xy()
    for i in range(nm):
        points[2 * i].set_data(x[i : i + 1], y[i : i + 1])
        xt, yt = s.trail(i)
        points[2 * i + 1].set_data(xt, yt)
    return points


anim = animation.FuncAnimation(
    fig, animate, init_func=init, frames=400, interval=20, blit=True
)


# plt.show()
plt.close()
anim
Contents