[交流] 【转帖】径向分布函数程序与简单说明 已有4人参与

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

Originally posted by ghcacj at 2010-09-13 21:38:22:

是伪码吧

Originally posted by zyj8119 at 2010-09-13 19:42:08:
径向分布函数g(r)代表了球壳内的平均数密度

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

Originally posted by yahoohoo at 2010-09-13 20:39:36:
这样的代码看着不够明了。。。。。。

伪代码：

for (int i = 0; i < TOTN - 1; ++i)
  for (int j = i + 1; j < TOTN; ++j) {
    double dij = sqrt( pow(Pos[0]-Pos[j][0], 2) + pow(Pos[1 ...