For formulation, see J, Kunes et al., Comp. Phys. Commun. 181, 1888 (2010)
Userguide for wien2wannier, https://wien2wannier.github.io
Userguide for wannier90, http://www.wannier.org

Step.1: prepare DFT data and input files

1. Prepare converged (non spin-polarized) DFT result with Wien2k package

2. Make “case.fermi” file with the Fermi energy (at the first line of the file) in the covered DFT result
grep :FER case.scf” finds Fermi energy (Ry unit) in the scf file

cat case.fermi
0.4782867816

Note for SOC. For the calculation with the spin-orbit coupling (SOC), copy case.fermi to case.fermiup, case.fermidn

cp case.fermi case.fermiup
cp case.fermi case.fermidn

3. Set k-grid in the whole (nonsymmetrized) BZ to be used in the wannierization (wannier90)
x kgen -fbz” with NO SHIFT generates the k-grid for the wannierization. Since the devision (ndiv for kx, ky, kz, see below) defines the maximum distance (nx, ny, nz in hr.dat file) of hopping integrals (matrix elements) computable in real space, the density of k-grid could affect the quality of wannierization.

x kgen -fbz
1  symmetry operations without inversion
NUMBER OF K-POINTS IN WHOLE CELL: (0 allows to specify 3 divisions of G)
500   ← total number of k points 
length of reciprocal lattice vectors:   1.379   1.379   1.379   7.937   7.937   7.937
Shift of k-mesh allowed. Do you want to shift: (0=no, 1=shift)
0     ← no shift 
343  k-points generated, ndiv=           7           7           7
STOP KGEN ENDS

4. Run “x lapw1″
Compute eigen energies and vectors on the new k-mesh created in Step.1
Thus one must rerun x lapw1 when the k-grid (x kgen -fbz) is changed.

## Important note when user-defined local matrix is used below ##
x lapw1 must be executed with “case.struct” with the original local matrix used in the self-consistent DFT calculation. Thus, when one has to rerun x lapw1 e.g.for a new k-grid, one must copy case.struct_orig, see blow, as case.struct!

Note for SOC. For the calculation with the spin-orbit coupling, to mimic the spin-polarized implementation, one needs to copy vector files and run x lapw1 for up spin and down spin, followed by x lapwso -up

cp case.vns case.vnsup
cp case.vns case.vnsdn
cp case.vsp case.vspup
cp case.vsp case.vspdn
x lapw1 -up
x lapw1 -dn
x lapwso -up

(Optional)
5. Modify local matrix in case.struct file
If the local matrix needs to be arranged to realize Wannier functions in different axis (i.e. define initial projections (in case.inwf) in different quantization axis) from the one used in Wien2k, one can provide a user-designed local matrix in case.struct file. Since the original case.struct file will be used later (x lapw1 etc), it must be saved. The following file names (orig, qtl) are used below.

cp case.struct case.struct_orig
cp case.struct case.struct_qtl   (←here user inserts his/her designed local matrix)
Example 
case.struct_orig
Ir         NPT=  781  R0=0.00000500 RMT=    1.9700   Z: 77.000  
LOCAL ROT MATRIX:    1.0000000 0.0000000 0.0000000
                     0.0000000 1.0000000 0.0000000
                     0.0000000 0.0000000 1.0000000

case.struct_qtl
Ir         NPT=  781  R0=0.00000500 RMT=    1.9700   Z: 77.000  
LOCAL ROT MATRIX:   -0.5665170-0.8240500 0.0000000
                     0.8240500-0.5665170 0.0000000
                     0.0000000 0.0000000 1.0000000

Step.2: prepare initial projection “case.inwf

1. Find indices of the bands to be included in Wannierization
The indices of the bands can be found in case.output file or alternatively by using “x findbands -emin ** -emax** -efermi **” (note. eV for emin and emax, Ry unit for efermi) command that will find and write the indices in case.outputfind file. With SOC, an argument “-so” is needed.

Example of case.outputfind file
# findbands: emin=  -7.000 eV emax=   1.500 eV efermi=  0.47729 Ry
# k-point first last #bands
...
Bloch bands in the interval
at all k : 6 12 7 (from index1 to index2, total no. of bands included in the given energy window at all k points)
at any k : 5 13 8 (from index1 to index2, total no. of bands included in the given energy window at any k points)

2. Define initial projections in case.inwf file
Write case.inwf file which defines initial projections (on atom, angular profile etc.) in the wanneirization.

A few tips (for beginners)
– atom index (#iat) is defined in case.struct file.
#LJMAX defines no. of the expansion order used to approximate exp(-ikb), usually 3~4.
– #WF is the number of the Wannier functions to be computed
– #min and #max band reads the index of min and max bands in the target window in the wannierization (#max-#min+1 ≧ #WF)

Example of case.inwf file (one transition metal (3d) and one oxygen (2p) in wannierization)
BOTH            #AMN, #MMN, or #BOTH
 5  12          #min band, #max band
 3   8          #LJMAX in exp(ibr) expansion, #WF
 2         # 1:d_x2-y2 (eg) = 1/{sqrt(2)} |l=2, m=2> + 1/{sqrt(2)} |l=2, m=-2>
 2 2 -2   0.70710678   0.00000000   # iat, #l, #m, #Re(coeff), #Im(coeff)
 2 2  2   0.70710678   0.00000000   → projection onto "atom=2, l=2, m=2 with (Re,Im)=(1/sqrt(2), 0)"
 1         # 2:d_z2 (eg) = |l=2, m=0>
 2 2  0   1.00000000   0.00000000
 2         # 3:d_xy (t2g) = -i/{sqrt(2)} |l=2, m=2> + i/{sqrt(2)} |l=2, m=-2> 
 2 2 -2   0.00000000   0.70710678
 2 2  2   0.00000000  -0.70710678
 2         # 4:d_yz (t2g) = i/{sqrt(2)} |l=2, m=1> + i/{sqrt(2)} |l=2, m=-1>
 2 2 -1   0.00000000   0.70710678
 2 2  1   0.00000000   0.70710678
 2         # 5:d_xz (t2g) = -1/{sqrt(2)} |l=2, m=1> + 1/{sqrt(2)} |l=2, m=-1>
 2 2 -1   0.70710678   0.00000000
 2 2  1  -0.70710678   0.00000000
 2         # 6:px
 1 1 -1   0.70710678   0.00000000
 1 1  1  -0.70710678   0.00000000
 2         # 7:py
 1 1 -1   0.00000000   0.70710678
 1 1  1   0.00000000   0.70710678
 1         # 8:pz
 1 1  0   1.00000000   0.00000000

When the spin-orbit coupling (SOC) is included in the wannierization (SCF with SOC must be done before exec)

Three inwf files are required: case.inwf, case.inwfup, case.inwfdn
1. case.inwf include full projections (i.e. both up- and down-spin channels)
2. case.inwfup defines projection onto the up-spin channel
3. case.inwfdn defines projection onto the down-spin channel

Below gives an example for d-only model (single transition metal in the unit cell), and compute wannier functions in cubic harmonics symmetry (x2-y2, z2, xy, yz, zx) with spin up or down. As seen in this example, case.inwfup and case.inwfdn contain the projection on up and down spinor (full projection is given in case.inwf). In this simplest example, each target wannier function is composed of one spinor (i.e. specific up or down). In practical purpose with SOC, one may want to target the wannier function with mixed spinor, for example, jj-basis (for Jeff basis in t2g materials or 4f rare-earth compounds). In such cases the target wannier function will have the (initial) projection onto both the up and down channel.

Example of case.inwf file with up (1~5) and down (6~10) spin  
...
 3   10          #LJMAX in exp(ibr) expansion, #WF
 2         # 1:d_x2-y2 (eg) / up-spin
 2 2 -2   0.70710678   0.00000000   
 2 2  2   0.70710678   0.00000000 
 1         # 2:d_z2 (eg) / up-spin
 1 2  0   1.00000000   0.00000000
 2         # 3:d_xy (t2g) / up-spin
 1 2 -2   0.00000000   0.70710678
 1 2  2   0.00000000  -0.70710678
 2         # 4:d_yz (t2g) / up-spin 
 1 2 -1   0.00000000   0.70710678
 1 2  1   0.00000000   0.70710678
 2         # 5:d_xz (t2g) / up-spin 
 1 2 -1   0.70710678   0.00000000
 1 2  1  -0.70710678   0.00000000
 2         # 6:d_x2-y2 (eg) / down-spin
 2 2 -2   0.70710678   0.00000000 
 2 2  2   0.70710678   0.00000000
1         # 7:d_z2 (eg) / down-spin
 1 2  0   1.00000000   0.00000000
 2         # 8:d_xy (t2g) / down-spin
 1 2 -2   0.00000000   0.70710678
 1 2  2   0.00000000  -0.70710678
 2         # 9:d_yz (t2g) / down-spin
 1 2 -1   0.00000000   0.70710678
 1 2  1   0.00000000   0.70710678
 2         # 10:d_xz (t2g) / down-spin
 1 2 -1   0.70710678   0.00000000
 1 2  1  -0.70710678   0.00000000
Corresponding case.inwfup file 
...
 3   10          #LJMAX in exp(ibr) expansion, #WF
 2         # 1:d_x2-y2 (eg) / up-spin
 2 2 -2   0.70710678   0.00000000   
 2 2  2   0.70710678   0.00000000 
 1         # 2:d_z2 (eg) / up-spin
 1 2  0   1.00000000   0.00000000
 2         # 3:d_xy (t2g) / up-spin
 1 2 -2   0.00000000   0.70710678
 1 2  2   0.00000000  -0.70710678
 2         # 4:d_yz (t2g) / up-spin 
 1 2 -1   0.00000000   0.70710678
 1 2  1   0.00000000   0.70710678
 2         # 5:d_xz (t2g) / up-spin 
 1 2 -1   0.70710678   0.00000000
 1 2  1  -0.70710678   0.00000000
0no projection on up-spin component
0   no projection on up-spin component
0   no projection on up-spin component
0   no projection on up-spin component
0   no projection on up-spin component
Corresponding case.inwfdn file  
...
 3   10          #LJMAX in exp(ibr) expansion, #WF
0   no projection on down-spin component
0   no projection on down-spin component
0no projection on down-spin component
0no projection on down-spin component
0no projection on down-spin component
 2         # 6:d_x2-y2 (eg) / down-spin
 2 2 -2   0.70710678   0.00000000 
 2 2  2   0.70710678   0.00000000
1         # 7:d_z2 (eg) / down-spin
 1 2  0   1.00000000   0.00000000
 2         # 8:d_xy (t2g) / down-spin
 1 2 -2   0.00000000   0.70710678
 1 2  2   0.00000000  -0.70710678
 2         # 9:d_yz (t2g) / down-spin
 1 2 -1   0.00000000   0.70710678
 1 2  1   0.00000000   0.70710678
 2         # 10:d_xz (t2g) / down-spin
 1 2 -1   0.70710678   0.00000000
 1 2  1  -0.70710678   0.00000000

Tics. opw_wannierbasis in opw2x package gives expansion coefficients in jj-basis. Note that the Clebsh-Goldan coefficients do not include the so-called Shortley-Condon phase. In the 4th column of the execution, see below, +1 and -1 mean up- and down-spin component, respectively.

opw_wannierbasis
enter L(ang)  : 2
 [1] C.F basis / [2] J-J basis 2
 [1] for wannier / [2]  for opw2x_spec 1  
 
  1 (J=5/2, Jz=5/2)
*at  2  2  1   1.00000000   0.00000000
  2 (J=5/2, Jz=3/2)
*at  2  1  1  -0.89442719   0.00000000
*at  2  2 -1  -0.44721360   0.00000000
  2 (J=5/2, Jz=1/2)
*at  2  0  1   0.77459667   0.00000000
*at  2  1 -1   0.63245553   0.00000000
  2 (J=5/2, Jz=-1/2)
*at  2 -1  1  -0.63245553   0.00000000
*at  2  0 -1  -0.77459667   0.00000000
  2 (J=5/2, Jz=-3/2)
*at  2 -2  1   0.44721360   0.00000000
*at  2 -1 -1   0.89442719   0.00000000
  1 (J=5/2, Jz=-5/2)
*at  2 -2 -1  -1.00000000   0.00000000
  2 (J=3/2, Jz=3/2)
*at  2  1  1   0.44721360   0.00000000
*at  2  2 -1  -0.89442719   0.00000000
  2 (J=3/2, Jz=1/2)
*at  2  0  1  -0.63245553   0.00000000
*at  2  1 -1   0.77459667   0.00000000
  2 (J=3/2, Jz=-1/2)
*at  2 -1  1   0.77459667   0.00000000
*at  2  0 -1  -0.63245553   0.00000000
  2 (J=3/2, Jz=-3/2) 
*at  2 -2  1  -0.89442719   0.00000000
*at  2 -1 -1   0.44721360   0.00000000 

Step.3: wannierization using wannier90 with the aid of wien2wannier

1. run “write_win”
This command makes “case.win” file using information written in case.inwf file created above. The case.win file is input for wannier90
Note for SOC. For wannierization with the spin-orbit coupling (SOC), execute with -up flag as “write_win -up”

A few tips.
– When num_bands (calculated with #min #max bands in case.inwf) > num_wann (#WF in case.inwf), the disentanglement option (dis_froz_min, dis_froz_max, dis_mix_ratio) can be activated. Some examples showing how the disentanglement option works (will) are provided later.
– In order to enforce the symmetry of the initial projections (given in case.inwf) in the resultant wannier functions, one can set num_itr=0, which skips the minimization (localization) process in constructing the wannier functions.
– In the latest version of wannier90, “write_hr” is used instead of “hr_plot” (old version).

case.win file
 num_bands       = 8
 num_wann        = 8               

!!! Disentanglement parameters !!!  option for disentanglement (when num_bands > num_wann case) 
!dis_froz_min     = 7.              dis_froz_min : bottom of the frozen energy window
!dis_froz_max     = 9.              dis_froz_max : top of the frozen energy window
!dis_mix_ratio    = 0.5             dis_mix_ratio : mixing ratio during minimization of Ω

!!! Iterations &c. !!!
 iprint                = 1
 num_iter              = 10000      max no.of iterations in minimization 
 num_print_cycles      =   100
 conv_window           = 3
!conv_tol              = 1e-10
 dis_num_iter          = 10000
!dis_conv_window       = 3
!dis_conv_tol          = 1e-10
!restart               = default ! wannierise ! plot ! transport

!!! Post-processing options !!!
 write_proj             = .true.
 write_xyz              = .true.
 translate_home_cell    = .true.
 write_hr               = .true.
!fermi_surface_plot     = .true.
!dos                    = .false.
!dos_project            =
!dos_kmesh              = 50 50 1
!dos_energy_min         = -1.8
!dos_energy_max         = +1.5

!!! Band structure !!!
!!! needs `kpoint_path' block
!bands_plot             = .true.
 bands_num_points       = 50
!bands_plot_format      = gnuplot xmgrace
!bands_plot_project     = 1-3
bands_plot_mode        = s-k | cut [Slater-Koster | truncate Hamiltonian] use s-k mode for standard usage 

2. run “x wannier90 -pp”
This command prepares a list of nn k-points in “case.nnkp” based on the “case.win” file created just above.

3. run “x w2w”
This command compute overlaps Mmn and initial projections Amn in “case.mmn” and “case.amn“, respectively.

Note for SOC. x w2w must be executed for both up- and down-spin channels with -so flag.

x w2w -so -up
x w2w -so -dn

4. run “x wannier90″
This runs wannierization using wannier90. With SOC, an argument with “-so” is needed.
Compare the band structure calculated with the Wannier-interporlated H(k), written in “case_band.dat”, and the one calculated by Wien2K.
GNUPLOT COMAND : p ‘case.spaghetti_ene’ u ($4/0.53):5, ‘case_band.dat’ w l
(Here conversion from Bohr-1 to Angstrom-1 unit, 1Bohr = 0.529177Angstrom)

Templates for various d, f materials