Wykonanie powierzchni energii potencjalnej dla układu bez relaksacji możemy zrobić przygotowując serie plików wsadowych, w którym każdy z nich jest jednym punktem na wykresie. Proces ten możemy z automatyzować i przygotować jeden plik wsadowy zawierający instrukcje postępowania. Jedna z metod zostanie przedstawiona na przykładzie cząsteczki H2O2[1].

W tym celu budujemy strukturę startową i zapisujemy ją w postaci Z-macierzy (opis Z-macierzy znajduje się na portalu molnet.eu):

H
O 1 B1
O 2 B2 1 A1
H 3 B3 2 A2 1 D1

B1 0.96000000
B2 1.32000000
B3 0.96000000
A1 109.50000032
A2 109.49999973
D1 0.00000000

 

(korzystając z programu molden uzyskujemy następująca formę Z-macierzy. Widzimy że różnica jest w nazywaniu zmienych-parametrów)

h
o 1 oh2
o 2 oo3 1 ooh3
h 3 ho4 2 hoo4 1 dih4

oh2 0.947000
oo3 1.480000
ooh3 109.471
ho4 0.947000
hoo4 109.471
dih4 180.000

Tak przygotowaną strukturę używamy do budowy pliku startowego. W linii definiującej polecenia # wpisujemy słowo kluczowe scan, oraz metodę i bazę lub samą metodę. Dodatkowo możemy dodać słowo kluczowe nosymm co pozwoli uniknąć pewnych błędów które kończą obliczenia. Parametr, który zmieniamy w tym przypadku – kąt dwuścienny powinien być zdefiniowany w Z-macierzy jako pojedyncza zmienna, tak jak jest to w przykładzie: D1. Najlepiej jest, aby zmiana tego parametru nie pociągała za sobą zmian innych parametrów, które ustawione są na „sztywno”. W tym celu dla bardziej skomplikowanych cząsteczek można przeprowadzić „przenumerowanie” atomów w Z-macierzy przykładowo za pomocą programu Molden.

Położenie drugiego atomu H zależy od odległości od drugiego atomu tlenu (H–O) oraz od kąta H–O–O i badanego kąta dwuściennego. Zmiana kąta dwuściennego nie powoduje zmiany poprzednich parametrów. Gdyby położenie atomu wodoru było definiowane odległością między atomami wodoru H∙∙∙∙∙∙∙H wówczas obliczenia zakończyłyby się błędem. Przy zmiennej D1 dopisujemy literę S (lub s, wielkość liter nie ma znaczenia), która oznacza, że dany parametr będzie zmieniany, następnie cyfrę, która informuje ile razy będzie parametr zmieniany (liczba kroków) oraz następną wartość definiująca jak wielka będzie zmiana (wielkość kroku). Przy tej cyfrze używamy znaków + i –, dzięki którym możemy określić, w którym kierunku zachodzi ta zmiana, nie użycie żadnego znaku powoduje użycie domyślnej wartości +. Wielkość kroku musi (!) być definiowana jako cyfra z kropką, liczby całkowite zapisujemy jako np.: 20.0. Użycie samej cyfry 20 spowoduję błąd. Uzyskujemy plik podobny do poniższego:

# scan am1 nosymm

Przyklad skanu bez relaksacji.

0 1
H
O 1 B1
O 2 B2 1 A1
H 3 B3 2 A2 1 D1

B1 0.96000000
B2 1.32000000
B3 0.96000000
A1 109.50000032
A2 109.49999973
D1 0.00000000 S 18 +20.0

Kąt dwuścienny zdefiniowany parametrem D1 będzie się zmieniał 18 razy po 20°. Wyniku tych obliczeń otrzymujemy plik z danymi, w którym znajduje się odpowiedni fragment na przy końcu pliku:

Summary of the potential surface scan:
N D1 SCF
---- --------- -----------
1 0.0000 -0.04440
2 20.0000 -0.04540
3 40.0000 -0.04785
4 60.0000 -0.05059
5 80.0000 -0.05267
6 100.0000 -0.05371
7 120.0000 -0.05380
8 140.0000 -0.05333
9 160.0000 -0.05277
10 180.0000 -0.05254
11 200.0000 -0.05277
12 220.0000 -0.05333
13 240.0000 -0.05380
14 260.0000 -0.05371
15 280.0000 -0.05267
16 300.0000 -0.05059
17 320.0000 -0.04785
18 340.0000 -0.04540
19 360.0000 -0.04440
---- --------- -----------

Wartości kąta znajdują się pod opisem D1, energia zaś w Hartree pod napisem SCF. Korzystając z tych wyników możemy sporządzić wykres:

Wykres_H2O2

Jak widzimy na wykresie otrzymane wartości kątów dla tej cząsteczki można powiązać ze sobą (τ’HOOH = 360°-τHOOH), umożliwia to skrócenie obliczeń o połowę obliczonych punktów. Ilość obliczonych punktów jest większa od zdefiniowanej liczby o jeden, ponieważ została uwzględniona pierwsza struktura. Z właściwości symetrii cząsteczek należy korzystać ,gdyż pozwalają one znacznie skrócić czas obliczeń. Na wykresie tym widoczne jest minimum energetyczne, w której cząsteczka jest w stanie równowagowym, oraz dwa maksima energetyczne reprezentujące stan przejściowy obrotu grupy –O–H.

Możemy także przeprowadzić obliczenia z optymalizacją wszystkich innych parametrów z wyłączeniem zmieniającego się. W tym celu w podanym przykładzie w linijce z # zamiast scan umieszczamy opt=z-matrix. Wynik naszych obliczeń znajdujemy na końcu pliku jako wypisanie zmiennych.

 

                     ----------------------------
! Optimized Parameters !
! (Angstroms and Degrees) !
---------------------- ----------------------
! Name Value Derivative information (Atomic Units) !
------------------------------------------------------------------------
! B1 0.978 -DE/DX = 0.0 !
! B2 1.3026 -DE/DX = 0.0 !
! B3 0.978 -DE/DX = 0.0 !
! A1 110.5471 -DE/DX = 0.0 !
! A2 110.5477 -DE/DX = 0.0 !
! D1 360.0 -DE/DX = 0.0 !
------------------------------------------------------------------------
GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad

Summary of Optimized Potential Surface Scan
1 2 3 4 5
EIGENVALUES -- -0.04517 -0.04626 -0.04898 -0.05213 -0.05459
B1 0.97814 0.97889 0.98107 0.98339 0.98473
B2 1.30265 1.30047 1.29566 1.29105 1.28947
B3 0.97814 0.97889 0.98107 0.98339 0.98473
A1 110.58674 110.37373 109.88186 109.20451 108.40023
A2 110.58673 110.37373 109.88185 109.20451 108.40021
D1 0.00000 20.00000 40.00000 60.00000 80.00000
.
.
.
.
.
16 17 18 19
EIGENVALUES -- -0.05213 -0.04898 -0.04626 -0.04518
B1 0.98366 0.98124 0.97887 0.97798
B2 1.29121 1.29594 1.30072 1.30261
B3 0.98363 0.98124 0.97887 0.97799
A1 109.17725 109.84687 110.35104 110.54706
A2 109.18035 109.85182 110.35419 110.54774
D1 300.00000 320.00000 340.00000 360.00000

Gdzie pod Name są nazwy definiowanych parametrów w macierzy-Z, cyfry od 1 do 19 to poszczególne kroki. Obliczona energia dla poszczególnych kroków podawana jest po słowie EIGENVALUES w jednostkach atomowych (Hartree).

Obliczenia dotyczące zmiany parametru geometrycznego z relaksacją możemy także wykonać stosując inną procedurę. Obliczenia te wykonujemy ze słowem kluczowym opt=ModRedundant lub z jego synonimem: opt=AddRedundand. Pod zdefiniowaniu Z-macierzy, lub układu kartezjańskiego pod nią, definiujemy zmieniający się parametr nie odwołując się do zmiennych Z-macierzy, dlatego można przeprowadzić to także dla układu kartezjańskiego. Definicję rozpoczynamy od określenia zmieniającego się parametru używając skrótów: B – dla wiązania (ang. bond), A – dla kąta (ang. angle), D – dla kąta torsyjnego (ang. dihedral angle), następnie podajemy, jakie atomy są opisywane parametrem. Po tym wpisujemy litere S, a następnie liczbę i rozmiar kroków podobnie jak wyżej. Przykłady podane są poniżej:

  # opt=modredundant am1 nosymm

Skan z relaksacja

0 1
H -1.68797589 0.12908249 0.90493540
O -1.36752130 0.12820513 0.00000000
O -0.04752130 0.12820513 0.00000000
H 0.27293328 0.12908249 0.90493541

D 1 2 3 4 S 9 20.0

w innym zapisie:

# opt=modredundant am1 nosymm

Skan z relaksacja

0 1
H
O 1 B1
O 2 B2 1 A1
H 3 B3 2 A2 1 D1

B1 0.96000000
B2 1.32000000
B3 0.96000000
A1 109.50000032
A2 109.49999973
D1 0.00000000

D 1 2 3 4 S 9 20.0

W ostatniej linii zapisujemy informację, że ma być wykonany skan (S) kąta dwuściennego (D) zbudowanego z atomów 1, 2, 3, 4 w 9 krokach zwiększając o 20° kąt. Wynik obliczeń znajduje się przy końcu pliku wynikowego:

                             ----------------------------
! Optimized Parameters !
! (Angstroms and Degrees) !
-------------------------- --------------------------
! Name Definition Value Derivative Info. !
--------------------------------------------------------------------------------
! R1 R(1,2) 0.9811 -DE/DX = 0.0 !
! R2 R(2,3) 1.3093 -DE/DX = 0.0 !
! R3 R(3,4) 0.9811 -DE/DX = 0.0 !
! A1 A(1,2,3) 104.4447 -DE/DX = 0.0 !
! A2 A(2,3,4) 104.4447 -DE/DX = 0.0 !
! D1 D(1,2,3,4) 180.0 -DE/DX = 0.0 !
--------------------------------------------------------------------------------
GradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGradGrad

Summary of Optimized Potential Surface Scan
1 2 3 4 5
EIGENVALUES -- -0.04518 -0.04626 -0.04898 -0.05213 -0.05459
R1 0.97795 0.97896 0.98114 0.98345 0.98475
R2 1.30256 1.30052 1.29579 1.29108 1.28942
R3 0.97795 0.97896 0.98114 0.98345 0.98475
A1 110.55266 110.35803 109.86820 109.20113 108.40458
A2 110.55266 110.35803 109.86820 109.20113 108.40458
D1 0.00000 20.00000 40.00000 60.00000 80.00000
6 7 8 9 10
EIGENVALUES -- -0.05591 -0.05632 -0.05632 -0.05623 -0.05620
R1 0.98465 0.98351 0.98220 0.98134 0.98109
R2 1.29201 1.29776 1.30385 1.30799 1.30935
R3 0.98465 0.98351 0.98220 0.98134 0.98109
A1 107.42643 106.32060 105.32093 104.65976 104.44470
A2 107.42643 106.32060 105.32093 104.65976 104.44470
D1 100.00000 120.00000 140.00000 160.00000 180.00000

Zdefiniowane parametry geometryczne w Z-macierzy zostały zastąpione przez nowe zmienne zdefiniowane w sekcji Optimized Parameters. Kąt, który zmienialiśmy znajduję się na końcu tych parametrów. Podobnie jak wcześniejszych obliczeniach podana jest wartość energii w kolumnach dla każdego kroku oraz wartość zmieniającego się parametru. Opcja opt=ModRedundant posiada więcej użytecznych opcji – informację o niej można znaleźć w manualu Gaussiana. Jednymi z takich opcji jest zamrożenie odpowiedniego kata, odległości między atomami, wykonujemy to poprzez napisanie po podaniu atomów litery F (ang. freeze – zamrożenie).

W jednym pliku możemy dokonywać zmiany kilku parametrów, tak jak zostało to przedstawione dla przykładu PES sporządzonego dla cząsteczki H-fosfonianu metylowego. Należy jednak pamiętać, że liczba punktów, a co za tym idzie czas obliczeń rośnie. Dla zmiany dwóch kątów dwuściennych w 36 krokach po 10° mamy już aż 1296 punktów (36 ∙ 36 = 1296).

 

[1] Nadtlenek wodoru (H2O2) używany jest popularnie jako roztwór 3% pod nazwą wody utlenionej jako środek odkażający. Roztwory 80% używane są jako utleniacze w paliwach rakietowych. Inną cząsteczka z tej rodziny jest H2O3, która powstaje w reakcji fotolizy H2O2 z O3, jest ona jeszcze silniejszym utleniaczem niż nadtlenek wodoru. Engdahl, A. Nelander, B. Science 295, 482–483, 2002.

Dalej