Molecular dynamic simulation using velocity verlet in Python

Question

I'm trying to implement a simple MD simulation in Python (I'm new to this),I'm using LJ potential and force equations along with Verlet method. I posted another question about the same code yesterday and finally made progress thanks to you! Here is my code:

def LJ_VF(r):
    #r = distance in Å
    #Returns V in (eV) and F in (eV/Å)
    V = 4 * epsilon * ( (sigma/r)**(12) - (sigma/r)**6 )
    F = 24 * epsilon * ( 2 * ((sigma**12)/(r**(13))) - ( (sigma**6)/(r**7) )) 
    return V , F

def velocity_verlet(x, v, f_old, f_new):                   #setting m=1 so that a=f
    x_new = x + v * dt + 0.5 * f_old * dt**2  
    v_new = v + 0.5 * (f_old + f_new) * dt  
    return x_new, v_new

Now to make sure it works I'm trying to use this on the simplest system, i.e 2 atoms separated by the distance 4 Å:

#Constants

epsilon = 0.0103     
sigma = 3.4  
m = 1.0
t0 = 0.0
v0 = 0.0
dt = 0.1 
N = 1000

def simulate_two_atoms(p1_x0, p1_v0, p2_x0, p2_v0):
    p1_x, p2_x = [p1_x0], [p2_x0]
    p1_v, p2_v = [p1_v0], [p2_v0]
    p1_F, p1_V, p2_F, p2_V = [], [], [], []

    r = abs(p2_x0 - p1_x0)
    V, F = LJ_VF(r)
    p1_F.append(F)
    p1_V.append(V)
    p2_F.append(-F)
    p2_V.append(V)

    for i in range(N - 1):
        r_new = abs(p2_x[-1] - p1_x[-1])  
        V_new, F_new = LJ_VF(r_new)  
        x1_new, v1_new = velocity_verlet(p1_x[-1], p1_v[-1], p1_F[-1], F_new)
        x2_new, v2_new = velocity_verlet(p2_x[-1], p2_v[-1], p2_F[-1], -F_new)

        p1_x.append(x1_new)
        p1_v.append(v1_new)
        p2_x.append(x2_new)
        p2_v.append(v2_new)
        
        p1_F.append(F_new)
        p2_F.append(-F_new)
        
        p1_V.append(V_new)
        p2_V.append(V_new)
    
    return np.array(p1_x), np.array(p1_v), np.array(p2_x), np.array(p2_v)

#Initial conditions

p1_x0 = 0.0
p1_v0 = 0.0  
p2_x0 = 4.0  
p2_v0 = 0.0 

#Plot 

p1_x, p1_v, p2_x, p2_v = simulate_two_atoms(p1_x0, p1_v0, p2_x0, p2_v0)
time = np.arange(N)  

plt.plot(time, p1_v, label="Particle 1", color="red")
plt.plot(time, p2_v, label="Particle 2", color="orange")
plt.xlabel("Time (t)")
plt.ylabel("v_x(t)")
plt.title("Velocities v_x(t)")
plt.legend()
plt.show()

Now when I'm plotting the velocities, I can see that they're increasing and hence the energies are not conserved while they should be!

I have been trying to find where I went wrong for a very long time now but not progressing in any level! Thought that a second eye might help! Any tips/advice?

Answer 1

(Your post reverted to the older code with the forces the wrong way round, so it doesn't produce your plot. Your plot comes from one with the forces corrected.)

You are slowly losing accuracy because you use F_new to calculate the velocity, but F_new was last calculated for positions x_old. Basically, F_new does not reflect latest position.

You should (in order):

• calculate new positions (using current position, current velocity, old force);
• calculate a new force (which depends on the new position);
• calculate a new velocity (which depends on the current velocity and the average of old and new forces).

The plot below is created by the code that follows. (As a matter of principle please make sure that you include a single, complete, runnable code, including imports.)

import numpy as np
import matplotlib.pyplot as plt

def LJ_VF(r):
    #r = distance in A
    #Returns V in (eV) and F in (eV/A)
    V = 4 * epsilon * ( (sigma/r)**(12) - (sigma/r)**6 )
    F = 24 * epsilon * ( 2 * ((sigma**12)/(r**(13))) - ( (sigma**6)/(r**7) )) 
    return V , F

def position_verlet(x, v, f_old):
    return x + v * dt + 0.5 * f_old * dt**2

def velocity_verlet(v, f_old, f_new):
    return v + 0.5 * (f_old + f_new) * dt


#Constants

epsilon = 0.0103     
sigma = 3.4  
m = 1.0
t0 = 0.0
v0 = 0.0
dt = 0.1 
N = 1000

def simulate_two_atoms(p1_x0, p1_v0, p2_x0, p2_v0):
    p1_x, p2_x = [p1_x0], [p2_x0]
    p1_v, p2_v = [p1_v0], [p2_v0]
    p1_F, p1_V, p2_F, p2_V = [], [], [], []

    r = p2_x0 - p1_x0
    V, F = LJ_VF(r)
    p1_F.append(-F)
    p1_V.append(V)
    p2_F.append(F)
    p2_V.append(V)

    for i in range(N - 1):
        x1_new = position_verlet(p1_x[-1], p1_v[-1], p1_F[-1] )      # update POSITION
        x2_new = position_verlet(p2_x[-1], p2_v[-1], p2_F[-1] )

        r_new = x2_new - x1_new
        V_new, F_new = LJ_VF(r_new)                                  # update FORCE

        v1_new = velocity_verlet(p1_v[-1], p1_F[-1], -F_new)         # update VELOCITY
        v2_new = velocity_verlet(p2_v[-1], p2_F[-1],  F_new)

        p1_x.append(x1_new)
        p2_x.append(x2_new)
        p1_v.append(v1_new)
        p2_v.append(v2_new)

        p1_F.append(-F_new)
        p2_F.append( F_new)
        p1_V.append( V_new)
        p2_V.append( V_new)
    
    return np.array(p1_x), np.array(p1_v), np.array(p2_x), np.array(p2_v)

#Initial conditions

p1_x0 = 0.0
p1_v0 = 0.0  
p2_x0 = 4.0  
p2_v0 = 0.0 

#Plot 

p1_x, p1_v, p2_x, p2_v = simulate_two_atoms(p1_x0, p1_v0, p2_x0, p2_v0)
time = np.arange(N)  

plt.plot(time, p1_v, label="Particle 1", color="red")
plt.plot(time, p2_v, label="Particle 2", color="orange")
plt.xlabel("Time (t)")
plt.ylabel("v_x(t)")
plt.title("Velocities v_x(t)")
plt.legend()
plt.show()