Python计算均方回转半径、质心、均方末端距

刚开始接触分子动力学（高分子相关），跑出数据后不会分析，虽然LAMMPS里compute命令能计算很多值，比如说Rg、RDF等等。但是按照导师的要求，还是需要从基本的公式，用编程的方法深入理解各项性质。我导主要用的是Fortran语言，而我考虑到目前python上手较为容易，且用的人也比较多，所以选择了学习Python语言。在看了老师的Fortran文件以及书里面关于回转半径、质心、末端距等高分子链基本性质的公式解释之后，开始编写我第一个Python程序。

实际上，此前我并没有接触过编程语言，所以这次用了一个星期的时间写出了这个程序（很多东西看着感觉还是有点繁琐）。希望各位大佬如果看到能帮忙指点指点！力求把程序改得更加简洁。

说明一下，我是用得Lammps进行的分子动力学模拟，得到的文件是dump后的轨迹文件。

以下是我的全部代码：

import numpy as np

np.set_printoptions(suppress=True)

global Rcmx,Rcmy,Rcmz,Rcm,Rg2x,Rg2y,Rg2z,Rg2,Re2x,Re2y,Re2z,Re2

### Set up some parameters related to the MD simulation
N =30                                # number of atoms
Nf = 10000                           # number of frames
dt = 0.007                           # timestep
fileName = "pe_30_k_5.record"        # read the dump lammpstrj file
sorted_fileName = "sorted_pe_pe_30_k_5.record"

### Sort
f = open(fileName, "r")
fs = open(sorted_fileName, "w")
for i in range(Nf):
    ratio = (i+1)/Nf *100
    ls = np.zeros((N, 5)) 
    print("Sorted Progress %s%%" %ratio)
    initial = i * (9+N)                         # The first 9 lines should be excluded
    for j in range(9):
        f.readline()
    for k in range(N):
        line = f.readline().split(" ")[0:5]     # 切片提取01234列的数据
        ls[k] = line
    idx = np.argsort(ls[:,0])
    sorted_ls = ls[idx]
    for row in sorted_ls:
        fs.write(",".join("%s" %id for id in row)+"\n")
    fs.close()
    fs = open("sorted_pe_pe_30_k_5.record", "a")
f.close()
fs.close()

### Read the atomic positon data from the lammps dump file
f = open(sorted_fileName, "r")
ls = []   
p_all = np.zeros((Nf,N,3))                      # Initialize the position data           
px = np.zeros((Nf,N,1))
py = np.zeros((Nf,N,1))
pz = np.zeros((Nf,N,1))
for i in range(Nf):
    ratio = (i+1)/Nf *100
    print("Read Data Progress %s%%" %ratio)
    for k in range(N):                          # 第i帧第k个原子
        line = f.readline().split(",")[2:5]     # 提取234列的数据
        line = [float(l) for l in line]         # 将数据转换为浮点数
        p_all[i, k] = line
        px[i, k] = line[0]
        py[i, k] = line[1]
        pz[i, k] = line[2]

f.close()

### pbc
def mic(dx, dy, dz):
    Lx = 45
    Ly = 45
    Lz = 45
    dx = dx - Lx * np.round(dx / Lx)
    dy = dy - Ly * np.round(dy / Ly)
    dz = dz - Lz * np.round(dz / Lz)

    return dx, dy, dz

for i in range(Nf):
    for k in range(N - 1):
        dx, dy, dz = mic(px[i, k + 1] - px[i, k], 
                         py[i, k + 1] - py[i, k], 
                         pz[i, k + 1] - pz[i, k])
        px[i, k + 1] = px[i, k] + dx
        py[i, k + 1] = py[i, k] + dy
        pz[i, k + 1] = pz[i, k] + dz

### Calculate Rcm from position data
def rcm():
    global Rcmx,Rcmy,Rcmz,Rcm
    Rcmx = np.zeros(Nf)
    Rcmy = np.zeros(Nf)
    Rcmz = np.zeros(Nf)

    for i in range(Nf):
        for k in range (N):
            Rcmx[i] += px[i, k]
            Rcmy[i] += py[i, k]
            Rcmz[i] += pz[i, k] 
        Rcmx[i] = Rcmx[i] / N
        Rcmy[i] = Rcmy[i] / N
        Rcmz[i] = Rcmz[i] / N
    
    Rcm = np.zeros((Nf,4))
    Rcm[:,0] = range(1,int(Nf+1))
    Rcm[:,1] = Rcmx[:]
    Rcm[:,2] = Rcmy[:]
    Rcm[:,3] = Rcmz[:]

    #print (Rcm)

    np.savetxt("rcm.txt",Rcm,fmt="%.5f")

    return Rcmx, Rcmy, Rcmz
rcm ()

### Calculate the mean-square radius of gyration from position data
def rg():
    global Rg2x,Rg2y,Rg2z,Rg2
    Rg2x = np.zeros(Nf)
    Rg2y = np.zeros(Nf)
    Rg2z = np.zeros(Nf)

    for i in range(Nf):
        for k in range(N):
            Rg2x[i] += (px[i,k]-Rcmx[i]) **2
            Rg2y[i] += (py[i,k]-Rcmy[i]) **2
            Rg2z[i] += (pz[i,k]-Rcmz[i]) **2
        Rg2x[i] = Rg2x[i] / N
        Rg2y[i] = Rg2y[i] / N
        Rg2z[i] = Rg2z[i] / N
    
    Rg2 = np.zeros((Nf,4))
    Rg2[:,0] = range(1,int(Nf+1))
    Rg2[:,1] = Rg2x[:]
    Rg2[:,2] = Rg2y[:]
    Rg2[:,3] = Rg2z[:]
   
    np.savetxt("rg2.txt",Rg2,fmt="%.5f")

    return Rg2x,Rg2y,Rg2z
rg ()

### Calculate the mean-square radius of gyration from position data
def re():
    global Re2x,Re2y,Re2z,Re2
    Re2x = np.zeros(Nf)
    Re2y = np.zeros(Nf)
    Re2z = np.zeros(Nf)

    for i in range(Nf):
        for k in range(N-1):
            Re2x[i] += (px[i,k+1]-px[i,k]) **2
            Re2y[i] += (py[i,k+1]-py[i,k]) **2
            Re2z[i] += (pz[i,k+1]-pz[i,k]) **2
    
    Re2 = np.zeros((Nf,4))
    Re2[:,0] = range(1,int(Nf+1))
    Re2[:,1] = Re2x[:]
    Re2[:,2] = Re2y[:]
    Re2[:,3] = Re2z[:]
   
    np.savetxt("re2.txt",Re2,fmt="%.5f")

    return Re2x,Re2y,Re2z
re ()