python - Lennard-Jones potential simulation -

import pygame import random import numpy np import matplotlib.pyplot 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 = 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))  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:     event in pygame.event.get():         if event.type == pygame.quit:             running = false     screen.fill(background_colour)      i, particle in enumerate(my_particles):         particle.move()         particle.bounce()         particle2 in my_particles[i+1:]:             collide(particle, particle2)         particle.display()      pygame.display.flip() pygame.quit() 

i wanted simulate particles lennard-jones potential. problem code not know how use verlet algorithm.

1. i not know should implement verlet algorithm; inside class or outside?
2. how can use velocity verlet algorithm in move method?
3. is implementation of verlet algorithm correct, or should use arrays saving results?
4. what else should change make work?

1. you can keep dynamical variables, position , velocity, inside classes, each class needs acceleration vector. verlet integrator has role of controller, acts outside on collection of particles. keep angle out of computations, forth , trigonometric functions , inverses not necessary. make position, velocity , acceleration 2d vectors.

2. one way implement velocity verlet variant (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);     = eval_a(x);     v += a*0.5*dt;     do_statistics(t,x,v);  

   which supposes vectorized variant. in framework, there iterations on particle collection include,

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

3. no, not have done nothing wrong since did not call verlet function update position , velocity. , no, strict vectorization not necessary, see above. implicit vectorization via particles array sufficient.