Section I: ab initio electronic structure methods Part II: introduction to the computational methodology

1
Section I ab initio electronic structure
methodsPart II introduction to the
computational methodology

• Instructor Marco Buongiorno Nardelli
• 516C Cox - tel. [phone] - e-mail
  mbnardelli_at_ncsu.edu
• Office hours T,Th, 200-300PM
• http//ermes.physics.ncsu.edu
• Suggested readings
• R. M. Martin, Electronic Structure, Cambridge,
  2004, http//www.electronicstructure.org
• W.E. Pickett, Pseudopotential methods in
  condensed matter applications, Comp. Phys. Rep.
  9, 115 (1989)
• J.I. Gersten, F.W. Smith, The Physics and
  Chemistry of Materials, Wiley, 2001
• http//www.pwscf.org

2
Electronic ground state

• Stable structure of solids are classified on the
  basis of their electronic ground state, which
  determines the minimum energy equilibrium
  structure, and thus the characteristics of the
  bonding between the nuclei

• Closed-shell systems rare gases and molecular
  crystals. They remain atom-like and tend to form
  close-packed solids
• Ionic systems compound formed by elements of
  different electronegativity. Charge transfer
  between the elements thus stabilizes structures
  via the strong Coulomb (electrical) interaction
  between ions
• Covalent bonding involves a complete change of
  the electronic states of the atoms with pair of
  electrons forming directional bonds
• Metals itinerant conduction electrons spread
  among the ion cores. Electron gas as electronic
  glue of the system

3
The many-body problem

• How do we solve for the electronic ground state?
  Solve a many-body problem the study of the
  effects of interaction between bodies, and the
  behavior of a many-body system
• The collection of nuclei and electrons in a piece
  of a material is a formidable many-body problem,
  because of the intricate motion of the particles
  in the many-body system
• Electronic structure methods deal with solving
  this formidable problem starting from the
  fundamental equation for a system of electrons
  (ri) and nuclei (RI)

4
Electronic structure methods

• In independent electron approximations, the
  electronic structure problem involves the
  solution of a Schroedinger-like equation for each
  of the electrons in the system
• In this formalism, the ground state energy is
  found populating the lowest eigenstates according
  to the Pauli exclusion principle
• Central equation in electronic structure theory.
  Depending on the level of approximation we find
  this equation all over
• Semi-empirical methods (empirical
  pseudopotentials, tight-binding)
• Density Functional Theory
• Hartree-Fock and beyond
• Mathematically speaking, we need to solve a
  generalized eigenvalue problem using efficient
  numerical algorithms

5
Towards Density Functional Theory

• The fundamental tenet of Density Functional
  Theory is that the complicated many-body
  electronic wavefunction ? can be substituted by a
  much simpler quantity, that is the electronic
  density
• This means that a scalar function of position,
  n(r), determines all the information in the
  many-body wavefunction for the ground state and
  in principle, for all excited states
• n(r) is a simple non-negative function subject to
  the particle conservation sum rule
• where N is the total number of electrons in the
  system

6
Kohn and Sham ansatz

• H-K theory is in principle exact (there are no
  approximations, only two elegant theorems) but
  impractical for any useful purposes
• Kohn-Sham ansatz replace a problem with another,
  that is the original many-body problem with an
  auxiliary independent-particle model
• Ansatz K-S assume that the ground state density
  of the original interacting system is equal to
  that of some chosen non-interacting system that
  is exactly soluble, with all the difficult part
  (exchange and correlation) included in some
  approximate functional of the density.
• Key steps
• Definition of the non-interacting auxiliary
  system
• The auxiliary Hamiltonian contains the usual
  kinetic energy term and a local effective
  potential acting on the electrons
• Actual calculations are performed on this
  auxiliary Hamiltonian
• through the solution of the corresponding
  Schroedinger equation for N independent electrons

7
Kohn and Sham ansatz
Non-interacting auxiliary particles in an
effective potential
Interacting electrons real potential

• The density of this auxiliary system is then
• The kinetic energy is the one for the independent
  particle system
• We define the classic electronic Coulomb energy
  (Hartree energy) as usual

8
Kohn and Sham equations

• Finally, we can rewrite the full H-K functional
  as
• All many body effects of exchange and correlation
  are included in Exc
• So far the theory is still exact, provided we can
  find an exact expression for the exchange and
  correlation term
• The minimization of this functional under the
  particle conservation constraint leads to a set
  of Schroedinger-like equations
• With an explicit effective potential

9
Kohn and Sham equations

• The great advantage of recasting the H-K
  functional in the K-S form is that separating the
  independent particle kinetic energy and the long
  range Hartree terms, the remaining exchange and
  correlation functional can be reasonably
  approximated as a local or nearly local
  functionals of the electron density
• Local Density Approximation (LDA) Excn is a
  sum of contribution from each point in space
  depending only upon the density at each point
  independent on other points
• where is the exchange and correlation
  energy per electron.
• is a universal functional of the
  density, so must be the same as for a homogeneous
  electron gas of given density n
• The theory of the homogeneous electron gas is
  well established and there are exact expression
  (analytical or numerical) for both exchange and
  correlation terms
• Exchange as
• Correlation from exact Monte Carlo calculations
  (Ceperley, Alder, 1980)

10
Kohn and Sham equations

• Finally, the set of K-S equations with LDA for
  exchange and correlation give us a formidable
  theoretical tool to study ground state properties
  of electronic systems
• Set of self-consistent equations that have to be
  solved simultaneously until convergence is
  achieved
• Note K-S eigenvalues and energies are
  interpreted as true electronic wavefunction and
  electronic energies (electronic states in
  molecules or bands in solids)
• Note K-S theory is a ground-state theory and as
  such is supposed to work well for ground state
  properties or small perturbations upon them
• Extremely successful in predicting materials
  properties - golden standard in research and
  industry

11
Overview of solid-state concepts

• In DFT calculations, one solves iteratively the
  Kohn-Sham equations for the electron density in
  the external potential of the ions
• We need to frame the K-S equations in a
  crystalline (or molecular or atomic) environment
  and be able to formulate the problem in a
  computationally tractable manner (a crystalline
  set-up is the most general for our purposes and
  for most available computer codes)
• A crystal is an ordered state of matter in which
  the position of the nuclei are repeated
  periodically in space
• Completely specified by the shape of one repeat
  unit (primitive unit cell) and type and position
  of nuclei in that unit (basis)
• The unit cell is repeated infinitely in space
  through a set of rules that describe the
  repetition (translations) - the Bravais lattice

12
Overview of solid-state concepts

• Simple Bravais lattice (cubic cell) primitive
  cell with one atom basis
• The mimimal primitive cell that can generate the
  whole periodic lattice is called the Wigner-Seitz
  cell, and depends on the particular symmetry of
  the crystal (lattice basis)
• Translations are defined as vectors T multiple of
  the primitive cell vectors
• Tla1ma2na3

Note A primitive cell is defined by 3 unit
vectors directed along the edges of the periodic
box
13
Overview of solid-state concepts

• Any function defined in a crystal has to satisfy
  the periodicity of the system
• This applies to all the quantities that we have
  introduced so far Hamiltonians, wavefunctions,
  electron densities, etc.
• It is well known that periodic functions can be
  represented by Fourier transform in terms of
  Fourier components at specific wavevectors q. We
  define a Fourier transform
• For periodic functions the above expression can
  be written as
• In a periodic crystal, Fourier components are
  restricted to those that are periodic within a
  large volume of crystal made of Ncell
• The set of Fourier components that satisfy the
  above condition form the reciprocal crystal
  lattice. The reciprocal lattice is itself a
  Bravais lattice whose axis bi are defined through
  the relation

14
Overview of solid-state concepts

• Since the reciprocal lattice is a Bravais
  lattice, it has a primitive cell. The minimal
  primitive cell of the reciprocal lattice is
  called Brillouin zone and its principal axis
  are given by
• The reciprocal space is comprised of a lattice of
  points defined by the reciprocal lattice vectors
  Glb1mb2nb3so that the Fourier transform of the
  periodic function can be written as

15
Overview of solid-state concepts

• Electrons in solids are subject to the external
  potential of the nuclei that has the periodicity
  of the crystal
• Independent electrons (such as the ones of DFT in
  the K-S approach) in a periodic potentials have a
  very important property as a general consequence
  of the periodicity of the potential
• Blochs theorem the eigenstates of the crystal
  Hamiltonian (the electrons wavefunctions) can be
  always chosen as the product of a plane wave
  times a function with the periodicity of the
  Bravais lattice
• with k within the primitive cell of the
  reciprocal lattice
• Note for the crystal we have r and T, for the
  reciprocal space we have k and G
• All possible eigenstates of the Hamiltonian are
  specified by k within the primitive cell of the
  reciprocal lattice. For each k there is a
  descrete set of eigenvalues that form what are
  the energy bands of the crystal. Usually they are
  calculated within the Brillouin zone and along
  specific directions in reciprocal space
• For a complete discussion on reciprocal space
  vectors and periodic boundary conditions see the
  Ashcroft and Mermin book or any other book on
  Solid-state theory

16
Example

• Brillouin zone and band structure for a cubic fcc
  solid Au

From http//cst-www.nrl.navy.mil/ElectronicStruct
ureDatabase/
17
An electronic structure code PWscf

• PWscf is a state-of-the-art software package that
  we will use as an introduction to the actual
  procedure of modeling the electronic structure of
  a solid
• http//www.pwscf.org

18
An electronic structure code PWscf
19
An electronic structure code PWscf

• Best first step READ THE USER MANUAL!

20
An electronic structure code PWscf

• Three basic preliminary steps. From the Users
  manual
• Installation The PWscf package can be downloaded
  from the http//www.pwscf.org site. Presently,
  only source files are provided. Some precompiled
  executables (binary files) are provided only for
  the GUI. Providing binaries for PWscf would
  require too much effort and would work only for a
  small number of machines anyway. Uncompress and
  unpack the code in an empty directory of your
  choice that will become the root directory of the
  distribution. On Linux machines, you may use
• tar -xvzf pw.2.0.tgz.
• On other Unix machines gzip -dc pw.2.0.tgz
  tar -xvf -
• Automatic configuration An experimental
  automatic configuration, using the GNU
  "configure" utility, is available (thanks to
  Gerardo Ballabio, CINECA). From the root
  directory, type
• ./configure
• The script will examine your hardware and
  software, generate dependencies needed by the
  Makefile's, produce suitable configuration files
  make.sys and make.rules. Presently it is expected
  to work for Linux PCs, IBM sp machines, SGI
  Origin, some HP-Compaq Alpha machines. For more
  details, read the INSTALL file.

21
An electronic structure code PWscf

• Compilation There are a few adjustable parameters
  in Modules/parameters.f90. The present values
  will work for most cases. All other variables are
  dynamically allocated you do not need to
  recompile your code for a different system. You
  can compile the following codes
• make pw produces PW/pw.x and PW/memory.x. pw.x
  calculates electronic structure, structural
  optimization, molecular dynamics, barriers with
  NEB. memory.x is an auxiliary program that checks
  the input of pw.x for correctness and yields a
  rough (under-)estimate of the required memory.
• make ph produces PH/ph.x. ph.x calculates phonon
  frequencies and displacement patterns, dielectric
  tensors, effective charges (uses data produced by
  pw.x).
• make d3 produces D3/d3.x d3.x calculates
  anharmonic phonon lifetimes (third-order
  derivatives of the energy), using data produced
  by pw.x and ph.x.
• make gamma produces Gamma/phcg.x. phcg.x is a
  version of ph.x that calculates phonons at 0
  using conjugate-gradient minimization of the
  density functional expanded to second-order. Only
  the ( 0 ) point is used for Brillouin zone
  integration. It is faster and takes less memory
  than ph.x, but does not (yet) support Ultrasoft
  pseudopotentials.
• make pp produces a variety of post-processing
  codes
• make tools produces utility programs, mostly for
  phonon calculations

22
An electronic structure code PWscf

• Running Pwscf for an electronic and ionic
  structure calculation
• Main information the input file

23
PWscf input file

• The input data is organized as several namelists,
  followed by other fields introduced by keywords.
  The namelists are
• CONTROL general variables controlling the run
• SYSTEM structural information on the system
  under investigation
• ELECTRONS electronic variables
  self-consistency, smearing
• IONS (optional) ionic variables relaxation,
  dynamics
• CELL (optional) variable-cell dynamics
• PHONON (optional) information needed to produce
  data for phonon calculations
• Optional namelist may be omitted if the
  calculation to be performed does not require
  them. This depends on the value of variable
  calculation in namelist CONTROL. Most variables
  in namelists have default values. Only the
  following variables in SYSTEM MUST be specified
  
• ibrav (integer) bravais-lattice index
• celldm (real, dimension 6) crystallographic
  constants
• nat (integer) number of atoms in the unit cell
• ntyp (integer) number of types of atoms in the
  unit cell
• ecutwfc (real) kinetic energy cutoff (Ry) for
  wavefunctions.

Description of all the input cards can be found
in the file INPUT_PW
24
PWscf input file

• Input file for an electronic structure
  calculation for a bulk Si crystal
• control
• scf self-consistent solution of the Kohn-Sham
  equations (ground state energy, electron
  densities) from an arbitrary initial density
  (from_scratch)
• prefix prefix for file names
• tstress,tprnfor compute also forces and
  stresses in the given geometry
• Specify working directories (optional)

25
Input parameters system

• The system namelist is the one where we specify
  the geometrical parameters of our simulation
  cell
• ibrav specifies the Bravais lattice type
• celldm specifies the crystallographic constants
  (the dimension of the simulation cell)

26
Input parameters system

• nat and ntyp specify the number of atoms in the
  simulation cell and how many different atomic
  species are in the system we want to study
• Together with the information contained under the
  keywords ATOMIC SPECIES and ATOMIC POSITIONS they
  determine the basis in the primitive cell
• ATOMIC SPECIES list of all the atomic species
  in the system with information on the atomic mass
  and a link to an external file that contains the
  information necessary to describe that particular
  atom - pseudopotential file
• ATOMIC POSITIONS list of all the different
  atoms in the simulation cell with the
  specification of their atomic coordinates

27
Input parameters the pseudopotential
28
Numerical solution plane waves

• Kohn-Sham equations are differential equations
  that have to be solved numerically
• To be tractable in a computer, the problem needs
  to be discretized via the introduction of a
  suitable representation of all the quantities
  involved
• Various discretization approeches. Most common
  are Plane Waves (PW) and real space grids.
• In periodic solids, plane waves of the form
  are most appropriate since they reflect the
  periodicity of the crystal and periodic functions
  can be expanded in the complete set of Fourier
  components through orthonormal PWs
• In Fourier space, the K-S equations become
• We need to compute the matrix elements of the
  effective Hamiltonian between plane waves

29
Numerical solution plane waves

• Kinetic energy becomes simply a sum over q
• The effective potential is periodic and can be
  expressed as a sum of Fourier components in terms
  of reciprocal lattice vectors
• Thus, the matrix elements of the potential are
  non-zero only if q and q differ by a reciprocal
  lattice vector, or alternatively, q kGm and q
  kGm
• The Kohn-Sham equations can be then written as
  matrix equations
• where
• We have effectively transformed a differential
  problem into one that we can solve using linear
  algebra algorithms!

30
Input parameters ecutwfc

• In this representation both the potentials and
  the Bloch functions, solution of the K-S problem,
  are expanded on a set of plane waves
• In principle, the plane waves basis set is
  infinite, since I have an infinite number of
  reciprocal lattice vectors (or, in other words,
  Fourier components).
• In practice, we need to limit ourselves to a
  finite basis set for the practical solution of
  the linear equations approximation!
• Remember in Fourier space, that is in reciprocal
  space in a crystal, small q components describe
  long-range features (wave-length), while large q
  components, describe short-range features. Very
  sharp oscillations, for instance, need to be
  described by a large number of plane waves with
  large G vectors
• Increasing the dimension of the basis (number of
  plane waves, so larger magnitudes of Gm) allows
  for a better description of short-range features
  in either the potential or the density.
• ecutwfc is the parameter that controls the number
  of PWs in the basis, so it affects directly the
  accuracy of the calculation convergence
  parameter
• Gmax2 ? ecutwfc

31
Input parameters electrons

• Kohn-Sham equations are always self-consistent
  equations the effective K-S potential depends on
  the electron density that is the solution of the
  K-S equations
• In reciprocal space the procedure becomes
• Iterative solution of self-consistent equations -
  often is a slow process if particular tricks are
  not used mixing schemes

where
and
32
Mixing schemes

• Key problem updating the potential and/or the
  density between successive iteration steps (loop
  in the solution of the K-S equations)
• A direct approach, where we start from an
  arbitrary density and use the solution of the K-S
  loop as input of the next will not work -
  instabilities in the solution of the minimization
  problem - the numerical procedure does not
  converges into the minimum

• Simplest approach to resolve the issue is linear
  mixing estimate an improved density input nini1
  at step I1 as the linear combinantion of input
  and output densities nini and nouti at step i

• mixing_mode particular algorithm to mix input
  and output density/potentials
• mixing_beta numerical value of ?
• conv_thr convergence threshold for
  self-consistency, estimated energy error lt
  conv_thr

33
Mixing schemes

• ? is the parameter that controls the rate of
  convergence
• Large ? output density from previous iteration
  weighs most - fast convergence (large density
  updates from one step to the other). Works well
  for strongly bound, rigid systems (insulators or
  semiconductors, regions around the ionic cores)
• Small ? output density weighs lest - slower
  convergence (smaller density updates from one
  step to the other) but more stable. Works well
  for soft systems such as metals, alloys,
  surfaces or open systems (lots of empty space).
• Various mixing scheme beyond linear mixing have
  been developed and are available in PWscf
  (Anderson, Broyden, etc.)

34
Symmetry and special points

• Crystals are very symmetric systems whose
  geometrical properties are best described by the
  specification of their space group
• Space group of a crystal is the group composed by
  the whole set of translations and point symmetry
  operations (rotations, inversions, reflections
  and combinations of the above) that leave the
  system invariant (including a particular set of
  operations that combines a rotation with a
  non-integer translation, or glide, of a fraction
  of a crystal translation vector, called
  non-symmorphic operations)
• Any function that has the full symmetry of the
  crystal is invariant upon any operation of the
  space group Sng(r) g(Snr)
• In particular, since the Hamiltonian is invariant
  upon any symmetry operation, any operation Sn
  leads to a new equation with r ? Snr and k ? Snk
• The new solution of the transformed equation is
  still an eigenfunction of the Hamiltonian with
  the same eigenvalue
• A high-symmetry k point is defined by the
  identity relation Snk ? k. Helpful in the
  classification of electronic states.
• One can define the Irreducible Brillouin Zone
  (IBZ), which is the smallest fraction of the
  Brillouin Zone that is sufficient to determine
  all the information on the electronic structure
  of the crystal. All the properties for k outside
  the IBZ are obtainable via symmetry operations

35
Symmetry and special points

• This concept becomes particularly significant
  when we have to compute properties that require
  an integration over the Brillouin Zone, such as
  energy or electron density.
• In general an integral over the BZ becomes a sum
  over a discrete set of states corresponding to
  different k points
• Two important issues
• To have an accurate numerical integration the
  discrete set of k values has to be dense enough
  to have a sufficient number of points in regions
  where the integrand varies rapidly. Crucial
  difference between metals and insulators. Since
  insulators have only filled bands, integrals can
  be computed using a few well-chosen points in the
  BZ special points
• Symmetry must be used to reduce the calculations
  to ones that involve only k-points comprised in
  the IBZ

36
Input parameters K-POINTS

• Set of special points in the BZ can be chosen for
  an efficient integration of smooth periodic
  functions
• Most general method has been proposed by
  Monkhorst and Pack (implemented in PWscf).
  Uniform sets of points are chosen in reciprocal
  space according to the formula
• They form a uniform mesh in reciprocal space
• Can be specified in an automatic way by the
  program
• 
  nk1,2,3 N1,2,3 and k1,2,3 0 or 1
• 
  0 centered grid, 1 shifted grid

• where bi are the reciprocal space basis vectors
  and Ni are the parameters that determine the mesh
  size (ni 1,2,,Ni)

37
Input parameters K-POINTS

• A sum over a uniform set of special k-points
  integrates exactly a periodic function that has
  Fourier components that extend only to NiTi in
  each direction
• This logic can be easily understood in one
  dimension
• the value of the following integral
• is given by the value of the integrand function
  f2(k) sin(k) at the mid-point
• k ?, f2(k ?) sin(?)0
• If one has two integrand functions,
  f2(k)A1sin(k) A2sin(2k), the integral is given
  by the sum over two points
• For functions with the symmetry of the full
  crystal symmetry group, the density of the
  special k-point mesh gives a direct measure of
  the accuracy with which we compute the integrals
  convergence parameter

38
Input parameters K-POINTS

• Integrals over the full Brillouin Zone (BZ) can
  be replaced by integrals over the first Brillouin
  Zone (IBZ). In this way we need to perform
  summations on a subset of the full k-point mesh
  in the BZ.
• NOTE All the k points of the full BZ can be
  obtained from the k points in the IBZ via the
  symmetry operations of the crystal space group.
• If we define the weight factor wk as the total
  number of distinguishable k points related by
  symmetry to a k point in the IBZ divided by the
  total number of points Nk.
• With this proviso, any sum over the BZ can be
  written as a sum over the IBZ that includes the
  appropriate weight for each point

39
A few notes on convergence

• The accuracy of an electronic structure
  calculation depends on various approximation
  factors, both physical and numerical
• Physical approximations
• description of nuclei/ions all-electron vs.
  pseudopotential, type of pseudopotential used
• choice of exchange and correlation functional
  (LDA vs. GGA. vs. hybrid functionals or mixed
  schemes (this is usually done in the
  pseudopotential input file)
• Numerical approximations
• The accuracy of the basis set
• In a plane wave basis, the extension of the
  reciprocal space mesh (Fourier basis)
• In a real space calculation, the density of the
  grid in real space
• The accuracy of the integrals in reciprocal
  space density of the special k-point mesh.
  Particularly important for metallic systems,
  where large sets of k-points have to be used for
  a consistent description of the system
• The accuracy of the self-consistent solution of
  the K-S equations mixing schemes and convergence
  threshold

40
Results

• Calculation for a crystal of Si
• Diamond structure
• Equilibrium lattice parameter
• FCC cell with two atoms in the basis - 48 point
  symmetry operation with non-symmorphic
  translations

a3
a2
a1
41
Results

• pw.x lt si.scf.in gt out
• Preliminaries
• First part geometry
• Warning on the existence of non-symmorphic point
  group operations

42
Results

• More on geometry and computational parameters
• lattice parameter and volume of the cell, atoms
  and species, kinetic energy cut-off (ecutwfc),
  convergence threshold and mixing parameters,
  exchange and correlation functional, real and
  reciprocal space basis vectors

43
Results

• Pseudopotential parameters
• Symmetry operations and atomic coordinates
• k-mesh for BZ integration

44
Results

• Iterate!

45
Results

• Convergence!
• Electronic energies (eigenvalues of the K-S
  equation) for the occupied bands of the solid
• Total energy ground state energy minimum of
  the K-S functional
• Individual contributions to the ground state
  energy (given in an alternative way, that groups
  together various contributions and uses the sum
  of the eigenvalues of the K-S equations)
• band energy
• 1-electron contribution
• Hartree energy
• x-c energy
• Ewald (ionic) energy
• Implicitly, we clearly know the ground state
  electron density

46
Results

• From the knowledge of the electron density, we
  can compute all the ground state properties of
  interest
• Forces (in this case are zero because of
  symmetry)
• Stress (zero, or negligible, since the
  calculation has been done at the equilibrium
  lattice parameter)

47
Results

• Collect run-time statistics and finish the run

48
Project

• Using the input file given on your home directory
  on pharos.physics.ncsu.edu do the following
  exercises
• Study the variation of the total ground state
  energy as a function of the convergence
  parameters
• ecutwfc 2,5,8,12,18,24,32 Ry with the given
  k-point mesh. Discuss the results and find the
  converged value for Etot.
• k-point mesh using the converged value for
  ecutwfc from above, use the option automatic
  and check convergence for various grid sizes both
  centered and not
• 1 1 1 0 0 0, 1 1 1 1 1 1, 2 2 2 0 0 0, 2 2 2 1 1
  1, etc.
• Discuss the results
• Displace one of the atoms in the cell by a small
  amount along the diagonal of the cell and repeat
  the study of point 1 above looking at the values
  of the forces and stresses. Discuss the results
• Using the ideal geometry, do a study of total
  ground state energy versus lattice parameter
  (volume). The minimum of the curve will give you
  the theoretical lattice parameter for Si
  corresponding to the pseudopotential used
  (GGA-PBE Si.pbe-rrkj.UPF)
• Repeat the study of point 3 above using the LDA
  pseudopotential (Si.pz-vbc.UPF). Discuss the
  results.