09 - Adsorption

This tutorial demonstrates calculating adsorption energies for atoms on metal surfaces, using H on Al(111) as an example.

Note: These calculations use reduced settings for demonstration purposes. Production calculations would use larger slabs, higher ENCUT, and relaxed surfaces.

import numpy as np
from ase.build import add_adsorbate, bulk, fcc111, molecule
from ase.constraints import FixAtoms

from vasp import Vasp

print("=" * 60)
print("Adsorption Energy Calculations")
print("=" * 60)
print()
print("Note: Using fast demo settings (small slab, frozen, low ENCUT)")
print()

============================================================
Adsorption Energy Calculations
============================================================

Note: Using fast demo settings (small slab, frozen, low ENCUT)

Part 1: Reference calculations

print("Part 1: Reference calculations")
print("-" * 40)
print()

# Gas-phase H atom
print("1a. Gas-phase H atom")
from ase import Atoms
h_atom = Atoms('H', positions=[[0, 0, 0]])
h_atom.center(vacuum=6.0)

calc_h = Vasp(
    label='results/adsorption/h_gas',
    atoms=h_atom,
    xc='PBE',
    encut=300,
    kpts=(1, 1, 1),
    ismear=0,
    sigma=0.05,
    ispin=2,
)

e_h = calc_h.potential_energy
print(f"  H energy: {e_h:.6f} eV")
print()

# Al bulk reference
print("1b. Al bulk reference")
al_bulk = bulk('Al', 'fcc', a=4.05)

calc_al_bulk = Vasp(
    label='results/adsorption/al_bulk',
    atoms=al_bulk,
    xc='PBE',
    encut=300,
    kpts=(10, 10, 10),
    ismear=1,
    sigma=0.1,
)

e_al_bulk = calc_al_bulk.potential_energy
e_al_per_atom = e_al_bulk / len(al_bulk)
print(f"  Al bulk energy per atom: {e_al_per_atom:.6f} eV")
print()

Part 1: Reference calculations
----------------------------------------

1a. Gas-phase H atom
  H energy: -1.111562 eV

1b. Al bulk reference
  Al bulk energy per atom: -3.744454 eV

Part 2: Clean Al(111) surface

We use a small 2x2 slab with 3 layers, all frozen (no relaxation). This is much faster but less accurate than a relaxed surface.

print("Part 2: Clean Al(111) surface")
print("-" * 40)
print()

# Create small Al(111) slab - 2x2, 3 layers, 10 A vacuum
slab = fcc111('Al', size=(2, 2, 3), vacuum=10.0, a=4.05)

# Freeze ALL atoms - no surface relaxation for speed
constraint = FixAtoms(indices=range(len(slab)))
slab.set_constraint(constraint)

print("  Slab: Al(111) 2x2, 3 layers (all frozen)")
print(f"  Atoms: {len(slab)}")
print()

# Static calculation only (no relaxation needed - all frozen)
calc_slab = Vasp(
    label='results/adsorption/al111_clean',
    atoms=slab,
    xc='PBE',
    encut=300,
    kpts=(4, 4, 1),
    ismear=1,
    sigma=0.1,
    ldipol=True,
    idipol=3,
)

e_slab = calc_slab.potential_energy
print(f"  Clean slab energy: {e_slab:.6f} eV")
print()

Part 2: Clean Al(111) surface
----------------------------------------

  Slab: Al(111) 2x2, 3 layers (all frozen)
  Atoms: 12

  Clean slab energy: -37.903796 eV

Part 3: H adsorption at different sites

We test two sites (ontop and fcc hollow) to demonstrate the workflow. Only the H atom position is relaxed; the slab remains frozen.

print("Part 3: H adsorption at different sites")
print("-" * 40)
print()

# Test just two sites for speed
adsorption_sites = ['ontop', 'fcc']
adsorption_results = {}

for site in adsorption_sites:
    print(f"  Calculating {site} site...")

    # Create fresh slab with H
    slab_h = fcc111('Al', size=(2, 2, 3), vacuum=10.0, a=4.05)

    # Add H atom at specified site
    add_adsorbate(slab_h, 'H', height=1.5, position=site)

    # Freeze slab atoms only (indices 0 to n_slab-1)
    # Let H (last atom) relax
    n_slab = 12  # 2x2x3 = 12 Al atoms
    constraint = FixAtoms(indices=range(n_slab))
    slab_h.set_constraint(constraint)

    calc_ads = Vasp(
        label=f'results/adsorption/al111_h_{site}',
        atoms=slab_h,
        xc='PBE',
        encut=300,
        kpts=(4, 4, 1),
        ismear=1,
        sigma=0.1,
        ldipol=True,
        idipol=3,
        ibrion=2,
        nsw=20,
        ediffg=-0.05,
    )

    e_ads = calc_ads.potential_energy

    # Calculate adsorption energy
    e_adsorption = e_ads - e_slab - e_h

    adsorption_results[site] = {
        'energy': e_ads,
        'e_ads': e_adsorption,
    }

    print(f"    Total energy: {e_ads:.6f} eV")
    print(f"    Adsorption energy: {e_adsorption:.3f} eV")
    print()

Part 3: H adsorption at different sites
----------------------------------------

  Calculating ontop site...
    Total energy: -44.989159 eV
    Adsorption energy: -5.974 eV

  Calculating fcc site...
    Total energy: -45.144757 eV
    Adsorption energy: -6.129 eV

Part 4: Analysis

print("Part 4: Summary of adsorption energies")
print("-" * 40)
print()

print(f"{'Site':<10} {'E_ads (eV)':<12} {'Stability':<15}")
print("-" * 37)

# Sort by adsorption energy (more negative = more stable)
sorted_sites = sorted(adsorption_results.items(), key=lambda x: x[1]['e_ads'])

for i, (site, data) in enumerate(sorted_sites):
    stability = "Most stable" if i == 0 else f"dE = {data['e_ads'] - sorted_sites[0][1]['e_ads']:.3f} eV"
    print(f"{site:<10} {data['e_ads']:<12.3f} {stability:<15}")

print()

Part 4: Summary of adsorption energies
----------------------------------------

Site       E_ads (eV)   Stability      
-------------------------------------
fcc        -6.129       Most stable    
ontop      -5.974       dE = 0.156 eV  

Summary

print("=" * 60)
print("Summary: H/Al(111) Adsorption")
print("=" * 60)
print()

most_stable = sorted_sites[0]
print(f"Most stable site: {most_stable[0]}")
print(f"Adsorption energy: {most_stable[1]['e_ads']:.3f} eV")
print()
print("Demo settings used (for speed):")
print("  - Small slab: 2x2, 3 layers")
print("  - Frozen substrate (no relaxation)")
print("  - Lower ENCUT (300 eV)")
print("  - Two sites only")
print()
print("For production calculations:")
print("  - Use 3x3 or 4x4 slab, 4+ layers")
print("  - Relax top 2 layers")
print("  - ENCUT 400+ eV")
print("  - Test all sites (ontop, bridge, fcc, hcp)")
print()
print("Key formula:")
print("  E_ads = E(slab+adsorbate) - E(slab) - E(adsorbate)")
print()
print("Next: Try 10_reactions/ for reaction energetics.")

============================================================
Summary: H/Al(111) Adsorption
============================================================

Most stable site: fcc
Adsorption energy: -6.129 eV

Demo settings used (for speed):
  - Small slab: 2x2, 3 layers
  - Frozen substrate (no relaxation)
  - Lower ENCUT (300 eV)
  - Two sites only

For production calculations:
  - Use 3x3 or 4x4 slab, 4+ layers
  - Relax top 2 layers
  - ENCUT 400+ eV
  - Test all sites (ontop, bridge, fcc, hcp)

Key formula:
  E_ads = E(slab+adsorbate) - E(slab) - E(adsorbate)

Next: Try 10_reactions/ for reaction energetics.