Systematic convergence for realistic projects From very fast to very accurate

1. Systematic convergence for realistic projects From very fast to very accurate Eduardo Anglada D. Sánchez-Portal Thanks to: J. Junquera, E. Artacho, S. Riikonen, S. García, José M. Soler, A. García and P. Ordejón

2. Basic strategy • Steps of a typical research project: • Exploratory-feasibility tests • Convergence tests • Converged calculations A fully converged calculation is impossible without convergence tests

3. Goal: Approx. parameter independence: • Monitor: • Convergence • CPU time & memory Convergence tests • Choose relevant magnitude(s) A of the problem (e.g. an energy barrier or a magnetic moment) • Choose set of qualitative and quantitative parameters xi (e.g. xc functional, number of k-points, etc)

4. Practical hints • Ask your objective: find the truth or publish a paper? • Do not try a converged calculation from the start • Start with minimum values of all xi • Do not assume convergence for any xi • Choose a simpler reference system for some tests • Take advantage of error cancellations (with care!) • Refrain from stopping tests when results are “good”

5. Almost complete parameter list for ab initio simulations • XC functional: LDA, GGAs • Pseudopotential • Method of generation • Number of valence states • Number of angular momenta • Core matching radii • Nonlinear core corrections • Basis set • Number of functions • LCAO: • Highest angular momentum • Number of zetas • Range • Shape (Optimized!) • Real space mesh cutoff (Vxc) •Spin polarization • Number of k-points •SCF convergence tolerance • Supercell size (solid & vacuum) •Geometry relaxation tolerance • Electronic temperature •O(N) Rc andminimization tolerance

6. Parameter interactions 2A/xixj0 • Number of k-points: • Supercell size • Geometry • Electronic temperature • Spin polarization • Harris vs SCF • Mesh cutoff: • Pseudopotential • Nonlinear core corrections • Basis set • Functional: GGA

7. Basis Sets

8. Convergence of the basis set size • Single- (minimal or SZ) • One single radial function per angular • momentum shell occupied in the free –atom Ways to improve the quality • Radialflexibilization: • Add more than one radial function with the same angular momentum: Multiple- • + Diffuse orbitals, ghosts, ... • Angular flexibilization: • Add shells ofdifferent atomic symmetryPolarization Golden Rule of Quantum Chemistry: a basis should be doubled before it’s polarized

9. Convergence of the basis set size: bulk Si Cohesion curves PW and NAO convergence

10. QUALITY OF THE RESULTS WITH A SZ BASIS? • not so bad as you might expect: • energetics changes considerably; however, energy differences might be reasonable enough • charge transfer and other basic chemistry is usually OK • (at least in simple systems) • if the geometry is set to the experiment, we typically have a nice band structure for occupied and lowest unoccupied bands • When can SZ basis sets be used: • long molecular dynamics simulations (once we have make ourselves sure that energetics is reasonable) • exploring very large systems and/or systems with many degrees of freedom (complicated energy landscape).

11. Example of SZ calculation: very long molecular dynamics “Atomic Layering at the liquid silicon surface: a first principles simulation” G. Fabricius, E. Artacho, D. Sánchez-Portal, P. Ordejón, D. Drabold, J. Soler Phys. Rev. B 60, R16283 (1999) Only the gamma k point was used in the simulations, since previous work found cell-size effects to be small. For the present calculation we used a mini- mal basis set of four orbitals and 3p for each Si atom, with a cutoff radius of 2.65 Ang. We have extensively checked the basis with static calculations of different crystalline Si phases and solid surfaces, and MD simulations of the bulk liquid. The energy differences between solid phases are described within 0.1 eV of other ab initio calculations. The diamond structure has the lowest energy, with a lattice parameter of 5.46 Å Adatom- and dimer-based and surface reconstructions found in other ab initio calculations are well reproduced, with geometries and relative energies changing less than 0.1 Ang and 0.15 eV when moving from a gamma-point calculation with a minimal basis set, to a converged k sampling with double- and polarization orbitals. The structural, electronic, and dynamical properties of l - Si are in good agreement with other ab initio calculations at the same density and temperature. The calculated diffusion constant is somewhat smaller than that obtained with a double- or polarized basis which is in agreement with other ab initio simulations. We interpret that the minimal basis overesti- mates the energies of the saddle point configurations occurring during diffusion, but we consider that this is not critical for the present application

12. Convergence with respect to the orbitals range • Remarks: • May require long radii (specially in molecules/surfaces) • Affects total energies, but energy differences converge better • For incomplete bases (SZ), increasing radii may not produce better properties. bulk Si equal s, p orbitals radii J. Soler et al, J. Phys: Condens. Matter, 14, 2745 (2002)

13. A single parameter for all cutoff radii E. Artacho et al. Phys. Stat. Solidi (b) 215, 809 (1999) Convergence vs Energy shift of Bond lengths Bond energies Range: Energy shift Reasonable values for practical calculations: ΔEPAO ~50-200 meV

14. Procedure to converge the basis Difference in energies involved in your problem? • SZ: (Energy shift) • Semiquantitative results and general trends • DZP: automatically generated(Split Valence and Peturbative polarization) • High quality for most of the systems. • Good valence: well converged results computational cost • ‘Standard’ • Variational optimization

15. Gold (111) surface: basis set completeness Special thanks to S. Garcíaand P. Ordejón

16. Gold (111): Wavefunctions -1.22 eV -0.32 eV -0.29 eV +3.38 eV

17. Gold (111): Bulk Siesta :DZP optimized Abinit Siesta: DZP optim. + diffuse orbitals

18. Au (111): Surface states Siesta :DZP Abinit Siesta: DZP + diffuse orbitals

19. Other Parameters

20. Real-space grid: Mesh cut-off Used to compute ρ(r) in order to calculate: -XC potential (non linear function of ρ(r) ) -Solve Poisson equation to get Hartree potential -Calculate three center integrals (difficult to tabulate and store) <φi(r-Ri) | Vlocal(r-Rk) | φj(r-Rj)> -IMPORTANTthis grid is NOT part of the basis set… is an AUXILIARY integration grid and, therefore, convergence of energy is not necessarily variational respect to its fineness. -Mesh cut-off: highest energy of PW that can be represented with such grid.

21. Convergence of energy with the grid: mesh cutoff • Important tips: • Convergence is rarely achieved for less than 100 Ry. • Values between 150 and 200 Ry provide good results in most cases • GGA and pseudo-core require larger values than other systems • To obtain very fine results use GridCellSampling • Filtering of orbitals and potentials coming soon (E. Anglada)

22. Egg box effect atom Energy of an isolated atom moving along the grid E(x) ΔE Grid points Δx • We know that ΔE goes to zero as Δx goes to zero, but • what about the ratio ΔE/Δx?: • Tipically covergence of forces is somewhat slowler than for the total energy • This has to be taken into account for very precise relaxations and phonon calculations. • Also important and related: tolerance in forces (for relaxations, etc) should not be smaller than tipical errors in the evaluation of forces.

23. k-sampling • -Only time reversal symmetry used in SIESTA (k=-k) • -Convergence in SIESTA not different from other codes: • Metals require a lot of k-point for perfect convergence • (explore the Diag.ParallelOverK parallel option) • Insulators require fewer k-points • -Gamma-only calculations should be reserved to really large simulation cells • -As usual, an incrementalprocedure might be the most intelligent approach: • Density matrix and geometry calculated with a “reasonable” • number of k-points should be close to the converged answer. • Might provide an excellent input for more refined calculations

24. K sampling: Al kgrid_cutoff (Moreno and Soler, PRB 45, 13891 (1992)): Automatic generation of integration grid kgrid_Monkhorst_Pack (Monkhorst and Pack, PRB 13, 5188 (1997)): Grid defined by hand

25. - | | Convergence of the density matrix DM.MixingWeight: αis not easy to guess, has to be small (0.1-0.3) for insulator and semiconductors, tipically much smaller for metals DM.NumberPulay (DM.NumberBroyden) : N such that is minimum N between 3 and 7 usually gives the best results

26. Convergence of the density matrix DM.Tolerance: you should stick to the default 10-4or use even smaller values …… … except in special situations: - Preliminary relaxations - Systems that resist complete convergence, but you are almost there - in particular if the Harris energy is very well converged - Warning: above 10-3 errors may be too large. - ALWAYS CHECK THAT THINGS MAKE SENSE.

27. A particular case: Determination of the Si(553)/Au structure More than 200 structures explored!! DM.Tolerance could be reduced Special thanks to: S. Riikonen and D. Sánchez-Portal

28. Si(553)/Au: Incremental approach to convergence • SZ basis set, 2x1 sampling, constraint relaxations, • slab with two silicon bulayers, DM.Tol=10-3 • only surface layer: first relax interlayer height, then relaxations with some • constraints to preserve model topology. • Selecting a subset with the most stable models • DZ basis set, 8x4 sampling, full relaxation, • slab with four silicon bilayers, DM.Tol=10-3 • Rescaling to match DZ bulk lattice parameter • iii) DZ basis set, 8x4 sampling, full relaxation, DM.Tol=10-4 • iv) DZP basis set, 8x4 sampling, full relaxation, DM.Tol=10-4 • Rescaling to match DZ bulk lattice parameter

29. Si(553)/Au:SZ energies might be a guide to select reasonable candidates, but… caution is needed!!!! SZ: 88 initial structures converge to the 41 most stable different models All converged to the same structure Already DZ recovers these as the most stable structures

30. Si(553)/Au:Finally we get quite accurate answer….

31. Conclusions • All the ab initio methods have approximations, and they must becarefully checked (convergence tests!!). • The best way to handle them is an incremental approach, start from a reasonable guess and converging the parameters that control the approximations. • In order to save computing time, use the lowest values of the parameters which provide reasonably converged results. • There are no “magic numbers”, the effect of the parameters is system dependent.