Tidying up the input, output and code for DOS and pDOS

docs/post-proc.rst

@@ -190,21 +190,30 @@ following parameters can be set:
 
 ::
 
-   Process.min_DOS_E real    (Ha, default lowest eigenvalue)
-   Process.max_DOS_E real    (Ha, default highest eigenvalue)
-   Process.sigma_DOS real    (Ha, default 0.001)
-   Process.n_DOS     integer (default 1001)
+   Process.min_DOS_E       real    (Ha, default lowest eigenvalue)
+   Process.max_DOS_E       real    (Ha, default highest eigenvalue)
+   Process.WFRangeRelative T/F
+   Process.sigma_DOS       real    (Ha, default 0.001)
+   Process.n_DOS           integer (default 1001)
 
 The limits for the DOS are set by the first two parameters (note that
-CONQUEST will output all eigenvalues, so the limits on these are set
-by the eigenspectrum).  The broadening applied to each state is set by
+CONQUEST will output all eigenvalues).  These limits can be given as
+absolute energies, or relative to the Fermi level, depending on whether
+``WFRangeRelative`` is set to ``F`` or ``T`` respectively.
+The broadening applied to each state is set by
 ``sigma_DOS``, while the number of bins is set by ``n_DOS``.  The
 integrated DOS is also calculated; the user can choose whether this
 is the total integrated DOS (i.e. from the lowest eigenvalue,
 regardless of the lower limit for DOS) or just the local integrated
 DOS (i.e. over the interval specified for the DOS) by setting
 ``Process.TotalIntegratedDOS`` to ``T`` or ``F``, respectively.
 
+If the user does not specify limits on the DOS, then the range will
+be expanded slightly so that the broadened peaks fit completely within
+the DOS range.  This behaviour can be controlled with an additional
+parameter ``Process.ExpandRange T/F`` (the default behaviour is false
+if the user sets any limits, otherwise true).
+
 We recommend that, for accurate DOS, CONQUEST should be run
 non-self-consistently with a very high k-point density, after reading
 in a well-converged input charge density: set ``minE.SelfConsistent
@@ -292,8 +301,8 @@ also be specified:
    Process.max_DOS_E  0.35
    Process.WFRangeRelative T
 
-where the final tag sets the minimum and maximum values relative to
-the Fermi level.
+where the final tag specifies that the limits on the DOS are given
+relative to the Fermi level.
 
 If you only want to produce pDOS for a few atoms, then you can set
 the  variable ``Process.n_atoms_pDOS`` and list the atoms you want

src/control.f90

@@ -152,7 +152,7 @@ subroutine control_run(fixed_potential, vary_mu, total_energy)
     use force_module,         only: tot_force, stress
     use minimise,             only: get_E_and_F
     use global_module,        only: runtype, flag_self_consistent, &
-                                    flag_out_wf, flag_write_DOS, wf_self_con, &
+                                    flag_out_wf, flag_write_projected_DOS, wf_self_con, &
                                     flag_opt_cell, optcell_method, min_layer, flag_DM_converged
     use input_module,         only: leqi
     use store_matrix,         only: dump_pos_and_matrices
@@ -180,12 +180,10 @@ subroutine control_run(fixed_potential, vary_mu, total_energy)
 !****lat>$
     flag_DM_converged = .false.
     if ( leqi(runtype,'static') ) then
-       !if(.NOT.flag_self_consistent.AND.(flag_out_wf.OR.flag_write_DOS)) return
        flag_ff = .true.
        flag_wf = .true.
-       if (flag_out_wf.OR.flag_write_DOS) then
+       if (flag_out_wf.OR.flag_write_projected_DOS) then
           ! This is done within get_E_and_F
-          !wf_self_con=.true.
           flag_ff = .false.
           flag_wf = .false.
        endif

src/global_module.f90

@@ -397,7 +397,7 @@ module global_module
   integer :: mx_temp_matrices                   ! Defaults to 100; used in mult_module (immi)
 
   ! DOS output (NB Maybe move these into DiagModule and revisit names)
-  logical :: flag_write_DOS, flag_write_projected_DOS, flag_normalise_pDOS, flag_pDOS_angmom, flag_pDOS_lm
+  logical :: flag_write_projected_DOS, flag_normalise_pDOS, flag_pDOS_angmom, flag_pDOS_lm
   real(double) :: E_DOS_min, E_DOS_max, sigma_DOS
   integer :: n_DOS
 

src/initial_read_module.f90

@@ -865,7 +865,7 @@ subroutine read_input(start, start_L, titles, vary_mu,&
          flag_propagateL,flag_dissipation,integratorXL, flag_FixCOM,   &
          flag_exx, exx_alpha, exx_scf, exx_scf_tol, exx_siter, exx_cutoff, &
          flag_out_wf,max_wf,out_wf,wf_self_con, flag_fire_qMD, &
-         flag_write_DOS, flag_write_projected_DOS, &
+         flag_write_projected_DOS, &
          E_wf_min, E_wf_max, flag_wf_range_Ef, &
          mx_temp_matrices, flag_neutral_atom, flag_diagonalisation, &
          flag_SpinDependentSF, flag_Multisite, flag_LFD, flag_SFcoeffReuse, &
@@ -1784,25 +1784,17 @@ subroutine read_input(start, start_L, titles, vary_mu,&
           call cq_abort("Won't output WFs for Order(N) or non-static runs")
        end if
     end if
-    ! DOS output
-    flag_write_DOS = fdf_boolean('IO.writeDOS',.false.)
-    if(flag_write_DOS) then
+    ! pDOS output
+    flag_write_projected_DOS = fdf_boolean('IO.write_proj_DOS',.false.)
+    if(flag_write_projected_DOS) then
        if(flag_diagonalisation) then
-          flag_write_projected_DOS = fdf_boolean('IO.write_proj_DOS',.false.)
-          if(flag_write_projected_DOS) then
-             flag_out_wf = .true.
-             E_wf_min = fdf_double('IO.min_wf_E',-BIG)
-             E_wf_max = fdf_double('IO.max_wf_E',BIG)
-          end if
-          ! Possibly needed to decide if MSSF needs dealing with
-          !flag_pDOS_angmom = fdf_boolean('IO.PDOS_Angmom',.false.)
-          !if (flag_pDOS_angmom .and. flag_basis_set==blips) then
-          !   flag_pDOS_angmom = .false.
-          !   if(inode==ionode) write(io_lun,'(2x,"Setting IO.PDOS_Angmom F as using blips")')
-          !endif
+          flag_out_wf = .true.
+          E_wf_min = fdf_double('IO.min_wf_E',-BIG)
+          E_wf_max = fdf_double('IO.max_wf_E',BIG)
+          flag_wf_range_Ef = fdf_boolean('IO.WFRangeRelative',.true.)
        else
-          flag_write_DOS = .false.
-          if(inode==ionode) write(io_lun,'(2x,"Setting IO.writeDOS F as solving O(N)")')
+          flag_write_projected_DOS = .false.
+          if(inode==ionode) write(io_lun,'(2x,"Setting IO.write_proj_DOS F as solving O(N)")')
        end if
     end if
 !!$

src/initialisation_module.f90

@@ -1111,8 +1111,7 @@ subroutine initial_H(start, start_L, find_chdens, fixed_potential, &
          MDinit_step,ni_in_cell,            &
          flag_dissipation,     &
          flag_propagateX, flag_propagateL, restart_X, &
-         flag_out_wf, wf_self_con, &
-         flag_write_DOS, flag_neutral_atom, &
+         flag_neutral_atom, &
          atomf, sf, flag_LFD, nspin_SF, flag_diagonalisation, &
          ne_in_cell, min_layer, flag_basis_set, PAOs
     use ion_electrostatic,   only: ewald, screened_ion_interaction
@@ -1401,77 +1400,6 @@ subroutine initial_H(start, start_L, find_chdens, fixed_potential, &
        call get_H_matrix(rebuild_KE_NL, fixed_potential, electrons, &
             density, maxngrid, level=backtrace_level)
     endif
-!    if ( flag_self_consistent ) then ! Vary only DM and charge density
-!       !
-!       if ( restart_DM ) then
-!          record  = .true.
-!          reset_L = .false.
-!          call new_SC_potl(record, sc_tolerance, reset_L, &
-!               fixed_potential, vary_mu, n_L_iterations,  &
-!               L_tolerance, total_energy, backtrace_level)
-!          !
-!       else
-!          if (flag_LFD .and. .not.read_option) then
-!             ! Hpao was already made in sub:initial_SFcoeff
-!             rebuild_KE_NL = .false. 
-!             call get_H_matrix(rebuild_KE_NL, fixed_potential, electrons, &
-!                  density, maxngrid, level=backtrace_level, build_AtomF_matrix=.false.)
-!          else
-!             rebuild_KE_NL = .true. 
-!             call get_H_matrix(rebuild_KE_NL, fixed_potential, electrons, &
-!                  density, maxngrid, level=backtrace_level)
-!          endif
-!          !
-!          electrons_tot = spin_factor * sum(electrons)
-!          !
-!          record  = .false.
-!          reset_L = .true.                
-!          call FindMinDM(n_L_iterations, vary_mu, L_tolerance, &
-!               reset_L, record, backtrace_level)
-!          !
-!          record  = .true.             
-!          reset_L = .false.
-!          call new_SC_potl(record, sc_tolerance, reset_L, &
-!               fixed_potential, vary_mu, n_L_iterations,  &
-!               L_tolerance, total_energy, backtrace_level)
-!          !
-!       end if
-!       !
-!    else ! Ab initio TB: vary only DM
-!
-!       rebuild_KE_NL = .true.
-!       !build_X = .false
-!       if (flag_LFD .and. .not.read_option) then
-!          ! Hpao was already made in sub:initial_SFcoeff
-!          rebuild_KE_NL = .false.
-!          call get_H_matrix(rebuild_KE_NL, fixed_potential, electrons, &
-!               density, maxngrid, level=backtrace_level, build_AtomF_matrix=.false.)
-!       else
-!          rebuild_KE_NL = .true.
-!          call get_H_matrix(rebuild_KE_NL, fixed_potential, electrons, &
-!               density, maxngrid, level=backtrace_level)
-!       endif
-!       electrons_tot = spin_factor * sum(electrons)
-!       if (flag_out_wf.OR.flag_write_DOS) then
-!          wf_self_con=.true.
-!       endif
-!
-!       if ( .not. restart_DM ) then
-!          record  = .false.   
-!          reset_L = .true.
-!          call FindMinDM(n_L_iterations, vary_mu, L_tolerance, &
-!               reset_L, record, backtrace_level)
-!       else
-!          record  = .false.
-!          reset_L = .false.
-!          call FindMinDM(n_L_iterations, vary_mu, L_tolerance, &
-!               reset_L, record, backtrace_level)
-!       end if
-!       if (flag_out_wf.OR.flag_write_DOS) then
-!          wf_self_con=.false.
-!       endif
-!       call get_energy(total_energy=total_energy,level=backtrace_level)
-!    end if
 !!$
 !!$
 !!$

tools/PostProcessing/local_module.f90

@@ -56,4 +56,5 @@ module local
   integer :: flag_smear_type, iMethfessel_Paxton
   integer :: n_atoms_pDOS
   integer, dimension(:), allocatable :: pDOS_atom_index
+  logical :: flag_expand_range
 end module local