径向分布函数matlab程序,【转帖】径向分布函数程序与简单说明 - 分子模拟 - 小木虫 - 学术 科研 互动社区...

CODE:

SUBROUTINE GR(NSWITCH)

IMPLICIT DOUBLE PRECISION(A-H,O-Z)

PARAMETER(NM=40000,PI=3.141592653589793D0,NHIS=100)

COMMON/LCS/X0(3,-2:2*NM),X(3,-2:2*NM,5),XIN(3,-2:2*NM),

$XX0(3,-2:2*NM),XX(3,-2:2*NM,5),XXIN(3,-2:2*NM)

COMMON/MOLEC/LPBC(3),MOLSP,MOLSA,NBX,NBY,NBZ,NPLA,LPBCSM,NC,NN,MC

COMMON/WALLS/HI(3,3),G(3,3),DH,AREA,VOLUME,SCM(3)

COMMON/PBCS/HALF,PBCX,PBCY,PBCZ

COMMON/GR_VAR/ NGR

DIMENSION H(3,3),GG(0:NHIS),R(0:NHIS)

EQUIVALENCE(X0(1,-2),H(1,1))

C   *****************************************************************

C      如何确定分子数密度：DEN_IDEAL

C      取分子总数作为模拟盒中的数密度，可保证采样分子总数=总分子数?

C====================================================================

C        　N1=MOLSP+1

C      N2=MOLSP+NC

DEN_IDEAL=MOLSP

G11=G(1,1)

G22=G(2,2)

G33=G(3,3)

G12D=G(1,2)+G(2,1)

G13D=G(1,3)+G(3,1)

G23D=G(2,3)+G(3,2)

IF(NSWITCH.EQ.0)THEN

NGR=0

DELR=HALF/NHIS

DO I=1,NHIS

GG(I)=0.D0

R(I)=0.D0

ENDDO

ELSE IF(NSWITCH.EQ.1)THEN

NGR=NGR+1

DO I=1,MOLSP-1

DO J=I+1,MOLSP

C====================================================================

C     USE PBC IN X DIRECTION:  SUITABLE FOR PBCX=1

C                              NOT GREAT PROBLEM FOR PBCX=0

C                              (THIS TIME USUALLY |DELTA X| < HALF)

C====================================================================

XIJ=X0(1,I)-X0(1,J)

IF(XIJ.GT.+HALF)XIJ=XIJ-PBCX

IF(XIJ.LT.-HALF)XIJ=XIJ+PBCX

YIJ=X0(2,I)-X0(2,J)

IF(YIJ.GT.+HALF)YIJ=YIJ-PBCY

IF(YIJ.LT.-HALF)YIJ=YIJ+PBCY

ZIJ=X0(3,I)-X0(3,J)

IF(ZIJ.GT.+HALF)ZIJ=ZIJ-PBCZ

IF(ZIJ.LT.-HALF)ZIJ=ZIJ+PBCZ

RSQ=XIJ*(G11*XIJ+G12D*YIJ+G13D*ZIJ)+

$      YIJ*(G22*YIJ+G23D*ZIJ)+G33*ZIJ*ZIJ

RRR=SQRT(RSQ)

RRR=RRR/H(1,1)

C====================================================================

C      以上用数组G和H的结果与下同

C      RRR=SQRT(XIJ**2+YIJ**2+ZIJ**2)

C      G11=H(1,1)**2

C====================================================================

IF(RRR.LT.HALF)THEN

IG=INT(RRR/DELR)

GG(IG)=GG(IG)+2

ENDIF

ENDDO

ENDDO

ELSE IF(NSWITCH.EQ.2)THEN

DO I=1,NHIS

R(I)=DELR*(I+0.5D0)

ENDDO

DO I=1,NHIS

VB=(4.D0/3.D0)*PI*(((I+1)**3-I**3)*(DELR**3))

GNID=VB*DEN_IDEAL

GG(I)=GG(I)/(NGR*MOLSP*GNID)

ENDDO

OPEN(UNIT=31,FILE="GR.DAT")

DO I=1,NHIS

WRITE(31,*)R(I),GG(I)

ENDDO

CLOSE(31)

ENDIF

RETURN

END