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.
在这里我想展示下我编写的用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势函数的参数以及模拟时长是符合现实物理意义的?

1 Like

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

你的最外层for loop的循环长度是一个ti.field定义的变量,一般来说并行的struct for需要循环次数是编译时常量比如max_num_particles,否则编译时可能会出问题。

对无量纲的LJ体系模拟,我们一般使用reduced units,见https://en.m.wikipedia.org/wiki/Lennard-Jones_potential#Dimensionless%20(reduced%20units)
一般来说 ,sigma和epsilon的值大概在原子单位附近,比如球形甲烷分子用LJ势能描述的sigma大约在0.373 nm, epsilon大约在148 K (再乘以Boltzmann常数)。对应地,LJ体系分子模拟的模拟体系大小一般在几到几百纳米不等。

我还注意到另外2点:

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

编程愉快~

谢谢!之后会尝试修改一下代码看看,另外我发现似乎taichi自带的mass_spring_game example在粒子多一些时运行一段时间也会卡住。还有就是我之后直接在程序运行时放置了粒子,程序能正常长时间运行,而如果是用鼠标放置粒子进去则会在一定时间内卡死。

目前还在想如何合理添加周期性边界条件。之后应该还会继续完善,希望能在之后有将真实的分子体系建模出来的能力 :grinning:

ps. taichi真的太香了