program tba
      implicit none
! input variables: energies of p and s orbitals, bond distance, 
      real ep(2), es(2), d, z
      real eh(2), v3, v2h, eb, ea, vsss, vppp, width_valence
      real width_conduction, cb_min, cb_max, vb_min, vb_max, gap
      print*, "This program calculates the band gap "
      print*, "for compound tetrahedral solids"
      print*, " "
      print*, "Give s and p energies of the first atom (in eV)"
      read*, es(1), ep(1)
      print*, "Give s and p energies of the second atom (in eV)"
      read*, es(2), ep(2)
      print*, "Give the bond distance (in Angstroms)"
      read*, d
      print*, " "
! 
! calculate the energies of the hybrids
      eh(1)=0.25*es(1)+0.75*ep(1)
      eh(2)=0.25*es(2)+0.75*ep(2)
! calculate v3
      v3=0.5*(eh(2)-eh(1))
! calculate v2h, using hbar^2 / (m a0^2) = 27.2 eV
      v2h=-3.22*27.2*(0.529/d)**2
! calculate the energies for the bonding and antibonding state of the AB
! molecule
      eb=0.5*(eh(2)+eh(1))-sqrt(v2h**2+v3**2)
      ea=0.5*(eh(2)+eh(1))+sqrt(v2h**2+v3**2)
      print '( "energy of sp3 orbital on atom 1: ", g8.2)', eh(1)
      print '( "energy of sp3 orbital on atom 2: ", g8.2)', eh(2)
      print '( "V3h, V2h = ", 2g8.2)', v3, v2h
      print '( "energy of bonding orbital:     " , g8.2)', eb
      print '( "energy of antibonding orbital: " , g8.2)', ea
      print*, " " 
! now the widths of the valence and the conduction band
      vsss=0.79*27.2*(0.529/d)**2
! for the p orbitals, we need to multiply by cos(pi-theta)=-cos(theta)=1/3
      vppp=-0.24*27.2*(0.529/d)**2
! z is the number of neighbors for each atom. Here z=4.
      z=4
      width_valence=2*z*abs(vsss)      
      width_conduction=2*z*abs(vppp)
! calculate the maxima and minima of the bands
      print '( "V2(v), V2(c) = ", 2g8.2)', vsss, vppp
      cb_min=ea-0.5*width_conduction      
      cb_max=ea+0.5*width_conduction      
      vb_min=eb-0.5*width_valence    
      vb_max=eb+0.5*width_valence      

      gap=cb_min-vb_max
! print results

      print '( "Valence band minimum    ", g8.2)', vb_min
      print '( "Valence band maximum    ", g8.2)', vb_max
      print '( "Conduction band minimum ", g8.2)', cb_min
      print '( "Conduction band maximum ", g8.2)', cb_max
      print*, " "
      if (gap .le. 0.0) then
          print*, "this is a metal"
      else
          print '( "Band gap: ", g8.2)', gap 
      endif
      end