Numpy计算平均平方位移（Mean Square Displacement）的Python和FFT

在物理学和化学中，Mean Square Displacement (MSD) 是一个被广泛应用的指标，用于衡量颗粒运动的随时间变化。MSD通常用来研究材料的动力学，催化剂的动力学特性，以及分子在液体中的扩散特性等。

在本文中，我们将介绍如何使用Python中的Numpy库，结合FFT快速傅里叶变换，来计算平均平方位移（Mean Square Displacement，MSD）。

阅读更多：Numpy 教程

什么是平均平方位移（Mean Square Displacement）

平均平方位移（Mean Square Displacement，MSD）是记录粒子在一段时间内运动路径长度的指标。MSD的定义为：

$MSD(t) = \frac{1}{N} \sum_{i=1}^{N} \left( r_i(t) – r_i(0) \right)^2$

其中， $r_i(t)$ 是粒子在时间t时的位置， $r_i(0)$ 是粒子在时间0时的位置， $N$ 是粒子数目。

在粒子采用布朗运动的情况下，MSD与时间t的关系为：

$MSD(t)=<R^2>=4Dt$

其中， $<R^2>$ 表示取所有的粒子平均后的平方位移， $D$ 称作扩散系数。

如何使用Python计算平均平方位移

首先，我们需要生成一个具有随机运动的粒子轨迹，用来模拟真实的运动数据。这个轨迹可以是一个二维的数组，其中每一行代表一个时间点的坐标。

在这里，我们使用Numpy的randn函数来生成一个二维的数组，用来表示粒子的随机起始位置。并使用cumsum函数求出每个时刻的位置。

import numpy as np

# set number of particle 'N' and time 'T'
N = 1000
T = 100

# initialize the position of the particle
pos = np.zeros((N, 2))
pos[0, :] = np.random.randn(2)

# genarate the simulation of Brownian motion
for i in range(1, N):
  pos[i, :] = pos[i - 1, :] + np.random.randn(2)

接下来，我们需要将坐标数据进行一些转化，才能够计算出MSD。

首先，我们需要将坐标数据沿着时间方向进行平移，然后计算出每个时间变化的时间差。这个时间差的取值范围可以设置为从1到N/2。

# Time shift for each particle
deltas = np.arange(1, N // 2)
shifts = np.arange(T - N // 2)

# Calculte the MSD for each time shift and all the particles
msds = np.zeros((shifts.size, deltas.size))
for i, delta in enumerate(deltas):
  for j, shift in enumerate(shifts):
    diffs = pos[shift + delta:, :] - pos[shift:-delta, :]
    sqdist = np.square(diffs).sum(axis=1)
    msds[j, i] = sqdist.mean()

最后，我们将MSD与时间t的关系用图表来呈现出来，以更好的观察粒子的扩散特性。

# plot the square of the distance traveled as a function of time
import matplotlib.pyplot as plt

for i, delta in enumerate(deltas):
  msd = msds[:, i]
  plt.plot(shifts, msd, 'o-', label=f' $\Delta t = {delta}$ ')
plt.xlabel('Time  $t$ ')
plt.ylabel('MSD')
plt.legend()
plt.show()