Numpy计算平均平方位移

在本文中，我们将介绍如何使用Python和Fast Fourier Transform（快速傅里叶变换）计算分子动力学中的平均平方位移（Mean Square Displacement，MSD）。平均平方位移是研究分子运动的重要物理量，它可以反映分子在空间中的运动轨迹和速度。

阅读更多：Numpy 教程

什么是平均平方位移？

平均平方位移是用来衡量分子在时间t内由于热运动导致的空间离散程度。在分子动力学中，我们通常使用以下公式来计算平均平方位移：

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

其中， $r_i(t)$ 是第 $i$ 个粒子在时间 $t$ 时的位置， $N$ 是粒子的数量。这个公式表示在 $t$ 时间内，每个粒子移动的距离的平方的平均值。

如何计算平均平方位移？

在计算平均平方位移时，我们通常使用快速傅里叶变换（FFT）来将时间域的数据转换成频域的数据。首先，我们需要从模拟数据中获取每个粒子的坐标数据，然后计算相邻两个瞬时的平移距离。接着，通过FFT计算出平移距离的功率谱密度函数。最后，将功率谱密度函数转换回时间域，并取平均值即可得到平均平方位移。

以下是利用Python和numpy计算平均平方位移的示例代码：

import numpy as np

# 获取坐标数据
coords = np.loadtxt('coords.txt')

# 计算平移距离
displacements = np.diff(coords, axis=0)

# 计算每个粒子的平均平方位移
msd = []
for i in range(1, len(displacements)):
    d = displacements[:i]
    msd.append(np.mean(np.sum(np.square(d), axis=1)))
msd = np.array(msd)

# 计算功率谱
psd = np.abs(np.fft.fft(msd))**2

# 取平均值
freq = np.fft.fftfreq(msd.size)
idx = np.argsort(freq)
freq = freq[idx]
psd = psd[idx]
mask = freq > 0
freq = freq[mask]
psd = psd[mask]
slope, intercept = np.polyfit(np.log10(freq), np.log10(psd), 1)

# 绘制结果图
import matplotlib.pyplot as plt
plt.loglog(freq, psd)
plt.loglog(freq, 10**(intercept) * freq**(slope), 'r')
plt.show()

print('Slope:', slope)