From 8bb6834f8724266e5fd045c58a3d6f3b3b708957 Mon Sep 17 00:00:00 2001 From: cophus Date: Fri, 27 Dec 2024 11:34:14 -0800 Subject: [PATCH] More content --- algorithms.md | 256 +++++++++++++++++++++++++- code.md | 87 ++++++++- figures/gaussian_plotting_example.jpg | Bin 0 -> 12759 bytes 3 files changed, 333 insertions(+), 10 deletions(-) create mode 100644 figures/gaussian_plotting_example.jpg diff --git a/algorithms.md b/algorithms.md index 9a53287..1baec9e 100644 --- a/algorithms.md +++ b/algorithms.md @@ -12,7 +12,7 @@ label : algorithms_page (numerical-solutions-of-the-schrodinger-equation)= -### Numerical Solutions of the Schrödinger Equation +## Numerical Solutions of the Schrödinger Equation As discussed in [](#physics_page), the [Schrödinger equation](wiki:Schrödinger_equation) typically cannot be solved analytically in complex systems. Therefore, in order to perform electron scattering simulations, we must calculate numerical solutions of [](#eq:Schrodinger_time) for electron waves. First, we define the {cite:t}`debroglie1925recherches` wavelength of a free electrons (corrected for relativistic effects) as @@ -48,7 +48,7 @@ where {math}`{\nabla_{xy}}^2 = \partial^2/\partial x^2 + \partial^2/\partial y^2 (multislice-method)= -### The Multislice Method +## The Multislice Method By a wide margin, the most common algorithm used for electron scattering simulations is the multislice method, first described by {cite:t}`cowley1957scattering`. Equation [](#eq:Shrodinger_electron) shows the overall numerical recipe we will use; when the wavefunction {math}`\psi_0(\bm{r})` is at position {math}`z_0`, we will evaluate the operators on the right hand side over a distance {math}`\Delta z` to calculate the new wavefunction {math}`\psi(\bm{r})` at position {math}`z_0 + \Delta z`. {cite:t}`kirkland2020` gives the formal operator solution to [](#eq:Shrodinger_electron) as @@ -97,7 +97,7 @@ Instead, we solve it numerically by using a split-step method, where we alternat The steps of the multislice method are detailed below. -#### 1 - Atomic Coordinates +### 1 - Atomic Coordinates We first generate a set of atomic coordinates for our desired sample. The atomic coordinates are placed in a *simulation cell*, a [rectangular cuboid](#wiki:Rectangular_cuboid) where all edges vectors are $90^\circ$ apart. @@ -106,7 +106,7 @@ For each atom, we define the $\bm{r}=(x,y,z)$ position, the atomic number, a the Ideally the atomic coordinates will be periodic in the $(x,y)$ plane, though this is not always possible. -#### 2 - Potential Slices +### 2 - Potential Slices Next, we calculate the potential $V(\bm{r})$ for the sample. We compute this potential numerically, either using the parameterization approach shown in [](#isolated-atomic-potentials), or using a DFT calculation as described in [](#dft-potentials). We divide the atomic potentials into *slices*, which are thin sections of the sample in the $(x,y)$ plane. @@ -119,14 +119,14 @@ We can also add additional electrostatic or electromagnetic fields to the potent The effect of both extrinsic and intrinsic magnetic fields can be calculated using the [Aharonov–Bohm equation](#wiki:Aharonov–Bohm_effect). -#### 3 - Initial Wavefunctions +### 3 - Initial Wavefunctions Next, we define the intitial condition of the electron beam wavefunction $\psi(\bm{r})$, described in Sections `insert sections`. In an ideal plane wave TEM or diffraction pattern simulation, we use only a single initial wavefunction. We can also include spatial coherence in a [TEM simulation](#tem_sims) by performing a multislice simulation where the initial probe is tilted to a range of incident probe angles, which are then summed incoherently to generate the simulation output. For a [STEM simulation](#stem_sims), we may need to calculate thousands or even millions of initial conditions for the electron probe, as each unique STEM probe position requires another simulation. -#### 4 - Transmission Operator +### 4 - Transmission Operator Following {cite:t}`kirkland2020`, if we assume a slice is infinitesimal thickness, we can set the ${\nabla_{xy}}^2$ term from [](#eq:Shrodinger_simple) to zero and obtain the solution ```{math} @@ -140,7 +140,7 @@ Following {cite:t}`kirkland2020`, if we assume a slice is infinitesimal thicknes We see from this expression that as the electron wavefuncton passes through a given slice, it will pick up a forward phase shift proportional to $V_{\Delta z}(\bm{r})$. This first Born approximation is quite accurate for high accelerating voltages, for small-to-intermediate atomic number species, and for thin slices. However we may require a more accurate expansion and / or numerical slicing of individual atomic potentals when using very low accelerating voltages or for calculating scattering from high atomic number species. -#### 5 - Propagation Operator +### 5 - Propagation Operator Next, we need to *propagate* the electron wave from one slice to the next by using the propagation operator in Equation [](#eq:Shrodinger_simple). We assume empty space between slices, setting $V(\bm{r})=0$ in [](#eq:Shrodinger_simple) to get @@ -234,13 +234,13 @@ We can now write the final propagation operator by combining ${k_x}^2+{k_y}^2=|\ If there are still remaining slices that the electron wave has not passed through, we alternate steps 4 and 5 until the prope wavefunction reaches the output surface of the sample, where it is referred to as the `exit wave`. -#### 6 - Transfer Function +### 6 - Transfer Function After we have calculated the exit wave, we then need to apply the effects of our microscope optics to this wave and reach the detector plane by using a microscope transfer function (MTF). The MTF could be very simple; for example in either a TEM diffraction simulation or a typical STEM simulation we assume the detector is placed at the far field limit, and therefore only need to Fourier transform the exit wave to reach the detector plane. For a TEM imaging simulation, we typically use a contrast transfer function (CTF) for the MTF. The CTF can include aplanatic [optical aberrations](wiki:Optical_aberration) such as defocus, spherical aberration, astigmatism, and higher order coherent wave aberrations. It can also include more complex optical affects such as field distortion, image rotation, or planatic aberrations, where the aberrations vary as a function of position. The CTF equations are described in `add section link`. -#### 7 - Detector Functions +### 7 - Detector Functions Finally, we convert from the complex wavefunction to a real-valued detector measurement. This intensity measurement may be performed in real space for near-field imaging giving $I(\bm{r})$, or in Fourier space for far-field diffraction space measurements giving $I(\bm{k})$. The measured intensity for a pixelated detector is just the magnitude squared of the wavefunction $|\psi(\bm{r})|^2$ or $|\psi(\bm{k})|^2$. To simulated an integrating detector intensity $I_D(\bm{k})$, such as those for BF or DF STEM measurements, we apply a detector function $D(\bm{k})$ to our measured intensity using the expression ```{math} @@ -267,3 +267,241 @@ I(\bm{R},n) ``` where $\theta$ is the annular coordinate and $k'$ is the radial coordinate for $\bm{k}$-space. + + +## Bloch Wave Simulations + +While the multislice method is widely used for electron scattering simulations due to its flexibility and scalability, the Bloch wave method provides an alternative approach that is especially efficient for small, periodic structures. +The Bloch wave formalism leverages the translational symmetry of the crystal lattice to express the electron wavefunction as a sum of periodic Bloch states, significantly reducing computational effort for crystalline samples. +This reduction in computational cost is possible because in a Bloch wave simulation, we only consider a small number of scattering vectors $k$, or equivalently to a small number of scattering angles $\alpha = \lambda k$ + +The Bloch wave method is particularly advantageous for periodic systems, as it reduces the computational domain to a single unit cell and directly incorporates crystal symmetry. However, it is less suited for non-periodic systems or those with large-scale defects, where the multislice method is more appropriate. +This approach requires careful numerical handling of eigenvalue decomposition and the summation over a sufficiently large number of reciprocal lattice vectors to ensure convergence. + + +### 1 - The Bloch Wave Ansatz + +When applying Equation [](#eq:Shrodinger_electron) to crystalline samples, $V(\bm{r})$ is periodic in the lateral plane. +This periodicity of the crystal allows the wavefunction inside it to be expressed as a sum of plane waves: +```{math} +\psi(\bm{r}) + = +\sum_n + \alpha_n b_n(\bm{k}_n, \bm{r}), +``` +where $b_n(\bm{k}_n, \bm{r})$ is the $n^{\rm{th}}$ wavefunction with the associated coefficient $\alpha_n$, and the set of these wavefunctions forms a complete basis. +In a Bloch wave expansion, we assume that each of these wavefunctions will satisfy the Schrödinger equation inside the specimen. +The set of Bloch wave coefficeints coefficient $\alpha_n$ can represent any arbitrary wavefunnction, but only one particular set of $\alpha_n$ values will match the incident wavefunction at the sample surface. + + + + +```{math} +\psi(\bm{r}) += +\sum_{\bm{g}} c_{\bm{g}}(z) \exp(2 \pi \ii (\bm{k_0} + \bm{g}) \cdot \bm{r}), +``` + +where $\bm{k_0}$ is the incident wavevector, $\bm{g}$ are the reciprocal lattice vectors, and $c_{\bm{g}}(z)$ are the expansion coefficients describing the contribution of each Bloch state. + +### 2 - Bloch Waves in the Schrödinger Equation + +Substituting the Bloch wave ansatz into Equation [](#eq:Shrodinger_electron) and projecting onto a specific plane wave $\exp(2 \pi \ii (\bm{k_0} + \bm{g}) \cdot \bm{r})$, we isolate the evolution equation for the coefficients $c_{\bm{g}}(z)$: + +```{math} +\frac{\partial}{\partial z} c_{\bm{g}}(z) += +-\ii \lambda \pi |\bm{g}_{xy} + \bm{k}_{xy}|^2 c_{\bm{g}}(z) ++ +\ii \sigma \sum_{\bm{g}'} V_{\bm{g} - \bm{g}'} c_{\bm{g}'}(z), +``` +where $V_{\bm{g} - \bm{g}'}$ are the Fourier coefficients of the crystal potential $V(\bm{r})$. The first term accounts for the propagation of each Bloch wave, while the second term represents the coupling between Bloch states due to the periodic crystal potential. +In this context, coupling refers to how the periodic crystal potential causes scattering between plane waves with different reciprocal lattice vectors $\bm{g}$ and $\bm{g}'$, redistributing the electron wave amplitude among the Bloch states. + +### 3 - Matrix Formulation + +The system of coupled differential equations for $c_{\bm{g}}(z)$ can be written in matrix form: + +```{math} +\frac{\partial}{\partial z} \bm{c}(z) += +\bm{H} \bm{c}(z), +``` +where $\bm{c}(z)$ is the vector of coefficients $[c_{\bm{g}1}(z), c_{\bm{g}2}(z), \dots]$, and $\bm{H}$ is the Hamiltonian matrix with elements: + +```{math} +H_{\bm{g}, \bm{g}'} += +-\ii \lambda \pi |\bm{g}_{xy} + \bm{k}_{xy}|^2 \delta_{\bm{g}, \bm{g}'} ++ +\ii \sigma V_{\bm{g} - \bm{g}'}. +``` + +### 4 - The Input Wavefunction + +The initial wavefunction $\Psi_0(\bm{r})$ is incorporated into the Bloch wave expansion by projecting it onto the reciprocal lattice plane waves: + +```{math} +c_{\bm{g}}(0) += +\int \Psi_0(\bm{r}) \exp(-2 \pi \ii \bm{g} \cdot \bm{r}) \, d\bm{r}. +``` +For an incident plane wave with wavevector $\bm{k}_0$, only the term corresponding to $\bm{g} = -\bm{k}_0$ is nonzero. For more complex wavefunctions, such as convergent STEM probes, the projection accounts for the full angular and spatial distribution of the initial beam. + +These coefficients form the initial condition for the matrix differential equation governing the propagation of the Bloch wave coefficients. + +### 5 - Propagation of the Wavefunction + +To propagate the wavefunction, we solve the matrix differential equation: + +```{math} +\bm{c}(z) += +\exp(\bm{H} z) \bm{c}_0, +``` +where $\bm{c}_0$ is the vector of initial coefficients from Section 4. The matrix exponential is computed via eigenvalue decomposition: + +```{math} +\bm{H} = \bm{U} \bm{\Lambda} \bm{U}^{-1}, +``` +where $\bm{\Lambda}$ is a diagonal matrix of eigenvalues $\lambda_i$, and $\bm{U}$ contains the corresponding eigenvectors. The matrix exponential becomes: + +```{math} +\exp(\bm{H} z) += +\bm{U} \exp(\bm{\Lambda} z) \bm{U}^{-1}, +``` +where $\exp(\bm{\Lambda} z)$ is a diagonal matrix with elements $\exp(\lambda_i z)$. The eigenvectors in $\bm{U}$ represent the Bloch waves, and the eigenvalues $\lambda_i$ describe their phase and amplitude evolution through the crystal. + +### 6 - The Output Wavefunction + +At the sample exit ($z = t$), the wavefunction is reconstructed from the Bloch wave coefficients: + +```{math} +\psi(\bm{r}, t) += +\sum_{\bm{g}} c_{\bm{g}}(t) \exp(2 \pi \ii (\bm{g} + \bm{k}) \cdot \bm{r}), +``` +This reconstruction combines all Bloch waves, each modulated by its respective coefficient $c_{\bm{g}}(t)$, to form the final wavefunction. This exit wavefunction can then be used to calculate quantities such as diffraction patterns or real-space images, depending on the simulation objectives. + + + + + + + + + + + + + + + + + +Substituting the Bloch wave ansatz into Equation [](#eq:Shrodinger_electron) and projecting onto a specific plane wave $\exp(2 \pi \ii (\bm{g} + \bm{k_0}) \cdot \bm{r})$, we isolate the coefficients $c_{\bm{g}}(z)$: +```{math} +\frac{\partial}{\partial z} c_{\bm{g}}(z) += +-\ii \lambda \pi |\bm{g}_{xy} + \bm{k}_{xy}|^2 c_{\bm{g}}(z) ++ +\ii \sigma \sum_{\bm{g}'} V_{\bm{g} - \bm{g}'} c_{\bm{g}'}(z). +``` +Here, $V_{\bm{g} - \bm{g}'}$ represents the Fourier coefficients of the potential $V(\bm{r})$, capturing the coupling between different Bloch states due to the crystal potential. +In this context, coupling refers to how the periodic crystal potential causes scattering between plane waves with different reciprocal lattice vectors $\bm{g}$ and $\bm{g}'$, redistributing the electron wave amplitude among the Bloch states. + +### 3 - Matrix Formulation + +The coupled differential equations for $c_{\bm{g}}(z)$ can be written in matrix form as: +```{math} +:label: eq:mat_diff +\frac{\partial}{\partial z} \bm{c}(z) += +\bm{H} \bm{c}(z), +``` +where $\bm{c}(z)$ is a vector of coefficients $c_{\bm{g}}(z)$, and the Hamiltonian matrix $\bm{H}$ has elements: +```{math} +H_{\bm{g}, \bm{g}'} += +-\ii \lambda \pi |\bm{g}_{xy} + \bm{k}_{xy}|^2 \delta_{\bm{g}, \bm{g}'} ++ +\ii \sigma V_{\bm{g} - \bm{g}'}. +``` +The first term describes the propagation of individual Bloch states, while the second term represents the coupling between states due to the crystal potential. + + +### 4 - Solution via Eigenvalue Decomposition + +The solution to Equation [](#eq:mat_diff) is given by: +```{math} +\bm{c}(z) += +\exp(\bm{H} z) \bm{c}_0, + +``` +where $\bm{c}_0$ represents the initial coefficients determined by the incident wavefunction. To compute $\exp(\bm{H} z)$ efficiently, we perform an eigenvalue decomposition: + +```{math} +\bm{H} = \bm{U} \bm{\Lambda} \bm{U}^{-1}, +``` +where $\bm{\Lambda}$ is a diagonal matrix of eigenvalues $\lambda_i$, and $\bm{U}$ contains the corresponding eigenvectors. The matrix exponential is then: + +```{math} +\exp(\bm{H} z) += +\bm{U} \exp(\bm{\Lambda} z) \bm{U}^{-1}. +``` + + +### 5 - The Initial Wavefunction + +The Bloch wave method requires incorporating the initial wavefunction $\Psi_0(\bm{r})$ into the expansion coefficients $c_{\bm{g}}(0)$ at the sample entrance. This is achieved by projecting $\Psi_0(\bm{r})$ onto the reciprocal lattice plane waves $\exp(2 \pi \ii \bm{g} \cdot \bm{r})$: + +```{math} +c_{\bm{g}}(0) += +\int \Psi_0(\bm{r}) \exp(-2 \pi \ii \bm{g} \cdot \bm{r}) \, d\bm{r}. +``` +For a plane wave incident along a specific wavevector $\bm{k}_0$, the initial coefficients are nonzero only for $\bm{g} = -\bm{k}0$, simplifying the expansion. For a more complex wavefunction, such as a convergent STEM probe, this projection must account for the entire angular and spatial distribution of the probe. Once $c{\bm{g}}(0)$ is determined, it serves as the input to the matrix formulation described in the next section. + +### 6 - Propagation of the Wavefunction + +Next, we propagate the electron wave from one plane to the next by solving the matrix differential equation for $c_{\bm{g}}(z)$: + +```{math} +\frac{\partial}{\partial z} \bm{c}(z) += +\bm{H} \bm{c}(z). +``` + +The solution is obtained using the matrix exponential: + +```{math} +\bm{c}(z) += +\exp(\bm{H} z) \bm{c}_0, +``` + +where $\bm{c}_0$ represents the initial coefficients determined in Section 5. The matrix exponential can be efficiently computed using eigenvalue decomposition of the Hamiltonian $\bm{H}$ as described earlier. + +This step iteratively applies the propagation operator until the wavefunction reaches the desired thickness $z = t$, ensuring the effects of both the crystal potential and free-space propagation are included. + + + + + + \ No newline at end of file diff --git a/code.md b/code.md index defff16..cde454b 100644 --- a/code.md +++ b/code.md @@ -4,7 +4,9 @@ numbering: enumerator: 1.%s --- -`Python` has inarguably become the leading programming language in the scientific community, continually gaining popularity for the past two decades. +## Overview + +`Python` has unarguably become the leading programming language in the scientific community, continually gaining popularity for the past two decades. This is due to `Python` being friendly for beginners, owing to its less strict language structure (dynamic typing, automatic memory management, non-compiled language) and the plethora of tutorials and examples. It also benefits from being open-source, with numerous high quality libraries covering a broad range of applications (numeric calculations, image analysis, machine learning, etc.), including a number of packages dedicated to S/TEM. The diverse and broad scope of the python ecosystem allows for interoperability between different domains, and file formats as well as easily re-purposing existing implementations algorithms to new domains. @@ -30,3 +32,86 @@ The main libraries used in the code examples of this text are: Dynamically interpreted languages including `Python` are particularly attractive for domain experts and scientists trying out new ideas, but the performance of the interpreter can be a barrier to high performance. However, by making appropriate use of `Python` open-source numerical libraries including NumPy, CuPy, and Numba, it is possible to write a purely `Python`-based code that performs as well or even better than prior codes based on traditional languages such C++ or Fortran. That is indeed what *ab*TEM has been able to achieve, which has made it one of the fastest-growing and popular tools for S/TEM simulations, and our choice for this text. + + +## Simple Python Examples + +Most python scripts begin with importing necessary libraries. +A *library* is a collection of pre-written code that provides tools, functions, and methods to perform specific tasks, such as mathematical operations, data visualization, or file manipulation, without requiring the user to write everything from scratch. +We can import libraries and alias them (meaning we substitute as a short form name to reduce typing) using this syntax: + + +```python +import numpy as np +import matplotlib.pyplot as plt +``` + +`Numpy` is a very useful library for computing mathematical functions. For example if we define a coordinate system $r=(x,y)$ where both $x$ and $y$ have a range -5 to +5 in steps of 0.1, and we want to evaluate the function $f(\vec{r}) = \exp(-|r|^2/(2 \sigma^2)$, we could type + +```python +sigma = 2.0 +x = np.arange(-5,5,0.1) +y = np.arange(-5,5,0.1) + +ya,xa = np.meshgrid(y,x) + +f = np.exp( + -(xa**2 + ya**2) / (2*sigma**2) +) +``` + +where we have set $\sigma=2$. We have used the function `np.arange` to calculate a range from [-5,+5) in steps of 0.1, which does not include the final point at +5, instead ending at $5 - 0.1 = 4.9$. We have then used the function `np.meshgrid` to expand our 1D $x$ and $y$ vectors into 2D arrays, where the $x$ coordinates are repeated along the $y$-axis and vice versa. Finally we compute the exponential function using `np.exp`. Note the indentation inside the exponential function, which we use to improve readability. + +One of the key features of 'Numpy' arrays is that they can be *broadcasted*, which means smaller arrays are automatically expanded in size during operations to match the shape of larger arrays, without creating unnecessary memory copies or requiring explicit loops. +We can therefore reshape $x$ and $y$ to be a row and column vector, allowing us to calculate the same expression as above without needing to call `np.meshgrid`: + +```python +sigma = 2.0 +x = np.arange(-5,5,0.1)[:,None] +y = np.arange(-5,5,0.1)[None,:] + +f = np.exp( + -(x**2 + y**2) / (2*sigma**2) +) +``` + +Here the modify the $x$ vector, by using square brackets $[,]$ to access indices of the $x$ array, and specifically use `[:,None]` where the first `:` means "all elements" and places these along the 1st dimension (array rows), and the `None` (or alternatively `np.newaxis`) expands the array to become 2D, by adding a new 2nd dimension (array columns). This changes the $x$ array shape `x.shape` from $(200,)$ to $(200,1)$. The similar indexing applied to $y$ changes its array shape `y.shape` from $(200,)$ to $(1,200)$. Finally, when using the addition operator, `Numpy` then *broadcasts* these two arrays from shapes $(200,1)$ to $(1,200)$ to $(200,200)$, which gives the compatible array shapes and allows them to be added together. + +We might also want to plot the results of our calculation. We can use the `Matplotlib` plotting library to plot this image, for example with the code which produces the result shown in [](#fig:gaussian_plotting_example): + +```python +plt.figure(figsize=(4,3)) +plt.imshow(f, + extent=(-5, 5, -5, 5), + origin='lower', + cmap='viridis', +) +plt.colorbar(label='Amplitude') +plt.title("2D Gaussian Function") +plt.xlabel("x") +plt.ylabel("y") +plt.show() +``` + +```{figure} figures/gaussian_plotting_example.jpg +:name: fig:gaussian_plotting_example +Example `python` plotting of a 2D Gaussian distribution. +``` + +Note that here we have moved the default origin for arrays from the upper left corner down to the lower left corner, using the argument `origin='lower'`. This change of origin is a common source of errors for `Python` newcomers, as in mathematics generally (and `Numpy` specifically) we typically place the origin in the upper left corner, while most plotting functions generally (and `Matplotlib` specifically) place it in the lower left. We should remember however that the directions are only display and plotting conventions; internally the entries in our 2D `Numpy` arrays are always self-consistent when being accessed using the `[index_0,index_1]` convention. + + + +## Interactive Python Notebooks + + +Interactive 'Python' notebooks, such as those provided by 'Jupyter', are a versatile tool for scientific computation and exploration. +They allow users to combine code, text, mathematics (using LaTeX), and visualizations in a single, interactive document. +This format is particularly advantageous for iterative workflows where the ability to test new ideas, visualize intermediate results, and document findings is critical. +Each section of code, called a *cell*, can be executed independently, which makes debugging and refining specific parts of a larger workflow straightforward. +The downside of this organizational scheme is that we must carefully check that our code works correctly when all cells are executed sequentially. + +Tools like *ab*TEM, interactive notebooks provide an excellent platform to run simulations, tweak their input parameters, and visualize the results. +For example, users can modify the defocus or aperture size in a simulated STEM image and then immediately see how the results change. +This interactive approach not only makes learning the software more intuitive but also facilitates efficient experimentation, enabling users to focus on understanding the science. +The ability to share notebooks further enhances collaboration and reproducibility, as researchers can publish notebooks alongside their papers, allowing readers to verify the results and modify them for their own needs. diff --git a/figures/gaussian_plotting_example.jpg b/figures/gaussian_plotting_example.jpg new file mode 100644 index 0000000000000000000000000000000000000000..80e52b310a75021c570b2fedbd649cb7dfcb3b3f GIT binary patch literal 12759 zcmc(FbzIcl)9BI-f^;gSgmehfNOvu;lyrB8uAp>Dmvpnhg0OT+qaqSZDIi@EOM}E+ zeENBR@4ffE|J+^9e%RSLXTB3>X3osHnZEfBxT7quBo9DAK>=7H-+-I109gPwCMFgp z1~wKJ77h+JEQQ4<*@(PMd z$~wAwKz##4kfoKijjf%%gQu6bkFTG9z^m7xVc`*xQHe>(DXD4c8JPuzMa3nhW#tv` z>gr()jZMuh-95d1{R4wTAI8Qf;FD9+Gqa0J%PXsE>l>R}h=aqUBYl>m$?5wMj^y* zl6p10e!`spYdE3h{0GXojQ3QuM`C(gJqVj<@Dd%WGrtD>bL!z8=4|srocOah=E*Tf zgChNQXj%K-umZmY{iGiEb9$PvIp-CFp#$bNLKOko(hr>+#P21VDo4UQu=eL`_pkeR zx0PR;ezKkOC$xM^%>Poh8v7o~xxVdUI|Xo*DQ;|8YrJ?oKm23Z?6{iHPi!VER{g<# zzX6fY&)vP^G0u|YN0~;QgN*QjT4}TUo6pqIcrC;2DQ!9Zze-WOpxN?m$$hQsKtbT? zY40tb(d1Hxy~1AkTpcj9>&QN((SXQ}y;N-3qG_dLHK=mtw6a{~vLqbi^Ri*7!P&s@ z5hahbUK}R>sysfLIy%>wK3SXsi7?@(C-=Lb7+Ge3d;K-(~pL;lD!ONq;Keof{N%pL-;}4-*x=STV^LKPSZfc1vA6U zJP{<1CxqsrQ$@~)h91D4ijhISZGq_4%9qSw1dQ&SS_&47`ZQGPKSa)(X^$u+xdl*M zH4uV3^U+JuGmMQS_WEgVQe^5>7Mi+^V!6ha`;{iL0|Ck8rZhfjD&u!`$5VHmWv)f2 zXRO%6KN$18e-e3~mA40?=)mesfI*&1_7Yf^Y!e=QGrCH!@=o}|k^H8n>oFDOgKBI@ zVNVOZB%#E9)~3@HT_{y+RGBst*F`6ue)t<)v9F40B=P|DPt9i)RAaRxF>>Or0oWoZz@;b4C zZ|q_bTYbbR-E7!cI&m5TC|Hq~)`7*2!a_b2vySa;4m*xlx!&H~OpkrS8O2Tp#XObV zm!GSjFklYl&&;m(vN7?Qq+WeyT!6d3rCpD6kK<2xvHpJQ1~8jhKkL#lDHKo_$+bvc zaez5}Cw+z4`9s+j4OT*w8h_Wq%;=e%WPQ-cnL?uONJ4zQYsA3SUK-U%006ew?yA!| zDHuy(pduSsa?ZUD;^6ccV6WAIMY2oFVmV6 z9lU-iG*-lyD9i@DcOU^v!V2NKhC}BTw!^m84aX;kzW&H3r^TYuBN!(-RT{y{@dFfE zFOA$)uwJY((rXy&&#fGCMU!@6zAH1%rFBXoP1UbYYg^v$a3s1ODR`x$r-tDL2_z8R zto)QV#2hbHh7yEsjh*+_YWFh{na^^!>5JgS9qyg4;M(>2tl@*Yb$@_+%w=osmiwNj zaH!xDI(^NSFM_M%-zQjq%JtPxS}Z;-0(=hlNO^#EZp<m@hVkj<47$BpWoJPPxnqM3|5MxHR-7%y5|Eab&Y~rM=hhVx4dwcT8Ap^;{x~ z!{KL-=&Q2=@j%r{CRnqOy{%i)tp7)Ljx==hDM#z93D2>-ru|3ltKkukxVWwQtAZMb zn3ZaGiu7XGbdnOOc{~S-@PZ>ml!wiF$<)AzSA-@Q)g$vZ`_N@}aO zEp?h|W0&g{t1w#RjRCt-Hs5;jy#c6hA{7h!A>ck(^FH__?czA&#nj641wLTXuFFm< zU7#kGr5vdy%KufhYbrSc;Pny`C(IYlLj7jO*8-iki)g$|3kWDTyct>)7ltS`|UG4_(is=T>Uk+Z(-v$3gML!xv zB?_!p!VT#b9f```@#QGuvodvAA9%o~VJ-**Eb!1BHU$P`JQ@kw4fAm`!kA?JaAd(3^1 zX8m=!s%bw~4Y)_zqL}AIQ#|mR8uKrD4=x4B7Orq>ZUDbM=)t4p zP~?ODj!5_qL^>oxZ~jC?t-nR&{Ygpx9!U)67e{8lI7<45B>WH{d(w__)fe4ZwMY%g zn*61qA|yxQ&9@wJ|HhG!oIsQWL=}pZyz9RtpZiaxLHgOS3%BZfkb?iGo*DK>g&F>) zp+39fW02)VJ&^;FSY*@Rh{f;k_3tgWlWc1}7jaUhkYD&{c=>ISxGo!Vbpud`A2>fJ z`Y8b(OsfD-BgJ{7su#u>V@}iC-%i;rtWc#-=e;y+w;hljmsK#WLGbz`nu=$cnCcG; zrN>%$?va4)*;eYqZQwvuoBj>JoD#V};(QzJ0+|k8n76&Fm!&ehcHM)tzUv23Hvo+~ zS0Z~>mCYDpVu;eqZsDZ08-T)h3!&q6vZ>862a;LLiJhvdY7bAk`W-&Y+HYD!o5pJ6 zwtc3gH-MU$+zTal(^v9ZQRPuR>@eR}+*+-agx-x=y|^TVhJIZC+&#nkX#1kCx!Io% zS0lwW&<@Hk4PD90m*Vj^P6t`WT|jGwd9Ka;1C$RdB_o19O^XE;jzN3ASCJ67ZSxKU zo*c~Z(aA3cuheLls1-|F$Sk9Yv-=0J6$}?i7q@3UPn?(Kyl<0W`hYp{<2~oC7d9*O z;dC71{QP401=EffPAsvpVGw5WD1!uQ%4LsN;wjlL=~{fM-<5A>nxw{Bq=JG_DM?dj zahPjGn3EEfkjp4Rk^x$Br2if$rU%y}>Mzxv`(jkeg<$1KMH%d4YItq}IfgozKbEVv zBsaDY=0D2LlFAIK6=1!i8et0j6x^8GBw>o_e;y|;GJmkGOacnF_T`jLa49?#D`9%< ztq(}}ve-!tIA30#x=c25{fW_Y11PkC@zcPF=&e~L0f9i6^BG~+WS zsND5@dzf>@F|s|;1ID3`C0fiYMt#ZF1zhXGc#I1u4OUI~kexyjAK@fwdfv16^ zjb382B2V&`q!>p7M=s{uppW*C@@}05o%<6UKMHMO1_xJV$UlmD4$>kvV2syTGHC^w zvW~J|L2`F;v<`bp2^mDrfiKwr##jroCldjggXAOMB&acwv5|JOk!MhmCLi#sF($}5 zZn(~zUO#Qux$Z_f!2%0t&c7km#%5jZ>Q4Gy?HVBHDC7zs^g2F_9+Suf-@#L1TMAI4I^OB^jzBM~ZUj zI*^O>?>B%y9J~!iz<&X8z<-x{zSQt9VH;%%0d%A4AC~LcXEZQYDMOSY+7&vy$zc_` zQOU>38QK}ID~Q>3gV{Ax?!`sFJ*bdEs-Y5lvs+Z3*E?$7gx&;nuyL?yoa>~^JFcu( zM<+0#zRKPwlS0E{%%*Bp>S+gOwsRRv7AxJ=!JHAtqCU7 z!lK$a#@YuYxjwnGl_Vah~K*z4DZYa^<6JZ|l;dQuvD5i&e|DXZp0u;L2}U z6w9?My$Ka)x9>zCl~}`68nH$D*5%m$jtz#H8^HC`IsNPI*xQ{}EoR%*u39UJ zDsI!3il|-IAO%kynO&%01i<|<3#Ji`QIgvi+*d3iZDLOTM)9UsG8##$ifP@W$0m z-4J3>-V}2$pEWyk?-=HtajfXU7@jS>`)p&Y-`Mx-K*ix_c{!gWeZxt&4)lx%2&@R} zS6ZT<2p&rFqf@fJM`a=dhRH`sqdUl$Zf0jP8kHrT0!YP~L%DWPs5T53IWU?W#-M43whoS4IM8yw0sNpKx{l3ko5tOH!D|n^ zP>8?Z*pRq6T0i|Z3sUYv!>Xh(pO79rpt^imO)U4z?mFW4zZ{eC;aWsxV0h&#B{EKI z)c1mFqE%PR?{C0F9|@9IGEOwO>6Njws7K1M7(A^kX8Hutl~Yc!EG^J`9oM|BJ&>zi zi#wo=-Ci*!n*-Un+ z(Qa6~9Dnx)F!kfsCIeag8tSJJQ;|1-?So%N88zC*>6eiN!`h1UN*#$4jN)0&4dh~W)^vo}-;<3ElQ;FU z9c7+Yyn7Vjmh!nm5SI|CrRqVLVV}2LP8C)1TRmSg_dehoU^qlsPJs+9>fKVFyHa+I zwGY`|E;fUU=2W%~d@N*jZ%UOo_qP@RI3tW%z zo)!}?7W}bcWs`6}>fU?RMBeV09Kf#9cYub2ZGi*aL{{tFm{U>v1Xn>rLp)7!N=VrV zTCR4W9uc7<|8}aj@kjccJ3A7L?9lYpe@r+SM1nkT0Ey99S0RfxfYCbr?5Y6n0m;AF z)zZVIvreeQvVbPM`GeyixHQUN_EqaIsgYt3}H=c~KnK;8$@ zRo2D|la0G1?_6}Mp4(u;egw_4W@kk5EtM;BD0^3sL|uURT141d`I18_njnd|kbBi< z{P29nqT8@eS>Nv?eBcr>E3%$8@vQ{)@d&$n4y%DUXWH+M3cr z?n#lwLYQN35T_>st;g*6nEN_O&GhgoJpGrQ=6L!qUisIz8+DJ`=l}4glTn(|KK#P0 z`H?$w@lEz>zo9$|)~f{2Jm;7o(cO&>qfPhomKWt@pUB5mi(M^agIA7Bx@vlS(Zr~K zJS-Exo4WrwDH;Ln7F^c z6Ibj`_swIjcT~R4#tv?tPSY7!qu{rBE8kvZvfga)p2{Q;yzcw zIxOzIim`gZ_>xzQyI_DJ+Io^R=wS~1yPw33UWt#jl~%Yau??az^l+s`&Kn6Ov`pas zr>_OF48+U?uaok&0<42;?Mj@bT6U9LQl_hWj|h@to-YI(Q3DiT>F;U7w`}^9Vbl9_A4W>mumerY^pZJ{zC?`c#g;YS819kFvM;PB{usvH$&M@yFGr z8XeKvUs(LVdR0_Z&Mj!H03~M<7Yl)Nc-VLC)nm+e<(&4ao z+oZFffJtpM$e*kst#o#}Yt~8Kq9N2$yI)EAu_mS)mIS<)8+2UL*{xtL8B9)6$)bO| zsjW5ZQ*1;bbfCB-v8}$ctqE^!(MCGYPl)X_#wwZ@1$6~gc%FmMC&aa@oe|RwbQ~a* zv3HN+l_!In$PNrc0&*G%3YCfVEy?6rA+@leY_Ie9vK~Q;L0WB|A_7@qvdB4a3mgpX zWW8qU@=uNdIES>9sg<_IgalUK;$RRrKfxrAfT_vV&5fW2-Yu(d-h#Q0%$md13H5oA^siF)M(dWr3n^i`;pDIn!K7GPqdLhaLvW0S2OxNk{-z#vUk;7AwuIjNxQ zPh8>a$gno5mZ`R18(Q{{4ec0Bwa={$33>pT3O=h`ULuMFZ*Sd(i}r&ytiAz;J8eqH z*scZM6cFrFSQ2v@GJSO{+xE~Cf%WG$JU2yKyqgfI9yUKjx&bQ%sBF&0&ze_}k^vC1tlU<&`t(8kPR%i*; zG1JmVac;E_6>?k|>pT{KS!_4sq(lJKIw?x-qPX^g5rVfRsh=8H=j&?A?EcmMi@!8+ z4&Udl!IgtDELE$EDSu<%SI5z@I+{&}K^la`t8%IJqB=9muM@)We|0y--i!-Cn8X4| zwdDTG2wF#sTq|=j{BfrADa|~R{8tY0X8^Ef@3pDOhhRGvMHb9b!>B6k96s3~t!1#8 z3Zj)+5?0yoNL4rpB{`CZc72<)a@J=ISgfKc{7_ao{Lyz4RHHQ3e8#dl>>$mp?zqld z&*imPF_TrRG*-==5=%ON?96%!P3c`U8OiTi(uT5#c_9t zMJH2PA#1*TShzia(ut#aeTkQ?fFn|LQm^M2!vAs4Floq4FC^cQK{5j3qnq6Vr&2RW zsVIsylam*cCrz0#mDRHn6bplhR8KQa=9AzM=y?K}1@(rV;z3dOJoUa!>^ih&+UBFn z!1v4moqy#V`|{AhGqf834@cB(l$ zJK)K6UiyZS%Xt-zd0*t-3&A7I?DA=claOt-_QtkG=Cagzy#yROdFei~Yo0C+tw*D# zp8iV&{$A}yeY_6RoL)OJiWF)GcZzEIxRY1k#x8!Q3d2SVeq~bK2w{ya&A2Sd@u!t@ zb}1zgt2)&L^92Hwz5^ea7aC8+I_ zjqf_0g{@ZClEeJc^nsagQx}(G7(=mRVN!G`mLVSQ6PeLVsQ#&DUa@bgswd3$xIa2e0A`JUb5ML5g(F zD+E-f+O-FIAGC}<#G zxq_@J;2~g!v0EX61`5^o{p(4fj8f93Ds0`_*NyX;VB4LT+p{v~L{u_XikAPLvX+`45RA$*J{V_=C zYWrOr?tr+etW}-bo-RugYb49CoT;bbAIupmfQyH(HD;mgtH)5u||wPH_Ak)qGjj6#BY;J+u$`>cX{oqtcOe`V2n z*XRL#L5It-@ev)!stGi`^C*-0g1W;OQ9IDtf|ASuNOF6?pP> zq~V*v)msW%PP!E9xo&>^g!&wH>|v*GegO->#>|CY97jl(OUwyk#G%#P*j%4>C8lwB zDs#G^Ml{Z|d!Arij7_MN`Cal7HrpRZr zk8c1&{&XeTOi#}+QVff}RY02PtKgn-95t0Y`Z8m%Cb%~TwhOK&&AY{AiHe>zon7oa zLfOTuuzZ%Q`6|w|JBH;sy@tjPz7%oHaV9Z+*$sr2PUfB?rUPF%a=(ai1&MgF&uld& zZoYE`*M3Xj8ZtShAR%sBNgcG9dGvNl`XAnt zi#C|^_Tb`T0@uC!_^O+nXKxzo1He3ZoNPi9>l2s~X7t=_;LSH91+Rz}Cq+4UWy4Ma zbLsA2h2nZO@3}s1Njt-GC$jJh{$mRXVYV9!?obo#~xze-MksG>ivMsByRJ^XWRX*i$Q!mJ6wl z=@K6GCl=UNEu&Rr6Ps*dImicu&Ns$r$fVq-t7wIZ%;x!2k%&|4_CL|{bqIt>4Xf_H zmTxYj9G{+0bx=>wk8xbM27m{%D{cVOTAj$=#UVoQ2%kA#PaY`9Cl-Me@!m`Wd9VMPlca~2Nh~K zGuXnEUk<5pt37S{!fMaX8W7ViKcZXUB@=NH5So)YqFDzo^3f?0$;tD{&&jtb@sUC0 z*X5k3gx{rBP_rTGblcM~-xoY0`hjd7{awT(ucy9%>rO8tGoe?Z48KlSe#p~R z`qGPAsV^*#$45#arCz;{EVYQ~9fbzXAOOot`VpPmZITxnWS)zX5&nnzvX|(dhzC8n z$Me5r{05pKM(6#gG~SEq`@5Wzr#Z$(U7Am>m$csd=q*E;y`V(`+e3R*=u_0shhM)= zSTU}*9_ZYiNjom2xReTchm5Lbf&OLLXvTx^Gbe@q_e`kt0i}y;* zF{@2XWE&vZ>wwsVFvpJ&#gxXj)bVtK;2gX-sSD!qDo=9Lf@Tf0kM3fkZ)@jFzHJO6 zTYw}(E-QSBh}R~cs5)$9el{d3G3BIJSO4c*O6+Z?Cq7fale)S0WPL;C=-P^8tM{yS zuu{8LmMS%|h8O6;D+OPX7W`Z7C>vQuzHR^h)gG`S2uMekkw=kbWJLlo>ZrE_Uep4r zRoW0eH~X!LLA&}YZ4Nz6I3tXo?P&v48v^Hn-wUuI33d}#q(_=UdG)9Q%6`$L6Q**+1H!B(Qo2`W%GEwn9*p0M7`>`&kje>$Sdt(BTL z&qF(`#jnk;Yh6wwe{Wc;HV?*=DPN8v&~3+>rC9&+R_}X%)4Q4YT{D^Mk4(P%m#tb11p&cFwKH82xFL7*Lt!C}622fm~PBYj~22CdTmKZx+6v%pepM^X&NfIr~OJ6mU>5BD4teJlGCsxlLkpV$0TC)`{xybXy zTU@6npkP4>@=oDUO6%Djb&|&;S#Q6mSx60y6hs|6_P3}gfj&@{G^3rm58scaL@wX^ zLOxrts=fEkkN~mbzGTK-QB+GO;%21{H^sppYAR9;JV7>TI%xxLLJ)gLN68Z}Y3<3? zQtb|_rY9M9)J>IN8J;*?^sd=@)GhQ)=V{wBJ^>k zyi8}tL|l7;yGn|A(EQF{?Iexg+ev=V{2PGn2(rQSJEYkijCw4BFhe~Scj$O_mO6cw zyTov@Lv)?gXQnarG579Q2GStu8-R&9vO4Go$;G*pO-yWyl%N)R@!%}*x>fC$EsY-^ ztoLA!^f3@2(;~xP#uSG%Ce(`hkf?{%mOKSf4*~zpKwDdgH$ztAJAW_6cXuS81@C8? z(fY-S8w%Ff;y9#U>wu^{VcB)r+rQfDAHh9OjKN=~D6a7|ZUA)ETY`(&59{8g_{rn& z_pv$7b}nXDOKvzK%T5@jH-NZFpR9epWzA0b4tS9zO7h2cal_mRYS5=1+U9|bfs0R; zIT?GP=MEXAllD341WyE!p>F^4spOF$^14OH_zhsd9=w>1yyfxdg5wu(o{L^0O3b$e zf!vq;eIr){?%?nB$jchDl8ZT9f5RL-uCPF2(qO)&ZjbD=g%oRp7dfJXbI&nzub=1& z5CunWP@vBR+?e0MTy}D80q}A_n&6+6y=#IGlG;OT8rF7O3%(xn}r;!;ni1`s1MWgw4rz`0;1arFx2!=dJ8Smn@nGk{i_E zvmFMb68c%JH!i0|bEZ8v0Q__0JrP-5$+Jh2mp<|TqWjz>Jt6tTAFUQ;ctBgu6BK?ZqM=T9ZtUBG#-EU=>QDyjwQ zU-GyAO=%iuw#b`0e+l=wL;7_sTAw|~Wd0eIyKuywSNxfh1%nMNWPSK2%6VIViqG<= zWqN(%!`83aNXHwS5;-R8;Gd(3=Y4iin=IEaEPN^+I7t~b`qud_^>N|j8M@~|RS2zN&77zsR;rCf5e?S!*XXN(ycb4PA|djq({X*&O@=XkrdjKkel-eMLilN|H? zs?#8y{WIkBB=Dvf>~{96Ydn=jS8@Z#M88m%zI>1SM%>))b~DJ??x%I*umCvOTIrKK z2n(ADvKMcV5&GXW<^OEd>$bPP8{@2Y-OYL8cgH_Tsj(VB^>`%*KsN#aTCa9SgARbT zZ~?_Cl@~h?#*2N?=3Z~%Tw+*mgTT>TEaUTb0@><#dj3VRHi^U3L6@=hl|`Yn!MW|P zrZusQN*%~s8nT?r3_h;Se(QhN!#dWQ9P=s!g(23kT4gcmTW60m+76EbF}W7@iL#l8 z6Q$)r%-zk;;%eUYzq(04`og8>L{iCybE^CPT4Cm7m2q!MEj~r*VbKk4YAQJ|%W|Z{ zo}*ZP^9dVUVOJ=^Nn0(7oY3$_(Q*>LdC~HQKOs~jN z$A0SNtWCU2Ss7;7wu*$42fFa^iRZN_<&fMO*)U@BFZTq#7iwl$6M|=l$pK!e9)`WUrrct}! zJ-btAuahO9gYe4$AZI2Ie@3vKk%&aD5Y|lqT4La}AWrNT908)Coz>tR_c`b%gFt6V z>WGmraD9EdIMju?HFa%GTsr-Oq2`lky+W@wQA;!c_;IE1jsGri{LkSvFan3Q@qQ3j zJoQuOy}af`;4|M>lxj^P0|?p`YOx1Z#@n-ceekiT@pQ+-HN>cx)2miGLt@|YzhLhC OpX~oX{LXMQ^S=P?(r|(R literal 0 HcmV?d00001