From 7d2310bc862543bf4461a8a5e9149c8950e4c278 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sun, 3 Jan 2021 18:45:58 +0100 Subject: [PATCH] Initial commit --- QMC.org | 788 ++++++++++++++++++++++++++++++++++++++++++ energy_hydrogen.f90 | 39 +++ energy_hydrogen.py | 23 ++ hydrogen.f90 | 25 ++ hydrogen.py | 13 + plot.png | Bin 0 -> 7938 bytes plot_hydrogen.f90 | 28 ++ plot_hydrogen.py | 21 ++ plot_py.png | Bin 0 -> 40693 bytes variance_hydrogen.f90 | 57 +++ variance_hydrogen.py | 36 ++ 11 files changed, 1030 insertions(+) create mode 100644 QMC.org create mode 100644 energy_hydrogen.f90 create mode 100644 energy_hydrogen.py create mode 100644 hydrogen.f90 create mode 100644 hydrogen.py create mode 100644 plot.png create mode 100644 plot_hydrogen.f90 create mode 100644 plot_hydrogen.py create mode 100644 plot_py.png create mode 100644 variance_hydrogen.f90 create mode 100644 variance_hydrogen.py diff --git a/QMC.org b/QMC.org new file mode 100644 index 0000000..c2b8394 --- /dev/null +++ b/QMC.org @@ -0,0 +1,788 @@ +#+TITLE: Quantum Monte Carlo +#+AUTHOR: Anthony Scemama, Claudia Filippi +#+SETUPFILE: https://fniessen.github.io/org-html-themes/org/theme-readtheorg.setup + +#+STARTUP: latexpreview +#+STARTUP: indent + + +* Introduction + +We propose different exercises to understand quantum Monte Carlo (QMC) +methods. In the first section, we propose to compute the energy of a +hydrogen atom using numerical integration. The goal of this section is +to introduce the /local energy/. +Then we introduce the variational Monte Carlo (VMC) method which +computes a statistical estimate of the expectation value of the energy +associated with a given wave function. +Finally, we introduce the diffusion Monte Carlo (DMC) method which +gives the exact energy of the H$_2$ molecule. + +Code examples will be given in Python and Fortran. Whatever language +can be chosen. +** Python + +** Fortran + +- 1.d0 +- external +- r(:) = 0.d0 +- a = (/ 0.1, 0.2 /) +- size(x) + + +* Numerical evaluation of the energy + + In this section we consider the Hydrogen atom with the following + wave function: + + $$ + \Psi(\mathbf{r}) = \exp(-a |\mathbf{r}|) + $$ + + We will first verify that $\Psi$ is an eigenfunction of the Hamiltonian + + $$ + \hat{H} = \hat{T} + \hat{V} = - \frac{1}{2} \Delta - \frac{1}{|\mathbf{r}|} + $$ + + when $a=1$, by checking that $\hat{H}\Psi(\mathbf{r}) = E\Psi(\mathbf{r})$ for + all $\mathbf{r}$: we will check that the local energy, defined as + + $$ + E_L(\mathbf{r}) = \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}, + $$ + + is constant. + + +** Local energy +:PROPERTIES: +:header-args:python: :tangle hydrogen.py +:header-args:f90: :tangle hydrogen.f90 +:END: +*** Write a function which computes the potential at $\mathbf{r}$ + The function accepts q 3-dimensional vector =r= as input arguments + and returns the potential. + + $\mathbf{r}=\sqrt{x^2 + y^2 + z^2})$, so + $$ + V(x,y,z) = -\frac{1}{\sqrt{x^2 + y^2 + z^2})$ + $$ + + #+BEGIN_SRC python +import numpy as np + +def potential(r): + return -1. / np.sqrt(np.dot(r,r)) + #+END_SRC + + #+BEGIN_SRC f90 +double precision function potential(r) + implicit none + double precision, intent(in) :: r(3) + potential = -1.d0 / dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) +end function potential + #+END_SRC + +*** Write a function which computes the wave function at $\mathbf{r}$ + The function accepts a scalar =a= and a 3-dimensional vector =r= as + input arguments, and returns a scalar. + + #+BEGIN_SRC python +def psi(a, r): + return np.exp(-a*np.sqrt(np.dot(r,r))) + #+END_SRC + + #+BEGIN_SRC f90 +double precision function psi(a, r) + implicit none + double precision, intent(in) :: a, r(3) + psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )) +end function psi + #+END_SRC + +*** Write a function which computes the local kinetic energy at $\mathbf{r}$ + The function accepts =a= and =r= as input arguments and returns the + local kinetic energy. + + The local kinetic energy is defined as $$-\frac{1}{2}\frac{\Delta \Psi}{\Psi}$$. + + $$ + \Psi(x,y,z) = \exp(-a\,\sqrt{x^2 + y^2 + z^2}). + $$ + + We differentiate $\Psi$ with respect to $x$: + + $$ + \frac{\partial \Psi}{\partial x} + = \frac{\partial \Psi}{\partial r} \frac{\partial r}{\partial x} + = - \frac{a\,x}{|\mathbf{r}|} \Psi(x,y,z) + $$ + + and we differentiate a second time: + + $$ + \frac{\partial^2 \Psi}{\partial x^2} = + \left( \frac{a^2\,x^2}{|\mathbf{r}|^2} - \frac{a(y^2+z^2)}{|\mathbf{r}|^{3}} \right) \Psi(x,y,z). + $$ + + The Laplacian operator $\Delta = \frac{\partial^2}{\partial x^2} + + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$ + applied to the wave function gives: + + $$ + \Delta \Psi (x,y,z) = \left(a^2 - \frac{2a}{\mathbf{|r|}} \right) \Psi(x,y,z) + $$ + + So the local kinetic energy is + $$ + -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (x,y,z) = -\frac{1}{2}\left(a^2 - \frac{2a}{\mathbf{|r|}} \right) + $$ + + #+BEGIN_SRC python +def kinetic(a,r): + return -0.5 * (a**2 - (2.*a)/np.sqrt(np.dot(r,r))) + #+END_SRC + + #+BEGIN_SRC f90 +double precision function kinetic(a,r) + implicit none + double precision, intent(in) :: a, r(3) + kinetic = -0.5d0 * (a*a - (2.d0*a) / & + dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) ) +end function kinetic + #+END_SRC + +*** Write a function which computes the local energy at $\mathbf{r}$ + The function accepts =x,y,z= as input arguments and returns the + local energy. + + $$ + E_L(x,y,z) = -\frac{1}{2} \frac{\Delta \Psi}{\Psi} (x,y,z) + V(x,y,z) + $$ + + #+BEGIN_SRC python +def e_loc(a,r): + return kinetic(a,r) + potential(r) + #+END_SRC + + #+BEGIN_SRC f90 +double precision function e_loc(a,r) + implicit none + double precision, intent(in) :: a, r(3) + double precision, external :: kinetic, potential + e_loc = kinetic(a,r) + potential(r) +end function e_loc + #+END_SRC + +** Plot the local energy along the x axis +:PROPERTIES: +:header-args:python: :tangle plot_hydrogen.py +:header-args:f90: :tangle plot_hydrogen.f90 +:END: +:LOGBOOK: +CLOCK: [2021-01-03 Sun 17:48] +:END: + +For multiple values of $a$ (0.1, 0.2, 0.5, 1., 1.5, 2.), plot the +local energy along the $x$ axis. + +#+begin_src python +import numpy as np +import matplotlib.pyplot as plt + +from hydrogen import e_loc + +x=np.linspace(-5,5) + +def make_array(a): + y=np.array([ e_loc(a, np.array([t,0.,0.]) ) for t in x]) + return y + +plt.figure(figsize=(10,5)) +for a in [0.1, 0.2, 0.5, 1., 1.5, 2.]: + y = make_array(a) + plt.plot(x,y,label=f"a={a}") + +plt.tight_layout() + +plt.legend() + +plt.savefig("plot_py.png") +#+end_src + +[[./plot_py.png]] + + +#+begin_src f90 +program plot + implicit none + double precision, external :: e_loc + + double precision :: x(50), energy, dx, r(3), a(6) + integer :: i, j + + a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) + + dx = 10.d0/(size(x)-1) + do i=1,size(x) + x(i) = -5.d0 + (i-1)*dx + end do + + r(:) = 0.d0 + + do j=1,size(a) + print *, '# a=', a(j) + do i=1,size(x) + r(1) = x(i) + energy = e_loc( a(j), r ) + print *, x(i), energy + end do + print *, '' + print *, '' + end do + +end program plot +#+end_src + +To compile and run: + +#+begin_src sh :exports both +gfortran hydrogen.f90 plot_hydrogen.f90 -o plot_hydrogen +./plot_hydrogen > data +#+end_src + +#+RESULTS: + +To plot the data using gnuplot" + +#+begin_src gnuplot :file plot.png :exports both +set grid +set xrange [-5:5] +set yrange [-2:1] +plot './data' index 0 using 1:2 with lines title 'a=0.1', \ + './data' index 1 using 1:2 with lines title 'a=0.2', \ + './data' index 2 using 1:2 with lines title 'a=0.5', \ + './data' index 3 using 1:2 with lines title 'a=1.0', \ + './data' index 4 using 1:2 with lines title 'a=1.5', \ + './data' index 5 using 1:2 with lines title 'a=2.0' +#+end_src + +#+RESULTS: +[[file:plot.png]] + +** Compute numerically the average energy +:PROPERTIES: +:header-args:python: :tangle energy_hydrogen.py +:header-args:f90: :tangle energy_hydrogen.f90 +:END: + + We want to compute + + \begin{eqnarray} + E & = & \frac{\langle \Psi| \hat{H} | \Psi\rangle}{\langle \Psi |\Psi \rangle} \\ + & = & \frac{\int \Psi(\mathbf{r})\, \hat{H} \Psi(\mathbf{r})\, d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} \\ + & = & \frac{\int \left[\Psi(\mathbf{r})\right]^2\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} + \end{eqnarray} + + If the space is discretized in small volume elements + $\delta x\, \delta y\, \delta z$, this last equation corresponds + to a weighted average of the local energy, where the weights are + the values of the square of the wave function at $(x,y,z)$ + multiplied by the volume element: + + $$ + E \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\; + w_i = \left[\Psi(\mathbf{r}_i)\right]^2 \delta x\, \delta y\, \delta z + $$ + + We now compute an numerical estimate of the energy in a grid of + $50\times50\times50$ points in the range $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$. + + Note: the energy is biased because: + - The energy is evaluated only inside the box + - The volume elements are not infinitely small + + #+BEGIN_SRC python :results output :exports both +import numpy as np +from hydrogen import e_loc, psi + +interval = np.linspace(-5,5,num=50) +delta = (interval[1]-interval[0])**3 + +r = np.array([0.,0.,0.]) + +for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: + E = 0. + norm = 0. + for x in interval: + r[0] = x + for y in interval: + r[1] = y + for z in interval: + r[2] = z + w = psi(a,r) + w = w * w * delta + E += w * e_loc(a,r) + norm += w + E = E / norm + print(f"a = {a} \t E = {E}") + + #+end_src + + #+RESULTS: + : a = 0.1 E = -0.24518438948809218 + : a = 0.2 E = -0.26966057967803525 + : a = 0.5 E = -0.3856357612517407 + : a = 0.9 E = -0.49435709786716214 + : a = 1.0 E = -0.5 + : a = 1.5 E = -0.39242967082602226 + : a = 2.0 E = -0.08086980667844901 + +#+begin_src f90 +program energy_hydrogen + implicit none + double precision, external :: e_loc, psi + double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm + integer :: i, k, l, j + + a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) + + dx = 10.d0/(size(x)-1) + do i=1,size(x) + x(i) = -5.d0 + (i-1)*dx + end do + + delta = dx**3 + + r(:) = 0.d0 + + do j=1,size(a) + energy = 0.d0 + norm = 0.d0 + do i=1,size(x) + r(1) = x(i) + do k=1,size(x) + r(2) = x(k) + do l=1,size(x) + r(3) = x(l) + w = psi(a(j),r) + w = w * w * delta + + energy = energy + w * e_loc(a(j), r) + norm = norm + w + end do + end do + end do + energy = energy / norm + print *, 'a = ', a(j), ' E = ', energy + end do + +end program energy_hydrogen +#+end_src + +To compile and run: + +#+begin_src sh :results output :exports both +gfortran hydrogen.f90 energy_hydrogen.f90 -o energy_hydrogen +./energy_hydrogen +#+end_src + +#+RESULTS: +: a = 0.10000000000000001 E = -0.24518438948809140 +: a = 0.20000000000000001 E = -0.26966057967803236 +: a = 0.50000000000000000 E = -0.38563576125173815 +: a = 1.0000000000000000 E = -0.50000000000000000 +: a = 1.5000000000000000 E = -0.39242967082602065 +: a = 2.0000000000000000 E = -8.0869806678448772E-002 + +** Compute the variance of the local energy +:PROPERTIES: +:header-args:python: :tangle variance_hydrogen.py +:header-args:f90: :tangle variance_hydrogen.f90 +:END: + + The variance of the local energy measures the intensity of the + fluctuations of the local energy around the average. If the local + energy is constant (i.e. $\Psi$ is an eigenfunction of $\hat{H}$) + the variance is zero. + + $$ + \sigma^2(E_L) = \frac{\int \left[\Psi(\mathbf{r})\right]^2\, \left[ + E_L(\mathbf{r}) - E \right]^2 \, d\mathbf{r}}{\int \left[\Psi(\mathbf{r}) \right]^2 d\mathbf{r}} + $$ + + Compute an numerical estimate of the variance of the local energy + in a grid of $50\times50\times50$ points in the range $(-5,-5,-5) \le \mathbf{r} \le (5,5,5)$. + + #+BEGIN_SRC python :results output :exports both +import numpy as np +from hydrogen import e_loc, psi + +interval = np.linspace(-5,5,num=50) +delta = (interval[1]-interval[0])**3 + +r = np.array([0.,0.,0.]) + +for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: + E = 0. + norm = 0. + for x in interval: + r[0] = x + for y in interval: + r[1] = y + for z in interval: + r[2] = z + w = psi(a, r) + w = w * w * delta + El = e_loc(a, r) + E += w * El + norm += w + E = E / norm + s2 = 0. + for x in interval: + r[0] = x + for y in interval: + r[1] = y + for z in interval: + r[2] = z + w = psi(a, r) + w = w * w * delta + El = e_loc(a, r) + s2 += w * (El - E)**2 + s2 = s2 / norm + print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}") + + #+end_src + + #+RESULTS: + : a = 0.1 E = -0.24518439 \sigma^2 = 0.02696522 + : a = 0.2 E = -0.26966058 \sigma^2 = 0.03719707 + : a = 0.5 E = -0.38563576 \sigma^2 = 0.05318597 + : a = 0.9 E = -0.49435710 \sigma^2 = 0.00577812 + : a = 1.0 E = -0.50000000 \sigma^2 = 0.00000000 + : a = 1.5 E = -0.39242967 \sigma^2 = 0.31449671 + : a = 2.0 E = -0.08086981 \sigma^2 = 1.80688143 + +#+begin_src f90 +program variance_hydrogen + implicit none + double precision, external :: e_loc, psi + double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm, s2 + integer :: i, k, l, j + + a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) + + dx = 10.d0/(size(x)-1) + do i=1,size(x) + x(i) = -5.d0 + (i-1)*dx + end do + + delta = dx**3 + + r(:) = 0.d0 + + do j=1,size(a) + energy = 0.d0 + norm = 0.d0 + do i=1,size(x) + r(1) = x(i) + do k=1,size(x) + r(2) = x(k) + do l=1,size(x) + r(3) = x(l) + w = psi(a(j),r) + w = w * w * delta + + energy = energy + w * e_loc(a(j), r) + norm = norm + w + end do + end do + end do + energy = energy / norm + + s2 = 0.d0 + norm = 0.d0 + do i=1,size(x) + r(1) = x(i) + do k=1,size(x) + r(2) = x(k) + do l=1,size(x) + r(3) = x(l) + w = psi(a(j),r) + w = w * w * delta + + s2 = s2 + w * ( e_loc(a(j), r) - energy )**2 + norm = norm + w + end do + end do + end do + s2 = s2 / norm + print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2 + end do + +end program variance_hydrogen +#+end_src + +To compile and run: + +#+begin_src sh :results output :exports both +gfortran hydrogen.f90 variance_hydrogen.f90 -o variance_hydrogen +./variance_hydrogen +#+end_src + +#+RESULTS: +: a = 0.10000000000000001 E = -0.24518438948809140 s2 = 2.6965218719733813E-002 +: a = 0.20000000000000001 E = -0.26966057967803236 s2 = 3.7197072370217653E-002 +: a = 0.50000000000000000 E = -0.38563576125173815 s2 = 5.3185967578488862E-002 +: a = 1.0000000000000000 E = -0.50000000000000000 s2 = 0.0000000000000000 +: a = 1.5000000000000000 E = -0.39242967082602065 s2 = 0.31449670909180444 +: a = 2.0000000000000000 E = -8.0869806678448772E-002 s2 = 1.8068814270851303 + + +* Variational Monte Carlo + + Instead of computing the average energy as a numerical integration + on a grid, we will do a Monte Carlo sampling, which is an extremely + efficient method to compute integrals in large dimensions. + + Moreover, a Monte Carlo sampling will alow us to remove the bias due + to the discretization of space, and compute a statistical confidence + interval. + +** Computation of the statistical error + +To compute the statistical error, you need to perform $M$ +independent Monte Carlo calculations. You will obtain $M$ different +estimates of the energy, which are expected to have a Gaussian +distribution by the central limit theorem. + +The estimate of the energy is + +$$ +E = \frac{1}{M} \sum_{i=1}^M E_M +$$ + +The variance of the average energies can be computed as + +$$ +\sigma^2 = \frac{1}{M-1} \sum_{i=1}^{M} (E_M - E)^2 +$$ + +And the confidence interval is given by + +$$ +E \pm \delta E, \text{ where } \delta E = \frac{\sigma}{\sqrt{M}} +$$ + +Write a function returning the average and statistical error of an +input array. + +#+BEGIN_SRC python :results output +def ave_error(arr): + M = len(arr) + assert (M>1) + average = sum(arr)/M + variance = 1./(M-1) * sum( [ (x - average)**2 for x in arr ] ) + return (average, sqrt(variance/M)) +#+END_SRC + +#+RESULTS: + +** Uniform sampling in the box + + Write a function to perform a Monte Carlo calculation of the + average energy. At every Monte Carlo step, + + - Draw 3 uniform random numbers in the interval $(-5,-5,-5) \le + (x,y,z) \le (5,5,5)$ + - Compute $\Psi^2 \times E_L$ at this point and accumulate the + result in E + - Compute $\Psi^2$ at this point and accumulate the result in N + + Once all the steps have been computed, return the average energy + computed on the Monte Carlo calculation. + + Then, write a loop to perform 30 Monte Carlo runs, and compute the + average energy and the associated statistical error. + + Compute the energy of the wave function with $a=0.9$. + + #+BEGIN_SRC python +def MonteCarlo(a, nmax): + E = 0. + N = 0. + for istep in range(nmax): + r = np.random.uniform(-5., 5., (3)) + w = psi(a,r) + w = w*w + N += w + E += w * e_loc(a,r) + return E/N + #+END_SRC + + #+BEGIN_SRC python +a = 0.9 +nmax = 100000 +X = [MonteCarlo(a,nmax) for i in range(30)] +E, deltaE = ave_error(X) +print(f"E = {E} +/- {deltaE}") + #+END_SRC + + #+RESULTS: + : E = -0.4952626284319677 +/- 0.0006877988969872546 + +** Gaussian sampling + + We will now improve the sampling and allow to sample in the whole + 3D space, correcting the bias related to the sampling in the box. + + Instead of drawing uniform random numbers, we will draw Gaussian + random numbers centered on 0 and with a variance of 1. Now the + equation for the energy is changed into + + \[ + E = \frac{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r})\right]^2}{P(\mathbf{r})}\, \frac{\hat{H} \Psi(\mathbf{r})}{\Psi(\mathbf{r})}\,d\mathbf{r}}{\int P(\mathbf{r}) \frac{\left[\Psi(\mathbf{r}) \right]^2}{P(\mathbf{r})} d\mathbf{r}} + \] + with + \[ + P(\mathbf{r}) = \frac{1}{(2 \pi)^{3/2}}\exp\left( -\frac{\mathbf{r}^2}{2} \right) + \] + + As the coordinates are drawn with probability $P(\mathbf{r})$, the + average energy can be computed as + + $$ + E \approx \frac{\sum_i w_i E_L(\mathbf{r}_i)}{\sum_i w_i}, \;\; + w_i = \frac{\left[\Psi(\mathbf{r}_i)\right]^2}{P(\mathbf{r})} \delta x\, \delta y\, \delta z + $$ + + #+BEGIN_SRC python +norm_gauss = 1./(2.*np.pi)**(1.5) +def gaussian(r): + return norm_gauss * np.exp(-np.dot(r,r)*0.5) + #+END_SRC + + #+RESULTS: + + #+BEGIN_SRC python +def MonteCarlo(a,nmax): + E = 0. + N = 0. + for istep in range(nmax): + r = np.random.normal(loc=0., scale=1.0, size=(3)) + w = psi(a,r) + w = w*w / gaussian(r) + N += w + E += w * e_loc(a,r) + return E/N + #+END_SRC + + #+RESULTS: + + #+BEGIN_SRC python :results output +a = 0.9 +nmax = 100000 +X = [MonteCarlo(a,nmax) for i in range(30)] +E, deltaE = ave_error(X) +print(f"E = {E} +/- {deltaE}") + #+END_SRC + + #+RESULTS: + : E = -0.4952488228427792 +/- 0.00011913174676540714 + +** Sampling with $\Psi^2$ + + We will now use the square of the wave function to make the sampling: + + \[ + P(\mathbf{r}) = \left[\Psi(\mathbf{r})\right]^2 + \] + + Now, the expression for the energy will be simplified to the + average of the local energies, each with a weight of 1. + + $$ + E \approx \frac{1}{M}\sum_{i=1}^M E_L(\mathbf{r}_i)} + $$ + + To generate the probability density $\Psi^2$, we can use a drifted + diffusion scheme: + + \[ + \mathbf{r}_{n+1} = \mathbf{r}_{n} + \tau \frac{\nabla + \Psi(r)}{\Psi(r)} + \eta \sqrt{\tau} + \] + + where $\eta$ is a normally-distributed Gaussian random number. + + + First, write a function to compute the drift vector $\frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})}$. + + #+BEGIN_SRC python +def drift(a,r): + ar_inv = -a/np.sqrt(np.dot(r,r)) + return r * ar_inv + #+END_SRC + + #+RESULTS: + + #+BEGIN_SRC python +def MonteCarlo(a,tau,nmax): + E = 0. + N = 0. + sq_tau = sqrt(tau) + r_old = np.random.normal(loc=0., scale=1.0, size=(3)) + d_old = drift(a,r_old) + d2_old = np.dot(d_old,d_old) + psi_old = psi(a,r_old) + for istep in range(nmax): + eta = np.random.normal(loc=0., scale=1.0, size=(3)) + r_new = r_old + tau * d_old + sq_tau * eta + d_new = drift(a,r_new) + d2_new = np.dot(d_new,d_new) + psi_new = psi(a,r_new) + # Metropolis + prod = np.dot((d_new + d_old), (r_new - r_old)) + argexpo = 0.5 * (d2_new - d2_old)*tau + prod + q = psi_new / psi_old + q = np.exp(-argexpo) * q*q + if np.random.uniform() < q: + r_old = r_new + d_old = d_new + d2_old = d2_new + psi_old = psi_new + N += 1. + E += e_loc(a,r_old) + return E/N + #+END_SRC + + #+RESULTS: + + #+BEGIN_SRC python :results output +nmax = 100000 +tau = 0.1 +X = [MonteCarlo(a,tau,nmax) for i in range(30)] +E, deltaE = ave_error(X) +print(f"E = {E} +/- {deltaE}") + #+END_SRC + + #+RESULTS: + : E = -0.4951783346213532 +/- 0.00022067316984271938 + + +* Diffusion Monte Carlo + + We will now consider the H_2 molecule in a minimal basis composed of the + $1s$ orbitals of the hydrogen atoms: + + $$ + \Psi(\mathbf{r}_1, \mathbf{r}_2) = + \exp(-(\mathbf{r}_1 - \mathbf{R}_A)) + + $$ + where $\mathbf{r}_1$ and $\mathbf{r}_2$ denote the electron + coordinates and \mathbf{R}_A$ and $\mathbf{R}_B$ the coordinates of + the nuclei. + + diff --git a/energy_hydrogen.f90 b/energy_hydrogen.f90 new file mode 100644 index 0000000..8b3fe7f --- /dev/null +++ b/energy_hydrogen.f90 @@ -0,0 +1,39 @@ +program energy_hydrogen + implicit none + double precision, external :: e_loc, psi + double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm + integer :: i, k, l, j + + a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) + + dx = 10.d0/(size(x)-1) + do i=1,size(x) + x(i) = -5.d0 + (i-1)*dx + end do + + delta = dx**3 + + r(:) = 0.d0 + + do j=1,size(a) + energy = 0.d0 + norm = 0.d0 + do i=1,size(x) + r(1) = x(i) + do k=1,size(x) + r(2) = x(k) + do l=1,size(x) + r(3) = x(l) + w = psi(a(j),r) + w = w * w * delta + + energy = energy + w * e_loc(a(j), r) + norm = norm + w + end do + end do + end do + energy = energy / norm + print *, 'a = ', a(j), ' E = ', energy + end do + +end program energy_hydrogen diff --git a/energy_hydrogen.py b/energy_hydrogen.py new file mode 100644 index 0000000..c6342b7 --- /dev/null +++ b/energy_hydrogen.py @@ -0,0 +1,23 @@ +import numpy as np +from hydrogen import e_loc, psi + +interval = np.linspace(-5,5,num=50) +delta = (interval[1]-interval[0])**3 + +r = np.array([0.,0.,0.]) + +for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: + E = 0. + norm = 0. + for x in interval: + r[0] = x + for y in interval: + r[1] = y + for z in interval: + r[2] = z + w = psi(a,r) + w = w * w * delta + E += w * e_loc(a,r) + norm += w + E = E / norm + print(f"a = {a} \t E = {E}") diff --git a/hydrogen.f90 b/hydrogen.f90 new file mode 100644 index 0000000..eccca10 --- /dev/null +++ b/hydrogen.f90 @@ -0,0 +1,25 @@ +double precision function potential(r) + implicit none + double precision, intent(in) :: r(3) + potential = -1.d0 / dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) +end function potential + +double precision function psi(a, r) + implicit none + double precision, intent(in) :: a, r(3) + psi = dexp(-a * dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) )) +end function psi + +double precision function kinetic(a,r) + implicit none + double precision, intent(in) :: a, r(3) + kinetic = -0.5d0 * (a*a - (2.d0*a) / & + dsqrt( r(1)*r(1) + r(2)*r(2) + r(3)*r(3) ) ) +end function kinetic + +double precision function e_loc(a,r) + implicit none + double precision, intent(in) :: a, r(3) + double precision, external :: kinetic, potential + e_loc = kinetic(a,r) + potential(r) +end function e_loc diff --git a/hydrogen.py b/hydrogen.py new file mode 100644 index 0000000..ff8fa96 --- /dev/null +++ b/hydrogen.py @@ -0,0 +1,13 @@ +import numpy as np + +def potential(r): + return -1. / np.sqrt(np.dot(r,r)) + +def psi(a, r): + return np.exp(-a*np.sqrt(np.dot(r,r))) + +def kinetic(a,r): + return -0.5 * (a**2 - (2.*a)/np.sqrt(np.dot(r,r))) + +def e_loc(a,r): + return kinetic(a,r) + potential(r) diff --git a/plot.png b/plot.png new file mode 100644 index 0000000000000000000000000000000000000000..3566478c10981d3a43b432c1af3e8f25dd0168ac GIT binary patch literal 7938 zcmZX32{=^$*ZAD)4r5R=jGZv{v8ycE#+IE4N!CeH3XyCT?nv2}N>M3O%2L@<_T^5p zeG!$4B2yHVP!h7td;9+0_y7Fg=bh)y%=0<-+_T^FIp=)RXwEhQcnLfPgAuT^wQ|K^ zKp2Ao*f=f}0&Z1*jV{b-j_%eR4hMxWnasa`2nGXl5DW)1H-}I$Jv|un;13)`76Vfe zn4?XHGcXt`jN!nXM9%PC3I%F@9V z^yw4IVnJ6ocP_VLH=e|)V^z(~ao%81ZkRyyH<6q{*38)V?>a+VJRB7T zgBeJ9_6_q_L=NWJa}MbeCy}-Q!yIeOANa(VE2kDdXSUWKYdFuhv^2(qWjS{-b4hb~ zn3vyUn2KI7Oh;fk2mUKRmBU0YI)@p~K@bkA4+j~jzq!C!VAc(= zIIIDr`V5E4pDlCWO1?;YYaj0MQ1>}n-J#*zTMsh*l3n^ z#06YiduBu;5c@6c`@%}+IhF{Sill8ui(gt^I~NXdy6ENQ?ub&EM$%yZySM_NeHw1S zm0x=cNliCL9xmkCqMY3ExOdt+T3(i1aIWFXjpoV5cJU^0zUSS~YTl`e;u;IG5twdrkK|KuOi`f ziiy8@MNtH~r0?L-Z1%#Hb-DA4Tl*l;0@?aOV|8)ko|KAkHI|Y^pTDx+aUo6%gz}1A zfmZLIIpa}FTyZMGs7b)(E9;lUyR+m2u1 z0X~&EBSTl#`^sh58u#3hAablC^CNCcNN=L@qGsKcB~_$a93V>l(X=ME@DgbJI9u~_ z*eHXp%0*-mwB7l^*M3wDBC%zX|2!N(rQ+?eR!}X{tjN(DXCB0Y_QH`N=gS@bYW^-( z8G2VkVbn=F#6eZ{Q;}T2-7*@2Tv<2aO~nNA);q}B zU5q9FUSqW?et#jAs+tNOZ4#GCu-{yB@ZaqG*TEZB)n#zM_s^}9$|GBQ1fVm-w-?Q> z#aKZ5EaHh$BNsXDwcXaXCd10~8P(Me)sNVDzRx!15TAJJcTa-0CqL-nM^)QHe4ubv zAB}eau=|ez$zmXPhB+ejL^6oT0thT|f?ogJUZl6|E7n+CZE#=IH#k|K_H-6KAC9sP zPi&0G7CPzs1oGOBYV!eHo$h?Gtr8sDrnb^4H@$vx5|>Ik`_+|oMh<+XLWe)FRq9ZO zI>gVkMg9T@Q@!+fbP*MH&(RKbvH{J}8@lfOlgAYF+I=19!s)dIiipKVKXn=Ic9~i} z9KUYcGB>$2ZL-q3{A19D^}EtENSz}`AOilEZS+t0oqu>hgeXPb169>cZH34y?PBUX zx91-~6_KN3Ou)+esjWxaaJNDVIt}T}k@(*mEl(f!KUCy|Yv&*6|8kzM#5-b6pmt3p z+yto1;B{t1Xw-)Enwe>(<#e3=kpjG~Csf>qUCZi|b{<+d<*@mEgT5XSWjRq{Ab`zt zErz}S^Dj>a*fWCD(ymS^<0s!Y=ct;){^)x_ii3Bu+n0E`tYKgM!F!_R;$h~dP zYPv=)wwIPp?Q@X_pVdvxe^Fh1WmD~jG?8Ol6^kTMT4pQtQ}g!7*K*yY4ey^ZOL^;G z%Uu!T^p$mA6GfSS^FX2ZPyg8|)zwhkk?WYy&ZqpDAy*R}Hh-C9Ec=WdtR8l!(3gfb zD_-*|ivJCb0A-ins+$s6{?Oj~_AkmpcR_N}e){;u_WPxOL!ay|@;@?V{Tu$=I9sc7 zIF&9^EzJkIrkf+z_+|o!(yk#@^H$J4FRC&jH`Nt5mP23OnfkaO55>gUQcCtAlbs`s zb2gFgmcx|Hk2uq)t}W`ug=O2kh~XRa;bo`n+k5~te@J=%Oc6`isCdHTo-RCc3l`+w zUDw~F$f>p7`@!w!UJ}i2jFNee-!H7!OsZ5G#ur_zx((Y?M8X>tIU`QDTGRmVRC+V~ zJHP_eMLiauTx+u+&+XuZnr>pGZc4IGv$+XZeuE-Ef_nStvlt1K028*qM;b!cx~5x7 zCMMi|5?Sg~VZDhdKQ?b(+zRYRf|8|e!X3e!iPkviI0%sd&~*=m&PALh4DE5;{Wq~G zh@6FAi457!QuxUYA3+HQ9GHjqjDt=?A^J%#)6J`ou+YY!RBn6;u7V z$L8&ru7Fz47MX2;Z>JUXfGH^_3ynQ0K!H3&M{dv_1N5SqE@^gwfm7lF+k|z`B*O!w zQ4ycnoUoCD(3)lvkrgS1ePq_c4GL!1J?KC;zfAD0%*5kvoU?DfX;fYxTU4m_xsEmc*d226tP2vUUzUmbs5Qip%bMn|7D?ll)52?@zWp(0@du_2bmCIZHYolZEKI*qe`UH&| zX%+y_`ks^%4bWn2kb{=_u^N&R&ETF52%5rYmry=bybA5gLcWd*eeNiW?*qz zlh}-ipm47~rKK&y*iR1v1aJoQo8lRwi?mkw^U}@O8sIKivM{&3=A|&mudnM_WkFtv z_D(wUWrqSAssC12Q|$~4v{6(iU{_$f>+_Y3jbf+sXNZm?t}@x|O#PkKqSuRdO*_Rq zS%;SheNKvs)eR(R8w2MvbetJ&0d7?VeElan$5e!KfrF1vtrYiMmFx2p`N6*#eXgcU z-RGrUGq}xT%QnAwt$jJuPN9Y3El*}J*x%B__?|22pIaxtIH>WCrtg6ith9^Z#em3F z)8}e#nfBknYwc`J-ept`;% zVHpi(jde3#QERi_&Ki{o(*m>__p%P0X-L1IM77XqQ(ZXxTQcZtw>#ZK{Yz=I0(<@B z(2DLG<1!~Ct~t46y*ObnsxzJ zTN`=fGeV;^R6g5G4wY4I`LoLxI_7;s>coV%%p1BU0~(b<@jjtkl)AbkrSTzB{sK7Qt<~e z*#=x1f4=y}UcL0K`)|?wk~7i75IPP-%*NR4CX>cosr3mw-~Sw1i}+qSeeg}YY2%GQ z%7?Q->mrRlCc97M8BR%Z=4Dk`;&rQqPb=$#UxLm&abqrYWX0nz`UMyPt$+3|lzle_ z2GOeehF|g?@oas6=(`JL1OvQ+DO0 z|02;?G<>Qdw(}I^*e*6|Mg6DKrVsqxk=|#VR9ZEtb+(yvZQ68Wm43Q&S#7*^;uX{j z-b@kBU$yS(iBhE>WrVIcCf#q37-Y<+rLRd1>0_HF>#tmDPa8W_aA0&*PopFA^Np`r z>47h%p zHSZnK@z^GD`&?gy@CasEPvc!o_o~g)Sk4%D=XSK9U(wBGvIO?smp(CFDMRi@bW0S_a>@E{a}%Nnt*mjm`SC4_W$0dq4F=p9d(x zy(w6&@AgzUM0P@@!C8ql!=_UUi-2r(aqJ-X#i4eb~3Ilh(c!k#}?S!cR~gb;Gm zgPC!L+M{wR%((JMj+CX^U!DD<&$w`f8d-e-Uw8)*!yh=A0GX67svGe^_4IpgO_3?T zgdW}o&-u3QmC_o?d=+7u+{;|=2?d1<$jke?-b1f!F_s9fOQ+p?xf4Zve>R*Wd2_?SGL%v0K=u~$mqf5cHdU-~o+PFwTX857IQTxza3L5+ZXg4bO4)-X zPe-Io|A|u;1HNpvaO1=$ok*^?e!3+WCr+i}4tYSEAvm%4n ze#9MGv+@&orCv>p=LNkTK~e+N9w0_-DT;CjXJyqwjZc=L~Z_J_y;K`(8H|qO@{{< z7NNViOaGGSLj4^Em!i?7cR4SJXWFNiJ&r=&5Wl3DdYviOQB?bz`>wiU!d45na^Yny zyQ*z_`+HBvYlh)7g%*zAyKqFd?*!#>{dK`1ur_M-4zF0U8#N5-qBg=`A8bck8*U%@ z4h8x|)Ohfj*2^&0;OJ_~dPLFVtj)+?88bVxWiifre8T!we@=d!acxoJpXEiPO-?b@ z>`sD0&YKzT9C*5IelQvc)~a4P7ho87>!Qc(PE~6c#;t(YyxLKC_4irJod)MecO>oP z2}$k6{pOFxm1!qi4(-n0i%c7jNUyzjWrv-sjSl7q3H<1CGoO0$yamd_Y68BR~cuB#))R9 zUWw0ggC!qCwEL*jsAoKlT=75I+Csn8!W znwlffKD3D!S5U~>bVKyEi-BFWzu@l*Lq*7ye+1AD{MG0Iq1-}NTL1&f7Nte$x54o9 zP7jehA??@HfpG&g{|-D*^bFFrf!+x=bfUwX8c!)tX9hp`RvF&CvX*+CBj5Ix>tqh(lR0ER1%Cu{Ny_U;*(m zcWTxsOJ3Fia>EXfWbw($KK7*dCtk(=`eU<`c0gYfC2vsxos^fQW2sc~En?B(@S89T zQj&7fK=tsAhPtnx5}Z*NiNxFPBKzy`qMSobY{WtXxP`F%e{6Osg^nLV-?z7Hp`o{hnNx?5UxeHzbnXyg_zH$I2X7;!3@uE66!o8g&@s0)Zvml2+MouglNZzT-v zQxA(CwWSX3p?ujbV~-s+4{SPA(GiIAjKW(|^=l@ZuEIMdQNaa9i6-}I@8?e8vRk<0 zJ{WiijEiBwru5E1%AG#CP2WYWB;{F~`%Fno72)vHg032erIq#}N8QB0w{{YPnU|0U zZfI&pLxuf%^n%_s5}-;fRA7ED_*)VYhCI|G^hQ0XeS0jCUGfGy2W<0M(LXhll#e&*Q4LZ!*U>t1~9uJzH2qHeQW40cWDuv}2J5u0R*5Q7s zid7xeyq>=NmWcS(;`}p$J$Z30TtTXp2oslHNmNa-UuG;ALM(o#^IcJkM*0$nMm(m&W5KA_B zgt}BtFOjKgK0Y`sU1{>?-Q8xNI6hmlj^>9xx}=~eODSEzy-EBQ73v$)5x~!6(Zc|< zPJ72B?;Bj{B!cLW)CFmrrv_aZv=kTWS7z##nf)7``Zf%aeKo-H&P7NpFe~ zFZGk+A!Du)lG{j)3c7o}gy*Stj; z$oSQg07R*xA<4i+FaN0|ubfLxaCCR4weSX!S|6 z?&CJOW6eh`=X9$mWY!tmMt zjO=!CP}c?K)o_TZRs}B6s>`(5qC-iAa4_k6cQiD7U_x$_64-o*s{d)1TRB|wa+Oe&5PcAUXzulZ_bV2%E9=u~)j>a2rpGO3i%+uxuFJR}r<=d)&KQ6k? z71(%b7?~J8zwyT@n;yn7@D`c^}v4Ch}gl3tXuF zP%U4Bul~iIWA4=DDn2Oiz_&1k{*>#u5$BP&qsJkm1^(Crz&nw70q7)MeA5E#dKTw> z*k+p@QN3nMMXlsv-5$Pz9n)WC4Jl~*aO8aYefDeHbawjO(+Cj2+LFR^{#N@j(aE-% z^n|~|JeP+yD}ERB$dxv2Q)Z;HG5?-fIhMnhCyC)$a>LL6 z89tgvucW5ZpG(RKpOlB(<(UZd-}6T@Id+PxV2^zuWu+~S zQvo0av>rsKX*}cteu$*{n-54^z{j?n%r6dPc&#}8<@+F^u<2Oh>}PHZsu>~oLd0*Q zaWTAmuC(76xQp!B)sAnvkeIZ~PYh4Gp^r^(#QKN`p>wscDYK{V(r%$AYtvZppUu3K zgE2uWH1okY2K>;mxl;1>^6*k0h*@D17si%1;tDudHUV0; z+B%o_3Bx5WZyD|p)s4Kb(+l{?ONnnO)L7GomX?C4T5J7%j{xpun!Bx~>5GcpW9RpY!54Y3e- zNAIx`PY89W3tBP-<_NS8E-qJWf$N;il|_rOqs5`rKhT?!IP$IvC6LkYqR zIl$1}`P~E0`|iDuWAFWc|9|`(Gxz{&ty%Yd)p?%RwVt5oDsm)OZ(W5zAS4R%vM(SI z+<6EDXZs31_)Q)g0T=j3)ai+)le+C2Cs!i}Gl;U0lbw~Vla+-rql=k?qlK-F0Ji`) zKPTf`Cnq~cF&-Z4e;>eY>tN2q^NE}YT!he0Uds^zAv40hak8W{Eg+C&V+C26m(axJ zF*lFf&NXKn!(Ah)yOLu%SFOJ-JpK40)%OMal`B{DU1ghRSwH#ELT{>7JLG2tQsrOK zf0)n2^y`@x5u4)3)gMv){n7A+XhDTy%#4*lo%CFnf=Yb4MZyg8!cn^keX*zc9IU@T z65&xdffz2ZMn09j_$un?a=HA^f6L-j2A+Sl1*SsuE`Di>_uPk@>HLfT65B&6^7Ahm zoExdZ7srM`uJ%5;IQA*z^8Cf2Om}@QIb9qYPX6%Um%Z};c-iWayn4F=&S1SNm++XF z03+&=ywHC>M+R;^Fgn_?zk%lK#I<1EKUy@yR@`)XX!sTQ<(MR$47-omLDd)xXzk*^EPF>$HyN3x^4b+QkS=-5lun= z>(c6t`0y`wXU6V}5=&4kdL*wUc6RGv9V+qI#rt0uOG~sy4h|0X4h}}G{0wpX`0?Y* z*RScFn95K8wT(2hsLAT;>O?yc3FZd$R1e%*Oib)=u{jw8-SU}s+x*|hHAl3Yiiz;B zv$wlAx0F0*!>do0gt77PFl?``KBHj0u+FCnlb*&YyJ$-XsGWDuW4AD{A%@9k5{Lwi;KfXGnIgO4_&oB4<^ZVykM}H z=U?ZW4kC-ph|H6mDc4{VN18S6LNNwjqA|J-gkh4qnkmD~kqVE#U5n9kp;3(CzTGmq ztMV@wFl|l?lj>+^O;Iv3x^wT|y+qI55v%d)gVa^$POdUu>rp|k)2$T8?LN7!Klxr? ziFm5N(AriPA}k#o9QM0yy&9=%*Dv4a;2_&IA-?smosPJ7v9b9{}-;z zoZQt>{8$QsV2X($UeOS!Q$l)vH&W2DQ}_3tcF;cAR(b-eqNFwTT$0 zsmUlR5{1h%r!C9>_$O#giG!LV=^T|O(0Kve3753<-c@=ZF^xYTkEiNq$ zaAE%Z`SbGCD?YQ1sBb?P1Jki5&+wLP)=C#!cP~D!Rh4nOoawz_39YVuoF5e!NaxZc zilNI;iu-L{Hr4{dppY~8=g<4bh18Qy_7iYxos?l7t)dTd^lo@Q9Uaf!Sy~G83=a)O zq@+mS*G9--EzoAfI@;Gmw0CUre6`f^M2F6EGg~RibvYHKDl1E>*M>0~3yU^e8NvP2 z+gtBXLR-E9M`{+C5<(D?d%92Ka2KX)p&Q20JO;`3r$MaAh8xq(xOwVuc5GS@?J9XN zh&xI4QEPX$lk`xEN7{*@bj>E++4jWnm9nfX{(B>8g28SHoll=Wv9YsH4e8jLp~~%F zYG^qAd8J_%CHjSX)!=k9rs7~O;b6C{t{Wq@Cz|ZO$i7;Ep{MuSpK^*zu&s?FpojY0 z#_hUIVkRwl>>Y`~&Gqkv{t(HPJwQ;QbDS1b^1(3mkVZEl~j#k%QE= zV(0Fq!P=^DxW9kH&v>w%Pn8B{R;$+LiaNQ|e*XLvCt$w$%i!dbx~k+(>Z&mCMDU?T zdF9hrZ7LU7SX2_**1Yz_qq&wAm7g0J)FA|I7>*|i1SL%@B!)b~EwkwO^2R!Vijo2C_RGv?llptWkGPp#YR;9!o_ zV?P_t)9WOVTE~&v3=i9OQ}*(bqJ3Z`zs4Xx>IITs`(X6d-YrtThX@Y#?(!6n*y!TH zjdpf)@b8{&NuNb73}iZXuxD<&eKw*735wMD56GpYL^i>Dcgj*T^IWUCK%#H^##el> z4z`=&Y+p`luL7}yygxoP)Dp$69xG@;fk}pO+t%z1@eaL^l{K}pP;|BDkxV}-=UUJG z3?r*h0oI6^g`6}HFr;T^2LG46LLkBy``|Kp@>(>qB#3>wMnOAq}e)a0r{?2rTNdNDL;L2(JhWCPq zJ~*!9U-D`B1{O)|pL2O03AaXr^fVKp$Ss<>+Adg^i6052!KH zgNvU6149TmG^2XpuOXYHg~3%u)C;ilj4kS=a$fDJGq-H05sx&f8vY zk3QtQeg&@$fk0@2ZBFO-vWtHS&%uH8z{~C4vl}3j^J^ASr0PUB-*!bcm`iQZRF8Mb zZc^zvk!u)uI=i3ANgtAjNuP?W)gE!@*X-U(IULOV3%5q%lvRTj$j{xd^y9O(7}lMC zCxea2myJ^0x9THRr4*cXANl z>T;;)TJCZ(ArV>Vjs7g`0N$jJ zbHQ3#Y8#$i!#J!B{N!E*yBomyOA7^Xfhp8^ZCZ-DUt#*zA77AaNA|h|5O&LjtrS=A zL^j**R&4c%@9CVP91+-TFf%j5>wW5Rx}7OaqG|^Fc99dO;Z7L(X}WoCS8nHfLD$K{ zt^W7&Q@#G1s4Gc4=4iDlOvrL*eFWhZ9n4)B=zWZ4Xw=-|)>@dxCjF5qxD%NR)!uFf zOO7eU0GcN1=u=?Z)Ks4w!)&YcDWBzVUNd^50AxID2FOrN@v7geoQ>9rePDe{`s@IK zBXT&O<_6MJ`9=$^7TiiJf7oNMYNGq^83_FzbeMF{$=(EF#5~1wl41VZd4$jKg0(O# z(#_}DdEsQV`dRCg2n&&RNIsFBS??Pk@A;~s-~KDVw&T-H4qoNRErc{fXAvY;684~? zqQaBgJlXBdp+1L3p6NpW51Ks#xD^ffK_XR5NLW~1SJz{nVU;lpn@DDOjqXK)E93YL z``W~^h(2j985-^zeV($q2^OwpioCYg)>G_hB;sQs!Balz8*+lQtxGB@;vdS@k$d&( zW@ii!D}C6WnzeN>?K(ZGSD%|)VJ5osg=&5 z0qOM><DM;Ml(lDKy8EOW*Nb0D9C5nJy^bL+-v;bdJS)sMUKD3-;-G! zU+qx*Jm84Yq#_UJ-1?Xly=)~3^SO2F*53jn+_EtjTE%EioXHAIawnT3qaKQiN55-x zwo{--dZXz&r+#($ggX{wVX@Bfu=fd-!SOT|f$rJ2b2L1|+b2YDI2n2!;f#5t0kAed zKXnU^vj8Ic`*)^i&z_a<4C&Aa+t2~DpzHkS6^$8bdhk(9-LaWX?ZKPclNH;o3WR59 zX0rPy*lI;U0E3r^-rQV*?c93;dc)>VH+=?jPnL@jb4vziVpO_zI3V01`zNP4%%9XL zU9oHYfc!1_mur6YMEg%Owg&lF#G^d-#&WbvKLK=cUjNFc{QYVtN-PGT`%?jKAWO#z z+tOnr;9zXcD~7LG+tX7rLD*Jh?d)`~E;B97Zz;bn*&f~UBZka-zcFBYgfvv7HdjWtXD@y%n#+RZ1LTd$+id5J%C1n$|pC%!YN?sl>>d?ew$QabDa-Whvq)mti9 zNb2V`Zm3Rlj$wXVPxaiM$?n~(ycZ-MW$+VQuU!A}zsFyhxihn__?R-uq?ejKU*iIw zBBF+6Fb>oO*08pSCbw#--g|}hmVS0sg_B`9h)83P?*N9;~&TKK^93mne?i2&7Z@2hu(jV_m@&g zbezy+epHfXBA730QM*?Y_UD&|RC^T`ue?J0V zD{{14JW}npB{Q}Ztz%;}T4Hg0jFHC7r-jV{kaeF*A9_c49E$eJ=VE971`ru%_Z~slah$O!Sk!L0^PPzutUtug_sHaSeU9d-=78_c>nY)C(&2NS}Fdc3v{v zT^^a4o7=6xLV*QXu5Lxd55<@_DDROXvmD(Dhq|5p-T*0CH#bqR9)9m56H&C*aV z7A($o##|3;J@V3LiyF_%&i)2M747ll$z?7sE?K7VkfJ;kO81<|^eeC$#k1^z}A#12X)6BOKvCf@a0d;Oq`XLYOj+|`dRgvmT953 zxu0o;9lC{v*ZUQD)C>$#Y((?(@-o0pu&ZYfzj}qk#>F+;!I7`D#^8#wUmh(DrnM=* z!+}#-Sq@w^#NPUGV#C;Ms4c(*U}v9W5*qz(ESauku?*Ye#r@2joW{I5pGgU zfwEg7=(GI&@gG?It;Z(qhzJ_S5@=(?8A32ZT9_Wx%pgdBtP)Vj7;Y}jxKvRZfMWdqoq z7pT}IzgZ>h&J##WIlKn;eoB#o%EE$$JPh64e?tiDnpkW``1n{T3=29%W&MZzaeR$V!D;GB9HwY6R;7O1?cRWX3l{D+s! zhH|uyXL-D*eiQLXM1{HZ5+mU&ITCKvr6ncLpFe--=qPY-aB!}Av=tnt@OVVXD+;8X z<{PLQk@r zfLpNp3W)&AQAY>129*fnKA!Nkp?-QubB;6}CeM>fh|x+mo;`g^*z@t;Jtp<Fr%D=fCO@4%I5_l&jLjodP@ zWa%6RniA8lGTTr6endRb1m_fskN2OA1L14(@gj~oiXuvcS1(sqCX)#bY=2RINzA5( zUHbcrJ9+9XC}gBV>2DOF{#?)s(arG!yHFahD%Z8S>nZzq0&luNLgB7E-o()Z2)SF! zJn4xDT?PliiuO+3O)n*Vd7|f z4^+o-pF9PY$VBHMttMkWEUxxCJDP9=DF7SD#P-&(iwzFBwTfrYz2|2A8S8!N-eFRQ zY}gzHw4$h}D7)Pe^J2HHIls!`!EE)4g^Yyq`R~DOY;4rMLRby}auFKKoVG9s3;@+L zl}mZ(;kr|yLe1sWMSHe?-|ep5y+Fqbgd@ffz?UhI-d@Yey{1IM$+GhPuXGdoe#E=k zj;~#D0s(-#Zx@<%M)&pi|JJdsahBhcq8>N${tjp*F;UrNSUD(tAeZG?5?2PwyZjBm zLK~_V`_&xbVSV$dh{XSzORJ@BT)pk2hgDYieH|HV4Aj*fv6sJyOEz%v{i^az?cT3g zm)kV~3Yo5Bw@~$VzvAHVaHdYruV25ifBrP?h+@~un_BZK%G0Bh^I61zXfg#^D(pBkGv@}-mo1LA7NC2=ngTXj<@wL+R*O0Xoum)A_j#?(w@s34n z6~EbmbsL-=YP{)=f5J;`+LI^>)F_YGcO2OCGy+%Fi3ZT-#OsmU!8x>%pZHHQ3!twO z09s)cG8!|ifD|rAKZ4;m`u!nf^Or{c8z7#Pl$NS#XtXZ=`GwUV4v%rkB)*V)Cwm<( z(19EBI2=fs`z$;@S=dUC28thjjJ`Yj`-g^ay&*Mi&nWk>n+1x>cK;Taf-6t;`XwL} zv>t$hj{r!p#|@LzihYlkS_TxhfS65F1RE7pJ|JmVj_dRs8;B#)1A-NRu4$FNy9tE0 z^HgfrND|-)&hC2C5!<q*!tcI)=c#G^FtXPf{J-Vo;8qzW*lS@VZNkukd1Hi7k3-*bAcL3*Hv9? z{RWJx94IvVz#h_Ui1O3JrGqIlZxOaV3b5qQV9BWO%kb1P3lOt`Z9Zl{f?=*80d>Xl z`6(wn*m17~@Jx~{=@%RhpWs^tGK5ffclUzX(e8?$wmB#k;oM`*HOz`Ql?lu$=ae6x za|y>?M~4ZN^t99(&3}j>q^#Cj@wj7OSEt9b&ZPrf zp=F#y;kpCm(*uNek@M1PYy!n%Jqda;tYhI|(_lGJV+LmHPIvOR0I1yDN;ym1SsHc( zoMChESNU9DW=wyD{9J=dMkTC+Fa=|?6d|LWBrI}I#c8EkWj zD1RstB7Ngh6xK~lQ~f0sPyi^5tsP*wdWs3dL7Uuuo17HCtPpmyb*7V&WRb}C%cZG) zk*PMerUce7S~hxVKFghJ;Llt|Qs<2Yd=1N=&7)MrQ)C~DUz6SYbR8%;ycc;eHV@DK z9<_gXsry?17s+_4(uUqzN_atqqH(GEjcSFTj zCqo2!^IlS>|<;qUE+rWh-{s zpkCnrrby{YANGYoK}lB?p_|)iWx#M499#uOLDyx(EWl$8_aC05z#7ESVdF#J(buq<#=w2kjv`gNS59_(!?t2^ckdZfQZ=B)6Ch{LEmT7r;p@5pA_l|%zX##;xz9+y2zgC@>Bm`{IM0RfO z?1$?rjMT)NAibC6BCjV%jNfeg288qeZ7;VF#mPJvR> zl_(l%aI&b8=)Rau;E63a>17)K6dSOCe=D6iARym?5`IUo1+%pQbd-3Gyz)y(xF}FB z51M(r8bC|xqhDU*xAR6Rzkf&+$i`*#xv;%A6)YeDQE=&1MmaAHnbC(QVT zkGQBlGec*}Y>n@&$hg@UJ&6q}K2&)&{FV89al}bSYjr6T$8*9t`Ze)y1v`L(SOLHk zR;8^Q&ebaRSLCUQ0&2oXP~2*FOY=8YY`tNz?~)YG+ho*0l5;=9IZWo_(o{%zo*YIK ze4p&lO?=!Y-N>A_N?3T0HWd>M*l)3Q)x?ksSb^V8zV}Hmna3pV@eHfv_)wa7b=o4n9c!`0QC=iitrU3MqaZEX)1GUy;BiMj$PW2*Ay55deQT_Ff zC?C!3hA1Vz>T0K=ckV^N(Ag^zu861a%HZ2NJuB1nVJ;LaZ&_m_VZDeO&{NrJn`qf0jvvbHb!|Jt?Zw8y+I^*gOJ0zbof_@-nkyt0ztwF z5}ET{g0DZ z=NE-Wjrg}{Y02h`iSaLF2b{i8m+3y0)H7y#xOFDz2|ICho+~SZ`}o9wFwRDHUdAC_4M?toUN{`j0b`%faea+_iC2%s`dd$YkN#4352Ty z!`E)|7t9$Km4~i4Qn-c*mF7u`4Iajb&S0=MxjfI+B~rWg!O^DXT*x(TbHl~qK02w< zg6+N+0jf_dzu8r+RE=&81Gf12^Cz&p&p@*$ChR98fMkJ^fpUK)ki6_#N zH7%$1w|_X+Rs7@IxA!Ep5A2ZEWrg4~Ajv?a4&f&lZ9YMGp9)|*Qsn^lW4mO#Isrf< z+5#H#S`*qHG~OUsr3gEKu)DhQPX)27D#wrrK9p_(8P=5!Hn#gZ!k#EBfA6?>bIo3dvHStsQp3U#fQw~;RVVQ_E?@2mUieR)WypGvs)sMb>3V<%XbfP)Ck{wK zBWN}w!DhOJ)|M_^)y*&pOF>PE-zs;%Pu4-XQtg?dzI(xF?dA!mbq{{r@vV@``Jj9f z3YY4v;2|9q#A3P^7tQ83e}29K&IXhs4M1;%>9$X>;F2N;4@;;3AW-B19B+8Q3=9YW z_0-Y5E;ZGsmBBj^$g~BZLsz@;+13sy2&QJAgwAD$w%FbWKHRH2;-m5`oW7#9Ue_yc zUuopcij7CWOGyT;ELH{9U6+1fF;)NAznjarNk6-HfMT#*)>S&J-wECc_ExK?3JD2|G7?--u_&4EomVw2wFF;A)O+H!?y5Mx2D-~D#2z+r{wQ$uak zykD5@7d9|{gEeDcNqdhJZctspqzfGVOaLWEe(fSy9RAfWJ*4rI`U`dobR8C%WJir_ zqWOsL>It7i-G3^VT90-3MRU}h>$*j7pDqc`H?T3HmjGj=>*CoL8gU)5CZMpND z%{J`$skF8wy6U|ns0N0iw!KTdsrK_v&b6-o3rLfeh%~FNV;@3oO!-NkqGnwj*KFx3 z6WH=z>49FQe^bL7=@?vrC~5Qa!kUJ9NYS)H^L+!8*Kzl1BF4(P0q%a*d`|;gyod=? z5DPSa`uf6Wx3k17jlg>fWTQy*&d9b$jK4Fc9+2)u7$M<%p)!^MMs5Ri{nu& zZjFU1L0>2yeTqVR6Y5U3U#rOrDbcX8goNy|{(hLxs*62Verrqpm6L-fPq4u6bd`u` zKzgkq>6^`%6~!p!HNmVH(2R3(@%rVK-3s&(piGoD5S;nC-NkpODI z>2JhqgAxHP?DP4i9_)ORfBG!yhwW2fbhRF3ieUiRZQkr`-U02gpNX&QX|2`-)O1Da zAOq-zXDk-B@0jaKtmK#Of(Fbc7>v60u%q`y1Xf*C>mqbX@bS{_4G;k|`1ZtjLE2|o z&I*KEqhLGxef)DOl`eb^8lDz{o3cXw-iBVGy0PpiMXuOikVZmXiD#dO`VN`$kd8J} z*BBeMq`B~%g!9a#)}?$cM*U@o#VYbgVP56U&LYCuqjnBbQVTueM)^9ttdV}0@Lpa; z?|DPIKl0)vt5m}2`2S0%Ez|3fs30E;s?MJ zk^#?tJK4jYyaLv-$?j(O5IN4a(%t9TQ?>K*e2Y9uU_PfHLq=j@3TWblvnE4*p%Og{ znJFqT7_4r17K1rgT#A53W%RQUT(E5xztXJ+Z+@j)O?%J97iU)8pY`D$LCuY@P8|2O zx$}k_lPB%2TQ7s>c3;Q|R&ChabifKLK+zeO$`>H^1EWsg^00T8HSKcKH8S5BgWuMS zc32PaU*o7V)|P^~ma13%_;`2DyvI3NV;X;+j4?L4=5~49f1aG{A@pFxr4*^_O@L`NJ9eM`28& zz&#pYy~cIwEGj?V>=a!{kIee{bFml&#K#IyuifUq1?;04&H~M}5~0&U%YO^ghBpCk ziH?pA0XJcqQwRVVvwb57?pJOfBH3&ZugyP|txD)KO!+=~h-DnI?j7%Zyg<{r$Yn)) zwf1P$tCLH=Iu<(~aIl5JPQuuoo$RRt3Tz8HMf#lvU;^jHjT`APU}%u!$qx9CqWk@$ zkNCqOjAe4Ff}KSp#og92Jk&$(H$QhaP_17JW3fQh*oMZfj$y|i(iMj;ZRdOK`~qUb z*(5OE8|-wVQ!WBy(%^g9YP2MKc>4kB1e0GOj-RWLL9Q$t_!2k5 z5Hv&XaT!ylj+uM1_Ih!s3;!JhZU?tMAD6=pM`1g`76Tc#*Civ`nQdziSU_1>+VNn# z#t7yhO|Y!rvY&AJkk(u39@zlhsjSs%`Q9kkQRw3}P`J;j{TK<1au4fxRjZZ3Sh!eS zEdf+n2vD9GZ(PN(>w7}QXY?BnJ5d9{>TLGE30duC!sYr02b%y>tp|;caxg#RUxgJE z?mpF#c1+D|O(Ao{lA7FcPqkjXkhhYAtFGX6bPNLZ@sPSvC#|xGhrt2 zdY_LzfBx(-MP!hgk>L-@IwiMm8<+>I2lA+Z?jK695PKB&^a*Pq?f3?Ud#|~ETyLdG z?E1c+UY$cZ&y9h9BS;ag^B$8?c0s|Jiy4|&o1iLY^clx>qErm zzN1Dnc{b&CE*H&X)N#YL_XPqw3Fj5g>AZo>fAxvDSl|<;hZ?V)>>O$a(<`IQt1j1A zIneuC^!CxEX@J4JW=&2D(3h3$4XAMrq{)iZtgPXZsr+BZ|* zffwV&303Exungoo5S(m-*qqL&%YTO;WkO(iXOUf!!QCKe`n=7}Kp!2;P zD((78Y2{YHxSq|sRTpP_FeS8b!4Pk}TjKDm5y(8RU7YNg76;IfZwI>){C7wGpRO0;+D zQ`Mzy|D+Y>Y0A?tG2cJg8N`rsGzw10zIEZ2<%FbssSVy=RB~5#7t1#V;4c1Zq0Z;MiBLd$9=8cGN}bk*sVYuR^<(Mf5EC;fig9zxtWAOBAIJo+!@q0>}g z_hk}WF?XXdcy)cDBAC4I5y2*2n1>^D-^5vAc43l7e4sil2X;dGH2vYB8HqU z8Fan$B}RoJ6N?7PFM!Bmx!&)sXNi>a(T6Hbx%Nu=G9&5LR~+haa^b>HR$dEVg5`2(7lrs>ob6Vo79^#8(kU-P(E z0gk*|0;~FE{=MAJ=#IaIu@2>JWo50LIADRc3odY@Fo)SXn|sQJY9kzUDPP4syP;85 z^VJmfTTpW5#!cM+yq#Y|Ggz!TfAcCkTn(0Q;aV&!!_mW8m;Y*Tiz#-g@!8+OHM_+e z>>v@7-UUu&az?@-pobV@;hr_o_mtm=ylW@i;8usj%Aq;+)5@=PEMM7rZU@w|?ow0@nt(#1le z{_gj04O+Pq7%;XGGMrNLx&uMBAIJT^8Y(ApTybP3TULD3ZK0v=O)++=RqGmFX zkZ|2&fM`N@4`yyq;-v9DRCTtu%?`e15o{DsAMG)140$bm=&HTAYg*_6P8|6sDQF&T zqXk#TpRm>W>?L|iLBiqG-Ol8OGv;=#mw!uL<@o(@xbtK778UVt1?8HO!X}uy_03?v0O*COONsU9(&-f}a=>r$%1Q&Guj%=#qh|nnsO~ zq7V>x+PTj?@-E7fOj=gX@%rW2n_WlsGz;iek=4Hbnz32WucP{F>jp&bDuUxnno%B4Y7i3#WBpYK1LDEJ-UkCOr%pon-Ug1$VsfJMd1y2%?_P-+MB*q@9Ds>@9Uz#6C|npgYfe%q$IC+$qA2CY=RdMjI#!bD_g3x_ zAo|cpD)67s=Mwaho9%87hee!*38_5jY#HlUm1+NOZ+2L{1QXWGmKUpejyRhQ(`?6! z#B`Z2qAt^5P(xebG^j+{+ER73+;UOXUCLm{=cSYp=r za3#9C`UE8asTLMeh0Yzk^pWZSPUEh41LKYmrMlMDnQ6bLG9-1R6equqN;!kn z+haUp?PQGdXHh;bB`&tE&Y(clIHP#GoqVhRgh?uGC67-1#P0v=i;cK5uDe04Tv|m) zb)#r~dFzmK_-1NY4;IO1AN&i+>vt-}Rr$2P+5%x(ZHHh8sPG4_*mzN8uUuOSDF9>s zxCQ&ZVQ4accdlXGVVR|I72Z^jXXAFS?RyXW^gp%ELJCYgP4dl%8x0kH-t-)EL=>1W zRW-iC#-hHyX>n6zR+Zn2^gt<=CS9!^C?42zO zY14|avk`5)M|9J#z9S8K!XW0`-g+Z4P)p)?NytjOG;`R2)pL0umYz)Kb^dqvm3X7P zA6RRj*NX@lN_$E(W966;&TK_yLWNl~2vb){d9`^T&5!yeT7 z1dw@|am&?a(s5X9{aqW{94Z#yW7l{onp~lXfJ}~5t=7X=rt*G*tA9E7rHCiWLRD_) zAUOzhjg_Z`SD-GwEZ2I{Rydc1Qh|51?NU#bG0v2S_k*8erV{#*WF9sTieo|uXp;{H z@bzI_Fkf#) zpCN8h*ax$}x|irU;&I8yAo`tHiPFGilxmmh6pcDW^Pv8Xez~psQ$;rwBl0y-j0YsS z>WGZ_dQi^bD^TdN&ON9-V$@Z*@btj@Qo6v`s&S7u1gfZ7(#|5$#qz5?$KH;%w;Ly$`^1mnZ*rk?8Zp3iF8VFdnYE==xc_&w>T) z(=)O71kIgwzmeQ5a8W&@iX0apG{FXzUz733yb6rLU zo3o8C)q&$TnCVa=P8U=lb1A1XS35WgBPdl*rmV-I@5CXj-ZN~JNsf5__WPGHd3PzJ z>~}RVQ$*0eg|*=UUKQ~@aM1P(i1jgZ$9%{ z_hVQewan+U9hZm(&j}?<)*w6&!LZi)MwoXZpyXTCpn;y~I_`#j=YS{BVl4srcbFf_ z_lTC%5k`}J)ZCs5iwInRT3Is`+TG$3+O36#t`-bbvo$@`pe1K+`pfdoxQLV5592kGkkt(VPVmZeZr!MYqT$XHPmyNK>GK@8SFNwPk4@=iLT9; zSh_en97yEg#}U^UPUjvyW1R)s{njMHz+|YM0pIXs>nweDOj#=rTf4r|L~5b(|Kkag zTdm+pjgdbu$^u|sPD0_)Z#FwyHgx;TqH|-iaY7<`ll_*j&^1uiFnVsB=ul3^ zp}3g}SVElO5gugq3}HD^c$&*O{2WZrO#-b)Ut9YVkOQ$YJJi~m>+9E=u-ibylStkC zFXDi@8?d6^^O{|f)3wFUkzzA}_lXjmXu;70iW{##9ITg7zIYiIxpF`B{`)U6j@$`| zDA9Y&Xv&?cxo{1!_n~tvr-x~4n};T6j701%SmZha*s)h#3K;JNotjVEe4Ry_g~pe` zWYX$RJW%7!qX>Y?!g|Qb@Xxr$(~qVoEl#5C;EzBfGDkkhKRxkUJ=xe>u(G~W2BjEJ zRQ(bdqy9iNjLdsWbMuqmNZwr5%oa>!tJ5*o-sAEE(WVH7jJL-gIi?YrHgz9cy|xsm z3*Vl5l3zm+W1NpHCsK&UVcbe$8;rEh2&5bkuel*$BB%T41JJ@!03*ODC(GvN6Ckm? zMwhU18w4wvi_S)=?}KrX9x%Aqj2+d)K1GpRHXi?Zku0uMI9c%(a%B%xOk7}s|N~tq~T9+CHYhvlcidQx!H}-jzY!44& z$G5izj%scTpd5P;4`oe7n#9H3wxbiE=Cj3+*i(~0olK)9Po-G3HiEMW-Q#aL9wcEj z^kMD&!C`m4pwZo#uDyL5Qq|u)D8kL%mQpyfknpdLbvtV_G12B#KfzN_1=a;z1r&z& z)R&rSTbE7WRULU<5|3kyQ!@N9tdL_WQL197MYt^6ucZ?-Sjgk? zLdG(yTVfI|T*IK1a7oPt2HP|>`93nZ*&4!?+H_XfFAOL`gMa?mgD%B$b@j&MgDp4k zFy}zfBp8U|>&|OA0aI7uSYc_+a}9=lTxK>aEp!_DWEg0E(f(oJKTFanJfCk&ecaJi5O7?Uk?3|O`MH0 z3KKu*21eLs)e(!Y3CdJFV+aJCC$`&}WJ60`D7FMm2;s0f!W_h~!%x^Fl zKQV?BKRbrlKTf-PY^aCm^1!wsb#}~^42H0UhV6Y`5+3?6vM|s$+ct)t^e7U9Aw<`9 zGKRb)^rbU~gIkwe^&de`sOcGErPiV9mIlMY0iHbvNhZmI@C84cvxz|Y-DN0@$h&TX z)k$y-U#a?=`Uc({-@?3csHmHS94A{Ny6-ICuYs2&`WUIEWbnL22eE3rwbW_^?FcG6 z^8;1IvGt|_dJv)8gBnn&?^V}%q-K-v^ zgeNt%{DyNePUVo8Pt+62ue?l+yQ*49T=1URkOu!jdXMN4Wsh5@{8!ba;;NN3M=58- zwwvGfR%J)TyO|~ZaXWy7;Mta=(iC(XG}c97K`PDO$F8xB0RSwvSaYT zJ`zH$Tvw`P=ky!3W0}1-o z8C3bL0mIu6CWPQ1lXjnlc4B~g(@q<%rt+7kOhz@5*_X}Se+&aP~LYC<%Q=bg~kcpiJiL(4EOf+}HmA(|9 znIj1)T$paovnY%baw2~-DfGMe4l$1n{ngaWzBbgtPaN@(zsBWKA)bm)w9@6!_Lp`s z%ao}Zof((lruTc5U5Fl3F_dkKx5h@SaZDblO2n+Z+g^HZ0{X{y+g9OYC9TLL?byItC^mN_23CS`CXMR+;K)5s&)XDgy%TQs;erdDGQ?G?sQx6q2&AeYxe z3!k`uV?22`I?zywaoN}}r@s?rdW>+!|h|*BV(i zri9KLhVXg{5Z+8~P-|5-+>Ic7snj`qHqpv>m<}3yJyoZ5FKLDE<(s5 zBSCYIKI~es>$Nb~YpSv@6r-R;kx>3;OmSCzxXgoedP>1DsAm*qNzT9Ffk*`kO=hbE zya)DH5r3IU> zlU=Xs7=9YWYHniRPNYzGo;fa!^N9IG>JrV^;u4>v5=U%wcO4g)IP!!%RrKLbQ$)m& z10z=Jf|?L>Is0V0^MmQ(tW3X3ggi?i%k`#6h41Zng|ZQvWhYQ~aRmt))?_6ygCS5q z#CG7~; zsY@&d~>@Q-=QbX`((bYcPkdCrDk&rn=-6mg2)2h@fShBovdWZI!DhPRBh; z6_&3@FIORDRKNHLm#QRn{@ru=O6c}@U&#P2UAcXSrVP6GS!~?hg-Tpl&=|=#y%J{J zoa9z+o>UTLcDDASTIz)6O@cMDdb2D5x{4?h*He@HvF{-DVSzmz>OtC?D3Kf&e4lV` z9`P>ynYMyhWRffeIk0-4Y*ddaYoBRUk1T7UDW{2ae8?LsS7}P+r%8-RlJe_n-`iED zgV^{uo(bOMZOsu;6AxbrlVX#Ae;#&eDUF9CEy*O{?~Q1C*6Ro1v?0c{)GqJXQdf2n zj+c}(Lxx8lRL4RSCkMAhEgOR^+&8H8i3b8&(T#$MrL9UNxZC^d?n%i|DO*TNYw+;) zwRP~b(lI@c)PNrX-Hh3|CZBtv&*L@VH-UDB*NhK}+TQQoNguXhT*8kVL+qi3g2iaS z6PpC!|IE?WP{^5}+3yf7TXwIML76$Ih3osO%kO6^tP3GJvS`BZg=`y8@gPmwAZU{o z1$hpX;B`@?Y0-6Y@`h(F#v@5zLX36P7(WGw&^Kx+-F@FGZ%9oD{*|E^koO74yBlQ0 zTu4HMl=n3H1-BlwlSOcx!E!Jz(tVjC^f3<>$J{{Zri zLD-2@w5|cS=qa7~isO)xfyhaQHmRe;;awxREgBqmc@JJz79?^RGjw7Jst+0MT~mp-ho=f;`hB zSalK{sxvG$q$2m-5B$tbzAzdp*WY@m^SwPh{PTn*63>b@#Ts2^?eSpWik6{=&vn8U zxi6oj$dzGjJemq6g$%Za6QYgn&XzjS9(h;YAa+oc!NuqH$ZRX*0g zUe~I$;wrw%$0yh3`n~P#HxQRVLJX0DuT;~`hiEymM93Pq=Q*Rk3nZrtBr{KQ7c$?@ zzDJ>WEvO}^nMt#cG@{RCS(D(kgfID8eW)Do^(1Iy65eOFJH{V^CU__3Hh(9+Xp_)WS?AbOt1xI6QBM4W&Ok-X zP?ccqG>>syG6miFjq&*;=jo-ZxF2ueUA^M=IN%8+Qb}8I`E8nuB9YMA1NGt`AEPVf z4?R?jd^X#6tlSsEn#AwUP~LRx`8r{C9HOO_Q@Wj5`YpNcWSSK#g%kUwg-XZl!td{{ zF3QIc&K_h*eE0ceaqyB}?cqU`hX;--rKO|?2*d8aH}T(<1?3Vfg5bS$mKj__aMj)w z{BfH1&}vCZPMa(|s{|i2a;3|iI7k;bmpEF}8*`j}?+a;(wM2B#$FTHSRwHr6eo~6k zOHC_`d>ry07+3HX3UFnw^aiB_-exqyS8(^6Br-aBL|x09!4pD$1x*47+XuqS#~)7U zA9<<1Te&y0b#2q_%AwLDDYEygM_yzy`RY=WuZZ&7-hRiX(Z@bE$OkUWMQ6TD{BI+o zlStRazbl9t6ONiFhu)GStDVts9? zJ|^UgiBrVm54v|B(UUubT&Zx1uyhhR0^eoZl)tv$d(=C>_|yW|p3yko=~i#})tS2) z3hAvbO@e{cEJ+zMP1jQgFs#!XV!;dF#(3xJGd?rlejU=pcrv{#R5NolRxg=S{jj!p zrIS_q`9m+p6yJxn>Mbd%Y6XGQhJ5;trmS;|Z-dl|OFxfV=MLIfrrE&6gP?PWE?&dU zzK?XCk7%mK-6CSnyci(pu=wVY7<~m586(Pf+yqSIVO!e_XECK`?pJ;8QgItMETp_a zTjpA6k3D1>LL*QgBJ7RsXUwSa3B`P3uecF~FXvQmIzri9_~sT9%p&x`z%9#fqC<}y zOW1EZmGx3sXWTa7sS*?QjIG;|CeZc1-R zBCbc#%5Ctf&>y~w*C_80`DPnd9jQv%B$SYU?B2`HFQ?6sxX zbz8T1iXU!d_Tj0^LfzWiFi;G6$5uZ^K>) z6o_f#O1-`<`aaZ$VW22$R1MM8GGufr%TLX+W_>STpv6+1!-DA0&8P)$&Kmjn$s%I> z4-H%Tz1s!i-mcTNy0@Oq^*-lsWgkSZspv1cJWATR3Pv7he6lZles~V%@kOZ>Ix~I< zsF=tlW`t6IUg3)un-+9x%CZf9F--9zkL#JO<9 z!yd6)Oe8%n>B)!kjXsrsZ&AOTd&z&tmMLnAkDY{< zN|MHg(hM|(qw%-g_PxCt@;=i-V!cOc=wK})H=?5Vw58$-vL;5t>949)m3|#+(abk% zbd%at^t+KBok+p6k?0bCT)G@;Z7zo-g_o@w!j@cIN>`YTR-`S5Mh?7;5mm>6wq-SRr!{aG79_Un?_`(nTKw4^3wEqYye}5|Ldzzs zzXoDj>D9i2^J5Y3QQN-C-5sTv`K*WQPjDX;RDBBmhowS}erm{b4U~LHfUtP{=Z%As z*JA^BOcgCzo%58yFWJ^Ba`Y1eR2iiMH}0PP&E`6ImIVLGze;J&!!`BZFrhG7#tdh^ z6K`kz;7ojyJTROwb?@QB#J*CFB;s>Bjk7uH!B?xlP`_MsBTZl5ImhNlPjIF*#a`2&>*bBm#$&U*pxK7jGY_%CV)ScwQK6eXDle8E5pK z&ofuQ6)MC@!lKc&m|Nie$8J^%c{bVvp$KgI-kMzBNRrfg7%PQB zqhz3tKgL~kKwenqW|utAx((x0pLh4C&*k5iz`dkss0?-O1?P8L zRLD5eo-91c$TqCudOONXD)Ef?+eFHF;$TssWp%MdU2=5K}-^#T^_* z<$B@zuD7n;U)X({M4uwU`P^MUtf=Ps2k~W^4wARu$ocH5USqO);T(6KH6O6HOajZo{Lf4BQP+CF|-mQ`_al*H@v7)q3A zz|qZQUu@C}HqWbn)`|m;2Rx277x&Tx4sv&63YVEgreIQd!zCkIuMyNWN@zUG)!{9d zris6~6z7r;_wQu?@EZ(HKka@-Alz}idG|O^>i#-&0h6@C3m-e*d{(iwVe*NNFyaM6o2O8orxhOA7M>&g?vH>+4@m^qif(~z zYJ!f&!zMcIckHU=>M98Bu7juz%874%EwIgKWRJ1UitbWWm#Ko%y;OPjT3_Dvypk^C}xK@g&W~s4Ko&cI6Y8k zp1jsH#TBv3%KgapPQ)&Yxm9#PhqNSn+!_N9DyUa?V{iTBiujTM@qDMG!Ieb6w{hVMP!a}LD}v8p zNx&Zep6Nviqk=iq?>#i3J47JSdb3cV@8DYEM)pRChb2I1S!lzx=jcmDblkf;e zoe+_LL|{JZ*Cg* zn9O~%XT?IHS`!X8q!P)LCoIFt4e>B1-%C9rw?nhZcALI|wOKV(v`5sF#F;ao21>!&|TSL(ZSEdSEo%YoqZBO&&+T~8+u!vfPnHY#ITosawSv{n`ye>cf`u6=&m z{A=-{BR{GQ(e%=Qm;2q?^_GLplQ2ISuI-Ba6u1M2+rc!Irpg(|h38i2x|_v4rUZ_| z4c}{PStUE&3w%GlkVFz8oSy_U9Ak9&qK`~>1&tF2pe@;X~; zs?-G@ON&Cx?u3jcHa*qZd1SMnRq!ZD)IUOO?`cS!>*P6(Lx*3Wyk5iML3FCknOO6k zDXs}Wsa@^4>-sAXjbWZNsPXytX%jzQ@~H-T(y07s+N-%m>gV-FMBm|wzZ*0mWQl^7 z{i#q~+Fzj)D*b0%=V9dD|NO<5ouk6=4N^{?hN`eEH;-6XuT|h1^afk}MhAVSp}U(x zf?-BdT7=64n{+bVE|>lr4pO|@Y$mqfhf6+*)fPV4)(n?ysWPuL`JLcMkGE~$YFjVS zYpq6TURH+NR7=P1IOW8aFv9;wy>>}D;ah25Zh?iV#0M>dftRx8=NuS2uhg8Nav(ZTZQ1#J3eXrW^1Hb7Fg5 zx5tkdzK%br)Tl8OvSrP2v!>=*mnm(e#*?Z=UL#9#<{&le+T4#m`50;FIIy(mHc##X z53^Cj6{y{D<4CZ7kR>lhRM%S9v;XXwR#`I}Z-E`ASAKfn+IOtf((MEaSe(kN_?*Us z8(;VEhff)!@>K&byv1Y~R+;2)A{U=s{G(Y4b-6Hri$XDa+$OuN=c z%1S0<$oy$#TFC!vJ>2{j%03jTaxyhI$Gz?LZLsR3FyFdhlAdoHeVpRR>FKNx!dc%eyg6nVcLe`c{O`*jQN=qKP54 zr0@>dyre{0=Su#T>%WnbXYx2@(~IMth|Ue<7~{9yV3SqyI#og?vo)ojDAsc4Y^ln< zvVr~u68G)n&Zk4a2XD*%d3-Ooz_Z+!S2#|_<%FB}k3i+GAd-95=4MDdAy*qy9X-E9 zT3jY;B(r~y+%t2tM}*xccLR16T)m3vQ!4Hu4TR5|VW`PgN8#iq1!t9qH}CGGuF%i< zQIgw9`Yyeu&jtaxNXgNIO(;_Aw@%JLYv?1vKtglffB0FjFdBnyEChGa2GQ`OD_f3% zcXV0Igz?`@hS8S6xofzjEycw%@9C@)8+WRd-?^Bv9HzB9YKg!rDzQ$ zl2csO{a~G$&n=2a$L%~7FO*a@%!4>_*0YUyJ+e{bh-Uq|~gVCDvAKYE!FKnb5I(XAP~6Ri(@cmjbkW zqN0KMm?k*Lstzgc@j_c#b+oLW)9@dbR%0Xcp*ng!%lfqw&Sy5&4Nj{_Ua}1D*qyy8 zzq$R`Y=e@?C{o7afDCil;GT6dJLs0vn%@{Sr+btj<&%yMPeYClx@ZT1r$6>5BC2aG zX0kA`?fZwF@YvmAuNW{0i4OY?enj4 zQ-o|}4h>4=TfUFzmy7q+0z*x8YAQ}%ZNiZHN%>xJA2*y3-P?Dh{2# zh>}vg;VFFE^Ls9WDn-oJ{zMz|?iq^@SL3-`VN&wqLk&TypIDpz#y*MRjcgu})#S zhlW{EBe0?7-yFODJ(;$wq-(9-#1HHirQ-I;-O4%b(OosW z!?h<9XEwG6qOqpj+BYVPY-ZMZxpnYGUTee-I@5)Y4ZZ8L7EWt!j2NTs+lzan%kyO{ z>J`Z3v(kfqVYjVm+5EOvn3vDdEHls%QZG;zz@kJOCzng_Ds}uA7Sn!WRLN$N1dS$*7Xw^v+mPWVC|slXdVi z^V34``9=T1s8Pn|6Wx2;l^Kqepz9Kz?%!;zwG($;PA#UKIv38Z#4i(@7=j_TZpAwg zToBZTCacl69mDwqgIri%m2oo{C1xiIZ?Ot<+o%ay5U3zIO;&4+W4SC{Mk<^sosX&b zhT!I!MLN%uwOTf?hV!B|8H%GdUl$b8iFQ>POq?wqi^2Pzb!&-w6t6^Q7#W4#&UJ-D z?laD`E_($1LGgC;!6wb?TIF&42W!XXr(T=ULZeo^P+vW;VqghtD4{HdZu89Y-ENdw z_fh^a&fzHww;~0&j|*j!uE?KZK3_Iy5mm zUF2B}$XFV!T5>MS*KZ!?Eg-e9=oT;#@nXOoA(7A1^-y#tct4N0n;iq@9$=52eI6KTV9QwaAKS|Pm`6?-S zyeodR$m5>Qv}W_BJgOJ=;Tp(NT^MtW6+&rdb{j}5E!Jle-KD&_f)`>Vy0rc^Wy*m{ z?0DTd`iLRr1n#j`IZ`TL06(tsxaYE#IGTqv9I4gm=}h`Ud18|$-vd49cPuhjY=%*R zT30o}rLq*oc$tc?Q4M zw(7$4Z0FUuHYO)#V%j4j&o=ibNZ+AEo~7Ob>3dI|lLn}?-s<-u;o(tT;?rYIdYsnC z6(VTA@d)8dF|ntt(LHx-Uo4i@HI*T&IiA7n2JsUNVJ)#U^qRqz-oA5MQVGq>jvT2B zS*WUZ-ykN2)7d<6OwSOXHGhtD(*8Gn%9eW~ELjb|_Bc^yqI0)W+y%QFZdwep)d;-6 zjUYTa+>Q{0ZVEJd*>dtL43(A!Sli~Ct*wrX9=ALmkAKEpO|Jgcc=gmcR>WK^Qi|r4 zjeSsUmgMV)<*k`RJ$7$gN8i0h4JL6Cc&y;Yxg%;Np02oChngSM;7zo}_3u)`_uZVu zbu`YDMO_h=ZR=y>@DSMkT$cd8Va6&b>fJokM}&g3(kQP3$R}nZZt%|H(gxx7_a-f< z3S+4bUm+6FwLW&t0scFntWm94xsyE*D|Dfm*`@QZg0gta?Swe+Sb1Qw&Ahcusl#@Zzk4%~^e(94E~o zEZr4s*wcn;5YaUD`bm?Dpb6dPIeXOyoJG8`ht~<+N(~>?P24hT4y@YiNUTYT zCJagZ*BX-iH@6+<;;gaEdbP~jXpOp{aR1xLXr&je`ZGt5;0=IPnN_AshZ#VO9cb77c}Ei|)ooV==@E`f zX#DviN4CGn-Nl&JwXPaeZtErlcEghSq~ zo1W``;&a!tI4^#5b;UmdflhozUFz?W54~bNvzcPbJ<&;rcThaF5=r8rIU#jl2C@>o z4BL*ZsE0+bg?dttD|3GI3*b5tp$a*?Gs`ma*#pHWzs}Uz?|kmT^9O; z?xSY9|7!V^$=a&MSm~_ws;jE`f|?y}qu>ukQNJCFfq@2S-{FmzrI{g+J4hS?cA+sD zxZ$*ArhS^l@uFa#3N-n^YD~^tXHCqMs44yZ=|`t}C15K|f^6W{4)Y42Nn5$W!%PcB%jG&k3yXns|JyapSyO;upIvyP|~Wqe;Sam=y(9w}xjANt8}hk$^sRYg-U)v#r_WCt&sF<8$K_1Z0!xhn)ae6ED)%IwPxDT#q%8%7GV&Y(_DVXO0`roAg zhtkthRhZ6?6ZMUmO^veRjH%f0Jb~7pZ|3p~KiKPh_)+!*WO9DRURc};^jmP_Hj3$^ z5m&dOV;Q-oh4JPQhjS$)NpmZO%G58Ye}Tx+J+SFTqvgg_?lq&{-e`(b_hLin@Wk}b z{5^B;x$PS0oZ{h*F>Q^?1fGC*3pq!^KdI38zk5ZqD-HG!pRmcObZt8I)|z*)d?>eK z8K0U$@~+Il-_Zojx2R}7P5NoCaNppqdg=_=2fhF3IsIiV)md>a6Q~3#+04tzOZ6Rw z?*r!iz(U#Qk>Ne+kJOp5HrZ_L@AIfJT*O}&rT~hsUOQ1Ada=#7=28yX$|hymJ=3bb zA#(gJe|trZyq@ynM&%+-m2Kl2J4En^{YA?O*5MMgw$0qk+e})lNpq_)pSu+vTwgz5JsQ+UZv>?0qCZx-C;1b64}X#pDn259OkDW!m4%HDDjCLOuBvut zE9AVc7;!=GqMy6*W&;I)lKbv?E8AfP>Pa^CJ#6fZfCw=euDtU_wvC-beZ>9tcd>_8 z4RsKZWw}$rzmZXz>RYBdvyn^PkL9joEK$(*tIE;e4jt{nO5v(Ky(8N}@;LtZJ*hm? z(U>v7^^dVbXw;HAybS~*v}`Zm-0B<_<|Kw+E?joZ#_o(CQ@n*Ll-YDG%&*i3A6#ocr5~ zsJC)cuR`1|avi)pW_8~<;z;b>cuia`^>u*IEV_7;Pb8ym262z$AyYy6kI6aZoi}fk z*ZQ6AC+q0?^@tyMM=>yn{@0z)hDW1brRPNg!X<5&ywjd8KYw3YIp|iRZjlb0!>lP# zu@(I0ANwPj((uT)Z}nq04?HTa&EjmYw+PkR_fW4}e49hZcp^KsJY7DJE{Sva$MGtp zX3Neb#N&xDPh8X<0iI*20ru8+O8VJp@-tyN(?z@)Or%^zUeHv;9wo2Wvd?M|YI~n@ zkMY~qA$oTI!tub7ZPuTEtEBL=EuuW|xQqg;psg<1nI~Dy=!^cB-VVE{ajJ2uRc?=w zy{Pnj{)WjKr}`unO}UB`A3HCjrK*&$xBBP5P^OFwbN@LdlXNg@`$C8@)P$0fTNc)t zX_aiFo0L9M+BydlbBM&&HRke)oJh9v3Z|iPwK(E`O#uOMnJbWW>(K2D=*ZT>ME9UB zDp9GWheIH6r?k9jK-^dInh2J_;)GrQ2tE(SNq!g*UgcHh(7 zli*;?q8d#Gv%TQtu%e#wfY2Qgm!RxJ*iWA%X6ZM1i^1cC4Mqv-bK6C{_P236DKc}{ zjeV^rB(-%C$bAqV^CcU3SPb?m8$}NY?+o)Jjs2#%-e%K939x?OZ%tq}>W-n?YnH`YjgEO4RuUp}@fd zJ>vMGhB%ZLtbJri73RpZPT$SwRkq%mD68BqrLHv28`!@H7mT(W^XoXcKv!5)b32a< zhYFo!aD$;SdMl=B1m%^u&)a?1J6Y7?@wA}v@aKP>V#5Hx@T7nC+7^+*%oMVtC5p@n zI>Sk=S8Mp$(1xR7QsiZ23U?N?R(tMu+DTIGDy6+a>d`z3a#TM zdES(yCB9Gh*OUUh2A6l=iGOBqyqA*tPzN8`EtF25mi&M-d)^?kv-S$41zu~(na2{+oj(*VV9lU{A@_jMS zzaRdYtIkS{KiOrsJ|guSmJC~)GpBwlSm@!*P!rZQG^>%Ly}w6(jJEa9%JcgH3;CnxmyHC?X-)%qq0O>Yz6jn@(ymdW5dC_4D-=88=P+A-d)Rgz@*m!of=ERo&kS(9`EtU(u~x za@ZIvJ5*`gI5d)!XNWy<@}~-5DXS zoJc}InZZ89#;Wo438V=6di%f3tF1rm!2lY{u*9Ln8Wtk7CMi91p4|;;0XS=klzPd1-yLBTA ztK8=H+d(Q%`NS8Uptj?9?(;KjY!6fiQSd1X&rh7<+or1w&T^fJ|MIe?+tX`gn)YM@ z4KRUCa&1CnwZ~s$mhZ5nLL?^NH-Q_R!%{?zmO z8&hXko(&Lv?-U)b=?%|%LFE34itLDsH*$~Rxyk!0_u(J_>;m}h&u+6)hDLmD9hhjN zn+TF*D5aIf9n_m>b7LsobQeh?C{k&Cw^{hUiWiy^3&W2XLAFvmX+QJy-DIa|P_=wz zMi^Ki)e#8oB)@sVvuwMX4f29Dx*^yeF1s8SHE z8hD-tufLyoq4n#lj+h-G_ZPF9`%)AxAk@{Zno}MgEO@X$w4H;#_(Ix}_~(Mqptg(B zag(Sa`uvsotdrOAGM4DrUN1;`dSC8caNDvLUr>Zs0?cb(;lj2;JTt-Zg|qmwnUe$TP$1UL9QAH^465Ta`d z>jZDi(b~0m7ut`r9*rkg-)VrEZQsp5pXqyZjg>KqdaqSjSmy7~zrtzR&L0ldKM2E~ zkaaB^EgeY9qE#KCznZzBIi~Zt9^N+xI&ZF4;~(gAUZ&+chhz-vwbsOJZmYV7w!&NB zQ~oP5kHm+54FDjNw?%9Jd zx!n}u*AlFl(2mA-o!2{J`}fGdpyrw&^9fwhu_RVjL${sU*_K9=aQAw7nlpa1h}R9U zSE)?LfgC28Y~O3d>ZG;RvO{XUJ;Od2j5wPAda=~n!GBBMsUuT7=iRcn&%cL;S~B5< z0^6gt5?u*N4)$5)x~~;0-84dl%K7zeTaEQZx2mXa9c(o4$C+MRrz1(Sc~P^neNc;rAerxy4lR+1T0*}O8A}$+TmOA z^hTV241sBGo99HFhw<9BfUTp-#eG+Jn{!A!9g^O8V1DJmW%Sd?$OgL|y1Z_adKNat zaMEw0vRX;@5~d&K4XElk+RkUamC-fjKE*5cO4EA&<#Ydn^K z>1#@UPZN$s_l5_#rJrf8XZf=B8DH1`$>!sVb8I6@5`pj_(nfgDyF{){HT{vY{cG}b zL_gGK4CW!sQoSi2M{N%3`kZjrvJJ&)XX3N1LrMRc)1aMj54&Z3Z5E+4JpyJ`=+%R{ z`|~j>Q2DTpDH+p~&KL9p10C_eDV#tD zMmC6&N$;X}1N+t&{jIj?tmoU>$Kgw%}Gyy3tqr z@iS2o_xofzEW{W8vS#&@$#3lsgv_mF(FRzoPFny#3Y^H+397|H+I9{#-z=VXTv@D$ z^awne{ObHuheR#wY2I{FHUgTO8N#drSISF|tm>w+Z(<;ke8`zxaJrdqrIBU;*!3T< z5w}fNPXiXnBKO>kP0UvSx?J*<7%XhO6X<$U;|NQ!eJzi4Hu8UdYmK@CluHNO-W^@XkxCAfE2Mso@r zK&8Md?giNsVc`k2TerVIa9Ld8K5C6aSmj$HU3#PrUpA+`ttLCS7xdf<L_Y}YS=%j`Z>TOarm#a)7T-$%sY78#ac**? zrZ=hFf8b*H&=mJRdg?kaJ%QMGdayamdf&(-AHKkus&Rdew$q6p8O)>n{ed``cVFQQ zf6}AYTvTMA9>Sz>L6WcCcdfoPG`Vv!C_v(^1SmV4z2Zm2=-I7Hh9v4rOFd~LG1#tu zG;OZ{_H;^y*U=vQDS8usK%PHM`)}e$<+!ltD!-XYw`QUrPzIzaAknm}#i?bHjekOI z6U7`qTRsMCRFO*{!cFc0p2NOMwHUlp<1QQZkPDu#@g#}zTj&B5fAH4#0I54Wf$u}h zi!|)^n>$M1OHYMln8{8_jZwg(zB}rjL=42}zcR zj#n5PUsr&S3vT~Y>jgS2mdGfiOk5otx4tzbz7q^n30>L$x7V10IYe9-2wBMc>1224 zLv23pdf`mnxAG#;KY=BWRQhg!c>w%kU`~~P3#_mgb-1gSGixC3Sk`k)>Fvhdcud9m zbh149(NmzP@KyAx)&F>i@`cL@ntbZ&3WLeDzs#OZyDEqspn5CvM%LP!7F9p8{jmxz zr<$T_p`xz}al?;~?|D!jA^6x*2hLEk#S7hRGmY&-IG7CLJSS2cs)6(Y3>+74Apf+r{>Z)UNrH#K?0h zPovQtPPct;D7o$rq={SRe(cAT>KlRakGK_Kt5U$68MkY@%f`*^6he5@1Q`-IY z`QO|LAMNWrcW6JrNDWs=5+f9KaiG%(Bo$~=ymML(F2Z^nRSJcQV3VAS|6z~WwPi9( zXhzdJEy-St+DF6`{IvsQD^`DJclT3j3&sRF(8h{7cwum)4pn0Q=DF@K zFD(5k<}}ZL)N{Xnrup=76T>s%fui3*9Un>S51I(~qT^SZN0#O%3Gu&Vm$rgi6aU(6 z-hzPcEPRq<+#4Ss|A>QQ%C=JcY()_u;FN)S*pP;Syj1IQyaLd-HMzee_hcp_IDzTC z_~C^1+)>2L-daKgo7eSyYVP*XGW+Gt;E!N|$gmEIkUp0R>rE6IowFJ(QLg;8b7Z-@ zwLO0|C~JW%U_77oTpdswJnE0G9ki`f^qeC-1weoWKy+IIoZBTpExkA#%n}UK<;k%(H$kEKrvM4!-h!&t->~qCecj;#gJk=z?nGfhjF8P(SRRZ1R4+Nqlu)5Ka2$%Y zX{`7LHqshDol|2UfvsD&**frOWP9c4TKUE|CKnwsJ~u69H!OaaI+M(S`As58p#)=!fFJT#RLkSg2PV!3tw*7a+| zhf|5+-#ku z+_~@6y*3}lC7EtN7>VvE@w4b|K-&6ces~yjGwULzYWSfz)Ds6Xt(%FN!VYefWlxud zT{jm7RTNEk7TPX(!+>)h7e|HBKWXob;|&HerV~Z?NvuChd(HXOeO}#MPsm?w;=lp*nE$D6qud9Cs~RSBjMO?EUIz{ zz#^9O*rxIh$3=P&e>}l{AbKuB`N@p+z9N~z>Tk_=7GWj9c?uZm#A42rxj z4w7U#WtqQLpMT$fn?O~U_yNX2Ove4d6`POGa&5${;OE)2>Ug8MHUtY6qp5xd-} zd2zm{{=>Oi=O2m3(E%lQu_gvkvVO(tYmCioa7K2waLxFRxE8hQtR#$9GwYe9he{ae z)R9rg->4*gcq38(Uu9qMuH0;3>oisZ#~E|wn>`Q}efp0!S`T6`PmxHTjoFPC$P{5O z>er_Ndrs(sej6EmmHp{kY^kN(4eQ%Sp6~j8JQ>93PncX$+TM7s?Zw^=gw+gd&(OSI z|M0@KHMLIb!tAwI^mm5D`}w(0KY^@sPF|r()2n%MrJnaX=qHW%v^HOawM&3wiVop% z?wppJJSP<7%@tdlyF4DnpC}t(Uxu67O`B=hH58I+3A_rL1q6$&E57Vq_vucD-y06* zS6zf@xidd2D2psk4LTCdu)NznIteogkFtbphn^=GZ81lYK8-SLyo0GMy1_6v%kB!9WGf|;ig%DDG|quAc}`=aMH zPQ-zId66P&s1MZ%+W`$$ssQO_^7m?(waY!c=^Rf8qt;YEP7o$O>#DkaCnXC|{Py3- zn%fl^zBF=+u>^9t&z@cJ6HQs}28>!q-zpak`zo}Xnxz`4j`bFyEw5s(67^xU8KxBW zy@9shxRz$Y@!Iw*|`0q*< zFO>k-REN8}GH+;%?h~LAQsyChBO4GLtz0)B6AtrX6{gK4XwzD zit3)3>N(uqseyPEc+mpmV_iLZI5&A!d>hQNrIPX$_jNbylE2!c!$sDIZ1|&Pi8e1# zZ2l*FI*ypp@gKDFZrodncTjrcURY*mT5u>{QS?I_oen;5(0M2?4z>eLlH`V5tM$4+ zGJF#UJe&VvDO<>0*5;=AAc^dhNQ9)l4gc=!Qdttcio;E@P|jVGaxy1VyZ?j7jMKdP z^BdQ<1z~&-%%~z3r9P(dw$WOtQlD6qQ@xZd0a=A$Vx)x$yFhb>pe(MglW1Y+bSihB zV1gP>Rn;(x`@ANP>`o%n(=pR%u2g9ZoAdL36i-!(gWEgq#|eZm-EM(5e#}hk;Drf; zYOW2(E-xg)Q4n}ZO4e#+yz(74$^tSBQtp3uGor|r+k}^=Py>Y-`{gxsKC2Py!y8dt zV&Oo<+B!I+`?-_THf&4Q+LOqJ01SusoZX$p2)X+#mrigq!z>ml)IBRuxd z$FPw`ugHOKzH#0x>{^9!IQa->WJ}sRH_@b^dCzG)Ozbq)D#1nnb3@rduMnBqmQ%!S zk}XGI#q^X5v;S}`efDRl1NPug&!MAq=!=4~j@Vza8I10reH?jZ_#d z`Ti@hO{L!qoIb5QqaXIg$c3D6+rf%8yd+#gKWme668-)6G%?TRMQ#6kc>PcrQ?Ea4 zO1hlv%Ti=?GXJQ5icWP`ccF+aVBzZFr9_s#=iC*x3Igosm$2#o0XF$k|LCxw^;=S) zawNsFW;hb!#4y=ivhn$ik2R@I4pld^kFAS4wxM!w9I+*|=|NZby}rV@@})lf%U3<2 zyIclWq~uqcd~5FDY5?ivi2XLQx2bdAn0W7gro~MQ7^oY3O?}n$(T|O{CO~NTo=}ns z{b=5WpMI4LINYXhhD9og4449Cc9d24u?_Q99` zWLuE-w%)a57jAeH#^pJvoPIfY^{@J|c{015k9o4P7fsY_v<%*OLtA7~)=IbYQy8+| zO=D+!Y;3C<8h!Y}3yko;ue6ZkVY%_srr2(jPQDG5a7C@vULTeEtv{i1W8<=C{s+b- zJAh)AjqIXB!e|mdEg|~qqzGy_K!6WLLLhdP@LK?^KS5hHCltQn?UA5%oJfpONwnV+rti`TU#D@rNm%B8~a|} z)#jy94Bk#q?1Op*Mf2f-x8Fd+Tvhv&KEh_f*R*?cc=0h|WAfG5b5N+a{L?DjE!f+k zyBsD+4CC5xJ#}x@U$(+i#;UsW7cdj|{zDb>KZtCDNqKPW<-wyIb$Um5a~X}8@ivaY z4hPn{_o;hM`A3B-4}JU}+%z7XaI#HL>n^JhW`-j>K*>)VI2`cbS^pS0vjMloIRo&> z!K~r=lC7Lezw$q!0Gs-1F>N6Eora6xku!_##>2CvIR?8FU=`p1s|1{iVQ8uj%{+}?{-%|gCJ!EtG;v@0*Xp=i2^1sKZybXOo~#d+;idqy}D zBdE+*PBR=RBTj8!l_Zf!F=D~ve25|SAFnS%Af5JsM!zxz(Av*YmJxD}jUy~4EUuOb zGCj#U`_ji;RgFDlt`G=IRKJhDfdOX|*IX^oy5(MkzJC3>1PG%=Mn*mo7dKcNDHs8Y zv?{7kArMCid6^`kE{A}{@Rq;~-ym-dZT?DwExj=D1w}V~OI%q0F0$WQXOc~Arb|7Q zpLNz93u5>4M?s4C!Km@-CP1YDqP%tci0yQ(9#Go-t?#-0idNJ^RZmav7-)Sb{``Oi ziJ|`R1x}?_hMkmc>%h4m^_!*cb-2?~d+f)hn;)@{N6E)(#a0N8HKli)rNhT93h^NI zKa}Zu2?+>*_Qk!EIG~T2Ua)`UDP%X}ljON01KZCr7;2AXbKaVH0SBObA+J+6f*d@^ zeYjQw+6>;CVVQ$azu3`nUTUX|eWn&+kXAHV6l5^yEfbzy_@ux?P4YQD#E{?-$uf{m zPy%}nwR;(yD?pcPWOQ@^rm2`Hl+o9xl&4qCcaU3Hs6#50~gISWcI#DUyTR3m@ohmuNLJpRpN z!RfRzF1R<3Xf2)f=9}C}*zi55*k-w#m&Tg&Seca#3<|l&vp%)Pt;A;-d++d?_avzq z8>hv^#~b$mA|XrpL>>?~G62e@QLXbbB|l<-Lh~AsqgCr~0rx^^u?<}C>rx$Uh&2-& zBW1-4-YFIy`=5DSyqVh=jw9odKbKhZtFSXMO=jN&_u`oKIR`pm?*S!VrPgeFu9;NK z_QQw%(a{(n1-u8O#Ww?!HPqA|-6hDbs}r}mcLs#bfp|=uQ*ENKYs=ZwdK5Y*Ki_86 z`!)o!^yl}>0^L9-uNsk^3D_6z<;^GM`iBMNdq&frf_j7t-XyahtuJtX9`3sn#~DOs z32x6k>!Yu)pTD!SqvL@BTWE_8m8wcgPiK2H@^!Nl5@_JtE!5g|UT~)D(o(f(4!wRL z;=;R9K^MxJ&9(!9q*t#$iT+IhuUZH`SD-zWBp`aT8N_Lv*Il6Q?6tp6PNw(b!mt7y zjbhBoY#RNSR@g-jC{+$YQSdA<4widUCB|!E)fnpIuXLsz+3MnA6@Z-=vY-3Nl>O_M zG|(_qcXTWUGK1bLK!g=Y{YGFU0TDpYGW#(t#78hSlS#vCk?n`dU3ac#KJ*W+E=?lq z**BHNx{rI1B$(T>fZMu@(9w!{YSsbeUobL;dV5ngCdyl_U(wUjs##j*1O^7a>j7fI zk1BOJ@(k+q8w}q9ao2$XMe)SnMPR;hfB53)tPER4tLMBMB)`!bSw5!GYgoV|^iU;G zsOw1k?RsQ~>8J(?1Y*HqiUk9*@BABj7>_Z)GoR!R*WZ|~GdMw`lDmM(Vs=3?8K*h> z+kkEbcVbcIpSvVT^+!D6q}3L5;Ax*M!g{eF_3M8)46Av*=`TOJop_<*x{xbO{uc|v z#fvQjfe;(%)gLJBwd`ekq-IGvo!e-n?-KmqKjY4zKglwZ^qb(G&D&)#VZOsJ0j6mE zgIxB7JVw*+_%QN}tqK_KtDcl{JVtLW!(d?uM1?U02QqWZebXv(Veyxa=`%e+WjECw z*D;N^PPk|A7arixT(sAmAP|~#V05EOo&*h2cf%3Ez1}dMUPhVY_a(aFB`KqpitNC# z-oE_c-w&q4$Wm;Nc+r48gC9)astFfWYTF4#c1+C2R0zTJzlKE%f#4T*yQ3v-70UC- zJ)EX|YqVvoY=J`&!gy$b&-XrZjw1*~g-~j)Cd|o%a}Y>{;H4Sf;qDEGIY(5>vn{KW zIQi3II-C6Jdg{NY=}RgxS0kMLhii{{7Ml>@^?m#`UmU`H=l=Q6idj^M?*~W%3ta1% zTg8SzA}-s1KVx@tXQ_gRsi0+-A!p=~!7Gwa#^n{sbf5R07q9>m*r{B#jO(i=1>M)F zl%yTUWnT^p5vGC2yr@)c8QnX;&`!5NC)!b13cDw>juVOwGE@Wm!D%{?Z{|srH1TH+l_JGb$(zd&43v& zS0DOE(5xXk^c_EYHNa;7Ij3&=_Mhy8CldE?lLkA9AoVm-n9=E*>o|>g{RNsTg}kxw z$|HsIVLO;s>1H;-vk6@GxCtWReA{h-Wq%NN?;nAOt&4LBvF(Q?Tu43EWy8Zto8ioE zdYn!J31>lR-q)9T9{+S+=sHgTA!3EO%02)+@FCtJW7aI_0_MDPOMG2G9Kxr@qSl5yc1GTp)z>+ z=Mrg&RC}3zKZgHvg?4TD=Cj9)B+4azo9JQb@FgU4p{p$pAqd9TtxKbrCQ8(C(N!wJ zhI>+UM0*_>*bXkIkJ)8&e%$Nr36Zo>)kx8yqD_L|NMl>3H8XYcDk}$_bM+AI+70Ku z$&Bsovy>m+ zmk|=yrQsI`1B>qRS=4BW{zc<35|Hyt(shFuy5*$bhti`a`&}G}&&mUge|fvRi`X%# zHi_$hu5Z0CYn7LxGaogWQ_ShmohS?Db`0fJqbOcm#aB8PALs7B3=I)huaVZkmk3WE zpQ{;*1%BfbW^BzE4)g7Hk8K8=(~I^ZD~nZQ;zflWvA{hpZQH5+LYyA1r3bTC7DzY^ z`f1+FX>PFQ3k2AAkif$>T$(4wnYhNe;_b7d^z9U$-U#q9w&e}8S0?%#pK&0DEmxnp z_e0pzSlM20!KkPwH8e!}oI_@0|2fi60GQ_MrB`Y0PC=8%GVg!?)F$?63dBW^DtDbZ zW`g3$!Az#V^hr7Fm1MnL;aXq@eADPpA)fwxIS}mFF7KhK_FW?o%(h0VggpJ_Lb~PP z!6Et;U;PeD{Q2Ka7ELNuwPqrDg5kYs^8vC{gQR*d%rTyu`*Qls<`PZ_@yNn^3)Xza z$oL{Ro~hnGOM?KDid~x2HIMxG^q+q!n=)5e%?)4<@izVYGT_~%78Q)~ot<|teM50s zd^wScY-I*$i}ox2u2IeCpQ61v!J|`P@aIWuJ@;PAD>tEN`DnnZn3S*wvmDN?oZTl8 z*?RuF9S>sI9c*onqK7Y*zVEXz~xr5gsKBPmhHF~R$KV`&uWbdl(sQ2dCfB7&p zOCKB$d|9VlJ9gE?hB^=tT|WWwC9R-?BnX`&jO1Vn?3|pOn{Wuk@Hm%!Gx4mczn=;H zvMwLr!@Q)KiXS*RJ4-8J|8u~&T?UeiAgUOUH-IkpKIlpoHvsyOL*RIamxt$v9|4$U zc$OWNwr;A7xP*)(5v9!ts;qSVL>L&;SP-9BnK6Q6Z%iCLoTrloBDO~c=rfeJK`9_ z{v-3ytL6pqK&MelX+xN478#vP$YCBzkp09Y&U*7VcSz)UOqhC5!r3g z7S2psYX*Yx4>mR-FckXa#tOvetGY@;PBFqOk>uVgh(3Af!cS?HWnk}cq(Da&lS%Jb zc-87mCFDwW;DW+KW}4oPjt({b|EIWXk7u%PtI5)Rz!s` z9q1{NQ)Uhuq7;d>2=8Q+6l0Wa4x^qP9i*Z#%Au!d%wb9_r|7wU>+^a4dH;St@Ac>I z`**s2zx#U~?rZngbGxjKonH{Asxk0B^7A2k`-Y^xEsQBX)(1$Mu043#Gg8+x@-oa1 zgUKq~6k%EoFqHTk$;B!$OV0tDqpI9=v;^DlRI!7Vf=nk?Uy}Z~X0n$(HPMqIuLUZ> zcn&o*v=BlMBl|2S_J&ooCnBz<4z#Qh2n79UE}n0ITbc_bx2NGl*#)M^ujE#}#gVRF z?kw_@CS==u7A>7T>7QV#F0t)HFc!#liXwHMQ&tTb@1k+<-5ZJLTuDwY6bLqkl0_Y| zR;wSA4zezhOI%zWe0R5bN?A4LS>ZZ*{naKG_eSF zbOVCEM3%5_k@(5`-9(MuG`}#%s}(7TrDeU zt2gJIjRUlMu9&2aE3VTy$Nsvye<0kj27##erFN;-c*E?m%=}D#Q1QCq&JXGs;_;uC z3JprsZ009lpi=Rv=;Ekf5VTv5EQR^`EgiQL=cMF)I&)S}gFP^;clEJC8+K^yIM*%($jq{qS@(u_!x+*Ac`J3+Br>mqjN z78V~6<08UWEOJ1JDh2z#YK`-xlcR189MB0KVY9ZkP*|Z@Tmo`Jnbt+cQqlscVv0VB z`7avO8lTN|T&VKw*A*p~QRs$U_;dUZuN-qN-JQDoUhiv8Dt=rR3xjkk+AJNc4%-#0 z)Md3b|D1F!ow^t`uvKgtUa;zUhttxr$CdeFRq&po1Uml!XBK5!6@dVUyJ|rDwnnygXuqyIswU3LV5)8iGkrjh&pHw zK$t414aoogmisbOea*eiso2SxC^k+B4f?ONw9Un8gag`(Un`%B!4~j%4ZVBN1&sY~ zzsx-dRJJ|i+P%d)r#EAlsdKqp^1ONT%++d#gPA*`9v-1mslEXLF2|1Lc1bWn6Lf>F z5p6p{&uBu(vy&DZGV=?F`Z)swzLg#utQ@l!J5JM-3U=j^7@ zcV2Dy+4j}e^83O5A6ekIfHYx{!N*`?EQrL5w;)88wF}5AQl1Sv#pr&Sw5{sdkQnl` z@1On2jmG|VJ5gK%Lag$d_n%gZU-TEEfCOZ>GTK)V)UmeFQZ-+ojWAcUx#PeQik1<9 zpS`_0yi7+&CzZ|S0{{6bs8E=HVMl0syqYljN0hj&?Zu$Hu%LhhC|xST5?TDnasMu; zXjp4(yx-35!f<8Z3V~4APSL{W0&8a-UTxd=FSViRv72LKl9NV(T-g@P9IXB^pU*d5 zrat)YBuA{ON6&$3pgo9kSq?9*TZWsMxt?Jk=WVpK1YTK61OT`kN;e#l5jEn_kJEf& zl2>diEOsv+m>B7an31<_BoGQq7_!qk%a-v)7Wxg(uNldQ${KyaR-nPpZ#rFP91TS1 z;?B+@lwqg>()Tgok9U8OZdNWb`Jj?wHCvG#-m6Wi2z)f`#5@Ni2_e(6@AJ?S>sg7lV$LDC`oTX6i(h3zmy6cWS@PIUm zKxFb#sVK+L%i{t4u)EOj4QF=DNFtHEqp0IoaT@P5H93i15LH0Q1nNL)E03ZZ*&AHR z2j)e{!@PW2ACoXi-%jd&n7JG0UhXf33X&@oKd#SdzIukHL|mJZUjT3W2cU@Mb>Sgy zOKa=zG}0e=;2PawOn>W#QiW2{67L+G3-1)2?g6P#_)8rQCvIvwe&Em_kN^Bwx*}@g z7-RCqkyaqSm6@tCu_y)e`N1^LK6|SVZA}Km;?J^#kZ9A6UoUYXlP`ahb=-jy4T-IG zin)1Rh!v-8aA>Hd?cpqy0Yt!WAp7Q5xD3-9k2VJwDGBt6(yrSWR^GR7q1(n@+Ymye z%xce)F3C~#ZZ$$L|NVPUN=*c{zRkq5+7*>%$)QxJ9-s7fbYvt5YEpb^;}=eWe}M0} zZJ@kwe6&AQw2B0YR=2X>uMAg@+HV{0JR4j&p}W?}Y{~pJW@Zj%HDfF|h1qXMfi&4+ zW@d&E{OJ? zQB&!OQu%H)hhYBaZF2gj@oBHWXjG~X&|pUZEoMtgOI~}sn=V(`0M~b9vfhOsEZH*y z%nH5fmVRQm2)yPt)8Ay*mb}|p4V5(8CO*W$quBvq0Ey#%N%|3}^{KhtNblVcnm+gO z+&M~=>ag|0@ACJtTi(9CyGxIedTVLLnuZ)hLlch*{P6Vm?{rcVYZdvIzEaA`lXqH+ z4&MUK6STei0|EkuPJ(;A$vJ*6I^1Z{T<7AZ<+RS+n#ghAz`*@bmvvucK}4r#w9%@> za98r(oE#!#VOe=MZyHE?d)qyhn^UZ_cvx$Cu{%! literal 0 HcmV?d00001 diff --git a/variance_hydrogen.f90 b/variance_hydrogen.f90 new file mode 100644 index 0000000..8187cfe --- /dev/null +++ b/variance_hydrogen.f90 @@ -0,0 +1,57 @@ +program variance_hydrogen + implicit none + double precision, external :: e_loc, psi + double precision :: x(50), w, delta, energy, dx, r(3), a(6), norm, s2 + integer :: i, k, l, j + + a = (/ 0.1d0, 0.2d0, 0.5d0, 1.d0, 1.5d0, 2.d0 /) + + dx = 10.d0/(size(x)-1) + do i=1,size(x) + x(i) = -5.d0 + (i-1)*dx + end do + + delta = dx**3 + + r(:) = 0.d0 + + do j=1,size(a) + energy = 0.d0 + norm = 0.d0 + do i=1,size(x) + r(1) = x(i) + do k=1,size(x) + r(2) = x(k) + do l=1,size(x) + r(3) = x(l) + w = psi(a(j),r) + w = w * w * delta + + energy = energy + w * e_loc(a(j), r) + norm = norm + w + end do + end do + end do + energy = energy / norm + + s2 = 0.d0 + norm = 0.d0 + do i=1,size(x) + r(1) = x(i) + do k=1,size(x) + r(2) = x(k) + do l=1,size(x) + r(3) = x(l) + w = psi(a(j),r) + w = w * w * delta + + s2 = s2 + w * ( e_loc(a(j), r) - energy )**2 + norm = norm + w + end do + end do + end do + s2 = s2 / norm + print *, 'a = ', a(j), ' E = ', energy, ' s2 = ', s2 + end do + +end program variance_hydrogen diff --git a/variance_hydrogen.py b/variance_hydrogen.py new file mode 100644 index 0000000..9fcf564 --- /dev/null +++ b/variance_hydrogen.py @@ -0,0 +1,36 @@ +import numpy as np +from hydrogen import e_loc, psi + +interval = np.linspace(-5,5,num=50) +delta = (interval[1]-interval[0])**3 + +r = np.array([0.,0.,0.]) + +for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: + E = 0. + norm = 0. + for x in interval: + r[0] = x + for y in interval: + r[1] = y + for z in interval: + r[2] = z + w = psi(a, r) + w = w * w * delta + El = e_loc(a, r) + E += w * El + norm += w + E = E / norm + s2 = 0. + for x in interval: + r[0] = x + for y in interval: + r[1] = y + for z in interval: + r[2] = z + w = psi(a, r) + w = w * w * delta + El = e_loc(a, r) + s2 += w * (El - E)**2 + s2 = s2 / norm + print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")