首先，因為我需要讀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(
)
徑向分布函式
rdf是徑向分布函式的radical distribution function的縮寫，指的是給定乙個空間，在此空間以乙個物件為中心，去尋找周圍物件的的概率。對於分子模擬的徑向分布函式實則也是求解粒子在週期性邊界盒子的區域密度和全域性密度的比值。區域密度實則就是每乙個球殼的數密度 球殼體積 全域性密...
python 用redis實現分布式鎖
在 redis 中設定值，預設，不存在則建立，存在則修改。引數 1.ex 過期時間 秒 這裡過期時間是3秒，3秒後p，鍵food的值就變成none import redis,time redis client redis.strictredis host localhost port 6379,db...
徑向基網路（RBF）實現函式插值（擬合）
訓練集與測試集 train x 1 0.1 7 train y sin 2.train x test x 0 0.01 8 test y sin 2.test x 引數初始化 rbf.inputsize size train x,1 輸入神經元的個數 rbf.hiddensize 60 隱層神經元的...