Clathrate_MBE_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clathrate cluster analysis\n",
    "\n",
    "This notebook contains analysis of the data used to produce results for methane in water clathrate paper.\n",
    "\n",
    "The ordering is as follows:\n",
    "* useful stuff is defined\n",
    "* data is loaded\n",
    "* analysis in order:\n",
    " * a\n",
    " * b\n",
    "\n",
    "Data are usually stored in arrays with indices used to describe basis set, functional, fragment index, use of singles corrections,  ...\n",
    "The indices of basis set or functional are stored in a python dictionary.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 315,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import stuff\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "#our little library with some useful functions (basis-set extrapolation, getting MBE contributions, ...)\n",
    "import mbe_lib as mbe\n",
    "\n",
    "#########################\n",
    "#Basis set dictionaries #\n",
    "#########################\n",
    "\n",
    "#standard basis sets\n",
    "basid={}\n",
    "basid['TZ']=0\n",
    "basid['QZ']=1\n",
    "basid['CBS']=2\n",
    "\n",
    "#QZ and 5Z basis sets\n",
    "bas5id={}\n",
    "bas5id['QZ']=0\n",
    "bas5id['5Z']=1\n",
    "bas5id['CBS']=2\n",
    "\n",
    "#DZ and TZ basis sets\n",
    "basDTid={}\n",
    "basDTid['DZ']=0\n",
    "basDTid['TZ']=1\n",
    "basDTid['CBS']=2\n",
    "\n",
    "#DZ, TZ, QZ basis sets\n",
    "bas24id={}\n",
    "bas24id['DZ']=0\n",
    "bas24id['TZ']=1\n",
    "bas24id['QZ']=2\n",
    "bas24id['CBS']=3\n",
    "\n",
    "#basis names for Marcin's code\n",
    "MM_basis={}\n",
    "MM_basis['small-basis']=0\n",
    "MM_basis['large-basis']=1\n",
    "MM_basis['CBS']=2\n",
    "\n",
    "###################\n",
    "# Unit conversion #\n",
    "###################\n",
    "#Hartree to kcal/mol (\"Marcin's units\")\n",
    "ha2kc=627.50960803059\n",
    "#Hartree to kJ/mol (\"Sirous' units\")\n",
    "ha2kj=2625.5002\n",
    "#Hartree to eV (\"Jiri's units\")\n",
    "ha2ev=27.211399\n",
    "\n",
    "#########################\n",
    "# DFT input eigenstates #\n",
    "#########################\n",
    "#for usual approximations\n",
    "methods={}\n",
    "methods['PBE']=0\n",
    "methods['PBE0']=1\n",
    "methods['SCAN']=2\n",
    "methods['SCAN0']=3\n",
    "methods['HF']=4\n",
    "\n",
    "#and for OEP\n",
    "oep_m={}\n",
    "oep_m['oepx']=0\n",
    "oep_m['oep2sc']=1\n",
    "oep_m['hf']=2\n",
    "\n",
    "##################\n",
    "# Arrays for MBE #\n",
    "##################\n",
    "#RPA: indices are methods; RSE (no and yes); MBE order 0th, 2-body, 3-body (1hb and 0hb), 4-body (2hb, 1hb, 0hb)\n",
    "RPA_mbe=np.zeros((5,2,7))\n",
    "#Moller-Plesset PT, MP2, MP3, MP4; MBE order\n",
    "MP_mbe=np.zeros((3,7))\n",
    "#CC with AVQZ and AV5Z: method HF+CABS, CCSD-F12, (T), HF+CABS+CCSD-F12(T); MBE order \n",
    "CC_mbeQ5=np.zeros((4,7))\n",
    "#CC with canonical CCSD(T)\n",
    "CCcan_mbeQ5=np.zeros((4,7))\n",
    "\n",
    "#interaction en ergies from supermolecular calculations in kcal/mol\n",
    "RPA_sup=np.array([[-3.73, -4.99],[-3.30, -4.03],  [-3.38, -5.11], [-3.17, -4.18],   [0.0 ,0.0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 316,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################\n",
    "# Trimers distance arrays set-up #\n",
    "##################################\n",
    "\n",
    "#The three-body and four-body contributions are divided according to the distance and number of hydrogen bonds\n",
    "#here we calculate the distances \n",
    "#There are four groups when ordering is done according to distance\n",
    "#There are two groups when orering is done according to number of hydrogen bonds\n",
    "\n",
    "#array storing distance between pair of water molecules\n",
    "dist=np.zeros((20,20))\n",
    "#array storing distance between water molecules in (indexed, ordered) trimers\n",
    "dist3_list=np.zeros((190))\n",
    "#number of fragments assigned in groups according to distance\n",
    "NN_num3=np.zeros((4))\n",
    "#distance group index for a given trimer\n",
    "trim_nn_case=np.zeros((190),dtype=np.int16)\n",
    "# arrays for transforming from pair of molecules (n,m) list into ordered trimers (1..190) list \n",
    "trim_idx=np.zeros((20,20),dtype=np.int32)\n",
    "trim_to2D=np.zeros((190,2),dtype=np.int32)\n",
    "\n",
    "#load the distances of water molecules, they are stored in the file\n",
    "idt=0\n",
    "with open('data/dist_trimer.dat','r') as f:\n",
    "    for i in range(19): #loop over first monomers 0 to 18\n",
    "        for j in range(i+1,20):  #loop over monomers i+1 to 19\n",
    "            #read the distance between molecules i and j\n",
    "            dist[i,j]=float(f.readline())\n",
    "            #put the conversion indices into right places\n",
    "            trim_idx[i,j]=idt\n",
    "            trim_to2D[idt,0]=i\n",
    "            trim_to2D[idt,1]=j\n",
    "            #assign the trimer to a distance group\n",
    "            if dist[i,j]<3.2:\n",
    "                trim_nn_case[idt]=0\n",
    "            elif dist[i,j]<5.5:\n",
    "                trim_nn_case[idt]=1\n",
    "            elif dist[i,j]<6.75:\n",
    "                trim_nn_case[idt]=2\n",
    "            else:\n",
    "                trim_nn_case[idt]=3\n",
    "            #increase the number of trimers in the assigned group\n",
    "            NN_num3[trim_nn_case[idt]]+=1\n",
    "            #go to next trimer\n",
    "            idt+=1\n",
    "            #possibly check that things are working fine\n",
    "            #if i==0:\n",
    "            #    print(i,j,dist[i,j])\n",
    "\n",
    "\n",
    "#############################\n",
    "# Tetramers distance set-up #\n",
    "#############################\n",
    "\n",
    "#assign tetramers into groups according to the number of hydrogen bonds in the fragment (0, 1, or 2)\n",
    "#distance in tetramers, stored as i,j,k and ordered tetramer index\n",
    "dist_tetr=np.zeros((20,20,20))\n",
    "dist_list=np.zeros((1140))\n",
    "#number of nearest neighbours in tetramer\n",
    "NN_num4=np.zeros((3))\n",
    "#array to store the group index (number of h-bonds)\n",
    "tetr_nn_num=np.zeros((1140),dtype=np.int16)\n",
    "#back and forth indices between [0,i,j,k] and ordered tetramer index\n",
    "tetr_idx=np.zeros((20,20,20),dtype=np.int32)\n",
    "tetr_to3D=np.zeros((1141,3),dtype=np.int32)\n",
    "idt=0\n",
    "for i in range(18): #loop over monomers 0 to 17\n",
    "    for j in range(i+1,19): #loop over i+1 to 18\n",
    "        for k in range(j+1,20): #loop over j+1 to 19\n",
    "            #calculate the total distance in the tetramer  (? maybe no reason for this...)\n",
    "            dist_tetr[i,j,k]=dist[i,j]+dist[i,k]+dist[j,k]\n",
    "            dist_list[idt]=dist_tetr[i,j,k]\n",
    "            #count the number of water nearest neighbours in the tetramer\n",
    "            #check the first dimer, then the second, then the last\n",
    "            if dist[i,j] <3.2:\n",
    "                tetr_nn_num[idt]+=1\n",
    "            if dist[i,k] <3.2:\n",
    "                tetr_nn_num[idt]+=1\n",
    "            if dist[j,k] <3.2:\n",
    "                tetr_nn_num[idt]+=1\n",
    "            #add one to the array counting the number of NNs\n",
    "            NN_num4[tetr_nn_num[idt]]+=1\n",
    "            #store the transition matrix between i,j,k and ordered index (not sure if needed)\n",
    "            tetr_idx[i,j,k]=idt\n",
    "            tetr_to3D[idt,0]=i\n",
    "            tetr_to3D[idt,1]=j\n",
    "            tetr_to3D[idt,2]=k\n",
    "            #go to next one\n",
    "            idt+=1\n",
    "            \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 317,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Few basic functions to process the data\n",
    "\n",
    "#calculate basic statistics (total, ME, MAE, RMSD from a reference set)\n",
    "#input: test   = data for which statistics is calculated\n",
    "#       ref    = reference data\n",
    "def calc_stat(test,ref):\n",
    "    if len(test) != len(ref):\n",
    "        print(len(test), len(ref), ' are different')\n",
    "    ME=0.\n",
    "    MAE=0.\n",
    "    RMSD=0.\n",
    "    #print ('Total ', np.sum(test), np.sum(ref))\n",
    "    AVG1=np.sum(test)/len(ref)\n",
    "    AVG2=np.sum(ref)/len(ref)\n",
    "    for i in range(len(test)):\n",
    "        ME+=test[i]-ref[i]\n",
    "        MAE+=abs(test[i]-ref[i])\n",
    "        RMSD+=(test[i]-ref[i])**2\n",
    "    ME=ME/len(ref)\n",
    "    MAE=MAE/len(ref)\n",
    "    RMSD=np.sqrt(RMSD/len(test))\n",
    "    \n",
    "    #print the statistics\n",
    "    print('Total  ME  MAE  RMSD')\n",
    "    print(np.sum(test), ME , MAE, RMSD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read data\n",
    "\n",
    "The total energies are read from output files or from data files into arrays with _orig_ in their name, the many-body interaction energy is then calculated using mbe library and stored in array with _val_ in its name."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-body energies\n",
    "\n",
    "Reference data\n",
    "\n",
    "* CC with accu=1e-16, orbital=1d-10, energy=1d-10, data/dimer/CC_orb10, AVTZ, AVQZ, AV5Z basis\n",
    " * stored as Molpro output files\n",
    " * read into HF_orig_o10, CCSD_orig_o10, CCSDT_orig_o10 arrays for AVTZ and AVQZ basis sets\n",
    " * read into HF_orig_o10Q5, CCSD_orig_o10Q5, CCSDT_orig_o10Q5 arrays for AVQZ and AV5Z basis sets\n",
    " \n",
    "* CC-F12 with accu=1e-16, orbital=1d-10, energy=1d-10, data/dimer/CCF12_orb10, AVTZ, AVQZ, AV5Z basis \n",
    " * stored as Molpro output files\n",
    " * read into CABS_orig_Fo10, CCSD_orig_Fo10, CCSDF12b_orig_Fo10, Ts_orig_Fo10, Tu_orig_Fo10 arrays for AVTZ and AVQZ basis sets\n",
    " * read into CABS_orig_Fo10Q5, CCSD_orig_Fo10Q5, CCSDF12b_orig_Fo10Q5, Ts_orig_Fo10Q5, Tu_orig_Fo10Q5 arrays for AVQZ and AV5Z basis sets\n",
    " \n",
    "CCSD-F12 using VTZ-F12 basis set as in Deible et al.\n",
    "* CC-F12 with orbital=1d-8, energy=1d-8, twoint=1.d-19, VTZ-F12 basis\n",
    " * stored as Molpro output files\n",
    " * read into CABS_orig_Fo8, CCSD_orig_Fo8, CCSDF12b_orig_Fo8, Ts_orig_Fo8, Tu_orig_Fo8 arrays for VTZ-F12 basis set\n",
    "\n",
    "MPn\n",
    "* with orbital=1d-8,energy=1d-8,twoint=1.d-19, data/dimer/MPn, AVTZ and AVQZ basis\n",
    " * stored as Molpro output files\n",
    "\n",
    "Notes:\n",
    "* Different twoint screening makes no significant difference to dimers(10^-6 kcal/mol for CC)\n",
    "* Orbital and energy convergence set to 1d-8 or to 1d-10 makes no significant difference to dimers (10^-5 kcal/mol for CC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [],
   "source": [
    "#####################\n",
    "# Canonical CCSD(T) #\n",
    "#####################\n",
    "\n",
    "# AVTZ and AVQZ basis set\n",
    "\n",
    "#HF_orig, indices: basis-set (AVTZ and AVQZ); dimer index; dimer and monomer energies\n",
    "#HF_val, indices: basis-set (AVTZ, AVQZ, CBS); dimer index\n",
    "HF2_orig_TQ=np.zeros((2,20,3),dtype=np.float64)\n",
    "HF2_val_TQ=np.zeros((3,20),dtype=np.float64)\n",
    "#CCSD correlation energy\n",
    "CCSD2_orig_TQ=np.zeros((2,20,3),dtype=np.float64)\n",
    "CCSD2_val_TQ=np.zeros((3,20),dtype=np.float64)\n",
    "#CCSD(T) correlation energy\n",
    "CCSDT2_orig_TQ=np.zeros((2,20,3),dtype=np.float64)\n",
    "CCSDT2_val_TQ=np.zeros((3,20),dtype=np.float64)\n",
    "#(T) correlation energy\n",
    "T2_val_TQ=np.zeros((3,20),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(20):\n",
    "        file='data/dimer/CC_orb10/'+basis+'/'+str(\"{:02d}\".format(i+1))+'_Dime_CC_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                #process the file\n",
    "                filelist=f.readlines()\n",
    "                # mbe.read functions process the file and return a tuple (or whatever that is) of lists\n",
    "                data=mbe.read_CC_F12_Molpro(filelist)\n",
    "                HF2_orig_TQ[basid[basis],i,:]=data[0]\n",
    "                HF2_val_TQ[basid[basis],i]=mbe.get_2body(data[0])*ha2kc\n",
    "                CCSD2_orig_TQ[basid[basis],i,:]=data[4]\n",
    "                CCSD2_val_TQ[basid[basis],i]=mbe.get_2body(data[4])*ha2kc\n",
    "                #read_MP4_Molpro actually looks only for \"Total correlation energy\"\n",
    "                data=mbe.read_MP4_Molpro(filelist)\n",
    "                CCSDT2_orig_TQ[basid[basis],i,:]=data\n",
    "                CCSDT2_val_TQ[basid[basis],i]=mbe.get_2body(data)*ha2kc \n",
    "                #store the data in (T)\n",
    "                T2_val_TQ[basid[basis],i]=CCSDT2_val_TQ[basid[basis],i]-CCSD2_val_TQ[basid[basis],i]\n",
    "                \n",
    "                \n",
    "for i in range(20):\n",
    "    CCSD2_val_TQ[basid['CBS'],i]=mbe.basis_extrap(CCSD2_val_TQ[basid['TZ'],i],CCSD2_val_TQ[basid['QZ'],i],3,4)\n",
    "    CCSDT2_val_TQ[basid['CBS'],i]=mbe.basis_extrap(CCSDT2_val_TQ[basid['TZ'],i],CCSDT2_val_TQ[basid['QZ'],i],3,4)\n",
    "    T2_val_TQ[basid['CBS'],i]=mbe.basis_extrap(T2_val_TQ[basid['TZ'],i],T2_val_TQ[basid['QZ'],i],3,4)\n",
    "    \n",
    "# AVQZ and AV5Z basis set\n",
    "\n",
    "#Data with QZ and 5Z basis sets and canonical CCSD(T)\n",
    "HF2_orig_Q5=np.zeros((2,20,3),dtype=np.float64)\n",
    "HF2_val_Q5=np.zeros((3,20),dtype=np.float64)\n",
    "#CCSD correlation energy\n",
    "CCSD2_orig_Q5=np.zeros((2,20,3),dtype=np.float64)\n",
    "CCSD2_val_Q5=np.zeros((3,20),dtype=np.float64)\n",
    "#CCSD(T) correlation energy\n",
    "CCSDT2_orig_Q5=np.zeros((2,20,3),dtype=np.float64)\n",
    "CCSDT2_val_Q5=np.zeros((3,20),dtype=np.float64)\n",
    "#T correlation\n",
    "T2_val_Q5=np.zeros((3,20),dtype=np.float64)\n",
    "\n",
    "\n",
    "for basis in ['QZ','5Z']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(20):\n",
    "        file='data/dimer/CC_orb10/'+basis+'/'+str(\"{:02d}\".format(i+1))+'_Dime_CC_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_CC_F12_Molpro(filelist)\n",
    "                HF2_orig_Q5[bas5id[basis],i,:]=data[0]\n",
    "                HF2_val_Q5[bas5id[basis],i]=mbe.get_2body(data[0])*ha2kc\n",
    "                CCSD2_orig_Q5[bas5id[basis],i,:]=data[4]\n",
    "                CCSD2_val_Q5[bas5id[basis],i]=mbe.get_2body(data[4])*ha2kc\n",
    "                data=mbe.read_MP4_Molpro(filelist)\n",
    "                CCSDT2_orig_Q5[bas5id[basis],i,:]=data\n",
    "                CCSDT2_val_Q5[bas5id[basis],i]=mbe.get_2body(data)*ha2kc \n",
    "                T2_val_Q5[bas5id[basis],i]=CCSDT2_val_Q5[bas5id[basis],i]-CCSD2_val_Q5[bas5id[basis],i]\n",
    "\n",
    "#extrapolate to CBS                \n",
    "for i in range(20):\n",
    "    CCSD2_val_Q5[bas5id['CBS'],i]=mbe.basis_extrap(CCSD2_val_Q5[bas5id['QZ'],i],CCSD2_val_Q5[bas5id['5Z'],i],4,5)\n",
    "    CCSDT2_val_Q5[bas5id['CBS'],i]=mbe.basis_extrap(CCSDT2_val_Q5[bas5id['QZ'],i],CCSDT2_val_Q5[bas5id['5Z'],i],4,5)\n",
    "    T2_val_Q5[bas5id['CBS'],i]=mbe.basis_extrap(T2_val_Q5[bas5id['QZ'],i],T2_val_Q5[bas5id['5Z'],i],4,5)\n",
    "    \n",
    "# store final CCSD and CCSD(T) values in additional arrays\n",
    "CCSD2_tot=np.zeros((20))\n",
    "CCSDT2_tot=np.zeros((20))\n",
    "for i in range(20):\n",
    "    CCSD2_tot[i]=CCSD2_val_Q5[bas5id['CBS'],i]+HF2_val_Q5[bas5id['5Z'],i]\n",
    "    CCSDT2_tot[i]=CCSDT2_val_Q5[bas5id['CBS'],i]+HF2_val_Q5[bas5id['5Z'],i]\n",
    "\n",
    "#check that we get -1.640 as in Table 1 for AVQZ->AV5Z extrapolated (T) contribution\n",
    "#print(np.sum(CCSDT_tot)-np.sum(CCSD_tot))\n",
    "\n",
    "#we now put the MBE two-body component into the CCcan_mbeQ5 array\n",
    "#first index is HF, CC, (T), total; second is order of MBE\n",
    "CCcan_mbeQ5[0,1]=np.sum(HF2_val_Q5[bas5id['5Z'],:])\n",
    "CCcan_mbeQ5[1,1]=np.sum(CCSD2_val_Q5[bas5id['CBS'],:])\n",
    "CCcan_mbeQ5[2,1]=np.sum(T2_val_Q5[bas5id['CBS'],:])\n",
    "CCcan_mbeQ5[3,1]=np.sum(CCSDT2_tot[:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 319,
   "metadata": {},
   "outputs": [],
   "source": [
    "#################\n",
    "# CCSD(T)-F12b reference data\n",
    "#################\n",
    "\n",
    "# AVTZ and AVQZ results first\n",
    "#HF+CABS data\n",
    "CABS2_orig_FTQ=np.zeros((2,20,3),dtype=np.float64)\n",
    "CABS2_val_FTQ=np.zeros((3,20),dtype=np.float64)\n",
    "#CCSD \n",
    "CCSD2_orig_FTQ=np.zeros((2,20,3),dtype=np.float64)\n",
    "CCSD2_val_FTQ=np.zeros((3,20),dtype=np.float64)\n",
    "#CCSD-F12b\n",
    "CCSDF12b2_orig_FTQ=np.zeros((2,20,3),dtype=np.float64)\n",
    "CCSDF12b2_val_FTQ=np.zeros((3,20),dtype=np.float64)\n",
    "#scaled Triples\n",
    "Ts2_orig_FTQ=np.zeros((2,20,3),dtype=np.float64)\n",
    "Ts2_val_FTQ=np.zeros((3,20),dtype=np.float64)\n",
    "#unscaled Triples\n",
    "Tu2_orig_FTQ=np.zeros((2,20,3),dtype=np.float64)\n",
    "Tu2_val_FTQ=np.zeros((3,20),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(20):\n",
    "        file='data/dimer/CCF12_orb10/'+basis+'/'+str(\"{:02d}\".format(i+1))+'_Dimer_F12_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_CC_F12_Molpro(filelist)\n",
    "                # return HF, CABS, MP2, MF12, CC, CCF12b, Ts, Tu, F12a, F12b\n",
    "                #sort out the data into what we need\n",
    "                CABS2_orig_FTQ[basid[basis],i,:]=data[1]\n",
    "                CCSD2_orig_FTQ[basid[basis],i,:]=data[4] \n",
    "                CCSDF12b2_orig_FTQ[basid[basis],i,:]=data[5]\n",
    "                Ts2_orig_FTQ[basid[basis],i,:]=data[6]\n",
    "                Tu2_orig_FTQ[basid[basis],i,:]=data[7]\n",
    "                #get the two-body interaction\n",
    "                CABS2_val_FTQ[basid[basis],i]=mbe.get_2body(data[1])*ha2kc\n",
    "                CCSD2_val_FTQ[basid[basis],i]=mbe.get_2body(data[4])*ha2kc\n",
    "                CCSDF12b2_val_FTQ[basid[basis],i]=mbe.get_2body(data[5])*ha2kc\n",
    "                Ts2_val_FTQ[basid[basis],i]=mbe.get_2body(data[6])*ha2kc\n",
    "                Tu2_val_FTQ[basid[basis],i]=mbe.get_2body(data[7])*ha2kc\n",
    "                \n",
    "#do basis-set extrapolation\n",
    "for i in range(20):\n",
    "    CCSD2_val_FTQ[basid['CBS'],i]=mbe.basis_extrap(CCSD2_val_FTQ[basid['TZ'],i],CCSD2_val_FTQ[basid['QZ'],i],3,4)\n",
    "    #Ts_val_FTQ[basid['CBS'],i]=mbe.basis_extrap(Ts_val_FTQ[basid['TZ'],i],Ts_val_FTQ[basid['QZ'],i],3,4)\n",
    "    Tu2_val_FTQ[basid['CBS'],i]=mbe.basis_extrap(Tu2_val_FTQ[basid['TZ'],i],Tu2_val_FTQ[basid['QZ'],i],3,4)\n",
    "    \n",
    "#store the details\n",
    "#CCSDF12_tot=np.zeros((20))\n",
    "#CCSDTF12_tot=np.zeros((20))\n",
    "#for i in range(20):\n",
    "#    CCSDF12_tot[i]=CCSD_val_Fo10[basid['CBS'],i]+CABS_val_Fo10[basid['QZ'],i]\n",
    "#    CCSDTF12_tot[i]=CCSDF12_tot[i]+Ts_val_Fo10[basid['CBS'],i]\n",
    "    \n",
    "#########################################################\n",
    "# The reference CCSD(T)-F12 data in QZ and 5Z basis-set #\n",
    "#########################################################\n",
    "#HF+CABS data\n",
    "CABS2_orig_FQ5=np.zeros((2,20,3),dtype=np.float64)\n",
    "CABS2_val_FQ5=np.zeros((3,20),dtype=np.float64)\n",
    "#CCSD \n",
    "CCSD2_orig_FQ5=np.zeros((2,20,3),dtype=np.float64)\n",
    "CCSD2_val_FQ5=np.zeros((3,20),dtype=np.float64)\n",
    "#CCSD-F12b\n",
    "CCSDF12b2_orig_FQ5=np.zeros((2,20,3),dtype=np.float64)\n",
    "CCSDF12b2_val_FQ5=np.zeros((3,20),dtype=np.float64)\n",
    "#scaled Triples\n",
    "Ts2_orig_FQ5=np.zeros((2,20,3),dtype=np.float64)\n",
    "Ts2_val_FQ5=np.zeros((3,20),dtype=np.float64)\n",
    "#unscaled Triples\n",
    "Tu2_orig_FQ5=np.zeros((2,20,3),dtype=np.float64)\n",
    "Tu2_val_FQ5=np.zeros((3,20),dtype=np.float64)\n",
    "\n",
    "for basis in ['QZ','5Z']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(20):\n",
    "        file='data/dimer/CCF12_orb10/'+basis+'/'+str(\"{:02d}\".format(i+1))+'_Dimer_F12_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_CC_F12_Molpro(filelist)\n",
    "                # return HF, CABS, MP2, MF12, CC, CCF12b, Ts, Tu, F12a, F12b\n",
    "                #sort out the data into what we need\n",
    "                CABS2_orig_FQ5[bas5id[basis],i,:]=data[1]\n",
    "                CCSD2_orig_FQ5[bas5id[basis],i,:]=data[4] \n",
    "                CCSDF12b2_orig_FQ5[bas5id[basis],i,:]=data[5]\n",
    "                Ts2_orig_FQ5[bas5id[basis],i,:]=data[6]\n",
    "                Tu2_orig_FQ5[bas5id[basis],i,:]=data[7]\n",
    "                ##print (data)\n",
    "                CABS2_val_FQ5[bas5id[basis],i]=mbe.get_2body(data[1])*ha2kc\n",
    "                CCSD2_val_FQ5[bas5id[basis],i]=mbe.get_2body(data[4])*ha2kc\n",
    "                CCSDF12b2_val_FQ5[bas5id[basis],i]=mbe.get_2body(data[5])*ha2kc\n",
    "                Ts2_val_FQ5[bas5id[basis],i]=mbe.get_2body(data[6])*ha2kc\n",
    "                Tu2_val_FQ5[bas5id[basis],i]=mbe.get_2body(data[7])*ha2kc\n",
    "                \n",
    "#do basis-set extrapolation for CCSD printed out by the code and for unscaled (T) contribution\n",
    "for i in range(20):\n",
    "    CCSD2_val_FQ5[bas5id['CBS'],i]=mbe.basis_extrap(CCSD2_val_FQ5[bas5id['QZ'],i],CCSD2_val_FQ5[bas5id['5Z'],i],4,5)\n",
    "    #Ts2_val_FQ5[bas5id['CBS'],i]=mbe.basis_extrap(Ts2_val_FQ5[bas5id['QZ'],i],Ts2_val_FQ5[bas5id['5Z'],i],4,5)\n",
    "    Tu2_val_FQ5[bas5id['CBS'],i]=mbe.basis_extrap(Tu2_val_FQ5[bas5id['QZ'],i],Tu2_val_FQ5[bas5id['5Z'],i],4,5)\n",
    "\n",
    "#store the totals\n",
    "CCSDF12b2_tot=np.zeros((20))\n",
    "CCSDTF12b2_tot=np.zeros((20))\n",
    "for i in range(20):\n",
    "    #this is HF+CABS+CCSD-F12b\n",
    "    CCSDF12b2_tot[i]=CCSDF12b2_val_FQ5[bas5id['5Z'],i]+CABS2_val_FQ5[bas5id['5Z'],i]\n",
    "    #this is HF+CABS+CCSD-F12b+(T)_scaled\n",
    "    CCSDTF12b2_tot[i]=CCSDF12b2_tot[i]+Ts2_val_FQ5[bas5id['5Z'],i]\n",
    "\n",
    "#store the total contributions to MBE, \n",
    "#first index is HF+CABS, CC-F12, (T), total; second is order of MBE\n",
    "CC_mbeQ5[0,1]=np.sum(CABS2_val_FQ5[bas5id['5Z'],:])\n",
    "CC_mbeQ5[1,1]=np.sum(CCSDF12b2_val_FQ5[bas5id['5Z'],:])\n",
    "CC_mbeQ5[2,1]=np.sum(Ts2_val_FQ5[bas5id['5Z'],:])\n",
    "CC_mbeQ5[3,1]=CC_mbeQ5[0,1]+CC_mbeQ5[1,1]+CC_mbeQ5[2,1]\n",
    "  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 320,
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################################\n",
    "# CCSD(T)-F12 with the VTZ-F12 basis set #\n",
    "##########################################\n",
    "\n",
    "#HF+CABS data\n",
    "CABS2_orig_Fo8=np.zeros((20,3),dtype=np.float64)\n",
    "CABS2_val_Fo8=np.zeros((20),dtype=np.float64)\n",
    "#CCSD \n",
    "CCSD2_orig_Fo8=np.zeros((20,3),dtype=np.float64)\n",
    "CCSD2_val_Fo8=np.zeros((20),dtype=np.float64)\n",
    "#CCSD-F12b\n",
    "CCSDF12b2_orig_Fo8=np.zeros((20,3),dtype=np.float64)\n",
    "CCSDF12b2_val_Fo8=np.zeros((20),dtype=np.float64)\n",
    "#scaled Triples\n",
    "Ts2_orig_Fo8=np.zeros((20,3),dtype=np.float64)\n",
    "Ts2_val_Fo8=np.zeros((20),dtype=np.float64)\n",
    "#unscaled Triples\n",
    "Tu2_orig_Fo8=np.zeros((20,3),dtype=np.float64)\n",
    "Tu2_val_Fo8=np.zeros((20),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(20):\n",
    "        file='data/dimer/CCF12_orb8/V'+basis+'/'+str(\"{:02d}\".format(i+1))+'_Dimer_F12_'+basis+'_F12.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_CC_F12_Molpro(filelist)\n",
    "                # return HF, CABS, MP2, MF12, CC, CCF12b, Ts, Tu, F12a, F12b\n",
    "                #sort out the data into what we need\n",
    "                CABS2_orig_Fo8[i,:]=data[1]\n",
    "                CCSD2_orig_Fo8[i,:]=data[4] \n",
    "                CCSDF12b2_orig_Fo8[i,:]=data[5]\n",
    "                Ts2_orig_Fo8[i,:]=data[6]\n",
    "                Tu2_orig_Fo8[i,:]=data[7]\n",
    "                #calculate data\n",
    "                CABS2_val_Fo8[i]=mbe.get_2body(data[1])*ha2kc\n",
    "                CCSD2_val_Fo8[i]=mbe.get_2body(data[4])*ha2kc\n",
    "                CCSDF12b2_val_Fo8[i]=mbe.get_2body(data[5])*ha2kc\n",
    "                Ts2_val_Fo8[i]=mbe.get_2body(data[6])*ha2kc\n",
    "                Tu2_val_Fo8[i]=mbe.get_2body(data[7])*ha2kc\n",
    "\n",
    "#print(np.sum(CABS_val_Fo8))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 321,
   "metadata": {},
   "outputs": [],
   "source": [
    "    \n",
    "########################################\n",
    "# Moller-Plesset perturbation theories #\n",
    "########################################\n",
    "\n",
    "#arrays to store original data, 2-body energies, and total final MP2 contributions\n",
    "MP2_orig=np.zeros((2,20,3),dtype=np.float64)\n",
    "MP2_val=np.zeros((3,20),dtype=np.float64)\n",
    "MP2_tot=np.zeros((20))\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(20):\n",
    "        file='data/dimer/MP2/'+basis+'/'+str(\"{:02d}\".format(i+1))+'_Dimer_MP2_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_MP2_Molpro(filelist)\n",
    "                MP2_orig[basid[basis],i,:]=data\n",
    "                MP2_val[basid[basis],i]=mbe.get_2body(data)*ha2kc\n",
    "\n",
    "# do CBS extrapolation\n",
    "for i in range(20):\n",
    "    MP2_val[basid['CBS'],i]=mbe.basis_extrap(MP2_val[basid['TZ'],i],MP2_val[basid['QZ'],i],3,4)\n",
    "    #use a sum of MP2/AVTZ->AVQZ and HF/AV5Z\n",
    "    MP2_tot[i]=MP2_val[basid['CBS'],i]+HF2_val_Q5[1,i]\n",
    "\n",
    "#MP2 corresponds to first index in MP_mbe array\n",
    "MP_mbe[0,1]=np.sum(MP2_tot[:])\n",
    "\n",
    "# do the same with MP3\n",
    "MP3_orig=np.zeros((2,20,3),dtype=np.float64)\n",
    "MP3_val=np.zeros((3,20),dtype=np.float64)\n",
    "MP3_tot=np.zeros((20))\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(20):\n",
    "        file='data/dimer/MP3/'+basis+'/'+str(\"{:02d}\".format(i+1))+'_Dimer_MP3_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_MP3_Molpro(filelist)\n",
    "                MP3_orig[basid[basis],i,:]=data\n",
    "                MP3_val[basid[basis],i]=mbe.get_2body(data)*ha2kc\n",
    "                \n",
    "for i in range(20):\n",
    "    MP3_val[basid['CBS'],i]=mbe.basis_extrap(MP3_val[basid['TZ'],i],MP3_val[basid['QZ'],i],3,4)\n",
    "    MP3_tot[i]=MP3_val[basid['CBS'],i]+HF2_val_Q5[1,i]\n",
    "\n",
    "#store the final MP3 value into MBE array\n",
    "MP_mbe[1,1]=np.sum(MP3_tot[:])\n",
    "\n",
    "# and finally MP4\n",
    "MP4_orig=np.zeros((2,20,3),dtype=np.float64)\n",
    "MP4_val=np.zeros((3,20),dtype=np.float64)\n",
    "MP4_tot=np.zeros((20))\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(20):\n",
    "        file='data/dimer/MP4/'+basis+'/'+str(\"{:02d}\".format(i+1))+'_Dimer_MP4_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_MP4_Molpro(filelist)\n",
    "                MP4_orig[basid[basis],i,:]=data\n",
    "                MP4_val[basid[basis],i]=mbe.get_2body(data)*ha2kc\n",
    "\n",
    "for i in range(20):\n",
    "    MP4_val[basid['CBS'],i]=mbe.basis_extrap(MP4_val[basid['TZ'],i],MP4_val[basid['QZ'],i],3,4)\n",
    "    MP4_tot[i]=MP4_val[basid['CBS'],i]+HF2_val_Q5[1,i]\n",
    "\n",
    "#store the final MP4 value into array for MBE \n",
    "MP_mbe[2,1]=np.sum(MP4_tot[:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 322,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read RPA data from output files\n",
    "#indices are functional, basis set, dimer, dimer/mono1/mono2\n",
    "DFT_orig=np.zeros((5,2,20,3),dtype=np.float64)\n",
    "#indices are functional, TZ/QZ/extrap, dimer\n",
    "DFT_val=np.zeros((5,3,20),dtype=np.float64)\n",
    "EXX_orig=np.zeros((5,2,20,3),dtype=np.float64)\n",
    "EXX_val=np.zeros((5,3,20),dtype=np.float64)\n",
    "RSE_orig=np.zeros((5,2,20,3),dtype=np.float64)\n",
    "RSE_val=np.zeros((5,3,20),dtype=np.float64)\n",
    "RPA_orig=np.zeros((5,2,20,3),dtype=np.float64)\n",
    "RPA_val=np.zeros((5,3,20),dtype=np.float64)\n",
    "\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0','HF']:\n",
    "    for basis in ['small-basis','large-basis']:\n",
    "        #loop over all the dimers\n",
    "        for i in range(20):\n",
    "            file='data/dimer/RPA/'+xc+'/'+basis+'/'+str(\"{:02d}\".format(i+1))+'.log'\n",
    "            if os.path.isfile(file):\n",
    "                with open(file,'r') as f:\n",
    "                    filelist=f.readlines()\n",
    "                    data=mbe.read_RPA_MMM(filelist)\n",
    "                    DFT_orig[methods[xc],MM_basis[basis],i,:]=data[0]\n",
    "                    DFT_val[methods[xc],MM_basis[basis],i]=mbe.get_2body(data[0])*ha2kc\n",
    "                    EXX_orig[methods[xc],MM_basis[basis],i,:]=data[1]\n",
    "                    EXX_val[methods[xc],MM_basis[basis],i]=mbe.get_2body(data[1])*ha2kc                    \n",
    "                    RSE_orig[methods[xc],MM_basis[basis],i,:]=data[2]\n",
    "                    RSE_val[methods[xc],MM_basis[basis],i]=mbe.get_2body(data[2])*ha2kc                    \n",
    "                    RPA_orig[methods[xc],MM_basis[basis],i,:]=data[3]\n",
    "                    RPA_val[methods[xc],MM_basis[basis],i]=mbe.get_2body(data[3])*ha2kc\n",
    "                    \n",
    "#extrapolate to CBS\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0','HF']:\n",
    "    for i in range(20):\n",
    "        RPA_val[methods[xc],MM_basis['CBS'],i]=mbe.basis_extrap(RPA_val[methods[xc],MM_basis['small-basis'],i],RPA_val[methods[xc],MM_basis['large-basis'],i],3,4)\n",
    "        \n",
    "#store the total RPA interaction energies\n",
    "#the indices are xc; rse or not; dimer idx\n",
    "RPA_tot=np.zeros((5,2,20),dtype=np.float64)\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0','HF']:\n",
    "    for rse in range(2):\n",
    "        for i in range(20):\n",
    "            RPA_tot[methods[xc],rse,i]=EXX_val[methods[xc],MM_basis['large-basis'],i]+RPA_val[methods[xc],MM_basis['CBS'],i]+rse*RSE_val[methods[xc],MM_basis['large-basis'],i]                    \n",
    "        RPA_mbe[methods[xc],rse,1]=np.sum(RPA_tot[methods[xc],rse,:])\n",
    "#print(RPA_mbe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 323,
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the oep RPA data\n",
    "\n",
    "RPA_oep=np.zeros((3,2,20),dtype=np.float64)\n",
    "\n",
    "for xc in ['oepx','oep2sc','hf']:\n",
    "    with open('data/dimer/RPA_oep/'+xc+'.dat') as f:\n",
    "        fli=f.readlines()\n",
    "    for i in range(20):\n",
    "        RPA_oep[oep_m[xc],1,i]=float(fli[i].split()[0])\n",
    "        RPA_oep[oep_m[xc],0,i]=float(fli[i].split()[1])\n",
    "        \n",
    "#print(RPA_oep[oep_m['oep2sc'],0,0:4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Three-body data\n",
    "\n",
    "We have\n",
    "* CCSD(T) with orbital=1d-8,energy=1d-8,twoint=1.d-19\n",
    "* CCSD(T)-F12 with orbital=1d-8,energy=1d-8,twoint=1.d-19, THRF12=1.d-12\n",
    "* MP2, MP3, MP4 with orbital=1d-8,energy=1d-8,twoint=1.d-19, no other setting, checked that differences to 10^-10 are on 4th place in kcal/mol\n",
    "* RPA with AVTZ and AVQZ basis, PBE, PBE0, SCAN, SCAN0, HF, ludicrous settings\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {},
   "outputs": [],
   "source": [
    "#################################################\n",
    "# Canonical CCSD(T) in AVTZ and AVQZ basis sets #\n",
    "#################################################\n",
    "\n",
    "#Hartree-Fock original data, indices: basis set; index; parts of fragment (trimer, monomers, dimers)\n",
    "HF3_orig=np.zeros((2,190,7),dtype=np.float64)\n",
    "#Hartree-Fock interaction energies, indices: basis set; index; three-body energy, two-body energies\n",
    "HF3_val=np.zeros((3,190,4),dtype=np.float64)\n",
    "#CCSD correlation energy\n",
    "CCSD3_orig=np.zeros((2,190,7),dtype=np.float64)\n",
    "CCSD3_val=np.zeros((3,190,4),dtype=np.float64)\n",
    "#CCSD(T) correlation energy\n",
    "CCSDT3_orig=np.zeros((2,190,7),dtype=np.float64)\n",
    "CCSDT3_val=np.zeros((3,190,4),dtype=np.float64)\n",
    "#(T) contributions\n",
    "T3_val=np.zeros((3,190,4),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(190):\n",
    "        file='data/trimer/CC_orb8/'+basis+'/'+str(\"{:03d}\".format(i+1))+'_CC_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_old_HF_Molpro(filelist)\n",
    "                HF3_orig[basid[basis],i,:]=data\n",
    "                inter=mbe.get_3body_312(data)\n",
    "                #loop over the trimer and three dimer energies and store them in the final array\n",
    "                for clus in range(4):\n",
    "                    HF3_val[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                    \n",
    "                data=mbe.read_CC_F12_Molpro(filelist)[4]\n",
    "                inter=mbe.get_3body_312(data)\n",
    "                CCSD3_orig[basid[basis],i,:]=data\n",
    "                for clus in range(4):\n",
    "                    CCSD3_val[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                    \n",
    "                data=mbe.read_MP4_Molpro(filelist)\n",
    "                CCSDT3_orig[basid[basis],i,:]=data\n",
    "                inter=mbe.get_3body_312(data)\n",
    "                for clus in range(4):\n",
    "                    CCSDT3_val[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                    T3_val[basid[basis],i,clus]=CCSDT3_val[basid[basis],i,clus]-CCSD3_val[basid[basis],i,clus]\n",
    "                    \n",
    "CCSD3_tot=np.zeros((190))\n",
    "CCSDT3_tot=np.zeros((190))\n",
    "for i in range(190):\n",
    "    for clus in range(4):\n",
    "        CCSD3_val[basid['CBS'],i,clus]=mbe.basis_extrap(CCSD3_val[basid['TZ'],i,clus],CCSD3_val[basid['QZ'],i,clus],3,4)\n",
    "        CCSDT3_val[basid['CBS'],i,clus]=mbe.basis_extrap(CCSDT3_val[basid['TZ'],i,clus],CCSDT3_val[basid['QZ'],i,clus],3,4)\n",
    "        T3_val[basid['CBS'],i,clus]=mbe.basis_extrap(T3_val[basid['TZ'],i,clus],T3_val[basid['QZ'],i,clus],3,4)\n",
    "    CCSD3_tot[i]=CCSD3_val[basid['CBS'],i,0]+HF3_val[basid['QZ'],i,0]\n",
    "    CCSDT3_tot[i]=CCSDT3_val[basid['CBS'],i,0]+HF3_val[basid['QZ'],i,0]\n",
    "\n",
    "#print(CCSD_val_o8[2,10:14])\n",
    "#print(CCSD3_val_o8[2,10:14,0])\n",
    "\n",
    "#HF+CABS data\n",
    "CABS3_orig_F=np.zeros((2,190,7),dtype=np.float64)\n",
    "CABS3_val_F=np.zeros((3,190,4),dtype=np.float64)\n",
    "#CCSD \n",
    "CCSD3_orig_F=np.zeros((2,190,7),dtype=np.float64)\n",
    "CCSD3_val_F=np.zeros((3,190,4),dtype=np.float64)\n",
    "#CCSD-F12b\n",
    "CCSDF12b3_orig_F=np.zeros((2,190,7),dtype=np.float64)\n",
    "CCSDF12b3_val_F=np.zeros((3,190,4),dtype=np.float64)\n",
    "#scaled Triples\n",
    "Ts3_orig_F=np.zeros((2,190,7),dtype=np.float64)\n",
    "Ts3_val_F=np.zeros((3,190,4),dtype=np.float64)\n",
    "#unscaled Triples\n",
    "Tu3_orig_F=np.zeros((2,190,7),dtype=np.float64)\n",
    "Tu3_val_F=np.zeros((3,190,4),dtype=np.float64)\n",
    "CCSDF12b3_tot_F=np.zeros((190))\n",
    "CCSDTF12b3_tot_F=np.zeros((190))\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(190):\n",
    "        file='data/trimer/CCF12_orb8/'+basis+'/'+str(\"{:03d}\".format(i+1))+'_F12_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_CC_F12_Molpro(filelist)\n",
    "                if len(data[1]) != 7:\n",
    "                    #things are OK now but this was not always the case\n",
    "                    print('CC 3-body output file corrupted !')\n",
    "                    print(basis, i, data, len(data[1]))\n",
    "                else:\n",
    "                    # return HF, CABS, MP2, MF12, CC, CCF12b, Ts, Tu, F12a, F12b\n",
    "                    #sort out the data into what we need\n",
    "                    CABS3_orig_F[basid[basis],i,:]=data[1]\n",
    "                    CCSD3_orig_F[basid[basis],i,:]=data[4] \n",
    "                    CCSDF12b3_orig_F[basid[basis],i,:]=data[5]\n",
    "                    Ts3_orig_F[basid[basis],i,:]=data[6]\n",
    "                    Tu3_orig_F[basid[basis],i,:]=data[7]\n",
    "                    ##print (data)\n",
    "                    inter=mbe.get_3body(data[1])\n",
    "                    for clus in range(4):\n",
    "                        CABS3_val_F[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                    inter=mbe.get_3body(data[4])\n",
    "                    for clus in range(4):\n",
    "                        CCSD3_val_F[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                    inter=mbe.get_3body(data[5])\n",
    "                    for clus in range(4):\n",
    "                        CCSDF12b3_val_F[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                    inter=mbe.get_3body(data[6])\n",
    "                    for clus in range(4):\n",
    "                        Ts3_val_F[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                    inter=mbe.get_3body(data[7])\n",
    "                    for clus in range(4):\n",
    "                        Tu3_val_F[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                \n",
    "for i in range(190):\n",
    "    for clus in range(4):\n",
    "        CCSD3_val_F[basid['CBS'],i,clus]=mbe.basis_extrap(CCSD3_val_F[basid['TZ'],i,clus],CCSD3_val_F[basid['QZ'],i,clus],3,4)\n",
    "        #Ts extrapolate to CBS? \n",
    "        #Ts3_val_Fo8[basid['CBS'],i,clus]=mbe.basis_extrap(Ts3_val_Fo8[basid['TZ'],i,clus],Ts3_val_Fo8[basid['QZ'],i,clus],3,4)\n",
    "        Tu3_val_F[basid['CBS'],i,clus]=mbe.basis_extrap(Tu3_val_F[basid['TZ'],i,clus],Tu3_val_F[basid['QZ'],i,clus],3,4)\n",
    "    CCSDF12b3_tot_F[i]=CCSDF12b3_val_F[basid['QZ'],i,0]+CABS3_val_F[basid['QZ'],i,0]\n",
    "    CCSDTF12b3_tot_F[i]=CCSDF12b3_tot_F[i]+Tu3_val_F[basid['CBS'],i,0]\n",
    "    \n",
    "#print(np.sum(CCSDTF123_tot_Fo8[:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 325,
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################\n",
    "# Moller-Plesset three-body #\n",
    "#############################\n",
    "\n",
    "#MP2\n",
    "#original data, indices: basis set; trimer index; fragment\n",
    "MP23_orig=np.zeros((2,190,7),dtype=np.float64)\n",
    "#interaction energies, indices: basis set; trimer index; fragment\n",
    "MP23_val=np.zeros((3,190,4),dtype=np.float64)\n",
    "\n",
    "#three-body nonadditive MP2 interaction energy, index: trimer index\n",
    "MP23_tot=np.zeros((190),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the trimer\n",
    "    for i in range(190):\n",
    "        file='data/trimer/MP2/'+basis+'/'+str(\"{:03d}\".format(i+1))+'_MP2_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_MP2_Molpro(filelist)\n",
    "                if len(data) != 7:\n",
    "                    print('MP3 3-body output file corrupted !')\n",
    "                    print(basis, i, data, len(data))\n",
    "                else:\n",
    "                    MP23_orig[basid[basis],i,:]=data\n",
    "                    inter=mbe.get_3body_312(data)\n",
    "                    for clus in range(4):\n",
    "                        MP23_val[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                \n",
    "for i in range(190):\n",
    "    for clus in range(4):\n",
    "        MP23_val[basid['CBS'],i,clus]=mbe.basis_extrap(MP23_val[basid['TZ'],i,clus],MP23_val[basid['QZ'],i,clus],3,4)    \n",
    "    MP23_tot[i]=MP23_val[basid['CBS'],i,0]+HF3_val[basid['QZ'],i,0]\n",
    "\n",
    "#MP3\n",
    "#arrays same as for MP2\n",
    "MP33_orig=np.zeros((2,190,7),dtype=np.float64)\n",
    "MP33_val=np.zeros((3,190,4),dtype=np.float64)\n",
    "MP33_tot=np.zeros((190),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(190):\n",
    "        file='data/trimer/MP3/'+basis+'/'+str(\"{:03d}\".format(i+1))+'_MP3_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_MP3_Molpro(filelist)\n",
    "                if len(data) != 7:\n",
    "                    print('Value corrupted !')\n",
    "                    print(basis, i, data, len(data))\n",
    "                else:\n",
    "                    MP33_orig[basid[basis],i,:]=data\n",
    "                    inter=mbe.get_3body_312(data)\n",
    "                    for clus in range(4):\n",
    "                        MP33_val[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                \n",
    "for i in range(190):\n",
    "    for clus in range(4):\n",
    "        MP33_val[basid['CBS'],i,clus]=mbe.basis_extrap(MP33_val[basid['TZ'],i,clus],MP33_val[basid['QZ'],i,clus],3,4) \n",
    "    MP33_tot[i]=MP33_val[basid['CBS'],i,0]+HF3_val[basid['QZ'],i,0]\n",
    "        \n",
    "#MP4\n",
    "#arrays same as for MP2\n",
    "MP43_orig=np.zeros((2,190,7),dtype=np.float64)\n",
    "MP43_val=np.zeros((3,190,4),dtype=np.float64)\n",
    "MP43_tot=np.zeros((190),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    #loop over all the dimers\n",
    "    for i in range(190):\n",
    "        file='data/trimer/MP4/'+basis+'/'+str(\"{:03d}\".format(i+1))+'_MP4_'+basis+'.out'\n",
    "        if os.path.isfile(file):\n",
    "            with open(file,'r') as f:\n",
    "                filelist=f.readlines()\n",
    "                data=mbe.read_MP4_Molpro(filelist)\n",
    "                if len(data) != 7:\n",
    "                    print('Value corrupted !')\n",
    "                    print(basis, i, data, len(data))\n",
    "                else:\n",
    "                    MP43_orig[basid[basis],i,:]=data\n",
    "                    inter=mbe.get_3body_312(data)\n",
    "                    for clus in range(4):\n",
    "                        MP43_val[basid[basis],i,clus]=inter[clus]*ha2kc\n",
    "                \n",
    "for i in range(190):\n",
    "    for clus in range(4):\n",
    "        MP43_val[basid['CBS'],i,clus]=mbe.basis_extrap(MP43_val[basid['TZ'],i,clus],MP43_val[basid['QZ'],i,clus],3,4)            \n",
    "    MP43_tot[i]=MP43_val[basid['CBS'],i,0]+HF3_val[basid['QZ'],i,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 361,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################\n",
    "# Read RPA three-body energies #\n",
    "################################\n",
    "\n",
    "#indices are functional; basis set; trimer index; fragment part\n",
    "DFT3_orig=np.zeros((5,2,190,7),dtype=np.float64)\n",
    "#indices are functional; basis set TZ/QZ/extrap; interaction energy index (three-body and three two-body)\n",
    "DFT3_val=np.zeros((5,3,190,4),dtype=np.float64)\n",
    "EXX3_orig=np.zeros((5,2,190,7),dtype=np.float64)\n",
    "EXX3_val=np.zeros((5,3,190,4),dtype=np.float64)\n",
    "RSE3_orig=np.zeros((5,2,190,7),dtype=np.float64)\n",
    "RSE3_val=np.zeros((5,3,190,4),dtype=np.float64)\n",
    "RPA3_orig=np.zeros((5,2,190,7),dtype=np.float64)\n",
    "RPA3_val=np.zeros((5,3,190,4),dtype=np.float64)\n",
    "\n",
    "#final non-additive three-body RPA energies\n",
    "#indices are functional; no RSE or yes; trimer index \n",
    "RPA3_tot=np.zeros((5,2,190),dtype=np.float64)\n",
    "\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0','HF']:\n",
    "    for basis in ['small-basis','large-basis']:\n",
    "        #loop over all the trimers\n",
    "        for i in range(190):\n",
    "            file='data/trimer/RPA/'+xc+'/'+basis+'/'+str(\"{:03d}\".format(i+1))+'.log'\n",
    "            if os.path.isfile(file):\n",
    "                with open(file,'r') as f:\n",
    "                    filelist=f.readlines()\n",
    "                    data=mbe.read_RPA_MMM(filelist)\n",
    "                    DFT3_orig[methods[xc],MM_basis[basis],i,:]=data[0]\n",
    "                    #print(data[0], 'RPA')\n",
    "                    inter=mbe.get_3body_312(data[0])\n",
    "                    for clus in range(4):\n",
    "                        DFT3_val[methods[xc],MM_basis[basis],i,clus]=inter[clus]*ha2kc\n",
    "                    EXX3_orig[methods[xc],MM_basis[basis],i,:]=data[1]\n",
    "                    inter=mbe.get_3body_312(data[1])\n",
    "                    for clus in range(4):\n",
    "                        EXX3_val[methods[xc],MM_basis[basis],i,clus]=inter[clus]*ha2kc   \n",
    "                    RSE3_orig[methods[xc],MM_basis[basis],i,:]=data[2]\n",
    "                    inter=mbe.get_3body_312(data[2])\n",
    "                    for clus in range(4):\n",
    "                        RSE3_val[methods[xc],MM_basis[basis],i,clus]=inter[clus]*ha2kc                   \n",
    "                    RPA3_orig[methods[xc],MM_basis[basis],i,:]=data[3]\n",
    "                    inter=mbe.get_3body_312(data[3])\n",
    "                    for clus in range(4):\n",
    "                        RPA3_val[methods[xc],MM_basis[basis],i,clus]=inter[clus]*ha2kc\n",
    "\n",
    "#print(np.sum(DFT3_val[0,0,:,0]))\n",
    "\n",
    "#do extrapolation of the correlation energy\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0','HF']:\n",
    "    for i in range(190):\n",
    "        for clus in range(4):\n",
    "            RPA3_val[methods[xc],MM_basis['CBS'],i,clus]=mbe.basis_extrap(RPA3_val[methods[xc],MM_basis['small-basis'],i,clus],RPA3_val[methods[xc],MM_basis['large-basis'],i,clus],3,4)\n",
    "\n",
    "#calculate the total RPA energy without and with singles corrections\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0','HF']:\n",
    "    for rse in range(2):\n",
    "        for i in range(190):\n",
    "            RPA3_tot[methods[xc],rse,i]=EXX3_val[methods[xc],MM_basis['large-basis'],i,0] +RPA3_val[methods[xc],MM_basis['CBS'],i,0]+rse*RSE3_val[methods[xc],MM_basis['large-basis'],i,0]\n",
    "            \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 365,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Divide into groups according to distance\n",
    "# CABS3_val_F[3,190,4]\n",
    "# CCSDF12b3_val_F\n",
    "# Tu_val_F\n",
    "# CCSDTF12b3_tot_F\n",
    "#\n",
    "# MP23_tot\n",
    "# MP33_tot\n",
    "# MP43_tot\n",
    "#\n",
    "# RPA4_tot\n",
    "\n",
    "#declare the arrays\n",
    "CABS3_val_F_NN=np.zeros((190,4))\n",
    "CCSDF12b3_val_F_NN=np.zeros((190,4))\n",
    "Tu3_val_F_NN=np.zeros((190,4))\n",
    "CCSDTF12b3_tot_F_NN=np.zeros((190,4))\n",
    "\n",
    "MP23_val_NN=np.zeros((190,4))\n",
    "MP33_val_NN=np.zeros((190,4))\n",
    "MP43_val_NN=np.zeros((190,4))\n",
    "\n",
    "RPA3_tot_NN=np.zeros((5,2,190,4))\n",
    "\n",
    "idx3_ok=np.zeros((4),dtype=np.int16)\n",
    "#distribute the data into arrays \n",
    "#masking did not work for some stuff\n",
    "for i in range(190):\n",
    "    nn=trim_nn_case[i]\n",
    "    \n",
    "    #first index is the basis set, for CABS and CCSDF12b we take directly the AVQZ value\n",
    "    #for Tu we take the AVTZ->AVQZ value\n",
    "    CABS3_val_F_NN[idx3_ok[nn],nn]=CABS3_val_F[1,i,0]\n",
    "    CCSDF12b3_val_F_NN[idx3_ok[nn],nn]=CCSDF12b3_val_F[1,i,0]\n",
    "    Tu3_val_F_NN[idx3_ok[nn],nn]=Tu3_val_F[2,i,0]\n",
    "    CCSDTF12b3_tot_F_NN[idx3_ok[nn],nn]=CCSDTF12b3_tot_F[i]\n",
    "        \n",
    "    MP23_val_NN[idx3_ok[nn],nn]=MP23_tot[i]\n",
    "    MP33_val_NN[idx3_ok[nn],nn]=MP33_tot[i]\n",
    "    MP43_val_NN[idx3_ok[nn],nn]=MP43_tot[i]\n",
    "    \n",
    "    for xc in range(5):\n",
    "        for rse in range(2):\n",
    "            RPA3_tot_NN[xc,rse,idx3_ok[nn],nn]=RPA3_tot[xc,rse,i]\n",
    "    idx3_ok[nn]+=1\n",
    "#print(np.sum(RPA3_tot_NN[0,0,:,0]))\n",
    "#print(np.sum(MP23_val_NN[:,0])) \n",
    "\n",
    "\n",
    "#store the contributions into MBE arrays\n",
    "#MP2 in 3b 1hb group contains only fragments with neighbouring molecules\n",
    "MP_mbe[0,2]=np.sum(MP23_val_NN[:,0])\n",
    "#MP2 in 3b 0hb group contains larger distances\n",
    "MP_mbe[0,3]=np.sum(MP23_val_NN[:,1:])\n",
    "MP_mbe[1,2]=np.sum(MP33_val_NN[:,0])  #same but for MP3\n",
    "MP_mbe[1,3]=np.sum(MP33_val_NN[:,1:])  \n",
    "MP_mbe[2,2]=np.sum(MP43_val_NN[:,0])  #same but for MP4\n",
    "MP_mbe[2,3]=np.sum(MP43_val_NN[:,1:])\n",
    "\n",
    "for xc in range(5):\n",
    "    for rse in range(2):\n",
    "        RPA_mbe[xc,rse,2]=np.sum(RPA3_tot_NN[xc,rse,:,0])\n",
    "        RPA_mbe[xc,rse,3]=np.sum(RPA3_tot_NN[xc,rse,:,1:])\n",
    "\n",
    "#store the total contributions to MBE, HF, CC-F12, (T), total\n",
    "#3B near\n",
    "CC_mbeQ5[0,2]=np.sum(CABS3_val_F_NN[:,0])\n",
    "CC_mbeQ5[1,2]=np.sum(CCSDF12b3_val_F_NN[:,0])\n",
    "CC_mbeQ5[2,2]=np.sum(Tu3_val_F_NN[:,0])\n",
    "CC_mbeQ5[3,2]=np.sum(CCSDTF12b3_tot_F_NN[:,0])\n",
    "#3B far\n",
    "CC_mbeQ5[0,3]=np.sum(CABS3_val_F_NN[:,1:])\n",
    "CC_mbeQ5[1,3]=np.sum(CCSDF12b3_val_F_NN[:,1:])\n",
    "CC_mbeQ5[2,3]=np.sum(Tu3_val_F_NN[:,1:])\n",
    "CC_mbeQ5[3,3]=np.sum(CCSDTF12b3_tot_F_NN[:,1:])\n",
    "\n",
    "#print(CC_mbeQ5)\n",
    "#print(RPA_mbe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Four-body interactions\n",
    "\n",
    "We have\n",
    "* CCSD(T)-F12, AVDZ and AVTZ basis set, orbital=1d-10, energy=1d-10, accu=1d-16, THRAOF12=1d-10\n",
    "* MP2 and MP3, AVTZ and AVQZ basis set, orbital=1d-10, energy=1d-10, accu=1d-16\n",
    "* RPA, AVTZ and AVQZ, PBE, PBE0, SCAN, SCAN0\n",
    "\n",
    "Results are stripped to datafiles to save space.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 328,
   "metadata": {},
   "outputs": [],
   "source": [
    "#we have a helper routine to read the stripped datafiles\n",
    "\n",
    "def readfile_4b_4123(path,data_orig,data_val,bas):\n",
    "    #file='data/tetramer/CC/AV'+basis+'_no_acc_in_HF/cabs_av'+basis.lower()+'_tot.dat'\n",
    "    with open(file,'r') as f:\n",
    "        filelist=f.readlines()\n",
    "    i=0\n",
    "    for line in filelist:\n",
    "        #print(len(line.split()), line.split()[0])\n",
    "        data_orig[bas,i,:]=line.split()[3:]\n",
    "        inter=mbe.get_4body_4123(data_orig[bas,i,:])\n",
    "        for clus in range(12):\n",
    "            data_val[bas,i,clus]=inter[clus]*ha2kc\n",
    "        i+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 329,
   "metadata": {},
   "outputs": [],
   "source": [
    "######\n",
    "# CC #\n",
    "######\n",
    "\n",
    "#data with accu=10e-16 in hf calculation and orbital=1.e-10\n",
    "#HF+CABS data\n",
    "HF4_orig_F=np.zeros((2,1140,15),dtype=np.float64)\n",
    "HF4_val_F=np.zeros((3,1140,12),dtype=np.float64)\n",
    "#HF+CABS data\n",
    "CABS4_orig_F=np.zeros((2,1140,15),dtype=np.float64)\n",
    "CABS4_val_F=np.zeros((3,1140,12),dtype=np.float64)\n",
    "#CCSD \n",
    "CCSD4_orig_F=np.zeros((2,1140,15),dtype=np.float64)\n",
    "CCSD4_val_F=np.zeros((3,1140,12),dtype=np.float64)\n",
    "#CCSD-F12b\n",
    "CCSDF12b4_orig_F=np.zeros((2,1140,15),dtype=np.float64)\n",
    "CCSDF12b4_val_F=np.zeros((3,1140,12),dtype=np.float64)\n",
    "#scaled Triples\n",
    "Ts4_orig_F=np.zeros((2,1140,15),dtype=np.float64)\n",
    "Ts4_val_F=np.zeros((3,1140,12),dtype=np.float64)\n",
    "#unscaled Triples\n",
    "Tu4_orig_F=np.zeros((2,1140,15),dtype=np.float64)\n",
    "Tu4_val_F=np.zeros((3,1140,12),dtype=np.float64)\n",
    "#MP2 correlation\n",
    "MP24_orig_F=np.zeros((2,1140,15),dtype=np.float64)\n",
    "MP24_val_F=np.zeros((3,1140,12),dtype=np.float64)\n",
    "#MP2-F12 correlation\n",
    "MP2F124_orig_F=np.zeros((2,1140,15),dtype=np.float64)\n",
    "MP2F124_val_F=np.zeros((3,1140,12),dtype=np.float64)\n",
    "\n",
    "for basis in ['DZ','TZ']:\n",
    "    file='data/tetramer/CC/AV'+basis+'/hf_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,HF4_orig_F,HF4_val_F,basDTid[basis])\n",
    "    \n",
    "    file='data/tetramer/CC/AV'+basis+'/cabs_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,CABS4_orig_F,CABS4_val_F,basDTid[basis])\n",
    "    \n",
    "    file='data/tetramer/CC/AV'+basis+'/cc_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,CCSD4_orig_F,CCSD4_val_F,basDTid[basis])\n",
    "    \n",
    "    file='data/tetramer/CC/AV'+basis+'/ccf12b_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,CCSDF12b4_orig_F,CCSDF12b4_val_F,basDTid[basis])\n",
    "    \n",
    "    file='data/tetramer/CC/AV'+basis+'/ts_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,Ts4_orig_F,Ts4_val_F,basDTid[basis])\n",
    "    \n",
    "    file='data/tetramer/CC/AV'+basis+'/tu_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,Tu4_orig_F,Tu4_val_F,basDTid[basis])\n",
    "     \n",
    "    file='data/tetramer/CC/AV'+basis+'/mp2_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,MP24_orig_F,MP24_val_F,basDTid[basis])\n",
    "    \n",
    "    file='data/tetramer/CC/AV'+basis+'/f12_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,MP2F124_orig_F,MP2F124_val_F,basDTid[basis])\n",
    "\n",
    "CCSDF12b4_tot_F=np.zeros((1140),dtype=np.float64)\n",
    "CCSDTF12b4_tot_F=np.zeros((1140),dtype=np.float64)\n",
    "\n",
    "for i in range(1140):\n",
    "    #for clus in range(4):\n",
    "    #    CCSD3_val_F[basid['CBS'],i,clus]=mbe.basis_extrap(CCSD3_val_F[basid['TZ'],i,clus],CCSD3_val_F[basid['QZ'],i,clus],3,4)\n",
    "    #    #Ts extrapolate to CBS? \n",
    "    #    #Ts3_val_Fo8[basid['CBS'],i,clus]=mbe.basis_extrap(Ts3_val_Fo8[basid['TZ'],i,clus],Ts3_val_Fo8[basid['QZ'],i,clus],3,4)\n",
    "    #    Tu3_val_F[basid['CBS'],i,clus]=mbe.basis_extrap(Tu3_val_F[basid['TZ'],i,clus],Tu3_val_F[basid['QZ'],i,clus],3,4)\n",
    "    CCSDF12b4_tot_F[i]=CCSDF12b4_val_F[basDTid['TZ'],i,0]+CABS4_val_F[basDTid['TZ'],i,0]\n",
    "    CCSDTF12b4_tot_F[i]=CCSDF12b4_tot_F[i]+Tu4_val_F[basDTid['TZ'],i,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 330,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5286505765036471\n"
     ]
    }
   ],
   "source": [
    "##################\n",
    "# Moller-Plesset #\n",
    "##################\n",
    "\n",
    "#this is without screening settings, same as CC calculation Fo10 for tetramers\n",
    "HF4mp2_orig=np.zeros((2,1140,15),dtype=np.float64)\n",
    "HF4mp2_val=np.zeros((3,1140,12),dtype=np.float64)\n",
    "MP24_orig=np.zeros((2,1140,15),dtype=np.float64)\n",
    "MP24_val=np.zeros((3,1140,12),dtype=np.float64)\n",
    "\n",
    "MP24_tot=np.zeros((1140),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    file='data/tetramer/MP2/AV'+basis+'_o10/hf_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,HF4mp2_orig,HF4mp2_val,basid[basis])\n",
    "    \n",
    "    file='data/tetramer/MP2/AV'+basis+'_o10/mp2_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,MP24_orig,MP24_val,basid[basis])\n",
    "    \n",
    "for i in range(1140):\n",
    "    for clus in range(12):\n",
    "        MP24_val[basid['CBS'],i,clus]=mbe.basis_extrap(MP24_val[basid['TZ'],i,clus],MP24_val[basid['QZ'],i,clus],3,4) \n",
    "    #Store the total MP2 4B contributions for individual tetramers\n",
    "    MP24_tot[i]=HF4mp2_val[basid['QZ'],i,0]+MP24_val[basid['CBS'],i,0]\n",
    "    \n",
    "\n",
    "#MP3 same_as_CC\n",
    "HF4_orig=np.zeros((2,1140,15),dtype=np.float64)\n",
    "HF4_val=np.zeros((3,1140,12),dtype=np.float64)\n",
    "MP34_orig=np.zeros((2,1140,15),dtype=np.float64)\n",
    "MP34_val=np.zeros((3,1140,12),dtype=np.float64)\n",
    "MP34_tot=np.zeros((1140),dtype=np.float64)\n",
    "\n",
    "for basis in ['TZ','QZ']:\n",
    "    file='data/tetramer/MP3/AV'+basis+'_o10/hf_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,HF4_orig,HF4_val,basid[basis])\n",
    "    \n",
    "    file='data/tetramer/MP3/AV'+basis+'_o10/mp3_av'+basis.lower()+'_tot.dat'\n",
    "    readfile_4b_4123(file,MP34_orig,MP34_val,basid[basis])\n",
    "    \n",
    "for i in range(1140):\n",
    "    for clus in range(12):\n",
    "        MP34_val[basid['CBS'],i,clus]=mbe.basis_extrap(MP34_val[basid['TZ'],i,clus],MP34_val[basid['QZ'],i,clus],3,4) \n",
    "    #Store the total MP3 4B contributions for individual tetramers\n",
    "    MP34_tot[i]=HF4_val[basid['QZ'],i,0]+MP34_val[basid['CBS'],i,0]\n",
    "print(np.sum(MP34_tot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 331,
   "metadata": {},
   "outputs": [],
   "source": [
    "##################### \n",
    "# RPA for four-body #\n",
    "#####################\n",
    "\n",
    "#indices are functional; basis set; tetramer index; fragment\n",
    "DFT4_orig=np.zeros((5,2,1140,15),dtype=np.float64)\n",
    "#indices are functional; TZ/QZ/extrap; interaction energies of various fragments starting with the 4-body\n",
    "DFT4_val=np.zeros((5,3,1140,12),dtype=np.float64)\n",
    "EXX4_orig=np.zeros((5,2,1140,15),dtype=np.float64)\n",
    "EXX4_val=np.zeros((5,3,1140,12),dtype=np.float64)\n",
    "RSE4_orig=np.zeros((5,2,1140,15),dtype=np.float64)\n",
    "RSE4_val=np.zeros((5,3,1140,12),dtype=np.float64)\n",
    "RPA4_orig=np.zeros((5,2,1140,15),dtype=np.float64)\n",
    "RPA4_val=np.zeros((5,3,1140,12),dtype=np.float64)\n",
    "\n",
    "RPA4_tot=np.zeros((5,2,1140),dtype=np.float64)\n",
    "\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0']:\n",
    "    for basis in ['small-basis','large-basis']:\n",
    "        #loop over all the trimers\n",
    "        for i in range(1140):\n",
    "            file='data/tetramer/RPA/'+xc+'/'+basis+'/'+str(\"{:04d}\".format(i+1))+'.log'\n",
    "            if os.path.isfile(file):\n",
    "                with open(file,'r') as f:\n",
    "                    filelist=f.readlines()\n",
    "                    data=mbe.read_RPA_MMM(filelist)\n",
    "                    if len(data[0])!=15 or len(data[1])!=15 or len(data[2])!=15 or len(data[3])!=15:\n",
    "                        print(xc, basis, i, ' corrupted')\n",
    "                        print(len(data[0]),len(data[1]),len(data[2]),len(data[3]))\n",
    "                    else:\n",
    "                        DFT4_orig[methods[xc],MM_basis[basis],i,:]=data[0]\n",
    "                        #print(data[0], 'RPA')\n",
    "                        inter=mbe.get_4body_MMM(data[0])\n",
    "                        for clus in range(12):\n",
    "                            DFT4_val[methods[xc],MM_basis[basis],i,clus]=inter[clus]*ha2kc\n",
    "                            \n",
    "                        EXX4_orig[methods[xc],MM_basis[basis],i,:]=data[1]\n",
    "                        inter=mbe.get_4body_MMM(data[1])\n",
    "                        for clus in range(12):\n",
    "                            EXX4_val[methods[xc],MM_basis[basis],i,clus]=inter[clus]*ha2kc   \n",
    "                            \n",
    "                        RSE4_orig[methods[xc],MM_basis[basis],i,:]=data[2]\n",
    "                        inter=mbe.get_4body_MMM(data[2])\n",
    "                        for clus in range(12):\n",
    "                            RSE4_val[methods[xc],MM_basis[basis],i,clus]=inter[clus]*ha2kc                   \n",
    "                            \n",
    "                        RPA4_orig[methods[xc],MM_basis[basis],i,:]=data[3]\n",
    "                        inter=mbe.get_4body_MMM(data[3])\n",
    "                        for clus in range(12):\n",
    "                            RPA4_val[methods[xc],MM_basis[basis],i,clus]=inter[clus]*ha2kc\n",
    "                \n",
    "\n",
    "#do CBS extrapolation\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0','HF']:\n",
    "    for i in range(1140):\n",
    "        for clus in range(12):\n",
    "            RPA4_val[methods[xc],MM_basis['CBS'],i,clus]=mbe.basis_extrap(RPA4_val[methods[xc],MM_basis['small-basis'],i,clus],RPA4_val[methods[xc],MM_basis['large-basis'],i,clus],3,4)\n",
    "        RPA4_tot[methods[xc],0,i]=RPA4_val[methods[xc],MM_basis['CBS'],i,0]+EXX4_val[methods[xc],MM_basis['large-basis'],i,0]\n",
    "        RPA4_tot[methods[xc],1,i]=RPA4_tot[methods[xc],0,i]+RSE4_val[methods[xc],MM_basis['large-basis'],i,0]\n",
    "        \n",
    "#print(np.sum(RPA4_val[0,2,:,0])+np.sum(EXX4_val[0,1,:,0]))            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Divide tetramers into groups according to distance\n",
    "\n",
    "# CABS4_val_F[3,1140,4]\n",
    "# CCSDF12b4_val_F\n",
    "# Tu4_val_F\n",
    "# CCSDTF12b4_tot_F\n",
    "#\n",
    "# MP24_tot\n",
    "# MP34_tot\n",
    "# \n",
    "#\n",
    "# RPA4_tot\n",
    "\n",
    "#declare the arrays\n",
    "CABS4_val_F_NN=np.zeros((1140,3))\n",
    "CCSDF12b4_val_F_NN=np.zeros((1140,3))\n",
    "Tu4_val_F_NN=np.zeros((1140,3))\n",
    "CCSDTF12b4_tot_F_NN=np.zeros((1140,3))\n",
    "\n",
    "MP24_val_NN=np.zeros((1140,3))\n",
    "MP34_val_NN=np.zeros((1140,3))\n",
    "#no MP43_val_NN=np.zeros((190,4))\n",
    "\n",
    "RPA4_tot_NN=np.zeros((5,2,1140,3))\n",
    "\n",
    "idx4_ok=np.zeros((3),dtype=np.int16)\n",
    "#distribute the data into arrays \n",
    "#masking did not work for some stuff\n",
    "for i in range(1140):\n",
    "    nn=2-tetr_nn_num[i]\n",
    "    \n",
    "    #                              we take the AVTZ value for ith tetramer \n",
    "    CABS4_val_F_NN[idx4_ok[nn],nn]=CABS4_val_F[1,i,0]\n",
    "    CCSDF12b4_val_F_NN[idx4_ok[nn],nn]=CCSDF12b4_val_F[1,i,0]\n",
    "    Tu4_val_F_NN[idx4_ok[nn],nn]=Tu4_val_F[1,i,0]\n",
    "    CCSDTF12b4_tot_F_NN[idx4_ok[nn],nn]=CCSDTF12b4_tot_F[i]\n",
    "        \n",
    "    MP24_val_NN[idx4_ok[nn],nn]=MP24_tot[i]\n",
    "    MP34_val_NN[idx4_ok[nn],nn]=MP34_tot[i]\n",
    "    #MP44_val_NN[idx3_ok[nn],nn]=MP43_tot[i]\n",
    "    \n",
    "    for xc in range(5):\n",
    "        for rse in range(2):\n",
    "            RPA4_tot_NN[xc,rse,idx4_ok[nn],nn]=RPA4_tot[xc,rse,i]\n",
    "    idx4_ok[nn]+=1\n",
    "    \n",
    "#stored the total contributions in the mbe arrays\n",
    "for nn in range(3):\n",
    "    CC_mbeQ5[0,4+nn]=np.sum(CABS4_val_F_NN[:,nn])\n",
    "    CC_mbeQ5[1,4+nn]=np.sum(CCSDF12b4_val_F_NN[:,nn])\n",
    "    CC_mbeQ5[2,4+nn]=np.sum(Tu4_val_F_NN[:,nn])\n",
    "    CC_mbeQ5[3,4+nn]=np.sum(CCSDTF12b4_tot_F_NN[:,nn])\n",
    "    MP_mbe[0,4+nn]=np.sum(MP24_val_NN[:,nn])\n",
    "    MP_mbe[1,4+nn]=np.sum(MP34_val_NN[:,nn])\n",
    "    for xc in range(5):\n",
    "        for rse in range(2):\n",
    "            RPA_mbe[xc,rse,4+nn]=np.sum(RPA4_tot_NN[xc,rse,:,nn])\n",
    "            \n",
    "#print(CC_mbeQ5)       \n",
    "#print(MP_mbe)\n",
    "#print(np.sum(MP_mbe[1,:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "### Scaling of triples in CCSD(T)-F12 \n",
    "\n",
    "For data in Tables 1, 2, and 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 333,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimers\n",
      "TZ\n",
      "-1.624379532737954 -1.6589257433826656\n",
      "QZ\n",
      "-1.637391434545468 -1.655427454087662\n",
      "5Z\n",
      "-1.637962775768483 -1.6485656806474744\n",
      "\n",
      "Trimers\n",
      "Scaled individually    Unscaled   Scalled with a common factor\n",
      "TZ\n",
      "0.046394688181151264 0.201295864293327 0.22115113733005973\n",
      "QZ\n",
      "0.11417086640714366 0.20154342620099752 0.21015706177525822\n",
      "\n",
      "Tetramers\n",
      "Scaled individually    Unscaled   Scalled with a common factor\n",
      "DZ\n",
      "-1.4795357199584864 0.019201716530162217 0.025365521081888463\n",
      "TZ\n",
      "0.8952191456874071 0.02098214910987935 0.023099723619988004\n"
     ]
    }
   ],
   "source": [
    "# Triples scaling in CC-F12\n",
    "\n",
    "#dimers\n",
    "Ts=np.zeros((3))\n",
    "print(\"Dimers\")\n",
    "for basis in ['TZ']:\n",
    "    print(basis)\n",
    "    tot_Ts=0.\n",
    "    tot_new_Ts=0.\n",
    "    for i in range(20):\n",
    "        #print(i)\n",
    "        #scale the (T) energies using the dimer scaling factor\n",
    "        for j in range(3):\n",
    "            Ts[j]=Tu2_orig_FTQ[basid[basis],i,j]*(Ts2_orig_FTQ[basid[basis],i,0]/Tu2_orig_FTQ[basid[basis],i,0])\n",
    "        #calculate the (T) contribution with common scaling and add to total\n",
    "        tot_new_Ts+=(Ts[0]-Ts[1]-Ts[2])*ha2kc\n",
    "\n",
    "        #just add the individually scaled values together\n",
    "        tot_Ts+=Ts2_val_FTQ[basid[basis],i]\n",
    "    #print the individual and total scaling values    \n",
    "    print(tot_Ts, tot_new_Ts)\n",
    "    \n",
    "for basis in ['QZ','5Z']:\n",
    "    print(basis)\n",
    "    tot_Ts=0.\n",
    "    tot_new_Ts=0.\n",
    "    for i in range(20):\n",
    "        #print(i)\n",
    "        #scale the (T) energies using the dimer scaling factor\n",
    "        for j in range(3):\n",
    "            Ts[j]=Tu2_orig_FQ5[bas5id[basis],i,j]*(Ts2_orig_FQ5[bas5id[basis],i,0]/Tu2_orig_FQ5[bas5id[basis],i,0])\n",
    "        #calculate the (T) contribution with common scaling and add to total\n",
    "        tot_new_Ts+=(Ts[0]-Ts[1]-Ts[2])*ha2kc\n",
    "\n",
    "        #just add the individually scaled values together\n",
    "        tot_Ts+=Ts2_val_FQ5[bas5id[basis],i]\n",
    "    #print the individual and total scaling values    \n",
    "    print(tot_Ts, tot_new_Ts)\n",
    "    #print('AVQZ->AV5Z estimate of CBS for unscaled T', np.sum(Tu2_val_FQ5[2,:]))\n",
    "    \n",
    "#trimers\n",
    "print(\"\\nTrimers\")\n",
    "print(\"Scaled individually    Unscaled   Scalled with a common factor\")\n",
    "Ts=np.zeros((7))\n",
    "for basis in ['TZ','QZ']:\n",
    "    print(basis)\n",
    "    tot_Ts=0.\n",
    "    tot_new_Ts=0.\n",
    "    tot_Tu=0.\n",
    "    for i in range(190):\n",
    "        #print(i)\n",
    "        for j in range(7):\n",
    "            Ts[j]=Tu3_orig_F[basid[basis],i,j]*(Ts3_orig_F[basid[basis],i,0]/Tu3_orig_F[basid[basis],i,0])\n",
    "        inter=mbe.get_3body(Ts)[0]\n",
    "        tot_Ts+=Ts3_val_F[basid[basis],i][0]\n",
    "        tot_Tu+=Tu3_val_F[basid[basis],i][0]\n",
    "        tot_new_Ts+=inter*ha2kc\n",
    "    print(tot_Ts, tot_Tu, tot_new_Ts)\n",
    "#OK, so in this case, triples scaled with a common factor help\n",
    "\n",
    "#tetramers\n",
    "Ts=np.zeros((15))\n",
    "print(\"\\nTetramers\")\n",
    "print(\"Scaled individually    Unscaled   Scalled with a common factor\")\n",
    "for basis in ['DZ','TZ']:\n",
    "    print(basis)\n",
    "    tot_Ts=0.\n",
    "    tot_new_Ts=0.\n",
    "    tot_Tu=0.\n",
    "    for i in range(1140):\n",
    "        #print(i)\n",
    "        for j in range(15):\n",
    "            Ts[j]=Tu4_orig_F[basDTid[basis],i,j]*(Ts4_orig_F[basDTid[basis],i,0]/Tu4_orig_F[basDTid[basis],i,0])\n",
    "        #print((Ts[0]-Ts[1]-Ts[2])*ha2kc)\n",
    "        #print(Ts_val_Fo10Q5[bas5id[basis],i])\n",
    "        inter=mbe.get_4body_4123(Ts)[0]\n",
    "        tot_Ts+=Ts4_val_F[basDTid[basis],i][0]\n",
    "        tot_Tu+=Tu4_val_F[basDTid[basis],i][0]\n",
    "        tot_new_Ts+=inter*ha2kc\n",
    "    print(tot_Ts, tot_Tu, tot_new_Ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table 1\n",
    "\n",
    "that is, most of it, the triples scaling part is right above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 371,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method       AVTZ       AVQZ     AVTZ->AVQZ      AV5Z      AVQZ->AV5Z\n",
      "HF     3.836182808812448 3.831864800986581  ---  3.8275206749278974  --- \n",
      "CCSD   -8.17484324532663 -8.38769524169429 -8.543019671476095 -8.449853540380532 -8.515068804575934\n",
      "(T)    -1.5279711327634202 -1.5972364370974481 -1.6477813889087658 -1.6182593279915878 -1.640316131552653\n",
      "HF+CABS     3.835462060060908 3.83049423466074  ---  3.827248685954868  --- \n",
      "CCSD-F12b   -8.4671999290372 -8.493662732726417  ---  -8.499713057375118  --- \n",
      "T_unscaled  -1.5126478713537377 -1.5893042902585113 -1.6452427581079405 -1.6136614093993906 -1.6392164196455592\n",
      "T_scaled  -1.6243795327379544 -1.637391434545468  ---  -1.6379627757684831  --- \n"
     ]
    }
   ],
   "source": [
    "#canonical CCSD(T)\n",
    "print('Method       AVTZ       AVQZ     AVTZ->AVQZ      AV5Z      AVQZ->AV5Z')\n",
    "print('HF    ', np.sum(HF2_val_TQ[0,:]), np.sum(HF2_val_TQ[1,:]), \" --- \" ,   np.sum(HF2_val_Q5[1,:]), \" --- \")\n",
    "print('CCSD  ', np.sum(CCSD2_val_TQ[0,:]), np.sum(CCSD2_val_TQ[1,:]), np.sum(CCSD2_val_TQ[2,:]), np.sum(CCSD2_val_Q5[1,:]), np.sum(CCSD2_val_Q5[2,:]))\n",
    "print('(T)   ', np.sum(T2_val_TQ[0,:]), np.sum(T2_val_TQ[1,:]), np.sum(T2_val_TQ[2,:]), np.sum(T2_val_Q5[1,:]), np.sum(T2_val_Q5[2,:]))\n",
    "\n",
    "#CCSD(T)-F12b\n",
    "print('HF+CABS    ', np.sum(CABS2_val_FTQ[0,:]), np.sum(CABS2_val_FTQ[1,:]), \" --- \" ,   np.sum(CABS2_val_FQ5[1,:]), \" --- \")\n",
    "print('CCSD-F12b  ', np.sum(CCSDF12b2_val_FTQ[0,:]), np.sum(CCSDF12b2_val_FTQ[1,:]), \" --- \" ,   np.sum(CCSDF12b2_val_FQ5[1,:]), \" --- \")\n",
    "print('T_unscaled ', np.sum(Tu2_val_FTQ[0,:]), np.sum(Tu2_val_FTQ[1,:]), np.sum(Tu2_val_FTQ[2,:]), np.sum(Tu2_val_FQ5[1,:]), np.sum(Tu2_val_FQ5[2,:]))\n",
    "print('T_scaled ', np.sum(Ts2_val_FTQ[0,:]), np.sum(Ts2_val_FTQ[1,:]), \" --- \", np.sum(Ts2_val_FQ5[1,:]), \" --- \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Total 2-body\n",
    "\n",
    "The total two-body contribution without and with the use of F12 is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 335,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total canonical CCSD(T):\n",
      "-6.32786426120069\n",
      "Total CCSD(T)-F12:\n",
      "-6.310427147188735\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print('Total canonical CCSD(T):')\n",
    "#print(np.sum(HF2_val_Q5[1,:])+np.sum(CCSD2_val_Q5[2,:])+np.sum(T2_val_Q5[2,:]))\n",
    "print(np.sum(CCSDT2_tot))\n",
    "print('Total CCSD(T)-F12:')\n",
    "#print(np.sum(CABS2_val_FQ5[1,:])+np.sum(CCSDF12b2_val_FQ5[1,:])+np.sum(Ts2_val_FQ5[1,:]))\n",
    "print(np.sum(CCSDTF12b2_tot))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 336,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HF+CABS:    3.8352517997130624\n",
      "CCSD-F12b:  -8.259842899392417\n",
      "T scaled:   -1.5436774296783413 T unscaled:  -1.427125071251579\n",
      "Total HF+CABS+CCSD-F12b+T unscaled:  -5.851716170930934\n"
     ]
    }
   ],
   "source": [
    "# \"Our value differs from $-5.85$~kcal/mol obtained by \n",
    "# Deible and co-workers in Ref.~\\citenum{deible2014theoretical}.\"\n",
    "\n",
    "#HF\n",
    "print('HF+CABS:   ', np.sum(CABS2_val_Fo8))\n",
    "#CCSD-F12b\n",
    "print('CCSD-F12b: ',np.sum(CCSDF12b2_val_Fo8))\n",
    "#T\n",
    "print('T scaled:  ', np.sum(Ts2_val_Fo8), 'T unscaled: ', np.sum(Tu2_val_Fo8))\n",
    "#total\n",
    "print('Total HF+CABS+CCSD-F12b+T unscaled: ', np.sum(CABS2_val_Fo8)+np.sum(CCSDF12b2_val_Fo8)+np.sum(Tu2_val_Fo8))\n",
    "#CABS2_val_Fo8[i]=mbe.get_2body(data[1])*ha2kc\n",
    "#CCSD2_val_Fo8[i]=mbe.get_2body(data[4])*ha2kc\n",
    "#CCSDF12b2_val_Fo8[i]=mbe.get_2body(data[5])*ha2kc\n",
    "#Ts2_val_Fo8[i]=mbe.get_2body(data[6])*ha2kc\n",
    "#Tu2_val_Fo8[i]=mbe.get_2body(data[7])*ha2kc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trimers, Table 2\n",
    "\n",
    "The data for (T) contribution with common scaling factor is in the first section of Analysis\n",
    "\n",
    "Table 2 results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 337,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method       AVTZ       AVQZ     AVTZ->AVQZ\n",
      "HF    -0.2727033005633549 -0.2728307859854403  --- \n",
      "CCSD  1.1116361822708614 1.1084722807095124 1.1061634876782573\n",
      "(T)   0.20515132647030618 0.2031772778798679 0.20173675593549398\n",
      "HF+CABS  -0.273348299683872 -0.2722324559092421  --- \n",
      "CCSD-F12  1.1055751918404624 1.105833678108463  --- \n",
      "T_unscaled  0.201295864293327 0.20154342620099752 0.20172407948497342\n",
      "T_scaled  0.04639468818115126 0.11417086640714369  --- \n"
     ]
    }
   ],
   "source": [
    "#canonical CCSD(T)\n",
    "print ('Method       AVTZ       AVQZ     AVTZ->AVQZ')\n",
    "print('HF   ', np.sum(HF3_val[0,:,0]), np.sum(HF3_val[1,:,0]), \" --- \" )\n",
    "print('CCSD ', np.sum(CCSD3_val[0,:,0]), np.sum(CCSD3_val[1,:,0]), np.sum(CCSD3_val[2,:,0]))\n",
    "print('(T)  ', np.sum(T3_val[0,:,0]), np.sum(T3_val[1,:,0]), np.sum(T3_val[2,:,0]))\n",
    "\n",
    "#CCSD(T)-F12b\n",
    "print('HF+CABS ', np.sum(CABS3_val_F[0,:,0]), np.sum(CABS3_val_F[1,:,0]), \" --- \" )\n",
    "print('CCSD-F12 ', np.sum(CCSDF12b3_val_F[0,:,0]), np.sum(CCSDF12b3_val_F[1,:,0]), \" --- \" )\n",
    "print('T_unscaled ', np.sum(Tu3_val_F[0,:,0]), np.sum(Tu3_val_F[1,:,0]), np.sum(Tu3_val_F[2,:,0]))\n",
    "print('T_scaled ', np.sum(Ts3_val_F[0,:,0]), np.sum(Ts3_val_F[1,:,0]), \" --- \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total contribution of three-body terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 338,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Three-body total:  1.0353253016841943\n"
     ]
    }
   ],
   "source": [
    "print('Three-body total: ', np.sum(CCSDTF12b3_tot_F))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tetramers, Table 3\n",
    "\n",
    "Only CCSD(T)-F12 data calculated, (T) contribution with a common scaling factor in first part of Analysis\n",
    "\n",
    "Table 3 results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 339,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method       AVDZ       AVTZ\n",
      "HF+CABS  0.5476723644914855 0.5521816614592011\n",
      "CCSD-F12  -0.010139570699090083 -0.01615748699134508\n",
      "T_unscaled  0.01920171653016218 0.020982149109879326\n",
      "T_scaled  -1.4795357199584835 0.8952191456874072\n"
     ]
    }
   ],
   "source": [
    "#CCSD(T)-F12b\n",
    "print ('Method       AVDZ       AVTZ' )\n",
    "print('HF+CABS ', np.sum(CABS4_val_F[0,:,0]), np.sum(CABS4_val_F[1,:,0]))\n",
    "print('CCSD-F12 ', np.sum(CCSDF12b4_val_F[0,:,0]), np.sum(CCSDF12b4_val_F[1,:,0]) )\n",
    "print('T_unscaled ', np.sum(Tu4_val_F[0,:,0]), np.sum(Tu4_val_F[1,:,0]) )\n",
    "print('T_scaled ', np.sum(Ts4_val_F[0,:,0]), np.sum(Ts4_val_F[1,:,0]) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moller-Plesset perturbation theory\n",
    "\n",
    "### Table 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 340,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method   2b |3b 1hb | 3b 0hb | 4b 2hb | 4b 1hb | 4b 0hb | Total  \n",
      "HF+CABS  3.827 -1.483 1.210 0.164 0.471 -0.083 4.107197891504827\n",
      "CCSD-F12  -8.500 0.913 0.193 -0.093 0.054 0.023 -7.410036866258\n",
      "(T)  -1.638 0.172 0.030 -0.015 0.027 0.009 -1.4152565471736307\n",
      "CC-F12 tot -6.310 -0.398 1.434 0.056 0.552 -0.051 -4.718095521926804\n",
      "MP 2 -5.192 -0.956 1.208 0.117 0.563 -0.067 -4.326238725786246\n",
      "MP 3 -5.452 -0.382 1.501 0.077 0.519 -0.067 -3.8038250474665607\n",
      "MP 4 -6.360 -0.433 1.409 0.000 0.000 0.000 -5.383625383417596\n"
     ]
    }
   ],
   "source": [
    "#we print the data stored in the mbe arrays\n",
    "print('Method   2b |3b 1hb | 3b 0hb | 4b 2hb | 4b 1hb | 4b 0hb | Total  ')\n",
    "print('HF+CABS ', \" \".join(\"{:.3f}\".format(CC_mbeQ5[0,ord]) for ord in range(1,7)), np.sum(CC_mbeQ5[0,:]) )\n",
    "print('CCSD-F12 ', \" \".join(\"{:.3f}\".format(CC_mbeQ5[1,ord]) for ord in range(1,7)), np.sum(CC_mbeQ5[1,:]) )\n",
    "print('(T) ', \" \".join(\"{:.3f}\".format(CC_mbeQ5[2,ord]) for ord in range(1,7)), np.sum(CC_mbeQ5[2,:]) )\n",
    "print('CC-F12 tot', \" \".join(\"{:.3f}\".format(CC_mbeQ5[3,ord]) for ord in range(1,7)), np.sum(CC_mbeQ5[3,:]) )\n",
    "#print('CC    tot',\" \".join(\"{:.3f}\".format(CCcan_mbe[3,ord]) for ord in range(1,7)) )\n",
    "\n",
    "for mpn in range(3):\n",
    "    print('MP',mpn+2, \" \".join(\"{:.3f}\".format(MP_mbe[mpn,ord]) for ord in range(1,7)), np.sum(MP_mbe[mpn,:]) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 341,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HF+CABS  0.5521816614592011\n",
      "MP2  0.6127398441951682\n",
      "MP3  0.5286505765036471\n",
      "CCSD(T)-F12  0.5570063235777354\n"
     ]
    }
   ],
   "source": [
    "# \"too repulsive total four-body contribution\"\n",
    "# Four-body contributions by different methods:\n",
    "# note that 4:7 is 4, 5, and 6\n",
    "print('HF+CABS ', np.sum(CC_mbeQ5[0,4:7]))\n",
    "print('MP2 ', np.sum(MP_mbe[0,4:7]))\n",
    "print('MP3 ', np.sum(MP_mbe[1,4:7]))\n",
    "print('CCSD(T)-F12 ', np.sum(CC_mbeQ5[3,4:7]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 342,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD4CAYAAAD2FnFTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dfXwU9bX/P99NNptsYjZIAIkKidXbVtRbMfhwFW+3iBVBBIohPsR629+NXXt/r+ulXH4CDWJCbiQ1Yq1XfO4PSVpIuWhTWn/1Iem1iNrEq1K0WjEBtGjFWrCAQpI9vz/mgdnZmdmZ3dnn83695pXs7DycmZ055/s953zPVxARGIZhGEaPJ90CMAzDMJkJGwiGYRjGEDYQDMMwjCFsIBiGYRhD2EAwDMMwhhSmWwA3qayspOrq6nSLwTAMk1W88sorHxPROP36nDIQ1dXVGBgYSLcYDMMwWYUQYo/RenYxMQzDMIawgWAYhmEMYQPBMAzDGMIGgmEYhjGEDQTDMAxjCBsIhkkh7S+0o2+oL2Jd31Af2l9oT5NEDGMOGwiGSSHTqqahbnOdaiT6hvpQt7kO06qmpVkyhokmp8ZBMEymE6wJonthN+o21yFUG8K6gXXoXtiNYE0w3aIxTBTcg2CYFBOsCSJUG0LL8y0I1YbYODAZCxsIhkkxfUN9WDewDk2XNmHdwLqomESmwvGT/IMNBMOkECXm0L2wG83BZtXdlA1GguMn+YcrBkIIcYUQ4m0hxC4hxG0G3/uEEJvk718WQlRrvlsmr39bCPF1zfoKIcRmIcRbQog/CCEuckNWhkkn/fv6I2IOSkyif19/miWLjTZ+srJvpWro2EWWu4hE56QWQhQA+COAmQDeB9AP4FoielOzzS0AziGi7wgh6gHMJ6JFQogzAfwUwPkAqgA8C+DviGhUCLEewG+J6BEhRBEAPxEdsJKltraWuFgfwySXlX0r0fJ8C5oubUJzsDnd4jAuIIR4hYhq9evd6EGcD2AXEQ0S0TEAGwFcrdvmagDr5f83A5ghhBDy+o1EdJSIhgDsAnC+EKIcwKUAHgUAIjoWyzgwDJN8sjV+wsSHGwbiZADvaT6/L68z3IaIRgAcBDDWYt/TAOwH8GMhxKtCiEeEEKVGJxdCNAohBoQQA/v373fhchiGMSKb4ydMfLhhIITBOr3fymwbs/WFAKYCWEdE5wI4DCAqtgEARPQQEdUSUe24cVHzXTAM4xLZHD9h4sONgXLvAzhV8/kUAPtMtnlfCFEIIADgE4t93wfwPhG9LK/fDBMDwTBMalh68dKodcGaIAepcxg3ehD9AM4QQtTIweR6AD26bXoAfFP+fyGAXpKi4z0A6uUspxoAZwD4HRF9COA9IcQX5X1mAHgTDMMwTMpIuAdBRCNCiH8B8GsABQAeI6I3hBDNAAaIqAdSsHmDEGIXpJ5DvbzvG0KIbkjKfwTAd4loVD70/wbQJRudQQD/lKisDMMwjH0STnPNJDjNlWEYxjnJTHNlGIZhchA2EAzDMIwhbCAYhmEYQ9hAMAzDMIawgWAYhmEMYQPBMAzDGMIGgslZeIIbhkkMNhBMzsIT3DBMYrhRi4lhMhLtBDeh2hDWDazjCW4YxgHcg2BymmBNEKHaEFqeb0GoNsTGgWEcwAaCyWl4ghuGiR82EEzOwhPcMExisIFgchae4IZhEoOruTIMw+Q5XM2VYZIAj7Vgchk2EAyTADzWgsll2EAweUMyWvvasRYr+1aqQXFOp2VyATYQjKtkssslWa19HmvB5CpsIBhXyWSXS7Ja+9qxFh0vduDuF++O+j4TDCTDOIVLbTCukunlLbSt/aZLm1wxDlpDU1FcgSVPLwEALL5occT3DJNtsIFgXMdtJewm+pHVwepgQvLpx1osvmgxAGBl30oc+PxAxhlIhnEEEeXMct555xFDtGbbGuod7I1Y1zvYS2u2rUnJ+XsHe6myvZKaepuosr0ySpZ0ocilyKP/7CZNvU2EVaCm3ibXj80wbgNggAx0KscgcpB0xgEyubxFqkZWc/0nJmcwshrZunAP4jjpasWnu/eSblLZS2EYt4BJD4JLbeQwK/tWqnGA5mBzusXJC9pfaMe0qmkRMYe+oT707+vH0ouXplEyhjHHrNQGG4gcRXH1ZGImEcMwmQXXYsojMjkOwDBM9sAGIgfhMtcMw7gBu5gYJg441sDkEuxiYhgXyeSSIgzjFmwgspRMLopnRLbJGwuu4srkA64YCCHEFUKIt4UQu4QQtxl87xNCbJK/f1kIUa35bpm8/m0hxNd1+xUIIV4VQmx1Q85cIttasNkmrx24iiuT8xgNjnCyACgA8C6A0wAUAXgdwJm6bW4B8ID8fz2ATfL/Z8rb+wDUyMcp0Oy3GMBPAGy1I0u+DZTL1JIWZmSbvLHIteth8hcksdTG+QB2EdEgER0DsBHA1bptrgawXv5/M4AZQgghr99IREeJaAjALvl4EEKcAmA2gEdckDGrsOuOybYWbLbJawWnEjP5gBsG4mQA72k+vy+vM9yGiEYAHAQwNsa+9wBYCiBsdXIhRKMQYkAIMbB///54r8ExyfSp23XHZFvNn2yT1wpOJWbyAqNuhZMFwDUAHtF8bgDwI902bwA4RfP5XUgG4j8B3KBZ/yiAbwCYA+B+ed1XkYEupmTX3Inlvsi2mj/ZJi/D5BNIoovpfQCnaj6fAmCf2TZCiEIAAQCfWOx7MYC5QojdkFxWXxNCdLogq2skO4slljsm21qw2SYvwzBwpQdRCGAQUpBZCVJP0W3zXUQGqbvl/6cgMkg9CE2QmjK4B6GQrLr/HABlkkm+V91lIkGyehAkxRT+BcCvAfxBVv5vCCGahRBz5c0eBTBWCLELUmbSbfK+bwDoBvAmgP8H4LtENJqoTKkiWT51DoAyySYX046ZJGBkNbJ1yZUYBLfuGCvcej64l8oogGeUc5dk+tSXXrw0KuYQrAlyjZ88wE52nFut/1xKO2aShJHVyNYl3wbKpRru2SQfuz1TN1r/bvYg+NnIbmDSg0i7UndzYQORXDhVNTXYVdyJJEi49VsqhkG7f+9gLzX2NPKzkUWwgWBcgf3WqSGW8k/0d3A7jqEYhkBbgEpWl1CgLcDPRhbBBoJxjWSl9jIS2TZIUitvyeoSfjayEDMDwUFqxhG5VC4jE7GT4pxpgw61wW4hBD8bOQTPKMfYRqu8gjXBqM9M4mTjTHV9Q32Yt2keRsIjKPQU4slFTwIAPxtZhNmMcmwgGNtko/JionHzd1QaCQu+tAD1Z9UDOG4YAPCzkSWwgWAYBoC7PUFuNOQGbCAYhlFZcc8buPOOcoQPngJP4H3cdvunaL11SrrFYtKEmYHgIDXD5BldXcA9K6YgfOBUgATCB07FPSumoKvLvXPk2hzkVuTytbKBYJgE6eoCqqsBj0f666aiTQYrVgBHjkSuO3JEWu8W+VQMMKev1Sj3NVsXHgfBpJrOTiK/XxpRpCxeL9HYsURCEE2eLG2jkAklKYQIR8irLEKEXT1PPg2qzPZrBY+DiCaXu4aMMW639o1a48PDwF/+IqndPXuAxsbj53GrtZnIdQTGf+pofbzkUzHAnL1WI6uRrYvTHkSmjUhlkotRa9/vj2zhO0WI6Ja40TJ58vF9Em1tJnodybgPRmR7q9oJ2X6t4FIbxjj9YTPBRcDEx+TJsZW3W8eMdt8c36d3sJdmrJ8Rd0kKN66js1Nygyn7jh0bv4Ho7JTOrXWp5VPjKxeulQ2EBU5qC+XCw5CvmLX2tcrbKUatcaNlfNURIpKel/K2cgq0BeJubdq5jlgNGbd6EWbHWbSyJ28aUrnQaGQDYUI8XcNs707mK8noQRBFtqDHjiUqKoo8vq9khMrqG6mpt0k1Dok0MKyuw275bbfuRbLuabLIBWWeDNhAGJBIbyATKpryw+6MVPnejVwuyvMyY/2MhH8zq+vQG4bytnLyt/qpvK084rxu9aaS0StLJuwBMIYNhAHxKthM6UHww+4cI+WdbJLxvFhdh/Z8/la/YUMmX3sQRJnz/mYSbCBcItOUMj/s6cOOsUnX86L0WEpWlxg+G8mOQaTC8CZCJngAMgk2EC6RiW4dfthTj13FmI7nxWhmNyPD5FZvKh29skTgRlU0bCByFH7Y00O6XCuxDI7yPDT2NBoGq/M9PpVpHoBMgQ1EDsIPe/pwGpxNxhzQRp8zsYdrRarlzbb7kyrYQOQga7atoeVrd0Z075ev3Rn1sPNL4T5OexC9g73kb/VTx/YO9XNleyV1bO+I20gku9e4ZtsatSeiPXdjT6Nrzw43cjIDNhA5iF0/OL+E7hNPcLZjeweJVYIatjRQaWsphbaGon4Xu4o33riTk8aCdlCfNm1WO47DDdhNmn7YQOQgTlqx/BK6h6JktcHZ8VVHaNHKnpj7NmxpIKwCnX3/2SRWiageRazfZc22NdSxvSPit3TSC3HaWNCOpdAGvfX3Qr+P0x4GJ1qkFzYQJmRbBoYWp35wfgndId4embLd9MemE1aBZj4+07HRVnohimHRf3Yiv93zKs+N0bPjRu+UGy/phw2EAdmaw63APYj04fR+amMOpa2lNPPxmSRWCZr5+EzCKlDDlgZbre5EexAKdhsL2h5EYXMhlbaWGsYk4n222P2ZGbCBMCAbR4FqsTJwevdHWX1jzr6E6eoFOumRaRW78nfuT+aqPQmnvYBEeoN2jZs+5qD0VhQjoX2OUhETYZIHGwgDsq2OjBFGytHIcPhKRqLKMeTCS5iuXmA8PTKtMlSU7czHZ5K/1U9zuuZE1UuyCh6nosVulMXUsb2DilcX04z1MyLGV5S2llLDloa4g+5MekmqgQBwBYC3AewCcJvB9z4Am+TvXwZQrflumbz+bQBfl9edCqAPwB8AvAHgX+3IkW89CDNy9bqMSMe1milZvTJVvjNSkmu2rVED1k29TdQ7KI1+buxpNDxHrHPbNRKzOmdF9VQ6tnfQrM5Z9i6eInsvyvnndM2hQFsgopekvZ5sIV97NEkzEAAKALwL4DQARQBeB3CmbptbADwg/18PYJP8/5ny9j4ANfJxCgBMBDBV3uYEAH/UH9NoybcYhBm50DOySyquVa80FHeRVmloffF2lLdRL8BOzyBRBZaogdHLqB2xrbijGrY0UHFLcUSPKFsUb77GRJJpIC4C8GvN52UAlum2+TWAi+T/CwF8DEDot9Vup9v/5wBmxpIl37KYzOAeRPzXqn0eKiYcoOVrd0YoCSNDoMWOkl++did5KvaSEGF1cKPel9+wpSHquG4PTnPqooqlPHsHe6lkdYlaJFDvLotH8aarllW+JXQk00AsBPCI5nMDgPt02+wEcIrm87sAKgHcB+AGzfpHASzU7VsNYC+AcpPzNwIYADAwadKkZN7DrCFXe0ZGuHmthrPDeQ+rRkJfAI/IWIFpXUdG5/D6jkXJu3ztTtXwNGxpiGuMhJZYijWeoLJyTO2xtZ8bexrV8uL+Vr+peyyeybmS0aI3ukeNPY3U2NMY5UbLtJ6O2yTTQFxjYCB+pNvmDQMDMRbAfxoYiG9oPpcBeAXAAjuy5NtAOStysWdkhlvXatYb8VTspabeJrV1rFWqeoWlHS1tpMjMzjG+6ojpcbQlOrTntVJaVoo13hay1Wx15W3lVNpaSuVt5YYz5ynEY5iS1aI3ukeBtoCt68g1ss7FBMArf15sVxY2EEwimMUzIMJqq9hISSmKxk7L3+ocZj0Ro+wgJ2MvjOIc8bTI9YZB26Oa3TU7IuagxCS0QepEFH2yBnnqZerY3kHFdd8iUbGHIMIkKvZQcd23MtIV5ibJNBCFAAblILMSpJ6i2+a7uiB1t/z/FF2QelAOUgsAjwO4x4ksbCCYRDBr3YuKPRHTdhopVbuxAzsxk85OqUcBjFJg/F+prL4xanCcXhmZZSd94YdfiFCsbgW5tT2qecu7qWLCATWmovTgtMd1wzAlIybQ2UkUGP9X9V5fOH+AfCUjEb+Nr2QkZhmVbA9uJ81ASMfGlXKm0bsAVsjrmgHMlf8vBvAzSOmsvwNwmmbfFfJ+bwOYJa+7BAAB2AHgNXm5MpYcbCCYRDCLQQQXPxxlGJav3SkrxeMDEe0osFgxk85OMlRQZfWNlrENoxIcWAUqbS11RbFq608pChWB3VT1tScI3sMxY0CJTu+bDMVrdK8hwnEnPWRzcDupBiJTFjYQTKIsWtlD46uORGUxaVvCi1b2mA5EtKPAjGImigI162GMmXAwprJXjIRS60lbFiNRxdo72Etl9Y1RAXYgfoVqh2S6bqRemvH9jk6bDts6ZrbWO2MDwTA2sNNijeUmikeBKecRJi1YYNSWsleMQ/U91a4r1jETDjpQqHGfxhXsGBaz3oLRUjHhQMxzcg8iwxc2EIwbxHrRkzU4r3ewlzwVew2PPb7qSNS2emWv7UE4re1kBycKNd1jbhIx9EZLKJT4+ZKBWz0sNhBMzmCW1uqmO8LKVeB0cJ4TueYt7yZ4D8X06etxowx4LMxcMnqD6caYGzd+y1iG3igeZGb8Y/UO7c7u6DZuGSY2EExOYBXkdetliUexWClFZZyAVQqo9rxn3txKomJPRFaQXjHpFeiszlkU2hqK2MZpjSUrlBiEUQA9FHJ/zI3Zb6lXwrHOFSsmoG9sWLnMrJ6ndA5OdcO1xQYiy8ingW5OsOP/d/Ky6O+ztuyF9nhGRsLu76OMGbAafKUfY6Cf6tNs7EWqXBqJzKIXL/rfcvnanY6UcDyK0+z5Coz/q+Ux0lneRp+qG4+uYAORRdhtjeSjEbHj/7ebSWJ0n72+Y7R87c6I7dzImukd7LUsQ6HvEShGQltW2+iYmRQUVZ5HiDCNrzriSnl57W9pNQpdf+x4DajxM3GUsOBay+cpXQUyzdKineoCNhBZhN3BVE5aU7liTNzsQaSy1dc7GF3IrrOTqHjshwTdADPFNWTH0CUjrTIe/7/VHCRuufrsZHglcg3aa9FPthXreUpXD8IsLqRPaogFG4gswk5rxMkDmUvF+9yMQSS71aetXRRRCnt1MRXXfSuq5ef3E113+1YSqwSFtoZiGrp4ehB2FGc8rW/TUeieEQJGo3oUsTCSwW6Gl1s4uQ+hUHKC9bEwzSyzOW5DPQ4biOzBjvJ3otxyrfy3VRaTk0ySZN8XRaEYTaYjKvYYv9iBIdU46BWTdkIiZZ0yL4XdVrpdpaesn7e8O6I0uZnCM60xFafCNDJky9fujBqo5/UdTZoS1g6aNEsWIDLPhoqVGusGbj3DbCCSQLLcNnZa/E4ejHyZQCget1uye1ZK7EFfcM/KXWKmmLTzUGjnuNYqezsuFDs9j85OopLyQ6QfKW12f+yOKUjU+BrVqUpG7MXs2TDK2Ep3gNqNZ5gNhMskW7nEMj5Ozp9rPQgz4rnOeIy8U/+2UZzAVKGWfGQ4X4S2AJ4bgWmr2IVhjaIY97Ozk8jrjW0g4mmUGMUEkp29ZeoyM3AjuXmt8eBGQ5UNhMtkgtK1+2DkUgzCilT1lJz4ps0Uutlv4g8cNn2uFCOUaGA6lpGJVaNIiOhnLxQiKiqKbSCcvB+dnURjxxocR4Qj3DduZJnpseMyU5aCgvTrgkRhA+EinZ3WL08mkitZTFa4YbTt3ic7LflYhsQoi8ncyIWpsr2SQltDpvNS2MGOcYtVUmPsWIOWs40yHIphsYNhZd04jxUPTspwAOGsb4CxgXCJWA9uNrUaco1Ee0pO94/Vko/lijL63qz17qnYS6GtoZgTEjk9p9anr2QaWSlHv9+kVW9zMUMvlx0FbaeAntU59PdGy6KVPQZuNmMjqL1v2doAYwPhErFenmx7MJyS6S9CIvKZTjlaMGoaMLbbkjdSTh3bOyIGzfUOGpezgPcQffGf76DS1lLTqUetpgM1GwlupOiV0hlGMYixY617ObEWq8aT3tiZB/GPL3ZLcJudQ/vZaER9WX2jmiwwvuoI+S58OCo+ZDQoLRtnl2MD4RJWL0emKUu3cdLCTtVLYjaTWjw1iOwoPl/JCBXXfSuiVIZRGqqCcs1m9Zhmd82OMjTHlVWYPBV76Yv/fIfhbHX68+hLdWhnwdMSqxdcPPZDKrzgQfIUjBIQJuEZobnX76U129ZQZyfJ650ZBzuNJ63bzmzMg12DY+cc2vtt9FwrZVf028YqaxHPOJJ0wwbCJTIhOJ0unFx7ql4SN6uY2vU7l437WFX82utr7GmkQFtALcLXO9hLpa2ldN3tW1U3jqjYQxffej8VtxSrRsbIVaUd51DZHj3ftRFa5aeU9XBSjfb4MhpVURbewzT3+r2W2U1mi8cj/VUmYFJQUnW1jQZl1rx5y7stj1lUFH+DTH+/rZ5r/bZ2s8jcyjZLFWwgXCJfMoKMcJollKqXxK15EGK1rPXXa2QElRHTyjVfd/vWqCk54T1EBQtvjHAJ6e+Rojy1bo4xEw6S95qbLCuaKgqtZHWJ4/ks1EUMm3xn5vaxP08ECj6nEyqOqteD2vtozISDxz8vuE4dM3JCxVHT4xR4R+N654zud6zEAGVbo3EnVs91Ns0uxwbCRTLdD58s4uk9peolUWZSm/7YdNNt7Li9TFMrNYu2tIORwrFTYA4I05gJBy1z+o0CpYXeUfIURCpwbamRQFuASlaXWFaCtbw+b/TguOQu+nNJKaym8RiL3yHe+a2tEgO02/pb/aYxILNzcQ8ig5ZcKbWRqTjtPSXrJdEr+o7tHeo0m2KVoNldsw0NgRKwVWIFer+94mO37kWESfg/jnCVKAZhxvoZaitTueZY6Z8FRZ9FBL4bexpVOczy682UpXJty9fujMhM0pYbMR/QFibh3y/VgjIrA5LixVMwSlXn7LQwWOGI5yyeRIHeQeM5xhOp6ssxiAxd2EAkH6fjBJwEcu2iPbZiHEpbSw0/K9srrT9tK7u0tVSVRzme3ThEoe/ziAKBDVsayNvsjYqHIDAUWxF6Rqlg4Y1qq9+uqytyGaWO7R0xjbjVCO6pD0zVuMUOOTx/spawpYEApEC+Xpk7xU2vQC5lMQnpu9ygtraWBgYG0i0GA6D9hXZMq5qGYE1Q/R8ANu7ciC1vbcGyS5ZhJDyCaVXTULe5Dt0LuxGsCdo+LgDUba5DwBfAnz79E278+xvx4FUPAgDufvFubNq5CYMHBhGqDWHdwDpccPIF+O3e3+LJRU+ib3cfWp5vAQBMPWkq/vibWvifvwf7PyiBk9dhzIRPUfC9L+DCky/EM4PP4CsnfQWvPv1lDD9zB+jgqUD5Xoz/Sj8+2j4LGC6NcTTpxJMnCxw6BPzlL/blAABRsRfl/+cc+H70AT7aVxL1/eTJwO7dgBBmRwijeHUpLqu5DATCcz+fAN9vOnDwowAA050iGDvWudxu4fcDDz0EXH99es6f7QghXiGi2qgvjKxGti7cg8hM9L0JJaisL2Dn9FiKa8doAh6iyPiH0nPwt/rJ3+qnopYiwioQFlwrTQgTV+t2VCrfXfct2S0zGtXaLSj6jOZev5eEx3lqqJMW9dR/+FiWwbi1rQTWlYyi6O/D1NjTSMWriwmroPaCZnzvUULB5zHlUnopzkYgu7vkQyZhsoBJD6Iw9baKyTeCNUF0L+xG3eY6tUV/wzk3YMOODWi6tMlWz0F/rPmb5uPY6DH4vX4UeiIf4/YX2lHoKcS6gXVourQJ6wbW4dUPXsWXK7+Ml/70EgDAV+CD1+PF8HP/geGjRXFdV2D8p9jQOQLvLx8BmRxj9FgxXuqrBIXttcJjY3ycV18cC6Kxpnt5PED97b9AODzH8BhEAvVn1ePx1x/HhSdfiLZtbXjtw9fw3NAoBBpg1bESnjAeesijtt4bG4EjRxxckkvs2UPQXltXF7BiBbB3LzBpEtDayj0Mp7CBYFJCsCaIUG0ILc+3oOGcBjy16ylVeQerg46MBAAcHj6MkfCIZGCqg6jbXKe6rQo9hVjy9BLcdfldWHzRYlQUV2DJ00vgER7JKISHcXT0KPxeP4Y/nRzX9XgKh/HZP/4bSn9zLw7HMDAf7fNBciG5ZSSioRiusdFRYFPLHKsjoH9fP351/a9Qt7kOX678Mjbs2IDS3+7H4VFvjJMLVfEqf1esAPbsJfjKDuHo38oQ+9rduT9dXZIMXV2RhmrPHumzVkbGBkbdimxd2MWUuWiDubHqCcWisaeRiq65iURAcuuMmXCQrrt9qxqYbuxppDldcyKyiUJbQ+T5xg1y0HiUEBgizzduoJLKP8fn0ij5SApCw4bryHRcgZ3FKkgb7/HMr6ezk6johIOa88Y+d0GBeXC3d7DX3j1yaVHcTPk8oDUewFlMTLrQGoF4JrrRj1guuuYmKvTp/OLeQxRc/HB0jGLBtXJpBCNlN6pZHCpUEZbKTsSMLUTHJZwq9OIT/uayItXLPEoXLfgfWrSyhyCcjpSOPaFQrPLhbi9E+TNJlluwgWBShj5lcNHKnoTS/vQGpmzcx4Yvf+m4/RHpr95rbnIQgJYNRclH0iLCUtC3ZL/h9mPHUrSRilCaYQNFHM/idg9CudbIYxYVEflKEunpHF+0M60JoQzMc3IN8V9zQYH0zHAPwhlsIJiUYJbDr1QCtYNV5dOm3iYLxTtKTb1NapbUmAkHnSmYwBBhFWhO1xwikgxTSd23qaDos8jtvIep5ASrcQLJUOpuLi6UzIixRD8D9o8txGhCvQ6z5zBfSuLEAxsIJiW4UQ7dbJCdUshNchkZnCMwRL4WHxWvLpaMhOOy1KNUcEeBOlht8mSSRkIrvQqMEgK7pVhGCv3q7hsCt7ZPzlLsH5ZGfOtKa9tZtD2EfC2JEw9sIJiUYEcp23lZtWU6yuob5d6AVGLZc/4D5CnStzBll05giLwLv0nL1+50VKoCIMkIKEFsnbL0+o5RwUIpwI5VyJhyFLm5jEb97nYC/dxDiJ+kGggAVwB4G8AuALcZfO8DsEn+/mUA1Zrvlsnr3wbwdbvHNFrYQKQfpwOlrIyFEmQuiPL1hwmVrxE8w2TszjlKnkKnrU87bqGwZEAWXCudP0Na3Lm2FJYeiCjwqDwHVr02uz0E7lUYkzQDAaAAwLsATgNQBOB1AGfqtrkFwGb9u9cAACAASURBVAPy//UANsn/nylv7wNQIx+nwM4xjRY2EOknHtdAoXeUSkuPfz6h4igFFz9Mle2V5u6ktCrnTI8xZNNicB8LPqeLb70/aj4M1N5HeiPhKxmhRSt7bD+bHJcwxsxAeFwYSnE+gF1ENEhExwBsBHC1bpurAayX/98MYIYQQsjrNxLRUSIaknsL59s8JpOBXH898O93/hEo+QtgOf72OCPDHhw+fPzz3w4Uoe+eBjR6f4OD+wMmeyVv0FlsRJrPn03onoGCo/CXH4EQhJLKj4Civ0XvMurDi/93Dq458xoseXoJrjnzGvS83QPv3H8DFtwABHYDCAOBPRBXNeLmm8oAAH1DfWh/od1UkhUrokd4HzkirWeMccNAnAzgPc3n9+V1htsQ0QiAgwDGWuxr55gAACFEoxBiQAgxsH///gQug3GL1lunoPeNHSir/w4C4w/ArqGIIOzDmiVT4to1f8j0m0M4b+7LGF/1GYAwysb9Bdd9/1l4l1Wh8Bv/hM9HPgeOnWC4Z/jAKXh0/We47LTL8OPXfozh8DCKCopQPPUJeL/3dyi8wwf8WzVwzk8ASMahbnOdWsjRiL17na0HpBHZ1dVSqZLqaulzPuGGgTBqSumfXLNtnK6PXkn0EBHVElHtuHHjLAVlUkewJoh/a5yAg7eMwbzlmy2qiJozOgpwSz17KTnhCPZcchU2bnsJvYP/jeJ//xKODB/BwZZdGN78Y9CBSTD/fQWo50E882QlLj71Ypw17iwEq4P41XW/wuzPuzBy9zvAqlF8/oO38C9rXrRVEXjSJGfrlXIde/ZIDimlXIcdI5ErhsUNA/E+gFM1n08BsM9sGyFEIYAAgE8s9rVzTCaD6RvqU4vlPTv4LIqKR9ItEpNCfCWjeHhdqVqkEQAuOXA/nvzBlcCRStgx/MNHixDYdj9e//PrWHTWIrz0p5fwyPojePIHs4GD1QA8wMHJePOxf8UlB+6PWc+rtVUqC67F75fWGxGvSyoRw5JxGAUmnCyQCv4NQgoyKwHlKbptvovIIHW3/P8URAapByEFqGMe02jhIHV60GeGLF+7Ux3H0NlJcU10z0usJUyofD3BOk9JkgvhiIGRSn0sBHY7Pp4Qx1OeQ1tDphMwaacHNXomFVmcZDHFW64jG0dxI8lprlcC+COkzKMV8rpmAHPl/4sB/AxSEPp3AE7T7LtC3u9tALOsjhlrYQORWszmb9ZO1xjX/AAiGwahpX/xlYxQKETyqO7My6zSzpVd3lZO8QwurJhwgIikY8zqnGU6hasQ4YhGiVW2kl0jEa+iz8Y6UEk1EJmysIFIjFgvzqKVPVQ+xl5tI+UlijUns9FLF9+Um/m5eAqkGlKl4/YTau8jj74siOPFXSMzebJUfTfQFrAsn+EpPEbwRJc00U4lunztTtMe0+TJx+t7WSl2J6mu8abFcg8iQxc2EPFjp9VV6LXfAlRaS3Zq6hi9dFpjlWkt40xdCn2fk+f8BxK4X05rNNk4T4yWPSD3gBZcSxffej95KvYSRPQ808vX7iR4D9t6fqxa8E6VdzwD67JxvAUbCMaSWC+OU1fR2LGafWL0ImK9OKkuF528JaxbzLaJ/xxl4z524OvXyFLyURwyHd/fbCpTBIbUEdFEkqLXGoHrbt8aMW+HUm9Luw8RUfHYD42P7xmJen6snuVUuX+ybcQ2GwjGklgvjtPCd16vPSWjlGe2orOToiuqZvUySvCaVJpVigImcOzrbt9q2tqOWMQwKQUIr7t9q2lpcztLURFRgb6H6T1M85Z3qz0Is7nJlcmjrOYqFxaxB6PnxawFn43un1RgZiDcSHNlcoBYOeJm3xvh8QDDw/q1xmmN0lgHa6r+oQ9F8/4FCOwBMn5wmB08wPAJiL4WAqZsirEvAQhDeIzvQ4k/jBf+72xguAQQI9L2/o9R4A1HH4cKoaSKdv/HLIijJ8ZzMQCAY8eAUe8nGF/1GYQgiIq9mHjtSoy/8Fk11XXjzo1Ydsky9O/rBwCMhEdw1+V3YSQ8gr6hPrRta8Ndl9+Fs8afpe5z8y9uRt9QHyZNMn5+jNZffz3w0EPA5MmAENLfhx6S1jtNdc17jKxGti7cg4gfOzEIO+4Pp8FlO77fsnEfk/eamyjQFqAzb24leJyXgc6axXvQsiXvDxymxp5GKim30UMAEbyHaN7y7oj7KQW2kyH/qJrS6m/1U2lrKQXaAmrvobGn0XR6Wf0cIJ2dsmtRdkXNvX6va379bHP/pAKwi4mJRawXx8oXrd3HrBuvd1MpL7heOSxa2RNd8E+eUnT52p0Ej9lMbu4puvQFxsNUUPqJpWxSCWz78unHCNjNLLO+PwbrA7vJ3+qnktUlEYZBG2OwM/e4YTDbe5jOuOJpaRyEQRCbSQw2EEzCxPLfao2D3hh4fcdo7vV7owbUKcZBqzzMZoJLrOVrpRT136W7Umwcc2RbLNoxAkTS2IK4ZSv5iDznryN4dTPqeQ/R+IZbCasQFWRu6m0yDDw7fc4QGKKGLQ1EFD2pFJMYZgaCYxCMbVpbgaKiyHVFRdJ6bXkBQHqllfpLkydLFV63//1U/LivD+Ew8OO+Pjw0/FVMq5qGYE1Q9Tmv7FuJv/65zPD84dF4H1fCjBkCoiAqMAIYlv6Kt/4ToaD4CFDysXzc+I5RNu6vpt9Zy2Z8zkmTBLoXdqu+//s6AvD69PfCRF4xCqVyKhZcj7Lvn4bwlSHgqn+OqKpacPUt+OgL98BX4ENJYQnuffle9A31RZRcWTewDn1DfRbyS5gWzzs4CVv+sAV9Q33qM6NcE5MkjKxGti7cg0gunZ3R2UmegmFavnanrd6FlK4qzQpXVt8Y1fpTWprmc0DEtxQVSW4rf8Cm3z7O1vXFt96vtpQXrexxNG7keEv8UwouftixG2h81RHyXfgweX2RAxnN/PTadFOpfIXZ+aRpWJWeAVaBCpsLIz4XtRRRaWsplbaW0pyuOdTY00jlbeURMQii463+xp7o314Z5EZk3oPwVOxV54ew6jkYzWmuPT4TDbgHwSTKihXR2Unh0ULceUc59u41boHu3Xu8d/HRvhIAHhz8qAJHn/gR9m2Xiqv1DfXh5l/crLY0R7/2fwxauPFz7Bjwy3X/gCMHi107ph5R8R52TlimtpTPv+IdFH8jBCEc9iSGS1F63s8ddUDGTPgU4Vsn4amffAE/frQIYyZ8CiCM8VWfqdk7elpvnYIf/OpnELcXYPraG+UMMQMC72GURuHRqIqR8AiqyqrUdeFwGF+t/iqag8146U8vof6sejy56ElMGTcFi6YsUovoKa1+AKjbXKf2JvSluo0yjeA9gttu/xSLL1ocs+cwrWqa5fEZBxhZjWxduAeRXMzHSkhBQ7MehFmLcHzVEeod7KVAW4DK28oj8uOx4DoaM+FgjKwbB61sMWoxj7RTv370fNUldd+mxp5GWrNtjdrK7djeQU7HNJRU/pkq2ytN4zD6c/tKRsh7zU3Usb0jIlOoY3uHZYtZac0rA9MuvvX+6LETXmmEs9Jr8LX4InoO/la/ur+vxRfxG8YiVvBajWfJPZx5y7sdPavxBMfzGXCQmkkUKzfSvOXdUYFLxb1hPshulAJtgSiXg6JkFQVnWGbBe0geVGZT+QZ2y0pQF1x1FLw+Lvf4qiMkhGTkFAVNJCmm0tZSCm0NOa5g6isZUYvbldU3RtVV8vuJQqHI9M/la3eqxfD8rX5LJa11vSj3ONAWoBnrZ1BleyXN+N6jsrxybSfZOCiLt9lLpa2lqsuptLWUegd7VSMx9YGpjp6nWMFrKyVvx43kNDiez7CBYBLGbKyEUt573vJu8lTsJSHCEWmyVlkpM9bPiHneNdvWqHEORSmX1H1bGhOhMxxFRdFxkkLf51Sw8EZJ0S24ljwVe6SWvUWpbOE9bGqAxlcdiZBPr5iUEcFn/PAMwoLrSBTFLhUiNCUj1mxbQ2fcewZhwbXkG/sBKamtwcUPR5xHq0D9rf6YylCb+aPvufUO9lLx6mLyt/qpY3uHGkPwNntp0t2TyLPKQ1gFmvn4TGrsaVRjDXO65lBpayk1bGlQDYad31PpZSnKX9/j0WcpxfuZexD2YAPBuILV3A9ExumHZnnt2jIMdtEfP7j4YRIVeyKM0nW3b5XcSXIre8b3Ho0KsHru8Fi4f8JS63nBdbaDvnpmPj6TsAp09v1nk/eam6QaSRil8jFHowyY13c0Imgf2hpSW+zlbeXqZyMFrLSSS1aXmCpDbWtbuX9TH5hKvhZfxLaNPY2qm0zp1Sluq/K2cjr9h6cTVkFV5qGtoahSGf5Wf8zf06zEhvJZL7P2tzczkE6MBxMNGwjGNk5GmtrNGDl+zMhBTk5fXq3yIpKUmpI9o7RMxSpBc7rmENFx5XPBwxdIRkE2FKWtpVR04j7Tno2y3Re+vVLNviob97E6dsPqnilup7PvP5uwChG5+2u2rYnI6BpfdUR1KyktaX+rn0JbQ1TeVq76/QubC9VrUujY3kHFq4sNB6Zp76d+neISUuSyQp95pNxPxQBe+PCFEdvFin0ov2GsHoRdjNxInMXkHDYQjC3M3EiLVva48tK58fLqXSWlraVqkbfC5kIqbimO8LWHtobo/IfOV33nnlUemt01W45tRA/4woJrae5P5lJxSzFhFdR4wuyu2VEtU0Xh6++Z13eMyuobqWFLQ1Tr2Oo+zOqcpa5XlJ/i3zcqdDena06UYTC6n9qgtFlBPKvfS3t8fWA6XjdOojECdiO5BxsIxhZWGUd65aj4q7WkqqWmVQ7alnZRc1FU7n2gLUClraVU3lauKkhfi48CbQE67dvfV7OlCireIyy4lmruqSFvs5d6B3sptDVE3mZvlKLXKkyre0YU6XqxYyCVFrYSeC5ZXaL2krQxB6f3Xt9zMOptWMmn3MuS1SVqQNysPLeT3y8e5c5uJHdhA8HYwqrst/alVjJ19K1au4FKN1BaoA1bGtQgraI4tcpnTteciOwebdaN1jVS2V6puoVmPj7T8DxGSs3OHANaJRtLsXVs74iIOWg/21XIyrm0PYDi1cVq7EFxV2m3iSVf72AvlawuUc9vVZ7bCjeUO7uR3MXMQPBAOSYCq7LfwZogQrUhtDzfggVfXoCfvfkzLLtkGeo21+HGJ27EkqeXoDnYrA6MSiZKCYeGcxrQuaMTAgJNlzah0FOI5v9uxqzTZ6Hl+RaEakOYPnk66qfUq/s9tespNJzTgDc/flMtQb3k6SU496RzsfOjnZh52kw8O/gs7n7x7ohSEU/teiriuMp1xiqVDgD9+/qjyorc+MSNmPPTOVjwpQUR+42ER3DhyRdiNDyKvt1SGeyOyzsQrA7iv/7wX7bKVkyrmoY5P52D3Qd2Y/6m+bjqp1ehqKAIF5xyAcIUxi/f+SWu+slV6rbagWQLvrRALXtSt7kO3Qu7EawJYuPOjSgqKELTpU249+V7ccd/3xFVnttOKY3+ff3qMQHEVTZj6cVLo56zYE0QSy9eavsYjA2MrEa2LtyDSByrst96t4DS6p7+2HTbQU830LY4lSC1Nl3T3+qn4tXFES393kFprEDx6mKa0zVHdeEE2gI0u2s2nXzXyRHX0LG9g4pbisnX4ovKtpn5+MyIbB2je6aMaVDo2N5B3maveiylJ3D6D09XXUf6npi2txBPq1uR9/Qfnq7KrbjJYgW49fEB/fmUeab1wXBuwWcnYBcTYxejLCYzBaVks0x/bHrKfMD6AV/awKziJ2/saYyQUzEm/lY/+Vp8EYrS2+xV3Up6JagEphVjGNoaotLW0oiaQEbjNIrrvhU1slpJCz3vwfPU84lVIsItpsQ6FPedYuRi1S8yujdExw1RoC0QZcRnrJ8RlSLbsb1DHY2tNbDs0slt2EAwCWGkIJT8fMUHbaeQWrJl1CtSxTAoikxRir4WH/lb/VTUXERYhQifvFGaqNKiN8piMjKeyiA0fXD75A6pp1JzT416z0pbS8nX4lN7YnpDZddHrzVGynUrxlAZk+Fr8dHsrtkRAWfFQDZsaYgqe8LB3/zAzEBwDIKxhd7n2zfUh/Wvr0fH5R14fP7j6F7YjbZtbRFTSqaaaVXTsOWtLREy1m2uQ/1Z9Vh68VL0DfWhf18//F4/iAhHho/gWPgYfAU+LL5oMQBjf3iwJogLT7kQG3ZswKzTZ0Xch3c/eVf1qWv99k8segL/esG/YsOODbjstMvQtq0NNz5xI/70tz9BQGDowBBmnT4L5550Lgo9hSjzluG3e3+L6ZOm44X3XsCyS5YBANpfaLfto59WNU39DeZvmo/HXnsMR0ePgkCYedpM/P6j34NAUvzhp1fhiUVP4Kav3ITh8DAKRAE27NiA6ZOm48lFTyYUH2ByByEZj9ygtraWBgYG0i1GXtD+QrsadFVQFHA6A4WKUQjVhrBuYJ0aDFXWdy/sxqsfvorvPf09dZ9Sbyl+ce0vDIPr7S+0q4Hvz4Y/w3B4GN+p/Q5Gw6PY+MZGCAg8segJBGuCWNm3Ei3Pt6Dp0iYEq4Oo21yHWafPQueOTkydOBWvfPAKzpt4Hl754BV4PV4Mh4dR6i3FjX9/I9YNrIOvwAdfoQ+3/+PtuOO/74g4ttPrLysqw+4Du+GBByXeEiy+aDHufvFujIZH8ZWTvoLXPnwN10y5Bp07OnHZaZfhmcFnMPWkqdj76d6IADKTHwghXiGi2qj1bCCYXOOyxy/Dc0PPoenSJjQHmwEAN//iZgBA/Vn1mLdpHo6OHMXR0aOYetJU7PrrLlNlfPeLd2PJ00tw1+V34dyTzsWVXVfi89HP4fV44ff61X20hunel+8FgdSW+C2/vAXrBtZhYtlEfHDoA4RqQ6j0V6Ll+RYICBR6CnHnZXfi3JPOxfxN8/HZyGcAgJv+/iY8eNWDhtdoZaB3frQTG3ZsQHVFNT757BNc/cWrsWHHBtVw9e/rx6Fjh9DyfAtmnjYTr374qmpQl12yDCPhEc4GyjPMDAS7mJispf2F9qi0yrtfvBvb9m6D3+vHD1/+oTqr2Za3tqD+rHr07+vHpZMuRXFhMZoubcLeT/fi0kmXYmLZRGzcuVE9jjJHxbODz+Kuy+9C27Y29O3ug5CnyRsOD6O2qhYbd27E3S/efTwdtDqIiWUTcWz0mHqcn735M8w8bSY+OPQBzpt4Hrp+34V7X74XTZc2ocBTgFEaxbknnYtgTRBzvzgXx0aP4axxZ5kaB8B8zoPdB3ajc0cnGs5pwKFjh3D92dern9cNrFP3VVKEnx18FssuWYbmYLPqJpxWNc3w3vYN9aH9hfbEfzgma+AeBJO1aN1GwZpgVGt//qb5ODZ6DEUFRVEtfa3raf6m+RgJj6DAU4AnFz0JAJi/aX5EL0BxH/kKfCjwFEgBPOEBgRAOh9E6oxXnnnQu5m2aBwGBlf+4EiPhEQBAoacQbdvaMOv0WdiwYwOKC4tx4zk3qr2Z0fAoCjwFuPqLV6NzRyduOOcGPLXrKSz40gLUn1Vv6sbTu9OWXbIMTX1NaAm2qC6lJU8vwXdqv4PqimpMq5oWcV39+/pV2bT3Qxmzob9P2s+5RKa6S1OJWQ8i7ZlHbi6cxZS5JCtNUj+6W1t+Qsnl15YUN5NDyevXlpHQZvIo5TyUjCRlXWlrKflb/WpJDKOxAdosoNldsyMyv5QU3akPTI0qg6GM07DKKNKOV9Bfm35eDaLjFVv112/0OyRaDiNb4LIdnObKpJlkvoRGRd/iUW7a4nj6AWKKYtWmkSqGRUmdNSqBoSjtWZ2zogbJzXx8plqcT5mdTX+PjMYkOL3GRIxzvky6ky/G0Aw2EEza0FcEVQrs6Vuy8WD0YsdjjLTjArQ9iFitcmWMhLfZa9iDUFBGNSsD4JQBhnN/MjemrFYG0M41xmuc801p5osxNIINBJMw8bZEjUo42JlYJhZmis/uqGPtd8oAMaWXoJ9tTWsQjEqNa91OZkZCGVh49v1nq8ZC7xbTy6q4vrQVXJUeiZFsVtdopez1v63ViHQjsn2kdb4ZQz1JMRAATgTwDIB35L9jTLb7przNOwC+qVl/HoDfA9gF4F4cD5r/AMBbAHYAeAJAhR152EAkl0TcRGYt9ERwSynpZ1DTlpxQXDxGrp0Z62dE1GoiInUaTv3MZ2u2raE129ao1WKnPzY9prxaw9U7eHyyHsW4OnXTWbWQ9cdSZpHTX3csA+T02cgEw8IxiOQZiHYAt8n/3wZgjcE2JwIYlP+Okf8fI3/3OwAXARAAngIwS15/OYBC+f81Rsc1WthAJJ94W1q9g5GlojP5JbRzjWbK1krZKAp++mPTDScR0mPkmjMqROjWNSXaio5n/0xQzplgpNJNsgzE2wAmyv9PBPC2wTbXAnhQ8/lBed1EAG+ZbadZPx9Alx152ECkhnh8tVpXidZFk8hLmMwX205r20wRGn1vZx5mu/LEuv/xuIuUfbTHjudexvNs5Lt7JxNIloE4oPv8V4NtlgD4vuZzk7yuFsCzmvXTAWw12P8XAG6wkKERwACAgUmTJiXtBjISmdRKdHpcuwbF6hrtnlOvKLVZTAod2ztoVucs29epBPf1hjbWfbHjLtLGT5p6m6KK9tkhEUWfzwHiTCBuAwHgWQA7DZarbRqIfzcwEN8DMM3AQPxCt+8KOQYhYslJ3INIOpnoZ3ailOzIH2sbO9fiZotYe359TMLq/juVQXtsrSFy6sJy8mwYudC0PR0mdWSdiwlSYPtFAH678rCBSC6Z6qt10vqMpThjXWOs7+0qSrv30mruC7N9FKzui1n6rjKmw6mLKZ5nQ9trUa7LqWFi3CFZBuIHuiB1u8E2JwIYkgPUY+T/T5S/6wdwoSZIfaW8/goAbwIY50QeNhD5Rzyt9UTcGXp3jaLUlFavUzdWsoKzduMk2vNrexCpigUo99PN2BTjnGQZiLEAnoOUvvqcRvHXAnhEs923IKWy7gLwT5r1tbK76l0A9+F4musuAO8BeE1eHrAjDxuIxIm3l5CO3kU8SjZR90+irhg3ZYl13Fj3Re/aSddEQRx/SD9JMRCZtrCBSJx4W7bJbhEb4dQouSVj76BUGsNqwJ9d2ZKhHJ3cF229qnS4D5NlJBlnsIFgbBPvS5vpL7tbvZzeweNjOkpWl1gGie0ExNN1vzLl/OnotTCRsIFgHBFvyzbX3QXaQGqsdFArBZxu5Zju8xNlbtJDPsIGgrFNrvYg3EAZ8KcP7pqlZpoZzHQrx3Sfn8ks2EAwtsimGEQ6cKJY88FgMrmBmYHgKUeZCPr39UfMGhasCaJ7YTf69/UnZb9ESMe0mEsvXho1o1qwJoilFy+NkEeZgW3ZJctQVlSG7oXdEVOExoKn/GQyAiOrka0L9yDyi0zrtWjPrww808tn14WTadfG5DYw6UHwnNRMVqOflzndcya7KU+mXRuTu5jNSc0uJiarCdYEEaoNoeX5FoRqQylToGYuoP59/VHyxOsuSte1MYwCGwgmq+kb6sO6gXVourQJ6wbW2fbxJ8q0qmkRMQWltV/oKYySx2zbaVXTLM+RrmtjGBUjv1O2LhyDyC/S7afXZykZxRy09YWcVlflGASTKsBZTEyukY7MKS16F9BIeMRUHqfuonRfG8MA4CA1w8SLkyAyB5yZTIaD1AzjIorC717YjeZgs+U4ByfbMkwmwQaCYeLAiQuI3UVMtsIuJoZhmDyHXUwMwzCMI9hAMLbh+kAMk1+wgWBsE++AL4ZhspPCdAvAZA9KcJXTNRkmP+AeBOMIrg/EMPkDGwjGEVwfiGHyBzYQjG14wBfD5BdsIBjb8IAvhskveKAcA0BKYZ1WNS0ipqDMb7D04qVplIxhmGTDA+UYSziFlWEYPZzmygDgFFaGYaLhHgSjwimsDMNoYQPBqHAKK8MwWthAMAA4hZVhmGjYQDAAOIWVYZhoOM2VYRgmz0lKmqsQ4kQhxDNCiHfkv2NMtvumvM07QohvatafJ4T4vRBilxDiXiGE0O23RAhBQojKRORkGIZhnJOoi+k2AM8R0RkAnpM/RyCEOBHA7QAuAHA+gNs1hmQdgEYAZ8jLFZr9TgUwE8DeBGVkGIZh4iBRA3E1gPXy/+sBzDPY5usAniGiT4jorwCeAXCFEGIigHIiepEkP9fjuv3XAlgKIHd8YAzDMFlEogZiAhF9AADy3/EG25wM4D3N5/fldSfL/+vXQwgxF8CfiOj1BOVjGIZh4iTmSGohxLMATjL4aoXNcwiDdWS2Xgjhl499ua2DC9EIyU2FSZMm2RSJYRiGiUVMA0FEl5l9J4T4sxBiIhF9ILuMPjLY7H0AX9V8PgXAb+T1p+jW7wPwBQA1AF6XY9anAPgfIcT5RPShgXwPAXhIlme/EGJPrGtyiUoAH6foXMmA5U8vLH96YfkjmWy0MqE0VyHEDwD8hYjuFELcBuBEIlqq2+ZEAK8AmCqv+h8A5xHRJ0KIfgD/G8DLAH4F4EdE9Cvd/rsB1BJRRv2YQogBo7SwbIHlTy8sf3ph+e2RaAziTgAzhRDvQMo4uhMAhBC1QohHAICIPgHQAqBfXprldQAQAvAIgF0A3gXwVILyMAzDMC6RUDVXIvoLgBkG6wcA/C/N58cAPGay3VkxzlGdiIwMwzBMfHCpjfh5KN0CJAjLn15Y/vTC8tsgp0ptMAzDMO7BPQiGYRjGEDYQDMMwjCFsICxwoRjh/xNCvC6EeEMI8YAQoiB10icmvxDCL4T4pRDiLVn+O1MpuyxDove/VQjxnhDiUOqkBoQQVwgh3paLUBrVJ/MJITbJ378shKjWfLdMXv+2EOLrqZRbPn9csgshxgoh+oQQh4QQ96Vabo188co/Uwjxilw89BUhxNdSLbssR7zyny+EeE1eXhdCzHdFICLixWQB0A7gNvn/2wCsMdjmRACD8t8x8v9j5O/K5b8CwH8BqM8W+QH4CzmZCQAAA1tJREFUAQTlbYoA/BbArGyRX/7uQgATARxKocwFkFK2T5Pv2+sAztRtcwuAB+T/6wFskv8/U97eB2mw6LsACrJE9lIAlwD4DoD7UvmcuCT/uQCq5P/PglTqJ5vk9wMolP9XBi0XJioT9yCsibsYIQAQ0afyNoWQfvBUZwTELT8RHSGiPgAgomOQBjieYrB/Mkn0/r9Ecq2wFHI+gF1ENCjft42QrkOL9ro2A5ghpLIBVwPYSERHiWgI0vig81MkN5CA7ER0mIi2Afg8deJGkYj8rxLRPnn9GwCKhRC+lEh9nETkP0JEI/L6Yrika9hAWJNIMUIAgBDi15Cs+d8g/aCpJGH5AUAIUQHgKkgl3VOJK/KnGDvyqNvIL/VBAGNt7ptMEpE9E3BL/m8AeJWIjiZJTjMSkl8IcYEQ4g0AvwfwHY3BiJuEBsrlAiJ5xQilf4i+LoQoBtAF4GuQWriukWz5hRCFAH4K4F4iGnQuYYyTJ1n+NGBHHkcFLBOWyD6JyJ4JJCy/EGIKgDWwWSzUZRKSn4heBjBFCPFlAOuFEE8RUUI9urw3EJS8YoTac3wuhOiB1D101UCkQP6HALxDRPe4IG4Uqbj/KeZ9AKdqPitFKI22eV82wAEAn9jcN5kkInsmkJD8QohTADwB4EYiejf54kbhyv0noj8IIQ5DiqUkNAczu5is6QGgZMV8E8DPDbb5NYDLhRBj5CybywH8WghRJis1pRV+JYC3UiCzlrjlBwAhxGpID+CtKZDViITkTxP9AM4QQtQIIYogBRJ7dNtor2shgF6Soos9AOrlTJUaSLMs/i5FcgOJyZ4JxC2/7Eb9JYBlRPRCyiSOJBH5a2Q9AyHEZABfBLA7YYlSHanPpgWSb+85AO/If0+U19cCeESz3bcgBRR3Afgned0E+QffASno9SO4kFWQQvlPgdR1/QOA1+Tlf2WL/PL6dkgtrrD8d1WK5L4SwB8hZaSskNc1A5gr/18M4GeyvL8DcJpm3xXyfm8jxVljLsi+G1Jr9pB8v8/MFvkBfB/AYc2z/hqA8Vkkf4OsZ16DlFAyzw15uNQGwzAMYwi7mBiGYRhD2EAwDMMwhrCBYBiGYQxhA8EwDMMYwgaCYRiGMYQNBMMwDGMIGwiGYRjGkP8PyNd2S9J0FtkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#HF and MP2 errors for all individual four-body terms\n",
    "\n",
    "plt.plot(CCSDTF12b4_tot_F[:],CABS4_val_F[1,:,0]-CCSDTF12b4_tot_F[:],'gx')\n",
    "plt.plot(CCSDTF12b4_tot_F[:],MP24_tot[:]-CCSDTF12b4_tot_F[:],'bo')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 343,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total MP2 error:  0.3918567961405577\n",
      "Total MP3 error:  0.9142704744602432\n"
     ]
    }
   ],
   "source": [
    "# Total MP2 and MP3 errors\n",
    "print('Total MP2 error: ', np.sum(MP_mbe[0,:])-np.sum(CC_mbeQ5[3,:]))\n",
    "print('Total MP3 error: ', np.sum(MP_mbe[1,:])-np.sum(CC_mbeQ5[3,:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimate of MP4 interaction energy:  -4.8314437219583946\n"
     ]
    }
   ],
   "source": [
    "# Estimate of MP4 interaction energy \n",
    "#2b and 3b terms from MP4, 4b terms from HF+CABS\n",
    "print('Estimate of MP4 interaction energy: ', np.sum(MP_mbe[2,1:4])+np.sum(CC_mbeQ5[0,4:7]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MP2 correlation energies for 3b subsets\n",
    "\n",
    "...and how they relate to ATM term.\n",
    "\n",
    "Here we just show in a graph how the MP2 errors for 3b are negative (too weak repulsion) for clusters\n",
    "with hydrogen bonded water molecules and are getting positive as the vertex angle increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 345,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Group  0 30 -0.5575934261112858\n",
      "Group  1 60 -0.2655433970311526\n",
      "Group  2 60 -0.0004569337433210441\n",
      "Group  3 40 0.040912644983374116\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAE/CAYAAACJnoCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3da2xc533n8e9/rrxTJEVdTFmRb4kbO43TMFayeeONr822kYsmjV9sqyxiOC02SLNdbOBsgFXWbgEn6G6KbBdNFUeo2l2kLtzFWmnTCrJTY3eBJDXt+Bo7keJYFi1ZokmK4v02/33BQ2dMz5BDnrmcc+b3AQacc+aZmefwSOfH5znPeY65OyIiIhJNqUZXQERERMpTUIuIiESYglpERCTCFNQiIiIRpqAWERGJMAW1iIhIhGUaXYFq2r59u+/bt6/R1RAREdmUJ5988g137y/1WqKCet++fQwNDTW6GiIiIptiZqfLvaaubxERkQhTUIuIiESYglpERCTCFNQiIiIRpqAWERGJMAW1iIhIhCmoRUREIkxBLSIiEmEKahERkQhTUIuIiERYoqYQFRGRePvkn38fgIJ7Tb/HADPjoc98qKbfUw0KahERaSh359LcEhMzi0zPL7FYcHwTQX1mbAaAy3vbNvW9KTOef22CnvYcPW1Z2nLRjMRo1kpERBJtan6JsakFxmYWGJ9ZYHl5JZg/f8s7N/1ZXz3+EgBfuP3aTb/39Yk5Xp+YAyCfTdHTlqO3feXRkk1v+vNqQUEtIiJ1MTm3yLmJOS5cmmducbnR1Xmb+cXCW4K7syXDjq4Wdne3NDS0FdQiIlJTk3OL/PT8JOPTi42uyqZMzi0xOTfFyyNT7Oxq4Z07O8ll6j8GW6O+RUSkZmYXlhl6ZTx2IV3MfaWLfOj02KbOnVeLWtQikiiro4ZLWVwubGE0sZFLG2ZW8tU4jBpupMn5RZYL9Q+3WpiZX2ZhuUA+U99ucAW1iCTaUsFZXCqwsE5IbzRqeBbIpFPkMimy6RSlI1tKicqArGpIpah7SIOCWkQS5q8+vZ+x6QXemJpndHqBxaXChu/ZzKhhM9jWlqO/I09vRy50fZNufnHj339cuMPCUqHu56kV1CISa4vLBS7OLDI+s8DY9AJTc0s1/T53GJ9eYHx6Ac5DLpOitz0X+WtxG2FxucDJC5ONrkbVuMNLr1/iPQPdZU+F1IL+RYlIrCwtF7g4u8j49EowT9Y4mDeysFQoeS1uT3uOvghdi1tvS8sFfvTqRWbmo3cZVhgXLs3zYmqSd1/WVbfvVFCLSOTNLixzYXKON6bmmZhdpBDh3tS11+K25dJs78zT35FnW1u2ri2xRvr5G9Ncmo3vSO/1nL04y/aOHDu6WuryfQpqEYmsC5NzvDo6w8WZ+B7wZxaWeXV0hldHZ8hnU+zpaePynlYy6WRfHZvkngQz6nqeWkEtIpH0k9cn3xyNnRTziwV+dmGKcxdn+cAVvWQTHNaX97aRz6Q4NTKVqO7vnvYsV+/opLs1W7fvrMq/EjO7w8x+YmanzOzeEq/nzeyh4PUfmtm+ote+GKz/iZndXulnikiyTc039txzLc0tLbO0nIxri9ezo6uFD13Zx/v2bmNXdwvpdDy7/fPZFHv72th/ZS/vf0dvXUMaqtCiNrM08N+BW4Fh4AkzO+buPy4q9mlg3N2vNrO7gK8AnzSzdwN3AdcBlwGPmtnqjOwbfaaIJNh1l3Xx/GsTse72LiWTNq67rJvWXHK7houZGX0defo68hQKztjMyqVzb0wuRHK+71UdLRm2d+Tp78zT1ZJp6NiCanR93wiccveXAczsr4EDQHGoHgC+HDx/GPhTW9nqA8Bfu/s88HMzOxV8HhV8pogkWEs2zfvf0cNrF2d5eWSahQquh44yM9jd3cpVO9obMmlGFKRSxvaOPNs78rBrpddkdGqeN6bmuTizSANm53xTOm30tuXY3pmP3Gj9agT1AHCmaHkY2F+ujLsvmdkE0Bes/8Ga9w4Ezzf6TBFJODNjT08bu7paeHVshtNjM2/eDjFO+jvzXLWjg468hgUV68hn6MhneEdfO4vLBcamFxiZnGdkar4u+zmfTbGjs4XtHTl62nKkUtHsmq/Gv5pSW7b2N1yuTLn1pc6dl9xrZnYPcA/A3r17y9dSRGIrk05xZX8He3raOD06zZnxmUhforWqpz3H1Ts66n5OM46y6RQ7u1rY2dVCoeCMTi/w+sQcI1NzVd3XmbSxq7uF3V2tdLU2tku7UtUI6mHg8qLlPcDZMmWGzSwDdANjG7x3o88EwN0PA4cBBgcH4/entohULJdJcc3OTi7vbePUhak3r1WOmvZ8hmt2dqx08cqmpVJGf+fK+eH5pQ6Gx2c5MzYTagBeay7Nvu3t7OpqIR3RlnM51Rj1/QRwjZldYWY5VgaHHVtT5hhwMHj+ceB7vnKvsGPAXcGo8CuAa4B/rvAzRaRJtWTTXD/QzeC+Htoj1J2cThnv3NnJB6/sVUhXST6T5qr+Dv7FVdvZ1b35CUbM4Mr+dj50ZR8D21pjF9JQhRZ1cM75s8BxIA0ccfcXzOw+YMjdjwHfAv4qGCw2xkrwEpT7G1YGiS0B/9bdlwFKfWbYuopIsmxry7H/il5+NjLF6dHGXnPd3Zblusu6NNd3jeQyKa4f6KYtl+blkemK3pNKwXsGttHfGe8/mqwRN8GulcHBQR8aGmp0NUSkAUYm53n+7MSWBiFt5u5ZpezpbeWdOzojOxgpTta7n/iqmYVl5pd+cWlXuduUtucz5DaYVCYq9xM3syfdfbDUa/rTT0QSob8zzwf29fLsmYvMLNTn+lwzeNeuTvb0lL6PtdRGWy5NJmVvjjB+16633yAjZSRm5jcFtYgkRkc+wweu6OXHZy8xMjlf0+/KZ1O8Z6CbbW26J3U1RaWFGyXJ+HNDRCSQTad47+XbeOfOTlI1OsL1d+bZf0WfQlrqQi1qEUmkvX1t9LRneW54ompd4akUXLOj823nQkVqSS1qEUmszpYsH7iil96O8C3fbCbF+/f2KqSl7tSiFlmj1KhTd6cQ8gKJlBmlJkHSObnayqZT3LBnGz86c5Hx6YUtfUY6bfzK3m10tmiGMak/BbVIBSbnl1heJ6nLXR5SLJNK0dmi/3KNkEoZv7ynmx+8PMr84ubno7zusi6FtDSMjhoia6xt4Z4enebk+al131Ppdbi/dFkXA9taw1VQtiSbTvHu3V386NWLm3rfZdta2dG5+RmxRKpFQS2yjuWC80oVZ7z6+cg0u7taNDFGDW00Ycbs4vLb7oO82iOy+gfXqnTK6GzJlrx70CqdupBaU1CLrOPM2AyLVbwP8tziMmcnZjVBRgO1ZtNvm63qusu6S5ZNma0b0iL1oKAWKWNucZlXRiubU3gzXh6ZZmdXS2JmTYoatXAlaXSkECmhUHBeOHsp1G31yllYKvDC2UskaZ59EakdBbXIGoWC8/zZiS1fylOJNybn+fE5hbWIbExd3yJFFpYKPPfaRcanF2v+XecuzrGwVOD6gW51g4tIWTo6iARGp+b54c9H6xLSv/jOBf7552NcnKld611E4k0taml6S8sFTl6Y4rXx2YZ8/+zCMkOvjLO3r42r+jtI69ItESmioJamdmFyjp+8Prml2aqq7dXRGUYm57l2Vyd9HflGV0dEIkJBLU1pbnGZn56f5MKl2t6zeLNmF5b50asX2dXdwjt3dpLL6OyUSLNTUEvTOXtxlp+en6zJpVfV8vrEHKPTC1y7q5OdXZq+UqSZKailaSwsFXjx3CVGJqPVii5ncanAc8MTjHTP865dnRoZLtKkFNTSFCZmF3l2+GIkzkVv1usTc1yaXeQ9e7p1ByeRJqQ/0SXxLkzO8eTpsViG9KqZhWWGTo8zOhWP3gARqR4FtSTahck5nhueoBDfjH7T8rLzzPBFhbVIk1FQS2JNzy/xwmuXSNIsnYUCPPfaxNtu0ygiyaWglsQ6dWGK5UKCUjqwtOz8bGSq0dUQkTpRUEsiFQrOGwnuIr4Qk5HrIhKegloSadk9UV3eay0vu+68JdIkQgW1mfWa2QkzOxn87ClT7mBQ5qSZHSxa/34ze87MTpnZ183MgvVfNrPXzOzp4PHRMPWU5pNNp2jLpRtdjZrpbMkQ/HcRkYQL26K+F3jM3a8BHguW38LMeoFDwH7gRuBQUaD/GXAPcE3wuKPorV9z9xuCx3dD1lOa0OW9bY2uQs3s7UvutonIW4UN6gPA0eD5UeDOEmVuB064+5i7jwMngDvMbDfQ5e7f95U+vL8s836RLdnT00pXa/ImCOlpz7FL04qKNI2wQb3T3c8BBD93lCgzAJwpWh4O1g0Ez9euX/VZM3vWzI6U61IXWY+Zcf1AF+l0crqIs5kU113WpW5vkSayYVCb2aNm9nyJx4EKv6PUEcXXWQ8rXeJXATcA54D/sk797jGzITMbGhkZqbBK0izachmuu6yr0dWoCjN4z0A3LdnknnsXkbfbcK5vd7+l3Gtmdt7Mdrv7uaAr+0KJYsPATUXLe4DHg/V71qw/G3zn+aLv+Cbwd+vU7zBwGGBwcFDDYOVtdnS28I6+RU6PzjS6KqFc2d9Bb3uu0dUQkToL2/V9DFgdxX0QeKREmePAbWbWE3Rh3wYcD7rKJ83sg8Fo799ZfX8Q+qt+A3g+ZD2lyV3V30FHS3zvQbOtLcs+DSATaUphg/oB4FYzOwncGixjZoNm9iCAu48B9wNPBI/7gnUAvwc8CJwCfgb8Q7D+q8FlW88C/xL4dyHrKU0ulTLeHdMu8FQKfmm3zkuLNKtQTQx3HwVuLrF+CLi7aPkIcKRMuetLrP/tMPUSKaWrJctATyuvjc82uiqbsre3jfZ8fHsDRCQczUwmTeXqHR3kMvH5Z9+STbOvr73R1RCRBorPEUukCrLpFL+0Ox5d4GZw3WVdZNL6byrSzHQEkKbT35nn6h0dja7Ght61q5MejfIWaXo68SVNad/2dlJm/PT8ZKOr8jZmcO3uLga2tTa6KiISAQpqaVp7+9poz6d54ewlFpYKja4OAPlsiusv61ZLWkTepK5vaWp9HXk+eGUfu7obP3f2Zdta+eCVfQppEXkLtail6eUyKa4f6GZgWys/OT/J1NxSXb+/qzXLu3Z20t2WvBuIiEh4CmqRQE97jv1X9HJ2Yo6XR6aYX6xtd3hrLs1V/R3s7MprMhMRKUtBLVLEzPiDh54GYG5xmbmlAit3YV3fmbGVecS/evylir6jNZsmH1zP/dBnPhSixiKSdApqkTJasmny2XRFQd090F3x55pZyVvHiYiUoqAWWUMtXBGJEo36FhERiTAFtYiISIQpqEVERCJMQS0iIhJhCmoREZEIU1CLiIhEmC7PqoFP/vn36/p9upxIRCS5FNRVVCg43395lIszC5t+76vBzFZ7e9s2/97RGfb2bf59IiISfQrqKnrt4iyzC8v8h9uv3fR7V6ee3Mp7XxmdZqCnlXRK812JiCSNzlFXSaHgvDI63ZDvXlgq8Nr4bEO+W0REaktBXSWvXZyt+d2W1vPK6DTLhY3npBYRkXhRUFdBI1vTq9SqFhFJJgV1FTS6Nb1KrWoRkeRRUIe0HIHW9KqFpQLD4zONroaIiFSRgjqksxFpTa96ZXRGrWoRkQRRUIcQhXPTay2qVS0ikiihgtrMes3shJmdDH72lCl3MChz0swOFq3/IzM7Y2ZTa8rnzewhMztlZj80s31h6lkrUTk3vdZptapFRBIjbIv6XuAxd78GeCxYfgsz6wUOAfuBG4FDRYH+nWDdWp8Gxt39auBrwFdC1rPqotiaXqUR4CIiyRE2qA8AR4PnR4E7S5S5HTjh7mPuPg6cAO4AcPcfuPu5DT73YeBmM4vUtFtnJ6LZml51emyaglrVIiKxFzaod64GbfBzR4kyA8CZouXhYN163nyPuy8BE0BfyLpWjbtzejTa54HnFwucnVCrWkQk7jac69vMHgV2lXjpSxV+R6mW8EZNvYrfY2b3APcA7N27t8IqhXP+0jyzC8t1+a4wTo/OMLCtlYh1RoiIyCZsGNTufku518zsvJntdvdzZrYbuFCi2DBwU9HyHuDxDb52GLgcGDazDNANjJWp32HgMMDg4GBd+nqjem56rdmFZS5MzrOzq6XRVRERkS0K2/V9DFgdxX0QeKREmePAbWbWEwwiuy1YV+nnfhz4nrtH4oTr2PQCU3NLja5GxaLeRS8iIusLG9QPALea2Ung1mAZMxs0swcB3H0MuB94InjcF6zDzL5qZsNAm5kNm9mXg8/9FtBnZqeAP6DEaPJGOR2T1vSqS7OLW7o/toiIREOo+1G7+yhwc4n1Q8DdRctHgCMlyn0B+EKJ9XPAJ8LUrRZmFpYYnYpf6J0Zm2VbW67R1RARkS3QzGSbMBzTa5MvTM4xtxj9wW8iIvJ2CuoKLRecsxfjGdTuxLbuIiLNTkFdoQuTcywtR2I825acvThHRMbjiYjIJiioK3T24lyjqxDK3OIy4zOLja6GiIhskoK6AnOLy4xPx28Q2VrnNFOZiEjsKKgrcOHSfKOrUBUXJuc1/7eISMwoqCtwfjLe3d6rlped0QT0DIiINBMF9Qbml5aZSNC53ZHJZPQOiIg0CwX1BuI4wcl63phSUIuIxImCegNjCesqXlgqMDUfn7nKRUSanYJ6A0kLaiARI9hFRJqFgnodMwtLLCwVGl2NqruYoHPuIiJJp6Bex6XZZHYRX5pTUIuIxIWCeh2TCQ202YVlFpeT11MgIpJECup1JHnQ1XSCt01EJEkU1OuYXUjurSGnE7xtIiJJoqAuw92ZW0pumOn+1CIi8aCgLmNhuUAhwadx5xcTvHEiIgmioC5jMcb3nq6EBpOJiMSDgrqMpYQH2VKSuwtERBJEQV3GcsJvB5nwv0NERBJDQV1GsmMaCp70LRQRSQYFdRnW6ArUWNK3T0QkKRTUZZglO8pSqWRvn4hIUiioy0gnPMiSvn0iIkmhoC4jm052kGUU1CIisaCgLiObTvavJpdJ9vaJiCRFqKO1mfWa2QkzOxn87ClT7mBQ5qSZHSxa/0dmdsbMptaU/5SZjZjZ08Hj7jD13IpsOpXo7uF8Jt3oKoiISAXCNqvuBR5z92uAx4LltzCzXuAQsB+4EThUFOjfCdaV8pC73xA8HgxZzy3JZ5Pb6mxJ8LaJiCRJ2KP1AeBo8PwocGeJMrcDJ9x9zN3HgRPAHQDu/gN3PxeyDjXTlss0ugo105ZN7raJiCRJ2KDeuRq0wc8dJcoMAGeKloeDdRv5TTN71sweNrPLQ9ZzS9pzye0ebssnd9tERJJkw2aVmT0K7Crx0pcq/I5SJ3o3mhbrO8C33X3ezH6Xldb6R8rU7x7gHoC9e/dWWKXKtOeT2erMZ1OJHywnIpIUGyaRu99S7jUzO29mu939nJntBi6UKDYM3FS0vAd4fIPvHC1a/CbwlXXKHgYOAwwODlZ1XszOlmQGdUdC/wAREUmisM2qY8DqKO6DwCMlyhwHbjOznmAQ2W3BurKC0F/1MeDFkPXckvZchlQCG55drdlGV0FERCoUNoYeAG41s5PArcEyZjZoZg8CuPsYcD/wRPC4L1iHmX3VzIaBNjMbNrMvB5/7OTN7wcyeAT4HfCpkPbcklTI6W5IXal0J3CYRkaQK1QcadFHfXGL9EHB30fIR4EiJcl8AvlBi/ReBL4apW7V0t2aZmFlsdDWqqlstahGR2Ehgx251bWtLVqi15dOalUxEJEZ0xN7AttZco6tQVT1tydoeEZGkU1BvIJdJ0ZGg0d8KahGReFFQV6C3PTnh1tOerK58EZGkU1BXICmt0I6WjG7GISISMwrqCvS0ZRNxPXVfgnoGRESaRQLip/Yy6VQiLmlKUhe+iEizUFBXqK893+gqhJJOWWK68EVEmomCukJ9HfEOuZ72HKlUqfujiIhIlCmoK9TZkiWfje+vS+enRUTiKb7J0wBxPscb9x4BEZFmpaDehO0d8TxP3ZZL05ZLzqQtIiLNREG9Cb3tOSyGp3n7YvoHhoiIKKg3JZtOxfJeznHushcRaXYK6k2KW+ilUisTtoiISDwpqDcpbqOnu1uzZNLazSIicaUj+CZ1tWRJp+Nzoro35hO1iIg0OwX1JqViNsOXur1FROJNQb0FvTEJ6nTK6GpRUIuIxJmCegu2xeSezt1tWU0bKiIScwrqLejMZ2JxnjpOXfQiIlKagnoLzCwWt73cFoM6iojI+hTUWxT1EDQjlpOziIjIWymotyjqLeqOfIa0zk+LiMSegnqLot5ajXr9RESkMgrqLcqmU7Tm0o2uRlkKahGRZFBQh9DZEt1bR0a5biIiUjkFdQidEZ1MxAw6dP9pEZFECBXUZtZrZifM7GTws6dMuYNBmZNmdjBY12Zmf29mL5nZC2b2QFH5vJk9ZGanzOyHZrYvTD1rpT0fza7vtlxGE52IiCRE2Bb1vcBj7n4N8Fiw/BZm1gscAvYDNwKHigL9j939WuB9wIfN7FeD9Z8Gxt39auBrwFdC1rMmOvLRbLVGtV4iIrJ5YYP6AHA0eH4UuLNEmduBE+4+5u7jwAngDnefcfd/AnD3BeApYE+Jz30YuNnMItdEbM2mSUXw5EFbRFv6IiKyeWFjZqe7nwMIfu4oUWYAOFO0PByse5OZbQN+nZVW+Vve4+5LwATQV6oCZnaPmQ2Z2dDIyEiITdk8M6MtgueC2yNYJxER2ZoNj+hm9iiwq8RLX6rwO0q1hL3o8zPAt4Gvu/vLlbznLSvdDwOHAQYHB0uWqaW2XJqpuaV6f+261KIWEUmODYPa3W8p95qZnTez3e5+zsx2AxdKFBsGbipa3gM8XrR8GDjp7n+y5j2XA8NBkHcDYxvVtRFas9ELxSjWSUREtiZs1/cx4GDw/CDwSIkyx4HbzKwnGER2W7AOM/tDVkL48+t87seB77l73VvLlYjapCeZtJFNR/DEuYiIbEnYI/oDwK1mdhK4NVjGzAbN7EEAdx8D7geeCB73ufuYme1hpfv83cBTZva0md0dfO63gD4zOwX8ASVGk0dF1FqvUauPiIiEE2rUkbuPAjeXWD8E3F20fAQ4sqbMMKXPRePuc8AnwtStXloiFoxRq4+IiISjPtKQ8plo/QoV1CIiyRKtlImhTDpFOh2dS7yj9oeDiIiEo6N6FUQpHPPZ6NRFRETC01G9CvKZ6HQ3R6kuIiISnoK6CqLUos5FqC4iIhKejupVEKVwzOkaahGRRNFRvQqiEo5mkI3QwDYREQkvGgkTc1FpUecyKSJ4kzEREQkhGgkTc1GZsjMq9RARkerRkb0KotKiVlCLiCSPjuxVEJVz1FEafS4iItWhI3sVRGUAl1rUIiLJoyN7FWTSKVIR+E1G5Q8GERGpngjESzLk0o2fESwq58pFRKR6dGSvkii0ZqNyrlxERKpHR/YqiUJrNgp1EBGR6tKRvUqiEJJRqIOIiFSXjuxVEoVuZ436FhFJHh3Zq6TRrdlUSkEtIpJEOrJXSaPvAx2FUeciIlJ9CuoqaXSLutHfLyIitaGje5U0Oig1faiISDLp6F4ljQ7KfFa7UkQkiXR0r5JsOkU61bhJTxp9jlxERGpDQV1FjWxVN7pFLyIitaGjexU1svtZQS0ikkw6uldRI7ufW7Lq+hYRSaJQQW1mvWZ2wsxOBj97ypQ7GJQ5aWYHg3VtZvb3ZvaSmb1gZg8Ulf+UmY2Y2dPB4+4w9ayXRoalglpEJJnCtqjvBR5z92uAx4LltzCzXuAQsB+4EThUFOh/7O7XAu8DPmxmv1r01ofc/Ybg8WDIetZFS4O6vrOZxg5kExGR2gmbLAeAo8Hzo8CdJcrcDpxw9zF3HwdOAHe4+4y7/xOAuy8ATwF7QtanoRrVqm3R+WkRkcQKe4Tf6e7nAIKfO0qUGQDOFC0PB+veZGbbgF9npVW+6jfN7Fkze9jMLi9XATO7x8yGzGxoZGRkq9tRFa2NCmp1e4uIJNaGQW1mj5rZ8yUeByr8jlJ9sl70+Rng28DX3f3lYPV3gH3u/svAo/yi1f72D3I/7O6D7j7Y399fYZVqo1GB2ZpTUIuIJFVmowLufku518zsvJntdvdzZrYbuFCi2DBwU9HyHuDxouXDwEl3/5Oi7xwtev2bwFc2qmcUpFNGLpNiYalQ1+9tVEteRERqL2zX9zHgYPD8IPBIiTLHgdvMrCcYRHZbsA4z+0OgG/h88RuC0F/1MeDFkPWsm0a0btX1LSKSXGGD+gHgVjM7CdwaLGNmg2b2IIC7jwH3A08Ej/vcfczM9gBfAt4NPLXmMqzPBZdsPQN8DvhUyHrWTSNat+r6FhFJrg27vtcTdFHfXGL9EHB30fIR4MiaMsOUPn+Nu38R+GKYujVKI0JTXd8iIsml63qqrN6hmc/qGmoRkSRTUFdZW51b1GpNi4gkm4K6yuo9sEvnp0VEkk1BXWUt2XRdu6LbcqGGGYiISMQpqGugnq1cdX2LiCSbgroG6hme6voWEUk2BXUN1HNAWb0Hr4mISH0pqGugXq3cTNrIprULRUSSTEf5GqhX17cGkomIJJ+CugbqFaDq9hYRST4FdQ20ZFOk6vCb1c04RESST0FdA2ZGS6b2IaoWtYhI8imoa6QeA8oU1CIiyaegrpF6BLW6vkVEkk9BXSO1HvmdSkE+o90nIpJ0OtLXSK1b1C3ZNGa6vaWISNIpqGuk1i1qzfEtItIcFNQ1Uuvzx5rjW0SkOSioaySbTpFJ165rWi1qEZHmoKCuoVq2qjXiW0SkOSioa6iWrV4FtYhIc1BQ11BtW9TadSIizUBH+xqqVYt65RpqtahFRJqBgrqGatXqrcc84iIiEg0K6hrK16hF3aJLs0REmoaCuobUohYRkbAU1DWUS9fmvtR5DSQTEWkaoY/4ZtZrZifM7GTws6dMuYNBmZNmdrBo/T+a2TNm9oKZfcPM0pv53Cir1X2pdWmWiEjzqEbT7F7gMXe/BngsWH4LM+sFDgH7gRuBQ0XB+1vu/l7geqAf+ESlnxsHtWj9tuiuWSIiTaMaR/wDwNHg+VHgzhJlbgdOuPuYu48DJ4A7ANz9UvZdsvMAAAjuSURBVFAmA+QA38TnRl4tLqNSi1pEpHlUI6h3uvs5gODnjhJlBoAzRcvDwToAzOw4cAGYBB7exOdiZveY2ZCZDY2MjITdlqqrxYAy3YdaRKR5VHTEN7NHzez5Eo8DFX5PqbtT+JtP3G8HdgN54CMVfubqew+7+6C7D/b392/mrXVR7RZ1Jm1k0gpqEZFmkamkkLvfUu41MztvZrvd/ZyZ7WalZbzWMHBT0fIe4PE13zFnZsdY6fI+AVTyuZFX7XPUmpFMRKS5VCNFjgGro7gPAo+UKHMcuM3MeoJBZLcBx82sIwhhzCwDfBR4aROfG3nVPp+sOb5FRJpLNY76DwC3mtlJ4NZgGTMbNLMHAdx9DLgfeCJ43BesaweOmdmzwDOstJq/sd7nxk21L89Si1pEpLlU1PW9HncfBW4usX4IuLto+QhwZE2Z88AHNvO5cZPLrEx6UihU5/PUohYRaS466tdBNVvVtZo/XEREoklBXQfVHFCmyU5ERJqLjvp1UM3zyprsRESkuSio66Ca4aqgFhFpLgrqOqjWALBM2kinSs0dIyIiSaWgroNqtYLVmhYRaT4K6jqoVsC2KqhFRJqOgroOqjVSWy1qEZHmo6Cug0w6RSYd/tyyJjsREWk+OvLXSTW6rdX1LSLSfBTUddKaCx+yLVX4DBERiRcFdZ1UozVc7Rt8iIhI9Cmo6yTsQLB02shp+lARkaajI3+dhO361vlpEZHmpKCuk7BB26bz0yIiTUlBXSdhg1otahGR5qSgrpNUykKdp67GqHEREYkfBXUdhQnbtlymijUREZG4UFDXUZjzzDpHLSLSnBTUdbTVsE2lIK9Ls0REmpKO/nW01e7r1mwGM92HWkSkGSmo62irLer2vLq9RUSalUYo1cAn//z7ZV+bml8quX54fBaArz3607e9lkun1p2V7KHPfGiTNRQRkbhQUNdZR770r/w9A911romIiMSBgroG1MIVEZFq0TlqERGRCFNQi4iIRFiooDazXjM7YWYng589ZcodDMqcNLODRev/0cyeMbMXzOwbZpYO1n/ZzF4zs6eDx0fD1FNERCSuwrao7wUec/drgMeC5bcws17gELAfuBE4VBTov+Xu7wWuB/qBTxS99WvufkPw+G7IeoqIiMRS2KA+ABwNnh8F7ixR5nbghLuPufs4cAK4A8DdLwVlMkAO8JD1ERERSZSwQb3T3c8BBD93lCgzAJwpWh4O1gFgZseBC8Ak8HBRuc+a2bNmdqRcl3rw/nvMbMjMhkZGRkJsioiISPRsGNRm9qiZPV/icaDC7yg19+WbLWd3vx3YDeSBjwSr/wy4CrgBOAf8l3If7u6H3X3Q3Qf7+/srrJKIiEg8bHgdtbvfUu41MztvZrvd/ZyZ7WalZbzWMHBT0fIe4PE13zFnZsdY6Uo/4e7ni77jm8DfbVRPERGRJArb9X0MWB3FfRB4pESZ48BtZtYTdGHfBhw3s44g3DGzDPBR4KVgeXfR+38DeD5kPUVERGIp7MxkDwB/Y2afBl4lGLVtZoPA77r73e4+Zmb3A08E77kvWLcTOGZmeSANfA/4RlDmq2Z2Aytd5K8AnwlZTxERkVgy9+QMtB4cHPShoaFGV0NERGRTzOxJdx8s+VqSgtrMRoDTja5HCNuBNxpdiRrS9sWbti/etH3R9g53LzkiOlFBHXdmNlTuL6ok0PbFm7Yv3rR98aW5vkVERCJMQS0iIhJhCupoOdzoCtSYti/etH3xpu2LKZ2jFhERiTC1qEVERCJMQV1nwU1GLphZydnWbMXXzexUcFOSX6l3HcOoYPtuMrOJonuN/6d61zEMM7vczP7JzF4M7qP++yXKxHYfVrh9sd2HZtZiZv9sZs8E2/efS5TJm9lDwf77oZntq39Nt6bC7fuUmY0U7b+7G1HXMMwsbWY/MrO3TS8d5/1XTtiZyWTz/gL4U+Avy7z+q8A1wWM/Kzco2V+XmlXHX7D+9gH8X3f/tfpUp+qWgH/v7k+ZWSfwpJmdcPcfF5WJ8z6sZPsgvvtwHviIu0+ZWRb4f2b2D+7+g6IynwbG3f1qM7sL+ArwyUZUdgsq2T6Ah9z9sw2oX7X8PvAi0FXitTjvv5LUoq4zd/8/wNg6RQ4Af+krfgBsWzP3eaRVsH2x5u7n3P2p4PkkKweLgTXFYrsPK9y+2Ar2yVSwmA0eawfqHACOBs8fBm42s1J3AYycCrcv1sxsD/CvgAfLFInt/itHQR09696/OyE+FHTN/YOZXdfoymxV0KX2PuCHa15KxD5cZ/sgxvsw6DZ9mpW7/Z1w97L7z92XgAmgr7613LoKtg/gN4PTMg+b2eV1rmJYfwJ8ASiUeT3W+68UBXX0rHv/7gR4ipWp8t4L/Dfgfze4PltiZh3A3wKfd/dLa18u8ZZY7cMNti/W+9Ddl939BlZuuXujmV2/pkis918F2/cdYJ+7/zLwKL9ofUaemf0acMHdn1yvWIl1sdl/pSioo2cYKP4Ldw9wtkF1qTp3v7TaNefu3wWyZra9wdXalODc398C/9Pd/1eJIrHehxttXxL2IYC7XwQeB+5Y89Kb+y+4BW83MTydU2773H3U3eeDxW8C769z1cL4MPAxM3sF+GvgI2b2P9aUScT+K6agjp5jwO8EI4c/CEy4+7lGV6pazGzX6vkiM7uRlX+Do42tVeWCun8LeNHd/2uZYrHdh5VsX5z3oZn1m9m24HkrcAvw0ppix4CDwfOPA9/zmEw4Ucn2rRkv8TFWxiHEgrt/0d33uPs+4C5W9s2/XlMstvuvHI36rjMz+zZwE7DdzIaBQ6wM+MDdvwF8F/gocAqYAf5NY2q6NRVs38eB3zOzJWAWuCtm/4k+DPw28FxwHhDgPwJ7IRH7sJLti/M+3A0cNbM0K39g/I27/52Z3QcMufsxVv5Q+SszO8VKS+yuxlV30yrZvs+Z2cdYGeE/BnyqYbWtkgTtv5I0M5mIiEiEqetbREQkwhTUIiIiEaagFhERiTAFtYiISIQpqEVERCJMQS0iIhJhCmoREZEIU1CLiIhE2P8HGYX40peQp+cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#calculate the MP2 errors wrt CCSD(T) for three-body groups\n",
    "#CABS3_val_F_NN=np.zeros((190,4))\n",
    "\n",
    "diff=np.zeros((190,4))\n",
    "for i in range(4):\n",
    "    for j in range(idx3_ok[i]):\n",
    "        diff[j,i]=MP23_val_NN[j,i]-CCSDTF12b3_tot_F_NN[j,i]\n",
    "    print('Group ', i, idx3_ok[i], np.sum(diff[:,i]) )\n",
    "    \n",
    "fig=plt.figure()\n",
    "ax2 = fig.add_axes([0,0,1,1])\n",
    "data=[diff[0:idx3_ok[0],0], diff[0:idx3_ok[1],1], diff[0:idx3_ok[2],2], diff[0:idx3_ok[3],3]]\n",
    "bp2=ax2.violinplot(data)\n",
    "plt.show()\n",
    "#data+=[MP23_val_NN[0:idx3_ok[0],0],diff3_mp2_NNord[0:idx3_ok[1],1],diff3_mp2_NNord[0:idx3_ok[2],2],diff3_mp2_NNord[0:idx3_ok[3],3]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two groups of dimers\n",
    "\n",
    "This is pyplot version of Figure 1, errors of MPn methods versus CCSD(T)-F12 reference data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 346,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAa2ElEQVR4nO3df5Dc9X3f8efrTjkhmYCxUCKQkA4V8oecMFK44kBtkzGNgbax7AmeCF8cmSGVKaHT1k0bUUxt0GgyeKZhkqBOpQlRiJGDCJ42moCjcUrN1LZMOFkyIIhiRZWE0MkRgRALRZx1eveP7/dgb7V7+93bX9/vd1+PmZvd/X6/u/e+z96+9rOf7+f7XUUEZmZWXgO9LsDMzDrLQW9mVnIOejOzknPQm5mVnIPezKzk5vS6gGoXX3xxDA8P97oMM7NC2b1792sRsbDWutwF/fDwMGNjY70uw8ysUCQdrrfOQzdmZiXnoDczKzkHvZlZyTnozcxKLlPQS7pJ0n5JByStr7H+w5K+K+mMpFuq1q2V9P30Z227Cjczs2waBr2kQWATcDOwArhV0oqqzY4AnwG+UnXf9wFfAD4AXAN8QdJFrZdtVgDbtsHwMAwMJJfbtvW6IutTWXr01wAHIuJgREwAjwGrKzeIiEMR8Txwtuq+NwJfj4jXI+IN4OvATW2o2yzftm2Ddevg8GGISC7XrXPYW09kCfrFwCsVt4+my7LIdF9J6ySNSRo7ceJExoc2y7F77oFTp6YvO3UqWW7WZVmCXjWWZT2Jfab7RsSWiBiJiJGFC2se2GVWLEeONLfcrIOyBP1R4LKK20uAYxkfv5X7mhXX0qXNLTfroCxB/xxwpaTLJQ0Ba4AdGR9/J/BRSRelO2E/mi4zK7eNG2H+/OnL5s9Plpt1WcOgj4gzwF0kAf0y8HhE7JN0v6SPAUj6p5KOAp8ENkval973dWADyZvFc8D96TKzchsdhS1bYNkykJLLLVuS5WZdprx9Z+zIyEj4pGZm9o5t25Kd2EeOJENfGzf6DbMGSbsjYqTWutydvdLM7B1T01SnZjBNTVMFh30TfAqEMvMBO1Z0nqbaFu7Rl5V7QlYGnqbaFu7Rl5V7QlYGnqbaFg76snJPyMrA01TbwkFfVu4JWRl4mmpbOOjLyj0hK4vRUTh0CM6eTS4d8k1z0JeVe0LW7zzr7B2edVNmo6MOdutPnnU2jXv0ZlY+nnU2jYPezMrHs86mcdCbWfl41tk0DnozKx/POpvGQW9m5eNZZ9N41o2ZlZNnnb3DPXozm1mW+eies55r7tGbWX1Z5qN7znru+RumzKy+4eEkuKstW5acjiDrNtZxM33DlIduzKy+LPPRPWc99xz0ZlZflvnonrOeew56M6svy3x0z1nPPQe9mdWXZT6656znnnfGmpmVgHfGmpn1sXIFvQ/a6Ay3q1mhleeAKR+00RluV7PCK88YvQ/a6Ay3q1kh9McYvQ/a6Ay3q1nhZQp6STdJ2i/pgKT1NdbPlbQ9Xf+spOF0+Y9JekTSC5JelnR3e8uv4IM2OsPtalZ4DYNe0iCwCbgZWAHcKmlF1Wa3A29ExBXAg8AD6fJPAnMj4meAq4HPTr0JtJ0P2ugMt6tZ4WXp0V8DHIiIgxExATwGrK7aZjXwSHr9CeAGSQICeI+kOcA8YAL4h7ZUXs0HbXSG29Ws8LLMulkMvFJx+yjwgXrbRMQZSW8CC0hCfzUwDswH/kNEvF79CyStA9YBLG1lSMBfNNAZblezQsvSo1eNZdVTdeptcw0wCVwKXA78R0nLz9kwYktEjETEyMKFCzOUZGZmWWUJ+qPAZRW3lwDH6m2TDtNcCLwOfAr484j4UUT8LfAtoOb0HzMz64wsQf8ccKWkyyUNAWuAHVXb7ADWptdvAZ6OZIL+EeAjSrwH+Dngr9pTupmZZdEw6CPiDHAXsBN4GXg8IvZJul/Sx9LNHgYWSDoAfA6YmoK5CTgfeJHkDWNrRDzf5r/BrPh8mgnroPIcGWtWVNWnmYBkCqtnN1kT+uPIWLOiuuee6SEPye177ulNPVY6DnqzXvNpJqzDHPRmvebTTFiHOejNes2nmbAOc9Cb9ZpPM2EdVp4vHjErMp9mwjrIPXozs5Jz0JuZlZyD3sys5Bz0ZmYl56A3Mys5B309PsmUmZWEp1fWUn2SqcOHk9vgKXBmVjju0dfik0yZWYk46GvxSabMrEQc9LX4JFNmViIO+lp8kikzKxEHfS0+yZSZlYhn3dTjk0yZWUm4R29mVnIOejOzknPQm5mVnIPezKzkHPRmZiXnoDczKzkHvZlZyTnozcxKzkFvZlZymYJe0k2S9ks6IGl9jfVzJW1P1z8rabhi3VWSdknaJ+kFSee1r3yzfBsfh+uvh+PHe12J9bOGQS9pENgE3AysAG6VtKJqs9uBNyLiCuBB4IH0vnOAR4E7IuL9wM8DP2pb9WY5t2EDfPObyaVZr2Tp0V8DHIiIgxExATwGrK7aZjXwSHr9CeAGSQI+CjwfEd8DiIi/i4jJ9pRulm/j47B1K5w9m1y6Vz97/mTUmixBvxh4peL20XRZzW0i4gzwJrAA+CkgJO2U9F1J/7nWL5C0TtKYpLETJ040+zeY5dKGDUnIA0xOulffCn8yak2WoFeNZZFxmznAB4HR9PITkm44Z8OILRExEhEjCxcuzFCSWb5N9eYnJpLbExPu1c+WPxm1LkvQHwUuq7i9BDhWb5t0XP5C4PV0+TMR8VpEnAKeAn621aLN8q6yNz/FvfrZ8Sej1mUJ+ueAKyVdLmkIWAPsqNpmB7A2vX4L8HREBLATuErS/PQN4HrgpfaUbpZfu3a925ufMjEB3/52b+opqlY/GXlsP9Ew6NMx97tIQvtl4PGI2CfpfkkfSzd7GFgg6QDwOWB9et83gN8mebPYC3w3Ip5s/59hli979kDEuT979vS6ss7pRKi2+snIY/sJJR3v/BgZGYmxsbFel2FmTbrzTti8Ge64AzZtas9jrloFe/eeu3zlysZvmuPjsHw5nD4N8+bBwYOwaFF76sojSbsjYqTWOh8ZW3L+6Grd0Kkdpq18MvLY/rsc9CXnj67WDXkLVc96ms5BX2KelmbtMtMnwzyGqmc9TeegL7G89bJsZnkeZpvpk2EeQ9WznqZz0JdUHntZNrO8DrM1+mSYx1Dtx1lPM3HQl1Qee1lWX56H2Rp9MnSo5p+DvqTy2Muy+vI6zOZPhuXgoC8p97KKI89h6k+G5eCgN+uxPIepPxmWw5xeF2DW7/Icpv4EWA4OerMec5hap3noxjLJ8xxvM5uZg94yyescbzNrzEFvDeV5jreZNeagt4byOsfbzLJx0NuM8jzH28yyKV3Qe6dhe+V5jreZZVO6oPdOw/bK8xxvM8umVF8l2G9fHWZmNqVvvkrQOw3NzM5VmqD3TkMzs9pKE/TeaWhmVltpgr4TOw09g8fMyqA0Qd+J8697Bo+ZlUFpgr7dfNi/mZWFg74Oz+Axs7Jw0NfgGTxmViYO+ho8g8esWDxxYmaZgl7STZL2SzogaX2N9XMlbU/XPytpuGr9UkknJf1Ge8ruLB/2b1Ysnjgxs4ZBL2kQ2ATcDKwAbpW0omqz24E3IuIK4EHggar1DwJfa73c7ujEDB4z6wxPnGgsS4/+GuBARByMiAngMWB11TargUfS608AN0gSgKSPAweBfe0p2czsXZ440ViWoF8MvFJx+2i6rOY2EXEGeBNYIOk9wG8C9830CyStkzQmaezEiRNZazezPueJE9lkCXrVWFZ9yst629wHPBgRJ2f6BRGxJSJGImJk4cKFGUoyM0t782cmpy2b/NGke/VVsgT9UeCyittLgGP1tpE0B7gQeB34APAlSYeAfw/8F0l3tVizmRkAu556nYkzg9OWTZwZ5NtPvt6jivJpToZtngOulHQ58CqwBvhU1TY7gLXALuAW4OlITnT/oakNJH0ROBkRD7WhbjMz9vCzwOEaa5YBh7pbTI41DPqIOJP2wncCg8AfRMQ+SfcDYxGxA3gY+LKkAyQ9+TWdLNrMDIAjR5pb3qdK9Q1TZtZnhofhcI0e/bJlcOhQt6vpqb75hikz6zMbN8L8+dOXzZ+fLLd3OOjNrGvafqqC0VHYsiXpwUvJ5ZYtyXJ7h4PerA/k5VwwHTlVwehoMkxz9mxy6ZA/h4Pe8mXbtmTcdWAgudy2LfNd8xJmeZSHc8H4VAW946C3/Ni2DdatS3auRSSX69ZlDvs8hNlsdPoNKi8B61MV9I6D3vLjnnvg1CnGWcT1fIPj/CScOpUsbyAvYTYbnX6DykPA+lQFveWgt/xI5z5v4F6+yQfZwL3Tls8kD2E2G51+g8pLwPo7HnrLQW/5sXQp4yxiK7dxlkG2clvSq1+6dMa75SXMZqPTb1B5CVh/x0NvOegtPzZuZP3Al3ibuQBMMsCGwfsazonOS5g1qxtvUHkJWH/HQ2856C03s1XGPzLKthgl0n/LCc5j68DtHL9h5ulyeQmzZnXjDcoBa+CgN/IzW2X9epiM6f+Sk5rTsK6ihllR36CseBz0fS5Ps1WefPLcZWUOvqK+Qc1aC8dIWGsc9H0uL7NVxsfhrbemL5s3L1le2uDrJy0eI2GtcdD3sTzNVinqDlXLKD1GYpqMx0hY6/om6Ltx9GEedmg2I0/h6vHqkvN543uqb4K+G0cf5mGHZjPyFK59N17db+odC9HgGAlrj74I+m4dfZiHHZrNcLha1/i88T3VF0G/4dN/zdnTbwMw+Y9vs+FX9rf38XOyQ9Mst3ze+J4q/VcJjj/0VZb/23/Baea9s2wepzj40NdY9Ou/1Prjj8Py5XD69LvL5s2Dgwdh0aKWH97MLJO+/irBDXef4iyatmySATasf6vOPZp8/Bzt0DQzq6X0Qb/r5E8zwXnTlk1wHt8++TPtefwc7dA0M6tlTq8L6LQ9yz5R/1viOdT643vHpZnlXOl79N7bb33Ppx7oe+UPeu/tt37mUw8YfTDrxqyvDQ/XH7o8dKjb1VgH9fWsG7O+5lMPGA769vAYaO8V9TnodN0+9YDhoG+dx0DbZ7ahV9TnoBt1ezKCAUREwx/gJmA/cABYX2P9XGB7uv5ZYDhd/gvAbuCF9PIjjX7X1VdfHYWybFmt08Ukyy27Rx+NmD9/ehvOn58sb6Soz0G36n700eQxpeQyS5ta4QBjUSdXG+6MlTQI/HUa2keB54BbI+Klim3uBK6KiDskrQE+ERG/LGkV8IOIOCbpp4GdEbF4pt9XuJ2xAwPJy7OadO4hs1ZfKzsNi/ocFLVuy6VWd8ZeAxyIiIMRMQE8Bqyu2mY18Eh6/QngBkmKiD0RcSxdvg84T9Lc5v+EHCv6GGhexrZb2WlY1OegqHVb4WQJ+sXAKxW3j6bLam4TEWeAN4EFVdv8ErAnIt6u/gWS1kkakzR24sSJrLXnQ5HHQPM0tt1K6BX1OShq3VY89cZ0pn6ATwK/X3H708DvVW2zD1hScftvgAUVt9+fLvsnjX5f4cboI4o7Bpqnse1Wxuin7l/E56CodVvu0OIY/bXAFyPixvT23ekbxG9VbLMz3WaXpDnAcWBhRISkJcDTwG0R8a1GbzyFG6MvsryNEW/blnyH6JEjSU9+40YfwWyWUatj9M8BV0q6XNIQsAbYUbXNDmBtev0W4Ok05N8LPAncnSXkrcvyNkY8OprseD17Nrl0yJu1RcOgj2TM/S5gJ/Ay8HhE7JN0v6SPpZs9DCyQdAD4HLA+XX4XcAVwr6S96c9PtP2vsNnxGLFZX/C5bvqdh0vMSmGmoZvSn4/eGhgddbCblZxPgWBmVnIOejOzknPQm5mVnIPezKzkHPRmZiXnoDczKzkHvZlZyTnozcxKzkFvZlZyDnozs5Jz0JuZlZyD3sys5Bz0ZmYl56A3Mys5B72ZWck56M3MSs5Bb2ZWcg56M7OSc9CbmZWcg97MrOQc9GZmJeegNzMrOQe9zWjV5lXc+eSdjP9wvNelmNksOehtRnuP7+XhPQ+z/HeXO/DNCspBbw1NTE5w+sxpB36H5PlTU5ba8lh/HmvqpdIEfTufWP+T1DYV+Jt3b2bNV9f0upzSyPOnpiy15bH+PNbUS6UJ+nY+sf4nqW1ocIh5c+Zxx9V3sP2W7b0up1Ty/KkpS22dqr+VTlee27TbMgW9pJsk7Zd0QNL6GuvnStqern9W0nDFurvT5fsl3di+0s/VzifW/yTvmgr4X1v1axz8dwfZ9C83sej8Rb0uq5Ty/KkpS23trr8dna48t2m3NAx6SYPAJuBmYAVwq6QVVZvdDrwREVcADwIPpPddAawB3g/cBPz39PE6qp1PbL//k6xctNIB30V5/tSUpbZO1N9qpyvPbdotioiZN5CuBb4YETemt+8GiIjfqthmZ7rNLklzgOPAQmB95baV29X7fSMjIzE2Ntb8H3KfuGjiXzMUyxFCEovOX8SyC5cxNDjU1GM9c/iZdx+3xcey/rR7fDcXzL0g0//MM4efye3/WZbaOll/5Wux0oVzL2TlopUt1Z1HKy69gC/84vtndV9JuyNipNa6ORnuvxh4peL2UeAD9baJiDOS3gQWpMu/U3XfxTUKXAesA1i6dGmGkmobHBhk4OxAW57YIv6TVGomaKz9Tk6c5K2Jtzh+8njD/6Hzh87P7XOVpbZu1V/9mpxJntu0F7IEvWosq/4YUG+bLPclIrYAWyDp0Weo6RwrF63kuiUD3Hv9r7Y8vLBq851ct+Q67r3+3sIOVei+6xiaHOLlNwa4beVt3Pvhe7nkxy/pdVl9Q/dd9871v2v4PFzb3eKakqW2ztWv+65jaHCIQQ0m7Zf5NZnnNu2+LEF/FLis4vYS4FidbY6mQzcXAq9nvG9b7Pnsnlw+Vi9NTE4A8PCeh9m6d6sDv0emnofNuzez78Q+nvlM7eEIO1fSgSt2pysPsgT9c8CVki4HXiXZufqpqm12AGuBXcAtwNMREZJ2AF+R9NvApcCVwF+2q3jLxkHTW9U9UsuuLJ2uXmsY9OmY+13ATmAQ+IOI2CfpfmAsInYADwNflnSApCe/Jr3vPkmPAy8BZ4Bfj4jJDv0tVoeDpjdmN+Rg1n5ZevRExFPAU1XL/mvF9dPAJ+vcdyOwsYUabZYcNL2TlyGHVZtXce2Saz1k1+dKc2SsTVe0+e+tHAGZx1NW7Pnsnobt3o2683KUdx6fo37ScB59t812Hr01L0+9Pd0nhgaHGFDzs4RauW8vdaNu3ffuxLdetlFRn6MimWkevXv0LSpyTyUvvb0prRwBWdRTVnSz7l4f5V3U56gMHPQtyltYNiuPL75WAqnXYTZb3ag7L6cCaOffWuSOVjc56Nsgj2HZrDwFZCuBlJcwa1Yn687bSena+bcWvaPVLZlm3RRZN8ehizxfPU9TMCfPTnLrVbfy+Q9/PnMgFXWGUafrzsvsH+jc3+oDAxsrfdDvPb6Xl0681JV/gDyFZVZ5CciVi1ay9/heACZjkj9+8Y/Zvm97pucsT2HWjG7UnZcDjrrxtxa5o9VppZ910+lZB1OzCfIQls1atXlVrgKy8rmaMqABPrj0g37RWk3Vr+8ivg7bpdWzV5ZGJ97xi9qbhPz09mop4qcj641+D/gs+iroOxEeeQ7LIvKL1ppR5I5WN/VF0Ds8isEvWmuWO1rZlD7oHR7F4RetWWeUPugdHmbW73zAlJlZyTnozcxKzkFvZlZyDnozs5Jz0JuZlVzuToEg6QRwuIWHuBh4rU3ltJPrao7rao7rak4Z61oWEQtrrchd0LdK0li98z30kutqjutqjutqTr/V5aEbM7OSc9CbmZVcGYN+S68LqMN1Ncd1Ncd1Naev6irdGL2ZmU1Xxh69mZlVcNCbmZVc4YJe0vskfV3S99PLi+psNylpb/qzo2L55ZKeTe+/XdJQN+tKt71A0quSHqpY9g1J+ytq/ol21NWm2q6W9IKkA5J+V9K53/nXobokLZO0O22TfZLuqFjXkTZrQ129bK+VknalNT0v6Zcr1v2hpP9X0V4rc1JXT1+Tkv5c0t9L+rOq5T1rrwZ1Nd9eEVGoH+BLwPr0+nrggTrbnayz/HFgTXr9fwD/ppt1pet/B/gK8FDFsm8AI71ssxlq+0vgWkDA14Cbu1UXMATMTa+fDxwCLu1km7Whrl62108BV6bXLwXGgfemt/8QuKVH7TVTXT19TQI3AL8I/FnV8p61V4O6mm6vtv4B3fgB9gOXpNcvAfbX2e6coE9feK8Bc9Lb1wI7u1zX1cBjwGfoXtDPurZ0+7+q2OZWYHM366rYfgFwhM4H/azrylN7pdt9ryJgOxVcs64rD6/JdP3P1wjUnrdXdV2zba/CDd0APxkR4wDpZb2P6+dJGpP0HUkfT5ctAP4+Is6kt48Ci7tVl6QB4L8B/6nOY2xNPyLe266P+22obTFJO03papultV0m6XngFZLez7GK1Z1os1bq6nl7VdR3Dcknj7+pWLwxHTp5UNLcHNTV09dkBj1vryqzaq9cfsOUpL8Aan3v3z1NPMzSiDgmaTnwtKQXgH+osV3m+aVtqOtO4KmIeKVGJo1GxKuSfhz4KvBp4I9yUFut8OxmmxERrwBXSboU+F+SnoiIH9BCm3WqLnLQXunjXAJ8GVgbEWfTxXcDx0lCdgvwm8D9vayrzptz19urjp63V62HrrGsYXvlMugj4p/XWyfpB5IuiYjx9J/mb+s8xrH08qCkbwCrSMLgvZLmpO+IS4Bjte7fobquBT4k6U6Scd0hSScjYn1EvJr+jh9K+gpwDU0EfadqIxmzX1KxXbfbrPKxjknaB3wIeKKVNutgXd+ix+0l6QLgSeDzEfGdisceT6++LWkr8Bs5qOs1evuanOmxe9pedcyqvYo4dLMDWJteXwv8afUGki6a+pgl6WLgnwEvRTKo9X+AW2a6f6fqiojRiFgaEcMk/zR/FBHrJc1J60TSjwH/CnixTXW1VFv6z/5DST+X9r5+tdb9O1WXpCWS5qXXLyJ5Lvd3uM1mXVcO2msI+J8kz9+fVK27JL0U8HG621416+r1a3ImvWyvembdXu3e0dDpH5Ixqv8NfD+9fF+6fAT4/fT6dcALJDt8XgBur7j/cpJZEQeAPyGdOdGNuqq2/wzv7vB8D7AbeB7YR9KLHuxmm9WrrWK7F0nGVB8iPaK6S8/lL6Tt8r30cl2n26yVunLQXr8C/AjYW/GzMl33dPp6eBF4FDg/J3X19DUJ/F/gBPCPJGPeN/a6vRrU1XR7+RQIZmYlV8ShGzMza4KD3sys5Bz0ZmYl56A3Mys5B72ZWck56M3MSs5Bb2ZWcv8fn8mQvpidUaUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(CCSDTF12b2_tot,MP2_tot-CCSDTF12b2_tot,'ro')\n",
    "plt.plot(CCSDTF12b2_tot,MP3_tot-CCSDTF12b2_tot,'b^')\n",
    "plt.plot(CCSDTF12b2_tot,MP4_tot-CCSDTF12b2_tot,'g>')\n",
    "\n",
    "plt.plot([-0.5,-0.1],[0.0,0.0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performance of RPA\n",
    "\n",
    "## Dimers\n",
    "\n",
    "### Table 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 347,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RPA( SCAN0 )      -4.350363816667299\n",
      "RPA( SCAN0 )+RSE  -5.360766165539584\n",
      "RPA( PBE0 )      -4.142345416671978\n",
      "RPA( PBE0 )+RSE  -5.40889445166321\n",
      "RPA( SCAN )      -4.432773911897005\n",
      "RPA( SCAN )+RSE  -6.4609842916652\n",
      "RPA( PBE )      -3.8739801026457625\n",
      "RPA( PBE )+RSE  -6.699209809368001\n",
      "RPA( HF )      -2.983009825884056\n",
      "RPA( HF )+RSE  -2.983009824555357\n",
      "RPA( oepx )      -2.7660169999999997\n",
      "RPA( oepx )+RSE  -4.470336\n",
      "RPA( oep2sc )      -3.939989\n",
      "RPA( oep2sc )+RSE  -6.372100999999998\n"
     ]
    }
   ],
   "source": [
    "# RPA energies\n",
    "for xc in ['SCAN0','PBE0','SCAN','PBE','HF']:\n",
    "    print('RPA(',xc, ')     ', np.sum(RPA_tot[methods[xc],0,:]))\n",
    "    print('RPA(',xc, ')+RSE ', np.sum(RPA_tot[methods[xc],1,:]))\n",
    "for xc in ['oepx','oep2sc']:\n",
    "    print('RPA(',xc, ')     ', np.sum(RPA_oep[oep_m[xc],0,:]))\n",
    "    print('RPA(',xc, ')+RSE ', np.sum(RPA_oep[oep_m[xc],1,:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 348,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RPAc( SCAN0 )      -9.741077975881769\n",
      "RPAc( PBE0 )      -9.766641242326015\n",
      "RPAc( SCAN )      -11.169213800958108\n",
      "RPAc( PBE )      -11.56452199536839\n",
      "RPAc( HF )      -6.814887207178057\n"
     ]
    }
   ],
   "source": [
    "# RPAc energies, use CBS extrapolated value\n",
    "for xc in ['SCAN0','PBE0','SCAN','PBE','HF']:\n",
    "    print('RPAc(',xc, ')     ', np.sum(RPA_val[methods[xc],2,:]))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 349,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HF( SCAN0 )      5.390714159214471\n",
      "HF( PBE0 )      5.624295825654036\n",
      "HF( SCAN )      6.736439889061102\n",
      "HF( PBE )      7.690541892722627\n",
      "HF( HF )      3.8318773812940012\n"
     ]
    }
   ],
   "source": [
    "# HF/EXX energies, use AVQZ value\n",
    "for xc in ['SCAN0','PBE0','SCAN','PBE','HF']:\n",
    "    print('HF(',xc, ')     ', np.sum(EXX_val[methods[xc],1,:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 350,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HF+RSE( SCAN0 )      4.380311810342186\n",
      "HF+RSE( PBE0 )      4.357746790662805\n",
      "HF+RSE( SCAN )      4.708229509292906\n",
      "HF+RSE( PBE )      4.865312186000388\n",
      "HF+RSE( HF )      3.831877382622701\n"
     ]
    }
   ],
   "source": [
    "# HF/EXX+RSE energies, AVQZ value\n",
    "for xc in ['SCAN0','PBE0','SCAN','PBE','HF']:\n",
    "    print('HF+RSE(',xc, ')     ', np.sum(EXX_val[methods[xc],1,:])+np.sum(RSE_val[methods[xc],1,:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trimers\n",
    "\n",
    "### Table 7\n",
    "\n",
    "Values used to produce Table 7. RPA3_tot has indices [method, RSE or not, fragment]. The \\_val arrays have indices [xc, basis set, fragment, subcluster]. The subcluster needs to be set to 0 to obtain the non-additive energy of the fragment, the other indices contain interaction energies of parts of the fragment (dimers, ...).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 351,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method   PBE | PBE0 | SCAN | SCAN0 | HF  \n",
      "RPA+RSE  0.9412 0.7406 0.4070 0.3988 0.5324\n",
      "RPA      -2.0088 -0.2112 0.5947 0.6409 0.5324\n",
      "DFT/HF   7.4506 4.3924 -0.5362 -0.1100 -0.2728\n",
      "RPAc     2.3397 1.4555 0.6395 0.6553 0.8052\n",
      "HF+RSE   -1.3985 -0.7149 -0.2325 -0.2565 -0.2728\n",
      "HF       -4.3485 -1.6668 -0.0448 -0.0145 -0.2728\n"
     ]
    }
   ],
   "source": [
    "print('Method   PBE | PBE0 | SCAN | SCAN0 | HF  ')\n",
    "print('RPA+RSE ', \" \".join(\"{:.4f}\".format(np.sum(RPA3_tot[ord,1,:])) for ord in range(0,5)) )\n",
    "print('RPA     ', \" \".join(\"{:.4f}\".format(np.sum(RPA3_tot[ord,0,:])) for ord in range(0,5)) )\n",
    "print('DFT/HF  ', \" \".join(\"{:.4f}\".format(np.sum(DFT3_val[ord,1,:,0])) for ord in range(0,5)) )\n",
    "print('RPAc    ', \" \".join(\"{:.4f}\".format(np.sum(RPA3_val[ord,2,:,0])) for ord in range(0,5)) )\n",
    "print('HF+RSE  ', \" \".join(\"{:.4f}\".format(np.sum(EXX3_val[ord,1,:,0])+np.sum(RSE3_val[ord,1,:,0])) for ord in range(0,5)) )\n",
    "print('HF      ', \" \".join(\"{:.4f}\".format(np.sum(EXX3_val[ord,1,:,0])) for ord in range(0,5)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table 8\n",
    "\n",
    "Which gives the RPA contributions divided into groups according to number of hydrogen bonds in fragment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 369,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method   2b |3b 1hb | 3b 0hb | 4b 2hb | 4b 1hb | 4b 0hb | Total  \n",
      "HF+CABS  3.827 -1.483 1.210 0.164 0.471 -0.083 4.107197891504827\n",
      "CC-F12 tot -6.310 -0.398 1.434 0.056 0.552 -0.051 -4.718095521926804\n",
      "RPA( PBE ) 0 RSE     & $ -3.874$ & $-1.849$ & $-0.160$ & $1.152$ & $1.551$ & $0.374 $ & $  -2.8057869236714215 $\\\\\n",
      "RPA( PBE ) 1 RSE     & $ -6.699$ & $-0.470$ & $1.411$ & $0.132$ & $0.705$ & $0.017 $ & $  -4.904339191808855 $\\\\\n",
      "RPA( PBE0 ) 0 RSE     & $ -4.142$ & $-1.083$ & $0.872$ & $0.416$ & $0.847$ & $0.055 $ & $  -3.0359686858161967 $\\\\\n",
      "RPA( PBE0 ) 1 RSE     & $ -5.409$ & $-0.608$ & $1.349$ & $0.090$ & $0.585$ & $-0.044 $ & $  -4.037806078802072 $\\\\\n",
      "RPA( SCAN ) 0 RSE     & $ -4.433$ & $-0.721$ & $1.316$ & $-0.113$ & $0.539$ & $-0.071 $ & $  -3.483615204389226 $\\\\\n",
      "RPA( SCAN ) 1 RSE     & $ -6.461$ & $-0.642$ & $1.049$ & $0.208$ & $0.725$ & $0.017 $ & $  -5.10396193430692 $\\\\\n",
      "RPA( SCAN0 ) 0 RSE     & $ -4.350$ & $-0.668$ & $1.308$ & $-0.007$ & $0.561$ & $-0.052 $ & $  -3.2077625180458496 $\\\\\n",
      "RPA( SCAN0 ) 1 RSE     & $ -5.361$ & $-0.729$ & $1.127$ & $0.162$ & $0.634$ & $-0.026 $ & $  -4.1923322854982255 $\\\\\n",
      "RPA( HF ) 0 RSE     & $ -2.983$ & $-0.770$ & $1.303$ & $0.000$ & $0.000$ & $0.000 $ & $  -2.4505968305070986 $\\\\\n",
      "RPA( HF ) 1 RSE     & $ -2.983$ & $-0.770$ & $1.303$ & $0.000$ & $0.000$ & $0.000 $ & $  -2.4505968528475357 $\\\\\n"
     ]
    }
   ],
   "source": [
    "#we print the data stored in the mbe arrays\n",
    "print('Method   2b |3b 1hb | 3b 0hb | 4b 2hb | 4b 1hb | 4b 0hb | Total  ')\n",
    "print('HF+CABS ', \" \".join(\"{:.3f}\".format(CC_mbeQ5[0,ord]) for ord in range(1,7)), np.sum(CC_mbeQ5[0,:]) )\n",
    "print('CC-F12 tot', \" \".join(\"{:.3f}\".format(CC_mbeQ5[3,ord]) for ord in range(1,7)), np.sum(CC_mbeQ5[3,:]) )\n",
    "for xc in ['PBE','PBE0','SCAN','SCAN0','HF']:\n",
    "    for rse in range(2):\n",
    "        print('RPA(',xc,')',rse,'RSE     & $', \"$ & $\".join(\"{:.3f}\".format((RPA_mbe[methods[xc],rse,ord])) for ord in range(1,7)) ,'$ & $ ',np.sum(RPA_mbe[methods[xc],rse,:]) ,'$\\\\\\\\')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tetramers and n>4-body clusters\n",
    "\n",
    "### Table 9\n",
    "\n",
    "First four-body contributions, then 5- and higher contributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 352,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method   PBE | PBE0 | SCAN | SCAN0 \n",
      "RPA+RSE  0.8537 0.6305 0.9501 0.7696\n",
      "RPA      3.0770 1.3176 0.3545 0.5017\n",
      "DFT/HF   -3.4636 -1.6772 1.9418 1.3465\n"
     ]
    }
   ],
   "source": [
    "print('Method   PBE | PBE0 | SCAN | SCAN0 ')\n",
    "print('RPA+RSE ', \" \".join(\"{:.4f}\".format(np.sum(RPA4_tot[ord,1,:])) for ord in range(0,4)) )\n",
    "print('RPA     ', \" \".join(\"{:.4f}\".format(np.sum(RPA4_tot[ord,0,:])) for ord in range(0,4)) )\n",
    "print('DFT/HF  ', \" \".join(\"{:.4f}\".format(np.sum(DFT4_val[ord,1,:,0])) for ord in range(0,4)) )\n",
    "\n",
    "#print(np.sum(RPA_mbe[0,0,4:7]))\n",
    "#print(RPA_mbe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 353,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Method   PBE | PBE0 | SCAN | SCAN0 \n",
      "RPA+RSE  -0.0857 0.0078 -0.0060 0.0123\n",
      "RPA      -0.9242 -0.2640 0.1036 0.0378\n",
      "[-3.73 -4.99]\n",
      "[[ 0.         -3.8739801  -1.0238456  -0.9849528   1.15174537  1.55075073\n",
      "   0.37449548]\n",
      " [ 0.         -6.69920981  0.35525191  0.58596007  0.13191714  0.70520967\n",
      "   0.01653183]]\n"
     ]
    }
   ],
   "source": [
    "print('Method   PBE | PBE0 | SCAN | SCAN0 ')\n",
    "print('RPA+RSE ', \" \".join(\"{:.4f}\".format(RPA_sup[ord,1]-np.sum(RPA_mbe[ord,1,:])) for ord in range(0,4)) )\n",
    "print('RPA     ', \" \".join(\"{:.4f}\".format(RPA_sup[ord,0]-np.sum(RPA_mbe[ord,0,:])) for ord in range(0,4)) )\n",
    "\n",
    "print(RPA_sup[0,:])\n",
    "print(RPA_mbe[0,:,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Incremental approach\n",
    "\n",
    "Here are the data obtained by using Eq. 11 from the manuscript. CCSD(T) interaction energies for small n-body contributions complemented with RPA terms for the higher-order contributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 366,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.         -5.40889445 -0.60817431  1.3487677   0.08986102  0.58456139\n",
      " -0.04392743]\n",
      "CC/RPA PBE        [-3.73       -6.16644704 -4.71567436 -3.12232334 -4.2177322  -5.21684588\n",
      " -5.6423086 ]\n",
      "CC/RPA PBE+RSE    [-4.99       -4.60121734 -4.52954217 -4.50710402 -4.58268464 -4.73625726\n",
      " -4.80375633]\n",
      "CC/RPA PBE0       [-3.3        -5.46808173 -4.78288125 -4.22151392 -4.58070137 -4.87618137\n",
      " -4.98212684]\n",
      "CC/RPA PBE0+RSE   [-4.03       -4.9315327  -4.72158391 -4.63680079 -4.67032529 -4.70324964\n",
      " -4.71028944]\n",
      "CC/RPA SCAN       [-3.38       -5.25765324 -4.9349942  -4.81699786 -4.64781495 -4.63487774\n",
      " -4.61448032]\n",
      "CC/RPA SCAN+RSE   [-5.11       -4.95944286 -4.71539922 -4.33106995 -4.48270528 -4.65641126\n",
      " -4.72413359]\n",
      "CC/RPA SCAN0      [-3.17       -5.13006333 -4.86074033 -4.735595   -4.67186205 -4.68162449\n",
      " -4.680333  ]\n",
      "CC/RPA SCAN0+RSE  [-4.18       -5.12966098 -4.79928427 -4.49315776 -4.5987332  -4.68106188\n",
      " -4.70576324]\n",
      "CC/MP2            [-4.31       -5.42880439 -4.87121096 -4.64612328 -4.70664676 -4.71830458\n",
      " -4.7018568 ]\n",
      "CC/HF             [ 4.11       -6.02767583 -4.94337364 -4.72011808 -4.82826619 -4.74779783\n",
      " -4.71529341]\n",
      "CC-F12 tot        0.000 -6.310 -6.709 -5.275 -5.219 -4.667 -4.718\n"
     ]
    }
   ],
   "source": [
    "#array to store RPA corrected MBE\n",
    "RPA_cortd=np.zeros((5,2,7))\n",
    "#array to store % errors\n",
    "RPA_cortd_pct=np.zeros((5,2,7))\n",
    "#array to store MP2 corrected MBE\n",
    "MP_cortd=np.zeros((7))\n",
    "HF_cortd=np.zeros((7))\n",
    "MP_cortd_pct=np.zeros((7))\n",
    "HF_cortd_pct=np.zeros((7))\n",
    "#start with the RPA supermolecular calculation\n",
    "RPA_act=RPA_sup\n",
    "CC_act=0.\n",
    "RPA_cortd[:,:,0]=RPA_act\n",
    "MP_cortd[0]=-4.31\n",
    "HF_cortd[0]=4.11\n",
    "\n",
    "#calculate correction scheme based on different inputs for the different orders (2b, 3b 1hb, 3b 0hb, ...)\n",
    "for order in range(1,7):\n",
    "    for xc in range(5):\n",
    "        for rse in range(2):\n",
    "            #take the previous energy as a reference and subtract the RPA part and add the CC part\n",
    "            RPA_cortd[xc,rse,order]=RPA_cortd[xc,rse,order-1]+CC_mbeQ5[3,order]-RPA_mbe[xc,rse,order]\n",
    "            #-4.71 is our reference value\n",
    "            RPA_cortd_pct[xc,rse,order]=(RPA_cortd[xc,rse,order]+4.71)/(-4.71)*100\n",
    "    #calculate MP2- and HF-based corrections\n",
    "    MP_cortd[order]=MP_cortd[order-1]+CC_mbeQ5[3,order]-MP_mbe[0,order]\n",
    "    HF_cortd[order]=HF_cortd[order-1]+CC_mbeQ5[3,order]-CC_mbeQ5[0,order]\n",
    "    MP_cortd_pct[order]=(MP_cortd[order]+4.71)/(-4.71)*100\n",
    "    HF_cortd_pct[order]=(HF_cortd[order]+4.71)/(-4.71)*100\n",
    "\n",
    "    \n",
    "print(RPA_mbe[1,1,:])\n",
    "    \n",
    "#print the results\n",
    "print('CC/RPA PBE       ', RPA_cortd[0,0,:])\n",
    "print('CC/RPA PBE+RSE   ', RPA_cortd[0,1,:])\n",
    "print('CC/RPA PBE0      ', RPA_cortd[1,0,:])\n",
    "print('CC/RPA PBE0+RSE  ', RPA_cortd[1,1,:])\n",
    "print('CC/RPA SCAN      ', RPA_cortd[2,0,:])\n",
    "print('CC/RPA SCAN+RSE  ', RPA_cortd[2,1,:])\n",
    "print('CC/RPA SCAN0     ', RPA_cortd[3,0,:])\n",
    "print('CC/RPA SCAN0+RSE ', RPA_cortd[3,1,:])\n",
    "#print('CC/RPA HF        ', RPA_cortd[4,0,:]) #we don't have the full cluster data\n",
    "print('CC/MP2           ', MP_cortd)\n",
    "print('CC/HF            ', HF_cortd)\n",
    "\n",
    "print('CC-F12 tot       ', \" \".join(\"{:.3f}\".format(np.sum(CC_mbeQ5[3,0:ord+1])) for ord in range(7)) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 370,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-4.721583905523378\n"
     ]
    }
   ],
   "source": [
    "#So for RPA(PBE0)+RSE the difference after adding 2b and 3b 1hb is \n",
    "print(RPA_cortd[1,1,2])"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}