Hybrid Quantum Mechanics / Molecular Mechanics (QM/MM) Approaches - - PowerPoint PPT Presentation

Hybrid Quantum Mechanics / Molecular Mechanics (QM/MM) Approaches - Partitioning the system: QM/MM Hamiltonian and calculation of forces Mauro Boero Institut de Physique et Chimie des Matériaux de Strasbourg University of Strasbourg - CNRS, F-67034 Strasbourg, France and @Institute of Materials and Systems for Sustainability, Nagoya University - Oshiyama Group, Nagoya Japan 1

Splitting of the system into a QM inner region and an outer MM region • The two subsystems are in strong interaction with each other • Hence the total energy (Hamiltonian) is not the simple (linear) sum of the energies of the two sub- parts • Coupling interaction must be accounted for • Particular care must be taken at the boundary especially if the border cuts through chemical bonds (more later) 2

Additive scheme: Total Hamiltonian   S     tot QM int QM MM MM H H H [ , { R } , { R } ] H J I S = QM+MM 3

MD Simulations: how to construct a Force Field? We need to consider all the relevant motions of the system that we want to study bond stretch   0 2 v ( r , k ) k ( r r ) 1 IJ IJ 2 IJ IJ IJ torsion v ( f IJKL ) intermolecular interactions bending intramolecular nonbonded      0 2 v ( ) 1 k ( )  IJK IJK 2 interactions 4

Typical form of the MM potential :    2    MM R MM 0 V ( ) k ( R R ) d Bond stretching I IJ I J IJ  I J    2     0 k ( R ) Bond bending  I  I J      f   k 1 cos( n ( R ) ) Torsion angles f I  I J 1 q I q   J Coulomb interaction   4 R R  non bond IJ 0 I J   6 12              Van der Waals     IJ IJ     IJ    R R R R         non bond IJ I J I J 5

MD Simulations: How do we get the parameters k IJ , k q , etc… for the force field ? From experiments molecular and bulk properties X-ray / neutron scattering structure factors isotopic substitutions, etc… From ab initio calculations molecular and cluster properties static geometry optimization (minima, saddle points) first principles methods 6

Subtractive scheme (e.g. IMOMM/ONIOM) (K. Morokuma et al. J. Comp. Chem . 16 , 1170 (1995)) 1. Compute H MM ( S ) 2. Compute H QM (QM) of the QM subsystem 3. Compute H MM (QM), i.e. MM calculation of the QM subsystem 4. Sum up terms 1 and 2 and subtract term 3 to get rid of the double counting H tot ( S ) = H MM ( S )+ H QM (QM)- H MM (QM) Warning: The coupling between QM and MM is driven by MM and this can give problems in Coulomb interactions between QM and MM, especially if moving (MD). 7

Typical form of the MM Hamiltonian :   MM 1  2    MM MM MM MM MM H M R V ( R ) I I I 2 I Classical kinetic term Parameterized form of the electrostatic interactions, i.e. an analytical function of the positions only ( not of the velocities ) keeping into account many-body effects up to the third (bending modes) or fourth (torsion modes) order. 8

QM driver: quantum electrons and classical nuclei and to compute forces from fundamental quantum mechanics y i ( x ) electron- Beside the atom-atom electron electron (ion-ion) interaction, - ion new ingredients must R I ion-ion be included electron - ion 9

Density Functional Theory: brief review • Define the electronic density  ( x ) as a superposition of single particle Kohn-Sham (KS) orbitals • Write the total energy functional as E [ y i , R I ] = E k + E H + E xc + E ps + E M i.e. sum of the interactions electron-electron + electron-ion + ion-ion after W. Kohn* and L. J. Sham, Phys. Rev. 140 , A1133 (1965) *Nobel Prize for Chemistry in 1998 10

Density Functional Theory: brief review • Electron-electron interaction: • E k = kinetic energy, E H = Coulomb interaction, E xc = exchange-correlation interaction (= what we do not know: The many-body interaction) 11

Density Functional Theory: brief review • Electron-ion interaction:    r ( x ) E ps   3 ps d x V x R I I the core-valence interaction is described by pseudopotentials • Ion-ion interaction:  1 Z I Z   E M J  2 R R  I J I J 12

(Not so) Typical form of the QM Lagrangean (Ok, and Hamiltonian): • In CPMD we generally use the extended MD Lagrangean where we add the electronic degrees of freedom y i • And a constraint to keep the orthogonality of the wavefunctions • And any external additional variables a q (e.g. thermostats, stress, etc.) 13

Car-Parrinello Molecular Dynamics 14

Car-Parrinello Molecular Dynamics • Solve the related Euler-Lagrange EOM   CP CP d L L   0 Electrons  y  y * * dt  i i   CP CP d L L   0 Ions    dt R R I I   CP CP d L L   0 external additional     dt  variables q q 15

Car-Parrinello Molecular Dynamics • Solve the related Euler-Lagrange EOM Temperature, pressure, etc. 16

Car-Parrinello Molecular Dynamics • It is of course straightforward to give a Hamiltonian, instead of a Lagrangean formulation, by a simple Legendre transform after defining the momenta   CP CP L L     y     y * * ( x )  ( x ) ( x )  ( x ) i i i i  y  y *  ( x )  ( x ) i i     CP p L M R  i I I R I  CP L        q q q    q 17

Car-Parrinello Hamiltonian for QM • so that the QM Hamiltonian reads   2   2 *   p  ( x ) ( x )       q    CP 3 tot I i i H d x E , { R }, I q   2 M 2 I i q I q       y y   3 * d x ( x ) ( x ) ij i j ij ij and the CPMD equations of motion will be given by the corresponding Hamilton equations for each set of variables and momenta. 18

CPMD dynamics: a little warning BO surface CP trajectory BO trajectory The difference between the CP trajectories R I CP (t) and the Born-Oppenheimer (BO) ones R I BO (t) is bound by BO (t)| < C  1/2 | R I CP (t) - R I           LUMO HOMO (C > 0) if 2 0 0 See F.A. Bornemann and C. Schütte, Numerische Mathematik vol. 78 , N. 3 , p. 359-376 (1998) 19

Practical implementation • To implement the CP-EOM numerically, the KS orbitals are generally expanded in plane waves y i ( x ) = S G c i ( G ) e i Gx • G are the reciprocal space vectors. The Hilbert space spanned by PWs is truncated to a suitable cut-off E cut such that G 2 /2 < E cut • The basis set accuracy can be systematically improved in a fully variational way • PWs are independent from atomic positions (no Pulay forces) • However they are not localized and can be inefficient, e.g., for small molecules or clusters in a large simulation cell 20

Plane wave expansion: y i ( x ) = S G c i ( G ) e i Gx For each electron i =1,…, N , G=1,…, M are the reciprocal space vectors. The Hilbert space spanned by PWs is truncated to a cut-off G cut 2 /2 < E cut R space a G space E 1 cut E 2 cut > E 1 cut 21

Practical implementation • Verlet’s algorithm on EOM gives ( m / D t 2 ) ・ [ | y i ( t + D t )> + | y i ( t - D t )> - 2 | y i ( t )>] = -( H CP + L ) | y i ( t )> or, in G -space,     D   D   c ( G , t t ) c ( G , t t ) 2 c ( G , t ) i i i D 2 t        CP G H G c ( G ) c ( G ) i ij j  G j • The ionic degrees of freedom R I ( t ) are updated at a rate (speed) D t while the electronic degrees of freedom | y i ( t )> are updated at a rate D t / m 1/2 ( D t ~ 5 a.u., m ~ 500 a.u.), hence they are much slower (adiabatic) with respect to the ions. • However this D t is driven by QM and is much shorter than  t in MM simulations 22

Typical form of the QM-MM Hamiltonian : Basically just the electrostatic interaction between the two subsystems, i.e.  QM MM ( x ) MM q Z         int 3 H ( x ), { R }{ q , R } q d x I J J I I I   x R | R R | QM    I 1 I 1 J 1 I I J QM atom-MM atom R I with charge q I QM electron-MM atom Particle-mesh interaction: r (x) We need good charge models for the MM part and a lot of computation 23

Total hamiltonian: scaling in a plane wave basis set   S     tot QM int QM MM MM H H H [ , { R } , { R } ] H J I Pure QM interaction N el x N G x N G Major cost = with non-bonding y j ( x ) = S G c j ( G ) e i Gx interactions* MM x MM QM/MM interface MM x N G * Ewald summations ( O ( MM log MM )) or spherical cut-off techniques ( O ( MM )) can reduce the MM computational cost. 24