Ugrás a tartalomhoz

A kvantumkémia alapjai és alkalmazása

Veszprémi Tamás, Fehér Miklós

Educatio Társadalmi Szolgáltató Nonprofit Kft.

16. fejezet - Molekuladinamikai szimulációk

16. fejezet - Molekuladinamikai szimulációk

A 7. fejezetben láttuk, hogy a számítástechnika rohamos fejlődésének, és a kvantumkémia új algoritmusainak köszönhetően ma már többszáz – maholnap többezer – atomos rendszereken is rutin számítások végezhetők. Ezzel kapcsolatban mindeddig gondosan elkerültünk egy alapvető kérdést. Vajon van-e egyáltalán értelme ilyen nagy rendszerek számításának? Egy nagymolekulában számos flexibilis kötés és forgási centrum lehetséges, ennek megfelelően a molekulától függően sokszáz, vagy akár milliónyi konformer is elképzelhető. Vajon melyik az „igazi”? Melyik tartozik a legstabilisabb minimumhoz, a globális minimumhoz? A sok közül melyek valósulnak meg a molekula folytonos mozgása során? Nyilvánvalóan egy konformer felbukkanása, illetve statisztikus súlya az anyagban az összes többi konformer relatív energiájától, a hőmérséklettől és egyéb termodinamikai állapotjelzőktől függ, ám az anyag minden egyes mérhető sajátsága az összes konformer tulajdonságából átlagolódik. Ennek alapján egyetlen konformer számítása, legyen az bármilyen pontos is, értéktelen, és semmitmondó. Komoly kihívás, melyre a molekuladinamika próbál válaszolni.

A molekuladinamika atomok és molekulák statisztikus mozgásának szimulálására alkalmazható módszer. Használata kiterjed kondenzált fázisú rendszerek, oldatok, makromolekulák, kolloidok és szervetlen molekulák modellezésére.

A molekuladinamikai szimulációkba a molekulák mozgását explicit módon veszik figyelembe. A szimuláció során rövid időlépésenként „mintát” veszünk és a részecskék helyzetét, valamint elmozdulását az előző lépésben kapott eredményből Newton klasszikus mozgásegyenlete alapján számítjuk. Először meghatározzuk az egyes atomokra ható erőket (a potenciális energia első deriváltjának –1-szerese), a gyorsulást, a végsebességeket, az egyes atomok helyzetét, és kiszámítjuk a rendszer potenciális és kinetikus energiáját. Ennek alapján meg tudjuk jósolni, hogy a rendszerrel mi fog történni egy kis időtartammal később. Végül is ezt a procedúrát ismételgetjük minden egyes időlépésben. A szimulációk alapvető paraméterei a kiindulási szerkezet, a rendszer energiája illetve hőmérséklete, az időlépések hossza, valamint a tanulmányozott időtartam.

A kiindulási szerkezet a végeredmény szempontjából alapvető fontosságú. Ehhez, fel is használnak minden lehetséges kísérleti és elméleti információt. A számítások általában egy optimált molekulaszerkezetből indulnak ki. Természetesen nem szükséges, hogy ez a rendszer globális energiaminimuma legyen, de mindenesetre reális, az anyagban gyakran előforduló, azaz nagy statisztikus súlyú szerkezetnek kell lennie.

Egy N atomból álló rendszer klasszikus mechanikai leírásához 6N változóra van szükség (3N helykoordináta és 3N sebesség vektor). A rendszer energiája a rendszer hőmérsékletétől függ, a hőmérséklet pedig a kinetikus energiával arányos. A szimuláció során a hőmérsékletet a rendszer összes atomjának kinetikus energiájából határozzák meg:

Ha a rendszer hőmérsékletét 0 K-nek tekintjük, ez egyenértékű azzal, hogy az atomok kiindulási sebessége zérus. Vagy ebből indulunk ki, vagy minden atomnak többé-kevésbé azonos energiát „osztunk ki”. Első lépésben az atomoknak kezdeti sebességeket adunk, és körülbelül száz lépésen keresztül hagyjuk a részecskéket bolyongani, azaz a rendszert relaxálódni, hogy elérje a Boltzmann-eloszlást. A rendszer „melegítése” a sebességvektorok folyamatos növelését jelenti. A leggyakrabban alkalmazott skálafaktor a következő alakú:

ahol T0 a környezet, vagy ha jobban tetszik, egy hipotetikus termosztát hőmérséklete, T a rendszer hőmérséklete, Δt az időlépés, és τ a csatolási paraméter (idő dimenziójú). τ határozza meg a rendszer és a környezete közti kapcsolat jellegét. Kis τ értékek esetén a rendszer gyakorlatilag felveszi a környezete hőmérsékletét. Ez hasznos lehet, ha egy adott hőmérsékletet szeretnénk gyorsan elérni, vagy ha minimális energiájú szerkezetet szimulálunk 0 K-en. Nagyobb értékeket (általában 0,1 ps felettieket) kell használni, ha termodinamikai tulajdonságokat határozunk meg.

Az egyes lépések hosszát úgy kell megválasztani, hogy ezek sokkal rövidebbek legyenek, mint a legrövidebb idő, mely alatt bármi érdekes, vagy lényeges lejátszódhatna a rendszerben. Ez általában a femtoszekundum nagyságrendjébe esik, ami azt jelenti, hogy mondjuk egy C–H vegyértékrezgés rezgési ideje alatt kb. 10 szimulációs lépést teszünk. Mindazonáltal az adott problémához illesztve nagyon gondosan kell megválasztani az időlépések hosszúságát. A nagyfrekvenciás mozgások (pl. egy kötés megnyúlása) leírásához rövid időintervallumokra van szükség, de ha a lépés túl rövid, csak erőforrást pocsékolunk. A rendszert célszerű kis lépésekben melegíteni, úgy, hogy a rendszer minden lépésben érje el az egyensúlyt. Az egyensúlyt gyakran alkalmazzák annak megállapítására, hogy a rendszer időbeli változása, fejlődése valóban egyensúlyi tulajdonságokat szimulál-e. Az egyensúly korrigálja az atomok sebességeit – a hibák a kiindulási rendszerben uralkodó feszültségekből és a kezdeti sebességek véletlen eloszlásából fakadnak. Az egyensúlyt egy sorozat dinamikai lépés végrehajtása során érjük el, rendszerint 5-100 ps szimulációs idővel. Természetesen minél hosszabb a szimulációs idő, annál jobban megközelíthetjük az egyensúlyt. Egyensúlynak azt az állapotot tekintjük, amikor a kinetikus energia (avagy a hőmérséklet), a teljes energia és néhány további tulajdonság ingadozása az egymás utáni lépések során elegendően kicsivé válik. A molekuladinamikai számítások során általános probléma, hogy a legtöbb fizikai és kémiai folyamat teljes modellezéséhez nanoszekundumokra vagy még hosszabb intervallumokra lenne szükség, de ehhez reménytelenül sok lépés kellene. Ez az, ami korlátozza a teljes szimulációs időt, vagy a tanulmányozható rendszerek méretét – és ez az, ami a módszer alkalmazhatóságának legnagyobb gátja.

A molekuladinamikai számítások előfeltétele, hogy ki tudjuk fejezni a molekulák között ható erőket az intermolekuláris (magok közötti) távolságok függvényeként. A klasszikus molekuladinamikai szimulációban empirikus (pl. az ún. Lennard–Jones potenciált,1 a molekulamechanikából átvett potenciálokat2 használnak számítási célra, a molekuladinamika kvantummechanikai leírása során viszont explicit kvantummechanikai potenciált3 alkalmaznak. Az utóbbi módszer használatához szükség van arra, hogy a mozgásegyenletek tartalmazzák rendszer hullámfüggvényét. Különösen előnyös kvantum-molekuladinamikát alkalmazni olyan esetekben, amikor kötések alakulhatnak ki és bomolhatnak fel, hiszen ezeket a folyamatokat csak kvantumkémiai számítások segítségével lehet egyértelműen vizsgálni.

A szimuláció során a rendszer teljes és kinetikus energiája is változhat. Amennyiben a rendszer egy valószínűtlen, nagyenergiájú állapotban van, kinetikus energiája nagy, ha viszont stabilis minimumban van, akkor kicsi, ekkor viszont potenciális energiája nő meg. A molekuladinamika alkalmas a rendszer termodinamikai tulajdonságainak meghatározására. Mivel a termodinamika egymással kölcsönhatásban lévő részecskéket tartalmazó rendszerekkel foglalkozik, ilyen esetekben a számítások során alapvető a sokaság (ensemble) fogalma. Formálisan úgy képzeljük, hogy egy sokaságban a rendszert sok azonos és egymástól független részecske építi fel, amelyek egy nagy térben helyezkednek el. A molekuladinamikában a sokaság két legfontosabb típusa a kanonikus és a mikrokanonikus sokaság. Mivel a rendszer kinetikus energiája megadja a hőmérsékletet, ezen alapul a molekuladinamika két alternatív formája: az egyikben a hőmérséklet, a másikban az energiát választjuk állandónak. A statisztikus mechanika nyelvére lefordítva a T = konstans feltétele az állapotok kanonikus sokaságát hozza létre, míg a konstans energiájú molekuladinamika mikrokanonikus sokaságot produkál. A kanonikus sokaságban minden egység azonos összetételű, térfogatú és hőmérsékletű (NVT). A mikrokanonikus sokaságban az állandó hőmérséklet helyett, minden egység azonos energiával rendelkezik (NVE). Ez utóbbi követelményt könnyebb kielégíteni. Ha viszont kanonikus sokaságot használunk, a hőmérséklet állandó értéken tartásához minden egyes időlépésben újra kell skálázni a hőmérsékletet a (16.2) egyenletnek megfelelően. Ugyanis ahogy a rendszer az egyensúly felé halad, a potenciális energia egyre nő, a kinetikus energia pedig csökken: a rendszer tehát hűl. Ha T = konstans, ilyenkor energiát kell a rendszerbe pumpálni. Ha viszont a rendszer valamely kisenergiájú, instabilis állapotba kerül, hűteni kell, azaz energiát elvonni.

A molekuladinamika fontos eszköz a folyadékok és oldatok leírásakor. Ilyenkor a rendszert rögzített méretű cellán belüli rögzített számú részecskével szimuláljuk. A kisszámú molekula miatt nagy nehézséget okozhatnak a felületi hatások a cella széleinél. Ezt megelőzhetjük periodikus határfeltételek alkalmazásával (l. 15. fejezet).

A kvantum-molekuladinamika az előző fejezetekben tárgyalt módszerek alternatívája lehet egy rendszer geometriájának meghatározásakor, vagy reakcióutak felderítésekor. Ugyanakkor jóval hatásosabb azoknál. Lehetővé teszi molekulák sokaságának vizsgálatát úgy, hogy a rendszernek nemcsak egy, hanem több kisenergiájú állapotát kapjuk meg, pl. feltérképezhetjük a rendszer összes konformerjét.Minthogy a rendszer fejlődése, változása során a fontosabb konformereket biztosan eléri, amennyiben az egyes konformerek statisztikus súlyát – az ott töltött időt – vizsgáljuk, könnyen és egyértelműen megkaphatjuk a rendszer legmélyebb, globális minimumát. A molekula viselkedését nemcsak elszigetelt körülmények között (pl. vákuumban) modellezhetjük, hanem különböző oldatokban is. Molekulákat ütköztetve képesek lehetünk valós kémiai reakciók lépéseinek vizsgálatára.

Persze a többlet információnak megvan az ára. Mivel lényegében minden időlépésben egy kvantumkémiai számítást kell elvégeznünk, a teljes időintervallumra ez akár 104105 számítást jelenthet. Ennek nyomán a kvantum-molekuladinamikai számítások, különösen azok, amelyek ab initio számításokat használnak, meglehetősen drágák, és csak a legegyszerűbb rendszerek esetében használhatók. A gyakorlatban legtöbbször a szemiempirikus módszerek alkalmazása a lehetséges alternatíva. E terület fejlődésének újabb állomása a sűrűségfunkcionál elmélet alkalmazása, amely jó minőségű szimulációt tesz elérhetővé.

16.1. A sósav disszociációja vizes oldatokban

R. Buesnel, I. H. Hillier és A. J. Masters munkája alapján (Chem. Phys. Letters, 247 , 391–394 (1995).

Minden kémikus számára elemi tudnivaló, hogy a sósav vizes közegben ionosan disszociál:

Az egyenlet azonban a fenti formájában nem tükrözi a valóságos folyamatot: egyértelmű bizonyítékaink vannak arra ugyanis, hogy egyetlen HCl-molekula egyetlen vízmolekulával való kölcsönhatása során nem disszociál. Azonnal felvetődik a kérdés: vajon hány vízmolekulának kell jelen lennie ahhoz, hogy a HCl disszociáljon? A kérdés nem csak teoretikus: számos, gyakorlati szempontból is fontos természeti jelenség megértésének alapkérdése ugyanez. Például hogyan viselkedik a sósav a poláros sztratoszférikus felhők felszínén úgy 180 K hőmérséklet környékén? Ha a klórgáz képződésének folyamatát – és ennek alapján a légköri ózonpaizs lebomlását – szeretnénk tisztázni a felső atmoszférában, ez alapvető kérdés. Más jellegű, ám lényegében ugyanazon kérdés: várható-e, hogy a sósav disszociál a jég felszínén? A vízmolekulák nem túl mozgékonyak a jégben, így ha nagyszámú molekulára van szükség, akkor valószínűtlen, hogy a vízmolekulák a jég belsejében szerepet játszhatnának a HCl disszociációjában.

Elvileg ilyen és hasonló problémák a szokásos kvantumkémiai módszerek segítségével is tanulmányozhatók, csupán ki kell számítanunk a rendszer teljes energiáját különböző disszociált és disszociálatlan elrendeződésekben. Ámde mely elrendeződéseket válasszuk ki, és honnan tudhatjuk, hogy kiválasztottuk-e az összes fontosat? A probléma könnyen megoldható egy-két vízmolekula esetén. Kitartó munkával még 3-4 vízmolekulából álló klaszter összes fontosabb elrendeződése is megtalálható, de nagyobb rendszerekben a dolog egyre bonyolultabbá válik. A molekuladinamika megelőzi a nehézséget azáltal, hogy statisztikusan vizsgál nagyszámú elrendeződést (melyek a hosszú szimulációs időben előfordulnak), és így alkalmasabb a lehetséges konformációk szisztematikus vizsgálatára.

A szimuláció minden egyes időlépésében kvantumkémiai számításokat kell végezni. Az itt ismertetendő munka szimulációs időtartama 100 ps volt, ami femtoszekundumos lépéseket figyelembe véve 100 000 számítást jelent! Logikus tehát, hogy a szerzők a viszonylag kis gépidő-igényű félempirikus AM1- és MNDO-szinteket választották. A félempirikus módszerek közelítéseiből adódóan persze nem várhatók precíz eredmények a HCl disszociációjára és a víz protonálódására. Például a víz AM1-szinten számított protonaffinitása 469 kJ/mol, szemben a kísérleti 706 kJ/mol értékekkel.4 MP2/6-311++G** szinten ez az adat jóval pontosabb, 722 kJ/mol. Ezért ahhoz, hogy megbízható adatokhoz jussunk, az AM1-eredményeket mindenképpen korrigálni kell.

Először egyedi szimulációkat hajtottak végre a protonált és nem-protonált H3O+Cl és H2OHCl rendszereken. Ezekben a számításokban az optimált ab initio egyensúlyi geometriákat használták és a molekulák szerkezetét a szimuláció során változatlannak vették. Ahhoz, hogy a folyamat energetikáját vizsgálhassuk, számítani kell a két rendszer (protonált és nem-protonált) energiakülönbségét. A H3O+Cl és H2OHCl formák relatív energiái 295 kJ/mol az AM1 és 194 kJ/mol az MP2 szinteken. A két adat különbségét tekinthetjük az AM1 szint hibájának.

A molekuladinamikai számítások stratégiája a következő volt: először megadott számú vízmolekulával szolvatált H2OHCl és H3O+Cl rendszereket vizsgáltak AM1-szinten. Ezután meghatározták az egyes klaszterek relatív energiáit, melyeket korrigálták az AM1-számítás hibájával. Ezt aztán megismételték különböző méretű víz klaszterekre a vízmolekulák számát egész 15-ig növelve.

Általában a molekulák kezdeti sebessége megadható pl. valamilyen sebességeloszlás függvény segítségével. Ebben a munkában a kezdeti sebességeket egyszerűen zérusnak vették, ami formálisan azt jelenti, hogy a rendszer 0 K-en van a szimuláció indításakor. Az első feladat a szimuláció kezdeti időszakában a molekulák kinetikus energiájának a növelése. A rendszert íly módon 1 ps alatt 0 K-ről 180K-re „melegítették” – ez éppen a sarki sztratoszférikus felhők hőmérsékletével egyezik. (Az 1 ps megfelel 1000 db 1 fs hosszúságú időlépésnek, azaz 1000 AM1-számításnak.) A „valódi” szimuláció csak ezután következik, és ennek során a rendszer termikus egyensúlyba jut. Ez 99 ps-os időtartamot, vagyis az összes atom helyzetének 99 000-szer történő kiszámítását jelentette. Minden szimuláció utolsó 90 000 lépésének eredményeit átlagolták.

 

16.1. táblázat. Átlagos kötési energia (kJ/mol) a HClH2O és H+ClH2O klaszterekben.

 

   

A molekuladinamikai számításokból nyert konformációk érdekes jellegzetességeket mutatnak. A H2OHCl rendszerekben a HCl-molekula mindig a klaszter szélén jelenik meg, a klóratom kifelé irányul a vízmolekulák közül. A H3O+Cl rendszerekben a környező vízmolekulák egy-egy hidrogénatomja a Cl, illetve a H3O+ oxigénje irányába mutat. A várakozásnak megfelelően a vízmolekulák mozgása sokkal koordináltabb az ionos klaszterekben, mivel az ion erősebb kölcsönhatásba lép az oldószermolekulákkal, mint a semleges molekula.

A klaszterek átlagos kötési energiája úgy definiálható, hogy az egyes klaszterek energiáját viszonyítjuk az izolált dimer és a vízmolekulák teljes energiájának az összegéhez. Az AM1-szimulációkból kapott értékek láthatók a 16.1. táblázatban a HClnH2O és H+ClnH2O rendszerek esetében (n=2-től 15-ig). Amint várható, a növekvő számú vízmolekula extra kölcsönhatásainak eredményeként a kötési energia a rendszer méretével nő. A kölcsönhatások nagyobbak az ionos klaszterekben annak megfelelően, hogy az ion-dipólus kölcsönhatások általában erősebbek, mint a dipólus-dipólus kölcsönhatások.

Talán még érdekesebb a kétfajta dimer relatív szolvatációs energiájának a változása. Ezt is bemutatjuk a 16.1. táblázatban, valamint grafikusan a 16.1. ábrán. Amint látható, a szolvatációs energia minden esetben nő a klaszter méretével, ám az ionos klaszteré a hosszútávú elektroszatikus kölcsönhatásoknak megfelelően erőteljesebben nő. A relatív szolvatációs energia a mérettel növekszik, de a növekedés egy határérték felé tart. A 16.1. ábrán az izolált dimerek MP2-számításokból kapott energiakülönbségét vízszintes vonallal ábrázoltuk. Látható, hogy a relatív szolvatációs energia ezt az energiakülönbséget 6 vízmolekulánál haladja meg. Ezek szerint a HCl ionizációja akkor válik fontossá, ha 6 vagy annál több vízmolekula veszi körül a sósavat. Az eredmény valószínűsíti, hogy a HCl a jégkristályok felszínén, valamint a felső atmoszférában is képes disszociációra.

 

16.1. ábra. Ionizált és nemionizált sósav AM1 molekuladinamikai szimulációkból kapott relatív szolvatációs energiái a H2O-molekula számának függvényében.