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:

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


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=())

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

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",
    while True:
        for e in gui.get_events(ti.GUI.PRESS):
            if e.key in (ti.GUI.ESCAPE, ti.GUI.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):

        # 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)


if __name__ == '__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?


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体系分子模拟的模拟体系大小一般在几到几百纳米不等。


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


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

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

ps. taichi真的太香了