结构因子（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