Week 5 Example 2 MD in a system of diatomic molecules in 2D

System Setup

%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=8

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

# Topology
M_gas_diatomic=np.array([[0, 1, 0, 0, 0, 0, 0, 0],
           [1, 0, 0, 0, 0, 0, 0, 0],
           [0, 0, 0, 1, 0, 0, 0, 0],
           [0, 0, 1, 0, 0, 0, 0, 0],
           [0, 0, 0, 0, 0, 1, 0, 0],
           [0, 0, 0, 0, 1, 0, 0, 0],
           [0, 0, 0, 0, 0, 0, 0, 1],
           [0, 0, 0, 0, 0, 0, 1, 0]]); 



######## "Force Field" Parameters #######
HS=1;   # Repulsive soft potential 
k=25.0; # Harmonic oscillator constant
req=1;  # Harmonic oscillator equilibrium distance
KAPPA=k*M_gas_diatomic
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      

        #Repulsive Wall + Harmonic potential
        dVdr=-12*HS/(np.power(r,13))+KAPPA*(r-req)

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

## 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 

## Initialise molecules in a reasonable initial configuration
for i in np.arange(0,N,2): 
    x0[i+1]=x0[i]+req/np.sqrt(2)
    y0[i+1]=y0[i]+req/np.sqrt(2)

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

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

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

Integration and Visualization

Initialise the system

## 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

## |------------------|
## |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[i,j],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;

|****** | 11%

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[4], line 28
     25     xnew[i]=verlet(x[i],xp[i],fx[i],m[i],dt) # new position (x-component)
     26     ynew[i]=verlet(y[i],yp[i],fy[i],m[i],dt); # new position (y-component)
---> 28 print_progress(timestep,nsteps)  
     30 # Reassign positions
     31 xp=x; yp=y; x=xnew+1-1; y=ynew+1-1;

Cell In[1], line 61, in print_progress(iteration, total, bar_length)
     58 arrow = '*' * int(round(bar_length * progress))
     59 spaces = ' ' * (bar_length - len(arrow))
---> 61 clear_output(wait=True)
     62 display(HTML(f"""
     63 <div style="color: blue;">
     64 |{arrow}{spaces}| {int(progress * 100)}%
     65 </div>
     66 """))

File ~/anaconda3/lib/python3.11/site-packages/IPython/core/display_functions.py:386, in clear_output(wait)
    384 from IPython.core.interactiveshell import InteractiveShell
    385 if InteractiveShell.initialized():
--> 386     InteractiveShell.instance().display_pub.clear_output(wait)
    387 else:
    388     print('\033[2K\r', end='')

File ~/anaconda3/lib/python3.11/site-packages/ipykernel/zmqshell.py:148, in ZMQDisplayPublisher.clear_output(self, wait)
    137 """Clear output associated with the current execution (cell).
    138 
    139 Parameters
   (...)
    145 
    146 """
    147 content = dict(wait=wait)
--> 148 self._flush_streams()
    149 assert self.session is not None
    150 msg = self.session.msg("clear_output", json_clean(content), parent=self.parent_header)

File ~/anaconda3/lib/python3.11/site-packages/ipykernel/zmqshell.py:67, in ZMQDisplayPublisher._flush_streams(self)
     65 """flush IO Streams prior to display"""
     66 sys.stdout.flush()
---> 67 sys.stderr.flush()

File ~/anaconda3/lib/python3.11/site-packages/ipykernel/iostream.py:578, in OutStream.flush(self)
    576     self.pub_thread.schedule(evt.set)
    577     # and give a timeout to avoid
--> 578     if not evt.wait(self.flush_timeout):
    579         # write directly to __stderr__ instead of warning because
    580         # if this is happening sys.stderr may be the problem.
    581         print("IOStream.flush timed out", file=sys.__stderr__)
    582 else:

File ~/anaconda3/lib/python3.11/threading.py:629, in Event.wait(self, timeout)
    627 signaled = self._flag
    628 if not signaled:
--> 629     signaled = self._cond.wait(timeout)
    630 return signaled

File ~/anaconda3/lib/python3.11/threading.py:331, in Condition.wait(self, timeout)
    329 else:
    330     if timeout > 0:
--> 331         gotit = waiter.acquire(True, timeout)
    332     else:
    333         gotit = waiter.acquire(False)

KeyboardInterrupt: 

Visualization of the trajectory

%%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=15, 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