[½»Á÷] ¡¾×ªÌû¡¿¾¶Ïò·Ö²¼º¯Êý³ÌÐòÓë¼òµ¥ËµÃ÷ ÒÑÓÐ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 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 ...

Originally posted by zyj8119 at 2010-09-13 21:21:47:

Õâ¸öÊÇCÓïÑԵİɣ¿

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

ÊÇαÂë°É

Originally posted by wangqingbo at 2010-09-16 18:41:17:
¿ÉÒÔÓÃÀ´Ö±½Ó¼ÆËãÂð ¶ÔÓÚÊäÈëµÄ²ÎÊýÊÇ Ê²Ã´ÒªÇó