用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 -= regionelif dr < -0.5*region:dr += region
def calRDF():global num_of_A,num_of_Bfor i in range(numatom-1):if type[i] == type_A:num_of_A += 1ix = 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 +=1jx = xu[j]jy = yu[j]jz = zu[j]#print(j)dx = jx-ixdy = jy-iydz = jz-izdx = 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] +=1density = 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)/3normFac = 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]*normFacprint >>f, '%g %g' % (i*dr,g[i])
f.close()

最后

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