Implementation of Stengel-Spaldin Functional in Wannierization

docs/docs/parameters/w90-wannierise-parameters.csv

@@ -21,3 +21,4 @@ slwf_num,I,"The number of objective WFs for selective localization"
 slwf_constrain,L,"Whether to constrain the centres of the objective WFs"
 slwf_lambda,R,"Value of the Lagrange multiplier for constraining the objective WFs"
 slwf_centres,P,"The centres to which the objective WFs are to be constrained"
+use_ss_functional, L,"If `true`, use Stengel-Spaldin functional and if `false`, use Marzari-Vanderbilt functional"

docs/docs/tutorials/tutorial_36.md

@@ -0,0 +1,113 @@
+# 34: 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)*
+
+- Input Files
+
+    - `silicon.scf` *The `pwscf` input file for ground
+        state calculation*
+
+    - `silicon.nscf` *The `pwscf` input file to obtain
+        Bloch states on a uniform grid*
+
+    - `silicon.pw2wan` *Input file for `pw2wannier90`*
+
+    - `silicon.win` *The `wannier90` input file*
+
+1. Run `pwscf` to obtain the ground state of silicon
+
+    ```bash title="Terminal"
+    pw.x < silicon.scf > scf.out
+    ```
+
+2. Run `pwscf` to obtain the Bloch states on a uniform
+    k-point grid. Note that we request the lower 4 (valence) bands
+
+    ```bash title="Terminal"
+    pw.x < silicon.nscf > nscf.out
+    ```
+
+3. Run `wannier90` to generate a list of the required overlaps (written
+    into the `silicon.nnkp` file).
+
+    ```bash title="Terminal"
+    wannier90.x -pp silicon
+    ```
+
+4. Run `pw2wannier90` to compute the overlap between Bloch states and
+    the projections for the starting guess (written in the `silicon.mmn`
+    and `silicon.amn` files).
+
+    ```bash title="Terminal"
+    pw2wannier90.x < silicon.pw2wan > pw2wan.out
+    ```
+
+5. Run `wannier90` to compute the MLWFs.
+
+    ```bash title="Terminal"
+    wannier90.x silicon
+    ```
+
+6. Plot the bandstructure by adding the following commands to the input
+    file `silicon.win`
+
+    ```vi title="Input file"
+    restart = plot
+
+    bands_plot = true
+    ```
+
+    and re-running `wannier90`. The files `silicon_band.dat` and
+    `silicon_band.gnu` are created. To plot the bandstructure using
+    gnuplot
+
+    ```bash title="Terminal"
+    gnuplot
+    ```
+
+    ```gnuplot title="Gnuplot shell"
+    load 'silicon_band.gnu'
+    ```
+
+    The k-point path for the bandstructure interpolation is set in the
+    `kpoint_path` block. Try plotting along different paths.
+
+## Stengel-Spaldin functional
+
+- Outline: *Obtain MLWFs for the valence and low-lying conduction-band
+    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>)*
+
+- Input Files
+
+    - `silicon.scf` *The `pwscf` input file for ground
+        state calculation*
+
+    - `silicon.nscf` *The `pwscf` input file to obtain
+        Bloch states on a uniform grid*
+
+    - `silicon.pw2wan` *Input file for `pw2wannier90`*
+
+    - `silicon.win` *The `wannier90` input file*
+
+- To run, the procedure is the same as the MV functional case.
+
+## Further ideas
+
+- Compare the spreads obtained by both functionals.
+  Increase the grid and find out that the two functionals converges to the same behavior.
+
+- Compare the Wannier-interpolated bandstructure with the both functionals.
+  Find out that the interpolated bandstructures are almost the same.

docs/docs/user_guide/wannier90/parameters.md

@@ -833,6 +833,19 @@ not possible if an explicit projection block is not defined).
 
 The default value is `false`.
 
+### `logical :: use_ss_functional`
+
+If `true`, use Stengel-Spaldin spread functional and if `false`,
+use Marzari-Vanderbilt spread functional.
+
+Both functionals converge to the same behavior if infinitely fine grid is used.
+Both of them are translationally invariant,
+but only Stengel-Spaldin functional is size consistent.
+
+For more information, refer to Phys. Rev. B 73, 075121.
+
+The default value is `false`.
+
 ### `integer :: num_guide_cycles`
 
 If `guiding_centres` is set to `true`, then the guiding centres are used

src/wannier90_readwrite.F90

@@ -713,6 +713,10 @@ subroutine w90_wannier90_readwrite_read_wannierise(settings, wann_control, num_w
                                    l_value=wann_control%guiding_centres%enable)
     if (allocated(error)) return
 
+    call w90_readwrite_get_keyword(settings, 'use_ss_functional', found, error, comm, &
+                                   l_value=wann_control%use_ss_functional)
+    if (allocated(error)) return
+
     call w90_readwrite_get_keyword(settings, 'num_guide_cycles', found, error, comm, &
                                    i_value=wann_control%guiding_centres%num_guide_cycles)
     if (allocated(error)) return

src/wannier90_types.F90

@@ -182,6 +182,7 @@ module w90_wannier90_types
     real(kind=dp) :: conv_tol = 1.0e-10_dp
     integer :: conv_window
     type(guiding_centres_type) :: guiding_centres
+    logical :: use_ss_functional = .false.
     real(kind=dp) :: fixed_step = -999.0_dp
     real(kind=dp) :: trial_step = 2.0_dp
     logical :: precond = .false.