diff --git a/package/AUTHORS b/package/AUTHORS index fface7c624..b79557763e 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -258,7 +258,7 @@ Chronological list of authors - Abdulrahman Elbanna - Tulga-Erdene Sodjargal - Gareth Elliott - + - Sirsha Ganguly External code ------------- diff --git a/package/MDAnalysis/analysis/msd.py b/package/MDAnalysis/analysis/msd.py index 8d0e7a39e1..85935519b4 100644 --- a/package/MDAnalysis/analysis/msd.py +++ b/package/MDAnalysis/analysis/msd.py @@ -1,3 +1,25 @@ +######## Documentation +# +# Issue: The default '_conclude_simple()' function previously calculated time-averaged mean squared displacement considering linear timespacing between the frames even when a dump file with non-linear timestep gap is provided +# +# Modified by Sirsha Ganguly +# +# Modification: +# '_conclude_modified()' is the new function added to calculate time-averaged mean squared displacement of non-linear dumps + # Note: This new function is generalized for non-linear dumps. + # Moreover, if you use this function to calculate for linear timedump this gives same result as '_conclude_simple()' which is implemented in the default version. + # So, I believe the algorithm in '_conclude_modified()' can potentially be used to calculate time-averaged msd in future. + # For now it only executes when the dump is non-linear. [ This is determined by the added function is_linear() ] + # It is tested with a .dump file of a Coarse Grained system of single type +# +# General Algorithm of '_conclude_modified()' to calculate time-averaged msd: + # It creates a dictionary/hash-map with the key being the time difference between the frames (Delta t) and value being a list containing the frame to frame MSD values [msd1, msd2, .....] for that particular Delta t + # Dictionary to collect MSDs: {Δt: [msd1, msd2, ...]} + # The reference frame is changed with each iteration and each time a new Key (i.e, Delta t) is found it is appended in the dictionary/hash-map + # If a duplicate Delta_t is found it is added to the value (msd list) of that specefic Key in the dictionary/hash-map + # Lastly in the dictionary/hash-map for each Delta_t the msd values inside the list are averaged +# + # -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- # vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 # @@ -263,7 +285,17 @@ del Doi + class EinsteinMSD(AnalysisBase): + + @staticmethod + def is_linear(sequence): # This function is to check if the timesteps in the array is linear or not! If it's linear default MDA code will carry-out other wise the modified code carries! + if len(sequence) < 2: + return True + step = sequence[1] - sequence[0] + return all(sequence[i+1] - sequence[i] == step for i in range(len(sequence) - 1)) + + r"""Class to calculate Mean Squared Displacement by the Einstein relation. Parameters @@ -386,7 +418,10 @@ def _conclude(self): if self.fft: self._conclude_fft() else: - self._conclude_simple() + if self.is_linear([ts.time for ts in self._trajectory]): + self._conclude_simple() + else: + self._conclude_modified() # New modified code def _conclude_simple(self): r"""Calculates the MSD via the simple "windowed" algorithm.""" @@ -421,3 +456,46 @@ def _conclude_fft(self): # with FFT, np.float64 bit prescision required. positions[:, n, :] ) self.results.timeseries = self.results.msds_by_particle.mean(axis=1) + + + def _conclude_modified(self): + from collections import defaultdict + print("The Dump file has non-linear time spacing") + + dump_times = [ts.time for ts in self._trajectory] + # print(dump_times) + n_frames = len(dump_times) + n_atoms = self.n_particles + positions = np.zeros((n_frames, n_atoms, 3)) + positions = self._position_array.astype(np.float64) + + msd_dict = defaultdict(list) # Dictionary to collect MSDs: {Δt: [msd1, msd2, ...]} + + # Looping over all the frames as if the referenced gets shifted frame to frame + for i in range(n_frames): + for j in range(i + 1, n_frames): + delta_t = dump_times[j] - dump_times[i] + + # Compute displacement and squared displacement + disp = positions[j] - positions[i] + squared_disp = np.sum(disp ** 2, axis=1) + msd = np.mean(squared_disp) + + # Store MSD under corresponding Δt + msd_dict[delta_t].append(msd) + + msd_dict[0] = [0] + + # Prepare averaged results + delta_t_values = [] + avg_msds = [] + for delta_t in sorted(msd_dict.keys()): + msd_list = msd_dict[delta_t] + avg_msd = np.mean(msd_list) + delta_t_values.append(delta_t) + avg_msds.append(avg_msd) + + self.results.timeseries = delta_t_values + self.results.msds_by_particle = avg_msds + + # self.results.timeseries = self.results.msds_by_particle.mean(axis=1) \ No newline at end of file