用Python实现径向分布函数（RDF）的计算

什么是RDF？

其实什么是径向分布函数（Radial distribution function）并不需要我去赘述，毕竟百度一搜一大堆。这里给出一个我第一接触RDF时看到的定义，图和公式给的也比较清楚，比较好理解：RDF定义再说明一点，代码所依靠的公式都是来自于剑桥大学的《THE ART OF MOLECULAR DYNAMICS SIMULATION》

思路

首先，因为我需要读lammps的输出的dump文件，所以readxyzufile()所实现的是读入dump文件的功能。然后进入正题，说白了，计算RDF就是计算粒子与粒子之间的距离，并按照距离的远近进行统计。在这里说明一下，为了能够使代码可以计算包含两种粒子的体系，在一些细节方面会有些小的变动。wrap(dr)实现的周期边界条件规范（wraparound）；calRDF()实现的是计算粒子与粒子之间的距离，统计每一层的包含的粒子数。
变量名称说明：
numatom:体系内的粒子总数。
nbins：计算RDF时，球层的总数。
nframe：计算的帧数。
rangemax：计算RDF的范围时0-rangemax。
type_A:粒子A在dump文件中的名称。A粒子为中心粒子。
region：盒子的边长。

代码实现

#calculate RDF
import numpy as np
from numpy import sqrt
def readxyzufile():

    line1 = f.readline()
    line2 = f.readline()
    line3 = f.readline()
    line4 = f.readline()
    natoms = int(line4)
    if natoms != numatom:
        raise ValueError('numatoms is wrong')
    line5 = f.readline()
    line6 = f.readline()
    line7 = f.readline()
    line8 = f.readline()
    line9 = f.readline()

    for d in range(natoms):
        line = f.readline()
        xyz = line.split()
        atomid = int(xyz[0])
        type[atomid-1] = int(xyz[1])
        xu[atomid-1] = float(xyz[2])
        yu[atomid-1] = float(xyz[3])
        zu[atomid-1] = float(xyz[4])
        #print(atomid,xu[atomid-1])
def wrap(dr):
    if dr > 0.5*region:
        dr -= region
    elif dr < -0.5*region:
        dr += region
def calRDF():
    global num_of_A,num_of_B
    for i in range(numatom-1):
        if type[i] == type_A:
            num_of_A += 1
            ix = xu[i]
            iy = yu[i]
            iz = zu[i]
            #print(i)
            for j in range(i+1,numatom):
                if type[j] == type_B:
                    num_of_B +=1
                    jx = xu[j]
                    jy = yu[j]
                    jz = zu[j]
                    #print(j)

                    dx = jx-ix
                    dy = jy-iy
                    dz = jz-iz

                    dx = wrap(dx)
                    dy = wrap(dy)
                    dz = wrap(dz)

                    dij = sqrt(dx*dx+dy*dy+dz*dz)
                    print(i,j,dx,dy,dz,dij)
                    if dij < rangemax:
                        layer = int(dij/dr)
                        #print(layer)
                        count[layer] +=1



density = 1
numatom = 64800
nbins = 50
nframe = 1
rangemax = 20.0
type_A = 1#RDF of atoms A with atom B
type_B = 2
region = 40.16
vol = pow(region,3)
filename = 'unwrap_90_ij400_xyz_2000wan.lammpstrj'
f = open(filename)
dr = rangemax/nbins
num_of_A = 0
num_of_B = 0
type = []
xu = []
yu = []
zu = []
count = []
g = []
for i in range(nbins):
    count.append(0)
    g.append(0)
for i in range(numatom):
    type.append(0)
    xu.append(0)
    yu.append(0)
    zu.append(0)

for i in range(nframe):
    readxyzufile()
    calRDF()
num_of_A = num_of_A/nframe
#density_A = num_of_A/vol
num_of_B = num_of_B/(nframe*num_of_A)
density_B = num_of_B/vol
f.close()
f = open('./rdf_result', 'w')
for i in range(nbins):
    #volume = 4*np.pi*(3*i*i*dr*dr*dr+3*i*dr*dr*dr+dr*dr*dr)/3
    normFac = vol/(4*np.pi*num_of_A*num_of_B*pow((i-0.5)*dr,2)*dr*nframe)
    #normFac = vol/(2*np.pi*numatom*numatom*pow((i-0.5)*dr,2)*dr*nframe)
    #print(density_B)
    print(count[i])
    print(normFac)
    g[i] = count[i]*normFac
    print >>f, '%g %g' % (i*dr,g[i])
f.close()

最后

第一次自己去实现物理量的计算，代码可能还有些问题和不规范的地方，可以指出来大家一起交流讨论。