CP2K notes

11 minute read

Published:

Content

Brief introduction of CP2K

another open sourced qunatum mechanical material science software. Slightly different to ABINIT and Quantum ESPRESSO:

  • very quick compared to ABINIT and QE
    • in particular computing large system (large period)
  • DFT using Gaussian and plane waves approaches GPW and GAPW (Augmented)
    • centered Gaussian orbitals
      • only few basis functions per atom are needed
      • foundation of wave function representation and Kohn-Sham matrix
    • plane waves as regular grids
      • electronic density
    • Kohn-Sham matrix and density matrix become sparse, allowing linear scaling method to perform density matrix optimization
    • memory friendly, efficient

Tutorial

official

Installation

found github and download source .tar.bz2

  • tar -xvf cp2k.tar.gz2

    1. Requirement

  • CMake
    • sudo apt install cmake
  • C, C++, and Fortran compilers
  • DBCSR
    • cd /cp2k/exts/dbcsr
    • mkdir build && cd build
    • cmake ..
    • make
    • sudo make install
  • OpenMP
  • BLAS and LAPACK
  • MPI
  • ScaLAPACK
    • download from official website
    • tar -zxvf scalapack.tar.gz
    • mkdir build && cd build
    • cmake .. -DCMAKE_Fortran_FLAGS=-fPIC -DCMAKE_C_FLAGS=-fPIC
    • cmake --build .
  • fftw3
    • sudo apt-get install libfftw3-dev

2. Using CMake

official

cd cp2k/
mkdir build/
cmake -S . -B build             # -S <SOURCE_DIR> -B <BUILD_DIR>
cmake --build build -- -j 16    # -build <BUILD_DIR> -- <BUILD_TOOL_OPTIONS>

not official

work and successful on 2025Jun15

cd cp2k/
mkdir build/
cd build/
cmake .. -DGMX_BUILD_OWN_FFTW=ON -DCMAKE_PREFIX_PATH="/home/yongnanli08/Software/scalapack-2.2.2/build/lib/"
cmake --build ..

remember to add path into $PATH

Use

Input files

.inp

main input file

&GLOBAL
  PROJECT Si_bulk8      ! project_name, also the prefix of output files
  RUN_TYPE ENERGY_FORCE ! static energy and force calculation
                        ! can be ENERGY: only energy calculation
                        ! GEO_OPT
  PRINT_LEVEL MEDIUM    ! verbosity
                        ! can be LOW
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep      ! GPW Gaussian and planewaves
  &SUBSYS
    &KIND Si
      ELEMENT	Si
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-PADE-q4
    &END KIND
    &CELL
      A     5.430697500    0.000000000    0.000000000
      B     0.000000000    5.430697500    0.000000000
      C     0.000000000    0.000000000    5.430697500
    &END CELL
    &COORD
      Si    0.000000000    0.000000000    0.000000000
      Si    0.000000000    2.715348700    2.715348700
      Si    2.715348700    2.715348700    0.000000000
      Si    2.715348700    0.000000000    2.715348700
      Si    4.073023100    1.357674400    4.073023100
      Si    1.357674400    1.357674400    1.357674400
      Si    1.357674400    4.073023100    4.073023100
      Si    4.073023100    4.073023100    1.357674400
      ! cartesian coordinates in Angstroms
      ! SCALED .TRUE. ! can use this to change to frational coordinates
    &END COORD
  &END SUBSYS
  &DFT
    BASIS_SET_FILE_NAME  BASIS_SET
    POTENTIAL_FILE_NAME  GTH_POTENTIALS
    &QS
    ! general control parameters used by QUICKSTEP
      EPS_DEFAULT 1.0E-10
    &END QS
    &MGRID
    ! integration grid used by QUICKSTEP
    ! QUICKSTEP multi-grid method for Gaussian functions
      NGRIDS 4          ! 4 levels of multi-grids
      CUTOFF 300        ! Ry, planewave cutoff of the finest grid
      REL_CUTOFF 60     ! Ry, grid spacing underneath Gaussian functions, finer than equivalent planewave cutoff 
    &END MGRID
    &XC
      &XC_FUNCTIONAL PADE ! consistent to the basis set and pseudopotential
      &END XC_FUNCTIONAL
    &END XC
    &SCF
      SCF_GUESS ATOMIC  ! how initial trial electron density function is to be generated, overlapping of atomic charge densities
      EPS_SCF 1.0E-7    ! tolerance of the charge density residual
      MAX_SCF 300       ! max scf steps
      ADDED_MOS 10      ! number of lowest empty molecular orbitals not to omit
      &SMEAR ON
        METHOD FERMI_DIRAC
        ELECTRONIC_TEMPERATURE [K] 300
      &END SMEAR
      &DIAGONALIZATION  ON
                        ! traditional diagonalization for ground state Kohn-Sham energy and electron density
                        ! alternative: orbital transform (OT) method
        ALGORITHM STANDARD ! lapack/scalapack for diagonalizatioin
      &END DIAGONALIZATION
      &MIXING  T
                        ! charge mixing in scf, only apply to traiditional diagonalization method
        METHOD BROYDEN_MIXING
        ! DIRECT_P_MIXING KERKER_MIXING PULAY_MIXING BROYDEN_MIXING MULTISECANT_MIXING
        ALPHA 0.4       ! mixing parameter
        NBROYDEN 8      ! number of histories to be used in mixing
      &END MIXING
    &END SCF
  &END DFT
  &PRINT
    &FORCES ON
    &END FORCES
  &END PRINT
&END FORCE_EVAL

BASIS_SET

basis sets

  • may be found in /cp2k/data
  • can change name and modify in .inp (BASIS_SET_FILE_NAME)

    GTH_POTENTIALS

    pseudopotentials

  • may be found in /cp2k/data
  • can change name and modify in .inp (BPOTENTIAL_FILE_NAME)

Run

mpirun cp2k.psmp .inp > .out &

Convergence

energy convergence: 1meV/atom, for large system 10meV/atom

CUTOFF

with PRINT_LEVEL MEDIUM

  • other than total energy convergence (~<10e-8 Ha?), also check multi-grid in .out
 -------------------------------------------------------------------------------
 ----                             MULTIGRID INFO                            ----
 -------------------------------------------------------------------------------
 count for grid        1:           3240          cutoff [a.u.]          150.00
 count for grid        2:          12976          cutoff [a.u.]           50.00
 count for grid        3:          13384          cutoff [a.u.]           16.67
 count for grid        4:           4488          cutoff [a.u.]            5.56
 total gridlevel count  :          34088

from 1 to 4, from fine to coarse

increase of CUTOFF, number of Gaussians assigned to finest grids decrease:

# Cutoff (Ry) | Total Energy (Ha) | NG on grid 1 | NG on grid 2 | NG on grid 3 | NG on grid 4
     50.00   -32.3795329864    5048    5432      16       0
    100.00   -32.3804557631    2720    5000    2760      16
    150.00   -32.3804554850    2032    3016    5432      16
    200.00   -32.3804554982    1880    2472    3384    2760
    250.00   -32.3804554859     264    4088    3384    2760
    300.00   -32.3804554843     264    2456    5000    2776
    350.00   -32.3804554846      56    1976    5688    2776
    400.00   -32.3804554851      56    1976    3016    5448
    450.00   -32.3804554851       0    2032    3016    5448
    500.00   -32.3804554850       0    2032    3016    5448
  • large CUTOFF lead to slow convergence in energy
  • reasonable distribution:
    • it is the lowest cutoff energy where the finest grid level is used, but at the same time with the majority of the Gaussians on the coarser grids.
      • above example, 250 is ok

REL_CUTOFF

similar to CUTOFF

# Rel Cutoff (Ry) | Total Energy (Ha) | NG on grid 1 | NG on grid 2 | NG on grid 3 | NG on grid 4
     10.00   -32.3902980020       0       0    2032    8464
     20.00   -32.3816384686       0     264    4088    6144
     30.00   -32.3805115576       0    2032    3016    5448
     40.00   -32.3805116025      56    1976    3016    5448
     50.00   -32.3804555002     264    2456    5000    2776
     60.00   -32.3804554859     264    4088    3384    2760
     70.00   -32.3804554859    1880    2472    3384    2760
     80.00   -32.3804554859    1880    2472    3384    2760
     90.00   -32.3804554848    2032    3016    5432      16
    100.00   -32.3804554848    2032    3016    5432      16
  • increase REL_CUTOFF, more Gaussians get mapped onto the finer grids
    • it is the lowest cutoff energy where the finest grid level is used, but at the same time with the majority of the Gaussians on the coarser grids.
      • above example, 60 is ok

Output files

-RESTART.wfn

final Kohn-Sham wavefunctions

-RESTART.wfn.bak-n

Kohn-Sham wavefunctions from the n previous SCF steps

  • 1
    • last SCF step

can start new SCF from the previous calculated wavefunctions by

  • SCF
    • SCF_GUESS RESTART
    • share the same PROJECT_NAME

Geometrical Optimization

other than GLOBAL and FORCE_EVAL, geometrical optimization needs an additional section MOTION

&MOTION
  &GEO_OPT
    TYPE MINIMIZATION
    MAX_DR    1.0E-03   ! convergence criterion, max geometry change in bohr
    MAX_FORCE 1.0E-03   ! convergence criterion, max force in hartree/bohr
    RMS_DR    1.0E-03   ! convergence criterion, RMS geometry change in bohr
    RMS_FORCE 1.0E-03   ! convergence criterion, RMS force in hartree/bohr
    MAX_ITER 200
    OPTIMIZER CG        ! BFGS LBFGS CG
    &CG
      MAX_STEEP_STEPS  0
      RESTART_LIMIT 9.0E-01
    &END CG
  &END GEO_OPT
  &CONSTRAINT
    &FIXED_ATOMS
    ! fix particular atoms
      COMPONENTS_TO_FIX XYZ
      LIST 1            ! indice in COORD
    &END FIXED_ATOMS
  &END CONSTRAINT
&END MOTION

for the example in website, not use PULAY_MIXING, otherwise abort

Example

Basis set & Pseudopotential

use this web to choose and download basis set & pseudopotential

Multiwfn

It seems that using Multiwfn can enhance the usage and efficiency of cp2k programming

Installation

Multiwfn
<enter input file name, can be qe, xyz, cif> (Press Enter can use GUI to choose input file)

Create cp2k input file

after showing a list of Main function menu

cp2k (to create cp2k input file)
<cp2k input file name>

before creating, generate basis set and pseudopotential first like this

Common setting

  • DFT
    • PBE, PBErev for surface, PBEsol for metals/oxides, PBE0 better but expensive
    • no need to define in Multiwfn cuz need to change after
  • basis set and pseudopotential
    • at least TZVP-MOLOPT-GTH (3) or TZV2P-MOLOPT-GTH (4), can try all electrons later
    • no need to define in Multiwfc cuz need to change after
  • diagonalization
    • or orbital transform (OT) later, can try later
      • OT is for large systems and large basis sets, 1k atoms (20k-30k basis functions)
  • mixing
    • Broyden slightly faster than Pulay
  • smearing
    • semiconductor not very big deal but important in metal
  • self-consistent continuum solvation (SCCS)
    • analytical expressions of the local solute-solvent interface functions determine the interface function and dielectric function values at a given real space position
  • k point
    • mostly primitive cell need test
    • large system can try gamma only or slightly increase and test
    • larger supercell, sparser k point grid
  • charge
    • other properties
  • Basis set and Pseudopotential
    • check cp2k-basis
    • select and download corresponding basis set & pseudopotential, than change inp with provided code

basis set, pseudopotential input geneartion

cp2k-basis

this website can choose and generate input of

  • basis set
  • pseduopotential
  • input file in cp2k

better to replace the input in Multiwfn, some of selection in Multiwfn cannot be directly found in these websites.

create 3 files:

  • INPUT
    • put in .inp
  • POTENTIAL
  • BASIS

remember to check

  • BASIS_SET_FILE_NAME BASIS
  • POTENTIAL_FILE_NAME POTENTIAL in Multiwfn generated .inp flie

Reference

basis set

in Gaussian but cp2k can take it as 对于非双杂化泛函的DFT计算,根据精度要求,一般问题基组选用建议如下。符号对应于同一档内基组尺寸关系。 垂死挣扎级别:STO-3G (极小基) 水深火热级别:3-21G (最烂2-zeta基组) 最低可接受级别: def2-SV(P) ≈ 6-31G* < def2-SVP ≈ 6-31G** ≈ pcseg-1 (2-zeta基组) 不错级别:6-311G** < def-TZVP (一般3-zeta基组) 理想级别: def2-TZVP < def2-TZVPP ≈ pcseg-2 (高档3-zeta基组) 无敌: def2-QZVP ≈ pcseg-3 (4-zeta基组。给DFT用很浪费)

对于后HF、双杂化泛函、多参考方法计算,根据精度要求,一般问题基组选用建议如下。 最低可接受级别: def-TZVP (一般3-zeta基组) 较好级别: cc-pVTZ ≈ def2-TZVPP (高档3-zeta基组) 高精度计算: cc-pVQZ ≈ def2-QZVPP (4-zeta基组) 无敌: cc-pV5Z (5-zeta基组。很浪费且极昂贵,一般通过TZ→QZ CBS外推来达到这个档次)

  • structural optimization (geometrical)
    • 低至6-31G*高至def-TZVP就够

DFT

  • PBE (Perdew-Burke-Ernzerhof) – Type: GGA (generalized gradient approximation) – Construction: parameter-free, enforces known exact constraints on exchange-correlation energy – Strengths: good balance of accuracy vs. cost for molecules and solids – Weaknesses: tends to overbind adsorbates on surfaces and slightly overestimate lattice constants

  • revPBE (revised PBE) – Type: GGA – What’s changed: exchange enhancement factor is “softened” (rises more steeply with density gradient) – Why it matters: gives more accurate chemisorption and surface energetics than PBE

  • PBEsol (“PBE for solids”) – Type: GGA – Re-tuning: restores the gradient-expansion approximation for slowly varying densities – Use case: equilibrium lattice constants, bulk moduli and surface energies of metals/oxides come out much closer to experiment than with PBE

  • PBE0 (aka PBE1PBE) – Type: hybrid GGA – Mixing: 25 % exact (Hartree–Fock) exchange + 75 % PBE exchange, with full PBE correlation – Trade-off: better atomization energies, reaction barriers and band-gaps than PBE, at ≈3× the computational cost of pure GGAs

  • B3LYP – Type: hybrid GGA – Mixing: Becke’s three-parameter fit of exact exchange + LYP correlation – Characteristics: heavily empirically parametrized for thermochemistry and kinetics of organic molecules – Limitations: suffers from delocalization error in charge-transfer, underestimates solid-state band-gaps and lattice constants

  • HSE06 (Heyd–Scuseria–Ernzerhof) – Type: screened hybrid – Split: short-range HF exchange (≈25 %) + long-range PBE exchange, full PBE correlation – Advantages: accurate band-gaps in semiconductors/insulators, much faster convergence in periodic solids vs. global hybrids like PBE0

— Beyond these, many researchers explore: – Dispersion-corrected GGAs (e.g., PBE-D3) for van der Waals systems – Meta-GGAs (e.g., SCAN) that incorporate kinetic-energy density – Range-separated hybrids (e.g., ωB97X) that tune the HF fraction with interelectronic distance