C -------------------------------------------------------------------------
C       USE THE FOLLOWING UNITS:                                 
C       ENERGY UNIT = eV  = 1.602E-19 JOULE
C       MASS   UNIT = AMU = 1.66054E-27 KG                      
C       LENGTH UNIT = ANGSTROM                                 
C       TIME   UNIT = 1.018043E-14 SECOND = 10.18043 FS
C       STRESS UNIT = 1.602192 MBAR
C -------------------------------------------------------------------------

      SUBROUTINE INITIALIZE_GENERAL_CONSTANTS
C     ---------------------------------------------------------
#     include "mulcp.fh"
      
      DC_I = DCMPLX(0.D0,1.D0)
      DC_ONE = DCMPLX(1.D0,0.D0)
      DC_ZERO = DCMPLX(0.D0,0.D0)
      
      UENERGY_IN_JOULE = 1.602192E-19
      ULENGTH_IN_METER = 1.E-10
      UMASS_IN_KG = 1.66053E-27
      UTIME_IN_SECOND = 
     A   ULENGTH_IN_METER*SQRT(UMASS_IN_KG/UENERGY_IN_JOULE)
      UFREQ_IN_HZ = 1./UTIME_IN_SECOND/2/PI
      UFREQ_IN_THZ = UFREQ_IN_HZ / 1E12

      AVO   = 6.0221367E+23
      BOLZ  = 1.380658E-23
      HBAR  = 1.054589E-34

C     RCUT = 3.00
      RCUT = 2.56
C     for static problems
      RLIST = RCUT      

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

      
      SUBROUTINE ASSIGN_DYNAMICS
C     ---------------------------------------------------------
#     include "mulcp.fh"
      DO 52 I=1,N     
      PMASS(I) = PM(NTYPE(I))
52    CONTINUE
      RETURN
      END     
C     ---------------------------------------------------------
      
      
      SUBROUTINE UPDATE_VERLET_LIST
C     ---------------------------------------------------------
C     UPDATE THE VERLET LIST. INPUT:SX,SY,SZ,H.
#     include "mulcp.fh"
      
      DO 57 IPT = 1,N
      INDEX(IPT) = 0
      IF (IPT.LE.2) THEN
      RRX(IPT,IPT) = ZERO
      RRY(IPT,IPT) = ZERO
      RRZ(IPT,IPT) = ZERO
      ENDIF
 57   CONTINUE
      
      DO 170 IPT = 1,N
      N1 = NTYPE(IPT)
      SXI = SX(IPT)
      SYI = SY(IPT)
      SZI = SZ(IPT)

      DO 270 JPT = IPT+1,N
      N2 = NTYPE(JPT)
      SXIJ  = SXI - SX(JPT)
      SYIJ  = SYI - SY(JPT)
      SZIJ  = SZI - SZ(JPT)
      IF(SXIJ.GT.HALF)  SXIJ=SXIJ-ONE
      IF(SXIJ.LT.-HALF) SXIJ=SXIJ+ONE
      IF(SYIJ.GT.HALF)  SYIJ=SYIJ-ONE
      IF(SYIJ.LT.-HALF) SYIJ=SYIJ+ONE
      IF(SZIJ.GT.HALF)  SZIJ=SZIJ-ONE
      IF(SZIJ.LT.-HALF) SZIJ=SZIJ+ONE
      RXX = H(1,1)*SXIJ + H(2,1)*SYIJ + H(3,1)*SZIJ
      RYY = H(1,2)*SXIJ + H(2,2)*SYIJ + H(3,2)*SZIJ
      RZZ = H(1,3)*SXIJ + H(2,3)*SYIJ + H(3,3)*SZIJ      
      R1  = SQRT(RXX*RXX+RYY*RYY+RZZ*RZZ)
      
C     STORE THE RELATIVE POSITIONS FOR PHONON DISPERSION CURVE
      IF (IPT.LE.2) THEN
      RRX(IPT,JPT) = RXX
      RRY(IPT,JPT) = RYY
      RRZ(IPT,JPT) = RZZ
      ENDIF

      IF (R1.LT.RLIST) THEN
      INDEX(IPT)=INDEX(IPT)+1
      LIST(IPT,INDEX(IPT))=JPT
      INDEX(JPT)=INDEX(JPT)+1
      LIST(JPT,INDEX(JPT))=IPT
      ENDIF
 270  CONTINUE
 170  CONTINUE

      RRX(2,1) = -RRX(1,2)
      RRY(2,1) = -RRY(1,2)
      RRZ(2,1) = -RRZ(1,2)
      RETURN
      END
C     ---------------------------------------------------------
     
      
      SUBROUTINE ASSIGN_3C_SIC
#     include "mulcp.fh"
C     --------------------------------------------------
C     ASSIGN DIAMOND CUBIC STRUCTURE TO SiC
      DOUBLE PRECISION PX(NM),PY(NM),PZ(NM)
      
      PX(1) = D_ZERO
      PY(1) = D_ZERO
      PZ(1) = D_ZERO
      PX(3) = D_HALF
      PY(3) = D_HALF
      PZ(3) = D_ZERO 
      PX(5) = D_ZERO
      PY(5) = D_HALF
      PZ(5) = D_HALF 
      PX(7) = D_HALF
      PY(7) = D_ZERO
      PZ(7) = D_HALF
      
      IP = 0
      DO 217 I=1,NC(1)
      DO 217 J=1,NC(2)
      DO 217 K=1,NC(3)
      DO 218 L=1,7,2
      SX(IP+L) = (PX(L)+I-7.D0/8.D0)/NC(1)
      SY(IP+L) = (PY(L)+J-7.D0/8.D0)/NC(2)
      SZ(IP+L) = (PZ(L)+K-7.D0/8.D0)/NC(3)
      NTYPE (IP+L) = 1
      SX(IP+L+1) = (PX(L)+I-5.D0/8.D0)/NC(1)
      SY(IP+L+1) = (PY(L)+J-5.D0/8.D0)/NC(2)
      SZ(IP+L+1) = (PZ(L)+K-5.D0/8.D0)/NC(3)
      NTYPE (IP+L+1) = 2
218   CONTINUE
      IP = IP+8
217   CONTINUE  
      
      IF (N.NE.IP) THEN
      PRINT *,'N IS NOT EQUAL TO',IP,'AS DEFINED BY NC(1,2,3).' 
      STOP
      ENDIF
      
      print *,'** n = ',n,sx(1),sy(1),sz(1),sx(2),sy(2),sz(2)
      
      RETURN 
      END
C     --------------------------------------------------
      
      
      SUBROUTINE SAMPLE_3C_SIC_BZ(NSAMPLE,IK,SKX,SKY,SKZ,WEIGHT)
C     ---------------------------------------------------------
C     Monte Carlo or special-point sampling of k space for 
C     the supercell
#     include "mulcp.fh"

      IF (NSAMPLE.EQ.1) THEN 
      SKX = FOURTH*2*PI
      SKY = FOURTH*2*PI
      SKZ = FOURTH*2*PI
      WEIGHT = ONE

      ELSEIF (NSAMPLE.EQ.4) THEN 
      IF (IK.EQ.1) THEN 
      SKX = 3*EIGHTH*2*PI
      SKY =   EIGHTH*2*PI
      SKZ =   EIGHTH*2*PI
      WEIGHT = 3*EIGHTH
      ELSEIF (IK.EQ.2) THEN 
      SKX = 3*EIGHTH*2*PI
      SKY = 3*EIGHTH*2*PI
      SKZ =   EIGHTH*2*PI
      WEIGHT = 3*EIGHTH
      ELSEIF (IK.EQ.3) THEN 
      SKX =   EIGHTH*2*PI
      SKY =   EIGHTH*2*PI
      SKZ =   EIGHTH*2*PI
      WEIGHT =   EIGHTH
      ELSEIF (IK.EQ.4) THEN 
      SKX = 3*EIGHTH*2*PI
      SKY = 3*EIGHTH*2*PI
      SKZ = 3*EIGHTH*2*PI
      WEIGHT =   EIGHTH
      ENDIF

      ELSE
      SKX = (RAN1(ISEED)-HALF)*2*PI
      SKY = (RAN1(ISEED)-HALF)*2*PI
      SKZ = (RAN1(ISEED)-HALF)*2*PI
      WEIGHT = ONE/NSAMPLE
      ENDIF
      RETURN 
      END
C     ---------------------------------------------------------
     
 
      
      SUBROUTINE ASSIGN_4H_SIC
#     include "mulcp.fh"
C     --------------------------------------------------
C     ASSIGN HEXAGONAL 4H STRUCTURE TO SiC
      DOUBLE PRECISION PX(NM),PY(NM),PZ(NM)

C     PLANE A
      PX(1) = D_ZERO
      PY(1) = D_ZERO
      PZ(1) = D_ZERO
      PX(3) = D_HALF
      PY(3) = D_HALF
      PZ(3) = D_ZERO 
c     PLANE B
      PX(5) = PX(1)+D_HALF
      PY(5) = PY(1)+D_SIXTH
      PZ(5) = PZ(1)+D_FOURTH
      PX(7) = PX(3)-D_HALF
      PY(7) = PY(3)+D_SIXTH
      PZ(7) = PZ(3)+D_FOURTH
c     PLANE C
      PX(9) = PX(1)
      PY(9) = PY(1)+D_THIRD
      PZ(9) = PZ(1)+D_HALF
      PX(11) = PX(3)
      PY(11) = PY(3)+D_THIRD
      PZ(11) = PZ(3)+D_HALF
c     PLANE B
      PX(13) = PX(5)
      PY(13) = PY(5)
      PZ(13) = PZ(5)+D_HALF
      PX(15) = PX(7)
      PY(15) = PY(7)
      PZ(15) = PZ(7)+D_HALF
      
      IP = 0
      DO 217 I=1,NC(1)
      DO 217 J=1,NC(2)
      DO 217 K=1,NC(3)
      DO 218 L=1,15,2
      SX(IP+L) = (PX(L)+I-1)/NC(1)
      SY(IP+L) = (PY(L)+J-1)/NC(2)
      SZ(IP+L) = (PZ(L)+K-1)/NC(3)
      NTYPE (IP+L) = 1
      SX(IP+L+1) = (PX(L)+I-1)/NC(1)
      SY(IP+L+1) = (PY(L)+J-1)/NC(2)
      SZ(IP+L+1) = (PZ(L)+K-1+3.D0/16.D0)/NC(3)
      NTYPE (IP+L+1) = 2
218   CONTINUE
      IP = IP+16
217   CONTINUE  
      
      IF (N.NE.IP) THEN
      PRINT *,'N IS NOT EQUAL TO',IP,'   AS DEFINED BY NC(1,2,3).' 
      STOP
      ENDIF

      RETURN 
      END
C     --------------------------------------------------
      
      
      SUBROUTINE SAMPLE_4H_SIC_BZ(NSAMPLE,IK,SKX,SKY,SKZ,WEIGHT)
C     ---------------------------------------------------------
C     Monte Carlo or special-point sampling of k space for 
C     the supercell
#     include "mulcp.fh"

      IF (NSAMPLE.EQ.1) THEN 
      SKX = 0.
      SKY = 0.
      SKZ = 0.
      WEIGHT = ONE

      ELSE
      SKX = (RAN1(ISEED)-HALF)*2*PI
      SKY = (RAN1(ISEED)-HALF)*2*PI
      SKZ = (RAN1(ISEED)-HALF)*2*PI
      WEIGHT = ONE/NSAMPLE
      ENDIF    
      RETURN 
      END     
C     ---------------------------------------------------------
      
      
      SUBROUTINE ASSIGN_HEX3C_SIC
#     include "mulcp.fh"
C     --------------------------------------------------
C     ASSIGN HEXAGONAL 3C STRUCTURE TO SiC
      DOUBLE PRECISION PX(NM),PY(NM),PZ(NM)

C     PLANE A
      PX(1) = D_ZERO
      PY(1) = D_ZERO
      PZ(1) = D_ZERO
      PX(3) = D_HALF      
      PY(3) = D_HALF
      PZ(3) = D_ZERO 
c     PLANE B
      PX(5) = PX(1)+D_HALF
      PY(5) = PY(1)+D_SIXTH
      PZ(5) = PZ(1)+D_THIRD
      PX(7) = PX(3)-D_HALF
      PY(7) = PY(3)+D_SIXTH
      PZ(7) = PZ(3)+D_THIRD
c     PLANE C
      PX(9) = PX(1)
      PY(9) = PY(1)+D_THIRD
      PZ(9) = PZ(1)+2*D_THIRD
      PX(11) = PX(3)
      PY(11) = PY(3)+D_THIRD
      PZ(11) = PZ(3)+2*D_THIRD
      
      IP = 0
      DO 217 I=1,NC(1)
      DO 217 J=1,NC(2)
      DO 217 K=1,NC(3)
      DO 218 L=1,11,2
      SX(IP+L) = (PX(L)+I-1)/NC(1)
      SY(IP+L) = (PY(L)+J-1)/NC(2)
      SZ(IP+L) = (PZ(L)+K-1)/NC(3)
      NTYPE (IP+L) = 1
      SX(IP+L+1) = (PX(L)+I-1)/NC(1)
      SY(IP+L+1) = (PY(L)+J-1)/NC(2)
      SZ(IP+L+1) = (PZ(L)+K-1+D_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

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

      
      
      SUBROUTINE SAMPLE_HEX3C_SIC_BZ(NSAMPLE,IK,SKX,SKY,SKZ,WEIGHT)
C     ---------------------------------------------------------
C     Monte Carlo or special-point sampling of k space for 
C     the supercell
#     include "mulcp.fh"

      IF (NSAMPLE.EQ.1) THEN 
      SKX = 0.
      SKY = 0.
      SKZ = 0.
      WEIGHT = ONE

      ELSE
      SKX = (RAN1(ISEED)-HALF)*2*PI
      SKY = (RAN1(ISEED)-HALF)*2*PI
      SKZ = (RAN1(ISEED)-HALF)*2*PI
      WEIGHT = ONE/NSAMPLE
      ENDIF
      RETURN 
      END     
C     ---------------------------------------------------------
      
      
      SUBROUTINE READ_CONFIG
C     ---------------------------------------------------------
#     include "mulcp.fh"
      DOUBLE PRECISION D_SX(NM),D_SY(NM),D_SZ(NM),D_H(3,3)

      OPEN (UNIT=LP_CONFIG, STATUS='UNKNOWN', FORM='UNFORMATTED',
     a      FILE=NAME_READ)
      READ (LP_CONFIG) N
      IF (N.GT.NM) THEN 
      PRINT*,'Incoming configuration N =',N,
     A       ' more than mulcp can handle.'
      PRINT*,'Increase _NM in mulcp.fh and recompile.'
      STOP
      ENDIF
      READ (LP_CONFIG) (NTYPE(I),I=1,N)
      READ (LP_CONFIG) ((D_H(I,J),J=1,3),I=1,3)
      DO I=1,3
      DO J=1,3
      H(I,J) = D_H(i,j)
      ENDDO
      ENDDO
      READ (LP_CONFIG) (D_SX(I),D_SY(I),D_SZ(I),I=1,N)
      DO I=1,N
      SX(I) = D_SX(I)
      SY(I) = D_SY(I)
      SZ(I) = D_SZ(I)
      ENDDO
c     ONLY NEED THE POSITION PART.
      READ (LP_CONFIG) (D_SX(i),D_SY(i),D_SZ(i),i=1,N)
      CLOSE (LP_CONFIG)
      RETURN
      END
C     ---------------------------------------------------------

      
      FUNCTION RAN1(IDUM)
C     ---------------------------------------------------------
      INTEGER IDUM,IA,IM,IQ,IR,NTAB,NDIV
      REAL RAN1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2E-7,RNMX=1.-EPS)
      INTEGER J,K,IV(NTAB),IY
      SAVE IV,IY
      DATA IV /NTAB*0/, IY /0/
      IF (IDUM.LE.0.OR.IY.EQ.0) THEN
        IDUM=MAX(-IDUM,1)
        DO 11 J=NTAB+8,1,-1
          K=IDUM/IQ
          IDUM=IA*(IDUM-K*IQ)-IR*K
          IF (IDUM.LT.0) IDUM=IDUM+IM
          IF (J.LE.NTAB) IV(J)=IDUM
11      CONTINUE
        IY=IV(1)
      ENDIF
      K=IDUM/IQ
      IDUM=IA*(IDUM-K*IQ)-IR*K
      IF (IDUM.LT.0) IDUM=IDUM+IM
      J=1+IY/NDIV
      IY=IV(J)
      IV(J)=IDUM
      RAN1=MIN(AM*IY,RNMX)
      RETURN
      END
C     ---------------------------------------------------------