From 045f91fa8ea6ed19903f7e8ebaca670cd4c82a4e Mon Sep 17 00:00:00 2001 From: Sirsha Ganguly Date: Fri, 13 Jun 2025 17:43:41 -0400 Subject: [PATCH 1/3] Inditial addition of code to calculate time-averaged msd for non-linear dump --- package/MDAnalysis/analysis/msd.py | 81 +++++++++++++++++++++++++++++- 1 file changed, 80 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/analysis/msd.py b/package/MDAnalysis/analysis/msd.py index 8d0e7a39e1..1746a6c15e 100644 --- a/package/MDAnalysis/analysis/msd.py +++ b/package/MDAnalysis/analysis/msd.py @@ -263,7 +263,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,9 +396,13 @@ 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): + print("The Dump file has linear time spacing") # optional comment I wrote to check if the is_linear flag working correctly r"""Calculates the MSD via the simple "windowed" algorithm.""" lagtimes = np.arange(1, self.n_frames) positions = self._position_array.astype(np.float64) @@ -421,3 +435,68 @@ def _conclude_fft(self): # with FFT, np.float64 bit prescision required. positions[:, n, :] ) self.results.timeseries = self.results.msds_by_particle.mean(axis=1) + + + + + + + # New function to calculate in the following + # + # Note: This 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 this following algorithm can potentially be used to calculate time-averaged msd in future. + # For now it only executes when the dump is non-linear. + # It is tested with a .dump file of a Coarse Grained system of single type + + def _conclude_modified(self): + from collections import defaultdict + import matplotlib.pyplot as plt + 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) + + print(msd_dict) + 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) + print(f"Δt = {delta_t}: Averaged MSD = {avg_msd} (from {len(msd_list)} pairs)") + + # Plot Δt vs averaged MSD # Optional + plt.figure(figsize=(8, 5)) + plt.plot(delta_t_values, avg_msds,marker='o', linestyle='-', label='Averaged MSD') + + plt.xlabel('Δt (Time Difference)') + plt.ylabel('Mean Squared Displacement (MSD)') + plt.title('MSD Averaged Over Frame Pairs with Same Δt') + plt.grid(True) + plt.legend() + plt.tight_layout() + plt.show() From 3462365b18f6dcd08c43b920a62020088b334d7f Mon Sep 17 00:00:00 2001 From: Sirsha Ganguly Date: Fri, 13 Jun 2025 19:23:43 -0400 Subject: [PATCH 2/3] updated_documentation_v1_plots_removed --- package/MDAnalysis/analysis/msd.py | 55 +++++++++++++++--------------- 1 file changed, 27 insertions(+), 28 deletions(-) diff --git a/package/MDAnalysis/analysis/msd.py b/package/MDAnalysis/analysis/msd.py index 1746a6c15e..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 # @@ -402,7 +424,6 @@ def _conclude(self): self._conclude_modified() # New modified code def _conclude_simple(self): - print("The Dump file has linear time spacing") # optional comment I wrote to check if the is_linear flag working correctly r"""Calculates the MSD via the simple "windowed" algorithm.""" lagtimes = np.arange(1, self.n_frames) positions = self._position_array.astype(np.float64) @@ -437,21 +458,8 @@ def _conclude_fft(self): # with FFT, np.float64 bit prescision required. self.results.timeseries = self.results.msds_by_particle.mean(axis=1) - - - - - # New function to calculate in the following - # - # Note: This 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 this following algorithm can potentially be used to calculate time-averaged msd in future. - # For now it only executes when the dump is non-linear. - # It is tested with a .dump file of a Coarse Grained system of single type - def _conclude_modified(self): from collections import defaultdict - import matplotlib.pyplot as plt print("The Dump file has non-linear time spacing") dump_times = [ts.time for ts in self._trajectory] @@ -476,7 +484,6 @@ def _conclude_modified(self): # Store MSD under corresponding Δt msd_dict[delta_t].append(msd) - print(msd_dict) msd_dict[0] = [0] # Prepare averaged results @@ -487,16 +494,8 @@ def _conclude_modified(self): avg_msd = np.mean(msd_list) delta_t_values.append(delta_t) avg_msds.append(avg_msd) - print(f"Δt = {delta_t}: Averaged MSD = {avg_msd} (from {len(msd_list)} pairs)") - - # Plot Δt vs averaged MSD # Optional - plt.figure(figsize=(8, 5)) - plt.plot(delta_t_values, avg_msds,marker='o', linestyle='-', label='Averaged MSD') - - plt.xlabel('Δt (Time Difference)') - plt.ylabel('Mean Squared Displacement (MSD)') - plt.title('MSD Averaged Over Frame Pairs with Same Δt') - plt.grid(True) - plt.legend() - plt.tight_layout() - plt.show() + + 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 From f6e5ddfbc6be4aa975a8872c768c1be988092785 Mon Sep 17 00:00:00 2001 From: Sirsha Ganguly Date: Sun, 15 Jun 2025 14:07:24 -0400 Subject: [PATCH 3/3] Added_name_in_the_package/AUTHORS --- package/AUTHORS | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 -------------