The following work can be helpful regarding the diffusion coefficient:

https://www.researchgate.net/publication/236586471_Silver-Sodium_Ion_Exchange_Dynamics_in_LTA_Zeolite_Membranes?_iepl%5BviewId%5D=ZyGeBQIoCaqM12lKur81n1HV&_iepl%5BprofilePublicationItemVariant%5D=default&_iepl%5Bcontexts%5D%5B0%5D=prfpi&_iepl%5BtargetEntityId%5D=PB%3A236586471&_iepl%5BinteractionType%5D=publicationTitle

Article Silver−Sodium Ion Exchange Dynamics in LTA Zeolite Membranes

I'm not calculating diffusion coefficient using Einstein relationship. I'm using MSD values instead. In this case, please let me know your views. I have attached MSD vs time plot.

You can calculate the diffusion coefficient (slope), as the plot divided into two domains, including stating and final MD time, so you can report the diffusion coefficient in two different MD times. Seems that after 100000, the diffusion is almost near zero.

In the beginning of simulation all atoms move out of their crystalline position - this is so called ballistic regime. It takes some time to establish thermal equilibrium (in this case, however, it takes an unusually large number of timesteps - my guess is that the timestep is very small). Therefore, the slope in the ballistic regime (initial stage) has nothing to do with diffusion. Eventually MSD becomes linear with time. In your case, as correctly noticed above - MSD becomes constant with zero slope and that means zero diffusion. it might be that you calculate MSD for all atoms while you are interested in other atoms. Create 2 kinds of atoms - even though interactions might be identical and follow MSD of those atoms you are interested in. Also, increase number of timesteps to see if MSD is indeed zero. Without knowing your problem in more details it is tough to describe all possible traps.

Dear Anatoly Belonoshko,

I equilibrated the system for 10 ps and then ran the MSD calculation for 20 ps with timestep being 0.0001 ps. After 10 ps of MSD calculation, there is still slope giving diffusivity values in the order of E-4 angstrom^2/ps.

As you suggested, I created 2 atom types, one for lattice atoms and another for interstitial atom. The calculated MSD value of an interstitial atom is very small and fluctuating a lot. So, I think its better to consider the MSD of the whole system.

My problem is to calculate diffusion coefficient of a self interstitial atom in Iron system at different temperature. Doing so, I can get Activation or Migration energy of that interstitial atom.

Hi,

   I believe if all of your atoms are iron and no distances are extremely short, you can safely use  the 10 times larger timestep, that is 0.001. Run your system to equilibrate for

50 ps and then 200 ps to get MSD. I strongly suspect you get zero diffusion. Eventually you get some deformed clusters of atoms where no atoms shall diffuse. This is at 800 K. If you go close to melting T, everything is possible but definitely not at a low T.

Dear Anatoly Belonoshko,

Yes, all atoms belong to Iron including the interstitial atom.

As  you suggested, I will equilibrate for 50 ps and run simulation for 200 ps calculating MSD values with timestep 0.001 at higher temperature (such as 1400 K) and get back to you.

Thanks a lot for your guidance.

Dear Anatoly Belonoshko,

As you suggested I have plotted MSD vs run time for both timestep 0.0001 and 0.001 ps. There is still a little value of slope.

In the plot: time = run time * timestep

Please let me know your feedback on this.

We are assuming that the longer you run MD the more equilibrated results are. In that spirit, if you take the last 50 ps of the run with 0.001 ps timestep the diffusion is actually zero. It is of course possible that atoms jumps from time to time but those are extremely rare events. To get statistics ove rare events one has to have extremely long runs, perhasp even 200 ps are not sufficient. Anyway, I see that the slope from 100 to 150 ps is different from the slope from 150 to 200 ps. That means simply that your results are not converged yet. Calculate please the slopes between 100 to 150 and from 150 to 200 ps. 

Dear Anatoly Belonoshko,

When I visualized the interstitial movement, I observed that, initially there was movement of interstitial atom from site to lattice. Later, there is no major change in its position. But you are asking me to calculate the slope (diffusion coefficient) at later or after longer runs where is there no much change in the position of the interstitial.

https://youtu.be/qkQkRaWEoAc

Hi,

   I believe you misunderstood me. I am not asking you anything. You asked an advice how to calculate diffusion - and I have given you such. Whether use it or not - is completely up to you. If you are interested in just relaxation of your initial configuration - which is what happens when interstitial atoms leave their places - you are done. If you are interested in the diffusion - you have to run the simulation until the diffusion becomes constant.

Dear Anatoly Belonoshko,

Sorry, we misunderstood each other. I mean to say that, you advised me to calculated the slope after longer run but what we observe from the simulation is the initial diffusion of the interstitial atom. I'm sorry that I communicated this wrongly.

Thanks for your kind advice and discussion. I understood that in regards to calculation of diffusion coefficient, the simulation should be run for longer time.

I did not understand what exactly you need, in following work the diffusions were determined based on MSD profiles...(Table 6)

https://www.researchgate.net/publication/236586471_Silver-Sodium_Ion_Exchange_Dynamics_in_LTA_Zeolite_Membranes?_iepl%5BviewId%5D=ZyGeBQIoCaqM12lKur81n1HV&_iepl%5BprofilePublicationItemVariant%5D=default&_iepl%5Bcontexts%5D%5B0%5D=prfpi&_iepl%5BtargetEntityId%5D=PB%3A236586471&_iepl%5BinteractionType%5D=publicationTitle 

Article Silver−Sodium Ion Exchange Dynamics in LTA Zeolite Membranes

Dear Ramin Ekhteiari Salmas,

I was able to plot out Mean square displacement vs time from LAMMPS MD simulations. The plot follows the above trend (image in the question). I'm unable to understand the profile and thus had confusion is selection of slope from that plot.

It's clear that the diffusion is high at the initial times of MD as expected , and then it came into an equilibrated state (diffusion of those atoms is zero). The MSD profiles show  you the dynamics of the atoms over the simulations, as the diffusion constant is a dynamic factor too....and you may analysis the diffusion degrees along the simulations...

Wanting to calculate the diffusion constant of those atoms into the...you may consider the slope between 0 to 60000.

Dear all,

I am trying to calculate MSD using he Diffusion Coefficient Tool for VMD (https://github.com/tonigi/vmd_diffusion_coefficient). However, the result plot of MSD is not linear (attached figure). I don't know why and this figure can not be used for linear fitting the diffusion coefficient with Einstein -Relation. My simulation was 1 M NaCl electrolyte solutions just in bulk state based on NAMD. Could anyone help me?

Hi

Looks like a zero diffusion case to me. In the beginning all atoms start moving out of their initial postions (so called ballistic regime) and that explains a parabolic growth of MSD up to 1 (picosecond? I cant see units on your graph). After that pretty much no increase of MSD - zero slope = zero diffudion coefficent. I am confused by the word 'solution' - in solution one would expect pretty considerable diffusion. It could be that the temperature is below the melting/solidification T of your model or concentration of electrolyte leads to freezing well above the melting T of the model.

The diffusion in crystalline systems is based on the atoms jumps between the lattice, vacancy and interstitial sites, which are rare event. Therefore, you need a long production simulation and only measuring the MSD of the atoms that leave their crystalline sites. To increase the chance to have jumps you have to use a large system and high temperature for the simulations.

Usually the MD simulations are done using an ideal monocrystal as starting structure. In reality the materials are policrystallines, with a high concentration of point and extended defects.

Thanks for your kind answers Dr. Viorel Chihaia. In metallic system with grain boundaries, dislocations and precipitates, the diffusivity of the interstitials is different in each cases. In an attempt to carry out diffusivity calculation at high temperature, the structural integrity of the grain boundaries, dislocations and precipitates seems lost which makes inaccurate diffusivity calculation. I'm trying to look into NEB or Dimer method to calculate migration energy. If you have any idea about the above techniques and codes relating it to calculation of migration energy of interstitial, I kindly request you to share with the research community.

Here is a list with some articles (by doi address) that deal with the electronic structure calculations for the T-dependent diffusion within crystals

10.1103/PhysRevB.77.134305 (2008- year)

10.1016/j.actamat.2009.05.006 (2009)

10.1016/j.actamat.2009.11.041 (2010)

10.1016/j.commatsci.2010.08.019 (2010)

10.1007/s11669-011-9854-5 (2011)

10.1016/j.actamat.2011.01.062 (2011)

10.1007/s11669-012-0185-y (2013)

10.1038/sdata.2016.54 (2016)

For an automated calculations look for MAST (https://pythonhosted.org/MAST/, https://github.com/uw-cmg/MAST or http://matmodel.engr.wisc.edu/software/) that works with VASP within pymatgen.

Thank you so much for your kind help Dr. Viorel Chihaia. The research articles and software were very helpful. I will look into MAST for my scientific problems.

Hi Kedharnath Arumugakani you can find it in "Helmut Mehrer Diffusion in Solids

Fundamentals, Methods, Materials,Diffusion-Controlled Processes"

according to this book, diffusion coefficient can be obtained from the slop of MSD curve.

Hi All,

How can one add atoms of carbon at interstitial points in bcc iron using LAMMPS or any other molecular dynamics software?