Px_IAPs_combined_CSP.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8ce686c9-8eab-4813-99c9-fb5adac4b195",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0.1+cpu\n"
     ]
    }
   ],
   "source": [
    "from pyxtal import pyxtal\n",
    "from pyxtal.lattice import Lattice\n",
    "from pymatgen.core import Structure\n",
    "from pymatgen.core import Composition\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import json\n",
    "from pyxtal.tolerance import Tol_matrix\n",
    "from datetime import datetime\n",
    "import os\n",
    "from pymatgen.analysis.structure_analyzer import SpacegroupAnalyzer\n",
    "from ase.io import read\n",
    "import matgl\n",
    "from matgl.ext.ase import M3GNetCalculator #, Relaxer\n",
    "import mace\n",
    "from mace.calculators import mace_mp\n",
    "from chgnet.model import CHGNet\n",
    "from chgnet.model import CHGNetCalculator\n",
    "\n",
    "from ase.spacegroup.symmetrize import FixSymmetry\n",
    "from ase.optimize import FIRE\n",
    "from ase.constraints import ExpCellFilter\n",
    "from ase import Atoms\n",
    "from ase.calculators.calculator import Calculator\n",
    "from ase.optimize.optimize import Optimizer\n",
    "from ase.optimize import FIRE, LBFGS, BFGSLineSearch\n",
    "from ase.io import read, write\n",
    "from ase.visualize import view\n",
    "from pymatgen.analysis import structure_matcher\n",
    "from pymatgen.io.ase import AseAtomsAdaptor\n",
    "import torch\n",
    "import re\n",
    "# from ase.visualize import view\n",
    "# import ase.optimize as opt\n",
    "print(torch.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6e3ca1d-a46a-4173-80ea-a819a8b25c7a",
   "metadata": {},
   "source": [
    "# Load IAPs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b88f4d28-d6b0-449a-8020-5b55d6fc4400",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using CHGNET IAPs \n",
      "\n",
      "CHGNet v0.3.0 initialized with 412,525 parameters\n",
      "CHGNet will run on cpu\n"
     ]
    }
   ],
   "source": [
    "# Use CHGNET\n",
    "print('Using CHGNET IAPs','\\n')\n",
    "calculator = CHGNetCalculator()\n",
    "\n",
    "\n",
    "# # Use MACE\n",
    "# print('Using MACE IAPs','\\n')\n",
    "# calculator = mace_mp(model=\"medium\", dispersion=False, default_dtype=\"float32\", device='cpu')\n",
    "\n",
    "\n",
    "\n",
    "# # Use M3GNET\n",
    "# print('Using M3GNET IAPs','\\n')\n",
    "# pot = matgl.load_model(\"M3GNet-MP-2021.2.8-PES\")\n",
    "# calculator = M3GNetCalculator(potential=pot,stress_weight=1 / 160.21766208)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b03e1793-e2c2-4217-8eb5-120082a97703",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": [
    "# Load data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1af859d-5301-496c-bb4c-f6d3e47b1dbc",
   "metadata": {},
   "source": [
    "# Modes\n",
    "#### 1. Composition only\n",
    "#### 2. += Sapce Group\n",
    "#### 3. += Wyckoff Sites\n",
    "#### 4. += volume/lattice parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "38bcd2eb-3e89-49ed-b88e-0d2710b2a7d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_formula = 'RbBiTe2'\n",
    "save_dir = 'temp_files/{}_example/'.format(test_formula)\n",
    "\n",
    "pymat_comp = Composition(test_formula)\n",
    "species = [i.value for i in pymat_comp.elements]\n",
    "test_ratio = np.array([pymat_comp.get_el_amt_dict()[i] for i in species])\n",
    "\n",
    "test_sg = 166\n",
    "# test_wyckoff = None\n",
    "test_wyckoff = [['3a'],['3b'],['6c']]\n",
    "test_volume = None\n",
    "# test_lattice = [5.6645, 10.715, 15.9032]\n",
    "test_lattice = None\n",
    "\n",
    "if test_wyckoff==None:\n",
    "    min_z = 1\n",
    "    max_z = int(40//np.sum(test_ratio))+1\n",
    "    if max_z>11:\n",
    "        max_z=11\n",
    "else:\n",
    "    tot_atom = np.sum([int(re.sub('[^0-9]','', j)) for i in test_wyckoff for j in i])\n",
    "    min_z = int(tot_atom/np.sum(test_ratio))\n",
    "    max_z = min_z+1\n",
    "\n",
    "# z=min_z\n",
    "# print(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c72c5b3c-e66b-4636-b829-c45e5b9c20ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating structures using Pyxtal for formula: RbBiTe2\n",
      "Space group: 166 | Wyckoff setting: [['3a'], ['3b'], ['6c']] | Lattic parameters: None\n",
      "Testing 'z' from 3 to 3\n"
     ]
    }
   ],
   "source": [
    "print('Generating structures using Pyxtal for formula: {}'.format(test_formula))\n",
    "print('Space group: {} | Wyckoff setting: {} | Lattic parameters: {}'.format(test_sg,test_wyckoff,test_lattice))\n",
    "print(\"Testing 'z' from {} to {}\".format(min_z, max_z-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2ae3e717-af5c-48a4-ab94-7b258b90b1b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def scale_volume(strc_in: Structure, pred_vol: float):\n",
    "#     lattice_constant = strc_in.lattice.abc[0]\n",
    "#     # factor set as the cubic root of the ratio of volume\n",
    "#     scaling_factor = (pred_vol/strc_in.volume)**(1/3)\n",
    "#     print(\"Scaling factor for abc is: \", scaling_factor)\n",
    "#     scaled_abc = lattice_constant * scaling_factor\n",
    "    scaled_structure = strc_in.copy()\n",
    "    scaled_structure.scale_lattice(pred_vol)\n",
    "    print(f\"Volume is scaled from {strc_in.volume:.4f} to {scaled_structure.volume:.4f}\")\n",
    "    return scaled_structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ccedf604-5d99-40b9-a503-af46ef88ff0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_pyxtal_structure(species,test_sg,test_ratio,z,test_wyckoff,test_lattice):\n",
    "\n",
    "    # read test input\n",
    "    test_ratio = test_ratio*z \n",
    "    test_angle = [90, 90, 90]\n",
    "\n",
    "    # temporary data processing\n",
    "    if not test_lattice==None:\n",
    "        if test_sg>=195:\n",
    "            test_lattice = [np.mean(test_lattice),np.mean(test_lattice),np.mean(test_lattice)]\n",
    "            ltype = 'cubic'\n",
    "        elif test_sg >= 143:\n",
    "            test_angle = [90, 90, 120]\n",
    "            test_lattice = [np.mean(test_lattice[:2]),np.mean(test_lattice[:2]),test_lattice[2]]\n",
    "            ltype = 'trigonal'\n",
    "            if test_sg>=168:\n",
    "                ltype = 'hexagonal'            \n",
    "        elif test_sg >=75:\n",
    "            test_lattice = [np.mean(test_lattice[:2]),np.mean(test_lattice[:2]),test_lattice[2]]\n",
    "            ltype = 'tetragonal'\n",
    "        else:\n",
    "            ltype = 'orthorhombic'\n",
    "    \n",
    "    # pyxtal crystal generation\n",
    "    inter_dis_mat = Tol_matrix(prototype=\"atomic\", factor=1.3)\n",
    "    lat_ = None\n",
    "    \n",
    "    if not test_lattice==None:\n",
    "        lat_ = Lattice.from_para(a=test_lattice[0],\n",
    "                                b=test_lattice[1],\n",
    "                                c=test_lattice[2],\n",
    "                                 alpha=test_angle[0],\n",
    "                                 beta=test_angle[1],\n",
    "                                 gamma=test_angle[2],\n",
    "                                ltype=ltype,\n",
    "                                )\n",
    "    structures = []\n",
    "    strc_pre = pyxtal()\n",
    "\n",
    "    try:\n",
    "        strc_pre.from_random(dim=3, group=test_sg, species=species,\n",
    "                                   lattice=lat_,numIons=test_ratio, sites=test_wyckoff, tm=inter_dis_mat)\n",
    "        print('Structure DoF = {}'.format(strc_pre.get_dof()))\n",
    "        strc_pre_pymatgen = strc_pre.to_pymatgen()\n",
    "        return strc_pre_pymatgen\n",
    "        \n",
    "    except Exception as e:\n",
    "        print('An exception occurred') \n",
    "        print(e.message)\n",
    "        return False\n",
    "    # strc_pre_set.get_dof()\n",
    "    # strc_pre_scaled = scale_volume(strc_pre_set.to_pymatgen(), test_volume)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1818f876-d080-4289-b709-7b420bea3cad",
   "metadata": {},
   "outputs": [],
   "source": [
    "# strc_pre_pymatgen = get_pyxtal_structure(species,test_sg,test_ratio,z,test_wyckoff,test_lattice)\n",
    "# sga = SpacegroupAnalyzer(strc_pre_pymatgen)\n",
    "# sga.get_symmetrized_structure()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f4aa5a4f-da87-47cb-8a9c-8e218d0244f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ase_relaxer(\n",
    "    atoms_in: Atoms,\n",
    "    calculator: Calculator,\n",
    "    optimizer: Optimizer = FIRE,\n",
    "    cell_filter = None,\n",
    "    fix_symmetry: bool = True,\n",
    "    fix_fractional: bool = False,\n",
    "    hydrostatic_strain: bool = False,\n",
    "    fmax: float = 0.05,\n",
    "    steps_limit: int = 500,\n",
    "    logfile: str = \"-\",\n",
    "    wdir: str = \"./\",\n",
    ") -> Atoms:\n",
    "    atoms = atoms_in.copy()\n",
    "    full_formula = atoms.get_chemical_formula(mode=\"metal\")\n",
    "    reduced_formula = atoms.get_chemical_formula(mode=\"metal\", empirical=True)\n",
    "    atoms.calc = calculator\n",
    "    E0 = atoms.get_potential_energy()\n",
    "    if fix_fractional:\n",
    "        atoms.set_constraint([FixAtoms(indices=[atom.index for atom in atoms])])\n",
    "    if fix_symmetry:\n",
    "        atoms.set_constraint([FixSymmetry(atoms)])\n",
    "    if cell_filter is not None:\n",
    "        target = cell_filter(atoms, hydrostatic_strain=hydrostatic_strain)\n",
    "    else:\n",
    "        target = atoms\n",
    "    now = datetime.now()\n",
    "    strnow = now.strftime(\"%Y%m%d%H%M%S\")\n",
    "    opt = optimizer(target,\n",
    "#                     maxstep = max_step,\n",
    "                    trajectory=f'{wdir}/{reduced_formula}_{full_formula}_{strnow}.traj',\n",
    "                    logfile=logfile,\n",
    "                   )\n",
    "    opt.run(fmax=fmax, steps=steps_limit)\n",
    "    if cell_filter is None:\n",
    "        write(filename=f'{wdir}/{reduced_formula}_{full_formula}_fix_cell_relaxed_{strnow}.cif',\n",
    "          images=atoms,\n",
    "         format=\"cif\",\n",
    "         )\n",
    "    else:\n",
    "        write(filename=f'{wdir}/{reduced_formula}_{full_formula}_relax_postitions_and_cell_{strnow}.cif',\n",
    "          images=atoms,\n",
    "         format=\"cif\",\n",
    "         )\n",
    "    cell_diff = (atoms.cell.cellpar() / atoms_in.cell.cellpar() - 1.0) * 100\n",
    "    E1 = atoms.get_potential_energy()\n",
    "    # print(\"Optimized Cell         :\", atoms.cell.cellpar())\n",
    "    # print(\"Optimized Cell diff (%):\", cell_diff)\n",
    "    # print(\"Scaled positions       :\\n\", atoms.get_scaled_positions())\n",
    "    # print(f\"Potential energy before opt: {E0:.4f} eV\")\n",
    "    # print(f\"Potential energy after  opt: {E1:.4f} eV\")\n",
    "    \n",
    "    return atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c72e1ce4-d4eb-4279-a8a8-53a9469484a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ase_iaps_relax(wdir, z,strc_pre_pymatgen,calculator):\n",
    "    # ase relaxation\n",
    "    atoms_pyxtal_scaled = AseAtomsAdaptor.get_atoms(strc_pre_pymatgen)\n",
    "    # wdir = save_dir+'z_{}'.format(z)\n",
    "    if not os.path.exists(wdir):\n",
    "        os.makedirs(wdir)\n",
    "    # fix cell relaxation\n",
    "    strc_relax_fixcell = ase_relaxer(\n",
    "        atoms_in=atoms_pyxtal_scaled,\n",
    "        calculator=calculator,\n",
    "        cell_filter=None,\n",
    "        logfile=wdir+\"/strc_volume-scaled_fix-cell_relaxed_atomic_postitions.log\",\n",
    "        wdir=wdir\n",
    "    )\n",
    "\n",
    "    # relax both cell and atoms\n",
    "    strc_relax_cell = ase_relaxer(\n",
    "        atoms_in=strc_relax_fixcell,\n",
    "        calculator=calculator,\n",
    "        cell_filter=ExpCellFilter,\n",
    "        fix_fractional=False,\n",
    "        logfile=wdir+\"/strc_volume-scaled_relax_postitions_and_cell.log\",\n",
    "        wdir=wdir\n",
    "    )\n",
    "    \n",
    "    return strc_relax_fixcell, strc_relax_cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ddd7a54f-82b5-4da4-bf87-009a35ffc886",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "62cdb99b-398e-4bd7-a639-508ebde23db7",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running z=1 for RbBiTe2\n",
      "An exception occurred\n",
      "Compoisition [1. 1. 2.] not compatible with symmetry 166\n",
      "-------------------------------------------------------------------\n",
      "Running z=2 for RbBiTe2\n",
      "An exception occurred\n",
      "Compoisition [2. 2. 4.] not compatible with symmetry 166\n",
      "-------------------------------------------------------------------\n",
      "Running z=3 for RbBiTe2\n",
      "Structure DoF = 3\n",
      "-------------------------------------------------------------------\n",
      "Running z=4 for RbBiTe2\n",
      "An exception occurred\n",
      "Compoisition [4. 4. 8.] not compatible with symmetry 166\n",
      "-------------------------------------------------------------------\n",
      "Running z=5 for RbBiTe2\n",
      "An exception occurred\n",
      "Compoisition [ 5.  5. 10.] not compatible with symmetry 166\n",
      "-------------------------------------------------------------------\n",
      "Running z=6 for RbBiTe2\n",
      "Structure DoF = 4\n",
      "-------------------------------------------------------------------\n",
      "Running z=7 for RbBiTe2\n",
      "An exception occurred\n",
      "Compoisition [ 7.  7. 14.] not compatible with symmetry 166\n",
      "-------------------------------------------------------------------\n",
      "Running z=8 for RbBiTe2\n",
      "An exception occurred\n",
      "Compoisition [ 8.  8. 16.] not compatible with symmetry 166\n",
      "-------------------------------------------------------------------\n",
      "Running z=9 for RbBiTe2\n",
      "Structure DoF = 4\n",
      "-------------------------------------------------------------------\n",
      "Running z=10 for RbBiTe2\n",
      "An exception occurred\n",
      "Compoisition [10. 10. 20.] not compatible with symmetry 166\n",
      "-------------------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "structures_result = []\n",
    "structures = []\n",
    "lat_volume = []\n",
    "all_predicted_volume = []\n",
    "all_volume = [] \n",
    "all_energy = [] \n",
    "ind_lst = []\n",
    "z_lst = []\n",
    "all_strc_relaxed = {}\n",
    "# all_vasp_volume = []\n",
    "# all_Dof = []\n",
    "# min_z = 10\n",
    "# max_z = 12\n",
    "\n",
    "for i in range(min_z, max_z):\n",
    "    z = int(i)\n",
    "    print('Running z={} for {}'.format(z,test_formula))\n",
    "\n",
    "    for ind in range(8):\n",
    "        # generate from pyxtal\n",
    "        strc_pre_pymatgen = get_pyxtal_structure(species,test_sg,test_ratio,z,test_wyckoff,test_lattice)\n",
    "        # if not test_volume==None:\n",
    "        #     strc_pre_pymatgen = scale_volume(strc_pre_pymatgen, test_volume)\n",
    "        \n",
    "        if not strc_pre_pymatgen:\n",
    "            print('-------------------------------------------------------------------')\n",
    "            break\n",
    "            \n",
    "        # IAPs relax\n",
    "        print('Relaxing structures using IPAs')\n",
    "        wdir = save_dir+'z_{}_{}'.format(z,ind)\n",
    "        strc_relax_fixcell, strc_relax_cell = ase_iaps_relax(wdir,z,strc_pre_pymatgen,calculator)\n",
    "        \n",
    "        # compare IAPs relaxed and DFT relaxed structure\n",
    "        ind_lst.append('z_{}_{}'.format(z,ind))\n",
    "        z_lst.append(z)\n",
    "        all_energy.append(strc_relax_cell.get_potential_energy()/len(strc_relax_cell))\n",
    "        strc_relaxed = AseAtomsAdaptor.get_structure(strc_relax_cell)\n",
    "        strc_relaxed = SpacegroupAnalyzer(structure=strc_relaxed #, symprec=0.1\n",
    "                      ).get_conventional_standard_structure()\n",
    "        all_volume.append(strc_relaxed.volume)\n",
    "        # reduced_formula = strc_relax_cell.get_chemical_formula(mode=\"metal\", empirical=True)\n",
    "        all_strc_relaxed['z_{}_{}'.format(z,ind)] = strc_relaxed\n",
    "    \n",
    "    if strc_pre_pymatgen:\n",
    "        print('-------------------------------------------------------------------')\n",
    "        continue\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7c47e4b4-605d-4be4-8d80-79b640d886d3",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ind</th>\n",
       "      <th>total_energy_per_atom</th>\n",
       "      <th>z_value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>z_3_0</td>\n",
       "      <td>-3.372469</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>z_9_0</td>\n",
       "      <td>-3.306879</td>\n",
       "      <td>9</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>z_6_0</td>\n",
       "      <td>-2.705953</td>\n",
       "      <td>6</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     ind  total_energy_per_atom  z_value\n",
       "0  z_3_0              -3.372469        3\n",
       "2  z_9_0              -3.306879        9\n",
       "1  z_6_0              -2.705953        6"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_result = pd.DataFrame({'ind':ind_lst,'total_energy_per_atom':all_energy,'z_value':z_lst})\n",
    "df_result.sort_values(by='total_energy_per_atom')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "b5c352d4-1ac5-4a4a-ad5b-5556d2bba0d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a8d36e9a-ba47-4a40-8ce6-8c699883373e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='z_value', ylabel='total_energy_per_atom'>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAGxCAYAAABhi7IUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA13UlEQVR4nO3de1yUdf7//+eIcpDDqIABCQpoopZmbq3QSpJW2kHtaC5pHr/ubgezNLEyKy2ybLdNP1luLmZ5WEtbtcOqrWUZ5Skt3YQVLFGQPBAMig7IzO+Pbs4vVsFhnGFmuh732+263Zjrel/XvGZuN5tn7/f7el8mu91uFwAAgEE183YBAAAA3kQYAgAAhkYYAgAAhkYYAgAAhkYYAgAAhkYYAgAAhkYYAgAAhkYYAgAAhtbc2wX4OpvNppKSEoWHh8tkMnm7HAAA4AS73a7KykrFxcWpWbOG+34IQ+dRUlKi+Ph4b5cBAABccODAAbVr167BNoSh8wgPD5f085cZERHh5WoAAIAzLBaL4uPjHb/jDSEMnceZobGIiAjCEAAAfsaZKS5MoAYAAIZGGAIAAIZGGAIAAIZGGAIAAIZGGAIAAIZGGAIAAIbmF2Hohx9+0JgxY5SYmKiQkBAlJydr+vTpqq6ubvAck8l0zu2dd95pwuoBAIAv84t1hvLy8mSz2fT666+rY8eO2r17t8aNG6cTJ05o9uzZ5zwnPj5ehw4dqrNv/vz5evHFFzVw4MCmKBsAAPgBk91ut3u7CFe8+OKLmjdvnvbt2+f0OT179tQVV1yhBQsWOH2OxWKR2WxWRUUFiy4CAOAnGvP77RfDZOdSUVGhNm3aON1++/bt2rlzp8aMGePBqgAAgL/xi2Gy/1VQUKA5c+bUO0R2LgsWLFCXLl2UlpbWYDur1Sqr1ep4bbFYXK4TgO+qqKrW0ePVspyqUURIC0WFBsrcMtDbZQHwAq/2DGVlZdU7yfnMlpeXV+ec4uJiDRgwQHfeeafGjRvn1PucPHlSS5YscapXKDs7W2az2bHxxHrg16ek/KTuX7pD/f68Ube+mqt+L23UA0t3qKT8pLdLA+AFXp0zdOTIER07dqzBNklJSQoM/Pn/1kpKStS3b1/17t1bCxcuVLNmzmW5t956S2PGjFFxcbGio6MbbHuunqH4+HjmDAG/EhVV1bp/6Q59vvfoWcfSO0VpzrCe9BABvwKNmTPk1WGy6Ojo84aTM4qLi5WRkaFevXopJyfH6SAk/TxENmjQIKfeKygoSEFBQU5fG4B/OXq8+pxBSJI+23tUR49XE4YAg/GLCdTFxcXq27evEhISNHv2bB05ckSlpaUqLS2t0yYlJUVbtmypc25BQYE+++wzjR07tqnLBuCDLKdqGjxeeZ7jAH59/GIC9fr161VQUKCCggK1a9euzrEzo3w1NTXKz89XVVVVneN///vf1a5dO11//fVNVi8A3xUR3KLB4+HnOQ7g18dv1xlqKqwzBPy6VFRV64GlO/QZc4aAXzVDrDMEAK4wtwzU87d3V3qnqDr70ztFadbt3QlCgAH5xTAZALhTXKsQzRnWU0ePV6vyVI3Cg1soKox1hgCjIgwBMCRzS8IPgJ8xTAYAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAyNMAQAAAzNL8LQDz/8oDFjxigxMVEhISFKTk7W9OnTVV1d3eB5paWlGj58uGJiYhQaGqorrrhCK1asaKKqAQCAP2ju7QKckZeXJ5vNptdff10dO3bU7t27NW7cOJ04cUKzZ8+u97wRI0aovLxcq1evVlRUlJYsWaK77rpL27ZtU8+ePZvwEwAAAF9lstvtdm8X4YoXX3xR8+bN0759++ptExYWpnnz5mn48OGOfZGRkZo1a5bGjh3r1PtYLBaZzWZVVFQoIiLigusGAACe15jfb78YJjuXiooKtWnTpsE2aWlp+sc//qGysjLZbDYtW7ZMp06dUt++fes9x2q1ymKx1NkAAMCvl1+GoYKCAs2ZM0fjx49vsN3y5ctVU1OjyMhIBQUFafz48XrvvffUsWPHes/Jzs6W2Wx2bPHx8e4uHwAA+BCvhqGsrCyZTKYGt7y8vDrnFBcXa8CAAbrzzjs1bty4Bq8/bdo0lZeX6+OPP9a2bdv08MMP66677tKuXbvqPWfq1KmqqKhwbAcOHHDLZwUAAL7Jq3OGjhw5omPHjjXYJikpSYGBgZKkkpIS9e3bV71799bChQvVrFn9Wa6wsNAx2bpbt26O/f3791fHjh312muvOVUjc4YAAPA/jfn99urdZNHR0YqOjnaqbXFxsTIyMtSrVy/l5OQ0GIQkqaqqSpLOahcQECCbzeZawQAA4FfHL+YMFRcXq2/fvkpISNDs2bN15MgRlZaWqrS0tE6blJQUbdmyRZKUkpKijh07avz48dqyZYsKCwv10ksvaf369RoyZIiXPgkAAPA1frHO0Pr161VQUKCCggK1a9euzrEzo3w1NTXKz8939Ai1aNFCH374obKysnTLLbfo+PHj6tixo958803deOONTf4ZAACAb/LbdYaaCnOGAADwP4ZYZwgAAMAdCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQCEMAAMDQmrt6YklJiTZt2qTDhw/LZrPVOfbggw9ecGEAAABNwaUwtHDhQo0fP16BgYGKjIyUyWRyHDOZTIQhAADgN0x2u93e2JPi4+P1hz/8QVOnTlWzZr/ukTaLxSKz2ayKigpFRER4uxwAAOCExvx+u5RkqqqqdPfddzdZEPrhhx80ZswYJSYmKiQkRMnJyZo+fbqqq6sbPK+wsFC33nqroqOjFRERobvuuks//vhjk9QMAAD8g0tpZsyYMXrnnXfcXUu98vLyZLPZ9Prrr+s///mP/vKXv+i1117TY489Vu85J06c0PXXXy+TyaQNGzboiy++UHV1tW655Zaz5jgBAADjcmmYrLa2VjfffLNOnjypyy67TC1atKhz/M9//rPbCqzPiy++qHnz5mnfvn3nPL5u3ToNHDhQP/30k6N7rKKiQq1bt9a6devUv39/p96HYTIAAPxPY36/XZpAnZ2drbVr16pz586SdNYE6qZQUVGhNm3a1HvcarXKZDIpKCjIsS84OFjNmjXTpk2bnA5DAADg182lMPTSSy/p73//u0aOHOnmcpxTUFCgOXPmaPbs2fW26d27t0JDQzVlyhQ999xzstvtysrKUm1trQ4dOlTveVarVVar1fHaYrG4tXYAAOBbXJozFBQUpKuvvvqC3zwrK0smk6nBLS8vr845xcXFGjBggO68806NGzeu3mtHR0frnXfe0Zo1axQWFiaz2azy8nJdccUVDU78zs7Oltlsdmzx8fEX/DkBAIDvcmnOUHZ2tg4dOqRXXnnlgt78yJEjOnbsWINtkpKSFBgYKOnnhR779u2r3r17a+HChU7fzXb06FE1b95crVq1UkxMjB555BFNnjz5nG3P1TMUHx/PnCEAAPyIx+cMbdmyRRs2bND777+vbt26nTWBeuXKlU5dJzo6WtHR0U61LS4uVkZGhnr16qWcnJxG3dYfFRUlSdqwYYMOHz6sQYMG1ds2KCiozjwjAADw6+ZSGGrVqpVuu+02d9dSr+LiYvXt21ft27fX7NmzdeTIEcexmJgYR5t+/fpp0aJFuuqqqyRJOTk56tKli6Kjo/Xll19qwoQJmjhxomPiNwAAgEthKCcnx911NGj9+vUqKChQQUGB2rVrV+fYmVG+mpoa5efnq6qqynEsPz9fU6dOVVlZmTp06KDHH39cEydObNLaAQCAb3NpztAZR44cUX5+viSpc+fOTg95+RPWGQIAwP94/HEcJ06c0OjRoxUbG6v09HSlp6crLi5OY8aMqdMzAwAA4OtcCkMPP/ywNm7cqDVr1qi8vFzl5eVatWqVNm7cqEceecTdNQIAAHiMS8NkUVFRevfdd9W3b986+z/55BPddddddSY4+zuGyQAA8D9N8tT6iy666Kz9bdu2ZZgMAAD4FZfCUGpqqqZPn65Tp0459p08eVJPP/20UlNT3VYcAACAp7l0a/3LL7+sAQMGqF27durRo4ck6ZtvvlFwcLDWrl3r1gIBAAA8yeVb66uqqrR48WLHs8O6dOmizMxMhYSEuLVAb2POEAAA/sfjj+P47LPPlJaWdtaDUk+fPq3PPvtM6enprlwWAACgybk0ZygjI0NlZWVn7a+oqFBGRsYFFwUAANBUXApDdrtdJpPprP3Hjh1TaGjoBRcFAADQVBo1THbm4awmk0kjR46s83T32tpaffvtt0pLS3NvhQAAAB7UqDBkNpsl/dwzFB4eXmeydGBgoHr37n3WPCIAAABf1qgwdOZp9R06dNCkSZMYEgMAAH7vgp5abwTcWg8AgP/x+K31kvTuu+9q+fLlKioqUnV1dZ1jX3/9tauXBQAAaFIu3U32yiuvaNSoUbrooou0Y8cOXXXVVYqMjNS+ffs0cOBAd9cIAADgMS6FoVdffVXz58/XnDlzFBgYqEcffVTr16/Xgw8+qIqKCnfXCAAA4DEuhaGioiLHLfQhISGqrKyUJA0fPlxLly51X3UAAAAe5lIYiomJcaxAnZCQoK+++kqS9P3334v52AAAwJ+4FIauvfZarV69WpI0atQoTZw4Udddd52GDh2qW2+91a0FAgAAeJJLt9bbbDbZbDY1b/7zzWjLli1Tbm6uOnXqpPHjxyswMNDthXoLt9YDAOB/GvP77dF1hv70pz/pmWeeUVRUlKfewuMIQwAA+J/G/H67NEzmrLffflsWi8WTbwEAAHBBPBqGmEwNAAB8nUfDEAAAgK8jDAEAAEMjDAEAAEMjDAEAAENrdBg6ffq0nnnmGR08ePC8be+55x5uRwcAAD7NpXWGwsPDtWvXLnXo0MEDJfkW1hkCAMD/eHydoWuvvVYbN250qTgAAABf0tyVkwYOHKisrCzt2rVLvXr1UmhoaJ3jgwYNcktxAAAAnubSMFmzZvV3KJlMJtXW1l5QUb6EYTIAAPxPY36/XeoZstlsLhUGAADgay741vpTp065ow4AAACvcCkM1dbWasaMGbr44osVFhamffv2SZKmTZumBQsWuLVAAAAAT3IpDD377LNauHChXnjhBQUGBjr2X3rppXrjjTfcVhwAAICnuRSGFi1apPnz5yszM1MBAQGO/T169FBeXp7bigMAAPA0l8JQcXGxOnbseNZ+m82mmpqaCy4KAACgqbgUhrp27arPP//8rP3vvvuuevbsecFFAQAANBWXwtCTTz6p+++/X7NmzZLNZtPKlSs1btw4Pfvss3ryySfdXaOknxdyTEhIUHBwsGJjYzV8+HCVlJQ0eM6pU6d03333KTIyUmFhYbr99tv1448/eqQ+AADgn1wKQ4MHD9aaNWv08ccfKzQ0VE8++aT27NmjNWvW6LrrrnN3jZKkjIwMLV++XPn5+VqxYoUKCwt1xx13NHjOxIkTtWbNGr3zzjvauHGjSkpKdNttt3mkPgAA4J9cWoHaF6xevVpDhgyR1WpVixYtzjpeUVGh6OhoLVmyxBGa8vLy1KVLF3355Zfq3bu3U+/DCtQAAPgfj69Afca2bdu0Z88eST/PI+rVq9eFXM5pZWVlWrx4sdLS0s4ZhCRp+/btqqmpUf/+/R37UlJSlJCQ0GAYslqtslqtjtcWi8W9xQMAAJ/i0jDZwYMH1adPH1111VWaMGGCJkyYoCuvvFK/+93vdPDgQXfX6DBlyhSFhoYqMjJSRUVFWrVqVb1tS0tLFRgYqFatWtXZf9FFF6m0tLTe87Kzs2U2mx1bfHy8u8oHAAA+yKUwNHbsWNXU1GjPnj0qKytTWVmZ9uzZI5vNprFjxzp9naysLJlMpga3X65bNHnyZO3YsUPr1q1TQECARowYIXeP8k2dOlUVFRWO7cCBA269PgAA8C0uDZNt3LhRubm56ty5s2Nf586dNWfOHPXp08fp6zzyyCMaOXJkg22SkpIcf0dFRSkqKkqXXHKJunTpovj4eH311VdKTU0967yYmBhVV1ervLy8Tu/Qjz/+qJiYmHrfLygoSEFBQU5/BgAA4N9cCkPx8fHnXFyxtrZWcXFxTl8nOjpa0dHRrpQgm80mSXXm9/xSr1691KJFC/373//W7bffLknKz89XUVHROcMTAAAwJpeGyV588UU98MAD2rZtm2Pftm3bNGHCBM2ePdttxZ2xefNmzZ07Vzt37tT+/fu1YcMGDRs2TMnJyY5gU1xcrJSUFG3ZskWSZDabNWbMGD388MP65JNPtH37do0aNUqpqalO30kGAAB+/Vy6tb5169aqqqrS6dOn1bz5z51LZ/4ODQ2t07asrOyCi9y1a5cmTJigb775RidOnFBsbKwGDBigJ554QhdffLEk6YcfflBiYqI++eQT9e3bV9LPiy4+8sgjWrp0qaxWq2644Qa9+uqrDQ6T/S9urQcAwP805vfbpTD05ptvOt323nvvbezlfQphCAAA/+PxdYacDTjPP//8WROYAQAAfIlLc4ac9dxzz7llmAwAAMBTPBqG/PRJHwAAwEA8GoYAAAB8HWEIAAAYGmEIAAAYGmEIAAAYmkfDUJ8+fRQSEuLJtwAAALggLoWha665RosWLdLJkycbbPfhhx8qNjbWpcIAAACagkthqGfPnpo0aZJiYmI0btw4ffXVV+6uCwAAoEm4FIZefvlllZSUKCcnR4cPH1Z6erq6du2q2bNn68cff3R3jQAAAB7j8pyh5s2b67bbbtOqVat08OBB/f73v9e0adMUHx+vIUOGaMOGDe6sEwAAwCMueAL1li1bNH36dL300ktq27atpk6dqqioKN18882aNGmSO2oEAADwGJeeWn/48GG99dZbysnJ0d69e3XLLbdo7NixuuGGG2QymSRJmzZt0oABA3T8+HG3F92UeGo9AAD+x+NPrW/Xrp2Sk5M1evRojRw5UtHR0We16d69u6688kpXLg8AANBkXApD//73v9WnT58G20REROiTTz5xqSgAAICm4tKcofMFIQAAAH/hUs9Qz549HXODfslkMik4OFgdO3bUyJEjlZGRccEFAgAAeJJLPUMDBgzQvn37FBoaqoyMDGVkZCgsLEyFhYW68sordejQIfXv31+rVq1yd70AAABu5VLP0NGjR/XII49o2rRpdfbPnDlT+/fv17p16zR9+nTNmDFDgwcPdkuhAAAAnuDSrfVms1nbt29Xx44d6+wvKChQr169VFFRoby8PF155ZWqrKx0W7HewK31AAD4n8b8frs0TBYcHKzc3Nyz9ufm5io4OFiSZLPZHH8DAAD4KpeGyR544AH94Q9/0Pbt2x1rCW3dulVvvPGGHnvsMUnS2rVrdfnll7utUAAAAE9waZhMkhYvXqy5c+cqPz9fktS5c2c98MAD+v3vfy9JOnnypOPuMn/GMBkAAP7HoytQnz59Ws8995xGjx6tzMzMetuFhIQ09tIAAABNrtFzhpo3b64XXnhBp0+f9kQ9AAAATcqlCdT9+vXTxo0b3V0LAABAk3NpAvXAgQOVlZWlXbt2qVevXgoNDa1zfNCgQW4pDgAAwNNcmkDdrFn9HUomk0m1tbUXVJQvYQI1AAD+x6MTqKWf1xACAAD4NXBpztAvnTp1yh11AAAAeIVLYai2tlYzZszQxRdfrLCwMO3bt0+SNG3aNC1YsMCtBQIAAHiSS2Ho2Wef1cKFC/XCCy8oMDDQsf/SSy/VG2+84bbiAAAAPM2lMLRo0SLNnz9fmZmZCggIcOzv0aOH8vLy3FYcAACAp7kUhoqLi896Yr3088TqmpqaCy4KAACgqbgUhrp27arPP//8rP3vvvuuevbsecFFAQAANBWXbq1/8sknde+996q4uFg2m00rV65Ufn6+Fi1apPfff9/dNQIAAHiMSz1DgwcP1po1a/Txxx8rNDRUTz75pPbs2aM1a9bouuuuc3eNAAAAHuPSCtRGwgrUAAD4n8b8fl/QoovV1dU6ePCgioqK6myeMGjQICUkJCg4OFixsbEaPny4SkpKGjxn/vz56tu3ryIiImQymVReXu6R2gAAgP9yKQzt3btXffr0UUhIiNq3b6/ExEQlJiaqQ4cOSkxMdHeNkqSMjAwtX75c+fn5WrFihQoLC3XHHXc0eE5VVZUGDBigxx57zCM1AQAA/+fSMNnVV1+t5s2bKysrS7GxsTKZTHWO9+jRw20F1mf16tUaMmSIrFarWrRo0WDbTz/9VBkZGfrpp5/UqlWrRr0Pw2QAAPgfjz+odefOndq+fbtSUlJcKvBClZWVafHixUpLSztvEGosq9Uqq9XqeG2xWNx6fQAA4FtcXmfo6NGj7q7lvKZMmaLQ0FBFRkaqqKhIq1atcvt7ZGdny2w2O7b4+Hi3vwcAAPAdLoWhWbNm6dFHH9Wnn36qY8eOyWKx1NmclZWVJZPJ1OD2y8d7TJ48WTt27NC6desUEBCgESNGyN03w02dOlUVFRWO7cCBA269PgAA8C0uzRlq1uznDPW/c4XsdrtMJpNqa2udus6RI0d07NixBtskJSXVeRjsGQcPHlR8fLxyc3OVmpra4DWYMwQAgLF4fM7QJ5984lJh/ys6OlrR0dEunWuz2SSpzvweAACAxnIpDF1zzTXurqNBmzdv1tatW/W73/1OrVu3VmFhoaZNm6bk5GRHr1BxcbH69eunRYsW6aqrrpIklZaWqrS0VAUFBZKkXbt2KTw8XAkJCWrTpk2TfgYAAOCbXF508fPPP9c999yjtLQ0FRcXS5Leeustbdq0yW3FndGyZUutXLlS/fr1U+fOnTVmzBh1795dGzduVFBQkCSppqZG+fn5qqqqcpz32muvqWfPnho3bpwkKT09XT179tTq1avdXiMAAPBPLs0ZWrFihYYPH67MzEy99dZb+u6775SUlKS5c+fqww8/1IcffuiJWr2COUMAAPgfjz+OY+bMmXrttdf0t7/9rc46P1dffbW+/vprVy4JAADgFS6Fofz8fKWnp5+132w28/wvAADgV1wKQzExMY5Jyb+0adMmJSUlXXBRAAAATcWlMDRu3DhNmDBBmzdvlslkUklJiRYvXqxJkybpj3/8o7trBAAA8BiXbq3PysqSzWZTv379VFVVpfT0dAUFBWnSpEl64IEH3F0jAACAx7h0N9kZ1dXVKigo0PHjx9W1a1eFhYXVOX7w4EHFxcU5Vqz2R9xNBgCA//H4CtRnBAYGqmvXrvUe79q1q3bu3Mk8IgAA4LM82mXj7oeoAgAAuJv/jl8BAAC4AWEIAAAYGmEIAAAYmkfDkMlk8uTlAQAALhgTqAEAgKFd0K315/Pdd98pLi7Ok28BAABwQZwOQ7fddpvTF125cqUkKT4+vvEVAQAANCGnw5DZbPZkHQAAAF7hdBjKycnxZB0AAABewa31AADA0FyeQP3uu+9q+fLlKioqUnV1dZ1jX3/99QUXBgAA0BRc6hl65ZVXNGrUKF100UXasWOHrrrqKkVGRmrfvn0aOHCgu2sEAADwGJfC0Kuvvqr58+drzpw5CgwM1KOPPqr169frwQcfVEVFhbtrBAAA8BiXwlBRUZHS0tIkSSEhIaqsrJQkDR8+XEuXLnVfdQAAAB7mUhiKiYlRWVmZJCkhIUFfffWVJOn7779n1WkAAOBXXApD1157rVavXi1JGjVqlCZOnKjrrrtOQ4cO1a233urWAgEAADzJZHehK8dms8lms6l5859vRlu2bJlyc3PVqVMnjR8/XoGBgW4v1FssFovMZrMqKioUERHh7XIAAIATGvP77VIYKioqUnx8/FlPpbfb7Tpw4IASEhIae0mfRRgCAMD/NOb326VhssTERB05cuSs/WVlZUpMTHTlkgAAAF7hUhiy2+1n9QpJ0vHjxxUcHHzBRQEAADSVRq1A/fDDD0uSTCaTpk2bppYtWzqO1dbWavPmzbr88svdWiAAAIAnNSoM7dixQ9LPPUO7du2qM1E6MDBQPXr00KRJk9xbIQAAgAc1Kgx98sknkn6+nf6vf/0rE4oBAIDfc+lBrTk5OY6/Dx48KElq166deyoCAABoQi5NoLbZbHrmmWdkNpvVvn17tW/fXq1atdKMGTNks9ncXSMAAIDHuNQz9Pjjj2vBggV6/vnndfXVV0uSNm3apKeeekqnTp3Ss88+69YiAQAAPMWlRRfj4uL02muvadCgQXX2r1q1Sn/6059UXFzstgK9jUUXAQDwPx5fdLGsrEwpKSln7U9JSXE8wBUAAMAfuBSGevTooblz5561f+7cuerRo8cFFwUAANBUXJoz9MILL+imm27Sxx9/rNTUVEnSl19+qQMHDujDDz90a4EAAACe5PKzyf773//q1ltvVXl5ucrLy3XbbbcpPz9f7du3d3eNAAAAHuPSBOqAgAAdOnRIbdu2rbP/2LFjatu2rWpra91W4BmDBg3Szp07dfjwYbVu3Vr9+/fXrFmzFBcXd872ZWVlmj59utatW6eioiJFR0dryJAhmjFjhsxms9PvywRqAAD8j8cnUNeXnzz5oNaMjAwtX75c+fn5WrFihQoLC3XHHXfU276kpEQlJSWaPXu2du/erYULF+pf//qXxowZ45H6AACAf2pUz9CZB7X+9a9/1bhx4875oNaAgAB98cUX7q/0f6xevVpDhgyR1WpVixYtnDrnnXfe0T333KMTJ06oeXPnpkvRMwQAgP9pzO+3Xz6otaysTIsXL1ZaWprTQUiS4wtxNggBAIBfP796UOuUKVM0d+5cVVVVqXfv3nr//fedPvfo0aOaMWOG/t//+38NtrNarbJarY7XFovF5XoBAIDvc2nOUE5OjluCUFZWlkwmU4NbXl6eo/3kyZO1Y8cOrVu3TgEBARoxYkS985d+yWKx6KabblLXrl311FNPNdg2OztbZrPZscXHx1/oxwQAAD7MpbvJ3OXIkSM6duxYg22SkpLqDMedcfDgQcXHxys3N9ex1tG5VFZW6oYbblDLli31/vvvn3eC97l6huLj45kzBACAH/HYnCF3i46OVnR0tEvn2mw2SaoTXP6XxWLRDTfcoKCgIK1evdqpO92CgoIUFBTkUk0AAMD/uDRM1tQ2b96suXPnaufOndq/f782bNigYcOGKTk52dErVFxcrJSUFG3ZskXSz0Ho+uuv14kTJ7RgwQJZLBaVlpaqtLTUI+sgAQAA/+QXt1W1bNlSK1eu1PTp03XixAnFxsZqwIABeuKJJxy9ODU1NcrPz1dVVZUk6euvv9bmzZslSR07dqxzve+//14dOnRo0s8AAAB8k1fnDPkD1hkCAMD/eHwFagAAgF8LwhAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0whAAADA0vwlDgwYNUkJCgoKDgxUbG6vhw4erpKSkwXPGjx+v5ORkhYSEKDo6WoMHD1ZeXl4TVQwAAPyB34ShjIwMLV++XPn5+VqxYoUKCwt1xx13NHhOr169lJOToz179mjt2rWy2+26/vrrVVtb20RVAwAAX2ey2+12bxfhitWrV2vIkCGyWq1q0aKFU+d8++236tGjhwoKCpScnOzUORaLRWazWRUVFYqIiLiQkgEAQBNpzO938yaqya3Kysq0ePFipaWlOR2ETpw4oZycHCUmJio+Pr7edlarVVar1fHaYrFccL0AAMB3+c0wmSRNmTJFoaGhioyMVFFRkVatWnXec1599VWFhYUpLCxMH330kdavX6/AwMB622dnZ8tsNju2hoITAABwTUVVtQoPH9eOop9UeOS4KqqqvVaLV4fJsrKyNGvWrAbb7NmzRykpKZKko0ePqqysTPv379fTTz8ts9ms999/XyaTqd7zKyoqdPjwYR06dEizZ89WcXGxvvjiCwUHB5+z/bl6huLj4xkmAwDATUrKT2rKim/1+d6jjn3pnaL0/O3dFdcqxC3v0ZhhMq+GoSNHjujYsWMNtklKSjpnT87BgwcVHx+v3NxcpaamOvV+1dXVat26td544w0NGzbMqXOYMwQAgPtUVFXr/qU76gShM9I7RWnOsJ4yt6x/BMdZfjNnKDo6WtHR0S6da7PZJKlOL8752O122e32Rp0DAADc5+jx6nMGIUn6bO9RHT1e7ZYw1Bh+MWdo8+bNmjt3rnbu3Kn9+/drw4YNGjZsmJKTkx29QsXFxUpJSdGWLVskSfv27VN2dra2b9+uoqIi5ebm6s4771RISIhuvPFGb34cAAAMy3KqpsHjlec57gl+EYZatmyplStXql+/furcubPGjBmj7t27a+PGjQoKCpIk1dTUKD8/X1VVVZKk4OBgff7557rxxhvVsWNHDR06VOHh4crNzVXbtm29+XEAADCsiOCG7wIPP89xT/DbdYaaCnOGAABwn4qqaj2wdIc+86E5Q37RMwQAAH4dzC0D9fzt3ZXeKarO/vROUZp1e/cmny8k+emiiwAAwH/FtQrRnGE9dfR4tSpP1Sg8uIWiwgK9EoQkwhAAAPACc0vvhZ//xTAZAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwNMIQAAAwtObeLsCoKqqqdfR4tSynahQR0kJRoYEytwz0dlkAABgOYcgLSspPasqKb/X53qOOfemdovT87d0V1yrEi5UBAGA8DJM1sYqq6rOCkCR9tveoslZ8q4qqai9VBgCAMRGGmtjR49VnBaEzPtt7VEePE4YAAGhKhKEmZjlV0+DxyvMcBwAA7kUYamIRwS0aPB5+nuMAAMC9CENNLCosUOmdos55LL1TlKLCuKMMAICmRBhqYuaWgXr+9u5nBaL0TlGadXt3bq8HAKCJcWu9F8S1CtGcYT119Hi1Kk/VKDy4haLCWGcIAABvIAx5ibkl4QcAAF/AMBkAADA0whAAADA0whAAADA0vwlDgwYNUkJCgoKDgxUbG6vhw4erpKTEqXPtdrsGDhwok8mkf/7zn54tFAAA+BW/CUMZGRlavny58vPztWLFChUWFuqOO+5w6tyXX35ZJpPJwxUCAAB/5Dd3k02cONHxd/v27ZWVlaUhQ4aopqZGLVrUv2rzzp079dJLL2nbtm2KjY1tilIBAIAf8ZueoV8qKyvT4sWLlZaW1mAQqqqq0u9//3v93//9n2JiYpqwQgAA4C/8KgxNmTJFoaGhioyMVFFRkVatWtVg+4kTJyotLU2DBw92+j2sVqssFkudDQAA/Hp5NQxlZWXJZDI1uOXl5TnaT548WTt27NC6desUEBCgESNGyG63n/Paq1ev1oYNG/Tyyy83qqbs7GyZzWbHFh8ffyEfEQAA+DiTvb400QSOHDmiY8eONdgmKSlJgYFnr9R88OBBxcfHKzc3V6mpqWcdf+ihh/TKK6+oWbP/P+/V1taqWbNm6tOnjz799NNzvp/VapXVanW8tlgsio+PV0VFhSIiIpz8ZAAAwJssFovMZrNTv99enUAdHR2t6Ohol8612WySVCe4/FJWVpbGjh1bZ99ll12mv/zlL7rlllvqvW5QUJCCgoIcr89kRYbLAADwH2d+t53p8/GLu8k2b96srVu36ne/+51at26twsJCTZs2TcnJyY5eoeLiYvXr10+LFi3SVVddpZiYmHNOmk5ISFBiYqLT711ZWSlJDJcBAOCHKisrZTabG2zjF2GoZcuWWrlypaZPn64TJ04oNjZWAwYM0BNPPOHoxampqVF+fr6qqqrc+t5xcXE6cOCAwsPD3b5W0ZkhuAMHDjAEdx58V87ju3Ie35Xz+K6cx3flPE9+V3a7XZWVlYqLiztvW78IQ5dddpk2bNjQYJsOHTqctyvMlelRzZo1U7t27Rp9XmNERETwD8ZJfFfO47tyHt+V8/iunMd35TxPfVfn6xE6w69urQcAAHA3whAAADA0wpAXBQUFafr06XXuXsO58V05j+/KeXxXzuO7ch7flfN85bvy6jpDAAAA3kbPEAAAMDTCEAAAMDTCEAAAMDTCUBObN2+eunfv7lhTITU1VR999JG3y/ILzz//vEwmkx566CFvl+KTnnrqqbMedJySkuLtsnxWcXGx7rnnHkVGRiokJESXXXaZtm3b5u2yfE6HDh3O+RDt++67z9ul+Zza2lpNmzZNiYmJCgkJUXJysmbMmOHSGndGUFlZqYceekjt27dXSEiI0tLStHXrVq/U4heLLv6atGvXTs8//7w6deoku92uN998U4MHD9aOHTvUrVs3b5fns7Zu3arXX39d3bt393YpPq1bt276+OOPHa+bN+ef+Ln89NNPuvrqq5WRkaGPPvpI0dHR2rt3r1q3bu3t0nzO1q1bVVtb63i9e/duXXfddbrzzju9WJVvmjVrlubNm6c333xT3bp107Zt2zRq1CiZzWY9+OCD3i7P54wdO1a7d+/WW2+9pbi4OL399tvq37+/vvvuO1188cVNWgt3k/mANm3a6MUXX9SYMWO8XYpPOn78uK644gq9+uqrmjlzpi6//HK9/PLL3i7L5zz11FP65z//qZ07d3q7FJ+XlZWlL774Qp9//rm3S/E7Dz30kN5//33t3bvX7Y8o8nc333yzLrroIi1YsMCx7/bbb1dISIjefvttL1bme06ePKnw8HCtWrVKN910k2N/r169NHDgQM2cObNJ62GYzItqa2u1bNkynThxwvHAWZztvvvu00033aT+/ft7uxSft3fvXsXFxSkpKUmZmZkqKirydkk+afXq1frNb36jO++8U23btlXPnj31t7/9zdtl+bzq6mq9/fbbGj16NEHoHNLS0vTvf/9b//3vfyVJ33zzjTZt2qSBAwd6uTLfc/r0adXW1io4OLjO/pCQEG3atKnJ66EP3Qt27dql1NRUnTp1SmFhYXrvvffUtWtXb5flk5YtW6avv/7aa+PI/uS3v/2tFi5cqM6dO+vQoUN6+umn1adPH+3evVvh4eHeLs+n7Nu3T/PmzdPDDz+sxx57TFu3btWDDz6owMBA3Xvvvd4uz2f985//VHl5uUaOHOntUnxSVlaWLBaLUlJSFBAQoNraWj377LPKzMz0dmk+Jzw8XKmpqZoxY4a6dOmiiy66SEuXLtWXX36pjh07Nn1BdjQ5q9Vq37t3r33btm32rKwse1RUlP0///mPt8vyOUVFRfa2bdvav/nmG8e+a665xj5hwgTvFeVHfvrpJ3tERIT9jTfe8HYpPqdFixb21NTUOvseeOABe+/evb1UkX+4/vrr7TfffLO3y/BZS5cutbdr186+dOlS+7fffmtftGiRvU2bNvaFCxd6uzSfVFBQYE9PT7dLsgcEBNivvPJKe2Zmpj0lJaXJa6FnyAsCAwMdybdXr17aunWr/vrXv+r111/3cmW+Zfv27Tp8+LCuuOIKx77a2lp99tlnmjt3rqxWqwICArxYoW9r1aqVLrnkEhUUFHi7FJ8TGxt7Vm9sly5dtGLFCi9V5Pv279+vjz/+WCtXrvR2KT5r8uTJysrK0t133y1Juuyyy7R//35lZ2fT43gOycnJ2rhxo06cOCGLxaLY2FgNHTpUSUlJTV4Lc4Z8gM1mk9Vq9XYZPqdfv37atWuXdu7c6dh+85vfKDMzUzt37iQIncfx48dVWFio2NhYb5fic66++mrl5+fX2fff//5X7du391JFvi8nJ0dt27atM9kVdVVVValZs7o/qwEBAbLZbF6qyD+EhoYqNjZWP/30k9auXavBgwc3eQ30DDWxqVOnauDAgUpISFBlZaWWLFmiTz/9VGvXrvV2aT4nPDxcl156aZ19oaGhioyMPGs/pEmTJumWW25R+/btVVJSounTpysgIEDDhg3zdmk+Z+LEiUpLS9Nzzz2nu+66S1u2bNH8+fM1f/58b5fmk2w2m3JycnTvvfeyXEMDbrnlFj377LNKSEhQt27dtGPHDv35z3/W6NGjvV2aT1q7dq3sdrs6d+6sgoICTZ48WSkpKRo1alTTF9PkA3MGN3r0aHv79u3tgYGB9ujoaHu/fv3s69at83ZZfoM5Q/UbOnSoPTY21h4YGGi/+OKL7UOHDrUXFBR4uyyftWbNGvull15qDwoKsqekpNjnz5/v7ZJ81tq1a+2S7Pn5+d4uxadZLBb7hAkT7AkJCfbg4GB7UlKS/fHHH7dbrVZvl+aT/vGPf9iTkpLsgYGB9piYGPt9991nLy8v90otrDMEAAAMjTlDAADA0AhDAADA0AhDAADA0AhDAADA0AhDAADA0AhDAADA0AhDAADA0AhDAADA0AhDAPA/nnrqKV1++eXeLgNAEyEMAQAAQyMMAQAAQyMMAfBJP/zwg0wm01lb37596z3HYrEoJCREH330UZ397733nsLDw1VVVSVJmjJlii655BK1bNlSSUlJmjZtmmpqauq9bt++ffXQQw/V2TdkyBCNHDnS8dpqtWrSpEm6+OKLFRoaqt/+9rf69NNPG/uxAXhBc28XAADnEh8fr0OHDjlel5aWqn///kpPT6/3nIiICN18881asmSJBg4c6Ni/ePFiDRkyRC1btpQkhYeHa+HChYqLi9OuXbs0btw4hYeH69FHH3W53vvvv1/fffedli1bpri4OL333nsaMGCAdu3apU6dOrl8XQCeRxgC4JMCAgIUExMjSTp16pSGDBmi1NRUPfXUUw2el5mZqeHDh6uqqkotW7aUxWLRBx98oPfee8/R5oknnnD83aFDB02aNEnLli1zOQwVFRUpJydHRUVFiouLkyRNmjRJ//rXv5STk6PnnnvOpesCaBqEIQA+b/To0aqsrNT69evVrFnDo/s33nijWrRoodWrV+vuu+/WihUrFBERof79+zva/OMf/9Arr7yiwsJCHT9+XKdPn1ZERITL9e3atUu1tbW65JJL6uy3Wq2KjIx0+boAmgZhCIBPmzlzptauXastW7YoPDz8vO0DAwN1xx13aMmSJbr77ru1ZMkSDR06VM2b//yfuy+//FKZmZl6+umndcMNN8hsNmvZsmV66aWX6r1ms2bNZLfb6+z75Ryj48ePKyAgQNu3b1dAQECddmFhYY35uAC8gDAEwGetWLFCzzzzjD766CMlJyc7fV5mZqauu+46/ec//9GGDRs0c+ZMx7Hc3Fy1b99ejz/+uGPf/v37G7xedHR0nflLtbW12r17tzIyMiRJPXv2VG1trQ4fPqw+ffo4XScA38DdZAB80u7duzVixAhNmTJF3bp1U2lpqUpLS1VWVnbec9PT0xUTE6PMzEwlJibqt7/9reNYp06dVFRUpGXLlqmwsFCvvPJKnflE53Lttdfqgw8+0AcffKC8vDz98Y9/VHl5ueP4JZdcoszMTI0YMUIrV67U999/ry1btig7O1sffPCBy98BgKZBGALgk7Zt26aqqirNnDlTsbGxju22224777kmk0nDhg3TN998o8zMzDrHBg0apIkTJ+r+++/X5ZdfrtzcXE2bNq3B640ePVr33nuvRowYoWuuuUZJSUmOXqEzcnJyNGLECD3yyCPq3LmzhgwZoq1btyohIaHxHx5AkzLZ/3cgHAAAwEDoGQIAAIZGGALgVwYOHKiwsLBzbqznA8AVDJMB8CvFxcU6efLkOY+1adNGbdq0aeKKAPg7whAAADA0hskAAIChEYYAAIChEYYAAIChEYYAAIChEYYAAIChEYYAAIChEYYAAIChEYYAAICh/X9WNdc7Rk8NvwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.scatterplot(data=df_result,x='z_value',y='total_energy_per_atom')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "127c0ea3-0759-4571-b677-7c2ed0ab9a6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_result.to_csv(save_dir+'/total_energy.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e58d796-3595-4449-847b-afabd9d4d3c4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}