Brillouin zone
The Brillouin zone objects are all special classes enabling easy manipulation of an underlying physical quantity.
Quite often a physical quantity will be required to be averaged, or calculated individually over a number of k-points. In this regard the Brillouin zone objects can help.
The BrillouinZone object allows direct looping of contained k-points while invoking particular methods from the contained object. This is best shown with an example:
>>> import sisl as si
>>> H = si.Hamiltonian(...)
>>> bz = si.BrillouinZone(H)
>>> bz.apply.array.eigh()
This will calculate eigenvalues for all k-points associated with the BrillouinZone and
return everything as an array. The dispatch property of
the BrillouinZone object has several use cases (here array is shown).
This may be extremely convenient when calculating band-structures:
>>> H = si.Hamiltonian(...)
>>> bs = si.BandStructure(H, [[0, 0, 0], [0.5, 0, 0]], 100)
>>> bs_eig = bs.apply.array.eigh()
>>> plt.plot(bs.lineark(), bs_eig)
and then you have all eigenvalues for all the k-points along the path.
Multiple quantities
Sometimes one may want to post-process the data for each k-point. As an example lets post-process the DOS on a per k-point basis while calculating the average:
>>> H = si.Hamiltonian(...)
>>> mp = si.MonkhorstPack(H, [10, 10, 10])
>>> E = np.linspace(-2, 2, 100)
>>> def wrap_DOS(eigenstate):
... # Calculate the DOS for the eigenstates
... DOS = eigenstate.DOS(E)
... # Calculate the velocity for the eigenstates
... v = eigenstate.velocity()
... V = (v ** 2).sum(1)
... return DOS.reshape(-1, 1) * v ** 2 / V.reshape(-1, 1)
>>> DOS = mp.apply.average.eigenstate(wrap=wrap_DOS, eta=True)
This will, calculate the Monkhorst pack k-averaged DOS split into 3 Cartesian
directions based on the eigenstates velocity direction. This method of manipulating
the result can be extremely powerful to calculate many quantities while running an
efficient BrillouinZone average. The eta flag will print, to stdout, a progress-bar.
The usage of the wrap method are also passed optional arguments, parent which is
H in the above example. k and weight are the current k-point and weight of the
corresponding k-point. An example could be to manipulate the DOS depending on the k-point and
weight:
>>> H = si.Hamiltonian(...)
>>> mp = si.MonkhorstPack(H, [10, 10, 10])
>>> E = np.linspace(-2, 2, 100)
>>> def wrap_DOS(eigenstate, k, weight):
... # Calculate the DOS for the eigenstates and weight by k_x and weight
... return eigenstate.DOS(E) * k[0] * weight
>>> DOS = mp.apply.sum.eigenstate(wrap=wrap_DOS, eta=True)
When using wrap to calculate more than one quantity per eigenstate it may be advantageous
to use oplist to handle cases of BrillouinZone.apply.average and BrillouinZone.apply.sum.
>>> H = si.Hamiltonian(...)
>>> mp = si.MonkhorstPack(H, [10, 10, 10])
>>> E = np.linspace(-2, 2, 100)
>>> def wrap_multiple(eigenstate):
... # Calculate DOS/PDOS for eigenstates
... DOS = eigenstate.DOS(E)
... PDOS = eigenstate.PDOS(E)
... # Calculate velocity for the eigenstates
... v = eigenstate.velocity()
... return si.oplist([DOS, PDOS, v])
>>> DOS, PDOS, v = mp.apply.average.eigenstate(wrap=wrap_multiple, eta=True)
Which does mathematical operations (averaging/summing) using oplist.
In some cases quantities are needed for all \(k\) points and in such cases
it may not always be that the returned quantities are commensurate.
Lets re-use the previous wrap_multiple function and try and return the
full quantity:
>>> DOS_PDOS_v = mp.apply.eigenstate(wrap=wrap_multiple, eta=True)
This will raise an error since wrap_multiple returns an oplist (same as a list)
and thus is unable to convert this into an equivalent numpy.ndarray. Additionally
this can not be merged together in a single numpy.ndarray since the shapes of the returned
quantities are not commensurate. One cannot concatenate the 3 different quantities.
To accomblish this one may use an zip flag where the two lines are equivalent:
>>> DOS, PDOS, v = mp.apply.array.renew(zip=True).eigenstate(wrap=wrap_multiple, eta=True)
>>> DOS, PDOS, v = mp.apply(zip=True).array.eigenstate(wrap=wrap_multiple, eta=True)
and the data is unpacked as wanted.
Parallel calculations
The apply method looping k-points may be explicitly parallelized.
To run parallel do:
>>> H = si.Hamiltonian(...)
>>> mp = si.MonkhorstPack(H, [10, 10, 10])
>>> with mp.apply.renew(pool=True) as par:
... par.array.eigh()
This requires you also have the package pathos available.
The above will run in parallel using a default number of processors
in priority:
Environment variable
SISL_NUM_PROCSReturn value of
os.cpu_count().
Note that this may interfere with BLAS implementation which defaults
to use all CPU’s for threading. The total processors/threads that will
be created is SISL_NUM_PROCS * OMP_NUM_THREADS. Try and ensure this is below
the actual core-count of your machine (or the number of requested cores in a
HPC environment).
Alternatively one can control the number of processors locally by doing:
>>> H = si.Hamiltonian(...)
>>> mp = si.MonkhorstPack(H, [10, 10, 10])
>>> with mp.apply.renew(pool=2) as par:
... par.eigh()
which will request 2 processors (regardless of core-count).
As a last resort you can pass your own Pool of workers that
will be used for the parallel processing.
>>> from multiprocessing import Pool
>>> pool = Pool(4)
>>> H = si.Hamiltonian(...)
>>> mp = si.MonkhorstPack(H, [10, 10, 10])
>>> with mp.apply.renew(pool=pool) as par:
... par.array.eigh()
The Pool should implement some standard methods that are
existing in the pathos enviroment such as Pool.restart and Pool.terminate
and imap and uimap methods. See the pathos documentation for detalis.
A class to construct Brillouin zone related quantities |
|
Create a Monkhorst-Pack grid for the Brillouin zone |
|
Create a path in the Brillouin zone for plotting band-structures etc. |
Adaptive integration
Some integrands need many k-points in a small part of the Brillouin zone only,
for example a Green function with a small broadening.
An adaptive Brillouin zone selects the k-points from the values of the integrand.
It has no k-points before the integration, so len(bz) and bz.k
raise a TypeError. Use integrate, or bz.apply.average:
>>> H = si.Hamiltonian(...)
>>> gk = si.GaussKronrodGrid(H, [15, 0, 0], atol=1e-6)
>>> E = 0.1 + 1e-3j
>>> def G(k):
... return np.linalg.inv(E * H.Sk(k, format="array") - H.Hk(k, format="array"))
>>> Gavg = gk.integrate(G)
>>> D = gk.apply.average.Pk(format="array")
>>> bz = gk.last_bz()
The result is the Brillouin zone average (as for MonkhorstPack).
When G returns a list of values with different shapes, the result is an oplist,
and the convergence includes all values.
gauss selects the rule of a box (3, 7 or 10 Gauss nodes, with \(2g+1\) k-points
along a direction). In 1 dimension, gauss=10 often needs fewer k-points.
In 3 dimensions, gauss=3 can halve boxes with a smaller max_k.
bz.apply.average accepts wrap as for BrillouinZone, and the convergence
is on the wrapped value. The self-energies RealSpaceSE and RealSpaceSI accept an
adaptive Brillouin zone as their bz option.
The other bz.apply modes (array, iter, sum, …) need fixed k-points and raise a TypeError.
last_bz returns a BrillouinZone with the k-points
and weights of the last integration. It is only created when you call it, and you can
use it with bz.apply to calculate more quantities on the same k-points.
Base class for Brillouin zone integrations with adaptive k-point grids |
|
Adaptive Gauss-Kronrod cubature of the Brillouin zone |