Graphene Susceptibility Code

c makes susceptibility from graphene bands.
      parameter(nq=300)
      parameter(ne=1000)
      dimension xx(ne)
      emax=7.0
      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 loop over momentum transfers.
      dq0=2.0*pi/3.0/4.0
      do 8 iq=1,4
      q0x=0.0
      q0y=iq*dq0
c      q0y=(iq+4)*dq0
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(q0x,q0y,qx,qy,eq,e1,e2,p)
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)*p
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
      write(*,*) x
c output result with hilbert transform.
      do 7  i=1,ne
      x=0.0
      do 6 ii=1,ne
      y=i+ii
      x=x-xx(ii)/y
      if(i.eq.ii) go to 6
      y=i-ii
      x=x+xx(ii)/y
    6 continue
      y=-xx(i)*pi
      e=(i-0.5)*de
      write(7,20) e,x,y
    7 continue
      write(7,21)
    8 continue
   20 format(f10.4,2e18.7)
   21 format(1h )
      stop
      end 
      subroutine bands(q0x,q0y,qx,qy,e,e1,e2,p)
      complex zm,zp,z1,z2,z3
      a=sqrt(3.0)/2.0
      b=1.0/2.0
c compute energy and wavefunction for lower band
      z1=cexp(cmplx(0.0,qx*a+qy*b))
      z2=cexp(cmplx(0.0,-qx*a+qy*b))
      z3=cexp(cmplx(0.0,-qy))
      zm=z1+z2+z3
      em=cabs(zm)
      em1=a*aimag(conjg(zm)*(z1-z2))/em
      em2=a*aimag(conjg(zm)*(z1-z3))/em
c compute energy and wavefunction for upper band.      
      qpx=qx+q0x
      qpy=qy+q0y
      z1=cexp(cmplx(0.0,qpx*a+qpy*b))
      z2=cexp(cmplx(0.0,-qpx*a+qpy*b))
      z3=cexp(cmplx(0.0,-qpy))
      zp=z1+z2+z3
      ep=cabs(zp)
      ep1=a*aimag(conjg(zp)*(z1-z2))/ep
      ep2=a*aimag(conjg(zp)*(z1-z3))/ep
c compute matrix element
      p=(1.0-real(conjg(zp)*zm)/em/ep)/2.0
c compute transition energy and gradient
      e=em+ep
      e1=em1+ep1
      e2=em2+ep2
      return
      end