Úlohy do cvičení z výpočetní chemie Tomáš Zelený
Petr Sklenovský 2008
Michal Otyepka
2
Poděkování Tato sbírka úloh vznikla za podpory grantu Fondu rozvoje vysokých škol FRVŠ 2261.
Tento text byl vysázen typografickým systémem LATEX 2ε .
i
ii
Seznam zkratek AO ATP BO BSSE CC CI CPU DFT DNA FCI (FullCI) GVB HDD HF IR MCSCF MD MM MO MP MP2 NMR OPT PDB PES QMC RAM RCSB RNA RTG SP STO UV VIS VMD vdW
atomový orbital adenosintrifosfát Born-Oppenheimerova aproximace bázová superpoziční chyba (Basis Set Superposition Error) Coupled Cluster theory Configuration Interaction procesor (Central Processor Unit) Density Functional Theory deoxyribonukleová kyselina metoda úplné konfigurační interakce Generalized Valence Bond pevný disk (Hard Disc Device) Hartree-Fockova metoda infračervené záření Multi-Configurations Self Consistent Field molekulová dynamika molekulová mechanika molekulový orbital Møller-Plessetova poruchová teorie Møller-Plessetova poruchová teorie druhého řádu nukleární magentická rezonance optimalizace Protein Data Bank hyperplocha potenciální energie (Potential Energy Surface) Quantum Monte carlo operační paměť (Random Access Memory) Research Collaboratory for Structural Bioinformatics ribonukleová kyselina rentgen single point Slater Type Orbital ultrafialové záření viditelné záření Visual Molecular Dynamics van der Waals iii
Obsah I
Teorie
2
1 Co je výpočetní chemie?
3
2 Metody výpočetní chemie 2.1 Ab initio výpočty . . . . . . . . . . . . . . . . . . . . . 2.2 Semiempirické metody . . . . . . . . . . . . . . . . . . 2.3 Molekulová mechanika . . . . . . . . . . . . . . . . . . 2.4 Molekulové simulace . . . . . . . . . . . . . . . . . . . 2.4.1 Molekulárně-mechanická molekulová dynamika 2.4.2 Protokol molekulárně-mechanické dynamiky . . 2.5 Statistická termodynamika . . . . . . . . . . . . . . .
II
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
Praktické úlohy
4 4 5 6 8 9 10 11
12
3 Úlohy z kvantové chemie 3.1 Stavba molekul . . . . . . . . . . . . . . . . . . 3.1.1 GaussView . . . . . . . . . . . . . . . 3.2 Optimalizace vs. Single Point . . . . . . . . . . 3.3 Sken povrchu potenciální energie . . . . . . . . 3.3.1 Nerelaxovaný . . . . . . . . . . . . . . . 3.3.2 Relaxovaný . . . . . . . . . . . . . . . . 3.4 Mezimolekulové komplexy . . . . . . . . . . . . 3.5 Výpočet frekvencí a termodynamických veličin 3.5.1 Výpočet frekvencí . . . . . . . . . . . . 3.5.2 Termodynamika . . . . . . . . . . . . . 3.5.3 Kinetika chemických rekcí . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
13 13 13 16 17 18 18 19 20 21 23 24
4 Úlohy z molekulové mechaniky 25 4.1 Vliv solventu na strukturu peptidu . . . . . . . . . . . . . . . 25 4.2 Sken povrchu potenciální energie podle jednoho parametru . 27 iv
5 Úlohy z molekulové simulace 29 5.1 Simulované žíhání . . . . . . . . . . . . . . . . . . . . . . . . 29 5.1.1 Nalezení globálního minima acetonu . . . . . . . . . . 29 5.1.2 Nalezení globálního minima peroxidu vodíku . . . . . 31 6 Úlohy z bioinformatiky 6.1 Databáze proteinů a nukleových kyselin RCSB 6.1.1 Stručná historie . . . . . . . . . . . . . . 6.1.2 Práce s databází . . . . . . . . . . . . . 6.2 Visual Molecular Dynamics . . . . . . . . . . .
III
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
Dodatky
7 Reprezentace molekul 7.1 Formáty souborů . . . . 7.1.1 Formát Z-matice 7.1.2 Formát XYZ . . 7.1.3 Formát PDB . .
32 32 32 32 35
40 . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
41 41 41 42 43
8 Povrchy potenciální energie
45
IV
47
Přílohy
9 Úlohy z kvantové chemie 9.1 Stavba molekul . . . . . . . . . . 9.1.1 GaussView . . . . . . . 9.2 Sken povrchu potenciální energie 9.2.1 GaussView . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
48 48 48 49 49
10 Úlohy z molekulové mechaniky 51 10.1 Sken povrchu potenciální energie . . . . . . . . . . . . . . . . 51 10.1.1 HyperChem . . . . . . . . . . . . . . . . . . . . . . . 51
v
vi
Úvod P
ředmět Molekulární modelování (KFC/MOMO) vyučovaný na Přírodovědecké fakultě Univerzity Palackého v Olomouci je zaměřen na získání základních znalostí z oblasti výpočetní chemie. Vzhledem ke skutečnosti, že je molekulární modelování dynamicky se rozvíjející vědní disciplínou, je nezbytné neustálé sledování aktuálních trendů a jejich včasné začlenění do výuky. Prioritou tohoto učebního textu je zařazení nových úloh do kurzu molekulárního modelování, které reflektují aktuální stav ve výpočetní chemii se zaměřením na výpočet termodynamických veličin aparátem kvantové mechaniky a statistické termodynamiky, strukturní databáze biomolekul (http://www.pdb.org), vizualizační software, základní principy molekulové mechaniky a dynamiky. Úlohy byly voleny tak, aby studentům byly nápomocny pro hlubší pochopení probírané problematiky a jejího významu v aplikacích. Studenti se samostatně zapojili do vyhledání potřebných informací nezbytných pro úspěšné řešení práce, vybudovali si strategii pro řešení problému za odborné pomoci lektora. Tento text je veřejně a zdarma zpřístupněn na univerzitní síti, konkrétně na webových stránkách Katedry fyzikální chemie (http://fch.upol.cz), kde je průběžně aktualizován.
J sme si vědomi, že mnohá zde uvedená tvrzení nejsou přesná a zda-
leka už ne úplná a vyčerpávající, dopustili jsme se mnoha zjednodušení, která mají za cíl spíše nastínit podstatu problému pro jednoduché pochopení než vyčerpávajícím způsobem dokázat všechna fakta. Berte proto tuto příručku jako úvod do problematiky teoretické chemie pro nezasvěcené a ne jako komplexní učebnici. Autoři
1
Část I
Teorie
2
Kapitola 1
Co je výpočetní chemie? Výpočetní chemie je rychle se rozvíjející vědní disciplína, která poskytuje fyzikálně-chemická data vypočtená pomocí teoretických metod jako je molekulární modelování, počítačové simulace systémů jako jsou biomolekuly (proteiny, DNA, RNA), polymery, drogy, anorganické a organické molekuly. Co nám poskytuje výpočetní chemie? • Energie a struktura molekul • Optimalizace geometrie vycházající z empirické představy • Energie a struktury tranzitních stavů • Vazebné energie • Reakční energie a termodynamické vlastnosti • Molekulové orbitaly • Atomové náboje a elektrostatický potenciál • Vibrační frekvence • Infračervená a Ramanova spektra molekul • NMR spektra • Magnetické vlastnosti • Polarizability a hyperpolarizability • Reakční cesty • Ionizační potenciály, elektronové a protonové afinity • Modelování excitovaných stavů 3
Kapitola 2
Metody výpočetní chemie Výpočetní chemie se skláda z teoretické (strukturální) části nazývané molekulární modelování a modelování procesů, známé pod pojmem molekulární simulace.
2.1
Ab initio výpočty
Latinský termín ab initio znamená „od počátkuÿ. Výpočty tohoto typu vycházejí přímo z teoretických principů (Schrödingerova rovnice) bez použití empirických dat. Někdy se jim také říká metody kvantově mechanické. Aproximace používané v ab initio metodách jsou často matematické, jako zjednodušení funkce nebo nalezení aproximativního řešení diferenciálních rovnic. Nejběžnějším typem ab initio výpočtu je Hartreeho Fockova metoda (HF). Nejvýznamnější aproximací HF metody je zanedbání (nezahrnutí) Coulombické elektron—elektronové repulze ve výpočtu, rovněž se lze setkat s formulací, že HF metoda nezahrnuje elekronovou korelaci. Metoda HF je variační, to znamená, že získané apoximativní energie jsou vyšší nebo rovny energii exaktní. Druhou aproximací v HF metodě je zavedení tzv. jednoelekronových funkcí, protože nelze pomocí HF metody exaktně řešit systémy s více než jedním elektronem. Často používanými funkcemi jsou lineární kombinace Slaterových orbitalů (e−ax ) nebo orbitalů, které jsou tvořeny 2 Gaussovou funkcí (e(−ax ) ), zkráceně STO a GTO. Vlnová funkce je tvořená lineární kombinací atomových orbitalů, resp. lineární kombinací bázových funkcí. Díky této aproximaci je Hartree Fockova energie často vyšší než energie exaktní. Popis orbitalů lze teoreticky vylepšovat až do nekonečna, použitím lineární kombinace „nekonečného množství bázových funkcíÿ, výpočtem energie s takovým množstvím bázových funkcí získáme energii, která odpovídá tzv. Hartee Fockově limitě. Pro označení bázových funkcí se často používá zkratek jako STO-3G, 6-31G*, 6-311++G**.
4
Většina ab initio výpočtů začíná výpočtem HF metodou, kterou následují tzv. post HF metody, které v různé míře zavádějí korekci pro explicitní elektron—elektronovou repulzi (zahrnutí elektronové korelace). K post HF metodám patří například Møller-Plessetova poruchová teorie (MPX, kde X označuje řád korekce), Generalized Valence Bond (GVB), multikonfigurační metoda self konzistentního pole (MCSCF – Multi-Cofigurations Self Consistent Field), metoda konfigurační interakce (CI – Configuration Interaction) a metoda vázaných klastrů (CC – Coupled Cluster theory). Metoda Quantum Monte Carlo (QMC) řeší mnoho nedostatků HF metody. Tato metoda bývá řešena mnoha různými způsoby, jmenovitě variační, difůzní a pomocí Greenových funkcí. Tyto metody pracují s explicitní elekronovou korelací přímo zahrnutou ve vlnové funkci použitím Monte Carlo integrací. Tyto výpočty jsou bohužel velice časově náročné. Alternativou ke klasickým ab initio metodám jsou metody nazývané Density Functional Theory (DFT) používající namísto vlnové funkce tzv. funkcionál hustoty, tedy je zde brán do úvahy vztah energie a elektronové hustoty. Závěrem je potřeba dodat, že ab initio výpočty nám poskytují velice kvalitní výsledky (především kvalitativně), které jsou bohužel vykoupeny obrovskou náročností ab initio metod na výpočetní techniku (počet a výkon CPU, RAM, HDD), to jsou limitující faktory. Těmi nejlepšími ab initio metodami lze v současné době počítat sytémy o velikosti řádově v desítkách atomů. Metoda DFT umožňuje pracovat se systémy, které jsou o něco větší.
2.2
Semiempirické metody
Mechanismus semiempirických výpočtů je v základu stejný jako u Hartreeho Fockovy (HF) metody s tím rozdílem, že některé části jsou aproximovány třeba empirickými parametry nebo úplně zanedbány jako například dvouelektronové integrály. Korekce těchto hrubých aproximací ev. zanedbání jsou prováděny parametrizací, fitováním tak, aby se výsledky co nejlépe shodovaly s experimentálními daty. Výhodou těchto metod je jejich výpočetní jednoduchost ve srovnání s ab initio metodami. Nevýhodou mohou být v některých případech zavádějící data při výpočtech molekul hodně se lišících od těch databázových, na kterých byla parametrizace dané semiempirické metody prováděna. Semiempirické metody lze s výhodou používat u organických molekul, kde je jen málo odlišných prvků a paramatrizace je tak snadná. Speciální případy, jako anorganické komplexy a jiné nestandardní případy je potřeba počítat přesnějšími metodami. 5
2.3
Molekulová mechanika
Chceme-li ab initio metody aplikovat na systémech obsahujících velké množství atomů (proteiny, nukleové kyseliny, polysacharidy), narazíme na překážku extrémních systémových nároků (potřeba astronomické velikosti RAM a rychlosti CPU). Z toho důvodu je třeba použít metod molekulové mechaniky (často označované jako metody silového pole nebo empirického potenciálu), které mají narozdíl od ab initio metod signifikantně menší systémové nároky. V molekulové mechanice (MM) popisujeme molekulu jako soubor bodových nábojů a vazeb mezi nimi, a proto můžeme využít k popisu klasickou newtonovskou mechaniku. Díky přístupu, který je podrobně vysvětlen v textu níže, je energie molekuly jednoduše dekomponovatelná na jednotlivé složky. Limitou v MM metodách je ztráta informace o distribuci elektronů, protože energie molekuly je vyjádřená pouze jako funkce vzájemné polohy jader (viz BO aproximace). Hlavním nedostatkem těchto metod je tedy znemožnění studia procesů, kde dochází ke změně elektronové distribuce, jako jsou například chemické reakce. Atomy v MM reprezentujeme jako sférické částice s určitým nábojem (lokalizovaným v jejich středu) a poloměrem (vdW poloměr). V místech, kde mezi sebou atomy tvoří kovalentní vazbu aplikujeme harmonickou funkci (potenciál), jejíž analytický tvar se volí podle pevnosti vazby. Na chemické vazby se tedy v molekulové mechanice nahlíží jako na pružiny s charakteristickou délkou a tuhostí. Model biatomické molekuly Biatomická molekula je dvojice atomů spojená kovalentní vazbou (např. molekula vodíku, oxid uhelnatý, molekula dusíku). V molekulové mechanice je modelem biatomické molekuly systém dvou koulí, které jsou spojené pružinou. Potenciálová funkce biatomické molekuly vypadá následovně: kb (r − r0 )2 , (2.1) 2 kde r0 je rovnovážná vzdálenost mezi atomy, r je okamžitá vzdálenost a kb je silová konstanta (pevnost) vazby. Potenciálová funkce biatomické molekuly je kvadratická funkce. Popis disociace atomů biatomické molekuly je prostřednictvím kvadratické funkce nemožný (energie roste s rostoucí vzdáleností do nekonečna). Na druhou stranu funkce dobře popisuje oblast rovnovážné meziatomové vzdálenosti (obrázek 2.1). E(R) =
6
Obrázek 2.1: Morseho a harmonický potenciál. Pro biatomickou molekulu (např. molekulu vodíku) představuje Morseho, resp. harmonický potenciál celkový povrch potenciální energie této molekuly. Jak je patrné z průběhu kvadratické funkce, aplikací harmonického potenciálu nebereme v úvahu disociaci molekul. Model molekuly vody Molekula vody obsahuje celkem tři atomy (dva vodíky a jeden kyslík) vzájemně spojené dvěma kovalentními vazbami (O—H). Obě vazby O—H stejně tak jako úhel H—O—H jsou popsány kvadratickým potenciálem. Silová konstanta pnutí vazby O—H a úhlu H—O—H představuje spolu s rovnovážnou hodnotou vzdálenosti H—O a úhlu H—O—H veškeré vazebné parametry modelu vody. Potenciálovou funkci molekuly vody vyjadřuje následující vztah:
E(R) =
2 X kb i=1
2
(ri − r0 )2 +
ka (θ − θ0 )2 , 2
(2.2)
kde r0 je rovnovážná vzdálenost mezi atomy O—H, r je okamžitá vzdálenost mezi atomy O—H, θ0 je rovnovážná hodnota úhlu H—O—H, θ je okamžitá hodnota úhlu H—O—H, kb je silová konstanta vazeb O—H a ka je silová konstanta deformace úhlu H—O—H. 7
Model molekuly peroxidu vodíku Molekula peroxidu vodíku obsahuje celkem čtyři atomy (dva atomy kyslíku a dva atomy vodíku). V systému jsou celkem tři vazby (dvě vazby O—H a jedna vazba O—O), dva vazebné úhly (H—O—O a O—O—H) a jedna torze (H—O—O—H). Pro konstrukci potenciálové funkce je nezbytné zahrnout jak vazebné termy (vazby, úhly a torze), tak nevazebné složky (elektrostatická a vdW interakce mezi vodíky). Silové pole musí obsahovat tudíž následující soubor parametrů: náboje a velikosti všech atomů, rovnovážné délky vazeb H—O, O—O a O—H, rovnovážné hodnoty velikosti úhlů H—O—O a O—O—H, rovnovážnou hodnotu dihedrálního úhlu H—O—O—H a parametry potřebné pro popis nevazebné interakce vodíků. Musí být zohledněn fakt, že v plynné fázi zaujímá molekula peroxidu vodíku synklinální konformaci s dihedrálním úhlem 111,5◦ , jejíž přičinou je odpuzování mezi vazbami O—H a volnými elektronovými páry na každém atomu kyslíku. Peroxid vodíku je nejmenší molekula, o níž je známo, že má bráněnou rotaci kolem jednoduché vazby s rotačními bariérami 4,62 kJ · mol−1 pro izomer trans a 29,45 kJ · mol−1 pro izomer cis, synklinální forma je stálá v kapalné fázi, kde je stabilizovaná vodíkovými vazbami. Funkce popisující rotaci vazby H—O—O—H se proto volí taková, aby energetické rozdíly mezi jednotlivými konformery peroxidu vodíku odpovídaly experimentálně naměřeným hodnotám. Pro lepší představu je na obrázku 2.2 zobrazena molekula peroxidu vodíku ve svých jednotlivých konformacích.
Obrázek 2.2: (a) Molekula peroxidu vodíku v konformaci, kdy jí přísluší minimální energie. (b) Molekula peroxidu vodíku v nezákrytové konformaci. (c) Molekula peroxidu vodíku v zákrytové konformaci.
2.4
Molekulové simulace
Molekulová dynamika se snaží zachytit časový vývoj systému mnoha částic řešením jejich pohybových rovnic. Mapuje se tak konfigurační prostor, což přináší cenné informace o studovaném systému. Existuje několik způsobů 8
jak výpočetně popsat časovou evoluci systému.
2.4.1
Molekulárně-mechanická molekulová dynamika
Nejpoužívanější a zároveň nejjednodušší je molekulárně-mechanická molekulová dynamika (MM/MD). Jedná se o přístup nesmírně vhodný pro studium časové evoluce biomolekul (proteinů, nukleových kyselin a polysacharidů) v případech, kde není zapotřebí znalosti elektronových stavů. Systém je v MM/MD definován polohou a hybností jednotlivých částic. Částice se nacházejí v potenciálovém poli, které je tvořeno superpozicí potenciálu každé částice. Pro provedení samotné simulace je potřeba znát jejich souřadnice (polohy) a počáteční rychlosti (vektory rychlostí), které jsou přiřazovány náhodně podle Maxwell-Boltzmannova rozdělení pro danou teplotu. Známe-li potenciál U , známe i síly působící na jednotlivé částice F, protože platí F=−
∂U , ∂r
(2.3)
kde r je polohový vektor částice. Z druhého Newtonova pohybového zákona lze pak vypočítat i vektory zrychlení částic a. a=
F . m
(2.4)
Časový vývoj systému je řešen numericky s integračním krokem, který se volí v rozmezí od 1 do 2 fs. Částice má v n-tém kroku polohu rn a vektor rychlosti vn . Její časový vývoj je určen silovým polem F(rn ) a poloha částice v dalším kroku je dána rovnicí rn+1 = rn + vn · ∆t
(2.5)
vn+1 = vn + an · ∆t,
(2.6)
a její rychlost rovnicí
kde ∆t je integrační krok, rn resp. rn+1 jsou polohové vektory částice v n-tém resp. v n + 1-vém kroku, vn resp. vn+1 jsou vektory rychlosti částice v n-tém resp. n + 1-vém kroku a an je vektor zrychlení částice v n-tém kroku simulace definovaný rovnicí 2.7 an =
F(rn ) , m
přičemž m představuje hmotnost částice. 9
(2.7)
Druhý Newtonův zákon – Jestliže na těleso působí síla, pak se těleso pohybuje se zrychlením, které je přímo úměrné působící síle a nepřímo úměrné hmotnosti tělesa.
2.4.2
Protokol molekulárně-mechanické dynamiky
Každý uživatel software pro molekulovou simulaci si v prvé řadě musí uvědomit, že data k nímž se dopracuje, jsou závislá na dvou základních faktorech, které významně mohou ovlivnit kvalitu (průběh) celé simulace. V prvé řadě se jedná o volbu počáteční struktury (startovní struktury). Často se pracuje s RTG strukturami biomolekul s vysokým rozlišením (pod 1 ˚ A), přičemž se struktura dále rafinuje optimalizačními agoritmy a krátkými simulacemi. Druhým prvkem je volba správného silového pole. Takového, aby co nejlépe fitovalo na data získaná z experimentálních měření studovaného systému (stabilita, flexibilita apod.). Víme-li například, že studovaná molekula tvoří za in vitro podmínek jasně definovanou, termodynamicky stabilní strukturu, je nepřípustné, aby během simulace došlo, díky volbě nekorektního silového pole, k jejímu kolapsu. Abychom dosáhli věrohodných dat, sproštěných artefaktů pramenících z neznalosti či nesprávného používání algoritmů pro molekulovou simulaci, doporučuje se postupovat pode níže uvedených kroků. 1. Volba systému – počátečních souřadnic (modelování, databáze) 2. Volba relevantního silového pole. 3. Optimalizace geometrie vodíků solutu. 4. Optimalizace protiiontů a vod. Krystalové vody zůstávají fixovány. 5. Relaxace protiiontů a molekul vody v „boxuÿ kvůli zvýšení hustoty soustavy – [N pT ] podmínky. 6. Optimalizace postranních řetězců proteinu/bází nukleových kyselin s postupným uvolňováním silové konstanty na atomech páteře (aplikace „restrainuÿ. 7. Optimalizace geometrie všech atomů systému. 8. Krátká dynamika celého systému (50–100 ps) za účelem zvýšení teploty (termalizační fáze). Teplota se zvyšuje lineárně a obvykle z 10 na 300 K. Počáteční rychlosti (např. pro 10 K) jsou jednotlivým atomům přiřazeny dle Maxwell-Boltzmannovy distribuce pro danou teplotu. 9. Produkční část molekulové dynamiky. Délka se volí obvykle taková, aby se stačil systém ekvilibrovat (došlo ke konvergenci energie a strukturních parametrů) a fázový prostor byl dostatečně „samplovánÿ. 10
2.5
Statistická termodynamika
Aparát statistické termodynamiky nám poskytuje spojení mezi mikrosvětem, který lze výborně popsat kvantovou mechanikou a makrosvětem, který popisujeme „klasickouÿ termodynamikou. Díky přiblížení, které nám poskytuje statistická termodynamika můžeme získat z výsledků kvantové mechaniky termodynamické veličiny jako je entalpie (H), Gibbsova energie (G), entropie (S), tepelné kapacity (C) apod. Jak již název napovídá, snažíme se brát termodynamicky pozorovaný soubor ze statistického hlediska, tedy že stav systému obsahujícího miliony molekul, který popisujeme termodynamickými veličinami je dán mnoha mikrostavy. Každý mikrostav popisuje stav jedné molekuly. Pro dobrý popis celého systému je důležité zahrnout všechny možné stavy a určit jejich význam (váha stavu) a pravděpodobnost výskytu daného stavu. Když tyto předpoklady statisticky zhodnotíme, dostaneme makroskopické (termodynamické) chování celého systému vycházející z pozorování mikrosvěta.
11
Část II
Praktické úlohy
12
Kapitola 3
Úlohy z kvantové chemie 3.1
Stavba molekul
Existuje mnoho programů vhodných pro stavbu molekul od velikosti dvouatomové molekuly až po DNA nebo proteiny. My budeme používat program GaussView, který zároveň tvoří grafické rozhraní (frontend) pro výpočetní program Gaussian. Pro stavbu jednoduchých molekul lze rovněž použít www.gaussian.com program Molden a pro složitější struktury Accelrys DS Visualizer. www.molden.de Oba tyto programy jsou narozdíl od GaussView dostupné pro nekomerční www.accelrys.com účely zdarma pro všechny běžné operační systémy.
3.1.1
Prostředí aplikace GaussView
Na obrázku 3.1 je zobrazen program GaussView po spuštění. Vidíme, že program se spouští ve dvou oknech. Prázdné modré okno je okno editační, do něj se umísťují atomy, ev. fragmenty. Druhé okno obsahuje několik lišt nástrojů, které slouží jak ke stavbě molukuly, tak k editaci. Další nástroje se používají pro zobrazení výsledků výpočtů (orbitaly, spektra). Součástí okna s nástroji je i náhled, který zobrazuje vybraný fragment (v případě obrázku 3.1 je to methyl). V horní liště najdeme několik nástrojů, které nám stavbu molekuly nebo komplexu zjednoduší, jako je periodická tabulka, knihovna funkčních skupin pro organické molekuly, knihovna organických cyklických a aromatických sloučenin a nástroj pro stavbu biomolekul (proteinů a aminokyselin) na obrázku 3.2. Další nástroje umožňují snadnou editaci geometrických parametrů jako délka vazby, úhel a dihedrální úhel, viz obrázek 3.3.
Cvičení 1. Nakreslete molekulu vody v symetrii C2v s délkou O—H vazby 0, 8 ˚ A a úhlem H—O—H 102◦ . 13
Obrázek 3.1: Program GaussView po spuštění. 2. Nakreslete molekuly 1,2-dichlorethanu s různými dihedrálními úhly Cl—C—C—Cl (viz níže) a spočtěte vždy single point HF/3-21G, získané SCF energie zapište do tabulky a vyneste do grafu: (a) −180◦ (b) −135◦ (c) −90◦ (d) −45◦ (e) 0◦ (f) 45◦ (g) 90◦ (h) 135◦ (i) 180◦ Řešení Úloha 2: Nastavení a změna dihedrálního úhlu Nakreslíme 1,2-dichlorethan, na panelu nástrojů vybereme editaci dihedrálního úhlu (viz obrázek 3.3, třetí ikona zleva) a postupně levým tlačítkem myši klikneme na atomy v pořadí Cl—C—C—Cl, objeví se nám okno, jako na obrázku 3.4. V polích Atom 1 a Atom 4 musí být zvoleno Rotate groups, aby se spolu s chlóry okolo vazby mezi uhlíky otáčely i vodíky. Pak již lze měnit dihedrální úhel buď posuvníkem, nebo přímo zadat hodnotu žádaného úhlu do textového pole. 14
Obrázek 3.2: Přehled knihoven a nástrojů pro stavbu molekuly.
Obrázek 3.3: Nástroje pro úpravu geometrie molekuly (zleva: vazby, úhly, dihedrální úhly, dotaz na parametry).
Spuštění single point výpočtu z GaussView Pokud máme nastaven požadovaný dihedrální úhel, můžeme spustit single point (viz kapitola 3.2) výpočet. Z nabídky programu GaussView vybereme Calculate → Gaussian, otevře se okno pro zadání parametrů výpočtu (obrázek 3.5). Na záložce Job Type zvolíme Energy a záložku Method upravíme dle obrázku 3.5. Kliknutím na tlačítko Submit budeme požádáni o uložení „inputuÿ pro program Gaussian. Místo uložení vstupu si zapamatujeme (vstup i výstup lze prozkoumat pomocí obyčejného textového editoru). Nakonec budeme dotázáni, zda chceme, aby se uložený vstup rovnou spustil, zvolíme OK. Po dokončení výpočtu nám program GaussView nabídne otevřít výsledky výpočtu, zvolíme OK. Na panelu programu GaussView vybereme Results 15
Obrázek 3.4: Nástroj pro změnu dihedrálního úhlu.
Obrázek 3.5: Okno dialogu pro spuštění výpočtu.
SCF – Self Consistent Field, metoda používaná pro výpočet energie Hartreeho–Fockovou aproximací.
→ Summary, odtud opíšeme hodnotu SCF energie. Další možností pro nalezení SCF energie je prohledání výstupu programu Gaussian. To provedeme pro všechny hodnoty dihedrálního úhlu. Po vynesení hodnot do grafu by měl výsledek vypadat jako na obrázku 3.6.
3.2
Optimalizace vs. Single Point
Optimalizace (OPT) a Single Point (SP) jsou dvě základní procedury, které je potřeba odlišit. SP je nejjednodušší typ výpočtu, kdy na fixované pozici jader spočteme energii systému, tedy nedochází během výpočtu ke změně geometrie. Během optimalizace se narozdíl od SP geometrie molekuly mění ve směru nejbližšího minima na povrchu potenciální energie (PES, viz kapitola 8). Optimalizace tedy vyžaduje provést několik po sobě jdoucích SP výpočtů s různými pozicemi jader, které jsou odvozovány na základě výpočtu gradientu (první derivace energie podle souřadnic). Výsledkem je geometrie s co nejnižší energií. 16
Obrázek 3.6: Závislost potenciální energie na dihedrálním úhlu dichlorethanu. Cvičení 1. Nakreslete molekulu vody v symetrii C2v s délkou O—H vazby 0, 8 ˚ A a úhlem H—O—H 102◦ . (a) Spočtěte single point metodou MP2/6-31G(d), poznamenejte si energii. (b) Proveďte optimalizaci toutéž metodou a porovnejte energie a geometrické parametry (délky vazeb, úhly) mezi počáteční a optimalizovanou geometrií.
3.3
Sken povrchu potenciální energie
Pomocí programu Gaussian lze prováděti dva druhy skenu povrchu potenciální energie (relaxovaný a nerelaxovaný). V obou případech si zvolíme jednu proměnnou (vazebná délka, úhel, torze) a přitom sledujeme jak se mění potenciální energie sledovaného systému. Nerelaxovaný sken provádí pouhou 17
změnu hodnoty proměnné a v každém bodě spočte hodnotu „Single Pointÿ energie. Relaxovaný sken navíc zafixuje proměnnou na její požadované hodnotě a zbytek molekuly optimalizuje.
3.3.1
Nerelaxovaný
Příklad vstupu # hf/631g(d) scan H2O Scan 0 1 O1 H2 H3 B1 B2 A1
O1 O1
B1 B2
0.94741 0.94741 70. 11
H2
A1
10.
V route (příkazové) části vstupu je příkaz sken, geometrie molekuly je definována formou Z-matice a jedna nebo více proměnných Z-matice, v našem případě úhel A1 70. 11 10. je definován pro provedení skenu, kde 70. udává počáteční hodnotu úhlu, 11 počet kroků a 10. velikost kroku. Program tedy provede nerelaxovaný sken povrchem potenciální energie podél změny úhlu v molekule vody a to v krocích 70◦ , 80◦ , 90◦ , 100◦ , 110◦ , 120◦ , 130◦ , 140◦ , 150◦ a 160◦ .
3.3.2
Relaxovaný
Příklad vstupu #opt=modredundant hf/631g(d) H2O Scan with optimization at each angle 0 1 O1 H2 H3 B1 B2 A1
O1 O1
B1 B2
H2
A1
0.94741 0.94741 70.
A 2 1 3 S 11
10.
V příkazové části je potřeba k příkazu opt přidat =modredundand a na konec vstupního souboru přidat řádek popisující sken, B – dékla vazby (bond), A – úhel (angle), D – dihedrální úhel. V našem případě jsme definovali úhel A 2 1 3 S 11 10., kde 2 1 3 jsou atomy, které úhel svírají, S 11 10. indikuje sken s jedenácti kroky, každý po 10◦ . 18
3.4
Mezimolekulové komplexy vázané slabými interakcemi
Slabými interakcemi rozumíme interakce nižších energií, než je disociační energie kovalentní vazby, proto se jim někdy říká nekovalentní nebo také van der Waalsovy interakce (vdW).
Uplatnění slabých interakcí v přírodě • Existence vody jako kapaliny je způsobena pouze slabými interakcemi (vodíkové vazby). • 3D struktura DNA je stabilizována vodíkovými vazbami a patrovou (stackovou) interakcí mezi bázemi nukleových kyselin, to dělá molekulu DNA dostatečně pevnou na to, aby byla stabilní. Zároveň musí být síla, která drží báze nukleových kyselin dostatečně malá na to, aby se mohla molekula DNA při replikaci rozplést. • Slabé interakce stabilizují sekundární, terciální a kvarterní strukturu proteinů. Jak je vidět z popisu některých vdW komplexů, mají tyto útvary velmi specifické chování, proto je třeba při výpočtu jejich vlastností brát tyto zvláštnosti v úvahu. Výpočet interakční energie mezi molekulami tvořícími vdW komplex je „ústředním tématemÿ výpočetní chemie. Nelze se omezit jen na tuto problematiku, protozě mnohem důležitější se jeví změna Gibbsovy energie formace komplexu. Při interpretaci experimentálních měření (především gas–phase) můžeme přispět výpočty infračervených a Ramanových spekter, rotačních konstant, dipólových a ostatních multipólových momentů, v některých případech i UV/VIS spekter. Při těchto výpočtech je třeba používat metody zahrnující dostatečnou míru elektronové korelace (MP2, CCSD(T), FullCI, . . . ). Při ignorování bázové superpoziční chyby (BSSE) můžeme dostat přeceněné interakční energie, pokud korekci této chyby zahrneme i do optimalizace nebo výpočtu frekvencí, dosáhneme rovněž lepších výsledků, které lze lépe porovnat s experimenty. Cvičení Na výpočtu interakční energie dimeru molekuly vody vázaného vodíkovou vazbou si ukážeme velikost BSSE, vliv BSSE na optimalizaci a infračervená spektra. 1. Optimalizujte dimer vody běžným způsobem a se započtením korekce na BSSE, výsledné geometrie porovnejte. 19
2. Vypočtěte interakční energii dimeru vody se zahrnutím korekce BSSE a bez, výsledky porovnejte. Diskutujte, jak se mění velikost BSSE v závislosti na velikosti použité báze. 3. Vypočtěte frekvence dimeru s korekcí na BSSE a bez korekce, pro pozorování červeného posunu ve spektru vypočtěte rovněž frekvence monomeru. Určete velikost červeného posunu. Diskutujte vliv BSSE korekce na infračervená spektra. 4. Vypočtěte ∆G formace dimeru vody.
Water IR spectrum MP2/6-31G(d) 400 monomer dimer dimer BSSE
INTENSITY
300
200
100
0
0
1000
2000 1/cm
3000
4000
Obrázek 3.7: Infračervené spektrum vody – monomeru, dimeru s BSSE a bez BSSE korekce (MP2/6-31G(d)).
3.5
Výpočet frekvencí a termodynamických veličin
Co můžeme zjistit z výpočtu frekvencí? • Predikce infračervených a Ramanových spekter molekul (frekvence a intenzity). • Výpočet silových konstant pro použití v následující optimalizaci. 20
• Výpočet ZPVE (Zero-Point Vibration Energy), korekcí na vnitřní energii a další termodynamické veličiny jako entalpie (H) a entropie (S). Výpočty energie a optimalizace geometrie ignorují vibrace molekul, tímto jsou tyto výpočty idealizované z pohledu reálného pohybu jader. Toto zjednodušení vychází z již zmíněné Born-Oppenheimerovy aproximace v kapitole ??. V reálu jsou jádra v molekule v neustálém pohybu. V rovnovážněm stavu jsou tyto vibrace pravidelné a předpovídatelné, každá molekula může být tak identifikována na základě jejího charakteristického spektra. Gaussian umožňuje výpočet vibračního spektra molekul v jejich základním, ale i excitovaném stavu. Kromě predikce intezit a poloh (frekvencí) spektrálních čar je program schopen popsat rovněž pohyby, které molekula podstupuje při jejich normálních vibračních módech. Molekulární frekvence závisí na druhých derivacích energie podle pozic jader v souřadnicích. Gaussian je rovněž schopen predikce vlastností molekul závislých na druhých a vyšších derivacích energie jako jsou polarizability a hyperpolarizability, ty závisejí na druhých derivacích s ohledem na elektrické pole, jejich hodnoty jsou automaticky zahrnuty v každém výpočtu frekvencí Hartreeho Fockovou metodou. Druhé derivace fekvence lze počítat dvěma způsoby a to analyticky pro metody HF, MP2, DFT a CASSCF nebo numericky, prakticky pro všechny ostatní metody. Numerický výpočet frekvencí je ale více náročný na výpočetní čas. Jak již bylo zmíněno v kapitole 2.5 lze pomocí aparátu statistické termodynamiky získat z ab initio výpočtů i termodynamické veličiny. Pro výpočet termodynamiky a frekvencí je potřeba mít dobře optimalizovanou geometrii systému, který hodláme studovat. To znamená, že všechny první derivace energie podle souřadnic musí být nulové (minimum, tranzitní stav nebo vyšší řády sedlových bodů). Proto obvykle před výpočet frekvencí řadíme optimalizaci struktury použitím kombinace hesel OPT a FREQ v route části vstupního souboru. Není vhodné počítat optimalizaci a frekvence jinou bází nebo ještě húře jinou metodou, výsledky pak mohou být zkresleny.
3.5.1
Výpočet frekvencí
Cvičení 1. Zadejte výpočet optimalizace a frekvencí pro molekulu formaldehydu Hartreeho-Fockovou metodou s bází 6-31G(d). %chk=form-opt-freq.chk %mem=500MB %nproc=1 # RHF/6-31G(d) OPT FREQ formaldehyde
21
0 1 C O H H
0.000000 0.000000 0.934017 -0.934017
0.000000 0.000000 0.000000 0.000000
-0.009504 1.201084 -0.591677 -0.591677
Výpočet frekvencí začíná výpočtem energie optimalizované struktury. Poté následuje vlastní výpočet frekvencí, predikce intenzit a Ramanových depolarizačních poměrů a rozptylu pro každou spektrální čáru: 1 B1 Frequencies -- 1334.7569 Red. masses -1.3684 Frc consts -1.4364 IR Inten -0.3791 Raman Activ -0.7551 Depolar (P) -0.7500 Depolar (U) -0.8571 Atom AN X Y Z 1 6 0.17 0.00 0.00 2 8 -0.04 0.00 0.00 3 1 -0.70 0.00 0.00 4 1 -0.70 0.00 0.00 4 A1 Frequencies -- 2026.7287 Red. masses -7.2617 Frc consts -17.5744 IR Inten -149.8183 Raman Activ -8.1462 Depolar (P) -0.3282 Depolar (U) -0.4942 Atom AN X Y Z 1 6 0.00 0.00 0.58 2 8 0.00 0.00 -0.41 3 1 0.00 -0.46 -0.19 4 1 0.00 0.46 -0.19
2 B2 1383.1886 1.3437 1.5147 23.0833 4.5181 0.7500 0.8571 X Y Z 0.00 0.15 0.00 0.00 -0.08 0.00 0.00 -0.25 -0.65 0.00 -0.25 0.65 5 A1 3163.5031 1.0490 6.1850 49.4458 137.4683 0.1830 0.3093 X Y Z 0.00 0.00 0.06 0.00 0.00 0.00 0.00 0.61 -0.35 0.00 -0.61 -0.35
3 A1 1678.8678 1.1038 1.8330 8.6215 12.8176 0.5911 0.7430 X Y Z 0.00 0.00 0.00 0.00 0.00 0.08 0.00 -0.35 -0.61 0.00 0.35 -0.61 6 B2 3236.2917 1.1206 6.9152 135.0162 58.0557 0.7500 0.8571 X Y Z 0.00 0.10 0.00 0.00 0.00 0.00 0.00 -0.60 0.37 0.00 -0.60 -0.37
Počet normálních vibračních módů odpovídá počtu stupňů volnosti (3N −6). Vizualizujte IR a Ramanova spektra a jednotlivé vibrační módy v programu GaussView a porovnejte jejich pohyby s hodnotami ve výstupu. Pokuste se jednotlivé vibrační módy pojmenovat (např. symetrický stretch, bend apod.). Vzhledem ke známé systematické chybě v Hartreeho Fockově metodě, která je způsobena především zanedbáním elektronové korelace si můžeme dovolit frekvence násobit empirickým škálovacím faktorem (0.8929), který hodnoty teoretických frekvencí přiblíží experimentálním hodnotám. Vypočtené hodnoty budou i po škálování odlišné asi o 15% od experimentálních vzhledem k použité menší bázi pro konstrukci atomových orbitalů. Čtvrtý vibrační mód o frekvenci 1809.7 cm−1 (po škálování) odpovídá streči karbonylové skupiny, který bývá často v IR spektrech využíván pro její detekci. Mód má experimentální hodnotu 1746 cm−1 , ke které bychom se více přiblížili použitím větší báze. 22
3.5.2
Termodynamika
Implicitně je nastaven výpočet termodynamiky na tlak 1 atm a teplotu 298.15 K. Výpočet termodynamiky v Gaussianu probíhá při každém výpočtu frekvencí. Termodynamické veličiny můžeme najít ve výstupu aplikace pod heslem Thermochemistry. ------------------- Thermochemistry ------------------Temperature 298.150 Kelvin. Pressure 1.00000 Atm. Atom 1 has atomic number 6 and mass 12.00000 Atom 2 has atomic number 8 and mass 15.99491 Atom 3 has atomic number 1 and mass 1.00783 Atom 4 has atomic number 1 and mass 1.00783 Molecular mass: 30.01056 amu.
Dále ve výstupu zde najdeme symetrii rotace, vibrační a rotační teploty, rotační konstanty a velikost ZPVE (energie nulového vibračního módu), vyjádřené vztahem 12 hν. Imaginární frekvence, se do výpočtu ZPVE nezahrnuje. Rotational Rotational Rotational Zero-point
symmetry number 2. temperatures (Kelvin) constants (GHZ): vibrational energy
Vibrational temperatures: (Kelvin) Zero-point correction=
14.08085 1.92668 1.69478 293.39764 40.14554 35.31359 76700.6 (Joules/Mol) 18.33188 (Kcal/Mol) 1920.42 1990.10 2415.51 2916.01 4551.57 4656.30 0.029214 (Hartree/Particle)
Všechny následující řádky zahrnují zero point energy. První dává korekci na celkovou vnitřní energii, Etot = Et + Er + Ev + Ee Thermal correction to Energy=
0.032067
Další dva řádky jsou: Hcorr = Etot + kB T
(3.1)
a kde Stot
Gcorr = Hcorr − T Stot , = St + Sr + Sv + Se
Thermal correction to Enthalpy= Thermal correction to Gibbs Free Energy=
(3.2)
0.033011 0.008256
Následující čtyři řádky rozpočítávají celkovou energii molekuly po aplikaci různých korekcí. Vzhledem k tomu, že E jsem již použil pro označení vnitřní energie (internal thermal energy), pro celkovou elektronickou energii použiji symbol ε0 . Sum Sum Sum Sum
of of of of
electronic electronic electronic electronic
and and and and
zero-point Energies= thermal Energies= thermal Enthalpies= thermal Free Energies=
-113.837117 -113.834264 -113.833320 -113.858074
Sum of electronic and zero-point Energies = ε0 + εZP E Sum of electronic and thermal Energies = ε0 + Etot Sum of electronic and thermal Enthalpies = ε0 + Hcorr Sum of electronic and thermal Free Energies = ε0 + Gcorr Následující tabulka rozepisuje příspěvky k vnitřní energii (Etot ), tepelné kapacitě ((Ctot ) a entropii ((Stot ): 23
E (Thermal) KCal/Mol 20.122 0.000 0.889 0.889 18.345
Total Electronic Translational Rotational Vibrational
CV Cal/Mol-Kelvin 6.256 0.000 2.981 2.981 0.294
S Cal/Mol-Kelvin 52.100 0.000 36.130 15.921 0.050
Poslední tabulka vypisuje konkrétní příspěvky k partiční funkci. Total Bot Total V=0 Vib (Bot) Vib (V=0) Electronic Translational Rotational
3.5.3
0.159345D-03 0.436200D+10 0.366482D-13 0.100322D+01 0.100000D+01 0.646199D+07 0.672854D+03
-3.797661 9.639685 -13.435948 0.001398 0.000000 6.810366 2.827921
-8.744437 22.196196 -30.937413 0.003219 0.000000 15.681448 6.511529
Kinetika chemických rekcí
Výsledky z výpočtu termodynamických veličin lze rovněž použít ke zkoumání kinetiky reakcí. Vezmeme-li v úvahu teorii tranzitního stavu (TS), pak budeme potřebovat spočítat geometrie a frekvence reaktantů, produktů a tranzitního stavu což bude nejobtížnější úkol. V teorii tranzitního stavu bereme v úvahu, že reaktanty vytvářejí nestabilní aktivovaný komplex ≈ tranzitní stav, který se samovolně rozkládá za vzniku produktu. Rychlostní konstanta rozkladu TS je dána kB T /h. Předpokládejme tedy, že běžná bimolekulární rekce probíhá způsobem: k1
k2
A + B AB ‡ → produkty, k−1
(3.3)
pak K‡ =
[AB ‡ ] [A][B]
k2 =
kB T . h
(3.4)
Rychlost reakce je v = k2 [AB ‡ ] = k2 K ‡ [A][B] = k[A][B],
(3.5)
kde k je rychlostní konstanta, tu je možné změřit i experimentálně. Z termochemického hlediska: ∆S ‡ =
∆H ‡ − ∆G‡ T
∆G‡ = ∆H ‡ − T ∆S ‡ = −RT ln K ‡ cn−1 , 0 kde n = 2 pro bimolekulární reakci. Pak při teplotě 298 K platí kB T ∆G‡ 6.2125 × 1012 ∆G‡ k= exp − = exp − , hc0 RT c0 RT
(3.6) (3.7)
(3.8)
kde c0 v jednotkách mol/l je faktor pro bimolekulární reakci. Při 1 atm, ideální plyn c0 = P/RT . 24
Kapitola 4
Úlohy z molekulové mechaniky 4.1
Vliv solventu na strukturu peptidu
Úkol: Zjistěte, jak se liší struktura tripeptidu THR–ILE–CYS (konformace β–skládaného listu), provede-li se optimalizace jeho geometrie ve vakuu a v prostředí explicitních vod. Použijte aplikaci HyperChem. Postup: 1. Zvolte silové pole AMBER: Setup → Molecular Mechanics. Zaškrtněte „AMBERÿ a v možnostech nastavení „Optionsÿ ponechejte vše, jak je přednastavené až na „Dielectric (Epsilon)ÿ, které zvolte „Distance dependentÿ. Nechte zahrnuty všechny energetické členy, což zkontrolujte v nabídce „Componentsÿ, kde musí být včechny energetické termy označené. Během editace struktury nezapomeňte opakovaně přepočítat atomové typy: Build → Calculate Types. 2. Vytvořte startovní strukturu tripeptidu THR–ILE–CYS. Níže jsou uvedeny dva nejschůdnější způsoby, jak k úkolu přistupovat: • „Naklikejteÿ peptid přímo v aplikaci HyperChem: Databases → Amino Acids. Nezapomeňte, že máte pracovat s peptidem v konformaci β–skládaného listu! • Pro zkušenější uživatele je dalším řešením stažení libovolného proteinu (PDB souboru) z databáze (např. www.pdb.org) a z něj „vystříhnoutÿ výše uvedenou aminokyselinovou sekvenci. Podmínkou samozřejmě je, že protein musí danou sekvenci v konformaci β–skládaného listu obsahovat. Úpravu staženého PDB souboru proveďte buď v textovém editoru, nebo v aplikacích jako jsou PyMOL či VMD. 25
3. Máte-li danou sekvenci postavenou/načtenou do aplikace HyperChem, vytvořte amfion tripeptidu (N-konec (THR) bude mít —NH+ 3 skupinu a C-konec (CYS) —COO− skupinu): Databases → Make Zwitterion. 4. Pokud máte amfion, uložte výchozí strukturu tripeptidu: File → Save As. Použijte formát PDB a nezapomeňte zaškrtnout položky „Hydrogensÿ a „Connectivityÿ. Soubor nazvěte „tri-start.pdbÿ. 5. Vypočítejte „Single Pointÿ energii a výsledek zapište do tabulky: Compute → Single Point. 6. Proveďte optimalizaci geometrie tripeptidu ve vakuu: Compute → Geometry Optimization. Zvolte optimalizační metodu „Conjugate gradientÿ a kriteria pro ukončení optimalizace nastavte následovně: RMS = 0.001 kcal/mol · ˚ A a 1000 optimalizačních kroků (1000 „maximum cyclesÿ). 7. Uložte strukturu zoptimalizovaného tripeptidu. Opět jako v předchozím případě. Soubor nazvěte „tri-optgas.pdbÿ. 8. Vytvořte v aplikaci HyperChem nový soubor: File → New. „Naklikejteÿ znovu strukturu tripeptidu, vyvořte amfion a silové pole nastavte opět na AMBER. Všechna nastavení ponechejte opět jako v předchozím případě až na položku „Dielectric (Epsilon)ÿ, kde označte „Constantÿ. 9. Vložte molekulu tripeptidu do periodického „boxuÿ vyplněného explicitními molekulami vody: Setup → Periodic Box. Ponechejte přednastavené dimenze kubického „boxuÿ a minimální vzdálenost mezi solutem a solventem nastavte na 2.3 ˚ A. 10. Proveďte optimalizaci geometrie tripeptidu v explicitním boxu, na molekuly vody aplikujte „constrainÿ. Metoda „Conjugate gradientÿ, RMS = 0.001 kcal/mol · ˚ A a (1000 optimalizačních kroků 1000 „maximum cyclesÿ). 11. Po skončení optimalizace odstraňte všechny explicitní vody. Postupujte následovně: Select → Molecules a pak vyselektujte molekulu solutu a následně proveďte Select → Complement Selection a smažte výběr. 12. Vypočítejte „Single Pointÿ energii struktury tripeptidu optimalizovaného v „boxuÿ vod. Nezapomeňte přepnout položku „Dielectric (Epsilon)ÿ na „Distance dependentÿ. Ostatní nastavení týkající se silového pole AMBER ponechejte. 13. Uložte strukturu tripeptidu. Soubor pojmenujte jako „tri-optwat.pdbÿ. 26
14. Otevřete soubory tri-optgas.pdb a tri-optwat.pdb v aplikaci VMD, proveďte přeložení obou systémů a spočítejte vzájemné RMS. V aplikaci VMD postupujte následovně: Jsou-li oba tripeptidy načteny, zvolte v okně VMD Main Extension → Analysis → RMSD Calculator. V okně „RMSD Toolÿ odznačte položku „Backbone onlyÿ a do pole nahoře napište „allÿ. Klepněte na tlačítko „Alignÿ, tím se provede přeložení obou systémů a tlačítkem „RMSDÿ proveďte vápočet hodnoty RMS. Vyhodnocení: Sestavte tabulku, do které uveďte energie příslušných konformerů tripeptidu THR–ILE–CYS. Do tabulky zapište rovněž RMS mezi strukturou tripeptidu získanou optimalizací ve vakuu a v prostředí explicitních vod. V aplikaci VMD vytvořte obrázky příslušných konformerů tripeptidu.
4.2
Sken povrchu potenciální energie podle jednoho parametru
Úkol 1: Spočítejte energetický profil rotace torze H—O—O—H v molekule peroxidu vodíku. Použijte aplikaci HyperChem. Postup: 1. Nakreslete strukturu peroxidu vodíku H—O—O—H (obrázek 2.2). Stačí spojit dva atomy kyslíku kovalentní vazbou a pak přejít na Build → Add H and Model Build. Dostanete tak strukturu peroxidu vodíku v nezákrytové konformaci s torzním úhlem H—O—O—H o velikosti 180◦ . 2. Nastavte silové pole MM+: Setup → Molecular Mechanics. Ve zobrazeném okně vyberte pole MM+, zahrňte všechny složky energie a v okně „Optionsÿ ponechejte vše implicitně až na nabídku „Electrostaticsÿ, kde označte „Bond dipolesÿ. 3. Spočítejte „Single Pointÿ energii pro strukturu v zákrytové konformaci (torze H—O—O—H = 180◦ ) a pak postupně měňte příslušný torzní úhel od 180 do −180◦ v krocích po −10◦ . Dihedrální úhel se nastaví tak, že se označí všechny atomy molekuly a pak se ručně nastaví velikost torze v okně „Set Bond Torsionÿ (Edit → Set Bond Torsion). Celkem provedete 37 „Single Pointÿ výpočtů. Pozn.: Aplikaci HyperChem je možné ovládat prostřednictvím maker tabulkového procesoru Excel. Právě pro případ výpočtu energetického profilu torzního úhlu je vhodné vytvořit jednoduché makro a práci si tak značně zjednodušit. V příloze skript je makro, které lze použít pro výpočet energetického profilu 27
torze H—O—O—H peroxidu vodíku. Vyhodnocení: Výsledky zapište do tabulky a sestrojte grafickou závislost MM+ energie na velikosti dihedrálního úhlu H—O—O—H. Porovnejte výsledky s kvantověmechanickým přístupem použitým v kapitole 9.2.1.
Úkol 2: Spočítejte energetický profil rotace torze Cl—C—C—Cl v molekule 1,2-dichlorethanu. Použijte aplikaci HyperChem. Postup: Postupujte naprosto stejně jako v úloze 1. Při výběru torze Cl—C—C—Cl je třeba postupně označovat jednotlivé atomy. Nesmí se označit celá molekula, protože pro nastavení velikosti torzního úhlu v okně „Set Bond Torsionÿ je potřeba mít vyselektované pouze čtyři atomy. Ve vašem případ jsou to atomy Cl, C, C a Cl. Vyhodnocení: Data vyneste do grafu závislosti MM+ energie na velikosti torze Cl—C— C—Cl. Srovnejte profil torze Cl—C—C—Cl s profilem dihedrálního úhlu H—O—O—H. Porovnejte energetické rozdíly jednotlivých konformací 1,2dichlorethanu s experimentem.
28
Kapitola 5
Úlohy z molekulové simulace 5.1
Simulované žíhání
Chceme-li získat globální energetické minimum studované molekuly nebo mít jistotu, že nás použitý optimalizační algoritmus do globálního minima skutečně zavedl, je výhodné nasadit molekulovou simulaci. Prostřednictvím molekulové simulace za zvýšené teploty dokážeme efektivně prohledávat konfigurační prostor systému (dostávat se přes energetické bariéry na PES), což vede k nalezení globálního minima. Během molekulové simulace se každý zvolený krok ukládají strukturní data (souřadnice molekuly), která jsou následně reoptimalizována. Konformaci, které z celého klastru uložených souřadnic molekuly přísluší nejnižší energie, označujeme jako globální minimum. Dalším řešením, jak získat potencionální globální minimum studovaného systému, je postupné snižování teploty simulace, tzv. „slow coolingÿ.
5.1.1
Nalezení globálního minima acetonu
Úkol: Nalezněte globální minimum acetonu metodou simulovaného žíhání v aplikaci HyperChem. 1. Nakreslete strukturu acetonu (CH3 —CO—CH3 ). Dvojnou vazbu mezi atomy kysíku a uhlíku vytvoříte tak, že na jednoduchou vazbu mezi danými atomy jednou klepnete. Pro získání věrohodnějšího modelu acetonu proveďte následující akci: Build → Add H and Model Build. 2. Vyberte silové pole MM+ (Setup → Molecular Mechanics Force Field ). Elektrostatickou složku celkové energie spočítejte přes parciální atomové náboje (v nabídce „Electrostaticsÿ označte „Bond dipolesÿ) a zahrňte všechny energetické složky. Vypočítejte „Single Pointÿ energii acetonu (silové pole MM+). 3. Modifikujte úhel C—C—C z výchozí hodnoty −120◦ na −140◦ . 29
4. Vyberte a pojmenujte strukturní parametry, které chcete monitorovat během simulace (Select → Name Selection). Určitě vyselektujte a pojmenujte (např. ccc) úhel C—C—C, abyste mohli sledovat jak se mění jeho velikost během simulace. 5. Nastavte simulaci tak, aby se provádělo simulované žíhání (obrázek 5.1). Heat time: Run time: Cool time: Step size: Starting temperature: Simulation temperature: Final temperature: Temperature step: In vacuo Constant temperature Data collection period : Screen refresh period : Playback : Restart:
2 ps 5 ps 2 ps 0.002 ps 0K 500 K 0K 10 K on off 1 time steps 1 data steps off off
Obrázek 5.1: Nastavení simulace pro provedení simulovaného žíhání v aplikaci HyperChem. 6. Nastavte cestu a název souboru pro kolekci strukturních dat (ukládání souřadnic): Molecular Dynamics Options → Snapshots. 7. Nastavte cestu a název souboru pro ukládání energetických termů a strukturních parametrů: Molecular Dynamics Options → Averages. 30
Do pole „Avg. and graphÿ přesuňte selekci „cccÿ a položku EKIN (kinetickou energii systému). 8. Odstartujte simulované žíhání: Molecular Dynamics Options → Proceed. 9. Po skončení MD simulace se podívejte na střední hodnotu sledovaného torzního úhlu C—C—C (selekce „cccÿ). Pro nalezení aritmetického průměru strukturního parametru je nutno přejít na Compute → Molecular Dynamics → Averages. 10. Spočítejte „Single Pointÿ energii acetonu po simulovaném žíhání. Opět ponechejte silové pole MM+ a výpočet elektrostatického energetického příspěvku přes vazebné dipóly. Vyhodnocení: Porovnejte energii struktury acetonu, ze které se startovalo simulované žíhání (struktura s C—C—C úhlem o velikosti −140◦ ) s energií acetonu získaného simulovaným žíháním. Zhotovte obrázky obou struktur. V aplikaci VMD proveďte přeložení příslušných molekul a spočítejte velikost RMS.
5.1.2
Nalezení globálního minima peroxidu vodíku
Úkol: Najděte globální minimum peroxidu vodíku pomocí simulovaného žíhání, máte-li jako startovní strukturu jeho zákrytovou konformaci (velikost torzního úhlu je 0◦ ). Opět jako v předchozí úloze pracujte v prostředí aplikace HyperChem.
31
Kapitola 6
Úlohy z bioinformatiky 6.1 www.pdb.org
Databáze proteinů a nukleových kyselin RCSB
Proteinová databanka obsahuje strukturu desítek tisíců biomolekul. Informace o struktuře je nejčastěji získána rentgenokrystalografickou analýzou (RTG) nebo nukleární magnetickou rezonancí (NMR). Jedná se o redundantní databázi.
6.1.1
Stručná historie
Databáze byla založena v roce 1971 Meyerem a Hamiltonem. V roce 1998 se stala součástí centra RCSB (Research Collabolatory for Structural Bioinformatics). Struktury biomolekul jsou volně stažitelné a lze je vizualizovat celou řadou programů (VMD, RasMol, PyMOL a další)
6.1.2
Práce s databází
Ovládání databáze je velmi intuitivní. Každé molekule přísluší čtyřmístný kód, který funguje jako adresa pro vyhledání systému v databázi. Celkem existují dvě základní možnosti, jak získat informace o biomolekule: • Zadáním kódu přímo na domovské stránce databanky (obrázek 6.1) nebo v databázi, která je s proteinovou databankou propojena (PDBsum, . . .). • Prostřednictvím vizualizačního programu (PyMOL, VMD, Chimera, MOE), který uživateli umožní vizualizovat biomolekulu a následně provést uložení strukturních dat. Neznáme-li PDB kód biomolekuly, lze vyhledávat podle názvu systému, autora, apod. (obrázek 6.1). Tento způsob hledání není tak jednoznačný, a proto je efektivnější použít rozšířeného hledání (obrázek 6.2). Uživatel zde může vyhledávat podle strukturních parametrů, které si předem nastaví 32
(identifikátor systému, počet a procentuální zastoupení sekundárních struktur, počet aminokyselin v biomolekule, metody použité pro shromáždění strukturních dat a další).
Obrázek 6.1: Hlavní stránka proteinové databanky. (a) pole pro vložení čtyřmístného PDB kódu, (b) odkaz sloužící pro přepnutí do režimu rozšířeného hledání. Příklad Na+ /K+ ATPasa - biologická funkce a struktura Souřadnice Na+ /K+ ATPasy byly získány krystalograficky s rozlišením 3, 5 ˚ A. + Jedná se o integrální membránový protein, jehož funkcí je transport tří Na iotů ven z buňky a import dvou K+ iontů do cytosolu. Skládá se ze tří podejdnotek, α, β a γ. V současnosti je kompletně vyřešena struktura celé α-podjednotky, která je dále složena ze tří domén (A, N a P). Cvičení 1. Prostřednictvím programů VMD nebo PyMOL určete přesný počet proteinů ve velké ribozomální podjednotce z „Haloarcula marismortui ÿ pod PDB kódem 1JJ2. Zjistěte, kolik reziduí má největší/nejdelší protein a kolik nejmenší/nejkratší. 2. Vyhledejte v databázi protein s největším počtem disulfidových můstků. Podmínkou je, že systém nesmí mít více jak sto aminokyselinových 33
Obrázek 6.2: Rozšířené vyhledávání v databázi PDB podle pěti parametrů: (a) typu experimentální metody, (b) typu biomolekuly obsažené ve struktuře, (c) počtu disulfidových můstků, (d) počtu reziduí, (e) počtu řetězců v systému. zbytků a jeden řetězec. Zjistěte, jakou biologickou funkci protein vykonává, jakou metodou byla strukturní data získána a vizualizujte disulfidové můstky. Diskutujte, do jaké míry jsou jednotlivé části systému disulfidovými můstky stabilizovány [2QSK, pět disulfidových můstků]. 3. Vizualizujte ve VMD protein s PDB kódem 1PGA a pomocí nástroje „Ramachandran Plotÿ zobrazte distribuci jednotlivých φ a ψ torzí. Rozhodněte, do jaké oblasti Ramachandranova diagramu spadají phi a psi torze α-šroubovice a β-skládaného listu. Porovnejte s hodnotami v literatuře. 4. Proteinový inhibitor p18INK4c (p18) kontroluje buněčný cyklus prostřednictvím specifické inhibice cyklin dependentních kináz. Podařilo se získat jak strukturu vázaného p18 (PDB ID 1G3N), tak strukturu nevázaného p18 (PDB ID 1IHB). Proveďte v programu PyMOL přeložení struktury vázané a nevázané p18 (přes Cα atomy) a diskutujte, jak se změní struktura p18, interaguje-li s kinázou.
34
6.2
Visual Molecular Dynamics
Visual Molecular Dynamics (VMD) je volně šířitelný program pro zobrazení struktury a dynamiky biomolekul. Součástí aplikace je celá řada nástrojů umožňujících analýzu struktury a dynamiky zejména proteinů. Nespornou výhodou VMD je podpora celé řady operačních systémů (MacOS-X, Unix, a Windows). V následujícím textu se budeme zabývat stručným popisem základních funkcí programu VMD ve verzi 1.8.6 používaným v prostředí Windows. Po startu aplikace (kliknutím na soubor vmd.exe nebo příslušného zástupce) se před uživatelem objeví tři okna – konzole aplikace VMD, „VMD Mainÿ nabídka umožňující interaktivní ovládání a okno „VMD 1.8.6 OpenGL Displayÿ vykreslující danou molekulu podle nastavených požadavků. 1. VMD konzole Prakticky celou aplikaci lze ovládat přímo z konzole (obrázek 6.3). Vyžaduje to však poměrně značné zkušenosti s programem VMD, a proto se v následujícím výkladu omezíme poze na interaktivní ovládání celé aplikace.
Obrázek 6.3: VMD konzole po spuštění aplikace VMD. V konzoli je vypsaná verze aplikace, primární citace na program VMD, konfigurace PC a grafické karty a cesta k doplňkovým aplikacím („pluginsÿ). Úplně na konci textu je příkazový řádek umožňující relativně rychlou obsluhu programu. 2. Nabídka „VMD Mainÿ V prvé řadě je zde uživateli umožněno načtení celé palety souborů uchovávajících strukturní data molekul (PDB, XYZ, AMBER Coordinates apod.). Molekulu lze načíst do aplikace dvojí cestou: - Přímo z pevného disku či připojeného datového média – File → New Molecule → Browse → Determine file type: → Load 35
- Z databází PDB nebo STING – Extensions → Data → PDB Database Query nebo STING Database Query Máme-li načtený systém (např. ve formátu PDB), zobrazí se v nabídce „VMD Mainÿ na tyrkysovém poli řádek informující o základních vlastnostech systému (VMD ID; sada čtyř písmen „T A D Fÿ pro snadnou manipulaci se systémem v případě, kdy je jich načteno více než jeden; název molekuly; počet atomů systému; počet snímků a položka „Volÿ podávající informace o počtu načtených volumetrických dat). Pro lepší představu viz obrázek 6.4.
Obrázek 6.4: Nabídka „VMD Mainÿ s načtenou strukturou B1 domény proteinu G (PDB ID 1PGA). Z okna „VMD Mainÿ lze v podstatě ovládat interaktivně veškeré funkce programu VMD. K nejdůležitějším částem přístupným z hlavní nabídky „VMD Mainÿ patří okno „Graphical Representationsÿ zajišťující plnou kontorlu nad způsobem zobrazení systému v okně „VMD . . . Displayÿ. V nabídce „Graphical Representationsÿ je k dispozici celá řada přednastavených režimů vykreslování systému. Zajímá-li nás jaké oblasti biomolekuly jsou přístupné molekulám solventu je vhodné vybrat vykreslování „Surfÿ. Chceme-li mít informaci o kompaktnosti proteinu, zvolíme metodu vykreslování „VDWÿ. Síť vodíkových vazeb lze jednoduše zobrazit selekcí „HBondsÿ. Metoda „Traceÿ vyselektuje Cα atomy proteinu a vzájemně je pospojuje virtuální vazbou (tzv. „α − traceÿ). Informaci o topologii systému a vzájemné orientaci proteinových sekundárních struktur poskytne režim vykreslení „Ribbonsÿ nebo „Cartoonÿ. Každému způsobu vykreslení lze přiřadit barvu přímo z palety nebo přednastavených barevných filtrů (např. „Nameÿ, „Typeÿ, „Chainÿ a další). Uživateli je v „Graphical Representationsÿ umožněno nastavit i texturu vykreslovacích režimů (viz Material ). Na obrázku 6.5 je nabídka „Graphical Representationsÿ ve 36
spojení s proteinem pod PDB ID 1PGA (B1 doména proteinu G). Nastavení je provedeno tak, abychom dostali vykreslení NewCartoon v barevném rozlišení sekundárních struktur a glutamové kyseliny společně s lysiny zobrazené na úrovni jednotlivých atomů (režim „CPKÿ) zbarvené podle druhu atomu (barevný filtr „Nameÿ). Výsledek nastavení reprezentace proteinu je na obrázku 6.6.
Obrázek 6.5: Nabídka „Graphical Representationsÿ pro B1 doménu proteinu G. Na obrázku vlevo je reprezentace s vykreslením „NewCartoonÿ a barvě „Structureÿ aplikovaná na všechny atomy. Na obrázku vpravo je všem glutamovým kyselinám a lysinům proteinu přiřazeno vykreslení „CPKÿ v barevném módu „Nameÿ.
37
3. Okno„VMD 1.8.6 OpenGL Displayÿ Zde lze manipulovat se zobrazeným systémem ve smyslu translačních a rotačních pohybů. Stisknutím klávesy „Rÿ se přepínáme do rotačního módu, stisknutím klávesy „Tÿ je aktivován translační mód a stisknutím klávesy „Sÿ se přepneme do režimu zmenšování/zvětšování („zoomÿ). Uživatel má rovněž možnost prostřednictvím numerických klásves „1-4ÿ aktivovat režim identifikace atomů, měření délky vazeb, vazebných úhlů a torzí. Na obrázku 6.6 je globulární B1 doména proteinu G (56 reziduí, α/β protein – β-list překrytý α-šroubovicí) vykreslovaná stylem „NewCartoonÿ v obarvení „Structureÿ s rozlišením jednotlivých atomů zobrazením postranních řetězců glutamové kyseliny a lysinu (metoda vykreslení „CPKÿ v barvě „Nameÿ.
Obrázek 6.6: Terciární struktura B1 domény proteinu G (PDB ID 1PGA). Reprezentace proteinu podle nastavení uvedeného v „Graphical Representationsÿ. Užitečné tipy a triky: Změna barvy pozadí Implicitně má VMD v okně „VMD . . . Displayÿ nastavenou černou barvu pozadí. Změnu na jakoukoliv jinou barvu lze provést následujícím způsobem – Graphics → Colors a v levém sloupci Categories vybrat Display, pak v prostředním sloupci Names zvolit Background a v pravém sloupci Colors aktivovat příslušnou barvu (třeba bílou). Barvu pozadí v okně „VMD . . . Displayÿ lze efektivněji změnit z konzole aplikace VMD příkazem color Display Background white
38
Vypnutí os kartézského souřadného systému Zobrazení os XYZ souřadného systému často zavazí grafickému zpracování vizualizovaného objektu (zhotovení rastrového obrázku, „renderingÿ). Zrušení os XYZ souřadného systému se provádí následovně – Display → Axes → Off. Trasování Potřebujeme-li získat obrázek systému v mimořádné kvalitě je vhodné provést trasování. Program VMD obsahuje modul, který je schopen sestavit popis scény pro kompilátory fotorealistických obrázků, jako jsou například POV-Ray, RenderMan, VRML a další. K modulu se dostaneme touto cestou – File → Render. Chceme-li dosáhnout ještě vyšší kvality zobrazení, nastavíme v nabídce „Graphical Representationsÿ rozlišení příslušného vykreslovacího stylu na co nejvyšší. Výrazné zvýšení hodnoty rozlišení jednotlivých stylů se nicméně nedoporučuje provádět na pomalejších PC!
39
Část III
Dodatky
40
Kapitola 7
Reprezentace molekul 7.1 7.1.1
Formáty souborů používané při reprezentaci molekul Formát Z-matice
Z-matice definuje vzájemnou orientaci atomů v molekule podle jejich přirozených vztahů, jako jsou vazebné délky, úhly a dihedrální úhly. Atom A
Vzdálenost
Úhel
B
B–A
C
C–B
C–B–A
D
D–C
D–C–B
Dihedr. úhel
D–C–B–A
První atom je umístěn do počátku. Pro druhý atom je třeba definovat vzdálenost od prvního atomu. Pro třetí atom je třeba přidat další upřesňující parametr — úhel. Polohu čtvrtého a každého následujícího atomu je třeba navíc definovat dihedrálním úhlem.
Tabulka 7.1: Z-matice O H 1 r1 H 1 r1 2 a1 r1 0.9584 a1 104.45
41
Obrázek 7.1: Geometrie molekuly vody.
7.1.2
Formát XYZ
Ve formátu XYZ jsou pozice atomů definovány v tzv. kartézském souředném systému, pomocí tří os X, Y a Z.
Obrázek 7.2: Kartézský souřadný systém.
8 C C H Cl H Cl
-4.3266 -2.7890 -2.4540 -2.1405 -4.7316 -4.9011
2.1938 0.0159 2.2132 0.0734 1.8293 1.0481 1.1943 -1.2158 2.8275 0.8184 0.5366 0.2265 42
H H
-4.6617 -2.4301
7.1.3
2.5771 -0.9595 3.2445 -0.0579
Formát PDB
• používaný pro reprezentaci biomolekul (proteiny, nukleové kyseliny) • podporuje většina programů pro počítačové modelování (VMD, PyMOL, YASARA, Chimera, . . . ) • polohy atomů v XYZ souřadnicích (jednotka ˚ A) • velmi striktní formátování: – prvních šest sloupců vyhrazeno pro „Record nameÿ (viz ATOM, ANISOU, TER, HETATM, MODEL, REMARK, atd.) – souřadnice atomů standardních reziduí uvedeny u ATOM – proteiny číslovány od N konce a nukleové kyseliny od 5’ konce – souřadnice molekul vody, iontů, ligandů, atd. uvedeny obvykle u HETATM – obsahuje-li biomolekula více řetězců, pak je každý řetězec oddělen prostřednictvím TER – každý PDB soubor je zakončen END
Obrázek 7.3: Obecný formát souboru PDB.
43
ATOM ATOM ATOM
9730 9731 9732
O H1 H2
WAT WAT WAT
3085 3085 3085
11.143 11.935 11.436
14.755 14.260 15.434
17.359 17.151 17.966
Obrázek 7.4: Molekula vody (AMBER, TIP3P) vizualizovaná programem YASARA s použitím souřadnic v PDB formátu uvedenými nad obrázkem.
44
Kapitola 8
Povrchy potenciální energie
Obrázek 8.1: Příklad jednoduchého povrchu potenciální energie. Povrchy potenciální energie (PES) nám umožňují grafické přiblížení vztahu geometrie daného systému a jeho energie. Složitost — dimenzionalita povrchů potenciální energie se zvyšuje zároveň se vzrůstajícím počtem atomů v pozorovaném systému. Složitost daného PES nám určuje počet stupňů volnosti, viz. tabulka 8.1.
45
Translace (x,y,z) Rotace (x,y,z) Vibrace Celkem
Atom 3 0 0 3
Lineární molekula 3 2 3N − 5 3N
Nelineární molekula 3 3 3N − 6 3N
Tabulka 8.1: Stupně volnosti, kde N je počet atomů v molekule a x,y,z jsou souřednice v kartézském systému.
46
Část IV
Přílohy
47
Kapitola 9
Úlohy z kvantové chemie 9.1 9.1.1
Stavba molekul Prostředí aplikace GaussView
Vstup pro program Gaussian, cvičení 2 %chk=dicleth.chk %mem=500MB %nproc=1 # hf/3-21g geom=connectivity dihedral 0 degree 0 1 C H H C H H Cl Cl B1 B2 B3 B4 B5 B6 B7 A1 A2 A3 A4 A5 A6 D1 D2 D3 D4 D5
1 1 1 4 4 4 1
B1 B2 B3 B4 B5 B6 B7
2 3 1 1 1 4
1.06999992 1.06999991 1.54000000 1.06999944 1.06999992 1.76000019 1.76000019 109.47126292 109.47124129 109.47125022 109.47118427 109.47122992 109.47122992 -119.99999401 120.00001377 0.00000001 -120.00001771 -119.99996852
48
A1 A2 A3 A4 A5 A6
2 3 3 3 5
D1 D2 D3 D4 D5
1 2 1.0 3 1.0 4 1.0 8 1.0 2 3 4 5 1.0 6 1.0 7 1.0 5 6 7 8
Jak je vidět ze vstupu, na prvních řádcích je definováno jméno souboru checkpoint (%chk=dicleth.chk), dále kolik si program může alokovat paměti (%mem=500MB) a na kolika procesorech bude výpočet probíhat (%nproc=1). Řádek uvozený symbolem # je nejdůležitější částí vstupu a je v něm uvedeno co a jak se bude počítat. V našem případě budeme počítat single point (bez parametru) metodou HF a bází 3-21G. Parametr geom=conectivity znamená, že na konci uvedené geometrie molekuly v Z-matici jsou uvedeny informace o tom, které atomy jsou spolu propojeny vazbami. Nejedná se o povinný parametr. Tento vstup byl generován programem GaussView.
9.2
Sken povrchu potenciální energie
9.2.1
Prostředí aplikace GaussView
Vstup pro program Gaussian %chk=hooh.chk %mem=6MW %nproc=1 # scan pbepbe/3-21g geom=connectivity nosymm Title Card Required 0 1 O H 1 B1 O 1 B2 2 A1 H 3 B3 1 A2 2 D1 B1 B2 B3 A1 A2 D1
0.96000000 1.32000000 0.96000000 109.50000006 109.50000006 180.0 36 10.
1 2 1.0 3 1.0 2 3 4 1.0 4
Opět jako v předchozím případě je definováno jméno souboru checkpoint (%chk=hooh.chk) Kolik si program může alokovat paměti (%mem=6MW) a na kolika procesorech bude výpočet probíhat (%nproc=1). Vstup informuje 49
program o tom, že se bude provádět scan metodou pbepbe a bází 3-21G. Sken peroxidu vodíku bude probíhat v celkovém počtu 36 kroků (od 180 do −180◦ ) s šířkou kroku 10◦ . Vstup byl generován programem GaussView.
50
Kapitola 10
Úlohy z molekulové mechaniky 10.1
Sken povrchu potenciální energie
10.1.1
Prostředí aplikace HyperChem
Makro pro výpočet profilu torze H—O—O—H v programu HyperChem 1. Zkopírujte soubor „Plot.xlmÿ z domovského adresáře aplikace HyperChem do svého pracovního adresáře. Cesta je následující: \$HYPERHOME\Samples\Scripts\ 2. V aplikaci Excel soubor upravte tak, aby obsah korespondoval s obsahem v tabulce 10.1. Za $PATH dosaďte celou cestu k souboru „chem.exeÿ: C:\Hyper80\Program\chem.exe 3. V programu Excel povolte makra a vytvořte v otevřeném souboru „Plot.xlmÿ nový list. 4. Do prvního sloupce v novém listu vepiště hodnoty úhlu H—O—O—H, pro které chcete vypočítat „Single Pointÿ energii. V sloupci bude celkem 37 položek (od 180 do −180◦ s velikostí kroku 10◦ ). 5. Označte buňku s velikostí úhlu 180◦ a spustěte aplikaci HyperChem. 6. Nakreslete strukturu H—O—O—H a zvolte metodu výpočtu potenciální energie příslušné molekuly. 7. Pokud je vše jak v aplikaci Excel, tak HyperChem nachystáno, spustěte makro (Compute.Results → Spustit). 51
8. Postupně se budou ke každému úhlu dopisovat energie H—O—O—H. Makro je ukončeno, jakmile se mu do cesty dostane prázdná buňka. Control-R Channel
NewChan
Comput.Results =OpenFile() =KDYŽ(JE.CHYBHODN(Channel)) =NÁVRAT() =KONEC.KDYŽ() =VYKONAT(Channel;”[query-response-has-tag(no)]”) =VYKONAT(Channel;”[menu-select-select-all]”) =DOKUD(NE(JE.PRÁZDNÉ(VÝBĚR()))) =VYKONAT(Channel;”[set-bond-torsion(”&VÝBĚR()&”)]”) =VYKONAT(Channel;”[do-single-point]”) =VZOREC.MATICE(POŽADAVEK(Channel;”total-energy”);”rc[1]”) =VYBRAT(”r[1]c”) =DALŠÍ() =KANÁL.ZAVŘÍT(Channel) =NÁVRAT() OpenFIle =KANÁL.OTEVŘÍT(”HyperChem”;”System”) =KDYŽ(JE.CHYBHODN(NewChan)) =KDYŽ(JE.CHYBHODN(SPUSTIT.PROGRAM(”$PATH”;1))) =NÁVRAT(NewChan) =KONEC.KDYŽ() =NÁVRAT(KANÁL.OTEVŘÍT(”HyperChem”;”System”)) =KONEC.KDYŽ() =NÁVRAT(NewChan)
Tabulka 10.1: Makro aplikace Excel pro sken torze H—O—O—H v programu HyperChem.
52
