From 11f68b9a43dba516afa781daee82cfe9de8fc141 Mon Sep 17 00:00:00 2001 From: Daniel Pozsar Date: Mon, 11 Nov 2024 15:18:08 +0100 Subject: [PATCH] =?UTF-8?q?valami=20elromlott=20menet=20k=C3=B6zben!!!?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .pickle | Bin 0 -> 9308 bytes grgoupy_comondor.py | 1850 +++++++++++++++++++++++++++++++++++++++ src/grogupy/__init__.py | 6 + src/grogupy/grogu.py | 24 +- src/grogupy/io.py | 2 +- test.ipynb | 6 +- 6 files changed, 1871 insertions(+), 17 deletions(-) create mode 100644 .pickle create mode 100644 grgoupy_comondor.py diff --git a/.pickle b/.pickle new file mode 100644 index 0000000000000000000000000000000000000000..04109891cc103ed168ae6827d4c7ed72b3fa6b26 GIT binary patch literal 9308 zcmcIq3p`X?+aI?IA*EDDrI#3D$fd|-kEmq2hzQXbW|(oA(U>8U3Y9cnXi~!HGSTfE zk|W(LU8tx~DXK%Fi$uC`(na5zv9~%yHRpT3Z~uPJ{P(k;wVw5?wg2l`Ywtbl<##XD z#fgu#d^}Gng2JYR(>OGCBwwTtV=#jlVKhF^)x;wbUNoUnm<(E21S>X@!ZvYYMKi-# z6l$ailNG`UF_}j*pHEvxGcyTiF3yZAf>aY5*uN$oGPi%-IB{d1^TIA%%wXs@~UPQHjT}!g|k85f0WP`pN zB>MQD4$%#pbwPst&kzhnera4d6!a;-4r<7*=>2pQB&%Tk`f-0mBafBB5cBfmTC-s- ztWQ`iO#fGX#47)z$^XwjVhv#pOV>xNZ{OFVacMmH&`27GPg3G3gi;wWQSf9}!>K4w z3t(~h^bx#ak%5#j3L8!>`fSG$)pKYOZZ!sF{7rEB!9t>{i@aoc3ZkwO^Zo|A5sqtbv~Kg4wI z;eva%Dqx2bE4Rhgtz`1X_EZa^XcW0T^>9islO_&;$zgC9w1J5)Pl2KFNie5`@kz2g zH5NO7!J&ka8B8i85EdIe8I)j(fOCyu@hu?@J=gFCFla+4)1eqaF@{2bVgkh!$_yxG zP((5EMny(2nB)M~>Yp^N0h@-?V!s=#Bs$#t0k;fA4-`F7T#lj_iry&tptu4>UljdN zT!|tXMSm12CdqG9HZSYhrrkA^oh;8S>Vt7Br}wZWq*O8BC= zE7MWD$M;%esU2OZj~HBQPYrk0%y-on-$x2>1dT~!2a9gs5FQ4&^RjRUr|14#1b6y< zUKhlZW@lFQLL8=B{kH((gh+*Tt`Jp?&D=LYToJ$OY15;(NKRPjxkiZfb^c}|6n8_0 z^wYN1KRRi*7JmxhwAqb(?otFdK7`6XpYTl|r9`*Oivl=*z3OCKG#IdB`u9VoLL!S6 zN=77y#te+%)1yU`6IWJ?PkF44h$p8wLvJO?(|IIW`nJJ?AR;G!izD6&y0wNf8;T7S zTPSu=?4ck~h)^7$h`LDMJJ=7#v<*6genzQH(+{8pYKp#-JFB;u;j= zP+W^*Jc{d3T#w=g6cbQPL{S`nj+kT?7|G}de@3v%In)uPbB(2^Ch6U1i>S5lwCYnf zxHGk%VDw{yheB?ES2=12S{3TwA|mkM;jxLIU7_d{*Bo<;sC_9}6 zu{P@j=_TZm!^%EpJH*>Naj!)v4*KNJTEgLu)hl(`A-Mg^RMd-H0S;H;A22&M{zfb=_u|&aW9G)DDFcs6U8hPvr!bF zn1kYe6c3>I9gi%fx!7CiBVcwO1PM5M3qd-^VB_e+Gk)a|VZ1^zVlUKIR+3iBYE1M6$?#^ZYw zh~gYbhw~l(-KqjlvAIthdIBrWjfNyW4XN&fpQsp|NxIN(k-&7r60PdE9_=VKG(o;IulAP}rgOME`yq&rQKxMA=QXOa~B%BP5hx#!` zn(eK7#X!IBY@(*wDd1f7X6rEh#8T^7!V~G!j)Pe{E-ddFky!dlX5D`#Lc7J~jurdh zb*IB_c_KUB1EGlYxlFB@zd^frycAr>j7azxBVIJ{fsPi*I&;w-otz8VeTdtj(7)#5 zw->Uxu{v{!H)k1FcwHGy93P!idnREv(f&5@(%X&@ud8!4ZUiD?YV%~AhVX=~-XVPD zGh_J(QQ@G%P2cTE*%hbM=mO|T9H5NLP&Eg#>k(CR%7L{uK(cjAC!AmV+n35YWmq(l z(Rnx8HCHV0EEBaQS|&pjWkmTZL`ej)A<7VjNV!rlF=Vu;s(4uAo}a9iJ1+l29%8xU z7d(9MpXTKGjY815;oQ|lTiubuJnmhY+$^9z$<4)ay&D*p(tf+4MF{qVWNljND*!v` z7ugZZ1mNsZ*B4KI&jQa1x)NI-3BiZLkC`#$LSS|`Q7vAUf@qjq=i{3)K~8=LvI+4A zTVh7Gg*9^r>3Nakg6-L6e(^d5N>Bhw5Qwf4KB;&GoA3x9Qs3G)MoiTLP0T%Lp zJ0-LC&o_hG$Koket}j*&o>_*)BN-ok`PaqW%bX-wDD!{HLQ>YQURg;NVkr=2Ln&FP zS2TQBb(0W`*mHA&wzWG*``o*$e5C-??>p?h=dc@MVSC|v^D`m%)HKc1{d^YCO!cR@ z6=s1b5(ucE3&7;@Eh)KAh2Vk1Q#FkWAy{ppcxzFqKM;ISm%sE7+7q)^->~#YwD7@G z{h)n$;bNVHx0#^gz?qv9A(|@Yc^fWCVwj+}5V|tuI-j$b5mStExlJQM0Dr)}fA<051 zbD`+9+V=;xJc|^E1G4`7qk(f#2v%N=;n%)_ zepdD9yh_bd*#0(q*Jtzyf$(~%ZyG)qG#@=^m!G!Qeq(FsZCq{_>bxWl(Q&aoa^D4a$XA9afcG43dt8!{=^rC=aqA=9X>Sr9JmF z$mT~2@%sqG+$j^MYIhhB^(JhbrsrcSO)l)Mn)rvo`{~5wB$karhf!JaG|jQIpeymc z{8cWj)0UsO-gzV<**>Oku9j~Y7H?68#UmNt=D7S<+AAbkC}l4Eo`ItCAF#0DsBP`n zV-hR!zt4qzk8ac-g)E%?hWf#a8;D%d8ozO}AQMbquRvK`830;)u5vCw7Fv9)YRUN^ z0Q8$5ywkr5!1e|0s`x#bfX-BNdsr+4Nu=|WD4jx3zOaKhS8Jnv{W0d)zcR8x!nC)O zK~r`qv8H#8A7t5+eCv9r!`Z+p`1$c(Xis`q=75L#dgD1$sc(gVK2kURC}iR6nUTpE zK5Oj_wP;^FAj?kdIdZ1kH^siNx;)qmvamJ$xOmz@o*`$kJw$%(hT*bf)jtO z8O*{%v5khS45V9;ee_)A83vZb`ku<%`(Jv&Nr!v4oFR2|=4kJo^}>jF*gEf|=U%)t zEZovadKu?sOx)UX&T`=kqq3#eV|15G%fd9h)FWDpvy?i}}_z!g2@NcG(A%R*pq zJGLZLmIB`0SyFVhQ2=J1e&S<8^heT*Vr%|{cK=Diesvw0VA-1T-%L6)!B_5^?r>Pn z>$SS-^Xo#eY31v6RX(|(?LuQs(G&vWXLfUR9yoQ>)ziX61rgo?DcH z9CEPvmSOQo##inC^}zK1?Ym^ z^o!HaFIBZ57G&J$^9ij38=u?OIY1W1dKn~n&M+XF_>Z3F)*#is&{M0;@~zy;2mJXO0H z;P}F`e3y62!I71JCER;pgcv@|2;RTW7bM{?hCI$#3UrVNF9#ntpgL1Iugz>ZaH_bV zkxlUfJ+-)(6Hfbr`-`X83}5MqFzLT-4thTW+;`A7Zc+9Is=Db~@(SKyLwoG++*1@# zsBM`1P9_itK5eRKl$i_?IZHl&w(ymFM*d{zt;n@|1H^B(e({Xlsia?{IPM`*WK}_$ z(ouup0?LLT;WmRNihWz(-gZPjFP@)K(kgy0_={&`c_jzsuWl_!(;B;r1VRh=e7Quk z&$|DWV#u%dh4D24!>6xXfn@ubo@rej<|F?SVDU)C*ScnzD`NjyJQUyF)JVy~FVeBl z$<_b)@cBCbEemb)J@cY=E&yR;wY+B*$|A}xuXD7B4!|VdrRuTlC~zoyWclB2&LIAh zVpdr89H5qZh}*H&5yY)y>L}B5Hv$aTdM01c@ zcpjqFOxd?hw+_MO_m%h1tEFM#WCj0?=RemXN33@e*0O7XMNxFa4&^gqV*HYYOINR? z&8t5@sC_J+I@SGGPX7RBEv+Usxz6wq>{ZH=V6xj6O79jceL@QL^3!uu~GDj zoD4rGi~cA>f&ZC_94OjT(tC175Kk$T44q?FyS=j2XX#w!jfa0(p|Rn&KpL6N>i?U~ zK+Rq!>f4ZjJ-kW5Y#NP84&pM!-r?6{_;o%Kwu3q&M)!8i@ECbe^zD?2U0bc9kJs3h Iz~vJD2f}Q#@&Et; literal 0 HcmV?d00001 diff --git a/grgoupy_comondor.py b/grgoupy_comondor.py new file mode 100644 index 0000000..cd9f7c5 --- /dev/null +++ b/grgoupy_comondor.py @@ -0,0 +1,1850 @@ +# Copyright (c) [2024] [] +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +import warnings +from argparse import ArgumentParser +from sys import getsizeof +from timeit import default_timer as timer + +# use numpy number of threads one +try: + from threadpoolctl import threadpool_info, threadpool_limits + + user_api = threadpool_info()["user_api"] + threadpool_limits(limits=1, user_api=user_api) +except: + print("Warning: threadpoolctl could not make numpy use single thread!") + +from pickle import dump, load + +import numpy as np +import sisl +from mpi4py import MPI +from numpy.linalg import inv +from scipy.special import roots_legendre + +try: + from tqdm import tqdm + + tqdm_imported = True +except: + print("Please install tqdm for nice progress bar.") + tqdm_imported = False + + +# 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]]) + + +# list of accepted input parameters +ACCEPTED_INPUTS = [ + "infile", + "outfile", + "scf_xcf_orientation", + "ref_xcf_orientations", + "kset", + "kdirs", + "ebot", + "eset", + "esetp", + "parallel_solver_for_Gk", + "padawan_mode", +] + +DEFAULT_ARGUMENTS = dict( + infile=None, + outfile=None, + scf_xcf_orientation=np.array([0, 0, 1]), + 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])]), + ], + kset=2, + kdirs="xyz", + ebot=None, + eset=42, + esetp=1000, + parallel_solver_for_Gk=False, + padawan_mode=True, +) + + +def commutator(a, b): + """Shorthand for commutator. + + Commutator of two matrices in the mathematical sense. + + Args: + a : np.array_like + The first matrix + b : np.array_like + The second matrix + + Returns: + np.array_like + The commutator of a and b + """ + + return a @ b - b @ a + + +# define some useful functions +def hsk(H, ss, sc_off, k=(0, 0, 0)): + """Speed up Hk and Sk generation. + + Calculates the Hamiltonian and the Overlap matrix at a given k point. It is faster that the sisl version. + + Args: + H : np.array_like + Hamiltonian in spin box form + ss : np.array_like + Overlap matrix in spin box form + sc_off : list + supercell indexes of the Hamiltonian + k : tuple, optional + The k point where the matrices are set up. Defaults to (0, 0, 0) + + Returns: + np.array_like + Hamiltonian at the given k point + np.array_like + Overlap matrix at the given k point + """ + + # this two conversion lines are from the sisl source + k = np.asarray(k, np.float64) + k.shape = (-1,) + + # this generates the list of phases + phases = np.exp(-1j * 2 * np.pi * k @ sc_off.T) + + # phases applied to the hamiltonian + HK = np.einsum("abc,a->bc", H, phases) + SK = np.einsum("abc,a->bc", ss, phases) + + return HK, SK + + +def make_kset(dirs="xyz", NUMK=20): + """Simple k-grid generator to sample the Brillouin zone. + + Depending on the value of the dirs + argument k sampling in 1,2 or 3 dimensions is generated. + If dirs argument does not contain either of x, y or z + a kset of a single k-pont at the origin is returned. The total number of k points is the NUMK**(dimensions) + + + Args: + dirs : str, optional + Directions of the k points in the Brillouin zone. They are the three lattice vectors. Defaults to "xyz" + NUMK : int, optional + The number of k points in a direction. Defaults to 20 + + Returns: + np.array_like + An array of k points that uniformly sample the Brillouin zone in the given directions + """ + + # if there is no xyz in dirs return the Gamma point + if not (sum([d in dirs for d in "xyz"])): + return np.array([[0, 0, 0]]) + + kran = len(dirs) * [np.linspace(0, 1, NUMK, endpoint=False)] + mg = np.meshgrid(*kran) + 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]] + kset = np.array([dirsdict[d] for d in "xyz"]).T + + return kset + + +def make_contour(emin=-20, emax=0.0, enum=42, p=150): + """A more sophisticated contour generator. + + Calculates the parameters for the complex contour integral. It uses the + Legendre-Gauss quadrature method. It returns a class that contains + the information for the contour integral. + + Args: + emin : int, optional + Energy minimum of the contour. Defaults to -20 + emax : float, optional + Energy maximum of the contour. Defaults to 0.0, so the Fermi level + enum : int, optional + Number of sample points along the contour. Defaults to 42 + p : int, optional + Shape parameter that describes the distribution of the sample points. Defaults to 150 + + Returns: + ccont + Contains all the information for the contour integral + """ + x, wl = roots_legendre(enum) + R = (emax - emin) / 2 # radius + z0 = (emax + emin) / 2 # center point + y1 = -np.log(1 + np.pi * p) # lower bound + y2 = 0 # upper bound + + y = (y2 - y1) / 2 * x + (y2 + y1) / 2 + phi = (np.exp(-y) - 1) / p # angle parameter + ze = z0 + R * np.exp(1j * phi) # complex points for path + we = -(y2 - y1) / 2 * np.exp(-y) / p * 1j * (ze - z0) * wl + + # just an empty container class + class ccont: + pass + + cont = ccont() + cont.R = R + cont.z0 = z0 + cont.ze = ze + cont.we = we + cont.enum = enum + + return cont + + +def tau_u(u): + """Pauli matrix in direction u. + + Returns the vector u in the basis of the Pauli matrices. + + Args: + u : list or np.array_like + The direction + + Returns: + np.array_like + Arbitrary direction in the base of the Pauli matrices + """ + + # u is force to be of unit length + u = u / np.linalg.norm(u) + + return u[0] * TAU_X + u[1] * TAU_Y + u[2] * TAU_Z + + +# +def crossM(u): + """Definition for the cross-product matrix. + + It acts as a cross product with vector u. + + Args: + u : list or np.array_like + The second vector in the cross product + + Returns: + np.array_like + The matrix that represents teh cross product with a vector + """ + return np.array([[0, -u[2], u[1]], [u[2], 0, -u[0]], [-u[1], u[0], 0]]) + + +def RotM(theta, u, eps=1e-10): + """Definition of rotation matrix with angle theta around direction u. + + + Args: + theta : float + The angle of rotation + u : np.array_like + The rotation axis + eps : float, optional + Cutoff for small elements in the resulting matrix. Defaults to 1e-10 + + Returns: + np.array_like + The rotation matrix + """ + + 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) + ) + + # kill off small numbers + M[abs(M) < eps] = 0.0 + return M + + +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. + + Args: + a : np.array_like + First vector + b : np.array_like + Second vector + eps : float, optional + Cutoff for small elements in the resulting matrix. Defaults to 1e-10 + + Returns: + np.array_like + The rotation matrix with the above property + """ + + v = np.cross(a, b) + c = a @ b + M = np.eye(3) + crossM(v) + crossM(v) @ crossM(v) / (1 + c) + + # kill off small numbers + M[abs(M) < eps] = 0.0 + return M + + +def read_siesta_emin(eigfile): + """It reads the lowest energy level from the siesta run. + + It uses the .EIG file from siesta that contains the eigenvalues. + + Args: + eigfile : str + The path to the .EIG file + + Returns: + float + The energy minimum + """ + + # read the file + eigs = eigSileSiesta(eigfile).read_data() + + return eigs.min() + + +def int_de_ke(traced, we): + """It numerically integrates the traced matrix. + + It is a wrapper from numpy.trapz and it contains the + relevant constants to calculate the energy integral from + equation 93 or 96. + + Args: + traced : np.array_like + The trace of a matrix or a matrix product + we : float + The weight of a point on the contour + + Returns: + float + The energy calculated from the integral formula + """ + + return np.trapz(-1 / np.pi * np.imag(traced * we)) + + +def onsite_projection(matrix, idx1, idx2): + """It produces the slices of a matrix for the on site projection. + + The slicing is along the last two axes as these contains the orbital indexing. + + Args: + matrix : (..., :, :) np.array_like + Some matrix + idx : np.array_like + The indexes of the orbitals + + Returns: + np.array_like + Reduced matrix based on the projection + """ + + return matrix[..., idx1, :][..., idx2] + + +def calc_Vu(H, Tu): + """Calculates the local perturbation in case of a spin rotation. + + Args: + H : (NO, NO) np.array_like + Hamiltonian + Tu : (NO, NO) array_like + Rotation around u + + Returns: + Vu1 : (NO, NO) np.array_like + First order perturbed matrix + Vu2 : (NO, NO) np.array_like + Second order perturbed matrix + """ + + Vu1 = 1j / 2 * commutator(H, Tu) # equation 100 + Vu2 = 1 / 8 * commutator(commutator(Tu, H), Tu) # equation 100 + + return Vu1, Vu2 + + +def build_hh_ss(dh): + """It builds the Hamiltonian and Overlap matrix from the sisl.dh class. + + It restructures the data in the SPIN BOX representation, where NS is + the number of supercells and NO is the number of orbitals. + + Args: + dh : sisl.physics.Hamiltonian + Hamiltonian read in by sisl + + Returns: + hh : (NS, NO, NO) np.array_like + Hamiltonian in SPIN BOX representation + ss : (NS, NO, NO) np.array_like + Overlap matrix in SPIN BOX representation + """ + + 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()) + ] + ) + + return hh, ss + + +def setup_pairs_and_magnetic_entities( + magnetic_entities, pairs, dh, simulation_parameters +): + """It creates the complete structure of the dictionaries and fills some basic data. + + It creates orbital indexes, spin box indexes, coordinates and tags for magnetic entities. + Furthermore it creates the structures for the energies, the perturbed potentials and + the Greens function calculation. It dose the same for the pairs. + + Args: + pairs : dict + Contains the initial pair information + magnetic_entities : dict + Contains the initial magnetic entity information + dh : sisl.physics.Hamiltonian + Hamiltonian read in by sisl + simulation_parameters : dict + A set of parameters from the simulation + + Returns: + pairs : dict + Contains the initial information and the complete structure + magnetic_entities : dict + Contains the initial information and the complete structure + """ + + # for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions + for mag_ent in magnetic_entities: + parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes + mag_ent["orbital_indices"] = parsed + mag_ent["spin_box_indices"] = blow_up_orbindx( + parsed + ) # calculate spin box indexes + # if orbital is not set use all + if "l" not in mag_ent.keys(): + mag_ent["l"] = "all" + + # tag creation for one atom + if isinstance(mag_ent["atom"], int): + mag_ent["tags"] = [ + f"[{mag_ent['atom']}]{dh.atoms[mag_ent['atom']].tag}({mag_ent['l']})" + ] + mag_ent["xyz"] = [dh.xyz[mag_ent["atom"]]] + # tag creation for more atoms + if isinstance(mag_ent["atom"], list): + mag_ent["tags"] = [] + mag_ent["xyz"] = [] + # iterate over atoms + for atom_idx in mag_ent["atom"]: + mag_ent["tags"].append( + f"[{atom_idx}]{dh.atoms[atom_idx].tag}({mag_ent['l']})" + ) + mag_ent["xyz"].append(dh.xyz[atom_idx]) + + # calculate size for Greens function generation + spin_box_shape = len(mag_ent["spin_box_indices"]) + + # we will store the second order energy derivations here + mag_ent["energies"] = [] + + # These will be the perturbed potentials from eq. 100 + mag_ent["Vu1"] = [] # so they are independent in memory + mag_ent["Vu2"] = [] + + mag_ent["Gii"] = [] # Greens function + mag_ent["Gii_tmp"] = [] # Greens function for parallelization + for _ in simulation_parameters["ref_xcf_orientations"]: + # Rotations for every quantization axis + mag_ent["Vu1"].append([]) + mag_ent["Vu2"].append([]) + # Greens functions for every quantization axis + mag_ent["Gii"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape, spin_box_shape), + dtype="complex128", + ) + ) + mag_ent["Gii_tmp"].append( + np.zeros( + (simulation_parameters["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: + # calculate distance + xyz_ai = magnetic_entities[pair["ai"]]["xyz"] + xyz_aj = magnetic_entities[pair["aj"]]["xyz"] + xyz_aj = xyz_aj + pair["Ruc"] @ simulation_parameters["cell"] + pair["dist"] = np.linalg.norm(xyz_ai - xyz_aj) + + # calculate size for Greens function generation + spin_box_shape_i = len(magnetic_entities[pair["ai"]]["spin_box_indices"]) + spin_box_shape_j = len(magnetic_entities[pair["aj"]]["spin_box_indices"]) + # tag generation + pair["tags"] = [] + for mag_ent in [magnetic_entities[pair["ai"]], magnetic_entities[pair["aj"]]]: + tag = "" + # get atoms of magnetic entity + atoms_idx = mag_ent["atom"] + orbitals = mag_ent["l"] + + # if magnetic entity contains one atoms + if isinstance(atoms_idx, int): + tag += f"[{atoms_idx}]{dh.atoms[atoms_idx].tag}({orbitals})" + + # if magnetic entity contains more than one atoms + if isinstance(atoms_idx, list): + # iterate over atoms + atom_group = "{" + for atom_idx in atoms_idx: + atom_group += f"[{atom_idx}]{dh.atoms[atom_idx].tag}({orbitals})--" + # end {} of the atoms in the magnetic entity + tag += atom_group[:-2] + "}" + pair["tags"].append(tag) + 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"] = [] + for _ in simulation_parameters["ref_xcf_orientations"]: + # Greens functions for every quantization axis + pair["Gij"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape_i, spin_box_shape_j), + dtype="complex128", + ) + ) + pair["Gij_tmp"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape_i, spin_box_shape_j), + dtype="complex128", + ) + ) + pair["Gji"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape_j, spin_box_shape_i), + dtype="complex128", + ) + ) + pair["Gji_tmp"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape_j, spin_box_shape_i), + dtype="complex128", + ) + ) + + return pairs, magnetic_entities + + +def parallel_Gk(HK, SK, eran, eset): + """Calculates the Greens function by inversion. + + It calculates the Greens function on all the energy levels at the same time. + + Args: + HK : (NO, NO), np.array_like + Hamiltonian at a given k point + SK : (NO, NO), np.array_like + Overlap Matrix at a given k point + eran : (eset) np.array_like + Energy sample along the contour + eset : int + Number of energy samples along the contour + + Returns: + Gk : (eset, NO, NO), np.array_like + Green's function at a given k point + """ + + # Calculates the Greens function on all the energy levels + return inv(SK * eran.reshape(eset, 1, 1) - HK) + + +def sequential_GK(HK, SK, eran, eset): + """Calculates the Greens function by inversion. + + It calculates sequentially over the energy levels. + + Args: + HK : (NO, NO), np.array_like + Hamiltonian at a given k point + SK : (NO, NO), np.array_like + Overlap Matrix at a given k point + eran : (eset) np.array_like + Energy sample along the contour + eset : int + Number of energy samples along the contour + + Returns: + Gk : (eset, NO, NO), np.array_like + Green's function at a given k point + """ + + # creates an empty holder + Gk = np.zeros(shape=(eset, HK.shape[0], HK.shape[1]), dtype="complex128") + # fills the holder sequentially by the Greens function on a given energy + for j in range(eset): + Gk[j] = inv(SK * eran[j] - HK) + + return Gk + + +def remove_clutter_for_save(pairs, magnetic_entities): + """Removes unimportant data from the dictionaries. + + It is used before saving to throw away data that + is not needed for post processing. + + Args: + pairs : dict + Contains all the pair information + magnetic_entities : dict + Contains all the magnetic entity information + + Returns: + pairs : dict + Contains all the reduced pair information + magnetic_entities : dict + Contains all the reduced magnetic entity information + """ + + # remove clutter from magnetic entities and pair information + for pair in pairs: + del pair["Gij"] + del pair["Gij_tmp"] + del pair["Gji"] + del pair["Gji_tmp"] + for mag_ent in magnetic_entities: + del mag_ent["Gii"] + del mag_ent["Gii_tmp"] + del mag_ent["Vu1"] + del mag_ent["Vu2"] + + return pairs, magnetic_entities + + +def blow_up_orbindx(orb_indices): + """Function to blow up orbital indices to make SPIN BOX indices. + + Args: + orb_indices : np.array_like + These are the indices in ORBITAL BOX + + Returns: + orb_indices : np.array_like + These are the indices in SPIN BOX + """ + + orb_indices = np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten() + + return orb_indices + + +def spin_tracer(M): + """Spin tracer utility. + + This takes 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. + + + Args: + M : np.array_like + Traceable matrix + + Returns: + dict + It contains the traced matrix with "x", "y", "z" and "c" + """ + M11 = M[0::2, 0::2] + M12 = M[0::2, 1::2] + M21 = M[1::2, 0::2] + M22 = M[1::2, 1::2] + + M_o = dict() + M_o["x"] = M12 + M21 + M_o["y"] = 1j * (M12 - M21) + M_o["z"] = M11 - M22 + M_o["c"] = M11 + M22 + + return M_o + + +def parse_magnetic_entity(dh, atom=None, l=None, **kwargs): + """Function to define orbital indexes of a given magnetic entity. + + Args: + dh : sisl.physics.Hamiltonian + Hamiltonian from sisl + atom : integer or list of integers, optional + Defining atom (or atoms) in the unit cell forming the magnetic entity. Defaults to None + l : integer, optional + Defining the angular momentum channel. Defaults to None + + Returns: + list + The orbital indexes of the given magnetic entity + """ + # 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. + + +def calculate_anisotropy_tensor(mag_ent): + """Calculates the renormalized anisotropy tensor from the energies. + + Args: + mag_ent : dict + An element from the magnetic entities + + Returns: + K : np.array_like + elements of the anisotropy tensor + consistency_check : float + Absolute value of the difference from the consistency check + """ + + # get the energies + energies = mag_ent["energies"] + + K = np.zeros((3, 3)) + + # calculate the diagonal tensor elements + K[0, 0] = energies[1, 1] - energies[1, 0] + K[1, 1] = energies[0, 1] - energies[0, 0] + K[2, 2] = 0 + + # calculate the off-diagonal tensor elements + K[0, 1] = (energies[2, 0] + energies[2, 1]) / 2 - energies[2, 2] + K[1, 0] = K[0, 1] + K[0, 2] = (energies[1, 0] + energies[1, 1]) / 2 - energies[1, 2] + K[2, 0] = K[0, 2] + K[1, 2] = (energies[0, 0] + energies[0, 1]) / 2 - energies[0, 2] + K[2, 1] = K[1, 2] + + # perform consistency check + calculated_diff = K[1, 1] - K[0, 0] + expected_diff = energies[2, 0] - energies[2, 1] + consistency_check = abs(calculated_diff - expected_diff) + + return K, consistency_check + + +def calculate_exchange_tensor(pair): + """Calculates the exchange tensor from the energies. + + It produces the isotropic exchange, the relevant elements + from the Dzyaloshinskii-Morilla (Dm) tensor, the symmetric-anisotropy + and the complete exchange tensor. + + Args: + pair : dict + An element from the pairs + + Returns: + J_iso : float + Isotropic exchange (Tr[J] / 3) + J_S : np.array_like + Symmetric-anisotropy (J_S = J - J_iso * I ––> Jxx, Jyy, Jxy, Jxz, Jyz) + D : np.array_like + DM elements (Dx, Dy, Dz) + J : np.array_like + Complete exchange tensor flattened (Jxx, Jxy, Jxz, Jyx, Jyy, Jyz, Jzx, Jzy, Jzz) + """ + + # energies from rotations + energies = pair["energies"] + # Initialize output arrays + J = np.zeros((3, 3)) + D = np.zeros(3) + + # J matrix calculations + # J(1,1) = mean([DEij(2,2,2), DEij(2,2,3)]) + J[0, 0] = np.mean([energies[1, 3], energies[2, 3]]) + + # J(1,2) = -mean([DEij(1,2,3), DEij(2,1,3)]) + J[0, 1] = -np.mean([energies[2, 1], energies[2, 2]]) + J[1, 0] = J[0, 1] + + # J(1,3) = -mean([DEij(1,2,2), DEij(2,1,2)]) + J[0, 2] = -np.mean([energies[1, 1], energies[1, 2]]) + J[2, 0] = J[0, 2] + + # J(2,2) = mean([DEij(2,2,1), DEij(1,1,3)]) + J[1, 1] = np.mean([energies[0, 3], energies[2, 0]]) + + # J(2,3) = -mean([DEij(1,2,1), DEij(2,1,1)]) + J[1, 2] = -np.mean([energies[0, 1], energies[0, 2]]) + J[2, 1] = J[1, 2] + + # J(3,3) = mean([DEij(1,1,1), DEij(1,1,2)]) + J[2, 2] = np.mean([energies[0, 0], energies[1, 0]]) + + # D vector calculations + # D(1) = mean([DEij(1,2,1), -DEij(2,1,1)]) + D[0] = np.mean([energies[0, 1], -energies[0, 2]]) + + # D(2) = mean([DEij(2,1,2), -DEij(1,2,2)]) + D[1] = np.mean([energies[1, 2], -energies[1, 1]]) + + # D(3) = mean([DEij(1,2,3), -DEij(2,1,3)]) + D[2] = np.mean([energies[2, 1], -energies[2, 2]]) + + J_iso = np.trace(J) / 3 + # based on the grogu output pdf + # traceless symmetric exchange matrix: + # Jxx, Jyy, Jxy, Jxz, Jyz + J_S = np.array([J[0, 0] - J_iso, J[1, 1] - J_iso, J[0, 1], J[0, 2], J[1, 2]]) + + return J_iso, J_S, D, J + + +def read_fdf(path): + """It reads the simulation parameters, magnetic entities and pairs from the fdf. + + Args: + path : string + The path to the .fdf file + + Returns: + fdf_arguments : dict + The read input arguments from the fdf file + magnetic_entities : list + It contains the dictionaries associated with the magnetic entities + pairs : dict + It contains the dictionaries associated with the pair information + """ + + # read fdf file + fdf = fdfSileSiesta(path) + fdf_arguments = dict() + + InputFile = fdf.get("InputFile") + if InputFile is not None: + fdf_arguments["infile"] = InputFile + + OutputFile = fdf.get("OutputFile") + if OutputFile is not None: + fdf_arguments["outfile"] = OutputFile + + ScfXcfOrientation = fdf.get("ScfXcfOrientation") + if ScfXcfOrientation is not None: + fdf_arguments["scf_xcf_orientation"] = np.array(ScfXcfOrientation) + + XCF_Rotation = fdf.get("XCF_Rotation") + if XCF_Rotation is not None: + rotations = [] + # iterate over rows + for rot in XCF_Rotation: + # convert row to dictionary + dat = np.array(rot.split()[:9], dtype=float) + o = dat[:3] + vw = dat[3:].reshape(2, 3) + rotations.append(dict(o=o, vw=vw)) + fdf_arguments["ref_xcf_orientations"] = rotations + + Kset = fdf.get("INTEGRAL.Kset") + if Kset is not None: + fdf_arguments["kset"] = Kset + + Kdirs = fdf.get("INTEGRAL.Kdirs") + if Kdirs is not None: + fdf_arguments["kdirs"] = Kdirs + + # This is permitted because it means automatic Ebot definition + fdf_arguments["ebot"] = fdf.get("INTEGRAL.Ebot") + + Eset = fdf.get("INTEGRAL.Eset") + if Eset is not None: + fdf_arguments["eset"] = Eset + + Esetp = fdf.get("INTEGRAL.Esetp") + if Esetp is not None: + fdf_arguments["esetp"] = Esetp + + ParallelSolver = fdf.get("GREEN.ParallelSolver") + if ParallelSolver is not None: + fdf_arguments["parallel_solver_for_Gk"] = ParallelSolver + + PadawanMode = fdf.get("PadawanMode") + if PadawanMode is not None: + fdf_arguments["padawan_mode"] = PadawanMode + + Pairs = fdf.get("Pairs") + if Pairs is not None: + pairs = [] + # iterate over rows + for fdf_pair in Pairs: + # convert data + dat = np.array(fdf_pair.split()[:5], dtype=int) + # create pair dictionary + my_pair = dict(ai=dat[0], aj=dat[1], Ruc=np.array(dat[2:])) + pairs.append(my_pair) + + MagneticEntities = fdf.get("MagneticEntities") + if MagneticEntities is not None: + magnetic_entities = [] + for mag_ent in MagneticEntities: + row = mag_ent.split() + dat = [] + for string in row: + if string.find("#") != -1: + break + dat.append(string) + if dat[0] in {"Cluster", "cluster"}: + magnetic_entities.append(dict(atom=[int(_) for _ in dat[1:]])) + continue + elif dat[0] in {"AtomShell", "Atomshell", "atomShell", "atomshell"}: + magnetic_entities.append( + dict(atom=int(dat[1]), l=[int(_) for _ in dat[2:]]) + ) + continue + elif dat[0] in {"AtomOrbital", "Atomorbital", "tomOrbital", "atomorbital"}: + magnetic_entities.append( + dict(atom=int(dat[1]), orb=[int(_) for _ in dat[2:]]) + ) + continue + elif dat[0] in {"Orbitals", "orbitals"}: + magnetic_entities.append(dict(orb=[int(_) for _ in dat[1:]])) + continue + else: + raise Exception("Unrecognizable magnetic entity in .fdf!") + + return fdf_arguments, magnetic_entities, pairs + + +def process_input_args( + default_arguments, + fdf_arguments, + command_line_arguments, + accepted_inputs=ACCEPTED_INPUTS, +): + """It returns the final simulation parameters based on the inputs. + + The merging is done in the order of priority: + 1. command line arguments + 2. fdf arguments + 3. default arguments + + Args: + default_arguments : dict + Default arguments from grogupy + fdf_arguments : dict + Arguments read from the fdf input file + command_line_arguments : dict + Arguments from the command line + + Returns: + dict + The final simulation parameters + """ + + # iterate over fdf_arguments and update default arguments + for key, value in fdf_arguments.values(): + if value is not None and key in accepted_inputs: + default_arguments[key] = value + + # iterate over command_line_arguments and update default arguments + for key, value in command_line_arguments.values(): + if value is not None and key in accepted_inputs: + default_arguments[key] = value + + return default_arguments + + +def save_pickle(outfile, data): + """Saves the data in the outfile with pickle. + + Args: + outfile : str + Path to outfile + data : dict + Contains the data + """ + + # save dictionary + with open(outfile, "wb") as output_file: + dump(data, output_file) + + +def load_pickle(infile): + """Loads the data from the infile with pickle. + + Args: + infile : str + Path to infile + + Returns: + data : dict + A dictionary of data + """ + + # open and read file + with open(infile, "wb") as input_file: + data = load(data, input_file) + + return data + + +def print_job_description(simulation_parameters): + """It prints the parameters and the description of the job. + + + Args: + simulation_parameters : dict + It contains the simulations parameters + """ + + print( + "================================================================================================================================================================" + ) + print("Input file: ") + print(simulation_parameters["infile"]) + print("Output file: ") + print(simulation_parameters["outfile"]) + print( + "Number of nodes in the parallel cluster: ", + simulation_parameters["parallel_size"], + ) + if simulation_parameters["parallel_solver_for_Gk"]: + print("solver used for Greens function calculation: parallel") + else: + print("solver used for Greens function calculation: sequential") + print( + "================================================================================================================================================================" + ) + print("Cell [Ang]: ") + print(simulation_parameters["cell"]) + print( + "================================================================================================================================================================" + ) + print("DFT axis: ") + print(simulation_parameters["scf_xcf_orientation"]) + print("Quantization axis and perpendicular rotation directions:") + for ref in simulation_parameters["ref_xcf_orientations"]: + print(ref["o"], " --» ", ref["vw"]) + print( + "================================================================================================================================================================" + ) + print("Parameters for the contour integral:") + print("Number of k points: ", simulation_parameters["kset"]) + print("k point directions: ", simulation_parameters["kdirs"]) + if simulation_parameters["automatic_ebot"]: + print( + "Ebot: ", + simulation_parameters["ebot"], + " WARNING: This was automatically determined!", + ) + else: + print("Ebot: ", simulation_parameters["ebot"]) + print("Eset: ", simulation_parameters["eset"]) + print("Esetp: ", simulation_parameters["esetp"]) + print( + "================================================================================================================================================================" + ) + + +def print_parameters(simulation_parameters): + """It prints the simulation parameters for the grogu out. + + Args: + simulation_parameters : dict + It contains the simulations parameters + """ + + print( + "================================================================================================================================================================" + ) + print("Input file: ") + print(simulation_parameters["infile"]) + print("Output file: ") + print(simulation_parameters["outfile"]) + print( + "Number of nodes in the parallel cluster: ", + simulation_parameters["parallel_size"], + ) + print( + "================================================================================================================================================================" + ) + print("Cell [Ang]: ") + print(simulation_parameters["cell"]) + print( + "================================================================================================================================================================" + ) + print("DFT axis: ") + print(simulation_parameters["scf_xcf_orientation"]) + print("Quantization axis and perpendicular rotation directions:") + for ref in simulation_parameters["ref_xcf_orientations"]: + print(ref["o"], " --» ", ref["vw"]) + print( + "================================================================================================================================================================" + ) + print("Parameters for the contour integral:") + print("Number of k points: ", simulation_parameters["kset"]) + print("k point directions: ", simulation_parameters["kdirs"]) + print("Ebot: ", simulation_parameters["ebot"]) + print("Eset: ", simulation_parameters["eset"]) + print("Esetp: ", simulation_parameters["esetp"]) + print( + "================================================================================================================================================================" + ) + + +def print_atoms_and_pairs(magnetic_entities, pairs): + """It prints the pair and magnetic entity information for the grogu out. + + Args: + magnetic_entities : dict + It contains the data on the magnetic entities + pairs : dict + It contains the data on the pairs + """ + + print("Atomic information: ") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print( + "[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz" + ) + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + # iterate over magnetic entities + for mag_ent in magnetic_entities: + # iterate over atoms + for tag, xyz in zip(mag_ent["tags"], mag_ent["xyz"]): + # coordinates and tag + print(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}") + print("") + + print( + "================================================================================================================================================================" + ) + print("Anisotropy [meV]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print("Magnetic entity x [Ang] y [Ang] z [Ang]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + # iterate over magnetic entities + for mag_ent in magnetic_entities: + # iterate over atoms + for tag, xyz in zip(mag_ent["tags"], mag_ent["xyz"]): + # coordinates and tag + print(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}") + print("Consistency check: ", mag_ent["K_consistency"]) + print("K: # Kxx, Kxy, Kxz, Kyx, Kyy, Kyz, Kzx, Kzy, Kzz") + print(mag_ent["K"]) + print("") + + print( + "================================================================================================================================================================" + ) + print("Exchange [meV]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print("Magnetic entity1 Magnetic entity2 [i j k] d [Ang]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + # iterate over pairs + for pair in pairs: + # print pair parameters + print( + f"{pair['tags'][0]} {pair['tags'][1]} {pair['Ruc']} d [Ang] {pair['dist']}" + ) + # print magnetic parameters + print("Isotropic: ", pair["J_iso"], " # Tr[J] / 3") + print("") + print("DMI: ", pair["D"], " # Dx, Dy, Dz") + print("") + print( + "Symmetric-anisotropy: ", + pair["J_S"], + " # J_S = J - J_iso * I ––> Jxx, Jyy, Jxy, Jxz, Jyz", + ) + print("") + print("J: # Jxx, Jxy, Jxz, Jyx, Jyy, Jyz, Jzx, Jzy, Jzz") + print(pair["J"]) + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + + print( + "================================================================================================================================================================" + ) + + +def print_runtime_information(times): + """It prints the runtime information for the grogu out. + + Args: + times : dict + It contains the runtime data + """ + + print("Runtime information: ") + print(f"Total runtime: {times['end_time'] - times['start_time']} s") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print(f"Initial setup: {times['setup_time'] - times['start_time']} s") + print( + f"Hamiltonian conversion and XC field extraction: {times['H_and_XCF_time'] - times['setup_time']:.3f} s" + ) + print( + f"Pair and site datastructure creatrions: {times['site_and_pair_dictionaries_time'] - times['H_and_XCF_time']:.3f} s" + ) + print( + f"k set cration and distribution: {times['k_set_time'] - times['site_and_pair_dictionaries_time']:.3f} s" + ) + print( + f"Rotating XC potential: {times['reference_rotations_time'] - times['k_set_time']:.3f} s" + ) + print( + f"Greens function inversion: {times['green_function_inversion_time'] - times['reference_rotations_time']:.3f} s" + ) + print( + f"Calculate energies and magnetic components: {times['end_time'] - times['green_function_inversion_time']:.3f} s" + ) + + +def parse_command_line(): + """This function can read input from the command line.""" + + parser = ArgumentParser() + + parser.add_argument( + "-i", + "--input", + dest="infile", + default=None, + type=str, + help="Input file name", + required=True, + ) + parser.add_argument( + "-o", + "--output", + dest="outfile", + default=None, + type=str, + help="Output file name", + ) + parser.add_argument( + "--kset", + dest="kset", + default=None, + type=int, + help="k-space resolution of calculation", + ) + parser.add_argument( + "--kdirs", + dest="kdirs", + default=None, + type=str, + help="Definition of k-space dimensionality", + ) + parser.add_argument( + "--ebot", + dest="ebot", + default=None, + type=float, + help="Bottom energy of the contour", + ) + parser.add_argument( + "--eset", + dest="eset", + default=None, + type=int, + help="Number of energy points on the contour", + ) + parser.add_argument( + "--eset-p", + dest="esetp", + default=None, + type=int, + help="Parameter tuning the distribution on the contour", + ) + + parser.add_argument( + "--parallel-green", + dest="parallel_solver_for_Gk", + default=None, + type=bool, + help="Whether to use the parallel or sequential solver for Greens function", + ) + parser.add_argument( + "--padawan-mode", + dest="padawan_mode", + default=None, + type=bool, + help="If it is on it turns on extra helpful information for new users", + ) + + cmd_line_args = parser.parse_args() + + return cmd_line_args + + +def main(simulation_parameters, magnetic_entities, pairs): + # runtime information + times = dict() + times["start_time"] = timer() + + # MPI parameters + comm = MPI.COMM_WORLD + size = comm.Get_size() + rank = comm.Get_rank() + root_node = 0 + + if rank == root_node: + # include parallel size in simulation parameters + simulation_parameters["parallel_size"] = size + + # check versions for debugging + try: + print("sisl version: ", sisl.__version__) + except: + print("sisl version unknown.") + try: + print("numpy version: ", np.__version__) + except: + print("numpy version unknown.") + + # rename outfile + if not simulation_parameters["outfile"].endswith(".pickle"): + simulation_parameters["outfile"] += ".pickle" + + # if ebot is not given put it 1 eV under the smallest energy + if simulation_parameters["ebot"] is None: + try: + eigfile = simulation_parameters["infile"][:-3] + "EIG" + simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 1 + simulation_parameters["automatic_ebot"] = True + except: + print("Could not determine ebot.") + print("Parameter was not given and .EIG file was not found.") + else: + simulation_parameters["automatic_ebot"] = False + + # add third direction for off-diagonal anisotropy + for orientation in simulation_parameters["ref_xcf_orientations"]: + o1 = orientation["vw"][0] + o2 = orientation["vw"][1] + orientation["vw"].append((o1 + o2) / np.sqrt(2)) + + # read sile + fdf = sisl.get_sile(simulation_parameters["infile"]) + + # read in hamiltonian + dh = fdf.read_hamiltonian() + + # read unit cell vectors + simulation_parameters["cell"] = fdf.read_geometry().cell + + # unit cell index + uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0]) + + if rank == root_node: + print("\n\n\n\n\n") + print( + "#################################################################### JOB INFORMATION ###########################################################################" + ) + print_job_description(simulation_parameters) + print( + "################################################################################################################################################################" + ) + print("\n\n\n\n\n") + + times["setup_time"] = timer() + print(f"Setup done. Elapsed time: {times['setup_time']} s") + print( + "================================================================================================================================================================" + ) + + NO = dh.no # shorthand for number of orbitals in the unit cell + + # reformat Hamltonian and Overlap matrix for manipulations + hh, ss = build_hh_ss(dh) + + # 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"] / 2 for f in traced]), + np.array([f["y"] / 2 for f in traced]), + np.array([f["z"] / 2 for f in traced]), + ] + ) + + # check if exchange field has scalar part + max_xcfs = abs(np.array(np.array([f["c"] / 2 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}" + ) + + if rank == root_node: + times["H_and_XCF_time"] = timer() + print( + f"Hamiltonian and exchange field rotated. Elapsed time: {times['H_and_XCF_time']} s" + ) + print( + "================================================================================================================================================================" + ) + + # initialize pairs and magnetic entities based on input information + pairs, magnetic_entities = setup_pairs_and_magnetic_entities( + magnetic_entities, pairs, dh, simulation_parameters + ) + + if rank == root_node: + times["site_and_pair_dictionaries_time"] = timer() + print( + f"Site and pair dictionaries created. Elapsed time: {times['site_and_pair_dictionaries_time']} s" + ) + print( + "================================================================================================================================================================" + ) + + # generate k space sampling + kset = make_kset( + dirs=simulation_parameters["kdirs"], NUMK=simulation_parameters["kset"] + ) + + # generate weights for k points + wkset = np.ones(len(kset)) / len(kset) + + # split the k points based on MPI size + kpcs = np.array_split(kset, size) + + # use progress bar if available + if rank == root_node and tqdm_imported: + kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop") + + if rank == root_node: + times["k_set_time"] = timer() + print(f"k set created. Elapsed time: {times['k_set_time']} s") + print( + "================================================================================================================================================================" + ) + + # this will contain the three Hamiltonian in the + # reference directions needed to calculate the energy + # variations upon rotation + hamiltonians = [] + + # iterate over the reference directions (quantization axes) + for i, orient in enumerate(simulation_parameters["ref_xcf_orientations"]): + # obtain rotated exchange field and Hamiltonian + R = RotMa2b(simulation_parameters["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 + + # store the relevant information of the Hamiltonian + hamiltonians.append(dict(orient=orient["o"], H=rot_H)) + + # these are the rotations (for now) perpendicular to the quantization axis + for u in orient["vw"]: + # section 2.H + Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) + Vu1, Vu2 = calc_Vu(rot_H_XCF_uc, Tu) + + for mag_ent in magnetic_entities: + idx = mag_ent["spin_box_indices"] + # fill up the perturbed potentials (for now) based on the on-site projections + mag_ent["Vu1"][i].append(onsite_projection(Vu1, idx, idx)) + mag_ent["Vu2"][i].append(onsite_projection(Vu2, idx, idx)) + + if rank == root_node: + times["reference_rotations_time"] = timer() + print( + f"Rotations done perpendicular to quantization axis. Elapsed time: {times['reference_rotations_time']} s" + ) + print( + "================================================================================================================================================================" + ) + + # provide helpful information to estimate the runtime and memory + # requirements of the Greens function calculations + if rank == root_node: + print("Starting matrix inversions.") + if simulation_parameters["padawan_mode"]: + print("Padawan mode: ") + print(f"Total number of k points: {kset.shape[0]}") + print( + f"Number of energy samples per k point: {simulation_parameters['eset']}" + ) + print(f"Total number of directions: {len(hamiltonians)}") + print( + f"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * simulation_parameters['eset']}" + ) + print( + f"The shape of the Hamiltonian and the Greens function is {NO}x{NO}={NO*NO}" + ) + # https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv + # memory is O(64 n**2) for complex matrices + memory_size = getsizeof(hamiltonians[0]["H"].base) / 1024 + print( + f"Memory taken by a single Hamiltonian is: {getsizeof(hamiltonians[0]['H'].base) / 1024} KB" + ) + print(f"Expected memory usage per matrix inversion: {memory_size * 32} KB") + print( + f"Expected memory usage per k point for parallel inversion: {memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32} KB" + ) + print( + f"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32 / 1024} MB" + ) + print( + "================================================================================================================================================================" + ) + + # MPI barrier + comm.Barrier() + + # make energy contour + cont = make_contour( + emin=simulation_parameters["ebot"], + enum=simulation_parameters["eset"], + p=simulation_parameters["esetp"], + ) + eran = cont.ze + + # sampling the integrand on the contour and the BZ + for k in kpcs[rank]: + # weight of k point in BZ integral + wk = wkset[rank] + + # iterate over reference directions + for i, hamiltonian_orientation in enumerate(hamiltonians): + # calculate Hamiltonian and Overlap matrix in a given k point + H = hamiltonian_orientation["H"] + HK, SK = hsk(H, ss, dh.sc_off, k) + + if simulation_parameters["parallel_solver_for_Gk"]: + Gk = parallel_Gk(HK, SK, eran, simulation_parameters["eset"]) + else: + # solve Greens function sequentially for the energies, because of memory bound + Gk = sequential_GK(HK, SK, eran, simulation_parameters["eset"]) + + # store the Greens function slice of the magnetic entities + for mag_ent in magnetic_entities: + idx = mag_ent["spin_box_indices"] + mag_ent["Gii_tmp"][i] += onsite_projection(Gk, idx, idx) * 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_indices"] + aj = magnetic_entities[pair["aj"]]["spin_box_indices"] + + # store the Greens function slice of the magnetic entities + pair["Gij_tmp"][i] += onsite_projection(Gk, ai, aj) * phase * wk + pair["Gji_tmp"][i] += onsite_projection(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) + + if rank == root_node: + times["green_function_inversion_time"] = timer() + print( + f"Calculated Greens functions. Elapsed time: {times['green_function_inversion_time']} s" + ) + print( + "================================================================================================================================================================" + ) + + 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 + ) # this is the on site projection + # evaluation of the contour integral + storage.append(int_de_ke(traced, cont.we)) + # fill up the magnetic entities dictionary with the energies + magnetic_entities[tracker]["energies"].append(storage) + # convert to np array + magnetic_entities[tracker]["energies"] = np.array( + magnetic_entities[tracker]["energies"] + ) + + # 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 + ) # this is the on site projection + # evaluation of the contour integral + storage.append(int_de_ke(traced, cont.we)) + # fill up the pairs dictionary with the energies + pairs[tracker]["energies"].append(storage) + # convert to np array + pairs[tracker]["energies"] = np.array(pairs[tracker]["energies"]) + + # calculate magnetic parameters + for mag_ent in magnetic_entities: + K, consistency = calculate_anisotropy_tensor(mag_ent) + mag_ent["K"] = K * sisl.unit_convert("eV", "meV") + mag_ent["K_consistency"] = consistency + + for pair in pairs: + J_iso, J_S, D, J = calculate_exchange_tensor(pair) + pair["J_iso"] = J_iso * sisl.unit_convert("eV", "meV") + pair["J_S"] = J_S * sisl.unit_convert("eV", "meV") + pair["D"] = D * sisl.unit_convert("eV", "meV") + pair["J"] = J * sisl.unit_convert("eV", "meV") + + times["end_time"] = timer() + print("\n\n\n\n\n") + print( + "##################################################################### GROGU OUTPUT #############################################################################" + ) + + print_parameters(simulation_parameters) + print_atoms_and_pairs(magnetic_entities, pairs) + print( + "################################################################################################################################################################" + ) + print_runtime_information(times) + print("") + + # remove unwanted stuff before saving + pairs, magnetic_entities = remove_clutter_for_save(pairs, magnetic_entities) + # create output dictionary with all the relevant data + results = dict( + parameters=simulation_parameters, + magnetic_entities=magnetic_entities, + pairs=pairs, + runtime=times, + ) + + # save results + save_pickle(simulation_parameters["outfile"], results) + + print("\n\n\n\n\n") + + +if __name__ == "__main__": + # loading parameters + # it is not clear how to give grogu.fdf path... + # command_line_arguments = parse_command_line() + # fdf_arguments, magnetic_entities, pairs = read_fdf(command_line_arguments["infile"]) + + # right now we do not use any of these input, but it shows + # the order of priority when processing arguments + default_arguments = False + fdf_arguments = False + command_line_arguments = False + # simulation_parameters = process_input_args( + # default_arguments, fdf_arguments, command_line_arguments, ACCEPTED_INPUTS + # ) + + #################################################################################################### + # This is the input file for now # + #################################################################################################### + magnetic_entities = [ + dict(atom=3, l=2), + dict(atom=4, l=2), + dict(atom=5, l=2), + ] + pairs = [ + dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([-1, -1, 0])), + dict(ai=1, aj=2, Ruc=np.array([-1, -1, 0])), + dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-2, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])), + ] + simulation_parameters = dict() + simulation_parameters["infile"] = "" + simulation_parameters["outfile"] = "./" + simulation_parameters["scf_xcf_orientation"] = np.array([0, 0, 1]) + simulation_parameters["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])]), + ] + simulation_parameters["kset"] = 0 + simulation_parameters["kdirs"] = "" + simulation_parameters["ebot"] = None + simulation_parameters["eset"] = 0 + simulation_parameters["esetp"] = 0 + simulation_parameters["parallel_solver_for_Gk"] = False + simulation_parameters["padawan_mode"] = True + #################################################################################################### + # This is the input file for now # + #################################################################################################### + + main(simulation_parameters, magnetic_entities, pairs) diff --git a/src/grogupy/__init__.py b/src/grogupy/__init__.py index 0cb25f8..691ef46 100644 --- a/src/grogupy/__init__.py +++ b/src/grogupy/__init__.py @@ -20,3 +20,9 @@ """Docstring in init. """ + +from grogupy.core import * +from grogupy.globals import * +from grogupy.io import * +from grogupy.magnetism import * +from grogupy.utilities import * diff --git a/src/grogupy/grogu.py b/src/grogupy/grogu.py index 105d245..225f3f3 100644 --- a/src/grogupy/grogu.py +++ b/src/grogupy/grogu.py @@ -48,11 +48,7 @@ except: tqdm_imported = False -from .core import * -from .globals import * -from .io import * -from .magnetism import * -from .utilities import * +from grogupy import * def parse_command_line(): @@ -166,7 +162,7 @@ def main(simulation_parameters, magnetic_entities, pairs): if simulation_parameters["ebot"] is None: try: eigfile = simulation_parameters["infile"][:-3] + "EIG" - simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 1 + simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 0.1 simulation_parameters["automatic_ebot"] = True except: print("Could not determine ebot.") @@ -347,7 +343,7 @@ def main(simulation_parameters, magnetic_entities, pairs): f"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * simulation_parameters['eset']}" ) print( - f"The shape of the Hamiltonian and the Greens function is {NO}x{NO}={NO*NO}" + f"The shape of the Hamiltonian and the Greens function is {2*NO}x{2*NO}={4*NO*NO}" ) # https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv # memory is O(64 n**2) for complex matrices @@ -357,7 +353,7 @@ def main(simulation_parameters, magnetic_entities, pairs): ) print(f"Expected memory usage per matrix inversion: {memory_size * 32} KB") print( - f"Expected memory usage per k point for parallel inversion: {memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32} KB" + f"Expected memory usage per k point for parallel inversion: {memory_size * 32 * len(hamiltonians) * simulation_parameters['eset']} KB" ) print( f"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32 / 1024} MB" @@ -549,7 +545,9 @@ if __name__ == "__main__": dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])), ] simulation_parameters = dict() - simulation_parameters["infile"] = "" + simulation_parameters["infile"] = ( + "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf" + ) simulation_parameters["outfile"] = "./" simulation_parameters["scf_xcf_orientation"] = np.array([0, 0, 1]) simulation_parameters["ref_xcf_orientations"] = [ @@ -557,11 +555,11 @@ if __name__ == "__main__": 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])]), ] - simulation_parameters["kset"] = 0 - simulation_parameters["kdirs"] = "" + simulation_parameters["kset"] = 10 + simulation_parameters["kdirs"] = "xy" simulation_parameters["ebot"] = None - simulation_parameters["eset"] = 0 - simulation_parameters["esetp"] = 0 + simulation_parameters["eset"] = 600 + simulation_parameters["esetp"] = 1000 simulation_parameters["parallel_solver_for_Gk"] = False simulation_parameters["padawan_mode"] = True #################################################################################################### diff --git a/src/grogupy/io.py b/src/grogupy/io.py index f290145..77fbc38 100644 --- a/src/grogupy/io.py +++ b/src/grogupy/io.py @@ -209,7 +209,7 @@ def load_pickle(infile): """ # open and read file - with open(infile, "wb") as input_file: + with open(infile, "rb") as input_file: data = load(input_file) return data diff --git a/test.ipynb b/test.ipynb index b7f8626..a68317f 100644 --- a/test.ipynb +++ b/test.ipynb @@ -16,7 +16,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -24,14 +24,14 @@ "output_type": "stream", "text": [ "0.14.3\n", - "1.24.4\n" + "numpy version unknown.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "[Mac:30692] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Mac.501/jf.0/6815744/sm_segment.Mac.501.680000.0 could be created.\n" + "[Daniels-MacBook-Air.local:40959] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Daniels-MacBook-Air.501/jf.0/3756457984/sm_segment.Daniels-MacBook-Air.501.dfe70000.0 could be created.\n" ] } ],