OrderN

Development of linear scaling electronic structure methods
Jinlong Yang Hefei National Laboratory for Physical Sciences at Microscale University of Science and Technology of China 2009.11.26
Collaborators
Dr. H. J. Xiang Dr. Z. Y. Li Prof. W. Z. Liang
1

Electronic structure calculations
DFT KS equation
Wavefunctions are orthogonal to each other, resulting in O(N3) scaling in traditional methods
Challenges
Many interesting systems:
Nanostructures Soft materials Quantum dots Biological systems Surfaces
First-principles modeling
Traditional DFT is successful, but O(N3) Need improved scaling
2

Introduction for O(N) methods
CPU load
Early 90’s ~N
3
~N
~ 100
N (# atoms)
Basis of O(N) methods
Kohn's “Nearsightedness” principles:
Basic strategies for O(N) scaling density matrix (DM) based localized orbitals based recursion method
W. Kohn, Phys. Rev. Lett. 76, 3168 (1996)
3

Density Matrix in DFT
Definition for DM: DFT energy from DM:
Total energy from DM
Properties of Density Matrix
charge conservation
idempotency commutation
Insulators at 0 K
4

Locality of Density Matrix
insulator:
metal:
DM based O(N) methods
Divide-and-conquer W. Yang, Phys. Rev. Lett. 66, 1438 (1991) Fermi operator expansion Goedecker et al., Phys. Rev. Lett. 73, 122(1994) DM minimization Li et al., Phys. Rev. B 47, 10891 (1993) M. S. Daw, Phys. Rev. B 47, 10895 (1993) DM purification PM: Palser et al., Phys. Rev. B 58, 12704 (1998) TC2: A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002) TRS4: Niklasson et al., J. Chem. Phys. 118, 8611 (2003)
S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999)
5

Divide and Conquer
↓
Divide
H
H
ρ
=
↓
ρ
+
↓ Approximate
ρ
Trace-Correcting Purification
ε0 and εN are the lowest and highest eigenvalues of H, respectively.
A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002)
6

Our recent works
Spin unrestricted linear scaling electronic structure theory O(N) calculation of maximally localized Wannier function with atomic basis set Linear scaling calculation of band edge states and doped semiconductors Linear-scaling treatment of electric field perturbation in solids Linear scaling phonon calculation using density matrix perturbation theory
Part I Spin unrestricted linear scaling electronic structure theory
7

Motivation
All previous O(N) methods were applied to spinrestricted systems Development of spintronics
SMM Mn84: 1032 atoms
Tasiopoulos et al., Angew. Chem. Int. Ed. 43, 2117 (2004)
Spin unrestricted O(N) method
General density matrix:
When the z component of total spin is definite,
R. McWeeny, Rev. Mod. Phys. 32, 335 (1960).
8

Predetermined spin multiplicity
Can use the same O(N) method for each spin component PSUTC2 algorithm:
similar formula for β
A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002)
Without predetermined spin multiplicity
We define new operators:
The goal is finding ρ, which can be achieved since H can be seen as a spin restricted Hamilton. matrix problem with dimension 2Nb can be reduced to two matrix problems with dimension Nb due to the block form Satisfy the ’Aufbau’ principle
9

Spin unrestricted TC2 method
Normalize Hα and Hβ using the same scaling factors Update ρα,n and ρβ,n in the same way depending on the sum of the traces
Implementation
SIESTA Kohn-Sham self-consistent density functional method in the local density (LDA/LSD) or generalized gradient (GGA) approximations Norm-conserving pseudopotentials in its fully nonlocal KB form numerical atomic orbitals: SZ, DZ, DZP Hamilton matrix H and overlap matrix S are sparse Geometry relaxation, fixed or variable cell Constant-temperature molecular dynamics Variable cell dynamics (Parrinello-Rahman)
10

PSUTC2 and SUTC2 in SIESTA
Multiatom BCSR sparse matrix Filter the sparse matrix Cholesky transform=>orthogonal basis
Blocked approximate inverse algorithm for Z
M. Benzi, R. Kouhia, and M. Tuma, Comput. Methods Appl. Mech. Eng. 190, 6533 (2001)
Tests
O2 molecule: By specifying the spin triplet state, using the PSUTC2 method,we can get the same energy for the triplet state as that from the diagonalization calculation. SUTC2 can also give the triplet ground state with the same energy without given the electronic occupation.
11

Performance
Applications
Increased interest in magnetism in metal-free systems A recent DFT calculation indicates that carbon doping induces spontaneous magnetization in BN nanotubes Chiral BN nanotubes have also been grown, however too large to be modeled using traditional methods
12

Carbon doped BN nanotubes
BN(7,6) nanotube: chiral angle 27.46° 508 atoms in the unit cell BN(5,5) nanotube: 500 atoms in the supercell
Density of states
DOS obtained by a non-SCF calculation using diagonalization
13

FM or AFM coupling?
AFM coupling is favorable
Perspective
Our methods are rather general: can be extended to other O(N) methods based on DM and localized orbitals applied to metallic systems when combined with wavelets basis Many other interesting magnetic systems: single-molecule magnets magnetic nanostructures other spintronic systems
14

Conclusion
Two methods are proposed to deal with systems with or without predetermined spin multiplicity respectively. Detail and implement the PSUTC2 and SUTC2 algorithms. Carbon doped BN(7,6) nantoube has similar magnetic properties as carbon doped BN(5,5) nantoube. FM coupling is unfavorable for carbon doped BN nantoubes. Much large magnetic systems can be studied using these new methods.
JCP 123, 124105(2005)
Part I I O(N) calculation of maximally localized Wannier function with atomic basis set
15

MLWF
Definition for Wannier function:
Spread function: MLWFs minimize the function
For isolated systems, MLWF is identical to Boys localized orbital
N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12 847 (1997)
Applications
Understand the nature of the chemical bond Describe local and global dielectric properties A basis set for linear-scaling approaches and for constructing model Hamiltionians Electron transport, excited electronic states, and many body correlations
N. Marzari et al., Psi-K Newsletter 57, 129 (2003)
16

The nature of the chemical bond
Local and global dielectric properties
Modern theory of polarization: relates the vector sum of the centers of the Wannier functions to the polarization of an insulating system. Born charges and their decomposition external fields on periodic systems
R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993)
17

MLWFs as a basis set: from linear-scaling to model Hamiltonians
Wannier-like localized orbitals based O(N) methods Develop model Hamiltonians Hubbard, t-J model, LDA+U, ballistic conductance O(N) QMC methods
BaTiO3
Calculation of MLWFs: large cubic supercell
Maximize the functional minimal spread
Coordinate of the Wannier-function center (WFC)
P. L. Silvestrelli et al., Solid State Commun. 107, 7 (1998)
18

Traditional calculating methods:
First O(N3) ground state calculations Jacobi rotation Include several Jacobi sweeps In each sweep, N(N-1)/2 Jacobi rotations In each Jacobi rotation, calculation amount is O(N) In total, O(N3) Unitary transformation Need O(N3) diagonalization
Plane wave basis set O(N3)
G. Berghold et al., Phys. Rev. B 61, 10040 (2000)
O(N) methods for calculating MLWFs
Traditional methods scale like O(N3) Our main standpoints: Use numerical atomic orbital basis instead of planewave basis Use linear scaling density matrix methods instead of O(N3) traditional KS/HF methods Employ the localization property of MLWF
19

Main procedure
Obtaining DM using O(N) Trace-Correcting (TC2) method Obtaining Nocc independent non-orthogonal localized orbitals (Nocc : number of occupied orbitals) Performing a modified Lowdin orthogonalization to obtain orthogonal localized orbitals Carrying out Jacobi sweeps to obtain MLWF
DM from O(N) methods
Any DM or localized orbitals O(N) methods can be used. We use the Trace-Correcting purification method
A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002)
O(N)
20