diff --git a/docs/groundstate.rst b/docs/groundstate.rst index 354f8db4..36073953 100644 --- a/docs/groundstate.rst +++ b/docs/groundstate.rst @@ -144,10 +144,22 @@ 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. +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 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 62479fd5..07405099 100644 --- a/docs/input_tags.rst +++ b/docs/input_tags.rst @@ -733,6 +733,16 @@ Diag.kT (*real*) *default*: 0.001 +Diag.IntegerOccs (*logical*) + Specifies integer occupation of eigenstates (overrides ``Diag.SmearingType``) + + *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 a8af360f..6ce78770 100644 --- a/src/DiagModule.f90 +++ b/src/DiagModule.f90 @@ -260,9 +260,13 @@ 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 + 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 integer :: nkp @@ -316,6 +320,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-6_double contains @@ -602,7 +607,49 @@ 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 + flag_gap = .false. call findFermi(electrons, evals, matrix_size, nkp, Efermi, occ) + ! Test for gap + if(flag_smear_type==0) then + if(flag_fix_spin_population) then ! Separate fermi levels + do spin=1,nspin + 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) + 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) + 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)).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 + 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) + 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; 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 ! ScaLAPACK distribution) allocate(expH(matrix_size,prim_size,nspin), STAT=stat) @@ -928,6 +975,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 @@ -2362,10 +2411,37 @@ 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)two*kT.and.(cbm(nspin)-vbm(nspin)>two*kT)) flag_gap = .true. + return + end if ! Integer occupations if (nspin == 2) then electrons_total = electrons(1) + electrons(2) else @@ -2373,8 +2449,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 @@ -2446,7 +2526,6 @@ subroutine findFermi_fixspin(electrons, eig, nbands, nkp, Ef, occ) integer, dimension(nspin) :: ne real(double) :: electrons_total, gaussian_width integer :: counter, ibrkt, lband, lkp, iband, ikp, spin - real(double), parameter :: tolElec = 1.0e-6_double if (nspin == 2) then electrons_total = electrons(1) + electrons(2) @@ -2549,7 +2628,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 @@ -2655,13 +2734,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) @@ -2718,7 +2797,6 @@ subroutine findFermi_varspin(electrons_total, eig, nbands, nkp, Ef, occ) real(double), dimension(nspin) :: lowEf, highEf, incEf real(double) :: gaussian_width, thisElec, lowElec, highElec integer :: counter, ne, ibrkt, lband, lkp, iband, ikp, spin, lspin - real(double), parameter :: tolElec = 1.0e-6_double ! Finding the correct bracket trapping Ef select case (flag_smear_type) @@ -2801,7 +2879,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 @@ -2819,7 +2897,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, 7) myid, lowElec, NElec_less ! find the lowest energy and start from there lband = 1 lkp = 1 @@ -2902,9 +2980,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 @@ -2975,8 +3057,9 @@ 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 if (nspin == 2 .and. present(spin)) then ss_start = spin @@ -2985,19 +3068,63 @@ 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 + cbm(ss) = BIG + gap_d(ss) = BIG kp: do ikp = 1, nkp band: do iband = 1, nbands select case (flag_smear_type) case (0) ! Fermi smearing - occu(iband,ikp,ss) = & - wtk(ikp) * fermi(ebands(iband,ikp,ss) - Ef(ss), kT) + locc = fermi(ebands(iband,ikp,ss) - Ef(ss), kT) + occu(iband,ikp,ss) = wtk(ikp) * locc + ! Gap detection + ! Check for VBM and allow for degenerate states + if(locc-half>-tolElec) then!locc>=half) then + if(vbm(ss)-tolElec) then!locc<=half) then + if(cbm(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) @@ -3035,6 +3162,81 @@ end subroutine occupy !!*** + ! ----------------------------------------------------------------------------- + ! Subroutine write_gaps + ! ----------------------------------------------------------------------------- + + !!****f* DiagModule/write_gaps * + !! + !! NAME + !! write_gaps + !! USAGE + !! write_gaps + !! PURPOSE + !! Writes out gaps (when there is one) + !! + !! INPUTS + !! + !! USES + !! units, global + !! AUTHOR + !! D.R.Bowler + !! CREATION DATE + !! 31/07/2026 + !! MODIFICATION HISTORY + !! + !! SOURCE + subroutine write_gaps(spin_ch) + + use units + use global_module, only: nspin + + implicit none + + ! Passed variables + integer, optional :: spin_ch + + ! 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 + write(io_lun,fmt='(4x,"Fermi level= ",f12.5," ",a2)') en_conv*Efermi(spin),en_units(energy_units) + 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 08b6065e..92284221 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 @@ -3228,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 + 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 @@ -3266,6 +3269,13 @@ subroutine readDiagInfo 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 + flag_adjust_Ef = fdf_boolean('Diag.AdjustEf',.true.) SmearingType = flag_smear_type iMethfessel_Paxton = fdf_integer('Diag.MPOrder',0) MPOrder = iMethfessel_Paxton @@ -3760,7 +3770,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 d119869c..ffdb9f5a 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>$ diff --git a/src/pao2blip.f90 b/src/pao2blip.f90 index 57ae0651..00ce0e16 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 diff --git a/testsuite/test_check_output.py b/testsuite/test_check_output.py index 3279b25d..796acc49 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(): diff --git a/tools/PostProcessing/process_module.f90 b/tools/PostProcessing/process_module.f90 index 483ad5cf..779951c3 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