Diffusion coefficient pro

Added in version 3.17.0.

This modifier computes the mean-squared displacement (MSD) of the particles as a function of lag time \(\tau\) from a simulation trajectory and derives the self-diffusion coefficient \(D\) via the Einstein relation

\({\displaystyle \mathrm{MSD}(\tau) = \left\langle \left| \vec{r}(t+\tau) - \vec{r}(t) \right|^2 \right\rangle = 2 n D \tau }\),

where \(n\) is the number of independent spatial directions (\(n = 3\) for the total MSD of a three-dimensional system, \(n = 2\) for a two-dimensional one). The linear fit is performed on the MSD itself, not on its square root, and includes a free intercept, \(\mathrm{MSD}(\tau) = 2 n D \tau + b\). The intercept \(b\) absorbs the ballistic/vibrational offset of the MSD at short lag times, which would otherwise bias the fitted slope; \(D\) is derived from the slope alone.

The modifier processes the entire input trajectory once, averaging the squared displacements over multiple time origins (windowed MSD). The result is cached, so stepping through the animation or changing the fit window does not retrigger the trajectory sweep. Changing any parameter that affects the displacement accumulation (time origins, selection, sampling), or any change to the upstream pipeline, restarts the computation.

How the windowed MSD is computed

In equilibrium, the MSD does not depend on the absolute time \(t\) at which a measurement starts, only on the lag time \(\tau\) separating the two configurations that are compared. Every pair of trajectory frames separated by \(\tau\) is therefore an equally valid sample of \(\mathrm{MSD}(\tau)\), and the windowed (also called multiple-origin or sliding-window) estimator averages over all of them:

\({\displaystyle \mathrm{MSD}(\tau) = \frac{1}{N_\mathrm{pairs}(\tau)} \sum_{t_0} \sum_{i} \left| \vec{r}_i(t_0+\tau) - \vec{r}_i(t_0) \right|^2 }\)

where the outer sum runs over all admissible time origins \(t_0\), the inner one over the particles contributing at that origin, and \(N_\mathrm{pairs}(\tau)\) counts all (particle, origin) pairs collected for that lag time. The average is thus taken over the pooled samples rather than over per-origin averages.

The sweep works in a single streaming pass over the trajectory, which never requires the full trajectory to be held in memory:

The frames selected by Sample every Nth frame are visited in order. Particle identity is established at the first visited frame and tracked by particle identifier from there on; particles that cannot be tracked throughout are dropped and counted in Diffusion.excluded_particles.
Every Origin stride-th sampled frame is registered as a new time origin. Registering means storing a snapshot of the unwrapped particle positions, the cell of that frame (used for the a/b/c decomposition), and, with Use only selected particles, the selection state.
At each sampled frame, the displacement of every tracked particle is computed against all time origins that are still live, and the squared displacement is added to the lag bin corresponding to the frame distance between that origin and the current frame.
An origin is retired, and its position snapshot released, as soon as the current frame is more than Maximum lag behind it. At any moment, the number of live origins is therefore at most Maximum lag / Origin stride (in sampled frames), which is what sets the memory footprint quoted under Parameters.

Two consequences of this scheme are worth keeping in mind when reading the resulting curve:

Statistics degrade towards long lag times: A lag of one sampled frame is sampled by every origin; a lag equal to the maximum lag is sampled by only a handful of origins near the beginning of the trajectory. The MSD curve therefore becomes progressively noisier towards the right end of the plot, which is why it can be advisable to limit Maximum lag (fraction) to e.g. half the trajectory and why the fit window should normally not be pushed into the last bins.
Samples are correlated, so the error is not \(1/\sqrt{N_\mathrm{orig}}\): Overlapping windows share trajectory segments and are not statistically independent. Reducing Origin stride to 1 maximizes the number of origins, but the accuracy gain saturates once the origin spacing exceeds the velocity correlation time of the system. A stride of a few frames usually costs little accuracy while cutting both memory and runtime proportionally.

Fixed reference mode is the degenerate case of this algorithm: a single origin is registered at Reference frame and every later frame contributes exactly one sample per particle to its own lag bin. Memory is then independent of the maximum lag, but each lag time rests on a single measurement.

Preparing the input trajectory

The modifier uses the input particle coordinates as-is. A few preparation steps, performed by standard modifiers inserted before this one in the pipeline, may be required:

Unwrapping

Displacements must be measured on continuous (unwrapped) particle trajectories. If your trajectory file stores coordinates wrapped back into the periodic simulation cell, insert the Unwrap trajectories modifier before this one.

Caution

Feeding wrapped coordinates directly into the analysis silently corrupts the MSD. Displacements are formed from the coordinates as they arrive, so with wrapped input no displacement can exceed the extent of the simulation cell: instead of growing linearly, the MSD saturates.

No unwrapping step is needed if the file already stores unwrapped coordinates (e.g. LAMMPS xu yu zu dump columns).

Constant simulation cell

The analysis requires a simulation cell whose size and shape do not change over the course of the trajectory. If the cell fluctuates (e.g. in an NPT simulation), the modifier aborts with an error message. Insert the Time averaging modifier, operating on the simulation cell, before this one to replace the fluctuating cell with its time average.

Caution

Usually diffusion coefficients are computed using a constant simulation cell size and shape. Please make sure you understand what you are doing before computing diffusion coefficients from a changing simulation cell.

Constant frame spacing

The physical time must advance by a constant amount from one sampled frame to the next, see Time axis below. A trajectory whose frames are spaced irregularly in time is rejected with an error message.

Dimensionality

Two-dimensional systems are fully supported, but the analysis takes the dimensionality from the simulation cell’s two-dimensional flag. The flag decides whether the out-of-plane contribution enters the total MSD and whether \(D\) is obtained as slope/\(2n\) with \(n = 2\) or \(n = 3\).

If the imported trajectory does not already declare the right dimensionality, set it with the Edit simulation cell modifier. As a safeguard, the modifier warns about a cell that is very flat but not marked two-dimensional. The dimensionality must not change over the course of the trajectory.

Center-of-mass drift correction (optional)

Thermostats and barostats can impart a net drift on the system, which contaminates the MSD at long lag times. Insert the Center of mass correction modifier (after the unwrapping step) to hold the center of mass in place at every frame.

Placing the fit window

The linear fit is restricted to a lag-time window, specified as fractions of the maximum lag time (Einstein fit interval group, Start/End). The window is displayed as a shaded band in the MSD plot of the modifier panel, and its edges can be dragged with the mouse to adjust the fit range interactively. Since the accumulated MSD data is cached, changing the fit range only redoes the inexpensive linear fit.

To help place the window, the plot can be switched to doubly logarithmic axes with the Log-log axes option below it. A power law \(\mathrm{MSD}(\tau) \propto \tau^{\alpha}\) appears as a straight line of slope \(\alpha\) there, which makes the dynamic regimes easy to tell apart: slope 2 at short lag times marks the ballistic regime, slope 1 the diffusive regime the Einstein fit requires, and a slope below 1 at long lag times indicates subdiffusive behavior (caging or confinement) or simply poor statistics. Details:

E.J. Maginn, R.A. Messerly, D.J. Carlson, D.R. Roe and J.R. Elliott,

Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics,

Living J. Comput. Mol. Sci. 1(1), 6324 (2019)

https://doi.org/10.33011/livecoms.1.1.6324

The Log-log axes option only changes the display of the plot.

Parameters

Time origins

Selects how the time origins of the displacement measurement are chosen.

In Windowed mode (the default), the squared displacements are averaged over all available time origins. This is the statistically preferable choice, because each lag time is sampled many times.

In Fixed reference mode, all displacements are measured relative to the single frame given by Reference frame. This is cheaper and closer to a textbook \(\left|\vec{r}(t) - \vec{r}(0)\right|^2\), but the statistics are much poorer.

Origin stride

Number of trajectory frames between successive time origins in windowed mode.

Note

The calculation keeps a snapshot of all particle positions in memory for every time origin within the maximum lag. Both the memory footprint and the computation time therefore scale as (number of particles) × (number of frames) × Maximum lag / Origin stride. For long trajectories of large systems this can amount to several gigabytes; the modifier issues a warning when it does. Increasing this value is the cheapest way to bring the cost down, at the price of averaging over fewer time origins.

Maximum lag (fraction)

The largest lag time for which the MSD is accumulated, as a fraction of the trajectory length. Lag times close to the trajectory length can only be sampled by very few time origins and are correspondingly noisy, so consider lowering this value, e.g. to 0.5. Lowering it also reduces the memory footprint and the computation time, see the note above.

Reference frame

The trajectory frame serving as the time origin in Fixed reference mode.

Sample every Nth frame

Processes only every Nth frame of the input trajectory. This coarsens the lag-time axis but speeds up the analysis proportionally.

Distance / Time

Declare the physical units of the input data, see below.

Source

Selects where the physical lag-time axis comes from, see below.

Timestep

The physical time interval between two successive trajectory frames, used when Source is Fixed timestep. Note that this is the spacing of the frames in the trajectory file, not the integration timestep of the MD simulation.

Attribute name

The global attribute providing the physical time of each frame, interpreted in the selected Time unit. Used when Source is Global attribute.

Use only selected particles

Restricts the MSD average to the currently selected particles. The selection state is evaluated at each time origin. Use this to obtain the diffusion coefficient of one species within a mixture.

Partition by particle type

Additionally computes a separate MSD curve and diffusion coefficient for each particle type, analogous to the partial radial distribution functions.

Caution

A particle is assigned to its type’s set once, at the first processed trajectory frame, and stays in that set for the whole analysis even if its type changes later on.

Partition by direction

Additionally resolves the full displacement/diffusion tensor, see below.

Einstein fit interval: Start / End

The lag-time window used by the linear fit, as fractions of the maximum lag time. The default range skips the short-time ballistic regime. See Placing the fit window above.

Time axis

The physical lag-time axis can be derived from a fixed timestep per frame, from a global attribute of the input trajectory (e.g. Time or Timestep), or simply measured in animation frames. In the first two cases the values are interpreted in the selected Time unit.

The modifier requires a linear time axis: the physical time must advance by the same, positive amount from one sampled frame to the next. Every lag bin then corresponds to exactly one physical lag time, namely an integer multiple of that step. The Fixed timestep and Animation frames sources satisfy this by construction; with Global attribute the spacing is taken from the first two sampled frames and verified at every subsequent frame.

Caution

If the time spacing is not constant, or if the time does not increase from one sampled frame to the next, the modifier aborts with an error message naming the two offending frames. Typical causes are a global attribute that is not a time axis, a trajectory concatenated from several restarted runs, and frames written at irregular intervals. Such a trajectory has to be resampled, split, or analysed with a different time source; the modifier does not attempt to compensate for a non-uniform frame spacing.

Units

The Distance and Time settings declare the physical units of the input data: distances in Å, nm or m, times in fs, ps or s. Both settings are purely declarative and no input values are rescaled. They state which units the input coordinates and the time axis are already expressed in, so that the axis labels of the MSD plot and the reported diffusion coefficient carry the right unit, and so that Diffusion.D_cm2_per_s can be computed from the correct conversion factor.

Data in reduced or otherwise unknown units can be marked Dimensionless. The results are then reported without a unit and no Diffusion.D_cm2_per_s attribute is generated. The same applies when the time axis is measured in animation frames, in which case the Time unit does not apply and is disabled.

Outputs

The modifier emits its results as a data table and a set of global attributes:

MSD data table

A data table named msd holds the MSD as a function of lag time. By default, it contains only the total MSD. When Partition by direction is enabled, additional columns hold the components of the symmetric displacement tensor \(\langle \Delta r_\alpha \Delta r_\beta \rangle\) in the order XX, YY, ZZ, XY, XZ, YZ, and – only if the simulation cell is not axis-aligned – the diagonal components a, b, c along the cell vectors (for an axis-aligned cell they would merely duplicate XX, YY, ZZ and are therefore omitted). When Partition by particle type is enabled, additional per-type columns are generated. Each data column is immediately followed by its fitted-line column (named e.g. XX (fit)), holding the Einstein fit evaluated over the entire lag range, so that the fit can be shown alongside the data in the same plot. Without any partitioning, the table thus consists of just two columns, named Total and Fit. The off-diagonal components XY, XZ, YZ quantify correlated motion between two axes; unlike the diagonal ones, they can be negative (anticorrelated motion) and vanish for isotropic diffusion.

The cell-vector components are obtained from the reduced (fractional) coordinates of the displacement vectors, i.e., from the decomposition \(\Delta\vec{u} = s_a \vec{a} + s_b \vec{b} + s_c \vec{c}\) evaluated with the origin frame’s cell. This yields independent components even for triclinic cells, unlike a naive projection onto the normalized cell vectors. Note that for non-orthogonal cells the a/b/c components do not sum to the total MSD (cross terms), whereas the XX, YY, ZZ components always do.

Two-dimensional systems (see the simulation cell’s dimensionality setting) are fully supported: the out-of-plane columns ZZ, XZ, YZ and c and the corresponding attribute components are omitted, and the total MSD contains only the in-plane contributions.

Diffusion coefficient

The attribute Diffusion.D holds the self-diffusion coefficient obtained from the linear fit, in units of Distance²/Time (per animation frame if the time axis is measured in frames). Diffusion.D_cm2_per_s provides the same value converted to cm²/s. The latter is only emitted if both units denote a physical unit and the time axis is not measured in animation frames.

Fit parameters

Diffusion.fit_slope and Diffusion.fit_intercept are the two parameters of the linear fit \(\mathrm{MSD}(\tau) = \mathrm{slope} \cdot \tau + b\), the intercept \(b\) being the MSD extrapolated to zero lag time. Diffusion.fit_r_squared is the coefficient of determination of the fit, and Diffusion.fit_lag_min and Diffusion.fit_lag_max report the lag-time window actually used by it.

Excluded particles

Diffusion.excluded_particles counts the particles that were excluded from the analysis because they could not be tracked across all trajectory frames. This attribute is only emitted if there are any.

Naming of partitioned results

With Partition by particle type enabled, every attribute listed above additionally comes in per-type variants such as Diffusion.D.Cu, following the naming convention of the partial radial distribution functions.

With Partition by direction enabled, the spatially resolved results are not emitted as one attribute per component. Instead, the components of a family are packed into a single multi-component attribute: Diffusion.D.tensor holds the six components of the symmetric diffusion tensor in the order (xx, yy, zz, xy, xz, yz), obtained as half the fitted slope of the corresponding displacement-tensor component, and Diffusion.D.abc holds the diagonal cell-vector components (only emitted for non-axis-aligned cells). Diffusion.D_cm2_per_s, Diffusion.fit_slope, Diffusion.fit_intercept and Diffusion.fit_r_squared come in the same .tensor/.abc variants. In two-dimensional systems the tensor attribute has the three in-plane components (xx, yy, xy) and the cell-vector attribute becomes the two-component Diffusion.D.ab. A packed attribute is omitted if the fit failed for any of its components.

In expression inputs and other places where attributes are addressed by component, the tensor components are available as Diffusion.D.tensor.XX through Diffusion.D.tensor.YZ (in 2D, the three components of the in-plane value are addressed generically as .X, .Y, .Z, meaning xx, yy, xy).