Merge pull request #533 from sjhong6230/wigner-seitz

Translationally-invariant computation for position-related matrix elements

docs/docs/parameters/postw90-global-parameters.csv

@@ -16,3 +16,5 @@ spin_moment,L,"Determines whether to evaluate the spin magnetic moment per cell"
 uHu_formatted,L,"Read a formatted `seedname.uHu` file"
 spn_formatted,L,"Read a formatted `seedname.spn` file"
 berry_curv_unit,S,"Unit of Berry curvature"
+transl_inv,L,"Use Marzari-Vanderbilt formula for wannier centres"
+transl_inv_full,L,"Use translationally-invariant formula for position-related matrix elements"

docs/docs/tutorials/tutorial_36.md

@@ -1,14 +1,12 @@
-# 34: Silicon — Valence and low-lying conduction states
+# 36: Silicon — Valence and low-lying conduction states
 
 ## Marzari-Vanderbilt functional
 
 - Outline: *Obtain MLWFs for the valence and low-lying conduction-band
     states of Si by minimizing the Marzari-Vanderbilt functional.
     Plot the interpolated bandstructure.*
 
-- Directory: `tutorials/tutorial34/Marzari-Vanderbilt`
-    *Files can be downloaded from
-    [here](https://github.com/wannier-developers/wannier90/tree/develop/tutorials/tutorial34)*
+- Directory: [`tutorial/tutorial36/`](https://github.com/wannier-developers/wannier90/tree/develop/tutorials/tutorial36)
 
 - Input Files
 
@@ -86,9 +84,7 @@
     states of Si by minimizing the Stengel-Spaldin functional.
     Plot the interpolated bandstructure.*
 
-- Directory: `tutorials/tutorial34/Stengel-Spaldin`
-  *Files can be downloaded from [here]
-  (<https://github.com/wannier-developers/wannier90/tree/develop/tutorials/tutorial34>)*
+- Directory: [`tutorial/tutorial36/`](https://github.com/wannier-developers/wannier90/tree/develop/tutorials/tutorial36)
 
 - Input Files
 

docs/docs/tutorials/tutorial_37.md

@@ -0,0 +1,70 @@
+# 37: Bcc Iron — Translationally-invariant Wannier Interpolation
+
+- Outline: *Perform Wannier interpolation of orbital magnetization
+    using translationally-invariant formulas.*
+
+- Directory: [`tutorial/tutorial37/`](https://github.com/wannier-developers/wannier90/tree/develop/tutorials/tutorial37)
+
+- Input Files
+
+    - `Fe.scf` *The `pwscf` input file for ground
+        state calculation*
+
+    - `Fe.nscf` *The `pwscf` input file to obtain
+        Bloch states on a uniform grid*
+
+    - `Fe.pw2wan` *Input file for `pw2wannier90`*
+
+    - `Fe.win` *The `wannier90` input file*
+
+1. For both directories, run `pwscf` to obtain the ground state of silicon
+
+    ```bash title="Terminal"
+    pw.x < Fe.scf > Fe.out
+    ```
+
+2. Run `pwscf` to obtain the Bloch states on a uniform
+    k-point grid.
+
+    ```bash title="Terminal"
+    pw.x < Fe.nscf > nscf.out
+    ```
+
+3. Run `wannier90` to generate a list of the required overlaps (written
+    into the `Fe.nnkp` file).
+
+    ```bash title="Terminal"
+    wannier90.x -pp Fe
+    ```
+
+4. Run `pw2wannier90` to compute the overlap between Bloch states and
+    the projections for the starting guess (written in the `Fe.mmn`,
+    `Fe.amn`, and `Fe.uHu` files).
+
+    ```bash title="Terminal"
+    pw2wannier90.x < Fe.pw2wan > Fe.out
+    ```
+
+5. Run `wannier90` to compute the MLWFs.
+
+    ```bash title="Terminal"
+    wannier90.x Fe
+    ```
+
+6. Run `postw90` to perform Wannier interpolation of orbital magnetization
+   with and without translationally-invariant formulas. You can use the option
+   `transl_inv_full`.
+
+    ```bash title="Terminal"
+    postw90.x Fe
+    ```
+
+7. Do step 6 varying the fermi energy and obtain a graph for the relation
+   between orbital magnetization and the fermi energy. Compare the results of
+   `without_translation` and `with_translation` directories for both
+   `transl_inv_full = F` and `transl_inv_full = T`.
+
+## Further ideas
+
+- Increase the grid and find out that the `transl_inv_full = T` case
+  converges faster.

docs/docs/user_guide/postw90/postw90params.md

@@ -421,6 +421,72 @@ correction, the convention of the Bloch sum (determined by
 
 The default value is `true`.
 
+### `logical :: transl_inv`
+
+If `transl_inv=true`, the Marzari-Vanderbilt formula is used
+for diagonal part of position matrix elements, or equivalently,
+for the wannier centres
+
+$$
+\mathbf{r}_i = - \frac{1}{N_{\mathbf{k}}}
+\sum_{\mathbf{k},\mathbf{b}} c_{\mathbf{b}}
+\mathbf{b} \, \text{Im} \ln M_{ii}^{(\mathbf{k,b})}
+$$
+
+instead of the naive finite difference formula.
+
+$$
+\mathbf{r}_i = \frac{i}{N_{\mathbf{k}}}
+\sum_{\mathbf{k},\mathbf{b}} c_{\mathbf{b}} \mathbf{b}
+\, M_{ii}^{(\mathbf{k,b})}
+$$
+
+where $M_{ij}^{(\mathbf{k,b})}$ is the overlap matrix in the Wannier gauge.
+
+$$
+M_{ij}^{(\mathbf{k,b})} = \sum_{m,n}
+U_{im}^{\dagger (\mathbf{k})} \langle u_{m\mathbf{k}}^{(\text{H})}
+| u_{n\mathbf{k+b}}^{(\text{H})} \rangle U_{nj}^{(\mathbf{k+b})}
+$$
+
+The default value is `false`.
+
+### `logical :: transl_inv_full`
+
+If `transl_inv_full=true`, the translationally-invariant finite difference
+formula is used for the computation of position and related matrix elements.
+For example, the naive finite difference formula used in the computation of
+off-diagonal position matrix elements is
+
+$$
+\mathbf{A}_{ij;\mathbf{R}} = \langle w_{i\mathbf{0}} | \hat{\mathbf{r}}
+| w_{j\mathbf{R}} \rangle = \frac{i}{N_{\mathbf{k}}} \sum_{\mathbf{k},
+\mathbf{b}} c_{\mathbf{b}} \mathbf{b} \, e^{-i\mathbf{k} \cdot \mathbf{R}}
+M_{ij}^{(\mathbf{k,b})} \, .
+$$
+
+For the `transl_inv_full=true`, the following formula is used
+
+$$
+\mathbf{A}_{ij;\mathbf{R}} = \langle w_{i\mathbf{0}} | \hat{\mathbf{r}}
+| w_{j\mathbf{R}} \rangle = \frac{i}{N_{\mathbf{k}}} \sum_{\mathbf{b}}
+c_{\mathbf{b}} \mathbf{b} \, e^{i\mathbf{b} \cdot
+\bar{\mathbf{r}}_{ij;\mathbf{R}}} \sum_\mathbf{k} e^{-i\mathbf{k}
+\cdot \mathbf{R}} M_{ij}^{(\mathbf{k,b})} +
+\bar{\mathbf{r}}_{ij;\mathbf{R}} \delta_{ij} \delta_\mathbf{0R}
+$$
+
+where $\bar{\mathbf{r}}_{ij;\mathbf{R}} =
+(\mathbf{r}_i + \mathbf{r}_j + \mathbf{R}) / 2$.
+The needed wannier centres are computed using
+the Marzari-Vanderbilt formula.
+
+This formalism is more accurate than the naive finite difference formula
+and converges faster. The formula is also manifestly translationally-invariant,
+that is, the results are the same if the system is translated by a whole.
+
+The default value is `false`.
+
 ## DOS
 
 Note that the behavior of the `dos` module is also influenced by the

src/comms.F90

@@ -154,6 +154,7 @@ module w90_comms
     module procedure comms_scatterv_real_2
     module procedure comms_scatterv_real_3
 !     module procedure comms_scatterv_cmplx
+    module procedure comms_scatterv_cmplx_3
     module procedure comms_scatterv_cmplx_4
   end interface comms_scatterv
 
@@ -218,6 +219,7 @@ module w90_comms
     module procedure comms_no_sync_scatterv_real_2
     module procedure comms_no_sync_scatterv_real_3
 !     module procedure comms_no_sync_scatterv_cmplx
+    module procedure comms_no_sync_scatterv_cmplx_3
     module procedure comms_no_sync_scatterv_cmplx_4
   end interface comms_no_sync_scatterv
 
@@ -1313,6 +1315,36 @@ subroutine comms_no_sync_scatterv_real_3(array, localcount, rootglobalarray, cou
 
   end subroutine comms_no_sync_scatterv_real_3
 
+  subroutine comms_no_sync_scatterv_cmplx_3(array, localcount, rootglobalarray, counts, displs, error, comm)
+    !! Scatter complex data from root node (array of rank 3)
+    implicit none
+
+    complex(kind=dp), intent(inout) :: array(:, :, :) !! local array for getting data
+    integer, intent(in) :: localcount !! localcount elements will be fetched from the root node
+    complex(kind=dp), intent(inout) :: rootglobalarray(:, :, :) !! array on the root node from which data will be sent
+    integer, intent(in) :: counts(0:) !! how data should be partitioned, see MPI documentation or function comms_array_split
+    integer, intent(in) :: displs(0:)
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
+
+#ifdef MPI
+    integer :: ierr
+
+    call mpi_scatterv(rootglobalarray, counts, displs, MPI_DOUBLE_COMPLEX, array, localcount, &
+                      MPI_DOUBLE_COMPLEX, root_id, comm%comm, ierr)
+
+    if (ierr .ne. MPI_SUCCESS) then
+      call set_base_error(error, 'Error in comms_scatterv_cmplx_3', code_mpi)
+      return
+    end if
+
+#else
+    !call zcopy(localcount, rootglobalarray, 1, array, 1)
+    array = rootglobalarray
+#endif
+
+  end subroutine comms_no_sync_scatterv_cmplx_3
+
   subroutine comms_no_sync_scatterv_cmplx_4(array, localcount, rootglobalarray, counts, displs, error, comm)
     !! Scatter complex data from root node (array of rank 4)
     implicit none
@@ -1860,6 +1892,25 @@ subroutine comms_scatterv_real_3(array, localcount, rootglobalarray, counts, dis
                                        error, comm)
   end subroutine comms_scatterv_real_3
 
+  subroutine comms_scatterv_cmplx_3(array, localcount, rootglobalarray, counts, displs, error, comm)
+    !! Scatter complex data from root node (array of rank 3)
+    implicit none
+
+    complex(kind=dp), intent(inout) :: array(:, :, :) !! local array for getting data
+    integer, intent(in) :: localcount !! localcount elements will be fetched from the root node
+    complex(kind=dp), intent(inout) :: rootglobalarray(:, :, :) !! array on the root node from which data will be sent
+    integer, intent(in) :: counts(0:) !! how data should be partitioned, see MPI documentation or function comms_array_split
+    integer, intent(in) :: displs(0:)
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
+
+    call comms_sync_error(comm, error, 0) ! sync error state across comm
+    if (allocated(error)) return
+
+    call comms_no_sync_scatterv_cmplx_3(array, localcount, rootglobalarray, counts, displs, &
+                                        error, comm)
+  end subroutine comms_scatterv_cmplx_3
+
   subroutine comms_scatterv_cmplx_4(array, localcount, rootglobalarray, counts, displs, error, comm)
     !! Scatter complex data from root node (array of rank 4)
     implicit none