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

欢迎大家来到IT世界,在知识的湖畔探索吧!

声明：本文转自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

欢迎大家来到IT世界,在知识的湖畔探索吧!

3.2 使用示例：

1) 查看帮助

欢迎大家来到IT世界,在知识的湖畔探索吧!(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

欢迎大家来到IT世界,在知识的湖畔探索吧!(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. 结果展示

欢迎大家来到IT世界,在知识的湖畔探索吧!(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

欢迎大家来到IT世界,在知识的湖畔探索吧!