• R/O
  • HTTP
  • SSH
  • HTTPS

python: Commit

libtetrabz python package


Commit MetaInfo

Revisionba17cd8b35cb58dc8ad61adfa5e66919477559bf (tree)
Zeit2021-11-10 02:00:31
AutorMitsuaki Kawamura <kawamitsuaki@gmai...>
CommiterMitsuaki Kawamura

Log Message

Add tutorial, manual, docstring

Ändern Zusammenfassung

Diff

--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
1+src/.ipynb_checkpoints/
2+src/sc.frmsf
3+dist/
4+*.pyc
5+doc/_build
6+src/libtetrabz.egg-info/
--- a/README.md
+++ b/README.md
@@ -1,9 +1,13 @@
11 # libtetrabz
22
3-Libtetrabz is a library which parform efficiently the Brillouin-zone
3+Libtetrabz is a library which perform efficiently the Brillouin-zone
44 integration in the electronic structure calculation in a solid by
55 using the tetrahedron method.
66
7+## Manual
8+
9+
10+
711 ## For citing libtetrabz
812
913 We would appreciate if you cite the following article in your
@@ -14,3 +18,7 @@ research with libtetrabz.
1418 M. Kawamura, Y. Gohda, and S. Tsuneyuki, Phys. Rev. B 89, 094515 (2014).
1519
1620 https://journals.aps.org/prb/abstract/10.1103/PhysRevB.89.094515
21+
22+## Git repository
23+
24+https://osdn.net/projects/libtetrabz/scm/git/python/
--- /dev/null
+++ b/doc/Makefile
@@ -0,0 +1,20 @@
1+# Minimal makefile for Sphinx documentation
2+#
3+
4+# You can set these variables from the command line, and also
5+# from the environment for the first two.
6+SPHINXOPTS ?=
7+SPHINXBUILD ?= sphinx-build
8+SOURCEDIR = .
9+BUILDDIR = _build
10+
11+# Put it first so that "make" without argument is like "make help".
12+help:
13+ @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
14+
15+.PHONY: help Makefile
16+
17+# Catch-all target: route all unknown targets to Sphinx using the new
18+# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS).
19+%: Makefile
20+ @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O)
--- /dev/null
+++ b/doc/api.rst
@@ -0,0 +1,517 @@
1+API usage
2+=========
3+
4+You can call a function in this library as follows:
5+
6+.. code-block:: python
7+
8+ import libtetrabz
9+..
10+
11+Total energy, charge density, occupations
12+-----------------------------------------
13+
14+.. math::
15+ \begin{align}
16+ \sum_{n}
17+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
18+ \theta(\varepsilon_{\rm F} -
19+ \varepsilon_{n k}) X_{n k}
20+ \end{align}
21+..
22+
23+.. code-block:: python
24+
25+ wght = libtetrabz.occ(bvec, eig)
26+..
27+
28+Parameters
29+
30+ .. code-block:: python
31+
32+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
33+ ..
34+
35+ Reciprocal lattice vectors (arbitrary unit).
36+ Because they are used to choose the direction of tetrahedra,
37+ only their ratio is used.
38+
39+ .. code-block:: python
40+
41+ eig = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
42+ ..
43+
44+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
45+ Where `ng0, ng1, ng2` are the number of grid for each direction of
46+ the reciprocal lattice vectors and `nb` is the number of bands.
47+
48+Return
49+
50+ .. code-block:: python
51+
52+ wght = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
53+ ..
54+
55+ The integration weights.
56+
57+Fermi energy and occupations
58+----------------------------
59+
60+.. math::
61+
62+ \begin{align}
63+ \sum_{n}
64+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
65+ \theta(\varepsilon_{\rm F} -
66+ \varepsilon_{n k}) X_{n k}
67+ \end{align}
68+
69+.. code-block:: python
70+
71+ ef, wght, iteration = libtetrabz.fermieng(bvec, eig, nelec)
72+..
73+
74+Parameters
75+
76+ .. code-block:: python
77+
78+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
79+ ..
80+
81+ Reciprocal lattice vectors (arbitrary unit).
82+ Because they are used to choose the direction of tetrahedra,
83+ only their ratio is used.
84+
85+ .. code-block:: python
86+
87+ eig = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
88+ ..
89+
90+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
91+ Where `ng0, ng1, ng2` are the number of grid for each direction of
92+ the reciprocal lattice vectors and `nb` is the number of bands.
93+
94+ .. code-block:: python
95+
96+ nelec = # float
97+ ..
98+
99+ The number of electrons.
100+
101+Return
102+
103+ .. code-block:: python
104+
105+ ef = # float
106+ ..
107+
108+ The Fermi energy. Unit is the same as `eig`.
109+
110+ .. code-block:: python
111+
112+ wght = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
113+ ..
114+
115+ The integration weights.
116+
117+Partial density of states
118+-------------------------
119+
120+.. math::
121+
122+ \begin{align}
123+ \sum_{n}
124+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
125+ \delta(\omega - \varepsilon_{n k})
126+ X_{n k}(\omega)
127+ \end{align}
128+
129+.. code-block:: python
130+
131+ wght = libtetrabz.occ(bvec, eig, e0)
132+..
133+
134+Parameters
135+
136+ .. code-block:: python
137+
138+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
139+ ..
140+
141+ Reciprocal lattice vectors (arbitrary unit).
142+ Because they are used to choose the direction of tetrahedra,
143+ only their ratio is used.
144+
145+ .. code-block:: python
146+
147+ eig = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
148+ ..
149+
150+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
151+ Where `ng0, ng1, ng2` are the number of grid for each direction of
152+ the reciprocal lattice vectors and `nb` is the number of bands.
153+
154+ .. code-block:: python
155+
156+ e0 = numpy.empty(ne, dtype=numpy.float_)
157+ ..
158+
159+ The energy point :math:`\omega` in the same unit of `eig`.
160+ Where `ne` is the number of energy points.
161+
162+Return
163+
164+ .. code-block:: python
165+
166+ wght = numpy.empty([ng0, ng1, ng2, nb, ne], dtype=numpy.float_)
167+ ..
168+
169+ The integration weights.
170+
171+Integrated density of states
172+----------------------------
173+
174+.. math::
175+
176+ \begin{align}
177+ \sum_{n}
178+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
179+ \theta(\omega - \varepsilon_{n k})
180+ X_{n k}(\omega)
181+ \end{align}
182+
183+.. code-block:: python
184+
185+ wght = libtetrabz.intdos(bvec, eig, e0)
186+..
187+
188+Parameters
189+
190+ .. code-block:: python
191+
192+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
193+ ..
194+
195+ Reciprocal lattice vectors (arbitrary unit).
196+ Because they are used to choose the direction of tetrahedra,
197+ only their ratio is used.
198+
199+ .. code-block:: python
200+
201+ eig = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
202+ ..
203+
204+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
205+ Where `ng0, ng1, ng2` are the number of grid for each direction of
206+ the reciprocal lattice vectors and `nb` is the number of bands.
207+
208+ .. code-block:: python
209+
210+ e0 = numpy.empty(ne, dtype=numpy.float_)
211+ ..
212+
213+ The energy point :math:`\omega` in the same unit of `eig`.
214+ Where `ne` is the number of energy points.
215+
216+Return
217+
218+ .. code-block:: python
219+
220+ wght = numpy.empty([ng0, ng1, ng2, nb, ne], dtype=numpy.float_)
221+ ..
222+
223+ The integration weights.
224+
225+Nesting function and Fr&oumlhlich parameter
226+-------------------------------------------
227+
228+.. math::
229+
230+ \begin{align}
231+ \sum_{n n'}
232+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
233+ \delta(\varepsilon_{\rm F} -
234+ \varepsilon_{n k}) \delta(\varepsilon_{\rm F} - \varepsilon'_{n' k})
235+ X_{n n' k}
236+ \end{align}
237+
238+.. code-block:: python
239+
240+ wght = libtetrabz.dbldelta(bvec, eig1, eig2)
241+..
242+
243+Parameters
244+
245+ .. code-block:: python
246+
247+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
248+ ..
249+
250+ Reciprocal lattice vectors (arbitrary unit).
251+ Because they are used to choose the direction of tetrahedra,
252+ only their ratio is used.
253+
254+ .. code-block:: python
255+
256+ eig1 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
257+ ..
258+
259+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
260+ Where `ng0, ng1, ng2` are the number of grid for each direction of
261+ the reciprocal lattice vectors and `nb` is the number of bands.
262+
263+ .. code-block:: python
264+
265+ eig2 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
266+ ..
267+
268+ Another energy eigenvalue :math:`\varepsilon'_{n k}` in the same unit of `eig1`.
269+ Typically it is assumed to be :math:`\varepsilon_{n k+q}`.
270+
271+Return
272+
273+ .. code-block:: python
274+
275+ wght = numpy.empty([ng0, ng1, ng2, nb, nb], dtype=numpy.float_)
276+ ..
277+
278+ The integration weights.
279+
280+A part of DFPT calculation
281+--------------------------
282+
283+.. math::
284+
285+ \begin{align}
286+ \sum_{n n'}
287+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
288+ \theta(\varepsilon_{\rm F} -
289+ \varepsilon_{n k}) \theta(\varepsilon_{n k} - \varepsilon'_{n' k})
290+ X_{n n' k}
291+ \end{align}
292+
293+.. code-block:: python
294+
295+ wght = libtetrabz.dblstep(bvec, eig1, eig2)
296+..
297+
298+Parameters
299+
300+ .. code-block:: python
301+
302+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
303+ ..
304+
305+ Reciprocal lattice vectors (arbitrary unit).
306+ Because they are used to choose the direction of tetrahedra,
307+ only their ratio is used.
308+
309+ .. code-block:: python
310+
311+ eig1 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
312+ ..
313+
314+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
315+ Where `ng0, ng1, ng2` are the number of grid for each direction of
316+ the reciprocal lattice vectors and `nb` is the number of bands.
317+
318+ .. code-block:: python
319+
320+ eig2 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
321+ ..
322+
323+ Another energy eigenvalue :math:`\varepsilon'_{n k}` in the same unit of `eig1`.
324+ Typically it is assumed to be :math:`\varepsilon_{n k+q}`.
325+
326+Return
327+
328+ .. code-block:: python
329+
330+ wght = numpy.empty([ng0, ng1, ng2, nb, nb], dtype=numpy.float_)
331+ ..
332+
333+ The integration weights.
334+
335+Static polarization function
336+----------------------------
337+
338+.. math::
339+
340+ \begin{align}
341+ \sum_{n n'}
342+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
343+ \frac{\theta(\varepsilon_{\rm F} - \varepsilon_{n k})
344+ \theta(\varepsilon'_{n' k} - \varepsilon_{\rm F})}
345+ {\varepsilon'_{n' k} - \varepsilon_{n k}}
346+ X_{n n' k}
347+ \end{align}
348+
349+.. code-block:: python
350+
351+ wght = libtetrabz.polstat(bvec, eig1, eig2)
352+..
353+
354+Parameters
355+
356+ .. code-block:: python
357+
358+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
359+ ..
360+
361+ Reciprocal lattice vectors (arbitrary unit).
362+ Because they are used to choose the direction of tetrahedra,
363+ only their ratio is used.
364+
365+ .. code-block:: python
366+
367+ eig1 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
368+ ..
369+
370+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
371+ Where `ng0, ng1, ng2` are the number of grid for each direction of
372+ the reciprocal lattice vectors and `nb` is the number of bands.
373+
374+ .. code-block:: python
375+
376+ eig2 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
377+ ..
378+
379+ Another energy eigenvalue :math:`\varepsilon'_{n k}` in the same unit of `eig1`.
380+ Typically it is assumed to be :math:`\varepsilon_{n k+q}`.
381+
382+Return
383+
384+ .. code-block:: python
385+
386+ wght = numpy.empty([ng0, ng1, ng2, nb, nb], dtype=numpy.float_)
387+ ..
388+
389+ The integration weights.
390+
391+Phonon linewidth
392+----------------
393+
394+.. math::
395+
396+ \begin{align}
397+ \sum_{n n'}
398+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
399+ \theta(\varepsilon_{\rm F} -
400+ \varepsilon_{n k}) \theta(\varepsilon'_{n' k} - \varepsilon_{\rm F})
401+ \delta(\varepsilon'_{n' k} - \varepsilon_{n k} - \omega)
402+ X_{n n' k}(\omega)
403+ \end{align}
404+
405+.. code-block:: python
406+
407+ wght = libtetrabz.fermigr(bvec, eig1, eig2, e0)
408+..
409+
410+Parameters
411+
412+ .. code-block:: python
413+
414+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
415+ ..
416+
417+ Reciprocal lattice vectors (arbitrary unit).
418+ Because they are used to choose the direction of tetrahedra,
419+ only their ratio is used.
420+
421+ .. code-block:: python
422+
423+ eig1 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
424+ ..
425+
426+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
427+ Where `ng0, ng1, ng2` are the number of grid for each direction of
428+ the reciprocal lattice vectors and `nb` is the number of bands.
429+
430+ .. code-block:: python
431+
432+ eig2 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
433+ ..
434+
435+ Another energy eigenvalue :math:`\varepsilon'_{n k}` in the same unit of `eig1`.
436+ Typically it is assumed to be :math:`\varepsilon_{n k+q}`.
437+
438+ .. code-block:: python
439+
440+ e0 = numpy.empty(ne, dtype=numpy.float_)
441+ ..
442+
443+ The energy point :math:`\omega` in the same unit of `eig`.
444+ Where `ne` is the number of energy points.
445+
446+Return
447+
448+ .. code-block:: python
449+
450+ wght = numpy.empty([ng0, ng1, ng2, nb, nb, ne], dtype=numpy.float_)
451+ ..
452+
453+ The integration weights.
454+
455+Polarization function (complex frequency)
456+-----------------------------------------
457+
458+.. math::
459+
460+ \begin{align}
461+ \sum_{n n'}
462+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
463+ \frac{\theta(\varepsilon_{\rm F} - \varepsilon_{n k})
464+ \theta(\varepsilon'_{n' k} - \varepsilon_{\rm F})}
465+ {\varepsilon'_{n' k} - \varepsilon_{n k} + i \omega}
466+ X_{n n' k}(\omega)
467+ \end{align}
468+
469+.. code-block:: python
470+
471+ wght = libtetrabz.polcmplex(bvec, eig1, eig2, e0)
472+..
473+
474+Parameters
475+
476+ .. code-block:: python
477+
478+ bvec = numpy.array([[b0x, b0y, b0z], [b1x, b1y, b1z], [b2x, b2y, b2z]])
479+ ..
480+
481+ Reciprocal lattice vectors (arbitrary unit).
482+ Because they are used to choose the direction of tetrahedra,
483+ only their ratio is used.
484+
485+ .. code-block:: python
486+
487+ eig1 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
488+ ..
489+
490+ The energy eigenvalue :math:`\varepsilon_{n k}` in arbitrary unit.
491+ Where `ng0, ng1, ng2` are the number of grid for each direction of
492+ the reciprocal lattice vectors and `nb` is the number of bands.
493+
494+ .. code-block:: python
495+
496+ eig2 = numpy.empty([ng0, ng1, ng2, nb], dtype=numpy.float_)
497+ ..
498+
499+ Another energy eigenvalue :math:`\varepsilon'_{n k}` in the same unit of `eig1`.
500+ Typically it is assumed to be :math:`\varepsilon_{n k+q}`.
501+
502+ .. code-block:: python
503+
504+ e0 = numpy.empty(ne, dtype=numpy.complex_)
505+ ..
506+
507+ The energy point :math:`\omega` in the same unit of `eig`.
508+ Where `ne` is the number of energy points.
509+
510+Return
511+
512+ .. code-block:: python
513+
514+ wght = numpy.empty([ng0, ng1, ng2, nb, nb, ne], dtype=numpy.complex)
515+ ..
516+
517+ The integration weights.
--- /dev/null
+++ b/doc/conf.py
@@ -0,0 +1,55 @@
1+# Configuration file for the Sphinx documentation builder.
2+#
3+# This file only contains a selection of the most common options. For a full
4+# list see the documentation:
5+# https://www.sphinx-doc.org/en/master/usage/configuration.html
6+
7+# -- Path setup --------------------------------------------------------------
8+
9+# If extensions (or modules to document with autodoc) are in another directory,
10+# add these directories to sys.path here. If the directory is relative to the
11+# documentation root, use os.path.abspath to make it absolute, like shown here.
12+#
13+import os
14+import sys
15+sys.path.insert(0, os.path.abspath('./'))
16+
17+
18+# -- Project information -----------------------------------------------------
19+
20+project = 'libtetrabz'
21+copyright = '2021, Mitsuaki Kawamura'
22+author = 'Mitsuaki Kawamura'
23+
24+# The full version, including alpha/beta/rc tags
25+release = '0.0.1'
26+
27+
28+# -- General configuration ---------------------------------------------------
29+
30+# Add any Sphinx extension module names here, as strings. They can be
31+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
32+# ones.
33+extensions = [
34+]
35+
36+# Add any paths that contain templates here, relative to this directory.
37+templates_path = ['_templates']
38+
39+# List of patterns, relative to source directory, that match files and
40+# directories to ignore when looking for source files.
41+# This pattern also affects html_static_path and html_extra_path.
42+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
43+
44+
45+# -- Options for HTML output -------------------------------------------------
46+
47+# The theme to use for HTML and HTML Help pages. See the documentation for
48+# a list of builtin themes.
49+#
50+html_theme = 'alabaster'
51+
52+# Add any paths that contain custom static files (such as style sheets) here,
53+# relative to this directory. They are copied after the builtin static files,
54+# so a file named "default.css" will overwrite the builtin "default.css".
55+html_static_path = ['_static']
\ No newline at end of file
--- /dev/null
+++ b/doc/index.rst
@@ -0,0 +1,14 @@
1+.. libtetrabz documentation master file, created by
2+ sphinx-quickstart on Mon Nov 8 16:12:42 2021.
3+ You can adapt this file completely to your liking, but it should at least
4+ contain the root `toctree` directive.
5+
6+Welcome to libtetrabz documentation!
7+====================================
8+
9+.. toctree::
10+ :maxdepth: 2
11+
12+ overview
13+ api
14+ ref
--- /dev/null
+++ b/doc/overview.rst
@@ -0,0 +1,94 @@
1+Introduction
2+============
3+
4+This document explains a tetrahedron method library ``libtetrabz``.
5+``libtetrabz`` is a library to calculate the total energy, the charge
6+density, partial density of states, response functions, etc. in a solid
7+by using the optimized tetrahedron method :ref:`[1] <ref>`.
8+Subroutines in this library receive the orbital (Kohn-Sham) energies as an input and
9+calculate weights :math:`w_{n n' k}` for integration such as
10+
11+.. math::
12+
13+ \begin{align}
14+ \sum_{n n'}
15+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
16+ F(\varepsilon_{n k}, \varepsilon_{n' k+q})X_{n n' k}
17+ = \sum_{n n'} \sum_{k}^{N_k} w_{n n' k} X_{n n' k}
18+ \end{align}
19+
20+``libtetrabz`` supports following Brillouin-zone integrations
21+
22+.. math::
23+
24+ \begin{align}
25+ \sum_{n}
26+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
27+ \theta(\varepsilon_{\rm F} - \varepsilon_{n k})
28+ X_{n k}
29+ \end{align}
30+
31+.. math::
32+
33+ \begin{align}
34+ \sum_{n}
35+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
36+ \delta(\omega - \varepsilon_{n k})
37+ X_{n k}(\omega)
38+ \end{align}
39+
40+.. math::
41+
42+ \begin{align}
43+ \sum_{n n'}
44+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
45+ \delta(\varepsilon_{\rm F} - \varepsilon_{n k})
46+ \delta(\varepsilon_{\rm F} - \varepsilon'_{n' k})
47+ X_{n n' k}
48+ \end{align}
49+
50+.. math::
51+
52+ \begin{align}
53+ \sum_{n n'}
54+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
55+ \theta(\varepsilon_{\rm F} - \varepsilon_{n k})
56+ \theta(\varepsilon_{n k} - \varepsilon'_{n' k})
57+ X_{n n' k}
58+ \end{align}
59+
60+.. math::
61+
62+ \begin{align}
63+ \sum_{n n'}
64+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
65+ \frac{
66+ \theta(\varepsilon_{\rm F} - \varepsilon_{n k})
67+ \theta(\varepsilon'_{n' k} - \varepsilon_{\rm F})}
68+ {\varepsilon'_{n' k} - \varepsilon_{n k}}
69+ X_{n n' k}
70+ \end{align}
71+
72+.. math::
73+
74+ \begin{align}
75+ \sum_{n n'}
76+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
77+ \theta(\varepsilon_{\rm F} - \varepsilon_{n k})
78+ \theta(\varepsilon'_{n' k} - \varepsilon_{\rm F})
79+ \delta(\varepsilon'_{n' k} - \varepsilon_{n k} - \omega)
80+ X_{n n' k}(\omega)
81+ \end{align}
82+
83+.. math::
84+
85+ \begin{align}
86+ \sum_{n n'}
87+ \int_{\rm BZ} \frac{d^3 k}{V_{\rm BZ}}
88+ \frac{
89+ \theta(\varepsilon_{\rm F} - \varepsilon_{n k})
90+ \theta(\varepsilon'_{n' k} - \varepsilon_{\rm F})}
91+ {\varepsilon'_{n' k} - \varepsilon_{n k} + i \omega}
92+ X_{n n' k}(\omega)
93+ \end{align}
94+
--- /dev/null
+++ b/doc/ref.rst
@@ -0,0 +1,6 @@
1+.. _ref:
2+
3+Reference
4+=========
5+
6+[1] `M. Kawamura, Y. Gohda, and S. Tsuneyuki, Phys. Rev. B 89, 094515 (2014). <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.89.094515>`_
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,7 +1,7 @@
11 [build-system]
22 requires = [
33 "setuptools>=42",
4- "wheel"
4+ "wheel",
55 "numpy"
66 ]
77 build-backend = "setuptools.build_meta"
--- a/src/libtetrabz/libtetrabz_common.py
+++ b/src/libtetrabz/libtetrabz_common.py
@@ -25,10 +25,11 @@ import numpy
2525
2626 def libtetrabz_initialize(ng, bvec):
2727 """
28- define shortest diagonal line & define type of tetragonal
2928 :param ng:
3029 :param bvec:
31- :return:
30+ :return: Weight
31+
32+ define shortest diagonal line & define type of tetragonal
3233 """
3334 bvec2 = numpy.empty([3, 3], dtype=numpy.float_)
3435 for i1 in range(3):
@@ -313,9 +314,12 @@ def libtetrabz_triangle_b1(e):
313314
314315
315316 def libtetrabz_triangle_b2(e):
316- """Cut triangle B2
317+ """
318+
317319 :param e:
318- :return:
320+ :returns:
321+ - x - first
322+ - y - second
319323 """
320324 a12 = (0.0 - e[2]) / (e[1] - e[2])
321325 a31 = (0.0 - e[1]) / (e[3] - e[1])
@@ -331,9 +335,13 @@ def libtetrabz_triangle_b2(e):
331335
332336
333337 def libtetrabz_triangle_c1(e):
334- """Cut triangle C1
335- :param e:
336- :return:
338+ """
339+ :param ndarray(4, float) e: energy eigenvalues
340+ :returns:
341+ Volume
342+ Vertex ratio
343+
344+ Cut triangle C1
337345 """
338346 a03 = (0.0 - e[3]) / (e[0] - e[3])
339347 a13 = (0.0 - e[3]) / (e[1] - e[3])
--- a/src/libtetrabz/libtetrabz_dos.py
+++ b/src/libtetrabz/libtetrabz_dos.py
@@ -26,11 +26,12 @@ from . import libtetrabz_common
2626
2727 def dos(bvec, eig, e0):
2828 """
29+ :param bvec: ndarray([3, 3], float) Reciprocal-lattice vectors
30+ :param eig: ndarray([ng0, ng1, ng2, nb], float) Energy eigenvalue
31+ :param e0: ndarray(ne, float) Energy grid
32+ :return : ndarray([ng0, ng1, ng2, nb, ne], float) Integration weight
33+
2934 Compute Dos : Delta(E - E1)
30- :param ndarray([3, 3], float) bvec: Reciprocal-lattice vectors
31- :param ndarray([ng0, ng1, ng2, nb], float) eig: Energy eigenvalue
32- :param ndarray(ne, float) e0: Energy grid
33- :return ndarray([ng0, ng1, ng2, nb, ne], float) wght: Weight
3435 """
3536 ng = numpy.array(eig.shape[0:3])
3637 nk = ng.prod(0)
--- /dev/null
+++ b/src/libtetrabz_tutorial.ipynb
@@ -0,0 +1,684 @@
1+{
2+ "cells": [
3+ {
4+ "cell_type": "markdown",
5+ "metadata": {
6+ "id": "pahorUv9vKeB"
7+ },
8+ "source": [
9+ "# LibTetraBZ tutorial\n",
10+ "This is tutorial for [LibtetraBZ](https://libtetrabz.osdn.jp/) which is a libraly to perform Brillouin-zone integrals with the optimized tetrahedron method [1].\n",
11+ "\n",
12+ "[1] \"Improved tetrahedron method for the Brillouin-zone integration applicable to response functions\", M. Kawamura, Y. Gohda, and S. Tsuneyuki, [Phys. Rev. B 89, 094515](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.89.094515) (2014).\n",
13+ "\n",
14+ "## 1, Preparation\n",
15+ "Before we start this tutorial, the nessesarry packages should be installed.\n",
16+ "LibTeteraBZ package is installed through the pip command."
17+ ]
18+ },
19+ {
20+ "cell_type": "code",
21+ "execution_count": null,
22+ "metadata": {
23+ "id": "J2SIkJdsDkU0"
24+ },
25+ "outputs": [],
26+ "source": [
27+ "!pip install libtetrabz"
28+ ]
29+ },
30+ {
31+ "cell_type": "markdown",
32+ "metadata": {
33+ "id": "dGU3CIgJipXD"
34+ },
35+ "source": [
36+ "## 2, Density of states of free electron\n",
37+ "Let us start the simplest case, namely the calculation of density of states of free electron system, $\\varepsilon_k = k^2/2$.\n",
38+ "$$\n",
39+ "D(\\varepsilon) = \\int d^3 k \\delta(\\epsilon - \\varepsilon_k) = \n",
40+ "4 \\pi \\int_0^\\infty dk k^2 \\delta \\left(\\epsilon - \\frac{k^2}{2}\\right)\n",
41+ "=4 \\pi \\sqrt{2 \\varepsilon}.\n",
42+ "$$\n",
43+ "As above, we already know the analytical result of $D(\\varepsilon)$.\n",
44+ "Therefore, it is appropriate to verify the numerical Brillouin-zone integration.\n",
45+ "\n",
46+ "### 2.1, Smearing method\n",
47+ "Since the delta function $\\delta(\\varepsilon)$ has infinity, in the simulation, it is replaced by a smeared function such as the gaussian,\n",
48+ "$$\n",
49+ "\\tilde{\\delta}(\\varepsilon)=\\frac{1}{\\sigma \\sqrt{\\pi}}\\exp\\left(-\\frac{\\varepsilon^2}{\\sigma^2}\\right)\n",
50+ "$$\n"
51+ ]
52+ },
53+ {
54+ "cell_type": "code",
55+ "execution_count": null,
56+ "metadata": {
57+ "id": "PMHCYZa6Dtkd"
58+ },
59+ "outputs": [],
60+ "source": [
61+ "import numpy"
62+ ]
63+ },
64+ {
65+ "cell_type": "code",
66+ "execution_count": null,
67+ "metadata": {
68+ "id": "D3sXKMQmsR_E"
69+ },
70+ "outputs": [],
71+ "source": [
72+ "k_max = 3\n",
73+ "e_max = k_max**2 * 0.5"
74+ ]
75+ },
76+ {
77+ "cell_type": "code",
78+ "execution_count": null,
79+ "metadata": {
80+ "id": "YTO3AWCNsb36"
81+ },
82+ "outputs": [],
83+ "source": [
84+ "bvec = 2.0 * numpy.array([[k_max, 0.0, 0.0],\n",
85+ " [0.0, k_max, 0.0],\n",
86+ " [0.0, 0.0, k_max]])"
87+ ]
88+ },
89+ {
90+ "cell_type": "code",
91+ "execution_count": null,
92+ "metadata": {
93+ "id": "gyBelNxVD4ld"
94+ },
95+ "outputs": [],
96+ "source": [
97+ "ng0 = 10\n",
98+ "nb = 1\n",
99+ "ng = numpy.array([ng0, ng0, ng0])\n",
100+ "eig = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)\n",
101+ "for i0 in range(ng[0]):\n",
102+ " for i1 in range(ng[1]):\n",
103+ " for i2 in range(ng[2]):\n",
104+ " ikvec = numpy.array([i0, i1, i2])\n",
105+ " for ii in range(3):\n",
106+ " if ikvec[ii] * 2 >= ng[ii]:\n",
107+ " ikvec[ii] += - ng[ii]\n",
108+ " kvec = ikvec[0:3] / ng[0:3]\n",
109+ " kvec[0:3] = kvec.dot(bvec)\n",
110+ " #\n",
111+ " eig[i0, i1, i2, 0] = 0.5 * kvec.dot(kvec)"
112+ ]
113+ },
114+ {
115+ "cell_type": "markdown",
116+ "metadata": {
117+ "id": "Bax0I1hrs1k4"
118+ },
119+ "source": [
120+ "$$\n",
121+ "D(\\varepsilon) = \\int d^3 k \\delta(\\varepsilon - \\varepsilon_k) \\approx\n",
122+ "\\sum_k \\frac{V_{BZ}}{N_k} \\tilde{\\delta}(\\varepsilon - \\varepsilon_k)\n",
123+ "$$"
124+ ]
125+ },
126+ {
127+ "cell_type": "code",
128+ "execution_count": null,
129+ "metadata": {},
130+ "outputs": [],
131+ "source": [
132+ "vbz = abs(numpy.linalg.det(bvec))"
133+ ]
134+ },
135+ {
136+ "cell_type": "markdown",
137+ "metadata": {},
138+ "source": [
139+ "$$\n",
140+ "\\sigma = \\left \\langle |\\nabla_k \\varepsilon_k| \\right \\rangle \\Delta k\n",
141+ "= \\left \\langle \\sqrt{2\\varepsilon_k} \\right \\rangle \\Delta k\n",
142+ "$$"
143+ ]
144+ },
145+ {
146+ "cell_type": "code",
147+ "execution_count": null,
148+ "metadata": {},
149+ "outputs": [],
150+ "source": [
151+ "sigma=numpy.sqrt(2.0*eig).mean()*2.0*k_max/ng0"
152+ ]
153+ },
154+ {
155+ "cell_type": "code",
156+ "execution_count": null,
157+ "metadata": {},
158+ "outputs": [],
159+ "source": [
160+ "e=numpy.linspace(0.0, e_max, 10)\n",
161+ "dos=numpy.empty(e.shape[0], dtype=numpy.float_)\n",
162+ "for ie in range(e.shape[0]):\n",
163+ " dos[ie] = numpy.exp(-(e[ie]-eig)**2/sigma**2).mean()*vbz/(sigma*numpy.sqrt(numpy.pi))"
164+ ]
165+ },
166+ {
167+ "cell_type": "code",
168+ "execution_count": null,
169+ "metadata": {},
170+ "outputs": [],
171+ "source": [
172+ "import matplotlib.pyplot as plt"
173+ ]
174+ },
175+ {
176+ "cell_type": "code",
177+ "execution_count": null,
178+ "metadata": {},
179+ "outputs": [],
180+ "source": [
181+ "e0=numpy.linspace(0.0, e_max, 100)\n",
182+ "dos0=4.0*numpy.pi*numpy.sqrt(2.0*e0)\n",
183+ "plt.plot(e0, dos0, label=\"Exact\")\n",
184+ "plt.plot(e, dos, label=\"Smearing\", linestyle=\"None\", marker=\"o\")\n",
185+ "plt.xlabel(\"Energy\")\n",
186+ "plt.ylabel(\"DOS\")\n",
187+ "plt.legend()\n",
188+ "plt.show()"
189+ ]
190+ },
191+ {
192+ "cell_type": "markdown",
193+ "metadata": {},
194+ "source": [
195+ "### 2.2, Tetrahedron method"
196+ ]
197+ },
198+ {
199+ "cell_type": "code",
200+ "execution_count": null,
201+ "metadata": {
202+ "id": "IjxVhoQXEAfW"
203+ },
204+ "outputs": [],
205+ "source": [
206+ "import libtetrabz"
207+ ]
208+ },
209+ {
210+ "cell_type": "code",
211+ "execution_count": null,
212+ "metadata": {},
213+ "outputs": [],
214+ "source": [
215+ "wght = libtetrabz.dos(bvec, eig, e)"
216+ ]
217+ },
218+ {
219+ "cell_type": "code",
220+ "execution_count": null,
221+ "metadata": {
222+ "id": "6Ne2FzwXEC4d"
223+ },
224+ "outputs": [],
225+ "source": [
226+ "dos_t = wght.sum(3).sum(2).sum(1).sum(0)*vbz"
227+ ]
228+ },
229+ {
230+ "cell_type": "code",
231+ "execution_count": null,
232+ "metadata": {},
233+ "outputs": [],
234+ "source": [
235+ "plt.plot(e0, dos0, label=\"Exact\")\n",
236+ "plt.plot(e, dos, label=\"Smearing\", linestyle=\"None\", marker=\"o\")\n",
237+ "plt.plot(e, dos_t, label=\"Tetrahedron\", linestyle=\"None\", marker=\"x\")\n",
238+ "plt.xlabel(\"Energy\")\n",
239+ "plt.ylabel(\"DOS\")\n",
240+ "plt.legend()\n",
241+ "plt.show()"
242+ ]
243+ },
244+ {
245+ "cell_type": "markdown",
246+ "metadata": {},
247+ "source": [
248+ "## 3, Lindhard function\n",
249+ "\n",
250+ "Next example is the calculation of polarization function. \n",
251+ "We again consider the simplest case, the free electrons system.\n",
252+ "$$\n",
253+ "\\chi(q) = \\int d^3 k \\frac{\\theta(\\varepsilon_{\\rm F} - \\varepsilon_k) - \\theta(\\varepsilon_{\\rm F} - \\varepsilon_{k+q})}{\\varepsilon_{k+q} - \\varepsilon_k}\n",
254+ "$$\n",
255+ "If we set $k_{\\rm F}=1$,\n",
256+ "$$\n",
257+ "\\chi(q) = 2 \\pi \\left(1 + \\frac{1 - q^2/4}{q} \\log\\left|\\frac{q+2}{q-2} \\right|\\right)\n",
258+ "$$"
259+ ]
260+ },
261+ {
262+ "cell_type": "code",
263+ "execution_count": null,
264+ "metadata": {
265+ "id": "lh1IHM8dENNW"
266+ },
267+ "outputs": [],
268+ "source": [
269+ "k_max = 1.5\n",
270+ "bvec = 2.0 * numpy.array([[k_max, 0.0, 0.0],\n",
271+ " [0.0, k_max, 0.0],\n",
272+ " [0.0, 0.0, k_max]])\n",
273+ "vbz = abs(numpy.linalg.det(bvec))"
274+ ]
275+ },
276+ {
277+ "cell_type": "code",
278+ "execution_count": null,
279+ "metadata": {
280+ "id": "9RBcyltzKpJD"
281+ },
282+ "outputs": [],
283+ "source": [
284+ "ng0 = 10\n",
285+ "nb = 1\n",
286+ "ng = numpy.array([ng0, ng0, ng0])\n",
287+ "eig = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)"
288+ ]
289+ },
290+ {
291+ "cell_type": "code",
292+ "execution_count": null,
293+ "metadata": {
294+ "id": "Ldzbq-9FKrBU"
295+ },
296+ "outputs": [],
297+ "source": [
298+ "ef = 0.5\n",
299+ "qmax = 4.0"
300+ ]
301+ },
302+ {
303+ "cell_type": "code",
304+ "execution_count": null,
305+ "metadata": {
306+ "id": "qLD0eBJ7KszC"
307+ },
308+ "outputs": [],
309+ "source": [
310+ "eig1 = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)\n",
311+ "eig2 = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)\n",
312+ "qx = numpy.arange(0.0, qmax, 0.4)\n",
313+ "chi_t = numpy.empty(qx.shape, dtype=numpy.float_)\n",
314+ "chi_s = numpy.empty(qx.shape, dtype=numpy.float_)"
315+ ]
316+ },
317+ {
318+ "cell_type": "code",
319+ "execution_count": null,
320+ "metadata": {},
321+ "outputs": [],
322+ "source": [
323+ "import math\n",
324+ "sigma=numpy.sqrt(2.0*(eig1+ef)).mean()*2.0*k_max/ng0\n",
325+ "for iq in range(qx.shape[0]):\n",
326+ " qvec = numpy.array([qx[iq], 0.0, 0.0])\n",
327+ " chi_s[iq] = 0.0\n",
328+ " for i0 in range(ng[0]):\n",
329+ " for i1 in range(ng[1]):\n",
330+ " for i2 in range(ng[2]):\n",
331+ " ikvec = numpy.array([i0, i1, i2])\n",
332+ " for ii in range(3):\n",
333+ " if ikvec[ii] * 2 >= ng[ii]:\n",
334+ " ikvec[ii] += - ng[ii]\n",
335+ " kvec = ikvec[0:3] / ng[0:3]\n",
336+ " kvec[0:3] = kvec.dot(bvec)\n",
337+ " eig1 = 0.5 * kvec.dot(kvec) - ef\n",
338+ " \n",
339+ " kvec[0:3] += qvec[0:3]\n",
340+ " eig2 = 0.5 * kvec.dot(kvec) - ef\n",
341+ " \n",
342+ " if abs(eig1 - eig2) < 1.0e-8:\n",
343+ " chi_s[iq] += numpy.exp(-eig1**2/sigma**2)/(sigma*numpy.sqrt(numpy.pi))\n",
344+ " else:\n",
345+ " chi_s[iq] += 0.5*(math.erfc(eig1/sigma) - math.erfc(eig2/sigma))/(eig2 - eig1)\n",
346+ " chi_s[iq] *= vbz / ng.prod(0)"
347+ ]
348+ },
349+ {
350+ "cell_type": "code",
351+ "execution_count": null,
352+ "metadata": {
353+ "id": "gs6Q1tFzKusG"
354+ },
355+ "outputs": [],
356+ "source": [
357+ "for iq in range(qx.shape[0]):\n",
358+ " qvec = numpy.array([qx[iq], 0.0, 0.0])\n",
359+ " for i0 in range(ng[0]):\n",
360+ " for i1 in range(ng[1]):\n",
361+ " for i2 in range(ng[2]):\n",
362+ " ikvec = numpy.array([i0, i1, i2])\n",
363+ " for ii in range(3):\n",
364+ " if ikvec[ii] * 2 >= ng[ii]:\n",
365+ " ikvec[ii] += - ng[ii]\n",
366+ " kvec = ikvec[0:3] / ng[0:3]\n",
367+ " kvec[0:3] = kvec.dot(bvec)\n",
368+ " eig1[i0, i1, i2, 0] = 0.5 * kvec.dot(kvec) - ef\n",
369+ " \n",
370+ " kvec[0:3] += qvec[0:3]\n",
371+ " eig2[i0, i1, i2, 0] = 0.5 * kvec.dot(kvec) - ef\n",
372+ " if iq == 0:\n",
373+ " e0 = numpy.array([0.0])\n",
374+ " wght = libtetrabz.dos(bvec, eig1, e0)\n",
375+ " chi_t[iq] = wght.sum() * vbz\n",
376+ " else:\n",
377+ " wght = libtetrabz.polstat(bvec, eig1, eig2)\n",
378+ " chi_t[iq] = wght.sum() * vbz * 2.0"
379+ ]
380+ },
381+ {
382+ "cell_type": "code",
383+ "execution_count": null,
384+ "metadata": {
385+ "id": "KQOcBHg5KzZq"
386+ },
387+ "outputs": [],
388+ "source": [
389+ "qx0 = numpy.arange(0.0, qmax, 0.011)\n",
390+ "chi0 = qx0.copy()\n",
391+ "for iq in range(qx0.shape[0]):\n",
392+ " if qx0[iq] < 1.0e-8:\n",
393+ " chi0[iq] = 4.0 * numpy.pi\n",
394+ " elif abs(qx0[iq] - 2.0) < 1.0e-8:\n",
395+ " chi0[iq] = 2.0 * numpy.pi\n",
396+ " else:\n",
397+ " chi0[iq] = 2.0 * numpy.pi*(1.0 + (1.0 - 0.25 * qx0[iq]**2) / qx0[iq] * numpy.log(numpy.abs((qx0[iq] + 2) / (qx0[iq] - 2))))"
398+ ]
399+ },
400+ {
401+ "cell_type": "code",
402+ "execution_count": null,
403+ "metadata": {
404+ "id": "IpGTnO8LK6lS"
405+ },
406+ "outputs": [],
407+ "source": [
408+ "plt.plot(qx0, chi0, label=\"Exact\")\n",
409+ "plt.plot(qx, chi_s, label=\"Smearing\", linestyle=\"None\", marker=\"o\")\n",
410+ "plt.plot(qx, chi_t, label=\"Tetrahedron\", linestyle=\"None\", marker=\"x\")\n",
411+ "plt.legend()\n",
412+ "plt.xlim(0)\n",
413+ "plt.ylim(0)\n",
414+ "plt.xlabel(\"$q / k_F$\")\n",
415+ "plt.ylabel(\"$\\chi(q)$\")\n",
416+ "plt.show()"
417+ ]
418+ },
419+ {
420+ "cell_type": "markdown",
421+ "metadata": {
422+ "id": "mo3GriM1K9UK"
423+ },
424+ "source": [
425+ "## 4, Multiband, low dimension, and partial DOS\n",
426+ "We caonsider the simple honeycomb lattice which has two orbitals (bands) in a unit cell."
427+ ]
428+ },
429+ {
430+ "cell_type": "code",
431+ "execution_count": null,
432+ "metadata": {},
433+ "outputs": [],
434+ "source": [
435+ "bvec = numpy.array([[1.0, 0.0, 0.0],\n",
436+ " [0.5, 0.5*numpy.sqrt(3.0), 0.0],\n",
437+ " [0.0, 0.0, 10.0]])"
438+ ]
439+ },
440+ {
441+ "cell_type": "code",
442+ "execution_count": null,
443+ "metadata": {},
444+ "outputs": [],
445+ "source": [
446+ "ham_indx = numpy.array([[0, 0, 0, 0, 1],\n",
447+ " [0, 0, 0, 1, 0],\n",
448+ " [-1, 0, 0, 0, 1],\n",
449+ " [1, 0, 0, 1, 0],\n",
450+ " [0, -1, 0, 0, 1],\n",
451+ " [0, 1, 0, 1, 0]])\n",
452+ "ham = numpy.array([-1.0, -1.0, -1.0, -1.0, -1.0, -1.0])"
453+ ]
454+ },
455+ {
456+ "cell_type": "code",
457+ "execution_count": null,
458+ "metadata": {},
459+ "outputs": [],
460+ "source": [
461+ "ng0 = 24\n",
462+ "nb = 2\n",
463+ "ng = numpy.array([ng0, ng0, 1])\n",
464+ "eig = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)\n",
465+ "proj = numpy.empty([ng[0], ng[1], ng[2], nb, nb], dtype=numpy.float_)\n",
466+ "wfc0 = numpy.array([[1.0, 1.0],\n",
467+ " [1.0, -1.0]]) / numpy.sqrt(2.0)\n",
468+ "for i0 in range(ng[0]):\n",
469+ " for i1 in range(ng[1]):\n",
470+ " for i2 in range(ng[2]):\n",
471+ " kvec = numpy.array([i0, i1, i2]) / ng[0:3]\n",
472+ " ham_k = numpy.zeros([nb, nb], dtype=numpy.complex_)\n",
473+ " for ih in range(ham.shape[0]):\n",
474+ " ham_k[ham_indx[ih, 3], ham_indx[ih, 4]] += ham[ih]*numpy.exp(2.0j*numpy.pi*kvec.dot(ham_indx[ih, 0:3]))\n",
475+ " eig[i0, i1, i2, :], ham_k = numpy.linalg.eigh(ham_k)\n",
476+ " ham_k[:, :] = wfc0.dot(ham_k)\n",
477+ " for ib in range(nb):\n",
478+ " for iwfc in range(2):\n",
479+ " proj[i0, i1, i2, ib, iwfc] = ham_k[iwfc, ib].real**2 + ham_k[iwfc, ib].imag**2"
480+ ]
481+ },
482+ {
483+ "cell_type": "code",
484+ "execution_count": null,
485+ "metadata": {
486+ "id": "2_6ggSgcD8ms"
487+ },
488+ "outputs": [],
489+ "source": [
490+ "e = numpy.linspace(-3, 3, 100)"
491+ ]
492+ },
493+ {
494+ "cell_type": "code",
495+ "execution_count": null,
496+ "metadata": {},
497+ "outputs": [],
498+ "source": [
499+ "wght = libtetrabz.dos(bvec, eig, e)"
500+ ]
501+ },
502+ {
503+ "cell_type": "code",
504+ "execution_count": null,
505+ "metadata": {},
506+ "outputs": [],
507+ "source": [
508+ "dos = wght.sum(3).sum(2).sum(1).sum(0)\n",
509+ "pdos = numpy.zeros([e.shape[0], 2], dtype=numpy.float_)\n",
510+ "for ie in range(e.shape[0]):\n",
511+ " for iwfc in range(2):\n",
512+ " pdos[ie, iwfc] = (wght[:,:,:,:,ie]*proj[:,:,:,:,iwfc]).sum()"
513+ ]
514+ },
515+ {
516+ "cell_type": "code",
517+ "execution_count": null,
518+ "metadata": {},
519+ "outputs": [],
520+ "source": [
521+ "plt.plot(e, dos, label=\"Total\")\n",
522+ "plt.plot(e, pdos[:, 0], label=\"Bonding\")\n",
523+ "plt.plot(e, pdos[:, 1], label=\"Anti\")\n",
524+ "plt.xlabel(\"Energy\")\n",
525+ "plt.ylabel(\"DOS\")\n",
526+ "plt.ylim(0)\n",
527+ "plt.legend()\n",
528+ "plt.show()"
529+ ]
530+ },
531+ {
532+ "cell_type": "markdown",
533+ "metadata": {},
534+ "source": [
535+ "## 5, Simple-cubic tight-binding model, Fermi energy and Fermi surface"
536+ ]
537+ },
538+ {
539+ "cell_type": "code",
540+ "execution_count": null,
541+ "metadata": {},
542+ "outputs": [],
543+ "source": [
544+ "bvec=numpy.array([[10.0,0.0,0.0],\n",
545+ " [0.0,10.0,0.0],\n",
546+ " [0.0,0.0,10.0]])"
547+ ]
548+ },
549+ {
550+ "cell_type": "code",
551+ "execution_count": null,
552+ "metadata": {},
553+ "outputs": [],
554+ "source": [
555+ "ng0 = 10\n",
556+ "nb = 1\n",
557+ "ng = numpy.array([ng0, ng0, ng0])\n",
558+ "eig = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)\n",
559+ "vf = numpy.empty([ng[0], ng[1], ng[2], nb], dtype=numpy.float_)\n",
560+ "for i0 in range(ng[0]):\n",
561+ " for i1 in range(ng[1]):\n",
562+ " for i2 in range(ng[2]):\n",
563+ " kvec = 2.0*numpy.pi*numpy.array([i0, i1, i2]) / ng[0:3]\n",
564+ " eig[i0, i1, i2, 0] = - numpy.cos(kvec).sum(0)\n",
565+ " vf[i0, i1, i2, 0] = numpy.sqrt((numpy.sin(kvec)**2).sum(0))"
566+ ]
567+ },
568+ {
569+ "cell_type": "code",
570+ "execution_count": null,
571+ "metadata": {},
572+ "outputs": [],
573+ "source": [
574+ "e=numpy.linspace(-3, 3, 100)"
575+ ]
576+ },
577+ {
578+ "cell_type": "code",
579+ "execution_count": null,
580+ "metadata": {},
581+ "outputs": [],
582+ "source": [
583+ "wght=libtetrabz.dos(bvec,eig,e)"
584+ ]
585+ },
586+ {
587+ "cell_type": "code",
588+ "execution_count": null,
589+ "metadata": {},
590+ "outputs": [],
591+ "source": [
592+ "dos=wght.sum(3).sum(2).sum(1).sum(0)"
593+ ]
594+ },
595+ {
596+ "cell_type": "code",
597+ "execution_count": null,
598+ "metadata": {},
599+ "outputs": [],
600+ "source": [
601+ "plt.plot(e, dos)\n",
602+ "plt.xlabel(\"Energy\")\n",
603+ "plt.ylabel(\"DOS\")\n",
604+ "plt.ylim(0)\n",
605+ "plt.show()"
606+ ]
607+ },
608+ {
609+ "cell_type": "code",
610+ "execution_count": null,
611+ "metadata": {},
612+ "outputs": [],
613+ "source": [
614+ "ef, wght, iteration=libtetrabz.fermieng(bvec, eig, 0.5)"
615+ ]
616+ },
617+ {
618+ "cell_type": "code",
619+ "execution_count": null,
620+ "metadata": {},
621+ "outputs": [],
622+ "source": [
623+ "wght.sum(), ef, iteration"
624+ ]
625+ },
626+ {
627+ "cell_type": "code",
628+ "execution_count": null,
629+ "metadata": {},
630+ "outputs": [],
631+ "source": [
632+ "with open(\"sc.frmsf\", \"w\") as f:\n",
633+ " print(ng[0], ng[1], ng[2], file=f)\n",
634+ " print(1, file=f)\n",
635+ " print(nb, file=f)\n",
636+ " for ii in range(3):\n",
637+ " print(bvec[ii, 0], bvec[ii, 1], bvec[ii, 2], file=f)\n",
638+ " for ib in range(nb):\n",
639+ " for i0 in range(ng[0]):\n",
640+ " for i1 in range(ng[1]):\n",
641+ " for i2 in range(ng[2]):\n",
642+ " print(eig[i0, i1, i2, ib]-ef, file=f)\n",
643+ " for ib in range(nb):\n",
644+ " for i0 in range(ng[0]):\n",
645+ " for i1 in range(ng[1]):\n",
646+ " for i2 in range(ng[2]):\n",
647+ " print(vf[i0, i1, i2, ib], file=f)"
648+ ]
649+ },
650+ {
651+ "cell_type": "code",
652+ "execution_count": null,
653+ "metadata": {},
654+ "outputs": [],
655+ "source": []
656+ }
657+ ],
658+ "metadata": {
659+ "colab": {
660+ "collapsed_sections": [],
661+ "name": "libtetrabz_tutorial.ipynb",
662+ "provenance": []
663+ },
664+ "kernelspec": {
665+ "display_name": "Python 3 (ipykernel)",
666+ "language": "python",
667+ "name": "python3"
668+ },
669+ "language_info": {
670+ "codemirror_mode": {
671+ "name": "ipython",
672+ "version": 3
673+ },
674+ "file_extension": ".py",
675+ "mimetype": "text/x-python",
676+ "name": "python",
677+ "nbconvert_exporter": "python",
678+ "pygments_lexer": "ipython3",
679+ "version": "3.7.8"
680+ }
681+ },
682+ "nbformat": 4,
683+ "nbformat_minor": 1
684+}
--- a/doc/tutorial.ipynb
+++ b/src/tutorial.ipynb
@@ -222,9 +222,17 @@
222222 ]
223223 },
224224 {
225+ "cell_type": "markdown",
226+ "id": "564dd3eb",
227+ "metadata": {},
228+ "source": [
229+ "$\\varepsilon$"
230+ ]
231+ },
232+ {
225233 "cell_type": "code",
226234 "execution_count": null,
227- "id": "11206e36",
235+ "id": "9fcff06e",
228236 "metadata": {},
229237 "outputs": [],
230238 "source": []
Show on old repository browser