Tutaj zebraliśmy wszystkie informacje dotyczące modelowania molekularnego, oprogramowania, zagadnień, sposobu rozwiązywania problemów - czyli wszystko to co może się przydać każdemu chemikowi, farmaceucie, biologowi czy pasjonatowi.
Przeprowadzenie obliczeń, interpretacja wyników
Energia elektronowa układu jest zawsze wyliczona podczas każdych obliczeń, służy między innymi w algorytmach optymalizacyjnych, w których poszukiwana jest najmniejsza energia. Aby obliczyć energię wystarczy policzyć energię dla konkretnego układu. W tym celu stosujemy słowo kluczowe SP (Single point) w linii ze znakiem # np:
# SP B3LYP/6-31G(d)
Można także przeprowadzić optymalizację, lub dowolne obliczenia. Po wykonaniu obliczeń energię cząsteczki możemy sprawdzić w pliku wynikowym szukając linii[1] zawierających słowo SCF Done:
W metodzie HF będzie to:
SCF Done: E(RHF) = -985.139230617 A.U. after 3 cycles
W metodzie B3LYP:
SCF Done: E(RB+HF-LYP) = -155.033805085 A.U. after 7 cycles
W metodzie MP2:
SCF Done: E(RHF) = -985.127777421 A.U. after 8 cycles
Oraz kilka lini poniżej w linii zawierającej:
E2 = -0.1992260933D+01 EUMP2 = -0.98712003835432D+03
Wszystkie energie są energiami elektronowymi podawanymi w jednosce atomowej Hartree[2], w metodzie MP2 podana jest energia wyliczona metodą HF (SCF Done), poprawka (E2) uwzględniająca korelacje między elektronami oraz energia wyliczona metodą Hartree Focka z poprawką (EUMP2), np.: -0.98712003835432D+03[3]. Te same wartości energii co w ostatnich liniach z napisem SCF Done możemy znaleźć na końcu pliku, na przykład dla optymalizacji prowadzoną w ćwiczeniu z cząsteczkami etanolu metodą B3LYP:
1\1\GINC-S15\FOpt\RB3LYP\6-31G(d)\C2H6O1\MAREKDOS\05-Sep-2009\0\\# OPT
B3LYP/6-31G(D)\\komentarz moj np.: moje pierwsza optymalizacja, wole
liczyc niz konsumowac ;)\\0,1\O,0.2191201171,-0.3795271757,-1.13665085
17\C,0.2373651866,-0.4111285632,0.2877256129\H,1.2715278049,-0.4276420
155,0.6701513276\C,-0.4795543291,0.8306124631,0.7901795859\H,-0.265415
0533,-1.3149963884,0.6701513276\H,0.6764276346,-1.1716070307,-1.457449
2431\H,-0.4924568475,0.8529602804,1.8857642295\H,0.0231001682,1.733736
3286,0.4285789895\H,-1.5130097881,0.8468628318,0.4285789895\\Version=I
A64L-G03RevC.02\State=1-A'\HF=-155.0338051\RMSD=6.252e-09\RMSF=9.767e-
05\Dipole=0.2646027,-0.4583053,0.3129192\PG=CS [SG(C2H2O1),X(H4)]\\@
Sacred cows make the best hamburger.
-- Mark Twain
Job cpu time: 0 days 0 hours 0 minutes 45.9 seconds.
File lengths (MBytes): RWF= 18 Int= 0 D2E= 0 Chk= 10 Scr= 1
Normal termination of Gaussian 03 at Sat Sep 5 11:10:34 2009.
Entalpia, entalpia swobodna (energia Gibbsa), entropia, energia elektronowa
Wyżej wymienione energie można obliczyć poprzez wyliczenie drgań cząsteczki. W tym celu w pliku tekstowym z obliczeniami w sekcji # umieszczamy napis freq. Obliczenia te można także połaczyć z wcześniejszą optymalizacja cząsteczki w tym celu dodajemy opcję opt np.:
# Opt Freq B3LYP/6-31G(d)
Obliczenia zostaną przeprowadzone w warunkach normalnych (temp. 298.15K, ciśnienie 1bar), oczywiście można zadać także inne warunki, o czym będzie później. W otrzymanym pliku z wynikami na początku sprawdzamy, czy cząsteczka znajduje się w minimum energetycznym. W tym celu szukamy pierwszej linii z napisem: Frequencies --
Low frequencies --- -9.7622 0.0006 0.0006 0.0006 9.5432 21.8606
Low frequencies --- 253.0301 307.2737 417.1536
Diagonal vibrational polarizability:
3.9168100 2.0707571 41.3969289
Harmonic frequencies (cm**-1), IR intensities (KM/Mole), Raman scattering
activities (A**4/AMU), depolarization ratios for plane and unpolarized
incident light, reduced masses (AMU), force constants (mDyne/A),
and normal coordinates:
1 2 3
A A A
Frequencies -- 253.0301 307.2615 417.1529
Red. masses -- 1.1497 1.0711 2.6404
Frc consts -- 0.0434 0.0596 0.2707
IR Inten -- 26.0210 105.2035 11.2676
Atom AN X Y Z X Y Z X Y Z
1 8 -0.00 0.00 0.08 -0.00 0.00 0.04 0.22 -0.06 0.00
2 6 0.00 -0.00 -0.05 -0.00 -0.00 0.06 -0.02 0.19 -0.00
3 1 0.03 0.08 -0.11 0.01 -0.02 0.08 -0.07 0.17 0.01
4 6 -0.00 0.00 -0.01 0.00 0.00 -0.03 -0.19 -0.05 -0.00
5 1 -0.03 -0.08 -0.11 -0.01 0.02 0.08 -0.07 0.17 -0.01
6 1 -0.00 0.00 -0.35 -0.00 0.00 -0.92 0.05 -0.35 -0.00
7 1 0.00 0.00 0.54 -0.00 -0.00 -0.25 0.09 -0.41 0.00
8 1 -0.29 -0.35 -0.25 0.15 0.11 0.04 -0.49 -0.10 -0.01
9 1 0.29 0.35 -0.25 -0.15 -0.11 0.04 -0.49 -0.10 0.01
4 5 6
A A A
Frequencies -- 830.4216 911.8792 1043.4321
Red. masses -- 1.0798 2.2036 2.1690
Frc consts -- 0.4387 1.0796 1.3914
IR Inten -- 0.0130 10.1043 55.6967
Atom AN X Y Z X Y Z X Y Z
1 8 0.00 0.00 0.02 0.14 0.09 -0.00 0.07 0.09 -0.00
Jeżeli zaznaczona cyfra nie jest ujemna, znaczy to, że struktura jest minimum energetycznym. Jeżeli wartości liczby jest ujemna, wówczas mamy do czynienia ze stanem przejściowym (maksimum energetyczne). Jeżeli poszukujemy struktury w minimum, musimy dokonać niewielkiej zmiany geometrii cząsteczki i potworzyć optymalizację oraz obliczenia częstości. Jeżeli wszystko jest w porządku szukamy w pliku tekstowym linii z napisem Zero-point correction= znajdujemy to w odpowiednim fragmencie pliku.
Vibrational temperatures: 364.05 442.08 600.19 1194.79 1311.99
(Kelvin) 1501.26 1617.76 1718.95 1859.58 1885.61
2054.69 2132.38 2178.18 2204.01 2246.67
4291.09 4324.44 4399.05 4500.17 4508.95
5394.86
Zero-point correction= 0.080327 (Hartree/Particle)
Thermal correction to Energy= 0.084586
Thermal correction to Enthalpy= 0.085530
Thermal correction to Gibbs Free Energy= 0.054970
Sum of electronic and zero-point Energies= -154.953473
Sum of electronic and thermal Energies= -154.949215
Sum of electronic and thermal Enthalpies= -154.948271
Sum of electronic and thermal Free Energies= -154.978831
E (Thermal) CV S
KCal/Mol Cal/Mol-Kelvin Cal/Mol-Kelvin
Total 53.078 13.368 64.320
Electronic 0.000 0.000 0.000
Translational 0.889 2.981 37.406
Rotational 0.889 2.981 22.318
Vibrational 51.301 7.407 4.596
Vibration 1 0.664 1.758 1.709
Vibration 2 0.697 1.660 1.377
Vibration 3 0.780 1.433 0.902
W tym fragmencie znajdują się najbardziej interesujące wartości. Dziś poznamy dwie wartości entalpię ΔH, oraz entalpię swobodną ΔG (energia Gibbsa), które znajdują się w określonych liniach tekstu:
ΔH = Sum of electronic and thermal Enthalpies = –154.948271
ΔG = Sum of electronic and thermal Free Energies = –154.978831
Mając reakcję np.: CH3OH + 3/2O2 → CO2 + 2H2O, możemy obliczyć ΔH, ΔG poszczególnych substratów oraz produktów, oraz korzystając z prawa Hessa wyliczyć energie spalania metanolu.
Ćwiczenie
Oblicz entalpię spalania metanolu (w kcal/mol) oraz podaj w komentarzach twoją wyliczoną wartość energii dla poszczególnych cząsteczek oraz dla całej reakcji. Dodatkowo podaj czas obliczeń. Wybierz jedną z podanych poniżej metod. Swoje wyniki umieść w komentarzu. Po wykonaniu ćwiczenia przez kilka osób można porównać dokładność metod oraz czas jaki jest potrzebny z przeprowadzenia obliczeń.
AM1
B3LYP/6-31G(d)
B3LYP/6-31G(d,p)
B3LYP/6-31++G(d,p)
MP2/6-31G(d)
MP2/6-31G(d,p)
*Pamiętaj aby w obliczeniach przeprowadzić także optymalizację najplepeij wykonaj ją razem z obliczeniami częstości #opt freq METODA/BAZA.
*aby przeliczyć energię z Hartree na kcal/mol wystarczy pomnożyc wartość energii przez 627.509.
[1] Jeżeli mamy w pliku kilka linii zawierające słowa SCF Done: wybieramy ostatnią linię zawierającą te słowa.
[2] Aby przeliczyć na kcal/mol mnożymy wartość uzyskaną w Hartree przez 627.509.
[3] D oznacza to samo co E, różnica polega na tym ze liczba przed została obliczona z podwójną precyzja (Double precision) czyli w naszym wypadku mamy -0.98712003835432 * 103 Hartree
Funkcje termodynamiczne ΔE, ΔH, ΔG
Obliczenie energii układu jest jednym z podstawowych obliczeń, energię wykorzystuje się do:
- porównywania konformerów, wybieranie konformeru o najmniejszej energii
- zestawiania diagramów energetycznych dla ścieżek reakcji chemicznych, dzięki temu można wybrać odpowiedni mechanizm
- obliczenia entalpii reakcji chemicznych
- spektroskopii
- obliczania energii oddziaływań, trwałości układów – korelacja z wynikami eksperymentalnymi
- Czytaj więcej: Funkcje termodynamiczne ΔE, ΔH, ΔG
Przeprowadzenie obliczeń optymalizacji
Po zbudowaniu cząsteczki oraz jej zapisu do pliku tekstowego, możemy zadać odpowiednie parametry w sekcji oznaczonej znakiem #. Dla osób nie pamiętających struktury pliku wsadowego, proszę zapoznać się z artykułem struktura pliku wsadowego Gaussiana. Aby wykonać optymalizację wpisujemy słowo kluczowe opt a następnie metodę i bazę np.:
# opt B3LYP/6-31G(d)
Korzystając z programu z interfejsem graficznym – Molden, możemy te opcje ustawić bezpośrednio z okna programu po kliknięciu przycisku Z-Mat Editor, a następnie zaznaczeniu na dole tego okna przycisku Gaussian, oraz wybranie Submit Job (ilustracja tego została przedstawiona w cześniejszych materiałach). Tak przygotowany plik, po ustawieniu odpowiednich parametrów dla pamięci (%mem=1500MB), procesora (%nproc=2) oraz check pointa (%chk=nazwa.chk) możemy skopiować na komputer ICM oraz wykonać obliczenia.
Ćwiczenie
Uprzejmie proszę o wykonanie tego ćwiczenia, ponieważ jest to etap obliczeń, którego nie można pominąć!
Proszę zbudować cząsteczkę etanolu, następnie zapisać ją w układzie kartezjanskim pod nazwą etanol_xyz.inp, oraz w postaci Z-macierzy pod nazwą etanol_zmat.inp[1] Następnie proszę przed wysłaniem obliczeń na komputery ICM otworzyć te pliki w edytorze tekstowym. Sugeruję na samym początku uruchomić program WinScp i zalogować się na komputer w Centrum obliczeniowym ICM, a następnie wybrać w oknie po prawej stronie katalog, w którym zapisany jest plik i otworzyć go do edycji (F4 Edit).
Proszę o dopisanie odpowiednich parametrów obliczeń takich jak chk, mem, nproc, #, komentarz, ładunek i multipletowość, oraz skasowanie jednej linii z formatu Kartezjańskiego mówiącej o liczbie atomów w cząsteczce, dodatkowo proszę na końcu pliku wprowadzić kilka pustych linii. Obliczenia proszę ustawić tak, aby przeprowadzić optymalizacje metodą B3LYP/6-31G(d)[2].
Pliki powinny mniej więcej wyglądać następująco:
etanol_xyz.inp:
%mem=1500mb
%nproc=2 %chk=etanol_xyz.chk
# opt B3LYP/6-31G(d)
komentarz moj np.: moje pierwsza optymalizacja, wole liczyć niz konsumować ;)
0 1
O 0.000000 0.000000 0.000000
C 0.000000 0.000000 1.400000
H 1.026720 0.000000 1.762996
C -0.683537 1.183920 1.883333
H -0.513360 -0.889165 1.763000
H 0.447834 -0.775672 -0.316667
H -0.683537 1.183920 2.972333
H -0.170177 2.073085 1.520333
H -1.710256 1.183920 1.520333
etanol_zmat.inp:
%mem=1500mb
%nproc=2
%chk=etanol_zmat.chk
# opt B3LYP/6-31G(d)
komentarz moj np.: moja druga optymalizacja
0 1
o
c 1 co2
h 2 hc3 1 hco3
c 2 cc4 1 cco4 3 dih4
h 2 hc5 1 hco5 3 dih5
h 1 ho6 2 hoc6 4 dih6
h 4 hc7 2 hcc7 1 dih7
h 4 hc8 2 hcc8 7 dih8
h 4 hc9 2 hcc9 7 dih9
co2 1.400000
hc3 1.089000
hco3 109.471
cc4 1.450000
cco4 109.471
dih4 120.000
hc5 1.089000
hco5 109.471
dih5 -120.000
ho6 0.950000
hoc6 109.471
dih6 180.000
hc7 1.089000
hcc7 109.471
dih7 180.000
hc8 1.089000
hcc8 109.471
dih8 120.000
hc9 1.089000
hcc9 109.471
dih9 240.000
[1] W programie Molden zapis dokonujemy poprzez wybranie opcji Write po skonstruowaniu cząsteczki. Przeważnie cząsteczka jest zapisywana w tym samym katalogu gdzie program został zainstalowany, chyba że została podana inna opcja.
[2] Nacisk na tekstowy sposób edycji plików wsadowych nie jest przypadkowy. Doświadczenie pokazało niejednokrotnie, że ustawienie w pliku tekstowym parametrów jest bardziej wydajne niż korzystanie z różnych opcji programów, w szczególności, że nie wszystkie metody i bazy są uwzględnione w opcjach programów graficznych.
Optymalizacja geometrii cząsteczki
Optymalizacja geometrii cząsteczki jest jednym z podstawowych obliczeń w modelowaniu molekularnym i jest wykorzystywana do:
– przeprowadzenia następnych obliczeń na bazie uzyskanej geometrii cząsteczki (energia, drgania, rozkład ładunku)
– porównywania uzyskanych parametrów geometrycznych (odległości atomów, długości wiązań, kąty) z wynikami eksperymentalnymi np: X-ray, NMR, IR itp.
– ilustracji konformerów, sposobu ułożenia grup w cząsteczce, sposobu oddziaływania
– prezentacji i ilustracja wyników badań
Znalezienie najlepszej geometrii (optymalnej) cząsteczki jest rzeczą bardzo ważną, ponieważ w dużej mierze to od niej zależą następne obliczenia jakie wykonujemy dla układu – np: wartość energii, drgania, rozkład ładunku.
Strona 6 z 23