New Castep Functionality

1. New Castep Functionality Linear response and non-local exchange-correlation functionals Stewart Clark University of Durham, UK

2. The Authors of Castep • Stewart Clark, Durham • Matt Probert, York • Chris Pickard, Cambridge • Matt Segal, Cambridge • Phil Hasnip, Cambridge • Keith Refson, RAL • Mike Payne, Cambridge

3. New Functionality • Linear response • Density functional perturbation theory • Atomic perturbations (phonons) • E-field perturbations (polarisabilities) • k-point perturbations (Born charges) • Non-local exchange-correlation • Summary of XC functionals • Why use non-local functionals? • Some examples

4. Density Functional Perturbation Theory • Based on compute how the total energy responds to a perturbation, usually of the DFT external potential v • Expand quantities (E, n, y, v) • Properties given by the derivatives

5. The Perturbations • Perturb the external potential (from the ionic cores and any external field): • Ionic positions  phonons • Cell vectors  elastic constants • Electric fields  dielectric response • Magnetic fields  NMR • But not only the potential, any perturbation to the Hamiltonian: • d/dk Born effective charges • d/d(PSP)  alchemical perturbation

6. Phonon Perturbations • For each atom i at a time, in direction a=x, y or z • This becomes perturbation denoted by l • The potential becomes a function of perturbation l • Take derivatives of the potential with respect to l • Hartree, xc: derivatives of potentials done by chain rule with respect to n and l So we need this E(2)

7. The expression for phonon-E(2) • Superscripts denote the order of the perturbation • E(2) given by 0th and 1st order wave functions and densities • This is a variational quantity – use conjugate gradients minimiser • Constraint: 1st order wave functions orthogonal to 0th order wave functions • This expression gives the electronic contribution

8. The Variational Calculation • E(2) is variational with respect to |y(1)> The plane-wave coefficients are varied to find the minimum E(2) under a perturbation • of a given ion i • in a given direction a • and for a given q • Analogous to standard total energy calculation • Based on a ground state (E(0)) calculation • Can be used for any q value

9. Sequence of calculation • Find electronic force constant matrix • Add in Ewald part • Repeat for a mesh of q And with Fourier interpolation: • Fourier transform to get F(R) • Fit and interpolate • Fourier transform and mass weight to get D at any q

10. Phonon LR: For and against • For • Fast, each wavevector component about the same as a single point energy calculation • No supercells requires • Arbitrary q • General formalism • Against • Details of implementation considerable

11. Symmetry Considerations • When perturbing the system the symmetry is broken • No time reversal symmetry • Implication is: k-point number increases • For example, Phonon-Si2 (Diamond): 6x6x6 MP set, 48 symmetry operations leads to SCF 28 k-points q=(0,0,0), 48 symmetry elements x-displacement leaves 12 elements 72 k-points for E(2) q=(1/2,0,0), 12 symmetry elements y-displacement leaves 4 elements 108 k-points for E(2) q=arbitrary, leaves only identity element and needs 216 k-points

12. So what can you calculate? …and phonon DOS • Phonon dispersion curves…

13. Thermodynamics • Phonon density of states • Debye temperature • Phase stability via • Entropic terms – derivatives of free energy • Vibrational specific heats

14. Electric Field Response • Bulk polarisability • Born effective charges • Phonon G-point LO/TO splitting • Dielectric permittivities • IR spectra • Raman Spectra

15. Why Electric Fields • Need Born effective charges to get LO-TO splitting: originates from finite dipole per unit cell • Example given later… • Found from d/dk calculation and similar cross-derivative expressions • Technical note: this means that all expressions for perturbed potentials different at zone centre than elsewhere

16. Examples of E-field Linear Response • Response of silicon to an electric field perturbation • Plot shows first order charge density Blue: where electrons are removed Yellow: where electrons go Acknowledgement: E-field work by my PhD student Paul Tulip

17. Phonon Dispersion of an Ionic Solid: G-point problems • At zone centre: Finite dipole (hence E-field) per unit cell caused by atomic displacements

18. Electronic Permittivity Tensor E-field on a polar system: NaCl • Without E-field: LO/TO Frequencies: 175 cm-1 • Born charge: Na Cl Ionic character as expected: Na+ and Cl- ions LO/TO splitting is 90cm-1: Smooth dispersion curve at the G-point • Polarisability Tensor Quantities given in atomic units

19. Summary of Density Functional Perturbation Theory • Phonon frequencies • Phonon DOS • Debye Temperate • PVT phase diagrams • Vibrational specific heats • Born effective charges • LO/TO Splitting • Bulk polarisabilities • Electric permittivities • First order charge densities (where electrons move from and to) • Etc… • Lots of new physics…and more planned in later releases

20. Non-local XC functionals • Basic background of XC interaction • Description within LDA and GGAs • Some non-local XC functional • Implementation within Castep • Model test cases

21. The exchange-correlation interaction Many body Hamiltonian – many body wavefunction Kohn-Sham Hamiltonian – single particle wavefunction Hartree requires self-interaction correction Single particle KE – not many body All goes in here

22. Approximations to XC • Local density approximation – only one • Generalised gradient approximations – lots • GGA’s + Laplacians – no real improvements • Meta-GGA’s (GGA + Laplacian + KE) – many recent papers – but nothing exciting yet. • Hybrid functionals (GGA/meta-GGA + some exact exchange from HF calculations) – currently favoured by the chemistry community. • Exc[n(r,r’)] – has been generally ignored recently by DFT community (although QMC and GW show it to have several interesting properties!).

23. New Functionals: The CPU cost • We aim to go beyond the local (LDA) or semi-local (GGA) approximations • Why? Cannot get any higher accuracy with these class of functionals • The computational cost is high: scaling will be O(N2) or O(N3) • This is because all pairs (or more) of electrons must be considered • That is, beyond the single particle model most DFT users are familiar with

24. First Class of New XC Functional: Exact and Screened Exchange • Exact Exchange (Hartree-Fock): • Screened Exchange: • And in terms of plane waves:

25. Why Exact/Screened Exchange? • Exact exchange gives us access to the common empirical “chemistry” GGA functionals such as B3LYP • Screened exchange can lead to accurate band gaps in semiconductors and insulators, so improved excitation energies and optical properties • The cost: scales O(N3), N=number of pl. waves • LDA/GGA is 0.1% of calculation, EXX is 99.9%

26. Silicon Band Structure • LDA band gap: 0.54eV • Exact Exchange band gap: 2.15eV • Screened Exchange band gap: 1.11eV • Experimental band gap: 1.12eV Extra cost can be worthwhile Some recent results (by my PhD student, Michael Gibson)

27. Another approach: Some theory on exchange-correlation holes… • The exact(!) XC energy within DFT can be written as • Relationship defined as the Coulomb energy between an electron and the XC hole nxc(r,r’) • XC hole is described in terms of the electron pair-correlation function Determines probability of finding an electron at position r’ given one exists at position r

28. Properties of the XC hole • Pauli exclusion principle: gxcobeys the sum rule • The size of the XC hole is exactly one electron – mathematically, this is the Pauli exclusion principle • The LDA and some GGAs obey this rule • For a universal XC function (applicable without bias), it must obey as many (all!) exact conditions as possible

29. The Weighted Density Method • In the WDA the pair-correlation function is approximated by • The weighted density is fixed at each point by enforcing the sum rule • This retains the non-locality of the function along with the Coulomb-like integral for Exc[n]

30. The XC potential • XC potential is determined in the usual manner (density derivative of XC energy) • We get 3 terms where

31. Example of WDA in action: An inhomogeneous electron gas • Investigate inhomogeneous systems by applying an external potential of the form • Very accurate quantum Monte Carlo results which to compare • Will have no pseudopotential effects • It’s inhomogeneous! • Given a converged plane wave basis set, we are testing the XC functional only – nothing else to consider

32. Investigate XC-hole shapes • XC holes for reference electron at a density maximum WDA results: acknowledgement to my PhD student Phil Rushton

33. But at the low density points… • Exchange-correlation hole at a density minimum

34. How accurate are the XC holes? • Compare to VMC data – M. Nekovee, W. M. C. Foulkes and R. J. Needs PRL 87, 036401 (2001)

35. Self Interaction Correction • H2+ molecule contains only one electron • HF describes it correctly, DFT with LDA fails • Self-interaction correction is the problem

36. Realistic Systems - Silicon • Generate self-consistent silicon charge density • Examine XC holes at various points Si [110] plane

37. XC holes in silicon Interstitial Region Bond Centre Region XC hole of electron moving along [100] direction

38. Some electronic properties • Band structure of Silicon – band gap opens

39. III-V semiconductor band gaps Band gaps in eV • The non-local potential opens up the band gap of some simple semiconductors

40. Timescales The following is my list of developments - other CGD members have other new physics and new improvements • First implementation of phonon linear response already in Castep v2.2 • Much faster (3-4 times faster) phonon linear response due to algorithmic improvements in Castep v3.0 • Phonon linear response calculations for metals under development. Aimed to be in v3.0 or v3.1 • Electric field response (Born charges, bulk polarisabilities, permittivity, LO/TO splitting) in Castep v3.0 • Exact and screened exchange: currently being tested and developed aimed at the v3.0 release • Weighted density method: currently working and tested – release schedule under discussion • Raman and IR intensities currently being implemented. Scientific evaluation still be performed before a possible release can be discussed