rdf分布函数（镜像分布函数rdf计算方法和程序）

rdf分布函数（镜像分布函数rdf计算方法和程序）选中第i个A原子坐标A[i 0:3] 计算所有B原子与A原子的距离dis=|B[: 0:3]-A[i 0:3]|2. 程序实现所以可得g(r)的定义式所以ρBg(r)就是以A为球心 B原子在r处的密度。

声明：本文转自cndaqiang的原创博文，本文只做学术交流，不做商业用途，原文请点击文末阅读原文。

1. 径向分布函数rdf定义

以A原子为中心 半径为r的球内 共计有NAB(r)个B原子

其中ρB是B原子的密度 为常数. 当r→∞应有关系

所以

可得g(r)的定义式

所以ρBg(r)就是以A为球心 B原子在r处的密度。

2. 程序实现

选中第i个A原子坐标A[i 0:3] 计算所有B原子与A原子的距离dis=|B[: 0:3]-A[i 0:3]|

计算在

内的B原子数ΔN(r)=N(r → r dr) 即dis[ (dis >= r) and (dis < r dr) ]的元素数量。

在程序中，可以使用整除向下取余的方式计算。

• 间距dr= dr 区间[minr maxr] 取点数ngrid=int((maxr-minr)/dr)

• r的采样网格点rj=0 ngrid=grid=[0 1 2 3 ... ngrid]*dr minr

• 计算Bj和Ai原子距离dis[j]=|B[j 0:3]-A[i 0:3]|所处区间j=floor(dis[j]-minr)/dr) 这里采用floor向下取整表示(dis >= r) and (dis < r dr)

  N(rj→rj dr)=N(rj→rj dr) 1

• 统计完所有距离Ai小于maxr的B原子后 计算g(rj)=N(rj→rj dr)/ρBΔV(rj)

• 可以外套一层循环以所有的A原子为中心计算一遍 最后处以A原子数平均

• 注意 在成键处的rb 因为键长固定 因此g(rB)随着dr反比变化

3. 代码下载：

https://github.com/cndaqiang/Radial_Distribution_Function_rdf

3.1 依赖环境：

python3Module: numpy matplotlib

3.2 使用示例：

1) 查看帮助

(python37) cndaqiang@mac Desktop$ ./rdf.pyUsage: ./rdf.py [ xxx.xyz | xxx.vasp ] [ none | a | a b c]

2) 计算晶格常数为30Ang的立方晶胞的rdf

(python37) cndaqiang@mac Desktop$ ./rdf.py 303030.xyz 30read from 303030.xyzXYZ: set cell toa: [30. 0. 0.]b: [ 0. 30. 0.]c: [ 0. 0. 30.]nstep ntyp nat 1 ['O' 'H'] [ 884 1768]Super cell [1. 1. 1.]['O-O' 'O-H' 'H-H']Save rdf to 303030.xyz.O-O.averagerdf.datSave rdf to 303030.xyz.O-H.averagerdf.datSave rdf to 303030.xyz.H-H.averagerdf.datSave png to 303030.xyz.averagerdf.png

3) 计算晶格常数a=3.967 b=8.121 c=8.320的长方型原胞rdf

(python37) cndaqiang@mac Desktop$ ./rdf.py 123.xyz 3.967 8.121 8.320read from 123.xyzXYZ: set cell toa: [3.967 0. 0. ]b: [0. 8.121 0. ]c: [0. 0. 8.32]nstep ntyp nat 1 ['O' 'H'] [12 24]Super cell [3. 1. 1.]['O-O' 'O-H' 'H-H']Save rdf to 123.xyz.O-O.averagerdf.datSave rdf to 123.xyz.O-H.averagerdf.datSave rdf to 123.xyz.H-H.averagerdf.datSave png to 123.xyz.averagerdf.png

4) vasp格式

(python37) cndaqiang@mac Desktop$ ./rdf.py 123.vaspread from 123.vaspnstep ntyp nat 1 ['H' 'O'] [24 12]Super cell [3. 1. 1.]['H-H' 'H-O' 'O-O']Save rdf to 123.vasp.H-H.averagerdf.datSave rdf to 123.vasp.H-O.averagerdf.datSave rdf to 123.vasp.O-O.averagerdf.datSave png to 123.vasp.averagerdf.png

4. 结果展示

(python37) cndaqiang@mac Desktop$ ./rdf.py cubic.xyz 12 12 12read from cubic.xyzXYZ: set cell toa: [12. 0. 0.]b: [ 0. 12. 0.]c: [ 0. 0. 12.]nstep ntyp nat 1 ['O' 'H'] [ 61 122]Super cell [2. 2. 2.]['O-O' 'O-H' 'H-H']Save rdf to cubic.xyz.O-O.averagerdf.datSave rdf to cubic.xyz.O-H.averagerdf.datSave rdf to cubic.xyz.H-H.averagerdf.datSave png to cubic.xyz.averagerdf.png

与VMD计算结果对比

参考：

https://zhuanlan.zhihu.com/p/178610319

http://bbs.keinsci.com/thread-15102-1-1.html