SelfDiffusion

class SelfDiffusion(md_trajectory, start_time=None, end_time=None, atom_selection=None, anisotropy=None, time_resolution=None, info_panel=None)

Constructor for the SelfDiffusion object.

Parameters:

  • md_trajectory (MDTrajectory | AtomicConfiguration) – The MDTrajectory or configuration to calculate the self diffusion for.

  • start_time (PhysicalQuantity of type time) – The start time.
    Default: 0.0 * fs.

  • end_time (PhysicalQuantity of type time) – The end time.
    Default: The last time frame.

  • atom_selection (PeriodicTableElement | str | list of ints) – Only include contributions from this selection. The atoms can be selected by element i.e. PeriodicTableElement, tag or a list of atomic indices.
    Default: all atoms.

  • anisotropy (list of type int | int | None) – The list of Cartesian directions (x=0, y=1, z=2) to calculate the anisotropic self diffusion in or a single Cartesian direction. By default an isotropic calculation is performed.
    Default: None.

  • time_resolution (PhysicalQuantity of type time) – The time interval between snapshots in the MD trajectory that are included in the analysis.

  • info_panel (InfoPanel (Plot2D)) – Info panel to show the calculation progress.
    Default: No info panel.

data()

Return the self-diffusion values.

times()

Return the time values.

Usage Examples

Load an MDTrajectory, and calculate the self-diffusion of all aluminum atoms.

# Load the trajectory
md_trajectory = nlread('alumina_trajectory.hdf5')[-1]

# Calculate the self diffusion
self_diffusion = SelfDiffusion(
    md_trajectory,
    atom_selection=Aluminum
)

# Get the times in ps and the MSD values in Ang**2.
times     = self_diffusion.times()
diffusion = self_diffusion.data()

# Plot the diffusion coefficient calculated with different time limits
model = Plot.PlotModel(ps, cm**2/Second)
model.title().setText('Diffusion coefficient of aluminum')
model.xAxis().setLabel('Time')
model.yAxis().setLabel(r'$D_{self}$')
model.legend().setVisible(False)

# Add diffusion integral
line = Plot.Line(times, diffusion)
line.setColor('blue')

# Add average
average = Plot.Average()
line.addItem(average)
max_time = times[-1]
average.setBounds(max_time*0.1, max_time*0.9)
average.setColor('red')

# Add line to the plot
model.addItem(line)

# Set limits and show the plot
model.setLimits()
Plot.show(model)

self_diffusion.py

Notes

The SelfDiffusion is calculated from the velocity auto-correlation function using the Green-Kubo expression:

\[D_{self} = \frac{1}{3} \int_0^\infty d\tau \langle \mathbf{v}_i(t) \mathbf{v}_i(t+\tau) \rangle\]

Here the ensemble average \(<...>\) averages over both all selected atoms \(i\) and time origins \(t\).

The diffusion coefficient can also be calculated with the MeanSquareDisplacement analysis by calculating the gradient of the mean-squared displacement. These two analyses provide complementary methods of calculating the diffusion coefficient.

When calculating quantities using Green-Kubo expressions, it is important that the trajectory both contains a short enough time step to be able to determine the decay in the auto-correlation function, as well as a long enough time scale to encompass the auto-correlation time length and provide enough time origins for accurate sampling. When calculating the self diffusion, ideally as the velocity auto-correlation decreases to zero the above integral should converge to a static value at longer time scales. Practically however, longer times in the auto-correlation function corresponds to fewer time origins in the ensemble average. This increases the sampling error and can lead to non-convergence in the velocity auto-correlation functions and subsequently the self diffusion coefficient.

The SelfDiffusion analysis method returns the calculated self diffusion coefficient using different time limits. The optimal value of the self diffusion coefficient is usually where the auto-correlation has decayed sufficiently for the integral to converge, and also there are enough time origins for accurate sampling. The auto-correlation time scale can vary significantly depending on the potential acting on the atoms. For instance single atoms on which only non-bonding forces are acting typically have a longer auto-correlation time compared to atoms that have strong molecular bonds. In the molecular case relatively short trajectory time steps of the order of 10fs may be required.

By default, all elements are taken into account, but a specified selection can be given as well. The atom_selection parameter accepts an element, a tag name, or a list of atom indices to select which velocities are included in the calculation. This can be useful for when only certain atoms should be included in the analysis, such as when there are constrained atoms in the system that excluded. It is also possible to select to only calculate diffusion in specific cartesian directions using the anisotropy argument. This can be useful in cases such as calculating the diffusion rate of atoms between two surfaces.