結構因子（structure factor）是表征材料對射線的散射能力，包括了原子序數，散射方向和原子位置的影響等等，反映了材料結構的平均信息。實驗中可以獲得結構因子，通過傅立葉變換獲得徑向分布函數（RDF）。分子模擬則相反，先獲得RDF，再通過變換獲得結構因子。
靜態結構因子( Static Structure Factor，SSF) 是一個連接實驗分析與模擬分析的重要參量，對於衍射分析來說，它表征的是材料對射線的散射能力，反映結構的平均信息。通過衍射數據計算結構因子的公式為［15］
統計發現，液態結構因子的第一峰高度在熔點附近都會達到一個同樣的高度，即 S( k) ≈2.8，這個規律叫做 Hansen-Verlet 凝固判據，是可以用來作為熔化和凝固的判據之一［3］。
……………………………………………………………………………………………
哪位能幫忙分析分析我計算出來的結構因子為什么在大q時不趨近於1？我調整過qx、qy、dq的值，可是沒什么作用！
do i = 1 , conf_2
mm=1
do j = 1 , 9
read(1,*)
enddo
do j = 1 , all_n
read(1,*) num_o , type_o , X_o , Y_o , Z_o
if (type_o==1)then !!!!!!!!水分子的坐標
xx(mm) = X_o
yy(mm) = Y_o
zz(mm) = Z_o
mm=mm+1 !!!!!!!水分子的數目
endif
end do
do nx=1,600
do ny=1,600
do nz=1,600
kr=sqrt((nx*qx)**2+(ny*qy)**2+(nz*qz)**2) !!!!!!!q的模
! if ((kr>9).and.(kr<62))then
bin=1+int(kr/0.5)
nhist(bin)=nhist(bin)+1
cossum=0.0
sinsum=0.0
mm1=0
do l = 1 , mm-1
rx=0.1*xx(l)
ry=0.1*yy(l)
rz=0.1*zz(l)
cossum=cossum+cos(nx*qx*rx+ny*qy*ry+nz*qz*rz)
sinsum=sinsum+sin(nx*qx*rx+ny*qy*ry+nz*qz*rz)
mm1=mm1+1
end do
hist(bin) =hist(bin) + (cossum**2+sinsum**2)/real(mm-1)
write(*,*) nhist(bin),hist(bin),mm1,mm,i
!endif
enddo
enddo
enddo
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do k = 1 ,600
if(hist(k)/=0) then
g(k) = real(hist(k))/real(nhist(k))
write(2,*) k*0.5,g(k)
else
g(k)=0
write(2,*) k*0.5,g(k)
endif
end do
……………………………………………………………………………………………
誰能幫忙看看我這段代碼是哪里錯了嗎？為什么我得到的S(k) 和文獻里的差距很大？謝謝
根據這個公式編的一段代碼但是得到的結果都不太對勁的樣子。
因為是正方box所以 倒空間的最小分度gx，gy，gz為：
gx=gy=gz=2*3.14159/Lx （這個大小合適嗎？？？）
gr=gx**2+gy**2+gz**2
sigma為reduced distance
do k=1,100
sinsum = 0.0
cossum = 0.0
param = 0.0
do i = 1, natom-1
do j = i+1, natom
rx=x(i)-x(j)
ry=y(i)-y(j)
rz=z(i)-z(j)
cossum = cossum + cos (k*gx*rx + k*gy*ry + k*gz*rz)
sinsum = sinsum + sin (k*gx*rx + k*gy*ry + k*gz*rz)
enddo
enddo
cossum=cossum/real(natom)
sinsum=sinsum/real(natom)
param = sqrt(cossum**2 + sinsum**2)
! param = param/real(natom)
hist(k)=param+1
enddo
open(200,file=filename2)
do k=1,100
write(200,"(f6.2,x,f11.7)" k*gr*sigma,hist(i)
enddo
close(200)
。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。。
你的問題出在q的處理上，在第二段做球平均以前，波矢q是矢量。所以你的計算中還需加q矢量的三重循環。我下面給你一個使用FFT算S(q)的例子
先給矩陣DT賦值，是你結構每點的form factor值,V是總原子數
CALL REALFT(DT,V,1)
DO I=1,V-1
KR=int(Rij) ！！！距離的bin size,可乘以常數控制曲線點數
SQ(KR)=SQ(KR)+ABS(DT(I+1))/KR/KR ！！結構因子
ENDDO
SUBROUTINE REALFT(DATA1,N,ISIGN)
REAL*8 WR,WI,WPI,WPR,WTEMP,THETA
DIMENSION DATA1(2*N) !USES FOUR1
INTEGER N,ISIGN,I1,I2,I3,I4,I
REAL C2,H2R,H2I,QRS,WIS,H1R,H1I
THETA=6.28318530717959D0/2.0D0/DFLOAT(N)
C1=0.5
IF (ISIGN==1) THEN
C2=-0.5
CALL FOUR1(DATA1,N,+1)
ELSE
C2=0.5
THETA=-THETA
ENDIF
WPR=-2.0D0*DSIN(0.5D0*THETA)**2
WPI=DSIN(THETA)
WR=1.0D0+WPR
WI=WPI
N2P3=2*N+3
DO I=2,N/2+1
I1=2*I-1
I2=I1+1
I3=N2P3-I2
I4=I3+1
WRS=SNGL(WR)
WIS=SNGL(WI)
H1R=C1*(DATA1(I1)+DATA1(I3))
H1I=C1*(DATA1(I2)-DATA1(I4))
H2R=-C2*(DATA1(I2)+DATA1(I4))
H2I=C2*(DATA1(I1)-DATA1(I3))
DATA1(I1)=H1R+WRS*H2R-WIS*H2I
DATA1(I2)=H1I+WRS*H2I+WIS*H2R
DATA1(I3)=H1R-WRS*H2R+WIS*H2I
DATA1(I4)=-H1I+WRS*H2I+WIS*H2R
WTEMP=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WTEMP*WPI+WI
END DO
IF (ISIGN==1) THEN
H1R=DATA1(1)
DATA1(1)=H1R+DATA1(2)
DATA1(2)=H1R-DATA1(2)
ELSE
H1R=DATA1(1)
DATA1(1)=C1*(H1R+DATA1(2))
DATA1(2)=C1*(H1R-DATA1(2))
CALL FOUR1(DATA1,N,-1)
ENDIF
END SUBROUTINE REALFT
SUBROUTINE FOUR1(DATA1,NN,ISIGN)
INTEGER ISIGN,NN
DIMENSION DATA1(2*NN)
INTEGER I,ISTEP,J,M,MMAX,N
REAL TEMPI,TEMPR
DOUBLE PRECISION THETA,WI,WPI,WPR,WR,WTEMP
N=2*NN
J=1
DO I=1,N,2
IF(J>I)THEN
TEMPR=DATA1(J)
TEMPI=DATA1(J+1)
DATA1(J)=DATA1(I)
DATA1(J+1)=DATA1(I+1)
DATA1(I)=TEMPR
DATA1(I+1)=TEMPI
ENDIF
M=N/2
DO WHILE((M>=2).AND.(J>M))
J=J-M
M=M/2
END DO
J=J+M
END DO
MMAX=2
DO WHILE(N>MMAX)
ISTEP=2*MMAX
THETA=6.28318530717959D0/(ISIGN*MMAX)
WPR=-2.D0*DSIN(0.5D0*THETA)**2
WPI=DSIN(THETA)
WR=1.D0
WI=0.D0
DO M=1,MMAX,2
DO I=M,N,ISTEP
J=I+MMAX
TEMPR=SNGL(WR)*DATA1(J)-SNGL(WI)*DATA1(J+1)
TEMPI=SNGL(WR)*DATA1(J+1)+SNGL(WI)*DATA1(J)
DATA1(J)=DATA1(I)-TEMPR
DATA1(J+1)=DATA1(I+1)-TEMPI
DATA1(I)=DATA1(I)+TEMPR
DATA1(I+1)=DATA1(I+1)+TEMPI
END DO
WTEMP=WR
WR=WR*WPR-WI*WPI+WR
WI=WI*WPR+WTEMP*WPI+WI
END DO
MMAX=ISTEP
END DO
END SUBROUTINE FOUR1
---------------------------------------------------------------------------------------------------
非常感謝！
今天又考慮了下代碼，發現少統計了好多對，而且確實也有您說的問題，我只考慮了qx=qy=qz的情況，按照您的意思修改如下：
!!!由公式可知，粒子i和j互換位置的話，就相當於在實部和虛部里考慮了rx=x(i)-x(j)和rx=x(j)-x(i)的正負項。y,z,類似
!!!而cos為偶函數，sin為奇函數，所以cossum相當於加兩倍，而sinsum直接正負抵消了
x方向： do nx = 1,50
y方向： do ny = 1, 50
z方向： do nz = 1, 50
kr = sqrt((nx*qx)**2+(ny*qy)**2+(nz*qz)**2)
bin = INT(kr/0.03) !0.03為dq, 波失q的最小分度
nhist(bin) = nhist(bin)+1 ！這個是對每個q的小窗口中計入q次數的統計
do i = 1, natom-1
do j = i+1, natom
考慮周期性條件后： rx = x(i)-x(j)
ry = y(i)-y(j)
rz = z(i)-z(j)
cossum = cossum+2*cos(nx*qx*rx+ny*qy*ry+nz*qz*rz)
enddo
enddo
cossum = 1+cossum/real(natom)
hist(bin) = hist(bin)+cossum ！hist(bin)為S(k)的統計分布
!!!疑問:這里需要對這些hist(bin)值再求平均嗎？就是用hist(bin)/nhis(bin)
z方向： enddo
y方向： enddo
x方向： enddo