C--------------------------------------------------------------------
C       PROGRAM LATTICE_DYNAMICS:
C       ORIGINAL VERSION BY LISA PORTER 
C       SENT TO JU LI ON MAY.24,1996. MINOR CHANGES MADE.
C       1. USE TERSOFF POTENTIAL
C       2. CALCULATE DYNAMICAL MATRIX FOR ARBITARY CONFIGURATION
C       3. JACOBI DIAGONALIZATION TO GET THE NORMAL MODES (DOS.OUT)
C       4. FOR PERFECT CRYSTAL: 
C          I) CALCULATE THE DENSITY OF STATES (DOS.OUT)
C         II) OUTPUT A STANDARDIZED MATRIX FILE DM.OUT FOR AUXILLARY
C             MATLAB PROGRAM DISPERSION.M TO GET THE PHONON DISPERSION 
C             CURVE
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 --------------------------------------------------------------------
 
      PROGRAM LATTICE_DYNAMICS
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NM=70,NX=10,PI=3.1415926536 )
      COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM),negno(nm,nm),
     A                NTYPE(NM),LIST(NM,NX),INDEX(NM),PMASS(NM),X(NM),
     A                Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM),
     A                DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE,
     A                ZOLD(NM),LIST2(NM,NX),INDEX2(NM),N,MODE,NEW,
     A                HSUM(3,3),RRX(2,NM),RRY(2,NM),RRZ(2,NM),WRITE_LP,
     A                POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0,
     A                VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE
      COMMON /BLOCK2/ PA(2), PB(2),PL1(2),PL2(2), PBT(2),PN(2),
     A                PC(2),PWR2(2),PD(2), PH(2), PR(2), PS(2), 
     A                C2D2(2), PWR(2),PWR1(2),CHIJ(NX),RLM1(NX),
     A                RLM2(NX),CA(NX),CB(NX),CR(NX),CS(NX),RLIST,
     A                CUBE2,CHI_SI_C,FCUBE2
      COMMON /BLOCK3/ DXX(nm,nm),DXY(nm,nm),DYX(nm,nm),DYY(nm,nm),
     a                DZZ(nm,nm),DXZ(nm,nm),DYZ(nm,nm),DZX(nm,nm),
     a                DZY(nm,nm),HESS(NM*3,NM*3),freq(NM*3),FREQ_K(12),
     a                DM(6,6,NM/2,4),MAX_DM(6,6),S_FREQ_K(NM*6),
     a                amp_si(12), amp_c(12)
      COMMON SEED

      LOGICAL WRITE_LP,UPDATE,NEW
      REAL KX,KY,KZ,DOS(2000),SAMPLE_K(1200000),SOLN(3*NM,3*NM),
     a     samp_si(1200000),dos_si(2000)
      
      DOUBLE PRECISION SX(NM),SY(NM),SZ(NM),S(NM),H(3,3),HH(3,3),DELTA
      DATA LP,LP_CONFIG,LP_DOS,LP_DM  /25,26,27,28/

      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)
      UTIME_IN_PS = UTIME_IN_SECOND/1.E-12 
      UTIME_IN_FS = UTIME_IN_SECOND/1.E-15
      USTRESS_IN_PASCAL = UENERGY_IN_JOULE/ULENGTH_IN_METER**3
      USTRESS_IN_BAR  = USTRESS_IN_PASCAL/1.E5
      USTRESS_IN_MBAR = USTRESS_IN_PASCAL/1.E11
      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
      PM(1) = 28.086
      PM(2) = 12.011 
      ISEED = 32471

C     SET OF PARAMTERS FOR Si. 
      PA(1)   = 1.8308E3
      PB(1)   = 4.7118E2
      PL1(1)  = 2.4799
      PL2(1)  = 1.7322
      PN(1)   = 7.8734E-1
      PBT(1)  = (1.1E-6)**PN(1)
      PC(1)   = (1.0039E5)**2.E0
      PD(1)   = (1.6217E1)**2.E0
      PH(1)   = -5.9825E-1
      PR(1)   = 2.70
      PS(1)   = 3.00
      C2D2(1) = 1.E0+PC(1)/PD(1)
      PWR(1)  = -1.E0/(2.E0*PN(1))
      PWR1(1) = PWR(1)-1.E0
      PWR2(1) = PN(1)-1.E0

C     SET OF PARAMTERS FOR C. 
      PA(2)   = 1.5448E3
      PB(2)   = 3.8963E2
      PL1(2)  = 3.4653
      PL2(2)  = 2.3064
      PN(2)   = 9.9054E-1
      PBT(2)  = (4.1612E-6)**PN(2)
      PC(2)   = (1.9981E4)**2.E0
      PD(2)   = (7.0340)**2.E0
      PH(2)   = -3.9953E-1
      PR(2)   = 1.80
      PS(2)   = 2.10
      C2D2(2) = 1.E0+PC(2)/PD(2)
      PWR(2)  = -1.E0/(2.E0*PN(2))
      PWR1(2) = PWR(2)-1.E0
      PWR2(2) = PN(2)-1.E0

      CHI_SI_C = 1.0086

      RCUT = 3.00
C     for static problems
      RLIST = RCUT

      UPDATE = .TRUE.

      READ *,LAST,CUBE,N,MODE,NPA,NSAMPLE
C     MODE=0, PERFECT CRYSTAL. MODE=1, CLUSTER.
      NC = NINT((N/8.)**(1.E0/3.E0))
      NOLD = N

      IF (LAST.EQ.1) THEN 
      OPEN (UNIT=LP_CONFIG, STATUS='UNKNOWN', FORM='UNFORMATTED',
     a      FILE="config")
       READ (LP_CONFIG) N,H,HH,HH,HH,HH,HH
       READ (LP_CONFIG) SX,SY,SZ
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) S,S,S
       READ (LP_CONFIG) NTYPE,DELTA
      ELSE
C     --------------------------------------------------
C     DIAMOND CUBIC STRUCTURE
       A = 1.E0/FLOAT(NC)
       X(1) = 0.E0
       Y(1) = 0.E0
       Z(1) = 0.E0
       X(3) = 1.E0/2.E0
       Y(3) = 1.E0/2.E0
       Z(3) = 0.E0 
       X(5) = 0.E0
       Y(5) = 1.E0/2.E0
       Z(5) = 1.E0/2.E0 
       X(7) = 1.E0/2.E0
       Y(7) = 0.E0
       Z(7) = 1.E0/2.E0

       IP = 0
       DO 217 I=1,NC
       DO 217 J=1,NC
       DO 217 K=1,NC
       DO 218 L=1,7,2
       SX(IP+L) = (X(L)+I-7.E0/8.E0)*A
       SY(IP+L) = (Y(L)+J-7.E0/8.E0)*A
       SZ(IP+L) = (Z(L)+K-7.E0/8.E0)*A
       NTYPE (IP+L) = 1
       SX(IP+L+1) = (X(L)+I-5.E0/8.E0)*A
       SY(IP+L+1) = (Y(L)+J-5.E0/8.E0)*A
       SZ(IP+L+1) = (Z(L)+K-5.E0/8.E0)*A
       NTYPE (IP+L+1) = 2
 218   CONTINUE
       IP = IP+8
 217   CONTINUE
C     --------------------------------------------------
       ENDIF

       CUBE = CUBE*NC
       CUBE2 = CUBE*0.5
       FCUBE2 = -CUBE2
       DO 52 I=1,N
       X(I) = CUBE*SX(I)
       Y(I) = CUBE*SY(I)
       Z(I) = CUBE*SZ(I)
       PMASS(I) = PM(NTYPE(I))
 52    CONTINUE

      IF (UPDATE) THEN
C    ---------------------------------------------------------  
C    UPDATE THE VERLET LIST

      DO 57 IPT = 1,N
      INDEX(IPT) = 0
      IF (IPT.LE.2) THEN
      RRX(IPT,IPT) = 0.
      RRY(IPT,IPT) = 0.
      RRZ(IPT,IPT) = 0.
      ENDIF
 57   CONTINUE
      
      DO 170 IPT = 1,N
      N1 = NTYPE(IPT)
      XI = X(IPT)
      YI = Y(IPT)
      ZI = Z(IPT)

      DO 270 JPT = IPT+1,N
      N2 = NTYPE(JPT)
      RXX  = X(JPT) - XI
      RYY  = Y(JPT) - YI
      RZZ  = Z(JPT) - ZI
      IF(RXX.GT.CUBE2)  RXX=RXX-CUBE
      IF(RXX.LT.FCUBE2) RXX=RXX+CUBE
      IF(RYY.GT.CUBE2)  RYY=RYY-CUBE
      IF(RYY.LT.FCUBE2) RYY=RYY+CUBE
      IF(RZZ.GT.CUBE2)  RZZ=RZZ-CUBE
      IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE
      R1  = SQRT(RXX*RXX+RYY*RYY+RZZ*RZZ)
      
C     STORE THE RELATIVE POSITION 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)
C     ---------------------------------------------------------
      ENDIF

      call tersoff_dynamical
      print *,' dynamical matrix construction complete.'

      if (mode.eq.0) then
C     ---------------------------------------------------------
C     GENERATE DOS FOR THE CONFIGURATION (K=0)
       call jacobi(hess,3*n,3*NM,freq,soln,nrot)
       do 78 i = 1,3*n
       if (freq(i).lt.0) then 
       print *,freq(i)
       freq(i) = 0.      
       sample_k(i) = 0.
       samp_si(i) = (soln(4,i)+soln(5,i)+soln(6,i))**2
       endif
       freq(i) = sqrt(freq(i))
       sample_k(i) = freq(i)
c      project on Si (111) direction for antisite-pair defect.
       samp_si(i) = (soln(4,i)+soln(5,i)+soln(6,i))**2
 78    continue   
c       call calc_dos(freq,3*n,dos,npa)
       call calc_ldos(sample_k,samp_si,3*n,npa)

      elseif (mode.eq.1) then 
C     ---------------------------------------------------------
C     MONTE CARLO SAMPLING OF K SPACE FOR THE SUPERCELL
       A = 2*PI/CUBE
       NCOUNT = 0
       IP = 0

 119   KX = A*(RAN1(ISEED)-0.5)
       KY = A*(RAN1(ISEED)-0.5)
       KZ = A*(RAN1(ISEED)-0.5)
       IF (NSAMPLE.EQ.1) THEN 
       KX = 0.
       KY = 0.
       KZ = 0.
       ENDIF
       CALL SUPERCELL(KX,KY,KZ) 
       NCOUNT = NCOUNT+1
       DO 29 I=1,N*6
 29    SAMPLE_K(I+IP) = S_FREQ_K(I)
       IP = IP+N*6
       IF (NCOUNT.GE.NSAMPLE) GOTO 209
       GOTO 119
 209   CALL CALC_DOS(SAMPLE_K,6*N*NCOUNT,DOS,NPA)

      else
C     ---------------------------------------------------------
C     MONTE CARLO SAMPLING OF K SPACE FOR PERFECT CRYSTAL
       A = 4*PI*NC/CUBE
       HA = A/2.
       NCOUNT = 0
       IP = 0

 110   KX = A*(RAN1(ISEED)-0.5)
       KY = A*(RAN1(ISEED)-0.5)
       KZ = A*(RAN1(ISEED)-0.5)
       PKX = ABS(KX)
       PKY = ABS(KY)
       PKZ = ABS(KZ)
       
      IF (((PKX+PKY+PKZ).GT.(HA*1.5)).OR.(PKX.GT.HA).OR.
     A   (PKY.GT.HA).OR.(PKZ.GT.HA)) GOTO 110
   
       CALL PHONON(KX,KY,KZ,0) 
       NCOUNT = NCOUNT+1
       DO 27 I=1,12
       SAMPLE_K(I+IP) = FREQ_K(I)
       samp_si(I+IP)  = amp_si(i)
       print *,i,freq_k(i),amp_si(i),1-amp_si(i)
       sum_si = sum_si + amp_si(i)
       sum_c = sum_c + 1 - amp_si(i)
 27    continue
       print *,sum_si, sum_c
       IP = IP+12
       IF (NCOUNT.GE.NSAMPLE) GOTO 200
       GOTO 110
c 200   CALL CALC_DOS(SAMPLE_K,12*NCOUNT,DOS,NPA)
  200   CALL CALC_LDOS(SAMPLE_K,samp_si,12*NCOUNT,NPA)
      ENDIF

      OPEN (UNIT=LP_DM, STATUS='UNKNOWN', FORM='FORMATTED',
     a      FILE="dm.out")
      DO 47 I=1,6
      DO 47 J=1,6
      WRITE (LP_DM,503) MAX_DM(I,J), CUBE/NC, 0, 0
 503  FORMAT (I10,1x,3(F10.4,1X))
      DO 47 K=1,MAX_DM(I,J)
      WRITE (LP_DM,505) 
     A   DM(I,J,K,1),DM(I,J,K,2),DM(I,J,K,3),DM(I,J,K,4)
 47   CONTINUE
 505  FORMAT (4(F10.4,1X))

c     Test the ordered phonon frequency calc.      
c      call phonon (0.1,0.2,0.3,1)
c      print *,freq_k(1)
c      print *,freq_k(2)
c      print *,freq_k(3)
c      print *,freq_k(4)
c      print *,freq_k(5)
c      print *,freq_k(6)
c      print *,' '
      STOP
      END

      SUBROUTINE TERSOFF_DYNAMICAL
      IMPLICIT REAL(A-H,O-Z)
      PARAMETER (NM=70,NX=10,PI=3.1415926536 )
      COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM),negno(nm,nm),
     A                NTYPE(NM),LIST(NM,NX),INDEX(NM),PMASS(NM),X(NM),
     A                Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM),
     A                DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE,
     A                ZOLD(NM),LIST2(NM,NX),INDEX2(NM),N,MODE,NEW,
     A                HSUM(3,3),RRX(2,NM),RRY(2,NM),RRZ(2,NM),WRITE_LP,
     A                POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0,
     A                VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE
      COMMON /BLOCK2/ PA(2), PB(2),PL1(2),PL2(2), PBT(2),PN(2),
     A                PC(2),PWR2(2),PD(2), PH(2), PR(2), PS(2), 
     A                C2D2(2), PWR(2),PWR1(2),CHIJ(NX),RLM1(NX),
     A                RLM2(NX),CA(NX),CB(NX),CR(NX),CS(NX),RLIST,
     A                CUBE2,CHI_SI_C,FCUBE2
      COMMON /BLOCK3/ DXX(nm,nm),DXY(nm,nm),DYX(nm,nm),DYY(nm,nm),
     a                DZZ(nm,nm),DXZ(nm,nm),DYZ(nm,nm),DZX(nm,nm),
     a                DZY(nm,nm),HESS(NM*3,NM*3),freq(NM*3),FREQ_K(12),
     a                DM(6,6,NM/2,4),MAX_DM(6,6),S_FREQ_K(NM*6),
     a                amp_si(12), amp_c(12)

      LOGICAL WRITE_LP

      DIMENSION ZET2(NM,NX,NX),ZET3(NM,NX,NX),FC(NM),FCP(NM),
     A  ZET12(NX,NX),ZET13(NX,NX),ZET22(NX,NX),ZET23(NX,NX),
     A  ZET33(NX,NX),RX(NM,NX),RY(NM,NX),RZ(NM,NX),R(NM,NX),
     A  R2(NM,NX),RXT(NX),RYT(NX),RZT(NX),RJK2(NX,NX),
     A  RJK(NX,NX),RJKX(NX,NX),RJKY(NX,NX),FF11(NM,NX),
     A  RJKZ(NX,NX),PRJKX(NX,NX),PRJKY(NX,NX),PRJKZ(NX,NX),
     A  ZET(NX),SZET1(NX),SZET11(NX),SFTM2X(NX),SFTM2Y(NX),
     A  SFTM2Z(NX),SFTM3X(NX),SFTM3Y(NX),SFTM3Z(NX),PF(NX),
     A  PF2(NX),PF3(NM,NX),FRST(NX),FRST2(NX),FF1(NM,NX),
     A  FF22(NM,NX,NX),FF23(NM,NX,NX),FF13(NM,NX,NX),FF12(NM,NX,NX),
     A  FF3(NM,NX,NX),FF2(NM,NX,NX),FF33(NM,NX,NX),FCPP(NM)

      real lloopxx,lloopxy,nloopxx,nloopxy,
     a     lloop2xx,lloop2xy,nloop2xx,nloop2xy,
     a     lloopyx,nloopyx,lloop2yx,
     a     nloop2yx,lloopyy,nloopyy,lloop2yy,
     a     nloop2yy,lloopzz,lloopxz,nloopzz,
     a     nloopxz,lloop2zz,lloop2xz,nloop2zz,
     a     nloop2xz,lloopyz, lloopzx, lloopzy,
     a     nloopyz, nloopzx, nloopzy, lloop2yz,
     a     lloop2zx, lloop2zy, nloop2yz,nloop2zx,
     a     nloop2zy,
     a     nloop2,lloop3xx,lloop3xy,lloop3xz,
     a     lloop3yx,lloop3yy,lloop3yz,lloop3zx,
     a     lloop3zy,lloop3zz,lloop4xx,lloop4xy,
     a     lloop4xz,lloop4yx,lloop4yy,lloop4yz,
     a     lloop4zx,lloop4zy,lloop4zz,
     a     nloop3xx,nloop3xy,nloop3xz,nloop3yx,
     a     nloop3yy,nloop3yz,nloop3zx,nloop3zy,
     a     nloop3zz,nloop4xx,nloop4xy,nloop4xz,
     a     nloop4yx,nloop4yy,nloop4yz,nloop4zx,
     a     nloop4zy,nloop4zz,
     a     lloop5xx,lloop5yy,lloop5zz,lloop5xy,
     a     lloop5xz,lloop5yx,lloop5yz,lloop5zx,
     a     lloop5zy

      DO 10 IPT=1,N
      FX(IPT) = 0.E0
      FY(IPT) = 0.E0
      FZ(IPT) = 0.E0
      EP(IPT) = 0.E0
      do 10 jpt=1,n
      negno(ipt,jpt) = 0
 10   CONTINUE

      DO 500 IPT = 1,N
      N1 = NTYPE(IPT)
      XI = X(IPT)
      YI = Y(IPT)
      ZI = Z(IPT)
      NEG = 0

      DO 600 J = 1,INDEX(IPT)
      JPT = LIST(IPT,J)
      N2  = NTYPE(JPT)
      RSMAX = PS(N1)*PS(N2)
      RXX  = X(JPT) - XI
      RYY  = Y(JPT) - YI
      RZZ  = Z(JPT) - ZI
      IF(RXX.GT.CUBE2)  RXX=RXX-CUBE
      IF(RXX.LT.FCUBE2) RXX=RXX+CUBE
      IF(RYY.GT.CUBE2)  RYY=RYY-CUBE
      IF(RYY.LT.FCUBE2) RYY=RYY+CUBE
      IF(RZZ.GT.CUBE2)  RZZ=RZZ-CUBE
      IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE
      RR = RXX*RXX+RYY*RYY+RZZ*RZZ
      IF(RR.GT.RSMAX) GOTO 600

      NEG = NEG + 1
      negno(ipt,jpt) = NEG
      LIST2(IPT,NEG) = JPT
      R1  = SQRT(RR)
      RX(IPT,NEG)  = RXX
      RY(IPT,NEG)  = RYY
      RZ(IPT,NEG)  = RZZ
      R2(IPT,NEG)  = RR
      R(IPT,NEG)   = R1
      RXT(NEG) = RXX/R1
      RYT(NEG) = RYY/R1
      RZT(NEG) = RZZ/R1

      CA(NEG) = SQRT(PA(N1)*PA(N2))
      CB(NEG) = SQRT(PB(N1)*PB(N2))
      CR(NEG) = SQRT(PR(N1)*PR(N2))
      CS(NEG) = SQRT(PS(N1)*PS(N2))

      if (n1.eq.1.and.n2.eq.1) then 
      cr(neg) = 1.6
      cs(neg) = 2.9
      elseif (n1.eq.2.and.n2.eq.2) then 
      cr(neg) = 1.8
      cs(neg) = 2.1
      else
      cr(neg) = 2.2
      cs(neg) = 2.5
      endif

      RLM1(NEG) = (PL1(N1)+PL1(N2))/2.E0
      RLM2(NEG) = (PL2(N1)+PL2(N2))/2.E0
      IF(N1.EQ.N2) THEN
      CHIJ(NEG) = 1.E0
      ELSE
      CHIJ(NEG) = 1.0109
      ENDIF

      IF(R1.LT.CR(NEG)) THEN
       FC(NEG)   = 1.E0
       FCP(NEG)  = 0.E0
       FCPP(NEG) = 0.E0
      ELSE
       FC(NEG)=0.5E0*(1.E0+COS(PI*(R1-CR(NEG))
     A    /(CS(NEG)-CR(NEG))))
       FCP(NEG)=-0.5E0*PI*SIN(PI*(R1-CR(NEG))
     A    /(CS(NEG)-CR(NEG)))/(CS(NEG)-CR(NEG))
       FCPP(NEG)=(0.5E0-FC(NEG))*PI*PI/(CS(NEG)-CR(NEG))
     A               /(CS(NEG)-CR(NEG))
      ENDIF
 600  CONTINUE
      INDEX2(IPT) = NEG

      DO 700 J = 1,NEG
      JPT = LIST2(IPT,J)
      N22 = NTYPE(JPT)

      ZET(J) = 0.E0
      SZET1(J)  = 0.E0
      SZET11(J) = 0.E0
      SFTM2X(J) = 0.E0
      SFTM2Y(J) = 0.E0
      SFTM2Z(J) = 0.E0
      SFTM3X(J) = 0.E0
      SFTM3Y(J) = 0.E0
      SFTM3Z(J) = 0.E0

      DO 800 K = 1,NEG
      IF(K.EQ.J) GOTO 800

      RMU = (RX(IPT,J)*RX(IPT,K)+RY(IPT,J)*RY(IPT,K)+
     a       RZ(IPT,J)*RZ(IPT,K))/R(IPT,J)/R(IPT,K)
      RJK2(J,K) = R2(IPT,J)+R2(IPT,K)
     a            -2.E0*RMU*R(IPT,J)*R(IPT,K)
      RJK(J,K)  = SQRT(RJK2(J,K))
      RJKX(J,K) = RX(IPT,J)-RX(IPT,K)
      RJKY(J,K) = RY(IPT,J)-RY(IPT,K)
      RJKZ(J,K) = RZ(IPT,J)-RZ(IPT,K)
      PRJKX(J,K) = RJKX(J,K)/RJK(J,K)
      PRJKY(J,K) = RJKY(J,K)/RJK(J,K)
      PRJKZ(J,K) = RJKZ(J,K)/RJK(J,K)

      RMU1 = RMU-PH(N1)
      RMU2 = RMU1*RMU1
      RMU3 = PD(N1)+RMU2
      GMU  = C2D2(N1)-PC(N1)/RMU3
      GMUP = 2.E0*RMU1*PC(N1)/RMU3**2

      COM1 = FCP(K)*GMU
      COM2 = FC(K)*GMUP
      COM3 = FCPP(K)*GMU
      ZET(J) = ZET(J)+FC(K)*GMU
      
      TIP1 = RMU/R(IPT,J)-1.E0/R(IPT,K)
      TIP2 = RMU/R(IPT,K)-1.E0/R(IPT,J)
      TIP3 = RJK(J,K)/R(IPT,J)/R(IPT,K)
      
      ZET1 = COM2*(-TIP1)
      SZET1(J)  = SZET1(J)+ZET1
      ZET2(IPT,J,K) = COM1+COM2*(-TIP2)
      ZET3(IPT,J,K) = COM2*(-TIP3)

      GMUPP = 2.E0*PC(N1)/RMU3**2-8.E0*PC(N1)*RMU2/RMU3**3
      COM4  = FCP(K)*GMUP
      COM5  = FC(K)*GMUPP

      TIP11 = 2.E0*RMU/R2(IPT,J)-1.E0/R(IPT,J)/R(IPT,K)
      TIP12 = RMU/R(IPT,J)/R(IPT,K)-1.E0/R2(IPT,K)-1.E0/R2(IPT,J)
      TIP13 = RJK(J,K)/R2(IPT,J)/R(IPT,K)
      TIP22 = 2.E0*RMU/R2(IPT,K)-1.E0/R(IPT,J)/R(IPT,K)
      TIP23 = RJK(J,K)/R2(IPT,K)/R(IPT,J)
      TIP33 = -1.E0/R(IPT,J)/R(IPT,K)

      ZET11 = COM5*TIP1*TIP1+COM2*TIP11
      SZET11(J) = SZET11(J) + ZET11
      ZET12(J,K) = COM4*(-TIP1)+COM5*TIP2*TIP1+COM2*TIP12
      ZET13(J,K) = COM5*TIP3*TIP1+COM2*TIP13
      ZET22(J,K) = COM3+2.E0*COM4*(-TIP2)+COM5*TIP2*TIP2+COM2*TIP22
      ZET23(J,K) = COM4*(-TIP3)+COM5*TIP2*TIP3+COM2*TIP23
      ZET33(J,K) = COM5*TIP3*TIP3+COM2*TIP33

      SFTM2X(J) = SFTM2X(J) + ZET2(IPT,J,K)*RXT(K)
      SFTM2Y(J) = SFTM2Y(J) + ZET2(IPT,J,K)*RYT(K)
      SFTM2Z(J) = SFTM2Z(J) + ZET2(IPT,J,K)*RZT(K)
      SFTM3X(J) = SFTM3X(J) + ZET3(IPT,J,K)*PRJKX(J,K)
      SFTM3Y(J) = SFTM3Y(J) + ZET3(IPT,J,K)*PRJKY(J,K)
      SFTM3Z(J) = SFTM3Z(J) + ZET3(IPT,J,K)*PRJKZ(J,K)

 800  CONTINUE
 700  CONTINUE

      DO 704 J=1,NEG
      JPT = LIST2(IPT,J)
      N22 = NTYPE(JPT)

      ELM2  = EXP(-RLM2(J)*R(IPT,J))
      ELM2B = ELM2*CB(J)
      ELM1  = EXP(-RLM1(J)*R(IPT,J))
      ELM1A = ELM1*CA(J)
      FCD1  = FCP(J)-RLM1(J)*FC(J)
      FCD2  = FCP(J)-RLM2(J)*FC(J)
      FCD1P = FCPP(J)-2.E0*RLM1(J)*FCP(J)+RLM1(J)**2.*FC(J)
      FCD2P = FCPP(J)-2.E0*RLM2(J)*FCP(J)+RLM2(J)**2.*FC(J)
 
      IF(ZET(J).NE.0.E0) THEN
       BEZET  = 1.E0+ PBT(N1)*(ZET(J)**PN(N1))
       BIJ    = (BEZET**PWR(N1))*CHIJ(J)
       BIJ1   = 0.5E0*PBT(N1)*(ZET(J)**PWR2(N1))
       BIJP   = BIJ1*(BEZET**PWR1(N1))*CHIJ(J)
       PF(J)  = FC(J)*ELM2B*BIJP
       PF2(J) = FCD2*ELM2B*BIJP
       PF3(IPT,J) = PF(J)*(PWR2(N1)/ZET(J)-
     A      (1.+2.*PN(N1))*BIJ1/BEZET)
      ELSE
       BIJ = 1.E0
       PF(J)  = 0.E0
       PF2(J) = 0.E0
       PF3(IPT,J) = 0.E0
      ENDIF

      FRST(J)  = FCD1*ELM1A - BIJ*FCD2*ELM2B
      FRST2(J) = FCD1P*ELM1A -BIJ*FCD2P*ELM2B

      FF1(IPT,J)  = FRST(J)+PF(J)*SZET1(J)
      FF11(IPT,J) = FRST2(J)+2.E0*PF2(J)*SZET1(J)
     A     +PF3(IPT,J)*SZET1(J)**2+PF(J)*SZET11(J)

C     CALCULATE POTENTIAL ENERGY PER PARTICLE
      VIJR = FC(J)*ELM1A
      VIJA = -FC(J)*BIJ*ELM2B
      VIJ  = VIJR+VIJA
      EP(IPT) = EP(IPT) + VIJ/4.E0
      EP(JPT) = EP(JPT) + VIJ/4.E0

C     THREE BODY FORCES AND THREE BODY STRESSES
      DO 840 K=1,NEG
      KPT = LIST2(IPT,K)
      IF (K.EQ.J) GOTO 840

      FF2(IPT,J,K) = PF(J)*ZET2(IPT,J,K)
      FF3(IPT,J,K) = PF(J)*ZET3(IPT,J,K)

      FF12(IPT,J,K) = PF2(J)*ZET2(IPT,J,K)+PF3(IPT,J)*SZET1(J)
     A                *ZET2(IPT,J,K)+PF(J)*ZET12(J,K)
      FF13(IPT,J,K) = PF2(J)*ZET3(IPT,J,K)+PF(J)*ZET13(J,K)
     A                +PF3(IPT,J)*SZET1(J)*ZET3(IPT,J,K)
      FF22(IPT,J,K) = PF3(IPT,J)*ZET2(IPT,J,K)**2+PF(J)*ZET22(J,K)
      FF23(IPT,J,K) = PF3(IPT,J)*ZET2(IPT,J,K)*ZET3(IPT,J,K)
     A                +PF(J)*ZET23(J,K)
      FF33(IPT,J,K) = PF3(IPT,J)*ZET3(IPT,J,K)**2+PF(J)*ZET33(J,K)

C     CALC FORCE ON PARTICLE K
      FKX = -(FF2(IPT,J,K)*RX(IPT,K)/R(IPT,K)+FF3(IPT,J,K)
     A        *(-PRJKX(J,K)))/2.E0
      FKY = -(FF2(IPT,J,K)*RY(IPT,K)/R(IPT,K)+FF3(IPT,J,K)
     A        *(-PRJKY(J,K)))/2.E0
      FKZ = -(FF2(IPT,J,K)*RZ(IPT,K)/R(IPT,K)+FF3(IPT,J,K)
     A        *(-PRJKZ(J,K)))/2.E0
      FX(KPT) = FX(KPT) + FKX
      FY(KPT) = FY(KPT) + FKY
      FZ(KPT) = FZ(KPT) + FKZ
 840  CONTINUE

C     CALC FORCE ON PARTICLE I
      FIX = ((FRST(J)+PF(J)*SZET1(J))*RXT(J)+PF(J)*SFTM2X(J))/2.E0
      FIY = ((FRST(J)+PF(J)*SZET1(J))*RYT(J)+PF(J)*SFTM2Y(J))/2.E0
      FIZ = ((FRST(J)+PF(J)*SZET1(J))*RZT(J)+PF(J)*SFTM2Z(J))/2.E0
      FX(IPT) = FX(IPT)+FIX
      FY(IPT) = FY(IPT)+FIY
      FZ(IPT) = FZ(IPT)+FIZ

C     CALC FORCE ON PARTICLE J
      FJX=-((FRST(J)+PF(J)*SZET1(J))*RXT(J)+SFTM3X(J)*PF(J))*0.5E0
      FJY=-((FRST(J)+PF(J)*SZET1(J))*RYT(J)+SFTM3Y(J)*PF(J))*0.5E0
      FJZ=-((FRST(J)+PF(J)*SZET1(J))*RZT(J)+SFTM3Z(J)*PF(J))*0.5E0
      FX(JPT) = FX(JPT)+FJX
      FY(JPT) = FY(JPT)+FJY
      FZ(JPT) = FZ(JPT)+FJZ
     
 704  CONTINUE
 500  CONTINUE

      POTE = 0.E0
      DO 400 I =1,N
      POTE  = POTE  + EP(I)
 400  CONTINUE
      print *, pote/n,n
      print *,' '

C     CALCULATE THE DYNAMICAL MATRIX
C     THE NOTATION FOLLOWS THE ANALYTIC DERIVATION BY L.J.PORTER
C     I.E, LLOOP AND LLOOP2 REFER TO THE FIRST SUM OVER L
C     AND THE SECOND SUM OVER L FOUND IN THOSE NOTES.
C
C     THE DERIVATIVE NOTATION IS A "HYBRID EINSTEIN NOTATION".
C     FOR EXAMPLE, DIJDIJVIL = THE SECOND DERIVATIVE OF VIL WITH
C     RESPECT TO RIJ. DIJDIKVIJ = THE SECOND DERIVATIVE OF VIJ WITH
C     RESPECT TO RIJ AND RIK. ETC

C     zero the dynamic matrix
      do 17 i = 1,n
      do 18 j = 1,n
      DXX(i,j) = 0.0
      DXY(i,j) = 0.0
      DXZ(i,j) = 0.0
      DYY(i,j) = 0.0
      DZZ(i,j) = 0.0
      DYX(i,j) = 0.0
      DYZ(i,j) = 0.0
      DZX(i,j) = 0.0
      DZY(i,j) = 0.0
 18   continue
 17   continue

C     MEIJIE DEFINED XIJ,YIJ,ZIJ AS XJ - XI, ETC
C     NOTE ALSO THAT BELOW, I HAVE QUANTITIES LIKE
C     RJNX = RX(I,N) - RX(I,J) BECAUSE AFTER THE
C     LOOPS BELOW, RX(I,J) WILL  = XI - XJ.

      do 21 ipt= 1,n
      do 22 j = 1,INDEX2(ipt)
      RX(ipt,j) = -RX(ipt,j)
      RY(ipt,j) = -RY(ipt,j)
      Rz(ipt,j) = -Rz(ipt,j)
 22   continue
 21   continue

      do 1002 ipt = 1,n
      do 1100 jpt = 1,n
      if (jpt .eq. ipt) goto 1100

      firstxx = 0.0
      firstyy = 0.0
      firstzz = 0.0
      firstxy = 0.0
      firstxz = 0.0
      firstyx = 0.0
      firstyz = 0.0
      firstzx = 0.0
      firstzy = 0.0

      lloopxx = 0.0
      lloopyy = 0.0
      lloopzz = 0.0
      lloopxy = 0.0
      lloopxz = 0.0
      lloopyx = 0.0
      lloopyz = 0.0
      lloopzx = 0.0
      lloopzy = 0.0
     
      lloop2xx = 0.0
      lloop2yy = 0.0
      lloop2zz = 0.0
      lloop2xy = 0.0
      lloop2xz = 0.0
      lloop2yx = 0.0
      lloop2yz = 0.0
      lloop2zx = 0.0
      lloop2zy = 0.0

      lloop5xx = 0.0
      lloop5yy = 0.0
      lloop5zz = 0.0
      lloop5xy = 0.0
      lloop5xz = 0.0
      lloop5yx = 0.0
      lloop5yz = 0.0
      lloop5zx = 0.0
      lloop5zy = 0.0

      negij = negno(ipt,jpt)
      negji = negno(jpt,ipt)
      if (negij.eq.0) goto 1800

      rij  = r(ipt,negij)
      rijx = RX(ipt,negij)
      rijy = RY(ipt,negij)
      rijz = RZ(ipt,negij)
      rij2 = r2(ipt,negij)
      rij3 = rij2*rij

      dijvij = ff1(ipt,negij)
      dijdijvij = ff11(ipt,negij)
      dijvji = ff1(jpt,negji)
      dijdijvji = ff11(jpt,negji)

      firstxx = (rijx*rijx/rij3 - 1.0/rij)*(dijvij + dijvji) +
     a  (rijx*(-rijx)/rij2)*dijdijvij + (rijx*(-rijx)/rij2)*dijdijvji
c
      firstyy = (rijy*rijy/rij3 - 1.0/rij)*(dijvij + dijvji) +
     a  (rijy*(-rijy)/rij2)*dijdijvij + (rijy*(-rijy)/rij2)*dijdijvji
c
      firstxy = (rijx*rijy/rij3)*(dijvij + dijvji) +
     a  (rijx*(-rijy)/rij2)*dijdijvij + (rijx*(-rijy)/rij2)*dijdijvji
c
      firstyx = (rijy*rijx/rij3)*(dijvij + dijvji) +
     a  (rijy*(-rijx)/rij2)*dijdijvij + (rijy*(-rijx)/rij2)*dijdijvji
c
      firstzz = (rijz*rijz/rij3 - 1.0/rij)*(dijvij + dijvji) +
     a  (rijz*(-rijz)/rij2)*dijdijvij + (rijz*(-rijz)/rij2)*dijdijvji
c
      firstxz = (rijx*rijz/rij3)*(dijvij + dijvji) +
     a  (rijx*(-rijz)/rij2)*dijdijvij + (rijx*(-rijz)/rij2)*dijdijvji
c
      firstyz = (rijy*rijz/rij3)*(dijvij + dijvji) +
     a  (rijy*(-rijz)/rij2)*dijdijvij + (rijy*(-rijz)/rij2)*dijdijvji
c
      firstzx = (rijz*rijx/rij3)*(dijvij + dijvji) +
     a  (rijz*(-rijx)/rij2)*dijdijvij + (rijz*(-rijx)/rij2)*dijdijvji
c
      firstzy = (rijz*rijy/rij3)*(dijvij + dijvji) +
     a  (rijz*(-rijy)/rij2)*dijdijvij + (rijz*(-rijy)/rij2)*dijdijvji

      do 1150 l = 1, INDEX2(ipt)
      if (l .eq. negij) goto 1150
      ril = r(ipt,l)
      ril2 = ril*ril
      rilx = RX(ipt,l)
      rily = RY(ipt,l)
      rilz = RZ(ipt,l)
      rmujl = (RX(ipt,negij)*RX(ipt,l) + RY(ipt,negij)*RY(ipt,l) + 
     a     RZ(ipt,negij)*RZ(ipt,l))/(rij*ril)
      rjl2 = rij2 + ril2 -2.0*rmujl*rij*ril
      rjlx = RX(ipt,l)-RX(ipt,negij)
      rjly = RY(ipt,l)-RY(ipt,negij)
      rjlz = RZ(ipt,l)-RZ(ipt,negij)
      rjl = sqrt(rjl2)

      nloopxx = 0.0
      nloopxy = 0.0
      nloopyx = 0.0
      nloopyy = 0.0
      nloopzz = 0.0
      nloopxz = 0.0
      nloopyz = 0.0
      nloopzx = 0.0
      nloopzy = 0.0
     
      do 1200 m = 1,INDEX2(ipt)
      if ((m .eq. negij).or.(m.eq.l)) goto 1200
      rinx = Rx(ipt,m)
      riny = RY(ipt,m)
      rinz = RZ(ipt,m)
      rin  = r(ipt,m)

      dindjlvij = pf3(ipt,negij)*zet2(ipt,negij,m)*zet3(ipt,negij,l)
      dildijvin = pf3(ipt,m)*zet2(ipt,m,negij)*zet2(ipt,m,l)
      dindjlvil = pf3(ipt,l)*zet2(ipt,l,m)*zet3(ipt,l,negij)

      nloopxx = nloopxx + (rilx*(-rijx)/(rij*ril))*dildijvin +
     a     (rinx*rjlx/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopyy = nloopyy + (rily*(-rijy)/(rij*ril))*dildijvin +
     a     (riny*rjly/(rin*rjl))*(dindjlvij + dindjlvil)
c     
      nloopxy = nloopxy + (rilx*(-rijy)/(rij*ril))*dildijvin +
     a     (rinx*rjly/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopyx = nloopyx + (rily*(-rijx)/(rij*ril))*dildijvin +
     a     (riny*rjlx/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopzz = nloopzz + (rilz*(-rijz)/(rij*ril))*dildijvin +
     a     (rinz*rjlz/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopxz = nloopxz + (rilx*(-rijz)/(rij*ril))*dildijvin +
     a     (rinx*rjlz/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopyz = nloopyz + (rily*(-rijz)/(rij*ril))*dildijvin +
     a     (riny*rjlz/(rin*rjl))*(dindjlvij + dindjlvil)
c
      nloopzx = nloopzx + (rilz*(-rijx)/(rij*ril))*dildijvin +
     a     (rinz*rjlx/(rin*rjl))*(dindjlvij + dindjlvil)
c     
      nloopzy = nloopzy + (rilz*(-rijy)/(rij*ril))*dildijvin +
     a     (rinz*rjly/(rin*rjl))*(dindjlvij + dindjlvil)

1200     continue

          dijvil = ff2(ipt,l,negij)
          dijdijvil = ff22(ipt,l,negij)
          dijdjlvij = ff13(ipt,negij,l)
          dijdjlvil = ff23(ipt,l,negij)
          dijdilvil = ff12(ipt,l,negij)
          dijdilvij = ff12(ipt,negij,l)
          djldilvij = ff23(ipt,negij,l)
          djldilvil = ff13(ipt,l,negij)

          lloopxx = lloopxx + (rijx*rijx/rij3 - 1.0/rij)*dijvil +
     a           (rijx*(-rijx)/rij2)*dijdijvil +
     a           (rijx*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilx*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilx*rjlx/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopxx
c
          lloopyy = lloopyy + (rijy*rijy/rij3 - 1.0/rij)*dijvil +
     a           (rijy*(-rijy)/rij2)*dijdijvil +
     a           (rijy*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rily*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rily*rjly/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopyy
c
          lloopxy = lloopxy + (rijx*rijy/rij3)*dijvil +
     a           (rijx*(-rijy)/rij2)*dijdijvil +
     a           (rijx*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilx*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilx*rjly/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopxy
c
          lloopyx = lloopyx + (rijy*rijx/rij3)*dijvil +
     a           (rijy*(-rijx)/rij2)*dijdijvil +
     a           (rijy*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rily*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rily*rjlx/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopyx
c
          lloopzz = lloopzz + (rijz*rijz/rij3 - 1.0/rij)*dijvil +
     a           (rijz*(-rijz)/rij2)*dijdijvil +
     a           (rijz*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilz*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilz*rjlz/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopzz
c
          lloopxz = lloopxz + (rijx*rijz/rij3)*dijvil +
     a           (rijx*(-rijz)/rij2)*dijdijvil +
     a           (rijx*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilx*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilx*rjlz/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopxz
c
          lloopyz = lloopyz + (rijy*rijz/rij3)*dijvil +
     a           (rijy*(-rijz)/rij2)*dijdijvil +
     a           (rijy*rjlz/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rily*(-rijz)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rily*rjlz/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopyz
c
          lloopzx = lloopzx + (rijz*rijx/rij3)*dijvil +
     a           (rijz*(-rijx)/rij2)*dijdijvil +
     a           (rijz*rjlx/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilz*(-rijx)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilz*rjlx/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopzx
c
          lloopzy = lloopzy + (rijz*rijy/rij3)*dijvil +
     a           (rijz*(-rijy)/rij2)*dijdijvil +
     a           (rijz*rjly/(rij*rjl))*(dijdjlvil + dijdjlvij) +
     a           (rilz*(-rijy)/(ril*rij))*(dijdilvil + dijdilvij) +
     a           (rilz*rjly/(ril*rjl))*(djldilvil + djldilvij) +
     a           nloopzy

 1150   continue

        do 1250 l = 1,INDEX2(jpt)
        lpt = LIST2(jpt,l)
        if(lpt .eq. ipt) goto 1250
        rjl = r(jpt,l)
        rjl2 = rjl*rjl
        rjlx = RX(jpt,l)
        rjly = RY(jpt,l)
        rjlz = RZ(jpt,l)
        rmuil = (RX(jpt,negji)*RX(jpt,l) + RY(jpt,negji)*RY(jpt,l) +
     a       RZ(jpt,negji)*RZ(jpt,l))/(rij*rjl)      
        ril2 = rij2 + rjl2 - 2.0*rmuil*rij*rjl
        ril = sqrt(ril2)
        rilx = RX(jpt,l)-RX(jpt,negji)
        rily = RY(jpt,l)-RY(jpt,negji)
        rilz = RZ(jpt,l)-RZ(jpt,negji)

        nloop2xx = 0.0
        nloop2xy = 0.0
        nloop2yx = 0.0
        nloop2yy = 0.0
        nloop2zz = 0.0
        nloop2xz = 0.0
        nloop2yz = 0.0
        nloop2zx = 0.0
        nloop2zy = 0.0
     
        do 1300 m = 1, INDEX2(jpt)
        npt = LIST2(jpt,m)
        if ((npt.eq.ipt) .or.(m.eq.l)) goto 1300
        rjn = r(jpt,m)
        rjn2 = rjn*rjn
        rmuin = (RX(jpt,negji)*RX(jpt,m) + RY(jpt,negji)*RY(jpt,m) +
     a       RZ(jpt,negji)*RZ(jpt,m))/(rij*rjn) 
        rin2 = rij2 + rjn2 - 2.0*rmuin*rij*rjn
        rin = sqrt(rin2)
        rinx = RX(jpt,m)-RX(jpt,negji)
        riny = RY(jpt,m)-RY(jpt,negji)
        rinz = RZ(jpt,m)-RZ(jpt,negji)

        djldijvjn = pf3(jpt,m)*zet2(jpt,m,negji)*zet2(jpt,m,l)
        djldinvji = pf3(jpt,negji)*zet2(jpt,negji,l)*zet3(jpt,negji,m)
        djldinvjn = pf3(jpt,m)*zet2(jpt,m,l)*zet3(jpt,m,negji)

        nloop2xx = nloop2xx + (rijx*rjlx/(rij*rjl))*djldijvjn 
     a       +(rinx*rjlx/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2yy = nloop2yy + (rijy*rjly/(rij*rjl))*djldijvjn 
     a       +(riny*rjly/(rin*rjl))*(djldinvji + djldinvjn)
c     
        nloop2xy = nloop2xy + (rijx*rjly/(rij*rjl))*djldijvjn 
     a       +(rinx*rjly/(rin*rjl))*(djldinvji + djldinvjn)
c     
        nloop2yx = nloop2yx + (rijy*rjlx/(rij*rjl))*djldijvjn 
     a       +(riny*rjlx/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2zz = nloop2zz + (rijz*rjlz/(rij*rjl))*djldijvjn 
     a       +(rinz*rjlz/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2xz = nloop2xz + (rijx*rjlz/(rij*rjl))*djldijvjn 
     a       +(rinx*rjlz/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2yz = nloop2yz + (rijy*rjlz/(rij*rjl))*djldijvjn 
     a       +(riny*rjlz/(rin*rjl))*(djldinvji + djldinvjn)
c
        nloop2zx = nloop2zx + (rijz*rjlx/(rij*rjl))*djldijvjn 
     a       +(rinz*rjlx/(rin*rjl))*(djldinvji + djldinvjn)
c     
        nloop2zy = nloop2zy + (rijz*rjly/(rij*rjl))*djldijvjn 
     a       +(rinz*rjly/(rin*rjl))*(djldinvji + djldinvjn)

 1300   continue

        dijvjl = ff2(jpt,l,negji)
        djldijvji = ff12(jpt,negji,l)
        dijdijvjl = ff22(jpt,l,negji)
        dildjlvjl = ff13(jpt,l,negji)
        dijdjlvjl = ff12(jpt,l,negji)
        dildjlvji = ff23(jpt,negji,l)
        dijdilvji = ff13(jpt,negji,l)
        dijdilvjl = ff23(jpt,l,negji)
c
        lloop2xx = lloop2xx +(rijx*rijx/rij3 - 1.0/rij)*dijvjl +
     a       (rijx*rjlx/(rij*rjl))*djldijvji +
     a       (rijx*(-rijx)/rij2)*dijdijvjl +
     a       (rilx*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijx*rjlx/(rij*rjl))*dijdjlvjl +
     a       (rilx*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2xx
c
        lloop2yy = lloop2yy +(rijy*rijy/rij3 - 1.0/rij)*dijvjl +
     a       (rijy*rjly/(rij*rjl))*djldijvji +
     a       (rijy*(-rijy)/rij2)*dijdijvjl +
     a       (rily*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijy*rjly/(rij*rjl))*dijdjlvjl +
     a       (rily*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2yy
c
        lloop2xy = lloop2xy +(rijx*rijy/rij3)*dijvjl +
     a       (rijx*rjly/(rij*rjl))*djldijvji +
     a       (rijx*(-rijy)/rij2)*dijdijvjl +
     a       (rilx*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijx*rjly/(rij*rjl))*dijdjlvjl +
     a       (rilx*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2xy
c     
        lloop2yx = lloop2yx +(rijy*rijx/rij3)*dijvjl +
     a       (rijy*rjlx/(rij*rjl))*djldijvji +
     a       (rijy*(-rijx)/rij2)*dijdijvjl +
     a       (rily*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijy*rjlx/(rij*rjl))*dijdjlvjl +
     a       (rily*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2yx
c
        lloop2zz = lloop2zz +(rijz*rijz/rij3 - 1.0/rij)*dijvjl +
     a       (rijz*rjlz/(rij*rjl))*djldijvji +
     a       (rijz*(-rijz)/rij2)*dijdijvjl +
     a       (rilz*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijz*rjlz/(rij*rjl))*dijdjlvjl +
     a       (rilz*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2zz
c
        lloop2xz = lloop2xz +(rijx*rijz/rij3)*dijvjl +
     a       (rijx*rjlz/(rij*rjl))*djldijvji +
     a       (rijx*(-rijz)/rij2)*dijdijvjl +
     a       (rilx*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijx*rjlz/(rij*rjl))*dijdjlvjl +
     a       (rilx*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2xz
c
        lloop2yz = lloop2yz +(rijy*rijz/rij3)*dijvjl +
     a       (rijy*rjlz/(rij*rjl))*djldijvji +
     a       (rijy*(-rijz)/rij2)*dijdijvjl +
     a       (rily*(-rijz)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijy*rjlz/(rij*rjl))*dijdjlvjl +
     a       (rily*rjlz/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2yz
c
        lloop2zx = lloop2zx +(rijz*rijx/rij3)*dijvjl +
     a       (rijz*rjlx/(rij*rjl))*djldijvji +
     a       (rijz*(-rijx)/rij2)*dijdijvjl +
     a       (rilz*(-rijx)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijz*rjlx/(rij*rjl))*dijdjlvjl +
     a       (rilz*rjlx/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2zx
c
        lloop2zy = lloop2zy +(rijz*rijy/rij3)*dijvjl +
     a       (rijz*rjly/(rij*rjl))*djldijvji +
     a       (rijz*(-rijy)/rij2)*dijdijvjl +
     a       (rilz*(-rijy)/(ril*rij))*(dijdilvji + dijdilvjl)+
     a       (rijz*rjly/(rij*rjl))*dijdjlvjl +
     a       (rilz*rjly/(ril*rjl))*(dildjlvji + dildjlvjl)+
     a       nloop2zy

 1250   continue

 1800  continue
       do 1700 l = 1,INDEX2(ipt)
       lpt = LIST2(ipt,l)
       if (lpt .eq. jpt) goto 1700
       if (negno(jpt,lpt) .eq. 0) goto 1700
       negli = negno(lpt,ipt)
       neglj = negno(lpt,jpt)

       ril = r(ipt,l)
       rilx = RX(ipt,l)
       rily = RY(ipt,l)
       rilz = RZ(ipt,l)
       ril2 = ril*ril
       rlj = r(lpt,neglj)
       rlj2 = rlj*rlj
       rljx = RX(lpt,neglj)
       rljy = RY(lpt,neglj)
       rljz = RZ(lpt,neglj)

c     The following is a calculation of the distance rij if j is not
c     a nearest neighbor of i. From above, l must be a neighbor of both
c     i and j, so we can get rij through l.

       if (negij .eq. 0) then
       rmuij=(RX(lpt,negli)*RX(lpt,neglj)+RY(lpt,negli)*RY(lpt,neglj)
     a        +RZ(lpt,negli)*RZ(lpt,neglj))/(ril*rlj)
        rij2 = ril2+rlj2-2.0*rmuij*ril*rlj
        rij = sqrt(rij2)
        rij3 = rij2*rij
        rijx = RX(lpt,neglj)-RX(lpt,negli) 
        rijy = RY(lpt,neglj)-RY(lpt,negli)
        rijz = RZ(lpt,neglj)-RZ(lpt,negli)
       endif 

       nloop3xx = 0.0
       nloop3xy = 0.0
       nloop3xz = 0.0
       nloop3yx = 0.0
       nloop3yy = 0.0
       nloop3yz = 0.0
       nloop3zx = 0.0
       nloop3zy = 0.0
       nloop3zz = 0.0
     
       do 1750 m = 1,INDEX2(lpt)
       npt = LIST2(lpt,m)
       if((npt.eq.ipt).or.(npt.eq.jpt)) goto 1750
       rln = r(lpt,m)
       rln2 = rln*rln
       rmuin = (RX(lpt,negli)*RX(lpt,m) + RY(lpt,negli)*RY(lpt,m) +
     a      RZ(lpt,negli)*RZ(lpt,m))/(ril*rln)
       rin2 = ril2 + rln2 - 2.0*rmuin*ril*rln
       rin = sqrt(rin2)
       rinx =  RX(lpt,m)-RX(lpt,negli)
       riny =  RY(lpt,m)-RY(lpt,negli)
       rinz =  RZ(lpt,m)-RZ(lpt,negli)
       rmujn = (RX(lpt,neglj)*RX(lpt,m) + RY(lpt,neglj)*RY(lpt,m) +
     a      RZ(lpt,neglj)*RZ(lpt,m))/(rlj*rln)           
       rjn2 = rln2 + rlj2 - 2.0*rmujn*rln*rlj
       rjn = sqrt(rjn2)
       rjnx = RX(lpt,m)-RX(lpt,neglj)
       rjny = RY(lpt,m)-RY(lpt,neglj)
       rjnz = RZ(lpt,m)-RZ(lpt,neglj)

       dindjlvli = pf3(lpt,negli)*zet2(lpt,negli,neglj)
     a      *zet3(lpt,negli,m)
       dindjlvln = pf3(lpt,m)*zet2(lpt,m,neglj)*zet3(lpt,m,negli)
       dildjnvlj= pf3(lpt,neglj)*zet2(lpt,neglj,negli)
     a      *zet3(lpt,neglj,m)
       dildjnvln = pf3(lpt,m)*zet2(lpt,m,negli)*zet3(lpt,m,neglj)
       dindijvli = pf3(lpt,negli)*zet3(lpt,negli,m)
     a      *zet3(lpt,negli,neglj)
       dildjlvln = pf3(lpt,m)*zet2(lpt,m,negli)*zet2(lpt,m,neglj)
       dijdjnvlj = pf3(lpt,neglj)*zet3(lpt,neglj,negli)
     a      *zet3(lpt,neglj,m)
       dindjnvln = pf3(lpt,m)*zet3(lpt,m,negli)*zet3(lpt,m,neglj)
     
       nloop3xx = (rinx*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a      (rilx*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a      (rinx*(-rijx)/(rin*rij))*dindijvli +
     a      (rilx*(-rljx)/(ril*rlj))*dildjlvln +
     a      (rijx*rjnx/(rij*rjn))*dijdjnvlj +
     a      (rinx*rjnx/(rin*rjn))*dindjnvln + nloop3xx
c
        nloop3yy = (riny*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a      (rily*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a      (riny*(-rijy)/(rin*rij))*dindijvli +
     a      (rily*(-rljy)/(ril*rlj))*dildjlvln +
     a      (rijy*rjny/(rij*rjn))*dijdjnvlj +
     a      (riny*rjny/(rin*rjn))*dindjnvln + nloop3yy
c
        nloop3zz = (rinz*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilz*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinz*(-rijz)/(rin*rij))*dindijvli +
     a       (rilz*(-rljz)/(ril*rlj))*dildjlvln +
     a       (rijz*rjnz/(rij*rjn))*dijdjnvlj +
     a       (rinz*rjnz/(rin*rjn))*dindjnvln + nloop3zz
c     
        nloop3xy = (rinx*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilx*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinx*(-rijy)/(rin*rij))*dindijvli +
     a       (rilx*(-rljy)/(ril*rlj))*dildjlvln +
     a       (rijx*rjny/(rij*rjn))*dijdjnvlj +
     a       (rinx*rjny/(rin*rjn))*dindjnvln + nloop3xy
c
        nloop3xz = (rinx*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilx*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinx*(-rijz)/(rin*rij))*dindijvli +
     a       (rilx*(-rljz)/(ril*rlj))*dildjlvln +
     a       (rijx*rjnz/(rij*rjn))*dijdjnvlj +
     a       (rinx*rjnz/(rin*rjn))*dindjnvln + nloop3xz
c
        nloop3yx = (riny*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rily*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (riny*(-rijx)/(rin*rij))*dindijvli +
     a       (rily*(-rljx)/(ril*rlj))*dildjlvln +
     a       (rijy*rjnx/(rij*rjn))*dijdjnvlj +
     a       (riny*rjnx/(rin*rjn))*dindjnvln + nloop3yx
c     
        nloop3yz = (riny*(-rljz)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rily*rjnz/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (riny*(-rijz)/(rin*rij))*dindijvli +
     a       (rily*(-rljz)/(ril*rlj))*dildjlvln +
     a       (rijy*rjnz/(rij*rjn))*dijdjnvlj +
     a       (riny*rjnz/(rin*rjn))*dindjnvln + nloop3yz
c     
        nloop3zx = (rinz*(-rljx)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilz*rjnx/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinz*(-rijx)/(rin*rij))*dindijvli +
     a       (rilz*(-rljx)/(ril*rlj))*dildjlvln +
     a       (rijz*rjnx/(rij*rjn))*dijdjnvlj +
     a       (rinz*rjnx/(rin*rjn))*dindjnvln + nloop3zx
c
        nloop3zy = (rinz*(-rljy)/(rin*rlj))*(dindjlvli+dindjlvln) +
     a       (rilz*rjny/(ril*rjn))*(dildjnvlj + dildjnvln) +
     a       (rinz*(-rijy)/(rin*rij))*dindijvli +
     a       (rilz*(-rljy)/(ril*rlj))*dildjlvln +
     a       (rijz*rjny/(rij*rjn))*dijdjnvlj +
     a       (rinz*rjny/(rin*rjn))*dindjnvln + nloop3zy
     
 1750    continue

         dijvli = ff3(lpt,negli,neglj)
         dijvlj = ff3(lpt,neglj,negli)
         dildjlvli = ff12(lpt,negli,neglj)
         dildjlvlj = ff12(lpt,neglj,negli)
         dildijvli = ff13(lpt,negli,neglj)
         dildijvlj = ff23(lpt,neglj,negli)
         dijdjlvli = ff23(lpt,negli,neglj)
         dijdjlvlj = ff13(lpt,neglj,negli)
         dijdijvli = ff33(lpt,negli,neglj)
         dijdijvlj = ff33(lpt,neglj,negli)

         lloop5xx = (rijx*rijx/rij3 - 1.0/rij)*(dijvli + dijvlj) +
     a             (rilx*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilx*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijx*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijx*(-rijx)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3xx + lloop5xx
c
         lloop5yy = (rijy*rijy/rij3 - 1.0/rij)*(dijvli + dijvlj) +
     a             (rily*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rily*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijy*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijy*(-rijy)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3yy + lloop5yy
c
         lloop5zz = (rijz*rijz/rij3 - 1.0/rij)*(dijvli + dijvlj) +
     a             (rilz*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilz*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijz*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijz*(-rijz)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3zz + lloop5zz
c
         lloop5xy = (rijx*rijy/rij3)*(dijvli + dijvlj) +
     a             (rilx*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilx*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijx*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijx*(-rijy)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3xy + lloop5xy
c
         lloop5xz = (rijx*rijz/rij3)*(dijvli + dijvlj) +
     a             (rilx*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilx*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijx*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijx*(-rijz)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3xz + lloop5xz
c
         lloop5yx = (rijy*rijx/rij3)*(dijvli + dijvlj) +
     a             (rily*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rily*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijy*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijy*(-rijx)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3yx + lloop5yx
c
         lloop5yz = (rijy*rijz/rij3)*(dijvli + dijvlj) +
     a             (rily*(-rljz)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rily*(-rijz)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijy*(-rljz)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijy*(-rijz)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3yz + lloop5yz
c
         lloop5zx = (rijz*rijx/rij3)*(dijvli + dijvlj) +
     a             (rilz*(-rljx)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilz*(-rijx)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijz*(-rljx)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijz*(-rijx)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3zx + lloop5zx
c
         lloop5zy = (rijz*rijy/rij3)*(dijvli + dijvlj) +
     a             (rilz*(-rljy)/(ril*rlj))*(dildjlvli + dildjlvlj) +
     a             (rilz*(-rijy)/(ril*rij))*(dildijvli + dildijvlj) +
     a             (rijz*(-rljy)/(rij*rlj))*(dijdjlvli + dijdjlvlj) +
     a             (rijz*(-rijy)/rij2)*(dijdijvli + dijdijvlj) +
     a             nloop3zy + lloop5zy
c
 1700  continue

        DXX(ipt,jpt) = 0.50*(firstxx + lloopxx + lloop2xx + lloop5xx)
        DXY(ipt,jpt) = 0.50*(firstxy + lloopxy + lloop2xy + lloop5xy)
        DYX(ipt,jpt) = 0.50*(firstyx + lloopyx + lloop2yx + lloop5yx)
        DYY(ipt,jpt) = 0.50*(firstyy + lloopyy + lloop2yy + lloop5yy)
        DZZ(ipt,jpt) = 0.50*(firstzz + lloopzz + lloop2zz + lloop5zz)
        DXZ(ipt,jpt) = 0.50*(firstxz + lloopxz + lloop2xz + lloop5xz)
        DYZ(ipt,jpt) = 0.50*(firstyz + lloopyz + lloop2yz + lloop5yz)
        DZX(ipt,jpt) = 0.50*(firstzx + lloopzx + lloop2zx + lloop5zx)
        DZY(ipt,jpt) = 0.50*(firstzy + lloopzy + lloop2zy + lloop5zy)
 1100  continue

C  THIS PART CALCULATES DXX(Ipt,I) DXY(Ipt,I) ETC, IE, J = I

       sumkxx = 0.0
       sumkyy = 0.0
       sumkzz = 0.0
       sumkxy = 0.0
       sumkxz = 0.0
       sumkyx = 0.0
       sumkyz = 0.0
       sumkzx = 0.0
       sumkzy = 0.0

       do 1500 k = 1,INDEX2(ipt)
       kpt = LIST2(ipt,k)
       negki = negno(kpt,ipt)
       rik = r(ipt,k)
       rik2 = rik*rik
       rik3 = rik2*rik
       rikx = RX(ipt,k)
       riky = RY(ipt,k)
       rikz = Rz(ipt,k)
       dikvik = ff1(ipt,k)
       dikvki = ff1(kpt,negki)
       dikdikvik = ff11(ipt,k)
       dikdikvki = ff11(kpt,negki)

       firstxx = (1./rik - rikx*rikx/rik3)*(dikvik + dikvki) +
     a      (rikx*rikx/rik2)*(dikdikvik + dikdikvki)
c
       firstyy = (1./rik - riky*riky/rik3)*(dikvik + dikvki) +
     a      (riky*riky/rik2)*(dikdikvik + dikdikvki)
c     
       firstzz = (1./rik - rikz*rikz/rik3)*(dikvik + dikvki) +
     a      (rikz*rikz/rik2)*(dikdikvik + dikdikvki)
c
       firstxy = -(rikx*riky/rik3)*(dikvik + dikvki) +
     a      (rikx*riky/rik2)*(dikdikvik + dikdikvki)
c     
       firstxz = -(rikx*rikz/rik3)*(dikvik + dikvki) +
     a      (rikx*rikz/rik2)*(dikdikvik + dikdikvki)
c
       firstyx = -(riky*rikx/rik3)*(dikvik + dikvki) +
     a      (riky*rikx/rik2)*(dikdikvik + dikdikvki)
c     
       firstyz = -(riky*rikz/rik3)*(dikvik + dikvki) +
     a      (riky*rikz/rik2)*(dikdikvik + dikdikvki)
c
       firstzx = -(rikz*rikx/rik3)*(dikvik + dikvki) +
     a      (rikz*rikx/rik2)*(dikdikvik + dikdikvki)
c
       firstzy = -(rikz*riky/rik3)*(dikvik + dikvki) +
     a      (rikz*riky/rik2)*(dikdikvik + dikdikvki)

      lloop3xx = 0.0
      lloop3yy = 0.0
      lloop3zz = 0.0
      lloop3xy = 0.0
      lloop3xz = 0.0
      lloop3yx = 0.0
      lloop3yz = 0.0
      lloop3zx = 0.0
      lloop3zy = 0.0

      do 1550 l = 1, INDEX2(ipt)
      if (l.eq.k) goto 1550
      ril = r(ipt,l)
      rilx = RX(ipt,l)
      rily = RY(ipt,l)
      rilz = RZ(ipt,l)
      nloop2 = 0.0
      do 1580 m = 1, INDEX2(ipt)
      if ((m.eq.k).or.(m.eq.l)) goto 1580
      dildikvin = pf3(ipt,m)*zet2(ipt,m,l)*zet2(ipt,m,k)
      nloop2 = nloop2 + dildikvin
 1580 continue

      dikvil = ff2(ipt,l,k)
      dikdikvil = ff22(ipt,l,k)
      dikdilvik = ff12(ipt,k,l)
      dikdilvil = ff12(ipt,l,k)
    

      lloop3xx = lloop3xx + (1./rik - rikx*rikx/rik3)*dikvil +
     a     (rikx*rikx/rik2)*dikdikvil +
     a     (rikx*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3yy = lloop3yy + (1./rik - riky*riky/rik3)*dikvil +
     a     (riky*riky/rik2)*dikdikvil +
     a     (riky*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3zz = lloop3zz + (1./rik - rikz*rikz/rik3)*dikvil +
     a     (rikz*rikz/rik2)*dikdikvil +
     a     (rikz*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3xy = lloop3xy  - (rikx*riky/rik3)*dikvil +
     a     (rikx*riky/rik2)*dikdikvil +
     a     (rikx*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3xz = lloop3xz  - (rikx*rikz/rik3)*dikvil +
     a     (rikx*rikz/rik2)*dikdikvil +
     a     (rikx*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3yx = lloop3yx  - (riky*rikx/rik3)*dikvil +
     a     (riky*rikx/rik2)*dikdikvil +
     a     (riky*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3yz = lloop3yz  - (riky*rikz/rik3)*dikvil +
     a     (riky*rikz/rik2)*dikdikvil +
     a     (riky*rilz/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3zx = lloop3zx  - (rikz*rikx/rik3)*dikvil +
     a     (rikz*rikx/rik2)*dikdikvil +
     a     (rikz*rilx/(rik*ril))*(dikdilvik + dikdilvil + nloop2)
c
      lloop3zy = lloop3zy  - (rikz*riky/rik3)*dikvil +
     a     (rikz*riky/rik2)*dikdikvil +
     a     (rikz*rily/(rik*ril))*(dikdilvik + dikdilvil + nloop2)

 1550 continue

      lloop4xx = 0.0
      lloop4yy = 0.0
      lloop4zz = 0.0
      lloop4xy = 0.0
      lloop4xz = 0.0      
      lloop4yx = 0.0
      lloop4yz = 0.0
      lloop4zx = 0.0
      lloop4zy = 0.0
      
      do 1600 l = 1, INDEX2(kpt)
      lpt = LIST2(kpt,l)
      if(lpt.eq.ipt) goto 1600
      rkl = r(kpt,l)
      rkl2 = rkl*rkl
      rklx = RX(kpt,l)
      rkly = RY(kpt,l)
      rklz = RZ(kpt,l)
      rmuil = (RX(kpt,negki)*RX(kpt,l) + RY(kpt,negki)*RY(kpt,l) +
     a     RZ(kpt,negki)*RZ(kpt,l))/(rik*rkl)
      ril2 = rik2 + rkl2 -2.0*rmuil*rik*rkl
      ril = sqrt(ril2)
      ril3 = ril2*ril
      rilx = -1.0*(RX(kpt,negki) - RX(kpt,l))
      rily = -1.0*(RY(kpt,negki) - RY(kpt,l))
      rilz = -1.0*(RZ(kpt,negki) - RZ(kpt,l))
      
      nloop4xx = 0.0
      nloop4yy = 0.0
      nloop4zz = 0.0
      nloop4xy = 0.0
      nloop4xz = 0.0
      nloop4yx = 0.0
      nloop4yz = 0.0
      nloop4zx = 0.0
      nloop4zy = 0.0
      
      do 1650 m = 1, INDEX2(kpt)
      npt = LIST2(kpt,m)
      if ((npt.eq.ipt) .or.(m.eq.l)) goto 1650
      rkn = r(kpt,m)
      rkn2 = rkn*rkn
      rmuin = (RX(kpt,negki)*RX(kpt,m) + RY(kpt,negki)*RY(kpt,m) +
     a     RZ(kpt,negki)*RZ(kpt,m))/(rik*rkn)
      rin2 = rik2 + rkn2 - 2.0*rmuin*rik*rkn
      rin = sqrt(rin2)
      rinx = RX(kpt,m)-RX(kpt,negki)
      riny = RY(kpt,m)-RY(kpt,negki)
      rinz = RZ(kpt,m)-RZ(kpt,negki)
      
      dildinvki = pf3(kpt,negki)*zet3(kpt,negki,l)*zet3(kpt,negki,m)
       
      nloop4xx = nloop4xx + (rilx*rinx/(ril*rin))*dildinvki
      nloop4yy = nloop4yy + (rily*riny/(ril*rin))*dildinvki
      nloop4zz = nloop4zz + (rilz*rinz/(ril*rin))*dildinvki
      nloop4xy = nloop4xy + (rilx*riny/(ril*rin))*dildinvki
      nloop4xz = nloop4xz + (rilx*rinz/(ril*rin))*dildinvki
      nloop4yx = nloop4yx + (rily*rinx/(ril*rin))*dildinvki
      nloop4yz = nloop4yz + (rily*rinz/(ril*rin))*dildinvki
      nloop4zx = nloop4zx + (rilz*rinx/(ril*rin))*dildinvki
      nloop4zy = nloop4zy + (rilz*riny/(ril*rin))*dildinvki    
 1650 continue

      dikvkl = ff2(kpt,l,negki)
      dilvki = ff3(kpt,negki,l)
      dilvkl = ff3(kpt,l,negki)
      dildilvkl = ff33(kpt,l,negki)
      dildikvki = ff13(kpt,negki,l)
      dildikvkl = ff23(kpt,l,negki)
      dildilvki = ff33(kpt,negki,l)
      dikdikvkl = ff22(kpt,l,negki)

      lloop4xx = (1./rik - rikx*rikx/rik3)*dikvkl +
     a               (1./ril -rilx*rilx/ril3)*(dilvki + dilvkl) +
     a   ((rilx*rikx + rikx*rilx)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikx*rikx/rik2)*dikdikvkl +
     a   (rilx*rilx/ril2)*(dildilvki + dildilvkl) +
     a    nloop4xx + lloop4xx
c
      lloop4yy = (1./rik - riky*riky/rik3)*dikvkl +
     a               (1./ril -rily*rily/ril3)*(dilvki + dilvkl) +
     a   ((rily*riky + riky*rily)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (riky*riky/rik2)*dikdikvkl +
     a   (rily*rily/ril2)*(dildilvki + dildilvkl) +
     a    nloop4yy + lloop4yy
c
      lloop4zz = (1./rik - rikz*rikz/rik3)*dikvkl +
     a               (1./ril -rilz*rilz/ril3)*(dilvki + dilvkl) +
     a   ((rilz*rikz + rikz*rilz)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikz*rikz/rik2)*dikdikvkl +
     a   (rilz*rilz/ril2)*(dildilvki + dildilvkl) +
     a    nloop4zz + lloop4zz
c
      lloop4xy = -(rikx*riky/rik3)*dikvkl 
     a              -(rilx*rily/ril3)*(dilvki + dilvkl) +
     a   ((rilx*riky + rikx*rily)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikx*riky/rik2)*dikdikvkl +
     a   (rilx*rily/ril2)*(dildilvki + dildilvkl) +
     a    nloop4xy + lloop4xy
c          
      lloop4xz = -(rikx*rikz/rik3)*dikvkl 
     a              -(rilx*rilz/ril3)*(dilvki + dilvkl) +
     a   ((rilx*rikz + rikx*rilz)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikx*rikz/rik2)*dikdikvkl +
     a   (rilx*rilz/ril2)*(dildilvki + dildilvkl) +
     a    nloop4xz + lloop4xz
c          
      lloop4yx = -(riky*rikx/rik3)*dikvkl 
     a              -(rily*rilx/ril3)*(dilvki + dilvkl) +
     a   ((rily*rikx + riky*rilx)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (riky*rikx/rik2)*dikdikvkl +
     a   (rily*rilx/ril2)*(dildilvki + dildilvkl) +
     a    nloop4yx + lloop4yx
c          
      lloop4yz = -(riky*rikz/rik3)*dikvkl 
     a              -(rily*rilz/ril3)*(dilvki + dilvkl) +
     a   ((rily*rikz + riky*rilz)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (riky*rikz/rik2)*dikdikvkl +
     a   (rily*rilz/ril2)*(dildilvki + dildilvkl) +
     a    nloop4yz + lloop4yz
c 
      lloop4zx = -(rikz*rikx/rik3)*dikvkl 
     a              -(rilz*rilx/ril3)*(dilvki + dilvkl) +
     a   ((rilz*rikx + rikz*rilx)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikz*rikx/rik2)*dikdikvkl +
     a   (rilz*rilx/ril2)*(dildilvki + dildilvkl) +
     a    nloop4zx + lloop4zx
c          
      lloop4zy = -(rikz*riky/rik3)*dikvkl 
     a              -(rilz*rily/ril3)*(dilvki + dilvkl) +
     a   ((rilz*riky + rikz*rily)/(rik*ril))*(dildikvki + dildikvkl) +
     a   (rikz*riky/rik2)*dikdikvkl +
     a   (rilz*rily/ril2)*(dildilvki + dildilvkl) +
     a    nloop4zy + lloop4zy
          
 1600 continue

      sumkxx = sumkxx + firstxx + lloop3xx + lloop4xx
      sumkyy = sumkyy + firstyy + lloop3yy + lloop4yy
      sumkzz = sumkzz + firstzz + lloop3zz + lloop4zz
      sumkxy = sumkxy + firstxy + lloop3xy + lloop4xy
      sumkxz = sumkxz + firstxz + lloop3xz + lloop4xz
      sumkyx = sumkyx + firstyx + lloop3yx + lloop4yx
      sumkyz = sumkyz + firstyz + lloop3yz + lloop4yz
      sumkzx = sumkzx + firstzx + lloop3zx + lloop4zx
      sumkzy = sumkzy + firstzy + lloop3zy + lloop4zy
 1500 continue

      DXX(ipt,ipt) = 0.50*sumkxx
      DYY(ipt,ipt) = 0.50*sumkyy
      DZZ(ipt,ipt) = 0.50*sumkzz
      DXY(ipt,ipt) = 0.50*sumkxy
      DXZ(ipt,ipt) = 0.50*sumkxz
      DYX(ipt,ipt) = 0.50*sumkyx
      DYZ(ipt,ipt) = 0.50*sumkyz
      DZX(ipt,ipt) = 0.50*sumkzx
      DZY(ipt,ipt) = 0.50*sumkzy
 1002 continue
      
C     CHECK RIGID TRANSLATION
      do 2002 ipt = 37,38
      sumDxx = 0.0
      sumDyy = 0.0
      sumDzz = 0.0
      sumDxy = 0.0
      sumDxz = 0.0
      sumDyx = 0.0
      sumDyz = 0.0
      sumDzx = 0.0
      sumDzy = 0.0

      do 2003 jpt = 1, n
      if(jpt.eq.ipt) goto 2003
      sumDxx = sumDxx + Dxx(ipt,jpt)
      sumDyy = sumDyy + Dyy(ipt,jpt)
      sumDzz = sumDzz + Dzz(ipt,jpt)
      sumDxy = sumDxy + Dxy(ipt,jpt)
      sumDxz = sumDxz + Dxz(ipt,jpt)
      sumDyx = sumDyx + Dyx(ipt,jpt)
      sumDyz = sumDyz + Dyz(ipt,jpt)
      sumDzx = sumDzx + Dzx(ipt,jpt)
      sumDzy = sumDzy + Dzy(ipt,jpt)
 2003 continue

      print *, 'i = ',ipt, '        (In eV/A^2)'
      print *,'Dxx: ', -1.0*sumDxx, Dxx(ipt,ipt)
      print *,'Dyy: ', -1.0*sumDyy, Dyy(ipt,ipt)
      print *,'Dzz: ', -1.0*sumDzz, Dzz(ipt,ipt)
      print *,'Dxy: ', -1.0*sumDxy, Dxy(ipt,ipt)
      print *,'Dxz: ', -1.0*sumDxz, Dxz(ipt,ipt)
      print *,'Dyx: ', -1.0*sumDyx, Dyx(ipt,ipt)
      print *,'Dyz: ', -1.0*sumDyz, Dyz(ipt,ipt)
      print *,'Dzx: ', -1.0*sumDzx, Dzx(ipt,ipt)
      print *,'Dzy: ', -1.0*sumDzy, Dzy(ipt,ipt)
      print *, '        '
 2002 continue

      print *,'Dxx(1,2) = ',Dxx(1,2),
     a        '             Dxy(1,2) = ',Dxy(1,2)
      print *,' '

C     THE FOLLOWING CALCULATES THE HESSIAN FROM
C     DXX(Ipt,Jpt),DXY(Ipt,Jpt) ETC. OBTAINED ABOVE. ALL IT
C     DOES IS FORM HESS(A,B) FROM THE FOLLOWING:
C     [X1,Y1,Z1,X2,Y2,Z2...XN,YN,ZN]^T*[X1,Y1,Z1,X2,Y2,Z2...XN,YN,ZN]
C     WHICH GIVES A 3N X 3N MATRIX. FOR EXAMPLE, 
C     HESS(1,1) = DXX(1,1), HESS(1,2) = DXY(1,1), HESS(2,9) = DYZ(1,3)

      do 28 ipt = 1,n
      do 28 jpt = 1,n
      ia = 3*(ipt-1)
      ib = 3*(jpt-1)
      smij = sqrt(pmass(ipt)*pmass(jpt))
      
      HESS(1+ia,1+ib) = Dxx(ipt,jpt)/smij
      HESS(1+ia,2+ib) = Dxy(ipt,jpt)/smij
      HESS(1+ia,3+ib) = Dxz(ipt,jpt)/smij
      HESS(2+ia,1+ib) = Dyx(ipt,jpt)/smij
      HESS(2+ia,2+ib) = Dyy(ipt,jpt)/smij
      HESS(2+ia,3+ib) = Dyz(ipt,jpt)/smij
      HESS(3+ia,1+ib) = Dzx(ipt,jpt)/smij
      HESS(3+ia,2+ib) = Dzy(ipt,jpt)/smij
      HESS(3+ia,3+ib) = Dzz(ipt,jpt)/smij
 28   continue

      DO 57 I=1,6
      IPT = 1+(I-1)/3
      DO 57 J=1,6
      JPT = 1+(J-1)/3
      MAXIJ = 0
      DO 96 IP=0,N-1,2
      IF (HESS(I,J+3*IP).NE.0) THEN
      MAXIJ = MAXIJ+1
      DM(I,J,MAXIJ,1) = RRX(IPT,JPT+IP)
      DM(I,J,MAXIJ,2) = RRY(IPT,JPT+IP)
      DM(I,J,MAXIJ,3) = RRZ(IPT,JPT+IP)
      DM(I,J,MAXIJ,4) = HESS(I,J+3*IP)
      ENDIF
 96   CONTINUE
      MAX_DM(I,J) = MAXIJ
 57   CONTINUE

c      read *,i,j
c      print *,'hess(i,j) = ',hess(i,j)
c      print *,'dxx12 = ', dxx(1,2)
c      print *,'dxx11 = ', dxx(1,1)

      return
      end

      SUBROUTINE PHONON(KX,KY,KZ,ORDERED)
      IMPLICIT REAL (A-H,O-Z)
      REAL KX,KY,KZ
      INTEGER ORDERED
      PARAMETER (NM=70,NX=10,PI=3.1415926536 )
      COMMON /BLOCK3/ DXX(nm,nm),DXY(nm,nm),DYX(nm,nm),DYY(nm,nm),
     a                DZZ(nm,nm),DXZ(nm,nm),DYZ(nm,nm),DZX(nm,nm),
     a                DZY(nm,nm),HESS(NM*3,NM*3),freq(NM*3),FREQ_K(12),
     a                DM(6,6,NM/2,4),MAX_DM(6,6),S_FREQ_K(NM*6),
     a                amp_si(12),amp_c(12)

      PARAMETER (TOL=1.E-4)
      REAL A(6,6),B(6,6),C(12,12),D(6),D_MODES(12,12)
      INTEGER NK(12)

      DO 82 I=1,6
      DO 82 J=I,6
      A(I,J) = 0.
      B(I,J) = 0.
      DO 23 M = 1,MAX_DM(I,J)
      PHASE = KX*DM(I,J,M,1)+KY*DM(I,J,M,2)+KZ*DM(I,J,M,3)
      A(I,J) = A(I,J) + DM(I,J,M,4)*COS(PHASE)
      B(I,J) = B(I,J) + DM(I,J,M,4)*SIN(PHASE)
 23   CONTINUE
      A(J,I) = A(I,J)
      B(J,I) = -B(I,J)

      C(I,J) = A(I,J)
      C(I+6,J+6) = A(I,J)
      C(I,J+6) = -B(I,J)
      C(I+6,J) = B(I,J)

      C(J,I) = A(J,I)
      C(J+6,I+6) = A(J,I)
      C(J,I+6) = -B(J,I)
      C(J+6,I) = B(J,I)
 82   CONTINUE

      CALL JACOBI(C,12,12,FREQ_K,D_MODES,NROT)
      DO 33 I=1,12
      IF (FREQ_K(I).LT.0) FREQ_K(I) = 0.
      FREQ_K(I) = SQRT(FREQ_K(I))
      amp_si(i) = d_modes(1,i)*d_modes(1,i)+d_modes(2,i)*d_modes(2,i)+
     a            d_modes(3,i)*d_modes(3,i)+d_modes(7,i)*d_modes(7,i)+
     a            d_modes(8,i)*d_modes(8,i)+d_modes(9,i)*d_modes(9,i)
      amp_c(i)  = 1-amp_si(i)
 33   CONTINUE

      IF (ORDERED.EQ.1) THEN
C     ----------------------------------------------------------------
C     WE SHALL PICK 6 INDEPENDENT NORMAL MODES OUT OF 12 (A HERMITIAN 
C     6*6 MATRIX IS EQUIVALENT TO A 12*12 REAL SYMMETRIC MATRIX. SEE
C     NUMERICAL RECIPES), AND ORDER THEM BY FREQUENCY IN ASCENDING 
C     ORDER.

      IP = 1
      NK(1) = 1
      DO 70 I=2,12
       DO 80 J = 1,IP
       SUM = 0.
       DO 90 L=1,6
 90    SUM = SUM+D_MODES(L,NK(J))*D_MODES(L+6,I)
     a          -D_MODES(L+6,NK(J))*D_MODES(L,I)
       IF (ABS(ABS(SUM)-1.).LE.TOL) GOTO 70
 80    CONTINUE

       J = 1
 92    IF ((J.GT.IP).OR.(FREQ_K(I).LT.FREQ_K(NK(J)))) GOTO 21
       J=J+1
       GOTO 92

 21    DO 100 L=IP,J,-1
 100   NK(L+1) = NK(L)
       NK(J) = I
       IP=IP+1
 70   CONTINUE

      IF (IP.NE.6) THEN 
      PRINT *,'SOMETHING IS WRONG WITH SORTING'
      STOP
      ENDIF

      DO 110 I=1,6
      D(I) = FREQ_K(NK(I))
c     print *, 'mode ', I, ' =', D(I)
      DO 130 J=1,6
      A(J,I) = D_MODES(J,NK(I))
      B(J,I) = D_MODES(J+6,NK(I))
c     A stores the real part of the 6 ordered modes
c     B stores the imaginary part of the 6 ordered modes
c     print *, A(J,I), '+', B(J,I),'i'
 130  CONTINUE
c     print *,' '
 110  CONTINUE

      DO 120 I=1,6
      FREQ_K(I) = D(I)
      FREQ_K(I+6) = 0
 120  CONTINUE
      ENDIF

      RETURN 
      END

      SUBROUTINE SUPERCELL(KX,KY,KZ)
      IMPLICIT REAL (A-H,O-Z)
      REAL KX,KY,KZ
      PARAMETER (NM=70,NX=10,PI=3.1415926536 )
      COMMON /BLOCK1/ FX(NM),FY(NM),FZ(NM),PM(2),EP(NM),negno(nm,nm),
     A                NTYPE(NM),LIST(NM,NX),INDEX(NM),PMASS(NM),X(NM),
     A                Y(NM),Z(NM),DX(NM),DY(NM),DZ(NM),DXOLD(NM),
     A                DYOLD(NM),DZOLD(NM),XOLD(NM),YOLD(NM),UPDATE,
     A                ZOLD(NM),LIST2(NM,NX),INDEX2(NM),N,MODE,NEW,
     A                HSUM(3,3),RRX(2,NM),RRY(2,NM),RRZ(2,NM),WRITE_LP,
     A                POUT(3,3),VIRIAL(3,3),CJX,CJY,CJZ,XF,CUBE0,
     A                VOLUME,VR,RCUT0,RCUT,POTE,C11,C12,C44,CUBE
      COMMON /BLOCK2/ PA(2), PB(2),PL1(2),PL2(2), PBT(2),PN(2),
     A                PC(2),PWR2(2),PD(2), PH(2), PR(2), PS(2), 
     A                C2D2(2), PWR(2),PWR1(2),CHIJ(NX),RLM1(NX),
     A                RLM2(NX),CA(NX),CB(NX),CR(NX),CS(NX),RLIST,
     A                CUBE2,CHI_SI_C,FCUBE2
      COMMON /BLOCK3/ DXX(nm,nm),DXY(nm,nm),DYX(nm,nm),DYY(nm,nm),
     a                DZZ(nm,nm),DXZ(nm,nm),DYZ(nm,nm),DZX(nm,nm),
     a                DZY(nm,nm),HESS(NM*3,NM*3),freq(NM*3),FREQ_K(12),
     a                DM(6,6,NM/2,4),MAX_DM(6,6),S_FREQ_K(NM*6),
     a                amp_si(12), amp_c(12)

      REAL A(3*NM,3*NM),B(3*NM,3*NM),C(6*NM,6*NM),S_MODES(6*NM,6*NM)

      N3 = 3*N
      DO 82 I=1,N3
      IPT = (I-1)/3+1
      XI = X(IPT)
      YI = Y(IPT)
      ZI = Z(IPT)
      DO 82 J=I,N3
      JPT = (J-1)/3+1
      A(I,J) = 0.
      B(I,J) = 0.
      RXX  = X(JPT) - XI
      RYY  = Y(JPT) - YI
      RZZ  = Z(JPT) - ZI
      IF(RXX.GT.CUBE2)  RXX=RXX-CUBE
      IF(RXX.LT.FCUBE2) RXX=RXX+CUBE
      IF(RYY.GT.CUBE2)  RYY=RYY-CUBE
      IF(RYY.LT.FCUBE2) RYY=RYY+CUBE
      IF(RZZ.GT.CUBE2)  RZZ=RZZ-CUBE
      IF(RZZ.LT.FCUBE2) RZZ=RZZ+CUBE
      PHASE = KX*RXX+KY*RYY+KZ*RZZ
      A(I,J) = HESS(I,J)*COS(PHASE)
      B(I,J) = HESS(I,J)*SIN(PHASE)
      A(J,I) = A(I,J)
      B(J,I) = -B(I,J)
      C(I,J) = A(I,J)
      C(I+N3,J+N3) = A(I,J)
      C(I,J+N3) = -B(I,J)
      C(I+N3,J) = B(I,J)
      C(J,I) = A(J,I)
      C(J+N3,I+N3) = A(J,I)
      C(J,I+N3) = -B(J,I)
      C(J+N3,I) = B(J,I)
 82   CONTINUE

      CALL JACOBI(C,6*N,6*NM,S_FREQ_K,S_MODES,NROT)
      DO 33 I=1,6*N
      IF (S_FREQ_K(I).LT.0) S_FREQ_K(I) = 0.
      S_FREQ_K(I) = SQRT(S_FREQ_K(I))
 33   CONTINUE

      RETURN 
      END

      SUBROUTINE jacobi(a,n,np,d,v,nrot)
      INTEGER n,np,nrot,NMAX
      REAL a(np,np),d(np),v(np,np)
      PARAMETER (NMAX=5000)
      INTEGER i,ip,iq,j
      REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
      do 12 ip=1,n
        do 11 iq=1,n
          v(ip,iq)=0.
11      continue
        v(ip,ip)=1.
12    continue
      do 13 ip=1,n
        b(ip)=a(ip,ip)
        d(ip)=b(ip)
        z(ip)=0.
13    continue
      nrot=0
      do 24 i=1,50
        sm=0.
        do 15 ip=1,n-1
          do 14 iq=ip+1,n
            sm=sm+abs(a(ip,iq))
14        continue
15      continue
c        if(sm.eq.0.)return
        if(sm.lt.1.e-7) return
        if(i.lt.4)then
          tresh=0.2*sm/n**2
        else
          tresh=0.
        endif
        do 22 ip=1,n-1
          do 21 iq=ip+1,n
            g=100.*abs(a(ip,iq))
            if((i.gt.4).and.(abs(d(ip))+
     *g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then
              a(ip,iq)=0.
            else if(abs(a(ip,iq)).gt.tresh)then
              h=d(iq)-d(ip)
              if(abs(h)+g.eq.abs(h))then
                t=a(ip,iq)/h
              else
                theta=0.5*h/a(ip,iq)
                t=1./(abs(theta)+sqrt(1.+theta**2))
                if(theta.lt.0.)t=-t
              endif
              c=1./sqrt(1+t**2)
              s=t*c
              tau=s/(1.+c)
              h=t*a(ip,iq)
              z(ip)=z(ip)-h
              z(iq)=z(iq)+h
              d(ip)=d(ip)-h
              d(iq)=d(iq)+h
              a(ip,iq)=0.
              do 16 j=1,ip-1
                g=a(j,ip)
                h=a(j,iq)
                a(j,ip)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
16            continue
              do 17 j=ip+1,iq-1
                g=a(ip,j)
                h=a(j,iq)
                a(ip,j)=g-s*(h+g*tau)
                a(j,iq)=h+s*(g-h*tau)
17            continue
              do 18 j=iq+1,n
                g=a(ip,j)
                h=a(iq,j)
                a(ip,j)=g-s*(h+g*tau)
                a(iq,j)=h+s*(g-h*tau)
18            continue
              do 19 j=1,n
                g=v(j,ip)
                h=v(j,iq)
                v(j,ip)=g-s*(h+g*tau)
                v(j,iq)=h+s*(g-h*tau)
19            continue
              nrot=nrot+1
            endif
21        continue
22      continue
        do 23 ip=1,n
          b(ip)=b(ip)+z(ip)
          d(ip)=b(ip)
          z(ip)=0.
23      continue
24    continue
      pause 'too many iterations in jacobi'
      return
      END

      SUBROUTINE CALC_DOS(FREQ,NF,DOS,NPA)
      IMPLICIT REAL (A-H,O-Z)
      REAL FREQ(NF),DOS(NPA)
      INTEGER NF,NPA
      DATA LP,LP_CONFIG,LP_DOS,LP_DM  /25,26,27,28/

      FMAX = 0.
      DO 71 I=1,NF
 71   IF (FREQ(I).GT.FMAX) FMAX = FREQ(I)

      PA = FMAX/NPA
      DO 72 I=1,NPA
 72   DOS(I) = 0.

      DO 73 I=1,NF
c     int(0.7) = 0
      IP = INT(FREQ(I)/PA)+1
      DOS(IP) = DOS(IP)+1
 73   CONTINUE

      OPEN (UNIT=LP_DOS,  STATUS='UNKNOWN', FORM='FORMATTED',
     a      FILE="dos.out")
      WRITE (LP_DOS,501) NF
      WRITE (LP_DOS,501) NPA
      WRITE (LP_DOS,502) PA
 502  FORMAT (F9.6)
      DO 74 IP = 1,NPA
 74   WRITE (LP_DOS,501) INT(DOS(IP))
 501  FORMAT (I8)

      RETURN 
      END

      SUBROUTINE CALC_LDOS(FREQ,SAMP_SI,NF,NPA)
      IMPLICIT REAL (A-H,O-Z)
      REAL FREQ(NF),SAMP_SI(NF)
      INTEGER NF,NPA
      DIMENSION DOS_SI(2000),DOS_C(2000)

      FMAX = 0.
      DO 71 I=1,NF
 71   IF (FREQ(I).GT.FMAX) FMAX = FREQ(I)

      print *,'maximum frequency = ',fmax
      
      
      PA = FMAX/NPA
      DO 72 I=1,NPA
      DOS_SI(I) = 0.
      DOS_C(I) = 0.
 72   CONTINUE

      DO 73 I=1,NF
c     int(0.7) = 0
      IP = INT(FREQ(I)/PA)+1
      DOS_SI(IP) = DOS_SI(IP)+SAMP_SI(I)
c      DOS_C(IP) = DOS_C(IP)+1.-SAMP_SI(I)
 73   CONTINUE

      LP_DOS_SI = 28
c      LP_DOS_C = 29

      OPEN (UNIT=LP_DOS_SI,  STATUS='UNKNOWN', FORM='FORMATTED',
     a      FILE="dos_old.out")
c      OPEN (UNIT=LP_DOS_C,  STATUS='UNKNOWN', FORM='FORMATTED',
c     a      FILE="dos_c.out")
      WRITE (LP_DOS_SI,501) NF
      WRITE (LP_DOS_SI,501) NPA+1
      WRITE (LP_DOS_SI,502) PA
c      WRITE (LP_DOS_C,501) NF
c      WRITE (LP_DOS_C,501) NPA
c      WRITE (LP_DOS_C,502) PA
      DO 74 IP = 1,NPA+1
      WRITE (LP_DOS_SI,503) DOS_SI(IP)
c      WRITE (LP_DOS_C,503)  DOS_C(IP)
 74   continue
 501  FORMAT (I8)
 502  FORMAT (F9.6)
 503  FORMAT (E13.7)

      RETURN 
      END

      FUNCTION ran1(idum)
      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

