Measure of how far particles (monomers, polymers, colloids) have moved in a given $\Delta t$.

environment.yml

name: hoomd-workshop
channels:
  - conda-forge
  - anaconda
dependencies:
  - python=3.11
  - hoomd=3.9.*=*cpu*
  - freud
  - matplotlib
  - fresnel
  - gsd
  - ipykernel

Calculation and plotting of mean squared displacement with the freud package - faster and more reliable code incuding unwrapping of trajectories:

import freud 
import matplotlib.pyplot as plt
import numpy as np 
import gsd, gsd.hoomd
 
 
 
fig = plt.figure()
ax = plt.axes()
 
 
trajectory = gsd.hoomd.open("KG_melt_N_300_period_1000.0.gsd",'r')
nparticles = len(trajectory[0].particles.position)
nframes = len(trajectory)
positions = np.zeros(shape=(nframes,nparticles,3))
images =  np.zeros(shape=(nframes,nparticles,3))
 
all_unwrapped_positions = []
all_unwrapped_coms = []
all_unwrapped_coms_particles = []
times = []
 
box  = trajectory[0].configuration.box[0:3]
bonds =trajectory[0].bonds.group 
polymer_delimiters = bonds[np.where((bonds[:-1,1]-bonds[1:,0])!=0)[0]+1][:,0]
time0 = trajectory[0].configuration.step
 
for frame in trajectory[10:]:
 
    unwrapped_positions = frame.particles.position + frame.particles.image*box
    unwrapped_positions = unwrapped_positions - np.mean(unwrapped_positions,axis=0)
    all_unwrapped_positions.append(unwrapped_positions)
    polymers = np.split(unwrapped_positions,polymer_delimiters)
    unwrapped_coms = []
 
    for p in polymers:
        com = np.mean(p,axis=0)
        unwrapped_coms.append(com)
 
    unwrapped_coms_tiled = np.repeat(unwrapped_coms,len(p),axis=0)
    all_unwrapped_coms_particles.append(unwrapped_positions-unwrapped_coms_tiled)
 
    all_unwrapped_coms.append(unwrapped_coms)
    times.append(frame.configuration.step-time0)
 
all_unwrapped_positions = np.array(all_unwrapped_positions)  # shape = (nframes,nparticles,3)
all_unwrapped_coms = np.array(all_unwrapped_coms)            # shape = (nframes,npolymers,3)
all_unwrapped_coms_particles = np.array(all_unwrapped_coms_particles)  # shape = (nframes,nparticles,3)
 
 
nframes= len(trajectory)
msd_calculator = freud.msd.MSD(mode='window')
msd_calculator.compute(all_unwrapped_positions)
result = msd_calculator.particle_msd
msd1 = np.nanmean(result,axis=1)
 
# average over innermost monomer
polymer_middles = np.array(polymer_delimiters-polymer_delimiters[0]/2.).astype(int)
msd1_middle = np.nanmean(result[:,polymer_middles],axis=1)
msd1_ends = np.nanmean(result[:,polymer_delimiters],axis=1)
 
 
msd_calculator.compute(all_unwrapped_coms_particles)
result = msd_calculator.particle_msd
msd2 = np.nanmean(result,axis=1)
 
msd_calculator.compute(all_unwrapped_coms)
result = msd_calculator.particle_msd
msd3 = np.nanmean(result,axis=1)
 
timestep = 0.005 # this needs to be the actual time between frames
period =1000
lagtimes = np.arange(len(msd1))*period*timestep # make the lag-time axis
 
    # plot the actual MSD
 
#for i in range(result.shape[1]):
#    ax.plot(lagtimes, result[:,i])
msd1[0]=0
msd2[0]=0
msd3[0]=0
msd1_middle[0]=0
msd1_ends[0]=0
 
msd1 = msd1[1:]
msd2 = msd2[1:]
msd3 = msd3[1:]
msd1_ends = msd1_ends[1:]
msd1_middle = msd1_middle[1:]
lagtimes = lagtimes[1:]
 
np.savetxt('freud2.txt', np.c_[lagtimes,msd1,msd2,msd3]) 
ax.plot(lagtimes, msd1)
ax.plot(lagtimes, msd2)
ax.plot(lagtimes, msd3)
 
ax.plot(lagtimes, msd1_middle)
ax.plot(lagtimes, msd1_ends)
 
 
ax.set_yscale('log')
ax.set_xscale('log')
 
plt.savefig('plot_KG_melt_MSD.png')
plt.show()

Straightforward calculation according to mathematical definition:

def msd_straight_forward(r):
    R"""
    Default way of calculating mean square displacement
 
    args: r (array) positions of particle/com over time. This only works
    if the timesteps are evenly spaced.
 
    returns: mean square displacement (array)
 
    """
    shifts = np.arange(len(r))
    msds = np.zeros(shifts.size)
 
    for i, shift in enumerate(shifts):
        diffs = r[:-shift if shift else None] - r[shift:]
        sqdist = np.square(diffs).sum(axis=1)
        msds[i] = sqdist.mean()
 
    return msds

More efficient implementation using Fourier transform:

def autocorrFFT(x):
    R"""
    caclulates autocorrelation function of signal
 
    args: x (array) signal
    returns: autocorrelation function (array)
 
    There are different convetions to define autocorrelation,
    convention A is the wikipedia one <https://en.wikipedia.org/wiki/Autocorrelation>
    convention B is multiplied by (N-m)
 
    The Wiener Khinchin theorem says that the  power spectral density (PSD)
    of a function is the Fourier transform of the autocorrelation. See
    <https://en.wikipedia.org/wiki/Wiener%E2%80%93Khinchin_theorem>
 
    So, compute PSD of x and back fourier transform it to get the
    autocorrelation (in convention B). It's the cyclic autocorrelation, so
    zero-padding is needed to get the non-cyclic autocorrelation.
 
    """
    N=len(x)
    F = np.fft.fft(x, n=2*N)  #2*N because of zero-padding
    PSD = F * F.conjugate()
    res = np.fft.ifft(PSD)
    res= (res[:N]).real   #now we have the autocorrelation in convention B
    n=N*np.ones(N)-np.arange(0,N) #divide res(m) by (N-m)
    return res/n #this is the autocorrelation in convention A
 
def msd_fft(r):
    R"""
    calculates mean square displacement
 
    args: r (array) positions of particle/com over time.
    returns: mean square displacement (array)
 
    Algorithm as described in section 4.2 and 4.1 of
    'nMoldyn - Interfacing spectroscopic experiments, molecular dynamics
    simulations and models for time correlation functions,
    V. Calandrini, E. Pellegrini, P. Calligari, K. Hinsen, G.R. Kneller
    <https://doi.org/10.1051/sfn/201112010>'
 
    This algorithm works only if the timesteps are evenly spaced.
    It is faster thant the straightforward implementation and yields the
    same results.
 
    """
    N=len(r)
    D=np.square(r).sum(axis=1)
    D=np.append(D,0)
    S2=sum([autocorrFFT(r[:, i]) for i in range(r.shape[1])])
    Q=2*D.sum()
    S1=np.zeros(N)
    for m in range(N):
        Q=Q-D[m-1]-D[N-m]
        S1[m]=Q/(N-m)
    return S1-2*S2

Read in the gsd file similarly to how it is described in 'analyzing gsd trajectories', then calculate msd. The difference is that the entire trajectory (all positions at all times) is needed for msd calculations. The following code calculates the mean-squared displacement for each monomer, each center of mass of each chain, and the relative msd of monomer to the center of mass of the respective chain.

parser = argparse.ArgumentParser(description='get center of mass msd from gsd file for a linear polymer chain configuration')
 
non_opt = parser.add_argument_group('mandatory arguments')
non_opt.add_argument('-i','--input', metavar="<gsd file>", type=str, required=True,dest='input_file',
    help="input gsd file")
 
parser.add_argument('-o','--output', metavar='<float>', type=str, required=True,dest='output_file',
                   help='output file, DEFAULT=""')
 
 
args = parser.parse_args()
output_file = args.output_file
input_file = args.input_file
 
 
in_file=gsd.fl.GSDFile(name=input_file, mode='rb')
N_frames=in_file.nframes
trajectory  = gsd.hoomd.HOOMDTrajectory(in_file)
 
# can be used to exclude frames, i.e. first ones, equilibration etc.
to_take_all=[np.arange(0,N_frames,1)]
 
# this can be simplified... 
for q in range(len(to_take_all)):
    frames_to_take=to_take_all[q]
    time=[]
    # array for all center of masses for all polymers at all times
    C=[[] for i in range(len(frames_to_take))]
    # array for all positions
    Pos=np.zeros((len(frames_to_take),trajectory[0].particles.N,3))
    Pos_minus_com=np.zeros((len(frames_to_take),trajectory[0].particles.N,3))
    image_checksum=0
 
    for k,j in enumerate(frames_to_take):
        sys.stdout.write("\rCurrent frame %3d/%3d"%(k,N_frames))
        sys.stdout.flush()
        # get positions and unwrap them
        b=trajectory[j].configuration.box[0:3]
        pos=trajectory[j].particles.position+trajectory[j].particles.image*b
        Pos[k]=pos
        # create checksum for images - if images are not written out, msd will be wrong
        # todo: find better way of doing this
        image_checksum += np.sum(trajectory[j].particles.image)
        # get bonds and figure out which particles are in a polymer
        polymers = connected_components(trajectory[j].bonds.group)
        pos_minus_com=[]
        # calculate center of mass for each polymer
        for pol in polymers:
            ids=np.unique(pol)
            com = np.average(pos[ids],axis=0)
            C[k].append(com)
            pos_minus_com.append(pos[ids] - com)
 
        Pos_minus_com[k] = np.vstack(pos_minus_com)
        time.append(trajectory[j].configuration.step)
 
    # sanity checks
    # check that images were written out into gsd file
    if  image_checksum==0:
        print("""WARNING: either no particle ever left the box or
        no images were written out in the gsd file - this can lead
        to a wrong msd calculation""")
 
    # check that timesteps are uniformly spaced - if not, this determination of msd won't work
    if len(np.unique(np.diff(time)))!=1:
        print("""ERROR: timesteps are not evenly spaced!
        mean square displacement is determined with fourier
        transform method, only works with evenly spaced timesteps""")
        exit(1)
 
    # acutal calculation of msd COM
    C=np.array(C)
    corr=np.zeros((C.shape[1],len(time)))
    i=0
    for r in np.rollaxis(C, 1):
        p=msd_fft(r)
        corr[i]=p
        sys.stdout.write("\rCurrent COM %3d/%3d    "%(i,C.shape[1]))
        sys.stdout.flush()
        i=i+1
    mean_x=np.mean(corr,axis=0)
    f1=open(output_file+'_'+str(q)+"_com.msd", 'w')
    f1.write("# mean square displacement center of mass polymer\n")
    f1.write("# timestep g_3(t) \n")
    i=0
    for i in range(len(mean_x)):
        f1.write("%s %f \n"%(str(time[i]-time[0]),mean_x[i]))
        i=i+1
 
    # particle msd relative to com
    C=np.array(Pos_minus_com)
    corr=np.zeros((C.shape[1],len(time)))
    i=0
    for r in np.rollaxis(C, 1):
        p=msd_fft(r)
        corr[i]=p
        sys.stdout.write("\rCurrent particle-com %3d/%3d    "%(i,C.shape[1]))
        sys.stdout.flush()
        i=i+1
    mean_x=np.mean(corr,axis=0)
    f1=open(output_file+'_'+str(q)+"_particle_com.msd", 'w')
    f1.write("# mean square displacement monomer-center of mass\n")
    f1.write("# timestep g_2(t) \n")
    i=0
    for i in range(len(mean_x)):
        f1.write("%s %f \n"%(str(time[i]-time[0]),mean_x[i]))
        i=i+1
 
    # particle msd
    C=np.array(Pos)
    corr=np.zeros((C.shape[1],len(time)))
    i=0
    for r in np.rollaxis(C, 1):
        p=msd_fft(r)
        corr[i]=p
        sys.stdout.write("\rCurrent particle %3d/%3d    "%(i,C.shape[1]))
        sys.stdout.flush()
        i=i+1
    mean_x=np.mean(corr,axis=0)
    f1=open(output_file+'_'+str(q)+"_particle.msd", 'w')
    f1.write("# mean square displacement of monomer\n")
    f1.write("# timestep g_1(t) \n")
    i=0
    for i in range(len(mean_x)):
        f1.write("%s %f \n"%(str(time[i]-time[0]),mean_x[i]))
        i=i+1
 
sys.stdout.write('DONE\n')
sys.stdout.flush()

The code that identifies the linear bonded sub-sections that make up a chain can be found here.

• Kremer, Kurt, and Gary S. Grest. “Dynamics of entangled linear polymer melts: A molecular‐dynamics simulation.” The Journal of Chemical Physics 92.8 (1990): 5057-5086. doi/10.1063/1.458541