Simula'ons of atomic scale systems using HPC Hannes Jnsson, Faculty - - PowerPoint PPT Presentation

Simula'ons of atomic scale systems using HPC

Hannes Jónsson, Faculty of Physical Sciences, University of Iceland

A few historical remarks … Computa(onal and theore(cal chemistry (electronic structure and rearrangements of atoms/spins)

• ApplicaFons to a few systems
• Development of algorithms and computer implementaFons

Computer clusters available to researchers at the University of Iceland

2003 Bjólfur (Beowoulf cluster, 130 computers) located in the basement of VR-III Funded in part by the Icelandic Research Fund (or rather its precursor). 2008 Sól (Sun system, 320 cores, each node 2x4 cores) located in the basement of VR-III Funded in part by the Icelandic Research Fund. 2009 Jötunn (IBM, 168 cores) located at Reiknistofnun in Taeknigardur Purchased by funds from Reiknistofnun. 2015 Garpur (IBM, 12 and 16 core nodes) located at Reiknistofnun in Neshagi Funded in part by the Icelandic Research Fund. 2011 Nordic HPC in Iceland (HP system, 12 core nodes) located at Thor Data Center Purchased by funds from Denmark, Norway, and Sweden. 2017 Garpur2 (??, 100 teraflop) located at Reiknistofnun in Neshagi Funded in part by the Icelandic Research Fund.

Examples of collabora'ons with computer scien'sts

at University of Washington in SeaAle: Ongoing collabora'on with computer scien'sts at Aalto University in Finland:

• `A Parallel ImplementaFon of the Car-Parrinello Method by Orbital DecomposiFon', J. Wiggs

and H. Jónsson, Computer Physics CommunicaFons, vol. 81, p. 1 (1994).

• `A Hybrid DecomposiFon Parallel ImplementaFon of the Car-Parrinello Method', J. Wiggs

and H. Jónsson, Computer Physics CommunicaFons, vol. 87, p. 319 (1995).

• `Dynamic-Domain-DecomposiFon Parallel Molecular Dynamics', S. Srivilliputhur, I. Ashok,
• H. Jónsson, G. Kalonji and J. Zarhojan, Comp. Phys. Commun., 102, 44 (1997).
• `Parallel Short-Range Molecular Dynamics Using Adhara RunFme System', - S. Srivilliputhur,
• I. Ashok, H. Jónsson, G. Kalonji and J. Zarhojan, Comp. Phys. Commun., 102, 28 (1997).
• 'Minimum energy path calculaFons with Gaussian process regression', O-P. KoisFnen,
• E. Maras, A. Vehtari, H. Jónsson. Nanosystems: Physics, Chemistry, MathemaFcs, 7, 925 (2016).
• 'Nudged ElasFc Band CalculaFons Accelerated with Gaussian Process Regression’ O-P. KoisFnen,
• F. Dagbjartsdójr, V. Ásgeirsson, A. Vehtari and H. Jónsson. J. Chem. Phys. (in press).
• 'Parallel implementaFon of gamma-point plane-wave DFT with exact exchange', E. J. Bylaska,
• K. Tsemekhman, S. B. Baden, J. H. Weare and H. Jónsson, J. Comp. Chem. 32, 5469 (2011).

at University of California San Diego:

ComputaFonal chemistry (reikniefnafræði)

Use basic laws of physics to calculate properFes of chemicals and materials Predict properFes of new systems and/or help interpret experiments on exisFng systems Tremendous progress in recent years

• improved approximaFons to the basic equaFons
• improved algorithms
• improved implementaFons
• increased computer power

Basic tasks in computaFonal chemistry

Typically work within the Born-Oppenheimer (adiabaFc) approximaFon. Two step process:

• 1. For fixed locaFon of the nuclei of the atoms, solve for the electronic

structure, i.e. find the distribuFon and energy of the electrons. Energy surface for the locaFon of the atoms and orientaFon of magneFc moment, E(R,ω). Need to solve many electron Schrödinger equaFon or equivalent formulaFon (most commonly now density funcFonal theory, DFT).

• 2. Move the atoms and rotate spins given the energy from step 1, i.e. on the energy surface.

Typically, use classical descripFon of the atoms, but someFmes need to include quantum mechanical descripFon (for example tunneling of atoms or spins). In principle, can solve Newton equaFon to calculate Fme evoluFon (trajectories), but the huge difference in Fme scale between atomic vibraFons and chemical reacFons

• r other atomic or spin rearrangements between stable states make it necessary to use

rate theory (such as transiFon state theory).

Self-interacFon error in standard implementaFons of DFT

In Kohn-Sham DFT, the energy is esFmated form

E[ρ(r)] =

N

X

i

Z d3r|rφi(r)|2 2 + Z d3rρ(r)vext(r) + 1 2 Z d3rd3r0 ρ(r)ρ(r0) |r r0| + Exc[ρ(r)]

Orbitals,

φi(r)

, introduced only to obtain an accurate esFmate of kineFc energy. Total electron density is a sum over orbital densiFes (Note: unitary invariance)

Ekin Eext

In terms of orbital densiFes, the KS esFmate of the electron Coulomb interacFon is

ρ(r) = X

i

ρi(r) = X

i

|φi(r)|2

It includes interacFon of the each orbital density with itself (i=j terms). A more accurate esFmate, which was used by Hartree, excludes the diagonal terms

EKS

C

EKS

C

= 1 2 X

i

X

j

Z d3r Z d3r0 ρi(r)ρj(r0) |r − r0| Coulomb self-interacFon correcFon

EC

H = 1

2 X

i

X

j6=i

Z d3r Z d3r0 ρi(r)ρj(r0) |r − r0| = EKS

C

− 1 2 X

i

Z d3r Z d3r0 ρi(r)ρi(r0) |r − r0|

Perdew-Zunger self-interacFon correcFon, PZ-SIC

should fix all approximaFons in the kineFc and Coulomb energy, but pracFcal funcFonals are approximate and a self-interacFon error remains.

Exc

EKS

C

can be corrected explicitly by subtracFng the diagonal terms, but

Exc

then also needs to be corrected. For a one-electron system, self-interacFon in a KS funcFonal can be removed by

EKSSIC[ρ1, . . . , ρN] = EKS[ρ] −

N

X

i

1 2 Z d3r Z d3r0 ρi(r)ρi(r0) |r − r0| − EKS

xc [ρi]

• Perdew and Zunger (PRB 1981)

proposed using this correcFon for many electron systems. PZ-SIC Note: A specific set of orbitals minimizes the total energy given by due to orbital density dependence (ODD). Convenient to work with two sets of orbitals

diagonalize the KS Hamiltonian minimize the ODD energy

The two sets are related by a unitary transformaFon, W

Canonical orbitals (give energy) OpFmal orbitals (give locaFon)

EKS−SIC Now need a new inner loop in SCF to find the W that minimizes EKS-SIC.

eigenstates of the Hamiltonian

• btained by minimizing the total energy

VariaFonal implementaFon of PZ-SIC

The energy is no longer unitary invariant with respect to the orbitals, now need an inner loop to find opFmal linear combinaFon

• f the orbitals at each SCF iteraFon.

Use algorithm developed in signal processing by Abrudan,Eriksson and Koivunen (Signal Processing 2009, 89, 1704). Same algorithm can be used to find opFmal local orbitals as post processing of Kohn-Sham or Hartree-Fock orbitals (S. Lehtola, E. Jónsson and HJ., JCTC 2013, 2014, 2017) Has been implemented in the GPAW code, a separate branch. Developed originally to use real space grid but can now also make use of atomic basis sets and plane waves. PAW representaFon of inner electrons. Elvar

Example applicaFon: Energy of two Mn atoms as a funcFon

• f the distance between them (Tushar, Aleksei, Elvar).

Two magneFc states: FerromagneFc (green) and anFferromagneFc (red)

Standard DFT methodology (PW91 funcFonal) fails, while self-interacFon corrected funcFonal gives results in close agreeement with high level quantum chemistry calculaFons and experiments.

Electronic structure calculaFons are computaFonally demanding

Can readily deal with 100 to 200 atoms with DFT. ComputaFonal effort scales as N3 (where N is number of electrons) Much larger systems can be simulated if an energy surface, E(R,ω), can be developed. Example: DescripFon of interacFon between H2O molecules using mulFpole expansion of the electrostaFcs, including dipole and octupole polarizaFon (SCME potenFal funcFon)

QM/MM approach: Divide the system up into a core region that is described by DFT, and an outer region that is described by SCME. ElectrostaFc embedding implemented in GPAW code Asmus and Elvar

Second step: How do the atoms move in chemical reac'ons (or other rearrangements such as diffusion)?

• A transiFon with an energy barrier of 0.5 eV and

a typical pre-exponenFal factor occurs 1000 Fmes per second at room temperature – fast on laboratory scale!

• A video of a direct classical dynamics simulaFon

where each vibraFon spans a second in the video would go on for more than 100 years in between such reacFve events – slow on atomic scale!

0.5 eV 1000/s Most interesting transitions are rare events (i.e., much slower than vibrations).

Time scale problem:

Typically there is a clear separation of time scales, and a statistical approach can be used

IdenFfy a 3N-1 dimensional dividing surface, that represents a boAleneck for going from the iniFal to a final state:

Transition State Theory (Wigner, Eyring 1930s)

3N-1 dimensional dividing surface, . Add thickness σ to define TransiFon State (TS)

The bottleneck can be due to an energy barrier and/or entropy barrier Initial state Final state

‡

σ

‡

HTST - Harmonic approximation to TST:

Good for solids at not too high T >kBT

R

1st order SP 2nd order SP

Works well when (1) energy of second order saddle points is significantly higher than kBT over the energy of first order saddle points, and (2) when the potenFal is smooth enough that a second order Taylor approximaFon to the PES is good enough in the region with large staFsFcal weight.

Energy ridge

Approximate the energy surface with second order Taylor expansions, (a) For reactant region expand around the local energy minimum, (b) For the transiFon state expand around the 1st order saddle point.

Nudged Elastic Band (NEB) Method

Initial state Final state

and images are distributed with springs Effective force on each image (R here denotes atom coordinates, was x before): where the perpendicular force is

tangent along current path

Initial guess

Converged MEP

Initial path (here linear)

! F

i s ≡ ki+1

! R

i+1 −

! R

i

( ) − ki

! R

i −

! R

i−1

( )

! F

i nudged = −

! ∇ V( ! R

i) ⊥ +

! F

i s ⋅ ˆ

τ

||

( )ˆ

τ

||

! ∇ V( ! R

i) ⊥ =

! ∇ V( ! R

i) −

! ∇ V( ! R

i)⋅ ˆ

τ

||

( )ˆ

τ

||

(Mills, Jónsson & Schenter, Surf. Sci. 1995; Henkelman & Jónsson, JCP. 2000)

Create several images (discretization pts.) Minimize all images simultaneously, in parallel. Typically 5 to 10 images. Estimate the tangent at each image using segment to adjacent image at higher energy.

Electrochemical reducFon of CO2

*CO + (H+ + e-) *COH *CO + (H+ + e-) *CHO *COH + (H+ + e-) *C + H2O *COH + (H+ + e-) *CHOH

S S IS IS S S

*CO + (H+ + e-) *COH *COH + (H+ + e-) *CHOH

(a) (b) (c)

• 20
• 15
• 10
• 5

5

• 1

1 2 3 4 5 6 7 8

Free energy / eV (H+ + e-) transferred

0.65 0.28 0.41 0.72 0.08 CO2 *COOH *CO + H2O *COH *CHOH *CH2 *CH3 CH4

U = 0 V No H2O U = -1.3 V With H2O

0.26 0.81 0.72 0.34 0.11 *CH2OH 0.18 *C + H2O *CHO CH2O(g) *CH3O *CH4 + *O *OH * + H2O 0.45 0.09 0.01 0.04

Tafel reaction Heyrovsky reaction

CH3OH *C2H4OH C2H4

1 2 3 4 5 6 7 8

(H+ + e-) transferred

C2H5 OH 10-4 10-2 100 102 104 10-4 10-2 100 102 104 log (i / mA cm-2) - Theory log (i / mA cm-2) - Experiment CH4 C2H4 C2H5OH CO CH3OH

The simulaFons explain why only copper electrodes give a range of products, including methane, ethene, and ethanol. Now being used to search for a beuer catalyst (with Javed Hussain and Egill Skúlason).

Methods for finding first order saddle points Two categories:

• A. Two point problem – both initial and final state minima

are known.

• B. One point problem – only initial state minimum is known.

R R P Easier, can use info about final state minimum to guide the search Harder, can only use local info about the energy surface

Find saddle points on a high-dimensional objective function surface

gradient of o the objecFve funcFon Transformed (effecFve) gradient TransformaFon

After transforming the gradient, a SP search is reformulated as a minimization using

• Geff. Any method that can find the zero of a gradient can the be applied to find

a first order saddle point (damped dynamics, conjugate gradients, L-BFGS …)

• 0.5
• 0.4
• 0.3
• 0.2
• 0.1

0.2 0.3 0.4 0.5 0.6 0.7

• 0.5
• 0.4
• 0.3
• 0.2
• 0.1

0.2 0.3 0.4 0.5 0.6 0.7

−G −Geff

• f the gradient

Minimum mode following method (G. Henkelman and HJ, 1999)

Transform the force by inverting the component along the minimum mode, , of the Hessian (the matrix of second derivatives) Within HTST, need to find all the relevant saddle points on the energy rim surrounding the current state minimum. The direction along the minimum mode is found by minimizing the energy of a dimer (two replicas of the system, ) or by using the Davidson algorithm. No need to construct the Hessian matrix. Use some minimization algorithm that only requires derivative of the objective function (not the objective function itself) and it will converge on a first order saddle point. The force projection locally transforms a first order saddle point to a minimum. Make a displacement

from the minimum, then follow effective gradient 1 2 3 4

Use random initial displacement and then climb up the PES

Liule or no bias from preconceived noFon of the mechanism, perhaps displace under-coordinated atoms and their neighbors. – Can discover unexpected mechanism and final state(s).

Two phases:

• 1. When lowest eigenvalue of H is posiFve,

move along minimum mode.

• 2. Axer lowest eigenvalue of H becomes negaFve, include also

force perpendicular to the minimum mode.

• No need to create a table of

transitions and their rates before the simulation as in regular KMC Adaptive Kinetic Monte Carlo G. Henkelman and HJ, JCP. (2001)

1. Find low energy saddle points near the current minimum using multiple SP searches started at random (Gaussian distributed) 2. Find the prefactor, ν, from normal mode analysis and calculate the rate of each process 3. Use a random number to pick one of the transitions according to the relative transition rates 4. Advance system to the final state of the chosen process (trajectory or slide down) 5. Increment time by an amount Δt = 1/Σ ri Repeat 1 - 5 until the desired time interval has been reached

T k E i

B i

e v r

Δ −

=

• Also, no need to assign atoms to lattice sites

(defects, glass, …)

EON Distributed soxware for long Fme scale simulaFons

• Distribute saddle point searches
• ver internet

– Communication builds on BOINC – Communication and computation are not entangled – Client runs as stand-alone

• Atomic Forces

– Can use DFT or potential functions for atomic forces

• Implementation

– Client side C++ – Server side Python

(Andreas Pedersen and H.J., ‘Distributed ImplementaFon of the AdapFve KineFc Monte Carlo Method‘, Mathematics and Computers in Simulation, 2009; Sam Chill et al., Modelling and Simula(on in Mat. Sci. & Eng., 22, 055002 (2014))

Several tricks to improve efficiency: Skipping path, systemaFc coarse graining, recycling, …

Structure of metal nanoparFcles

Long Fme scale simulated annealing, move from one configuraFon of atoms to another by finding saddle points on the energy surface. EON soxware for distributed compuFng.

• 0.6

0.4 0.0 Energy [eV] 2 4 6 Reaction coordinate [A]

• 1.0

Atom index [#] 20 40 60 80 100 120 140 0.0 1.0 2.0 3.0 4.0 R e a c t a n t Atomic displacement [A]

• Saddle

Product S a d d l e P r

• d

u c t

• 1.0

1.0 0.0

• 2.0

Energy [eV] 10 20 30 Reaction coordinate [A]

• R

e a c t a n t S a d d P r

• A hole in the center of a 561 atom

Au cluster forms as a row of atoms moves out and an addiFonal atom enters the outermost shell TransformaFon from FCC to icosahedral str.

Coarse Graining Problems arise when two or

more states are connected by low energy barriers while barriers leading out of this group of states are several kBT higher (as in regular KMC simulations)

Coarse graining

– Visited minima are assigned a reference energy value – At each visit the reference value is increased – A composite state is created when reference value exceeds the saddle point energy – Composite states can merge

Accurate estimate of time spent in coarse grained state and selection of exit state by using absorbing Markov chain theory (Novotny: http://arxiv.org/abs/cond-mat/ 0109182)

energy

• 5
• 4
• 3
• 2
• 1

1 2 3 4 5

5 10 15 20 25 30 35 1 2 3 4 5 6 7 8 9 x 10

7

d, number of dimensions

• No. of function evaluations

GOUST

• 4.5 × 102 · d 2.1

GAmi n

• 1.2 × 104 · d 2.6

GOUST (global optimization using saddle traversals) Results for Lunacek function

24

5 10 15 20 25 30 35

d, number of dimensions

5 10 15 20 25 30 35 1 2 3 4 5 6 7 8

d, number of dimensions Log of No. of function evaluations

GAmi n GOUST

• M. Lunacek, D. Whitley and A. Suuon,

Parallel Problem Solving from Nature, Springer Verlag, 2008,

• pp. 498

f(x) = min d X

i

(xi − µ1)2, d + s

d

X

i

(xi − µ2)2 ! + 10

d

X

i

(1 − cos 2π(xi − µ1))

µ1 = 2.5, µ2 = −2.5, s = 0.7

Extend transiFon state theory to esFmate lifeFme of magneFc states

(Bessarab, Uzdin, Jónsson 2012).

Thermal Stability of MagneFc States

Apply to studies of various magneFc systems, including non-collinear. Example quesFon: How can localized magneFc states be made stable enough, and how can they be manipulated to store and read informaFon?

E, ✏

π 2π θ2 π 2π θ3 1 2 3 4 5 6 P SP AP

0.0 0.4 0.8 1.2 1.6 2.0 ⇡/2 ⇡ 3⇡/2

E, ✏

q

✓2

2 + ✓2 3

P SP AP

MagneFc Tip / Surface InteracFon

(b)

V

(1)

V

(2)

(a)

80 60 40 20 20 40 0.2 0.4 0.6 0.8 1 Energy (meV) Reaction coordinate Only tip B = 0 T B = 0.4 T

B

23 meV 33 meV

10-7 100 107 5 10 15 20 25 30 35 Lifetime (s) N 10 100 1000 26 28 30 P state AP state Only tip B = 45 mT Expt. B = 0 T

Aleksei

Quantum mechanical description of the atoms

Use statistical Feynman path integrals:

classical stat. mech. qantum stat. mech.

ks = mP ✓2πkBT h ◆2

The statistical mechanical partition function of a quantum particle is mathematically the same as the partition function of a closed chain

• f replicas connected by temperature dependent springs, ring polymer.

Should take the limit as P goes to infinity, i.e. infinitely many replicas of the system. The distribution of replicas represents quantum delocalization. They get pulled inn to a point as mass or T become large. V (~ r) → V eff(~ r1,~ r2, . . . ,~ rP ) =

P

X

i

• ks(~

ri − ~ ri−1)2 + V (~ ri)/P

For a spin of length s, with orientaFon , the acFon is

Quantum mechanical descrip'on of magne'c moments Onset temperature for tunneling

S(θ, φ) = Z β/2

−β/2

dτ h −is(1 − cos θ) ˙ φ + U(θ, φ) i

where U(θ, φ) is the potenFal energy surface. gives classical eqns. of moFon, Landau-Lifshitz eqns, in imaginary Fme.

(θ, φ)

δS = 0

δ2S = Z β/2

−β/2

dτ h −2isδθδ ˙ φ sin θ +

• aδ2θ + 2bδθδφ + cδ2φ

i

where At the first order saddle point of the energy surface, a ≡ ∂2U(θ†, φ†) ∂θ2

• θ†, φ† , c ≡ ∂2U(θ†, φ†)

∂φ2

• θ†, φ† , b ≡ ∂2U(θ†, φ†)

∂θ∂φ

• θ†, φ† .

θ†, φ†

The onset temperature for tunneling is the temperature at which a second eigenvalue becomes negaFve. This gives

Tc = √ b2 − ac 2π s kB sin θ†

Can be applied to any system, described for example by DFT or Alexander-Anderson SCF model.

Tunneling of magne'c moments in molecular magnets

• S. M. J. Aubin et. al., J. Am. Chem. Soc. 120, 4991 (1998)

MagneFc group in monomer

Sergei

Methodology development:

• Find opFmal path between given endpoints (NEB method)
• Find the global minimum or maximum of a funcFon (GOUST)

Other applicaFons:

• Development of models for

geothermal reservoirs, here Laugarnes area (previously Hellisheiði area) with Manuel Plasencia

• Calculate the fastest path for a seismic wave
• Calculate the path of an electromagneFc signal in the ionosphere

Similar methodology can be used in many contexts

More approximate descripFons to simulate larger/more complex systems

Non-collinear Alexander Anderson model for magneFc states of transiFon metal systems Analogous to Fght binding approximaFon for electronic structure.