C     ---------------------------------------------------
C          SPC/E ICE Ih COHESIVE ENERGY CALCULATION
C                                  JAN.12,1995 LIJU
C     ---------------------------------------------------

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 576)
      PARAMETER (ID  = 20000)
      COMMON /EWALD/ ALPHA,ERR(ID),CLI(ID),DERR,A(9,9,9)
      COMMON /POS/ X(NPD),Y(NPD),Z(NPD),SX(NPD),SY(NPD),SZ(NPD)
      COMMON /PARM/ NM,NP,ITYPE(NPD),Q(NPD)
      COMMON /ENERGY/ POTE,PLJ12,PLJ6,ELEC
      COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2
      COMMON /CONTROL/ SEED,VOLUME,XLN,YLN,ZLN,CUBE,R,AA,IB
      COMMON /DIRECTION/ MTYPE(192),VECTOR(2,3,2)
      COMMON /CONST1/ SQ2,SQ3,SQ32
      REAL MO,MH

      AVO=6.0221367E+23
      BOLZ=1.380658E-23
      AMU=1.6605E-24

      NM = 192
      NP = 3*NM

      QO    = -0.8476
      QH    = 0.4238
      QOO   = QO*QO
      QOH   = QO*QH
      QHH   = QH*QH
      MO = 15.9994
      MH = 1.00797
      AMASS = MO+2*MH
      BOND = 1.0
      ANGLE = 1.910633236
      HHDIST= 1.632993162
      PI = 3.141592654
      P2 = 2.*PI
      SQ2 = SQRT(2.)
      SQ3 = SQRT(3.)
      SQ32 = SQ3/2.

      SEED  = 0.2573

      VECTOR(1,1,1) = 0.
      VECTOR(1,1,2) = 2*SQ2/3.

      VECTOR(1,2,1) = SQ2/SQ3
      VECTOR(1,2,2) = -SQ2/3.

      VECTOR(1,3,1) = -SQ2/SQ3
      VECTOR(1,3,2) = -SQ2/3.

      VECTOR(2,1,1) = -SQ2/SQ3
      VECTOR(2,1,2) = SQ2/3.

      VECTOR(2,2,1) = SQ2/SQ3
      VECTOR(2,2,2) = SQ2/3.

      VECTOR(2,3,1) = 0.
      VECTOR(2,3,2) = -2*SQ2/3.

      DO 1100 I = 1,NP-2,3
      ITYPE(I) = 0
      ITYPE(I+1)=1
      ITYPE(I+2)=1
      Q(I) = QO
      Q(I+1) = QH
      Q(I+2) = QH
1100  CONTINUE

      OPEN (UNIT=25, FORM='FORMATTED',FILE="data.out")
      LP = 25
      WWW = 5.6

      DO 310 L = 1,200

      DENSITY = 0.900+L*0.001

      TOTALMASS = NM*AMASS
      VOLUME = TOTALMASS*1.66054/DENSITY
      AA = (VOLUME/48./SQ2)**(1./3.)
      R  = SQ3/SQ2/2.*AA
      XLN = 4*AA
      YLN = 3*SQ3*AA
      ZLN = 4*SQ2/SQ3*AA
      CUBE = VOLUME**(1./3.)
     
      ALPHA = WWW/CUBE
      DERR = ALPHA*CUBE/2.*1.05
      DERR = DERR/ID
      C = DERR
      WW = 2./SQRT(PI)

      DO 500 I=1,ID
      ERR(I) = ERFC(C)
      CLI(I) = WW*EXP(-C*C)
      C = C+DERR
500   CONTINUE

      DO 600 I = -4,4
      DO 600 J = -4,4
      DO 600 K = -4,4
      QX = I*P2/XLN
      QY = J*P2/YLN
      QZ = K*P2/ZLN
      Q2 = QX*QX+QY*QY+QZ*QZ
      IF (Q2.EQ.0) GOTO 600
      C1 = EXP(-Q2/ALPHA/ALPHA/4.)
      A(I+5,J+5,K+5) = 4.*PI*C1/Q2/VOLUME
600   CONTINUE

      CALL ICE_Ih

      CALL INTERACT
      POTE=POTE*2.3070795E-18*AVO/NM/1000.
      POTE=POTE+5.22-2./3.*PI*2616.91*NM/VOLUME/(CUBE/2.)**3

      WRITE(LP,278) DENSITY,POTE
 278  FORMAT(2(2X,F13.7))

 310  CONTINUE

      STOP
      END

C     --------------------------------------------------------
      SUBROUTINE RANDOM
      IMPLICIT REAL(A-H,O-Z)
      COMMON /CONTROL/ SEED,VOLUME,XLN,YLN,ZLN,CUBE,R,AA,IB
      DATA K,J,M,RM / 5701,3612,566927,566927.0 /
      IX=INT(SEED*RM)
      KR = J*IX + K
      IRAND = MOD(KR,M)
      SEED = (FLOAT(IRAND) + 0.5) / RM
      RETURN
      END
C     ---------------------------------------------------------


C     --------------------------------------------------------------
      SUBROUTINE INTERACT

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 576)
      PARAMETER (ID  = 20000)
      COMMON /EWALD/ ALPHA,ERR(ID),CLI(ID),DERR,A(9,9,9)
      COMMON /POS/ X(NPD),Y(NPD),Z(NPD),SX(NPD),SY(NPD),SZ(NPD)
      COMMON /PARM/ NM,NP,ITYPE(NPD),Q(NPD)
      COMMON /ENERGY/ POTE,PLJ12,PLJ6,ELEC
      COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2
      COMMON /CONTROL/ SEED,VOLUME,XLN,YLN,ZLN,CUBE,R,AA,IB

      DIMENSION SINN(NPD,9,9,9),COSN(NPD,9,9,9)

      ELEC  = 0.
      PLJ6  = 0.
      PLJ12 = 0.
      W = 2.*ALPHA/SQRT(PI)
      RCUT = CUBE/2.
      FROGX = P2/XLN
      FROGY = P2/YLN
      FROGZ = P2/ZLN
      XLN2 = XLN/2.
      YLN2 = YLN/2.
      ZLN2 = ZLN/2.
      XNL2 = -XLN2
      YNL2 = -YLN2
      ZNL2 = -ZLN2

      DO 805 I = 1,NP
      ELEC = ELEC-Q(I)*Q(I)*W/2.
805   CONTINUE

      DO 1000 I=1,NP
      SINX = SIN(SX(I))
      COSX = COS(SX(I))
      SINY = SIN(SY(I))
      COSY = COS(SY(I))
      SINZ = SIN(SZ(I))
      COSZ = COS(SZ(I))
 
      SINN(I,1,1,1) = SIN(-4.*SX(I)-4.*SY(I)-4.*SZ(I))
      COSN(I,1,1,1) = COS(-4.*SX(I)-4.*SY(I)-4.*SZ(I))

      DO 1100 J=2,9
      SINN(I,J,1,1) = SINN(I,J-1,1,1)*COSX+COSN(I,J-1,1,1)*SINX
      COSN(I,J,1,1) = COSN(I,J-1,1,1)*COSX-SINN(I,J-1,1,1)*SINX
1100  CONTINUE

      DO 1200 J=1,9
      DO 1200 K=2,9
      SINN(I,J,K,1) = SINN(I,J,K-1,1)*COSY+COSN(I,J,K-1,1)*SINY
      COSN(I,J,K,1) = COSN(I,J,K-1,1)*COSY-SINN(I,J,K-1,1)*SINY
1200  CONTINUE

      DO 1300 J=1,9
      DO 1300 K=1,9
      DO 1300 L=2,9
      SINN(I,J,K,L) = SINN(I,J,K,L-1)*COSZ+COSN(I,J,K,L-1)*SINZ
      COSN(I,J,K,L) = COSN(I,J,K,L-1)*COSZ-SINN(I,J,K,L-1)*SINZ
1300  CONTINUE
1000  CONTINUE

      DO 1400 J=1,9
      DO 1400 K=1,9
      DO 1400 L=1,9

      QX = (J-5)*FROGX
      QY = (K-5)*FROGY
      QZ = (L-5)*FROGZ
      Q2 = QX*QX+QY*QY+QZ*QZ
      IF (Q2.EQ.0) GOTO 1400
 
      COSSUM = 0.
      SINSUM = 0.

      DO 1500 I=1,NP
      COSSUM = COSSUM + Q(I)*COSN(I,J,K,L)
      SINSUM = SINSUM + Q(I)*SINN(I,J,K,L)
1500  CONTINUE

      ELEC = ELEC+A(J,K,L)*(COSSUM*COSSUM+SINSUM*SINSUM)/2.
1400  CONTINUE

      DO 522 I=1,NP-1
      DO 5520 J=I+1,NP

      RX = X(I)-X(J)
      RY = Y(I)-Y(J)
      RZ = Z(I)-Z(J)
      IF (RX.LT.XNL2)  RX=RX+XLN
      IF (RX.GT.XLN2)  RX=RX-XLN
      IF (RY.LT.YNL2)  RY=RY+YLN
      IF (RY.GT.YLN2)  RY=RY-YLN
      IF (RZ.LT.ZNL2)  RZ=RZ+ZLN
      IF (RZ.GT.ZLN2)  RZ=RZ-ZLN

      R2 = RX*RX+RY*RY+RZ*RZ
      R1 = SQRT(R2)
      IF (R1.GT.RCUT) GOTO 5520

      R3 = R2*R1
      IJT = ITYPE(I)+ITYPE(J)+1

C     CALCULATE EWARD FACTOR IN REAL SPACE SUMMATION
      PRODUCT = ALPHA*R1
      I0 = PRODUCT/DERR
      XX = I0*DERR
      ARGIN = PRODUCT-XX
      GACTOR = ERR(I0)-CLI(I0)*ARGIN*(1-XX*ARGIN)

C     IF THEY ARE THE SAME MOLECULE
      IF ((J-1)/3.EQ.(I-1)/3) GACTOR = GACTOR-1.
      
      IF ( IJT .EQ. 1 ) THEN
       R6  = R3*R3
       R12 = R6*R6
       E1  = 1895.2940/R12
       E2  = 1.8835408/R6
       PLJ12 = PLJ12+E1
       PLJ6  = PLJ6-E2
       ELEC  = ELEC+GACTOR*QOO/R1

      ELSEIF ( IJT .EQ. 2 ) THEN
       ELEC = ELEC+GACTOR*QOH/R1

      ELSEIF ( IJT .EQ. 3 ) THEN
       ELEC = ELEC+GACTOR*QHH/R1

      ENDIF

5520  CONTINUE
522   CONTINUE

C     CALCULATE TOTAL POTENTIAL ENERGIES
      POTE = PLJ12+PLJ6+ELEC

      RETURN
      END
C     ----------------------------------------------------------------

C     ---------------------------------------------------------------
      BLOCK DATA
      COMMON /DIRECTION/ MTYPE(192),VECTOR(2,3,2)

      DATA MTYPE/ 3,-1,1,-2,-1,3,-3,2,-1,1,-2,-2,1,3,2,-3,2,-3,
     1            -3,3,-2,1,-1,2,
     2            1,2,-2,-2,3,3,-3,2,-1,-2,1,-2,-1,1,-3,2,2,-1,
     3            3,-3,3,-3,1,-1,
     4            -1,-3,3,1,-2,-3,2,-1,1,3,-2,1,3,-2,1,-1,-2,2,
     5            -1,2,-3,3,-3,2,
     6            -1,1,-3,3,-2,2,2,1,3,-1,-2,3,3,-2,-1,-3,1,-1,
     7            -3,-3,2,-2,2,1,
     8            3,-3,2,-2,2,-3,-3,-3,-1,1,1,-2,-2,3,2,1,-1,1,
     9            3,3,-2,2,-1,-1,
     A            -3,-1,3,-2,-1,2,-1,2,-1,-3,1,1,-2,1,-3,2,2,-2,
     B            -2,1,3,-3,3,3,
     C            2,2,-1,2,1,-2,2,-2,3,3,-3,-3,3,-1,3,-3,-1,1,1,
     D            -1,-3,1,-2,-2,
     E            -3,2,-3,2,3,-1,1,-2,2,-1,3,1,-1,-3,-2,2,-3,1,3,
     F            -1,1,-2,3,-2/
      END
C     ----------------------------------------------------------------

C     ----------------------------------------------------------------
      SUBROUTINE ICE_Ih

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 576)
      COMMON /POS/ X(NPD),Y(NPD),Z(NPD),SX(NPD),SY(NPD),SZ(NPD)
      COMMON /PARM/ NM,NP,ITYPE(NPD),Q(NPD)
      COMMON /CONST/ QO,QH,QOO,QOH,QHH,MO,MH,BOND,HHDIST,ANGLE,PI,P2
      COMMON /CONTROL/ SEED,VOLUME,XLN,YLN,ZLN,CUBE,R,AA,IB
      COMMON /DIRECTION/ MTYPE(192),VECTOR(2,3,2)
      COMMON /CONST1/ SQ2,SQ3,SQ32

      BASEX = 0.
      BASEY = 0.
      BASEZ = ZLN
      IP    = 1
      IM    = 1

      DO 100 I=1,8

      IH = MOD(I-1,4)+1
      IF (IH.EQ.1.OR.IH.EQ.4) THEN
      IQ = 1
      ELSE 
      IQ = 2
      ENDIF
      IF (IH.EQ.1.OR.IH.EQ.3) THEN
      PH = 1.
      ELSE 
      PH = -1.
      ENDIF

      DO 200 J=1,6
      DO 200 K=1,4

      IF (IH.EQ.1.OR.IH.EQ.4) THEN
      X(IP) = BASEX+(K-1)*AA+(1-MOD(J,2))*0.5*AA
      ELSE
      X(IP) = BASEX+(K-1)*AA-(1-MOD(J,2))*0.5*AA
      ENDIF

      Y(IP) = BASEY+(J-1)*SQ32*AA
      Z(IP) = BASEZ
      
      MT = MTYPE(IM)
   
      IF (MT.LT.0) THEN
      X(IP+1) = X(IP) 
      Y(IP+1) = Y(IP)
      Z(IP+1) = Z(IP)+PH*BOND
      MT = -MT
      X(IP+2) = X(IP)+VECTOR(IQ,MT,1)*BOND
      Y(IP+2) = Y(IP)+VECTOR(IQ,MT,2)*BOND
      Z(IP+2) = Z(IP)-PH/3.*BOND
     
      ELSE
      MT1 = MOD(MT,3)+1
      X(IP+1) = X(IP)+VECTOR(IQ,MT1,1)*BOND 
      Y(IP+1) = Y(IP)+VECTOR(IQ,MT1,2)*BOND
      Z(IP+1) = Z(IP)-PH/3.*BOND
      MT2 = MOD(MT+1,3)+1
      X(IP+2) = X(IP)+VECTOR(IQ,MT2,1)*BOND
      Y(IP+2) = Y(IP)+VECTOR(IQ,MT2,2)*BOND
      Z(IP+2) = Z(IP)-PH/3.*BOND
      ENDIF

      IP = IP+3
      IM = IM+1
200   CONTINUE

      IF (IH.EQ.1) THEN
      BASEX = BASEX+SQ2/SQ3*R
      BASEY = BASEY-SQ2/3.*R
      BASEZ = BASEZ-1./3.*R
      
      ELSEIF (IH.EQ.3) THEN
      BASEX = BASEX-SQ2/SQ3*R
      BASEY = BASEY+SQ2/3.*R
      BASEZ = BASEZ-1./3.*R

      ELSE
      BASEZ = BASEZ-R
      ENDIF

100   CONTINUE
           
      DO 400 I=1,NP
      IF (X(I).GT.XLN) X(I) = X(I)-XLN     
      IF (X(I).LT.0.)  X(I) = X(I)+XLN     
      IF (Y(I).GT.YLN) Y(I) = Y(I)-YLN     
      IF (Y(I).LT.0.)  Y(I) = Y(I)+YLN     
      IF (Z(I).GT.ZLN) Z(I) = Z(I)-ZLN     
      IF (Z(I).LT.0.)  Z(I) = Z(I)+ZLN     
      SX(I) = X(I)*P2/XLN
      SY(I) = Y(I)*P2/YLN
      SZ(I) = Z(I)*P2/ZLN
400   CONTINUE

      RETURN
      END
C     ---------------------------------------------------------------