diff --git a/.gitignore b/.gitignore index deae2a3..b3eaab6 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,6 @@ tmp* # Mac stuff *.DS_Store* + +# PYPI test tokens +.pypirc diff --git a/README.md b/README.md index fb52fed..f9e167c 100644 --- a/README.md +++ b/README.md @@ -1,34 +1,71 @@ # Relativistic magnetic interactions from non-orthogonal basis sets +More on the theoretical background can be seen on [arXiv](https://arxiv.org/abs/2309.02558). + # TODO -- Definition of magnetic entities: - * Through simple sequence o forbitals in the unit cell - * Through atom specification - * Through atom and orbital specification -- Separation of TR and TRB components of the Hamiltonian, Identification of the exchange field. -- Definition of commutator expressions, old projection matrix elements -- Efficient calculation of Green's functions -- Calculation of energy and momentum resolved derivatives -- Parallel BZ and serial energy integral +- Parallel or serial energy integral to reduce memory overhead +- Check the symmetrization of the Hamiltonian and overlap matrix to make them hermitian +- Check if exchange field has scalar part +- Add more tests +- Run tests on different magnetic materials and compare it to Grogu Matlab. +- Ehm MAKE IT WORK SOMEHOW :'( # Building wheel See detailed documentation on [PYPI](https://packaging.python.org/en/latest/tutorials/packaging-projects/). +- First you need some API Tokens for Test PYPI ,to be able to upload. You can read about it [here](https://test.pypi.org/help/#apitoken). I own the current project, so you have to contact me. + + Use the following commands for a quick setup: - Build wheel + ``` python -m build +``` + +- Push to PYPI test repository +``` +python -m twine upload --repository testpypi dist/* ``` -Build wheel: +# Usage -Push to pypi(testpypi for beginners): python3 -m twine upload --repository testpypi dist/* +## For end users +Download and install from [PYPI](https://test.pypi.org/project/grogu-magn/) test repository. -Ehhez kellenek tokenek: -You will be prompted for a username and password. For the username, use __token__. For the password, use the token value, including the pypi- prefix. +``` +python3 -m pip install --index-url https://test.pypi.org/simple/ grogu_magn +``` +## For developers -Végfelhasználóknak (egyelőre testpypi): python3 -m pip install --index-url https://test.pypi.org/simple/ example-package-YOUR-USERNAME-HERE +- Clone repository from Gitea + +``` +git clone https://gitea.vo.elte.hu/et209d/grogu.git +``` + +- Create .venv (for example with VsCode) + * Use python 3.9.6 + * install dependencies from: + * requirements.txt + * requirements-dev.txt + * /docs/requirements.txt + +- Install and run pre-commit + +``` +pre-commit install +pre-commit run --all-files +``` + +To build the documentation navigate to the `docs/source` folder and run `make clean` and `make html`. After this the html page can be found in `docs/source/_build/html`. + +``` +cd docs/source +make clean +make html +``` diff --git a/docs/source/_build/doctrees/environment.pickle b/docs/source/_build/doctrees/environment.pickle index 74e6848..2b1ac65 100644 Binary files a/docs/source/_build/doctrees/environment.pickle and b/docs/source/_build/doctrees/environment.pickle differ diff --git a/docs/source/_build/doctrees/grogu.doctree b/docs/source/_build/doctrees/grogu.doctree index a099c16..a4f2663 100644 Binary files a/docs/source/_build/doctrees/grogu.doctree and b/docs/source/_build/doctrees/grogu.doctree differ diff --git a/docs/source/_build/doctrees/index.doctree b/docs/source/_build/doctrees/index.doctree index e84e375..1baf1f4 100644 Binary files a/docs/source/_build/doctrees/index.doctree and b/docs/source/_build/doctrees/index.doctree differ diff --git a/docs/source/_build/doctrees/modules.doctree b/docs/source/_build/doctrees/modules.doctree index 93edda6..4a65c60 100644 Binary files a/docs/source/_build/doctrees/modules.doctree and b/docs/source/_build/doctrees/modules.doctree differ diff --git a/docs/source/_build/html/.buildinfo b/docs/source/_build/html/.buildinfo index ed21c1e..7580025 100644 --- a/docs/source/_build/html/.buildinfo +++ b/docs/source/_build/html/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 -# This file records the configuration used when building these files. When it is not found, a full rebuild will be done. -config: 0f8ff1598a5b65221fc2ef8f29fdbeac +# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. +config: d8f86c3651bae6bbb2cfa348920c9196 tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/docs/source/_build/html/_modules/grogu/useful.html b/docs/source/_build/html/_modules/grogu/useful.html index 15af420..8407880 100644 --- a/docs/source/_build/html/_modules/grogu/useful.html +++ b/docs/source/_build/html/_modules/grogu/useful.html @@ -13,7 +13,7 @@ - + @@ -89,15 +89,22 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. +from itertools import permutations, product + import numpy as np from scipy.special import roots_legendre -from itertools import permutations, product + +# Pauli matrices +tau_x = np.array([[0, 1], [1, 0]]) +tau_y = np.array([[0, -1j], [1j, 0]]) +tau_z = np.array([[1, 0], [0, -1]]) +tau_0 = np.array([[1, 0], [0, 1]]) # define some useful functions
[docs] -def hsk(dh, k=(0, 0, 0)): +def hsk(H, ss, sc_off, k=(0, 0, 0)): """ One way to speed up Hk and Sk generation """ @@ -105,13 +112,12 @@ k.shape = (-1,) # are from the sisl source # this generates the list of phases - phases = np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell), dh.sc.sc_off.T)) + phases = np.exp(-1j * 2 * np.pi * k @ sc_off.T) - HKU = np.einsum("abc,c->ab", dh.hup, phases) - HKD = np.einsum("abc,c->ab", dh.hdo, phases) - SK = np.einsum("abc,c->ab", dh.sov, phases) + HK = np.einsum("abc,a->bc", H, phases) + SK = np.einsum("abc,a->bc", ss, phases) - return HKU, HKD, SK
+ return HK, SK @@ -175,80 +181,190 @@ -
-[docs] -def make_atran(nauc, dirs="xyz", dist=1): +
+[docs] +def commutator(a, b): + "Shorthand for commutator" + return a @ b - b @ a
+ + + +
+[docs] +def tau_u(u): """ - Simple pair generator. Depending on the value of the dirs - argument sampling in 1,2 or 3 dimensions is generated. - If dirs argument does not contain either of x,y or z - a single pair is returend. + Pauli matrix in direction u. """ - if not (sum([d in dirs for d in "xyz"])): - return (0, 0, [1, 0, 0]) + u = u / np.linalg.norm(u) # u is force to be of unit length + return u[0] * tau_x + u[1] * tau_y + u[2] * tau_z
- dran = len(dirs) * [np.arange(-dist, dist + 1)] - mg = np.meshgrid(*dran) - dirsdict = dict() - for d in enumerate(dirs): - dirsdict[d[1]] = mg[d[0]].flatten() - for d in "xyz": - if not (d in dirs): - dirsdict[d] = 0 * dirsdict[dirs[0]] - ucran = np.array([dirsdict[d] for d in "xyz"]).T - atran = [] - for i, j in list(product(range(nauc), repeat=2)): - for u in ucran: - if (abs(i - j) + sum(abs(u))) > 0: - atran.append((i, j, list(u))) +# +
+[docs] +def crossM(u): + """ + Definition for the cross-product matrix. + Acting as a cross product with vector u. + """ + return np.array([[0, -u[2], u[1]], [u[2], 0, -u[0]], [-u[1], u[0], 0]])
+ - return atran
+
+[docs] +def RotM(theta, u, eps=1e-10): + """ + Definition of rotation matrix with angle theta around direction u. + """ + u = u / np.linalg.norm(u) + + M = ( + np.cos(theta) * np.eye(3) + + np.sin(theta) * crossM(u) + + (1 - np.cos(theta)) * np.outer(u, u) + ) + M[abs(M) < eps] = 0.0 # kill off small numbers + return M
-
-[docs] -def add(x, y): - """The sum of two numbers for testing. - This function adds to numbers together. I only created this for testing documentation and examples. +
+[docs] +def RotMa2b(a, b, eps=1e-10): + """ + Definition of rotation matrix rotating unit vector a to unit vector b. + Function returns array R such that R@a = b holds. + """ + v = np.cross(a, b) + c = a @ b + M = np.eye(3) + crossM(v) + crossM(v) @ crossM(v) / (1 + c) + M[abs(M) < eps] = 0.0 # kill off small numbers + return M
- Parameters - ---------- - x : float - First number - y : float - Second number added to `x` - Returns - ------- - sum : int - The sum of the inputs - See Also - -------- - numpy.add : Adds more than two numbers. +
+[docs] +def spin_tracer(M): + """ + Spin tracer utility. + This akes an operator with the orbital-spin sequence: + orbital 1 up, + orbital 1 down, + orbital 2 up, + orbital 2 down, + that is in the SPIN-BOX representation, + and extracts orbital dependent Pauli traces. + """ - Notes - ----- - We can create some latex notes here [1]_ : + M11 = M[0::2, 0::2] + M12 = M[0::2, 1::2] + M21 = M[1::2, 0::2] + M22 = M[1::2, 1::2] - .. math:: a + b = c + M_o = dict() + M_o["x"] = M12 + M21 + M_o["y"] = 1j * (M12 - M21) + M_o["z"] = M11 - M22 + M_o["c"] = M11 + M22 - References - ---------- - .. [1] https://numpydoc.readthedocs.io/en/latest/format.html + return M_o
- Examples - -------- - >>> add(1, 2) - 3 +
+[docs] +def parse_magnetic_entity(dh, atom=None, l=None, **kwargs): + """ + Function to define orbital indeces of a given magnetic entity. + dh: a sisl Hamiltonian object + atom: an integer or list of integers, defining atom (or atoms) in the unicell forming the magnetic entity + l: integer, defining the angular momentum channel + """ + # case where we deal with more than one atom defining the magnetic entity + if type(atom) == list: + dat = [] + for a in atom: + a_orb_idx = dh.geometry.a2o(a, all=True) + if ( + type(l) == int + ): # if specified we restrict to given l angular momentum channel inside each atom + a_orb_idx = a_orb_idx[[o.l == l for o in dh.geometry.atoms[a].orbitals]] + dat.append(a_orb_idx) + orbital_indeces = np.hstack(dat) + # case where we deal with a singel atom magnetic entity + elif type(atom) == int: + orbital_indeces = dh.geometry.a2o(atom, all=True) + if ( + type(l) == int + ): # if specified we restrict to given l angular momentum channel + orbital_indeces = orbital_indeces[ + [o.l == l for o in dh.geometry.atoms[atom].orbitals] + ] + + return orbital_indeces # numpy array containing integers labeling orbitals associated to a magnetic entity.
+ + + +
+[docs] +def blow_up_orbindx(orb_indices): + """ + Function to blow up orbital indeces to make SPIN BOX indices. """ - return x + y
+ return np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten()
+ + + +
+[docs] +def calculate_exchange_tensor(pair): + o1, o2, o3 = pair["energies"] # o1=x, o2=y, o3=z + # dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]), + # dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]), + # dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]), + + J_ii = np.array([o2[-1], o3[0], o1[0]]) # xx, yy, zz + J_S = -0.5 * np.array([o3[1] + o3[2], o2[1] + o2[1], o1[1] + o1[2]]) # yz, zx, xy + D = 0.5 * np.array([o1[1] - o1[2], o2[2] - o2[1], o3[1] - o3[2]]) # x, y, z + return J_ii.sum() / 3, D, np.concatenate([J_ii[:2] - J_ii.sum() / 3, J_S]).flatten()
+ + + + diff --git a/docs/source/_build/html/_modules/index.html b/docs/source/_build/html/_modules/index.html index 85e47e5..a7e8580 100644 --- a/docs/source/_build/html/_modules/index.html +++ b/docs/source/_build/html/_modules/index.html @@ -13,7 +13,7 @@ - + diff --git a/docs/source/_build/html/_static/basic.css b/docs/source/_build/html/_static/basic.css index 93a4776..91cf78f 100644 --- a/docs/source/_build/html/_static/basic.css +++ b/docs/source/_build/html/_static/basic.css @@ -1,5 +1,12 @@ /* + * basic.css + * ~~~~~~~~~ + * * Sphinx stylesheet -- basic theme. + * + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. + * :license: BSD, see LICENSE for details. + * */ /* -- main layout ----------------------------------------------------------- */ @@ -108,11 +115,15 @@ img { /* -- search page ----------------------------------------------------------- */ ul.search { - margin-top: 10px; + margin: 10px 0 0 20px; + padding: 0; } ul.search li { - padding: 5px 0; + padding: 5px 0 5px 20px; + background-image: url(file.png); + background-repeat: no-repeat; + background-position: 0 7px; } ul.search li a { diff --git a/docs/source/_build/html/_static/doctools.js b/docs/source/_build/html/_static/doctools.js index 0398ebb..4d67807 100644 --- a/docs/source/_build/html/_static/doctools.js +++ b/docs/source/_build/html/_static/doctools.js @@ -1,5 +1,12 @@ /* + * doctools.js + * ~~~~~~~~~~~ + * * Base JavaScript utilities for all Sphinx HTML documentation. + * + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. + * :license: BSD, see LICENSE for details. + * */ "use strict"; diff --git a/docs/source/_build/html/_static/language_data.js b/docs/source/_build/html/_static/language_data.js index a5ea78e..434cd3d 100644 --- a/docs/source/_build/html/_static/language_data.js +++ b/docs/source/_build/html/_static/language_data.js @@ -1,6 +1,13 @@ /* + * language_data.js + * ~~~~~~~~~~~~~~~~ + * * This script contains the language-specific data used by searchtools.js, * namely the list of stopwords, stemmer, scorer and splitter. + * + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. + * :license: BSD, see LICENSE for details. + * */ var stopwords = ["a", "and", "are", "as", "at", "be", "but", "by", "for", "if", "in", "into", "is", "it", "near", "no", "not", "of", "on", "or", "such", "that", "the", "their", "then", "there", "these", "they", "this", "to", "was", "will", "with"]; diff --git a/docs/source/_build/html/_static/searchtools.js b/docs/source/_build/html/_static/searchtools.js index 2c774d1..b08d58c 100644 --- a/docs/source/_build/html/_static/searchtools.js +++ b/docs/source/_build/html/_static/searchtools.js @@ -1,5 +1,12 @@ /* + * searchtools.js + * ~~~~~~~~~~~~~~~~ + * * Sphinx JavaScript utilities for the full-text search. + * + * :copyright: Copyright 2007-2024 by the Sphinx team, see AUTHORS. + * :license: BSD, see LICENSE for details. + * */ "use strict"; @@ -13,7 +20,7 @@ if (typeof Scorer === "undefined") { // and returns the new score. /* score: result => { - const [docname, title, anchor, descr, score, filename, kind] = result + const [docname, title, anchor, descr, score, filename] = result return score }, */ @@ -40,14 +47,6 @@ if (typeof Scorer === "undefined") { }; } -// Global search result kind enum, used by themes to style search results. -class SearchResultKind { - static get index() { return "index"; } - static get object() { return "object"; } - static get text() { return "text"; } - static get title() { return "title"; } -} - const _removeChildren = (element) => { while (element && element.lastChild) element.removeChild(element.lastChild); }; @@ -65,13 +64,9 @@ const _displayItem = (item, searchTerms, highlightTerms) => { const showSearchSummary = DOCUMENTATION_OPTIONS.SHOW_SEARCH_SUMMARY; const contentRoot = document.documentElement.dataset.content_root; - const [docName, title, anchor, descr, score, _filename, kind] = item; + const [docName, title, anchor, descr, score, _filename] = item; let listItem = document.createElement("li"); - // Add a class representing the item's type: - // can be used by a theme's CSS selector for styling - // See SearchResultKind for the class names. - listItem.classList.add(`kind-${kind}`); let requestUrl; let linkUrl; if (docBuilder === "dirhtml") { @@ -120,10 +115,8 @@ const _finishSearch = (resultCount) => { "Your search did not match any documents. Please make sure that all words are spelled correctly and that you've selected enough categories." ); else - Search.status.innerText = Documentation.ngettext( - "Search finished, found one page matching the search query.", - "Search finished, found ${resultCount} pages matching the search query.", - resultCount, + Search.status.innerText = _( + "Search finished, found ${resultCount} page(s) matching the search query." ).replace('${resultCount}', resultCount); }; const _displayNextItem = ( @@ -145,7 +138,7 @@ const _displayNextItem = ( else _finishSearch(resultCount); }; // Helper function used by query() to order search results. -// Each input is an array of [docname, title, anchor, descr, score, filename, kind]. +// Each input is an array of [docname, title, anchor, descr, score, filename]. // Order the results by score (in opposite order of appearance, since the // `_displayNextItem` function uses pop() to retrieve items) and then alphabetically. const _orderResultsByScoreThenName = (a, b) => { @@ -255,7 +248,6 @@ const Search = { searchSummary.classList.add("search-summary"); searchSummary.innerText = ""; const searchList = document.createElement("ul"); - searchList.setAttribute("role", "list"); searchList.classList.add("search"); const out = document.getElementById("search-results"); @@ -326,7 +318,7 @@ const Search = { const indexEntries = Search._index.indexentries; // Collect multiple result groups to be sorted separately and then ordered. - // Each is an array of [docname, title, anchor, descr, score, filename, kind]. + // Each is an array of [docname, title, anchor, descr, score, filename]. const normalResults = []; const nonMainIndexResults = []; @@ -345,7 +337,6 @@ const Search = { null, score + boost, filenames[file], - SearchResultKind.title, ]); } } @@ -363,7 +354,6 @@ const Search = { null, score, filenames[file], - SearchResultKind.index, ]; if (isMain) { normalResults.push(result); @@ -485,7 +475,6 @@ const Search = { descr, score, filenames[match[0]], - SearchResultKind.object, ]); }; Object.keys(objects).forEach((prefix) => @@ -596,7 +585,6 @@ const Search = { null, score, filenames[file], - SearchResultKind.text, ]); } return results; diff --git a/docs/source/_build/html/genindex.html b/docs/source/_build/html/genindex.html index 112bda8..585a91c 100644 --- a/docs/source/_build/html/genindex.html +++ b/docs/source/_build/html/genindex.html @@ -13,7 +13,7 @@ - + @@ -71,16 +71,35 @@

Index

- A + B + | C | G | H | M + | P + | R + | S + | T
-

A

+

B

+
+ +

C

+ + +
@@ -132,8 +151,6 @@

M

+

P

+ + + +
+ +

R

+ + + +
+ +

S

+ + +
+ +

T

+ + +
+ diff --git a/docs/source/_build/html/grogu.html b/docs/source/_build/html/grogu.html index 23c92ec..ab23e48 100644 --- a/docs/source/_build/html/grogu.html +++ b/docs/source/_build/html/grogu.html @@ -14,9 +14,8 @@ - + - @@ -50,11 +49,19 @@
  • grogu.example module
  • grogu.jij module
  • grogu.useful module
  • Module contents
    • @@ -103,62 +110,46 @@

    grogu.useful module

    -
    -grogu.useful.add(x, y)[source]
    -

    The sum of two numbers for testing.

    -

    This function adds to numbers together. I only created this for testing documentation and examples.

    -
    -
    Parameters:
    -
      -

    • x (float) – First number

      • -

    • y (float) – Second number added to x

      • -
    -
    -
    Returns:
    -

    sum – The sum of the inputs

    -
    -
    Return type:
    -

    int

    -
    -
    -
    -

    See also

    -
    -
    numpy.add

    Adds more than two numbers.

    -
    -
    -
    -

    Notes

    -

    We can create some latex notes here [1] :

    -
    -\[a + b = c\]
    -

    References

    - -

    Examples

    -
    >>> add(1, 2)
-3
-
    -
    +
    +grogu.useful.RotM(theta, u, eps=1e-10)[source]
    +

    Definition of rotation matrix with angle theta around direction u.

    -
    -grogu.useful.hsk(dh, k=(0, 0, 0))[source]
    -

    One way to speed up Hk and Sk generation

    +
    +grogu.useful.RotMa2b(a, b, eps=1e-10)[source]
    +

    Definition of rotation matrix rotating unit vector a to unit vector b. +Function returns array R such that R@a = b holds.

    -
    -grogu.useful.make_atran(nauc, dirs='xyz', dist=1)[source]
    -

    Simple pair generator. Depending on the value of the dirs -argument sampling in 1,2 or 3 dimensions is generated. -If dirs argument does not contain either of x,y or z -a single pair is returend.

    +
    +grogu.useful.blow_up_orbindx(orb_indices)[source]
    +

    Function to blow up orbital indeces to make SPIN BOX indices.

    +
    + +
    +
    +grogu.useful.calculate_exchange_tensor(pair)[source]
    +
    + +
    +
    +grogu.useful.commutator(a, b)[source]
    +

    Shorthand for commutator

    +
    + +
    +
    +grogu.useful.crossM(u)[source]
    +

    Definition for the cross-product matrix. +Acting as a cross product with vector u.

    +
    + +
    +
    +grogu.useful.hsk(H, ss, sc_off, k=(0, 0, 0))[source]
    +

    One way to speed up Hk and Sk generation

    @@ -176,6 +167,39 @@ If dirs argument does not contain either of x,y or z a kset of a single k-pont at the origin is returend.

    +
    +
    +grogu.useful.parse_magnetic_entity(dh, atom=None, l=None, **kwargs)[source]
    +

    Function to define orbital indeces of a given magnetic entity. +dh: a sisl Hamiltonian object +atom: an integer or list of integers, defining atom (or atoms) in the unicell forming the magnetic entity +l: integer, defining the angular momentum channel

    +
    + +
    +
    +grogu.useful.print_atomic_indices(pair, magnetic_entities, dh)[source]
    +
    + +
    +
    +grogu.useful.spin_tracer(M)[source]
    +

    Spin tracer utility. +This akes an operator with the orbital-spin sequence: +orbital 1 up, +orbital 1 down, +orbital 2 up, +orbital 2 down, +that is in the SPIN-BOX representation, +and extracts orbital dependent Pauli traces.

    +
    + +
    +
    +grogu.useful.tau_u(u)[source]
    +

    Pauli matrix in direction u.

    +
    +

    Module contents

    diff --git a/docs/source/_build/html/index.html b/docs/source/_build/html/index.html index 6708b9c..1c2feba 100644 --- a/docs/source/_build/html/index.html +++ b/docs/source/_build/html/index.html @@ -14,9 +14,8 @@ - + - diff --git a/docs/source/_build/html/modules.html b/docs/source/_build/html/modules.html index 68a6a30..d7f7f4a 100644 --- a/docs/source/_build/html/modules.html +++ b/docs/source/_build/html/modules.html @@ -14,9 +14,8 @@ - + - @@ -84,11 +83,19 @@
  • grogu.example module
  • grogu.jij module
  • grogu.useful module
  • Module contents
    • diff --git a/docs/source/_build/html/objects.inv b/docs/source/_build/html/objects.inv index 45e7134..9bfc70c 100644 --- a/docs/source/_build/html/objects.inv +++ b/docs/source/_build/html/objects.inv @@ -2,6 +2,6 @@ # Project: Grogu # Version: # The remainder of this file is compressed using zlib. -xڕM -@ =E@z$Ķv~3MՊv^TV~` UVVujZ`nD_ltj&B)5>ڷ3h-084sMš`kk=2tV Ձ].9}L:-VTGr^BR&0TC|٠G{!İw$>N6:QYcϑ; - \ No newline at end of file +xڕMn F>Ruvt)jxb'fPuH);0 b~H",{dPaMv#[TYMIq2_f)J|:kֺ[Watv\a*0*`8܀ё55j2TA[P*اXM\XRGXl +B +k.ʐdz,h3tM/x_뽺@RP Z@iyt?Eq7 \ No newline at end of file diff --git a/docs/source/_build/html/py-modindex.html b/docs/source/_build/html/py-modindex.html index afa884e..f2c3645 100644 --- a/docs/source/_build/html/py-modindex.html +++ b/docs/source/_build/html/py-modindex.html @@ -13,7 +13,7 @@ - + diff --git a/docs/source/_build/html/search.html b/docs/source/_build/html/search.html index 61399b5..e8b7a3c 100644 --- a/docs/source/_build/html/search.html +++ b/docs/source/_build/html/search.html @@ -14,7 +14,7 @@ - + diff --git a/docs/source/_build/html/searchindex.js b/docs/source/_build/html/searchindex.js index 3d96987..3ca980e 100644 --- a/docs/source/_build/html/searchindex.js +++ b/docs/source/_build/html/searchindex.js @@ -1 +1 @@ -Search.setIndex({"alltitles": {"Contents:": [[1, null]], "Indices and tables": [[1, "indices-and-tables"]], "Module contents": [[0, "module-grogu"]], "Submodules": [[0, "submodules"]], "asd documentation": [[1, null]], "grogu package": [[0, null]], "grogu.example module": [[0, "module-grogu.example"]], "grogu.jij module": [[0, "module-grogu.jij"]], "grogu.useful module": [[0, "module-grogu.useful"]], "src": [[2, null]]}, "docnames": ["grogu", "index", "modules"], "envversion": {"sphinx": 64, "sphinx.domains.c": 3, "sphinx.domains.changeset": 1, "sphinx.domains.citation": 1, "sphinx.domains.cpp": 9, "sphinx.domains.index": 1, "sphinx.domains.javascript": 3, "sphinx.domains.math": 2, "sphinx.domains.python": 4, "sphinx.domains.rst": 2, "sphinx.domains.std": 2, "sphinx.ext.viewcode": 1}, "filenames": ["grogu.rst", "index.rst", "modules.rst"], "indexentries": {"add() (in module grogu.useful)": [[0, "grogu.useful.add", false]], "grogu": [[0, "module-grogu", false]], "grogu.example": [[0, "module-grogu.example", false]], "grogu.jij": [[0, "module-grogu.jij", false]], "grogu.useful": [[0, "module-grogu.useful", false]], "hsk() (in module grogu.useful)": [[0, "grogu.useful.hsk", false]], "make_atran() (in module grogu.useful)": [[0, "grogu.useful.make_atran", false]], "make_contour() (in module grogu.useful)": [[0, "grogu.useful.make_contour", false]], "make_kset() (in module grogu.useful)": [[0, "grogu.useful.make_kset", false]], "module": [[0, "module-grogu", false], [0, "module-grogu.example", false], [0, "module-grogu.jij", false], [0, "module-grogu.useful", false]]}, "objects": {"": [[0, 0, 0, "-", "grogu"]], "grogu": [[0, 0, 0, "-", "example"], [0, 0, 0, "-", "jij"], [0, 0, 0, "-", "useful"]], "grogu.useful": [[0, 1, 1, "", "add"], [0, 1, 1, "", "hsk"], [0, 1, 1, "", "make_atran"], [0, 1, 1, "", "make_contour"], [0, 1, 1, "", "make_kset"]]}, "objnames": {"0": ["py", "module", "Python module"], "1": ["py", "function", "Python function"]}, "objtypes": {"0": "py:module", "1": "py:function"}, "terms": {"0": 0, "1": 0, "150": 0, "2": 0, "20": 0, "3": 0, "42": 0, "A": 0, "If": 0, "One": 0, "The": 0, "ad": 0, "add": [0, 1, 2], "argument": 0, "b": 0, "c": 0, "can": 0, "contain": 0, "content": 2, "contour": 0, "creat": 0, "depend": 0, "detail": 1, "dh": 0, "dimens": 0, "dir": 0, "dist": 0, "document": 0, "doe": 0, "either": 0, "emax": 0, "emin": 0, "en": 0, "enum": 0, "exampl": 2, "first": 0, "float": 0, "format": 0, "function": 0, "gener": 0, "grid": 0, "grogu": [1, 2], "here": 0, "hk": 0, "hsk": [0, 2], "html": 0, "http": 0, "i": 0, "index": 1, "input": 0, "int": 0, "io": 0, "jij": 2, "k": 0, "kset": 0, "latest": 0, "latex": 0, "make_atran": [0, 2], "make_contour": [0, 2], "make_kset": [0, 2], "modul": [1, 2], "more": 0, "nauc": 0, "note": 0, "number": 0, "numk": 0, "numpi": 0, "numpydoc": 0, "onli": 0, "origin": 0, "p": 0, "packag": [1, 2], "page": 1, "pair": 0, "paramet": 0, "pont": 0, "readthedoc": 0, "refer": 0, "restructuredtext": 1, "returend": 0, "return": 0, "sampl": 0, "search": 1, "second": 0, "see": 1, "simpl": 0, "singl": 0, "sk": 0, "some": 0, "sophist": 0, "sourc": 0, "speed": 0, "src": 1, "submodul": 2, "sum": 0, "syntax": 1, "test": 0, "than": 0, "thi": 0, "togeth": 0, "two": 0, "type": 0, "up": 0, "us": [1, 2], "valu": 0, "wai": 0, "we": 0, "x": 0, "xyz": 0, "y": 0, "your": 1, "z": 0}, "titles": ["grogu package", "asd documentation", "src"], "titleterms": {"asd": 1, "content": [0, 1], "document": 1, "exampl": 0, "grogu": 0, "indic": 1, "jij": 0, "modul": 0, "packag": 0, "src": 2, "submodul": 0, "tabl": 1, "us": 0}}) +Search.setIndex({"alltitles": {"Contents:": [[1, null]], "Indices and tables": [[1, "indices-and-tables"]], "Module contents": [[0, "module-grogu"]], "Submodules": [[0, "submodules"]], "asd documentation": [[1, null]], "grogu package": [[0, null]], "grogu.example module": [[0, "module-grogu.example"]], "grogu.jij module": [[0, "module-grogu.jij"]], "grogu.useful module": [[0, "module-grogu.useful"]], "src": [[2, null]]}, "docnames": ["grogu", "index", "modules"], "envversion": {"sphinx": 62, "sphinx.domains.c": 3, "sphinx.domains.changeset": 1, "sphinx.domains.citation": 1, "sphinx.domains.cpp": 9, "sphinx.domains.index": 1, "sphinx.domains.javascript": 3, "sphinx.domains.math": 2, "sphinx.domains.python": 4, "sphinx.domains.rst": 2, "sphinx.domains.std": 2, "sphinx.ext.viewcode": 1}, "filenames": ["grogu.rst", "index.rst", "modules.rst"], "indexentries": {"blow_up_orbindx() (in module grogu.useful)": [[0, "grogu.useful.blow_up_orbindx", false]], "calculate_exchange_tensor() (in module grogu.useful)": [[0, "grogu.useful.calculate_exchange_tensor", false]], "commutator() (in module grogu.useful)": [[0, "grogu.useful.commutator", false]], "crossm() (in module grogu.useful)": [[0, "grogu.useful.crossM", false]], "grogu": [[0, "module-grogu", false]], "grogu.example": [[0, "module-grogu.example", false]], "grogu.jij": [[0, "module-grogu.jij", false]], "grogu.useful": [[0, "module-grogu.useful", false]], "hsk() (in module grogu.useful)": [[0, "grogu.useful.hsk", false]], "make_contour() (in module grogu.useful)": [[0, "grogu.useful.make_contour", false]], "make_kset() (in module grogu.useful)": [[0, "grogu.useful.make_kset", false]], "module": [[0, "module-grogu", false], [0, "module-grogu.example", false], [0, "module-grogu.jij", false], [0, "module-grogu.useful", false]], "parse_magnetic_entity() (in module grogu.useful)": [[0, "grogu.useful.parse_magnetic_entity", false]], "print_atomic_indices() (in module grogu.useful)": [[0, "grogu.useful.print_atomic_indices", false]], "rotm() (in module grogu.useful)": [[0, "grogu.useful.RotM", false]], "rotma2b() (in module grogu.useful)": [[0, "grogu.useful.RotMa2b", false]], "spin_tracer() (in module grogu.useful)": [[0, "grogu.useful.spin_tracer", false]], "tau_u() (in module grogu.useful)": [[0, "grogu.useful.tau_u", false]]}, "objects": {"": [[0, 0, 0, "-", "grogu"]], "grogu": [[0, 0, 0, "-", "example"], [0, 0, 0, "-", "jij"], [0, 0, 0, "-", "useful"]], "grogu.useful": [[0, 1, 1, "", "RotM"], [0, 1, 1, "", "RotMa2b"], [0, 1, 1, "", "blow_up_orbindx"], [0, 1, 1, "", "calculate_exchange_tensor"], [0, 1, 1, "", "commutator"], [0, 1, 1, "", "crossM"], [0, 1, 1, "", "hsk"], [0, 1, 1, "", "make_contour"], [0, 1, 1, "", "make_kset"], [0, 1, 1, "", "parse_magnetic_entity"], [0, 1, 1, "", "print_atomic_indices"], [0, 1, 1, "", "spin_tracer"], [0, 1, 1, "", "tau_u"]]}, "objnames": {"0": ["py", "module", "Python module"], "1": ["py", "function", "Python function"]}, "objtypes": {"0": "py:module", "1": "py:function"}, "terms": {"0": 0, "1": 0, "10": 0, "150": 0, "1e": 0, "2": 0, "20": 0, "3": 0, "42": 0, "A": 0, "If": 0, "One": 0, "act": 0, "add": 1, "ak": 0, "an": 0, "angl": 0, "angular": 0, "argument": 0, "around": 0, "arrai": 0, "atom": 0, "b": 0, "blow": 0, "blow_up_orbindx": [0, 2], "box": 0, "calculate_exchange_tensor": [0, 2], "channel": 0, "commut": [0, 2], "contain": 0, "content": 2, "contour": 0, "cross": 0, "crossm": [0, 2], "defin": 0, "definit": 0, "depend": 0, "detail": 1, "dh": 0, "dimens": 0, "dir": 0, "direct": 0, "doe": 0, "down": 0, "either": 0, "emax": 0, "emin": 0, "entiti": 0, "enum": 0, "ep": 0, "exampl": 2, "extract": 0, "form": 0, "function": 0, "gener": 0, "given": 0, "grid": 0, "grogu": [1, 2], "h": 0, "hamiltonian": 0, "hk": 0, "hold": 0, "hsk": [0, 2], "i": 0, "indec": 0, "index": 1, "indic": 0, "integ": 0, "jij": 2, "k": 0, "kset": 0, "kwarg": 0, "l": 0, "list": 0, "m": 0, "magnet": 0, "magnetic_ent": 0, "make": 0, "make_contour": [0, 2], "make_kset": [0, 2], "matrix": 0, "modul": [1, 2], "momentum": 0, "more": 0, "none": 0, "numk": 0, "object": 0, "oper": 0, "orb_indic": 0, "orbit": 0, "origin": 0, "p": 0, "packag": [1, 2], "page": 1, "pair": 0, "parse_magnetic_ent": [0, 2], "pauli": 0, "pont": 0, "print_atomic_indic": [0, 2], "product": 0, "r": 0, "represent": 0, "restructuredtext": 1, "returend": 0, "return": 0, "rotat": 0, "rotm": [0, 2], "rotma2b": [0, 2], "sampl": 0, "sc_off": 0, "search": 1, "see": 1, "sequenc": 0, "shorthand": 0, "simpl": 0, "singl": 0, "sisl": 0, "sk": 0, "sophist": 0, "sourc": 0, "speed": 0, "spin": 0, "spin_trac": [0, 2], "src": 1, "ss": 0, "submodul": 2, "syntax": 1, "tau_u": [0, 2], "theta": 0, "thi": 0, "trace": 0, "tracer": 0, "u": 0, "unicel": 0, "unit": 0, "up": 0, "us": [1, 2], "util": 0, "valu": 0, "vector": 0, "wai": 0, "x": 0, "xyz": 0, "y": 0, "your": 1, "z": 0}, "titles": ["grogu package", "asd documentation", "src"], "titleterms": {"asd": 1, "content": [0, 1], "document": 1, "exampl": 0, "grogu": 0, "indic": 1, "jij": 0, "modul": 0, "packag": 0, "src": 2, "submodul": 0, "tabl": 1, "us": 0}}) diff --git a/local_operator.ipynb b/local_operator.ipynb index 421af10..968674a 100644 --- a/local_operator.ipynb +++ b/local_operator.ipynb @@ -220,7 +220,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.7" + "version": "3.9.6" } }, "nbformat": 4, diff --git a/pyproject.toml b/pyproject.toml index 6be180a..187cb8e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,48 +3,52 @@ requires = ["hatchling"] build-backend = "hatchling.build" [project] -name = "grogu" -version = "0.0.1" +name = "grogu_magn" +version = "0.0.3" authors = [ - { name="Laszlo Oroszlany", email="laszlo.oroszlany@ttk.elte.hu" }, { name="Daniel Pozsar", email="danielpozsar@student.elte.hu" }, + { name="Laszlo Oroszlany", email="laszlo.oroszlany@ttk.elte.hu" }, ] description = "Relativistic magnetic interactions from non-orthogonal basis sets" readme = "README.md" kewords = [ "DFT", "physics", - "grogu" + "grogu", + "magnetic interactions", ] requires-python = ">=3.9" dependencies = [ - "numpy", + "numpy==1.24.4", "scipy", - "sisl", + "sisl==0.14.3", + "netcdf4==1.6.2", "openmpi", "mpi4py", - "argparse", ] + classifiers = [ "Development Status :: 3 - Alpha", "Programming Language :: Python", "Programming Language :: Python :: 3.9", - "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3 :: Only", "Topic :: Scientific/Engineering :: Physics", "License :: OSI Approved :: MIT License", "Operating System :: OS Independent", ] +[project.scripts] +grogu = "grogu:main" + [project.urls] -Homepage = "https://github.com/pypa/sampleproject" -Documentation = "https://readthedocs.org" -Repository = "https://github.com/me/spam.git" -Issues = "https://github.com/pypa/sampleproject/issues" +Homepage = "https://gitea.vo.elte.hu/et209d/grogu" +Documentation = "https://gitea.vo.elte.hu/et209d/grogu" +Repository = "https://gitea.vo.elte.hu/et209d/grogu" +Issues = "https://gitea.vo.elte.hu/et209d/grogu" [tool.pytest.ini_options] pythonpath = [ - "src/grogu/" + "src/grogu_magn/" ] addopts = [ "--import-mode=importlib", diff --git a/requirements-dev.txt b/requirements-dev.txt index 58f3049..1d5eb40 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -14,6 +14,5 @@ pytest pytest-randomly pytest-cov isort -sphinx tqdm matplotlib diff --git a/requirements.txt b/requirements.txt index 4bd7910..8020276 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,4 @@ # Package dependencies -argparse numpy==1.24.4 scipy sisl==0.14.3 diff --git a/src/grogu/__init__.py b/src/grogu_magn/__init__.py similarity index 100% rename from src/grogu/__init__.py rename to src/grogu_magn/__init__.py diff --git a/src/grogu/example.py b/src/grogu_magn/example.py similarity index 100% rename from src/grogu/example.py rename to src/grogu_magn/example.py diff --git a/src/grogu_magn/grogu.py b/src/grogu_magn/grogu.py new file mode 100644 index 0000000..536734b --- /dev/null +++ b/src/grogu_magn/grogu.py @@ -0,0 +1,526 @@ +from useful import * + + +def main(): + import os + from sys import stdout + from timeit import default_timer as timer + + from tqdm import tqdm + + os.environ["OMP_NUM_THREADS"] = "1" # export OMP_NUM_THREADS=4 + os.environ["OPENBLAS_NUM_THREADS"] = "1" # export OPENBLAS_NUM_THREADS=4 + os.environ["MKL_NUM_THREADS"] = "1" # export MKL_NUM_THREADS=6 + os.environ["VECLIB_MAXIMUM_THREADS"] = "1" # export VECLIB_MAXIMUM_THREADS=4 + os.environ["NUMEXPR_NUM_THREADS"] = "1" # export NUMEXPR_NUM_THREADS=6 + + import warnings + + import numpy as np + import sisl + from mpi4py import MPI + from numpy.linalg import inv + + start_time = timer() + + # this cell mimicks an input file + fdf = sisl.get_sile("../../lat3_791/Fe3GeTe2.fdf") + # this information needs to be given at the input!! + scf_xcf_orientation = np.array([0, 0, 1]) # z + # list of reference directions for around which we calculate the derivatives + # o is the quantization axis, v and w are two axes perpendicular to it + # at this moment the user has to supply o,v,w on the input. + # we can have some default for this + ref_xcf_orientations = [ + dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]), + dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]), + dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]), + ] + + # human readable definition of magnetic entities + # magnetic_entities = [ + # dict(atom=0, ), + # dict(atom=1, ), + # dict(atom=2, ), + # dict(atom=3, l=2), + # dict(atom=4, l=2), + # dict(atom=5, l=2), + # ] + # pairs = [ + # dict(ai=3, aj=4, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV + # dict(ai=3, aj=5, Ruc=np.array([0, 0, 0])), # these should all be around -41.9 in the isotropic part + # dict(ai=4, aj=5, Ruc=np.array([0, 0, 0])), + # dict(ai=3, aj=0, Ruc=np.array([0, 0, 0])), + # dict(ai=3, aj=1, Ruc=np.array([0, 0, 0])), + # dict(ai=3, aj=2, Ruc=np.array([0, 0, 0])), + # ] + magnetic_entities = [ + dict(atom=3, l=2), + dict(atom=4, l=2), + dict(atom=5, l=2), + ] + + # pair information + pairs = [ + dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV + dict( + ai=0, aj=2, Ruc=np.array([0, 0, 0]) + ), # these should all be around -41.9 in the isotropic part + dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=1, Ruc=np.array([-1, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])), + dict(ai=0, aj=1, Ruc=np.array([1, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([1, 0, 0])), + dict(ai=0, aj=1, Ruc=np.array([0, -1, 0])), + dict(ai=0, aj=2, Ruc=np.array([0, -1, 0])), + dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])), + dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])), + dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])), + ] + + # Brilloun zone sampling and Green function contour integral + kset = 20 + kdirs = "xy" + ebot = -30 + eset = 50 + esetp = 1000 + + # MPI parameters + comm = MPI.COMM_WORLD + size = comm.Get_size() + rank = comm.Get_rank() + root_node = 0 + if rank == root_node: + print("Number of nodes in the parallel cluster: ", size) + + simulation_parameters = dict( + path="Not yet specified.", + scf_xcf_orientation=scf_xcf_orientation, + ref_xcf_orientations=ref_xcf_orientations, + kset=kset, + kdirs=kdirs, + ebot=ebot, + eset=eset, + esetp=esetp, + parallel_size=size, + ) + + # digestion of the input + # read in hamiltonian + dh = fdf.read_hamiltonian() + try: + simulation_parameters["geom"] = fdf.read_geometry() + except: + print("Error reading geometry.") + + # unit cell index + uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0]) + + setup_time = timer() + + NO = dh.no # shorthand for number of orbitals in the unit cell + + # preprocessing Hamiltonian and overlap matrix elements + h11 = dh.tocsr(dh.M11r) + h11 += dh.tocsr(dh.M11i) * 1.0j + h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") + + h22 = dh.tocsr(dh.M22r) + h22 += dh.tocsr(dh.M22i) * 1.0j + h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") + + h12 = dh.tocsr(dh.M12r) + h12 += dh.tocsr(dh.M12i) * 1.0j + h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") + + h21 = dh.tocsr(dh.M21r) + h21 += dh.tocsr(dh.M21i) * 1.0j + h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") + + sov = ( + dh.tocsr(dh.S_idx) + .toarray() + .reshape(NO, dh.n_s, NO) + .transpose(0, 2, 1) + .astype("complex128") + ) + + # Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation + U = np.vstack( + [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])] + ) + # This is the permutation that transforms ud1ud2 to u12d12 + # That is this transforms FROM SPIN BOX to ORBITAL BOX => U + # the inverse transformation is U.T u12d12 to ud1ud2 + # That is FROM ORBITAL BOX to SPIN BOX => U.T + + # From now on everything is in SPIN BOX!! + hh, ss = np.array( + [ + U.T + @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) + @ U + for i in range(dh.lattice.nsc.prod()) + ] + ), np.array( + [ + U.T + @ np.block( + [[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]] + ) + @ U + for i in range(dh.lattice.nsc.prod()) + ] + ) + + # symmetrizing Hamiltonian and overlap matrix to make them hermitian + for i in range(dh.lattice.sc_off.shape[0]): + j = dh.lattice.sc_index(-dh.lattice.sc_off[i]) + h1, h1d = hh[i], hh[j] + hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2 + s1, s1d = ss[i], ss[j] + ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2 + + # identifying TRS and TRB parts of the Hamiltonian + TAUY = np.kron(np.eye(NO), tau_y) + hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())]) + hTRS = (hh + hTR) / 2 + hTRB = (hh - hTR) / 2 + + # extracting the exchange field + traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77 + XCF = np.array( + [ + np.array([f["x"] for f in traced]), + np.array([f["y"] for f in traced]), + np.array([f["z"] for f in traced]), + ] + ) # equation 77 + + # Check if exchange field has scalar part + max_xcfs = abs(np.array(np.array([f["c"] for f in traced]))).max() + if max_xcfs > 1e-12: + warnings.warn( + f"Exchange field has non negligible scalar part. Largest value is {max_xcfs}" + ) + + H_and_XCF_time = timer() + + # for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions + for i, mag_ent in enumerate(magnetic_entities): + parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes + magnetic_entities[i]["orbital_indeces"] = parsed + magnetic_entities[i]["spin_box_indeces"] = blow_up_orbindx( + parsed + ) # calculate spin box indexes + spin_box_shape = len( + mag_ent["spin_box_indeces"] + ) # calculate size for Greens function generation + + mag_ent["energies"] = ( + [] + ) # we will store the second order energy derivations here + + mag_ent["Gii"] = [] # Greens function + mag_ent["Gii_tmp"] = [] # Greens function for parallelization + mag_ent["Vu1"] = [ + list([]) for _ in range(len(ref_xcf_orientations)) + ] # These will be the perturbed potentials from eq. 100 + mag_ent["Vu2"] = [list([]) for _ in range(len(ref_xcf_orientations))] + for i in ref_xcf_orientations: + mag_ent["Gii"].append( + np.zeros((eset, spin_box_shape, spin_box_shape), dtype="complex128") + ) # Greens functions for every quantization axis + mag_ent["Gii_tmp"].append( + np.zeros((eset, spin_box_shape, spin_box_shape), dtype="complex128") + ) + + # for every site we have to store 2x3 Greens function (and the associated _tmp-s) + # in the 3 reference directions, because G_ij and G_ji are both needed + for pair in pairs: + spin_box_shape_i, spin_box_shape_j = len( + magnetic_entities[pair["ai"]]["spin_box_indeces"] + ), len( + magnetic_entities[pair["aj"]]["spin_box_indeces"] + ) # calculate size for Greens function generation + + pair["energies"] = [] # we will store the second order energy derivations here + + pair["Gij"] = [] # Greens function + pair["Gji"] = [] + pair["Gij_tmp"] = [] # Greens function for parallelization + pair["Gji_tmp"] = [] + + pair["Vij"] = [ + list([]) for _ in range(len(ref_xcf_orientations)) + ] # These will be the perturbed potentials from eq. 100 + pair["Vji"] = [list([]) for _ in range(len(ref_xcf_orientations))] + + for i in ref_xcf_orientations: + pair["Gij"].append( + np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype="complex128") + ) + pair["Gij_tmp"].append( + np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype="complex128") + ) # Greens functions for every quantization axis + pair["Gji"].append( + np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128") + ) + pair["Gji_tmp"].append( + np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128") + ) + + site_and_pair_dictionaries_time = timer() + + kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling + wkset = np.ones(len(kset)) / len(kset) # generate weights for k points + kpcs = np.array_split(kset, size) # split the k points based on MPI size + kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop", file=stdout) + + k_set_time = timer() + + # this will contain all the data needed to calculate the energy variations upon rotation + hamiltonians = [] + + # iterate over the reference directions (quantization axes) + for i, orient in enumerate(ref_xcf_orientations): + # obtain rotated exchange field + R = RotMa2b(scf_xcf_orientation, orient["o"]) + rot_XCF = np.einsum("ij,jklm->iklm", R, XCF) + rot_H_XCF = sum( + [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])] + ) + rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx] + + # obtain total Hamiltonian with the rotated exchange field + rot_H = hTRS + rot_H_XCF # equation 76 + + hamiltonians.append( + dict(orient=orient["o"], H=rot_H, rotations=[]) + ) # store orientation and rotated Hamiltonian + + for u in orient[ + "vw" + ]: # these are the infinitezimal rotations (for now) perpendicular to the quantization axis + Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H + + Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100 + Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100 + + for mag_ent in magnetic_entities: + mag_ent["Vu1"][i].append( + Vu1[:, mag_ent["spin_box_indeces"]][mag_ent["spin_box_indeces"], :] + ) # fill up the perturbed potentials (for now) based on the on-site projections + mag_ent["Vu2"][i].append( + Vu2[:, mag_ent["spin_box_indeces"]][mag_ent["spin_box_indeces"], :] + ) + + for pair in pairs: + ai = magnetic_entities[pair["ai"]][ + "spin_box_indeces" + ] # get the pair orbital sizes from the magnetic entities + aj = magnetic_entities[pair["aj"]]["spin_box_indeces"] + pair["Vij"][i].append( + Vu1[:, ai][aj, :] + ) # fill up the perturbed potentials (for now) based on the on-site projections + pair["Vji"][i].append(Vu1[:, aj][ai, :]) + + reference_rotations_time = timer() + + if rank == root_node: + print("Number of magnetic entities being calculated: ", len(magnetic_entities)) + print( + "We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site." + ) + print(f"The shape of the Hamiltonian and the Greens function is {NO}x{NO}.") + comm.Barrier() + # ---------------------------------------------------------------------- + + # make energy contour + # we are working in eV now ! + # and sisil shifts E_F to 0 ! + cont = make_contour(emin=ebot, enum=eset, p=esetp) + eran = cont.ze + + # ---------------------------------------------------------------------- + # sampling the integrand on the contour and the BZ + for k in kpcs[rank]: + wk = wkset[rank] # weight of k point in BZ integral + for i, hamiltonian_orientation in enumerate( + hamiltonians + ): # iterate over reference directions + # calculate Greens function + H = hamiltonian_orientation["H"] + HK, SK = hsk(H, ss, dh.sc_off, k) + Gk = inv(SK * eran.reshape(eset, 1, 1) - HK) + + # store the Greens function slice of the magnetic entities (for now) based on the on-site projections + for mag_ent in magnetic_entities: + mag_ent["Gii_tmp"][i] += ( + Gk[:, mag_ent["spin_box_indeces"]][..., mag_ent["spin_box_indeces"]] + * wk + ) + + for pair in pairs: + # add phase shift based on the cell difference + phase = np.exp(1j * 2 * np.pi * k @ pair["Ruc"].T) + + # get the pair orbital sizes from the magnetic entities + ai = magnetic_entities[pair["ai"]]["spin_box_indeces"] + aj = magnetic_entities[pair["aj"]]["spin_box_indeces"] + + # store the Greens function slice of the magnetic entities (for now) based on the on-site projections + pair["Gij_tmp"][i] += Gk[:, ai][..., aj] * phase * wk + pair["Gji_tmp"][i] += Gk[:, aj][..., ai] * phase * wk + + # summ reduce partial results of mpi nodes + for i in range(len(hamiltonians)): + for mag_ent in magnetic_entities: + comm.Reduce(mag_ent["Gii_tmp"][i], mag_ent["Gii"][i], root=root_node) + + for pair in pairs: + comm.Reduce(pair["Gij_tmp"][i], pair["Gij"][i], root=root_node) + comm.Reduce(pair["Gji_tmp"][i], pair["Gji"][i], root=root_node) + + green_function_inversion_time = timer() + + if rank == root_node: + # iterate over the magnetic entities + for tracker, mag_ent in enumerate(magnetic_entities): + # iterate over the quantization axes + for i, Gii in enumerate(mag_ent["Gii"]): + storage = [] + # iterate over the first and second order local perturbations + for Vu1, Vu2 in zip(mag_ent["Vu1"][i], mag_ent["Vu2"][i]): + # The Szunyogh-Lichtenstein formula + traced = np.trace( + (Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2 + ) + # evaluation of the contour integral + storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we))) + + # fill up the magnetic entities dictionary with the energies + mag_ent["energies"].append(storage) + + # iterate over the pairs + for tracker, pair in enumerate(pairs): + # iterate over the quantization axes + for i, (Gij, Gji) in enumerate(zip(pair["Gij"], pair["Gji"])): + site_i = magnetic_entities[pair["ai"]] + site_j = magnetic_entities[pair["aj"]] + + storage = [] + # iterate over the first order local perturbations in all possible orientations for the two sites + for Vui in site_i["Vu1"][i]: + for Vuj in site_j["Vu1"][i]: + # The Szunyogh-Lichtenstein formula + traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2) + # evaluation of the contour integral + storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we))) + + # fill up the pairs dictionary with the energies + pairs[tracker]["energies"].append(storage) + + end_time = timer() + + print( + "############################### GROGU OUTPUT ###################################" + ) + print( + "================================================================================" + ) + print("Input file: ") + print(simulation_parameters["path"]) + print( + "Number of nodes in the parallel cluster: ", + simulation_parameters["parallel_size"], + ) + print( + "================================================================================" + ) + try: + print("Cell [Ang]: ") + print(simulation_parameters["geom"].cell) + except: + print("Geometry could not be read.") + print( + "================================================================================" + ) + print("DFT axis: ") + print(simulation_parameters["scf_xcf_orientation"]) + print("Quantization axis and perpendicular rotation directions:") + for ref in ref_xcf_orientations: + print(ref["o"], " --» ", ref["vw"]) + print( + "================================================================================" + ) + print("number of k points: ", simulation_parameters["kset"]) + print("k point directions: ", simulation_parameters["kdirs"]) + print( + "================================================================================" + ) + print("Parameters for the contour integral:") + print("Ebot: ", simulation_parameters["ebot"]) + print("Eset: ", simulation_parameters["eset"]) + print("Esetp: ", simulation_parameters["esetp"]) + print( + "================================================================================" + ) + print("Atomic informations: ") + print("") + print("") + print("Not yet specified.") + print("") + print("") + print( + "================================================================================" + ) + print("Exchange [meV]") + print( + "--------------------------------------------------------------------------------" + ) + print("Atom1 Atom2 [i j k] d [Ang]") + print( + "--------------------------------------------------------------------------------" + ) + for pair in pairs: + J_iso, J_S, D = calculate_exchange_tensor(pair) + J_iso = J_iso * sisl.unit_convert("eV", "meV") + J_S = J_S * sisl.unit_convert("eV", "meV") + D = D * sisl.unit_convert("eV", "meV") + + print(print_atomic_indices(pair, magnetic_entities, dh)) + print("Isotropic: ", J_iso) + print("DMI: ", D) + print("Symmetric-anisotropy: ", J_S) + print("") + + print( + "================================================================================" + ) + print("Runtime information: ") + print("Total runtime: ", end_time - start_time) + print( + "--------------------------------------------------------------------------------" + ) + print("Initial setup: ", setup_time - start_time) + print( + f"Hamiltonian conversion and XC field extraction: {H_and_XCF_time - setup_time:.3f} s" + ) + print( + f"Pair and site datastructure creatrions: {site_and_pair_dictionaries_time - H_and_XCF_time:.3f} s" + ) + print( + f"k set cration and distribution: {k_set_time - site_and_pair_dictionaries_time:.3f} s" + ) + print(f"Rotating XC potential: {reference_rotations_time - k_set_time:.3f} s") + print( + f"Greens function inversion: {green_function_inversion_time - reference_rotations_time:.3f} s" + ) + print( + f"Calculate energies and magnetic components: {end_time - green_function_inversion_time:.3f} s" + ) + + +if __name__ == "__main__": + main() diff --git a/src/grogu/jij.py b/src/grogu_magn/jij.py similarity index 99% rename from src/grogu/jij.py rename to src/grogu_magn/jij.py index 3de5d32..c91c289 100644 --- a/src/grogu/jij.py +++ b/src/grogu_magn/jij.py @@ -31,7 +31,7 @@ if __name__ == "__main__": from mpi4py import MPI from scipy.special import roots_legendre - from grogu.useful import hsk, make_atran, make_contour, make_kset + from grogu_magn.useful import hsk, make_atran, make_contour, make_kset start = timer() diff --git a/src/grogu/useful.py b/src/grogu_magn/useful.py similarity index 100% rename from src/grogu/useful.py rename to src/grogu_magn/useful.py diff --git a/test.ipynb b/test.ipynb index 7b1f469..07c01cd 100644 --- a/test.ipynb +++ b/test.ipynb @@ -9,7 +9,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "[Daniels-Air:17787] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Daniels-Air.501/jf.0/2288320512/sm_segment.Daniels-Air.501.88650000.0 could be created.\n" + "[Daniels-Air:31350] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Daniels-Air.501/jf.0/1950875648/sm_segment.Daniels-Air.501.74480000.0 could be created.\n" ] } ], @@ -27,7 +27,7 @@ "\n", "import numpy as np\n", "import sisl\n", - "from grogu.useful import *\n", + "from grogu_magn.useful import *\n", "from mpi4py import MPI\n", "from numpy.linalg import inv\n", "import warnings\n", @@ -375,20 +375,7 @@ "Number of magnetic entities being calculated: 3\n", "We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.\n", "The shape of the Hamiltonian and the Greens function is 84x84.\n", - "k loop: 27%|██▋ | 107/400 [00:29<01:20, 3.62it/s]\n" - ] - }, - { - "ename": "KeyboardInterrupt", - "evalue": "", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[7], line 26\u001b[0m\n\u001b[1;32m 24\u001b[0m H \u001b[38;5;241m=\u001b[39m hamiltonian_orientation[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mH\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[1;32m 25\u001b[0m HK, SK \u001b[38;5;241m=\u001b[39m hsk(H, ss, dh\u001b[38;5;241m.\u001b[39msc_off, k)\n\u001b[0;32m---> 26\u001b[0m Gk \u001b[38;5;241m=\u001b[39m \u001b[43minv\u001b[49m\u001b[43m(\u001b[49m\u001b[43mSK\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43meran\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mreshape\u001b[49m\u001b[43m(\u001b[49m\u001b[43meset\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mHK\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 28\u001b[0m \u001b[38;5;66;03m# store the Greens function slice of the magnetic entities (for now) based on the on-site projections\u001b[39;00m\n\u001b[1;32m 29\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m mag_ent \u001b[38;5;129;01min\u001b[39;00m magnetic_entities:\n", - "File \u001b[0;32m<__array_function__ internals>:200\u001b[0m, in \u001b[0;36minv\u001b[0;34m(*args, **kwargs)\u001b[0m\n", - "File \u001b[0;32m~/Downloads/grogu/.venv/lib/python3.9/site-packages/numpy/linalg/linalg.py:538\u001b[0m, in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m 536\u001b[0m signature \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mD->D\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m isComplexType(t) \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124md->d\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 537\u001b[0m extobj \u001b[38;5;241m=\u001b[39m get_linalg_error_extobj(_raise_linalgerror_singular)\n\u001b[0;32m--> 538\u001b[0m ainv \u001b[38;5;241m=\u001b[39m \u001b[43m_umath_linalg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43minv\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msignature\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msignature\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mextobj\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mextobj\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 539\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m wrap(ainv\u001b[38;5;241m.\u001b[39mastype(result_t, copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m))\n", - "\u001b[0;31mKeyboardInterrupt\u001b[0m: " + "k loop: 100%|██████████| 400/400 [01:46<00:00, 3.75it/s]\n" ] } ], @@ -453,9 +440,75 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "############################### GROGU OUTPUT ###################################\n", + "================================================================================\n", + "Input file: \n", + "Not yet specified.\n", + "Number of nodes in the parallel cluster: 1\n", + "================================================================================\n", + "Cell [Ang]: \n", + "[[ 3.79100000e+00 0.00000000e+00 0.00000000e+00]\n", + " [-1.89550000e+00 3.28310231e+00 0.00000000e+00]\n", + " [ 1.25954923e-15 2.18160327e-15 2.05700000e+01]]\n", + "================================================================================\n", + "DFT axis: \n", + "[0 0 1]\n", + "Quantization axis and perpendicular rotation directions:\n", + "[1 0 0] --» [array([0, 1, 0]), array([0, 0, 1])]\n", + "[0 1 0] --» [array([1, 0, 0]), array([0, 0, 1])]\n", + "[0 0 1] --» [array([1, 0, 0]), array([0, 1, 0])]\n", + "================================================================================\n", + "number of k points: 20\n", + "k point directions: xy\n", + "================================================================================\n", + "Parameters for the contour integral:\n", + "Ebot: -30\n", + "Eset: 50\n", + "Esetp: 10000\n", + "================================================================================\n", + "Atomic informations: \n", + "\n", + "\n", + "Not yet specified.\n", + "\n", + "\n", + "================================================================================\n", + "Exchange [meV]\n", + "--------------------------------------------------------------------------------\n", + "Atom1 Atom2 [i j k] d [Ang]\n", + "--------------------------------------------------------------------------------\n", + "[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n", + "Isotropic: -60.89330922309355\n", + "DMI: [-5.95741551e+00 7.27737654e+00 6.90431275e-04 -8.11057566e-04\n", + " -5.49031203e-06]\n", + "Symmetric-anisotropy: [-9.32966923e-01 -8.92579299e-04 -2.04258659e-06]\n", + "\n", + "[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n", + "Isotropic: -60.55651225519789\n", + "DMI: [-0.34565796 0.65919757 0.07106945 -6.23149871 -0.0424978 ]\n", + "Symmetric-anisotropy: [ 3.78506176e+00 -6.13838308e+00 3.59037036e-03]\n", + "\n", + "================================================================================\n", + "Runtime information: \n", + "Total runtime: 107.576339458\n", + "--------------------------------------------------------------------------------\n", + "Initial setup: 0.18185883300000016\n", + "Hamiltonian conversion and XC field extraction: 0.600 s\n", + "Pair and site datastructure creatrions: 0.010 s\n", + "k set cration and distribution: 0.016 s\n", + "Rotating XC potential: 0.230 s\n", + "Greens function inversion: 106.512 s\n", + "Calculate energies and magnetic components: 0.027 s\n" + ] + } + ], "source": [ "if rank == root_node:\n", " # iterate over the magnetic entities\n", @@ -565,7 +618,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -575,16 +628,27 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "dh.geometry.plot(axes=\"xy\")" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -593,9 +657,27 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "xyz[-3:]: red, green, blue\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABNoAAAHACAYAAAB0/gUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCsElEQVR4nO3df5hWdZ0//ucMyAwWM2rK8GMnf5b4E0iTMFu1Rln1YmU/5aLtAktqW9l+VLYfUiqWJqlpbEayaUbq9vG35ie9IGUjM1lNlNZMTRSElEH9qjOICjJzf//g47QToAOcmXtmeDyu61x5n/u87/v1PkznNfdzzn1ORalUKgUAAAAA2CqV5S4AAAAAAHoDQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEAB+pa7gO6otbU1zz//fAYMGJCKiopylwPQ45VKpaxatSpDhgxJZaW/8egzAMXSZzak1wAUq6O9RtC2Ec8//3zq6+vLXQZAr7N8+fL81V/9VbnLKDt9BqBz6DN/ptcAdI536zWCto0YMGBAkvU7r6ampszVAPR8zc3Nqa+vbzu+buv0GYBi6TMb0msAitXRXiNo24i3T62uqanRlAAK5Ksr6+kzAJ1Dn/kzvQagc7xbr3EBAwAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAAChA33IX0KusWJHMnp388Y/JgAHJCSckhx2WVFSUuzIAeoM330xuvDG59971veXww5NPfSqpri53ZQAA0C2tW5fccUcyd+76//7wh5N/+If1sU1nELQV5fLLkzPPTEqlpLLyz+uOOCK5/faktrac1QHQ0/32t8lxxyUvvpj0/X/t+6qrki99KbnrruRDHypvfQAA0M0sXpyMGZM888yff4X+8Y+TL385ueWW5Oiji39PXx0twq23Jv/7fyctLUlr6/qIdN269c/9+tfJ+PHlrQ+Anm3FiuSoo5KXX17/+H/2mZdeShoakhdeKF99AADQzbz+evLxjyfPPrv+8du/QpdK658bOzb5wx+Kf19B29YqlZJvfGPTXw9taVl/fuKiRV1aFgC9yKxZyWuvre8pf6mlJWlqWn92GwAAkCS54YZk+fKN/wrd2rp++bd/K/59BW1b67nnkv/+7/WB26b07bv+66MAsCVuumnjvyG8rbV1/TYAAECS9V8+rHyH1GvduvWXPy6aoG1rvf76u29TUZG88Ubn1wJA77R6dTHbAADANmL16vV/j34nb75Z/PsK2rZWfX3ynve88zZvvZXsv3/X1ANA7zNy5J+v3roxffuu3wYAAEiSDB/+zr9CV1Z2TlQjaNta/fsnJ5+c9Omz8ecrKpIddkhOOKFLywKgF/nCF/5884ONWbcu+fznu64eAADo5v75n9/5V+jW1uSLXyz+fQVtRfjmN5NhwzYM2/r0Wb/8x38k1dXlqQ2Anu+oo5LTTlv/3//z5jtvX3RiypTkiCO6vCwAAOiuhg1LLr54/X//z2u1VVSsX/7X/0r+8R+Lf19BWxFqa5Pf/CY566xkp53Wr6usXH+v2PvvT449trz1AdCzVVQkl1+e/PjHyb77/nn9/vsn116bfOc75asNAAC6qS9/OfnZz5KPfOTP63bfff3dRm+8cdNfTtwaFaXSO90uc9vU3Nyc2traNDU1paamZvMGt7YmTU3J9tsnVVWdUyBAD7NVx9VeaKv3R3Pz+vBtwIDiiwPogfSZDdknAO2tXr3+q6Q1Ne2/JNJRHT2uvsNl4dgilZXJjjuWuwoAejMfmAAAYLO8230si+KrowAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFCAsgZt9957b8aOHZshQ4akoqIit99+e7vnb7311hx99NF53/vel4qKiixatOhdX3P27NmpqKhot1RXV3fOBADo1vQZAACgK5U1aFu9enWGDx+emTNnbvL5ww47LBdddNFmvW5NTU1WrFjRtjz77LNFlAtAD6PPAAAAXalvOd/8mGOOyTHHHLPJ5ydMmJAkWbp06Wa9bkVFRQYNGrQ1pQHQC+gzAABAV+qV12h77bXXsuuuu6a+vj7HH398HnvssXfcfs2aNWlubm63AMCm6DMAAMDG9Lqgbe+9987VV1+dn/3sZ7nuuuvS2tqaQw89NH/60582OWb69Ompra1tW+rr67uwYgB6En0GAADYlF4XtI0ePToTJ07MiBEjcvjhh+fWW2/NLrvskn//93/f5JipU6emqampbVm+fHkXVgxAT6LPAAAAm1LWa7R1he222y4jR47M4sWLN7lNVVVVqqqqurAqAHoLfQYAAHhbrzuj7S+1tLTk0UcfzeDBg8tdCgC9kD4DAAC8raxntL322mvtzgBYsmRJFi1alJ122invf//78/LLL2fZsmV5/vnnkyRPPvlkkmTQoEFtd3ubOHFihg4dmunTpydJvvnNb+YjH/lI9tprr7z66qu55JJL8uyzz+aUU07p4tkBUG76DAAA0JXKGrQ99NBDOfLII9seT5kyJUkyadKkzJ49O3fccUcmT57c9vyJJ56YJJk2bVrOO++8JMmyZctSWfnnE/NeeeWVnHrqqWlsbMyOO+6Ygw46KPfff3/23XffLpgRAN2JPgMAAHSlilKpVCp3Ed1Nc3Nzamtr09TUlJqamnKXA9DjOa62Z38AFKu7H1fvvffeXHLJJVm4cGFWrFiR2267LePGjXvHMfPnz8+UKVPy2GOPpb6+PmeffXb+6Z/+qcPv2d33CUBP09Hjaq+/RhsAAEA5rV69OsOHD8/MmTM7tP2SJUty3HHH5cgjj8yiRYtyxhln5JRTTsncuXM7uVIAtlavv+soAABAOR1zzDE55phjOrz9rFmzsvvuu+fSSy9Nkuyzzz6577778t3vfjdjxozprDIBKIAz2gAAALqRBQsWpKGhod26MWPGZMGCBZscs2bNmjQ3N7dbAOh6gjYAAIBupLGxMXV1de3W1dXVpbm5OW+88cZGx0yfPj21tbVtS319fVeUCsBfELQBAAD0cFOnTk1TU1Pbsnz58nKXBLBNco02AACAbmTQoEFZuXJlu3UrV65MTU1N+vfvv9ExVVVVqaqq6oryAHgHzmgDAADoRkaPHp158+a1W3f33Xdn9OjRZaoIgI4StAEAAHSi1157LYsWLcqiRYuSJEuWLMmiRYuybNmyJOu/9jlx4sS27T/3uc/lmWeeyVe+8pU88cQT+cEPfpAbb7wxZ555ZjnKB2AzCNoAAAA60UMPPZSRI0dm5MiRSZIpU6Zk5MiROffcc5MkK1asaAvdkmT33XfPnXfembvvvjvDhw/PpZdemquuuipjxowpS/0AdJxrtAEAAHSiI444IqVSaZPPz549e6NjHnnkkU6sCoDO4Iw2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAApQ1aLv33nszduzYDBkyJBUVFbn99tvbPX/rrbfm6KOPzvve975UVFRk0aJFHXrdm266KcOGDUt1dXUOOOCA3HXXXcUXD0C3p88AAABdqaxB2+rVqzN8+PDMnDlzk88fdthhueiiizr8mvfff39OOumknHzyyXnkkUcybty4jBs3Lr///e+LKhuAHkKfAQAAulJFqVQqlbuIJKmoqMhtt92WcePGbfDc0qVLs/vuu+eRRx7JiBEj3vF1xo8fn9WrV+fnP/9527qPfOQjGTFiRGbNmtWhWpqbm1NbW5umpqbU1NRszjQA2IjucFzVZwB6L8fVDdknAMXq6HG1112jbcGCBWloaGi3bsyYMVmwYMEmx6xZsybNzc3tFgDYGH0GAADYlF4XtDU2Nqaurq7durq6ujQ2Nm5yzPTp01NbW9u21NfXd3aZAPRQ+gwAALApvS5o2xJTp05NU1NT27J8+fJylwRAL6LPAADAtqFvuQso2qBBg7Jy5cp261auXJlBgwZtckxVVVWqqqo6uzQAegF9BgAA2JRed0bb6NGjM2/evHbr7r777owePbpMFQHQm+gzAADAppT1jLbXXnstixcvbnu8ZMmSLFq0KDvttFPe//735+WXX86yZcvy/PPPJ0mefPLJJOvPJnj7zIGJEydm6NChmT59epLk9NNPz+GHH55LL700xx13XK6//vo89NBD+eEPf9jFswOg3PQZAACgK5X1jLaHHnooI0eOzMiRI5MkU6ZMyciRI3PuuecmSe64446MHDkyxx13XJLkxBNPzMiRIzNr1qy211i2bFlWrFjR9vjQQw/NT3/60/zwhz/M8OHDc/PNN+f222/P/vvv34UzA6A70GcAAICuVFEqlUrlLqK7aW5uTm1tbZqamlJTU1PucgB6PMfV9uwPgGL1hOPqzJkzc8kll6SxsTHDhw/P5ZdfnkMOOWST28+YMSNXXHFFli1blp133jmf+tSnMn369FRXV3fo/XrCPgHoSTp6XO1112gDAADoTm644YZMmTIl06ZNy8MPP5zhw4dnzJgxeeGFFza6/U9/+tOcddZZmTZtWh5//PH86Ec/yg033JCvfe1rXVw5AJtL0AYAANCJLrvsspx66qmZPHly9t1338yaNSvbb799rr766o1uf//99+ejH/1oPv3pT2e33XbL0UcfnZNOOikPPvhgF1cOwOYStAEAAHSStWvXZuHChWloaGhbV1lZmYaGhixYsGCjYw499NAsXLiwLVh75plnctddd+XYY4/d5PusWbMmzc3N7RYAul5Z7zoKAADQm7300ktpaWlJXV1du/V1dXV54oknNjrm05/+dF566aUcdthhKZVKWbduXT73uc+941dHp0+fnm984xuF1g7A5nNGGwAAQDcyf/78XHjhhfnBD36Qhx9+OLfeemvuvPPOnH/++ZscM3Xq1DQ1NbUty5cv78KKAXibM9oAAAA6yc4775w+ffpk5cqV7davXLkygwYN2uiYc845JxMmTMgpp5ySJDnggAOyevXqfPazn83Xv/71VFZueL5EVVVVqqqqip8AAJvFGW0AAACdpF+/fjnooIMyb968tnWtra2ZN29eRo8evdExr7/++gZhWp8+fZIkpVKp84oFYKs5ow0AAKATTZkyJZMmTcrBBx+cQw45JDNmzMjq1aszefLkJMnEiRMzdOjQTJ8+PUkyduzYXHbZZRk5cmRGjRqVxYsX55xzzsnYsWPbAjcAuidBGwAAQCcaP358XnzxxZx77rlpbGzMiBEjMmfOnLYbJCxbtqzdGWxnn312KioqcvbZZ+e5557LLrvskrFjx+Zb3/pWuaYAQAdVlJx7vIHm5ubU1tamqakpNTU15S4HoMdzXG3P/gAoluPqhuwTgGJ19LjqGm0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFKGvQdu+992bs2LEZMmRIKioqcvvtt7d7vlQq5dxzz83gwYPTv3//NDQ05KmnnnrH1zzvvPNSUVHRbhk2bFgnzgKA7kyvAQAAukpZg7bVq1dn+PDhmTlz5kafv/jii/O9730vs2bNygMPPJD3vOc9GTNmTN588813fN399tsvK1asaFvuu+++zigfgB5ArwEAALpK33K++THHHJNjjjlmo8+VSqXMmDEjZ599do4//vgkyTXXXJO6urrcfvvtOfHEEzf5un379s2gQYM6pWYAeha9BgAA6Crd9hptS5YsSWNjYxoaGtrW1dbWZtSoUVmwYME7jn3qqacyZMiQ7LHHHvmHf/iHLFu27B23X7NmTZqbm9stAPR+XdVr9BkAANg2dNugrbGxMUlSV1fXbn1dXV3bcxszatSozJ49O3PmzMkVV1yRJUuW5GMf+1hWrVq1yTHTp09PbW1t21JfX1/MJADo1rqq1+gzAACwbei2QduWOuaYY3LCCSfkwAMPzJgxY3LXXXfl1VdfzY033rjJMVOnTk1TU1Pbsnz58i6sGICeZnN7jT4DAADbhm4btL193ZuVK1e2W79y5crNuibODjvskA9+8INZvHjxJrepqqpKTU1NuwWA3q+reo0+AwAA24ZuG7TtvvvuGTRoUObNm9e2rrm5OQ888EBGjx7d4dd57bXX8vTTT2fw4MGdUSYAPZheAwAAFKmsQdtrr72WRYsWZdGiRUnWX5R60aJFWbZsWSoqKnLGGWfkggsuyB133JFHH300EydOzJAhQzJu3Li21/jEJz6R73//+22Pv/SlL+VXv/pVli5dmvvvvz9/93d/lz59+uSkk07q4tkB0B3oNQAAQFfpW843f+ihh3LkkUe2PZ4yZUqSZNKkSZk9e3a+8pWvZPXq1fnsZz+bV199NYcddljmzJmT6urqtjFPP/10XnrppbbHf/rTn3LSSSfl//v//r/ssssuOeyww/Jf//Vf2WWXXbpuYgB0G3oNAADQVSpKpVKp3EV0N83NzamtrU1TU5Pr6AAUwHG1PfsDoFg94bg6c+bMXHLJJWlsbMzw4cNz+eWX55BDDtnk9q+++mq+/vWv59Zbb83LL7+cXXfdNTNmzMixxx7boffrCfsEoCfp6HG1rGe0AQAA9HY33HBDpkyZklmzZmXUqFGZMWNGxowZkyeffDIDBw7cYPu1a9fmqKOOysCBA3PzzTdn6NChefbZZ7PDDjt0ffEAbBZBGwAAQCe67LLLcuqpp2by5MlJklmzZuXOO+/M1VdfnbPOOmuD7a+++uq8/PLLuf/++7PddtslSXbbbbeuLBmALdRt7zoKAADQ061duzYLFy5MQ0ND27rKyso0NDRkwYIFGx1zxx13ZPTo0TnttNNSV1eX/fffPxdeeGFaWlo2+T5r1qxJc3NzuwWAridoAwAA6CQvvfRSWlpaUldX1259XV1dGhsbNzrmmWeeyc0335yWlpbcddddOeecc3LppZfmggsu2OT7TJ8+PbW1tW1LfX19ofMAoGMEbQAAAN1Ia2trBg4cmB/+8Ic56KCDMn78+Hz961/PrFmzNjlm6tSpaWpqaluWL1/ehRUD8DbXaAMAAOgkO++8c/r06ZOVK1e2W79y5coMGjRoo2MGDx6c7bbbLn369Glbt88++6SxsTFr165Nv379NhhTVVWVqqqqYosHYLM5ow0AAKCT9OvXLwcddFDmzZvXtq61tTXz5s3L6NGjNzrmox/9aBYvXpzW1ta2dX/84x8zePDgjYZsAHQfgjYAAIBONGXKlFx55ZX5yU9+kscffzyf//zns3r16ra7kE6cODFTp05t2/7zn/98Xn755Zx++un54x//mDvvvDMXXnhhTjvttHJNAYAO8tVRAACATjR+/Pi8+OKLOffcc9PY2JgRI0Zkzpw5bTdIWLZsWSor/3wORH19febOnZszzzwzBx54YIYOHZrTTz89X/3qV8s1BQA6qKJUKpXKXUR309zcnNra2jQ1NaWmpqbc5QD0eI6r7dkfAMVyXN2QfQJQrI4eV311FAAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACtC33AUA0L09+WTyk58kzz+fDBqUTJiQ7LdfuasCoLd45ZXkmmuS3/0uqapKxo5NxoxJ+vQpd2UAsPkEbQBsVGtrcvrpyfe/n/Ttm5RKSUVFctFFycknJ7NmrV8PAFvqlluSf/zHZM2aPwdrs2Yl+++fzJmTDB1a3voAYHP56igAG/Wtb60P2ZJk3bqkpWX9/ybJ1VcnX/96+WoDoOd74IFk/Pj1IVuptL7HvN1nnngiOfroPz8GgJ5ii4K2j3/84/nGN76xwfpXXnklH//4x7e6KADK6/XXk0su2fTzpVLyve8lTU2d8/76DEDv9+1vrz9TulTa8Ll165I//CG5666ur+ttkyZNyr333lu+AgDokbYoaJs/f36+//3vZ9y4cVm9enXb+rVr1+ZXv/pVYcUBUB733pusWvXO27z5ZnL33Z3z/voMQO/W0pL83//7zmes9e2b3HZb19X0l5qamtLQ0JAPfOADufDCC/Pcc8+VrxgAeowt/uroPffck8bGxnzkIx/J0qVLCywJgHJ7/fWObffGG51Xgz4D0Hu9fUmCd9La2rl95t3cfvvtee655/L5z38+N9xwQ3bbbbccc8wxufnmm/PWW2+VrzAAurUtDtoGDx6cX/3qVznggAPy4Q9/OPPnzy+wLADKaf/9i91uS+gzAL1XVVWy++7rvzr6Tjqzz3TELrvskilTpuR3v/tdHnjggey1116ZMGFChgwZkjPPPDNPPfVUeQsEoNvZoqCt4v91xKqqqvz0pz/N6aefnr/5m7/JD37wg0KLA6A8PvjB5PDD/3wHuL/Up09y0EHJyJGd8/76DEDv9y//8s7PV1Ymn/lM19TyblasWJG77747d999d/r06ZNjjz02jz76aPbdd99897vfLXd5AHQjfbdkUOkvrlh69tlnZ5999smkSZMKKQqA8rvqqmT06OTVV9tfQ6dv3+S9701+8pPOe299BqD3O+205Oc/T+bPX/810bf16bP+a6VXXJEMGVK28vLWW2/ljjvuyI9//OP84he/yIEHHpgzzjgjn/70p1NTU5Mkue222/KZz3wmZ555ZvkKBaBb2aKgbcmSJdlll13arfvkJz+ZYcOG5aGHHiqkMADKa6+9kocfTqZPT2bPXn+dnKqq5B//Mfna15I99ui899ZnAHq/fv3W31X03/4tufzy5E9/Wr/+r/86mTo1Oeqo8tY3ePDgtLa25qSTTsqDDz6YESNGbLDNkUcemR122KHLawOg+6oo/eVpA6S5uTm1tbVpampq+2sVwLZs3bqkuTkZMCDZbrvNH++42p79AdBeqZQ0Na0P37bffvPHd8Zx9dprr80JJ5yQ6urqQl6vq+k1AMXq6HF1i85oA2Db0rdvstNO5a4CgN6qoiLpbieGTZgwodwlANADbfFdRwEAAACAPxO0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABShr0Hbvvfdm7NixGTJkSCoqKnL77be3e75UKuXcc8/N4MGD079//zQ0NOSpp55619edOXNmdtttt1RXV2fUqFF58MEHO2kGAHR3eg0AANBVyhq0rV69OsOHD8/MmTM3+vzFF1+c733ve5k1a1YeeOCBvOc978mYMWPy5ptvbvI1b7jhhkyZMiXTpk3Lww8/nOHDh2fMmDF54YUXOmsaAHRjeg0AANBVKkqlUqncRSRJRUVFbrvttowbNy7J+jMMhgwZkn/913/Nl770pSRJU1NT6urqMnv27Jx44okbfZ1Ro0blwx/+cL7//e8nSVpbW1NfX59/+Zd/yVlnndWhWpqbm1NbW5umpqbU1NRs/eQAtnHd5bjaXXpNd9kfAL2F4+qG7BOAYnX0uNptr9G2ZMmSNDY2pqGhoW1dbW1tRo0alQULFmx0zNq1a7Nw4cJ2YyorK9PQ0LDJMQBsu/QaAACgSH3LXcCmNDY2Jknq6urara+rq2t77i+99NJLaWlp2eiYJ554YpPvtWbNmqxZs6btcXNz85aWDUAP0lW9Rp8BAIBtQ7c9o60rTZ8+PbW1tW1LfX19uUsCoBfRZwAAYNvQbYO2QYMGJUlWrlzZbv3KlSvbnvtLO++8c/r06bNZY5Jk6tSpaWpqaluWL1++ldUD0BN0Va/RZwAAYNvQbYO23XffPYMGDcq8efPa1jU3N+eBBx7I6NGjNzqmX79+Oeigg9qNaW1tzbx58zY5JkmqqqpSU1PTbgGg9+uqXqPPAADAtqGs12h77bXXsnjx4rbHS5YsyaJFi7LTTjvl/e9/f84444xccMEF+cAHPpDdd98955xzToYMGdJ2t7gk+cQnPpG/+7u/yxe/+MUkyZQpUzJp0qQcfPDBOeSQQzJjxoysXr06kydP7urpAdAN6DUAAEBXKWvQ9tBDD+XII49sezxlypQkyaRJkzJ79ux85StfyerVq/PZz342r776ag477LDMmTMn1dXVbWOefvrpvPTSS22Px48fnxdffDHnnntuGhsbM2LEiMyZM2eDi1YDsG3QawAAgK5SUSqVSuUuortpbm5ObW1tmpqafL0HoACOq+3ZHwDFclzdkH0CUKyOHle77TXaAAAAeouZM2dmt912S3V1dUaNGpUHH3ywQ+Ouv/76VFRUtLukAQDdl6ANAACgE91www2ZMmVKpk2blocffjjDhw/PmDFj8sILL7zjuKVLl+ZLX/pSPvaxj3VRpQBsLUEbAABAJ7rsssty6qmnZvLkydl3330za9asbL/99rn66qs3OaalpSX/8A//kG984xvZY489urBaALaGoA0AAKCTrF27NgsXLkxDQ0PbusrKyjQ0NGTBggWbHPfNb34zAwcOzMknn9yh91mzZk2am5vbLQB0PUEbAABAJ3nppZfS0tKywZ2p6+rq0tjYuNEx9913X370ox/lyiuv7PD7TJ8+PbW1tW1LfX39VtUNwJYRtAEAAHQTq1atyoQJE3LllVdm55137vC4qVOnpqmpqW1Zvnx5J1YJwKb0LXcBAAAAvdXOO++cPn36ZOXKle3Wr1y5MoMGDdpg+6effjpLly7N2LFj29a1trYmSfr27Zsnn3wye+655wbjqqqqUlVVVXD1AGwuZ7QBAAB0kn79+uWggw7KvHnz2ta1trZm3rx5GT169AbbDxs2LI8++mgWLVrUtvzt3/5tjjzyyCxatMhXQgG6OWe0AQAAdKIpU6Zk0qRJOfjgg3PIIYdkxowZWb16dSZPnpwkmThxYoYOHZrp06enuro6+++/f7vxO+ywQ5JssB6A7kfQBgAA0InGjx+fF198Meeee24aGxszYsSIzJkzp+0GCcuWLUtlpS8bAfQGFaVSqVTuIrqb5ubm1NbWpqmpKTU1NeUuB6DHc1xtz/4AKJbj6obsE4BidfS46s8mAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABSg2wdtq1atyhlnnJFdd901/fv3z6GHHprf/va3m9x+/vz5qaio2GBpbGzswqoB6En0GgAAoAh9y13AuznllFPy+9//Ptdee22GDBmS6667Lg0NDfnDH/6QoUOHbnLck08+mZqamrbHAwcO7IpyAeiB9BoAAKAI3fqMtjfeeCO33HJLLr744vz1X/919tprr5x33nnZa6+9csUVV7zj2IEDB2bQoEFtS2Vlt54qAGWi1wAAAEXp1p8I1q1bl5aWllRXV7db379//9x3333vOHbEiBEZPHhwjjrqqPzmN795x23XrFmT5ubmdgsA24au6DX6DAAAbBu6ddA2YMCAjB49Oueff36ef/75tLS05LrrrsuCBQuyYsWKjY4ZPHhwZs2alVtuuSW33HJL6uvrc8QRR+Thhx/e5PtMnz49tbW1bUt9fX1nTQmAbqYreo0+AwAA24aKUqlUKncR7+Tpp5/OZz7zmdx7773p06dPPvShD+WDH/xgFi5cmMcff7xDr3H44Yfn/e9/f6699tqNPr9mzZqsWbOm7XFzc3Pq6+vT1NTU7to7AGyZ5ubm1NbWdtvjamf3Gn0GoHN19z5TDvYJQLE6elzt1me0Jcmee+6ZX/3qV3nttdeyfPnyPPjgg3nrrbeyxx57dPg1DjnkkCxevHiTz1dVVaWmpqbdAsC2o7N7jT4DAADbhm4ftL3tPe95TwYPHpxXXnklc+fOzfHHH9/hsYsWLcrgwYM7sToAegO9BgAA2Bp9y13Au5k7d25KpVL23nvvLF68OF/+8pczbNiwTJ48OUkyderUPPfcc7nmmmuSJDNmzMjuu++e/fbbL2+++Wauuuqq/Od//md+8YtflHMaAHRjeg0AAFCEbh+0NTU1ZerUqfnTn/6UnXbaKZ/85CfzrW99K9ttt12SZMWKFVm2bFnb9mvXrs2//uu/5rnnnsv222+fAw88MPfcc0+OPPLIck0BgG5OrwEAAIrQ7W+GUA4uHApQLMfV9uwPgGI5rm7IPgEoVq+5GQIAAEBPN3PmzOy2226prq7OqFGj8uCDD25y2yuvvDIf+9jHsuOOO2bHHXdMQ0PDO24PQPchaAMAAOhEN9xwQ6ZMmZJp06bl4YcfzvDhwzNmzJi88MILG91+/vz5Oemkk/LLX/4yCxYsSH19fY4++ug899xzXVw5AJtL0AYAANCJLrvsspx66qmZPHly9t1338yaNSvbb799rr766o1u/x//8R/5whe+kBEjRmTYsGG56qqr0tramnnz5nVx5QBsLkEbAABAJ1m7dm0WLlyYhoaGtnWVlZVpaGjIggULOvQar7/+et56663stNNOm9xmzZo1aW5ubrcA0PUEbQAAAJ3kpZdeSktLS+rq6tqtr6urS2NjY4de46tf/WqGDBnSLqz7S9OnT09tbW3bUl9fv1V1A7BlBG0AAADd1Le//e1cf/31ue2221JdXb3J7aZOnZqmpqa2Zfny5V1YJQBv61vuAgAAAHqrnXfeOX369MnKlSvbrV+5cmUGDRr0jmO/853v5Nvf/nbuueeeHHjgge+4bVVVVaqqqra6XgC2jjPaAAAAOkm/fv1y0EEHtbuRwds3Nhg9evQmx1188cU5//zzM2fOnBx88MFdUSoABXBGGwAAQCeaMmVKJk2alIMPPjiHHHJIZsyYkdWrV2fy5MlJkokTJ2bo0KGZPn16kuSiiy7Kueeem5/+9KfZbbfd2q7l9t73vjfvfe97yzYPAN6doA0AAKATjR8/Pi+++GLOPffcNDY2ZsSIEZkzZ07bDRKWLVuWyso/f9noiiuuyNq1a/OpT32q3etMmzYt5513XleWDsBmErQBAAB0si9+8Yv54he/uNHn5s+f3+7x0qVLO78gADqFa7QBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUIBuH7StWrUqZ5xxRnbdddf0798/hx56aH7729++45j58+fnQx/6UKqqqrLXXntl9uzZXVMsAD2SXgMAABSh2wdtp5xySu6+++5ce+21efTRR3P00UenoaEhzz333Ea3X7JkSY477rgceeSRWbRoUc4444yccsopmTt3bhdXDkBPodcAAABFqCiVSqVyF7Epb7zxRgYMGJCf/exnOe6449rWH3TQQTnmmGNywQUXbDDmq1/9au688878/ve/b1t34okn5tVXX82cOXM69L7Nzc2pra1NU1NTampqtn4iANu47nxcLUev6c77A6AnclzdkH0CUKyOHle79Rlt69atS0tLS6qrq9ut79+/f+67776NjlmwYEEaGhrarRszZkwWLFiwyfdZs2ZNmpub2y0AbBu6otfoMwAAsG3o1kHbgAEDMnr06Jx//vl5/vnn09LSkuuuuy4LFizIihUrNjqmsbExdXV17dbV1dWlubk5b7zxxkbHTJ8+PbW1tW1LfX194XMBoHvqil6jzwAAwLahWwdtSXLttdemVCpl6NChqaqqyve+972cdNJJqawsrvSpU6emqampbVm+fHlhrw1A99fZvUafAQCAbUPfchfwbvbcc8/86le/yurVq9Pc3JzBgwdn/Pjx2WOPPTa6/aBBg7Jy5cp261auXJmampr0799/o2OqqqpSVVVVeO0A9Ayd3Wv0GQAA2DZ0+zPa3vae97wngwcPziuvvJK5c+fm+OOP3+h2o0ePzrx589qtu/vuuzN69OiuKBOAHkyvAQAAtka3D9rmzp2bOXPmZMmSJbn77rtz5JFHZtiwYZk8eXKS9V/HmThxYtv2n/vc5/LMM8/kK1/5Sp544on84Ac/yI033pgzzzyzXFMAoJvTawAAgCJ0+6Ctqakpp512WoYNG5aJEyfmsMMOy9y5c7PddtslSVasWJFly5a1bb/77rvnzjvvzN13353hw4fn0ksvzVVXXZUxY8aUawoAdHN6DQAAUISKUqlUKncR3U1zc3Nqa2vT1NSUmpqaDo9bsWpFZi+anT++/McM6DcgJ+x7Qg57/2GpqKjoxGoBur8tPa72Vlu6P95c92ZufOzG3PvsvalIRQ7f7fB8at9PpbpvdSdWC9D96TMb2pJ9Umptzb3/9/LcvOBHeW3dGxm24175p09fkrrd9+/kagG6v44eV7v9GW09xeUPXJ7679bn7F+enev++7pc8dAV+evZf52PX/PxNL3ZVO7yAOjhfvvcb/P+774/k26flJ/87ieZ/bvZmXDbhOw2Y7c8vOLhcpcHwLuYOXNmdtttt1RXV2fUqFF58MEH33H7m266KcOGDUt1dXUOOOCA3HXXXZ1a3ysrluSvp+yQIxadkVn9Hs112y/O19bOyV/9+IDMuuzTnfreAL2JoK0Atz5+a/73nP+dllJLWkutWde6Luta1yVJfv3srzP+5vFlrhCAnmzFqhU56tqj8vIbLydJuz7z0usvpeGahryw+oVylgjAO7jhhhsyZcqUTJs2LQ8//HCGDx+eMWPG5IUXNn7svv/++3PSSSfl5JNPziOPPJJx48Zl3Lhx+f3vf98p9ZVaW/PJiz6UBbWrkiTr+qxfWivX/+/nV/2f/Py6czvlvQF6G0HbViqVSvnGr76Rimz866EtpZbMfXpuFjUu6trCAOg1Zj00K6+tfS0tpZYNnmsptaRpTVOueviqMlQGQEdcdtllOfXUUzN58uTsu+++mTVrVrbffvtcffXVG93+3/7t3/I3f/M3+fKXv5x99tkn559/fj70oQ/l+9//fqfU99t7fpJf7vhqWjbx6bCyNblg4WWd8t4AvY2gbSs9t+q5/PfK/04pm77UXd+Kvrn9idu7rigAepWb/nDTRkO2t7WWWnPTYzd1YUUAdNTatWuzcOHCNDQ0tK2rrKxMQ0NDFixYsNExCxYsaLd9kowZM2aT2yfJmjVr0tzc3G7pqNt/fWX6brrNpLUyeWCH1Xlh6WMdfk2AbZWgbSu9/tbr77pNRUVF3njrjS6oBoDeaPVbqwvZBoCu99JLL6WlpSV1dXXt1tfV1aWxsXGjYxobGzdr+ySZPn16amtr25b6+voO1/j6ujc28f2c9t5Y/WqHXxNgWyVo20r1NfV5z3bvecdt3mp9K/sPdKceALbMyEEj07ey7yaf71vZNyMHjezCigDobqZOnZqmpqa2Zfny5R0ee8Dg4XnrXT4Z1q5JBu8xfCurBOj9BG1bqf92/XPyyJPTp6LPRp+vSEV2qN4hJ+x3QhdXBkBv8YUPf6Ht5gcbs651XT7/4c93YUUAdNTOO++cPn36ZOXKle3Wr1y5MoMGDdromEGDBm3W9klSVVWVmpqadktHnTjh4gxYm1Rs4mo4fVqTUysOTr/+7+3wawJsqwRtBfjmkd/MsJ2HbRC29anokz6VffIf/+s/Ut23ukzVAdDTHbXHUTntw6clSbub71T+vzY+5SNTcsRuR5SjNADeRb9+/XLQQQdl3rx5betaW1szb968jB49eqNjRo8e3W77JLn77rs3uf3Wes+OA3PdPl9LZSnp8xfXauvTmuzfXJ1zp/ysU94boLcRtBWgtro2v/nMb3LWYWdlp/47JUkqKyozdu+xuf8z9+fYDxxb5goB6MkqKipy+TGX58fH/zj77rJv2/r96/bPtX93bb5z9HfKWB0A72bKlCm58sor85Of/CSPP/54Pv/5z2f16tWZPHlykmTixImZOnVq2/ann3565syZk0svvTRPPPFEzjvvvDz00EP54he/2Gk1/u2Eb+W+0Vfm2FV1qWxdv27nNyrytXwsvz776Qx435BOe2+A3qSiVCpt+naZ26jm5ubU1tamqalps065Ttbf+a3pzaZsv932qepb1UkVAvQsW3Nc7Y22dn80r2lORSoyoGpAJ1QH0PP0hD7z/e9/P5dcckkaGxszYsSIfO9738uoUaOSJEcccUR22223zJ49u237m266KWeffXaWLl2aD3zgA7n44otz7LEd/wP+1uyTN197NW+seiW1A+tT2WfT1wgF2JZ09LgqaNuIntCoAXoSx9X27A+AYjmubsg+AShWR4+rvjoKAAAAAAUQtAEAAABAAQRtAAAAAFAAQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEAB+pa7gO6oVColSZqbm8tcCUDv8Pbx9O3j67ZOnwEolj6zIb0GoFgd7TWCto1YtWpVkqS+vr7MlQD0LqtWrUptbW25yyg7fQagc+gzf6bXAHSOd+s1FSV/9tlAa2trnn/++QwYMCAVFRWbPb65uTn19fVZvnx5ampqOqHCbYP9WAz7sRj249YplUpZtWpVhgwZkspKVy3QZ96dOfZ8vX1+iTl2J/rMhram1/SUf/fuzn4shv1YDPtx63W01zijbSMqKyvzV3/1V1v9OjU1NX6AC2A/FsN+LIb9uOWcYfBn+kzHmWPP19vnl5hjd6HPtFdEr+kJ/+49gf1YDPuxGPbj1ulIr/HnHgAAAAAogKANAAAAAAogaOsEVVVVmTZtWqqqqspdSo9mPxbDfiyG/Uh3si38PJpjz9fb55eYI72Xf/di2I/FsB+LYT92HTdDAAAAAIACOKMNAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgbQvNnDkzu+22W6qrqzNq1Kg8+OCD77j9TTfdlGHDhqW6ujoHHHBA7rrrri6qtHvbnP04e/bsVFRUtFuqq6u7sNru5957783YsWMzZMiQVFRU5Pbbb3/XMfPnz8+HPvShVFVVZa+99srs2bM7vc7ubnP34/z58zf4WayoqEhjY2PXFMw2YVvoM5szxyuvvDIf+9jHsuOOO2bHHXdMQ0PDu+6T7mBz/x3fdv3116eioiLjxo3r3AK30ubO79VXX81pp52WwYMHp6qqKh/84Ae7/c/q5s5xxowZ2XvvvdO/f//U19fnzDPPzJtvvtlF1W4+v0tsu7aFPtMVfJ7Zeo5DxfCZpvsQtG2BG264IVOmTMm0adPy8MMPZ/jw4RkzZkxeeOGFjW5///3356STTsrJJ5+cRx55JOPGjcu4cePy+9//vosr7142dz8mSU1NTVasWNG2PPvss11YcfezevXqDB8+PDNnzuzQ9kuWLMlxxx2XI488MosWLcoZZ5yRU045JXPnzu3kSru3zd2Pb3vyySfb/TwOHDiwkypkW7Mt9JnNneP8+fNz0kkn5Ze//GUWLFiQ+vr6HH300Xnuuee6uPKO25I+lyRLly7Nl770pXzsYx/rokq3zObOb+3atTnqqKOydOnS3HzzzXnyySdz5ZVXZujQoV1cecdt7hx/+tOf5qyzzsq0adPy+OOP50c/+lFuuOGGfO1rX+viyjvO7xLbpm2hz3QFn2eK4ThUDJ9pupESm+2QQw4pnXbaaW2PW1paSkOGDClNnz59o9v//d//fem4445rt27UqFGlf/7nf+7UOru7zd2PP/7xj0u1tbVdVF3Pk6R02223veM2X/nKV0r77bdfu3Xjx48vjRkzphMr61k6sh9/+ctflpKUXnnllS6piW3PttBnNneOf2ndunWlAQMGlH7yk590VolbbUvmuG7dutKhhx5auuqqq0qTJk0qHX/88V1Q6ZbZ3PldccUVpT322KO0du3aripxq23uHE877bTSxz/+8XbrpkyZUvroRz/aqXUWxe8S245toc90BZ9niuc4VAyfacrLGW2bae3atVm4cGEaGhra1lVWVqahoSELFizY6JgFCxa02z5JxowZs8nttwVbsh+T5LXXXsuuu+6a+vr6HH/88Xnssce6otxew89isUaMGJHBgwfnqKOOym9+85tyl0MvsS30mS3tAf/T66+/nrfeeis77bRTZ5W5VbZ0jt/85jczcODAnHzyyV1R5hbbkvndcccdGT16dE477bTU1dVl//33z4UXXpiWlpauKnuzbMkcDz300CxcuLDtq2PPPPNM7rrrrhx77LFdUnNX6GnHGza0LfSZruDzTPn4eSyWzzTFE7RtppdeeiktLS2pq6trt76urm6T32VubGzcrO23BVuyH/fee+9cffXV+dnPfpbrrrsura2tOfTQQ/OnP/2pK0ruFTb1s9jc3Jw33nijTFX1PIMHD86sWbNyyy235JZbbkl9fX2OOOKIPPzww+UujV5gW+gzWzLHv/TVr341Q4YM2eAX7e5iS+Z433335Uc/+lGuvPLKrihxq2zJ/J555pncfPPNaWlpyV133ZVzzjknl156aS644IKuKHmzbckcP/3pT+eb3/xmDjvssGy33XbZc889c8QRR3Trr45uLr9L9HzbQp/pCj7PlI/jUDF8puk8fctdAHTU6NGjM3r06LbHhx56aPbZZ5/8+7//e84///wyVsa2Zu+9987ee+/d9vjQQw/N008/ne9+97u59tpry1gZbBu+/e1v5/rrr8/8+fN7zUWkV61alQkTJuTKK6/MzjvvXO5yOkVra2sGDhyYH/7wh+nTp08OOuigPPfcc7nkkksybdq0cpdXiPnz5+fCCy/MD37wg4waNSqLFy/O6aefnvPPPz/nnHNOucsDysznGboTn2k6j6BtM+28887p06dPVq5c2W79ypUrM2jQoI2OGTRo0GZtvy3Ykv34l7bbbruMHDkyixcv7owSe6VN/SzW1NSkf//+ZaqqdzjkkENy3333lbsMeoFtoc9sTQ/4zne+k29/+9u55557cuCBB3ZmmVtlc+f49NNPZ+nSpRk7dmzbutbW1iRJ37598+STT2bPPffs3KI3w5b8Gw4ePDjbbbdd+vTp07Zun332SWNjY9auXZt+/fp1as2ba0vmeM4552TChAk55ZRTkiQHHHBAVq9enc9+9rP5+te/nsrKnv9lEr9L9HzbQp/pCj7PlI/jUOfxmaYYPb/bd7F+/frloIMOyrx589rWtba2Zt68ee3+OvE/jR49ut32SXL33XdvcvttwZbsx7/U0tKSRx99NIMHD+6sMnsdP4udZ9GiRX4WKcS20Ge2tAdcfPHFOf/88zNnzpwcfPDBXVHqFtvcOQ4bNiyPPvpoFi1a1Lb87d/+bdsd1err67uy/He1Jf+GH/3oR7N48eK2ADFJ/vjHP2bw4MHdLmRLtmyOr7/++gZh2tvBYqlU6rxiu1BPO96woW2hz3QFn2fKx89j5/GZpiDlvhtDT3T99deXqqqqSrNnzy794Q9/KH32s58t7bDDDqXGxsZSqVQqTZgwoXTWWWe1bf+b3/ym1Ldv39J3vvOd0uOPP16aNm1aabvttis9+uij5ZpCt7C5+/Eb3/hGae7cuaWnn366tHDhwtKJJ55Yqq6uLj322GPlmkLZrVq1qvTII4+UHnnkkVKS0mWXXVZ65JFHSs8++2ypVCqVzjrrrNKECRPatn/mmWdK22+/fenLX/5y6fHHHy/NnDmz1KdPn9KcOXPKNYVuYXP343e/+93S7bffXnrqqadKjz76aOn0008vVVZWlu65555yTYFeZlvoM5s7x29/+9ulfv36lW6++ebSihUr2pZVq1aVawrvanPn+Je6+11HN3d+y5YtKw0YMKD0xS9+sfTkk0+Wfv7zn5cGDhxYuuCCC8o1hXe1uXOcNm1aacCAAaX/83/+T+mZZ54p/eIXvyjtueeepb//+78v1xTeld8ltk3bQp/pCj7PFMNxqBg+03QfgrYtdPnll5fe//73l/r161c65JBDSv/1X//V9tzhhx9emjRpUrvtb7zxxtIHP/jBUr9+/Ur77bdf6c477+ziirunzdmPZ5xxRtu2dXV1pWOPPbb08MMPl6Hq7uPtWzL/5fL2fps0aVLp8MMP32DMiBEjSv369SvtsccepR//+MddXnd3s7n78aKLLirtueeeperq6tJOO+1UOuKII0r/+Z//WZ7i6bW2hT6zOXPcddddN/r/02nTpnV94Zthc/8d/6fuHrSVSps/v/vvv780atSoUlVVVWmPPfYofetb3yqtW7eui6vePJszx7feeqt03nnntfWI+vr60he+8IXSK6+80vWFd5DfJbZd20Kf6Qo+z2w9x6Fi+EzTfVSUSr3kPHYAAAAAKCPXaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNuhhXnzxxQwaNCgXXnhh27r7778//fr1y7x588pYGQC9xTXXXJP3ve99WbNmTbv148aNy4QJE8pUFQC9hc809GYVpVKpVO4igM1z1113Zdy4cbn//vuz9957Z8SIETn++ONz2WWXlbs0AHqBN954I4MHD86VV16ZE044IUnywgsvZOjQofnFL36RI488sswVAtDT+UxDbyVogx7qtNNOyz333JODDz44jz76aH7729+mqqqq3GUB0Et84QtfyNKlS3PXXXclSS677LLMnDkzixcvTkVFRZmrA6A38JmG3kjQBj3UG2+8kf333z/Lly/PwoULc8ABB5S7JAB6kUceeSQf/vCH8+yzz2bo0KE58MADc8IJJ+Scc84pd2kA9BI+09AbuUYb9FBPP/10nn/++bS2tmbp0qXlLgeAXmbkyJEZPnx4rrnmmixcuDCPPfZY/umf/qncZQHQi/hMQ2/kjDbogdauXZtDDjkkI0aMyN57750ZM2bk0UcfzcCBA8tdGgC9yBVXXJEZM2bkqKOOylNPPZW5c+eWuyQAegmfaeitBG3QA335y1/OzTffnN/97nd573vfm8MPPzy1tbX5+c9/Xu7SAOhFmpqaMmTIkKxbty7XXHNNxo8fX+6SAOglfKaht/LVUehh5s+fnxkzZuTaa69NTU1NKisrc+211+bXv/51rrjiinKXB0AvUltbm09+8pN573vfm3HjxpW7HAB6CZ9p6M2c0QYAwCZ94hOfyH777Zfvfe975S4FAKDbE7QBALCBV155JfPnz8+nPvWp/OEPf8jee+9d7pIAALq9vuUuAACA7mfkyJF55ZVXctFFFwnZAAA6yBltAAAAAFAAN0MAAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAvz/e/gHozxT/KgAAAAASUVORK5CYII=", + "text/plain": [ + "
    " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "plt.figure(figsize=(15,5))\n", "plt.subplot(131)\n", @@ -615,9 +697,18 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-0.06089330922309355\n", + "-0.06055651225519789\n" + ] + } + ], "source": [ "print(calculate_exchange_tensor(pairs[0])[0]) # isotropic should be -82 meV\n", "print(calculate_exchange_tensor(pairs[1])[0]) # these should all be around -41.9 in the isotropic part\n", @@ -628,9 +719,18 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, - "outputs": [], + "outputs": [ + { + "ename": "SyntaxError", + "evalue": "invalid syntax (659628047.py, line 1)", + "output_type": "error", + "traceback": [ + "\u001b[0;36m Cell \u001b[0;32mIn[14], line 1\u001b[0;36m\u001b[0m\n\u001b[0;31m These are reasonably converged:\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" + ] + } + ], "source": [ "These are reasonably converged:\n", "\n", @@ -647,10 +747,6 @@ "metadata": {}, "outputs": [], "source": [ - "# symmetrizing Hamiltonian and overlap matrix to make them hermitian \n", - "# Check if exchange field has scalar part\n", - "# parallel over integrals\n", - "\n", "========================================\n", " \n", "Atom Angstrom\n", diff --git a/test.py b/test.py deleted file mode 100644 index 766b8cd..0000000 --- a/test.py +++ /dev/null @@ -1,515 +0,0 @@ -import os -from sys import stdout -from timeit import default_timer as timer - -from tqdm import tqdm - -os.environ["OMP_NUM_THREADS"] = "1" # export OMP_NUM_THREADS=4 -os.environ["OPENBLAS_NUM_THREADS"] = "1" # export OPENBLAS_NUM_THREADS=4 -os.environ["MKL_NUM_THREADS"] = "1" # export MKL_NUM_THREADS=6 -os.environ["VECLIB_MAXIMUM_THREADS"] = "1" # export VECLIB_MAXIMUM_THREADS=4 -os.environ["NUMEXPR_NUM_THREADS"] = "1" # export NUMEXPR_NUM_THREADS=6 - -import warnings - -import numpy as np -import sisl -from mpi4py import MPI -from numpy.linalg import inv - -from grogu.useful import * - -start_time = timer() - -# this cell mimicks an input file -fdf = sisl.get_sile("./lat3_791/Fe3GeTe2.fdf") -# this information needs to be given at the input!! -scf_xcf_orientation = np.array([0, 0, 1]) # z -# list of reference directions for around which we calculate the derivatives -# o is the quantization axis, v and w are two axes perpendicular to it -# at this moment the user has to supply o,v,w on the input. -# we can have some default for this -ref_xcf_orientations = [ - dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]), - dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]), - dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]), -] - -# human readable definition of magnetic entities -# magnetic_entities = [ -# dict(atom=0, ), -# dict(atom=1, ), -# dict(atom=2, ), -# dict(atom=3, l=2), -# dict(atom=4, l=2), -# dict(atom=5, l=2), -# ] -# pairs = [ -# dict(ai=3, aj=4, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV -# dict(ai=3, aj=5, Ruc=np.array([0, 0, 0])), # these should all be around -41.9 in the isotropic part -# dict(ai=4, aj=5, Ruc=np.array([0, 0, 0])), -# dict(ai=3, aj=0, Ruc=np.array([0, 0, 0])), -# dict(ai=3, aj=1, Ruc=np.array([0, 0, 0])), -# dict(ai=3, aj=2, Ruc=np.array([0, 0, 0])), -# ] -magnetic_entities = [ - dict(atom=3, l=2), - dict(atom=4, l=2), - dict(atom=5, l=2), -] - -# pair information -pairs = [ - dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV - dict( - ai=0, aj=2, Ruc=np.array([0, 0, 0]) - ), # these should all be around -41.9 in the isotropic part - dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])), - dict(ai=0, aj=1, Ruc=np.array([-1, 0, 0])), - dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])), - dict(ai=0, aj=1, Ruc=np.array([1, 0, 0])), - dict(ai=0, aj=2, Ruc=np.array([1, 0, 0])), - dict(ai=0, aj=1, Ruc=np.array([0, -1, 0])), - dict(ai=0, aj=2, Ruc=np.array([0, -1, 0])), - dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])), - dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])), - dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])), -] - -# Brilloun zone sampling and Green function contour integral -kset = 20 -kdirs = "xy" -ebot = -30 -eset = 50 -esetp = 1000 - - -# MPI parameters -comm = MPI.COMM_WORLD -size = comm.Get_size() -rank = comm.Get_rank() -root_node = 0 -if rank == root_node: - print("Number of nodes in the parallel cluster: ", size) - -simulation_parameters = dict( - path="Not yet specified.", - scf_xcf_orientation=scf_xcf_orientation, - ref_xcf_orientations=ref_xcf_orientations, - kset=kset, - kdirs=kdirs, - ebot=ebot, - eset=eset, - esetp=esetp, - parallel_size=size, -) - -# digestion of the input -# read in hamiltonian -dh = fdf.read_hamiltonian() -try: - simulation_parameters["geom"] = fdf.read_geometry() -except: - print("Error reading geometry.") - -# unit cell index -uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0]) - -setup_time = timer() - -NO = dh.no # shorthand for number of orbitals in the unit cell - -# preprocessing Hamiltonian and overlap matrix elements -h11 = dh.tocsr(dh.M11r) -h11 += dh.tocsr(dh.M11i) * 1.0j -h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") - -h22 = dh.tocsr(dh.M22r) -h22 += dh.tocsr(dh.M22i) * 1.0j -h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") - -h12 = dh.tocsr(dh.M12r) -h12 += dh.tocsr(dh.M12i) * 1.0j -h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") - -h21 = dh.tocsr(dh.M21r) -h21 += dh.tocsr(dh.M21i) * 1.0j -h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") - -sov = ( - dh.tocsr(dh.S_idx) - .toarray() - .reshape(NO, dh.n_s, NO) - .transpose(0, 2, 1) - .astype("complex128") -) - - -# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation -U = np.vstack( - [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])] -) -# This is the permutation that transforms ud1ud2 to u12d12 -# That is this transforms FROM SPIN BOX to ORBITAL BOX => U -# the inverse transformation is U.T u12d12 to ud1ud2 -# That is FROM ORBITAL BOX to SPIN BOX => U.T - -# From now on everything is in SPIN BOX!! -hh, ss = np.array( - [ - U.T @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) @ U - for i in range(dh.lattice.nsc.prod()) - ] -), np.array( - [ - U.T - @ np.block([[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]]) - @ U - for i in range(dh.lattice.nsc.prod()) - ] -) - - -# symmetrizing Hamiltonian and overlap matrix to make them hermitian -for i in range(dh.lattice.sc_off.shape[0]): - j = dh.lattice.sc_index(-dh.lattice.sc_off[i]) - h1, h1d = hh[i], hh[j] - hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2 - s1, s1d = ss[i], ss[j] - ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2 - -# identifying TRS and TRB parts of the Hamiltonian -TAUY = np.kron(np.eye(NO), tau_y) -hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())]) -hTRS = (hh + hTR) / 2 -hTRB = (hh - hTR) / 2 - -# extracting the exchange field -traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77 -XCF = np.array( - [ - np.array([f["x"] for f in traced]), - np.array([f["y"] for f in traced]), - np.array([f["z"] for f in traced]), - ] -) # equation 77 - -# Check if exchange field has scalar part -max_xcfs = abs(np.array(np.array([f["c"] for f in traced]))).max() -if max_xcfs > 1e-12: - warnings.warn( - f"Exchange field has non negligible scalar part. Largest value is {max_xcfs}" - ) - -H_and_XCF_time = timer() - -# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions -for i, mag_ent in enumerate(magnetic_entities): - parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes - magnetic_entities[i]["orbital_indeces"] = parsed - magnetic_entities[i]["spin_box_indeces"] = blow_up_orbindx( - parsed - ) # calculate spin box indexes - spin_box_shape = len( - mag_ent["spin_box_indeces"] - ) # calculate size for Greens function generation - - mag_ent["energies"] = [] # we will store the second order energy derivations here - - mag_ent["Gii"] = [] # Greens function - mag_ent["Gii_tmp"] = [] # Greens function for parallelization - mag_ent["Vu1"] = [ - list([]) for _ in range(len(ref_xcf_orientations)) - ] # These will be the perturbed potentials from eq. 100 - mag_ent["Vu2"] = [list([]) for _ in range(len(ref_xcf_orientations))] - for i in ref_xcf_orientations: - mag_ent["Gii"].append( - np.zeros((eset, spin_box_shape, spin_box_shape), dtype="complex128") - ) # Greens functions for every quantization axis - mag_ent["Gii_tmp"].append( - np.zeros((eset, spin_box_shape, spin_box_shape), dtype="complex128") - ) - -# for every site we have to store 2x3 Greens function (and the associated _tmp-s) -# in the 3 reference directions, because G_ij and G_ji are both needed -for pair in pairs: - spin_box_shape_i, spin_box_shape_j = len( - magnetic_entities[pair["ai"]]["spin_box_indeces"] - ), len( - magnetic_entities[pair["aj"]]["spin_box_indeces"] - ) # calculate size for Greens function generation - - pair["energies"] = [] # we will store the second order energy derivations here - - pair["Gij"] = [] # Greens function - pair["Gji"] = [] - pair["Gij_tmp"] = [] # Greens function for parallelization - pair["Gji_tmp"] = [] - - pair["Vij"] = [ - list([]) for _ in range(len(ref_xcf_orientations)) - ] # These will be the perturbed potentials from eq. 100 - pair["Vji"] = [list([]) for _ in range(len(ref_xcf_orientations))] - - for i in ref_xcf_orientations: - pair["Gij"].append( - np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype="complex128") - ) - pair["Gij_tmp"].append( - np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype="complex128") - ) # Greens functions for every quantization axis - pair["Gji"].append( - np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128") - ) - pair["Gji_tmp"].append( - np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128") - ) - -site_and_pair_dictionaries_time = timer() - -kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling -wkset = np.ones(len(kset)) / len(kset) # generate weights for k points -kpcs = np.array_split(kset, size) # split the k points based on MPI size -kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop", file=stdout) - -k_set_time = timer() - -# this will contain all the data needed to calculate the energy variations upon rotation -hamiltonians = [] - -# iterate over the reference directions (quantization axes) -for i, orient in enumerate(ref_xcf_orientations): - # obtain rotated exchange field - R = RotMa2b(scf_xcf_orientation, orient["o"]) - rot_XCF = np.einsum("ij,jklm->iklm", R, XCF) - rot_H_XCF = sum( - [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])] - ) - rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx] - - # obtain total Hamiltonian with the rotated exchange field - rot_H = hTRS + rot_H_XCF # equation 76 - - hamiltonians.append( - dict(orient=orient["o"], H=rot_H, rotations=[]) - ) # store orientation and rotated Hamiltonian - - for u in orient[ - "vw" - ]: # these are the infinitezimal rotations (for now) perpendicular to the quantization axis - Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H - - Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100 - Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100 - - for mag_ent in magnetic_entities: - mag_ent["Vu1"][i].append( - Vu1[:, mag_ent["spin_box_indeces"]][mag_ent["spin_box_indeces"], :] - ) # fill up the perturbed potentials (for now) based on the on-site projections - mag_ent["Vu2"][i].append( - Vu2[:, mag_ent["spin_box_indeces"]][mag_ent["spin_box_indeces"], :] - ) - - for pair in pairs: - ai = magnetic_entities[pair["ai"]][ - "spin_box_indeces" - ] # get the pair orbital sizes from the magnetic entities - aj = magnetic_entities[pair["aj"]]["spin_box_indeces"] - pair["Vij"][i].append( - Vu1[:, ai][aj, :] - ) # fill up the perturbed potentials (for now) based on the on-site projections - pair["Vji"][i].append(Vu1[:, aj][ai, :]) - -reference_rotations_time = timer() - -if rank == root_node: - print("Number of magnetic entities being calculated: ", len(magnetic_entities)) - print( - "We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site." - ) - print(f"The shape of the Hamiltonian and the Greens function is {NO}x{NO}.") -comm.Barrier() -# ---------------------------------------------------------------------- - -# make energy contour -# we are working in eV now ! -# and sisil shifts E_F to 0 ! -cont = make_contour(emin=ebot, enum=eset, p=esetp) -eran = cont.ze - -# ---------------------------------------------------------------------- -# sampling the integrand on the contour and the BZ -for k in kpcs[rank]: - wk = wkset[rank] # weight of k point in BZ integral - for i, hamiltonian_orientation in enumerate( - hamiltonians - ): # iterate over reference directions - # calculate Greens function - H = hamiltonian_orientation["H"] - HK, SK = hsk(H, ss, dh.sc_off, k) - Gk = inv(SK * eran.reshape(eset, 1, 1) - HK) - - # store the Greens function slice of the magnetic entities (for now) based on the on-site projections - for mag_ent in magnetic_entities: - mag_ent["Gii_tmp"][i] += ( - Gk[:, mag_ent["spin_box_indeces"]][..., mag_ent["spin_box_indeces"]] - * wk - ) - - for pair in pairs: - # add phase shift based on the cell difference - phase = np.exp(1j * 2 * np.pi * k @ pair["Ruc"].T) - - # get the pair orbital sizes from the magnetic entities - ai = magnetic_entities[pair["ai"]]["spin_box_indeces"] - aj = magnetic_entities[pair["aj"]]["spin_box_indeces"] - - # store the Greens function slice of the magnetic entities (for now) based on the on-site projections - pair["Gij_tmp"][i] += Gk[:, ai][..., aj] * phase * wk - pair["Gji_tmp"][i] += Gk[:, aj][..., ai] * phase * wk - -# summ reduce partial results of mpi nodes -for i in range(len(hamiltonians)): - for mag_ent in magnetic_entities: - comm.Reduce(mag_ent["Gii_tmp"][i], mag_ent["Gii"][i], root=root_node) - - for pair in pairs: - comm.Reduce(pair["Gij_tmp"][i], pair["Gij"][i], root=root_node) - comm.Reduce(pair["Gji_tmp"][i], pair["Gji"][i], root=root_node) - -green_function_inversion_time = timer() - -if rank == root_node: - # iterate over the magnetic entities - for tracker, mag_ent in enumerate(magnetic_entities): - # iterate over the quantization axes - for i, Gii in enumerate(mag_ent["Gii"]): - storage = [] - # iterate over the first and second order local perturbations - for Vu1, Vu2 in zip(mag_ent["Vu1"][i], mag_ent["Vu2"][i]): - # The Szunyogh-Lichtenstein formula - traced = np.trace((Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2) - # evaluation of the contour integral - storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we))) - - # fill up the magnetic entities dictionary with the energies - mag_ent["energies"].append(storage) - - # iterate over the pairs - for tracker, pair in enumerate(pairs): - # iterate over the quantization axes - for i, (Gij, Gji) in enumerate(zip(pair["Gij"], pair["Gji"])): - site_i = magnetic_entities[pair["ai"]] - site_j = magnetic_entities[pair["aj"]] - - storage = [] - # iterate over the first order local perturbations in all possible orientations for the two sites - for Vui in site_i["Vu1"][i]: - for Vuj in site_j["Vu1"][i]: - # The Szunyogh-Lichtenstein formula - traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2) - # evaluation of the contour integral - storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we))) - - # fill up the pairs dictionary with the energies - pairs[tracker]["energies"].append(storage) - - end_time = timer() - - print( - "############################### GROGU OUTPUT ###################################" - ) - print( - "================================================================================" - ) - print("Input file: ") - print(simulation_parameters["path"]) - print( - "Number of nodes in the parallel cluster: ", - simulation_parameters["parallel_size"], - ) - print( - "================================================================================" - ) - try: - print("Cell [Ang]: ") - print(simulation_parameters["geom"].cell) - except: - print("Geometry could not be read.") - print( - "================================================================================" - ) - print("DFT axis: ") - print(simulation_parameters["scf_xcf_orientation"]) - print("Quantization axis and perpendicular rotation directions:") - for ref in ref_xcf_orientations: - print(ref["o"], " --» ", ref["vw"]) - print( - "================================================================================" - ) - print("number of k points: ", simulation_parameters["kset"]) - print("k point directions: ", simulation_parameters["kdirs"]) - print( - "================================================================================" - ) - print("Parameters for the contour integral:") - print("Ebot: ", simulation_parameters["ebot"]) - print("Eset: ", simulation_parameters["eset"]) - print("Esetp: ", simulation_parameters["esetp"]) - print( - "================================================================================" - ) - print("Atomic informations: ") - print("") - print("") - print("Not yet specified.") - print("") - print("") - print( - "================================================================================" - ) - print("Exchange [meV]") - print( - "--------------------------------------------------------------------------------" - ) - print("Atom1 Atom2 [i j k] d [Ang]") - print( - "--------------------------------------------------------------------------------" - ) - for pair in pairs: - J_iso, J_S, D = calculate_exchange_tensor(pair) - J_iso = J_iso * sisl.unit_convert("eV", "meV") - J_S = J_S * sisl.unit_convert("eV", "meV") - D = D * sisl.unit_convert("eV", "meV") - - print(print_atomic_indices(pair, magnetic_entities, dh)) - print("Isotropic: ", J_iso) - print("DMI: ", D) - print("Symmetric-anisotropy: ", J_S) - print("") - - print( - "================================================================================" - ) - print("Runtime information: ") - print("Total runtime: ", end_time - start_time) - print( - "--------------------------------------------------------------------------------" - ) - print("Initial setup: ", setup_time - start_time) - print( - f"Hamiltonian conversion and XC field extraction: {H_and_XCF_time - setup_time:.3f} s" - ) - print( - f"Pair and site datastructure creatrions: {site_and_pair_dictionaries_time - H_and_XCF_time:.3f} s" - ) - print( - f"k set cration and distribution: {k_set_time - site_and_pair_dictionaries_time:.3f} s" - ) - print(f"Rotating XC potential: {reference_rotations_time - k_set_time:.3f} s") - print( - f"Greens function inversion: {green_function_inversion_time - reference_rotations_time:.3f} s" - ) - print( - f"Calculate energies and magnetic components: {end_time - green_function_inversion_time:.3f} s" - )