CP2K 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 bulk
• poisson_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)