Skip to content

Support for lag-time based calculations #34

Description

@amruthesht

Feature Request

Many scientifically relevant analyses require repeated access to data over a range of time points or simulation frames. The mathematical formulation behind most such analyses are time correlation-based computations (eg. autocorrelation). Typical examples of scientifically relevant quantities are MSD (mean-squared displacement), VACF (velocity autocorrelation) etc.

However, as described above, these analyses require repeated access to data over a range of time values. A MDA data structure like Universe would thus need access to data points at previous times—a limitation inherent to streaming.

Describe the solution you'd like

The solution would ideally require the Universe to store simulation frames from previous time points up to a certain user-defined number temporarily. This could be implemented via the BufferedTrajectory. However, unlike current batching strategy this would have to be implemented as FIFO Buffer with LIFO-esque consumption.

An example of the kind of analysis described above with rudimentary sample solution code, inspired by MDAnalysis' present mean-squared displacement routine.

import MDAnalysis as mda
u = mda.Universe("topol.tpr", "imd://localhost:8889")

if not isinstance(u.trajectory, BufferedTrajectory):
    u.trajectory = BufferedTrajectory(u.trajectory, n_lags+1) 
    # n_lags is the number of lag times for the calculation of MSD
    # could be done by making _maxlen initailization parameter in the class

msd_dict = collections.defaultdict(list)
msds_by_particle_dict = collections.defaultdict(list)

msd_dict[0] = [0]
msds_by_particle_dict[0.0] = [np.zeros(u.atoms.n_atoms)]

for ts in u.trajectory: # might wanna implemenet as seperate get_next_frame 
    # deque should already behave FIFO-esque with removing old frames when maxlen is reached
    for i in range(1, len(u.trajectory)): # BufferedTrajectory needs a `len` method
        delta_t = u.trajectory[0].time - u.trajectory[-i].time
        # actual calculation - would be different for VACF
        disp = u.trajectory[0].position - u.trajectory[-i].position
        squared_disp = np.sum(disp**2, axis=1)
        msd = np.mean(squared_disp)
        msd_dict[delta_t].append(msd)
        msds_by_particle_dict[delta_t].append(squared_disp)
    
    delta_t_values = sorted(msd_dict.keys())
    avg_msds = [np.mean(msd_dict[dt]) for dt in delta_t_values]
    msds_by_particle_array = np.zeros((len(delta_t_values), u.atoms.n_atoms))
    for idx, dt in enumerate(delta_t_values):
        arr = np.vstack(msds_by_particle_dict[dt])
        msds_by_particle_array[idx, :] = np.mean(arr, axis=0)

    # this averaging can be further optimized by reusing means calculated at previous ts 
    # and calculating new means using msd values from newly arrived ts at relevant delta_t_values
    
    # return avg_msds, delta_t_values, msds_by_particle_array

EDIT: Updated sample code to provided updating average results when a new frame is read into trajectory

Additional context

It might be useful to add buffered trajectory support with calculations modified in MDAnalysis to suit buffered and non-buffered universes in the future.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions