Electronic optimization

Georg KRESSE

Institut für Materialphysik and Center for Computational Material Science

Universität Wien, Sensengasse 8, A-1090 Wien, Austria

b-initio

ackage
ienna imulation

G. K RESSE , E LECTRONIC O PTIMISATION Page 1

 Overview

Determination of the electronic grounstate

– general strategies
– strategy adopted in VASP
iterative matrix diagonalization and mixing
– how to overcome slow convergence

molecular dynamics
the algorithms are particularly well suited for molecular dynamics

G. K RESSE , E LECTRONIC O PTIMISATION Page 2

 Density functional theory according to Kohn-Sham

density and kinetic energy:

sum of one electron charge densities and kinetic energies

Ne 2
ρtot r ∑ ψn r ρion r
 
2
2 Ne number of electrons









n 1


h̄2 Ne 2 1 ρtot r ρtot r 3 3

2me n∑1



ψn r ∇2 ψn r d 3 r Exc ρ r


2 d rd r





 

 2 r r

 



 
 
 





 LDA/GGA





kinetic energy electrost. energy

KS-functional has a (the) minimum at the electronic groundstate

G. K RESSE , E LECTRONIC O PTIMISATION Page 3

 Numerical determination of the Kohn-Sham groundstate

direct minimization of the DFT functional (Car-Parrinello, modern)

start with a set of wavefunctions ψn r n 1 Ne 2 (random numbers) and minimizes the

 







value of the functional (iteration)

h̄2 2
Gradient: Fn r ∇ V eff r ψn r εn ψn r





2me







iteration – self consistency (old fashioned)
start with a trial density ρ, set up the Schrödinger equation, and solve the Schrödinger equation
to obtain wavefunctions ψn r



h̄2 2
∇ V eff r ρ r ψn r ε n ψn r n 1 Ne 2



 
2me














!
as a result one obtains a new charge density ρ r ∑ n ψn r 2 and a new Schrödinger
 



equation iteration
"

G. K RESSE , E LECTRONIC O PTIMISATION Page 4

 disordered diamond, insulator disordered fcc Fe, metal
energy 2
0
0 n=4
-2
-2 direct

log 10 E-E 0
DCC DCC BA DCC BA @?? BA @?? >= @?? >= <;; >= <;; <;; :9 :9 877 :9 877 6877 6 6 43 43 21 43 21 0/ 21 0/ 0/ .- .- ,++ .- ,++ *) ,++ *) ('' *) ('' &%('' &% $## &% $## $##
D D BA D BA @ BA @ >= @ >= < >= < < :9 :9 8 :9 8 568 56 56 43 43 21 43 21 0/ 21 0/ 0/ .- .- , .- , *) , *) ( *) ( &%( &% $ &% $ $

log 10 E-E 0
-4 direct Car−Parrinello
Car−Parrinello -4
-6 -6 n=1
n=1,2,4,8 n=1 n=8
555

L -8 -8
n=8
1 self.consistent 1 self.consistent
0 0
|

log 10 |F-F 0 |
0

-1 -1
log 10 |F-F

-2 -2
n=1
-3 -3
-4 -4
0 5 10 15 20 0 10 20 30 40
iteration forces iteration

G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).

G. K RESSE , E LECTRONIC O PTIMISATION Page 5

 Direct minimization (not supported by vasp.4.5)

preconditioned conjugate gradient algorithm was applied

h̄2 2
Gradient: Fn r ∇ V eff r ψn r εn ψn r





2me







the main troubles are
– to keep the set of wavefunctions ψn r n 1 Ne 2 orthogonal


  





!

– “sub-space” rotation
at the end one aims to have a set of wavefunction ψn r n 1 Nbands that


E

  






diagonalize the Hamiltonian

ψn H ψ m δnm ε̄n
F

G


for metals, this condition is difficult to achieve with direct algorithms
E

in metals, actually this optimisation subproblem leads to a linear slowdown

E

with the longest dimension of the (super)cell

G. K RESSE , E LECTRONIC O PTIMISATION Page 6

 Selfconsistency Scheme

trial-charge ρin and trial-wavevectors ψn

HH

HH
M M I
set up Hamiltonian H ρin two subproblems

N



optimization of ψn and ρin



iterative refinements of wavefunctions ψn
I refinement of density:

N


DIIS algorithm
I

new charge density ρout  ∑n fn ψn r 2 P. Pulay, Chem. Phys. Lett. 73,




393 (1980).
refinement of density ρin ρout new ρin
I

refinement of wavefunctions:

N
"


DIIS or Davidson algorithm

no
I JK L KJ I

∆E
J K

K J

Ebreak
J K

K J
KJ
JK

calculate forces, update ions

G. K RESSE , E LECTRONIC O PTIMISATION Page 7

 ALGO flag

ALGO determines how the wavefunctions are optimized

all algorithms are fully parallel for any data distribution scheme
– ALGO= Normal (default): blocked Davidson algorithm
– ALGO= Very Fast: DIIS algorithm
– ALGO= Fast: 5 initial steps blocked Davidson, afterwards DIIS algorithm
after ions are moved: 1 Davidson step, afterwards again DIIS

RMM-DIIS is 1.5 to 2 times faster, but Davidson is more stable

ALGO= Fast is a very reasonable compromise, and should be specified for system
with more than 10-20 atoms

generally the user can not influence the behavior of these algorithms (delicately
optimized black box routines)

G. K RESSE , E LECTRONIC O PTIMISATION Page 8

 OSZICAR and OUTCAR files

POSCAR, INCAR and KPOINTS ok, starting setup

WARNING: wrap around errors must be expected
prediction of wavefunctions initialized
entering main loop
N E dE d eps ncg rms rms(c)
DAV: 1 0.483949E+03 0.48395E+03 -0.27256E+04 96 0.166E+03
DAV: 2 0.183581E+01 -0.48211E+03 -0.47364E+03 96 0.375E+02
DAV: 3 -0.340781E+02 -0.35914E+02 -0.35238E+02 96 0.129E+02
DAV: 4 -0.346106E+02 -0.53249E+00 -0.53100E+00 112 0.158E+01
DAV: 5 -0.346158E+02 -0.52250E-02 -0.52249E-02 96 0.121E+00 0.198E+01
RMM: 6 -0.286642E+02 0.59517E+01 -0.50136E+01 96 0.584E+01 0.450E+00
RMM: 7 -0.277225E+02 0.94166E+00 -0.47253E+00 96 0.192E+01 0.432E+00

initial charge corresponds to the charge of isolated overlapping atoms (POTCAR)

DAV: blocked Davidson algorithm RMM: RMM-DIIS was used
ALGO=F: 5 initial steps blocked Davidson, than RMM-DIIS
4 steps charge fixed, than charge is updated (rms(c) column)
G. K RESSE , E LECTRONIC O PTIMISATION Page 9
 OSZICAR file

N iteration count
E total energy
dE change of total energy
d eps change of the eigenvalues (fixed potential)
ncg number of optimisation steps Hψ
rms total initial residual vector ∑nk wk fnk H εnk ψnk


rms(c) charge density residual vector

G. K RESSE , E LECTRONIC O PTIMISATION Page 10

 OUTCAR file

initial steps (delay no charge update)

cpu time wall clock time
POTLOK: VPU time 0.04: CPU time 0.04 local potential
SETDIJ: VPU time 0.08: CPU time 0.08 set PAW strength coefficients
EDDAV : VPU time 0.94: CPU time 0.94 blocked Davidson
DOS : VPU time 0.00: CPU time 0.00 new density of states
----------------------------------------
LOOP: VPU time 1.06: CPU time 1.06

charge update:
cpu time wall clock time
POTLOK: VPU time 0.04: CPU time 0.04 new local potential
SETDIJ: VPU time 0.09: CPU time 0.09 set PAW strength coefficients
EDDIAG: VPU time 0.14: CPU time 0.14 sub-space rotation
RMM-DIIS: VPU time 0.77: CPU time 0.77 RMM-DIIS step (wavefunc.)
ORTHCH: VPU time 0.01: CPU time 0.02 orthogonalisation
DOS : VPU time -0.01: CPU time 0.00 new density of states
CHARGE: VPU time 0.07: CPU time 0.07 new charge
MIXING: VPU time 0.01: CPU time 0.01 mixing of charge

G. K RESSE , E LECTRONIC O PTIMISATION Page 11

What have all iterative matrix diagonalisation schemes in common ?

one usually starts with a set of trial vectors (wavefunctions) representing the filled
states and a few empty one electron states

ψn n 1 Nbands


 





these are initialized using a random number generator
then the wavefunctions are improved by adding to each a certain amount of the
residual vector
the residual vector is defined as
R ψn H εapp S ψn εapp ψn H ψ n






G


G

G




adding a small amount of the residual vector

ψn ψn λR ψn
O





is in the spirit of the steepest descent approach (usually termed “Jacobi relaxation”)

G. K RESSE , E LECTRONIC O PTIMISATION Page 12

 Iterative matrix diagonalization based on the DIIS algorithm

for our case we need a rather specialized eigenvalue solver

– it should be capable of doing only little work

– efficiency and parallelization are important issues
two step procedure
– start with a set of trial vectors (wavefunctions) representing the electrons
ψn n 1 Nbands (initialized with random numbers)

 




– apply Raighly Ritz optimization in the space spanned by all bands (“sub-space”
rotation)
transform: ψn n 1 Nbands so that

  



ψn H ψ m G δnm ε̄n
F



– then optimize each vector individually ψn n 1 Nbands two or three times


  




G. K RESSE , E LECTRONIC O PTIMISATION Page 13
 “In space” optimization EDDIAG

a set of vectors, that represent the valence electrons ψn n 1 Nbands


 





Raighly Ritz optimization in the space spanned by these vectors (subspace)
search for the unitary matrix Ū such that the wavefunctions ψn


ψn ∑ Ūmn ψm



m

fulfill ψ n H ψm εm δnm
 F

G


this requires the calculation of the subspace matrix H̄

ψn H ψ m H̄mn ψn S ψ m δmn always holds

G

G
F

F






and it’s diagonalisation 
2
the setup of the matrix scales like Nbands NFFT (worst scaling part of VASP)
in the parallel version, communication is required, but modest
worse is the fact that the diagonalisation of H̄mn is not done in parallel

G. K RESSE , E LECTRONIC O PTIMISATION Page 14

 Iterative matrix diagonalization based on the DIIS algorithm

“out of space optimization” EDDRMM

– minimize norm of residual vector using the DIIS method

R ψn H εapp S ψn R ψn R ψn minimal

O G


G


G
F







– each vector is optimized individually (without regard to any other vector)

– easy to implement on parallel computers since each processor handles a subset of
the vectors (no communication required, NPAR=number of proc.)

scaling is propotional to Nbands NFFT with a large prefactor

dominates the compute time for medium to large problems

orthogonalization of wavefunctions ORTHCH

G. K RESSE , E LECTRONIC O PTIMISATION Page 15

 Problem of the DIIS algorithm

eigenstates can be missed for large systems

and there is no clear way to tell when this happens
– in the “best” case no convergence
– but convergence might also slows down after reaching a precision of 10 2 or 10 3

P
– in the worst case, you might not notice anything

in these cases, switch to blocked Davidson (manual contains a number of tricks how you
might be able to use the DIIS algorithm even when it initially fails)

things are not that bad

if the Davidson algorithm is used for the first steps, there is practically no danger of
missing eigenstates

G. K RESSE , E LECTRONIC O PTIMISATION Page 16

 VASP.4.5: new blocked Davidson algorithm

combines “in space” and “out of space” optimization

selects a subset of all bands ψn n 1 Nbands ψk k n1 n2


  

 
"










– optimize this subset by adding the orthogonalized residual vector to the presently
considered subspace
ψk H εapp S ψk k n1 n2


  






– apply Raighly Ritz optimization in the space spanned by these vectors
(“sub-space” rotation in a 2(n2-n1+1) dim. space)
– add additional residuals calculated from the yet optimized bands (“sub-space”
rotation in a 3(n2-n1+1) dim. space)
approximately a factor of 1.5-2 slower than RMM-DIIS, but always stable

available in parallel for any data distribution

G. K RESSE , E LECTRONIC O PTIMISATION Page 17

 charge density mixing (RMM-DIIS)

VASP aims at the minimization of the norm of residual vector

R ρin ρout ρin ρin R ρin min




"





with ρout r ∑occupied wk fnk ψnk r 2
Q

Q




DIIS algorithm is used for the optimization of the norm of the residual vector
linearization of R ρin around ρsc (linear response theory)



Rρ Jρ ρsc






with the charge dielectric function J 
J 1 χ U




 4πe2
q2

leads to
R ρin ρout ρin ρin J ρin ρsc









G. K RESSE , E LECTRONIC O PTIMISATION Page 18

 Divergence of the dielectric function

eigenvalue spectrum of J determines convergence

J 1 χ U


 
4πe2
q2

“broader” eigenvalue spectrum slower convergence

"
for insulators and semi-conductors the width of the eigenvalue spectrum is constant
and system size independent !
for metals the eigenvalue spectrum diverges, its width is proportional to the square of
the longest dimension of the cell:
– short wavelength limit J 1 (no screening)
R

– long wavelength limit J 1 q2 ∝ L2 (metallic screening)

R
!

complete screening in metals causes slow convergence to the groundstate (charge

sloshing)

G. K RESSE , E LECTRONIC O PTIMISATION Page 19

 VASP charge density mixer

0.4
VASP uses a model dielectric function
AMIX
which is a good initial approximation 0.3
for most systems
0.2

J
1 q2
J G1q max
STU WUT 
V V
P

q2
R

STU
X



0.1 AMIN
defaults:
0
AMIX=0.4 ; AMIN=0.1 ; 0 1 2 3 4
BMIX=1.0 2
G (1/A )

this is combined with a convergence accelerator

the initial guess for the dielectric matrix is improved using information accumulated
in each electronic (mixing) step
direct inversion in the iterative subspace (DIIS)

G. K RESSE , E LECTRONIC O PTIMISATION Page 20

 How can you tune VASP to achieve faster convergence

try linear mixing (AMIX=0.1-0.2, BMIX=0.0001)

1
J G1q A


R
VASP also gives information on how good the initial mixing parameters are
allow VASP to run until selfconsistency is achieved and search for the last occurrence
of
eigenvalues of (default mixing * dielectric matrix)
average eigenvalue GAMMA= 2.200
– for linear mixing (e.g. AMIX=0.1 ; BMIX=0.0001) the optimal AMIX is given by
the present AMIX GAMMA
Y

– Kerker like mixing (model dielectric matrix):

GAMMA larger 1 decrease BMIX
E

GAMMA smaller 1 increase BMIX

E

G. K RESSE , E LECTRONIC O PTIMISATION Page 21

 What to do when electronic convergence fails
fails to converge
fails to converge
use Davidson (ALGO=N) ICHARG=12 (no charge update)

converges
converges
play with mixing parameters
converges ICHARG=2
use this setting AMIX=0.1 ; BMIX=0.01

fails to converge

converges
increase BMIX
BMIX=3.0 ; AMIN=0.01

fails to converge
bug report
after positions have been checked

G. K RESSE , E LECTRONIC O PTIMISATION Page 22

 ab initio Molecular dynamics

CP approach exact KS−groundstate

elegant large timestep

simple to implement

problematic for metals, since

direct minimization selfconsistency cycle
electrons must decouple from ionic
degrees of freedom
not the case for metals
problematic for metals very stable
large memory requirements
small timestep efficient for insulators
damped second order and metals
(Tassone, Mauri, Car)
conjugate gradient
(Arias, Payne, Joannopoulos)
RMM−DIIS
(Hutter, Lüthi, Parrinello)

G. K RESSE , E LECTRONIC O PTIMISATION Page 23

 Selfconsistency cycle is very well suited for MDs

MD on the Born Oppenheimer surface (exact KS-groundstate)

selfconsistency cycle determines the dielectric matrix

first time step is rather expensive
but since the dielectric matrix changes only little when ions are moved, the method
becomes very fast in successive steps

wavefunctions and charges etc. are “forward” extrapolated between time-steps

all this makes an extremely efficient scheme that is competitive with the so called
“Car-Parrinello” scheme for insulators
for metals, our scheme is generally much more robust and efficient than the
Car-Parrinello scheme
to select this feature in VASP, set MAXMIX in the INCAR file

G. K RESSE , E LECTRONIC O PTIMISATION Page 24

 Using MAXMIX

usually VASP resets the dielectric matrix to it’s default after moving the ions
but if the ions move only a little bit one can bypass this reset
– definitely a good option for molecular dynamics
– damped molecular dynamics (optimisation)
– works also well during relaxations, if the forces are not large ( 0.5 eV/Å)

L
you need to specify MAXMIX in the INCAR file
set MAXMIX to roughly three times the number of iterations in the first ionic step
the resulting speedups can be substantial (a factor 2 to 3 less electronic steps for each
ionic step)

G. K RESSE , E LECTRONIC O PTIMISATION Page 25

 Using Molecular dynamics

a simple INCAR file

ENMAX = 250 ; LREAL = A # electronic degrees
ALGO = V # very fast algorithm
MAXMIX = 80 # mixing

IBRION = 0 # MD
NSW = 1000 # number ofMD steps
POTIM = 3.0 # time step
TEBEG = 1500 ; TEEND = 500 # target temperature 1500-500 K
SMASS = -1 ; NBLOCK = 50 # scale velocities every 50 steps
SMASS = 2 # use a Nose Hoover thermostat
SMASS = -3 # micro canonical

G. K RESSE , E LECTRONIC O PTIMISATION Page 26

 Using Molecular dynamics

timestep POTIM, depends on the vibrational frequencies and the required energy
conservation
as a rule of thumb: increase POTIM until 3 electronic minisation steps are required per
timestep
another rule of thumb:
H 0.5 fs
increase by 1 fs for each row
Li-F 1 fs
SMASS controls the MD simulation
– SMASS=-3 micro canonical ensemble
– for equilibration and simulated annealing SMASS = -1 ; NBLOCK = 50-100
microcanonical MD, and every NBLOCK steps the kinetic energy is scaled to
meet the requied temperature criterion
– for positive values a Nose Hoover thermostat is introduced

G. K RESSE , E LECTRONIC O PTIMISATION Page 27