[PPT] - Ionic optimisation Georg KRESSE Institut f ur Materialphysik and PowerPoint Presentation

SLIDE 1

Ionic optimisation

Georg KRESSE

Institut f¨ ur Materialphysik and Center for Computational Material Science Universit¨ at Wien, Sensengasse 8, A-1090 Wien, Austria

ienna imulation ackage b-initio

G. KRESSE, IONIC OPTIMISATION

Page 1

SLIDE 2

Overview

the mathematical problem

– minimisation of functions – rule of the Hessian matrix – how to overcome slow convergence

the three implemented algorithms

– Quasi-Newton (DIIS) – conjugate gradient (CG) – damped MD strength, weaknesses

a little bit on molecular dynamics
G. KRESSE, IONIC OPTIMISATION

Page 2

SLIDE 3

The mathematical problem

search for the local minimum of a function f

✁ ✂

x

✄

for simplicity we will consider a simple quadratic function f

✁ ✂

x

✄ ☎

a

✆ ✂

b

✂

x

✆

1 2

✂

xB

✂

x

☎

¯ a

✆

1 2

✁ ✂

x

✝ ✂

x0

✄

B

✁ ✂

x

✝ ✂

x0

✄ ✞

where B is the Hessian matrix Bij

☎

∂2 f ∂xi∂xj

✟

for a stationary point, one requires

✂

g

✁ ✂

x

✄ ☎

∂f ∂

✂

x

☎

B

✁ ✂

x

✝ ✂

x0

✄

gi

✁ ✂

x

✄ ☎

∂f ∂xi

☎

∑

j

Bij

✁

xj

✝

x0

j

✄

at the minimum the Hessian matrix must be additionally positive definite

G. KRESSE, IONIC OPTIMISATION

Page 3

SLIDE 4

The Newton algorithm

educational example

start with an arbitrary start point

✂

x1

calculate the gradient

✂

g

✁ ✂

x1

✄

multiply with the inverse of the Hessian matrix and perform a step

✂

x2

☎ ✂

x1

✝

B

✠

1

✂

g

✁ ✂

x1

✄

by inserting

✂

g

✁ ✂

x1

✄ ☎

∂f ∂

✡

x

☎

B

✁ ✂

x1

✝ ✂

x0

✄

, one immediately recognises that

✂

x2

☎ ✂

x0 hence one can find the minimum in one step

in practice, the calculation of B is not possible in a reasonable time-span, and one

needs to approximate B by some reasonable approximation

G. KRESSE, IONIC OPTIMISATION

Page 4

SLIDE 5

Steepest descent

approximate B by the largest eigenvalue of the Hessian matrix

☛

steepest descent algorithm (Jacobi algorithm for linear equations)

1. initial guess

✂

x1

2. calculate the gradient

✂

g

✁ ✂

x1

✄

3. make a step into the direction of the steepest descent

✂

x2

☎ ✂

x1

✝

1

☞

Γmax

✁

B

✄ ✂

g

✁ ✂

x1

✄

4. repeat step 2 and 3 until convergence is reached

for functions with long steep valleys convergence can be very slow

Γ Γ min

max

G. KRESSE, IONIC OPTIMISATION

Page 5

SLIDE 6

Speed of convergence

how many steps are required to converge to a predefined accuracy

assume that B is diagonal, and start from

✂

x1

☎ ✂

x0

✆

1 B

☎ ✌ ✌

Γ1

✟ ✟ ✟

Γn

✍ ✍ ✂

x1

☎ ✂

x0

✆ ✌ ✌

1

✟ ✟ ✟

1

✍ ✍

with Γ1

✎

Γ2

✎

Γ3

✟ ✟ ✟

gradient

✂

g

✁ ✂

x1

✄

and

✂

x2 after steepest descent step are:

✂

g

✁ ✂

x1

✄ ☎

B

✁ ✂

x1

✝ ✂

x0

✄ ☎ ✌ ✌

Γ1

✟ ✟ ✟

Γn

✍ ✍ ✂

x2

☎ ✂

x1

✝

1 Γn

✂

g

✁ ✂

x1

✄ ☎ ✂

x0

✆ ✌ ✌

1

✝

Γ1

☞

Γn

✟ ✟ ✟

1

✝

Γn

☞

Γn

✍ ✍

G. KRESSE, IONIC OPTIMISATION

Page 6

SLIDE 7

Convergence

the error reduction is given by

✂

x2

☎ ✂

x0

✆ ✌ ✌

1

✝

Γ1

☞

Γn

✟ ✟ ✟

1

✝

Γn

☞

Γn

✍ ✍

1−Γ/Γmax Γ Γ Γ Γ Γ Γ 4 5 3 2 1 1 1−2Γ/Γmax −1

– the error is reduced for each component

– in the high frequency component the error vanishes after on step – for the low frequency component the reduction is smallest

G. KRESSE, IONIC OPTIMISATION

Page 7

SLIDE 8

the derivation is also true for non-diagonal matrices

in this case, the eigenvalues of the Hessian matrix are relevant

for ionic relaxation, the eigenvalues of the Hessian matrix correspond to the

vibrational frequencies of the system the highest frequency mode determines the maximum stable step-width (“hard modes limit the step-size”) but the soft modes converge slowest

to reduce the error in all components to a predefined fraction ε,

k iterations are required 1

✝

Γmin Γmax

k

✏

ε kln 1

✝

Γmin Γmax

✏

lnε

✝

k Γmin Γmax

✏

lnε

✑

k

✏ ✝ ✁

lnε

✄

Γmax Γmin k ∝ Γmax Γmin

G. KRESSE, IONIC OPTIMISATION

Page 8

SLIDE 9

Pre-conditioning

if an approximation of the inverse Hessian matrix is know P

✏

B

✠

1,

the convergence speed can be much improved

✂

xN

✒

1

☎ ✂

xN

✝

λP

✂

g

✁ ✂

xN

✄ ✟

in this case the convergence speed depends on the eigenvalue spectrum of

PB

for P

☎

B

✠

1, the Newton algorithm is obtained

G. KRESSE, IONIC OPTIMISATION

Page 9

SLIDE 10

Variable-metric schemes, Quasi-Newton scheme

variable-metric schemes maintain an iteration history they construct an implicit or explicit approximation of the inverse Hessian matrix B

✠

1 approx

✟

search directions are given by B

✠

1 approx

✂

g

✁ ✂

x

✄ ✟

the asymptotic convergence rate is give by number of iterations ∝ Γmax Γmin

G. KRESSE, IONIC OPTIMISATION

Page 10

SLIDE 11

Simple Quasi-Newton scheme, DIIS

direct inversion in the iterative subspace (DIIS)

set of points

✓ ✂

xi

✔

i

☎

1

✞ ✟ ✟ ✟ ✞

N

✕

and

✓ ✂

gi

✔

i

☎

1

✞ ✟ ✟ ✟ ✞

N

✕

search for a linear combination of xi which minimises the gradient, under the

constraint

∑

i

αi

☎

1

✞

✂

g

✁

∑

i

αi

✂

xi

✄ ☎

B ∑

i

αi

✂

xi

✝ ✂

x0

☎

B ∑

i

αi

✂

xi

✝

∑

i

αi

✂

x0

☎

∑

i

αiB

✁ ✂

xi

✝ ✂

x0

✄ ☎

∑

i

αi

✂

gi

✟

gradient is linear in it’s arguments for a quadratic function

G. KRESSE, IONIC OPTIMISATION

Page 11

SLIDE 12

Full DIIS algorithm

1. single initial point

✂

x1

2. gradient

✂

g1

☎ ✂

g

✁ ✂

x1

✄

, move along gradient (steepest descent)

✂

x2

☎ ✂

x1

✝

λ

✂

g1

3. calculate new gradient

✂

g2

☎ ✂

g

✁ ✂

x2

✄

4. search in the space spanned by

✓ ✂

gi

✔

i

☎

1

✞ ✟ ✟ ✟ ✞

N

✕

for the minimal gradient

✂

gopt

☎

∑αi

✂

gi

✟

and calculate the corresponding position

✂

xopt

☎

∑αi

✂

xi

5. Construct a new point

✂

x3 by moving from

✂

xopt along

✂

gopt

✂

x3

☎ ✂

xopt

✝

λ

✂

gopt

G. KRESSE, IONIC OPTIMISATION

Page 12

SLIDE 13

1. steepest descent step from

✖

x0 to

✖

x1 (arrows correspond to gradients

✖

g0 and

✖

g1)

2. gradient along indicated red line is now know, determine optimal position

✖

x1

pt
3. another steepest descent step form

✖

x1

pt along

✖

gopt

✗ ✖

g

✘ ✖

x1

pt

✙

4. calculate gradient x2

✚

now the gradient is known in the entire 2 dimensional space (linearity condition) and the function can be minimised exactly

pt

x1 x0

1

x2 x x0 a x + a x, a +a =1

0 0 1 1 0 1

G. KRESSE, IONIC OPTIMISATION

Page 13

SLIDE 14

Conjugate gradient

first step is a steepest descent step with line minimisation search directions are “conjugated” to the previous search directions

1. gradient at the current position

✂

g

✁ ✂

xN

✄

2. conjugate this gradient to the previous search direction using:

✂

sN

☎ ✂

g

✁ ✂

xN

✄ ✆

γ

✂

sN

✠

1 γ

☎ ✁ ✂

g

✁ ✂

xN

✄ ✝ ✂

g

✁ ✂

xN

✠

1

✄ ✄✜✛ ✂

g

✁ ✂

xN

✄ ✁ ✂

g

✁ ✂

xN

✠

1

✄ ✄✜✛ ✂

g

✁ ✂

xN

✠

1

✄

3. line minimisation along this search direction

✂

sN

4. continue with step 1), if the gradient is not sufficiently small.

the search directions satisfy:

✂

sNB

✂

sM

☎

δNM

✢

N

✞

M the conjugate gradient algorithm finds the minimum of a quadratic function with k degrees of freedom in k

✆

1 steps exactly

G. KRESSE, IONIC OPTIMISATION

Page 14

SLIDE 15

1. steepest descent step from

✖

x0, search for minimum along

✖

g0 by performing several trial steps (crosses, at least one triastep is required)

✣ ✖

x1

2. determine new gradient

✖

g1

✗ ✖

g

✘ ✖

x1

✙

and conjugate it to get

✖

s1 (green arrow) for 2d-functions the gradient points now directly to the minimum

3. minimisation along search direction

✖

s1

1

x x0 x2

1

x

1

x x0 s1

G. KRESSE, IONIC OPTIMISATION

Page 15

SLIDE 16

Asymptotic convergence rate

asymptotic convergence rate is the convergence behaviour for the case that the

degrees of freedom are much large than the number of steps e.g. 100 degrees of freedom but you perform only 10-20 steps

how quickly, do the forces decrease?
this depends entirely on the eigenvalue spectrum of the Hessian matrix:

– steepest descent: Γmax

☞

Γmin steps are required to reduce the forces to a fraction ε – DIIS, CG, damped MD: Γmax

☞

Γmin steps are required to reduce the forces to a fraction ε Γmax

✞

Γmin are the maximum and minimal eigenvalue of the Hessian matrix

G. KRESSE, IONIC OPTIMISATION

Page 16

SLIDE 17

Damped molecular dynamics

instead of using a fancy minimisation algorithms it is possible to treat the minimisation problem using a simple “simulated annealing algorithm”

regard the positions as dynamic degrees of freedom
the forces serve as accelerations and an additional friction term is introduced
equation of motion (

✂

x are the positions) ¨

✂

x

☎ ✝

2

✤

α

✂

g

✁ ✂

x

✄ ✝

µ˙

✂

x

✞

using a velocity Verlet algorithm this becomes

✂

vN

✒

1

✥

2

☎ ✦ ✁

1

✝

µ

☞

2

✄ ✂

vN

✠

1

✥

2

✝

2

✤

α

✂

FN

✧ ☞ ✁

1

✆

µ

☞

2

✄ ✂

xN

✒

1

☎ ✂

xN

✒

1

✆ ✂

vN

✒

1

✥

2 for µ

☎

2, this is equivalent to a simple steepest descent step

G. KRESSE, IONIC OPTIMISATION

Page 17

SLIDE 18

behaves like a rolling ball with a friction

it will accelerate initially, and then deaccelerate when close to the minimum

if the optimal friction is chosen the ball will glide right away into the minimum
for a too small friction it will overshoot the minimum and accelerate back
for a tool large friction relaxation will also slow down (behaves like a steepest

descent)

x0

G. KRESSE, IONIC OPTIMISATION

Page 18

SLIDE 19

Algorithms implemented in VASP

additional flags termination DISS

IBRION =1 POTIM, NFREE EDIFFG

CG

IBRION =2 POTIM EDIFFG

damped MD

IBRION =3 POTIM, SMASS EDIFFG POTIM determines generally the step size

for the CG gradient algorithm, where line minisations are performed, this is the size of the very first trial step

EDIFFG determines when to terminate relaxation

positive values: energy change between steps must be less than EDIFFG negative values:

✔ ✂

Fi

✔ ✎ ✔✩★✪ ✫✬ ✬ ✭ ✔ ✢

i

☎

1

✞

Nions

G. KRESSE, IONIC OPTIMISATION

Page 19

SLIDE 20

DIIS

POTIM determines the step size in the steepest descent steps

no line minisations are performed !!

NFREE determines how many ionic steps are stored in the iteration history

set of points

✮ ✖

xi

✯

i

✗

1

✰✲✱ ✱ ✱ ✰

N

✳

and

✮ ✖

gi

✯

i

✗

1

✰✲✱ ✱ ✱ ✰

N

✳

searches for a linear combination of xi, that minimises the gradient NFREE is the maximum N

for complex problems NFREE can be large (i.e. 10-20)
for small problems, it is advisable to count the degrees of freedom carefully

(symmetry inequivalent degrees of freedom)

if NFREE is not specified, VASP will try to determine a reasonable value, but

usually the convergence is then slower

G. KRESSE, IONIC OPTIMISATION

Page 20

SLIDE 21

CG

the only required parameter is POTIM

this parameter is used to parameterise, how large the trial steps are

CG requires a line minisations along the search direction

1

x x0 x0

1

x

trial 1

x

trial 2

x

this is done using a variant of Brent’s algorithm – trial step along search direction (conjg. gradient scaled by POTIM) – quadratic or cubic interpolation using energies and forces at

✂

x0 and

✂

x1 allows to determine the approximate minimum – continue minimisation as long as approximate minimum is not accurate enough

G. KRESSE, IONIC OPTIMISATION

Page 21

SLIDE 22

Damped MD

two parameters POTIM and SMASS

✂

vN

✒

1

✥

2

☎ ✦ ✁

1

✝

µ

☞

2

✄ ✂

vN

✠

1

✥

2

✝

2

✤

α

✂

FN

✧ ☞ ✁

1

✆

µ

☞

2

✄ ✂

xN

✒

1

☎ ✂

xN

✒

1

✆ ✂

vN

✒

1

✥

2 α ∝ POTIM and µ ∝ SMASS

POTIM must be as large as possible, but without leading to divergence

and SMASS must be set to µ

✏

2 Γmin

☞

Γmax, where Γmin and Γmax are the minimal und maximal eigenvalues of the Hessian matrix

a practicle optimisation procedure:

– set SMASS=0.5-1 and use a small POTIM of 0.05-0.1 – increase POTIM by 20 % until the relaxation runs diverge – fix POTIM to the largest value for which convergence was achieved – try a set of different SMASS until convergence is fastest (or stick to SMASS=0.5-1.0)

G. KRESSE, IONIC OPTIMISATION

Page 22

SLIDE 23

Damped MD — QUICKMIN

alternatively do not specify SMASS (or set SMASS

✎

0) this select an algorithm sometimes called QUICKMIN

QUICKMIN

✂

vnew

☎ ✴ ✵ ✶

α

✂

F

✆ ✂

F

✁ ✂

vold

✛ ✂

F

✄ ☞ ✔ ✔ ✂

F

✔ ✔

for

✂

vold

✛ ✂

F

✷

α

✂

F else – if the forces are antiparallel to the velocities, quench the velocities to zero and restart – otherwise increase the “speed” and make the velocities parallel to the present forces

I have not often used this algorithm, but it is supposed to be very efficient
G. KRESSE, IONIC OPTIMISATION

Page 23

SLIDE 24

Damped MD — QUICKMIN

my experience is that damped MD (as implemented in VASP) is faster than QUICKMIN but it requires less playing around defective ZnO surface: 96 atoms are allowed to move! relaxation after a finite temperature MD at 1000 K

20 40 60 80 steps

6
4
2

2 log(E-E0) damped: SMASS=0.4 quickmin

G. KRESSE, IONIC OPTIMISATION

Page 24

SLIDE 25

Why so many algorithms :-(... decision chart

close to minimum 1−3 degrees of freedom Really, this is too complicated no yes no very broad vib. spectrum >20 degrees of freedom yes yes no no yes DIIS damped MD or QUICKMIN CG

G. KRESSE, IONIC OPTIMISATION

Page 25

SLIDE 26

Two cases where the DIIS has huge troubles

X0 X1 DIIS is dead, since it consideres the forces only in cartesian coordinates the Hessian matrix changes force increases along the search direction it will move uphill instead of down when the octahedron rotates! rigid unit modes i.e. in perovskites (rotation) molecular systems (rotation)

G. KRESSE, IONIC OPTIMISATION

Page 26

SLIDE 27

How bad can it get

the convergence speed depends on the eigenvalue spectrum of the Hessian matrix

– larger systems (thicker slabs) are more problematic (acoustic modes are very soft) – molecular system are terrible (week intermolecular and strong intramolecular forces) – rigid unit modes and rotational modes can be exceedingly soft the spectrum can vary over three orders of magnitudes

✑

100 or even more steps might be required ionic relaxation can be painful

to model the behaviour of the soft modes, you need very accurate forces since
therwise the soft modes are hidden by the noise in the forces

EDIFF must be set to very small values (10

✠

6) if soft modes exist

G. KRESSE, IONIC OPTIMISATION

Page 27

SLIDE 28

Electronic optimization

Georg KRESSE

Institut f¨ ur Materialphysik and Center for Computational Material Science Universit¨ at Wien, Sensengasse 8, A-1090 Wien, Austria

ienna imulation ackage b-initio

G. KRESSE, ELECTRONIC OPTIMISATION

Page 1

SLIDE 29

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. KRESSE, ELECTRONIC OPTIMISATION

Page 2

SLIDE 30

Density functional theory according to Kohn-Sham

density and kinetic energy: sum of one electron charge densities and kinetic energies ρtot

✁

r

✂ ✄

2 Ne

☎

2

∑

n

✆

1

✝

ψn

✁

r

✂ ✝

2

✞

ρion

✁

r

✂ ✟

Ne number of electrons

✠

¯ h2 2me 2

Ne

☎

2

∑

n

✆

1 ψn

✁

r

✂☛✡

∇2ψn

✁

r

✂

d3r

☞ ✌✍ ✎

kinetic energy

✞

1 2 ρtot

✁

r

✂

ρtot

✁

r

✏ ✂ ✝

r

✠

r

✏ ✝

d3rd3r

✏ ☞ ✌ ✍ ✎

electrost. energy

✞

Exc

✑

ρ

✁

r

✂ ✒ ☞ ✌ ✍ ✎

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

G. KRESSE, ELECTRONIC OPTIMISATION

Page 3

SLIDE 31

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)

Gradient: Fn

✁

r

✂ ✄ ✠

¯ h2 2me ∇2

✞

V eff

✁

r

✟ ✢

ψn

✁

r

✏ ✂ ✣ ✂ ✠

εn ψn

✁

r

✂

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

✕ ✠

¯ h2 2me ∇2

✞

V eff

✁

r

✟ ✢

ρ

✁

r

✏ ✂ ✣ ✂

ψn

✁

r

✂ ✄

εnψn

✁

r

✂

n

✄

1

✟✥✤ ✤ ✤ ✟

Ne

✦

2 as a result one obtains a new charge density ρ

✔

r

✕ ✗

∑n

✖

ψn

✔

r

✕✖

2 and a new Schr¨

dinger

equation

✧

iteration

G. KRESSE, ELECTRONIC OPTIMISATION

Page 4

SLIDE 32

iteration log 10 E-E0

8
6
4
2

2 log 10 |F-F 0 |

4
3
2
1

1 10 20 30 40 iteration log 10 E-E0

8
6
4
2

log 10 |F-F 0 |

4
3
2
1

1 5 10 15 20

★ ★ ★ ★ ★ ★ ✩ ✩ ✩ ✩ ✩ ✩ ✪ ✪ ✪ ✪ ✪ ✪ ✫ ✫ ✫ ✫ ✫ ✫ ✬ ✬ ✬ ✬ ✬ ✬ ✭ ✭ ✭ ✭ ✭ ✭ ✮ ✮ ✮ ✮ ✮ ✮ ✯ ✯ ✯ ✯ ✯ ✯ ✰ ✰ ✰ ✰ ✰ ✰ ✱ ✱ ✱ ✱ ✱ ✱ ✲ ✲ ✲ ✲ ✲ ✲ ✳ ✳ ✳ ✳ ✳ ✳ ✴ ✴ ✴ ✴ ✴ ✴ ✵ ✵ ✵ ✵ ✵ ✵ ✶ ✶ ✶ ✶ ✶ ✶ ✷ ✷ ✷ ✷ ✷ ✷ ✸ ✸ ✸ ✸ ✸ ✸ ✹ ✹ ✹ ✹ ✹ ✹ ✺ ✺ ✺ ✺ ✺ ✺ ✻ ✻ ✻ ✻ ✻ ✻ ✼ ✼ ✼ ✼ ✼ ✼ ✽ ✽ ✽ ✽ ✽ ✽ ✾ ✾ ✾ ✾ ✾ ✾ ✿ ✿ ✿ ✿ ✿ ✿ ❀ ❀ ❀ ❀ ❀ ❀ ❁ ❁ ❁ ❁ ❁ ❁ ❂ ❂ ❂ ❂ ❂ ❂ ❃ ❃ ❃ ❃ ❃ ❃ ❄ ❄ ❄ ❄ ❄ ❄ ❅ ❅ ❅ ❅ ❅ ❅ ❆ ❆ ❆ ❆ ❆ ❆ ❇ ❇ ❇ ❇ ❇ ❇ ❈ ❈ ❈ ❈ ❈ ❈ ❉ ❉ ❉ ❉ ❉ ❉

n=8 n=1 n=1 n=4

self.consistent

n=8 n=1 L n=1,2,4,8

forces self.consistent Car−Parrinello direct Car−Parrinello direct disordered fcc Fe, metal disordered diamond, insulator energy

G. Kresse and J. Furthm¨

uller, Phys. Rev. B 54, 11169 (1996).

G. KRESSE, ELECTRONIC OPTIMISATION

Page 5

SLIDE 33

Direct minimization (not supported by vasp.4.5)

preconditioned conjugate gradient algorithm

was applied Gradient: Fn

✁

r

✂ ✄ ✠

¯ h2 2me ∇2

✞

V eff

✁

r

✟ ✢

ψn

✁

r

✏ ✂ ✣ ✂ ✠

εn ψn

✁

r

✂

the main troubles are

– to keep the set of wavefunctions

✢

ψn

✁

r

✂ ✝

n

✄

1

✟ ✤ ✤ ✟

Ne

✦

2

✣

rthogonal

– “sub-space” rotation

❊

at the end one aims to have a set of wavefunction

✢

ψn

✁

r

✂ ✝

n

✄

1

✟ ✤ ✤ ✟

Nbands

✣

that diagonalize the Hamiltonian

❋

ψn

✝

H

✝

ψm

✄

δnm¯ εn

❊

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

❊

in metals, actually this optimisation subproblem leads to a linear slowdown with the longest dimension of the (super)cell

G. KRESSE, ELECTRONIC OPTIMISATION

Page 6

SLIDE 34

Selfconsistency Scheme

trial-charge ρin and trial-wavevectors ψn

❍ ❍ ❍ ❍ ■

set up Hamiltonian H

✁

ρin

✂ ■

iterative refinements of wavefunctions

✢

ψn

✣ ■

new charge density ρout

✄

∑n fn

✝

ψn

✁

r

✂ ✝

2

■

refinement of density ρin

✟

ρout

✧

new ρin

■ ❏ ❏ ❏ ❏ ❏ ❑ ❑ ❑ ❑ ❑ ❑ ❑ ❑ ❑ ❑ ❏ ❏ ❏ ❏ ❏

no ∆E

▲

Ebreak calculate forces, update ions

▼ ■ ▼ ◆

two subproblems

ptimization of

✓

ψn

✜

and ρin

◆

refinement of density: DIIS algorithm

P. Pulay, Chem. Phys. Lett. 73,

393 (1980).

◆

refinement of wavefunctions: DIIS or Davidson algorithm

G. KRESSE, ELECTRONIC OPTIMISATION

Page 7

SLIDE 35

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
ptimized black box routines)
G. KRESSE, ELECTRONIC OPTIMISATION

Page 8

SLIDE 36

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. KRESSE, ELECTRONIC OPTIMISATION

Page 9

SLIDE 37

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. KRESSE, ELECTRONIC OPTIMISATION

Page 10

SLIDE 38

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

rthogonalisation

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. KRESSE, ELECTRONIC OPTIMISATION

Page 11

SLIDE 39

What have all iterative matrix diagonalisation schemes in common ?

ne 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

✠

εappS

✂ ✝

ψn

εapp

✄ ❋

ψn

✝

H

✝

ψn

adding a small amount of the residual vector

ψn

❖

ψn

✞

λR

✁

ψn

✂

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

G. KRESSE, ELECTRONIC OPTIMISATION

Page 12

SLIDE 40

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

✄

δnm¯ εn – then optimize each vector individually

✢

ψn

✝

n

✄

1

✟ ✤ ✤ ✟

Nbands

✣

two or three times

G. KRESSE, ELECTRONIC OPTIMISATION

Page 13

SLIDE 41

“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 ¯ U such that the wavefunctions ψ

✏

n

ψ

✏

n

✄

∑

m

¯ Umnψm fulfill

❋

ψ

✏

n

✝

H

✝

ψ

✏

m

✄

εmδnm this requires the calculation of the subspace matrix ¯ H

❋

ψn

✝

H

✝

ψm

✄

¯ Hmn

✁ ❋

ψn

✝

S

✝

ψm

✄

δmn always holds

✂

and it’s diagonalisation

the setup of the matrix scales like N2

bandsNFFT (worst scaling part of VASP)

in the parallel version, communication is required, but modest worse is the fact that the diagonalisation of ¯ Hmn is not done in parallel

G. KRESSE, ELECTRONIC OPTIMISATION

Page 14

SLIDE 42

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

✠

εappS

✂ ✝

ψn

❋

R

✁

ψn

✂ ✝

R

✁

ψn

✂

❖

minimal – 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 NbandsNFFT with a large prefactor

dominates the compute time for medium to large problems

rthogonalization of wavefunctions

ORTHCH

G. KRESSE, ELECTRONIC OPTIMISATION

Page 15

SLIDE 43

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

P

2 or 10

P

3 – 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. KRESSE, ELECTRONIC OPTIMISATION

Page 16

SLIDE 44

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

✠

εappS

✂

ψ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. KRESSE, ELECTRONIC OPTIMISATION

Page 17

SLIDE 45

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

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. KRESSE, ELECTRONIC OPTIMISATION

Page 18

SLIDE 46

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) – long wavelength limit J

❘

1

✦

q2 ∝ L2 (metallic screening) complete screening in metals causes slow convergence to the groundstate (charge sloshing)

G. KRESSE, ELECTRONIC OPTIMISATION

Page 19

SLIDE 47

VASP charge density mixer

VASP uses a model dielectric function

which is a good initial approximation for most systems J

P

1

❘

G1

q

✄

max

✁

q2

❙❚❯ ❱

q2

✞ ❲ ❚❯ ❱ ✟ ❙❚❯ ❳ ✂

defaults:

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

2)

0.1 0.2 0.3 0.4 J AMIN AMIX

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. KRESSE, ELECTRONIC OPTIMISATION

Page 20

SLIDE 48

How can you tune VASP to achieve faster convergence

try

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

P

1

❘

G1

q

✄

A

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

f

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

– Kerker like mixing (model dielectric matrix):

❊

GAMMA larger 1

❖

decrease BMIX

❊

GAMMA smaller 1

❖

increase BMIX

G. KRESSE, ELECTRONIC OPTIMISATION

Page 21

SLIDE 49

What to do when electronic convergence fails

use Davidson (ALGO=N) use this setting fails to converge converges ICHARG=12 (no charge update) ICHARG=2 AMIX=0.1 ; BMIX=0.01 converges increase BMIX BMIX=3.0 ; AMIN=0.01 converges fails to converge fails to converge bug report after positions have been checked fails to converge play with mixing parameters converges

G. KRESSE, ELECTRONIC OPTIMISATION

Page 22

SLIDE 50

ab initio Molecular dynamics

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

exact KS−groundstate direct minimization selfconsistency cycle CP approach

G. KRESSE, ELECTRONIC OPTIMISATION

Page 23

SLIDE 51

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,

ur 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. KRESSE, ELECTRONIC OPTIMISATION

Page 24

SLIDE 52

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/ ˚ A)

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. KRESSE, ELECTRONIC OPTIMISATION

Page 25

SLIDE 53

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. KRESSE, ELECTRONIC OPTIMISATION

Page 26

SLIDE 54

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 Li-F 1 fs increase by 1 fs for each row

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. KRESSE, ELECTRONIC OPTIMISATION

Page 27