python - Lennard-Jones 势模拟

Question

import pygame
import random
import numpy as np
import matplotlib.pyplot as plt
import math

number_of_particles = 70
my_particles = []
background_colour = (255,255,255)
width, height = 500, 500
sigma = 1
e = 1
dt = 0.1
v = 0
a = 0
r = 1

def r(p1,p2):
    dx = p1.x - p2.x
    dy = p1.y - p2.y
    angle = 0.5 * math.pi - math.atan2(dy, dx)
    dist = np.hypot(dx, dy)
    return dist

def collide(p1, p2):
    dx = p1.x - p2.x
    dy = p1.y - p2.y

    dist = np.hypot(dx, dy)
    if dist < (p1.size + p2.size):
        tangent = math.atan2(dy, dx)
        angle = 0.5 * np.pi + tangent

        angle1 = 2*tangent - p1.angle
        angle2 = 2*tangent - p2.angle
        speed1 = p2.speed
        speed2 = p1.speed
        (p1.angle, p1.speed) = (angle1, speed1)
        (p2.angle, p2.speed) = (angle2, speed2)

        overlap = 0.5*(p1.size + p2.size - dist+1)
        p1.x += np.sin(angle) * overlap
        p1.y -= np.cos(angle) * overlap
        p2.x -= np.sin(angle) * overlap
        p2.y += np.cos(angle) * overlap
def LJ(r):
    return -24*e*((2/r*(sigma/r)**12)-1/r*(sigma/r)**6)

def verlet():
    a1 = -LJ(r(p1,p2))
    r = r + dt*v+0.5*dt**2*a1
    a2 = -LJ(r(p1,p2))
    v = v + 0.5*dt*(a1+a2)
    return r, v
class Particle():
    def __init__(self, (x, y), size):
        self.x = x
        self.y = y
        self.size = size
        self.colour = (0, 0, 255)
        self.thickness = 1
        self.speed = 0
        self.angle = 0

    def display(self):
        pygame.draw.circle(screen, self.colour, (int(self.x), int(self.y)), self.size, self.thickness)

    def move(self):
        self.x += np.sin(self.angle) 
        self.y -= np.cos(self.angle) 

    def bounce(self):
        if self.x > width - self.size:
            self.x = 2*(width - self.size) - self.x
            self.angle = - self.angle

        elif self.x < self.size:
            self.x = 2*self.size - self.x
            self.angle = - self.angle

        if self.y > height - self.size:
            self.y = 2*(height - self.size) - self.y
            self.angle = np.pi - self.angle

        elif self.y < self.size:
            self.y = 2*self.size - self.y
            self.angle = np.pi - self.angle             

screen = pygame.display.set_mode((width, height))

for n in range(number_of_particles):
    x = random.randint(15, width-15)
    y = random.randint(15, height-15)

    particle = Particle((x, y), 15)
    particle.speed = random.random()
    particle.angle = random.uniform(0, np.pi*2)

    my_particles.append(particle)

running = True
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
    screen.fill(background_colour)

    for i, particle in enumerate(my_particles):
        particle.move()
        particle.bounce()
        for particle2 in my_particles[i+1:]:
            collide(particle, particle2)
        particle.display()

    pygame.display.flip()
pygame.quit()

我想通过伦纳德-琼斯势来模拟粒子。我对这段代码的问题是我不知道如何使用 Verlet 算法。

1. 我不知道应该在哪里实现 Verlet 算法；课内还是课外？
2. 如何在 move 方法中使用 Verlet 算法中的速度？
3. 我对 Verlet 算法的实现是否正确，还是应该使用数组来保存结果？
4. 我还应该更改哪些内容才能使其正常工作？

Answer 1

1. 您可以将动态变量、位置和速度保留在类实例内，但是每个类都需要一个加速度向量来累积力贡献。 Verlet 积分器具有 Controller 的作用，它从外部作用于所有粒子的集合。将角度排除在计算之外，不需要用三角函数及其反函数来回计算。使位置、速度和加速度都是二维向量。

2. 实现速度 Verlet 变体的一种方法是(参见 https://stackoverflow.com/tags/verlet-integration/info)

   verlet_step:
       v += a*0.5*dt;
       x += v*dt; t += dt;
       do_collisions(t,x,v,dt);
       a = eval_a(x);
       v += a*0.5*dt;
       do_statistics(t,x,v); 
   

   它假设一个矢量化变体。在您的框架中，将包含一些对粒子集合的迭代，

   verlet_step:
       for p in particles:
           p.v += p.a*0.5*dt; p.x += p.v*dt; 
       t += dt;
       for i, p1 in enumerate(particles):
           for p2 in particles[i+1:]:
               collide(p1,p2);
       for i, p1 in enumerate(particles):
           for p2 in particles[i+1:]:
               apply_LJ_forces(p1,p2);
       for p in particles:
           p.v += p.a*0.5*dt;
       do_statistics(t,x,v); 
   

3. 不，您不可能没有做错任何事，因为您实际上没有调用 Verlet 函数来更新位置和速度。不，严格的矢量化是没有必要的，见上文。通过粒子数组进行隐式矢量化就足够了。如果您想与使用相同模型提供 ODE 函数的标准积分器(例如 scipy.integrate 中的积分器)的结果进行比较，则只需要完全矢量化。

带有一些附加组件但没有冲突的代码，具有去奇异化的潜力

import pygame
import random
import numpy as np
import matplotlib.pyplot as plt
import math

background_colour = (255,255,255)
width, height = 500, 500
aafac = 2 # anti-aliasing factor screen to off-screen image

number_of_particles = 50
my_particles = []

sigma = 10
sigma2 = sigma*sigma
e = 5
dt = 0.1 # simulation time interval between frames
timesteps = 10 # intermediate invisible steps of length dt/timesteps  

def LJ_force(p1,p2):
    rx = p1.x - p2.x
    ry = p1.y - p2.y

    r2 = rx*rx+ry*ry

    r2s = r2/sigma2+1
    r6s = r2s*r2s*r2s
    f = 24*e*( 2/(r6s*r6s) - 1/(r6s) )

    p1.ax += f*(rx/r2)
    p1.ay += f*(ry/r2)
    p2.ax -= f*(rx/r2)
    p2.ay -= f*(ry/r2)

def Verlet_step(particles, h):
    for p in particles:
        p.verlet1_update_vx(h); 
        p.bounce()
    #t += h;
    for i, p1 in enumerate(particles):
        for p2 in particles[i+1:]:
            LJ_force(p1,p2);
    for p in particles:
        p.verlet2_update_v(h);

class Particle():
    def __init__(self, (x, y), (vx,vy), size):
        self.x = x
        self.y = y
        self.vx = vx
        self.vy = vy
        self.size = size
        self.colour = (0, 0, 255)
        self.thickness = 2
        self.ax = 0
        self.ay = 0

    def verlet1_update_vx(self,h):
        self.vx += self.ax*h/2
        self.vy += self.ay*h/2
        self.x += self.vx*h
        self.y += self.vy*h
        self.ax = 0
        self.ay = 0

    def verlet2_update_v(self,h):
        self.vx += self.ax*h/2
        self.vy += self.ay*h/2

    def display(self,screen, aa):
        pygame.draw.circle(screen, self.colour, (int(aa*self.x+0.5), int(aa*self.y+0.5)), aa*self.size, aa*self.thickness)

    def bounce(self):
        if self.x > width - self.size:
            self.x = 2*(width - self.size) - self.x
            self.vx = - self.vx

        elif self.x < self.size:
            self.x = 2*self.size - self.x
            self.vx = - self.vx

        if self.y > height - self.size:
            self.y = 2*(height - self.size) - self.y
            self.vy = - self.vy

        elif self.y < self.size:
            self.y = 2*self.size - self.y
            self.vy = - self.vy            

#------------ end class particle ------------
#------------ start main program ------------

for n in range(number_of_particles):
    x = 1.0*random.randint(15, width-15)
    y = 1.0*random.randint(15, height-15)
    vx, vy = 0., 0.
    for k in range(6):
        vx += random.randint(-10, 10)/2.
        vy += random.randint(-10, 10)/2.

    particle = Particle((x, y),(vx,vy), 10)

    my_particles.append(particle)

#--------- pygame event loop ----------
screen = pygame.display.set_mode((width, height))
offscreen = pygame.Surface((aafac*width, aafac*height))

running = True
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
    offscreen.fill(background_colour)

    for k in range(timesteps):
        Verlet_step(my_particles, dt/timesteps)

    for particle in my_particles:
        particle.display(offscreen, aafac)

    pygame.transform.smoothscale(offscreen, (width,height), screen)
    pygame.display.flip()
pygame.quit()