Lennard-Jones势模拟 [英] Lennard-Jones potential simulation

问题描述

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

我想用Lennard-Jones势来模拟粒子.我的这段代码的问题是我不知道如何使用Verlet算法.

I wanted to simulate particles by Lennard-Jones potential. My problem with this code is that I do not know how to use the Verlet algorithm.

1. 我不知道应该在哪里实现Verlet算法；在班级内部还是外部?
2. 如何在move方法中使用Verlet算法中的速度?
3. 我的Verlet算法实现是否正确，还是应该使用数组保存结果?
4. 要使其正常工作，我还需要更改什么?

1. I do not know where I should implement the Verlet algorithm; inside the class or outside?
2. How can I use velocity from the Verlet algorithm in the move method?
3. Is my implementation of the Verlet algorithm correct, or should I use arrays for saving results?
4. What else should I change to make it work?

推荐答案

1. 您可以在类实例中保留动态变量，位置和速度，但是每个类都需要一个加速度矢量来累积力的贡献. Verlet积分器具有控制器的作用，它从外部作用于所有粒子的收集.不要将角度放在计算之外，不要使用三角函数来回计算角度，也不必使用它们的逆函数.将位置，速度和加速度设为所有2D向量.

1. You can keep the dynamical variables, position and velocity, inside the class instances, however then each class needs an acceleration vector to accumulate the force contributions. The Verlet integrator has the role of a controller, it acts from outside on the collection of all particles. Keep the angle out of the computations, the forth and back with trigonometric functions and their inverses is not necessary. Make position, velocity and acceleration all 2D vectors.

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

One way to implement the velocity Verlet variant is (see 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); 

，它假定为向量化变体.在您的框架中，粒子集合将包含一些迭代，

which supposes a vectorized variant. In your framework, there would be some iterations over the particle collection to include,

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

不，您没有做错任何事情，因为您实际上没有调用Verlet函数来更新位置和速度.不，没有必要进行严格的矢量化处理，请参见上文.通过particles数组进行隐式矢量化就足够了.如果要与scipy.integrate等标准积分器的结果进行比较，并使用同一模型提供ODE功能，则只需要完整的向量化即可.

No, you could not have done nothing wrong since you did not actually call the Verlet function to update position and velocity. And no, a strict vectorization is not necessary, see above. The implicit vectorization via the particles array is sufficient. You would only need a full vectorization if you wanted to compare with the results of a standard integrator like those in scipy.integrate using the same model to provide the ODE function.

具有一些附加功能但没有冲突的代码，具有潜在的奇异功能

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

这篇关于Lennard-Jones势模拟的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！