Two ways to do it, not any of them in LAMMPS. 

First, you can transfer your trajectory to VMD, which has a diffusion coefficient calculator tool.

Or you can do what it does, manually, using a simple code. What you need to do is just average the msd in sample bins, like for 1ns, 2ns, 5ns, 10ns etc. If you have n long run, the largest bin will be n/2 for sure. It will give you a smooth MSD behaviour, from which you can calculate diffusivity.

Thank you for your answer.

Using VMD is very interesting. I don't know such function of VMD.

I will check it.

Could you tell me more detail about latter one?

I know how to output the position data in LAMMPS calculation, it is just like

dump pos Ow custom 10 pos.dump id type xu yu zu

And I can understand your answer. If MSD is sampled during 100ns by 1ns interval, MSD in the first 1ns is the average from 100 samples. And next MSD, time = 2ns is the average from 99 samples. Thus, calculated MSD until the half of sampling time (50ns, in this case) is good.

Unfortunately, I'm not familiar with programing. Could you tell me how to extract each position data from dump file or sample code?

You have nothing to do with positions, msd file is good enough. 

All you need is a simple code that substracts (i+n)-i for whole array (which is from 1:(i+n)) and average the values.

Thank you so much. I just got it!

I will try it.

Thanks again!

I tried to calculate the MSD using msd file from LAMMPS but I think something wrong my calculation. Could you give me some advices?

MSD calcualtion in LAMMPS was performed by;

compute msdOw Ow msd com yes

fix msdOw Ow ave/time 10 1 10 c_msdOw[4] file wat.msd

As a result of LAMMPS calculation, MSD was calculated each 10 time steps.

When obtained MSD data from LAMMPS are defined as 0, a1, a2, a3, a4,,,,an, then I performed MSD calculation as ((a1-0)+(a2-a1)...)/n.

But this calculation can be written as an/n. Other terms between a1 and an-1 are offset each other.

I think my calculation is wrong.I'm very sorry I didn't understand your comments perfectly.

Could you give me more comments?

Thank you,

You are right for the smallest bin, which is the case you shown. So this is your average value with 10 timesteps bin.

Then you have, ((a2-0)+(a3-a1)+(a4-a2)...) 

goes on like this. 

I found my matlab code for a 40 ns run with 20 bins which may be helpful for you

n=length(A)

for k=1:20

k1=100*k;

b=0;

for i=1:1:(n-k1);

k2=k1+i;

b(i)=A(k2,2)-A(i,2);

end

%b(b

Thank you for your comment.

Unfortunately, I'm not matlab user. So I have little confidence about second loop;

for i = 1:1:(n-k1)

.

.

end

%b(b

That is commented out (%).

average of b is taken, then it is written as c(k,2).

Well... I see.

Thank you for your comments.

I'm interested in calculating the diffusion coefficient of self interstitial atom in iron system at 800 K. I got a plot like (attached) this 'MSD vs time' using LAMMPS. The slope (diffusion coefficient) is changing with time. which part of the curve to be considered to calculate the slope (diffusion coefficient) ?

Kedharnath Arumugakani

Dear,

Can I ask you how did you obtain your results in LAMMPS.

My MSD figure has a lot of fluctuations and goes up and down.

Dear Ensieh Yousefi,

I got the plots as attached using the LAMMPS commands:

#Heating

velocity all create 800.0 4928459 dist gaussian sum yes units box

fix 1 all npt temp 800.0 800.0 1 iso 0.0 0.0 10.0 dilate all drag 1

#MSD calculation

fix 1 all nve

compute 3 all msd com yes

variable msd equal c_3[4]

I'm not sure whether this is the correct script but the plots seems to be ok.

Dear Ensieh Yousefi,

I think these might solve the problem (but still I'm not sure)- have less number of atoms which will reduce the fluctuations due to other atoms in the systems during displacements.

How to calculate diffusion coefficient in VMD?