WannierBerri advanced tutorial - Fermi surface and Fermi sea calculation of the Berry curvature dipole
author: Jae-Mo Lihm (jaemo.lihm@gmail.com)
The Berry curvature dipole (BCD) is the dipole of the Berry curvature in the reciprocal space. BCD has attracted considerable attention, partially because it contributes to the nonlinear Hall effect.
In WannierBerri, there are two ways to compute this quantity: the Fermi surface integral and the Fermi sea integral:
\[D^{\rm surf.}_{ab} = -\int_{\rm BZ} d\mathbf{k} \sum_n \left[ \frac{\partial}{\partial k^a} f_0 (\varepsilon_{n\mathbf{k}}) \right] \Omega^b_{n\mathbf{k}}\]
\[D^{\rm sea.}_{ab} = \int_{\rm BZ} d\mathbf{k} \sum_n f_0 (\varepsilon_{n\mathbf{k}}) \left(\frac{\partial}{\partial k^a} \Omega^b_{n\mathbf{k}} \right)\]
Formally, one can show that the two methods give identical results using partial integration:
\[D^{\rm sea.}_{ab} - D^{\rm surf.}_{ab}
= \int_{\rm BZ} d\mathbf{k} \frac{\partial}{\partial k^a} \left[ \sum_n f_0 (\varepsilon_{n\mathbf{k}}) \Omega^b_{n\mathbf{k}} \right]
= 0.\]
The integral is zero because the Brillouin zone has no boundary.
However, in practice, there are numerical differences between the two methods. In this tutorial, we compare the Fermi sea and Fermi surface integral methods for the BCD using the chiral tight-binding model and the ab-initio tight-binding model of ferroelectric GeTe.
[1]:
# Preliminary (Do only once)
%load_ext autoreload
%autoreload 2
# Set environment variables - not mandatory but recommended
import os
os.environ['OPENBLAS_NUM_THREADS'] = '1'
os.environ['MKL_NUM_THREADS'] = '1'
import numpy as np
import scipy
import matplotlib.pyplot as plt
import wannierberri as wberri
from wannierberri import calculators as calc
# 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
parallel = wberri.Parallel(num_cpus=1, progress_step_percent=10)
2023-08-12 02:17:11,717 INFO worker.py:1519 -- Started a local Ray instance. View the dashboard at 127.0.0.1:8267
Model 1: Chiral tight-binding model
Model, band structure
The model is taken from T. Yoda et al, “Orbital Edelstein Effect as a Condensed-Matter Analog of Solenoids”, Nano Lett. 18, 2, 916–920 (2018) https://arxiv.org/abs/1706.07702
The system is a honeycomb lattice with two orbitals per unit cell with an AA stacking (i.e. an AA stacked 2d hexagonal boron nitride). The original model contains in-plance nearest-neighbor hopping (hop1), vertical hopping (hopz_vert), and chiral inter-layer hopping (hopz_left and hopz_right). Here, we additionally add onsite energy difference (delta) and in-plane next-nearest-neighbor hopping (hop2 and phi). Nonzero phi breaks the time-reversal symmetry.
When we set phi=0, the model has time-reversal symmetry and C3z symmetry. In total, there are 6 symmetry operations.
[2]:
import wannierberri.models
# Initiate a tight-binding model (wannierberri.models.Chiral uses PythTB)
model_Chiral_left = wannierberri.models.Chiral(delta=2., hop1=1., hop2=0.3, phi=0, hopz_left=0.2, hopz_right=0.0, hopz_vert=0)
# Initialize WannierBerri system object
system = wberri.System_PythTB(model_Chiral_left, use_wcc_phase=True)
print("=" * 40)
# Set symmetry from the generators.
system.set_symmetry(["TimeReversal", "C3z"])
print("Number of symmetry operations: ", system.symgroup.size)
R=0 found at position(s) [[10]]
NOT using ws_dist
Number of wannier functions: 2
Number of R points: 21
Recommended size of FFT grid [3 3 3]
Reading the system from PythTB finished successfully
========================================
Number of symmetry operations: 6
Now, we plot the band structure along a k-point path that passes the high-symmetry k points. We use the plot_path_fat method to plot the band structure.
The energy can be accessed via path_result.get_data("Energy", iband). The valence band maximum is at the Gamma point, while the conduction band minimum on the K-H line.
Figure taken from T. Yoda et al, Nano Lett. 18, 2, 916–920 (2018).
[6]:
path = wberri.Path(
system,
k_nodes=[
[2/3, 1/3, 0.], # K
[0., 0., 0.], # Gamma
[1/2, 0., 0.], # M
[2/3, 1/3, 0.], # K
[2/3, 1/3, 0.5], # H
[0., 0., 0.5], # A
[1/2, 0., 0.5], # L
[2/3, 1/3, 0.5], # H
],
labels=["K","$\Gamma$","M","K", "H", "A", "L", "H"],
length=50
)
calculators = {
"tabulate": calc.TabulatorAll(
{"Energy": calc.tabulate.Energy()},
ibands=np.arange(system.num_wann),
mode="path",
)
}
path_result = wberri.run(
system,
grid=path,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)
path_result = path_result.results["tabulate"]
calculator not described
Calculation along a path - checking calculators for compatibility
tabulate <wannierberri.calculators.TabulatorAll object at 0x7f304d5baeb0>
All calculators are compatible
Symmetrization switched off for Path
The set of k points is a Path() with 184 points and labels {0: 'K', 33: '$\\Gamma$', 62: 'M', 79: 'K', 104: 'H', 137: 'A', 166: 'L', 183: 'H'}
WARNING : symmetry is not used for a tabulation along path
generating K_list
Done
Done, sum of weights:184.0
symgroup : None
processing 184 K points : using 1 processes.
# K-points calculated Wall time (sec) Est. remaining (sec)
18 0.2 1.8
36 0.2 0.8
54 0.3 0.7
72 0.4 0.6
90 0.5 0.5
108 0.6 0.4
126 0.7 0.3
144 0.7 0.2
162 0.8 0.1
180 0.9 0.0
time for processing 184 K-points on 1 processes: 0.9362 ; per K-point 0.0051 ; proc-sec per K-point 0.0051
time1 = 0.021979331970214844
Totally processed 184 K-points
[7]:
e_vbm = np.amax(path_result.get_data("Energy", 0))
e_cbm = np.amin(path_result.get_data("Energy", 1))
print(f"Valence band maximum : {e_vbm:6.3f}")
print(f"Conduction band minimum: {e_cbm:6.3f}")
print("VBM k point: ", path_result.kpoints[np.argmax(path_result.get_data("Energy", 0))])
print("CBM k point: ", path_result.kpoints[np.argmin(path_result.get_data("Energy", 1))])
fig = path_result.plot_path_fat(path, close_fig=False, show_fig=False)
ax = fig.get_axes()[0]
ax.axhline(e_vbm, c="k", ls="--")
ax.axhline(e_cbm, c="k", ls="--")
plt.show(fig)
Valence band maximum : -0.606
Conduction band minimum: -0.099
VBM k point: [0. 0. 0.]
CBM k point: [0.66666667 0.33333333 0.16 ]
Berry curvature dipole
[8]:
from wannierberri.smoother import FermiDiracSmoother
efermi = np.linspace(-2.0, 1.0, 101, True)
# efermi = np.linspace(e_vbm - 0.5, e_cbm + 0.5, 301, True)
kwargs = dict(
Efermi=efermi,
kwargs_formula={"external_terms": False},
smoother=FermiDiracSmoother(efermi, T_Kelvin=1200, maxdE=8)
)
calculators = {
'berry_dipole_fermi_sea': calc.static.BerryDipole_FermiSea(**kwargs),
'berry_dipole_fermi_surface': calc.static.BerryDipole_FermiSurf(**kwargs),
'dos': calc.static.DOS(Efermi=efermi,),
}
Berry curvature dipole (dimensionless)
| With Fermi sea integral. Eq(29) in `Ref <https://www.nature.com/articles/s41524-021-00498-5>`_
| Output: :math:`D_{\beta\delta} = \int [dk] \partial_\beta \Omega_\delta f`
Berry curvature dipole (dimensionless)
| With Fermi surface integral. Eq(8) in `Ref <https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.115.216806>`_
| Output: :math:`D_{\beta\delta} = -\int [dk] v_\beta \Omega_\delta f'`
Density of states
Compare calculations using the full k-point grid and the irreducible grid. This can be done using the use_irred_kpt argument of wberri.run.
You can find the following two lines in the output: K_list contains 1000 Irreducible points(100.0%) out of initial 10x10x10=1000 grid K_list contains 172 Irreducible points(17.2%) out of initial 10x10x10=1000 grid
These lines show that using symmetry can speed up calculation by 1000 / 172 ~ 6 times. This speedup corresponds to the 6 symmetry operations. Check that the BCD calculated from the two grids are identical.
[159]:
grid = wberri.Grid(system, NK=20, NKFFT=2)
print("=" * 40)
print("Calculation using the full k-point grid")
result_full_grid = wberri.run(
system,
grid=grid,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
use_irred_kpt = False,
)
print("=" * 40)
print("Calculation using the irreducible k-point grid")
result_irr_grid = wberri.run(
system,
grid=grid,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
use_irred_kpt = True,
)
determining grids from NK=20 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=2 (<class 'NoneType'>)
The grids were set to NKdiv=[10 10 10], NKFFT=[2 2 2], NKtot=[20 20 20]
========================================
Calculation using the full k-point grid
Grid is regular
The set of k points is a Grid() with NKdiv=[10 10 10], NKFFT=[2 2 2], NKtot=[20 20 20]
generating K_list
Done in 0.02189946174621582 s
Done in 0.02212214469909668 s
K_list contains 1000 Irreducible points(100.0%) out of initial 10x10x10=1000 grid
Done, sum of weights:1.0000000000000007
symgroup : <wannierberri.symmetry.Group object at 0x7f0be010e7f0>
processing 1000 K points : using 1 processes.
# K-points calculated Wall time (sec) Est. remaining (sec)
300 3.5 8.2
600 6.3 4.2
900 9.0 1.0
time for processing 1000 K-points on 1 processes: 16.7419 ; per K-point 0.0167 ; proc-sec per K-point 0.0167
time1 = 1.0496785640716553
Totally processed 1000 K-points
========================================
Calculation using the irreducible k-point grid
Grid is regular
The set of k points is a Grid() with NKdiv=[10 10 10], NKFFT=[2 2 2], NKtot=[20 20 20]
generating K_list
Done in 0.015077590942382812 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.11146044731140137 s
Done in 0.11165475845336914 s
K_list contains 172 Irreducible points(17.2%) out of initial 10x10x10=1000 grid
Done, sum of weights:1.0000000000000007
symgroup : <wannierberri.symmetry.Group object at 0x7f0be010e7f0>
processing 172 K points : using 1 processes.
# K-points calculated Wall time (sec) Est. remaining (sec)
52 0.8 1.7
104 1.3 0.9
156 2.0 0.2
time for processing 172 K-points on 1 processes: 3.4148 ; per K-point 0.0199 ; proc-sec per K-point 0.0199
time1 = 0.20007848739624023
Totally processed 172 K-points
[164]:
dir_string = ["D_xx", "D_yy", "D_zz"]
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for label, res in zip(["Full grid", "Irr. grid"], [result_full_grid, result_irr_grid]):
data_fsurf = res.results['berry_dipole_fermi_surface'].data
data_fsea = res.results['berry_dipole_fermi_sea'].data
for i in range(3):
if label == "Full grid":
fmt = f"C{i}-"
else:
fmt = "k--"
axes[0].plot(efermi, data_fsurf[:, i, i], fmt, label=label + " " + dir_string[i])
axes[1].plot(efermi, data_fsea[:, i, i], fmt, label=label + " " + dir_string[i])
axes[0].set_title("Fermi surface")
axes[1].set_title("Fermi sea")
for ax in axes:
ax.axhline(0, c="k", lw=1)
ax.set_xlabel("Fermi level (eV)")
ax.set_ylabel("Berry curvature dipole")
axes[1].legend()
plt.tight_layout()
plt.show()
From now on, we use the irreducible k-point grid, which is the default in wberri.run.
[166]:
grid = wberri.Grid(system, NK=20, NKFFT=2)
result = wberri.run(
system,
grid=grid,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)
determining grids from NK=20 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=2 (<class 'NoneType'>)
The grids were set to NKdiv=[10 10 10], NKFFT=[2 2 2], NKtot=[20 20 20]
Grid is regular
The set of k points is a Grid() with NKdiv=[10 10 10], NKFFT=[2 2 2], NKtot=[20 20 20]
generating K_list
Done in 0.018593788146972656 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.11524176597595215 s
Done in 0.11572003364562988 s
K_list contains 172 Irreducible points(17.2%) out of initial 10x10x10=1000 grid
Done, sum of weights:1.0000000000000007
symgroup : <wannierberri.symmetry.Group object at 0x7f0be010e7f0>
processing 172 K points : using 1 processes.
# K-points calculated Wall time (sec) Est. remaining (sec)
52 0.6 1.3
104 1.2 0.8
156 1.7 0.2
time for processing 172 K-points on 1 processes: 2.8150 ; per K-point 0.0164 ; proc-sec per K-point 0.0164
time1 = 0.2164595127105713
Totally processed 172 K-points
[176]:
data_fsurf = result.results['berry_dipole_fermi_surface'].data
data_fsea = result.results['berry_dipole_fermi_sea'].data
print("Shape of data: ", data_fsurf.shape)
print("Maximum value for each components (Fermi surface): ")
print(np.linalg.norm(data_fsurf, axis=0))
print("Maximum value for each components (Fermi sea): ")
print(np.linalg.norm(data_fsea, axis=0))
Shape of data: (101, 3, 3)
Maximum value for each components (Fermi surface):
[[4.79826162e-03 3.60947837e-19 2.68729962e-18]
[5.02055851e-19 4.79826162e-03 5.92648672e-19]
[3.87336309e-19 7.92158054e-20 1.07935885e-02]]
Maximum value for each components (Fermi sea):
[[4.65044414e-03 1.03203601e-19 3.69577356e-18]
[4.43330684e-19 4.65044414e-03 7.09898574e-19]
[3.16277692e-19 6.54682410e-20 9.30088828e-03]]
[177]:
dir_string = ["D_xx", "D_yy", "D_zz"]
for i in range(3):
plt.plot(efermi, data_fsurf[:, i, i], c=f"C{i}", label="Fermi surface " + dir_string[i])
plt.plot(efermi, data_fsea[:, i, i], c=f"C{i}", label="Fermi sea " + dir_string[i], ls="--")
# plt.plot(efermi, result.results['dos'].data, "k-")
plt.axhline(0, c="k", lw=1)
plt.xlabel("Fermi level")
plt.ylabel("Berry curvature dipole")
plt.axvline(e_vbm, c="k", lw=1, ls="--")
plt.axvline(e_cbm, c="k", lw=1, ls="--")
plt.legend()
plt.show()
The trace of the BCD tensor should be zero due to a topological reason (see footnote 2 of S.Tsirkin et al, PRB 97, 035158 (2018) for details).
The Fermi sea formula for the BCD tensor is explicitly traceless at every k point: this fact can be proven from the definition of the Fermi sea formula using
\[\sum_{a} \frac{\partial}{\partial k^a} \Omega^a_{n\mathbf{k}}
= \sum_{a} \epsilon_{abc} \frac{\partial}{\partial k^a} \frac{\partial}{\partial k^b} A^c_{n\mathbf{k}}
= 0.\]
In contrast, tracelessness is not explicit in the Fermi surface formula. The trace of the BCD tensor becomes zero only after summing over the full Brillouin zone due to the cancellation of contributions from all the k points.
Therefore, in numerical calculations using a finite number of k points, only the Fermi sea formula, not the Fermi surface formula, gives a zero trace. The trace will converge to 0 in the Fermi surface case as one makes the k-point grid denser.
[216]:
data_fsurf_trace = np.trace(data_fsurf, axis1=1, axis2=2)
data_fsea_trace = np.trace(data_fsea, axis1=1, axis2=2)
plt.plot(efermi, data_fsurf_trace, label="Fermi surface")
plt.plot(efermi, data_fsea_trace, label="Fermi sea")
plt.axhline(0, c="k", lw=1, zorder=1)
plt.xlabel("Fermi level")
plt.ylabel("Trace(D)")
plt.axvline(e_vbm, c="k", lw=1, ls="--")
plt.axvline(e_cbm, c="k", lw=1, ls="--")
plt.legend()
plt.show()
Try to calculate the Berry curvature dipole using the two methods for various grid sizes and see how they converge. Also, check that the trace of Berry curvature dipole indeed converges to 0 for both methods.
(NK=120 will be enough to see convergence.)
[178]:
data_fsurf_all = {}
data_fsea_all = {}
[198]:
for nk in [20, 40, 60, 80]:
grid = wberri.Grid(system, NK=nk)
result = wberri.run(
system,
grid=grid,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)
data_fsurf_all[nk] = result.results['berry_dipole_fermi_surface'].data
data_fsea_all[nk] = result.results['berry_dipole_fermi_sea'].data
determining grids from NK=120 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=None (<class 'NoneType'>)
Minimal symmetric FFT grid : [3 3 3]
The grids were set to NKdiv=[40 40 40], NKFFT=[3 3 3], NKtot=[120 120 120]
Grid is regular
The set of k points is a Grid() with NKdiv=[40 40 40], NKFFT=[3 3 3], NKtot=[120 120 120]
generating K_list
Done in 1.3817126750946045 s
excluding symmetry-equivalent K-points from initial grid
Done in 6.342855215072632 s
Done in 6.358941316604614 s
K_list contains 10682 Irreducible points(16.69%) out of initial 40x40x40=64000 grid
Done, sum of weights:1.0000000000000344
symgroup : <wannierberri.symmetry.Group object at 0x7f0be010e7f0>
processing 10682 K points : using 1 processes.
# K-points calculated Wall time (sec) Est. remaining (sec)
3205 60.7 141.7
6378 120.8 81.5
9583 176.0 20.2
time for processing 10682 K-points on 1 processes: 271.4751 ; per K-point 0.0254 ; proc-sec per K-point 0.0254
time1 = 11.60431694984436
Totally processed 10682 K-points
[208]:
ic = 0
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
nklist = sorted(list(data_fsurf_all.keys()))
for nk in nklist:
c = f"C{ic}"
axes[0].plot(efermi, data_fsurf_all[nk][:, 0, 0], c=c, label=f"NK={nk}")
axes[1].plot(efermi, data_fsea_all[nk][:, 0, 0], c=c, label=f"NK={nk}")
ic += 1
axes[0].set_title("Fermi surface")
axes[1].set_title("Fermi sea")
for ax in axes:
ax.axhline(0, c="k", lw=1)
ax.set_ylim([-1.5e-3, 1.5e-3])
ax.set_xlabel("Fermi level (eV)")
ax.set_ylabel("Berry curvature dipole")
axes[1].legend()
plt.tight_layout()
plt.show()
[200]:
ic = 0
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for nk in sorted(list(data_fsurf_all.keys())):
c = f"C{ic}"
axes[0].plot(efermi, np.trace(data_fsurf_all[nk], axis1=1, axis2=2), c=c, label=f"NK={nk}")
axes[1].plot(efermi, np.trace(data_fsea_all[nk], axis1=1, axis2=2), c=c, label=f"NK={nk}")
ic += 1
axes[0].set_title("Fermi surface")
axes[1].set_title("Fermi sea")
for ax in axes:
ax.axhline(0, c="k", lw=1)
ax.set_ylim([-1.5e-3, 1.5e-3])
ax.set_xlabel("Fermi level (eV)")
ax.set_ylabel("Trace of Berry curvature dipole")
axes[1].legend()
plt.tight_layout()
plt.show()
[206]:
nklist = sorted(list(data_fsurf_all.keys()))
plt.plot(nklist, [np.linalg.norm(np.trace(data_fsurf_all[nk], axis1=1, axis2=2)) for nk in nklist], "o-", label="Fermi surface")
plt.plot(nklist, [np.linalg.norm(np.trace(data_fsea_all[nk], axis1=1, axis2=2)) for nk in nklist], "o-", label="Fermi sea")
plt.legend()
plt.xlabel("NK")
plt.ylabel("Trace(D)")
# plt.yscale("log")
plt.show()
Model 2: Ferroelectric germanium telluride (α-GeTe)
The second model is an ab-initio tighit-binding model of ferroelectric GeTe. This material has a rhombohedrally distorted rocksalt structure where the Te atom is shifted from the center of the unit cell along the [111] direction.
Figure taken from H. Wang et al, npj Computational Materials, 6, 7 (2020).
We load the system from a Wannier90 output. We set symmetry using the set_symmetry_from_structure method, which calls spglib to automatically determine the symmetry of the system.
We also symmetrize the system. See data_GeTe/GeTe.win file and check that the initial projections are correct. For details, refer to the symmetrization tutorial.
[5]:
system_GeTe = wberri.System_w90("data_GeTe/GeTe", berry=True)
system_GeTe.set_structure([[0., 0., 0.], [0.4688262167, 0.4688262167, 0.4688262167]], ["Ge", "Te"])
system_GeTe.set_symmetry_from_structure()
system_GeTe.symmetrize(
proj=["Ge:s", "Ge:p", "Te:s", "Te:p"],
atom_name=["Ge", "Te"],
positions=[[0., 0., 0.], [0.4688262167, 0.4688262167, 0.4688262167]],
soc=False,
DFT_code='qe'
)
using fortio to read
Reading restart information from file ../data_GeTe/GeTe.chk :
Time to read .chk : 0.13194584846496582
Time for MMN.__init__() : 0.3212873935699463 , read : 0.31705737113952637 , headstring 0.004230022430419922
time for FFT_q_to_R : 0.010440587997436523 s
using ws_distance
irvec_new_all shape (37,)
using ws_dist for Ham_R
using ws_dist for AA_R
Number of wannier functions: 8
Number of R points: 63
Recommended size of FFT grid [3 3 3]
Real-space lattice:
[[ 2.11676981 -1.22211762 3.64031036]
[ 0. 2.44423524 3.64031036]
[-2.11676981 -1.22211762 3.64031036]]
Wannier atoms info
(1, 'Ge', [0.0, 0.0, 0.0], ['s', 'p'], [[0], [1, 2, 3]])
(2, 'Te', [0.4688262167, 0.4688262167, 0.4688262167], ['s', 'p'], [[4], [5, 6, 7]])
[get_spacegroup]
Spacegroup is R3m (160).
--------------- 1 ---------------
det = 0.9999999999999999
rotation: cart:
[ 1 0 0] [1.00 -0.00 -0.00]
[ 0 1 0] [-0.00 1.00 -0.00]
[ 0 0 1] [-0.00 0.00 1.00]
translation:
( 0.00000 0.00000 0.00000) ( 0.00000 0.00000 0.00000)
--------------- 2 ---------------
det = 1.0
rotation: cart:
[ 0 1 0] [-0.50 0.87 0.00]
[ 0 0 1] [-0.87 -0.50 -0.00]
[ 1 0 0] [-0.00 0.00 1.00]
translation:
( 0.00000 0.00000 0.00000) ( 0.00000 0.00000 0.00000)
--------------- 3 ---------------
det = 1.0
rotation: cart:
[ 0 0 1] [-0.50 -0.87 -0.00]
[ 1 0 0] [0.87 -0.50 -0.00]
[ 0 1 0] [-0.00 0.00 1.00]
translation:
( 0.00000 0.00000 0.00000) ( 0.00000 0.00000 0.00000)
--------------- 4 ---------------
det = -0.9999999999999999
rotation: cart:
[ 0 0 1] [-1.00 0.00 0.00]
[ 0 1 0] [-0.00 1.00 -0.00]
[ 1 0 0] [-0.00 0.00 1.00]
translation:
( 0.00000 0.00000 0.00000) ( 0.00000 0.00000 0.00000)
--------------- 5 ---------------
det = -1.0
rotation: cart:
[ 1 0 0] [0.50 -0.87 -0.00]
[ 0 0 1] [-0.87 -0.50 -0.00]
[ 0 1 0] [-0.00 0.00 1.00]
translation:
( 0.00000 0.00000 0.00000) ( 0.00000 0.00000 0.00000)
--------------- 6 ---------------
det = -1.0
rotation: cart:
[ 0 1 0] [0.50 0.87 0.00]
[ 1 0 0] [0.87 -0.50 -0.00]
[ 0 0 1] [-0.00 0.00 1.00]
translation:
( 0.00000 0.00000 0.00000) ( 0.00000 0.00000 0.00000)
##########################
Symmetrizing Start
rot = 1
rot = 2
rot = 3
rot = 4
rot = 5
rot = 6
number of symmetry oprations == 6
nRvec_add = 0
Symmetrizing Finished
Testing AA with diag = True
[0,0,0]
[-0. -0. -0. -0. -0. -0. -0. -0.]
[ 0. -0. 0. 0. 0. 0. 0. 0.]
==============================================
[ 0. 0. -0.1036 0.1036 0. 0. 0.1009 -0.1009]
[-0. -0. -0.1114 0.1191 0. 0. 0.0812 -0.0865]
==============================================
[ 0.0013 0.1771 -0.0008 -0.0008 5.1191 5.0203 5.1347 5.1347]
[ 0.0013 0.1771 -0.0009 -0.0007 5.1191 5.0203 5.1347 5.1346]
==============================================
[1,0,0]
[-0.0005 0.0209 0.0003 0.0115 0.0004 -0.0113 -0.0016 -0.0073]
[-0.0005 0.0209 -0.0009 0.0121 0.0004 -0.0113 -0.0009 -0.0074]
==============================================
[ 0.0003 -0.0121 -0.0056 -0.0012 -0.0002 0.0065 0.0037 0.0014]
[ 0.0003 -0.0121 -0.0064 -0.0009 -0.0002 0.0065 0.0045 0.0009]
==============================================
[ 0.0009 -0.0123 0.0021 -0.0057 -0.0007 0.0073 0.0005 0.0039]
[ 0.0009 -0.0123 0.0025 -0.0059 -0.0007 0.0073 -0.0003 0.0046]
==============================================
[0,1,0]
[-0. -0. -0. -0. 0. 0. 0. 0.]
[ 0. 0. -0. 0. 0. -0. -0. 0.]
==============================================
[-0.0005 0.0242 0.0144 -0.0007 0.0005 -0.0131 -0.0089 -0.0013]
[-0.0005 0.0242 0.0156 -0.0013 0.0005 -0.0131 -0.0092 -0.0016]
==============================================
[ 0.0009 -0.0123 -0.0095 0.006 -0.0007 0.0073 0.0057 -0.0013]
[ 0.0009 -0.0124 -0.0103 0.0065 -0.0007 0.0073 0.0071 -0.0026]
==============================================
[0,0,1]
[ 0.0005 -0.0209 -0.0003 -0.0115 -0.0004 0.0113 0.0016 0.0073]
[ 0.0005 -0.0209 0.0009 -0.0121 -0.0004 0.0113 0.0009 0.0074]
==============================================
[ 0.0003 -0.0121 -0.0056 -0.0012 -0.0002 0.0065 0.0037 0.0014]
[ 0.0003 -0.0121 -0.0064 -0.0009 -0.0002 0.0065 0.0045 0.0009]
==============================================
[ 0.0009 -0.0123 0.0021 -0.0057 -0.0007 0.0073 0.0005 0.0039]
[ 0.0009 -0.0123 0.0025 -0.0059 -0.0007 0.0073 -0.0003 0.0046]
==============================================
[6]:
path = wberri.Path(
system_GeTe,
k_nodes=[
[0.00000, 0.00000, 0.00000], # GAMMA
[0.50000, 0.00000, 0.00000], # L
[0.62910, 0.62910, 0.24180], # U
[0.50000, 0.50000, 0.00000], # X
[0.00000, 0.00000, 0.00000], # GAMMA
[0.50000, 0.50000, 0.50000], # Z
[0.62910, 0.62910, 0.24180], # U
],
labels=["$\Gamma$","L","U","X", "$\Gamma$", "Z", "U"],
length=50
)
calculators = {
"tabulate": calc.TabulatorAll(
{"Energy": calc.tabulate.Energy()},
ibands=np.arange(system_GeTe.num_wann),
mode="path",
)
}
path_result = wberri.run(
system_GeTe,
grid=path,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)
path_result = path_result.results["tabulate"]
calculator not described
Calculation along a path - checking calculators for compatibility
tabulate <wannierberri.calculators.TabulatorAll object at 0x7f535c4f24f0>
All calculators are compatible
Symmetrization switched off for Path
The set of k points is a Path() with 39 points and labels {0: '$\\Gamma$', 7: 'L', 15: 'U', 18: 'X', 26: '$\\Gamma$', 33: 'Z', 38: 'U'}
WARNING : symmetry is not used for a tabulation along path
generating K_list
Done
Done, sum of weights:39.0
symgroup : None
processing 39 K points : using 1 processes.
# K-points calculated Wall time (sec) Est. remaining (sec)
4 1.2 10.5
8 1.2 4.8
12 1.3 2.9
16 1.3 1.9
20 1.4 1.3
24 1.4 0.9
28 1.5 0.6
32 1.5 0.3
36 1.5 0.1
time for processing 39 K-points on 1 processes: 1.6132 ; per K-point 0.0414 ; proc-sec per K-point 0.0414
time1 = 0.03001260757446289
Totally processed 39 K-points
[7]:
# Calculate the VBM and CBM energies using a 30*30*30 grid.
calculators = {
"tabulate": calc.TabulatorAll(
{"Energy": calc.tabulate.Energy()},
ibands=np.arange(system_GeTe.num_wann),
mode="grid",
)
}
grid_energy_result = wberri.run(
system_GeTe,
grid=wberri.Grid(system_GeTe, NK=30),
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)
grid_energy_result = grid_energy_result.results["tabulate"]
e_vbm_GeTe = np.amax(grid_energy_result.get_data("Energy", 4))
e_cbm_GeTe = np.amin(grid_energy_result.get_data("Energy", 5))
print(f"Valence band maximum : {e_vbm_GeTe:6.3f}")
print(f"Conduction band minimum: {e_cbm_GeTe:6.3f}")
print("VBM k point: ", grid_energy_result.kpoints[np.argmax(grid_energy_result.get_data("Energy", 0))])
print("CBM k point: ", grid_energy_result.kpoints[np.argmin(grid_energy_result.get_data("Energy", 1))])
calculator not described
determining grids from NK=30 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=None (<class 'NoneType'>)
Minimal symmetric FFT grid : [3 3 3]
The grids were set to NKdiv=[10 10 10], NKFFT=[3 3 3], NKtot=[30 30 30]
Grid is regular
The set of k points is a Grid() with NKdiv=[10 10 10], NKFFT=[3 3 3], NKtot=[30 30 30]
generating K_list
Done in 0.020868301391601562 s
excluding symmetry-equivalent K-points from initial grid
Done in 0.20849919319152832 s
Done in 0.2090928554534912 s
K_list contains 116 Irreducible points(11.6%) out of initial 10x10x10=1000 grid
Done, sum of weights:1.0000000000000007
symgroup : <wannierberri.symmetry.Group object at 0x7f535ed4b640>
processing 116 K points : using 1 processes.
# K-points calculated Wall time (sec) Est. remaining (sec)
12 0.5 4.3
24 0.9 3.5
36 1.3 2.9
48 1.6 2.3
60 2.0 1.8
72 2.3 1.4
84 2.6 1.0
96 2.9 0.6
108 3.2 0.2
time for processing 116 K-points on 1 processes: 5.1376 ; per K-point 0.0443 ; proc-sec per K-point 0.0443
time1 = 0.1799149513244629
setting the grid
setting new kpoints
finding equivalent kpoints
collecting
collecting: to_grid : 0.2142021656036377
collecting: TABresult : 0.0028984546661376953
collecting - OK : 0.21726107597351074 (0.00016045570373535156)
using a pool of 4 processes to write txt frmsf of 54000 points per process
using a pool of 4 processes to write txt frmsf of 54000 points per process
using a pool of 4 processes to write txt frmsf of 54000 points per process
Totally processed 116 K-points
Valence band maximum : 9.226
Conduction band minimum: 9.836
VBM k point: [0.2 0.5 0.76666667]
CBM k point: [0. 0. 0.5]
[8]:
fig = path_result.plot_path_fat(path, close_fig=False, show_fig=False)
ax = fig.get_axes()[0]
ax.axhline(e_vbm_GeTe, c="k", ls="--")
ax.axhline(e_cbm_GeTe, c="k", ls="--")
plt.show(fig)
We compute the BCD tensor for the hole-doped case. (The electron-doped case is much harder to converge due to band crossings near the CBM.)
Try to calculate BCD for various grids and compare the convergence of the two methods. (NK=120 will be enough to see convergence.)
[22]:
data_GeTe_fsurf_all = {}
data_GeTe_fsea_all = {}
[10]:
from wannierberri import calculators as calc
from wannierberri.smoother import FermiDiracSmoother
efermi = np.linspace(e_vbm_GeTe - 1.5, e_cbm_GeTe, 301, True)
smoother = FermiDiracSmoother(efermi, 300.0)
kwargs = dict(
Efermi=efermi,
smoother=smoother,
)
calculators = {
'berry_dipole_fermi_sea': calc.static.BerryDipole_FermiSea(**kwargs),
'berry_dipole_fermi_surface': calc.static.BerryDipole_FermiSurf(**kwargs),
}
Berry curvature dipole (dimensionless)
| With Fermi sea integral. Eq(29) in `Ref <https://www.nature.com/articles/s41524-021-00498-5>`_
| Output: :math:`D_{\beta\delta} = \int [dk] \partial_\beta \Omega_\delta f`
Berry curvature dipole (dimensionless)
| With Fermi surface integral. Eq(8) in `Ref <https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.115.216806>`_
| Output: :math:`D_{\beta\delta} = -\int [dk] v_\beta \Omega_\delta f'`
[37]:
for nk in [20, 30, 40, 60, 80, 120]:
grid = wberri.Grid(system_GeTe, NK=nk)
result = wberri.run(
system_GeTe,
grid=grid,
calculators=calculators,
parallel=parallel,
print_Kpoints = False,
)
data_GeTe_fsurf_all[nk] = result.results['berry_dipole_fermi_surface'].dataSmooth
data_GeTe_fsea_all[nk] = result.results['berry_dipole_fermi_sea'].dataSmooth
determining grids from NK=80 (<class 'int'>), NKdiv=None (<class 'NoneType'>), NKFFT=None (<class 'NoneType'>)
Minimal symmetric FFT grid : [3 3 3]
The grids were set to NKdiv=[20 20 20], NKFFT=[4 4 4], NKtot=[80 80 80]
Grid is regular
The set of k points is a Grid() with NKdiv=[20 20 20], NKFFT=[4 4 4], NKtot=[80 80 80]
generating K_list
Done in 0.18623900413513184 s
excluding symmetry-equivalent K-points from initial grid
Done in 1.052736520767212 s
Done in 1.0575153827667236 s
K_list contains 781 Irreducible points(9.76%) out of initial 20x20x20=8000 grid
Done, sum of weights:0.999999999999985
symgroup : <wannierberri.symmetry.Group object at 0x7f535ed4b640>
processing 781 K points : using 1 processes.
# K-points calculated Wall time (sec) Est. remaining (sec)
78 27.0 242.9
156 48.1 192.7
234 70.9 165.6
312 99.5 149.6
390 121.4 121.7
468 141.7 94.8
546 162.2 69.8
624 182.4 45.9
702 204.4 23.0
780 224.4 0.3
time for processing 781 K-points on 1 processes: 240.9678 ; per K-point 0.3085 ; proc-sec per K-point 0.3085
time1 = 1.158562183380127
Totally processed 781 K-points
[34]:
nk = 20
print("Shape of data: ", data_GeTe_fsurf_all[nk].shape)
print("Maximum value for each components (Fermi surface): ")
print(np.amax(abs(data_GeTe_fsurf_all[nk]), axis=0))
print("Maximum value for each components (Fermi sea): ")
print(np.amax(abs(data_GeTe_fsea_all[nk]), axis=0))
Shape of data: (301, 3, 3)
Maximum value for each components (Fermi surface):
[[2.75821033e-18 1.78581033e-01 5.36068541e-15]
[1.78581033e-01 9.57931302e-18 1.09150709e-17]
[5.36012901e-15 1.09047721e-17 6.93608374e-18]]
Maximum value for each components (Fermi sea):
[[4.05362734e-19 5.76069248e-02 1.72783080e-15]
[5.76069248e-02 3.37626401e-18 1.77990429e-18]
[1.72920817e-15 1.38700899e-18 8.71348106e-19]]
[43]:
ic = 0
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for nk in sorted(list(data_GeTe_fsurf_all.keys())):
c = f"C{ic}"
axes[0].plot(efermi, data_GeTe_fsurf_all[nk][:, 0, 1], c=c, label=f"NK={nk}")
axes[1].plot(efermi, data_GeTe_fsea_all[nk][:, 0, 1], c=c, label=f"NK={nk}")
ic += 1
axes[0].set_title("Fermi surface")
axes[1].set_title("Fermi sea")
for ax in axes:
ax.axhline(0, c="k", lw=1)
ax.axvline(e_vbm_GeTe, c="k", lw=1, ls="--")
ax.axvline(e_cbm_GeTe, c="k", lw=1, ls="--")
ax.set_ylim([-0.05, 0.05])
ax.set_xlabel("Fermi level (eV)")
ax.set_ylabel("Berry curvature dipole")
axes[1].legend()
plt.tight_layout()
plt.show()
Let us focus on the undoped case: Fermi level between the VBM and the CBM. Then, $:nbsphinx-math:frac{partial}{partial k^a} f_0 (\varepsilon_{n:nbsphinx-math:mathbf{k}}) =0 $ holds so that the BCD is always 0. For the Fermi surface method, the integrand is 0 at every k point, so the integral is also 0. However, for the Fermi sea method, the integrand is not 0, and the integral becomes 0 only due to a calculation.
[44]:
ic = 0
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for nk in sorted(list(data_GeTe_fsurf_all.keys())):
c = f"C{ic}"
axes[0].plot(efermi, data_GeTe_fsurf_all[nk][:, 0, 1], c=c, label=f"NK={nk}")
axes[1].plot(efermi, data_GeTe_fsea_all[nk][:, 0, 1], c=c, label=f"NK={nk}")
ic += 1
axes[0].set_title("Fermi surface")
axes[1].set_title("Fermi sea")
for ax in axes:
ax.axhline(0, c="k", lw=1)
ax.axvline(e_vbm_GeTe, c="k", lw=1, ls="--")
ax.axvline(e_cbm_GeTe, c="k", lw=1, ls="--")
ax.set_xlabel("Fermi level (eV)")
ax.set_ylabel("Berry curvature dipole")
ax.set_ylim(np.array([-0.01, 0.03]) * 1)
ax.set_xlim([e_vbm_GeTe - 0.1, e_cbm_GeTe + 0.1])
axes[1].legend()
fig.suptitle("zoom around the band gap")
plt.tight_layout()
plt.show()
[47]:
# Choose the Fermi level index at the middle of the band gap.
ifermi = np.argmin(np.abs(efermi - (e_vbm_GeTe + e_cbm_GeTe) / 2))
nklist = sorted(list(data_GeTe_fsurf_all.keys()))
print("D_xy for the undoped case")
print("NK : ", nklist)
print("Fermi sea: ", [data_GeTe_fsea_all[nk][ifermi, 0, 1] for nk in nklist])
print("Fermi surface: ", [data_GeTe_fsurf_all[nk][ifermi, 0, 1] for nk in nklist])
plt.plot(nklist, [abs(data_GeTe_fsurf_all[NK][ifermi, 0, 1]) for NK in nklist], "o-")
plt.plot(nklist, [abs(data_GeTe_fsea_all[NK][ifermi, 0, 1]) for NK in nklist], "o-")
# plt.yscale("log") # Also try logscale
plt.xlabel("Nk")
plt.ylabel("BCD for the undoped case")
plt.show()
D_xy for the undoped case
NK : [20, 30, 40, 60, 80, 120]
Fermi sea: [0.005090029617682658, -0.001636927143447324, -7.100518562204471e-06, -3.727637728706892e-07, 5.031723634526716e-09, 1.0901526500907337e-12]
Fermi surface: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Further questions
If you are interested, try to answer the following questions: - What happens if one uses an un-symmetrized model? Does Fermi sea BCD with the Fermi energy inside the gap converge to 0? - What happens if one uses the tetrahedron method? - What happens for the zero-temperature case (you may use .data instead of .dataSmooth)?