DFT Calculations on Periodic Systems Herbert Früchtl

DFT Calculations on Periodic
Systems
Herbert Früchtl
Outline
• Theory
–
–
–
–
DFT on periodic systems
Basis Sets
Pseudopotentials
Functionals
• CASTEP, VASP and SIESTA input syntax and
capabilities
• Steps in a typical surface science project
• Visualisation of results
• How to do a calculation on wardlaw
• Links to more help
Density Functional Theory
Hohenberg-Kohn Theorem:
All ground state properties can be determined as a
functional of the electron density
Kohn-Sham equations
Can be solved self-consistently, using a basis set to
express ψi(r), analogous to Hartree-Fock-Roothaan
method
Periodic Systems
Bloch Theorem



ik  r
 nk (r )  e unk (r )



unk (r )  unk (r  R )
Eigenfunctions of a periodic
potential can be expressed as
plane wave x periodic function
• After Fourier transformation,
only need to solve for 1st
Brillouin zone
• Plane waves are natural (but
not the only) choice of basis
set
Basis Sets
Plane Waves
 Good at describing periodic systems
 Straight-forward way of extending
basis
 No geometry dependence (BSSE)
 Mathematically simple
 Not good at describing “rough”
features
 Does not take advantage of geometric
knowledge (e.g. empty space)
 Need many functions to converge
result
Atom Centred
 Good at describing core
 Correspond to atomic orbitals
 Relatively few functions
necessary
 Better scaling for large systems
due to integral screening
 Not orthogonal
 BSSE and linear dependencies
Different types available:
Gaussian, Slater, numerical
Mixed
Atom centred functions in non-overlapping spheres
Plane waves in interstitial regions
Many flavours: PAW, FLAPW, GPW, …
Pseudopotentials
Problems with core electrons
• wave function has many nodes
need large basis set to describe
• high kinetic energy
 relativistic effects
But: they are barely changed by
bonding
Replace core with potential
that “looks like a core” to
valence electrons outside
the core area
Pseudopotential Gobbledigook
Small/large core
How many electrons are removed?
Hardness
A measure for the cutoff-radius (smaller=harder)
US-PP
Ultrasoft potentials: very fast
Transferability
Does it work for different chemical environments?
Norm-conserving
US-PP’s are not norm-conserving
PAW: Projector Augmented Wave
Construct wavefunction in core area from
pseudopotential wavefunction
All-electron wavefunction (frozen core) based on US-PP calculation
• Little difference for structural properties, but better electronic
description (charges, DOS, …)
• Recommended method in VASP
• GIPAW (Gauge Invariant PAW) used for NMR calculations in CASTEP
Density Functionals
Less choice than for molecules
• LDA
• GGA
– PBE
– PW91
Very similar in cost
and accuracy
• Hybrid functionals
– now available in CASTEP and VASP
– very expensive
• VdW correction (Grimme’s DFT-D2) available
• A few meta-GGA functionals appearing…
K-Points
Band Structure of fcc Ag
Need to sample Brilloiun zone
– Monkhorst-Pack grid: equallyspaced points in BZ
– For band structure, specify
explicit k-points along
high-symmetry axes in BZ
BZ describes interaction with remote
cells in real space.
Need more k-points for
– small unit cells
– metals
– magnetic materials
Programs
(marked in blue if available in St Andrews )
CASTEP
CASTEP
Cambridge Serial Total Energy Package
CPMD
Plane wave pseudopotential code
Abinit
Fairly high accuracy
Dacapo
CRYSTAL
VASP
Vienna Ab-initio Simulation Package
VASP
Similar to CASTEP in theory
PWscf (QUANTUMESPRESSO)
Faster, but not as robust and user friendly
SIESTA
No NMR properties
QuickStep (CP2K)
NWChem
SIESTA
Conquest
Spanish Initiative for Electronic Simulations
ASW
with Thousands of Atoms
Numerical atomic orbital basis + pseudopotentials
Wien2k
Good scaling towards large numbers of atoms
More at http://www.psi-k.org/codes.shtml
CASTEP input and output
Need (at least) two input files:
<project>.cell
geometry, k-points
<project>.param
everything else
Output will be appended to
<project>.castep
Pseudopotential files
created automatically by default (Vanderbilt US-PP),
unless you want others
CASTEP .cell File: Geometry
Free-form and case
insensitive
! starts comment
Multiline inputs in %block /
%endblock
Cell dimensions can be
given as Cartesian
vectors or as lengths and
angles
Atomic positions in
fractional or absolute
(Cartesian) coordinates
CASTEP .cell File: K-points and Constraints
K-points
kpoints_mp_grid <nx> <ny> <nz> Number of points in xyz directions
kpoints_mp_spacing <R>
Max distance between k-points in Å-1
%block KPOINTS_LIST
<x1> <y1> <z1> <weight1> Explicit list (Σ(weights)=1)
<x2> <y2> <z2> <weight2>
…
Only one of these!
%endblock KPOINTS_LIST
Simple geometric constraints
fix_all_cell = true
Default is not to freeze
fix_all_ions = true
anything
More complex linear constraints possible using blocks
CELL_CONSTRAINTS and IONIC_CONSTRAINTS
CASTEP .param File
Format:
<parameter> : <value> [units]
Important parameters:
task : GeometryOptimization
or SinglePoint, MolecularDynamics
xc_functional : PBE
or LDA, PW91, …
cut_off_energy : 400 eV
accuracy of basis set
max_SCF_cycles : 100
geom_max_iter : 100
fix_occupancy : false
Necessary for metals
smearing_width : 0.1 eV
Only compare total energies
obtained with same smearing!
Order of input is irrelevant
VASP Files
Need (at least) four input files:
POSCAR geometry
POTCAR pseudopotentials and method
KPOINTS grid of points in Brillouin zone
INCAR
everything else
Files MUST have these names
Output:
OSZICAR short overview
OUTCAR human readable output
Several other files for analysis or restart
POSCAR file
Title (ignored)
Lattice constant
lattice vectors and cartesian
coordinates will be multiplied by this
Lattice vectors
Number of atoms of each species
(here: 7 atoms of 3 species)
“selective dynamics”
Optional. Necessary if atoms to be
frozen
Direct (fractional) / cartesian
type of atom positions
Atomic coordinates and constraints:
F=frozen (actually FALSE…), T=free
Optionally: velocities
(in case of dynamics)
Information on atom types comes
purely from POTCAR file!
POTCAR File
Pseudopotential (including electron
configuration and functional to use) for
each atom species in same order as in
POSCAR
• Created by concatenating atomic PP files:
uncompress -c /usr/local/VASP/potpaw_PBE/Rh/POTCAR.Z >POTCAR
uncompress -c /usr/local/VASP/potpaw_PBE/C/POTCAR.Z >>POTCAR
uncompress -c /usr/local/VASP/potpaw_PBE/O/POTCAR.Z >>POTCAR
Overwrite if it exists
Append
PBE GGA functional and corresponding
pseudopotential with PAW
KPOINTS (Regular Grid)
Automatic mesh
0
Gamma
4 4 4
0. 0. 0.
Title (comment)
Number of points (0 = automatic)
Centering (Gamma/Monkhorst)
Points along reciprocal vectors
Shift from centre (optional)
VASP will take care of symmetry
KPOINTS for Band Structure
Band structure
10 ! 10 points per line
Line-mode
cartesian
0
0
0 ! gamma
0
1
0 ! X
0
1
0.5 1
0
0
! X
! W
0.5 1
0
1
0
0
! W
! gamma
Best to run a calculation with full K-grid first, and then
calculate babnd structure at unchanged density
(ICHARG=11 in INCAR)
Example INCAR File
general:
SYSTEM = Cu 111 surface
ISTART = 1
ICHARG = 2
ENMAX = 270
ISMEAR = 2
SIGMA = 0.2
ALGO = Very_fast
geometry optimisation:
NSW = 100
POTIM = 0.2
IBRION = 2
Title
Start / restart
initial guess charges
energy cutoff in eV
Initial Fermi smearing (for metals)
Smearing width
fast or robust?
max. number of steps
scaling for forces
optimisation/dynamics algorithm
Many more options. See manual!
SIESTA Input
Files
FDF (Flexible Data Format) file
Pseudopotential file for each element (download from
http://www.uam.es/departamentos/ciencias/fismateriac/siesta/
or create your own (ask me…))
FDF format
<parameter> <value> [units]
Case insensitive; random order
Characters -_. ignored
Lattice-Constant = lattice_constant = LatticeConstant
# starts a comment
Logical values: T = true = .true. = yes
Complex data structures in blocks (like CASTEP)
%block <label>
…
%endblock <label>
SIESTA: A Simple Input
Name of project
Root for file names
Atoms and elements
Scaling factor for further geometry
information
Unit cell (also possible to specify through
LatticeVectors)
Atomic coordinates (various methods,
including Z-matrix)
Monkhorst-Pack grid (one of several
ways of specifying it)
Functional (LDA, GGA, hybrid)
Geometry optimisation
Properties to Calculate
• Total energy
• Equilibrium Structures
with and without optimising the unit cell
• Dynamics
• Band Structures
• Density of States
Local DOS (projected on axis or atom)
• Mulliken Charges (even with CASTEP…)
• Phonon Spectrum
• NMR: chemical shift, EFG (only CASTEP)
How To Simulate Adsorption
on surfaces or microporous structures
Solid
Molecule
• Periodic
• Fractional
coordinates
• Plane Waves
• Band structure
• Isolated (vacuum or
solvent)
• Z-Matrix
• LCAO
• Energy levels
?
Cluster Calculations
Treat limited part of solid as a
molecule
• Use any molecular program
• In case of molecular crystal, cap
dangling bonds with hydrogen
• Freeze all or most of the solid at
bulk geometry
• Can not determine bulk
geometry
• For metals, consider using
ECPs (=pseudopotentials)
Periodic Surface Model
Use sufficiently large unit cell
– Vacuum slab between bulk
layers
– Potential self-interaction
between adsorbate
molecules
– Need to freeze lower
layer(s) of bulk
Example: Adsorption Energy
Steps
• Optimise bulk structure
–
–
–
–
Create unit cell geometry
Converge energy cutoff to a few meV/atom
Converge k-points to a few meV/atom
Optimise cell size and ionic positions
A
• Relax surface
– Create surface cell with vacuum slab large enough to accommodate
adsorbate and enough layers to get bulk properties
– Re-converge k-points (need less now)
B
– Freeze cell dimensions and lowest layer(s)
– Optimise geometry (1)
• Optimise complete system
– Add adsorbate molecule
– Optimise geometry (3)
C
• Optimise adsorbate
– Remove adsorbent
– Optimise adsorbate in the same cell with the same parameters (2)
• Calculate adsorption energy
– Difference in final energies: ΔE = E(3) – (E(1) + E(2))
A
Create Unit Cell
• Convert existing structure file
– Scripts on wardlaw in
/usr/local/progs/CASTEP-6.0/scripts
– cryst (syntax like babel) translates
between fdf, cell, poscar, …
– Most viewers allow saving in different
formats
• Build from scratch
– Fairly easy for basic crystals in fractional
coordinates
– Atoms at top and rear are not part of the
unit cell!
• Use primitive cell if possible
– This is not what we did in the fcc
examples…
B
Create Surface Slab
• Use a calculator and a bit of
trigonometry…
– Make sure you use the lattice constant from
bulk optimisation
• Use tetr (non-graphical)
– Tool to create or modify inputs for VASP,
SIESTA and (less so) CASTEP
– Can be used to freeze bottom layers
• Use gdis (graphical)
– No support for VASP, but easy to convert
from cif or fdf formats
– Less powerful than tetr, but easier to use
Example: (1 1 5) fcc-Surface with gdis
C
Add Molecule to Surface
gdis viewer can add, delete
or move atoms
Select “Editing” in Model menu
Add atom:
–Select “Add atoms”
–Specify element
–Left-click in model window
Move atom
–Select “Normal Mode”
–Select atom
–Ctrl-middle-mouse-button to
move atom
Or: change coordinates
C
More Complex Molecules
Combine Cartesian coordinates
Use a spreadsheet to move molecule
Use Z-Matrix input for SIESTA
• Define molecule by distances, angles and
dihedrals
• Can combine Cartesian and Z-Matrix coordinates
Gdis can copy between different files
• Open slab model (e.g. fdf file)
• Open molecule file (e.g. xyz or Gaussian log)
• Select molecule, Edit -> Copy
• Switch to slab window, Edit -> Paste
• Move molecule with ctrl-middle-mousebutton
Visualising Periodic Structures
Want to see
• Geometry
• Electron density
• Orbitals
• Spectra
• DOS
Available programs
• Gdis
Blue: available on wardlaw
• XCrysDen
• Mercury
• VESTA
• Discovery Studio Visualizer
• Materials Studio
Gdis
• Fairly good rendering
• Can read .castep and
write .cell or .fdf files
(problems reading .cell)
• Display animation of MD
or optimisation
• Calculate diffraction
spectra (X-Ray, neutron,
electron)
• The only reasonable
builder (except Materials
Studio)
Mercury
• Part of Cambridge
Structural Database
• Very good rendering
• Use castep2cif.py
script to create cif file
• Can calculate X-ray
pattern
XCrysDen
• Good rendering
• Very good at
displaying molecular
surfaces (density, …)
• Scripts to convert
different formats to .xsf
in
/usr/local/progs/CASTEP6.0/scripts
/usr/local/progs/siesta3.0/Utils
How to do a Calculation
• Get an account on wardlaw
– Email me!
[email protected]
– Look at introduction lecture or ask me
• Create input files
– Templates available.
• Create job script
– Templates available…
• Submit batch jobs
– See introduction lecture (sbatch, squeue, scancel…)
• Sample job scripts available in
/usr/local/examples
Links
This lecture and previous ones:
Available soon!
CASTEP
http://www.tcm.phy.cam.ac.uk/castep/
SIESTA:
http://http://icmab.cat/leem/siesta//
CPMD (with a very good tutorial)
http://www.cpmd.org/
VASP (also with good tutorial)
http://www.vasp.at/
Other Lectures and Seminars
Lectures on
 Molecular Quantum Chemistry (Monday)
 DFT on Periodic Systems (Tuesday)
Computation in Chemistry Seminars
See Chemistry Newsletter
Sign up to mailing list with Tanja van Mourik
Postgrad course on computational chemistry
Email your name and supervisor to Tanja van Mourik
[email protected]
Not only for postgrads!