Tetrahedron method to calculate spin Hall conductivity

src/Makefile.2

@@ -171,7 +171,10 @@ kpath.o: $(POSTDIR)kpath.F90 comms.o types.o constants.o io.o spin.o berry.o pos
 kslice.o: $(POSTDIR)kslice.F90 comms.o types.o constants.o io.o spin.o berry.o postw90_types.o ws_distance.o
 	$(COMPILER) $(POSTOPTS) $(FCOPTS) -c $(POSTDIR)kslice.F90
 
-berry.o: $(POSTDIR)berry.F90 comms.o types.o constants.o utility.o postw90_common.o get_oper.o io.o spin.o wan_ham.o postw90_types.o ws_distance.o
+tetrahedron.o: $(POSTDIR)tetrahedron.F90 constants.o utility.o
+	$(COMPILER) $(POSTOPTS) $(FCOPTS) -c $(POSTDIR)tetrahedron.F90
+
+berry.o: $(POSTDIR)berry.F90 comms.o types.o constants.o utility.o postw90_common.o get_oper.o io.o spin.o wan_ham.o postw90_types.o ws_distance.o tetrahedron.o
 	$(COMPILER) $(POSTOPTS) $(FCOPTS) -c $(POSTDIR)berry.F90
 
 gyrotropic.o: $(POSTDIR)gyrotropic.F90 comms.o types.o constants.o utility.o postw90_common.o get_oper.o io.o spin.o wan_ham.o berry.o postw90_types.o ws_distance.o

src/Makefile.header

@@ -5,7 +5,7 @@ OBJS  =  constants.o io.o utility.o types.o hamiltonian.o overlap.o kmesh.o dise
 OBJS_LIBV1 = wannier_lib.o
 
 OBJS_POST  = ws_distance.o types.o kmesh.o io.o comms.o utility.o get_oper.o constants.o \
-        postw90_common.o wan_ham.o spin.o dos.o berry.o gyrotropic.o kpath.o kslice.o \
+        postw90_common.o wan_ham.o spin.o dos.o tetrahedron.o berry.o gyrotropic.o kpath.o kslice.o \
         boltzwann.o geninterp.o postw90_types.o postw90_readwrite.o readwrite.o \
 	error.o error_base.o
 

src/postw90/postw90_common.F90

@@ -506,6 +506,12 @@ subroutine pw90common_wanint_w90_wannier90_readwrite_dist(print_output, ws_regio
     if (allocated(error)) return
     call comms_bcast(pw90_berry%sc_use_eta_corr, 1, error, comm)
     if (allocated(error)) return
+    call comms_bcast(pw90_berry%tetrahedron_method, 1, error, comm)
+    if (allocated(error)) return
+    call comms_bcast(pw90_berry%tetrahedron_higher_correction, 1, error, comm)
+    if (allocated(error)) return
+    call comms_bcast(pw90_berry%tetrahedron_cutoff, 1, error, comm)
+    if (allocated(error)) return
 ! ----------------------------------------------
 !
 ! New input variables in development

src/postw90/postw90_readwrite.F90

@@ -836,6 +836,24 @@ subroutine w90_wannier90_readwrite_read_berry(pw90_calculation, pw90_berry, pw90
                                           r_value=pw90_berry%kdotp_kpoint)
     if (allocated(error)) return
 
+    pw90_berry%tetrahedron_method = .false.
+    call w90_readwrite_get_keyword('tetrahedron_method', found, error, comm, &
+                                   l_value=pw90_berry%tetrahedron_method)
+    if (pw90_berry%tetrahedron_method .and. index(pw90_berry%task, 'shc') == 0) &
+      call set_error_input(error, 'Error: tetrahedron_method only implemented in the shc routine', comm)
+
+    pw90_berry%tetrahedron_higher_correction = .true.
+    call w90_readwrite_get_keyword('tetrahedron_higher_correction', found, error, comm, &
+                                   l_value=pw90_berry%tetrahedron_method)
+    if (.not. pw90_berry%tetrahedron_higher_correction) call set_error_input &
+      (error, 'Error: Set tetrahedron_higher_correction = .true., tetrahedron_method works only with correction', comm)
+
+    pw90_berry%tetrahedron_cutoff = 1.e-4_dp
+    call w90_readwrite_get_keyword('tetrahedron_cutoff', found, error, comm, &
+                                   r_value=pw90_berry%tetrahedron_cutoff)
+    if (pw90_berry%tetrahedron_cutoff <= 0._dp) call set_error_input &
+      (error, 'Error: tetrahedron_cutoff must be greater than zero', comm)
+
     kdotp_num_bands = 0
     call w90_readwrite_get_keyword('kdotp_num_bands', found, error, comm, i_value=kdotp_num_bands)
     if (allocated(error)) return
@@ -2055,6 +2073,15 @@ subroutine w90_postw90_readwrite_write(print_output, w90_system, fermi_energy_li
         write (stdout, '(1x,a21,5x,a47,4x,a1)') '|  Smearing Function ', &
           trim(w90_readwrite_get_smearing_type(pw90_dos%smearing%type_index)), '|'
       endif
+      if (pw90_berry%tetrahedron_method) then
+        write (stdout, '(1x,a46,10x,a8,13x,a1)') '|  Tetrahedron method (PRB 106, 075125)                        :', '       T', '|'
+        write (stdout, '(1x,a46,10x,a8,13x,L8)') '|  with higher-order correction                                :', &
+          pw90_berry%tetrahedron_higher_correction, '|'
+        write (stdout, '(1x,a46,10x,E8.3,13x,a1)') '|  Tetrahedron cutoff                        :', &
+          pw90_berry%tetrahedron_cutoff, '|'
+      else
+        write (stdout, '(1x,a46,10x,a8,13x,a1)') '|  Tetrahedron method                        :', '       F', '|'
+      endif
       if (pw90_extra_io%global_kmesh%mesh(1) == pw90_dos%kmesh%mesh(1) .and. &
           pw90_extra_io%global_kmesh%mesh(2) == pw90_dos%kmesh%mesh(2) .and. &
           pw90_extra_io%global_kmesh%mesh(3) == pw90_dos%kmesh%mesh(3)) then

src/postw90/postw90_types.F90

@@ -165,8 +165,19 @@ module w90_postw90_types
     integer :: kubo_nfreq
     complex(kind=dp), allocatable :: kubo_freq_list(:)
     real(kind=dp) :: kubo_eigval_max
+    ! + Module tetrahedron
+    !type(pw90_tetrahedron_type) :: tetrahedron
+    logical :: tetrahedron_method
+    logical :: tetrahedron_higher_correction
+    real(kind=dp) :: tetrahedron_cutoff
   end type pw90_berry_mod_type
 
+  !type pw90_tetrahedron_type
+  !  logical :: method
+  !  logical :: correction
+  !  real(kind=dp) :: cutoff
+  !end type pw90_tetrahedron_type
+
   ! spin Hall conductivity (postw90 - common, get_oper, berry, kpath)
   type pw90_spin_hall_type
     !!==================================================