Brownian Dynamics Simulation using Periodic Boundary
Condition and a Spatially-Periodic Potential
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

This example was inspired from the Matlab program available at the following MIT OpenCourseWare web page:

http://ocw.mit.edu/OcwWeb/Chemical-Engineering/10-34Fall-2005/CourseHome/index.htm

Wait until system equilibrates then perform Brownian Simutation (this computation takes several minutes)

<< Statistics`ContinuousDistributions`

Di = 1;

Δt = 0.001;

Vf = 0;

Ea = 1;

\(For[j = 1, 
      j < 1001, {i = 1, 
        While[i < 1501, {x[0] = 0, 
            x[i] = x[
                  i - 1] + \((Vf - 
                      2\ Pi\ Ea\ Sin[x[i - 1]\ Pi]\ Cos[
                          x[i - 1]\ Pi])\)\ Δt + \(\@\(2\ Di\ Δt\)\) Random[NormalDistribution[0, 1]], 
            If[Mod[x[i], 1] ≠ 1, x[i] = Mod[x[i], 1]], i = i + 1}], 
        For[k = 1, 
          k < 1001, {x[k - 1] = x[i - 1], 
            x[k] = x[
                  k - 1] + \((Vf - 
                      2\ Pi\ Ea\ Sin[x[k - 1]\ Pi]\ Cos[
                          x[k - 1]\ Pi])\)\ Δt + \(\@\(2\ Di\ Δt\)\) Random[NormalDistribution[0, 1]], 
            If[Mod[x[k], 1] ≠ 1, x[k] = Mod[x[k], 1]], 
            X[j, k] = x[k], k = k + 1}], 
        If[Mod[j, 10] == \ 0, 
          Print[j, " ", i - 1, " ", k - 1]], \(j++\)}]\)

Probability distribution function compared to Boltzmann distribution

<< Statistics`DataManipulation`

tbl = Table[X[j, k], {j, 1, 1000}, {k, 1, 1000}];

t = BinCounts[Flatten[tbl], {0, 1, 0.001}]

midpoints = Table[i + 0.0005, {i, 0, 1 - 0.001, 0.001}]

su = Sum[t[[i]], {i, 1, Length[t]}]*0.001

plt1 = 
      ListPlot[Table[ {midpoints[[i]], t[[i]]/su}, {i, 1, 
            Length[t]}], PlotJoined ->  True,  
        PlotStyle  ->  RGBColor[1, 0, 0]]

kbT = 1;

plt2 = 
      Plot[Exp[-Ea Sin[Pi x]^2/kbT]/
          NIntegrate[Exp[-Ea Sin[Pi u]^2/kbT], {u, 0, 1}], {x, 0, 1}, 
        PlotStyle ->  RGBColor[0, 0, 1]]

Show[plt1, plt2]