Systematic convergence for realistic projects

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
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 and minimization
  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
Radial flexibilization Add more than one radial
function with the same angular momentum
Multiple-? Diffuse orbitals, ghosts, ...
Angular flexibilization Add shells of different
atomic symmetry Polarization
Golden Rule of Quantum Chemistry a basis should
be doubled before its 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 A 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
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ía and 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)
ltfi(r-Ri) Vlocal(r-Rk) fj(r-Rj)gt -IMPORTANT
this 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
?x
Grid points

• 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 incremental procedure 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
a 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-4 or 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.Tol10-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.Tol10-3
• Rescaling to match DZ bulk lattice parameter
• iii) DZ basis set, 8x4 sampling, full relaxation,
  DM.Tol10-4
• iv) DZP basis set, 8x4 sampling, full relaxation,
  DM.Tol10-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 be carefully 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.