Software cheatsheet

The DFT chapters tell you what the equations are. This page tells you which program to run when you actually want a number.

The cheatsheet is organised by use case, not by history. Each section opens with a summary table and follows with a short profile of every code in the category: name, license, language, repository, the DFT methods it implements, the basis sets / boundary conditions it supports, the DFT Notes chapters it is most relevant to, and a “when to use” recommendation.

How to read this page. Jump to the section that matches the boundary condition of your system (molecule vs. solid vs. surface) and its physics (ground-state DFT vs. TDDFT vs. high-accuracy wavefunction vs. excited-state $GW$/BSE). The chapter cross-references in each profile are the most useful pointers: chapter numbers tell you which theoretical machinery the code is exposing to you.

Table of contents

  1. Molecular / cluster codes (Gaussian basis)
  2. Plane-wave / PAW codes
  3. All-electron codes (full potential)
  4. Real-space grid codes
  5. Wavelet codes
  6. Coupled-cluster / high-accuracy wavefunction
  7. Codes for special systems
  8. Visualisation
  9. Workflow managers

1. Molecular / cluster codes (Gaussian basis)

For finite systems (molecules, clusters, bio-organic fragments, transition-metal complexes). Common thread: a Gaussian-type orbital (GTO) basis (chapter 06 §6.3–6.6) centred on the atoms, with the Fock or Kohn–Sham matrix built from analytical electron-repulsion integrals (chapter 03 §3.6.3).

Code License Language Best at
PySCF Apache-2.0 Python (C ext) Open-source QC, embedding, method development
Psi4 LGPL-3.0 Python / C++ Open-source thermochemistry, DFT, CC
ORCA Free for academics C++ Main-group, TM spectroscopy, multi-reference
Q-Chem Commercial C++ / Fortran High-throughput, large molecules, excited states
Gaussian Commercial Fortran General-purpose, broad method coverage

1.1 PySCF

1.2 Psi4

1.3 ORCA

1.4 Q-Chem

1.5 Gaussian

2. Plane-wave / PAW codes

For periodic solids: crystals, surfaces, defects, interfaces, liquids (with PBC). The common thread is a plane-wave basis (chapter 06 §6.7, chapter 07) with Bloch’s theorem for periodicity (chapter 07 §7.3) and pseudopotentials or PAW to handle the core electrons (chapter 08). The plane-wave basis is position-independent, which eliminates the Pulay force (chapter 04 §4.7.3).

Code License Language Best at
Quantum ESPRESSO GPL Fortran Open-source periodic DFT, phonons, TDDFT
VASP Commercial Fortran Industry / academic standard for solids
CASTEP Commercial (UK) Fortran UK academic, phonons, NMR, TDDFT
ABINIT GPL Fortran Open-source, response properties, GW
CP2K GPL Fortran Mixed Gaussian / plane wave, AIMD, large cells
WIEN2k Commercial (ac.) Fortran Full-potential LAPW (see §3)

2.1 Quantum ESPRESSO

2.2 VASP

2.3 CASTEP

2.4 ABINIT

2.5 CP2K

2.6 WIEN2k

See section 3 (all-electron codes). WIEN2k is the most widely-used full-potential LAPW code and is listed in this table because of its prominence in solid-state DFT; it is the entry that breaks the “plane-wave + pseudopotential” mould.

3. All-electron codes (full potential)

Codes that make no pseudopotential approximation — the core electrons are computed explicitly, on the same footing as the valence. The price is computational cost: the core wavefunctions are sharp and require a very high plane-wave cutoff or an alternative basis (Muffin-tin orbitals in LAPW). The reward is the correct treatment of magnetic core states, hyperfine parameters, NMR contact shifts, Mössbauer isomer shifts, spin–orbit coupling, and $f$-electron localisation (chapter 08 §8.13–8.14).

Code License Language Basis Best at
WIEN2k Commercial (ac.) Fortran APW + LAPW + lo Solids, $f$-electron systems, hyperfine
Elk GPL Fortran FP-LAPW Open-source all-electron DFT, response
FLEUR GPL (BSD-style) Fortran FP-LAPW + lo Open-source, magnetism, MAE, NRIXS

3.1 WIEN2k

3.2 Elk

3.3 FLEUR

4. Real-space grid codes

Codes that discretise all of real space (not just the atoms) on a finite-difference, finite-element, or multigrid grid. The basis is position-independent (no atom-centred functions), so the Pulay force vanishes (chapter 04 §4.7.3) and the representation is exact in the limit of an infinitely-fine grid. Real-space grids are the natural setting for time-dependent problems: the wavefunction can be propagated in time with a leap-frog, Crank– Nicolson, or exponential-midpoint scheme, and there is no wrap-around artefact from periodic images of a localised excitation.

Code License Language Best at
Octopus GPL Fortran + C TDDFT, atoms in strong fields, large clusters
GPAW GPL Python + C PAW on a grid, ASE integration, ML hooks
SIESTA GPL Fortran NAO basis, linear-scaling DFT, large systems

4.1 Octopus

4.2 GPAW

4.3 SIESTA

5. Wavelet codes

Wavelet bases — Daubechies wavelets, in particular — are the next step up from finite-difference real-space grids: they combine the locality of the real-space grid (so the Pulay force is zero) with the systematic improvability of a basis-set expansion (chapter 06 §6.14). The convergence of the total energy with the wavelet resolution is exponential, rather than the polynomial convergence of finite-difference / finite-element methods.

5.1 BigDFT

6. Coupled-cluster / high-accuracy wavefunction

Codes that implement wavefunction-based post-HF methods to high accuracy: MP2, MP3, CCSD, CCSD(T), FCI, EOM-CC, explicitly correlated (F12) methods, multireference methods. The hierarchy is described in chapter 02 §2.3; DFT is the cheap alternative, but for small molecules where accuracy is paramount, CCSD(T) in a quadruple-zeta basis is the gold standard (chapter 02 §2.4). The “all in Hartree” cost of CCSD(T) is $K^7$, so these codes are aimed at small-to-medium molecules.

Code License Language Best at
MRCC Free Fortran / C High-accuracy CC, F12, arbitrary-order CC
CFOUR GPL Fortran / C Open-source CC, EOM-CC, MR-CC

6.1 MRCC

6.2 CFOUR

7. Codes for special systems

Codes that occupy a niche. Either they implement a basis that is unique (CRYSTAL with Gaussians for solids, FHI-aims with NAOs of tier-4 quality, exciting with FP-LAPW for excited states), or they are layered on top of a more general code as a post-processor (WEST, Yambo for $G_0W_0$ and BSE).

Code License Language Best at
CRYSTAL Free (ac.) Fortran Gaussian basis for crystals, phonons
FHI-aims Academic (ac.) Fortran All-electron NAO, numerical basis, hybrid DFT
exciting GPL Fortran FP-LAPW excited states, BSE, $GW$
WEST GPL Fortran $G_0W_0$ + BSE on top of PW / PAW codes
Yambo GPL Fortran $G_0W_0$ + BSE on top of PW / PAW codes

7.1 CRYSTAL

7.2 FHI-aims

7.4 WEST

7.5 Yambo

8. Visualisation

Codes for rendering the outputs of electronic- structure calculations. None of these are themselves electronic-structure codes; they read the outputs of the production codes and turn them into pictures. The DFT chapters that produce the data are 04 (charge densities, Kohn–Sham orbitals), 05 (spin densities), 06 (basis-set isosurfaces), 07 (band structures, Fermi surfaces), 08 (PAW augmentation densities), and 10/11 (response properties, optical absorption).

Code License Language Best at
VESTA Freeware C Crystal structures, isosurfaces, charge
OVITO Commercial C++ / Python Molecular-dynamics visualisation
ParaView BSD C++ / Python General 3-D scientific visualisation
XCrySDen GPL C / Tcl/Tk Band structures, Fermi surfaces

8.1 VESTA

8.2 OVITO

8.3 ParaView

8.4 XCrySDen

9. Workflow managers

Codes for organising electronic-structure calculations. A workflow manager accepts a high-level specification (“compute the band structure of MoS₂, then the phonons, then the GW correction”) and turns it into a sequence of job submissions, with dependency tracking, error recovery, and provenance. Workflow managers are not electronic-structure codes; they sit on top of the production codes from sections 1–7. | Code | License | Language | Best at | |:———-|:———|:———–|:————————————————-| | ASE | LGPL | Python | Lightweight scripting, structure handling, IO | | AiiDA | MIT | Python | Workflow + provenance, sharing, materials cloud | | Fireworks | BSD | Python | Lightweight workflow + database | | atomate | MIT | Python | Pre-built workflows (MP-style) for VASP, Q-Chem |

9.1 ASE (Atomic Simulation Environment)

9.2 AiiDA

9.3 Fireworks

9.4 atomate

10. Code cheatsheets (input files)

Companion to the use-case profiles in §§1–9. The profiles tell you which code to reach for; the cheatsheets below tell you what to type into it. Every entry shows a minimal input that produces a sensible result on a trivial system (a hydrogen molecule or a primitive cell), an annotated version of the same input for a slightly harder calculation (a band structure, a spin-polarised system, an MD step), a short list of what to tweak for X, a list of common pitfalls*, and a cross-reference to the chapter where the underlying theory is derived.

The cheatsheet is not a replacement for the manual of any of these codes. It is the first thing to look at when starting a new project, before reading the manual. The “what to tweak” lists are the answers to the questions you will ask on the third day of the project; the “common pitfalls” lists are the answers to the questions you will ask on the third week.

Table of contents

  1. Quantum ESPRESSO — Si band structure
  2. VASP — Si band structure with HSE06
  3. GPAW — Si band structure in Python
  4. SIESTA — Si band structure
  5. CP2K — water in a box and NVT MD
  6. FLEUR — Fe bcc, spin-polarised, with SOC

10.1 Quantum ESPRESSO

Quantum ESPRESSO is a Fortran 2003 plane-wave DFT code (§2.1). The cheatsheet below covers the standard ‘pw.x’ (ground-state) and ‘bands.x’ (post-processing) executables. The input format is name-list Fortran — every section starts with ‘&NAME’ and ends with ‘/`.

10.1.1 Minimal input — H₂ SCF

The smallest input that produces a sensible single-point energy on a hydrogen molecule:

&CONTROL
    calculation = 'scf'         ! 'scf' | 'relax' | 'vc-relax' | 'bands' | 'nscf'
    prefix      = 'h2'          ! base name for output files
    outdir      = './tmp/'      ! scratch directory (must exist)
    pseudo_dir  = './pseudo/'   ! directory containing *.UPF pseudopotentials
    verbosity   = 'low'         ! 'low' | 'high'
/
&SYSTEM
    ibrav = 0                    ! 0 = free-format lattice (ATOMIC_POSITIONS below)
    nat   = 2
    ntyp  = 1
    ecutwfc = 30.0               ! plane-wave cutoff for wavefunctions, Ry
    ecutrho = 240.0              ! charge-density cutoff, Ry
                                 !   >= 4*ecutwfc for NC-PP
                                 !   >= 8-12 * ecutwfc for USPP / PAW
/
&ELECTRONS
    conv_thr = 1.0d-8            ! SCF convergence threshold, Ry
/
ATOMIC_SPECIES
  H  1.008  H.upf                ! element, atomic mass, pseudopotential file
ATOMIC_POSITIONS angstrom
  H  0.0  0.0  0.0
  H  0.0  0.0  0.74              ! bond length 0.74 Å
K_POINTS gamma                    ! Gamma-only: 1 k-point

To run: ‘pw.x < h2.scf.in > h2.scf.out’. The total energy (in Rydberg) is in `h2.scf.out’ after ‘! total energy’.

10.1.2 Annotated example — Si band structure

A band structure of bulk silicon along the standard Γ–X–W–K–Γ–L path. Three calculations: (i) self-consistent ground state on a uniform Monkhorst–Pack mesh; (ii) non-self-consistent calculation along the k-path using the converged charge density; (iii) post-processing with `bands.x’ to extract the band energies.

Step 1 — SCF (si.scf.in):

&CONTROL
    calculation = 'scf'
    prefix      = 'si'
    outdir      = './tmp/'
    pseudo_dir  = './pseudo/'
/
&SYSTEM
    ibrav     = 0
    nat       = 2
    ntyp      = 1
    ecutwfc   = 40.0            ! converged at 40 Ry for Si.pbe-n-rrkjus_psl.0.1.UPF
    ecutrho   = 320.0           ! 8x ecutwfc for ultrasoft / PAW
    occupations = 'fixed'       ! semiconductor: fixed occupations
/
&ELECTRONS
    conv_thr = 1.0d-10          ! tighter than the H2 default
/
ATOMIC_SPECIES
  Si  28.086  Si.pbe-n-rrkjus_psl.0.1.UPF
ATOMIC_POSITIONS crystal
  Si  0.00  0.00  0.00
  Si  0.25  0.25  0.25         ! diamond structure, fractional
CELL_PARAMETERS angstrom
    0.00   2.715  2.715
    2.715  0.00   2.715
    2.715  2.715  0.00
K_POINTS automatic              ! Monkhorst-Pack mesh
  8 8 8  0 0 0

Step 2 — bands (si.bands.in): same as ‘si.scf.in’ but with ‘calculation = ‘bands’, an added nbnd`, and the k-path:

&CONTROL
    calculation = 'bands'          ! only this line differs from scf
    ...
/
&SYSTEM
    ...
    nbnd      = 12              ! >= number of electrons / 2
    occupations = 'fixed'
/
&ELECTRONS
    conv_thr = 1.0d-10
/
ATOMIC_SPECIES
  Si  28.086  Si.pbe-n-rrkjus_psl.0.1.UPF
ATOMIC_POSITIONS crystal
  Si  0.00  0.00  0.00
  Si  0.25  0.25  0.25
CELL_PARAMETERS angstrom
    0.00   2.715  2.715
    2.715  0.00   2.715
    2.715  2.715  0.00
K_POINTS tpiba_b               ! high-symmetry path; cartesian, 2*pi/bohr
  6
  0.0   0.0   0.0    30         ! Gamma
  0.5   0.0   0.5    30         ! X
  0.5   0.25  0.75   30         ! W
  0.375 0.375 0.375  30         ! K
  0.0   0.0   0.0    30         ! Gamma
  0.5   0.5   0.5    30         ! L

The ‘6’ is the number of segments, the ‘30’ is the number of k-points per segment. Run: ‘pw.x < si.scf.in > si.scf.out && pw.x < si.bands.in > si.bands.out’.

Step 3 — post-processing (si.bands.post.in): extract band energies:

&BANDS
    prefix   = 'si'
    outdir   = './tmp/'
    filband  = 'si.bands.dat'   ! output: k-index, E_1, E_2, ..., E_nbands
    lsym     = .true.           ! enforce the symmetries of the path
/

Run: ‘bands.x < si.bands.post.in > si.bands.post.out’. The file `si.bands.dat’ has one row per k-point and one column per band.

10.1.3 What to tweak for X

10.1.4 Common pitfalls

10.1.5 Cross-references

10.2 VASP

VASP is the Fortran PAW plane-wave code (§2.2). The cheatsheet below covers the four input files — INCAR, KPOINTS, POSCAR, POTCAR. The “POSCAR” is the structure, the “INCAR” the run-time parameters, the “KPOINTS” the k-point sampling, and the “POTCAR” the PAW pseudopotential concatenation (one entry per species, in POSCAR order).

10.2.1 Minimal input — Si SCF

The smallest input that produces a sensible single-point energy on bulk silicon:

# INCAR
ENCUT  = 250       # plane-wave cutoff, eV (typically 1.3 x max ENMAX in POTCAR)
EDIFF  = 1.0e-5    # SCF convergence, eV
ISMEAR = 0         # 0 = Gaussian (semiconductors and insulators)
SIGMA  = 0.05      # smearing width, eV
IBRION = -1        # -1 = no ionic update (static calculation)

# KPOINTS
Automatic mesh
0              # 0 = auto: Monkhorst-Pack with explicit offsets
Gamma          # Gamma-centred (use Monkhorst for non-Gamma-centred)
8  8  8        # mesh
0. 0. 0. # shift

# POSCAR
Si
 5.43           # cubic lattice constant, A
  1.0  0.0  0.0
  0.0  1.0  0.0
  0.0  0.0  1.0
Si
 2              # number of atoms
Cartesian
 0.00  0.00  0.00
 1.3575 1.3575 1.3575   # 1/4 of the cube diagonal

# POTCAR -- built by concatenating per-element POTCAR files from
# the VASP PAW library, in the order the species appear in POSCAR:
$ cat $VASP_PP_PATH/potcar/Si/PBE/Si/POTCAR > POTCAR

To run: ‘mpirun -n 8 vasp_std’. The total energy (in eV) is on the last “free energy” line of OUTCAR.

10.2.2 Annotated example — Si band structure with HSE06

A band structure of bulk silicon using the HSE06 range-separated hybrid (chapter 05 §5.5). Three steps: (i) standard SCF at PBE; (ii) HSE06 SCF on the same mesh (re-using the PBE charge density as the initial guess); (iii) HSE06 non-self-consistent calculation along the k-path. The HSE06 step is $\sim 100\times$ slower than PBE; SCF convergence is harder.

Step 1 — PBE SCF (INCAR.PBE): add LWAVE’ and LCHARG’ to the minimal input above (to write the wavefunction and charge for the next step):

LWAVE  = .TRUE.
LCHARG = .TRUE.

Step 2 — HSE06 SCF (INCAR.HSE): same POTCAR, KPOINTS, POSCAR. INCAR changes:

LHFCALC = .TRUE.        # turn on hybrid functional
HFSCREEN = 0.2          # HSE06 screening, 1/Angstrom
ALGO = Damped           # robust for HSE; 'All' is faster when stable
TIME = 0.4              # damped/All mixing parameter
PRECFOCK = Normal       # 'Fast' (FFT2) or 'Accurate' for the final run
EDIFF  = 1.0e-5         # tighter (1.0e-6) for production band gaps

The HSE06 SCF starts from the PBE charge density in CHGCAR’ if it is present. Run: ‘mpirun -n 8 vasp_std’.

Step 3 — HSE06 bands (INCAR.bands): the k-path is specified in KPOINTS’ as a “line-mode” path:

# KPOINTS (line mode)
K-Path                  # header
6                       # number of high-symmetry segments
0.0  0.0  0.0   ! Gamma
0.5  0.0  0.5   ! X
0.5  0.25 0.75  ! W
0.375 0.375 0.375 ! K
0.0  0.0  0.0   ! Gamma
0.5  0.5  0.5   ! L

The bands INCAR has ICHARG = 11' (read the converged charge from CHGCAR, do not update) and LWAVE = .FALSE.’ (we only need the eigenvalues):

LHFCALC = .TRUE.
HFSCREEN = 0.2
ALGO = Normal           # or 'Eigenvalue'
ICHARG = 11              # read CHGCAR, do not update
LWAVE  = .FALSE.
LCHARG = .FALSE.
EDIFF  = 1.0e-6
NBANDS = 16              # >= number of electrons / 2

To run: ‘mpirun -n 8 vasp_std’. The band energies are in EIGENVAL. Plot them with p4vasp, sumo, or a custom script.

10.2.3 What to tweak for X

10.2.4 Common pitfalls

10.2.5 Cross-references

10.3 GPAW

GPAW is a Python / C PAW code with a real-space grid (§4.2). The cheatsheet below uses the ASE interface (§9.1) — the same script works in serial, in parallel (with ‘mpiexec gpaw python script.py’), and on a GPU (with the `gpaw.cuda’ build).

10.3.1 Minimal script — H₂ SCF

The smallest script that produces a sensible single-point energy on a hydrogen molecule:

# h2.py
from ase import Atoms
from gpaw import GPAW, PW

h2 = Atoms('H2', positions=[(0, 0, 0), (0, 0, 0.74)])
h2.center(vacuum=4.0)             # pad the cell to avoid image interactions

calc = GPAW(
    mode=PW(300),                 # plane-wave cutoff, eV
    xc='PBE',
    kpts=(1, 1, 1),               # Gamma-only for a molecule
    txt='h2.txt',                 # log file
)
h2.calc = calc
energy = h2.get_potential_energy()
print(f'H2 energy: {energy:.4f} eV')

Run: ‘gpaw python h2.py’. The total energy is in eV (the GPAW output is in eV by default, unlike Quantum ESPRESSO which uses Rydberg).

10.3.2 Annotated example — Si band structure in 20 lines

A band structure of bulk silicon with a plane-wave PAW calculation:

# si_bands.py
import numpy as np
from ase.build import bulk
from ase.dft.kpoints import bandpath
from gpaw import GPAW, PW

# 1. Build the structure
si = bulk('Si', 'diamond', a=5.43)

# 2. Standard SCF on a uniform mesh
si.calc = GPAW(
    mode=PW(300),
    xc='PBE',
    kpts=(8, 8, 8),
    convergence={'energy': 1.0e-6},
    txt='si.scf.txt',
)
si.get_potential_energy()                   # runs the SCF
si.calc.write('si.gpw', mode='all')         # save everything for re-use

# 3. Non-self-consistent band structure along a path
calc = GPAW('si.gpw').fixed_density()       # reload, hold density fixed
kpts, x, X = bandpath('GXWKL', si.cell, npoints=60)
calc.set(kpts=kpts)                        # re-attach the calculator
si.calc = calc
si.get_potential_energy()                   # NSCF run

# 4. Extract the band structure
bs = calc.band_structure()
bs.plot(filename='si.bands.png', emax=15)  # one line for the plot

The ‘mode=’all’’ write in step 2 saves bot’h’ the wavefunctions and the density. The fixed_density()' load in step 3 re-uses the SCF density for a non-self-consistent NSCF run along the k-path — exactly the same two-step workflow as QE / VASP, but expressed in Python and re-runnable from a Jupyter notebook. The bs.plot()’ call in step 4 produces a publication-quality band-structure plot with one line of code.

10.3.3 What to tweak for X

10.3.4 Common pitfalls

10.3.5 Cross-references

10.4 SIESTA

SIESTA is a Fortran NAO basis code (§4.3). The cheatsheet below uses the fdf (flexible data format) input — a Fortran-NAMELIST-inspired free format. The pseudopotential is a separate file (.psf' or .vps`).

10.4.1 Minimal input — H₂ SCF

The smallest input that produces a sensible single-point energy on a hydrogen molecule:

# h2.fdf
SystemName          h2
SystemLabel         h2
LatticeConstant     1.0 Ang
%block LatticeVectors
    10.0  0.0   0.0
     0.0 10.0   0.0
     0.0  0.0  10.0
%endblock LatticeVectors
NumberOfAtoms        2
NumberOfSpecies      1
%block ChemicalSpeciesLabel
   1  1  H
%endblock ChemicalSpeciesLabel
%block AtomicCoordinatesAndAtomicNumbers
    0.0  0.0  0.0   1
    0.0  0.0  0.74  1
%endblock AtomicCoordinatesAndAtomicNumbers
PAO.BasisSize        DZP            # SZ | DZ | DZP | TZP
PAO.EnergyShift      0.05 eV        # energy shift of the radial functions
SplitNorm            0.15
MeshCutoff           100 Ry
XC.functional        GGA            # PBE
XC.authors           PBE
SolutionMethod       diagon
MaxSCFIterations     100
DM.MixingWeight      0.1
DM.NumberPulaySteps  5
DM.Tolerance         1.0d-4

The H₂ SCF total energy is in `h2.out’ after ‘siesta: Total energy =’.

10.4.2 Annotated example — Si band structure

A band structure of bulk silicon. The key SIESTA-specific parameters are the PAO basis (PAO.BasisSize, PAO.EnergyShift, SplitNorm) and the MeshCutoff (the real-space grid for the Poisson solver and the local part of the Hamiltonian).

# si.fdf
SystemName          si
SystemLabel         si
LatticeConstant     1.0 Ang
%block LatticeVectors
     0.0   2.715  2.715
     2.715  0.0   2.715
     2.715  2.715  0.0
%endblock LatticeVectors
NumberOfAtoms        2
NumberOfSpecies      1
%block ChemicalSpeciesLabel
   1  14  Si
%endblock ChemicalSpeciesLabel
%block AtomicCoordinatesAndAtomicNumbers
    0.00  0.00  0.00  14
    0.25  0.25  0.25  14
%endblock AtomicCoordinatesAndAtomicNumbers
%block kgrid_Monkhorst_Pack
    8   0   0   0.0
    0   8   0   0.0
    0   0   8   0.0
%endblock kgrid_Monkhorst_Pack
PAO.BasisSize        DZP
PAO.EnergyShift      0.05 eV
SplitNorm            0.15
MeshCutoff           200 Ry
XC.functional        GGA
XC.authors           PBE
SolutionMethod       diagon
MaxSCFIterations     200
DM.MixingWeight      0.1
DM.NumberPulaySteps  5
DM.Tolerance         1.0d-5

For a band structure, set ‘%block BandLines’ instea’d’ of the Monkhorst–Pack kgrid_Monkhorst_Pack:

BandPointsScale      pi/a
WriteBands           .true.
%block BandLines
   1   0.000  0.000  0.000  Gamma
  30   0.500  0.000  0.500  X
  30   0.500  0.250  0.750  W
  30   0.375  0.375  0.375  K
  30   0.000  0.000  0.000  Gamma
  30   0.500  0.500  0.500  L
%endblock BandLines

The first column of BandLines’ is the number of k-points per segment; the last column is a label. SIESTA writes the bands to SystemLabel.bands.

10.4.3 What to tweak for X

10.4.4 Common pitfalls

10.4.5 Cross-references

10.5 CP2K

CP2K is a Fortran mixed Gaussian / plane-wave (GPW) code (§2.5). The cheatsheet below covers the standard ‘cp2k.psmp’ (or ‘cp2k.popt) executable. The input format is a *structured* Fortran-NAMELIST with explicit section keywords (&GLOBAL, &SUBSYS, &DFT, …) and matching &END’ lines. Indentation is optional but is the standard convention.

10.5.1 Minimal input — water in a box

The smallest input that produces a sensible geometry optimisation of a water molecule in a cubic box:

! h2o.inp
&GLOBAL PROJECT h2o RUN_TYPE GEO_OPT PRINT_LEVEL MEDIUM &END GLOBAL

&MOTION
  &GEO_OPT TYPE MINIMIZATION OPTIMIZER BFGS MAX_ITER 100 &END GEO_OPT
&END MOTION

&FORCE_EVAL
  METHOD QUICKSTEP
  &DFT
    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    POTENTIAL_FILE_NAME  GTH_POTENTIALS
    CHARGE 0  MULTIPLICITY 1
    &QS EPS_DEFAULT 1.0E-12 &END QS
    &POISSON PERIODIC NONE PSOLVER MT &END POISSON
    &MGRID CUTOFF 300 REL_CUTOFF 50 &END MGRID
    &SCF MAX_SCF 50 EPS_SCF 1.0E-6 SCF_GUESS ATOMIC &END SCF
    &XC
      &XC_FUNCTIONAL PBE &END XC_FUNCTIONAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL ABC 6.0 6.0 6.0 PERIODIC NONE &END CELL
    &COORD
      O   0.000  0.000  0.000
      H   0.957  0.000  0.000
      H  -0.240  0.927  0.000
    &END COORD
    &KIND O BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q6 &END KIND
    &KIND H BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q1 &END KIND
  &END SUBSYS
&END FORCE_EVAL

To run: ‘mpiexec -n 4 cp2k.psmp -i h2o.inp -o h2o.out’. The geometry optimisation converges in ~5 steps; the total energy is in `h2o.out’ after ‘ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.)’.

10.5.2 Annotated example — NVT molecular dynamics of water

A short Born–Oppenheimer MD of a 32-water box at 300 K, the canonical “CP2K can do this” test case. The &FORCE_EVAL' block from the previous example is re-used almost unchanged; only the &MOTION’ block changes to MD.

! md32h2o.inp
&GLOBAL
  PROJECT  md32h2o
  RUN_TYPE MD
  PRINT_LEVEL LOW
&END GLOBAL

&MOTION
  &MD
    ENSEMBLE NVT
    STEPS 1000
    TIMESTEP 0.5
    TEMPERATURE 300
    &THERMOSTAT REGION MASSIVE
      &NOSE LENGTH 3 YOSHIDA 3 TIMECON 1000 &END NOSE
    &END THERMOSTAT
  &END MD
  &PRINT
    &TRAJECTORY &EACH MD 1 &END EACH &END TRAJECTORY
    &RESTART   &EACH MD 50 &END EACH &END RESTART
  &END PRINT
&END MOTION

! &FORCE_EVAL re-uses the water molecule block from h2o.inp
! above, with the following changes:
&FORCE_EVAL
  METHOD QUICKSTEP
  &DFT
    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    POTENTIAL_FILE_NAME  GTH_POTENTIALS
    &POISSON PERIODIC XYZ PSOLVER PERIODIC &END POISSON
    &MGRID CUTOFF 400 REL_CUTOFF 60 &END MGRID
    &SCF MAX_SCF 20 EPS_SCF 1.0E-6 SCF_GUESS RESTART
      &OT MINIMIZER DIIS PRECONDITIONER FULL_SINGLE_INVERSE &END OT
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE &END XC_FUNCTIONAL
      &VDW_POTENTIAL POTENTIAL_TYPE PAIR_POTENTIAL
        &PAIR_POTENTIAL PARAMETER_FILE_NAME dftd3.dat
          TYPE DFTD3 REFERENCE_FUNCTIONAL PBE
        &END PAIR_POTENTIAL
      &END VDW_POTENTIAL
    &END XC
  &END DFT
  &SUBSYS
    &CELL ABC 9.852 9.852 9.852 PERIODIC XYZ &END CELL
    &TOPOLOGY COORD_FILE_NAME water32.xyz COORD_FILE_FORMAT XYZ
      &CENTER_COORDINATES &END CENTER_COORDINATES
    &END TOPOLOGY
    &KIND O BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q6 &END KIND
    &KIND H BASIS_SET DZVP-MOLOPT-GTH POTENTIAL GTH-PBE-q1 &END KIND
  &END SUBSYS
&END FORCE_EVAL

The water32.xyz' is a standard XYZ with 32 waters (96 atoms). D3 dispersion ([chapter 05](/DFT-notes/dft-notes/chapter-05/) §5.7) is enabled with '&VDW_POTENTIAL ... TYPE DFTD3'. The OT solver (&OT`) is the recommended SCF solver for large systems: it is much more parallel-scalable than the default diagonalisation, at the cost of requiring a localised initial guess.

10.5.3 What to tweak for X

10.5.4 Common pitfalls

10.5.5 Cross-references

10.6 FLEUR

FLEUR is a Fortran FP-LAPW code (§3.3). The cheatsheet below covers the ‘inp’ file (the main input) and the ‘inpgen’ generator. FLEUR uses a two-step workflow: inpgen’ reads a structure (CIF or template) and produces the ‘inp’ file, then ‘fleur’ reads the inp’ and runs.

10.6.1 Minimal input — Fe bcc, spin-polarised

The smallest input that produces a sensible spin-polarised ground state for bcc iron. The structure is taken from Fe.cif; the generator produces the inp’ file.

Step 1 — generate the input:

# in the directory containing Fe.cif:
cat > inp.txt <<EOF
&input
  cartesian=F
  symrel=T
  film=T              ! 3D periodic (set film=F for 2D)
  rmt= 2.0            ! muffin-tin radius, bohr
  lmax= 8
  lnonsph= 4
  kcrel= 0            ! 0 = scalar-relativistic,
                      ! 1 = full-relativistic (with SOC)
/
&lattice
  ext=T               ! read lattice from external file (CIF)
/
&atom
  id= 26.1
  bmu= 2.2            ! initial magnetic moment, mu_B
  lo= 3s 3p           ! local orbitals for semicore states
  econfig= 3d6 4s2
/
EOF
inpgen < inp.txt

The ‘inpgen’ input is read from ‘inp.txt; the output 'inp' is written. 'bmu=2.2' sets the initial magnetic moment; lo=3s 3p’ adds local orbitals for the 3s and 3p semicore states (necessary for an accurate all-electron treatment of 3$d$ transition metals).

Step 2 — run FLEUR:

fleur_MPI

The total energy is in ‘out’ (or ‘out.xml`) after ‘Total energy :’.

10.6.2 Annotated example — Fe bcc with spin–orbit coupling

The same system as above, but with spin–orbit coupling (necessary for a correct magnetic anisotropy energy; chapter 04 §4.9). Re-run ‘inpgen’ with ‘kcrel=1’ (full-relativistic) and add the magnetisation direction ‘cth=30, 30, 30’:

sed -i 's/kcrel=  0/kcrel=  1/' inp.txt   # or edit by hand
cat >> inp.txt <<EOF
&input
  cth= 30.0, 30.0, 30.0   ! magnetisation direction (theta, phi, alpha)
  swsp= 0                 ! 0 = ferromagnet, 1 = spin-spiral
/
EOF
inpgen < inp.txt
fleur_MPI

The MAE is the difference in total energy between two magnetisation directions (e.g. cth=0, 0, 0' — z axis — versus cth=90, 0, 0’ — x axis). With SOC, the k-mesh must be denser (an MAE of 0.1 meV/atom is the typical target), and lmax=8' should be increased to lmax=10’ for production.

10.6.3 What to tweak for X

10.6.4 Common pitfalls

10.6.5 Cross-references

Back to the chapter index or jump to chapter 06 for the theory of basis sets that makes the “Gaussian vs. plane wave vs. real space” split sensible.