      program chadicohen
      parameter (n=30)
      complex g0,g0s,g1,g1s,g2,g2s,g3,g3s
      complex h(8,8),zero
      real kx,ky,kz,kxp,kyp,kzp
      real kxpoint(100),kypoint(100),kzpoint(100),xkk(100)
      real hr(8,8),hi(8,8),w(8),ww(3*n+1,8),wx(3*n+1)
      real zr(8,8),zi(8,8),aa(8),bb(8),cc(2,8)
      pi=3.1415926535897931
      zero=complex(0.0,0.0)
      open(1,file="chadicoh.dat")
      write(*,*) "zadaj krystal podla cisla"
      write(*,*) "1: C"
      write(*,*) "2: Si"
      write(*,*) "3: Ge"
      write(*,*) "4: GaAs"
      write(*,*) "5: ZnSe"
      read(*,*) nn
c     References:
c     D.J. Chadi and M.L. Cohen, "Tight Binding Calculations
c          of the Valence Bands of Diamond and Zincblende Crystals."
c     Phys. Stat. Sol. (b) 68, 405 (1975) 
c     Define a set of k points along the Brillouin zone boundary from
c     L (0.5,0.5,0.5) to Gamma (0,0,0) to X (1.,0.,0.). Space the points
c     evenly, based on the scaling parameter n.
c      Tight binding parameters; these are in eV:
c     structure == 'C':
      if (nn.EQ.1) then
      e_s_c = 0.0 
      e_s_a = 0.0     ! Arbitrary;
      e_p_c = 7.4 - e_s_c
      e_p_a = 7.40 - e_s_c
      v_ss = -15.2
      v_sc_p = 10.25
      v_sa_p = 10.25
      v_xx = 3.0
      v_xy = 8.30
c     structure == 'Si':
      else if (nn.EQ.2) then 
      e_s_c = 0.0
      e_s_a = 0.0     ! Arbitrary
      e_p_c = 7.20 - e_s_c
      e_p_a = 7.20 - e_s_c
      v_ss = -8.13
      v_sc_p = 5.88
      v_sa_p = 5.88
      v_xx = 3.17
      v_xy = 7.51
c     structure == 'Ge':
      else if (nn.EQ.3) then 
      e_s_c = 0.0
      e_s_a = 0.0     ! Arbitrary
      e_p_c = 8.41 - e_s_c
      e_p_a = 8.41 - e_s_c
      v_ss = -6.78
      v_sc_p = 5.31
      v_sa_p = 5.31
      v_xx = 2.62
      v_xy = 6.82
c     structure == 'GaAs':
      else if (nn.EQ.4) then
      e_s_c = -6.01
      e_s_a = -4.79
      e_p_c = 0.19
      e_p_a = 4.59
      v_ss = -7.00
      v_sc_p = 7.28
      v_sa_p = 3.70
      v_xx = 0.93
      v_xy = 4.72
c     structure == 'ZnSe':
      else if (nn.EQ.5) then
      e_s_c = -8.92
      e_s_a = -0.28
      e_p_c = 0.12
      e_p_a = 7.42
      v_ss = -6.14
      v_sc_p = 5.47
      v_sa_p = 4.73
      v_xx = 0.96
      v_xy = 4.38
      end if
c
      kx=0.5   ! Start at the L point (1/2,1/2,1/2)
      ky=0.5
      kz=0.5
      step = 0.5/float(n)
      xk=0.0
      do i=1,3*n+1
c         write(*,*) i,kx,ky,kz
         kxp=kx*pi/2.0
         kyp=ky*pi/2.0
         kzp=kz*pi/2.0
         xkk(i)=xk
         kxpoint(i)=kx
         kypoint(i)=ky
         kzpoint(i)=kz
         g0_r = cos(kxp)*cos(kyp)*cos(kzp)
         g0_i = -sin(kxp)*sin(kyp)*sin(kzp)
         g1_r = -cos(kxp)*sin(kyp)*sin(kzp)
         g1_i = sin(kxp)*cos(kyp)*cos(kzp)
         g2_r = -sin(kxp)*cos(kyp)*sin(kzp)
         g2_i = cos(kxp)*sin(kyp)*cos(kzp)
         g3_r = -sin(kxp)*sin(kyp)*cos(kzp)
         g3_i = cos(kxp)*cos(kyp)*sin(kzp)
c     "s" stands for "star": the complex conjugate
         g0 = complex(g0_r,g0_i)
         g0s = conjg(g0)
         g1 = complex(g1_r,g1_i)
         g1s = conjg(g1)     
         g2 = complex(g2_r,g2_i)
         g2s = conjg(g2)
         g3 = complex(g3_r,g3_i)
         g3s = conjg(g3)
c
         do k=1,8
            do l=1,8
               h(k,l)=zero
            end do
         end do 
c
c     Make the diagonal elements
          H(1,1) = e_s_c
          H(2,2) = e_s_a
          H(3,3) = e_p_c
          H(4,4) = H(3,3)
          H(5,5) = H(3,3)
          H(6,6) = e_p_a
          H(7,7) = H(6,6)
          H(8,8) = H(6,6) 
c     Make the off-diagonal parts
          H(2,1) = v_ss*g0s
          H(1,2) = v_ss*g0

          H(3,2) = -v_sa_p*g1
          H(2,3) = -v_sa_p*g1s
          H(4,2) = -v_sa_p*g2
          H(2,4) = -v_sa_p*g2s
          H(5,2) = -v_sa_p*g3
          H(2,5) = -v_sa_p*g3s

          H(6,1) = v_sc_p*g1s
          H(1,6) = v_sc_p*g1
          H(7,1) = v_sc_p*g2s
          H(1,7) = v_sc_p*g2
          H(8,1) = v_sc_p*g3s
          H(1,8) = v_sc_p*g3

          H(6,3) = v_xx*g0s
          H(3,6) = v_xx*g0
          H(7,3) = v_xy*g3s
          H(3,7) = v_xy*g3
          H(8,3) = v_xy*g2s
          H(3,8) = v_xy*g2

          H(6,4) = v_xy*g3s
          H(4,6) = v_xy*g3
          H(7,4) = v_xx*g0s
          H(4,7) = v_xx*g0
          H(8,4) = v_xy*g1s
          H(4,8) = v_xy*g1

          H(6,5) = v_xy*g2s
          H(5,6) = v_xy*g2
          H(7,5) = v_xy*g1s
          H(5,7) = v_xy*g1
          H(8,5) = v_xx*g0s
          H(5,8) = v_xx*g0
c
          do k=1,8
             do l=1,8
                HR(k,l)=real(H(k,l))
                HI(k,l)=aimag(H(k,l))
             end do
          end do
c
          call CH(8,8,hr,hi,w,0,zr,zi,aa,bb,cc,ierr) ! main procedure
c          write(*,*) ierr
          i1=2
c          go to 20
  10      continue
               if (w(i1-1).GT.w(i1)) then
                  p=w(i1-1)
                  w(i1-1)=w(i1)
                  w(i1)=p
                  i1=2
                  go to 10
               end if
               i1=i1+1
               if (i1.LE.8) go to 10
          do k=1,8
               ww(i,k)=w(k)
          end do
  20      continue
          write(1,*) (w(k),k=1,6)
c
          if (i.LE.n) then
            kx=kx-step  ! Move to the Gamma point (0,0,0)
            ky=ky-step
            kz=kz-step
            xk=xk+step
          else if ((i.GT.n).AND.(i.LT.2*n)) then
            kx=kx+2.0*step ! Now go to the X point (1,0,0)
            ky=0.0
            kz=0.0
            xk=xk+step
          else if (i.EQ.2*n) then 
            kx=1.0             ! Jump to the U,K point
            ky=1.0
            kz=0.0
            xk=xk+step
          else if (i.GT.2*n) then
            kx=kx-2.0*step     ! Now go back to Gamma
            ky=ky-2.0*step
            kz=0.0
            xk=xk+step
          end if
      end do
      close(1)
c
c     stop
      call metafl('cons')
      CALL DISINI
      CALL NAME('k','X')
      CALL NAME('E(k) (eV)','Y')
      call title
      CALL HSYMBL(12)
      CALL GRAF(0.05,1.5,1.8,0.2,-15.,16.,-15.,2.0)
      do i=1,8
        do j=1,3*n+1
           wx(j)=ww(j,i)
        end do
        CALL COLOR('RED')
        CALL curve(xkk,wx,3*n+1)
        call color('GREEN')
        do k=1,3*n+1
           call rlsymb(21,xkk(k),wx(k))
        end do
      end do
      call rline(0.5,-15.0,0.5,16.0)
c      call rline(1.0-step,-15.0,1.0-step,16.0)
      call rline(1.0,-15.0,1.0,16.0)
      call rlmess("L",0.05,-15.5)
      call rlmess("GAMMA",0.4,-15.5)
      call rlmess("X",1.0,-15.5)
      call rlmess("GAMMA",1.4,-15.5)
      call disfin
      end