C --------------------------------------------------------------------
C     PROGRAM STACKING:
C     TO CREATE A STACKING FAULT IN HEX3C-SIC STRUCTURE
C --------------------------------------------------------------------
      
      PROGRAM STACKING
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NM=8*6**3)
      PARAMETER (PI=3.14159265358979323846D0)
      PARAMETER (ZERO=0.D0,ONE=1.D0,HALF=1.D0/2.D0)
      PARAMETER (THIRD=1.D0/3.D0,FOURTH=1.D0/4.D0,SIXTH=1.D0/6.D0)
      COMMON /FLOAT/ SX(NM),SY(NM),SZ(NM),X(NM),Y(NM),Z(NM),H(3,3)
      COMMON /INT/ N,NTYPE(NM),NC(3)
      DATA LP_CONFIG /25/

C     Number of particles in reference state
      N = 1008
      
C     Cubic lattice constant:
      A = 4.28533
      
C     Bond length:
      B = A * SQRT(3.) /4

C     Triangular lattice constant:
      D = SQRT(B**2-(B/3)**2)*3/2/(SQRT(3.)/2)
      
      H(1,1) = D
      H(1,2) = 0
      H(1,3) = 0
      H(2,1) = 0
      H(2,2) = D*SQRT(3.)
      H(2,3) = 0
      H(3,1) = 0
      H(3,2) = 0
      H(3,3) = D*SQRT(6.)
      
C     Translation of the unit cell in three dimensions, NC(1,2,3):
      NC(1) = 7
      NC(2) = 4
      NC(3) = 3     

      DO 549 I=1,3
      DO 549 J=1,3
 549  H(I,J) = NC(I)*H(I,J)

C     --------------------------------------------------
C     ASSIGN HEXAGONAL 3C STRUCTURE TO SiC

C     PLANE A
      X(1) = ZERO
      Y(1) = ZERO
      Z(1) = ZERO
      X(3) = HALF
      Y(3) = HALF
      Z(3) = ZERO 
C     PLANE B
      X(5) = X(1)+HALF
      Y(5) = Y(1)+SIXTH
      Z(5) = Z(1)+THIRD
      X(7) = X(3)-HALF
      Y(7) = Y(3)+SIXTH
      Z(7) = Z(3)+THIRD
C     PLANE C
      X(9) = X(1)
      Y(9) = Y(1)+THIRD
      Z(9) = Z(1)+2*THIRD
      X(11) = X(3)
      Y(11) = Y(3)+THIRD
      Z(11) = Z(3)+2*THIRD
      
      IP = 0
      DO 217 K=1,NC(3)
      DO 217 I=1,NC(1)
      DO 217 J=1,NC(2)
      DO 218 L=1,11,2
      SX(IP+L) = (X(L)+I-1)/NC(1)
      SY(IP+L) = (Y(L)+J-1)/NC(2)
      SZ(IP+L) = (Z(L)+K-1)/NC(3)
      NTYPE (IP+L) = 1
      SX(IP+L+1) = (X(L)+I-1)/NC(1)
      SY(IP+L+1) = (Y(L)+J-1)/NC(2)
      SZ(IP+L+1) = (Z(L)+K-1+FOURTH)/NC(3)
      NTYPE (IP+L+1) = 2
218   CONTINUE
      IP = IP+12
217   CONTINUE  
      
      IF (N.NE.IP) THEN
      PRINT *,'N IS NOT EQUAL TO',IP,'   AS DEFINED BY NC(1,2,3).' 
      STOP
      ENDIF

C     Perform squashing of the top C layer:
      H(3,3) = H(3,3)/NC(3)/3*(3*NC(3)-1)
      IPP = 0
      DO IP = 1,N
      SZ(IP) = SZ(IP)*3*NC(3)/(3*NC(3)-1)
      IF (SZ(IP).LT.ONE) THEN 
      IPP = IPP+1
      SX(IPP) = SX(IP)
      SY(IPP) = SY(IP)
      SZ(IPP) = SZ(IP)
      ENDIF
      ENDDO
      IF (IPP.NE.N/NC(3)/3*(3*NC(3)-1)) PRINT *, 'SOMETHING IS WRONG.'
      N = IPP
      
      NSI = 0
      DO I = 1,N
      IF (NTYPE(I).EQ.1) NSI = NSI+1
      ENDDO
      PRINT *, 'There are',NSI,' Si and',N-NSI,' C atoms in cell.' 
      
      WRITE (*,'("               |",3f9.4," |")') H(1,1),H(1,2),H(1,3)
      WRITE (*,'(" H(angstrom) = |",3F9.4," |")') H(2,1),H(2,2),H(2,3)
      WRITE (*,'("               |",3F9.4," |")') H(3,1),H(3,2),H(3,3)

C     WRITE ON CONFIGURATION FILE
      OPEN (UNIT=LP_CONFIG, STATUS='UNKNOWN', FORM='UNFORMATTED',
     a      FILE='stack1008')
      WRITE (LP_CONFIG) N
      WRITE (LP_CONFIG) (NTYPE(I),I=1,N)
      WRITE (LP_CONFIG) ((H(I,J),J=1,3),I=1,3)
      WRITE (LP_CONFIG) (SX(I),SY(I),SZ(I),I=1,N)
      WRITE (LP_CONFIG) (0.,0.,0.,I=1,N)
      CLOSE (LP_CONFIG)
      
      STOP
      END