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

1

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

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

3

S = QM+MM

Additive scheme: Total Hamiltonian

 

MM MM QM int QM tot

] } { , } { , [ H H H H

I J

   S R R 

4

MD Simulations: how to construct a Force Field?

We need to consider all the relevant motions of the system that we want to study intermolecular interactions intramolecular nonbonded interactions torsion bond stretch bending

2 2 1

) ( ) , (

IJ IJ IJ IJ IJ

r r k k r v  

2 2 1

) ( ) (   



 

IJK IJK

k v

v(fIJKL)

5

Typical form of the MM potential:  





  

J I IJ J I IJ I

d k V

2 MM

) ( ) ( R R RMM

 





 

J I I

k

2

) (  



R

 





  

J I I

n k ) ) ( cos( 1  f

f

R





 

IJ bond non J I J Iq

q R R 4 1 





                            

IJ bond non J I IJ J I IJ IJ 12 6

R R R R   

Bond stretching Bond bending Torsion angles Coulomb interaction Van der Waals

6

MD Simulations: How do we get the parameters kIJ, kq, 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

7

Subtractive scheme (e.g. IMOMM/ONIOM)

(K. Morokuma et al. J. Comp. Chem. 16, 1170 (1995))

• 1. Compute HMM(S)
• 2. Compute HQM(QM) of the QM subsystem
• 3. Compute HMM(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 Htot(S) = HMM(S)+ HQM(QM)-HMM(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).

8

Typical form of the MM Hamiltonian:

 

) ( 2 1

MM MM MM 2 MM MM MM I I I I

V M H R R



  

Classical kinetic term Parameterized form of the electrostatic interactions, i.e. an analytical function of the positions only (not

• f the velocities) keeping into account many-body

effects up to the third (bending modes) or fourth (torsion modes) order.

9

QM driver:

quantum electrons and classical nuclei and to compute forces from fundamental quantum mechanics electron- electron ion-ion electron

• ion

electron

• ion

yi(x) RI

Beside the atom-atom (ion-ion) interaction, new ingredients must be included

10

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[yi,RI] = Ek + EH + Exc + Eps + EM

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

11

Density Functional Theory: brief review

• Electron-electron interaction:
• Ek = kinetic energy, EH = Coulomb interaction,

Exc = exchange-correlation interaction (= what we do not know: The many-body interaction)

12

Density Functional Theory: brief review

• Electron-ion interaction:

the core-valence interaction is described by pseudopotentials

• Ion-ion interaction:

 



 

I I ps

V x d R x

3

Eps r (x)





  

J I J I J I Z

Z R R 2 1

EM

13

(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 yi

• And a constraint to keep the orthogonality of the wavefunctions
• And any external additional variables aq (e.g. thermostats, stress,

etc.)

Car-Parrinello Molecular Dynamics

14

15

Car-Parrinello Molecular Dynamics

• Solve the related Euler-Lagrange EOM

* CP * CP

 

i i

dt d y   y   L L 

CP CP

     

I I

dt d R R L L 

CP CP

     

q q

dt d   L L 

Electrons Ions external additional variables

• Solve the related

Euler-Lagrange EOM

Temperature, pressure, etc.

16

Car-Parrinello Molecular Dynamics

17

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

x x x x x x

i i i i i i

y  y    y  y            L L

I I i

M

I

R p

R





  

CP

L

q q q CP q

           L

18

Car-Parrinello Hamiltonian for QM

• so that the QM Hamiltonian reads

 

q I tot q q q i i i I I I

E x d M H        }, { , 2 ) ( ) ( 2

2 * 3 2 CP

R x x p    

  

 

 

 

 

ij j i ij ij

x d  y y  ) ( ) (

* 3

x x

and the CPMD equations of motion will be given by the corresponding Hamilton equations for each set of variables and momenta.

19

BO surface

CP trajectory BO trajectory The difference between the CP trajectories RI

CP(t) and the

Born-Oppenheimer (BO) ones RI

BO(t) is bound by

| RI

CP(t) - RI BO(t)| < C 1/2

(C > 0) if

 

2        

HOMO LUMO

See F.A. Bornemann and C. Schütte, Numerische Mathematik vol.78,

• N. 3, p. 359-376 (1998)

CPMD dynamics: a little warning

20

Practical implementation

• To implement the CP-EOM numerically, the KS orbitals are

generally expanded in plane waves

yi(x) = SG ci(G) eiGx

• G are the reciprocal space vectors. The Hilbert space spanned

by PWs is truncated to a suitable cut-off Ecut such that G2/2 < Ecut

• 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

21

R space a G space E2

cut > E1 cut

Plane wave expansion: yi(x) = SG ci(G) eiGx

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 Gcut

2/2 < Ecut

E1

cut

22

Practical implementation

• Verlet’s algorithm on EOM gives

(m/Dt2) ・[ |yi(t+Dt)> + |yi(t-Dt)> - 2 |yi(t)>] = -(H CP+L) |yi(t)>

• r, in G-space,
• The ionic degrees of freedom RI(t) are updated at a rate (speed)

Dt while the electronic degrees of freedom |yi(t)> are updated at a rate Dt/m1/2 (Dt ~ 5 a.u., m ~ 500 a.u.), hence they are much slower (adiabatic) with respect to the ions.

• However this Dt is driven by QM and is much shorter than  t in

MM simulations

 

 

       D   D  D

 j j ij i CP i i i

c c H t c t t c t t c t ) ( ) ( ) , ( 2 ) , ( ) , (

2

G G G G G G G

G



23

Typical form of the QM-MM Hamiltonian:

Basically just the electrostatic interaction between the two subsystems, i.e.

 

  

  

   

MM 1 QM 1 QM 3 1 int

| | ) ( } , }{ { ), (

I J J I J I I MM I I I I J

Z q x d q q H R R R x x R R x  

r(x) RI with

charge qI

QM electron-MM atom Particle-mesh interaction: We need good charge models for the MM part and a lot of computation QM atom-MM atom

24

Total hamiltonian: scaling in a plane wave basis set

Pure QM interaction Nel x NG x NG with yj(x) = SG cj(G) eiGx QM/MM interface MM x NG Major cost = non-bonding interactions* MM x MM

 

MM MM QM int QM tot

] } { , } { , [ H H H H

I J

   S R R 

*Ewald summations (O(MM log MM)) or spherical cut-off techniques

(O(MM)) can reduce the MM computational cost.

25

Total Hamiltonian: scaling in a localized basis set

Pure QM interaction Nel x Nbasis-set x Nbasis-set QM/MM interface MM x Nbasis-set Major cost = non-bonding interactions MM x MM

 

MM MM QM int QM tot

] } { , } { }, [{ H H H H

I J i

   S R R y







el

N k I k k i i

c

1

}) { ; ( ) ( R x x f y

 

I k k z k y k x k I k

r r r r N

z y x

R x r R x      

2

exp }) { , (  f

e.g. in the case of Gaussian basis-set

26

…and now the forces: CPMD as the QM driver

• the inner part -
• Euler-Lagrange EOM for electrons, ions & Co. inside the QM

region + interaction with MM region

+ ???? + ????

New force components ?

• QM atoms interact in the usual DFT-based

way inside the QM box and also with the

• uter MM atoms (warning: No periodic

boundary conditions here !)

• Electrons interact in the usual DFT-based

way (Kohn-Sham potential) but also with the outer MM electrostatic charges (again, no PBC allowed here)

27

28

MM interaction: forces on atoms

   

???? } { }, {   

I MM I I MM

V q H

I

R R

R

Remember that computing this gradient means evaluating the derivatives with respect to each atom position RI of all the five components of V MM:

• stretching
• bending
• torsion
• Coulomb electrostatics
• van der Waals – Lennard-Jones
• + contribution from the QM atoms !

Forces acting on the MM atoms:

29

QM/MM electrostatic interaction (1): atoms

  



 

MM 1 QM QM MM int

} , { , } , {

I J J I J I J J I I

Z q Z q H R R R R

MM

• QM

where MM = number of classical atoms, QM = number of quantum atoms, qI = MM charge, ZJ = QM charge. We have then two forces components: From MM to QM From QM to MM

 

) 1 (

int MM 1 3 int MM QM J I J I J I I J J

q Z H F R R R R R       



 

 

) 2 (

int QM 1 3 int MM QM I J J I J I J I I

Z q H F R R R R R       



 

30

QM/MM electrostatic interaction (2): electrons

 

 

 

  QM 3 MM 1 int

) ( } , { ), (

I I I I I e

x d q q H R x x R x

MM

 

Functional form (MM = number of classical atoms): Potential acting on the QM wave functions yi(x): Forces acting on the MM charged atoms:

 

x R x x

MM int MM 1 int

) ( V q H

I I I e

   

 

  

 

) 3 ( ) (

int 3 QM 3 int I I I I I e

x d q H F R x R x x R

MM

       







Expensive if MM is large and/or QM is too accurate !

Equations of motion: what we have in the code(s)

31

 

) 3 ( ) 2 (

int int I I I MM I I

V M

I

F F R R

R

     

 

x x

i

V y ) (

int



MM QM

) 1 (

int J DFT J J

E M

J

F R

R

    

32

Car-Parrinello & Born-Oppenheimer MD their respective (dis)advantages - I

CPMD BO Conservation of constants of motion Good Convergence dependent Electronic optimization Not needed Needed Hamiltonian diagonalization Not needed Needed Integration step Small Large Minimum of the BO surface Approximate Exact

• CP-like yi propagation Stability
• Large Dt with no SCF cycle Efficiency
• Dynamics ~ true BO Accuracy
• BO deviation under control Min. Error

33

Car-Parrinello & Born-Oppenheimer MD their respective (dis)advantages - II

BO:

• Diagonalization / minimization of EDFT are required in BO
• Hellmann-Feynman forces are just one component, Pulay

forces exist if local basis set are used

• Residual forces components appear unless strong convergence

criteria are imposed. Dr = r  rBO CP:

• Only Hellmann-Feynman forces (at least in PWs)
• No residual SCF forces because of the fully Lagrangean

(variational) equations of motion

BO H BO BO xc extra

I

V V x d r           r D  r D  r  r  

R

F ] [ ] [

3