Week 5 exercise: Simple MD simulations

In this session we will implememnt simple MD simulations of increasing complexity. Using this notebook as a starting point (also accessible on Google colab):

• Implement a dynamic simulation of an ideal gas in 2D.

• Introduce an interaction \(\Gamma_{ii}(r)=r^{-12}\) potential to model a system of soft spheres, in a box defined by harminic repulsive walls.

• Introduce a Lennard Jones interaction potential.

• Observe what happens by performing simulations at different values of \(\epsilon\)

• Introduce an Harmonic potential to model diatomic molecules.

• If you are inclined to adventure, introduce a three body angle potential to model non-linear, triatomic pseudo molecules.

Model of a Fluid Phase in a 2D box: A template

Work with the function forceij to introduce different interaction potentials:

Ideal Gas:

\(\Gamma_{ij}=0\)

Soft sphere:

\(\Gamma_{ij}=r_{ij}^{-12}\)

Soft sphere + harmonic:

\(\Gamma_{ij}=k_{ij}(r_{ij}-r_{eq})^2+r_{ij}^{-12}\)

Lennard Jones:

\(\Gamma_{ij}=4\epsilon\left[r_{ij}^{-12}-r_{ij}^{-6}\right]\)

Define Functions

%matplotlib inline
## Import libraries to plot and do math
import matplotlib.pyplot as plt 
import numpy as np
from IPython.display import clear_output, display, HTML

## Useful functions
verlet=lambda r, r_past, force, mass, dt:  2*r-r_past+(dt**2)*force/mass
forcebox=lambda x, boxx,boxk: np.greater(np.abs(x),boxx)*(-boxk)*x


#Define the system's box
boxx=10 #x dimension of the simulation' box
boxy=10 #y dimension of the simulation' box
boxk=1  #k constant for harmonic repulsive force

#Number of particles
N=50

#mass of the particles
m=np.ones(N)

######## "Force Field" Parameters #######
HS=1; # Repulsive soft potential 
k=25.0; # Harmonic oscillator constant
req=1; # Harmonic oscillator equilibrium distance
KAPPA=k*np.zeros([N,N])
epsilon=0;


##Use this function to implement different potentials
def forceij(xi,xj,yi,yj,HS,KAPPA,req,epsilon): 
        r=np.sqrt((xi-xj)**2+(yi-yj)**2); #Distance      

        #Ideal Gas
        dVdr=0 
        #Repulsive Wall 
        # dVdr=-12*HS/(np.power(r,13))
        #Repulsive Wall + Harmonic potential
        #... ... ...
        #Lennard Jones Potential
        #... .... .... ... ....         
        cx=-dVdr*((xi-xj))/r;  #Pairwise component of the force in x
        cy=-dVdr*((yi-yj))/r;  #Pairwise component of the force in y
             
        return [cx,cy]


def print_progress(iteration, total, bar_length=50):
    progress = (iteration / total)
    arrow = '*' * int(round(bar_length * progress))
    spaces = ' ' * (bar_length - len(arrow))
    
    clear_output(wait=True)
    display(HTML(f"""
    <div style="color: blue;">
    |{arrow}{spaces}| {int(progress * 100)}%
    </div>
    """))

System Setup

## Set the initial Conditions
# Random initial positions
x0=(np.random.rand(N)*2*boxx)-(boxx);     #Initial position in x
y0=(np.random.rand(N)*2*boxy)-(boxy);     #Initial position in y 

# Random initial velocities
v0=(np.random.rand(2,N)-0.5); # Initial random velocitites

## Define the timestep and the total time
dt=0.01; # Timestep
total_time=25;  # Total simulation time
nsteps=int(total_time/dt); # Total number of steps

## Initialise vectors 
time=np.zeros(nsteps)

## Compute a trajectory with the Verlet Algorithm
# Initialise positions at t-dt
xp=x0;
yp=y0;

# Position at time t
x=xp+v0[0,:]*dt;
y=yp+v0[1,:]*dt;

# Position at time t+dt
xnew=np.zeros(N);
ynew=np.zeros(N);

# time
time=np.arange(0,nsteps);
time[0]=0;
time[1]=time[0]+dt;

## Initialize verctors for plotting 
xx=np.zeros((np.size(time),N));xx[0]=x0
yy=np.zeros((np.size(time),N));yy[0]=y0

## |------------------|
## |Compute trajectory|
## |------------------|
for timestep in np.arange(1,nsteps): #Cycle over timesteps
    timestep=int(timestep)           #Make sure timestep is an integer
    
    # Initialise force vectors
    fx=np.zeros(N);  
    fy=np.zeros(N); 
    
    # Cycle over all particles
    for i in np.arange(0,N):
        fx[i]+=forcebox(x[i],boxx,boxk)
        fy[i]+=forcebox(y[i],boxy,boxk)
        for j in np.arange(i+1,N):
            
            [cx,cy]=forceij(x[i],x[j],y[i],y[j],HS,KAPPA,req,epsilon)
            
            fx[i]=fx[i]+cx;      #update total x-component of the force on particle i
            fx[j]=fx[j]-cx;      #update total x-component of the force on particle j 

            fy[i]=fy[i]+cy;      #update total y-component of the force on particle i
            fy[j]=fy[j]-cy;      #update total y-component of the force on particle j 

        xnew[i]=verlet(x[i],xp[i],fx[i],m[i],dt) # new position (x-component)
        ynew[i]=verlet(y[i],yp[i],fy[i],m[i],dt); # new position (y-component)

    print_progress(timestep,nsteps)  

    # Reassign positions
    xp=x; yp=y; x=xnew+1-1; y=ynew+1-1;

    ## Store trajectory for animation 
    xx[timestep]=x;
    yy[timestep]=y;

|**************************************************| 99%

Visualization

%%capture
## Display the trajectory
%matplotlib inline
from matplotlib.animation import FuncAnimation
from matplotlib import animation, rc
from IPython.display import HTML

fig, ax = plt.subplots(figsize=(5, 5))
line, = ax.plot([]) 
ax.set_xlim(-boxx, boxx)
ax.set_ylim(-boxy, boxy)
line, = ax.plot([], [], lw=2, marker='o', markersize=10, markerfacecolor=(0.8, 1.0, 0.8, 0.5),
             markeredgewidth=1,  markeredgecolor=(0, 0, 0, .5), linestyle=' ',color='red')
# initialization function: plot the background of each frame
def init():
    line.set_data([], [])
    return (line,)

def animate(frame_num):
    x=xx[frame_num,:]
    y=yy[frame_num,:]
    line.set_data((x, y))
    return (line,)

# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=np.arange(1,int(nsteps),100), interval=50);

HTML(anim.to_jshtml())

Once Loop Reflect