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

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[1], line 1
----> 1 get_ipython().run_line_magic('matplotlib', 'inline')
      2 ## Import libraries to plot and do math
      3 import matplotlib.pyplot as plt 

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/IPython/core/interactiveshell.py:2486, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2484     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2485 with self.builtin_trap:
-> 2486     result = fn(*args, **kwargs)
   2488 # The code below prevents the output from being displayed
   2489 # when using magics with decorator @output_can_be_silenced
   2490 # when the last Python token in the expression is a ';'.
   2491 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/IPython/core/magics/pylab.py:103, in PylabMagics.matplotlib(self, line)
     98     print(
     99         "Available matplotlib backends: %s"
    100         % _list_matplotlib_backends_and_gui_loops()
    101     )
    102 else:
--> 103     gui, backend = self.shell.enable_matplotlib(args.gui)
    104     self._show_matplotlib_backend(args.gui, backend)

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/IPython/core/interactiveshell.py:3758, in InteractiveShell.enable_matplotlib(self, gui)
   3755     import matplotlib_inline.backend_inline
   3757 from IPython.core import pylabtools as pt
-> 3758 gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
   3760 if gui != None:
   3761     # If we have our first gui selection, store it
   3762     if self.pylab_gui_select is None:

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/IPython/core/pylabtools.py:338, in find_gui_and_backend(gui, gui_select)
    321 def find_gui_and_backend(gui=None, gui_select=None):
    322     """Given a gui string return the gui and mpl backend.
    323 
    324     Parameters
   (...)    335     'WXAgg','Qt4Agg','module://matplotlib_inline.backend_inline','agg').
    336     """
--> 338     import matplotlib
    340     if _matplotlib_manages_backends():
    341         backend_registry = matplotlib.backends.registry.backend_registry

ModuleNotFoundError: No module named 'matplotlib'

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

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

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