VASP Geometry Optimization

Set up and submit VASP geometry optimizations. Three main scenarios with different parameter requirements.

Discussion Checkpoints

🔴 Must discuss with user:

• Functional (METAGGA/GGA/+U) — PBE vs SCAN vs PBE+U fundamentally changes energetics; wrong functional invalidates the entire study
• ISPIN — must be 2 for magnetic systems (Fe, Co, Ni, Mn oxides, NRR substrates); default ISPIN=1 gives wrong energies for magnetic materials
• Structure source — bulk from Materials Project vs user-uploaded CIF vs previous optimization; wrong starting structure wastes all compute

🟡 Recommend confirming:

• ENCUT (default: 520) — increase to 600+ for accurate equation of state or when comparing across different compositions
• EDIFFG (default: -0.02 eV/A) — tighten to -0.01 for frequency calculations downstream; loosen to -0.05 for quick screening
• k-points — must be converged for the system; small unit cells need denser meshes
• Selective dynamics / frozen layers (default: freeze_layers=2 for slabs) — adjust based on slab thickness and whether subsurface relaxation matters
• ISIF (default: 2) — must be 3 for bulk relaxation, 2 for slabs; wrong ISIF is a common mistake

🟢 Safe defaults:

• EDIFF = 1E-5
• ISMEAR = 0, SIGMA = 0.05
• NSW = 200
• IBRION = 2 (conjugate gradient)
• PREC = Accurate
• NCORE = 4

Scenario 1: Bulk Relaxation

Full cell + ionic relaxation. Use ISIF=3 to allow cell shape and volume to change.

from catgo.workflow import Workflow
from catgo.workflow.builtins import geo_opt

wf = Workflow("Bulk TiO2 relaxation")
struct = wf.add_task("structure_input", structure=bulk_json)
opt = wf.add_task(geo_opt, structure=struct.output.structure,
                  ISIF=3,        # Relax cell shape + volume + ions
                  EDIFFG=-0.02,  # Force convergence (eV/A)
                  system_name="bulk_TiO2")
wf.submit()

Key parameters:

• ISIF=3 — relax ions + cell shape + cell volume
• EDIFFG=-0.02 — converge when max force < 0.02 eV/A (negative = force criterion)
• NSW=200 — max ionic steps (default, usually converges in 50-100)

MCP equivalent:

catgo_workflow_engine(action="create", params={"name": "Bulk TiO2"})

catgo_workflow_engine(action="add_task", params={
  "workflow_id": "wf_xxx",
  "task_type": "structure_input",
  "name": "input",
  "structure": "<bulk_json>"
})

catgo_workflow_engine(action="add_task", params={
  "workflow_id": "wf_xxx",
  "task_type": "geo_opt",
  "software": "vasp",
  "structure": "{{t_001.output.structure}}",
  "ISIF": 3,
  "system_name": "bulk_TiO2"
})

catgo_workflow_engine(action="submit", params={"workflow_id": "wf_xxx"})

Scenario 2: Slab Relaxation (Clean Surface)

Fixed cell, relax only ions. Bottom layers frozen to mimic bulk.

wf = Workflow("RuO2(110) slab")
struct = wf.add_task("structure_input", structure=slab_json)
opt = wf.add_task(geo_opt, structure=struct.output.structure,
                  ISIF=2,              # Fix cell, relax ions only
                  selective_dynamics=True,
                  freeze_layers=2,     # Freeze bottom 2 layers
                  system_name="clean_slab")
wf.submit()

Key parameters:

• ISIF=2 — MANDATORY for slabs. Fixes cell shape and volume
• freeze_layers=2 — freeze bottom N layers (sorted by z-coordinate)
• selective_dynamics=True — enable per-atom freeze in POSCAR
• Vacuum: ensure >= 15 A in z-direction to avoid periodic image interaction

ISIF reference:

ISIF	Ions	Cell shape	Cell volume	Use case
2	Yes	No	No	Slabs, adsorbates
3	Yes	Yes	Yes	Bulk relaxation
4	Yes	Yes	No	Bulk at fixed volume

Scenario 3: Adsorbate on Slab

Same as slab relaxation, but the structure has an adsorbate. Adsorbate atoms are always free.

wf = Workflow("OH on RuO2(110)")
struct = wf.add_task("structure_input", structure=adsorbate_slab_json)
opt = wf.add_task(geo_opt, structure=struct.output.structure,
                  ISIF=2,
                  selective_dynamics=True,
                  freeze_layers=2,
                  EDIFFG=-0.02,
                  system_name="*OH")
wf.submit()

Guidelines for adsorbates:

• freeze_layers only affects the slab — adsorbate atoms above the surface are always relaxed
• Use system_name="*OH" convention (asterisk = adsorbed species)
• For weak adsorbates (CO2, H2O physisorption), add IVDW=11 for DFT-D3

Confirmation Gate

By default, HPC tasks (including VASP relaxation) pause at PENDING_REVIEW after local preprocessing completes. This lets users verify the structure, frozen layers, and VASP parameters in the task detail panel before committing HPC resources. Click "Confirm & Submit" in the UI, or use wf.submit(auto_submit=True) to bypass the gate.

Visual Verification Steps

Before submitting, verify the structure in the viewer:

# Step 1: Check current structure in viewer
catgo_view(action="get_state")
# Verify: correct composition, reasonable cell parameters, vacuum > 15 A for slabs

# Step 2: For slabs, verify frozen atoms
catgo_view(action="get_state")
# Check that bottom-layer atoms will be frozen

# Step 3: After optimization completes, push result to viewer
catgo_workflow_engine(action="get_result", params={"task_id": "t_opt"})
catgo_view(action="push", params={"structure": "<optimized_structure_json>"})

Convergence Monitoring

# Check if optimization is converged
catgo_analyze(action="convergence", params={"task_id": "t_opt"})
# Returns: energy vs step, max force vs step, whether converged

# Check remaining forces
catgo_analyze(action="forces", params={"task_id": "t_opt"})
# Returns: per-atom forces, max force, force on each species

Common Chain: Relaxation then Frequency

For thermodynamic properties (Gibbs energy, ZPE), chain relaxation with frequency:

from catgo.workflow.builtins import geo_opt, freq, gibbs_energy

wf = Workflow("OH adsorption thermodynamics")
struct = wf.add_task("structure_input", structure=slab_oh_json)

opt = wf.add_task(geo_opt, structure=struct.output.structure,
                  ISIF=2, freeze_layers=2, system_name="*OH")

frq = wf.add_task(freq, structure=opt.output.structure,
                  freeze_mode="layers", freeze_layers=4,
                  system_name="*OH")

gib = wf.add_task(gibbs_energy,
                  energy=opt.output.energy,
                  frequencies=frq.output.frequencies,
                  phase="adsorbed", system_name="*OH")

wf.submit()

Troubleshooting

Problem	Fix
Forces not converging	Increase NSW (e.g., 400), or loosen EDIFFG to -0.03
Atoms escaping into vacuum	Check initial adsorbate placement, reduce POTIM to 0.3
SCF not converging	Switch ALGO=All, increase NELM=400, try AMIX=0.1
Cell shape changing for slab	Ensure ISIF=2, not ISIF=3
Wrong energy (magnetic system)	Set ISPIN=2, provide MAGMOM

Parameter Defaults (inherited from config)

These are applied automatically unless overridden:

• ENCUT=520, EDIFF=1e-5, PREC=Accurate
• IBRION=2 (conjugate gradient), NSW=200
• EDIFFG=-0.02, ISIF=2
• ISMEAR=0, SIGMA=0.05 (Gaussian smearing)
• NCORE=4, LREAL=Auto