Improved the accuracy of delta sum close to the edge of the Fermi window by including energy eigenvalues on-the-fly of the tetrahedral vertices with energies lying outside the chosen Fermi-window and added a test for it

fpm.toml

@@ -83,3 +83,8 @@ main = "check_interactions_symmetries.f90"
 name = "bte_regression"
 source-dir="test"
 main = "bte_regression.f90"
+
+[[test]]
+name = "fermi_window_adjustment"
+source-dir="test"
+main = "check_fermi_window_adjustment.f90"

src/delta.f90

@@ -29,7 +29,7 @@ module delta
   public form_tetrahedra_3d, fill_tetrahedra_3d, &
        form_triangles, fill_triangles, real_tetra, &
        delta_fn, get_delta_fn_pointer, &
-       delta_fn_triang, delta_fn_tetra
+       delta_fn_triang, delta_fn_tetra, fermi_window_adjusted_fill_tetrahedra_3d
 
   abstract interface
      pure real(r64) function delta_fn(e, ik, ib, mesh, map, count, evals)
@@ -586,11 +586,25 @@ subroutine form_tetrahedra_3d(nk, mesh, tetra, tetracount, tetramap, &
              kk = scvol_vertices(aux,3)
              aux = mux_vector([ii, jj, kk], mesh, 1_i64)
              tmp = aux !Guaranteed to be > 0
-             if(blocks) then
-                !Which point in indexlist does aux correspond to?
-                call binsearch(indexlist, aux, tmp) !tmp < 0 if search fails.
-             end if
-             tetra(count, tl) = tmp
+
+
+	    ! If energy-restricted blocks are being used, find the index in indexlist
+	    ! The updated snippet keeps track of vertices whose eigenvalues lie outside the chosen window around fermi energy
+
+	    if(blocks) then
+		call binsearch(indexlist, aux, tmp)   ! Binary search in indexlist
+
+
+		if (tmp < 0) then
+		    tetra(count, tl) = -aux            ! If vertex is outside Fermi window, save negative index (for the exceptional vertices)
+		else
+		    tetra(count, tl) = tmp             ! Save the index from indexlist
+		end if
+
+	    else
+	    tetra(count, tl) = aux                ! Save the multiplexed index
+	    end if
+
 
              if(tmp > 0) then
                 !Save the mapping of a wave vector index to a (tetrahedron, vertex)
@@ -644,6 +658,77 @@ subroutine fill_tetrahedra_3d(tetra, evals, tetra_evals)
        end do
     end do
   end subroutine fill_tetrahedra_3d
+
+
+subroutine fermi_window_adjusted_fill_tetrahedra_3d(wann, crys, wave_vector_mesh, tetra, evals, tetra_evals)
+
+    type(wannier), intent(in) :: wann
+    type(crystal), intent(in) :: crys
+    integer(i64), intent(in) :: tetra(:,:)
+    real(r64), intent(in) :: evals(:,:)
+    integer(i64), intent(in) :: wave_vector_mesh(1, 3)
+
+    real(r64), allocatable, intent(out) :: tetra_evals(:,:,:)
+
+    integer(i64) :: iv, it, ib, numbands, aux, numtetra
+    integer(i64) :: kpoint(3)
+
+    real(r64), allocatable :: energies(:,:)
+    real(r64) :: kvecs(1,3)
+
+    ! Calculate dimensions
+    numtetra = size(tetra, 1)
+    numbands = size(evals, 2)
+
+    print *, "Number of tetrahedra (numtetra) =", numtetra
+    print *, "Number of bands (numbands) =", numbands
+
+	  print *, "tetra shape:", shape(tetra)
+	  print *, "evals shape:", shape(evals)
+	  print *, "tetra_evals shape:", shape(tetra_evals)
+
+    ! Allocate arrays with correct dimensions
+   allocate(tetra_evals(numtetra, numbands, 4))
+   allocate(energies(1, wann%numwannbands))
+
+    ! Initialize tetra_evals and energies
+    tetra_evals(:,:,:) = 0.0_r64
+    energies(:, :) = 0.0_r64
+
+    ! Loop over all tetrahedra and vertices
+    do it = 1, numtetra
+       do iv = 1, 4
+          aux = tetra(it, iv)
+
+          !print *, "Accessing tetra array with index it=", it, ", iv=", iv
+
+	  !print *, "aux = ", aux
+
+          if (aux < 0) then
+             ! Vertex outside Fermi window, calculate on-the-fly
+             call demux_vector(-aux, kpoint, wave_vector_mesh, 1_i64)
+             kvecs(1, :) = real(kpoint) / real(wave_vector_mesh(1, :))
+
+             call wann%el_wann(crys, 1_i64, kvecs, energies)
+             tetra_evals(it, :, iv) = energies(1, :)
+          else
+             ! Use pre-calculated eigenvalue
+             tetra_evals(it, :, iv) = evals(aux, :)
+          end if
+       end do
+    end do
+
+    ! Sort eigenvalues for each tetrahedron vertex
+    do it = 1, numtetra
+       do ib = 1, numbands
+          call sort(tetra_evals(it, ib, :))
+       end do
+    end do
+    deallocate(energies)
+end subroutine fermi_window_adjusted_fill_tetrahedra_3d
+
+
+
 
   subroutine form_triangles(nk, mesh, triang, triangcount, triangmap, &
        blocks, indexlist)

src/electron.f90

@@ -25,7 +25,7 @@ module electron_module
   use wannier_module, only: wannier
   use crystal_module, only: crystal, calculate_wavevectors_full
   use symmetry_module, only: symmetry, find_irred_wedge, create_fbz2ibz_map
-  use delta, only: form_tetrahedra_3d, fill_tetrahedra_3d, form_triangles, &
+  use delta, only: form_tetrahedra_3d, fermi_window_adjusted_fill_tetrahedra_3d, form_triangles, &
        fill_triangles
 
   implicit none
@@ -584,7 +584,9 @@ subroutine calculate_electrons(self, wann, crys, sym, num)
        call form_tetrahedra_3d(self%nwv, self%wvmesh, &
             self%simplicial_complex, self%simplex_count, &
             self%simplex_map, .true., self%indexlist)
-       call fill_tetrahedra_3d(self%simplicial_complex, self%ens, self%simplex_evals)
+       ! call fill_tetrahedra_3d(self%simplicial_complex, self%ens, self%simplex_evals)
+       call fermi_window_adjusted_fill_tetrahedra_3d(wann, crys, self%wvmesh, &             ! --> here's the fermi-window adjusted call
+       	    self%simplicial_complex, self%ens, self%simplex_evals)
     else
        call print_message("Calculating electron mesh triangles...")
        call form_triangles(self%nwv, self%wvmesh, &

test/3C-SiC/fpm_run_fermi_window_adj_caf.sh

@@ -0,0 +1,23 @@
+#!/bin/bash
+
+# Download data from public repo
+curl -X GET "https://nomad-lab.eu/prod/v1/api/v1/uploads/5b-c4JYRSG6qitMFjFUUzQ/raw/?offset=0&length=-1&compress=true" -o 3C-SiC.zip
+
+#Extract data
+unzip 3C-SiC.zip
+rm 3C-SiC.zip
+
+workdir="./3C-SiC_test_output"
+datadir="./3C-SiC/elphbolt_input_data"
+inputdir="./test/3C-SiC"
+
+mkdir -p $workdir
+cp $datadir/* $workdir
+cp $inputdir/input.nml $workdir
+rm -r 3C-SiC
+cd $workdir
+
+#gcc+opencoarrays
+cafrun -n 2 ../build/caf_*/test/fermi_window_adjustment | tee 3C_SiC_test.output
+
+cd ..

test/check_fermi_window_adjustment.f90

@@ -0,0 +1,125 @@
+program check_fermi_window_adjustment
+
+  use iso_fortran_env, only : r64 => real64, i64 => int64
+  use testify_m, only : testify
+  use misc, only: print_message, subtitle, timer, exit_with_message
+  use numerics_module, only: numerics
+  use crystal_module, only: crystal
+  use symmetry_module, only: symmetry
+  use electron_module, only: electron
+  use phonon_module, only: phonon
+  use wannier_module, only: wannier
+  use bte_module, only: bte
+  use bz_sums, only: calculate_dos, calculate_qTF, calculate_el_dos_fermi, calculate_el_Ws
+  use interactions, only: calculate_gReq, calculate_gkRp, calculate_3ph_interaction, &
+       calculate_eph_interaction_ibzq, calculate_eph_interaction_ibzk, &
+       calculate_echimp_interaction_ibzk, calculate_bound_scatt_rates
+
+  implicit none
+
+  integer :: itest
+  integer, parameter :: num_tests = 4
+  type(testify) :: test_array(num_tests), tests_all
+
+  type(numerics) :: num
+  type(crystal) :: crys
+  type(symmetry) :: sym
+  type(wannier) :: wann
+  type(electron) :: el
+  type(phonon) :: ph
+  type(bte) :: bt
+  type(timer) :: t_all, t_event
+
+  !character(:), allocatable :: datalink, curl_arg, workdir, datadir, inputdir
+
+  if(this_image() == 1) then
+     write(*, '(A)')  'Regression test on 3C-SiC'
+     write(*, '(A, I5)') 'Number of coarray images = ', num_images()
+  end if
+
+  !Test counter
+  itest = 0
+
+  call t_all%start_timer('elphbolt: Fermi-window adjustment')
+
+  call t_event%start_timer('Initialization')
+
+  !Set up crystal
+  call crys%initialize
+
+  !Set up numerics data
+  call num%initialize(crys)
+
+  !Calculate crystal and BZ symmetries
+  call sym%calculate_symmetries(crys, num%qmesh)
+
+  !Test symmetries
+  if(this_image() == 1) then
+     itest = itest + 1
+     test_array(itest) = testify("number of symmetries")
+     call test_array(itest)%assert(sym%nsymm, 24_i64)
+
+     itest = itest + 1
+     test_array(itest) = testify("symmetry group")
+     call test_array(itest)%assert(sym%international, "F-43m")
+  end if
+  sync all
+  !!
+
+  if(num%onlyebte .or. num%drag .or. num%phe &
+       .or. num%plot_along_path .or. num%runlevel == 3) then
+     !Read EPW Wannier data
+     call wann%read(num)
+
+     !Calculate electrons
+     call el%initialize(wann, crys, sym, num)
+
+     !Test electron energies
+     if(this_image() == 1) then
+        itest = itest + 1
+        test_array(itest) = testify("electron energies")
+        call test_array(itest)%assert( &
+             [transpose(el%ens_irred)], &
+             [0.1091276329E+02_r64, 0.1459479591E+02_r64, 0.2329596906E+02_r64, 0.2329596906E+02_r64, &
+              0.1096209018E+02_r64, 0.1442980731E+02_r64, 0.2334947484E+02_r64, 0.2340698517E+02_r64, &
+              0.1079030242E+02_r64, 0.1369997715E+02_r64, 0.2354625098E+02_r64, 0.2354625098E+02_r64, &
+              0.1107938027E+02_r64, 0.1467078718E+02_r64, 0.2311065489E+02_r64, 0.2358101936E+02_r64], &
+             tol = 1.0e-8_r64)
+     end if
+     !!
+  end if
+
+  call t_event%end_timer('Initialization')
+
+  call t_event%start_timer('Density of states and one-particle scattering rates')
+
+  call subtitle("Calculating density of states...")
+  if(num%onlyebte .or. num%drag) then
+     !Calculate electron density of states
+     call calculate_dos(el, num%tetrahedra)
+
+
+     !Test electron density of states
+     if(this_image() == 1) then
+        itest = itest + 1
+        test_array(itest) = testify("electron DOS")
+        call test_array(itest)%assert( &
+             [transpose(el%dos)], &
+             [0.6944753689E-02_r64, 0.7200285677E-02_r64, 0.2139544764E-01_r64, 0.2139544797E-01_r64, &
+              0.2264851372E-01_r64, 0.3001614228E-02_r64, 0.1790594732E-01_r64, 0.2216848884E-01_r64, &
+              0.0000000000E+00_r64, 0.0000000000E+00_r64, 0.6902017130E-02_r64, 0.6902016684E-02_r64, &
+              0.9485621719E-02_r64, 0.2888870927E-02_r64, 0.9858175722E-02_r64, 0.1918518823E-02_r64], &
+             tol = 1.0e-11_r64)
+     end if
+    end if
+     !!
+
+  if(this_image() == 1) then
+     tests_all = testify(test_array)
+     call tests_all%report
+
+     if(tests_all%get_status() .eqv. .false.) error stop -1
+  end if
+  sync all
+  call t_all%end_timer('elphbolt: Fermi-window adjustment')
+end program check_fermi_window_adjustment