D G Kanhere, Center for Simulation and modeling Pune University, Mastani School- IISER Pune , July 2014 What is Molecular Dynamics (An Overview)

Ab Initio MD Born Oppenheimer Dynamics Car Parrinello Dynamics General Comment on How to set up dynamics run Thermodynamics of clusters Multiple histogram methods Results of clusters As Examples

System Many particles interacting with each other Solid/liquid : Box of Volume V, Temperature T and Pressure P: Periodic BC Finite size systems clusters :Free BC Given the forces acting on all the ions,

and initial state at time t=0, Compute Trajectories of all the particles as a function of time t, by using Newton's laws: Essentially exact We have entire phase space trajectories! Therefore all the information to compute various statistical quantities from microscopic description. Molecular Dynamics

Consider N particles interacting VIA Lennard Jones Potential (N ~ 10-13) For given co-ordinates {Ri}write down total energy : V ( R1,R2,) Compute Force on atom i. Write a subroutine which takes coordinates as input and returns total potential energy and forces on each of

the atoms. You are now ready to do simulation experiments after one more step. ow to integrate Newtons equations xI (t ) d 2 RI MI FI 2

dt t t (1) RI (t t ) RI (t ) vI (t )t 1 FI t 2 O(t 3 ) 2 MI

(2) RI (t t ) RI (t ) vI (t )t 1 FI t 2 O(t 3 ) 2 MI (1) (2) RI (t t ) RI (t t ) 2 RI (t ) FI t 2 O(t 4 ) MI

Verlets algorithm FI RI (t t ) 2 RI (t ) RI (t t ) t 2 O(t 4 ) MI typical molecular dynamics protocol Initialize: select starting atomic positions and velocities as close as possible to thermal equilibrium

Integrate: compute all forces and determine new positions using Verlets algorithm Equilibrate: let the system loose memory of the initial configurations, i.e. let it reach thermal equilibrium Average: accumulate averages of observables of interest (A) T 1 A A( R (t ), P(t )) dt T 0

A( R, P ) exp( E ) dRdP A exp( E )dRdP Average in Molecular Dynamics = Average in Statistical Mechanics

Molecular Dynamics Sum of pair interactions? How to treat chemical complexity? Where are the electrons? Empirical potentials Starting point: Pair potentials (LennardJones, Born-Mayer, Coulomb, etc) R0 12

4 R R0 R 6 Ze 2

e R , etc... , R + three-body corrections (3) ( I , J , K ) V0 cos 1 / 3 IJK

I K R RI RJ J + density dependent terms (embedded atom models) + atomic distorsion terms (includes polarization)

+ charge transfer terms Parameters are determined from empirical data, such as experimental EOS, vibrations, phase diagrams, dynamical properties, etc. Quantum simulations: The standard model Molecular dynamics for atoms

Ma = F = -dE/dR Schroedinger equation for electrons HE R. Cohen e--e- interactions: Density Functional Theory e--nuclei interactions: Pseudopotentials

Ab-initio molecular dynamics = Classical molecular dynamics in t potential energy surface generated by th electrons in their quantum ground state here are the electrons? The adiabatic approximation d 2 RI dE ({R}) MI FI 2

dt dRI Electrons respond much faster than nuclei to external forces, because of their lighter mass, therefore: Electrons are always in their instantaneous quantum mechanical ground state, for each given {R}. E R e H e R e Consequence:

Every MD step requires the calculation of the QM ground state of the electrons ee - e- eeNuclei

e- The electronic Hamiltonian He depends parametrically on the nuclear positions {R} Density functional theory Electron-electron (many-body) interaction

Density-functional theory [W. Kohn et al. 1964-1965] states that the e-e interaction can be written as a one-electron effective term: j e2 r e(r')

dr'Vxc (r e(ri )) ri r' ri rj However, the exact functional form of Vxc is not (yet) known. Current approximations to the exact Vxc go under different names (LDA, BLYP, BP, GGA, hybrid, et al). While GGA and hybrid functionals provide slightly better results than other approximations, the choice of the Vxc is often made in such a way to improve agreement with exps (ab fine??).

b-initio forces: the Hellmann-Feynman theorem d 0 H e ({R}) 0 d 2 RI dE ({R}) MI FI 2 dt dRI

dRI d 0 H e ({R}) 0 dRI 0 dH e ({R}) d 0

d 0 0 H e ({R}) 0 0 H e ({R}) dRI dRI dRI 0 dH e ({R}) d 0

d 0 0 E ({R}) 0 E ({R}) 0 dRI dRI dRI d 0 0 dH e ({R}) 0 0 E ({R})

dRI dRI dH e ({R}) 0 0 dRI Forces can be calculated without recalculating the ground state wave function for small atomic displacements

The Born Oppenheimer Dynamics How to keep electron in the ground state (1) d 2 RI dE ({R}) MI F

I dt 2 dRI E ({R}) 0 H e ({R}) 0 If we expand wavefunctions in a basis set cii i H e ci c*j i H e j

i, j finding the ground state is equivalent to minimizing a quadratic form in the {c}s (variational principle). So, standard minimization schemes can be used, e.g. steepest descent:

H e H e Carr- Parrinello dynamics How to keep electron in the ground state (2) d 2 RI dE ({R}) MI

F I dt 2 dRI E ({R}) 0 H e ({R}) 0 H e ({R})

down to the minimum for each {R} H e ({R}) ??? dE ({R})

M I RI dRI {R} The Car-Parrinello algorithm Car-Parrinello Molecular Dynamics (i) integrate the equations of motion on the (long) time scale set by the nuclear motion but nevertheless

(ii) take intrinsically advantage of the smooth time evolution of the dynamically evolving electronic subsystem as much as possible. The second point allows to circumvent explicit diagonalization or minimization to solve the electronic structure problem for the next molecular dynamics step as it is done in Born-Oppenheimer Molecular dynamics. CarParrinello molecular dynamics is an efficient method to satisfy requirement (ii) in a numerically stable fashion and makes an acceptable compromise concerning the length of the time step (i).

In CPMD a two-component quantum / classical problem is mapped onto a two-component purely classical problem with two separate energy scales at the expense of loosing the explicit time dependence of the quantum subsystem dynamics. Now, in classical mechanics the force on the nuclei is obtained from the derivative of a Lagrangian with respect to the nuclear positions. This suggests that a functional derivative with respect to the orbitals, which are interpreted as classical fields, might yield the force on the orbitals, given a suitable Lagrangian. In addition, possible constraints within the set of orbitals have to be imposed, such as e.g.

orthonormality (or generalized orthonormality conditions that include an overlap matrix). Car-Parrinello Lagrangian Car and Parrinello (1985) have postulated the following Lagrangian LCP 1 1

M R R i i | i EDFT { i } ij ( i | i ij ) R 2 i 2 ij Kinetic Energy Potential Energy Constrain ts

where i are fictitious masses and i are classical fields. Classical action needs to be minimized which results in equation of motions S L CP (t )dt d dL CP (t ) dL CP (t ) dt dR dR

d dL CP (t ) dL CP (t ) dt d i ( r, t ) d i ( r, t ) Car-Parrinello Equations of Motions Car and Parrinello (1985) have derived equations of motions 2 EDFT [{ i },{R}] Z Re

M R R r ( r, t ) dr R | r R (t ) | EDFT [{ i },{R}] ij j ( r, t ) i i ( r , t )

i ( r, t ) j H DFT i ( r, t ) ij j ( r, t ) j which are obviously transformed back to BornOppenheimer molecular dynamics if fictitious masses for the electrons i At the eqilibrium, there are no forces on electrons, therefore Ground state of density functional theory is reached with the

Why does the Car-Parrinello Method works? Conserved Energy in CPMD is not a physical energy but supplemented with a small fictitious kinetic term 1 E { } E phys M R R DFT i 2

R 1 1 Econs M R R i i | i EDFT [{ i }{R}] E phys T fict R 2 i 2 1 T fict i i | i i 2

Various Energies extracted from CPMD for a model system. EDFT Tfict The fictitious kinetic energy of the electrons is found to perform bound oscillations around a constant, i.e. the electrons do not heat up systematically in the presence of the nuclei; note that Tfict is a measure for deviations from the exact Born-Oppenheimer surface. Closer inspection shows actually two time scales of oscillations: the one visible in the Figure stems from the drag exerted by the moving nuclei on the electrons and is

the mirror image of the EDFT fluctuations. As a result the physical energy (the sum of the nuclear kinetic energy and the electronic total energy which serves as the potential energy for the nuclei) is essentially constant on the relevant energy and time scales. Given the adiabatic separation and the stability of the propagation, the central question remains if the forces acting on the nuclei are actually the correct" ones in Car-Parrinello molecular dynamics. As a reference serve the forces obtained from full self-consistent minimizations of the electronic energy at each time step, i.e. Born-Oppenheimer molecular dynamics with extremely well converged

wavefunctions. How to control adiabaticity? Since the electronic degrees of freedom are described by much heavier masses than the electronic masses, time step to perform CPMD simulations needs not to be too small as compared to Ehrenfest molecular dynamics. For a system with a gap in the spectrum, the lowest possible frequency of fictitious electronic oscillations 2 min ~ E gap

min E gap ~ 1/ 2

To guarantee adiabatic separation this frequency should be much larger than the typical phonon energy and/or the gap in the spectrum which would make sure that the electrons follow the nuclei adiabatically. Hence fictitious mass . At the same time, small fictitious mass would imply smaller and smaller time step because maximum fictitious electronic frequency is proportional to the plane-wave cutoff energy 2 max ~ Ecut

max E ~ cut 1/ 2 tmax ~

E cut 1/ 2 As a result a compromise fictitious mass needs to be found in CPMD simulations. For metals gap is zero and zero frequency fictitious electronic modes occur in the spectrum overlapping with the phonon spectrum. Thus, a well-controlled Born-Oppenheimer approach can only be recommended

P Method as dynamical solution of DFT equatio CP Method invented a new way to solve Kohn-Sham equations alternative to diagonalization. Consider variational principle which can be used to find an upper bound for the lowest eigenstates of the hamiltonian | H | cnn

using basis set expansion n c n n | H | m cm

nm 2 | c | n 1 nm CPMD offers a way to determine the coefficients without reduction to the eigenvalue problem.

orn-Oppenheimer versus- Car-Parrinello AIMD Born-Oppenheimer Car-Parrinello rn-Oppenheimer versus- Car-Parrinello forces rn-Oppenheimer versus- Car-Parrinello: summary

from M Sprik b-initio Molecular Dynamics: bibliography R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985) M. Payne, M. Teter, D. Allan, T. Arias, J. Joannopoulos, Rev. Mod. Phys. 64, 1045 (1992). D. Marx, J. Hutter, "Ab Initio Molecular Dynamics: Theory and Implementation", in "Modern Methods and Algorithms of Quantum Chemistry" (p. 301-449), Editor: J. Grotendorst, (NIC, FZ Jlich 2000) http://www.theochem.ruhr-uni-bochum.de/research/marx/cprev.en.html

R. Rousseau and S. Scandolo, Car-Parrinello Molecular Dynamics, in Encyclopedia of Condensed Matter Physics, edited by G. Bassani, G. Liedl, and P. Wyder, Elsevier, Amsterdam (2005) Finite Temperature Properties of finite size systems D G Kanhere Pune University Lindemann Criteria

Mean Square Displacements Extracting ionic density of States Constant Temperature Constant Volume Phase space trajectories The Multiple Histogram Method To Extract Density of States

Sodium Clusters Simplest of atomic clusters Jellium model works ( Depends) Nice delocalized charge density

Magic Clusters at N=8,20,40,58,92,138, Icosahedra for N=13,55,147, Chacko et al., Phys. Rev. B 71,155407 (2005) Lee et al., J. Chem. Phys. 123,164310 (2005) Lee et al., Phys. Rev. B (2007) Shahab et al., Phys Rev B ( 2007) Gazi et al J Chem Phys ( 2008) The heat Capacities N=40 : peaked

N=50 : flat N=55 : very sharp N=58 : peaked but broad BUT Highest melting Point ( Tm = 375 K ) Sodium Clusters

Simplest of atomic clusters Jellium model works ( Depends) Nice delocalized charge density Magic Clusters at N=8,20,40,58,92,138, Icosahedra for N=13,55,147, Chacko et al., Phys. Rev. B 71,155407 (2005) Lee et al., J. Chem. Phys. 123,164310 (2005)

Lee et al., Phys. Rev. B (2007) Shahab et al., Phys Rev B ( 2007) Gazi et al J Chem Phys ( 2008) Application II Thermodynamics Size sensitive specific heats Magic melters Higher than Bulk Melting temperatures

The Case of Gallium and Aluminum clusters Gallium and aluminum clusters show extreme size sensitivity Gallium clusters melt at Higher than their

Bulk Tm Magic Melters ----- Non - Melter Gallium Clusters: Heat Capacity Experimental Data from Indiana group * N=30 Flat N=31 Peaked * Tm Shifts by 200K @ N=45

And by 350K across the series Calculations Density functional with LDA/GGA Born Oppenheimer MD

Delta T 100 au Total simulation time per temp 100 ps or more Soft pseudo potentials, plane waves Extensive search for equilibrium geometries. Simulated annealing, Basin Hopping The Heat Capacity : Ga30- Ga31 Joshi et al. Phys. Rev. Lett. 96,135703 (2006)

The MSD : Ga30 and Ga31 Finite Size Effect : Amorphous Clusters In amorphous cluster each atom may have different environment. Atoms may be bonded with the rest of the cluster with different strength. When heated, they will begin diffusive motion at different temperatures This may result in continuous phase change.

No sharp peak in specific heat, very broad transition. Cluster with large island having local order will show a peak in the heat capacity ..Most of the atoms melt together. Thermodynamics: Ga27Si3 Ghazi and Kanhere J Phys Chem (2012) Ga13 Ga12C and Al13 - Al12C Ga13 is Decahedra Ga12C is a Perfect

Icosahedra Calculated Specific Heat Ga Specific Heat