Skip to content

Dynamics

Dynamic analysis routines for trajectory data.

compute_D(time, msd, fit_from=0, dim=3)

Compute the diffusion coefficient from an MSD curve.

Parameters:

Name Type Description Default
time ndarray

Time values corresponding to the MSD entries.

 required
msd ndarray

Mean-squared displacement values.

 required
fit_from int

Index at which to start the linear regression. Default is 0.

 0
dim int

Dimensionality of the system (e.g., 3 for 3D). Default is 3.

 3

Returns:

Type Description
float

Estimated diffusion coefficient in units of msd/time.
Source code in m2dtools/basic/dynamics.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def compute_D(time, msd, fit_from=0, dim=3):
    """Compute the diffusion coefficient from an MSD curve.

    Parameters
    ----------
    time : np.ndarray
        Time values corresponding to the MSD entries.
    msd : np.ndarray
        Mean-squared displacement values.
    fit_from : int, optional
        Index at which to start the linear regression. Default is ``0``.
    dim : int, optional
        Dimensionality of the system (e.g., ``3`` for 3D). Default is ``3``.

    Returns
    -------
    float
        Estimated diffusion coefficient in units of ``msd/time``.
    """
    coeffs = np.polyfit(time[fit_from:], msd[fit_from:], 1)
    D = coeffs[0] / 2 / dim # in unit of msd / time
    return D

compute_SISF(positions, box_lengths, lag_array, k, num_vectors=100, n_repeat=100)

Compute the self-intermediate scattering function (SISF).

Parameters:

Name Type Description Default
positions ndarray

Atom positions with shape (n_frames, n_atoms, 3).

 required
box_lengths ndarray or list

Simulation box lengths per frame with shape (n_frames, 3).

 required
lag_array ndarray or list

Lag times (in frames) at which to evaluate the SISF.

 required
k float

Magnitude of the wave vector |k|.

 required
num_vectors int

Number of random k-vectors to sample. Default is 100.

 100
n_repeat int

Number of time origins to average over. Default is 100.

 100

Returns:

Type Description
ndarray

Self-intermediate scattering values for each lag time.
Source code in m2dtools/basic/dynamics.py 
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def compute_SISF(positions, box_lengths, lag_array, k, num_vectors=100, n_repeat=100):
    """Compute the self-intermediate scattering function (SISF).

    Parameters
    ----------
    positions : np.ndarray
        Atom positions with shape ``(n_frames, n_atoms, 3)``.
    box_lengths : np.ndarray or list
        Simulation box lengths per frame with shape ``(n_frames, 3)``.
    lag_array : np.ndarray or list
        Lag times (in frames) at which to evaluate the SISF.
    k : float
        Magnitude of the wave vector ``|k|``.
    num_vectors : int, optional
        Number of random ``k``-vectors to sample. Default is ``100``.
    n_repeat : int, optional
        Number of time origins to average over. Default is ``100``.

    Returns
    -------
    np.ndarray
        Self-intermediate scattering values for each lag time.
    """
    n_frames, n_atoms, _ = positions.shape

    # Unwrap positions considering periodic boundary conditions
    unwrapped_pos = np.zeros_like(positions)
    unwrapped_pos[0] = positions[0]

    for t in range(1, n_frames):
        box_length_tmp = (box_lengths[t] + box_lengths[t - 1]) / 2
        delta = positions[t] - unwrapped_pos[t - 1]
        delta -= box_length_tmp * np.round(delta / box_length_tmp)
        unwrapped_pos[t] = unwrapped_pos[t - 1] + delta

    sisf = np.zeros(len(lag_array))
    vectors=random_vectors(k, num_vectors)

    for i, lag in enumerate(lag_array):
        if len(positions)-lag < n_repeat:
            displacements = unwrapped_pos[lag:] - unwrapped_pos[:-lag]
            #compute cos(k*r) for each random vector and average
            cos_kr = np.cos(np.einsum('ij,tkj->tki', vectors, displacements))
        else:
            random_indices = np.random.choice(len(positions)-lag, n_repeat)
            displacements = unwrapped_pos[lag + random_indices] - unwrapped_pos[random_indices]
            cos_kr = np.cos(np.einsum('ij,tkj->tki', vectors, displacements))
        sisf[i] = np.mean(cos_kr)
    return sisf

compute_SISF_zresolved(positions, box_lengths, lag_array, k, z_ranges, num_vectors=100, n_repeat=100)

Layer-resolved SISF with z-bin assignment performed at every time origin t0.

Parameters:

Name Type Description Default
z_ranges list of tuple

List of (zmin, zmax) tuples defining the layers along z-axis.

 required
Other
required

Returns:

Type Description
ndarray

Self-intermediate scattering values for each layer and lag time.
Source code in m2dtools/basic/dynamics.py 
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
def compute_SISF_zresolved(positions, box_lengths, lag_array, k,
                            z_ranges,
                            num_vectors=100, n_repeat=100) -> np.ndarray:
    """
    Layer-resolved SISF with z-bin assignment performed at every time origin t0.

    Parameters
    ----------
    z_ranges : list of tuple
        List of (zmin, zmax) tuples defining the layers along z-axis.
    Other parameters are the same as in `compute_SISF`.

    Returns
    -------
    np.ndarray
        Self-intermediate scattering values for each layer and lag time.
    """

    n_frames, n_atoms, _ = positions.shape
    n_layers = len(z_ranges)

    # ---- unwrap ----
    unwrapped = np.zeros_like(positions)
    unwrapped[0] = positions[0]

    for t in range(1, n_frames):
        box_tmp = 0.5 * (box_lengths[t] + box_lengths[t-1])
        delta = positions[t] - unwrapped[t-1]
        delta -= box_tmp * np.round(delta / box_tmp)
        unwrapped[t] = unwrapped[t-1] + delta

    # ---- random k vectors ----
    vectors = random_vectors_2D(k, num_vectors)

    # ---- output ----
    sisf_layers = np.zeros((n_layers, len(lag_array)))

    # ---- compute for each lag ----
    for i, lag in enumerate(lag_array):

        max_origin = n_frames - lag
        if max_origin < n_repeat:
            origins = np.arange(max_origin)
        else:
            origins = np.random.choice(max_origin, n_repeat, replace=False)

        # accumulator for each layer
        accum = [ [] for _ in range(n_layers) ]

        for t0 in origins:
            z_now = unwrapped[t0, :, 2]   # z positions at time origin

            # assign atoms at this time origin
            atom_layers = []
            for (zmin, zmax) in z_ranges:
                atom_layers.append(np.where((z_now >= zmin) & (z_now < zmax))[0])

            # compute displacements from t0 to t0+lag
            disp = unwrapped[t0 + lag] - unwrapped[t0]   # (N,3)

            for L, atoms_L in enumerate(atom_layers):
                if len(atoms_L) == 0:
                    continue

                d = disp[atoms_L]           # (n_atoms_L, 3)

                cos_kr = np.cos(np.einsum("ij,nj->ni", vectors, d))
                accum[L].append(np.mean(cos_kr))

        # average for each layer
        for L in range(n_layers):
            if len(accum[L]) > 0:
                sisf_layers[L, i] = np.mean(accum[L])

    return sisf_layers

compute_msd_pbc(positions, box_lengths, lag_array)

Compute mean-squared displacement with periodic boundaries.

Parameters:

Name Type Description Default
positions ndarray

Atom positions with shape (n_frames, n_atoms, 3).

 required
box_lengths ndarray or list

Simulation box lengths per frame with shape (n_frames, 3).

 required
lag_array ndarray or list

Lag times (in frames) at which to evaluate the MSD.

 required

Returns:

Type Description
ndarray

Mean-squared displacement values for each lag time.
Source code in m2dtools/basic/dynamics.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def compute_msd_pbc(positions, box_lengths, lag_array):
    """Compute mean-squared displacement with periodic boundaries.

    Parameters
    ----------
    positions : np.ndarray
        Atom positions with shape ``(n_frames, n_atoms, 3)``.
    box_lengths : np.ndarray or list
        Simulation box lengths per frame with shape ``(n_frames, 3)``.
    lag_array : np.ndarray or list
        Lag times (in frames) at which to evaluate the MSD.

    Returns
    -------
    np.ndarray
        Mean-squared displacement values for each lag time.
    """
    n_frames, n_atoms, _ = positions.shape

    # Unwrap positions considering periodic boundary conditions
    unwrapped_pos = np.zeros_like(positions)
    unwrapped_pos[0] = positions[0]

    for t in range(1, n_frames):
        box_length_tmp = (box_lengths[t] + box_lengths[t - 1]) / 2
        delta = positions[t] - unwrapped_pos[t - 1]
        delta -= box_length_tmp * np.round(delta / box_length_tmp)
        unwrapped_pos[t] = unwrapped_pos[t - 1] + delta

    msd = np.zeros(len(lag_array))

    for i, lag in enumerate(lag_array):
        displacements = unwrapped_pos[lag:] - unwrapped_pos[:-lag]
        squared_displacements = np.sum(displacements**2, axis=2)
        msd[i] = np.mean(squared_displacements)

    return msd