任何人都可以给我一个关于使用 python 计算均方位移的建议吗？

Question

我正在研究分子动力学模拟。由于身边有很多好的MD包，所以我从来没有自己写过自己的代码。但是我面临着我需要编写自己的代码来计算均方位移。以下是我所做的程序。

1. 读取轨迹 a) 读取框尺寸 b) 读取原子的坐标 c) 解开原子的坐标

2. 计算 MSD 我已经编写了用于计算 MSD 的代码，如下所示。但是，它与 lammps 内置的 MSD 计算结果不匹配。谁能给我一个建议，出了什么问题？提前致谢。

    msd = zeros((Nsteps, Natoms), float)
    msd_x = zeros((Nsteps, Natoms), float)
    msd_y = zeros((Nsteps, Natoms), float)
    msd_z = zeros((Nsteps, Natoms), float)
    msdTotX = zeros(Nsteps, float)
    msdTotY = zeros(Nsteps, float)
    msdTotZ = zeros(Nsteps, float)
    msdTot = zeros(Nsteps, float)
    cnt = zeros(Nsteps, int)
    skip = 10
    for i in range(Natoms):
        if i % 10 == 0:
            print(f"Calculating MSD of {i+1}/{Natoms} is done")
        for j in range(0,Nsteps-1,skip):
            xj = crd_unwrap[j,i,0]
            yj = crd_unwrap[j,i,1]
            l = 0
            for k in range(j+1,Nsteps, skip):
                xk = crd_unwrap[k,i,0]
                yk = crd_unwrap[k,i,1]
                dx = xj - xk
                dy = yj - yk
                msd_x[k-j, i] = msd_x[k-j, i] + dx**2
                msd_y[k-j, i] = msd_y[k-j, i] + dy**2
                msd[k-j,i] = msd[k-j,i] + dx**2 + dy**2
                cnt[k-j] += 1
        for j in range(1,Nsteps-1, skip):
            N = cnt[j] 
            N = float(N)
            msdTot[j] += msd[j, i]/N
            msdTotX[j] += msd_x[j, i]/N
            msdTotY[j] += msd_y[j, i]/N
            if msd[j,i] == 0:  
                continue
   
    out = open('mymsd.txt','w')
    for i in range(1, Nsteps-1, skip):
        xmsd = msdTotX[i]/Natoms
        ymsd = msdTotY[i]/Natoms
        msd_tot = msdTot[i]/Natoms
        out.write('%f %f %f %f\n'%(i+1, xmsd, ymsd, msd_tot))
    out.close()
    print("Done")
   

Answer 1

我会将此作为评论，但没有足够的声誉。你可能想看看这篇文章： Mean square displacement python