C     ---------------------------------------------------
C     SPC/E DIAMOND CUBIC ICE COHESIVE ENERGY CALCULATION
C                                        JAN.6,1995 LIJU
C     ---------------------------------------------------

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 648)
      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,CUBE,FROG,IB
      COMMON /DIRECTION/ VECTOR(4,3),MTYPE(192)
      REAL MO,MH

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

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

      SEED  = 0.2573

      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

      RIVER = BOND/(3.**0.5)
      DO 300 I=1,4
      DO 300 J=1,3
      VECTOR(I,J)=VECTOR(I,J)*RIVER
300   CONTINUE

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

      READ *,DENSITY
      TOTALMASS = NM*AMASS
      VOLUME = TOTALMASS*1.66054/DENSITY
      CUBE = VOLUME**(1./3.)
      FROG = P2/CUBE

      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*FROG
      QY = J*FROG
      QZ = K*FROG
      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

      DO 310 I=1,12

      PRINT *,'TYPE',I
      IB = 1+(I-1)*16	
      CALL SEAMLESS
      CALL INTERACT
      POTE=POTE*2.3070795E-18*AVO/NM/1000.
      POTE=POTE+5.22-2./3.*PI*2616.91*NM/VOLUME/(CUBE/2.)**3

      print *,'POTE=',POTE
      PRINT *,' '
310   CONTINUE

      STOP
      END

C     --------------------------------------------------------
      SUBROUTINE RANDOM
      IMPLICIT REAL(A-H,O-Z)
      COMMON /CONTROL/ SEED,VOLUME,CUBE,FROG,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     ---------------------------------------------------------

      BLOCK DATA
      COMMON /DIRECTION/ VECTOR(4,3),MTYPE(192)
      DATA VECTOR/ 1,1,-1,-1,
     Y             1,-1,-1,1,
     Z             1,-1,1,-1/
      DATA MTYPE/ 4,2,4,2,3,1,3,1,-2,-1,-2,-1,-3,-4,-3,-4, 
     2            3,2,4,1,3,2,4,1,-1,-2,-4,-3,-2,-1,-4,-3, 
     3            4,3,4,3,2,1,2,1,-1,-3,-1,-3,-4,-2,-4,-2, 
     4            1,4,2,3,3,2,4,1,-2,-4,-1,-3,-3,-1,-4,-2, 
     5            1,3,4,2,2,4,3,1,-2,-3,-1,-4,-4,-1,-3,-2, 
     6            4,3,1,2,3,4,2,1,-1,-4,-3,-2,-4,-1,-3,-2, 
     7		  1,2,4,3,2,1,3,4,-2,-3,-1,-4,-3,-2,-4,-1, 
     8            4,2,1,3,3,1,2,4,-1,-4,-3,-2,-3,-2,-4,-1, 
     9            3,2,4,1,1,4,2,3,-1,-3,-4,-2,-4,-2,-3,-1, 
     A            1,2,1,2,3,4,3,4,-2,-4,-4,-2,-3,-1,-3,-1, 
     B            1,4,2,3,1,4,2,3,-3,-4,-1,-2,-4,-3,-2,-1, 
     C            1,3,1,3,2,4,2,4,-3,-4,-4,-3,-2,-1,-2,-1/
      END

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

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 648)
      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,CUBE,FROG,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.
      
      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)*FROG
      QY = (K-5)*FROG
      QZ = (L-5)*FROG
      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.-RCUT) RX=RX+CUBE
      IF (RX.GT.RCUT)  RX=RX-CUBE
      IF (RY.LT.-RCUT) RY=RY+CUBE
      IF (RY.GT.RCUT)  RY=RY-CUBE
      IF (RZ.LT.-RCUT) RZ=RZ+CUBE
      IF (RZ.GT.RCUT)  RZ=RZ-CUBE

      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     ----------------------------------------------------------------     
      SUBROUTINE DIAMOND

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 648)
      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,CUBE,FROG,IB
      COMMON /DIRECTION/ VECTOR(4,3),MTYPE(192)

      TA  = TAN(PI-ANGLE)

      SX(1) = 0.
      SY(1) = 0.
      SZ(1) = 0.

      SX(4) =  1./2.
      SY(4) =  1./2.
      SZ(4) =  0.

      SX(7) =  0.
      SY(7) =  1./2.
      SZ(7) =  1./2. 

      SX(10) =  1./2.
      SY(10) =  0.
      SZ(10) =  1./2.

      DO 100 I = 13,22,3
      SX(I) = SX(I-12) + 1./4.
      SY(I) = SY(I-12) + 1./4.
      SZ(I) = SZ(I-12) + 1./4.
100   CONTINUE

      NX = NINT((NM/8)**(1./3.))
      NY = NX
      NZ = NX
      A  = CUBE/NX
      M  = 0
      DO 200 I = 0,NX-1
      DO 200 J = 0,NY-1
      DO 200 K = 0,NZ-1
      DO 300 L = 1,22,3
      X(L+M) = (SX(L)+I+1./8.)*A
      Y(L+M) = (SY(L)+J+1./8.)*A
      Z(L+M) = (SZ(L)+K+1./8.)*A
300   CONTINUE
      M = M + 24
200   CONTINUE

      DO 400 I = 2,NP-1,3
      J=I+1
      
      CALL RANDOM
      XX1 = SEED-0.5
      CALL RANDOM
      YY1 = SEED-0.5
      CALL RANDOM
      ZZ1 = SEED-0.5
      R = SQRT(XX1*XX1+YY1*YY1+ZZ1*ZZ1)
      XX1 = XX1/R
      YY1 = YY1/R
      ZZ1 = ZZ1/R

      CALL RANDOM
      XX2 = SEED-0.5
      CALL RANDOM
      YY2 = SEED-0.5
      CALL RANDOM
      ZZ2 = SEED-0.5
      R1 = SQRT(XX2*XX2+YY2*YY2+ZZ2*ZZ2)
      XX2 = XX2/R1
      YY2 = YY2/R1
      ZZ2 = ZZ2/R1

      X(I) = XX1 + XX2
      Y(I) = YY1 + YY2
      Z(I) = ZZ1 + ZZ2
      R  = SQRT(X(I)*X(I)+Y(I)*Y(I)+Z(I)*Z(I))

      DXX = XX1 - XX2
      DYY = YY1 - YY2
      DZZ = ZZ1 - ZZ2
      R1 = SQRT(DXX*DXX+DYY*DYY+DZZ*DZZ)

      X(J) = -(X(I) + TA*R*DXX/R1)
      Y(J) = -(Y(I) + TA*R*DYY/R1)
      Z(J) = -(Z(I) + TA*R*DZZ/R1)
      R1 = SQRT(X(J)*X(J)+Y(J)*Y(J)+Z(J)*Z(J))

      X(I) = X(I)*BOND/R  + X(I-1)
      Y(I) = Y(I)*BOND/R  + Y(I-1)
      Z(I) = Z(I)*BOND/R  + Z(I-1)
      X(J) = X(J)*BOND/R1 + X(I-1)
      Y(J) = Y(J)*BOND/R1 + Y(I-1)
      Z(J) = Z(J)*BOND/R1 + Z(I-1)
400   CONTINUE

      DO 500 I= 1,NP
      IF(X(I).LT.0)     X(I) = X(I)+CUBE
      IF(X(I).GT.CUBE)  X(I) = X(I)-CUBE
      IF(Y(I).LT.0)     Y(I) = Y(I)+CUBE
      IF(Y(I).GT.CUBE)  Y(I) = Y(I)-CUBE
      IF(Z(I).LT.0)     Z(I) = Z(I)+CUBE
      IF(Z(I).GT.CUBE)  Z(I) = Z(I)-CUBE
      SX(I) = P2*X(I)/CUBE
      SY(I) = P2*Y(I)/CUBE
      SZ(I) = P2*Z(I)/CUBE
500   CONTINUE
      
      RETURN
      END
C     ----------------------------------------------------------------      

C     ----------------------------------------------------------------     
      SUBROUTINE SEAMLESS

      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NPD = 648)
      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,CUBE,FROG,IB
      COMMON /DIRECTION/ VECTOR(4,3),MTYPE(192)

      NX = NINT((NM/8)**(1./3.))
      NY = NX
      NZ = NX
      A  = CUBE/NX

      SX(1) = 0.
      SY(1) = 0.
      SZ(1) = 0.

      SX(7) =  1./2.
      SY(7) =  1./2.
      SZ(7) =  0.

      SX(10) =  0.
      SY(10) =  1./2.
      SZ(10) =  1./2. 

      SX(4) =  1./2.
      SY(4) =  0.
      SZ(4) =  1./2.

      DO 550 I=13,22,3
      SX(I) = SX(I-12)+0.25
      SY(I) = SY(I-12)+0.25
      SZ(I) = SZ(I-12)+0.25
550   CONTINUE

      IO = -2
      DO 500 I = 1,24

      IF (MOD(I,3).EQ.1) THEN
      SX(I) = SX(I)*A
      SY(I) = SY(I)*A
      SZ(I) = SZ(I)*A
      IO = IO+3

      ELSE
      IC = MTYPE(IB)
      IB = IB+1
      IF (IC.LT.0) THEN
      IC = -IC
      SX(I) = SX(IO)-VECTOR(IC,1)
      SY(I) = SY(IO)-VECTOR(IC,2)
      SZ(I) = SZ(IO)-VECTOR(IC,3)
      ELSE
      SX(I) = SX(IO)+VECTOR(IC,1)
      SY(I) = SY(IO)+VECTOR(IC,2)
      SZ(I) = SZ(IO)+VECTOR(IC,3)
      ENDIF
      ENDIF

500   CONTINUE

      M  = 0
      DO 200 I = 0,NX-1
      DO 200 J = 0,NY-1
      DO 200 K = 0,NZ-1
      DO 300 L = 1,24
      X(L+M) = SX(L)+I*A
      Y(L+M) = SY(L)+J*A
      Z(L+M) = SZ(L)+K*A
300   CONTINUE
      M = M + 24
200   CONTINUE

      Q11 = 0.
      Q22 = 0.
      Q33 = 0.
      Q12 = 0.
      Q13 = 0.
      Q23 = 0.
      xx =  0.
      yy =  0.
      zz =  0.

      DO 320  N=1,NP
      R2 = X(N)*X(N)+Y(N)*Y(N)+Z(N)*Z(N)
      Q11 = Q11+(3*X(N)*X(N)-R2)*Q(N)
      Q22 = Q22+(3*Y(N)*Y(N)-R2)*Q(N)
      Q33 = Q33+(3*Z(N)*Z(N)-R2)*Q(N)
      Q12 = Q12+3*X(N)*Y(N)*Q(N)
      Q13 = Q13+3*X(N)*Z(N)*Q(N)
      Q23 = Q23+3*Y(N)*Z(N)*Q(N)
      xx = xx+q(N)*x(n)
      yy = yy+q(N)*y(n)
      zz = zz+q(N)*z(n)
 320  continue

      Q11 = Q11/NM
      Q22 = Q22/NM
      Q33 = Q33/NM
      Q12 = Q12/NM
      Q13 = Q13/NM
      Q23 = Q23/NM
      PRINT *,Q11,Q12,Q13
      PRINT *,Q12,Q22,Q23
      PRINT *,Q13,Q23,Q33
      print *,xx,yy,zz

      DO 700 I= 1,NP
      IF(X(I).LT.0)     X(I) = X(I)+CUBE
      IF(X(I).GT.CUBE)  X(I) = X(I)-CUBE
      IF(Y(I).LT.0)     Y(I) = Y(I)+CUBE
      IF(Y(I).GT.CUBE)  Y(I) = Y(I)-CUBE
      IF(Z(I).LT.0)     Z(I) = Z(I)+CUBE
      IF(Z(I).GT.CUBE)  Z(I) = Z(I)-CUBE
      SX(I) = X(I)*FROG
      SY(I) = Y(I)*FROG
      SZ(I) = Z(I)*FROG
700   CONTINUE
      
      RETURN
      END
C     ---------------------------------------------------