Applications of Density-Functional Theory: Structure Optimization, Phase Transitions, and Phonons

1. Applications of Density-Functional Theory: Structure Optimization, Phase Transitions, and Phonons Christian Ratsch UCLA, Department of Mathematics In previous talks, we have learned how to calculate the ground state energy, and the forces between atoms. Now, we will discuss some important concepts and applications of what we can do with this.

2. Outline • Structure Optimization • Optimize bond length • Optimize atomic structure of a cluster or molecule • Optimize structure of a surface • Phase transitions • Phonons • Structural and vibrational properties of metal clusters • Dynamics on surfaces • Molecular dynamics • Use transition state theory: DFT can be used to calculate energy barriers, prefactors.

3. F F F F • Algorithms for structure optimization: • Damped Newton dynamics (option: NEWT) • Parameters l and m need to be optimized. • Quasi-Newton structure optimization (BFGS scheme; option: PORT) • No parameters needed • In principle, only one iteration needed if Hessian is known, and fully harmonic. • In practise, a few iterations and Hessian updates are needed. • Nevertheless, this is typically the recommended option. 1.80 Å Hessian: My results for vanadium: lbond = 1.80 Å (experimental value: 1.77 Å) Structure optimization; example: Vanadium dimer • Put atoms “anywhere” • Calculate forces • Forces will move atoms toward configuration with lowest energy (forces = 0)

4. E=0.4 eV E=0.8 eV Bigger vanadium clusters: V8+ Different start geometries lead to very different structures E=0 eV E=1.8 eV • Each structure is in an energetic local minimum (i.e., forces are zero). • But which one is the global minimum? • Finding the global minimum is a challenging task • Sometimes, good intuition is all we need • But even for O(10) atoms, “intuition” is often not good enough. • There are many strategies to find global minima.

5. Surface relaxation on a clean Al(111) surface Relaxation obtained with DFT (DMol3 code) (Jörg Behler, Ph.D. thesis) D12 D23 • Top layer relaxes outward • Second layer relaxes inward (maybe)

6. In - terminated As - terminated [110] [110] Which reconstruction ? Use density-functional theory (DFT) Things are more complicated on semiconductor surfaces • Semiconductor surfaces reconstruct. Example InAs(100) • Surface reconstruction is important for evolution of surface morphology which influences device properties • RHEED experiments show transition of symmetry from (4x2) to (2x4)

7. Computation details of the DFT calculations • Computer code used: fhi98md • Norm-conserving pseudopotentials • Plane-wave basis set with Ecut = 12 Ryd • k summation: 64 k points per (1x1) cell • Local-density approximation (LDA) for exchange-correlation • Supercell with surface on one side, pseudo-hydrogen on the other side • Damped Newton dynamics to optimize atomic structure

8. Possible structures b(2x4) a(2x4) z(4x2) b2(2x4) a2(2x4) a3(2x4) b3(2x4) We also considered the corresponding (4x2) structures (which are rotated by 90o, and In and As atoms are interchanged)

9. b2 (4x2) a2 (2x4) a3 (2x4) b2 (2x4) z(4x2) Typical experimental regime c (4x4) Phase diagram for InAs(001) C. Ratsch et al., Phys. Rev. B 62, R7719 (2000).

10. Low As pressure High As pressure a2(2x4) b2(2x4) Predictions confirmed by STM images Barvosa-Carter, Ross, Ratsch, Grosse, Owen, Zinck, Surf. Sci. 499, L129 (2002)

11. Without pressure, structure B is stable • With pressure, eventually structure A becomes stable • Minimize Gibbs free energy (at T=0): Pressure for phase transition is determined by: Structure B V1 V2 Phase transitions Yesterday, we have learned how to calculate the lattice constant, by calculating E(V). Etot Structure A lattice constant volume

12. Classical example: Phase transition of silicon Si has a (cubic) diamond structure, which is semiconducting The pressure of phase transition has been computed from DFT to be 99 kbar (experimental value: 125 kbar) Under pressure, there is a phase transition to the tetragonal b-tin structure, which is metalic M.T. Yin, and M.L. Cohen, PRL 45, 1004 (1980)

13. Ying and Cohen, PRL, 1980 Ying and Cohen, PRB, 1982 Historical remark • In the original paper, the energies for bcc and fcc were not fully converged • Luckily, this “did not matter” (for the phase transition) • Nevertheless, these calculations are considered one of the first successes of (the predictive power) of DFT calculations.

14. n-2 n-1 n n+1 n u(na): displacement of atom n Equation of Motion: Mass of atom Spring constant Assume solution of form: Lattice vibrations: A 1-dimensional monatomic chain of atoms Periodic boundary condition requires: Upon substitution, we get solution Dispersion curve for a monatomic chain

15. n-1 n n+1 n u1(na): displacement of atom n,1 u2(na): displacement of atom n,2 Spring K Spring G Lattice vibrations of a chain with 2 ions per primitive cell Coupled equations of motion: Solution

16. Because for small k, which is characteristic of sound waves Dispersion relation for the diatomic linear chain There are N values of k: • For each k, there are 2 solutions, leading to a total of 2N normal modes. • The normal modes are also called “phonons”, in analogy to the term “photons”, since the energy of the N elastic modes are quantized as Optical branch Because long wavelength modes can interact with electromagnetic radiation Acoustic branch

17. Analysis essentially the same • For each k, there are 3p normal modes • The lowest 3 branches are acoustic • The remaining 3(p-1) branches are optical • A “real” phonon spectrum might look slightly different; • The reason is that interactions beyond nearest neighbors are not included, the potential might not be harmonic, there are electron-phonon coupling, etc. • More in the talk by Claudia Ambrosch-Draxl Lattice vibrations in 3D with p ions per unit cell

18. How can we calculate phonon spectrum? • Frozen phonon calculation • This is what you will do this afternoon. • DFT perturbation theory • Molecular Dynamics • Do an MD simulation for a sufficiently long time • “Measure” the time of vibrations; for example for dimer, this is obvious • For bigger systems, one needs to do Fourier analysis to do this • Very expensive

19. Calculate dynamical matrix D • in principle, this is done by displacing each atom in each direction, and get the forces acting on all other atoms: • eigenvalues l are the frequencies • in practise, one exploits the symmetry of the system (need group theory) • Repeat for several k p/3a p/2a p/a Frozen phonon calculations • Choose a supercell that corresponds to the inverse of wave vector k • More details in the presentation by Mahboubeh Hortamani

20. Structural and vibrational properties of small vanadium clusters • Why do we care about small metal clusters? • Many catalytic converters are based on clusters • Clusters will play a role in nano-electronics (quantum dots) • Importance in Bio-Chemistry • Small clusters (consisting of a few atoms) are the smallest nano-particles! This work was also motivated by interesting experimental results by A. Fielicke, G.v.-Helden, and G. Meijers (all FHI Berlin)

21. 14 3 4 15 5 16 6 17 7 18 8 19 9 20 10 21 11 22 12 23 13 Spectra for VxAry+ • Each cluster has an individual signature • V13+ is the only structure with peaks that are beyond 400 cm-1 • Beginning at size 20, the spectra look “similar”. This suggests a bulk-like structure

22. Tunable free electron Laser (FELIX) Example: Excitation of V7Ar1 and V7Ar2 at 313 cm-1 Experimental setup using a tunable free electron laser Laser Beam: clusters are formed, Ar attaches Mass-Spectrometer Gas flow (~1% Ar in He) metal-rod

23. DFT calculations for small metal clusters • Computer Code used: DMol3 • GGA for Exchange-Correlation (PBE); but we also tested and compared LDA, RPBE • We tested a large number of possible atomic structures, and spin states. • All atomic structures are fully relaxed. • Determine the energetically most preferred structures • Calculate the vibrational spectra with DFT (by diagonalizing force constant matrix, which is obtained by displacing each atom in all directions) • Calculate the IR intensities from derivative of the dipol moment • What can we learn from these calculations? • Confirm the observed spectra • Determine the structure of the clusters • Is the spectrum the result of one or several isomers?

24. Structure determination for V8+ experiment theory E=0 E=0.4eV E=0.8eV E=1.8eV

25. Structure determination for V9+ experiment theory E=0 S=1 E=0.01eV S=0 E=0.06eV S=1 E=0.08eV S=0

26. Niobium • Open issues: • Sometimes neutral and cationic niobium have similar spectrum, sometimes they are very different • Cationic Nb is sometimes like cationic V, sometimes different.

27. Calculated Spectra of lowest energy structure neutral cationic Niobium 7 Experimental Spectra neutral cationic

28. Calculated Spectra of lowest energy structure neutral cationic Niobium 6 Experimental Spectra neutral cationic

29. But there is the timescale problem: Even with big computers, we can’t describe dynamics beyond microseconds, for up to ~ 106 atoms 10-13 s More about molecular dynamics in the talk by Karsten Reuter. One way to “speed up dynamics” is to use transition state theory (TST) Molecular dynamics Once we have the forces, we can solve the equation of motion for a large number of atoms, describing the dynamics of a system of interest.

30. Attempt frequency (using harmonic approximations) Energy barrier normal mode frequencies at adsorption and transition site , • Finding transition state is a big challenge • Sometimes, intuition is enough • Often, sophisticated schemes are needed (nudged elastic band method, dimer method, …. ) Transition state theory (TST) to calculate microscopic rate parameters Transition state theory (Vineyard, 1957):

31. Ag/Pt(111) Ag/ 1ML Ag/Pt(111) Ag/Ag(111) o T = 65 K, Coverage = 0.12 ML 100 A Nucleation Theory: N ~ (D/F)-1/3 Model system: Ag/Ag(111) and Ag/Pt(111) (Brune et al, Phys. Rev. B 52, 14380 (1995))

32. Tensile strain: less diffusion • Compressive strain: increased diffusion Results and comparison: Diffusion barrier Lowered diffusion barrier for Ag on Ag/Pt(111) is mainly an effect of strain. Ratsch et al., Phys. Rev. B 55, 6750 (1997)

33. Attempt frequency • Calculate force constant matrix by displacing each atom in x,y,z-direction, and by calculating the forces that act on all atoms • Eigenvalues of this matrix are the normal mode frequencies • Important question: how many degrees of freedom need to be included? , Convergence test for Ag/Ag(111) only adatom 2x2 cell adatom and top layer, 2x2 cell adatom and 2 layers, 4x4 cell # degrees of freedom 3 15 99 Prefactor G0 (THz) 1.55 0.82 0.71 How to calculate the prefactor Ratsch et al, Phys. Rev. B 58, 13163 (1998)

34. Results and comparison: Prefactor Compensation effect: Higher barrier Higher prefactor • Prefactor is O(1 THz) for all systems • Compensation effect is not confirmed • Explanation for experimental result: • Long range interactions are important for systems with small barrier (Fichthorn and Scheffler, PRL, 2000; Bogicevic et al., PRL, 2000), • Simple nucleation theory does not apply any longer but can be modified (Venables, Brune, PRB 2002)

35. Conclusion and summary • DFT calculations can be used to optimize the atomic structure of a system • DFT calculations can be used to calculate the pressure of a phase transition. This will be part of the exercises this afternoon! • DFT calculations can be used to calculate a phonon spectrum. This will also be part of the exercises this afternoon. • DFT calculations can be used to obtain structural and vibrational properties of clusters • DFT calculations can be used to obtain the relevant microscopic parameters that describe the dynamics on surfaces.