# Examples

**Note**: The following exampole is written with an outdated (but still working so far) inteface using ithe integrate() and tabulate() functions.
This is done only for demonstaration purposes. Updated tutorials in Jupyter notebooks are available here

## AHC of bcc Fe

First let’s consider a simple example - we will calculate a anomalous Hall conductivity of bcc iron, and also visualize the Berry curvature over the Fermi surface. This is equivalent to example18 of Wannier90.

We assume that we have performed the ab-initio calculations and constructed the MLWFs for bcc iron with magnetization along [001] direction (\(z\)). For information how to do it, please refer to the manuals of the abinitio codes and Wannier90. Also a nicely ilustrated tutorial is published on the WanierToools website.

After that, the calculation is performed by the following short python script. First, we import the needed packages

```
import wannierberri as wberri
import numpy as np
```

Then we read the information about the system. Either we read the information abot Wanier functions:

```
system=wberri.System_w90('Fe',berry=True)
```

thie reads the files `’Fe.chk’`

, `’Fe.eig’`

, `’Fe.mmn’`

. [6]
Alternatively, we can read all information from a file `Fe_tb.dat`

,
which is also written by Wannier90, or maybe composed by user from any
tight-binding model. This can be done by the following line:

```
system=wberri.System_tb('Fe_tb.dat',berry=True)
```

Next, we define the symmetries of the system that we wish to take into account:

```
generators = ['Inversion','C4z','TimeReversal*C2x']
system.set_symmetry(generators)
```

In the ab initio calculation we have specified the magnetization along
the \(z\) axis, hence the symmetries that are preserved are
inversion \(\cal I\), 4-fold rotation around the \(z\) axis
\(C_{4z}\),and a combination of time-reversal \(\cal T\) and
2-fold rotation around the \(x\) axis \(C_{2x}\). Here we need
only the generators of the symmetry group. The other symmetries will be
automatically obtained by taking products of these generators. E.g. the
mirror \(M_z=(C_{4z})^2\cdot \cal I\). Next we need to set the
\({\bf k}\)-point grid. Most conveniently it can be done by setting
the `’length’`

parameter.

```
grid=wberri.Grid(system,length=100)
```

This will guarantee the grid to be consistent with the symmetries, and the spacing of \(k\)-points will be \(\Delta k\approx \frac{2\pi}{\rm length}\). In this particular case a grid of \(52\times52\times52\) will be generated, but that is dependent on the size of the unit cell.

To determine all parameters of parallel execution we define an object
of class `Parallel`

```
parallel=wberri.Parallel(method="ray",num_cpus=num_proc)
```

Or, if for a serial execution just set paralle=None.

Next, we want to integrate the Berry curvature to get the AHC. This is
done by the `wberri.integrate`

method.

```
wberri.integrate(system, grid,
Efermi=np.linspace(12.,13.,1001),
smearEf=10, # 10K
quantities=["ahc","dos","cumdos"],
parallel = parallel,
adpt_num_iter=10,
fout_name="Fe")
```

and in addition to AHC we evaluate the density of states (DOS) and
cumulative DOS. We consider Fermi level as a free parameter, and we scan
over a set of Fermi levels from 12 to 13 eV with a step of 1 meV. To
avoid too strong jittering of the curve, we introduce a small smearing
over the Fermi level corresponding to temperature 10K (\(\sim1\)
meV). It is known, that in BZ integrations, some \({\bf k}\) points
give huge spites in the \(E_F\)-resolved curves. This is especially
strong for Berry curvature, which diverges near band degeneracies and
avoided crossings, that fall close to the Fermi level. To make the
calculation more precise around those points, and adaptive recursive
refinement algorithm is used, and we set the number of iterations to
10 [7]. The integration is done in parallel by means of the
`multiprocessing`

module and the parameter `’numproc’`

specifies
that a `Pool`

of 16 worker processes is used. Both the smeared and
unsmeared result is written to the files, e.g. `Fe-???_iter-????.dat`

.
In particular, from the cumulative dos (`Fe-cumdos_iter-????.dat`

) we
can find the precise position of the Fermi level \(E_F=12.610\) eV —
the energy at which the cumulative DOS reaches 8 electrons per unit
cell. This is much more accurate, then the result of the `postw90.x`

which
is evaluated from a coarse abinitio grid.

Next, it is instructive to plot the AHC after each iteration (Fig. 2). One can see that already after a few iterations most of the chaotic peaks are removed, and we can get a reasonably sooth curve already starting from a not very dense grid.

## Tabulating Berry curvature

Now we wish to visualize the Berry curvature to see, from which parts of the BZ mostly contribute for the AHC. For that purpose we employ the following method:

```
wberri.tabulate(system, grid,
quantities=["berry"],
frmsf_name="Fe",
parallel = parallel,
ibands=np.arange(4,10),
Ef0=12.6)
```

Which produce files `Fe_berry-x.frmsf`

, `Fe_berry-y.frmsf`

,
`Fe_berry-z.frmsf`

, containing the Energies and Berry curvature of
bands `4,5,6,7,8,9`

[8] over. The format is chosen such that the
files can be directly passed to the `FermiSurfer`

visualization
tool(“Fermisurfer Visualization Tool,” n.d.; Kawamura 2019) [9]
However, the hotspots of Berry curvature usually present tiny areas with
huge magnitude of \({\cal O}\). Hence to get a smoother picture we
apply a logarithmic scale as

with \(x_0=1\) ???. Now we can use the `FermiSurfer`

to produce
Fig. 1

It is often more convinient to get `TABresult`

object, and saving the large formatted files only when needed, e.g.

```
tab_result = wberri.tabulate(system,
grid,
quantities=["berry"],
)
```

Ths object can be pickled (saved to disk):

```
pickle.dump(tab_result,open("Te_berry.pickle","wb")) #Then the result is saved in the file "Te_berry.pickle" now.
```

and then loaded in another script:

```
tab_result=pickle.load(open("Te_berry.pickle","rb"))
```

and then get the FermiSurfer files only for specific components and bands, e.g.:

```
open("Te-berry-x.frmsf","w").write(tab_result.fermiSurfer(quantity='berry',component='x',efermi=Ef0,npar=num_proc,iband=np.arange(14,18))) #iband - counts from zero
open("Te-E.frmsf","w").write(tab_result.fermiSurfer(efermi=Ef0,npar=num_proc))
```

One can also extract the data as arrays like this:

```
Energy_CB=tab_result.get_data(iband=18,quantity='Energy')
berry=tab_result.get_data(iband=18,quantity='berry',component='z')
```

to further process in any other way.

## Plotting lines and planes from the 3D grid(tab_plot)

Once we have the `TABresult`

object, we can plot the band structure with a quantity on a path or a plane of k-points

Path:

```
python3 -m wannierberri.utils.tab_plot tab_result.pickle type=Line quantity=True kpath=0,0,0,0,0,60 namelist=G,Z qtype=berry component=z
```

plane:

```
python3 -m wannierberri.utils.tab_plot tab_result.pickle type=Plane quantity=True Efermi=-0.5 vec1=1,0,0 vec2=0,1,0 qtype=berry component=z
```

For more comments on the parameters please check

```
python3 -m wannierberri.utils.tab_plot -h
```

## Plotting along lines

Since version 0.8.2 it is also possible to tabulate bands and other properties (e.g. spin or berry curvature) along any arbitrary pathes in the BZ. This feature is not well tested yet, so please try it and report any problems you meet. An example for valence band of tellurium along the K-H-K line is given below.

```
path=wberri.Path(system,
k_nodes=[[1./3,1./3,0],[1./3,1./3,0.5],None,[1./3,1./3,0.5],[1./3,1./3,1],],
labels=["K1","H","H","K2"],
length=1000)
path_result=wberri.tabulate(system,
grid=path,
quantities=['berry'],
)
```

now plot the result

```
path_result.plot_path_fat( path,
quantity='berry',
component='z',
save_file="Te-berry-VB.pdf",
Eshift=0,
Emin=5.65, Emax=5.85,
iband=None,
mode="fatband",
fatfactor=20,
cut_k=True
)
```

Which should produce a figure like this:

or you may get the data to plot in whatever way you like

```
k=path.getKline()
E=path_result.get_data(quantity='Energy',iband=(16,17,18,19))
curv=path_result.get_data(quantity='berry',iband=(16,17,18,19),component="z")
```

## Optical conductivity

The Kubo formula (2) for the (interband) optical conductivity may be evaluated by adding ‘opt_conductivity’ to the list of quantities to integrate, for more details see the example below.

The implementation is based on the one in postw90. Thus, with identical input and equivalent parameters, it reproduces the results from postw90. Note, however, that when using the full power of wannier-berri (symmetries, adaptive refinement etc.) small deviations are to be expected. Please refer to sec-benchmark.

In contrast to the other quantities currently implemented, this is an optical quantity and therefore requires a list of frequencies rather than a list of Fermi energies; the relevant argument name for wberri.integrate() is omega and its values are expected to be in eV. Additionally, there are several options (following the usual syntax of wannier-berri) that can be specified:

‘mu’: chemical potential in units of eV

‘kBT’: temperature in units of eV/kB (can also be 0)

‘smr_fixed_width’: fixed smearing parameter in units of eV

‘smr_type’: analytical form of the broadened delta function (must be one of ‘Lorentzian’ or ‘Gaussian’)

Here smearing refers to the approximation of the delta function in the Kubo formula.

An example call might look as follows (with appropriate initialization):

```
wberri.integrate(system,
grid = grid,
omega = np.linspace(0., 7., 701),
Efermi=np.linspace(12.,13.,21),
smearEf = 100,
quantities = [ 'opt_conductivity' ],
adpt_num_iter = 10,
fout_name = 'Fe',
restart = False,
parameters = { 'smr_fixed_width': 0.01, 'smr_type':'Gaussian' }
)
```

However, note that smearEf parameter has a sence only when many Fermi levels are considered, (when difference between Fermi levels is smaller then smearEf in eV) otherwise ` ‘kBT’:<value in eV> ` shoould be added to the parameters dictionary. But do not use both kBT and smearEf at the same time.

## Spin Hall conductivity

Utilizing the similar formula as the optical conductivity, the intrinsic spin Hall conductivity(SHC) can be calculated. Two methods are available:
`SHCryoo`

and `SHCqiao`

. The former requires `.sHu`

and `.sIu`

files from pw2wannier90.x(see mmn2uHu),
while the latter does not and instead uses an approximation \(\mathbf{1}=\sum_{l\in \it{ab\,initio}} \vert u_{l{\bf q}}\rangle\langle u_{l{\bf q}}\vert\).

There are more additional quantities than optical conductivity that can be specified. Note that if one of these is not specified, the module will calculate all the 27 components of SHC.

‘shc_alpha’: direction of spin current (1, 2, 3)

‘shc_beta’: direction of applied electric field (1, 2, 3)

‘shc_gamma’: direction of spin polarization (1, 2, 3)

You can check which component to calculate by the code below, considering the point group your material belongs to.

```
wberri.symmetry.Group(['Inversion', 'C4x', 'C4y', 'C4z']).get_symmetric_components(3, False, False)
```

You can also set symmetry before integrating, but note that specification of `shc_alpha`

, `shc_beta`

, and `shc_gamma`

may not work with symmetry when the point group does not belong to the Laue group cubic I, so it is not recommended to use specification and symmetry at the same time.

```
# with symmetry
system=wberri.System_w90(seedname='pt', SHCryoo=True, SHCqiao=True, use_ws=True, transl_inv=False)
generators=[SYM.Inversion, SYM.C4z, SYM.C4x, SYM.C4y]
system.set_symmetry(generators)
wberri.integrate(system,
grid = grid,
omega = np.linspace(0., 7., 701),
Efermi=np.linspace(12.,13.,21),
smearEf = 100,
quantities = ['SHC_ryoo', 'SHC_qiao'],
adpt_num_iter = 10,
fout_name = 'pt',
restart = False,
parameters = { 'smr_fixed_width': 0.01, 'smr_type':'Gaussian'}
)
```

```
# with specification of shc_alpha, shc_beta, and shc_gamma
system=wberri.System_w90(seedname='pt', SHCryoo=True, SHCqiao=True, use_ws=True, transl_inv=False)
wberri.integrate(system,
grid = grid,
omega = np.linspace(0., 7., 701),
Efermi=np.linspace(12.,13.,21),
smearEf = 100,
quantities = ['SHC_ryoo', 'SHC_qiao'],
adpt_num_iter = 10,
fout_name = 'pt',
restart = False,
parameters = { 'smr_fixed_width': 0.01, 'smr_type':'Gaussian', 'shc_alpha':1, 'shc_beta':2, 'shc_gamma':3 }
)
```