From 4c59fd606f76e775e813439ea4205caf63291d65 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 13:43:14 +0200 Subject: [PATCH 01/24] Skip sites check because it cannot be None --- src/gemdat/path.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index 4cc8a344..0b099ece 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -52,9 +52,6 @@ def wrap(self, dims: tuple[int, int, int]): F: np.ndarray Grid in which the path sites will be wrapped """ - if self.sites is None: - raise ValueError('Voxel coordinates of the path are required.') - X, Y, Z = dims self.sites = [(x % X, y % Y, z % Z) for x, y, z in self.sites] @@ -102,15 +99,11 @@ def path_over_structure( @property def start_site(self) -> tuple[int, int, int]: """Return first site.""" - if self.sites is None: - raise ValueError('Voxel coordinates of the path are required.') return self.sites[0] @property def stop_site(self) -> tuple[int, int, int]: """Return stop site.""" - if self.sites is None: - raise ValueError('Voxel coordinates of the path are required.') return self.sites[-1] From 40fb2539f20e56532024a464c3ab1676fd1665ac Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 14:35:53 +0200 Subject: [PATCH 02/24] Match volume with Path vasp_full_path is calculated from vasp_path_vol, whereas the test is performed using vasp_full_path. There is a mismatch in the resolution between the two Volumes. This seems like a bug in the tests. I simplify the code to use a single volume instead. --- .../baseline_images/plot_test/path_energy.png | Bin 32850 -> 33744 bytes tests/integration/conftest.py | 9 +------ tests/integration/path_test.py | 24 +++++++++--------- tests/integration/plot_test.py | 6 ++--- 4 files changed, 16 insertions(+), 23 deletions(-) diff --git a/tests/integration/baseline_images/plot_test/path_energy.png b/tests/integration/baseline_images/plot_test/path_energy.png index 771a8af8c5cc2e787a43cfcfd4017704d1899c26..08494d4b0903128eeede8249277b8adee3188537 100644 GIT binary patch literal 33744 zcmbrmcQl;c+crEHy%W9n-h1ysFp20yuTe)AMDL>a7=%O(q687WM;Bd0526K$-tukt z{d?Z^KF_nhKfd*?#j?zrxvqWfeeQD~=W!h8#6E$k;$l%^K_C!Zbv5Ou5C}2?0zuNn zKm)JDlJBp9KN4OlhF&^uc3!@g9<~rIOE0*yo0s!*D|#PW56|aru0p)xyaL?x4qjex zPf0#Lm;dtyUN;YWzR5~rX>butxSEkC1cGmQ|BF;6SNa?R371h~B#@P_GEQLp=2(SoQT&Eo|2bEYif)a(X)XK1s*f0) z%EmT?gy-w2pP7j*jJ?6+jcJ>Q2XLEk%XNa4lZ^R!+@=vvIW2zgP6?-Z-{6TyFbx z{i{X5QC zZj^E)TR>^n8o6_|E8JF&6~b-$`9a6u+d~ow9~G5I#(^2CM{R5;519|w)K;=mTIU!q z1K7QnI3iQ$>WSD-DGGmxs2*Y}9{${ssdbv;tu$_^&}xbdGwI5C{Td5GF5@p;qM5_E zR~sK6e}8p^dp|`T9Wi9tZv*!`yazXq`a0LGUHiu;f2H#Ch~=0cGSoZQf~n`m%)$7@ z)l$Vze*5bq7t1oGw{HEEQTR#I*2bn^Wjs@|sHmW%^uOmS63TWz`2^5Hc1#X&Y%Hvg zcnpkK*w_tt48OiE285xf2b%Z)SqZPUe<=MofhlsCSlNF4SSBSo%=7n;LE~rw!sCnn z4P0;+Kl#CPTj{?tKmF}`|E!2V1$XyirClKgm;KMOp3?vSTGsz-sQ=$0#gr+oXqMsr zDP=4sB_Uz-F>1&DZHD!M679>4xGzzW%<3tuOJIH9RpKq`7e3-8f#CWHy$?LE$hp0~ zRw~!dTNv9NC5|G59+};&7TLgMDbBmaH|-~gxK$j^a$;y9*DSYr;&mSE!N)i)2mwSB z3(ftpI%cNq#hc7ELU3nO`TLkPJw9!8$$o(2b2> z>*(OAU!XveuNzklSB5HlR470X{O7Yv^ccvyTC5RQR|H6YqjngI?COsQlo?l#{f>_I z7ryw4!~0gdz@Pj7j{KyJqcdynw}hg>tD)f@JEb`~(_CD%XTmP8hX@IwUBt9W=mwix zjJ|zX&K!ILi|F`P8ZRU~3@SVInwyy=MZ)gm;bnCq~%BYc^8uPDCn^W>J@7_5IpLRI910ccnHs)t(|e`H?y|QP8crCkp4&@)Nzo zAcD$_F}-$k{R!D>;gR1JWWmn;nV6pbU_|Or`h<%syu2vv)cNiR1%9^B>C%Msp$+W` z9h(9(brCyb-!`^WO|pO6Or+rbr-;%P6XZLB20Dac5~zM_sH>-Qq4otcqT*rn!U_&c z(htKUZxsHSZ>!3x$PV+4#E#Z85evem$L{S{ zMV>c>3?RPLzu>Sz#Hl-n`A`hZFuBdg9_>RFm~k457`ds%3E26f{dS58(D2bRd8Q;Q zBccMjnnC_~Kt|?!VXB~j#QHeDQX2$`1B@SIe3B~4p z!XphPAM4Jq|U7uD&CU-m@pVx)DW1MR*-j@Kt%k51VH2v1|bpof6Hj z_i0eP+i#yd!h+!Lp? zVMazqvU79qAGI_%%~@TY?8vld`h%@MK~wV07$Q%vt8UuujJE8DQuYDeZUmy3tw0PV zRBb&%PA(WXc7o7L&OUbuyKLrs)__ALjwom7DX2a>% zLucL_k^af**uOg7cm6yp67^4T_slOD57*#X*yO6cdibU$`P3Vvo5#pVqIp3}-tuJ% ziYSZ_(_6lB4IFYv>z(A$T6ASEj37exay&TX_PT|)R3rJaW$T}JnnRJD zll&{gh#|bqP`@99t6x<^wQ{ZGXcY){|8!$_?v$cz7tmcvq_0g=6Cot!AE;C*(POys z)5>5eJlG^eK;VaMBjk7eR12#LTwY6r5qM?QXtN*e7e;?MPuA6kS`UOoq0J~VNrK~E zN17I&HeZJ75)ZwaT$&vGx&MzgNp)`mLs92jBTk>7@HA~Wl3_iX8Z2si_xz`ZbR@}zT>s8}fU};TO2qy0NP(w^2}4U@b{08Vpl4x-*sUW$A}mB5 zV(JSmM3=oIzxdVi3#>4EPq_rMxP;Apo+SV82gNG7juvhhDRvkKM;)#2KW`{JJ?42Px}+7+f#ts1F!_nn20N5;(yK>!Jq}Ah=3W`{iG=x zfOvoM1cF2P;LgdRt+)h?F1b0LbZrnp?JP+CO^}x{L<^EXG}!f@<3-UV{R6tGY z3YeiT;(&AHfnkNBZh#>wYdhO~hmwqkH3wg4#W@~x$yxkwfOMcSLKVsav7 zEEH&7s*KjLlYHiTd%fWSA&$c9%g+%SB#Nc?Tkj(a9noC$61t$ zRQ8g&#*w6VfHoYykpS@BzT)+;`8^;I?kRf`MeNY_E0YJMx&?1Nl1y>F21)#nrplNK zf${@k5rO*iMOT!OtuJOyO(?2K$O#F~R(VfIhTy+Py23-o(Je3F^dcmmHt<&6g01(O zF&GJtFD)wG>uACo*Wfj=r*(q1<4Od|HaC0b@-95Nu%rirH+gl4*B~L)0g;#&({5*HQgJM_(_BRF*LVCli zj<44};H1}@n96UCeP0)X~`=U}OZy z>qB~$khX*${O9(wU>SZm{ZTJ*MWhOWF&OOHRE!ur3+;KAZbRBI^e;E2 zN}pDW+y6H>;9T<(tX){DHjGF8v`xEucO5GCmx5eat?dFs5k<410i!Hq4f=&Lmj2@w z2p}C867OfV{EY^4H;RfdwTO3RWf&?dDl|d9tDT`WX2Wf(2mxy}P@fcelkXBcDz*uA zT;dYvK9zdg=`($Vb=M)I(!WpP@6#}N?{Mr@H6 zncem?^9Hywa9_d8gXxr<^Dy2McmF)$lMi!NJv6{B#zA^q4&W7a)n#PGfxkTY~z_z+Y)atyEmCYX@k9cUDuZU12k{O+=oC9B3etEf!V z8CpnKk=qUqa4ChdVQfA~mMN)`i4ik^CL$v8p7nwQY@Z{bY)Q(3h08hDzfjKhk$iro zi?N1Hkn7DEH%0^KT}aEq0(Q<24z{}5Oiu;p=9t3$?xu}Zfs;Of*dvfYo&39BnSjz7 zY^%H{S#WW~%c@%KTTf!;qT}^=`G;urN>8z#FA|)(SAfcXnV@eTT1pT3vxRinjve?X zp|zGUA}ag4>j6j_HEso0+pIBAMQ2yEyE|XI6W0b#U&e2x)%Oqhwvn{#eSDSg=I6ca zyxvrny8G>6?bR|P@-mB_4;J3QS5|l z#q|I*a#OaBou-IK%@&GYc%-ZQg@{zaA$DRKYONG58s<4c-)Y0xRnDYL`Fa}lf7DB^ zVvUe5Je=JCM-0WnzKOoOAC=P8(|8p)a^>1*U$@4WG4x)E<=#o3S!ECpnfPrmFkxf& zWoLt*z%PV2p#cu_ zN{2xB;g+Uh_8A_SJBzN4Nri6bTig1FO$a>XIV)?jzkm}Vs5egi1rDRA_%~1IYTN3F z+~}lZcs4h`nAR#TXi^Xow}`s-_FB{o53r(!&lkTd_3{z@&qS!9D1}Y@@eHs!bMD+m z(XWcTH3Ykem1h<#d}YcfyV2j9pH_C1LrHT48970n8aE#OWW@_P#HQeJ^BKkvyaJpa z%Q3bfBjn_Cj*6iSXCi57{iaT%7nUbNyiT7+B0?bk+!w#^^uW_%DG;oRYUjA^IfQQrt8AOxQ(AT%$5Cs^q z>~TL2&YiM2n2PL1rpwv!YPH?!yF%jAPu4#!hBB#!WRil|^a_EpJ&s|_eVXN&H#f+J z(s6+bQmH&W{faxLqKc-m5uC61?x-+)_Hw`LH{aKs0R~lEV)}~r8(arxFE~@mNq_hLk8f)Ue}CP z)F}fDum3+VWdm-Z_>_p3PMm-^3bzLUky_Vt*Li2rO&DzJu)=zG2XkW}N!2A?XN<4o zMsQ5}&q$eBZt@q~*y$Y4an#?xLk8&V7hs>K*u0{)D8}w?Lx0$SDHXIna~@f&N4I*; zfz>(mlQDZfkdqC4FMTAPP{49;i3&T324!nB{9-_^A}hA!G_zvr?dtayU%rQl)UkKB z*Q#EtzAG=uPESv@(KD0TS;Uc6mjglv(?Smxb1s-MXe=wMRHnY9d_$$j5Ekzh;Naje ztn=YnrVzt>EF#J%!GVcuWJE6HLOCn*JEjgqhUlY5A+IS{Q#Gbt0@NPcW77bvMk#Z@ zqEiKT5gL4L!bK=ksI;gswib-^3pcUnMEs5qjpk0@`{mjbC$FnF6%V{fIz^27lny_~0As})B2pGqHpsI89)3T5P=b1F(_`ciQc_wtI8TqA zn0BvO@Dka4=wVT}wS|OS{Ta76;<}PLPGYYU{h-17b?atfA$?6Z1zdW&AUHAj@AZmf zucCO~Z(~O9erL3;5u~`dkeAnxkk%cHUjonDW}kMiK`5zma-m%`q@~2|Bbpy^sQOz; zoE$&m#r}-_?XkM5Z~$ck8@5=r>vG1x+_F0YlU&dSgEsI;;pV6Z-?2Bt@8VWoC9;kJ z2KjWV@8>Kk9oR>5PYsy|$=)ZO_oRGLb9;Jc-=T!4CHh%7qZ+%D4GXW;nJ&7G^`o}d zde+j=o^g2g$w~Rt6j=8=cm{!;{U|By&+hvV{l7PdQ*v@*Rq0+nkdTo0bkvP~eDi1g zZFMy)ft)xz96jf!t|8oCnb06d^ULin*FCKi5RGpgV9QbhD#@)SI10e7QlEz4$kPLZ z9J6(T!nkyJ9l}!biebTEeS0d6O&5Kn7JIk_}awA)^_b`?}MtlyFiAWwF&@$KGL-oNL!zlvY{WVKu!7p8lHXVwVq+0 zuYGa@CFFRh4nNy@>q0L|<(UGGemYsYN)i$WT3PD5S2$+?4&T&UNBeAgt$w)b)GJU>ckC#2bZdxnS8xj~=DRiL|G*_eSCA zoe2M*^okXx5ZLCDTeQs5qtcQ(TAL~cCuZJSeO>vZ>#zH8dEwj~Z8=O>yL)p8 zUGTI%{L$7SU=6Ues~R6bkuEPU4eRZ2D=I2#{7)T!{P?lT7kr8RFoBZA_f>4@)Zp@s zx_Vrm;scd;7l9bIRa_9Ys_J~CkpU1A4SUyOWKkq-2Un_zj%0Rv8} z8Q(p%`Ju^0$*qI#l$H45=B6c}1p(9dRZ?2f@|ARUy9QveN(&rmP`_e28ByGnmt(vu ztv3f_D9DovC1gsBm*x4nQjjw?6`(>=Gz6cWo3lRMogo)Xu)M=4S2WCl6mV|ISo{Z1EaLuz7ma3|fTf zB1n;1?fF-0*LvJHDvja>2r!*8`33+n3>KTx5u7ji(xKtgz6uLNZatkfVP|LW>`$Q5*47SX(a2mqnR;x% zG<#9poYS_ZM}QE`#ohBc><8zC0D}aD&ZylFT)_3r1wP3=f59@Xbh4B8OZ~FpyJGCE?>BX~NNv z0kQV&41>tD-jmHK;7MUZr79qJlH_J3v6lV^4SApqH}hMIyx7WhgMWKoszPS^3g|8W z(2_sCR}`kAZ)4)nU=0^4Y3_Tl&qEdD+XmmYFt91g;+XN_@bDaI^w1kyH`MVx8iIQ$ zEk^inAZu9)1Zbr4X7{1YBDnz92&aWVm05qU1h!~YNx%p@su6ffdqV@memXo z;UuVy{tA+5gwuyMQgnRU2zrccad*b*wliZ$=gae@=-00)TmRl&*P&C+iRV@7mNpho zCo@1G^}~I0Wd2fiBl!t00?y>Qxs*cL6(PgUu8uSw#@!3a#PXsr4=cc+c6J#LpMjl; zFXsSuW?Vxh&S+0mZLD=H63YBvHjo18f=7C3_`;@K?e*G~q76^Op1g{2PWkRCRszez zitY)O2{^Op3PYu#?KtkIs(pFwBP=ZZh0^Bx+fFRC_iR7p z_r8y7G~!FYy4gD%HZdh7h^@H9+fkx&nTJgEEdbonfistaLLQblo5IgdRZU5VId-}? z_qMngl}gODw|P51-mJ~v+wbP%6fAsv+es*tty`G|l=eSY=oH?QkS3zuhvO|^|IbSGAKD1u#l0B54GKrdb;WuysSr2^-JqFI+LLbOGkrad7Ua` z$6q@zCmI<@V%7fHi9ONys>a5Y5CngvWrHzNNS$~}QFmq`FCp||-=|s#L1dSOt;M9O z2+4i4NmA#YN@nfCgZ%vYZd0yTv{zMTmJ=mf*D@VoVTpwhkerS`YMw%1D?Ejv&o?qS zc07rUftD`04z$0#o|v4R+}vpKJ&dBqAeZziJ`ocN`lG9I%UT``zBAW z)o%nU+3zbR1ijdAc6RvPU#~g4ZA_EcY&?(vWc?l5x4L?yctGl*8x9ZAQ{p)6a1N&A z;ZjSN1wB!A@Y#oe`omssN-{MOl*OFN%;NQRbu}KFk6KSAHOQs?_|5!xSoW7%$gKNc zb^ZMMcJ0*vc)xy}lf)^osCL!4cF}T68r#CF#nQo1QNu?xFem33_77JV1O!o#(Z|Xd zTeCU@z?B6g3Xi305D5Tix1&dW&NTYvaAx7o$&|>eOlxJrhqM~IkQTf6xKF0o zxfmF@xHY=^J~I5!hYzEaMKoXL+Tqci6%uf z!j)cswTpeQ`$vbMOm_&o=!+nN%)HW@*;EOXbQ9LU3?Hc6oPG%IyE@%tl7zE=fECfW zitRoB`UG$Um1BPz8R~~4qx7i~8ecij>mm{yGqS(0r5Fpra2D1K0CnHFwiIIbr)__PF;Pun7ptV<~GR; z2?_3m*xgl!!WXLl*MfZ)C{vPuRsfZqH0a*3~q{i4f{?cV*9U z3iEUS=KH+Pdeh3D3@S&Ciq9u0lgBj46i~w9RSk_o(_3f!Y_TG{+KCi&EynV5OW9ow z?ikxoyi4Sgwn7~z4$chawDNET$-_OK;`u-s9NfrW8^xKI=WCsy+?L&&v!1VN%7+Wl zv&wI!4L4*SxQowgF*ReCB{t~&j0VInk*H^1-*#^z&HlCu*>ft@$AI;|zXv8?SpBI| z&VahNr!7marWoa7&_K&&fQi0nzjv|fhgKr6FThmbvyqT?|D7@TFHY-xb9-YfmbljO z*vhZXOrm=W4+!;;Bw7TZa2KaW*Sv`j4f9}7=f-?7@15QoQ<&w?{}=Z{T^Z64(=2)7hCI?U`SY~<%?!E&J~-A0 z_R}*GGW!euQrmv$NPOECV7>f`BKxUs#6d=e-{Fy^^M}gb^f~hZVC;y|_*x#1!zTX{ zK{?T=)K==loQPUV)IG6ryjt=7*^g?&P(YvgU9T#DYF71E9LXn6F$1cv&7N$3H@!mI z4De{Iq*g%Y?#U2S#)4@`ffZdIT^Z#iyUEF^!)>{!t2-(ML-JH6 zL@mPD;nIgZA;7U}pNciefAKKW^BCIX#KRJWYb~b^uSlHP*>6VRLCVYJA7S~1>2O$H z^oF6{>y_$V7SZFe>Rt4|I5@f{(dp6;Rv5S|DnjUjR)Z}tEGWh-2+FEwB~w4K6crar z`Tpt!smFcMrOEQgmy}0+>Ev*$Smv(xWTgw6l3iwH;g{DjgvJpwJYO}-&U=^f<&@dg zjREr$5oS`d>i$)H4P=H^=t~fD&EYJJO+nWFxNwx{6 z;U9LC(^}FCXnz1K&2uYVdyUr(V1L-Rikn}1LCl#=rEEkaHI z+S6(P&!RVG!z$4pNS5o$x_??$!?0I9nrp23u!KcL*GqC;!vMk?$q_mRd<@3}E5uVX zJR2Lhb;R(7|8VluGRSWf6yISL!+C3J<0#w zy@_axK&CJk2Mruv*1)LeLKd_UZmN7siE z7NBq-VJ<0;D)@V|ZCGoCibpFQL*>>=2$(qu*N$7C`2(NtLDvUmWo5}+V;A8{U6@BA zF;jVM1Ru;0CA*Fop&|P+2qC-i{(B@xf#YJX%NTNwgC;E+DDH|7(U|@7JfXai6mr2= zE=4V`f;ucrSYg4#r8GQINK-laSp|qLeoq^V()khp@&1&+3W&jzOJ}5ougLvY^FmVX zsW#p3Xf^{ng|8Uf@k31jtUT5L`s@RK!5SdL0m{Z=XR`RyLF3|4Kb0GR%hwcuNMbXn zJVyNxmi9J&^r_L{l?YEJG3|zTo2>~|$ViwF3k#~oY;|8iWaJkrae}{DDWRB3I%6YAcF%JxllQMYd*H4!^%=grYqTYe@GQ0yJCN$&gOi6hCw0q+J zlE-nVqK47&dLw?=?C&g5?fLwQ&>G!)4q|&%zPf$eKVoEsgvF|h0tcnDc?LRk13h%xcRN2nI}oYn}yA!n%KLMB$t@$`Wf@~GkX9o zIX}>-+|e9OG^~}SI61SMxmCk-J{qp>trWcPDy-IAVF@%c;ujU}d}5wNuWG~}KO8Cl zC~jQ!Q6 z8d2*DcVkYjuYOlMrm|aZhkb)7tNXlMTw>lYUXTMn!&4w*E(0srCc@FOqG~gssyJc+ z2ZvHQOoxTL>Tg~yQ?t=Ff6qeH6ebQ;$<2=o6#SaCJl|!cX>#hT3v2zsDV@Kpoy%N$ zdfW^TT$A{>)ilty($g(0_Y@(7^2gv%bljWG7TQfA_*w-E3X%mi4^`$oG(;gBjX^2u zG{{}`g#?!n4o!tOGtH{w;uC%h9X-9Dsgk9OQ6qYKE(e&8rqa_;Kt15R*jVl%E3y_e z{3O8Psi{cn&_Bgz><>E`35uUn(T}^zXZ<(gd9BvnzA}1$O9lX8aw-%-y*bXp<|9R% ztD5LQ6YG5R@g?W_cYh5HjW5B0vIRr}{d7QXNUJk$NPLGPZ7Z%dGsAyz%YZxml7eF7 zChu*&&Ko^=MfzF<{Tp6B{PM*{_qkz8E>L0Jop@FCVBiY|8{h&~(-;jef_%rGEgC@k z^WSnukfT91;GP&Na2mkN)ZT1EZxNEV?|J)=Z)biqHPn-KE1d>&HR+jbc)oIP*!KdAd z_eU$B;!KjgV)xmf2f9&{sz)4g_P;IxqA7t-Jv&gCr1~B^tkb@!XJQih?P5Yme-%Nd z>%Y70nXPza|3lz8vu#q|#XFSPrcX`{n|4`z0*EOgf{@xqeL$-2I=*So&mFZqz9+E; zcw4xCliAs$$lWpr%h%G|_RTz!#XuhRF!Ls<2@p%peLwO-_8P~*aWZKwAmE*g5<>UN zs`xfFwEBmQ@M7t``M}U{e$CfuuDb7JXDT4jMfrJU(CwAwREd^ho4<(c&5_c*r5l@! zw|AJY1E+Kfvwq>#$%&Wq_fLdwt;c9SzrM$UBq*1v!;FTo6EA*R6Ny4-UtOu(X{S_^ z%LfOMY(N2H;Ym?sFy=0!6GCgmmwj@8cw{Z@gA>O7yJPmL;^)laVG&%WWy065EvhSD zjS#7a4d9fX1Zmf2{f*-?G%rqF7@(1wPna7IuKq5O4t$eW=i7&EL@3fO1#xQdYLh0s zHviM^R}?}pze=~2_Ioq?@rckk1XA68CnW6uinZ<<48IKU1(^bf zO4;HB14NH<+Rv4jKu)_mGvQ4*dRypHRcz)jUhLHDJs?(yGyB^Eqlnx@J+IWvp2z_g{;ML_ ziX@|p+4=|(CIm2M7-zeEeXGgEb#()pmf*5!d?}(`<&){3m2dmVl!PfQ<7;0=l6{*s z-eu{@$>>u?!Nuc<)2maC6YW`;MEsLG`9rj@^7539?13HY2JRF|YeA&*CG^%yR-nlJrVThW zds7W70HvaP<{H#E3kTuP*=WdZh;K2x=w`iEXS-W1sS*+PCyuMSbH3L4o|?uaW?jV( zzCfW3zH$HlVml8mPR$O`J~+wnAQu-KNwsb(vd_G}s>DI7O%O1jE_}*&fO1$L9-Ecz z;fUim&_8?jG^2^;XRSyEr(b=(=v zo-gEwqAWblhn8lx-H5G>2>z*BnX9X{Mjt*a2&iry&`T2^hBs;(W4)s2iHa7!sss^; z6&7;78~{eyLw5uQX7#m5AoX(r2qlpw9Pbj4g3$`_+k6iu?}zHq0=&DA&84!{;heB( zAjl9KK;+C2T3VXZ4l~9X>Cl_I^Ul{Iop-mtX?zFpfaQ7N_%v>5tPwZzsCt^GxOgOo zUc!hD$bo>Wu7R%s^%XWfH@ER6g{kbrhk=&E_IS_JK!dlW2w0iMk;NTvUdKSH|3ft8 z@X%bG(&mas^*E(&M;5{f-ZZ;44K$y-D zr=yX+n1OzsrX{&dU~6_qQv;Gi_L{Twt=GlXxqZuw_bXnF-G8@S6sU@CP|%cv36`PT zt6kp$^v?FszrSFsu|kEPI5-b{t3?std_L6P2w(fwx^$J+Qv1UtrSSE@Lhi~g*fL6xX_ zEx2R$T0|FQe7@FsFF7#))7$9j;Ubxu2Gbwk$QsE!?iPw+dWy$kI$Rp>s)5&TOZ#fB zqu&QXBXwd9L%-NN&j_Dti$xyR4-M>t# zaeQ1Fi>aT7BR}g0ROLGe1%G(2#LuI(_C?XUx*|+YC~Zl>G3LDW55fKnc4L)^H;1@Z zm90!wWPon+sxF##L_;O_{&+;PjlEEx<6%apZBC0cHDl zfKd)gq+oJ;X56JC`V8-p2JodU{B`9M=ldi-!w<&h(?T<67EtxHwpR|u##QS1i|=Rs zq6n7`HZCr3k)AokGOsyp_bpV`uI{^6n4$K*WA@{1B>(_#?MqdYBMq#VyjW(}pNbNB z9{n;3*I1R%H@o7!W#xyYTeOO{%a+Wr`u2RIUq`w-l#&m!W%G=fmtNTaxL)$zx1d-& z@v>dlp=zPI_$5{KHaS{82Xnddg$j!%8380xYKslMw1e2>(&|uf+k`Rm=qk00G~=!*g)E z8{6~Wxcp;8@nZD#7jr65<`W(3K9z$*jxH0iVe<8kn!RE56%|+pW-nL-1qI<m!>6Frx5o@7}P#@CK)Z4^xg8;8X-jLWlMkBu5oq4{;Ta&4Tyq#h3?iQ z;eDURgd2Y8IGK5|INm?{rRZB0ee6Oc6vN^8*M@Lj))x`ZsMHUcI1EA^8Nt1Aj%9uniU^AiqZb;E*Nb%5Ynbg8^%r?0@!@Vbr-9LQx>JfPS=y8?Y zM8Gk|!<9Df9!l0V=()ZYiXP(Dfhy{Xb9i|?+kx3T=RZ986cC{@Pp^}PZbyKwmOC27 zw1mS56otRx^cX)^I)X<+`)DLVd&IOZG&N0x(A_SJYyVCQb5sa8`6+fahifbyXhm_# zCAr==*@~ZH88HpU6;|VIF{!RUhX1~i0i8DdKHN19+T${>fQ>!iyuV;IU3R!3=))I< z0DKr`g6La-C_2B<2O3D7SwJskm3B|#~|11k!ThSM$mIArR@_rwriXs8j> z^OQBOX2wRc-?p{_UFF*0O$(cf0E{QsXSjNA^mjGBY^g3ONoZ&|+FLO0TYKZ8BxvPw ztN4nmYj7&F@Cooi9PN~L2H#EH{_#1mDrh-o6|d5_JAQYLQ_AkcBQD~2e`x0eF!^b( zx&=N!<7@Hl{uXTMJpu^CLNkfd$pE{+4e zZ6n~Q7aErizvol~&VTh~@|YYg`L0+Frn3B;`#={LA1~m(uFSGNE$2^+pv>;b#iz6D z5lg`u)33FkEK5o==vshpV<)Yg&TPhpz-i1wPb z{g8NbvQVTq-ehL{{=PA($pro7$W+)h>aJ_;QExY@_~L~G zCl}W-=$6G75fQlnUKjW0mG1}EDS$^XgHLoYyq-nzc+1I@q(YC*Uf2avVIkb@ImlAu zvRLBHY{_sD$nP(GH4cws6bfOp>NnPBU5Su<({?HKq2$oi{c^MOd1-CzX*FTPPQet_Di&az~7k9c*&6S z&_FkHC*9o)id5V#VT~Z4OJhFnDq1tRmSjdTmh2rUio1DV^ioYAuy?FK}#nl zAXBy0g!Z;GMY3XKAPqKxLfLal%Q$+9XVO?>=^G=`$l;g+a~B7SqNflyI{pxs_|z0q zy^%?mgmZS#qlr0Fbm!LO{M#m_rj)cArwjCqOhGYdG+AD_{^ z*UNx2k|KP+yQ71IC=tWJ8m12S2{pU&`Nn_=iytn|Uwu_5yd+v)Ehy~|Z&cNUJ-kIF zqxas^|Aktr{q1e_c3kWr6o*1FOb3m7bCV!U;RR|by@LtoOizS=fF|o|`*n!aFUZ}E z{c5^SqFkdjI{d3jAm#jU=+YRGimMp$tp7Wyu7&p7(%Th11xEi;=Y=}daoKBbyD>!U z%j+`+0t|Z4A9Ydc*`EgKXz^^h7eE)atF6?#o*GdccP2;ZgDkBMznApqKl%gM+K^HC zqVEBI81Wl*@0{c%B6UTLuJU@bAV{Z{n1K1v0Ik$rOD3ts`L7F}Qj~#3zLRJAhw5zN&oLqt2?$ab&p*n{kplWkyj`n_=u2 zM@@LV$1(Q%hVl;{fB_W!6T&3+VI!HTY22%uYmyxg3zk2s5J7~X{P-6fuTxV~>NXF0gglZ!z=D{-6&0Cbuf^SdnyXI3{*n@8Yk#JxBD!gTvumK+Jy33rEPdjqCjQ*Wc)b= zcQxShuj&6Tj6aac01bk~*rXPq#SZu?fm2<;ZdZ52^dm-Nb>)*2MS6(S3yDOt7sWRN zV`Fwy!IV@;%{ZGq3ejyhME7k1RVL9AUG${D>G|LL#rLjnNwH*1m~Kyqoi_HWmxk9M zx>t$lj1Lgt+XI4+dx;-IAo?H75@y~T6dX}25U9Z$2ett;_j+Nl2jm$^ub6{%QA!dS zgWR0}9ANR~UQS`MjepMus0xRt$AMob11HrTXr;;TIZPgL5{|-(ghydh|4Nx7fI;ZC z;ICE=DXt@m0`B-nyO7~$-dVRPRHej#TfKcNly6uA@eH~^O8a*EPya!L;nT2UeA$@D z;mqeG`5r$PZcHi3cdKzH*C5M5CJo7qoaNiw+iyV&@nZ0qhE^Z+ZBZCS$DTIa&c0{$ zyhm5(?L{*f8pgYK^x3Z|XZCI>D0Xu@;u(kGyKmm%DTJ``CArXGp*19yXuN7@IM#kT zK4Oz26se(i&w@(;^U^e@Xy^vxYQ8?<8-vQKTl+X{n8?>$mZwraF5hEf8Dz!oyBHcQ z!-ojVBl3k1OIl^Xhzgpp7n`b)32ng{&V`0YU8c)zTK7Nuq8zEU-)^+i432!?_vn(d z2r1gG)72MTBRzQX`YVR{-@6Jd7_fQ)DgQHjIj{)`ZDbfs?6n+D`1Mm5M|&cJa1Uy` zT53=w(A`P+*5J-W>cfPb|Gfrj;`!_#FjC8Vc$mC@Rq4PY?)`0;FOKnPplrmI4Jm~A zI*>p`;4UhzzCLHAjFy83HRw40SHV;}rRP?TB~9ah@3x|vDPbrW|HpbPcXNaSXf^kP zhB;&iKpptGqpW)s!xTbx7DMSAbu^E?L60D+uuaTtin8WnBSRK;0@is^AJFv*8%0uu zpY(0gp7J`7%TV4;>(0wwgh+L92#fXnJY^)nr+~@3cTpTNcJ6xBg)*>M{nL!xK!}(s zDGz{A5J2z09D#9S^I8iE_JQ}M>>eTLF#={_Svnyve6}7;kcqJWH>KLkg?tZ)aOh;I z0r%baz1j2lyyI`%K`nqwAG89AszC01qY%*Ec z%ZMcZ%fpbth${_V1m=Ox#Bi3l$mCE#oy zI2^GhrGtRk(+W6a92^{yMbu2+P6=W`us|e4oUT4Gd2d2#{hDpuH6h z<-ht{hmX>gu3(cUboSza)M~k;6n*u$t2yY_|K7m2(VrlG{@wHI`?TWT=cRdQ1>lR;=iJO47e~N$J`E+XS|t9)7qrvK_leQQM-4p04K52B+nY7trr4& zP*Uy9$UHmbfq9uB2sI-X8UReeK?5d#qu_tAA*K@x4nKkfE zLd-rGlpL)?N{O@i753PRO8Q(4zQ+pznGe%kpA_sQep;!h5Ug1y|KnBX$NVg)%YVkg zT^2>BBffgRC3Rsr^#gqY`&R}I2jx5*n0uP6+%ytIk443pn#(b9_5Pp6-aDS^|Ns9t zj&aH+d#@7;p`2utm0c*=WGm@VM)t~)c>O0u%WNl{ipL&!0*_g?Y4pI)!`=lA)1 zzJGkbUHYSoi;kz~cs%cq`|WnS-mkMl3)5r%`@9jNO1DT%1WiGM0FFc5gxB*%?4N;R zJmT?;d94dunM;PNZ>WI!XSSX^m5iAMWQ58aTv4PHp;rJGR|Q0&^mu^3aD?z$&-CV1 zAPRliIR?cIv(og5JfyG&@kN}luD;lGdjf&Jo|RaRh?LKYp<_SJyi##)o$(Y%NFpiU ze13hqiUq3$YI+$LI_TluriXvENh6!iS~p`SsaYcz5*AWYq|yPmyzS-Pp9&(| z(r_K?1Q&exn-WS6(UVjd)``^IaE= zDT9ezcmj$%PTqk&-MFJN7g1l|50+Ge5NTEfw=dOWxf4zvobYV#{FVqlI-o$it7fdf z2qHLyRlT@S8bTn*;C_~(~=q2c%68)52qrn+^fQ@L7>m?vC~CiS%ZKLyw)k zN7frzwz>J6ihjnsRE76PUhoO}DOQg^MieTR4+UE^R6XMN6VxuC7kK9zF#5KHEV7*n zMNceOqe?>bSZ9TrHMay`QWE>1Ga<~>65AI@tB^K8oI9Ih*qfBq9fxd$kVR^avoV^p zo%-?bNq{_5_K=7-l+3q#m>bnp2h56yqYVvyq?fRVR06?ec zJSYmD0w8AuoJthrq-j~En3NeDt)$q&lS128dB2&4wWVWG?$pZ~6aHepo&nL+Hr$Gn zvTY5it}fB#QtlB0(+-5yiC2sT|FmO_!(&-|=gKA1$l6vorMBi({Ij zfv=r8m2CxL&9%!I zxvYgLUQR~W5Ty-=L+(#u^}q_sy0^eO)XYvBB7MQEyd=*YTKAx}IUELT6rE~6{?H`+ zI@Yuw({rA3V~IMm&EDhFhQgB8ToRAl?Qh}^Z~riR@85wPk&oy7gpRqRy>$6$Bjk1A z0il5o!cQ|_9Y2%lG}C>cV2Z1|%X9nBRXPT&7Ci;L1QZ*x?5Te*sgX+ivyc%7WK=x7 z^}wx7CEQRZbF6O1rh6a#>`j^*q|UW$lvPT*aJKQpq2$D}!;tlSRt^^P zg_EBPuNZ_&!-)>^3t5jsJ@)_fswVu1JuCI=>HTeJ+I>BNJ$?G({Qi|Q8kjkrCT};b zG>tEIXyQY}be!z3NP>5~xVV_mczcqWImBIQ`ic-_6ucc&uLj~|XjWxXnSqVL2>(+d z25`@@1nbp)`e8)xl9)m*V2U6gQge9yjRumOAP{<@?T9% zn(ioD=;Dwk#RQ)`F?NVOmADaSSxpT_5_q7d$*zRl`(A2c3?UV0!Cy5uk1aMV^z`;S z*jYtm@I<1U7*wu|+x;6HP(PIo-cN*mEX!wH!2ic zeV3w^HswTLzLd+Ok<84GrTR@)_~oO#67mnMdXRk3mZ}$e0?!3VOeT(5WQF}v5~kav zUd&iR4+7{$p@v}gc}uye;c>RNm+nSLe_biOzqzzEe*vSfi@k@}EA-bYKb})CNz0&PVZbXHTB<1Y~ zbw1GS*UD^8PnT}>O4h8qN7rFnQ#lV^#CeT+lo9R$o~QwLJrXJh`*e%&^_~l~Ag&qL zxeL*7qJ!Gn)csa)6ogbnIYcwMGG$+@Fw5{owMr56o2@*)XNW>DxWiVE)r|VXZ((kE zv=C(d2r*4@Fl@Wft2J7fMEw&UCjnB}3w}B1pst!3$g*ePBijx@xF1}Le6Vg4bBV&9 zU=ya&)tx)8_pAv$QPP)MHiJ88dPY`MM5G=YxwHvJHSCXOt6#N0>C2^tYuyBT&QEi8 zVb3#)8I`nvBA=L(6AxMrHE7BbQ=fBQt77*6z~d&9F#{2aFB6F#W$(w`aXD$i`-<+* z2&UTfg*!y1dIb@LCMw5|gdB1x#l z<9oU0wsG$z6N^o_tj{7YVS8%#xHXJ7%6?TeV(%P_7HUPFG&PDh#yp2*2O6dyz&{z2 zZs~)(0VNbtX*OnC52F|h7%H$$iBN|PR%Q~-k0F?(0%^$y za}BX}=w^V+B`cRm0|LnAdpLD;_Gs@9N|k$+KFOa<>EiHyi)YhbCCrm&uUFK6omDNn zA+h;d>=_yA4qG(sM#9*^ENV;I6nEZ^(k$a;EOdlBc-u!-NY)HyW18wH^ zvB*;mtW+dmC;%b`)-7EA%)vNIm(i@%5lO&02XHzBj0zbcEVkJwVCh;>neo|~x8*0# zWWn2Pp4g7mqr@wkL%UEe~Rz=7~~|=kqPbkNgy>yvNA~r)@VFmcDri8@~=JCC<>4Ol{Pkv zL^`khtKEj6XY)Da>E%JCv9kYEg%Va&1}1;HWq)6I9t#1b0f|CVLdMXg*U)vmEL|yY zY=phIDy8YWO`V&8w}FBUZWpXEc!KdYM&fNGJik;{0lW4?whf+(@4)9K0?WA|4H-&D zc!!)S;VwplX;9JqVVK!aL4DWz{xm)!*3|TfZoUewfBAX&(|65ovwc>;CM6u34mD>A z{t%OaZp4#3n5dm%$A{AHh`o%---GRy*49*NdS%W)D;v*>>mhe16be^D*4rV-fjxf@ z^WDNL4|F>(kDq?4F;)j4ge*NYn;of=XQY%&p7v(MCVXfAo>T;%NHVrPeeU6Xw#MPZ1&!&%LLkbkR|H~H4cY2N?- zvu{fhr10H5^e^hJQp6|V=}`B32khe+FHLpwdE(>Kp?<42a5;iI#gmn2^g8FB=oXPK zG)LibvG{>}iPWp)(9FT{R$3;jV2o~59@>&cIyrp=YA1AC?YUq)5Y`i{d2-kutv2RA`L?TTVrp3`TBWk{>p3zWiGqRx{MPPz@b0?Xo=D+EjVlxZ zThbSAmLscN*`kwO+%yn;ja|5Qp30ltWB@lFr$)l!98C zO~p|ev0iU)g1I7}AO8uJ;a1rP7Z@|0RZbE?Vz{*!De@_qtdp3x5JE0`zCQt7pezH+qFhyG(?*F7r4pGKPlT->LwR? z(2Y4Gq$h+fd| z)BL#j!8uC0f7_K!iIAnfcwE(oEAZ zJ+2q@=k(;eSvwLT`n_aEWP^g66gKZ?pbalS?m?d>3b*K4Lw%e+>qm3WpTu!vOvl&< z5qpi)g)zh=3zndw^!M?$ja~eKoSd%QUaP(X&;`@m8FWG(nQ?=39?aM)O5GhZS!qru zX;A%h5LMlglBU%)HT}KX>ufdck8^W%Zk(Cl4oW1%{s(%+0Bw8V}}S9US-^a z4<~L|e@x)pncz7(>C5YOiKIkoeQ@#2>lzm^U9r<>d&|4dMwk;QFDf>YSk}BEg(Ur- zQYn-njOcRNxfy-y-&*)b7lIW12jc8oh^&XvMn>Q1x;+5}gul3DsE1=~Z#qCy$4@Q0 z!?`=@6vw1b>j(lYq4bhp`zYYO-Yx6)>%nxMRiVW3OLS%TXAMk32dMQ>nQMLTT8xGK zq*4YTy`-9qgaD>}^o!e&({%y&_Q{y9I?BfUZ%|5BUljGAnMfq4PK(0iN$mu^q9tZD{JB&Yz?k) zR;ETG9M#moh%ii<*-f6_`4DQa{jh1lMm(p?$csR87|)to`d=Pa>aGRzeQ#!Qx#Ac2q7TvKihvt$OL=#l?keBM8Hz za>|ZIdRsbXn-AZ3HdLR@He&HpG!_-D`V5}NM(8McHx}rA z$IssmV^l?A)g5G+ev_d<_P>5ZOfiqDsH(vTR@OCd zdC2Cm)B$?-cu)-oM>*NWrTNzr5aPR|Bqnp}rJ0ttiDxp!Jnwlj#|P)7*RE(ojfGG; zc{A?>c%)0VjgdGLOMoTXhiT27fo-?r@z>@2y9tVXtadG85V#KB|8JIJvH8jUits?e( zsdU)x?lo~0B4p2F0>FC7OHRLKo+V28y-l2*{-;nw<$sRHkCdF~d1Nft#uaxd#;ju- zZci8rBJTY~hVmN{*kYyB9M>KdF^Z;sy)c}6!^kFdH2z3N#=2wvwe2rzStl{M-wWLW zIJVIhnPTx8p`)5k?DwI|+sQ2gX_08{Aq(or5N1Cx%#=?%QSZ8cC<}U#( zVpENOrT3|zfMNsJ6ufm<-(P^&^z;djp+4 z>v|<`rNk!AMpk|%H!w6Yu@J^NYvsO416H{j#BMB%RQ_eHQs>XTmk7u!Li(X)XGBbY zOv9x$oK*kmy;kYbF6~N$(#Lxr!7CWK zpCm^_kd}aleyn^1pFp|Awt7vEhCl=R##d4D9IU|?VOb@5&3oTB)bg&BvsbA29ggVh zm^9cHK^=Io-s$Sn(2sSp2CJFj1!9ce1q3q5_ zbCakxEN7j{F-_m6d(Zo?#Luw;9)tDZ@cn5JX-gCxY4DpBBHgT-1yqIqDq0&jiq*+* zYK^VkUh8wVT~*G;+aCois)d|GLhLtx0eeOJ*6nh(i9CVH zwxb`57cH_}avth4QkHFJCb9J>_}VE5Qd57psHJr*8#t}y6fHZ7{+QjC;K!2x6tUXT}qv5A+(kslUA& zZ7T=MFoRMqZgS;s^GU$hf6BZ%_Vs7VVzs%Pqm<>x`H`_KmMy(QciXw1UYiYntPBG2weYRHkJ_#b!0*D|c#e{?l1I+J=ml52o1U z-vUX?PZ5w(g!E{VNo7L5eKn=^@-png7g_aC_OcUMTT{!XMuDjns!2{3A7Y}9U3|V- z2@8$jt6El22vu^elmePHm;lysV%10c6NBbyh*ev~LUK%;C8WuRFX3aR3Sw;Mi8dSH zzmZ>>i!a!TzBJR;CIBo#XFt*}n7!r9{^tGwE}3hN63=6$93=|cIv=#N^~twmLIi8L z>yHSU>r*?6uM&Z$2WZ;_`u;=#jok%l&Sbucr?UByLERNIhgY)y>fdeY{wz1bl4p0Gww@)~m7sjo}N^ z&YL8KHKn6?WcEfKp7NuV=Zzp1#b+@|TlfIX97WF8B^y34s&cwVXKXF+$r zEu5T{|0B7o?NdPtFqy+6Z}aD*zIEvfgC^=AZE_mBCahFRolG69O?9&kuDib_26fdC z4PWEVo@Lz{8s!Iz1pdhd7K+0)7R2T`G~>%%Up2D**e}T+h4VurHq*4gtn7JPkO8Z2 z^P*IJZ8`^Kcb+Sasdco~tJ}q9d)pLi zYwtu~Zu`%lMnh~AG@OHtHjWh_o2;#hsA(NbNMp{#Ga{`gcDznfZ*;S_BTrW7N$!n=cnHbqKCej8iXfX6;$jsACi^#Tjv4-Wa4>d?{tCep2-PgFEnnZnH&^7c!!E?H)|Ydt!{8oT6Qdqk_u`m{mH$IqMqwwe=qc5T zJ0BO^`JBp~(!NH2i%Aw9XaH1Ea6v&);aA-^Z|HQNhc(}vNc|7+l|fCTa*8B~MF0IV z(JA&b`~k1<@U(?UlLX!*DrtZ9BKlM!!$V^(Y8$oFIA@s(KESB_B>bwm9~j(y)%gj* z{THPU;q406Zgs|>nl@DWo0#aQ4}e>fUj$h@lV9{#gX(mh`yz)T{sC?!;qaA`(?lzq zTZ?n#1;T@a0Iul$>BxC*f2%wvWlV>vuP?~5XYQOLHH#w&6 z)xkxuaz0kksi?jzMp%m-9Ij0*SJ|IL?CwG*AE-<3CC(;MplTxXBZ0(_=Sqek5>XsQ zhN`-$R81{U;xdeZGGFx(Ab^!Fa>$;YHoq0#r0pf$h1m_yDpl*l(o|^nNQjBKuWuWP zp?*pYkWI9cT1r6LmcDraXe@wJ)lfRw0jlt z19=)hqrY_}7FO|~s&N}{@wcqF6sA`yVUKGt%qxi_Rz%zJ+l9S&e};<*~_W@_cwalm`KdmV|I3`D0DG4Thd#eTmGPx_Ty+d z>Hm)pCdf$#pR~0HjY^0~Qovcwufc0( z*=;DnM4Pt92&S0uJC&=Zj1`<#QGo4pS9Z5ii5i~5@x`lWNpkHF`Yjo1;?M>M{uS3B zME6uk0ozLeL4GPR(TLjSSaK>!?C1>6`JG&`?w-emr2EkE{-3J5=uhD=yT3V-|{e*hpO{j*$W{Ub4uIk+)LKMs2yo(0Y<1!2NgnEpu0;& z^X3j)J{1ZlmbB|ztKs`x6Vs+J_9I+K%$Z|q{mkbZZ@xPx;jhIHz^Md%3GPkEgVeSH z`OH>@6yPX)IksPZFlh1{0NGfIi5wm@t3c|o-O=j|*nw6}46MDC4M@a0-?tlTI%eVH z*%+5r&}Q7P2mje0+;@LZY5$1l7z1`LvHF4ZT-&3EKVGK8QAbEnN$9=O1(rfaMpuDi zmjnYK&jiuFHV58Y_vXKqWxKkJ;Fa0*Ld?LT=cwx8;r-c5ZgA`LUWo=}eI@mfizGhZ zxj*GX-A}ML&`@+vwm+@|VizeOahu*AaeLTp{Svx~)Jxay2N@5%Wnq}^QBq=I$h&r! zC#RBO1AKn?V_{QQGyfR&$s6#PGCt zvlB90=SgtOC>$G&0gofASM!}Tjg8+rD_R91m5tbrS!9?UNA-U^b%oDuKozMO&VMA- zS9Vk9{);Fb!j9UP7Ybc4BE!4MG~ykKT&1Re6!x)v2$Es2A%uNx7x}Y)Ra1cqURrSr zeI0LHB@<^a;4Z-Ko?OL_O^=;b>!jdN8Nq*+o@S$Mz`#ty@x0^*Vj$;gIzWl!p&=>w zXGdUBZ<#n}gho4(){z@x_ay#%yiRYzVOY?(D<5u9+O}P4zDV-MCoG)49dT+6j;o_Iv1uN+dzVGt8Tuc3MA&*L+2@MlLd znzhX$Yw28HO)Za=)X4D2gk>EQh4h?UXD5a_~BcBK6jhP z%g6g;`WM{G4|V6xQpF$m38Ph5`LzMUf%v$Qh0nr(VV+Eft>^tzmvnHWy!Ed8?;ic@ zWHK_aNM9~d`V&^rXCPV_5iD_mgFuV1NMZCYS52b)9@Ec{)Z>?Li>DZN#lL)A@BajR zOkv!V-w(17C_WI~K>k?So$lbY5(8f&1GXa(%IKi(@V(;_19D7D&YwC$c%v=SUO?NF z1t48xt{nNNRp{MuWGqEv92_#XPSyaW2>OW#H(4X8f8 z9y52>tO-?l{`4|&ZyKL6g`Fd51nsswLYW&I7WRd*3yW+#^yvynwJ6hL*jFcnx>CVl21t9ZwaoI-HamYr}=rt7h){>gs-+rqtDi&6OZp!2ZA+ zV$HdAb)6A}(uo#@yU^!obmc4EldGto_c~8^{l?c#PbyWXk*c(;;%ZqtU=eQMy|eCs za5?9F@7^04HI~RGAKYloGBTu}AXN{%|G@1Lb#8szb?M@>cJGd^xsU_D>NuP2C2^i} zS|Q0WabP*@_Mu)#veT8(w=Hq9BtUCAeR_>M6O3K_(Hg!-(EH}0{o?XXj(r5$DOc#Q z22ok4v2ST{i^BJAPD#32K$7XAy81^&JrTIuJVJcm(PQ&9%)%oXGtM5z)CGCWPb#=N zi%8mpHMYlw6$l#jp>QHRU<;$2f9wCxGCUKVCSH;#NUN8vH5XhZ<8b)PBHR$_k|)$c z@7`Mzi8smE5y}Q;T!A3Mchlg!WJqMf#{CvxAYr)W5He@=!;@z|&K6c~<=b~oUU0DB z!c0b1v+40p5jJ938e-z%N~fter?Uz@PS(G3tf@ijAOQ?l-EM6?h+RUUO!(>ZYWH2OKXY%zb2{CL61)^* zs9rF;ap33ar5zRr!d}`7Q6q<2AAmemk%>jdW_OKtN6@(0|1J<)ThAbc@(m?#=icCK z9ob8(@W4&$_F++h-hCK+;$XGWCB4QO51pjMYX&5@)yN*{KF)>XZ6z>0eo)#rbTQv0 zINnN1T)&%vO#S3)eA-y+APa`YpU?>3lSTgN?XTu5Oy0fQw`Z0wout;6dAA~<#~3W| zeSKBiz8>fwXTJf!Kyt(Sr4C!pb2Wn~b&qBZWa!h1o zns_7nK@TQeCNRBG1-$0awj>?(_j;Gr^l33%Nx|Z2E+xORPY3j;lJSDd{e`rceq9(T z!rYgMcBukbV#2Rw{j{%ejD7P%cW<8ePU!meRIv9oPYAPiCSK+*)sS$Yg!DKaWy}4^ zkH~no`8jJB7Z?C?Jl(P*z*F>P=B?ZGxm2h;7SPLJ5j;I~qA!gvgz`g)Nmn<%Th{hM z&lROrJ~Rcj%_2+^uoZYy%X%2_2M=blS1MmnZ`iAJqm zC9VKB+8O*Tm-aO?^X7&(tb@I@eU3M7D0l}G4TE&Ml`qE=`qKV)ej-t58DTmzr z5L%qG6tPM;EQ?5eRa3V!k1RadQu}QoQs~DWZzhFcbJOk9gfCuEIlzZ($lmt{lKaNKp_yO+_=A8UEN=KkJ zF$_L=kXG;?GtkVMiPB8|5uW?I^Qo_b1QZ8^vjeYBXuNxWmeZcUdS5q!%F$k^Q;_%gMI`sqIy=Kd28aJclGRvZ+U%~! zvJYK!pTW6BMNIBNxf&_5hn13i@cy2C`l`sI`>HV`H1|BthNUq-wh$Kw93ym=`J}Pe z^wNWWH3bh3Pxxc68K5@8dBe^HddUTsoD$YR=A&HSW-Rg_L3IbO-edLeDo&FBW8r}C8Q8NjAl(gp8<$dtj#V2Ngg=i;&D zUFdE{B*S`5S&|ioaI`Bp!P1(Z5!ELBc>sPu-Urvq9n(3aMB73hxldUfGiNU#^>$&d z-p-1O^-%yE?BT~A^0|59`01TxQBmWmC?9l`M%{vx_2kEu(@EZysxUvunCD&?X4(CJ z&sE|a3MM5ka2*V>B_o3yfIesB(d6A``=d0Z*o?0n*~8}YUxs5X zj{ZpaKq2KZv^F# zZ($-=TDtUDZ?Z9OWNGcJFgUJ*QpLdJ)DQ@7${{P)H!m2?SeXcWqY4HXQixxy0aonK zH@DCD?`?2fWtJ5;!9KX^$0Fg(y#5~3{&t^cm#9?HhsjG18~KmMb&Km~l*ZKZB7_CS zS?p!)C!@%YOBNO&l8Go^HUIb74PPYXLMTH^tEtatb>$>aaNG*91+s~29hvvi`P}mI z{pq!Hu-syZEG2qf$@w`M3^99uG=?h3JQ#(`?w_TXwtUU9?pjMc(R{#QM&z=xPZ)}0 zXc@=_Ukuj%GB_C1Z+`vjw;|2iMQoWgnRk3Xe#W zFYJcv1BMqJVx+y9u>SiMT}sFv#>9>$!f!t$CRScvj69-2)QWG`97Mr?9}w$riHI*5`F;Xb834o4%&o2 zN-AB9eAm6C{L^^=PXInMTgWOAL>OUO)}Y2wsxf0T6`+->C zkT{UhXTYn#HPe52v+N!JG2EvnTE`MqgQef6NefZ`f$F=YYck{BPuaFc*dikaPgz}9 z_fkl>E#F20>@$a#B-E3iN*8SJngQtQdc~ z1VzAl=EN^N%mfrlcVaiKQ62lh_MJ@+cx^XphC}M|Bikux+k9mjbJ>Sux{U-S^eR;O zDf6`pIKG5Z45o|0xptTKJ-4!0(?Y>hb~O6quEb!I~U*OmoJ$lS0 zi_}wg6^to>(Tp7(J2g6RU$~syk-bZ@Tz3d_o8J-HI7_0rwe35C) zyovqorBE2vgXeUKL)xq zNV)$RJN;-ym6OJLx~D9MTb4}^K+v4cdnb4)F>eo~qIB~wgEu?dZ<=e8w3ejh^>h~w zJQL1eMJ~~CIlob7NKP^yWxw<>mN3|KxhY~ejoGF6b^B0&QguqY0E~DxPt80IZ9I9 zC6eS@S&87O;4p0F*=>!O7#pnIWg-f5ePDwl>Drw|>5#~k;0eIreplT9hsjkA zj`tKKF`(f^0}2q@$_EL0&#<$zy zkkhlvL0FEKhh~1Zo(bCo$OpLdFLmZ#!YCfeL+izfiy=?o@7$7s#Sa#DrB=NNAH-hY z4IgAjm}PeE{^?KuHAM9oV}w9}GA^ zbY(cPo;|jo{eaxYz2X}my`1pBi=7y_eA*#I9Q~Q%^BdsGVbwIKqDkbXj3Z3k zF8x`50@|{|Z{$oPt!>uf@JlAn&a2Kaz=98fbmGpmSx2%Om;mrSV)NP(*>opxURA(H z^RU4a=i77s$I(VVo$qF!+$nw<4KOG#u1}M$PtQ+}3sWkeEhy;+boTpXj&m>+S0!{X z_DuitnLKzZlP)j97<4R&wGZp-FKV=o#6jnEb&c1D@0ms@!`OMjlY-}vAz!BS zIo_Qs5xdWaCYG8Iri8igA^6}W`CBS^95|6ewr#?AE8mNv2G_0?cu=~<|N6x!pF6oP zm3_p0UqO0&_b#W%MmDcp>$-0~D*ory8FfsDb6XSGdw-A`&}WZ&uZH~J|3?y=>?JrL zy0ESL)A}68VC=|$9ZBd`X5zbtzPR0K7<%#V`!fUXOv3j&>FjrF#7IlM^P#R$#Cmvm z5C!GvJ`TF%hz+Jy+__k7%L8ET<;L|#F73%X016V@KSNBYQX#tZiodtK2%Cn8V=s}y7Z`3@Dh ziUh~=^C4y2o$%y;z8Fbc%frt9&nIk9BjH^^__wk?VVwKlGZPKNrvLfYt4z?U`{!)s p6NYpC{Z=Lzit_hsh`f_Wv^+bznacjvzlq?Frkb8=sj_v*{|Bt$X+NlAAjND2&!bcdh_ z?2G65y?eiVuf4y&J{K<58fWhNn)^EAIFIAFjCllCA;hD@gFqmJYO0F55C}RF0zuQq z!2++u(Cn^)U-!L~488Q+?7Vy}J!~ObmR<;FH!o*LD;6JH4^KxoS5YCD&>aC52QM## zr=+m3%l~|Xkei3S@Ru@jX>bu-gsPDz1VU_i`xmW5w%8E@329eTl+*V^?&Uo9qxgB& z-(EY&#WBpahQ0<5=TTNtQX-R+>v_oYma1Ht)TMp~-{SEjil^0-ZZGjD1Y@mqpe|e5 z55peDkqnKtun6(onI&Y!MXNDFpSRgK@6X1~uFBrE6*}0>9X~k~*tK1MvpW*FfjkdP zV2Lfj31|7|=MhfC3>pbk?)ImX5EuUM%kny^2(o`)jFk^h{?A(}$dLcO7FqBH#qsYI zC|*Fh|M@U9pN9j6`|rnZuB=}SU9lg}%1hw;J-mrf^tjO&wPXB9vDK45h zm<9uz==9=&E9OQxCjQ}iw%apCtbebzis{yR6xBc5NQsQVB1DJ$U6*L@M0G~kF76-y z9KK^aL^0p6p&k&>%(h7QczrxCp`Yow=8lMbEco-QeK@F0wLkh>`;*y>v(Tgs7Ma!2 z-?}@V4>zv(&rA!QqMz4YUy(0YTlW*S{$9`CYaS*4hBuaxK?CUkx40Ah*KfYzsi7gy z&M%Xu&4M_JdDoKTNfqYR*HZfvlCR!;2P=2DY5CS4>lHX~_1rt&h?vKp=nMHWQn>yB zw%2mI?#jJppGK9P71>o8x3=GQEv_ffdxWI_Uk@}Ai5eLnpOu7mY2e}E)e*72)7dt} zJzL?=iD*^XIf$LPB=+9iU^W5oV~}t~&ajf0$#c`Y{6F9HpY56czitm}r`o?;&V%>= z*N*=mm-LF~i8YRX2=3yg(D+Y)UBr&sb@SKSTlY76zl8ad#jmc?h!9lBIP)}!&$s~!R3B{ zcUdp?a6xb&$2ss3KHkpDI4Vb4CJ;rdShw-I8OyjcWBF=o>sNhXG5KYK>thKP5TylA z|G5{vqE6)Cs>8Umk21u<16GpwlI5k(W?%>QZ7fd^hO>0Iigiu>RDe#YZcSO%&estM+T1>e$O;;9OLRQKeQi^E0xXljy?QUI?iw*&_zihc4-%hli>KWQ?PUdS%+rZE|64RoF+OH|Et}j%O zCe`RF+}Z~BVc{%wJ*^&`>+g!t-@hkOi4T91OC(_U*4vPk@I2h}_wi320fF5^li=V7 zeev|gI*MSWucvUQzc!5`)aNQ%h`jZRSRAhdkqlxuTInVQR-xbxeS^szv6LOBTMzwU zk<dGC)0BS>C?#i>BUw1A3_&9#>^KPSzs^At1YVXm+`+`f#Kh#`;n_U= z7E8%5G9SlBR=x_NEM$77rr#CH{YqKC^%#$UTbpgKv^YbR6NjzQsF{LBhHhZU>ISZs zPr!}+&n|{cPgnH|sc)0Le#fEz;6{36@e{lI{M{EA>^SmzL`={IH|+h-=OdnjgQ25X zV&c!JspW1v1|u8p*#2FliA5kKPXb2Z?0 z3DLw(WMrKzx2bboQjYmbnMHhNCxO9#lNH%|s{2m%8rJoi>9Gjd53&3v+KYbKP+q;z zfdSF+;&)XlWXP<40^v$WuX|{#?Nu>8u}bf-u!4i5+@CU|hfW_oj=(NI$O|4CBRu~$ zb>z?6t+oTZ_fKK`(;cz(+_&y@DWGUS4Nzmj19)k>>VN`8=@)MIVcD>?)@kO`7jbD3rp(-Si z+XAL|CVs3UUL_mnkD~NERT0||9g2FQFqVteO`V4nwxmeNqw|)PPM@={s?WE1X+C~j z%q#x!n#Jp(i{RlPw?NB$s+O@j9^yL`as2LCJ}KiZ4(O5*&nB!K$&CkfxN{8w3aZsy zByPqXU!jrQnGzC39#prXMx_YNqu+l|U>j*qCw3B54%6d=y4niZPMsa)to>kE@eDD% zQ^nmjcL!@_!&%mTT%Bi5Af=mj&YR%$k{b*^f3S^(-$gQ5)y|9XN>KmbiY&Eww$=^_ zuDnzrgVPov$eYgN+4k>a%)uQOk3UwGf56*udUpHk$8D>T+7DfFO?L@_$_<3_)Y#=` zyP}hkw**&;xxAcle};t%C-S%b&Hpoz2QjInB=Me_B|oH&_~U;LH1mc3dLQ^v@6C<1q^g zc)CkjIbin{D3mf7i2RKI2tZ!?z)JGNa7;NgILy~ed?R>bH|M{6`RQ(705NEiv5(zrNt}?r{>5L-8r;s%K74_5b3CggnwiZx_#G;~o{lDFd zik@;;A!$9FROu+y>sa!AenU&=R>}+RZ3A|e=M>T3Flo_o7())a4GT;jYs%M`fbKpS zy$plzfBIp~0f7E<{A zI|3s#ad6Pk{86YsbKt`VNN6YiK`X}T9{%6$Cv8EmZ!g<-fc&R>aY4V5~8 znuX-xYKP<8d&jGQhY!a=&(hPqg-qHM)EG}2E|vv2N6)(oM6L+zZXS9vh?&owI*qu_<}C*6iYLAH>5k^KgBf++_{k) z2|2v@(a{+J8P(^Ok6$qMsn+Q%D7h~fj|-=>S0h=dwZs_K{4*u6h6f@Hgtzn$HN^Tl z`V(`))>AT5f|v<;IgNK?FkL-MhiCGJL7(|9*>MW+LG2f6A)t|lRqcU?CaQ`!9@?53 zp&HCVK@HP9&s{D0iM)TX;RxysoGvkrmX<;u+IojiWPRX6)sFvl`#1VP6MQl>Z1hWa z`_aKjqk9J?HE0TQk$8M|uCs7$Lfx} z%S9S-q-c%ej+RbUlGZm|P%Tp*Ei7-(=A3DJA)A^C9qfUfIG;Oysc zdSO*nRk!Ew-@j?VWT8x2s;Gubo^3iB8~+l*WU@oYKLP_gNGFG7rl3&#j zWRX6gTOfUzL;@S77d;u>wo1cQ{cHe+F>La%jo#qM}GX ztZVqt04l%h(^7MXG%~FOHltf` z{~o2h2L3CR*Cn@~?Z|qu988}9Bi&ql(zb_(Zp7{g#f39DI@>sX1$S3d>oE1Lu95U! zChd%kwHD+Yh)ZS>`@T~d_WAkG+X_^cCg#gj+ygbc@j757E%moNX!NH%wZ`oF5#VcHe92$VEhA6yb#tYZ6H zv*QR(@^o)G(9=4-izAogiS_#q6h7(g?YBNaI_E6@(^Ri75z zd^Len*Pe|nRE62LfsB%}CK6oZ%XCk(LA?Xr_9C^F9}E}H$f`T_PN-_$IU^-I4$OAd zO?4Db+dI{?5BG;s!D}Cpo1o9$?vRAK`kXaOIJpIW#q(b_Q_(o`2M!j0mf?mp3~Fj3 zNvO3fqRg68uUBJ3U7}OP{Mo$h=-*T4YmZ8L!>i%kkqsMA$p4xI`9Gbk*$g2*00g1<*5*UL0Bzs8Qfe?MT_YB;Y+OG+X& zQ{GKeJhkTpF1h>Yi=HGW=$?E4FiO0OVwKyx{=9T2ND_sBxV6zRae4&-l9N=`0)^CTwmZ?Zw7ZNP|V4yPy;S86}?p@IY{Wh9xJ5HA}v2 zGHMDUVPxqXBrh0fp8mF(!>6KdV_lo!ewX@+haC(}eE)cY6_u>Lwj?B1@u#OkcVR3v+JXb(BE@|e4|MZ&~) z$s4wApMLWE;B$1GaNAtvcoc_+5u3N&0Exw#D(Zc@VMFr%DemWd_o&)?517VEc#Qi8 zc}3U9dfy2?W!YazsJ5{6BGq(mOUsG>s%|&*GF0SzH|9P!93pBjANb4Y?Czx4vaWYC zCGz?Ces4z@#^hrEZkK&tP*h)^aytgTS_;{}T@*B_1{cjhvdCO8g9$Y$wp`((dw=kJ zmBO2ymlXblo}P565F5_LX*t&uV|ew=a|qn4qslO}P$kT%LdSPrK^~oi2~G>A+%iC* zcbH!vQ@^6SH&jeGg|?5)$+^q&=)1jykAR2+2}pWSg08?9B6p}&E+F9kl!1xnNk<>X zB*Q|)L}0lAFIQk|wY>CU5g-^C)bVY{!Z2bl)%@fHb8GrZ#i8NE7>H_EbX8!E7ba;| z*eCLS0u8ghFX0UNwg1?T8JbT~U%!A4uTlP$1_TyKtd2ZjHn9`p%6#Mj*hUA zr-{MZfAx&<49gdA;UKMV*V=N>reJ2qbocJv!_z74AH{jWy6hDbh3_#-N*;-W`R{xg zKYRX_y=yD!|4vFdE{k+>(7hnR+>YjoReyY*xaC}8Ph=?6-!pWK}|C+eyX*bx$txOm4_U0bzmfN`NeGI z#SgQ#w!Y)NtsnUNXZpoaKl9PWkA%XS8nTd!Zss>BDUpkA?Tjy8ywJ}gq~H`xH1WF~ zwT#VI_8&(#^n>chRBb47vsFyFBrA9EDrHqub6N~w_{cP~&n0>}Dhx#8-`J$+Rn@J= z6H#k<1#EjmTOG6~zFFe$@E&N5bxcj;ZHLoX!P6CP z|9cmUkdDLmRlb}W?&b0MX8xSeKV)Pmx05N_!@BNH?}n_6Ih+?0er50ZUr+0UcqZkY zbie++n>e$6H2!OE_h|cLeO%aWICs=eCosvx9v^$8q^Dy*y1G9fCvi&_)^U=}#fL9B zL~;*up;ytzZ^NG-{-@s+a}_^&fidNaJB09v?d_^TeL0VWa^E}E->@@1=C}QxW_E@l!ShxF)S1mG2oh{tA&a54`%27d0 z4df4QC1x>C>|C8AIhS_st9S1>zIl$z!2AwbnEcmBcNaeqn+N~qWD3}1?dk31rFA>h zHU7}bNft`$x$*t$uCucU&4OqSTPPtGhS~5~pY^y*u^&O(;n429-?UYyVR7<**73Sd z45z$Y(~Imc0tF#Ql7YD^`%i$y%7Z;j?uRRsw*}0&+(7!M2aI*vDuCX*KC292m6c`Q znJHILQo>vx&5p(wZhaJQ7?D0zJWp+g+3ENn5;jwk(tmQYk@YT_e`e;Ks5c)NKQU{` zJilRtW+8DZ04ZRG)Hm3`DW3X_DyzUUep(+FcR2Ez(Rql3I6N8jgu+G4bb1L_fKyHo zy_8Q-rNz1Oa&loGj7I)Xr2gy<58Q}W>O!9c2_et%q!NG-V7=U#Re9?+{$%n zK$t2G8!s@I#i69!Ef~wU@O1I;pk}Ppf=9lZb1p;gh(S@XvrT*{Ciebt5;>i}Y=l5B z^H2;=mVI}YEsbp9?cO5)-`Gn3<26^=MHe+0DXB|=>MaE9aNa&|wSLyDTZK7j`a00x z|Nd)6sdblckziQaNpZ@yDWIdJeW*u1D;j+^vo=jz8tx?|Y|9m%+ZaRW%k1@q_MW4J*# z_%uXkq54R-t0!+EO_!E-w|ewTtjg)k8uU$kav9?8%$L;lIx*w!L2@d55d_N)ENyct z`YIdB9kRPT+QTRuE_y%n@$c6Rlafhq`?))f%*_*U$+=kg`SG9s+U?w(ua#;I-pKK_ zT5NEE`5(V4lf#nXGh)km!2uf{nqK*oUMF$i)SvSh@(scw%pS4G#a-i`L`x)!w7DLm zP2x=Ax3nDDj{m7J|5(DaCPl#1~trPfJ9MA$=mNL9A7B3 zd(}ZAjJs8lQ{!Fj=2CWQ0qiCybJT8xRtzVb5GxSf;v|YNqt*h=*NU8lo_oHu)N-el zjI8eUyPWlDeVPQO7ZcO(ym_D~GX%PgZ#VzqCyegslj>mdw~h*Y-;+^w#^swU?Xp(Z z)@Z0QyZC_tRmd%!#@k9F_kvf!l@yxxbIk;SZecX1A!Nz+`p224%?}RI`uc}~t$&8K zn%^$b#q_-THT|_(b|1=}*?!W$HC6U*OyfSIMGUoZa6aAx+_AMao0Qbl{HCV)y`2`H zJ)q}oZZ!L4Z)Wu=`*cmqBhdelA)A^sIiBb{FvdR(w=|IPMcBk}On#gBref<&1Ve<% z=tk)=2cA)a4s0H-XB;d?3-3KrL%e2|Wf~nF%@B8?Gd4Dsw1Ejy;oFU5gpslDviDYma6b|xg9IU(ML}U{%zW;<_wL!N(B^7)bP;Kp99)MkDcA_2%OaGW z^Tj{+MH{=1?eKK$Ygh_}&;H?~M0PQYOCQvSkB!D`7@;+*e5+v4o!yz9)1i>u9@-?|6Abc-n=7Z@!R92l&?Hu^v{YvH)7epLrVfp zAqtWgh;9!9Ev;39 zKM-|!`N<=MWYQ(dbN^e@#m{u}z~#f*@dt5uuSIpfy-npziWGwcrmW zhh42^Hslq*$qoMHq!5j69y|#5=#IwtXi&<4gu~p?5wJ$V(&&5QY_HA_EJZvIS)QMA zz{trUg^5sZVbH19)*8*-KHlTy4W070tD()u0YDBDpXzO0-Sh_Gc%Ez?)ZX`H!LRUI z72=uSqn9?<&BH zem6{5qM*JMCpEJ?DovqQOtghnwi)lJep6UK_%PwJZ%Co^ch4=Pn>#P?`&5*Szc-oC z^la@tL1R34XowTXU<5fVIkCJ1BGmY>C8N2p;QNRNrnaII+Ce+9N#Gq3C#u74;tp(P zkwud5;p#N9z4ZeFez>V$e&1vy4;4Onic5s*H~IC0K+vO4ve(_OxF+0xV+y+KzH-I1 z&wCJ5ABs5ql){N-_UGt}SxSWbsj)G(b4yLGad z_{8%)L2&`Vtmg+SmC6kHxS*P|W6ZTYS)`G+tVQd8ox@wIqFBWVX%(I=DS;dxw}#wQ z2KOqrX(4e;{5cQDL(mJAd8zZ5O3H1i!!G>F+3(ZKS%`jB8*zy4{c?3DDG412qLN1> zxx9Q7N<>wJfN5fPxZ(@?j3?9VCO^Fkgz8foKZTiTW1U<;v{$j)7 zb{AP;V~=4T*a;u;_wV@L8+Z+$E-1oM)wC%{R$y4`EoBrG>=>7Qp1#UKnlIQ|xsn1H zkKWZ2O(s6oG*#op4>mR;KYSfW8MOTX{bLVHV6;vXTqIH%23aDcbg8;wB>m8kVR-25SvHxMX90`@R>F9BnHMYNu zPtGO<*ax3WS?m~~it|EQ6)YB;$j3GFUQZesHQtfA73&X2PPYI#gBCr)2`UR?X{iW| z0cWH7F3GN;K$+C%mc$OpJ0agrOY6O)dBhXd|Rb zI2NJ{_|(FEU}CQVQhi3mD7>_&=#n5Us?nIe{&er7yBKQ}fJI9l9Ir8pmt6agU*K1} z?HKgn-ik3%S{b(p_fL?(qO&*{wR#Jw@Yi8EfkpI9+vP}IcLv7A!r?7wn;D) z$hbM3I^K=rh6$HSnLxSo?_v}A`HGT}clUkw;DYV>NZBpBWA{}Wzi@f{*wDJ(3fLGt zZokqullI!gTYXPOkm{8Sd3I>uM$GmzgZ$OcmN*u)@Z#igvnZC%VgjHtE*=eX=b_pz z_qzcNks8Y@hUz6#V`GmPcT}*owO_<<%4$DxXf`$5zxks@ygOHg<|*AO4vY&Qoo3PT zsRo+0wcI5$b;SBPz3sikjKBYITzIqal8qes>D!m9=l65O>jGEhrHypM1JlpilSC*4 zsYXJ+s}iP6Oy z*4W^hg(Y4}~-k*iR^Y)oeF(x3>AKC?I@G(!Gogr=i508cYt%;Rd7y8phbm`?>i?X|X@_ zMuoY#v7?EZ`@Z+aO-@fow(K|K4-O7yNP6>uhbr|#o2i8WUq0D%;S19oWajZa)L>%cQotYWBgr5N3g&Kw~ZAdwE{Quv;ov+ z3z_pp>%fS!_>b$$gZ?MxK`(&$3G35Ju+)RUhYBXazkl9#|9?0IadAFUCq$=4bhDx| z)>F}004O=^5gVPnoleC8PW;0~117r6wJ--h?^7rD&b63|3N75)`yMn}Vds%Nb{r<` zK>9!`hO=`pg4>opLiu%%GG-*}lZ_Zgskw0zVIUZ`WG7p&vvm4j99l{PmV!pifh0}X z{Ppgx^GKd7^ZAl5^MgN|>?R(AoIF3ONnc=t3M75iHZ(NXoHZ~|YyRV2H?N?^oBQ4W zeG<*gY-Hud^@sdg`-xe}sp|)u_?SME=%wcQj1Q_4nOMWN&1%Hvuc3~Y{kJ#Kd@p5d z>PX&zY6d^Ntn0AsMDZElaY12qvAuTLa+~SWiJ}gv=MHL1a+cBNx`GFwx#FpqTk_<4 zPt4|QCG+*4*>*Mw32IO-Oh~{plt@{N_Dt`qEca z?yYv-gGS)~L+sAonGYwV{kr-sLa*RHFK$6+D9%SppYjM?%%(YS6rJ#lZR`2+bDPXHqGbSr{X|2A z`7xph+;6fzQgge7>xJzWJ>`N!bA!@s$tqglIYA#VH>?F0@6fGoE~ZO<{rbHF?N`YO zL2=~dvlr}~n3L-EXn6s#_YMO8@JhMyjP5U$B^4K+WkW72^q*jVtz*qC(*EuBQ8eJx zq0F$li@zEd@R^RPnwqElRd0*xQc`K9Z4c_gxjVPK?%vr_Ydw;5DJhWu6~&rBz%&wh zzEHottF7ggvudZ z{BGut$vRh{fm+p8z;{$_B8s>49&F5HXJ^Y?yFV%}@=ZyplH&Cp6*2Cbs>|<=HY%M< z=Fx`R736$^d*2g9hIiCZI?Wz`O1t;Sivu?J*HlS>s)1r~c8E55{lLL7th(8l5#SYp z4Qg}9(*~aQ$^y1Dou7e(zl;z_!t?Z7U!f+o0+vrnM%7{IyKBHuSawC?npWJy#m83w zWP&^B9VpQ86urE>9-;}Ew=n=(lQDOAe7wKW{RjF&MbY6*zEEOfClonVdKFW^PYe$& zH`0W!i+104ObWPcS;a!dON0hphQZcLsxa1aRT7 z)-`nHr=doMt%GwtH~3nh#Z8&1(|I9gO+h>_cIxswBkJRouu?K~KY9FAFEeSX>Gn)D zGm8!)L(b37H-aEJ8k0I;0jUe7vKRBt%8rh_4}Pzab8vH8e|e{N0jl{eDO^}mf&(~; z6I*KgTmS7VX4waCT@Wg2xQMa)#wzaOSAVMAlj`e7Obj)efiTZ+dXlm#?VprtHO*Iw z(eaYBOwDm`0Y!?`Rnrir9h=sVoSp4P9CIMimMzGaEUGVj719Bsl;jAJMq{zbKkV#2N>2MHZV`W7Hgx)Qq6uGzf z$$F|tt8JG}40XR7_pthPZP}~xme`o1zqxf{ZFiA)#%%r^cgbj^s0?eBKL_z?1<5>R zdz(yP-*S)4jdRfgjOOJb=k}&73AAGkehkmXXxopJf1T*vy)@MCQ5H83DD}8}gjblD z;|sg2fTC+p6$g?bKR&5c*4L*f(#-UjadLNe_q$k)u6CGYgLK^T(r6f1sGZ98H101J z`K2*SGDiY1f5ZNHFi1k(XyLq6Ac7(MA4AS4L!74nDsck|`W?YzsJ4V`=Rjw~t}@VA z2Pcg;UJj`Pu{ex|G=CXRhy}vY=+hriX=)P+*evEmJ{3xM%e?_4i&)SQpr%Ysmu$Dc zT8N;khK2b!IWYkha;-xBO4aJMga!xXWr96pPOA;9g&f(M*{&fgF<7G}czlxiqQ5xr zLyP(Qi9th+n6w^67&}{uyQ<^gYQ-c~HA^S_ot>ZYyb!iBgR`Nn3JQvnH2gt{1lj zc_ELVpI^Z($3#rOur5XBIAaxn+_J>hQN+StM-2IP&omO*N4rpF0bB1J zTbQLtM`DsZr}z76!?~lM0q)6sLH-LWh63{|#LGbze`(I4uT}P&y1>h1t zxL~QNsj;XwJp7^|j0W*JSdOo$xl6$Oz$CH5TSE<33iWj?++U+T7-&>Frw>4|=@|$M zKWRq_-En6-T)C;)(g-MSNGpH8{p?D}9;igRYNn2XjE-Be9J??wWR@)q4Yl8r>F3sB zHms~e+OqmftV)3xY7|wQ^W6$;>5Mc7VV(|c9<-PR$FMzE=K8c)c6B2!r0*LLH`kK2 zzkfbo?ta6nP!)9P0|Yue$eZM34A9BrBHnu%2&tJ8=%m01v74B@l8E>t_Uhif1zcb+3higYNE1#! z4>Z?@5-@JJX!B`vaMwC2Ftze~Yk=vI>*ed2XcRP*NhZKz40)f2msif(ntgkFdtz?x zU-OvTP=f9BY zPen~k_qaA-5m@Z_p0U8p1CNwW0e@FsUJeP}+;r&a@9%zFuhy755rHcEXq}s6yc&Pj8T3VFOE_J8nfg6|P?r@J6h%HU>D(VXGc%|x z|7segM}|~pYK!Z9+1thS4AUxUED`mpq%4x1SM7P_(1jnW6FaOoF6?3L2^gbo)2l~H za}!gb#+iX1w6gSh5$V_kTmbE5z-{Qc;wZg0UGLoz(>2EqatH}{X3r*4n=(on_kB2%aQ z +%LGjh}@|w*TTHWp~eUF;{LN=IM!ecQ))9~7JMxea#h%qO=-m5n;E0*c$r)P}H zhwc|5Pj8qu2~$4)Zixz?0!W}Mf$>N2JEE$Sy|uv+wyjR=pu98#{T}CeGP0ByZI`5E z;Fz73#PfG>v0YFYLiF|3IK9V;R|qb)zGbh0qDm(F_b>u^1eIv@Ezxke?=tX|7>;XJ zMCj-evo_Qd)3LwQk;rbW@?9!`-f2r73_-6RKi-9E4@6Jr?^j|MM@}Q&xl9o; z7QU4kX?_zn-@mkvu)Qi$6iJ=t=SPeF@d^bM+tEKaLKr-G@}nr*Ro=>q&BMb3xDXRQ zwfg|A=sTFilkauBmA_>QBXD7YT!MZ(oxH5TUb$43|1%|j3^+*z;;y%fZxJB63$#LiMdf2urR^u7P$Lp=6SazE7q^$@k59Y$?_Iy!tlKx#L{ ztQjuK`2rU%BVn!%n0++9H%j2~{MzcH-OZ3~(v2Jt8 z&lp!x$2`2Y$tat^$cXmo#DPl_uQW1;0c2G}T~!>7aCoh+8xRXl9TW8Rx_70#9Lm{j z%LeGw>;$NQO%&r32_C-nd#`xRVlWkmvY`t~FYC+v<~SS{?&xQJev`SkIRE486Mph! zDp2vuP?LqmDr7KFZj^_=GO7^Mj@@&bLEZKu@b|R}Fc)t4^0}M1_fk4IJ1aKL)HeFB zZR#{UNVM1Bz*U1e9z7S}x3-!}kJB*DYM>qaz?Vk2{vHHK&7bwr1^BCDY{xRrgB|1& zi)|SI9CmA)kJ3?ZX*obtvvVq4Vk!87(IQJDH)Qz|?+KXwa*xug2_r4;hRAxoG|gsX z8y^Ic*9I zef8jY$K^U_H*Qs_f|gcC$aSGNM_axEh_z}JVi*~@8e&L`1P{Is(2n#OMOBDwm1~ag z0dd1w_CrMcCjlWl(P`7cAuC#w)bYm;aN%I$Z5F~>V1^+e;Vf?@65gYiPbQ;Y#aVM| zea3;8kFNM<+?~!m8MAh(nJJr;BwLwSHPm==*y(gy|Toi;VaN&)(?yynT-X1x^rh?cZFQyTSpB=z;@Qm z?OK3M7f^9X23Q9ZSc%6y+W~TI8O*@noB+=FT(tp;!(U-Xh(lpe(rprss{ z$HPy66Qf6gSSIpTc_QbA_oN9vePm|#nmgAIo!ZQg)o1@(3~pCIt~4%#RFS61QfC~ty*(G)27NzHxbAh&b1xHJY*ob2>)tNc&~Yj4+`acj zH3u2ptmD4|d%(d2a46^+bK`6l1GKXmq^XX6k78?Xh0AhR=cro8^j9?T#rbf(tdC>0 z0G}wuQLUcMn)fe%eqQZ=V)qy1c!q=<)At`gbPNq&-m&iEAuYb1EJW!`bvxGAyx6`5 z2Z{FRqWz9hUYy*~sl{jd2I*N*Ua{esi|65&0N$kmH5P5h&^H=*NzWc10H7Z(WA zOy5tdgmg5JaIT_L3m~FauF&#x*B?FseL&Rx{o|%~xup{YkJtq6pFm><^Yh z=L3H^(c!?`kzOId-x+w+OAQOW5HXgz2Tn7F%3f`0*FdSo7#@t&NCXTc&^0wNC3VN1 zg6`S5>hxPc;1wui1E38(43fzMd|G+$_GEJzH3Uh8Qx?)`awMV_64o!HlePOIadg(v(`*Ia`)=?qgcaV3V zY4;mNJ1!|AV5`^C1#WLtASw*|5cER1(+N_|@F~IoU?+S7)i=s)4xciQ#wgV@b@4kz zdG*?$D<|0P!2GF^l1xtXa6cGW%QTDwTlT8wgMtSho!ms7{{(l0NWHk@IpZtwx$a=X9x_ zj)6hc%!~mnXN)LB(q0-DNs0{MD=jXnA&<%0fnks$je)|o**Tv{?`rXj)6{R%BMjJo zL9}O=!a2QZ-z_tXp7n4dYa@)LBv4dBuC|MW2xMWM9bIe3!7>U+A0Hj;tXdUKt* z(>ZHi%SqfR0DBL%!iAw`F_S)myE8L)EUoZ}4~U{b#GJC(iH8BOuB17#mNn?3&HuG)AM&~PTM`Hx zG9-v0EUt@vMU$ohr48=9F=8&Nm9%5#!IlZRP`#ci&T_`90>i^dBB)@$8vNUWEiEftHF{0{bVx z#LegX;c{Bcrgk0yyo*25p1OAKeOh_l8G+5dprvwj0CccolVqZxdxL)0>0z$t?&e>uoW|^#aKl8Ugkb6jKctr z1Bqg?_C4l*Kdqaa&Q;fdj!lfNZla;9I~VV$M@~Y?dGkw#1sA@QI^9`um26PiR|w}` z*VlgtVL|6HX{AYSaPHhTZL#pN|6J&M>KyO2N%Ayz`cvh`xs5H_^Yaf>Nr1itjC((S zWh*-5ma_q^z-Q(;RF0){E4I5+F~&?3nVP!0#!OYT+lCIrx2Y&#QO(V0&2KeXy81l~ zj#o%6)jR#BMBck6XIcH-q2>1-0|P$3{}Qfp1$aAGkYhfiRqkK;M!wZz1KN0^a#KLi zmXi3r`Zi?P#5?sWkrrTdgIe0L44SXMc&YypwL<>jTASvQLtc4HvKs#*k^yT>n87M< zbnfWN`>6biG&wul(&p)*;-o)-f7Q{^nf#gOQC~cD2v@>Q)7Yqi}DXY3(RI5^$zu&-ho+8fgx2YC5Q?h62I9SL$oJu*pkq4fQuG zT;gb4=>O|q{1mg~QyJ=}@WAsyJ^924o4+%4?70a0fQ&uF$yP@H6CL@V&kk(&xJ0gsC6?YuQxKy*ESa2NRT8~W-`-{z#Atx-i3v8|P<`;2F zE4afG)!UV&n?m9Mya0?9D-=qkT+tZ0V?&^3a^k3O~dx0N&>Fl{l zd7@C;+{e)hlO{4?I|GMm->1)-wsR5QXCQK3sKb@LIwjZD)x9-ip+km-hryX=TSo8u z|I(0V4SS_e()0^>`xGfqR^Wy#F)zxJqg5{O^WcQ)P3+o51KnT(FgB^Y}KXOdq(4xhQBH&&qOQfJc z0)Ti&nj!bk7bOTt$v>26>PBp63!3yEl4laO11~Jex!z86E9S>zVMJ*oOFR_$5|g@S z$Ug1`lkxogS$19^)tQ;n(|q)JTmM`QF*7@B4G8Ojp`or^aJnA2(TFNn{5Nvu8eHh% zxzfB&Z<>`W(^NY$N6P4SJcftRsa-)uO;4|-tqT|$T0&q{x#fK@pfMvQuog6R5H96+ zcr50yrs{mkf-X)&(0H-nM7j8>;AKG#=oY4ksM+a`i7y3<_mH>V^^f!Oa~2|XD{C(; zR%-w2OgC+=hwA!4z|;P%cuum8ilfZTK)f;S2fv|!hcx`;${0mJ8(soWBcld@g@FkB z9#IqUe^w@@UG#stCLCW+kVIe)4nPoS zY2io(=vN2*bBqu$#4;0d&@u^{L#vzXM$ZMY-e2WBDz9%(6#eMr&96Z@pQEPpImsXy zqjtB)#-A6#mCWq8Tmjbue*MzSp=(W34Yk_|!YFFYcwE=2_Pu1YcWa6c{_$(<_0=W2 zxHuJLjBi#^SsC5Jhf$H9lk=fj?u^cbwlY@L^G`LGIMYD?mdPo{EUUyOBSTO4D8c*D4wdBeZh{trIh}6P zcIXj6+t3D=A5NyaT7=}?G64>ZF#?TQ?fgoWXFG|&g%Etc{n4=1o!xAq2J^=EXR05U zfM&ZDxkb3Bw24u|ND9TnRZijx8h-vsGG^k*(YF*bWMtr|yzP|Q@vNH_ojMB0fOgF__bNH{myia$0XyUO>|rVVrpHdIf>oy{hD;x51D%tK0` z)}EVmIOw}F18>z*Mb^#RTp5@`?wUpGCqDJLC+HK_iActem~+R4)n9|vKm5gY~O@vvVXw~2>mykp@z3izSy+IHDGXWAMz=RYG1V5JBiei0u zR_{-0@Udcei#2|EWJQVe;wSYAl>6yX;{Jc2JiR!CDXqT*UX{_-GJ#Bs z!r^R-D|p%SUSkqwgOern`F|Al-SJfS@BjARdt|Td&9O2g;S?1iQg*~aMpR@|G9x;a zo$OuqiDZxL3P)uNg^Y~fb-KUz=ll8n{`uXH9*^$3$K7$xdB5M+>$;xTbJVmOf#$s8 zx)OlWOsse& zTDq86==yKFYDNTW-6 zc^;pdIfwu3=C02KO}9L_*~j**Wl!gY<4z=rpL}5rHl^&!Ei&>b!p@aj(}Ac+br=!O zr=BAu3D`Y;gIIXmUyK!O(%VQqhbL}YSg>M%dOiPM^hkaCsu}6z<(LOYd$JE{-%%BU zToqoR&emK%P}(N8V9i}&VWHu^&4l9)Xi4D69dWqHYRODgSlPf4I%3FyHjwPi#LUkXa6HPeLWy zi^3JnHu8S{JzF<= z|Ilajc-?uy-tN zo+LF?Rxtn(^yTmR)`+EizJ~RlmBI_vkMcjUBO}sg*-Fi=T8o@>axKHsx#iH*)LJ%~ z!R$z_pH&UlLw&lFgu43s&+6z9cr+{*{9b$!w64+wcuN&=p&J0Ztcq+>IcCBaVmkR$ z)L^$DL`iF*)b0u#!}wKAkofHHuPxSkWctaTV;`_Xn5*)o?FW4g+Q?QECi(E-$gNy6 zx3TIVTel>yumJ1s;gmee1UOhejc;=sGwsMHeT910MbQNPtW!v`z&Znpa$-YTTi!SB_`qn)>qnL zrP&N)JmbYqj-uM9Ev18-nTf&7@A)`t<&=gaVo>EnYmXRkFPP4< zU0XwhBTtG}R&D{d0h_Qc#!OrK_gE&hcllh&>zn=O&O0EUG66wmjp)1M7Rj!IQ&qRF z+ql-~JEAku@@O-s-r1_acLqG#CDxFJ$evD zB`Ru+`XwWnD}J~P_^1iW)^ExYc5WEu*?WXocg$O0nQyO*Q7#tWP3sYW{e7S6BvIiP zcH~!Y#VvYH5Rzlu60?%H23sr zNg|RDPEsar z0i#f6+Q zOh9EyT6*9O(A;iNQ@?8bvmT-|bk?%Dgv#ZX*>#lI+I#jJwMnGTGA~KDb{CxdfC<~y zI_T_JDHvy_Yl}5VI{;AyPM%}8GFpV=kJa7&966k8@U0!bpl|V-Pb{mQMCl7#5{L>N z6zQlyR$UVEy{TPXobc?o$|))Hi)sje*Uzh8Tzmg3Q?~uK%uH2H$`uNNZWKXmD(=1} zCnw)qs~)TfK6%jeXT321NNqncq&zR1nodbKZH6+)uEc*db5b5DPV0Qw6)jlt{TgSH zUX&`z=~g1cDk2jp>=34hTEALPW3T>L9*j$HOjJigAikgeXxEQ=r@6H|Ie11_&sXu* z##@6_wl~u7wDr|r{n7LAU8sD{>sXG+P<#ys1TELi_+yHJ$5i7~#V@hjK{pKj{FJK3 zS~$ljG2=2scNBtKH+AmTS=yF8#;$#~P}$o3WMy7P8=BPc507$m@H8IqDe2DsEZO-~ zMW6RiBd{gexO9jA=d}N@jJF{2MNpZ9tim;aF0+M!vwguQkL;CaC3=+{=iBhdv&(*a zm!vQ$X_wYMy}Q>FJ@?%(oJUPRRGv*H{n<3#Ap5}An+uL>8@v534&5ACa%j8?{|(zd zv*ofFy`70tY4N4Oq-L5O}STCc`l1J@9#)?aQTDh7{gYMwEr47^AE;i*N>b7mT!vWF_5k?4%VM1KNFhm^# z`6pm|*x=`oE^T&00|%|f_Qj});G@;+_8T+JuPg$TBqwU>SBGSR}$F_fnTs;3P0J#eHJYg~Is067X4wqGJD%3F!U zD*X5qv2^suRK(uzQs`%ooljOugV%~-UXMk&D)o<)Vj35yX#D7#OF?Ck|K z{r=<@9UYA*9`zegS=3@WdzHB+9h>!9OsQ6gjkpj!eh!ojX8Eo{@@RJUV-s&I9c*33 z8vSkzjuv(G_H~qh{hB~(H?lKfKSRpX<{&4^36@4aVq`3g_myip-)Ad|^%Ns_W;U9_ z3Mi$t>^;g3y*91PCtkxjmF^p+olfQKDAel_*OUs zIg@x1N}8Plhu6&A#ctN@nbVSv@#o1!CA+)qgu5jrA9@y-P({a@l__Lql|6>?br;p3 z8y$^nQn$Y}TW4~S{KPGSMm{!37VX+LFcr$nTP|Dky^(<0VC=51=TaRU(Kj3C5W{pW zFHeS38Qhr84ZzT>(IbRl7AyWEh!8=dqr=0AzT6oxnQZj;{PSi0QX%rj^2hwYdeiP? z5dl?It3*3rkjqxjlXt$*2I6q8*auaXq8sOmKGvi}qy%v})0wdU3ctT}5$A*zeyK2J z)-vIFixH4$>ii`e1%55Qps)57-EH@%haq14cfJ4+*tJ5bYcsbmNoJPAeB_x$|a)%0uN9pDNangksvt zdcmpZMck}(+K+H~CF;-WoB zoU7IiltL|0=hK^idD~4^wAEcPC!u0z2r$k9#iC+PSNdFyRk|ISH|B0&8ycbw8Dqt80_<9J%wq7I)inm9VNl#7VzFC3O?Pvc02?;P}~-YmHkcB+C;0l5X;2WN(SWsZiZvwMgu8ht@FcNLcgT8 zRk;xqyXY(p{Oa^ss={k6&6bv0tc3T#g!BJ~^E^HQz>7N$4wo+k>pbfdXbX*8a1us^ zDz=k%*~g7ecHP@e@hCZeAx*iRT`ZKbq=b9tTRfSXm(JTv7-j8+iR92`1;svvfdMvL zUynCVP5(F#T2{`Rpk3=uW@2D!Bs(o=PYC@pXAB?koaQYXN5 zSv~_Qugc3kJz4lwO<@yjD2Qm(cZE&)+x}T`v1iVC9Ddpy_JiFId1Gbn@8^4~kv{}L zADK)p>E-3@`&y*wmr_JCI(+!0X+?b&`Ab^3OEU+)y}nV2XW!{k?HLmePBEAC`EJx$F{iVfbuspq z+uWJ3G8m>8B3_(N>#gorf;N?znHh&ig1Z{$y-eleUmK|}2Z0{S-*lfLcC~x_G<7s|qx1Ns1 zU_)UsbWh!F>1kLv(Kv6uR^&)w%>{6b0Fm_g{ENU+)@(OJ8Ec@4u(xPcq0a*xbxug~ z?YLLuMMEyUH49sEQ@RxK;)&6UWb0@JJ|(cG2-uUsqQG9{ME|NT;mUf8b#;J(bTKPy zwKjuQ+9(M z`DCkuye>hN2o`i!_V!+4YSwRM)A;0Ia7aug!7PVY<`Bey{|Tf0hg<6-0dQ>!q+5mx ztZkt7$zj;FxzQNP_=OE`N~D$==B>{n64wB6Wx2V!lK(I+AnmyRfcya_9o;v6eTA=m zsozpNj$0hL;G7K>B-8`VdgOyhKg9A~0!e?HqrEC^9ahi|#*Zdj_ybOIl zE8-rAkT-}k}Hbfa<%vVyW3#-9+< ztzRgZaP_6F!ofl)Dy}f45gi^;X;#|`6r1$S$uT-KNqc%{wpQYCYql`L$J?%?SU0am z#s;a=u{{lp)&(wq-V^B5+1nZ}6;`lo9e2h(yZeDFBNi#yyYe+&vuaJ4c3fMI0%?DR z8GuV*RjtHaBq5flA;tfdr_r4qPSMFFnyJNgM^SbVBbVSPSp(MEr4B zdl}U6>NNV4`qgw_G+)x(mwI3E&cH6F?`oNo^}fXo_`W)&Re@2Q+pon+ybu*{D-+{6 z@32oSB$7HOv|lQ0_=2k(cfFFZR3T+jf;brkzk;lFMl2>}C*sGw1s0 zqL=JYZq(sZ?=u_J^HDsn>s8C=e{nR=4 zLwxndg#>2=v9=bNUh1{|c&mf+9WL8f(V=GCyKi0p)zjw7%hm=z$qKuN{D0P240%ab z-cNm^fV;zzpIBa0^d9+u&p@JcMU1#c5J61LCo-fdQ({`$#M~_2=IZs6gOXkH&=8;p zuOh#1oYC1pEr;R}e3|?+hWQQ{^jYc8?Ye*s$As@g-)Ke{%KgjQ^tEEE| zl2EllC+sloS~2v1PgF|(4E#3bB7?16b1!+Yt`aA3Ia;f4wznswS1MUKt1Ln9=Gll7 zpBwR~@J~HSkSroBDalgs0MsSiXI)7EcAw#CbsiD0shM`^c$p|5 z>-kCHnXz!GF40SO%jU&_%r^^@vhG9)8w?Hy0)Y z_^W3P$zUZ}FJ~{8zPU1c9O=8iACPY${meudt-vAc{#6h6HtScd!$Pw(rq;}_6Bmwj zhrh^RAL!5n?NtSGhfeXn|NE?DLR}+DV~{*`-e@P7wKJjKx)b+^@}9LSOgwYPmenDa z)o*Mjrftnn27^&+%B_jcxd8h6&uGs^Q2uOtYLWA@faIp7Wt6eBGCUxT2rhDRy4)NG zidZa-#g@NhvKKS61bt2(Hs$vXg`O%y1~h(BlaI(kc*I;=h9cr*rTA@A!s$@W%oT$t z6+0<_5_TVIBp`j4PuO1z`QjB;S7Db%NDHPL=jqY4wQEjV>+jc9LneImf-LWk2wOrQ zQAzJw?JVHb%QE`KpJ#>V9-Z_H)v>JQAH?l$LpoSk-aVLq>Jk^f_**B3?3fH$3)2@trqX5G0S*eLm6BgBG#A18%Xc6pHxCJCAg%N_ z$w!Z-s2j~DE&<6doyd79pVdJSI_84%mUl6-_BAwT?;oZs<$(d(1zL%A(YvT$cf3aN zf{-*h5-j0dYAk!5a4Oh%Dm><<#*5w`3_hIEH@?Alz{;dWMQiZ0=jo4|)$-S3OhU*9 zI`w(+2m^D+&I1%|AZ!p8A}ZxHF&9ndSob0kQ+78$zv?!8TD%T0p$fqaohcH?10YH| za?H_j%`PtPijU9h@v^`R&Hw* zm2(D)pF8wHun@ToEh^j|!LMHNV#r_SdtOKP^dvd=2{|I(&coUb*Q;Kr9&3NXeRhH8 zIHY5D_QM>z{KfrWy{Mi!r4a|gvcDt=BC(`nK5&P`PCy%FLdO>J2y9a0JpL& zk)(&!`XA4u3i2zFYTdCPz3X*`lj;jPz_~PWbK5&FC8<6t%YB8 zUgwyAYD#ci`Yk1#8Dt-Lzxi}|IjEM{qxn$xH5fDDlBL3?V>s92IxF&pevkYh@6MF7 z*IL2zi&gIIR!K>hpP1nj&@t#2GfUn-&`Rm9PQbF2%A?OiBM7vcu`#Vn$r8^6Fw3NK z22a|y`<1h;k)#8A)`8FG)E}Uy6J7ZbLA$?a4{Z)kOd}cbu#aQ3MB%EJ~SGW#N*P{fWjcm?zy zk1?RlMqU_au6%D_Carv~-0rH7g9)29=o%`2u^W0x+(4AwNp!O`w;PGij4b59&Z09yX3B>M(g> z{nd*#?^&Us5sS1D<+kFgXo&t;tv!z_)7;KZ^x?kqi4q`Ia9hQIR3F@;+Y}9+U|WHs zaQ#vksN>mMy@j}+D^b@jW{rB{|VGGZCmu+TQ#_uS?4G+v(ce6CdkOjYG?_ zO%F0*nbrW;tY@ zlJ#qM62o`Ggus4cy#1VLO@e1?c3PWOnEV1c%<`jF8N9hsQMT2UE(X!dWJANecB2ev zu$4G42))F4lY(;Yar7oaLm;oK%Z?P?K714;G*h-E1 z(ouQGC6ke04Ol+*8Uf=Agp>8RpEuwa!OGUBVACrjwbCSDh+}|`Pbu|EbEpOk!=H+f zp(6%xKi{0+215&)vq5rlld5MtA1BqW06D>EYw$2r4()-vvfCEdwKo@>a>4zl=0Wl$ zJ)8oaC7^+s4*Ibt)f_RQB%1xR-C^<_2|(gcoO5?~$owXV2;&5r)>3njZll^Xv!Pbt zCG!a@R6*e%576ZZA~c4cC5KfM2#Kr&hURq`&}S%O zqF;+e&&{ERc2(~T^>7&8x+Q>N|I*<;Lj_D_cL9T2TogR#@|}Q#-B$057C}}t8j2OY zcwe+NXkcXgI2a6dON1u=F}9YtKLw(>uXeUb!MiaBr*YgNAb7g>bJBCvC_XF7`eD%v zL=W*tQykoy+?xgGPzShb=?*r;g^+!5`ktK}J~{a`exrQOQ)YU+_U^($c@f| zRmD+{_|ew4JSvfif32^Z|0f>ee`zbwRi}2HTR1=G7Fk`DH0~Km2cGMLhue($tju$r z|1lD6cd8a`-THdd@5hzZQ!BUPn0UYZLH)nF)XQ=0S4|%-@jTGlxA8A(xW3-qwe;Xq zg5Vd>2nY*TJ&FOb`Y-Gn0~}mQ*YDAjG)~t)mw{&??jIlR&EuT20j{g&&`<+u3itIu zHx?#3>G|fULSN4O%G_w%LkeQO&2WeV0q+R)10Yok#L|_RN8?f8>?9=ramaYO{UL=T za@1=4s~Ccmwi9f~D!u1>XgrIy&#UF z*B=;iv~2LMVmmvFUnfD>lGeAOY&t6|Iyfq3cvwqI+lF(Ybe>A?4hM^F0;m4a)3KOF zrxC6PWCH*M)|uuqQAA^7ESv^A!DyqiA>|`6&TTasOgA)-6mSB-Yk0T(UJ**D&~Z&? z%}Ja-t);E49aRsqXM$vr3+sNDE}QKp{q_#LQ!}2u8KFE9+0!%d_}ysI_MKTdSm;SO zS*6XVE<89ILZ^KHtU8cuU@6byi15`t7h}@qjB%~_Ob=|0r=z;zf0Ii{>yyFt(67X% z*d6Z{B~m+~)GL~(kX0`wO`zu{CZ^K9mL@_o4|t)Jg}ptKhF@GROH29h-8x2cqiqDp z0A|<|#@(qHwWR&VG`7zYiRVxxntwf8@3?ZZl;+)j){ZGFycItph$-TheXfLIXD0qU zX5q#~5$xblfm|xzSq*l?Y({H@UgMGyZ9V#hMS-#s%~_>y|eP; z2R9lXPKaKyQ|4%_p|0kP>$yVQ(>t~6-G=8_w*g#&75YXm0S^7gA?LPg2K5;nY@UwY zrAOo57%M{$$}e6+gcNTGSP)Hr0DcdGdmbYWP~9oM`xcFSpz4s9V~~uAp&-@YY^xw= zytsr>a?izt%rT0H)%(;+CJI(ng?u!7VbJ0iN#v}omy*J}g4I??m}16VPj6vtV&onWQ=Wl=(}5j%i;coEyhG3Kmjme6ud=>PauGt<-wfFcsCy!hcMj!; zJ8G>=gcZPn6pG)-)BQ$4y}NWCEzfXt?VmTA14X%0XS{x|#S}b}(z)O{?giV6Awt_g zuSi7&-aZd!H*oHE{G?#s|EKljk_S)DD==nvSiTJ3HgRP{M1@U|Bv@pZ>drl1Cvv{Y z&*-U+EoCiqh|zPSZqAVviBhDwE2HNTC5H~@-ZVpTeD<&{%H#fu^Q7_?tUzUZ|9MqL zg3N*j@Ub4ZYB=}%_m`2okq?fZdD{4FO_l%r2pYFiR*}1?13f4#WSU+C8Z@Vt4!|i| za62SzVdFc-utxsAYK1*1AcsBf*P)|0Q;&X@!Si4=HCPXced%9)W+Ukf7L{c+Fxt z_3KSuScN=E4gNfn@brBN!Tt2{3HUf3xHf3}v}m?>L;bn*wt-bJ#MAij9P947d=L)0 z+HnrE$Pp0^fV>VCILX#i&G8h*lZI%`7%%$zrg=W9uyT!5^@9{lw)fv41?yQ<{SY-;AK7q=tS+yuFa0gY8dHT0m) z-%OEXZ2|E)jGgX%`D(Y^9=v5`D+N$PU){jIi)tSJn6452E+G3j2OLd2N*xB7FbnYLeF9uYDGuyyn za}!#1^r#~->j^lCDfzyVhfVWRj7I_W*{h{4F>w*q&KH=&`0qS)NZkJU!Ck}oB`8lT zD?eFKR##2Rf&t-m9&iuQS?i2*0kL`}Wc>(`8%_A4n=JVa?|{b{?LY`@*}PhxB*Vqa zTIxKUcBwVZB3YzsvwG`>IZM!;#kW6*1rb7d>)XS*>G3RcoJm?f z;Hcw_&`586#1EWnIv5_8Os_pgb(0u^rFnSPCO-rRDU-)CuC4_*DI>}Nv zc)%OT8E#x>PSZcT0`G1KRr?1oPjf~;R`*%7udbAhRHqoIQ4m&|PFLkoBboe>GDn41 zG`B7hBnZ!6Aiq#VCT&K>mlk(HjNzRRHLmVmT_1Qz>oRqt_UK@z;Qkhme4q{YTpQg} zUrW1Q=6m`8O1`DGXnS`qxC*8;->H9cNvt**yZwAToJfKC)MR#v6oeayATL7Qq)loq zXjBX+EtTOuV=UyG<2Inl!zLeTDfm`t^x)m{$XY{37?z6Gc>-cO@Ofk5M(c=KEid#A z&yTM;|Z9z`$ZW=lA;Aj=6T1m)?PPZ;+f~;`U!PR4hL~`$Eq%NYyp% zhyZ^wA#<>8@FJ#%J5GjMB{ce@8J^Z&rTrT5-ybB7netnINhLBYY5$>$T$_^hVE3*q z(EQ=?)`tV}R7_08b-#3|7|#dT2e$9j&4eivVPN5QmUwGcekhxsOhmf%6gc(wdJ~Yr z1fQ*Z(|#(hh!ev;Jl;HGY%dwJF2goyH4uDMY&1a_&%|2o;cx5V0k0l1(%mpp-kNp6 zD^CPwnZUx4!u=o65`BMfn0l@X`}j#`_akW!&r>-`s@ZBpnB4uKaNs%J!lb*gmtmEXkckI$dIRL!R8|vQS4I}ef47Z*>J)|4}lVCxF zW45dasYdG_0f?R12eLFCF1eJp^OJ4-#J+c}v6_3iCOijuAbn{*-3+rPBaa?g8fLK) zROnOO#jsmnB{SHVL8t^&z*oV8^WYK^%v zbfq}U_EM|)&x-*+u_++A5xG6T4fAM?(*+vA?B831B{c($*fsa#nNK}c>?w085y`AJ z;p0m#5Ps!^dujgdNKy>Ev0;g~Hgw_2=pd=*l`W|3T^2*62#A9Lmj?X%lM@5Fh5qPA zz0CM1wqvr=DvIHh-dR~w-~;|q9t<7H1gTy78qDH7VZ&s9m@Qrc-b|JY7s}uZ84@!W zQdtnGI!Zp9euiuB_h1Z)(u1<3b_if#>M)$ba~}Tr1%l*Mc=`}}$3H>~!v67xw|8?*}>w?W1Ls8?!f;&~f7#~%4$weY04(U7fi?`+k%z`Yx( zYfI}KV3Qyh`?0~02EtT?%#Zg3AZbI!iy!ahM;py>&Nr;+O9jT5y`SG4XFax+XWOn? z7)&pFj|h)< zsLk3F-^9eGHLV`_Qza5@9;7QP_YSVQPz1>j%-pFR*w#;{GSMplOx*a_kYubND8G!0zQ&+#uQASDGUJP+* z9`$ZqRu})OyN0cmimUh4$dOZYMEJlo+q@&jj1Lnw+2tF&0fotmNtLf-<U|LmL=$hwkpBQor`rvQ1^`=S zy@s9QkFhWkgv$cvXfB?%C{?wTpsoMPcaI52jFKC4RB@!(EV z!qbQy65T)ky8oTJQ9~8a+thv*rC+0wgcT&v2c_Ir@`W!!R!P`V5gSG6|Evdz`QXXY z%ATpROp3QX`<@*G%={^tbmM3BH;Kntnd}^HiT$r;2UG#!4*(E9{cf@@6<7W&pQ*D0 z!Gy&yQ&ftCfg!}U`7&*<_#qt{&=HbY<$lV7I&W!zb_AC!#FH6=swD?vSZ=2~{2RKV z?o=)4WodcWDT&PHoWuM7>{3v$AowV&zZ4H7MBB|tS-q>JK-X*6zEn-7B9k5xv_U-c z;Q78)hOx8lx7r)$_{AxKFYxVUv)~eVlkvVB&d-}*{|i~!kP&wLUxVRJV*kq6D$ua* zswAXzJjHlaXn$bKHUg0l4teVC*1GzvPL>|iu1kwcg^Xq9dNC0O-Vmw^qrL`1&WLs- zx^1lRg1tb_yw88+0<+KJa(<`@5GD73>6+;@0PJZD^Z5W{#GDMI$uJ#Q!D%9ZE} zt0LXZgtR_mC4TQY9qewa<~M5FYgE^K;tH~TV|Z-Q+L?1y#pMPI8oPFpr@!WJZ&zoq zj&o84R0aosQMsLwdV@sh%Oex=?nlWIiDaRuD1~<)Zwb*Nb=V6okpA8vJ#&ISYAlMm z+xiPVx843o-i6rtsts4w`diXDU*A{ObyU(2?)jsi^4q~Srfc+aiS+5VCfYM0-M-UA z^j7ECkr{G&bm40J0uMVn0p>0O2ui1~D%=`;4jWsilcnVE$UlF~i+N(2F@`-u`s z>y5MkDj!4JmxtW(83f=%I}Jn)76PlF3dk-hvT^r|><}JCT5*WmC@3aE(APT+QBg4L z!Rrc3z8L-4)0y!ei!@2>Bfj5D_<|yCTvsOcT;iB6vm*i4ED2i2*P|EeO-;1m7}mHA z$u1#UxUkzV;I`}jg$xf4iWELojylmxn?0}ymL)oqAHR-$4%0#Dr=;?)t~c`O z>wit#ov}aTN;FIi_+zrSPh#UOW+g%#AK@(VEUB_K-D~^msm;EV<0v>9Zn_2?% zXSD@;_z6`7{^-9^AAw7BodjI^sCZx1pV99__b+Ic-nhfv)kC01|EhRv#kh60d)^Y zt`JlupPIyFEe5b*H2zJ+?V`-%PXA<8K2zp~Vc#FQVC)X0&+nQ3LrKqG`JM|$4$Wgs z=)zyu%sz{~8O1pQNy63V&*Pv@t4ISP@_*1Yl~A>Y+eFG1d4+@|x|@^(+$@y1TlY!J zebZm-jV}4Vis5NXzp=t(I#(JRWc&Z&Rr?b1hcf(wn4QAhR_56EnUv*7Un(j337?b+ zjUxxZEujt4#q*Fx2Dz24TnyO)cY;uiMA;lOU!G+wb*3c5{H6P1ZIf>Q6Tb0O;8#*f z3AP`l4{q*q!6!lJ^;w>L&t_t0_mXiwdR4vpzt;lWALj~y`1CV|Dw%(S7c=!yP!Kk) z{FJzXdnMVvWo&Maa)RS-pCmnaj}QrV5a-h*9wr-H7t>+==hG;{rd4rpKI1P58MxoTmXjF*k-m;ZF@hGki1VBO2(N0^g#)ypb{##&4)hRka5e{u?}lzQr`f(sGq*!k6p zXkM+fzu&^}e&=YwFG28G#>I@@9Vx&#$hmUy)8{SC1%kqUvnNAwFD-3j(H3>Qyyw*s ze5@jT*5KPdtF4wojEN5ZHVSts&Xg$^{Q3PO*N0de``F+I;Kl=l5Vg&8B#@KaJ>9FPk>gq zFnBeZt zhW0$S!FG_-(Ayr(fW70PjVC+-duNCkr|*si>L`aMs}KGY++I^C9&dGFNe5NtxuxFt zpFdW3?;8+XkCH~5d@5$e9i703{`W10cfLXD2lB}d7Z0LV+oy(x@bKOUO4GW+C;*(8 zYjwAIRJ@5^JUA$F|D237jdYw99!>;N4*=sjt1_Xfv9T28KrZz9-urW|cVjB;*-oU{ zq8RtZU-ybe@9GqN3oGO5(V3v0N9^{9d&RXsJOWeXvd`*#(cU%p#Uak@Of*PVF0IP^ z^8ung0&wOVz)TFdm-(4$5u_JXcyAz%)3r?TpKmk@QrsX>q28wcZs=!7%lrE^)1N}k zDnJ@H?!T5eak105fB*j`+T*f+|NgK49hrX2!x$<0bp6?pD<1q~cowNwrehcSe*mPC B|D*r_ diff --git a/tests/integration/conftest.py b/tests/integration/conftest.py index bc621138..2a74c89f 100644 --- a/tests/integration/conftest.py +++ b/tests/integration/conftest.py @@ -101,13 +101,6 @@ def vasp_shape_data(vasp_traj): return shapes -@pytest.fixture(scope='module') -def vasp_full_vol(vasp_full_traj): - trajectory = vasp_full_traj - diff_trajectory = trajectory.filter('Li') - return trajectory_to_volume(trajectory=diff_trajectory, resolution=0.3) - - @pytest.fixture(scope='module') def vasp_path_vol(vasp_full_traj): trajectory = vasp_full_traj @@ -116,7 +109,7 @@ def vasp_path_vol(vasp_full_traj): @pytest.fixture(scope='module') -def vasp_full_path(vasp_path_vol): +def vasp_path(vasp_path_vol): peaks = vasp_path_vol.find_peaks() F = vasp_path_vol.get_free_energy(temperature=650.0) path = find_best_perc_path(F, diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index 32105dd9..349aed5a 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -10,8 +10,8 @@ @pytest.vaspxml_available # type: ignore -def test_fractional_coordinates(vasp_path_vol, vasp_full_path): - frac_sites = vasp_path_vol.voxel_to_frac_coords(vasp_full_path.sites) +def test_fractional_coordinates(vasp_path_vol, vasp_path): + frac_sites = vasp_path_vol.voxel_to_frac_coords(vasp_path.sites) assert isclose(frac_sites[-1][0], 0.4107142857142857) assert isclose(frac_sites[19][1], 0.6785714285714286) @@ -39,25 +39,25 @@ def test_optimal_path(vasp_F_graph, start, stop, method, expected): @pytest.vaspxml_available # type: ignore -def test_find_best_perc_path(vasp_full_path): - assert isclose(vasp_full_path.total_energy, 11.488013690080908) - assert vasp_full_path.start_site == (11, 9, 6) +def test_find_best_perc_path(vasp_path): + assert isclose(vasp_path.total_energy, 11.488013690080908) + assert vasp_path.start_site == (11, 9, 6) @pytest.vaspxml_available # type: ignore -def test_nearest_structure_reference(vasp_full_vol, vasp_full_path): +def test_nearest_structure_reference(vasp_path_vol, vasp_path): structure = load_known_material('argyrodite') - nearest_structure_label, nearest_structure_coord = vasp_full_path.path_over_structure( - structure, vasp_full_vol) + nearest_structure_label, nearest_structure_coord = vasp_path.path_over_structure( + structure, vasp_path_vol) assert nearest_structure_label[0] == '48h' assert nearest_structure_label[10] == '48h' - assert isclose(nearest_structure_coord[0][0], 0.2381760000000003) - assert isclose(nearest_structure_coord[0][1], 1.8160919999999998) - assert isclose(nearest_structure_coord[20][2], 3.145908) - assert isclose(nearest_structure_coord[-1][-1], 1.816092) + assert isclose(nearest_structure_coord[0][0], 3.145908) + assert isclose(nearest_structure_coord[0][1], 8.107908) + assert isclose(nearest_structure_coord[20][2], 5.200176) + assert isclose(nearest_structure_coord[-1][-1], 5.200176) @pytest.vaspxml_available # type: ignore diff --git a/tests/integration/plot_test.py b/tests/integration/plot_test.py index 272cad7d..27c8dfb5 100644 --- a/tests/integration/plot_test.py +++ b/tests/integration/plot_test.py @@ -81,10 +81,10 @@ def test_msd_per_element(vasp_traj): @image_comparison2(baseline_images=['path_energy']) -def test_path_energy(vasp_full_vol, vasp_full_path): +def test_path_energy(vasp_path_vol, vasp_path): structure = load_known_material('argyrodite') - plots.energy_along_path(paths=vasp_full_path, - volume=vasp_full_vol, + plots.energy_along_path(paths=vasp_path, + volume=vasp_path_vol, structure=structure) From 34c514ccf422ae0cee56e9174805adde07530a86 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 14:58:34 +0200 Subject: [PATCH 03/24] Refactor Pathway --- src/gemdat/path.py | 56 ++++++++++++++++++++------- src/gemdat/plots/matplotlib/_paths.py | 2 +- tests/integration/path_test.py | 6 ++- 3 files changed, 46 insertions(+), 18 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index 0b099ece..1c81b25c 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -25,16 +25,21 @@ class Pathway: List of voxel coordinates of the sites defining the path energy: list[float] List of the energy along the path + dims: [int, int, int] | None + Voxel dimensions of bounding box. If set (usually to `Volume.dims`), + enable some site transformations. """ sites: list[tuple[int, int, int]] energy: list[float] + dims: tuple[int, int, int] | None = None def __repr__(self): s = ( f'Path: {self.start_site} -> {self.stop_site}', f'Steps: {len(self.sites)}', f'Total energy: {self.total_energy:.3f} eV', + f'Dimensions: {self.dims}', ) return '\n'.join(s) @@ -44,21 +49,43 @@ def total_energy(self): """Return total energy for path.""" return sum(self.energy) - def wrap(self, dims: tuple[int, int, int]): - """Wrap path in periodic boundary conditions in-place. + def wrapped_sites(self) -> list[tuple[int, int, int]]: + """Wrap sites to bounding box. + + Returns + ------- + np.ndarray + Voxel coordinates wrapped to bounding box. + """ + if not self.dims: + raise AttributeError( + f'Dimensions are needed for this method {self.dims=}') + xdim, xdim, xdim = self.dims + return [(x % xdim, y % xdim, z % xdim) for x, y, z in self.sites] + + def frac_sites(self, wrapped: bool = False) -> np.ndarray: + """Return fractional sites. Parameters ---------- - F: np.ndarray - Grid in which the path sites will be wrapped + wrapped : bool + If True, wrap coordinates to bounding box using + `.wrapped_sites()`. + + Returns + ------- + np.ndarray + Fractional coordinates for sites """ - X, Y, Z = dims - self.sites = [(x % X, y % Y, z % Z) for x, y, z in self.sites] + if not self.dims: + raise AttributeError( + f'Dimensions are needed for this method {self.dims=}') + sites = self.wrapped_sites() if wrapped else self.sites + return (np.array(sites) + 0.5) / np.array(self.dims) def path_over_structure( self, structure: Structure, - volume: Volume, ) -> tuple[list[str], list[np.ndarray]]: """Find the nearest site of the structure to the path sites. @@ -66,8 +93,6 @@ def path_over_structure( ---------- structure : Structure Reference structure - volume : Volume - Volume object that contains the information about the nearest sites of the structure Returns ------- @@ -76,7 +101,8 @@ def path_over_structure( nearest_structure_coord: list[np.ndarray] List of cartesian coordinates of the closest site of the reference structure """ - frac_sites = volume.voxel_to_frac_coords(np.array(self.sites)) + frac_sites = np.array(self.frac_sites(wrapped=True)) + nearest_structure_tree, nearest_structure_map = nearest_structure_reference( structure) @@ -94,6 +120,8 @@ def path_over_structure( for index in nearest_structure_indices ] + # TODO: Return as dict? + return nearest_structure_label, nearest_structure_coord @property @@ -420,8 +448,6 @@ def find_best_perc_path(F: Volume, best_percolating_path: Pathway Optimal path that percolates the graph in the specified directions """ - xyz_real = F.dims - # Find percolation using virtual images along the required dimensions if not any([percolate_x, percolate_y, percolate_z]): raise ValueError('percolation is not defined') @@ -438,7 +464,7 @@ def find_best_perc_path(F: Volume, # reaching the percolating image image = tuple( x * px - for x, px in zip(xyz_real, (percolate_x, percolate_y, percolate_z))) + for x, px in zip(F.dims, (percolate_x, percolate_y, percolate_z))) # Find the lowest cost path that percolates along the x dimension best_cost = float('inf') @@ -466,7 +492,7 @@ def find_best_perc_path(F: Volume, best_path = path if best_path: - # Before returning, wrap the path in the original volume - best_path.wrap(F.dims) + # Before returning, set dimensions of original volume + best_path.dims = F.dims return best_path diff --git a/src/gemdat/plots/matplotlib/_paths.py b/src/gemdat/plots/matplotlib/_paths.py index 903e1b11..acf00b1c 100644 --- a/src/gemdat/plots/matplotlib/_paths.py +++ b/src/gemdat/plots/matplotlib/_paths.py @@ -46,7 +46,7 @@ def energy_along_path( ax.set(ylabel='Free energy [eV]') nearest_structure_label, nearest_structure_coord = path.path_over_structure( - structure, volume) + structure) # Create costum labels for the x axis to avoid consecutive repetitions site_xlabel = [] diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index 349aed5a..22e69c35 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -11,7 +11,7 @@ @pytest.vaspxml_available # type: ignore def test_fractional_coordinates(vasp_path_vol, vasp_path): - frac_sites = vasp_path_vol.voxel_to_frac_coords(vasp_path.sites) + frac_sites = vasp_path.frac_sites(wrapped=True) assert isclose(frac_sites[-1][0], 0.4107142857142857) assert isclose(frac_sites[19][1], 0.6785714285714286) @@ -48,8 +48,10 @@ def test_find_best_perc_path(vasp_path): def test_nearest_structure_reference(vasp_path_vol, vasp_path): structure = load_known_material('argyrodite') + assert vasp_path_vol.dims == vasp_path.dims + nearest_structure_label, nearest_structure_coord = vasp_path.path_over_structure( - structure, vasp_path_vol) + structure) assert nearest_structure_label[0] == '48h' assert nearest_structure_label[10] == '48h' From 061365acc772f02dd6156a89afb5aebb5d3bcae4 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 15:23:33 +0200 Subject: [PATCH 04/24] Add method to calculate path length --- src/gemdat/path.py | 26 +++++++++++++++++++++++++- tests/integration/path_test.py | 8 ++++++++ 2 files changed, 33 insertions(+), 1 deletion(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index 1c81b25c..b5a744fe 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -4,16 +4,21 @@ from __future__ import annotations from dataclasses import dataclass -from typing import Literal +from itertools import pairwise +from typing import TYPE_CHECKING, Literal import networkx as nx import numpy as np from pymatgen.core import Structure +from pymatgen.core.units import FloatWithUnit from gemdat.volume import Volume from .utils import nearest_structure_reference +if TYPE_CHECKING: + from pymatgen.core import Lattice + @dataclass class Pathway: @@ -49,6 +54,25 @@ def total_energy(self): """Return total energy for path.""" return sum(self.energy) + def total_length(self, lattice: Lattice) -> FloatWithUnit: + """Return total length of pathway in Ã…ngstrom. + + Parameters + ---------- + lattice : Lattice + Lattice parameters + + Returns + ------- + length : FloatWithUnit + Total distance in Ã…ngstrom + """ + length = 0 + for a, b in pairwise(self.frac_sites()): + dist, _ = lattice.get_distance_and_image(a, b) + length += dist + return FloatWithUnit(length, 'ang') + def wrapped_sites(self) -> list[tuple[int, int, int]]: """Wrap sites to bounding box. diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index 22e69c35..3db92584 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -44,6 +44,14 @@ def test_find_best_perc_path(vasp_path): assert vasp_path.start_site == (11, 9, 6) +@pytest.vaspxml_available # type: ignore +def test_path_length(vasp_path_vol, vasp_path): + structure = load_known_material('argyrodite') + length = vasp_path.total_length(structure.lattice) + assert length == 20.10367043479136 + assert str(length.unit) == 'ang' + + @pytest.vaspxml_available # type: ignore def test_nearest_structure_reference(vasp_path_vol, vasp_path): structure = load_known_material('argyrodite') From bb028c32fab6ed6b1d99fa4a33529f55e7bd3b98 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 15:43:41 +0200 Subject: [PATCH 05/24] Return PeriodicSite from path_over_structure --- src/gemdat/path.py | 23 +++++++---------------- tests/integration/path_test.py | 16 ++++++++-------- 2 files changed, 15 insertions(+), 24 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index b5a744fe..d1521799 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -17,7 +17,7 @@ from .utils import nearest_structure_reference if TYPE_CHECKING: - from pymatgen.core import Lattice + from pymatgen.core import Lattice, PeriodicSite @dataclass @@ -110,7 +110,7 @@ def frac_sites(self, wrapped: bool = False) -> np.ndarray: def path_over_structure( self, structure: Structure, - ) -> tuple[list[str], list[np.ndarray]]: + ) -> list[PeriodicSite]: """Find the nearest site of the structure to the path sites. Parameters @@ -120,10 +120,8 @@ def path_over_structure( Returns ------- - nearest_structure_label: list[str] - List of the label of the closest site of the reference structure - nearest_structure_coord: list[np.ndarray] - List of cartesian coordinates of the closest site of the reference structure + nearest_sites: list[PeriodicSite] + List closest sites of the reference structure """ frac_sites = np.array(self.frac_sites(wrapped=True)) @@ -135,18 +133,11 @@ def path_over_structure( nearest_structure_tree.query(site)[1] for site in frac_sites ] # and use it to get its label and coordinates - nearest_structure_label = [ - structure.labels[nearest_structure_map[index]] + nearest_sites = [ + structure[nearest_structure_map[index]] for index in nearest_structure_indices ] - nearest_structure_coord = [ - structure.cart_coords[nearest_structure_map[index]] - for index in nearest_structure_indices - ] - - # TODO: Return as dict? - - return nearest_structure_label, nearest_structure_coord + return nearest_sites @property def start_site(self) -> tuple[int, int, int]: diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index 3db92584..b9238ec5 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -58,16 +58,16 @@ def test_nearest_structure_reference(vasp_path_vol, vasp_path): assert vasp_path_vol.dims == vasp_path.dims - nearest_structure_label, nearest_structure_coord = vasp_path.path_over_structure( - structure) + nearest_sites = vasp_path.path_over_structure(structure) - assert nearest_structure_label[0] == '48h' - assert nearest_structure_label[10] == '48h' + assert all(s.label == '48h' for s in nearest_sites) - assert isclose(nearest_structure_coord[0][0], 3.145908) - assert isclose(nearest_structure_coord[0][1], 8.107908) - assert isclose(nearest_structure_coord[20][2], 5.200176) - assert isclose(nearest_structure_coord[-1][-1], 5.200176) + np.testing.assert_allclose(nearest_sites[0].coords, + [3.145908, 8.107908, 5.200176]) + np.testing.assert_allclose(nearest_sites[20].coords, + [1.816092, 6.778092, 5.200176]) + np.testing.assert_allclose(nearest_sites[-1].coords, + [3.145908, 8.107908, 5.200176]) @pytest.vaspxml_available # type: ignore From d1faa1e4d34e7e60b0e39df1d436a5556109ce10 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 15:45:30 +0200 Subject: [PATCH 06/24] Use np.testing.assert_allclose --- tests/integration/path_test.py | 14 ++++++-------- tests/rotations_test.py | 10 +++++----- 2 files changed, 11 insertions(+), 13 deletions(-) diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index b9238ec5..39be8df2 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -4,6 +4,7 @@ import numpy as np import pytest +from numpy.testing import assert_allclose from gemdat.io import load_known_material from gemdat.path import multiple_paths, optimal_path @@ -58,16 +59,13 @@ def test_nearest_structure_reference(vasp_path_vol, vasp_path): assert vasp_path_vol.dims == vasp_path.dims - nearest_sites = vasp_path.path_over_structure(structure) + nearest = vasp_path.path_over_structure(structure) - assert all(s.label == '48h' for s in nearest_sites) + assert all(s.label == '48h' for s in nearest) - np.testing.assert_allclose(nearest_sites[0].coords, - [3.145908, 8.107908, 5.200176]) - np.testing.assert_allclose(nearest_sites[20].coords, - [1.816092, 6.778092, 5.200176]) - np.testing.assert_allclose(nearest_sites[-1].coords, - [3.145908, 8.107908, 5.200176]) + assert_allclose(nearest[0].coords, [3.145908, 8.107908, 5.200176]) + assert_allclose(nearest[20].coords, [1.816092, 6.778092, 5.200176]) + assert_allclose(nearest[-1].coords, [3.145908, 8.107908, 5.200176]) @pytest.vaspxml_available # type: ignore diff --git a/tests/rotations_test.py b/tests/rotations_test.py index 6de6de6c..4a72c1c5 100644 --- a/tests/rotations_test.py +++ b/tests/rotations_test.py @@ -3,6 +3,7 @@ from math import isclose import numpy as np +from numpy.testing import assert_allclose from pymatgen.core import Species from gemdat.rotations import ( @@ -32,8 +33,8 @@ def test_normalize(trajectory): in_vectors=np.array([[1, 2, 2], [2, 2, 1]], dtype=float)) ret = orientations.normalize() - assert np.allclose( - ret.vectors, np.array([[1 / 3, 2 / 3, 2 / 3], [2 / 3, 2 / 3, 1 / 3]])) + assert_allclose(ret.vectors, + np.array([[1 / 3, 2 / 3, 2 / 3], [2 / 3, 2 / 3, 1 / 3]])) def test_conventional(trajectory): @@ -46,8 +47,7 @@ def test_conventional(trajectory): dtype=float)) matrix = np.eye(3) * [1, 2, 3] ret = orientations.transform(matrix=matrix) - assert np.allclose(ret.vectors, np.array([[1, 0, 0], [0, 2, 0], [0, 0, - 3]])) + assert_allclose(ret.vectors, np.array([[1, 0, 0], [0, 2, 0], [0, 0, 3]])) def test_symmetrize(trajectory): @@ -59,7 +59,7 @@ def test_symmetrize(trajectory): dtype=float)) sym_ops = np.array([[0, -1, 0], [1, 0, 0], [0, 0, -1]]) ret = orientations.symmetrize(sym_ops=sym_ops) - assert np.allclose( + assert_allclose( ret.vectors, np.array([[[0., 1., 0.], [-1., 0., 0.], [0., 0., -1.]], [[0., 1., 0.], [-1., 0., 0.], [0., 0., -1.]]])) From 6ada1a2702a9455cf654c79c37833a0acdb33a8e Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 16:22:57 +0200 Subject: [PATCH 07/24] Fix bug --- src/gemdat/path.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index d1521799..f24e90d9 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -84,7 +84,7 @@ def wrapped_sites(self) -> list[tuple[int, int, int]]: if not self.dims: raise AttributeError( f'Dimensions are needed for this method {self.dims=}') - xdim, xdim, xdim = self.dims + xdim, ydim, zdim = self.dims return [(x % xdim, y % xdim, z % xdim) for x, y, z in self.sites] def frac_sites(self, wrapped: bool = False) -> np.ndarray: From 8c8316b6fcdc269259e7b173a69355eba6581fd8 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 16:25:54 +0200 Subject: [PATCH 08/24] Update plot to work with PeriodicSites --- src/gemdat/plots/matplotlib/_paths.py | 31 +++++++++++++-------------- tests/integration/plot_test.py | 4 +--- 2 files changed, 16 insertions(+), 19 deletions(-) diff --git a/src/gemdat/plots/matplotlib/_paths.py b/src/gemdat/plots/matplotlib/_paths.py index acf00b1c..baae932e 100644 --- a/src/gemdat/plots/matplotlib/_paths.py +++ b/src/gemdat/plots/matplotlib/_paths.py @@ -4,14 +4,12 @@ from pymatgen.core import Structure from gemdat.path import Pathway -from gemdat.volume import Volume def energy_along_path( *, paths: Pathway | list[Pathway], structure: Structure, - volume: Volume, ) -> plt.Figure: """Plot energy along specified path. @@ -21,8 +19,6 @@ def energy_along_path( Pathway object containing the energy along the path, or list of Pathways structure : Structure Structure object to get the site information - volume : Volume - Volume object to get the site information Returns ------- @@ -35,41 +31,43 @@ def energy_along_path( else: path = paths - if path.energy is None: - raise ValueError('Pathway does not contain energy data') - if path.sites is None: - raise ValueError('Pathway does not contain site data') - fig, ax = plt.subplots(figsize=(8, 4)) ax.plot(path.energy, marker='o', color='r', label='Optimal path') ax.set(ylabel='Free energy [eV]') - nearest_structure_label, nearest_structure_coord = path.path_over_structure( - structure) + nearest_sites = path.path_over_structure(structure) # Create costum labels for the x axis to avoid consecutive repetitions site_xlabel = [] sitecoord_xlabel = [] - for i in range(len(path.sites)): - coord = nearest_structure_coord[i] + + prev = nearest_sites[0] + for i, site in enumerate(nearest_sites): # only non repeated labels will get an entry - if (coord != nearest_structure_coord[i - 1]).any() or i == 0: - sitecoord_xlabel.append(', '.join([f'{val:.1f}' for val in coord])) - site_xlabel.append(f'{nearest_structure_label[i]}') + if (site.coords != prev.coords).any() or i == 0: + sitecoord_xlabel.append(', '.join(f'{val:.1f}' + for val in site.coords)) + site_xlabel.append(site.label) else: sitecoord_xlabel.append('') site_xlabel.append('') + + prev = site + non_empty_ticks = [ i for i, label in enumerate(sitecoord_xlabel) if label != '' ] + extra_ticks = non_empty_ticks.copy() extra_ticks.append(ax.get_xlim()[1]) centered_ticks = [(extra_ticks[i] + extra_ticks[i + 1]) / 2 for i in range(len(extra_ticks) - 1)] + ax.set_xticks(centered_ticks) ax.set_xticklabels([sitecoord_xlabel[i] for i in non_empty_ticks], rotation=45) + # Change background color alternatively for different sites for i in range(0, len(non_empty_ticks), 2): if i + 1 < len(non_empty_ticks): @@ -82,6 +80,7 @@ def energy_along_path( max(non_empty_ticks), facecolor='lightgray', edgecolor='none') + # and add on top the site labels ax_up = ax.twiny() ax_up.set_xlim(ax.get_xlim()) diff --git a/tests/integration/plot_test.py b/tests/integration/plot_test.py index 27c8dfb5..af78e535 100644 --- a/tests/integration/plot_test.py +++ b/tests/integration/plot_test.py @@ -83,9 +83,7 @@ def test_msd_per_element(vasp_traj): @image_comparison2(baseline_images=['path_energy']) def test_path_energy(vasp_path_vol, vasp_path): structure = load_known_material('argyrodite') - plots.energy_along_path(paths=vasp_path, - volume=vasp_path_vol, - structure=structure) + plots.energy_along_path(paths=vasp_path, structure=structure) @image_comparison2(baseline_images=['rectilinear']) From d0073f05e1a07c9f555ef0ce581ea362af6289b4 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 15 Apr 2024 16:34:13 +0200 Subject: [PATCH 09/24] Add shortcuts to path plots as methods --- src/gemdat/path.py | 10 ++++++++++ src/gemdat/plots/matplotlib/_paths.py | 23 ++++++++--------------- tests/integration/plot_test.py | 4 ++-- 3 files changed, 20 insertions(+), 17 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index f24e90d9..00a5b307 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -149,6 +149,16 @@ def stop_site(self) -> tuple[int, int, int]: """Return stop site.""" return self.sites[-1] + def plot_energy_along_path(self, **kwargs): + """See [gemdat.plots.energy_along_path][] for more info.""" + from gemdat import plots + return plots.energy_along_path(path=self, **kwargs) + + def plot_path_on_grid(self, **kwargs): + """See [gemdat.plots.path_on_grid][] for more info.""" + from gemdat import plots + return plots.path_on_grid(path=self, **kwargs) + def free_energy_graph(F: np.ndarray | Volume, max_energy_threshold: float = 1e20, diff --git a/src/gemdat/plots/matplotlib/_paths.py b/src/gemdat/plots/matplotlib/_paths.py index baae932e..911e69e7 100644 --- a/src/gemdat/plots/matplotlib/_paths.py +++ b/src/gemdat/plots/matplotlib/_paths.py @@ -7,29 +7,27 @@ def energy_along_path( + path: Pathway, *, - paths: Pathway | list[Pathway], structure: Structure, + other_paths: list[Pathway] | None = None, ) -> plt.Figure: """Plot energy along specified path. Parameters ---------- - paths : Pathway | list[Pathway] + path : Pathway | list[Pathway] Pathway object containing the energy along the path, or list of Pathways structure : Structure Structure object to get the site information + other_paths : Pathway | list[Pathway] + Optional list of alternative paths to plot Returns ------- fig : matplotlib.figure.Figure Output figure """ - # The first Pathway in paths is assumed to be the optimal one - if isinstance(paths, list): - path = paths[0] - else: - path = paths fig, ax = plt.subplots(figsize=(8, 4)) @@ -90,8 +88,8 @@ def energy_along_path( ax_up.get_yaxis().set_visible(False) # If available, plot the other pathways - if isinstance(paths, list): - for idx, path in enumerate(paths[1:]): + if other_paths: + for idx, path in enumerate(other_paths): if path.energy is None: raise ValueError('Pathway does not contain energy data') ax.plot(range(len(path.energy)), @@ -102,7 +100,7 @@ def energy_along_path( return fig -def path_on_grid(*, path: Pathway) -> plt.Figure: +def path_on_grid(path: Pathway) -> plt.Figure: """Plot the 3d coordinates of the points that define a path. Parameters @@ -115,11 +113,6 @@ def path_on_grid(*, path: Pathway) -> plt.Figure: fig : matplotlib.figure.Figure Output figure """ - if path.energy is None: - raise ValueError('Pathway does not contain energy data') - if path.sites is None: - raise ValueError('Pathway does not contain site data') - # Create a colormap to visualize the path colormap = plt.get_cmap() normalize = plt.Normalize(0, len(path.energy)) diff --git a/tests/integration/plot_test.py b/tests/integration/plot_test.py index af78e535..07ea09dd 100644 --- a/tests/integration/plot_test.py +++ b/tests/integration/plot_test.py @@ -81,9 +81,9 @@ def test_msd_per_element(vasp_traj): @image_comparison2(baseline_images=['path_energy']) -def test_path_energy(vasp_path_vol, vasp_path): +def test_path_energy(vasp_path): structure = load_known_material('argyrodite') - plots.energy_along_path(paths=vasp_path, structure=structure) + plots.energy_along_path(path=vasp_path, structure=structure) @image_comparison2(baseline_images=['rectilinear']) From 291fae621044410647911d3ee5fc301840d1385b Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 11:29:51 +0200 Subject: [PATCH 10/24] Optimize api for generating Pathways from Volumes --- src/gemdat/path.py | 7 +++-- src/gemdat/volume.py | 63 ++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 66 insertions(+), 4 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index 00a5b307..1993d885 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -353,9 +353,9 @@ def optimal_path( ---------- F_graph : nx.Graph Graph of the free energy - start : tuple + start : Collection Coordinates of the starting point - stop: tuple + stop: Collection Coordinates of the stoping point method : str Method used to calculate the shortest path. Options are: @@ -380,6 +380,9 @@ def optimal_path( if method in ('dijkstra-exp', 'minmax-energy', 'simple'): method = 'dijkstra' + start = tuple(start) + stop = tuple(stop) + optimal_path = nx.shortest_path(F_graph, source=start, target=stop, diff --git a/src/gemdat/volume.py b/src/gemdat/volume.py index 5b2cb71e..431a5730 100644 --- a/src/gemdat/volume.py +++ b/src/gemdat/volume.py @@ -2,7 +2,7 @@ from __future__ import annotations -from dataclasses import dataclass +from dataclasses import dataclass, field from pathlib import Path from typing import TYPE_CHECKING, Any, Optional @@ -18,6 +18,7 @@ from .segmentation import watershed_pbc if TYPE_CHECKING: + import networkx as nx from pymatgen.core import Lattice, PeriodicSite from skimage.measure._regionprops import RegionProperties @@ -42,11 +43,69 @@ class Volume: data: np.ndarray lattice: Lattice label: str = 'volume' - units: Unit | None = None + units: Unit = field(default_factory=lambda: Unit('')) def __post_init__(self): self.dims = self.data.shape + def __repr__(self): + units = f' ({self.units})' if self.units else '' + + def to_str(x): + return f'{x:>10.6f}' + + abc = ' '.join(to_str(i) for i in self.lattice.abc) + angles = ' '.join(to_str(i) for i in self.lattice.angles) + pbc = ' '.join(str(p).rjust(10) for p in self.lattice.pbc) + + s = ( + f'Data: {self.data.shape}', + f'Lattice, abc: {abc}', + f' angles: {angles}', + f' pbc: {pbc}', + f'Label: {self.label}{units}', + f'Dimensions: {self.dims}', + ) + return '\n'.join(s) + + def free_energy_graph(self, **kwargs) -> nx.Graph: + """Compute the graph of the free energy for networkx functions. + + See [gemdat.path.free_energy_graph][] for more info. + """ + from .path import free_energy_graph + return free_energy_graph(self.data, **kwargs) + + def optimal_path(self, *args, F_graph: nx.Graph | None = None, **kwargs): + """Calculate the shortest cost-effective path using the desired method. + + Parameters + ---------- + F_graph : Graph | None + Optionally, define your own free energy graph. Otherwise, + it will be calculated on the fly using default parameters. + **kwargs: + These parameters are passed to [gemdat.path.optimal_path][]. + See [gemdat.path.optimal_path][] for more info. + + Returns + ------- + path : Pathway + Voxel coordinates and energy of optimal path from start to stop. + """ + if not self.label == 'free_energy': + return TypeError( + f'Cannot calculate path for class of type `{self.label}`') + + from .path import optimal_path + + if not F_graph: + F_graph = self.free_energy_graph(max_energy_threshold=1e7) + + path = optimal_path(F_graph, *args, **kwargs) + path.dims = self.dims + return path + def normalized(self) -> np.ndarray: """Return normalized data.""" return self.data / self.data.max() From 1d41a252f0864ce52c5b9adfc28ed61e3edba4b1 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 11:40:45 +0200 Subject: [PATCH 11/24] Subclass to FreeEnergyVolume --- src/gemdat/path.py | 12 +++---- src/gemdat/volume.py | 79 +++++++++++++++++++++----------------------- 2 files changed, 44 insertions(+), 47 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index 1993d885..756abf85 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -12,7 +12,7 @@ from pymatgen.core import Structure from pymatgen.core.units import FloatWithUnit -from gemdat.volume import Volume +from gemdat.volume import FreeEnergyVolume from .utils import nearest_structure_reference @@ -160,14 +160,14 @@ def plot_path_on_grid(self, **kwargs): return plots.path_on_grid(path=self, **kwargs) -def free_energy_graph(F: np.ndarray | Volume, +def free_energy_graph(F: np.ndarray | FreeEnergyVolume, max_energy_threshold: float = 1e20, diagonal: bool = True) -> nx.Graph: """Compute the graph of the free energy for networkx functions. Parameters ---------- - F : np.ndarray | Volume + F : np.ndarray | FreeEnergyVolume Free energy on the 3d grid max_energy_threshold : float, optional Maximum energy threshold for the path to be considered valid @@ -194,7 +194,7 @@ def free_energy_graph(F: np.ndarray | Volume, G = nx.Graph() - data = F.data if isinstance(F, Volume) else F + data = F.data if isinstance(F, FreeEnergyVolume) else F for index, Fi in np.ndenumerate(data): if 0 <= Fi < max_energy_threshold: @@ -451,7 +451,7 @@ def _optimal_path_minmax_energy( return optimal_path -def find_best_perc_path(F: Volume, +def find_best_perc_path(F: FreeEnergyVolume, peaks: np.ndarray, percolate_x: bool = True, percolate_y: bool = False, @@ -460,7 +460,7 @@ def find_best_perc_path(F: Volume, Parameters ---------- - F : Volume + F : FreeEnergyVolume Energy grid that will be used to calculate the shortest path peaks : np.ndarray List of the peaks that correspond to high probability regions diff --git a/src/gemdat/volume.py b/src/gemdat/volume.py index 431a5730..08c81cfa 100644 --- a/src/gemdat/volume.py +++ b/src/gemdat/volume.py @@ -68,44 +68,6 @@ def to_str(x): ) return '\n'.join(s) - def free_energy_graph(self, **kwargs) -> nx.Graph: - """Compute the graph of the free energy for networkx functions. - - See [gemdat.path.free_energy_graph][] for more info. - """ - from .path import free_energy_graph - return free_energy_graph(self.data, **kwargs) - - def optimal_path(self, *args, F_graph: nx.Graph | None = None, **kwargs): - """Calculate the shortest cost-effective path using the desired method. - - Parameters - ---------- - F_graph : Graph | None - Optionally, define your own free energy graph. Otherwise, - it will be calculated on the fly using default parameters. - **kwargs: - These parameters are passed to [gemdat.path.optimal_path][]. - See [gemdat.path.optimal_path][] for more info. - - Returns - ------- - path : Pathway - Voxel coordinates and energy of optimal path from start to stop. - """ - if not self.label == 'free_energy': - return TypeError( - f'Cannot calculate path for class of type `{self.label}`') - - from .path import optimal_path - - if not F_graph: - F_graph = self.free_energy_graph(max_energy_threshold=1e7) - - path = optimal_path(F_graph, *args, **kwargs) - path.dims = self.dims - return path - def normalized(self) -> np.ndarray: """Return normalized data.""" return self.data / self.data.max() @@ -404,11 +366,9 @@ def get_free_energy( free_energy = -temperature * physical_constants[ 'Boltzmann constant in eV/K'][0] * np.log(prob) - return Volume( + return FreeEnergyVolume( data=np.nan_to_num(free_energy), lattice=self.lattice, - label='free_energy', - units=Unit('eV K-1'), ) def plot_3d(self, **kwargs): @@ -417,6 +377,43 @@ def plot_3d(self, **kwargs): return plots.plot_3d(volume=self, **kwargs) +class FreeEnergyVolume(Volume): + + def free_energy_graph(self, **kwargs) -> nx.Graph: + """Compute the graph of the free energy for networkx functions. + + See [gemdat.path.free_energy_graph][] for more info. + """ + from .path import free_energy_graph + return free_energy_graph(self.data, **kwargs) + + def optimal_path(self, *args, F_graph: nx.Graph | None = None, **kwargs): + """Calculate the shortest cost-effective path using the desired method. + + Parameters + ---------- + F_graph : Graph | None + Optionally, define your own free energy graph. Otherwise, + it will be calculated on the fly using default parameters. + **kwargs: + These parameters are passed to [gemdat.path.optimal_path][]. + See [gemdat.path.optimal_path][] for more info. + + Returns + ------- + path : Pathway + Voxel coordinates and energy of optimal path from start to stop. + """ + from .path import optimal_path + + if not F_graph: + F_graph = self.free_energy_graph(max_energy_threshold=1e7) + + path = optimal_path(F_graph, *args, **kwargs) + path.dims = self.dims + return path + + def trajectory_to_volume( trajectory: Trajectory, resolution: float = 0.2, From 68fc29078b2873f16e0ef809026b976bdf93a3e8 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 13:09:25 +0200 Subject: [PATCH 12/24] Remove redundant date.shape --- src/gemdat/volume.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/gemdat/volume.py b/src/gemdat/volume.py index 08c81cfa..dbfbc385 100644 --- a/src/gemdat/volume.py +++ b/src/gemdat/volume.py @@ -59,12 +59,11 @@ def to_str(x): pbc = ' '.join(str(p).rjust(10) for p in self.lattice.pbc) s = ( - f'Data: {self.data.shape}', + f'Data: {self.dims}', f'Lattice, abc: {abc}', f' angles: {angles}', f' pbc: {pbc}', f'Label: {self.label}{units}', - f'Dimensions: {self.dims}', ) return '\n'.join(s) From e65be6905c821dc75348dd076240bf2051297d4d Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 13:43:41 +0200 Subject: [PATCH 13/24] Add shortcut for multiple/percolating_paths on volume --- src/gemdat/path.py | 36 ++++++++++++++++++------------ src/gemdat/plots/plotly/_plot3d.py | 3 ++- src/gemdat/volume.py | 33 +++++++++++++++++++++++++-- tests/integration/path_test.py | 2 +- 4 files changed, 56 insertions(+), 18 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index 756abf85..87bb794c 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -3,6 +3,7 @@ from __future__ import annotations +from collections.abc import Collection from dataclasses import dataclass from itertools import pairwise from typing import TYPE_CHECKING, Literal @@ -278,15 +279,15 @@ def _paths_too_similar(path: list, list_of_paths: list, def multiple_paths( - *, F_graph: nx.Graph, - start: tuple, - stop: tuple, + *, + start: Collection, + stop: Collection, method: _PATHFINDING_METHODS = 'dijkstra', n_paths: int = 3, min_diff: float = 0.15, ) -> list[Pathway]: - """ Calculate the Np shortest paths between two sites on the graph. + """Calculate the n_paths shortest paths between two sites on the graph. This procedure is based the algorithm by Jin Y. Yen (https://doi.org/10.1287/mnsc.17.11.712) and its implementation in NetworkX. Only paths that are different by at least min_diff are considered. @@ -294,9 +295,9 @@ def multiple_paths( ---------- F_graph : nx.Graph Graph of the free energy - start : tuple + start : Collection Coordinates of the starting point - stop: tuple + stop: Collection Coordinates of the stopping point method : str Method used to calculate the shortest path. Options are: @@ -315,9 +316,11 @@ def multiple_paths( list_of_paths: list[Pathway] List of the n_paths shortest paths between the start and stop sites """ + start = tuple(start) + stop = tuple(stop) # First compute the optimal path - best_path = optimal_path(F_graph, start, stop, method) + best_path = optimal_path(F_graph, start=start, stop=stop, method=method) list_of_paths = [best_path] @@ -343,8 +346,9 @@ def multiple_paths( def optimal_path( F_graph: nx.Graph, - start: tuple, - stop: tuple, + *, + start: Collection, + stop: Collection, method: _PATHFINDING_METHODS = 'dijkstra', ) -> Pathway: """Calculate the shortest cost-effective path using the desired method. @@ -390,8 +394,10 @@ def optimal_path( method=method) if method == 'minmax-energy': - optimal_path = _optimal_path_minmax_energy(F_graph, start, stop, - optimal_path) + optimal_path = _optimal_path_minmax_energy(F_graph, + start=start, + stop=stop, + optimal_path=optimal_path) elif method not in ('dijkstra', 'bellman-ford', 'dijkstra-exp'): raise ValueError(f'Unknown method {method}') @@ -402,6 +408,7 @@ def optimal_path( def _optimal_path_minmax_energy( F_graph: nx.Graph, + *, start: tuple[int, int, int], stop: tuple[int, int, int], optimal_path: list, @@ -452,11 +459,12 @@ def _optimal_path_minmax_energy( def find_best_perc_path(F: FreeEnergyVolume, + *, peaks: np.ndarray, percolate_x: bool = True, percolate_y: bool = False, percolate_z: bool = False) -> Pathway | None: - """Calculate the best percolating path. + """Calculate the optimal percolating path. Parameters ---------- @@ -507,8 +515,8 @@ def find_best_perc_path(F: FreeEnergyVolume, try: path = optimal_path( F_graph, - tuple(start_point), - tuple(stop_point), + start=start_point, + stop=stop_point, ) except nx.NetworkXNoPath: continue diff --git a/src/gemdat/plots/plotly/_plot3d.py b/src/gemdat/plots/plotly/_plot3d.py index 999a9cb6..6f2b6974 100644 --- a/src/gemdat/plots/plotly/_plot3d.py +++ b/src/gemdat/plots/plotly/_plot3d.py @@ -1,5 +1,6 @@ from __future__ import annotations +from collections.abc import Collection from typing import TYPE_CHECKING, Sequence import numpy as np @@ -196,7 +197,7 @@ def plot_paths( fig : go.Figure Plotly figure to add paths too """ - if isinstance(paths, list): + if isinstance(paths, Collection): optimal_path = paths[0] else: optimal_path = paths diff --git a/src/gemdat/volume.py b/src/gemdat/volume.py index dbfbc385..d45c59ca 100644 --- a/src/gemdat/volume.py +++ b/src/gemdat/volume.py @@ -22,6 +22,7 @@ from pymatgen.core import Lattice, PeriodicSite from skimage.measure._regionprops import RegionProperties + from gemdat.path import Pathway from gemdat.trajectory import Trajectory @@ -386,7 +387,9 @@ def free_energy_graph(self, **kwargs) -> nx.Graph: from .path import free_energy_graph return free_energy_graph(self.data, **kwargs) - def optimal_path(self, *args, F_graph: nx.Graph | None = None, **kwargs): + def optimal_path(self, + F_graph: nx.Graph | None = None, + **kwargs) -> Pathway: """Calculate the shortest cost-effective path using the desired method. Parameters @@ -408,10 +411,36 @@ def optimal_path(self, *args, F_graph: nx.Graph | None = None, **kwargs): if not F_graph: F_graph = self.free_energy_graph(max_energy_threshold=1e7) - path = optimal_path(F_graph, *args, **kwargs) + path = optimal_path(F_graph, **kwargs) path.dims = self.dims return path + def multiple_paths(self, + F_graph: nx.Graph | None = None, + **kwargs) -> list[Pathway]: + """Calculate the n_paths shortest paths between two sites on the graph. + + See [gemdat.path.multiple_paths][] for more info. + """ + from .path import multiple_paths + + if not F_graph: + F_graph = self.free_energy_graph(max_energy_threshold=1e7) + + paths = multiple_paths(F_graph, **kwargs) + + for path in paths: + path.dims = self.dims + return paths + + def optimal_percolating_path(self, **kwargs) -> Pathway | None: + """Calculate the optimal percolating path. + + See [gemdat.path.find_best_perc_path][] for more info. + """ + from .path import find_best_perc_path + return find_best_perc_path(self, **kwargs) + def trajectory_to_volume( trajectory: Trajectory, diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index 39be8df2..a29b80e0 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -35,7 +35,7 @@ def test_fractional_coordinates(vasp_path_vol, vasp_path): @pytest.vaspxml_available # type: ignore @pytest.mark.parametrize('start,stop,method,expected', TEST_DATA) def test_optimal_path(vasp_F_graph, start, stop, method, expected): - path = optimal_path(vasp_F_graph, start, stop, method) + path = optimal_path(vasp_F_graph, start=start, stop=stop, method=method) assert isclose(path.total_energy, expected) From 334df248d8f5e67bedec1c316b2c45c3bcbc5a22 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 14:30:44 +0200 Subject: [PATCH 14/24] Update function names --- src/gemdat/path.py | 4 ++-- src/gemdat/volume.py | 14 +++++++------- tests/integration/conftest.py | 4 ++-- tests/integration/path_test.py | 8 ++++---- 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index 87bb794c..b0de7e6f 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -278,7 +278,7 @@ def _paths_too_similar(path: list, list_of_paths: list, return False -def multiple_paths( +def optimal_n_paths( F_graph: nx.Graph, *, start: Collection, @@ -458,7 +458,7 @@ def _optimal_path_minmax_energy( return optimal_path -def find_best_perc_path(F: FreeEnergyVolume, +def optimal_percolating_path(F: FreeEnergyVolume, *, peaks: np.ndarray, percolate_x: bool = True, diff --git a/src/gemdat/volume.py b/src/gemdat/volume.py index d45c59ca..61f977cd 100644 --- a/src/gemdat/volume.py +++ b/src/gemdat/volume.py @@ -415,19 +415,19 @@ def optimal_path(self, path.dims = self.dims return path - def multiple_paths(self, + def optimal_n_paths(self, F_graph: nx.Graph | None = None, **kwargs) -> list[Pathway]: """Calculate the n_paths shortest paths between two sites on the graph. - See [gemdat.path.multiple_paths][] for more info. + See [gemdat.path.optimal_n_paths][] for more info. """ - from .path import multiple_paths + from .path import optimal_n_paths if not F_graph: F_graph = self.free_energy_graph(max_energy_threshold=1e7) - paths = multiple_paths(F_graph, **kwargs) + paths = optimal_n_paths(F_graph, **kwargs) for path in paths: path.dims = self.dims @@ -436,10 +436,10 @@ def multiple_paths(self, def optimal_percolating_path(self, **kwargs) -> Pathway | None: """Calculate the optimal percolating path. - See [gemdat.path.find_best_perc_path][] for more info. + See [gemdat.path.optimal_percolating_path][] for more info. """ - from .path import find_best_perc_path - return find_best_perc_path(self, **kwargs) + from .path import optimal_percolating_path + return optimal_percolating_path(self, **kwargs) def trajectory_to_volume( diff --git a/tests/integration/conftest.py b/tests/integration/conftest.py index 2a74c89f..d937f7a2 100644 --- a/tests/integration/conftest.py +++ b/tests/integration/conftest.py @@ -6,7 +6,7 @@ from gemdat.io import load_known_material from gemdat.jumps import Jumps -from gemdat.path import find_best_perc_path, free_energy_graph +from gemdat.path import optimal_percolating_path, free_energy_graph from gemdat.rdf import radial_distribution from gemdat.rotations import Orientations from gemdat.shape import ShapeAnalyzer @@ -112,7 +112,7 @@ def vasp_path_vol(vasp_full_traj): def vasp_path(vasp_path_vol): peaks = vasp_path_vol.find_peaks() F = vasp_path_vol.get_free_energy(temperature=650.0) - path = find_best_perc_path(F, + path = optimal_percolating_path(F, peaks=peaks, percolate_x=True, percolate_y=False, diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index a29b80e0..def16eb2 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -7,7 +7,7 @@ from numpy.testing import assert_allclose from gemdat.io import load_known_material -from gemdat.path import multiple_paths, optimal_path +from gemdat.path import optimal_n_paths, optimal_path @pytest.vaspxml_available # type: ignore @@ -40,7 +40,7 @@ def test_optimal_path(vasp_F_graph, start, stop, method, expected): @pytest.vaspxml_available # type: ignore -def test_find_best_perc_path(vasp_path): +def test_optimal_percolating_path(vasp_path): assert isclose(vasp_path.total_energy, 11.488013690080908) assert vasp_path.start_site == (11, 9, 6) @@ -69,8 +69,8 @@ def test_nearest_structure_reference(vasp_path_vol, vasp_path): @pytest.vaspxml_available # type: ignore -def test_multiple_paths(vasp_F_graph): - paths = multiple_paths(F_graph=vasp_F_graph, +def test_optimal_n_paths(vasp_F_graph): + paths = optimal_n_paths(F_graph=vasp_F_graph, start=(10, 4, 13), stop=(21, 3, 10), n_paths=3, From 2c044c4cd4bfca05f85eee2812667cb3a94318e6 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 14:31:16 +0200 Subject: [PATCH 15/24] Simplify percolate argument --- src/gemdat/path.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index b0de7e6f..b2383875 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -458,12 +458,12 @@ def _optimal_path_minmax_energy( return optimal_path -def optimal_percolating_path(F: FreeEnergyVolume, - *, - peaks: np.ndarray, - percolate_x: bool = True, - percolate_y: bool = False, - percolate_z: bool = False) -> Pathway | None: +def optimal_percolating_path( + F: FreeEnergyVolume, + *, + peaks: np.ndarray, + percolate: str, +) -> Pathway | None: """Calculate the optimal percolating path. Parameters @@ -472,19 +472,20 @@ def optimal_percolating_path(F: FreeEnergyVolume, Energy grid that will be used to calculate the shortest path peaks : np.ndarray List of the peaks that correspond to high probability regions - percolate_x : bool - If True, consider paths that percolate along the x dimension - percolate_y : bool - If True, consider paths that percolate along the y dimension - percolate_z : bool - If True, consider paths that percolate along the z dimension + percolate : str + Directions to percolate, e.g. 'x' to consider paths that + percolate along the x dimension, 'yz' for the y/z dimension, + or any other combinition of 'x', 'y', and 'z'. Returns ------- best_percolating_path: Pathway Optimal path that percolates the graph in the specified directions """ - # Find percolation using virtual images along the required dimensions + percolate_x = 'x' in percolate + percolate_y = 'y' in percolate + percolate_z = 'z' in percolate + if not any([percolate_x, percolate_y, percolate_z]): raise ValueError('percolation is not defined') From b5dd84623d24d3f7bf8d87cdd26c65b7af9dd08f Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 14:31:49 +0200 Subject: [PATCH 16/24] Fix plot and get rid of Volume requirement for paths --- src/gemdat/plots/plotly/_plot3d.py | 63 ++++++++++++++---------------- 1 file changed, 30 insertions(+), 33 deletions(-) diff --git a/src/gemdat/plots/plotly/_plot3d.py b/src/gemdat/plots/plotly/_plot3d.py index 6f2b6974..6ed98365 100644 --- a/src/gemdat/plots/plotly/_plot3d.py +++ b/src/gemdat/plots/plotly/_plot3d.py @@ -181,28 +181,29 @@ def plot_volume( def plot_paths( - paths: Pathway | list[Pathway], + paths: Pathway | Collection[Pathway], *, - volume: Volume, + lattice: Lattice, fig: go.Figure, ): """Ploth paths over free energy. Arguments --------- - paths : Pathway | list[Pathway] + paths : Pathway | Collection[Pathway] Pathway object containing the energy along the path, or list of Pathways - volume : Volume - Input volume to create the landscape + lattice : Lattice + Input lattice for coordinate transformations fig : go.Figure Plotly figure to add paths too """ - if isinstance(paths, Collection): - optimal_path = paths[0] - else: - optimal_path = paths + if not isinstance(paths, Collection): + paths = (paths, ) + + path, *other_paths = paths - x_path, y_path, z_path = volume.voxel_to_cart_coords(optimal_path.sites).T + x_path, y_path, z_path = lattice.get_cartesian_coords( + path.frac_sites(wrapped=True)).T fig.add_trace( go.Scatter3d( @@ -221,25 +222,25 @@ def plot_paths( )) # If available, plot the other pathways - if isinstance(paths, list): - for idx, path in enumerate(paths[1:]): - x_path, y_path, z_path = volume.voxel_to_cart_coords(path.sites).T + for idx, path in enumerate(other_paths): + x_path, y_path, z_path = lattice.get_cartesian_coords( + path.frac_sites(wrapped=True)).T - fig.add_trace( - go.Scatter3d( - x=x_path, - y=y_path, - z=z_path, - mode='markers+lines', - line={'width': 3}, - marker={ - 'size': 5, - #'color': color, - 'symbol': 'circle', - 'opacity': 0.9 - }, - name=f'Alternative {idx+1}', - )) + fig.add_trace( + go.Scatter3d( + x=x_path, + y=y_path, + z=z_path, + mode='markers+lines', + line={'width': 3}, + marker={ + 'size': 5, + #'color': color, + 'symbol': 'circle', + 'opacity': 0.9 + }, + name=f'Alternative {idx+1}', + )) def plot_jumps(jumps: Jumps, *, fig: go.Figure): @@ -374,11 +375,7 @@ def plot_3d(*, plot_structure(structure=structure, lattice=lattice, fig=fig) if paths: - # TODO: Does this need the volume? - # Revise after https://github.com/GEMDAT-repos/GEMDAT/pull/282 - if not volume: - raise NotImplementedError - plot_paths(paths=paths, volume=volume, fig=fig) + plot_paths(paths=paths, lattice=lattice, fig=fig) if jumps: plot_jumps(jumps=jumps, fig=fig) From cea568e0820ce42ed343fdc62b3931aa66d7f123 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 14:43:37 +0200 Subject: [PATCH 17/24] Fix tests --- tests/integration/conftest.py | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/tests/integration/conftest.py b/tests/integration/conftest.py index d937f7a2..ff2a2dd7 100644 --- a/tests/integration/conftest.py +++ b/tests/integration/conftest.py @@ -6,7 +6,6 @@ from gemdat.io import load_known_material from gemdat.jumps import Jumps -from gemdat.path import optimal_percolating_path, free_energy_graph from gemdat.rdf import radial_distribution from gemdat.rotations import Orientations from gemdat.shape import ShapeAnalyzer @@ -111,21 +110,15 @@ def vasp_path_vol(vasp_full_traj): @pytest.fixture(scope='module') def vasp_path(vasp_path_vol): peaks = vasp_path_vol.find_peaks() - F = vasp_path_vol.get_free_energy(temperature=650.0) - path = optimal_percolating_path(F, - peaks=peaks, - percolate_x=True, - percolate_y=False, - percolate_z=False) + free_energy = vasp_path_vol.get_free_energy(temperature=650.0) + path = free_energy.optimal_percolating_path(peaks=peaks, percolate='x') return path @pytest.fixture(scope='module') def vasp_F_graph(vasp_path_vol): - F = vasp_path_vol.get_free_energy(temperature=650.0) - F_graph = free_energy_graph(F, max_energy_threshold=1e7, diagonal=True) - - return F_graph + free_energy = vasp_path_vol.get_free_energy(temperature=650.0) + return free_energy.free_energy_graph() @pytest.fixture(scope='module') From ce915e947ac151d7b6fea959524296faf74c23bc Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 15:12:25 +0200 Subject: [PATCH 18/24] Refactor plotting code --- src/gemdat/plots/plotly/_plot3d.py | 34 ++++++++++++------------------ 1 file changed, 13 insertions(+), 21 deletions(-) diff --git a/src/gemdat/plots/plotly/_plot3d.py b/src/gemdat/plots/plotly/_plot3d.py index 6ed98365..64999492 100644 --- a/src/gemdat/plots/plotly/_plot3d.py +++ b/src/gemdat/plots/plotly/_plot3d.py @@ -205,24 +205,16 @@ def plot_paths( x_path, y_path, z_path = lattice.get_cartesian_coords( path.frac_sites(wrapped=True)).T - fig.add_trace( - go.Scatter3d( - x=x_path, - y=y_path, - z=z_path, - mode='markers+lines', - line={'width': 3}, - marker={ - 'size': 6, - 'color': 'teal', - 'symbol': 'circle', - 'opacity': 0.9 - }, - name='Optimal path', - )) - - # If available, plot the other pathways - for idx, path in enumerate(other_paths): + for idx, path in enumerate(paths): + if idx == 0: + name = 'Optimal path' + size = 6 + color = 'teal' + else: + name = f'Alternative {idx+1}' + size = 5 + color = None + x_path, y_path, z_path = lattice.get_cartesian_coords( path.frac_sites(wrapped=True)).T @@ -234,12 +226,12 @@ def plot_paths( mode='markers+lines', line={'width': 3}, marker={ - 'size': 5, - #'color': color, + 'size': size, + 'color': color, 'symbol': 'circle', 'opacity': 0.9 }, - name=f'Alternative {idx+1}', + name=name, )) From 0e892c7451639720c0c57459dd00ad3c22d5c92f Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Tue, 16 Apr 2024 15:24:30 +0200 Subject: [PATCH 19/24] Yapf it --- src/gemdat/volume.py | 4 ++-- tests/integration/path_test.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/gemdat/volume.py b/src/gemdat/volume.py index 61f977cd..518c79dc 100644 --- a/src/gemdat/volume.py +++ b/src/gemdat/volume.py @@ -416,8 +416,8 @@ def optimal_path(self, return path def optimal_n_paths(self, - F_graph: nx.Graph | None = None, - **kwargs) -> list[Pathway]: + F_graph: nx.Graph | None = None, + **kwargs) -> list[Pathway]: """Calculate the n_paths shortest paths between two sites on the graph. See [gemdat.path.optimal_n_paths][] for more info. diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index def16eb2..4d87847f 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -71,10 +71,10 @@ def test_nearest_structure_reference(vasp_path_vol, vasp_path): @pytest.vaspxml_available # type: ignore def test_optimal_n_paths(vasp_F_graph): paths = optimal_n_paths(F_graph=vasp_F_graph, - start=(10, 4, 13), - stop=(21, 3, 10), - n_paths=3, - min_diff=0.1) + start=(10, 4, 13), + stop=(21, 3, 10), + n_paths=3, + min_diff=0.1) assert len(paths) == 3 assert isclose(sum(paths[-1].energy), 5.351758190646607) From 06e00d2520bffa9a37d302e465c3b3751b7eb68e Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 22 Apr 2024 10:33:45 +0200 Subject: [PATCH 20/24] Update src/gemdat/plots/matplotlib/_paths.py Co-authored-by: SCiarella <58949181+SCiarella@users.noreply.github.com> --- src/gemdat/plots/matplotlib/_paths.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gemdat/plots/matplotlib/_paths.py b/src/gemdat/plots/matplotlib/_paths.py index 911e69e7..97025b4c 100644 --- a/src/gemdat/plots/matplotlib/_paths.py +++ b/src/gemdat/plots/matplotlib/_paths.py @@ -16,8 +16,8 @@ def energy_along_path( Parameters ---------- - path : Pathway | list[Pathway] - Pathway object containing the energy along the path, or list of Pathways + path : Pathway + Pathway object containing the energy along the path structure : Structure Structure object to get the site information other_paths : Pathway | list[Pathway] From 03a2ea3f00ee0ad391624788be616ad08f2c6f2e Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 22 Apr 2024 11:06:28 +0200 Subject: [PATCH 21/24] Always wrap for frac_sites --- src/gemdat/path.py | 14 +++++--------- src/gemdat/plots/plotly/_plot3d.py | 5 ++--- tests/integration/path_test.py | 2 +- 3 files changed, 8 insertions(+), 13 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index b2383875..b8669e69 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -88,14 +88,10 @@ def wrapped_sites(self) -> list[tuple[int, int, int]]: xdim, ydim, zdim = self.dims return [(x % xdim, y % xdim, z % xdim) for x, y, z in self.sites] - def frac_sites(self, wrapped: bool = False) -> np.ndarray: - """Return fractional sites. + def frac_sites(self) -> np.ndarray: + """Return fractional site coordinates. - Parameters - ---------- - wrapped : bool - If True, wrap coordinates to bounding box using - `.wrapped_sites()`. + Note that these wrap around the periodic boundary conditions. Returns ------- @@ -105,7 +101,7 @@ def frac_sites(self, wrapped: bool = False) -> np.ndarray: if not self.dims: raise AttributeError( f'Dimensions are needed for this method {self.dims=}') - sites = self.wrapped_sites() if wrapped else self.sites + sites = self.wrapped_sites() return (np.array(sites) + 0.5) / np.array(self.dims) def path_over_structure( @@ -124,7 +120,7 @@ def path_over_structure( nearest_sites: list[PeriodicSite] List closest sites of the reference structure """ - frac_sites = np.array(self.frac_sites(wrapped=True)) + frac_sites = np.array(self.frac_sites()) nearest_structure_tree, nearest_structure_map = nearest_structure_reference( structure) diff --git a/src/gemdat/plots/plotly/_plot3d.py b/src/gemdat/plots/plotly/_plot3d.py index 64999492..a919d0c9 100644 --- a/src/gemdat/plots/plotly/_plot3d.py +++ b/src/gemdat/plots/plotly/_plot3d.py @@ -202,8 +202,7 @@ def plot_paths( path, *other_paths = paths - x_path, y_path, z_path = lattice.get_cartesian_coords( - path.frac_sites(wrapped=True)).T + x_path, y_path, z_path = lattice.get_cartesian_coords(path.frac_sites()).T for idx, path in enumerate(paths): if idx == 0: @@ -216,7 +215,7 @@ def plot_paths( color = None x_path, y_path, z_path = lattice.get_cartesian_coords( - path.frac_sites(wrapped=True)).T + path.frac_sites()).T fig.add_trace( go.Scatter3d( diff --git a/tests/integration/path_test.py b/tests/integration/path_test.py index 4d87847f..380c90a7 100644 --- a/tests/integration/path_test.py +++ b/tests/integration/path_test.py @@ -12,7 +12,7 @@ @pytest.vaspxml_available # type: ignore def test_fractional_coordinates(vasp_path_vol, vasp_path): - frac_sites = vasp_path.frac_sites(wrapped=True) + frac_sites = vasp_path.frac_sites() assert isclose(frac_sites[-1][0], 0.4107142857142857) assert isclose(frac_sites[19][1], 0.6785714285714286) From b76956e3cf186614aa0b6b657726241bab767519 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 22 Apr 2024 11:06:41 +0200 Subject: [PATCH 22/24] Rename file for consistency --- tests/integration/{rotations_test.py => orientations_test.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename tests/integration/{rotations_test.py => orientations_test.py} (100%) diff --git a/tests/integration/rotations_test.py b/tests/integration/orientations_test.py similarity index 100% rename from tests/integration/rotations_test.py rename to tests/integration/orientations_test.py From d8fc0992d1c67def42453b501471fd3bd574f729 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 22 Apr 2024 11:14:31 +0200 Subject: [PATCH 23/24] Clean percolate function --- src/gemdat/path.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/src/gemdat/path.py b/src/gemdat/path.py index b8669e69..6e2665e5 100644 --- a/src/gemdat/path.py +++ b/src/gemdat/path.py @@ -478,33 +478,25 @@ def optimal_percolating_path( best_percolating_path: Pathway Optimal path that percolates the graph in the specified directions """ - percolate_x = 'x' in percolate - percolate_y = 'y' in percolate - percolate_z = 'z' in percolate + percolate_xyz = np.array([dim in percolate for dim in 'xyz']) - if not any([percolate_x, percolate_y, percolate_z]): + if not percolate_xyz.any(): raise ValueError('percolation is not defined') # Tile the grind in the percolation directions - F_data_periodic = np.tile( - F.data, (1 + percolate_x, 1 + percolate_y, 1 + percolate_z)) + F_data_periodic = np.tile(F.data, tuple(1 + percolate_xyz)) # Get F on a graph - F_graph = free_energy_graph(F_data_periodic, - max_energy_threshold=1e7, - diagonal=True) + F_graph = free_energy_graph(F_data_periodic, max_energy_threshold=1e7) # reaching the percolating image - image = tuple( - x * px - for x, px in zip(F.dims, (percolate_x, percolate_y, percolate_z))) + image = F.dims * percolate_xyz # Find the lowest cost path that percolates along the x dimension best_cost = float('inf') best_path = None for start_point in peaks: - # Get the stop point which is a periodic image of the peak stop_point = start_point + image From 12559b4f1f2e1543f5366e44f147f74425ad78a4 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Mon, 22 Apr 2024 11:22:03 +0200 Subject: [PATCH 24/24] Remove whitespace --- src/gemdat/plots/matplotlib/_paths.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gemdat/plots/matplotlib/_paths.py b/src/gemdat/plots/matplotlib/_paths.py index 97025b4c..bb40fa0a 100644 --- a/src/gemdat/plots/matplotlib/_paths.py +++ b/src/gemdat/plots/matplotlib/_paths.py @@ -16,7 +16,7 @@ def energy_along_path( Parameters ---------- - path : Pathway + path : Pathway Pathway object containing the energy along the path structure : Structure Structure object to get the site information