{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Convergence classes\n", "\n", "BigDFT uses wavelets on a discretization grid to represent the wavefunction and density of the electronic system. The key parameters of a BigDFT calculation therefore are the grid spacing (noted `hgrids`) and the grid extension (noted `rmult`). The convergence of a given BigDFT calculation with respect to those two parameters can easily be studied via two classes: [HgridsConvergence](https://mmoriniere.gitlab.io/MyBigDFT/hgrids_convergence.html) and [RmultConvergence](https://mmoriniere.gitlab.io/MyBigDFT/rmult_convergence.html). Both will be used here in the case of the N$_2$ molecule with the LDA exchange-correlation potential." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import the relevant modules, classes and constant.\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from mybigdft import Posinp, Atom, InputParams, Job\n", "from mybigdft.globals import EV_TO_HA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## rmult convergence\n", " \n", "The `RmultConvergence` class allows to run all the necessary calculations to determine the smallest grid extensions which must be used so that the energy error compared to the reference calculation (*i.e.*, with the largest `rmult` considered) lie within the required precision per atom. Using the smallest `rmult` as possible allows to save computational time because there are less grid points (or degrees of freedom). \n", "\n", "Note that there actually are two grids centered on the atoms of the system: a fine one with a shorter extension and a coarse one with a longer extension. They respectively correspond to the two numbers in the `reference` list given when initializing the `RmultConvergence` instance. These numbers will be decreased by 1 until the resulting energy lies outside the given `precision_per_atom` (here, 10 meV). The maximal number of jobs to be run is fixed to 6." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from mybigdft.workflows.convergences import RmultConvergence\n", "\n", "atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.1])]\n", "pos = Posinp(atoms, units=\"angstroem\", boundary_conditions=\"free\")\n", "base = Job(posinp=pos, run_dir=\"N2/rmult_convergence\", name=\"N2\")\n", "rmc = RmultConvergence(base, reference=[9, 11], delta=[-1, -1], n_jobs=6,\n", " precision_per_atom=0.01*EV_TO_HA)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The jobs in the queue have an extra attribute `param` which returns the value of rmult:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[9.0, 11.0]\n", "[8.0, 10.0]\n", "[7.0, 9.0]\n", "[6.0, 8.0]\n", "[5.0, 7.0]\n", "[4.0, 6.0]\n" ] } ], "source": [ "for job in rmc.queue:\n", " print(job.param)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can run the calculations as usual:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/rmult_convergence/9.0_11.0\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/rmult_convergence/8.0_10.0\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/rmult_convergence/7.0_9.0\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/rmult_convergence/6.0_8.0\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/rmult_convergence/5.0_7.0\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/rmult_convergence/4.0_6.0\n", "Logfile log-N2.yaml already exists!\n", "\n" ] } ], "source": [ "rmc.run(nmpi=6, nomp=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Once the calculations are performed, a converged job is determined. It corresponds to the one with the smallest `rmult` so that the total energy of the system is below the convergence threshold, defined by the sum of the energy found for the largest `rmult` and the user defined precision per atom. It can be accessed via the `converged` attribute:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[5.0, 7.0]" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rmc.converged.param" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result of the previous cell shows that `rmult = [5.0, 7.0]` gives converged results for the given precision per atom. To make sure that everything ran as expected, all the relevant data can be accessed via the `summary` method:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requested precision per atom: 3.67e-04 (Ha)\n", "----------------------------------------------------\n", " rmult precision_per_atom (Ha) is_converged \n", "----------------------------------------------------\n", " [9.0, 11.0] 0.00e+00 True \n", " [8.0, 10.0] 1.09e-06 True \n", " [7.0, 9.0] 3.71e-06 True \n", " [6.0, 8.0] 4.96e-06 True \n", " [5.0, 7.0] 2.75e-04 True \n", " [4.0, 6.0] 2.31e-03 False \n" ] } ], "source": [ "rmc.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the convergence\n", "\n", "The figure below shows the convergence of the total energy of the N$_2$ system with respect to the rmult. The threshold is also shown for comparison." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAERCAYAAABCcWF4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlYVfW+x/H3VpxnygE3Jgo4MSrgkJqobNEy1HIqj2JplJ2T5/YUDqdTqc81zc5tVDMsE+/pqNkhqCwkVNJKRVQqswxNb7JFUsAZB3DdP7zuGyGCwGInfF7P4/Owfvu31vr+kvbHNf2WxTAMAxEREZPUcnYBIiJSvSloRETEVAoaERExlYJGRERMpaARERFTKWhERMRUCpoyWLduHT4+PtSqVYu0tLQS+7322mv4+vri4+PDq6++6mhPT0+nd+/eBAYGEhwcTGpqKgCnTp3i3nvvJSAgAB8fH959991Sa3n44Ydp1aoVvr6+FR+YiEgVUNCUga+vL3Fxcdx1110l9tm7dy/Lly8nNTWVb775hk8++YQDBw4AMGPGDJ5//nnS09OZN28eM2bMAGDJkiV069aNb775hpSUFJ566ikuXbp0w1omT55MYmJi5Q1ORMRkCpoy6Nq1K507d75hnx9++IFevXrRsGFDXFxcGDBgAHFxcQBYLBZOnz4NXD2Kadu2raP9zJkzGIbB2bNncXV1xcXFBYCXXnqJkJAQ/P39ef755x37ueuuu3B1dTVjmCIipnBxdgHVha+vL8888ww5OTk0aNCATz/9lODgYABeffVVwsPDefrpp7ly5Qpff/01AH/5y1+IiIigbdu2nDlzhrVr11KrVi2SkpLIyMggNTUVwzCIiIhgy5YtNzyiEhH5o1LQ/J+wsDCOHTtWrH3+/PmMGDGi1PW7du3KzJkzGTJkCI0aNSIwMJDatWsD8Oabb/LKK69w//338/777zNlyhSSk5PZsGEDgYGBbNq0iYMHD2Kz2ejfvz9JSUkkJSXRvXt3AM6ePUtGRoaCRkRuSQqa/5OcnFzhbUyZMoUpU6YA8Le//Q13d3cAYmNjee211wAYM2YMU6dOBeDdd99l1qxZWCwWvLy86NChAz/++COGYTB79mweffTRCtckIuJsukZTiX799VcAfvnlF+Li4njwwQcBaNu2LV988QUAmzZtwtvbG4A77riDjRs3ApCdnc3+/fvp2LEj4eHhrFixgrNnzwJgt9sd2xYRueUYTvTZZ58ZnTp1Mjw9PY0FCxYU+/zChQvG2LFjDU9PT6Nnz57GoUOHHJ+98MILhqenp9GpUycjMTGxzNssj7i4OMNqtRp169Y1WrVqZQwZMsQwDMOw2+3GsGHDHP369etndO3a1fD39zeSk5Md7Vu3bjV69Ohh+Pv7Gz179jTS0tIc69tsNsPX19fw8fEx/vu//9uxzquvvmr4+voavr6+Ru/evY0DBw4YhmEY48ePN9q0aWO4uLgYVqvVePvttytljCIiZnFa0BQUFBgdO3Y0Dh48aFy8eNHw9/c3vv/++yJ9lixZYjz66KOGYRjG6tWrjbFjxxqGYRjff/+94e/vb1y4cMH4+eefjY4dOxoFBQVl2qaIiFQtp506S01NxcvLi44dO1K3bl3Gjx9PQkJCkT4JCQlERkYCMHr0aDZu3IhhGCQkJDB+/Hjq1atHhw4d8PLyIjU1tUzbFBGRquW0mwHsdjvt2rVzLLu7u7Njx44S+7i4uNCsWTNycnKw2+307t27yLp2ux2g1G1eExMTQ0xMDAA//vgjXbp0qZyBiYjUEIcPH+bEiROl9quxd51FRUURFRUFQHBw8A2nlhERkeKuPStYGqedOrNarRw5csSxnJmZidVqLbFPQUEBp06d4rbbbitx3bJsU0REqpbTgiYkJISMjAwOHTrEpUuXWLNmDREREUX6REREEBsbC8AHH3zAoEGDsFgsREREsGbNGi5evMihQ4fIyMigZ8+eZdqmiIhULaedOnNxcWHx4sWEh4dTWFjIww8/jI+PD8899xzBwcFEREQwZcoUJk6ciJeXF66urqxZswYAHx8fxo4dS7du3XBxcWHJkiWOp/Cvt00REXEei2EYhrOLcDZdo5Fb3eXLl8nMzOTChQvOLkWqofr16+Pu7k6dOnWKtJf1u7PG3gwgUp1kZmbSpEkTPDw8sFgszi5HqhHDMMjJySEzM5MOHTqUaxsKmnKK32PnpQ37OXoyn7bNGxAd3pmR3XXjgTjHhQsXFDJiCovFwm233cbx48fLvQ0FTTnE77EzO+478i8XAmA/mc/suO8AFDbiNAoZMUtFf7c0qWY5vLRhvyNkrsm/XMhLG/Y7qSIRkT8uBU05HD2Zf1PtIjVB7dq1CQwMJCAggB49ejhe8Hf06FFGjx5d6voeHh74+fnh5+dHt27d+Pvf/+64ueHw4cM0aNCAwMBAx5+33nrL8XPdunXx8/MjMDCQWbNmmTrOa1auXMlf/vIXAOLj49m3b1+V7PdWpFNn5dC2eQPs1wmVts0bOKEakZtnxjXGBg0akJ6eDsCGDRuYPXs2X3zxBW3btuWDDz4o0zY2b97M7bffztmzZ4mKiuLRRx91PEvn6enp2P41197Z5OHh4Vj3RgoKChyvS69M8fHxDB8+nG7dulX6tqsDHdGUQ3R4ZxrUqV2krUGd2kSHd3ZSRSJld+0ao/1kPgb/f40xfo+90vZx+vRpWrRoAVw9GvH19QXg/PnzjmfgRo0aRa9eva57e2zjxo1ZtmwZ8fHx5ObmVqiWlStXEhERwaBBgxg8eDAAL730EiEhIfj7+/P8888DcO7cOe655x4CAgLw9fVl7dq1wNUQuzafV1paGqGhoUW2//XXX/PRRx8RHR1NYGAgBw8e5PXXX6dbt274+/szfvz4CtVfHeiIphyu/ctPd53JrehG1xgr8jucn59PYGAgFy5cICsri02bNhXrs3TpUlq0aMG+ffvYu3cvgYGBJW6vadOmdOjQgYyMDFq3bs3Bgwcd/fv27cuSJUvKXNvu3bv59ttvcXV1JSkpiYyMDFJTUzEMg4iICLZs2cLx48dp27Yt69evB+DUqVNl2vadd95JREQEw4cPd5wiXLhwIYcOHaJevXqcPHmyzHVWVwqachrZ3apgkVuSWdcYf3vqbNu2bUyaNIm9e/cW6fPll1/y17/+FQBfX1/8/f1vuM3fPk9+vVNnZWWz2XB1dQUgKSmJpKQkunfvDsDZs2fJyMigf//+PPXUU8ycOZPhw4fTv3//cu0LwN/fnwkTJjBy5EhGjhxZ7u1UFzp1JlLDlHQtsTKvMfbp04cTJ05U6NmLM2fOcPjwYTp16lTheho1auT42TAMZs+eTXp6Ounp6Rw4cIApU6bQqVMndu/ejZ+fH3//+9+ZN28ecHW6rCtXrgCUeeaF9evX8+c//5ndu3cTEhJCQUFBhcdwK1PQiNQwVXGN8ccff6SwsJDbbrutSHvfvn15//33Adi3bx/ffffdddc/e/Ysjz/+OCNHjnRc66ks4eHhrFixgrNnzwJX33v166+/cvToURo2bMif/vQnoqOj2b17N3D1Gs2uXbsA+Pe//33dbTZp0oQzZ84AcOXKFY4cOcLAgQN58cUXOXXqlGNfNZVOnYnUMGZdY7x2jQauHjXExsY6Jru95vHHHycyMpJu3brRpUsXfHx8aNasmePzgQMHYhgGV65cYdSoUTz77LMVqul6hgwZwg8//ECfPn2Aqzce/POf/+TAgQNER0dTq1Yt6tSpw5tvvgnA888/z5QpU3j22WeL3Qhwzfjx43nkkUd4/fXXWbNmDVOmTOHUqVMYhsH06dNp3rx5pY/jVqJJNdGkmnLr++GHH+jatauzyyhVYWEhly9fpn79+hw8eJCwsDD2799P3bp1nV2alOJ6v2OaVFNE/nDOnz/PwIEDuXz5MoZhsHTpUoVMDaCgEZEq06RJE509qIF0M4CIiJhKQSMiIqZS0IiIiKkUNCIiYioFjYhUimPHjjF+/Hg8PT0JCgri7rvv5qeffnJ2WU5z8uRJli5d6lhOSUlh+PDhlb6f376uoKx+O1Hob82ZM4d//OMflVWag4JGRCrMMAxGjRpFaGgoBw8eZNeuXSxYsIDs7Owqr6WwsLD0TlXg90FTVn+U+iuTgkZEKmzz5s3UqVOHxx57zNEWEBBA//79MQyD6OhofH198fPzc0y/n5KSQmhoKKNHj6ZLly5MmDABwzBITExkzJgxju389kggKSmJPn360KNHD8aMGeOY2sXDw4OZM2fSo0cP1q1bx86dO/H39ycwMNCxb7j6JR4dHe14RcBbb711w1oAdu7cyZ133klAQAA9e/bkzJkzJW7nt2bNmuWYcTo6Ohq4OrXO9fbx+/oPHjzI0KFDCQoKon///vz4448ArFu3Dl9fXwICArjrrrsc+zp69ChDhw7F29ubGTNmONpXr16Nn58fvr6+zJw587p/d/Pnz6dTp07069eP/fvNeUuwnqMRqW4+mwXHrj+HWLm18YNhC0v8eO/evQQFBV33s7i4ONLT0/nmm284ceIEISEhji/JPXv28P3339O2bVv69u3LV199RVhYGFFRUZw7d45GjRqxdu1axo8fz4kTJ/jP//xPkpOTadSoES+++CIvv/wyzz33HAC33XabY34yX19fli9fTp8+fYq8cfOdd96hWbNm7Ny5k4sXL9K3b1+GDBlSYi09e/Zk3LhxrF27lpCQEE6fPk2DBg1K3E6HDh0c+1q4cCF79+51zDidkpJy3X3069evWP2DBw9m2bJleHt7s2PHDh5//HE2bdrEvHnz2LBhA1artcjrB9LT09mzZw/16tWjc+fOPPHEE9SuXZuZM2eya9cuWrRowZAhQ4iPjy8ym/SuXbtYs2YN6enpFBQU0KNHjxL/HitCQSMipvryyy954IEHqF27Nq1bt2bAgAHs3LmTpk2b0rNnT9zd3QEIDAzk8OHD9OvXj6FDh/Lxxx8zevRo1q9fz6JFi/jiiy/Yt28fffv2BeDSpUuO+coAxo0bB1w9ZXXmzBnHZw8++CCffPIJcPWI6Ntvv3W88fPUqVNkZGRQt27d69bSrFkz3NzcCAkJAa6+I+dG2/lt0FxPSeP9bf1nz57l66+/LnJUd/HiReDqpKSTJ09m7Nix3HfffY7PBw8e7Jgzrlu3bvzP//wPOTk5hIaG0rJlSwAmTJjAli1bigTN1q1bGTVqFA0bNgQgIiLihvWXl4JGpLq5wZGHWXx8fMr8uubfqlevnuPn2rVrO6bTHz9+PIsXL8bV1ZXg4GCaNGmCYRjYbDZWr1593W399lUAJTEMgzfeeIPw8PAi7SkpKSXWcjPbKc2N9nGt/itXrtC8efPrvntn2bJl7Nixg/Xr1xMUFOSYVfpmancGXaMRkQobNGgQFy9eJCYmxtH27bffsnXrVvr378/atWspLCzk+PHjbNmyhZ49e95wewMGDGD37t0sX77c8Srk3r1789VXX3HgwAHg6quXr3dXW/PmzWnSpAk7duwAYM2aNY7PwsPDefPNN7l8+TIAP/30E+fOnSuxjs6dO5OVlcXOnTuBq+/IKSgoKNN2fvvqgJtx7c2i69atA66G2jfffAPAwYMH6dWrF/PmzaNly5YcOXKkxO307NmTL774ghMnTlBYWMjq1asZMGBAkT533XUX8fHx5Ofnc+bMGT7++OObrrcsFDQiUmEWi4UPP/yQ5ORkPD098fHxYfbs2bRp04ZRo0bh7+9PQEAAgwYNYtGiRbRp0+aG26tduzbDhw/ns88+c9wI0LJlS1auXMkDDzyAv78/ffr0cVwk/7133nmHRx55hMDAQM6dO+c4rTR16lS6detGjx498PX15dFHH73hv/7r1q3L2rVreeKJJwgICMBms3HhwoUybee2226jb9+++Pr6Om4GKKv33nuPd955h4CAAHx8fEhISAAgOjracXH/2g0KJXFzc2PhwoUMHDiQgIAAgoKCGDFiRJE+PXr0YNy4cQQEBDBs2DDHKcLKptcEoNcEyK3vVnlNQFU5e/YsjRs3Bq5elM/KyuK1115zclW3Nr0mQETkN9avX8+CBQsoKCigffv2rFy50tkl1WgKGhGpdsaNG+e4i0ucT9doRKoJnQUXs1T0d0tBI1IN1K9fn5ycHIWNVDrDMMjJyaF+/frl3oZOnYlUA+7u7mRmZnL8+HFnlyLVUP369R0PmpaHgkakGqhTp06pT6WLOItTTp3l5uZis9nw9vbGZrORl5d33X6xsbF4e3vj7e1NbGyso33Xrl34+fnh5eXF9OnTHacL1q1bh4+PD7Vq1dLtyiIifxBOCZqFCxcyePBgMjIyGDx4MAsXFp8yIzc3l7lz57Jjxw5SU1OZO3euI5CmTZvG8uXLycjIICMjg8TERODqRHpxcXFFZjUVERHnckrQJCQkEBkZCUBkZCTx8fHF+mzYsAGbzYarqystWrTAZrORmJhIVlYWp0+fpnfv3lgsFiZNmuRYv2vXrnTu3LlKxyIiIjfmlGs02dnZuLm5AdCmTZvrvhzJbrfTrl07x7K7uzt2ux273V7kotS19psVExPjmJdJF1BFRMxjWtCEhYVx7NixYu3z588vsmyxWLBYLGaVUaKoqCiioqKAq9MoiIiIOUwLmuTk5BI/a926NVlZWbi5uZGVlUWrVq2K9bFaraSkpDiWMzMzCQ0NxWq1kpmZWaTdarVWau0iIlJ5nHKNJiIiwnEXWWxsbLEZReHqdN5JSUnk5eWRl5dHUlIS4eHhuLm50bRpU7Zv345hGKxateq664uIyB+DU4Jm1qxZfP7553h7e5OcnOx41WpaWhpTp04FwNXVlWeffZaQkBBCQkJ47rnncHV1BWDp0qVMnToVLy8vPD09GTZsGAAffvgh7u7ubNu2jXvuueemX0okIiKVT68JQK8JEBEpj7J+d2quMxERMZWCRkRETKWgERERUyloRETEVAoaERExlYJGRERMpaARERFTKWhERMRUChoRETGVgkZEREyloBEREVMpaERExFQKGhERMZWCRkRETKWgERERUyloRETEVAoaERExlYJGRERMpaARERFTKWhERMRUChoRETGVgkZEREyloBEREVMpaERExFQKGhERMZWCRkRETKWgERERUyloRETEVAoaERExlYJGRERMpaARERFTKWhERMRUChoRETGVy40+zMzMZM2aNWzdupWjR4/SoEEDfH19ueeeexg2bBi1aimnRETkxkpMioceeoiHH36YunXrMnPmTFavXs3SpUsJCwsjMTGRfv36sWXLlnLtNDc3F5vNhre3Nzabjby8vOv2i42NxdvbG29vb2JjYx3tu3btws/PDy8vL6ZPn45hGABER0fTpUsX/P39GTVqFCdPnixXfSIiUomMEnz33XclfWQYhmFcvHjRyMjIuGGfkkRHRxsLFiwwDMMwFixYYMyYMaNYn5ycHKNDhw5GTk6OkZuba3To0MHIzc01DMMwQkJCjG3bthlXrlwxhg4danz66aeGYRjGhg0bjMuXLxuGYRgzZsy47navJygoqFzjEBGpycr63VniEY2vr+8NA6pu3bp4eXmVK9wSEhKIjIwEIDIykvj4+GJ9NmzYgM1mw9XVlRYtWmCz2UhMTCQrK4vTp0/Tu3dvLBYLkyZNcqw/ZMgQXFyung3s3bs3mZmZ5apPREQqzw2v0QBkZGQwe/Zs9u3bx4ULFxztP//8c7l3mp2djZubGwBt2rQhOzu7WB+73U67du0cy+7u7tjtdux2O+7u7sXaf2/FihWMGzeuxBpiYmKIiYkB4Pjx4+Uei4iI3FipQfPQQw8xd+5cnnzySTZv3sy7777LlStXSt1wWFgYx44dK9Y+f/78IssWiwWLxXITJZdu/vz5uLi4MGHChBL7REVFERUVBUBwcHCl7l9ERP5fqUGTn5/P4MGDMQyD9u3bM2fOHIKCgpg3b94N10tOTi7xs9atW5OVlYWbmxtZWVm0atWqWB+r1UpKSopjOTMzk9DQUKxWa5FTYpmZmVitVsfyypUr+eSTT9i4cWOlB5iIiNy8Uu9PrlevHleuXMHb25vFixfz4Ycfcvbs2QrtNCIiwnEXWWxsLCNGjCjWJzw8nKSkJPLy8sjLyyMpKYnw8HDc3Nxo2rQp27dvxzAMVq1a5Vg/MTGRRYsW8dFHH9GwYcMK1SgiIpWktLsFUlNTjTNnzhhHjhwxJk+ebIwaNcrYtm1bhe5UOHHihDFo0CDDy8vLGDx4sJGTk2MYhmHs3LnTmDJliqPfO++8Y3h6ehqenp7GihUrHO07d+40fHx8jI4dOxp//vOfjStXrhiGYRienp6Gu7u7ERAQYAQEBBiPPvpomerRXWciIjevrN+dFsP4v4dQarDg4GDS0tKcXYaIyC2lrN+dJV6juffee294jeOjjz4qX2UiIlKjlBg0Tz/9NACGYfDII4/w9ttvV1lRIiJSfZQYNAMGDHD83Lhx4yLLIiIiZVWmWTF1m7CIiJRXiUc0ubm5jp8LCwvJy8vjt/cNuLq6mluZiIhUCyUGTVBQEBaLxREuPXr0cHxmsVgqNAWNiIjUHCUGzaFDh6qyDhERqaZKvEZz+PDhG65oGIZmRxYRkVKVeEQTHR3NlStXGDFiBEFBQbRs2ZILFy5w4MABNm/ezMaNG5k7d26RmZRFRER+r8SgWbduHfv27eO9995jxYoVZGVl0bBhQ7p27crdd9/NM888Q/369auyVhERuQXdcPbmbt26FZvWX0RE5GaU6TkaERGR8lLQiIiIqRQ0IiJiqlKD5r777mP9+vVlen2ziIjI75UaNI8//jj/+te/8Pb2ZtasWezfv78q6hIRkWqi1KAJCwvjvffeY/fu3Xh4eBAWFsadd97Ju+++y+XLl6uiRhERuYWV6RpNTk4OK1eu5O2336Z79+789a9/Zffu3dhsNrPrExGRW9wNn6MBGDVqFPv372fixIl8/PHHuLm5ATBu3DiCg4NNL1BERG5tpQbN9OnTGThw4HU/K8u7okVEpGYrNWjy8vKIi4sr0tasWTP8/Pxo1aqVaYWJiEj1UGrQvPPOO2zbts1xVJOSkkJQUBCHDh3iueeeY+LEiaYXKSIit65Sg+by5cv88MMPtG7dGoDs7GwmTZrEjh07uOuuuxQ0IiJyQ6XedZaZmekIGYBWrVpx5MgRXF1dqVOnjqnFiYjIra/UI5rQ0FCGDx/OmDFjAPj3v/9NaGgo586do3nz5qYXKCIit7ZSg2bJkiXExcXx5ZdfAjBp0iTuv/9+LBYLmzdvNr1AERG5td0waAoLCwkLC2Pz5s3cf//9VVWTiIhUIze8RlO7dm1q1arFqVOnqqoeERGpZko9dda4cWP8/Pyw2Ww0atTI0f7666+bWpiIiFQPpQbNfffdx3333VcVtYiISDVUatBERkaSn5/PL7/8QufOnauiJhERqUZKfY7m448/JjAwkKFDhwKQnp5ORESE6YWJiEj1UGrQzJkzh9TUVMczM4GBgfz888+mFyYiItVDqUFTp04dmjVrVnSlWmV6jY2IiEjpQePj48O//vUvCgsLycjI4IknnuDOO++s0E5zc3Ox2Wx4e3tjs9nIy8u7br/Y2Fi8vb3x9vYmNjbW0b5r1y78/Pzw8vJi+vTpGIYBwLPPPou/vz+BgYEMGTKEo0ePVqhOERGpuFKD5o033uD777+nXr16PPDAAzRt2pRXX321QjtduHAhgwcPJiMjg8GDB7Nw4cJifXJzc5k7dy47duwgNTWVuXPnOgJp2rRpLF++nIyMDDIyMkhMTAQgOjqab7/9lvT0dIYPH868efMqVKeIiFRcqUHTsGFD5s+fz86dO0lLS2P+/PnUr1+/QjtNSEggMjISuHpXW3x8fLE+GzZswGaz4erqSosWLbDZbCQmJpKVlcXp06fp3bs3FouFSZMmOdZv2rSpY/1z585hsVgqVKeIiFRcqbc3//TTT/zjH//g8OHDFBQUONo3bdpU7p1mZ2c7Xgndpk0bsrOzi/Wx2+20a9fOsezu7o7dbsdut+Pu7l6s/ZpnnnmGVatW0axZsxvOxRYTE0NMTAwAx48fL/dYRETkxkoNmjFjxvDYY48xdepUateuXeYNh4WFcezYsWLt8+fPL7JssVgq9chj/vz5zJ8/nwULFrB48WLmzp173X5RUVFERUUBEBwcXGn7FxGRokoNGhcXF6ZNm3bTG05OTi7xs9atW5OVlYWbmxtZWVnXfSW01WolJSXFsZyZmUloaChWq5XMzMwi7Vartdj6EyZM4O677y4xaEREpGqUeo3m3nvvZenSpWRlZZGbm+v4UxERERGOu8hiY2MZMWJEsT7h4eEkJSWRl5dHXl4eSUlJhIeH4+bmRtOmTdm+fTuGYbBq1SrH+hkZGY71ExIS6NKlS4XqFBGRirMY1+4NLkGHDh2Kr2SxVOihzZycHMaOHcsvv/xC+/btef/993F1dSUtLY1ly5bx9ttvA7BixQpeeOEF4Oq1l4ceegiAtLQ0Jk+eTH5+PsOGDeONN97AYrFw//33s3//fmrVqkX79u1ZtmzZdY92fi84OJi0tLRyj0dEpCYq63dnqUFTEyhoRERuXlm/O0s8dbZo0SLHz+vWrSvy2d/+9rcKlCYiIjVJiUGzZs0ax88LFiwo8tm1ByRFRERKU2LQ/PaM2u/Prulsm4iIlFWJQfPbZ1t+/5yLnrgXEZGyKvE5mm+++YamTZtiGAb5+fmO6V0Mw+DChQtVVqCIiNzaSgyawsLCqqxDRESqKb1YRkRETKWgERERUyloRETEVAoaERExlYJGRERMpaARERFTKWhERMRUChoRETGVgkZEREyloBEREVMpaERExFQKGhERMZWCRkRETKWgERERUyloRETEVAoaERExlYJGRERMpaARERFTKWhERMRUChoRETGVgkZEREyloBEREVMpaERExFQKGhERMZWCRkRETKWgERERUyloRETEVAoaERExlVOCJjc3F5vNhre3Nzabjby8vOv2i42NxdvbG29vb2JjYx3tu3btws/PDy8vL6ZPn45hGEXW+6//+i8sFgsnTpzcBOnXAAAMo0lEQVQwdRwiIlI6pwTNwoULGTx4MBkZGQwePJiFCxcW65Obm8vcuXPZsWMHqampzJ071xFI06ZNY/ny5WRkZJCRkUFiYqJjvSNHjpCUlMQdd9xRZeMREZGSOSVoEhISiIyMBCAyMpL4+PhifTZs2IDNZsPV1ZUWLVpgs9lITEwkKyuL06dP07t3bywWC5MmTSqy/pNPPsmiRYuwWCxVNh4RESmZizN2mp2djZubGwBt2rQhOzu7WB+73U67du0cy+7u7tjtdux2O+7u7sXa4WqAWa1WAgICSq0hJiaGmJgYAI4fP16h8YiISMlMC5qwsDCOHTtWrH3+/PlFli0WS6UcfZw/f54XXniBpKSkMvWPiooiKioKgODg4ArvX0RErs+0oElOTi7xs9atW5OVlYWbmxtZWVm0atWqWB+r1UpKSopjOTMzk9DQUKxWK5mZmUXarVYrBw8e5NChQ46jmczMTHr06EFqaipt2rSpvIGJiMhNccqps4iICGJjY5k1axaxsbGMGDGiWJ/w8HD+9re/OW4ASEpKYsGCBbi6utK0aVO2b99Or169WLVqFU888QR+fn78+uuvjvU9PDxIS0vj9ttvN28gn82CY9+Zt30REbO18YNhxW/IqkxOuRlg1qxZfP7553h7e5OcnMysWbMASEtLY+rUqQC4urry7LPPEhISQkhICM899xyurq4ALF26lKlTp+Ll5YWnpyfDhg1zxjBERKQMLMbvH0KpgYKDg0lLS3N2GSIit5SyfndqZgARETGVgkZEREyloBEREVM55a4zuTXF77Hz0ob9HD2ZT9vmDYgO78zI7lZnlyUif3AKGimT+D12Zsd9R/7lQgDsJ/OZHXf11m6FjYjciE6dSZm8tGG/I2Suyb9cyEsb9jupIhG5VShopEyOnsy/qXYRkWsUNFImbZs3uKl2EZFrFDRSJtHhnWlQp3aRtgZ1ahMd3tlJFYnIrUI3A0iZXLvgr7vORORmKWikzEZ2typYROSm6dSZiIiYSkEjIiKmUtCIiIipFDQiImIqBY2IiJhKQSMiIqZS0IiIiKkUNCIiYioFjYiImEpBIyIiplLQiIiIqRQ0IiJiKgWNiIiYSkEjIiKmUtCIiIip9D4akRuI32PXy95EKkhBI1KC+D12Zsd9R/7lQgDsJ/OZHfcdQLUOm5oYrhqzuWPWqTOREry0Yb8jZK7Jv1zISxv2O6ki810LV/vJfAz+P1zj99idXZppNGbzx6ygESnB0ZP5N9VeHdTEcNWYrzJzzAoakRK0bd7gptqrg5oYrhpz6e0VpaARKUF0eGca1KldpK1BndpEh3d2UkXmq4nhqjGX3l5RChqREozsbmXBfX5YmzfAAlibN2DBfX7V+iJxTQxXjfkqM8fslKDJzc3FZrPh7e2NzWYjLy/vuv1iY2Px9vbG29ub2NhYR/uuXbvw8/PDy8uL6dOnYxgGAHPmzMFqtRIYGEhgYCCffvpplYxHqq+R3a18NWsQhxbew1ezBlXrkIGaGa4as/ljthjXvqWr0IwZM3B1dWXWrFksXLiQvLw8XnzxxSJ9cnNzCQ4OJi0tDYvFQlBQELt27aJFixb07NmT119/nV69enH33Xczffp0hg0bxpw5c2jcuDFPP/30TdVzbT8iIlJ2Zf3udMoRTUJCApGRkQBERkYSHx9frM+GDRuw2Wy4urrSokULbDYbiYmJZGVlcfr0aXr37o3FYmHSpEnXXV9ERP4YnBI02dnZuLm5AdCmTRuys7OL9bHb7bRr186x7O7ujt1ux2634+7uXqz9msWLF+Pv78/DDz9c4ik5ERGpOqYFTVhYGL6+vsX+JCQkFOlnsViwWCyVss9p06Zx8OBB0tPTcXNz46mnniqxb0xMDMHBwQQHB3P8+PFK2b+IiBRn2hQ0ycnJJX7WunVrsrKycHNzIysri1atWhXrY7VaSUlJcSxnZmYSGhqK1WolMzOzSLvVanVs95pHHnmE4cOHl1hDVFQUUVFRwNXzjCIiYg6nnDqLiIhw3EUWGxvLiBEjivUJDw8nKSmJvLw88vLySEpKIjw8HDc3N5o2bcr27dsxDINVq1Y51s/KynKs/+GHH+Lr61s1AxIRkRI55a6znJwcxo4dyy+//EL79u15//33cXV1JS0tjWXLlvH2228DsGLFCl544QUAnnnmGR566CEA0tLSmDx5Mvn5+QwbNow33ngDi8XCxIkTSU9Px2Kx4OHhwVtvveW4FnQjt99+Ox4eHuUay/Hjx2nZsmW51r1Vacw1g8ZcM1RkzIcPH+bEiROl9nNK0FQnNfHWaI25ZtCYa4aqGLNmBhAREVMpaERExFS158yZM8fZRdzqgoKCnF1CldOYawaNuWYwe8y6RiMiIqbSqTMRETGVgkZEREyloKmAwsJCunfvfsMZCKobDw8P/Pz8CAwMrBEzKpw8eZLRo0fTpUsXunbtyrZt25xdkqn279/veM1GYGAgTZs25dVXX3V2WaZ75ZVX8PHxwdfXlwceeIALFy44uyTTvfbaa/j6+uLj42P637Gu0VTAyy+/TFpaGqdPn+aTTz5xdjlVwsPDg7S0NG6//XZnl1IlIiMj6d+/P1OnTuXSpUucP3+e5s2bO7usKlFYWIjVamXHjh20b9/e2eWYxm63069fP/bt20eDBg0YO3Ysd999N5MnT3Z2aabZu3cv48ePJzU1lbp16zJ06FCWLVuGl5eXKfvTEU05ZWZmsn79eqZOnersUsQkp06dYsuWLUyZMgWAunXr1piQAdi4cSOenp7VOmSuKSgoID8/n4KCAs6fP0/btm2dXZKpfvjhB3r16kXDhg1xcXFhwIABxMXFmbY/BU05/cd//AeLFi2iVq2a9Z/QYrEwZMgQgoKCiImJcXY5pjp06BAtW7bkoYceonv37kydOpVz5845u6wqs2bNGh544AFnl2E6q9XK008/zR133IGbmxvNmjVjyJAhzi7LVL6+vmzdupWcnBzOnz/Pp59+ypEjR0zbX836lqwkn3zyCa1ataqR99t/+eWX7N69m88++4wlS5awZcsWZ5dkmoKCAnbv3s20adPYs2cPjRo1YuHChc4uq0pcunSJjz76iDFjxji7FNPl5eWRkJDAoUOHOHr0KOfOneOf//yns8syVdeuXZk5cyZDhgxh6NChBAYGUrt2bdP2p6Aph6+++oqPPvoIDw8Pxo8fz6ZNm/jTn/7k7LKqxLVXMrRq1YpRo0aRmprq5IrM4+7ujru7O7169QJg9OjR7N6928lVVY3PPvuMHj16FHn1RnWVnJxMhw4daNmyJXXq1OG+++7j66+/dnZZppsyZQq7du1iy5YttGjRgk6dOpm2LwVNOSxYsIDMzEwOHz7MmjVrGDRoULX/FxDAuXPnOHPmjOPnpKSkav0qhjZt2tCuXTv2798PXL1m0a1bNydXVTVWr15dI06bAdxxxx1s376d8+fPYxgGGzdupGvXrs4uy3S//vorAL/88gtxcXE8+OCDpu3LtBefSfWTnZ3NqFGjgKunlR588EGGDh3q5KrM9cYbbzBhwgQuXbpEx44deffdd51dkunOnTvH559/zltvveXsUqpEr169GD16ND169MDFxYXu3bs7XopYnd1///3k5ORQp04dlixZYuqNLrq9WURETKVTZyIiYioFjYiImEpBIyIiplLQiIiIqRQ0IiJiKgWNSBV4/fXX6dq1KxMmTKjU7aakpDhmD09JSakRDxrKrUfP0YhUgaVLl5KcnIy7u7ujraCgABeXyvtfMCUlhcaNG3PnnXdW2jZFKoOOaERM9thjj/Hzzz8zbNgwmjVrxsSJE+nbty8TJ05k5cqVjBw5EpvNhoeHB4sXL+bll1+me/fu9O7dm9zcXABCQ0NJS0sD4MSJE3h4eBTZx+HDh1m2bBmvvPIKgYGBbN26taqHKVIiBY2IyZYtW0bbtm3ZvHkzTz75JPv27SM5OZnVq1cDV98NEhcXx86dO3nmmWdo2LAhe/bsoU+fPqxatapM+/Dw8OCxxx7jySefJD09nf79+5s5JJGboqARqWIRERE0aNDAsTxw4ECaNGlCy5YtadasGffeey8Afn5+HD582ElVilQeBY1IFWvUqFGR5Xr16jl+rlWrlmO5Vq1aFBQUAODi4sKVK1cAasRrhqV6UdCI3AI8PDzYtWsXAB988MF1+zRp0sQxu7bIH4mCRuQW8PTTT/Pmm2/SvXt3Tpw4cd0+9957Lx9++KFuBpA/HM3eLCIiptIRjYiImEpBIyIiplLQiIiIqRQ0IiJiKgWNiIiYSkEjIiKmUtCIiIip/hfTctv3XiSGKQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "energies = [job.logfile.energy\n", " for job in rmc.queue if job.is_completed]\n", "rmults = [job.param\n", " for job in rmc.queue if job.is_completed]\n", "frmult = [rmult[0] for rmult in rmults]\n", "n_at = len(rmc.queue[0].posinp)\n", "threshold = min(energies) + n_at*rmc.precision_per_atom\n", "\n", "fig=plt.figure()\n", "fig.patch.set_facecolor('white') # When dark background\n", "plt.plot(frmult, energies,\n", " label=\"BigDFT results\", marker=\"o\", linestyle='')\n", "plt.plot(frmult, [threshold]*len(energies),\n", " label=\"Convergence threshold\")\n", "plt.xlabel(\"frmult\")\n", "plt.ylabel(\"Energy (Ha)\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remember that there are actually two grids: one known as the fine grid, the other being the coarse one. The former is more expensive (from a computational point of view) than the latter, and should therefore have a smaller extension. One can perform a more refined rmult convergence by testing the convergence with respect to that fine grid extensions: this will be done after the convergence with respect to `hgrids`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## hgrids convergence\n", "\n", "The `HgridsConvergence` class allows to run all the necessary calculations to determine the largest `hgrids` which must be used so that the energy error compared to the reference calculation (with the lowest `hgrids` considered) lie within the required precision per atom. Using the largest `hgrids` as possible allows to save computational time because there are less grid points (or degrees of freedom)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "from mybigdft.workflows.convergences import HgridsConvergence\n", "\n", "atoms = [Atom('N', [0, 0, 0]), Atom('N', [0, 0, 1.1])]\n", "pos = Posinp(atoms, units=\"angstroem\", boundary_conditions=\"free\")\n", "inp = InputParams({\"dft\": {\"rmult\": rmc.converged.param}})\n", "run_dir = \"N2/hgrids_convergence/rm_{:.1f}_{:.1f}\".format(*rmc.converged.param)\n", "base = Job(inputparams=inp, posinp=pos, run_dir=run_dir, name=\"N2\")\n", "hgc = HgridsConvergence(base, reference=0.24, delta=0.02,\n", " precision_per_atom=0.01*EV_TO_HA)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The jobs in the queue have an extra attribute `hgrids`:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.24, 0.24, 0.24]\n", "[0.26, 0.26, 0.26]\n", "[0.28, 0.28, 0.28]\n", "[0.3, 0.3, 0.3]\n", "[0.32, 0.32, 0.32]\n", "[0.34, 0.34, 0.34]\n", "[0.36, 0.36, 0.36]\n", "[0.38, 0.38, 0.38]\n", "[0.4, 0.4, 0.4]\n", "[0.42, 0.42, 0.42]\n" ] } ], "source": [ "for job in hgc.queue:\n", " print(job.param)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can run the calculations as usual:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/0.24_0.24_0.24\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/0.26_0.26_0.26\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/0.28_0.28_0.28\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/0.30_0.30_0.30\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/0.32_0.32_0.32\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/0.34_0.34_0.34\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/0.36_0.36_0.36\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/0.38_0.38_0.38\n", "Logfile log-N2.yaml already exists!\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/workflows/workflow.py:114: UserWarning: Some jobs of the workflow were not run.\n", " UserWarning)\n" ] } ], "source": [ "hgc.run(nmpi=6, nomp=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that not all calculations were run: this is because the job with `hgrids` = 0.38 gave an energy that was above the requested precision per atom.\n", "\n", "Once the calculations are performed, a converged job is determined. It corresponds to the one with the largest `hgrids` so that the total energy of the system is below the convergence threshold. It can be accessed via the `converged` attribute. Its `param` attribute can in turns be accessed:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.36, 0.36, 0.36]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hgc.converged.param" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result of the previous cell shows that `hgrids = 0.36` gives converged results for the given precision per atom. To make sure that everything ran as expected, all the relevant data can be accessed as follows:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requested precision per atom: 3.67e-04 (Ha)\n", "-----------------------------------------------------------\n", " hgrids precision_per_atom (Ha) is_converged \n", "-----------------------------------------------------------\n", " [0.24, 0.24, 0.24] 0.00e+00 True \n", " [0.26, 0.26, 0.26] 3.89e-07 True \n", " [0.28, 0.28, 0.28] 2.63e-05 True \n", " [0.30, 0.30, 0.30] 4.53e-05 True \n", " [0.32, 0.32, 0.32] 1.23e-04 True \n", " [0.34, 0.34, 0.34] 1.88e-04 True \n", " [0.36, 0.36, 0.36] 3.46e-04 True \n", " [0.38, 0.38, 0.38] 5.86e-04 False \n" ] } ], "source": [ "hgc.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It shows that everything ran as expected: the calculation with `hgrids = 0.36` has the largest `hgrids` while still having an energy below the requested threshold." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the convergence\n", "\n", "The figure below shows the convergence of the total energy of the N$_2$ system with respect to `hgrids`. The threshold value is also shown for comparison." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEVCAYAAACmMTGfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XlcVPX+P/DXiBoqoWKiMLigjAQMiwKGmoriOOBVVFyyWxcsjdar9z4S1Ns1s/KKWrlUaqgp2qqlUOkFLiGmuSAaZlmKCynjSArI5mAMfH5/8PX8pGEVhsPyej4ePB4zn/mcc16HcN6dz/mccxRCCAEiIqIm1k7uAERE1DaxABERkSxYgIiISBYsQEREJAsWICIikgULEBERyYIFqAF2794NNzc3tGvXDmlpadX2W7duHdRqNdzc3LB27VqpPT09HX5+fvDy8oKPjw9SU1MBAPn5+Zg0aRI8PT3h5uaGbdu21Zrl6aefhq2tLdRqdcN3jIioCbAANYBarcaePXswatSoavv89NNP2Lx5M1JTU3H69Gl88803uHDhAgAgMjISS5cuRXp6Ol5//XVERkYCAN5//324urri9OnTSElJwcsvv4w//vijxiyzZ89GfHx84+0cEZGZsQA1gIuLC5ydnWvs88svv+CRRx5B586d0b59e4wePRp79uwBACgUChQUFACoOOqxt7eX2gsLCyGEQFFREWxsbNC+fXsAwOrVq+Hr6wsPDw8sXbpU2s6oUaNgY2Njjt0kIjKL9nIHaO3UajVeeeUV5OTkoFOnTti/fz98fHwAAGvXroVWq8WCBQtQXl6OI0eOAABeeuklBAcHw97eHoWFhfj888/Rrl07JCYmIiMjA6mpqRBCIDg4GN99912NR2BERM0VC1Atxo0bh+vXr5u0L1++HJMnT651eRcXFyxcuBDjx49Hly5d4OXlBQsLCwDAxo0bsWbNGkybNg27du3CnDlzkJSUhISEBHh5eSE5ORkXL16ERqPByJEjkZiYiMTERAwePBgAUFRUhIyMDBYgImqRWIBqkZSU1OB1zJkzB3PmzAEA/Otf/4KDgwMAICYmBuvWrQMAzJgxA3PnzgUAbNu2DYsWLYJCoYCTkxMcHR3x66+/QgiBxYsX49lnn21wJiIiufEcUBP4/fffAQBXrlzBnj178Ne//hUAYG9vj4MHDwIAkpOToVKpAAB9+/bFt99+CwDIzs7GuXPnMGDAAGi1Wnz44YcoKioCAOh0OmndREQtjqD7tmfPHqFUKkXHjh2Fra2tGD9+vBBCCJ1OJ4KCgqR+jz76qHBxcREeHh4iKSlJaj906JAYMmSI8PDwEEOHDhVpaWnS8hqNRqjVauHm5iZ27twpLbN27VqhVquFWq0Wfn5+4sKFC0IIIWbNmiV69+4t2rdvL5RKpdiyZUtT/AqIiO6bQgg+joGIiJoeh+CIiEgWnIRQg4ceegj9+/eXOwYRUYuSmZmJmzdv1tqPBagG/fv3r/EWO0REZOrutY614RAcERHJggWIiIhkwQJERESy4DmgeiotLUVWVhZKSkrkjkKtkKWlJRwcHNChQwe5oxCZHQtQPWVlZeHBBx9E//79oVAo5I5DrYgQAjk5OcjKyoKjo6PccYjMjgWonkpKSlh8yCwUCgV69OiBGzduyB2F2rDYH3RYnXAO124ZYN+tEyK0zpgyWGmWbbEA3QcWHzIX/m2RnGJ/0GHxnjMwlJYBAHS3DFi85wwAmKUIcRICEREBAFYnnJOKz12G0jKsTjhnlu2xALVAFhYW8PLygqenJ4YMGSI9yO7atWuYPn16rcv3798f7u7ucHd3h6urK/79739LkyoyMzPRqVMneHl5ST8ffPCB9Lpjx45wd3eHl5cXFi1aZNb9vGv79u146aWXAACxsbE4e/Zsk2yXqK25dstQr/aG4hCcmZljPLVTp05IT08HACQkJGDx4sU4ePAg7O3t8cUXX9RpHQcOHMBDDz2EoqIihIeH49lnn0VMTAwAYODAgdL677r7DKL+/ftLy9bEaDRKjxFvTLGxsZg4cSJcXV0bfd1EbZ19t07QVVFs7Lt1Msv2eARkRnfHU3W3DBD4/+OpsT/oGm0bBQUF6N69O4CKoxe1Wg0AuH37NmbOnAlXV1dMnToVjzzySJW3FbKyssKmTZsQGxuL3NzcBmXZvn07goODMXbsWAQEBAAAVq9eDV9fX3h4eGDp0qUAgOLiYvzlL3+Bp6cn1Go1Pv/8cwAVxe3u/aPS0tLg7+9faf1HjhzBV199hYiICHh5eeHixYtYv349XF1d4eHhgVmzZjUoP1FbF6F1RqcOFpXaOnWwQITW2Szb4xGQGdU0ntqQoyCDwQAvLy+UlJRAr9cjOTnZpM+GDRvQvXt3nD17Fj/99BO8vLyqXZ+1tTUcHR2RkZGBXr164eLFi1L/ESNG4P33369ztlOnTuHHH3+EjY0NEhMTkZGRgdTUVAghEBwcjO+++w43btyAvb099u3bBwDIz8+v07qHDx+O4OBgTJw4URpqjIqKwuXLl/HAAw/g1q1bdc5JRKbufi9xFlwrYK7x1HuH4I4ePYrQ0FD89NNPlfocPnwY8+fPBwCo1Wp4eHjUuM57HwtV1RBcXWk0GtjY2AAAEhMTkZiYiMGDBwMAioqKkJGRgZEjR+Lll1/GwoULMXHiRIwcOfK+tgUAHh4eeOKJJzBlyhRMmTLlvtdDRBWmDFaareD8GYfgzKi6cdPGHE8dNmwYbt682aBrRwoLC5GZmYlBgwY1OE+XLl2k10IILF68GOnp6UhPT8eFCxcwZ84cDBo0CKdOnYK7uzv+/e9/4/XXXwcAtG/fHuXl5QBQ5ztN7Nu3Dy+++CJOnToFX19fGI3GBu8DETUNFiAzaorx1F9//RVlZWXo0aNHpfYRI0Zg165dAICzZ8/izJkzVS5fVFSEF154AVOmTJHOJTUWrVaLDz/8EEVFRQAAnU6H33//HdeuXUPnzp3x5JNPIiIiAqdOnQJQcQ7o5MmTAIAvv/yyynU++OCDKCwsBACUl5fj6tWrGDNmDFauXIn8/HxpW0TU/HEIzozMNZ569xwQUHGUERMTAwuLyoXuhRdeQFhYGFxdXfHwww/Dzc0NXbt2lT4fM2YMhBAoLy/H1KlTsWTJkgZlqsr48ePxyy+/YNiwYQAqJjx89NFHuHDhAiIiItCuXTt06NABGzduBAAsXboUc+bMwZIlS0wmINw1a9YsPPPMM1i/fj0+++wzzJkzB/n5+RBCYN68eejWrVuj7wcRmYdC3Dv4T5X4+PiYzBz75Zdf4OLiIlOiuisrK0NpaSksLS1x8eJFjBs3DufOnUPHjh3ljka1aCl/Y0TVqeq7syo8Amqlbt++jTFjxqC0tBRCCGzYsIHFh4iaFRagVurBBx/k48SJqFmTZRJCbm4uNBoNVCoVNBoN8vLyquwXExMDlUoFlUolXaUPACdPnoS7uzucnJwwb948aQrxY489Jt0ypn///tJ5ktLSUoSFhcHd3R0uLi5YsWKF+XeSiIhqJEsBioqKQkBAADIyMhAQEICoqCiTPrm5uVi2bBmOHz+O1NRULFu2TCpUzz//PDZv3oyMjAxkZGQgPj4eAPD5559LU36nTZuGkJAQAMDu3btx584dnDlzBidPnsQHH3yAzMzMJttfIiIyJUsBiouLQ1hYGAAgLCwMsbGxJn0SEhKkixq7d+8OjUaD+Ph46PV6FBQUwM/PDwqFAqGhoSbLCyGwa9cuPP744wAqbnFfXFwMo9EIg8GAjh07wtra2vw7SkRE1ZKlAGVnZ8POzg4A0Lt3b2RnZ5v00el06NOnj/TewcEBOp0OOp0ODg4OJu33OnToEHr16gWVSgUAmD59Orp06QI7Ozv07dsXCxYskK7WJyIieZitAI0bNw5qtdrkJy4urlI/hULR6A/h+vTTT6WjHwBITU2FhYUFrl27hsuXL+Ptt9/GpUuXqlw2OjoaPj4+8PHxabZPprx+/TpmzZqFgQMHwtvbGxMmTMD58+fljiWbW7duYcOGDdL7lJQUTJw4sdG3c+9jIerq3hus3uu1117DW2+91VjRiFoks82CS0pKqvazXr16Qa/Xw87ODnq9Hra2tiZ9lEolUlJSpPdZWVnw9/eHUqlEVlZWpXal8v9f2Gk0GrFnzx7pinoA+OSTTxAYGIgOHTrA1tYWI0aMQFpaGgYMGGCy3fDwcISHhwOomMve3AghMHXqVISFheGzzz4DAJw+fRrZ2dmNciud+igrKzO5AFYOdwvQCy+8UK/lmkt+orZKliG44OBgaVZbTEwMJk+ebNJHq9UiMTEReXl5yMvLQ2JiIrRaLezs7GBtbY1jx45BCIEdO3ZUWj4pKQkPP/xwpWG6vn37SneMLi4uxrFjx/Dwww+beS/N48CBA+jQoQOee+45qc3T0xMjR46EEAIRERFQq9Vwd3eXHnOQkpICf39/TJ8+HQ8//DCeeOIJCCEQHx+PGTNmSOu598ghMTERw4YNw5AhQzBjxgzpFjf9+/fHwoULMWTIEOzevRsnTpyAh4cHvLy8pG0DFV/uERER0qMYPvjggxqzAMCJEycwfPhweHp6YujQoSgsLKx2PfdatGiRdAfviIgIABW3GKpqG3/Of/HiRQQGBsLb2xsjR47Er7/+CqBi4oparYanpydGjRolbevatWsIDAyESqVCZGSk1P7pp5/C3d0darUaCxcurPK/3fLlyzFo0CA8+uijOHfOPE+YJGpRhAxu3rwpxo4dK5ycnERAQIDIyckRQghx4sQJMWfOHKnf1q1bxcCBA8XAgQPFhx9+KLWfOHFCuLm5iQEDBogXX3xRlJeXS5+FhYWJjRs3VtpeYWGhmD59unB1dRUuLi5i1apVdcrp7e1t0nb27Nn//2b/QiE+nNC4P/sX1php3bp14h//+EeVn33xxRdi3Lhxwmg0iuvXr4s+ffqIa9euiQMHDghra2tx9epVUVZWJvz8/MShQ4dEaWmp6NOnjygqKhJCCPHcc8+JnTt3ihs3boiRI0dK7VFRUWLZsmVCCCH69esnVq5cKW3Tzc1NHDlyRAghxMKFC4Wbm5sQQogPPvhAvPHGG0IIIUpKSoS3t7e4dOlStVnu3LkjHB0dRWpqqhBCiPz8fFFaWlrteu51+fJlabtCiGq3UVX+sWPHivPnzwshhDh27JgYM2aMEEIItVotsrKyhBBC5OXlCSGE2LZtm3B0dBS3bt0SBoNB9O3bV1y5ckXodDrRp08f8fvvv4vS0lIxZswYsXfvXml7N27cEGlpaUKtVovi4mKRn58vBg4cKFavXl3lf8dKf2NELVBV351VkeVC1B49euDbb781affx8cGWLVuk908//TSefvrpKvv9+fEDd23fvt2kzcrKCrt3777/wC3E4cOH8fjjj8PCwgK9evXC6NGjceLECVhbW2Po0KHSUaGXlxcyMzPx6KOPIjAwEF9//TWmT5+Offv2YdWqVTh48CDOnj2LESNGAAD++OMP6X5uQMX1VkDF0FdhYaH02V//+ld88803ACqOoH788UfpCa35+fnIyMhAx44dq8zStWtX2NnZwdfXFwCkWYrVrcfR0bHG30V1+3tv/qKiIhw5cqTSUeCdO3cAVNzMdfbs2Zg5c6Y0nR8AAgICpHvqubq64rfffkNOTg78/f3Rs2dPAMATTzyB7777rtLjIQ4dOoSpU6eic+fOACpGAYjaOt4JoSGCTK9fMjc3N7c6P3b7Xg888ID02sLCQnpswaxZs/Dee+/BxsYGPj4+ePDBByGEgEajwaefflrluu595EJ1hBB49913odVqK7WnpKRUm6U+66lNTdu4m7+8vBzdunWr8tlHmzZtwvHjx7Fv3z54e3tL5xTrk52IasbHMbQwY8eOxZ07dxAdHS21/fjjjzh06BBGjhyJzz//HGVlZbhx4wa+++47DB06tMb1jR49GqdOncLmzZulR1r7+fnh+++/x4ULFwBUnDerapZdt27d8OCDD+L48eMAIE2KACrO4W3cuBGlpaUAgPPnz6O4uLjaHM7OztDr9Thx4gSAimcUGY3GOq3n3kc01MfdJ8HePToWQuD06dMAgIsXL+KRRx7B66+/jp49e+Lq1avVrmfo0KE4ePAgbt68ibKyMnz66acYPXp0pT6jRo1CbGwsDAYDCgsL8fXXX9c7L1FrwwLUwigUCuzduxdJSUkYOHAg3NzcsHjxYvTu3RtTp06Fh4cHPD09MXbsWKxatQq9e/eucX0WFhaYOHEi/vvf/0oTEHr27Int27fj8ccfh4eHB4YNGyadnP+zrVu34plnnoGXlxeKi4ul4am5c+fC1dUVQ4YMgVqtxrPPPlvj0ULHjh3x+eef4+9//zs8PT2h0WhQUlJSp/X06NEDI0aMgFqtliYh1NXHH3+MrVu3wtPTE25ubtJlAhEREdKkgrsTI6pjZ2eHqKgojBkzBp6envD29jaZWDNkyBA89thj8PT0RFBQkDTUSNSW8XEMNWjJj2NoKkVFRbCysgJQcYslvV6PdevWyZyqZePfGLV0fBwDNYl9+/ZhxYoVMBqN6NevX5WTQIiIqsICRA3y2GOPSbPKiIjqg+eA7gNHLclc+LdFbQkLUD1ZWloiJyeHXxTU6IQQyMnJgaWlpdxRiJoEh+DqycHBAVlZWc32RqXUsllaWla6jRRRa8YCVE8dOnSo9Sp8IiKqHYfgiIhIFixAREQkCxYgIiKSBQsQERHJggWIiIhkwQJERESyYAEiIiJZsAAREZEsWICIiEgWLEBERCQLFiAiIpIFCxAREcmCBYiIiGTBAkRERLJgASIiIlnIUoByc3Oh0WigUqmg0WiQl5dXZb+YmBioVCqoVCrExMRI7SdPnoS7uzucnJwwb9486emk6enp8PPzg5eXF3x8fJCamgqg4kmT8+bNg5OTEzw8PHDq1Cnz7yQREdVMyCAiIkKsWLFCCCHEihUrRGRkpEmfnJwc4ejoKHJyckRubq5wdHQUubm5QgghfH19xdGjR0V5ebkIDAwU+/fvF0IIodFopNf79u0To0ePll4HBgaK8vJycfToUTF06NA65fT29m7orhIRtTl1/e6U5QgoLi4OYWFhAICwsDDExsaa9ElISIBGo4GNjQ26d+8OjUaD+Ph46PV6FBQUwM/PDwqFAqGhodLyCoUCBQUFAID8/HzY29tL2wsNDYVCoYCfnx9u3boFvV7fRHtLRERVkeWR3NnZ2bCzswMA9O7dG9nZ2SZ9dDod+vTpI713cHCATqeDTqeDg4ODSTsArF27FlqtFgsWLEB5eTmOHDlS47ruZrhXdHQ0oqOjAQA3btxohL0lIqKqmK0AjRs3DtevXzdpX758eaX3CoUCCoWiUba5ceNGrFmzBtOmTcOuXbswZ84cJCUl1Wsd4eHhCA8PBwD4+Pg0Si4iIjJltgJU0xd/r169oNfrYWdnB71eD1tbW5M+SqUSKSkp0vusrCz4+/tDqVQiKyurUrtSqQRQMWlh3bp1AIAZM2Zg7ty50rquXr1a5TJERCQPWc4BBQcHS7PaYmJiMHnyZJM+Wq0WiYmJyMvLQ15eHhITE6HVamFnZwdra2scO3YMQgjs2LFDWt7e3h4HDx4EACQnJ0OlUknb27FjB4QQOHbsGLp27Vrl8BsRETUdWc4BLVq0CDNnzsTWrVvRr18/7Nq1CwCQlpaGTZs2YcuWLbCxscGSJUvg6+sLAHj11VdhY2MDANiwYQNmz54Ng8GAoKAgBAUFAQA2b96M+fPnw2g0wtLSUjqXM2HCBOzfvx9OTk7o3Lkztm3bJsNeExHRvRRC/N9FNGTCx8cHaWlpcscgImpR6vrdyTshEBGRLFiAiIhIFixAREQkCxYgIiKSBQsQERHJggWIiIhkwQJERESyYAEiIiJZsAAREZEsWICIiEgWLEBERCQLFiAiIpIFCxAREcmCBYiIiGTBAkRERLJgASIiIlmwABERkSxYgIiISBbta/owKysLn332GQ4dOoRr166hU6dOUKvV+Mtf/oKgoCC0a8f6RURE96faAvTUU09Bp9Nh4sSJWLhwIWxtbVFSUoLz588jPj4ey5cvR1RUFEaNGtWUeVuG/y4Crp+ROwUR0f3r7Q4ERZl1E9UWoJdffhlqtdqkXa1WIyQkBH/88QeuXLli1nBERNR6VVuAqio+9+rYsSOcnJwaPVCrYOb/ayAiag1qPAcEABkZGVi8eDHOnj2LkpISqf3SpUtmDUZERK1brbMInnrqKTz//PNo3749Dhw4gNDQUDz55JNNkY2IiFqxWguQwWBAQEAAhBDo168fXnvtNezbt68pshERtXixP+gwIioZjov2YURUMmJ/0MkdqdmotQA98MADKC8vh0qlwnvvvYe9e/eiqKioQRvNzc2FRqOBSqWCRqNBXl5elf1iYmKgUqmgUqkQExMjtZ88eRLu7u5wcnLCvHnzIIQAAKSnp8PPzw9eXl7w8fFBamoqAODjjz+Gh4cH3N3dMXz4cJw+fbpB+YmI6iL2Bx0W7zkD3S0DBADdLQMW7znDIvR/ai1A69atw+3bt7F+/XqcPHkSO3furFQM7kdUVBQCAgKQkZGBgIAAREWZnrTPzc3FsmXLcPz4caSmpmLZsmVSoXr++eexefNmZGRkICMjA/Hx8QCAyMhILF26FOnp6Xj99dcRGRkJAHB0dMTBgwdx5swZLFmyBOHh4Q3KT0RUF6sTzsFQWlapzVBahtUJ52RK1LzUOgnB19cXAGBlZYVt27Y1ykbj4uKQkpICAAgLC4O/vz9WrlxZqU9CQgI0Gg1sbGwAABqNBvHx8fD390dBQQH8/PwAAKGhoYiNjUVQUBAUCgUKCgoAAPn5+bC3twcADB8+XFqvn58fsrKyGmU/iIhqcu2WoV7tbU21BWjSpElQKBTVLvjVV1/d90azs7NhZ2cHAOjduzeys7NN+uh0OvTp00d67+DgAJ1OB51OBwcHB5N2AFi7di20Wi0WLFiA8vJyHDlyxGS9W7duRVBQULXZoqOjER0dDQC4cePG/e0gEREA+26doKui2Nh36yRDmuan2gK0YMECAIAQAs888wy2bNlSrxWPGzcO169fN2lfvnx5pfcKhaLGQlcfGzduxJo1azBt2jTs2rULc+bMQVJSkvT5gQMHsHXrVhw+fLjadYSHh0tDdD4+Po2Si4japgitMxbvOVNpGK5TBwtEaJ1lTNV8VFuARo8eLb22srKq9L4u7v3i/7NevXpBr9fDzs4Oer0etra2Jn2USqU0TAdU3JfO398fSqWy0hBaVlYWlEolgIpJC+vWrQMAzJgxA3PnzpX6/fjjj5g7dy7++9//okePHvXaFyKi+zFlcMV30+qEc7h2ywD7bp0QoXWW2tu6Ot1NtLGOUO4KDg6WJjLExMRg8uTJJn20Wi0SExORl5eHvLw8JCYmQqvVws7ODtbW1jh27BiEENixY4e0vL29PQ4ePAgASE5OhkqlAgBcuXIFISEh2LlzJwYNGtSo+0JEVJMpg5X4ftFYXI76C75fNJbF5x7VHgHl5uZKr8vKypCXlydNdwYgTQ64H4sWLcLMmTOxdetW9OvXD7t27QIApKWlYdOmTdiyZQtsbGywZMkSaRLEq6++Km1zw4YNmD17NgwGA4KCgqRzOps3b8b8+fNhNBphaWkpnct5/fXXkZOTgxdeeKFip9u3R1pa2n3nJyKihlOIe6vKPRwdHaFQKFDVxwqFok3cisfHx4eFioionur63VntEdDly5cbNRAREdG9qj0HlJmZWeOCQgheT0NERPet2iOgiIgIlJeXY/LkyfD29kbPnj1RUlKCCxcu4MCBA/j222+xbNmyStfkEBER1VW1BWj37t04e/YsPv74Y3z44YfQ6/Xo3LkzXFxcMGHCBLzyyiuwtLRsyqxERNSK1HgrHldXV5MLR4mIiBpDna4DIiIiamwsQEREJAsWICIikkWtBSgkJAT79u1DeXl5U+QhIqI2otYC9MILL+CTTz6BSqXCokWLcO4cH6REREQNV2sBGjduHD7++GOcOnUK/fv3x7hx4zB8+HBs27YNpaWlTZGRiIhaoTqdA8rJycH27duxZcsWDB48GPPnz8epU6eg0WjMnY+IiFqpWh/JPXXqVJw7dw5/+9vf8PXXX0tPMn3sscf4wDYiIrpvtRagefPmYcyYMVV+xjtFExHR/aq1AOXl5WHPnj2V2rp27Qp3d/cqn2RKRERUF7UWoK1bt+Lo0aPSUVBKSgq8vb1x+fJlvPrqq/jb3/5m9pBERNT61FqASktL8csvv6BXr14AgOzsbISGhuL48eMYNWoUCxAREd2XWmfBZWVlScUHAGxtbXH16lXY2NigQ4cOZg1HREStV61HQP7+/pg4cSJmzJgBAPjyyy/h7++P4uJidOvWzewBiYiodaq1AL3//vvYs2cPDh8+DAAIDQ3FtGnToFAocODAAbMHJCKi1qnGAlRWVoZx48bhwIEDmDZtWlNlIiKiNqDGc0AWFhZo164d8vPzmyoPERG1EbUOwVlZWcHd3R0ajQZdunSR2tevX2/WYERE1LrVWoBCQkIQEhLSFFmIiKgNqbUAhYWFwWAw4MqVK3B2dm6KTERE1AbUeh3Q119/DS8vLwQGBgIA0tPTERwc3KCN5ubmQqPRQKVSQaPRIC8vr8p+MTExUKlUUKlUiImJkdpPnjwJd3d3ODk5Yd68eRBCSNn8/Pzg5eUFHx8fpKamVlrfiRMn0L59e3zxxRcNyk9ERI1A1GLIkCHi1q1bwsvLS2pzc3OrbbEaRUREiBUrVgghhFixYoWIjIw06ZOTkyMcHR1FTk6OyM3NFY6OjiI3N1cIIYSvr684evSoKC8vF4GBgWL//v1CCCE0Go30et++fWL06NHS+oxGoxgzZowICgoSu3fvrlNOb2/vhuwmEVGbVNfvzlqPgDp06ICuXbtWamvXrk6PEapWXFwcwsLCAFQM8cXGxpr0SUhIgEajgY2NDbp37w6NRoP4+Hjo9XoUFBTAz88PCoUCoaGh0vIKhQIFBQUAgPz8fNjb20vre/fddzFt2jTeQJWIqJmo9RyQm5t2aS1jAAAXmElEQVQbPvnkE5SVlSEjIwPr16/H8OHDG7TR7Oxs6blCvXv3RnZ2tkkfnU6HPn36SO8dHByg0+mg0+ng4OBg0g4Aa9euhVarxYIFC1BeXo4jR45I69q7dy8OHDiAEydO1JgtOjoa0dHRAIAbN240aD+JiKh6tR7KvPvuu/j555/xwAMP4PHHH4e1tTXWrl1b64rHjRsHtVpt8hMXF1epn0KhgEKhuP89uMfGjRuxZs0aXL16FWvWrMGcOXMAAP/4xz+wcuXKOh25hYeHIy0tDWlpaejZs2ej5CIiIlO1HgF17twZy5cvx/Lly+u14qSkpGo/69WrF/R6Pezs7KDX66scFlMqlUhJSZHeZ2Vlwd/fH0qlEllZWZXalUolgIpJC+vWrQMAzJgxA3PnzgVQ8eC8WbNmAQBu3ryJ/fv3o3379pgyZUq99omIiBpPrYcE58+fR3h4OMaPH4+xY8dKPw0RHBwszWqLiYnB5MmTTfpotVokJiYiLy8PeXl5SExMhFarhZ2dHaytrXHs2DEIIbBjxw5peXt7exw8eBAAkJycDJVKBQC4fPkyMjMzkZmZienTp2PDhg0sPkREMqv1CGjGjBl47rnnMHfuXFhYWDTKRhctWoSZM2di69at6NevH3bt2gWg4khl06ZN2LJlC2xsbLBkyRL4+voCAF599VXY2NgAADZs2IDZs2fDYDAgKCgIQUFBAIDNmzdj/vz5MBqNsLS0lM7lEBFR86MQ4v8uoqmGt7c3Tp482VR5mhUfHx+kpaXJHYOIqEWp63dnrUNwkyZNwoYNG6DX65Gbmyv9EBERNUStQ3B3z9WsXr1aalMoFLh06ZL5UhERUatXawG6fPlyU+QgIqI2ptohuFWrVkmvd+/eXemzf/3rX+ZLREREbUK1Beizzz6TXq9YsaLSZ/Hx8eZLREREbUK1BejeyXF/nihXy8Q5IiKiWlVbgO69Pc6fb5XTWLfOISKitqvaSQinT5+GtbU1hBAwGAywtrYGUHH0U1JS0mQBiYiodaq2AJWVlTVlDiIiamMa9mAfIiKi+1TrdUBERM1N7A86rE44h2u3DLDv1gkRWmdMGayUOxbVEwsQEbUosT/osHjPGRhKK04T6G4ZsHjPGQBgEWphOARHRC3K6oRzUvG5y1BahtUJ52RKRPeLBYiIWpRrtwz1aqfmiwWIiFoU+26d6tVOzRcLEBG1KBFaZ3TqUPnhmJ06WCBC6yxTIrpfnIRARC3K3YkGnAXX8rEAEVGLM2WwkgWnFeAQHBERyYIFiIiIZMECREREsmABIiIiWbAAERGRLFiAiIhIFrIUoNzcXGg0GqhUKmg0GuTl5VXZLyYmBiqVCiqVCjExMVL7yZMn4e7uDicnJ8ybN096RHh6ejr8/Pzg5eUFHx8fpKamSsukpKTAy8sLbm5uGD16tHl3kIiIaidkEBERIVasWCGEEGLFihUiMjLSpE9OTo5wdHQUOTk5Ijc3Vzg6Oorc3FwhhBC+vr7i6NGjory8XAQGBor9+/cLIYTQaDTS63379onRo0cLIYTIy8sTLi4u4rfffhNCCJGdnV2nnN7e3g3aTyKitqiu352yHAHFxcUhLCwMABAWFobY2FiTPgkJCdBoNLCxsUH37t2h0WgQHx8PvV6PgoIC+Pn5QaFQIDQ0VFpeoVCgoKAAAJCfnw97e3sAwCeffIKQkBD07dsXAGBra9sUu0lERDWQ5U4I2dnZsLOzAwD07t0b2dnZJn10Oh369OkjvXdwcIBOp4NOp4ODg4NJOwCsXbsWWq0WCxYsQHl5OY4cOQIAOH/+PEpLS+Hv74/CwkLMnz8foaGhVWaLjo5GdHQ0AODGjRuNs8NERGTCbAVo3LhxuH79ukn78uXLK71XKBRQKBSNss2NGzdizZo1mDZtGnbt2oU5c+YgKSkJRqMRJ0+exLfffguDwYBhw4bBz88PgwYNMllHeHg4wsPDAQA+Pj6NkouIiEyZrQAlJSVV+1mvXr2g1+thZ2cHvV5f5ZCYUqlESkqK9D4rKwv+/v5QKpXIysqq1K5UVtwTKiYmBuvWrQMAzJgxA3PnzgVQcZTUo0cPdOnSBV26dMGoUaNw+vTpKgsQERE1DVnOAQUHB0uz2mJiYjB58mSTPlqtFomJicjLy0NeXh4SExOh1WphZ2cHa2trHDt2DEII7NixQ1re3t4eBw8eBAAkJydDpVIBACZPnozDhw/DaDTi9u3bOH78OFxcXJpob4mIqCqynANatGgRZs6cia1bt6Jfv37YtWsXACAtLQ2bNm3Cli1bYGNjgyVLlsDX1xcA8Oqrr8LGxgYAsGHDBsyePRsGgwFBQUEICgoCAGzevBnz58+H0WiEpaWldC7HxcUFgYGB8PDwQLt27TB37lyo1WoZ9pyIiO5SCPF/F9GQCR8fH6Slpckdg8jsYn/Q8fk61Gjq+t3J5wERtXGxP+iweM8ZGErLAAC6WwYs3nMGAFiEyKx4Kx6iNm51wjmp+NxlKC3D6oRzMiWitoIFiKiNu3bLUK92osbCAkTUxtl361SvdqLGwgJE1MZFaJ3RqYNFpbZOHSwQoXWWKRG1FZyEQNTG3Z1owFlw1NRYgIgIUwYrWXCoyXEIjoiIZMECREREsmABIiIiWbAAERGRLFiAiIhIFixAREQkCxYgIiKSBQsQERHJggWIiIhkwQJERESyYAEiIiJZsAAREZEsWICIiEgWLEBERCQLFiAiIpIFCxAREcmCBYiIiGTBAkRERLKQpQDl5uZCo9FApVJBo9EgLy+vyn4xMTFQqVRQqVSIiYmR2k+ePAl3d3c4OTlh3rx5EEIAANLT0+Hn5wcvLy/4+PggNTUVAJCfn49JkybB09MTbm5u2LZtm/l3ktq82B90GBGVDMdF+zAiKhmxP+jkjkTUrMhSgKKiohAQEICMjAwEBAQgKirKpE9ubi6WLVuG48ePIzU1FcuWLZMK1fPPP4/NmzcjIyMDGRkZiI+PBwBERkZi6dKlSE9Px+uvv47IyEgAwPvvvw9XV1ecPn0aKSkpePnll/HHH3803Q5TmxP7gw6L95yB7pYBAoDulgGL95xhESK6hywFKC4uDmFhYQCAsLAwxMbGmvRJSEiARqOBjY0NunfvDo1Gg/j4eOj1ehQUFMDPzw8KhQKhoaHS8gqFAgUFBQAqjnrs7e2l9sLCQgghUFRUBBsbG7Rv376J9pbaotUJ52AoLavUZigtw+qEczIlImp+ZPkWzs7Ohp2dHQCgd+/eyM7ONumj0+nQp08f6b2DgwN0Oh10Oh0cHBxM2gFg7dq10Gq1WLBgAcrLy3HkyBEAwEsvvYTg4GDY29ujsLAQn3/+Odq1q7r2RkdHIzo6GgBw48aNxtlhanOu3TLUq52oLTLbEdC4ceOgVqtNfuLi4ir1UygUUCgUjbLNjRs3Ys2aNbh69SrWrFmDOXPmAKg4mvLy8sK1a9eQnp6Ol156STpS+rPw8HCkpaUhLS0NPXv2bJRc1PbYd+tUr3aitshsBSgpKQk//fSTyc/kyZPRq1cv6PV6AIBer4etra3J8kqlElevXpXeZ2VlQalUQqlUIisry6QdqJi0EBISAgCYMWOGNAlh27ZtCAkJgUKhgJOTExwdHfHrr7+aa9eJEKF1RqcOFpXaOnWwQITWWaZERM2PLOeAgoODpVltMTExmDx5skkfrVaLxMRE5OXlIS8vD4mJidBqtbCzs4O1tTWOHTsGIQR27NghLW9vb4+DBw8CAJKTk6FSqQAAffv2xbfffgugYvjv3LlzGDBgQFPsKjWiljSrbMpgJVaEuEPZrRMUAJTdOmFFiDumDFbKHY2o2VCIu3OYm1BOTg5mzpyJK1euoF+/fti1axdsbGyQlpaGTZs2YcuWLQCADz/8EP/5z38AAK+88gqeeuopAEBaWhpmz54Ng8GAoKAgvPvuu1AoFDh8+DDmz58Po9EIS0tLbNiwAd7e3rh27Rpmz54NvV4PIQQWLVqEJ598stacPj4+SEtLM98vgurs7qyye0/sd+pgwS91omaort+dshSgloIFqPkYEZUMXRUn8JXdOuH7RWNlSERE1anrdyfvhEAtAmeVEbU+LEDUInBWGVHrwwJELQJnlRG1PrwdALUIdycarE44h2u3DLDv1gkRWmdOQCBqwViAqMWYMljJgkPUinAIjoiIZMEjoEYW+4OuRQ0TtbS8RNR6sAA1oj9fLHn3FvwAmuWXekvLS0StC4fgGlFLuwV/S8tLRK0LC1AjamkXS7a0vETUurAANaKWdrFkS8tLRK0LC1AjamkXS7a0vETUunASQiNqaRdLtrS8RNS6sAA1spZ2sWRLy0tErQeH4IiISBYsQEREJAsWICIikgULEBERyYIFiIiIZKEQQgi5QzRXDz30EPr3739fy964cQM9e/Zs3EBm1JLytqSsQMvK25KyAi0rb0vKCjQsb2ZmJm7evFlrPxYgM/Hx8UFaWprcMeqsJeVtSVmBlpW3JWUFWlbelpQVaJq8HIIjIiJZsAAREZEsLF577bXX5A7RWnl7e8sdoV5aUt6WlBVoWXlbUlagZeVtSVkB8+flOSAiIpIFh+CIiEgWLEBERCQLFqD7EB8fD2dnZzg5OSEqKsrk83feeQeurq7w8PBAQEAAfvvtt0qfFxQUwMHBAS+99FKzznrlyhWMHz8eLi4ucHV1RWZmZrPOGxkZCTc3N7i4uGDevHkw9+hybVk3bdoEd3d3eHl54dFHH8XZs2elz1asWAEnJyc4OzsjISHBrDkbmvd///sfvL294e7uDm9vbyQnJzfbrHdduXIFVlZWeOutt8yetaF5f/zxRwwbNgxubm5wd3dHSUlJs8xaWlqKsLAwuLu7w8XFBStWrGh4GEH1YjQaxYABA8TFixfFnTt3hIeHh/j5558r9UlOThbFxcVCCCE2bNggZs6cWenzefPmiccff1y8+OKLzTrr6NGjRWJiohBCiMLCQqlfc8z7/fffi+HDhwuj0SiMRqPw8/MTBw4ckDVrfn6+9DouLk5otVohhBA///yz8PDwECUlJeLSpUtiwIABwmg0mi1rQ/OeOnVK6HQ6IYQQZ86cEfb29s02613Tpk0T06dPF6tXrzZr1obmLS0tFe7u7iI9PV0IIcTNmzfN+rfQkKwff/yxeOyxx4QQQhQXF4t+/fqJy5cvNygPj4DqKTU1FU5OThgwYAA6duyIWbNmIS4urlKfMWPGoHPnzgAAPz8/ZGVlSZ+dPHkS2dnZGD9+fLPOevbsWRiNRmg0GgCAlZWV1K855lUoFCgpKcEff/yBO3fuoLS0FL169ZI1q7W1tfS6uLgYCoUCABAXF4dZs2bhgQcegKOjI5ycnJCammq2rA3NO3jwYNjb2wMA3NzcYDAYcOfOnWaZFQBiY2Ph6OgINzc3s2VsrLyJiYnw8PCAp6cnAKBHjx6wsKj8lOLmklWhUKC4uBhGoxEGgwEdO3as1Pd+sADVk06nQ58+faT3Dg4O0Ol01fbfunUrgoKCAADl5eV4+eWXm2xYoCFZz58/j27duiEkJASDBw9GREQEysrKmm3eYcOGYcyYMbCzs4OdnR20Wi1cXFxkz/r+++9j4MCBiIyMxPr16+u1bHPJe68vv/wSQ4YMwQMPPNAssxYVFWHlypVYunSp2fI1Zt7z589DoVBAq9ViyJAhWLVqVbPNOn36dHTp0gV2dnbo27cvFixYABsbmwblYQEyo48++ghpaWmIiIgAAGzYsAETJkyAg4ODzMlM/Tmr0WjEoUOH8NZbb+HEiRO4dOkStm/fLm/Ie/w574ULF/DLL78gKysLOp0OycnJOHTokMwpgRdffBEXL17EypUr8eabb8odp1Y15f3555+xcOFCfPDBBzKlq6yqrK+99hr++c9/wsrKSuZ0pqrKazQacfjwYXz88cc4fPgw9u7di2+//VbmpFVnTU1NhYWFBa5du4bLly/j7bffxqVLlxq0HRagelIqlbh69ar0PisrC0ql6SOtk5KSsHz5cnz11VfS/y0ePXoU7733Hvr3748FCxZgx44dWLRoUbPM6uDgAC8vLwwYMADt27fHlClTcOrUKbNlbWjevXv3ws/PD1ZWVrCyskJQUBCOHj0qe9a7Zs2ahdjY2PtatjE0JO/d/lOnTsWOHTswcODAZpv1+PHjiIyMRP/+/bF27Vr85z//wXvvvdds8zo4OGDUqFF46KGH0LlzZ0yYMMGs/84akvWTTz5BYGAgOnToAFtbW4wYMaLh94pr0BmkNqi0tFQ4OjqKS5cuSSfxfvrpp0p9Tp06JQYMGCDOnz9f7Xq2bdtm9kkIDclqNBqFh4eH+P3334UQQsyePVu89957zTbvZ599JgICAkRpaan4448/xNixY8VXX30la9Z7M3711VfC29tbCCHETz/9VGkSgqOjo9knITQkb15envDw8BBffvmlWTM2RtZ7LV26tEkmITQkb25urhg8eLAoLi4WpaWlIiAgQHzzzTfNMmtUVJSYPXu2EEKIoqIi4eLiIk6fPt2gPCxA92Hfvn1CpVKJAQMGiDfffFMIIcSSJUtEXFycEEKIgIAAYWtrKzw9PYWnp6eYNGmSyTqaogA1NGtiYqJwd3cXarVahIWFiTt37jTbvEajUYSHh4uHH35YuLi4iH/+85+yZ503b55wdXUVnp6ewt/fv9I/9DfffFMMGDBADBo0SOzfv9/sWRuS94033hCdO3eWfueenp4iOzu7WWa9V1MVoIbm3blzp3B1dRVubm4iIiKi2WYtLCwU06dPF66ursLFxUWsWrWqwVl4Kx4iIpIFzwEREZEsWICIiEgWLEBERCQLFiAiIpIFCxAREcmCBYiIiGTBAkRERLJgASK6D5mZmVCr1Q1ax/Dhw6tsf+211+p9w9rmeO8zotq0lzsAUVsjKu5AgiNHjsgdhUhWPAIiuk9lZWV45pln4ObmhvHjx8NgMAAA3njjDTg7O+PRRx/F448/jrfeeguZmZlwdnZGaGgo1Go1rl69WumoZfny5Rg0aBAeffRRnDt3rtptTpkyBd7e3nBzc0N0dHSVfd555x2o1Wqo1WqsXbsWQMURm4uLS53z/tmYMWPwv//9DwDw73//G3//+9/v75dGdK8G38yHqA26fPmysLCwED/88IMQQogZM2aInTt3itTUVOHp6SkMBoMoKCgQTk5OYvXq1eLy5ctCoVCIo0ePSuvo0qWLEEKItLQ0oVarRXFxscjPzxcDBw6s9h5mOTk5Qgghbt++Ldzc3MTNmzerXFdRUZEoLCwUrq6u4tSpU/XO+2cHDx4Uo0ePFh999JGYMGGC2W+eSm0Dh+CI7pOjoyO8vLwAAN7e3sjMzMTNmzcxefJkWFpawtLSEpMmTZL69+vXD35+fibrOXToEKZOnSo96TU4OLjaba5fvx579+4FAFy9ehUZGRno0aOH9Pnhw4cxdepUdOnSBQAQEhKCQ4cOITg4uN557zVq1CgIIfDOO+8gJSXFrE/tpLaDQ3BE9+nep4JaWFjAaDTW2P9uUbhfKSkpSEpKwtGjR3H69GkMHjwYJSUldV6+vnnvdebMGej1enTs2BEPPvhgvXITVYcFiKgRjRgxAl9//TVKSkpQVFSEb775ptZlRo0ahdjYWBgMBhQWFuLrr7+usl9+fj66d++Ozp0749dff8WxY8dM+owcORKxsbG4ffs2iouLsXfvXowcObJBefV6PZ544gnExcXBysoK8fHxte4TUV1wCI6oEfn6+iI4OBgeHh7o1asX3N3d0bVr1xqXGTJkCB577DF4enrC1tYWvr6+VfYLDAzEpk2b4OLiAmdn5yqH84YMGYLZs2dj6NChAIC5c+di8ODByMzMvK+8t2/fRkhICN5++224uLhgyZIlWLhwIQIDA+v4GyGqHp8HRNTIioqKYGVlhdu3b2PUqFGIjo7GkCFD5I5VrZaWl1oPHgERNbLw8HCcPXsWJSUlCAsLa/Zf5i0tL7UePAIiIiJZcBICERHJggWIiIhkwQJERESyYAEiIiJZsAAREZEsWICIiEgWLEBERCSL/wcVVJOFzVeELwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "energies = [job.logfile.energy\n", " for job in hgc.queue if job.is_completed]\n", "hgrids = [job.param \n", " for job in hgc.queue if job.is_completed]\n", "simple_hgrids = [hgrid[0] for hgrid in hgrids]\n", "n_at = len(hgc.queue[0].posinp)\n", "threshold = min(energies) + n_at*hgc.precision_per_atom\n", "\n", "fig=plt.figure()\n", "fig.patch.set_facecolor('white') # When dark background\n", "plt.plot(simple_hgrids, energies,\n", " label=\"BigDFT results\", marker=\"o\", linestyle='')\n", "plt.plot(simple_hgrids, [threshold]*len(energies),\n", " label=\"Convergence threshold\")\n", "plt.xlabel(\"hgrid along $x$\")\n", "plt.ylabel(\"Energy (Ha)\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## frmult convergence\n", "\n", "We started the convergence study with a coarse grid having a larger extension than the fine grid. This means that the coarse grid extension is converged, but there might still be some room to use a smaller fine grid. This is the reason why the convergence study ends by looking only at the fine grid extension multiplying factor `frmult` (note that the delta value for `crmult` is 0):" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "[frm, crm] = rmc.converged.param\n", "frmc = RmultConvergence(base, reference=[frm+2, crm], delta=[-1, 0],\n", " n_jobs=5, precision_per_atom=0.01*EV_TO_HA)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/7.0_7.0\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/6.0_7.0\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/5.0_7.0\n", "Logfile log-N2.yaml already exists!\n", "\n", "/Users/maximemoriniere/Documents/Python/MyBigDFT/doc/source/notebooks/N2/hgrids_convergence/rm_5.0_7.0/4.0_7.0\n", "Logfile log-N2.yaml already exists!\n", "\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/Users/maximemoriniere/Documents/Python/MyBigDFT/mybigdft/workflows/workflow.py:114: UserWarning: Some jobs of the workflow were not run.\n", " UserWarning)\n" ] } ], "source": [ "frmc.run(nmpi=6, nomp=3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nothing changed: `frmult = 5.0` gives converged results for the given precision:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requested precision per atom: 3.67e-04 (Ha)\n", "---------------------------------------------------\n", " rmult precision_per_atom (Ha) is_converged \n", "---------------------------------------------------\n", " [7.0, 7.0] 0.00e+00 True \n", " [6.0, 7.0] 1.17e-05 True \n", " [5.0, 7.0] 1.55e-04 True \n", " [4.0, 7.0] 1.93e-03 False \n" ] } ], "source": [ "frmc.summary()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAERCAYAAAA9oHOJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X9UVVX+//HnVcw0A8VU8GJqcjPl11WgqMZS8YqYoZaVjRM0H4myb+NMa4akabSxz2JkbOYz1ZSZVgrzIyozsSwgStImlcCoSRu7mZT3eiUV/JmW4vn+4cfzkfipAkfx9ViLtTj77rPPe3fsvjnn7LO3zTAMAxERkTbWweoARETkwqQEJCIillACEhERSygBiYiIJZSARETEEkpAIiJiCSWgs/Dqq68SFhZGhw4dKC0tbbDek08+SXh4OGFhYTzxxBNmeXl5OXFxcTidTmJiYigpKQFg37593HzzzURFRREWFsaSJUuajOW//uu/6N27N+Hh4WffMRGRNqAEdBbCw8NZvnw5N9xwQ4N1PvvsMxYvXkxJSQmffPIJb775Jl9++SUADz30EI8++ijl5eU89thjPPTQQwA888wzDB06lE8++YTi4mJ+/etf88MPPzQay913301+fn7LdU5EpJUpAZ2FIUOGMHjw4EbrfP7551xzzTV07doVPz8/brzxRpYvXw6AzWZj//79wImrnr59+5rlBw4cwDAMDh48SGBgIH5+fgA8/vjjxMbGEhkZyaOPPmoe54YbbiAwMLA1uiki0ir8rA6gvQsPD+eRRx5hz549dOnShbfeeouYmBgAnnjiCRISEvjNb37D8ePH+fDDDwF44IEHSEpKom/fvhw4cICXX36ZDh06UFhYiNvtpqSkBMMwSEpKYs2aNY1egYmInKuUgJowZswYdu7cWac8MzOTiRMnNrn/kCFDmDVrFmPHjuWSSy7B6XTSsWNHAJ599ln+8pe/cOutt/LKK68wffp0ioqKKCgowOl08t5777F161ZcLhcjRoygsLCQwsJChg0bBsDBgwdxu91KQCJyXlICakJRUdFZtzF9+nSmT58OwG9/+1tCQkIAyM7O5sknnwTgtttuIzU1FYAlS5aQkZGBzWYjNDSUgQMH8p///AfDMHj44Ye59957zzomERGr6RlQG/j2228B+Oabb1i+fDk//elPAejbty/vv/8+AO+99x4OhwOAyy+/nHfffReAyspKtmzZwhVXXEFCQgIvvvgiBw8eBMDr9Zpti4icdww5Y8uXLzfsdrtx0UUXGb179zbGjh1rGIZheL1eIzEx0az3k5/8xBgyZIgRGRlpFBUVmeVr1641hg8fbkRGRhpXX321UVpaau7vcrmM8PBwIywszPjb3/5m7vPEE08Y4eHhRnh4uBEXF2d8+eWXhmEYxtSpU42goCDDz8/PsNvtxvPPP98W/wlERM6YzTC0HIOIiLQ93YITERFLaBBCIy677DIGDBhgdRgiIueViooKdu/e3WQ9JaBGDBgwoNEpdkREpK6T7zo2RbfgRETEEkpAIiJiCSUgERGxhJ4BibRzR48exePxcOTIEatDkXbm4osvJiQkhE6dOp3R/kpAIu2cx+Ph0ksvZcCAAdhsNqvDkXbCMAz27NmDx+Nh4MCBZ9SGElArWPGxl8cLtrBj72H6du9CesJgJg2zWx2WXKCOHDmi5CMtzmaz0bNnT3bt2nXGbSgBtbAVH3t5ePm/OXy0BgDv3sM8vPzfAEpCYhklH2kNZ/vvSoMQWtjjBVvM5HPS4aM1PF6wxaKIRETOTUpALWzH3sOnVS5yIejYsSNOp5OoqCiGDx9uLr64Y8cOpkyZ0uT+AwYMICIigoiICIYOHcrvfvc7c1BFRUUFXbp0wel0mj/PPfec+ftFF11EREQETqeTjIyMVu3nSUuXLuWBBx4AYMWKFWzevLlNjnu+0S24Fta3exe89SSbvt27WBCNyOlrjWeYXbp0oby8HICCggIefvhh3n//ffr27cuyZcua1cbq1au57LLLOHjwIGlpadx7771kZ2cDMGjQILP9k06umzVgwABz38YcO3YMP7+W/0pcsWIFEyZMYOjQoS3e9vlOV0AtLD1hMF06daxV1qVTR9ITBlsUkUjznXyG6d17GIP/e4a54mNvix1j//799OjRAzhx9RIeHg7Ad999x+23387QoUOZPHky11xzTb1TYXXr1o2FCxeyYsUKqqqqziqWpUuXkpSUxOjRo4mPjwfg8ccfJzY2lsjISB599FEADh06xE033URUVBTh4eG8/PLLwInkdnLOs9LSUkaOHFmr/Q8//JCVK1eSnp6O0+lk69atPPXUUwwdOpTIyEimTp16VvGf73QF1MJO/qWoUXByPmrsGebZ/Bs+fPgwTqeTI0eO4PP5eO+99+rUWbBgAT169GDz5s189tlnOJ3OBtvz9/dn4MCBuN1u+vTpw9atW836119/Pc8880yzY9u4cSOffvopgYGBFBYW4na7KSkpwTAMkpKSWLNmDbt27aJv376sWrUKgH379jWr7euuu46kpCQmTJhg3mrMyspi27ZtdO7cmb179zY7zvZICagVTBpmV8KR81JrPcM89RbcunXrSE5O5rPPPqtV54MPPuCXv/wlAOHh4URGRjba5qlLmdV3C665XC4XgYGBABQWFlJYWMiwYcMAOHjwIG63mxEjRvDrX/+aWbNmMWHCBEaMGHFGxwKIjIxk2rRpTJo0iUmTJp1xO+2BbsGJiKmhZ5Ut+Qzz2muvZffu3Wf1/siBAweoqKjgyiuvPOt4LrnkEvN3wzB4+OGHKS8vp7y8nC+//JLp06dz5ZVXsnHjRiIiIvjd737HY489BoCfnx/Hjx8HaPZME6tWreL//b//x8aNG4mNjeXYsWNn3YfzlRKQiJja4hnmf/7zH2pqaujZs2et8uuvv55XXnkFgM2bN/Pvf/+73v0PHjzI/fffz6RJk8xnSS0lISGBF198kYMHDwLg9Xr59ttv2bFjB127duVnP/sZ6enpbNy4ETjxDKisrAyA1157rd42L730Ug4cOADA8ePH2b59O6NGjeKPf/wj+/btM491IdItOBExtdYzzJPPgODEVUZ2djYdO9ZOdPfffz8pKSkMHTqUq666irCwMAICAszPR40ahWEYHD9+nMmTJzN79uyziqk+Y8eO5fPPP+faa68FTgx4+Pvf/86XX35Jeno6HTp0oFOnTjz77LMAPProo0yfPp3Zs2fXGYBw0tSpU7nnnnt46qmnyM3NZfr06ezbtw/DMJg5cybdu3dv8X6cL2zGqTdSpZaYmBgtSCfnvc8//5whQ4ZYHUaTampqOHr0KBdffDFbt25lzJgxbNmyhYsuusjq0KQR9f37au53p66AROSc8N133zFq1CiOHj2KYRgsWLBAyaeds+QZUFVVFS6XC4fDgcvlorq6ut562dnZOBwOHA6H+cIZQFlZGREREYSGhjJz5kxzNMyrr75KWFgYHTp0qJN9582bR2hoKIMHD6agoKD1OiciZ+TSSy+ltLSUTz75hE8//ZTExESrQ5JWZkkCysrKIj4+HrfbTXx8PFlZWXXqVFVVMXfuXDZs2EBJSQlz5841E9WMGTNYvHgxbrcbt9tNfn4+cGLo5vLly7nhhhtqtbV582Zyc3PZtGkT+fn53H///dTU1NQ5poiItB1LElBeXh4pKSkApKSksGLFijp1CgoKzPH5PXr0wOVykZ+fj8/nY//+/cTFxWGz2UhOTjb3HzJkCIMH1x2tk5eXx9SpU+ncuTMDBw4kNDSUkpKS1u2kiIg0ypIEVFlZSXBwMABBQUFUVlbWqeP1eunXr5+5HRISgtfrxev1EhISUqe8MQ21JSIi1mm1QQhjxoxh586ddcozMzNrbdtstnNqrZJFixaxaNEigLN6UU5ERBrXaldARUVFfPbZZ3V+Jk6cSJ8+ffD5fAD4fD569+5dZ3+73c727dvNbY/Hg91ux2634/F46pQ3pqG26pOWlkZpaSmlpaX06tXrtPosIvXbuXMnU6dOZdCgQURHRzN+/Hi++OILq8OyzN69e1mwYIG5XVxczIQJE1r8OKcuC9Fcp06weqrf//73/OlPf2qp0ACLbsElJSWZo9qys7OZOHFinToJCQkUFhZSXV1NdXU1hYWFJCQkEBwcjL+/P+vXr8cwDHJycurd/8fHy83N5fvvv2fbtm243W6uvvrqVumbiNRmGAaTJ09m5MiRbN26lbKyMubNm1fvrffWdq4MPvpxAmqucyX+lmJJAsrIyOCdd97B4XBQVFRkLhJVWlpKamoqAIGBgcyePZvY2FhiY2OZM2eOOWHgggULSE1NJTQ0lEGDBpnDNV9//XVCQkJYt24dN910EwkJCQCEhYWZ07yPGzeOZ555ps5b2CLSOlavXk2nTp247777zLKoqChGjBiBYRikp6cTHh5ORESEucxBcXExI0eOZMqUKVx11VVMmzYNwzDIz8/ntttuM9s59cqhsLCQa6+9luHDh3PbbbeZU9wMGDCAWbNmMXz4cF599VU++ugjIiMjcTqd5rHhxJd7enq6uRTDc88912gsAB999BHXXXcdUVFRXH311Rw4cKDBdk6VkZFhzuCdnp4OnJhiqL5j/Dj+rVu3Mm7cOKKjoxkxYgT/+c9/gBOvoYSHhxMVFVVrJPCOHTsYN24cDoeDhx56yCx/6aWXiIiIIDw8nFmzZtV77jIzM7nyyiv5yU9+wpYtrbCqsyENio6OtjoEkbO2efPm/9t4a5ZhvDi+ZX/emtXo8Z988knjV7/6Vb2fLVu2zBgzZoxx7NgxY+fOnUa/fv2MHTt2GKtXrzb8/f2N7du3GzU1NUZcXJyxdu1a4+jRo0a/fv2MgwcPGoZhGPfdd5/xt7/9zdi1a5cxYsQIszwrK8uYO3euYRiG0b9/f+OPf/yjecywsDDjww8/NAzDMGbNmmWEhYUZhmEYzz33nPHf//3fhmEYxpEjR4zo6Gjjq6++ajCW77//3hg4cKBRUlJiGIZh7Nu3zzh69GiD7Zxq27Zt5nENw2jwGPXFP3r0aOOLL74wDMMw1q9fb4waNcowDMMIDw83PB6PYRiGUV1dbRiGYSxZssQYOHCgsXfvXuPw4cPG5ZdfbnzzzTeG1+s1+vXrZ3z77bfG0aNHjVGjRhmvv/66ebxdu3YZpaWlRnh4uHHo0CFj3759xqBBg4zHH3+8zjms9e/rfzX3u1MzIYiIZT744APuvPNOOnbsSJ8+fbjxxhv56KOP8Pf35+qrrzZHvDqdTioqKvjJT37CuHHjeOONN5gyZQqrVq1i/vz5vP/++2zevJnrr78egB9++MGczw3gjjvuAE7c+jpw4ID52U9/+lPefPNN4MQV1Keffmqu0Lpv3z7cbjcXXXRRvbEEBAQQHBxMbGwscGKNosbaGThwYKP/LRrq76nxHzx4kA8//LDWVeD3338PnJjM9e677+b222/nlltuMT+Pj48359QbOnQoX3/9NXv27GHkyJHmc+5p06axZs2aWstDrF27lsmTJ9O1a1fgxKOMlqYEJHIhSaz70ndrCwsLa/ay26fq3Lmz+XvHjh3NZQumTp3K008/TWBgIDExMVx66aUYhoHL5eKll16qt61Tl1xoiGEY/PWvfzVv3Z9UXFzcYCyn005TGjvGyfiPHz9O9+7d6137aOHChWzYsIFVq1YRHR1tztJ9OrG3NS3HICKtavTo0Xz//ffm6w0An376KWvXrmXEiBG8/PLL1NTUsGvXLtasWdPkAKEbb7yRjRs3snjxYnNJ67i4OP71r3/x5ZdfAieW0K5vlF337t259NJL2bBhAwC5ubnmZwkJCTz77LMcPXoUgC+++IJDhw41GMfgwYPx+Xx89NFHwIk1io4dO9asdk5douF0nFwJ9tVXXwVOJLtPPvkEgK1bt3LNNdfw2GOP0atXr1ojf3/s6quv5v3332f37t3U1NTw0ksvceONN9aqc8MNN7BixQoOHz7MgQMHeOONN0473qYoAYlIq7LZbLz++usUFRUxaNAgwsLCePjhhwkKCmLy5MlERkYSFRXF6NGjmT9/PkFBQY2217FjRyZMmMDbb79tDkDo1asXS5cu5c477yQyMpJrr73WfDj/Yy+88AL33HMPTqeTQ4cOmbenUlNTGTp0KMOHDyc8PJx777230auFiy66iJdffplf/OIXREVF4XK5OHLkSLPa6dmzJ9dffz3h4eHmIITm+sc//sELL7xAVFQUYWFh5OXlAZCenm4OKjg5MKIhwcHBZGVlMWrUKKKiooiOjq4zmnj48OHccccdREVFkZiYaN5qbElajqERWo5B2oPzZTmGtnLw4EG6desGnJiX0ufz8eSTT1oc1flLyzGIiDTTqlWrmDdvHseOHaN///4sXbrU6pAuWEpAInJBueOOO8xRZWItPQMSuQDoTru0hrP9d6UEJNLOXXzxxezZs0dJSFqUYRjs2bOHiy+++Izb0C04kXYuJCQEj8ej2d2lxV188cW1lsc5XUpAIu1cp06dmnwLX8QKugUnIiKWUAISERFLKAGJiIgllIBERMQSSkAiImIJJSAREbGEJQmoqqoKl8uFw+HA5XJRXV1db73s7GwcDgcOh4Ps7GyzvKysjIiICEJDQ5k5c6b5gt2rr75KWFgYHTp0qDURXkVFBV26dMHpdOJ0OmstDSwiItawJAFlZWURHx+P2+0mPj6erKy6i2RVVVUxd+5cNmzYQElJCXPnzjUT1YwZM1i8eDFutxu3201+fj4A4eHhLF++vNZ66CcNGjSI8vJyysvLWbhwYet2UEREmmRJAsrLyyMlJQWAlJQUVqxYUadOQUEBLpeLwMBAevTogcvlIj8/H5/Px/79+4mLi8Nms5GcnGzuP2TIEAYPHtymfRERkTNjSQKqrKwkODgYgKCgICorK+vU8Xq99OvXz9wOCQnB6/Xi9XprTf1wsrwp27ZtY9iwYdx4442sXbu2wXqLFi0iJiaGmJgYTV0iItKKWm0qnjFjxrBz58465ZmZmbW2bTYbNputtcIATqz+980339CzZ0/KysqYNGkSmzZtwt/fv07dtLQ00tLSgBOLKomISOtotQRUVFTU4Gd9+vTB5/MRHByMz+ejd+/ederY7XaKi4vNbY/Hw8iRI7Hb7Xg8nlrldru90Vg6d+5M586dAYiOjmbQoEF88cUXSjAiIhay5BZcUlKSOaotOzu7zlrkAAkJCRQWFlJdXU11dTWFhYUkJCQQHByMv78/69evxzAMcnJy6t3/VLt27aKmpgaAr776CrfbzRVXXNHyHRMRkWazJAFlZGTwzjvv4HA4KCoqIiMjA4DS0lJSU1MBCAwMZPbs2cTGxhIbG8ucOXMIDAwEYMGCBaSmphIaGsqgQYNITEwE4PXXXyckJIR169Zx0003kZCQAMCaNWuIjIzE6XQyZcoUFi5caLYlIiLWsBlapapBMTExtd4nEhGRpjX3u1MzIYiIiCWUgERExBJKQCIiYgklIBERsYQSkIiIWEIJSERELKEEJCIillACEhERSygBiYiIJZSARETEEkpAIiJiCSUgERGxhBKQiIhYQglIREQsoQQkIiKWUAISERFLKAGJiIglLElAVVVVuFwuHA4HLpeL6urqeutlZ2fjcDhwOBxkZ2eb5WVlZURERBAaGsrMmTM5uahreno6V111FZGRkUyePJm9e/ea+8ybN4/Q0FAGDx5MQUFB63ZQRESaZEkCysrKIj4+HrfbTXx8PFlZWXXqVFVVMXfuXDZs2EBJSQlz5841E9WMGTNYvHgxbrcbt9tNfn4+AC6Xi88++4xPP/2UK6+8knnz5gGwefNmcnNz2bRpE/n5+dx///3U1NS0XYdFRKQOSxJQXl4eKSkpAKSkpLBixYo6dQoKCnC5XAQGBtKjRw9cLhf5+fn4fD72799PXFwcNpuN5ORkc/+xY8fi5+cHQFxcHB6Pxzze1KlT6dy5MwMHDiQ0NJSSkpI26q2IiNTHkgRUWVlJcHAwAEFBQVRWVtap4/V66devn7kdEhKC1+vF6/USEhJSp/zHXnzxRRITExttqz6LFi0iJiaGmJgYdu3adWYdFBGRJvm1VsNjxoxh586ddcozMzNrbdtsNmw2W4seOzMzEz8/P6ZNm3ba+6alpZGWlgZATExMi8YlIiL/p9USUFFRUYOf9enTB5/PR3BwMD6fj969e9epY7fbKS4uNrc9Hg8jR47Ebrebt9ZOltvtdnN76dKlvPnmm7z77rtmYrPb7Wzfvr3BfUREpO1ZcgsuKSnJHNWWnZ3NxIkT69RJSEigsLCQ6upqqqurKSwsJCEhgeDgYPz9/Vm/fj2GYZCTk2Pun5+fz/z581m5ciVdu3atdbzc3Fy+//57tm3bhtvt5uqrr26bzoqISL0sSUAZGRm88847OBwOioqKyMjIAKC0tJTU1FQAAgMDmT17NrGxscTGxjJnzhwCAwMBWLBgAampqYSGhjJo0CDzWc8DDzzAgQMHcLlcOJ1O7rvvPgDCwsK4/fbbGTp0KOPGjeOZZ56hY8eOFvRcREROshknX6KROmJiYigtLbU6DBGR80pzvzs1E4KIiFii0UEIHo+H3Nxc1q5dy44dO+jSpQvh4eHcdNNNJCYm0qGD8peIiJyZBhPQz3/+c7xeLxMmTGDWrFn07t2bI0eO8MUXX5Cfn09mZiZZWVnccMMNbRmviIi0Ew0moF//+teEh4fXKQ8PD+eWW27hhx9+4JtvvmnV4EREpP1q8B5afcnnVBdddBGhoaEtHpCIiFwYmnwR1e128/DDD7N582aOHDliln/11VetGpiIiLRvTY4i+PnPf86MGTPw8/Nj9erVJCcn87Of/awtYhMRkXasyQR0+PBh4uPjMQyD/v378/vf/55Vq1a1RWwiItKONXkLrnPnzhw/fhyHw8HTTz+N3W7n4MGDbRGbiIi0Y01eAT355JN89913PPXUU5SVlfG3v/2t1uqkIiIiZ6LJK6DY2FgAunXrxpIlS1o9IBERuTA0mIBuvvnmRtfpWblyZasEJCIiF4YGE9BvfvMbAAzD4J577uH5559vs6BERKT9azAB3Xjjjebv3bp1q7UtIiJytpo1m2hLL5ktIiLS4BVQVVWV+XtNTQ3V1dWcunTQycXhREREzkSDV0DR0dHExMQQHR3N/v37GT58ONHR0Wb52aiqqsLlcuFwOHC5XFRXV9dbLzs7G4fDgcPhqDX0u6ysjIiICEJDQ5k5c6aZGNPT07nqqquIjIxk8uTJ7N27F4CKigq6dOmC0+mstVKqiIhYyLBAenq6MW/ePMMwDGPevHnGQw89VKfOnj17jIEDBxp79uwxqqqqjIEDBxpVVVWGYRhGbGyssW7dOuP48ePGuHHjjLfeesswDMMoKCgwjh49ahiGYTz00ENmu9u2bTPCwsJOO87o6Ogz6p+IyIWsud+dDV4BVVRUNJW48Hg8Z5T08vLySElJASAlJYUVK1bUqVNQUIDL5SIwMJAePXrgcrnIz8/H5/Oxf/9+4uLisNlsJCcnm/uPHTsWP78TdxXj4uLOOD4REWl9DSag9PR0br31VnJycti0aRPffvst33zzDe+99x6zZ8/m+uuv5/PPPz+jg1ZWVhIcHAxAUFAQlZWVdep4vV769etnboeEhOD1evF6vYSEhNQp/7EXX3yRxMREc3vbtm0MGzaMG2+8kbVr1zYY26JFi4iJiSEmJoZdu3adUf9ERKRpDQ5CePXVV9m8eTP/+Mc/ePHFF/H5fHTt2pUhQ4Ywfvx4HnnkES6++OIGGx4zZgw7d+6sU56ZmVlr22aztfgou8zMTPz8/Jg2bRoAwcHBfPPNN/Ts2ZOysjImTZrEpk2b8Pf3r7NvWloaaWlpAGf9rEtERBrW6FQ8Q4cOrZMwmquoqKjBz/r06YPP5yM4OBifz0fv3r3r1LHb7RQXF5vbHo+HkSNHYrfba91a83g82O12c3vp0qW8+eabvPvuu2Zi69y5M507dwZODK4YNGgQX3zxhRKMiIiFmvUeUEtLSkoyR7VlZ2czceLEOnUSEhIoLCykurqa6upqCgsLSUhIIDg4GH9/f9avX49hGOTk5Jj75+fnM3/+fFauXEnXrl3Ntnbt2kVNTQ1wYiE9t9vNFVdc0QY9FRGRhliSgDIyMnjnnXdwOBwUFRWRkZEBQGlpKampqcCJ94xmz55NbGwssbGxzJkzx3z3aMGCBaSmphIaGsqgQYPMZz0PPPAABw4cwOVy1RpuvWbNGiIjI3E6nUyZMoWFCxfqPSYREYvZDOOUt0ullpiYGEpLS60OQ0TkvNLc784mr4BuueUWVq1axfHjx1skMBEREWhGArr//vv55z//icPhICMjgy1btrRFXCIi0s41mYDGjBnDP/7xDzZu3MiAAQMYM2YM1113HUuWLOHo0aNtEaOIiLRDzRqEsGfPHpYuXcrzzz/PsGHD+OUvf8nGjRtxuVytHZ+IiLRTTS7JPXnyZLZs2cJdd93FG2+8Yc5gcMcdd+g9GhEROWNNJqCZM2cyatSoej/TCDERETlTTSag6upqli9fXqssICCAiIiIemcwEBERaY4mE9ALL7zAunXrzKug4uJioqOj2bZtG3PmzOGuu+5q9SBFRKT9aTIBHT16lM8//5w+ffoAJ2ayTk5OZsOGDdxwww1KQCIickaaHAXn8XjM5APQu3dvtm/fTmBgIJ06dWrV4EREpP1q8gpo5MiRTJgwgdtuuw2A1157jZEjR3Lo0CG6d+/e6gGKiEj71GQCeuaZZ1i+fDkffPABAMnJydx6663YbDZWr17d6gGKiEj71GgCqqmpYcyYMaxevZpbb721rWISEZELQKPPgDp27EiHDh3Yt29fW8UjIiIXiCZvwXXr1o2IiAhcLheXXHKJWf7UU0+1amAiItK+NZmAbrnlFm655Za2iEVERC4gTSaglJQUDh8+zDfffMPgwYPbIiYREbkANPke0BtvvIHT6WTcuHEAlJeXk5SUdNYHrqqqwuVy4XA4cLlcVFdX11svOzsbh8OBw+EgOzvbLC8rKyMiIoLQ0FBmzpzJyYVdZ89QgTydAAAWGUlEQVSebS6/PXbsWHbs2AGAYRjMnDmT0NBQIiMj2bhx41n3QUREzoLRhOHDhxt79+41nE6nWRYWFtbUbk1KT0835s2bZxiGYcybN8946KGH6tTZs2ePMXDgQGPPnj1GVVWVMXDgQKOqqsowDMOIjY011q1bZxw/ftwYN26c8dZbbxmGYRj79u0z93/yySeNe++91zAMw1i1apUxbtw44/jx48a6deuMq6++uskYo6Ojz7qfIiIXmuZ+dzZ5BdSpUycCAgJqlXXo0KxlhBqVl5dHSkoKcOI234oVK+rUKSgowOVyERgYSI8ePXC5XOTn5+Pz+di/fz9xcXHYbDaSk5PN/f39/c39Dx06hM1mM4+XnJyMzWYjLi6OvXv34vP5zrofIiJyZpp8BhQWFsY///lPampqcLvdPPXUU1x33XVnfeDKykpzbaGgoCAqKyvr1PF6vfTr18/cDgkJwev14vV6CQkJqVN+0iOPPEJOTg4BAQHmy7INtXUyhpMWLVrEokWLANi1a9dZ91NEROrX5KXMX//6VzZt2kTnzp2588478ff354knnmhW42PGjCE8PLzOT15eXq16NpvNvFJpCZmZmWzfvp1p06bx9NNPn9a+aWlplJaWUlpaSq9evVosJhERqa3JK6CuXbuSmZlJZmbmaTdeVFTU4Gd9+vTB5/MRHByMz+erd20hu91OcXGxue3xeBg5ciR2ux2Px1Or3G6319l/2rRpjB8/nrlz52K329m+fXuT+4iISNto8groiy++IC0tjbFjxzJ69Gjz52wlJSWZo9qys7OZOHFinToJCQkUFhZSXV1NdXU1hYWFJCQkEBwcjL+/P+vXr8cwDHJycsz93W63uX9eXh5XXXWVebycnBwMw2D9+vUEBATUuf0mIiJtp8kroNtuu4377ruP1NRUOnbs2GIHzsjI4Pbbb+eFF16gf//+vPLKK8CJZb4XLlzI888/T2BgILNnzyY2NhaAOXPmEBgYCMCCBQu4++67OXz4MImJiSQmJprtbtmyhQ4dOtC/f38WLlwIwPjx43nrrbcIDQ2la9euLFmypMX6IiIip89mGP/7Ak0DoqOjKSsra6t4zikxMTGUlpZaHYaIyHmlud+dTd6Cu/nmm1mwYAE+n4+qqirzR0RE5Gw0eQvu5HOaxx9/3Cyz2Wx89dVXrReViIi0e00moG3btrVFHCIicoFp8Bbc/Pnzzd9fffXVWp/99re/bb2IRETkgtBgAsrNzTV/nzdvXq3P8vPzWy8iERG5IDSYgE4dHPfjgXJNDJwTERFpUoMJ6NSpcX48TU5LTpsjIiIXpgYHIXzyySf4+/tjGAaHDx82Z5k2DIMjR460WYAiItI+NZiAampq2jIOERG5wJz9wj4iIiJnQAlIREQsoQQkIiKWUAISERFLKAGJiIgllIBERMQSSkAiImIJSxJQVVUVLpcLh8OBy+Wiurq63nrZ2dk4HA4cDoe5LARAWVkZERERhIaGMnPmTHNqoNmzZxMZGYnT6WTs2LHs2LEDgOLiYgICAnA6nTidTh577LHW76SIiDTKkgSUlZVFfHw8breb+Ph4srKy6tSpqqpi7ty5bNiwgZKSEubOnWsmqhkzZrB48WLcbjdut9ucHDU9PZ1PP/2U8vJyJkyYUCvRjBgxgvLycsrLy5kzZ07bdFRERBpkSQLKy8sjJSUFgJSUFFasWFGnTkFBAS6Xi8DAQHr06IHL5SI/Px+fz8f+/fuJi4vDZrORnJxs7n9yuiCAQ4cOac46EZFzWJML0rWGyspKgoODAQgKCqKysrJOHa/XS79+/cztkJAQvF4vXq+XkJCQOuUnPfLII+Tk5BAQEMDq1avN8nXr1hEVFUXfvn3505/+RFhYWL2xLVq0iEWLFgGwa9eus+uoiIg0qNWugMaMGUN4eHidn7y8vFr1bDZbi16pZGZmsn37dqZNm8bTTz8NwPDhw/n666/55JNP+MUvfsGkSZMa3D8tLY3S0lJKS0vp1atXi8UlIiK1tVoCKioq4rPPPqvzM3HiRPr06YPP5wPA5/PRu3fvOvvb7Xa2b99ubns8Hux2O3a7HY/HU6f8x6ZNm8Zrr70GnLg1161bNwDGjx/P0aNH2b17d4v2V0RETo8lz4CSkpLMUW3Z2dlMnDixTp2EhAQKCwuprq6murqawsJCEhISCA4Oxt/fn/Xr12MYBjk5Oeb+brfb3D8vL4+rrroKgJ07d5oj5UpKSjh+/Dg9e/Zs7W6KiEgjLHkGlJGRwe23384LL7xA//79eeWVVwAoLS1l4cKFPP/88wQGBjJ79mxiY2MBmDNnDoGBgQAsWLCAu+++m8OHD5OYmEhiYqLZ7pYtW+jQoQP9+/dn4cKFACxbtoxnn30WPz8/unTpQm5urgYoiIhYzGZofe0GxcTEUFpaanUYIiLnleZ+d2omBBERsYQSkIiIWEIJSERELKEEJCIillACEhERSygBiYiIJZSARETEEkpAIiJiCSUgERGxhBKQiIhYQglIREQsoQQkIiKWsGQ27AvC2xmw899WRyEicmaCIiAxq1UPoSsgERGxhK6AWksr/+UgInK+0xWQiIhYQglIREQsYVkCqqqqwuVy4XA4cLlcVFdX11svOzsbh8OBw+EgOzvbLC8rKyMiIoLQ0FBmzpzJjxd2/fOf/4zNZmP37t0AGIbBzJkzCQ0NJTIyko0bN7Ze50REpEmWJaCsrCzi4+Nxu93Ex8eTlVX3mUlVVRVz585lw4YNlJSUMHfuXDNRzZgxg8WLF+N2u3G73eTn55v7bd++ncLCQi6//HKz7O233zbrLlq0iBkzZrR+J0VEpEGWJaC8vDxSUlIASElJYcWKFXXqFBQU4HK5CAwMpEePHrhcLvLz8/H5fOzfv5+4uDhsNhvJycm19n/wwQeZP38+Nput1vGSk5Ox2WzExcWxd+9efD5f63dURETqZVkCqqysJDg4GICgoCAqKyvr1PF6vfTr18/cDgkJwev14vV6CQkJqVMOJxKN3W4nKiqqWW392KJFi4iJiSEmJoZdu3adXSdFRKRBrToMe8yYMezcubNOeWZmZq1tm81W62rlTH333Xf84Q9/oLCw8IzbSEtLIy0tDYCYmJizjklEROrXqgmoqKiowc/69OmDz+cjODgYn89H796969Sx2+0UFxeb2x6Ph5EjR2K32/F4PLXK7XY7W7duZdu2bebVj8fjYfjw4ZSUlGC329m+fXudfURExBqW3YJLSkoyR7VlZ2czceLEOnUSEhIoLCykurqa6upqCgsLSUhIIDg4GH9/f9avX49hGOTk5DBx4kQiIiL49ttvqaiooKKigpCQEDZu3EhQUBBJSUnk5ORgGAbr168nICDAvAUoIiJtz7IElJGRwTvvvIPD4aCoqIiMjAwASktLSU1NBSAwMJDZs2cTGxtLbGwsc+bMITAwEIAFCxaQmppKaGgogwYNIjExsdHjjR8/niuuuILQ0FDuueceFixY0LodFBGRRtmMH79AI6aYmBhKS0utDkNE5LzS3O9OzYQgIiKWUAISERFLKAGJiIgllIBERMQSSkAiImIJJSAREbGEVkSVC8KKj708XrCFHXsP07d7F9ITBjNpmGbCELGSEpC0eys+9vLw8n9z+GgNAN69h3l4+b8BlIRELKRbcNLuPV6wxUw+Jx0+WsPjBVssikhEQAlILgA79h4+rXIRaRtKQNLu9e3e5bTKRaRtKAFJu5eeMJgunTrWKuvSqSPpCYMtikhEQIMQ5AJwcqCBRsGJnFuUgOSCMGmYXQlH5ByjW3AiImIJJSAREbGEJQmoqqoKl8uFw+HA5XJRXV1db73s7GwcDgcOh8NcvhugrKyMiIgIQkNDmTlzJj9eU+/Pf/4zNpuN3bt3A1BcXExAQABOpxOn08ljjz3Wep0TEZFmsSQBZWVlER8fj9vtJj4+nqysrDp1qqqqmDt3Lhs2bKCkpIS5c+eaiWrGjBksXrwYt9uN2+0mPz/f3G/79u0UFhZy+eWX12pvxIgRlJeXU15ezpw5c1q3gyIi0iRLElBeXh4pKSkApKSksGLFijp1CgoKcLlcBAYG0qNHD1wuF/n5+fh8Pvbv309cXBw2m43k5ORa+z/44IPMnz8fm83WZv0REZHTZ0kCqqysJDg4GICgoCAqKyvr1PF6vfTr18/cDgkJwev14vV6CQkJqVMOJxKb3W4nKiqqTnvr1q0jKiqKxMRENm3a1GBsixYtIiYmhpiYGHbt2nXGfRQRkca12jDsMWPGsHPnzjrlmZmZtbZtNluLXK189913/OEPf6CwsLDOZ8OHD+frr7+mW7duvPXWW0yaNAm3211vO2lpaaSlpQEQExNz1nGJiEj9Wi0BFRUVNfhZnz598Pl8BAcH4/P56N27d506drud4uJic9vj8TBy5Ejsdjsej6dWud1uZ+vWrWzbts28+vF4PAwfPpySkhKCgoLM+uPHj+f+++9n9+7dXHbZZS3QUxE5E1oiQyy5BZeUlGSOasvOzmbixIl16iQkJFBYWEh1dTXV1dUUFhaSkJBAcHAw/v7+rF+/HsMwyMnJYeLEiURERPDtt99SUVFBRUUFISEhbNy4kaCgIHbu3GmOlCspKeH48eP07NmzTfssIv/n5BIZ3r2HMfi/JTJWfOy1OrQL3oqPvVyf9R4DM1ZxfdZ7rXpOLElAGRkZvPPOOzgcDoqKisjIyACgtLSU1NRUAAIDA5k9ezaxsbHExsYyZ84cAgMDAViwYAGpqamEhoYyaNAgEhMTGz3esmXLCA8PJyoqipkzZ5Kbm6tBCiIW0hIZ56a2/sPAZvz4JRoxxcTEUFpaanUYIu3OwIxV1PfFYwO2Zd3U1uHI/7o+6z289SxTYu/ehX9ljG52O8397tRMCCLS5rRExrmprdfOUgISkTanJTLOTW39h4ESkIi0uUnD7My7JQJ79y7YOHGLZ94tERoFZ7G2/sNAyzGIiCW0RMa5p63XzlICEhERU1v+YaBbcCIiYgklIBERsYQSkIiIWEIJSERELKEEJCIiltBUPI247LLLGDBgwBnvv2vXLnr16tVyAVmkvfQD1JdzUXvpB6gvJ1VUVLB79+4m6ykBtaL2Mpdce+kHqC/novbSD1BfTpduwYmIiCWUgERExBIdf//73//e6iDas+joaKtDaBHtpR+gvpyL2ks/QH05HXoGJCIiltAtOBERsYQSkIiIWEIJqAXU1NQwbNgwJkyYUOez77//njvuuIPQ0FCuueYaKioq2j7A09BYX5YuXUqvXr1wOp04nU6ef/55CyJsngEDBhAREYHT6SQmJqbO54ZhMHPmTEJDQ4mMjGTjxo0WRNk8TfWluLiYgIAA87w89thjFkTZtL179zJlyhSuuuoqhgwZwrp162p9fj6dk6b6cr6cky1btpgxOp1O/P39eeKJJ2rVac3zouUYWsCTTz7JkCFD2L9/f53PXnjhBXr06MGXX35Jbm4us2bN4uWXX7YgyuZprC8Ad9xxB08//XQbR3VmVq9ezWWXXVbvZ2+//TZutxu3282GDRuYMWMGGzZsaOMIm6+xvgCMGDGCN998sw0jOn2//OUvGTduHMuWLeOHH37gu+++q/X5+XROmuoLnB/nZPDgwZSXlwMn/vi02+1Mnjy5Vp3WPC+6AjpLHo+HVatWkZqaWu/neXl5pKSkADBlyhTeffddztVxH031pT3Jy8sjOTkZm81GXFwce/fuxefzWR1Wu7Vv3z7WrFnD9OnTAbjooovo3r17rTrnyzlpTl/OR++++y6DBg2if//+tcpb87woAZ2lX/3qV8yfP58OHer/T+n1eunXrx8Afn5+BAQEsGfPnrYMsdma6gvAa6+9RmRkJFOmTGH79u1tGN3psdlsjB07lujoaBYtWlTn81PPC0BISAher7ctQ2y2pvoCsG7dOqKiokhMTGTTpk1tHGHTtm3bRq9evfj5z3/OsGHDSE1N5dChQ7XqnC/npDl9gXP/nPxYbm4ud955Z53y1jwvSkBn4c0336R3797tYtx/c/py8803U1FRwaefforL5TKv7M5FH3zwARs3buTtt9/mmWeeYc2aNVaHdMaa6svw4cP5+uuv+eSTT/jFL37BpEmTLIq0YceOHWPjxo3MmDGDjz/+mEsuuYSsrCyrwzojzenL+XBOTvXDDz+wcuVKbrvttjY9rhLQWfjXv/7FypUrGTBgAFOnTuW9997jZz/7Wa06drvdvFI4duwY+/bto2fPnlaE26jm9KVnz5507twZgNTUVMrKyqwItVns9hNLCvfu3ZvJkydTUlJS5/NTr+A8Ho+5z7mmqb74+/vTrVs3AMaPH8/Ro0ebNRFkWwoJCSEkJIRrrrkGOHE7+scPs8+Xc9KcvpwP5+RUb7/9NsOHD6dPnz51PmvN86IEdBbmzZuHx+OhoqKC3NxcRo8ezd///vdadZKSksjOzgZg2bJljB49GpvNZkW4jWpOX06977ty5UqGDBnS1mE2y6FDhzhw4ID5e2FhIeHh4bXqJCUlkZOTg2EYrF+/noCAAIKDg60It1HN6cvOnTvN54olJSUcP378nPsjJygoiH79+rFlyxbgxPOGoUOH1qpzvpyT5vTlfDgnp3rppZfqvf0GrXteNAquFcyZM4eYmBiSkpKYPn06d911F6GhoQQGBpKbm2t1eKfl1L489dRTrFy5Ej8/PwIDA1m6dKnV4dWrsrLSHMlz7NgxfvrTnzJu3DgWLlwIwH333cf48eN56623CA0NpWvXrixZssTKkBvUnL4sW7aMZ599Fj8/P7p06UJubu45+UfOX//6V6ZNm8YPP/zAFVdcwZIlS87LcwJN9+V8OSdw4g+bd955h+eee84sa6vzoql4RETEEroFJyIillACEhERSygBiYiIJZSARETEEkpAIiJiCSUgEYs99dRTDBkyhGnTprVou8XFxeas5sXFxXz44Yct2r7I2dJ7QCIWW7BgAUVFRYSEhJhlx44dw8+v5f73LC4uplu3blx33XUt1qbI2dIVkIiF7rvvPr766isSExMJCAjgrrvu4vrrr+euu+5i6dKlTJo0CZfLxYABA3j66af5n//5H4YNG0ZcXBxVVVUAjBw5ktLSUgB2797NgAEDah2joqKChQsX8pe//AWn08natWvbupsi9VICErHQwoUL6du3L6tXr+bBBx9k8+bNFBUV8dJLLwHw2WefsXz5cj766CMeeeQRunbtyscff8y1115LTk5Os44xYMAA7rvvPh588EHKy8sZMWJEa3ZJpNmUgETOIUlJSXTp0sXcHjVqFJdeeim9evUiICCAm2++GYCIiIhzfnVdkaYoAYmcQy655JJa2ydnHwfo0KGDud2hQweOHTsGnFhn6vjx4wAcOXKkjSIVOXtKQCLnuQEDBphLYyxbtqzeOpdeeqk5q7bIuUIJSOQ895vf/IZnn32WYcOGNbjmzM0338zrr7+uQQhyTtFs2CIiYgldAYmIiCWUgERExBJKQCIiYgklIBERsYQSkIiIWEIJSERELKEEJCIilvj/2IoaJLsfGRYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "energies = [job.logfile.energy \n", " for job in frmc.queue if job.is_completed]\n", "rmults = [job.param\n", " for job in frmc.queue if job.is_completed]\n", "frmult = [rmult[0] for rmult in rmults]\n", "n_at = len(frmc.queue[0].posinp)\n", "threshold = min(energies) + n_at*frmc.precision_per_atom\n", "fig=plt.figure()\n", "fig.patch.set_facecolor('white') # When dark background\n", "plt.plot(frmult, energies,\n", " label=\"BigDFT results\", marker=\"o\", linestyle='')\n", "plt.plot(frmult, [threshold]*len(energies),\n", " label=\"Convergence threshold\")\n", "plt.xlabel(\"frmult\")\n", "plt.ylabel(\"Energy (Ha)\")\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusion\n", "\n", "The convergence with a precision per atom of 10 meV for the N$_2$ molecule is found using `hgrids = 0.36` and `rmult = [5.0, 7.0]` as input parameters for a BigDFT calculation." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.0" } }, "nbformat": 4, "nbformat_minor": 2 }