diff --git a/docs/post-proc.rst b/docs/post-proc.rst index 70e2e0d86..0770e0006 100644 --- a/docs/post-proc.rst +++ b/docs/post-proc.rst @@ -190,14 +190,17 @@ 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, @@ -205,6 +208,12 @@ 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 diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index fc5e02291..a8af360fb 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -3797,288 +3797,6 @@ subroutine buildK(range, matA, occs, kps, weight, localEig, overlapEig, matSij, end subroutine buildK !!*** - !!****f* DiagModule/accumulate_DOS - !! - !! NAME - !! accumulate_DOS - !! USAGE - !! - !! PURPOSE - !! Accumulates DOS - !! INPUTS - !! - !! USES - !! - !! AUTHOR - !! D. R. Bowler - !! CREATION DATE - !! 2016 ? - !! MODIFICATION HISTORY - !! 2017/11/01 18:00 nakata - !! Introduced PDOS with MSSFs, projecting on neighbor atoms with global ID. - !! Added spinSF to specify the spin of the MSSF coefficients. - !! SpinSF is not used for primitive PAOs. - !! 2017/11/13 18:15 nakata - !! Introduced the normalization of each eigenstate of PDOS. - !! Added optional argument weight_pDOS, the normalisation weight. - !! 2018/09/19 18:30 nakata - !! Introduced orbital angular momentum resolved DOS. - !! Added optional projDOS_angmom, l-projected PDOS. - !! 2018/10/22 14:18 dave & jsb - !! Adding (l,m) projection for pDOS - !! 2018/10/30 11:43 dave - !! Implementing semi-core exclusion given right flags - !! (NB at present, semi-core states are not flagged but will be !) - !! 2018/11/02 16:30 nakata - !! Bug fix: changed atom_spec to neigh_species for semicore of neighbour atoms - !! SOURCE - !! -!**! subroutine accumulate_DOS(weight,eval,evec,DOS,spinSF,projDOS,projDOS_angmom,weight_pDOS) -!**! -!**! use datatypes -!**! use numbers, only: half, zero -!**! use global_module, only: n_DOS, E_DOS_max, E_DOS_min, flag_write_DOS, sigma_DOS, flag_write_projected_DOS, & -!**! sf, atomf, id_glob, species_glob, flag_normalise_pDOS, flag_pDOS_angmom, flag_pDOS_lm, & -!**! flag_SpinDependentSF -!**! use ScalapackFormat, only: matrix_size -!**! use species_module, only: nsf_species, natomf_species -!**! use group_module, only: parts -!**! use primary_module, only: bundle -!**! use cover_module, only: BCS_parts -!**! use matrix_data, only: mat, halo, SFcoeff_range -!**! use mult_module, only: matSFcoeff, matrix_pos, mat_p -!**! use GenComms, only: cq_abort -!**! use pao_format -!**! -!**! implicit none -!**! -!**! ! Passed variables -!**! real(double) :: weight -!**! complex(double_cplx), dimension(:,:), intent(in) :: evec -!**! real(double), dimension(:) :: eval -!**! real(double), dimension(n_DOS) :: DOS -!**! integer, intent(in) :: spinSF ! used only if atomf/=sf -!**! real(double), OPTIONAL, dimension(:,:) :: projDOS -!**! real(double), OPTIONAL, dimension(:,:,:,:) :: projDOS_angmom ! Dimensions are bin, atom, l, m -!**! real(double), OPTIONAL, dimension(:) :: weight_pDOS -!**! -!**! ! Local variables -!**! integer :: iwf, n_band, n_min, n_max, i, acc, spin_SF, & -!**! atom, isf1, nsf1, atom_spec, l1, nacz1, m1, & -!**! neigh_global_num, iatomf2, natomf2, neigh_species, l2, nacz2, m2, & -!**! atom_num, gcspart, neigh_global_part, j_in_halo, wheremat -!**! integer :: iprim, part, memb, neigh, ist -!**! real(double) :: Ebin, a, fac, fac1, fac2, val -!**! real(double), dimension(6,13) :: fac_angmom, fac2_angmom ! up to h-orbital and 2l+1 -!**! real(double), dimension(n_DOS) :: tmp -!**! -!**! if(present(projDOS).AND.(.NOT.flag_write_projected_DOS)) call cq_abort("Called pDOS without flag") -!**! if(present(projDOS_angmom).AND.(.NOT.flag_pDOS_angmom)) call cq_abort("Called pDOS_angmom without flag") -!**! if(present(weight_pDOS) .AND.(.NOT.flag_normalise_pDOS)) call cq_abort("Normalised pDOS without flag") -!**! ! --------------- -!**! ! DOS calculation -!**! ! --------------- -!**! ! Now accumulate DOS for this band -!**! do iwf=1,matrix_size ! Effectively all bands -!**! tmp = zero -!**! n_band = floor((eval(iwf) - E_DOS_min)/dE_DOS) + 1 -!**! n_min = n_band - n_DOS_wid -!**! if(n_min<1) n_min = 1 -!**! n_max = n_band + n_DOS_wid -!**! if(n_max>n_DOS) n_max = n_DOS -!**! do i = n_min, n_max -!**! Ebin = real(i-1,double)*dE_DOS + E_DOS_min -!**! a = (Ebin-eval(iwf))/sigma_DOS -!**! tmp(i) = weight*pf_DOS*exp(-half*a*a) -!**! DOS(i) = DOS(i) + tmp(i) -!**! end do -!**! -!**! ! Having found DOS, we now project onto atoms -!**! if(flag_write_projected_DOS) then -!**! if (atomf == sf) then -!**! acc = 0 -!**! do atom=1,bundle%n_prim -!**! atom_spec = bundle%species(atom) -!**! fac = zero -!**! if (.not.flag_pDOS_angmom) then -!**! do isf1 = 1,nsf_species(atom_spec) -!**! fac = fac + real(evec(iwf,acc+isf1)*conjg(evec(iwf,acc+isf1)),double) -!**! end do -!**! else if(flag_pDOS_lm) then -!**! fac_angmom(:,:) = zero -!**! isf1 = 0 -!**! do l1 = 0, pao(atom_spec)%greatest_angmom -!**! do nacz1 = 1, pao(atom_spec)%angmom(l1)%n_zeta_in_angmom -!**! if((pao(atom_spec)%angmom(l1)%semicore(nacz1)==0) .OR. & -!**! (flag_pDOS_include_semicore)) then -!**! do m1 = -l1,l1 -!**! isf1 = isf1 + 1 -!**! fac = fac + real(evec(iwf,acc+isf1)*conjg(evec(iwf,acc+isf1)),double) -!**! ! l, m so shift m1 by l1+1 so it runs from 1 to 2*l1+1 -!**! fac_angmom(l1+1,m1+l1+1) = fac_angmom(l1+1,m1+l1+1) + & -!**! real(evec(iwf,acc+isf1)*conjg(evec(iwf,acc+isf1)),double) -!**! enddo -!**! end if -!**! enddo -!**! enddo -!**! else -!**! fac_angmom(:,:) = zero -!**! isf1 = 0 -!**! do l1 = 0, pao(atom_spec)%greatest_angmom -!**! do nacz1 = 1, pao(atom_spec)%angmom(l1)%n_zeta_in_angmom -!**! if((pao(atom_spec)%angmom(l1)%semicore(nacz1)==0) .OR. & -!**! (flag_pDOS_include_semicore)) then -!**! do m1 = -l1,l1 -!**! isf1 = isf1 + 1 -!**! fac = fac + real(evec(iwf,acc+isf1)*conjg(evec(iwf,acc+isf1)),double) -!**! fac_angmom(l1+1,1) = fac_angmom(l1+1,1) & -!**! + real(evec(iwf,acc+isf1)*conjg(evec(iwf,acc+isf1)),double) -!**! enddo -!**! end if -!**! enddo -!**! enddo -!**! if (isf1.ne.nsf_species(atom_spec)) call cq_abort("Error in NSF in the PDOS calculation.") -!**! endif -!**! if (flag_normalise_pDOS) then -!**! fac = fac / weight_pDOS(iwf) -!**! if (flag_pDOS_angmom) fac_angmom(:,:) = fac_angmom(:,:) / weight_pDOS(iwf) -!**! endif -!**! do i=n_min,n_max -!**! projDOS(i,atom) = projDOS(i,atom) + tmp(i)*fac -!**! end do -!**! if (flag_pDOS_angmom) then -!**! if(flag_pDOS_lm) then -!**! do l1 = 0, pao(atom_spec)%greatest_angmom -!**! do m1=-l1,l1 -!**! do i=n_min,n_max -!**! projDOS_angmom(i,atom,l1+1,m1+l1+1) = projDOS_angmom(i,atom,l1+1,m1+l1+1) + & -!**! tmp(i)*fac_angmom(l1+1,m1+l1+1) -!**! end do -!**! end do -!**! end do -!**! else -!**! do l1 = 0, pao(atom_spec)%greatest_angmom -!**! do i=n_min,n_max -!**! projDOS_angmom(i,atom,l1+1,1) = projDOS_angmom(i,atom,l1+1,1) + tmp(i)*fac_angmom(l1+1,1) -!**! end do -!**! end do -!**! end if -!**! endif -!**! acc = acc + nsf_species(atom_spec) -!**! end do ! atom -!**! else -!**! spin_SF = 1 -!**! if (flag_SpinDependentSF) spin_SF = spinSF -!**! acc = 0 -!**! iprim = 0 -!**! do part = 1,bundle%groups_on_node ! Loop over primary set partitions -!**! if(bundle%nm_nodgroup(part)>0) then ! If there are atoms in partition -!**! do memb = 1,bundle%nm_nodgroup(part) ! Loop over primary atoms -!**! atom_num = bundle%nm_nodbeg(part)+memb-1 -!**! iprim=iprim+1 -!**! nsf1 = nsf_species(bundle%species(atom_num)) ! = mat(part,SFcoeff_range)%ndimi(memb) -!**! do neigh = 1, mat(part,SFcoeff_range)%n_nab(memb) ! Loop over neighbours of atom -!**! fac = zero -!**! if (flag_pDOS_angmom) fac_angmom(:,:) = zero -!**! ist = mat(part,SFcoeff_range)%i_acc(memb)+neigh-1 -!**! gcspart = BCS_parts%icover_ibeg(mat(part,SFcoeff_range)%i_part(ist))+ & -!**! mat(part,SFcoeff_range)%i_seq(ist)-1 -!**! neigh_global_part = BCS_parts%lab_cell(mat(part,SFcoeff_range)%i_part(ist)) -!**! neigh_global_num = id_glob(parts%icell_beg(neigh_global_part)+ & -!**! mat(part,SFcoeff_range)%i_seq(ist)-1) -!**! neigh_species = species_glob(neigh_global_num) -!**! j_in_halo = halo(SFcoeff_range)%i_halo(gcspart) -!**! natomf2 = natomf_species(neigh_species) -!**! ! Now loop over support functions and atomf (basically PAOs) -!**! do isf1 = 1, nsf1 -!**! fac1 = real(evec(iwf,acc+isf1)*conjg(evec(iwf,acc+isf1)),double) -!**! fac2 = zero -!**! if (.not.flag_pDOS_angmom) then -!**! do iatomf2 = 1, natomf2 -!**! wheremat = matrix_pos(matSFcoeff(spin_SF),iprim,j_in_halo,isf1,iatomf2) -!**! val = mat_p(matSFcoeff(spin_SF))%matrix(wheremat) -!**! fac2 = fac2 + val*val -!**! enddo -!**! fac = fac + fac1 * fac2 -!**! else if(flag_pDOS_lm) then ! m and l resolved -!**! fac2_angmom(:,:) = zero -!**! iatomf2 = 0 -!**! do l2 = 0, pao(neigh_species)%greatest_angmom -!**! do nacz2 = 1, pao(neigh_species)%angmom(l2)%n_zeta_in_angmom -!**! if((pao(neigh_species)%angmom(l2)%semicore(nacz2)==0) .OR. & -!**! (flag_pDOS_include_semicore)) then -!**! do m2 = -l2,l2 -!**! iatomf2 = iatomf2 + 1 -!**! wheremat = matrix_pos(matSFcoeff(spin_SF),iprim,j_in_halo,isf1,iatomf2) -!**! val = mat_p(matSFcoeff(spin_SF))%matrix(wheremat) -!**! fac2 = fac2 + val*val -!**! fac2_angmom(l2+1,m2+l2+1) = fac2_angmom(l2+1,m2+l2+1) + val*val -!**! enddo ! m2 -!**! end if -!**! enddo ! nacz2 -!**! enddo ! l2 -!**! fac = fac + fac1 * fac2 -!**! fac_angmom(:,:) = fac_angmom(:,:) + fac1 * fac2_angmom(:,:) -!**! else ! NOT m resolved -!**! fac2_angmom(:,:) = zero -!**! iatomf2 = 0 -!**! do l2 = 0, pao(neigh_species)%greatest_angmom -!**! do nacz2 = 1, pao(neigh_species)%angmom(l2)%n_zeta_in_angmom -!**! if((pao(neigh_species)%angmom(l2)%semicore(nacz2)==0) .OR. & -!**! (flag_pDOS_include_semicore)) then -!**! do m2 = -l2,l2 -!**! iatomf2 = iatomf2 + 1 -!**! wheremat = matrix_pos(matSFcoeff(spin_SF),iprim,j_in_halo,isf1,iatomf2) -!**! val = mat_p(matSFcoeff(spin_SF))%matrix(wheremat) -!**! fac2 = fac2 + val*val -!**! fac2_angmom(l2+1,1) = fac2_angmom(l2+1,1) + val*val -!**! enddo ! m2 -!**! end if -!**! enddo ! nacz2 -!**! enddo ! l2 -!**! fac = fac + fac1 * fac2 -!**! fac_angmom(:,1) = fac_angmom(:,1) + fac1 * fac2_angmom(:,1) -!**! if (iatomf2.ne.natomf2) call cq_abort("Error in NATOMF in the PDOS calculation.") -!**! endif ! flag_pDOS_angmom -!**! enddo ! isf1 -!**! if (flag_normalise_pDOS) then -!**! fac = fac / weight_pDOS(iwf) -!**! if (flag_pDOS_angmom) fac_angmom(:,:) = fac_angmom(:,:) / weight_pDOS(iwf) -!**! endif -!**! ! project on the neighbour atom -!**! do i=n_min,n_max -!**! projDOS(i,neigh_global_num) = projDOS(i,neigh_global_num) + tmp(i)*fac -!**! end do -!**! if (flag_pDOS_angmom) then -!**! if(flag_pDOS_lm) then ! (l,m) resolved -!**! do l2 = 0, pao(neigh_species)%greatest_angmom -!**! do m2 = -l2,l2 -!**! do i=n_min,n_max -!**! projDOS_angmom(i,neigh_global_num,l2+1,m2+l2+1) = & -!**! projDOS_angmom(i,neigh_global_num,l2+1,m2+l2+1) + tmp(i)*fac_angmom(l2+1,l2+m2+1) -!**! end do ! i -!**! end do ! m2 -!**! end do ! l2 -!**! else ! l resolved only -!**! do l2 = 0, pao(neigh_species)%greatest_angmom -!**! do i=n_min,n_max -!**! projDOS_angmom(i,neigh_global_num,l2+1,1) = projDOS_angmom(i,neigh_global_num,l2+1,1) & -!**! + tmp(i)*fac_angmom(l2+1,1) -!**! end do ! i -!**! end do ! l2 -!**! end if -!**! endif ! flag_pDOS_angmom -!**! end do ! neigh -!**! acc = acc + nsf1 -!**! end do ! memb -!**! end if ! nm_nodgroup -!**! end do ! part -!**! endif ! atomf -!**! end if ! flag_write_projected_DOS -!**! end do ! iwf -!**! end subroutine accumulate_DOS -!**! !!*** - !!****f* DiagModule/write_wavefn_coeffs !! !! NAME diff --git a/src/control.f90 b/src/control.f90 index 2a3ae7b0e..a39f48f7e 100644 --- a/src/control.f90 +++ b/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 diff --git a/src/global_module.f90 b/src/global_module.f90 index a4eef0650..dc4b4eed4 100644 --- a/src/global_module.f90 +++ b/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 diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index 70c67fe20..305883f74 100644 --- a/src/initial_read_module.f90 +++ b/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 !!$ diff --git a/src/initialisation_module.f90 b/src/initialisation_module.f90 index e42198254..4d8934cc2 100644 --- a/src/initialisation_module.f90 +++ b/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 !!$ !!$ !!$ diff --git a/tools/PostProcessing/local_module.f90 b/tools/PostProcessing/local_module.f90 index 7588d4be2..7ec43a618 100644 --- a/tools/PostProcessing/local_module.f90 +++ b/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 diff --git a/tools/PostProcessing/process_module.f90 b/tools/PostProcessing/process_module.f90 index 80b413b3c..81e1249cc 100644 --- a/tools/PostProcessing/process_module.f90 +++ b/tools/PostProcessing/process_module.f90 @@ -1,6 +1,11 @@ module process + use datatypes + implicit none + + ! Maximum possible number of spin components for simplicity + real(double), dimension(4) :: range_offset contains @@ -290,7 +295,8 @@ subroutine process_dos use datatypes use numbers, ONLY: zero, RD_ERR, twopi, half, one, two, four, six - use local, ONLY: eigenvalues, n_bands_total, nkp, wtk, efermi, flag_total_iDOS, flag_procwf_range_Ef + use local, ONLY: eigenvalues, n_bands_total, nkp, wtk, efermi, & + flag_total_iDOS, flag_procwf_range_Ef, flag_expand_range use read, ONLY: read_eigenvalues, read_psi_coeffs use global_module, ONLY: nspin, n_DOS, E_DOS_min, E_DOS_max, sigma_DOS use units, ONLY: HaToeV @@ -299,7 +305,7 @@ subroutine process_dos ! Local variables integer :: i_band, i_kp, i_spin, n_DOS_wid, n_band, n_min, n_max, i - real(double) :: Ebin, dE_DOS, a, pf_DOS, spin_fac + real(double) :: Ebin, dE_DOS, a, pf_DOS, spin_fac, peak_width real(double), dimension(nspin) :: total_electrons, iDOS_low real(double), dimension(:,:), allocatable :: total_DOS, iDOS real(double), dimension(:,:), allocatable :: occ @@ -315,77 +321,68 @@ subroutine process_dos allocate(total_DOS(n_DOS,nspin),iDOS(n_DOS,nspin),occ(n_bands_total,nkp)) total_DOS = zero iDOS = zero - ! Set limits on DOS output + ! Define peak width + peak_width = six*sigma_DOS + ! Offset for the energy range used for DOS display + range_offset = zero + do i_spin = 1, nspin + if(flag_procwf_range_Ef) range_offset(i_spin) = efermi(i_spin) + end do + ! Set limits and broaden energy range if needed if(abs(E_DOS_min)=E_DOS_min .and. & - eigenvalues(i_band, i_kp, i_spin)<=E_DOS_max) then - n_band = floor((eigenvalues(i_band, i_kp, i_spin) - E_DOS_min)/dE_DOS) + 1 + if(eigenvalues(i_band, i_kp, i_spin)>=E_DOS_min + range_offset(i_spin) .and. & + eigenvalues(i_band, i_kp, i_spin)<=E_DOS_max + range_offset(i_spin)) then + n_band = floor((eigenvalues(i_band, i_kp, i_spin) - (E_DOS_min + range_offset(i_spin)))/dE_DOS) + 1 n_min = n_band - n_DOS_wid if(n_min<1) n_min = 1 n_max = n_band + n_DOS_wid if(n_max>n_DOS) n_max = n_DOS do i = n_min, n_max - Ebin = real(i-1,double)*dE_DOS + E_DOS_min + Ebin = real(i-1,double)*dE_DOS + E_DOS_min + range_offset(i_spin) a = (Ebin-eigenvalues(i_band, i_kp, i_spin))/sigma_DOS total_DOS(i,i_spin) = total_DOS(i,i_spin) + wtk(i_kp)*pf_DOS*exp(-half*a*a) total_electrons(i_spin) = total_electrons(i_spin) + occ(i_band,i_kp)*wtk(i_kp)*pf_DOS*exp(-half*a*a) iDOS(i,i_spin) = iDOS(i,i_spin) + wtk(i_kp)*pf_DOS*exp(-half*a*a) end do - else if(eigenvalues(i_band, i_kp, i_spin)E_DOS_min .and. & - eigenvalues(i_band, i_kp, i_spin)E_DOS_min + range_offset(i_spin) .and. & + eigenvalues(i_band, i_kp, i_spin)n_DOS) n_max = n_DOS do i = n_min, n_max - Ebin = real(i-1,double)*dE_DOS + E_DOS_min + Ebin = real(i-1,double)*dE_DOS + E_DOS_min + range_offset(i_spin) a = (Ebin-eigenvalues(i_band, i_kp, i_spin))/sigma_DOS do i_atom = 1, n_atoms_pDOS i_spec = species_glob(pDOS_atom_index(i_atom)) @@ -603,10 +591,6 @@ subroutine process_pdos end if end do end do - if(flag_procwf_range_Ef) then - E_DOS_min = E_DOS_min - efermi(i_spin) - E_DOS_max = E_DOS_max - efermi(i_spin) - end if end do ! do i_spin = 1, n_spin ! Include spin factor and convert Ha to eV if(flag_l_resolved .and. flag_lm_resolved) then @@ -621,9 +605,11 @@ subroutine process_pdos pDOS = pDOS*spin_fac/HaToeV total_electrons = total_electrons*dE_DOS*spin_fac end if + write(*,fmt='(2x,"Integration of DOS in terms of absolute energies")') if(nspin==1) then - write(*,fmt='(2x,"Results of integrating pDOS between ",f11.3," and ",f11.3," Ha (electrons per atom).")') & - E_DOS_min, E_DOS_max + !write(*,fmt='(2x,"Results of integrating pDOS between ",f11.3," and ",f11.3," Ha (electrons per atom).")') & + write(*,fmt='(2x,"Results of integrating pDOS between ",f11.3," and Ef (electrons per atom).")') & + E_DOS_min + range_offset(1) if(flag_l_resolved) then ! l and l-m write(*,fmt='(4x," Atom Total l=0 l=1 l=2")') write(fmt_DOS,*) max_l + 2 ! Number of columns @@ -653,8 +639,8 @@ subroutine process_pdos end if else if(flag_l_resolved) then - write(*,fmt='(2x,"Results of integrating pDOS between ",f11.3," and ",f11.3," Ha (electrons per atom).")') & - E_DOS_min, E_DOS_max + write(*,fmt='(2x,"Results of integrating pDOS between ",f11.3," and Ef (electrons per atom).")') & + E_DOS_min + range_offset(1) write(*,fmt='(4x," Atom Spin Up l=0 l=1 l=2 Spin Down l=0 l=1 l=2")') write(fmt_DOS,*) 2*(max_l + 2) ! Number of columns fmt_DOS = '(4x,i7,x,'//trim(adjustl(fmt_DOS))//'f11.3)' @@ -677,8 +663,8 @@ subroutine process_pdos fmt_DOS = '('//trim(adjustl(fmt_DOS))//'f12.5)' end if else - write(*,fmt='(2x,"Results of integrating pDOS between ",f11.3," and ",f11.3," Ha (electrons per atom).")') & - E_DOS_min, E_DOS_max + write(*,fmt='(2x,"Results of integrating pDOS between ",f11.3," and Ef (electrons per atom).")') & + E_DOS_min + range_offset(1) write(*,fmt='(4x," Atom Spin Up Spin Down")') do i_atom = 1, n_atoms_pDOS write(*,fmt='(4x,i7,x,2f11.3)') pDOS_atom_index(i_atom), total_electrons(i_atom,1), total_electrons(i_atom,2) @@ -702,50 +688,26 @@ subroutine process_pdos write(17,fmt='("# Spin ",I1)') i_spin write(17,fmt='("# Original Fermi-level: ",f12.5," eV")') HaToeV*efermi(i_spin) write(17,fmt='("# DOS shifted relative to Fermi-level")') - if(flag_procwf_range_Ef) then - if(flag_l_resolved .and. flag_lm_resolved) then - write(17,fmt='("# Energy(eV) pDOS(/eV)")') - write(17,fmt='("# Total l=0 l=1 l=2")') - do i=1, n_DOS - write(17,fmt=fmt_dos) HaToeV*(E_DOS_min + dE_DOS*real(i-1,double)), pDOS(i_atom,i,i_spin), & - ((pDOS_lm(i_m,i_l,i_atom,i,i_spin),i_m=-i_l,i_l),i_l=0,pao(i_spec)%greatest_angmom) - end do - else if(flag_l_resolved) then - write(17,fmt='("# Energy(eV) pDOS(/eV)")') - write(17,fmt='("# Total l=0 l=1 l=2")') - do i=1, n_DOS - write(17,fmt=fmt_dos) HaToeV*(E_DOS_min + dE_DOS*real(i-1,double)), pDOS(i_atom,i,i_spin), & - pDOS_l(0:pao(i_spec)%greatest_angmom,i_atom,i,i_spin) - end do - else - write(17,fmt='("# Energy(eV) pDOS(/eV)")') - do i=1, n_DOS - write(17,fmt='(2f12.5)') HaToeV*(E_DOS_min + dE_DOS*real(i-1,double)), & - pDOS(i_atom,i,i_spin) - end do - end if + if(flag_l_resolved .and. flag_lm_resolved) then + write(17,fmt='("# Energy(eV) pDOS(/eV)")') + write(17,fmt='("# Total l=0 l=1 l=2")') + do i=1, n_DOS + write(17,fmt=fmt_dos) HaToeV*(E_DOS_min + range_offset(i_spin) - efermi(i_spin) + dE_DOS*real(i-1,double)), pDOS(i_atom,i,i_spin), & + ((pDOS_lm(i_m,i_l,i_atom,i,i_spin),i_m=-i_l,i_l),i_l=0,pao(i_spec)%greatest_angmom) + end do + else if(flag_l_resolved) then + write(17,fmt='("# Energy(eV) pDOS(/eV)")') + write(17,fmt='("# Total l=0 l=1 l=2")') + do i=1, n_DOS + write(17,fmt=fmt_dos) HaToeV*(E_DOS_min + range_offset(i_spin) - efermi(i_spin) + dE_DOS*real(i-1,double)), pDOS(i_atom,i,i_spin), & + pDOS_l(0:pao(i_spec)%greatest_angmom,i_atom,i,i_spin) + end do else - if(flag_l_resolved .and. flag_lm_resolved) then - write(17,fmt='("# Energy(eV) pDOS(/eV)")') - write(17,fmt='("# Total l=0 l=1 l=2")') - do i=1, n_DOS - write(17,fmt=fmt_dos) HaToeV*(E_DOS_min + dE_DOS*real(i-1,double)-efermi(i_spin)), pDOS(i_atom,i,i_spin), & - ((pDOS_lm(i_m,i_l,i_atom,i,i_spin),i_m=-i_l,i_l),i_l=0,pao(i_spec)%greatest_angmom) - end do - else if(flag_l_resolved) then - write(17,fmt='("# Energy(eV) pDOS(/eV)")') - write(17,fmt='("# Total l=0 l=1 l=2")') - do i=1, n_DOS - write(17,fmt=fmt_dos) HaToeV*(E_DOS_min + dE_DOS*real(i-1,double)-efermi(i_spin)), pDOS(i_atom,i,i_spin), & - pDOS_l(0:pao(i_spec)%greatest_angmom,i_atom,i,i_spin) - end do - else - write(17,fmt='("# Energy(eV) pDOS(/eV)")') - do i=1, n_DOS - write(17,fmt='(2f12.5)') HaToeV*(E_DOS_min + dE_DOS*real(i-1,double)-efermi(i_spin)), & - pDOS(i_atom,i,i_spin) - end do - end if + write(17,fmt='("# Energy(eV) pDOS(/eV)")') + do i=1, n_DOS + write(17,fmt='(2f12.5)') HaToeV*(E_DOS_min + range_offset(i_spin) - efermi(i_spin) + dE_DOS*real(i-1,double)), & + pDOS(i_atom,i,i_spin) + end do end if write(17,fmt='("&")') end do diff --git a/tools/PostProcessing/read_module.f90 b/tools/PostProcessing/read_module.f90 index 7fa2631a1..274b04433 100644 --- a/tools/PostProcessing/read_module.f90 +++ b/tools/PostProcessing/read_module.f90 @@ -231,7 +231,7 @@ subroutine read_input E_procwf_min = fdf_double('Process.min_wf_E',E_wf_min) E_procwf_max = fdf_double('Process.max_wf_E',E_wf_max) ! Is the range relative to Ef (T) or absolute (F) - flag_procwf_range_Ef = fdf_boolean('Process.WFRangeRelative',.true.) + flag_procwf_range_Ef = fdf_boolean('Process.WFRangeRelative',flag_wf_range_Ef) n_bands_process = fdf_integer('Process.noWF',0) if(n_bands_process>0) then allocate(band_proc_no(n_bands_process)) @@ -269,12 +269,29 @@ subroutine read_input flag_outputWF_real = .false. if (leqi(job,'ban')) flag_outputWF_real = fdf_boolean('Process.outputWF_real',.false.) ! DOS - ! Add flag for window relative to Fermi level - E_DOS_min = fdf_double('Process.min_DOS_E',E_wf_min) - E_DOS_max = fdf_double('Process.max_DOS_E',E_wf_max) - sigma_DOS = fdf_double('Process.sigma_DOS',0.001_double) ! Better than adaptive - n_DOS = fdf_integer('Process.n_DOS',1001) - flag_total_iDOS = fdf_boolean('Process.TotalIntegratedDOS',.false.) + flag_wf_range_Ef = fdf_boolean('IO.WFRangeRelative',.true.) + flag_procwf_range_Ef = fdf_boolean('Process.WFRangeRelative',flag_wf_range_Ef) + if(i_job==6.or.i_job==7) then + ! Add flag for window relative to Fermi level + E_DOS_min = fdf_double('Process.min_DOS_E',E_wf_min) + E_DOS_max = fdf_double('Process.max_DOS_E',E_wf_max) + ! ExpandRange allows for the broadening of peaks in energy range + ! If user did not set limits, default to expanding range + if(((abs(E_DOS_min-E_wf_min)>1e-8_double).or.(abs(E_DOS_max-E_wf_max)>1e-8_double)).or. & + (abs(E_wf_min)>1e-8_double.or.abs(E_wf_max)>1e-8_double)) then ! At least one set + flag_expand_range = fdf_boolean('Process.ExpandRange',.false.) + flag_wf_range_Ef = fdf_boolean('IO.WFRangeRelative',.true.) + flag_procwf_range_Ef = fdf_boolean('Process.WFRangeRelative',flag_wf_range_Ef) + else ! Defaults used + flag_expand_range = fdf_boolean('Process.ExpandRange',.true.) + ! If the limits are set automatically then don't treat limits as relative + flag_wf_range_Ef = fdf_boolean('IO.WFRangeRelative',.false.) + flag_procwf_range_Ef = fdf_boolean('Process.WFRangeRelative',flag_wf_range_Ef) + end if + sigma_DOS = fdf_double('Process.sigma_DOS',0.001_double) ! Better than adaptive + n_DOS = fdf_integer('Process.n_DOS',1001) + flag_total_iDOS = fdf_boolean('Process.TotalIntegratedDOS',.false.) + end if if(i_job==7) then ! If no limits specified, cover whole range if(abs(E_wf_max-E_wf_min)<1e-8_double) then @@ -282,8 +299,6 @@ subroutine read_input E_wf_max = BIG end if flag_wf_range = .true. - flag_wf_range_Ef = fdf_boolean('IO.WFRangeRelative',.true.) - flag_procwf_range_Ef = fdf_boolean('Process.WFRangeRelative',.true.) flag_l_resolved = fdf_boolean('Process.pDOS_l_resolved',.false.) flag_lm_resolved = fdf_boolean('Process.pDOS_lm_resolved',.false.) if(flag_lm_resolved .and. (.not.flag_l_resolved)) flag_l_resolved = .true.