CSC 2002 Conference, Vancouver, BC, 2 June 2002

Ab Initio

Molecular Dynamics with a Continuum Solvation Model

Hans Martin Senn, Tom Ziegler

University of Calgary

Outline

 Background 

Car –Parrinello ab initio molecular dynamics



The Projector-Augmented Wave (PAW) method



Continuum solvation / The Conductor-like Screening Model (COSMO)

 Continuum solvation within

ab initio

MD 

Surface with smooth analytic derivatives



Surface charges as dynamic variables

 Tests and first applications 

Hydration energies



Conformers of glycine



CO insertion into Ir –CH 3

 

Ab Initio

Molecular Dynamics



Ab initio

molecular dynamics

Nuclei

Classical mechanics

Electrons

Density-functional theory R. Car, M. Parrinello, Phys. Rev. Lett. 1985 ,

55,

2471

Ab initio molecular dynamics

Forces on nuclei derived from instantaneous electronic structure  Car –Parrinello scheme for

ab initio

MD 

Fictitious dynamics of the wavefunctions



Simultaneous treatment of atomic and electronic structure

+ 

T



V

 1 2 

J M J

Ý

J

 

i m c

Ý

i

Ý

i



E

 {

R I

},{

c i

}   

i

,

j

 DFT total energy

c

nuclear kin. energy fictitious kin. energy of the wavefunctions

i c j



d ij



L

orthonormality constraint

ji M J R J m c

Ý

j

  

E

 {

R I

},{

c



R J i

}   

dE

 {

R I

},{

c d c j i

}   

i c i L ij



Friction dynamics for minimizations



The Projector-Augmented Wave (PAW) Method

 Basis set: Augmented plane waves P.E. Blöchl, Phys. Rev. B 1994 ,

50,

17953 

Spatial decomposition of the wavefunction according to its shape

 Smooth plane-wave expansion far from nuclei  Rapidly oscillating augmentation functions around the nuclei

c

  

k f k

˜ (

r

)

f f

˜

k k

( (

r r

) )   

G

x k c

(

G

)e i

G r

Y l

,

m

˜

k

(

r

) ˜

k

 

k

˜

k

˜

k

Plane-wave expansion Numerical , atomic -orbital-like functions Smoothened atomic orbitals Projector functions 

Integrations

for total energy and matrix elements decompose accordingly  Smooth parts in Fourier representation  One-centre parts on radial grids in spherical harmonics expansions  All-electron method 

Frozen-core approximation



Continuum Solvation Models

  Implicit electrostatic solvation  

“Solvent” = homogeneous, isotropic dielectric medium Fully characterized by relative permittivity (

 r

)

Electrostatic interactions 1.

Solute polarizes the dielectric (  reaction field) 2.

3.

4.

Solute charge density interacts with polarized medium Solute electronic and atomic structure adapt (“back-polarization”) Repeat to self-consistency!

 Electrostatic free energy 

For a linear dielectric:

G

diel 

E

int 

E

pol  1 2

E

int

E

pol > 0

E

int < 0 ∆

E

en > 0

 

The Solvation Energy

 Non-electrostatic contributions 

Cavity formation



Dispersion and exchange-repulsion



Approximately proportional to cavity surface area

np 

b



gA

 Free energy of solvation

G

sol  sol

E

en 

G

diel 

G

np  solv

G



G

sol 

E

g en 

G

diel 

G

np 

DE

en 

Neglect changes in thermal motion upon solvation

   

The Conductor-Like Screening Model (COSMO)

 Main approximations in COSMO

1.

Dielectric with infinite permittivity (



r

 ∞)  No volume polarization; polarization expressed as surface charge density

only

 Vanishing potential at the cavity boundary

G

 diel    

V d

3

r



S d

2

s

r

(

r

)

s

r



s

(

s

)  1 2 

S d

2

s



S d

2

s



s

(

s

)

s

(

s

 )

s

f



(

s

)



dG

 diel

ds

s

 

V d

3

r

r

(

r

)

r



s

 

S d

2

s



s

(

s



)

s

  

0

2.

Recover true dielectric behaviour by scaling

  Derived for rigid multipoles in spherical cavity

f

=

f

(  r ) = (  r – 1) / (  r +

x

) [

x

= 0, 1/2 for monopole, dipole]

G

COSMO diel 

fG

 diel

s

(

s

) 

fs

(

s

)

3.

Discretize cavity surface and surface charge

 Surface segments carrying point charges

q i



G

COSMO diel    

i q i



V d

3

r r

r

(

r

)



s

i

Energy expression is variational in

q i



1 2

f



i

,

j

(

j



i

)

s

i q i q j



s

j



1.07 4

p

2

f



i q i

2

a i

 

Surface Charges as Dynamic Variables

 Extended Lagrangean including solvation

+



+

CP  

i m q



G

COSMO diel 

{

R

J

},{

q j

}

 

G

np 

{

R

J

}

 

Equations of motion for charges

m q

Ý

 

+



q i



m q

Ý

q



Additional forces acting on the nuclei

 Surface segments and charges implicitly depend on atomic positions

via

the surface construction  Forces obtained as analytic derivatives 

Energy-conserving dynamics requires smooth derivatives

 Segments/charges must not vanish or be created during the simulation  Each segment/charge must be unambiguously assigned to an atom  Conventional schemes for building the cavity are not compatible with this!



Surface Construction

 A surface construction compatible with AIMD 1.

Surface = union of atomic spheres 2.

Each sphere is triangulated 3.

All segments and charges are retained Each segment/charge moves

rigidly

along with its atom 4.

A

switching function

effectively removes charges lying within another sphere from the energy expression COSMO

G

diel  

i q i Q i



V d

3

r r

r

(

r

) 

s

i

    1 2

f



i

,

j

(

j



i

)

q i Q i q

s

i



j Q j

s

j

 1.07

f

Q

i

= 1 if the charge is exposed on the molecular surface Q

i

= 0 if the charge is buried Smooth transition between “on” and “off”

p



i q i

2

Q i

2

a i

  

Switching Functions

 Building the switching function 

Rectangular pulse

h

(

d

,

d

0

)



1



exp

    2

n

 

,

 Centred at

d

0 =

R

solv –

c

 Vanishes smoothly at

d

0

n



N



Atomic switching function

u J

,

i

 

h

 0 

R

J



s

i

,

d

0

J

 for

R

J

for

R

J

 

s

i

s

i



d J

0 

d J

0 h

d

0

R

solv 

Total switching function

Q i

 

J

(

J



I

(

i

))

u J

,

i R

A solv  Charge is smoothly switched off if it lies within any one of the other spheres

R

B solv

n n

= 1 = 10

n

= 24

d



Rounding the Edges…

 Behaviour of hidden charges 

Magnitude of switched-off charges physically undetermined

 Problematic if charge becomes again exposed 

Restoring potential

 Added to the Lagrangean 

E

rest  

i k

(1



Q i

)

q i

2  Modified Coulomb potential at short range 

Charges can collide at seams



Modified Coulomb potential

  True 1/

r

behaviour replaced by linear continuation for |

s

i

–

s

j

| <

R ij

c

R ij

c derived from average of self-interaction potentials

R ij

c pot.

|

s

i

–

s

j

|

 

Representation of the Solute Density

 Decoupling of periodic images 

Plane-wave-based methods always create periodic images

 Artificial in molecular calculations 

Electrostatic decoupling scheme in PAW

 Model density  Superposition of atom-centred Gaussians  Coefficients determined by fit to true density  Reproduces long-range behaviour

ˆ (

r

)

 

I

,

m Q I

,

m g I

,

m

(

r

)

 Surface charges couple to model density 

Not instrumental



Computationally convenient

 3 –5 Gaussians per atom  Fitting procedure relays forces acting on the wave functions due to the charges

E

s– m  

i

,

I

,

m q i Q i Q I

,

m

R

I



s

i

erf



R

I



s

i r I

c ,

m



Parameterization and Validation

 Parameterization 

Free energies of hydration

 Parameters: Solvation radii

R I

solv , non-polar parameters b , g  Fit set: Neutral organic molecules containing C, H, O , N CH 4 MUE: 3.2 kJ/mol H H O H H O O H O H O O O O H O H O H O H 2 N N H O O H O N N H 3 N H O O N H O N H N H 2 O N H 2 N O 2 N H 4.2 kJ/mol 5.2 kJ/mol  Mean unsigned error (MUE) over whole set:

4.1 kJ/mol

 On par with similar DFT/COSMO results

Conformers of Glycine

 H-bonding patterns

A B



Relative stabilities

(kJ/mol)

Gas p hase

Method

PAW

O ther pure DFT Hybrid DFT Expt.

C

0

0 0 0

A

4.1

1.1…4.6

–1.6…–1.1

–5.9

B

9.6

6.8

4.7…5.6

C Z Aque ous ph ase

Method

PAW /COSMO

O ther pure DFT Hybrid DFT Expt.

Z

0

0 0

A

49

36…39 13 30…32

B

55



Hydration energies

(kJ/mol) Method

PAW /COSMO

O ther pure DFT Hybrid DFT Expt.

A

(g) 

A

(aq) ²

–57

–54 –49…–42 hyd

G

A

(g) 

Z

(aq)

–105

–91 –61…–57 –80 

Zwitterion in solution

 Most stable conformer  Not present in the gas phase

Energy Conservation

 Conservation of total energy 

Glycine without and with continuum solvation

  ∆

t

= 7 a.u.,

T

≈ 300 K (no thermostats, no friction) Drift in total energy: solvation off solvation on +2.78   1.43   5

E

h /ps per atom  5

E

h /ps per atom  Energy is conserved

A First Application

 Carbonylation of methanol 

“Monsanto” / ”Cativa” acetic-acid process

  M = Rh, Ir

Via

[MI 3 (CH 3 )(CO) 2 ] –  [MI 3 (COCH 3 )(CO)] – MeOH cat. [MI 2 (CO) 2 ] – , HI, H 2 O CO 

AIMD simulation of CO insertion step

 Solvent: iodomethane;

T

= 300 K   Slow-growth thermodynamic integration along

d

(C Me –C CO ) ∆

A

‡ = 111 (128) kJ/mol ∆

E

‡ = 126 (155) kJ/mol; ∆

S

‡ = 60 (91) J/(K mol) 

Dissociation of I – trans to incipient acyl!

 Only in solution at finite temperature  Enthalpy of Ir –I bond traded for entropy of free (solvated) I – MeCOOH RC

Summary

 COSMO continuum solvation in AIMD 

Surface charges as dynamic variables



Surface construction

 Switching function smoothly disables unexposed charges  Smooth analytic derivatives wrt. atomic positions 

Energy-conserving



Solvation energies

 On par with other DFT/COSMO implementations 

Modelling of finite-temperature and solvation effects

Acknowledgments



Peter E. Blöchl,

Institute of Theoretical Physics, Clausthal University of Technology, Germany 

Peter M. Margl,

Corporate R&D Computing, Modeling and Information Sciences, Dow 

Rochus Schmid,

Institute of Inorganic Chemistry, Munich University of Technology, Germany $

Swiss National Science Foundation

$

NSERC

$

Petroleum Research Fund