Higher-order finite difference

docs/docs/user_guide/wannier90/parameters.md

@@ -77,6 +77,8 @@ and the the `wannier90` Tutorial.
 |   skip_B1_tests   |  L   | Check the condition B1 of Ref [@marzari-prb97].                                                   |
 |      nnkpts       |  I   | Explicit list of nearest-neighbour k-points.                                     |
 |     kmesh_tol     |  R   | The tolerance to control if two kpoint belong to the same shell                  |
+|   higher_order_n  |  I   | The order of higher-order finite difference to get b-vectors and weights         |
+| higher_order_nearest_shells |  L   | Use the b-vectors on the nearest shells                                |
 <!-- markdownlint-enable MD013 -->
 
 `seedname.win` file keywords defining the system. Argument types are
@@ -496,6 +498,19 @@ less than `kmesh_tol`. Units are Ang.
 
 The default value is 0.000001 Ang.
 
+### `integer :: higher_order_n`
+
+Specifies the order of finite-difference(FD) formula. The default is 1.
+Assuming the first-order FD uses #1, #3 shells for neighbors,
+then the n-th-order FD also includes 2, 3, ..., n times #1, #3 shells.
+A very high value of n may require larger `search_shells`.
+
+### `logical :: higher_order_nearest_shells`
+
+Instead of 2, 3, ..., n times the first-order shells, search the
+nearest shells from the origin to satisfy the higher-order version of
+the condition mentioned above. Not extensively tested.
+
 ## Projection
 
 The projections block defines a set of localised functions used to

src/library_interface.F90

@@ -809,7 +809,7 @@ subroutine w90_set_option_real(common_data, keyword, rval)
 
   subroutine w90_create_kmesh(common_data, istdout, istderr, ierr)
     use w90_error_base, only: w90_error_type
-    use w90_kmesh, only: kmesh_get
+    use w90_kmesh, only: kmesh_get, kmesh_sort
 
     implicit none
 
@@ -826,6 +826,11 @@ subroutine w90_create_kmesh(common_data, istdout, istderr, ierr)
     call kmesh_get(common_data%kmesh_input, common_data%kmesh_info, common_data%print_output, &
                    common_data%kpt_latt, common_data%real_lattice, common_data%num_kpts, &
                    common_data%gamma_only, istdout, common_data%timer, error, common_data%comm)
+
+    if (.NOT. common_data%gamma_only) then
+      call kmesh_sort(common_data%kmesh_info, common_data%num_kpts, error, common_data%comm)
+    endif
+
     if (allocated(error)) then
       call prterr(error, ierr, istdout, istderr, common_data%comm)
       return

src/readwrite.F90

@@ -729,6 +729,7 @@ subroutine w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
 
     integer :: itmp, ierr
     logical :: found
+    integer :: n
 
     call w90_readwrite_get_keyword(settings, 'search_shells', found, error, comm, &
                                    i_value=kmesh_input%search_shells)
@@ -737,6 +738,32 @@ subroutine w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
       call set_error_input(error, 'Error: search_shells must be positive', comm)
       return
     endif
+    call w90_readwrite_get_keyword(settings, 'search_supcell_size', found, error, comm, &
+                                   i_value=kmesh_input%search_supcell_size)
+    if (allocated(error)) return
+    if (kmesh_input%search_supcell_size < 0) then
+      call set_error_input(error, 'Error: search_supcell_size must be positive', comm)
+      return
+    endif
+    call w90_readwrite_get_keyword(settings, 'higher_order_n', found, error, comm, &
+                                   i_value=kmesh_input%higher_order_n)
+    if (allocated(error)) return
+    if (kmesh_input%higher_order_n < 0) then
+      call set_error_input(error, 'Error: higher_order_n must be positive', comm)
+      return
+    endif
+
+    n = kmesh_input%higher_order_n
+    kmesh_input%max_shells_h = n*(4*n**2 + 15*n + 17)/6
+    kmesh_input%max_shells_aux = kmesh_input%max_shells_h
+    kmesh_input%num_nnmax_h = 2*kmesh_input%max_shells_h
+
+    call w90_readwrite_get_keyword(settings, 'higher_order_nearest_shells', found, error, comm, &
+                                   l_value=kmesh_input%higher_order_nearest_shells)
+    if (allocated(error)) return
+    if (.not. kmesh_input%higher_order_nearest_shells) then
+      kmesh_input%max_shells_aux = 6
+    endif
 
     call w90_readwrite_get_keyword(settings, 'kmesh_tol', found, error, comm, &
                                    r_value=kmesh_input%tol)
@@ -750,8 +777,8 @@ subroutine w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
                                         .true., error, comm)
     if (allocated(error)) return
     if (found) then
-      if (kmesh_input%num_shells < 0 .or. kmesh_input%num_shells > max_shells) then
-        call set_error_input(error, 'Error: number of shell in shell_list must be between zero and six', comm)
+      if (kmesh_input%num_shells < 0 .or. kmesh_input%num_shells > kmesh_input%max_shells_h) then
+        call set_error_input(error, 'Error: number of shell in shell_list must be between zero and kmesh_input%max_shells_h', comm)
         return
       endif
       !if (allocated(kmesh_input%shell_list)) deallocate (kmesh_input%shell_list)
@@ -770,7 +797,7 @@ subroutine w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
     else
       !if (allocated(kmesh_input%shell_list)) deallocate (kmesh_input%shell_list)
       ! this is the default allocation of the shell_list--used by kmesh_shell_automatic()
-      allocate (kmesh_input%shell_list(max_shells), stat=ierr)
+      allocate (kmesh_input%shell_list(kmesh_input%max_shells_h), stat=ierr)
       if (ierr /= 0) then
         call set_error_alloc(error, 'Error allocating shell_list in w90_readwrite_read_kmesh_data', comm)
         return
@@ -986,6 +1013,8 @@ subroutine w90_readwrite_clear_keywords(settings, comm)
     call w90_readwrite_get_keyword(settings, 'fixed_step', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'gamma_only', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'guiding_centres', found, error, comm)
+    call w90_readwrite_get_keyword(settings, 'higher_order_n', found, error, comm)
+    call w90_readwrite_get_keyword(settings, 'higher_order_nearest_shells', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'hr_cutoff', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'hr_plot', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'iprint', found, error, comm)
@@ -1003,6 +1032,7 @@ subroutine w90_readwrite_clear_keywords(settings, comm)
     call w90_readwrite_get_keyword(settings, 'num_shells', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'num_valence_bands', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'num_wann', found, error, comm)
+    call w90_readwrite_get_keyword(settings, 'use_ss_functional', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'one_dim_axis', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'optimisation', found, error, comm)
     call w90_readwrite_get_keyword(settings, 'postproc_setup', found, error, comm)

src/types.F90

@@ -121,10 +121,18 @@ module w90_types
     !!==================================================
     integer :: num_shells = 0
     !! no longer an input keyword
+    !! higher-order finite difference
+    integer :: max_shells_h ! \sum ((3, 2i)) (combination with rep.) = n(4n*2 + 15n + 17)/6
+    integer :: max_shells_aux ! = 6 for higher-order simple algorithm, n(4n*2 + 15n + 17)/6 for search algorithm
+    integer :: num_nnmax_h ! 2*n(4n*2 + 15n + 17)/6
+    integer :: higher_order_n = 1
+    logical :: higher_order_nearest_shells = .false. ! experimental feature
+    !! A simpler algorithm of determining bvectors
     logical :: skip_B1_tests = .false.
     !! do not check the B1 condition
     integer, allocatable :: shell_list(:)
     integer :: search_shells = 36
+    integer :: search_supcell_size = 5 !Size of supercell (of recip cell) in which to search for k-point shells
     real(kind=dp) :: tol = 0.000001_dp
   end type kmesh_input_type