Molekulární modelování a simulace
1/35 pch01
Evropský sociální fond þPraha & EU: Investujeme do va¹í budoucnostiÿ Inovace pøedmìtu Poèítaèová chemie je podporována projektem CHEMnote (Inovace bakaláøského studijního programu Chemie { moderní vzdìlávání podpoøené pou¾itím notebookù { CZ.2.17/3.1.00/33248) v rámci Operaèního programu PRAHA { ADAPTABILITA.
Základní prvky modelování
2/35 pch01
? elementární èástice + gravitace: þteorie v¹ehoÿ { temná hmota. . . (známé) elementární èástice: standardní model { atomová jádra. . . jádra + elektrony + fotony: QED { pøesná spektroskopie. . . jádra + elektrony: Schrödingerova rovnice { (praktická) chemie vlastnosti malých molekul, spektra, kinetika, fotochemie. . . Atomy { klasické molekulové modelování∗; kvantová jádra: PI Hrubozrnné (coarse-grained) modely: mezo/nanoskopická ¹kála více atomù = 1 objekt (surfaktant = hlavièka + ocásek, èlánek polymerního øetìzce. . . ) Mikroskopická ¹kála: disperzní systémy, sypké materiály
P
Materiál jako kontinuum: parciální diferenciální rovnice gravitace: prostoroèas multiscale modeling: QM/MM, . . .
∗ pøíp. pomocná centra / vìt¹í skupiny (-CH ) 3
Schrödingerova rovnice aneb teorie ve¹keré chemie
3/35 pch01
V principu (QED) umíme spoèítat bod varu vody na 12 desetinných míst z velièin e, h, me, mO, mH, c. V realitì dosáhneme po mnohaleté práci chybu nìkolik % Pøímo celkem dobøe umíme popsat malé systémy: tvar a energie (malých) molekul, spektra, kinetika K pro reakce v plynné fázi leccos umí fyzika i pro krystaly
(Hyper)plocha potenciální energie
[plot/rcoord.sh]
4/35 pch01
Jádra jsou mnohem tì¾¹í ne¾ elektrony ⇒ elektronové pohyby jsou mnohem rychlej¹í (tzv. Bornova{Oppenheimerova aproximace) potential energy surface (PES) energie jako funkce souøadnic poloh v¹ech atomových jader Epot(~r1, ~r2, . . . , ~rN)
Reakce probíhá cestou nejmen¹ího odporu = pøes sedlový bod (pøesnìji: v jeho blízkosti) =
tranzitní stav
credits: http://www.ucl.ac.uk/ ucecmst/publications.html, http://theory.cm.utexas.edu/henkelman/research/ltd/
Jak získám PES?
5/35 pch01
z kvantových výpoètù (Schrödingerova rovnice +) aproximujeme vzorcem (þsilové poleÿ, þpotenciálÿ, þmodelÿ, . . . ): silové pole: (force eld) Epot = souèet mnoha èlenù, èlen = funkèní tvar + parametry pro rùzné atomy/skupiny kombinace { QM/MM metody (quantum mech./molec. mech.) PES a modelování v chemii
pou¾iju klasickou mechaniku: na statické výpoèty (minimum energie, potenciál v okolí aj.) na výpoèet vývoje systému v èase (molekulová dynamika) na výpoèet termodynamických velièin vzorkováním (Monte Carlo, MD) pou¾iju kvantovou mechaniku (na jádra): metoda dráhového integrálu (PI MC, PI MD) pou¾iju klasickou mechaniku s kvantovými korekcemi kombinace silové pole + klasická mechanika = þmolekulová mechanikaÿ (MM); v u¾¹ím smyslu nezahrnuje MC a MD
6/35 pch01
Modelování v chemii: dìlba práce
kvantová mechanika
jádra elektrony
makroskopické vlastnosti
sta m tis e kin cha tick sim et nik á ula ická a ce te or ie
(kvantové simulace)
-
vá a to ik an an kv ech m
@ @
@
@ @
@
@ @
@
@ R @
molekulový model
7/35 pch01
Pøed r. 1930
sta m tis e kin cha tick et nik á ick a á te or ie
makroskopické vlastnosti
molekulový model
8/35 pch01
Molekulová mechanika { statický pohled
Energie jako funkce souøadnic (hyperplocha potenciální energie, PES).
Aproximujeme funkcí zvanou silové pole (FF). Minimalizace energie (T = 0), þoptimalizace strukturyÿ Re nement { zpøesnìní struktury (z rozptylových experimentù) Biochemie: tvar molekul (klíè + zámek), síly (hydro lní/hydrofobní. . .)
. . . ale co
pohyb?
9/35 pch01
Co je to pohyb?
þSkuteènýÿ pohyb molekul (tekutiny, . . . ) v èase Soubor v¹ech mo¾ných kon gurací (molekul) zprùmìrovaný v èase:
→
Statistická termodynamika se systematicky zabývá výpoètem velièin (bod
varu, a nita ligandu k receptoru) na základì pøedstavy (makro)stavu systému jako þprùmìruÿ v¹ech mo¾ných kon gurací
Molekulové simulace
10/35 pch01
molekulová dynamika (MD) èasový vývoj systému slo¾eného z mnoha molekul pohyb ka¾dého atomu je urèen silami, které na nìj v ka¾dém okam¾iku pùsobí metoda Monte Carlo (MC); pøesnìji Metropolisova metoda a varianty posloupnost kon gurací systému se generuje pomocí náhodných èísel provedeme náhodný pohyb molekuly a rozhodneme se, zda jej pøijmeme { tak, aby pravdìpodobnosti výskytu kon gurací molekul byly stejné jako v realitì kinetické Monte Carlo simulovaný dìj je rozdìlen na elementární události (napø. adsorpce atomu na rostoucím krystalu, reakce na katalyzátoru) událost, ke které dojde, vybíráme podle známé pravdìpodobnosti kvantové simulace { MD, MC
Co simulujeme
Kapaliny: vliv struktury na vlastnosti (anomálie vody), roztoky fázové rovnováhy, rozpustnost povrchy a rozhraní, surfaktanty Pevné látky: struktura krystalù, materiály (poruchy) adsorpce (zeolity) Biochemie: proteiny, nukleové kyseliny, iontové kanály, lipidické membrány Nanoobjekty: micely, polymery, samoskladba (coarse-grained modely, møí¾ky) Podobnými metodami lze studovat: sypké materiály, rùzné minimalizace (MC), ¹íøení epidemií atd.
11/35 pch01
[uvodsim/blend.sh]
Optimalizace struktury (molekulová mechanika)
←
twist (skew) boat zkøí¾ená vanièka experiment: 45 kJ/mol model: 53 kJ/mol
→
chair ¾idlièka experiment: 28 kJ/mol model: 26 kJ/mol
12/35 pch01
[cd /home/jiri/tex/talks/letniskola; cytox.sh]
Pøíklad { elektrosprej Cytochromu C
Yi Mao, J. Woenckhaus, J. Kolafa, M.A. Ratner, M.F. Jarrold Elektrosprej: rozpra¹ování nabitých èástic Mìøí se úèinný prùøez
13/35 pch01
Pøíklad { elektrosprej Cytochromu C
14/35 pch01
Pøíklad { voda
10000 molekul 300 K periodické ve smìrech x,y adhezivní podlo¾ka neadhezivní poklièka
[cd /home/jiri/tex/talks/letniskola; showvid.sh]
15/35 pch01
SIMOLANT
[cd uvodsim; simolant -g +100+50 -T.1]
Vlastnosti: 2D þatomyÿ (potenciál Lennard-Jonesova typu) odpudivé/pøita¾livé stìny, gravitace MC i MD konstantní energie i termostat Jevy: kondenzace plynu zmrznutí kapky poruchy krystalu kapilární deprese a elevace plyn v gravitaèním poli nukleace Chcete si sami nainstalovat? Staèí Google SIMOLANT. . .
16/35 pch01
[show/janus.sh]
Self-assembly (samoskladba)
17/35 pch01
Supramolekulární chemie: skládání molekul pomocí (zpravidla) nekovalentních sil (van der Waals, vodíkové vazby) do strukturovaných celkù Ukázka: dvoufunkèní èástice v roztoku ≈ Janus particles
Janus
Janus
Iapetus credit: wikipedie, www.nasa.gov/mission pages/cassini
Ukázka: + ètyøfunkèní èástice
credit: Atwood et al., Science 309, 2037 (2005)
Jak dostat minimum energie
Na10Cl10 rychlé chlazení pomalé chlazení
↓
→
[uvodsim/min.sh]
18/35 pch01
[uvodsim/salesman.sh 100]
Simulované ¾íhání (simulated annealing)
19/35 pch01
Hledáme globální minimum funkce (þenergieÿ) s mnoha lokálními minimy Zaèneme nìjakou ¹patnou kon gurací (napø. náhodnou) Navrhneme vhodné zmìny kon gurace Ai → Aj Aplikujeme Metropolisovu metodu za sni¾ující se þteplotyÿ T
Pøíklad: Problém obchodního cestujícího (traveling salesman) 50 mìst náhodnì ve ètverci 1 × 1 Kon gurace = poøadí mìst þenergieÿ = délka cesty Zmìna kon gurace = zámìna 2 náhodnì zvolených mìst
T = 0.1 → 10−5 l = 5.37
T =0 l = 7.93
Genetické algoritmy
+
Hledáme maximum funkce zvané zde þ tnessÿ 0 kon gurace → jedinec → genom = chromozom = seznam 1 alel (þsloupec èíselÿ) 2 3 Zvolíme (napø. náhodnou) poèáteèní populaci 4 Generování následující generace: 5 { vyhodíme nejhor¹í èást populace 6 { èást jedincù zkopírujeme s mutací 7 { (nejvìt¹í) èást jedincù získáme køí¾ením (crossing over) ze dvou rodièù 8 Èísla se kódují Grayovým kódem (po sobì jdoucí pøirozená èísla 9 10 se li¹í v jenom bitu) 11 Aplikace: 12 { logistika, ekonomie, øízení prùmyslových procesù 13 { biochemie { protein folding aj. 14 { elektronika { návrh obvodù, tvar antény { vývoj algoritmù 15
20/35 pch01
0000 0001 0011 0010 0110 0111 0101 0100 1100 1101 1111 1110 1010 1011 1001 1000
[. ./simul/rayleigh/show.sh]
(Plateauova-)Rayleighova nestabilita
Èúrek vody se rozpadá na kapky. Nestabilita pro kr < 1 (pro poruchu ∝ sin(kz)), max. nestabilita pro kr = ln 2.
21/35 pch01
Nukleace pøi supersonické expanzi
[show/supexp.sh]
22/35 pch01
Vodní pára o tlaku cca 5 bar se pou¹tí velmi úzkým otvorem pøes trysku do vakua a adiabaticky se ochlazuje pod bod mrazu. Lze tak studovat napø. chem. reakce ve stratosféøe.
credit: M. Fárník
Otázka: Jaký je tvar, velikost a struktura klastrù ledu?
Tání nanoèástic
kroupa z 600 molekul vody (led Ih) ohøívání èas simulace = 5 ns tento model vody taje pøi 253 K nanoèástice taje pøi ni¾¹í teplotì
[show/kroupa.sh]
23/35 pch01
24/35 pch01
Síly mezi molekulami
Londonovy (disperzní) síly pro vìt¹í vzdálenosti: model uktuující dipól{ 2 E / (kJ mol-1)
{indukovaný dipól elst. pole E ∝ 1/r3 indukovaný dipól µind ∝ E energie u(r) ∝ µE ∝ 1/r6 (v¾dy záporná) Odpuzování na krat¹ích vzdálenostech: u(r) ∝ e−const r Celkem: u(r) = Ae−Br − C/r6 Aproximace odpudivých sil:
Ar...Ar
0
Ae−Br → A 0/r12
Lennard-Jonesùv potenciál : u(r) = 4
σ 12 r
−
σ 6 r
Tyto síly jsou souèástí interakcí mezi
0
0.2
0.4
0.6
r / nm
v¹emi atomy a molekulami
0.8
1
25/35 pch01
Elektrické síly
náboj{náboj (ionty) 1
qi qj U= 4π0 rij
parciální náboje: takové náboje na atomových jádrech, aby se to chovalo stejnì jako skuteèné nábojové rozlo¾ení dipólový moment ~µ =
X
qi~ri
i
polarizovatelnost (el. pole indukuje dipól) ~µind = α~E
26/35 pch01
Elektrické síly
náboj{náboj (ionty) 1
qi qj U= 4π0 rij
parciální náboje: takové náboje na atomových jádrech, aby se to chovalo stejnì jako skuteèné nábojové rozlo¾ení dipólový moment ~µ =
X
qi~ri
i
polarizovatelnost (el. pole indukuje dipól) ~µind = α~E
Silové pole
27/35 pch01
Molekulový model èi silové pole (force eld) je matematický zápis energie molekuly nebo souboru molekul jako funkce souøadnic atomù, ~ri, i = 1, . . . , N. malé: tuhá◦ tìlesa { rotace (voda 25 C: vibruje 0.05 % molekul) velké: mnoho èlenù vazebné síly: vibrace vazeb (1{2): U = K(r − r0)2 lze nahradit pevnou vazbou vibrace úhlù 1 4 P torze (1{4) a torzní potenciál: n Kn cos(nφ) φ 3 2 \improper torsion" (dr¾í >C=O v rovinì) nevazebné síly (èást. 1{4, 1{dále): Lennard-Jones, náboje A v¹echny pøíspìvky seèteme = aproximace párové aditivity Noo, ideálnì pøesná není, øeknìme na 90 %
Konstrukce silových polí
28/35 pch01
geometrie: spektroskopie, difrakce, kvantové výpoèty vazebné síly: kvantové výpoèty, spektroskopie Lennard-Jones σ: experimentální hustota, struktura (difrakce) Lennard-Jones : výparná entalpie parciální náboje: { dipólové momenty: spektroskopie, permitivita { kvantové výpoèty (Mulliken, CHELPG = CHarges from Electrostatic Potentials using a Grid based method) a/nebo: struktura klastrù (z kvantových výpoètù)
29/35 pch01
Molekulová dynamika
tuhé koule ap. { nárazy þklasickáÿ MD { integrace pohybových rovnic Brownovská (stochastická) dynamika { MD + náhodné síly Síly: N) ∂U(~ r ~fi = − ∂~ri
Pøíklad: U=
i = 1, . . . , N X
u(rij)
i<j
⇒ ~fi =
N X j=1 j6=i
~fji ≡
du(rji) ~rji − drji rji j=1
N X
j6=i
30/35 pch01
Newtonovy rovnice ~fi ~ ri = , mi
Verletova metoda Taylor
i = 1, . . . , N
⇒ ~ri(t − h) − 2~ri(t) + ~ri(t + h) 2 + O(h ) ~ri(t) = 2 h
Verletova metoda: ~ri(t + h) = 2~ri(t) − ~ri(t − h) + h
~
2 fi (t)
mi
Poèáteèní podmínky: 2~ h fi(t0) ~ri(t0 − h) = ~ri(t0) − h~r_ i(t0) + + O(h3) 2 mi
31/35 pch01
Metoda leap-frog
rychlost = dráha (zmìna polohy) za jednotku èasu (h)
zrychlení = zmìna rychlosti za jednotku èasu ~f ~v(t + h/2) − ~v(t − h/2) a ~ (t) =
h
=
m
⇒ ~v(t + h/2) = ~v(t − h/2) + a ~ (t)h ~r(t + h) = ~r(t) + ~v(t + h/2)h
opakujeme s
t := t + h
Verlet a leap-frog jsou ekvivalentní
credit: http://www.anagrammer.com/scrabble/leapfrog
~r(t + h) − ~r(t) ~v(t + h/2) = h
Pøíklad: dráha planety
WWW verze: mujweb.cz/kolafa/planet.html
[uvodsim/verlet.sh]
32/35 pch01
33/35 pch01
Teplota
Ekvipartièní princip: h m2 v2ixi = 12 kBT Poèet stupòù volnosti: f = 3N − fzachování ≈ 3N Ve standardním (mikrokanonickém) MD platí Epot + Ekin = const mìøíme: * Tkin =
Ekin
1k f 2 B
+ = hTkini
⇒
teplotu
34/35 pch01
Konstantní teplota v MD *
Pøibli¾né metody:
Tkin =
Ekin
1k f 2 B
+ = hTkini
pøe¹kálování rychlostí: ~r_ i,new = ~r_ i(T/Tkin)1/2 Berendsen (friction): ~r_ i,new = ~r_ i(T/Tkin)q, q < 1/2 Pøesný kanonický soubor = þsystém ve styku s termostatemÿ Maxwell-Boltzmann/Andersen: systematicky/náhodnì vybereme novou rychlost podle Maxwellova-Boltzmannova rozlo¾ení Langevin (Brownovská dynamika): aplikujeme malé náhodné síly a zároveò chladíme (þtøeníÿ) Nosé-Hoover pøidáme dal¹í stupeò volnosti (þpíst entropieÿ) Canonical Sampling through the Velocity Space (pøe¹kálování se stochastickým q) Konstantní tlak
. . . podobnì, ale trochu víc matematiky . . .
The End
[showvid /home/jiri/macsimus/ray/dogrun/dogrun.vid]
35/35 pch01
