Real-space multigrid methods for DFT and TDDFT:

1
Real-space multigrid methods for DFT and TDDFT
Tuomas Torsti CSC The finnish IT center for
Science Laboratory of Physics, Helsinki
University of Technology
http//www.csc.fi/physics/mika
2
Acknowledgements

• For Funding
• CSC The finnish IT center for Science
• COMP, Helsinki University of Technology
• For advice
• Martti Puska (COMP)
• Risto Nieminen (COMP)
• Janne Ignatius (CSC)
• For collaboration using MIKA/cyl2
• Bo Hellsing (Chalmers)
• Vanja Lindberg (Växjö, Chalmers)
• Nerea Zabala (San Sebastian)
• Eduardo Ogando (Bilbao)
• Paula Havu (COMP)
• Tero Hakala (COMP)

• For development of RQMG
• Mika Heiskanen (then COMP)
• For collaboration in development of MIKA/rspace
• Sampsa Riikonen (now San Sebastian)
• Ville Lehtola (COMP)
• Kaarle Ritvanen (COMP)
• For work done with MIKA/RS2Dot
• Henri Saarikoski
• Esa Räsänen
• For response iterations
• Eckhardt Krotscheck (Linz)
• Michael Aichinger (Linz)
• For work done with MIKA/doppler
• Ilja Makkonen

3
Motivation for using real-space grids

• With uniform grids the control of the basis set
  is simple Only one parameter (the grid spacing
  h)
• Flexible choice of boundary conditions cluster,
  wire, surface, bulk.
• cluster
• wire
• surface
• bulk
• ...
• Parallelization using domain decomposition
• It is possible to use nonuniform grids to refine
  the mesh close to atomic nuclei or hard
  pseudopotential, and/or to push the vacuum
  boundary far away in cluster calculations
• adaptive grids
• composite grids
• finite elements
• Multigrid techniques can be used to obtain
  optimal scaling for PDEs
• Natural framework for Order-N (localized orbitals
  required)

4
Multigrid methods
A. Brandt. Math. Comput. 31, 333 (1977)., T. L.
Beck. Rev. Mod. Phys. 72, 1041 (2000). W. L.
Briggs et al., A Multigrid Tutorial, Second
Edition. (SIAM 2000).
As a simple example, take the Poisson equation
Simple relaxation schemes (e.g. the Gauss-Seidel
method) efficiently remove the short wavelength
components of the residual
(they are good smoothers), while critical slowing
down occurs for the long wavelength components.
Solution treat long wavelength components of V
on a coarse grid
The idea can be applied recursively (V-cycle).
Linear scaling with problem size can be acchieved
with the full-multigrid method.
5
Classification of MG-methods for the eigenproblem

• Steepest descent (or CG or RMM-DIIS) with
  MG-preconditioning
  e.g.
  Bernholc et al., Phys. Rev. B 54 14362 (1996)
• Full approximation storage
  A. Brandt et al.
  SIAM J. Sci. Comput. 4, 244 (1983)
  J. Wang and T. L. Beck , J. Chem. Phys. 112,
  9223 (2000)
• Rayleigh Quotient Multigrid method (RQMG)
  J. Mandel and S. F. Cormick,
  J. Comput. Phys. 80, 442 (1989). M. Heiskanen et
  al., Phys. Rev. B 63, 245106, (2001).

6
Rayleigh quotient multigrid method
J. Mandel and S. F. Cormick, J. Comput. Phys. 80,
442 (1989). M. Heiskanen et al., Phys. Rev. B 63,
245106, (2001).

• Discretized Schrödinger equation

• With search vector d vary a to minimize the
  Rayleigh quotient

• Coordinate relaxation choose a coordinate
  vector de.
• RQMG method on coarse grids minimize the
  fine grid RQ with

• The fine grid Rayleigh quotient can be
  evaluated entirely on the coarse grid

• If eigenpairs other than the lowest one are
  required, add a penalty functional to take care
  of the orthogonality requirement

7
Rayleigh quotient multigrid method (continued)

• Galerkin conditions should hold

• In the original implementation, approximated by
  discretization coarse grid approximation (DCA).
  In MIKA/rspace 1.0 also the Galerkin case
  implemented

• Can we get rid of the penalty functional by
  minimizing the residual norm instead of the
  Rayleigh Quotient (In analogy with the familiar
  RMM-DIIS method) ?

8
Response iteration method full response
J. Auer and E. Krotscheck, Comp. Phys. Comm. 151
(2003), 265-271

• Newton-Raphson method for the equation

• Full response equation (needs unoccupied states)
  (solve with CG or GMRES)

where
9
Response iteration method collective
approximation
J. Auer and E. Krotscheck, Comp. Phys. Comm. 151
(2003), 265-271

• requires only occupied states
• implemented in MIKA/cyl2 and MIKA/RS2Dot

10
MIKA/rspace 1.0

• Parallelized over k-points and real-space domains
• Periodic and cluster boundary-conditions
  implemented
• Norm-concerving nonlocal pseudopotentials of the
  Kleynman-Bylander form (usually Troullier-Martins
  pseudopotentials are used), double-grid method
• Hellman-Feynman Forces
• Structural optimization with the BFGS-method (two
  implementations)
• Mixing schemes
• Pulay
• Broyden
• GR-Pulay (D. R Bowler and M. J. Gillan. Chem.
  Phys. Lett. 325, 473 (2000) ),
• screened Coulomb interaction (M. Manninen et
  al., Phys. Rev. B 12, 4012 (1975). )
• Pulay-Kerker (Note rough Fourier components
  obtained using a MG-technique)
• Pulay-Kerker with metric (motivated by Kresse and
  Furthmuller, PRB 54, 11169).

11
MIKA/rspace (future)

• Mixed boundary conditions for surface
  computations
• Spin-dependent version of the code
• Alternative MG-solver (e.g. RMM-DIIS with
  MG-preconditioning)
• PBE (Perdew, Burke, Ernzerhof) GGA correction
  already implemented, and will be included in the
  next release
• Response iterations (already implemented in other
  MIKA-codes, 3D subroutines from prof. Krotscheck
  available)
• Build an interface to Octopus for time-dependent
  calculations

12
Double grid method for nonlocal pseudopotentials
T. Ono and K. Hirose, PRL 82, 5016 (1999)

• Replaces the fourier filtering of
  pseudopotentials of Briggs et al.
• The idea should be understood as a general
  scheme to transfer a function from a fine grid to
  a coarse grid, and is in fact equivalent to the
  MG restriction operation.
• Should be applied also to the local part, and
  compensating gaussian charges (all functions that
  are transferred from a radial grid to the
  computational grid)
• Thanks to J. J. Mortensen (CAMP, DTU) who
  implemented this in grid-based PAW.
  

13
All-electron finite-element calculations with
ELMER

• These are outside the scope of the MIKA-project,
  but demonstrated the capabilities of CSCs ELMER
  package.

14
Vortex clusters in quantum dots
Left SDFT density of 24-electron QD at 5T
showing 14 vortice Right CSDFT density and
currents at the edge of the QD.

• Saarikoski et al. Phys. Rev. Lett (2004)
  (cond-mat/0402514)
• Exact diagonalization and DFT (both CSDFT and
  SDFT) give corresponding results limitations
  and differences of the methods discussed.
• Finding the vortex solution in DFT requires high
  numerical accuracy. Our real-space implementation
  is superior to existing plane-wave schemes in
  describing the vanishing density at the vortex
  core

15
Conductance oscillations in metallic nanocontacts
P. Havu et al., Phys. Rev. B, 66, 075401 (2002).

• We model a chain of N Na atoms between two
  conical stabilized jellium leads
• Since only one channel contributes to the
  conductance, and because of the mirror symmetry,
  the Friedel sum rule can be applied for the
  conductance

• We observe the even-odd behaviour of the
  conductance as the function of N
• In addition, the important role of the leads is
  manifested as an additional oscillation as a
  function of the cone opening angle

16
Ultimate jellium model for breaking nanowires
E. Ogando et al., Phys. Rev. B 67, 075417 (2003).

• Ultimate jellium is a locally neutral model, the
  compensating background charge density equals
  the electron density at every point.
• The shape of the system in the central part is
  free to vary to minimize the total energy.
• The shape of the leads is frozen to the uniform
  wire solution.
• In the beginning of the elongation, classical
  catenoid shape is observed
• Quantum mechanical shell structure in
  cylindrical symmetry -gt cylinders with magic
  radii.
• Quantum mechanical shell structure in sperical
  symmetry -gt Cluster derived structures (CDS)
• Oscillation of elongation force

17
Model study of adsorbed metallic quantum dots Na
on Cu(111)
T. Torsti et al., Phys. Rev. B 66, 235420 (2002)

• Roughly hexagonal islands are observed to form on
  the second monolayer of Na grown on Cu(111)
• Bandgap at Fermi level in Cu for electrons moving
  in the (111) direction gt quantum well states
• We developed a two-jellium model to fit the
  bottoms of two surface state bands
• The infinite monolayer is described with as a
  hexagonal lattice of circles, the k-space is
  sampled with two points.
• In the largest system studied, 2400 states are
  solved the code is parallelized over the 652
  different values of (m,k). This is also a
  demanding test for the charge density (or
  potential) mixing.
• The local density of states is calculated at a
  realistic STM-tip distance (15 a.u.) above the
  surface and compared with measured differential
  conductance

18
Quantum corrals (Tero Hakala, M.Sc. project)

• We use a pseudopotential (E. Ogando et al.
  submitted to PRB, cond-mat/0310533) for the
  Cu(111) surface
• A ring of 45 Pb atoms on both sides of a Cu(111)
  slab with 5 atomic layers and radius 60 bohr a
  localized surface state observed within the
  corral
• The total system size was 3272 electrons and
  required about 2000 SCF-iterations to converge
  (about 1 day with 8 processor in the IBM server
  cluster of CSC).

19
Quantum corrals (continued)

• Charge transfer in a corral with 8 Pb-atoms on
  both sides of a Cu(111)-slab with15 atomic
  layers. This transfer is due to the equilibration
  of chemical potentials between Pb and Cu.
• It has been observed also in 1D-calculations of
  Pb on top of Cu(111) by Ogando et al.

20
Partial list of publications related to MIKA
Numerical methods M. Heiskanen, T. Torsti, M.J.
Puska, and R.M. Nieminen, Multigrid method for
electronic structure calculations, Phys. Rev. B
63, 245106 (2001). T. Torsti, M. Heiskanen, M.J.
Puska, and R.M. Nieminen, MIKA a multigrid-based
program package for electronic structure
calculations, Int. J. Quantum Chem. 91, 171-176
(2003). T. Torsti, Real-Space Electronic
Structure Calculations for Nanoscale Systems, CSC
Research Reports R01/03 (Ph. D.
-thesis). Applications to axially symmetric
model systems P. Havu, T. Torsti, M.J. Puska,
and R.M. Nieminen, Conductance oscillations in
metallic nanocontacts, Phys. Rev. B 66, 075401
(2002). T. Torsti, V. Lindberg, M. J. Puska, and
B. Hellsing Model study of adsorbed metallic
quantum dots Na on Cu(111) Physical Review B 66,
235420 . E. Ogando, T. Torsti, N. Zabala, and M.
J. Puska, Electronic resonance states in
metallic nanowires ... simulated with the
ultimate jellium model, Phys. Rev. B. 67, 075417
T. Torsti, Real-Space Electronic Structure
Calculations for Nanoscale Systems, CSC Research
Reports R01/03 (Ph. D. -thesis) Applications to
quantum dots in 2DEG Saarikoski, H. , Harju, A. ,
Puska, M. J., Nieminen, R. M., Vortex Clusters in
Quantum Dots, Submitted to Physical Review
Letters on 19.2.2004 Harju, A., Räsänen, E.,
Saarikoski, H., Puska, M.J., Nieminen, R.M., and
Niemelä, K., Broken symmetry in
density-functional theory Analysis and cure,
Submitted to Physical Review B on 3.2.2004
Räsänen, E., Harju, A., Puska, M. J., and
Nieminen, R. M., Rectangular quantum dots in high
magnetic fields, Submitted to Physical Review B
on 27.11.2003. Räsänen, E., Puska, M.J., and
Nieminen, R.M., Maximum-density-droplet formation
in hard-wall quantum dots, Submitted to Physica E
on 9.6.2003. Räsänen, E., Saarikoski, H.,
Stavrou, V. N., Harju, A., Puska, M.J., and
Nieminen, R.M., Electronic structure of
rectangular quantum dots, Physical Review B 67,
235307 (2003) . Saarikoski, H., Räsänen,
E.,Siljamäki, S., Harju, A., Puska, M.J.,
Nieminen, R.M., Testing of two-dimensional local
approximations in the current-spin and
spin-density-functional theories, Physical Review
B 67, 205327 (2003) . Räsänen, E., Saarikoski,
H., Puska, M. J., and Nieminen, R. M., Wigner
molecules in polygonal quantum dots A
density-functional study, Physical Review B 67 ,
035326 (2003) . Saarikoski, H., Räsänen, E.,
Siljamäki S., Harju A., Puska, M.J., and
Nieminen, R.M., Electronic properties of model
quantum-dot structures in zero and finite
magnetic fields, European Physical Journal B 26 ,
241-252 (2002) . Applications of the RQMG method
to one-dimensional problems Engström, K.,
Kinaret, J., Puska, M.J., and Saarikoski, H.,
Influence of Electron-Electron Interactions on
Supercurrent in SNS structures, Low Temperature
Physics 29, 546 (2003). Ogando,E. Zabala,N.,
Chulkov,E.V., Puska,M.J., Quantum size effects in
Pb islands on Cu(111) Electronic-structure
calculations, Submitted to Phys. Rev. B on
22.10.2003
21
Summary

• MIKA (Multigrid Instead of the K-spAce) is a
  collection of programs that solve the Kohn-Sham
  equations of DFT in one, two and three
  dimensional cartesian coordinate systems or in
  axial symmetry
• The core numerical method is the Rayleigh
  quotient multigrid method for the eigenproblem
• No TDDFT yet, but this has a high priority as a
  future development.
• MIKA / rspace 1.0 was released on 2.9.2004. Along
  with the other codes, it is licensed with the
  GPL, and available from http//www.csc.fi/physics/
  mika