#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 20 10:49:24 2020

@author: bagger
"""
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.patches import Ellipse
from scipy.stats import norm, chi2
import numpy.linalg as linalg
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import NullFormatter
from scipy.interpolate import splrep, splev
from scipy.optimize import minimize
from matplotlib.ticker import MultipleLocator
import matplotlib.patches as mpatches
import csv
from numpy import genfromtxt
import pandas as pd
#from pandas.tools.plotting import scatter_matrix
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from matplotlib.artist import setp
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import seaborn as sns
from os import path
import ase.db
from ase.db import connect
from scipy.optimize import fsolve

########################################################
### Defining covariance function and ellipse.
########################################################
def cov_ellipse(cov, q=None, nsig=None, **kwargs):

    if q is not None:
        q = np.asarray(q)
    elif nsig is not None:
        q = 2 * norm.cdf(nsig) - 1
    else:
        raise ValueError('One of `q` and `nsig` should be specified.')
    r2 = chi2.ppf(q, 2)
    val, vec = np.linalg.eigh(cov)
    width, height = 2 * np.sqrt(val[:, None] * r2)
    rotation = np.degrees(np.arctan2(*vec[::-1, 0]))

    return width, height, rotation
    
def plot_point_cov(points, nstd=2, ax=None, **kwargs):
    pos = points.mean(axis=0)
    cov = np.cov(points, rowvar=False)
    return cov, pos

#############################################################
#### load in data
##################################################################
with open('../DFT_data/Molecules//H2_BEEF-vdw_energy.txt', 'r') as f:
    EH2mBEEF = [line.rstrip('\n') for line in f]
    data = [float(i) for i in EH2mBEEF]
    EH2=np.asarray(data)
with open('../DFT_data/Molecules//CO2_BEEF-vdw_energy.txt', 'r') as f:
    EH2mBEEF = [line.rstrip('\n') for line in f]
    data = [float(i) for i in EH2mBEEF]
    ECO2=np.asarray(data)

#############################################################
#### load in Metal data
##################################################################
dict_Metals={}
d=ase.db.connect('Metals_large.db')
names=['Ag','Au','Co','Cu','Ni','Pd','Pt','Rh','Ru','Cd','Fe','In','Ir','Pb','Sn','Tl']
# Key words for database:
# Metal 'Cu','Ag','Au','Ir','Pt','Pd','Rh','Sn','In','Pb','Cd','Tl','Hg','Ru', 'Zn', 'Ni', 'Fe', 'Ti'
# ads: 'COOH', 'H','CO', 'Clean', 'HCOO', 'CH3O', 'CH2OH'
for metal in ['Cu','Ag','Au','Ir','Pt','Pd','Rh','Sn','In','Pb','Cd','Tl','Ru', 'Zn', 'Ni', 'Fe', 'Ti']:
    dict_Metals[metal+'_slab']=d.get(Metal=metal, ads='Clean').data.Beef
    dict_Metals[metal+'_COOH']=d.get(Metal=metal, ads='COOH').data.Beef
    dict_Metals[metal+'_H']=d.get(Metal=metal, ads='H').data.Beef
    dict_Metals[metal+'_HCOO']=d.get(Metal=metal, ads='HCOO').data.Beef 



#############################################################
#### load in molecule metal data
##################################################################
conMol=ase.db.connect('Molecules.db')
ECO2_metal= conMol.get(Molecule='CO2').data.Beef
EH2_metal= conMol.get(Molecule='H2').data.Beef


#############################################################
#### load in molecule new gpaw version  data
##################################################################
MolBEEF=ase.db.connect('molecules_BEEF_vdW.db')
ECO2=np.asarray(MolBEEF.get(mol='CO2').data.Beef)
EH2=np.asarray(MolBEEF.get(mol='H2').data.Beef)

#############################################################
#### load in BEEF MNC data 
##################################################################
names_MNC=['Co','Fe','Ir','Mn','Mo','Pt','Rh','Ru']
dict_MNC={}     
dMNC=ase.db.connect('MNC_structures.db')
# Key words for database:
# metal: 'Co','Fe','Ir','Mn','Mo','Pt','Rh','Ru'
# ads: 'slab', 'CO', 'H', 'N2'
# note that there are several structures for each, as spin 0,1,2,3,4 has been tested.
for metal in names_MNC:
    for ads in ['slab','N2','CO','H','COOH','HCOO']:
        energies = []
        ids=[]
        for row in dMNC.select(metal=metal, ads=ads, sort='energy'):
            energies.append(row.energy)
            ids.append(row.id)
        dict_MNC[metal+'_'+ads]=np.asarray(dMNC.get(id=ids[0]).data.Beef)

#############################################################
#### load in BEEF XHg4 data 
##################################################################
names_XHg=['Cu','Pt','Pd']
dict_XHg={}
d1=ase.db.connect('XHg4.db')
for metal in names_XHg:
    for ads1 in ['slab','H','COOH','HCOO']:
        dict_XHg[metal+'_'+ads1]=d1.get(metal=metal+'Hg4', ads=ads1).data.Beef

#############################################################
#### load in BEEF Au111 + single Pt atom data 
##################################################################
dict_AuPt={}
d2=ase.db.connect('Au111_XPt.db')
dict_AuPt['Au_1Pt_slab']=d2.get(metal='Au_1Pt.gpw').data.Beef
dict_AuPt['Au_1Pt_H']=d2.get(metal='Au_1Pt_H.gpw').data.Beef
dict_AuPt['Au_1Pt_COOH']=d2.get(metal='Au_1Pt_COOH.gpw').data.Beef
dict_AuPt['Au_1Pt_HCOO_1']=d2.get(metal='Au_1Pt_HCOO_1.gpw').data.Beef
dict_AuPt['Au_1Pt_HCOO_2']=d2.get(metal='Au_1Pt_HCOO_2.gpw').data.Beef

# Creation of references:
metal='Au'
val_Au_H=(dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-0.5*np.asarray(EH2_metal)).mean()
val_Au_COOH=(dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-ECO2_metal-0.5*np.asarray(EH2_metal)).mean()
val_Au_HCOO=(dict_Metals[metal+'_HCOO']-dict_Metals[metal+'_slab']-ECO2_metal-0.5*np.asarray(EH2_metal)).mean()

metal='Fe'
val_FeNC_H=(dict_MNC[metal+'_H']-dict_MNC[metal+'_slab']-0.5*np.asarray(EH2)).mean()
val_FeNC_COOH=(dict_MNC[metal+'_COOH']-dict_MNC[metal+'_slab']-ECO2-0.5*np.asarray(EH2)).mean()
val_FeNC_HCOO=(dict_MNC[metal+'_HCOO']-dict_MNC[metal+'_slab']-ECO2-0.5*np.asarray(EH2)).mean()

metal='Pt'
val_PtHg4_H=(dict_XHg[metal+'_H']-dict_XHg[metal+'_slab']-0.5*np.asarray(EH2)).mean()
val_PtHg4_COOH=(dict_XHg[metal+'_COOH']-dict_XHg[metal+'_slab']-ECO2-0.5*np.asarray(EH2)).mean()
val_PtHg4_HCOO=(dict_XHg[metal+'_HCOO']-dict_XHg[metal+'_slab']-ECO2-0.5*np.asarray(EH2)).mean()

fonttype='sansserif'
L=3; size1=28; size2=20; size3=24;A1=0.5;
s1=15; s2=15; s3=15; sdot=150; sdot2=150;
GiveColor='b'

##############################################################################
# Figure COOH vs H 
##############################################################################
fig = plt.figure(figsize=(7,7));
plt.locator_params(axis='x',nbins=5);plt.grid(True)
plt.xticks(fontsize = size3); plt.yticks(fontsize = size3)
ax=fig.gca()
for metal in ['Ag','Cu','Ni','Pd','Pt','Rh','Cd','Fe','In','Ir','Pb','Sn','Tl']:
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    if metal=='Fe':
        p1=plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P')    
    else:
        plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P') 
    if metal=='Sn':
        ax.text(np.asarray(Exval).mean()-0.15, np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)
    elif metal=='Pt':
        ax.annotate('Pt', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.2, np.asarray(Eyval).mean()),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
    else:
        ax.text(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), metal,color='black',alpha=A1, fontsize=size2)

for metal in ['Au']:
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-0.5*np.asarray(EH2_metal)
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-ECO2_metal-0.5*np.asarray(EH2_metal)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P') 
    ax.annotate('Metals', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.4, np.asarray(Eyval).mean()+0.1),color='blue',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
    ax.text(np.asarray(Exval).mean()-0.00, np.asarray(Eyval).mean(), metal,color='black',alpha=A1, fontsize=size2)

MetalmeanCOOH=[]; PorpyrinmeanCOOH=[]; MetalmeanH=[]; PorpyrinmeanH=[]
for metal in ['Ag','Cu','Ni','Pd','Pt','Rh','Fe','Ir','Au']:        
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    cov1, pos1=plot_point_cov(np.transpose([Exval,Eyval]))
    MetalmeanH.append(pos1[0])  
    MetalmeanCOOH.append(pos1[1])
s=1
def f(i):
    """Fixed slope 1-deg polynomial residuals"""
    return ((np.asarray(MetalmeanCOOH) - (s*np.asarray(MetalmeanH) + i))**2).sum()
metalnumber=fsolve(f, x0=1)
plt.plot(np.asarray(MetalmeanH), np.asarray(MetalmeanH)+fsolve(f, x0=1), '--b')
plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)
angle = 45
th2 = plt.text(0-0.05, fsolve(f, x0=1)+0.02-0.05, '$\Delta$E$_{^*\!\mathrm{COOH}}$=$\Delta$E$_{^*\!\mathrm{H}}$+0.29', fontsize=16,
               rotation=angle, rotation_mode='anchor')

for metal in ['Co','Ir','Mn','Rh']:
    Exval=dict_MNC[metal+'_H']-dict_MNC[metal+'_slab']-(dict_MNC['Fe_H']-dict_MNC['Fe_slab'])+val_FeNC_H
    Eyval=dict_MNC[metal+'_COOH']-dict_MNC[metal+'_slab']-(dict_MNC['Fe_COOH']-dict_MNC['Fe_slab'])+val_FeNC_COOH
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='o', color='k',s=sdot) 
    if metal=='Mn':
        ax.annotate('MNC', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()+0.2, np.asarray(Eyval).mean()-0.1),color='black',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
    if metal=='Rh' or metal=='Co':
        ax.text(np.asarray(Exval).mean()-0.08, np.asarray(Eyval).mean()+0.02, metal,color='black',alpha=A1, fontsize=size2)
    else:
        ax.text(np.asarray(Exval).mean()-0.05, np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)
    
for metal in ['Fe']:
    Exval=dict_MNC[metal+'_H']-dict_MNC[metal+'_slab']-0.5*np.asarray(EH2)
    Eyval=dict_MNC[metal+'_COOH']-dict_MNC[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='o', color='k',s=sdot) 
    ax.text(np.asarray(Exval).mean()-0.05, np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)

for metal in ['Pd']:
    Exval=dict_XHg[metal+'_H']-dict_XHg[metal+'_slab']-(dict_XHg['Pt_H']-dict_XHg['Pt_slab'])+val_PtHg4_H
    Eyval=dict_XHg[metal+'_COOH']-dict_XHg[metal+'_slab']-(dict_XHg['Pt_COOH']-dict_XHg['Pt_slab'])+val_PtHg4_COOH
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='s', color='m',s=sdot) 
    ax.annotate('PdHg$_{4}$', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()+0.05, np.asarray(Eyval).mean()-0.01),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))

for metal in ['Pt']:
    Exval=dict_XHg[metal+'_H']-dict_XHg[metal+'_slab']-0.5*np.asarray(EH2)
    Eyval=dict_XHg[metal+'_COOH']-dict_XHg[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='s', color='g',s=sdot) 
    ax.annotate('PtHg$_{4}$', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()+0.3, np.asarray(Eyval).mean()-0.01),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))

for metal in ['Au_1Pt']:
    Exval=dict_AuPt[metal+'_H']-dict_AuPt[metal+'_slab']-0.5*np.asarray(EH2)
    Eyval=dict_AuPt[metal+'_COOH']-dict_AuPt[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='^', color='c',s=sdot) 
    ax.annotate('Pt$_1$Au(111)', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()+0.3, np.asarray(Eyval).mean()+0.005),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))


plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)
plt.ylabel('$\Delta E_{~^*\!\mathrm{COOH}}$ [eV]',fontsize=size1)
plt.plot([-2,3],[-2,3],'k--',lw=2)
plt.plot([-0.08429,-0.08429],[-5.25,3.5],'k',lw=2)
plt.plot([-1.0,3.0],[0.12,0.12],'k',lw=2)
ax.annotate('$\Delta E^{0}_{FAOR}$', xy=(0.75,0.12), xytext=(0.72, 0.0),color='black',alpha=1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
ax.annotate('$\Delta E_{H_{UPD}}$', xy=(-0.08429,0.8), xytext=(-0.45, 0.9),color='black',alpha=1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
plt.xlim([-0.5,1.0])
plt.ylim([-0.5,1.0])
fig.text(0.0, 0.9, 'a)', ha='center',weight = 'bold', fontsize=size3+10)
plt.savefig('COOHvsH.png',dpi=400, bbox_inches='tight')


##############################################################################
# Figure COOH vs H as a function of potential U = - 0.12 V vs RHE
##############################################################################
fig = plt.figure(figsize=(7,7));
plt.locator_params(axis='x',nbins=5);plt.grid(True)
plt.xticks(fontsize = size3); plt.yticks(fontsize = size3)
ax=fig.gca()

for metal in ['Pt']:
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c='b',s=sdot2, marker='P')    

    ax.annotate('Pt', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.1, np.asarray(Eyval).mean()),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
MetalmeanCOOH=[]; PorpyrinmeanCOOH=[]; MetalmeanH=[]; PorpyrinmeanH=[]
for metal in ['Ag','Cu','Ni','Pd','Pt','Rh','Fe','Ir','Au']:        
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    cov1, pos1=plot_point_cov(np.transpose([Exval,Eyval]))
    MetalmeanH.append(pos1[0])  
    MetalmeanCOOH.append(pos1[1])
s=1
def f(i):
    """Fixed slope 1-deg polynomial residuals"""
    return ((np.asarray(MetalmeanCOOH) - (s*np.asarray(MetalmeanH) + i))**2).sum()
from scipy.optimize import fsolve
metalnumber=fsolve(f, x0=1)
plt.plot(np.asarray(MetalmeanH), np.asarray(MetalmeanH)+fsolve(f, x0=1), '--b')
plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)
V=-0.12
plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)
plt.ylabel('$\Delta E_{~^*\!\mathrm{COOH}}$ [eV]',fontsize=size1)
plt.plot([-2,3],[-2,3],'k--',lw=2)
plt.plot([0.12+(0.12+V)-0.29,3.0],[0.12+(0.12+V),0.12+(0.12+V)],'r',lw=2)
plt.plot([-0.08429-V,3.0],[0.12-(0.12+V),0.12-(0.12+V)],'r',lw=2)
plt.plot([-0.08429-V,-0.08429-V],[-5.25,0.29+-0.08429-V],'r',lw=2)
plt.annotate(s='', xy=(-0.08429+0.12-0.3-0.1,-0.1), xytext=(-0.08429-V,-0.1),color='red', arrowprops=dict(arrowstyle='<->',color='red'),size=16)
ax.text(-0.3, 0.0-0.07, '$\eta_{^*\!\mathrm{H}}$',color='red',alpha=1, fontsize=size2)
plt.annotate(s='', xy=(0.4,0.12+(0.12+V)), xytext=(0.4,0.4+0.29), arrowprops=dict(arrowstyle='<->', color='red'),size=16)
ax.text(0.42, 0.67, '$\eta_{^*\!\mathrm{COOH}}$',color='red',alpha=1, fontsize=size2)
ax.text(0.45,0.12+(0.12+V)+0.05,'$\Delta E^{0}_{FAOR}$+%s eV' % np.round(((0.12+V)),2),color='black',alpha=1, fontsize=size2-2)
ax.text(0.45,0.12-(0.12+V)-0.1,'$\Delta E^{0}_{FAOR}$%s eV' % np.round(-(0.12+V),2),color='black',alpha=1, fontsize=size2-2)
ax.annotate('$\Delta E_{H_{UPD}}$+%s eV' % np.round(-V,2), xy=(-0.08429-V,-0.45), xytext=(0.1, -0.38),color='black',alpha=1, fontsize=size2-2,arrowprops=dict(arrowstyle="-"))
plt.annotate('@ %s V$_{\mathrm{RHE}}$' % np.round(V,2), xy=(0.03, 0.9), xycoords='axes fraction', fontsize=size1-4,
    bbox=dict(boxstyle='round',facecolor='white', alpha=1.0),
    horizontalalignment='left', verticalalignment='bottom')
xval=-0.08429-V; yval=-0.08429+0.29-V;
plt.scatter(xval, yval, c=GiveColor,s=sdot2, marker='o') 
ax.annotate('$^*\!\mathrm{COOH}$ limit', xy=(xval, yval), xytext=(-0.4, 0.45),color='black',alpha=1, fontsize=16,arrowprops=dict(arrowstyle="-"))
xval=0.12+(0.12+V)-0.29; yval=0.12+(0.12+V);
plt.scatter(xval, yval, c=GiveColor,s=sdot2, marker='o') 
ax.annotate('$^*\!\mathrm{H}$ limit', xy=(xval, yval), xytext=(-0.45, 0.2),color='black',alpha=1, fontsize=16,arrowprops=dict(arrowstyle="-"))
ax.annotate('Metals', xy=(0.2,0.2+0.29), xytext=(0.2-0.15, 0.2+0.29+0.2),color='blue',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
ax.annotate('Single-site', xy=(0.8,0.8), xytext=(0.3, 0.9),color='black',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
xval=0.12; yval=0.12
plt.scatter(xval, yval, c='r',s=sdot2+100,alpha=0.5, marker='o') 
plt.xlim([-0.5,1.0])
plt.ylim([-0.5,1.0])
fig.text(0.0, 0.9, 'a)', ha='center',weight = 'bold', fontsize=size3+10)
plt.savefig('COOHvsH_potential_%s.png' % V,dpi=400, bbox_inches='tight')


##############################################################################
# Figure COOH vs H as a function of potential U = 0.0 V vs RHE
##############################################################################
fig = plt.figure(figsize=(7,7));
plt.locator_params(axis='x',nbins=5);plt.grid(True)
plt.xticks(fontsize = size3); plt.yticks(fontsize = size3)
ax=fig.gca()

for metal in ['Pt']:
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c='b',s=sdot2, marker='P')    
    ax.annotate('Pt', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.1, np.asarray(Eyval).mean()),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))

MetalmeanCOOH=[]; PorpyrinmeanCOOH=[]; MetalmeanH=[]; PorpyrinmeanH=[]
for metal in ['Ag','Cu','Ni','Pd','Pt','Rh','Fe','Ir','Au']:        
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    cov1, pos1=plot_point_cov(np.transpose([Exval,Eyval]))
    MetalmeanH.append(pos1[0])  
    MetalmeanCOOH.append(pos1[1])
s=1
def f(i):
    """Fixed slope 1-deg polynomial residuals"""
    return ((np.asarray(MetalmeanCOOH) - (s*np.asarray(MetalmeanH) + i))**2).sum()
from scipy.optimize import fsolve
metalnumber=fsolve(f, x0=1)
plt.plot(np.asarray(MetalmeanH), np.asarray(MetalmeanH)+fsolve(f, x0=1), '--b')

plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)

V=0.00
plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)
plt.ylabel('$\Delta E_{~^*\!\mathrm{COOH}}$ [eV]',fontsize=size1)
plt.plot([-2,3],[-2,3],'k--',lw=2)
plt.plot([0.12+(0.12+V)-0.29,3.0],[0.12+(0.12+V),0.12+(0.12+V)],'r',lw=2)
plt.plot([-0.08429-V,3.0],[0.12-(0.12+V),0.12-(0.12+V)],'r',lw=2)
plt.plot([-0.08429-V,-0.08429-V],[-5.25,0.29+-0.08429-V],'r',lw=2)
plt.annotate(s='', xy=(-0.08429+0.12-0.3-0.1,-0.1), xytext=(-0.08429-V,-0.1),color='red', arrowprops=dict(arrowstyle='<->',color='red'),size=16)
ax.text(-0.3, 0.0-0.07, '$\eta_{^*\!\mathrm{H}}$',color='red',alpha=1, fontsize=size2)
plt.annotate(s='', xy=(0.4,0.12+(0.12+V)), xytext=(0.4,0.4+0.29), arrowprops=dict(arrowstyle='<->', color='red'),size=16)
ax.text(0.42, 0.67, '$\eta_{^*\!\mathrm{COOH}}$',color='red',alpha=1, fontsize=size2)
ax.text(0.45,0.12+(0.12+V)+0.05,'$\Delta E^{0}_{FAOR}$+%s eV' % np.round(((0.12+V)),2),color='black',alpha=1, fontsize=size2-2)
ax.text(0.45,0.12-(0.12+V)-0.1,'$\Delta E^{0}_{FAOR}$%s eV' % np.round(-(0.12+V),2),color='black',alpha=1, fontsize=size2-2)
ax.annotate('$\Delta E_{H_{UPD}}$+%s eV' % np.round(-V,2), xy=(-0.08429-V,-0.45), xytext=(0.1, -0.38),color='black',alpha=1, fontsize=size2-2,arrowprops=dict(arrowstyle="-"))
plt.annotate('@ %s V$_{\mathrm{RHE}}$' % np.round(V,2), xy=(0.03, 0.9), xycoords='axes fraction', fontsize=size1-4,
    bbox=dict(boxstyle='round',facecolor='white', alpha=1.0),
    horizontalalignment='left', verticalalignment='bottom')
xval=0.12+(0.12+V)-0.29; yval=0.12+(0.12+V);
plt.scatter(xval, yval, c=GiveColor,s=sdot2, marker='o') 
ax.annotate('$^*\!\mathrm{COOH}$ limit', xy=(xval, yval), xytext=(-0.4, 0.45),color='black',alpha=1, fontsize=16,arrowprops=dict(arrowstyle="-"))
xval=-0.08429-V; yval=-0.08429+0.29-V;
plt.scatter(xval, yval, c=GiveColor,s=sdot2, marker='o') 
ax.annotate('$^*\!\mathrm{H}$ limit', xy=(xval, yval), xytext=(-0.45, 0.2),color='black',alpha=1, fontsize=16,arrowprops=dict(arrowstyle="-"))
ax.annotate('Metals', xy=(0.2,0.2+0.29), xytext=(0.2-0.15, 0.2+0.29+0.2),color='blue',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
ax.annotate('Single-site', xy=(0.8,0.8), xytext=(0.3, 0.9),color='black',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
left, bottom, width, height = (-0.08429-V, 0.12-(0.12+V), 2, 2*(0.12+V))
rect=mpatches.Rectangle((left,bottom),width,height, 
                        fill=True,
                        color="red",
                       #linewidth=2,
                       facecolor="red", alpha=0.5
                       )
plt.gca().add_patch(rect)
plt.text(0.2,0.1,'Active Region',color='black',alpha=0.8, fontsize=size1)
plt.xlim([-0.5,1.0])
plt.ylim([-0.5,1.0])
fig.text(0.0, 0.9, 'b)', ha='center',weight = 'bold', fontsize=size3+10)
plt.savefig('COOHvsH_potential_%s.png' % V,dpi=400, bbox_inches='tight')

##############################################################################
# Figure COOH vs H as a function of potential U = 0.17 V vs RHE
##############################################################################
fig = plt.figure(figsize=(7,7));
plt.locator_params(axis='x',nbins=5);plt.grid(True)
plt.xticks(fontsize = size3); plt.yticks(fontsize = size3)
ax=fig.gca()
for metal in ['Pt']:
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c='b',s=sdot2, marker='P')    
    ax.annotate('Pt', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.1, np.asarray(Eyval).mean()),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
MetalmeanCOOH=[]; PorpyrinmeanCOOH=[]; MetalmeanH=[]; PorpyrinmeanH=[]
for metal in ['Ag','Cu','Ni','Pd','Pt','Rh','Fe','Ir','Au']:        
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    cov1, pos1=plot_point_cov(np.transpose([Exval,Eyval]))
    MetalmeanH.append(pos1[0])  
    MetalmeanCOOH.append(pos1[1])
s=1
def f(i):
    """Fixed slope 1-deg polynomial residuals"""
    return ((np.asarray(MetalmeanCOOH) - (s*np.asarray(MetalmeanH) + i))**2).sum()
from scipy.optimize import fsolve
metalnumber=fsolve(f, x0=1)
plt.plot(np.asarray(MetalmeanH), np.asarray(MetalmeanH)+fsolve(f, x0=1), '--b')
plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)
V=0.17
plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)
plt.ylabel('$\Delta E_{~^*\!\mathrm{COOH}}$ [eV]',fontsize=size1)
plt.plot([-2,3],[-2,3],'k--',lw=2)
plt.plot([0.12+(0.12+V)-0.29,3.0],[0.12+(0.12+V),0.12+(0.12+V)],'r',lw=2)
plt.plot([-0.08429-V,3.0],[0.12-(0.12+V),0.12-(0.12+V)],'r',lw=2)
plt.plot([-0.08429-V,-0.08429-V],[-5.25,0.29+-0.08429-V],'r',lw=2)
plt.annotate(s='', xy=(-0.08429+0.12-0.3-0.1,-0.1), xytext=(-0.08429-V,-0.1),color='red', arrowprops=dict(arrowstyle='<->',color='red'),size=16)
ax.text(-0.499, 0.0-0.07, '$\eta_{^*\!\mathrm{H}}$',color='red',alpha=1, fontsize=size2)
plt.annotate(s='', xy=(0.4,0.12+(0.12+V)), xytext=(0.4,0.4+0.29), arrowprops=dict(arrowstyle='<->', color='red'),size=16)
ax.text(0.42, 0.67, '$\eta_{^*\!\mathrm{COOH}}$',color='red',alpha=1, fontsize=size2)
ax.text(0.45,0.12+(0.12+V)+0.05,'$\Delta E^{0}_{FAOR}$+%s eV' % np.round(((0.12+V)),2),color='black',alpha=1, fontsize=size2-2)
ax.text(0.45,0.12-(0.12+V)-0.1,'$\Delta E^{0}_{FAOR}$%s eV' % np.round(-(0.12+V),2),color='black',alpha=1, fontsize=size2-2)
ax.annotate('$\Delta E_{H_{UPD}}$+%s eV' % np.round(-V,2), xy=(-0.08429-V,-0.45), xytext=(0.1, -0.38),color='black',alpha=1, fontsize=size2-2,arrowprops=dict(arrowstyle="-"))
plt.annotate('@ %s V$_{\mathrm{RHE}}$' % np.round(V,2), xy=(0.03, 0.9), xycoords='axes fraction', fontsize=size1-4,
    bbox=dict(boxstyle='round',facecolor='white', alpha=1.0),
    horizontalalignment='left', verticalalignment='bottom')
xval=0.12+(0.12+V)-0.29; yval=0.12+(0.12+V);
plt.scatter(xval, yval, c=GiveColor,s=sdot2, marker='o') 
ax.annotate('$^*\!\mathrm{COOH}$ limit', xy=(xval, yval), xytext=(-0.4, 0.45),color='black',alpha=1, fontsize=16,arrowprops=dict(arrowstyle="-"))
xval=-0.08429-V; yval=-0.08429+0.29-V;
plt.scatter(xval, yval, c=GiveColor,s=sdot2, marker='o') 
ax.annotate('$^*\!\mathrm{H}$ limit', xy=(xval, yval), xytext=(-0.45, 0.2),color='black',alpha=1, fontsize=16,arrowprops=dict(arrowstyle="-"))
ax.annotate('Metals', xy=(0.2,0.2+0.29), xytext=(0.2-0.15, 0.2+0.29+0.2),color='blue',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
ax.annotate('Single-site', xy=(0.8,0.8), xytext=(0.3, 0.9),color='black',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
left, bottom, width, height = (-0.08429-V, 0.12-(0.12+V), 2, 2*(0.12+V))
rect=mpatches.Rectangle((left,bottom),width,height, 
                        fill=True,
                        color="red",
                       #linewidth=2,
                       facecolor="red", alpha=0.5
                       )
plt.gca().add_patch(rect)
plt.text(0.2,0.1,'Active Region',color='black',alpha=0.8, fontsize=size1)
plt.xlim([-0.5,1.0])
plt.ylim([-0.5,1.0])
fig.text(0.0, 0.9, 'c)', ha='center',weight = 'bold', fontsize=size3+10)
plt.savefig('COOHvsH_potential_%s.png' % V,dpi=400, bbox_inches='tight')


##############################################################################
# Figure HCOO vs H
##############################################################################
fig = plt.figure(figsize=(7,7));
plt.locator_params(axis='x',nbins=5);plt.grid(True)
plt.xticks(fontsize = size3); plt.yticks(fontsize = size3)
ax=fig.gca()
for metal in ['Ag','Cu','Ni','Pd','Pt','Rh','Cd','Fe','In','Ir','Pb','Sn','Tl']:   
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H
    Eyval=dict_Metals[metal+'_HCOO']-dict_Metals[metal+'_slab']-(dict_Metals['Au_HCOO']-dict_Metals['Au_slab'])+val_Au_HCOO
    if metal=='Fe':
       plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P')    
    else:
       plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P') 
    if metal=='Sn':
        ax.text(np.asarray(Exval).mean()-0.15, np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)
    elif metal=='Pt':
        ax.annotate('Pt', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.15, np.asarray(Eyval).mean()),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
    elif metal=='Pd':
        ax.annotate('Pd', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.15, np.asarray(Eyval).mean()-0.1),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
    else:
        ax.text(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), metal,color='black',alpha=A1, fontsize=size2)

for metal in ['Au']:
    Exval=dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-0.5*np.asarray(EH2_metal)
    Eyval=dict_Metals[metal+'_HCOO']-dict_Metals[metal+'_slab']-ECO2_metal-0.5*np.asarray(EH2_metal)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P') 
    ax.annotate('Metals', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.4, np.asarray(Eyval).mean()+0.1),color='blue',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
    ax.text(np.asarray(Exval).mean()-0.00, np.asarray(Eyval).mean(), metal,color='black',alpha=A1, fontsize=size2)

for metal in ['Co','Ir','Mn','Rh']:
    Exval=dict_MNC[metal+'_H']-dict_MNC[metal+'_slab']-(dict_MNC['Fe_H']-dict_MNC['Fe_slab'])+val_FeNC_H
    Eyval=dict_MNC[metal+'_HCOO']-dict_MNC[metal+'_slab']-(dict_MNC['Fe_HCOO']-dict_MNC['Fe_slab'])+val_FeNC_HCOO
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='o', color='k',s=sdot) 
    if metal=='Mn':
        ax.annotate('MNC', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.2, np.asarray(Eyval).mean()-0.3),color='black',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
    if metal=='Rh' or metal=='Ir' or metal=='Co':
        ax.text(np.asarray(Exval).mean()-0.1, np.asarray(Eyval).mean()+0.02, metal,color='black',alpha=A1, fontsize=size2)
    else:
        ax.text(np.asarray(Exval).mean()+0.02, np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)
    
for metal in ['Fe']:
    Exval=dict_MNC[metal+'_H']-dict_MNC[metal+'_slab']-0.5*np.asarray(EH2)
    Eyval=dict_MNC[metal+'_HCOO']-dict_MNC[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='o', color='k',s=sdot) 
    ax.text(np.asarray(Exval).mean()-0.05, np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)

for metal in ['Pd']:
    Exval=dict_XHg[metal+'_H']-dict_XHg[metal+'_slab']-(dict_XHg['Pt_H']-dict_XHg['Pt_slab'])+val_PtHg4_H
    Eyval=dict_XHg[metal+'_HCOO']-dict_XHg[metal+'_slab']-(dict_XHg['Pt_HCOO']-dict_XHg['Pt_slab'])+val_PtHg4_HCOO
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='s', color='m',s=sdot) 
    ax.annotate('PdHg$_{4}$', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()+0.1, np.asarray(Eyval).mean()+0.1),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))

for metal in ['Pt']:
    Exval=dict_XHg[metal+'_H']-dict_XHg[metal+'_slab']-0.5*np.asarray(EH2)
    Eyval=dict_XHg[metal+'_HCOO']-dict_XHg[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='s', color='g',s=sdot) 
    ax.annotate('PtHg$_{4}$', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()+0.2, np.asarray(Eyval).mean()-0.2),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))

for metal in ['Au_1Pt']:
    Exval=dict_AuPt[metal+'_H']-dict_AuPt[metal+'_slab']-0.5*np.asarray(EH2)
    Eyval=dict_AuPt[metal+'_HCOO_1']-dict_AuPt[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='^', color='c',s=sdot) 
    ax.annotate('Pt$_1$Au(111)', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.35, np.asarray(Eyval).mean()+0.05),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))

plt.xlabel('$\Delta E_{~^*\!\mathrm{H}}$ [eV]',fontsize=size1)
plt.ylabel('$\Delta E_{~^*\!\mathrm{HCOO}}$ [eV]',fontsize=size1)
plt.plot([-0.08429,-0.08429],[-5.25,3.5],'k',lw=2)
plt.plot([-1.0,3.0],[0.12,0.12],'k',lw=2)
ax.annotate('$\Delta E^{0}_{FAOR}$', xy=(0.75,0.12), xytext=(0.72, 0.2),color='black',alpha=1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
ax.annotate('$\Delta E_{H_{UPD}}$', xy=(-0.08429,0.6), xytext=(-0.45, 0.8),color='black',alpha=1, fontsize=size2,arrowprops=dict(arrowstyle="-"))

plt.xlim([-0.5,1.0])
plt.ylim([-1.5,1.0])
fig.text(0.0, 0.9, 'b)', ha='center',weight = 'bold', fontsize=size3+10)
plt.savefig('HCOOvsH.png',dpi=400, bbox_inches='tight')




##############################################################################
# Figure HCOO vs H
##############################################################################
fig = plt.figure(figsize=(7,7));
plt.locator_params(axis='x',nbins=5);plt.grid(True)
plt.xticks(fontsize = size3); plt.yticks(fontsize = size3)
ax=fig.gca()
for metal in ['Ag','Cu','Ni','Pd','Pt','Rh','Cd','In','Ir','Pb','Sn','Tl']:
    Exval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH
    Eyval=dict_Metals[metal+'_HCOO']-dict_Metals[metal+'_slab']-(dict_Metals['Au_HCOO']-dict_Metals['Au_slab'])+val_Au_HCOO
    if metal=='Ag':
        plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P')   
    else:
        plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P') 
    if metal=='Sn':
        ax.text(np.asarray(Exval).mean()-0.15, np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)
    else:
        ax.text(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), metal,color='black',alpha=A1, fontsize=size2)

for metal in ['Au']:
    Exval=dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-ECO2_metal-0.5*np.asarray(EH2_metal)
    Eyval=dict_Metals[metal+'_HCOO']-dict_Metals[metal+'_slab']-ECO2_metal-0.5*np.asarray(EH2_metal)
    ax.annotate('Metals', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.5, np.asarray(Eyval).mean()+0.3),color='blue',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"), fontname='tex')
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), c=GiveColor,s=sdot2, marker='P') 
    ax.text(np.asarray(Exval).mean()-0.00, np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)


for metal in ['Co','Ir','Mn','Rh','Ru']:
    Exval=dict_MNC[metal+'_COOH']-dict_MNC[metal+'_slab']-(dict_MNC['Fe_COOH']-dict_MNC['Fe_slab'])+val_FeNC_COOH
    Eyval=dict_MNC[metal+'_HCOO']-dict_MNC[metal+'_slab']-(dict_MNC['Fe_HCOO']-dict_MNC['Fe_slab'])+val_FeNC_HCOO
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='o', color='k',s=sdot) 
    if metal=='Mo':
        ax.text(np.asarray(Exval).mean()-0.1, np.asarray(Eyval).mean()+0.05, metal,color='black',alpha=A1, fontsize=size2)
    if metal=='Ir':
        ax.text(np.asarray(Exval).mean(), np.asarray(Eyval).mean()-0.1, metal,color='black',alpha=A1, fontsize=size2)
        ax.annotate('MNC', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()+0.1, np.asarray(Eyval).mean()+0.2),color='black',alpha=0.8, fontsize=size1,arrowprops=dict(arrowstyle="-"))
    else:
        ax.text(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), metal,color='black',alpha=A1, fontsize=size2)
    
for metal in ['Fe']:
    Exval=dict_MNC[metal+'_COOH']-dict_MNC[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    Eyval=dict_MNC[metal+'_HCOO']-dict_MNC[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='o', color='k',s=sdot) 
    ax.text(np.asarray(Exval).mean(), np.asarray(Eyval).mean(), metal,color='black',alpha=A1, fontsize=size2)

for metal in ['Pd']:
    Exval=dict_XHg[metal+'_COOH']-dict_XHg[metal+'_slab']-(dict_XHg['Pt_COOH']-dict_XHg['Pt_slab'])+val_PtHg4_COOH
    Eyval=dict_XHg[metal+'_HCOO']-dict_XHg[metal+'_slab']-(dict_XHg['Pt_HCOO']-dict_XHg['Pt_slab'])+val_PtHg4_HCOO
    ax.annotate('PdHg$_{4}$', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()+0.2, np.asarray(Eyval).mean()+0.1),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='s', color='m',s=sdot) 

for metal in ['Pt']:
    Exval=dict_XHg[metal+'_COOH']-dict_XHg[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    Eyval=dict_XHg[metal+'_HCOO']-dict_XHg[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    ax.annotate('PtHg$_{4}$', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.3, np.asarray(Eyval).mean()-0.1),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='s', color='g',s=sdot) 

for metal in ['Au_1Pt']:
    Exval=dict_AuPt[metal+'_COOH']-dict_AuPt[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    Eyval=dict_AuPt[metal+'_HCOO_1']-dict_AuPt[metal+'_slab']-ECO2-0.5*np.asarray(EH2)
    ax.annotate('Pt$_1$Au(111)', xy=(np.asarray(Exval).mean(), np.asarray(Eyval).mean()), xytext=(np.asarray(Exval).mean()-0.5, np.asarray(Eyval).mean()),color='black',alpha=A1, fontsize=size2,arrowprops=dict(arrowstyle="-"))
    plt.scatter(np.asarray(Exval).mean(), np.asarray(Eyval).mean(),marker='^', color='c',s=sdot) 

plt.ylabel('$\Delta E_{~^*\!\mathrm{OOCH}}$ [eV]',fontsize=size1)
plt.xlabel('$\Delta E_{~^*\!\mathrm{COOH}}$ [eV]',fontsize=size1)
plt.plot([-2,3],[-2,3],'k--',lw=2)
plt.xlim([-0.6,1.0])
plt.ylim([-0.6,1.0])
fig.text(0.0, 0.9, 'c)', ha='center',weight = 'bold', fontsize=size3+10)
plt.savefig('HCOOvsCOOH.png',dpi=400, bbox_inches='tight')


##############################################################################
# Free energy diagram at U = -O.12, 0, 0.17  (V vs RHE)
##############################################################################
metal='Pt'
ECOOH_Pt=(dict_Metals[metal+'_COOH']-dict_Metals[metal+'_slab']-(dict_Metals['Au_COOH']-dict_Metals['Au_slab'])+val_Au_COOH).mean()
EH_Pt=(dict_Metals[metal+'_H']-dict_Metals[metal+'_slab']-(dict_Metals['Au_H']-dict_Metals['Au_slab'])+val_Au_H).mean()+0.08429
L1=1
L2=4
#V=-0.12

for V in [-0.12, 0.0,0.17]:
    fig = plt.figure(figsize=(7,7));plt.grid(True)
    plt.xticks(fontsize = size3); plt.yticks(fontsize = size3)
    plt.subplot(2, 1, 1)
    if V==-0.12:
        plt.annotate('FAOR@%s V$_{\mathrm{RHE}}$' % np.round(V,2), xy=(0.03, 0.1), xycoords='axes fraction', fontsize=size1-4,
            bbox=dict(boxstyle='round',facecolor='white', alpha=1.0),
            horizontalalignment='left', verticalalignment='bottom')
    else:
        plt.annotate('FAOR@%s V$_{\mathrm{RHE}}$' % np.round(V,2), xy=(0.03, 0.8), xycoords='axes fraction', fontsize=size1-4,
             bbox=dict(boxstyle='round',facecolor='white', alpha=1.0),
             horizontalalignment='left', verticalalignment='bottom')
    dx=0.2
    xval1=1; yval1=0
    plt.plot([xval1-dx,xval1+dx],[yval1,yval1],'k-',linewidth=L)
    xval2=2; yval21=ECOOH_Pt-0.12-(0.12+V)
    plt.plot([xval2-dx,xval2+dx],[yval21,yval21],'b--',linewidth=L2)
    xval=2; yval22=-0.08429+0.29-0.12-(0.12+V)
    plt.plot([xval1+dx,xval2-dx],[yval1,yval21],'b', linestyle='dotted',linewidth=L1)
    plt.plot([xval1+dx,xval2-dx],[yval1,yval22],'k', linestyle='dotted',linewidth=L1)
    plt.plot([xval2-dx,xval2+dx],[yval22,yval22],'k--',linewidth=L2)
    plt.text(xval2-dx,yval22+0.015,'Metal',fontsize=size2)
    plt.text(xval2-0.5*dx,yval21+0.015,'Pt',color='b',fontsize=size2)
    dx=0.2
    xval3=3; yval3=-0.24-2*V
    plt.plot([xval2+dx,xval3-dx],[yval21,yval3],'b', linestyle='dotted',linewidth=L1)
    plt.plot([xval2+dx,xval3-dx],[yval22,yval3],'k', linestyle='dotted',linewidth=L1)
    plt.plot([xval3-dx,xval3+dx],[yval3,yval3],'k-',linewidth=L)
    plt.xticks([1,2,3], ['HCOOH', '$^*\!$COOH','CO$_2$'],rotation=0, fontsize = size3); plt.yticks(fontsize = size3)
    plt.ylim([-0.6,0.3])
    plt.ylabel('$\Delta G$ [eV]',fontsize=size1)
    
    plt.subplot(2, 1, 2)
    if V==-0.12 or V==0.0:
        plt.annotate('HOR@%s V$_{\mathrm{RHE}}$' % np.round(V,2), xy=(0.03, 0.1), xycoords='axes fraction', fontsize=size1-4,
            bbox=dict(boxstyle='round',facecolor='white', alpha=1.0),
            horizontalalignment='left', verticalalignment='bottom')
    elif V==0.17:
        plt.annotate('HOR@%s V$_{\mathrm{RHE}}$' % np.round(V,2), xy=(0.03, 0.8), xycoords='axes fraction', fontsize=size1-4,
            bbox=dict(boxstyle='round',facecolor='white', alpha=1.0),
            horizontalalignment='left', verticalalignment='bottom')
    xval1=1; yval1=0
    plt.plot([xval1-dx,xval1+dx],[yval1,yval1],'k-',linewidth=L)
    xval2=2; yval21=EH_Pt-V
    plt.plot([xval2-dx,xval2+dx],[yval21,yval21],'b--',linewidth=L2)
    xval2=2; yval22=-V
    plt.plot([xval2-dx,xval2+dx],[yval22,yval22],'k--',linewidth=L2)
    plt.plot([xval1+dx,xval2-dx],[yval1,yval21],'b', linestyle='dotted',linewidth=L1)
    plt.plot([xval1+dx,xval2-dx],[yval1,yval22],'k', linestyle='dotted',linewidth=L1)
    plt.text(xval2-dx,yval22+0.015,'Metal',fontsize=size2)
    plt.text(xval2-0.5*dx,yval21+0.015,'Pt',color='b',fontsize=size2)
    dx=0.2
    xval3=3; yval3=-2*V
    plt.plot([xval2+dx,xval3-dx],[yval21,yval3],'b', linestyle='dotted',linewidth=L1)
    plt.plot([xval2+dx,xval3-dx],[yval22,yval3],'k', linestyle='dotted',linewidth=L1)
    plt.plot([xval3-dx,xval3+dx],[yval3,yval3],'k-',linewidth=L)
    plt.xticks([1,2,3], ['1/2H$_2$', '$^*\!$H','H$^+$+e$^-$'],rotation=0, fontsize = size3); plt.yticks(fontsize = size3)
    plt.ylim([-0.6,0.3])
    plt.ylabel('$\Delta G$ [eV]',fontsize=size1)
    if V==-0.12:
        fig.text(0.0, 0.9, 'd)', ha='center',weight = 'bold', fontsize=size3+10)
    elif V==0.0:
        fig.text(0.0, 0.9, 'e)', ha='center',weight = 'bold', fontsize=size3+10)
    elif V==0.17:
        fig.text(0.0, 0.9, 'f)', ha='center',weight = 'bold', fontsize=size3+10)
    plt.savefig('FreeEnergyDiagram_%s.png' % V,dpi=400,bbox_inches='tight')
