有人能给我一个关于使用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