Large scale electronic structure: linear scaling

1
Large scale electronic structure linear scaling?

• In this lecture well review traditional
  plane-wave electronic structure, grid methods,
  other basis sets, and how one can achieve
  near-linear scaling for electronic structure
• First, we go back to discussions of the density
  matrix (DM), as that quantity will play a key role

2
Martin p. 60 --gt Density Matrix

• Recall
• Helmholtz free energy
• Equilibrium DM
• In basis of hamiltonian eigenstates

3
Grand canonical ensemble

• Number of particles allowed to vary the chemical
  potential is assigned to yield correct average
  number of particles
• Averages

4
Non-interacting particles

• This is case for Hartree and KS-DFT electrons
  dont directly interact, they move in an
  effective potential (which is created by
  distribution of all the electrons)

5
Finite T non-interacting

• Expectation values are sums over many-particle
  states each with occupation numbers
  for independent particle states with energy
• For electrons occupancies are 0 or 1 and we have
• Then

6
Non-interacting (contd)

• Average energy at some T
• Single-body DM
• In position representation

7
Deriving above equations
Exercise Derive the Fermi-Dirac distribution for
non-interacting particles from the general
definition of the DM using the
fact that the sum over many-body states can be
reduced to a sum over all possible occupation
numbers for each of the independent particle
states, subject to the conditions that each
occupation number can be either 0 or 1 and the
sum of all the occupanices gives the total number
of electrons. (Martin exercise 3.7, p. 71)
Exercise Following Exercise above, show
that simplifies to
for any operator in the independent particle
approximation (Martin ex. 3.8 p. 71)
Exercise show that in the limit of T--gt0, the
above sums truncate sharply at the highest
occupied state (hint what does the FD
distribution look like at 0 K?
8
Plane waves p. 236 Martin

• Plane waves are completely nonlocal, but they
  also have desirable properties, esp that they
  allow for use of FFTs
• Expand wavefunctions in PW basis

9
PWs contd

• Reciprocal lattice vectors
• Then
• Note the KE is diagonal
• Need pseudopotentials, need huge number of PWs to
  represent core states

10
Grids

• Martin, p. 248. Beck, Rev. Mod. Phys. 72, 1041
  (2000).
• Finite differences and finite elements
• Finite differences Taylor expansion of function
  about a central point (not a basis set expansion,
  so variational theorem not satisfied)
• Finite elements localized basis set
  representation

11
Finite differences

• Taylor expansion
• Add and rearrange to get

Exercise show this
12
FDs contd

• The above approximation is accurate to 2nd order
  in the grid spacing
• We can go to higher orders with higher order
  Taylor expansions
• Note the representation is near-local in space
  only neighboring grid points are needed to
  calculate the action of the kinetic energy
  operator on the wavefunction.
• Localization important for linear scaling

Exercise see the Physica Status Solidi B issue
-- vol. 243, 5, April 2006 for extensive
discussion of real-space methods
13
Intro to solving PDEs with FDs

• Functional derivatives, an example, the Poisson
  equation (1D here)
• Where does this come from?
• A functional maps a function onto a number, as we
  see above

14
Functionals contd

• S above is an action, and it is a functional.
  We want to minimize that functional wrt
  variations in the potential until we reach the
  minimum. Minus the functional derivative is a
  force which drives the action downhill.
• To see how the functional derivative is obtained,
  we discretize the equations with FDs. Lets then
  write out the action in a FD representation.

15
FD action

• FD action in 1D
• Functional derivative as a limit
• What we get
• In the limit of grid spacing going to zero
• We see that when this force goes to zero we
  have a solution to the Poisson equation

Exercise confirm the above equations
16
Iterative scheme to drive the action downhill

• Write out a steepest descent equation for the
  dynamics of the potential
• If this equation is iterated, the action will
  gradually be pushed downward until it reaches a
  minimum. What does the update equation then look
  like?

Note this may seem esoteric, but it turns out to
be highly practical -- well derive the typical
update equations below
17
Weighted Jacobi updates

• We discretize this equation in time and space.
  The time step size will be called tau, and the
  time step number indicated by t. The spatial
  grid points are labeled with i.
• Now introduce the parameter
• We get

18
Iterative updates contd

• It turns out that the limit of stability of this
  weighted Jacobi scheme is Well
  show that below.
• One minor change we can make is to say that, once
  we have already updated the i-1 point, we can use
  the new value

This is called Successive Over-relaxation (SOR)
19
Iterative updates contd

• If we make omega1, then we have
• This is called Gauss-Seidel iteration, which
  actually is useful in multigrid methods
• We can view all of these update equations as a
  matrix times a vector

Exercise be able to derive all of these update
schemes
20
Spectral props of update eqns

• Lets go back and analyze the spectral properties
  of the weighted Jacobi update equation
• The update matrix is highly banded, having terms
  only on and next to the diagonal
• For the present argument (stability), the charge
  density does not play a role, so well assume no
  charge for now (which means were solving the
  Laplace eqn.)

21
Spectral props of update eqns

• Let the update matrix act on a plane wave
• One can always diagonalize a symmetric matrix, so
  application of the update matrix involves
  repeated multiplication by the eigenvalues. Thus
  if the eigenvalues are greater than 1, the
  process will explode.

22
Spectral props of update eqns

• The eigenvalues
• For the magnitude of the eigenvalues to be less
  than or equal to 1, need
• Also, notice that for small k (long wavelength)
• Thus for long wavelengths the eigenvalues
  approach 1

See Beck, Rev Mod Phys v72, p1041 (2000)
Exercise derive the above expression for the
eigenvalues and discuss how this explains why a
solver on a large domain may slow down in the
solution process
23
Action for Schrodinger eq.

• Solving the Schrodinger equation is quite
  similar, with some differences
• Differences since the wavefunctions can in
  general be complex, when we take the functional
  derivatives, the derivative wrt is
  independent of the derivative wrt
• Also there is a contraint for wave function
  normalization (or orthonormality more generally).
  See my review for more discussion.

24
Multigrid methods

• Vcycle

Exercise discuss how going to multiple levels
might eliminate the slowing down in the solution
process
25
MG methods contd

• Multigrid was developed to overcome the slowing
  down that we found above from the spectral
  analysis of the weighted Jacobi update matrix.
• Developed by Brandt, Hackbusch, and others in
  1970s.
• Now applied to many problems in computational
  science
• Very powerful numerical methods
• See my RMP review and Martins book

Exercise why do you think these grid methods
might be good for linear scaling applications?
26
Linear scaling

• Martin p. 450
• Beck, Rev. Mod. Phys. 72, 1041 (2000) --gtreal
  space
• Goedecker, Rev. Mod. Phys. 71, 1085 (1999). --gt
  linear scaling algorithms
• Physical Status Solidi issue vol. 243, 5, April
  2006 for other algorithms (ONETEP, MIKA,
  CONQUEST, etc)

27
Linear scaling contd

• LS methods rely somehow on localization or
  nearsightedness (Walter Kohn)
• This is a consequence of the decay of the
  one-particle, off-diagonal density matrix in real
  space (discussed above)

28
LS contd

• If there is a band gap (HOMO/LUMO gap in chemical
  terms), then the DM decays exponentially
• Also, even a metal at finite T yields exponential
  decay
• Thus at some distance we can truncate the DM in
  space without substantial loss of information
• We are interested in integrated quantities like
  the electron density and the total energy

29
LS contd

• Represent the orbitals as orthogonal or
  non-orthogonal Wannier functions
• Rely on sparsity of Hamiltonian and overlap
  matrices
• Build Hamiltonian
• Solve the equations
• A wide range of methods has been developed for
  the above. See Martins excellent review in Ch.
  23 of his book. Some of the new codes are
  ONETEP, SIESTA, and CONQUEST

30
Basis sets for LS

• SIESTA numerical, localized atomic orbitals
• ONETEP localized functions built from periodic
  plane waves
• CONQUEST finite elements
• MIKA and other MG codes finite differences (not
  a basis)

Exercise again see Physical Status Solidi issue
vol. 243, 5, April 2006 for algorithm
discussions
31
Electrostatics in PBC

• In periodic boundaries we need a method to handle
  the long-ranged charge-charge interactions. This
  is done with the Ewald method.
• The Ewald method has been converted to linear
  scaling in recent particle mesh Ewald
  formulations T. Darden, D. York, and L. Pederson,
  J. Chem. Phys. 98 (1993) 10089.
• See our book (Ch. 5) for discussion of Ewalds
  method, and its importance for computation of
  ion solvation free energies.