From 0ce04d34a7ad9ecdd599771082ec5d9941624b9c Mon Sep 17 00:00:00 2001 From: aichhorn Date: Tue, 5 May 2020 10:29:13 +0200 Subject: [PATCH] finished SOC tutorial --- doc/tutorials/images_scripts/Sr2MgOsO6_SOC.py | 112 ++++++++++++++++++ .../images_scripts/Sr2MgOsO6_SOC_Sigmas.png | Bin 0 -> 26720 bytes doc/tutorials/sr2mgoso6_nosoc.rst | 2 +- doc/tutorials/sr2mgoso6_soc.rst | 101 +++++++++------- 4 files changed, 170 insertions(+), 45 deletions(-) create mode 100644 doc/tutorials/images_scripts/Sr2MgOsO6_SOC.py create mode 100644 doc/tutorials/images_scripts/Sr2MgOsO6_SOC_Sigmas.png diff --git a/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.py b/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.py new file mode 100644 index 00000000..fe07a9a7 --- /dev/null +++ b/doc/tutorials/images_scripts/Sr2MgOsO6_SOC.py @@ -0,0 +1,112 @@ +# Import the modules: +from triqs_dft_tools.sumk_dft import * +from pytriqs.gf import * +from pytriqs.archive import HDFArchive +from pytriqs.operators.util import * +from pytriqs.operators.util.U_matrix import * +from triqs_cthyb import * +import pytriqs.utility.mpi as mpi + +# Init the SumK class: +filename = 'Sr2MgOsO6_SOC.h5' +SK = SumkDFT(hdf_file=filename,use_dft_blocks=True) + +# Find diagonal local basis set: +SK.calculate_diagonalization_matrix(prop_to_be_diagonal='eal',calc_in_solver_blocks=True) + +########################### +# Now we pick the orbitals: +# BE CAREFUL: THIS NEEDS TO BE DONE PROPERLY +# AND IS DIFFERENT FORM CASE TO CASE! +SK.block_structure.pick_gf_struct_solver([{'ud_0': [0,1,2],'ud_1': [0,1,2]}]) +########################### + +# Now we set up the U matrix, first in cubic (Wien2k) convention: +U = 2.0 +J = 0.2 +U_sph = U_matrix(l=2, U_int=U, J_hund=J) +U_sph = np.kron(np.reshape(np.eye(2),(1,2,1,2)),np.kron(np.reshape(np.eye(2),(2,1,2,1)),U_sph)) +U_mat = transform_U_matrix(U_sph, SK.T[0].conjugate()) + +# Now we set up the interaction Hamiltonian +h_sumk = h_int_slater(['ud'], range(10), U_mat, off_diag=True, complex=True) +# And convert it to the solver structure +h_int = SK.block_structure.convert_operator(h_sumk) + +# Solver Init: +beta = 40.0 +S = Solver(beta=beta, gf_struct=SK.block_structure.gf_struct_solver_list[0]) + +# Solver parameters: +p = {} +# solver +p["random_seed"] = 123 * mpi.rank + 567 +p["length_cycle"] = 200 +p["n_warmup_cycles"] = 100000 +p["n_cycles"] = 3000000 +# tail fit +p["perform_tail_fit"] = True +p["fit_max_moment"] = 4 +p["fit_min_w"] = 4.0 +p["fit_max_w"] = 8.0 + +# double counting correction: +dc_type = 0 # FLL +# DMFT loops: +n_loops = 1 + +#for first iteration, add the outout group: +if mpi.is_master_node(): + with HDFArchive(filename) as ar: + if (not ar.is_group('dmft_output')): + ar.create_group('dmft_output') + +for iteration_number in range(1,n_loops+1): + + mpi.report("Iteration = %s"%iteration_number) + + SK.set_Sigma([ S.Sigma_iw ]) # put Sigma into the SumK class + chemical_potential = SK.calc_mu( precision = 0.01 ) # find the chemical potential for given density + S.G_iw << SK.extract_G_loc()[0] + + if (iteration_number==1): + # Put Hartree energy on Re Sigma + dm = S.G_iw.density() + SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) + S.Sigma_iw << SK.block_structure.convert_matrix(SK.dc_imp[0],space_from='sumk',space_to='solver')['ud_0'][0,0] + + mpi.report("Orbital densities of local Green function :") + for s,gf in S.G_iw: + mpi.report("Orbital %s: %s"%(s,dm[s].real)) + mpi.report("Total charge of Gloc : %.6f"%S.G_iw.total_density().real) + + # Calculate new G0_iw to input into the solver: + S.G0_iw << S.Sigma_iw + inverse(S.G_iw) + S.G0_iw << inverse(S.G0_iw) + + # Solve the impurity problem: + S.solve(h_int=h_int, **p) + + # Solved. Now do post-solution stuff: + dm = S.G_iw.density() + mpi.report("Orbital densities of impurity Green function:") + for s,gf in S.G_iw: + mpi.report("Orbital %s: %s"%(s,dm[s].real)) + mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density().real) + + # Write the final Sigma and G to the hdf5 archive: + if mpi.is_master_node(): + with HDFArchive(filename) as ar: + ar['dmft_output']['iterations'] = iteration_number + ar['dmft_output']['G_0'] = S.G0_iw + ar['dmft_output']['G_tau'] = S.G_tau + ar['dmft_output']['G_iw'] = S.G_iw + ar['dmft_output']['Sigma_iw_%s'%iteration_number] = S.Sigma_iw + + # Set the new double counting: + dm = S.G_iw.density() # compute the density matrix of the impurity problem + SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) + + # Save stuff into the user_data group of hdf5 archive in case of rerun: + SK.save(['chemical_potential','dc_imp','dc_energ']) + diff --git a/doc/tutorials/images_scripts/Sr2MgOsO6_SOC_Sigmas.png b/doc/tutorials/images_scripts/Sr2MgOsO6_SOC_Sigmas.png new file mode 100644 index 0000000000000000000000000000000000000000..7147695d860c912bf8d0f8476fe615f706ceb75b GIT binary patch literal 26720 zcmeFZbySt_+AcalK)R6z>2B#1X%G;lq`SKtq@_fpq@<)<>7Kv@Y3Y(ir9tUB&-Ay} z`o6WkHO^jRpK-?iXX{jk&UeoDsrR|B`?{|CiB?yY$3iDVhd>}$iV89s5C|d@_(P70 z0{(~lhXe-rf#fEw_yQID@<%m~0{=#HR?u~WKrry(e-OSd1zv;y6n2-@aewY)>F#CX zY5{RDad)Sp8Q_>h+q{APb>?e6X@!o~Gp7jQbcT5(0>OoBT-geb~9 ze&L_or6W{zsFNbS+!P@%p)dn_kb)H^x`*m>>qJh+J!TEjY&{u*=!rfBB zZN#sSn9phCBGZjf8Wn*S{7P3>GJxRW;f<@+2-ARn>zcSp%gD&oc5j76MMas152JxM z2hcMJz1MTrmT44^jojE5w&=ygD4Lp@G*Ta9OBYQj)*5c5v4@i~%8$)ZxHVw^FK<%)zuP*U z0^zw@3kz!7rQzkphj0ho`kPx>(Fq9Do2htu3TTze{Ca-SB(FsNqqbM4eW=9%qA2gL59+ zxvoll&n^d3%Da@5h)2;FOuHef{-{L7-H|?k9JVQKlqF|Dv=Y|zc zO=;25(Vh1fg@2s*+Rs)au&}UX^%WEpJmljefK=N|FfMv7IAl6W*LMFk)QM?rdh&%D zt8s$gI)DG1FnL6Fv=WL#%gKq&)%F|TU4#;Jx!V1w2D^SE?z3ml4$sbD5m@Bx&qNb@ zxnbr9X(q1p=rmz(nxpoAxQlFTZ7CU`2&0x9_fRAiazi&zFW)m|)1>^T%lLnLh5i5f z!MgWJ@86TrGci3Wo1zsKCjI>Rvq2w)m^W?`<5OvUeTu$#%9uEE!S{7av@vG!1K)U{ z8`eX%s=U0sZ{EENd-V#Do`HdejxK+Y1P2GF_3WFUkmoK*z5U!PrW}8K4-b#6@|6C( zG_zW-A#PzZaI@WAdv?QC5=gc8k*T?bg$F4CB_&pc76t?r3rp76_#rDB+psQ2qAbr} zQ^wMg@!gkw!P$bsLK+SZ zEXaq6pT2By#iPf+e)0JJ{=s?12Q%kNGa*y)bx^H$SeRZnd;It@GaK6z1qIYEU%udy zkj(BiMns?trSaTj)2{5n5?tgyw2emTw8oE$kC$_FbZkffUx5SvilO&zZ{NOsH8i9G zkplO%9DYA#5o=flGw@mpmyLN~vcK50RUTigHXqBlkyX_z>6zX2*2@j|?#s*!#pquP zc|2_@v!JNRY-hS+%=Ytay#x8zi3tP<*zXjcK1GK3{~nj0{cJ}(U7;6we(v1^gRQT0 zg>Gr9x`&U`J69KDJMD|b$=;=gx>Whc3E!Pl-PLWpY!*7Fo)qvZ__sC zgP~XqQw+_(Mf_g+{xD}k^IYvQ(P?HyU*Us$kcT4u7=kjTJS-YE@-YP zMMw4P+p8+1gy9Nlx}Ytn8TO1MCHJ?eW8^D#E@+06=TDxWaL{qfk_WHas+A}+(9&X< zxk=ZK+j5(ZrrX+!J`34n;>`Qhjts#iUkf$5qs8xt!j|dw2z2&bW5|t)q9Dy6c>CgJ zQ@~zQGo&hB%UJs{1LsbhvzKSLxdVhs?d#g~vk=r{m!xJDTH0C~0VZ>|IEwbre&6fNWU_2IxR{=2D$&V{Gp zqJ80Om0{fU5m#>5AV_0Hm1huqF2YbjEN&@}oR6`(eHS|Ddk=IdBH@&X=FFet9lEj7 zZ*?(AI+W*Tue{%3;o+q}ogR+9OO7b`7(E>Q4otI#Mgoo-=gPzr3u}A+hrFExqqnYv zmI*lPJ_=D${du4xNZ!p5NGaLL7>c6)(8I}DTdC77eE;e*#Q0=^&ZRjtzP^Z+>>iUb zEb49_&=Yea(346>juJ=P-u(HrR$QGjEh1oXjL*u zM@xI3X=khF?n_G2d&ah&WFd*Rt}w2I6qB3mU;W}n7~0ABe6WbMcJj(6-Q1Fn3;?lY zBzm~<_CQI|Wv?75YvW%>JZSAWLzI~HL(q(?AgMg<{M7(6-bccU(6_;kj*7BZ1YQ23 z7)?XPS?svsNU4l3{qVq_V>-YnibHZY!{r$HYES-ZyOkcS2K9An5S~!=@>)t zT`x1+gfGWf>4!Ei9;Pf3hLFn_3wI?0M5cT?94*<9`FrlY(`eq5v5ZduhE~vX2o1zZ zo9MwZNxbZ$#}gSDPib)3{2=nlu!A+<`<2~-cg+-YtmU@d*>zJo5&`s7ZT#q191hII zjbYEaBo5t^Z-V&n@SfeD3tBJ|$J66r>NT8*2#|HwVCD9Dg#uYw--<9vN~Tw{a%rNG3ecpvz=a@F zYr86wG%EbT{WbNJ^2}k@3!zKG@m%(|IvN|8_gc(3HH8(&tPzR+rr_^)0fd~YN@Ct_ zNMP|n8!nQ)cc+<)g-p`GrV*^WY{}2hub`?Lp;ad9>dI?tXXi(;Fg6B>lMA&NOgXWm zY~LRNEmwzfFwF+#K424|hJyWMpM>~zL%VpWQN=Abgqtg|#U16#17XHmMu1-8bcYq? zgK0a0l z;nHT}BbK$ety=2iHoc<4;V78IV#55OZ>z^0(z(y4gYDPUTaj@%YTY>pR7@`6r^8yH zoe%!jPFY!GX=!N)q^ql|xVZQSZ=X5uYPwZ+SE**1ufezX{B|H#f|t6k?wm`<{Gi3y z-n5MkpXuPZtJv)-g35VS#rdi|)~(+uNl*BtlEHgujb?mwIUI^l?1ashJTJ;?ZEYRA z0a0_;T`yPN-LEZ0aM_=rVL4Fr5!utK-v1VoMfSD*n%C4F45~)GkStDybvnL_Nendw z&ks6~tgXhrJ69iijRVC>v@IE|)4RC1AQAUz z2`C163_T2GwLKc@vX081&|MS!ESv`75Tsy&*LXIb9M7 zHp}WZbmCrB(^EBhQ7y=usz2((0P;2M*m8_Q@`2 zcpQk`&w{`>@J=oeNjS1@4E~Oy3>^=%pr(+rF~GR>TJ;XQY*IGB2g7Ahn_;v^gpQmS z+Un9Q9XK3~i&*9FOK~S`5H2}P8UuFdLQW9l`W@$nQwShqC_XfIYYeC~1)PORf?nBO zm((Nct1O9WsX{7$rf2tk06$a6x8jz?2v$-9Gppx@Jh#_VEZ~{#t?Qs*@DPS~a4h5HN)s}6^_?<*ilB&F3@@z%4IoOmY--v72 z&w&sp8?+^ZOKL5&my1KR6J`^=2M1tHvfk%8J|slDaY3AIRNoxM2FO6?3TKDMMWn9q zMCTtmDK=zYA2e7UE*oX{fUWyma}c7K&s%;lfibmZ986MI1GDI?1;VD^2k9)%54hA! zw`Or$>x2^FgPBM)Yv6`Y^ds&iTKsH@lOt3)BY}o<@_yA4!Kg4JXpEMr*K##wQB8;S zf)5iMX_u1gCGz;Ozfko&s|Q5H99GL!%FU6p>B6w+oVHFV&Z%MN*9LO@ zy|Oz*SG*BSx-V&|W0r-vVs?U5k&A;9hMBSDCg|S@kt@W=mOk!eU z=)wZ|>guXGv-hipZ>Tit1_p6!Yo?$B3I^CAr2AZ!I)Q1cw{I{HKlaS1mfEts>w{tP2-<))#liN(b349 z9BRwa?6|NnBzB{A3c386h0O}xIuwyCZ174>P8P4!CZ7%L87VnNar<(x=52HD_V|S5 z+SZ+$sxk43IOG(E+YRK(t7jc|;*byjL@xDE39?vOT1uBqffp?wU)}mGh_I)PP4g5Q zn`!n>USr-Jr*?$FTQHJ2_w)AR*Z2GnA3Do@o1vJZ0(w^GAWuwtWP`u9w#Mzbt7E^| zShLP!XVHi;9s^Jc9x?IErka*k%Hhe$-lOT!(W=2;lDFq(tLy7@BbmZ<^YimPyK{BD z^Plbf^eWQ0O{5bN(p4B;e%<-g?fYFt8m3NY5Br9uoo+WspG8*WHP=j;4)Vi})?tO5 zkG*|;-=AD7CNV~T`JzMeHw==9p$kz9W8miiqlqE*K#(zu%PEDOofSJ77>nXx=`<63O;Pw7fYJ8zue?j z1ajm(kF23;f6PA291#ztjql$vXu<%+g$FQBUVXh5>i8KZ#`ulj%#7HH@A&`-w7&lA z%*?1378a;zXpePubsMg@_{p%=*4Ob!NzovUzeYrlKpTyJel5J`SPWM3zS7vr$sRIu znUh>p)m3n+$5d&R>gtPLhSnB(yxVL}ORmwch8Mp{;!|IGDnGL?#sN#HQq9b3e=$1M zD3CNTFi^;K19fXSBa+9g&-i3}%G%KpAi4)~u_WU^g#om2_VAFEl7a-EJ{@J#Du+OZ zvP6joQ#e9#^CljYR<&QE4-539J|&ABO{x zQ=Sr6v02O%%zIKA8u*_3i*tn$14&HlA@`|*yNrXCraJQP(oBA0rC`d--+9l+j*Y1! zp`hf|*5ZRkT3R$p4$W62y?y(&W=4bgpvTbiAR?i>zJGgiqJTa-VqbR5UoR^`Mm-bx?s0_74m^wzXw7 z47xetT=zQo76(>|*CHniOINiC6aryQMsII#FRrZI|DBl=k1#K(uB@a4fjo)8d~9cD zm-;{}J|*^j5fTKKu7? z_n{0yoIi9c4t0Ej#3d@^O-4pc!#7_HyP3-Hq$C2ZikIQw89=4sHt9wZ`#tu=l9-BW z1#Rh&=G;7pC-5ewA7K3IMZ9mfxxBC{>oY$CBctc;3>A+l(#NM6 zC=e1pOH9w**>MYhU0u=$Y)ad4GQ@lLp1pjD{`sXtjfGPBdRG{!qZdDvfKjD-mN(sG zfyADalZ_499gYF;BwFDy%z<>e|8O*YDiiRc{=bF=4fbb#y6d!fI>7e8zDyzEKbxuP zVR!0{m3?)FN5Nr3GG?=xaUyqpnRAt^&V1N$rWV{1lhkYHPG4DBNqMKmnoQY{ti*tc z@+9?We}BKYs_LMT^1*|MygXX#M`R?JD3dcYUN*el-Q6M%*j${PeH;yoUtL-X3z7Ht zw+{NQSW4Vf>dj&>&mu?}-DX*#SXfwTJZ30bmHN@-j*Yn!X7J|DEq(R~wCzj7Z7hCw zK6ASkD5;Ga0+KHzra|(<`D<_eqAePalm{S!0KVI6@%HD(?@u!XR_hlWL#1=26cy1x z!UTbU_h2(dL2>i*FDIv`*LJEp`T&`W@aZlmIW?6|Ol)zhDjcOK#Jo}PuN2^u<~GA+ zQ;4}v)y#Tj%YJvBYz#Me&^QM_$wuro3m=_0%ekuH3Q)7kW{gi(eQ!V6)Ae`#fTYGK z>y{mztdh=)d>@-47#9%&r-ol>Ytyl^qJt;Q78accs6P=#t&6Lxb&$Hr06rR@fZ*t> zGM9%+{P5r)6tIBCQ;Iq=lKgKl@9hy79kgFL!Y6CY`qBAzSrz0H_3_UBe*DUc@j?d) z**lcCtMG-_;|mjK#NDWadjhHUo{Tr{p4r=G>_G2b8yt4hH(o zkn&v5jrP>yoMPmUKN`cYP{od{PTx6#Q4@nap=~P-43z=U1kP-9!A7{_ktrNmrAvkb zHIE3|p3*kjn^x3UFt#hGn+LVDWPH3w}3}_9pk9|m|B0Qc;`)=pZ?@Z2D zX~x$c--MWx#m3+Bivvkxw#Dm!AI7d%U#>Pt29m8D!7KvZ*9?5S-@iwOV-Ua0V+Xnb z40GYckrLwUsX^=EO7EFn+5?ci9|e=lRqwSHfRhvJud00m5dj>!JGS5g+bp2*ALu?v zV1LM>!wXeTVPADrtQGbfs(y0)PE3peUtW19v1XH=Hqw_ySfF@fJ@w`MVT)S)N*&Z9 z?E!h%3~_8^UDhMUybJ=dXQMsIH>LE{L(X^H&|@g(Sx%F1eS*^uL*``k^7o?HX5Icr zJrVbWo=pd$NCFTXD(lg>esZ@GLuUbJE#Np%|HbK1?b6~RfJ4JLWfJ((ADWxV0G|t( z25kUOpqPU=e2YOT{WVYQg%KcL#7wPi3X!@_yM0Sedm&{?(eDR+Q4x@?8cm+e)9X{O za(SfR+26;IGu!MWRO21i<`60!yd!}t6L0TY253FxTV-Q_N)&rl>IAu!mDL|)hB6+^1BFWyHbO7l5FtD?5<23NkGFY29$+23 zYj%tb045gM-vl&5>5tFfw6(T!+fFi<4E{(0fi}jdX_5@ZE$eS0+7I5VZbMhGv^o>4 zq3S)XEqk4qcfu)qnMu#`oocY{d=mHb;V7S|-7OqxuTGum7U57UeQ6>3 zM>j!X+MoXKq9EFqTEu~8${{i$3& z9vT{2ny}l)$wK|-1_qQ{BUv%qBxNHrPj?d4t0&RmWo>G+2L)QQZSFc-J#Q+Qe8ySF7Uus9nZUJ(NLdJLyMapNX zVDN>XX;hk*dU~=B4i!B_C)9lQlwjm$l_SE9<0~Fxlcp@CB>k5{$H~c= zCK)7ATUS@eSr!8bh*)Y#D#*&-UR+ul$-8&&6!bOEUdO~7U3|}>NIq z2ww}85g_2>hrfO$HJ09jG5^z<0eENL>+`a=Yysg52^xX5(^xTXZBG8F3>^c*Wu+5A zB~$3xwnp>GMv~f^TF@oIJEgfBFE3Z1nqg9k5Q4|+v9vUz#`EXnlaurwx;K58dQ^P& zdZ%e7!NT=OP3Ai72XC^h?d^Lqh20M4?W($9F!^oQ=JVaUQJ{KQ6T_y+F;RAxJE#-H zj#05Ggdl2aYPttIo`gfyo7~W^Ch;6kI7w4LhV&M3X?A;kG20#>y5DgZSTeYsB9|W`m;dv;!V1VEvmHSl zJ(}>b)&*ld_dEh&Yxb>ql1av>>rYudy$rAjajizP_tK&~e*bs|KHY838x8iO*G$>o z(%TxB_ON)$+UzifFkHYi)d~3=7XSnaaF_1&X+uYHV;&gYib&FWEWd$j+@Tu&ei~?P z3u}#1VK=!`Ysk-uBB=GFwXXZrGWz3CbWP{japM`y zR_9dJaVxUlav%W+B^f!nwI#3R_ z(FFLKGLw=(c|+V1j+YGoyBA>gayD}|8UUm&b6%6*KgNZEZpCho5`rMs9UsG_bKf+I zeSKY$~x?7yiKty@tFGA0QBCAXvwe3t}@cFf^V_>mAK>+``rFgm;Db z9x;IB>!esf$siCl;Y}Z#okh(bk5ZGkKuYB_kalz9g9}uC2W`ZW*p%IEXWwFqvVGB6 zlXp4OWdNi3#5msJN&U<;8Y-%kw|C>D#hXxOc6J#}&E(;k!GB1uW^?r(AZrCsHW_XP z5Xn+~d%_DdX92x@)xMSX&pXn8DxJDNb&yM%WEgu$w~$l~QvBb5dhTiegp7xf4cKxz z)-epUnTyp31^&g%tOCvi3PQMFVLK5Ky+Acy=CLeYeDd^qOuhv!B3xT~`pou6_swdl zhV*7ZSduhC5-#`q`8ydTZ&T9(Nilu1y%jltzQP2?-XK5UWJ)F zS=nDh60~AshGbUQ$^`n*jcDtmYcAi<49yTDOjbR>< z^edB1!hIRfa++F%R#Or~8ZFEp8Y<#7PTB@{F&lHI8$Uu5Vj%?d_y=`n=kiZfEyrv? zx25`f1C{|1im~sTzsgvth6Po=j<-Lh(B9n%WJ36QDe&jZ~gpoPz3x!@y-v*z6 z-Yjx=*tOSh_h~Jjjo>hs+VmN10&|rOV#|nx}u6pNpGCg8O zyGYYoPv-2L@pg;+DP!rf+_q@awbjKV*&Fv0i`GvU%(sUZXxFR12q;I-Ts3oRrp9d2 zm5_0N>_zlxJ%O3!ziIums3!ql8dtcpO?C|9Jrspe;UrGMe zp@PFfuApKR+8d!&-t0%wE~F758RUU(*C86?p$^N@MAXwZvABqLP7@|@L45guGJ(|W30x?; zAG3YfFZt_mApvK#4(1e5R@!!ob4Endk%vcG2kT@;f{)H(_seUaeWAtl!}Z$8GV!B> z<;Ae|LM)*d`jK?KCvXiCNV`v*&W});jl8Tdwf&^{dMIhy_fd}=mZjZ2+K1G+w9gyM zP^Pqq>oN$y1}0A(KnUmHVOou4$>YF$J)4bW$k;pd3HzT8j|o#+I0?MpSZ!~f&P+zZ~!A9$6jlB z#oGE7|Gt+M69sD}G*Xzj`a+#0MU6Td^jjpGSaJv_sMQWLrF<3K312ypt#ixJu# z65-K!w~b6LBya{ILx@5;$}!_UQ(;~}%!t}nE6GztVp~qSjx~JF7e;*-?-}7wSJGx` z^}0F3`-ND8lza`=mG7&uK&m|aCV z^m@>STHEnY9zeu%H^}4}>hi&4>N_!w^xhdKmc}F?TdV z8`p>`dht6YAyF2fWrl5k**oo6xZ5juA+=c$tT7&#qfrtJcJe}eb=LUq@-CpBK_Gzf+e(p} zwEE+0YV=}dV#v-Gz)zAxT2V8~l1a-EeQOqa81p0d@YfHYnXwo?Gr07h5abOBpL&`m zXhI?vd9wE-=T?K6@bJXhF((p#@ydGzb|Ln(0tMIrs;+rHGz~{h3lO}VkG19@$?f3% zUi5hwscjsoE%(L=Bs~K8KqEm4efDAR%Z6oW&CMI<(4R20QAQdOMyZV0U`5n46-366 zcv!e~hY^H3eg5u2UVw_xfUUIs%a>$dzkY>40By;&jK}06yTR=4wBTqUIm(mR=JFd; zjV!mXBgpM}4(!UWc_mrbMd!2&0j#6d^q=Ee<&}4VB02IL9z1Ph0h0|bmV!(sCML|R ztlNzS`uefe)zyx$e^3ukDpIcW6hhJ5R8Heb>nnkCis*d+nscMbu|RI!*Bo+hXQlrA;?bhDl11a zr|@+N|4cHYJfP3tW@W|o#t`3wJm4~nEiW(UG|S7+&z-VGfaK=pqGDpoSXdMamod~b z|IO$nXEBn*#>$F|kBtoWC741t>%|In#&|HNDAyJ^Xr6AkeUjKl3l*UqF`_7V}K zL09|l&C*R=A6fD<{rvfJw%LR2(@PFOTt@(A6C&m4$O+~R0YV+{8_#oZJ{+*+G!hb2 zM@L6O-iIG9k_mKC-BeuF$H!jnhBBNlF-5v8lx8&sZVGz;4!xOetU?72hZ%oPkXp!r zlvsm!oH9#wW}Q(b@AfS#aC@B;)GH>4Ta_K%q_g@5eI&DD%YZ-twgX-vGt)51%7jC> z+6gvO&_dLJ;Ur!89e@~@Kas`uvz7NYAScfTG^4|m$R8rDY?gMP;N9xH@up;qQ>SC> zxv%;qo9Poi77CogDBsZ6Bxo1+z+2wZUW=^bw|3;~o>L47+lg2yDH(5MVTG(eYF8IN zQlgcS>;>?wUm7@QcE$BV_n5d$VI?s_MggP{T)dTNKm16OS^)?F-PED$raag9=R+oN zHRe=Y#cbd)ZGBd8rx&tjW4U0rW?@ep-K37M8LeVGwEmC2l`|5=U~PD8)nsXD^qxz(>{s+R^_7;ML4`GZJ=81HeUxOiX^% zo5%~drhw=Rfz;4QcwY-}fL$J!cpG|WMvDPF(dM%JPkpV8&hf9ZbKqG)$zX1vh-M%f zTs91T3 z|G-8}aNpb)pwE-uD<4r+5-_AcG-OON;0ucZm$g=b2J0BghvpKZ;D*UfbsSc*mxn)d zHO#5mdlTRKG81;5tlq+{Q|FsMa4}CT%@oa*!Gx6e|NNq9R7QS+uPrS(GilqIlD~O= z7J|Pk@f?{ZxQ5KAm75L)2=uhGsd@@g1;FFo|N4dM9szwBxgJAh7^Kr2f7Y#@><*ZE zQsoE77cF}{ZmGO?1vvI&ab!2txk{AjX~C9P%bKi#P>5LUxRye|l2zr{{jt^bnyreY z@Py;g^ZDvX821nStUUQ_(u(RP?8_c(WFJZ~diV$f8m1%6* z%RL@I6oLdPJ2Fh65@avB>GKnSyR!eGVjK=wygj;_v1ULswfsrcgyjj_R z3hx8~WMOR2{hFK`;*fK*54ojNa>=NSeSj>V~(t!Vvhy4^}o`DOaF?C&^jU1%FEl+@u)X zWZz*~>Fc7U43HsicSgnmf<~tn1okyodiboV1Z;jf0n1CZf4AszOW@ish(}%_`jkh8It-G?JC5D>9hZE4jkUgtt3%9Rcvizw$2SH~a%G;3vds=f`Ij%rT zjP3;ga&P#T?{23gu~MlnS`l-7HeMC=8wbncqmkMx&PAi`w^98O@KKplYsgeT*v?p} z{*guT;?bwIoDu>alXKh>Wm;PD=b%6!G5@)rfsN|>#?6h=6ahMbTnzR*KuW6oZMbU= zXx6%qLgLWIANF#9A=uh1eHEiLz`+EGWemY_GS=)x3Hjs=B_u zXez`)biVM)8!Uk2W_{YXsDGxCnGkj5#N&-u`vWcu>fV~KD2qbUaVAg14uK2C|D7Ht z#Mwn-tvz~XRj@&hrXTin%{E#=33St5YHw4d*lXw7xt{E-{pKE=tUdFCt)KFA0bjBA z&v+ewmjX-xx6QgvW@si84%OzmuK-(+TDR%pclF$ahor$Y6)xd>!&D1Z&+vr3kEU=>4ACT{SI50&ed zDNO3JBmP5a4N*`=PN&;D9W=?90k}aJc*-8Ku%IC$BP%E=Nz2K}C72ocAISniqc?L3 zPdA*yl>)3DfQM0OhFGTd?Ps;aIAC9D{4W(^$Iq1js9Lk)`L`ZzQBE!{93mp3Gy#YA z!ew-|=z1GEdz-JcLeZXHRz{ABl(zoj_uiZ^ODo5c+8$2SNL?R;R|NExZx!U_A=1#$ zK!Snf@#4h`cTzLn^ur;(QRAsn%{ES_f2t`cg4Yzt#!qE~g;zm^1pWfv3JhI&?z0MU zbzlSI@=6HR7_eK>A%OSyw+M~r|FC#38HVi)3oiHUH7sMt$xYbt*x1-8C@bHC{BUnJ zgb|!*{v&?tZJuvO&jOz*a8~Ym-JD9HhU|tG!bK1DiOJUtf^sWAZN@>YX!TQmBLbDF zpWG?x`=7G#pFkfwS7(LOezh;!da;!Oh|@5`s>syT8UC|qeBGDs?vGsE+*$&L9tdL} zx`Gl3xbmfO%CnJK)+H3jHO-;e^MiYYr-r_SPsLo>@HdxfwrU}jVfe(O5&~8N0 z*H@(8aS0tPE{BDB8sI1UqjXMAzK9`IFn6f$cae8nlhRWKE)|!%`8K#vn^r^~K{>kk z%8$tq)+d8_&CV#FkIO4w|Ey@@q!BHtxTE09gxZYg)@BvA)hiS>jWv72$ighvj%sYz z1iVwO-0;UtW^eN~nP^A5s)rdHRX~s=Vyi(7j|2Xz0G|*E+FG%?txVi!%O&D_SC=A; zNu5dLIC)#OR6!kHV&nw0wPaw26aaQF*+fQZ32CUChsw@l_$nw zo;W(GvXbrQ3;Q^Y`cSZxz?8jux@Vxez*J?8yNYd2OAeSOOi3K8%=TGk4;Jy}a={h2 zq1QYJ?bfXsy(5b1wP;_2g24kkEZ#qZnyda$WQ-h{C>1!h^fGXg*P7JgKh#V?hH<^b zJnaanT7_T9f5%N6hWCrc;SI5b_zO_?V&DEd9M^kiq3%3~gDH_j0%=r*JbDLule;k$ z9<4wN|7K>9esU}h6tJa#S^Ir}l8?saBoT3TX4px&E?u*|U*oJPtO#$0YbLI;`IGU- z-yp&M>AUkjG&{`Smi-VC7I_iN4x+8hq-AP50a%CtbQGjNUwF=En-(6SpdC>?z$wAzIt=%#e%3kQ6wi$hklEXW5!b88X%(?# zWM<&NrSH?@Z-CU4(hZ0io(TX%rdp6Y>Z&b{SeDq=tsN%XZ4CHodqBxNTrgXHk0I?> z!T|y#sKZHY(V$T@GHN3lK9zHq&ZhZ_D`z=WFtAol%lmM0Bs}@nqyrpk9clAIR>E{> z$;sD#zRLYp+21gZC+cJ`8F#l3fUWF^4s&aOb<^Cr1X=RODaWQt# zf?=@kQLgtzyii$HLR3HgIU8Xj(x;Gb%#E|CNdvy}yMPBcoa_rlJ@)K9u!^gn_l=5T zkJD2%FR!S;0(y+KP0h_S%F{2y!&o53jEszEq`E!@vjb2{9w^p;UW=F4#^2qp(E}cN zH8hg3d&c1XSt#(9`3)Qzcm0tpE!}9%z&9P9p2`8P-oDb9%}N>-h~@c>jU*uF)kt-F zs;y>bMaz& z?iSoh=t5hM7Jh!s9t;C^FTXGFMxQN9xm`a8=@Ka^DLg(tevAOz1wdmQ-)n?d;anVV z3c0Kyxw^W-Wv&t1`$V@=zelCMak04W?(D}tq}JX@)$`OS$3^5ar<7f?5I1bI%4Ao4 zI*NIg&}Qz+ap zMaQKSlGk0xB;om!9ER%aCI4c=9BLM5((`;D&ZNg-IlX@ht$_OZ-?CfTR36}Q0FodS z*sc=ayb0bMOaE{)Sp8u(O3IaN%XxwvV!X?K(9umxRAY@hBPO?N0)dd02+T6MKtkFjN|g2S@ku|qE-z=lzCG)>580jnjIj<&hy>IM z;9z|HM?k*M?Yf}=iWWWQ#l%6$A|W9mW^w%Wwg23vx2+1Um_g>7fDi6cRVB{@R+>1DN06Jn# zu;b!s)nGQa$t&SP;V{r@WrJ6Tzkfg4xd8QV_<*ScV(&k<2eP)T<^L0lU zP(XnXj0J9i14?i&AVGR-%wV8=YIm{edBMkzr_Ce#-x};T6~9nOeqIeGxjB9?UM~2% zx!lf2_^O!XcGqexKY)1XWqcu4A0`G(5AeRT0p~u@7xjy-FfcGc)q!i(&ZiQkiM$hHU_ee!PX|oB{NIbqQ~*o8y}b<_ z`q5yD;9~{KHn++>)qMd-6Ed7P9=GU=xcq%lo-69L-!sV7fhj2|Ng?ct1YC1a;Nq$F zJ+XZxC`iY~hCwIi{w?`UP7aQ?wl)&dC+aU`R8-gyP}=m(I5jCL6jZzdl@rl=`6ZmW z0SdiJwiqdn*dYWH<=zX{K}uvUi+C%Bu+^k{46{C^_#G6KVjnEZ?^ea&^)Nt{S#ezA zuCuF6>3VAdN#nIpYLM^7q?VwVnwkQ&d-BeQowpq=&{5L`5_VE-G_;?yPk;PFIe#c1 zFl`Grb8Y>y^h-QWtOo>N;m0h6ln&E-k|=%~D^ zDmEZZjxKkqyhuBM#tWqScB36&69l}jWMG%=BQhY^L#*#iUcY|b4F-0|k4BQtBHoEC zR}NVFYjILAe_Y(L6nDQKhB+@3OBLL)_U-yoWO48=W&kgn+1dWm0|Q3MnEN<-R*luw zxWHu82_#z4{W`gbppkMT$thd@b^swc&vx8J$_-lT*FAv;6`NYZexwvdTv(r&frqCi zKjYt;%F=t!wC@9#-ltE5e_6;1hiX6(65u(H9*-bD`t+#_R^epzoKV^rUZ!zX1&t!x z4xjq9ZXr4Pm+>s62R_R7q>le2{1>eU+L@n!_jaJ%=F`@mi4S5+{t~~SREUofnmg$2 z<#iUA6UOu}^}B!bi0$8S7))FCMoK$zn+zr%OT=-NHpLdEH+gT|+@O-Lg?Leob+sfs z9X7TC9K42y*!)*8`gTSt80`T*}i|0&DO@YetQzDL%-@3WL=!?m?5eK6$){T66#p4$O zYZ28P#G-gi+GkIHbG;D}#%E=^Mb(>-;>CnuJYerEP)u^(o@4%x9$;DvCpP@yktuFD`;6x=2bqF5GtV=pY`=YB zAbOm=H-p=Zi!4Bpg7c^;%*viXF|YQl!7fQ}(rN#R8;0A+MuQ2e#)O{Fg-s3xpq_ZZ zni--&rj<~^q00*eVRBULS1iL9+HF3+A#TY{b5&H68R9H3;oFa_dMJQ><)=q1Zc2-4 zT!n#yI+ZUM@t`1O-|WnXmHy*mY|?ui83iSP|)FKE-x=)uxh(z43Ld*tQLi6A&(_a z&zgZ1-X#oEl+3d6jy=&#fQk3x|8#Yw^$yRWB`mIy5bv7-IDp5ytMHch8b*?{K#d6E*4mGxq-aLIElP_uX z>ST?`aw}DMIZUma$!%mHM5~PdZnDmGQwjK@sqmdo>bw*KQ>+GEmhReX=tX16)l?;$Iqp8Og6!0#zB625#hXIa-qi zy2nC;6ETPjpLz(MQt*j+q?s&Kr($AsQK%sco9lRLGia4f%>wQRe%1zX2@Gcl(y_28 zZr4PBnhzMj2XK!YkPsn|?4Vm8Q2C+1Hn6(@&{ix^V;)z{(20nURkfYIpQL+A!Y9_o z1Fi3+(95CV9WpXf16I#k+7>hTF(}r)W`8@H~=FmtsKF($9;I zVgx2_;LB1IPazOJ2K3uX)Aww0S0fb_OyE6acpqhL3XZx!8akRZ6RsfuhjDPHtAO%X zIS&th@Ff_=hp)8J!==}LM)mc9BS?TtD|E9*HB;1U#J^+@;HM7&bwHd!vIeh70CfO5 zwdQDW1-0C;?d6we(Q?4r_4f~@sW7^Gcys|fU4J6okKJ7vCnqjoFN^wYHXlr|0FBltpC= z&Y!Iy;=Xjg_hS#8wP<)d?R@xRPqr3kmRS)`1jJ4N>`rBRBwNVh9 z;&LIp;D>UUaBQjHx)}YHEsNPUYql8FEtYZ>kmjM`bS2nj=h6|Bn9kegMkGfy1Y#7)Ov)%)%CQNk4SibTX4m%ifAP# zx#Gt(IRZopd1b#mqAUqLE8yOfFFhWrtjH*_PmE4bPd&Q~>R#ckevf!=^k~)HgFKcI zXgg`E2O()vYV7GmVa{>n-TaP?i{p_qk|lDdiKEr)S*GP}7mW7J{OB|imAKL?lf+MF z{Yoi_xx5}O&QRjDJEBfBb`5G`30e^?X1@3_PCh}iEA zDs{ftq|sh)v$kY-mv#5Jywej+Wt`i+eyhg2GisAx4$o$t?dZYH$H-OA2lQlxd@HGx z2|^Q%reildu5j;~Tsgnuhr*7k<4J z1uEc(5a=p2@@OLGmJp0qL}3PF?fa51H3RhO3)3<8w<67+eUfyHKb8+Js=FYvl4sBm(DWWdo|Ii&v>Bx zA$_pw>f}hm*>ILy+eCiH z?DV+dTX5ArmaCUm*US=WHm5D{#TUQpWFlHR|L}>_Z+Zk;!GrLmer`@Mcsqqk-*&GH zzNAHb?;n1TB&Ok<#c=5|92-9D?YP@RJPpzbxuMoFwlPI*zQj5&E+<_`8ScS3h@%*@C@$vaIb9P~kf3W_Z`4aj{93j34Sm1f2mjSIo)5FNII3 zyaTy_yU#_Kg_&8(-ku#~{0Im}>dPiSqh!q9ms0q(2-6%4IK^K!YLrM7Cy+}AI9{XoydrJy*ZCmZ(G?YMcn`5n zzwg8a_(#}de?z+!Bqu+eu5`*Ld#~s4qfRxhh2Q_acKeQ7s3hVR1_YYrU^H>xdt&;N zbPsZlKXUrtjoFk-st4x;+A%sFyTR@k4xd+lLgTM!+wuUD&{}UpV`&#Kbavq@wB9?# z;_Yf9h8>_0CwQwHIwEl|Fp1}W*#kKUBo9%{c_Xr6?d%{q|5}Mye^N=PEZJFuCGz@v zg*s+RVpjvVHenKOFp0n2Bdm5J)9Nvjr@`aqm&X+(?9un+2vnpyjq$lj!kdsz$YZ(F z+FQjNu}Qe!qs&BCEjrGMKa27! zpOsA=9vwYaR>nX?j0kMn1I!Sdih(@h=hqb>8*iLbjyITOMe>ciW*yL>yf9fimY!I9ED3Dls{eFxLBP^K7Ch6WU6#F4dpXdQvUE zm-WJodArQJKX<#eNnuqkzxP_-n^TU6?t7wjzsJ(Tf(B3;@UwoZZx$5s)nwwO8{ifIA(YtpkQ>Z0AJC|C}2f9S~Wr2v{6ILqJwQ!9zDVAE-AVlb-GGCt&q* z36xIR0^;vJyzp<_x(AdCCW1qlj6k6RAb0=+1Orf{I`GD~r?0PjbQC8Nhq}IRS4Qh^ zzL!MiF#W*rjTzu<0@b=o;(tHO3Vzf8YF`G^K0PA{l)5JW3C{SdYU97I%zJi1uHwCt z?4Q{*$#}y2D6SfKgvX0BU!W5Ty03Jo&M1kaKc!;_$GbgZE&P9~JM(C$`?rq|l4yjK zEE!wUElVQ%l68uZeJO)RD6)+$OOvw8(k)D6tSJ%UCU{8i^Y#~dQBsze zA-rVHDc3i(cdrR947?JWUWxLw7F7%0ji{#hs~|Q%#Jg2ST8507s(awxutYHqDtkS zFXR65oDIwO@iFDA@8h<#mKi)g5$0ws6YMltSmkxUCQR~CYlfWMR{NSg@(pCVx0t6l zipSiESzh73aA@_R-|^&qeETdbR-QZ0YiSZ>gi|UEdIqCkBw~zTduCyGCpauXUYz4j zK@)vMbQ!dQ(-a`oV1o-4U>gEcKVOjBXE_%sE68pT(qJ~6?TgtQ&Qv4c6#S6*i!+HqTQ z?+X><^l=!JkvZbZ)_ZMsXbH5WiGJBP3QJs36S_lCsc2Xj9cwLw_@whvNxd_Me#u+Y zGo9|eQ9bj~k%({SnKbCoY|N$FcZBAPkwsI!z1YdVwc^aZvt}Sj#z`B^xaUOw`a<={ zYFgt4{+M>5XGV4r>vgj?@7{fJz>Q}c@(S0|&}ZPvUG+{~qgG;J{TA{3Jd^L6{yyUB zuXCaRYZ8`jKDc@!jCYF1#MjqdDq!$QHl`)xuT^wys^oCS9}T*9=8FT4+MP9rg4|IJV|6&nN4esGI*-|CQN>s`=IuPl zFrJq1O!RqE9^t0d21-@TM<%DZ0g=y-oYnA z)<+5~RAmcWv+u@R;?nQb^Sd5Rd|(~TZvU4W#r+ovOvSBjoL`Kzq#8o0!O5vc~Jf-eEfx@#MUswPK3D7x4z5SyNr zMFY}YXfvTr3vsu)-vLesgn{vADHc&#OLtn1D*DeUvT1EVZB7}RdKAe_Z7t6ADTstQm#sy6{>Eg7luR<`TQ3nk) zIIJVe%A~h8`f4!6Ibhj>u@y`yz;OHa9_>N!1)`(12$ZSbX-_b_JehOUP#-)661xE! zIfrrxUDc(tG%(>Dd&x6MEu^q+#RY-t>hw2tJXS?n%(ZD)+39>^DUEp{!OM{U>Ue0% z)_&sC&9OC<$m15nvlipqWXJB)LcocW^_n??vKeoRM!f4`fvV6&OxeYNM=L#@7m|1M zI<2bRUZ%;FJPtvGm`}Qo-3{2UUu6-TWcuv9;0Gg{&HWE#*J6d|cGIsln|LQM*gz>O zJ`=Xf9HE}uTS)bZoM&t0a2 zUZ=21&FewzYnrb#ygL7lA~R0o?({(AUZoZDx9LcW4BUZ5&xg)9&YSHX9;-t><3A?a zDD%TGmWcy(=UAvX;hdDzUZ=mMCCE_jht3@s%{8~SQlD9S&-AbqI{UnGa<}`p;gdh; z-@GW`(^Tttb~I%7k08U`TGz9a{aZ=l91l(&S0QrjPHi&;nD!~pE!~#6@ybb;0fNP6 zPqAc@MB7Ci?#uI7S~|KcknezG)zQ)M{N>B<=qylR> z)(ln&;9X*x6Qslr9i_W}(Z&q-2dBbDTws?7goaBkWhej&A_qqY@zT+&9VSPN%Q7>n zm<*^ezB@q!51v5csIz7fHvh--+k@GR4VY7jBY>m-TU0&~F^49}zB2d=MI z13m+>Hm|Tc5f zzl8tI!kOL0NP){?nfvZXS>-;H3<)ZvHOIxc%rgtlMP6q9-S57HwS5nS1K?=B;kB0h z^MeXP=;836>fv3{Y4O=IUlc8j^>Pl^@Js%5Nk~4i5Ij3OL++=wE=DnpHrd zs^|5>7m+-(b=Zfu9~hy?u-s$^{2j+$s6Y;*ec~ zV!4P1^OTBOavPvN#?vo)s`iet+lj`g4kU;@sq7Djw2Yh()YR$V@bpf>Mq}YF%ePTay1ou??zzvORy$FhCFxMPF^_RN~#9)-< z$KF6sBCDu~1y)_D9YMqB-fu!~79;>EX=#2bsUs*DoZo9|YO-~@bdr_qH}gm3V*Nq9 z@^E$F=5K`AyN@63S{yaToBq*-KYwfK_HUsBgyGlv-0z*W!@pC1qa-*)CU%NFwNBp6p@u>2RkuBHdR-r0SH;>%k{PK*rm}h zp0NZooT{1Gv7dfQK=^Vk*MHRf3FsmQ=H{|vcqIDUlXBuE+FGI!AsB4pfq?-4WYmEY zOvSk)ZRe-Z@-P!@BiLhiB;xYV`~yp6 zzjX76s+&pTaH^inj1EbIRdCO>fDqv7e3u)BnRQ^j71w(JTkrqDLuR5MigT#876qWb z?!G>ul~dl7o}QzYFWL0;^bBLwY(jsrHz+Qa~_7z~&A>thHv_Q^?6HG}pU89wdh+hwocxx2JL1y2>7Dbf& zB-k6swRCheHbq01>`EY?NJw?)|*+k6G zOmboo?NBb38Q@K9Rl9&L;ZIIZ27Ly+t)Zun572*nClbsW0BIQdQidJPFE4NTC@Cpw zXXlQDL+5e#ovqujL7wTQ9)?CJ-i46em6DWnj35#;O|peQ-Q8Yo1Xnn6WpSG@3)Zl- z{6zX}7dnBLw=b`*#ywFEv^>D4fyIVF9F!#TU|d5(*w+b=j-qdu>HT|&!0x{L6i}1e z+6;T$tKs!dpT!oTG#fbPvl$B)ml{3OvqeN4KoJ2KKBCA06$=-Pgtef34#6-I3M;+^ zBSS-KpnC1ww-1F%6tTE=#~FB8zS}E(VYi-vZ@Z|dXtENq!ay;2;nJlBZlA%Rg}?>> zGwkTgc|zWq7GPzE4mvjAfgBm*oM4MGuk>9Cr%UMs!m+5RXnDZ;S-?V`lv06FeXEG-|^EA`%~sVM{&AfVhAL__D+rL>B!sDVp>;p*gP*z~Zx{y?42HAdt>Fv|cCQeF8J z<{-hD>w_Q#{}^Zt0C|SZFEP5x%_s@Km7tiIn8^D75^~i8=`$kRMbI`1R(er zR4w#tuv>o=K|KinI%OX4Mtu8z;f~*IAIhs0AeTTtipRqmhrUpp0UG`@Q`AK$!^#3PtxAxkms1 zb^yeX;<~!BQrFSZF%~j#?Oisp+7M-Ek z>BIwjf%;0BGeV~z2mENU?Ya?#E_w3qkwvAiL?>D4Am2We4G?X<&CJx_XHtd=@(5uX zC|S~l(3AXU3>ONBwS3D`hKm>PtGt9NpuVHymNrv8|Hg|JJ9ZSt?;GEsN1N!{Sd@5r z&rO7tfx52}otTiYWX@FlxW!IZMTG~^1pyB0A-D+P9rlzM&ZJTI)DHqd(|0zKnSaQJ zbQr2^bjN~pN64jo5~SvAuoE%Mf8zp{8r@Sj99yQHdZmX5?Pw;iugxq=%W1m-?OuDx z9rb7#irS?ms2LO22D^B>y+0UyW+t(; ziuDf)F&;@5O&t*Mh7H(3uzEb_U>sd_VHPj@UrR zz}F0^XuR4qj1YVSCyOq|6St(@dQB zy0$b9IGK~z@c1IalwR&z4g+0XMhK+$Ak_in-WK#*hOwUNFdai=a)GNuP-w__uFRO^ zU;X1u~{dUtxq=;y;o7JJbV&KD4AnvNzx2K4MIL%XtzC;^nC> zIuu}T!I41)bK1lD=Zay=Z*4CSK{}Z8F>CO$zrTX=+6ayFKrt<9WNr>x69P-dtZ&g* z0KIbM*PjR?Mj8+?-V@pX9dJ#{R4WQzykJ1KKA0~=R!#|Q08He;@YbZvW#Q|qT}m3) zuM2{F6Ir_mbPYo2xB2-*zkGVYuIxc=+7YPDeSq2wcg3Z`1_$%Xey}aly1Kf;_b};7 z+~k7ld2X!I4|ay6L+8jtET#@@7tnF%n++U_xyMn4u{W1*FY;;@q!gbPg5N*S&6OA~ zc==Kj)XBlD(JnSNC=mXdzjon9iiH(ffmRsOLwA6ggChd8vm;+y;@9SVYEG}13%tKs z3XuKDG(QAp2Jd!3v8O?T(%XpsOjf$i4wRpZ(i~;$24Z-2c6Q{p5#SB358&;cN+XJ` zm{VZHT?8%oBClD~Tr9_H&=u(S(%I7UA*K{DF6b#DW4? zNAYo>|N2Byk;{-g!kxil$4mQg@qXFx0YCfAotXt6#`^3mz*{ev+WE922T5I5qSIh50=MK!mhoDBR(-E%|BlMDPs;!Kb-3LZ?y- z^n#2J({Tw2AqXfFF4$+#oes75RP#?kI2%zx4UzBV%@4=Asm@f%VLg3)o*%w|^n;UI z_u$|Iz&l-u#Y!Q-Z{(@Stwa0{K(xu6{h!e$>> print mat[0]['down'] - [[ 1. +0.j 0. +0.j 0. +0.j 0. +0.j 0. +0.j ] - [ 0. +0.j -0.985+0.j 0. -0.173j 0. +0.j 0. +0.j ] - [ 0. +0.j 0.173+0.j 0. -0.985j 0. +0.j 0. +0.j ] - [ 0. +0.j 0. +0.j 0. +0.j 1. +0.j 0. +0.j ] - [ 0. +0.j 0. +0.j 0. +0.j 0. +0.j 1. +0.j ]] - >>> - -This transformation is already stored in the SK.block_structure class. The next step is actually not needed for a DMFT calculation, but it is good to do this check to see what the transformation does to the local Hamiltonian. We can calculate it before rotation, rotate it, and look at the 2x2 block with and without off-diagonals:: +This transformation is stored in the SK.block_structure class and can be used to transform Greens function or operators. The next step is actually not needed for a DMFT calculation, but it is good to do this check to see what the transformation does to the local Hamiltonian. Note that in case of SOC, although spin is not a good quantum number any more, there are two blocks of size 5x5, each corresponding to negative/positive :math:`m_j` values:: eal = SK.eff_atomic_levels() eal2 = SK.block_structure.convert_matrix(eal[0],space_from='sumk', space_to='solver') - print eal[0]['up'][1:3,1:3] # prints the 2x2 block with offiagonals - [[ 0.391+0.j -0. -0.815j] - [-0. +0.815j 4.884-0.j ]] + print eal2['ud_0'].diagonal().real. # diagonal of the first 5x5 block + [0.071 0.131 0.608 4.572 5.128] - print eal2['up_1'] # prints the 2x2 block after rotation - [[0.247-0.j 0. +0.j] - [0. -0.j 5.028+0.j]] + print eal2['ud_1'].diagonal().real. # diagonal of the second 5x5 block + [0.071 0.131 0.608 4.572 5.128] -So the local Hamiltonian has been diagonalised. From the other entries we can see that the *up_0* block and the [1,1] entry of the *up_1* block correspond to :math:`e_g`-like orbitals, and the others are the -:math:`t_{2g}` orbitals that we want to keep. So we pick the according orbitals in the block structure:: +We see that the orbitals are ordered from low to high energy, which makes picking the orbitals around the Fermi level quite easy. We just take indices 0 to 2:: - SK.block_structure.pick_gf_struct_solver([{'up_1': [0],'up_2': [0],'up_3': [0],'down_1': [0],'down_2': [0],'down_3': [0]}]) + SK.block_structure.pick_gf_struct_solver([{'ud_0': [0,1,2],'ud_1': [0,1,2]}]) We can now look at the final result:: - print SK.block_structure.convert_matrix(eal[0],space_from='sumk',space_to='solver') - {'up_2': array([[0.156-0.j]]), 'up_3': array([[0.156-0.j]]), 'up_1': array([[0.247-0.j]]), 'down_3': array([[0.156-0.j]]), 'down_2': array([[0.156-0.j]]), 'down_1': array([[0.247-0.j]])} - -We see that we arrived at a structure with 3 orbitals per spin only, and blocks of size 1x1. + eal3 = SK.block_structure.convert_matrix(eal[0],space_from='sumk', space_to='solver') + + print eal3['ud_0'] + [[ 0.071+0.j 0. +0.j -0. -0.j] + [-0. -0.j 0.131+0.j -0. -0.j] + [-0. +0.j 0. +0.j 0.608-0.j]] + + print eal3['ud_1'] + [[ 0.071+0.j -0. -0.j 0. +0.j] + [-0. +0.j 0.131-0.j 0. +0.j] + [ 0. -0.j 0. -0.j 0.608+0.j]] + +We see that we arrived at a diagonal structure with two blocks of size 3x3, and we have picked the orbitals around the Fermi level. The interaction Hamiltonian --------------------------- -We now set up the interaction Hamiltonian. Since we want to rotate the interaction matrix into the local basis, we are using the Slater convention for it:: +We now set up the interaction Hamiltonian. Since we want to rotate the interaction matrix into the local basis, we are using the Slater convention for it. We use *l=2* for *d* orbitals. Also, for SOC calculations, we need to inflate the resulting matrix to size 10x10:: from pytriqs.operators.util import * from pytriqs.operators.util.U_matrix import * U = 2.0 J = 0.2 - U_mat = U_matrix(l=2,U_int=U,J_hund=J,basis='other', T=SK.T[0].conjugate()) + U_sph = U_matrix(l=2, U_int=U, J_hund=J) + U_sph = np.kron(np.reshape(np.eye(2),(1,2,1,2)),np.kron(np.reshape(np.eye(2),(2,1,2,1)),U_sph)) # inflating the matrix + U_mat = transform_U_matrix(U_sph, SK.T[0].conjugate()) + -In the last line we use the Wien2k convention to write the U matrix in the cubic harmonics. Next, we want to rotate it, and pick out only the relevant :math:`t_{2g}` orbitals. At this step we need to give the indices of the wanted orbitals in the *sumk* basis:: +In the last line we use the Wien2k convention to write the U matrix in the cubic harmonics. Next, we want to set up a Hamiltonian and rotate it into the *solver* basis:: - U_trans = transform_U_matrix(U_mat, SK.block_structure.transformation[0]['up'].conjugate()) - indices_to_pick_sumk = [1,3,4] # these are the relevant indices in the sumk basis - U_red = subarray(U_trans,len(U_trans.shape)*[indices_to_pick_sumk]) - h_int = h_int_slater(['up','down'], indices_to_pick_sumk, U_red, map_operator_structure=SK.block_structure.sumk_to_solver[0], off_diag=False, complex=False) + h_sumk = h_int_slater(['ud'], range(2*(2*l+1)), U_mat, off_diag=True, complex=True) + h_int = SK.block_structure.convert_operator(h_sumk) + +Note that we needed to set up the interaction Hamiltonian first for the full set of *d* orbitals. The :meth:`convert_operator` method then takes care of rotating and picking the relevant orbitals. Now we have the interaction Hamiltonian for the solver, which we set up next:: @@ -157,16 +154,16 @@ Now we have the interaction Hamiltonian for the solver, which we set up next:: # tail fit p["perform_tail_fit"] = True p["fit_max_moment"] = 4 - p["fit_min_n"] = 40 - p["fit_max_n"] = 100 + p["fit_min_w"] = 4.0 + p["fit_max_w"] = 8.0 The DMFT loop with automatic basis rotations -------------------------------------------- -After these initialisation steps, the formal DMFT cycle is very similar to a calculation without these basis rotations, since these rotations are done automatically, once the :class:`BlockStructure` property *transformation* is set, see :ref:`basisrotation`. +After these initialisation steps, the formal DMFT cycle is very similar to a calculation without SOC, since the rotations are done automatically, once the :class:`BlockStructure` property *transformation* is set, see :ref:`basisrotation`. -The DMFT loop itself looks very much the same as in :ref:`SrVO3`:: +The DMFT loop itself looks very much the same as in :ref:`SrVO3` or :ref:`Sr2MgOsO6_noSOC`:: # double counting correction: dc_type = 0 # FLL @@ -177,7 +174,6 @@ The DMFT loop itself looks very much the same as in :ref:`SrVO3`:: mpi.report("Iteration = %s"%iteration_number) - SK.symm_deg_gf(S.Sigma_iw) # symmetrizing Sigma SK.set_Sigma([ S.Sigma_iw ]) # put Sigma into the SumK class chemical_potential = SK.calc_mu( precision = 0.01 ) # find the chemical potential for given density S.G_iw << SK.extract_G_loc()[0] @@ -186,7 +182,7 @@ The DMFT loop itself looks very much the same as in :ref:`SrVO3`:: # Put Hartree energy on Re Sigma dm = S.G_iw.density() SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type) - S.Sigma_iw << SK.block_structure.convert_matrix(SK.dc_imp[0],space_from='sumk',space_to='solver')['up_1'][0,0] + S.Sigma_iw << SK.block_structure.convert_matrix(SK.dc_imp[0],space_from='sumk',space_to='solver')['ud_0'][0,0] # Calculate new G0_iw to input into the solver: S.G0_iw << S.Sigma_iw + inverse(S.G_iw) @@ -203,12 +199,29 @@ The DMFT loop itself looks very much the same as in :ref:`SrVO3`:: # Save stuff into the user_data group of hdf5 archive in case of rerun: SK.save(['chemical_potential','dc_imp','dc_energ']) -The only difference to the other example is in the initialisation of the real part of the self energy. We cannot just take an element of the *dc_imp* array, since this array is stored in the *sumk* structure. Therefore, we first need to transform this matrix into *solver* space, and then take the appropriate matrix element. After the first iteration (here done with 24e6 MC sweeps), you should get self energies like this: +The only difference to the other example is in the initialisation of the real part of the self energy. We cannot just take an element of the *dc_imp* array, since this array is stored in the *sumk* structure. Therefore, we first need to transform this matrix into *solver* space, and then take the appropriate matrix element. After the first iteration (here done with 18e6 MC sweeps), you should get this output at the end of the run:: -.. image:: images_scripts/Sr2MgOsO6_noSOC_Sigmas.png + Total number of measures: 18000000 + Average sign: (0.884535,-4.11253e-06) + Orbital densities of impurity Green function: + Orbital ud_0: + [[ 5.20045070e-01 -8.24863778e-10 5.95348202e-12] + [-8.24863778e-10 4.30734642e-01 -1.29359496e-03] + [ 5.95348202e-12 -1.29359496e-03 4.80477133e-02]] + Orbital ud_1: + [[ 5.24181422e-01 2.22991244e-09 -8.16290063e-10] + [ 2.22991244e-09 4.30431196e-01 2.19004569e-03] + [-8.16290063e-10 2.19004569e-03 4.77161009e-02]] + Total charge of impurity problem : 2.001156 + +We see that there is a small sign problem due to the off-diagonal elements in the hybridisation function. However, this sign problem is treatable, since we have rotated into the local basis where these off-diagonal elements are minimised. + +The imaginary part of the self energy matrix of the *ud_0* block looks like this: + +.. image:: images_scripts/Sr2MgOsO6_SOC_Sigmas.png :width: 600 :align: center -The two :math:`d_{xz}` and :math:`d_{yz}` orbitals are degenerate (blocks *up_2* and *up_3*), whereas the :math:`d_{xy}`-like orbital is different. +Plotted on the same scale, the off-diagonal elements are very small, only the *(1,2)* and *(2,1)* elements are visibly different from zero. -A complete python script for this tutorial, including some more input/output, is available (:download:`Sr2MgOsO6_noSOC.py `). When running the script, you will encounter warnings during the transformation from the *sumk* to the *solver* basis. These warnings just reflect that the off-diagonal elements of the full Greens function are not zero at all frequencies, although the local Hamiltonian is. In that sense, we still do an approximation when restricting ourselves to the :math:`t_{2g}`-like orbitals. \ No newline at end of file +A complete python script for this tutorial, including some more input/output, is available (:download:`Sr2MgOsO6_SOC.py `). When running the script, you will encounter warnings during the transformation from the *sumk* to the *solver* basis. These warnings just reflect that the off-diagonal elements of the full Greens function are not zero at all frequencies, although the local Hamiltonian is. In that sense, we still do an approximation when restricting ourselves to the low-energy subset. \ No newline at end of file