MeanSquareDisplacement

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

Constructor for the MeanSquareDisplacement object.

Parameters:

  • md_trajectory (MDTrajectory | AtomicConfiguration) – The MDTrajectory or configuration to calculate the mean-square-displacement 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 mean square displacement in or a single Cartesian direction. By default an isotropic calculation is performed.
    Default: None.

  • use_com (bool) – Whether or not displacement is calculated from molecular center of mass.
    Default: False.

  • 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 mean-square-displacement values.

times()

Return the time values.

Usage Examples

Load an MDTrajectory, and calculate the mean-square-displacement (MSD) of all aluminum atoms. Estimate the diffusion coefficient from the slope of the MSD-curve, according to \(MSD(t)=6 Dt\):

md_trajectory = nlread('alumina_trajectory.nc')[-1]

msd = MeanSquareDisplacement(md_trajectory, atom_selection=Aluminum)

# Get the times in ps and the MSD values in Ang**2.
t        = msd.times().inUnitsOf(ps)
msd_data = msd.data().inUnitsOf(Angstrom**2)

# Plot the data using pylab.
import pylab

pylab.plot(t, msd_data, label='MSD of aluminum')
pylab.xlabel('t (ps)')
pylab.ylabel('MSD(t) (Ang**2)')
pylab.legend()

pylab.show()

# Fit the slope of the MSD to estimate the diffusion coefficient.
# If you discover non-linear behavior at small times, discard this initial part in the fit.
a = numpy.polyfit(t[5:], msd_data[5:], deg=1)

# Calculate the diffusion coefficient in Ang**2/ps.
diffusion_coefficient = a[0]/6.0

mean_square_displacement.py

Notes

The MeanSquareDisplacement is calculated as:

\[MSD(t) = \left \langle \left [ \mathbf{r}(t) - \mathbf{r}(0)\right ]^2 \right\rangle.\]

In practice, the average \(<...>\) runs over all selected atoms \(i\) in the trajectory, and an additional average over simulation time is carried out to improve the statistical sampling. That means for a given time difference \(t\) all image pairs that are separated by \(t\) are taken into account in the average, as

\[MSD(t) = \frac{1}{(t_{max} - t)N_{atoms}} \sum_{i=1}^{N_{atoms}} \sum_{t'=0}^{t_{max}-t} \left [ \mathbf{r}_i(t' + t) - \mathbf{r}_i(t')\right ]^2 .\]

Note, that this requires a system which is equilibrated, i.e. its macroscopic properties do not change during the simulation.

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 indices to select atoms for the velocity distribution. This can be useful, e.g. in the presence of constraints as constrained atoms should be excluded in this analysis.