""" Created on February 8 2018 This script plots the formation energy. @author: Alexander Bagger """ import os import numpy as np import matplotlib.pylab as plt from tempfile import TemporaryFile import sys from matplotlib.patches import Ellipse from scipy.stats import norm, chi2 import matplotlib.patches as patches import ase.db from ase.db import connect EH2=-6.72136774092 # 0.17 grid EH2O=-14.2249062056 # 0.17 grid ECO=-14.8523695434 # 0.17 grid ECO2=-23.1040825464 # Energies needed for calculating G, data obtained from ChemCatChem Noerskov 2014 H2zpl=0.28; H2ts=-0.40; H2cp=0.09; G1_H2=H2zpl+H2ts+H2cp; COzpl=0.13; COts=-0.61; COcp=0.09; G1_CO=COzpl+COts+COcp; CO2zpl=0.31; CO2ts=-0.66; CO2cp=0.10; G1_CO2=CO2zpl+CO2ts+CO2cp; A_COOHzpl=0.62; A_COOHts=-0.178; A_COOHcp=0.10; AG1_COOH=A_COOHzpl+A_COOHts+A_COOHcp; A_COzpl=0.19; A_COts=-0.16; A_COcp=0.08; AG1_CO=A_COzpl+A_COts+A_COcp; # Water and RPBE functional corrections: CO2coor=0.45 WaterCoorOH=0.5 WaterCoorROH=0.25 WaterCoorCO=0.1 # Defining COOH* and CO* adsorbate energies from the following functions. def fCOOH(ECOOH, EBULK, V): E=ECOOH-EBULK-ECO2-G1_CO2-0.5*EH2-0.5*G1_H2+AG1_COOH-CO2coor-WaterCoorROH-WaterCoorCO+V return E def fCO(EnergyCO, EBULK, V): E=EnergyCO-EBULK-ECO-G1_CO+AG1_CO-WaterCoorCO+2*V #E=EnergyCO-EBULK-ECO-G_CO+AG_CO-WaterCoorCO+2*V return E # Load in data from database file dict={} d=ase.db.connect('CO2RR_project.db') # Key words for database: # facet 'Ag111', 'Ag211', 'Ag211K', 'Ag211A', 'Ag111A', 'Cu111', 'Cu211', 'Au111', 'Au211', 'Ag111_O', 'Ag2O' # ads: 'slab', 'CO','COOH' # field: '0', '2', '4', '6' V=0.0 for facet in ['Ag111', 'Ag211', 'Ag211K', 'Ag211A', 'Ag111A', 'Cu111', 'Cu211', 'Au111', 'Au211', 'Ag111_O', 'Ag2O']: if facet=='Ag111' or facet=='Ag211' or facet=='Ag211K' or facet=='Ag211A' or facet=='Ag111A': for k in [0,-2,-4,-6]: dict[facet+'_CO_'+str(k)]=fCO(d.get(facet=facet, ads='CO', field=k).energy,d.get(facet=facet, ads='slab', field=k).energy,V) dict[facet+'_COOH_'+str(k)]=fCOOH(d.get(facet=facet, ads='COOH', field=k).energy,d.get(facet=facet, ads='slab', field=k).energy,V) else: dict[facet+'_CO_0']=fCO(d.get(facet=facet, ads='CO', field=0).energy,d.get(facet=facet, ads='slab', field=0).energy,V) dict[facet+'_COOH_0']=fCOOH(d.get(facet=facet, ads='COOH', field=0).energy,d.get(facet=facet, ads='slab', field=0).energy,V) E_names=['111','211','Kink','Kink add','Add','Cu111','Cu211','Au111','Au211'] listx=[dict['Ag111_CO_0'],dict['Ag211_CO_0'],dict['Ag211K_CO_0'],dict['Ag211A_CO_0'],dict['Ag111A_CO_0'],dict['Cu111_CO_0'],dict['Cu211_CO_0'],dict['Au111_CO_0'],dict['Au211_CO_0']] listy=[dict['Ag111_COOH_0'],dict['Ag211_COOH_0'],dict['Ag211K_COOH_0'],dict['Ag211A_COOH_0'],dict['Ag111A_COOH_0'],dict['Cu111_COOH_0'],dict['Cu211_COOH_0'],dict['Au111_COOH_0'],dict['Au211_COOH_0']] fit = np.polyfit(listx,listy,1) fit_fn = np.poly1d(fit) #fit2 = np.polyfit(G_CO[-3:],G_COOH[-3:],1) #fit_fn2 = np.poly1d(fit2) L=4 L1=2 s3=50 plt.figure() # at applied potential plt.plot([-0.22, 1.0],[-fit_fn(-0.22)+0.6,-fit_fn(1.0)+0.6],'r-',linewidth=L) plt.plot([-0.57-0.02, -0.2-0.02],[-0.57+0.12-0.02,-0.2+0.12-0.02],'r-',linewidth=L,label='With -0.6 V vs RHE') plt.plot([-1.5, -0.575-0.02],[-1.5+0.12,-0.575+0.12-0.02],'b-',linewidth=L,label='Without potential') plt.plot([-0.575, 1.0],[-fit_fn(-0.575),-fit_fn(1.0)],'b-',linewidth=L) plt.plot([-1.5, 1.0],[0.12,0.12],'r--',linewidth=L1, label='Thermodynamics') plt.plot([0.0, 0.0],[-2,1.2],'k--',linewidth=L1)#, label='CO(g) -> CO*') # points. ax=plt.gca() plt.scatter(listx[-4:],-1*np.asarray(listy[-4:]), marker="o",s=s3 , c='k') plt.scatter(listx[:-4],-1*np.asarray(listy[:-4]), marker="s",s=s3 , c='k') plt.scatter(listx[:-4],-1*np.asarray(listy[:-4])+0.6, marker="s",s=s3 , c='k') # Adding names ax.text(listx[5]-0.32,-1*listy[5]-0.06, r'%s' % E_names[5],fontsize=18) ax.text(listx[6]-0.32,-1*listy[6]-0.06, r'%s' % E_names[6],fontsize=18) ax.text(listx[7]+0.05,-1*listy[7], r'%s' % E_names[7],fontsize=18) ax.text(listx[8]-0.32,-1*listy[8]-0.22, r'%s' % E_names[8],fontsize=18) ax.arrow(listx[8], -1*listy[8]-0.19, 0, 0.1, width=0.001, head_width=0.03, head_length=0.05, fc='k', ec='k') # Oxide points plt.scatter(dict['Ag2O_CO_0'],-1*dict['Ag2O_COOH_0']+0.6, marker="s",s=s3 , c='r') plt.scatter(dict['Ag111_O_CO_0'],-1*dict['Ag111_O_COOH_0']+0.6, marker="s",s=s3 , c='r') ellip=Ellipse(xy=[-0.15, -0.38+0.63], width=0.1, height=0.3, angle=55, color='r', alpha=0.1) ax.add_patch(ellip) ax.text(-0.7,-0.28+0.63, r'Ag Subsurface',fontsize=18) ax.text(-0.7,-0.37+0.63, r'Oxygen',fontsize=18) # Adding electric field dots listx2=[dict['Ag111_CO_-2'],dict['Ag211_CO_-2'],dict['Ag211K_CO_-2'],dict['Ag211A_CO_-2'],dict['Ag111A_CO_-2']] listx4=[dict['Ag111_CO_-4'],dict['Ag211_CO_-4'],dict['Ag211K_CO_-4'],dict['Ag211A_CO_-4'],dict['Ag111A_CO_-4']] listx6=[dict['Ag111_CO_-6'],dict['Ag211_CO_-6'],dict['Ag211K_CO_-6'],dict['Ag211A_CO_-6'],dict['Ag111A_CO_-6']] listy2=[dict['Ag111_COOH_-2'],dict['Ag211_COOH_-2'],dict['Ag211K_COOH_-2'],dict['Ag211A_COOH_-2'],dict['Ag111A_COOH_-2']] listy4=[dict['Ag111_COOH_-4'],dict['Ag211_COOH_-4'],dict['Ag211K_COOH_-4'],dict['Ag211A_COOH_-4'],dict['Ag111A_COOH_-4']] listy6=[dict['Ag111_COOH_-6'],dict['Ag211_COOH_-6'],dict['Ag211K_COOH_-6'],dict['Ag211A_COOH_-6'],dict['Ag111A_COOH_-6']] plt.scatter(listx2,-1*np.asarray(listy2)+0.6, marker="p",s=s3, c='y') plt.scatter(listx4,-1*np.asarray(listy4)+0.6, marker="p",s=s3, c='y') plt.scatter(listx6,-1*np.asarray(listy6)+0.6, marker="p",s=s3, c='y') fitE2 = np.polyfit(listx2,listy2,1) fit_fnE2 = np.poly1d(fitE2) plt.plot([0,0.45],[-fit_fnE2(0)+0.6,-fit_fnE2(0.45)+0.6],'y--',linewidth=1) fitE4 = np.polyfit(listx4,listy4,1) fit_fnE4 = np.poly1d(fitE4) plt.plot([0,0.45],[-fit_fnE4(0)+0.6,-fit_fnE4(0.45)+0.6],'y--',linewidth=1) fitE6 = np.polyfit(listx6,listy6,1) fit_fnE6 = np.poly1d(fitE6) plt.plot([0,0.45],[-fit_fnE6(0)+0.6,-fit_fnE6(0.45)+0.6],'y--',linewidth=1) ax.arrow(0.40, -0.5+0.2, 0, 0.22, width=0.0015, head_width=0.04, head_length=0.06, fc='y', ec='y') ax.text(0.45,-0.3+0.2, r'-0.6 V/$\AA$',fontsize=16) ax.text(0.45,-0.4+0.2, r'-0.4 V/$\AA$',fontsize=16) ax.text(0.45,-0.5+0.2, r'-0.2 V/$\AA$',fontsize=16) plt.legend(loc=1, fontsize=11.5) plt.xlim([-1.0,0.8]) plt.ylim([-1.2,0.6]) plt.xlabel('$\Delta$G$_{CO}$ (eV)',fontsize=20) plt.ylabel('-$\Delta$ G (eV)',fontsize=20) plt.xticks(fontsize=16) plt.yticks(fontsize=16) xmin, xmax = -1, 0 ymin, ymax = -1.2, 0.12 X = [[0.0, 0.0], [1.6, 1.6]] ax.imshow(X, interpolation='bicubic',cmap='Greys', extent=(xmin, xmax, ymin, ymax), alpha=0.3) ax.text(-0.95,-1.15, r'Beyond CO (g)',fontsize=20) plt.savefig('CO2RR_vulcano_plot.png',bbox_inches='tight',dpi=300)