【转帖】径向分布函数程序与简单说明 (小木虫)

径向分布函数g(r)代表了球壳内的平均数密度

为离中心分子距离为r，体积为 的球壳内的瞬时分子数。
具体参见李如生，《平衡和非平衡统计力学》科学出版社：1995

 

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=================================