Graphene Density of States Code

c makes density of states from graphene bands.
      parameter(nq=300)
      parameter(ne=1000)
      dimension xx(ne)
      emax=3.5
      eps=0.0000001
      pi=3.1415927
c determine energy increment.
      x=ne
      de=emax/x
c determine momentum increment.
      x=nq
      dq=4.0*pi/3.0/x
c calculate sin and cos once.
      s=sqrt(3.0)/2.0
      c=1.0/2.0
c open output file.
      open(7,file="output.dat",status="new")
c zero out accumulator.
      do 1 i=1,ne
      xx(i)=0.0
    1 continue
c begin looping over brillouin zone.
      do 4 j=1,nq
      do 4 k=1,nq
      qx=(j-0.5)*dq+(k-0.5)*dq*c
      qy=(k-0.5)*dq*s
      call bands(qx,qy,eq,e1,e2)
c order the gradients.
      e1=abs(e1)+eps
      e2=abs(e2)+eps
c calculate minimum and maximum energies
      eq=eq-(e1+e2)*dq/2.0
      im=eq/de-2.0
      if(im.lt.1) im=1
      ip=im+(e1+e2)*dq/de+4.0
      if(ip.gt.ne) ip=ne
c loop over energies in this range.
      do 3 i=im,ip
c for this energy calculate the bin edges.
      em=(i-1.0)*de
      ep=em+de
c now calculate two areas.
      dq1=(em-eq)/e1
      dq2=(em-eq)/e2
      am=dq1*dq2
      if(dq1.gt.dq) am=am-(dq1-dq)**2*dq2/dq1
      if(dq2.gt.dq) am=am-(dq2-dq)**2*dq1/dq2
      if((dq1*dq2).gt.((dq1+dq2)*dq)) am=2.0*dq**2
      if(dq1.lt.0.0.or.dq2.lt.0.0) am=0.0
      dq1=(ep-eq)/e1
      dq2=(ep-eq)/e2
      ap=dq1*dq2
      if(dq1.gt.dq) ap=ap-(dq1-dq)**2*dq2/dq1
      if(dq2.gt.dq) ap=ap-(dq2-dq)**2*dq1/dq2
      if((dq1*dq2).gt.((dq1+dq2)*dq)) ap=2.0*dq**2
      if(dq1.lt.0.0.or.dq2.lt.0.0) ap=0.0
c the relevant area is the difference of these.
      xx(i)=xx(i)+abs(ap-am)
c end of loop on energy bins.
    3 continue    
c end of loop over momenta.
    4 continue
c output results.
      x=0.0
      y=(3.0/4.0/pi)**2/2.0
      do 5 i=1,ne
      xx(i)=xx(i)*y
      x=x+xx(i)
      xx(i)=xx(i)/de
    5 continue
      do 6  i=1,ne
      e=(i-0.5-ne)*de
      write(7,20) e,xx(ne-i+1)
    6 continue
      do 7  i=1,ne
      e=(i-0.5)*de
      write(7,20) e,xx(i)
    7 continue
   20 format(f10.4,e18.7)
      write(*,*) x
      stop
      end 
      subroutine bands(qx,qy,e,e1,e2)
      complex z,z1,z2,z3
      a=sqrt(3.0)/2.0
      b=1.0/2.0
      z1=cexp(cmplx(0.0,qx*a+qy*b))
      z2=cexp(cmplx(0.0,-qx*a+qy*b))
      z3=cexp(cmplx(0.0,-qy))
      z=z1+z2+z3
      e=cabs(z)
      e1=a*aimag(conjg(z)*(z1-z2))/e
      e2=a*aimag(conjg(z)*(z1-z3))/e      
      return
      end