# A simple L-J potential for simulating atomic interaction (with bug)

Hi, I am a student and recently started learning physics simulation with a strong interest.

Herein, I would like to show my simple code simulating the atomic interaction with Lennard-Jones (L-J) potential. And would like to ask some questions about the physical model.

Firstly, the L-J potential is a potential used to describe the interaction between two atoms and is widely used in molecular dynamics simulations. It has the following form and appearance:
L-J势函数是一种用于描述两原子间相互作用的势，被广泛用于分子动力学模拟之中。他有着如下的形式:

V(r) = 4\epsilon[(\frac{\sigma}{r})^{12} - (\frac{\sigma}{r})^{6}]

in which \sigma measures the location of the point where the potential energy is zero and \epsilon measures the depth of the potential well.

I obtained the force field by finding the gradient of this potential function and updated the velocity of each particle by the distance of different particles within a cutoff radius. Here I do not use any high-end operations, but simply a modification of the spring-mass system.

\pmb{F}(\pmb{r})=24[2\frac{\sigma^{12}}{r^{13}}- \frac{\sigma^{6}}{r^{7}}]\hat{\pmb{r}}

So, without other words, let me show this rough code:

import taichi as ti

ti.init()

max_num_particles = 1024
particle_mass = 1
dt = 1e-3
substeps = 10

sigma = 0.05
sigma2 = sigma**2
sigma6 = sigma2**3
sigma12 = sigma6**2
potential_energy = 0.1
cutoff_r = 0.5
k = 0.001

x = ti.Vector.field(2, dtype=ti.f32, shape=max_num_particles)
v = ti.Vector.field(2, dtype=ti.f32, shape=max_num_particles)
f = ti.Vector.field(2, dtype=ti.f32, shape=max_num_particles)
charge = ti.field(dtype=ti.f32, shape=max_num_particles)

num_particles = ti.field(dtype=ti.i32, shape=())

@ti.kernel
def add_ion(x_pos: ti.f32, y_pos: ti.f32, pos_neg: ti.f32):
new_particle_id = num_particles[None]

x[new_particle_id] = ti.Vector([x_pos, y_pos])
charge[new_particle_id] = pos_neg

num_particles[None] += 1

@ti.kernel
def substep():
n = num_particles[None]
for i in range(n):
f[i] = ti.Vector([0,0])
for j in range(n):
if i != j:
x_ij = x[i] - x[j]
if x_ij.norm() <= cutoff_r:
c_ij = charge[i] * charge[j]
delta_x = x_ij.norm()
x2 = delta_x**2
x6 = x2**3
x7 = x6 * delta_x
x13 = x6**2 * delta_x
d = x_ij.normalized()
f[i] += (24 * potential_energy * (2 * sigma12/x13 - sigma6/x7)) * d
# f[i] += -k * c_ij/ x_ij.norm()**2 * d
# Boundary conditions
for i in range(n):
v[i] += dt * f[i] / particle_mass
x[i] += dt * v[i]
for d in ti.static(range(2)):
if x[i][d] < 0:
x[i][d] = 0
v[i][d] = 0
if x[i][d] > 1:
x[i][d] = 1
v[i][d] = 0

def main():
gui = ti.GUI("Ionic Lattice Simulation",
res=(512,512),
background_color=0x262629)
while True:
for e in gui.get_events(ti.GUI.PRESS):
if e.key in (ti.GUI.ESCAPE, ti.GUI.EXIT):
exit()
elif e.key == ti.GUI.LMB:
elif e.key == ti.GUI.RMB:

for step in range(substeps):
substep()

# draw ions
X = x.to_numpy()
for i in range(num_particles[None]):
c = 0x3CC1FF if charge[i] < 0 else 0xFF623F

gui.show()

if __name__ == '__main__':
main()


There are actually two types of particles added to this code to model the attraction between positive and negative ions later on, but the Coulomb field is not implemented yet, just coloured differently.

The result is as follow

But I have a problem with this simulation. When there are a lot of particles, the simulation will suddenly get stuck and not move after a few dozen seconds, and I would like to ask what is happening.

Also, I would like to ask if there is any simple way to measure the parameters of the L-J potential function and the duration of the simulation that is realistic in a physical sense?

1 Like

I assume that you’re more comfortable with Chinese, if not, let me know.

1. 你的截断半径是sigma的10倍，一般如果使用带有cutoff的LJ势的时候截断半径会在sigma的2.5到4倍之间
2. 体系用了slip boundary condition，热力学上这会使体系的等效温度或能量不断降低，最终趋于接近固体的凝聚态，假如你的目标是分子模拟的话，多数情况下会用周期性边界条件。

ps. taichi真的太香了