<jupyter_start><jupyter_text>Stochastic Processes: Data Analysis and Computer Simulation  Brownian motion 3: data analyses and applications -Interacting Brownian particles-  Note 1- So far we have successfully performed simulations for independent Brownian particles and confirmed that the simulation results are consistent with the corresponding theoretical predictions.- In this lesson, we modify our simulations to calculate the diffusive dynamics of interacting Brownian particles for which fully analytical predictions are not available.- Theoretical calculations based on such models, despite their apparent simplicity, can be carried out only by the use of computer simulations.   Necessary changes for interacting Brownian particles Periodic boundary conditions![](fig/PBC.png)  Note 2- Because we already have a simulation code for independent Brownian particles, modifying it for interacting Brownian particles is rather simple.- There are two main changes that we must implement, namely the inter-particle forces and the use of periodic boundary conditions.- The former is obvious, if we want to consider interacting particles, we need to specify and calculate their interactions.- The latter, while not strictly necessary, is often used in order to approximate infinite systems and avoid boundary effects in our results.- The idea of periodic boundaries is quite simple, imagine that your finite-sized simulation box (colored in orange in the figure) is periodically repeated throughout space.- The unit cell is usually a cubic box, which can be used to generate a perfect three-dimensional tiling that fills all of space, as illustrated in the figure.- The orange unit cell, is surrounded by infinitely many blue periodic copies of itself. - Note that all the boxes are exactly the same, which one you call the unit cell does not matter. - In addition, since they are all the same, we only need to consider what happens in one of them. - We just have to be careful to properly handle particle exchanges between cells.- When a particle, shown here in blue as an example, passes through one of the boundaries of the unit cell, it enters the neighboring cell. - This means that it re-appears on the opposite side of the unit cell with the same velocity.  Inter-particle interaction![](fig/soft_core.png)$$U=\sum_{r_{ij}<2\sigma}\phi(r_{ij})=\sum_{r_{ij}<2\sigma}4\epsilon\left(\frac{\sigma}{r_{ij}} \right)^{12},\ \ \ \ \mathbf{F}_{i}=-\frac{dU}{d \mathbf{r}_{ij}}=\sum_{r_{ij}<2\sigma}48\epsilon\left(\frac{\sigma}{r_{ij}} \right)^{12} \frac{\mathbf{r}_{ij}}{r^2_{ij}}\tag{I1, I2}$$$$\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i,\ \ \ r_{ij}=|\mathbf{r}_{ij}|$$  Note 3- To introduce inter-particle interactions, a well defined potential function must be defined for the system.- Here we assume a purely repulsive soft-core potential to represent the excluded volume of each Brownian particle with a diameter equal to $\sigma$.- The total potential energy of the system is defined by Eq.(I1), where the sum is taken over particle pairs separated by a distance less than 2 $\sigma$ only.- The inter-particle force acting the i-th particle is then given by Eq.(I2).  Simulation code for interacting Brownian particles with animation Import libraries<jupyter_code>% matplotlib nbagg
import numpy as np # import numpy library as np
import matplotlib.pyplot as plt # import pyplot library as plt 
import matplotlib.animation as animation # import animation modules from matplotlib
from   mpl_toolkits.mplot3d import Axes3D # import Axes3D from mpl_toolkits.mplot3d
from   numpy import newaxis # import newaxis used for inter-particle force calculations
plt.style.use('ggplot') # use "ggplot" style for graphs<jupyter_output><empty_output><jupyter_text>Note 4- As always, we begin by importing the necessary numerical and plotting libraries.- Compared to the previous code example, we import an additional module `newaxis` from the `numpy` library.- The use of the `newaxis` module allows us to compute inter-particle interactions efficiently in Python by avoiding unnecessary for loops.  Define `init` function for `FuncAnimation`<jupyter_code>def init():
    global R,V,W,F,Rs,Vs,Ws,time # define global variables
    initconf()          # create random particle configuration without overlaps
    V[:,:] = 0.0        # initialize all the variables to zero
    W[:,:] = 0.0        # initialize all the variables to zero
    F[:,:] = 0.0        # initialize all the variables to zero
    Rs[:,:,:] = 0.0     # initialize all the variables to zero
    Rs[0,:,:] = R[:,:]  # store initial particle positions in Rs
    Vs[:,:,:] = 0.0     # initialize all the variables to zero
    Ws[:,:,:] = 0.0     # initialize all the variables to zero
    time[:]   = 0.0     # initialize all the variables to zero
    title.set_text(r'') # empty title
    line.set_data([],[]) # set line data to show the trajectory of particle n in 2d (x,y)
    line.set_3d_properties([]) # add z-data separately for 3d plot
    particles.set_data([],[]) # set position current (x,y) position data for all particles
    particles.set_3d_properties([]) # add current z data of particles to get 3d plot
    return particles,title,line # return listed objects that will be drawn by FuncAnimation<jupyter_output><empty_output><jupyter_text>Note 5- Compared to the previous code example, a few modifications have been done for the `init` function.- The array F, which will be used to store the current values of the inter-particle force, is added to the list of global variables and initialized by setting all the elements to zero.- Because particles should not overlap with each other due to the inter-particle interaction, we prepare the initial particle configuration using uniform random numbers by calling a function `initconf` that will be define later.- The remaining parts of the initialization code are unchanged.  Define `animate` function for `FuncAnimation`<jupyter_code>def animate(i):
    global R,V,W,F,Rs,Vs,Ws,time # define global variables
    time[i]=i*dt # store time in each step in an array time
    particleforces() # compute inter-particle force F by examining all nump^2 pairs
    W = std*np.random.randn(nump,dim) # generate an array of random forces accordingly to Eqs.(F10) and (F11)
    V = V*(1-zeta/m*dt)+F/m*dt+W/m # update velocity via Eq.(F9) with inter-particle force F
    R = R+V*dt # update position via Eq.(F5)
    Rs[i,:,:]=R # accumulate particle positions at each step in an array Rs
    Vs[i,:,:]=V # accumulate particle velocitys at each step in an array Vs
    Ws[i,:,:]=W # accumulate random forces at each step in an array Ws
    title.set_text(r"t = "+str(time[i])+"/"+str((nums-1)*dt))  # set the title to display the current time
    line.set_data(Rs[:i+1,n,0],Rs[:i+1,n,1]) # set the line in 2D (x,y)
    line.set_3d_properties(Rs[:i+1,n,2]) # add z axis to set the line in 3D
    particles.set_data(pbc(R[:,0],box[0]), pbc(R[:,1],box[1])) # set the current position of all the particles in 2d (x,y) with PBC
    particles.set_3d_properties(pbc(R[:,2],box[2])) # add z axis to set the particle in 3D with PBC
    return particles,title,line # return listed objects that will be drawn by FuncAnimation<jupyter_output><empty_output><jupyter_text>Note 6- In the `animate` function, we compute the inter-particle force F on all particles by calling the `particleforce` function in the 5-th line. - This function will be defined later.- In the 6-th line, we modified the velocity update operation to include the inter-particle force F just computed above.- In the 14th and 15th lines, the object `particles`, which has to be drawn in the animation, is updated to use the current particle positions under periodic boundary conditions. - This is done by calling the `pbc` function on the x,y,and z particle positions.- You might wonder why we don't use the pbc function directly when updating the positions, in line 7 of the code. - The reason for this is to simplify the analysis that we will perform later. - In particular, calculating the mean-squared displacements becomes more difficult if we use the periodically wrapped positions.- As written, our code keeps track of the particles that were initially in the center unit cell. - While the positions stored in R  and Rs might go outside of the center cell, if we know the position of one of the copies, we know the position of all of the copies.- We just have to be careful when computing interactions between particles, to make sure that we properly take this into account.  Newly defined functions for interacting Brownian particles<jupyter_code>def pbc(r, lbox): # enforce Periodic Boundary Conditions for all positions
    return np.fmod(r+lbox,lbox)
def distance(r1,r2,lbox): # Compute distance vector R2 - R1 with PBC
    return r2-r1-np.around((r2-r1)/lbox)*lbox
def fij(r2,rij): # calculate Fij=dU/drij
#    f=-24*eps*(2*(r2/sig**2)**(-6)-(r2/sig**2)**(-3))/r2*rij # Lennard-Jones potential
    f=-48*eps*((r2/sig**2)**(-6))/r2*rij # soft-core potential
    return f 
def particleforces(): # compute inter-particle force F by examining all nump^2 pairs
    global F
    F[:,:] = 0.0
    for n in range(nump): # repeat below for all particles
        rij = distance(R[n,:], R, box) # distance vectors rij=R_i-R_j for all i (1 <= i <= nump)
        r2  = np.linalg.norm(rij, axis=1)**2 # square distance rij**2
        nei = (r2 < (2.0*sig)**2)  # list neighbor particles of j
        nei[n] = False # ignore self pair (i=j)
        F[n,:] = np.sum(fij(r2[nei, newaxis], rij[nei,:]), axis=0) # total force on particle j
def initconf():  # create random particle configuration without overlapping
    global R,V,W,F,Rs,Vs,Ws,time
    for n in range(nump): # repeat below from n=0 to nump-1
        nn=0 # set overlap true to perform while loop below for the n-th particle 
        while nn == 0: # repeat the loop below while overlap is true (nn==0)
            R[n,:]=np.random.rand(dim)*box # generate a position candidate for n using uniform random number.
            nn = 1 # initialize overlap as false
            for l in range(n): # examine overlap generated positions (from l=0 to n-1)
                rij = distance(R[n,:],R[l,:],box) # calculate distance vector rij=R_l-R_n
                r2  = np.linalg.norm(rij)**2 # calculate the squared distance rij**2
                if r2 < (0.90*sig)**2: # check if the distance is smaller than threshold
                # Yes -> perform below (nn=0) -> repeat while loop, No (nn=1) -> exit while loop
                    nn = 0 # set overlap true<jupyter_output><empty_output><jupyter_text>Note 7- To modify our previous code for the present case of interacting Brownian particles, we have defined 5 new functions as shown here.- The `pbc` function applies the periodic boundary conditions to a vector 'r', assuming the size of the unit cell is 'box'. - If any of the components of 'r' lie outside the unit cell, they are wrapped back to be inside the proper domain.- As written here, the function works if r is a vector or a matrix (vector of vectors).- The `distance` function computes the distance vector from r1 to r2.- Given python's advanced array capabilities, the same function can be used to compute the distance between two particles, or between one particle and a set of particles.- In the first case r1 and r2 are both 3d vectors, and the result is the vector r12. - In the second case r1 is a 3d vector and r2 would be a $n\times3$ matrix, and the result is a matrix, with row 'i' containing the distance vector between particle 1 and particle 'i'.- If you want to write efficient python code, you should learn how these array broadcasting rules work.- More information can be found online at the SciPy.org website [https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html].- The `fij` function together with the `particleforce` function computes the inter-particle force F acting on each particle by examining all nump^2 particle pairs.- The `initconf` function is used to create a random particle configuration using uniform random numbers.- To avoid the occurrence of particle overlap, the distance between all existing particles is examined. If an overlap exists, the particle position is rejected and a new random position is tested.  Set parameters and initialize variables<jupyter_code>dim  = 3    # system dimension (x,y,z)
nump = 100  # number of interacting Brownian particles to simulate 
nums = 4096 # number of simulation steps
dt   = 0.01 # set time increment, \Delta t
zeta = 1.0  # set friction constant, \zeta
m    = 1.0  # set particle mass, m
kBT  = 1.0  # set thermal energy, k_B T
std  = np.sqrt(2*kBT*zeta*dt) # calculate std for \Delta W via Eq.(F11)
sig  = 1.0  # unit of length of inter-particle potential
eps  = 1.0  # unit of energy inter-particle potential
vf   = 0.001 ##volume fraction of particles < 0.45 
boxl = np.power(nump*np.pi/6/vf,1/3) # calculate the side length of unit cell
print('Volume fraction =',vf,'  boxl =',boxl) # print vf and boxl
box  = np.array([boxl,boxl,boxl])*sig # set array box[dim]
np.random.seed(0) # initialize random number generator with a seed=0
R = np.zeros([nump,dim]) # array to store current positions and set initial condition Eq.(F12)
V = np.zeros([nump,dim]) # array to store current velocities and set initial condition Eq.(F12)
W = np.zeros([nump,dim]) # array to store current random forcces
F = np.zeros([nump,dim]) # rray to store current particle orcces
Rs = np.zeros([nums,nump,dim]) # array to store positions at all steps
Vs = np.zeros([nums,nump,dim]) # array to store velocities at all steps
Ws = np.zeros([nums,nump,dim]) # array to store random forces at all steps
time = np.zeros([nums]) # an array to store time at all steps<jupyter_output>Volume fraction = 0.001   boxl = 37.4110192682<jupyter_text>Note 8- Here we set the system parameters and initialize all the particle variables.- Compared to the previous code example, we introduced three new parameters `sig`, `eps`, and `vf`. - `sig` and `eps` are the units of length and energy appearing in Eq.(I1), which we set to be both equal to one.- `vf` represents the volume fraction of the particles $=\pi\sigma^3 / 6V_{cell}$.- Using the parameters defined here, we perform a 4096 step simulation for 100 interacting Brownian particles with diameter one, in a cubic box of side length 37 $\sigma$, which gives rise to a particle volume fraction of 0.001.- This means that only $0.1\%$ of the space in the unit cell is occupied by particles.  Perform and animate the simulation using `FuncAnimation`<jupyter_code>fig = plt.figure(figsize=(10,10)) # set fig with its size 10 x 10 inch
ax = fig.add_subplot(111,projection='3d') # creates an additional axis to the standard 2D axes
ax.set_xlim(0.0,box[0]) # set x-range
ax.set_ylim(0.0,box[1]) # set y-range
ax.set_zlim(0.0,box[2]) # set z-range
ax.set_xlabel(r"x",fontsize=20) # set x-lavel
ax.set_ylabel(r"y",fontsize=20) # set y-lavel
ax.set_zlabel(r"z",fontsize=20) # set z-lavel
ax.view_init(elev=12,azim=120)  # set view point
particles, = ax.plot([],[],[],linestyle='None',color='r',marker='o',ms=250/box[0],alpha=0.5) # define object particles
title = ax.text(0.,0.,0.,r'',transform=ax.transAxes,va='center') # define object title
line, = ax.plot([],[],[],'b',lw=2,alpha=0.8) # define object line
n = 0  # trajectry line is plotted for the n-th particle
anim = animation.FuncAnimation(fig,func=animate,init_func=init,
            frames=nums,interval=5,blit=True,repeat=False)
## If you have ffmpeg installed on your machine 
## you can save the animation by uncomment the last line
## You may install ffmpeg by typing the following command in command prompt
## conda install -c menpo ffmpeg
## 
#anim.save('movie.mp4',fps=20,dpi=400)<jupyter_output><empty_output><jupyter_text>Note 9- Now we can run the simulation, visualize the results, and store the trajectory data.- On a standard pc it may take a few minutes until the simulation finishes, but please be patient to complete it.- The trajectory data generated here is needed in the next step.  The mean square displacement and the diffusion constant<jupyter_code># mean square displacement vs time
msd = np.zeros([nums])
for i in range(nums):
    for n in range(nump):
        msd[i]=msd[i]+np.linalg.norm(Rs[i,n,:]-Rs[0,n,:])**2 # (R(t) - R(0))^2
    msd[i] = msd[i]/nump
dmsd = np.trapz(msd, dx=dt)/(3*(nums*dt)**2)
print('D =',kBT/zeta,'(Theoretical)')
print('D =',dmsd,'(Simulation via MSD)')
print('Volume fraction =',vf) ### print vf
fig, ax = plt.subplots(figsize=(7.5,7.5))
ax.set_xlabel(r"$t$", fontsize=20)
ax.set_ylabel(r"mean square displacement", fontsize=16)
ax.plot(time,6*kBT/zeta*time,'r',lw=6,label=r'$6Dt={6k_BT t}/{\zeta}$')
ax.plot(time,msd,'b',lw=4, label=r'$\langle R^2(t)\rangle$')
ax.legend(fontsize=16, loc=4)
plt.show()<jupyter_output>D = 1.0 (Theoretical)
D = 0.885552952684 (Simulation via MSD)
Volume fraction = 0.001
{'fraction_non_alphanumeric': 0.0993154726779181, 'fraction_numerical': 0.01917856721350171, 'mean_word_length': 4.259776536312849, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 30, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '597916', 'license': 'MIT', 'n_tokens_mistral': 5147, 'n_tokens_neox': 4856, 'n_words': 2178, 'path': 'edx-stochastic-data-analysis/downloaded_files/05/.ipynb_checkpoints/SP053-checkpoint.ipynb', 'repo_name': 'mirandagil/extra-courses', 'size': '16946'}