Cross-posted on MMSE (Matter Modeling Stack Exchange).
The following are the movements of the center of mass of a polymer chain over time in a monte carlo simulation.
# cm-X cm-Y cm-Z
9.507 21.232 9.910
9.092 18.427 7.308
7.994 14.856 4.675
0.533 12.138 -1.163
-0.056 6.735 -5.470
-2.138 4.950 -7.280
0.736 13.076 -10.012
9.516 14.611 -12.470
0.886 17.235 -10.714
-8.954 21.381 -11.735
-7.457 16.042 -16.710
-0.842 13.984 -17.874
-2.799 13.195 -14.247
0.738 14.792 -11.774
5.760 9.218 -10.846
8.406 13.735 -15.946
6.848 12.679 -9.621
13.569 14.036 -6.348
8.665 12.724 -8.238
... ... ... ... ...
I want to draw a log-log plot of MSD (mean square displacement) versus t of a movement of the polymer chain to examine the diffusion behavior.
The problem I am facing is that the plot should be a near-straight line, whereas my plot is not.
What am I doing incorrectly?
Without smoothing:
With smooting:
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
# Function to compute moving average
def moving_average(x, w):
return np.convolve(x, np.ones(w), 'valid') / w
# Read data from file, skipping the first row (header)
data = np.loadtxt('cm.dat', skiprows=1)
# Wrap CM values if they exceed -1000 or +1000
data = np.mod(data + 1000, 2000) - 1000 # performs the wrap operation
# Initialize reference point
x0, y0, z0 = data[0]
# Compute squared displacement for each time step
SD = [(x - x0)**2 + (y - y0)**2 + (z - z0)**2 for x, y, z in data]
# Compute the cumulative average of SD to get MSD at each time step
MSD = np.cumsum(SD) / np.arange(1, len(SD) + 1)
# Size of the moving window for smoothing
window_size = 100
# Apply moving average smoothing to MSD
MSD_smoothed = moving_average(MSD, window_size)
# Generate time steps for smoothed MSD
t_smoothed = np.arange(window_size, len(MSD) + 1)
# Create a log-log plot of smoothed MSD versus t
plt.figure(figsize=(8, 6))
plt.loglog(t_smoothed, MSD_smoothed, marker='o')
plt.title('Smoothed Mean Squared Displacement vs Time')
plt.xlabel('Time step')
plt.ylabel('Smoothed MSD')
plt.grid(True, which="both", ls="--")
plt.savefig('msd_plot_smoothed.png')
