用Python實現徑向分布函式（RDF）的計算

首先，因為我需要讀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
defreadxyzufile()
:    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])
defwrap
(dr)
:if dr >
0.5*region:
dr -= region
elif dr <
-0.5
*region:
dr += region
defcalrdf()
:global num_of_a,num_of_b
for i in
range
(numatom-1)
:iftype
[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 =
1numatom =
64800
nbins =
50nframe =
1rangemax =
20.0
type_a =
1#rdf of atoms a with atom b
type_b =
2region =
40.16
vol =
pow(region,3)
filename =
'unwrap_90_ij400_xyz_2000wan.lammpstrj'
f =open
(filename)
dr = rangemax/nbins
num_of_a =
0num_of_b =
0type=[
]xu =
yu =
zu =
count =
g =[
]for i in
range
(nbins):0
)0)for i in
range
(numatom)
:type0)
0)0)
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(
)

第一次自己去實現物理量的計算，**可能還有些問題和不規範的地方，可以指出來大家一起交流討論。

用Python實現徑向分布函式（RDF）的計算

徑向分布函式

python 用redis實現分布式鎖

徑向基網路（RBF）實現函式插值（擬合）

用Python實現徑向分布函式（RDF）的計算

徑向分布函式

python 用redis實現分布式鎖

徑向基網路（RBF）實現函式插值（擬合）

相關推薦