14 - Van der Waals Corrections

This tutorial demonstrates dispersion-corrected DFT calculations using graphite as an example, where van der Waals interactions are essential.

Note: Using reduced settings for demonstration. Production calculations would use higher ENCUT and denser k-points.

import numpy as np
from ase import Atoms

from vasp import Vasp
from vasp.parameters import get_vdw_params

print("=" * 60)
print("Van der Waals Corrections")
print("=" * 60)
print()

============================================================
Van der Waals Corrections
============================================================

Part 1: Graphite structure

print("Part 1: Graphite structure")
print("-" * 40)
print()

# Create graphite (AB stacking)
# Hexagonal cell with 4 atoms
a = 2.46   # In-plane lattice constant
c = 6.71   # Interlayer spacing (experimental)

graphite = Atoms(
    symbols=['C', 'C', 'C', 'C'],
    scaled_positions=[
        [0.0, 0.0, 0.0],
        [1/3, 2/3, 0.0],
        [0.0, 0.0, 0.5],
        [2/3, 1/3, 0.5],
    ],
    cell=[
        [a, 0, 0],
        [-a/2, a*np.sqrt(3)/2, 0],
        [0, 0, c],
    ],
    pbc=True,
)

print("Structure: Graphite (AB stacking)")
print(f"  In-plane lattice: a = {a} A")
print(f"  Interlayer spacing: c/2 = {c/2:.2f} A")
print("  Experimental interlayer: 3.35 A")
print()

Part 1: Graphite structure
----------------------------------------

Structure: Graphite (AB stacking)
  In-plane lattice: a = 2.46 A
  Interlayer spacing: c/2 = 3.35 A
  Experimental interlayer: 3.35 A

Part 2: Standard PBE (no vdW)

print("Part 2: Standard PBE (no dispersion)")
print("-" * 40)
print()

calc_pbe = Vasp(
    label='results/vdw/graphite_pbe',
    atoms=graphite.copy(),
    xc='PBE',
    encut=400,  # Carbon needs decent cutoff
    kpts=(8, 8, 4),  # Reduced for speed
    ismear=0,
    sigma=0.05,
)

e_pbe = calc_pbe.potential_energy
print(f"PBE energy: {e_pbe:.6f} eV")
print()

# Scan interlayer distance with PBE (fewer points)
print("Interlayer binding scan (PBE):")
c_values = [6.0, 6.71, 8.0]
pbe_energies = []

for c_test in c_values:
    graphite_test = graphite.copy()
    cell = graphite_test.get_cell()
    cell[2, 2] = c_test
    graphite_test.set_cell(cell, scale_atoms=True)

    calc = Vasp(
        label=f'results/vdw/graphite_pbe_c{c_test:.1f}',
        atoms=graphite_test,
        xc='PBE',
        encut=400,
        kpts=(8, 8, 4),
        ismear=0,
        sigma=0.05,
    )
    e = calc.potential_energy
    pbe_energies.append(e)
    print(f"  c = {c_test:.2f} A: E = {e:.4f} eV")

print()
print("PBE shows no binding minimum (energy decreases with c)")
print("This is WRONG - graphite layers should bind!")
print()

Part 2: Standard PBE (no dispersion)
----------------------------------------

PBE energy: -36.872474 eV

Interlayer binding scan (PBE):
  c = 6.00 A: E = -36.7092 eV
  c = 6.71 A: E = -36.8725 eV
  c = 8.00 A: E = -36.9350 eV

PBE shows no binding minimum (energy decreases with c)
This is WRONG - graphite layers should bind!

Part 3: D3-BJ correction

print("Part 3: PBE-D3(BJ) (Grimme's D3 with BJ damping)")
print("-" * 40)
print()

d3bj_params = get_vdw_params('d3bj')
print("D3-BJ parameters:")
print("  IVDW = 12 (D3 with BJ damping)")
print("  Empirical C6 coefficients with geometry-dependent terms")
print()

d3bj_energies = []
for c_test in c_values:
    graphite_test = graphite.copy()
    cell = graphite_test.get_cell()
    cell[2, 2] = c_test
    graphite_test.set_cell(cell, scale_atoms=True)

    calc = Vasp(
        label=f'results/vdw/graphite_d3bj_c{c_test:.1f}',
        atoms=graphite_test,
        xc='PBE',
        encut=400,
        kpts=(8, 8, 4),
        ismear=0,
        sigma=0.05,
        **d3bj_params,
    )
    e = calc.potential_energy
    d3bj_energies.append(e)
    print(f"  c = {c_test:.2f} A: E = {e:.4f} eV")

# Find minimum
min_idx = np.argmin(d3bj_energies)
print()
print(f"D3-BJ minimum near c = {c_values[min_idx]:.2f} A")
print("Experimental: c = 6.71 A (interlayer = 3.35 A)")
print()

Part 3: PBE-D3(BJ) (Grimme's D3 with BJ damping)
----------------------------------------

D3-BJ parameters:
  IVDW = 12 (D3 with BJ damping)
  Empirical C6 coefficients with geometry-dependent terms

  c = 6.00 A: E = -37.4198 eV
  c = 6.71 A: E = -37.4981 eV
  c = 8.00 A: E = -37.4444 eV

D3-BJ minimum near c = 6.71 A
Experimental: c = 6.71 A (interlayer = 3.35 A)

Part 4: Compare dispersion methods

print("Part 4: Comparison of vdW methods")
print("-" * 40)
print()

# Test at experimental c
graphite_exp = graphite.copy()

# Test just three methods for speed
methods = [
    ('None', {}),
    ('d3bj', get_vdw_params('d3bj')),
    ('ts', get_vdw_params('ts')),
]

print("Comparison at experimental geometry (c = 6.71 A):")
print()
print(f"  {'Method':<10} {'Energy (eV)':<15} {'IVDW':<10}")
print("  " + "-" * 35)

for name, params in methods:
    calc = Vasp(
        label=f'results/vdw/graphite_{name}',
        atoms=graphite_exp.copy(),
        xc='PBE',
        encut=400,
        kpts=(8, 8, 4),
        ismear=0,
        sigma=0.05,
        **params,
    )
    e = calc.potential_energy
    ivdw = params.get('ivdw', 0)
    print(f"  {name:<10} {e:<15.4f} {ivdw:<10}")

print()

Part 4: Comparison of vdW methods
----------------------------------------

Comparison at experimental geometry (c = 6.71 A):

  Method     Energy (eV)     IVDW      
  -----------------------------------
  None       -36.8725        0         
  d3bj       -37.4981        12        
  ts         -37.4459        2         

Part 5: Binding energy calculation

print("Part 5: Graphite binding energy")
print("-" * 40)
print()

# Calculate isolated graphene layer with PBE (no vdW)
graphene = Atoms(
    symbols=['C', 'C'],
    scaled_positions=[
        [0.0, 0.0, 0.5],
        [1/3, 2/3, 0.5],
    ],
    cell=[
        [a, 0, 0],
        [-a/2, a*np.sqrt(3)/2, 0],
        [0, 0, 15.0],  # Large vacuum
    ],
    pbc=True,
)

calc_graphene_pbe = Vasp(
    label='results/vdw/graphene_pbe',
    atoms=graphene.copy(),
    xc='PBE',
    encut=400,
    kpts=(8, 8, 1),
    ismear=0,
    sigma=0.05,
)
e_graphene_pbe = calc_graphene_pbe.potential_energy

# Calculate graphene with D3-BJ
calc_graphene_d3 = Vasp(
    label='results/vdw/graphene_d3bj',
    atoms=graphene.copy(),
    xc='PBE',
    encut=400,
    kpts=(8, 8, 1),
    ismear=0,
    sigma=0.05,
    **d3bj_params,
)
e_graphene_d3 = calc_graphene_d3.potential_energy

# Binding energy per atom (graphite at experimental c = 6.71)
# PBE binding energy
e_graphite_pbe = pbe_energies[1]  # c = 6.71
e_binding_pbe = (e_graphite_pbe - 2 * e_graphene_pbe) / 4

# D3-BJ binding energy
e_graphite_d3 = d3bj_energies[1]  # c = 6.71
e_binding_d3 = (e_graphite_d3 - 2 * e_graphene_d3) / 4

print("Interlayer binding energy (per C atom):")
print()
print(f"  PBE (no vdW): {e_binding_pbe*1000:.1f} meV")
print(f"  PBE-D3(BJ):   {e_binding_d3*1000:.1f} meV")
print("  Experimental: ~52 meV")
print()
print("Note: PBE shows weak/no binding - this is the vdW problem!")
print()

Part 5: Graphite binding energy
----------------------------------------

Interlayer binding energy (per C atom):

  PBE (no vdW): 15.3 meV
  PBE-D3(BJ):   -52.3 meV
  Experimental: ~52 meV

Note: PBE shows weak/no binding - this is the vdW problem!

Summary

print("=" * 60)
print("Summary")
print("=" * 60)
print()

print("Available dispersion corrections:")
print()
print("  Empirical (pairwise):")
print("    - D2: Original Grimme, fixed C6")
print("    - D3: Geometry-dependent C6")
print("    - D3-BJ: D3 with Becke-Johnson damping (recommended)")
print("    - TS: Tkatchenko-Scheffler, uses electron density")
print()
print("  Nonlocal (density-based):")
print("    - vdW-DF family: Physical but expensive")
print("    - SCAN+rVV10: Meta-GGA with nonlocal correlation")
print()

print("Recommendations:")
print("  - D3-BJ: Good accuracy, low cost (default choice)")
print("  - TS-SCS: Better for metals and surfaces")
print("  - vdW-DF: When electron density matters")
print()
print("Demo settings used:")
print("  - ENCUT = 400 eV (adequate for C)")
print("  - Reduced k-points (8x8x4)")
print("  - Fewer c values scanned")
print()
print("Next: Try 15_workflows/ for automated calculation pipelines.")

============================================================
Summary
============================================================

Available dispersion corrections:

  Empirical (pairwise):
    - D2: Original Grimme, fixed C6
    - D3: Geometry-dependent C6
    - D3-BJ: D3 with Becke-Johnson damping (recommended)
    - TS: Tkatchenko-Scheffler, uses electron density

  Nonlocal (density-based):
    - vdW-DF family: Physical but expensive
    - SCAN+rVV10: Meta-GGA with nonlocal correlation

Recommendations:
  - D3-BJ: Good accuracy, low cost (default choice)
  - TS-SCS: Better for metals and surfaces
  - vdW-DF: When electron density matters

Demo settings used:
  - ENCUT = 400 eV (adequate for C)
  - Reduced k-points (8x8x4)
  - Fewer c values scanned

Next: Try 15_workflows/ for automated calculation pipelines.