First-principles molecular dynamics P. Giannozzi Dept of Chemistry, - PowerPoint PPT Presentation

First-principles molecular dynamics P. Giannozzi Dept of Chemistry, Physics, Environment, University of Udine, Italy XIX School of Pure and Applied Biophysics, Venice, 26 January 2015 Outline 1. First-principles molecular dynamics: when and what? 2. A quick look at the theoretical basis, with emphasis on Density-Functional Theory 3. Born-Oppenheimer and Car-Parrinello Molecular Dynamics 4. Technicalites you should be aware on: plane waves, pseudopotentials 5. Software and needed hardware – Typeset by Foil T EX –

1. What is First-Principles Molecular Dynamics? In First-Principles Molecular Dynamics (FPMD) the forces acting on nuclei are computed from the electronic states , i.e. by solving, in principle exactly , the many-body Schr¨ odinger equation for electrons (nuclei are approximated by classical particles). FPMD is the method of choice whenever “force fields” do not give a sufficiently accurate or reliable description of the system under study, or of a specific phenomenon: typically, whenever bonds are broken or formed, or in presence of complex bonding, e.g.: transition metals. With FPMD, all techniques used in classical (with force fields) MD can be used; moreover , one has access to the electronic structure and charge density, can in principle compute excitation spectra: however , the heavy computational load of FPMD, in spite of available efficient implementations, limits its range of application to systems containing O (1000) atoms max, for rather short times (tens of ps).

An example of application Below, the simplest model for the myoglobin active site: iron-porphyrin- imidazole complex. Yellow: Fe. Dark Gray: C. Blue: N. Light gray: H. To the right: extended model of the myoglobin active site (including the 13 residues in a 8A sphere centered on Fe, terminated with NH 2 groups, containing 332 atoms, 902 electrons). Red: O, all other atoms as in the above picture.

The need for first-principles MD • High-level theoretical description needed, due to the presence of transition metal atoms • Complex energy landscape : many possible spin states and local minima, not easy to guess from static (single-point) calculations – Molecular Dynamics may be better suited. Interest: respiration, photosynthesis, enzymatic catalysis... Note however that sizable systems (hundreds of atoms in this case) are needed in order to accurately describe the effect of surrounding environment: efficient first-principles techniques are needed (P. Giannozzi, F. de Angelis, and R. Car, J. Chem. Phys. 120 , 5903 (2004); D. Marx and J. Hutter, Modern Methods and Algorithms of Quantum Chemistry , John von Neumann Institute for Computing, J¨ ulich, NIC Series, Vol.3, pp. 329-477 (2000) http://juser.fz-juelich.de/record/152459/files/FZJ-2014-02061.pdf )

2. The basics: Born-Oppenheimer approximation The behavior of a system of interacting electrons r and nuclei R is determined by the solutions of the time-dependent Schr¨ odinger equation :   h∂ ˆ � h 2 � h 2 ∂ 2 ∂ 2 Φ( r , R ; t ) ¯ ¯  ˆ  − i ¯ = − + V ( r , R ) Φ( r , R ; t ) ∂ r 2 ∂ R 2 ∂t 2 M µ 2 m µ i µ i where V ( r , R ) is the potential describing the coulombian interactions: � � � Z µ Z ν e 2 Z µ e 2 e 2 V ( r , R ) = | R µ − R ν | − | r i − R µ | + | r i − r j | µ>ν i,µ i>j ≡ V nn ( R ) + V ne ( r , R ) + V ee ( r ) Born-Oppenheimer (or adiabatic ) approximation (valid for M µ >> m ): Φ( r , R ; t ) ≃ Φ( R )Ψ( r | R ) e − i ˆ ˆ Et/ ¯ h NB: in order to keep the notation simple: r ≡ ( r 1 , .., r N ) , R ≡ ( R 1 , .., R n ) ; spin and relativistic effects are not included

Potential Energy Surface The Born-Oppenheimer approximation allows to split the original problem into an electronic problem depending upon nuclear positions : � � � h 2 ∂ 2 ¯ − Ψ( r | R ) = E ( R )Ψ( r | R ) + V ( r , R ) ∂ r 2 2 m i i and a nuclear problem under an effective interatomic potential , determined by the electrons:   � h 2 ∂ 2 ¯  Φ( R ) = ˆ  − + E ( R ) E Φ( R ) ∂ R 2 2 M µ µ µ E ( R ) defines the Potential Energy Surface (PES). The PES determines the motion of atomic nuclei and the equilibrium geometry: the forces F µ on nuclei, defined as F µ = − ∂E ( R ) , ∂ R µ vanish at equilibrium: F µ = 0 .

Solving the electronic problem: Hartree-Fock Let us write the many-body electronic wave-function as a Slater determinant : � � � � ψ 1 (1) . . . ψ 1 ( N ) � � � � 1 . . . � � √ Ψ( r 1 , r 2 , ... r N ) = � � . . . N ! � � � � ψ N (1) . . . ψ N ( N ) for N electrons, where ( i ) ≡ ( r i , σ i ) , σ i is the spin of the i − th electron, the ψ i states are all different. Minimization of E yields the Hartree-Fock equations . For the “restricted” case (pairs of states with opposite spin): h 2 − ¯ 2 m ∇ 2 ψ i ( r ) + V ( r ) ψ i ( r ) + V H ( r ) ψ i ( r ) + ( ˆ V x ψ i )( r ) = ǫ i ψ i ( r ) where i = 1 , ..., N/ 2 , V ( r ) is the external (nuclear) potential acting on each electron: � Z µ e 2 V ( r ) = − | r − R µ | , µ

V H is the Hartree potential: � � e 2 | ψ j ( r ′ ) | 2 | r − r ′ | d r ′ V H ( r ) = 2 j ˆ V x is the (nonlocal) exchange potential: � � e 2 ( ˆ | r − r ′ | ψ ∗ j ( r ′ ) ψ i ( r ′ ) d r ′ . V x ψ i )( r ) = − ψ j ( r ) j These are integro-differential equations that must be solved self-consistently . In practice, HF is not a very good approximation. More accurate results require perturbative corrections (Møller-Plesset) or to add more Slater determinants (Configuration Interaction) in order to find more of the correlation energy , defined as “the difference between the true and the Hartree-Fock energy”. Solving HF and post-HF equations is the “traditional” approach of Quantum Chemistry. An alternative approach, based on the charge density , is provided by Density-Functional Theory (DFT)

Solving the electronic problem: Hohenberg-Kohn theorem Let us introduce the ground-state charge density n ( r ) . For N electrons: � | Ψ( r , r 2 , ... r N ) | 2 d r 2 ...d r N . n ( r ) = N The Hohenberg-Kohn theorem (1964) can be demonstrated: there is a unique potential V ( r , R ) having n ( r ) as ground-state charge density. Consequences: • The energy can be written as a functional of n ( r ) : � E [ n ( r )] = F [ n ( r )] + n ( r ) V ( r ) d r where F [ n ( r )] is a universal functional of the density, V ( r ) is the external (nuclear) potential acting on each electron: � Z µ e 2 V ( r ) = − | r − R µ | . µ • E [ n ( r )] is minimized by the ground-state charge density n ( r ) . Note: the energy E defined above does not include nuclear-nuclear interaction terms

Density-Functional Theory: Kohn-Sham approach Let us introduce the orbitals ψ i for an auxiliary set of non-interacting electrons whose charge density is the same as that of the true system: � | ψ i ( r ) | 2 , n ( r ) = � ψ i | ψ j � = δ ij i Let us rewrite the energy functional in a more manageable way as: � E = T s [ n ( r )] + E H [ n ( r )] + E xc [ n ( r )] + n ( r ) V ( r ) d r where T s [ n ( r )] is the kinetic energy of the non-interacting electrons: � � h 2 T s [ n ( r )] = − ¯ i ( r ) ∇ 2 ψ i ( r ) d r , ψ ∗ 2 m i E H [ n ( r )] is the Hartree energy , due to electrostatic interactions: � n ( r ) n ( r ′ ) E H [ n ( r )] = e 2 | r − r ′ | d r d r ′ , 2

E xc [ n ( r )] is called exchange-correlation energy (a reminiscence from the Hartree- Fock theory) and includes all the remaining (unknown!) energy terms. Minimization of the energy with respect to ψ i yields the Kohn-Sham (KS) equations: � � h 2 − ¯ 2 m ∇ 2 + V ( r ) + V H ( r ) + V xc ( r ) ψ i ( r ) = ǫ i ψ i ( r ) , � �� � H KS where the Hartree and exchange-correlation potentials: � n ( r ′ ) V H ( r ) = δE H [ n ( r )] V xc ( r ) = δE xc [ n ( r )] = e 2 | r − r ′ | d r ′ , δn ( r ) δn ( r ) depend self-consistently upon the ψ i via the charge density. The energy can be rewritten in an alternative form using the KS eigenvalues ǫ i : � � ǫ i − E H [ n ( r )] − E = n ( r ) V xc ( r ) d r + E xc [ n ( r )] i

Exchange-correlation functionals: simple approximations What is E xc [ n ( r )] ? Viable approximations needed to turn DFT into a useful tool. The first, ”historical” approach (1965) is the Local Density Approximation (LDA): replace the energy functional with a function of the local density n ( r ) , � � � V xc ( r ) = ǫ xc ( n ( r )) + n ( r ) dǫ xc ( n ) � E xc = n ( r ) ǫ xc ( n ( r )) d r , � dn n = n ( r ) where ǫ xc ( n ) is calculated for the homogeneous electron gas of uniform density n (e.g. using Quantum Monte Carlo) and parameterized. Generalized Gradient Approximation (GGA): A more recent class of functionals depending upon the local density and its local gradient |∇ n ( r ) | , � E xc = n ( r ) ǫ GGA ( n ( r ) , |∇ n ( r ) | ) d r There are many flavors of GGA, yielding similar (but slightly different) results. GGA is the ”standard” functional in most FPMD calculations, with excellent price-to- performance ratio, but some noticeable shortcomings.