Python-Particle motion simulation
Thing you want to do
Simulate the movement of atoms using Lennard-Jones potential as follows
ideal
http://www.engineering-eye.com/AUTODYN/function/index.html
Potential and force between atoms
Lennard Jones Potential
Lennard-Jones potential when the distance between two particles is $ r_ {ij} $:
\begin{eqnarray}
U(r_{ij})&=&4ε((\frac{σ}{r_{ij}})^{12}-(\frac{σ}{r_{ij}})^{6})\\
&=&4ε(\frac{σ}{r_{ij}})^{6}((\frac{σ}{r_{ij}})^{6}-1)
\end{eqnarray}
For $ N $ particles
\begin{eqnarray}
U(\vec{x_1},...,\vec{x_N})&=&4ε\sum_{i=1}^{N}\sum_{j=i+1}^{N}((\frac{σ}{r_{ij}})^{12}-(\frac{σ}{r_{ij}})^{6})\\
\end{eqnarray}
The force acting on the atom $ i $ under the Lennard-Jones potential is due to the gradient at the interatomic distance $ r $.
\begin{eqnarray}
\vec{F}_i&=&-∇U(\vec{x_1},...,\vec{x_N})\\
&=&24ε\sum_{j=1,j≠i}^{N}\frac{1}{r_{ij}^2}(\frac{σ}{r_{ij}})^{6}(1-2(\frac{σ}{r_{ij}})^{6})\vec{r}_{ij}\\
\end{eqnarray}
https://ja.wikipedia.org/wiki/%E3%83%AC%E3%83%8A%E3%83%BC%E3%83%89-%E3%82%B8%E3%83%A7%E3%83%BC%E3%83%B3%E3%82%BA%E3%83%BB%E3%83%9D%E3%83%86%E3%83%B3%E3%82%B7%E3%83%A3%E3%83%AB
Cutoff distance and cell split
To reduce the amount of calculation, a cutoff radius of $ r_ {cut} $ is provided, and forces from particles (on the grid) farther than $ r_ {cut} $ are ignored as their contribution is small.
Divide the space into cubes of $ r_ {cut} × r_ {cut} × r_ {cut} $ and manage which particles belong to each grid.
For particles in a particular grid, only the force from particles in adjacent grids is calculated.
Calculation of particle position and velocity
The equation of motion of the second-order differential equation is transformed into the first-order simultaneous differential equation and calculated by the Runge-Kutta method.
\begin{eqnarray}
m\frac{d^2\vec{r}}{dt^2}=\vec{F}\\
\end{eqnarray}
↓
\begin{eqnarray}
\frac{d\vec{r}}{dt}&=&\vec{f}(\vec{r},\vec{v},t)=\vec{v}\\
\frac{d\vec{v}}{dt}&=&\vec{g}(\vec{r},\vec{v},t)=\frac{\vec{F}(\vec{r},\vec{v},t)}{m}\\
\end{eqnarray}
Update by Runge-Kutta method
\begin{eqnarray}
\vec{k_1} &=& dt×\vec{f}(\vec{r},\vec{v},t)=dt×\vec{v}\\
\vec{l_1} &=& dt×\vec{g}(\vec{r},\vec{v},t)=dt×\frac{\vec{F}(\vec{r},\vec{v},t)}{m}\\
\vec{k_2} &=& dt×\vec{f}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})=dt×(\vec{v}+\frac{\vec{l_1}}{2})\\
\vec{l_2} &=& dt×\vec{g}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})=dt×\frac{\vec{F}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})}{m}\\
\vec{k_3} &=& dt×\vec{f}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})=dt×(\vec{v}+\frac{\vec{l_2}}{2})\\
\vec{l_3} &=& dt×\vec{g}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})=dt×\frac{\vec{F}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})}{m}\\
\vec{k_4} &=& dt×\vec{f}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)=dt×(\vec{v}+\vec{l_3})\\
\vec{l_4} &=& dt×\vec{g}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)=dt×\frac{\vec{F}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)}{m}\\
\vec{r}(t+dt)&=&\vec{r}(t)+\frac{\vec{k_1}+2\vec{k_2}+2\vec{k_3}+\vec{k_4}}{6}\\
\vec{v}(t+dt)&=&\vec{v}(t)+\frac{\vec{l_1}+2\vec{l_2}+2\vec{l_3}+\vec{l_4}}{6}\\
\end{eqnarray}
Position $ \ vec {r} (t + dt) $ at $ t + dt $ ・ Velocity $ \ vec {v} (t + dt) $ is calculated sequentially by the above formula.
Implementation
Base class
import numpy as np
from abc import abstractmethod
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
import random
import math
import datetime
class Particle():
def __init__(self, pos = [0.0, 0.0, 0.0], vel = [0.0, 0.0, 0.0], force = [0.0, 0.0, 0.0], mass = 1.0, type = -1):
self.pos = np.array(pos)
self.vel = np.array(vel)
self.force = np.array(force)
self.mass = mass
self.type = type
#The force acting between particles under Lennard-Jones potential
@staticmethod
def force(pos1, pos2, eps = 1.0, sigma = 1.0):
#return np.array([0, 0, 0])
r2 = ((pos2 - pos1)**2).sum()
if abs(r2) < 0.00001:
return np.array([0.0, 0.0, 0.0])
sig6_div_r6 = (sigma**2 / r2)**3
#print("force : ", 24 * eps * (1/r2) * sig6_div_r6 * (1 - 2 * sig6_div_r6) * (pos2 - pos1))
return 24 * eps * (1/r2) * sig6_div_r6 * (1 - 2 * sig6_div_r6) * (pos2 - pos1)
class Calculater():
#The position of the particle after dt(pos)And speed(vel)Is updated by the Runge-Kutta method
@staticmethod
def runge_kutta(particle, adjuscent_particles, ps, dt, t, eps, sigma):
k1 = dt * particle.vel
l1 = 0
for p in adjuscent_particles:
l1 += dt * Particle.force(pos1 = particle.pos, pos2 = p.pos, eps = eps, sigma = sigma) / particle.mass
l1 += dt * ps.force(particle.pos, particle.vel, particle, t)
k2 = dt * (particle.vel + k1 / 2)
l2 = 0
for p in adjuscent_particles:
l2 = dt * Particle.force(particle.pos + l1/2, p.pos, eps, sigma) / particle.mass
l2 += dt * ps.force(particle.pos + l1/2, particle.vel + k1/2, particle, t + dt/2)
k3 = dt * (particle.vel + k2 / 2)
l3 = 0
for p in adjuscent_particles:
l3 = dt * Particle.force(particle.pos + l2/2, p.pos, eps, sigma) / particle.mass
l3 += dt * ps.force(particle.pos + l2/2, particle.vel + k2/2, particle, t + dt/2)
k4 = dt * (particle.vel + k3)
l4 = 0
for p in adjuscent_particles:
l4 = dt * Particle.force(particle.pos + l3, p.pos, eps, sigma) / particle.mass
l4 += dt * ps.force(particle.pos + l3, particle.vel + k3, particle, t + dt)
particle.pos += (k1 + 2*k2 + 2*k3 + k4)/6
particle.vel += (l1 + 2*l2 + 2*l3 + l4)/6
#The position of the particle after dt(pos)And speed(vel)Is updated by the Euler method
@staticmethod
def euler(particle, adjuscent_particles, ps, dt, t, eps, sigma):
f = 0
for p in adjuscent_particles:
f += dt * Particle.force(pos1 = particle.pos, pos2 = p.pos, eps = eps, sigma = sigma) / particle.mass
f += dt * ps.force(particle.pos, particle.vel, particle, t)
particle.pos += dt * particle.vel
particle.vel += dt * f
class ParticleSystem():
#
def __init__(self, cutoff_r, NX, NY, NZ, e = 0.9, mass = 1, eps = 1, sigma = 1):
self.cutoff_r = cutoff_r
self.NX = NX
self.NY = NY
self.NZ = NZ
self.X_MAX = cutoff_r * NX
self.Y_MAX = cutoff_r * NY
self.Z_MAX = cutoff_r * NZ
self.e = e
self.mass = mass
self.eps = eps
self.sigma = sigma
self.time = 0.0
self.prev_timestamp = 0.0
self.fps = 0.1
self.particles = []
self.snapshots = []
self.fig = plt.figure()
self.ax = self.fig.add_subplot(111, projection='3d')
self.gif_output_mode = False
self.init_particles()
self.update_cells()
self.setup_particle_type_dict()
#Find the force acting on all particles and the position / velocity after dt.
def update(self, dt, calc_method):
self.time += dt
#update_forces_on_particle()
self.update_pos_and_vel(dt, calc_method)
self.update_cells()
'''
#Calculate the force acting on all particles.
#Abolished because it is unlikely to be usable by the Runge-Kutta method.
def update_forces_on_particle():
for particle in self.particles:
particle.force = 0
#Calculate the intermolecular force acting on particles
for p in get_particles_adjuscent_cell(particle):
f = Particle.force(particle.pos, p.pos)
particle.force += f
particle.force += self.force(particle)
'''
#Calc the position and time after dt of all particles_Calculate with method.
def update_pos_and_vel(self, dt, calc_method):
for particle in self.particles:
calc_method(particle = particle, adjuscent_particles = self.get_particles_adjuscent_cell(particle), ps = self, dt = dt, t = self.time, eps = self.eps, sigma = self.sigma)
#Make it bounce at both ends of the space.(Bounce coefficient e)
if particle.pos[0] < 0:
#particle.pos[0] += 2*(0 - particle.pos[0])
particle.pos[0] = 0.0
particle.vel[0] = -self.e * particle.vel[0]
if particle.pos[0] > self.X_MAX:
#particle.pos[0] -= 2*(particle.pos[0] - self.X_MAX)
particle.pos[0] = self.X_MAX - 0.0001
particle.vel[0] = -self.e * particle.vel[0]
if particle.pos[1] < 0:
#particle.pos[1] += 2*(0 - particle.pos[1])
particle.pos[1] = 0.0
particle.vel[1] = -self.e * particle.vel[1]
if particle.pos[1] > self.Y_MAX:
#particle.pos[1] -= 2*(particle.pos[1] - self.Y_MAX)
particle.pos[1] = self.Y_MAX - 0.0001
particle.vel[1] = -self.e * particle.vel[1]
if particle.pos[2] < 0:
#particle.pos[2] += 2*(0 - particle.pos[2])
particle.pos[2] = 0.0
particle.vel[2] = -self.e * particle.vel[2]
if particle.pos[2] > self.Z_MAX:
particle.pos[2] = self.Z_MAX - 0.0001
#particle.pos[2] -= 2*(particle.pos[2] - self.Z_MAX)
particle.vel[2] = -self.e * particle.vel[2]
if self.gif_output_mode is True:
self.take_snapshot(save_to_snapshots = True)
#Get a list of particles contained in cells containing particle cells and 26 adjacent cells
def get_particles_adjuscent_cell(self, particle):
particles = []
index = self.get_cellindex_from_pos(particle.pos)
for i in range(index[0] - 1, index[0] + 1):
for j in range(index[1] - 1, index[1] + 1):
for k in range(index[2] - 1, index[2] + 1):
try:
for p in self.cell[i][j][k] if self.cell[i][j][k] is not None else []:
if p is not particle:
particles.append(p)
except IndexError:
continue
return particles
#Calculate the cell index from the position of the particles.
def get_cellindex_from_pos(self, pos):
return [int(pos[0] / self.cutoff_r), int(pos[1] / self.cutoff_r), int(pos[2] / self.cutoff_r)]
#Update which cell each particle belongs to
def update_cells(self):
self.cell = [[[ [] for i in range(self.NX)] for j in range(self.NY)] for k in range(self.NZ)]
for p in self.particles:
try:
index = self.get_cellindex_from_pos(p.pos)
self.cell[index[0]][index[1]][index[2]].append(p)
except Exception as e:
print("Exception : ", e)
print("pos : ", p.pos)
#input()
#Get particle list
def get_particles(self):
return self.particles
#Get time
def get_time(self):
return self.time
#Take a snapshot of the state of the particles. It is not displayed.
def take_snapshot(self, save_to_snapshots = False):
if self.time - self.prev_timestamp > self.fps:
self.ax.set_title("Time : {}".format(self.time))
self.ax.set_xlim(0, self.X_MAX)
self.ax.set_ylim(0, self.Y_MAX)
self.ax.set_zlim(0, self.Z_MAX)
scats = []
for type, particles in self.particle_dict.items():
x_list = []
y_list = []
z_list = []
for p in particles[0]:
x_list.append(p.pos[0])
y_list.append(p.pos[1])
z_list.append(p.pos[2])
scats.append(self.ax.scatter(x_list, y_list, z_list, c=particles[1]))
if save_to_snapshots is True:
self.snapshots.append(scats)
print(len(self.snapshots))
self.prev_timestamp = self.time
#Display the state of particles
def show_snapshot(self):
if self.time - self.prev_timestamp > 0.1:
self.fig = plt.figure()
self.ax = self.fig.add_subplot(111, projection='3d')
self.take_snapshot()
#plt.savefig('box_gravity_time-{}.png'.format(self.time))
plt.show()
self.prev_timestamp = self.time
#Write mode on to gif file
def start_output_gif(self, fps = 0.1):
self.gif_output_mode = True
self.snapshots = []
self.take_snapshot(save_to_snapshots = True)
#gif file write mode off. Write to file
def stop_output_gif(self, filename = "hoge.gif"):
print("stop_output_gif : ", len(self.snapshots))
self.gif_output_mode = False
ani = animation.ArtistAnimation(self.fig, self.snapshots)
ani.save(filename, writer='imagemagick')
self.snapshots = []
#Make a dictionary for each type of particle.
def setup_particle_type_dict(self):
self.particle_dict = {}
for p in self.particles:
if str(p.type) not in self.particle_dict:
#Assign a random color for each particle type
self.particle_dict[str(p.type)] = [[], tuple([random.random() for _ in range(3)])]
self.particle_dict[str(p.type)][0].append(p)
#Calculate the kinetic energy of all particles
def get_kinetic_energy(self):
return sum([(p.mass*np.array(p.vel)**2).sum()/2 for p in self.particles])
#Particle initialization
@abstractmethod
def init_particles(self):
raise NotImplementedError()
#Power to work as a system(Gravity etc.)
@abstractmethod
def force(self, pos, vel, particle, t):
raise NotImplementedError()
Box-shaped particle army system under gravity
#Box-shaped particle army system under gravity
class BoxGravitySystem(ParticleSystem):
def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
#Generates a box-shaped 10x10x10 particle army at the bottom. The initial speed is set to go diagonally upward.
def init_particles(self):
for x in np.linspace(0.1*self.X_MAX, 0.2*self.X_MAX, 10):
for y in np.linspace(0.1*self.Y_MAX, 0.2*self.Y_MAX, 10):
for z in np.linspace(0.1*self.Z_MAX, 0.2*self.Z_MAX, 10):
self.particles.append(Particle(pos = [x, y, z], vel=[0.1*self.X_MAX, 0.05*self.Y_MAX, 0.5*self.Z_MAX], mass = self.mass))
#gravity
def force(self, pos, vel, particle, t):
return np.array([0.0, 0.0, -particle.mass * 9.8])
if __name__ == '__main__':
cutoff_r = 1.0
e = 0.1
mass = 1.0
eps = 1
sigma = 0.5
dt = 0.001
system = BoxGravitySystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
time = 0
system.start_output_gif(fps = 0.1)
while time <= 5:
system.update(dt = dt, calc_method = Calculater.runge_kutta)
time = system.get_time()
print("Time : ", time)
#system.show_snapshot()
system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))
Large mass (no interatomic interaction)
A system that shoots bullets into the wall
#A system that shoots bullets into the wall
class BulletWallSystem(ParticleSystem):
def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
#Sphere with initial velocity(bullet)And set up a fixed wall
def init_particles(self):
#ball(bullet)Generate
bullet_center = [0.35*self.X_MAX, 0.5*self.Y_MAX, 0.5*self.Z_MAX]
for i in range(200):
r = 0.05 * self.X_MAX * (random.random() - 0.5) * 2
phi = 2 * np.pi * random.random()
theta = np.pi * random.random()
self.particles.append(Particle(pos = [bullet_center[0] + r*np.sin(theta)*np.cos(phi), bullet_center[1] + r*np.sin(theta)*np.sin(phi), bullet_center[2] + r*np.cos(theta)], vel=[0.2*self.X_MAX, 0, 0], mass = self.mass, type = 1))
#Wall generation
for x in np.linspace(0.49*self.X_MAX, 0.50*self.X_MAX, 2):
for y in np.linspace(0.3*self.Y_MAX, 0.7*self.Y_MAX, 30):
for z in np.linspace(0.3*self.Z_MAX, 0.7*self.Z_MAX, 30):
self.particles.append(Particle(pos = [x, y, z], vel=[0.0, 0.0, 0.0], mass = self.mass, type = 2))
#No external force
def force(self, pos, vel, particle, t):
return np.array([0.0, 0.0, 0])
if __name__ == '__main__':
cutoff_r = 1.0
e = 0.1
mass = 1.0
eps = 1
sigma = 0.5
dt = 0.001
system = BulletWallSystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
time = 0
system.start_output_gif(fps = 0.1)
while time <= 5:
system.update(dt = dt, calc_method = Calculater.runge_kutta)
time = system.get_time()
print("Time : ", time)
#system.show_snapshot()
system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))
cutoffr：1.0、ε：0.001、σ：0.1、dt：0.002
cutoffr：1.0、ε：0.02、sigma：0.5
cutoffr：1.0、ε：0.02、sigma：1.3、dt-0.001
cutoffr：1.0、ε：0.3、σ：1.3、dt：0.01
cutoffr：1.0、ε：0.1、σ：1.2、dt：0.01
cutoffr：1.0、ε：0.1、σ：1.5、dt：0.001
cutoffr：1.0、ε:1、σ：1.5、dt：0.001
A system that receives a force that rotates like a tornado
#A system that receives a force that rotates like a tornado
class TornadeSystem(ParticleSystem):
def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
#Create 3000 particles at random positions
def init_particles(self):
for _ in range(3000):
x = self.X_MAX * random.random()
y = self.Y_MAX * random.random()
z = self.Z_MAX * random.random()
self.particles.append(Particle(pos = [x, y, z], vel=[0.0, 0.0, 0.0], mass = self.mass, type = 1))
#A force that rotates around the z-axis(rot(r)∝[0,0,1]Power to become) × z
def force(self, pos, vel, particle, t):
return 0.05 * pos[2] * np.array([-(pos[1] - self.Y_MAX/2), pos[0] - self.X_MAX/2 , 0.0])# - 0.5 * np.array(vel)
if __name__ == '__main__':
cutoff_r = 1.0
e = 0.1
mass = 1.0
eps = 1
sigma = 0.5
dt = 0.001
system = TornadeSystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
time = 0
system.start_output_gif(fps = 0.1)
while time <= 5:
system.update(dt = dt, calc_method = Calculater.runge_kutta)
time = system.get_time()
print("Time : ", time)
#system.show_snapshot()
system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))
Failure work
Summary
Adjusting the parameters doesn't work very well
ideal
reality
