Ferromagnetism and Antiferromagnetism in 2D Hubbard model.

      program dos
c     density of states
      implicit none

c     This program will caculate the density of states for the half filled
c     2-D Hubbard model.
c     write Hubbard model...
      
      integer ne, n, i, j, k, im, ip


c     i,j,k,im,ip: auxiliary integers
c     ne: the energy bins
c     n: the mommentum bins
      parameter(ne=200)
      parameter(n=1000)

      real pi, d(ne), emax, de, dk, ek, ex, ey, kx, ky
      real em, ep, e0, a, b, c, dx, dy

c     dx,dy: the sides of the triangle in the bins
c     a,b,c: auxiliary real numbers
c     em,ep: min and max energies for each bin...
c     e0 the energy on lower left the corner...

      
c     kx,ky: the x and y components of the wave vector..
c     ek: the dispersion relation for energy
c     ex: it's derivative with respect to kx
c     ey: it's derivative with respect to ky
c     de: energy binning
c     dk: momentum binning
c     emax: the maximun energy we will go
c     d(ne): the density of states

c     pi=3.14...
      pi=acos(-1.)

c     defining maximum energy.
      emax=6.

c     establishing energy and momentum binning.
      de=emax/real(ne)
      dk=2*pi/real(n)

c     opening output file.
      open(1, file="output.dat")

c     setting d(i) zero.
      do i=1,ne
         d(i)=0.
      enddo

c     loop in mommenta.

      do i=1,n
         do j=1,n

            kx= (real(i)-0.5)*dk
            ky= (real(j)-0.5)*dk
c     E(k)=-2t(cos(kx)+cos(ky))
            ek=-2.0*(cos(kx)+cos(ky))
            ex=2.0*abs(sin(kx))
            ey=2.0*abs(sin(ky))
            
c     displace energy to the lower left corner.
            a=(ex+ey)*dk/2.0
            e0=ek-a
            im=int(e0/de)-4
            ip=im+2*int(a/de)+8

            if (im .lt. 1) then
               im=1
               
            elseif (im .gt. ne) then
               im=ne

            endif

            if (ip .lt. 1) then
               ip=1

            elseif (im .gt. ne) then
               ip=ne

            endif

c     loop over the bins
            do k=im,ip
               
               em=real(k-1)*de
               ep=real(k)*de

c     compute the fist area (b)
               dx=(em-e0)/ex
               dy=(em-e0)/ey
         
               b=0.0
               
               if (dx .lt. 0.0 .or. dy .lt. 0.0) then
                  goto 2

               endif

                  b=dx*dy
                  
               if(dx .gt. dk) then
                  b=b-(dx-dk)**2*dy/dx

               elseif(dy .gt. dk) then
                  b=b-(dy-dk)**2*dx/dy
                  
               elseif(dx*dy .lt. dk*(dx+dy)) then
                  b=2.0*dk**2
                     
               endif

 2             continue

c     compute the second area (c)
               dx=(ep-e0)/ex
               dy=(ep-e0)/ey
         
               c=0.0
               
               if (dx .lt. 0.0 .or. dy .lt. 0.0) then
                  goto 3

               endif

                  c=dx*dy
                  
               if(dx .gt. dk) then
                  c=c-(dx-dk)**2*dy/dx

               elseif(dy .gt. dk) then
                  c=c-(dy-dk)**2*dx/dy
                  
               elseif(dx*dy .lt. dk*(dx+dy)) then
                  c=2.0*dk**2

               endif

 3             continue

c     increment energy bins

               d(k)=d(k)+abs(b-c)

            enddo
   

         enddo
      enddo

c     normalize the result.

      a=0.5/(2*pi)**2
      c=0.
      
      do i=1,ne
         d(i)=a*d(i)
      enddo

c     write the output.

      do i=1,ne
         
         b=(real(i)-0.5)*de
         c=c+d(i)

         write(1,20) b, d(i)

      enddo
      print*, b

 20   format(f10.4, e18.7)

      stop
      end