Graphene Band Structure Code

c makes graphene bands.
      pi=3.1415927
      qk=4.0*pi/3.0/sqrt(3.0)
      qm=2.0*pi/3
c open output files.
      open(7,file="out1.dat",status="new")
      open(8,file="out2.dat",status="new")
c loop from gamma to k.
      dqx=qk/20.0
      dqy=0.0
      dq=sqrt(dqx**2+dqy**2)
      qx=-dqx
      qy=-dqy
      q=-dq
      do 1 i=1,21      
      qx=qx+dqx
      qy=qy+dqy
      q=q+dq
      e=ep(qx,qy)
      write(7,20) q,e,-e
    1 continue
c make vertical line at k.
      e=3.2
      write(8,21) q,e
      e=-e
      write(8,21) q,e
      write(8,22)
c loop from k to m.
      dqx=(sqrt(3.0)/2.0*qm-qk)/20.0
      dqy=0.5*qm/20.0
      dq=sqrt(dqx**2+dqy**2)
      do 2 i=1,20
      qx=qx+dqx
      qy=qy+dqy
      q=q+dq
      e=ep(qx,qy)
      write(7,20) q,e,-e
    2 continue
c make vertcal line at m.
      e=3.2
      write(8,21) q,e
      e=-e
      write(8,21) q,e
      write(8,22)
c loop from m to gamma.
      dqx=-sqrt(3.0)/2.0*qm/20.0
      dqy=-0.5*qm/20.0
      dq=sqrt(dqx**2+dqy**2)
      do 3 i=1,20
      qx=qx+dqx
      qy=qy+dqy
      q=q+dq
      e=ep(qx,qy)
      write(7,20) q,e,-e
    3 continue
   20 format(f10.4,2e18.7)
   21 format(2f10.4)
   22 format(1h )
      stop
      end 
      function ep(qx,qy)
      complex z
      a=sqrt(3.0)/2.0
      b=1.0/2.0
      z=0.0
      z=z+cexp(cmplx(0.0,qx*a+qy*b))
      z=z+cexp(cmplx(0.0,-qx*a+qy*b))
      z=z+cexp(cmplx(0.0,-qy))
      ep=cabs(z)      
      return
      end