From 31c0b24317222e6e9fdcf41927ae0ec60eb3ea1a Mon Sep 17 00:00:00 2001 From: wiederm Date: Thu, 29 Oct 2020 13:30:29 +0100 Subject: [PATCH 01/11] updating scripts to run sampling for tautomers given two input smiles --- data/test_data/exp_results.pickle | Bin 0 -> 261676 bytes neutromeratio/analysis.py | 101 +++++++++++++++++- neutromeratio/parameter_gradients.py | 9 +- neutromeratio/plotting.py | 6 +- neutromeratio/tests/test_neutromeratio.py | 56 +++++++++- neutromeratio/utils.py | 13 +++ ...ing.py => Generate_equilibrium_samples.py} | 0 ...ilibrium_samples_for_new_tautomer_pairs.py | 98 +++++++++++++++++ .../generate_equilibrium_samples_in_vacuum.sh | 30 ++++++ ...amples_in_vacuum_for_new_tautomer_pairs.sh | 30 ++++++ ...gradient-opt.py => parameter_opt-kfold.py} | 0 11 files changed, 333 insertions(+), 10 deletions(-) create mode 100644 data/test_data/exp_results.pickle rename scripts/{Generate_equilibrium_sampling.py => Generate_equilibrium_samples.py} (100%) create mode 100644 scripts/Generate_equilibrium_samples_for_new_tautomer_pairs.py create mode 100644 scripts/generate_equilibrium_samples_in_vacuum.sh create mode 100644 scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh rename scripts/{parameter_gradient-opt.py => parameter_opt-kfold.py} (100%) diff --git a/data/test_data/exp_results.pickle b/data/test_data/exp_results.pickle new file mode 100644 index 0000000000000000000000000000000000000000..7b2ec7c500bdd39e772a3b8d5d63d1f3b07dd949 GIT binary patch literal 261676 zcmd3P1)Lnk_H}Ug;1+y?WoOobkR=ODg6sl80!aw41PHQWZg3df-QArexN8=7cY?dS z>vzuWx?NS>vmx-_`+a}j>rS6rw|i#ioVs;Oy8FOsD&<}^LucyMtJjb<*56>gUB>KO zEKgM_PhB%~8rVjT-nKkV=Qd-9^5~?N(T6>@K5s9yPW+!_eNiP)A##e_LDswjwtEzu3RMP@XYy zEt>00L+7Jw|F%Lq9NKCU-(p|KY~6?_Hjy)V;uV`k#JH8MR$`)}gba zT>psTD#bNM76+DR8#*maU0d%msyutIIO%t_E%Ui`P1}Qi2pSF4hv;|hlP`$ zOn-;)C&RC>Ncbx(TI!xU^D8VC{tESLjI1A6Ufll*eKLQ_68=|MGW`{n3V(&Acfj92 z8-54~W6jMr8nfrfLTznrd6`OiS^pa>r@ujs{suE9N*H%Yxo`G2SYAK9qqES}+1b@u zD6f$Dvscvbp3vXFt(~N=6#ls@m%69D>;4P=*y@v%^nOF9!*^?LDD)XvuJwOsVd(6^ zS8kU=vDjQL`aiZV{bTFHKej>5yLatWZd5ZAEIW-YHx1nkR~yzbXqCcL8~IsfBMx!`HgEU56J%XtwZO)h4rT| zcJjX>30NaTz`&vN(2Kh3tWmGNokiuLE#zU%Qul0o&3M>M6VFR|Sj+KH1P^O_9!hy0 zTFQew4}()4)(Lr7H|1fykcah$9$TFU{Io@~(X12{YLkCR)}U2du*Yv+{9A@Pag~~$ zuTt#mt`m{h9+1}&lea;hykdF7N_its-Vh~kUOWvkK7scP*)?U+D=&Wh0 zX)m_*?cdQCm8b^$LSKZ+j`C(1Du)j3lbnMnS_2BPEAT9GThQ|niaN@hhrDi48u0bV zC05$AaAL~qu%R>HGR+Njeb!iUV0pM_c0`WZy4v!Vp4pKpvs;DCZk;l_O~`EL(5LY+ zs2z3Wuzi(!G6B{K*n`d*R3(zrj+#_duf`F5`=ge%jOYtq3&n)+dQx4egCnbILsC{4 zQRu5RvB(n&BTzf}F{rJ)nq9dq`IvTIOKWBU7lV-=sK5nYMbF1oQLbN&U$fx8bNX$x z&L^K#-Zl`&s8}4^<;77~-o8@a!HZ+GietyHW?yx4d8e#6c20vT1c5reOGXG?KB!U# z460*7(swO&FZ{x1XV1UOq?GistIYMN%;G5@o2R@`-py0KdrJ8pA?16fl zJ}CRsA3Stk*0Et-D~u>0qB3A&ua&ZVNQQ<(hc2kT?JJDXn%L$8tX9OsLQW1ZP5knN z@)G0M9h-7;L{bqOYx``y;x^?YJu64$St*u}_N*L}vT|(5%5fA9uuC6D`dyE_K=Pbodm#nuLmt>=46FUV6`SH93w zdQnR0#UZ7aXu`Ty`O=Wk%Z46JLUXx`RWml3Ml733jdcMTj@Lr5YF3k0Kz|txinvTt ze~nUWi;LyU1LCfTiMujSTx0pFN_m_oZoCpVAtY{Gzw*^t;;u2`jJ9G|SD}1uhPdl| zcLT)v-UjmXLizfTy&Fp1bAI~bz7PMjZ_3_{%6fA{OP_(OwXe5&p|&Y$BX9CN-kj&L zzI==4@z#{b+d>|1S2MHFJ3=P!9C|{QN&L*(=DCXbH|Y=JD$~~L`Ke3Pf(`goqu=JFGs z&?i$up9%?mIwka(kkDs`jyFQtnb0Hy3~9pJn1o49LKsh7Cj#LUKYa&SO0p4z*iK|j zOm3~Fp?8v0)jD#PToh;Hsw7UY^!HcgnY=9WAxAy(qC)w(fb8dEvR}xNT~}Xzu~L4? zll`)i{R+sQ)O#>`iOqj~c449XYL@KRlDez)kvVI-wvUrCB)^_`CDj&giPsw;vu~EV z7hGt&xeqz9O@Dum5xkW!+uW=n_ifMYJJ~)izWuwN*!NOm-w%oXV24YCU*tn|+GXp# z%O8a_e>`*tM>De&Z2(cFnv)vEG?bN?ts(PKIN*m)GgUH8F4$^TJ1NygdQ$n5JP{4$ zPb=lmJQ1HO5nm_~>o_7-Eq|FM;wvK}VFAe_dkkM^i1?-|5#NSHd{^pT{PI`tS?}ef zyGJ6vPl;$lzv2f^#E&^58p=O;B7RPZ_$4Ib*NBMUQX+m2iTEQ=gjTEMw@2{AFPx^3 zIU=k>uSC>`>@~B8;VY=J*EFL1XP$`0@?RC{g<++X3z-UtbjZ{=h!`0x)4-G=MW%(# z4JBzF8FbW1k4(pMj^^!VqiQmJdJ1Ywl6DsQcA^`YR2rFq<$i2UW(3Y4(kv@_V; z(g&o;%rw~en#@9!87Hc1lv&M*BW*Gp5r-8q*e0_R@MDwlHkkvsbimIr;6GojWKKH? zHB08QlaiqsnVa}RoWKl+6vhc=yN{-w|0utVvYi?PS?!UgSxu78SiiHD0$>(ZpbT27 z$f0?J%!6;byY|I`k+gd04Rb~2#olXTJ|NY^{5a6Iai(NYZNk2?08Cj;EC`!YgI+$n z_OcKSu?Q9pd^idqi_j3Tz9?|1b-82vKk;eJUrG|{sVqjXsqyLP6_Wr6cUb6 zw)$FG!km=wC5epir3eDTm$pKrIzt}CAwYLctniXxn3 zIbMRpgF3b*HH6{tn)HPwM@fS$4|7FUz}{1`B9KzD5)QPQuAT(W3ucm)Vaig{54NGR z;IowyM@ub@F)f9RpHfq#F`%Xn7}V^))|S`2oBW)lrk>WczNEY~@MMmJ2Qq7rMq25k zv@{VJEzJZ0EiJ(Qy8icD#ZF>9T-8o;bgf3Lqoi(d1+i2|SHKqsV%V%Pa*|4yx;XaG z0#@}6$*dnC{qX_0fZEU~tHWNA0oZ%8T7i_THE^J0t>wtNT?WFGC94fKl`_d<#-YiC ztVwfB*IG1jXt@o!=&-IqBiP&TEXdll2aJ_~O9#S_|D}~_5WN^eTMF7&kij&h%&kKd z)+h|M<1bv-oRqrth>W`R2?FZcfvZv1VJDHg4eTUG-G;=&@vW2|%~6;qwhbz#c?o(B z*1n~QhNDmWhm3sMq3qJJ?q~=AaTn>0Y=lqA^VcLpV6Vu=*n9pq0aE@p#ev2@n^fqO z&0xy%$9``W{x+vM=5LF@vz=*;<}@;l=77H8z@_f#KRfl+*1gV1=^H^WPhXpCNkdBC zNTTetLvxF4WnRkP)=r3Sw_QNksY!3RPF?%RPKy}h{{p23rtxmyI`wU!N$-YbGa)`tS%he z_FR@}4!9KHJeOnXWn2ziu}yZPA?0#+qAZs!wX%nKDVKW^8JBwz1YGWIh2(M{D|B4$ zOB~j@Bzu^A9r8793YpxGpW>OU5JpVy4@;iO7C8X+iX4c&XYwE*W%6JgL`;4mhrpC& z@=(}}$z-O2@vq<*94PC(~vTG z3{j5Bf*fmJ%H(lG#^mt?0h1?KA(=eU3LTRt5$BoIsE3>|GZu=U0L>vg=xhqpSc~8~ z`8s(rKgBb73Sq?Lsj%cynMOGc=8Bw-y(jVvASLol94L{NT^=N7!IUNPY}gVana1G% zE}J=r_L#_X^RAlEoJV`W_O2UZCt6<61 z%Enq52YW@vWADkF0HkDIjf1FGJ}KA0lqK_8*sNBD3Cwl0$7Eg~xP}$@2HFEMZv-w) z{I&4ZwaYBBU`pmq^s;I>aK$y{W*X9Jc?(glS{CJ2^HNrCBQjQRCkR-*!wSjjomS{r zy^DB@fYpB%U=#1jS7N!FpW%tThcF^?A}o0#>!lm^irkC6C-OcZCGvh8fXGR`OD<>q zOdfzKOXP#Fr9^h2R=P~|A=+aiA5LAhp+pZc$$B24Js|Q?U@Rz!XDQ?{dL=|Q^}hy*z17%HvZ+#^cij0gulBhr>7Wta&qAea^fcvCk9lXvFHY1E-0T zmNV^gpytNGH_1X>+|_Xc4s^S?i_?;MT%A?UXGAaHGjqYUNnV7xA}?X@X?z(-X?z6- zQE=T{UWF-3<7==Dof&@)r;pnBUqL3(7E}0o>V`oc=k^9|0ex=*ms*!xW04tm8~A8i z58tAf>%hOwBmamDQ5JoleQ&@5}v8gDZ!CsNi zvG){y0i+aui36>P)_c5GzJe)B;n%QfhZ}vuPE^D(SXt8{-_RP9_$^J0d8~eHL_3mx z=6AFQM1Buk8u08MM?JQ{-P@-`{y;CShsB0I%yWLEA!YL?qL5A0L;Qh1o0GEn3z4z; zD?z~KZ@|9!^}putb`q=OA9j*s`A_2QqUxwK#!*ci!*x?#$c1zqtw4I4Y#oLhvqR^j zEjovyVg1BOZR|iMvmWvnJ}*}lo23^NUXiJ=_dHGwq&!Z8gQzNwl4)Ve@;DuAL+4Ja zVht-IrkC7Qu1wFen9Uhznv?EHmGsUtvLxVhCg4)*Qr{M~zG%`BDW5ab3mZq)X5n#4 z<*Y;@m8_TzGMhOmm9rBWm2(gTRL*IIteA6Ip`&bW;!{bPub9@|4{M~3V|IcR&glos z0J1a0m<{Z?%DTF$b6UDc3)j;H1QAx|Oj!vm44X;-YpuBUG#O#9qqB(br9br8d z$Ke-aSs;SNflJ-D`dTihLrVYq8#f@ z(qdl9`YJ@m`l#JEIS?_O!j`h`vkC`m%{Xh{Cq<)ORj&~3^9I(h4;fO((RBwHRWnLUlFShF*gI18;%;cxocQ=zu z1|fhwXbJ?dCvd5C*?phge!JEyQUUBmFAAWkcEvu0LW}H8Ln?rMh_b0*qwH%=%Kv^u z#(#w%;D3KBB>x9kq2vER;*FJm(1YI_={9@9**@Ao=$r+;vV(YaPuszS5p9RSlBKOs zD~G~dk;Aa}v>gtlv>kzih_*xINSLy;9R*vG$t6o3kj!IobcfE4AI-A3E+3P+cW{A3 z;@yV1Jxr7z%d&vTkP}>K?Wqsr_T2eJ!PLvJlbuDs*d8r((Br-W%MG(kg954pV zXmJ`CPmgn{0sj6PXq!oyE=h6AveUII(IRVh^4JtZ7fAKo`m;Jxuf}N3YijUj zmkuB18iC}Z7>To$P|47mt|HW{A=970Z|TKwHDM%%YhcNXp;4}dy&~6P@5OLEkc#03 z96$_{dbhdp(!J$In6hHH2{v7>k2sDMq3ELZ#iU`)MQCtzWKKjT&w z2XeR#xYU{*J-wY?+Q2syb=A}zG^A>{lPHW0BMTTlmAlMI#c(%~iQyiCKnxRstBu-p z+esX?xz|o|dby8yoa)6b@MN}61;lf~N+zYDH9Q>@=csGi@Z%PlKnT#28w`LVv=33` zswTL9a#gt>-#IUi7I^^niadzD7so?DDvpP7ppnA*{cp%4FlEK@C~V695WP1M4S+mm zCf~?EPLtka%7ibuH?<&7&>HCAN#NYVQ+bMBs)NRQu055fX-IYO3{f^J)Hlkr=A}Y- zj>v@YJV79Y7l<)B)~iikv?53MOT?R#?ja1%sKadLyI9>_1h?RjjLM*%gG6R*mJ+zw zkoN3#SO<}Sfu3POUgk&my8Q}aRJUJ+CD+y$dK}`}_35 z#`gaM9;e)XNEFs_4(GMVN9Lr|eoSQ4enJpX`>7SO{r}7g9l4(qUq*7Z8Am*J_p<$c zv#r_r3J?Q-C#b^dzi9H`2jYi0`;8t=+w0B+_$SBKhhja<0qOpx;sSc zIvp_{$d295ED7ZB3vg-T*Jn)q;kN4^sT6W?-mml`pDm3I@*6>_hu?{Eb-hXcFfSFv zpLPK*mHmq#5JWFrInDs4f+d;7MitleE*i15daYF-gS|U5OVbFD)O1^AuDvyv1(=(w zVbf*j z8n5lfJ4EeSSP(NgD@`1zQ6`z)o34BAO%RYdFYu1pkUbxrj>j2QLU#4;Gvm76=E3`BHA6Gfs|R!T z!cn0utB*(Glis{u$kWa^KR>^hr+op!i1r0x$t#9fQqpUPuqo+%2?ElW2hNec0-f|0_nzBn=3bc;JYn7usC6Gq zaZYA(Ex5l&794=|;-_sJUkYM%Uogc12|u;j^ak~Lwk$XeKY^4A7Z@=G{~n(LQk5KLL}2g9Zu zCuuq6)ny$s`Lujpn!>c4_xQAq!%{u*C5*xr2LH{%JJ7E+nJa0w>=RyN(aFf6$)b5Ra~xgv*R?-@G`NEtgE2g;bO>0MWjfGNw^k+7wsVO+Pw|Iqv8!PtTv z#p1ZK9!=8%_yIdON$wocCmL=i%Q5r_m^>D^)Vjo3d%Q8>%Xd;HkE0j5Grf>XAJ5~I z$PKfX7p;kp1LSt z1uTdKb0JM>0*>mA)*yK>hFrv=Kq(gk=WjT@gkIQC=bPkGf>bJ(5#{Q9i(GDAs+B8< zu&GwABnY%}6>zT3kE4@5FC*7BTb=iLbJE%JoFGFXFED0v2Kg;im zGM?Yb6Fz}3BK&Gt@`TsPHLzFYTI@aH*8wTv*W&<$PwH*=f3GVyz?3EYM%a{aJg%k- zkB#E`MzrQbu=ImQis9#SpoLg#qQa0++g%$xNQyMlWo=TAFL+ zc7l}hJBY&QfRS>{gx_gS%J^MG*p%_R2?ECN0ZthAzve_cfurQzc7lr-_Y$Y8qjU)y zZ;b0w#yHr9r!au*LVgmr8#&p&U;c(ca zw8z{&MiYm_+?qIY%c<@m@;FNZVxIsmbuUqS?ZsbjQC2S4em_YsY)b4?1Szpk6J?35 zmuJjLiG7v`n-cpRK|t*Dz*%BnuoFn^i*|yIhZDzZ$ax2&r;;p;(cy6&991ewZqmr7 zjHZv%NEkIC@aztr9TNF6J~K~bgS-NJMP9|;6ZsmD5;+M65s?Grb(pe5z5!cEBo?nI zk#Eu-6ZsZR86t-$k#DmkAo3mHQg@%s(zSQ#g-wZkk02%TeWEOp4f26GDUlx%VN)VM zA_$257&uGhCw2mf{M1fxM1DqGjYzv2hg7g();d|k?lo>d29f?AEUlW3NLbk|zdmg_6QXap?fi|y}mDR~NFlBlC7Pg+UvM~kuj>R#Z z-_sOlWn=I+1X{mkQ()nAB0Rxx|gP%pol zld}375jJJ@cY=V`Kdg{h*`HSE82pRaWo0(qlw@T(_MuMzoU*JeT34J*PJ`Pqj%r<( zz^n|D9n8um^_rDU;pa;){9ZvLJr!ZpNKXw*PBry14a^mp7JILn>3~!<)8imi z&G<<&158=f%m`cZL=fIoVClH#@<`DNG84;V+00B6pGTJ_^Yn~#GYee;-OLJHYRxp;>`n+bd*6!bl-2!IFzJ z4YD%K73qh)7f3CT3Z#I8FwRVvMv5?H1yToF^*ECZ`P8#KmPi9lVVqIVD9$v}C6GxI zaH(~^%%c~Y>E+@~3y-lS0$GJ98)q72RdZ5(tVU$|=uZ&nV|6QJoEcz+F3z+P|Ci%T z7yAQ@K&BUhzl$?Hon$iM7(%AC>M_o&!9U8YXdq#vqBdA^ai*bG)`T5{`q+CJtqr6y zD&at7o!WGu*&>qWY-PBbhjJt7lh^$9@pqcf7OWiYO zR%Eu*%Y~LY4J{qCU`v#<0Z}%zv=n7S^HS+-L}bz#LJ&x2W8hq9*@RAUQc7ZU0D)m7 z$u0dIsxHiE~-6 z%-h9O(J2IhicYmcs^~N;bSgTX z_)&&~|3-HtFNox+u#?LYCVS4m(a;(Ele~t`B#bn47A(0Y1vlcJ4Lg?ZWA8?>X-^mAUT7viqj(Wb3?G~^P8Z}}Olu&OOMo#Cbyw#%YxEhO zhMr65m1GBvjk+TDGFq@Dg1MY1lO5C-8{`V}QpH?JWQw_pAW+OWE2Lt^TcK0T1ma;- z3~!4}`yA+fYHpb1Zb@U3%nT$$R=PwR?S(cd0BlIOZr?Q_1Uk9J7Nj6o^Fus`*APY= zUJFZ}!=hXVdqu9t-g9^ZkaBn<4k8Zcm78G7a(FXrIvz z-;c!hOPRfm9s#qr1D9Ge&-j%)=oK-0CoR|#X73`(G21A2o0l?s50NoDksx5U+X~6- zy;kU$y^r{MGTRA$nD6mlH=&q3Dik=OStk{k9V)LbPpGDWq{onGV9GCx`*)M{f^yND zd8_EfJ&Bllxu0L%>*WE$NG}h(75txqlN#xqk}>5%=TeZJ4s$zXO{e`Q~bC{67_eyvyQP0q@b|p7V}3 z+$3kjQjquQ6X@Ur;8JVm6;kpcy&@faL<_b=2OksVbkHoHn3w9{QzFyBX9R%`KDR>Z z;0r5sI{1?K66&A>q$7|+3+y9usG2Lau`b^hY(pCqgm(1P0Rgun&;=NaE4mcnZ>~qF z=z5aYK~j_mm%D%Btii*pfVvBi-aYr<9-R5-a5wnij+F%R9fqvodP=m41O*egiHI z$ll-gJG~;c{6P!0L@j?3<u6r5Ap}R4r2x1ZtVu3aOT9tk9`tTH+fMwItl> zw4;dy(_2+U!wdoefB_|m9C zL9??k*3BF=F&}9gl5Jz$NpsR8P|RGwrHS8UhPdaZSEQJEXh{{*n<%T8x>}jnyi_ss z5m|hhpCC}o0#-=HENF#JF$)o&LB+70!N3AXjO=E=WI($MOfpM}#gPOS-Z@;umd4#5 zRI~=iu0$c|7zScVA_QCv%=fg$ID0u+m|qZ^Q^g{Lkt!C2C9jGmSq%1yERMZbMIRtl z#S%D(RB@av2~$=TOTp&aks-Y2R|h@%1_D7ou{3>RaV$gA;`q)T{lnAZkiuMKSvm$9 zSq`|=Jx%5zrZx16G}4!rR3pn1+?BP&}WHPX)toknVj-=Idw zbF%iN#w8)pMWPN18Xz7vm(i*`()5+HBWa9D`Y?qgg|S&TCcspa%IIHxSYB@~@VH(O zKM?Z62d06?eWDcjXJT_oD-uRZtAi!4v}UP?y&?_Rd!;o3snVKo0Hsap-RkZ@Zjfe} zvPx@#P4ipLYl54WR$)ObuT^Q{S}N_?vO^Tyw6q!v14;D<&Ub@WrxyoXaYJzvJGdNS ztH;g716YtMtCc7m7{S14bF-{rPAaT{L?)~@f|BClnGIwaz3(XL6ln(vgtMnZIno=~frwN2X0D^J72}m4mY@$b2;VF( zffgAIdqvj4-b-LzAeF#+I0%Ey_;qD{n6eURhs{NeI2FlCqJxgHBsQRF8T?{Bq%j$h zY)JP&BpU&j24tT-|d91Z6~5IDq+wcPe9(;oK$SP5SiGz2!g0I#tJF6U9HfG ztxP;cgZ>zcam|=(^*B-o5_?!Y!lFc2NhqL#kNOAe$ZRnI|0xmZ;%2b-75=MW;3&!Gf? zd=9fh%I9z^bn-cZ_;`GNwZcP)k%KE8#SvmUWgQaRF^$|Tu8@+ykxTN znH(;=LwElbiI z=-?z4$2vHfCPtG?lEVkDxjO|;p-Z5LQ-SlV&`zV5(?g4n<)2PVs)sX(a(ZZxGtEo& za2Ap2;cS9H59e4R^>D5gIz606{74Ua$RsmK7Go_1DNM2%#DS}nVa78?Y^?d$C*qo9sAdHlBAuM?%73Ct>D{?XRUP+e#sgf?m zL8zn&hsk9yWtDU}Y!*u_lJNbVas>-w4P8kS@A*{?VVAk*_bQeL0vZR5hj$%+{R73v zhNjKYczPuQO5Wlr6KF{VbTv^cppc{$x24H5x znx_N{RFlfzQ|}y{^4+-l&_2!_a0mzSBpFHMRgz@59f>j1V>h-Kj|QVT)4r14_Rn<( zMTKa+f~0~D{>q5UJ17YsS@tML8OE+awkk#CENv@67GhA zp@F;2_> zd{eO?515xS{~(bu{}4gI{KHmA<{z;_$NZzjn_K8m=GlR%HU)$%*6u(^`PI%&GWw$3 z9JCtl4nx#wB_2;Q;zwX}3_ea6G57>5c?KKgN!Tm$6!xCMr-78gXK(-peW$lYo`osP z;B&AUgUA+Hi=Q`>XYd7@QU+0r_14fAX$=^B2{@k&yi70hrO!Km#gMKFdX*?lBhk5S zX^_{Uvbr>SlFW zCEvkEBzS3ZO(%hyM+S?t9+F1M%`?f1tlh-2?3 zG^KeQSo7`nrz{Gn`wY00d;00;^m5eI%NK@}x-W^c)HTXi=A_hpO=Q%4Ll98+EpV2) z@9YFp_r0CqsQZEV*r+a&I<1S!7E&*!|9=~gL6U2!dx*l}bx>jT{=aCB&keO99?OsT zR@ojPZs+_7=8F7`y(jw@ASL@(97N%DHu(*vEZM)q7WDu-!MM8%m`lI@V0p~>pEUK< z1Ju*yFS-QW_rh_hb^cd-uTa}BeV@Hd1*?Xkv;0hR_M5%iFk^7fFPY+K&asVL=SK>Uq@!=3zmsqTV!B8Qx?P zC-brpTdz@Fqh$NNWj-9^Wl@DOy>^&%tg_7!g!$>^Vn;z1Fr=|#L84GQBMba17cwW+ z&cZ~doka++nRXV%aW-}=W+yOqEN&+_J@p~(sY{sY0%ZDa2q+wdwE5iscsR2?9L-z; z-zyhA3bG{36t$J(vV<=OTNph!{u4%z8kWbr_obU4tyeIe?@vxKo&h#)b*o9T|c@K3)4WdGEpc729l!mGba^8 zEs=?#K!DA}P_#k@k~%AN#Mcu~K9C@EILhOW*8iuWq`&_<(!f8;2a-m@D3CP4l9N?Y znqjU;3-(@Cs{pC2R>eUitD9vtn6k3!4_ii7h$6|ODcC5j&hl7R18C}5Rt!R|birn_ zS_8+W*6Gi2}NkY0INHLI-JSeVLcO`@!<>SQf*QdzA{WU?v|U^7__vO>yg zuoXI4twa1MoG~=r!~Jd5c!cf$kgB3@BkS@{@}gRgFcQ`Ju;f}N^tIYyuSf^>UR4_a zsj4=_L8z(;^Tf+*94tZ!_RVdkY$8BSzU89{)}q_U+IQYs^@&`D)0;`Jz% zbnqbQLMBQ`=0Fg6#%2cx=*6`FoB;DqHGI$MMh!o z+1d_B+1eflAzS1Bk{w{mvNakut%q%XN`m!j3^ zSe^C^;?ZN*hWP0#dq{I~lL?0iI75-6{U2YJlU+&w1wt9#vr_itH}^`}i!f5k-mv7A zQYZVsUXgvV_e$9hNR?8-0hBVS_dqwP-YNUTlvT;gLX?i`5cP_Sr!ZB zAeuN_h$W`H!C*)~+<(+QJe~78nB~|^Cx_s;)H>5b+s%E*iEU{!aVWh~og8LJb#gdS z*e)RV#B)lHFejDCkwhkwqX@8>OpeBJwfXd8>?G!&$J$9wKF1MXh>IlLt8QGKhmg3E zP{)jX=QZwbIyEH&$2tmVLrsVc@nbG@euH5nxgEN|dzKo)o8>2aWjyo=-!!j=dN~31 zikyhOSHnp_s)mzspmD?QVLnGrfhntoQ(;S<1kKrSBn{ZvVz04Wccf2(K8;?nE>5SZ z4}QMn=4W?e7((piGw6xU6mlkx^S3phMXyvLXB$$5oI{k2BzS7{x#p!3IgiLBay|hz zlgI^D$VhUb6*?(gMEoY45lQM+(d9SnyCR^ZlapMuc>hcN| z6;^^)UlaDvKP_6#!zhDj2hE|&twYdEXE=_@Wkc~QG`X06r;j_A5Jqw5QdsgzYmm!e zugK-td!=0gq)NLI2Vn~}VJ5i>rmWJ&!KTkq(PoLS3vA;gg6Mbb7j#_-N{>N`9fiflUi)J^!Iv_^-0amwW2m6W3gfsg`~) zxy7VJa#L*jE>L%}T?oW zHb6eaFVI4@f0aepjQZDbTxy*+^E|Xk^m5b}<#j_!{ToC%>g(i9^HS>HA~NdVCctLY zzhi}@{#`3{)W1i35S)9g_Gmexwm!*IgDT%M?U@ZC!XuT6T-r}X`v33`R|Bw z%-75J=B3R4KxEASNPx|l|H%r;{LfbCnE!?N6gX=F0g5#jOw&kCQ1j+?CX30)Cb+87 zQAx(UAP0AY5-!w^&Pk1Oi+%M8ZQhZj*L4%uK)R+AX(lJN&ELiWOugPROy5g&uOr{4ewa&BP0@3h#3hJKh)WU#LR`uUDa56%&@PS z$I`e*_hbQ=6LUOU$d6%itS;K$N|)wjBiyq5le~_WBaC!Z14~{<1?dZWMV80j>u3ca z)zOMLh;;O%tOQe5M=QgY)e&1eZSkn1esqX+R7(?Y#YAM&HqSYSEnk5?fs%^ArPkRq zk8`i1ms3)s)EiPIH4x>L)Fh4OrAlfdG9@(=1WIbLLMmw$D|AX)l~^ZoTr~OHRu6GX zC6tVs=(=C@dlFk!5wJz=kLn z-fL+MAl1@996(Ec?a!*x22)l`Yr^Kn=CG6`7@I>cc`cU3GFqFaWNfYjV{`58IyRR) zY8(QT63YV#4Fb-m%!BF0KM~1OvY=}nTGHfcU82nRTwPJtGbdHi`b4Ipc7i}f9l+Jb z=QgmDI6k+bo#fQD5%KwqQ^xqtQ7& z@zSZ)wVQ52(nY^Sb-cfK2)<`t4|TFJ>=oGrd#{I0fm9Eh;UJ7H6P6Rc>n^K@&0*6~ zeU6zx5dH=bOYo$RBguNmMVehqa zFpz5H5FBWmvWA=QDaoNQWwmk`Y(p2tZ^lF;8Q@2+E7!anPPbSzN6^$KIukE=uGY{T zNzXt%M**YzmAPg2XnIBJIfj;0J;xGd8=3|=&YV;{#}k=&P9O-xb0ToHhUO$Yi4D!k zc9Ij-Da1cfR0t*oODkM@%TaW-C#uoR*7dd~*CyGm2~Kh@vvrZmb>VQ=fDcYn{?Tsm zm#X}m_6u_s`qd}@B?@vXet^9C8s#+DD{?yaUVUc(srt^uL8QKQN^{@su>N& z|IT4qEV*-ON-~-u)r58)%LAdE4~*y39QW;ZcMR^IHenagD-zm;w4_42h$t(xM!DFW zRA`qFnb0mJ2!wVSa5bS_ZYNP_SJ+8TXjc;d7c&|gEOZ&okiS>Nlw~wm;d|!w&?Mtv zugG}py&fh2sUEJzLDX{fm1|(i>fu`0s%JD1#C0r;1#vx1meF7`Aa1p8U}+$S8-eo~ z%}w-*>UkcVl?wo+At$0PI%)bc2B{)vc>(JNBR$Ozhb7U7O3CApm>_BL9NtMcrLifQyeP=W{F**0eL@&TrB7kW>7u??K7+j? zpJVTJ@dc3T;!7Mvy4YI2f+?$uuVM53tgJA;VQDOkZ)pnqSy^p-M~6Tg-vi^-=+8fK z+`AL%lt}i$f1p>SjUQ=Aweb^CP8*H#vw5jDejze#{7MjL<2Ng$Hh#B4r;R^|GySZ| zE9CDwMCh2I!(|fj{j5|e!zUw^zrDCNga!GNU)yWtFTzMGy@2vsDacf?S7d7Jy;i0H zQmssjgGeiz$aFAewK6?yLuWxZLf5x+4eJ_&!J!#wj3qK7O}1pV9r{4Tn2EMP6*B|p zBgQQB!r#=ZslE;YV@;WrhEx=@5#>bDB(s~Biee5T6UCeafhgv(LW*K;D|Dimhd5m| zt8(Cr*z7U7=pg9~XBQvC=WeM@*rgX30{`vEQ7{6d3*sDSy?Is7*}R03@=7+r^3t;a#TM$S&TL=dcXFJKlFl9Mg1U97yvy)xu5*ONKQ8W2_Ef=GS z%z~=+(votv6j6?|W?9<2l(S`sjI(74 z0?wARLULAPg^shn#GB&tb-)CiB}4RD856b;^bzj4l#Rn}P4=iVxVEDafEB!PGy};c z7=v>gonhg$kSx#Zd-hf!#OB#s5h%}IovZ|VMOMb%v)2zu*{j6?*qhYb?w^`L3NU5a zE5fEVsZEDRai?sX)X^4KxO$rS8$0^?aRF<;Hff+eV6YLm)Vh4(scV;6WWh8PH_;25 zXRw);l))CFOgdPP*I=8I61Xan5x5#bKwy91YU$wWb`sOU0d|rjwv~8WBbGC{!JUA=aH%^xSGFF# z7r!6J3ohVx!ar<+Jc`2J45E+{r5(I2+WQDAbL#)uTwlQ(1v4%-a z{@A*FR4bM~z7)SW9!Z%zsFw2+%%P(sN}KeweSD zn-XI4b#pVI9IJR|3B4+^Irg5_Er68OVK|6b#k>McSyo5D7O+aPNYs`z$E=PFJ&{Sb z_f^lWXbvdd8o1P*eVf)c^m3HeE2VYPNn1+kwnRBf8)cMvDW%&H8Kv731eESzg`{+} z6*@|HB;Lv>#i$|Y7d)j2y$r0zrK8qS>Q4ZVXd0a!lw2iQPh=J?JMoh|pF0y`^L*|C zl;^Wvx?smUFR=G~?h2%QmT{o5c2F{keBpE=FlG513tLk6SW6K=+Yn8)?(N2sn9|*8 z(j{H#v!+o`c?LDG2TKEH_XI9=&$ic$hs`wcyfp6aMK7{gz-vilZ$mnTLbj*d0fw6uct3y0jFEBAV<+FrvGS~^huwL)8=U}r|=h% zW9S$t;aK4Oo2-tbS05;$seyMD^0w9F9@e_LRr}PhG_+66hJ9Dn-0FCGrZxWrqEIEg zjkQHiG$&QbNkpcSlL-QqoMMHn_@`Q-Q^RRm$OT;_f>k#O-|p1}Zm89kC!kY}d`pvN zeXLzFlA^=&u6N@m2RQGbi}&3!a2>)=<`0|)T6CP9oP|X;2qZncs&R5>K|e^t<5xEN zd4hlXRpoU47(OK zWki`KrvdL@y4<`}c~=mb@~$KZly{XCQhDR7&?#>`@e%lZP?(EDolM^-5Na(mIi6?+ z%AA5IG@SZyk5JLzLs6<|AVY>CNY8@M1b!jU{?&xoJp0!G<=HRDwXj#@I_y3B*8?g0 zH{c**{{y)ZrY!q6!KPD23{xrlsSf03md7Hvg{Fn^1*21vnQ@bY+)AfF3%3F1wQxJV zuz4-qK})KIJBe~yXpy_jOSNz}k!j%`f!#C_)Ib=0IBpk^55_dYFs5K? zGpXOCpIK$v@AE22B8xq>gAR-{$rH|fh?~+6i5fMPlRe*>$u7Q+-YSWI@v%wUhTB2- zXL%*vM~KZU>3*QRl8W*G>=k(sd#|L2fK*8j;~-Mf)$#~TStUIRTeVjz=$cSG6^g~N zkRGS0>MIq}up&>;C(zN8!1-4yJVh^TUPn*UlIrLgqO6Ya5YT7MOLg=dk?H7pfB0VlkM3mgOIz{hM@Pg@KO%$6Xgds!==#hy`(lmxIW}rBd(IMzD5(j6o zR1`QJ|KcM|953f}|}EAkrlUL2EvR2;A4K*eEiQrKJGfGI1E zH(^Vjn1C>Y|AlW-Z~}RYF0nw~rm4p_DR_~*L%%>I?*ixFr0^cSuz8WZPfIG24~TLi zDaeQBr6T!=$VBonK_HS(tdI@Qr&j3n@)_~XiCz-&*u8Oy&(7d_vTss=Y+S125{)Fl zxTn8@=p@ig(S?&7nF7P~k;=zfrzD#AyCGFhyt}R-pY!W`{d_@)&Fkk&puB$SeUe!j&)q@U&FJD9Ti`5v~knQE(9yD$E~HecA{#Pb8oWAXe*lQmQ6DVb*K zCprbn`5Cx0@!rfE6@H-?Hm{stX-Spy8&OUXVm8l8w2YRhc14~{j4Kgk4 z6`2luua)V6R4X&!AkxadG9ye`t;__QvXLaTc+L-Q#YWG0X48hmxy(z2F*lJ3V;+J)7`=gW!kCv% zyHOa-VVPiRpBH(FPCusW`4BXcL>9q9XQ4IcY_<_GJCPhu!IKx<$ub{5-LpGCVZ`nN zu;kfolm%hOYs0bk>@EzX>@I?Xklk_j%Azo3*tiW&SHL)ULq=}Va z$!nrXR)!sK3di1Sq83OsQNTf{iScVm5vHsr>R{7uX~J;wZZ0$qi6rV-iY<^t15LaW zDLf%kL?b-{MKl4I2E3M;-e@*2r-&BfR1vEX+vBy{wI<0|RRhu&4ezIZsLmgBxVb1@Vj|E z3?z*7&<0CR4^0iSChU0SH}+l+YXhksN;ru0u&4}zDXWLUuvOJV7n%ZDho#s8J*-Pp zk9y$h|Mln*=wW@}(!^&oi9)-1K@ZJn4LXQZJ#0Xf)dTLw+t9pJ4;vAg9)=JEdf3h z4y0-rfrC&D>El|YPG%XTF2zOExCAOhgpo&i5e7mtN zy{H042^nQb2jjOR%4woWwl^=;#12HJiO~dsCUyjld!{?l<6?Xunv5yx5Ks?k>If<1 zdZuGpsWm`M(e*M(WBQO4w|;@W416LxEwVGetCz+ugi)aAf+Z&n+!Zkf_KNI^y_ZHA zNTne-h@^3ujD;yHjoo0=!SS>m2&Iv}1ql+_ou1eNiR?ksQawwgGd&;iV^2BB!a)%Y&3~*Tz0o@PIM)SEM(TBX0HVUf>3KsxzLP%BK}B{D|4c8kg9#&% z9Rf>UWJNg?cD$e*doQxXfmCEi;2;v&d2%F7S&#?T49fG*eqWn4(pqAL%myl|ND+JZ%z=8 zH*c4|Od!4n-%dNV?ituBC^qK=F<9SUgbR6P@tb>YG@EEs(~|kvG;xdhc6U6v*LL;Cz4PHhNM2cqduY6Q4;5_<0#68nmGT zvy-ED9CKv6z*PacM|gotTA9JPfp5)PZ2pC`yC=P&)j$+LQ}s_SCin6idO_bu7zK>` zVacljQ!KD!iUoVGhKGPu4G-fWQo~?*1g5MS9)-;s1*Zl-XP1tJK1MHWfiNDAPE=u3 zIq?a41sZu0xHRCMOdsbddPN#}nwC@}&k*G_(kRcGmuloWBGbt81c63gutI9&MJsd~ zd5QS@Dbz>}G66RnYW*9REx6?XE>*#Z8IuNfk1=NEMowLN!Ngee3jN0h+W%Hj|C5iF zm-*Lw{k=jM>F-ro^7^Zj*I=*6B<#KZUI$YBy@7*Jf8*YgH(|=^?=9Gp_4&A$jLVVQ z+lu;tIXsIHtC~9GZ8~5Jr1uU@^TO1|bJ_|;4PSV`9*%LwyDSbg_a1O*V)iB2@6#*N z+y}Izn){F_r@1Ei$h=f@9}}78J|PG+_o)?9bDvqE)70tJnHd=b=WI1J@#H6 zGXSYPX2d}hR(6+}V9Ls4X4sO|u4u!$I=tzrSYLdK&12r~TwAg|X{4Vm-Nq&;mX<*%*sL4uyAZ9M^%^T6F z$Z0a%&d(3`+%7;Eal0Tad2Wlc5bPCM7<-P5k$YvrYCeoO26()5a$U51vF z*=31x%r?t%=B3Qm5E-+52?A!92hK6O0-d%uW|>rg4R%?QtRYD+bX6kT3SDZL(96*@ za>@~OB~Xb=B81jq5Uf0F_)qJSbkC6#RR^?}6F#G-)HW?;Oe{+E@` zoOr+tDW@&;ia1?`mXy;~iE^B_$ZF=Loc1R&PFE)gI2{0-MQT{NwNgTu!eTM04OH4?7Sp>N6i|jfTVAu5m~!F?h3L7zlvvj zAYsII8!UOY>t#*YF^z(~XM1fRWxIrfi0zwY5KLLN2g8=E_KesDo17nChov#&>(Zq6 zZ~8inaS6I7`QdtW2zXx~xHR#P)@7eqCalYeVx=-ZvsL z-iHtbyl)Jg<9!o4otE+r9(64-iBEb2T(s#_pj8&t88;Zfi@}`a@y>-J;;+t>7c8UG zwwx`z-q?Wdh|N6=Q=$f(l4-E1 zpsnc;s9_smsNti>e|+z!<5D$r(koKKwzQ;b7)6v*LqWDPFIB_#M5cxv2m&>X2F|Hr zM>^@^8w=B*7HM^dwXZcwJoZ|vli?{HcTU=I3^h?Kpk2F`&F zJi-OWsgk}XQp4cNh$Ts2!G(Fn^un?ezqr@Q&V-Rpc7Y|YlSb)+9e285?{%^(km{t2 zgGeV&i@=oC$ynHuW?W+iT6;V;scort*^LFUKz66et-c8x2)ggVvOpVq0+%NKa^Ecd z*Li>WR2zHIE7Hc^w4~bDhbX5Fyd>PbR2%ydnKmi}fj0IB&S~QSI#n!EFgR#D$fj); zEZpQ;5+xc(EJaS!{$AWg&@ke1{h)QwrQ|Nmq!*9_`87P#2N6b09}G*L=_WY@cHH2C zy=VF`AZ7Y+97Ig-DM!GRW%@|is$}g)(H!%9G)+EhXVUI^|HrT-VEI^JWbNDT@Yn-4 zJ)5$89K8aT+vIp!QkG93%CTH0Cz_YCd=indd@@16@+rVMmQSVA{~&AE`wEhzkD*0t zZmvgm8b91~dpcpn?HRD-xows+VXw$p*n4iz22yU%!9m3BzH%;1S#HmREuXcYPisu= z1vFXK&g9ORy^sX~vljs)YybYTp`Y&ePRi`X^op3hgqD=qONnyK*2`t)rOaMVWXxVc z5HNcsaE{rl==9&p+B4aD(z<8*VB`3Kp3U)u5t|cW$+OuaSHoVBYq0lhUJImbUWbE- z%@gE$n6hl%09*B}{YI9=jNU|3nzd`XZk*oCqJYy|fRVK~oN>}Y=PjLbdMmvmPH&?n z<@9!<9H$L(hj}TdcM=(=cM$}f-VL1N^d36>kFxg3r05!?mZe%e?E8T$+Ta$W~75=J_B36{JL@NN&-EAk5VUI(uNsSaMlL8OBRWfDwT9lQ?P zY2jU3QZ2kkl+!|! zyl-Brg%5~K3m*~$TKEV!r-hH{^nW2^M{3SN01oy0jGd#dJ!I^k@QZt$d`cMU@?7dFD1X7)Rg@Z^Z&&t;@Wp(lmY*jP%Z&?rv=8i8DnQ32V^Fe1VqmaTpBR`>(}=lmAY!<$ECf@Q;e}yK({atk-5~lREQkqSH1rKq_r+Kc z@Vhv0Y2x*nWdMEX74f?SEh)cC5@q=<*2+@mrTi{UWc)5e5b(P!aE{;Q=(MBp%j&0H z? zkVnJl1s27u79;N{eXnCtz-v8lX~1pL{n1^kjY#Wl1HB?%8)->-Z6eC?T99V*QeIn# zjMr5N0$x`I&hfe$oenl$|KXqrS9s}=dG%qHY>*D!t3~Pl{4SpF)d?fM2f&i&yHQ$U zugDtMd%gz(Dc@~4i1=PX)`Tg`_gb(`b`X7SI>fA(XzDqIAH?E-`N6=YiI-+NWfvSoN~x1B zn6gS41DnsyeLKO?^j%pNi=-U6$Em!qEYQbT;L?C>2Cy5wB7N*mORA4Oh;sU0@t#Kp*=6=k&2Joeng8gjGCFwl%U>tMxsPzfZIRN&G9EiQ=`ye3Y`(PY^?@7JgGsQ0WOAdi4%lDzMX|89Bs<0eJ7fc++f|%^X zLtpr%!?Z`RAYl4P;8OR@ncEzXqF2In0XI1#mp_`8l<8xL;+5ic_$!Y!CnfqgA|v{E zf`I4~fTNY-Cz>}`il1cOZaLn`#5+&{oGvA4RdPW0gbY5WDU1TB{n?X{0tU`i$8S%; zC+3N4lv81^$Z6PnB2NcWBG15q5?OMu82L@kgegnpS+J$uGCtCl5Ajnb&t^%?j=dvWA@;u;D_e@9kJ#@q)Mg8HqChmNCMO0ouOG@R1L|G~ujWGqt@h@a&E#sBvb0_UTh;n`Eem2+uSPoMIfd2B^-d-tADVEG|nO~!<41= z71-?7yaHa6+t!YS?8W}=cu$HxS*9qj(j%t$HJa?PRmn}B?#+|!Me^L%R+LHf38;P@ zxYRoTo`aV=Y4w9rs^6ejLUpmwXCR(Q{w6Ic)o&4HGpM>odE2~{>vxEZ>vstPuHUml zE^B|^3LUi{5O=d3jc&IAFC4b!L+X;$8-v;40)1=TS;?lh$WV&Ji#JT*C33nh0jprp z5M$4To&RvrBmTi{n26B+`Zem0eUJX>ofyhL5a-6RxH*|vaf9gicG?%vZ&{QN`Iq=i z>?1<{u<3KqKZYf*q9*wS_KJLpy;sp^K&qn8aS%1&@5vW1WmWVgY}!ZDo3x6E-MB9z z6wz07h(+`@P3*DroUURQpAH(x?i&^d`uP?(e`Du&^m6*a;|8Fg?`cW(^8-;PR5Ufo zkLINE`H9Hn^D{vppI?BZQ1PpIGgSO$-cC=y6W>%#Pq-Bq^dY9}%{3I0PLSYU2}+66 zBL7B9TxUyZ*hBOo$u7NeC!}k4|4!+Zx&%U*0l3sURc5euMtV7+ z@V>q^WhPouq0CH_4K(!)GK+bsQf4JGrOZYUC}nnFD1{cMhB@eQS2Z;-x_}{-1BLMa zGq@xQ>SS<9R_kfh!k0pTf{0{J{uMs5%taVQmbqcc%ce!UMaE7BKxFPi0nR5UB#K$|0b*!T>xB1~D)tOT1+YS0|DXWOHd=@DzD zUzL;kw}nkoEjkxX5Z($vwCN~WGD(OfzuQ zB(<0~o1|6D+XbOji65w@CTo&3lhF}IWGw%ydws0`PVm_Bt5pnJZvS8DSp0cku=c-jV|WBunz?MKwSZaMPwL(FgpbgID$S zdjrC#es2g%o~U}+2=nplo4+u^Eto*HC(*ekLF_MWHFK+4mOI8dH!0H{NDf+@?>&akBw zs?$}dU1*FeR2NM~6&K)vsxh<$RP72}YR$gCqfD=eDxoE%YAjKfss`E3oRq5FiHxc} z2m-411g=KaUUm|x+S^WYRP96DUWKY8t4jKZ-~g2l93+*_bCp!OYGmz;>t^59fLl=a zgSjFV>^)ce11VPr;2;Wq3(J8pWw|;CHr68EtEO{hdQ;26w8mr|5;?+4-{Ns7tpQz! z0hd}A+Vjen)<0m5v}PSnFGh_*ZF8TsG)?E|$DucG@@LCoGL|zA6>Yg?8;MVKu<%n$P!}+TlXi15@ktj!G zz1(D8O61K%M&vC70g<;_A&I=r3LTNR6YH9#bm}QtdqB=mnUovulLy?yvy2)a-7J?* zeCLVO(Nna8!DP4G!O!&E-bom7dlxKuZtLZ4*eh}m_MY2`K+0`54#IFc{w%o{rYyJj z!B&;q`&kfk`#|c;g$=BpZuaXz76jZr1YBy(J~HoNdIj9BDUZ;Sa{DMzj@t%#%)FG_ z$BB&FCkO&=pR_`9`;-+rZl5OpkGOS{9X)Z&%9`Q!8Gfec_F2M++vi}(bK4-#!;Ud? z>^-+H0x7pI;UMI8!c_7yOj&MUfz7y$>f);`h`D_&_04cQi3I_-uLI-VMI+xn`l}np zrvdm4dZpa5$M7aCDYtJC<+#Q8xOpkJ?+_Wc?-B&uzGsCxK5m7M+Yg9SZgqw+Sz(e; zrg2xr6evSih zJE?b@>slNoU%-^*_Dk6G2j!~-Src7I{)+aP+OJbrQtNB&H?#-LehXY0@M~s8!FTjZ zm~E=>LuS9HC1v&pqA(j8S*XX~_@g-~u|E+Ru|E?8#Qp*t4HNuo-a2w^-j3(svao?ccY)v}1PCr0xVQ}N?(Po3-4op1-QC^cJLgv2>Z#QUmQ4u~yQR45brL3|EGLb8Q41H`Cu2ecb8_g~ z*m)4a+t6~rlcjnoQewnO|bo_?S`zo=BV8Wvq7ddCxtl;mBE~r zB7!*`bWh!o>Fpvioxv{hOlPFt%9!dZksfMYH!qPtemRaxSNBzNji;XK+?nu+WuA>P zGwc`&#COgYH@u(=hTJ#L~lRFgNQtAUn) zwUH{TO6==7?-3Do2YT4E4MGSa)~mPc0I9+ubYD= zZWPJ-gRZsr&i0~SkzN`u4H&NAURf*A62fI=sytj! zN3UXD8ZN6+S-7l55rs>qB{E#PEYXL{>eS&qedNzf%b?Il=KVk}I=^6vBg%I-7U&n3 za8C{0e?v!@G$u;7@M9sNLxp=pbvz@CV~hv{_RM4rexZW?nv@Csp|F(cS7a^N>#{bE z1^soPl>WLnN$9UDEL(H>>%*pubgtqzTOiH#Uzt9}RT8VVKs4_i*T1 z`}}u~{AStS!xRq3eK(?)=gpU`cgP4@0`HBfO1v9n6Z2Bun^GC?krWZ{%`B0;H@8I3 zdkgCI@!j#~_HO}o1H56>u$$u?0eBd%*D($q492{m8{$B|1v6&(z?Q@ITk<0cf?H80 z1hf9VI2HuAgHnRq<0K*2B;7FO1V_Q<*l=x;<=q`z0bH1ch)+GL$io*6jqOfSc=GEnpFF0=%WyHb@nVp+C%DaYNZjN={@5yw3( zksSB3M9*<=>aC3!CDxu9`7j@2XP{V?xCBaA?e)9lPWFqgUA}a6N)v z5!Vhml9s^rD5?_ICOO)?l#6_@ZiQ7C}PUJ@xTu-7*xSkA4nQNn*0()If#j)Ue8kBN99VZFb7C8f^oa>pe z#r3V`SY+IaoW+!s?b$R%^=-oU941D5&xJ1M+we~nokMV072BU>Rj}=B=&J1a7XOa^<@`bg z{VOOF`lDef({GY1VXwjYd8Vv*JfG4 zk8Hjqufvp+egigrikjf?{8O~n_Dv?GwBHKe>Kw1-eVa)U?02ARWB$x4>v!qJr7l&x z6>`&lBjr6B0_*pw@LxQFf*<;UIVtK7sSNc;6cOr=p?ew;`@}BdK+30fkq?2-s1IQX zfKLI+{Vi`mv3Oofqx5|*2vmn3a?^v%bLDpVv78Gnz-*3 z0;JOEmi1TrVo`QQeqd6B{72|o`;>ETG1Gy^b%gx=lX+Dd8c~1C&(wkVFI4z1CXQeF zt2rt0->8iE?-UX7KcIUe{-<3;;(ysip7;bPP}V!p5!YQNOSO2Yy}6!0aMocJw~D&^ zZ&qozhQL4GLQ8{+Oo%U2#@!?n!Cse%aV&680;RYo#fjn`?DO~RG8s%c?#W@((!}+T zu_$T^T2jVS(qt=bx&Rt*r=l%_J2iBzeWA+qYg~E8Ujgnk<^{M-fIBU9fIA%({_+Xj zCYj!x6z&XE26sk^2<}YKZtX@I7jkpDp^AM58krigYl5z|_qt@xjh|mPDdP`A6^p=M1iQXW)i1e1Ri#)x-)aO|aOCK|L z-Bs3%s3rSEEIIc{ROja^?G`u2B1eT%7n)VQ=(8l4e(`>YkmGc6`%vOti|fu-9cd9E+IffYO*4f)kC2A--HaQK16c!9?s{|emyyj0i*-?J15Z3AUn*z2+#jzu7B52b<7jS~$7Yg}(AqhQJd zVF%dky{UMdt*9qTk0}}qRkGgsL_}GW~2gXOb z6qCzti|WGH6K+DI!pK3bfqoa?CV*1k9$p-;2X4{Wz}M;n20gUuHe} zBk1MBtU*^S9Z5?Fv!kf;FssPX=A>bE43&l1u@q659S5C++41JhFgwA#eVCm{eOpgq z7PTwbY5^?&*aM&;p17%o>)!K-K6|;j;R6KEd!tlJqBt5Lu$N-k^WKOEz96OA-=gH7 z?ky+bCzL~_DksBUms4;oLgiE_4VBYyqM>4&s(d7;!<2{08L*|Hg6>2bDreFq4VANK zibDm@Wyhg%Hq)a}IS0BncA}%d-QxC@289~yTzYY=u+q?^i`$m&Ty>dBwV}255IK)& zA!yF0%7dmV7nqX<&4pAJG#61sL31&55;T{XH-qL<^Y%e=8TBpy95fix`~Qm>lqi^= z^@Yg444KRE1Ii)OAXmU%m(e&DA#){^hRjts(U7s{CMTAwVah}18rUpk_;yU|%5XeL zAL_5jwRA~C<~o|+|>>8(cBO zhY(oaAfkL`0rMQbemP*8Ywa`ly5y#*r*~hGt9jm} zmoHGXx=_LKB7J-jSiL# zZmSuTC@3yalyHd_%Z9k!i^egM{KHEqfRZ+euIThnIiUfJi}`*nv_#?-Y%;XXO&vbQ zw<*JJmQP^E8hRWH*q=cu?9Xwcux+O4G5G?f9QK#689S{wI(0V~Mg0}iQq*75#OY(d z5llC>XK>55zF|rP`djGos?P7|q@r zfd19I3Hmql_Mm^K-n%Ez#nbR03Q&aI!MbD+;=r|3#{rqWkg(4STMVIB;-8_NnV!lPu%O&3ryo3TlUm0Ft3LVw?v zpstDM=Yl!jzy@K@okSGb?rBC~_P5D8gueCY?ujB_7R7fe!>&jl*sw={6=^ajMchnfh+8Nk#I4W?;x_Xp;w8=7BVLMn?Y~6~gtn+Q zuUK8>Jcw1VdtjKgtJTwRmZfp~GOel%f*pI2;8@Ta45hT%aiX-Wh20>_!j#im4z@sx z<24oOFjGNm2u%?!eJnRI5QJs5;G?sJhJCqgtK1XH+($%iQ*F9`ueLRgyJu`!cEqSrc~bMuKC3YAq;5 zwKh%^l{I(Xk#%6oQLPJGK!s8is930G<>&e|MX21|?gq3*s5XQy?{YWHynqVxgEE{t zpxTHkN7W!B%t=vgOl7Dxp@>jz3Z0-DY2HM&nR$Cun^X6U$_mc$M%8Jp+TxM?EpYoX zsz%uoc5Fz3V}WXGC`GjmPLf>JOSXk6N3|Vn0TpY!ifnJDB3E_O6i^Wq%fL~zMyPgx zM(;H1b?j(fm4+7VGbchFQ0+vOqiU3$%}G)1LS?9SrHD}N2A!bV-Moov5A*h@_N49^ zRg|m#rU2B@PG%%62J-@NFWkP2s!8^S9s7~sSfJVmN>S~LlLXb`vL8%2s{LUr+9s}L zQf1ycfC*{a*?Q_?aV{yw{5XNB8k9EWdJ2Hq^k!;alda4dkI2&KSJ!ifU6 z8MQ6tWSDZ`r@%H&h^I0!g?<`MVKm8&5}<#fn@Vz7)9Fl)xSs)yy@o$K>6P}~PYvA9 zq!*{^E7hik-kbE>R7aB<5CeRQ_bjG{7&x0MkAY@6$DA|<&ZV*#IFBNVf%Bo07`VW^ z83PxZw~v8~s87@w2x}G7M#=c#^Bu&VS~u>wA;0Z{e-c6aq#90|<^4beA}VSKF%&Mw zmn(-tt6TzmT`tA32!+d_G!!n!iH3spFK>}6V9G;bG;F$W6q~ACugHhJuB0&yfvafZ zYAMb4Zp`3n+9K@NK-czz5C0diUrR4`E89?^4gu-wXb7aQr^-pU$_?hEq;I4$(l=2= zq;G~!NZ(@KBz>!Sd(yX2uZpiz-sh%UTc60GZlbgYjINAUw+^}~D7H^;$6d>a+T;${ z>vAWK1){s46w%!{QAD;X=3*HGQ;ujXY>KF(8!JmIc*%98OYWgDMRYGs{7dpNO!ofW zW7zl698tO-x;FTO=ld;v&gWwTr3dI0Q+kkwKoKl-SY)(q)5h|nfC`Cl+G3bQS zZ+_;CcBd*z1=UnV8ahi6)kg z$Y35{XC++lGE*anuRzxZ|Byj^m0mH#*Jub3U#H3-;@7@mP73i&DuehIMFjC}=mg?B z=1qw2nzx7e9`(6tR#S#<4hKuVJww8+a4NT@hgIM<`)`%*xQAr*hFX21^T_1{<9S1Q zA78W_2Nn4M_PTtCV-W`*L1`R(jFY6U-AO)yDUX9sVbj(Q3)oSF{~0YQ{Lg6$r7L$( ztjHI%MaaK|E^nXk6}@8QU(*nfe?yfcugJINq{zRcGUVS=M96=DPLTg--bDVBd3)qP zQ|l^JWvfC$D}!jzGvpv!A2)c)%;~w26^m9ly9cfmlM>@dDtO`24r1i^Ys+w%4xdqG z)*#cvUY8kgESSv*rOal+i88a*7GKEBFy+iCu;Bn(Em|ke%7ipy%|;Uj z;IY+Wh2<(IXSmxL$3&CEOg>fu6E&`<-7sZKkw85Mgqz_Cv$G)(|9J_%>IWES8l;h$w zxfjKIj%2AQ`!Oj(*&n(#=J#w@z9s0Dpd3I;Ksk^qN7*P9b5fL5Dnr>o5ut2^PEa

p+FCA99TzvR1HKbdDENeRmOvCV;cqj#hkwGTJ6-gS{?G z;#i@LqwNJRlkzSS_ z5$WZi%X5?+^omQ{Av6Tg%Twi`n`8xZQqVOj1HB?e1bQXt1oX=0P0*{Dw+Fo{^_cMi z%?6OJmMaDtQ2e4;wly!V=H!A+h6dHY42FMke-sE=ysU;FQjV8q>4d#5T{sr;vO1K; z%NjUI;^iq>6Q(>~hQemmFd9jyhtW=2i^eoY)~1Q@uybx_7pflj;;MCMjsjy{=-QZk znLj7KePD3uZilQ-LkNrwsPe#QmJQ8G17jGK1;%iSC@?mHP6A_uc{4CJHg6von^1>p zIFtaP&IMhrS3v2_XIv2mXxwMh{cf2Q78(FSNl}>5PZurzP z%U0PP_PXqWW5IGyC}p`9PL!n$F8Ucsf;1 zu}#h}C#85Ml~Fv4BBFRUbVBhQ^Crb}&D&Evk2)?(t#!-2omez#j_k^+w4yZSngMoG z0dEwYdSZ4yKBL5}T9FH2ugirv7R)Y!Qf3$9M48#@kS4hVrkvTOu(|#}SfPu8!ci`x zEd_RYaKiwBce{eNh}UT7@(lEq^dcV2m9@!Lrn)!Tesj5+DuadJb&WYGtZS(Z)^!vS ztm~l@ST~qAVclrn9@b6N8?pG*o{1HWEH7E1fdNZWREYeC6`JnQl@9yZ0^Yyzmz(i9 zWnfjg1@^k!iemxnHYf#lJ5CgsE#IC}?tm!=b|-8Fur3957i}r9yMr5nIkz#iMZm^F zmtSpj54{4gdrc{@`>1kYRk`1s6xah)2JAtK2-ri=39yIFo4_71Zx8HIY8{MqgF3zk z$p?mQyr!yI3!nnMza8G`!CfB1r<8Fu$m6ir$W^b4`nZ0S=p4nT}tN&9b8}HmzILh0&bD2}4yaRh(-o>%t z^d6LQdLJjs$%b}*k`G|YIeiFQ!3p`LA|IKl;PkPXl#_aWLQ}-)Q|Q{*ds|k(Jp%O?bU8y!fT140M8NMJJ0UG0H%&y9!)ld@%}HTRLS?WfrHEin2A#l~ z+`I{E3iI}`rlj7+!%7Sb&@sNE zWLgSyCYp4_6}yfF-$6rVW+p^5XMwJb?URjv&q^FUzlY;DJ z*KfZ24l*Z21adCu1mxW2O_1}Lw+A^d^}z(mo>tRvZYZTw4Q!+jK)Vl_(Dk)>Ojdok zWgFnrci4tQF(_lZ|Dwh#koOtSTMWW>#Mt9^K75@rW9)7XdtK(ov0%IalrmlrCkf-# zWFeSx#tXyNGvh8Uu8mzZd%N&r^fJbaQv}BSsB*@Q(%+nv z@e+1D84sX{7!QO_7+1`jjH~AD88=Y>_l)CV^Z$@Bv%NIp>y#NcNfYdKX~wZ&+ybSH zTXCX{Z6l%?r46Q>@shA9M~=a`w?i#urlP)F+Dvf4^cBWlWDrde)McP+?F+BC+n1e- zyr3MpQfn~1CPHq(8?f7X9+)mml{0OU<;+Q$c2HRXJcJ@*x;%8g3DF960kPKX0*`e? z>Y8O;KRCkT)Ma38d!h{H-&{L(Y>T4I4W+U5!D(X&EGyxLWiZXMGVFC(1;+x+s!$4M zHJm6I+uf$Gbi$N_>4GhQS+KJ&M#WY)QvqfTGXYGW4!Cv6nlwdVhC`Xh*xq6B26I<3{dq4H1qt`kdoBZD@LEb}3pkNyDY&@4M*eq3dNF*4hc5 zQ^MDXMYpimWfvR^e7iy^zTI%5_?GbxLQN~X!<6IO12+BBY1@;fn7%b{PnuJHd&Qo} zSY3V9b8ng>fOY6v`y$yQ#eL{C4FGPx$#$cJlHNG#Y5~OPw`M{>xOf91l)?|q>VB*xEJPOu#sh9Gs{BMWRW|H zJ{JL2ATr%wg>8@z!>5!PHp}6#*X0Nt3x-ERDZ`_1l4QXT{UCza^AbB2?k~|+LO40_IkCh8x%1K@bn`>NxB>H#k;gKXO zx{H{W!n~L!E+Fn0(9vfAyl`|E4O9M?Fg;>=DRg-q&1L4*(5y%=rw&N3pvsZ9$Y^s? zq*qcI(yJ&Uq*p`tlq;{Xi%9lbyU3Hhj`{?v{rGIC&z8tjE_eC@)vfjo1DhO|R-lT@ zXV+HiRJ5UIqBk1A+=IeT~nbF(JEimN)a4T%CI7a2Db*WaTx6z!kzdd*Yb>Vpj%@OiDq01{8@1mDx z<4P-+jLY3L1m=(V%*^JA7p zu44u2S8}IyCOe15!M8u7%-Z8+@&G@z*pBl-N-kFkZ+m$NmQrqRYLJIvugfDi7BC-$ zQkak7M02xE|BaT%Vaj1X0UO(dm{xZqJJ>$2V<<|ME_sqkDbS~A((UKzJ%FvFVME;u z+`8mxCPrAFfi6!1K1(nDJiMEWkAQT@bF>7m&r`)*M^u{dTVF6I<@zF(aeav*;`%al zPkH$jyNG%DRlCTee~tPSd@V3)(8oN8of)@FpBYCr>asLPxltA)>+-Fq8XP_ZF!aNo zJxgSc;%HA~T(@FocHzW(VN^lm#sW>VIm&t#6uv}%9pA4U15NS<>~(n)$07#ag3=gx z8z&kA)|H%I-hnBPfp=lk7+3-2IpzboSf|pr&kB84koTCB#=!eDu|S8nM!*2p(F_H5 z$p=i2!r((_G`23j{rWj(Jzt@*fB6x;n7bQWTQmwjrXfVZCscV9G|8vtq*3r0l|{kl z6j2m>0i9IRUz#`Ze`Vet|JT&V6aNZIRcl*D0C3w3A=hb)>k%IOc3g?36pq>nEp8K` zMg@&gC)Xt7olvr)z;LTx5o#*Atwek;`37IDOukvZg}pA{;aHIW9!kmofRluLFZmIs zocvF)X_J=|R7PHYrZt8B3r+ZM0yw77QaQr;6_r;@;XLW+p@M9>NEiOrk5Coyl&ds1p0DB{~Cs?6nzAuzm{`KReM}6XzK1bEtZraF z)bb-yKr{O@w<7=BvnVVWptw>*%Wjs|ZbjY-$lqFYw`eth1wh2Bagy8~LBaw)cXk1! z9Xo!q78|~#2wR$IgQM6p&l(%m;W7|Er5rnLQh~iLRUC`hX@Jt$X~c=fjy*HFqcp*k z$4)bBT1%kOSJVxu(+NNxY6;XBY#oe$Ex|%3y}`e)C6>q6 zDKTyyC@a8Tml}=*;}xNl@k%&R#@4GpR91#5XS@n*Hbh!5UX|9A@oF@YF=pU1#+^)x z7G3RZ`jozGnI=E*IkgY2lB{1 zcz}CkMK?-k!M_O|Q_P#vG{62?vTLg&Q7YnLVXoA!^6HUvjj(S9U4HH3=JetuW2Ld7 z!L5|3Xt#a~CItLjQswv?Wh--1{998Q{%t5C{M$k&__s4};@{r9J^pU$8{F??2JKv2 zoBjmfhb6o<{xD(#rM&^`38HYuYJ68spi!^wcvcft7W=3mwD{M~Kv)0qzBTvC)q2O> z5SujH5>j?9@D=n^^%I86D17;H@HEK|u-9cr9E;!)C=H&SaH7Fu0|O_^&M@V{vkPqY zP-+El0}c;l=wQ#TOid$ZH=1;N2wf%+J@$c@aO}>+C`R^xE)V(aX6j5mG2c3k*{^rflIKaGpXdFoGo?%58Ao~Jh#rM=hko2%c z1nLHILVp7P@X3D7`C1Hc2vAE5|AeG}_|6q98eE{whe7JQV2fB`zD ze*u#t))zvT-@S2>dHFH@i>U+aOQ>?zZE~qODeKFqjP>Od5$h|UdtyD>E+XqI?IO?m zD(VCNcdXGMAna@`A7k~1lQMRC>)0Zj6mm7bONn<2mS@9Wmuqn>cwYylysyWJ^0qep z_HqMEIqw@`vt}Hh;Tc+NG8{+bIL-JE$_!_?35>lajuR%1GZ$5s@AP-4p4tb`eS6V;6bS_fqdq(wr-G#k{7BG>WNH^Cth}{k*8qv1`*hX%xGqOabCioKet{5aA}K9(WGRsphkhfdhI;9#6+B)L+;HwJw`VXr9Wf9^CK759h_eoCj5sHp%VQcnHJIi3A(&H(a-cs=>I}Xp#LjXPQOupGbg40JC)J@ zgCe5;Cv;Er|FVlne**Mu^de7xLh3V(KKXLa?{8zk_Y2%`;@`$VxZU_hg1a^PjcFuI zgfCfI1c0^jFxO=g9E(7h6iNeOGMpd~uKHnetjXH#*HbS)Q6`5e4}>XTbG&$5w;)d@g;bVwT0zM z&EqVrQ7(k{&qL^QiFswfqJtF$TMllnNh=zh_gsLu8lWn(@td(>waI3iZN8bzPRYy? z)(*@8OQ|quZINEE*JVx|ix`*-N@HMdoM>TSk5=`Sd0@(8U|!flk3Y05bSt@v%*Vtu z271$kOyxI{!!Cku;el8+!9O{)K1=VX!b&90piB zfZw`^IcXRyN@Zcthaw7tzR<}C!D8mkLSS+8_HogVT009`@@gr}(%O$5IhY=x%l3Y= zNA4Ero)_k5VKj!Qf#p_t@wrvDO6gmK5UUWVFn36SW=K`~w_)PGi|5N@;bW(_$z~k z=1}HgU)C0He;662hd^&b@x zXmtKV8PV-0tx)p(>q6pGe9>|!w8&|&*X48^i%>WNN<-mHoMBl}n=WaY$pbcw>@eCXQXce2sV3+NSx!-X`2aJYyn4~G`H z*qk&RE}^n;xRfFahs&UoaJbyO84g#Nw-1NW)Yp{D3D5R_T~K5J@n4h_SK^D8!=qKM zg1s(R<5+~pHBcHJ*WyIO!}gkdM6QD=50C3%>$&BMlESJ7xq)73l-x)Ymq4(@h+3`| zF=1E3irhrkC|YiYu8qlG?RyKo;%K>*h7c{cQRUInDz}@HM#~*k7A<#DMA32=bP_Fh zn>V9njCuQL8B6^mqa_~INDB|!l+brcdTd5w9z(Y>(Bgv)Lv2B5^VQcHJ+&jzGPTyL zYk%@GR6kMC^7mKcLrVeZ+PwexWwox7d+__qLDwet!d{pAa4drEekcvP2XLZ6XUi?u zkOyJPgYF^NlqGj7M1R>gSRST1jj~5*vMtqgQ4Qm(B9GD@h16rvwe}fC%-?^v%5fS2 zoL+gHUU5i0K|=_sC#mv~YLlnTNki&sDhsJ+D58*h7CH&3=ggZS^}Koeka~eS-n+d+ z%N$?+_^I)~Z-KZq*NhnUN&EyjyE5o1DQ|3{j!(&lI*f6A<%q}DyGSij|O?)oRsDVR7UediiqY%mdKHhk1f%| z`-J+jv^?menbPHGH#v%Y%aKi;)6#!@nc(s`pI*Qvt`F-wlRg`}uCAX6@eAj0mn8G^ zhQ5#M%M5>aEn`B*M0>S1XwviT-Iu7!r~FaHpvh;H$)L&Su;iscGv0{z1?+YC62~H% zzJk(d`Wh!%3fQ8`x_kpu9!=lE=C{7YVmoY?=!=2xXip>Pdz$<#nX6#1BeLtV{J^9r zc7BAe4gNEGZ~Ra6ic5l@X$Z0N3soFDBP#rrznYWA&Tmu}JHJyzvGWIXQWE@W-Yg0J zGH)MD6X1OY?y1+b%i@ayEf(B#*Ld_uqHGVw?iS4(DXh>A=_mx%1t z(B&uQrlD6%c3K(&+3BcqvK5)$oRsVgR7Q42iiqq?&g+f+-J=xnV2b;^5voHxF%TOw3CYcZ)-l(!FvABcJop7)3#EXbete8?ww# zFE2@1*XOTv;Mb2j zyDY^D)`kQP6?w+kkx3Z>9F*)}6%OGBpZSkp9iP4A3b%5xA?^5LW$bukC+u}u4#xs} z2b5wTf)mAV^A_L9@-XGtSAfl)Yd{^0ZR>zkYP6=zS4h#8j)TZy1ezvD)ds) zd?hRAC0C^-kY0@{C*2^O=A@*%sEqXL6cOn)pcB$-nm0)gHE&OPE$a2smYi*y3LuUDe0nM$CW07lySA77rv`H$1c%yI!O z&I`7(4L+cZt4X$ny)N6~Sm4?oN^y1LB*8VKjDji0wF7J!t{u%(;1ZewE@UsowG*uo zuAQN4V`upEg}om7X)pcbvXi?Dy&_y4vMVhC*KSlft|r;toD|m{REBF$iU`+U&Zqoxl&yF;seUKnq@!O>#{$N1+D|26xV?`QCzlYV^KK>rX1J7 zuod;7)r6QpKZFTsnRqBo{#h^9gVEc_4`Wh9?r`YZ*x5$BeZ*JSj}GLHpqD4dw~=?q zk+cMIM^WYEn&oJ7QgX*o8M$LAB67z;C*+PdZ<0H~ygj)SsgL;EjHV4wKYQxCmcC`_ zDn&m{QhpnG%CxwRYY^Mg9-EN0so;oCdD(_n10d&!A-*G}oVqsq0w*7W(eH{^|3ZnNp^AL+Ur-?S7B4b5@| z>~$H9V-XEkLTNNyg%gbi8;;mbu7)X(hHGHcGPpw@+u(4-fDXBq)-(pL3yyrk6Hj~h z9oE?)*V7u2zX7_|p1(NjMtW(MZ>qQ@8aL4r(BDjzp~vsL#heuVtyG5oHi`)S?a&GO zJItHt?=)|Z{x0h9egGY#<_bLbq=i0n&COtyZABm72;Fm!n+=@EJjCXS?V*Gd)+%AQMFhMEr4jHFP7ndbep{pDWtj2^ zcm+1Mf{Hy+%qL*eEv~^sd$%uJQC`cd^h<-`HJawaUk*)aJh#(X<&_-KE>0`*I-R1R zcmulpZ04Kva%0o(IgYnz2r=n{JhN%t>S7T`G%-_b8&6c;6B^HvNGm`Vjb# z`Zq1WL$lY7M`Ku-Yb5bPB(BNR(BOlT43Ub@?Obf78pG0v0<^Vsm>U0Dua(PcxNu9o zP#NgIyfg2DK??o1SZ>7*G1uZC*mQ7;-Co{+3pL%QQRBFQ8#(uF1NW8vHRJxu?ixh6 ze8fLOF_!%?C9`Q*?E48UrE;hVgW#~&rYNfGvN53KR11%xe zex!oNhJyevRY2&Et=!ij>k-M;3N ziDAk?P6As1$$A83h{Bna2`R|Qf-ga`u<>tRoSX>}$SI&}?F&q_%p6-ExQ71O*>Rha zUiiurRM^j6uHO)S?!VH;*S-&tshAp|PED0VZIEfqNuf?lWl*Q1h@egnoj{$zya{zi z^Y&0@qSkFWbEvMA&;b$@5uq0y148>Cr>bR4;q4G-jk}rgnPsevG7Ic=nH9$Z>ugYp zb#|O2Sm&2HV9K%ff-PWGQofBpClgYva|PcF>)cF;u+9Tr#yT&(G=CB+^XGiD1eCq0 za+HlSzd0$&1*i<=f)o+Tg`g9Z3!67lE@Ivu<)YL*qAWpS8z6=9AENAo&n%;ClD@Fl zWicEJl#4?t%6>Relr~g7t@MW}N4W%S0VV7CPLFZ`6H=4|gD)Eg7~ZrTT46$jvI>o% zm2CY$1HD{P+tQ-BvXPd6vWY53*(A;8q$pdc3}q`tgt84fLAj)P6XjCo?NKgGoffrP zx9cQ=tGl(>bmLs(CK+w9%jYPVTdmmD21C)|8iY?QgKU;%V6V$y91D=`PzrKcoFpKp zmgQi|L3Y5VAkq55NJ}@4oZGF+5GJG`m#2wiEpAB~`sv-+X{8&N$qG!0P}ZPpW2eg| zfLElKcDV=g9hb|>N;CwLD^ukpn`ISqQj)7u8OhZsB9fiZ`R?#tb^&3oZWnl%Yfx{- zdJ_#LHm~pxYEQSL9@^~g85nCuVW%x!;en}sLhjcrN$cUVxB}L>d9LW#uVnPF42eI&XBVhglWWoZV zUSE|h@i}FHt+EyDb=exn0^l}K3UFJTC_ozsTUNG%DF?VcY{f`sm1CJ5Rg|ijGU#Sf zif|N7`T!GJJl$1JLR8gv2PQ-acZ9~^*T8Eo_= z(p^hjHfGU{OFAoZ22&!)XF}sa#H?|27QM7_)KulP?hrYfmVo&jsthxJ>AB{ln9rj! z%;!@?m@j}%Fkfii#C(x?d(0P8@7+^Q23DO+!qIg#35HLhabSZz7|qTuhv$g7^n!4~ z&{K8U*s*o@vzaQC^KuEkMHz6DTnc+#F2k_^d^wZ?z5*u-*w$*lBBNo-0bdE5u8-;D z)GyaeVT+loXiUjoO%wm}D9Aa$fIY}LopBA#5!P#=Yo%98Uq>&8wQ-QD~YW76ABt{39$VE$2T?ETLMUnT= z;?*rrK?xQJ!#D+Wj;~&>n!Ujm@9UdyrlUw1;R3&>p7Bp*71R z=A_Ubr7~!bQAE%lhfbh9VcvxHq>GdyC4ry-g8udj~q<_O5x8+k58ixxG)VE7Y0Ud=s|E zRffz_Y|0`NIdH{1A1D?zp1sB~p2hnBpHe2*Dj&jLmyd8P$bAf@!SX3g zIl0eZb2vKs;Mr>%hW(u86x%ok#wHDTg)zRudXDJiu(wCZsuqHW5t; zG?b6Jj6jFtdxa;h9!n@nL&3T;X%gEkdK1Z`^Q z1lly_O=#1aw}&1PVXYW^v{rv>t|W#v^O;F;y&c#%GkE;hlvuz+RUbaV((C z1f|es#z~UX2FNTh<?e? zzt+Dpxy?y09VTx=g9MpPhsm3;!vQniTucoN=cdXS;^)p|PRejzDq}bwMZ~Z-bi#0c z^CrUu%-b_ukUB0e!`w$y0{-U^UHK{t;Zw>KE3z=`by)<*g5si3O0f@4lA7)r=?ha% zaWUAmr`IvmhSs>IuuB$aLdvlpP5LkvcD#*{_oqGLxCC@<%n#ZA0R!me^IxSmmw^tX zC9tbd~&ck$AaP-P)c!4oG3+`*_c&^!jw~73pQU)>VwynzAh`SO?wJ)9h&%O=6rKe z_^wNH#E`4oYGY^2mY=OpuY}A5QsoRAWSBWA!{JoMa3hL{;Rxu2;l}1ohMSnT zXSgYKT#9nwRPpJogc}jp>_x8rc_!a_3H4^$>MWPJG7_IsCf6vN!CsfmaV*Gf0j1=& z#EFu#2U>^9Rxst{wuY^Ez=gx771@Til-jm5`3GFwE(>_TWjh)pX4^xTA7$;PS764P zZxk(o*$z}Wvqss`oRpbR8MB=zB4#^7C(L#+Z!+7}ygjqss8?m_$=ovtB^}B5K{Oi> zDJntGN;Uy2x(_T~QBvwCx9pC)mKin49X(9Yz%ZBsl${r+6T7M zAAesOQ%L*KG!=eXSUO*k{b`FB9ROV(_B@bYff3hr9z;uEbTC!Us7VenCuMXfl`%Su zB4Ttnbi(Ke^CqJs&D%3NirO*K^>{A#;O}i+xA<>QLNQ4lP3dqQa*mQ|)y8!;N8_$# zM$K{z>~%R7$AZyuP|E0doG2q3;d(<(fGKBmB5Xai6i=c(#dI=FRsw44Q^!kAp*bRT zDs;KueHy(IQm4}rNS#5IlWLYT%}Gg}MP;PUrie(L1D%jM*StyUJoEOX&Zo{x!13yL zyD11QFHr)rykyZAaJhWxSI=I6Psvx|Vfo^PFxTZG91Cz4Ln*jRaH8O>+r7A43R4d5 zGT3Yei5q5j;So?>i&>S+nUu1-f~N9@+33oRW=e$bN@%?0d*gFI-m2jt{bRFBa236D zwH-Fh=B~+C(-7!gL&c{lwjZ_KkagD_wc7};(!ADu6y9}I2Jd=`2;L3Q3A`K4oA7Qj zZx8Qg>iyHA(B>D=hUD9g{8S)ke{7x?+0d_ccWWpypxk32$!q#BqZZQ(5aAj=PtIA! zQDiT-;H#7fH_NRs*X1@G3&OWUDd9VCqJ*uFzM$L*Q%?9U*kZz%6*9tiGbtrJh9<7d z?b4-{T$Iaa{$0OsEE6Nh_du6liE=N!V#xQ=5Fp=Al|gQn2h2wyKS*VeAEJmLKMdW| z# z)5|k3<+PuLEmVS-6~a)OJV#rK`FWaHjuP`ud=~8mS|YYDLf6K;f6zN0%`oKKkbz&K zmt)(Y8Te%y0^3)p;^G&}WgD90RdZ6fuTdG?*C`^nZ&)Ii%f4xe9^zZnzU`t5*pP$t zE$)63OH9)H?FKnh;GfP#k;oMHMfciEhWJlJEZ}YCj`qOtP`KJl%Pv3>5gsV;7 zgSjs6<5=MO07`Lvh?4}@{PGb@Ij)aki*d32`3Y?)u1{%7+Mgiy8I2LK&!KB$9=zd( z-b06-8;E^DFD}4oY$GuDiXzAEOC|+&Us1*RvQt`Gh1S_jH6FxTZ*91D2AK`Ff7agxBB zPyT=@hxaFJ8tomOtU8g~v?cQw%_+YLu!)ND>*&U#M_2^Kftm?vjQC9iT^sXM_s^4# zxMR)0Z(@4IIc^df0<}r0GHR_dnfWNS$*GLm6ciD)DJ_xIrm{p&ZE9*=Bg1hMblKQK ziK@sB*UZE}rNK^r0sDY`d~?+y(+NM0&J;8;+c z5lSh}gcGG`d*l2fGsBcqoCP+kAEPb5F`O;4GASiF8%-R9VZj-_Hez<#BbakQ*V<=W z^v`D}I^~7{vlqP_%tkFf=cFOPoQoD3s#vgOeoF9x8od%5g6STMxKdPcF`+6n8(G4EIh&%c(!@5$+|RYlDBD z=#O3Aon-rfdjP#++yiL{xGPjS?l!5Klj3fmGTe<65$+~SB<^NQ^tfB7J#JkP1Hzz! zI!`ksZ~}J#O##kqmhufQ<$9DYS9b@^?XDhG-|(PP&$pCTevtxq8)brfNmxp_8)PY% z>#{VC1@1vmihCKHB)AWj!7$~x+hOYgH|yDDnUvyQE_fq%v&~x&?_g3ydC1LHh{8JY1QrQXCPv^lfUXVx zFl%OQNUs?DFd72*;Z!+z?CW4o3VsBYf!~-S0>6nR68xr?=)sSqKA9yut6j3?O1P+^ zb_27>;=VJAl^LO|$HGSD`kvq0!h3N+9K(oXyZ?9#79Z*Ub`{x-U#|#;%_);$*aDV3 z7}~IO7UsHag<}y6TSI9uY=aX7!)vqH)4>PIwlL+vupMkY?~{od9=V+V!*0*?N`s=C zrXF^99~EzSg<;@P^o!zS2k7$aVRxjL#)seGU1$jLu@hC?j2ltG@88*cG(2{pvhdiI zA_|Y)pnK{R?`{{dQ@n>=hJRgL9w6`xDha@&rwgfi0OI9cHPF#&wsN^h0FZ9 zcwg9XAO-jE^@E$5LuBRqCu!aEZ+6Pe_jDf5ujpxXyZ zgKl4(XbI%olvj?F{b0(2ZhzR;!+Hrk-R!no)R)t;*^j5wyW{{Sr_pvGP52Lc_wtd} z4s5H|frWyD^aPt5$U$_7!t7vZJj3+P{lCxi-M~=)96~P_X06)5Jd}nIW`|K_#ZYrq z4mT%_vm>Z1&W@yr;_N6(WHEHKCHiOJyThO$IoeBJT?^r>e4^EN* z!}?{I@&LOJw)N&=p@MZKJx593K||$!x}_2JKzt>Op`JtRLApgD_7HTfeZEUxxqH=@ zRth2ZFuh!eRpb$h5Mqx~!$Tj(i9-Ka7<&XxxwqsxECoqp8F|Lu@h9dS`=DuK$qWp@+Q6b>#75*{d!jiR^%;OLVUeV zm9)OA_^I!hkA~N~R94@;M-heB`_Rd};|J!=X7`8Y?c?nu>QgM<5D^`&*`iC(v?k(f zxHuimMt3^IlJwtvDTdy=ulW4EuyPG?>*tcSLFLDM#&Iff~ZJq;lO zexOPt0KfQ0^U(+J_3HHz8T*}OPmU;HNa<4;YMM? z=n8>ZzuOl*zK-KFiDVo`?D(K6XH*(cW>`C@xbXlT5kOop)ZBZzxmyg))&h5pWX1tM zQme^ve5Cw=?_Y|L2Kf`_y8MM>5g`-I;Ui>1oM?nte|>+M2&OzjCWb8x-f&uFP-g^t z5~ijhGAT{ElN9Dw`oI~}D0V`3B>lS7wxhMIz2E;t%`Z`yAYnUaPO98*zc!OeM@PP>a0_ZTDt1UF+BxQ}LgZr*~}bG#sA=rTX|_|flP*$_&GdFZ9^o7|d|d1(mn=c7v79u3WH(%YPr{`^!%e*ub! z{(_dsw#P!2=qWEuy^>Mp%46-W;)f{buKwUuv`Jn2gnNA6cM^MOC-7_EzzS} zih5l`=WE;XAsAf$0p!yBfCA(o$^_&xu#`YH%3zr5(vD*Ra#<(^xg1VnkXKGG9Wdn} zhrp&Qy|}enk_VS(N=kACn%ri7Wt250MkrT=#^zYrRM|@OiczjiLqNF-RgMyG*EAtR8%nmKiniGJ5Y4an2e51M;|j1tDHE`3!BPU-Bx}Q5 zmvwL~z^)6WVAsP*40iM*vOY{X*bQI{jqkoG+6|eOq8&z)AIa?KTgs=ynH<602)foj zO}6ZB1ifOo8`BWrZbFsAt;(k6q;N-48Qje%BDkAdBH?afi5~8j)LSNSwMpe_!oL4p z;ws(vl7GOv6+f`RyESEkcNd;q?Mwb~HIXTQ=gDINN zJvWostny>ddMo}-d_@~Z3sL_v3jW1SWgmXQA}sc$Ou}M6SW02hBKyN!mjiGt!s0+E z4U2+=hMsKZqW&-3up+qFQm$Gx5-83q_{7pGTfI?M7S@t#Et>W%Vn16abHf|K-{h$ z%*6Mr{Q=@PJbAFc)pQVTpK4*zU_`D^0gfCF9MTEs%i5q(K<4(Tm?8{b+Eg<{F&nmHO zye+T<&k>+&^@MF4yQr2+6QPT~N#a%=exraS7*D~c37GG1CvtvKhngV z4m#+qr80`-0qz?4iK!9&pP_5*3st6HVz{moK=d zrw+JhpvnsF7Mam}l=e(iMtf$8i1sX&$hFE^2cY`)OV7 zxrZOm_9zrO3yzr7T)WE;H-iX&x78Un?os4|L;vb6aq#X(d?aT$t; z;$Y~6V!L^h;FnX0VlMPbsbvxyBS%r9Gm!8gzN@?@oFd#ZCK-lr9gB7hk(pRd|@DXd0s;1rOg@?Kv-6 zW_3j}@4;aSy9y&2Y<(dwaxgEi6|1s7taaG{#{%MpPzrGvPLf*jNEr@O4sj#c+;fOl zE26#|L1RjA8T}d=|b6#pIp?E+fydBWH&5%&cqs( zQ83qK2OJBiJ3=W`!HMQfn{hu_c7iE~x-)EU&tPmxqeZ4KHFse`3U*hT`1j}IZ}=aU zC}RoyZnQ>ZcZV)7f!~8(nlW1|m}MU+d(sfl?nM=&9npXvxVQNz+B%h?-G?GVyDxN- zG50fXX3YJ~+p|7^dP~n*yEa@(sN+RGS1BbFpFI?m1Ra(i>q-eYkF-bVxW3&rm)+dF zTW)?JJ}@WPfOXw4*X3Xw3xbD0DZxW=qS?|OVLVd~gDEF?IBeSg&|RWIrb;$^WkMC* z)gzdiVmvZjfO@j<$Ijur;3%d>IFE*|wNIFBMsy6l$fDZXB*#(&l*dtJofZ7ltuSIqr7~hhP+rLau6SwfB{`dtXptk@Td2yAF7H5-5Rua zit)f|=&FJwQHSBDe{nwy9a@p??qZlLSV-p!T+P!!tea4UZfh|d?5}4H?4rt->KwyBjJY>!@p^_$xL&U0_b6yzMVZjP z8kU@Pvs?pfU9QEkpnV;b(!L%iG40V)%MCE)v~PqhOyrVwg!?9@rMPeI!QFV@!nBC@ ztJHQ#y0eN^%HcB!sAKk z+SrMX{&tJoR~i(;<0<=$hDuYv-c5L$o;*!ch>&NfvRa{~O`bI;4Uy-lEJU8Ch(hEA zOJs<=Xo)^VUZNiB{v0fn{R%2HkGc5O^@DVAHYzxatu_p8GIq3R z4|$$-Ki6A@5BP8H?}8Ig*R#^n53LNBm-$PIpn8Qe3946N$qSVhc@5UOypCfLQ*S_N zOudPdIHpEVEN{V-$JE=fjiU;Ahe>G+y-SlW4jF)9Ps9or=&~Fdpu;ThF*%Bx_n~WJ zKCL`?&BAl{3UTv+c{T6_j3ebk>JT&^QDwzSgM4g08Z)0zSmnE2mRBN*b$77DOFam#JX6$-sWjUr$2yA;&FqD-iN4NFeFRlb3>F5luZA`)Q_+<#$u#Qj(3+Tf36>_>MlJv?y# z&AdGK->C!lKd3V9jq<1YDEGgpjQa$0`VSo?TW`5_w@gUQL~@_V5Z~Er!v2<>4EbyH%7ojgTgsBz4Gv(hO4`AuX`^mER}^qBSh7Oi06` zjV8D98@}l(mL-`M#llk1*k@(IT~FwL(lH?xmZq1Bg{rQ57(_#eg=MJnSg6WibJAF7 zr?OaBmLiIU#`D#1?-if6!t1ONnkh0sxalSSA)$p{2Xh( z)<{R3TuIQ$

=qO_+h!hfJ*y@E5s4V0C&#;Me#+=H3F{s-o{6M+v*TL5zz9=iJr@ zL9cSexTx5Ig&nA0Y&~{&cQ;^ngNfbP-QAtPwZ1d6r_RO8Bk%M7zn|yb?zLw2KKruP zd?)v8;L3nIl2*gI%w>c(izJ4UD1={+C>6e5)_0jy_zj4h@EZ~Y!ZXTgB8A`BC0gM( zA->)TZ-aj=#o5E)Hn1E^@*`2Un-L61MZHN|VMMo0z0lbkt66-PVpBf2H_I@>$SlLb zWM-+C%|KUWa~!={wg6JIY>AVo?zl*{0+pI&1UPLVaBUOTL~yL4j24BCHEzws*ecr; z`86vt4#CT^HvhIv3|z7uFvcdY+v%|GuWnX^hwYimxWp``*nvdhk{yXsmo&&uE|a=s zXCmj4T?hh~jC6@~$*wNZx@0%vv$fWM4T=UM&A_l=R!0=`a;_z4BLRjmIW3FY1`AWs zII!c2tT13Cwnxm$xMJm$z0TZ5$3fkXu+_lVj5#XCo7tW$W{7s@6M8!cVPuCJZS?tGbu`%|a^p$S7;s9m~ zEK&uo3`~bB4rDIMTi;kO2N4wZIG8B02X-+%#AQ-@9O|y;@bzH?fjtffj$0B(Fvm^m ziD>Y^wku%^^bfg7K{z@QHzjQJz+0Y^F)T5eH#TNX#v_j8vwOE3MHso|XfT;u>g5>F zRXG+%@0R0$)Gf#3By!7XassH-EhmE0dZTDT+MrN$Nls#p*eEA=c##_qKZQ90kDLme z&sd$tT*e~}aymibku!)=k2J}dE|YrXEO$LUayCKWk#k%kL&&)<(T0%oh;_U(9&3v1 z!gB6n3VSMm9gI^a%}vAt9TPL==-@WS3h{WkoX_X-R=9vLvciR6GAlI7MWCy4F^=8} zmjJ01F2zY`g=>zH%Rr@8xE!1g3Fs=jHY{7`(B>6nWBIQn#gCvOV_!Q$(U+^-uVO-= z`qjXdflDR5(b3FhRBw_o1cmBjiBi>Dd{j0C*}OeScH+zPrX zx8djoy&Xsey#psIsGG0aLGA>V3VIhfzig5X1Gn&VSQ*R5%@r`0qOqHSd=6CnoR_y%QcKjGXWonA8dNt@1eNsyuV&6o5;@^2c^XveglE9nml||o0W{Swk(cCIro=9IF3PO=M&N=GC3&7Hf&4E3 z;~Nr74%p+^xz1L#IFI!rbE*6-{j~2mkwhW?OGJtMcx(7&mr3P+g~-YODnTIsYc7#x z{B@UT%lI3_dpIZ9L19$w>=U8^z>l$@$Dl(%E!F~@hu|KJ{A{CHgvszw^HTlp6|x~S zdkP%5(!Qg)XGz}VuX~Z-B8)_S8%#3D+|VHJfUU~AIC_oW15%CO$4L~v2gnDYQjI?Z zr~k_MzVQ(Gh$zmL!^90b48XTppw;K@ zP9wM6WgGdz<#Dfev`b@IzI1t1&##D_o?jCLdVb>)spq#Y(dzjfaUZ8AghY(xJ4Hnw zpYOA1W2qrCQm^}XpbX9)F(W~$akgD)2PNI@_n*0VTRpEDV zI`qPUTK3s>^A+=Xhd-DWYw{;4zHNjae!L&86a)rS262;Rv=x)vD@~FP^6FGesAPDqb&?Qpegi5JI7B=PRD1gKQvCBf-$VJt-(KvG907Pk*6el>!dGH{%go1!d5GEjDDV7!e| z{qXwlE*VwC)n!~R6IW$f;zHHsh!Rz?*4yP#Rr?V+RZ9ecs%4i*RqI`%RkeY5uu~Np zvV+3Fssgb>!%js(faR8LI!xB#Yk<JQE;1)*5=tw<(TYCx2U)!#reflMm_<9o)lp4qtN^luh2t?Y7Hnc9d8 znN}f6WNMUET^^Nb5RsFqLJ-KbnoFcigI%JPX?5a3kxWp@v^cnxR9jLH>dty@ftJ-! z^Vrg~2Jh!JT9Yu+Xe}_QMlI3~wkkt#^ct-Vq#CV*lTf2;wv%;1r5X(dX9@}nlPyr| zk&Ihj@9&k{iZD`c z1ejF0R@oXXzTb$WS8iJ%RcVqmyLLZEMSq=o1`rsgN`VWse2ET2|!DM3x970Mz=)pR1nDdJq zNk0N02Q%aHPuU$U|_t5j#2_uD10+aMy8v8ZK$)Kxp3XWdlQ-M_C({Q2^yDbtnl+!__5}yIC*onxj zM`>^OJNTc;vr>lpA?Fp%UtZh zG~kPc{4(fyBnq+5CyF~TP3SLQ;4-P&7ZN$OFCqxkz8HANiBntBs(cAEUF%AExD3V4 zDi^P8p&aWf6s|*O|3=rr7i$sNSwWY@&K#I|M@%u;!9rcCfT4mq)5pQp9DOijg=QIw zzPZ3XZ98>axs;FX?Qt1lWRJ_iWcFy3D?nG}N*ujCt^!hfT#b{^9%DX~(V$X$i~%#V#E9fMBX#)LW`}WJVOb+~U9dKpb_g8G%b<>G+6#lrLxx7ExN6*$UQ)6kb7|w8RRUv4^(Q9ao{?R2IH9; z`{RC6cI&em{&;}Nfj=Guu8dFLl6i=^oIh4xrcE9uQTXE#qQoDl5+}G!>W@c>oIf5T z2>kIlaORIEm??@MTv0To%6SkzrmQ`lk`NKaS=O*5+xmC-`-Fmw&yt-g|_R_1*`Q zS+7Mt09}<2arD;v2uQ8>F;3KaZq@D0@(HNadY^*Z2y2Q>5SM{E?0B2HgdiuE)4e#QEE7>Ja?JYaNiL*!+lQ> z814t)ctY$)=6J@rEpFlHsy)_iv`hF6sgYZ3MatFVVQww;j<+16Fa}RLm5ZXlmff zfS*1;{ndf{oTZ|3Uac#0ac8k+jauv+HMG!AZPJZNg`c_;rG9Fa9xjvmX&NHur)dcS zKTYQn*+rV(CECXG48%{-TDFU1#|>@3;>uc#@Y}UDw%HQ9NSD|oqsX?kOUAaKwkE|CuE=@PAjdJ!Mx zf>SZCQ1r6x#E+@ptTBwA%m@j@kH=Tcq6>9MuoR`G^4+X zv@v5~tyO^Y1%9hCmuYF&Yr|&{iNaeIqSRX@SaDR;z)n$m6@J7$RN1S&P<&fq!~>ULqq*q0+oaZI8UzTB0017Ge2Tp2%ovM6?U z=5oGdEi5DoUydS5eOWJixJ>HHJ&Bwz_aX>9-lDhJ}|t$7fTTJvC>MAi4hatNr@numf5`rXI{d~PG zx7vA|p7r%B^lVsmqn>qZ_gO(Zn@2OfaON>YsWThpSeHqic^r{*=J5o9Gf!}dbmoaJ z(K_=a;-|EI=mQqQP0-myPci6iqU-s8(b?o40~`f2-AZ&>@LjJ)i@0_M4}^!nnaH_!64l1L%2}XNcbyHcqYj60nw-Nd zvBAzI^-l&nkC_64oe!L^55Iu9IxrYR--S#s40aJwYOp4`*kw|KT|(pxb}2z%u*+N` z4R*OpvdU^0(2$!O438H1zu*jOO-*fltb zBG@>&7F6o7>%fIAKKjb^9!gjc>*Cn;%oW@022xSYI?Oltl5ln-^9F9a2{?bP{buH> z;kI>bIJ<@Eh1+f=O5N5Xx4BH}w%duE+wLF;+;*o+q}%RtiPmj*E1qNnkRG!0$_Qi# zXKZ}=R=*wXw>=ZvPeD9F7Q`RB##YzOnlL!|<7!;1fP2;H9a1;+VnZ%Db*hQcZQyH!2I7i^?Cl2c&ENV7&-MtFqudQdt)K_q#PIce7I7VItl{)n`aCQqA`;}I+2^3~Mhsf(pjNSQ02l;i^s5hAy*zqmk z%J?Pk?cILoC3~w~IN9(vbMen@z;1La$vcFF``#r=+=r#<@3~CszW0fo`#vBD-1nhN zr29T{iPn7|6Q9a9AJ;*$ZO5?EiowBl%(9~e#XPERpcvDyY6 z`Y-QPER1$5EI#40c^iC67}?-6FzGBoSw06_l`nAgF8C5iUGNo7!aDVuUh*}l)CJ#w z(+FV)^~|8ZuT8&YdTfO6NO9(~#&DjUGx(la0z3QwTp76F{14wh;MwO2JN(F8*_^>o zWD8gP?DF_I0}K@Z;_|2|ekF3I_>CYi#qTbWruf4pT2uT<{M2N01Q>SyO-}%q|I6WM zKEXdsgUMg~v%H`FCXD>l1>5_iVFlanOa;0sQ{(6j)fGq$)eR?MSQ)d1bO)6hst34@ zdg7nduB_mM64z0*V`+t(A4j0Ux{7I-Gq%>Wr11alV($H4*X}R0tGuRTw!mJ~16Rh+ zFt~Y*N3Qy&u-6RCWkx@mw2wI>iNaYk5vB2_L1uQD)LOF;Icv>I5LjzA;CO<3cINn* zekwu(C&CqB;+iLX-5cSA(S}ilv4HbV3_`jxLYrw^?O?`*XsC1c8WZhWu{rn;;GH)oVdT8Iz+}#Al(|7y zWgZ;8^X3Io=go(csOj~w%nvGc-U8t4i^g^8IJWW8L%^4%+tqLjGDB>*g@S9^O@0n- zm>sog(uJ8Ju-YQP`GUfp%w??Bs#fboqOjVcM5)zUrMJtZR$Gk7S#5EGz-mhX$5vaC zIc}(BHTuC=N>|y~C1kGo4@05(Z)U`?ID=yg9BYOA9C2Xv7QLmT6{WfeIsW;6I#$8S zQN16Z0v@q*v97kOLpJEUA(f}yV(m((*oFh~P%W!#U+HoT24mwf}FiMa4bGg15oLQ_dv zTpqPYE0MFu@&tiBR&a@|r{R1*hO7BSI3#?!;R38(JnQ& z!N=A?wkv1)mY6PO%bJxOb=MVJ|0lm9EAnsh9vVOxd1xS*)I$xj64bth0_wq;=MHiPkzpiSrUg9<{(WM9u=+5(F054!CmAFZeC|U+^!wy}O8?V%Wi5 zWUa9y@y58omUpuj53N%ZGjx>6G*?W(K=g=u-TZhLY@rpGfs< zEXmHGtFjA@UcZq*s^6|S(L(Iz>6Vt=K&ATa4o?5g+Rb_xY<~&~g=}1aN0H*TH0@}J zQa<;l!5&NqRNNC7-?jN{;{FqU-oH?BFXqyQft7A=5`}d85XB9H=2D~V>oTc!`w=E1LYQF&b>IIx=|o|jkZr~iDGK>xTzoq@+W;s zJcuwVi3fwplq}03pz)Dp9KDi<0jZLQ<0Mpa^yhK}s8q=#!I@=l46W301lED$`zA** zAr|v!Qo78IVj^3h=P^tR^gI^0GJfu4x#V%o73z6Bi9*j4h*CY93K3i zpyw&Tv7V4cHKXMoA{t(P-FSLG}my}oAyslMmnMD=xVkj*3Kf=cy051js!+C64= z?xWH6eA2PH7m(t%7dn)IXp7z>kg=evat)4N$!mdB$?I?uDLGWG z2bC&$131>Y=-dx64JTV#-FG@iBDfBb8<{Ou_9jw2-Ib>7x{G`f=EXNNXQ1>gz?A{% z?)A4aSE%%DBnqW(CrXuWkvm)_Rr*dMr}SL}fzo#a$4cMB9A8ptZAocehuK}!jcAl} zj>EnUfdZ}u&>otZF=2bU9v)EPVG;9X1P$JF1^e>6fhq8&^G3rIt*(Y*C4EP~YBsJy z6^zGe>(u@GfL8 z+XF!Aw+C?&)ePUtL!eT>Jq*s)@pW#k!w`9dY;3Oyq^LKC3AB->!v&9$4!rajFlvO& z-#PZ{o5rg29MgWBxk4{JL89=|lSHYPTIDI1Nxk$mk@M0s1c8^H1&+P+9CK*YU~LWU z*#xgMi`yZ*X@Q&CUPpUO(@eW%wuV>3Ax5Z-5FOIvDj!6W`2ReA(5v_YVWi@VU@{e( zWFlxxQ{d=Td>Kepd<7?wioeUNpi&iI1E-@p#a@ylyqYZXuQNj|=o_TCZ*g%=C@e2} zlQ{xy-vZ9t;BPZmsO>u>3T@vdO0+Fw#=>P%ZQmzy+I~O~X!{{>tnEk4ahQ2@mTvn& z6jSOE0U%bfL9ywA>+X@JjR%Q~-GkZKjy5etxeoP4LjroE%&G^dKAfcGx44CT%!l#1 ze?l1P{wbJD_h$JFbX7jb(d+&Nkm~*=PD0(sG|N|@Qr*7>r!S@qf!=j(JNH2m=Fl9f z^KY0GEB-Ahw(PX0#q74Xh3wxkEztUV;L7+}k`H73z+9o$Kawc4{)s5nx+FilOse%S zL{95p2?DKu1CF)+ojFFS)~5Gtim|ECcSecSm=B>`+dU}FTAfp4s0IabrP+BHEkA8B zsH!N&MZc=mwiUD`fZc@u;E#K4|0ImG{R>Q{ZHxR3x++~34z-;MNVT0BC!w~Zzm%?^ zQf<3|)5zQw_P)B4jb-gY%4`Z6_r9iKLZIoizGINW`yY9#E->^MbP?`hYC+k&Q*1pA`SOIql=C*9Axi5-tdw zOSll(T*8IjJ^i#`bGcp?Azet=lgLTfiy)A2QQ%C%-pr&mE9F3}EoxTOnEzO_+H15h zYWc%O|4_3o#vk+oE>0K;xCEF~Kumxx3A!qEIC=s50I7gW;UpArOsgynDiv@UaOMM) z+H>}=uWAQtS*FH9E=P*1N$hY|T`{3&2DbV#Jy5hCa2{w&WOGH!yf`hC4N_0KP_%)_ zDcVR7DB1*^Dca0Ty6J+pe;_D3L#RmEpF->mUs5%rmBB;_ve}su)RUMhadDW&U?;hw z5iQ;qcen4p{^MY9!zEoyyycOFZ$r)~aN4V+`WOsWBR*Lb@=< zYDCT$g9!p-tPY$RV-03Xi~-Sf>;p$MpD|5lT3GB!X(TA|jb{7mpJ`>+G>h4+#lK?8 zn%T{@=&bG}ueFL8Zo651jkZTu^hfLA^ecVpnWHN_X>J0dUg&GCi450Qi zEw!S@?K<4W9gO}c=j?bDa&lfx>j`MxhMSq*Gj3)&@v;RU$4k8>VI=idU^1!eWd!J| zY>lIrdK(~>dRv@?F<{IbvK^>Y>g~Z5Qv3Ry%I(0USn3@~*@4xtl4tgvm=*}VGjL`6 z9LbhEyO7O=9?6SSq45Q6(uL5w5jmlECkTWV;7sUI%yd{yp*03TkWTlqSnEeZPi`-3 z4?c{Sc~8Pf=DomVGB?QHpsTVEj$Y<{fmG)Ga1xd1gJgeDsmup}Gs5YR9BO_YLf3(~ zD$`<}4mYkxrVnCrAoan(`J&WA$mUWX%8OH}TjVg(h17==IjN5z2&6s|IFtG) zW;(@5jq=A&zoEEs{LIX0bEu4g-$rG|d=zg}vN`sOBs0>*h;E0{A85Ihv2=g+GQc68=~)nedHr9O$YXkE0j<1RxduM4UwBdoDQ% zR4V+*;Jom00Y8OlvGAvoGV|CFK4`|XgrCO5K=#vt^Adgr*0yTFgzNJ>6aRE*-c@`jiA1*tZ^sRA5{bg z=T-IMMshwM%!_^jVI=y6U^3C062I9Z{Wju>2D;Aq`wJFCVjKq z47w_};OM2l6-cGO4JT1~?<%*0N~OO8oR>Z>?{_jSmi{hMTHdL9Qr_=oVj%oIz?JcH zCObgfOEwq&K3<#(Uy^a83*pBTIpOao2!wwCI1~OsW-7vbqL+5`OdR*E`c4)2`#Kzf9_K<;OOE93sy;-XKsub)tq^XJIsazD?DQ@P9X0_j5T7m1wQ z6A1#jUjokLewmqc=#nr?mFX*d;k)W@E$pmec*Foc3-#B_f-vm-Y-@-{$j=Rg- zpi)8K0p|qurC8o29Siy%DQ=nzL4yJ;?~@PY`~WyFy&saz<@|^jr*hWI$D|86KOu5* zeo7F?`5AB~=jY6{yQ{68DwMO%2!>*7nrCJ(h3zrhdP6PUeRzk$73{6A(F zZ2F=XaRJ#>1iZ*DAe)-_NMO!@>*NI0=4`f3bC{*^5LtIfVLOUW2-@hT4^BPHoCiVk zi={NYhn8J7oeF zM&#t4ogk2V4&bZ^&&f>M;HDHBjM0|1m9D9PjXHy%FW!V0gs9yIN3*&mX_DE2Rbva@ z2F15{Nz1BS<31?rS&K~O;)8il&P^B<@Oi*w!k1)T&{dfaM=$*RKq~wKIEkY5JhC9D zRQQF!6~cQJbv1%4%*0suMMxFG`wQGaYEPyI!uJBsds>T<&4usHi&NnnWiir)@QV{U z;g=u?gkKUk6TXg_obW84rmUjeqBLuHtf3T&a_!ESPHDPzO&mC#)NQnYchL7LNgqC# z7k(+iNcg3}WWtwa8PHW(7Dq4qazHA4Uz|iG{srj=Diyv2&Q#MV^2Ml%lu5^u*OStA z7lu@#6{ZbL3B+y$u8jL->Oc2WVlPh+h`j=E zCiZ`rDQtH!c0+Pbs!U$HOD{D3=xn>IKY!c{x*}mD=m0R8p!G5kbX8Wu(F?jVkP6y{ zldvR@{##Z7l?u8lI45Y(?ixfo7PLZ2+g(GVcGqg;133o+=aV3-lg;H^gBPcAHp`l% z3pv*!a&op41ab}m&g5L1nZkBgA!prWwz~>BJKOGBhd=GbT$eBsb10Zh%m!HxbXC^J z(TlkOkczn>PNLHJsqkg|RLqURxz2D{J~bRq1PL{8YP2m)b80B6E(%}jZBSjRqD8UM>+5Ics-*j%g5HN7ST z=31kk+wjr6@Y@nb!fyvA6TVTl2VIpNaP-3O2&BUAgp;rokC{$(29*lG3pg*lFUOQ| zB-3KycO_*D@#t+EH{LiRn0oa}oM z1hVf9oXNfqGZmwr#o$FDyO!$O74IZQJ>&bZN;ZYOgPIqrXQQ6`^1;05`w>Q>?++#u zy-5xLU6m@1Ui1TjRP=*z5|#F^!5In&fVrM5TNhxd&8gl6%3WtB$B##gTCYb=7ndDJs zawgH|kq5mOAA`VrITnG_r31yWiaup<_Xm%hk^z!(#O@%6{;zn%V^+kR3#*|UC$xpcg= zNnYkTW!Lf2SBT>A(iVBuWzq4{*NB`mUMC2g@dmKBod2hM(_O^z(zo12cD(d$;=?t( zpb|it*Gk@um+G*o9Tbl~b6jJ*v?&=cHElVz@i6p?YJ=San1l@h`qEg%y3@UJP?mS_ zG-)7c=_l`kuF88jdXe7;QjtHvNfZcXk`FG~e6_FG8Yl1-FZ-BF+|1C4=;$n=DcDBNw z>~;P@mH!=o(7Wqh zET+}{JLy=&KS<@TZLTcxCsP6u{{pUzyFY1m|4lX*u}e?;0B%mBU79iRF%{`T#Hop# zh+PQ+5xW6rB6epcC!!ziG|NV{M6N7TsNJnMk{}*L?Icm1OjpEu0>!9&HKY!4RxBy`!;(}l@ z5$k0k&{bI&M=#KXvlZqPEg&df9vJ*R(o@n$*xv=_85VI*ijFqxnYQUYC- zGLBx*dLR|F0Vh$Z>?Mt$QbC)*6@vOPENLbm3)(`8Eg%XSE&^y}LLlezz?evy@5)~ z2EAgHRpl(MYh(LGizFH_P?d}EijfXBcG#@8>#SGi4|_4&2qQ690h5V|kFtZV${-xQ zm=z!ub2XeqrL(^b29=7rIyf^TIbPIKsFF(_%r_B&9>q$Q~C_Ua=OF16kXF zD+AJv4u+7;WnG&Wr?R%nI;0C(*CleY4kZX=T@N^ub$wSE8x5|k06`NzBMmSWXIP`w;^4~ zzAcfHeLI3c_U(Z)*>_;3le%;$JH*#fs-m|$sr7i1klhq=gaj0F*IDg&0WL(}kq_oY z--$31eP=M4=*_YV=&Fpw(Tlz-kcz$=PNLGjzU&Sv6`7iW49q3y?}=Ssep&$Br3JD$q}GZ0gnWi4n+@= zqsYbr9!<&$h_30NqCI>J69NH`1U$C1qiJf0V)0+!_j(uII05;*}+A_xRL88{R0 z6lQV)CY3HqDOy!EhN7uO&E~L*J*n34sr)%F;Aw=BfTx4W1ZN0jYpz z<0LA9Bjg-VsetE#OSY4UE8X*$A=dJIQuBAAzn4XvS$cSln)d?c31qzxxH2G}BD#od zF6+g-IF+?tE+JjWdMS~U^)iA$*2{r2S+8KGBf98wqd64rCBaIULTd=dk>fkix#QQ_ zvb`j9<`=hUD_T_;8=U0pY*+GeyzEyIMzUWGCLK3ymC<1F

SKXa^=ojVs@sqC zUQ?B8L8Yo+2Tp&hO?M?mm3cjxSmqnTOhZs$66AYVyVR_s7-1ZNgAkSO=?hF+w52_up21Cy4`)_yV$G`^jHqgUyE zAXVuBoP@!3^u6*Rs8pqgz?qoIWw6+}*X$-!mkH@nUH$O%jg zbbA!IGJgJh55>}{?uBlTF&F=~lJ2e8uUQ@^Qz-WYQCcb+{G!Z--v(SiO5A*$mL?%2&i;$jAPSl6Ah{>!y~> zwXM-^zxWAXbkPcz=lEP+`R55EUHaS|ziy}SY{RsL0Q zrZLPekozK}*_81$rp4mFPHN6BxHY;!@QPl5;M3MtM$vwQse$@$0$0W_n7ps>7TG)$ zzs-wN(HrF*(sYK3{w|Rd{XK#}^!I@?(LZ3O3%l5IYL>@hydF`whSX&>f)$VB?+D%; zk=GeSu5;EA3}|MHG;hcE6Eq~-`x!5UlRd^cc0Yq6R@<&{*>;6LQViTj1OT-;vEt@I5b1P0%DikSiA->(%n2%u z1apD2CdeYe+)Rs2Fb^qXf*O%vUZw^nm=Cxzeu3m;5%ZJHO|SqjPEF7v3z9BOun>_m z!NLTA2^Im)Owf~=3KJwfE*lA0J)k{Mv*6m;kw*f%Az+Qs&kmDGVsP0$CN_PWfny>NS?qG7lc z^TcLYn$$vw2x-GG+?J*!%P?Euie-WG)!56C&0W!#7pJahm42iPSCojHE6N0cE9!wW zS2Qrw$9Q!6HWO=n^?=S_pnm7}Gl&>ktLp+I*u(VVN@h%4KkEikzYNmpwneFv%}7%j z<;_Urtj<4~R#9wnHu4Yl_G%)G?9~h=vlrfC09}<<9KF4k2U2^jfRoT(qsPgAK&AHT z56;@F*u)==4zxD>R!1u`Pi(INq!zIDY70Hp?$@FWWTwDeD*@;IkCn;h=4#``iMg=j z*D9n7bFE6`%r%H0FjoaQGuLX&^kHnSI`|1XFw8Kbz&K`r&e@29LtHjfwW8vDnl?T# zw-RiHWP2#4|KglzeTh9VbdbpRQ5c)}!@A0(SohSy#qnlHYzEFIQ3IWHX>d4X=5Vir%ea~KWz$}`Dqw4 z-9wIsj*46Bjrp|qn)!o(*u+209+YA-~8Loe%Q%vsx1er zb|jmdWG7ynnxrf{lP*lM3z0L)NP@s5y8>q>*^QYdWfBAeous#I9yE4Lqo<~hTKw@+ zzJ?%nlrY`U>?QwAy3>7tw1eAg8;HC1aq1aY^BX0yJO2i65+RICG73y) z61>X+x+;6(=uNU0keXy~oP;JBvx4jcDmBTz;5M4q1O}ZQjzS5i{`O;z*eLsxn%`bj zBb@NE1DGc;OBFcZs_{Uwxmgb4#i?29p`W;v7~Fw0@UnOP2Jrbk^wLE*M< z8FARd1bj15%o41XmFX3`!fq=ZmU#-zh~jo}Do^$x?810YP;Kgs!y8bp_ zbuxhkzh*~Mj^JP9t#c${WSyhHB-Ux_hj&9jSLGNSy>*TSQtKRtlh8V&UzOuQrPetC zoVQLIVCaDpnIpE&Nu&zvzqO4taBzo zV4bsoGwYnqOp~yVh8Z1X(e4ymGbYa18X}gJ(E)0jID=})OlxMt=z%?1C7EdBC!)2f zu0QNrRq3y2xM9szKx&$saT1zl^aQyDRBD=A!Fki9@kVZAa%`I0NfoBagN@w5 z41sCx1V%Ti^z`+MFF2oS%@)mFWOLKp&5KjhG|4@r3)9?7{;{~~Xi2M8n6JP0N;4c>$SU6qG%^rm?PNKG>VC!uM^tRs(tN=@?^ zIJ2rrx2K9)H}W_WW6L~2YA(hZZc>$v_ddzwz%Wk%SH@45G=ZNcn;YgCUYr`HS)L_b z80I-5XPD;+0>iujoEhdtW*XyKHHb3VPF3hsYkS6aqgng-28V6)vNdBju|f@|^FL^z zdWYyuIfk#1L?ag~F2<9ntz_#v^Bct^5GtlNcGI-k(r>s-$@2m(ud44hfw6K1Mu2@M-i!w%g+ps=AEiDw!z=>2MM zh}~`UhB~tomT4B#Vz&L9-d39>9A&YsnjyHYe9C9`miUY?vc%_LGD|ed7oe;1C63+_ zUjeBlzQ##ti81TTH=t5Wd<)KeaGaBZepehfKU6BqcgzwS<9kvy2ZD>;l(2L_eqgr1 zAU^_E#!s8d5pQ(Y=K3#DV+Gq6Z(hXE<eznoi=s;^L zy8}ndZ2Wt?b!I1wtTP9g%sS07C+Moog`>C5+(2rbd2kY1XUyg@FR0Wy^MSLguJh^W z`I#7-W&u)K8F8POtbMW|lLNaf1Y8;4Bl*VU!en#1EW(RZyOg9S>B26(h@4#(B?#=& z8#uGeV$2j)M)pl%PCc0&ou)+l#EjogsrJd@d}eQnB?u!+ED0vFM2pmcu1X&qy(N|c zQcEn2lh6{Q@0VpjrIuJ0oU=sMK3R^5u_gMFvX-dPKIzBgz!D|ke5$BSHn&7QFHS8{ zmIl&=B^rsGC7K8VOEd##mS|z7nwHSwJ%!pQg(XTlgoe8!D(J+v3nP=$K56AMdrK@& z7+GQkFqtJ<!Wi`@;IR+CsbF5Afm}3p#%p7Yn z)1$6GgAi}xcCp0}qlPY2N_rk?Wx%SSbECBJZ%n!25;sb_7XKn|op!>=Izzyu)+zOq zwLw>99UQ%N)&)}Q48=(lO`ekVK&94MADp*N8c*nf4VWXg&W5B4>*PU&muN`L}p+gTAUb$4)jOn_Fj7UYuH|L57hotTUX*S!XkXz&e`)XV%$*nI>T!jV@Dk zls31H2{uz|ly*!0Mcz7F5k}S-0VcCfNwx-Em2GhJ*4Y+Ft+O3YqNsD6Y!526&JN&o zk(5~jKbul|KEjG&R8Hs-nP3-gwp5Zm`IO!Sdl5z^*c(h{g0k!bx+?qP z=uNO6keXnBoJ8^C894w{YJw^_Yl37I%V0T>X|V|oBE^+*G=be9%P(6wn5lsU4gt<% z#-U_$3mnFaQwuc7;iL--96{tPa3n!sfun#k3mnZ%=UNM}Ct+8QSR=%X<3p`Ee^=wJ zrODYq%Z{F*I@7*)2ftO8WB82T0LKzW1~?8(W`KG*9&}Ysz|kAvL?AW5NjQn3z$J1r zsMG+bfYZ9bL;=4EQVrhK@7?s(YBCFasmH# zZ?y{vBdc8mCbL?DTnxG@m*D8Fb}5ir?J}H1vE~W698_wxE5KQ+C9y`mr~9&9$vm;u zt|GNy5pzO+1xx0xX0E_!qk;2r+A(BvqmAXosnJ^G8q$T)t|f9tyN)0*+V#Mh(QaU- zZ>-U@XQIpXo!i{PgIY%6=1Iky%^0a?Ey5nP?OZb%dRxEPn#!tff*7hQjP4JCDWSct zc606Q8n&r@gKQr!w7AWSByL>Uw#m~!pe#4?FZVXPi7>L+&0sQ{HOeiZt8y!j-e$J} zsm*T3NfgnZkUK!7HoFs?cCbxE!=_F_D7%Y^vA^yn#c9N>VWay&mE;~~2&{E4a3$L< z|30$0wZ`${)LN}Fo^)ZY`-z;j9v}#;^&oI&t%sON7l31su4c=|e8w3jVmSKWHEp7= zN0;PbKBqUrBZQF=CV<8UYytf+nv2gy0F1SB4>k_2m%|t z44m2E6=re`T&=y#+Sp0AIN*m4S~ruoTUck?H?Q(3y$N0;j7;!4n9Kyt@&@RtyosYX z!COFTg12!JMT3d*4ye=w?}7`Wf$k$#baUQgVr+uP1c41c2hMEp1v5=D8bq6l{qLedw6j%7zT|UyBYZ^| z8R2U%nGss#8_-qx7DsP{?|{?@-{T~T1~14Dpi(3J2(ETCz>Z`;F)=p6&!l`bFgud@ zXz&Y@0~`DboJWJ-$mTZqofoGzD9azD3mg1N5WWNFd9tWj%A&V z22=4Vy$Pl!j7-oKOlE>s=>{4f(!tT2pa+neU>cl6(cp8L7F242>A;ycz06)_@q&@* znHsxb22vdLFgu!MGcq$WJ+Q${z?E$8{+Y?__r*tijv_EtE2f%bqlMB|qyf_|@QRGMyuPPcFdIecL_x+hnC zOo6wC=i(E4OUz9eSz;bAsU^z&WM0tt#14+$67vJ8B^JO*6c^r>1wo~jSO{EUiF5!} z7G`Q}iA6{`OXP#6(v#_dC3*o@#?O1=jSrO{-$<3@$mybFb4&E*#i=D4WHHi(B^D>@ zYAvw@L12j`fip|gF_X1~+3C(V7F2OJb4!b`BWZgCwkvD=VS6w&UK&IT)NZ$yzBuSs znbb|&|LvsNiJ{%oA9n5jg)2Koi;A8qK^owML_~kS?sz zO4QX_V|jwW8Y=*2*7y%IJ<4cdclYREd$nSC!&<}k{cNXGI~kLnkD0_CTvKNGs?Nt? zgRcJki@bGKB#f*x08D0`vJ3=`&-~!%t+O(aTBi*sQ8am1Rsoe-XH{_CI%zne2L>@m zY@G_J!a8|K;bp5aPhg$FzU_OR>yA>UkO_&y&W>Zo% z=VQ26SV@L4Juu90;L5mnlg*DeBbys$b6%VprdhTiT^MFdqOR63TM-0?83CLbW@~10 z4b+;W7U4#3|GREVxZ_w!w&7EH6KqQunP59GnF$(Xd(im65RTpiI|8W*cEU*%KPJk~ zpi&d;0Pofdn=)m#vAD6_?0IE>UHtjL31=4!P<4rj)|7DoVA#-*FQ z9!WO0#ZkOCwMDBOO}enfF+^RhEsiA!Y;hcLW{cyQ=^pd&9IU=KQ+&Jg) z;>0-E_w9Vrg>fz*>S~R1AwgiAi-0rZT+B?hjZI5XmjU{P#@7c@S`Fhgvav832TO6E*_-E$3d1YWroxH5k6)%O0XeTj)G zF?+$+F;`dIr=huCt|y=sbl<}ph*H1Q%Z)CRZlrb-QCI7in+XEH+ycD)<;vLI4sT_S z3CPa2rOt!t1}O+Ks+^swja#Q>vxI~1GZb5FnGj(YeAp0SxvOJV zzNK^V9C9zH)HnBmQ(JJJZKQsQYn-F@C@+hD@}0LCCF4nP`%0Zu>ECWPSiGNX;Eo4? zu_4806ZfC+^Ztc99%QbrxVE9WK^`I~-0?6`xGU!7P584Oae34n6NtK6b395AnBy^E zAEo?1;p6TienRI7cainalf)Zl{?U~^T1mMDBW7Y2FV*NAY_N9;v-Zch96BBWk~|H%D$n5P<$4xK<$4Y$Q8E8fo(Gl6^#ZtJb4m%WC!ILrMDeY z7b3sQr}B<_k1#5&RYUsG1oNefu^2YF0++T1RRo*)qa2jHwl@FO!FTU&I*QXP-5ea|S48|fgHZEkiD#aRuQfeRaz z1G}Io>q0-=9_2ca{3kw>m;7hKNb+C6WRf??ub`{)8;)M`-+@%}KX4LCK4w1o6I3es zU*NQPhs}Xdzt z6-Yff8&1@d?v000GCQc$lXHO6m*d*n@)x)HJ_Cn5=VYeXk8_b?R|VH4?>ykeJiV-s z?%_K(GX{2?2e>l6d-8qHd6~-us8$V7^O5N0EIB_>Y{`~>Wm&*wQcEsKU|QS!rPZ5Zst#w-FChH1SGt*PuU+vl)b*Wl;y7154=t`8DnJbf`7rkQ9 zW5>P`TDjF@Sx_HliaobfaUCx(4JU+AEDr+0)%C1Cmu9BGZ_5Bz#{G5D?!CsYF-q;h z?#{BzrG9H^=(*K0Tgq}Ix;ek~B})B}60*xbl4bu$QM=c@zgTWkT;VNg5O%qD%RfIqFPkDI#ro~Dk z+>LRJKV|O5EasOmFJOF(K9+LP4>1&SVXR%qW^ZmQ@Rs5F7>)QPeVDQ|@vrmNY9@@V z)dD87R)e&HuFCQ_dTXr!q}KWmPQog7^hMGiRBEji!I>qhun4=*`t7qayuo0E|7h=M z0CUCO8dzM2y48BpIC8{@dN#E*n61QIfw@)&u8i-M4CS{mm)5Ni4tSs z1*AbPlNzf+X8p{X^wRMH6#s!LPHD@T(N{$|b*lj^0 zN$aS3ItQwofWc)42r%`9$K240*}4KQ)~*wO%(z&;T8vv`@D9JgfTiJ3S4n%@X3ta4 z-0Aa#&mn8_FZAYGi!d@*JDAK|EiwdjRo2GQn`<2)HP^a0fw}xQYv+)mpi*WOdRE@Dr7 zOLvjA+*ZV|s^wsp{<Q@1RSr7tQbuer z4lu2s#5}gSW843BI>5BWQ2F_R;g0TC!y&jOEz!?91otE?J~rlVDtq%`eYDw!Fp4(& zf=R8>Sd#rfS7m=3y)_O1QfpLkqS0oc{X%w04g{52;~;Rl?9tC66(I*&AI#*~9*2

V@Kp+#-iEPhgS5fGY#j52_r_TsmgdqLGHHq&M8L{#Q_I@`%x~B8s%sr zXOv?I0;3!YY(@?IzwvSI0tTDo-38V>ClDWH{=#5#sxSG#&ikMoc4AJ zC!hk0;1(K0Ow_R7V766~zHG?w%3na_l)sQ5Q2rv9$PjU{OSCFqLj0tv zY`%!I8fq4nP}7w0U)mU9{{4&BQmRt@d^-{B30k$ev;^@6}6_9%AYMg|%%-FkSG^o@|W58)6gx}1n z=?eDp$1+!Jr)z=>`5mpA*UGib6*%fT;L5;h-#dKcGxu1daMblKmvz()#NC{uZX`+^ z)hIW)JnE>MiJYTuAqX6Gt4pM#ZgYv&QMVJnMMtrcfLB~aV=ZMrtzv5$Tg6z5*!e;X z0?<;rJvCa+VSmgus%^W5uL){5+ja4a-!m(C`{UP*#isPn${qYWy~pk(j68N1n6$oX zZj`%0SLGfYy~*wcQj^_>lQ3Y7ZIf}JQj?7brwe9GQwD>S?f6F9P=v}MazE2zYdt`U z3uc10%urm=t`|JW)WA{?0apgRmTVaLFmtgb(^%@K4dO>gbaS4XK$Ns(@IBN=T_!cv zV?@qWj}ru@dIC5b)_#(iu5@N%7l-u^TdU{?u`h+VV_Lv)5RF#iJQGcTY_`FpwTphc zjl)f*I%#!&>EM=HNomK5_v~x7ALKm|()b^HPy7NpQ5+wo30sHDQ+!@;il+%9Q#=DE z4LmLKEZC|%hokqz^FZo}7jO~`}L5*dcxb%{0%y+^zVN5MJ#-QPzcsD;?V75uEyZ2XDGy-}sH6Aqhc976Ag zbsOIv<0ZEGGj&|qG>-AYv+c@j9uv-NEoft#X&fP%u}@_xL3j&KVFJ`@Xqb=S{nM-H zE%Pqx|H)I!`}~LUmi>S*vh0Uok~r4fuSGrrU6qe<^q&0$NIm;0PQo}gW+(X!RO;E! z!TC;B+$m!$`-17QO}`{HA3S7ZSzTXNmVLwK`d7>nnDcAkd>H&2=9&dJYOF8yTv=yp z&3xUrOz7rZ`W;ac*jk$8dzVRV`U8=(>5l|~O@9K8hr@qnj&b+}ts=F9l=aIHRZL)` zHDQY2&1h#@Yed`1a-prst*^1Ghily);f?q2OCxIjHnz0;PTmjUQxL!K5Aojll`!(o zZ(x#ehcz$1gRRORIC|^+38dEf3n#G7>vOJZyKE23-=I?KbXgp~r8|VRl`tz}XZFz= zXxA^BjmNS%9GM$$~_2XlX6C%0e!aI%Q#ZJtuz`Aqbq(6F3eny_n+<=LYV>p$#VJ zhiXvUNT}6JU||~DE>ztH8PrC;PAr<(N`|Yo;4j{yp0PtKR^t{Mcp98(I{I~OToqe1 z=m@}+{e4-K{|epluVBhgZ^9@pEe0lyOAWF(*s3gnqc`7@Kx)1^oM>EH#m1$xqz|an zd`p2V7N=wBkdI75v5aVGrpAt2hLkQlLaR&Tl6~>2&&ZN2%f!HA%K=vg{Fuz>_GPX# zCiNp(c&+5J_?XlnWtT;rR!`(&QUgKYv_{}s>B8H}BcaaTEEyNeljBFgU8N&d; z_Bqc@g$4L~5c9#`FrEVAiZlfjfFJL-qI$h6ceont)!c6~H?!#Ht@<{EB z4dex}K~^Hwr^Y?Z>{+RgMmNm6;g)6Y-+-MHc}ZZKHsIVgt1y>#>6*+O$EqX>%M2n) z!c$X&R9q%C&1yu>G=m8O)2t30ho?1|qs6UFDWmFPc;b2!geMJ1XkIh6uw{<_nprNS z*$!6;sx$DQ*=_p$jXGV1m=+Eaus%YSr_qd+oUxjjr)7f&O?5ZFipXb{xSB2;)(qHY zQ^(IX8IMMr+jv`yfTJ7X*Ra@bhzYg&{lQ`k%|1Amuk4IE*Zj5ku+p0R7wleywzUYO z(AEwnbu+%5J_Ix-<#6}E&29!|5g7_9b@O`Qim@$xXIgzcl1(WWvPO(5 z$@-ZJ3!em!fJNHWWjrK0>yaWJBf*oVyWlWnlW<(T$l)o!i`KX6Hy0 z&fSzKaW2M8hq+AZ+~GvdxtkFL&fVN4vTM18OSEBaOX7F%*Qgb>idDC@X|IxfJExa^ zbO7#Sk~w#Ct#Z95<}+x9pK%{F$T z>;Nh?+m7ILJc9$N!H5K7XXd$`m^t>`&ZIQ-8qfJL4IN$0t+xv=2&^{}xH2&PUf{0G z6+L}hSZ`0@IB4$09Ca>erhe0LifCkl;Sa<3wv}Iu zC)fzrNlzLxx|M*oBJHdv|BGju9N%{QY}U7X^Ixoc;oE%(Bj4@|CT;&U%YI-nm58Ia z?Eye)+bT{X+g>6Ef=X?B5IFNLBdcdD*JHKb$8<^$W?JmhLyGLk>T+`y3bpv5ObZNo z7;t6a)F*tm?OlWWtADt>>Ts8!(g}#Ls^e$KjLlTU@YyE`e>m6WvT~nCT*!SsQ7U&? zE^wJt?hA>W+!ql9a$oEcDfcBV(aL=(@zIom7uF~6^`BkR|Ba7M&m>s zV5g18j*>B;QU{C$$M#8$#(i2ru~B{vv&W9Omedk8?f^riO(WjoI_3{7ay@X~uD^k~ zw2EpeYwPGn5`{-@B8olI+FF*IT_*L&Ekw>Ew-N*%xeYk3qHbr7Mps3#>IsG&(2{lY zG+{eQr*SkJ;<+~JQ8m#Uump$SJ8m}F_L8mV(1}2G$DRTIgZ?Bg_RS`9$6_0f4j--4 zPe%Xo(QJ>jX8Yt>WRnJDX0Eerf04!_34Wx zoNdF*Vq;$XAae%BeF(TR?!AX!T;!^48j7myVdgT%ZR*MQrXC?t7hh>>A0u+UeVich?GwQ18ki^D1suS5%3WaN-P4L`YcpeGOvxD=HowhFY}gvApAMay zX(F|tI?=Y7tS^mO(T~-6B0K6~)2xH7Yjl8Bms7aac^IT%Et$?^!HFh^nVN!wbnkq2oD3K4*s5e_xQ|y5qKXE&{CMRWvK|b7fec zVcMGilDPu!eFcoqlC0Hbx!!fZ7vB4txs3N()O+8MD7^PAQCfGm$agN2n(uoeXTBc@ z0`vU{oGr8XiJ2b6U!e6cw1+B;4K`G8Vvs}rG>IzIL5M*HJ4tS&gPyJ;*Z66NFz784 zdx{{+z=V*>AHdHPli+D!)Z6hVU>*OWVV_~9O4hc+aQT^kOn3Y-X3*jn!l*L+6-?SE zXqDf{ZO!6u<#=OkZ!y98h;L5-m1~;$q$W`AI@nn9N%Nk<=;=&jU5~apykcC_(HO9h3 z&KQdj^zg>$=@My-UM|rZV^QMblUNM$Fa*L1M+ne;M89wfIHcxX|9h+S<}-V%EJhev zWpOa6RZ6l1*s3gvqqj;OkXoe=PQtopY;Rc#RBDx_!P)m$=#yf=(hUolL8N7vD|XAW z#f4$a!5!SNyjDQVF;@?7n7+W3arY-LRP}SYtYJ#Tg<;A>sbLzW-epq5G!Qw%G!pdi zhG}w%G)%Kgw1#ORej1sR;Utw4?Ns-+vT2pLHi~U7*yc0(W}K9s(nC-F5lri9(vQFy zCIr2b9j0pK-{noUJYi(26~H7dnijmP_8-vr&Ki#1RVxCis|Mf%uJZl6yJaA#)Kx2i zvkhxDOYF)%HcrsOL0``jvNH33nA+NYLX zwFpLon2U{@rV^H1p?_B)QJ8BrqNribe}lm;k9uo$BIm6&2zuatww78GIBi(3*DB5I21@txE@X-6Fw;GgGx=f0l311)^VnV9Gh@M=8H|p@v)7t@wb{VX<4&>*Br zz)57p39=)o)QCHQOO1$F2rtG)+?n}eBkn@V8BsfS#*J8cQBI9Gl6iZ0Bkl@}wS>uw zKf5uPF=D;O!`(>~MiipNh-De&@~9E_AaX|Blc0w;;$Fb15%+c%(1`oE3#<|MRh%2q z3TI=D4M4_^cI-eSU=%meP19Z^7hj_@$v!1wQ760ji4oy1`${R_2UCx$FKa)L8)0KS z#8c?cmy8cJYIgc5WIz0Sq;0*j><_jo2jJ*^SOroa9*C2ut+%Bd1S<96!Qiypr&TL^ z1v*cNvC{~bz7w>>qjJT~e$w}{>V_Ha0isXe^o4hOCb_~xvsKiYb&ql?gY zgv(`HcSjNzhC7NVZQZrV(JqrZ?ieEHxMK-=c*h;*64|;t-X&Uhoj`mzGVA;gduVp8 zk6G%&5eD=T8O7La;>sSoN(U?=P#ddaJ3l~waK${x;xf&kx6sUR8Z*01Ky+@wjNI#8 zZn(715_&hnc!-OkFRYQ@y}v< zEdJS~YKxCAd7Q)49$x%&fhz-^Pnuunxm;HK^N9=bFCa?AZoHTw@hL2^{#R`pT>)S1z{xqm0;2;sJ=arEks22%CM z;6&A*1}jAO{_U*M@5@+FsruJ|Q+2Tk5Vuyr^vJc2@>>XAM~cCJD9blCZefih*OTqx zg}wne-w^#q=F)nlzCr7in@AK&-%P|ej@??&b}9A4_64`NTq^agL{9442zq#_ZwJms z5$|B8jge{7p=b}V>B)IAG+<0q6{WS-PJx_9V=Q1^#Cf}u05U^(G?zgITLs?9pYXcf zMHuOJH<(PfX1ND+Rqn;n>vkWI>NXB1k#2X(cu=Ws_k+_HMyzhdl*9u}h>Pokq`0V; z>iK4Z5%D1=_CTWTEIkZd8K16kdxW|8@88r=mI(xffR7So0=CLyE|&`UxVr&Yjz2-r z19!7sr6++i0iR-~6Okzuu;>m=-qwMkv~~oUYKFZS(~aiK2e1n#uxxP;$9LsvK9^Vh z8Nx{QXTfBux5{&%tMWXKUiBA%RP`5e5~=>HOazsx{t~#Nw3-b+Uv`u)t*?;srPXfu z`6}5SUg+0=E92)%R&T$~Tt?^yd4r%3`c0x#XuS9Tmdm9=zwK^7q2D3s;e~z|I1~Ck zX3{M^^OC9?{dKaSvVC8(oW9Q=@QQsv7%BE4n6z?jkdMIPoiZG~V4nb~V4vb7tX#+3 zE}wx)1^XPF)|6Z~*oRB7Fa?8Vl37#r1ykd)`ekxK30b*9)m-ouQ+s#~zXq<1`(ulX zKH0v0g3884Sl_r@jg7hp>s#VN#_x!dsud#}-@8nz;txbF-u_6?!>jm{OJvphvrDv9 z>o3G7BQt^}hE>r#Hv{fOSag%D$Kr|$wCA;r#RXZ!DBl`@wRPoIK%-^v?}NJ*H$n?M zxTwy|jpKeAy>U&k#7BPRlX)BbMi|-PcQ9$eZRjU|fUe4)IC>xa1*AUs8z)i0eNehA znfqWWa7DNdw|a{bS-tnrKlgmEd*Of{%wAPw#HPQal3zf76ER8G*A>JQFh=Pu-JJ3^}y`|F5Do-yc%>(dL$! z`BYx>SqLM|X9bfM{YIG$Y*l8*(MvuDkV-x$PD07Y441h;rA2>kaJoS^-fc&h$_41oS4pDGC87)~IT=OmfzY@Ac5j%)2_wTS z0VWL+jrFo5=&ID==pEAsNFB2jPQuD!>>;u=sMImbfHR9Y>)M96p>Pb1o9)JAStiGJ zS&o!`w#aNTHmt@LV||$?FibySycw7nro>#+K(nT%23@9CCQ&%1o+ynPO^wpva;aq+ ziJWDc2m;GAyF^;1#U)zHv=TpKEJLdpd!VV#bwj0~YFOJ*0YX>G`$J)HN zt)jrH4-gi+0Bc1>W!0mI^=uSHQB=fkQLwvZ?d}%4yE}j1@0po9Gk4w!{{HZKzh}<9 z@3NmWbEj`s^=sEYKHndSQl51-yE6ITheOQb1-^IuyeVzYn@LEx7;%ce87)L~nIfdy zM2ubA|KXuzUH-#-Jgr9=#nbw*O|YnR%vB&98Vk3v?MU-(wH_Z4r{F)ingA1FtS^5uxw11Ae=UVP8TMZ z=vBpV!Zw>6Td*KAor+0I!fG5b_ zFy*1|1Dp9&h%W=;0fv$Eh>L;_n!4jpXop*~A|M=T-Iq>51ndW0+VI!h=VImX^8wCGhV9W1JQWtO#YA@hJRK*ep>KBEST{@EF69O z908^Ab0ki*(6D`kQ@hJiFy-UtXxO$}34bwNH;fqSaszuV3pj?}aV#B6)5`7!c#X9B z*&uBuCOK|p#}MQ=x(87;20Goy{&;#-;igRmvsKCoG-OdVmMT}S6bkinqVv+YI*H1~ z)i{bEu1F z*CZD?FAd&{sa)`0LJF;c6G?b}6t4Y-Rxy!*T_m(?`l=$|zE< zge8fT2Du96vP{9zN6J(vjg+f#5=F|Dat%!RNVyiaT%-VZt&0&(XH($U(Jzja>uE|O zrGwY%>H)p*26_dNawBvaDL2upDu|SuX~-hw7OEstn&npKrIB(Qm5Y?yDS}A3!zD6O z?sSPZQtqOz94R4K{#g(l4ye-gWdLF*kCrdSeCUNYw4eH7oY>YJgN?DHMOFTGayOsT zN6I~vQKZ}pOA;xKav#iPnTDf}lujs(l>2d_k>Zx#jgkjo%16qBuxX_lYdbph^91-B z7g-nOAtuGa@o?VTh6R?BJi??PA|8cKSJXX5uPVOy8tSy{c$|hTBA%d1BBDi}bY2<} zPf@vuc$y-Jh-X|PBjQO)x!w&9JZFCo6zkd^gChSJ1r#A{0oJ4qkF;9+t; zay#(P#MQcgI70x}*%)GfZc(1+GkDrBP)4+0ge9TfA}_&QmX~q#v|oWz+OOgyqCKm; z22-B)>#&(Q^4Przn>lfrnIc42eX&=YyutJs`E;5V$DgcKFM(ri8=BUv@U)>_Rg^dB z5P*LRI&D#{phU?yQUJCwwDhK`piU9nFE|K6ra)}oF$JA#M zyta+mBCldoaQhr~Cq-~ig-4hDB_Yu{FBTM zQy%dgu;p^SL2PsVoJ@~7&qb5zJt^1>kIc^ej{gc$#pHl>HFUZ$WDUKl$hz6=$2~U< z8S8ndV%GQv&+B}Y^?X#0_52h8>jj`Iv0l(!MAi$ri!AGfsZXEzlE4i8C@)_UAk%hH z61cA?<@3Jj!fDm}_7RQN}{H`6m zu`B^oJ_eSAO=BQz0JtbX)qg3b$1%{o(j_LhYaMWDrUyaL1G>~XU+-0(TKcld3XKE0 z%g{@MprsC-^L=Gm8nPf*jw)a9z}(c)f(dUxWk zTxCORAtub;_^cJ$^);?m$QJ|)1Xt@{mA~hU@>u{0PU`ZGkWp$X# zvIdTx@tRP|crBcSZ}+JqWo?-9jMsrpGgM!}o87{dp0vhf3z4HmDn(iYo^{Y^rmm+~ z$g_cljAtWN!m~x1oR{)!rgA)6C<31Axr)RPPo00}UfMcg)_p{Q^P@CF z;VoUk%k`zx4R8N~^x{1|vkfS#LuMPok}#{2-Y}P?500K$UnphPij#=hMA-h6RC#7twKvFl zDYG^z$80b~z-%*@NM@V6M9XXoYGtMqMrX*Zf>*+9OWxBn+lsO}WVSUd3A1|H2IjJC zi=$_@9h5TL9w!mAyJZKM^2~<7rY}bj3$|=?OolQ!hBvJ86)ubJ$mD?HPSEMpn4Rer zQrv}xjN-0T3B`i6J1?cU87M@u{Qp~Sj%-9}RA6NkhB6)tn{V|Uqs z^a>3gcol4Ov_qOqIk%Q4Vol8W)FBxwtrtB8ZE_T_WS6 z>=JEU96@~`OJuw)QJ-yFLt9p+0~vfwSh&L04N#d0>3q>P>&8_TE)A~KfU9(=9LZnz z;zSEXJJs3vGjkLD?gau%LB9~szRd=kuMIU7eG8RtN0WSol=EfZ{Qd&&uN9!&YjI3G6M zM_A_~aL&Ryt#Sd4aWGs+Q#BxBh81Tl43>*%3E*E0UFxiTY1|h>)*6&$`Ag`)Xai)9E7#%8V4kDc**c5FE7ZPGux zj_+}p54ni{Lx$o1(6a?lenTW|Ds4G+eSzG^U-HrUjtG#9pH}X^gYiTbZ8) zMZdRc3{brT4OBZF^=6&)$jZQX>19wg$$JzTs`sf9RE_e1^HNkFx*HJHM-&07kD+5! zpU`8+8KTnsh3wS}`D>=LSKu$7@^?J8&nT-yY@frDU~7>tU@ps-IC^YfK`FMcaS~y> zQoeyHkL_F70&JK(2OQtg7-RcB^HXfb?*|$KY(GMmI+vOKo)>o*Gd9Ea6TJ+!X8D;S z!}beRf~`q@bzX|?H+KVK`<)`d_6Kx~?N54ii489*1vRs!Vwe8%7k|fN`l#P$y? zd2DE8&$1e^&5EPPHXD>;n;j<+wufa7nDW@>ge`cNalwVmMOzIXz?X_unVX|lkZRfj zXf@F3ocy`zrIk(nW?5gEhlUK?yi^IeW|_}*v*;!uie37i1c^lk%8o8jV;ED2K{ z)l#qps5moCy3-pJ)OztqykO+OPczHnPe4aKNitf8feHMHB#mC64p31IWNW3Oy%-Y3q^owUFbM3tw#@C zQ0nImatzBfk-wa{g!{A*6u=TG>^=tUW2stsUCpw-ESFqfq_jviVc zD23J+ClRzB(h5@^+D5RMPEx*Go7mt<$4yg`bz%nM<4VCq2q7C?k|x!jhn@m#ttf%hou0l-ocl z%58BHp&tcA&K;LOFyci;`PDaGpOWVaiY@1Sp3=rzm%% zSBP>a8ZwkSQza;ivWxRll)F+nlDKe^A z`C(>jA`|NN46quEx!mYx8n9g_yYnYK%{?e1ntQ^M&}@(qFqdU796imwp_Jx6I8mDJ zd(cy4Busgl9kA(QEv^YdHy1uR!BzVE(p(ea+>fROVd5v{BcA&+CEzL0rD^}XefQ{= zANuHDpAEMI=oRujkcN!sD5`{KogC!6l;>zF$Max{faf95an?MP9tSv>C`*F{KmIc% zgvd2uV2s|;j5M&@{OK;t`bjw5#N|%pQM)7%7iaTPn}^9^{Be)w5z{d=aj6~~fy|$#tUaFQ z0O$$OrOpK=zkc7qSNmr`$I>eVdLj)O(37YVp!G7&c`4A7sT|N#C<374p<|$@(xaU~ z!3;IaYznA(H=2oiwkKP2sQCz={qh2$Akdl z`OxW*)dlnlAznyB2Js@Q1Y(0+?7S4>B~%XLr4#|g%UmLdCN6i0wg!3y^&XZd8KU6S z!q$)^VI4TGSlA*22r(&Z`8o|=q)}NrnLp^6UP&1-y$Y6uX|qg$xhzw0^h~dYQl{77 zB+8WA%C#`%nO+B*$!(?*#QjjPqs{d+*F;cnpowLzT^@rEpRqCo=Np+4P`wGd)LETd z6LT}YLaMjWkWsyrDxumaw>dARdOMY)dIv>7^-kzGYu-f4g(X4RBKN_JNx3+Bl$}tD z@_w8|+4B&20H!?32Vv7Wd4t<)fw30*8F89F2|dK5nh58^H0h2;Sew#HXVz>T(npvY zAbk`XvpV{X{<>|Im$TCLF?xkaAEzNh`UF*iv`L|&@s|y=`qUY zPPRe70CTam;5E~0lreDdgD09(bzmrCq-4#v#d)Q9e=0rm`s90*fVGbR?_dE9mk|wS8f*wfrH~k9*S&9h($nMak z&Utg27%fe&5M&P;GLXwqB_Ip3tn*Tk%TYOy%Toj(S8$1}_O9p>ZMAnL>f!dqIzy)Y z(MjOD7`j5LUSM5;=*s*#PjnT^i0G=YBt+|_7Ur_7hNCCCI+PM!11C{UX_~AFQ=aHr zu<7a)Kdsp;iDa9U^C{M*LkxEvnz->uc1_$F$jf@tC*WLwPMbGHdWD?pXvjF%Qze{> z(%`(5b0d}G+(Z#@ZiY@eky_|9#{2_V*~N^jV9Bm&xv~Pq+LBgp&0+-<&F-V|D0A6Y z!PjlLkd4h(?q1?*$h?IzHw{>A4*9J^?M-A|K91+V9%aOTeOMCy*d_#KY!iZ`=f5G8 z^6!n4Fh@_BB7I=W^Y06r=IAyyWZ8^m|57WhG4zdSVvc6HjGhSXG1NUlm=KWO1iIAO z?^(JLg~kA{qNf2xFZoeXeZN_tZ&M|vPdKza~#Tr{`QV+1jCF@d&u zxKJ8M?ZUOFf6$^vDyp~uYp-z|kC~F%clcVcc6yHM`5nxk^e8u@j8JY4OM(*Hd%%qC zJ#h3Ww}MiXTjL}|IdxUp2Btj9ZDI2$4a|aUM{A67dzviDsHboTCIl#lK$kif&P|3M zO0N**Fd8zHJ5nVm>t!eBr6_l%awvD92vF_{9iwcgM}~5SMRKm&mnb#w<-2|Vtt(5MDUZ3pswcpNQ#vI9n zm}Un}+%Zmv(#%Xq-B5w`rhSz#9?q8O$)g&Knabo(du?l#xN-eg5#k}olAeY zRpTAM-I$foC%79lH}o_?Q0Up%tOKlLnUqDsiBw4>G|5TMOCw<%m5YRvDS}8i#U(Nl z#=Art38zw@Tsac7{Iv-2l6C7$85MTYs+ke1B?R#0KNDvw*}RRum+ead_hH8~ql;Ij XF;rngh(T+ literal 0 HcmV?d00001 diff --git a/neutromeratio/analysis.py b/neutromeratio/analysis.py index a9efc2ab..46b97e73 100644 --- a/neutromeratio/analysis.py +++ b/neutromeratio/analysis.py @@ -30,7 +30,7 @@ ) from .tautomers import Tautomer from .parameter_gradients import FreeEnergyCalculator -from .utils import generate_tautomer_class_stereobond_aware +from .utils import generate_tautomer_class_stereobond_aware, generate_new_tautomer_pair from .ani import ( ANI1_force_and_energy, AlchemicalANI2x, @@ -557,6 +557,105 @@ def setup_alchemical_system_and_energy_function( return energy_function, tautomer, flipped +def setup_new_alchemical_system_and_energy_function( + env: str, + ANImodel: ANI, + checkpoint_file: str = "", + base_path: str = None, + diameter: int = -1, +): + + import os + + if not ( + issubclass(ANImodel, (AlchemicalANI2x, AlchemicalANI1ccx, AlchemicalANI1x)) + ): + raise RuntimeError("Only Alchemical ANI objects allowed! Aborting.") + + ####################### + + #################### + # Set up the system, set the restraints + tautomers = generate_new_tautomer_pair(t1_smiles, t2_smiles) + tautomer.perform_tautomer_transformation() + + # if base_path is defined write out the topology + if base_path: + base_path = os.path.abspath(base_path) + logger.debug(base_path) + if not os.path.exists(base_path): + os.makedirs(base_path) + + if env == "droplet": + if diameter == -1: + raise RuntimeError("Droplet is not specified. Aborting.") + # for droplet topology is written in every case + m = tautomer.add_droplet( + tautomer.hybrid_topology, + tautomer.get_hybrid_coordinates(), + diameter=diameter * unit.angstrom, + restrain_hydrogen_bonds=True, + restrain_hydrogen_angles=False, + top_file=f"{base_path}/{name}_in_droplet.pdb", + ) + else: + if base_path: + # for vacuum only if base_path is defined + pdb_filepath = f"{base_path}/{name}.pdb" + try: + traj = md.load(pdb_filepath) + except OSError: + coordinates = tautomer.get_hybrid_coordinates() + traj = md.Trajectory( + coordinates.value_in_unit(unit.nanometer), tautomer.hybrid_topology + ) + traj.save_pdb(pdb_filepath) + tautomer.set_hybrid_coordinates(traj.xyz[0] * unit.nanometer) + + # define the alchemical atoms + alchemical_atoms = [ + tautomer.hybrid_hydrogen_idx_at_lambda_1, + tautomer.hybrid_hydrogen_idx_at_lambda_0, + ] + + model = ANImodel(alchemical_atoms=alchemical_atoms).to(device) + # if specified, load nn parameters + if checkpoint_file: + logger.debug("Loading nn parameters ...") + model.load_nn_parameters(checkpoint_file) + + # setup energy function + if env == "vacuum": + energy_function = ANI1_force_and_energy( + model=model, + atoms=tautomer.hybrid_atoms, + mol=None, + ) + else: + energy_function = ANI1_force_and_energy( + model=model, + atoms=tautomer.ligand_in_water_atoms, + mol=None, + ) + + # add restraints + for r in tautomer.ligand_restraints: + energy_function.add_restraint_to_lambda_protocol(r) + for r in tautomer.hybrid_ligand_restraints: + energy_function.add_restraint_to_lambda_protocol(r) + + if env == "droplet": + tautomer.add_COM_for_hybrid_ligand( + np.array([diameter / 2, diameter / 2, diameter / 2]) * unit.angstrom + ) + for r in tautomer.solvent_restraints: + energy_function.add_restraint_to_lambda_protocol(r) + for r in tautomer.com_restraints: + energy_function.add_restraint_to_lambda_protocol(r) + + return energy_function, tautomer, flipped + + def _error(x: np.ndarray, y: np.ndarray): """ Simple error """ return x - y diff --git a/neutromeratio/parameter_gradients.py b/neutromeratio/parameter_gradients.py index ad9392af..a74c3fb1 100644 --- a/neutromeratio/parameter_gradients.py +++ b/neutromeratio/parameter_gradients.py @@ -37,7 +37,6 @@ def __init__( bulk_energy_calculation: bool, potential_energy_trajs: list, lambdas: list, - n_atoms: int, max_snapshots_per_window: int = 200, pickle_path: str = "", ): @@ -1048,14 +1047,18 @@ def parse_lambda_from_dcd_filename(dcd_filename): assert len(lambdas) == len(energies) assert len(lambdas) == len(md_trajs) + if env == "vacuum": + pickle_path = f"{data_path}/{name}/{name}_{ANImodel.name}_{max_snapshots_per_window}_{len(tautomer.hybrid_atoms)}_atoms.pickle" + else: + pickle_path = f"{data_path}/{name}/{name}_{ANImodel.name}_{max_snapshots_per_window}_{diameter}A_{len(tautomer.ligand_in_water_atoms)}_atoms.pickle" + # calculate free energy in kT fec = FreeEnergyCalculator( ani_model=energy_function, md_trajs=md_trajs, potential_energy_trajs=energies, lambdas=lambdas, - pickle_path=f"{data_path}/{name}/{name}_{ANImodel.name}_{max_snapshots_per_window}.pickle", - n_atoms=len(tautomer.hybrid_atoms), + pickle_path=pickle_path, bulk_energy_calculation=bulk_energy_calculation, max_snapshots_per_window=max_snapshots_per_window, ) diff --git a/neutromeratio/plotting.py b/neutromeratio/plotting.py index 2b003df4..f10551c8 100644 --- a/neutromeratio/plotting.py +++ b/neutromeratio/plotting.py @@ -55,11 +55,13 @@ def plot_correlation_analysis( Y, yerr=error, mfc="blue", - ms=3, + mec="blue", + ms=4, fmt="o", capthick=2, + capsize=2, alpha=0.6, - ecolor="r", + ecolor="red", ) else: diff --git a/neutromeratio/tests/test_neutromeratio.py b/neutromeratio/tests/test_neutromeratio.py index d2f0e792..e00f75b5 100644 --- a/neutromeratio/tests/test_neutromeratio.py +++ b/neutromeratio/tests/test_neutromeratio.py @@ -1273,6 +1273,57 @@ def test_setup_mbar(): ) assert np.isclose(-15.506057798507124, fec._end_state_free_energy_difference[0]) + for testdir in [ + f"data/test_data/droplet/{name}", + f"data/test_data/vacuum/{name}", + ]: + # remove the pickle files + for item in os.listdir(testdir): + if item.endswith(".pickle"): + os.remove(os.path.join(testdir, item)) + + +def test_setup_mbar_test_pickle_files(): + # test the setup mbar function, write out the pickle file and test that everything works + from ..parameter_gradients import setup_mbar + from ..ani import AlchemicalANI2x, AlchemicalANI1x, AlchemicalANI1ccx + + test = os.listdir("data/test_data/vacuum") + + name = "molDWRow_298" + # remove the pickle files + for testdir in [ + f"data/test_data/droplet/{name}", + f"data/test_data/vacuum/{name}", + ]: + # remove the pickle files + for item in os.listdir(testdir): + if item.endswith(".pickle"): + os.remove(os.path.join(testdir, item)) + + # vacuum + fec = setup_mbar( + name, + env="vacuum", + data_path="data/test_data/vacuum", + ANImodel=AlchemicalANI1ccx, + bulk_energy_calculation=False, + max_snapshots_per_window=20, + ) + assert np.isclose(-3.2194223855155357, fec._end_state_free_energy_difference[0]) + + # vacuum + fec = setup_mbar( + name, + env="vacuum", + data_path="data/test_data/vacuum", + ANImodel=AlchemicalANI1ccx, + bulk_energy_calculation=False, + max_snapshots_per_window=20, + ) + assert np.isclose(-3.2194223855155357, fec._end_state_free_energy_difference[0]) + + def test_change_stereobond(): from ..utils import change_only_stereobond, get_nr_of_stereobonds @@ -1578,9 +1629,7 @@ def test_orca_input_generation(): qm_results = pickle.load(open("data/results/QM/qm_results.pickle", "rb")) mol = qm_results[name][t1_smiles]["vac"][0] - orca_input = qmorca.generate_orca_script_for_solvation_free_energy( - mol, 0 - ) + orca_input = qmorca.generate_orca_script_for_solvation_free_energy(mol, 0) print(orca_input) @@ -1785,7 +1834,6 @@ def test_parameter_gradient(): potential_energy_trajs=potential_energy_trajs, lambdas=lambdas, bulk_energy_calculation=False, - n_atoms=len(tautomer.hybrid_atoms), max_snapshots_per_window=10, ) diff --git a/neutromeratio/utils.py b/neutromeratio/utils.py index 389312e7..cd8bab92 100644 --- a/neutromeratio/utils.py +++ b/neutromeratio/utils.py @@ -129,6 +129,19 @@ def get_stereotag_of_stereobonds(smiles: str) -> int: return stereo_tag +def generate_new_tautomer_pair(name: str, t1_smiles: str, t2_smiles: str): + + from neutromeratio import Tautomer + + return Tautomer( + name=name, + initial_state_mol=generate_rdkit_mol(t1_smiles), + final_state_mol=generate_rdkit_mol(t2_smiles), + nr_of_conformations=1, + enforceChirality=True, + ) + + def generate_tautomer_class_stereobond_aware( name: str, t1_smiles: str, diff --git a/scripts/Generate_equilibrium_sampling.py b/scripts/Generate_equilibrium_samples.py similarity index 100% rename from scripts/Generate_equilibrium_sampling.py rename to scripts/Generate_equilibrium_samples.py diff --git a/scripts/Generate_equilibrium_samples_for_new_tautomer_pairs.py b/scripts/Generate_equilibrium_samples_for_new_tautomer_pairs.py new file mode 100644 index 00000000..19f4e51e --- /dev/null +++ b/scripts/Generate_equilibrium_samples_for_new_tautomer_pairs.py @@ -0,0 +1,98 @@ +from simtk import unit +import numpy as np +import mdtraj as md +import neutromeratio +from neutromeratio.constants import _get_names, device, platform, kT +import sys +import os +import torch +from neutromeratio.analysis import setup_alchemical_system_and_energy_function +from neutromeratio.ani import AlchemicalANI2x + +n_steps = int(sys.argv[2]) +base_path = str(sys.argv[3]) # where to write the results +env = str(sys.argv[4]) +model_name = str(sys.argv[5]) + + +if model_name == 'ANI2x': + model = neutromeratio.ani.AlchemicalANI2x + print(f'Using {model_name}.') +elif model_name == 'ANI1ccx': + model = neutromeratio.ani.AlchemicalANI1ccx + print(f'Using {model_name}.') +elif model_name == 'ANI1x': + model = neutromeratio.ani.AlchemicalANI1x + print(f'Using {model_name}.') +else: + raise RuntimeError(f'Unknown model name: {model_name}') + + + +if not (env in ('droplet','vacuum')): + raise RuntimeError('Env must be `droplet` or `vacuum`. Aborting.') +# diameter +if env == 'droplet': + diameter_in_angstrom = int(sys.argv[6]) + if not diameter_in_angstrom or diameter_in_angstrom < 1: + raise RuntimeError('Diameter must be above 1 Angstrom') +else: + diameter_in_angstrom = -1 + +protocol = [] +names = neutromeratio.parameter_gradients._get_names() +for name in names: + for lamb in np.linspace(0, 1, 11): + protocol.append((name, np.round(lamb, 2))) + +name, lambda_value = protocol[idx-1] +print(name) +print(lambda_value) + +base_path = f"{base_path}/{name}" +os.makedirs(base_path, exist_ok=True) + +energy_function, tautomer, flipped = setup_alchemical_system_and_energy_function( + name=name, + ANImodel=model, + env=env, + diameter=diameter_in_angstrom, + base_path=base_path) + +energy_and_force = lambda x: energy_function.calculate_force(x, lambda_value) + +if env == 'droplet': + x0 = tautomer.get_ligand_in_water_coordinates() + langevin = neutromeratio.LangevinDynamics(atoms=tautomer.ligand_in_water_atoms, + energy_and_force=energy_and_force) +elif env == 'vacuum': + x0 = tautomer.get_hybrid_coordinates() + langevin = neutromeratio.LangevinDynamics(atoms=tautomer.hybrid_atoms, + energy_and_force=energy_and_force) +else: + raise RuntimeError() + +torch.set_num_threads(1) + + +x0, e_history = energy_function.minimize(x0, maxiter=5000, lambda_value=lambda_value) +equilibrium_samples, energies, restraint_contribution = langevin.run_dynamics(x0, + n_steps=n_steps, + stepsize=0.5*unit.femtosecond, + progress_bar=True) + +# save equilibrium energy values +for global_list, poperty_name in zip([energies, restraint_contribution], ['energy', 'restraint_energy_contribution']): + f = open(f"{base_path}/{name}_lambda_{lambda_value:0.4f}_{poperty_name}_in_{env}.csv", 'w+') + for e in global_list[::20]: + e_unitless = e / kT + f.write('{}\n'.format(e_unitless)) + f.close() + +equilibrium_samples = [x[0].value_in_unit(unit.nanometer) for x in equilibrium_samples] +if env == 'vacuum': + ani_traj = md.Trajectory(equilibrium_samples[::20], tautomer.hybrid_topology) +else: + ani_traj = md.Trajectory(equilibrium_samples[::20], tautomer.ligand_in_water_topology) + +ani_traj.save(f"{base_path}/{name}_lambda_{lambda_value:0.4f}_in_{env}.dcd", force_overwrite=True) diff --git a/scripts/generate_equilibrium_samples_in_vacuum.sh b/scripts/generate_equilibrium_samples_in_vacuum.sh new file mode 100644 index 00000000..07cb3ca7 --- /dev/null +++ b/scripts/generate_equilibrium_samples_in_vacuum.sh @@ -0,0 +1,30 @@ +######################################## +# idx has to be provided +idx=${1} +######################################## +# idx is an integer between 0 and 5000 +# idx is used to index in this code: + +#names = neutromeratio.parameter_gradients._get_names() +#for name in names: +# for lamb in np.linspace(0, 1, 11): +# protocol.append((name, np.round(lamb, 2))) +#name, lambda_value = protocol[idx-1] +######################################## + +n_steps=400000 # nr of steps +env='vacuum' +potential_name='ANI1ccx' # which potential do you want to use? (ANI1ccx, ANI1x, ANI2x) +echo 'Using potential ' ${potential_name} + +hostname +echo ${idx} +echo ${n_steps} +echo ${env} + +#conda activate conda.env # fire up your conda environment + +base_path="./" # where do you want to save the ouput files +mkdir -p ${base_path} + +python Generate_equilibrium_sampling.py ${idx} ${n_steps} ${base_path} ${env} ${potential_name} \ No newline at end of file diff --git a/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh b/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh new file mode 100644 index 00000000..07cb3ca7 --- /dev/null +++ b/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh @@ -0,0 +1,30 @@ +######################################## +# idx has to be provided +idx=${1} +######################################## +# idx is an integer between 0 and 5000 +# idx is used to index in this code: + +#names = neutromeratio.parameter_gradients._get_names() +#for name in names: +# for lamb in np.linspace(0, 1, 11): +# protocol.append((name, np.round(lamb, 2))) +#name, lambda_value = protocol[idx-1] +######################################## + +n_steps=400000 # nr of steps +env='vacuum' +potential_name='ANI1ccx' # which potential do you want to use? (ANI1ccx, ANI1x, ANI2x) +echo 'Using potential ' ${potential_name} + +hostname +echo ${idx} +echo ${n_steps} +echo ${env} + +#conda activate conda.env # fire up your conda environment + +base_path="./" # where do you want to save the ouput files +mkdir -p ${base_path} + +python Generate_equilibrium_sampling.py ${idx} ${n_steps} ${base_path} ${env} ${potential_name} \ No newline at end of file diff --git a/scripts/parameter_gradient-opt.py b/scripts/parameter_opt-kfold.py similarity index 100% rename from scripts/parameter_gradient-opt.py rename to scripts/parameter_opt-kfold.py From 5dc66476fdbf9e0ce34c2b0b15da381985fdf4f1 Mon Sep 17 00:00:00 2001 From: wiederm Date: Thu, 29 Oct 2020 16:46:59 +0100 Subject: [PATCH 02/11] providing scripts to generate equilibrium samples with alchemical ANI potential for two tautomer pairs defined by SMILES strings --- neutromeratio/analysis.py | 59 +++++---- neutromeratio/utils.py | 2 +- scripts/Analyse_equilibrium_samples.py | 69 +++++----- scripts/Generate_equilibrium_samples.py | 83 +++++++----- ...ilibrium_samples_for_new_tautomer_pairs.py | 123 ++++++++++-------- .../generate_equilibrium_samples_in_vacuum.sh | 4 +- ...amples_in_vacuum_for_new_tautomer_pairs.sh | 35 ++--- 7 files changed, 205 insertions(+), 170 deletions(-) rename scripts/{ => backup}/generate_equilibrium_samples_in_vacuum.sh (82%) diff --git a/neutromeratio/analysis.py b/neutromeratio/analysis.py index 46b97e73..b3101e14 100644 --- a/neutromeratio/analysis.py +++ b/neutromeratio/analysis.py @@ -1,44 +1,46 @@ import copy import logging import pickle +from glob import glob import matplotlib.pyplot as plt +import mdtraj as md import networkx as nx import numpy as np import pandas as pd +import pkg_resources +import scipy.stats as scs import seaborn as sns +import torch +import torchani from rdkit import Chem, Geometry -from rdkit.Chem import AllChem +from rdkit.Chem import AllChem, rdFMCS from scipy.special import logsumexp from simtk import unit -import mdtraj as md -import torchani -import torch -from rdkit.Chem import rdFMCS -import pkg_resources -import scipy.stats as scs -from .constants import ( - num_threads, - kT, +from neutromeratio.ani import ( + ANI, + AlchemicalANI1ccx, + AlchemicalANI1x, + AlchemicalANI2x, + ANI1_force_and_energy, +) +from neutromeratio.constants import ( + device, + exclude_set_ANI, gas_constant, - temperature, + kT, mols_with_charge, - exclude_set_ANI, multiple_stereobonds, - device, + num_threads, + temperature, ) -from .tautomers import Tautomer -from .parameter_gradients import FreeEnergyCalculator -from .utils import generate_tautomer_class_stereobond_aware, generate_new_tautomer_pair -from .ani import ( - ANI1_force_and_energy, - AlchemicalANI2x, - AlchemicalANI1ccx, - ANI, - AlchemicalANI1x, +from neutromeratio.parameter_gradients import FreeEnergyCalculator +from neutromeratio.tautomers import Tautomer +from neutromeratio.utils import ( + generate_new_tautomer_pair, + generate_tautomer_class_stereobond_aware, ) -from glob import glob logger = logging.getLogger(__name__) @@ -306,7 +308,7 @@ def compare_confomer_generator_and_trajectory_minimum_structures( mol.RemoveAllConformers() # generate energy function, use atom symbols of rdkti mol - from .ani import ANI1ccx, ANI1_force_and_energy + from .ani import ANI1_force_and_energy, ANI1ccx model = ANI1ccx() energy_function = ANI1_force_and_energy( @@ -558,11 +560,14 @@ def setup_alchemical_system_and_energy_function( def setup_new_alchemical_system_and_energy_function( + name: str, + t1_smiles: str, + t2_smiles: str, env: str, ANImodel: ANI, - checkpoint_file: str = "", base_path: str = None, diameter: int = -1, + checkpoint_file: str = "", ): import os @@ -576,7 +581,7 @@ def setup_new_alchemical_system_and_energy_function( #################### # Set up the system, set the restraints - tautomers = generate_new_tautomer_pair(t1_smiles, t2_smiles) + tautomer = generate_new_tautomer_pair(name, t1_smiles, t2_smiles) tautomer.perform_tautomer_transformation() # if base_path is defined write out the topology @@ -653,7 +658,7 @@ def setup_new_alchemical_system_and_energy_function( for r in tautomer.com_restraints: energy_function.add_restraint_to_lambda_protocol(r) - return energy_function, tautomer, flipped + return energy_function, tautomer def _error(x: np.ndarray, y: np.ndarray): diff --git a/neutromeratio/utils.py b/neutromeratio/utils.py index cd8bab92..07f1bf6e 100644 --- a/neutromeratio/utils.py +++ b/neutromeratio/utils.py @@ -9,7 +9,7 @@ from rdkit.Chem import AllChem from simtk import unit -from .constants import kT +from neutromeratio.constants import kT logger = logging.getLogger(__name__) diff --git a/scripts/Analyse_equilibrium_samples.py b/scripts/Analyse_equilibrium_samples.py index 1e907793..ac986df8 100644 --- a/scripts/Analyse_equilibrium_samples.py +++ b/scripts/Analyse_equilibrium_samples.py @@ -4,10 +4,17 @@ import pickle import sys import torch -from neutromeratio.parameter_gradients import FreeEnergyCalculator, setup_mbar -from neutromeratio.constants import kT, device, exclude_set_ANI, mols_with_charge, _get_names +from neutromeratio.parameter_gradients import setup_mbar +from neutromeratio.constants import ( + kT, + device, + exclude_set_ANI, + mols_with_charge, + _get_names, +) from glob import glob from neutromeratio.ani import AlchemicalANI1ccx, AlchemicalANI2x + ####################### ####################### # job idx @@ -19,52 +26,54 @@ ####################### ####################### # read in exp results, smiles and names -exp_results = pickle.load(open('../data/exp_results.pickle', 'rb')) +exp_results = pickle.load(open("../data/exp_results.pickle", "rb")) # name of the system names = _get_names() torch.set_num_threads(4) -name = names[idx-1] -print(name) -print(base_path) -print(potential_name) +name = names[idx - 1] +print(f"Analysing samples for tautomer pair: {name}") +print(f"Saving results in: {base_path}") +print(f"Using potential: {potential_name}") -if potential_name == 'ANI1ccx': +if potential_name == "ANI1ccx": AlchemicalANI = AlchemicalANI1ccx -elif potential_name == 'ANI2x': +elif potential_name == "ANI2x": AlchemicalANI = AlchemicalANI2x else: - raise RuntimeError('Potential needs to be either ANI1ccx or ANI2x') - + raise RuntimeError("Potential needs to be either ANI1ccx or ANI2x") -if env == 'droplet': - print('Env: droplet.') - fec = setup_mbar(name=name, - data_path=base_path, - ANImodel=AlchemicalANI, - bulk_energy_calculation=False, - env='droplet', - max_snapshots_per_window=300, - diameter=18) -elif env == 'vacuum': - print('Env: vacuum.') - fec = setup_mbar(name=name, - data_path=base_path, - ANImodel=AlchemicalANI, - bulk_energy_calculation=True, - env='vacuum', - max_snapshots_per_window=300 +if env == "droplet": + print("Simulating in environment: droplet.") + fec = setup_mbar( + name=name, + data_path=base_path, + ANImodel=AlchemicalANI, + bulk_energy_calculation=False, + env="droplet", + max_snapshots_per_window=300, + diameter=18, + ) +elif env == "vacuum": + print("Simulating in environment: vacuum.") + fec = setup_mbar( + name=name, + data_path=base_path, + ANImodel=AlchemicalANI, + bulk_energy_calculation=True, + env="vacuum", + max_snapshots_per_window=300, ) else: - raise RuntimeError('No env specified. Aborting.') + raise RuntimeError("No env specified. Aborting.") DeltaF_ji, dDeltaF_ji = fec._end_state_free_energy_difference if fec.flipped: DeltaF_ji *= -1 print(DeltaF_ji) -f = open(f"{base_path}/{potential_name}_{env}_rfe_results_in_kT_300snapshots.csv", 'a+') +f = open(f"{base_path}/{potential_name}_{env}_rfe_results_in_kT_300snapshots.csv", "a+") f.write(f"{name}, {DeltaF_ji}, {dDeltaF_ji}\n") f.close() diff --git a/scripts/Generate_equilibrium_samples.py b/scripts/Generate_equilibrium_samples.py index cabb6e89..bf9e1f8d 100644 --- a/scripts/Generate_equilibrium_samples.py +++ b/scripts/Generate_equilibrium_samples.py @@ -11,32 +11,31 @@ idx = int(sys.argv[1]) n_steps = int(sys.argv[2]) -base_path = str(sys.argv[3]) # where to write the results +base_path = str(sys.argv[3]) # where to write the results env = str(sys.argv[4]) model_name = str(sys.argv[5]) -if model_name == 'ANI2x': +if model_name == "ANI2x": model = neutromeratio.ani.AlchemicalANI2x - print(f'Using {model_name}.') -elif model_name == 'ANI1ccx': + print(f"Using {model_name}.") +elif model_name == "ANI1ccx": model = neutromeratio.ani.AlchemicalANI1ccx - print(f'Using {model_name}.') -elif model_name == 'ANI1x': + print(f"Using {model_name}.") +elif model_name == "ANI1x": model = neutromeratio.ani.AlchemicalANI1x - print(f'Using {model_name}.') + print(f"Using {model_name}.") else: - raise RuntimeError(f'Unknown model name: {model_name}') + raise RuntimeError(f"Unknown model name: {model_name}") - -if not (env in ('droplet','vacuum')): - raise RuntimeError('Env must be `droplet` or `vacuum`. Aborting.') +if not (env in ("droplet", "vacuum")): + raise RuntimeError("Env must be `droplet` or `vacuum`. Aborting.") # diameter -if env == 'droplet': +if env == "droplet": diameter_in_angstrom = int(sys.argv[6]) if not diameter_in_angstrom or diameter_in_angstrom < 1: - raise RuntimeError('Diameter must be above 1 Angstrom') + raise RuntimeError("Diameter must be above 1 Angstrom") else: diameter_in_angstrom = -1 @@ -46,9 +45,10 @@ for lamb in np.linspace(0, 1, 11): protocol.append((name, np.round(lamb, 2))) -name, lambda_value = protocol[idx-1] -print(name) -print(lambda_value) +name, lambda_value = protocol[idx - 1] +print(f"Generating samples for tautomer pair: {name}") +print(f"Saving results in: {base_path}") +print(f"Using potential: {model_name}") base_path = f"{base_path}/{name}" os.makedirs(base_path, exist_ok=True) @@ -58,18 +58,21 @@ ANImodel=model, env=env, diameter=diameter_in_angstrom, - base_path=base_path) + base_path=base_path, +) energy_and_force = lambda x: energy_function.calculate_force(x, lambda_value) -if env == 'droplet': +if env == "droplet": x0 = tautomer.get_ligand_in_water_coordinates() - langevin = neutromeratio.LangevinDynamics(atoms=tautomer.ligand_in_water_atoms, - energy_and_force=energy_and_force) -elif env == 'vacuum': + langevin = neutromeratio.LangevinDynamics( + atoms=tautomer.ligand_in_water_atoms, energy_and_force=energy_and_force + ) +elif env == "vacuum": x0 = tautomer.get_hybrid_coordinates() - langevin = neutromeratio.LangevinDynamics(atoms=tautomer.hybrid_atoms, - energy_and_force=energy_and_force) + langevin = neutromeratio.LangevinDynamics( + atoms=tautomer.hybrid_atoms, energy_and_force=energy_and_force + ) else: raise RuntimeError() @@ -77,23 +80,31 @@ x0, e_history = energy_function.minimize(x0, maxiter=5000, lambda_value=lambda_value) -equilibrium_samples, energies, restraint_contribution = langevin.run_dynamics(x0, - n_steps=n_steps, - stepsize=0.5*unit.femtosecond, - progress_bar=True) - -# save equilibrium energy values -for global_list, poperty_name in zip([energies, restraint_contribution], ['energy', 'restraint_energy_contribution']): - f = open(f"{base_path}/{name}_lambda_{lambda_value:0.4f}_{poperty_name}_in_{env}.csv", 'w+') +equilibrium_samples, energies, restraint_contribution = langevin.run_dynamics( + x0, n_steps=n_steps, stepsize=0.5 * unit.femtosecond, progress_bar=True +) + +# save equilibrium energy values +for global_list, poperty_name in zip( + [energies, restraint_contribution], ["energy", "restraint_energy_contribution"] +): + f = open( + f"{base_path}/{name}_lambda_{lambda_value:0.4f}_{poperty_name}_in_{env}.csv", + "w+", + ) for e in global_list[::20]: e_unitless = e / kT - f.write('{}\n'.format(e_unitless)) - f.close() + f.write("{}\n".format(e_unitless)) + f.close() equilibrium_samples = [x[0].value_in_unit(unit.nanometer) for x in equilibrium_samples] -if env == 'vacuum': +if env == "vacuum": ani_traj = md.Trajectory(equilibrium_samples[::20], tautomer.hybrid_topology) else: - ani_traj = md.Trajectory(equilibrium_samples[::20], tautomer.ligand_in_water_topology) + ani_traj = md.Trajectory( + equilibrium_samples[::20], tautomer.ligand_in_water_topology + ) -ani_traj.save(f"{base_path}/{name}_lambda_{lambda_value:0.4f}_in_{env}.dcd", force_overwrite=True) +ani_traj.save( + f"{base_path}/{name}_lambda_{lambda_value:0.4f}_in_{env}.dcd", force_overwrite=True +) diff --git a/scripts/Generate_equilibrium_samples_for_new_tautomer_pairs.py b/scripts/Generate_equilibrium_samples_for_new_tautomer_pairs.py index 19f4e51e..f0b811bd 100644 --- a/scripts/Generate_equilibrium_samples_for_new_tautomer_pairs.py +++ b/scripts/Generate_equilibrium_samples_for_new_tautomer_pairs.py @@ -6,93 +6,112 @@ import sys import os import torch -from neutromeratio.analysis import setup_alchemical_system_and_energy_function +from neutromeratio.analysis import setup_new_alchemical_system_and_energy_function from neutromeratio.ani import AlchemicalANI2x -n_steps = int(sys.argv[2]) -base_path = str(sys.argv[3]) # where to write the results -env = str(sys.argv[4]) -model_name = str(sys.argv[5]) - - -if model_name == 'ANI2x': +smiles1 = str(sys.argv[1]) +smiles2 = str(sys.argv[2]) +name = str(sys.argv[3]) +lambda_value = float(sys.argv[4]) +n_steps = int(sys.argv[5]) +base_path = str(sys.argv[6]) # where to write the results +env = str(sys.argv[7]) +model_name = str(sys.argv[8]) + +print(f"Generating samples for tautomer pair: {name}") +print(f"With SMILES1: {smiles1}") +print(f"With SMILES2: {smiles2}") +print(f"Lambda coupling parameter: {lambda_value} (between 0. and 1.)") +print(f"Saving results in: {base_path}/{name}") +print(f"Using potential: {model_name}") +print(f"Simulating in environment: {env}") + +if model_name == "ANI2x": model = neutromeratio.ani.AlchemicalANI2x - print(f'Using {model_name}.') -elif model_name == 'ANI1ccx': + print(f"Using {model_name}.") +elif model_name == "ANI1ccx": model = neutromeratio.ani.AlchemicalANI1ccx - print(f'Using {model_name}.') -elif model_name == 'ANI1x': + print(f"Using {model_name}.") +elif model_name == "ANI1x": model = neutromeratio.ani.AlchemicalANI1x - print(f'Using {model_name}.') + print(f"Using {model_name}.") else: - raise RuntimeError(f'Unknown model name: {model_name}') + raise RuntimeError(f"Unknown model name: {model_name}") - -if not (env in ('droplet','vacuum')): - raise RuntimeError('Env must be `droplet` or `vacuum`. Aborting.') -# diameter -if env == 'droplet': +# INPUT CHECKS +if not (env in ("droplet", "vacuum")): + raise RuntimeError("Env must be `droplet` or `vacuum`. Aborting.") +if env == "droplet": diameter_in_angstrom = int(sys.argv[6]) if not diameter_in_angstrom or diameter_in_angstrom < 1: - raise RuntimeError('Diameter must be above 1 Angstrom') + raise RuntimeError("Diameter must be above 1 Angstrom") else: diameter_in_angstrom = -1 -protocol = [] -names = neutromeratio.parameter_gradients._get_names() -for name in names: - for lamb in np.linspace(0, 1, 11): - protocol.append((name, np.round(lamb, 2))) - -name, lambda_value = protocol[idx-1] -print(name) -print(lambda_value) - +# Generating base path base_path = f"{base_path}/{name}" os.makedirs(base_path, exist_ok=True) -energy_function, tautomer, flipped = setup_alchemical_system_and_energy_function( +# setting up energy function +energy_function, tautomer = setup_new_alchemical_system_and_energy_function( name=name, - ANImodel=model, + t1_smiles=smiles1, + t2_smiles=smiles2, env=env, + ANImodel=model, diameter=diameter_in_angstrom, - base_path=base_path) + base_path=base_path, +) energy_and_force = lambda x: energy_function.calculate_force(x, lambda_value) -if env == 'droplet': +if env == "droplet": x0 = tautomer.get_ligand_in_water_coordinates() - langevin = neutromeratio.LangevinDynamics(atoms=tautomer.ligand_in_water_atoms, - energy_and_force=energy_and_force) -elif env == 'vacuum': + langevin = neutromeratio.LangevinDynamics( + atoms=tautomer.ligand_in_water_atoms, energy_and_force=energy_and_force + ) +elif env == "vacuum": x0 = tautomer.get_hybrid_coordinates() - langevin = neutromeratio.LangevinDynamics(atoms=tautomer.hybrid_atoms, - energy_and_force=energy_and_force) + langevin = neutromeratio.LangevinDynamics( + atoms=tautomer.hybrid_atoms, energy_and_force=energy_and_force + ) else: raise RuntimeError() +# pytorch number of threads to be used torch.set_num_threads(1) - +# start minimization x0, e_history = energy_function.minimize(x0, maxiter=5000, lambda_value=lambda_value) -equilibrium_samples, energies, restraint_contribution = langevin.run_dynamics(x0, - n_steps=n_steps, - stepsize=0.5*unit.femtosecond, - progress_bar=True) -# save equilibrium energy values -for global_list, poperty_name in zip([energies, restraint_contribution], ['energy', 'restraint_energy_contribution']): - f = open(f"{base_path}/{name}_lambda_{lambda_value:0.4f}_{poperty_name}_in_{env}.csv", 'w+') +# start sampling +equilibrium_samples, energies, restraint_contribution = langevin.run_dynamics( + x0, n_steps=n_steps, stepsize=0.5 * unit.femtosecond, progress_bar=True +) + +# save equilibrium energy values +for global_list, poperty_name in zip( + [energies, restraint_contribution], ["energy", "restraint_energy_contribution"] +): + f = open( + f"{base_path}/{name}_lambda_{lambda_value:0.4f}_{poperty_name}_in_{env}.csv", + "w+", + ) for e in global_list[::20]: e_unitless = e / kT - f.write('{}\n'.format(e_unitless)) - f.close() + f.write("{}\n".format(e_unitless)) + f.close() +# save trajectory equilibrium_samples = [x[0].value_in_unit(unit.nanometer) for x in equilibrium_samples] -if env == 'vacuum': +if env == "vacuum": ani_traj = md.Trajectory(equilibrium_samples[::20], tautomer.hybrid_topology) else: - ani_traj = md.Trajectory(equilibrium_samples[::20], tautomer.ligand_in_water_topology) + ani_traj = md.Trajectory( + equilibrium_samples[::20], tautomer.ligand_in_water_topology + ) -ani_traj.save(f"{base_path}/{name}_lambda_{lambda_value:0.4f}_in_{env}.dcd", force_overwrite=True) +ani_traj.save( + f"{base_path}/{name}_lambda_{lambda_value:0.4f}_in_{env}.dcd", force_overwrite=True +) diff --git a/scripts/generate_equilibrium_samples_in_vacuum.sh b/scripts/backup/generate_equilibrium_samples_in_vacuum.sh similarity index 82% rename from scripts/generate_equilibrium_samples_in_vacuum.sh rename to scripts/backup/generate_equilibrium_samples_in_vacuum.sh index 07cb3ca7..3ba3a6e1 100644 --- a/scripts/generate_equilibrium_samples_in_vacuum.sh +++ b/scripts/backup/generate_equilibrium_samples_in_vacuum.sh @@ -22,9 +22,7 @@ echo ${idx} echo ${n_steps} echo ${env} -#conda activate conda.env # fire up your conda environment - base_path="./" # where do you want to save the ouput files mkdir -p ${base_path} -python Generate_equilibrium_sampling.py ${idx} ${n_steps} ${base_path} ${env} ${potential_name} \ No newline at end of file +python Generate_equilibrium_samples.py ${idx} ${n_steps} ${base_path} ${env} ${potential_name} \ No newline at end of file diff --git a/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh b/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh index 07cb3ca7..96d09aa3 100644 --- a/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh +++ b/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh @@ -1,30 +1,23 @@ -######################################## -# idx has to be provided -idx=${1} -######################################## -# idx is an integer between 0 and 5000 -# idx is used to index in this code: -#names = neutromeratio.parameter_gradients._get_names() -#for name in names: -# for lamb in np.linspace(0, 1, 11): -# protocol.append((name, np.round(lamb, 2))) -#name, lambda_value = protocol[idx-1] +# User specified runtime variables ######################################## - -n_steps=400000 # nr of steps -env='vacuum' +SMILES1='OC1=CC=C2C=CC=CC2=N1' +SMILES2='O=C1NC2=C(C=CC=C2)C=C1' +lambda_value=0.0 # alchemical coupling parameter (between 0. and 1.) --> we recommend using 11 lambda_values +name='SAMP2Lmol4' potential_name='ANI1ccx' # which potential do you want to use? (ANI1ccx, ANI1x, ANI2x) -echo 'Using potential ' ${potential_name} +n_steps=200 # nr of steps (dt = 0.5fs) +######################################## +######################################## -hostname -echo ${idx} -echo ${n_steps} -echo ${env} -#conda activate conda.env # fire up your conda environment +env='vacuum' +echo 'Using potential ' ${potential_name} +echo 'Nr of steps : ' ${n_steps} +echo 'Simulating in : ' ${env} +echo 'Lambda value: ' ${lambda_value} base_path="./" # where do you want to save the ouput files mkdir -p ${base_path} -python Generate_equilibrium_sampling.py ${idx} ${n_steps} ${base_path} ${env} ${potential_name} \ No newline at end of file +python Generate_equilibrium_samples_for_new_tautomer_pairs.py ${SMILES1} ${SMILES2} ${name} ${lambda_value} ${n_steps} ${base_path} ${env} ${potential_name} \ No newline at end of file From 26d706ce52b30103117f7bfe1235a0c8dbf62ead Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Thu, 29 Oct 2020 14:07:32 -0400 Subject: [PATCH 03/11] parameters --> parameter_path in logger --- neutromeratio/ani.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/neutromeratio/ani.py b/neutromeratio/ani.py index 8c1f36fe..d6850f30 100644 --- a/neutromeratio/ani.py +++ b/neutromeratio/ani.py @@ -87,7 +87,7 @@ def load_nn_parameters( else: self.tweaked_neural_network.load_state_dict(parameters) else: - logger.info(f"Parameter file {parameters} does not exist.") + logger.info(f"Parameter file {parameter_path} does not exist.") def _from_neurochem_resources(self, info_file_path, periodic_table_index): ( From 688c6aa81a4541cceb16c62ff1974670e054b516 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Thu, 29 Oct 2020 14:10:05 -0400 Subject: [PATCH 04/11] float32 --> float64 in ani.py TODO: set default precision throughout? --- neutromeratio/ani.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neutromeratio/ani.py b/neutromeratio/ani.py index d6850f30..f45cb0a4 100644 --- a/neutromeratio/ani.py +++ b/neutromeratio/ani.py @@ -775,7 +775,7 @@ def calculate_force( raise RuntimeError(f"Shape of coordinates: {x.shape} is wrong. Aborting.") coordinates = torch.tensor( - x, requires_grad=True, device=self.device, dtype=torch.float32 + x, requires_grad=True, device=self.device, dtype=torch.float64 ) energy_in_kT, restraint_energy_contribution = self._calculate_energy( @@ -921,7 +921,7 @@ def calculate_energy( coordinate_list.value_in_unit(unit.nanometer), requires_grad=requires_grad_wrt_coordinates, device=self.device, - dtype=torch.float32, + dtype=torch.float64, ) logger.debug(f"coordinates tensor: {coordinates.size()}") From 2d06a7570f0783c80ee1e9ee3e17eae23edaa234 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Thu, 29 Oct 2020 14:14:38 -0400 Subject: [PATCH 05/11] add docstring to generate_new_tautomer_pair --- neutromeratio/utils.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/neutromeratio/utils.py b/neutromeratio/utils.py index 07f1bf6e..0e340df1 100644 --- a/neutromeratio/utils.py +++ b/neutromeratio/utils.py @@ -8,6 +8,7 @@ from rdkit import Chem from rdkit.Chem import AllChem from simtk import unit +from neutromeratio import Tautomer from neutromeratio.constants import kT @@ -129,10 +130,9 @@ def get_stereotag_of_stereobonds(smiles: str) -> int: return stereo_tag -def generate_new_tautomer_pair(name: str, t1_smiles: str, t2_smiles: str): - - from neutromeratio import Tautomer - +def generate_new_tautomer_pair(name: str, t1_smiles: str, t2_smiles: str) -> Tautomer: + """Constructs and returns a Tautomer pair object, generating RDKit mols + from t1 and t2 smiles""" return Tautomer( name=name, initial_state_mol=generate_rdkit_mol(t1_smiles), From 3735132e77b8befbb2c5f9abcda607f1a8555fea Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Thu, 29 Oct 2020 14:21:41 -0400 Subject: [PATCH 06/11] update documentation and type hints in geenerate_tautomer_class_stereobond_aware --- neutromeratio/utils.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/neutromeratio/utils.py b/neutromeratio/utils.py index 0e340df1..fb6b84a7 100644 --- a/neutromeratio/utils.py +++ b/neutromeratio/utils.py @@ -1,6 +1,7 @@ import logging import os import random +from typing import List, Tuple, Any, Union import mdtraj as md import nglview @@ -10,6 +11,8 @@ from simtk import unit from neutromeratio import Tautomer +StereoBondType = Union[Any, str] + from neutromeratio.constants import kT logger = logging.getLogger(__name__) @@ -133,6 +136,7 @@ def get_stereotag_of_stereobonds(smiles: str) -> int: def generate_new_tautomer_pair(name: str, t1_smiles: str, t2_smiles: str) -> Tautomer: """Constructs and returns a Tautomer pair object, generating RDKit mols from t1 and t2 smiles""" + # TOOD: should this also accept nr_of_conformations, enforceChirality? return Tautomer( name=name, initial_state_mol=generate_rdkit_mol(t1_smiles), @@ -148,7 +152,7 @@ def generate_tautomer_class_stereobond_aware( t2_smiles: str, nr_of_conformations: int = 1, enforceChirality=True, -) -> list: +) -> Tuple[StereoBondType, List[Tautomer], bool]: """ If a stereobond is present in the tautomer pair we need to transform from the molecule with the stereobond (e.g. enol) to the tautomer without the stereobond (e.g. keto). This @@ -167,6 +171,8 @@ def generate_tautomer_class_stereobond_aware( nr of conformations Returns ---------- + stereobond_type: StereoBondType + one of [None, "generic", "Imine"] tautomers: list a list of tautomer(s) flipped: bool @@ -174,7 +180,6 @@ def generate_tautomer_class_stereobond_aware( """ tautomers = [] - from neutromeratio import Tautomer flipped = False stereobond_type = None From feaf360dbf9bb33b1e77f0102ad3fb3e9cd98109 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Thu, 29 Oct 2020 14:42:00 -0400 Subject: [PATCH 07/11] start to refactor sterebond handler --- neutromeratio/utils.py | 77 ++++++++++++++++++++++++++++++------------ 1 file changed, 55 insertions(+), 22 deletions(-) diff --git a/neutromeratio/utils.py b/neutromeratio/utils.py index fb6b84a7..5f0b6ffe 100644 --- a/neutromeratio/utils.py +++ b/neutromeratio/utils.py @@ -159,6 +159,7 @@ def generate_tautomer_class_stereobond_aware( function makes sure that this happens and returns a list of tautomers. If there is no stereobond present in either tautomer, the list contains only one tautomer. If there is a stereobond present the list contains two tautomers (one with cis, one with trans configuration). + Parameters ---------- name: str @@ -169,6 +170,7 @@ def generate_tautomer_class_stereobond_aware( a SMILES string nr_of_conformations: int nr of conformations + Returns ---------- stereobond_type: StereoBondType @@ -195,11 +197,18 @@ def generate_tautomer_class_stereobond_aware( t1_kappa_1 = change_stereobond_in_imine_to_trans( generate_rdkit_mol(t1_smiles) ) + + t1_a = t1_kappa_0 + t2_a = generate_rdkit_mol(t2_smiles) + + t1_b = t1_kappa_1 + t2_b = generate_rdkit_mol(t2_smiles) + tautomers.append( Tautomer( name=name, - initial_state_mol=t1_kappa_0, - final_state_mol=generate_rdkit_mol(t2_smiles), + initial_state_mol=t1_a, + final_state_mol= t2_a, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) @@ -207,8 +216,8 @@ def generate_tautomer_class_stereobond_aware( tautomers.append( Tautomer( name=name, - initial_state_mol=t1_kappa_1, - final_state_mol=generate_rdkit_mol(t2_smiles), + initial_state_mol=t1_b, + final_state_mol=t2_b, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) @@ -222,11 +231,18 @@ def generate_tautomer_class_stereobond_aware( t2_kappa_1 = change_stereobond_in_imine_to_trans( generate_rdkit_mol(t2_smiles) ) + + t1_a = t2_kappa_0 + t2_a = generate_rdkit_mol(t1_smiles) + + t1_b = t2_kappa_1 + t2_b = generate_rdkit_mol(t1_smiles) + tautomers.append( Tautomer( name=name, - initial_state_mol=t2_kappa_0, - final_state_mol=generate_rdkit_mol(t1_smiles), + initial_state_mol=t1_a, + final_state_mol=t2_a, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) @@ -234,8 +250,8 @@ def generate_tautomer_class_stereobond_aware( tautomers.append( Tautomer( name=name, - initial_state_mol=t2_kappa_1, - final_state_mol=generate_rdkit_mol(t1_smiles), + initial_state_mol=t1_b, + final_state_mol=t2_b, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) @@ -245,6 +261,9 @@ def generate_tautomer_class_stereobond_aware( raise RuntimeError("Stereobonds present in both tautomers ... aborting!") elif get_nr_of_stereobonds(t1_smiles) == get_nr_of_stereobonds(t2_smiles): + t1 = generate_rdkit_mol(t1_smiles) + t2 = generate_rdkit_mol(t2_smiles) + if ( get_nr_of_stereobonds(t1_smiles) == 0 and get_nr_of_stereobonds(t2_smiles) == 0 @@ -255,8 +274,8 @@ def generate_tautomer_class_stereobond_aware( tautomers.append( Tautomer( name=name, - initial_state_mol=generate_rdkit_mol(t1_smiles), - final_state_mol=generate_rdkit_mol(t2_smiles), + initial_state_mol=t1, + final_state_mol=t2, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) @@ -266,8 +285,8 @@ def generate_tautomer_class_stereobond_aware( tautomers.append( Tautomer( name=name, - initial_state_mol=generate_rdkit_mol(t1_smiles), - final_state_mol=generate_rdkit_mol(t2_smiles), + initial_state_mol=t1, + final_state_mol=t2, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) @@ -275,18 +294,24 @@ def generate_tautomer_class_stereobond_aware( else: # stereobonds on both endstates - # we need to add a torsion bias to make sure that the lambda protocol stopp at the correct torsion - + # we need to add a torsion bias to make sure that the lambda protocol stops at the correct torsion raise RuntimeError("Two stereobonds ... aborting") elif get_nr_of_stereobonds(t1_smiles) > get_nr_of_stereobonds(t2_smiles): stereobond_type = "generic" t1_smiles_kappa_0 = t1_smiles t1_smiles_kappa_1 = change_only_stereobond(t1_smiles) + + t1_a = generate_rdkit_mol(t1_smiles_kappa_0) + t2_a = generate_rdkit_mol(t2_smiles) + + t1_b = generate_rdkit_mol(t1_smiles_kappa_1) + t2_b = generate_rdkit_mol(t2_smiles) + tautomers.append( Tautomer( name=name, - initial_state_mol=generate_rdkit_mol(t1_smiles_kappa_0), - final_state_mol=generate_rdkit_mol(t2_smiles), + initial_state_mol=t1_a, + final_state_mol=t2_a, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) @@ -294,8 +319,8 @@ def generate_tautomer_class_stereobond_aware( tautomers.append( Tautomer( name=name, - initial_state_mol=generate_rdkit_mol(t1_smiles_kappa_1), - final_state_mol=generate_rdkit_mol(t2_smiles), + initial_state_mol=t1_b, + final_state_mol=t2_b, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) @@ -307,20 +332,28 @@ def generate_tautomer_class_stereobond_aware( t1_smiles_kappa_0 = t2_smiles t1_smiles_kappa_1 = change_only_stereobond(t2_smiles) t2_smiles = t1_smiles + + t1_a = generate_rdkit_mol(t1_smiles_kappa_0) + t2_a = generate_rdkit_mol(t2_smiles) + + t1_b = generate_rdkit_mol(t1_smiles_kappa_1) + t2_b = generate_rdkit_mol(t2_smiles) + tautomers.append( Tautomer( name=name, - initial_state_mol=generate_rdkit_mol(t1_smiles_kappa_0), - final_state_mol=generate_rdkit_mol(t2_smiles), + initial_state_mol=t1_a, + final_state_mol=t2_a, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) ) + tautomers.append( Tautomer( name=name, - initial_state_mol=generate_rdkit_mol(t1_smiles_kappa_1), - final_state_mol=generate_rdkit_mol(t2_smiles), + initial_state_mol=t1_b, + final_state_mol=t2_b, nr_of_conformations=nr_of_conformations, enforceChirality=enforceChirality, ) From b1a970b2f5c7122cf5dbae6f7aec58d6ce8f86fc Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Thu, 29 Oct 2020 14:49:02 -0400 Subject: [PATCH 08/11] reduce code duplication in stereobond handler --- neutromeratio/utils.py | 112 ++++++++--------------------------------- 1 file changed, 21 insertions(+), 91 deletions(-) diff --git a/neutromeratio/utils.py b/neutromeratio/utils.py index 5f0b6ffe..49ad5233 100644 --- a/neutromeratio/utils.py +++ b/neutromeratio/utils.py @@ -179,6 +179,11 @@ def generate_tautomer_class_stereobond_aware( a list of tautomer(s) flipped: bool to indicate that t1/t2 SMILES have been exchanged + + + Notes + ----- + one branch (for detecting stereobonds in heterocycles) depends on name """ tautomers = [] @@ -186,6 +191,10 @@ def generate_tautomer_class_stereobond_aware( flipped = False stereobond_type = None + def _tautomer(t1, t2): + """Tautomer constructor, with name, nr_of_conformations, enforceChirality applied""" + return Tautomer(name, t1, t2, nr_of_conformations, enforceChirality) + if flag_unspec_stereo(t1_smiles) or flag_unspec_stereo(t2_smiles): logger.debug("Imines present ... switching to imine generation.") stereobond_type = "Imine" @@ -204,24 +213,9 @@ def generate_tautomer_class_stereobond_aware( t1_b = t1_kappa_1 t2_b = generate_rdkit_mol(t2_smiles) - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1_a, - final_state_mol= t2_a, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1_b, - final_state_mol=t2_b, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) + tautomers.append(_tautomer(t1_a, t2_a)) + tautomers.append(_tautomer(t1_b, t2_b)) + elif flag_unspec_stereo(t2_smiles): flipped = True @@ -238,24 +232,8 @@ def generate_tautomer_class_stereobond_aware( t1_b = t2_kappa_1 t2_b = generate_rdkit_mol(t1_smiles) - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1_a, - final_state_mol=t2_a, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1_b, - final_state_mol=t2_b, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) + tautomers.append(_tautomer(t1_a, t2_a)) + tautomers.append(_tautomer(t1_b, t2_b)) else: raise RuntimeError("Stereobonds present in both tautomers ... aborting!") @@ -271,26 +249,11 @@ def generate_tautomer_class_stereobond_aware( # no stereobond -- normal protocol # generate both rdkit mol logger.debug("No stereobonds ...") - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1, - final_state_mol=t2, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) + tautomers.append(_tautomer(t1, t2)) + elif name == "molDWRow_1636": logger.debug("molDWRow_1636 -- stereobonds in hetereocycle ...") - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1, - final_state_mol=t2, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) + tautomers.append(_tautomer(t1, t2)) else: # stereobonds on both endstates @@ -307,24 +270,8 @@ def generate_tautomer_class_stereobond_aware( t1_b = generate_rdkit_mol(t1_smiles_kappa_1) t2_b = generate_rdkit_mol(t2_smiles) - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1_a, - final_state_mol=t2_a, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1_b, - final_state_mol=t2_b, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) + tautomers.append(_tautomer(t1_a, t2_a)) + tautomers.append(_tautomer(t1_b, t2_b)) elif get_nr_of_stereobonds(t1_smiles) < get_nr_of_stereobonds(t2_smiles): stereobond_type = "generic" @@ -339,25 +286,8 @@ def generate_tautomer_class_stereobond_aware( t1_b = generate_rdkit_mol(t1_smiles_kappa_1) t2_b = generate_rdkit_mol(t2_smiles) - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1_a, - final_state_mol=t2_a, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) - - tautomers.append( - Tautomer( - name=name, - initial_state_mol=t1_b, - final_state_mol=t2_b, - nr_of_conformations=nr_of_conformations, - enforceChirality=enforceChirality, - ) - ) + tautomers.append(_tautomer(t1_a, t2_a)) + tautomers.append(_tautomer(t1_b, t2_b)) else: raise RuntimeError("Stereobonds present in both tautomers ... aborting!") From 1168b955009487093ba4b00ecf1d070c662c483e Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Thu, 29 Oct 2020 16:03:55 -0400 Subject: [PATCH 09/11] revert Tautomer cyclic-import --- neutromeratio/utils.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/neutromeratio/utils.py b/neutromeratio/utils.py index 49ad5233..26b31fbc 100644 --- a/neutromeratio/utils.py +++ b/neutromeratio/utils.py @@ -9,7 +9,6 @@ from rdkit import Chem from rdkit.Chem import AllChem from simtk import unit -from neutromeratio import Tautomer StereoBondType = Union[Any, str] @@ -133,10 +132,11 @@ def get_stereotag_of_stereobonds(smiles: str) -> int: return stereo_tag -def generate_new_tautomer_pair(name: str, t1_smiles: str, t2_smiles: str) -> Tautomer: +def generate_new_tautomer_pair(name: str, t1_smiles: str, t2_smiles: str): """Constructs and returns a Tautomer pair object, generating RDKit mols from t1 and t2 smiles""" # TOOD: should this also accept nr_of_conformations, enforceChirality? + from neutromeratio.tautomers import Tautomer return Tautomer( name=name, initial_state_mol=generate_rdkit_mol(t1_smiles), @@ -152,7 +152,7 @@ def generate_tautomer_class_stereobond_aware( t2_smiles: str, nr_of_conformations: int = 1, enforceChirality=True, -) -> Tuple[StereoBondType, List[Tautomer], bool]: +) -> Tuple[StereoBondType, List, bool]: """ If a stereobond is present in the tautomer pair we need to transform from the molecule with the stereobond (e.g. enol) to the tautomer without the stereobond (e.g. keto). This @@ -185,6 +185,7 @@ def generate_tautomer_class_stereobond_aware( ----- one branch (for detecting stereobonds in heterocycles) depends on name """ + from neutromeratio.tautomers import Tautomer tautomers = [] From b2692190681288746cdfbc04f8af8649ef29c5d8 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Thu, 29 Oct 2020 16:11:38 -0400 Subject: [PATCH 10/11] revert float32 --> float64 --- neutromeratio/ani.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neutromeratio/ani.py b/neutromeratio/ani.py index f45cb0a4..d6850f30 100644 --- a/neutromeratio/ani.py +++ b/neutromeratio/ani.py @@ -775,7 +775,7 @@ def calculate_force( raise RuntimeError(f"Shape of coordinates: {x.shape} is wrong. Aborting.") coordinates = torch.tensor( - x, requires_grad=True, device=self.device, dtype=torch.float64 + x, requires_grad=True, device=self.device, dtype=torch.float32 ) energy_in_kT, restraint_energy_contribution = self._calculate_energy( @@ -921,7 +921,7 @@ def calculate_energy( coordinate_list.value_in_unit(unit.nanometer), requires_grad=requires_grad_wrt_coordinates, device=self.device, - dtype=torch.float64, + dtype=torch.float32, ) logger.debug(f"coordinates tensor: {coordinates.size()}") From d8bc9d70e49da2df2df024343d1ae904dc0b26d4 Mon Sep 17 00:00:00 2001 From: wiederm Date: Fri, 30 Oct 2020 19:09:02 +0100 Subject: [PATCH 11/11] adding script to generate equilibrium samples, analyse dG and reweight based on optimized parameters --- neutromeratio/parameter_gradients.py | 105 +++++++++++++++++- neutromeratio/tests/test_neutromeratio.py | 5 +- .../Calculate_dG_for_new_tautomer_pairs.py | 82 ++++++++++++++ ...amples_in_vacuum_for_new_tautomer_pairs.sh | 23 ---- .../generate_samples_and_analyse_results.sh | 22 ++++ 5 files changed, 209 insertions(+), 28 deletions(-) create mode 100644 scripts/Calculate_dG_for_new_tautomer_pairs.py delete mode 100644 scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh create mode 100644 scripts/generate_samples_and_analyse_results.sh diff --git a/neutromeratio/parameter_gradients.py b/neutromeratio/parameter_gradients.py index a74c3fb1..9bcfdc95 100644 --- a/neutromeratio/parameter_gradients.py +++ b/neutromeratio/parameter_gradients.py @@ -104,9 +104,9 @@ def get_mix(lambda0, lambda1, lam=0.0): f"There are {len(snapshots)} snapshots per lambda state (max: {max_snapshots_per_window}). Aborting." ) - # test that we have not less than 80% of max_snapshots_per_window + # test that we have not less than 60% of max_snapshots_per_window if max_snapshots_per_window != -1 and len(snapshots) < ( - int(max_snapshots_per_window * 0.8) + int(max_snapshots_per_window * 0.6) ): raise RuntimeError( f"There are only {len(snapshots)} snapshots per lambda state. Aborting." @@ -125,7 +125,7 @@ def get_mix(lambda0, lambda1, lam=0.0): snapshots.extend(ani_trajs[lam]) logger.debug(f"Snapshots per lambda {lam}: {len(ani_trajs[lam])}") - if len(snapshots) < 100: + if len(snapshots) < 300: logger.critical( f"Total number of snapshots is {len(snapshots)} -- is this enough?" ) @@ -1065,3 +1065,102 @@ def parse_lambda_from_dcd_filename(dcd_filename): fec.flipped = flipped return fec + + +def setup_mbar_for_new_tautomer_pairs( + name: str, + t1_smiles: str, + t2_smiles: str, + max_snapshots_per_window: int, + ANImodel: ANI, + bulk_energy_calculation: bool, + env: str = "vacuum", + checkpoint_file: str = "", + data_path: str = "../data/", + diameter: int = -1, +): + + from neutromeratio.analysis import setup_new_alchemical_system_and_energy_function + import os + + if not (env == "vacuum" or env == "droplet"): + raise RuntimeError("Only keyword vacuum or droplet are allowed as environment.") + if env == "droplet" and diameter == -1: + raise RuntimeError("Something went wrong.") + + def parse_lambda_from_dcd_filename(dcd_filename): + """parsed the dcd filename + + Arguments: + dcd_filename {str} -- how is the dcd file called? + + Returns: + [float] -- lambda value + """ + l = dcd_filename[: dcd_filename.find(f"_energy_in_{env}")].split("_") + lam = l[-3] + return float(lam) + + data_path = os.path.abspath(data_path) + if not os.path.exists(data_path): + raise RuntimeError(f"{data_path} does not exist!") + + ####################### + (energy_function, tautomer,) = setup_new_alchemical_system_and_energy_function( + name=name, + t1_smiles=t1_smiles, + t2_smiles=t2_smiles, + ANImodel=ANImodel, + checkpoint_file=checkpoint_file, + env=env, + diameter=diameter, + base_path=f"{data_path}/{name}/", + ) + # and lambda values in list + dcds = glob(f"{data_path}/{name}/*.dcd") + + lambdas = [] + md_trajs = [] + energies = [] + + # read in all the frames from the trajectories + if env == "vacuum": + top = tautomer.hybrid_topology + else: + top = f"{data_path}/{name}/{name}_in_droplet.pdb" + + for dcd_filename in dcds: + lam = parse_lambda_from_dcd_filename(dcd_filename) + lambdas.append(lam) + traj = md.load_dcd(dcd_filename, top=top) + logger.debug(f"Nr of frames in trajectory: {len(traj)}") + md_trajs.append(traj) + f = open( + f"{data_path}/{name}/{name}_lambda_{lam:0.4f}_energy_in_{env}.csv", "r" + ) + energies.append(np.array([float(e) * kT for e in f])) + f.close() + + if len(lambdas) < 5: + raise RuntimeError(f"Below 5 lambda states for {name}") + + assert len(lambdas) == len(energies) + assert len(lambdas) == len(md_trajs) + + if env == "vacuum": + pickle_path = f"{data_path}/{name}/{name}_{ANImodel.name}_{max_snapshots_per_window}_{len(tautomer.hybrid_atoms)}_atoms.pickle" + else: + pickle_path = f"{data_path}/{name}/{name}_{ANImodel.name}_{max_snapshots_per_window}_{diameter}A_{len(tautomer.ligand_in_water_atoms)}_atoms.pickle" + + # calculate free energy in kT + fec = FreeEnergyCalculator( + ani_model=energy_function, + md_trajs=md_trajs, + potential_energy_trajs=energies, + lambdas=lambdas, + pickle_path=pickle_path, + bulk_energy_calculation=bulk_energy_calculation, + max_snapshots_per_window=max_snapshots_per_window, + ) + + return fec diff --git a/neutromeratio/tests/test_neutromeratio.py b/neutromeratio/tests/test_neutromeratio.py index e00f75b5..2e9c405c 100644 --- a/neutromeratio/tests/test_neutromeratio.py +++ b/neutromeratio/tests/test_neutromeratio.py @@ -1324,7 +1324,6 @@ def test_setup_mbar_test_pickle_files(): assert np.isclose(-3.2194223855155357, fec._end_state_free_energy_difference[0]) - def test_change_stereobond(): from ..utils import change_only_stereobond, get_nr_of_stereobonds @@ -1748,8 +1747,10 @@ def test_load_parameters(): ) ): # set tweaked parameters - print(model_name) + model._reset_parameters() + model_instance = model([0, 0]) + model_instance._reset_parameters() # initial parameters params1 = list(model_instance.original_neural_network.parameters())[6][ 0 diff --git a/scripts/Calculate_dG_for_new_tautomer_pairs.py b/scripts/Calculate_dG_for_new_tautomer_pairs.py new file mode 100644 index 00000000..8d4ba258 --- /dev/null +++ b/scripts/Calculate_dG_for_new_tautomer_pairs.py @@ -0,0 +1,82 @@ +import neutromeratio +from simtk import unit +import numpy as np +import pickle +import sys +import torch +from neutromeratio.parameter_gradients import setup_mbar_for_new_tautomer_pairs +from neutromeratio.constants import ( + kT, + device, + exclude_set_ANI, + mols_with_charge, + _get_names, +) +from glob import glob +from neutromeratio.ani import AlchemicalANI1ccx, AlchemicalANI2x + +####################### +####################### +smiles1 = str(sys.argv[1]) +smiles2 = str(sys.argv[2]) +name = str(sys.argv[3]) +base_path = str(sys.argv[4]) +env = str(sys.argv[5]) +potential_name = str(sys.argv[6]) +####################### +####################### +torch.set_num_threads(4) + +print(f"Analysing samples for tautomer pair: {name}") +print(f"With SMILES1: {smiles1}") +print(f"With SMILES2: {smiles2}") +print(f"Saving results in: {base_path}") +print(f"Using potential: {potential_name}") + +if potential_name == "ANI1ccx": + AlchemicalANI = AlchemicalANI1ccx +elif potential_name == "ANI2x": + AlchemicalANI = AlchemicalANI2x +else: + raise RuntimeError("Potential needs to be either ANI1ccx or ANI2x") + +if env == "droplet": + print("Simulating in environment: droplet.") + fec = setup_mbar_for_new_tautomer_pairs( + name=name, + t1_smiles=smiles1, + t2_smiles=smiles2, + data_path=base_path, + ANImodel=AlchemicalANI, + bulk_energy_calculation=False, + env="droplet", + max_snapshots_per_window=300, + diameter=-1, + ) +elif env == "vacuum": + print("Simulating in environment: vacuum.") + fec = setup_mbar_for_new_tautomer_pairs( + name=name, + t1_smiles=smiles1, + t2_smiles=smiles2, + data_path=base_path, + ANImodel=AlchemicalANI, + bulk_energy_calculation=True, + env="vacuum", + max_snapshots_per_window=300, + checkpoint_file="../data/retraining/parameters_ANI1ccx_vacuum_best.pt", + ) +else: + raise RuntimeError("No env specified. Aborting.") + +DeltaF_ji, dDeltaF_ji = fec._end_state_free_energy_difference +print(f"dG obtained from vacuum simulations: {DeltaF_ji} with {dDeltaF_ji} [kT]") + +if potential_name == "ANI1ccx": + DeltaF_ji = fec._compute_free_energy_difference() + print( + f"dG obtained with optimized parameters from vacuum simulations: {DeltaF_ji} [kT]" + ) + +else: + print("We only provide optimized parameters for ANI1ccx.") diff --git a/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh b/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh deleted file mode 100644 index 96d09aa3..00000000 --- a/scripts/generate_equilibrium_samples_in_vacuum_for_new_tautomer_pairs.sh +++ /dev/null @@ -1,23 +0,0 @@ - -# User specified runtime variables -######################################## -SMILES1='OC1=CC=C2C=CC=CC2=N1' -SMILES2='O=C1NC2=C(C=CC=C2)C=C1' -lambda_value=0.0 # alchemical coupling parameter (between 0. and 1.) --> we recommend using 11 lambda_values -name='SAMP2Lmol4' -potential_name='ANI1ccx' # which potential do you want to use? (ANI1ccx, ANI1x, ANI2x) -n_steps=200 # nr of steps (dt = 0.5fs) -######################################## -######################################## - - -env='vacuum' -echo 'Using potential ' ${potential_name} -echo 'Nr of steps : ' ${n_steps} -echo 'Simulating in : ' ${env} -echo 'Lambda value: ' ${lambda_value} - -base_path="./" # where do you want to save the ouput files -mkdir -p ${base_path} - -python Generate_equilibrium_samples_for_new_tautomer_pairs.py ${SMILES1} ${SMILES2} ${name} ${lambda_value} ${n_steps} ${base_path} ${env} ${potential_name} \ No newline at end of file diff --git a/scripts/generate_samples_and_analyse_results.sh b/scripts/generate_samples_and_analyse_results.sh new file mode 100644 index 00000000..0137559a --- /dev/null +++ b/scripts/generate_samples_and_analyse_results.sh @@ -0,0 +1,22 @@ +# This is a simple wrapper script that generates equilibrium samples for a given tautomer set +# and estimates the free energy difference with ANI1cccx and the tautomer optimized parameter set +################################################ +SMILES1='OC1=CC=C2C=CC=CC2=N1' +SMILES2='O=C1NC2=C(C=CC=C2)C=C1' +name='test_mol' # defines where the output directory name +base_path="./" # where do you want to save the ouput files -> the ouput will have the form ${base_path}/${name} +potential_name='ANI1ccx' # which potential do you want to use? (ANI1ccx, ANI1x, ANI2x) +n_steps=10000 # nr of steps (dt = 0.5fs) +env='vacuum' +################################################ +echo 'Using potential ' ${potential_name} +echo 'Nr of steps : ' ${n_steps} +echo 'Simulating in : ' ${env} +echo 'Lambda value: ' ${lambda_value} +mkdir -p ${base_path} + +for lambda_value in $(seq 0 0.2 1); # using 6 lambda states +do python Generate_equilibrium_samples_for_new_tautomer_pairs.py ${SMILES1} ${SMILES2} ${name} ${lambda_value} ${n_steps} ${base_path} ${env} ${potential_name}; +done + +python Calculate_dG_for_new_tautomer_pairs.py ${SMILES1} ${SMILES2} ${name} ${base_path} ${env} ${potential_name};