** Reference to paper:

https://pubs.acs.org/doi/10.1021/acs.jpclett.0c02452

** DOI:

10.1021/acs.jpclett.0c02452

** Title:

Density Functional Theory for Molecule–Metal Surface Reactions: When Does the Generalized Gradient Approximation Get It Right, and What to Do If It Does Not

** Authors:

Gerrits, Nick; Smeets, Egidius W.F.; Vuckovic, Stefan; Powell, Andrew D.; Doblhoff-Dier, Katharina; Kroes, Geert-Jan

** Contact e-mail:

n.gerrits@lic.leidenuniv.nl
g.j.kroes@chem.leidenuniv.nl

** Abstract:

While density functional theory (DFT) is perhaps the most used electronic structure theory in chemistry, many of its practical aspects remain poorly understood.
For instance, DFT at the generalized gradient approximation (GGA) tends to fail miserably at describing gas-phase reaction barriers, while it performs surprisingly well for many molecule–metal surface reactions.
GGA-DFT also fails for many systems in the latter category, and up to now it has not been clear when one may expect it to work.
We show that GGA-DFT tends to work if the difference between the work function of the metal and the molecule’s electron affinity is greater than ∼7 eV and to fail if this difference is smaller, with sticking of O2 on Al(111) being a spectacular example.
Using dynamics calculations we show that, for this system, the DFT problem may be solved as done for gas-phase reactions, i.e., by resorting to hybrid functionals, but using screening at long-range to obtain a correct description of the metal.
Our results suggest the GGA error in the O2 + Al(111) barrier height to be functional driven.
Our results also suggest the possibility to compute potential energy surfaces for the difficult-to-treat systems with computationally cheap nonself-consistent calculations in which a hybrid functional is applied to a GGA density.

** Description per file:

The program used for DFT is VASP v5.4.4 with a special modification (made by Egidius Smeets) in order to use the MS-RPBEl functional and with VTST (http://theory.cm.utexas.edu/vtsttools/index.html).
Molecular dynamics (MD) are performed with an inhouse developed MD code (QCTraj).
All figures are made with matplotlib 3.0, except the table of content figure which was made with Blender 1.79 and annoted with matplotlib and GIMP.

** Folder DFT_inputs:
KPOINTS - k-point grid employed
INCAR - The input files for the VASP calculations. *_1 is for a non-spin-polarized calculation, whereas *_2 is for spin-polarized calculations restarting from the *_1 calculation. HSEatRPBE is used to perform an NSCF calculation.
The POTCAR is not included due to licensing, where the standard PAW VASP 5.4 data set is used.

** Folder CRP:
pes - Program used to generator POSCARs for the construction of the CRP PES. I did not write the program nor have the source code. The original was written by Egidius Smeets and is likely to be present somewhere in his pubs folder (https://pubs.tc.lic.leidenuniv.nl/public/guido.smeets/). This program essentially makes POSCARs with spherical coordinates and center of mass direct coordinates as an input for a homogeneous diatomic molecule. It implicitly assumes an (111) surface with an angle of 120 degrees between the u and v lattice vectors (NB: typically people use 60 degrees, as such, careful with rotation of the surface sites).

Parameters employed for the cuts:
# [ "u vector, v vector, theta, phi" ]
location    = [ "0 0 0 0", "0 0 90 0", "0 0 90 30", "0.5 0 0 0", "0.5 0 90 0", "0.5 0 90 60", "0.5 0 90 90", "0.66666666 0.33333333 0 0", "0.66666666 0.33333333 45 30", "0.66666666 0.33333333 45 210", "0.66666666 0.33333333 90 0", "0.66666666 0.33333333 90 30", "0.33333333 0.16666666 0 0", "0.33333333 0.16666666 45 30", "0.33333333 0.16666666 45 120", "0.33333333 0.16666666 45 210", "0.33333333 0.16666666 90 30", "0.33333333 0.16666666 90 120", "0.33333333 -0.33333333 0 0", "0.33333333 -0.33333333 45 150", "0.33333333 -0.33333333 45 330", "0.33333333 -0.33333333 90 0", "0.33333333 -0.33333333 90 330", "0.16666666 -0.16666666 0 0", "0.16666666 -0.16666666 45 150", "0.16666666 -0.16666666 45 240", "0.16666666 -0.16666666 45 330", "0.16666666 -0.16666666 90 240", "0.16666666 -0.16666666 90 330", "0 0 90 0", "0 0 90 0" ];
R           = [ 1.0, 1.1, 1.15, 1.175, 1.2, 1.225, 1.25, 1.3, 1.4, 1.5, 1.6 ];#11
Z           = [ 1.00, 1.50, 2.00, 2.25, 2.50, 2.75, 3.00, 3.25, 3.50 ];#9

Zgasphase   = [ 5.00 ]
Rgasphase   = [ 1.100, 1.125, 1.150, 1.175, 1.200, 1.210, 1.225, 1.250, 1.275, 1.300, 1.350, 1.400, 1.450, 1.500, 1.600 ]

Zextrapol   = [ 5.00, 4.75, 4.50, 4.25, 4.00, 3.75, 3.50, 3.25, 3.00, 2.75, 2.50, 2.25, 2.00 ]
Rextrapol   = [ 1.21 ]

The "pes" program is run as follows (read it as a python syntax):
"./pes 4.022386413 2 4 10 fix 0.00000 2.35576845951271 4.6033263509659 6.95593267263204 Molecule {:s} {:1.5f} {:1.5f} > POSCAR".format(location, Z, r)
The first argument is the lattice constant, second is the super cell size (assuming Nx = Ny), third is the number of layers, fourth is the vacuum space, fifth is (as far as I'm aware) always 'fix', after that the positions of the layers (i.e., it takes care of interlayer distances this way), followed by 'Molecule' or 'Atom' for the 6D and 3D PESs, respectively.
The final arguments are the position of the molecule in u, v, theta, phi, Z, and r.

CRProutine.tar.gz - Contains the MS-RPBEl and HSE03-1/3X PESs and an example of how to retrieve energies and forces with an ASE calculator. This tarball is also available with the SI at the publisher's site.

** Folder Figure01:
plotworkfunction.py - Generates the figure.
workfunction.pdf

** Folder Figure02:
plotsticking.py - Generates the figure.
reactionprobability.pdf

** Folder Figure03:
plotreactioncoordinate.py - Generates the figure.
reactioncoordinate.pdf

** Folder FigureS01:
draw_elbow.py - Generates the figure.
elbow.pdf
energy.dat - Contains the data of the figure.
phi_scan.py - Generates the MEP shown in the figure.

** Folder FigureS02:
Same as FigureS01 but for HSE03-1/3X instead of MS-RPBEl.

** Folder FigureS03:
magneticmoment.pdf
plotmagneticmoment.py - Generates the figure.

** Folder FigureS04:
plotsticking_rotational.py - Generates the figure
reactionprobability_rotational.pdf

** Folder FigureS05:
plotsticking_angle.py - Generates the figure
reactionprobability_angle.pdf