%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

Morse potential

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

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

Force¶

F = -V.diff(r)
F

Plotting¶

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