import numpy as np

E_O = -2.06642499326
E_CO = -14.60836

Nvalues = [13, 55, 147, 309, 561, 923]

E_Pt = [ -46.90724,
         -242.91267,
         -695.69851,
         -1511.54142648, # -1511.54751
         -2801.00825,
         -4670.75541
         ]

E_O_Pt = [ -53.73204,
           -248.67198,
           -701.23476,
           -1517.09939,
           -2806.4281758,
           -4676.22316
           ]

E_CO_Pt = [-63.83768, -259.20228, -711.45056, -1527.32629, -2816.69506]

de = np.array(E_O_Pt) - np.array(E_Pt) - E_O

descriptions = ['O fcc', 'CO 111 top']
e_mol = [E_O, E_CO]


import pylab as pl
for i, energies in enumerate([E_O_Pt, E_CO_Pt]):
    E = [energy - E_Pt[j] - e_mol[i] for j, energy in enumerate(energies)]
    pl.plot(E, lw=2, label=descriptions[i])
pl.xticks(range(len(Nvalues)), map(str, Nvalues))
pl.xlabel('Number of atoms in cluster')
pl.ylabel('Adsorption energy [eV]')
pl.show()
#for i, N in enumerate(Nvalues):
#    dE = E_O_Pt[i] - E_Pt[i] - E_O
#    print N, dE