From 9cc88f79f4c9cbff256b8366b9728e9aef86abc8 Mon Sep 17 00:00:00 2001
From: Jesse Huhtala <jesse.j.huhtala@utu.fi>
Date: Mon, 30 Nov 2020 00:40:22 +0200
Subject: [PATCH] attempted implementation of another functional

---
 basis.f90    |   2 +-
 energies.f90 |  82 +++++++++++++++------------------------------------
 minim        | Bin 48104 -> 43904 bytes
 minim.f90    |  13 ++++----
 optimize.f90 |  21 ++++++++-----
 5 files changed, 45 insertions(+), 73 deletions(-)

diff --git a/basis.f90 b/basis.f90
index de1400e..bc4a8ae 100644
--- a/basis.f90
+++ b/basis.f90
@@ -32,7 +32,7 @@ function create_basis(nelectrons,nfuns,ngrid,L) result(bas)
 		x(i) = (i) * ep
 	end do
     do i=1, nfuns
-        f(:,i) = sqrt(2/L) * sin(i*pi*x(:)/L)
+        f(:,i) = sqrt(2/L) * sin(i*pi*x(:)/L)!+1.0_dp
         fder(:,i) = sqrt(2/L) * i*pi/L *cos(i*pi*x(:)/L)
     end do
     bas%nfuns = nfuns
diff --git a/energies.f90 b/energies.f90
index 8571d0a..a8a1a81 100644
--- a/energies.f90
+++ b/energies.f90
@@ -55,28 +55,7 @@ function get_exchange_energy(bas,nx) result(ex)
 
     
 end function
-function get_tf_pot(bas,alpha,sn,vec,mu) result(k)
-    type(basis_t)           :: bas
-    real(dp)                :: k,sn(:),vec(:),alpha,beta, kindens(size(bas%x)), deldens(size(bas%x)),mu,const
-    real(dp)                :: psi(size(bas%x)), pi
-    pi = 4.0_dp*atan(1.0_dp)
-    psi = expand_psi_ofdft(bas,vec(:)+alpha*sn(:))
 
-    k = sum( abs(2.0_dp*psi(:)*(0.5_dp * psi(:)**2.0_dp+ pi**2.0_dp/2.0_dp *psi(:)**4.0_dp-mu) )**2.0_dp*bas%dx)
-    !print *, k 
-end function
-function get_tf_pot_grad(bas,vec,mu) result(grad)
-    type(basis_t)           :: bas
-    real(dp)                :: grad(bas%nfuns), vec(:), mu, psi(size(bas%x)),pi
-    integer                 :: i
-    pi = 4.0_dp*atan(1.0_dp)
-    psi = expand_psi_ofdft(bas,vec)
-    grad = 0.0_dp
-    do i=1, bas%nfuns
-        grad(i) = grad(i)+sum(bas%dx*(4.0_dp*bas%f(:,i)*(0.5_dp * psi(:)**2.0_dp+ pi**2.0_dp/2.0_dp *psi(:)**4.0_dp-mu)))
-        grad(i) = grad(i)+sum(bas%dx*(4.0_dp*psi(:)*(bas%f(:,i)*psi(:)+2.0_dp*pi**2.0_dp *bas%f(:,i)*psi(:)**3.0_dp)) )
-    end do
-end function
 function get_nagy_total_energy(bas,alpha,sn,vec,mu) result(k)
     type(basis_t)           :: bas
     real(dp)                :: k, g(size(bas%x)), beta(size(bas%x)),pi,alpha,mu,rho(size(bas%x)), drho(size(bas%x))
@@ -104,38 +83,31 @@ function get_nagy_total_energy(bas,alpha,sn,vec,mu) result(k)
     g(:) = beta(:)/rho(:)
 
     drho(:) = derivative(bas,rho)
-    nonloc(:) =  2*g(:)*drho(:)
-    !goto 10
+    nonloc(:) = sum(derivative(bas,psi(:)**2.0_dp)*psi(:)**2.0_dp * bas%dx)/psi(:)**2.0_dp
     
-    do i=size(bas%x)/2+1, size(bas%x)
-        do j=1, size(bas%x)-i+1
-            xint(i) = xint(i) - g(i+j-1)*drho(i+j-1)*psi(i+j-1)**4.0_dp*dx 
-            !xint(i)  = xint(i) - 2.0_dp*g(i+j-1)*drho(i+j-1)*rho(i+j-1)**2.0_dp
-        end do
-    end do
-    xint(size(bas%x)/2) = xint(size(bas%x)/2+1)
-    do i=1, size(bas%x)/2-1
-        xint(i) = xint( size(bas%x)-i+1)
-    end do
-    !10 continue
-    xint(:) = xint(:)/(rho**2.0_dp)
-    !goto 20
+    
+    !goto 10
+
+    goto 20
     open(12, file="gnonloc.dt")
     
     do i=1,size(bas%x)
-        write(12,*) bas%x(i), xint(i),g(i),drho(i),psi(i)**4
+        write(12,*) bas%x(i), nonloc(i),g(i),drho(i),psi(i)**2.0_dp
     end do
     close(12)
-    !20 continue
+    stop
+    20 continue
     k = sum(( abs( derivative(bas,psi(:)**2.0_dp) ) )**2.0_dp/(psi(:)**2.0_dp)*dx)*1.0_dp/(8.0_dp)
     k = k+mupar*sum(psi(:)**2.0_dp*dx) - 2.0_dp * sum(g(:)*psi(:)**2.0_dp*dx) &
-        +sum(0.25_dp*psi(:)**4.0_dp*dx)-mu*sum(psi(:)**2.0_dp*dx)+ sum(xint(:)*dx)
-    !print *, k
+        +sum(0.25_dp*psi(:)**4.0_dp*dx)-mu*sum(psi(:)**2.0_dp*dx)-sum(nonloc(:)*bas%dx)
+    !print *, "here"
+    !print *, k 
+    !stop
 end function
 function get_nagy_gradient(bas,vec,mu) result(grad)
     type(basis_t)           :: bas
     real(dp)                :: vec(:),mu,psi(size(bas%x)),grad(bas%nfuns),xint(size(bas%x)),drho(size(bas%x)),mupar
-    real(dp)                :: nonloc(size(bas%x)),rho(size(bas%x)),beta(size(bas%x)),dx,g(size(bas%x)),pi
+    real(dp)                :: nonloc(size(bas%x),bas%nfuns),rho(size(bas%x)),beta(size(bas%x)),dx,g(size(bas%x)),pi
     integer                 :: i,j,k,z
     !Beta term for non-interacting system
     beta = 0.0_dp
@@ -160,26 +132,17 @@ function get_nagy_gradient(bas,vec,mu) result(grad)
     g(:) = beta(:)/rho(:)
     
     drho(:) = derivative(bas,rho)
-    nonloc(:) =  2*g(:)*drho(:)
+    nonloc = 0.0_dp
+    do i=1, bas%nfuns
+        nonloc(:,i) = nonloc(:,i)+ 2.0_dp*sum((bas%fder(:,i)*psi(:)+&
+                      bas%f(:,i)*derivative(bas,psi(:)) )*psi(:)**2.0_dp*bas%dx)/psi(:)**2.0_dp
+        nonloc(:,i) = nonloc(:,i) + sum(bas%f(:,i)*derivative(bas,psi(:)**2.0_dp)*bas%dx)/psi(:)**2.0_dp
+        nonloc(:,i) = nonloc(:,i) - 2.0_dp*bas%f(:,i)*sum(psi(:)**2.0_dp*derivative(bas,psi(:)**2.0_dp)*bas%dx)/psi(:)**3.0_dp
+    end do
     !grad(:) = get_vw_gradient
     !goto 10
-    do k=1, bas%nfuns
-        do i=size(bas%x)/2+1, size(bas%x)
-            do j=1, size(bas%x)-i+1
-                !print *, i
-                xint(i) = xint(i) - 4.0_dp*g(i+j-1)*drho(i+j-1)*bas%f(i+j-1,k)*psi(i+j-1)**3.0_dp*dx
-            end do
-        end do
-        xint(size(bas%x)/2) = xint(size(bas%x)/2+1)
-        do z=1,size(bas%x)/2-1
-            xint(z) = xint( size(bas%x)-z+1)
-        end do
-        grad(k) = grad(k) + sum(xint(:)/(rho(:)**2.0_dp) *dx)
-        xint = 0.0_dp
-        
-        !print *, k
-    end do
-    !10 continue
+
+    !10 continu
     do i=1, bas%nfuns
         !print *, i
         grad(i) = grad(i) + sum(2.0_dp*mupar * bas%f(:,i) *psi(:)* dx)-sum(4.0_dp*g(:)*bas%f(:,i)*psi(:)*dx)
@@ -191,6 +154,7 @@ function get_nagy_gradient(bas,vec,mu) result(grad)
         grad(i) = grad(i) - 2.0_dp/(8.0_dp) *&
          sum(bas%dx * bas%f(:,i)*abs(derivative(bas,psi**2.0_dp))**2.0_dp/(psi**3.0_dp))
 
+        grad(i) = grad(i) - sum(nonloc(:,i)*bas%dx)
     end do
     !print *, "we here?"
 end function
diff --git a/minim b/minim
index 5b3f8192f8a88289efbd7abc5a889dc21a92bf31..fd105a2c92a3c97ec465f96acd9ed19a9f6e8084 100755
GIT binary patch
delta 16805
zcmaFyovGnE(*zABhb0rW?$x*EF@OP#W)J{#85kKD7+65e1`xr(z<@?ua6`mlG_pQ6
z1_p+VUWk0gTvU1kGlc&DY&gT@c19sVsL~3k(u%p0S24yZo@j=sM5pgSgc<Hk)?^BR
z=|xyGxr@n5<ot0>+qyRCSF2oe+_&7S3{Oq4m)rb;Nsz4`Y6l0@kN{x_qX12u0V;lC
z0Yo0A9>)KHCjI~xJPT3PCqT{NfQr|{3@VU=D7-NrLW9B#5+o4*10-<~h!B{3fg~;l
z7J(2Ski^9yLSXU-lDGr|Lp@j!L`Xn=z|Oz`j%<(^2rD3oGl2x4SOZBM8i!yR10-=)
zun2^(KoSRe4kB95z~F!+0kRw-!NA~wBn}HDkX!(gI44K|iX)K3xu9YoDgj9xnn*#y
z3>hGCkpDrJ@_>aPL;*;Gfq?-e&I=I&lNCtfd|(j>(SRh*4-o>B9Z2HPTm%-MfFvF!
z2o{16Gms>NAVOer0g^a0{eZ<+Ac>2DMIgk61)EQCJF(RFGcb6x9w=e@f5D^q2*+Wt
z<bTsg{R|BMRiE@TF!0MeF#J~q@iRbjFCYB>|Np=0t$qfE3{ZA|c>&D71mc5Y{^bEM
z{}hN1it(2l!2ClXJ}9ZaTma_p0`WmHU;lCfSYQ)K02Jjf8^HWkAU-H?UlxG*i$Hu(
z48KeO^JjthpqPCb0On5u@j)^A(gDox0`WmH`O*N)Zvyc_5%^L8%&!9Rc^Md@7+wm1
z1&Tldpcs3}0On_b_@J12`QabP$4MYQD284>0P~|jd{D%^ya470f%u>Zd3gZL_v)XV
zB50lPw~&EhKgb_1n*aU(@6pZ5Q^dgF)A`)vIBQ!!C_ot)JbJgN=rA%ccpPWl3>EU}
zZBa3R2>pNxdGzuwEM#Cf-U3z&@~KC!=;Xr5{{^KL`WYBr<bjOu6`d{vGSF546s^6c
zGgTNEUa(|Mwi2>2`j!on;jIT#wl6@A6piX-VE8h@qt~{v4<yN21)?l}l!$tCvzqlz
zULhpI0TSso_3E2^Mo5OSck)XiIjueh28L4OP>*KYLp=-(3?7|NOL#rHZI|^jF#H$Y
z(8Iv+BKiOS|HoM0^lmm3)@0JlVPNn$_<&_6C_Nm_j%2OMVPNn)_=x#vwSZ_{4g-Tv
zFV9y{P#momkS)uZ+#wpqm^1mls03rq<R7B4Y9Q6UJl8xAK4IMjN<v3}3+R5%W?(q_
zTR;^`%YM$@Y%M0uA(+U(V0gfz*YsH;$Un9V5;v#IRWmatOn$5w%$P9QT&amMe)BG+
zau%-$P&B*(DYi9;WMFt9{qp~Rk6v4!NRSCH)`EGaXCSIPdTn_?YQDb!MWKj_N3ZFo
z2#`9DUfa(Ro2%75nHVD`AJw#Enierh%YZ2-VzZIfNhYS1p_BjW$T3X}oh+*>&lob<
zR#%QGGIVpSt`{rQU;oJ~jgB#e`fn~YZeV0O?K}Ci$pfZ?KAW$aK44)I_t?DMYByuO
zW0+&8V~A(xpU_~B&d(l=Z+3uEMQ4x73<d^<6pwBeRSj?qw5Y5ANpw4DcyzL;c(fki
zZ}AmnVDR{U!=uwh#lWMR$<^?i<6cmX1qn!abURsiXge8rG#_9D321n9J3Dx2J6nJR
zJo!D2`6U17bWy42@NNC?)A_@vvqnV$B<#_7(6jTnM`w#l2gu%o51BnW&w6zFs1$he
zyS(-2eE(wS<^TUZI^Tno9AIEzc(D?}KLO&;y!`*aXLpH;21o^{0`zG8=E?8!)uVHZ
z$_fSs2G4F26%hLZ$o%?bV;}zY|2#XtcyyMi2!I{)-;>|tpU1)H%%1!%e_wdRjCWC~
zcwq<Ocy=CtVGLn{EZG6FRvE-^u2E57@a#MW61~B|z))xBk?f+vVff9XdyfjpwLZN%
zM?opkr}Nc|iI@KWuLnth^!s$aGW_<U{?h;dV3tqk_ZLN%{{Q!Au2B(S;P0CWaxB=W
zx1gZ1@afJ`sc-~^LX3)oPp=8o(6>IFpFvUzKHViM1+Uro<ry41kNdX1t>g9Ryk+<e
z6n@f|{{MIEb!YTp{OV{~f28!jM|TT2U|c%?dUS(~>#k9;@aUeR0uCXl*@g%B<y$~5
zXXrfWaU4_-Gk{{2@ss7*(zhPnE-DHhy;D>YKxw4&o=@i<6$LOeM8(15c#Db#BLle5
z1F=EH4yc?3g;7L3NEOKc42;a6jNYPR0TStKQK<mY-7YE`AUzI@3=CjB$6Hi97$Hpq
zkaz%$-Pxit1FV5X1tea8q`m=(JpoC52S|N=w~LAd$R%KxcY{)LXNwA`;_>Khf%@S%
zDE)f8XuA%|o6Z~_-N6D!e+#%=5IFifCg%c2w=>Jp--0?SoyQKolV?13@n7=~M*hBU
zpa|;@_Hf(<$`VI!gg^1<ZU!akqrU|v|6pKX@ac|GQSfo<ZULu@qpu?_dGy+VqtS!e
zIlzP2x#H+=0oh4TpkjvSkVkTe3WsOsNsoh%nA0X4IT)>V%J9I^&ygpOexAJ6UZS4y
z=<mn_M}J2iJaVu)f`#G8!D>Mk29M?g93Ga3O7A*0A7b=m{OMwOt#+eNcZiCD<4%yP
ze7ajyKw+Q8FYf{>Y<ycU)hSOn`dfenYFxwtY{uPlY(B{70W$7f?J2Os4C@bo90p3P
z3Ji7T9^E})PkD52Q2}Sb?mb|8_~k)7P_YWKAdZFo=x;$5Hkd14cy#-yD0p<<14SS_
z!F0E%fO0iR44!bhp{WX<D!S`yR6IaR;Ax|?M<oIjzo3E)M0fXqT>y${aKZ$68#T!t
zZ&5*{rXEO2J`PC=pz0cAIy_aulIQUj6-er;2UV^uDxf3<3RsY0SW<g&AD)ft;z4<{
z!VsK=1;8Z(Gz&|3bmyoj_;kJoDJX#D=QLQB^8jUFpUzL941E8E-?{()9h;9c`Y^tB
zwfs~1$fvueUPS<ubtSr)J6lwyFfcH*9^mh*;$vX&NcK_T02S#BzMX$Oj=QJ`fZf-5
z&ZDzOWdkUU9(=~^!S8y^v-2yc1PlPB^$19QZhiqRxevau2MK3@geyFHGk$pXvY3NJ
z5?~Ij_vm&G@abFv$uZr<6`&A+#i37kH>imS5`smcPq#BTzTh#~>7v2`GF}1X2#M|x
zl?)G15$@CJqhjIHy&F_+`+)12&K5|-!Xwb9d-DZQ2Lcqe7NDrcSww7s6cNW;RB%Tm
z_IL!9d##uFTg-VF81{n#&~YD#e4&08R3y}>IMhRvEvVoCDfZ|NQSkt|2vn@Rcz6M%
zY!5h19eo*o37p-I{uT^T;dpWM45W}(VDRV-<?!izaP+r82q-LkR1!diLY)>O*FZ9K
z^C1p!W`7-b;OOrNryoau3p&;R@UT2sde@`d*~6#vJ1A;6KrsShNO*KRSAZA>pzzfI
zISQ0BUgSb8)?o112XZpF{Dl_PAigK076#{*<B-N8xCHLp0x6=8x2V8l?|40Cu?#Y!
zLe=m<njS2X964Amz`_D@JGgja2W5Ou%cG^wJ-VwM!0`!6BJj9@XD$W?5AA9TkQ{pc
z0@>lwe4Npvp5Nn#2iOR(H(>cQ`M*!E|9>B7A>iA3lE3o=Co~FQ<edh^3OJs?MGhkF
zYE&FR!32&wjq{)y`Tq+Auxe-=ehfbe>a`sGE$E}d@q+dA|9Wsm4Ua=mNtST*w}6id
zN1YnKJgCG&jl0iCafcCg2C%3D#UVKAEKs9v>nTW@K|~$2`~l@SP?(Suci5wjzvBiM
zs7`>?4S1r?MMcA>v;I3Mz(6(d3n`G}AjuzhuwscT<ZOu?^~W)bRgf27OgamRn68um
z|9jMTLt@7804RSQ{VmW98a)8T5yIt1e@C)_(;TRpJo-A~0x0c00B1vq7p@>v;1zU;
z3b;)GF17gAyQm0wFoVh`aO<uv9mKJ4HT)0O$^j~;4C;M)O^!i=%A=e0q7?(fk;ZBP
z7Mmjvf99Qdk$Lj}|D(SJPJqSij{c4S)zmI397li09XztITENGa<LKv@14kBC3%b~H
zcs3vB@U{F~y3L~*R9t&7=cp8b;$OfLQnd5;&iwWNe>byl=X0O#B`P3CcGg3rK&Ac!
z4h9BTi_WL>FR081cW^-rQ27YzDKUd$3#1XG5ZuZGm#w`mDi0VK7(kMs0vB5F`0%e!
zQL*r0P7VNj6;jEAN*a&uU{D){kH2pwXaLER@d~Ke3joD{U$4#sP>T^(-`g^Pg5;oc
zqzfy<(cd5vl=(q1q5_HuP$T6gsKp6VZ+PJ7!|>anEEE6<(WAcwOea_}FdRJ{$-;Ex
zU^PfPnj!T^A4l8(6@gcdzK%HQ(_N#I;M@7mv-1teF-Jegd_MYHutY@w>LZVCRzr|^
z;1bxg`3Q$+^Dz!kLUB>~0g16)$N&Ea+2r;E;zn?p405Tj<@3@5U~d{8I0`NXTo@RR
z9IUSvz~Zyd5hp-CdwBG9L<p#YZUJ|yJQ2lh=NfRG04o1{I=4WYyvS`1?P_pCUjbD1
zVJjDp{uW>WB@KgmkP2Ao2ufh!Qt^0;3al04)42v*`yKrqdGhG%n3G3;3vz(VOAmgJ
zAK><ZNAiD<&OIs;pepp>b7l{Imme=ek3ot^2M<u2ev1mI-3sckfKmr2m4H*r@fMXA
zpaCpcE9nCRqOA_fh2UD_c#8@LtnGcg1>6b&*PvIx`OXK_-)IJz=F?pq;L{xqN>{5v
znmzc}9|Xn!6mtd!59Vf&#YcY&=BP+`bngZwV4uzy6^^5)(<9${fZ73|W_7iPN9Q4E
zV;j`Mg%~?U<psFW-2=88WULS42jAAW{C!TWpbX=q65zr3$%pZVFXKgzUYqMjZ-zhd
z>@_)Vc)+_?<U|^%RUE_0Fj>#oRQp8fA#k7{y&b8_YR<rL^mc^8(cgj!KAp#3Ts;B`
zG3N@O?rIO8?&8UP#&Y#O-6@b>46OA73OZ0$>~>Mf07V?49z~SJpvDxao(9Q-+NGc}
z57g;yW`fmUg&?a?D(s^mgFzKGqR#c`uCAxBZpB@-rvC@I5>)E@wm#wS6l8=s?fPL*
zC_t;$dT5FU`M{$axo(|!98`qAf6;mP|9_uuXkIw(q5?`W46s&%Pv`p=8DR5VR182P
zGv1~Q3`Z|V!VBc1zXc>fohndrIr=!}Mx7HVb}Ekk3_k!WL<;Ib-4Jk7w>N;pqnG8B
z2je-9UXf2Q3y!{y<bae$5eGrp^*zWCkP!x;ehes)feK;^Mo{TzarC#Kga_jZkPN6Z
z=fUWsQUGpTS#ANhyZL+0GJz7|Vvy|1#$W&cPpAhCB|x$mC{OJ>gcLS7gQnsbBxnl2
zs$EnpJowjxyN|GN@?<>b!(5`m0V>wP;pN|J^4-<&KfgSKM>ngW38V%J`O9!bWMQ>{
z%U=d?Z6yFH-ew<~{Lxg3mFJ&l^YO`oW}@}Jme1>!yv+ai|3AEltWk-0@$MkV6^}q^
z5Ye=Nr!vC>;B4V?6jY>xdWwQBM?uXTMi0vib<bWV{{8=d0^IN5mOQvnnRyT+06@hD
zC=((lj^LvpyZ^uN1}k?_5damph)}2j6#<Y+->2J0#R8TBJQ<IJ;*P_Ixkkmnt5@U)
zs2XwskL9r(0u|=RK?6`ne@7lV`Z?l+r{%%gr!Q~Re+79QGQ0swq7g5?9l&rTDE%^k
zWD$+S4WKO89ik!ts%AX8p(Pl&%K{onfHq1z_}3o-r9_RRzXcY9QYqMt;MnwK{O8?k
z^06MQ1l-18c2SY=>=iixayG<4po)p-2&h+n=;-H2nPUh~fV%&`BMyVxGruE`fa)g?
z%M*3yK!ZOKFI_%^JmmnXAdVgk|9|vw1Sl0S_=3t2hu{DIzZCiX|NjKf&iCL^DA0J-
ziwFBb{)TkdUtHb~iGtH$!5Wo{7YD(tJ>Xvbi!D$NsDJ)qIhX?(B~0<qh723>w>|;2
zC!k|r9*svp*&sUBF~%|0G0rjmWzqis|I;QMJy{?A?da!-YoHi@S@HwqY*0f0>~alo
ze+DwT^CB5+3A9|nm3PQ!TY}ObD3OC|OHeewxP1uX^UM1n!2#+Dfa3UN(bxa~L0Jww
z=I?QwMFmugy<Gb0|Nj@$_rV?R7!Da$4Gs2ad;_ZYe7Y+&K>gPMpYBKtkIotu4xjEK
z4NzZKz@ytRr5<^l6;#O@9sq?D186vs-{psocBFxC>o*^Mmk+LnPe8gnJbS$vL3N)7
zXlOLTqnisvCwO%GJ9udOTX-}dWc1|sIMd0a;*osPqxB?6R>Fth^P%U#C(J(lE)RS<
zpT7{_2lC`|kItj@FL+>V7Zr;a|Mx;ATvQBRd_r&(UOa_xJUfrRxCLQ)bRK?jelIAL
zd{iVrV_F)Z;?TjP^<-U$Pp`^ZP;i2VE#8A&>(luN?35QaVBL>EGWTET*MkkHQE~9;
zd}R3Tg#whT0O1OPxp3pVOH?>~Izf%%&K~f%B`7{BV66{O!U2_s9^Gyn;3$Q)M!LfU
zd^+cVTLGYwTmU2q8si2<Ah`JiD(OKSkOXL`7Gw&jGkyWA6x2!viGXB5IR&Kh=x+hO
zqpxE^E^r{Fx#N&VG|1ne>KIfHfM`%A@ol}$-#Ho90R<Px5}@D^sDIJF2NHx}yTA=I
z25^}bq9Wnhc?4tts0ITM9DrK493H(jpo+&wMFKQV36cXF2X6Xwf(A`HYg7b$I{&>W
z1E=KYFZlNS|L@rC&H^fkK-Fc4iicz85yyjX<Q*B0)H`1M10J=!0v;a^_i)?|_TPo@
zx1QbYATNLlq4}Vc37QjGtjoY~^l`*(*f^CZvwMIivwMX{ugM{gWDymQ&hMb%z#|8P
zwNArEDI?Du{Tz4t=;sJfa;#GWb<IFU8arq>;As8t$b(0JM}V8qf-bC}N(4M`dC#%=
zFryFScUQ|hwVS{LFW?X`JOD~wil8<b3mez~aF;5Ag&8!U$PB8?K^6G}@EF-QN6SC8
zXF<*~JOGM#1CY}|`lUen!B)oALv+H2BA<h*5KvtoqT=DX2NadyVTJAx(4+!1(Sed7
zXb=t5tOV5x0-%_b044qctO*j-R04%0C;?4?CCcNFE&*s1oB<>OPlKIvAft1~AtSV~
z)D26QAj3c@3Zw^H$^?mn+P$Dy1m|;bw1Y%p`2&=)RXjlEbTBZ~gDppQ6-YBk3&>HR
zF!1TE{0}OTI`@FH9e-yk8??a(Nqdfaz&bj8R3u&;*$GL{65y}^>4sPco*QvdDS$P;
zKrzCxi;-dS$0%`7Clc%sSpEW84)O&^yGOU1gNL>oc>2hL-{X)+GSYyLPjCExP#ysV
ziEry|SUv|8IUf9-e_4>783BqXj~7B<2ZFj~&>#kd0t47NAd^5j3KXjvFHY?N6^7qm
z9NGc#k;l=K;opvex|0&1t`JD!(bEFwz}}0P03H%b0EY+EqArlWWDylm>hn?Y0NJ1b
zQ4C6m;6b7UP;yk_uLq67gEAU8QAQp(3LYgw8X<q|)A=09XwWPRs6+I^W5<6`#DL=1
z10LSsV#7tH08unSZ3PXkz|tRVnhn%bWbk8v&#}M?Mo<jEqQ62Nnu!iRV&<0zO))_G
zyx^uE3kxjUJp*}>V;99yS`Lp=P~RAo&3s$G@prbuIzUJv@xpXFC~{C@6V#q8saN4Z
zjm@{aK(YD!#nWvd4L&L!kg<k6;M@iYAlRUTU<@d?`=|&&EPAnN8+a5aMn&W3<4DsZ
z8XyZGvLz}W9<3+&`;1r^7<`gVRD3(%gN!};TcAEfCBmzh1w8o-DjGEWdO^cwFObFr
zA>~)(A<zhhr{xjQXlRIvf#CsA2NV>e8qgLwa>Do>$-)E*Zg4a9cEkhF$j5!(&Uav!
z3LL!{{`Tl^!5Ec@qxHW9Q&bWlKJw^hy{QhWR0MTYU;`d*(D9B4u-_tC4}cX!f?C_Z
z1$0!7UXFO-Yx%4c)cORqazBQjJo;MzWbkhRxVa#Qv4aM!;6C{raRM|>^30>V9yH1Y
z?TCTrBtdQ94WKr#N9PGo#&h7}{OIXO0f=uR4tZK0D}8_TainZM$R)QUo`6c?fTO<!
zK?XgJhyi7w8Wo45k7KSJ{T#`0^mW8l=x_#)$_wGGpqb(OFSxgYBKiFb=B@RhSguh~
zc=2-!gk|vJJ(Okf;wgyL`ToV-E#QddQ31{Bsc?XX;y~kKpxO?pC0oSAz~Iq%1k{>C
zv|taz`#0dAj{tRWAmz&nkm227m`M^nMJDV3IrjOB2(Shp6%W|>7jl|}rU6imdHl%1
z>dEyH;`PuJbFB9IOQ*m8|D(Dj0p^m37rbDD!Cgg2WdbQqK`!fT0jGONL3eaJ$RXcf
z?A;6ssuC3rP(cSNwIU$p9jKgxm3cn(-7zX2M}G^<1H~pd=fG!!Kq&&+D@1fwen%dL
z^tcY!-ha94*Z=>XooB#}TyV&FfEuQtkduJ5PQ$^@05x!qyQoxv+zhtF16G5u8y*19
zk4?_;lC1Z%JYM_kW#Lbd%WG63z~+KY1DD4cU@J60)qwygm_T9U37X$O_>B3*olPK3
zpezip2tYYIL?r=Sf<i*b@BpZ3D*(!mpvIGrN(R_58lZgqIg-i7hw;<V-w|v+jQ_w>
z8%IEWOi#<hwU1wlef$3(l$DR(7Pxctw!oRAw*}rDy)AI$CDY&k|6dqv`u~6OgoJGA
zgUSpHM}LO1a;h*e9EMenM}G^jzEGYlmgvE#yg5G6l94H7!{k{>c1$80Cf`a*Ve(r)
z*&+Ej)6q?nIaBhPLN-mVN~vO0nEW**gDHRQ<oHx+rl7TxYf^)l<kn8UkZPs&XwCos
zhL^yT11>5F9<2vTr9FC0g=9eE{2slw#&Qe{FYMM#Hb@I}_g?e=Kd20Q`RL#O|E`Ar
zJ$g;&NrM)&`~oj%*(1xq@Z$LD|NmbItpTlS5nU_`F~oG1G-$%w_Ap3+*_zEK(&jQU
zy_cDso}t3UC=2qVKJQ7H$x}0?F=@(7*3DGmvX%jhvNFqTj?O&C$S65kD?5tGV)^8@
zY&jJTR#1tm4(jq~cyxxSRDdE?0cBZ=#N-p%Dyj_2L5cT0sEp-!@p&1T16rl=0>(V{
z;_kA^f;n2YpfNXo&lesCpMu6XUxQXe%v%O3H^7<`UQC4O096~Hn(sv;jDPe+*|N!X
zIbykK5|A-+P-D!u^+}zCM>mTKXo&}C_@Y+?Jcs|{=2B3y@&PT*Sfj!KTHA2EMdbt&
zbQwf9>qH4qXubkfuosvR74HKWyIV#DRF4WwzF#fH<){EE1~WW5=S==pEydH>qVfY&
z8g;v<Busu?q*~wEqH+Kv0PZ+|o2s2HDmOqPoh>R7Aj^B9DnOHZpcN0G89ipmI*J!9
ztN;J^XuVyhf$#^1Pp=BtA3mLrUL08h3Qlle?fr}GOF;gTQ2{kRKz90czVe(blqci?
zwhlalkK)l16>ygZ?#IqKP`85ya60FxSb$vD*#ZfoZWon`$(O5Dxk1Cl%#g5}e7DYI
z@`FMa!~YjNI&XONM*jybKI-G<2jx;06;PJ+IQpxO*GUYN@F1ZN89bTHlCLES4yfZT
zDhLO(sN9%*p<3A*RHS%x+o*sGi=#g^eN+r!g~N+Ui~s*W`cnfW1oaJAb>|^avN*sr
zd11abBgf>k`Qnq`7b#9o&Xt-hUm#xZ0P@T67O<ZnP2nE!DyIJ>DjGiEr7YbgppkA+
zg;ew%6pnkq?mGHgutf#b=mpnv2`|bPf$|jVei2X-d4Ie`MS+=tAq_N=0*lTapk_R1
zF(b%50?Z5yFO(Mj*ZU9Z3xdqq1?s8#fK~^AS0aI0b1#e*{r?XdOU(iq-py()GWl^q
zQavaN5iLYeV1gS=;K+np268z#{ef5Fz>=^}ugYCWn7|?wJZT1v(C!iy1MglLbx;Qh
zl=yu*UxF$b)TD5{MP>5!YDJXDeqCfRS)oX@9^_<@Z$O<dkRL&P6A&NN9ROwW0-w$o
zzMX%4JHLa&EWoq#*m1}zCeYF(aI*rOWqdk6fwCCP)gU%JC4km}b^CyZP+L?G`3{sW
zj|a4<ykLT450LHPL<Guhpmk&56x$6=*B{`jdLT&{M^>E7TH{v#f)iBGG~1{!)baUr
zKJw{&=hONBMaq0oN?=g|P0@g&1LWkLpr8e3<n9_!gP{kSn|mO+8EiXxUIQ)cserq-
za}UHl-C+NK(?2r<2gBrx)fz&O6bcFp4p=68U1T8Z`2Ugz$Z`CApw-!s`7KD(5!B?U
zDwY$z2b!Y;MU7AA572NJzkI{wHN}?oM-EoUfIBU)&byE0gVLv<8212=_`G}ZbS}h`
z9?fq!d^<lPEjj4+Q4#RyTmvo}K~4fi53IWjN`|1M4DuBy8G(Wy6jLB|p!fmH9dA)k
zfk(q}NYMz%eJv{B9NF13S))Xram(a@63u#~yx5}xPF3L5M2L*Z_X(8wY|Hr>Kucjk
zIj<6&^Y%c3em^4T?NI^cJRfipM9RMzFHX+^d3FnUwhx*WzrT1k2bBNnzrT<J2ZiZ-
z(CQq}+!-jU!8U9GFSv_{f<}TGJE(*UQQ?5J!5qL{#R%}g$!EwEXar~(l^~=m01ENv
zpc;b%6yl(1mKUZVXY{hJ;{z2V-$Ao(%nZ&9ph0C&c!Hwd1Em}Rh42eOkf!b};Gh97
zh6SxeJTn{Q_7|Y>S$=t!$$Vwf3I#Ba9sMoX1D-a1@d2#r`wI<_PbYhpDKL6W&MMo@
zx}K4N!Dq5oxx8>OFSxG1$=^O5BnDbFI60<V+sYIq{yWA+CE)1K5LR~(@2CEYGEh^?
z_9Hi_&ms__5^?mW{xAM^bI=53H;W1=%^YBwe6U>D8<b2SrBvq@NGx@O3L{V<1WGCJ
z%nz?OYd}d3nhUx?6G`w2V)DgmX%(dE2wdhNr!G*ifs*HBP)XHSp<e&|2e|jJ2ON36
zo&P|czu$s3DldaTg*Lb(0xjfw0P@Dm@JmO3>fhvV%lrol&=0<y|3I0fz_$}r%7Z3O
zyFmkRoku--S>w4OYmq<`vY=(<pg03{f<1a!13^-qkmZS!n=3`@1we%!Xc21E3{cAB
zQSs>nw=Cd32bW!NPlD1_j0$+_5tN(2xd>dOfRZCpfd(n*Km`c6n8Z@FLCSt1klR3!
z019nLU`_s8?Pu^JcpkWA!QbZvsu`1~sH|ZCHxs*kR1!e5a-aqlXcPq6y_tNVN}6%O
z<U3Vz!WT+EgBG+2_<%zOl+XGn^HxjMujPcay+F}`nR<_Yjy!PmL-@&~pS8SMj{c7M
zc@z{AM}G@KV&}#GY5)K40_|M@)s4NTRh$eAFIG?g{~wx`--8y)gND;yuuTWW&-)j4
z(?Q|s?63=zNj<ckQ!Ky(qp)=spk_X3FbcHpViIV2xw8f|2o6qOIVuXBM?tf3y{tbu
zKsgW;#UPJN1ZVaZ6@(XH(F@LiG1DeT)yg^LOar;#8O&ubV!-_Opw$X5CQSSP-|}v4
z4`?J<0=#k%oQgbpS<5)Ur4}e9O_{v2R<0hT4J7ID;?q=cBlpEC5Z%28yl&*hT`;Fc
zCE>+oD9ho+$*CZlzrQ#LrvJa#36<e^u^!B7QITK;Ee(gy@$fN%b0~|7N8=GtgB_`R
zTtC61m-Qq&$otR&h*0^X1GWRL{CPeF6b&UR@InWa&_GQAP_(?z1UvNmi^Qq_|AUIC
zm%gBWAh@9^49P>FIt7%<LBpNk#-BGxqetuQN*(a1@{5}w3s72mkoNM+jE|s(VEBQT
zJb(ZH2Q@R#voSEd{PW@e{}-iG{{IJM@5%Y~vW)DLr`9XkoB%gQ6u?u@$6Zt;z_#?}
zs2KS4f)Yk=iHgF@E(Qh$aK{xS3GSXMc=WOcvrT?ouT<{=%EY@S!&BrxQ1|2(f)8H6
zeQGi!CYh#yOCA;#P@xGLi#_^NpVgO*fuSrIJdyz_K7CXaj{c4~2wJhi0cwPG-aGm`
zl63>9TLm6VbWzcGq4nnfe^70H^mo00kBY|8pZXX0+uZ(x3>WlKk$CC&ALLO6k6zY^
zprsO>?~l8vICz0mF=%2MsSW|n7yX^||NkyX)AYSZuW1h}Xo)_=#n8n3?gb}=|H2VO
zceki`FhLU3WbQ^OtL>9O)`6;7hR^(RzZ}DTdTpcS85kVHJUg#BhI)2>bqsOr{1fWY
ztLr1rz!2=C`4`lxcyWH>=J>{ojGHESFilRFw2lk5H}HY|WP{1F(M1{gy1AKonYnr?
zB@C%KsmUcp`FRQ{nQ3XMMX7nosSLTL3VHeEpjcV)W%}&$$3N})tA1=>I(3!BteGIN
z)84Ab;89<ts~t!U1B0hSwa7a6hgB{1lUGcZVYHcibh3FpXd9pu0|P@u&;S3RrM4@2
z|Nno%#=tP6@Be?$fYOZq|NlWtKYmR3{~xqoCt~vd|3?@Z7&4~(|Ifk1z_4Q4|Nj9@
z3=9#o{{LS9;?Mg3e-9G_L&co`|3L-DiN*i_M}P{`<^TV$U}j+WvFiW-AIuC49H7w>
z7KU(!6+8a_pTWYw5OMDR{}(I_3_EW9|8D_`w%h;z`>--F+_?S!e*!B5!;c65|AUgj
zjfemLUtncm*zxrLe+4!MhKe8m|4(6KV6gc0|Gxk`1H+16|Nlp@PkuB-O=?E(|Nr_R
zzq_ytfG8CP1_l!b1_p=j$x2hzm?HWndrkFY%IKdwYpNAf&&0{sraCc2Oqwh;&5bo+
z%K!hXC+ALAncO)|U2enl|Ns3!COGm5G&2=&@o_kEgHrwjMh1o_(<h&s=EamTbF$QQ
zKc<WYle4D#G0j*!dDnD5rXR~DbItH$+OckO)C@m1jxGQHuN0juH(h-4+*w@uUVI;z
znRkeA@mV<X894H3IPobs@ku!G2{`d_xO0OFX9fla9u@|MjFbQWC-6*`o30#sfQk7&
zCsrMfARQoI#;`CjR9yZ4zZ2wM1_tnk(gGF+hA&8BJ`4;DEi4QS1=s%nmtvhPH(h43
z(JU>kD=Z8Q4!8gRzYfw0_A1kME<O$q?jSx6kd`l?y+(Kb|G&Wu^@Ki`3o8S|g!}*h
z?-l}?65-5e(ai3{_kfEziiPh5R|ek)E~c$~FSr~*LXVmGUT`t(1BuOriK((<^^FsE
z00RTV4ps(+J0Ji5-^mEIk85)FY)v*EHU@?rZIfrsR$?<@V_;DDJo(UUC(ar+1_qDM
z|NpaYe#mXjG&y7L4x@@XkVhC8Bq0<7Y<H_Pln*LlK}=Bk2hjy{CkM@|XI!xP(mZ`e
z#<`OnnZzbbEZE4n0L*%@K$5Xx^REScjP)>sL7HHgnSmd%yB8+Nz`(%FAix0Idy6E@
zz|0`XumaS!MHU3B7lJz!S&ErKm;t)I7Bu^hkYZ*KVSw%VMTmkd7G;2Ki$xU|V_*O^
zWRL~H7K<|+z$z}mpa9)>3(+#UZE=&-0_eVC8HgBzF?a*41cSrE$x2JgC7~KX{Th&h
zYOpNG2a`7~(Px}H`O%U-CX?>TMN9Pr(L-n3<c&)eJfS;@L7G8W6>O3aEC@hCU~C2z
zhs7;e5Y#3GFWv+Bbh6a4IKhTYh{SHtenikr^<+n8@yQF8aT-D$3Xx`DU}a@sP+-6e
ziC2sa41(Y&o_uSWo|HokL<4wR3IhX!4if`|5CboRK>TE-<$6*FQXtacoW#JuP{Rb$
zEza;@!Q`yv<*+@#5H$=8!psoO3PBK&1W@~mfq}tdvf>JHPk*TF1v7~34A3S~1_p*{
zYzz!SOneOJiLnT(zQ6{eekF8|>*UFeE5r?#L1hm>cSFZOTjG16;tg>SU9jTg0a(|6
zsH}iCM0OQu<01nC11BrQ<q6PW1(&v<;Fv79Qk>Cs^1_wkl8Gz~41!GJ3=(J|oePp>
zte;%CQk-!Pn0*XHdfsAzm^ncmVplQLE@oCxVTKtqhKwL_2?mcv5cM!pU^3$>aYn<*
zhO5LGBS55N4cG(;43~6(W!FIc!T>EbQlSRznLKfoIO9FAn!l`|c;shz0IuyB7(kON
zAZ1FEKdusIbOe!-fuML{;AJ?V2a;i6;Dr`51(Ow5iyQVpWecFTKqHf34pe*r)Y=<R
zU#|z5$#?)v-Upe$_zpx)-msd5QE0LwldzCEJ0!do#6#?pV%^NOMu&NGz<Ohr&5Jg>
zF>ZdoMUk0_HFq-C4m~bdYN@E3?6kvN0O}7|)8oYC&K>&j3<=AacP5|Qp)Uk0tYGv6
zw#i&O^*N#WG^2O2(@t|3ub_8w6M|RKJNeK~b57784FdzijsD3@yUaOZCC!A%P6*zC
z36q<4nR7x@ip0dphj!VsfchVk3z>x{2e7hCmf5Y!0TN(fVAwF(bGJF?gV_*oM^5hC
RZO&!67{X&_oXoXV4ggmQ&x!y5

delta 18163
zcmZp;&h+9t(*zBs6<a21-K*bP$N&Z~nn3`}Wng4rU|<0;Cx8eB1_m@*f(IfFqmlKo
zF)%P>%!SBjtVN|aFhlqcz=ktSZf6t{get9oDy>*Mc@<-v;*ZG?mFV;vh%m#O$(l?7
zFue$CCU-G;i4+->1*Xqi@x{2Bd%~T3*%kk9ZROqkf=Q6A9%=^%)Q|vS2%`W^oB=9+
zVgp1TrXI%sfhPU{7CakK)h9sB;ed+Q!wf2rgDAYQ9zuh{3=$*|{sSa&5r`0&e1Rk`
z1{Q%3ACSbwAwppC2a>o114BJn5JX5oeZbDZ0FG>s7zisMi8Fx&pjZP*92$pU83QD7
zR<H<!us{+Ac@834&%oe-BmuG<BEi7mfg}zKC6HVIk~k+w0E#1!#JQkiASwY#oEszn
z#Tg)RkpDrJLQ^kTssJRxz`y_!=LL&EhzcZeK8O&QY(Ns{2a7<64kU4CE&@wTKoSoV
z1Peik8AuXB5Fs$R07)F0e!$`@ki<p7A`oK3hRvtAomlGUGcb6x9w=e@f5D^q2*+Wt
z<bTsg^BEZat3H{}z`!r>!0=xc#Loc9y?pTh|NsB0x8^f2WPq~!%L`!sB@iDJ^Dhs8
z`KLg9P>jFa0OlV8@j*%T<pMB&7l;pv`TCa=zyg~<0-z{=*#PFR0`WnC`?3JcUj*WV
zV)$hOm_G}|2gU5m05E?Nh!2X<mkwZl7l;pv$(IITeiMifioll&V15;d&&$9N#qd%9
zEKmdz0L9o#1~5Mh#0SOH%MbrRK28GhK{5350hk{J;)5dQ<pnT52*d|P$jbv@zSsQ8
zDT3Asx=joW`$7JA(fsfKe~)ffn<fSZpU&qV$64Fvg94O+!J~JJiVh<KgU4~!%}^nq
z-WC-Dh|mwHkVh}?#YP5(<1JvdAfI~liXLv9{9jO7VLk)Hi#(9=y`tsrAOmd$K+)Q3
z>gmnE@WLc_vXzjHky<`ThIf4)h_ZbFa-?X~Tn2_O6Fhot8|Q&US*t*l<&P3kk8W18
zxsz82$#8%~dQH9NO+F(e!#H>HOCdR}c?=8;rN*Hi&9;Z;FfcH9bUrQN_2{-;HkX0n
zzvzZJ3=A)l|NsAgjP=dj&4$97OnNm83?2s`u<Qh-hojk%tW`A(44wxcF(0iK5Us0W
zVDRbX`3ee(qtybkWi^vKM8g<sCf^s8V62(^LsV7`q`H^qn&-hMth+!-=;&_&-Otqw
z3`c(ps6uJk&()i)#iTg|3mF&;4|w#NJ}U(I$96&C=5)DgX2ycaj}?O%3nrT@H8JLI
z-lbH|;*|l4hF2iPwg#CD3@@Z#{{QdMYwME<GU3HqFwgW1M3qOcEe}Y|_ZOfj6jAZ$
zHQkf}Qs>cY`#EECwYn!0W5(p8nwCuKGA3ymFtucCHqtuD#B?xq@?RY}rj@CaWp(8l
zQzqN$$}ttDZjROUVr3Fdn7q>H7*lTi=0f8JMyAKHlRuk0V7eKz`Ksvy7AF0O&D*VZ
zGuC@_UT_R^40Q|%_2~Q??9uq<1}If@_NahzwMVy?hDYZXl@kmM3@IMnCaN0X*w~_S
zf{}s2@S8{LN&c25eozFR_xOI_qq9bZ!=u~6aUUq(dUX1zXn1r7Sa@g$7=WY%Ji7OQ
zP4eix;L&`5(WBEvMWWuL+ugxK+ug$OKUj?4<G4@q3y&Edtp`Bz79RYbKRi4CcpiMh
z?91=+%%}7Di=CJM{|BY-&K8vm3=9k}RwDQ}K>V4PL17rH0kh7t^SDR1FM~(xZBKre
zdoJDPV9h)#FA6UI|6dR0w5U7)=}dr`bo50ijOn6M@xlXSVsnj(1%pR7$WxszDhhQ$
zAoD>E@JzmH_|2zR<tiwJ`E)*dvGWqh0W~TWKAn#Yzr9!m<tjk9v%%awU>kip>p?0y
z!9pODyFpQJau;Nh;Q`myZ*{^GKp}M;RPZr)Fn)5hJXZS7quWKr!lQSJN&`5m)~Glz
zGBEgbL!!x}cMq8F(HWxR;c>i0#e<OnT;71#0gMa`pd!wrvqhx>q<IIZ=3!uD21Q|u
zN(9K9&Tk;u86eqi7ZnMRwgi~A<1H$nG8r7b$6Hhi7$KE3NW%`O1_6+U2}l|iAh9<f
zX#hLEyWT~G!=tlDWd%rk=N^>?Z~|*l>0o4Fcp-e_|9{VJHxAG4FoC1L1zau&9Q_@W
za)HCK+l}SuZ^0B5j?QBT-^nu`yZEp92P1zUBgp;TF)9|0J3&r7dL#UaM|Tf6&G~ds
z0jHLuzXfDNf*BY*x*Iz{RylSDdK~>NF!2Y7(`$1S6!IR-jsZtsM_lq?cC7H|<vHY;
zeAVONBWBOe6a4ZZ_c0tf7_D{C@W9dEkq3_co~-L2QO^KjM;<(KusVW;>BzxqK^7*D
z<^voamWN93IyN6<^kn?uVtKB1qfd8;iiIO6v>1H4eN-wY7#{F#y;P^2#xLIjvK&<C
zurNZ5i#Px_juC7eBh0vaj-VvL_|wJmSS`qHh6g<A4?be{04LBoEs&2v$<m{{M5O{$
zw;ug105UL+g$-;DJItKt9^Ex66`<&XCz0+J6;KueiNTVJPj`um1xNs%2)dyO3!V}>
z_tZm@A}k4ZBjOR1G(C2JJO)Z+KAl@2>8=x!*g)pO5)nua6f7VyQ0#;1NDv#ImO5Kh
zz^U?hiwbg@0;QaKP;`S7!xGqw%*+4(`?h|ovjhbP$WNeLdcdQ53pn3_qLW{q;mE;i
z0T%Y7zav0V2?`@Y7B*PuJoe~r0jHd9hRzn10}Kodtq1t~zH)++4x|{{1@hH#7ZnXg
zQ2P~>|3LP5beCIz%e2lCl>$%^mf@lK0a~;j>O9oxq7rcMl|17ij~SiUF8=HFW9;0b
za)g0_p*QS*XOGGZP%=zD=aKvglqtX|wA<eUR0M-s@*dsg0Uo^~;Nrofdp}sZGejlA
zqrO{41*G8vNPFiN6>w^1{O8fRN96?rsQOd+*u4ep36IVg6$6i6k-s3vgE$Tzol{g+
zfE-ZA*tth#0my{z8WjbPUY6tFe9geX;ML3W7NXuoC4pbwg@M5%c@NmF-99P-9-Tkx
zLFKzd=XuTJoi!=~mIwKJK<#@_;`0EPA>BSI3f(a(0zSPuM}PeP|NrQ3fslhgJQxpq
z9DK|SF|eC;r9Y^~6%<i9vanjfMVH~o!=HI4UMQafD>t1DmK8X0^tXVB%8^FzNEZDg
zjn$DXdLXjiMTf()`6$QH-;oE8{*H7};W+v_=7Nvq>(Uia2b8ENH0P)YFdqFa7;@|f
zI6gqtPq&MTK=Y6P{H>w?L1ng&iUKIwK>^!bqY}Uj$~p%@BB0Pe`dg6YV4a>v_Y@US
zf(E5khHgmY)$_}P@~zVkkLFh#M}G^t90Xfod93~=$bJV%-aPtSK!U%|oCB1ELR2aY
z4;+0NehE~gD1h>60<>@e*~xVDcO=NCEZ|Ik^miQGR}mLHERU7m2Ag{Hb;M=EZywzt
zD)kAVysU8ax1hia(=-47H-l0o$k2#x7ZnW;#u60=ALb?CXtr(v$1#7aEvPxyU80fz
z_K0WaTMtND2DxsJ3MeJ{bna0Bc^u?pP>S;DUINY$9<ARzIzM~xyB-525Kxk>2f4?>
zqxl86Qg`q;_`)8P^b<UKS@b|^z^SUYMFrHQ1Mxut3<_3I;!M+nm3$uj>p`yV_E9nL
z1QlK$j2GYmC&}OUo{fRQgYm=BkKreMdTow_%3lH19li_<M-J9UN3yUVIanR(!pZ>3
zz|b`CIpTze<>AsBAP;~%0ts=1g`oDGhc5%e(c6*WG8?NQuR%e6xb*VT&k?6V8AigV
z^Q}+kZ%}>*dF$xI@Y_f0KgYZ}`dctZMIa59njiXrQ?rbUNApn*kLJT1M}G^r{(z+E
z(9@ta9qIf7>dxN+GAc(ON8IqVd{=tJqdP`L0^-sMkOB(RV^B9eIQlpuMn%A*+ouKG
zT-yiAq@MdhB&4)m11^_AMR*1%GC{?sk2a`WIu5B}G0Ip_**aO#OQarDs)7ntP~qs=
ze2fvfaP{H$_;&Pn<R#B!a0z1ps=5z8XFmEk=GM{Qf)bFl(AlB_N?0Jw;KA?m<He_w
zpd1#W;^EWz4CELIkhbnUD&Vxx2~8WvTU0ne?Ul(#4K?bIL%Jv+XS^sp3)1bQqT$&M
zX{>p6Ly9fm?i!T<-|iTd3J?DE2SG9F<_$^`)dFk4$<4ES3pl}pqesHCdk>@;a69s?
zXLk<BF)}J1jE6v`f*j)04RHvlx!?(MiBESAq@TpU{-6)z4@i3C@0$i{F){w~VZ7nX
zc=71X@FyO<HrG6QOO7%a9`NoJK{OPq1!7nkChHlS*54>S2nrd4qqief*LZ;vT?8n@
zD)@9Bf01$m)X;HJsQ|kP>?KgtF9G|&qmxAi+}e@wfJZIJBnMEOz^Z9bQ4eaqfyCf7
zHF9&U9@^yVc2P+HDFRiypoG$TfWIY$k%3`9IGFz&|NkE(3u+O7WIa2<c^p)WgA=xA
zw~LAc*cBe&#0E-cAoZ|H+EW|UH0|7iSyO{z6PBnzIzTlw$Odpx%<pmBBl(9%uls*c
z=>SSIpyEToqxArPrzG5Q2_VOHPXV_PeLA0kOaV0%z@~K9sCa<%f?Aa?#7~2I2;X1u
zAOHU!luto|^&ZDrR6t6h#jH=~_ZJ_IffALAiUBCuy!HUaNhFU7s3?s%`dfg(16qvw
zsAwGh9dq#LZ$SnR#vdS;8Cdp!Yg7K7Rz?N}&+Zx(2T+M+0ZT@n-65cQu|_2WB<JDL
zS)yW6@bb_vkaduX^TkxK-61OASi}<Ch6niN8z!%bm8|!$yioV-Wg)1zKKPs&6gZF!
zWe+xih~W5k5)vHmj)LqjQPBVuV*#H0>$iY2A1s7CIzNC?I5Ue1EU0{XYmPG18@~1E
zeC^Q<D%>V`bhFNMhh#62d!Quha+d*ArK$)(N;aNjpps1>L<Q7dzj5@pK#GdQ(btiO
zjvS1R<hk$Be3<d*XHcU#?(osy5jvn`t)t>$dAD{ssMHI184B`bHz-0trRw7dklFPN
zzMwSv>CgZFFAx3s{~uDq`gB9;CZEo?p!iaFAq@6hjY`D}P7tf}{R>7A-3{qzz4&qj
z6bAoayg332gFWD8?u+|iP7Am<o8qC(0`A_HM!-6D9*swCfZBG^v5qm0v5s+$@h{gM
z`Trl3-9a_?(Z>;YUJ8Tb)J4VN=*952M;}Mr1(kp=g?@mtk&6lksDK6cMP8PD`2YV!
z{E`3v(-8F!sDb1d?$LPy+5-o5zav;c!v-#(wquQog-2&ey$T1oKb_*yEdzF-N9P<B
z4;D~09jW2b%?avfM|gC5Sa`JF_TYE9=b`Ok;L+>I2oe_nCDsIwZYK>;f12Orm@6nU
zT|liICXiGHxTUJ?ZUM@Lp8Ot1J6%*9;GJ&+AAZjlp7jTxGW+nmyz}XN|Kj3dP$<6l
z>^%137=#Hj#pA_J2*<PY$cuFlCP+Ns#X=b8=!<C(rf28z7u|>d{|6PM2H@^A*o_)>
z{-7ew!0?+-=XbEnK}|(acJ-<6eD&hjA&|YVK(hZ|ygdYp!WtC~pUzi?-(K8<aswdT
z^I)zIs8$4ZTsqgNd;nE6$6Hi>fGQSHGw*nd3Ih|gBYeC?MF7U`cH;n*DH2T3UTt?6
zs54uyq5u~tfpk_gm_TD`Eh-K$oyS{L0-$V2KNVz71WeoU7L^2W2OiWS5a5IKVZq*o
z7VDiYDjFaoLCs2dXPg5h(%GU??*U@8p5*T=;er*@J|IsXZ&9&;8td74095!HFo8x>
zLCteeC)lI61?+d9P8SsppWZbpAZd@z<DjUC@acTz(RuHM^9fK136k$T{^Huf|Npz4
zSy21b??FSL9^JtnjyorxvzM$7Q857ZjY0kBZfDRC7_86i9N@w14DK&`B!{Rtcy^wE
z_p+k3P8uFK`Z@CC(a({mkA99jb@X$@Nl>p>9h3t=#Rv;(2O82-h4;mfI?(l?9<+~&
z0i*{FDGNY7=#zMR(4e-VAPc<bi_(JzyVmdk$gz+_S*HZ8BESv<IVBF<xD;ez@o27x
zbm^ahqR+#p^A)Hi0}j`2P`?T}mBISW@TAsVqY?p17x47f35g#>_ZXD3!3m+WMFkxF
z$6Hk3squIVC?VE^n$Qd&g|I{jN&}!`24plmL3XyNfRhbU!UV~qCrNOfe;$+_TW{By
z!$P9lN5uowc?bZ73Aow^hXf0#Em(i_cO+&{8B|h1$MHZViAU!la8mZ*_qfn`p)*7!
z;ou8-#tRpJbh@ZyG#_BJJjCAv>cskVu7TD-IVv8Vmpzg{XdXWLTY%Tj5mW}2s8o1%
ze(Q`;Q2+%pXy7GHkG~!|N&*Uh4ZmKSub}Q^cZ`YxsM8qo;fD|7XHcJV0*LEk3#vcB
zb&;Tt3df6#{m7Nmk;ZBP78_9S(Z!a*v-voOkLB;9zXg0$K)uMr;Fh@rs25qf1#DMy
zj*0@~k%jfuf*~J%bi1f1fLg`iu4B=EP#!8#sQ}eUAZPk?L%ijY?4ptZ3P_FS8kGca
zXYyFD!+(z%pq}KhIunq50;vC|(OsgF;n6E{5L~%K67|vF0@%8f^+=JV0g4}3nuYci
zj{X+lfs7t~jynLd2i$`MjY8e`=w&UiX8?`wgZeGippHIhNCiaqvL@Su@-n1O1r;UW
z5u_aqj0}vRDGiVar~~EESq~cGJq{Te0ZkDcZ&5h_N=2aIQ|Y~sg4qH*GJ;&Sc<_6G
zMm`ikLnxs395_ThI&)M6d^(@)2I&U{gNILNiAsWJ=QmJe!@{RCMJ2<hvj&t3j(c{V
znc&gOy3wv4lrr8Qhm^f_Ah&lyib0PZprD627L+Al+}Q(iPpE)T=X;NCR}K%xTOPe8
z7j}ULKwo^=1FBK}fjOX>?E#qo4bHy+=0Ah;4}i^bQK<)+21+fUc;x^~do~~A@N7N;
z>fpKD_yKO?fCmOa?Hs`ta1Z)L6WF--FLv$u|KBo3MS;H;G-(9tm?ijhr+`ukC@jF^
zaowO0?qxk;3y&(0$3QUv>1l&=7R;ZZNO@5Pb(;Ywi^Gc><eCk$9066W3?97^|M8SQ
zt+)9*7c#@@oC~`_ZJI<#Jqpqet#m-uD}3-06e`f@hsF^ov_bBI)sG-HdNBkt&7<4R
z0bEl}PKcDO_vnoWcLl(OSmzuSaA5~>Pc9S8J-Q(Gbcd)!fT9<vE(ALP6kOmm@O&RA
zb3!Y)Tf0CBVh^|?1h>~fU3P_|mj$kZ>LGADmgB|NUH|{rg9fNTgGiC8_iY#$AbAX2
zSx4~qg@gKj$u6KO?k%Vt=yC7?vq!Ip2B==O@aTL28hJN>7zpY^8~~LkpaDh~P*WH*
zm~g1{DQM8f!telS@Bq;q1T`I4!4+XWsLKzkZoryB8Xrg80A=F4pmwqWDCanUo1wo2
zc-d`0%?iN~m4u_e1!Gh)AVX2EprS8DMF2GVaP&90_>0>GD*hr_?^!c2cv#+rmT;gZ
zI;gP@D)(+jJSaU3%6#<>pivf3E;Rs$M8eVE0x>EX;4;RQ0XotK^V8>u6P}hwOCN&*
zUIP@YKHWK>mg@&k#;?$t^yqJa8Wjof$jsqLhaX2TN4x-Od2#f0Bx{E?1H;kB5#S-M
z0&p)~;As8t2p5$IP{SV7^E!O=b0o*n-w{WS{*HMFjs=bv!aG6v_x}s-o&W!PbiRMV
z45C3p0WW^;0LeUm@o@(zvTIa4UOWeje1CC&2Pgmpa#RG4UXB3GJ9M{z8`AKmKY#0V
zP+PI_2vY>hWWE@&1n|Jc4zNLl5?T=`hj&90S}IrzG@;!Le*zk-5=>E1IC@*)3b;ZA
z4e}j7`g`((81Z^g^g;%7YOlSF{0nil0k|sgIQX0y>;@%}fzZ|kD7PTxJ;MX=+y_cP
zN5KgQJkoyvTo+p&u6_Q}=nuqL$Y38hK?s1F*&H4R>mgQvYRwmWwnLODfO09QVF9uW
z6t59b3qg6x@Bk?FW`Ig;P?-m6Y9D+K8O1&NIg-hP@z>Gc5o{iep#GsC3)B$@YahSV
z`}+SsD8(OrEO7YfV}a909}7G_`nX=;^~;BU|Nnm>u>JpkP;Uq}9~lf9WNdr`8hP~S
z4$=TQS0lxvn+Mtn>H!aJ8=mxNy~N+r0Lqh{*F3)81NF5Pe7X}XK-o~jM?1m5r#oN6
zr(3|MyV$|A*F}ZVqg&s%^(`p3a)1IkjKQT_*@xfdg^zZzh2ejn&JYy=P~Gms@A2KW
z^Sf{IH{aH8pqk0Sli%|nJVqvO0|ixxiiAgRjY@(?=RJ^t2B6s;aEB&ga-W8*Zs!!p
zfNAd(@E8<miUU-LgGO32z-_k97I04#l+zR@7ix%sMoB>9P~a-Dvjse23yKR+IS&>)
z4jJi#m8zh52^bqxWrNiCbV6)<!M^AJf3IE}=+ONu(1hZXI&pCN1I-}8CKZ3VSUxCy
z*Dc(7iNCM@0w@E3293RXO~4u*_kzqk`deT(sO1R?-H4;V1#3X9BmoPL?q-l{!2%i{
zy*$wAGv@$MWI_4^-QFIa$@hJFRSxz#F?PDB_<;Os;L$5{1MCn`(dqyiJ1+$lSS9r;
z3jBSI;Ay-P6$9VS=b-cgQP|CT(+phF#i&R?s@&ropmu54kE5>z4uK1A=pcmQR&YKC
zH~b=5KZA6=j=XgAw}9z%2wU_zhz1pD5=UQ0TsU&jJJS8f(bq8{D)kyiUq@az`Z_X1
zh2!Ysh+96EpG(`0K8~nSQ2;sSX803OooxX2umz-{=F!b+2sRFs!BHIsa$Y>T^T4IF
z#Eaf7U`NKNNF4nQp49+3_i=qBt1rl~$C0<dgC`0{Uq^~QHU;?`G_w%_YL^}Ym)smj
zZ$~`vw7gf^aP)RWii*Y2--4+OKFld90glZ_89kVDR3cogb5wlzTjziWRa(H=1zyyZ
zI)F-71CPcdp!@|N;i!*)S+*HA#&Yyz__w2{BT`fpUM~VAS@61+U7+Ofvf|JG|1WGn
zDiPJ!1W@LM3}t}R4APPU&>U2Eg@#Xek%mX7i3+%$_vkLLK*`RaMF#bTCq2639ei8=
zdvtU7^1J-?&@M3WZT;rU@AB0{JKn<Zzel$aC~--6^tv&EQlkK<!3D~D8Xo)}_d%K7
zpfg5A#UuGPcs#`OHz>z?^1J+fv2PP7$C{{sM!zaNJCA@gIUqHEj=QLMM1XVUaTk>U
zP_8=qTYwLwDgsml3V?zj0VE7!fF?f~Kn(CWD!4$k0MCnT0Z$IW@+??6sFgiA+(;Qb
z*ADVOf2S^dP?H<vD$sB|sHkONVDRib3TnK7+I*lk1bCVs+{}FkDrg*DNP{c0_b)bY
z1o;Lu)Ad^b6!|saA-BDtObSk&5+1F$`TJZ!**V!oMdj$h@c;EbosW<H7PL|E?A18}
zpLFSF%`^rzalkVa!HCq#0U6Odvk{z5LHPlcPG3hJf^`2QS*wg07><69<OS2eBMu(@
z9peIWL?kGI!Wsw2qkQ$HlR@h@es2J|9MXk=r$_$QRSck1F2JG!8g^q~V1TE|_`}cw
zphiUl9+}s`hWMy}*C&9=AW-mm_1Zw@HM>Jp43Gj6RAZT__=WYdKvn^PdU>GS3!QNP
z9LeOvc<AWw2u2UapP-u4<tS(tiP6LISnZ>iwg3MA2Ri^1gy7(NY5xEJ{}-wopyHDo
z;<BYr8Zs~({Ta@xX~e*A7#6lie+#g37)@r2_h2;K>>qE*$mFtaa#w;K6U(~ErxH?_
z?AA^;NIcH8b;IPJiTO+}8zyHZRWa&KewviQ6u)}1f3h@_)9T4N$w5qft0x~wwo<#a
z>i>VkOQ7{EouE<t)&r%|9=)baS`6SNNw(HH3=A*yR!vq&33Rt!1*vOaKKl3nzpLSY
zk6zP0O;F$I7kEL^5p4#B7u#3<|Nnw%6=+qG=xS|<A*NlLppKyJX^;Z7RhxID%w=Ta
z*Pa}nroyGD4f3Nt?@O)8t!dMkJhdjvrmJv8Yk@^sm9;i|r=MeFRG%!B8O5ZrbaGjy
zoJs~O1A|YpIw+%RfbwPrIPa!-bVJ5kK`A|A@{UXu)elQRt+w}|+{W?ZF@yseE4%??
z9(!?i$z;YXt$N7(r3b&~3y*_OL9L_Lp!Q$)5^y^ZWFV-nd{GV22jYWj+!xs}{?Qi+
zAf?A46)z~&fb8_}=w{7OV_*QO2i2**txxJ4Ji1v_48MU|vko4;GLX8h)BhHzySR5T
zI9IUlRRe|T3(wvVSq4ZV0ZkG>*L{2jO>u(y+MvNV2T%(IG`+<PU5?WYnw0H)HThhj
zq#bDb{&<Uu1~YWz5!9yzsqcn#A3^i~VEx@CDxe`r6$7}!IVu&*lTYQ!hJjX8fV6;S
zSwVE?aZnlL0kaJ>5X6ia^GrZuXTS^tt>}TVK}(!oc(0s1D_gF<dl4vZ9vOZE1$5ma
zQ0|pc0gcHpFoS&m%5x7W7{OCPpy?r>&Lt`mASZM~vnL145isWqz-;fFqhbM830iQ_
z3Gxdh3Lu_!QON-LVX|ws=;R4G-1Q(1HV=U;i~tz}8fb+X4VwFcc?QX2|1WrS-tg#+
z{tud4>U$4byVlL30&2^89Q{?tYpDWC{UYGF2G78Nd<Lq<S=f(GzMdn`0&?!;2^A)j
z^>St8qZa=EfApt@j|#Zw4e};Pxku+AP%>g*o}8JhT@QCKD2{s}(?%eV1RyySYixA3
zsDOh3+}{TWE>aYOSDG}x;Q*yw3$TTAR2rZ`PyzFHXOBulCIbV*{}L4iAJD7{s2Z~X
z)%yxByuO2?agPcpDvtgZgyaU$$_15#7j_FkMFQ(~B~a{w##%t`O9od$Eh^xE29LFX
z*3yDJ3i8Q|SMw*=<VlKw^WLxdp!D+og%QY{ZdPff$*c2{>LH2cB~pdk4IWN!fg~2F
z=RocR&me1n0tpnPpc>hyR|Pa)2p;zY#Q>-Sq5v8M1^LyxS4JJ=DFNhc0*(=o=fE}a
z@fMZI_X-tEL1w}ugxD-G`C-0HJ*bKS4XA^XIw;T!KpfDbK@cawr}Kp;XjQ;pa0mwY
zc7FHlJa!z?d<PBDgNp`_W*ZfTIzFGyM;@K;d^-QXV4eqx0~QsJ&NV8K(D{TktJ>`Y
znrzqt2_5s!Es&50FU^7FBY5z1hp1$Lq5>Q?ojs6*1Im7r&lMWPgJeJfMV9ZOiOKQ*
zC9wPX`#@`R4G)0E2f)5M`g^i;k(|hV&>9O+Ao+Cu01b+Os=&!UMV9rTRskz$rIjFb
zZHABKgVLv<`10`SeB{&l?nUbyXfQOt;qdMJ=&=)&iNFaS<WJBTG02}F$AjV%lpH{8
zPzeK$TUb8n>`_sG<^WJR3YI(GGC87Hv>p@y;M~!<MFkZ7ApM|J3o{FprorkvdsIMa
z9V8Db6~SUig$GC+6p)}MFjA%pQE`CgtRC=iHpp`dFZO)`<sw^mIfi-$kDZ|Q94Ny;
zi>P97ZrTHmQ=iW7s5vkLw3G|vz8BfEK|!+xJcbA!YXr?&x6KA+ukSDRfhBF1%R(lj
zKx=70AqQFs0O|@xL_s6$8EAJ0D93Ps+Ce@l4&YJl2vF<t$iZUJl9{8wBOqh?;3bQo
z2zmzc5C<r!g66hgT$u&(3u~<`s6_Y<%{Z{pS8#};ln9`5>%|6;p6(Krf?X1zHPs-O
z=YfPeUqI4M0|Ue4btTdY2{4Zx{Vmu7UY+-1!YqgnPk_`<ep#Zx=rftCbT{jEP-}Mb
zsZx1iXBlv#<0gOmd~j#cMMYxrpHgirR**QT9U5@-X9%k}i1$<fMH#4%Wjj#{+`tS`
zi8%UG{}+F|J;=0f78OtmV_=?~SSIWUN>q>%2%g_UK$E2<h!P&0X5jU6jfwy$ae)#A
zs9(!5`COr@FDylZy3wE{1Fx+?t8tm3)ju>4VMz>>>|ktA(gP(iklupH>g9U%`+tI)
z8<4dtzMcO(JKy_sz6H1Q4PLUcGB6-|VW5e*6(FC!48L^rr~Xa;wnR|>uUFPt3KV%C
zd^`U^t!@POHb4e}o62(_Spy!CpkY<0Ph4N5Pyhda7ig0KNVM0KRSL8^cLpeNf+83a
zzwciNP6x%?`xh&xgJeL%BE77SBtcdB9F-RglTVe)*MruTwSp2qXh@oofdMp-4IQZl
zCEgcPrh&AAx`?kp-UsD;aEQZNLLjxUnzMTjcvJ<X4wUo3!3QoMQJPX)R3Nbfu59XC
zz$GnGB?N1Nbi(SPEzr^vy{UyO@}LzL$Q3vv5n718Xs-vAb)b&f5*E-R$|)*q7{KlB
z?jCR+2dz&9HCI6sJ4YeQcEPg<uw|(iN?(H3vx1Zy{TY71qw|ypBW$hF(cckF9?gfq
zD~>;dhKWHl70{^e1ci9FGlxfgcQAOc_kzID&oMq1I6Qh;MI}HH)PffKpi~XY#^B`|
zAHq)_{jBB9a`bo1Pe`IW`dbi^7+;u81r-fjz^kgDIr;sIK#*X!i3)N!uTeR`#9#xS
z@&`||pJ0OIgrmO&CWCsdARBLhrn;LUt56w0A$kEis}8E#e@wnutvvaE78~Q1$^6+;
zlg}6OOrBezUC-b52BfVUQcQvh7ms8}z<P9k@aPq}1}f<w>v_P#>7cIZp`+lnfrm<;
zgG1`*Z$TavP*?i^6KExU3TRZ%^rRSQjRkZ?qJl@~ec#Rp|1W?h*i>L11(l!jpcO+U
zDhVJH5k=VB$&jp|I%V>ON{M<<jDbgE`8^JSCc!}MB%e+j6`yW!36i5C(Rq|#p24G+
zRZt9+O+l?LPz}%xia5|@E~sYP0jfkm83q(opvZaQG8wdPEF2Wmp55-?8LsBTjEGsT
z5U?wr?E=*R44{4sNP|!3d(d>%i;Btr|6AUz?eXchQ2`D5OMvRF1dm?Ueo^q$6DX_p
zLhT1tE|AEC*$-;_yl|Zi%B`R=_S=)d{jV37L3B5CjngqGX98&C<Hb%euSUh;#rjF0
z!1?}SIhg+cVm?%c<Hb}kt3@S(6|~_5x*~_aRf-W@bi1g42CP7RN~FO7=qOw->rD|*
z$YHIW;=w@&s(;}@#{;$rQF?uUacUySGbJkU3JH_|LEZ!xOd$23I?QF#|No#g_HqKK
z1r04dCMV~Z#DSD~wBD{Xg)F7p12PWU;R8#7);KaVFuZL305!q+FR1$ve&FT0-~a!E
zdLK8185myjfffq+PeflRG}*pJmQw^A47{8olVfX?OfG<1o1jI6ppI5=j*7w2&jJU1
zdTUf9e0ob%6kdV~T1es-1g&iDWepadyt_uJ-W@bs+cn|;e`q@X2kO|ZLh!*Iy;&0=
zF?wkNxb9(50oACf!r*$B)l`^)p)43Q-vA!J3{g=y`a9wfxRK%0`O>5F-qGKYtP_Ml
z5%L`DN{ts>uOTB6M}OA~fL7G~)W5;s<_2C2F6g5o@zU}C|NoFIFi{9pBE3HjnqmW|
zKhT;&q!@HjF?eyd|NsA8kO+P6(QDcu#K7>P65?X0r{2A|*$?r7AV?fE*TV!!RFiMl
zN?En{gRBGPD2C7ca=#qIeR^%Rbr={N!#q2$Ifi<6esv6S?EDkz(W@(?!@v;iqxsjP
z`3=X5`F)%9>n<{G{=hB3f*>Y;;9jt)!GlR6H7~U&Ju|gfFU`_`Aw4fYFDE})FQtTG
zvPaW8r3d!+T=$){*(>Y-p^Y6t34?)w0qhJAds4G(ba6?3K|yK?Ls3S)Zf<5?W-ds7
zZmB|EemO`9GG6gz`t0(@KkfOeer#Vlb(KZ6$U66jRW0^bJqC~ZDqZa$YCWOWcpjKs
z(JaF#GI?sVc|B;qvJ3+QgGSB&{|bx@3>CHi|Ic7&V92QZ{~xsQHKYFjf6%hT6Al0W
zgVqgeH2?p9g^_{5qUHa8(3+Hrw*UV@tC}>r{{P>=#K2I|_5c4FP&wT5|39c1HDmJs
z|DcA9#Pt9FcQ7+BfZB2#EDQ`c7XAMp!NO3_P_g3w{}n6@3>iEA|Np_lz_0?eT#uE3
zq2uKL|1($_7-pRL|NjLm1B1rd|NsB6GBE5o`~SZH8w10QOaK2TurV;~xcvWr0~-Uw
zjO+jZUtnWksCfPVzY9A9gT>qb{|~S;Fsyj{|33!@Bm`hvl^Lso7#J%A7^Qh8?`q+b
zim3hnUmv8xg&nl-k^!_n)C9z@p8TjqjY*<zGFPh~lScjIs8%Z`pT^14TAi3AnkL_B
zbz>E1`Tu|QWZia^$)0WMas}=G|NDV7I`RoLGZk?0aX4~&FfcGYU}Rue(muJb&5KE+
zbMmb=KPHU{leOCYm?9=mu4?yV+A(eNsdhi6ig}ZzI{es9EdKxhA^&8(b{QTQK8t2{
zcfJWM%r*Rz{o17Xd>9xQctD}G`Tu_h-pP9H%B&MOnWuA4W^0$@jA3D5h&b^7|CPzQ
z?MmEUd>fdUTbQ}{EF33CwM*5nVPRl+arpm#ZqV*XCq99GCRaX*K4vFAg<cj9K8+q$
zM?Ql#Hkge}%-&2Ot$YTKd>T%C3Ql|yPJ9ARd>n3|y{@2O`U8rrWB>ny8WkWvgTw?_
z85mTKBg8<?RbgddSaAG5^M5JU$@ki1CI@wDX_c@tFch5m|NlBjBRH^_wsY}ucyI^t
zae(wpVP#-2IQ#$q4d%&t9g35WcdBwdVP#-AaN+;|i^3pNBAh`X<HNUro4JXFZv}S-
z-v(}`t$ZuE9YI3e%zP`jnf8Ii=EB4taA5V0J9hvB149HG1B1?!|NnO~PR{F4p4{7|
z&b9`W#!4se>QZ96!p6W*@O1K{E+<YIb_RwSPyhdC-TaU{nsIVQ?+&9CpawJp1A`>U
zFQBsppqurjp?qlCV~~RK3)W5!>Z@nmu=!G-J|pAW$&O57lO-l>WZVE|J(wWLIAQaz
z34M(9FoU7`z!WnBKjN?hun2@;W)NV29Yp{Ufso7$f($RFK-kD6GlLKV>{J3&abX7N
zX#=3j8X?8ZAi@AU0s$e)z`)EP$^bhh099O!fdRC46j_h~q+Xoi09J7c1_kIb1Q0Eg
z+a@<jU4S0504kzD!p7k903;Y3HcnQWQZ5PA0O}=yWI>CW7#To5n7nC<KI7WSkEZl7
ziOiW?G*wR!J#;Qi-Z)jk6B<__%^<7_Hc1E;1Rx<WHiL?z2bm{$BO}PClclD`2~Nlc
zt7f<e-4}m!vLmzj<OS0>4WSN&NHZ|VvobI!FkptnD@FzeL2wjLzBNrxssI`X;HAS1
z3=BF<3=BdHybKD7la;3HNnJ>Tm<P^93=9l4Od#Fj3=cL;&YE5hJAnYAhJitt8KSu%
z79s*J2N)O_bfMx2pp3`Bz~BHqWgrMDZXf|sZx0&5Vqjp%oUAxQ+;9UM1A`C~9|L+K
z?Sbkluz~1$2-?}sz`(F}^2QnBj0Yw+&JdTp4plQD0ip(0g8Txhk(6YCxKltIBD)Gy
zUNJB*C{KPkLtN4gD(j#Ck#z#~Mj03w;wKBv6ql@KVPFts5@(P=3-5NY>`bU$1?UkJ
zOyF2$VAuwdWqbf8zf4{@Q`}RQm4QJ3OZ*rzg2W{lJT^gWhml&W5I4gP)Chu_>^+%r
zmN;V(h?JZLHd6w_ZHqv%hR2{L7(h#tRH)I{pyIF-IY1NhAgLc9U6SH#pcv<8fF1Aw
zs<A<`Hj@v|5@(DCli8C$&Jt(rn%p>B+;ABv));sh7C>zQ9gF}{zYQwBAsxhHV7LJd
z^OIl`Zh^_4AQL8Un9ahdFxinwm{Do+#@XUpf$R{UZ%Bmbl!6*r3>9Ax2@$`-y!qB_
z9cH$YwD^MjlF5PdB{y@=*JGI+G+&q-A{L)sl$f&FVzE2p=7uE-%uJ>QlNT-5<ANoy
z6^)Z`EjJec^_V~=Lg^oq6<6rP^ENDxznPr4LSG10?!jmS_Q{J@=yO7|e8$|#w^o?L
zcm;DOORa?SD&|g(T4~M++PuiXz;I*!<XJ1tIbns?g2}fKyaNj+ORX~JgeEtMg_EOJ
y*|RKUU|^VB$Sgd0!zxV<5C_x*nEY^+Ii~=qpU=R+Fln;lYICMVn<q!DRs#SI2erQd

diff --git a/minim.f90 b/minim.f90
index a3edb23..538d4aa 100644
--- a/minim.f90
+++ b/minim.f90
@@ -33,12 +33,13 @@ L = 8.0_dp
 !g_ptr => get_vw_gradient
 !f_ptr =>get_nagy_total_energy
 !g_ptr => get_nagy_gradient
-!f_ptr => get_vw_total_energy
-!g_ptr => get_vw_gradient
-f_ptr => get_tf_pot
-g_ptr => get_tf_pot_grad
-bas = basis_t(nel,20,2000,L)
-
+f_ptr => get_vw_total_energy
+g_ptr => get_vw_gradient
+!f_ptr => get_tf_pot
+!g_ptr => get_tf_pot_grad
+bas = basis_t(nel,60,2000,L)
+!f_ptr => get_thomas_fermi_total_energy
+!g_ptr => get_thomas_fermi_gradient
 !call deq_optimize(bas,f_ptr,g_ptr)
 call optimizer(bas,f_ptr,g_ptr)
 !L = get_nagy_total_energy(bas,0.0_dp,0.0_dp)
diff --git a/optimize.f90 b/optimize.f90
index d33b8ae..72db14e 100644
--- a/optimize.f90
+++ b/optimize.f90
@@ -56,10 +56,11 @@ subroutine optimizer(bas,defunc,gradfunc)
     procedure(func)   :: defunc
     procedure(gs)     :: gradfunc
     real(dp)          :: mu, a,b,c, optres
-    a = -50.0_dp
+    a = 0.0_dp
     b = 50.0_dp
     c = 25.0_dp
     optres = 10.0_dp
+    !goto 10
     do while(abs(optres) > 1e-2_dp)
         c = (a+b)/2.0_dp
         optres = deq_optimize(bas,defunc,gradfunc,c)
@@ -71,23 +72,25 @@ subroutine optimizer(bas,defunc,gradfunc)
         print *, a,b
     end do
     print *,"mu now", c
-    
+    !10 continue
+    !optres = deq_optimize(bas,defunc,gradfunc,a)
     
 end subroutine
 function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron)
     type(Basis_t)     :: bas
     real(dp)          :: grad(bas%nfuns),gradprev(bas%nfuns),b,stp,etp
     real(dp)          :: sn(bas%nfuns),snprev(bas%nfuns), curr(bas%nfuns),prev(bas%nfuns),gamma,mu, psi(size(bas%x))
-    real(dp)          :: rho(size(bas%x)), prevenergy,currenergy,nelectron
+    real(dp)          :: rho(size(bas%x)), prevenergy,currenergy,nelectron,pi
     integer           :: i,j
     procedure(func)   :: defunc
     procedure(gs)     :: gradfunc
     !mu = 8.7_dp
+    pi = 4.0_dp*atan(1.0_dp)
     stp = 0.0_dp
     etp = 10.0_dp
-    
     prev(:) = 0.0_dp
-    prev(:) = 0.01_dp
+    prev(5) = 0.1_dp
+    !pi = defunc(bas,0.0_dp,prev,prev,mu)
     psi(:) = expand_psi_ofdft(bas,prev)
     !print *, (bas%x(2)-bas%x(1))*sum(psi(:)**2.0_dp)
     curr(:) = 0.0_dp
@@ -104,6 +107,9 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron)
   
     prevenergy = defunc(bas,stp,sn(:),curr(:),mu)
     do i=1, 100000
+        psi(:) = expand_psi_ofdft(bas,curr)
+        !print *, sum(psi(:)),"lol"
+        !print *, "mi"
         gradprev(:) = grad(:)
         grad(:)     = gradfunc(bas,curr,mu)
         b           = beta(grad,gradprev)
@@ -126,7 +132,7 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron)
         print *, i, sum(grad*grad)
         if(sum(grad(:)*grad(:)) < 1e-8_dp .or. &
           (abs(defunc(bas,stp,sn(:),curr(:),mu)-prevenergy )<1e-8_dp .and.sum(grad(:)*grad(:)) < 2.0_dp ) ) then
-            print *, "stopped",i
+            !print *, "stopped",i
             
             
             
@@ -134,10 +140,11 @@ function deq_optimize(bas,defunc,gradfunc,mu) result(nelectron)
             psi = expand_psi_ofdft(bas,curr(:))
             rho = psi**2.0_dp
             nelectron= sum(rho)*(bas%x(2)-bas%x(1))-bas%nelectrons
-            print *, nelectron, defunc(bas,0.0_dp,sn,curr,mu)+mu*sum(bas%dx*psi**2)
+            print *, "electron difference",nelectron
             do j=1, size(bas%x)
                 write(12, *) bas%x(j), rho(j)
             end do
+            !print *, curr
             close(12)
             exit
         end if
-- 
GitLab