Python计算绘制粒子轨迹MSD异常问题求助

问题描述 #

需要计算并绘制粒子轨迹的MSD，以此判断粒子是纯扩散、受约束还是受外力作用，但当前Python代码生成的MSD曲线始终为水平线，且不同数据集结果相同。

用户原始代码：

import numpy as np
import matplotlib.pyplot as plt

# Loading data from the txt file
data = np.loadtxt('/Users/arthurcoet/Documents/PhD COËT/CiNAM/Python/Exercices/Dtxt1.txt')

# Retrieving data columns
temps = data[:, 0]
indices = data[:, 1]
positions_x = data[:, 2]
positions_y = data[:, 3]

# Calculation of the MSD for each particle
particle_indices = np.unique(indices)
msd_particles = []

for particle_index in particle_indices:
    indices_particle = np.where(indices == particle_index)[0]
    x_particle = positions_x[indices_particle]
    y_particle = positions_y[indices_particle]
    displacements = np.sqrt((x_particle - x_particle[0])**2 + (y_particle - y_particle[0])**2)
    msd_particle = np.mean(displacements ** 2)
    msd_particles.append(msd_particle)

# Calculation of average MSD
msd_mean = np.mean(msd_particles)

# Plot of average MSD against time
plt.plot(temps, msd_mean * np.ones_like(temps), '-')
plt.xlabel('Temps')
plt.ylabel('MSD')
plt.title('MSD moyen en fonction du temps')
plt.show()

数据格式：时间、粒子索引、X、Y，示例数据：

2.000000000000000042e-02    1.000000000000000000e+00    1.240807597782861649e+00    -3.785534313700077980e-01
2.000000000000000042e-02    2.000000000000000000e+00    -7.309480801080029400e-01   -7.170078639286830979e-01
2.000000000000000042e-02    3.000000000000000000e+00    8.098986703909789586e-01    5.815296326847829711e-01
2.000000000000000042e-02    4.000000000000000000e+00    2.750517661172900349e-03    5.522082763915119319e-01
2.000000000000000042e-02    5.000000000000000000e+00    9.127289184940016176e-01    5.289936832099206843e-01
2.000000000000000042e-02    6.000000000000000000e+00    1.340115717383320693e+00    -2.206167086146532119e-01
2.000000000000000042e-02    7.000000000000000000e+00    -9.347632084034850353e-01   -1.031270972401856945e+00
2.000000000000000042e-02    8.000000000000000000e+00    -1.564283425436366892e+00   9.906037061246882880e-01
2.000000000000000042e-02    9.000000000000000000e+00    -1.092736507472809315e-01   -7.019693105410058642e-01
2.000000000000000042e-02    1.000000000000000000e+01    1.079635002230163246e-02    1.936345833496776025e+00

粒子轨迹图：

错误分析 #

1. MSD计算逻辑错误：原代码对每个粒子计算的是「从初始位置到所有时间点位移平方的平均值」，再取所有粒子的平均得到一个单一数值，导致绘图时只能生成水平线。正确的MSD应该是随时间变化的函数，即每个时间点对应所有粒子在该时刻相对于初始位置的位移平方的平均值。
2. 未按时间排序轨迹：原代码提取粒子轨迹时未确保时间序列是递增的，可能导致位移计算混乱。
3. 绘图逻辑错误：用所有时间点乘以同一个均值，必然生成水平线。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt

# 加载数据
data = np.loadtxt('/Users/arthurcoet/Documents/PhD COËT/CiNAM/Python/Exercices/Dtxt1.txt')

# 提取列数据
temps = data[:, 0]
indices = data[:, 1]
positions_x = data[:, 2]
positions_y = data[:, 3]

# 获取所有唯一粒子索引和时间点
particle_indices = np.unique(indices)
unique_times = np.sort(np.unique(temps))
num_times = len(unique_times)

# 初始化MSD数组，每个时间点对应一个MSD值
msd_over_time = np.zeros(num_times)

for i, t in enumerate(unique_times):
    # 存储当前时间点所有粒子的位移平方
    displacements_sq = []
    for pid in particle_indices:
        # 获取该粒子的所有数据行
        particle_mask = (indices == pid)
        particle_data = data[particle_mask]
        # 按时间排序确保轨迹顺序正确
        particle_data_sorted = particle_data[np.argsort(particle_data[:, 0])]
        # 找到当前时间点的位置
        t_mask = particle_data_sorted[:, 0] == t
        if np.any(t_mask):
            # 初始位置（第一个时间点的位置）
            x0 = particle_data_sorted[0, 2]
            y0 = particle_data_sorted[0, 3]
            # 当前时间点的位置
            x_t = particle_data_sorted[t_mask, 2][0]
            y_t = particle_data_sorted[t_mask, 3][0]
            # 计算位移平方
            disp_sq = (x_t - x0)**2 + (y_t - y0)**2
            displacements_sq.append(disp_sq)
    # 计算当前时间点的平均MSD
    if displacements_sq:
        msd_over_time[i] = np.mean(displacements_sq)

# 绘制MSD随时间变化的曲线
plt.plot(unique_times, msd_over_time, '-')
plt.xlabel('Temps')
plt.ylabel('MSD')
plt.title('MSD moyen en fonction du temps')
plt.grid(True)
plt.show()

结果说明 #

修正后的代码会生成随时间变化的MSD曲线：

• 若曲线呈线性增长（MSD ∝ t），说明粒子是纯扩散运动；
• 若曲线趋于平稳，说明粒子受约束；
• 若曲线呈超线性增长（MSD ∝ t^α, α>1），说明粒子受外力作用。

内容的提问来源于stack exchange，提问作者ArthurCo