Linear Scaling ‘Order-N’ Methods in Electronic Structure Theory Richard M. Martin University of Illinois Acknowledgements: Pablo Ordejon Matthew Grumbach Daniel Sanchez-Portal David Drabold Uwe Stephan Satoshi Itoh Thanks to: Jose Soler, Emilio Artacho, Giulia Galli, ...

Comp. Mat. Science School 2001 1

Linear Scaling ‘Order-N’ Methods and Car-Parrinello Simulations • Fundamental Issues of locality in quantum mechanics • Paradigm for view of electronic properties • Practical Algorithms • Results Comp. Mat. Science School 2001 2

Locality in Quantum Mechanics • V. Heine (Sol. St. Phys. Vol. 35, 1980) “Throwing out k-space” Based on ideas of Friedel (1954) , . . .

• Many properties of electrons in any region are independent of distant regions • Walter Kohn “Nearsightness” Comp. Mat. Science School 2001 3

Locality in Quantum Mechanics • Which properties of electrons are independent of distant regions?

• Total integrated quantities Density, Forces on atoms, . . .

• Coulomb Forces are long range but they can be handled in O(N) fashion just as in classical systems Comp. Mat. Science School 2001 4

Non-Locality in Quantum Mechanics • Which properties of electrons are non-local ?

• Individual Eigenstates in crystals • Sharp features of the Fermi surface at low T • Electrical Conductivity at T=0 Metals vs insulators : distinguished by delocalization of eigenstates at the Fermi energy (metals) vs localization of the entire many-electron system (insulators) • Approach in the Order-N methods: Identify localized and delocalized aspects Comp. Mat. Science School 2001 5

Density Matrix I • Key property that describes the range of the non locality is the density matrix r (r,r’) • In an insulator r (r,r’) • In a metal r (r,r’) is exponentially localized decays as a power law at T = 0, exponentially for T > 0.

(Goedecker, Ismail-Beigi) • For non-interacting Bosons or Fermions, Landau and Lifshitz show that the correlation function is uniquely related to the square of r (r,r’) g(r,r’) • Thus correlation lengths and the density matrix generally become shorter range at high T Comp. Mat. Science School 2001 6

Density Matrix II • Key property that describes the range of the non locality is the density matrix r (r,r’) • Definition: r (r,r’) = S i f i * (r) f i (r’) • Can be localized even if each f i * (r) is not!

Atom positions r’ r fixed at r =0 Comp. Mat. Science School 2001 7

Toward Working Algorithms I (My own personal view) Heine and Haydock laid the groundwork - but it was applied only to limited Hamiltonians, ….

1985 - Car-Parrinello Methods changed the picture Key quantity is the total energy E[{ f i }] which does not require eigenstates occupied states - only traces over the the { f i } can be linear combinations of eigenstates Comp. Mat. Science School 2001 8

Toward Working Algorithms II How can we use the advantages of the Car-Parrinello and the local approaches?

1992 - Galli and Parrinello pointed out the key idea to make a Car-Parrinello algorithm that takes advantage of the locality Require that the states in localized.

Note this does not require a localized basis it may be very convenient, but a localized basis is not essential to construct localized states (example: sum of plane waves can be localized) Comp. Mat. Science School 2001 9

Toward Working Algorithms III What are localized combinations of the eigenfunctions?

Extended Bloch Eigenfunctions Wannier Functions (generalized)! Wannier Functions span the same space as the eigenstates - all traces are the same Wannier Functions One localized Wannier Ftn centered on each site Comp. Mat. Science School 2001 10

Toward Working Algorithms IV Can work with either localized Wannier functions w i (r) Functions of one variable But not unique or localized density matrix r (r,r’) = S i f i * (r) f i (r’) = S i w i * (r) w i (r’) Functions of two variables - more complex But unique Comp. Mat. Science School 2001 11

Linear Scaling ‘Order-N’ Methods • Computational complexity ~ N = number of atoms (Current methods scale as N 2 or N 3 ) • Intrinsically Parallel • “Divide and Conquer” • Green’s Functions • Fermi Operator Expansion • Density matrix “purification” • Generalized Wannier Functions • Spectral “Telescoping” (Review by S. Goedecker in Rev Mod Phys) Comp. Mat. Science School 2001 12

Divide and Conquer (Yang, 1991) • Divide System into (Overlapping) Spatial Regions.

Solve each region in terms only of its neighbors.

(Terminate regions suitably) • Use standard methods for each region • Sum charge densities to get total density, Coulomb terms Comp. Mat. Science School 2001 13

Expansion of the Fermi function • Sankey, et al (1994); Goedecker, Colombo (1994); Wang et al (1995) • Explicit T nonzero • Projection into the occupied Subspace • Multiply trial function f by “ Fermi operator ”: F • Localized f f = [(H - E F )/K B T +1] -1 f leads to localized projection since the Fermi operator (density matrix) is localized • Accomplish by expanding F in power series in H operator Comp. Mat. Science School 2001 14

Density Matrix “Purification” • Li, Nunes, Vanderbilt (1993); Daw (1993) Hernandez, Gillan (1995) • Idea: A density matrix r at T=0 has eigenvalues = occupation = 0 or 1 • Suppose we have an approximate r that does not have this property • The relation r n+1 = 3 ( r n ) 2 - 2 ( r n ) 3 always produces a new matrix with eigenvalues closer to 0 or 1.

Comp. Mat. Science School 2001 15

Density Matrix “Purification” • The relation r n+1 = 3 ( r n ) 2 - 2 ( r n ) 3 always produces a new matrix with eigenvalues closer to 0 or 1.

3 x 2 - 2 x 3 1 1 x Comp. Mat. Science School 2001 instability 16

Generalized Wannier Functions • Divide System into (Overlapping) Spatial Regions.

• Require each Wannier function to be non-zero only in a given region • Solve for the functions in each region requiring each to be orthogonal to the neighboring functions • New functional invented to allow direct minimization without explicitly requiring orthogonalization • Mauri, et al.; Ordejon, et al; 1993; Stechel, et al 1994; Kim et al 1995 Comp. Mat. Science School 2001 17

Generalized Wannier Functions • Factorization of the density matrix r (r,r’) = S i w i * (r ) w i (r’) • Can chose localized Wannier functions (really linear combinations of Wannier functions) • Minimize functional: Overlap matrix E = Tr [ (2 - S) H] • Since this is a variational functional , the Car Parrinello method can be used to use one calculation as the input to the next • Mauri, et al.; Ordejon, et al; 1993; Stechel, et al 1994; Kim et al 1995 Comp. Mat. Science School 2001 18

Functional (2-S)(H - E F ) • Minimization leads to orthonormal filled orbitals focres empty orbitals to have zero amplitude • Each matrix element (S and H) contains two factors of the wavefunction - amplitude ~ x.

• For occupied states (eigenvalues below E F ) - ( 2 x 2 - x 4 ) 1 x Minimum at zero for empty states above E F 1 Comp. Mat. Science School 2001 Minimum for normalized wavefunction (x = 1) 19

Example of Our work Prediction of Shapes of Giant Fullerenes S. Itoh, P. Ordejon, D. A. Drabold and R. M. Martin, Phys Rev B 53, 2132 (1996).

See also C. Xu and G. Scuceria, Chem. Phys. Lett. 262, 219 (1996).

Comp. Mat. Science School 2001 20

Wannier Function in a-Si

U. Stephan

Comp. Mat. Science School 2001 21

Combination of O(N) Methods

Comp. Mat. Science School 2001 22

Collision of C 60 Buckyballs on Diamond Galli and Mauri, PRL 73, 3471 (1994) Comp. Mat. Science School 2001 23

Deposition of C 28 Buckyballs on Diamond •

Simulations with ~ 5000 atoms, TB Hamiltonian from Xu, et al. ( A. Canning, G.~Galli and J .Kim, Phys.Rev.Lett. 78, 4442 (1997).

Comp. Mat. Science School 2001 24

Example of DFT Simulation (not order N)

•

Daniel Sanchez-Portal (Phys. Rev. Lett. 1999)

•

Simulation of a gold nanowire pulled between two gold tips

•

Full DFT simulation

•

Explanation for very puzzling experiment! Thermal motion of the atoms makes some appear sharp, others weak in electron microscope

Comp. Mat. Science School 2001 25

Simulations of DNA with the SIESTA code • • •

Machado, Ordejon, Artacho, Sanchez-Portal, Soler (preprint) Self-Consistent Local Orbital O(N) Code Relaxation - ~15-60 min/step (~ 1 day with diagonalization) Iso-density surfaces

Comp. Mat. Science School 2001 26

HOMO and LUMO in DNA (SIESTA code) •

Eigenstates found by N 3 method after relaxation

•

Could be O(N) for each state

Comp. Mat. Science School 2001 27

O(N) Simulation of Magnets at T > 0

• • • • Collaboration at ORNL, Ames, Brookhaven Snapshot of magnetic order in a finite temperature simulation of paramagnetic Fe. These calculations represent significant progress towards the goal of full implementation of a first principles theory of the finite temperature and non equilibrium properties of magnetic materials. Record setting performance for large unit cell models (up to 1024-atoms) led to the award of the 1998 Gordon Bell prize.

The calculations that were the basis for the award were performed using the locally self-consistent multiple scattering method, which is an O(N) Density Functional method •

Web Site:

http://oldpc.ms.ornl.gov/~gms/MShome.html

Comp. Mat. Science School 2001 28

FUTURE! ---- Biological Systems

• Examples of oriented pigment molecules that today are being simulated by empirical potentials Comp. Mat. Science School 2001 29

Conclusions

• It is possible to treat many thousands of atoms in a full simulation - on a workstation with approximate methods - intrinsically parallel for a supercomputer • Why treat many thousands of atoms?

• Large scale structures in materials - defects, boundaries, ….

• Biological molecules • The ideas are also relevant to understanding even small systems Comp. Mat. Science School 2001 30