cp2k-geo-opt
CP2K geometry optimization. Handles bulk, slab, and molecular systems with GPW method. Efficient for large systems (200+ atoms).
SKILL.md.claude/skills/cp2k-relax/SKILL.md
GitHubCP2K Geometry Optimization
Set up and submit CP2K geometry optimizations using the Gaussian and Plane-Wave (GPW) method. CP2K is the preferred code for systems larger than ~200 atoms where VASP becomes memory-limited.
Scenario 1: Bulk Optimization
Full cell and ionic relaxation for periodic bulk systems.
from catgo.workflow import Workflow
wf = Workflow("CP2K bulk MgO")
struct = wf.add_task("structure_input", structure=bulk_json)
opt = wf.add_task("geo_opt",
structure=struct.output.structure,
software="cp2k",
cell_opt=True, # Relax cell + ions (like ISIF=3 in VASP)
cutoff=600, # Ry
rel_cutoff=60, # Ry
basis_set="DZVP-MOLOPT-SR-GTH",
xc_functional="PBE",
max_iter=200, # Max geo_opt steps
eps_geo=3e-4, # Force convergence (Hartree/Bohr)
system_name="bulk_MgO")
wf.submit()
MCP equivalent:
catgo_workflow_v2(action="create", params={"name": "CP2K bulk MgO"})
catgo_workflow_v2(action="add_task", params={
"workflow_id": "wf_xxx",
"task_type": "structure_input",
"structure": "<bulk_json>"
})
catgo_workflow_v2(action="add_task", params={
"workflow_id": "wf_xxx",
"task_type": "geo_opt",
"software": "cp2k",
"structure": "{{t_001.output.structure}}",
"cell_opt": true,
"cutoff": 600,
"basis_set": "DZVP-MOLOPT-SR-GTH",
"system_name": "bulk_MgO"
})
catgo_workflow_v2(action="submit", params={"workflow_id": "wf_xxx"})
Scenario 2: Slab Optimization
Fixed cell, ionic relaxation with frozen bottom layers. Analogous to VASP ISIF=2.
wf = Workflow("CP2K TiO2 slab")
struct = wf.add_task("structure_input", structure=slab_json)
opt = wf.add_task("geo_opt",
structure=struct.output.structure,
software="cp2k",
cell_opt=False, # Fix cell (slab)
cutoff=600,
basis_set="DZVP-MOLOPT-SR-GTH",
freeze_layers=2, # Freeze bottom 2 layers
poisson_solver="MT", # Martyna-Tuckerman for slab geometry
system_name="TiO2_slab")
wf.submit()
Slab-specific settings:
cell_opt=False— mandatory for slabs (equivalent to ISIF=2 in VASP)freeze_layers=2— freeze bottom layers to mimic bulkpoisson_solver="MT"— Martyna-Tuckerman solver handles the vacuum correctly for 2D-periodic systems. Use "PERIODIC" for bulk (3D-periodic) and "MT" or "WAVELET" for slabs
Scenario 3: Adsorbate on Slab
Same as slab, with adsorbate atoms free to relax:
opt = wf.add_task("geo_opt",
structure=adsorbate_slab_json,
software="cp2k",
cell_opt=False,
freeze_layers=2,
cutoff=600,
vdw_method="DFTD3", # Dispersion for adsorption
poisson_solver="MT",
system_name="*OH_on_TiO2")
Scenario 4: Large System (500+ atoms)
CP2K's GPW method with OT (Orbital Transformation) SCF solver scales linearly for large systems:
opt = wf.add_task("geo_opt",
structure=large_system_json,
software="cp2k",
cutoff=400, # Lower cutoff acceptable for screening
basis_set="SZV-MOLOPT-SR-GTH", # Minimal basis for speed
ot_minimizer="DIIS", # OT method for large systems
ot_preconditioner="FULL_ALL",
eps_scf=1e-5,
system_name="large_system")
OT vs diagonalization:
- OT: O(N) scaling, no HOMO-LUMO gap requirement, default for > 100 atoms
- Diagonalization: O(N^3) scaling, needed for metallic systems (zero gap)
For metals: OT does not work for metallic systems (zero band gap). Use Fermi-Dirac smearing with diagonalization:
opt = wf.add_task("geo_opt",
structure=metal_json,
software="cp2k",
scf_method="diag", # Standard diagonalization
smearing_method="FERMI_DIRAC",
electronic_temperature=300, # K
system_name="metal")
Key Parameters
| Parameter | Default | Purpose |
|---|---|---|
| cutoff | 600 Ry | PW cutoff for density grid |
| rel_cutoff | 60 Ry | Multi-grid relative cutoff |
| basis_set | DZVP-MOLOPT-SR-GTH | Gaussian basis set |
| xc_functional | PBE | Exchange-correlation functional |
| max_iter | 200 | Max geometry optimization steps |
| eps_geo | 3e-4 | Force convergence (Hartree/Bohr, ~ 0.015 eV/A) |
| eps_scf | 1e-6 | SCF convergence (Hartree) |
| cell_opt | False | Whether to optimize cell parameters |
| freeze_layers | 0 | Number of bottom layers to freeze |
| vdw_method | None | Dispersion correction ("DFTD3", "DFTD3(BJ)") |
| poisson_solver | PERIODIC | Poisson solver ("PERIODIC", "MT", "WAVELET") |
Convergence Monitoring
catgo_workflow_v2(action="status", params={"workflow_id": "wf_xxx"})
catgo_analyze(action="convergence", params={"task_id": "t_opt"})
# Returns: energy vs step, max force vs step
catgo_analyze(action="forces", params={"task_id": "t_opt"})
Chain: CP2K Optimization then Frequency
wf = Workflow("CP2K opt + freq")
struct = wf.add_task("structure_input", structure=slab_oh_json)
opt = wf.add_task("geo_opt", structure=struct.output.structure,
software="cp2k", freeze_layers=2,
system_name="*OH")
frq = wf.add_task("freq", structure=opt.output.structure,
software="cp2k",
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()
Output
The geo_opt task produces:
output.structure— optimized structure (pymatgen dict as JSON string)output.energy— total DFT energy in eV
Troubleshooting
| Problem | Fix |
|---|---|
| SCF not converging | Use OT method, increase scf_max_iter, reduce mixing |
| Energy oscillations | Increase cutoff (try 800 Ry), check rel_cutoff |
| Forces not converging | Loosen eps_geo, increase max_iter |
| OT fails for metal | Switch to diagonalization with Fermi smearing |
| Memory error | Reduce cutoff, use SZV basis, increase nodes |
| Missing basis for element | Check CP2K basis set library, may need to download |
| Poisson solver error for slab | Use poisson_solver="MT" instead of "PERIODIC" |
Unit Conversions
CP2K uses atomic units internally. CatGo converts automatically, but for reference:
- 1 Hartree = 27.2114 eV
- 1 Bohr = 0.529177 A
- Force: 1 Ha/Bohr = 51.422 eV/A
- eps_geo=3e-4 Ha/Bohr corresponds to ~0.015 eV/A (comparable to VASP EDIFFG=-0.02)