Self-diffusion coefficient

mean_square_displacement(
    sim::Simulation,
    selection::AbstractVector{<:PDBTools.Atom};
    natomspermol::Integer,
    maxdelta::Integer = length(sim) ÷ 10,
    unwrap::Bool = true,
    parallel::Bool = true,
    show_progress::Bool = true,
)

Computes the mean square displacement (MSD), as a function of the time lag delta (in frames), of the centers of mass of the molecules defined by selection.

selection is expected to be the concatenation of n identical molecules, each with natomspermol atoms (the same convention used, for instance, by occupancy).

Trajectory files typically store coordinates wrapped back into the simulation box by periodic boundary conditions, which would introduce spurious jumps in a molecule's displacement whenever it crosses a periodic boundary between two frames. To avoid this, each molecule's center of mass is reconstructed into a continuous ("unwrapped") trajectory: at each frame, the periodic image closest to that same molecule's (already unwrapped) position at the previous frame is chosen. Since consecutive frames are assumed to be much closer in time than the time it takes a molecule to diffuse across half a box length, this reconstruction is unambiguous. If that assumption is violated, that is, if some molecule is reconstructed as moving by about half a unit cell vector between two consecutive frames, a warning is emitted (once per call), since the resulting MSD is then likely meaningless.

Some trajectories are already stored with continuous (unwrapped) coordinates. For those, set unwrap=false: the centers of mass are then used exactly as they are read, and neither the reconstruction nor the associated warning is performed.

Returns an OffsetArray with indices 0:maxdelta, in squared length units (typically Ų, matching the units of the input coordinates). The value at delta is the average, over all molecules and all pairs of frames separated by delta steps, of the squared displacement of the molecule's center of mass.

Use self_diffusion_coefficient to estimate the self-diffusion coefficient from the linear (diffusive) regime of the resulting curve.

Arguments

• sim::Simulation: The Simulation object.
• selection::AbstractVector{<:PDBTools.Atom}: The atoms of the n molecules considered, concatenated in sequence.

Optional keyword arguments

• natomspermol::Integer: Number of atoms of each molecule. Required.
• maxdelta::Integer: The maximum delta-step to be considered. Defaults to length(sim) ÷ 10.
• unwrap::Bool: Defines if the centers of mass are unwrapped before the displacements are computed. Defaults to true. Set it to false for trajectories whose coordinates are already continuous, so that they are used exactly as they are stored.
• parallel::Bool: Defines if the unwrapping and the averaging over time lags are run in parallel. Defaults to true. Requires starting Julia with multi-threading.
• show_progress::Bool: Show progress bar. Defaults to true.

Example

julia> using MolSimToolkit, PDBTools, MolSimToolkit.Testing

julia> sim = Simulation(Testing.namd2_pdb, Testing.namd2_traj);

julia> tmao = select(get_atoms(sim), "resname TMAO");

julia> msd = mean_square_displacement(sim, tmao; natomspermol=14, maxdelta=4, show_progress=false);

julia> msd[0]
0.0

julia> msd[4]
878.4761526314077

self_diffusion_coefficient(
    msd::AbstractVector;
    dt::Real = 1.0,
    dim::Integer = 3,
    mindelta::Integer = 1,
    maxdelta::Integer = lastindex(msd),
)

Estimates the self-diffusion coefficient from a mean-square-displacement curve msd, as computed by mean_square_displacement, using the Einstein relation MSD(t) = 2 * dim * D * t.

The coefficient is obtained as the slope of an ordinary least-squares fit of msd[delta] against t = delta * dt, for delta in mindelta:maxdelta, divided by 2 * dim.

Arguments

• msd::AbstractVector: The mean-square-displacement curve, as computed by mean_square_displacement.

Optional keyword arguments

• dt::Real: The time interval between consecutive frames of the trajectory used to compute msd, in whatever time unit is desired for the output (e.g. ps). Used to convert delta (in frames) into that time unit. Defaults to 1, that is, delta is used directly, and the resulting units of D are squared-length per frame.
• dim::Integer: The dimensionality of the diffusion process. Defaults to 3 (three-dimensional diffusion).
• mindelta::Integer: The smallest delta (in frames) included in the fit. Defaults to 1, excluding the trivial delta = 0 point.
• maxdelta::Integer: The largest delta (in frames) included in the fit. Defaults to lastindex(msd), that is, the whole curve.

Example

julia> using MolSimToolkit, OffsetArrays

julia> msd = OffsetArray([0.0, 2.0, 4.0, 6.0, 8.0], 0:4); # a purely diffusive MSD, dt=1, dim=3

julia> self_diffusion_coefficient(msd)
0.3333333333333334