issues with calculating diffusion coefficient when jumps

Hi all.

I have previously done diffusion of hydrogen in iron tracking the MSD with no issues. I am now moving on to diffusion in SiO2, where the diffusion occurs via jumps. In this case it is not as simple as just the slope of MSD over the entire time. Instead, several small intervals are averaged.

While I am getting a value, it seems significantly lower than other values in literature, making me wondering if I am missing something. My procedure is:

  • create amorphous SiO2 (check density), add hydrogen, equilibrate the structure in an NPT at desired temperature and 1 atm, then compute the MSD in an NVT for 4ps (keeping the temperature at the desired value)
  • I do see several jumps. however, what I see is a large jump in MSD right when it starts, making me think that I should be adding something to minimize my hydrogen before running the NPT and tracking the MSD. Perhaps an NVE limit?
  • for diffusion of hydrogen in iron I didn’t need this, and did not observe this big jump, but maybe it is different for this system

Is there something in my methodology that seems wrong? Any advice would be great.


When you say “lower than other values in the literature” are you referring to other MD simulations? If yes, are you using the same interatomic potential, the same protocols for buidling the initial structures, the same simulation duration, etc.? If the answer to any of these questions is no, then there is no reason to expect agreement between your estimate and theirs.

Here are a few more comments:

  1. 4 ps is a very shrot simulation. See what happens when your run for 10x longer, 100x longer

  2. The large jump in MSD indicates that something is badly wrong. You said that you had already equilibrated the system at the desired temperature, but apparently that is not true.

  3. Things like NVE limit and even barostats and thermostats tend to mask underlying problems. You should only use them to bring the system in to equilibrium, and then run NVE, or maybe NVT dynamics.