stochastic_processes_data_analysis/10_Interacting_Brownian_Particles.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Interacting Brownian particles\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Necessary changes for interacting Brownian particles: periodic boundary conditions\n",
    "\n",
    "![](images/PBC.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Inter-particle interaction: repulsive soft-core potential\n",
    "\n",
    "![](images/soft_core.png)\n",
    "\n",
    "$$\n",
    "U=\\sum_{r_{ij}<2\\sigma}\\phi(r_{ij})=\\sum_{r_{ij}<2\\sigma}4\\epsilon\\left(\\frac{\\sigma}{r_{ij}} \\right)^{12},\\ \\ \\ \\ \n",
    "\\mathbf{F}_{i}=-\\frac{\\partial U}{\\partial \\mathbf{r}_{ij}}\n",
    "=\\sum_{r_{ij}<2\\sigma}48\\epsilon\\left(\\frac{\\sigma}{r_{ij}} \\right)^{12} \\frac{\\mathbf{r}_{ij}}{r^2_{ij}}\\tag{I1, I2}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{r}_{ij} = \\mathbf{r}_j - \\mathbf{r}_i,\\ \\ \\ r_{ij}=|\\mathbf{r}_{ij}|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simulation code for interacting Brownian particles with animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-03-08T11:38:24.873975",
     "start_time": "2017-03-08T11:38:24.865189"
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "%matplotlib nbagg\n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt \n",
    "import matplotlib.animation as animation # import animation modules from matplotlib\n",
    "from   mpl_toolkits.mplot3d import Axes3D # import Axes3D from mpl_toolkits.mplot3d\n",
    "from   numpy import newaxis # import newaxis used for inter-particle force calculations\n",
    "plt.style.use('ggplot') # use \"ggplot\" style for graphs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Define `init` function for `FuncAnimation`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-03-08T11:33:20.948957",
     "start_time": "2017-03-08T11:33:20.940296"
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def init():\n",
    "    global R,V,W,F,Rs,Vs,Ws,time # define global variables\n",
    "    initconf()          # create random particle configuration without overlaps\n",
    "    V[:,:] = 0.0        # initialize all the variables to zero\n",
    "    W[:,:] = 0.0        # initialize all the variables to zero\n",
    "    F[:,:] = 0.0        # initialize all the variables to zero\n",
    "    Rs[:,:,:] = 0.0     # initialize all the variables to zero\n",
    "    Rs[0,:,:] = R[:,:]  # store initial particle positions in Rs\n",
    "    Vs[:,:,:] = 0.0     # initialize all the variables to zero\n",
    "    Ws[:,:,:] = 0.0     # initialize all the variables to zero\n",
    "    time[:]   = 0.0     # initialize all the variables to zero\n",
    "    title.set_text(r'') # empty title\n",
    "    line.set_data([],[]) # set line data to show the trajectory of particle n in 2d (x,y)\n",
    "    line.set_3d_properties([]) # add z-data separately for 3d plot\n",
    "    particles.set_data([],[]) # set position current (x,y) position data for all particles\n",
    "    particles.set_3d_properties([]) # add current z data of particles to get 3d plot\n",
    "    return particles,title,line # return listed objects that will be drawn by FuncAnimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Define `animate` function for `FuncAnimation`\n",
    "\n",
    "The Euler method used in previous `animate` function is replaced with the Leap-Frog method in the following new `animate` function because of the numerical instability taking place with the former method for the present case of interacting Brownian particles. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-03-08T11:33:23.887393",
     "start_time": "2017-03-08T11:33:23.876963"
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def animate(i):\n",
    "    global R,V,W,F,Rs,Vs,Ws,time # define global variables\n",
    "    time[i] = i*dt # store time in each step in an array time\n",
    "    particleforces() # compute inter-particle force F by examining all nump^2 pairs\n",
    "    W = std*np.random.randn(nump,dim) # generate an array of random forces accordingly to Eqs.(F10) and (F11)\n",
    "    \n",
    "#    V = V*(1-zeta/m*dt)+F/m*dt+W/m # update velocity via Eq.(I11) using the L-F method\n",
    "\n",
    "    V = (V*(1-zeta/2/m*dt) + F/m*dt+W/m) / (1+zeta/2/m*dt) # update velocity via Eq.(I11) using the L-F method\n",
    "    R = R + V*dt # update position via Eq.(I12) using the L-F method\n",
    "    Rs[i,:,:] = R # accumulate particle positions at each step in an array Rs\n",
    "    Vs[i,:,:] = V # accumulate particle velocitys at each step in an array Vs\n",
    "    Ws[i,:,:] = W # accumulate random forces at each step in an array Ws\n",
    "    title.set_text(r\"t = \"+str(time[i])+\"/\"+str((nums-1)*dt))  # set the title to display the current time\n",
    "    line.set_data(Rs[:i+1,n,0], Rs[:i+1,n,1]) # set the line in 2D (x,y)\n",
    "    line.set_3d_properties(Rs[:i+1,n,2]) # add z axis to set the line in 3D\n",
    "    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\n",
    "    particles.set_3d_properties(pbc(R[:,2], box[2])) # add z axis to set the particle in 3D with PBC\n",
    "    return particles,title,line # return listed objects that will be drawn by FuncAnimation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Newly defined functions for interacting Brownian particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-03-08T11:33:29.453575",
     "start_time": "2017-03-08T11:33:29.433781"
    },
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def pbc(r, lbox): # enforce Periodic Boundary Conditions for all positions\n",
    "    return np.fmod(r + lbox,lbox)\n",
    "\n",
    "def distance(r1,r2,lbox): # Compute distance vector R2 - R1 with PBC\n",
    "    return r2-r1-np.around((r2-r1)/lbox)*lbox\n",
    "\n",
    "def fij(r2,rij): # calculate Fij=dU/drij\n",
    "#    f=-24*eps*(2*(r2/sig**2)**(-6)-(r2/sig**2)**(-3))/r2*rij # Lennard-Jones potential\n",
    "    f = -48*eps*((r2/sig**2)**(-6))/r2*rij # soft-core potential\n",
    "    return f \n",
    "\n",
    "def particleforces(): # compute inter-particle force F by examining all nump^2 pairs\n",
    "    global F\n",
    "    F[:,:] = 0.0\n",
    "    for n in range(nump): # repeat below for all particles\n",
    "        rij = distance(R[n,:], R, box) # distance vectors rij=R_i-R_j for all i (1 <= i <= nump)\n",
    "        r2  = np.linalg.norm(rij, axis = 1)**2 # square distance rij**2\n",
    "        nei = (r2 < (2.0*sig)**2)  # list neighbor particles of j\n",
    "        nei[n] = False # ignore self pair (i=j)\n",
    "        F[n,:] = np.sum(fij(r2[nei, newaxis], rij[nei,:]), axis=0) # total force on particle j\n",
    "        \n",
    "def initconf():  # create random particle configuration without overlapping\n",
    "    global R,V,W,F,Rs,Vs,Ws,time\n",
    "    for n in range(nump): # repeat below from n=0 to nump-1\n",
    "        nn = 0 # set overlap true to perform while loop below for the n-th particle \n",
    "        while nn == 0: # repeat the loop below while overlap is true (nn==0)\n",
    "            R[n,:] = np.random.rand(dim)*box # generate a position candidate for n using uniform random number.\n",
    "            nn = 1 # initialize overlap as false\n",
    "            for l in range(n): # examine overlap generated positions (from l=0 to n-1)\n",
    "                rij = distance(R[n,:], R[l,:], box) # calculate distance vector rij=R_l-R_n\n",
    "                r2  = np.linalg.norm(rij)**2 # calculate the squared distance rij**2\n",
    "                if r2 < (0.90*sig)**2: # check if the distance is smaller than threshold\n",
    "                # Yes -> perform below (nn=0) -> repeat while loop, No (nn=1) -> exit while loop\n",
    "                    nn = 0 # set overlap true \n",
    "                    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Set parameters and initialize variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-03-08T12:32:34.139411",
     "start_time": "2017-03-08T12:32:34.125859"
    },
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dim  = 3    # system dimension (x,y,z)\n",
    "nump = 100  # number of interacting Brownian particles to simulate \n",
    "nums = 4096 # number of simulation steps\n",
    "dt   = 0.01 # set time increment, \\Delta t\n",
    "zeta = 1.0  # set friction constant, \\zeta\n",
    "m    = 1.0  # set particle mass, m\n",
    "kBT  = 1.0  # set thermal energy, k_B T\n",
    "std  = np.sqrt(2*kBT*zeta*dt) # calculate std for \\Delta W via Eq.(F11)\n",
    "sig  = 1.0  # unit of length of inter-particle potential\n",
    "eps  = 1.0  # unit of energy inter-particle potential\n",
    "vf   = 0.2 # volume fraction of particles < 0.45 \n",
    "boxl = np.power(nump*np.pi/6/vf,1/3) # calculate the side length of unit cell\n",
    "print('Volume fraction =',vf,'  boxl =',boxl) # print vf and boxl\n",
    "box  = np.array([boxl, boxl, boxl])*sig # set array box[dim]\n",
    "np.random.seed(0) # initialize random number generator with a seed=0\n",
    "R = np.zeros([nump,dim]) # array to store current positions and set initial condition Eq.(F12)\n",
    "V = np.zeros([nump,dim]) # array to store current velocities and set initial condition Eq.(F12)\n",
    "W = np.zeros([nump,dim]) # array to store current random forcces\n",
    "F = np.zeros([nump,dim]) # rray to store current particle orcces\n",
    "Rs = np.zeros([nums,nump,dim]) # array to store positions at all steps\n",
    "Vs = np.zeros([nums,nump,dim]) # array to store velocities at all steps\n",
    "Ws = np.zeros([nums,nump,dim]) # array to store random forces at all steps\n",
    "time = np.zeros([nums]) # an array to store time at all steps\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Perform and animate the simulation using `FuncAnimation` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-03-08T12:36:39.814630",
     "start_time": "2017-03-08T12:36:39.774689"
    },
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize = (10,10)) # set fig with its size 10 x 10 inch\n",
    "ax = fig.add_subplot(111, projection = '3d') # creates an additional axis to the standard 2D axes\n",
    "ax.set_xlim(0.0, box[0]) # set x-range\n",
    "ax.set_ylim(0.0, box[1]) # set y-range\n",
    "ax.set_zlim(0.0, box[2]) # set z-range\n",
    "ax.set_xlabel(r\"x\", fontsize = 20) # set x-lavel\n",
    "ax.set_ylabel(r\"y\", fontsize = 20) # set y-lavel\n",
    "ax.set_zlabel(r\"z\", fontsize = 20) # set z-lavel\n",
    "ax.view_init(elev = 12, azim = 120)  # set view point\n",
    "particles, = ax.plot([], [], [], linestyle = 'None', color = 'r', marker = 'o', ms = 250/box[0], alpha = 0.5) # define object particles\n",
    "title = ax.text(0., 0., 0., r'', transform = ax.transAxes, va = 'center') # define object title\n",
    "line, = ax.plot([], [], [], 'b', lw = 2, alpha = 0.8) # define object line\n",
    "n = 0  # trajectry line is plotted for the n-th particle\n",
    "anim = animation.FuncAnimation(fig, func = animate, init_func = init,\n",
    "            frames = nums,interval = 5, blit = True, repeat = False)\n",
    "\n",
    "## If you have ffmpeg installed on your machine \n",
    "## you can save the animation by uncomment the last line\n",
    "## You may install ffmpeg by typing the following command in command prompt\n",
    "## conda install -c menpo ffmpeg\n",
    "## \n",
    "#anim.save('movie.mp4',fps=20,dpi=400)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The mean square displacement and the diffusion constant "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-03-08T12:44:23.453755",
     "start_time": "2017-03-08T12:44:21.408854"
    }
   },
   "outputs": [],
   "source": [
    "# mean square displacement vs time\n",
    "msd = np.zeros([nums])\n",
    "\n",
    "for i in range(nums):\n",
    "    for n in range(nump):\n",
    "        msd[i] = msd[i] + np.linalg.norm(Rs[i,n,:] - Rs[0,n,:])**2 # (R(t) - R(0))^2\n",
    "    msd[i] = msd[i]/nump\n",
    "    \n",
    "dmsd = np.trapz(msd, dx = dt)/(3*(nums*dt)**2)\n",
    "print('D_0 =', kBT/zeta, '(Theoretical)')\n",
    "print('D =', dmsd, '(Simulation via MSD)')\n",
    "print('Volume fraction =', vf) ### print vf\n",
    "fig, ax = plt.subplots(figsize = (7.5,7.5))\n",
    "ax.set_xlabel(r\"$t$\", fontsize =20)\n",
    "ax.set_ylabel(r\"mean square displacement\", fontsize = 16)\n",
    "ax.plot(time, 6*kBT/zeta*time, 'r', lw = 6, label = r'$6Dt={6k_BT t}/{\\zeta}$')\n",
    "ax.plot(time, msd, 'b', lw = 4, label = r'$\\langle R^2(t)\\rangle$')\n",
    "ax.legend(fontsize = 16, loc = 4)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Velocity auto-correlation function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute self-correlation of vector v\n",
    "\n",
    "def auto_correlate(v):\n",
    "    # np.correlate computes C_{v}[k] = sum_n v[n+k] * v[n]\n",
    "    corr = np.correlate(v, v, mode = \"full\") # correlate returns even array [0:2*nums-1] centered at nums-1\n",
    "    return corr[len(v)-1:]/len(v) # take positive values and normalize by number of points\n",
    "\n",
    "corr = np.zeros([nums])\n",
    "for n in range(nump):\n",
    "    for d in range(dim):\n",
    "        corr = corr + auto_correlate(Vs[:,n,d]) # correlation of d-component of velocity for particle n\n",
    "corr = corr/nump #average over all particles\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (7.5,7.5))\n",
    "ax.plot(time, dim*kBT/m*np.exp(-zeta/m*time), 'r', lw = 6, label = r'$(3k_BT/m)\\exp(-\\zeta t/m)$')\n",
    "ax.plot(time, corr, 'b', lw = 3, label = r'$\\langle\\mathbf{V}(t)\\cdot \\mathbf{V}(0)\\rangle$')\n",
    "ax.set_xlabel(r\"tsize = 20\")\n",
    "ax.set_ylabel(r\"auto correlation $\\varphi_V(t)$\", fontsize = 20)\n",
    "ax.set_xlim(-1,10)\n",
    "ax.set_ylim(-0.5,3.5)\n",
    "ax.legend(fontsize = 16)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power spectrum of particle velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# return power spectrum for positive frequencies of even signal v\n",
    "from numpy import fft\n",
    "\n",
    "def psd(v,dt):  \n",
    "    vw = fft.fft(v)*dt # V(w) with zero-frequency component at vw(0)\n",
    "    return np.abs(vw[:nums//2])**2/(nums*dt) # S_V for w > 0 (Eq. (G9))\n",
    "\n",
    "Sw   = np.zeros([nums//2])\n",
    "for n in range(nump):\n",
    "    for d in range(dim):\n",
    "        Sw = Sw + psd(Vs[:,n,d],dt) # power spectrum of d-component of velocity for particle n\n",
    "Sw = Sw/nump\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (7.5,7.5))\n",
    "gamma = zeta/m\n",
    "omega = fft.fftfreq(nums, d = dt)[:nums//2]*2.0*np.pi\n",
    "\n",
    "ax.plot(omega, (6.0*kBT/m)*gamma/(omega**2 + gamma**2), 'r', lw = 6, label = r'$6k_BT\\gamma_c/(m(\\omega^2+\\gamma_c^2))$')\n",
    "ax.plot(omega, Sw, 'b', lw = 3, label = r'$|\\mathbf{V}_\\tau(\\omega)|^2 / \\tau$')\n",
    "\n",
    "ax.set_xlabel(r\"angular frequency $\\omega$\", fontsize = 16)\n",
    "ax.set_ylabel(r\"spectral density $S_V(\\omega)$\", fontsize = 16)\n",
    "ax.legend(fontsize = 16)\n",
    "plt.xlim(-1, 10)\n",
    "plt.ylim(-1, 8)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "# References\n",
    "\n",
    "- The SciPy.org website, https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html\n",
    "\n",
    "- Listed in the major unsolved problems in physics, https://en.wikipedia.org/wiki/List_of_unsolved_problems_in_physics#Condensed_matter_physics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of simulation method\n",
    "\n",
    "### Original differential equation\n",
    "$i$: particle index\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{d\\mathbf{R}_i(t)}{dt}=\\mathbf{V}_i(t)\\tag{I3}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "m\\frac{d\\mathbf{V}_i(t)}{dt}=\\color{black}{-\\zeta\\mathbf{V}_i(t)}+\\color{black}{\\mathbf{F}_i(t)} +\\color{black}{\\mathbf{F}_R(t)} \n",
    "\\tag{I4}\n",
    "\\end{equation}\n",
    "$\\hspace{80mm}$with\n",
    "\\begin{equation}\n",
    "\\langle \\mathbf{F}_R(t)\\rangle=\\mathbf{0} \n",
    "\\tag{I5}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle \\mathbf{F}_R(t)\\mathbf{F}_R(0)\\rangle = {2k_B T\\zeta}\\mathbf{I}\\delta(t)  \n",
    "\\tag{I6}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Euler method\n",
    "$i$: index of particles,\n",
    "$n$: time step\n",
    "$$\n",
    "\\mathbf{V}_{i,n+1}\n",
    "=\\left(1-\\frac{\\zeta}{m}\\Delta t\\right)\\mathbf{V}_{i,n}  + \\frac{1}{m}\\mathbf{F}_{i,n}\\Delta t + \\frac{1}{m} {\\Delta \\mathbf{W}_n} \n",
    "\\tag{I7}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{R}_{i,n+1}=\\mathbf{R}_{i,n}+\\mathbf{V}_{i,n} \\Delta t \\hspace{15mm}\\tag{I8}\n",
    "$$\n",
    "$\\hspace{80mm}$with\n",
    "\\begin{equation}\n",
    "\\langle \\Delta \\mathbf{W}_n\\rangle=\\mathbf{0} \n",
    "\\tag{I9}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle \\Delta \\mathbf{W}_n\\Delta \\mathbf{W}_l\\rangle = {2k_B T\\zeta}\\Delta t\\mathbf{I}\\delta_{nl}  \n",
    "\\tag{I10}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Leap-Frog method\n",
    "$i$: index of particles,\n",
    "$n$: time step\n",
    "$$\n",
    "\\mathbf{V}_{i,n+\\frac{1}{2}}\n",
    "=\\mathbf{V}_{i,n-\\frac{1}{2}}-\\frac{\\zeta}{m}\\Delta t\\mathbf{V}_{i,n}  + \\frac{1}{m}\\mathbf{F}_{i,n}\\Delta t + \\frac{1}{m} {\\Delta \\mathbf{W}_n} \\\\\n",
    "\\hspace{45mm}=\\mathbf{V}_{i,n-\\frac{1}{2}}-\\frac{\\zeta}{m}\\Delta t\\frac{1}{2}\\left(\\mathbf{V}_{i,n+\\frac{1}{2}}+\\mathbf{V}_{i,n-\\frac{1}{2}}\\right)  + \\frac{1}{m}\\mathbf{F}_{i,n}\\Delta t + \\frac{1}{m} {\\Delta \\mathbf{W}_n}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\therefore\\ \\ \\mathbf{V}_{i,n+\\frac{1}{2}}\n",
    "=\\left(\\left(1-\\frac{\\zeta}{2m}\\Delta t\\right)\\mathbf{V}_{i,n-\\frac{1}{2}} + \\frac{1}{m}\\mathbf{F}_{i,n}\\Delta t + \\frac{1}{m} {\\Delta \\mathbf{W}_n}\\right)\\left({1+\\frac{\\zeta}{2m}\\Delta t}\\right)^{-1} \\tag{I11} \n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{R}_{i,n+1}=\\mathbf{R}_{i,n}+\\mathbf{V}_{i,n+\\frac{1}{2}} \\Delta t \\hspace{15mm}\\tag{I12}\n",
    "$$\n",
    "$\\hspace{80mm}$with\n",
    "\\begin{equation}\n",
    "\\langle \\Delta \\mathbf{W}_n\\rangle=\\mathbf{0} \n",
    "\\tag{I13}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle \\Delta \\mathbf{W}_n\\Delta \\mathbf{W}_l\\rangle = {2k_B T\\zeta}\\Delta t\\mathbf{I}\\delta_{nl}  \n",
    "\\tag{I14}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2nd order Runge-Kutta method\n",
    "$i$: index of particles,\n",
    "$n$: time step\n",
    "\n",
    "$$\n",
    "\\mathbf{V}'_{i,n+\\frac{1}{2}}\n",
    "=\\mathbf{V}_{i,n}-\\frac{\\zeta}{m}\\frac{\\Delta t}{2}\\mathbf{V}_{i,n}  + \\frac{1}{m}\\mathbf{F}_{i,n}\\frac{\\Delta t}{2} \n",
    "=\\left(1-\\frac{\\zeta}{m}\\frac{\\Delta t}{2}\\right)\\mathbf{V}_{i,n}  \n",
    "+ \\frac{1}{m}\\mathbf{F}_{i,n}\\frac{\\Delta t}{2} \\tag{I15}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{R}'_{i,n+\\frac{1}{2}}=\\mathbf{R}_{i,n}+\\mathbf{V}_{i,n} \\frac{\\Delta t}{2} \\hspace{15mm}\\tag{I16}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{V}_{n,i+1}\n",
    "=\\mathbf{V}_{i,n}-\\frac{\\zeta}{m}\\Delta t\\mathbf{V}'_{n,i+\\frac{1}{2}} + \\frac{1}{m}\\Delta t\\mathbf{F}_{i,n+\\frac{1}{2}} + \\frac{1}{m} {\\Delta \\mathbf{W}_i} \n",
    "\\tag{I17}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbf{R}_{n,i+1}=\\mathbf{R}_{n,i}+\\mathbf{V}'_{n,i+\\frac{1}{2}} \\Delta t \\hspace{15mm}\\tag{I18}\n",
    "$$\n",
    "\n",
    "$\\hspace{80mm}$with\n",
    "\\begin{equation}\n",
    "\\langle \\Delta \\mathbf{W}_n\\rangle=\\mathbf{0} \n",
    "\\tag{I19}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\langle \\Delta \\mathbf{W}_n\\Delta \\mathbf{W}_l\\rangle = {2k_B T\\zeta}\\Delta t\\mathbf{I}\\delta_{nl}  \n",
    "\\tag{I20}\n",
    "\\end{equation}\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "457px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}