From 2e6f4f6aab2278038dbaa96418738e387488eee4 Mon Sep 17 00:00:00 2001 From: mailhexu Date: Wed, 22 May 2024 22:19:20 +0200 Subject: [PATCH 1/7] adding splitOSC --- TB2J/abacus/splitSOC.py | 33 +++++++++++++++++++++++++++++++++ TB2J/mathutils/__init__.py | 1 + TB2J/mathutils/lowdin.py | 12 ++++++++++++ TB2J/mathutils/rotate_spin.py | 28 ++++++++++++++++++++++++++++ TB2J/pauli.py | 19 +++++++++++++++++++ 5 files changed, 93 insertions(+) create mode 100644 TB2J/abacus/splitSOC.py create mode 100644 TB2J/mathutils/__init__.py create mode 100644 TB2J/mathutils/lowdin.py create mode 100644 TB2J/mathutils/rotate_spin.py diff --git a/TB2J/abacus/splitSOC.py b/TB2J/abacus/splitSOC.py new file mode 100644 index 0000000..676471d --- /dev/null +++ b/TB2J/abacus/splitSOC.py @@ -0,0 +1,33 @@ +from TB2J.abacus.abacus_wrapper import AbacusWrapper, AbacusParser +from TB2J.mathutils.rotate_spin import rotate_Matrix_from_z_to_axis + + +class AbacusSplitSOCWrapper(AbacusWrapper): + """ + Abacus wrapper with Hamiltonian split to SOC and non-SOC parts + """ + + def __init__(self, **kwargs): + HR_soc = kwargs.pop("HR_soc", None) + super().__init__(**kwargs) + self.HR_soc = HR_soc + + def rotate_Hxc(self, axis): + """ + Rotate SOC part of Hamiltonian + """ + for iR, R in enumerate(self.Rlist): + self.HR[iR] = rotate_Matrix_from_z_to_axis(self.HR[iR], axis) + + +class AbacusSplitSOCParser: + """ + Abacus parser with Hamiltonian split to SOC and non-SOC parts + """ + + def __init__(self, outpath_nosoc=None, outpath_soc=None, binary=False): + self.outpath_nosoc = outpath_nosoc + self.outpath_soc = outpath_soc + self.binary = binary + parser_nosoc = AbacusParser(outpath=outpath_nosoc, binary=binary) + parser_soc = AbacusParser(outpath=outpath_soc, binary=binary) diff --git a/TB2J/mathutils/__init__.py b/TB2J/mathutils/__init__.py new file mode 100644 index 0000000..4d8f75e --- /dev/null +++ b/TB2J/mathutils/__init__.py @@ -0,0 +1 @@ +from .lowdin import Lowdin diff --git a/TB2J/mathutils/lowdin.py b/TB2J/mathutils/lowdin.py new file mode 100644 index 0000000..bee9ccc --- /dev/null +++ b/TB2J/mathutils/lowdin.py @@ -0,0 +1,12 @@ +import numpy as np +from scipy.linalg import inv, eigh + + +def Lowdin(S): + """ + Calculate S^(-1/2). + Which is used in lowind's symmetric orthonormalization. + psi_prime = S^(-1/2) psi + """ + eigval, eigvec = eigh(S) + return eigvec @ np.diag(np.sqrt(1.0 / eigval)) @ (eigvec.T.conj()) diff --git a/TB2J/mathutils/rotate_spin.py b/TB2J/mathutils/rotate_spin.py new file mode 100644 index 0000000..08fc70d --- /dev/null +++ b/TB2J/mathutils/rotate_spin.py @@ -0,0 +1,28 @@ +import numpy as np +from TB2J.pauli import pauli_block_all, s0, s1, s2, s3, gather_pauli_blocks + + +def rotate_Matrix_from_z_to_axis(M, axis, normalize=True): + """ + Given a spinor matrix M, rotate it from z-axis to axis. + The spinor matrix M is a 2x2 matrix, which can be decomposed as I, x, y, z components using Pauli matrices. + """ + MI, Mx, My, Mz = pauli_block_all(M) + axis = axis / np.linalg.norm(axis) + # M_new = s0* MI + Mz * (axis[0] * s1 + axis[1] * s2 + axis[2] * s3) *2 + M_new = gather_pauli_blocks( + MI, Mz * axis[0] * 2, Mz * axis[1] * 2, Mz * axis[2] * 2 + ) + return M_new + + +def test_rotate_Matrix_from_z_to_axis(): + M = np.array([[1.1, 0], [0, 0.9]]) + print(pauli_block_all(M)) + Mnew = rotate_Matrix_from_z_to_axis(M, [1, 1, 1]) + print(pauli_block_all(Mnew)) + print(Mnew) + + +if __name__ == "__main__": + test_rotate_Matrix_from_z_to_axis() diff --git a/TB2J/pauli.py b/TB2J/pauli.py index 0733919..1aa192d 100644 --- a/TB2J/pauli.py +++ b/TB2J/pauli.py @@ -143,7 +143,26 @@ def pauli_block_all(M): return MI, Mx, My, Mz +def gather_pauli_blocks(MI, Mx, My, Mz): + """ + Gather the I, x, y, z component of a matrix. + """ + return np.kron(MI, s0) + np.kron(Mx, s1) + np.kron(My, s2) + np.kron(Mz, s3) + + +def test_gather_pauli_blocks(): + M = np.random.rand(4, 4) + MI, Mx, My, Mz = pauli_block_all(M) + M2 = gather_pauli_blocks(MI, Mx, My, Mz) + print(M) + print(M2) + assert np.allclose(M, M2) + + def op_norm(M): + """ + Return the operator norm of a matrix. + """ return max(svd(M)[1]) From 9def387472901c69465a5507929d63ea86e1042c Mon Sep 17 00:00:00 2001 From: hexu Date: Thu, 23 May 2024 14:13:38 +0200 Subject: [PATCH 2/7] ongoing: rotate spin for splitSOC --- TB2J/abacus/abacus_wrapper.py | 18 ++++- TB2J/abacus/splitSOC.py | 123 ++++++++++++++++++++++++++++++++-- TB2J/green.py | 15 +---- TB2J/mathutils/fermi.py | 22 ++++++ TB2J/mathutils/rotate_spin.py | 13 +++- 5 files changed, 169 insertions(+), 22 deletions(-) create mode 100644 TB2J/mathutils/fermi.py diff --git a/TB2J/abacus/abacus_wrapper.py b/TB2J/abacus/abacus_wrapper.py index f7b0171..033401c 100644 --- a/TB2J/abacus/abacus_wrapper.py +++ b/TB2J/abacus/abacus_wrapper.py @@ -19,7 +19,7 @@ def __init__(self, HR, SR, Rlist, nbasis, nspin=1): self.R2kfactor = -2j * np.pi self.is_orthogonal = False self._name = "ABACUS" - self.HR = HR + self._HR = HR self.SR = SR self.Rlist = Rlist self.nbasis = nbasis @@ -27,6 +27,14 @@ def __init__(self, HR, SR, Rlist, nbasis, nspin=1): self.norb = nbasis * nspin self._build_Rdict() + @property + def HR(self): + return self._HR + + @HR.setter + def set_HR(self, HR): + self._HR = HR + def _build_Rdict(self): if hasattr(self, "Rdict"): pass @@ -74,6 +82,14 @@ def solve(self, k, convention=2): Hk, Sk = self.gen_ham(k, convention=convention) return eigh(Hk, Sk) + def solve_all(self, kpts, convention=2): + nk = len(kpts) + evals = np.zeros((nk, self.nbasis), dtype=float) + evecs = np.zeros((nk, self.nbasis, self.nbasis), dtype=complex) + for ik, k in enumerate(kpts): + evals[ik], evecs[ik] = self.solve(k, convention=convention) + return evals, evecs + def HSE_k(self, kpt, convention=2): H, S = self.gen_ham(tuple(kpt), convention=convention) evals, evecs = eigh(H, S) diff --git a/TB2J/abacus/splitSOC.py b/TB2J/abacus/splitSOC.py index 676471d..5521bcd 100644 --- a/TB2J/abacus/splitSOC.py +++ b/TB2J/abacus/splitSOC.py @@ -1,5 +1,16 @@ +import numpy as np from TB2J.abacus.abacus_wrapper import AbacusWrapper, AbacusParser from TB2J.mathutils.rotate_spin import rotate_Matrix_from_z_to_axis +from TB2J.kpoints import monkhorst_pack +from TB2J.mathutils.fermi import fermi +from copy import deepcopy +from scipy.spatial.transform import Rotation +import matplotlib.pyplot as plt + +# TODO List: +# - [x] Add the class AbacusSplitSOCWrapper +# - [x] Add the function to rotate the XC part +# - [x] Compute the band energy at arbitrary class AbacusSplitSOCWrapper(AbacusWrapper): @@ -7,17 +18,66 @@ class AbacusSplitSOCWrapper(AbacusWrapper): Abacus wrapper with Hamiltonian split to SOC and non-SOC parts """ - def __init__(self, **kwargs): + def __init__(self, *args, **kwargs): HR_soc = kwargs.pop("HR_soc", None) - super().__init__(**kwargs) + super().__init__(*args, **kwargs) + self._HR_copy = deepcopy(self._HR) self.HR_soc = HR_soc + @property + def HR(self): + return self._HR + self.HR_soc + def rotate_Hxc(self, axis): """ Rotate SOC part of Hamiltonian """ + # print("Before rotation:") + # print(self._HR[0][:2, :2]) for iR, R in enumerate(self.Rlist): - self.HR[iR] = rotate_Matrix_from_z_to_axis(self.HR[iR], axis) + self._HR[iR] = rotate_Matrix_from_z_to_axis(self._HR_copy[iR], axis) + + # print("After rotation:") + # print(self._HR[0][:2, :2]) + + +class RotateHam: + def __init__(self, model, kmesh, gamma=True): + self.model = model + self.kpts = monkhorst_pack(kmesh, gamma_center=gamma) + self.kweights = np.ones(len(self.kpts), dtype=float) / len(self.kpts) + + def get_band_energy(self, dm=False): + evals, evecs = self.model.solve_all(self.kpts) + eband = np.sum( + evals + * fermi(evals, self.model.efermi, width=0.01) + * self.kweights[:, np.newaxis] + ) + if dm: + density_matrix = self.model.get_density_matrix(evecs) + return eband + + def get_band_energy_vs_theta( + self, + angle_range=(0, np.pi * 2), + rotation_axis="y", + initial_direction=(0, 0, 1), + npoints=21, + ): + thetas = np.linspace(*angle_range) + es = [] + for theta in thetas: + axis = Rotation.from_euler(rotation_axis, theta).apply(initial_direction) + self.model.rotate_Hxc(axis) + e = self.get_band_energy() + es.append(e) + return thetas, es + + +def get_model_energy(model, kmesh, gamma=True): + ham = RotateHam(model, kmesh, gamma=gamma) + return ham.get_band_energy() class AbacusSplitSOCParser: @@ -29,5 +89,58 @@ def __init__(self, outpath_nosoc=None, outpath_soc=None, binary=False): self.outpath_nosoc = outpath_nosoc self.outpath_soc = outpath_soc self.binary = binary - parser_nosoc = AbacusParser(outpath=outpath_nosoc, binary=binary) - parser_soc = AbacusParser(outpath=outpath_soc, binary=binary) + self.parser_nosoc = AbacusParser(outpath=outpath_nosoc, binary=binary) + self.parser_soc = AbacusParser(outpath=outpath_soc, binary=binary) + spin1 = self.parser_nosoc.read_spin() + spin2 = self.parser_soc.read_spin() + if spin1 != "noncollinear" or spin2 != "noncollinear": + raise ValueError("Spin should be noncollinear") + + def parse(self): + nbasis, Rlist, HR, SR = self.parser_nosoc.Read_HSR_noncollinear() + nbasis2, Rlist2, HR2, SR2 = self.parser_soc.Read_HSR_noncollinear() + # print(HR[0]) + HR_soc = HR2 - HR + model = AbacusSplitSOCWrapper(HR, SR, Rlist, nbasis, nspin=2, HR_soc=HR_soc) + model.efermi = self.parser_nosoc.efermi + model.basis = self.parser_nosoc.basis + model.atoms = self.parser_nosoc.atoms + return model + + +def test_AbacusSplitSOCWrapper(): + path = "/Users/hexu/projects/2D_Fe/2D_Fe" + outpath_nosoc = f"{path}/Fe_SOC0/OUT.ABACUS" + outpath_soc = f"{path}/Fe_SOC1/OUT.ABACUS" + parser = AbacusSplitSOCParser( + outpath_nosoc=outpath_nosoc, outpath_soc=outpath_soc, binary=False + ) + model = parser.parse() + kmesh = [4, 4, 1] + e_z = get_model_energy(model, kmesh=kmesh, gamma=True) + print(e_z) + + model.rotate_Hxc([0, 0, 1]) + e_z = get_model_energy(model, kmesh=kmesh, gamma=True) + print(e_z) + + model.rotate_Hxc([1, 0, 0]) + e_x = get_model_energy(model, kmesh=kmesh, gamma=True) + print(e_x) + + r = RotateHam(model, kmesh) + # thetas, es = r.get_band_energy_vs_theta(angle_range=(0, np.pi*2), rotation_axis="z", initial_direction=(1,0,0), npoints=21) + thetas, es = r.get_band_energy_vs_theta( + angle_range=(0, np.pi * 2), + rotation_axis="y", + initial_direction=(0, 0, 1), + npoints=21, + ) + + plt.plot(thetas / np.pi, es, marker="o") + plt.savefig("E_along_z_x_z.png") + plt.show() + + +if __name__ == "__main__": + test_AbacusSplitSOCWrapper() diff --git a/TB2J/green.py b/TB2J/green.py index d5d0884..7862311 100644 --- a/TB2J/green.py +++ b/TB2J/green.py @@ -7,6 +7,8 @@ from pathos.multiprocessing import ProcessPool import sys import pickle +import warnings +from TB2J.mathutils.fermi import fermi MAX_EXP_ARGUMENT = np.log(sys.float_info.max) @@ -26,19 +28,6 @@ def eigen_to_G(evals, evecs, efermi, energy): ) -def fermi(e, mu, width=0.01): - """ - the fermi function. - .. math:: - f=\\frac{1}{\exp((e-\mu)/width)+1} - - :param e,mu,width: e,\mu,width - """ - - x = (e - mu) / width - return np.where(x < MAX_EXP_ARGUMENT, 1 / (1.0 + np.exp(x)), 0.0) - - def find_energy_ingap(evals, rbound, gap=4.0): """ find a energy inside a gap below rbound (right bound), diff --git a/TB2J/mathutils/fermi.py b/TB2J/mathutils/fermi.py new file mode 100644 index 0000000..66d60ec --- /dev/null +++ b/TB2J/mathutils/fermi.py @@ -0,0 +1,22 @@ +import numpy as np +import warnings +import sys + +MAX_EXP_ARGUMENT = np.log(sys.float_info.max) + + +def fermi(e, mu, width=0.01): + """ + the fermi function. + .. math:: + f=\\frac{1}{\exp((e-\mu)/width)+1} + + :param e,mu,width: e,\mu,width + """ + x = (e - mu) / width + # disable overflow warning + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ret = np.where(x < MAX_EXP_ARGUMENT, 1 / (1.0 + np.exp(x)), 0.0) + + return ret diff --git a/TB2J/mathutils/rotate_spin.py b/TB2J/mathutils/rotate_spin.py index 08fc70d..af0c5b2 100644 --- a/TB2J/mathutils/rotate_spin.py +++ b/TB2J/mathutils/rotate_spin.py @@ -10,9 +10,7 @@ def rotate_Matrix_from_z_to_axis(M, axis, normalize=True): MI, Mx, My, Mz = pauli_block_all(M) axis = axis / np.linalg.norm(axis) # M_new = s0* MI + Mz * (axis[0] * s1 + axis[1] * s2 + axis[2] * s3) *2 - M_new = gather_pauli_blocks( - MI, Mz * axis[0] * 2, Mz * axis[1] * 2, Mz * axis[2] * 2 - ) + M_new = gather_pauli_blocks(MI, Mz * axis[0], Mz * axis[1], Mz * axis[2]) return M_new @@ -23,6 +21,15 @@ def test_rotate_Matrix_from_z_to_axis(): print(pauli_block_all(Mnew)) print(Mnew) + M = np.array( + [ + [-9.90532976e-06 + 0.0j, 0.00000000e00 + 0.0j], + [0.00000000e00 + 0.0j, -9.88431291e-06 + 0.0j], + ] + ) + print(M) + print(rotate_Matrix_from_z_to_axis(M, [0, 0, 1])) + if __name__ == "__main__": test_rotate_Matrix_from_z_to_axis() From 82b09ef563173564248ae52b22f872dedd9dddbd Mon Sep 17 00:00:00 2001 From: mailhexu Date: Thu, 23 May 2024 17:01:42 +0200 Subject: [PATCH 3/7] updates --- TB2J/abacus/E_along_z_x_z.png | Bin 0 -> 28080 bytes TB2J/abacus/splitSOC.py | 85 +++++++++++++++++++++++++------- TB2J/mathutils/kR_convert.py | 90 ++++++++++++++++++++++++++++++++++ 3 files changed, 158 insertions(+), 17 deletions(-) create mode 100644 TB2J/abacus/E_along_z_x_z.png create mode 100644 TB2J/mathutils/kR_convert.py diff --git a/TB2J/abacus/E_along_z_x_z.png b/TB2J/abacus/E_along_z_x_z.png new file mode 100644 index 0000000000000000000000000000000000000000..9e41eb571669f133bb2686344cb7b9e51231741b GIT binary patch literal 28080 zcmd?RXH*s2wk=%boO6Z+f}((+ph$)#prQyWIg647$w_ipfPjdIh=PC$NX|J+P!vQ# zP%=nHg5>BM<=*Gq`|dk0eD9rqpRMh-b=6$8s%9E}^xnsexOzp6ikzJsK@h47>R4?A z!4E(XyfIQ@_(Za2=m-2y+V#Ai>ov#Qt{&#jw~))`u1@xjuJ$$-yzaN0U2GihoDw}F zdh!IXwX3U>%Q-PIhkrH@b#%58lV!=Lg`1E$sq4ES2#q=VA6}M1rVWBPa9zMEUH5$U zbJ$Don$@AyiVWG_V^TtV@*s{ele9 z8iPB17t0LttfR^dF4>(gwP3{Naf}gTY63#I@THe|)NL*A@;?h1N?my0H~m9aX!Pb| zm$e_yX5vVLcmr_oCpYZ+B0l<;EEjQvo{o-ApN&`%1AmUUAnI5w*76oB5`ng(B;bXs zaS7 -AWpbUVkSe-Du2XM9H{mK77f4T7rGc7F=fwgcoEVI2)CE1xfkjlW7Xsa|k+ChJA7 z5=Pk+O2P7}#={Z)@^Fnud)IRraj#!Lqz1WUJ!y~Y=D+R9P#d3{ix~ZVxP^;}2`eda z$yvyJ`O@r1M_i@vj>P`K!E&_;%pnP}eKI}59G$$N@Nkl)nq~Rgy`>s!YwOCNWtb5BjySUKN(+6?Od5^8O^Ry=%*QNX^YEp}{`SO|%gTWWl zOu^dQpP=TJ#T~iyxOB(JAdi_of{;9PB7Z9ZCg}$nA4KLG36F1pc}rN+@Gw_6HCOZd zgyVE2eYm&W18<644*XTtZ=zCQQnxZ5&Z8`v zIO*uW$KmYk?DXr$g@JG1%*UHUAH1{NUjL#=H7M)3!Yg>?8Sc#+;j?GYvdWn{IEb{x za7B%cnc06Y&7SPOiX+V3OvhoTinnCagQDml8}(pZhSTL-00I?os3nOfUze;jOx6ZgcHbX zJeJ#fv#*I;yd^~kth~HD%z->807m|qoDVY}A79wix#!fdlsn$NQ?jtY3=It}t*3^d z{lhATKZkTCinRw3QG7Sf)GaW{6XlV#Wn^b(e~^$6`r^e4slm=E=Z-1o^`C<*8EUbr z6?cXU@9nwEuD)ehSzX0L)|W>|3+uJBLozdkic3oHkc5N;^mSG;6%Mjn%N=DKG&8@9 zqk}iqnHcxrLF0-V7TZfCj<_SQw8>s+d+9o8K0Aw-eVtrPtYGjCKPf(;lA^IO7h<+J zP=o~BzfZY7+e3g{H86;rn6Ox1{KkNsv>!N$zS+IqbsLwOkK@bxJ8J_Kjy$rmvP-t9 z=-yPjDl~TGWskMfj}AJGr>DvC^YcGhpVq5(J;v=l%jdh=%7O$84H;>q$g*)r**A8? z@!y^MeB^Cg8``fS+S%AESEzOJ3_@SN6nOdS6Qvn5x*#A$IreDQl3#mLA= zSoQGmu(_3$(EA297!^0!O~f7JSFs2cgTk()}lUZC#zXc{3sBIoFP?ey5f&D)b$1e|Q4@85zMKgprpY zaq1c!YYHy%|Ls@1JD% z7fCEJuH*A2AU)&3Cw)SG8yBcyEy@1gB(vc``4T(R`fSIirv=@mv99582i6*4yj_lZ zq9X|R;xd(*Op>|&<$L^ovWC{ZI)ZUyW04nD;$+;l%dov#p8L0HcT z^o32{S?%;Z(X&eJusWO2fgQUPt!)AX&zf%2;LA}`&%JQh54t2~S*iFE*)RKM_%RY{ z=_HMxKWp!=#`?d7-Iy=oA`Dm{jlw>ysHmtzd3Ii2>#wQKqFbHJy?uS=^S#;IKjJmV ze*UCF^cwu+F#P<64u*Pqq+()Xh$8F=k3|f_PFZy|c6ZbG?+uO?nl%MckE&tl+oKf* zFCs%Tqg08P-dxqbF?yRLf^d8CrM#0%spN^{1wG>pWGqf=-A6`KI$A1(evOUWhhJzi zaVm_naB-IM7Zk*l2<>V%(#WmT{#sjGoy$yB&iukBAwkctdcSe{WAdQ=S~G>%7zD2N z-Hiq4(;~}v^y;s)1&$oSLu_qrU-a2Ve?OsDh?JF;scUEuB9^yrTZ}aZdcvMrL z+g=$*0=Bn3HB#ldq@|@-HZ~l3RCD^2v#(n~P<|q5$BgjG$#KD63;PVKoVRO~6zrwv zx0Z(K7#Mb_N9p zk1sFNM5}(WO{U_IYJB(Z9W85e>Ry`z(g=}caEO>rKg_htw6&Cmrj-5AlAi1~<)%2f zKQhMe5mH^zF5)w9VA!oA*!pB7#q;O*VWl7Yi@1h+KI9G08oWJZRIrwqbmGRKApDjc z&yV)@_G3#+F+o9({V$xUC@LulR!%!welLXXrD79Tdz#II(48WJYk5=_$lZ!Z5D!eTw%3$O&{>J{NTHXOLu zEWpm5G?8rx1EzlIQcGJj$Fi4UPfw5M+Q&08FN87kr0DWwyv{qPEH1@ef8XFaw?sCf zW)8v1=STHUE`_bd#}4i+%nZlC{-7YZevLh)`&}3D@53EVDWl(Bo~EjF%%h_uG&HN{ zQ_911$r<9THRzE*_r1v-OP!}Ko-i!`;M7oC%jxIm z_j5^zE=kq?qLj|t=grYp;<4ne?F!s%INpRD^Q73JqH8K&KXz~Jwh+J*^O%g@C03KT zw0fB0;_2TljnryC!M(lt5XrhOoOXR=_oY>tL3DAk#LL&O+3^jOmGJ?LsYbK+*(QJb zbamkScMFIh5NQ#_($bPwMux+-?7Q{w`}Y^^hpU;O6Z7O@>eoa+{`svWvsLPjrK2Ma zU>n0F7P=%I6Mab?UzhjRddv}odU8IbzD!rr`}!W6C2*9>%CK`#o3<^;bg@2OAR0@k zy!Pai=DQR{%>4EA6arR2f&l!?d(zLVU%8^Htxf;*L4sPEIEu#>e9! z&Q+XRT3UML_H3I!2g7@NuVYnJmGg8eD)G-sNijozd;WSSWGK6})~yhGDN%I1!2jU* z=PzGY-wCCfEqu*`5Sl&b@eH?QkYUouG5Xv=TriavNxbyIBgP~G8&q#1(tgH)G2oJJ z4mJz2q@*l7T13j>dkCOC6i|H0)NHT3J6U-9BYVpv;4U^cwygZapFigMa&&J_e@KE} zYYZV{c%hw5%fixcF(@uBPD(@Lc(;iZAtWZAK=;BFIiHh&)oXtr?&;OI|8$7OQ^!6p zcwJA|?Vl8xYLJwzA{fWp_+{Vw_Z`OBxuN-l=G3JI+|%8*E-`qZ$Fl@Vh7vLS*j}%T z!V9z;ty&OQh-EmAaZ~&Z6nUq6M>?RSs>(~z23|(BE?{q;NpG*$$+L3hGy~?^%TGn>1q$^I%yB%?) zE3V8n*mTk8c1H-_&1Y!i=D0fR!$M6t&)eLP+aKa8l6MyHK0cUoe2Cuo=$p4eZtps( zbO=1GWODC$pAFmbV5m8fE*jPDcSn)Vq!ww=>@4<6V;CX{oyRqEg9oJ+|{F zx;CY@+xZFZKbDMU@_Z^6i0w^-)VL$}?W<{Ms40f%w|*6mT}rSM2{y^#)%$qJpfzI# z$#ThEVYR5=f=1RwL-cLi3r8KIne_S_@8U;fO+wNfrw=b4ZB?(XZ4`BhH1DtA^|UmI+xqpFpZp zWihyBd;P1V-VauPu2){UiZm=T`*S&y>S9@n3XUwk2hFUN(gTsk+;Qh&ROJG&0>!(; zaw(%vBH`;R?j`Nj9du4^g`u%Sbn_$JI&6ms>t;i0#c_GlllUWdHqYwnMX=OsoWB=r zviJ+v(F|NcoHe0RW2EC1WtVKiE_Vsvdb5x|;`+^PVb0Z^AUdb=c)mGhGM%>P`2Mwg z={`*W)6(@1*1X>0`@!5+6nyD2t?6!xFteB zqi#=!Gk+PLD0vN+F4~oZsw0Oo^gRwY5rOn8aZj{D^3+_AMq8_$*`s+MnK1}-<1rI0 ztfkgjh7P$g^$Bc6W@Rz;m%{^!$m0turPV?y4Gk4Hx_|X(J+-3PIk!PxFa4O}RgT<* zDJKmLi%-!dehE8SzcpbqjD)xn?Mi=vW;%RdX(V4yT8@yqT!OfP8V$-Yh+8|yog~zv zU~sn=`h|)%=f)ololZPS*Au-wB~@rC z__;=8-Q4-@Un{qYE@3dI+Myc+UGfX@e=qOqxgCE^#IJYjxk8z17%k74+Wo61kum7f zabNS$4ko%dX?S-`d{&?2gr+H%nsG2U!`?AEX(`3Uht;(m;&`xyM^qEb9hOcJ&a<<( zibmi*FhV~jy_H#IgSi?IUcxovV^7+BVS!LPn{Zi9igl^WUAVzIX)fp9u! zWT8qYulb%3sofheJ8JB8!`5om*{!jrgqQ%apuRoPSe1j(#bT#%Xk(4~a^ls@wW8CT zx=7=f`4msy0mC0}a=U-fI)308Wn0(w+NvTbRn|9JKYW(h2LM{_HsfE)tFVY=WSYO5(bNL_5Bv6?!H*sIfm32V2YR2$A=>3O&Jib=t8x0cFYF6Q90(JGtz9s$LXJ2wR61N#)h+ zHm;_SPvorxK3lqM35f$XoJN8?FzzP7b~P=X7(T4LYdLmr?PR(T1l{ph%;n%_Zy%To z-i}>q5AnNLTQo@tX$dZ59kYMSfOgU0>=B+g$(xXo499D3nJVAD%T#OhbM<@6LE6&K zrWLE(zCU6-w^sSCR)3$=zRM&(?0;eLAx_;A=s?&0;>etpyAwo6R?=l%BLC85dWj8ORdg#VBcL zXeui!d4bgl4Gj%z?=LW=iNJoI&&$dpMIcg!AG`KiG5dNmeA%{)UEy#$>E_x_m>wUl7jEO6^*zvKHenH+TYFqpVpR7QW7dvy76P?A7W_&0Iatpow6#l z>gGn1aT03I#*%LN!8^%lvOUxFe4yNcV{>^l7V_1Q@Nf?9R5>3h$XufDT+6)D`Zj`L zd}@jaIqNiryw*qx0xFDOQj#fN@N#ntjl%uO_Sj@s|A^|yWzDI{#I@GKlDG+;Sjrzp z4a|y|xH0%^0Iafunf!_HFGT3X#>A0g4#}Q8ejIjqu*ZMo2vJKo4GHqLz5V9y`W$dI zfdvJQQBpu@U3^8Ke8#a2SOf%_o0}sgBMVf`Xae@mu(Zm1L%7;wIR*vNsR|7ry@@GV z1HP45O#*vWlegRXpE6Zz59Qm(Q`FOAeW`schEw{pZ8Bi0S=VzRp`mP#iAYFV zMn`$rqhI$)0hbaGA0N+mvRt*GKO;RIe`#sSqv{Oy>twLp%+>2!;Rub8tU^)6z|<$M z-Ijy8o-ZwAp;K;pbu+aCQ?zlUSgHUbBBHQV%j9nZ1C2n}5fSD6KCN~2Dk38zbCi)h zV+z<4hw_Msh^C<`m+{$w&4dV)BNIM;e06&=wxa3^5J!-lN8iDIop#MW~kdi@~e+ z5@1|CY%COxK)w`cf*fxCw3&5G>5RoypN^*_tn=M*hHUo5u$%Xv1mISsaD~2#3IPzh zY1-Ka4w>57S4~Z00d@oOghWkeFxR+-StXqMJaD8rgM-DlD40c!a1MjzYzMhRPk{~n zLlEsj!bBNpDG4OR_G06aeeYM4qMK4-eCIoy>!eB)Hm!_WR~IFGy~pUow%^gs&V9Jc z8Ic=s0Usb}h_=^43f})4)W5TkQ&8W0s#`W##IW3~?LJ#YZvV&#_vYqig;8u)PQ>MR zZSO0uSuY7ZWy}q!NDi(%-y^r??$QYP(@Lzd)vCMimjVfGGFkvWC2HGIMqG4B53^>F zu#gb5@1-_+rPg$S>VX3d0gBhG>zQ=e@c~`UBAwVG8jN#@xgkG zi;KIO(HyyAB3YcBA*2C#2xww!KFg85x1yL<~f_L;0$;>dV`X z?&CFX39ZEyy~_k@*QM4~f&4*Oy7DsoQ?Mixw2H# z47J_YTyfd_FLMv#;sSwhZ7KfWgM3)U#E^sN)q0Bi`ub8RVsm+qb&yKn!Tla-@riUGwA7c=Z*XZN;!)(Pp(J-0X_dgOQVQok2x- zzV{-k&&C@ODr1dRV7k5Gc`8cnY4MWB@61aX;2Kz|tDk)1-xp6yrT;zT#ayGwGWMyJ z9!RQ=m(NCESS=u4l-ca98ogw>-uSVJSWJ^5th$1v#gt`=-{bA?y#wCLa|6604ViQF;s59+Kv6%+qJ4M~`p$4~?~>qr+jpEzqC|d!!p=z?F18h>3=X0Fkm5oxFIABvf$T z7EJ~^d=zde_8&F#HPyVo?6V+Umu>g-Q$NGxu4j6^R_s>3DJ2$bOagd^?b(N8EglWi zN|?Yj5?pRP_Lrg9SPX;L*@PCa&k@!1@T>R;X;46gdnO*@r%)T z2_BIY_>a{su1o}pyNUlPBRRf$Q^Zv3%Tu8tXO`ayKkHj0IiBTG3h>yh`MDfY(9Ey z9^_|IzHJm3+boXr@k>^#_DarDZEfJMTpY5g*TkJ*CE$%?>O7GjbohF0mxBmTzm;<< z+Am(;C^T8@_t;)YM=|>A;-*zMqtMxnmm*8@z>i^e?tXGRLE^>6EwVTwfSSb?dlHfN zmpTRLYnkL}OT^jrokIe!KfJRRynQQq07;^-AyYxKY2XGK?(#e4niCjX0s>ww zmV|=5CI~CEcl0)QV@@q5@x@o_7{^9_CvLe{ht$%0JY7D6IR)8V@A+HIx=LcKVs?y# zmQ1=)*J#~rl3H(6&i!oL;)r8okV76nerz#P>jk8%KKpz?y7%-dtl*U!(gH2!c@USh?K)IQcDd%(_ zrr9~lxyv)G>PN=*qB}GqxRIHU$>G4Wn{}ti4V2ljP;;IQuu8pGJ}d<^R4efKMGk{f zi=z!ZpalabbVgPlG@2hXGaSROuRKTfoMymt)`quk-SV8z&OGV$%QAkECr;0=E?u%( z?5}+u7=-h0aMcn#Rs88KBy6FLLaMuK2db8@LAw|lvKzJe{BI<+Ipn0Qx;Rl~P$G$P z)NfDZ3MgWO8;OaDDx!%hIeh@ue^PHg>nH`?tYa+HMY-kpzn;8A0G$ zK70~bd9#YO*EnYR@%>P2b7yQvE%4aRG%r*$0)Z+HQ9wv3D1xsGx5n|S892m0eMzuyM&%gNbUg0Ss2aCj|%=T~+Y3OpJ|Xp>@tSt9SAtKf_& zmpI||Laet=j`c9>8uZj^1>3>6+40tFBN1xj6X;^L4b8$U1YFRnSu7cX^2BL7X z!{XxNi>pWPKL(Mm;+vqLAhF9#S7^qhBT3`hO!rODbhm$XikaWM>0-F?wXZL&wsvS% z)iUPL(f>;>4`7<6;cr)d*miuclo>SdBK1d3eo$1wPPYo+Il)&QOhS5mFKKF?X}wnQ zqR09Hqyt%5M_Y>9+Eh_$+_8Q~5VR4G!NEvS){R8Z4DIjlYiLRI$Vn$%3wud%q%iJc zp_A{6fvE=xk?bs_5T?-A6;|_EAjbA)-zivFpcXN$=cHghVeRqUq|RIEm3A56&a^bV zSFc`O&Ay&jel7cYRzZQ9zJ8LHn3=h`d5*)^eB-Dzm2k}I)4hecfaBF-xm#B!+fj7} z2{0)L=OryHDr(YHI0zx(_I%M=u9I&RajWw&+Y$`@DZ+=}Y?co8hG|XLXWNdSyuj(V zbY&a${#+0!mPzb`>a0O47cOQj~BTu=wQ{TDRm&9w$)R?=IZRVgvhqb%? ztqdOwTK~$?Eol3YbPCl0$ruEoXZb`$^IOYrT3FC>bB{~}(WU2SWz{l|zMRHKOe#}h zJLXPl7-2ZLY)c*a-tK&U81CsmLFScofAyBkJu!37?mOGNq2eIy627p}c(4CEQwc;< zPd^LCT+gF=Km&O0hq0|A^gVT7zRniTIY}X#hjCtU>H=XVjhDm4GHP{%oX`KEF5JkH zV~kX9858V!)tG5N=+EPzbxVKOWHP%j6%<)b+`=P_6%d7RM4<~YdxoS0IsM%6NP_cs z5eYHv?6*!gc4sam0mT#vTmlLbDjj2{UH}m`%oOA4souF$67ivFIc&Qb;&l#_8&2l& zaezNMu@FpyBSAzH-|z3=~OHnlnW`V`}z*d!-x`4s@u zK;q%x{<3%AvLU7jjiN6;Ob8FjhWKRg^5#@lZxCZ|^OjVSpGn78TcIe7Wa$*N&doCflw*LXU_3((aP(r1@!3>_*6WNO`I<)S2jpp8`7w zLN;&w{$qRTCPz|$#!92jIue@iWtDgxayZ2Q>N(HHHB z;bK7S#`S@Lte6)NR!fufc%=zOSH*QpfS!ywhJp~W?=_lUcfvlukB$_&Ljc+Ds(aRA z9&s${V#M+01bMc@01d%*q8Uv@RzZ5L#&|MC#dQ^oeKD8j(ad z@e#z1=C*RBmpw@2Y}`t?6K@aKbnm8R(Ti?EtLm6<{0C$9LM@@`7Z(*gZY(OS47Zh_MHak;iJ=0jV$qLQv1Awh?%jzh? zdvGzpw2Otm?u8M0FMPEBt(VxmtRK*LyH4hk^_Fc{$(VmvWC|Gdu8!y^g%SK-{sx!7xfGCf~#ez#Cn=E{$>taP0pwzz6}SOTXiR#$G#vJ-`DE{5&x*8on8yKuU$ zfho(eck&VyNtmF7B-%>nn*U^D_u@fDex8@pT6DPd$v{L^!d9CFpt69#3}ciYcG!((W}V;O;|?S#sL2{!#rm`b@`K1#d9<%u?*dlpIR+`*x6z8BeND@I!K%XX@<2#7|-h0Q{uO7 z9_0*GnDnj<$nw+$8hdgjJ{&BP{#im7%Li`JxgA3vM8MQAIN@$4?~v2{mzq_!MLfP7 zO#H?B@&tcztO_JpE)+c92TR5H^xDXu{GxAy5i1b%Vd4#Q=e9KOyxyPtfZyVb>B4c8 zWPLX`+Tl&;eQt6(YGmgh4c-oaqwYFmw8HlvkSaGa%Vp;;apob1@))S9joo-!BonVf zv9~R^vVB2@X#`FBcZwWiU7eBVC;pt%pbGca$*$DIKD?@vacZ{(&_Bu+e+c0v)k}Ot zgy@|VL?(BI*$zA2?D(9WOj@|cn}6wj2(Ocwb$SfO8#y5&(pu}aHa0y?9I6}$kl?zG zPWj;L@^a=V7O^n!&lwpR(ZXim+^pO(xzJw_j>Xt)=v=>UsaKGnPhC*IE9||w2#zuz zLPEkZkUGzNtEVrU(c6@K_MxY_OEnUJVHC_P{U9{XmM>QM3V^Q%imE;Co`2g|AUwwBdQOt_m{TBu?H(Q@qU?5I{i z7`B*r8?o~1mqhLska9Osp*SK!NmZ4&es}&jArVoMy=lbvwVjQT1~u5QTUt$qEN~w} zGNa`D&h)aN+}z{JR$4bcdZoAW8>OYE!C=O7?AS3@dEc`Lue>}r=&NdA&+@%@&#<(| zwhY{s1GQeKZhw4Ags5w3s;R4!o@S>4JJxXeM<8Q*>U=!jzI|(UxW5xN1I#bAsBtx+ zhsRdU8&ot0r5NbdJ9qAQlu=<`*z7mlbzDZtkR^p^HGk9kRNT|&U40hQ{xr?1!c9$y zJOLM%c-=2pj(hGiSCL%J=uu3&$_m6k5TQ!K2@_x=%PYTd`Eux$Ocg!7C+7MfNQ3n; zt;gEHp~Sp}_^h-v8$LhcyEW{AzRif&443Qj2-m^xMoa4NeXHJHrx#lDzF_R=$znMv z{OS34?5rR`O{cg-BNvr79tVvT>;zhgyz{qi086T>@D4v9I+|)|co+!MPq#V#HAT8T zNVbtEYTVouOfo(;7VvYhBBQdB)qlH%=c=}LklVu7vg+qxvjCUSAMiHXb%(wz2Xfhh zhsm$2A~`vZI!j`wk2%{T4@)i8g`y$c(|G={6Zn9)W4c^7D^uKYN&3AhH=x_w!hUMee>szWGq5Szb z7FLEMWFZ(7dkT^TfLBq?z`&rU;x8o2W}HBQ{D*_WZ4m5`8yg!(j~_o}V4|ltV4GZK z*KYy&87UQ&q`}hUWH8tgAUWZ992yd$NK5$R$B)w8A+WH9=;TNafi=xp<~JwSq5AG@ zVrKh#$l~-I4UQBs=B9}>&T!I)&88&%5HkC>0N09vu7RSPXwnC3yScTM2yt2( zk_P`1SROO_X38rn!t3M;x?aX)R9z2R4M4$R8h<#%8(8{eE6^sP&+K;b|)gIJ9k9E`V;r;871Muw~hnRkt!V78Ml%-`p6lDgSMv3PceKkY>O_F|qiq1eHDxe{fKssGVy5YJ8|{~w@NV+Xx*?+x1z z&!iauy*@cX)a)Cu2D|c<+c~nK8~92x8X#P1)#J-W~s27~A-)NyOkB_#(C@Ol!I< zel}^Kl2No)Ug|GJtjNb$=^%*ZH+!}fYtLmgKsrSv*g4v$dqZYR1parrj;^5A`CPuo0+aOuv%RalRH5gDi zS%uE&H~;|W^q0zP_r+E*$e4|rrn>YxMGsfp#Uq}(-eP0l30UG>>sfaUx+bu$LTF9_*9q105N^g+vqRRIxMuN<{quJ*rB)AUfwjby!#Xj*aS;+C2KX`yOe$zi6FpO=W| zWu~5=phjv5eOe&$y?c}bp_Qy>fux=75O+>->&_L=LtwN5#Hep6`L9KkV|#_s_g?7? zwJ;LS4~nXbHka_h%(K#^$&Fb!1_NE`w4z)DVNIV9ZTD$tZ);%XH31Wy8zqswUs2Omu%?U7qsv}VS(WaOj7b%@^QRTn(@t7KNMh&% zQ@UVWdyi~rSDS(8Lzp&qO6{*ZrNQK|m+>Y-C|wfi<7`e;9=L?!QIeL-GV%IN$K8!IQ_ocrtxVfP zHcv-yzWhfzQP4gKW{*T5)l>NClb{aNC5WPP`^K(k6-t7W;uo%-L?4L?SUQrNe87vegd&_q9k+MT)2493GO@H@Di7?R%F3Q|#=lW=Odj zPkL5)U+D%<$8P_-%m6G$n#!4 za2Wg-1`h|B6pot5k`bS$u0Td15sRRVcE@M3pn#2V2|PsX%2HHMo$MbTjWUkyiCi-e z!|t)w)q0BUUrth8KOU>l5~>Surbpy*0O(wbo3X1V1d( z{adz$=IX2<8G!GDo{NCDN6x_G-q;MzTyAXI>^N+sy~0%)A{gZwf&lYEcZ?mFuU%ev50*cBfE3-Wk&W{NWK!Js(? zgXYQy@C~A>#N45Fzn1YR$REd%07wESm{qV9uc}y9a(>uf?>~csC4$4&k)H(#8SK64 z7cMldjJMz+zYmuEQ500WHRAIR|1bvP8>;LR(8wYPTK=M}tn9SCVt$`njvXqAoLu*R z3hV>1gW==D1D{~nEpYbY;o)t6Nj)S)NGU1bx^&BhSB&}%KC00$?%15HI$GeMs7WI< zG}?U$GSXQo4f4&Ph8l+AOF64QS7kqNO^E^XmQGvte+6F6zvk+P-2(T8ZoW|y;3z&) zd~ynk)s22r)Ud!+>d~8kPR|Sa$Ww2|muF(7E@5w95KI_5=JsW$dV%JaGHr-qDewH# zejxS#61?7)w{DH^%%nQHpsrwmOJT{Wsj1nWLZ`}>yq8CF zufk4bY4and9LM@~J$&#SHJOUedx055^uLg;{{^D{uT8$Hof#P!;_i!j`%_~6cnBr8 zENfR+7h1G&xKZHW)ZU)?U&&U3w8e+EiAEWg(?D76GYokiMP_D#$%3FDJ&MX(=syl&b#BDIsml= z;uN9hS#yg9!!2iqxSgkrKSEQuufURRgq4fph-?A;g1`nluyeq4`(#@Nk`zM?C+hn< zd20jD2wbyvr}ie3?gpanp{=9SW`c)R3nzBHtHv^4G4RGYsSDXPPFA3NtbXLO-$#_^ zI{W=1Kuj%`)qxdmr{ef=U#%nZdsFm4KEgeEk3n?1qXHPyWXKsGxtjB3GW^=k1Tz-P z)%>xEIoPKI4{`i?rTuU1txEV07gi0LJQcV+C5Bl3A~T7we;f$!l3s=Fc?YIC7E=Rg zsd%7oM09Z3>b(>>KA6M|9Tl<6!Ej9nz4pS@MtU;?*1@h##B*tU1ugjQ-B5TF4L zc(6YrZx+R`IGKs=Yr(u8rRUvkBsX6(MM7$I0Z>NSnEj?ERAu+bRc*R`` zJQE7Y*xE~txEwwSfLG~K&${|1)NdcNe@QRwx}jF;n1$G)ntNyBU_8Xj@y^^Q@u$_QutKWe$A?@WG?s{hu4URvYA$?P6pfHGJSB>;z zpEpbcSELn=whZJY4NEgkVP`-%*c33>uN_bZa%~<#6hAa+kOqD7BfwC!WN`FQold;? zz~m~PnjstT=rD~j^gRqfznUFCZ{Q~lfCqXZof-sM*WhlN+js~d46vCC0pJi24vpu+ zgl!;=foa(&Sy(NaQ%UagCHb{~Bxh%o0*(~N$$`AEG3e+#asz^IZy?O#d_(`E6hVjv zxpsQR=kpNzIw;pyom=1fT_jch;z6jT3mr!IqE;u|L9GB8`y__7%-xYOAM61HK99{0 z%eo2tJD#4;82y`h4!W`Nkh}L}0Z!mPbcEEW%uY@Py{{&&Wk?P5b%=&dZ1@zhVzIaB z5vxfBG8P0%1Kt#EcxytT8nBo2luM!{2)h4*BZpV&2>4EMOMr#-%V5<&{&lNT3_&fS zzng_H!SGSNT%28H0n4_W(qqWEat*+Jd53$XOpErx4I`0@a#%}0@CChSn2QP;Ww4xM zMc%!+QSPi<@kaxQC+(w+lO|tx4ZChpt09cnAk>jAjy$}vcdl&C>iP7tdlELNqL~DY zKSxmuvd&+ak}&B<1>?YIdos1q(*}Tw1Beppf=F6;=*r%K%}BrHU`4Tz zg%A%X?TMO2m)M~LUQN6|=Ntc&0P%EJ7K6otxm0ZJ`^BTJJVW|ub}mB|dwLBNpzGKK zPfN_Rud{wuiHQqbf^8AYZ0RcIfmQ&3wU~p(6lOL2JhaS&iUQX#kjCK$R`pe#)5#Q>Td z7Kc$4j`qfbx57u5V-|>Dg5!f{+}=_f8x2}T?!gp2yCe6pryU<6{ML{EpYyXBZra+~ zs9X;oM%RU}RH&J>pyANAwCi2GkhtHzCsYliwRg+4YgfXTd6E)SyBNouFGw<8DfhJZ zq=8_Hq1fDZQt04vh`XU3EoVjHI2_P%ANT*=P++339|0vaP2kvsqOc~E6`Y@s1kP4D z(+(6MphTf8niq)U2(Wmf#o3;d(Pz+7miQ&H2sK(lU+)(-17SNEKF#l9TFmHshF!`l zn&LAK@KCvu5r$Ak05=OX4L-tm97c6LVqNHKR3M4DM{r~U3UsO(~%>bABE{^IIAzq z`cVXB$p!qWjl&10?U@r$VfDXu6(C4ZP|(=;cp$iId;9w{UcV+rt98i8loS;M#;b$? z1^^^&aYx-5P^b*$@gE6&E?F45)NP4b=y@D9{)x(})%aq%B>RW;OGtp`mW!EnG3Q^k zah_gYY$7E^MS(RnGUVjs-|cT&S&@V)o0kp^4i0{5>+53zqbq7G7@wUb1@;S!!IRLR zK>9|T%YYG9q`YTe&kl-(ubP2(P`Z6~0xS(k&cfQvnJdBIg!;af>0MU@on> z0*DDLrdL6dD6RdU+YQpvJ_4B+mykddtxAx6-Bu5f<{1F}1TNDb0--AI?v&HoL3xVL z)cXV|=}H>+^YW5bc|in$(26xP>PtOc+%Z|1G?vZREaS<&oF|e$)1-*Xehm2l#@;^$ z*X-IUZ^&7IH&nlT8H}_?5njlqfSiV75^T#Z!w2>b4n-czCV=b^1PZMmJbzAAW!Z87nXcFp0-nOQgB9g3v907JRvbOG(;MqL<^u@$erQpvfUxDx)xTkxz+aRFaV4Z$BLF7b&=#k*tGE{2(BxhCM)K^>TIDM6bBF}3VsuM480tbEep z`{@$~rKEs+)h}Lbet6{4NwC&z@9yFw=otwA8?!4^oSsI^-ak5q)`WwZ{yyj?0H=w; z3X2*wUS5bt5Y)2r{5e;8db-o}2R5|U9Zc(|=g0XwnUfzz7v<0w9(ZyiW&s5pnfLfH z@sN(gKJAD!0!H~yYFtMGUD=iOyU3aUq2TR)+r*c@;j&wq`M3GL9C4Rv4d_&}<`6RQ znc$+LPRVk^i#(U{m;ehccy!Uc=GU)Zv^+dCJUl$eWu;%4sea0Eqy^YMioh#fyZPO& zf;0HR36;MHnDMrzora5Ds`Z;o-FQbZw%udgcU1rO7y2zf(0@o8)VO1}qzyCVhqH94 zfe4TfcZLJbeS|uXCwuyxX(9Cnr`bWO##{$eCRc;U36j{++n)&3Ey+;?y3=j~0TaNZ z);wQ&z8LUwpOfo7x9e1jcF&l)7fRuipkseV>xa~`dk0aA~mjS$Ur z_0yol29P};SZ%=<39p+KFOFy@0%s`yy}x&F>t(eN7;P)8m9V|5)zX9JhqEb2z;D4q zmIohb2zcF&g^NBOdz&CJFsNSS#EB05DmI!Xdr9YuDm@I~S5;Od`_hhMGeA** zaA~enrT4G`+|VB%*m_Qz`mTvFV5fuN>F)P-zB=oh%~&o$mh!xDL;*uT>cZgD0@_(x z)DiVM-^+l@DB)h6*7d*6Sm@Dg3w!p7(M1%1z(3wnO@PSX7UNB#BUYTJzcOQ6BsMAh zl74I_q7IGM9b=(;dy2w%^q-@g&gr7&4g_u#WqsJWJI+^+79UozW70=~vA#KNW{NZl zMSufw4aha-WnXimiEc)4x;pvups>T%X7%Ok<<~q{SjLpKLOc3sp}SaCE=-&UKGQAW zOwjdw6m;||KX}h{ZP!PI_Gx_y?k&(>%px@}*+&nak|25#*#+=Dx-y4%8zzlhLcg3| zl^UQ1FTfhDs1zRk-0Bl(c)*SZdv~A}Pg$2$_8D16`xyXUfd4sAg(V6#n%A5t)<%sO zu4n#n=0Lxd?)-8~s$WFMzDQCpNy9U-HmIv#`57Urenb(HpAX#QV(vDUx1f%VS;*Az zCgUM8Rzi@pf9WaiLrOz6f91F7|s@|fIz9cCDYn7c8r80U^tP-t9wyz z(asSnOt_1dcLTp~pr4HE@P!Uk5f(RTTNXo&EBc(q!Dj13!v7G$r4KmhlFt0?ExZ?g z;AN}e=dYDHd7jtu?qS^)=Q{pLi5DD?uaGH&kEB61qQe57H`c_Z_N?$ zX9oYmCi3un)z#~Pj8ov5xYl*n(5U4)p=3brk-UZ;2dE8p( z82zbyFh`i#=qaEOZ2sh>%=KUmVI^g6I&HB;KpOgF#EKF=LA+znd6K<W*0!ft7|8+!_iQ0e6s+w-S+i$_Gg1=9{6jPuxn> z;oA*AjWoJ2cq=Vg9tW!z_O*~lacjw>4<0fSt=TR%Q-CowLc-7P2)+N&X`8mXLl2ri>UU!#Cjtl<1`u@uOns8FDNml@)VRF}<>lp_ z+}+zgrKuDFc@w|LKb!=9bPZo2I;bc@kTWE(hZb|kP2!J(9CAfpWV&78=HeT|x34Xu zjG5iU|1W4_T^$_@C`eh^+G+ zfj(}|H?GmBZs@K5i!y<#{QBCewG+Sb1u;2VS~)6i_cH#W%PlL-Uj@&D-N*ljV<%3> zaLYx~@vEX0`2g9F0LV1Zr#kC@C<}+>9N2zlJ9WEpaF8PvP91XSGc_@xt}>#|)1f$X zrlMv?=tVMTwch=zR4$vM{had|Q2xiK3`6=nYBczD`(#DO-{|?jJ6WRCca1>tUmYz$ zfRK`tqvt$XSXiJ{O;4Y)gAs#QK!6Y!Y=^HIaC*wp(ulZ>nn&xmD_<@Caa{<=|5t?=Z=Cw~x>uQ*$rIUlh@_+>s-A;-e;4R_&&UdlSnoKyNl2M7M3*W6fzKPKK5~8MD)A zWj+OzK%y%f@Ox_(s-2;5i?g#{5Ef8E=?Is5?BjnUx)v`TZ17j$_S-hk1#VMQZsU55 zWPwo?edlA5Q$uhxj1}o3uXkO*8IKeqM3(tGL3~i%ot_|Emfc#MjK#=o%GtZGCEzERrn*Y9to(J{GxTEpV#?(FB49ldhGS( zWLiMWz!0G73hS9WoFKp4^fZD4bfB|nL~n9(fXk->%-!;{ZIv+F#kY1DKVPReEi;^f*(jAH9*6T=`%sCj zuRAGgo#&FjP|NEs-Ajb~y9RAb6o8O9&!QTzjHz;xoPd(h`u#Zl-8l`sl;J1j;IQ2k zO(s?Fpr&xC@Z6UQ2CnGNZ^teD-ZS&!^#wxB%AVsSbXtZGzS{=mxw9NUP38^T0P13n z!FzypnN=3|lf8yI*&DhS9bJ3i=rXo=NJkw0TX$z3P37A6@f|{KW1)~CN=T@PB2wZ| zk<%oEO@<~rge`O0C`C>wB@`7UQzUFtR8oYLdB{9vN`{Qj_jbGHXa+J;QtBTKJ#nBoUw!)rvANt6sQCoF zpI!RYY_yXIbRmKW*4bInzF1p|os#WZ{+TTM;QWiXpJnD2_luPr=KvcD!1(Xf_3KT) zJH6kHcrae=*g?NnKK46&Z0c`$A5^U;Jj`KTz1?R6A5#t#7v+xa*)efnMms;VOYWHV zAcgeEfFD`Z(*+IR4>>v{UrybO4H7xlV>*+U?sN_fQG$Bl?Pqq=t+t#{C{$X=`oHeY0x zfRApmO!)RrYdoP7d0K5Re>>w(+j!~g{rFAV58}5vX!z?ejm+F~07J_KzO9XE6P8k> z2ku##|0)Yw*=TFuE7@$l6@{J>Z2{Yy=kxn1j{L{2X5B zNhIm2YiE#Z9*Eg1sSUrH) zs=Rc~>U#}Xe5o+!B%wrCl>?I#KQhFpJkHxQNk|+T=xe)Hu%DBOq&;k~$3c)S;7JWO zP2ftFC6Ea`e3s8E$&}M3Mi;aqJ>qitp)c4nHN-*sR>92* zX4vu0tAM>2Dm*6NKG4TTqB++Jdirb$9b>skWKfEPLY)zM^I%*xGs$d?>L*p+Fk^v` zVY+DZ8k{mgZzDth42MRFcuphU6+xIrsB5p?o;x>lYW}UB_PRe}ubOl=pXUuCbVD7J zTa>ZOxzZH!*i~X|o21#{QWbbCDj{*^4&mPmWa`9ESd$)kxF$UlUC_&_w8~!W?*#g+ zl~}O{+tGxnJe3lOSwc&q4LE#YrYEv%NFc(Fm~o;` z9K22R?Y#4`(c!!v&ST4Nx?lwvnJP(@Bh^%7Q80WiA^D5(h6#>_sR`xZX`P-EmusEa zIG;YL9`N&qwNucm_g}+|!Odt6$L7xE^I#$otPXD<3vQ@IBw|gkvJd*949Y~h{Vuw& z>knB7o+07?ERy(?`ew}HNi6K7#Zbn(==cykjo8MYlXvGFV+xp{YoX15v`%`KmGEtZ!Oz-n)>8GMw z!?~AQUr^uVIGI$W^&w)T1Nfp`Fu9Y3UUH z6ln>GejaE+_fdtoZ#c7>va#+uU3)k8hk&zQ|I+jFdJ?hnu?ewVVn6S1%^JWS+IFNy z`CF^Urr`3l@Y^c6uK;)Yu@7zb8#~q%pgeTg)R{Hqn%Wfzxx$`F9JrkB0YHVVg+E3 zZ-%`m)d>6;L7-`8v*N zoGQN*jaDRMue<4SZGERp7FQ@0`31VHMsa5fgG~({VE4r_HK)ltvE2WjQ==SrVtmSL z17sIaJM^{|XcFBB&#CjIT?h&YKhgdBM5Qs)9x?gD!?I~`%9q_!CX+z04S|}}UP*}1 zqYxgjAYMtLI;yFuNh(A}rbrvHyszma3J~Y!v%_6?EOa0T$7)}@V}7x;pV<=p!|3BN zOb{1kH8fgw6d9o^$^!?dM8Ecia21%H)3vp%#xnE>iGkDVK&R0rK`{<0pmFbnK z|GxR}e=XBLJOOxF=DEgzP%uL>`;!tA_qe-fS0!zH`c>JgDj=@XOx7kPRQj~b{VmBGd6h{5e0+t30rT#Emt8XidIsD*JcJvN4oYm^ z%xmHOYyBD#k+$~sOHh%_judmJrk(14CodxdeE>Ptdv3bXVBs%GiK+8&?BkEMHR7A& z*vT4(?_QJNqLxST>SX5_8UoO-B3}7 z$RMB*9Etre$x%|<)(fL&Y56*bZ_Ot;deZ(&E@T}`nlA%oYTmH-IxJb%$XB4Lqby`6 zMoc+K0KLWXSCjMmk6zJ;&gY-Nc(uk=4klM89#eo;2OwL)O;mgPU8hdT9F{$jI0$8L z|G)sU0^2WSV@SX33BU%FHehb-k_{KS4MjSqn(OyMRIj0A_8rcj;UWkhBxtpJ$1f6Q zTNvMH-FNLw)L-T0?fgE-_+&Gh#QCEsVb!iH?S25gGhkmsB$1bp4)lm^+cxj2Q5+(o zHZsvj?LaZ2wFF82IJpdHIBO8#qt+l$4Oz=ZfCh?Xwyo8M99u9b+OZ$zwV*l07DzsN~+VfoBF;z!^7nS30Z`I(x`|D2KBabE14v(@Wr5>L^ni)=lI=8utqfZ za&`8(3d1tusd#_*Uf|Vi_rnp^$3Bd_u?yUs$`D42n=a_@1%XE(l;7zZhN+1i+kGbs zZ}spf1u`-<)+9x>vXKZ$tw4Vu%7(v{%D&y`7+)pL9Q%Cqi_ddW8>oO$Xez71IgPWpT( zyj~U6Gl}yrd!T2qpdDQbMXFSKs)QQG@qq3Q*D=Q?Ap+xk_3CnxFI154!{I!uih7oi0Q}A>28Ypa zEN4OrCV%@o3|b^-5aP_KU%%Qw%Z8F|pz>TK^H5CJK>y79xpdDHUSiokZiatHMMra% zh!%Q?hwI~j_TznETAR)HOB+pSY@benHOSSgY)>7BIAE>;6w!-cL1n?64M)`* zzv*xD4N`9s0#kN&c3*#gW>8P*unD==wt@{tfOBD^((-4e0C1`kO1|VVU7(6GF=LMFp(VJKi*9iIJ{ygEqoSR-| z1>Ro4(;_JK)TDLRj)1>SP4hvpyw$c$@wv4@N*j*v5d;y%;A!K(fms;b)VNMc%Bn7v zlKSeEZj6Q(5mN{ZH1+i~+#0$d*3sUsX!C(jI$WP9v?m?yC=v5fGU>Ul(qca<#LC7N zd0wHxHGbw!L_}S`N5+E(h4wq0-8qLZESSuA>OsG>TCp%SWP6@BLX_2;tys)XoB;KV zt}Eun2bAfnvv#LScJC`9)h)BKwoW@+za@8mdO+3n`|1D=!Apn0IG$QilxCKec3vf}AiQQnb4cw}> z?4IbAELuE2FB0O4%Wkb+wIej38!V|*1y4)|kgD^aXzrVDe*W;G)_E5f{`{$Sa`6lF z)lG<@6gzUhFmUg)9GZ7%$=`XEotHczlx&`7Li{AKC!8P429?jU?yAf8@?Z?iLUPPc zTx`bKiSPAA07G+54v)8%-Ey!KGpz<0yRBTMRjJ7bxeR^no`4Hx6p>EX32JM;7Za%s zHPUm**2yVt&AHx}FMn++zU3K^!kpslle@0pnY@ZiU;?G9$!Hz9hjYuA%A4)24DqF-L)JDn>j_6Jk{9(K?m4pVV z^}<~1LR%y?|Fuf1Co>kY$!B4XcriFl6WA5jn*|C@fzO)x1tH4wzj$n)z(fPP(PB}x0Fj3&MT@g`pT^%Jk zi8Nhs-Up2@P(=6$HCvMDD+yUy{cFEdO%ERAfvSsRvgDeeKb+K7e(~aNshbZ|Q!DM=f#D_+nsfMvF`%+A zX^?qJ)-Bzq9E-a1xyj4rnsYx_vT}y*qaIIAitfiLzv7P`75&Ub+PRU>-|8obT$+8wH z5Be_3a=1-(47U^`?yBowPQR?0CJ|W7;%PrWKjcX|13<5Lz`1+Qq5ANU?Sk)XIe?b3WOk0)usOyCF)4x(PbzH;aLhK z8n|s?x}v;H%*=FT)ak7S#CKg?U5n6z-b4_Y9QgV93F+eAz4VbD&4uGoGD5iYWCTWg zGDt)#rzZwC67RF=`Q~zsd;lbi@R&)P<bXk^e zSK*M>)MQ?aFsp59w4(^o!%|WaUVY+gG`-C*N#}()`=wN?x63Jnoqb&4M9w;R%9r{4 z`e<<(imx^dh7mZqve)zyKAZL+phH(mf zc*uznCu`;8_K*Dh$p`aQnfJoc6!T4<7*SJ%Q;&}R@iIIWu?WiO*ktuYU=ZMz<5Lc( zq8^3tjW|@_H!oOaPSb*AVYii4Y}HIsPEJE=$YU|`CQfF{4Zntc7B&LXL?%`b+lh|N zBP<%J>KHqM**U^V1z`9*v?%bt6crVvXjw>V0{~s+qFG8U8OKoTQYw}Dx~0V&N`M!D z|5OPwK0i~dFM#@PImx7nxPj1#KFrB6p}||e^zrX9K+`Wrt&LA$lgYXUznAq<-#hjx zPi*t|HYpHXn;bzCb4OExRoR` zs^>hgigX;_WG*hQOAH3%JkXP5?B7k9XYihB&eg@uVA0jpC0IUxu4Byvr*BoyvyE%W zo33c{u)IFm98ARbCtZ<0{m-H5|Hn~W^O}4tM$+4Bpq;$i^3p#5L6(zb literal 0 HcmV?d00001 diff --git a/TB2J/abacus/splitSOC.py b/TB2J/abacus/splitSOC.py index 5521bcd..f8c9329 100644 --- a/TB2J/abacus/splitSOC.py +++ b/TB2J/abacus/splitSOC.py @@ -3,9 +3,11 @@ from TB2J.mathutils.rotate_spin import rotate_Matrix_from_z_to_axis from TB2J.kpoints import monkhorst_pack from TB2J.mathutils.fermi import fermi +from TB2J.mathutils.kR_convert import k_to_R, R_to_k from copy import deepcopy from scipy.spatial.transform import Rotation import matplotlib.pyplot as plt +from pathlib import Path # TODO List: # - [x] Add the class AbacusSplitSOCWrapper @@ -28,17 +30,41 @@ def __init__(self, *args, **kwargs): def HR(self): return self._HR + self.HR_soc - def rotate_Hxc(self, axis): + def rotate_HR_xc(self, axis): """ Rotate SOC part of Hamiltonian """ - # print("Before rotation:") - # print(self._HR[0][:2, :2]) for iR, R in enumerate(self.Rlist): self._HR[iR] = rotate_Matrix_from_z_to_axis(self._HR_copy[iR], axis) - # print("After rotation:") - # print(self._HR[0][:2, :2]) + def rotate_Hk_xc(self, axis): + """ + Rotate SOC part of Hamiltonian + """ + for ik in range(len(self._Hk)): + self._Hk[ik] = rotate_Matrix_from_z_to_axis(self._Hk_copy[ik], axis) + + def get_density_matrix(self, kpts, kweights): + evals, evecs = self.solve_all(kpts) + nkpt = len(kpts) + rho = np.zeros((nkpt, self.nbasis, self.nbasis), dtype=complex) + for ik, k in enumerate(kpts): + rho[ik] = ( + evecs[ik] + * fermi(evals[ik], self.efermi, width=0.05) + @ evecs[ik].T.conj() + * kweights[ik] + ) + print(np.trace(np.sum(rho, axis=0))) + return rho + + def rotate_DM(self, rho, axis): + """ + Rotate the density matrix + """ + for ik in range(len(rho)): + rho[ik] = rotate_Matrix_from_z_to_axis(rho[ik], axis) + return rho class RotateHam: @@ -51,11 +77,28 @@ def get_band_energy(self, dm=False): evals, evecs = self.model.solve_all(self.kpts) eband = np.sum( evals - * fermi(evals, self.model.efermi, width=0.01) + * fermi(evals, self.model.efermi, width=0.05) * self.kweights[:, np.newaxis] ) if dm: density_matrix = self.model.get_density_matrix(evecs) + return eband, density_matrix + else: + return eband + + def calc_ref(self): + # calculate the Hk_ref, Sk_ref, Hk_soc_ref, and rho_ref + self.rho_ref = self.model.get_density_matrix(self.kpts, self.kweights) + self.Hk_xc_ref = R_to_k(self.kpts, self.model.Rlist, self.model.HR) + self.Hk_soc_ref = R_to_k(self.kpts, self.model.Rlist, self.model.HR_soc) + + def get_band_energy_from_rho(self, axis): + eband = 0.0 + for ik, k in enumerate(self.kpts): + rho = rotate_Matrix_from_z_to_axis(self.rho_ref[ik], axis) + Hk_xc = rotate_Matrix_from_z_to_axis(self.Hk_xc_ref[ik], -axis) + Hk_soc = self.Hk_soc_ref[ik] + eband += np.trace(rho @ (Hk_xc + Hk_soc)) return eband def get_band_energy_vs_theta( @@ -65,14 +108,20 @@ def get_band_energy_vs_theta( initial_direction=(0, 0, 1), npoints=21, ): - thetas = np.linspace(*angle_range) es = [] + es2 = [] + # e,rho = self.model.get_band_energy(dm=True) + self.calc_ref() + thetas = np.linspace(*angle_range) for theta in thetas: axis = Rotation.from_euler(rotation_axis, theta).apply(initial_direction) - self.model.rotate_Hxc(axis) + self.model.rotate_HR_xc(axis) e = self.get_band_energy() + e2 = self.get_band_energy_from_rho(axis) es.append(e) - return thetas, es + es2.append(e2) + print(f"{e=}, {e2=}") + return thetas, es, es2 def get_model_energy(model, kmesh, gamma=True): @@ -102,16 +151,17 @@ def parse(self): # print(HR[0]) HR_soc = HR2 - HR model = AbacusSplitSOCWrapper(HR, SR, Rlist, nbasis, nspin=2, HR_soc=HR_soc) - model.efermi = self.parser_nosoc.efermi + model.efermi = self.parser_soc.efermi model.basis = self.parser_nosoc.basis model.atoms = self.parser_nosoc.atoms return model def test_AbacusSplitSOCWrapper(): - path = "/Users/hexu/projects/2D_Fe/2D_Fe" - outpath_nosoc = f"{path}/Fe_SOC0/OUT.ABACUS" - outpath_soc = f"{path}/Fe_SOC1/OUT.ABACUS" + # path = Path("~/projects/2D_Fe").expanduser() + path = Path("/home/hexu/projects/TB2Jflows/examples/2D_Fe") + outpath_nosoc = f"{path}/Fe_soc0/OUT.ABACUS" + outpath_soc = f"{path}/Fe_soc1_nscf/OUT.ABACUS" parser = AbacusSplitSOCParser( outpath_nosoc=outpath_nosoc, outpath_soc=outpath_soc, binary=False ) @@ -120,24 +170,25 @@ def test_AbacusSplitSOCWrapper(): e_z = get_model_energy(model, kmesh=kmesh, gamma=True) print(e_z) - model.rotate_Hxc([0, 0, 1]) + model.rotate_HR_xc([0, 0, 1]) e_z = get_model_energy(model, kmesh=kmesh, gamma=True) print(e_z) - model.rotate_Hxc([1, 0, 0]) + model.rotate_HR_xc([1, 0, 0]) e_x = get_model_energy(model, kmesh=kmesh, gamma=True) print(e_x) r = RotateHam(model, kmesh) # thetas, es = r.get_band_energy_vs_theta(angle_range=(0, np.pi*2), rotation_axis="z", initial_direction=(1,0,0), npoints=21) - thetas, es = r.get_band_energy_vs_theta( + thetas, es, es2 = r.get_band_energy_vs_theta( angle_range=(0, np.pi * 2), rotation_axis="y", initial_direction=(0, 0, 1), npoints=21, ) - plt.plot(thetas / np.pi, es, marker="o") + # plt.plot(thetas / np.pi, es, marker="o") + plt.plot(thetas / np.pi, es2, marker=".") plt.savefig("E_along_z_x_z.png") plt.show() diff --git a/TB2J/mathutils/kR_convert.py b/TB2J/mathutils/kR_convert.py new file mode 100644 index 0000000..06aa8d0 --- /dev/null +++ b/TB2J/mathutils/kR_convert.py @@ -0,0 +1,90 @@ +import numpy as np + + +def HR_to_k(HR, Rlist, kpts): + # Hk[k,:,:] = sum_R (H[R] exp(i2pi k.R)) + phase = np.exp(2.0j * np.pi * np.tensordot(kpts, Rlist, axes=([1], [1]))) + Hk = np.einsum("rlm, kr -> klm", HR, phase) + return Hk + + +def Hk_to_R(Hk, Rlist, kpts, kweights): + phase = np.exp(-2.0j * np.pi * np.tensordot(kpts, Rlist, axes=([1], [1]))) + HR = np.einsum("klm, kr, k->rlm", Hk, phase, kweights) + return HR + + +def k_to_R(kpts, Rlist, Mk, kweights=None): + """ + Transform k-space wavefunctions to real space. + params: + kpts: k-points + Rlist: list of R vectors + Mk: matrix of shape [nkpt, n1, n2] in k-space. + + return: + MR: matrix of shape [nR, n1, n2], the matrix in R-space. + + """ + nkpt, n1, n2 = Mk.shape + if kweights is None: + kweights = np.ones(nkpt, dtype=float) / nkpt + phase = np.exp(-2.0j * np.pi * np.tensordot(kpts, Rlist, axes=([1], [1]))) + MR = np.einsum("klm, kr, k -> rlm", Mk, phase, kweights) + return MR + + # nkpt, n1, n2 = Mk.shape + # nR = Rlist.shape[0] + # MR = np.zeros((nR, n1, n2), dtype=complex) + # if kweights is None: + # kweights = np.ones(nkpt, dtype=float)/nkpt + # for iR, R in enumerate(Rlist): + # for ik in range(nkpt): + # MR[iR] += Mk[ik] * np.exp(-2.0j*np.pi * np.dot(kpts[ik], R)) * kweights[ik] + # return MR + + +def R_to_k(kpts, Rlist, MR): + """ + Transform real-space wavefunctions to k-space. + params: + kpts: k-points + Rlist: list of R vectors + MR: matrix of shape [nR, n1, n2] in R-space. + + return: + Mk: matrix of shape [nkpt, n1, n2], the matrix in k-space. + + """ + phase = np.exp(2.0 * np.pi * 1j * np.tensordot(kpts, Rlist, axes=([1], [1]))) + Mk = np.einsum("rlm, kr -> klm", MR, phase) + + nkpt, n1, n2 = Mk.shape + nR = Rlist.shape[0] + Mk = np.zeros((nkpt, n1, n2), dtype=complex) + for iR, R in enumerate(Rlist): + for ik in range(nkpt): + Mk[ik] += MR[iR] * np.exp(2.0 * np.pi * 1j * np.dot(kpts[ik], R)) + return Mk + + +def R_to_onek(kpt, Rlist, MR): + """ + Transform real-space wavefunctions to k-space. + params: + kpt: k-point + Rlist: list of R vectors + MR: matrix of shape [nR, n1, n2] in R-space. + + return: + Mk: matrix of shape [n1, n2], the matrix in k-space. + + """ + phase = np.exp(2.0j * np.pi * np.dot(Rlist, kpt)) + Mk = np.einsum("rlm, r -> lm", MR, phase) + return Mk + # n1, n2 = MR.shape[1:] + # Mk = np.zeros((n1, n2), dtype=complex) + # for iR, R in enumerate(Rlist): + # Mk += MR[iR] * np.exp(2.0j*np.pi * np.dot(kpt, R)) + # return Mk From c1df33711951735b9dcb348c2f7e139b42bc2589 Mon Sep 17 00:00:00 2001 From: mailhexu Date: Thu, 23 May 2024 20:02:49 +0200 Subject: [PATCH 4/7] update --- TB2J/abacus/splitSOC.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/TB2J/abacus/splitSOC.py b/TB2J/abacus/splitSOC.py index f8c9329..cf4e6a0 100644 --- a/TB2J/abacus/splitSOC.py +++ b/TB2J/abacus/splitSOC.py @@ -47,6 +47,12 @@ def rotate_Hk_xc(self, axis): def get_density_matrix(self, kpts, kweights): evals, evecs = self.solve_all(kpts) nkpt = len(kpts) + rho = np.einsum( + "kb,kib,kjb->kij", + fermi(evals, self.efermi, width=0.05), + evecs, + evecs.conj(), + ) rho = np.zeros((nkpt, self.nbasis, self.nbasis), dtype=complex) for ik, k in enumerate(kpts): rho[ik] = ( @@ -96,7 +102,7 @@ def get_band_energy_from_rho(self, axis): eband = 0.0 for ik, k in enumerate(self.kpts): rho = rotate_Matrix_from_z_to_axis(self.rho_ref[ik], axis) - Hk_xc = rotate_Matrix_from_z_to_axis(self.Hk_xc_ref[ik], -axis) + Hk_xc = rotate_Matrix_from_z_to_axis(self.Hk_xc_ref[ik], axis) Hk_soc = self.Hk_soc_ref[ik] eband += np.trace(rho @ (Hk_xc + Hk_soc)) return eband From fc96517fd74f54f373db32c8592307f9c18553fc Mon Sep 17 00:00:00 2001 From: hexu Date: Fri, 24 May 2024 10:58:56 +0200 Subject: [PATCH 5/7] update --- TB2J/abacus/.gitignore | 1 + TB2J/abacus/E_along_z_x_z.png | Bin 28080 -> 0 bytes TB2J/abacus/abacus_wrapper.py | 4 +- TB2J/abacus/splitSOC.py | 127 +++++++++++++++++++++-------- TB2J/abacus/test_density_matrix.py | 38 +++++++++ TB2J/mathutils/kR_convert.py | 12 +-- 6 files changed, 140 insertions(+), 42 deletions(-) delete mode 100644 TB2J/abacus/E_along_z_x_z.png create mode 100644 TB2J/abacus/test_density_matrix.py diff --git a/TB2J/abacus/.gitignore b/TB2J/abacus/.gitignore index 75caa95..5efca39 100644 --- a/TB2J/abacus/.gitignore +++ b/TB2J/abacus/.gitignore @@ -1 +1,2 @@ TB2J_results/ +*.png diff --git a/TB2J/abacus/E_along_z_x_z.png b/TB2J/abacus/E_along_z_x_z.png deleted file mode 100644 index 9e41eb571669f133bb2686344cb7b9e51231741b..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 28080 zcmd?RXH*s2wk=%boO6Z+f}((+ph$)#prQyWIg647$w_ipfPjdIh=PC$NX|J+P!vQ# zP%=nHg5>BM<=*Gq`|dk0eD9rqpRMh-b=6$8s%9E}^xnsexOzp6ikzJsK@h47>R4?A z!4E(XyfIQ@_(Za2=m-2y+V#Ai>ov#Qt{&#jw~))`u1@xjuJ$$-yzaN0U2GihoDw}F zdh!IXwX3U>%Q-PIhkrH@b#%58lV!=Lg`1E$sq4ES2#q=VA6}M1rVWBPa9zMEUH5$U zbJ$Don$@AyiVWG_V^TtV@*s{ele9 z8iPB17t0LttfR^dF4>(gwP3{Naf}gTY63#I@THe|)NL*A@;?h1N?my0H~m9aX!Pb| zm$e_yX5vVLcmr_oCpYZ+B0l<;EEjQvo{o-ApN&`%1AmUUAnI5w*76oB5`ng(B;bXs zaS7 -AWpbUVkSe-Du2XM9H{mK77f4T7rGc7F=fwgcoEVI2)CE1xfkjlW7Xsa|k+ChJA7 z5=Pk+O2P7}#={Z)@^Fnud)IRraj#!Lqz1WUJ!y~Y=D+R9P#d3{ix~ZVxP^;}2`eda z$yvyJ`O@r1M_i@vj>P`K!E&_;%pnP}eKI}59G$$N@Nkl)nq~Rgy`>s!YwOCNWtb5BjySUKN(+6?Od5^8O^Ry=%*QNX^YEp}{`SO|%gTWWl zOu^dQpP=TJ#T~iyxOB(JAdi_of{;9PB7Z9ZCg}$nA4KLG36F1pc}rN+@Gw_6HCOZd zgyVE2eYm&W18<644*XTtZ=zCQQnxZ5&Z8`v zIO*uW$KmYk?DXr$g@JG1%*UHUAH1{NUjL#=H7M)3!Yg>?8Sc#+;j?GYvdWn{IEb{x za7B%cnc06Y&7SPOiX+V3OvhoTinnCagQDml8}(pZhSTL-00I?os3nOfUze;jOx6ZgcHbX zJeJ#fv#*I;yd^~kth~HD%z->807m|qoDVY}A79wix#!fdlsn$NQ?jtY3=It}t*3^d z{lhATKZkTCinRw3QG7Sf)GaW{6XlV#Wn^b(e~^$6`r^e4slm=E=Z-1o^`C<*8EUbr z6?cXU@9nwEuD)ehSzX0L)|W>|3+uJBLozdkic3oHkc5N;^mSG;6%Mjn%N=DKG&8@9 zqk}iqnHcxrLF0-V7TZfCj<_SQw8>s+d+9o8K0Aw-eVtrPtYGjCKPf(;lA^IO7h<+J zP=o~BzfZY7+e3g{H86;rn6Ox1{KkNsv>!N$zS+IqbsLwOkK@bxJ8J_Kjy$rmvP-t9 z=-yPjDl~TGWskMfj}AJGr>DvC^YcGhpVq5(J;v=l%jdh=%7O$84H;>q$g*)r**A8? z@!y^MeB^Cg8``fS+S%AESEzOJ3_@SN6nOdS6Qvn5x*#A$IreDQl3#mLA= zSoQGmu(_3$(EA297!^0!O~f7JSFs2cgTk()}lUZC#zXc{3sBIoFP?ey5f&D)b$1e|Q4@85zMKgprpY zaq1c!YYHy%|Ls@1JD% z7fCEJuH*A2AU)&3Cw)SG8yBcyEy@1gB(vc``4T(R`fSIirv=@mv99582i6*4yj_lZ zq9X|R;xd(*Op>|&<$L^ovWC{ZI)ZUyW04nD;$+;l%dov#p8L0HcT z^o32{S?%;Z(X&eJusWO2fgQUPt!)AX&zf%2;LA}`&%JQh54t2~S*iFE*)RKM_%RY{ z=_HMxKWp!=#`?d7-Iy=oA`Dm{jlw>ysHmtzd3Ii2>#wQKqFbHJy?uS=^S#;IKjJmV ze*UCF^cwu+F#P<64u*Pqq+()Xh$8F=k3|f_PFZy|c6ZbG?+uO?nl%MckE&tl+oKf* zFCs%Tqg08P-dxqbF?yRLf^d8CrM#0%spN^{1wG>pWGqf=-A6`KI$A1(evOUWhhJzi zaVm_naB-IM7Zk*l2<>V%(#WmT{#sjGoy$yB&iukBAwkctdcSe{WAdQ=S~G>%7zD2N z-Hiq4(;~}v^y;s)1&$oSLu_qrU-a2Ve?OsDh?JF;scUEuB9^yrTZ}aZdcvMrL z+g=$*0=Bn3HB#ldq@|@-HZ~l3RCD^2v#(n~P<|q5$BgjG$#KD63;PVKoVRO~6zrwv zx0Z(K7#Mb_N9p zk1sFNM5}(WO{U_IYJB(Z9W85e>Ry`z(g=}caEO>rKg_htw6&Cmrj-5AlAi1~<)%2f zKQhMe5mH^zF5)w9VA!oA*!pB7#q;O*VWl7Yi@1h+KI9G08oWJZRIrwqbmGRKApDjc z&yV)@_G3#+F+o9({V$xUC@LulR!%!welLXXrD79Tdz#II(48WJYk5=_$lZ!Z5D!eTw%3$O&{>J{NTHXOLu zEWpm5G?8rx1EzlIQcGJj$Fi4UPfw5M+Q&08FN87kr0DWwyv{qPEH1@ef8XFaw?sCf zW)8v1=STHUE`_bd#}4i+%nZlC{-7YZevLh)`&}3D@53EVDWl(Bo~EjF%%h_uG&HN{ zQ_911$r<9THRzE*_r1v-OP!}Ko-i!`;M7oC%jxIm z_j5^zE=kq?qLj|t=grYp;<4ne?F!s%INpRD^Q73JqH8K&KXz~Jwh+J*^O%g@C03KT zw0fB0;_2TljnryC!M(lt5XrhOoOXR=_oY>tL3DAk#LL&O+3^jOmGJ?LsYbK+*(QJb zbamkScMFIh5NQ#_($bPwMux+-?7Q{w`}Y^^hpU;O6Z7O@>eoa+{`svWvsLPjrK2Ma zU>n0F7P=%I6Mab?UzhjRddv}odU8IbzD!rr`}!W6C2*9>%CK`#o3<^;bg@2OAR0@k zy!Pai=DQR{%>4EA6arR2f&l!?d(zLVU%8^Htxf;*L4sPEIEu#>e9! z&Q+XRT3UML_H3I!2g7@NuVYnJmGg8eD)G-sNijozd;WSSWGK6})~yhGDN%I1!2jU* z=PzGY-wCCfEqu*`5Sl&b@eH?QkYUouG5Xv=TriavNxbyIBgP~G8&q#1(tgH)G2oJJ z4mJz2q@*l7T13j>dkCOC6i|H0)NHT3J6U-9BYVpv;4U^cwygZapFigMa&&J_e@KE} zYYZV{c%hw5%fixcF(@uBPD(@Lc(;iZAtWZAK=;BFIiHh&)oXtr?&;OI|8$7OQ^!6p zcwJA|?Vl8xYLJwzA{fWp_+{Vw_Z`OBxuN-l=G3JI+|%8*E-`qZ$Fl@Vh7vLS*j}%T z!V9z;ty&OQh-EmAaZ~&Z6nUq6M>?RSs>(~z23|(BE?{q;NpG*$$+L3hGy~?^%TGn>1q$^I%yB%?) zE3V8n*mTk8c1H-_&1Y!i=D0fR!$M6t&)eLP+aKa8l6MyHK0cUoe2Cuo=$p4eZtps( zbO=1GWODC$pAFmbV5m8fE*jPDcSn)Vq!ww=>@4<6V;CX{oyRqEg9oJ+|{F zx;CY@+xZFZKbDMU@_Z^6i0w^-)VL$}?W<{Ms40f%w|*6mT}rSM2{y^#)%$qJpfzI# z$#ThEVYR5=f=1RwL-cLi3r8KIne_S_@8U;fO+wNfrw=b4ZB?(XZ4`BhH1DtA^|UmI+xqpFpZp zWihyBd;P1V-VauPu2){UiZm=T`*S&y>S9@n3XUwk2hFUN(gTsk+;Qh&ROJG&0>!(; zaw(%vBH`;R?j`Nj9du4^g`u%Sbn_$JI&6ms>t;i0#c_GlllUWdHqYwnMX=OsoWB=r zviJ+v(F|NcoHe0RW2EC1WtVKiE_Vsvdb5x|;`+^PVb0Z^AUdb=c)mGhGM%>P`2Mwg z={`*W)6(@1*1X>0`@!5+6nyD2t?6!xFteB zqi#=!Gk+PLD0vN+F4~oZsw0Oo^gRwY5rOn8aZj{D^3+_AMq8_$*`s+MnK1}-<1rI0 ztfkgjh7P$g^$Bc6W@Rz;m%{^!$m0turPV?y4Gk4Hx_|X(J+-3PIk!PxFa4O}RgT<* zDJKmLi%-!dehE8SzcpbqjD)xn?Mi=vW;%RdX(V4yT8@yqT!OfP8V$-Yh+8|yog~zv zU~sn=`h|)%=f)ololZPS*Au-wB~@rC z__;=8-Q4-@Un{qYE@3dI+Myc+UGfX@e=qOqxgCE^#IJYjxk8z17%k74+Wo61kum7f zabNS$4ko%dX?S-`d{&?2gr+H%nsG2U!`?AEX(`3Uht;(m;&`xyM^qEb9hOcJ&a<<( zibmi*FhV~jy_H#IgSi?IUcxovV^7+BVS!LPn{Zi9igl^WUAVzIX)fp9u! zWT8qYulb%3sofheJ8JB8!`5om*{!jrgqQ%apuRoPSe1j(#bT#%Xk(4~a^ls@wW8CT zx=7=f`4msy0mC0}a=U-fI)308Wn0(w+NvTbRn|9JKYW(h2LM{_HsfE)tFVY=WSYO5(bNL_5Bv6?!H*sIfm32V2YR2$A=>3O&Jib=t8x0cFYF6Q90(JGtz9s$LXJ2wR61N#)h+ zHm;_SPvorxK3lqM35f$XoJN8?FzzP7b~P=X7(T4LYdLmr?PR(T1l{ph%;n%_Zy%To z-i}>q5AnNLTQo@tX$dZ59kYMSfOgU0>=B+g$(xXo499D3nJVAD%T#OhbM<@6LE6&K zrWLE(zCU6-w^sSCR)3$=zRM&(?0;eLAx_;A=s?&0;>etpyAwo6R?=l%BLC85dWj8ORdg#VBcL zXeui!d4bgl4Gj%z?=LW=iNJoI&&$dpMIcg!AG`KiG5dNmeA%{)UEy#$>E_x_m>wUl7jEO6^*zvKHenH+TYFqpVpR7QW7dvy76P?A7W_&0Iatpow6#l z>gGn1aT03I#*%LN!8^%lvOUxFe4yNcV{>^l7V_1Q@Nf?9R5>3h$XufDT+6)D`Zj`L zd}@jaIqNiryw*qx0xFDOQj#fN@N#ntjl%uO_Sj@s|A^|yWzDI{#I@GKlDG+;Sjrzp z4a|y|xH0%^0Iafunf!_HFGT3X#>A0g4#}Q8ejIjqu*ZMo2vJKo4GHqLz5V9y`W$dI zfdvJQQBpu@U3^8Ke8#a2SOf%_o0}sgBMVf`Xae@mu(Zm1L%7;wIR*vNsR|7ry@@GV z1HP45O#*vWlegRXpE6Zz59Qm(Q`FOAeW`schEw{pZ8Bi0S=VzRp`mP#iAYFV zMn`$rqhI$)0hbaGA0N+mvRt*GKO;RIe`#sSqv{Oy>twLp%+>2!;Rub8tU^)6z|<$M z-Ijy8o-ZwAp;K;pbu+aCQ?zlUSgHUbBBHQV%j9nZ1C2n}5fSD6KCN~2Dk38zbCi)h zV+z<4hw_Msh^C<`m+{$w&4dV)BNIM;e06&=wxa3^5J!-lN8iDIop#MW~kdi@~e+ z5@1|CY%COxK)w`cf*fxCw3&5G>5RoypN^*_tn=M*hHUo5u$%Xv1mISsaD~2#3IPzh zY1-Ka4w>57S4~Z00d@oOghWkeFxR+-StXqMJaD8rgM-DlD40c!a1MjzYzMhRPk{~n zLlEsj!bBNpDG4OR_G06aeeYM4qMK4-eCIoy>!eB)Hm!_WR~IFGy~pUow%^gs&V9Jc z8Ic=s0Usb}h_=^43f})4)W5TkQ&8W0s#`W##IW3~?LJ#YZvV&#_vYqig;8u)PQ>MR zZSO0uSuY7ZWy}q!NDi(%-y^r??$QYP(@Lzd)vCMimjVfGGFkvWC2HGIMqG4B53^>F zu#gb5@1-_+rPg$S>VX3d0gBhG>zQ=e@c~`UBAwVG8jN#@xgkG zi;KIO(HyyAB3YcBA*2C#2xww!KFg85x1yL<~f_L;0$;>dV`X z?&CFX39ZEyy~_k@*QM4~f&4*Oy7DsoQ?Mixw2H# z47J_YTyfd_FLMv#;sSwhZ7KfWgM3)U#E^sN)q0Bi`ub8RVsm+qb&yKn!Tla-@riUGwA7c=Z*XZN;!)(Pp(J-0X_dgOQVQok2x- zzV{-k&&C@ODr1dRV7k5Gc`8cnY4MWB@61aX;2Kz|tDk)1-xp6yrT;zT#ayGwGWMyJ z9!RQ=m(NCESS=u4l-ca98ogw>-uSVJSWJ^5th$1v#gt`=-{bA?y#wCLa|6604ViQF;s59+Kv6%+qJ4M~`p$4~?~>qr+jpEzqC|d!!p=z?F18h>3=X0Fkm5oxFIABvf$T z7EJ~^d=zde_8&F#HPyVo?6V+Umu>g-Q$NGxu4j6^R_s>3DJ2$bOagd^?b(N8EglWi zN|?Yj5?pRP_Lrg9SPX;L*@PCa&k@!1@T>R;X;46gdnO*@r%)T z2_BIY_>a{su1o}pyNUlPBRRf$Q^Zv3%Tu8tXO`ayKkHj0IiBTG3h>yh`MDfY(9Ey z9^_|IzHJm3+boXr@k>^#_DarDZEfJMTpY5g*TkJ*CE$%?>O7GjbohF0mxBmTzm;<< z+Am(;C^T8@_t;)YM=|>A;-*zMqtMxnmm*8@z>i^e?tXGRLE^>6EwVTwfSSb?dlHfN zmpTRLYnkL}OT^jrokIe!KfJRRynQQq07;^-AyYxKY2XGK?(#e4niCjX0s>ww zmV|=5CI~CEcl0)QV@@q5@x@o_7{^9_CvLe{ht$%0JY7D6IR)8V@A+HIx=LcKVs?y# zmQ1=)*J#~rl3H(6&i!oL;)r8okV76nerz#P>jk8%KKpz?y7%-dtl*U!(gH2!c@USh?K)IQcDd%(_ zrr9~lxyv)G>PN=*qB}GqxRIHU$>G4Wn{}ti4V2ljP;;IQuu8pGJ}d<^R4efKMGk{f zi=z!ZpalabbVgPlG@2hXGaSROuRKTfoMymt)`quk-SV8z&OGV$%QAkECr;0=E?u%( z?5}+u7=-h0aMcn#Rs88KBy6FLLaMuK2db8@LAw|lvKzJe{BI<+Ipn0Qx;Rl~P$G$P z)NfDZ3MgWO8;OaDDx!%hIeh@ue^PHg>nH`?tYa+HMY-kpzn;8A0G$ zK70~bd9#YO*EnYR@%>P2b7yQvE%4aRG%r*$0)Z+HQ9wv3D1xsGx5n|S892m0eMzuyM&%gNbUg0Ss2aCj|%=T~+Y3OpJ|Xp>@tSt9SAtKf_& zmpI||Laet=j`c9>8uZj^1>3>6+40tFBN1xj6X;^L4b8$U1YFRnSu7cX^2BL7X z!{XxNi>pWPKL(Mm;+vqLAhF9#S7^qhBT3`hO!rODbhm$XikaWM>0-F?wXZL&wsvS% z)iUPL(f>;>4`7<6;cr)d*miuclo>SdBK1d3eo$1wPPYo+Il)&QOhS5mFKKF?X}wnQ zqR09Hqyt%5M_Y>9+Eh_$+_8Q~5VR4G!NEvS){R8Z4DIjlYiLRI$Vn$%3wud%q%iJc zp_A{6fvE=xk?bs_5T?-A6;|_EAjbA)-zivFpcXN$=cHghVeRqUq|RIEm3A56&a^bV zSFc`O&Ay&jel7cYRzZQ9zJ8LHn3=h`d5*)^eB-Dzm2k}I)4hecfaBF-xm#B!+fj7} z2{0)L=OryHDr(YHI0zx(_I%M=u9I&RajWw&+Y$`@DZ+=}Y?co8hG|XLXWNdSyuj(V zbY&a${#+0!mPzb`>a0O47cOQj~BTu=wQ{TDRm&9w$)R?=IZRVgvhqb%? ztqdOwTK~$?Eol3YbPCl0$ruEoXZb`$^IOYrT3FC>bB{~}(WU2SWz{l|zMRHKOe#}h zJLXPl7-2ZLY)c*a-tK&U81CsmLFScofAyBkJu!37?mOGNq2eIy627p}c(4CEQwc;< zPd^LCT+gF=Km&O0hq0|A^gVT7zRniTIY}X#hjCtU>H=XVjhDm4GHP{%oX`KEF5JkH zV~kX9858V!)tG5N=+EPzbxVKOWHP%j6%<)b+`=P_6%d7RM4<~YdxoS0IsM%6NP_cs z5eYHv?6*!gc4sam0mT#vTmlLbDjj2{UH}m`%oOA4souF$67ivFIc&Qb;&l#_8&2l& zaezNMu@FpyBSAzH-|z3=~OHnlnW`V`}z*d!-x`4s@u zK;q%x{<3%AvLU7jjiN6;Ob8FjhWKRg^5#@lZxCZ|^OjVSpGn78TcIe7Wa$*N&doCflw*LXU_3((aP(r1@!3>_*6WNO`I<)S2jpp8`7w zLN;&w{$qRTCPz|$#!92jIue@iWtDgxayZ2Q>N(HHHB z;bK7S#`S@Lte6)NR!fufc%=zOSH*QpfS!ywhJp~W?=_lUcfvlukB$_&Ljc+Ds(aRA z9&s${V#M+01bMc@01d%*q8Uv@RzZ5L#&|MC#dQ^oeKD8j(ad z@e#z1=C*RBmpw@2Y}`t?6K@aKbnm8R(Ti?EtLm6<{0C$9LM@@`7Z(*gZY(OS47Zh_MHak;iJ=0jV$qLQv1Awh?%jzh? zdvGzpw2Otm?u8M0FMPEBt(VxmtRK*LyH4hk^_Fc{$(VmvWC|Gdu8!y^g%SK-{sx!7xfGCf~#ez#Cn=E{$>taP0pwzz6}SOTXiR#$G#vJ-`DE{5&x*8on8yKuU$ zfho(eck&VyNtmF7B-%>nn*U^D_u@fDex8@pT6DPd$v{L^!d9CFpt69#3}ciYcG!((W}V;O;|?S#sL2{!#rm`b@`K1#d9<%u?*dlpIR+`*x6z8BeND@I!K%XX@<2#7|-h0Q{uO7 z9_0*GnDnj<$nw+$8hdgjJ{&BP{#im7%Li`JxgA3vM8MQAIN@$4?~v2{mzq_!MLfP7 zO#H?B@&tcztO_JpE)+c92TR5H^xDXu{GxAy5i1b%Vd4#Q=e9KOyxyPtfZyVb>B4c8 zWPLX`+Tl&;eQt6(YGmgh4c-oaqwYFmw8HlvkSaGa%Vp;;apob1@))S9joo-!BonVf zv9~R^vVB2@X#`FBcZwWiU7eBVC;pt%pbGca$*$DIKD?@vacZ{(&_Bu+e+c0v)k}Ot zgy@|VL?(BI*$zA2?D(9WOj@|cn}6wj2(Ocwb$SfO8#y5&(pu}aHa0y?9I6}$kl?zG zPWj;L@^a=V7O^n!&lwpR(ZXim+^pO(xzJw_j>Xt)=v=>UsaKGnPhC*IE9||w2#zuz zLPEkZkUGzNtEVrU(c6@K_MxY_OEnUJVHC_P{U9{XmM>QM3V^Q%imE;Co`2g|AUwwBdQOt_m{TBu?H(Q@qU?5I{i z7`B*r8?o~1mqhLska9Osp*SK!NmZ4&es}&jArVoMy=lbvwVjQT1~u5QTUt$qEN~w} zGNa`D&h)aN+}z{JR$4bcdZoAW8>OYE!C=O7?AS3@dEc`Lue>}r=&NdA&+@%@&#<(| zwhY{s1GQeKZhw4Ags5w3s;R4!o@S>4JJxXeM<8Q*>U=!jzI|(UxW5xN1I#bAsBtx+ zhsRdU8&ot0r5NbdJ9qAQlu=<`*z7mlbzDZtkR^p^HGk9kRNT|&U40hQ{xr?1!c9$y zJOLM%c-=2pj(hGiSCL%J=uu3&$_m6k5TQ!K2@_x=%PYTd`Eux$Ocg!7C+7MfNQ3n; zt;gEHp~Sp}_^h-v8$LhcyEW{AzRif&443Qj2-m^xMoa4NeXHJHrx#lDzF_R=$znMv z{OS34?5rR`O{cg-BNvr79tVvT>;zhgyz{qi086T>@D4v9I+|)|co+!MPq#V#HAT8T zNVbtEYTVouOfo(;7VvYhBBQdB)qlH%=c=}LklVu7vg+qxvjCUSAMiHXb%(wz2Xfhh zhsm$2A~`vZI!j`wk2%{T4@)i8g`y$c(|G={6Zn9)W4c^7D^uKYN&3AhH=x_w!hUMee>szWGq5Szb z7FLEMWFZ(7dkT^TfLBq?z`&rU;x8o2W}HBQ{D*_WZ4m5`8yg!(j~_o}V4|ltV4GZK z*KYy&87UQ&q`}hUWH8tgAUWZ992yd$NK5$R$B)w8A+WH9=;TNafi=xp<~JwSq5AG@ zVrKh#$l~-I4UQBs=B9}>&T!I)&88&%5HkC>0N09vu7RSPXwnC3yScTM2yt2( zk_P`1SROO_X38rn!t3M;x?aX)R9z2R4M4$R8h<#%8(8{eE6^sP&+K;b|)gIJ9k9E`V;r;871Muw~hnRkt!V78Ml%-`p6lDgSMv3PceKkY>O_F|qiq1eHDxe{fKssGVy5YJ8|{~w@NV+Xx*?+x1z z&!iauy*@cX)a)Cu2D|c<+c~nK8~92x8X#P1)#J-W~s27~A-)NyOkB_#(C@Ol!I< zel}^Kl2No)Ug|GJtjNb$=^%*ZH+!}fYtLmgKsrSv*g4v$dqZYR1parrj;^5A`CPuo0+aOuv%RalRH5gDi zS%uE&H~;|W^q0zP_r+E*$e4|rrn>YxMGsfp#Uq}(-eP0l30UG>>sfaUx+bu$LTF9_*9q105N^g+vqRRIxMuN<{quJ*rB)AUfwjby!#Xj*aS;+C2KX`yOe$zi6FpO=W| zWu~5=phjv5eOe&$y?c}bp_Qy>fux=75O+>->&_L=LtwN5#Hep6`L9KkV|#_s_g?7? zwJ;LS4~nXbHka_h%(K#^$&Fb!1_NE`w4z)DVNIV9ZTD$tZ);%XH31Wy8zqswUs2Omu%?U7qsv}VS(WaOj7b%@^QRTn(@t7KNMh&% zQ@UVWdyi~rSDS(8Lzp&qO6{*ZrNQK|m+>Y-C|wfi<7`e;9=L?!QIeL-GV%IN$K8!IQ_ocrtxVfP zHcv-yzWhfzQP4gKW{*T5)l>NClb{aNC5WPP`^K(k6-t7W;uo%-L?4L?SUQrNe87vegd&_q9k+MT)2493GO@H@Di7?R%F3Q|#=lW=Odj zPkL5)U+D%<$8P_-%m6G$n#!4 za2Wg-1`h|B6pot5k`bS$u0Td15sRRVcE@M3pn#2V2|PsX%2HHMo$MbTjWUkyiCi-e z!|t)w)q0BUUrth8KOU>l5~>Surbpy*0O(wbo3X1V1d( z{adz$=IX2<8G!GDo{NCDN6x_G-q;MzTyAXI>^N+sy~0%)A{gZwf&lYEcZ?mFuU%ev50*cBfE3-Wk&W{NWK!Js(? zgXYQy@C~A>#N45Fzn1YR$REd%07wESm{qV9uc}y9a(>uf?>~csC4$4&k)H(#8SK64 z7cMldjJMz+zYmuEQ500WHRAIR|1bvP8>;LR(8wYPTK=M}tn9SCVt$`njvXqAoLu*R z3hV>1gW==D1D{~nEpYbY;o)t6Nj)S)NGU1bx^&BhSB&}%KC00$?%15HI$GeMs7WI< zG}?U$GSXQo4f4&Ph8l+AOF64QS7kqNO^E^XmQGvte+6F6zvk+P-2(T8ZoW|y;3z&) zd~ynk)s22r)Ud!+>d~8kPR|Sa$Ww2|muF(7E@5w95KI_5=JsW$dV%JaGHr-qDewH# zejxS#61?7)w{DH^%%nQHpsrwmOJT{Wsj1nWLZ`}>yq8CF zufk4bY4and9LM@~J$&#SHJOUedx055^uLg;{{^D{uT8$Hof#P!;_i!j`%_~6cnBr8 zENfR+7h1G&xKZHW)ZU)?U&&U3w8e+EiAEWg(?D76GYokiMP_D#$%3FDJ&MX(=syl&b#BDIsml= z;uN9hS#yg9!!2iqxSgkrKSEQuufURRgq4fph-?A;g1`nluyeq4`(#@Nk`zM?C+hn< zd20jD2wbyvr}ie3?gpanp{=9SW`c)R3nzBHtHv^4G4RGYsSDXPPFA3NtbXLO-$#_^ zI{W=1Kuj%`)qxdmr{ef=U#%nZdsFm4KEgeEk3n?1qXHPyWXKsGxtjB3GW^=k1Tz-P z)%>xEIoPKI4{`i?rTuU1txEV07gi0LJQcV+C5Bl3A~T7we;f$!l3s=Fc?YIC7E=Rg zsd%7oM09Z3>b(>>KA6M|9Tl<6!Ej9nz4pS@MtU;?*1@h##B*tU1ugjQ-B5TF4L zc(6YrZx+R`IGKs=Yr(u8rRUvkBsX6(MM7$I0Z>NSnEj?ERAu+bRc*R`` zJQE7Y*xE~txEwwSfLG~K&${|1)NdcNe@QRwx}jF;n1$G)ntNyBU_8Xj@y^^Q@u$_QutKWe$A?@WG?s{hu4URvYA$?P6pfHGJSB>;z zpEpbcSELn=whZJY4NEgkVP`-%*c33>uN_bZa%~<#6hAa+kOqD7BfwC!WN`FQold;? zz~m~PnjstT=rD~j^gRqfznUFCZ{Q~lfCqXZof-sM*WhlN+js~d46vCC0pJi24vpu+ zgl!;=foa(&Sy(NaQ%UagCHb{~Bxh%o0*(~N$$`AEG3e+#asz^IZy?O#d_(`E6hVjv zxpsQR=kpNzIw;pyom=1fT_jch;z6jT3mr!IqE;u|L9GB8`y__7%-xYOAM61HK99{0 z%eo2tJD#4;82y`h4!W`Nkh}L}0Z!mPbcEEW%uY@Py{{&&Wk?P5b%=&dZ1@zhVzIaB z5vxfBG8P0%1Kt#EcxytT8nBo2luM!{2)h4*BZpV&2>4EMOMr#-%V5<&{&lNT3_&fS zzng_H!SGSNT%28H0n4_W(qqWEat*+Jd53$XOpErx4I`0@a#%}0@CChSn2QP;Ww4xM zMc%!+QSPi<@kaxQC+(w+lO|tx4ZChpt09cnAk>jAjy$}vcdl&C>iP7tdlELNqL~DY zKSxmuvd&+ak}&B<1>?YIdos1q(*}Tw1Beppf=F6;=*r%K%}BrHU`4Tz zg%A%X?TMO2m)M~LUQN6|=Ntc&0P%EJ7K6otxm0ZJ`^BTJJVW|ub}mB|dwLBNpzGKK zPfN_Rud{wuiHQqbf^8AYZ0RcIfmQ&3wU~p(6lOL2JhaS&iUQX#kjCK$R`pe#)5#Q>Td z7Kc$4j`qfbx57u5V-|>Dg5!f{+}=_f8x2}T?!gp2yCe6pryU<6{ML{EpYyXBZra+~ zs9X;oM%RU}RH&J>pyANAwCi2GkhtHzCsYliwRg+4YgfXTd6E)SyBNouFGw<8DfhJZ zq=8_Hq1fDZQt04vh`XU3EoVjHI2_P%ANT*=P++339|0vaP2kvsqOc~E6`Y@s1kP4D z(+(6MphTf8niq)U2(Wmf#o3;d(Pz+7miQ&H2sK(lU+)(-17SNEKF#l9TFmHshF!`l zn&LAK@KCvu5r$Ak05=OX4L-tm97c6LVqNHKR3M4DM{r~U3UsO(~%>bABE{^IIAzq z`cVXB$p!qWjl&10?U@r$VfDXu6(C4ZP|(=;cp$iId;9w{UcV+rt98i8loS;M#;b$? z1^^^&aYx-5P^b*$@gE6&E?F45)NP4b=y@D9{)x(})%aq%B>RW;OGtp`mW!EnG3Q^k zah_gYY$7E^MS(RnGUVjs-|cT&S&@V)o0kp^4i0{5>+53zqbq7G7@wUb1@;S!!IRLR zK>9|T%YYG9q`YTe&kl-(ubP2(P`Z6~0xS(k&cfQvnJdBIg!;af>0MU@on> z0*DDLrdL6dD6RdU+YQpvJ_4B+mykddtxAx6-Bu5f<{1F}1TNDb0--AI?v&HoL3xVL z)cXV|=}H>+^YW5bc|in$(26xP>PtOc+%Z|1G?vZREaS<&oF|e$)1-*Xehm2l#@;^$ z*X-IUZ^&7IH&nlT8H}_?5njlqfSiV75^T#Z!w2>b4n-czCV=b^1PZMmJbzAAW!Z87nXcFp0-nOQgB9g3v907JRvbOG(;MqL<^u@$erQpvfUxDx)xTkxz+aRFaV4Z$BLF7b&=#k*tGE{2(BxhCM)K^>TIDM6bBF}3VsuM480tbEep z`{@$~rKEs+)h}Lbet6{4NwC&z@9yFw=otwA8?!4^oSsI^-ak5q)`WwZ{yyj?0H=w; z3X2*wUS5bt5Y)2r{5e;8db-o}2R5|U9Zc(|=g0XwnUfzz7v<0w9(ZyiW&s5pnfLfH z@sN(gKJAD!0!H~yYFtMGUD=iOyU3aUq2TR)+r*c@;j&wq`M3GL9C4Rv4d_&}<`6RQ znc$+LPRVk^i#(U{m;ehccy!Uc=GU)Zv^+dCJUl$eWu;%4sea0Eqy^YMioh#fyZPO& zf;0HR36;MHnDMrzora5Ds`Z;o-FQbZw%udgcU1rO7y2zf(0@o8)VO1}qzyCVhqH94 zfe4TfcZLJbeS|uXCwuyxX(9Cnr`bWO##{$eCRc;U36j{++n)&3Ey+;?y3=j~0TaNZ z);wQ&z8LUwpOfo7x9e1jcF&l)7fRuipkseV>xa~`dk0aA~mjS$Ur z_0yol29P};SZ%=<39p+KFOFy@0%s`yy}x&F>t(eN7;P)8m9V|5)zX9JhqEb2z;D4q zmIohb2zcF&g^NBOdz&CJFsNSS#EB05DmI!Xdr9YuDm@I~S5;Od`_hhMGeA** zaA~enrT4G`+|VB%*m_Qz`mTvFV5fuN>F)P-zB=oh%~&o$mh!xDL;*uT>cZgD0@_(x z)DiVM-^+l@DB)h6*7d*6Sm@Dg3w!p7(M1%1z(3wnO@PSX7UNB#BUYTJzcOQ6BsMAh zl74I_q7IGM9b=(;dy2w%^q-@g&gr7&4g_u#WqsJWJI+^+79UozW70=~vA#KNW{NZl zMSufw4aha-WnXimiEc)4x;pvups>T%X7%Ok<<~q{SjLpKLOc3sp}SaCE=-&UKGQAW zOwjdw6m;||KX}h{ZP!PI_Gx_y?k&(>%px@}*+&nak|25#*#+=Dx-y4%8zzlhLcg3| zl^UQ1FTfhDs1zRk-0Bl(c)*SZdv~A}Pg$2$_8D16`xyXUfd4sAg(V6#n%A5t)<%sO zu4n#n=0Lxd?)-8~s$WFMzDQCpNy9U-HmIv#`57Urenb(HpAX#QV(vDUx1f%VS;*Az zCgUM8Rzi@pf9WaiLrOz6f91F7|s@|fIz9cCDYn7c8r80U^tP-t9wyz z(asSnOt_1dcLTp~pr4HE@P!Uk5f(RTTNXo&EBc(q!Dj13!v7G$r4KmhlFt0?ExZ?g z;AN}e=dYDHd7jtu?qS^)=Q{pLi5DD?uaGH&kEB61qQe57H`c_Z_N?$ zX9oYmCi3un)z#~Pj8ov5xYl*n(5U4)p=3brk-UZ;2dE8p( z82zbyFh`i#=qaEOZ2sh>%=KUmVI^g6I&HB;KpOgF#EKF=LA+znd6K<W*0!ft7|8+!_iQ0e6s+w-S+i$_Gg1=9{6jPuxn> z;oA*AjWoJ2cq=Vg9tW!z_O*~lacjw>4<0fSt=TR%Q-CowLc-7P2)+N&X`8mXLl2ri>UU!#Cjtl<1`u@uOns8FDNml@)VRF}<>lp_ z+}+zgrKuDFc@w|LKb!=9bPZo2I;bc@kTWE(hZb|kP2!J(9CAfpWV&78=HeT|x34Xu zjG5iU|1W4_T^$_@C`eh^+G+ zfj(}|H?GmBZs@K5i!y<#{QBCewG+Sb1u;2VS~)6i_cH#W%PlL-Uj@&D-N*ljV<%3> zaLYx~@vEX0`2g9F0LV1Zr#kC@C<}+>9N2zlJ9WEpaF8PvP91XSGc_@xt}>#|)1f$X zrlMv?=tVMTwch=zR4$vM{had|Q2xiK3`6=nYBczD`(#DO-{|?jJ6WRCca1>tUmYz$ zfRK`tqvt$XSXiJ{O;4Y)gAs#QK!6Y!Y=^HIaC*wp(ulZ>nn&xmD_<@Caa{<=|5t?=Z=Cw~x>uQ*$rIUlh@_+>s-A;-e;4R_&&UdlSnoKyNl2M7M3*W6fzKPKK5~8MD)A zWj+OzK%y%f@Ox_(s-2;5i?g#{5Ef8E=?Is5?BjnUx)v`TZ17j$_S-hk1#VMQZsU55 zWPwo?edlA5Q$uhxj1}o3uXkO*8IKeqM3(tGL3~i%ot_|Emfc#MjK#=o%GtZGCEzERrn*Y9to(J{GxTEpV#?(FB49ldhGS( zWLiMWz!0G73hS9WoFKp4^fZD4bfB|nL~n9(fXk->%-!;{ZIv+F#kY1DKVPReEi;^f*(jAH9*6T=`%sCj zuRAGgo#&FjP|NEs-Ajb~y9RAb6o8O9&!QTzjHz;xoPd(h`u#Zl-8l`sl;J1j;IQ2k zO(s?Fpr&xC@Z6UQ2CnGNZ^teD-ZS&!^#wxB%AVsSbXtZGzS{=mxw9NUP38^T0P13n z!FzypnN=3|lf8yI*&DhS9bJ3i=rXo=NJkw0TX$z3P37A6@f|{KW1)~CN=T@PB2wZ| zk<%oEO@<~rge`O0C`C>wB@`7UQzUFtR8oYLdB{9vN`{Qj_jbGHXa+J;QtBTKJ#nBoUw!)rvANt6sQCoF zpI!RYY_yXIbRmKW*4bInzF1p|os#WZ{+TTM;QWiXpJnD2_luPr=KvcD!1(Xf_3KT) zJH6kHcrae=*g?NnKK46&Z0c`$A5^U;Jj`KTz1?R6A5#t#7v+xa*)efnMms;VOYWHV zAcgeEfFD`Z(*+IR4>>v{UrybO4H7xlV>*+U?sN_fQG$Bl?Pqq=t+t#{C{$X=`oHeY0x zfRApmO!)RrYdoP7d0K5Re>>w(+j!~g{rFAV58}5vX!z?ejm+F~07J_KzO9XE6P8k> z2ku##|0)Yw*=TFuE7@$l6@{J>Z2{Yy=kxn1j{L{2X5B zNhIm2YiE#Z9*Eg1sSUrH) zs=Rc~>U#}Xe5o+!B%wrCl>?I#KQhFpJkHxQNk|+T=xe)Hu%DBOq&;k~$3c)S;7JWO zP2ftFC6Ea`e3s8E$&}M3Mi;aqJ>qitp)c4nHN-*sR>92* zX4vu0tAM>2Dm*6NKG4TTqB++Jdirb$9b>skWKfEPLY)zM^I%*xGs$d?>L*p+Fk^v` zVY+DZ8k{mgZzDth42MRFcuphU6+xIrsB5p?o;x>lYW}UB_PRe}ubOl=pXUuCbVD7J zTa>ZOxzZH!*i~X|o21#{QWbbCDj{*^4&mPmWa`9ESd$)kxF$UlUC_&_w8~!W?*#g+ zl~}O{+tGxnJe3lOSwc&q4LE#YrYEv%NFc(Fm~o;` z9K22R?Y#4`(c!!v&ST4Nx?lwvnJP(@Bh^%7Q80WiA^D5(h6#>_sR`xZX`P-EmusEa zIG;YL9`N&qwNucm_g}+|!Odt6$L7xE^I#$otPXD<3vQ@IBw|gkvJd*949Y~h{Vuw& z>knB7o+07?ERy(?`ew}HNi6K7#Zbn(==cykjo8MYlXvGFV+xp{YoX15v`%`KmGEtZ!Oz-n)>8GMw z!?~AQUr^uVIGI$W^&w)T1Nfp`Fu9Y3UUH z6ln>GejaE+_fdtoZ#c7>va#+uU3)k8hk&zQ|I+jFdJ?hnu?ewVVn6S1%^JWS+IFNy z`CF^Urr`3l@Y^c6uK;)Yu@7zb8#~q%pgeTg)R{Hqn%Wfzxx$`F9JrkB0YHVVg+E3 zZ-%`m)d>6;L7-`8v*N zoGQN*jaDRMue<4SZGERp7FQ@0`31VHMsa5fgG~({VE4r_HK)ltvE2WjQ==SrVtmSL z17sIaJM^{|XcFBB&#CjIT?h&YKhgdBM5Qs)9x?gD!?I~`%9q_!CX+z04S|}}UP*}1 zqYxgjAYMtLI;yFuNh(A}rbrvHyszma3J~Y!v%_6?EOa0T$7)}@V}7x;pV<=p!|3BN zOb{1kH8fgw6d9o^$^!?dM8Ecia21%H)3vp%#xnE>iGkDVK&R0rK`{<0pmFbnK z|GxR}e=XBLJOOxF=DEgzP%uL>`;!tA_qe-fS0!zH`c>JgDj=@XOx7kPRQj~b{VmBGd6h{5e0+t30rT#Emt8XidIsD*JcJvN4oYm^ z%xmHOYyBD#k+$~sOHh%_judmJrk(14CodxdeE>Ptdv3bXVBs%GiK+8&?BkEMHR7A& z*vT4(?_QJNqLxST>SX5_8UoO-B3}7 z$RMB*9Etre$x%|<)(fL&Y56*bZ_Ot;deZ(&E@T}`nlA%oYTmH-IxJb%$XB4Lqby`6 zMoc+K0KLWXSCjMmk6zJ;&gY-Nc(uk=4klM89#eo;2OwL)O;mgPU8hdT9F{$jI0$8L z|G)sU0^2WSV@SX33BU%FHehb-k_{KS4MjSqn(OyMRIj0A_8rcj;UWkhBxtpJ$1f6Q zTNvMH-FNLw)L-T0?fgE-_+&Gh#QCEsVb!iH?S25gGhkmsB$1bp4)lm^+cxj2Q5+(o zHZsvj?LaZ2wFF82IJpdHIBO8#qt+l$4Oz=ZfCh?Xwyo8M99u9b+OZ$zwV*l07DzsN~+VfoBF;z!^7nS30Z`I(x`|D2KBabE14v(@Wr5>L^ni)=lI=8utqfZ za&`8(3d1tusd#_*Uf|Vi_rnp^$3Bd_u?yUs$`D42n=a_@1%XE(l;7zZhN+1i+kGbs zZ}spf1u`-<)+9x>vXKZ$tw4Vu%7(v{%D&y`7+)pL9Q%Cqi_ddW8>oO$Xez71IgPWpT( zyj~U6Gl}yrd!T2qpdDQbMXFSKs)QQG@qq3Q*D=Q?Ap+xk_3CnxFI154!{I!uih7oi0Q}A>28Ypa zEN4OrCV%@o3|b^-5aP_KU%%Qw%Z8F|pz>TK^H5CJK>y79xpdDHUSiokZiatHMMra% zh!%Q?hwI~j_TznETAR)HOB+pSY@benHOSSgY)>7BIAE>;6w!-cL1n?64M)`* zzv*xD4N`9s0#kN&c3*#gW>8P*unD==wt@{tfOBD^((-4e0C1`kO1|VVU7(6GF=LMFp(VJKi*9iIJ{ygEqoSR-| z1>Ro4(;_JK)TDLRj)1>SP4hvpyw$c$@wv4@N*j*v5d;y%;A!K(fms;b)VNMc%Bn7v zlKSeEZj6Q(5mN{ZH1+i~+#0$d*3sUsX!C(jI$WP9v?m?yC=v5fGU>Ul(qca<#LC7N zd0wHxHGbw!L_}S`N5+E(h4wq0-8qLZESSuA>OsG>TCp%SWP6@BLX_2;tys)XoB;KV zt}Eun2bAfnvv#LScJC`9)h)BKwoW@+za@8mdO+3n`|1D=!Apn0IG$QilxCKec3vf}AiQQnb4cw}> z?4IbAELuE2FB0O4%Wkb+wIej38!V|*1y4)|kgD^aXzrVDe*W;G)_E5f{`{$Sa`6lF z)lG<@6gzUhFmUg)9GZ7%$=`XEotHczlx&`7Li{AKC!8P429?jU?yAf8@?Z?iLUPPc zTx`bKiSPAA07G+54v)8%-Ey!KGpz<0yRBTMRjJ7bxeR^no`4Hx6p>EX32JM;7Za%s zHPUm**2yVt&AHx}FMn++zU3K^!kpslle@0pnY@ZiU;?G9$!Hz9hjYuA%A4)24DqF-L)JDn>j_6Jk{9(K?m4pVV z^}<~1LR%y?|Fuf1Co>kY$!B4XcriFl6WA5jn*|C@fzO)x1tH4wzj$n)z(fPP(PB}x0Fj3&MT@g`pT^%Jk zi8Nhs-Up2@P(=6$HCvMDD+yUy{cFEdO%ERAfvSsRvgDeeKb+K7e(~aNshbZ|Q!DM=f#D_+nsfMvF`%+A zX^?qJ)-Bzq9E-a1xyj4rnsYx_vT}y*qaIIAitfiLzv7P`75&Ub+PRU>-|8obT$+8wH z5Be_3a=1-(47U^`?yBowPQR?0CJ|W7;%PrWKjcX|13<5Lz`1+Qq5ANU?Sk)XIe?b3WOk0)usOyCF)4x(PbzH;aLhK z8n|s?x}v;H%*=FT)ak7S#CKg?U5n6z-b4_Y9QgV93F+eAz4VbD&4uGoGD5iYWCTWg zGDt)#rzZwC67RF=`Q~zsd;lbi@R&)P<bXk^e zSK*M>)MQ?aFsp59w4(^o!%|WaUVY+gG`-C*N#}()`=wN?x63Jnoqb&4M9w;R%9r{4 z`e<<(imx^dh7mZqve)zyKAZL+phH(mf zc*uznCu`;8_K*Dh$p`aQnfJoc6!T4<7*SJ%Q;&}R@iIIWu?WiO*ktuYU=ZMz<5Lc( zq8^3tjW|@_H!oOaPSb*AVYii4Y}HIsPEJE=$YU|`CQfF{4Zntc7B&LXL?%`b+lh|N zBP<%J>KHqM**U^V1z`9*v?%bt6crVvXjw>V0{~s+qFG8U8OKoTQYw}Dx~0V&N`M!D z|5OPwK0i~dFM#@PImx7nxPj1#KFrB6p}||e^zrX9K+`Wrt&LA$lgYXUznAq<-#hjx zPi*t|HYpHXn;bzCb4OExRoR` zs^>hgigX;_WG*hQOAH3%JkXP5?B7k9XYihB&eg@uVA0jpC0IUxu4Byvr*BoyvyE%W zo33c{u)IFm98ARbCtZ<0{m-H5|Hn~W^O}4tM$+4Bpq;$i^3p#5L6(zb diff --git a/TB2J/abacus/abacus_wrapper.py b/TB2J/abacus/abacus_wrapper.py index 033401c..10dde17 100644 --- a/TB2J/abacus/abacus_wrapper.py +++ b/TB2J/abacus/abacus_wrapper.py @@ -16,7 +16,7 @@ class AbacusWrapper(AbstractTB): def __init__(self, HR, SR, Rlist, nbasis, nspin=1): - self.R2kfactor = -2j * np.pi + self.R2kfactor = 2j * np.pi self.is_orthogonal = False self._name = "ABACUS" self._HR = HR @@ -70,6 +70,8 @@ def gen_ham(self, k, convention=2): S = self.SR[iR] * phase # Sk += S + S.conjugate().T Sk += S + # Hk = (Hk + Hk.conj().T)/2 + # Sk = (Sk + Sk.conj().T)/2 elif convention == 1: # TODO: implement the first convention (the r convention) raise NotImplementedError("convention 1 is not implemented yet.") diff --git a/TB2J/abacus/splitSOC.py b/TB2J/abacus/splitSOC.py index cf4e6a0..933dcc3 100644 --- a/TB2J/abacus/splitSOC.py +++ b/TB2J/abacus/splitSOC.py @@ -4,6 +4,7 @@ from TB2J.kpoints import monkhorst_pack from TB2J.mathutils.fermi import fermi from TB2J.mathutils.kR_convert import k_to_R, R_to_k +from scipy.linalg import eigh from copy import deepcopy from scipy.spatial.transform import Rotation import matplotlib.pyplot as plt @@ -25,10 +26,11 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) self._HR_copy = deepcopy(self._HR) self.HR_soc = HR_soc + self.soc_lambda = 0.1 @property def HR(self): - return self._HR + self.HR_soc + return self._HR + self.HR_soc * self.soc_lambda def rotate_HR_xc(self, axis): """ @@ -44,25 +46,28 @@ def rotate_Hk_xc(self, axis): for ik in range(len(self._Hk)): self._Hk[ik] = rotate_Matrix_from_z_to_axis(self._Hk_copy[ik], axis) - def get_density_matrix(self, kpts, kweights): - evals, evecs = self.solve_all(kpts) - nkpt = len(kpts) - rho = np.einsum( - "kb,kib,kjb->kij", - fermi(evals, self.efermi, width=0.05), - evecs, - evecs.conj(), - ) - rho = np.zeros((nkpt, self.nbasis, self.nbasis), dtype=complex) - for ik, k in enumerate(kpts): - rho[ik] = ( - evecs[ik] - * fermi(evals[ik], self.efermi, width=0.05) - @ evecs[ik].T.conj() - * kweights[ik] + def get_density_matrix(self, kpts): + rho = np.zeros((len(kpts), self.nbasis, self.nbasis), dtype=complex) + for ik, kpt in enumerate(kpts): + Hk, Sk = self.gen_ham(kpt) + evals, evecs = eigh(Hk, Sk) + rho[ik] = np.einsum( + "ib, b, jb -> ij", + evecs, + fermi(evals, self.efermi, width=0.05), + evecs.conj(), ) - print(np.trace(np.sum(rho, axis=0))) return rho + # rho = np.zeros((nkpt, self.nbasis, self.nbasis), dtype=complex) + # for ik, k in enumerate(kpts): + # rho[ik] = ( + # evecs[ik] + # * fermi(evals[ik], self.efermi, width=0.05) + # @ evecs[ik].T.conj() + # * kweights[ik] + # ) + # print(np.trace(np.sum(rho, axis=0))) + # return rho def rotate_DM(self, rho, axis): """ @@ -79,6 +84,20 @@ def __init__(self, model, kmesh, gamma=True): self.kpts = monkhorst_pack(kmesh, gamma_center=gamma) self.kweights = np.ones(len(self.kpts), dtype=float) / len(self.kpts) + def get_band_energy2(self): + for ik, kpt in enumerate(self.kpts): + Hk, Sk = self.model.gen_ham(kpt) + evals, evecs = eigh(Hk, Sk) + rho = np.einsum( + "ib, b, jb -> ij", + evecs, + fermi(evals, self.model.efermi, width=0.05), + evecs.conj(), + ) + eband1 = np.sum(evals * fermi(evals, self.model.efermi, width=0.05)) + eband2 = np.trace(Hk @ rho) + print(eband1, eband2) + def get_band_energy(self, dm=False): evals, evecs = self.model.solve_all(self.kpts) eband = np.sum( @@ -94,9 +113,24 @@ def get_band_energy(self, dm=False): def calc_ref(self): # calculate the Hk_ref, Sk_ref, Hk_soc_ref, and rho_ref - self.rho_ref = self.model.get_density_matrix(self.kpts, self.kweights) - self.Hk_xc_ref = R_to_k(self.kpts, self.model.Rlist, self.model.HR) + self.Sk_ref = R_to_k(self.kpts, self.model.Rlist, self.model.SR) + self.Hk_xc_ref = R_to_k(self.kpts, self.model.Rlist, self.model._HR_copy) self.Hk_soc_ref = R_to_k(self.kpts, self.model.Rlist, self.model.HR_soc) + self.rho_ref = np.zeros( + (len(self.kpts), self.model.nbasis, self.model.nbasis), dtype=complex + ) + print(f"{self.Hk_xc_ref[0][:4,0:4].real=}") + print(f"{self.Sk_ref[0][:4,0:4].real=}") + for ik, kpt in enumerate(self.kpts): + # evals, evecs = eigh(self.Hk_xc_ref[ik]+self.Hk_soc_ref[ik], self.Sk_ref[ik]) + evals, evecs = eigh(self.Hk_xc_ref[ik], self.Sk_ref[ik]) + self.rho_ref[ik] = np.einsum( + "ib, b, jb -> ij", + evecs, + fermi(evals, self.model.efermi, width=0.05), + evecs.conj(), + ) # @self.Sk_ref[ik] + print(f"{self.rho_ref[0][:4,0:4].real=}") def get_band_energy_from_rho(self, axis): eband = 0.0 @@ -104,7 +138,23 @@ def get_band_energy_from_rho(self, axis): rho = rotate_Matrix_from_z_to_axis(self.rho_ref[ik], axis) Hk_xc = rotate_Matrix_from_z_to_axis(self.Hk_xc_ref[ik], axis) Hk_soc = self.Hk_soc_ref[ik] - eband += np.trace(rho @ (Hk_xc + Hk_soc)) + Htot = Hk_xc + Hk_soc * self.model.soc_lambda + Sk = self.Sk_ref[ik] + # evals, evecs = eigh(Htot, Sk) + # rho2= np.einsum("ib, b, jb -> ij", evecs, fermi(evals, self.model.efermi, width=0.05), evecs.conj()) + if ik == 0 and False: + print(f"{evecs[:4,0:4].real=}") + print(f"{evals[:4]=}") + print(f"{Hk_xc[:4,0:4].real=}") + print(f"{Htot[:4,0:4].real=}") + print(f"{Sk[:4,0:4].real=}") + print(f"{rho[:4,0:4].real=}") + print(f"{rho2[:4,0:4].real=}") + # eband1 = np.sum(evals * fermi(evals, self.model.efermi, width=0.05)) + # eband2 = np.trace(Htot @ rho2).real + # eband3 = np.trace(Htot @ rho).real + # print(eband1, eband2, eband3) + eband += np.trace(Hk_soc @ rho) * self.kweights[ik] * self.model.soc_lambda return eband def get_band_energy_vs_theta( @@ -118,12 +168,15 @@ def get_band_energy_vs_theta( es2 = [] # e,rho = self.model.get_band_energy(dm=True) self.calc_ref() - thetas = np.linspace(*angle_range) + thetas = np.linspace(*angle_range, npoints) for theta in thetas: axis = Rotation.from_euler(rotation_axis, theta).apply(initial_direction) self.model.rotate_HR_xc(axis) + # self.get_band_energy2() e = self.get_band_energy() + # e=0 e2 = self.get_band_energy_from_rho(axis) + # e2=0 es.append(e) es2.append(e2) print(f"{e=}, {e2=}") @@ -164,8 +217,8 @@ def parse(self): def test_AbacusSplitSOCWrapper(): - # path = Path("~/projects/2D_Fe").expanduser() - path = Path("/home/hexu/projects/TB2Jflows/examples/2D_Fe") + path = Path("~/projects/2D_Fe").expanduser() + # path = Path("~/projects/TB2Jflows/examples/2D_Fe") outpath_nosoc = f"{path}/Fe_soc0/OUT.ABACUS" outpath_soc = f"{path}/Fe_soc1_nscf/OUT.ABACUS" parser = AbacusSplitSOCParser( @@ -173,16 +226,16 @@ def test_AbacusSplitSOCWrapper(): ) model = parser.parse() kmesh = [4, 4, 1] - e_z = get_model_energy(model, kmesh=kmesh, gamma=True) - print(e_z) + # e_z = get_model_energy(model, kmesh=kmesh, gamma=True) + # print(e_z) - model.rotate_HR_xc([0, 0, 1]) - e_z = get_model_energy(model, kmesh=kmesh, gamma=True) - print(e_z) + # model.rotate_HR_xc([0, 0, 1]) + # e_z = get_model_energy(model, kmesh=kmesh, gamma=True) + # print(e_z) - model.rotate_HR_xc([1, 0, 0]) - e_x = get_model_energy(model, kmesh=kmesh, gamma=True) - print(e_x) + # model.rotate_HR_xc([1, 0, 0]) + # e_x = get_model_energy(model, kmesh=kmesh, gamma=True) + # print(e_x) r = RotateHam(model, kmesh) # thetas, es = r.get_band_energy_vs_theta(angle_range=(0, np.pi*2), rotation_axis="z", initial_direction=(1,0,0), npoints=21) @@ -190,11 +243,15 @@ def test_AbacusSplitSOCWrapper(): angle_range=(0, np.pi * 2), rotation_axis="y", initial_direction=(0, 0, 1), - npoints=21, + npoints=11, ) + # print the table of thetas and es, es2 + print("theta, e, e2") + for theta, e, e2 in zip(thetas, es, es2): + print(f"{theta=}, {e=}, {e2=}") - # plt.plot(thetas / np.pi, es, marker="o") - plt.plot(thetas / np.pi, es2, marker=".") + plt.plot(thetas / np.pi, es - es[0], marker="o") + plt.plot(thetas / np.pi, es2 - es2[0], marker=".") plt.savefig("E_along_z_x_z.png") plt.show() diff --git a/TB2J/abacus/test_density_matrix.py b/TB2J/abacus/test_density_matrix.py new file mode 100644 index 0000000..589d901 --- /dev/null +++ b/TB2J/abacus/test_density_matrix.py @@ -0,0 +1,38 @@ +from scipy.linalg import eigh +import numpy as np + + +def gen_random_hermitean_matrix(n): + A = np.random.rand(n, n) + 1j * np.random.rand(n, n) + return A + A.conj().T + + +def gen_overlap_matrix(n): + A = np.random.rand(n, n) + 1j * np.random.rand(n, n) + return np.dot(A, A.conj().T) + + +def fermi_function(x, ef, beta): + return 1.0 / (np.exp(beta * (x - ef)) + 1) + + +def test(): + n = 10 + A = gen_random_hermitean_matrix(n) + S = gen_overlap_matrix(n) + beta = 0.1 + ef = 0 + + evals, evecs = eigh(A, S) + + etot = np.sum(evals * fermi_function(evals, ef, beta)) + + rho = np.einsum("ib,b,jb->ij", evecs, fermi_function(evals, ef, beta), evecs.conj()) + + etot2 = np.trace(np.dot(A, rho)) + + print(etot, etot2) + + +if __name__ == "__main__": + test() diff --git a/TB2J/mathutils/kR_convert.py b/TB2J/mathutils/kR_convert.py index 06aa8d0..eaa2501 100644 --- a/TB2J/mathutils/kR_convert.py +++ b/TB2J/mathutils/kR_convert.py @@ -59,12 +59,12 @@ def R_to_k(kpts, Rlist, MR): phase = np.exp(2.0 * np.pi * 1j * np.tensordot(kpts, Rlist, axes=([1], [1]))) Mk = np.einsum("rlm, kr -> klm", MR, phase) - nkpt, n1, n2 = Mk.shape - nR = Rlist.shape[0] - Mk = np.zeros((nkpt, n1, n2), dtype=complex) - for iR, R in enumerate(Rlist): - for ik in range(nkpt): - Mk[ik] += MR[iR] * np.exp(2.0 * np.pi * 1j * np.dot(kpts[ik], R)) + # nkpt, n1, n2 = Mk.shape + # nR = Rlist.shape[0] + # Mk = np.zeros((nkpt, n1, n2), dtype=complex) + # for iR, R in enumerate(Rlist): + # for ik in range(nkpt): + # Mk[ik] += MR[iR] * np.exp(2.0 * np.pi * 1j * np.dot(kpts[ik], R)) return Mk From addd42988e2317bf7b57fe840a769138e6ab7ab5 Mon Sep 17 00:00:00 2001 From: mailhexu Date: Sun, 26 May 2024 00:36:07 +0200 Subject: [PATCH 6/7] fix splitSOC occupation --- TB2J/abacus/occupations.py | 278 +++++++++++++++++++++++++++++++++++++ TB2J/abacus/splitSOC.py | 88 ++++++++---- 2 files changed, 337 insertions(+), 29 deletions(-) create mode 100644 TB2J/abacus/occupations.py diff --git a/TB2J/abacus/occupations.py b/TB2J/abacus/occupations.py new file mode 100644 index 0000000..d279198 --- /dev/null +++ b/TB2J/abacus/occupations.py @@ -0,0 +1,278 @@ +""" +This file is stolen from the hotbit programm, with some modification. +""" + +import numpy as np +from scipy.optimize import brentq, brent +import sys + +MAX_EXP_ARGUMENT = np.log(sys.float_info.max) +from ase.dft.dos import DOS +from scipy import interpolate, optimize, integrate + +# import numba +import math + +# from numba import float64, int32 + + +# @numba.vectorize(nopython=True) +# def myfermi(e, mu, width, nspin): +# x = (e - mu) / width +# if x < -10: +# ret = 2.0 / nspin +# elif x > 10: +# ret = 0.0 +# else: +# ret = 2.0 / nspin / (math.exp(x) + 1) +# return ret + + +def myfermi(e, mu, width, nspin): + x = (e - mu) / width + return np.where(x < 10, 2.0 / (nspin * (np.exp(x) + 1.0)), 0.0) + + +class Occupations(object): + def __init__(self, nel, width, wk, nspin=1): + """ + Initialize parameters for occupations. + :param nel: Number of electrons + :param width: Fermi-broadening + :param wk: k-point weights. eg. If only gamma, [1.0] + :param nspin(optional): number of spin, if spin=1 multiplicity=2 else, multiplicity=1. + """ + self.nel = nel + self.width = width + self.wk = wk + self.nk = len(wk) + self.nspin = nspin + + def get_mu(self): + """Return the Fermi-level (or chemical potential).""" + return self.mu + + def fermi(self, mu): + """ + Occupy states with given chemical potential. + Occupations are 0...2; without k-point weights + """ + return myfermi(self.e, mu, self.width, self.nspin) + + def root_function(self, mu): + """This function is exactly zero when mu is right.""" + f = self.fermi(mu) + return np.einsum("i, ij->", self.wk, f) - self.nel + + def occupy(self, e, xtol=1e-11): + """ + Calculate occupation numbers with given Fermi-broadening. + + @param e: e[ind_k,ind_orb] energy of k-point, state a + Note added by hexu: With spin=2,e[k,a,sigma], it also work. only the *2 should be removed. + @param wk: wk[:] weights for k-points + @param width: The Fermi-broadening + + Returns: fermi[ind_k, ind_orb] + """ + self.e = e + eflat = e.flatten() + ind = np.argsort(eflat) + e_sorted = eflat[ind] + if self.nspin == 1: + m = 2 + elif self.nspin == 2: + m = 1 + n_sorted = (self.wk[:, None, None] * np.ones_like(e) * m).flatten()[ind] + + sum = n_sorted.cumsum() + if self.nel < sum[0]: + ifermi = 0 + elif self.nel > sum[-1]: + raise ("number of electrons larger than number of orbital*spin") + else: + ifermi = np.searchsorted(sum, self.nel) + try: + if ifermi == 0: + elo = e_sorted[0] + else: + elo = e_sorted[ifermi - 1] + if ifermi == len(e_sorted) - 1: + ehi = e_sorted[-1] + else: + ehi = e_sorted[ifermi + 1] + guess = e_sorted[ifermi] + dmu = np.max((self.width, guess - elo, ehi - guess)) + mu = brentq(self.root_function, guess - dmu, guess + dmu, xtol=xtol) + # mu = brent( + # self.root_function, + # brack=(guess - elo, guess, guess + dmu), + # tol=xtol) + except: + # probably a bad guess + dmu = self.width + if self.nel < 1e-3: + mu = min(e_sorted) - dmu * 20 + elif self.nel - sum[-1] > -1e-3: + mu = max(e_sorted) + dmu * 20 + else: + # mu = brent( + # self.root_function, + # brack=(e_sorted[0] - dmu * 10, + # guess, + # e_sorted[-1] + dmu * 10), + # tol=xtol) + mu = brentq( + self.root_function, + e_sorted[0] - dmu * 20, + e_sorted[-1] + dmu * 20, + xtol=xtol, + ) + + if np.abs(self.root_function(mu)) > xtol * 1e4: + # raise RuntimeError( + # 'Fermi level could not be assigned reliably. Has the system fragmented?' + # ) + print( + "Fermi level could not be assigned reliably. Has the system fragmented?" + ) + + f = self.fermi(mu) + # rho=(self.eigenvecs*f).dot(self.eigenvecs.transpose()) + + self.mu, self.f = mu, f + return f + + def plot(self): + import pylab as pl + + for ik in range(self.nk): + pl.plot(self.e[ik, :], self.f[ik, :]) + pl.scatter(self.e[ik, :], self.f[ik, :]) + pl.title("occupations") + pl.xlabel("energy (Ha)") + pl.ylabel("occupation") + pl.show() + + +class GaussOccupations(Occupations): + def get_mu(self): + return self.mu + + def delta(self, energy): + """Return a delta-function centered at 'energy'.""" + x = -(((self.e - energy) / self.width) ** 2) + return np.exp(x) / (np.sqrt(np.pi) * self.width) + + def get_dos(self, npts=500): + eflat = self.e.flatten() + ind = np.argsort(eflat) + e_sorted = eflat[ind] + if self.nspin == 1: + m = 2 + elif self.nspin == 2: + m = 1 + n_sorted = (self.wk * np.ones_like(e) * m).flatten()[ind] + dos = np.zeros(npts) + for w, e_n in zip(self.w_k, self.e_skn[0]): + for e in e_n: + dos += w * self.delta(e) + + def root_function(self, mu): + pass + + # @profile + def occupy(self, e, xtol=1e-8, guess=0.0): + self.e = e + dos = myDOS(kweights=self.wk, eigenvalues=e, width=self.width, npts=501) + edos = dos.get_energies() + d = dos.get_dos() + idos = integrate.cumtrapz(d, edos, initial=0) - self.nel + # f_idos = interpolate.interp1d(edos, idos) + # ret = optimize.fmin(f_idos, x0=edos[400], xtol=xtol, disp=True) + ifermi = np.searchsorted(idos, 0.0) + # self.mu = ret[0] + self.mu = edos[ifermi] + self.f = self.fermi(self.mu) + return self.f + + +class myDOS(DOS): + def __init__( + self, kweights, eigenvalues, nspin=1, width=0.1, window=None, npts=1001 + ): + """Electronic Density Of States object. + + calc: calculator object + Any ASE compliant calculator object. + width: float + Width of guassian smearing. Use width=0.0 for linear tetrahedron + interpolation. + window: tuple of two float + Use ``window=(emin, emax)``. If not specified, a window + big enough to hold all the eigenvalues will be used. + npts: int + Number of points. + + """ + self.npts = npts + self.width = width + # self.w_k = calc.get_k_point_weights() + self.w_k = kweights + self.nspins = nspin + # self.e_skn = np.array([[calc.get_eigenvalues(kpt=k, spin=s) + # for k in range(len(self.w_k))] + # for s in range(self.nspins)]) + # self.e_skn -= calc.get_fermi_level() + self.e_skn = np.array([eigenvalues.T]) # eigenvalues: iband, ikpt + + if window is None: + emin = None + emax = None + else: + emin, emax = window + + if emin is None: + emin = self.e_skn.min() - 10 * self.width + if emax is None: + emax = self.e_skn.max() + 10 * self.width + + self.energies = np.linspace(emin, emax, npts) + + if width == 0.0: + bzkpts = calc.get_bz_k_points() + size, offset = get_monkhorst_pack_size_and_offset(bzkpts) + bz2ibz = calc.get_bz_to_ibz_map() + shape = (self.nspins,) + tuple(size) + (-1,) + self.e_skn = self.e_skn[:, bz2ibz].reshape(shape) + self.cell = calc.atoms.cell + + def get_idos(self): + e, d = self.get_dos() + return np.trapz(d, e) + + def delta(self, energy): + """Return a delta-function centered at 'energy'.""" + x = -(((self.energies - energy) / self.width) ** 2) + return np.exp(x) / (np.sqrt(np.pi) * self.width) + + def get_dos(self, spin=None): + """Get array of DOS values. + + The *spin* argument can be 0 or 1 (spin up or down) - if not + specified, the total DOS is returned. + """ + + if spin is None: + if self.nspins == 2: + # Spin-polarized calculation, but no spin specified - + # return the total DOS: + return self.get_dos(spin=0) + self.get_dos(spin=1) + else: + spin = 0 + + dos = np.zeros(self.npts) + for w, e_n in zip(self.w_k, self.e_skn[spin]): + for e in e_n: + dos += w * self.delta(e) + return dos diff --git a/TB2J/abacus/splitSOC.py b/TB2J/abacus/splitSOC.py index 933dcc3..2d57f6c 100644 --- a/TB2J/abacus/splitSOC.py +++ b/TB2J/abacus/splitSOC.py @@ -9,6 +9,7 @@ from scipy.spatial.transform import Rotation import matplotlib.pyplot as plt from pathlib import Path +from TB2J.abacus.occupations import Occupations # TODO List: # - [x] Add the class AbacusSplitSOCWrapper @@ -16,6 +17,17 @@ # - [x] Compute the band energy at arbitrary +def get_occupation(evals, kweights, nel, width=0.1): + occ = Occupations(nel=nel, width=width, wk=kweights, nspin=2) + return occ.occupy(evals) + + +def get_density_matrix(evals=None, evecs=None, kweights=None, nel=None, width=0.1): + occ = get_occupation(evals, kweights, nel, width=width) + rho = np.einsum("kib, kb, kjb -> kij", evecs, occ, evecs.conj()) + return rho + + class AbacusSplitSOCWrapper(AbacusWrapper): """ Abacus wrapper with Hamiltonian split to SOC and non-SOC parts @@ -23,10 +35,14 @@ class AbacusSplitSOCWrapper(AbacusWrapper): def __init__(self, *args, **kwargs): HR_soc = kwargs.pop("HR_soc", None) + # nbasis = HR_soc.shape[1] + # kwargs["nbasis"] = nbasis super().__init__(*args, **kwargs) self._HR_copy = deepcopy(self._HR) self.HR_soc = HR_soc - self.soc_lambda = 0.1 + self.soc_lambda = 1.0 + self.nel = 16 + self.width = 0.1 @property def HR(self): @@ -46,17 +62,22 @@ def rotate_Hk_xc(self, axis): for ik in range(len(self._Hk)): self._Hk[ik] = rotate_Matrix_from_z_to_axis(self._Hk_copy[ik], axis) - def get_density_matrix(self, kpts): + def get_density_matrix(self, kpts, kweights=None): rho = np.zeros((len(kpts), self.nbasis, self.nbasis), dtype=complex) - for ik, kpt in enumerate(kpts): - Hk, Sk = self.gen_ham(kpt) - evals, evecs = eigh(Hk, Sk) - rho[ik] = np.einsum( - "ib, b, jb -> ij", - evecs, - fermi(evals, self.efermi, width=0.05), - evecs.conj(), - ) + evals, evecs = self.solve_all(kpts) + # occ = Occupations(self.efermi, width=self.width, wk=self.nel, nspin=1) + occ = get_occupation(evals, kweights, self.nel, width=self.width) + rho = np.einsum("kib, kb, kjb -> kij", evecs, occ, evecs.conj()) + + # for ik, kpt in enumerate(kpts): + # Hk, Sk = self.gen_ham(kpt) + # evals, evecs = eigh(Hk, Sk) + # rho[ik] = np.einsum( + # "ib, b, jb -> ij", + # evecs, + # fermi(evals, self.efermi, width=0.05), + # evecs.conj(), + # ) return rho # rho = np.zeros((nkpt, self.nbasis, self.nbasis), dtype=complex) # for ik, k in enumerate(kpts): @@ -100,11 +121,11 @@ def get_band_energy2(self): def get_band_energy(self, dm=False): evals, evecs = self.model.solve_all(self.kpts) - eband = np.sum( - evals - * fermi(evals, self.model.efermi, width=0.05) - * self.kweights[:, np.newaxis] + occ = get_occupation( + evals, self.kweights, self.model.nel, width=self.model.width ) + eband = np.sum(evals * occ * self.kweights[:, np.newaxis]) + # * fermi(evals, self.model.efermi, width=0.05) if dm: density_matrix = self.model.get_density_matrix(evecs) return eband, density_matrix @@ -119,18 +140,23 @@ def calc_ref(self): self.rho_ref = np.zeros( (len(self.kpts), self.model.nbasis, self.model.nbasis), dtype=complex ) - print(f"{self.Hk_xc_ref[0][:4,0:4].real=}") - print(f"{self.Sk_ref[0][:4,0:4].real=}") + + evals = np.zeros((len(self.kpts), self.model.nbasis), dtype=float) + evecs = np.zeros( + (len(self.kpts), self.model.nbasis, self.model.nbasis), dtype=complex + ) + for ik, kpt in enumerate(self.kpts): # evals, evecs = eigh(self.Hk_xc_ref[ik]+self.Hk_soc_ref[ik], self.Sk_ref[ik]) - evals, evecs = eigh(self.Hk_xc_ref[ik], self.Sk_ref[ik]) - self.rho_ref[ik] = np.einsum( - "ib, b, jb -> ij", - evecs, - fermi(evals, self.model.efermi, width=0.05), - evecs.conj(), - ) # @self.Sk_ref[ik] - print(f"{self.rho_ref[0][:4,0:4].real=}") + evals[ik], evecs[ik] = eigh(self.Hk_xc_ref[ik], self.Sk_ref[ik]) + print(f"{evals.shape=}, {evecs.shape=}") + print(f" {self.kweights=}, {self.model.nel=}, {self.model.width=} ") + occ = get_occupation( + evals, self.kweights, self.model.nel, width=self.model.width + ) + # occ = fermi(evals, self.model.efermi, width=self.model.width) + self.rho_ref = np.einsum("kib, kb, kjb -> kij", evecs, occ, evecs.conj()) + print(f"{self.rho_ref[0][:4, :4].real}") def get_band_energy_from_rho(self, axis): eband = 0.0 @@ -154,7 +180,11 @@ def get_band_energy_from_rho(self, axis): # eband2 = np.trace(Htot @ rho2).real # eband3 = np.trace(Htot @ rho).real # print(eband1, eband2, eband3) - eband += np.trace(Hk_soc @ rho) * self.kweights[ik] * self.model.soc_lambda + print(rho[:4, :4].real) + e_soc = np.trace(Hk_soc @ rho) * self.kweights[ik] * self.model.soc_lambda + + eband += e_soc + print(eband) return eband def get_band_energy_vs_theta( @@ -217,15 +247,15 @@ def parse(self): def test_AbacusSplitSOCWrapper(): - path = Path("~/projects/2D_Fe").expanduser() - # path = Path("~/projects/TB2Jflows/examples/2D_Fe") + # path = Path("~/projects/2D_Fe").expanduser() + path = Path("~/projects/TB2Jflows/examples/2D_Fe").expanduser() outpath_nosoc = f"{path}/Fe_soc0/OUT.ABACUS" outpath_soc = f"{path}/Fe_soc1_nscf/OUT.ABACUS" parser = AbacusSplitSOCParser( outpath_nosoc=outpath_nosoc, outpath_soc=outpath_soc, binary=False ) model = parser.parse() - kmesh = [4, 4, 1] + kmesh = [6, 6, 1] # e_z = get_model_energy(model, kmesh=kmesh, gamma=True) # print(e_z) From 3115c8915633351631d0df813907e644170f1f26 Mon Sep 17 00:00:00 2001 From: mailhexu Date: Tue, 4 Jun 2024 16:08:30 +0200 Subject: [PATCH 7/7] add MAE.py to abacus --- TB2J/abacus/{splitSOC.py => MAE.py} | 143 +++++++++++++++++----------- 1 file changed, 86 insertions(+), 57 deletions(-) rename TB2J/abacus/{splitSOC.py => MAE.py} (70%) diff --git a/TB2J/abacus/splitSOC.py b/TB2J/abacus/MAE.py similarity index 70% rename from TB2J/abacus/splitSOC.py rename to TB2J/abacus/MAE.py index 2d57f6c..1e6713c 100644 --- a/TB2J/abacus/splitSOC.py +++ b/TB2J/abacus/MAE.py @@ -28,6 +28,19 @@ def get_density_matrix(evals=None, evecs=None, kweights=None, nel=None, width=0. return rho +def spherical_to_cartesian(theta, phi, normalize=True): + """ + Convert spherical coordinates to cartesian + """ + x = np.sin(theta) * np.cos(phi) + y = np.sin(theta) * np.sin(phi) + z = np.cos(theta) + vec = np.array([x, y, z]) + if normalize: + vec = vec / np.linalg.norm(vec) + return vec + + class AbacusSplitSOCWrapper(AbacusWrapper): """ Abacus wrapper with Hamiltonian split to SOC and non-SOC parts @@ -65,30 +78,11 @@ def rotate_Hk_xc(self, axis): def get_density_matrix(self, kpts, kweights=None): rho = np.zeros((len(kpts), self.nbasis, self.nbasis), dtype=complex) evals, evecs = self.solve_all(kpts) - # occ = Occupations(self.efermi, width=self.width, wk=self.nel, nspin=1) occ = get_occupation(evals, kweights, self.nel, width=self.width) - rho = np.einsum("kib, kb, kjb -> kij", evecs, occ, evecs.conj()) - - # for ik, kpt in enumerate(kpts): - # Hk, Sk = self.gen_ham(kpt) - # evals, evecs = eigh(Hk, Sk) - # rho[ik] = np.einsum( - # "ib, b, jb -> ij", - # evecs, - # fermi(evals, self.efermi, width=0.05), - # evecs.conj(), - # ) + rho = np.einsum( + "kib, kb, kjb -> kij", evecs, occ, evecs.conj() + ) # should multiply S to the the real DM. return rho - # rho = np.zeros((nkpt, self.nbasis, self.nbasis), dtype=complex) - # for ik, k in enumerate(kpts): - # rho[ik] = ( - # evecs[ik] - # * fermi(evals[ik], self.efermi, width=0.05) - # @ evecs[ik].T.conj() - # * kweights[ik] - # ) - # print(np.trace(np.sum(rho, axis=0))) - # return rho def rotate_DM(self, rho, axis): """ @@ -149,16 +143,16 @@ def calc_ref(self): for ik, kpt in enumerate(self.kpts): # evals, evecs = eigh(self.Hk_xc_ref[ik]+self.Hk_soc_ref[ik], self.Sk_ref[ik]) evals[ik], evecs[ik] = eigh(self.Hk_xc_ref[ik], self.Sk_ref[ik]) - print(f"{evals.shape=}, {evecs.shape=}") - print(f" {self.kweights=}, {self.model.nel=}, {self.model.width=} ") occ = get_occupation( evals, self.kweights, self.model.nel, width=self.model.width ) # occ = fermi(evals, self.model.efermi, width=self.model.width) self.rho_ref = np.einsum("kib, kb, kjb -> kij", evecs, occ, evecs.conj()) - print(f"{self.rho_ref[0][:4, :4].real}") def get_band_energy_from_rho(self, axis): + """ + This is wrong!! Should use second order perturbation theory to get the band energy instead. + """ eband = 0.0 for ik, k in enumerate(self.kpts): rho = rotate_Matrix_from_z_to_axis(self.rho_ref[ik], axis) @@ -180,37 +174,28 @@ def get_band_energy_from_rho(self, axis): # eband2 = np.trace(Htot @ rho2).real # eband3 = np.trace(Htot @ rho).real # print(eband1, eband2, eband3) - print(rho[:4, :4].real) e_soc = np.trace(Hk_soc @ rho) * self.kweights[ik] * self.model.soc_lambda - eband += e_soc - print(eband) return eband - def get_band_energy_vs_theta( + def get_band_energy_vs_angles( self, - angle_range=(0, np.pi * 2), - rotation_axis="y", - initial_direction=(0, 0, 1), - npoints=21, + thetas, + phis, ): es = [] - es2 = [] + # es2 = [] # e,rho = self.model.get_band_energy(dm=True) - self.calc_ref() + # self.calc_ref() thetas = np.linspace(*angle_range, npoints) - for theta in thetas: - axis = Rotation.from_euler(rotation_axis, theta).apply(initial_direction) + for i, theta, phi in thetas: + axis = spherical_to_cartesian(theta, phi) self.model.rotate_HR_xc(axis) # self.get_band_energy2() e = self.get_band_energy() - # e=0 - e2 = self.get_band_energy_from_rho(axis) - # e2=0 es.append(e) - es2.append(e2) - print(f"{e=}, {e2=}") - return thetas, es, es2 + # es2.append(e2) + return es def get_model_energy(model, kmesh, gamma=True): @@ -246,26 +231,39 @@ def parse(self): return model +def abacus_get_MAE( + path_nosoc, path_soc, kmesh, thetas, psis, gamma=True, outfile="MAE.txt" +): + """Get MAE from Abacus with magnetic force theorem. Two calculations are needed. First we do an calculation with SOC but the soc_lambda is set to 0. Save the density. The next calculatin we start with the density from the first calculation and set the SOC prefactor to 1. With the information from the two calcualtions, we can get the band energy with magnetic moments in the direction, specified in two list, thetas, and phis.""" + parser = AbacusSplitSOCParser( + outpath_nosoc=path_nosoc, outpath_soc=path_soc, binary=False + ) + model = parser.parse() + ham = RotateHam(model, kmesh, gamma=gamma) + es = [] + for theta, psi in zip(thetas, psis): + axis = spherical_to_cartesian(theta, psi) + model.rotate_HR_xc(axis) + e = ham.get_band_energy() + es.append(ham.get_band_energy()) + if outfile: + with open(outfile, "w") as f: + f.write("theta, psi, energy\n") + for theta, psi, e in zip(thetas, psis, es): + f.write(f"{theta}, {psi}, {e}\n") + return es + + def test_AbacusSplitSOCWrapper(): # path = Path("~/projects/2D_Fe").expanduser() - path = Path("~/projects/TB2Jflows/examples/2D_Fe").expanduser() - outpath_nosoc = f"{path}/Fe_soc0/OUT.ABACUS" - outpath_soc = f"{path}/Fe_soc1_nscf/OUT.ABACUS" + path = Path("~/projects/TB2Jflows/examples/2D_Fe/Fe_z").expanduser() + outpath_nosoc = f"{path}/soc0/OUT.ABACUS" + outpath_soc = f"{path}/soc1/OUT.ABACUS" parser = AbacusSplitSOCParser( outpath_nosoc=outpath_nosoc, outpath_soc=outpath_soc, binary=False ) model = parser.parse() kmesh = [6, 6, 1] - # e_z = get_model_energy(model, kmesh=kmesh, gamma=True) - # print(e_z) - - # model.rotate_HR_xc([0, 0, 1]) - # e_z = get_model_energy(model, kmesh=kmesh, gamma=True) - # print(e_z) - - # model.rotate_HR_xc([1, 0, 0]) - # e_x = get_model_energy(model, kmesh=kmesh, gamma=True) - # print(e_x) r = RotateHam(model, kmesh) # thetas, es = r.get_band_energy_vs_theta(angle_range=(0, np.pi*2), rotation_axis="z", initial_direction=(1,0,0), npoints=21) @@ -286,5 +284,36 @@ def test_AbacusSplitSOCWrapper(): plt.show() +def abacus_get_MAE_cli(): + import argparse + + parser = argparse.ArgumentParser( + description="Get MAE from Abacus with magnetic force theorem. Two calculations are needed. First we do an calculation with SOC but the soc_lambda is set to 0. Save the density. The next calculatin we start with the density from the first calculation and set the SOC prefactor to 1. With the information from the two calcualtions, we can get the band energy with magnetic moments in the direction, specified in two list, thetas, and phis. " + ) + parser.add_argument("path_nosoc", type=str, help="Path to the calculation with ") + parser.add_argument("path_soc", type=str, help="Path to the SOC calculation") + parser.add_argument("thetas", type=float, nargs="+", help="Thetas") + parser.add_argument("psis", type=float, nargs="+", help="Phis") + parser.add_argument("kmesh", type=int, nargs=3, help="K-mesh") + parser.add_argument( + "--gamma", action="store_true", help="Use Gamma centered kpoints" + ) + parser.add_argument( + "--outfile", + type=str, + help="The angles and the energey will be saved in this file.", + ) + args = parser.parse_args() + abacus_get_MAE( + args.path_nosoc, + args.path_soc, + args.kmesh, + args.thetas, + args.psis, + gamma=args.gamma, + outfile=args.outfile, + ) + + if __name__ == "__main__": - test_AbacusSplitSOCWrapper() + abacus_get_MAE_cli()