Obliczanie energii oddziaływania

Pierwszym etapem obliczeń energii oraz sposobów oddziaływania jest optymalizacja kompleksu cząsteczek. Należy tak zbudować cząsteczki i tak ustawić je koło siebie, żeby jak najlepiej odzwierciedlały możliwe oddziaływania. Następnie taki układ poddajemy procesowi optymalizacji. Czasami do końca nie wiemy, w którym miejscu następuje oddziaływanie wówczas próbujemy różnych sposobów.

Jeżeli otrzymaliśmy zoptymalizowany kompleks dwóch cząsteczek, wówczas przystępujemy do analizy wiązań wodorowych (odległości, kąty) lub do klasyfikacji występujących oddziaływań. Tą analizę opieramy na strukturze kompleksu. Jednym z ważniejszych parametrów są obliczenia energii. Załóżmy, że nasz kompleks składa się z dwóch fragmentów: A i B. Jeżeli znamy energie zoptymalizowanych osobno fragmentów A i B oraz energię kompleksu, wówczas możemy obliczyć energię dysocjacji kompleksu. Liczona jest ona jako różnica energii kompleksu AB i osobno fragmentów A i B. Można to wyrazić wzorem:

ΔEdysocjacji = EAB – EA’ – EB'

 

Kolejną energię, którą możemy obliczyć to energia oddziaływania. Energia oddziaływania jest różnicą energii kompleksu EAB i sumy energii poszczególnych fragmentów wchodzących w skład tego kompleksu.

Grafika_oddziaływań

 

ΔEoddzial = EAB – (EA + EB)

Energie poszczególnych fragmetów A i Bliczone dla geometrii takiej samej jaką posiadają te fragmenty w kompleksie. Należy wspomnieć, że energia kompleksu jest liczona z zastosowaniem bazy AB, zaś energia poszczególnych fragmentów – w bazach tych fragmentów (EA w bazie A i EB w bazie B, rysunek powyżej). W takim podejściu obliczona energia zawiera błąd związany z rozmiarem bazy funkcyjnej dla podukładów A i B. Rozszerzenie bazy dla podukładów A oraz B o bazę funkcyjną całego kompleksu AB wiąże się z obniżeniem energii tych fragmentów i lepsze zbliżenie do wartości oczekiwanych. Błąd ten nazywa się błędem superpozycji bazy BSSE. Błąd superpozycji bazy, zgodnie ze schamatem zaproponowanym przez Boysa i Bernardiego, jest to różnica sumy energii substratów liczonych w pełnej bazie AB i sumy energii substratów liczonych osobno: EA w bazie A i EB w bazie B.

BSSE = (EAB + EAB) – (EA + EB)

Chcąc obliczyć energię oddziaływania kompleksu z uwzględnieniem BSSE, od energii kompleksu AB liczonego w bazie AB trzeba odjąć sumę energii podukładu A liczonego w bazie AB i podukładu B liczonego w bazie AB rysunek poniżej:

 

Grafika_oddziaływań

 

ΔEoddzialBSSE = EAB – (EA(AB) + EB(AB))

Błąd superpozycji bazy można policzyć osobno, korzystając z wcześniej zoptymalizowanej struktury kompleksu, lub uwzględnić BSSE podczas procesu optymalizacji kompleksu. Dla porównania przedstawione są energie oddziaływania etylenu z wodorem liczone w różny sposób metodą MP2/aug-cc-pVTZ. Energia oddziaływania bez BSSE wynosi 0.48 kcal/mol, energia z uwzględnieniem BSSE -0.35, zaś energia z uwzględnieniem BSSE podczas procesu optymalizacji -0.36. Widzimy że ostanie przypadki są najbardziej zbliżone do siebie i odpowiadają wartości eksperymentalnej.

Ukazany powyżej sposób obliczeń nie jest trudny do przeprowadzenia. Do obliczeń wybierzmy kompleks chlorowodoru z cząsteczką wody (H2O:HCl). W pierwszym etapie optymalizujemy kompleks oraz obliczamy częstości, żeby sprawdzić czy kompleks jest w minimum energetycznym. Plik obliczeń (h2ohcl.inp) wygląda następująco:

%chk=h2ohcl.chk
%nproc=2
%mem=1500mb
# opt freq b3lyp/6-31G(d)

kompleks H2O z HCl

0 1
O 0.448141 0.316244 0.198285
H 0.480086 -0.077493 1.083718
H 1.167086 -0.111305 -0.292018
H -1.108258 -0.165684 -0.494586
Cl -2.260912 -0.606682 -0.977408

Następnie sprawdzamy częstości, jeżeli wszystkie są  dodatnie, to przechodzimy dalej do obliczeń. Za pomocą programu Molden wczytujemy plik z wynikami h2ohcl.log. Przechodzimy do ostatniego kroku optymalizacji.

 

Grafika_oddziaływań

Następnie zapisujemy strukturę w postaci kartezjańskiej. Zapisany plik h2ohcl_BSSE.inp wygląda następująco:

5
scf done: -537.218903
O 0.448728 0.324467 0.209110
H 0.497205 -0.076345 1.091713
H 1.136485 -0.127585 -0.305288
H -1.114323 -0.146537 -0.489077
Cl -2.241953 -0.618921 -0.988466

Następnie dokonujemy modyfikacji pliku. W pierwszym etapie nagłówek. W linii ze znakiem # wstawiamy opcje counterpoise=2. Ładunek i multipletowość – tutaj  nie mamy już dwóch cyfr, (0 1) tylko musimy podać trzy pary cyfr. Pierwsza para informuje nas o ładunku i multipletowości kompleksu, druga para o ładunku i multipletowości pierwszego fragmentu, druga para informuje nas o drugim fragmencie. W układzie kartezjańskim dodajemy kolejną kolumnę informującą jakie atomy wchodzą w skład poszczególnych fragmentów – 1 oznacza pierwszy fragment, 2 oznacza drugi fragment (numerację atomów można podejrzeć graficznie programem Molden). Końcowy rezultat modyfikacji pliku powinien wyglądać następująco:

%chk=h2ohcl_BSSE.chk
%nproc=2
%mem=1500mb
# counterpoise=2 b3lyp/6-31G(d)

kompleks H2o z HCl obliczenie bledu superpozycji bazy

0 1, 0 1, 0 1
O 0.448728 0.324467 0.209110 1
H 0.497205 -0.076345 1.091713 1
H 1.136485 -0.127585 -0.305288 1
H -1.114323 -0.146537 -0.489077 2
Cl -2.241953 -0.618921 -0.988466 2

W opcjach można użyć także słowa kluczowego opt, ale wydłuża to obliczenia, a uzyskane rezultaty są zbliżone to uzyskanych powyższą metodą. Po skończeniu obliczeń otwieramy plik w formie tekstowej i szukamy następujących informacji: SCF Done.

Pierwsza linia SCF Done reprezentuje energię całej cząsteczki.

There are    40 symmetry adapted basis functions of A   symmetry.
Integral buffers will be 131072 words long.
Raffenetti 2 integral format.
Two-electron integral symmetry is turned on.
40 basis functions, 92 primitive gaussians, 40 cartesian basis functions
14 alpha electrons 14 beta electrons
nuclear repulsion energy 47.1715240149 Hartrees.
NAtoms= 5 NActive= 5 NUniq= 5 SFac= 1.00D+00 NAtFMM= 60 Big=F
One-electron integrals computed using PRISM.
NBasis= 40 RedAO= T NBF= 40
NBsUse= 40 1.00D-06 NBFU= 40
Harris functional with IExCor= 402 diagonalized for initial guess.
ExpMin= 1.43D-01 ExpMax= 2.52D+04 ExpMxC= 3.78D+03 IAcc=1 IRadAn= 1 AccDes= 1.00D-06
HarFok: IExCor= 402 AccDes= 1.00D-06 IRadAn= 1 IDoV=1
ScaDFX= 1.000000 1.000000 1.000000 1.000000
Initial guess orbital symmetries:
Occupied (A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A)
(A) (A)
Virtual (A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A)
(A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A)
(A) (A)
The electronic state of the initial guess is 1-A.
Warning! Cutoffs for single-point calculations used.
Requested convergence on RMS density matrix=1.00D-04 within 128 cycles.
Requested convergence on MAX density matrix=1.00D-02.
Requested convergence on energy=5.00D-05.
No special actions if energy rises.
Keep R1 integrals in memory in canonical form, NReq= 1255471.
SCF Done: E(RB+HF-LYP) = -537.218903306 A.U. after 6 cycles
Convg = 0.5621D-05 -V/T = 2.0038
S**2 = 0.0000

 

Druga linia SCF Done reprezentuje fragment pierwszy w bazie dimer, należy zwrócić uwagę na fragment poprzedzający SCF Done, kilka lini wyżej jest użyty skrót DCBS – oznacza to, że obliczenia wykonywane są w bazie dimeru.

 Counterpoise: doing DCBS calculation for fragment   1 NewBq=T
Basis read from rwf: (6D, 7F)
There are 40 symmetry adapted basis functions of A symmetry.
Integral buffers will be 131072 words long.
Raffenetti 2 integral format.
Two-electron integral symmetry is turned on.
40 basis functions, 92 primitive gaussians, 40 cartesian basis functions
5 alpha electrons 5 beta electrons
nuclear repulsion energy 9.0679564666 Hartrees.
NAtoms= 5 NActive= 5 NUniq= 5 SFac= 1.00D+00 NAtFMM= 60 Big=F
One-electron integrals computed using PRISM.
NBasis= 40 RedAO= T NBF= 40
NBsUse= 40 1.00D-06 NBFU= 40
Harris functional with IExCor= 402 diagonalized for initial guess.
ExpMin= 1.43D-01 ExpMax= 2.52D+04 ExpMxC= 3.78D+03 IAcc=1 IRadAn=
1 AccDes= 1.00D-06
HarFok: IExCor= 402 AccDes= 1.00D-06 IRadAn= 1 IDoV=1
ScaDFX= 1.000000 1.000000 1.000000 1.000000
Initial guess orbital symmetries:
Occupied (A) (A) (A) (A) (A)
Virtual (A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A)
(A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A)
(A) (A) (A) (A) (A) (A) (A) (A) (A) (A) (A)
The electronic state of the initial guess is 1-A.
Warning! Cutoffs for single-point calculations used.
Requested convergence on RMS density matrix=1.00D-04 within 128 cycles.
Requested convergence on MAX density matrix=1.00D-02.
Requested convergence on energy=5.00D-05.
No special actions if energy rises.
Keep R1 integrals in memory in canonical form, NReq= 1255471.
SCF Done: E(RB+HF-LYP) = -76.4113995017 A.U. after 5 cycles
Convg = 0.4853D-05 -V/T = 2.0083
S**2 = 0.0000

 

Trzecia linia SCF Done odpowiada energii drugiego fragmentu w bazie dimeru.

Counterpoise: doing DCBS calculation for fragment 2 NewBq=T
.
.
SCF Done: E(RB+HF-LYP) = -460.795541716 A.U. after 4 cycles
Convg = 0.1717D-04 -V/T = 2.0031
S**2 = 0.0000

Czwarta linia SCF Done informuje o energii fragmentu pierwszego obliczony w swojej bazie (monomeru). Użyty jest skrót MCBS.

Counterpoise: doing MCBS calculation for fragment 1
.
.
Keep R1 integrals in memory in canonical form, NReq= 885565.
SCF Done: E(RB+HF-LYP) = -76.4089088740 A.U. after 5 cycles

Piąta linia SCF Done informuje nas o energii drugiego fragmentu w swojej bazie.

Counterpoise: doing MCBS calculation for fragment 2
.
.
SCF Done: E(RB+HF-LYP) = -460.795172630 A.U. after 4 cycles

Pod koniec możemy znaleźć energie układu z korektą BSSE, oraz błąd superpozycji bazy. Wszystkie energie podane są w jednostkach Hartree.

Counterpoise: corrected energy = -537.216043592858
Counterpoise: BSSE energy = 0.002859713284

Mając te wszystkie dane można policzyć energię oddziaływania z uwzględnieniem superpozycji bazy:

Eoddziaływania = EAB – EA(DCBS) – EB(DCBS)

Możemy także obliczyć sam bład superpozycji bazy:

BSSE = EA(DCBS) + EB(DCBS) – (EA(MCBS) + EB(MCBS))

Energia wyliczona powinna się zgadzać z wyliczoną wartością corrected energy.

Obliczona energia oddziaływania dla naszego przykładu H2O:HCl wynosi -0.01196 Hartree, czyli 7.50 kcal/mol, błąd superpozycji bazy wynosi 1.79 kcal/mol.

Przeprowadzając obliczenia dla kilkunastu przypadków, osoby prowadzące obliczenia ułatwiają sobie analizę danych poprzez pisanie własnych skryptów (malych programów). Postaramy się w przyszłości o przekazanie kompletu skryptów, co powinno znacznie ułatwić obliczenia.

Ćwiczenie

Oblicz energię oddziaływania, z uwzględnieniem BSSE, dla dimeru wody. Podaj w komentarzach twoją wyliczoną wartość w kcal/mol. Wybierz jedną z podanych poniżej metod. Swoje wyniki umieść w komentarzu. Dodatkowo możesz napisać czas jaki był potrzeby do obliczeń BSSE. 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)

A4 Infociacho