python - Lennard-Jones potential simulation

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

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. 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?

Answer 1

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.

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

3. 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.

Code with some add-ons but without collisions, desingularized potential

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