MD Example 1 in a system of LJ particles 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=15 #x dimension of the simulation' box
boxy=15 #y dimension of the simulation' box
boxk=0.1 #k constant for harmonic repulsive force
#Number of particles
N=20
#mass of the particles
m=np.ones(N)
# Topology
M=np.array([[0, 1, 0, 0, 1, 0, 0, 0, 0, 0],
[1, 0, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 1, 1, 0, 1, 0],
[1, 0, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 0, 1, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 0, 0, 1],
[0, 0, 0, 1, 0, 0, 0, 0, 0, 1]]);
M=np.array([[0, 1, 0, 0, 0, 0, 0, 0, 0, 1],
[1, 0, 1, 0, 0, 0, 0, 0, 0, 1],
[0, 1, 0, 1, 0, 0, 0, 0, 0, 1],
[0, 0, 1, 0, 1, 0, 0, 0, 0, 1],
[1, 0, 0, 1, 0, 1, 0, 0, 0, 1],
[0, 0, 0, 0, 1, 0, 1, 0, 0, 1],
[0, 0, 0, 0, 0, 1, 0, 1, 0, 1],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]);
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=np.zeros([N,N]) #k*M
epsilon=0.1;
distance=lambda xi,xj,yi,yj : np.sqrt((xi-xj)**2+(yi-yj)**2)
##Use this function to implement different potentials
def forceij(xi,xj,yi,yj,r,HS,KAPPA,req,epsilon):
#Lennard Jones Potential
dVdr=-48*epsilon*(1/(np.power(r,13))-0.5/(np.power(r,7)))
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.randint(-2*boxx, 2*boxx, N)/2; #Initial position in x
y0=np.random.randint(-2*boxy, 2*boxy, N)/2; #Initial position in y
# Random initial velocities
v0=0.01*(np.random.rand(2,N)-0.5); # Initial random velocitites
## Define the timestep and the total time
dt=0.05; # Timestep
nsteps=50000; # Total number of steps
total_time=dt*nsteps; # Total simulation time
## 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):
rij = distance(x[i],x[j],y[i],y[j])
if rij <= 2 :
[cx,cy]=forceij(x[i],x[j],y[i],y[j],rij,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;
|*********** | 21%
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[4], line 29
26 xnew[i]=verlet(x[i],xp[i],fx[i],m[i],dt) # new position (x-component)
27 ynew[i]=verlet(y[i],yp[i],fy[i],m[i],dt); # new position (y-component)
---> 29 print_progress(timestep,nsteps)
31 # Reassign positions
32 xp=x; yp=y; x=xnew+1-1; y=ynew+1-1;
Cell In[1], line 82, in print_progress(iteration, total, bar_length)
79 arrow = '*' * int(round(bar_length * progress))
80 spaces = ' ' * (bar_length - len(arrow))
---> 82 clear_output(wait=True)
83 display(HTML(f"""
84 <div style="color: blue;">
85 |{arrow}{spaces}| {int(progress * 100)}%
86 </div>
87 """))
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:66, in ZMQDisplayPublisher._flush_streams(self)
64 def _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:573, in OutStream.flush(self)
562 """trigger actual zmq send
563
564 send will happen in the background thread
565 """
566 if (
567 self.pub_thread
568 and self.pub_thread.thread is not None
(...)
571 ):
572 # request flush on the background thread
--> 573 self.pub_thread.schedule(self._flush)
574 # wait for flush to actually get through, if we can.
575 evt = threading.Event()
File ~/anaconda3/lib/python3.11/site-packages/ipykernel/iostream.py:266, in IOPubThread.schedule(self, f)
264 self._events.append(f)
265 # wake event thread (message content is ignored)
--> 266 self._event_pipe.send(b"")
267 else:
268 f()
File ~/anaconda3/lib/python3.11/site-packages/zmq/sugar/socket.py:696, in Socket.send(self, data, flags, copy, track, routing_id, group)
689 data = zmq.Frame(
690 data,
691 track=track,
692 copy=copy or None,
693 copy_threshold=self.copy_threshold,
694 )
695 data.group = group
--> 696 return super().send(data, flags=flags, copy=copy, track=track)
File zmq/backend/cython/socket.pyx:742, in zmq.backend.cython.socket.Socket.send()
File zmq/backend/cython/socket.pyx:789, in zmq.backend.cython.socket.Socket.send()
File zmq/backend/cython/socket.pyx:250, in zmq.backend.cython.socket._send_copy()
File ~/anaconda3/lib/python3.11/site-packages/zmq/backend/cython/checkrc.pxd:13, in zmq.backend.cython.checkrc._check_rc()
KeyboardInterrupt:
Visualization of the trajectory#
%%capture
## Display the trajectory
%matplotlib inline
import matplotlib
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);
#matplotlib.rcParams['animation.embed_limit'] = 2**128
HTML(anim.to_jshtml())