Tukukkk

I have a particle trajectory in 1D, j= and for the time=np.arange(0, 10 + dt, dt) where dt is the time step. I have calculate the MSD according to this article.

I have searched in google as well as here for 1d MSD in python but did not find any suitable one as my python knowledge is very beginner level. I have written one code and it is working without any error but I am not sure that it represent the same thing according to the given article. Here is my code,

j_i = np.array(j)

MSD=

diff_i=

tau_i=

for l in range(0,len(time)):

    tau=l*dt

    tau_i.append(tau)

    for i in range(0,(len(time)-l)):

        diff=(j_i[l+i]-j_i[i])**2



        diff_i.append(diff)



    MSD_j=np.sum(diff_i)/np.max(time)

    MSD.append(MSD_j) 

Can anyone please check verify the code and give suggestion if it is wrong.

The code is mostly correct, here is a modified version where:

• I simplified some expressions (e.g. range)


• I corrected the average, directly using np.mean because the MSD is a squared displacement [L^2], not a ratio [L^2] / [T].

Final code:

j_i    = np.array(j)

MSD    = 

diff_i = 

tau_i  = 



for l in range(len(time)):

    tau = l*dt

    tau_i.append(tau)



    for i in range(len(time)-l):

        diff = (j_i[l+i]-j_i[i])**2

        diff_i.append(diff)



    MSD_j = np.mean(diff_i)

    MSD.append(MSD_j)

EDIT: I realized I forgot to mention it because I was focusing on the code, but the ensemble average denoted by <.> in the paper should, as the name implies, be performed over several particles, preferentially comparing the initial position of each particle with its new position after a time tau, and not as you did with a kind of time-running average

EDIT 2: here is a code that shows how to do a proper ensemble average to implement exactly the formula in the article

js     = # an array of shape (N, M), with N the number of particles and

         # M the number of time points

MSD_i  = np.zeros((N, M))

taus   = 



for l in range(len(time)):

    taus.append(l*dt)  # store the values of tau



    # compute all squared displacements at current tau

    MSD_i[:, l] = np.square(js[:, 0] - js[:, l])



# then, compute the ensemble average for each tau (over the N particles)

MSD = np.mean(MSD_i, axis=0)

And now you can plot MSD versus taus and Bob's your uncle