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.
在这里我想展示下我编写的用L-J势函数模拟原子间相互作用简陋代码，并想在此问些问题。
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.
其中 \sigma 衡量势能为0点的位置， \epsilon 衡量势阱深度。
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:
add_ion(e.pos[0], e.pos[1], 1)
elif e.key == ti.GUI.RMB:
add_ion(e.pos[0], e.pos[1], -1)
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.circle(X[i], color=c, radius=5)
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?
同时我还想问一下，是否有什么简单方法来衡量L-J势函数的参数以及模拟时长是符合现实物理意义的？
