Modelowanie Ab-initio Modelowanie Ab-initio

(1)

Modelowanie molekularne

metodami chemii kwantowej

Dr hab. Artur Michalak Zakład Chemii Teoretycznej

Wydział Chemii UJ

http://www.chemia.uj.edu.pl/~michalak/mmod2007/

Wykład 13

Ck08

(2)

• Podstawowe idee i metody chemii kwantowej:

Funkcja falowa, gęstość elektronowa; równanie Schrodingera; Teoria Funkcjonałów Gęstości (DFT);

przyblienie Borna-Oppenheimera, zasada wariacyjna w mechanice kwantowej i w DFT, przyblienie

jednoelektronowe; metoda HF; korelacja elektronowa; metody korelacyjne oparte na funkcji falowej; metoda Kohna-Shama

• Dane do obliczeń kwantowo-chemicznych; GAMESS:

Geometria czasteczki; macierz Z; bazy funkcyjne

w obliczeniach ab initio ; input/output programu GAMESS

• Struktura geometryczna układów molekularnych:

Optymalizacja geometrii; optymalizacja z wiazami; analiza konformacyjna; problem minimum globalnego

• Struktura elektronowa układów molekularnych:

Orbitale molekularne, orbitale KS; wiazanie chemiczne; gęstość rónicowa; orbitale zlokalizowane; analiza populacyjna; analiza rzędów wiązań

• Analiza wibracyjna; Wielkości termodynamiczne;

Reaktywność chemiczna:

Analiza wibracyjna; wielkosci termodynamiczne; modelowanie reakcji chemicznych;

optymalizacja geometrii stanu przejściowego, IRC; indeksy reaktywności chemicznej, molekularny potencjał elektrostatyczny, funkcja Fukui’ego i teoria orbitali granicznych;

jedno- i dwu-reagentowe indeksy reaktywności

• Inne zagadnienia:

Dynamika molekularna (MD); Metody hybrydowe QM/MM; modelowanie wielkich układów; efety rozpuszczalnika; modelowanie w katalizie homo- i heterogenicznej; oddziaływania międzycząsteczkowe, i. in.

(3)

Modelowanie Ab-initio Modelowanie Ab-initio

PES

Podejście Statyczne

elektrony – kwant.; jądra zamrożone

Dynamika kwantowa

elektrony i jądra – mech. kwantowa

(4)

Dynamika Molekularna

elektrony – mech.kwant.; jadra – mech. klas.

Dynamika kwantowa

elektrony i jądra – mech. kwantowa

Modelowanie Ab-initio

(5)

Dynamika molekularna (MD) Dynamika molekularna (MD)

Ab Initio MD

potencjał ab initio wyznaczany w trakcie symulacji („on-the-fly”)

Model-Potential Classical MD

eg. MM-MD

Classical Trajectory Calculations

pre-determined PES

≠

≠ ≠

≠ ≠ ≠ ≠ ≠

(6)

Dynamika molekularna (MD) Dynamika molekularna (MD)

Ab Initio MD

potencjał ab initio wyznaczany w trakcie symulacji („on-the-fly”)

Model-Potential Classical MD

eg. MM-MD

Classical Trajectory Calculations

pre-determined PES

≠

≠ ≠

≠

first-principle MD

≠ ≠ ≠ ≠

quantum-chemical MD on-the-fly MD

direct MD

potential-free MD Hellmann-Feynman MD

Car-Parinello MD

(7)

From: Marx, D.; Hutter, J. in “Modern Methods and Algorithms of Quantum Chemistry”, Grotrndorst, J. (Ed.), NIC Series, Vol. 1, John von Neumann Institute of Computing: Jülich, 2000, pp. 301-449.

Dynamika molekularna (MD)

(8)

Mechanika kwantowa

„kwantowe” jądra i elektrony

Dynamika molekularna

(„kwantowe” elektrony; „klasyczne” jądra) Lata 30 XX wieku:

Erhenfest,

Born-Oppenheimer

Dynamika molekularna (MD)

(9)

Mechanika kwantowa

„kwantowe” jądra i elektrony

Dynamika molekularna

(„kwantowe” elektrony; „klasyczne” jądra) Lata 30 XX wieku:

Erhenfest,

Born-Oppenheimer

Dynamika molekularna (MD) Dynamika molekularna (MD)

Erhenfest MD: Born-Opprnheimer MD:

(10)

Dynamika molekularna (MD) Dynamika molekularna (MD)

Born-Opprnheimer MD:

krok czasowy (‘timestep’):

rzędu 100 a.u.

(11)

Car-Parinello Molecular Dynamics (CP-MD) Car-Parinello Molecular Dynamics (CP-MD)

Równania ruchu dla jąder i f. falowej

Car, R.; Parinello, M. Phys. Rev.Lett. 1985, 55, 2471.

Lagrangian:

s constraint )

; 2 (

1 2

) 1 ,

( = + ² + ₀ Ψ +

= ^L ^R ^R

∑ ∑

^M^R ^E ^R

L

i

i i

α

ψ α

ψ

µ _& _& ^&

&

en. kinet.

ruchu jąder

Potencjał dla ruchu jąder

} s constraint {

)

; (

} s constraint {

)

; (

0 0

E i E R R

M

i i

a

ψ δψ

ψ δ µ

α α

α

∂ + ∂ Ψ

−

=

∂ + ∂ Ψ

−∇

=

R R

&

(12)

Lagrangian:

s constraint )

; 2 (

1 2

) 1 ,

( = + ² + ₀ Ψ +

= ^L ^R ^R

∑ ∑

^M^R ^E ^R

L

i

i i

α

ψ α

ψ

µ _& _& ^&

&

en. kinet.

f. falowej

} s constraint {

)

; (

} s constraint {

)

; (

0 0

E i E R R

M

i i

a

ψ δψ

ψ δ µ

α α

α

∂ + ∂ Ψ

−

=

∂ + ∂ Ψ

−∇

=

R R

&

en. kinet.

ruchu jąder

(13)

Lagrangian:

s constraint )

; 2 (

1 2

) 1 ,

( = + ² + ₀ Ψ +

= ^L ^R ^R

∑ ∑

^M^R ^E ^R

L

i

i i

α

ψ α

ψ

µ _& _& ^&

&

en. kinet.

f. falowej

} s constraint {

)

; (

} s constraint {

)

; (

0 0

E i E R R

M

i i

a

ψ δψ

ψ δ µ

α α

α

∂ + ∂ Ψ

−

=

∂ + ∂ Ψ

−∇

=

R R

&

en. kinet.

ruchu jąder Uwaga! Nie mylić z T_el

(14)

Lagrangian:

s constraint )

; 2 (

1 2

) 1 ,

( = + ² + ₀ Ψ +

= ^L ^R ^R

∑ ∑

^M^R ^E ^R

L

i

i i

α

ψ α

ψ

µ _& _& ^&

&

en. kinet.

f. falowej

} s constraint {

)

; (

} s constraint {

)

; (

0 0

E i E R R

M

i i

a

ψ δψ

ψ δ µ

α α

α

∂ + ∂ Ψ

−

=

∂ + ∂ Ψ

−∇

=

R R

&

en. kinet.

ruchu jąder

Fictitious mass, fictitious kinetic energy wave-function mass, kinetic energy

(15)

CP MD: BO MD:

(16)

Car-Parinello Molecular Dynamics Car-Parinello Molecular Dynamics

CP MD: BO MD:

(17)

rzędu 100 a.u.

rzędu 10 a.u.

(18)

Zasada zachowania energii dla CP-MD:

s constraint )

; 2 (

1 2

) 1 ,

( = + ² + ₀ Ψ +

= ^L ^R ^R

∑ ∑

^M^R ^E ^R

L

i

i i

α

ψ α

ψ

µ _& _& ^&

&

(19)

Porównanie BO-MD i CP-MD. Etylen, symulacje z tej samej geometrii startowej;

program CPMD program; pseudopotencjał Troulliera-Martinsa ; timestep 4 a.u., masa ff 400 a.u., energia cut-off 70 Ry, kom.el. 12 A x 12 A x12 A).

(20)

Siły w dynamice molekularnej Siły w dynamice molekularnej

Analityczne gradienty:

(21)

Twierdzenie Hellmanna-Feynmana :

Prawdziwe dla : - dokładnej funkcji falowej,

- funkcji wariacyjnej w kompletnej bazie

Analityczne gradienty:

(22)

Sily Hellmanna-Feynmana

Incomplete-basis-set force

„sily Pulay”

Poprawka ‘Non-self-consistency’

(23)

Bazy funkcyjne Bazy funkcyjne

Typu Slatera:

Typu Gaussa:

Fale płaskie (plane waves):

(24)

Typu Slatera:

Typu Gaussa:

Centrowane na

atomach

‘origin-less’:

nie ma sil Pulay’a

(25)

Typu Slatera:

Typu Gaussa:

Centrowane na

atomach

‘origin-less’:

- nie ma sil Pulay’a - Periodyczność (!)

(26)

Fale plaskie: energia obcięcia ( cut-off energy) Fale plaskie: energia obcięcia ( cut-off energy)

Orbitale KS :

Gęstość:

Obcięcie

Liczba fal plaskich dla zalożonej energii E_cut

(27)

Fale płaskie Fale płaskie

Problemy z obszarem w pobliżu rdzenia

• Pseudopotenciały

• ‘Augmentation schemes’

I II Podejście I

‘muffin-tin’:

I sfery atomowe II obszary

międzyatomowe

(28)

Etylen

Butadien

CPMD program; Troullier-Martins pseudopotentials, time step of 4 a.u., fictitious mass 400 a.u., unit cell 12 A x 12 A x12 A

(29)

Etylen

Butadien

CPMD program; Troullier-Martins pseudopotentials, time step of 4 a.u., fictitious mass 400 a.u., unit cell 12 A x 12 A x12 A

(30)

Zbieżność i wymagana E_cut-offzależy od metody reprezentacji rdzenia

pseudopotencjały Troulliera-Martinsa : 60-100 Ry pseudopotencjały Vanderbilta (‘ultgasoft’): 20-40 Ry Pseudopotencjały Goedeckera: 100-200 Ry

metoda PAW 20-40 Ry

(31)

Termostaty Termostaty

Nose-Hoover chain thermostat:

Jądrowe równania ruchu:

Elektronowe równania ruchu:

(32)

Symulacja MD – typowe elementy

(33)

• dynamika z wiązami,

(jądrowymi lub elektronowymi)

• ‘umbrella sampling’; ‘bias potentials’

Modelowanie reakcji chemicznych Modelowanie reakcji chemicznych

E

RC

Skala czasowa symulacji ab initio MD – max. kilka ns:

brak szans na spontaniczne zajście reakcji

(34)

TS

min.

• Założona współrzędna reakcji

• dynamika z więzami dla punktów na ściezce

• całkowanie termodynamiczne (thermodynamic integration) Energia swobodna reakcji chemicznych

Energia swobodna reakcji chemicznych

∆A = F

_i _λ

∆ λ

_i

i n_po_{int s}

∑

(35)

TS

min.

• Założona współrzędna reakcji

• dynamika z więzami dla punktów na ściezce

• całkowanie termodynamiczne (thermodynamic integration) Energia swobodna reakcji chemicznych

Energia swobodna reakcji chemicznych

Każdy punkt wymaga:

1) Wstępnego uzbieżnienia f.falowej 2) podgrzania do zalożonej T

3) równowagowania termicznego

(36)

TS

min.

• współrzędna reakcji λ zmieniana w sposób ciągły

Metoda powolnego wzrostu (slow-growth)

(37)

TS

min.

• współrzędna reakcji λ zmieniana w sposób ciągły

Metoda powolnego wzrostu (slow-growth) Metoda powolnego wzrostu (slow-growth)

Tylko pierwszy punkt wymaga:

1) wstępnego uzbieżnienia f.falowej 2) podgrzania do zalożonej T

3) równowagowania termicznego

(38)

Typowy problem – histereza w profilach energii swobodnej

∆A

RC

forward sampling backward sampling

Metoda powolnego wzrostu (slow-growth)

(39)

Wybór współrzędnej reakcji Wybór współrzędnej reakcji

TS

min.

Przekrój w kierunku prostopadłym do ścieżki:

zmienia się od punktu do punktu

Kierunek prostopadły do scieżki

(40)

TS

min.

Przekrój w kierunku prostopadłym do ścieżki:

łagodne przejście od punktu do punktu

Kierunek prostopadły do scieżki

Wybór współrzędnej reakcji Wybór współrzędnej reakcji

(41)

MD wzdłuż IRP MD wzdłuż IRP

A. Michalak, T. Ziegler „First-principle Molecular Dynamics along Intrinsic Reaction Paths”, J. Phys Chem. A 105, 2001, 4333-4343.

• wstępne wyznaczenie IRC

• MD z więzem zamrażającym ruch wzdłuż IRC

(42)

HCN TS CNH

IRP:

HCN →→→ CNH→ HCN →→→ CNH→

(43)

HCN →→→→ CNH HCN →→→→ CNH

MD wzdłuż IRP (300K)

MD z więzem

R_NH -R_CH= const.

IRC (T=0K)

(44)

Trajektoria wodoru

HCN →→→ CNH→ HCN →→→ CNH→

MD z więzem

R_NH -R_CH= const.

(45)

HCN →→→ CNH→ HCN →→→ CNH→

Trajektoria wodoru

MD z więzem

R_NH -R_CH= const.

(46)

HCN →→→ CNH→ HCN →→→ CNH→

Trajektoria wodoru

MD z więzem

R_NH -R_CH= const.

(47)

cyclobutene TS gauche-butadiene Conrotatory ring opening of cyclobutene

Conrotatory ring opening of cyclobutene

(48)

cyclobutene TS gauche-butadiene

IRP:

Conrotatory ring opening of cyclobutene Conrotatory ring opening of cyclobutene

(49)

(50)

Cl^- + CH₃Cl TS ^Cl-CH3 + Cl^- Prototype SN2 reaction : Cl^- + CH₃Cl →→→→ CH₃Cl + Cl^- Prototype SN2 reaction : Cl^- + CH₃Cl →→→→ CH₃Cl + Cl^-

(51)

Cl^- + CH₃Cl TS ^Cl-CH3 + Cl^-

IRP ( T = 0 K ):

Prototype SN2 reaction : Cl^- + CH₃Cl →→→→ CH₃Cl + Cl^- Prototype SN2 reaction : Cl^- + CH₃Cl →→→→ CH₃Cl + Cl^-

(52)

0 1 2 3 s[amu-1 bohr]

0 30 60 90 120 150 180 Angle

Cl1-C-Cl2 IRC Cl1-C-H

0 1 2 3 s [amu-1 bohr]

2 3 4 5 R [A]

IRC

C-Cl2 Cl1-C Cl1 - Cl2

0 1 2 3 s[amu-1 bohr]

-5 -4 -3 -2 -1 0 E [kcal/mol]

∆ _IRC

∆ G

TS E

vdW complex

Prototype SN2 reaction : Cl^- + CH₃Cl →→→→ CH₃Cl + Cl^- Prototype SN2 reaction : Cl^- + CH₃Cl →→→→ CH₃Cl + Cl^-

(53)

Cl -CH₂-CH=CH₂ TS ^CH2=CH-CH₂-Cl

IRP (TS → R):

CH₂=CH-CH₂Cl isomerization CH₂=CH-CH₂Cl isomerization

(54)

Cl-CH₂-CH=CH₂ Cl-CH₂-CH=CH₂

(55)

TS

conf. 2 (gauche) conf. 1 (cis)

(56)

TS

IRP (T = 0 K)

(57)

TS

IRP (T = 0 K ) T = 300 K

(58)

0 2 4 6 8 s [amu-1 bohr]

-40 -30 -20 -10 0 E [kcal/mol]

E G

∆ IRC

∆

0 2 4 6 8 s [amu-1 bohr]

1 2 3 4 R [A]

IRC Cl-C3 Cl-C1 C1-C2 C2-C3

0 2 4 6 8 s [amu-1 bohr]

0 30 60 90 120 Angle

Cl-C1-C2-C3 Cl-C1-C2 IRC

C1-C2-C3

TS

cis-

CH₂=CH-CH₂Cl isomerization CH₂=CH-CH₂Cl isomerization

(59)

Cdn.