Tutorial-04 - Dr. Nabil Khossossi

Prerequisites & Setup

Before you begin

Required Background

Completion of Tutorial 01 (ASE & DFT Fundamentals)
Completion of Tutorial 02 (Python for Materials Science)
Completion of Tutorial 03 (Crystal Structures & Symmetry)
Access to VASP (license required) or alternative DFT code
Basic understanding of quantum mechanics and solid-state physics

Alternative DFT Codes Supported

While this tutorial focuses on VASP, the concepts and workflows can be adapted to:

Quantum ESPRESSO (Open Source) GPAW (ASE-native) SIESTA / FHI-aims CP2K (AIMD)

Note: ASE provides calculator interfaces for most DFT codes, making the Python workflows largely transferable.

Environment Setup

bash

# Clone the tutorial repository
git clone https://github.com/nabkh/Tutorial-04-VASP-Catalysis.git
cd Tutorial-04-VASP-Catalysis

# Create conda environment
conda env create -f environment.yml
conda activate vasp-catalysis

# Verify installation
python -c "import ase; import pymatgen; print('Ready!')"

VASP Input Files Mastery

1.5 hours | Foundation

Learning Objectives

Understand the role and syntax of all four VASP input files
Configure INCAR parameters for different calculation types
Generate POSCAR files from ASE and pymatgen structures
Select appropriate KPOINTS meshes for surfaces vs bulk
Choose correct pseudopotentials (POTCAR) for your elements

Open Jupyter Notebook

1.1 The Four Essential Files

Every VASP calculation requires exactly four input files. Understanding their interplay is crucial for successful simulations.

File	Purpose	Key Parameters
INCAR	Calculation settings	ENCUT, EDIFF, ISMEAR, IBRION, NSW
POSCAR	Atomic structure	Lattice vectors, atomic positions
KPOINTS	k-point sampling	Mesh density, symmetry
POTCAR	Pseudopotentials	PAW potentials for each element

1.2 INCAR Deep Dive: Catalysis Settings

The INCAR file controls how VASP performs the calculation. For surface catalysis, specific settings are critical.

INCAR - Surface Relaxation

# System identification
SYSTEM = Pt111_surface_relaxation

# Electronic convergence
ENCUT  = 450        # Plane-wave cutoff (eV) - check POTCAR ENMAX
EDIFF  = 1E-6       # Energy convergence (eV)
PREC   = Accurate   # Precision level
ALGO   = Fast       # Electronic minimization algorithm

# Smearing for metals
ISMEAR = 1          # Methfessel-Paxton (metals/surfaces)
SIGMA  = 0.1        # Smearing width (eV)

# Ionic relaxation
IBRION = 2          # Conjugate gradient
NSW    = 200        # Maximum ionic steps
EDIFFG = -0.02      # Force convergence (eV/A) - negative = forces
ISIF   = 2          # Relax ions only, fix cell

# Surface-specific settings
LDIPOL = .TRUE.     # Dipole corrections
IDIPOL = 3          # Dipole in z-direction (perpendicular to surface)
LVHAR  = .TRUE.     # Write electrostatic potential

# Parallelization
NCORE  = 4          # Cores per orbital band
LREAL  = Auto       # Real-space projection (large systems)

Pro Tip: ENCUT Selection

Always check the ENMAX values in your POTCAR files. Use ENCUT = 1.3 * max(ENMAX) for production calculations. For Pt, ENMAX is typically ~230 eV, so ENCUT = 400-450 eV is appropriate.

1.3 Generating POSCAR with ASE

ASE provides an elegant interface for creating and manipulating structures that can be directly exported to VASP format.

Python - ASE to POSCAR

from ase.build import fcc111, add_adsorbate
from ase.io import write
from ase.constraints import FixAtoms

# Create Pt(111) surface: 3x3 supercell, 4 layers
slab = fcc111('Pt', size=(3, 3, 4), vacuum=15.0)

# Fix bottom 2 layers (bulk-like region)
z_positions = slab.positions[:, 2]
z_min = z_positions.min()
z_threshold = z_min + 2.5  # ~2 layer thickness

constraint = FixAtoms(indices=[atom.index for atom in slab
                               if atom.position[2] < z_threshold])
slab.set_constraint(constraint)

# Export to VASP format
write('POSCAR', slab, format='vasp', vasp5=True)

print(f"Created slab with {len(slab)} atoms")
print(f"Fixed atoms: {len(constraint.index)}")
print(f"Free atoms: {len(slab) - len(constraint.index)}")

1.4 KPOINTS Strategy for Surfaces

Surface calculations require different k-point strategies than bulk calculations due to the vacuum region.

KPOINTS - Gamma-centered for surfaces

Automatic mesh
0              # Auto-generate
Gamma          # Gamma-centered (better for hexagonal surfaces)
4 4 1          # Dense in-plane, single point perpendicular
0 0 0          # No shift

Critical: k-points in z-direction

For slab calculations with vacuum, always use 1 k-point in the z-direction (perpendicular to the surface). The vacuum region breaks periodicity, making additional k-points physically meaningless and computationally wasteful.

Alternative Codes: Input File Equivalents

The ASE calculator interface makes it easy to adapt these workflows:

QE: pw.x uses &CONTROL, &SYSTEM namelist GPAW: Direct Python calculator

Surface Slab Models

2 hours | Core Methodology

Learning Objectives

Convert Miller indices to physical surface terminations
Build symmetric and asymmetric slab models
Determine optimal vacuum thickness to eliminate spurious interactions
Test convergence with respect to slab thickness
Apply dipole corrections for polar/asymmetric surfaces

Open Jupyter Notebook: Slab Generation

2.1 Miller Indices and Surface Termination

Miller indices (hkl) define the crystallographic plane exposed at the surface. For FCC metals like Pt:

Surface	Coordination	Surface Energy	Catalytic Sites
(111)	CN = 9	Lowest	Top, Bridge, FCC-hollow, HCP-hollow
(100)	CN = 8	Medium	Top, Bridge, Hollow
(110)	CN = 7	Higher	Top, Short-bridge, Long-bridge
(211)	Stepped	Highest	Step edge, Terrace sites

Python - Comparing Surface Facets

from ase.build import fcc111, fcc100, fcc110
from ase.visualize import view
import matplotlib.pyplot as plt

# Create different Pt surface facets
surfaces = {
    'Pt(111)': fcc111('Pt', size=(4,4,4), vacuum=10),
    'Pt(100)': fcc100('Pt', size=(4,4,4), vacuum=10),
    'Pt(110)': fcc110('Pt', size=(4,2,4), vacuum=10),
}

# Calculate surface atom density
for name, slab in surfaces.items():
    cell = slab.get_cell()
    area = np.linalg.norm(np.cross(cell[0], cell[1]))
    print(f"{name}: Area = {area:.2f} A^2")

2.2 Surface Energy Calculation

γ = (E_slab − N × E_bulk) / (2 × A)

Surface energy (eV/Å²) | N = number of atoms, A = surface area

The slab must be thick enough that the center atoms exhibit bulk-like behavior. This requires systematic testing.

Python - Convergence Testing

from ase.build import fcc111
from ase.calculators.vasp import Vasp
import numpy as np

# Test slab thickness from 3 to 7 layers
layers_range = range(3, 8)
surface_energies = []

for n_layers in layers_range:
    # Create slab
    slab = fcc111('Pt', size=(2,2,n_layers), vacuum=15.0)

    # Setup VASP calculator
    calc = Vasp(
        encut=450,
        kpts=(6,6,1),
        xc='PBE',
        ismear=1,
        sigma=0.1,
        ediff=1e-6,
        ibrion=2,
        nsw=100,
        ediffg=-0.02,
        ldipol=True,
        idipol=3,
        directory=f'slab_{n_layers}L'
    )

    slab.calc = calc
    E_slab = slab.get_potential_energy()

    # Calculate surface energy
    n_atoms = len(slab)
    A = np.linalg.norm(np.cross(slab.cell[0], slab.cell[1]))
    E_bulk_per_atom = -6.05  # Pt bulk energy (eV/atom)

    gamma = (E_slab - n_atoms * E_bulk_per_atom) / (2 * A)
    surface_energies.append(gamma)
    print(f"{n_layers} layers: gamma = {gamma:.4f} eV/A^2")

Convergence Criterion

Surface energy is converged when adding more layers changes gamma by less than 0.001 eV/A² (approximately 1 meV/A²). For Pt(111), typically 4-5 layers are sufficient.

2.3 Dipole Corrections

Asymmetric slabs (e.g., with adsorbates on one side) create artificial dipole interactions across the periodic boundary. VASP can correct for this.

INCAR - Dipole Correction Settings

# Dipole corrections for asymmetric slabs
LDIPOL = .TRUE.      # Enable dipole corrections
IDIPOL = 3           # Direction: 1=x, 2=y, 3=z
DIPOL  = 0.5 0.5 0.5 # Center of cell (fractional coordinates)

# For accurate workfunction calculations
LVHAR  = .TRUE.      # Write LOCPOT (electrostatic potential)

Adsorption & Binding Sites

2 hours | Catalysis Core

Learning Objectives

Enumerate all symmetry-unique adsorption sites on surfaces
Calculate adsorption energies with proper reference states
Understand coverage effects and lateral interactions
Build and analyze CO adsorption on Pt(111) - the classic benchmark

Open Jupyter Notebook: Adsorption Calculations

3.1 Adsorption Sites on Pt(111)

The Pt(111) surface has four distinct high-symmetry adsorption sites:

Site	Description	Coordination	Typical Adsorbates
Top (atop)	Directly above surface atom	1	CO, NO, H
Bridge	Between two surface atoms	2	O, H
FCC hollow	Three-fold site, no atom below	3	O, S, H
HCP hollow	Three-fold site, atom directly below	3	O, S, H

Python - Add CO to Pt(111)

from ase.build import fcc111, add_adsorbate
from ase import Atoms
from ase.visualize import view

# Create clean Pt(111) surface
slab = fcc111('Pt', size=(3,3,4), vacuum=15.0)

# Create CO molecule (C-O bond length ~1.13 A)
# add_adsorbate places the first atom at the specified height
# CO adsorbs via C atom, so C should be closer to surface

# Top site adsorption
add_adsorbate(slab, 'C', height=1.85, position='ontop')
add_adsorbate(slab, 'O', height=1.85+1.13, position='ontop')

# Alternative: Create CO as molecule and add
from ase.build import molecule
co = molecule('CO')

# View the structure
view(slab)

3.2 Adsorption Energy Calculation

The adsorption energy quantifies the binding strength between adsorbate and surface.

E_ads = E_{slab+adsorbate} − E_{clean slab} − E_{adsorbate(gas)}

Adsorption energy (eV) | Negative values = favorable (exothermic) binding

Python - Calculate E_ads for CO/Pt(111)

from ase.calculators.vasp import Vasp

def calculate_adsorption_energy(slab_with_ads, clean_slab, adsorbate_gas):
    """
    Calculate adsorption energy.

    Parameters:
    -----------
    slab_with_ads : ASE Atoms
        Optimized slab with adsorbate
    clean_slab : ASE Atoms
        Optimized clean slab
    adsorbate_gas : ASE Atoms
        Optimized gas-phase adsorbate molecule

    Returns:
    --------
    E_ads : float
        Adsorption energy in eV
    """
    E_total = slab_with_ads.get_potential_energy()
    E_slab = clean_slab.get_potential_energy()
    E_gas = adsorbate_gas.get_potential_energy()

    E_ads = E_total - E_slab - E_gas
    return E_ads

# Example usage
E_ads_CO = calculate_adsorption_energy(pt111_co, pt111_clean, co_gas)
print(f"CO adsorption energy on Pt(111): {E_ads_CO:.3f} eV")
# Expected: ~ -1.5 to -1.8 eV for top site (PBE)

The CO Puzzle

Standard DFT (PBE/PW91) predicts CO prefers the hollow site on Pt(111), contradicting experiments showing top-site preference. This is the famous "CO puzzle" - solved by using higher-level functionals (BEEF-vdW, RPBE) or many-body corrections.

Spin-Polarized DFT

1 hour | Magnetic Systems

Learning Objectives

Recognize when spin-polarized calculations are necessary
Configure ISPIN and MAGMOM parameters correctly
Understand ferromagnetic (FM) vs antiferromagnetic (AFM) ordering
Apply spin-polarization to Fe, Co, Ni-based catalysts

Open Jupyter Notebook

4.1 When is Spin-Polarization Needed?

System Type	ISPIN Setting	Examples
Non-magnetic metals	ISPIN = 1	Pt, Au, Ag, Cu, Pd
Ferromagnetic metals	ISPIN = 2	Fe, Co, Ni
Antiferromagnetic oxides	ISPIN = 2	NiO, FeO, MnO
Paramagnetic molecules	ISPIN = 2	O₂, NO, radicals

INCAR - Spin-Polarized Fe(110)

SYSTEM = Fe110_ferromagnetic

# Enable spin polarization
ISPIN  = 2          # Spin-polarized calculation
MAGMOM = 36*2.5     # Initial magnetic moment: 36 Fe atoms x 2.5 mu_B each

# For AFM ordering, specify alternating moments
# MAGMOM = 18*2.5 18*-2.5  # Half up, half down

# Standard settings
ENCUT  = 450
EDIFF  = 1E-6
ISMEAR = 1
SIGMA  = 0.1

4.2 NiFe Alloy for HER - Research Application

NiFe alloys are promising earth-abundant catalysts for hydrogen evolution. Both elements are magnetic, requiring careful treatment.

Python - NiFe Alloy Surface

from ase.build import fcc111
from ase import Atoms
import numpy as np

# Create Ni(111) slab as template
slab = fcc111('Ni', size=(3,3,4), vacuum=15.0, a=3.52)

# Substitute 25% Ni with Fe (random or ordered)
n_substitute = len(slab) // 4
fe_indices = np.random.choice(len(slab), n_substitute, replace=False)

for idx in fe_indices:
    slab[idx].symbol = 'Fe'

# Set initial magnetic moments
magmoms = []
for atom in slab:
    if atom.symbol == 'Ni':
        magmoms.append(0.6)   # Ni: ~0.6 mu_B
    else:
        magmoms.append(2.5)   # Fe: ~2.5 mu_B

slab.set_initial_magnetic_moments(magmoms)
print(f"NiFe alloy: {sum(s=='Ni' for s in slab.symbols)} Ni, {sum(s=='Fe' for s in slab.symbols)} Fe")

Research Reference

For detailed analysis of NiFe catalysts for hydrogen evolution, see our published work on ASCICat and related systems.

DFT+U for Transition Metals

1.5 hours | Strongly Correlated Systems

Learning Objectives

Understand why standard DFT fails for transition metal oxides
Apply Hubbard U correction to localized d-electrons
Choose appropriate U values from literature or linear response
Model battery cathode materials (LiCoO₂, NMC)

Open Jupyter Notebook

5.1 The Self-Interaction Problem

Standard DFT suffers from self-interaction error, causing excessive delocalization of d-electrons. DFT+U adds an on-site Coulomb interaction to correct this.

E_DFT+U = E_DFT + (U − J)/2 × Σ_σ [Tr(n^σ) − Tr(n^σ × n^σ)]

Dudarev formulation | U = Coulomb parameter, J = exchange (often set to 0)

INCAR - DFT+U for CoO₂

SYSTEM = LiCoO2_battery_cathode

# DFT+U settings (Dudarev approach)
LDAU     = .TRUE.       # Enable DFT+U
LDAUTYPE = 2            # Dudarev: simplified rotationally invariant
LDAUL    = -1 2 -1      # L value: Li(-1=off), Co(2=d), O(-1=off)
LDAUU    = 0.0 3.3 0.0  # U values (eV)
LDAUJ    = 0.0 0.0 0.0  # J values (0 for Dudarev)
LMAXMIX  = 4            # For d-electrons

# Recommended U values for common TM oxides:
# Co: 3.3-3.5 eV, Fe: 4.0-4.3 eV, Ni: 6.0-6.4 eV, Mn: 3.9-4.0 eV
# Ti: 3.0-4.0 eV, V: 3.0-3.5 eV, Cu: 4.0-6.0 eV

Choosing U Values

Common approaches: (1) Use values from Materials Project database, (2) Fit to experimental band gap or formation energy, (3) Calculate from linear response theory. Be consistent within your study.

5.2 Transition Metal Oxide Surfaces

Python - TiO₂ (110) with DFT+U

from ase.build import surface
from ase.spacegroup import crystal
from ase.calculators.vasp import Vasp

# Build rutile TiO2
tio2 = crystal(
    ['Ti', 'O'],
    basis=[(0, 0, 0), (0.305, 0.305, 0)],
    spacegroup=136,  # P4_2/mnm
    cellpar=[4.594, 4.594, 2.958, 90, 90, 90]
)

# Create (110) surface
tio2_110 = surface(tio2, (1,1,0), layers=4, vacuum=15)

# Setup DFT+U calculator
calc = Vasp(
    xc='PBE',
    encut=450,
    kpts=(4,4,1),
    ldau=True,
    ldautype=2,
    ldau_luj={'Ti': {'L': 2, 'U': 3.5, 'J': 0},
              'O':  {'L': -1, 'U': 0, 'J': 0}},
    lmaxmix=4,
    ismear=0,  # Gaussian for semiconductors
    sigma=0.05
)

Transition States & Reaction Barriers

2.5 hours | Advanced Methods

Learning Objectives

Understand the concept of minimum energy path (MEP) and saddle points
Apply Nudged Elastic Band (NEB) method for barrier calculations
Use Climbing Image NEB (CI-NEB) for accurate transition states
Model CO oxidation on Pt(111) as a benchmark reaction

Open Jupyter Notebook: NEB & Dimer Methods

6.1 Nudged Elastic Band (NEB) Method

NEB finds the minimum energy path between two known states (reactants and products) by optimizing a chain of images connected by springs.

E_a = E_TS − E_IS

Activation barrier (eV) | TS = transition state, IS = initial state

Python - NEB Setup with ASE

from ase.neb import NEB
from ase.optimize import BFGS
from ase.io import read, write
import matplotlib.pyplot as plt

# Load optimized initial and final states
initial = read('initial_state.vasp')  # CO + O on surface (separate)
final = read('final_state.vasp')      # CO2 on surface

# Create intermediate images via interpolation
n_images = 7  # 5 intermediate + 2 endpoints
images = [initial.copy() for _ in range(n_images)]
images[-1] = final.copy()

# Setup NEB with climbing image
neb = NEB(images, climb=True, k=0.5)  # k = spring constant
neb.interpolate()  # Linear interpolation

# Attach VASP calculators to intermediate images
for i, image in enumerate(images[1:-1]):
    calc = Vasp(
        encut=400,
        kpts=(4,4,1),
        ibrion=-1,  # No ionic relaxation (ASE controls)
        nsw=0,
        directory=f'neb_image_{i+1:02d}'
    )
    image.calc = calc

# Optimize the NEB path
optimizer = BFGS(neb, trajectory='neb.traj')
optimizer.run(fmax=0.05)

# Extract barrier
energies = [img.get_potential_energy() for img in images]
barrier = max(energies) - energies[0]
print(f"Activation barrier: {barrier:.3f} eV")

6.2 VASP Native NEB (VTST Tools)

For production calculations, VASP's built-in NEB (via VTST add-on) is often more efficient.

INCAR - VASP NEB Settings

SYSTEM = CO_oxidation_NEB

# NEB-specific settings
IMAGES = 5          # Number of intermediate images
SPRING = -5.0       # Spring constant (eV/A^2)
LCLIMB = .TRUE.     # Climbing image for accurate TS

# Ionic relaxation (NEB-optimized)
IBRION = 3          # Quick-min (better for NEB)
POTIM  = 0.0        # Automatic timestep
NSW    = 300        # Maximum steps
EDIFFG = -0.03      # Force convergence

# Standard settings
ENCUT  = 400
EDIFF  = 1E-5
ISMEAR = 1
SIGMA  = 0.1

6.3 CO Oxidation on Pt(111): Complete Example

CO oxidation (CO + O* -> CO₂) is a benchmark reaction for heterogeneous catalysis.

Expected Results for CO + O -> CO₂ on Pt(111)

Activation barrier: ~0.6-0.9 eV (PBE)
Reaction energy: ~-1.5 eV (exothermic)
TS geometry: CO tilted toward O, C-O bond elongated

Computational Hydrogen Electrode

1.5 hours | Electrochemistry

Learning Objectives

Apply the Computational Hydrogen Electrode (CHE) model
Calculate free energy diagrams for HER and OER
Include zero-point energy and entropy corrections
Construct volcano plots for catalyst screening

Open Jupyter Notebook: CHE Model

7.1 The Norskov CHE Framework

The CHE model connects DFT energies to electrochemical potentials through a simple but powerful approximation.

G(H⁺ + e⁻) = ½ G(H₂)

At U = 0 V vs. RHE and pH = 0

ΔG_H* = ΔE_H + ΔZPE − TΔS

Free energy of H* intermediate | ZPE ≈ 0.04 eV, TΔS ≈ −0.20 eV at 298 K

Practical Approximation

For hydrogen adsorption: ΔG_H* ≈ ΔE_H + 0.24 eV

This allows calculating electrochemical energetics from gas-phase DFT!

Python - HER Free Energy Diagram

import numpy as np
import matplotlib.pyplot as plt

def calculate_her_free_energy(dE_H, U=0.0, pH=0):
    """
    Calculate HER free energy at given potential and pH.

    HER: H+ + e- -> H* (adsorption step)
         H* + H+ + e- -> H2 (desorption/Heyrovsky step)

    Parameters:
    -----------
    dE_H : float
        Hydrogen adsorption energy from DFT (eV)
    U : float
        Electrode potential vs RHE (V)
    pH : float
        Solution pH

    Returns:
    --------
    dG_H : float
        Free energy of H* intermediate
    """
    # Zero-point energy correction (from vibrations)
    ZPE_H_ads = 0.17   # H on metal surface (eV)
    ZPE_H2_gas = 0.27  # H2 molecule (eV)
    dZPE = ZPE_H_ads - 0.5 * ZPE_H2_gas

    # Entropy correction at 298 K
    T = 298.15  # K
    S_H2_gas = 0.41  # eV at 298 K, 1 bar
    TdS = -0.5 * T * S_H2_gas / 298.15  # ~ -0.20 eV

    # Free energy at U=0, pH=0
    dG_H = dE_H + dZPE - TdS  # ~ dE_H + 0.24 eV

    # Potential and pH corrections
    dG_H -= U  # Shift with applied potential
    dG_H += 0.059 * pH  # pH correction (kT*ln(10)*pH)

    return dG_H

# Example: Calculate for Pt(111)
dE_H_Pt = -0.35  # DFT adsorption energy (eV)
dG_H_Pt = calculate_her_free_energy(dE_H_Pt, U=0, pH=0)
print(f"Pt(111) dG_H* = {dG_H_Pt:.3f} eV")  # ~ -0.1 eV (near optimal!)

7.2 Volcano Plots

The Sabatier principle: optimal catalysts have intermediate binding strength. Too weak = poor activation; too strong = surface poisoning.

Python - HER Volcano Plot

import matplotlib.pyplot as plt
import numpy as np

# Sabatier volcano for HER
dG_range = np.linspace(-0.8, 0.8, 100)

# Activity limited by stronger of two steps
# Step 1: H+ + e- -> H* (limited when dG_H > 0)
# Step 2: H* -> 1/2 H2 (limited when dG_H < 0)
activity = -np.abs(dG_range)  # Log(exchange current) ~ -|dG_H|

# Experimental data points
catalysts = {
    'Pt': (-0.09, -0.1),
    'Pd': (-0.15, -0.2),
    'Ni': (-0.25, -0.5),
    'Au': (0.10, -0.6),
    'MoS2': (0.08, -0.3),
    'NiFe': (-0.12, -0.25),  # Our research!
}

plt.figure(figsize=(8, 6))
plt.plot(dG_range, activity, 'k-', lw=2, label='Volcano')

for name, (dG, log_i0) in catalysts.items():
    plt.scatter(dG, log_i0, s=100, label=name, zorder=5)

plt.xlabel(r'$\Delta G_{H*}$ (eV)', fontsize=12)
plt.ylabel(r'log($i_0$) / Activity', fontsize=12)
plt.legend()
plt.title('HER Volcano Plot')
plt.tight_layout()
plt.savefig('her_volcano.png', dpi=150)

Best Practices & Workflows

1 hour | Professional Skills

Learning Objectives

Implement systematic convergence testing protocols
Write efficient HPC submission scripts
Ensure reproducibility through documentation and version control
Avoid common pitfalls in computational catalysis

Open Convergence Checklist

8.1 Convergence Testing Protocol

Parameter	Test Range	Convergence Target
ENCUT	300-600 eV (50 eV steps)	Energy change < 1 meV/atom
KPOINTS	2x2x1 to 8x8x1	Energy change < 5 meV
Slab layers	3-7 layers	Surface energy < 1 meV/A²
Vacuum	10-20 A	Work function < 10 meV
Supercell size	2x2 to 4x4	Adsorption energy < 20 meV

8.2 HPC Submission Script

bash - SLURM Job Script

#!/bin/bash
#SBATCH --job-name=vasp_surf
#SBATCH --nodes=2
#SBATCH --ntasks-per-node=48
#SBATCH --time=24:00:00
#SBATCH --partition=compute
#SBATCH --account=your_project

module load vasp/6.3.0
module load intel/2022

# Run VASP
srun vasp_std > vasp.out

# Post-processing
python analyze_output.py

8.3 Common Pitfalls to Avoid

Critical Errors in Catalysis DFT

Wrong reference state: Always use consistent gas-phase references (e.g., O₂/2 for O, H₂/2 for H)
Missing ZPE corrections: Essential for accurate free energies, especially for H-containing species
Unconverged k-points: Surface calculations need dense in-plane meshes
Frozen bottom layers: Always fix bulk-like layers to prevent artificial slab bending
Insufficient vacuum: Minimum 12-15 A to avoid image interactions

Dr. Nabil Khossossi

Prerequisites & Setup

Required Background

Alternative DFT Codes Supported

Environment Setup

VASP Input Files Mastery

Learning Objectives

1.1 The Four Essential Files

1.2 INCAR Deep Dive: Catalysis Settings

Pro Tip: ENCUT Selection

1.3 Generating POSCAR with ASE

1.4 KPOINTS Strategy for Surfaces

Critical: k-points in z-direction

Alternative Codes: Input File Equivalents

Surface Slab Models

Learning Objectives

2.1 Miller Indices and Surface Termination

2.2 Surface Energy Calculation

Convergence Criterion

2.3 Dipole Corrections

Adsorption & Binding Sites

Learning Objectives

3.1 Adsorption Sites on Pt(111)

3.2 Adsorption Energy Calculation

The CO Puzzle

Spin-Polarized DFT

Learning Objectives

4.1 When is Spin-Polarization Needed?

4.2 NiFe Alloy for HER - Research Application

Research Reference

DFT+U for Transition Metals

Learning Objectives

5.1 The Self-Interaction Problem

Choosing U Values

5.2 Transition Metal Oxide Surfaces

Transition States & Reaction Barriers

Learning Objectives

6.1 Nudged Elastic Band (NEB) Method

6.2 VASP Native NEB (VTST Tools)

6.3 CO Oxidation on Pt(111): Complete Example

Expected Results for CO + O -> CO2 on Pt(111)

Computational Hydrogen Electrode

Learning Objectives

7.1 The Norskov CHE Framework

Practical Approximation

7.2 Volcano Plots

Best Practices & Workflows

Learning Objectives

8.1 Convergence Testing Protocol

8.2 HPC Submission Script

8.3 Common Pitfalls to Avoid

Critical Errors in Catalysis DFT

Expected Results for CO + O -> CO₂ on Pt(111)