notebooks/collect_vacuum_phase_benchmark26.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Collect vaccum phase DLPNO-CCSD(T)/DFT-level ion-pair energies and compare with ISAPT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pandas import DataFrame as df\n",
    "import scipy as sc\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import r2_score\n",
    "au2kcal = 627.50957099203276"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Globals\n",
    "pattern='FINAL SINGLE POINT ENERGY'\n",
    "position=4\n",
    "\n",
    "mols=[\n",
    "    'C1L1A1',\n",
    "    'C1L1A2',\n",
    "    'C1L1A3',\n",
    "#   'C2L1A1', \n",
    "    'C2L1A2',\n",
    "    'C2L1A3',\n",
    "    'C3L1A1',\n",
    "    'C3L1A2',\n",
    "    'C3L1A3',\n",
    "    'C1L2A1',\n",
    "    'C1L2A2',\n",
    "    'C1L2A3',\n",
    "    'C2L2A1',\n",
    "    'C2L2A2',\n",
    "    'C2L2A3',\n",
    "    'C3L2A1',\n",
    "    'C3L2A2',\n",
    "    'C3L2A3',\n",
    "    'C1L3A1',\n",
    "    'C1L3A2',\n",
    "    'C1L3A3',\n",
    "    'C2L3A1',\n",
    "    'C2L3A2',\n",
    "    'C2L3A3',\n",
    "    'C3L3A1',\n",
    "    'C3L3A2',\n",
    "    'C3L3A3'\n",
    "]\n",
    "\n",
    "Ctypes=np.array([\n",
    "    [1,1],\n",
    "    [1,2],\n",
    "    [1,3],\n",
    "#   [2,1],\n",
    "    [2,2],\n",
    "    [2,3],\n",
    "    [3,1],\n",
    "    [3,2],\n",
    "    [3,3],\n",
    "    [1,1],\n",
    "    [1,2],\n",
    "    [1,3],\n",
    "    [2,1],\n",
    "    [2,2],\n",
    "    [2,3],\n",
    "    [3,1],\n",
    "    [3,2],\n",
    "    [3,3],\n",
    "    [1,1],\n",
    "    [1,2],\n",
    "    [1,3],\n",
    "    [2,1],\n",
    "    [2,2],\n",
    "    [2,3],\n",
    "    [3,1],\n",
    "    [3,2],\n",
    "    [3,3]\n",
    "])\n",
    "\n",
    "DFT=np.array([\n",
    "    'DLPNOCCSDT',\n",
    "    'B3LYPD3BJ_def2SVPD',\n",
    "    'BLYPD3BJ_def2SVPD',\n",
    "    'CAMB3LYP_def2SVPD',\n",
    "    'LC_BLYP_def2SVPD',\n",
    "    'LC_PBE_def2SVPD',\n",
    "    'M062XD3Zero_def2SVPD',\n",
    "    'PBE0D3BJ_def2SVPD',\n",
    "    'PBED3BJ_def2SVPD',\n",
    "    'TPSSD3BJ_def2SVPD',\n",
    "    'wB97XD3_def2SVP',\n",
    "    'wB97XD3_def2TZVP',\n",
    "    'wB97XD3_def2SVPD',\n",
    "    'wB97XD3_def2TZVPD'\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_E(filename):\n",
    "    file = open(filename, \"r\")\n",
    "    iline=0\n",
    "    for line in file:\n",
    "        if re.search(pattern, line):\n",
    "            iline=iline+1\n",
    "            eline=line\n",
    "    if (iline < 1):\n",
    "        print('ERROR: Can\\'t find <',pattern,'> in ',filename)\n",
    "    file.close()\n",
    "    E=float( eline.split()[position] )\n",
    "    return E\n",
    "\n",
    "def calc_E_SB(mol,Ctype):\n",
    "    \n",
    "    dire='../'+directory+'/'+mol+'/'\n",
    "    \n",
    "    file_A_pos=dire+'A_pos.out'\n",
    "    file_A_neg=dire+'A_neg.out'\n",
    "    \n",
    "    E_A_pos=extract_E(file_A_pos)\n",
    "    E_A_neg=extract_E(file_A_neg)\n",
    "    A=E_A_pos+E_A_neg\n",
    "    \n",
    "    file_B=dire+'B.out'\n",
    "\n",
    "    B=extract_E(file_B)\n",
    "\n",
    "    N_D=np.prod(Ctype)\n",
    "    E=np.zeros(N_D)\n",
    "    \n",
    "    i_D=0\n",
    "    \n",
    "    for i_Cm in range(Ctype[1]):\n",
    "        for i_Cp in range(Ctype[0]):\n",
    "        \n",
    "            file_C_pos=dire+'C'+str(i_Cp+1)+'_pos.out'\n",
    "            file_C_neg=dire+'C'+str(i_Cm+1)+'_neg.out'\n",
    "        \n",
    "            E_C_pos=extract_E(file_C_pos)\n",
    "            E_C_neg=extract_E(file_C_neg)        \n",
    "            C=E_C_pos+E_C_neg\n",
    "        \n",
    "            file_D=dire+'D'+str(i_D+1)+'.out'\n",
    "        \n",
    "            D=extract_E(file_D)\n",
    "        \n",
    "            E[i_D]=((B - D) - (A - C))*au2kcal\n",
    "            \n",
    "            i_D=i_D+1\n",
    "        \n",
    "    return E\n",
    "\n",
    "def make_csv(filename,Nrows,Ncols,Header,Values):\n",
    "    \n",
    "    file = open(filename, \"w\")\n",
    "\n",
    "    for icol in range(Ncols):\n",
    "        if icol == Ncols-1:\n",
    "            file.write(Header[icol]+'\\n')\n",
    "        else:\n",
    "            file.write(Header[icol]+',')\n",
    "\n",
    "    for irow in range(Nrows):\n",
    "        for icol in range(Ncols):\n",
    "            if icol == Ncols-1:\n",
    "                file.write(str(Values[icol,irow])+'\\n')\n",
    "            else:\n",
    "                file.write(str(Values[icol,irow])+',')\n",
    "    file.close() \n",
    "    \n",
    "    return 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ion-pair energies estimated with terminating group -CH$_3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "NDFT=len(DFT)\n",
    "Nmol=len(mols)\n",
    "\n",
    "Eip=np.zeros([NDFT,Nmol])\n",
    "Eip_fluct=np.zeros([NDFT,Nmol])\n",
    "\n",
    "for iDFT in range(NDFT):\n",
    "\n",
    "    directory='IO_files_benchmark26/'+DFT[iDFT]+'/NN1'\n",
    "    \n",
    "    for imol in range(Nmol):\n",
    "        mol=mols[imol]\n",
    "        Ctype=Ctypes[imol]\n",
    "        E=calc_E_SB(mol,Ctype)\n",
    "        Eip[iDFT,imol]=np.mean(E)\n",
    "        Eip_fluct[iDFT,imol]=np.std(E)\n",
    "\n",
    "tmp=make_csv('../csv_files/benchmark26_term_group_1.csv',Nmol,NDFT,DFT,Eip)\n",
    "tmp=make_csv('../csv_files/benchmark26_term_group_1_fluct.csv',Nmol,NDFT,DFT,Eip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ion-pair energies estimated with terminating group -CH$_2$CH$_3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eip=np.zeros([NDFT,Nmol])\n",
    "Eip_fluct=np.zeros([NDFT,Nmol])\n",
    "\n",
    "for iDFT in range(NDFT):\n",
    "\n",
    "    directory='IO_files_benchmark26/'+DFT[iDFT]+'/NN2'\n",
    "    \n",
    "    for imol in range(Nmol):\n",
    "        mol=mols[imol]\n",
    "        Ctype=Ctypes[imol]\n",
    "        E=calc_E_SB(mol,Ctype)\n",
    "        Eip[iDFT,imol]=np.mean(E)\n",
    "        Eip_fluct[iDFT,imol]=np.std(E)\n",
    "\n",
    "tmp=make_csv('../csv_files/benchmark26_term_group_2.csv',Nmol,NDFT,DFT,Eip)\n",
    "tmp=make_csv('../csv_files/benchmark26_term_group_2_fluct.csv',Nmol,NDFT,DFT,Eip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ion-pair energies estimated with terminating group -CH$_2$CH$_2$CH$_3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Eip=np.zeros([NDFT,Nmol])\n",
    "Eip_fluct=np.zeros([NDFT,Nmol])\n",
    "\n",
    "for iDFT in range(NDFT):\n",
    "\n",
    "    directory='IO_files_benchmark26/'+DFT[iDFT]+'/NN3'\n",
    "    \n",
    "    for imol in range(Nmol):\n",
    "        mol=mols[imol]\n",
    "        Ctype=Ctypes[imol]\n",
    "        E=calc_E_SB(mol,Ctype)\n",
    "        Eip[iDFT,imol]=np.mean(E)\n",
    "        Eip_fluct[iDFT,imol]=np.std(E)\n",
    "\n",
    "tmp=make_csv('../csv_files/benchmark26_term_group_3.csv',Nmol,NDFT,DFT,Eip)\n",
    "tmp=make_csv('../csv_files/benchmark26_term_group_3_fluct.csv',Nmol,NDFT,DFT,Eip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error metrics w.r.t. ISAPT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#benchmark26_ISAPT=pd.read_csv('../csv_files/benchmark26_ISAPT_Psi4_1.4.csv',header=0)\n",
    "benchmark26_ISAPT=pd.read_csv('../csv_files/benchmark26_ISAPT_junccpvdz.csv',header=0)\n",
    "benchmark26_NN1=pd.read_csv('../csv_files/benchmark26_term_group_1.csv',header=0)\n",
    "benchmark26_NN2=pd.read_csv('../csv_files/benchmark26_term_group_2.csv',header=0)\n",
    "benchmark26_NN3=pd.read_csv('../csv_files/benchmark26_term_group_3.csv',header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                        MAD-CH3    STD-CH3   MPAD-CH3  MAD-CH2CH3  STD-CH2CH3  \\\n",
      "DLPNOCCSDT            28.528980  31.141644  24.603804   23.791993   29.376541   \n",
      "B3LYPD3BJ_def2SVPD    28.648859  31.670699  24.726923   24.035482   29.841745   \n",
      "BLYPD3BJ_def2SVPD     28.381362  31.550243  24.519415   23.581342   29.693238   \n",
      "CAMB3LYP_def2SVPD     28.459683  31.611785  24.559248   24.038311   29.841327   \n",
      "LC_BLYP_def2SVPD      28.273705  31.483325  24.396842   23.878983   29.726607   \n",
      "LC_PBE_def2SVPD       29.305121  31.791713  25.261871   25.055548   30.017057   \n",
      "M062XD3Zero_def2SVPD  29.049511  31.630669  24.980936   24.473769   29.746132   \n",
      "PBE0D3BJ_def2SVPD     29.133940  31.739473  25.123956   24.534155   29.882380   \n",
      "PBED3BJ_def2SVPD      28.859709  31.577874  24.914151   23.970250   29.676830   \n",
      "TPSSD3BJ_def2SVPD     28.887597  31.577270  24.932768   24.108214   29.720558   \n",
      "wB97XD3_def2SVP       37.195439  33.366940  31.998489   30.821205   30.832675   \n",
      "wB97XD3_def2TZVP      31.150965  31.963580  26.828929   25.826068   29.922650   \n",
      "wB97XD3_def2SVPD      28.873536  31.672669  24.913666   24.322714   29.911260   \n",
      "wB97XD3_def2TZVPD     28.725031  31.406411  24.848099   24.186586   29.678818   \n",
      "\n",
      "                      MPAD-CH2CH3  MAD-CH2CH2CH3  STD-CH2CH2CH3  \\\n",
      "DLPNOCCSDT              20.738663      22.272933      29.169834   \n",
      "B3LYPD3BJ_def2SVPD      20.963397      22.567824      29.578933   \n",
      "BLYPD3BJ_def2SVPD       20.605754      22.415050      29.422307   \n",
      "CAMB3LYP_def2SVPD       20.940956      22.569052      29.587722   \n",
      "LC_BLYP_def2SVPD        20.787332      22.632296      29.471135   \n",
      "LC_PBE_def2SVPD         21.763707      23.506255      29.726066   \n",
      "M062XD3Zero_def2SVPD    21.247862      22.610905      29.487768   \n",
      "PBE0D3BJ_def2SVPD       21.365097      22.872665      29.594495   \n",
      "PBED3BJ_def2SVPD        20.923140      22.699076      29.394555   \n",
      "TPSSD3BJ_def2SVPD       21.036067      22.546568      29.438226   \n",
      "wB97XD3_def2SVP         26.835606      28.014552      30.494182   \n",
      "wB97XD3_def2TZVP        22.485810      23.742324      29.661911   \n",
      "wB97XD3_def2SVPD        21.193019      22.737548      29.634648   \n",
      "wB97XD3_def2TZVPD       21.108437      22.573988      29.455429   \n",
      "\n",
      "                      MPAD-CH2CH2CH3  \n",
      "DLPNOCCSDT                 19.374545  \n",
      "B3LYPD3BJ_def2SVPD         19.637786  \n",
      "BLYPD3BJ_def2SVPD          19.501388  \n",
      "CAMB3LYP_def2SVPD          19.621520  \n",
      "LC_BLYP_def2SVPD           19.640757  \n",
      "LC_PBE_def2SVPD            20.391406  \n",
      "M062XD3Zero_def2SVPD       19.635078  \n",
      "PBE0D3BJ_def2SVPD          19.896641  \n",
      "PBED3BJ_def2SVPD           19.740614  \n",
      "TPSSD3BJ_def2SVPD          19.630984  \n",
      "wB97XD3_def2SVP            24.427282  \n",
      "wB97XD3_def2TZVP           20.667868  \n",
      "wB97XD3_def2SVPD           19.777256  \n",
      "wB97XD3_def2TZVPD          19.657414  \n"
     ]
    }
   ],
   "source": [
    "methods=benchmark26_NN2.columns[:]\n",
    "\n",
    "diff_NN1=benchmark26_NN1.copy()\n",
    "diff_NN2=benchmark26_NN2.copy()\n",
    "diff_NN3=benchmark26_NN3.copy()\n",
    "\n",
    "for imethod in methods:\n",
    "    ref=benchmark26_ISAPT.ISAPT.values[:].copy()\n",
    "    diff_NN1[imethod]=df(benchmark26_ISAPT.ISAPT.values[:]-diff_NN1[imethod].values,columns=[[imethod]])\n",
    "    diff_NN2[imethod]=df(benchmark26_ISAPT.ISAPT.values[:]-diff_NN2[imethod].values,columns=[[imethod]])\n",
    "    diff_NN3[imethod]=df(benchmark26_ISAPT.ISAPT.values[:]-diff_NN3[imethod].values,columns=[[imethod]])\n",
    "    \n",
    "rdiff_NN1=benchmark26_NN1.copy()\n",
    "rdiff_NN2=benchmark26_NN2.copy()\n",
    "rdiff_NN3=benchmark26_NN3.copy()\n",
    "\n",
    "for imethod in methods:\n",
    "    ref=benchmark26_ISAPT.ISAPT.values[:].copy()\n",
    "    rdiff_NN1[imethod]=df(diff_NN1[imethod].values/ref,columns=[[imethod]])\n",
    "    rdiff_NN2[imethod]=df(diff_NN2[imethod].values/ref,columns=[[imethod]])\n",
    "    rdiff_NN3[imethod]=df(diff_NN3[imethod].values/ref,columns=[[imethod]])\n",
    "\n",
    "metrics=df()\n",
    "metrics['MAD-CH3']=np.mean(diff_NN1.abs())\n",
    "metrics['STD-CH3']=diff_NN1.std(ddof=0)\n",
    "metrics['MPAD-CH3']=np.mean(rdiff_NN1.abs()*100)\n",
    "\n",
    "metrics['MAD-CH2CH3']=np.mean(diff_NN2.abs())\n",
    "metrics['STD-CH2CH3']=diff_NN2.std(ddof=0)\n",
    "metrics['MPAD-CH2CH3']=np.mean(rdiff_NN2.abs()*100)\n",
    "\n",
    "metrics['MAD-CH2CH2CH3']=np.mean(diff_NN3.abs())\n",
    "metrics['STD-CH2CH2CH3']=diff_NN3.std(ddof=0)\n",
    "metrics['MPAD-CH2CH2CH3']=np.mean(rdiff_NN3.abs()*100)\n",
    "\n",
    "print(metrics)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scatterplot for ISAPT vs $\\omega$B97XD3/def2SVPD with terminal group -CH$_2$CH$_2$CH$_3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGBCAYAAAAgxmXKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABTeElEQVR4nO3deXxU1f3/8dfHsG9CDA27QOICCiJEDCIUaREQRW2rVUCktaJfWW2t+m1ttf4A26rIphaUFmW1tF+tCyDSuqBCS6DgAipEVgUEkoCsIvn8/rh3wmQyy00ya/J5Ph73QebeM3fO3JB5z7n33HNEVTHGGGPi7YxEV8AYY0z1ZAFkjDEmISyAjDHGJIQFkDHGmISwADLGGJMQFkDGGGMSwgLIGGNMQqR0AIlIMxGZKyJ7ROSoiGwQkaEBZZq4ZQ66y1wRaZygKhtjjHGldAABzwMdgGuBC93Hc0Wkt1+ZBUBXYIC7dAXmxrmexhhjAkgqj4QgIoeBMar6F79124HpqvqYiHQANgKXq+p77vbLgZXA+ar6aah9Z2RkaNu2bWNaf2OMSRV79+6lqKiIw4cP71fVptHYZ41o7CSB3gVuFJGXgULgGqApsMLd3gM4DLzv95z3gCPAZUDIAGrbti15eXmxqLMxxqQMVWXKlCm8/PLL3Hzzzdxxxx3bo7XvVD8FdyOgwH7gBDAfuFlV17vbmwH71K+Z5/78lbutFBEZKSJ5IpK3b9++WNfdGGOSWmD43H777VHdf9IFkIhMEBGNsPRxi08AMoDvAznAo8DzInJRRV5bVWepao6q5jRtGpUWpjHGpKRg4SMiUX2NZDwFNwWYF6HMDhHJAsYAXVR1g7t+g4j0ctf/DNgDNBUR8bWCxDmC33G3GWOMCRCP8IEkDCBV3Y9zSi0sEann/ngqYNMpTrfsVgENcK4F+a4D9QDqU/q6kDHGGOIXPpCEAVQOnwBbgKdE5B7gAHAd0A+nWzaquklElgEzRWSk+7yZwKvhesAZY0x1FM/wgSS8BuSVqp4ErgL2Aa8AHwDDgZ+o6it+RYcAG4DX3WUDcEt8a2uMMckt3uEDqd0CQlU3Az+MUKYQGBafGhljTOpJRPhACreAjPHJz89n9NjxpGdkkpaWRnpGJqPHjic/Pz/RVTMm6SUqfMACyKS4pUuX0jUnl1WbhX5jljBiym76jVnCqs3C+R07ISIRl8bpjRL9NoxJiESGD6T4UDyxlJOTozYSQnLLz8+na04u371tHpntLymzfe/na1gx+3p+8u9jpGeF3s+EGs4fojHVSUXDR0TWqmpONOpgLSCTsp6YOp3s3GFBwwcgs/0lnJv7E/Jm1IxzzYxJbolu+fhYAJmUtWDBQrJzw/cvOS/3p3y0IC1ONTIm+SVL+IAFkElhBwv30yC9ddgyDdJbcazoRJxqZExyS6bwAQsgk8LObJLB4YKdYcscLthF3ca141QjY5JXsoUPWACZFDZkyM1sWR1+2MBPV/+ZC4cEjtZkTPWSjOEDFkAmhd09bgxbVs9j7+drgm7f+/kaPlv9F3JGn4xzzYxJHskaPpDiIyGY6i0rK4tFC57npiHDyM4dSnbuLTRIb8Xhgl1sWT2XDf+cxQ3/F74LtjFVWTKHD1gLyKS4gQMHsi5vNT3PFVZMH8Scu1uyYvogep4rnDp5jOyBia6hMYmR7OEDdiNqSHYjauprnN6Ig4VfRyx3ZpOGFBUcikONjImPWIZPNG9EtVNwpsqyUDHVUSq0fHzsFJwxxlQRqRQ+YAFkjDFVQqqFD1gAGWNMykvF8AELIGOMSWmpGj5gAWSMMSkrlcMHLICMMSYlpXr4gAWQMcaknKoQPmABZIwxKaWqhA9YABljTMqoSuEDFkDGGJMSqlr4gAWQMcYkvaoYPmABZIwxSa2qhg9YABljTNKqyuEDKR5AIpIlIi+KyD4ROSQifxWRzIAyTURkrogcdJe5ItI4QVU2xhhPqnr4QAoHkIjUB5YDAvQFegK1gFdExP99LQC6AgPcpSswN761NcYY76pD+EBqzwfUE2gH5KhqIYCI3AoU4gTSChHpgBM6l6vqKrfMHcBKETlPVT9NTNWNMSa46hI+kMItIKA2oMBxv3XHgWLgcvdxD+Aw8L5fmfeAI8BlcaijMcZ4Vp3CB1I7gFbjhMujIlLfPSX3GJAGNHfLNAP2qd+84+7PX7nbShGRkSKSJyJ5+/bti/kbMMYYn+oWPpCEASQiE0REIyx9VHUfcAMwEPgaOAg0BtbhtILKTVVnqWqOquY0bdo0Wm/JGGPCqo7hA8l5DWgKMC9CmR0AqrocyBKRDOBbVS0SkT3A5265PUBTERFfK0ic3+p33G3GGJNQ1TV8IAkDSFX3A/sr8BxEpC9OuLzsbloFNMC5FuS7DtQDqE/p60LGGBN31Tl8IAkDqDxE5CfAJzjXdHoAU4EnfL3bVHWTiCwDZorISPdpM4FXrQecMSaRqnv4QIoHEHAe8AiQDmwDJgJPBJQZAkwHXncfvwyMjlP9jDGmDAsfR0oHkKreD9wfoUwhMCw+NTLGmPAsfE5Lul5wxpjo2r17N+efk82ePdbvJtEsfEqzADKmivvjxIl8tWsnf5w4IdFVqdYsfMqyADKmCtu9ezfPPTeHf/bvxnPPzbFWUIJY+ARnAWRMFfbHiRO5Nas5F2c0Ynj75tYKSgALn9AsgIyponytn3s7tgLg3o6trBUUZxY+4VkAGVNF+Vo/zevVBqB5vdrWCoojC5/ILICMqYICWz8+1gqKDwsfbyyAjKmCAls/PtYKij0LH+8sgIypYkK1fnysFRQ7Fj7lYwFkTBUTqvXjY62g2LDwKT/xm6vN+MnJydG8vLxEV8OYcmvfpg1bd+6MWK5d69Z8vmNHHGpU9VWn8BGRtaqaE419pfRYcMaYsixU4qs6hU+02Sk4Y4ypIAufyrEAMsaYCrDwqTwLIGOMKScLn+iwADLGmHKw8ImekJ0QROTlCuzvdlXdW4n6GGNM0rLwia5wveCuBv4KHPO4ryFA/UrXyBhjkpCFT/RF6oY9VlW/8rIjEflRFOpjjDFJx8InNsJdA7oCKCjHvgYCX1SuOsaYZJCfn8/oseNJz8gkLS2N9IxMRo8dT35+fqKrFncWPrETMoBU9W1V/dbrjlT1XVU9EZ1qmcponN4IEYm4NE5vlOiqmiS0dOlSuubksmqz0G/MEkZM2U2/MUtYtVnompPL0qVLE13FuLHwia2QQ/GISLrXnahqeVpKKSGVh+IRER7w8NVhQg3nD8wYn/z8fLrm5PLd2+aR2f6SMtv3fr6Gt2cPY13earKyshJQw/ix8AkumkPxhDsFtx/YF2HxlTHGVAFPTJ1Odu6woOEDkNn+ErJzhzJ12ow41yy+LHziI1wnhCviVgtjTFJYsGAh/cYsCVsmO/cW5k8fxLSpT8SpVvFl4RM/IQNIVd+OZ0WMMYl3sHA/DdJbhy3TIL0VRUX741Sj+LLwiS/PIyGISG0R+amIPCYij4rICBEJPuFIFIjISBF5U0SKRERFpG2QMk1EZK6IHHSXuSLSOKBMJxF5W0SOicgXIvJbsf9RxgR1ZpMMDheEn8rhcMEuGjfOiFON4sfCJ/48BZCIdAQ2A5OBS4FcYArwmYh0iFHd6gHLgYfClFkAdAUGuEtXYK5vo4g0At4A9gKXAOOAXwI/j0mNjUlxQ4bczJbV88KW2bJ6LkOH3hynGsWHhU9ieJ0PaCrwX+AWVT0EJR/u83CCqH+0K6aqU9zXCdrbwg2+AcDlqrrKXXcHsFJEzlPVT4GhOEF2q6oeAz4SkfOBn4vIZLUuYMaUcve4MczNyaXlBQNC9oLbsno+f52xOgG1iw0Ln8TxegquJ/ArX/gAuD//Grg8FhXzoAdwGHjfb917wBHgMr8yK93w8XkdaAG0jUMdjUkpWVlZLFrwvNPV+pWHObRvK8WnTnJo31bWvfIwb88exqIFz1eZLtgWPonlNYCOA42DrD/T3ZYIzYB9/q0Y9+ev3G2+MoGDo+7121aKe90pT0Ty9u2z3uWmeho4cCDr8lbT81xhxfRBzLm7JSumD6LnucK6vNUMHDgw0VWMCgufxPN6Cu4V4BkRuR3wtb17ADMBz6Nmi8gEnFZTOFeo6lte9xlNqjoLmAXOjaiJqEM0nNmkIRNqfO2pnDHBZGVlMW3qE9bV2sSU1wAaBzwHrAROuevOwAmf8eV4vSk4143C8Tqh/R6gqYiIrxXk9m77jrvNVyYz4HmZftuqpKKCQ5ELGZMAjdMbcbDQ25ejWP0/tvBJHp4CSFWLgGtF5BzgfHf1JlXdUp4XU9X9OKMnRMMqoAFOS8x3HagHzpQQ7/uV+YOI1FFV36nCfsCXwLYo1cMY49HBwq89DhMVOaQqwsInuXhtAQGgqptxumPHnIg0w7lOc667qqN7j88OVS1Q1U0isgyYKSIj3TIzgVfdHnDgdNN+EJjjnv47F7gf+J31gDOmerHwST6eA0hErscZnuc7BHReUNUbo1wvgDtxwsPnNfffnwBz3J+HANNxeraBc0pwtF+9DopIP+BJIA8oBB7HuZ/JGFNNWPgkJ08BJCKPA2Nwujnv5fR1oJhR1YcIfxMqqloIDItQ5kOgd9QqZoxJKRY+yctrC+hW4AZV/UcsK2OMMdFk4ZPcvN4HdBT4JJYVMcaYaLLwSX5eA+j3wL0iUq5OC8YYkwgWPqnBa6A8A1wDfCEinwEn/Teqat9oV8wYYyrCwid1eA2gP+GM+bYMpxOCdWE2xiQdC5/U4jWAfgxcr6pvxLIyxpiqLZbDRFn4pB6vAbQP+CKWFTHGVH02vI7x57UTwoPAwyLSIJaVMcaY8rLwSV1eW0C/xJk/Z6+I7KBsJ4TOUa6XMcZEZOGT2rwG0N9iWgtjjAngdeTsmrXTmDHtKQufFCThxuQUkSbucDfVTk5Ojubl5SW6GsZUWyLiceRsKC4utvCJExFZq6o50dhXpGtAe0TkXyIyTkTOjsYLGmNMtFn4pKZIAXQ2sAi4EvhERNaLyO9EpGvsq2aMMaYqCxtAqrpHVWep6iCgKTABaA+sEJHtIjJdRL4vImnxqKwxxpiqw2s3bFT1sKr+TVVvwQmjnwLFwLPAPhEZGqM6GmOMqYIqNLioqp4C/uku40SkC1AzivUyxhhTxYUMoHJc51FV/W+U6mOMMaaaCNcCysMZdDRS9xIF7BqQMcaYcgkXQO3iVgtjjPET7v5EU3WEDCBV3R7PihhjDJweXqdm7TQm1DgVsXxFRs42ycFzLzgAEWkhIrki0tt/iVXljDHVi//YbjOmPUVxcTGqGnaJ1QjbJvY89YITkRbAAqA3p68L+beR7RqQMaZSbGDR6sdrC2gKcAroCBwFegE3AJuAATGpmTGm2rDwqZ683gf0XWCQqn4iIgrsU9X3ROQE8P8AmynVGFMhFj7Vl9cWUF1gv/tzAfAd9+eNgM0FZIypEAuf6s1rAH0CnO/+vB640x0dexQ2VbcxpgIsfIzXAJoKNHN/fhhndOzPgbuAX8WgXojISBF5U0SKRERFpG2QMr8WkfdE5Ih7ajDYftqIyCtumf0iMk1EasWizsYYbyx8DHi8BqSq8/1+XueGwfnADlXdH/KJlVMPWA78A3giRJnawP8BbxEkCN1Rul8DDuB0nDgLeA6nF9+YqNfYGBORhY/x8doNuxZwhqoeB1DVo8A6EakjIrVU9ZtoV0xVp7ivHXLmPVX9rVvmRyGKXAlcAJytqjvdsvcCz4rIr1XVbiAwJo4sfIw/r6fgFuOcbgt0J/DX6FUn6noAm3zh43odp+XULbCwe9ovT0Ty9u3bF686GlMtWPiYQF4DqCfO6bBAbwCXRa86UdcM2Buwbj/OPU3NAgu7k+/lqGpO06ZN41E/Y6oFCx8TjNcAqgd8G2R9MeB5ICYRmeB2KAi39PG6P2NM8rPwMaF4DaAPgJuDrB8CfFSO15sCdIiw/Kcc+4tkD5AZsC4DZ+igPVF8HWOqlPz8fEaPHU96RiZpaWmkZ2Qyeux48vPzy7UfCx8TjteREB4G/iEi2cC/3HXfwxmO53qvL+b2mItVr7lgVgEPiEgrVd3lrusHnADWxrEexqSMpUuXctOQ4WTnDqPfmCU0SG/N4YKdrFo9j7k5uSxa8DwDBw6MuB8LHxOJ127YS0TkGuABYJq7+r/AYFVdGouKiUgznOs057qrOopIY5yu3wVumTZAOtDWfdzFLbtFVQ/jXLf6GHheRH6B0w37UeAZ6wFnTFn5+fncNGQ4371tHpntLylZ36hpO7pe8xtaXjCAm4YMY13earKyskLux8LHeOF5OgZVXaaql6tqfXe5XFWXikjNGNXtTpyQ892D9Jr7eLBfmYfddY+6j//rLjlunU8Bg3AGUH0PeAH4O3BPjOpsEmz37t2cf042e/bYGdaKeGLqdLJzh5UKH3+Z7S8hO3coU6fNCLkPCx/jlacAEpH/F2J9LZwP9KhT1YdUVYIsc/zKjAhR5i2/MjtU9WpVraeqZ6nqWFU9EYs6m8T748SJfLVrJ3+cOCHRVUlJCxYsJDt3WNgy2bm3MH/+wqDbLHxMeXhtAd0mImP9V7jh839Am6jXypgK2L17N889N4d/9u/Gc8/NsVZQBRws3E+D9NZhyzRIb0VRUdlLuRY+pry8BtBA4CERGQol4fMi0BqnM4IxCffHiRO5Nas5F2c0Ynj75tYKqoAzm2RwuGBn2DKHC3bRuHFGqXUWPqYiPAWQqm4ArgWeEpEf4rR8WgJ9VfVADOtnjCe+1s+9HVsBcG/HVtYKqoAhQ25my+p5YctsWT2XoUNP35Vh4WMqqjydEFbi3PezEAsfk2R8rZ/m9WoD0LxebWsFVcDd48awZfU89n6+Juj2vZ+vYcvq+YwbOxqw8DGVI6pBZzFARF4O8ZwcnKkYCnwrVHVwiLIpKycnR/Py8hJdDePB7t27ueDcc/j4mktKAghg99ETXPjqGj7+bAvNmpUZecmEcPo+oKFk595Cg/RWHC7YxZbVc9myen7JfUAWPtWTiKxV1ZCDRJdHuBbQgRDL68DmgHXGJExg68fHWkEVM3DgQNblrabnucKK6YOYc3dLVkwfRM9zhXV5qy18TNSEbAFVd9YCSg2hWj8l260VFHWVCZ9Dhw7x1VdfcfLkyRjX0lREzZo1+c53vkOjRo1ClolmC8jrUDzGJKVQrR8f/1bQ5Omhb5403lQ2fPbu3UvLli2pW7eutZiSjKpy7NgxvvjiC4CwIRQt4a4BLQFuVtWDnnYk8nfgLlUNnP4gJVkLKDW0b9OGrTvDdxsGaNe6NZ/v2BGHGlVdlT3ttmXLFlq0aEG9evViWEtTWUePHuXLL78kOzs76PZ4tYD6A83caa296AfUr3yVjPHOQiU+onHN5+TJk9StWzdGNTTRUrdu3bidIg0XQAJsjEstjDFJK5odDuy0W/KL5+8oXABdUYH9fVHRihhjko/1djOxFDKAVPXteFbEGJNcLHxMrHkeCcEYU31Y+HhTWFhIZmZmuWeKjYYbbriBxx9/PO6vG00WQCkmWlMlGxOKhY93kyZN4qqrrio1Od/evXsZN24cWVlZ1K5dm5YtWzJw4ECWLFlSUmbEiBFcffXVZfaXl5eHiLBt27aIr/3b3/6WiRMncvCgp47KSckCKIUsXbqUrjm5rNos9BuzhBFTdtNvzBJWbRa65uSydGlMJqc11YiFj3dHjx7l2Wef5bbbbitZt23bNrp27crrr7/OI488wgcffMCKFSsYNGgQd955Z1Rfv1OnTrRv355588IPHpvM7EbUFBGtqZJN1dc4vREHC7+OWO7MJg0pKjg9M30qhc/u3bu5oncv3lr5bsJGuFiyZAkiQs+ePUvW3XXXXYDTkmnQoEHJ+g4dOjBsWPiJ/ipi8ODBLFy4kFGjRkV93/EQtgUkImFnphKROiIyPLpVMsFEY6pkUz0cLPyaB74l4uIfUqkUPpAcM9+uXLmSbt26lRyngoICli1bxqhRo0qFj0/jxo2jXofu3bvzn//8h2PHjkV93/EQ6RTcNhF5UURC3WB6JvCXKNfJBFHZqZKNCSXVwidZZr7dvn07LVq0KHm8ZcsWVJUOHTp4ev6yZcto0KBBqaV3794l23fu3EmfPn3o2LEjnTt3ZvHixWX20aJFC06ePMmXX34JQFFRETk5OXTp0oULL7yQZ555ppLvMrYiBZAAlwD/FpH2caiPCaEyUyUbE0qqhQ/Ef+bb+fPnlwqJlStXAnDs2DHq1KlTUq68Azv37t2b9evXl1oWLFhQsr1GjRpMmTKFjRs3snz5csaPH8+RI0dK7cM3soSvBdSwYUPeeecd1q9fz7///W8mTZrEgQMHQr6HRIt0DUhxptyeCqwRkR+r6orYV8sE8k2V3Khpu5Blgk2VbEw4qRY+vtbPx9c4p6Lv7diKC5+bw72/fiBm14IGDx7MpZdeWvK4ZcuWAGRkZFBYWFiy/pxzzkFE2LRpE9dff33E/darV6/MeGtFRUUlPzdv3pzmzZsD0KxZMzIyMigoKKB+/dMnpAoKnGnZmjZtCkBaWlrJWHsnTpxAVVHVkO8h0by0gAqBgcBsYImI3B3zWpkyKjJVsjGRpFL4QGJmvm3YsCHZ2dkli6/VcfHFF7Nx4+nRytLT0+nfvz8zZszg8OHDZfbjHy7ltXbtWk6dOkXr1qXPgnz00Ue0bNmSzMzMUq9z0UUX0apVK375y1+SkZER8j0kmqdu2Oq4F7gVmCAic0SkVmyrZvyVd6pkY7xIpfDxtX7u7diq1Pp7O7ZKyLWg/v37s2nTJg4cOD0n55NPPomqkpOTw+LFi/n000/55JNPePrpp+ncuXOFXqegoIDhw4cza9asMttWrlxJ//79S61r3LgxGzZsYOvWrSxYsIC9e5N3goJy3QekqguBXkAf4B0gOdpx1UBWVhaLFjzP27OHse6Vhzm0byvFp05yaN9W1r3yMG/PHsaiBc9bF2xTLqkSPpB8M9926tSJ7t27s2jRopJ17du3Z926dfTr14/77ruPzp0707dvX15++eWgARLJiRMnuO6667j//vu57LLLSm07fvw4L774IrfffnvQ52ZmZnLRRRclzfWeYMLOiCoixUAzVf0qYH1T4G9AJ+BMVfU6ZUPKSNb5gPLz85k6bQbz5y+kqGg/jRtnMHTozYwbO9rCxwDOaMYPfBu53IQa5b9wXhmbNm3y3EMsULLOfLts2TLGjRvHxo0bSUuL7segqjJkyBDOO+88HnrooTLbn3zySf7xj3+wfPnyknV79+6lXr16NGzYkIMHD9KzZ08WLlxIp06dyvXa4X5X0ZwPKFILaDtwKnClqu7D6ZywEOc6UdSJyEgReVNEikRERaRtwPa2IjJbRD4XkWPuv4+ISN2Acm1E5BUROSIi+0VkWiqfPszKymLa1Cc4sH8Pp779lgP79zBt6hMWPqZKK8/Mt/E0YMAARo0axa5du6K+7/fee48XXniBl156iS5dutClSxc+/PDDku01a9Zk+vTppZ6zfft2evXqxUUXXUSvXr0YM2ZMucMnniK1gIYDL6jqifhVqeS1xwN1gWPAE0A7Vd3mt30AcBNOCG4GOgCzgNdUdaRbJg1YDxwAfg6cBTwH/F1Vx4R7/WRtARkTSUVHQoi1yrSAbObb+IpXCyhSN+y/AMuAryKUizpVnQIgIkHfqKouw6mbz+ciMhH4f8BId92VwAXA2aq6093fvcCzIvJrVY3fX58xcZKqw+uEY6FSNXnphp1KGuF0G/fpAWzyhY/rdaA20C2eFTMm3qpK+Jiqy0svuPhdpawEETkbuAd4ym91MyCwD+J+nOtaZa5Uuted8kQkb9++fTGrqzH+du/ezfnnZEe1G7GFj0kFXgLoRRH5V7jF64uJyAS3Q0G4pU9534SIZOKcjnsD53pRhajqLFXNUdUc353FxsRatAfWtPAxqcLLdAyfAkej9HpTgEiTV5TrZK+INAP+BXwE3KKle1XsAXoGPCUDSHO3GZNQ/gNrfj8KQ8pY+JhU4iWA/jfwPqCKUtX9OKfAokJEmgNvAh8DN6tq4N0Pq4AHRKSVqvr6SfYDTgBro1UPYyoq2MCak6dXbEoNCx+TaiKdgkvY9R8RaSYiXYBz3VUdRaSLiKS721sAb+O0ZMYDGe5zmrndrwGW44TT8yJysYh8H3gUeMZ6wJlECxxapjJDylj4mFSUzL3g7gT+C8x3H7/mPh7sPr4SOAf4Ls5pu91+S2sAVT0FDMI5hfge8ALwd5zOCsYkVLQG1rTwMakqUgCNxLkRNO5U9SFVlSDLHHf7nBDbxf+GVVXdoapXq2o9VT1LVccm4sZaY/xFa2BNCx+TyiIF0DTgCxF5RkRy41EhY6qDaAysaeFjUl2kAGoB/Arnps33ReRjEblbRM6KfdWMqZpCtX58vLSCLHxMVRA2gFS1SFVnqGpXIAd4C/gNTqvoryLSP9zzjTFlVXZgTQuf5FFYWEhmZib5+flxf+0bbriBxx9/PO6vG1W+KVu9LjjD2AwBVgDfAlvLu49UWLp166bGxEK71q0Vp4dp2KVd69ZlnltcXKyTJ0/WPn366MyZM7W4uDgB76BiNm7cmOgqRN0999yjI0aMKLVuz549OnbsWG3fvr3WqlVLW7RooQMGDNDXXnutpMytt96qgwYNKrO/NWvWKKBbt26N+NoffPCBNmnSRIuKiir9PgKF+10BeRqlz1kv9wEFBtYJEVkHdHWX5pUNQWOqk4oOrKnVuOWTjCN8Hz16lGeffZZXXnmlZN22bdvo2bMnDRs25JFHHuGiiy6iuLiYf/7zn9x5553siOKgqp06daJ9+/bMmzePUaNGRW2/8eQ5gESkAc70Bz8FLsUZIeERnOkNjDExVJ3DB+Bg4dceJ9mLHFLRsmTJEkSEnj1PD7Zy1113AZCXl0eDBg1K1nfo0IFhw4ZFvQ6DBw9m4cKFKRtAEceCE5FeIvIXnPtrpuAEz3dVtaOqPqpRGiXBGBNcdQ+fZLVy5Uq6detW8rsoKChg2bJljBo1qlT4+DRu3DjqdejevTv/+c9/OHYsIXfLVFrYABKRz3A6HlyAc/Nmc1X9iaq+G4e6GVPtWfgkr+3bt9OiRYuSx1u2bEFVPU+6t2zZMho0aFBq6d27d8n2nTt30qdPHzp27Ejnzp1ZvHhxmX20aNGCkydP8uWXX5as27p1K1dccQUdO3akU6dOHDlypBLvMrYitYCWAherandVnamq8WvfGlPNWfgkh/nz55cKiZUrVwJw7Ngx6tSpU1JOw8wuHUzv3r1Zv359qWXBggUl22vUqMGUKVPYuHEjy5cvZ/z48WXCpG7duiV18RkxYgQPP/wwGzdu5O2336Z27doh30Oihb0GpKrj/B+LSBbO1NfgTPQW/76HxniUn5/PE1Ons2DBQg4W7ufMJhkMGXIzd48bQ1ZWVqKrF5aFT/IYPHgwl156acnjli1bApCRkUFh4en5L8855xxEhE2bNnH99ddH3G+9evXIzs4uta6oqKjk5+bNm9O8udPHq1mzZmRkZFBQUED9+vVLyhQUFADgmz7m448/pmbNmvTq1QuA9PT0sO8h0bzMB4SInCUiLwGbgZfc5TMR+YfdlGqS0dKlS+mak8uqzUK/MUsYMWU3/cYsYdVmoWtOLkuXLk10FUOy8EkuDRs2JDs7u2TxtTouvvhiNm7cWFIuPT2d/v37M2PGDA4fPlxmP/7hUl5r167l1KlTtG7dutT6jz76iJYtW5KZmQnA5s2badCgAddccw1du3Zl0qRJYd9DonkKIOBZIBvoBdRxl95AO+CZ2FTNmIrJz8/npiHD+e5t8+h6zW9o1LQdZ6TVoFHTdnS95jd897Z53DRkeEJuHozEwid19O/fn02bNnHgwIGSdU8++SSqSk5ODosXL+bTTz/lk08+4emnn6Zz584Vep2CggKGDx/OrFmzymxbuXIl/fufHg/g22+/ZeXKlTz11FOsWrWKN954gzfeeKNCrxsPXgOoP3C7qr6nqt+6y3vAHe42Y5LGE1Onk507jMz2lwTdntn+ErJzhzJ1WsXm3YkVC5/U0qlTJ7p3786iRYtK1rVv355169bRr18/7rvvPjp37kzfvn15+eWXgwZIJCdOnOC6667j/vvv57LLLiu17fjx47z44ovcfvvtJetatmxJTk4OrVu3pnbt2lx11VWsX7++wu8x1rwG0D4gWFeKo8CBIOuNSZgFCxaSnRv+novs3FuYP39hnGoUmYVPanrwwQeZNm0ap06dKlnXvHlzpk+fzueff86JEyf48ssvWbp0KQMGDCgpM2fOHF599dUy+8vJyUFVadu2LarKiBEj6Nu3L7fcckuZsrNnz+bSSy8lN/f0ONGXXHIJX331FYWFhRQXF/POO+947pWXCF5vRH0YmCIit6jqFwAi0hJ43N1mTNI4WLifBumtw5ZpkN6KoqKoTc5bKRY+kZ3ZpKGnm0zPbNIwDrU5bcCAAYwaNYpdu3Zx9tlnR3Xf7733Hi+88AKdO3fmpZdeAmDu3Ll06tQJgJo1azJ9+vRSz6lRowaTJk2id+/eqCpXXnklV199dVTrFU3ipeugiHwItMW59vOFu7olcBzY6l9WVSt2ojPJ5OTkaF5eXqKrYSogPSOTfmOW0Khpu5BlDu3byorpgziwv/yzj0ZTdQqfTZs2JfW3cXNauN+ViKxV1ZxovI7XFtDfovFixsTDkCE3s2q10wEhlC2r5zJ06M1xrFVZ1Sl8jAnGUwCp6u9iXRFjouXucWOYm5NLywsGBO2IsPfzNWxZPZ+/zlidgNo5LHyM8d4JoYSIPCUiGbGojDHRkJWVxaIFz/P27GGse+VhDu3bSvGpkxzat5V1rzzM27OHsWjB8xW6GTU/P5/RY8eTnpFJWloa6RmZjB47vlxdui18jHGUO4CAYUCjaFfEmGgaOHAg6/JW0/NcYcX0Qcy5uyUrpg+i57nCurzVDBw4sNz7jMbNrRY+xpzmqRNCqSeIfA1cpKqfx6ZKycE6IRh/+fn5dM3J5bu3zQt5Wu/t2cNYl7c6ZMuqPOGTjPPfVJZ1QkgdydYJwZhqrTw3t06b+kSZ7eVt+STj/DfGRFu5T8GpasOq3vqprMbpjRCRiEvjdDuTmSoqc3OrnXYzJjhPLSARuQH4RlX/EbD+WqCmqlo3bT/27bXqqejNrRY+xoTmtQX0EM5Np4GOuNuMqdLObJLB4YKdYcscLthF48anO4ha+BgTntcAao8zFXegLe62qBORkSLypogUiYiKSNuA7WeIyMsiskNEjovIbhGZ5w4R5F+ujYi8IiJHRGS/iEwTkVqxqLOpuoYMuZktq+eFLeN/c6uFjzGReQ2gQuCcIOvPBWJ1HqkesJzwLax/ATcC5wE/xAnDF30bRSQNeA1oiDOVxM3Aj3DGsDPGs7vHjWHL6nns/XxN0O2+m1vHjR1t4VONFBYWkpmZmZCpPW644QYefzy1P8q8BtA/gCdE5FzfChE5D5iMMzld1KnqFFV9BHg3xPZit8xqVd2uqu8DvwcuERHfPLlXAhcAt6jqOlV9A7gXuF1ErAeA8czrza3t27e38KlGJk2axFVXXVWq6/3evXsZN24cWVlZ1K5dm5YtWzJw4ECWLFlSUmbEiBFBBwnNy8tDRNi2bVvE1/7tb3/LxIkTOXjwYFTeSyJ4DaD7gIPARhHZKSI7gY+BQ8AvY1W58hCRdGAo8G9V9V2v6oEzdbj/yfvXgdpAt2i9duDd8Wk16/D6uJoUJN98Z6YSIt3cOmDAAAufGIvGSBTRcvToUZ599lluu+22knXbtm2ja9euvP766zzyyCN88MEHrFixgkGDBnHnnXdG9fU7depE+/btmTcv/KnhZOZ1LLhDQE8R6Qd0cVf/F/inlvdO1igTkT8Ao3FO2a0G/L9WNAP2BjxlP3DK3Ra4r5HASIA2bdp4ev2lS5dy05DhZOcOo9+YJTRIb83hgp18tvov/OXSP3Pt3GNkl/+me5OksrKymDb1iTL3+thpt9gL9be2avU85ubksmjB8xUa4aKilixZgojQs2fPknV33XUX4LRkGjRoULK+Q4cODBsWvht/RQwePJiFCxcyatSoqO87Hsp1I6p7CqvC87uKyATg1xGKXaGqb5Vjt48Cs4GzgQeBeSIysCLBqKqzgFngjIQQqbz/1M/+Nyg2atqOnGsepvUF1/CPW67nJ/8+Rnr5hx0zKSIW4ZOs898kSri/ta7X/IaWFwzgpiHhR6KItpUrV9KtW7eS33VBQQHLli1jwoQJpcLHp3HjxlGvQ/fu3ZkwYQLHjh2jbt26Ud9/rHk6BScifxGRXwRZ/3MRebYcrzcF6BBh+U859oeq7lfVz9xwvAlnivDL3c17gMyAp2QAae62SvFyd/y5uT8hb0bNyr6USVKxavkUFRxCVSMuqTIMT2Ul4zTr27dvp0WLFiWPt2zZgqp6Hm5o2bJlNGjQoNTSu3fvku07d+6kT58+dOzYkc6dO7N48eIy+2jRogUnT57kyy+/LFn32GOPccEFF3DhhRcm/ek5r9eABuL0OAv0L+Aqry/mhsUnEZajXvcXhO/91Hb/XQV0EJFWfmX6ASeAtZV4HcDb3fHn5f6UtX9KY0INyizV5dtrVWWn3eInkdOsz58/v1RIrFy5EoBjx45Rp06dknLlPenSu3dv1q9fX2pZsGBByfYaNWowZcoUNm7cyPLlyxk/fjxHjhwptQ9fq+fYsWMAfPjhhyxYsIC1a9eyZs0aZsyYQVFRUcj3kGheT8E1Bg4HWX8ESI9abfyISDOc6zS+nncdRaQxsENVC0SkB9AVp5dcEZAF/D9gG6d7zi3H6SzxvNuCOwvnlN0z7nWtSvF6d7wWnyz3f06T3Cx84iuR06wPHjyYSy+9tORxy5bOrYYZGRkUFhaWrD/nnHMQETZt2sT1118fcb/16tUjOzu71LqioqKSn5s3b07z5s0BaNasGRkZGRQUFFC/fv2SMgUFBQA0bdoUcAYR7dGjR0kwXnTRRSxbtizke0g0ry2gzwje0hmEczNqLNyJ09Fhvvv4NffxYPfxMZx7ev6Fc5PsbOADoJevF5yqnnLreBR4D3gB+DtwTzQqWJG7403qs/CJv0T+rTVs2JDs7OySxdfquPjii9m4cWNJufT0dPr378+MGTM4fLjs93X/cCmvtWvXcurUKVq3Lh3CH330ES1btiQz07nScOGFF/LWW29RVFREYWEhb731Fl988UXI95BoXgPoceD3IjJRRL7nLpOAScTopk5VfUhVJcgyx92+XlWvUNWzVLWOqrZT1f9R1V0B+9mhqleraj237FhVPRGNOpb37niT+ix8EiMZ/9b69+/Ppk2bOHDgQMm6J598ElUlJyeHxYsX8+mnn/LJJ5/w9NNP07lz5wq9TkFBAcOHD2fWrFlltq1cuZL+/fuXPO7YsSNjx46lb9++/OAHPyA3N5e0tLQKvW48eAogVX0OGA8Mx+kF9wZwC/BzVf1LzGqX5Mpzd7xJfRY+iZOMf2udOnWie/fuLFq0qGRd+/btWbduHf369eO+++6jc+fO9O3bl5dffjlogERy4sQJrrvuOu6//34uu+yyUtuOHz/Oiy++yO23315q/R133MG6det48803qVmzJuecE2wQmyThpaeN/wI0BZqW93mptnTr1k29WLJkiTZqnKFdB4zTGx/K059O26M3PpSnXQeM00aNM3TJkiWe9mOSW3FxsU6ePFn79OmjM2fO1OLi4kRXKeVs3LixUs9Pxr+1pUuX6rnnnqvffvtt1PddXFysN910kz744INBt8+YMUP79etXZv3evXtVVfWTTz7RTp066cmTJ8v92uF+V0CeRulztiIBdD/QOFoVSNbFawCpqm7ZskXHjB2v6Wdl6hlpaZp+VqaOGTtet2zZ4nkfJnlZ+ERHZQNINTn/1qZOnarbtm2L+n5XrlypIqIXXXRRyfLBBx+UbJ85c6Z+8sknZZ6Xm5urHTp0UHdW5wq9drwCqCJTch8CumgVn5TOpuQ2YKfdosmm5E4d8ZqSu9wzogL212eqBQsfY2KrIgFkTJVn4WNM7JVrLDhXR+DLiKWMSVEWPsbER7kDSEtPbWBMlWLhY0z8eB2MtFhEToVYjojIBhEZG+vKGhNLFj6xV95OTyb+4vk78toCGo0zNfaLwL/ddZcC1wF/AFrjjJSgqjo9ynU0JuYsfGKvZs2aHDt2jHr16iW6KiaMY8eOUbNmfEbw9xpA/YH/VdXZfuv+LCL/AQar6rUi8ikwBrAAMinFwic+vvOd7/DFF1/QsmVL6tata8c4yagqx44d44svvigZWy7WvAbQ94Ay8wEBbwPT3J/fACZHo1LGxIuFT/w0atQIgC+//JKTJ08muDYmmJo1a5KZmVnyu4o1rwF0AOd022MB66/DmeIaoAFwMCq1MiYOLHzir1GjRnH7cDPJz2sA/Q54RkT6cnrG0kuAKwHfSHj9cFpExiQ9Cx9jEs9TAKnqn0VkE841Ht98PJ/gzL2z2i0T2DoyJilZ+BiTHDzfB6Sqq3CmuDYmZVn4GJM8PAWQiISddltVC6JTHWNix8LHmOTitQW0Hwh3d1LyTrlnDBY+xiQjrwF0RcDjmsDFwP8AD0S1RsZEmYWPMcnJayeEYL3bVojI58DPgAVRrZUxUWLhY0zyKveEdKWeLJIFfKCq9aNXpeRgE9KlnsbpjThY+HXEcmc2aUhRwaE41MiYqieaE9JVZDoGXyUaAOMBGx3bJIWDhV/zwLeRy02oETmkjDGx57UX3NeU7oQgQD3gCDA0BvUyxhhTxZVnNGx/xcA+4N+qWhjdKhljjKkOPM0HpKrPBSxzVXWZhU/Vlp+fz+ix40nPyCQtLY30jExGjx1Pfn5+oqtmjKkCPAWQPxFpJiJt/JdYVMwk1tKlS+mak8uqzUK/MUsYMWU3/cYsYdVmoWtOLkuXLq3QfmMVajbRmTGpx+uMqGeKyHMicgz4AtgasESdiIwUkTdFpEhEVETahilbx52VVUUkJ2BbGxF5xZ25db+ITBORWrGoc1WRn5/PTUOG893b5tH1mt/QqGk7zkirQaOm7eh6zW/47m3zuGnI8HKHRqxCzdfV2hiTWry2gB4DLsKZfuE4MAT4JbAL+HFMauZ0cliOMxNrJI+5dSlFRNKA14CGQC/gZuBHwONRq2UV9MTU6WTnDiOz/SVBt2e2v4Ts3KFMnTbD8z5jFWr+9/kYY1KL1wAaCIxR1deBU8BaVZ0M3A/cEYuKqeoUVX0EeDdcORG5FmekhnuCbL4SuAC4RVXXqeobwL3A7SJik5KEsGDBQrJzh4Utk517C/PnL/S8z1iEWuBNpsaY1OI1gBoD292fDwJnuT+vAi6Lcp08E5FWwNM4LbJjQYr0ADapqv+9Sq8DtYFusa9hajpYuJ8G6a3DlmmQ3oqiov1hy/iLdqgFG+HAGJNavHbDzgfaAzuATcBNIvIf4AdAQkbCdk+vzQceV9UNIa4RNQP2Bqzbj9OKaxZknyOBkQBt2lTfvhVnNsngcMFOGjVtF7LM4YJdNG6c4Xmf0Qy1UMPrnNmkoaebTM9s0tBzvY0xseO1BTQH6Oz+/Huc027fAI8Cf/D6YiIywe0oEG7p43F3v3LrMNnr60eiqrNUNUdVc5o2bRqt3aacIUNuZsvqeWHLbFk9l6FDvZ/28oVaOF5CLdzYbkUFh1DViIsNw2NMcvA6GOkTfj//S0TOB3KAzar6YTlebwoQ/pPNaWV58T2cjgUnAwaXXC0iL6jqUGAP0DPgeRk400fs8fg61c7d48YwNyeXlhcMCHrNZu/na9iyej5/nbHa8z6HDLmZVaudDgihRAo1G1jUmKolbACJyGTgReA9VS32rVfVHXgPihKquh/nFFg0/ATwHwS1Bc71naHAe+66VcADItJKVX295PoBJ4C1UapHlZOVlcWiBc9z05BhZOcOJTv3Fhqkt+JwwS62rJ7LltXzWbTgebKysjzvs7KhZuFjTNUTqQVUF1gE1BKR14CXgNdVNdgF/6gSkWY412nOdVd1FJHGwA5VLVDVrQHlD7s/5vuFzXLgY+B5EfkFTueJR4FnVNXOw4QxcOBA1uWtZuq0GcyfPoiiov00bpzB0KE389cZq8sVPlC5ULPwMaaK8nLOHOgOTAQ+whmA9GXgp0BTL8+vyIJz/48GWUaEKN/W3Z4TsL4N8CpwFDgATANqR3r9bt26qYm+LVu26Jix4zX9rEw9Iy1N08/K1DFjx+uWLVuCli8uLtbJkydrnz59dObMmVpcXBznGhtj/AF5GqXP+XLPByQi2cC17nIpsAanZbRQVb8ofwQmJ5sPKPHUWj7GJJ1ozgdU7rHgVHWLqj6uqr2BlsCfgctxRhkwJiosfIyp+iJ1QjgFNFfVr4JtV6dTwZ/dxZiosPAxpnqI1AKyv3oTVxY+xlQfFZ6S25iKapzeiIOFkUcs+O+GtYwcOTIONTLGJIKXALpRRMJ2WVbV56NUH1MNHCz8mge+jVzOy7A6xpjU5SWAfo/TvTkUBSyAjDHGlIuXAGofqhOCMcYYU1GROiHYPMfGGGNiwnrBGWNMELt37+b8c7LZs8fGLY6VSAH0HMEnejPGmCrtjxMn8tWunfxx4oREV6XKihRAtwF3ish7IvIfEZkkInXiUTFTNZV36CdjEmH37t0899wc/tm/G889N8daQTESKYDuw+kFdwTYDfwcZzBPY8rNd5OpMcnujxMncmtWcy7OaMTw9s2tFRQjkQJoBDBGVa9U1WuB64DhYremm3LyH+HAmGTma/3c27EVAPd2bGWtoBiJ1A37bJypDHxex+mY0AKoMiNfm9gKHF7nvxvWerrJ9MwmDeNQO2NK87V+mterDUDzerVLWkGTp89IcO2qlkgtoFr4dUJw54L4Bqgdy0qZqiPY2G5FBYc8zRVSVGBzBpr4Cmz9+FgrKDa83Ij6iIgc9XtcC3hQRA76Vqjq2KjXzKQ8G1jUpJrA1o+PtYJiI+yEdCLyFh5uRlXVK6JYp6RgE9JVjoWPSTW7d+/mgnPP4eNrLikTQAC7j57gwlfX8PFnW2jWrFkCapgcojkhXdgWkKr2icaLmOrFwsekolCtHx9rBUVfuafkLvVkkXbAbar6QPSqlBysBVQxFj4mVbVv04atO3dGLNeudWs+37EjDjVKTnFrAYV48VrAD3FuUr0C2AFUuQAy5WfhY1JZdQ6VRInUC66EiHQSkWnAHmAekA/0UNV2saqcSR0WPsaY8gobQCLSUETuEJE1wPtAI+DHQDEwVVX/E4c6miRn4WOMqYhIp+C+AFYBU4H/U9WjgH24mBIWPsaYiop0Cu4EcKa72CCkphQLH2NMZUQKoBbA48A1wBci8rKI3IBNVFftWfgYYyorbACp6klVXayqA4DzgbXAozin7n4nIoNEpGYsKiYiI0XkTREpEhEVkbZBymxzt/kvvw8o00ZEXhGRIyKyX0SmuT35TAVZ+BhTcTbR3Wmee8Gp6nZV/R3QDrgKpxX0d2BfjOpWD1gOPBSh3MNAc7+lZNx0EUkDXgMaAr2Am4Ef4bTqTAVY+BhTOTbR3WmeA8hHHctU9UagNfC76FcLVHWKqj4CvBuh6NequsdvOey37UrgAuAWVV2nqm8A9wK3i0ijWNS7KrPwMaZybKK70sodQFByWqsjsF9Vn4hyncrrHhE5ICLrReTXAafXegCbVNX/9ubXcUbz7hbXWqY4Cx9jKs8muist0n1APxaR/wlY9zSwFfgQ+FBEWsawfpFMwzmtdgUwA7gbeMpvezNgb8Bz9gOn3G2luNed8kQkb9++WJ1ZTD0WPsZUnk10V1akFtAYnJtOARCR7wN3AL8FbgDSgN94fTERmRCk00Dg0sfr/lR1sqq+qaofqOqzwF3AbSJyltd9BOxvlqrmqGpO06ZNK7KLKsfCx5joCDfRXXUVKYDOA1b7Pb4WWK6qE1X1/4Bf4Fxn8WoK0CHCUpnRFf7t/pvt/rsHyAwok4ETnFX6a0c0etpY+BgTHTbRXXCRAqgBUOj3+DLgX36PPybIqaxQVHW/qn4SYTkaeU8hdXH/3e3+uwroICL+v/V+ODfYrq3E6yS9yva0sfAxJnq8THRXHUUKoF04vchwe411At7z234WcDjI8ypNRJqJSBfgXHdVRxHpIiLp7vYeInK3u66diNyIc/3nZVX1DWu7HCcknxeRi91TiI8Cz6hqlZ3vubI9bSx8jImeUK0fn+rcCooUQIuBaSLyU+BZnJaF/ym5HOCTGNXtTuC/wHz38Wvu48Hu4xM4A6O+BWzEuR/oGZxOCQCo6ilgEHAUJzhfwLl36Z4Y1TkpVKanjYWPMdFVnonuqptIU3LXBWbiDMWzBxipqiv9tr8JLFPVP8S6ovGWqhPSBU4rXJ5phC18jIm+qjbRXTyn5D4GDA+z/YpoVMJET7ieNuGmEbbwMSY2UiFUEqVCN6Ka5FTRnjYWPsaYRAjbAhKRl73sRFUHRy5lYs1LT5vAVpCFjzEmUSK1gA54XEyCVaSnjYWPMSaRIl0D+km8KmIqpzw9bSZPn2HhY4xJuLC94KqzVOsFV56eNvnbt1v4GGMqJG694Ezq8NrTxlo+xphkYb3gqhELH2NMMrEAqiYsfIwxycYCqBqw8DHGJCMLoCrOwsckUjSmBTFVlwVQFWbhY6KhMiFS2WlBTNVmAVRFWfiYaKloiFR2WhBT9VkAVUEWPiZaKhMilZkWxFQPFkBVjIWPiaaKhkjg0FDVedI1E5oFUBVi4WOiqTIhEm5aEGN8LICqCAsfUxE7d+6kT58+dOzYkc6dO7N48eKSbRUNkYpOC2KqHxsLLoRUGgvOwsdU1O7du9m7dy9dunRhz549dOvWjc8++4xDhw6Vmlm3pLyHGXbvHj0a3n2dJ3Kyym7Ly0d69Q87OaJJbtEcC85aQCnOwsdURvPmzenSpQsAzZo1IyMjg4KCAk9zSwVTkWlBTPVlAZTCLHxMNK1du5ZTp05Ro0aNCodIeaYFCfTUU0/Rrl076tSpQ7du3Vi5cqWnekd63pNPPknnzp1p1KgRjRo1okePHrz22mvlLmNiQFVtCbJ069ZNk1lxcbFOnjxZ+/TpozNnztTi4uJEV8mkoJMnT6qq6oEDB7Rjx4763nvv6fhRo3T8Rdmqt/UPuYy/KFvvHj2qzP7atW6tQMSlXevWpZ63aNEirVGjhs6aNUs3btyoo0eP1vr16+v27dvD1t/L81566SVdsmSJbt68WT/99FP91a9+pTVq1NANGzaUq4xxAHkapc/ZhH/QJ+uSzAFk4WMqYufOnQrowoUL9YorrtDatWvr7Nmz9fjx49qrVy99/vnnVbXiIVIZ3bt315/97Gel1mVnZ+v9998fk+c1adJE//SnP1W6THUUzQCyU3ApRtVOu5mK2bBhAwB/+MMfuOeee/j444+59tprGTFiBH379uWWW24BnLmlAj8oJk6cSP369UstXxUU0KBBg5LF6ymzQN988w1r167lyiuvLLX+yiuv5P3334/q806dOsWiRYs4fPgwl112WYXLmOiwCelSiIWPqYz169dTp04dFi9eTHZ2NgDvvvsuL7zwAp07d+all14CYO7cuXTq1KnUc++8805uvPHGsPtv2bJlheq1f/9+Tp06RWZmZqn1mZmZrFixIirP+/DDD+nRowfHjx+nQYMGvPjii2Xeo5cyJrosgFKEhY+prA0bNnDVVVeVhA/A5ZdfTnFxccTnpqenk56eXuk6zJ8/nzvuuKPk8dKlS8nKKttdO9rOO+881q9fz8GDB/nb3/7GrbfeyltvvcWFF15YrjImuiyAUoCFj4mG9evXM3r06Ao9d9KkSUyaNClsmaVLl9KrV6+wZQYPHsyll15a8rhly5akpaWRlpbG3r17S5Xdu3dvyHuNADIyMjw/r1atWiXB261bN9asWcMTTzzB7Nmzy1XGRFfSBpCIjARuBi4GzgTaqeq2IOX6Aw8BFwHfAOtUta/f9jbAk0Bf4BiwALhHVb+J8VuICgsfEw1HjhwhPz+frl27Vuj50ToF17BhQxo2bFhmfbdu3XjjjTe44YYbSta98cYb/PCHPwy5r1q1alXoeQDFxcWcOHGi0mVMJUWrN0O0F2A88L/uvwq0DVLmOqAQuAs4D+gADPXbngZ8CLwFdAX6AV8C0yO9fjL0grPebiZa3n//fT3jjDP066+/TnRVglq0aJHWrFlTn3nmGd24caOOHTtW69evr9u2bSspM336dD3vvPPK/bz77rtP33nnHd26dat+8MEHev/996uI6JIlS8pVxjioTt2wgZxgAeSGyw7g9jDPHQgUA6391g0DjgONwr1uogPIwsdE09NPP13mwzvZPPnkk3r22WdrrVq1tGvXrvr222+X2v7ggw+q8525fM+79dZbtU2bNlqrVi1t2rSpfu9739Nly5aVu4xxRDOAkn4sOBHJAdYQcApORLoD/wZGAOOAFsAHwH2q+l+3zMPAD1X1Ar/nNQW+Avqq6puhXjeRY8Gp2mk3Y0xysrHgHO3dfx8GJgGDgF3AWyLS3N3WDNgb8Lz9wCl3WykiMlJE8kQkb9++fbGpdQQWPsaY6iKuASQiE0REIyx9PO7OV/eJqvo3VV0LjAQOAsMrUj9VnaWqOaqa07Rp04rsolIsfIwx1Um8e8FNAeZFKLPD4752u/9u9K1Q1W9FZDPQxl21B+gZ8LwMnOtHSTUcr4WPMaa6iWsAqep+nFNg0bAWOIHT++1dABE5A8gCXnfLrAIeEJFWqrrLXdfPfd7aKNWj0ix8jDHVUTLfB9QM5zrNue6qjiLSGNihqgWqekhE/gT8TkR2AduA0UATYK77nOXAx8DzIvIL4CzgUeAZVT0UtzcThoWPMaa6StoAAu4EHvR77Juc4yfAHPfnX+LcfPocUA9YB1yhqrsBVPWUiAwCngLew7kRdb77vISz8DHGVGdJ3w07UWLdDdvCxxiTiqwbdoqz8DHGGAuguLPwMcYYhwVQHFn4GGPMaRZAcWLhY4wxpVkAxYGFjzHGlGUBFGMWPsYYE5wFUAxZ+BhjTGgWQDFi4WOMMeFZAMWAhY8xxkRmARRlFj7GGOONBVAUWfgYY4x3FkBRYuFjjDHlYwEUBRY+xhhTfhZAlWThY4wxFWMBVAkWPsYYU3EWQBVk4WOMMZVjAVQBFj7GGFN5FkDlZOFjjDHRYQFUDhY+xhgTPRZAHln4GGNMdFkAeWDhY4wx0WcBFIGFjzHGxIYFUBgWPsYYEzsWQGFY+BhjTOxYAIWwd+9eCx9jjIkhC6AQioqKLHyMMSaGRFUTXYekJCL7gO0JrkYGsD/BdUh2dowis2MUmR2jyHzH6GxVbRqNHVoAJTERyVPVnETXI5nZMYrMjlFkdowii8UxslNwxhhjEsICyBhjTEJYACW3WYmuQAqwYxSZHaPI7BhFFvVjZNeAjDHGJIS1gIwxxiSEBZAxxpiEsABKAiIyUkTeFJEiEVERaRuiXH8RWSUiR92y/wrY3kZEXhGRIyKyX0SmiUituLyJGPN6jNyydURkg1suJ2BblTxGkY6PiLQVkdki8rmIHHP/fURE6gaUq5LHB7z9HxKRJiIyV0QOustcEWkcUKaTiLztHscvROS3UoXvVheRZu5x2ON+9mwQkaEBZSIet2AsgJJDPWA58FCoAiJyHbAImAtcDPQAZvttTwNeAxoCvYCbgR8Bj8eozvEW8Rj5eQzYFbiyih+jSMfnfCAN+B/gAmAMMByY6itQxY8PePs/tADoCgxwl644f3MAiEgj4A1gL3AJMA74JfDzmNQ4OTwPdACuBS50H88Vkd5+ZcIet5BU1ZYkWYAcQIG2AevTgB3A7WGeOxAoBlr7rRsGHAcaJfq9xfoY+W2/FvjY/YNRIKc6HaNIxyeg7F3Agep0fMIdI7//Mz391l3urjvPffw/wCGgrl+ZB4AvcDt1VbUFOAz8JGDdduAer8ct1GItoNTQDWgNfCMi69ym8HIRudivTA9gk6ru9Fv3OlDbfX6VJyKtgKeBIcCxIEWq/TEK0Ago9Htc3Y9PD5wP2/f91r0HHAEu8yuzUlX9/3+9DrQA2sahjonwLnCjiJwlImeIyLVAU2CFu93LcQvKAig1tHf/fRiYBAzCOcX0log0d7c1wzkt4G8/cMrdVqW5p4/mA4+r6oYQxar1MfInImcD9wBP+a2u7senGbBP3a/wAO7PX3H6/Qc7Rnv9tlVFN+K0ZvYDJ3D+zm5W1fXudi/HLSgLoBgRkQnuhc5wSx+Pu/P9niaq6t9UdS0wEjiIcx4/JUX5GP0K+AaYHLMKx1mUj4//fjOBZTjXMp6IcrXjKlbHqKor53GbgDMQ6fdxTl8+CjwvIhdVth41KrsDE9IUYF6EMjs87mu3++9G3wpV/VZENgNt3FV7gJ4Bz8vAuX60x+PrxNsUoneMvodz4fxkQIek1SLygqoOJfWO0RSid3wAp0cT8C/gI+AW/2+tpN7xgegeoz1AUxER33Fxe7d9h9Pvfw+QGfC8TL9tqWIKHo6biGThdFjp4ndmYYOI9HLX/wxvxy0oC6AYUdX9RG9497U4Td/zcM7HIiJnAFk4558BVgEPiEgrVfX1AOvnPm9tlOoRVVE+Rj8B6vs9boFzbIbinI+GFDtGUT4+uKdr38TppHGzqn4bUCSljg9E/RitAhrgXNPwXc/ogfP/6n2/Mn8QkTqqetxd1w/4EtgWpXrEnNfjJiL13B9PBWw6xekzM16OW8iK2JL4XibNgC44F88VuMp9nO5XZgrOdZ/+OEE0HecUXHN3exrwIc6324txmstfANMT/f7idYwCyrelbC+4KnuMIh0fnED+DHgLp0NLM78lraofH6//h4Cl7jHo4S4fAq/4bT8T51v9IpwuyT/A6RX3i0S/vxgds5rAZuAdoDvOl95f4PSWvMbrcQu5/0S/QVsUnPsSNMgyIuA/wh/d//yH3A+SrgH7aQO8ChwFDgDTgNqJfn/xOkYB5csEUFU+RpGODzAixPZS3ZGr6vHx+n8IaIJzauqQu8wDGgfsp5P7gXwc5/T4g1TRLtju+z0H+DtOZ4sjwAbg1oAyEY9bsMUGIzXGGJMQ1gvOGGNMQlgAGWOMSQgLIGOMMQlhAWSMMSYhLICMMcYkhAWQMcaYhLAAMsYYkxAWQMaUg4i8KiJzylH+DBGZKSIHquLAmCLylt/glbl+62bE+HXniMirsXwNr6/tPvYdgx8lok6pygLIeBbwh6biTNn8qoicH1CuoYhMEZHt4kxb/L6IXBJQZluIEXiXisg7IvJKQPl6IvKpiPwpSF1OishX4ky3PEpEagY8d5SIfCAih9xllYgMCvEe/yQiU6JywBxX4YxTdw3QHHhfRP5XRNa4ddknzhTYFwbUo6mIPOUepxMisldE/iki/dztL4vIP0O8hw7ucbkyyHH6XEQeE5H6fuU9H8sQ/uK+t4SPFxeD358X43DevyknCyBTXitw/tiaA1cCdYEXA8o8izNm3a04w5YsB1aISEu/Mpf47ac5zhS+ijPG1gjgChH5qV/5P+CMVfaLIHVp69blFeB3wEr/D1icMfTuc18jB2ess5dEpLN/pd0RfAcDL3k4Dl5lA7tV9X1V3aOq3wB9cObhuQzoC3yLc3zS/Z73d5yxt24DzgWuxhlv6yx3+2ycY9Q2yGvehjNjpW/CMN9xao8ze+ddONOW+/N6LIM56r63kxHKxVSMfn8RqepBVU2lkbCTR6LHGbIl8QvOnDB5wBlBtuUBU9yf5wCvBmy/Gic46rqP6+J8oF4bUG4tMCFMHX4NFPnt506cwVbPxplq4Vvgcr/yZerirr8QZ16g30V4zwXAHQHruuOMf1bDfVzPfZ3DOONg/QpnnLQ57nYB7gXycWZg/RAYFlBH/zHHtoWoSwOc0YWvcR83dst/P0z9a+CMQ/a7gPU13br+Nszv7BmcUKz0scQZk3BGuHXu768IuNNvneB8mdiMM9r2LuARd9sAYCXObK0FOKOadwh4jWDvq+T359bhaeBxdx/7cFoqtYEn3frswJmSwn8ftXEG/t2LM9bbav//dxGOlwI/SvTfcyot1gKq5kTkPJx5PX6pqsVBimzCGRk52HMbAj8GPtTTUxTXwGmpHA8ofgxnnvhg+xGcb+3zfPtR1T/h/PHPxTnFM1lV3430flT1I5zJ1n4Y4rXSROQmnA/9wKHirwNe09PTFDyGM9T+D3E+RC8GevuVn+DWexTQEXgEmOl3em8cziy2u3BaF6VOQ/ppiHM2wjc99mF3GSwidUK8z2+B54AR7tQcPtfgzOHzlxCvBc7vIuKptUjH0gv3msiLwEj3d+ozCfgNzjG7ALgB8E0FXh8nBLrjtBYPAq+ISK0IL3cdpX9/Q4GvgUuB37v7fAlnVPAcnOP3rJyeVRicAX9/DPwU5/f9IbAsoIyJlkQnoC2JXXC+Da8Ns30q8Jn78xyclojvA1JxvkVeGPCc93G+wbbECaNhON/wPw3xGle6+7ooYH07nGHfNxMwIjMhvoW6236Pc1rIf10nt87f4nz7HRTkeR8DP3B/boDzzXyo3/YG7nPn4HxIHgN6BexjCrDE7/E9hGj5+JX5K/Bf3GkR3HU/xPnmfhxnvpXHgEsDnneOe9yu9Fv3GrA01HHC+VDfD7xQ0WMZsP0tQrSAOD1r75UB2xu47+vOUPsNKF/f/f8TtgUc8Pt7C1jlt01wWkEv+62ridPC+5Hf63wDDPcrk4bTwp0Q7rXd9dYCKudiLaBqzP3m/CPgb37rnhCRMX7FGuIMwe7zDs4cKl1wPsz+CSwXkdZ+ZW7BCY5dOB/iY4GF7rpgbgfW6OkZF31+ivMh3wonjDy/NZwPA3+funW+FOfUzHP+F/5FJBvnGolvgr8soBbOhz8AqnoY5xsxOC2eOjjfjg/7FuB/3Od6q6jIZJyW4Q9VtWTSL1X9O84cPtfgXPu5DGd211/5ldkMvI1znBCRFjjX3mYHvMwAt36+MHsHp9XrqYqUPZZeXIdzqmuAqi4P2NYR51RXqE4UWSKyQETyReQQzumwMzg9+2+w5wT+/gA+8P2gTkJ8xenfH+pcsyrEmbkTnN9bTU5PYIj7O1nl1tlEmQVQ9dYO53rDh37rbsSZC8bnIvymAsf5NrzFXdbgTMnbCOfbLgCqmq+q38X5pttaVbvj/GF/HlgBEfkOcC1OS8x//SXA/TgB+QZOYKR5fF8dA19LVb9x67xWVf8XWA/c7VfkOuCfquoftuH4/nau4XQgd8E5nXSllx2IyBPAzUBfVS1zbFT1uKq+oaoPq+plOMHyUMCpqNnAdW4HhhE4raZ/BOzK96XhPKCOqv5AVb/yUkeCHEuPNuBco7rNPcVaHq8CTYE7cL4wXIzTcg13Cu46yv7+AjtFaIh1Xj4Hbd6aGLAAqt6auP8eBhDnHpUWOKchEJFzcD64Anu5+VOclk29MhtUj6jqbhFpgvPNPPCDEZwPzRM4LSTc160DPI9zsX8pTrhl41zwD8tt1QzAr1UXwhk438J9rqV076l8nA+rXL9918e5MA9OKJ8AzvYLZN+y3UM9p3I6fD6JVN7vNWvgtLx8/oZzOmsYTkvoeS3bG833pWF7kG3h6uj1WAazFef6zZXArIAQ2oRz7L4X5DXPAs4HJqnqClXdhNMKrxHh9QJ/fxWRj/N/v6dffdJwZvjcGOpJpuIi/VJN1bYDJzyGiEgRzuyXrwBXi8gGnK7CH1A6gGqLSDP35ybAaJyWTsl9OyLSH+cD/hOc4HjU/bnUhXH3Q+lnwCL39JbPIzgfsj8HUNU9IjIKpxX0sqp+HFCXM3C+MX8Pp6faWvy6GYvI73GujezE+TAbgvPhOMjd3hQnaEpuIlTVwyIyG/iDiOwDvgR+i3NNAFX9WkQeAx5z38c77nHIBYpVdVbQI+683pM4pymvAwr9judh93XPAhYDf8Y5/l/jXDS/F+db/iG/eh4TkQU4s302oezpN688HcvyUNXPReQKnOsxM0XkDnV87QbwIyJyAufYnQV0A2biXKO6XUR24lxHfBSnBRRUsN9fBet7RESexvmd78cJ0buBTJy/BRNtib4IZUtiF5wPtUKc8+OjcO45+QznNNw/gOZ+ZedQulvxIeA/ONcv/Pd5I863yRM4p2FmAGcGee0r3P1091vXG+fDpk+Q8otxuoXXCKjLtzgfWm/hBGKtgOfNwbkv5oT7PlcA/f22/xR4P8jr1cdpiR12n/cbynbDHsPp1tA+nNOF/fz2UaYTQsAx9F8ecrfXxukltsb93RzF6YgxGUgPUk/fPVTvBdk2hxAdDEL8XsMeyyDPfYvI3bCzcMJ/FpTMwnwGzinWz3FaHTuBie62vsBHOC27j3Baz4cpPXV2yfsK9vsLUa+PfMfYb90eYLTfY/9u2CewbtgxXWxKblPticg/cD68/5jouqQaEXkL+EhVRyewDknx+xMRBW5Q1YqcsqyW7BqQMU6vp4URS5lQRrq97ELd5xRrCf39ucP/HI5c0gSyFpAxpsLc4ZXqug93quqJRNYnEdyenI3ch7vVe0/Kas8CyBhjTELYKThjjDEJYQFkjDEmISyAjDHGJIQFkDHGmISwADLGGJMQFkDGGGMSwgLIGGNMQvx/mMM7M4wjDeEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams.update({'font.size': 14})\n",
    "\n",
    "colors=['tomato','lawngreen','cornflowerblue']\n",
    "legend=['-(CH$_2$)$_3$-','-(CH$_2$)$_6$-',\"-(CH$_2$)$_9$-\"]\n",
    "markers=['^','s','o']\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6,6))\n",
    "\n",
    "lb=np.array([0,8,17],dtype=int)\n",
    "ub=np.array([8,17,26],dtype=int)\n",
    "\n",
    "for i in range(3):\n",
    "    xs=benchmark26_NN3.wB97XD3_def2SVPD.values[lb[i]:ub[i]]\n",
    "    ys=benchmark26_ISAPT.ISAPT.values[lb[i]:ub[i]]\n",
    "    ax.scatter(xs,ys,c=colors[i],marker=markers[i],edgecolors='black',s=100,alpha=1.0,linewidths=1.0)\n",
    "\n",
    "plt.legend(legend)   \n",
    "\n",
    "lims = [\n",
    "    np.min([ax.get_xlim(), ax.get_ylim()]),  # min of both axes\n",
    "    np.max([ax.get_xlim(), ax.get_ylim()]),  # max of both axes\n",
    "]\n",
    "\n",
    "# now plot both limits against eachother\n",
    "ax.plot(lims, lims, 'k-', alpha=0.75, zorder=0)\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlim(lims)\n",
    "ax.set_ylim(lims)\n",
    "\n",
    "linreg = sc.stats.linregress(benchmark26_NN3.wB97XD3_def2SVPD.values, benchmark26_ISAPT.ISAPT.values)\n",
    "plt.text(-120,-160, r'$r^2=$ {0:5.3f}'.format(linreg.rvalue), color='k')\n",
    "\n",
    "plt.xlabel('$\\omega$B97XD3/def2SVPD [kcal/mol]')\n",
    "plt.ylabel('ISAPT/aug-cc-pVTZ [kcal/mol]')\n",
    "plt.savefig('ISAPT_vs_wB97XD2_NN3.png')  \n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Breakdown of benchmark results with terminal group -CH$_2$CH$_2$CH$_3$ for L1, L2 and L3 linkers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "ISAPT=pd.read_csv('../csv_files/benchmark26_ISAPT_Psi4_1.1.csv',header=0)\n",
    "#ISAPT=pd.read_csv('../csv_files/benchmark26_ISAPT_Psi4_1.4_junccpvdz.csv',header=0)\n",
    "NN=pd.read_csv('../csv_files/benchmark26_term_group_1.csv',header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0\n",
    "f=8\n",
    "mDLPNO=df(ISAPT.values[i:f,0]-NN.DLPNOCCSDT.values[i:f],columns=[['DLPNO']])\n",
    "mDFT_SVP=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVP.values[i:f],columns=[['wB97XD3_SVP']])\n",
    "mDFT_TZVP=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVP.values[i:f],columns=[['wB97XD3_TZVP']])\n",
    "mDFT_SVPD=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVPD.values[i:f],columns=[['wB97XD3_SVPD']])\n",
    "mDFT_TZVPD=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVPD.values[i:f],columns=[['wB97XD3_TZVPD']])\n",
    "\n",
    "benchmark26=pd.concat([mDLPNO,mDFT_SVP,mDFT_TZVP,mDFT_SVPD,mDFT_TZVPD], axis=1)\n",
    "\n",
    "pDLPNO=df(np.abs(100*(ISAPT.values[i:f,0]-NN.DLPNOCCSDT.values[i:f])/ISAPT.values[i:f,0]),columns=[['DLPNO']])\n",
    "pDFT_SVP=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVP.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_SVP']])\n",
    "pDFT_TZVP=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVP.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_TZVP']])\n",
    "pDFT_SVPD=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVPD.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_SVPD']])\n",
    "pDFT_TZVPD=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVPD.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_TZVPD']])\n",
    "\n",
    "benchmark26p=pd.concat([pDLPNO,pDFT_SVP,pDFT_TZVP,pDFT_SVPD,pDFT_TZVPD], axis=1)\n",
    "\n",
    "mad=np.mean(benchmark26.abs())\n",
    "std=benchmark26.std(ddof=0)\n",
    "mpad=np.mean(benchmark26p)\n",
    "metrics=df()\n",
    "metrics['MAD-L1']=mad\n",
    "metrics['STD-L1']=std\n",
    "metrics['MPAD-L1']=mpad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0+8\n",
    "f=8+9\n",
    "mDLPNO=df(ISAPT.values[i:f,0]-NN.DLPNOCCSDT.values[i:f],columns=[['DLPNO']])\n",
    "mDFT_SVP=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVP.values[i:f],columns=[['wB97XD3_SVP']])\n",
    "mDFT_TZVP=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVP.values[i:f],columns=[['wB97XD3_TZVP']])\n",
    "mDFT_SVPD=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVPD.values[i:f],columns=[['wB97XD3_SVPD']])\n",
    "mDFT_TZVPD=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVPD.values[i:f],columns=[['wB97XD3_TZVPD']])\n",
    "\n",
    "benchmark26=pd.concat([mDLPNO,mDFT_SVP,mDFT_TZVP,mDFT_SVPD,mDFT_TZVPD], axis=1)\n",
    "\n",
    "pDLPNO=df(np.abs(100*(ISAPT.values[i:f,0]-NN.DLPNOCCSDT.values[i:f])/ISAPT.values[i:f,0]),columns=[['DLPNO']])\n",
    "pDFT_SVP=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVP.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_SVP']])\n",
    "pDFT_TZVP=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVP.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_TZVP']])\n",
    "pDFT_SVPD=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVPD.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_SVPD']])\n",
    "pDFT_TZVPD=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVPD.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_TZVPD']])\n",
    "\n",
    "benchmark26p=pd.concat([pDLPNO,pDFT_SVP,pDFT_TZVP,pDFT_SVPD,pDFT_TZVPD], axis=1)\n",
    "\n",
    "mad=np.mean(benchmark26.abs())\n",
    "std=benchmark26.std(ddof=0)\n",
    "mpad=np.mean(benchmark26p)\n",
    "#metrics_TZVPD2=df()\n",
    "metrics['MAD-L2']=mad\n",
    "metrics['STD-L2']=std\n",
    "metrics['MPAD-L2']=mpad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=0+8+9\n",
    "f=8+9+9\n",
    "mDLPNO=df(ISAPT.values[i:f,0]-NN.DLPNOCCSDT.values[i:f],columns=[['DLPNO']])\n",
    "mDFT_SVP=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVP.values[i:f],columns=[['wB97XD3_SVP']])\n",
    "mDFT_TZVP=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVP.values[i:f],columns=[['wB97XD3_TZVP']])\n",
    "mDFT_SVPD=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVPD.values[i:f],columns=[['wB97XD3_SVPD']])\n",
    "mDFT_TZVPD=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVPD.values[i:f],columns=[['wB97XD3_TZVPD']])\n",
    "\n",
    "benchmark26=pd.concat([mDLPNO,mDFT_SVP,mDFT_TZVP,mDFT_SVPD,mDFT_TZVPD], axis=1)\n",
    "\n",
    "pDLPNO=df(np.abs(100*(ISAPT.values[i:f,0]-NN.DLPNOCCSDT.values[i:f])/ISAPT.values[i:f,0]),columns=[['DLPNO']])\n",
    "pDFT_SVP=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVP.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_SVP']])\n",
    "pDFT_TZVP=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVP.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_TZVP']])\n",
    "pDFT_SVPD=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVPD.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_SVPD']])\n",
    "pDFT_TZVPD=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVPD.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_TZVPD']])\n",
    "\n",
    "benchmark26p=pd.concat([pDLPNO,pDFT_SVP,pDFT_TZVP,pDFT_SVPD,pDFT_TZVPD], axis=1)\n",
    "\n",
    "mad=np.mean(benchmark26.abs())\n",
    "std=benchmark26.std(ddof=0)\n",
    "mpad=np.mean(benchmark26p)\n",
    "#metrics=df()\n",
    "metrics['MAD-L3']=mad\n",
    "metrics['STD-L3']=std\n",
    "metrics['MPAD-L3']=mpad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                  MAD-L1     STD-L1    MPAD-L1     MAD-L2    STD-L2  \\\n",
      "DLPNO          16.082499  10.334850  14.872596   9.523370  6.345837   \n",
      "wB97XD3_SVP    24.359147  11.090008  21.964101  20.043715  9.277303   \n",
      "wB97XD3_TZVP   18.159709  10.497811  16.626754  12.375565  7.342746   \n",
      "wB97XD3_SVPD   15.965609  10.731412  14.716597   9.781178  6.809893   \n",
      "wB97XD3_TZVPD  16.071649  10.575202  14.912407   9.745563  6.503269   \n",
      "\n",
      "                 MPAD-L2     MAD-L3    STD-L3    MPAD-L3    MAD-all   STD-all  \\\n",
      "DLPNO           7.075165   9.741644  5.482192   6.893689  11.617120  8.112891   \n",
      "wB97XD3_SVP    14.892804  21.677590  8.986185  15.435650  21.937112  9.932630   \n",
      "wB97XD3_TZVP    9.160682  13.362595  6.765574   9.457714  14.496966  8.634106   \n",
      "wB97XD3_SVPD    7.254612  10.186142  5.823995   7.185915  11.824260  8.395911   \n",
      "wB97XD3_TZVPD   7.269574  10.139866  5.389075   7.191606  11.828540  8.179234   \n",
      "\n",
      "                MPAD-all  \n",
      "DLPNO           9.411556  \n",
      "wB97XD3_SVP    17.256496  \n",
      "wB97XD3_TZVP   11.560754  \n",
      "wB97XD3_SVPD    9.526828  \n",
      "wB97XD3_TZVPD   9.594226  \n"
     ]
    }
   ],
   "source": [
    "i=0\n",
    "f=8+9+9\n",
    "mDLPNO=df(ISAPT.values[i:f,0]-NN.DLPNOCCSDT.values[i:f],columns=[['DLPNO']])\n",
    "mDFT_SVP=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVP.values[i:f],columns=[['wB97XD3_SVP']])\n",
    "mDFT_TZVP=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVP.values[i:f],columns=[['wB97XD3_TZVP']])\n",
    "mDFT_SVPD=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVPD.values[i:f],columns=[['wB97XD3_SVPD']])\n",
    "mDFT_TZVPD=df(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVPD.values[i:f],columns=[['wB97XD3_TZVPD']])\n",
    "\n",
    "benchmark26=pd.concat([mDLPNO,mDFT_SVP,mDFT_TZVP,mDFT_SVPD,mDFT_TZVPD], axis=1)\n",
    "\n",
    "pDLPNO=df(np.abs(100*(ISAPT.values[i:f,0]-NN.DLPNOCCSDT.values[i:f])/ISAPT.values[i:f,0]),columns=[['DLPNO']])\n",
    "pDFT_SVP=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVP.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_SVP']])\n",
    "pDFT_TZVP=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVP.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_TZVP']])\n",
    "pDFT_SVPD=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2SVPD.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_SVPD']])\n",
    "pDFT_TZVPD=df(np.abs(100*(ISAPT.values[i:f,0]-NN.wB97XD3_def2TZVPD.values[i:f])/ISAPT.values[i:f,0]),columns=[['wB97XD3_TZVPD']])\n",
    "\n",
    "benchmark26p=pd.concat([pDLPNO,pDFT_SVP,pDFT_TZVP,pDFT_SVPD,pDFT_TZVPD], axis=1)\n",
    "\n",
    "mad=np.mean(benchmark26.abs())\n",
    "std=benchmark26.std(ddof=0)\n",
    "mpad=np.mean(benchmark26p)\n",
    "#metrics=df()\n",
    "metrics['MAD-all']=mad\n",
    "metrics['STD-all']=std\n",
    "metrics['MPAD-all']=mpad\n",
    "print(metrics)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 16.08 ( 10.33) &  14.87 &  9.52 (  6.35) &   7.08 &  9.74 (  5.48) &   6.89 & 11.62 (  8.11) &   9.41 &\n",
      "\n",
      " 24.36 ( 11.09) &  21.96 & 20.04 (  9.28) &  14.89 & 21.68 (  8.99) &  15.44 & 21.94 (  9.93) &  17.26 &\n",
      "\n",
      " 18.16 ( 10.50) &  16.63 & 12.38 (  7.34) &   9.16 & 13.36 (  6.77) &   9.46 & 14.50 (  8.63) &  11.56 &\n",
      "\n",
      " 15.97 ( 10.73) &  14.72 &  9.78 (  6.81) &   7.25 & 10.19 (  5.82) &   7.19 & 11.82 (  8.40) &   9.53 &\n",
      "\n",
      " 16.07 ( 10.58) &  14.91 &  9.75 (  6.50) &   7.27 & 10.14 (  5.39) &   7.19 & 11.83 (  8.18) &   9.59 &\n",
      "\n"
     ]
    }
   ],
   "source": [
    "Nmethods=len(metrics)\n",
    "for imethod in range(Nmethods):\n",
    "    a=metrics[imethod:imethod+1]\n",
    "    b=np.array(a)\n",
    "    Ncols=b.shape[1]\n",
    "    for icols in range(0,Ncols,3):\n",
    "        print('{:6.2f} ({:6.2f}) & {:6.2f} &'.format(b[0,icols],b[0,icols+1],b[0,icols+2]),end=\"\")\n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-CH3\n",
    "```\n",
    " 16.08 ( 10.33) &  14.87 &  9.52 (  6.35) &   7.08 &  9.74 (  5.48) &   6.89 & 11.62 (  8.11) &   9.41 &\n",
    " 24.36 ( 11.09) &  21.96 & 20.04 (  9.28) &  14.89 & 21.68 (  8.99) &  15.44 & 21.94 (  9.93) &  17.26 &\n",
    " 18.16 ( 10.50) &  16.63 & 12.38 (  7.34) &   9.16 & 13.36 (  6.77) &   9.46 & 14.50 (  8.63) &  11.56 &\n",
    " 15.97 ( 10.73) &  14.72 &  9.78 (  6.81) &   7.25 & 10.19 (  5.82) &   7.19 & 11.82 (  8.40) &   9.53 &\n",
    " 16.07 ( 10.58) &  14.91 &  9.75 (  6.50) &   7.27 & 10.14 (  5.39) &   7.19 & 11.83 (  8.18) &   9.59 &\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-CH2CH3\n",
    "```\n",
    " 12.34 ( 10.45) &  11.89 &  4.88 (  4.77) &   3.84 &  4.31 (  2.72) &   3.17 &  6.98 (  7.28) &   6.09 &\n",
    " 17.52 (  9.78) &  16.55 & 12.82 (  5.47) &   9.80 & 14.53 (  3.93) &  10.52 & 14.86 (  6.98) &  12.13 &\n",
    " 13.04 ( 10.12) &  12.58 &  6.83 (  4.92) &   5.29 &  7.33 (  3.02) &   5.32 &  8.92 (  7.05) &   7.54 &\n",
    " 12.43 ( 10.55) &  11.89 &  5.41 (  5.19) &   4.21 &  4.83 (  3.25) &   3.52 &  7.37 (  7.40) &   6.33 &\n",
    " 12.47 ( 10.53) &  12.00 &  5.34 (  5.11) &   4.19 &  4.91 (  2.83) &   3.60 &  7.38 (  7.34) &   6.39 &\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-CH2CH2CH3\n",
    "```\n",
    " 10.93 ( 10.87) &  10.51 &  4.21 (  4.76) &   3.32 &  2.73 (  2.49) &   2.07 &  5.77 (  7.25) &   5.10 &\n",
    " 13.09 ( 10.02) &  12.79 &  9.48 (  5.08) &   7.30 & 11.41 (  2.91) &   8.29 & 11.26 (  6.70) &   9.33 &\n",
    " 11.13 ( 10.50) &  10.80 &  4.97 (  4.80) &   3.88 &  4.83 (  2.44) &   3.53 &  6.81 (  6.94) &   5.89 &\n",
    " 10.71 ( 10.86) &  10.32 &  4.57 (  5.07) &   3.60 &  3.25 (  2.57) &   2.40 &  6.00 (  7.24) &   5.25 &\n",
    " 10.87 ( 10.88) &  10.47 &  4.55 (  5.05) &   3.59 &  2.97 (  2.52) &   2.22 &  5.95 (  7.26) &   5.23 &\n",
    "```\n",
    "\n",
    "junccpvdz\n",
    "\n",
    "-CH3\n",
    "```\n",
    " 19.24 ( 10.35) &  18.39 & 12.64 (  6.40) &   9.75 & 13.73 (  5.32) &  10.15 & 15.05 (  8.06) &  12.55 &\n",
    " 27.92 ( 11.14) &  26.06 & 23.71 (  9.35) &  18.22 & 25.67 (  8.86) &  18.94 & 25.68 (  9.93) &  20.88 &\n",
    " 21.72 ( 10.53) &  20.54 & 16.04 (  7.42) &  12.33 & 17.35 (  6.62) &  12.79 & 18.24 (  8.60) &  15.01 &\n",
    " 19.04 ( 10.76) &  18.16 & 12.89 (  6.89) &   9.93 & 14.17 (  5.66) &  10.45 & 15.23 (  8.36) &  12.64 &\n",
    " 19.28 ( 10.59) &  18.48 & 12.97 (  6.59) &  10.04 & 14.13 (  5.22) &  10.45 & 15.31 (  8.14) &  12.78 &\n",
    "\n",
    "```\n",
    "-CH2CH3\n",
    "```\n",
    " 13.98 ( 10.41) &  14.10 &  7.25 (  4.83) &   5.86 &  8.29 (  2.60) &   6.32 &  9.68 (  7.21) &   8.55 &\n",
    " 21.08 (  9.77) &  20.48 & 16.49 (  5.58) &  13.00 & 18.52 (  3.82) &  13.89 & 18.60 (  6.98) &  15.61 &\n",
    " 15.83 ( 10.10) &  15.73 &  9.93 (  5.01) &   7.90 & 11.32 (  2.88) &   8.53 & 12.23 (  7.00) &  10.53 &\n",
    " 14.06 ( 10.51) &  14.09 &  7.83 (  5.27) &   6.29 &  8.81 (  3.11) &   6.67 & 10.09 (  7.35) &   8.82 &\n",
    " 14.17 ( 10.49) &  14.27 &  7.67 (  5.19) &   6.19 &  8.90 (  2.70) &   6.76 & 10.10 (  7.29) &   8.88 &\n",
    "```\n",
    "-CH2CH2CH3\n",
    "```\n",
    " 12.30 ( 10.82) &  12.49 &  5.72 (  4.80) &   4.68 &  6.22 (  2.40) &   4.79 &  7.92 (  7.19) &   7.12 &\n",
    " 16.55 ( 10.00) &  16.52 & 13.15 (  5.22) &  10.43 & 15.39 (  2.83) &  11.59 & 14.97 (  6.71) &  12.70 &\n",
    " 12.79 ( 10.45) &  12.98 &  7.33 (  4.89) &   5.90 &  8.81 (  2.33) &   6.68 &  9.52 (  6.90) &   8.35 &\n",
    " 12.33 ( 10.81) &  12.46 &  6.23 (  5.14) &   5.05 &  7.00 (  2.45) &   5.34 &  8.37 (  7.19) &   7.43 &\n",
    " 12.32 ( 10.82) &  12.50 &  6.12 (  5.12) &   4.98 &  6.77 (  2.42) &   5.19 &  8.25 (  7.21) &   7.37 &\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}