program qdos3

c

c     Performs a "quick" calculation of the electronic density

c      of states (DOS) of a first-principles or tight-binding

c      calculation, using the Fermi broadening term so that the

c      "smeared" DOS is that same as that used to calculate the

c      total energy.

c

c     The eigenvalues, k-points, and weights are read from a QLMT-like

c      file (QLMT files are produced by the LAPW and tight-binding

c      "static" codes), with a user-supplied Fermi Temperature.

c      The output can be directed to any file.  At the moment we only

c      calculate the total DOS, though the angular momentum DOS are

c      rather easy to calculate with this scheme.

c

c     We'll also assume that this is a spin-restricted calculation,

c      i.e., two electrons per state

c

c     The only other assumptions we make about the QLMT files are that

c      it represents only one window and that the eigenvalues are

c      ordered from lowest to highest within each k-point.  This

c      will speed some of the calculations, as we won't have to

c      consider large eigenvalues.

c

c                                     --  mjm  14 October 1997

c

c     Since it costs next to nothing to calculation rho'(mu), we'll add

c      that in.                       --  mjm  15 October 1997

c

c     The number of atoms is in the QLMT file, so use that, don't

c      ask for it anymore.            --  mjm  11 February 1999

c

c     In this version, the code is rewritten as a filter.  The QLMT

c      file is read from standard input, the DOS is written on

c      standard output, and the command line is used to transmit

c      the Fermi temperature, starting and stopping ranges, and

c      and options (see the POD section below for more details).

c                                     --  mjm  20 August 2001

c=======================================================================

cpod

cpod 

cpod This is a man page written in Perl's Plain Old Documentation (POD)

cpod format.  To construct a man page, from this documentation, run the

cpod following commands:

cpod

cpod $ grep '^cpod' qdos3.f | sed "s/^cpod//" > /tmp/qdos3.perlman

cpod $ pod2man /tmp/qdos3.perlman > qdos3.1

cpod $ rm /tmp/qdos3.perlman

cpod 

cpod qdos3.1 is in standard manual page format, and can be placed in

cpod any convenient man/man1 subdirectory.  To read the man page, use

cpod the command

cpod

cpod $ nroff -man qdos3.1 | less -imeX

cpod

cpod Users of non-Unix-like systems should consult their OS manuals to

cpod determine the appropriate commands equivalent to the above.

cpod

cpod

cpod=head1 NAME

cpod

cpodI<qdos3> -- Quick and dirty Electronic Density of States

cpod

cpod=head1 SYNOPSIS

cpod

cpodI<qdos3> [B<-h>] [B<-r>] [B<-v>] Temperature E(start) E(stop) E(delta) < QLMT > dos

cpod

cpod=head1 DESCRIPTION

cpod

cpodA filter which determines the quick and dirty approximation to the

cpodelectronic Density of States (DOS).  Each eigenvalue's contribution

cpodto the DOS at a given energy is determined by using Fermi

cpodbroadening.

cpod

cpodIf the atomic and angular momentum decomposition of the state is

cpodgiven, the output will include a similar decomposition at every

cpodenergy level.

cpod

cpodNote that this is a filter.  Diagnostic messages are written to

cpodSTDERR if the B<-v> option is used.

cpod

cpod=head1 ARGUMENTS

cpod

cpodAll arguments must have the same units of energy as is used in the

cpodQLMT file.

cpod

cpodB<Temperature> is that used to evaluate the Fermi function.  It may

cpodor not have anything to do with a physical temperature.

cpod

cpodB<E(start)> is the first energy to be plotted.

cpod

cpodB<E(stop)> is I<near> the last energy to be plotted.  The exact

cpodvalue of the last energy will depend on E(start), E(delta), and the

cpodway the computer rounds.

cpod

cpod=head1 OPTIONS

cpod

cpodB<-h> Print a brief help message.

cpod

cpodB<-r> Determine the smallest and largest eigenvalue in the file,

cpodprinting them on one line on standard output.  In this case it is

cpodnot necessary to define the temperature or energy range.

cpod

cpodNote that B<-r> suspends reading of the arguments.  Thus if you

cpodwant the program to be verbose you must specify B<-v> F<before> B<-r>.

cpod

cpodB<-v> Be verbal -- print information about the program on standard

cpoderror.

cpod

cpod=head1 NO WARRANTY

cpod

cpodThis software and any accompanying documentation are released "as

cpodis." The U.S. Government makes no warranty of any kind, expressed

cpodor implied, concerning this software and any accompanying

cpoddocumentation, including without limitation, any warranties of

cpodmerchantability or fitness for a particular purpose. In no event

cpodwill the U.S. Government be liable for any damages, including any

cpodlost profits, lost savings, or other incidental or consequential

cpoddamages arising out of the use, or inability of use, of this

cpodsoftware or any accompanying documentation, even if informed in

cpodadvance of the possibility of such damages.

cpod

cpod=head1 AUTHOR

cpod

cpod    Michael Mehl

cpod    Center for Computational Materials Science

cpod    Code 6393

cpod    Naval Research Laboratory

cpod    Washington DC 20375-5000

cpod

cpod    email: mehl@dave.nrl.navy.mil

cpod

cpod=cut

cpod

cpod

c=======================================================================

c

c     Theory:

c

c     The following theory is written for a Fermi-style smearing

c      function, where the Fermi temperature is T and b = 1/T:

c

c     f(x) = 1/(1+exp(b*x))

c

c     The formalism can be easily generalized to any appropriate

c      f(x).  Note that for the Fermi function,

c

c     f'(x) = -b exp(b*x)/(1+exp(b*x))^2 = -b f(x) f(-x)

c           = -b exp(b*x) f(x)^2

c

c     We assume that the k-points in the QLMT file give an adequate

c      sampling of the Brillouin zone for the given temperature T.

c      Let w(k) be the weight associated with the k-th k-point, and

c      eig(i,k) be the ith eigenvalue at the point.  Further assume

c      that there are nk k-points and neig(k) eigenvalues at the

c      k-th k-point.

c

c     Then occupation number for a given fermi energy mu is then

c

c     N(mu) = Sum[w(k) Sum[f(eig(i,k)-mu),{i,1,neig(k)}],{k,1,nk}]

c

c     Since

c

c     N(mu) = Integral[rho(e),{e,-oo,mu}],

c

c     where rho(e) is the electronic density of states at the

c      energy e.  Thus

c

c     rho(e) = N'(mu)

c            = -Sum[w(k) Sum[f'(eig(i,k)-mu),{i,1,neig(k)}],{k,1,nk}]

c

c     For the Fermi distribution we have

c

c     rho(e) = b*Sum[w(k) Sum[f(eig(i,k)-mu)*f(mu-eig(i,k)),

c                                        {i,1,neig(k)}],{k,1,nk}]

c

c     From the above, calculating rho'(e) is rather trivial.  Note

c      that f''(x) = f'(x) [ 1 - 2 e^x f(x)] and we're done.

c

c     Note that we can also calculate the band energy up to mu:

c

c     E(mu) = Integral[e*rho(e),{e,-oo,mu}]

c

c=======================================================================

c

c     Start with a clean slate:

c

      implicit none

c

c=======================================================================

c

c     maxk is the maximum number of k-points in the calculation

c     maxeig is the maximum number of eigenvalues at a given k-point

c     matom is the maximum number of atom types in the QLMT file

c

      integer maxk,maxeig,matom

c

      parameter (maxk = 3000)

      parameter (matom = 20)

      parameter (maxeig = 9*matom)

c

c=======================================================================

c

c     Real scalars and arrays

c

c     farg(i) is the ith numeric argument on the command line.  These

c      will be assigned in the order:

c

c     t     = farg(1)

c     mumin = farg(2)

c     mumax = farg(3)

c     mudel = farg(4)

c

      real*8 farg(4)

c

c     t is the Fermi temperature

c     b = 1/t  -- Note that t=0 will produce a generic large beta

c     mu is the current chemical potential

c     mumin, mumax, mudel are the minimum, maximum, and step size

c      in mu.

c     t20 = 20*t

c     temp is a temporary holding variable

c

      real*8 t,b,mu,mumin,mumax,mudel,t20,temp

c

c     wtsum is the sum of the k-point weights

c     nsum will be the total occupation of the QLMT file.  That is,

c      lim (e -> oo) rho(e) = nsum

c     esum will be the total energy summation

c     elow is the lowest eigenvalue in the file

c     ehigh is the highest eigenvalue in the file

c     elocal is an energy summation variable over a single k-point

c

      real*8 wtsum,nsum,esum,elow,ehigh,elocal

c

c     w(k) is the weight at the k-th k-point

c     eig(j,k) is the jth eigenvalue at the k-th k-point

c

      real*8 w(maxk),eig(maxeig,maxk)

c

c     lowe = mu - 20*t is the lowest energy where f(e-mu) is

c      appreciably different from 1

c     highe = mu + 20*t is the highest energy where f(e-mu) is

c      appreciably different from 0

c

      real*8 lowe,highe

c

c     nstart(k) is the total weight of electrons at energies

c      less than lowe for the k-th kpoint

c

c     estart(k) is the sum of all energies < lowe for the k-th

c      point

c

      real*8 nstart(maxk),estart(maxk)

c

c     n = N(mu)

c     rho = rho(mu)

c     rhop = rho'(mu)

c     e = e(mu)

c     nlocal, rlocal, and rplocal are similar to elocal

c

      real*8 n,rho,rhop,e,nlocal,rlocal,rplocal

c

c     f is the local value of f(e-mu), fp = - f'(e-mu), and

c     fpp = f''(e-mu)

c

      real*8 f,fp,fpp

c

c

c     awght(l,i,j,k) is the weight inside the muffin tin of

c      of orbitals with angular momentum l-1 on atom i, at the

c      jth state at the kth k-kpoint.

c

c     pdos(l,i) is the partial dos of angular momentum l-1

c      on atom i.  pdosl(l,i) is the local dos for the current

c      state.

c

      real*8 awght(4,matom,maxeig,maxk),pdos(4,matom),pdosl(4,matom)

c

c=======================================================================

c

c     Integers

c

c     natom is the number of atoms represented in the QLMT file

c      this may be zero.

c     nk is the number of k-points read in

c

      integer natom,nk

c

c     neig(k) is the number of eigenvalues at the k-th k-point

c

      integer neig(maxk)

c

c     istart(k) is the first-eigenvalue we examine at the k-th

c      k-point.  This assumes that all eigenvalues below istart

c      are fully occupied.

c     The loop actually starts at ik, since istart(k) may

c      change during execution

c

      integer istart(maxk),ik

c

c     nargs is the actual number of arguments in the argument list:

c

      integer nargs

c

c     Do loop and other counting variables

c

      integer i,j,k,l

c

c=======================================================================

c

c     Logical variables:

c

c     frange is true if the "-r" option is invoked, i.e., we want to

c      find the smallest and largest eigenvalues.  

c

      logical frange

c

c     if verbose is true ("-v") we print extra information and

c      diagnostics to standard error.  Otherwise the program is

c      relatively silent, except for the DOS, which goes to standard

c      output.

c

      logical verbose

c

c=======================================================================

c

c     Strings and Things

c

c     fout holds the output string format

c

      character*11 fout

c

c     Defined flags (actual definitions in character statements below):

c

      character*2 hflag,rflag,vflag

c

c     arg(i) is a string containing the ith command line argument

c

      character*80 arg(10)

c

c=======================================================================

c

c     Real numbers we'd like to keep

c

      real*8 zero,one,two,twenty

c

      parameter (zero = 0d0)

      parameter (one = 1d0)

      parameter (two = 2d0)

      parameter (twenty = 2d1)

c

      real*8 tol,small,big

c

      parameter (tol = 1d-10)

      parameter (small = 1d-10)

      parameter (big = 1d20)

c

c=======================================================================

c

      data hflag/"-h"/

      data rflag/"-r"/

      data vflag/"-v"/

c

c=======================================================================

c

c     Start execution

c

c     Default behavior is a quiet program which prints the full DOS

c

      verbose = .false.

      frange = .false.

c

c     How many arguments?

c

      call args(nargs,arg)

c

c     Scroll through the arguments.  Use "j" as the flag for the

c      numeric arguments.

c

      j = 1

c

      do i = 1,nargs

c

c         Check for verbosity:

c

         if (arg(i).eq.vflag) then

            verbose = .true.

            write(0,*) 'Writing diagnostics to standard error'

c

c        Check for a "usage" request

c

         else if (arg(i).eq.hflag) then

            write(0,1013)

 1013       format ('qdos [-h] [-r] [-v] Temperature E(start) ',

     $           'E(stop) E(delta) < QLMT > dos')

            call exit

c

c        If "rflag" is defined, then we're only looking for the

c         minimum and maximum eigenvalues:

c

         else if (arg(i).eq.rflag) then

            frange = .true.

c

c           Make the program happy by giving it a fake temperature

c            and energy range:

c

            j = 5

            farg(1) = 5d-2

            farg(2) = -1d1

            farg(3) = 1d1

            farg(4) = 1d0

c

c           Forgive us, for we are about to goto:

c

            go to 100

         else

c

c           The argument is assumed to be a number

c

            if(j.lt.5) then

               read(arg(i),'(f80.70)') farg(j)

               j = j + 1

            end if

         end if

      end do

c

c     If j isn't 5, then we don't have enough arguments:

c

      if(j.lt.5) then

         write(0,1013)

         call exit(1)

      end if

c

c     The only goto to this statement is from the "rflag"

c      section of the arguments, see above

c

 100  continue

c

c     Assign the arguments as needed:

c

      t = farg(1)

      mumin = farg(2)

      mumax = farg(3) + 1d-2*abs(farg(4))

      mudel = farg(4)

c

      if (verbose) then

         write(0,*) 'Fermi temperature = ',t

         write(0,*) 'DOS printed from ',mumin,' to ',mumax

         write(0,*) 'with step size ',mudel

      end if

c

c

c     Make a default size for t, then define b

c

      t = max(t,small)

      t20 = twenty*t

      b = one/t

c

c     Get the range of the calculation

c     make sure mumin < mumax:

c

      mudel = abs(mudel)

      if(mumin.gt.mumax) then

         temp = mumax

         mumax = mumin

         mumin = temp

      end if

c

c     Now read in the information in the QLMT file.  Which is

c      presumed to be on standard input.  We also assume

c      only one window, so the only thing worth reading in the

c      first line is the number of atoms.  Note that there

c      are 3 dummy integer variables ahead of it:

c

      read(*,'(15x,i5)') natom

c

c     Check atom size:

c

      if(natom.gt.matom) then

         write(0,*) 'QLMT file has',natom,' atoms'

         write(0,*) 'Program dimensioned for ',matom,' atoms'

         write(0,*) 'Redimension matom in qdos.f'

         call exit(1)

      end if

c

c     nk is the number of k-points we'll need:

c

      read(*,*) nk

c

c     Make sure we can handle the number of k-points:

c

      if(nk.gt.maxk) then

         write(0,*) 'QLMT file has',nk,' k-points'

         write(0,*) 'Program dimensioned for ',maxk,' k-points'

         write(0,*) 'Redimension maxk in qdos.f'

         call exit(2)

      end if

c

      if (verbose) write(0,'(''Reading '',i5,'' k-points'')') nk

c

c     Now read the eigenvalues at each k-point.

c

c     Keep some auxillary terms:

c

c     wtsum is the sum of the k-point weights.  It should sum to

c      1, but we can renormalize

c

      wtsum = zero

c

c     nsum will be the total occupation of the QLMT file.  That is,

c      lim (e -> oo) rho(e) = nsum

c

      nsum = zero

c

c     Similarly, esum will be the total energy summation

c

      esum = zero

c

c     elow is the lowest eigenvalue in the file

c     ehigh is the highest eigenvalue in the file

c

      elow = big

      ehigh = -big

c

      do k = 1,nk

         read(*,'(40x,i10,f20.10)') neig(k),w(k)

         wtsum = wtsum + w(k)

         nsum = nsum + dble(neig(k))*w(k)

c

         elocal = zero

c

         do j = 1,neig(k)

            read(*,*) eig(j,k)

            elow = min(elow,eig(j,k))

            ehigh = max(ehigh,eig(j,k))

            elocal = elocal + eig(j,k)

c

c           Read the partial occupancies for each atom:

c

            do i = 1,natom

               read(*,'(4e15.6)') (awght(l,i,j,k),l=1,4)

ctemp

c$$$               write(0,'(3i5,4f10.6)') i,j,k,(awght(l,i,j,k),l=1,4)

cend temp

            end do

         end do

         esum = w(k)*elocal + esum

      end do

c

c     If the "-r" option is invoked, print out and stop here:

c

      if (frange) then

         write(*,*) 'Minimum eigenvalue = ',elow,

     $        ' ; Maximum eigenvalue = ',ehigh

c

c        Standard exit

c

         call exit

      end if

ctemp

c$$$      write(0,*) (k,w(k),(eig(j,k),j=1,neig(k)),k=1,nk)

cend temp

c

c     We really want elow and ehigh adjusted by the 20*t rule:

c

      elow = elow - t20

      ehigh = ehigh + t20

c

c     Check the weight:

c

      if(abs(wtsum-one).gt.tol) then

         if (verbose) write(0,'(/''WARNING:  Found sum of weights = '',

     $        f20.10/''This will be renormalized to 1''/)') wtsum

         temp = one/wtsum

         do k = 1,nk

            w(k) = w(k)*temp

         end do

         nsum = nsum*temp

         esum = esum*temp

c$$$         wtsum = one

      end if

c

c     This is spin-restricted, and the simplest way to handle it

c      is to multiply all weights by two

c

      do k = 1,nk

         w(k) = two*w(k)

      end do

      nsum = two*nsum

      esum = two*esum

c

      if (verbose) write(0,'(/''Found '',f10.5,

     $     '' electrons with e = '',f15.5)') nsum,esum

c

c     We get quite a bit of speedup if we assume that mumin < mumax

c      Then once an eigenvalue is below lowe for one value of mu,

c      it is less than lowe for all succeding values of mu

c

c     Thus all k-points below istart(k) are fully occupied,

c      contributing nstart(k) electrons to the calculation,

c      with energy estart(k)

c

c     Enter the starting values of these variables

c

      do k = 1,nk

         istart(k) = 1

         nstart(k) = zero

         estart(k) = zero

      end do

c

c     Start going through the energies

c

c     The output will be

c

c     mu,n(mu),rho(mu),rho'(mu),e(mu)

c

c     Set up output format

c

      write(fout,'(''('',i4.4,''f12.5)'')') 4*natom+3

c$$$      write(0,'(a11)') fout

c

      do mu = mumin,mumax,mudel

c

c        Some calculations are easy

c

         if(mu.lt.elow) then

            write(*,fout) mu,zero,zero,((zero,l=1,4),j=1,natom)

         else if(mu.gt.ehigh) then

            write(*,fout) mu,nsum,zero,((zero,l=1,4),j=1,natom)

         else

            lowe = mu - t20

            highe = mu + t20

c

c           Here we go

c

            n = zero

            rho = zero

            rhop = zero

            e = zero

            do i = 1,natom

               do l = 1,4

                  pdos(l,i) = zero

               end do

            end do

c

            do k = 1,nk

c

c              Start with the stuff currently below istart(k)

c

               nlocal = nstart(k)

               elocal = estart(k)

               rlocal = zero

               rplocal = zero

               do j = 1,natom

                  do l = 1,4

                     pdosl(l,j) = zero

                  end do

               end do

c

c              Since we may want to change istart(k) below,

c               start at ik instead

c

               ik = istart(k)

c

c              Note that for large mu this loop may not execute

c

               do i = ik,neig(k)

c

c                 Is this eigenvalue now to small?  Then

c                  add it to the starting values

c

                  if(eig(i,k).lt.lowe) then

                     nlocal = nlocal + one

                     elocal = elocal + eig(i,k)

c

c                    Update starting values

c

                     istart(k) = istart(k) + 1

                     nstart(k) = nlocal

                     estart(k) = elocal

                  else if(eig(i,k).gt.highe) then

c

c                    We're done for this k-point.  Get out:

c

                     go to 1000

                  else

c

c                    Evaluate f and its derivative

c

                     temp = exp(b*(eig(i,k)-mu))

                     f = one/(one+temp)

c

c                    Note that we store -f' rather than f',

c                     and that we're leaving out the multiplications

c                     by b and b^2 until the end

c

                     fp = (temp*f)*f

c

c                    Thus we have to change the sign on fpp:

c

                     fpp = fp*(two*(temp*f)-one)

c

c                    Increment the counts

c

                     nlocal = nlocal + f

                     elocal = elocal + eig(i,k)*f

                     rlocal = rlocal + fp

                     rplocal = rplocal + fpp

                     do j = 1,natom

                        do l = 1,4

                           pdosl(l,j) = pdosl(l,j)+fp*awght(l,j,i,k)

                        end do

                     end do

                  end if

               end do

c

c              Now update according to the weight of this k-point

c

 1000          n = n + w(k)*nlocal

               e = e + w(k)*elocal

               rho = rho + w(k)*rlocal

               rhop = rhop + w(k)*rplocal

               do j = 1,natom

                  do l = 1,4

                     pdos(l,j) = pdos(l,j) + w(k)*pdosl(l,j)

                  end do

               end do

            end do

c

c           Write it out, including the factors of b and b^2:

c

            write(*,fout) mu,n,b*rho,((b*pdos(l,j),l=1,4),j=1,natom)

ctemp

c$$$            write(0,105) mu,n,b*rho,b2*rhop,e

cend temp

         end if

      end do

      end

c_______________________________________________________________________

c

      subroutine  args (nargs,arg)

c     Reads unix command line and parses it into arguments

C     Should work on any Unix computer

      character*(*) arg(1)

c     iargc reads the number of arguements

      nargs=iargc()

      do i = 1,nargs

c        getarg reads an individual argument

         call getarg(i,arg(i))

      end do

      return

      end