How are diffusion coefficients calculated?

Question

To support my laboratory result, I'm searching for an MD/statistical mechanics approach to predict translational diffusion coefficients in the liquid phase.

I'm aiming to describe small peptides and small organic molecules and in future proteins and protein aggregates.

Answer 1

There are several different methods to calculate diffusion coefficients depending on the problem and some nuances to consider.

Mean square displacement

The most common way to calculate diffusion coefficients in homogeneous media is to calculate the mean square displacement. This is what the Gromacs tool described in the other answer does. However, these calculations are relatively easy to do on your own (using a script in Python, MathLab, or whatever language you like or even a spreadsheet). The general formula is:

$$D = \frac{1}{2} \lim_{t\rightarrow \infty} \left< \left[ x(t)-x(0) \right]^2 \right>/t$$

Here, the mean $\left< \cdots \right>$ is over many different trajectories (probably at least 20). These trajectories can be from multiple simulations with different initial conditions, multiple molecules in the same simulation, or even different time intervals of a single molecule in a single simulation.

You can shift your time origin to $t'$ to make use of different intervals. If your system is three-dimensional and isotropic (all directions are the same), you can use the vector displacement along all three axes to get:

$$D = \frac{1}{6} \lim_{t\rightarrow \infty} \left< \left| \mathbf{r}(t - t')-\mathbf{r}(t') \right|^2 \right>/t$$

The plot below from Mark et al. shows some mean squared displacements for the TIP3P water molecule as a function of time, $\mathsf{MSD}(t) = \left< \left| \mathbf{r}(t - t')-\mathbf{r}(t') \right|^2 \right>$. Here the mean is calculated over all 901 water molecules in the system. The different curves represent different time origins. From here, you can calculate the slope of the curve using a least squares fit (I'm not sure this is the best way in terms of error analysis, but this is often how its done). In the graph below, the mean square displacement is roughly 300 Ų in 100 ps, so $D \approx \frac{1}{6} 300~Å^2/(100~\mathrm{ps}) = 0.5~Å^2/\mathrm{ps}$.

Given that you have a limited simulation time, you usually have to choose between having more subtrajectories to get better statistics or having a longer time ($t$) for each subtrajectory. More subtrajectories give you better statistics, but if your trajectories are too short, they will not yet be in the diffusive regime ($\left<x^2 \right> \sim t$). As shown in the MSD graph below from Flenner et al. for lateral diffusion of lipids in bilayer membrane, different regimes can exist at different timescales. At short times motion is ballistic (the particles move a near constant velocity and $\left<x^2 \right> \sim t^2$). At longer times, there is an intermediate subdiffusive regime ($\left<x^2 \right> \sim t^{0.7}$). Standard diffusion ($\left<x^2 \right> \sim t$) only appears after many nanoseconds.

Finite-size effects

The diffusion coefficient that you measure will be biased somewhat due to the small size of the simulation system, which leads to hydrodynamic interactions between particles and their periodic images. How small is too small? For the periodic boundary conditions typically used in molecular dynamics simulations, you can use the correction of Yeh and Hummer:

$$D_\mathrm{corrected} = D_\mathrm{PBC} + 2.84 k_\mathrm{B}T/(6 \pi \eta L),$$

where $D_\mathrm{PBC}$ is the diffusion coefficient calculated in the simuilation, $k_\mathrm{B}$ is Boltzmann's constant, $T$ is the temperature, $\eta$ is the shear viscosity of the solvent, and $L$ is the dimension of the cubic box. With modern computers, you can usually make your systems big enough to where this correction is negligible.

Velocity autocorrelation

While the MSD method is the most common for homogeneous systems, there are other methods, such as the velocity autocorrelation method:

$$D = \int_0^\infty \left< v(0) v(t) \right> \, \mathrm{d}t.$$

For isotropic system in 3D, you can write,

$$D = \frac{1}{3} \int_0^\infty \left< \mathbf{v}(0) \cdot \mathbf{v}(t) \right> \, \mathrm{d}t.$$

Inhomogeneous systems

Things get much more difficult if your system is not homogeneous in the direction that you need to calculate the diffusion coefficient. In this case, the diffusion coefficient likely varies with position, $D(x)$. Furthermore, in inhomogeneous systems the particles typically experience systematic forces in addition to the random forces of the medium. The effect of these position-dependent forces, $F(x)$, also must be considered.

To calculate the diffusion coefficients along the coordinate, you can restrain your particle harmonically to different positions. For the case of a particle in an approximately harmonic well, an unbiased diffusion coefficient can be calculated based on the Generalized Langevin Equation:

$$D(z) = \frac{\left<\delta z^2 \right>^2}{\int_0^\infty \left< \delta z(t) \delta z(0) \right> \, \mathrm{d}t},$$ where $\delta z(t) = z(t) - \left<z\right>$ is the displacement from the mean position in the well.

This method is commonly used for inhomogeneous systems in the context of umbrella sampling calculations (which also obtain the free energy and its gradient, the mean force $F(x) = -\nabla A(x)$). however, more robust methods based on Bayesian inference have been proposed and can obtain $D(x)$ in the case of arbitrary position-dependent forces. See Hummer 2005 and others.

Answer 2

I decided to make my comment an answer...

Diffusion coefficients are simple to calculate in molecular dynamics. GROMACS has built in software specifically for this calculation.

GROMACS diffusion coefficient manual

You simply run your "normal" NVT or NPT simulation, and then call this function

gmx msd -f *.trr

adding additional options as you see fit, and replacing the wildcard * with the actual file name if more than one .trr file exists in the folder. GROMACS will ask you what molecules or groups you would like to know the diffusion coefficient for, and there you go.

This will probably take less than 1 hour of calculation time, perhaps a couple if you really want to equilibrate and sample well (note: time depends on cpu or gpu ratio to number of atoms - longer times may occur if there are few computational resources thrown at the simulation). The amount of time taken by gmx msd is in the seconds to minutes... it is quite fast.