DIY tutorial

As we learned from the basic tutorial, in WannierBerri calculations are done by `Calculator <https://wannier-berri.org/docs/calculators.html#wannierberri.calculators.Calculator>`__s. In this tutorial we will se how to create our own calculators.

As an example let’s consider anomalous Nernst conductivity (Xiao et al. 2006 ). In the zero-temperature limit \(\alpha_{\alpha\beta}^{\rm ANE}\) may be obtained from AHC \(\sigma_{\alpha\beta}(\epsilon)^{\rm AHE}\) evaluated over a dense grid of Fermi levels \(\epsilon\)

\[\alpha_{\alpha\beta}^{\rm ANE}=-\frac{1}{e}\int d\varepsilon \frac{\partial f}{\partial\varepsilon}\sigma_{\alpha\beta}(\varepsilon)\frac{\varepsilon-\mu}{T}, \label{eq:ane} \tag{1}\]

where \(f(\varepsilon)=1/\left(1+e^\frac{\varepsilon-\mu}{k_{\rm B}T}\right)\) and the AHC is defined as a Fermi-sea integral of Berry curvature

\begin{equation} \sigma^{\rm AHC}_{xy} = -\frac{e^2}{\hbar} \sum_n \int\frac{d\mathbf{k}}{(2\pi)^3} \Omega^n_\gamma f(\epsilon_n) \label{eq:ahc}\tag{2} \end{equation}

In the zero-temperature limit it reduced to

\[\alpha_{\alpha\beta}^{\rm ANE} \propto \frac{\partial \sigma_{\alpha\beta}(\varepsilon)}{\partial \varepsilon} \propto \int\frac{d\mathbf{k}}{(2\pi)^3} \Omega^n_\gamma f'(\epsilon_n) \tag{3}\label{eq:ane-ahc}\]

Where we omit the dimensional factor for simplicity.

Thus, there are two ways of calculating ANE :

  • via eq (\ref{eq:ane-ahc})

  • via eq (\ref{eq:ane})

let’s try eq (\ref{eq:ane-ahc}) first, because AHC is already implemented

Warning : this is an advanced tutorial. You something may not work from first try, you need to work with documentation to resolve the problem, but do not hesitate to ask the TA if you cannot.

[ ]:

# Preliminary

# Set environment variables - not mandatory but recommended if you use Parallel()
#import os
#os.environ['OPENBLAS_NUM_THREADS'] = '1'
#os.environ['MKL_NUM_THREADS'] = '1'


import wannierberri as wberri
import wannierberri.models
print (f"Using WannierBerri version {wberri.__version__}")
import numpy as np
import scipy
import matplotlib.pyplot as plt
from termcolor import cprint

#  This block is needed if you run this cell for a second time
#  because one cannot initiate two parallel environments at a time
try:
    parallel.shutdown()
except NameError:
    pass

# Chiose one of the options:

parallel = wberri.Parallel(num_cpus=2)
#parallel = wberri.Parallel()  # automatic detection
#parallel = wberri.Serial()


Using WannierBerri version 1.2.1
initializing ray with  {'num_cpus': 2}

2025-05-14 11:55:56,422 INFO worker.py:1879 -- Started a local Ray instance. View the dashboard at http://127.0.0.1:8265 
(raylet) [2025-05-14 15:04:56,665 E 18313 18313] (raylet) node_manager.cc:3287: 2 Workers (tasks / actors) killed due to memory pressure (OOM), 0 Workers crashed due to other reasons at node (ID: 41f5f8eb248e8ca703d3f988466034c2421af859bb67197d42528a9a, IP: 128.179.188.185) over the last time period. To see more information about the Workers killed on this node, use `ray logs raylet.out -ip 128.179.188.185`
(raylet)
(raylet) Refer to the documentation on how to address the out of memory issue: https://docs.ray.io/en/latest/ray-core/scheduling/ray-oom-prevention.html. Consider provisioning more memory on this node or reducing task parallelism by requesting more CPUs per task. To adjust the kill threshold, set the environment variable `RAY_memory_usage_threshold` when starting Ray. To disable worker killing, set the environment variable `RAY_memory_monitor_refresh_ms` to zero.

[2]:

# Instead of Wannier functions, we work with a Hlldane model here https://wannier-berri.org/docs/models.html

model=wberri.models.Haldane_ptb()
system = wberri.System_PythTB(model)

number of wannier functions: 2
R=0 found at position(s) [[3]]
shape of Ham_R = (7, 2, 2)
Real-space lattice:
 [[1.        0.        0.       ]
 [0.5       0.8660254 0.       ]
 [0.        0.        1.       ]]
Number of wannier functions: 2
Number of R points: 7
Recommended size of FFT grid [3 3 1]
Reading the system from PythTB finished successfully

anomalous Hall conductivity

Now, let’s evaluate the AHC on a grid of Ef-points, and then take a derivtive

[3]:

calculators = {}
Efermi = np.linspace(-3,3,21)
# Set a grid
grid = wberri.Grid(system, length=30 )   # length [ Ang] ~= 2*pi/dk

calculators ["ahc"] = wberri.calculators.static.AHC(Efermi=Efermi,
                                                    tetra=True,
                                                    constant_factor=1,
                                                    kwargs_formula={"external_terms":False})

result_run_ahc = wberri.run(system,
            grid=grid,
            calculators = calculators,
            parallel=parallel,
            adpt_num_iter=0,
            fout_name='Fe',
            restart=False,
            file_Klist="Klist_ahc.pickle",  # needed to restart a calculation in future
            print_Kpoints=False
            )



Minimal symmetric FFT grid :  [3 3 1]
Starting run()
Using the follwing calculators :
############################################################

 'ahc'  :  <wannierberri.calculators.static.AHC object at 0x74552017eb70>  : Anomalous Hall conductivity (:math:`s^3 \cdot A^2 / (kg \cdot m^3) = S/m`)

        | With Fermi sea integral Eq(11) in `Ref <https://www.nature.com/articles/s41524-021-00498-5>`__
        | Output: :math:`O = -e^2/\hbar \int [dk] \Omega f`
        | Instruction: :math:`j_\alpha = \sigma_{\alpha\beta} E_\beta = \epsilon_{\alpha\beta\delta} O_\delta E_\beta`
############################################################
Calculation on  grid - checking calculators for compatibility
ahc <wannierberri.calculators.static.AHC object at 0x74552017eb70>
All calculators are compatible
Grid is regular
The set of k points is a Grid() with NKdiv=[7 7 1], NKFFT=[5 5 1], NKtot=[35 35  1]
generating K_list
Done in 0.00014543533325195312 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.0006546974182128906 s
K_list contains 49 Irreducible points(100.0%) out of initial 7x7x1=49 grid
Done, sum of weights:1.0000000000000007
processing 49 K points : using  2 processes.
# K-points calculated  Wall time (sec)  Est. remaining (sec)
time for processing     49 K-points on   2 processes:     2.4343 ; per K-point          0.0497 ; proc-sec per K-point          0.0994
time1 =  0.0007510185241699219
Totally processed 49 K-points
run() finished

[4]:

!ls

Fe-ahc_iter-0000.dat  Fe-ane_iter-0000.npz   tutorial-wb-DIY-solution.ipynb
Fe-ahc_iter-0000.npz  Klist_ahc.pickle
Fe-ane_iter-0000.dat  tutorial-wb-DIY.ipynb

[5]:

# taking derivative and plotting

ef = result_run_ahc.results["ahc"].Energies[0]
ahc_z = result_run_ahc.results["ahc"].data[:,2]

# take the derivative
d_ef=ef[1]-ef[0]
d_ahc_z = (ahc_z[1:]-ahc_z[:-1])/d_ef
efnew = (ef[1:]+ef[:-1])/2

plt.plot(efnew, d_ahc_z)





[5]:

[<matplotlib.lines.Line2D at 0x7454fe78dee0>]
../../_images/tutorials_4_DIY_tutorial-wb-DIY-solution_6_1.png

Look into documentation

Look into documentation to see how the AHC calculator is defined : look for “Calculator” on https://docs.wannier-berri.org and press the “source” link.

One can see that it is based on the StaticCalculator, and only redefines the __init__ method. Namely, it redefines formula that stands under the integral (Berry curvature) and fder=0 meaning that the formula is weighted by the 0th derivative of the Fermi distribution.

copy the definition of AHC calculator and redefine the and below and define the init class

(hint : leave the factor the same, for simplicity) (hint : formula shouldbe imported from wannierberri.formula.covariant.DerOmega )

[6]:

from wannierberri.calculators.static import StaticCalculator
from wannierberri.formula.covariant import Omega
from wannierberri import factors


class ANC(StaticCalculator):
    r"""Anomalous Nernst conductivity

        | Output: :math:`O = \int [dk] \Omega f' `
        """

    def __init__(self,
                 constant_factor=1, # Figure out the pre-factor yourself
                 **kwargs):
        """describe input parameters here"""
        self.Formula = Omega
        self.fder = 1
        super().__init__(constant_factor=constant_factor, **kwargs)

Evaluate ANE using the new calculator and plot the result

[7]:

# insert the needed code below

calculators = {}
# Set a grid
grid = wberri.Grid(system, length=50 )   # length [ Ang] ~= 2*pi/dk

calculators ["ane"] = ANC(Efermi=Efermi,tetra=True)

result_run_ane = wberri.run(system,
            grid=grid,
            calculators = calculators,
            parallel=parallel,
            adpt_num_iter=0,
            fout_name='Fe',
            restart=False,
           # file_Klist="Klist_ahc.pickle"  # needed to restart a calculation in future
            )


Minimal symmetric FFT grid :  [3 3 1]
Starting run()
Using the follwing calculators :
############################################################

 'ane'  :  <__main__.ANC object at 0x7455200d9fa0>  : Anomalous Nernst conductivity

        | Output: :math:`O = \int [dk] \Omega f' `

############################################################
Calculation on  grid - checking calculators for compatibility
ane <__main__.ANC object at 0x7455200d9fa0>
All calculators are compatible
Grid is regular
The set of k points is a Grid() with NKdiv=[19 19  1], NKFFT=[3 3 1], NKtot=[57 57  1]
generating K_list
Done in 0.0009851455688476562 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.003620147705078125 s
K_list contains 361 Irreducible points(100.0%) out of initial 19x19x1=361 grid
Done, sum of weights:0.999999999999993
processing 361 K points : using  2 processes.
# K-points calculated  Wall time (sec)  Est. remaining (sec)

/home/stepan/github/wannier-berri/wannierberri/grid/grid.py:223: UserWarning:  the requested k-grid [58 58 50] was adjusted to [57 57 50].
  warnings.warn(f" the requested k-grid {NK} was adjusted to {NKFFT * NKdiv}. ")

time for processing    361 K-points on   2 processes:     0.6829 ; per K-point          0.0019 ; proc-sec per K-point          0.0038
time1 =  0.0052797794342041016
Totally processed 361 K-points
run() finished

[8]:


ef = result_run_ane.results["ane"].Energies[0]
ane_z = result_run_ane.results["ane"].dataSmooth[:,2]


plt.scatter(efnew,d_ahc_z)
plt.plot(ef, ane_z)
plt.show()
../../_images/tutorials_4_DIY_tutorial-wb-DIY-solution_11_0.png

Questions:

  1. Make the k-points or Efermi denser. Will the agreement improve?

  2. we calculated the AHC and ANE with zero refinement iterations, and the results matched well (or at least they should). If we calculate them separately, but with some refinement iterations, will the results match?

  3. What if we use refinement itrations and run the two calculations together, in one run() call ?

FormulaProduct

Look at the definition of erryDipole_FermiSurf

https://wannier-berri.org/_modules/wannierberri/calculators/static.html

It evaluates

\begin{equation} D_{ab} =\sum_n \int\frac{d\mathbf{k}}{(2\pi)^3} v_a\Omega^n_b f(\epsilon_n) \label{eq:ahc}\tag{2} \end{equation}

also look at definition of the Feormula https://wannier-berri.org/docs/formula.html

you can see that it is base on a FormulaProduct. Essentially, by analogy you may define any tensor, e.g.

\[\sum_n \int\frac{d\mathbf{k}}{(2\pi)^3} v_a\Omega^n_b v_c f''(\epsilon_n)\]

i.e. product of several quantities weighted by rge second derivative of Fermi distribution. You may write any analogous formula, if you want. Try to define corresponding FormulaProduct and StaticCalculator classes, and try to evaluate them.

[ ]:



[ ]:

Problem :

try to prove (analytically and numerically) that

\[\sum_n \int\frac{d\mathbf{k}}{(2\pi)^3} \partial_a \partial_b \partial_c \epsilon_{n\mathbf{k}} f(\epsilon_{n\mathbf{k}}) \propto \sum_n \int\frac{d\mathbf{k}}{(2\pi)^3} \partial_a \epsilon_{n\mathbf{k}} \partial_b \epsilon_{n\mathbf{k}} \partial_c \epsilon_{n\mathbf{k}} f''(\epsilon_{n\mathbf{k}})\]

and tell me what is the factor missing?

** Have fun ! **

[ ]: