用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()

最后

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

  • 15
    点赞
  • 38
    收藏
    觉得还不错? 一键收藏
  • 5
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值