From 9e654897cc8187ce54fcb9bf58d4d14560d25322 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 27 Oct 2026 17:26:24 +0100 Subject: [PATCH] [wien2k] Rotate symmetry matrices into local solver frames --- .repo-hygiene-allow | 1 - python/triqs_dftkit/wien2k/_symmetry.py | 37 +++++++++++++ python/triqs_dftkit/wien2k/converter.py | 19 +++++-- test/python/wien2k/CMakeLists.txt | 6 +++ test/python/wien2k/wien2k_convert.ref.h5 | Bin 716936 -> 716936 bytes test/python/wien2k/wien2k_soc_convert.ref.h5 | Bin 717882 -> 722219 bytes test/python/wien2k/wien2k_symmetry_frames.py | 53 +++++++++++++++++++ 7 files changed, 112 insertions(+), 4 deletions(-) delete mode 100644 .repo-hygiene-allow create mode 100644 python/triqs_dftkit/wien2k/_symmetry.py create mode 100644 test/python/wien2k/wien2k_symmetry_frames.py diff --git a/.repo-hygiene-allow b/.repo-hygiene-allow deleted file mode 100644 index 86f9197..0000000 --- a/.repo-hygiene-allow +++ /dev/null @@ -1 +0,0 @@ -test/python/wien2k/CaOs2_pcubic.symqmc.gz Native dmftproj cubic-p regression output from the published CaOs2 inputs retained in this repository. diff --git a/python/triqs_dftkit/wien2k/_symmetry.py b/python/triqs_dftkit/wien2k/_symmetry.py new file mode 100644 index 0000000..5c140ef --- /dev/null +++ b/python/triqs_dftkit/wien2k/_symmetry.py @@ -0,0 +1,37 @@ +################################################################################ +# TRIQS: a Toolbox for Research in Interacting Quantum Systems +# (GPL-3.0-or-later) +################################################################################ + +"""Small, NumPy-only helpers for Wien2k symmetry conversion.""" + +import numpy + + +def rotate_to_local_frames(mat, perm, time_inv, orbits, rot_mat): + """Rotate global-frame symmetry matrices into each shell's solver frame. + + A symmetry maps source shell ``i`` to ``j``. In local frames its matrix is + ``R_j^dag M_ji R_i``. For an antiunitary operation Symmetry.symmetrize + conjugates the source quantity, so the source rotation is conjugated too. + """ + lookup = {} + for iorb, orbit in enumerate(orbits): + key = (orbit['atom'], orbit['sort'], orbit['l'], orbit['dim']) + if key in lookup: + raise ValueError('duplicate shell identity while mapping symmetry frames') + lookup[key] = iorb + + for isym, atom_perm in enumerate(perm): + targets = [] + for orbit in orbits: + key = (atom_perm[orbit['atom'] - 1], orbit['sort'], + orbit['l'], orbit['dim']) + if key not in lookup: + raise ValueError('symmetry operation maps to an unknown shell') + targets.append(lookup[key]) + for iorb, jorb in enumerate(targets): + source = rot_mat[iorb].conjugate() if time_inv[isym] else rot_mat[iorb] + mat[isym][iorb] = (rot_mat[jorb].conjugate().transpose() + @ mat[isym][iorb] @ source) + return mat diff --git a/python/triqs_dftkit/wien2k/converter.py b/python/triqs_dftkit/wien2k/converter.py index 37e25f3..d0ac4f6 100644 --- a/python/triqs_dftkit/wien2k/converter.py +++ b/python/triqs_dftkit/wien2k/converter.py @@ -25,6 +25,8 @@ from types import * import numpy + +from ._symmetry import rotate_to_local_frames from h5 import * from ..converter_tools import * import os.path @@ -364,7 +366,8 @@ def convert_dft_input(self): # Symmetries are used, so now convert symmetry information for # *correlated* orbitals: self.convert_symmetry_input(orbits=self.corr_shells, symm_file=self.symmcorr_file, - symm_subgrp=self.symmcorr_subgrp, SO=self.SO, SP=self.SP) + symm_subgrp=self.symmcorr_subgrp, SO=self.SO, SP=self.SP, + rot_mat=rot_mat) self.convert_misc_input() def convert_parproj_input(self): @@ -470,7 +473,8 @@ def convert_parproj_input(self): # Symmetries are used, so now convert symmetry information for *all* # orbitals: self.convert_symmetry_input(orbits=self.shells, symm_file=self.symmpar_file, - symm_subgrp=self.symmpar_subgrp, SO=self.SO, SP=self.SP) + symm_subgrp=self.symmpar_subgrp, SO=self.SO, SP=self.SP, + rot_mat=rot_mat_all) def convert_bands_input(self): """ @@ -810,7 +814,7 @@ def convert_transport_input(self): for it in things_to_save: ar[self.transp_subgrp][it] = locals()[it] - def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP): + def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP, rot_mat=None): """ Reads and stores symmetrisation data from symm_file, which can be is case.sympar or case.symqmc. @@ -828,6 +832,10 @@ def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP): Is spin-orbit coupling considered? SP : integer Is the system spin-polarised? + rot_mat : list of ndarray, optional + The local-frame rotation per orbit (rot_mat / rot_mat_all). When + given, the symmetry matrices are rotated into that local solver + frame. """ @@ -889,6 +897,11 @@ def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP): R.close() # Reading done! + # Symmetrization acts in each shell's local solver frame; antiunitary + # operations also conjugate the source-side rotation. + if rot_mat is not None: + rotate_to_local_frames(mat, perm, time_inv, orbits, rot_mat) + # Save it to the HDF: with HDFArchive(self.hdf_file, 'a') as ar: if not (symm_subgrp in ar): diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt index f94fffa..6c1b428 100644 --- a/test/python/wien2k/CMakeLists.txt +++ b/test/python/wien2k/CMakeLists.txt @@ -91,3 +91,9 @@ add_test(NAME Py_wien2k_outband_python WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) set_property(TEST Py_wien2k_outband_python APPEND PROPERTY ENVIRONMENT PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD}) + +add_test(NAME Py_wien2k_symmetry_frames + COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_symmetry_frames.py + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +set_property(TEST Py_wien2k_symmetry_frames APPEND PROPERTY ENVIRONMENT + PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD}) diff --git a/test/python/wien2k/wien2k_convert.ref.h5 b/test/python/wien2k/wien2k_convert.ref.h5 index f837d1785b41e90c6883e74731748c4a43426d36..870b04b5f614756ba9490bad1bfa767e9f1c80bf 100644 GIT binary patch delta 43251 zcmeHQd3;sH^}ity5<~%og6wEP6cEIX8|L*wRop9D_o|?EL8<~SRNO9AtXivvMP9HL zf(znO#Tp40=Xt@EDgrLJfwBgoqC8N7tohCIJ@ek1mrD{B+dqE!e3JX}X6DSfGv~~l z`OcX$_pj)7|B7y-OLOLY?sq!ooHs_CkmQc_vwY{*^Qx-Whdn9=e*Rs4!-SRWAl6e+EpUN#$3^$1p z36;&%KI1Y zv^Yv-`#}>(Tw^tH<{PB)_jt@LUT|ZF+4@;!UVw7^p1uCWBRiTDrk789f5$oD>3XZ) z{n$xF5kNXMv3a9(f>fD&wmx<9fj9K%8Oe~S=TyZp2pEDyWblejtKtKi#=uaPkRe^k z(0yZ;he({b%8gRc!2vRvqZGows^YS&g|R6E#BXukvjaNHLj89jnLA}Rz4-QKYJ#)N)sln52>8ph4j5vRy7Ry zLc41>4pCwG+`^$NYKm?KAST)S4>n zcFFu>-k;U75OmT0;yqCnpqF-l*cB{)$@1JUcn)=M{GK`F- z3+WB0l8Uxr6O~=SkG(*>YEp1q@T|^7K1Ow|Xoo)zH3lJmJQuIANiJ?kG zr|hmSxgVWdlkrjQI4ucMWqUYG7npH+{Gy9ml_hV_O*SN+vSeu8$hY5{QhtCbLN7*7 z;;A2udHSB75D~sY8~ITUSP&ljXH`)mzWlS?@nze_W#eXQ=N=%8V#pe;9bF6!}Sv%_K;_tu{l7K*!s4MMMA$W%W{>mgs(( zp+tnMU!zm%-d&=rxv5rb^xasZ*P_dO1+S#uO@AVxtiI3oxXaRH4s-Kfmg=+$GF~Xm z5^7%qPmu-E%!rIe+FMqHUwbn`TM@f5w06VJLi>w%+g(Fz&qiv0eN97k)<_*C)g-ia zva?Q4>s1>0h)`OR`j}~GH8FMum8<2DPg6(Mke^1UHd9bz$}gv&vN6+g=bu9+R@Ruf zy#|!56h~8fUExj5UW#`6=8!|<7eE}u`IuF{Zsyjnmj(ZasBaS ze`--bWoiHSWivkf)q01Qjcr$d7k20`WE~e6>o~rsMZsC|4()SKy|1WCeev5s{N|Mv zb!phA$akv5&!C}Zmg23Jtxi{e+@_(q>^KZcAsNgumJ2yuTj) z-Ved|lCi2OeA_#`>5rdpVypk@9evioa{jaU_Vhz(wVw1k=FM-s?Y(l=>|LD?42ADJ znX_O0gYR~c&M7^!yq^qF0srK&FJAcYxpH35apu&&g&2r4#Yfe$kOq%G^%}M4wwd2$ zY4ka7E&cJ?lN-fuxV%fFQ%>7%T9nFNoB#O7AAQ_t_R#e5zDU&@JIk^N{)~||${M{; zjY3V$WYMnbf?;!dH8$GZUyo9;sPbC{eZRU5o+Del_S(6fx1-M;Cu&+eA#)ftEktij zLp%N7c2J8E_e^a};F)vyE3lYt;%|T-VWMH2u8K_{NT@`+E zT>Fh~fg(`CPywe3ECv+gi+U$H315R=txn0^_!G;2YQKJ`{h^XM4FsXVq-7;+|E1@> zk=Lh^9pumZ+ePm;-()yiI-~37&m1tuMF#$E{ORNRF8ql(xfuj_p00qa3YKhQLTDYy zfa2FldZaxoCLZy}t)r4IY0n0&O(vf&)uF1YVWO^TL=~ZY4BEq}0?j~Nw#b~vZ%+oQ znrtodv--S_YGi~CT4es$3_8x|03+I0YmWPl(^8=jt&G}ij~_VZ`nAoCS~;F+P!bN&w!o&T?FPTswg-n#+O6A!^O`k6T{Yp{49LtFTPggV?Z&t!KtQb5u{9Prx zD_Fd4*5yuOr{H+m!0$<0!uL8{+_~R4C-%;38@1wZo29&5g8 zl$Us3iS-ourX`d9kyE%r4@A75PN&&1YvZ}sq|@NH^6|5GR@S@LMLMT+U1bUHCqI0A3DLLGld`jXr*!E;cZYn3o@Vg# zt_lt_zTx^LT@#5CD$Rl`*8%>Vlg_tu^BS^pk6u`kcM?C|)cH6%vfDW`t& z_7{_VN>AJ}qwBD;U$5ogj6MH@r&qlk^p^}dMSN-EZ+G5Sv4(Fd^Pf0;>+fJn>(XIn``>Tj{LkNB)bmS{b4qm_n#nO^oAPGD_EUM)r&qA6 z+)fWEkn2>Yw$>PA4>e5zh$5d9tBfDgbJ+5lNE-!Ba%Q(+uKyDJ=D+(9H<$s z{r0yenDo9C9nB22=RgSUS8MSNV(y@m7l)p6xx&bCJ#|{rgj%U-QbNDDPLr z*)UxJ`&cg>)PR!4K&a};RwIvOgP!wVc%|%`v9EK|7W$r^Wi5IiH|&@&NMAlG`8Tr> z<-~K4_9x`6tVxHm>2e*=QjVhuw=w@~DQ^-(d5n})+Hp=DQ|wpkuxLwr6alJ$)G%?WQP~u+D%GCQm+K>?=3q#auZQ zj(X0_0sqwB!}=>H842G%dH%+s3nq+rsUR(aECL}a;FM2D{fvj--{GF87%1k_svCW0 z{@==GuWH!kU|-HTJT`e%N}n$yy&YuKPKla$$(K2v+Y9TK|bvX$y}>H*Dij%!GJ1TsYvD%%Y=+HoI({Rz41`0H5%6=rSjgN?<=JW zkRkS|#^u=ZOai1-3`&&{5XkU1KiFxMYP4K-3`%K9;?U7Y?te)oA2#L2%KiVC7n54H ztk*gtoZ$~!UcaU=FUEaqHQA?P>9h-f`%v(@&X5>})udXDrU-VS1URF)dREY7rfQTH;{$; zR{W`G{OFC_?)_@!!1Wuu9>exCQv{{Ys(dd0{zePjS(<<;A0>p^0kUa?AGUaZHf-wo zT~F+2hWd41-Mz76VYlO{LJ)}b>!I@GXGrDj7qHV+KH2`aTrp(Lsm+^ZxRMhQ`6L%k z`mWnQTAX@WX&b{Zm9MsfRC#&;RU9cPKWfKso6J_~mS4#~P#|v}hTde@RnFqmUS{7*a;+-?4(Aajn%~M_api=1y?<;`t;!4unvlwnN*DgQ}W1SXWkUM$ieFqo2NM{-}%qq_{{noFV z48(w}pmb>%#0NJk7;O9%^!ud#YPZ8C0`Bb7gg!%vwD)2m6(ylHWyxbv6m3kOWzqNx zo%5gYX%-G$y<{#W2C0|P0(03$bQ#`|%!KfT+gA+w@%#q@q*Dvp<*_6$Js`~eV2Ip4O4$pz7BmF(`VHs>$>Joml(DS3>%GQ z^n(GzwlQkLx*+sszLusJa77%*%l1p>nk zfMF+qVViQmuqizT44VPNcKb_)oFWDcn-PYsmRwa{k_dUK^Rr=tEOGyT3|qBm+8;Wi zMPJ!)#GYci)U%*j zncE|O%E4okdH!^aE-M4z*&KMb8J^^-N5atrR|C&hf-NTiz}YV0-GH-AdBE8Tz}XCN zwodmQVqmokJex}+O77YK&rSf(j>&HYo=pND5?l>9+f0Fgvt6Vo0B5Uo0cTTrPr}*0 zz_Vk(vtzPifoHqGvsrC5Y7gLS>GGQG$5DGcTegYg>^_|BXC7s)GoHPtaJIjDcy48C{+w7^e#%di9e$dHFGx##Uwq#of@LziB9?wfZ6NN0!t=fWqR z9JpGg4?c|u2>-6nzQ3vS+ErZd@qfdfR%$bhAEFuHyVN3IWzv84cw$PUB|jdGV;Oyw zFF$S7r!vVs6<-|mS4|?nzUT24|2n;MY7g|+I0R6`u8?=Eb^hb$jsID)e9{i%6L7WO zTQ&6%3&Nt&F7Gv+uBHN*m^l4UjCIzGn(rxnl4dO8e1;r(eP_<)Tb2jQu~(^=*lWS5 z=v&cw%Vw1E?68++opkXyO6N--+jz{| zO`O-0=I)+B%7Mq#I>}H5!`hWf%x|sYcY{;#s%J$p6mRhf%vi zXJB)+n^q$=4X4%TbkVy-M)+(kTCERJGdrK6A@v_?f)Ev9#B&^5jV{ugt62+^h8KFfhNQ*eM$Ox^ppr(}PC3--pPfKahnSw?#gc`Z+ zAeEFsSC9-rAOWVl$=?Y;j!djcL<1N1g& z=Csmp3(ZTsj@ztEvDj~mJ5F!2vJCpd+adf!GeC{usR9B@29`}wv+_rTKZk^}XIoG; z4LKp8?7)J$IiPF~D7#QX*?!TTDJ&d~Nep0g-g_WiH5i;4-F1b!cJ+7HSau2nqawQx zW&1Mv2|2wwPOY$PQ7b{wq*%lVJ!=hR`zriOiUErN!HA}XB$h3ZWEW8OObcq(!bX^r zBRpL+lpTr7`4PZe2u2bXR?aL2j??)Fftld~%7y}JUw5g%1VsbNrgS~JiwqvHY!_HI zh6sa&6M=+vvA($i43|1!6*YGZB!%XTfS8{>j-k_Q&9 z+=bgzazNQJK-p5AC6o=MI6{kj>TvO_EtXx48p&gUMn;HaZJ=y_H?VA?zc-+4Uu{W2 zjUwKlshQ8TdhR}yO`MYucU$7i9iaotEciJ?jU+s6X-P)9BxXSo$DXCsi$&37O$HO8%yJC)NH2wWLsDAQ;-g`H$1oJ#0&Z_y*WTS zqjcF=Y321kkS{iZr~=Kv6xo9Py1tQ&+D|pR8H%3xQ%B`7`!`bbe7Vg*gy^Y5YUD(Q zk)kPp@}K+HmJ0@UXhEWl3ATz?)+j#K<|9&cUqGQb|l?@&$t!Fr< z3TOifu|GI{PfxqzsJU&;DSq`8fIw&m7twn%KO9F|MBAKrsxB)Ykn(h~lmm>W2E^4p zGxu9nhBN{m?Xyw^Ge_o+ZY_;B66_hlV7-INPSA3Q#?;eJ&xf2?ah~2mg_x-&x|G8; z2amX{{=^dkq>~uhC&k(D&5vJv5*tf?2#nHOsPbxucEr3Z3zQ#%G~NsW-R((YAgh(_SwrXF zBgLuyG`fz~(<)MNe)F#bW_R0pn9<-1w{L&_@4b#Oii7fzibLglJf@8|^p)e>^!4l0 zj(YIdMo~tVEI#$j1+7d)^v1ODSyNmh$`8fCVCpzjexx5!0coBcw`lNFXMObMyLrdk ziZkn|7tg)#j;4yy(xntOF&Wl0>-4fTOSMEA0i!rnLA+>eeoJ2s_QKJIuhZLcJYUP8 zHQT)9K0HM4ph5=CihbFj`;R^5``amM!w93zONhfLE(DdPTAe9yt8i$4E()uQ2W)klU3yoVA z7|nbhEHu}$;EW28g{PM{ex;fEV4=x|l@{8*LJ7<;t^}jHz%Rd`g;wlKj}zS&{|yVR zvU9a2plzYWdmUBK%D2oV)j`6XGCf<|vx9)^sq0%Gv1azS+Y_)o;+-Uvv*6Rlhn8CHF8~>|xmGtA{-dhCK}0qoAR} zfUnrYpgj&FP_ze<4!$-&b%;SuJoRy!iP*)GAD+lUIs3?YUjHd!MrNb*g zU1>ZGYuF5wRf_RxO90fB19jE(dO&Y4&`NwW&C_+ z0tfmkIFez&-cxU;LdCR4%cN+~a92}gK`LQ;snCn9#Gbl&6DmLzP~<`SUT$5d&Qc`3 zYWCFCR=Cbz!shCuVToxtxkZaR+Eed`46%vW6A%Rjfo+&eOSpZ~(qbTt&p9af!Lrxh z$`^rJi2hMOKAH3~Pli%o`6T9)9^Va3{ z%_9L0Pw9N|^XE;C2S`s8Up8;Yf=kJ7Y#w-XoV{(F94 z7zhz3L1IWUPC|e<2^7Wva-nxB#rcI;T~D^?R;KcP!vu~9al}brh?5}4FBvDnmHJd0 zjoPg+338yt#-I^s&0$S(6827*1Q~54Yx~;PsA}UR2rfzB_aIDyuWtg-hnJ%?%EZl~ zmXStedc>J({Jx}j3MmAmFbD{fKvR?O`x>C)Bm@*k0dW$fkfp06X_WE%k{Xr9Kp0L&QsPDV_r2B`7r%@e)kF z9+$r>!zD0;OVER!?-aSzR0!VVl3G!09a1X<-$Bz07khp3j)L=0Mc%FsilqmPko0{S z9f(feiH&|4yax=#quf@A%dg7?tS>RGXM7=rf%2;RdeHV=aLkc@eo6Fg^jr`9T;TuI4v$^0I^>j1NN?A9)1KmvMXyZ>xj`AK#t!T$AC^2s}cTpXklu zjh8ruWck(A;uDtM*Zo@MpnRMkeWkA=n;@P~1+y6N1V7?uuUmXZf0k)+%J!n%3X2Tx z9cUqvW(vNY+JCd7VRb?@EOHM?mDHcMA!eWnn@bsNg^FU}(9CF9=8cFB1r?)9i2kCox`nJu))ahGG-ubeLM~p9>yU{4rWE~BQRp3ZG3A}67aiHdL;a=T*-p) zyP|`^9DY}1?|JaMqI3qos|5V6l)ehLD;I88F}PhRHxl3JWf<_J@Kz(fA+;qv>*ja0 zOSdau(cGtgSE2^8j&Q=?go-V2sg5?+?{pJVPaosyotIeDJu@m*mn+_koqX|M^f8{^ zd9mj*qEeaio}}-K3Sufu!=ZVx4>B0ICy@W{gBjcmFSMcyRa2z>X^Uu4i`9B372<=k zbg<&=*zt5`sB%pbc0MYP8Qg1Ip^x$C)wM{)cmr{CoG1cRK~1;4P=%UqDxQwMV5?;; z@>Zv)DCLNw!$k-X_T5`MT|1tRt}yAzm`n0hjoLaFp?n`k2oRc+-bXkcpF)061|lHA zdwD9xj;BL$bj0#$(N_;w1rjj%E$-p=@y%dCL0r#~%#$yPuV-T+DVPK>d z1K|Y0kUV}ZkK!(W;cavL`uT-7Y<}<^6U9m&gu3;L*(`6NI-kvuf$Rux%EjKWyuRaH zq;pC~3n}07s#|B{{Pn9;K1K}ad~al)_01A>ePc*Val>Apru3N%%&{vjti|mk@-Dzu zB+>Q-B=JQ}3yS^cWRnQ*qdsa^_cnHW@132dPN25HuFB|V+zK7Pk|Uj(BK&D;vxBa? zWdfyBc}zOA%Uu+JK9z@BiCya1ngNfOF}IHkjjcMYfV17jaO7o2>#tLeACn+PzH|Ms$XrcLD;XXCRvM_&t5 zB~E|#X-c2dI;R%#{rgbQg!UCP@utlSCi>KNF9tXy&zq{?+Nr3BN{mL$=4)_BUI*7M zH%(cu>al^BlwA=Z9r#|kn&ORvgK){DL--l=_tSA8MIA57GZ1AnS{B@|PPY8ajVxX>Zvj6j`Fm}AzJ;ifgu#}4o0 zB=pc$$#bN`wFvku9^bUvZ+y@+_v*ji*}0N!@11#|=!S3KK7VeHhbWzYiUUPW5bya* zEqlCD(Dy|JdhA@ZNgt393@FaT!b9cDHp8=|r;ZC-Pt#Xcc=bDRxHTkko`OG8(a5Rb zudQj(#%1I+ZsK;4gq@Gd(;(6}Vp>ECI1UJy|A|0v06#Td6G!=NT;^Ovfuw0Mt4X7- zH_-%>p|pmt2^z+1T!sKn0eJqn;O>nsMJIn2_zZ>x4?t(&6TnFsC_1^QsN~KD8w1PE zC)pqd65k3&(n4^B{V{8m;53n9cr?YqEn&MkEgAwG6(oufwb5e0yHHHnQcLPNquukG zP7k&h-{c>(-)ik!w@9z92Y*Qxy>R_jQ^m!b`6%mjYxkuc4q5T{TeJ)J#hckdINd7W zND@$QCLsV_f%bSExS2nuG7M!UN_`%Q+Ah#t?rSRQZyis!rmK-WBc3Hg8*wSwrCN)O z?MGEgn@^r@HPUUX5pE5t$+`?TS+}3ri&+?#rXFJ8oE4jVsi(1;E;pnCwi;cakImBd zUsel{acSRLr0FHYl_FNS_*rkq;%8}U^h)l%8}w1MffZ~{rW%?30&QGY@VQ7&)1u|4 zoUnnk$Sh_J1FHTiy1OlP%Jm3D^gRp*0wvD0>RLLXoDA!_MzhL$hhq)sx73R?^vz1_ zx=rU;qj9X!0<6(uSfeS6NUYJ!ET}>$h&342Xa*C)8cjzI7Hc@H(Nwx9Z_Yo3v@#uF-rU2{ap_BBQ6j#2QW3@0dke#I75u zk-oAsb(0juI>j0-u$E9tMucS+VvPnvrkr63CS^dKZWZwvb$CXs(P$Beme7a9+P7B^zXZ8)V+hO?5osYRYr0*b-Q= z`P^c_lV>BUeA4!$H+1v~XB@U{jSiC9u3_R0IX$bbyg4_lSyif|)_;(fcvDBUPZ{-3 z9_|+}>PP9ZOLbIxHdxPn;>ASyKl-HsD$wyMSLG$PWik-Te`}*uyEQdQ*@{xCn=oa3 zrGXD4^HLn;PixBu%LYE^Vol^@x>V|~)G1#cYZ_Vo+++$w5m*FdrKX!ZsY)C8=+kYA z|3Q=A?4Urx9%ZB_uJ+8jFOHN-1Qg!hp zi*>A1J_#8#1IZ8sC{i&6RgIBS2|`VuZCnNISWPnn1ANE;A9BEdM$1g#1s`)5IEsvF zZxQf-A9CP_P*cNznQ}VFHxWx)tgUeX)vS&AQWq_Sh`QmG9C0U#g_3@ShfK!S)!CUkIl0U=osld%#2K?N2> zMZ&JO*#rc|eq7LA#~37xiWtys{)W42n1gVqJ;R`=Q`C)!WiJdWXxc`!N-_la#4T_L zjetfQ{Y(JxiNGal1P|a7rCZ<6(iR~Moi?BISJ z&8eB!r~nO4_EoYXUF?jn40OKf+>Y-f_oqHqc?anx zbw=||rr;TCbDRF5``y85Is^MiUw6ix9paU$hGCF&oa27@`O8;V|HL&!vPNB!QgW|g zAU@)-`Um+Xr- z@`%jPWG2Kc8B!#k!!s4f94h&nRxEZf^yVR5UU+G!9}#^Gk++DDOn-BnVE~-Zfnqq$ z^j}@t=YkhX`KuH5FZpy(g-S14OMq80EEqna{wTMCB`5v$tj|B6X-0eRpgX5L_2wv7 zr;BTlIy`!Gr%(NVs|;xbyq9zZ?4AdgG`?(@u0UKqVHl7@gSf2(V_h=Dbtw$jr7$@% zi0e{ZqfH1eM_QtW78w{QN{w(|3JAnGu1n##E=9d@8Uvoxp}K03FRt4}7wC?Vf2+GM zMSc%8^5p0ApFOxPMa#dI@?GASvX8Dyf&SHeUy4>wnr719UVg6Ox|D#nd2Q&x9VRzB zMm>ivNQSsC1wKh;SI9sIbYlaZg}`+wRNk|p10iJcx|D!cI}$n=_)dPOvMpl@j&Gw6 zm7Hmuc0*j3!f;)RmThzSrl6&-ONrsS6r)tQF2zN90@tN*T$iFmRkjnZOW_J$6W66A zD2k#f4t1lOe`Bq}YgOEKDl>rxVQT?!|Ta9xUv>r!I4E=6GvX#_mFk3}5{ z@0ypubtzzgwp#UJcxhk3z~eCR7^Wod!@{Ghy2yefPQ9qBykOyRSa^X3Un&?ZJkxSg z=cR=y8bG!|DF2r~Fz~Xku-*F>o*cNN3_PKK``p6QFI>HWC-r|{TX;mjIs?xa79Q-G z>C9dl!NN0r3B!!R!Xu&eDD){TJPr$w!@{HTVKplwUs`y?0M{)wCqRplg~wpwsdQnW z8R3M5$6?_ahG5|ZWC6m!Q|W2t!C++ANg@=$Fd0<mNcGNLuN8vrxhpdlz8vQ?@TbAG--KkPwp=nROEI0V6AfA1Hv~_g5PW&n%P(fdIgu zE_f3N1rZ=MN)LUieE)R=P1lj_K{b>tuNe|CV~hPX*NB0e89Wy_JQt!Sr;y%=p1ZX> zBnn5ucD3JhinU{01WYf4K1fk&X+2Hti}n92bm?X#XF7 zOaRXX2G0df2OLHBNY4c~!vBXDfG{nGmQ3exT;OnAVA8}Bjtem@E5^B8-+AndOOavX zMeYt+1VVT&FgPxN0RvSW>jgCn-I^(9dXp><8ZyZLJ8c8?00B&>~ zUUH;H#{UOibd1u|v`Bi5Nd6>DlZ0C1!BP4g&L8Zsb^g#BmcIUHLKA+bXCy$PuSTn@2Y_DNBAo2 ziK_yKs{)=q>v&IlYN>z5u-DOfQXRv+hND^>|5oB{(Oa*!ED=yG|aClEM@iz;N`AS%z&ueb6? z{_Mmo29zg(h!#XYWILZruIEEc@I9EptN+4R0*{nyjAeW*EAzBdTz~6npI!#Y;96h+7mD1CA4# z10vQvkeZcgbGda(qG?IoqUc(;ETSu^1?QRrBDh6Sr>72L=*R&QeNb;bdSQTD6zP(~ z#7H$4J)#*NJ#40ROQNwLaf_mB-Lh!P;}%6lU*4iGR_}-HjA9`Zh(=Q}qvDm+U`igSF z{nTJQmALn>KJ(xG6CFi&crf1gTZgWB_?7ZA_C|IAZQ91zps26BQMiu6t&4eACcP@4xIH-uLRp$-WIY7XD>_#ow4APUWi% zXT!uNf5<{L-@g;g>Fmk9fjV z;!7n~m{*<6XwmjvUR{`fOwSbf+0dJ#TsM?o+0(W%F@9U}Dz{@WzM3Uh$&8&$trjm` zijDtn@)pRpsGb<`vs?`p8gi+`eMXFmXAgiOo=OD3>{s|Iq3{^LHw&QTD_>%$(>n#^ z5TR7)=X&H&wfB^aQwi>oqLKKoh7`wi1bT;iIEX7X|G2XAk}V@OC}y5yu#%mnTD>ENev-?It)ke3 zw@%r&_OwYWlTp;vkQDS-Z{LWjhVu~F^w|(Zw{fi_ zY`jdLdQ=N!;{#^g$LqleUi%ESuHp++Zj>mc&2isK(xh9>1NplK+kSMtXKGufc;i#Q zoywQ*Y%6yRs6IC@;PpLyW0tk=Ou2^!J=}JC z%i2n!Ee-N(>Tn?L$^FZE^AG>!eDBZjDj26p9RQ?7bU=oYY9gFe(&Eqk()MLT`HxJr zO?N=F#&4#Eew{p02!y*CF3%ANNf8WORuQ(Zr3bRwfR3^&R2F4mZRfCI)%(dj@T5;g zhBckUFn;pAOH~-KK5#easzXo>Use7}SZ39tcT7bo42No!69!<0a1}e?C!9gKwru_D zRXjJl+sA`uzwnL*fCKcdT?ha}rf0G5yuI2x`&%*H(K3-4Ti7DY*%~X~xPr0sfPx8p zu}fl07aE?h&8%YPNpA}X?yOv0mYlaypKoomEt@e zuJDZ+pd=^~tYnN$>&6}IH#vo$x9v?tk|t;=qptKdmEds^*@gtbPsGVcXg^67%MhQr zAYSg!qU91=DxWjmw$~|fEhSy8_bET3~S)4F|+ zfh_K4)URr|E-Ms3(D0@H1R0tmJG!TAbHmy6sqfA|K!TwIq-uOlNo<~)qX5z(je0zf=LpZ zCREizAih8wS1|*)uT@_o_DY=*!9Wfbe{JfO?_a<9E!&06Kh^IYGaJItV*)AEBg~Ml;KzuM2A(VYY>|nl0^+u-|>~ zLl+?u6?n~|haeiMFXne&Br58cg7uwPOc2Nk2C_;y)F z@<<=7n?W5+)!{x$AGj^+YwC zM7_b~gz-4!?v>*oXmQT{>sw+j=%xSGf74GeHp zYYcPs+}dmQJUHAlpfx05-@y{42GbTrG@|LZ_f+1c?+wnCVtmw06iawTz|<;GutAzy z-fX6pQ9{F)e!nuCa-HXNJU)V6MD&<-)lC-JR+d?9k&T_)=N|f8;=cUk%N>IkT8u9Z zq(baDSQ~fy%g1En&Q2GKvBJ5C<&CdZAxHf&T z|6Ll)Xaj7ZjfQZ&+>9~OMchl(7~(NVRDGT0n{Tl1xBOKZ;G`jq-+A`Hn?RH5!cP~C zx*(^BEF#_kXn@HTfp)W`@38o5OTRS?IX0Un@K~Fn`#I2MEpP95 zVEa8UKbx5=Jo4Ow`Nv;eK0-ZNHr{YbWpYrC9yV9ExFx)*EyDs@m=ZLcJiP7&(?7B=Py43e1em4`-qih$8m?RIdp>AlJd80L|cV?4=um8rj zop=s*RikDU4@+Wb9+hM{e#e})&tiLc{l5(#Om#lfXF0*cQN)Ndn<}&YBKj>#K&wL% za*u6Mxp#=c2<1lHMG2+O##+$Us-4UmBX7MxURM^pj}Ld+*HTNywkA^DD5BR|KjI#3 zPxP!Q+dc$_>cLjJ2DUD-5+!Vn2yCUzOgLjJyw-}bWvwj3!EYQ5ynszgW&*ZA9)W8c z_mKztC~8m}nTHL^J5mM!S3^R`DsGzSS-=mX3V@Z(I@#6zkQF^J$oJW`SuKPuLM@$tFp=wy-fo^*E}YCv6o& zHY*MWFhXPnwP|?k4PfmF&#{>B1|_Kv!kY+vQu-hW2m>_%#hN|XbsV<@vM!hXy|zM1Uit()1vfpNU)Db$iS>^82h2k8ypjrENgpiQq&{Gt3+SYkcs!M4#~e~0cvj5J zE=NsUJ!Y~$o=!#S15c4z4vjLXKbZECBao1=SsxJV3TZ7?+O>tPqNnkJ4FmbqwTWH* zLOEKE#SnKN3t^FmQ@c4|b38coUCx738v3eyO#$V>QS;2}Nt>ZTs%KJ=-a(>U*eJTI z#whw6%+S;l$?BZ$(D`M81Wlw&gp3v;{sx;`l?ReI$s&-}Rpj!5$U6tT>wB{9p5cqnb@2hs!*%-rkacze zE|{QfU|$~hT%Q&_@5R%S8YS{NMI!HQ8>)gD*DDZJ+7vc05S6H=)_AiRzX|#>7F*K3 zKYk!_sBhO|nsgc!w`lCdzAyjgkG}J4yT$kC{^py>*zb3b03N!JH%xTsIbZd$SVIWo z=C{OG19J>Wdem2`#6%}S)IBy{hlLZEqgDOAS&0&AE*8T^KiqFP{rMk%X7ZOU=D0pB zyWfkaCAGrGfAN8iwHL7g2i0DSK~P6MyVpRi`LLL?!6cpvY$qJ$I`*lB%D^}iAfWEf z7p*k~0nK6}&nEALK#YN#LKJsdxdbG!IL3E&uEuPRT41(Mc~t@Af?J+YlOYgagj0R5 zL9XHc4o1qlqi>WQ&G-W8-QvD((Jr7VHg3m(BH`-{ij*%jS7pqAA{mvHZ&vsbPFJDo zZ_hhQ8+2Gl*^<7_tG+(A!QrhQl(T^5s%vv5VmyLKlp(XmC1bCyJVz63IZkI4kYJMu z1>5B!U24MHSnbKTsQMVOgMp06CUl=j%yb}|z(p0Ac*4v@DXqk#l~`o1z#1#J&2`&F zO`I7+eVIA4sWvc>;f3T`&ceur=1jltpigQp!>($k`LAEbqX+s3j7= zr!0};cc(s@Jd!{AL#-?V8qnkxVbj`i@aE(M1Mq?+oFaJO^fKqmR9*JU+g7b-p*fuuS+G>A{=0aSEjqtd z(%%{dpOEz3?$_Ej<)42sF>Rx;Rx=mHo86EO2JD0<7Q_fBV@)r@Iz+PAMnTA#CklhE zq!py5_J`J}&q)S~EsDZO*Br$~)V?IK#ysR0wo#tk>W0*nT5H^OeYHlt3yF5(NHG*K zof%>vaFoT6x9pshl_plb9qGKJo_C}lG0b~MT~M7^Y?@=<7K?u?@FTS`Bz|}A_d1T` z?>wa)3p1C)T(}%91+ldRLzW0^me8BHE@-iqsvLx7AIQaG@14nbPuLNQQ-Gvq4^VRg zZ5elJ;(|^2Zwt=f>f@WVHY@*hw z1Rt2o2nMnds$q>BuuLIN-GgVU@2J{A-2k~*?7eKZb1fD>w&K+m|D!W!6A59k?TQkB z2kU7{ke)Q?=6ByHmCI4z=BOg-+vSaHL}o2YZ9q zW(5PK2y0@20{}GWXmh5XZxv{WAW~_Z2@4T0 zTV9iV&>g-zzHzLq4x(MFc6L2UM-nKh@3P-@$v*GQJ0;#5mgn???kj7O54&%yNgj1C z>q$P~UUGMQz*52r4Bw_7HN~YV*;N$kAk4kyAWU;!o_kS`qGkx))sx&YR=y~ItT%b@ zg;>7OjENvePP&N+z1?gVifpw?=vo5JTtRkmgLVL0RgsgQK?4kexq@II*IYX=@)l01 z$iaYcR8hh`uE-K{k}t{*c$pHJmLxmCWTYJ^pc!?onwR0TnKC$SNmGLz@P>GN1eD%} z%w`~S2aWTrVG$jQ4rirUv_v_+R04I(E7g{0s~B>CRe@2waU{ejggRp*?4Y5v2Q+bJ ze1?0@TQhIodr#T1T`q$bW(Wkp5JZp>4AhPJ;56-ZHBTza0d!HArv%1qfieeOY7b)7 z72h7_DfBkmj!i0^rj7V$U5oLE^oA`5Ta;g^WGGdkmYMl{pNDyVN}m-9*GHKUZ}+q4 zP;@vDfFrjckSK}P7O22{q(ocAkiK@S;yJIV}UUN{0C$J*&A zFvKEgFA4_g#(Xf|yRPPGL^UyguEl$c6{!V%q(lde@TH*-i4^C;b-zzds=oZ8?Lvx0%T)gWm4>$Lpy z4SPODhsza4P57X?BGDJzlY1u+0qP*+uXIgn4_wBs11-qf};~VE752mnyM-{(MNr zk`)PKsl8#5MMw;532KW>A4aqHpA(z%KR=w@8QxC~`W%fni$6fpDl^k}LXZN+$~Z z7oQJVYO(Yk7QdXs4NG63p-FIbLjM^=-3+mpEP$~Z0|E^8f^faH%HlNOnQXtEj$ljW zMTCNZT!el|A8bdj$+huOP9UQ?wB}>TQ5f@&R>WX4 z#1IcQWDns5SZBQWDd@byYnB~;fZxo{^dJBO)-!lcFc5|HpmIg~0{#B<5&tt#3BN}T z0-7`Uy*e}sK)bQ3qhD1jM+5J`;J0bwhg1o6A?Xj)5FHbKd(M%W zJ9KF4BuduFE?S9{sE}~e7_!zmtuWzt6H;#l1EGCj;(0l@iH&6op1G56oZ)=79XiyD z`KHg@ontFB(4pvXK(bzm+ALM%#GSxA%w$NSLc&dBNMhbbB8=fi8RF;bUM?BELZild z$h%)pMa8;U@2Q^iC%j2jb>(Jk{cdFEH!uu3n*HBJ1W&?lVFU&% zQ3`sx+B+E1>Q#IRXZv6M!4>A@8WeF4*NHU9A{n~WF1f`U$yGHWhTPj;S2MK>Jzf6K aURSfT=FivPP~UuH8a`b7F#(Z+sH6?c17VNBfA#uD9c!4 zEQzv9i7^IQ!&oz6@Vm!y&Uu~loPN*mJg0Mhe>{KqWA3@fea+{4eXs9zy{~KX`sANn z+Y&AZljr91kq?lEWXVHsh(m%Sp3A~^Par9l7of9zWT~tbJmkdpyxYhcsfY}6tkvWaL-DEj*i!l(JYi%7!QGI zO#;T0MC0Q+%j4hQ9YVh;9R+>|J*XQbfMcZHA?Rroz{Bf^xd-17obcuh!Zl4jvaNPT zd#qV7K!RHnA2Rl0+Xj@es?j$-o=FnL78 z6h5ijx?kuJEh!f$dI!%V*3tR6Hf~1-%z{R972Ai+F6c#d0QFih&+Gz6FPJaSb~tWP zu3?f$yp5+D?NTg$r;%r^X1@5sk+$ei&W7kIBsyeDKz? zQIwq?(0*=6GRIH#zAX5u9*XR`2_Xms)%lUVl6Ttgt|-9F(bY!h=ZQqRNt9B}n$cm$ z=(VC7Wasmn@K=D?LRZKa~OF3aI zWRvb&+IhB9B+gBtIUnr&(kd0cU!PhxHb||Os>D07lJbSQ;u2#d$Ui)V=9Ak)K)K>M z)_q>K4p8Od;_`sF81I-Lnx7vY5)I@{KW#!m#}WmPH$7ct)ndIX=jC%N*{v-8t%P-E zevAS?Suf0f|8c#S%~n?#1-t-SD5}(sbl3WLXpf31_tsQ44_7|7PO=B&tQD6JX41-m zS77NI&-}n@kW1AmCF?*!&J5a-_96H;uXNUVhg$;^4Edq2gO4j{PK7}+xeXoThD3AB z44PaLkwSd}ErBfVUtJQ^07d>_jJQJ~=Rw=KA!Z!^Oe_NZ({%ZJu#`|Mttc;k6x+#4 z*z?M4l)jHNY&FVLrJmEOX54y}$#2HIRS@|RmleVgGN08PwPqG~xA`~=OA{q@#>JAQ z*6EX!sS#@~?!QoSi@{a&@#HMIkQdQOR;U3LBygVhwg!3}1zh0@Ix)_eo1Z6XM5^^| z%D?4C(`0^mmmB*?s+LMfR}~?HTR~zAb3+&ZL*mdA^7Pk>HdZ3a9U+zKimX**O;!;> zJt{NBomg$=l>nTwH*M8XL)>LErcW-}f|=BaeMn^Od&f#bSaL22bCK>Ie&bOx9*d5d zf$#B7%BRHBpd;Li24!{pp+WSMZGul24}VlM77yh-m|uvn${F5RFp2o`0f!k}E&0a= z4a?k1uo%l}g}LAgTD2aj8?l%{DZ^TVg$tLs{6|i7s~+PyvgXz72s>N2nk+ua%TLzp zHTq z5wCQ{*FF0?ntXdC)OHinov?F~&A0<}j!`grb#pfN_2P*0UPxMsbkj}Mb z=m^VQSD2UaTk5%;(@NjwInK{hDTI(l2enhLsDv2*%qNkWGe7d0|A<~WPqhsk(*_cR zo%35Mq(Z?$hZcT^{INY@z^-z*<+Gs^qCFX%FNN?2IB_tU`P?_J9QwJ9Zp3OaqInC~ zO6qZ!ow^fD-5^((9)yu>RCg2t-_j~2XlxO;X5P!2DH`2;p8IO=8=>%#4!_r@X#CLy z=_6X4fA(l-@7m~~IL74NZi;!|WO6s9pl__Phx5v=UpPh-`!%_rcRjT^=|0}HTg4@F zF%;L2X^@TT5n_`1tm6d94b&m+N_;vyEP*juG$v7{ZZ}np<&w=hj3OgTkX2+$C72uZ zBuCr_>M}4CSu+4;02nNUn6p+tZ-hfaJ1@q(9uO+4WguayZ)K`m>E~b(_t# z^SaFH0kP=44)17@Zf(**peGfz^10PwT?v|WfR81ax92{4&1q;Me)N=2EG7eU>Z>L_ zADBOUA8t=5+g3t8$|{;Tbvm_YMM`3U;$EMT=Ox|PXrgrw%6)ANW3sdrQH?)Y$F|`Y z463hVXY<7{+B>*!KG=x{FDtrOckmr%h)mwSFB>|6Rth=_ex$v_kF(J_4>8ima6U&_ z?G2b}{K$x3bCA6EKKca=ct02eut&?1r>ccg!Zg1QE<6xYZT2REOfN?2kmK4Byp&Qc zSUWdViK8c_QWy4{fK%yoXEB`pR0F7HDa3diI3-l7_1o=V{DO0pmBF`>YJi;`C|(j&Z{eWH>O-LK2IuS-D=@i66S=j zfV^^^6s!(#B}bSEXhf@K+1Y&^^1UtlvC(GisChGzqvZwaapTvQm0&QPuO(H+zEdu|$0B-{UwmMX3NczOdCE~(nOgi#&;T35c$zeY##$&|c6is< z&n;9~H*wus1AMbVP-i+NY6g7M0^8w?sVc~6P*J9CD;|eswZfI)1uuQ?TblKYeQuit zKe4h<(Ph(p@ZaL9#LYTj@{$V=m}QS7V_hV|Jg6*H{_$x;c?RShz|*XC#!<27W)*vEIp z%_TJ!4456E&sU^oN!ajnyf2@-6DuF>DS~zz#)>4ixMp=b9O?S%8&KO{tJXLH=zrbT z3Uf=FJ!&o5Ta7}Abm%Lv0_M%l;vMAC6v@xnQ1I3TFVZXAKo`E)~cR95|DyN5sEtDpw~@596Dl0QVP1gaK31di=|GSgDC~ID1(BQB82vW z8a5WpB=?JnriQgJeS1_qJlfyRr;kB6>&(rS;+^S7#-$<^$OqQJ1jv&j@I92UbFfHm z2#ph1^&Mbf3iwMLPipTBtVk2r&2caOs z;T9wC)R=yvaLQ%f14(1A%UZ2`1Pue~wni`k;5(-B(JrW-t#eCGv>G4Dy*_J9jLSr- z-0_z*=o{&S%Ad_2e!PanQbq_cIS8bbqo(HN->TWLsOHmOypSej9(R&$*aWB_aIbxA zl`B)h3}~L2j8-Gg49M(Ua>HSttEGJno?_OOa|%(4d$amAUUN%3et}n8^=u8ygOiCw zDP}XmEe&lTnbq%5`97}bYYQK^wY`tl;DLK*PVZ#C>b9s=jN_%O zA^+T)YdW#|W~26SH_#JHszfhlauU5mk>D}>);Gd*3e6RWdhwi-=vf_%21_Z@a=CJQ zxm+3Unj5H|CqYh&sMp;H%7y9N%_k($bt~o`JsNLgQ)|=pe^XE?=Mn84HP%$=*5H%8 zR31W<}xG9x$Iwwe69Q7=VM&>!J7F6KR(F3WNbnKvtX7J>ke%| z=+S%k)V_JF>!?qsV84Z!X)X>74XZ>f!DC@xDg1%Ae=DC-kcD#cSMiHXnXP7E-gv2F zf8J*8ejlT0S-Q?~iS;N2DunCI>hlsO-`PH2C|H0P*iymrC6r^J@sfed0Ag>vuMvEG6yAki8JmR7`ht4kV-C6_;2g|{nWJa%@W zIwVK2e3YF=njs4uoqvK z%NN5u7##TR3jfUoDK}V=E-83DAB^Pk3pVynE!hZ|Ey%c3(*OO%Deq!+8Zb+lt$+=3 zLv$ z*lzOSPA3h^R5jTaT(-^_)41Ju`z8;k`QZRr4f(|EnXx^r{MX4Dq_(+GvChx$-Vb=- zH9Ou|y{ONRPY{;O7Pe`ft%-hJ$v$HFkFdf1b_@i|g4llQkAc-p!1L z_x3~)?hdrvJGA}F%!gvUk@LDpW;BWB&2-AxL3=q69O4@MRyD0|Fj>%%)jX%+zHxfc zpm7k@eEZ=?QQa2iS=&3EZMz8Wf#B61?`H7zL3m8T6mcs9;=XY+&3n)wG_*J9^^7C4 zT4~xFe2p59AC;)6*>sOKXl;9kh1s5QEs}hdh%6@C)4|1Ga8L?Q8%N!&$Zh(tZ6bkw zzIXTUR+T$y?siO|BtQ_*w59uZb}<0~WW3q^j02YhwfDU$dt-8`?pa;&`xAytZ_Lc@ zoz-Jb!py=T97F7E)tmQLaR$}fQ@>cM>=lOgJ!Tyn@ZK((iruPHXXsedkv#6iGMiFW z^t!lp=y;(_^bInxf9$m<;QZu#TH>{|mBxb!sp4J)+KsAbA3giU?X+s<=Pf@EfDEWz3t`$7s$9Y@7&=^Z&u@rr)?gqvbWUX(EQ!HZ#Uj$wYU^AY|Gs9 z^c*Awja{3@y}LgBL_U-)G{ti+E8NBX^L=U{#rIVi zcV+c3o(OzxrsTBE&C%T2HJfU5CHGBiy9b2KrOevUK5dARmAv?Nt~f)hga}#f~-XK}`sb)F0Yz8zkZ@TN~cnO=hN^@K+5Z4ky*Q57cT zVTVi+lh?gpwQk3w%}@*ROg@J$viJnz z9f&ekxIv{B_B(7GjOb(>{|X-9^p{<9pv~8Irqsh`^%c87E&fQNi}Dp_W@4V}4cXv? z9G6e=U2YA;^|>3ex$J2i&Zxf92t3IaPshTR*VS+4a~yBS$hJc2>zvv!<8bNO%qU#~ zm@cGDwNJRp3ZZvUW&FIJj4k@#-bXMzNpn6fqg6T_3c77n?(H?= zGY%OJX(l4USi9cR)(q(9WjJLy;&IjvLK|1-3t@=B?8S^#_sv=x9@Kc&&3Vr7-ZKw5 z$N2r8pyVD31Y@W@r=?E?(@tV72&ITlT@)?LgyIH{IX)@FvwmAZytZf^glkt%9FcMf zv*I&1G+@CiP1(y{Su4GQ7qI)730CQerH`lZ^YZ|{R`IKra7BT;{6Mx(>%c0`ksS=+ z>YBG6g%0JzTc%QZSvlHtzs@$0`mLPwYrOzRGiJj_W~uO1h21c0A$aEBT*Z|SFDuyp z6@TScQWs?wHYjrDsqPVK-9!=v)as&ff- zc*B259EU`c+P`NSXYvPm+<%_NkufF6a`Kn;pgX%S2+x>Z7QRh!g&^>cFp*!72=I?E zkzbGq=1(yZ(R!qGkjeHH;2$^m#{>SsDJI)_Bd#tsGmZtq{&h2hUBSE9%y=Lk3t|*b ztAhfzI!2dT7&$cF|F10!b|3Fz3*!MS6^8l^So=pPIPw<^PLZADEnRA8Y<>ybk0OtZ zA$Z9aya>!c-8+CRb~KU@h@}Ey6`L0<5R`w*&i-D|Ouj{ficqA45g=Sm=2WepAS;$s zqkw%45XmpJ6B_wfP||u3dML_c6@)L(0=%dAt@-)Zq(hurIo-DhxBGcfH>oIUsLt(3$cGq)yW`L=OqvEf=xw!M{qe7c)vOlVSuvGTv(1@ z5XjHnlD`IJzvM-1qXc>)`j(jf6-)4;z`j3v=s(q>AxAKwIpobZKy4NP{}V^hOp0JE z!kcFi!AD);TmTXLS;!)O-+i4Dk%lWH zm3-!=vQ52Rz})n%FQZO*`6h?G_o?EuG=vL^5TA51=IBFkA?7XLCqXTShef(KQ|?4h z77KRzkUcze?wm`x-_hk5K3T3+6Htx|d-Kxbxi&zG@Y>RL7*$3L1vcK^vGbD;IWO(t zOgo}UwUeqb6;BC|L#Xfa>oTx=BU&|MTXDOtfxwK9Q_*sM&Hha49(pj)0m;I-@Q=Pa z-BdY!q{e?+uQSe}(5)?Q43ja<&Pt5a+B5oLZF94QB0HRDEt>cGnmdBOX`sjIMVFNUAUBKCE?2wYO6$7=t^c zue{=VeR#M)@XNB~_iavvbxI6YR2UKWa{N-wl<3*KFu%mym3a9IYXqNd#}(|v0E9fDuJ=F0PX!&#Jc1Q;@pO~j$xEVA=M_V?W+M9c8DJ7YV4 zxa{)P&eEA@^c9UTUvYax!*&qoCQXzNjRXd1VR2gTG4U$u^TsoCkQ>4mdC81obtkkR{L_`sMc;{ui|W2DEQ1Ml6ALCm3=Wv~NKm|9Q}E z`u#vVdD{TYlMJ%(46?s4(rtl%e_L`f45m^ zbq%RZsgXf~eZocZEiHlsbB6zSs0dwBl#*VAjD;={sw~*@RcbAg zxB;vs{)0~Bf8K}$e(ps6+K|7mqX#Us^nTsZ^ML%FAR72E{W}u%!E#h{reG5(-2(Y#{UTNja zWZ2c-x6k(`&7n_)lPU2+S-o^!-1(##@f;EF9xi&hm9&gep&GC}hxK?j|Cx_BpM}ZZ zZgWgkot#SQa!)FFUBN7N#eG!&V8j~q$rBc7vts#{#6Gqn|CH9A?AV(0)%1(}!}IWQ z(oMk9IMat6_|{%nsI&ym#BO1l=cnAEdtOZs9k=1}DLC#X{xLmJ&S?74b8yP&u!h`d z*@{azMK7R;$sHJ0x;C5`1tR~lf=0r+HHTd_tMTr3Q|Q%_N$l)3x+K`NWE|sZerJ8)Iz-zLS?}2!i@(jF>i8Wx|=>jI$JV^XwXeQgnCyHSOP85YrT3y1)QHN44 zkur2R6L5IHD8`izSPee^QkRMv*W?>r!1SMm89>)R3p0SYKM6BDAb$y~0R9ZD_zPJD z@MmDf|Cv>IKnfzI5AOS6b>IU(M2(HWZG-4NA$;=Na|eMjtElH^1Y(P&?xb&gVv6Pw z&J2F>v`i{g!u|5uGDXWesTaGW?H`s|wZLkrPF6$h@|Hx_r`B=$+vdheM{j}SbA9+M z$&v}D7R5(VeZwfprnT+GI$KFCX6sUNhbKPsmE~}seak4HyC6Fzow6%0*5rnlx6MdC z;fuJ%J@II7pqYK--bvKWY0K`|2PKAl$3KKz4(T24FEzmW!{xl2?{bw0wzZICo9LB& zB=L7=_*JQ%dj^IIFY|2qax&lg$C~XDJt&76nI`?L+Ri;bgIZF1I?x!SX@?p6`7pJ8 zc%Qy%C)+7R!Ae5U+#{9FU3XOySuJO%&nZbSaYIv=YyosV?W6AW8GO5p=Xytt`dGnU zN_vx^XGV!e-Q=g5&_49+rvh=RVT103jc)T3-~GSqg3i$kxof-Hwhy(;I1g=-a=gFm z(dLmMAxD={u@@(#T0V`Cu208rKPfslKKFhovUyD;$ssjIaJ=oQc#XMK4B+#&Xz1N5 zniOIGL_7oIES%W3M=V^e>1vOjUu*SAx8oV@A5X<931)lx_C>wz3OljRQUTqtc`Kwp z-ozbli=|B2xC~rw-70YSeDAjXfQP4dje(3yM{HiAmNH)q=NAy*Kc?>)G%-*Ljb&le7rdqS5(=~)i&6#|wYPIu|n9o1^<&-0Zj--3+=k%A+My#VW@5Y@t)DCJhXE+!f zc$Ls4Wq{Lyr+0U&$Y)OEOM)}8?G^Z_bXRVD+eUs3I}+_i6a#4I{K{nuy{KU&_`?=G zM^C5PgQ~7!7MT-DdcSZ4(g1aPO$Vi!imoUNBb#{if&d*Fk-89ZZIZxs$8I z+sdx+yl_?}Tvgjft8F5-&YofaMQvMB^SH_zcYN~<Oj2Zhl#zpXV zM_BQx>Q&%xU|8T}0G8`S!R4dB^sRrDM*Qh4;x9@g{&W@r{55IBpUxtHzb1_Uvxq;6 z@ZS#bz;}=E-wyD=caQMj4)DNtkMQ3P@W6ME@ZS#bz;}=E-wyD=caQMj4)DNtkMQ3P R@W6ME@ZS#bg8w4I{}1lTQxO0F delta 4143 zcmds(Yg80v7{{4+_A)a}EP@ItcGGlqMHbxyp&)oLigodp$*z~&f|TS8+JzuNK?^t` znVj|bill&-gGyk8_##@YVR=E_ympA7g_N)km?m^#mEO4#Wn5e<)W_k~TwS2PBN0EJSeM8u{SVBDO z+(od6H4o++3PAgtc zpM8q-ggw6MQ1oFuxlIXP6sqT}CKs|9(yTXb&~j_cyL-)D{#T}Gxi2$Sr<55<-YKuN zxbNpK?K^m}d2Dk=WS75f#I~01`I(Je%cfO5^VjuNzFQa&vs`C&-y59cn4bP!d(&fD zE-w6L^w)jyN}X>Erw;PTj*kz&5|ndEpJGoqFr&3At2WLjFZ-Cc@pj44CA}*TpJ@Md zisdWD7}+@@^U~t^X`Iz?s+mnIIpTlye(%eM>gJT&Hyf_C%{_k6?-Cf9P!T3?ik{Z z*TKA8hC-d$cuZi#B%noph5Qi3qCg81M!<}nWU`=_7M!J67?DMq2s6a7T=Z)x*&)`H z7Tl6za9%B`M!^;)n1CM~mr9YBi<-^C;uJs;%fwpZEUuMo9_AR@HcNx`T66zksal41 z#{&kd>St1`mUhe!gu=wad3~e`jVxzoN=qn;P{>P2JBF6#1O|90gsdn6rD+s+g$WA~ zEt=V%;doL(p&C2$7=Ybl@h!xJ=F9yFEqaAFzay>UD$KD=Rp3jBFBQJj$TF46e_vR_ zUS}wN3EzJRA6#cBaUB<28<>GGVBC@ReN7Yf5IVRzC=ne>b4@hAut$mHFtvx0$svrP zu@^5O_=(n8rkaIYK6aieO1Q z6ijbYiIlLfgK#dT_FHr?>FVMf@kuWv}C% z6BweV1RqO8+8*J?3;@zeTkvD?ao|%Zx#L|(uU83tOI1h`#)vpVkiI|8ORo->*x`nW zd59e%F4OXd-0J zQ!gXIsbc=4(eGydfdmtZu_eOtA%~aCqv2W3Oje7#j67nO)lck9A&}I1hpg5)n6zHq)$ diff --git a/test/python/wien2k/wien2k_symmetry_frames.py b/test/python/wien2k/wien2k_symmetry_frames.py new file mode 100644 index 0000000..fac6f17 --- /dev/null +++ b/test/python/wien2k/wien2k_symmetry_frames.py @@ -0,0 +1,53 @@ +################################################################################ +# TRIQS: a Toolbox for Research in Interacting Quantum Systems +# (GPL-3.0-or-later) +################################################################################ + +import numpy as np + +from triqs_dftkit.wien2k._symmetry import rotate_to_local_frames + + +# Two angular momenta on each atom make an atom/sort-only lookup ambiguous. +orbits = [ + dict(atom=1, sort=1, l=0, dim=1), + dict(atom=1, sort=1, l=1, dim=2), + dict(atom=2, sort=1, l=0, dim=1), + dict(atom=2, sort=1, l=1, dim=2), +] +r0 = np.array([[1j]]) +r1 = np.array([[0, 1], [1j, 0]], dtype=complex) +r2 = np.array([[1.]]) +r3 = np.array([[1, 1j], [1j, 1]], dtype=complex) / np.sqrt(2) +rotations = [r0, r1, r2, r3] +permutations = [[2, 1], [1, 2]] +rng = np.random.default_rng(1729) + + +def random_matrix(dim): + return rng.normal(size=(dim, dim)) + 1j * rng.normal(size=(dim, dim)) + + +def action(matrix, density, antiunitary): + source = density.conj() if antiunitary else density + return matrix @ source @ matrix.conj().T + +for antiunitary in (False, True): + matrices = [[np.linalg.qr(random_matrix(orbit['dim']))[0] for orbit in orbits] + for _ in permutations] + global_matrices = [[matrix.copy() for matrix in row] for row in matrices] + rotate_to_local_frames(matrices, permutations, + [antiunitary, antiunitary], orbits, rotations) + # Independent covariance oracle: acting in the global frame and then + # changing frame must equal changing frame first and acting locally. + for isym, targets in enumerate(((2, 3, 0, 1), (0, 1, 2, 3))): + for iorb, jorb in enumerate(targets): + seed = random_matrix(orbits[iorb]['dim']) + density = seed + seed.conj().T + local_density = rotations[iorb].conj().T @ density @ rotations[iorb] + global_result = action(global_matrices[isym][iorb], density, antiunitary) + expected = rotations[jorb].conj().T @ global_result @ rotations[jorb] + actual = action(matrices[isym][iorb], local_density, antiunitary) + np.testing.assert_allclose(actual, expected, atol=1e-13, rtol=0) + +print('wien2k_symmetry_frames: ok')