From 2239f94edd784cc182f425956c1d750f4ba16dd3 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Tue, 14 Jul 2026 13:42:26 +0100 Subject: [PATCH 1/9] Updates to occupancy calculations Added flag to allow user to specify integer occupancies Track valence band maximum and conduction band minimum in occupy (except for Methfessel-Paxton) and write out gap when present, specifying direct or indirect and location --- src/DiagModule.f90 | 115 ++++++++++++++++++++++++++++++++---- src/initial_read_module.f90 | 3 +- src/pao2blip.f90 | 2 +- 3 files changed, 108 insertions(+), 12 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index a8af360fb..4993dfd71 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -260,9 +260,11 @@ module DiagModule type(DistributeData) :: DistribH, DistribS ! Fermi Energy - real(double), dimension(2) :: Efermi + real(double), dimension(2) :: Efermi, vbm, cbm, gap_d + integer, dimension(2) :: vbm_k, cbm_k, gap_k integer, dimension(2) :: band_ef - + logical :: flag_integer_occ + ! K-point data - here so that reading of k-points can take place in ! different routine to FindEvals integer :: nkp @@ -316,6 +318,7 @@ module DiagModule ! that Methfessel-Paxton approximation may casue the bracket search ! algorithm to fail.) integer :: max_brkt_iterations + real(double), parameter :: tolElec = 1.0e-8_double contains @@ -602,7 +605,38 @@ subroutine FindEvals(electrons) end if ! Find Fermi level, given the eigenvalues at all k-points (in w) ! if (me < proc_rows*proc_cols) then + vbm = -BIG + cbm = BIG + vbm_k = 0 + cbm_k = 0 + gap_k = 0 + gap_d = BIG call findFermi(electrons, evals, matrix_size, nkp, Efermi, occ) + ! Test for gap + if(flag_smear_type==0.or.flag_integer_occ) then + if((cbm(1)-vbm(1))>two*kT.and.(cbm(nspin)-vbm(nspin)>two*kT)) then + if(inode==ionode) then + do spin=1,nspin + write(io_lun,fmt='(4x,"Spin ",i1," gap found. VBM=",f12.5,"Ha CBM=",f12.5,"Ha Gap ",f12.5,"Ha")') & + spin,vbm(spin),cbm(spin),cbm(spin)-vbm(spin) + if(vbm_k(spin)==cbm_k(spin)) then + write(io_lun,fmt='(4x,"Direct gap found at ",3f8.4)') kk(:,gap_k(spin)) + else + write(io_lun,fmt='(4x,"Indirect gap found; smallest direct gap is ",f12.5," Ha at ",3f6.2)') & + gap_d(spin),kk(:,gap_k(spin)) + write(io_lun,fmt='(4x,"VBM is at ",3f8.4)') kk(:,vbm_k(spin)) + write(io_lun,fmt='(4x,"CBM is at ",3f8.4)') kk(:,cbm_k(spin)) + end if + end do + end if + ! Adjust Fermi level to lie at mid-gap and revisit occupancies + do spin=1,nspin + Efermi(spin) = half*(cbm(spin)+vbm(spin)) + locc(spin) = electrons(spin) + end do + call occupy(occ, evals, Efermi, locc, matrix_size, nkp) + end if + end if ! Allocate space to expand eigenvectors into (i.e. when reversing ! ScaLAPACK distribution) allocate(expH(matrix_size,prim_size,nspin), STAT=stat) @@ -2362,10 +2396,36 @@ subroutine findFermi(electrons, eig, nbands, nkp, Ef, occ) ! local variables real(double) :: electrons_total - real(double), parameter :: tolElec = 1.0e-6_double real(double) :: locals_occ, localt_occ - integer :: ikp, i, ispin - + integer :: ikp, i, spin, ne + + ! For integer occupancies + if(flag_integer_occ) then + occ = zero + gap_d = BIG + do spin=1,nspin + ne = int(electrons(spin)) + electrons_total = zero + vbm(spin) = maxval(eig(ne,:,spin)) + cbm(spin) = minval(eig(ne+1,:,spin)) + vbm_k(spin) = maxloc(eig(ne,:,spin),dim=1) + cbm_k(spin) = minloc(eig(ne+1,:,spin),dim=1) + do ikp=1,nkp + occ(1:ne,ikp,spin) = wtk(ikp) + if(eig(ne+1,ikp,spin) - eig(ne,ikp,spin)half) then + if(vbm(ss)ebands(iband,ikp,ss)) then + cbm_k(ss) = ikp + cbm(ss) = ebands(iband,ikp,ss) + if(ebands(iband,ikp,ss) - ebands(iband-1,ikp,ss)half) then + ! if(vbm(ss)ebands(iband,ikp,ss)) then + ! cbm_k(ss) = ikp + ! cbm(ss) = ebands(iband,ikp,ss) + ! gap_d(ss) = min(gap_d(ss),ebands(iband,ikp,ss)-ebands(iband-1,ikp,ss)) + ! end if + !end if case default call cq_abort ("FindEvals: Smearing flag not recognised",& flag_smear_type) diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index 08b6065ec..a804887a7 100644 --- a/src/initial_read_module.f90 +++ b/src/initial_read_module.f90 @@ -3228,7 +3228,7 @@ subroutine readDiagInfo use DiagModule, only: nkp, kk, wtk, kT, maxefermi, & flag_smear_type, iMethfessel_Paxton, & max_brkt_iterations, gaussian_height, & - finess, NElec_less + finess, NElec_less, flag_integer_occ use energy, only: SmearingType, MPOrder use memory_module, only: reg_alloc_mem, reg_dealloc_mem, & type_dbl @@ -3263,6 +3263,7 @@ subroutine readDiagInfo ! Read Control Flags associated to diagonalisation method maxefermi = fdf_integer('Diag.MaxEfIter',50) + flag_integer_occ = fdf_boolean('Diag.IntegerOccs',.false.) kT = fdf_double('Diag.kT',0.001_double) ! Method to approximate step function for occupation number flag_smear_type = fdf_integer('Diag.SmearingType',0) diff --git a/src/pao2blip.f90 b/src/pao2blip.f90 index 57ae06513..00ce0e168 100644 --- a/src/pao2blip.f90 +++ b/src/pao2blip.f90 @@ -806,7 +806,7 @@ subroutine blips_in_star(inode,ionode,sym_type,& scal_prod_sym(n) = scal_prod(m)/coeff(n) end do if((inode == ionode).and.(iprint_basis >= 4)) then - write(unit=io_lun,fmt='(/6x" symmetric scalar product:"/)') + write(unit=io_lun,fmt='(/6x," symmetric scalar product:"/)') do n = 1, n_star_in_sphere write(unit=io_lun,fmt='(6x,3x,i5,3x,e15.6)') n, scal_prod_sym(n) end do From 34687ac806ea83adf4f111bdc1eea9a4d4c9e582 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Fri, 31 Jul 2026 14:43:40 +0100 Subject: [PATCH 2/9] Tidied output and added documentation --- docs/groundstate.rst | 14 +++-- docs/input_tags.rst | 5 ++ src/DiagModule.f90 | 105 +++++++++++++++++++++++++----------- src/initial_read_module.f90 | 2 +- src/minimise.f90 | 5 +- 5 files changed, 96 insertions(+), 35 deletions(-) diff --git a/docs/groundstate.rst b/docs/groundstate.rst index 354f8db46..5fe1e8c23 100644 --- a/docs/groundstate.rst +++ b/docs/groundstate.rst @@ -144,10 +144,18 @@ Go to :ref:`top `. Electronic occupation smearing ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The occupation numbers of the eigenstates are slightly smeared near +By default, the occupation numbers of the eigenstates are slightly smeared near the Fermi level, following common practice. The default smearing type -is Fermi-Dirac smearing with a temperature (in Hartrees) set with the -flag ``Diag.kT`` which defaults to 0.001Ha. +is Fermi-Dirac smearing (``Diag.SmearingType 0``) with a temperature (in Hartrees) set with the +flag ``Diag.kT`` which defaults to 0.001Ha. If you know that you have an insulator with a +finite gap, you can force integer occupancies with the flag: + + :: + + Diag.IntegerOccs T + +This is false by default, and will lead to significant problems if the system is metallic +at any point during the self-consistent ground state search, so should be considered carefully. The Methfessel-Paxton approach :cite:`g-Methfessel:1989ny` to occupations allows much higher smearing temperatures with minimal effect on the free energy (and diff --git a/docs/input_tags.rst b/docs/input_tags.rst index 62479fd57..41946f908 100644 --- a/docs/input_tags.rst +++ b/docs/input_tags.rst @@ -733,6 +733,11 @@ Diag.kT (*real*) *default*: 0.001 +Diag.IntegerOccs (*logical*) + Specifies integer occupation of eigenstates (overrides ``Diag.SmearingType``) + + *default*: False + Diag.MPOrder (*integer*) Order of Bessel function approximation to delta-function used in Methfessel-Paxton smearing diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index 4993dfd71..b9227ad48 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -263,7 +263,7 @@ module DiagModule real(double), dimension(2) :: Efermi, vbm, cbm, gap_d integer, dimension(2) :: vbm_k, cbm_k, gap_k integer, dimension(2) :: band_ef - logical :: flag_integer_occ + logical :: flag_integer_occ, flag_gap ! K-point data - here so that reading of k-points can take place in ! different routine to FindEvals @@ -611,31 +611,17 @@ subroutine FindEvals(electrons) cbm_k = 0 gap_k = 0 gap_d = BIG + flag_gap = .false. call findFermi(electrons, evals, matrix_size, nkp, Efermi, occ) ! Test for gap - if(flag_smear_type==0.or.flag_integer_occ) then - if((cbm(1)-vbm(1))>two*kT.and.(cbm(nspin)-vbm(nspin)>two*kT)) then - if(inode==ionode) then - do spin=1,nspin - write(io_lun,fmt='(4x,"Spin ",i1," gap found. VBM=",f12.5,"Ha CBM=",f12.5,"Ha Gap ",f12.5,"Ha")') & - spin,vbm(spin),cbm(spin),cbm(spin)-vbm(spin) - if(vbm_k(spin)==cbm_k(spin)) then - write(io_lun,fmt='(4x,"Direct gap found at ",3f8.4)') kk(:,gap_k(spin)) - else - write(io_lun,fmt='(4x,"Indirect gap found; smallest direct gap is ",f12.5," Ha at ",3f6.2)') & - gap_d(spin),kk(:,gap_k(spin)) - write(io_lun,fmt='(4x,"VBM is at ",3f8.4)') kk(:,vbm_k(spin)) - write(io_lun,fmt='(4x,"CBM is at ",3f8.4)') kk(:,cbm_k(spin)) - end if - end do - end if - ! Adjust Fermi level to lie at mid-gap and revisit occupancies - do spin=1,nspin - Efermi(spin) = half*(cbm(spin)+vbm(spin)) - locc(spin) = electrons(spin) - end do - call occupy(occ, evals, Efermi, locc, matrix_size, nkp) - end if + if((flag_smear_type==0.or.flag_integer_occ).and.flag_gap) then + if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps + ! Adjust Fermi level to lie at mid-gap and revisit occupancies + do spin=1,nspin + Efermi(spin) = half*(cbm(spin)+vbm(spin)) + locc(spin) = electrons(spin) + end do + call occupy(occ, evals, Efermi, locc, matrix_size, nkp) end if ! Allocate space to expand eigenvectors into (i.e. when reversing ! ScaLAPACK distribution) @@ -2424,6 +2410,7 @@ subroutine findFermi(electrons, eig, nbands, nkp, Ef, occ) end if Ef(spin) = half*(vbm(spin) + cbm(spin)) end do + if((cbm(1)-vbm(1))>two*kT.and.(cbm(nspin)-vbm(nspin)>two*kT)) flag_gap = .true. return end if ! Integer occupations if (nspin == 2) then @@ -2452,6 +2439,7 @@ subroutine findFermi(electrons, eig, nbands, nkp, Ef, occ) end do end do end if + if((cbm(1)-vbm(1))>two*kT.and.(cbm(nspin)-vbm(nspin)>two*kT)) flag_gap = .true. return end subroutine findFermi !!***** @@ -3071,16 +3059,12 @@ subroutine occupy(occu, ebands, Ef, electrons, nbands, nkp, spin) gap_d(ss) = ebands(iband,ikp,ss) - ebands(iband-1,ikp,ss) gap_k(ss) = ikp end if - !gap_d(ss) = min(gap_d(ss),ebands(iband,ikp,ss)-ebands(iband-1,ikp,ss)) end if end if case (1) ! Methfessel Paxton smearing locc = MP_step(ebands(iband,ikp,ss) - Ef(ss), iMethfessel_Paxton, kT) - occu(iband,ikp,ss) = & - wtk(ikp) * locc - !wtk(ikp) * MP_step(ebands(iband,ikp,ss) - Ef(ss), & - !iMethfessel_Paxton, kT) - ! It really doesn't make sense to look for vbm and cbm with MP but this is how + occu(iband,ikp,ss) = wtk(ikp) * locc + ! It really doesn't make sense to look for vbm and cbm with MP but this is how we could !if(locc>half) then ! if(vbm(ss)1) then + write(io_lun,fmt='(4x,"Spin ",i1," gap found. VBM=",f12.5," ",a2," CBM=", & + f12.5," ",a2," Gap ",f12.5," ",a2)') & + spin,en_conv*vbm(spin),en_units(energy_units), & + en_conv*cbm(spin),en_units(energy_units),& + en_conv*(cbm(spin)-vbm(spin)),en_units(energy_units) + else + write(io_lun,fmt='(4x,"Gap found. VBM=",f12.5," ",a2," CBM=",f12.5," ",a2," Gap ",f12.5," ",a2)') & + en_conv*vbm(spin),en_units(energy_units),& + en_conv*cbm(spin),en_units(energy_units),& + en_conv*(cbm(spin)-vbm(spin)),en_units(energy_units) + end if + if(vbm_k(spin)==cbm_k(spin)) then + write(io_lun,fmt='(4x,"Direct gap found at ",3f8.4," (1/a0)")') kk(:,gap_k(spin)) + else + write(io_lun,fmt='(4x,"Indirect gap found; smallest direct gap is ",f12.5," ",a2," at ",3f6.2," (1/a0)")') & + en_conv*gap_d(spin),en_units(energy_units),kk(:,gap_k(spin)) + write(io_lun,fmt='(4x,"VBM is at ",3f8.4," (1/a0)")') kk(:,vbm_k(spin)) + write(io_lun,fmt='(4x,"CBM is at ",3f8.4," (1/a0)")') kk(:,cbm_k(spin)) + end if + end do + end if + return + end subroutine write_gaps + !!*** + ! ----------------------------------------------------------------------------- ! Function fermi ! ----------------------------------------------------------------------------- diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index a804887a7..3b0b80bb5 100644 --- a/src/initial_read_module.f90 +++ b/src/initial_read_module.f90 @@ -3761,7 +3761,7 @@ subroutine readDiagInfo 3 format(8x,'Proc grid (row, col): ',2i5) 4 format(/8x,'***WARNING***',/,2x,& 'No Kpoints block found - defaulting to Gamma point') -51 format(/8x,i4,' symmetry inequivalent Kpoints in Cartesian form (1/A): ') +51 format(/8x,i4,' symmetry inequivalent Kpoints in Cartesian form (1/a0): ') 7 format(8x,' All ',i4,' Kpoints in fractional coordinates: ') 9 format(/8x,'***WARNING***',/,8x,& 'Specified Kpoint mesh shift in fractional coords >= 1.0.') diff --git a/src/minimise.f90 b/src/minimise.f90 index d119869c6..ffdb9f5ac 100644 --- a/src/minimise.f90 +++ b/src/minimise.f90 @@ -167,7 +167,7 @@ subroutine get_E_and_F(fixed_potential, vary_mu, total_energy, & use multisiteSF_module,only: flag_LFD_nonSCF, flag_mix_LFD_SCF use units use io_module, only: return_prefix - use DiagModule, only: nkp + use DiagModule, only: nkp, write_gaps use H_matrix_module, only: flag_write_locps, locps_output, get_H_matrix, get_occ_matrix use maxima_module, only: maxngrid use polarisation, only: get_polarisation @@ -372,6 +372,9 @@ subroutine get_E_and_F(fixed_potential, vary_mu, total_energy, & ! Get occupation matrix for output density matrix if(flag_DFTplusU) call get_occ_matrix + ! Output gap if appropriate + if(inode==ionode.and.iprint_minE + min_layer >= 0) call write_gaps + !****lat<$ call final_energy(nkp,backtrace_level) !****lat>$ From 78559c232c2b09277af2eb883365899d54308202 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Fri, 31 Jul 2026 16:01:37 +0100 Subject: [PATCH 3/9] Update to account for degenerate levels Test 007 identified a problem when there were degenerate levels which has now been fixed --- src/DiagModule.f90 | 13 ++++++++----- testsuite/test_check_output.py | 5 ++++- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index b9227ad48..ec6eda614 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -3043,15 +3043,18 @@ subroutine occupy(occu, ebands, Ef, electrons, nbands, nkp, spin) select case (flag_smear_type) case (0) ! Fermi smearing locc = fermi(ebands(iband,ikp,ss) - Ef(ss), kT) - occu(iband,ikp,ss) = & - wtk(ikp) * locc - !wtk(ikp) * fermi(ebands(iband,ikp,ss) - Ef(ss), kT) - if(locc>half) then + occu(iband,ikp,ss) = wtk(ikp) * locc + ! Gap detection + ! Check for VBM and allow for degenerate states + if(locc-half>-very_small) then!locc>=half) then if(vbm(ss)-very_small) then!locc<=half) then if(cbm(ss)>ebands(iband,ikp,ss)) then cbm_k(ss) = ikp cbm(ss) = ebands(iband,ikp,ss) diff --git a/testsuite/test_check_output.py b/testsuite/test_check_output.py index 3279b25d8..796acc496 100644 --- a/testsuite/test_check_output.py +++ b/testsuite/test_check_output.py @@ -50,7 +50,10 @@ def precision(key='_'): return 1e-4 ''' - return 1e-4 + if(key == 'Total stress'): + return 1e-3 + else: + return 1e-4 @pytest.fixture def testsuite_directory(): From bc89e3eb906c6d86b25e5bc5f2e15746ad3ae76c Mon Sep 17 00:00:00 2001 From: David Bowler Date: Wed, 5 Aug 2026 17:21:37 +0100 Subject: [PATCH 4/9] Updated occupancies and post-processing Occupancies now correct for all spin cases, and made zero of energy for fixed spin common between two channels --- src/DiagModule.f90 | 117 +++++++++++++++++------- src/initial_read_module.f90 | 12 ++- tools/PostProcessing/process_module.f90 | 11 ++- tools/PostProcessing/read_module.f90 | 21 ++++- 4 files changed, 118 insertions(+), 43 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index ec6eda614..37c638b16 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -263,7 +263,8 @@ module DiagModule real(double), dimension(2) :: Efermi, vbm, cbm, gap_d integer, dimension(2) :: vbm_k, cbm_k, gap_k integer, dimension(2) :: band_ef - logical :: flag_integer_occ, flag_gap + logical :: flag_integer_occ + logical, dimension(2) :: flag_gap ! K-point data - here so that reading of k-points can take place in ! different routine to FindEvals @@ -614,14 +615,35 @@ subroutine FindEvals(electrons) flag_gap = .false. call findFermi(electrons, evals, matrix_size, nkp, Efermi, occ) ! Test for gap - if((flag_smear_type==0.or.flag_integer_occ).and.flag_gap) then + if(flag_smear_type==0) then + if(flag_fix_spin_population) then ! Separate fermi levels + do spin=1,nspin + if(flag_gap(spin)) then + Efermi(spin) = half*(cbm(spin)+vbm(spin)) + locc(spin) = electrons(spin) + if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps(spin_ch=spin) + call occupy(occ, evals, Efermi, locc, matrix_size, nkp, spin=spin) + if(inode==ionode.and.iprint_DM + min_layer >= 2) & + write(io_lun,'(4x, "Adjusted Fermi level for spin ", i2, " is ", f12.5)') & + spin,Efermi(spin) + end if + end do + else + if(flag_gap(1) .or. flag_gap(2)) then ! Fermi levels should be same + if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps + ! Adjust Fermi level to lie at mid-gap and revisit occupancies + do spin=1,nspin + Efermi(spin) = half*(minval(cbm)+maxval(vbm)) + locc(spin) = electrons(spin) + if(inode==ionode.and.iprint_DM + min_layer >= 2) & + write(io_lun,'(4x, "Adjusted Fermi level for spin ", i2, " is ", f12.5)') & + spin,Efermi(spin) + end do + call occupy(occ, evals, Efermi, locc, matrix_size, nkp) + end if + end if + else if(flag_integer_occ) then ! There has to be a gap! if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps - ! Adjust Fermi level to lie at mid-gap and revisit occupancies - do spin=1,nspin - Efermi(spin) = half*(cbm(spin)+vbm(spin)) - locc(spin) = electrons(spin) - end do - call occupy(occ, evals, Efermi, locc, matrix_size, nkp) end if ! Allocate space to expand eigenvectors into (i.e. when reversing ! ScaLAPACK distribution) @@ -948,6 +970,8 @@ subroutine FindEvals(electrons) do j = 1, matrix_size ! Calculate entropic contribution to electronic energy select case (flag_smear_type) + case (-1) ! Integer occupancies + entropy = zero case (0) ! Fermi smearing if (occ(j,kp,spin) > RD_ERR .and. & (wtk(kp) - occ(j,kp,spin)) > RD_ERR) then @@ -2420,8 +2444,12 @@ subroutine findFermi(electrons, eig, nbands, nkp, Ef, occ) end if if (flag_fix_spin_population .or. nspin == 1) then call findFermi_fixspin(electrons, eig, nbands, nkp, Ef, occ) + do spin=1,nspin + if((cbm(spin)-vbm(spin))>two*kT) flag_gap(spin) = .true. + end do else call findFermi_varspin(electrons_total, eig, nbands, nkp, Ef, occ) + if((cbm(1)-vbm(1))>two*kT.and.(cbm(nspin)-vbm(nspin)>two*kT)) flag_gap = .true. end if if(flag_DeltaSCF.AND.flag_excite) then if(nspin==1) then @@ -2439,7 +2467,6 @@ subroutine findFermi(electrons, eig, nbands, nkp, Ef, occ) end do end do end if - if((cbm(1)-vbm(1))>two*kT.and.(cbm(nspin)-vbm(nspin)>two*kT)) flag_gap = .true. return end subroutine findFermi !!***** @@ -3021,7 +3048,7 @@ subroutine occupy(occu, ebands, Ef, electrons, nbands, nkp, spin) integer, optional, intent(in) :: spin ! local variables - integer :: ikp, iband, ss + integer :: ikp, iband, ss, ne integer :: ss_start, ss_end real(double) :: locc @@ -3032,7 +3059,16 @@ subroutine occupy(occu, ebands, Ef, electrons, nbands, nkp, spin) ss_start = 1 ss_end = nspin end if - + if(flag_integer_occ) then + do ss=ss_start,ss_end + occ(:,:,ss) = zero + ne = int(electrons(ss)) + do ikp=1,nkp + occ(1:ne,ikp,ss) = wtk(ikp) + end do + end do + return + end if electrons = zero labspin: do ss = ss_start, ss_end vbm(ss) = -BIG @@ -3141,36 +3177,49 @@ end subroutine occupy !! MODIFICATION HISTORY !! !! SOURCE - subroutine write_gaps + subroutine write_gaps(spin_ch) use units use global_module, only: nspin implicit none - integer :: spin + ! Passed variables + integer, optional :: spin_ch - if((flag_smear_type==0.or.flag_integer_occ).and.flag_gap.and.myid==0) then - do spin=1,nspin - if(nspin>1) then - write(io_lun,fmt='(4x,"Spin ",i1," gap found. VBM=",f12.5," ",a2," CBM=", & - f12.5," ",a2," Gap ",f12.5," ",a2)') & - spin,en_conv*vbm(spin),en_units(energy_units), & - en_conv*cbm(spin),en_units(energy_units),& - en_conv*(cbm(spin)-vbm(spin)),en_units(energy_units) - else - write(io_lun,fmt='(4x,"Gap found. VBM=",f12.5," ",a2," CBM=",f12.5," ",a2," Gap ",f12.5," ",a2)') & - en_conv*vbm(spin),en_units(energy_units),& - en_conv*cbm(spin),en_units(energy_units),& - en_conv*(cbm(spin)-vbm(spin)),en_units(energy_units) - end if - if(vbm_k(spin)==cbm_k(spin)) then - write(io_lun,fmt='(4x,"Direct gap found at ",3f8.4," (1/a0)")') kk(:,gap_k(spin)) - else - write(io_lun,fmt='(4x,"Indirect gap found; smallest direct gap is ",f12.5," ",a2," at ",3f6.2," (1/a0)")') & - en_conv*gap_d(spin),en_units(energy_units),kk(:,gap_k(spin)) - write(io_lun,fmt='(4x,"VBM is at ",3f8.4," (1/a0)")') kk(:,vbm_k(spin)) - write(io_lun,fmt='(4x,"CBM is at ",3f8.4," (1/a0)")') kk(:,cbm_k(spin)) + ! Local variables + integer :: spin, spin_st, spin_end + + if(present(spin_ch)) then + spin_st = spin_ch + spin_end = spin_ch + else + spin_st = 1 + spin_end = nspin + end if + if((flag_smear_type==0.or.flag_integer_occ).and.myid==0) then + do spin=spin_st,spin_end + if(flag_gap(spin)) then + if(nspin>1) then!spin_end-spin_st>0) then + write(io_lun,fmt='(4x,"Spin ",i1," gap found. VBM=",f12.5," ",a2," CBM=", & + f12.5," ",a2," Gap ",f12.5," ",a2)') & + spin,en_conv*vbm(spin),en_units(energy_units), & + en_conv*cbm(spin),en_units(energy_units),& + en_conv*(cbm(spin)-vbm(spin)),en_units(energy_units) + else + write(io_lun,fmt='(4x,"Gap found. VBM=",f12.5," ",a2," CBM=",f12.5," ",a2," Gap ",f12.5," ",a2)') & + en_conv*vbm(spin),en_units(energy_units),& + en_conv*cbm(spin),en_units(energy_units),& + en_conv*(cbm(spin)-vbm(spin)),en_units(energy_units) + end if + if(vbm_k(spin)==cbm_k(spin)) then + write(io_lun,fmt='(4x,"Direct gap found at ",3f8.4," (1/a0)")') kk(:,gap_k(spin)) + else + write(io_lun,fmt='(4x,"Indirect gap found; smallest direct gap is ",f12.5," ",a2," at ",3f6.2," (1/a0)")') & + en_conv*gap_d(spin),en_units(energy_units),kk(:,gap_k(spin)) + write(io_lun,fmt='(4x,"VBM is at ",3f8.4," (1/a0)")') kk(:,vbm_k(spin)) + write(io_lun,fmt='(4x,"CBM is at ",3f8.4," (1/a0)")') kk(:,cbm_k(spin)) + end if end if end do end if diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index 3b0b80bb5..4b4a623d2 100644 --- a/src/initial_read_module.f90 +++ b/src/initial_read_module.f90 @@ -2987,6 +2987,8 @@ subroutine write_info(titles, mu, vary_mu, HNL_fac, NODES) write(io_lun,fmt='(6x,"Solving for the K matrix using ",a16)') 'diagonalisation ' if(iprint_init>0) then select case (flag_smear_type) + case(-1) + write(io_lun,'(6x,"Using integer occupancies")') case (0) write(io_lun,'(6x,"Using Fermi-Dirac smearing")') case (1) @@ -3219,7 +3221,8 @@ subroutine readDiagInfo use functions, only: is_prime use global_module, only: iprint_init, rcellx, rcelly, rcellz, & area_general, ni_in_cell, numprocs, & - species_glob, io_lun, io_ase, ase_file, write_ase, flag_calc_pol + species_glob, io_lun, io_ase, ase_file, write_ase, flag_calc_pol, & + flag_fix_spin_population, nspin use numbers, only: zero, one, two, pi, RD_ERR, half use GenComms, only: cq_abort, cq_warn, gcopy use input_module @@ -3263,10 +3266,15 @@ subroutine readDiagInfo ! Read Control Flags associated to diagonalisation method maxefermi = fdf_integer('Diag.MaxEfIter',50) - flag_integer_occ = fdf_boolean('Diag.IntegerOccs',.false.) kT = fdf_double('Diag.kT',0.001_double) ! Method to approximate step function for occupation number flag_smear_type = fdf_integer('Diag.SmearingType',0) + flag_integer_occ = fdf_boolean('Diag.IntegerOccs',.false.) + if(flag_integer_occ) then + if(nspin>1.and.(.not.flag_fix_spin_population)) & + call cq_abort("Cannot have free spin population with integer occupancies") + flag_smear_type = -1 + end if SmearingType = flag_smear_type iMethfessel_Paxton = fdf_integer('Diag.MPOrder',0) MPOrder = iMethfessel_Paxton diff --git a/tools/PostProcessing/process_module.f90 b/tools/PostProcessing/process_module.f90 index 483ad5cfd..779951c35 100644 --- a/tools/PostProcessing/process_module.f90 +++ b/tools/PostProcessing/process_module.f90 @@ -298,7 +298,8 @@ subroutine process_dos 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 global_module, ONLY: nspin, n_DOS, E_DOS_min, E_DOS_max, sigma_DOS, & + flag_fix_spin_population use units, ONLY: HaToeV implicit none @@ -326,7 +327,13 @@ subroutine process_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) + if(flag_procwf_range_Ef) then + if(nspin>1.and.flag_fix_spin_population) then + range_offset(i_spin) = efermi(1) + else + range_offset(i_spin) = efermi(i_spin) + end if + end if end do ! Set limits and broaden energy range if needed if(abs(E_DOS_min)** Added the optional filename + if ( species_from_files ) then + read (unit=input_array(block_start+i-1),fmt=*) & + j, mass(j), & + species_label(j), & + species_file(j) + else + read (unit=input_array(block_start+i-1),fmt=*) & + j, mass(j), & + species_label(j) + end if type_species(j)=j end do call fdf_endblock From f7bc3869f84fa0e8a3d62e658ae136495c79907b Mon Sep 17 00:00:00 2001 From: David Bowler Date: Thu, 6 Aug 2026 10:46:30 +0100 Subject: [PATCH 5/9] Update to documentation to describe gap output --- docs/groundstate.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/docs/groundstate.rst b/docs/groundstate.rst index 5fe1e8c23..88d44e46b 100644 --- a/docs/groundstate.rst +++ b/docs/groundstate.rst @@ -156,6 +156,9 @@ finite gap, you can force integer occupancies with the flag: This is false by default, and will lead to significant problems if the system is metallic at any point during the self-consistent ground state search, so should be considered carefully. +In both of these cases, if there is gap in the band structure, it is written out along with +the location in the Brillouin zone where it is found (though this result will depend very +sensitively on the k-point sampling used). The Methfessel-Paxton approach :cite:`g-Methfessel:1989ny` to occupations allows much higher smearing temperatures with minimal effect on the free energy (and From 8d7f3a967158647d27764ec61c12b5999a3a51c9 Mon Sep 17 00:00:00 2001 From: David Bowler Date: Thu, 20 Aug 2026 12:07:16 +0100 Subject: [PATCH 6/9] Tweaks to Methfessel-Paxton Mainly just making sure that the starting low electron number is sensible, and tweaking formatting --- src/DiagModule.f90 | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index 37c638b16..e38757022 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -319,7 +319,7 @@ module DiagModule ! that Methfessel-Paxton approximation may casue the bracket search ! algorithm to fail.) integer :: max_brkt_iterations - real(double), parameter :: tolElec = 1.0e-8_double + real(double), parameter :: tolElec = 1.0e-6_double contains @@ -2623,7 +2623,7 @@ subroutine findFermi_fixspin(electrons, eig, nbands, nkp, Ef, occ) ! Fill the bands for the first (electrons-NElec_less) electrons if (NElec_less >= electrons(spin)) then if (inode == ionode) write (io_lun, 7) myid, spin - NELec_less = electrons(spin) + NElec_less = electrons(spin)-one end if thisElec(spin) = zero band1: do iband = 1, nbands @@ -2729,13 +2729,13 @@ subroutine findFermi_fixspin(electrons, eig, nbands, nkp, Ef, occ) 5 format(10x, 'Proc: ', i5, ' findFermi_fixspin: level, Ne: ', 2f12.5) 6 format(10x, 'Proc: ', i5, ' findFermi_fixspin: found upper bound', f12.5) 7 format(10x, 'Proc: ', i5, ' findFermi_fixspin: Warning! Diag.NElecLess >= & - &total number of electrons for spin channel ', i2, & - ' setting it equal to number of electrons, but this is slow & + &total number of electrons for spin channel ', i2,/ & + 12x,' setting it equal to number of electrons, but this is slow & &and you may want to change it to something smaller.') 8 format(10x, 'Proc: ', i5, ' findFermi_fixspin: Warning! the & &calculated number of electrons (',f12.5, & - ') > electron_number (for spin ', i2, ' ) - 1.0. May be you & - &should increase the value of Diag.NElecLess (at the moment =& + ') >'/12x,'electron_number (for spin ', i2, ' ) - 1.0. May be you & + &should increase the value of Diag.NElecLess '/12x,'(at the moment =& & ',f12.5,')') 10 format(10x, 'Fermi level is ', f12.5) 11 format(10x, 'Fermi level for spin ', i2, ' is ', f12.5) @@ -2874,7 +2874,7 @@ subroutine findFermi_varspin(electrons_total, eig, nbands, nkp, Ef, occ) ! Fill the bands for the first (electrons_toal - NElec_less) electrons if (NElec_less >= electrons_total) then if (inode == ionode) write (io_lun, 6) - NELec_less = electrons_total + NELec_less = electrons_total - one end if thisElec = zero band1 : do iband = 1, nbands @@ -2892,7 +2892,7 @@ subroutine findFermi_varspin(electrons_total, eig, nbands, nkp, Ef, occ) lowElec = spin_factor * sum(electrons(:)) ! check if we indeed have a good lower bound if ((electrons_total - lowElec) < two) then - if (inode == ionode) write (io_lun, 8) lowElec, NElec_less + if (inode == ionode) write (io_lun, 8) lowElec ! find the lowest energy and start from there lband = 1 lkp = 1 @@ -3221,6 +3221,7 @@ subroutine write_gaps(spin_ch) write(io_lun,fmt='(4x,"CBM is at ",3f8.4," (1/a0)")') kk(:,cbm_k(spin)) end if end if + write(io_lun,fmt='(4x,"Fermi level= ",f12.5," ",a2)') en_conv*Efermi(spin),en_units(energy_units) end do end if return From f7e39ca834565f8f4f68f7be7dbc5aa6ae83a83c Mon Sep 17 00:00:00 2001 From: David Bowler Date: Thu, 20 Aug 2026 12:28:28 +0100 Subject: [PATCH 7/9] Further tweaks to MP format --- src/DiagModule.f90 | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index e38757022..e6980ec01 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -2892,7 +2892,7 @@ subroutine findFermi_varspin(electrons_total, eig, nbands, nkp, Ef, occ) lowElec = spin_factor * sum(electrons(:)) ! check if we indeed have a good lower bound if ((electrons_total - lowElec) < two) then - if (inode == ionode) write (io_lun, 8) lowElec + if (inode == ionode) write (io_lun, 7) myid, lowElec, NElec_less ! find the lowest energy and start from there lband = 1 lkp = 1 @@ -2975,9 +2975,13 @@ subroutine findFermi_varspin(electrons_total, eig, nbands, nkp, Ef, occ) 4 format(10x, 'Proc: ', i5, ' findFermi_varspin: found upper bound ', f12.5) 5 format(10x, 'Proc: ', i5, ' bracketed Ef: ', 2f12.5) 6 format(10x, 'In findFermi, Warning! Diag.NElecLess >= total number & - &of electrons, setting it equal to number of electrons, but & - &this is slow and you may want to change it to something & - &smaller.') + &of electrons, setting it equal to number of electrons, '/12x, & + 'but this is slow and you may want to change it to something smaller.') +7 format(10x, 'Proc: ', i5, ' findFermi_varspin: Warning! the & + &calculated number of electrons (',f12.5, & + ') > '/12x,'electron_number - 2.0. May be you & + &should increase the value of Diag.NElecLess (at the moment =& + & ',f12.5,')') 8 format(10x, 'Fermi level is ', f12.5) end subroutine findFermi_varspin From 4f608494f803dddbb9ff0d20a5dc905d87855c8a Mon Sep 17 00:00:00 2001 From: David Bowler Date: Thu, 20 Aug 2026 13:20:19 +0100 Subject: [PATCH 8/9] Bug fix: use tolElec when deciding on VBM vs CBM I had used very_small to decide whether an energy level was VBM or CBM (difference to one half) but this broke when tolElec was bigger than very_small. Only affected test007 but broke that comprehensively. --- src/DiagModule.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index e6980ec01..563df1a51 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -3086,7 +3086,7 @@ subroutine occupy(occu, ebands, Ef, electrons, nbands, nkp, spin) occu(iband,ikp,ss) = wtk(ikp) * locc ! Gap detection ! Check for VBM and allow for degenerate states - if(locc-half>-very_small) then!locc>=half) then + if(locc-half>-tolElec) then!locc>=half) then if(vbm(ss)-very_small) then!locc<=half) then + if(half-locc>-tolElec) then!locc<=half) then if(cbm(ss)>ebands(iband,ikp,ss)) then cbm_k(ss) = ikp cbm(ss) = ebands(iband,ikp,ss) From c8e90c051aca0211aef76a7b8aed9ca2776ff5ba Mon Sep 17 00:00:00 2001 From: David Bowler Date: Thu, 20 Aug 2026 13:46:36 +0100 Subject: [PATCH 9/9] Added flag to allow Fermi level adjustment to be turned off This is unlikely to be needed unless for debugging (test_007 made me think of it) but it's included for completeness. Small updates to documentation. --- docs/groundstate.rst | 3 ++- docs/input_tags.rst | 5 +++++ src/DiagModule.f90 | 11 ++++++++--- src/initial_read_module.f90 | 3 ++- 4 files changed, 17 insertions(+), 5 deletions(-) diff --git a/docs/groundstate.rst b/docs/groundstate.rst index 88d44e46b..360739535 100644 --- a/docs/groundstate.rst +++ b/docs/groundstate.rst @@ -158,7 +158,8 @@ This is false by default, and will lead to significant problems if the system is at any point during the self-consistent ground state search, so should be considered carefully. In both of these cases, if there is gap in the band structure, it is written out along with the location in the Brillouin zone where it is found (though this result will depend very -sensitively on the k-point sampling used). +sensitively on the k-point sampling used). The Fermi level will be placed at the mid-point +between the valence band maximum and the conduction band minimum. The Methfessel-Paxton approach :cite:`g-Methfessel:1989ny` to occupations allows much higher smearing temperatures with minimal effect on the free energy (and diff --git a/docs/input_tags.rst b/docs/input_tags.rst index 41946f908..074050995 100644 --- a/docs/input_tags.rst +++ b/docs/input_tags.rst @@ -738,6 +738,11 @@ Diag.IntegerOccs (*logical*) *default*: False +Diag.AdjustEf (*logical*) + Turns off the automatic adjustment of Ef to be mid-gap (unlikely to be needed). + + *default*: True + Diag.MPOrder (*integer*) Order of Bessel function approximation to delta-function used in Methfessel-Paxton smearing diff --git a/src/DiagModule.f90 b/src/DiagModule.f90 index 563df1a51..6ce787702 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -265,6 +265,7 @@ module DiagModule integer, dimension(2) :: band_ef logical :: flag_integer_occ logical, dimension(2) :: flag_gap + logical :: flag_adjust_Ef ! K-point data - here so that reading of k-points can take place in ! different routine to FindEvals @@ -618,7 +619,7 @@ subroutine FindEvals(electrons) if(flag_smear_type==0) then if(flag_fix_spin_population) then ! Separate fermi levels do spin=1,nspin - if(flag_gap(spin)) then + if(flag_gap(spin).and.flag_adjust_Ef) then Efermi(spin) = half*(cbm(spin)+vbm(spin)) locc(spin) = electrons(spin) if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps(spin_ch=spin) @@ -626,10 +627,12 @@ subroutine FindEvals(electrons) if(inode==ionode.and.iprint_DM + min_layer >= 2) & write(io_lun,'(4x, "Adjusted Fermi level for spin ", i2, " is ", f12.5)') & spin,Efermi(spin) + else if(flag_gap(spin)) then + if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps(spin_ch=spin) end if end do else - if(flag_gap(1) .or. flag_gap(2)) then ! Fermi levels should be same + if((flag_gap(1) .or. flag_gap(2)).and.flag_adjust_Ef) then ! Fermi levels should be same if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps ! Adjust Fermi level to lie at mid-gap and revisit occupancies do spin=1,nspin @@ -640,9 +643,11 @@ subroutine FindEvals(electrons) spin,Efermi(spin) end do call occupy(occ, evals, Efermi, locc, matrix_size, nkp) + else if(flag_gap(1).or.flag_gap(2)) then + if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps end if end if - else if(flag_integer_occ) then ! There has to be a gap! + else if(flag_integer_occ) then ! There has to be a gap; Ef set to mid-gap in occupy if(inode==ionode.and.iprint_DM + min_layer >= 2) call write_gaps end if ! Allocate space to expand eigenvectors into (i.e. when reversing diff --git a/src/initial_read_module.f90 b/src/initial_read_module.f90 index 4b4a623d2..922842214 100644 --- a/src/initial_read_module.f90 +++ b/src/initial_read_module.f90 @@ -3231,7 +3231,7 @@ subroutine readDiagInfo use DiagModule, only: nkp, kk, wtk, kT, maxefermi, & flag_smear_type, iMethfessel_Paxton, & max_brkt_iterations, gaussian_height, & - finess, NElec_less, flag_integer_occ + finess, NElec_less, flag_integer_occ, flag_adjust_Ef use energy, only: SmearingType, MPOrder use memory_module, only: reg_alloc_mem, reg_dealloc_mem, & type_dbl @@ -3275,6 +3275,7 @@ subroutine readDiagInfo call cq_abort("Cannot have free spin population with integer occupancies") flag_smear_type = -1 end if + flag_adjust_Ef = fdf_boolean('Diag.AdjustEf',.true.) SmearingType = flag_smear_type iMethfessel_Paxton = fdf_integer('Diag.MPOrder',0) MPOrder = iMethfessel_Paxton