From d61b55f0a430e01fde0003a83174b7834d1d1a6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Malte=20Sch=C3=BCler?= Date: Fri, 19 Jul 2019 13:33:11 +0200 Subject: [PATCH] uncompleted work on csc NiO tutorial --- .../images_scripts/NiO_local_lattice_GF.py | 2 +- doc/tutorials/images_scripts/maxent.py | 2 +- doc/tutorials/images_scripts/nio.py | 2 +- doc/tutorials/images_scripts/nio_Aw.png | Bin 25036 -> 25166 bytes doc/tutorials/images_scripts/nio_csc.py | 201 ++++++++++++++++++ doc/tutorials/images_scripts/nio_csc.py.rst | 9 + doc/tutorials/images_scripts/plo.cfg | 1 + doc/tutorials/nio_csc.rst | 29 ++- 8 files changed, 241 insertions(+), 5 deletions(-) create mode 100644 doc/tutorials/images_scripts/nio_csc.py create mode 100644 doc/tutorials/images_scripts/nio_csc.py.rst diff --git a/doc/tutorials/images_scripts/NiO_local_lattice_GF.py b/doc/tutorials/images_scripts/NiO_local_lattice_GF.py index 02dbe051..21cf2eb6 100644 --- a/doc/tutorials/images_scripts/NiO_local_lattice_GF.py +++ b/doc/tutorials/images_scripts/NiO_local_lattice_GF.py @@ -11,7 +11,7 @@ from triqs_cthyb import * import warnings warnings.filterwarnings("ignore", category=FutureWarning) -filename = 'vasp' +filename = 'nio' SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False) beta = 5.0 diff --git a/doc/tutorials/images_scripts/maxent.py b/doc/tutorials/images_scripts/maxent.py index 8795990f..841d8067 100644 --- a/doc/tutorials/images_scripts/maxent.py +++ b/doc/tutorials/images_scripts/maxent.py @@ -2,7 +2,7 @@ from pytriqs.gf import * from pytriqs.archive import * from triqs_maxent import * -filename = 'vasp' +filename = 'nio' ar = HDFArchive(filename+'.h5','a') if 'iteration_count' in ar['DMFT_results']: diff --git a/doc/tutorials/images_scripts/nio.py b/doc/tutorials/images_scripts/nio.py index 27f47f6a..0195e797 100644 --- a/doc/tutorials/images_scripts/nio.py +++ b/doc/tutorials/images_scripts/nio.py @@ -15,7 +15,7 @@ import triqs_dft_tools.version as dft_tools_version import warnings warnings.filterwarnings("ignore", category=FutureWarning) -filename = 'vasp' +filename = 'nio' SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False) diff --git a/doc/tutorials/images_scripts/nio_Aw.png b/doc/tutorials/images_scripts/nio_Aw.png index 35b9b0c3b524177cabd36dbd207f720a7f6fdd05..45b62fd0be9fd10e1411d19123c4445eaf4ef6dd 100644 GIT binary patch literal 25166 zcmc%wbyU<}_%;d;!_d+V0z(TD0@4aYhm>?TA|YMUAl-tLfGD7-bax{mQUcQ5-CggV z@9#Y8UFW>#dDr>xu+}UV<2auk_rCAzy6$U7s;S7~VpCy5AP`&yd1-YB1SJ8yI55${ zJG(Wj?cf*8RZ>9{6a4sMT1J4sv7F@fTp?c zB|HD4+8VBj{C`s^KP9AJs*y;%OetluWVM{PX5m{u@9s@!vSdk5|C8d>Qfybuv_7+* zDXO^kJllv_Au$92C%}Yp<(!?P!=OQkCqjf!2Ji;7hT#9D-$PeEs6l~awqF|?GY<}& zG})6|Y8b(YjHMEpK7RZtuc#~$F7z=E6AIvKgR;t~=FUb?QLF)8F!RVDxY$qpeRCT_ZF3?xn{GN^GmpY=>i zP3?}M?X*fYk~K7Zcz+&n-wl@G5gXf~eAGt-LBdll zNh88Mq`_D03X0ZJI8G=d4$-i&~XIi($3MMOU zk~*$Sj2cpUdgP~HoU{==psL4-cJ40FTS{mx3@@ga&npyOT&1fepOUh82WTUUS6kbO9*0t(^dsBhHoTZgM)G? z@FImn;;VX@!)!fm_q^ZvEG?ZaIhUclghWusH3=!{7xf=PP78?mWrJ#av$a3z(sFXd zG&KD)b*_^Qo+Z`+SlmQv(43LYp5MP2#Kl{>4{vUKzo}+BoPVU`d*yU<;asR)u&XN} zff5-Rx#1?#Fb>1On7%TOyweGmkTNnG%#Ovqo9q@2GcH*_?54{vE*@B8qzky=@j2U}=jTtJnKAtR`!{;IRW_j-6k6)Cs%Y%Bh-U1)$&5!Q_UHNe zjFhu8PhToGhUwwbFEH`n)puvpj``odg*EKeM+XN(?~b$Y4<;<)ox5o~FSd2Sw!Eur zS|4nO9%8w!4=~fyw>gk~mki)4e3>hY;?#yl{#pG;+H)TPNJ-0`$Ybq7I3HgU7BNey zRUcWQ_G3D~llQw*;K&u%E@3}FN6nRuU5XW)GhOLG@`=tc+{- zukcGwPH;@jgX%ZauU_6i-hKgU7zmQ>2Ow>!$UI$#AHT6y0weDWatEIp14t$KkfTKHF8w z6dL~GnO9XsS~2;j)U3n4n~d97352ky=N{9eM-j=7b*8YrvtSzoip=39RKEMN6K!EFvoSb$gwRL&GOZ& zqgGV>gQaNh`kOx!6S|Nvd^#0L$>0umkTK|bdwc(6iIN90%EguXU3)ZNZ9Oips6g>r z3?-)!bja-UIY38QZzsZp;cdL$zWH`MQsHEF`JtOGAXCT0FIg%a@9pZH>a#GuD+EB@&)4v#^mS~=RcA^2h!t!v-Z z)1H_?X=rG6=3BB03ppThbmFv6XktTgUcf0h9+AnX?WofIA?q=<`kSY~p`GD+^eASU5c)qFlVJGp;o@dMQ*HmM35JZ~V+Uk0G3a+kvaGV}+=IVB< zrzaXbxv$Un*Q=);I!vFyqzx$O;m{zQXVtq;`2OuKXevkV)6EU{Ls6daztcMLJ?00| zRochJRHWU^i%<`Rx92 zka6g?atB<#Xz%EF<4hEglL-Blz9vUSO&u8%llSF|1PLkj!}4n{e!->#W=J%b;_5J3 zl5a!Mo0ia~4701O(re1giKhHkf(y>4Q`6II92}OT-(R8}b_u4jadIN#es{X&IG66W zXYXV4osZAat%@#~2{DBz$IN|v@m|fouZFBtL40>PI}d>$>~*lHAj^!dR2Ge>$HAP> z{!3lmp1s+Iy2UVR3Sk!}u+2y3|3U&WVFRVt(zr7GMV%pNd;(lC0u`Y}Cgj&1#EK83 z?Gj{Q=;y(y#YhvLlIE%T3#$<`x!jWJ2#`F(VKQFdrN{lr&dk4)JU;qQKp7 zUqjqN%^bUT^tG-_-%1f-&yLsnm%;zd2bB`++i`rr4Wnt;8gBD5k)T$5-!V4$h%By! zZ{y9Hk>T{#wDJT@8o#X_N_atM7KAxK>5;s*bdP-(e^%k6cq-K}+3xw(Wok8)G5kFd z=U`^By?y>{>TTu+KMZ$QGb|kG1vxOK0j%bqH&G`07yS6UZzHlBuZiRS(YjwvIlkE9 z(JSl>!$=Z%7Gm~AMAF@b-&N1gKwQU%z}1sH1;q0`096!6qt)Xs`}ybQ$%PxEBJ2k2 zm&?uKB7V0zGUeA>?e9C(erSHzWLFOhg@}1rOd9H<56S(W$01K?jmJbGeDHQo>`K6T zYT${8%Tq_AOEPv%U*kLBU$FyLFKHgVpO#fp-nKQ($CQsWfl2GagIf`p5vD!IFUy`~ zTK9G`$3)Y1$E;97Zi*HhRg3~oDrW?4Z#s)}3N%}=BaHJ7v7u0Ucrd-Pf(DkHZT#n# zHfdSB)PbSvoUn!L?)yls@-t`tr_`d4ek3nHTxSlMj(Z&BHeX*VvW*BnVTcuM-Ed@B zRFeUZl&{G3!hlR3PD*cEk&iO>rV`}fei_j1E~~fjGY%*mxr)HWf2I3+oTP4aQSm`1 zoy5r=l-Da<5f?p$LoZnItNCob`}4S#>%E5OA1S#EYZE}VFmW8*M=pU43CiV5eVsvg zk^YNUhAfbIO}~-#Q?x9oi?Od){BjgdL1p!rX9Z<_Lus(&>VCgiv5YXr@b>AmyNl(~ zu`yvoizgHj(b30?#ZXk_Y@>V4mrSc>TDB7NUU+9_P_e3dU9Y*#LtXO9H4sGJ*gX-i z^$;IpXhyGdSX0KZ$-X8>(z+cTzqu`g8C&el)B)62>anXIZkOq72FiSy)1qWfZZ4&m zuMqNF%gD%7AA@W(^L9N;T1JM)_rhtf>D<^K)IQgK{&BA@okJl>1WYO=xE}gPo9|`Y z%NR_$;pLgp3oew+zF;TIppAD=J41T8ijL6f8Y2zr`=#*z_WcTZ;)N!i*gCqDof=J+ zE1SDm)=g`%+!gyFC8hp!ySloX*LvV#dwcuE$%yLZPQ@T54$gA7_#F>R5>|rTTJzo2 zcH_-apD#d=01_YVHSE`)J~N@$Psa@7n9sVz3yevNG4AkW`K2UT@K1eATgJwfjct+5 zW8Zl+2n|X=lapINxsp-bI2$Ti5DxzZrOe~vaDDtNWhL8s-UP2+EtN~J)D$H>J^k51 zD{AA#QZ$$NQBlpDWRdnUsEZSKoXz0rU^&WM*Hv6xT_#zbi6WXCvL5&fKr!oMIwq?_U!=H98XV`vKrK)Kq`=I~EF?6P$A01o$c4hr zOE6qq+)OdQk2u_oLmsiRxn;3sxSAno3>4*$4&!6e9yF};nez%ctp(2fvsBtYW*zLY zkV(pY;IZy0M`cN_ZB)=JCAx=dMfK$$gnewG1`#Z2`J&Zi(Bqw6#Va1ixo3MV_x>QT zm;P317bJt4WE^mNYP^n2bo1G7*f6ddTLLQ)M@L8YNtNa0%p_$<$Vy8 z7=e=#BrAF_UM}+o>Md>mD_SR}(RV`<6s{{4GU zW-2`Qji-^&`u2L)pgTyi6|3i8kwG0egnCWh{GiNP1I+m+TReNoy!#C+0|UeN@835L zi7faA-{6h)ID`?JF-H8I9*m>6>ZLY1YEpI*MAZakftJ7`0dN;Pd6i&wH&$Do;1riby+rKGiJT%T#{~<3g z4@Q_do%vlHrS)ycfNQO68RF;k&mVUCS<5E6oDHQ9_24h*%W?Wx$>h&f8iz)38JFN< ztBc`u9+OohL`3cR+5U>Qw(R4^KAfAIn>S13EuV{u`UO&NK-F9AqIR?d*mh-Qg@BB# ztMTnI^s2)N@e85ELJXkaa`W}xHj?Yqs3w)w)umTfKII_8K^}^mgAPN0XYt0VV`XEy ziYt4%WQXz0zdJYH2|+A!c3^@cSXkKf!W`-4YrJt|G|TOsI6j%&JKdqH({S$^FibWk zc0?kYG~MEcVjIM3HCe}-=EaK_+p7#F<(3L9&hJA)V8(u@TH)(e`kY{eEA<3H;GLYf z6}Q}Zr=IUVhW>{BC@$_YFC(C%n_u?w^?e31mq$rDKu#d;u=ktJX^}|JyXg~>Iqak@ zlabnu7gt@YG8~)J z*)+aAaa~iI^S^!rK;C?Qyvx{6t*1|)P8b82LOa{i+9IFW85yR6ve^Mep&$n}-Wu(1 z-}dwqcimx|Q6^6Tm1orR=)(9 zX0~0oTs(%i0bCJXd+LDTsqHnLF@q>t-asl(!e<9)4o;A^^aINXWGJp@2kh;?sg*|P z!9}v-z~MhSKK6MyI5=p1b+V~43*-#Ehe8N*^BlqIM-;F}dj+Ehnr>7^*>uyX_}(;Kr(9LJ%ATpU3GSdVaErmMe`XI%O?Ev>DSl9D>s`ZL(LxL#cz ztwK1A8`D-+EDbCw!fCxDkSpuF+%-Hg(St;tL8usBXOq^qmv)#p5uA6!@z>aJeg@J& zSq;Y1-P7t1CLo{F*oohv7;$k_$mEP!(m$Pu4p5<)rD9``uyqT{#u7r&t91etkiSr` zA`CDtz}qLP-%tZN0761eo~!LisRE>q8IS2V05lU(Qud;e8@7UzG?mL7xG*Cb84wfq zaCTq+-9Ln{{Ih%mRgF2Jmk$HkVWv3ng*9+=IL&1f2j_Ljb~F&V*EA|oV}pEkxjFQs1y0*7DT;$R~5rKyl#JELizL7YQE z80VrUjps`wQR}wgmM!}o&Orbs^qm5LO{^&D)nYV8Mp=<6V`qyj z+h=HqhDmd~;W?cLuzH3gW91b}>{*b{FBYDL zvQ^fVv|sH?t?=RG`3M@LEw&R~P$JQZjB}E`%p>hkzC4+{$>Lfftnmj=*3q4Rf>DRBU`iLrx&-FMx za1S$uL}`QY01no7{}x7LVu4U4cCmm;>v!_$s(r%0#{JsSsBU}6_f-1d_wAXz zHNXQL?8#Ap?K=Z#dLk0A(nAQKEKrM_uO_1=U`J{CIzTR^;Qheg>v+f5E_+f@F8G^& zZ1?C9ikwv501Jn1T#(l+^rbm`>vOTiJ1<%!Y5yI&q6)N#0i-l(Y)pGT)|)dvFd9$iIM zdNSQEO+UKMRk#QJTh}e|I_Lq}-OakV#Z#(uK}Xtb|0_07C>rnYt_?jVt(aI@k2iAT z&)zlxRcl@RUU=UBnzO-UcWFR$f8}bsgpTMCA2}yU5ZIm^2^3NfM5^!-)X>I6ciXYY z*tj@M(1v$);osVWMYvf<*y-C_poG)` z+Fl0$Zf}1-Gb7`>#>QuKyK~o=5ew$-|-}pRm0Exm3-(>W4t_r4fJ9E6*6MThW) zb?(UbVUOZBr{AdiiUZyB;be@uO_^vc1QQZS0;gtOm|K*KOUiD5n9OG&nvk)XwZ$|v$j=NJ~-(vpApor1CUf;6^$M zdamgs(45**{YFyMGPex=oHAB8tPw3oyMU7476`VaA5v4VeMF@`0ZCt7L*r#Z)3znH zE=q4--*PW&b{-HQk(!5D2l}gl%*lOs$EH&J7PA_wp(V3@ib47(Esub)y?>KCbUGzOCV|7#bF zIcMx?#@rWy)f_!{O54KONzQ# ziBRI4ds`j;@`ol24ULe5Wc1!#qLny7;^cIDyuxi$J%!Ui7HHj|YXAp=l5&QCa3g@2 zrLTt>XuUV+W#eggJ7%$Ol;8Mp?xuKXewW8|FHxB79pI^E6tQ?zj42=S_x-p*@`eD8 zG@7mQpzFIHDB*+Wpe6A?aGDVS>Fz@yI)m=R?m<1j-9G|q>ORn_Xljy?2>eo8VlHKl zB{KbMnTb`DGofs_JhWgmSx72vV490Vl&ul_UG3Kf(*t5BuEMP?TL8ol4-cQmb$54b z=;_&MlUQa^LCRH z6PMmU$GMRL5ebPji*ztD(I1!1m<(1aU=daD2>4e${&#C7-}~zq^MY|hU5}qbq#-yj z7UZA

_Hna)(dK1pqz;vG=PQG z+0@k3p%D+vKX`>Uko63YC_`c*aT%uwWNNuGSu}{P`mKN%G`@doT!Y~7(}m{{8U6RD z;Y~M7E35a}*?<-J`-=fhZNJLDJp}9FGp__x47>?r9V2&A*R?;3twE?Bq4YvRX%P_- zwA^hAUvp%t5i`n*oEszw$4jf)EQ+Ee)~K7tgYD&lN{SMl2$(<90tZ1L(7OR*tx3qr zQU=H(L(uVCZEf0WUn**;*P(f0Qj%VS2PdE%F)=Z;<8m@GC?_IZ6dxzd;ix1D&WDSE zGkB{`yJqWXd3+Ij{NyQ#UDNqEFEAro&tZ9i96%@f_>t;gp`MhC3>qaRrTuhuEC56Z zTFw1k>=$hdh{QY;oiB7nnj8sSC@7EKu)A-0Govq=kf%fhVfv~fu>UU-{PC2I@!Ma= z)I>0**l)7>e_XJOMI}z~qQ_F&MVmDcvLOv??MHgd^+x}r?h()}qrIbIzT;66fAH>K z9Z%&f+fSQcmSyIcf1fbuHyf8BGD#B>9Jn|Nv|Qfe--5hxMe<*crUkiBB6g`+@qs&c z=((=-jfJl-#n!-ZrJ~&Y-*4(nKLRAwz|v_~kT1cb&|ExNjR|ZTp)Cvm-F>HgP5dj@ z8g1xHAJk^|_18Jw`Tgg?RYn1_yT*hG&FsS%I@}y>Ly{n+tMLu7pIGFtVxJt>W677Y zr0Gg=EM!u@DAN9dpLqXas_}i~7+6do``*+4&o`5Lr=k)4TdpON*O6U43n#*`Ij+Ka z!5X~rzdj*PT3^;7X}10_;h3o1C$j?WZ@>@}%Bo?l^!Y z)*8Q%3~BBb*5NCAAi{mtQU3KX%2+qJ6*<14zqWhWE)C1=>e!gw-Qnl9qnoxoU(AEN zm?IpbgHFtX+$trk@r8m~b&x8c{)Y#J2cMl7+=|cnG@u*H%N;&hu2pGjgarG#a7Dp6 zvEU7OGywHd0H>pKDx#+}pqwx!8}Y-fb9LdXoU6Kpx6#2%?MJo*v!8GZd=B%REU^{WL zp@2uJ-4yw_`1FR3_py`7N`+h~hm!Dy5=A$HAB{38TjBr93jl3dEOTfCLgLO;Rdh}c zgT8@5WNhqeFb)?I@HH9F-{Zodcd^GyI@d;)y6s4AUSb zr6{-Bu9V1CX&y5j`Y!;F4bY%7!9neiDblX6>?K0-bWq!TG6kt(ktnEU&O3T%#~!Fz zwSbABSDFPuJ{RT2PMztSnvz8{U4#c4aB-9fXjBqQV+@7EP^87^{uL295{BvIWqQhjp=|6sm`jyogz2`X}9XHx2B2e zZhJp`_%1dsgCR#;j3XzRllLQUuu9C?ERrGu-g=N02z13V|66Ysr8ER!WXNhh`275# zT-?LMQAQY)x>^|S02^xIOkSqbkQKuFj0&`VS=kuag68rMvGFaxzp*vqD1=g9Va>lR zphq{ygRTIB|L*J4Yw>>$RJxl@)d2IQ&PXPf`f*!HxKk&>7Sa6JeE^CgADzI55lt0+ zWk|l%@6h<~4knG2g9CaN2G5M;DQIa`+}xfH4-W&)^~;#8g1r2}i3Mz3XzyR(XjPh! z7GpOy3u^-$jw|C=d~7RdY;B~QgD1UU@r(t2KIg*?3L{ui9%MJH!;-hmc8pymTde=? z_WJoLMG(ZP1%0J~chZlao#jIh0fuAGtkY!Z>UR3ZP;!@@?b%HQTCf;cnP-`D%SGVMC1r9((McAx*uXJvJj$9DA5iT86oJv|WUeT!2_n*@j?&|$~K#bpM) z-j*2b2mR2AVFmhuTk7Ms|ebH@KDjA9jlWKC(YQ_$QqYX{Y%0>zBtkKD9v1 zSb^|?U{^75Eqf~=rSan0J@UY&9n z)qk++OL_n4(*W>`$SWwA0S(Gyd=Q)jG#s)7@HSu_wSXR{J$uWW?RS>W>d!9`IQ%?t zPZ2j_jQvo#i~~N&7%uFtEU$*?947SH@v(GdP#P!EtLX4(EY7SdW?Kv^gZhEDEXadG z5znsyo5Y?3U477Owj3`uLRw?e(|a?8T{D}*do9i9oBd`y=l!0?t*@<>gNC+$Nh(t$ zX!kh*36?hC#<{^RMwU@FclgZJ!XkI%n)rgw(+&}*#XmASClR625fU<|fl0%|>%ye< zMvH`>j7Qj`e3m~74eIicNq)KQ=EI$#Zd4#e#yLIPssKB)qsLg5m9<50Rd0 zy;3C6C!br8lZi8TTs2l1&33Al&3~Dz9i$cPkraxF!MQY>BTlYeY90jmS!<1#65{SC zwWXGsB2H+UpI$*ddoK#ukB|}_X#4_;0pVkvV#||_VUVF3#;3Fjw9A0+VR|q&Ha0IW z&&zTiy1(pz+cTcOvKFmpGt6Zq$YuKw(!CWYYzX9N%wZcF0S_-R;=Z-^EpD6_wWa@^ za~Gt5yCMco*z2bsjsMANX=!m6d7q2-1)g8Qar?XQ2CdB}ITFC~BKJ83Mc6i#+eAe{Vg7kANwvf5@9mP7 zxqT798nm^^E-x>^{;sd{`CYl?pR*{a#|wpT z18TzIq;>Y;T81O9c_*fTV{5Rv@*YUY^2GA{zMNQGtffY-(v0^8&A#$!61pFEOd^RmdfAf!r<)?_u>^Y#i>ePmfu{ zX9-qV$x{^3>Ybhl^f(iIp@Y&b4vQpw|Nh8eI}scqfj3i~KppM@l5nYRquQ6pLhEX4 zE_ors+59?_Ln;8Y*^(+zuNtB(PCg1cu{e!GW0GD+BTr!psu~GheH2T+qQqhGPJWg2 zf+XH^c@{XZ1m}JE0AI16seKQ`zT?wz6Q3)dU{{iizz|VxxY1#}eCjeCHgT9TWI=wH zP^9oIACIbWSnZ`v9%kI*VTsXT%$s02j)cE|)sd1t(oCeUuU`xLl{^{!H@Eu%>!2o{ zcY;#-`2|9qQ9Y}9tzvl$vJ&Gw=!fRQW9FqMq`bxg)|=7K%e#c%c(ccDmj!;S2Z!5A zz{IRt$$Eg>Zs2k~d+>`?_ylM>uCA_vu?`da5Mz43oM93ib^8Z>&lG9YJpg%PXIMBr zs8ukk6qN;Q4dOeq@CM&k02^v^egKmeB@Qy7C=K+zSPbVdsE$GE?2n6Up|P?f_cG8R zYTZII#SBFhz*H98`0}Y^H{Oj}#=8~WE!*>2@S$9s{i*R7^!dtI@-?HOkefW^_4h{y zsHkKKAcGB$j!J`bL2l$m7y!D&!mevEtY<6)C!c>-VONaBXcsZlkJhoJMx8nje_s?J zPijz#%!0NCv0sx#v{Abblj2CzFj*T2OCl+kX)6?D$qqowL9cE>r{#(isWAXEkZzGS z=tJ1p60uaZ@`0wqN~vX~pC8qJ7b{xUUcupBLkS7XFBQqqDA0id{W6~{L9dR`F%#^P zU_g~_0=vsgB_$;Q7GJ*AGP4uL|bas zDP3w@JmtZ0hx=l6HH(Rs|JZh;&RUM=$w;Kmf84Nj*`^V)p5;l4G#P#U`n3@d=sB;D zg@@8*2=$JvzSb%CLaL-+ilEuj)fsjg5JYU4=!*prjq*GZ$=8K{QeJj{i9pKh^u~k) zXtA;#{NsilKKo7U{oR-Ygo;!~6$mjAU#dVl?s-5Z3pp+zy$F5UI^1tC(y{b^uDR1O!P$9;wafM1G0fxWan$BC)gj9Vx_6 zthdgaP$!Kicdm(;qJ$C3_F3`r-t0J$$Q7l&TekO}8Nf`Krbx1c_A3^I zK*Zgw&U`ArI?yOUHQzaYMt;9B@VhV>W*0nZQk9+Ee6C?Eek1ZGMH}*Gahk{7QG%Bs zbD*7w3~=%PCSn*FaFaiBmBvT!&^oprC;#Au0)N`^-) zcUm^Te`b92MZ|Do_;l)a^Brp2=MmH0g!!Jw2J9IA zQr6oD$O-5R%sn|+$A+oEVbTN#nf5|hJr-8iOGPD?U+TDgU{qBB!98HgiF*Ew+>UH? z7r0qPs6oX-5rEM7O1EOm17#i;Cm9$+*CfXYW+S(zDFgUic#>Glm5a5D)LicESX&R%qNCGh?!%jlb? z2p$!{0VZhSd0QkxN$hHaC8TX*e3 zU5Qt(ScqAaOo2@Rq&~5$4aWbTaOG-NmBHCqrgP`yJR&pDk9|mno`r6gALzvtu@{AA z9vuXtXo33wUV|Kc4wyJ7OAn~m>bMG__ome!+1{|DnYyjKZm_#6Yl)d6Bkp~P?lKsE zzIx_uhWU4jgJ2LdBFP;N2jk8rc1qqD8ynN*AX{4YX(j{)3=O@_%ohuQZvbMRUgx?# z_1y6(6%SQI7z)3DKz)l@kPe#-`xmLL*3bL|bbk#VGAS(m8nh%IFXEHU8M`>=d!HG~ z+>l7699x-YfQgIkejhgQFDUT8ZD`elV_T~xBr)~=_~G(@HgeAGeXOz@2VGhlpRair z+#Bp08SPia7vW9!_t_+wi(XoMUqo=nRgdrQyPqgyY>$?hAygV@oE~;giVzThR`eif z2mfCbO6%iM4b?;7>^%fV3f{Lzmf_plMfx&QQtx)P)YX^XH{m=}HOu9I$Ia+~GIm=@ z!f>K8s_2zo_;0XIq9KVSX(qj<@onIb#^TZcV0kwzGY}Y-fj6dSaB$!GnYeiN=19J1 z;~r4aSHL6IZqMggMMWD1t%0MNn2>OInHf__<)0|tXYx!3^YJoUP4-%M?|FWPW`>G8 zDTukqO5Mh37fN$_gBHgB#E=uu|1)V4-caO!e?5N;oF7kIR^&5IiO9%gb#>F`{qfsc zg|!171eU$D@X4C4T^54Wt#>Dh%mn;s5YYcgX~aFKfM4<3b(mTiHKrtzgF^yB3kM|RxHvfk$blW)_gVrk$P{*NE19GnfJiG9ID6rZTch_~0##UGF zdfvS@EU~+7Ki~h{75oyL!7vf>nsy@u@e%_hcBnc60**l&uo&-rnyz!Tj?a#cifU}3 z4;o~1VDGZzHyfhz?w2shHi5K~CK08lk!;5Ap00%$dY)vif~yvyD6bu;UO(teR&=X_ z(rQ#Z6m@2CU?6x97zmCC47+>Vl6(xflD=yh-D~217)=V;cxp>J<2vaT@TP<{5oR?H zIf%5vewA$aZgf=bY3VgPbV-Ru-l<0SQ>i3O44ErUok{M!t~zod9JQt}i)AQN*O)Fk zrq?{CH=zJJ2NDoSwC3o-XP4<3;&m#vb@?TzNIeH} z6DBYw3gc6HoGFIs^>86iGcc%7oxkA0q}fYa^FeKbaabF`zu`R&@Y=_KkQseg8~+yB z^e7cAkKY3uLRMfO^wj-tpx6=-@q|&BKj6`PT6vfsZy9Yl zwz?@z1%1@fV;g#edQ_u84jG!F?`HT|^et$1I@wrvs6Uh$7!%Gx%Nm;fdnMcy<`x?i3Ap2AF@p>jwC-wCTsT)cm84rewFmiOl-}2|=^_!8-a+t*m#}gN0_6 zw>Gl~^`)V1fBvk`H^c%Ixi#TL2uHg}49&vEoyJx-BNR=Vus9FQxDd0WKZ!%8oHVbo zQde}9T`MxAzJ7SXq50WE{aFVI3L$_v*T^t>c&!qE6IjQe9xaEzT&;H&u7C!szK~;a z^N7IPjCmL0arc(Ls@?|Itpvr{oGSljNQFG*2i6}{4s#OBD`b)XT&mA}q`HF!m*fu< zpMOfK{MpLLeSOs7ay!DNTFR`=CCN_8x?+3npY8eHt}QNRPiV+*Z=fN2A)O(C<{fev zHt!Wtyp~wS=~@0L%(sXqFg==K4iLBfw(Im=8L=XGv?d^E$7@)Gf-~OD6KF|KFt`K5pk;W$w_>p-k_l)Uja-ZV?8~G&OV>0 zQIK0uDl!-`#W;mJ9PiI)sNT)oZ!N9h+%z(>KDyDo6EUZPQQC(mQ^*q>m7(B7RPe*G zs}c2}_%s4s95gSiJJq@yFe?|X#UL!<(V(=sV@o|m;$B}OJFw%d>uhGD=~81_V^EJj z`k%iE`LKh>0DwC6y)gM}UjD>l$BD^^5Rzoha_MPtYgXl^2ha7!j#z3}AD((6Vk#A7AKbQ&rYv1bC+BBZYham)o5g<0 z-JsTU5Jcy>t-SHvrz{?Tmc2oobP!#j10#RS!5f9=f^RN^!vy#=G;%Nc{^AI0*3bp^ zmL*|%oNCmR2%o54ZihVfX=)2c@8cn;(4wuG-f`fQ(XkGCgai(Zm3MQ%@U*;~t2J|Y zC5Sw`UuBMe7mZ+!V;r7tmZ0>+4W1q2UANf=b?o=ft2TZc-&D)Gv9>Dj0X`3h7oMJ+ z2~UOx00Kcu1;3|Fl1sPfFv@h9U&Te6jr8k*7zbRnS(Oc4bW!i4Y>0K4v^hq$2V`RfV76Z;QhV`?|QE2MZw68?}O>N04hcy5X zk=w_cs^XqgJ3uP19-;80EA33z#Lvvk(0025_j8j^P4_2YU$L>b-#H}mA}i1`@5^V^ zOE?z>8XP6zn}w`tYbW+2i84BFnB+V#R3ZFs=Z$;`|9_e@bi z;c%h}O-*~gYtUq_FLxTxqFVk>2E|}YfNNw&I_-7A*~rxJ?KNL#^S5N|n2EM3szNsn}ADIO(6r^NUm?o8@lF|dJ*l(#- zFSWGB`q(KkK`V*Dc*zV~MbLSEzGXRM$s1BWlB0J%(i;zRVnBiHDe~hHT^P4gR-e+4(9@J8Y$`1gK6cDyiZP+0X=dY zGLBLj;OAW5DkwQy9)QKYe$AJ3e&}p`)kshw!@!~2*L2KUy3wvVs`4OML@yvQ372*F zXVssFJQee&0x`;JhcucriQT^iJ|JH3SN`iMrmK?m-^^XPq%d=K@skFjv-yl&S!_-x zR{GMvCjf$N6=bpZ?#ZP{elJ(nK2eWH zS&AcNBsh^mfZJuilONZurWgQCl?B`^1GgMC3=RJrt@felHhS>@NeJBVBuhZu2e(Pf z{YUM zU}opTcK5|Naz&yKyYGU5%z9Gp=km-P{yLbRU5o6)hKfX4+N$8}-zGI9c*)4tV$`Wq z8$-S$4OQWnYP1LJoThn`y8BmmF4()B(SjDV{=YY!9bU1&cL0bXo7$ZIsF1 z(4hu;ID4|aRhIAr_fUX5jOqVRUe*d3Gmv-X6El*Vq4{O{%y-v*%+#I_m*v0SbO}7u z?0%jayB=o4_L{#KqeYef)6(nNmrvPoUDemzvYu46oQ$99I`!&dEjw0 zhgC-D`&jCR-Z0SW9>O!!REa{7%ZOolQCBy-!n)q7%t24tV`t2K4!hGsgx4)Q8DJF)z zQ3m<(47{Cwi;z-9?TkP}Vm{KlyB8OC`=8y4V`;H80`StK;WMt3Ud`sZ-!rQVQ!If< z^V;>6B?_m{2#4uQMSdCpA*levmeMN_Vf0aXE%%uNivGYNoHQ`4)YSh|QrDn$43jy8 zOwRL6+gAbf5*ELwAI~|K@fXs#vnjP`cxd1*ag}$pTcVWdpMDM_i>7WWB_U4{``e_< zV75zm6p;SEEYKxX+9wh8v`ul%*0ctnf zv7=Loar6gzHV1n3_GabaPJsjot{=2DFht`t5BD7(PdMjcC&^R`xaspB^AXivEnKCFz}Zdui1OE!d?plj{}d5?2jkT$-H;XvGWKqx_R&VOI z+tq)-#ZXoJIU33U`gj%to)Y0v5t=18c(A|!`*OpH^vzna!duXU(G}thEy!oW7iu1U zx5jsL6H{WGpkW(SREFjDUDNWH68&w|`t$Qbh-4W;UkUvQB6`KVKWeG@?*{xeWOtXC z#T^S6HOxW+A2<&@rJQ_Q85JL|eP|=W-``a%SpLBby*M+oQU)}-P{40ZPz2g*RL z$4ph+J~1u6Gx3)#j7KwAl;ZB2=Lt-a_^N(A+0B>4%*ss-+BNQ#hcds&eP14;4bpn0 zi1O|))E^h)%oP8I#xuYy6NVkC5B&5{ucNmwztzT-MA`mu*h^m|rY&T){^R8TfwnvI zcqt{Jws(ZB1o?US+V|pD)&*dtql9!H!tegz7k7Z33v!0`20_Gx?DHAl1nRMU=3Twje-T<(B--gv)Digf36 zaQ$HS!q{T12kNq9kqY2Z(MNILI}DHHT__xw^6k2f`*uD|#>CRTNDzrA|Hdo!_&TO0 zu2}{c(BpkCF{iO>4oAPL5>gXwpH#r!5!1s_VnC#gRy%*L{N$)IR@RdDNoo~0jq|xd z4e0wQQ47axT9~Jay%GY4WyY3siG-8irER{(QSwsIs5?CVolT$VZ(NV*fDinshf8hb zC>QQ)CBr|%@5GzHMbh6}4Td89ocUOBPU?sO7qceU9h#5*{V_2Vy@iAowQikX9~M3T zy>~QQI8_zZVhk?xkK!fL$qv4s9?qU&YoN?>GVoVF!U~wU#ju^W8~)+WkyZJDnK+ez zh~%%u?8HVHm}g)A`uv}NxdxN90gT{2Wk-7ex8?25sbv>#Y^8mWU)Y_=A-rwWY*kl_ zIQ%!Yu;y?N{%Y()F2}sy2~i(b7Nu^qAHzz@u`G`{KeDmC|5PwiDi=7k>U@ZxS4&m*VEW>*YwwA#lAR_1|4@8{BcRB3%SsOdv(Z$Nk1(~hCD9nM9Ht$HAKmO0FdtF}3%OoAP?^t-yyk_Is&XyOL*;Gt6JFc8nH_yC z!hRn7N`p^)^aj%Um65B8GZH^`9Sq~sy#JL!wQo%*r2^khz|weVd#$0+kz*a%dUyP! z-^!PdQZqH%V%XH!D7wCP=6jjmtZh{5){%KKd52x3jt32c0djXgd}$|cA>hVKh40ak z>)QP_NEZa$re%siaC&Jtt>~=$IFz3=^Lt`wPHi=5!Mn9KLT}D>Rv9OMb0!dyLgH+} zJY{AZm15ZP0$C})mU>#O0_qa|aR@#!2|mqqLGAo1>em;whE02W^b&!W86Qu(ziEWQ z5fvz%yZk^_Wjwp{(E9(=*;z(K6~23WKtvEF1Obs0K~YI*q!E>pQo2OCM7mKz1Obt5 zL_$D@?o=A2yTJjZYluPWz2|?eBdy-rWSlGWKH4p>h@7K`AY-L_I4hCfx!PYh7w%<0k6Cqy1X{H?x|Q)7)k4 z78Q5IVop&n#P=SFTk|SZs(OdUExonNEf3@0{D4;0da+!ZdsZJ6ggBk99bk~lmRuOtd zc&|Hc-H7}%I^0xdIJDZDhELgI4+=dvo_@QL}tdi?hli#UvPsNDT z%II!sotj75f_!1h@GE9+a#mfY5Lw-7H5?5OXf597M_WR=9&yWa<|(y%&)Otsn1Tm- zWH=uqoYc!-3i;RzP>sw4AL6v(UZbJtEICI%5U{I>hg%e_z0}7);)?=rbWEH6wkrMX zM($ZoFLuoQU@Xidh74N{t}}h>k6XoYUAi1k4{>eJeBuww98_$&BPVoYVP`OiHt!cHrj+5FFF0! z!Y|E9C87*T^4IrilVZ5K=XaY6WzO7>477IoFNcjXm>nG~zHT!^qawsjb@v3{^Zgvk zP+a`b#O*01Bz|J{wN^m=ynGVKa=pLM5&)Qyr`mYL78SA+g<~NDmzb~6Ro&8b5IjXL^O!%=s;@fe4 z9xRqFaH7wpW-8oj?J-yD#ZLJQQqQXYrV|KV)wut2YUsS7t|mMzy_QPwELoACbHOCC z@UtSW%v1G7#_QL}?%n=%aUnQJpKJw3XmRnf21ViUH`V)Tysy{>E~xJ4U&*igu%OwZ zEvMrbAQ5$0>T|4mpKc`L8O^)T0=vaE=qYj5@xq0G*My0PhF&2D=fik2(-gPj_}QLE zy$>ITxruAGZ;(|~H0aXH6|)qX3{X;|&GM{qmbf0z&aiSk*i5j~(G|HSfxIMl zBV31_)8D~Owovd;vB`?Q^q?@EUNkgurR9pid)d?sC7~vB#b-{0y(Xe<^IIgSyrbgP znHLI`CcpS5^<73?eTv?{z3(qg80H+_WNs}pYoJF^>UChBW7uGr^4vWo-p{lcRYSMV zIrH4egsFOfaEFMA<91BV`4J(TS+C)SJ$>z1&V2dyr0dq~EUM%8osNQbv!$^32MlE= zSt}!6%wG@#!%;l4_8QSrCfre}fkfG#4DIKm`du1Mcm*e186!HZlb-EW)1`0nrpTDx zeD!E@lE{W#j;;KB*7+rarlOmk5YrOr$n63n^-pGXuC^2-hCb;JaFX&@LvCPZ+Wb#} zb(S2t&pFkHU)(TNxh{CiQszTB9UX0VEfYFUAR{1*Xy+26&9#LyGBG_Il>ScVVMdYe zmZE!jOUzIueE%0q_8?d8YxVr_ zOp5%?gdN|sARNQcOKUc*%XVw|q?A^?U)-SeiB4f7dxlbP8Tw30e-AXQE5bUleL9#h`ATop|nKjs#&nur&v|3Um5 zuXvWXKJl<&v+NFK;6?u(aq{yI6>n!=p^Lg-C$$^;ba?C2>|8sPnMOtP1obfCELfhP zpTAjlQFiz#s$4pB@P=iH5!LDvQ{PJAykf;b5Q;OqxkHpcqSA#_t9g1ZoGhK$x}kvU zi5uth*V8*`8 za-2jI+oNQ?kGc$0;fV-9i$ANK)~H3@9rKRFA@8a0ePja5aG;DZ1`rv`ZjY9KzW%ti zu_gwCFJ*{fQH$Gg?Vp=zr5KtN{6?^8QhP{JUv5$-x9GIE&N=@4Al{k>{V=EW=IHXN zg7+imW8>d+JXBVOA{^aQuLGN@r(a0ZnFf1qP{jT%i*K0IXq0(_>yIELEgur=(3t#Ti$RPn3e~-re9!DEikEHiH==rZ=zA;09HRA#m5D)Ke_@ zwyJW4O&e-b4b<1O`(|UAC0q4k9G(*|%!ZllY!Rv`JYVf5V2>Xy-?_1aK z{`CB6%ZE@x##=p=*L}8rw#KCeowXV^=ebF}B6@H}sZYCeD}3r7!5VQs^RWVdDy4U4 zc0lSia@wX!Zxr>U-2>?#(%a~rn%&C*qvmosBOer8BFLC!xj+Nq4UU1z0sN~p%|khA z&lDAHVtgRWsIIQQ9ouly6776`dPq)B-vjxkg33Aoo&g_%;Fy1)^KRr(b)2+|vQ5EK z`suCV1It{xVfD$poZNyxo*QmG!zn1w&>i=z)hpTeo~zi{9ql~1@ZbXV4bd+}UTt5@ z<)h5RN(bxm&Xs-~T@BlN)MP%LT4A~&jdp7IeduLECq7f0TKh1QmJst?)gy4*+nm5f z_7mgwN#y#8*&t&xqU58|95JgYlcP29 zo2WqT8j73LwvK%S#WXHN2P8NRhNxljvC7Pb@bSub-@*+BJtC{->&7d|va{Zq^03m( zO6xmuZ(aBPb@0b}mE*>;<%s8v*EbsJVgzjo!-Iwk2y%&7qvx^j;_uZ zMDKO&-lvdzhg|y}m>I>xs>c?kDodXnc$oIBFB(^)KFbZcX#Uhe-Yy}Wr+oU(8L>N_ z(|X)+d|QnzmBi}W+DiYoOIr=+^^+{$TSBNDFxy73;Ge)gSYRNI4R-+PrDuV)5`fO) zTUy#&@Z?^pl!Kf@yKr0U^U;%_C$3>&Q7QRm#%K5D?7@nUfN}@k8aGV z=6HcupEFrw>qqwHz_C_YWT5P8<=YkhGr#6i{BBG+a7lYyjedBQRCO#cs&(E3usIk$i-N%8ZeO8vvxc64fCS5Y(0G?>Jt3_F*lDt zA81pNFz~<;P0qyhAC!v>K2}7b8xhX<`*S6Y2o6lCHv_vtX3!g9LWfnndygeh8EbfV z15#Dv5RteH2{Vbxn93QyKubQ-vg=>UUWTQo{Ky-&DJ31+q^^4Bu^B({b>GnxvEm$C z*png}GPrUZ9{RJ3^6-EO@e5@v8ynk-@)_s|S635=e?Z)oeW~V)u~~T0(U@s?%szif zhr(x}cD{>=euN|h*=7>YZ{ad8-{_0{o_6ZP{vd=k#OpbU}3tgQwdb-s~80 zqzw``H-RS(Q6V|^4wJ@h*8WIiD^F78Yc0PM!Ts#p+dgJ_g6ePH6(||ZIci;HSuBfE zA5wksw=w5tLr~r-t-8}mg&vQSNe{d~Yv! zVn%vD44N`1#d~353qI0B5`M^1SzVL2s1Qll)5+WK48H~;uq($~jj!3s7CMF%hZPmYRB{`T&+>GNn7^vyweybh`iNr6y)U z2y;wkxG5>ot8>HWFRksyxTk+*M+Vufh(BhA3EVm}YR}jn=j+)K&=arH2|H(pNoxq| z7F@1DqTbv~s-7v-9lV&>x<(cnYx?O>aVqn@mP z8p3ibck}1VjC^P+Lh@=ab4+t69;ch%3z1bDx6hY%aI&s>J*qP_agk~G*Z8Hz&rEFj z1;2@^^$jA&6BPfn;oummm@+ZFip ziDr1WeJO(Z34hl)M1o2aW@l$p`}&kCo)Iv~v^my1B=v=uY7E$TBI4rO+as7k#^eH) znb%n1&?E4?0W4JM{+b=E8IZXy1Y`sYte#=o#EU^t1l_OZt$;u2pTH6J3+;O=gF-Hw zw15F*0J#Y3l+Np8lGZ=C3v&t!Ne5X-1ZV*L#j0XJ<%WW=F)Jt~fo8N4!IvRmvVA*^ zylJOV8PWERfT=lC80|W`F(KhkDy7SC+n0dJe>(_|l*ylzekCS!d1hqn%EC&JBa?rn zJ}>=W;crJexCF{t7E{EHkw5*V2<%03fBa~ERiryCoIyp|n#>D<9s)wbmF<#CG*60B zg`MQw-{{G%X4$BU+eH4em>wPqRSoX!lmnuFI#AkGsz(A5*&&lE=7qGR<4Lb`QWe0y zEOF@h+N1o$JW;gy?a*Z@0;Wf|aA3#%pLjw%zrHq)$i5D}v%|!~a?!w`>XxfePGb49 z(^AFNY^Uk3wi#r4*q1I`iy@UF2!DzolOmWWBV_!qPidA0+yEQ12X2U*n))dOg|V&~ z5QW$th0%-3y?x8`?%g{$6JV0;RH_P)aH#;$2qdYG!|6p~CCZBXxseeQIQ|$VB_$!m zm*Tb*Z@TaDZ30r<0Q%fvWb_400whK=NDkH5u@W;9>}vrlF;8d(k;MYP0b2_!cY&sm zm5b|hR1^jHA7&nB?jwhW)S(i>SN9d(tB@=lO@2qIaxo__54Dn=+G<&Lje?>%icL+> z>tE?;=Utp;X~-1vRU=CvUx-ymfU~9ptC5iBtXt?*Q&IVe*+IcZmVl7`VbG8QoHE#p zCO95cxeSetKD3$OrR2~Z`{S8C*l~g6kp!3^$H&Jrva``YFEMOKrk*r+vFz7_mlAwY zzaY=p-qB$MYvKwCydOd5h2%&WeW9W8&~aHA*692GoB&|v2kAX5Lk^6JZ?~!O31$a5 zNiy9Iz0(yWf-r{9u;EAI4}PG^RMykm=)_Wcz*AsRZP+x}O1UASm!y>KaydsP+YJ>e zpniUKs+THYJN5kmWk_y*zEfsqR#q;AKr`3eH7NSmAfOi2^2#En-nfm}geQIo*6?*L zU|s*k+1XicZS7Rznmz!elQR+ z#%k4j3MFKf33(kY#@cp)O~4+mLb2lqJU@T=w{~lea|7}B!e;g>M7A?*Ck9mM?q(vp zs#m_p)KdvGd3#8$UZ!N%ZicEj2KO`9(e5SQ3ulN}El{9@i=VkBy4hAOQjQ+O_XTkb zG1)NM4zQP^VtP$W_>NhCZMn>Q{=!Lg*?hXR#k6Fnh035Z6C_ba6_W*B)fSw2!>~&8 zJY+WHrdxxkQOP(073-|79(DiJG%pWKKMAl2P_-UZPV-HfmH^iBTuUo{QtPRO1slYS zca`RA2L|5y3B#1C-XGB!Mum$bX^Eed9D%G#9qeqm;GT^>06m91bG4Tli@c6M~A0JD@^T&&8L8PzneOh`y*IaO)D3X(ZXU~_Vx zr;1quKdPs@YsD!lsu$mFTIPl}jZRF=gDIW<{{4G(Z*Ol6s_r!y=6bVNCGz!ZUM?;y zds8TrU}2T%MBUXvGy*95vI-)-d;I(_p?2qh;H0YB zbjBLW+3=sk*9Yt4a%?CK%>9DebZh>^-NQo-FruhD0nfvPl*coxwpf7gHBtb{F9Wdn zIvLql+zS_m7dj$e7#V$qO6r2Tj_nUC0WAKcznHIK+n}oqZqWayb^NCR{XY#d0&(lV z-?W|HaDKXVg;tPMZ&bvo4&)Ptol$Jb)k;iA?1$;^0WefzJRG5|Da?Nh7fa8|TFO1t2s-$TxFsV3qc#w_I*d>Rc%Oe+mVJqiUhFLKIJAYT zsH4}0jnykCILaSS&Yrab;tE~U4DR<$fG))z4rFq0xL&uDLI->*FYlnUIqfow&1X4_ zGTdb3<|eC|?ty+?;O-M(P6D00)UV~${3US=4!8vIKN%K*^#kh?hM^4i5%X|vA$VU} z9R}iQ>}AFPq`q4~s@Mo%_zf;Y1Xz8rK1S@ek&mzz33vu^0cb}AG@+fdARU&-!1^2q%h|IGAywg<$#HML-37ZUf1Z*5O0mA z)ApR?OQh~In5b3>{{{2v*1D4Ug2pawO&F>H0-k`RMOp3~?ZL$ca?!tX#IvRvA*IWy zs;W9x>-NyVYyUcH%BK;IQER;ew7~$3x3*9RYK2=g; z1s@iW4xodcTH1A)cA~1Ps!8Il0+*mC2-~poK@GHp0sJ(< z2kdVIn8)0KBp+C=fU`%2Z6%S?|7Sq(mHVBX9X9+OZ8+mI#EZ5V%z~OdX?>1}6STk8 zaJpePJ%F{3K74L}+Tf`NcaLClkCL99iV9i7 z@$}YhSAv2s5C3x^9r54){%^VXU)7lZ;Z3{ec7ND;%l*^m;kO?U&mJquAf*i7{~s!_ BKdAr! literal 25036 zcmce7bySpZyDnV=GD!E3(jg^1gmj~H2ny0&(jX}!rG$i10@B?LQj*e*bc0COdFHqG zxA!^ge0#0)*I}($;5hTX?-TcZU-xxgPo##LJRS}O4iXX)o}z-R77`M2EO>EXVSsn` zYu7u#FEm$a#b;RH#~15W1o$_$lY+i05)!@{;)R?qnP&~&6n1;6=cesw>E>zXVu56D z=H_JQ=w@ea4)d^ZakX}IfIa4V%*Dq6vvPBD65;0l?+Q!Adkjc&bseTN=7P80)Kr=fKl^UBC~z`#`vE%5;Sjenq8kP>C>n0$#}D-?kp=)3QJo)=j1Tnoz=GT`d&CGDJey~e~;U?gN_okd?Q6{ zM~Kbcxc8W=e(j-&*9!V|&EB4!L){8SBDd+H?LDcfZ`>z#U6+d^vnqpP$s{%~Pn!Z# zOfYcavNE1UpF{DDpD8bw{FnUB%N>^_T5kRfnDiLA4hnR2cN371{J^C3g`1nx*Uq~@ z@Fgx_6wyCc+*OzfnG^c%8V4=2OG-XmTzEZW`*gmP3`T{99nUyAKHgPk)~Qj5Nz#qQ zMyjK#ilDzdW-cKip`fV9&Fg1h!lTQZli?ImPcBL!hqdYQmxyU{b8|atT)G`v z;T=fg-aPK1J-@mFqkK`Op{3Pvw9<*`_TlS8XsBy3l45+C0ft~|G?;l!DfuEp!u+WPvDGj?)$zfMJ%`T0MNjEtDfWoKm} zHdOElieugGTefK3ykBvtsRQ3sQ@!n$s;mjMzo|HT{R9gs?0@q_Fm1QuJ+LXLckqLiroVs)q_WQRv36~L8Thp5py^2?ne)o-sEvd7! z1_zCY5+7*9`{ilA`8Mwl~FOPq9Q&lEDX+V(Mw23 zD5IqGZg_Y0pD1hn8X>Qx_0Y`BEX!&NQz2f~)|UPG^XGwSxjV{;4>*|H)gG<2MkR~W zePZunLiZMai&8jYLmG2&aRKLkt?-hE45|Kxj+M1Q-E$I--N3;S1e<>2w@JQu`|YD? zb9wkJQwQ}bj?zyC4i1hsc40}4yL3Nl+lMIh(7GT-#io0ER-}ywl#r!qnbaAI<>z=b zA#?qb7rRUPSP}4LjHM(>Va)|;tBB8aPSnKF`@Xu1WLp;4Xvk9k{P+E|Ax7K#6IEy~ioviKw^yfKU0rD2m6w6ig)Nvn=yD zcQ|aprlm)dg*n9uR(?vD7(_C9Odt01!U;)jstUtJAB8r1YWL06c3wegsb@n>`wJLr z-|hP*PQx0flR=@AZc@`S*G)CK_hjuqQ}|ff*(rPDZ^NL__BU_)Y|L%9bxedWa8CD! zT5;2y>&7BfgF_-T;R=Ht_E_OA9^ zTte41vs$?}~Cg~vac>RQ7{FVEKCDovkZE1FTU zN*^Vm*>c2Y{;GA%ieKz>TzlSSfOFpGu$czWY_{1)D9OZ|HC4!QqXWY2Hp|*nT}>R- zEp23!<~IkzCBLA+Y#@=-s@i~r+r*BW2$x?P4uZop(Z4U~Gb))hZgLTG3(IA&;PrP@{>|xx?m~+tNQOwq$HzH2IZqZ2Gs9`e|JNw|Z!0XDu6FJiQQwFB?Ec}dq>43I8BR<^V{^4-I#{iZ|M-ldP;CX9D zw7_4)QeO_3-brg}#;Fh8hfCfB)h_v`o-YN^Hk;tL#WBcZVPmhvD>pc}$;Hs9+SwJU zPcCzHpb+4`Q4`EqlO(>ED%R5WlA<%HPCSok9&Xw3dsg!-{^7yR2aR@{geOKbzfh2S zm9X-PGNXS9aHz= z*ey1tL=rxltrfLRq~TB<4-Lc$6H6Vo!}tm-4m2a^L?WXV$~d%Xed)RGP(tZ$+`9SX z`7i0l?l7rBI(EAG7I6lN_xn97D#g8BE|{N`)@vwP$@I$2(S?PD-&2cD7O#Fx4xgA6 zJ*ae9SIY4HqbPB?M%b9DWHw*t3<^V;=b?#j`rBix$ln1y|jsr zOCI_=7D#*Fr`N}XKmJ|NR9vwdH%6YW(whl=zBs+{sk-6)AjawIw%L`X-Hc5mx6!9R z(?yd6Y$xm7jX=f%`#47OF1hmc@YpnqC-+4s?Ohu~j9AB)`WJ1Cri`<>jG0 zc$lBbyl)#18ln~&JuAUO9W?B6;bO>WFxq1w6X04G39?@oo}oG^Y9-NLj*TlL>%)zl zFTNKz(zw zDfx(xFTwiPH*{}n5Y#fM&Ch+1XMR6P!p0^dAKXI?E-7 z%?`>m7h#xve80U7V7%!h3NKScpBo}EW!Lu8@>~BN>0RC?AXoRcZ#~HlPd4|1J0o%S-rtii?lFl75LfkuffMHuOjiRr7SXQB+a9Nj06LEuHR8Wos39iw~(*5~c;Pl;d?jU&qq@4fB97(*0c7j= z;9oDFu@SSfvWnB7v>$uVgbIlKc+>v;9Nxf}HWA3ubx}zo8ySQJk%gk6eF_9Y0q)-k z>Yk6?7hzXsQI9<*ZzES%J^+7kX(iHVZ@1+478<9fl!?)Vqrb_0pIvM7G9^*asS;x? zsVkT^W7Y|ukW)-w6=cTIyepC;B~j2ATaxX_R06p0-+|UT=T*H{e{sYXb!@#A1USDy zJI`e=`9xGHp)X|AYVR0|%L%eD8cdiYMZq}d*XmcMC=|m zKm}MeKP+S^@tS3uazL3d9%5q<$es&;JTD4rFEP8$-|4&IfV)b0ITJH8I{6q8^)9mS zV8ASCscbX4L6hrArbETGi~Tu@2`nru%L|@M7JJj>ejWoxExy8|qoZf?-rfzmQ`}f$ zd6ep?vE=D9oBjsTiQ`Al;5E|T>3Mv2aRj7c2KyI|Wj?%_!vrJ|OQXA|Vz+;PjZcZ* zdkT`_GBPn;IzMlh`o93ObC^M~<)5`aKLeO@B1dL<`RE-BD=Ya(1A4o>V!}Dj5`eNz+b_w2sCl#pK z^}=iI$p4(0fHg`5B+6h{NEa9EE+QydWYvAaSa1Jr*~5h~9CbX%&uUsO%=D9cb$JjA z;c019G=6_IP5e%u*K%JDi5`B{DSWg3Nq6ty0${f2bH;M|Q8xlP43F})KbE!~Ym zz{L@8^7ogz>%?>2<2F?MI={UPaUW3DpTxz96N=2v%?UgIMvIP)P80X3={ZmQr?(?b z9TZ%yzJB+W5E}S2iliLi((qFInTp(=Lx4y*z^jN>hlz}X|1*VkMemG~gFY<1psAla zYl2yG@6|&y&%MwP|8jVk0@Rq2*$^eDp$1bWRWKS%!7(Tk=b~aeMs|Xf_4W0LHElLg zR#B!7B>3XH|Lz`L?ADi%@)cZaH=3M3_!C zy}<07*QiVpi{>@v2W6R%V>7E!Gdi=kVhqOA&*%nQGz;cE7|>DKb&6n*9)*{`>^W<7 z6+(3Ff_BqzEv@)p$+j*oj~hJ?zkU6Rgj6_z$gzEW+v4wFdStQz&L}~)EG7AWS*FG+ zWdt%GXutV%(Q`D-Bh5b~$+nVdd~Q~}@t}5Umt=tx>Jx5;trD$nHir z5(xoogoy~8nC>!TNm4Lge4$0Tt>w`iu^HbILI2jhx0(W5fMCnG86S|1jx9Kf9d{8r zG$wX*bWAt6vopxOU;aVq$S5kBYI^r4sj9kqhBF?tpRuWqEh%qKHsZc~QJ0lPtF#=% z*=h0b2xTI3`6$@u=~j|+#KuCR@YFTb&oDrhsF9H5d;7di_TU#1tlkn1<-!}NS30O< z=d+1%aifKE=K!iMIvRnz720~U&aGyOkBx0M_*uYX_lIR=TzdMTs?fsPS;K88XpLfJ z2x6-u_>*GFAqvV#foTQKgDA*Yzih19q@((#NEB+xuzH&`fV=bw_s?pP>m3D zdNM5*-uF3WjQlpy13j#hN`~^3c;7J8ruy{-t;;czc0@Yd{>b;`N&R(Qr5#a)hd`8< z1OEBivzi9ld!%|2MiXvQs7}?hI^z342a*Wxw_BT*FQjV8A^;M+16}hA1zY5hR;C>8W;_RmfCc;4$5=2Pf=O0-zs{c`qEP@7 zA9{{_EJIN`!6+H-wj0fe*0uTZ<42x=JKr!o$}qJ&1_p+=2pVN}b{5Xfxr_)9@>$W@ z0-?T2UJN4H!y-BL=&|gnhD01Y{H*Ko&W%K&Ru-GaGj`5*Ve>3|5 z2lQ4gukt*fg4Y}kG%c1V8-pJ|eJXd`*2?RVIGcV6xV9sZkK-~j)b#YoB(GcVCkq{0eMWrma;d3cKuJe5_BQYy2%e?6f-FhkT-_%hVT-5 zEZ(}m=bE~@J&(Tlb{q!W`{0IUyzSG{(n>~9cY=0^?nG_a6O3cg@7(zAv^2nspCUyv zK;nGfHKwyK0Oo8?$j1ubfWrH(=HB?We@o1YZHW<3G)V@KqyB>C|En5UU+r zG)SN>PumH8pc4Lla(ans7K9as`wn6YL&?(PX0u0Sn1@1JUzrKxX5p(rx}fGSrZ^i> zJnVJJ#J%|eNXyEi2OMjdNhz&_Z}RlC-ou9vjaq8x1>42GO;RKFGA^pS8FH7Eg>Hy} zp}J1+E*D4Zc+g*B+1k$*q`IExI#7VYja(Se;{khG(sE^k;C~(cPf}@|WBdD+7JW}t zK~|W`BLnbS+r(t3%&fwq5fDA&HUckuKAc`&-3B?CVS~F^zs&kt@dM@S6y<)_GL8~O zX*V0_W#0JSN;nhAPI(Xt4=h(mW1tR8-it9OM1p!(?!K!7_|S`s)tLILokB%r<<;}n zd%rI|yQkj%$5Ezb0r&o(Jp5=rUmO5Qlb`W7Eyxr+vayVjmBp`vg7jhE&a9BeMb5Q2 zKk=I!nNf9VLd9jL*2gGbITW?GCG)<#0K^fQfDIAA5$`Dk-XU6L1dd7<^?XQ5I<B z9OyDha@yJAoczLKGpi^3Wb0qu!oLSSNZmbPF^-@pmh-TvGVopd;X))EE=7q+&hhHi ztD{-lQoxJ`&Dedbt4rR@@IlbkNgrB4_yV>pySO;I_5St{;Lj_s?mGq+^jDzmo)5Uc zQbH&e;V#dUB}2a3Vb+|pfQa4?Ew@iQm~%R%82KdQ_HW3hvbJV7+Y zF_>9p%nx_sx05}5QbZ70SG6>mSLWC#vKlzO-=7yIgQ$RzaKJz*oG8>y$jBIia5pXg zB?F~?zTWlC&2bO8uuBfDE;|+m##m^rXr9enEJp>UG69n=HvXxzC}AVI%Lh0P=zd_q zYVeZapZGS9U~MRSRC!bLqL1$lSDF&wVrc=jhB~!OQep9Dt^G zEo4%3Np@TSdH|wb8Q^Z+bCDYUvpx{{;X{q&*sovB6_!%VM`3|!MR6~&lsF?Axi~=1 zp3z{<<0@(=MAiIWi{8?0dOsxfI*5^x@$0CSw>~9c8k&n91_CnM>FTew|II4xTv=wud7?Lv~ zG6w2b)hH~|{%A`89PC@U;wbepO$yx&sDSTD21 z!_1xJ;e!rX1HiP*%wgZzgey5YIP7|>POz8(u?`xJ@88+X&CN%wj8b2(vH~oQ6h+Ft zQ8VYj0BXj~jW5tTI^%f^&R@rj)Re4W^Chpnz!fMQH z_iIBj1zA2iePiGcrf-!6CN=Q1dNpc-Z>P?Z9J|YEtIxle5)%+yd;tAG>$678+(04{ z5-mf+feMR({QUeE9U<6&+=&B16a0WXTL;trLT>r>cf-PRMM3Ftx-Lp1_5+u$K8ni+ z5|0ERANo0aqAeLlGe^0oOI5h}FKw86Vfa~v&)BrIwH1|=I?T(_yoF^*igk8f2W042 z(9vFuWJjHY@IhAFGt56G8+`{f@zW-%BM+u zXgCw;PaRl*LN90cTXsU%^Hb;p$3F7(qRUEUk6!FU)Rqmq?3-SUtVYvIGX-Wq36wsX z43s%pNHpxG=)FT9JFtaxSwqF^jP-&u!$E=l+Cz!##uV!LD(l7~jUsM#b-rX-8##`V zbj%Ag9H1f*L{;~0O6d?)B8VMwW#yAI%2bB?E_Z0n>%!PZ(WTxZF0yGCHUGm8PNI+( z_+BOUaDy1xIgg}^2K6`;eQb0>*sKg=6kP@%lt8H%Pa8N@W7)j*>n!~3w5S6J^7bzj zaw_Yl|9qU+U%5DIcK(L!Zx?MD^M1>^g}WTQ{{hnF$w2V}hFp(t!+~uW6)0 zFpEV`r28S2N}RpX{|BRdmd?QQO|Wc$B^NIlrIi@$aM1~UV9dp#44L+Q_Ld&{1P$GD+wtFR&NG-u>rLeo(SD$q)CV zfL@zuUUeN0@AxhH!{2QvqZENd{>B~AOh*g=c)g`?eD#9bnHlClLU=VEYEdff2}RW` zkY!GU!TuR1#_W%r1@m>4;=QwaSa&)4ZyV-Bbsw`CW4#9)P-u{w6dZ~|Pd1LrkKt0$ zfL^FdKl_!5wPJIjnI2YxF^ZQ?5>ek7Hn_1K=`<+;z5uu8w|thb@Id)&i|QTvV~Oh8K*mref-IbN^Gc10 z*?S`KXoFB4218yGZst93C&1aGp&@ogl(`=x*~1Z=DhFC_MJ>?jM~HM~L_QyT_Ipft z|F*12d{PGu?6(NSezW*`P^k1$6t9?s--b9lfl9YLPyQRK#_Y45F<}|MiLvG-SOWkC zH!+A0cT4mO&hwo=*Q~emd$0C>i=cID=qO?Z>7!sHGfZMPd}8HQ^2cC+(vcCXJo0#i z&Fv|is=kzn2-0OinppY{JAF__7cp^2h>0U{w0L2WCnvJ-+AGH{=%c`=dN(zYuDzAd z{oihx{i#IUS0SQS%?;#bW~|32|MCJL&;$^atnKZYczNTZxLcAzCj0CCsD2lv9f{0} zU?pfw4v9JEDIkRGAziw+WmkAvQT1CW_k1sq-|QCGdL7hnP6HK4QBe^jJTo4${KCT3 zu1LbIQ}2zT^nthd5(;v1VM3D}P$X1WvkAde7F#RSD-vioN`cVCFFS{WT8pNI(Di*X z0KiGn=ulnEj9vPeiVUjX89Fw9`jw{+gx)?tK}@&!i2~gnFgX$z%cy`3c=P5>6Ck1% z4QiLgdd#aee(9ikSq4#YBoxQQVXw$;r{^Js+jd6&>On3+A3JU7@7RIzLZyy)VUV@R zM;XfidAMJBTUoJyTfDydyDfIPmXIdmZUajq#M(=<7mo@FQ7f(tXm9LnS3940+L>-h zi_^&vrafgMZ>`O``uoRez2lS$lvW?6hMqI0!5Y-2j@yf}tt>XgL`N3{`p3qm``K#S zzW#pc(y0WXLW7A;JGj2PK1@hTl2=i|U0GQ{bOJy`X!PZ5U(WWyFLL$KtM3R>gREsq z%&u z!S_RG%WkXN{^zGD==3{qG^vOwG)8**;OB)pybkkEKwI@=P-y8V*c`npo#6`-Q~x7g za=Hh52_P}5ipolIfHDDpK?}gZtLw9UR6oOD*i!wBgn zLc8~CdRhU{BO>V^Y)hLHfnenXdVbM)&&xmSvhwnCw15TRDpiM4;pSkh&Hcte;-Vqg z9lImfB+u6{3a9ZslHzK;aRHr+*IJxB zPwF6tenkja+X|p*)ed_)ym%I$Bk?vDs7iZtwK1iwcb}K;@6I_)nv#L|{AAvL(*3gN^6z^5_wkH>2VoJ@`**K_3hqhcH_}4F&zJ*m$fakU=jbJhkng zDdB47Za_Q)N~i$|fomlPj6iX*xdpzARH_g2IX|%29?RWZY)%CVQ(Sz!G$7l72!Rs* zfZF9!3`{w2@d(<_j=#JH_KxJo7Vic&;~23M`quv31-Hv&7sMqFd1H1zcKVcs!vJgt zlwkrg{RK9-V*<0PIWR*Zgl(j~rPeQ?3vfMb zJ~ywNDKju`WiV{=;;miqK%||u&CN=O1(A1AQTlaGbV#7TX85Gh5C~ zFm(o4JpyibP17yyu`z@!I9*&SW< z#lp$i{ep)Fgzlo(Djp65(tLS_CMEsGRd~7O6!~Z{MWjvqoY)P!x9>-S-@z0y>ai=7 zy5TqMKc)g(0qp$2!9kVJ{E;A3P;Up)5d7D6Om2u&WsVzYKIVD6JZczf&N#Qw0WyFq5&0SS6_O%jDM*Gumuw z3|99xn<;Yf!;FA6!1Bpr138T9M%F3bf%J(mQao65K|x-whbTLn8YdDTj=CEYLN6I6 zHdG-w*@l9dpmYXDAlH`xoiYUw?!ZN~(y&{CSZTA)FwhGF^G9BOe#HJ@<{t#b z3!1r42us*u%XLZ5E{JV}m<`^n)^%hLY9riHJ*{%@O+Z^Vj!qb^_V1FYrLM&<-zNI>6?7Sqp#m^}$&dM6eq?}mxs_!8M z1;rEB4P1MBd(gb`*o?CPZ4YGWlOG1?D2>-_eX(mtka==b>ZCl<4T~S&>gO(?g{c8Y1|n_& zd%w0J;&=;#*1gq|-EzQ~ta(`xs6w5-d|dsU?UOP8neM)4ow{fH(PzqmS@lOIpaVvc zTs?JYLI}^bi{9)&(FQti*i7YX;H_Emx_%sdLIjG5Xp!KZtMpJC8`)}oAkrGsCfMNm zJpsQAqA8*cY92NvMox?lb3p&&a-T$E$HufKOHGFn?g0Y>gSzc+sbzW`H@AlYuWv4n z7wbacNN)9IVWEv&FwA}m>&Na;`w^EFp7R~aAtvrzZ8%%2GsWsic)m5SVY5|Gokr*b z`UtKgQPdiEN&N=Jj~H5YO-%yZ$)e1rrZi9yyk6R)NO2R9CO*!aCvFkb!fP8JJ@(jb z(`N2chnJ3OJAU7k?Zghh_G-@c2{;)c2E^wPA1X<#y5#LgbO^~2Ac9=?tEaL?z38l1 zpikNfNhe9=zau64x3?Ll*$NoA#JeDEI&hSd5}pa+8@F@)6T{)ikX`^ctZZ!P!6qHx zYJ9)>B_LyG$HuWp6}AP;SRTmE)e&Xi3(`?@|J=lElQMc+P%BX;7!x;AwH%y!Fgsw5 zMw=HEw{hd@!2$F6=5jr$AU}UuxAltnu<7K1wY4=M=Ro$A)lT5aLto6SX(maS2rSkibP;%F@AM3c?0se+>{T{zm{_7idKPpeHrnU22m zUKLrDJo>73T#NI&?xl@odn5{yyij|xd=bep`dF?bR1@`8^t$C>2*NM|2QCPp3fFLR z(LKR;7D6ptFt5fODgBcb2)tqb{D0(;Zm!o_eVPy4NfXB+Yb7DFxhUyaIgp!hbu>f^ z>(q%SbhHvqd-VtGI<#EO#6ogCQ|i*pZaXFcgW;~KpxncZfog_Q^cSAp4r)Rf2gjcm z=Wff4@Ul=t*kp3w|HptOh=F}C6@;fQUzx;R$-}^u)dWp^9=~%VAVnt&`cgc6$nH@J z-yx|w^TzF(fQdw?a0HnGU7n1;Zna~lYg4sx+`9DQ;4I1;$aN8pC32qM>cb`$87L%BFRvCTX*1tZx>y7AN zP9pTjZp#*>vFJV^9l5(8*#Sy=43x|wHq zLdX<&{w_dA$IUEOrS$Fm#rm4dYY(zCZS~S>ldJEm6OWHqt#_?n>JEI-2v59*MH~d9 zXC11LSULV3szwi%Ky;HtM3kZnx^WCw)MHhrR(Lmkjh|`$va-nxZyjDL8^{Nl% z$vn=atwBzx1sUI1kBdl}LL}UxTF7g27Db8{(YfNxM+W$osXwT>xB6X%ETry+S!^)S z&#sfjj23;z;roK{f<|;X2L6wspBchMo9z~Ss!GM}|get`f1GN4!=ptSt(9uF;5}-RR zhci^2oy$9KYk_o(0P3KIv=60=k^=z(VDX%yqMyU;q!&24oU!LEtjNoG^_|zzpDy0* zk<}YI5{8%SaiU6^g;SD5JVXUKAbDUO6aZMZ>8Jzp;5-0i zU;Qd80)#P)Qs_gHsXuw@;E&f|{5Cc}L43iMDX-<-9zbDb6LR@l6$hBui8Ls`4z>y{6W3Cu(|Cu`i0)CAM0F0hT$ z6|fLVnH^Ac|Kdy&eC*WK)V@_zjODtD{T-13y(;g^AGqPA+@>vQTU%Dw4C3N`ZM!5F z4|4U*ET4VTu6?H67#eW!$vDxK`%qGQAZe?bf-!1ONsA*7c}0yIP@&~j9`!9~D70V% zSxyI}p+E|}thzr6mplR}-~=pD_LQjLWfF(J^vjpox0+g7zds2;d^=fJS&<%R^Q9l0 z#57w}eJ*~lg!?EUC%{HBEpzI7x&9jWpLdqrIe!`mWvAp()4vb`chN<1ES!7U2s8sD zR{YIH(Ex&*==9Yiyk`HK^-pAh`yV58Tm2KXgUq@hkRPd99I6$Y8jML&=D3Ej^p1=& z8_~5Ue}~F`cpwB|H^yI07bb5*R9P>_2k+ocT!6{_H|7ns7{huOCUvo6G}?e0UVwjF zOun_nu#iOXAz!v0vOUqBQ&LGzOl8Jgq%Oc8eX&afE!BhR2yFxfA`+iNhEZyCK zOsvEpVO96qven{#)p#B+3pIvTfiIy60IMQ051Zw)n0-<8@ju;hHqzs4{(LX(V2K-K ziot5gtb#GfUr`D-`T85r$6SOXxUK-x1n5^O4c8(Ql33pImwR&n{&wF}qbw&5$i=2W zX=1B|b2eDUzjTBoGHgU4vEYvej&{14bYCIi8|fh;mlH5rQp)y*2;weZe+*thgKu`nSGs!%+w2LvIs5?IIN{cvg%657=H1F zkRg@n*%#fc@gRAsQ7W3highP?(o7wASVb6DnfsL4? z?zron>QXi5hMG?ySjHL_Cde)~t3P#-s75`p72e84k8hbN6XxP479Ow%mJOM_C>F>t z)`mGDissjKBV^VB+w~dMKQ;Fm0v{qqVC*ly3*&H3K(q4rIp(-v)NI7|a`==U2D59g ziz)RjSJH1w4Q$&jPxl&7f*^=F=k!31B!A1OcysXf(5h0eo}c~0mvub+qaCV|?H2l5 zEw@ZFm^JqHO@siMvbbno?sNaidODcdaUc$3DhEAEPz-_+4n) z%>Ch6_OJXmHbYO559@Yghmb_bC1EOg$j?Fw)}~coM%h-{c>@cAC`kl!&Or2p(fW)# zItuAOR~jN49+?~sGs%9WmGE=*o(*OBU^|;EohStFvAxcd)7+br-tCZZ8zijfeYW;T z{4{2sX5{mtYSycZ2={;}X!e5RvkdfkxFQ+@+MfP%8JTLA%&#>J=`pRPvW1$* zc^y%swgSi|DTHdibds-$1DU%Vg7V|EZ73&ftp0pW*3g3m!{o3cU@WZP6F6gy!`DR* zmxh`Gqr&7dT@?LXE41NqCfI>ek023BJ4)-h-Ke&0MS3s?oQU{HCIDd&=6|T$G$Z9x z;5zkDX497MaeZAj;p^5jkJ5b~7eK=VQSjP)!Tw3rp>ix==WXBHZ|^)v6I1sM=tYe6 zG9nV;%k+r7EFW1OLZZwU&l$26yHeD^idij!Rarvx!4ad7+?Fu+!pqIQ-M!RK5l-k+I0c4XPy`2u z)Dci0p!~juM#IgAdrz~uF_vEw`a_qA>)Q&7UF;a>(Z~Mmr4l3ardsUwxsp!8Ib&_g zX69Q<;L9k8fwjl=Gaq>E+Zu2(OXIq1&MHPGpfs_OJ}Hy(r}4Cpai}*-=i;bTbT-4| zi70X9Ljalh=ReU2E^tAt1Zo@2_1l*oWR(v}dxQ}vPDt5VacswJ2bmZwRZU8bhW{BR zU7B6ZEgKSv6-*F&TBce2WAY#}(2GoV%Bqo>@ZFP$HhIK|KWYnb-uN!4egCXt6id~2 znl!**^oU6fXUJc7LHMXdjF;z-s?ANJv+qvl)vBw9OQPsKYcN(vad4Hdsul5#e4*ApNc z2L^5N6zEm~Ito-k#sADL`4fBuAIuE$$_s0TQ11PlhQ6-4j|ze}Y+t&1HO`*zqhd(i ztoF*~OE8|dG;F z63EE=tyitT&wU>)a5wtsr1>>zn@Ta5`%7X#%U($dLB#^XC&vG$huw$)vn@K2x@-7j3->VHmIFf=sRX3@`8RJpvy|ii9w>b0$SE57(+QYn4_XRu)Qg z!J7Ej$Kafm0V=9T0r-8Ao7}-po1q_8(owRedikV2i(t^hXwXLYQCxo9z69wG>wD%H zID04Vb9OW$gMT_L;ku#Xu{XU^y_>AUtiqgd>zO_J_ODXVKw6Q5kZ@;G~LlZW}HxgGX2 zp7E|F@8fR9{kzUjs*JK*UN74%=7sHFl@cwi;`J2RmUcqTu@*0S+FRQX-2V(u-IXtK zpB`ykLJSMz!j;G;-1Lb`UQ`FFs@a|a^OY97_)`?Nc*;|H=wBHP?1+wfQ?7W%H`UB9 z{+m8#G7+tXf2*14T5$9oKt0^gn;CDy>hb41pbV|p?lAfIAqW)km)Z4f=VB7waA zSNCY8bp1=f5ha;;_L2;4dS~f)x;A#|y#4DeJ=*0Ih;ngQ@=z)R0Vcsk;MWRk!v7Q0 z!m;r1h>suA(kVW?7tv_12=G$}2Y{f}@0N_j6HQeq6(fFJ-u6#c6eAJgavx1;8E~3FGB77; zFxY?Z9>8?6K+$C>v0snFZ+W?^S&zCmy@sW;@Ub(bs^t+j#h_e^_%$ceuQ2tKOA+I5 z3oLoN2fikEet&UUe5uDvq9^u0E3q_w`8-d}Hqnf;pbarXz(55JY^*kLnC=caSOT*`m~82`NAmA!@rixI!ukak^EG3Ijxys#Eca% z7};9egf`z;P?^VmsrG(D>so1N&#&7tla(=uxhmiJbEy+u(#_94on>Ez2${Em9fn#U ztsAjV6x)!H@ddGvBVcC#3t$fshyMJj0%V8I?`SyGI+XVk8ttxMs8?S_8ySnU@ur2B z_HBH3R&-zzOD?vo`Dll{{BA|}$x?YP$&g+6uep_<5A>ff?!p(tT&4fUUZVDlOg9HSNX<0l8UN+lJ~j9cGOS`R>_XvbY>hmDj$mzkX;U* zNoN~Gk2`d&mmUQX5wJ17NB@K|OS$&UB)2^azKjm?=~ef%!%DNG(&fpS_D6wUQyUm^dGWUg>F3Y1NV(#;R{qv&QQvD_mYE6H%`UrMw!rrVJYx7+W5HJ+ zW+k5W)7RQ;YZO>HmiAbLe+12AG0kszsX@&d96(PUAux`but!Z~>z>%)^-|)fP2$4S z$hgcKn^^^i7Pk;_)d_82!Mx16fFM~i8R4eVz7-D$s%a=m4>QGCEln_1EKfmg$Y?E( zsjjvf&3~pBXJu3BDg$qV%3>-|a6@`$=#zGI(NJ{3&9)*z-7OsJD&Ah3YZ1J;m zL2FHp?^EX@|6sCfKw%r;>tV@POkJFlWihrJYRvV0%Pm=qmYL+ij_5|2T9CbfVI=>S zRapv8jZo~`4mQTvj~vlHQt44i3vi0|pRr8aJ@k3EdpReRi$Knn+L%Z-#75#0DW)oY zi~+Bj&Kr=E4VT**Z*JSZz)@Xy3N;l4dpQCJ{Lx>&#*b`Kjq|<~=dW>xk`pMs9hB8> zq&(?sUg&5K3;Lz%P#4bwX`AAdrzhs zUGk+Zm5U>G*U(NU2||Y+dh6rdr&>AJQ>E=-pBo7O618EFD%dKp7VK0Y=#KOmyapj-kSzea%sl!pDN=| zr4543#pU5y_C#*I2fl%P=#gksH z#KN9Ou7e@ZaDFh?Jrg*1{Ed0~{=rC4zjzm z1@F(CgMeQkH)<|zdM}mN{&fQ2m@G)#brAUiIokxsIT3@H^Kw+*jZ|KsFViCyEHo21 zRxeC0p0S@|cG_^htwg(ePei|`{QCyYwbeKq7db)aUr2nNxgGtgc6QTdhE|#W%GRT! z!9pir*Qv9M(F!Y-2J*iy%B~3$UAs?9QJt5%{g{;mDMFmV*B0!&?bUMG0`>w;CVxflpPb?EsD9w<70j5$YFqo* zH~(fdeGPwr_QD(eEebF@MldRZ7!WoZmDLQ51K`Hjbg{S3bzT#R&W3nZiyKUUHk>bu%`9JDqlCg~E+8&CMWLU_1fN~+i8@&M9&f-AmAH2xGh zBg;K0n_^9yjR1vJv8yUk({{($&ic+RpJ|_X_Z2~`jxzghtLgC7;nc6mZI7s5J;!gW zXr0V%LV@iE)KsNZrJ3jTQKjsyHIg`EI9~ksn(>{@LY?x^G%RFnSUjWVFT0Ib?tG8F z8$QvQ=*m=Ux3a2!B4#OTM3yl2vR5&N?veOsNcdBq;<8fp?=ij3(Wr<0eC@5h@Vo$i4=c<@^Ubz9E@l4Fp)t#5sZ>k82maU>fehy{j zr07>@%%b*8$hD#&6r6o0n-is?$<{0+OM$9+*bsq5=~Ig~H#guvCh$+onRX+2ClPq- zD0R75vVTQZx0M&`gPEh^wEZmH2LJd+rVx@B><~HJR8}VLB4S8Fk6Ae3JDV)3Jr@xn z+d(1JFNUKCU){I&u;q$+a!9RfuD|BEzg@RZbko6sGJw*+o%W!AZ+8t@zeggp27{+- z+BBXqvN+f9jmL0hnw1!G8?a8zfANvz71r*5YX#0 zuRhxPAYZPeV`(vjGgMwJA?<8dow|0(a!P@>vB7rbzBsd>t9TZ8o>A$e*WkyBr2E=P z)xm{}Qr%$3mNK?kG&a7O;1EhseD7-ucwYj`(mEVIK^+4cz^hKL{vM(gyj!%WT4QU zYuNqm?QH_TmK<|f0y!-WS*AuTV|{mbxw^#)%iOP% z(<-C7E#~AS!ib}vq^(zT6I9-A_BM8!6(hMO;y!V{{`Eb^H z*ZFW*tVJwlp4xlg`+4^7x~R}8Uf57_jgCvg5E!sbt52}`oNlI89CiGfXJ2shGR8gD zc#fB+t=@mD8n55>7g_kF8n5F!2w zgY?;MS1wbLXfTk}c32-mShu)n-HIjZK* zcOL>jqR+*h6cT1>;>Bg+#FTl2H6ftPAfKt#7~U=nX)@IdI7rY{Q;Y8!LwF>5BAVhd z_Y`a9qQ(EFcAf<>V#b%NY)Q_AQbyq0JG=0YUl?SjD5tjx7{x`6r8F`23bN{P-n-nHVOm z5cYsnHFTIkg&oK3u8a9Cni}@%m5Y;_|Kh1a=|zIi zfTyJ45tTtCRnt@WyGY&BkT&FJsW_%dY`G;tL@njgy`!e5t<gB5)lV}y1fWJ_0jmjs!u<@C>qm+_7$!2epK+V4j!9_02%(nhug^vUy_9bMWPBx$tb1FQsOd zVTf4|GwN;Yez07}RE+3g%{}Y4Y~k2~Cp$y}%ob0$3f@0&zBP4`=3n#pxaU+rS;RdU z9J=xeTvYeIZ@ixvZ(8t4=u%aAl++4Kiz&iSbeAp<#|QqhD7(yhl~`5X++pj?k%!wz zX7`e!5Z>0eF4oJL;VX037-AnOsA;NAZNDpzwx4zUyIjq8AeoAE(Jmb1N~(;2HBxu0 zRk^joNyV9_i%$f-A8%^3y_NZPc$m-f!s%1m>3wbHDAr#=X@B?zlb7$zucCy=G;cP4 zR-Sij3-$^pnBKUiblHyCC)Bdh=N);zzj74Mv-J;E0k*`yRb42rlnkU@)G#~sU7zsO zrd~SvyM4gL)#)um*ZeD(M78fHeQAi^o8@8U#(SJANK`^IB_|=p-)|EaT(0Xcl&Z~O3 zKlsvh_HVvFFhk11ftLMOwwZJM*KPh z)qUEwC&`v9-^>l^i=n$asn^g<<_*L}nN0o*YT-@T^t{^OC4^^+{}RQyt*?E=zwbL- zX`J+4yPO(io$x>}Ed(34x+Ns?_Ag<>+*oaL0`If+5-HW_NqNIFG|H~m8RMt-shnKu z^pe+B<(zYK|AFh&kvJ6t(~qfBtX&bJG9~ zuHcaD%lNq{`NyQ%v@#2ai@kT-UhXBAVs28lY`pzEgJM1Nn?Y?qwGB7%Syp*$Mq`)a z7G1qFnbD&Y6y6;4%S;!mw`S&jkaLcL(LE?0b`-5Qhc5o2Jy2S9UfmI2d@DRk4zDXD zdEkVwZj)~7gl?7EfqEDEZUe0^%fN)Q2cw&n68$ z`a}8V;lCeE6-!MA>SM>K(0emxdahGqqVMZF{-&R}A*utS>CQEiWE&G<{%MMN^kNod z^X}g4oJz$y(~<0xTD0{sOqlTb$K9~Xl3ZL?5M94Vf98lj8|FdgA5g0fxZOHpYg7DS zIbfOCE)n;`z}>FwU)=;(h(XJmxTW&7mARuAcq3(I#5g6+_wV12VL|r!>QJ#t!r=Cm zE}R(ve!x9}fkj`O5C%eyqZNS?du(e)zUyAboFt#6`=l`p{*Ls6Vl4N!YEw4$?;EXn zIc_tOCy7h99WZz~t}rgue~O^j8LxPITP~;liK2L$qUfj;Dg`SS`c~inO_bf#4b$hu zyaPuT6fRE-InI5)cYPUg>z3_&iM`graA`H^VsG?0NDiE0QOSjZjmaJnlVzB%9gjDr7P&&Y>mVcA9l?}(his+9Eo`rTb8 zS}wh>C}#0bA`aKZT{q@IO&|bkv}6qYS*oWHp%t+H$)4i5A&R5UT{b7X!5spbKHlw z&63F|xsoKr$?2JZm6NURElIaDv*Go%!%fNh_yf+|PNeuiF}@BbMSNGQ%#|APwGBoO z%~O}9{nHX=Z+<^qh7%r~B_K{#kUfUkxm-5#{!oBT4l#=s;q;dJwj*0)08y=9U-rC2prV z^*NDMzGJ;Dp%>lnw8&O9e>G}5q%XYkvp+4%1p@HJ@p;z^V%8@3BU^cYdFTDpu0oOxwCJu~(5vy+IwI zJrln^Jzx5<8MAF|#Ukv6ufN(HnEamJBhf%YI^CB89Z+qV( zY)>xVW@culUN+oTu!--bF_G6_TIqdot?n z4WF@(QwNsZwP8#YF>zqv)xz&2e)nwPep^8J%$b@&SpUAZT)sda%ja{t;}g3y z8eatDM|&u<3KGW-d8B`AvlQ}Pq92nOw7wd2>@~8viP12*GQ}9?r5oAn+TB28!aev3 z^I~{gR}8d6cBRwws`e6PmL+1cvG$kok0eWazEQ;jk?RKoI) z(sb?)9=faA-3hMr zWK3~63%6v)P^KhfR&Df;PIqGBLmYy_Npr6*Ddg zES|RTZ;?F&>HtG+T(BT(F6StnzhvrsqbkhhPr_PHlgU(8wW$Bfd){C|)5KirlYDmB zyEhH%UR^OwNAcR;c^mi^azGz`^VtY}=$@PWz1{qdijk<;mJV?oWBTYN zvz84r@0U9Q^zD0sk#_fAw=Ju_>3^7CO=oUw{QG{`YhK=t29Gj-O$v_?9to{M|06>@ zYlpoa<=ekj$k#n~aXl2iic&L-9hutc$tX>C%|4Ahnnh!3Z2P>EeWnOhCqJ-U?pS6^XI|`gR%qznAmI#i zx%$FoP3*T5b6rW2|N8swk28;SHy6k1{!B#u?30M8`jLTmNtM7qM>Kp*P*oJwCm1M0 zpemS~{tWHW_>;MlAU(m9>#EtM;wEt;|EmUqafu!IRgVYKhP6tEHS9_=CykK}#^FZQ z>6vD_LTeo;qmlkZ;<|Ft-XpgY;8_76fJ3VB zBzX@kilg&qkc$=Xvg3~@At!fX2%wZTOqFuCD)-p4!IeXx({VxV(uV!3fH`v+d?E7- zyhz7y)pu=)`3X|$ke!s2wEpjGsKoKsw_^xin@%-)1J=B^XhFi(%3X{kiV$^NwBMdp zz(rgkTJPZCP=9v%PyWG!bLx6}!OeXhJJhVKNPYdQO^4lfdeyI(!RwQio9nv5K*<(B zxv>)6`^Igw_>?J_Qnqrr=fo)0mVio@#KAEF?|J}b3-@Jd*kg1yQYgbyUs)ic4w-e4 zt+hMD95Kw?J#|^y4>(i2Y8Fmv+2BMlU(#4$!y4=qHd@X(8R2r~X=!ONGs|DNy}#LX zX1p<7`+iV$Da%Lf=H5UVV?JXxH9)w)3CCgB#@Q^d`87X%x;5P_TuJVp(O{-2oE&|$ zk~E0Kj29;|N@A_p!j@%8*v7Aexq&UfpdM@|pVO-^I#ghHcCO#PO=@RXbKOlaKe0#Y z#7XT>u90{zxfHQvDDQ4s3>&^xA*krB#)XkPwHe#lB0Vhk}I#CxzrLrA{FY%nE@tGFrTw{rmGX z0{8_2fryHV5^u|?qX>8pE$MM+M*tBEki|@90Y+|yhqO7E(p;P#HrFN20VHZ6(sBJEkIF{7}&A#I)UmkX; zS9|~_O4+amUXYi!^_)ub6~yIjfN1BjE-Ndu0oOx=c19zw@iU8-1ToxlRDm#3y)Fi zk*?K%s2=z|9o29Wht^hk*s{P&0QIP8Xf)K-`z~^#n6yRr|?%iZsq_;Lwirp z#`mbbqPDgNkXwxxdYPbn1QI?Gr)72C*aAB-ARah4ICz}R`Qsm|sMK?>7Gz|c`&ttY zhL@cD{2$pqa8;qVvlOzDUYlTtyM&ixU&_F>WT~&DL<*IY=TJXK`OzZ@{l>9RSt!-7g-P>Evy4dXfX(8Z4C&pbYHw*%yd&*FD%MOknIQtsb6g{F+?5qT=bC z7E^nBbEx_FpYd?2(t3-#T1Mt;>pyiyIgmbkZ79 z)qTQ6TVG#4dSX1@F}98Z$JSSeEJ^!mY0IN^9bJ$yiU1hJ0wGPFhIfkc@_m58xx3eN zK6wZ!dysFY+)h%|*MFt4)RU}*_}fe|2m*!JZVOuhU+@OrVP}^IN0rC`PKRtQuUj0p z3`^cUmGCi^L*)(V9FH~Jw!}TyCBNEl0`?V4eL1;rF9(84R%`XkOasgHi~m6Tf4NwF zKO;9cS9}0kvczaG4|Lk4-c(&6_H4j__6IgufJXlotwehrmS3Qt&;?3M3mB9?$*z+c z&CSh!Vs(vIzDIesLGt*qni^^p?T~qaj7%F+&3A6zRD>b)c6Lnp$csBXJZh@({0hvH z9@@|yk(XE`^ok8@%HxzLu|I$QjOn-12KcVFqq9>Bym^}V`1l1dKI8f@yhYC@Q5u9k zGeA)~2-}I(L=o^*j1*0}FEPg^CKiH2B$tSsDeunRyPh@kl6G+H9WR6m35dzjYbLNF zpX?5a)f^%UmLkxb<$|%cfD@Tg^Z>;)vVs2+8G$#O1KpIcQrvLm@AjMxuxlO9`Samo zRRNYXHoIm96#x2A_vy*MWcQvF$wJ&%w6L_43wMCII1}^hPHgd7xbJc1`wt(qp>^8q zF8B8s2-pe#JJ_^4@Wk0=vt#CN9m4?r7}^c)%>x3?!;X*D1@1?n_IUxw>5*`DcD8{= zS76fQ)6v_D>`oHZqM@brDEjb|u@qMz+Ku%7M}cAYgW6#qdLCCMUq$}UYN&7yisk=T zPvl7`j^?A0AmvpU-g{-V2`3_6&{n6qL@U@9_i|!MFufH)de?tu6wv^^o{RP(*i^jg zuv(g$!1VwFSq^M0cw}V3pzBKaB$!+YkE%ruhb#iyrlh8}q3OfQ%DUK_+RJEHFXL%6 zl%11P4zYJHFE74$6$$Y50E1Ut7s@Rk1UuhP2$W?LOdlMKJzzj9?Ykdla`H#FM*;z6 z;do&?JJU>*;J>^0!n z;KP_zNfflD1@4On&ZlsFgRD~^oTJ)nfD*%G0NM{;dwcs=)3S)j$eAw3>96>}qKO+S zLInijjI+|e<;-xLzyuf4hMWJj@tSDiNjLxQL}7xhOxBZNn8!vYB>X8Khw4KIjMAnT zLI^iYvbM3A1*#W4vrdz7_rDGg!cU$XH9^n@@la2<2ept6YrpS7-W)uFGtlio;pL1* z0dr8KQ>Z)&4H&1X3#H{Yg(Fl+7K_CxTtFowhQ*Pnz7W#_GZ2Bgkx@6|bekfsg2K&8b<28yA#P*9+H-s0um+})w^`)&%j);Ndiu(Qys@ph!oPPEXW!P~Hr z7Z%{K1VMN?0CNz~f?Svoj2$mO);O);MB2iiV`ZI^aTh8tn;G|z3G6Nz)>mPc1O8!k zNX1L|NYkUrM@S3aV8vmMRjPhP^t0<_9^_J0b_w)tj#J0V_Cd_!geciczMULQ zg}0z3l;B?z5aRwvz{D$8qEL7`0`h7Z&{V;Yyy-LoNC?pD7>MB|TjRL--V{k3g96RV z0cWT;@1(tcbn}l;e1u;egO;2{xOg9BlI8Za8OlUY>nT zXMwmENLy6m`2)3b)c~8VbDofOlJAiY*h&NM&}9Y&n6_*V?a|WEY!4aODQ!`*a9I9l rviaX%{xfs^zYi=hhyCAQv~Wfs>iT`kJ)Z|NKpT7NnfNtjI diff --git a/doc/tutorials/images_scripts/nio_csc.py b/doc/tutorials/images_scripts/nio_csc.py new file mode 100644 index 00000000..993293da --- /dev/null +++ b/doc/tutorials/images_scripts/nio_csc.py @@ -0,0 +1,201 @@ +from itertools import * +import numpy as np +import pytriqs.utility.mpi as mpi +from pytriqs.archive import * +from pytriqs.gf import * +import sys, pytriqs.version as triqs_version +from triqs_dft_tools.sumk_dft import * +from triqs_dft_tools.sumk_dft_tools import * +from pytriqs.operators.util.hamiltonians import * +from pytriqs.operators.util.U_matrix import * +from triqs_cthyb import * +import triqs_cthyb.version as cthyb_version +import triqs_dft_tools.version as dft_tools_version +from triqs_dft_tools.converters.vasp_converter import * + + +import warnings +warnings.filterwarnings("ignore", category=FutureWarning) + + +def dmft_cycle(): + filename = 'nio' + + Converter = VaspConverter(filename=filename) + Converter.convert_dft_input() + + SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False) + + beta = 5.0 + + Sigma = SK.block_structure.create_gf(beta=beta) + SK.put_Sigma([Sigma]) + G = SK.extract_G_loc() + SK.analyse_block_structure_from_gf(G, threshold = 1e-2) + for i_sh in range(len(SK.deg_shells)): + num_block_deg_orbs = len(SK.deg_shells[i_sh]) + mpi.report('found {0:d} blocks of degenerate orbitals in shell {1:d}'.format(num_block_deg_orbs, i_sh)) + for iblock in range(num_block_deg_orbs): + mpi.report('block {0:d} consists of orbitals:'.format(iblock)) + for keys in SK.deg_shells[i_sh][iblock].keys(): + mpi.report(' '+keys) + + # Setup CTQMC Solver + + n_orb = SK.corr_shells[0]['dim'] + spin_names = ['up','down'] + orb_names = [i for i in range(0,n_orb)] + + #gf_struct = set_operator_structure(spin_names, orb_names, orb_hyb) + gf_struct = SK.gf_struct_solver[0] + mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver) + mpi.report('GF struct sumk: %s'%SK.gf_struct_sumk) + mpi.report('GF struct solver: %s'%SK.gf_struct_solver) + + S = Solver(beta=beta, gf_struct=gf_struct) + + # Construct the Hamiltonian and save it in Hamiltonian_store.txt + H = Operator() + U = 8.0 + J = 1.0 + + + U_sph = U_matrix(l=2, U_int=U, J_hund=J) + U_cubic = transform_U_matrix(U_sph, spherical_to_cubic(l=2, convention='')) + Umat, Upmat = reduce_4index_to_2index(U_cubic) + + H = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat) + + # Print some information on the master node + mpi.report('Greens function structure is: %s '%gf_struct) + mpi.report('U Matrix set to:\n%s'%Umat) + mpi.report('Up Matrix set to:\n%s'%Upmat) + + # Parameters for the CTQMC Solver + p = {} + p["max_time"] = -1 + p["random_name"] = "" + p["random_seed"] = 123 * mpi.rank + 567 + p["length_cycle"] = 100 + p["n_warmup_cycles"] = 2000 + p["n_cycles"] = 20000 + p["fit_max_moment"] = 4 + p["fit_min_n"] = 30 + p["fit_max_n"] = 50 + p["perform_tail_fit"] = True + + # Double Counting: 0 FLL, 1 Held, 2 AMF + DC_type = 0 + DC_value = 59.0 + + # Prepare hdf file and and check for previous iterations + n_iterations = 1 + + iteration_offset = 0 + if mpi.is_master_node(): + ar = HDFArchive(filename+'.h5','a') + if not 'DMFT_results' in ar: ar.create_group('DMFT_results') + if not 'Iterations' in ar['DMFT_results']: ar['DMFT_results'].create_group('Iterations') + if not 'DMFT_input' in ar: ar.create_group('DMFT_input') + if not 'Iterations' in ar['DMFT_input']: ar['DMFT_input'].create_group('Iterations') + if not 'code_versions' in ar['DMFT_input']: ar['DMFT_input'].create_group('code_versio\ + ns') + ar['DMFT_input']['code_versions']["triqs_version"] = triqs_version.version + ar['DMFT_input']['code_versions']["triqs_git"] = triqs_version.git_hash + ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version + ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.cthyb_hash + ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version + ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.dft_tools_hash + if 'iteration_count' in ar['DMFT_results']: + iteration_offset = ar['DMFT_results']['iteration_count']+1 + S.Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_it'+str(iteration_offset-1)] + SK.dc_imp = ar['DMFT_results']['Iterations']['dc_imp'+str(iteration_offset-1)] + SK.dc_energ = ar['DMFT_results']['Iterations']['dc_energ'+str(iteration_offset-1)] + SK.chemical_potential = ar['DMFT_results']['Iterations']['chemical_potential'+str(iteration_offset-1)].real + ar['DMFT_input']["dmft_script_it"+str(iteration_offset)] = open(sys.argv[0]).read() + iteration_offset = mpi.bcast(iteration_offset) + S.Sigma_iw = mpi.bcast(S.Sigma_iw) + SK.dc_imp = mpi.bcast(SK.dc_imp) + SK.dc_energ = mpi.bcast(SK.dc_energ) + SK.chemical_potential = mpi.bcast(SK.chemical_potential) + + # Calc the first G0 + SK.symm_deg_gf(S.Sigma_iw,orb=0) + SK.put_Sigma(Sigma_imp = [S.Sigma_iw]) + SK.calc_mu(precision=0.01) + S.G_iw << SK.extract_G_loc()[0] + SK.symm_deg_gf(S.G_iw, orb=0) + + #Init the DC term and the self-energy if no previous iteration was found + if iteration_offset == 0: + dm = S.G_iw.density() + SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value) + S.Sigma_iw << SK.dc_imp[0]['up'][0,0] + + mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset)) + + + + # The infamous DMFT self consistency cycle + for it in range(iteration_offset, iteration_offset + n_iterations): + mpi.report('Doing iteration: %s'%it) + + # Get G0 + S.G0_iw << inverse(S.Sigma_iw + inverse(S.G_iw)) + # Solve the impurity problem + S.solve(h_int = H, **p) + if mpi.is_master_node(): + ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)] = p + ar['DMFT_results']['Iterations']['Gimp_it'+str(it)] = S.G_iw + ar['DMFT_results']['Iterations']['Gtau_it'+str(it)] = S.G_tau + ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)] = S.Sigma_iw + # Calculate double counting + dm = S.G_iw.density() + SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value) + # Get new G + SK.symm_deg_gf(S.Sigma_iw,orb=0) + SK.put_Sigma(Sigma_imp=[S.Sigma_iw]) + SK.calc_mu(precision=0.01) + S.G_iw << SK.extract_G_loc()[0] + + # print densities + for sig,gf in S.G_iw: + mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0])) + mpi.report('Total charge of Gloc : %.6f'%S.G_iw.total_density()) + + if mpi.is_master_node(): + ar['DMFT_results']['iteration_count'] = it + ar['DMFT_results']['Iterations']['Sigma_it'+str(it)] = S.Sigma_iw + ar['DMFT_results']['Iterations']['Gloc_it'+str(it)] = S.G_iw + ar['DMFT_results']['Iterations']['G0loc_it'+str(it)] = S.G0_iw + ar['DMFT_results']['Iterations']['dc_imp'+str(it)] = SK.dc_imp + ar['DMFT_results']['Iterations']['dc_energ'+str(it)] = SK.dc_energ + ar['DMFT_results']['Iterations']['chemical_potential'+str(it)] = SK.chemical_potential + + + + + if mpi.is_master_node(): + print 'calculating mu...' + SK.chemical_potential = SK.calc_mu( precision = 0.000001 ) + + if mpi.is_master_node(): + print 'calculating GAMMA' + SK.calc_density_correction(dm_type='vasp') + + if mpi.is_master_node(): + print 'calculating energy corrections' + + correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density() + + 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,use_dc_value=DC_value) + dc_energ = SK.dc_energ[0] + + if mpi.is_master_node(): + ar['DMFT_results']['Iterations']['corr_energy_it'+str(it)] = correnerg + ar['DMFT_results']['Iterations']['dc_energy_it'+str(it)] = dc_energ + + if mpi.is_master_node(): del ar + + return correnerg, dc_energ diff --git a/doc/tutorials/images_scripts/nio_csc.py.rst b/doc/tutorials/images_scripts/nio_csc.py.rst new file mode 100644 index 00000000..cd339bfe --- /dev/null +++ b/doc/tutorials/images_scripts/nio_csc.py.rst @@ -0,0 +1,9 @@ +.. _nio_csc.py: + +nio_csc.py +------------- + +Download :download:`nio_csc.py <./nio_csc.py>`. + +.. literalinclude:: nio_csc.py + :language: python diff --git a/doc/tutorials/images_scripts/plo.cfg b/doc/tutorials/images_scripts/plo.cfg index 91c72333..c370bf17 100644 --- a/doc/tutorials/images_scripts/plo.cfg +++ b/doc/tutorials/images_scripts/plo.cfg @@ -7,6 +7,7 @@ SHELLS = 1 2 EWINDOW = -9 2 NORMION = False NORMALIZE = True +BANDS = 2 10 [Shell 1] # Ni d shell LSHELL = 2 diff --git a/doc/tutorials/nio_csc.rst b/doc/tutorials/nio_csc.rst index 1192a3bc..23f60f43 100644 --- a/doc/tutorials/nio_csc.rst +++ b/doc/tutorials/nio_csc.rst @@ -3,7 +3,8 @@ DFT and projections ================================================== -We will perform charge self-consitent DFT+DMFT calcluations for the charge-transfer insulator NiO. We start from scratch and provide all necessary input files to do the calculations: First for doing a single-shot calculation and then for charge-selfconsistency. +We will perform DFT+DMFT calcluations for the charge-transfer insulator NiO. We start from scratch and provide all necessary input files to do the calculations: First for doing a single-shot calculation. +.. (and then for charge-selfconsistency). VASP setup ------------------------------- @@ -62,7 +63,7 @@ For sensible results run this script in parallel on at least 20 cores. As a quic Local lattice Green's function for all projected orbitals ---------------------- -We calculate the local lattice Green's function - now also for the uncorrelated orbitals, i.e., the O p states, for what we use the script :ref:`NiO_local_lattice_GF.py`. The result is saved in the h5 file as `G_latt_orb_it`, where `n_it>` is the number of the last DMFT iteration. +We calculate the local lattice Green's function - now also for the uncorrelated orbitals, i.e., the O p states, for what we use the script :ref:`NiO_local_lattice_GF.py`. The result is saved in the h5 file as `G_latt_orb_it`, where `` is the number of the last DMFT iteration. Spectral function on real axis: MaxEnt ---------------------- @@ -72,3 +73,27 @@ To compare to results from literature we make use of the `maxent triqs applicati .. image:: images_scripts/nio_Aw.png :width: 400 :align: center + +.. Charge self-consistent DMFT +.. ================================================== + + +.. In this part we will perform charge self-consistent DMFT calculations. To do so we have to adapt the VASP `INCAR` such that :program:`VASP` reads the updated charge density after each step. We add the lines + +.. ICHARG = 5 +.. NELM = 1000 +.. NELMIN = 1000 + +.. which makes VASP wait after each step of its iterative diagonalization until the file vasp.lock is created. It then reads the update of the charge density in the file `GAMMA`. It is terminated by an external script after a desired amount of steps, such that we deactivate all automatic stoping criterion by setting the number of steps to a very high number. + +.. We take the respective converged DFT and DMFT calculations from before as a starting point. I.e., we copy the `CHGCAR` and `nio.h5` together with the other :program:`VASP` input files and :file:`plo.cfg` in a new directory. We use a script called :program:`vasp_dmft` to invoke :program:`VASP` in the background and start the DMFT calculation together with :program:`plovasp` and the converter. This script assumes that the dmft sript contains a function `dmft_cycle()` and also the conversion from text files to the h5 file. Additionally we have to add a few lines to calculate the density correction and calculate the correlation energy. We adapt the script straightforardly and call it :ref:`nio_csc.py`. The most important new lines are + +.. SK.chemical_potential = SK.calc_mu( precision = 0.000001 ) +.. SK.calc_density_correction(dm_type='vasp') +.. correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density() + +.. where the chemical potential is determined to a greater precision than before, the correction to the dft density matrix is calculated and output to the file :file:`GAMMA`. The correlation energy is calculated via Migdal-Galitzki formula. We also slightly increase the tolerance for the detection of blocks since the DFT calculation now includes some QMC noise. + +.. We can start the whole machinery by excectuing + +.. vasp_dmft -n -i -p nio_csc.py