1
Contents
1. Kivonat
3
2. Bevezetés
5
3. Káoszelmélet [1, 2]
6
4. A Bloch-egyenlet iteratív megoldása
10
4.1. Az iterációs séma
10
4.2. Ljapunov-exponens számítás
12
4.3. Példák
14
4.3.1. A számítás kiindulási paraméterei
14
4.3.2. Hidrogén molekula
14
4.3.3. Hélium atom
19
5. A Hartree-Fock s¶r¶ségmátrix iterációja
23
5.1. Az iterációs séma
23
5.2. Ljapunov-exponens számítás
26
5.3. Példák
29
5.3.1. Számolás ab initio Hartree-Fock szinten
29
5.3.2. A Hückel-szint¶ számítás kiindulási paraméterei
30
5.3.3. Butadién
31
5.3.4. Benzol és más molekulák
37
5.3.5. Számítás lineáris poliénre
39
5.3.6. Számítás nanocs®re
40
6. Összefoglalás
42
Hivatkozások
43
2 .
3
1. KIVONAT A dolgozatban két iteratív eljárást tanulmányozunk: az N -elektron Schrödingeregyenlettel ekvivalens Bloch-egyenlet megoldását, és egy Hartree-Fock szinten alkalmazható algoritmust a s¶r¶ségmátrix közvetlen meghatározására. A két eljárásban közös, hogy olcsón tudunk energiát számolni, mert egyikben sincs szükség mátrixok diagonalizálására. A dolgozat els® felében bevezetjük a hullámoperátort, majd ennek segítségével levezetjük a Bloch-egyenletet. A Bloch-egyenlet rekurzív alakját iterálva diagonalizálás nélkül számolhatunk energiát. Az iterációs paraméter (η ) változtatásával azt tapasztaljuk, hogy kaotikus eredmények léphetnek fel. Levezetünk egy összefüggést, melynek segítségével a Blochegyenlet Ljapunov-exponensét számolni lehet, és megvizsgáljuk a H2 molekula és a He atom esetén, hogy mennyire megbízhatóan jósolja az iteráció instabillá válását. A dolgozat második felében a s¶r¶ségmátrix (P) Bloch-típusú iterációjával foglalkozunk. A P-iteráció meg®rzi az idempotenciát és a spurt. A P-iterációt a dolgozatban f®ként Hückelközelítésben használjuk, de az algoritmus minden további nélkül általánosítható ab initio Hartree-Fock, vagy Kohn-Sham (DFT) szintre is [3]. Ennek az iterációnak is vannak kaotikus megoldásai. Legrészletesebben a butadién példáján tanulmányozzuk a kaotikus iterációk tulajdonságait: nemcsak az iterációs-, hanem a kon gurációs-, illetve a fázistérben is ábrázoljuk a kapott eredményeket. A Bloch-egyenlettel analóg levezetést végzünk a Ljapunovexponens számolásra, majd különböz® molekulákra teszteljük ennek predikciós képességét. Végzünk P-iterációs számolást egy 60 vízmolekulából álló klaszterre, egy lineáris fémes poliénre, valamint egy királis nanocs®re.
4 .
5
2. BEVEZETÉS ˆ Hamilton-operáA kvantumkémia egyik legfontosabb feladata, hogy egy molekula H torának sajátérték-problémája, a
ˆ i (r1 , r2 , . . . rN ) = Ei Ψi (r1 , r2 , . . . rN ) HΨ Schrödinger-egyenlet megoldásával megadja a molekula energiaszintjeit.
Ennek a
dierenciál-egyenletrendszernek megoldásával kapjuk a Ψi (r1 , r2 , . . . rN ) N -elektron hullámfüggvényeket, valamint az Ei energiákat. Egyelektron esetben a Schrödinger-egyenlet néha analitikusan megoldható.
Két vagy több elektronból álló rendszer esetén a Hamilton-
operátor sajátérték problémáját csak közelítések alkalmazásával, numerikusan tudjuk megoldani. Ebben a dolgozatban kétfajta iterációs eljárást vizsgálunk, amelyek a Schrödinger-egyenlet közelít® megoldását segítik el®. Az els® módszer a Schrödinger-egyenlettel elvben ekvivalens Bloch-egyenlet iterálásán alapul. A másik eljárás egyelektron-modellben használható, és idempotens egyrészecske s¶r¶ségmátrixok diagonalizálás-mentes megkeresésére irányul. Mindkét iterációs összefüggés tartalmaz egy tetsz®legesen megválasztható úgynevezett iterációs együtthatót (η ). Azt tapasztaltuk, hogy ezen együttható változtatásával az iteráció jellege (és eredménye) változik. Az η paraméter kis értékei esetén konvergál az iteráció, majd
η -t növelve két xpont közötti oszcillációt, majd bifurkációkat találunk. Tovább növelve az η paramétert, bizonyos molekulák esetén az iteráció kaotikussá válik. Végül, még nagyobb η -k esetén az eljárás divergál. Energiaszámítás szempontjából azok az η értékek megfelel®ek, melyek eredményeképpen az iteráció a lehet® leggyorsabban a legkisebb energiájú xpontba konvergál. Ebb®l számítható ugyanis a molekula alapállapotának energiája. Feladatom egy olyan eljárás kidolgozása volt, melynek segítségével el lehet dönteni, hogy egy adott η paraméter esetén az iteráció konvergál-e.
6
3. KÁOSZELMÉLET [1, 2] A kaotikus iterációk megértéséhez segít, ha megpróbáljuk az iterációk során tapasztalt jelenségeket a zikai káosz tapasztalatain keresztül értelmezni. Az ehhez szükséges legfontosabb fogalmakat ebben a fejezetben példák segítségével de niáljuk. Ha egy homorú felületre helyezünk egy golyót, akkor a golyó a felület alján kerül nyugalmi állapotba. A felület alja kitüntetett pont, mivel ott a golyó nyugalomban, egyensúlyban van; az ilyen pontokat xpontoknak hívjuk. Ha a golyót kicsit kitérítjük az egyensúlyi helyzetb®l, visszatér oda. Az ilyen típusú xpontot stabil xpontnak nevezzük. Analóg módon, azt az értéket, amelybe az iteráció konvergál, xpontnak nevezzük. Ha az iteráció kiindulási paraméterét kicsit megváltoztatva az eljárás ugyanahhoz az értékhez konvergál, az iteráció xpontját stabil xpontnak nevezzük. Helyezzük a golyót egy domború felületre. A golyó csak a felület tetején lehet nyugalmi állapotban, ezért ez a pont ezen felület egyetlen xpontja. Tudjuk azonban, hogy a golyó nem jut el minden esetben a xpontba, csak ha megfelel® helyr®l, megfelel® kezd®sebességgel indítjuk. Ha a golyót a xpontból kis távolságra elvisszük, a golyó nem tér vissza a xpontba, hanem egyre jobban eltávolodik onnan, legurul a lejt®n. Az ilyen típusú xpontokat instabil
xpontoknak nevezzük. Ennek megfelel®en, ha az iterációt a xpont környékér®l indítjuk, és az iteráció divergál, vagy más xpontba megy be, akkor a xpontot szintén instabilnak nevezzük. Egy acélgolyót tartalmazó inga stabil xpontja a függ®leges, kitérés nélküli állapot. Közelítsünk az acélgolyóhoz két mágnessel jobbról, illetve balról. Ha a két mágnes végtelen távol van a golyóhoz képest, a hatásuk nem érvényesül: a függ®leges kitérés nélküli állapot stabil marad. Ha a mágnesek elég közel kerülnek, a vonzás miatt két új xpont alakul ki a függ®leges állapothoz képest jobbra, illetve balra. A keletkez® két xpont stabil lesz, az eredetileg stabil állapot instabil xpont lesz. Ezt a jelenséget bifurkációnak nevezzük. Az, hogy a két xpont közül a golyó melyikbe megy, a kezdeti feltételeken, és véletlenszer¶, el®re nem kiszámítható hatásokon múlik. Amikor azonban a golyó elindult az egyik mágnes felé a mozgása dinamikailag jól leírható, jósolható, tehát nem kaotikus. Ennek analógjakénti azt, ha az iteráció két xpont között ugrál (oszcillál), szintén bifurkációnak nevezzük. Ha az iteráció eredményét ábrázoljuk a lépések számának a függvényében, jellegzetes villaszer¶
7 görbét kapunk. Ilyenre mutatunk példát kés®bb a 3. ábrán. A kaotikus mozgás kialakulásának szükséges és elégséges feltétele, hogy a rendszer végtelen sok instabil állapoton menjen keresztül. A káosz kialakulásának folyamata egy megfelel® paraméter változtatásával a következ®: az eredetileg stabil xpont instabillá válik, új xpontok jelennek meg, kialakulnak a bifurkációs villák. A paraméter továbi növelésével a mozgás periodicitása megsz¶nik, kaotikussá válik. A káosznak három alapvet® tulajdonsága van:
• nem ismétli önmagát • nem jelezhet® el®re, mert érzékeny a kezd®feltételekre, amelyeket soha nem ismerünk teljesen pontosan
• a visszatérési szabály bonyolult geometriájú: a hely - sebesség ábrázolásban (fázistér) szabályos szerkezet jelenik meg. Iterációk esetén hasonló a helyzet: az iterációs paraméter kis értékei esetén konvergens iteráció a paraméter növelésének hatására oszcillálva konvergál, majd további növelés után a xpont elveszti stabilitását, újabb xpontok alakulnak ki, az iteráció oszcillál, majd több bifurkációs villa alakul ki. Ha a paraméter értékét még tovább növeljük, a periodikusság megsz¶nik, az iteráció kaotikussá válhat, végül divergens lesz. A kaotikus iterációk egyik jól ismert példája a logisztikus-leképezés 1 :
x(n+1) = ηx(n) (1 − x(n) )
(1)
Ha η 6= 1 érték mellett megoldjuk az x = ηx(1 − x) egyenletet, azt kapjuk, hogy a leképezésnek két xpontja van az x = 0, illetve az x = 1 − 1/η helyen. Az x = 0 xpont instabil, ha
η > 1. A másik, x = 1 − 1/η xpont stabil 3 > η > 1 választás esetén. Ha η ertékét a (3,4) intervallumból választjuk, bifurkációt tapasztalunk, majd η = 4 esetben az iteráció kaotikus [2]. A kaotikus mozgás (x, t), vagy (v, t) térben ábrázolása nem alkalmas a mozgás a áttekintésére, mivel a periodicitás hiánya miatt mindig számíthatunk újabb viselkedésformákra. Tekintsünk egy N elemb®l álló rendszert. Ennek a rendszernek az (r1 , r2 , . . . , rN , v 1 , v 2 , . . . , v N ) 1
Fels® indexben az iterációs lépés sorszámát szerepeltetjük.
8 koordináták által kifeszített tere a fázistér.
Mivel ri = (xi , yi , zi ) és hasonlóképpen
v i = (vxi , vyi , vzi ), a fázistér 6N dimenziós. A fázistérben egy pont jeleníti meg a rendszer mozgásállapotát. A pont annak megfelel®en vándorol a fázistérben, ahogy a rendszer mozog. A mozgás fázistérbeli pályáját trajektóriának hívjuk. A trajektóriák a fázistérben nem metszhetik egymást. A fázistérben ábrázolva a kaotikus mozgást az ábrának szerkezete lesz. A fázistér azon pontjait, ahová a trajektóriák tartanak, azaz a fázistér vonzó pontjait
attraktoroknak nevezzük. Az attraktorokat az odatartó trajektóriák alakja alapján nevezzük el (például spirális-, hiperbolikus attraktor). Gyakorlatban a fázistérr®l metszeteket készítünk az ábrázolás lehet®vé tétele érdekében. Az iterációk esetén a (K, N ) teret nevezzük iterációs térnek (K az a mennyiség, amit az iteráció hivatott kiszámolni, N az iteráció sorszáma). Az iterációs térben ábrázolva egy kaotikus iterációt a valós térben ábrázolt mozgáshoz hasonló struktúrálatlan ábrát kapunk. Ennek kiküszöbölésére mátrix-iterációkra két, a zikai fázistérrel analóg teret de niálunk: i Ha ábrázoljuk az iterált K mátrix valamely elemét egy t®le független másik elem függvényében az iteráció során, az ábra szintén valamilyen geometriai alakzatot vesz. Ebben a dolgozatban ezt a teret kon gurációs térnek nevezzük. ii Ha a K mátrix valamely elemének változását ábrázoljuk az elem értékének függvényében, akkor a zikai fázistérre emlékeztet® teret kapunk, hiszen ez esetben egy változó "deriváltját" ábrázoljuk a változó függvényében. Ezért a jelen munkában ezt a teret nevezzük "fázistérnek". Mozgások esetén lehet®ség van xpontok stabilitásának vizsgálatára, a mozgás kaotikusságának el®rejelzésére. A Ljapunov-exponens (λ) a kaotikus mozgást végz®, eredetileg széttartó trajektóriák id®ben exponenciális széttartásának mértéke, azaz az el®rejelezhet®ség mér®száma. Ha λ < 0, a mozgás nem kaotikus, létezik valamilyen xpont; ha λ>0, nem létezik stabil xpont. Komplex λ értékek esetén a Ljapunov-exponens valós részének el®jele határozza meg a mozgás stabilitását, a fent leírtaknak megfelel®en. Iterációk esetén is lehet de niálni Ljapunov-exponenst, amely itt az iterációs tér közeli pontjainak távolodását méri. A távolság növekedésének mér®száma a K (N ) −K (f ixpont) ∼ eλN összefüggéssel de niált Ljapunov-exponens. A konvergencia el®rejelzése szempontjából hasonlóan m¶ködik, mint a
9 mozgásokra de niált Ljapunov-exponens: ha λ<0, az iterációnak létezik stabil xpontja, ha
λ>0 az iterációnak nem létezik xpontja. Ha λ = λ1 + iλ2 és λ2 6= 0, azaz a Ljapunovλn
exponens komplex, akkor λ komplex alakját beírva a de niáló egyenletbe, az e
= eλ1 n eiλ2 n
egyenletet kapjuk. A λ1 < 0 esetben az iteráció konvergál, ellenkez® esetben instabilitás lép fel. A képzetes rész oszcilláló eredményt ad, mely nem az iteráció stabilitásáról, hanem a xpont típusáról ad információt. A Bloch-egyenlet, és a s¶r¶ségmátrix-iteráció esetén az itt bevezetett fogalmakat fogjuk vizsgálni, illetve a Ljapunov-exponens segítségével el®rejelezzük az iteráció instabilitásának megjelenését.
10
4. A BLOCH-EGYENLET ITERATÍV MEGOLDÁSA 4.1. Az iterációs séma Egy molekula Schrödinger-egyenletét szeretnénk megoldani:
ˆ = EΨ HΨ
(2)
Legyen |Φi egy tetsz®leges 1-re normált függvény, amelyr®l csak azt kötjük ki, hogy ne legyen szingulárisan rossz közelítése Ψ-nek, azaz hΨ|Φi 6= 0. A szokásoktól eltér®en most a (2) egyenletben szerepl® egzakt Ψ hullámfüggvényt nem 1-re normáljuk, hanem kényelmi okokból a hΦ|Ψi = 1 úgynevezett közbüls® normálást írjuk el® [4].
ˆ operátort, amely Φ és Ψ között teremt kapcsolatot! Az Ω ˆ operátort De niáljunk egy Ω nevezzük hullámoperátornak:
ˆ Ψ = ΩΦ
(3)
ˆ =1−Ω ˆ Q
(4)
Vezessük be a
operátort is.
ˆ = |ΨihΦ| alakban. Ha ezt behelyettesítjük a (3) A hullámoperátor formálisan felírható Ω de niáló egyenletbe, akkor a Ψ = |ΨihΦ|Φi összefüggést kapjuk, amely Φ normáltsága miatt
ˆ hullámoperátor kivetíti egy Φ próbafüggvényazonosság. Azt mondhatjuk tehát, hogy az Ω b®l az egzakt hullámfüggvényt.
ˆ idempotens: Könnyen látható, hogy Ω ˆ 2 = |ΨihΦ|ΨihΦ| = Ω, ˆ Ω ˆ2 = Ω ˆ , a hullámoperátor mégsem nevezhet® szigorú értelemben mivel hΦ|Ψi = 1. Bár Ω ˆ + = |ΦihΨ| 6= Ω ˆ , azaz a hullámoperátor nem hermitikus. projektornak, mert Ω ˆ operátorra fennállnak Az idempotencia fontos következménye, hogy a (4) alatt de niált Q a
ˆΩ ˆ =0 Q
ˆQ ˆ=0 Ω
ˆ és Ω ˆ ortogonális idempotens operátorok. összefüggések, tehát Q
11
ˆ operátort meghatározó egyenletet nevezzük. Látni fogjuk, hogy Bloch-egyenletnek az Ω ennek egy lehetséges formája:
ˆH ˆΩ ˆ = 0, Q
(5)
ˆΩ ˆ −Ω ˆH ˆΩ ˆ = 0 alakban is írható. Löwdin ezt az egyenletet nemlineáris ami (4) miatt a H Schrödinger egyenletnek nevezi [5], de felfogható a Kvasni£ka-Lindgren-féle Bloch-egyenlet speciális esetének is [6, 7]. Az (5) egyenlet bizonyítása céljából a (2) egyenletet szorozzuk meg balról hΦ|-vel. Ekkor a
ˆ hΦ|H|Ψi =E
(6)
egyenletet kapjuk, mivel a közbüls® normálás miatt hΦ|Ψi = 1. Helyettesítsük vissza a (6) egyenletet a (2) egyenletbe:
ˆ ˆ H|Ψi = |ΨihΦ|H|Ψi
(7)
A kapott egyenletet szorozzuk meg jobbról hΦ|-vel:
ˆ |ΨihΦ| = |ΨihΦ| H ˆ |ΨihΦ| H | {z }
| {z }
| {z }
ˆ Ω
ˆ Ω
ˆ Ω
ˆ = |ΨihΦ|, a fenti egyenlet ekvivalens az (5) egyenlettel, ami bizonyítandó Tekintve, hogy Ω volt. Írjuk át az (5) egyenletet rekurzív formára, így egy iterációs képlethez jutunk:
ˆ (n+1) = Ω ˆ (n) + η Q ˆ (n) H ˆΩ ˆ (n) Ω
(8)
Itt η egy tetsz®leges paraméter, amelynek alkalmas megválasztása befolyásolja az iteráció konvergenciáját. Ennek az iterációs sémának fontos tulajdonsága, hogy meg®rzi a hullámo³
ˆ (n) perátor idempotenciáját. Felhasználva, hogy Ω ³
ˆ (n+1) Ω
´2
³
ˆΩ ˆ (n) ˆ (n) + η Q ˆ (n) H = Ω
´³
´2
ˆ (n) =Ω ´
ˆΩ ˆ (n) = ˆ (n) + η Q ˆ (n) H Ω
ˆ (n) + η Q ˆ (n) H ˆΩ ˆ (n) + η Ω ˆ (n) Q ˆ (n) H ˆΩ ˆ (n) + η 2 Q ˆ (n) H ˆΩ ˆ (n) Q ˆ (n) H ˆΩ ˆ (n) = = Ω |
{z
}
|
0 (n+1) ˆ
{z
}
0
ˆ (n) + η Q ˆ (n) H ˆΩ ˆ (n) = Ω = Ω
Fennáll továbbá, hogy az iteráció meg®rzi a hullámoperátor spurját is:
³
´
ˆ (n) ˆ (n) + η Tr ˆΩ ˆ (n) Q ˆΩ ˆ (n) = TrΩ ˆ (n+1) = TrΩ ˆ (n) + η Tr Q ˆ (n) H ˆ (n) , TrΩ H = TrΩ |
{z 0
}
12 hiszen a mátrixszorzat spurja invariáns a mátrixok sorrendjére. Az energiaszámításhoz vessük egybe a (2) és a (3) egyenleteket:
ˆ Ω|Φi ˆ ˆ H = E Ω|Φi Szorozzuk meg balról az egyenlet mindkét oldalát a hΦ| vektorral, ekkor az energiára a következ® összefüggést kapjuk:
ˆ Ω|Φi, ˆ E = hΦ|H
(9)
ˆ mivel hΦ|Ω|Ψi = hΦ|Ψi = 1. Vegyük észre a hasonlóságot az (1) logisztikus leképezés, és a (8) egyenlet közötti analógiát. Ez a nagyfokú hasonlóság okozhatja a Bloch-egyenelet iterációjának kaotikusságát.
4.2. Ljapunov-exponens számítás Az iteráció stabilitásának vizsgálatához Ljapunov-exponenst kell számolni. Ehhez tekintsük a (8) egyenlet mátrix reprezentációját:
Ω(n+1) = Ω(n) + ηQ(n) HΩ(n)
(10)
Az iteráció n-edik lépésében kapott Ω(n) mátrixot írjuk fel Ω(n) = Ω + ξ (n) alakban, ahol
Ω a xpont, ξ (n) pedig a xponttól való kis eltérés. Ekkor a (10) egyenlet a következ®képp alakul: ³
´
³
´
ξ (n+1) = ξ (n) + η 1 − Ω − ξ (n) H Ω + ξ (n) = ³
= ξ (n) + η HΩ + Hξ (n) − ΩHΩ − ξ (n) HΩ − ΩHξ (n) + O(2)
´
(11)
A ξ -ben kvadratikus tagot elhanyagoljuk, azaz linearizáljuk az egyenletet. Mivel Ω a xpont, HΩ = ΩHΩ. Vegyük fel a ξ (n) korrekciós mátrixot ξ (n) = ξeλn alakban, ahol
λ a Ljapunov-exponens. Ha ezt alkalmazzuk a ξ (n+1) mátrixra, akkor a ξ (n+1) = ξ (n) eλ összefüggést kapjuk. A kapott eredményt írjuk be a (11) egyenletbe:
ξeλ = ξ + η (Hξ − ξHΩ − ΩHξ) Vezessük be az m := eλ jelölést. Így a probléma egy négyindexes mátrix sajátérték egyenleteként fogható fel:
Aξ k = mk ξ k ,
(12)
13 ahol k a sajátvektorokat indexeli. Írjuk ki a mátrixszorzásokat: X
k k Aµν,λσ ξλσ = mk ξµν
λσ
(Aξ)µν = ξµν + η
à X
Hµλ ξλν −
λ
X
ξµλ Hλσ Ωσν −
λσ
X
!
Ωµλ Hλσ ξσν
λσ
Indexcserék után a fenti egyenlet helyett az alábbi írható:
(Aξ)µν =
X λσ
=
X
δµλ δνσ ξλσ + η "
à X
δσν Hµλ ξλσ −
λσ
X
δµλ ξλσ Hστ Ωτ ν −
σλτ
δµλ δνσ + ηδσν Hµλ − η
X
δµλ Hστ Ωτ ν − η
τ
σλ
X λστ
X
!
δνσ Ωµτ Hτ λ ξλσ
=
#
δνσ Ωµτ Hτ λ ξλσ
τ
Az A négyindexes mátrix elemeit tehát a következ® összefüggés adja meg:
Aµν,λσ = δµλ δνσ + ηδσν Hµλ − ηδµλ (HΩ)σν − ηδσν (ΩH)µλ Az A mátrix nem szimmetrikus, de a tapasztalat szerint sajátértékei valósak, így általánosított sajátvektor keres® eljárásokkal könnyen diagonalizálható. Erre a feladatra a WilkinsonReinsch ismertette eljárást használtuk [8]. Ha megoldjuk az A mátrix sajátérték problémáját, az eredményül kapott m sajátértékb®l a Ljapunov exponenst a λ = ln |m| összefüggéssel számolhatjuk. Ha ugyanis m < 0, a Ljapunov-exponens a
λ = ln |m| + iπ összefüggéssel adható meg. A fenti egyenlet triviális átalakítás után a
eλ = eiπ eln|m| alakba írható. Tudjuk, hogy eiπ = −1, tehát a − |m| = m összefüggésre jutottunk. Így végeredményül azt kaptuk, hogy Re λ = ln |m|
(13)
Az iteráció stabilitásához az ln |m| < 0 feltételnek teljesülnie kell, ami |m| < 1 esetben következik be.
14
4.3. Példák 4.3.1. A számítás kiindulási paraméterei
Minden iteráció elején kell adni egy kiindulási értéket az iterálandó mennyiségnek.
ˆ (0) -lal. A konkrét számításokhoz szükségünk van az Ω ˆ (0) kiindulási hullámoJelöljük ezt Ω perátor mátrixreprezentációjára. Ehhez a következ® meggondolásokat tegyük. Vegyük fel az
ˆ (0) operátort Ω ˆ (0) = |ΦihΦ| alakban, ahol Φ egyszer¶en az els® egységvektor (ez felel meg Ω például a Hartree-Fock hullámfüggvénynek a determinánsok bázisán egy CI számításban):
Φ :=
1
0
0 .. .
0 Egy operátor mátrixának (i, j) eleme egy vektorrendszer i-edik és j -edik elemével vett skalárszorzatával egyenl®. Ennek, illetve Φ fenti megválasztásának gyelembevételével a hullámoperátor mátrixelemeit a következ® meggondolással számíthatjuk ki: (0) ˆ Ωij = hi|Ω|ji = hi|ΦihΦ|ji = hi|1ih1|ji = δij δi1
A kiindulási mátrixot tehát
1 0 (0) Ω = .. .
0 0 ... 0 0 .. .
0 ... .. . ...
0 .. .
0 0 0 ... 0
alakban vettük fel. Az energiát a (9) egyenlet alapján a következ®képp számíthatjuk:
E=
X
H1k Ωk1 = (HΩ)11
k
4.3.2. Hidrogén molekula
Ahhoz, hogy számításokat tudjunk végezni, szükségünk van a vizsgálandó molekula Hamilton operátorának mátrixára is. Ezt egy 3-21G bázison történt ab initio kon gurációs
15 kölcsönhatás (CI) számítás során kaptuk, és a Mungauss program Budapest verziójával [9, 10] számoltuk ki. Az I. táblázatban összefoglaltuk a H2 molekula Bloch-egyenlettel történ® iterációjának eredményét η függvényében. Látható, ha η abszolút értéke egy kritikus határon túl n®, az iteráció nem konvergál, hanem egyre kevésbé lesz stabil, végül divergál. Nem érdemes azonban túl kis abszolút érték¶ η -t sem választani, ekkor ugyanis az iteráció csak lassan konvergál. Jelen példában az η = −0, 1 esetben az iteráció 150 lépés alatt, míg η = −0, 5 választás esetén 40 lépés alatt konvergált. Az iteráció η = −0, 600 érték körül a megszokottól eltér® viselkedést mutat (1. ábra). Ha
η − 0, 601 és −0, 604 közé esik az iteráció egyre lassabban konvergál, és oszcillálva közelíti meg a xpontot. Az η = −0, 605 értéknél az iteráció látszólag nagyon gyorsan, 40 lépés alatt bekonvergál. A 400. lépés után azonban a 8. értékes jegybe bizonytalanság kerül, az iteráció kicsit oszcillál (2. ábra). Az η = −0, 607 esetben a kezeti látszólagos konvergencia után egyre távolodnak egymástól az energia értékek (2. ábra). Mivel az energiát 6 tizedes pontossággal elég számolnunk, ezek az η értékek tulajdonképpen gyakorlati szempontból kedvez®ek. Az η paraméter abszolút értékét tovább növelve az iteráció energia számításra alkalmatlanná válik a bifurkáció és a káosz megjelenése miatt (I. táblázat). Az iteráció η < 0 esetben a legalcsonyabb CI energiához, azaz az alapállapot energiájához konvergált, η > 0 esetben pedig a CI számítás legmagasabb szintjéhez tartozó energiát kaptuk meg. Mivel általában az alapállapotot számítjuk, legtöbbször negatív η értékeket adunk meg. Az A mátrix diagonalizálását elég csak a független paraméterekre kell elvégezni. A független paraméterek számát az Ω mátrix független elemeinek száma adja. Az Ω elemeire megszorítást jelent, hogy Ω idempotens és a nyoma konstans (esetünkben 1). Ezen túl megszorítást jelent Ω-ra az is, hogy egyértelm¶en származtatható kell legyen az M részecskés rendszer közelít®, illetve egzakt |Φi és hΨ| hullámfüggvényeib®l. Ez utóbbi feltétel az úgynevezett M-reprezentábilitási feltétel. Mivel ezen feltételekb®l következ® összefüggés az Ω mátrixelemei között általában nem ismeretes, a redundanciák kiküszöbölésével empírikusan próbálkoztunk az alábbi módokon:
16 1. csak a µ < ν és λ < σ alsó háromszöget vesszük gyelembe 2. csak az alsó háromszöget és a f®átlót (µ ≤ ν és λ ≤ σ ) vesszük gyelembe Ezen kívül foglalkoztunk a teljes n2 ∗ n2 -es A mátrixszal, azaz 3. teljes mátrixot azaz ∀µ, ν és ∀λ, σ értéket gyelembe vesszünk
I. táblázat. Energia számítás a H2 molekulára a Bloch-egyenlet segítségével Az iteráció jellegét az
iterációs térben ábrázolt pontok mintázata alapján állapítottuk meg. Oszcillációról egy bifurkációs villa megjelenésekor beszélünk. Kaotikus iteráció esetén z¶rzavaros mintázatot látunk, a bifurkációs villák nem azonosíthatók.
η
E/Eh
jelenség
-0,100 -0,024769473
konvergál 158 lépésben
-0,500 -0,024769473
konvergál 44 lépésben
-0,600 -0,024774342
konvergál 500 lépésben
-0,605
-
oszcilláció
-0,700
-
oszcilláció
-0,770
-
bifurkáció (4 villa)
-0,776
-
bifurkáció (7 villa)
-0,777
-
káosz
-0,800
-
káosz
-0,900
-
káosz
0,100 3,28484366
konvergál 178 lépésben
0,500 3,28484366
konvergál 43 lépésben
0,700
-
két xpont között oszcillál
0,770
-
bifurkáció (4 villa)
0,776
-
bifurkáció (7 villa)
0,777
-
káosz
0,800
-
káosz
0,900
-
divergál
17 1. ábra. A H2 molekula Bloch egyenlettel történ® iterációjával nyert energiák konvergens η -k esetén. A görbék mellett feltüntettük a hozzájuk tartozó η értékeket, N az iterációs lépés sorszáma. -0.018
-0.02 -0.604 -0.022 -0.601
E/Eh
-0.024 -0.605 -0.026 -0.602 -0.028
-0.03
-0.032 100
200
300
400
500
600
N
A következ®kben a fontos mennyiségek (m, λ, η ) számozása ennek a számozásnak megfelel®en történik. A fenti három eset mindegyikében azt tapasztaltuk, hogy az A spektrumában megjelenik a degeneráció, és a legnagyobb sajátértékek 1,0000-gyel egyenl®ek. Konvergens iterációt akkor kapunk, ha a Ljapunov-exponens valós része negatív ami az
m ∈ (−1, 1) tartományban teljesül. Ez a feltétel az A összes sajátértékére fenn kell álljon. Példáinkban a legnagyobb sajátérték egyszer sem vált 1-nél nagyobbá. A legkisebb sajátérték és a hozzá tartozó Reλ = ln |m| számot a II. táblázatban foglaltuk össze, különböz® η iterációs paraméterekre. A II. táblázatból kit¶nik, hogy m1 azon η értéknél lép ki el®ször a megadott intervallumból, amikor az iteráció egy kicsit elkezdett oszcillálni. (I. táblázat). A másik két módszer csak nagyobb abszolút érték¶ η esetén jelez. A Bloch-egyenlet lehet®séget ad arra is, hogy egyszerre több sajátvektort meghatározzunk. Ehhez az Ω(n) mátrix megfelel® oszlopába a gerjeszetett állapot próbavektorát kell
18 2. ábra. A H2 molekula Bloch egyenlettel történ® iterációja. A görbék mellett feltüntettük a hozzájuk tartozó η értékeket, N az iterációs lépés sorszáma. A három görbe közül az η = −0, 605 eset látszólag konvergens iterációt mutat; η = −0, 605, illetve -0,607 esetén nemkonvergens iterációt látunk. -0.02476944
-0.02476945
-0.02476946
-0,606
-0,605
E/Eh
-0.02476947
-0.02476948
-0.02476949
-0,607
-0.02476950
-0.02476951 700
720
740
760
780
800
N
II. táblázat. Az A mátrix sajátértékei és a Ljapunov-exponensek valós része a H2 molekulára. Az mi , λi mennyiségek i indexe a szövegben de niált esetek sorszáma η
m1
λ1
m2
λ2
m3
λ3
-0,100 0,66904 -0,37048 0,67152 -0,39821 0,71632 -0,33362 -0,500 -0,65481 -0,42341 -0,64242 -0,44251 -0,41838 -0,87136 -0,605 -1,00232 0,00231 -0,98733 -0,01275 -0,71624 -0,33373 -0,630 -1,08506 0,08163 -1,06945 0,06714 -0,78716 -0,23932 -0,705 -1,33280 0,28728 -1,31581 0,27445 -0,99991 -0,00009
írni. Ha például a
0
1 Φ1 := 0 . ..
0
19 (0)
az els® gerjesztett állapot próbavektora, akkor az Ω22 elem nem 0-nak, hanem 1-nek adódik. Az els® gerjesztett állapot energiáját az E2 =
P k
H2k Ωk2 képlettel tudjuk számolni. Hason-
lóan, ha Ω(0) f®átlójában az (i, i) elemet 1-re cserélem akkor az i-edik állapot energiáját is ki lehet számolni. A gyakorlatban azonban ezen cserék után az iteráció többnyire hamis xpontokhoz konvergál, így hamis energiaértékeket kapunk. Ezért a Bloch-egyenlet iterálása ebben a formában nem bizonyult jónak gerjesztett állapotok energiáinak számítására.
4.3.3. Hélium atom
A He atom esetén 6-311G bázison történt a CI számítás. Az iteráció eredményeit az η változtatásával a III. táblázat foglalja össze, és a 3., 4. ábrák illusztrálják. Hasonlóan a H2 -molekulán végzett számításokhoz, ebben a példában is kisebb abszolút érték¶ η -k esetén konvergál az iteráció, majd növelve a paramétert ismét instabillá válik a megoldás: bifurkáció, s®t káosz is fellép. Számítások szempontjából lényeges különbség azonban, hogy az iteráció η < −0, 155 esetben nem konvergál, tehát lényegesen kisebb intervallumból válogathatunk alkalmas η -kat. III. táblázat. Energia számítás a He atomra a Bloch-egyenlet segítségével. η
E/Eh
jelenség
-0,100 -0,0165229 konvergál 119 lépésben -0,155 -0,0165229 konvergál 260 lépésben -0,156
-
oszcilláció
-0,190
-
oszcilláció
-0,198
-
bifurkáció (4 villa)
-0,200
-
káosz
-0,240
-
divergál
0,100 12,8368518
konvergál
A 3. ábrán az η paraméter 4 különböz® értéke esetén ábrázoltuk a számított energiát az iteráció számának (N ) függvényében. Az iteráció konvergál η = −0, 155 esetben, majd η -t csökkentve egyre jobban távolodik a xponttól, és egyre bonyolutabb oszcillációt végez. A
20 4. ábra a káosz kiszámíthatatlanságát mutatja: η értékét mindössze 0,01-dal változtattuk meg és az iteráció látványosan más értékeket adott. A Ljapunov-exponens számítás eredménye is lényegesen különbözik a H2 molekulára végzett számolások alapján várttól (IV. táblázat). Míg az el®z® példában m1 jóslata vált be nagyon jó közelítéssel, most az m2 , illetve m3 módszer bizonyult megbízhatóbbnak (ezek ezred pontossággal megjósolták az oszcilláció kritikus η értékét). IV. táblázat. Az A mátrix sajátértékei és a Ljapunov-exponensek valós része a He atomra. η
m1
λ1
m2
λ2
m3
λ3
-0,100 0,23753 -1,43746 -0,28369 -1,25987 -0,28534 -1,25407 -0,130 0,00880 -4,73300 -0,66879 -0,40228 -0,67094 -0,39908 -0,156 -0,18945 -1,66363 -1,00255 0,00255 -1,00513 0,00512 -0,263 -1,00528 0,00527 -2,37609 0,86546 -2,38044 0,86728
A következ®kben látni fogjuk, hogy a kaotikus iterációkra jellemz® struktrálatlan ábrák, mint a 4a, 4b, a kon gurációs, illteve a fázistérben milyen struktúrált formát öltenek. Ezt azonban már nem a Bloch-egyenlet példáján, hanem egy gyakorlati szempontól fontosabb esetben, az els®rend¶ s¶r¶ségmátrix iterációs formuláján mutatjuk be.
21 3. ábra. A hélium atom Bloch egyenlettel történ® iterációjának instabillá válása az η iterációs paraméter függvényében.
7 6 -0,190
5 4 3
E/Eh
2 -0,157
1 -0,155
0 -1 -0,198 -2 -3 -4 0
50
100
150
200
250
300
350
400
450
N
4. ábra. A kaotikus iteráció ábrázolása η =-0,21 (a) illetve η =-0,22 (b) értékek esetén. 10
12
8
10 8
6 6 E/Eh
E/Eh
4 4
2 2
0 0
-2
-2 -4
-4 0
(a)
100
200
300
400 N
500
600
700
800
0
100
(b)
200
300
400 N
500
600
700
800
22 .
23
5. A HARTREE-FOCK SRSÉGMÁTRIX ITERÁCIÓJA 5.1. Az iterációs séma A Hartree-Fock-Roothan-elméletben a φi molekulapályákat a molekulát felépít® atomok atompályáinak lineárkombinációjával kapjuk meg:
φi =
X
ciµ χµ ,
µ
ahol χµ az atompálya. A ciµ koe cienseket a Roothan-egyenletek segítségével kaphatjuk meg:
Fci = ²i Sci
(14)
Itt F a Fock-mátrix. Az S átfedési mátrixra azért van szükség, mert az atompályák bázisa nem ortogonális (ebben a fejezetben végig spinpályákkal dolgozunk). A Fock-mátrix elemei a hµν egy-, és a (µν|λσ) kételektron-integrálokból az alábbi módon épülnek fel [4]:
Fµν = hµν +
X
Pµν [(µν|λσ) − (µσ|λν)] ,
(15)
µν
ahol a P egyrészecske s¶r¶ségmátrixot a ci koe ciensek segítségével építhetjük fel:
Pµν =
occ X
ciµ ciν
(16)
i
A P mátrix legfontosabb tulajdonságai:
P+ = P
(17)
Tr (PS) = M
(PS)2 = PS ahol M a molekula elektronjainak a száma. Ha a bázis ortonormált (S = 1), amit a továbbiakban - az általánosság megsértése nélkül - felteszünk, akkor a (17) tulajdonságok a
P+ = P, Tr P = M, P2 = P
(18)
képletekre egyszer¶södnek. Ezek szerint a P mátrix hermitikus, idempotens, tehát projektor, és spurja megadja az elektronok számát. A P mátrix projekciós jellege annak felel meg,
24 hogy - mint ezt a (16) egyenlet mutatja - P a betöltött molekulapályák alterére vetít. Ennek megfelel®en sajátértékei az 1 és 0 számok, betöltött pályára 1, virtuálisokra 0. A molekula energiája, mint a P mátrix funkcionálja az alábbi képlettel adható meg [11]:
1 E = Tr [(H + F) P] , 2
(19)
ahol H elemei a hµν egyelektron integrálok . Ez Hückel-típusú rendszerekre, amelyekben
(µν|λσ) = 0, az E = Tr (HP)
(20)
alakra egyszer¶södik. A P mátrix elemeit nemcsak a (16) egyenlet segítségével, hanem több más közvetlen eljárással is meg lehet határozni [1221]. Ezek lényege, hogy elkerüljük a (14) sajátértékegyenlet megoldását, ami nagy bázison rendkívül költséges, hiszen a mátrix méretének harmadik hatványával arányos számítási munkát igényel. A P mátrixot a ci koe ciensek nélkül lényegében háromféle eljárással lehet meghatározni: 1. Az F mátrix G(z) = (z − F)−1 Green-függvényének a
P=
1 I G(z)dz 2πi C
kontúrintegráljával, ahol a C az úgynevezett Coulson-kontúr a betöltött pólusokat veszi körbe [12, 13]. Ennek a módszernek hátránya, hogy a (z − F) mátrix invertálását igényli. 2. A (19) (vagy(20)) energiaképletek közvetlen minimalizálásával [15, 16]. Ekkor a (17) vagy (18) tulajdonságok, mint variációs mellékfeltételek veend®k vigyelembe. 3. Iterációs eljárásokkal. Ezek közül a McWeeny-féle puri kációs eljárást [22], és az ezzel gyakorlatilag ekvivalens NémethScuseria-féle el®jel-mátrix technikát [21] említjük meg. Ezek legtöbbször gyorsan konvergálnak, de több esetben (például fémes szerkezetekre) nem alkalmazhatóak. Ebben a dolgozatban egy nemrégen javasolt, a 3. csoportba tartozó iterációs technikát vizsgálunk, amelynek lényege a következ® (a képleteket az egyszer¶ség kedvéért egyelektron
25 (Hückel) szinten adjuk meg, de az eljárás minden további nélkül általánosítható ab initio Hartree-Fock szintre is): Az ortogonális bázison reprezentált egzakt s¶r¶ségmátrix kielégíti a
[H, P] = 0 egyenletet [22, 23], amit jobbról P -vel szorozva a HP − PHP = 0 Bloch-típusú egyenletet kapjuk. Vagyis (21)
QHP = 0,
ahol Q = 1 − P az úgynevezett lyuks¶r¶ség-mátrix. A de nicióból látszik, hogy QP = 0. A Bloch-egyenletnél tárgyaltaknak megfelel®en a P mátrixra a ³
P(n+1) = P(n) + η Q(n) HP(n)
´
(22)
iteratív összefüggést írhatjuk fel. A korábbihoz hasonlóan látható, hogy a fenti képlettel végrehajtott iteráció meg®rzi P nyomát:
³
´
´
³
³
´
³
´
³
Tr P(n+1) = Tr P(n) + η Tr Q(n) HP(n) = Tr P(n) ³
´
´
mivel Tr Q(n) HP(n) = Tr P(n) Q(n) H = 0. Az iteráció meg®rzi a P idempotenciáját is: ³
P(n+1)
´2
³
= P(n)
´2
+ ηP(n) Q(n) HP(n) + ηQ(n) HP(n) P(n) + ηQ(n) HP(n) Q(n) HP(n) =
= P(n) + ηQ(n) HP(n) = P(n+1) Az idempotencián és a nyomon kívül az iterációnak a hermiticitást is meg kell ®riznie, ez fontos különbség az Ω-iterációhoz képest. Ennek vizsgálatához vegyünk fel egy közelít® P(0) s¶r¶ségmátrixot és a hozzá tartozó Q(0) lyuk-mátrixot. Be fogjuk látni, hogy a (22) iteráció nem konvergálhat. Végezzünk indirekt bizonyítást! Tegyük fel, hogy az iteráció az egzakt
P -t adja eredményül. Ekkor P -t a P = P(0) + Q(0) hef f P(0)
(23)
alakban írhatjuk, ahol hef f az iterációk során felgy¶lt közbens® eredményeket foglalja magában. Szorozzuk meg a (23) egyenletet balról, majd jobbról P(0) -lal. Ekkor a P(0) P = P(0) illetve a PP(0) = P egyenleteket kapjuk. Ezen egyenleteket és P(0) illetve P hermiticitását
26 kihasználva a következ®ket írhatjuk: P(0) P = P és P(0) P = P(0) , amib®l az következik, hogy P(0) = P. Ezzel ellentmondásra jutottunk a kiindulási P(0) közelít® jellegét illet®en. Az ellentmondás azzal oldható fel, hogy a (22) iteráció nem ®rzi meg a hermiticitást a (23) alakban. Ezért tehát P a (23) alakban nem lehet egzakt. A problémán a következ®képp segíthetünk: Vegyük észre, hogy a P(n+1) = P(n) + ηP(n) HQ(n) képlettel végzett iteráció ugyanúgy meg®rzi P nyomát és idempotenciáját, mint (22). Ezek alapján a s¶r¶ségmátrixot az alábbi dupla-iterációval számíthatjuk:
P(n+1) = P(n)
+ ηQ(n) HP(n)
(24)
P(n+2) = P(n+1) + ηP(n+1) HQ(n+1) Felhasználva, hogy Q(n+1) = 1 − P(n+1) = 1 − P(n) − ηQ(n) HP(n) = Q(n) − ηQ(n) HP(n) , a (24) két egyenletét kombinálva a ³
´
³
P(n+2) = P(n) + ηQ(n) HP(n) + η P(n) + ηQ(n) HP(n) H Q(n) − ηQ(n) HP(n)
´
(25)
összefüggést kapjuk. Mivel (24) egyenlet mindkét lépése meg®rzi a spurt és az idempotenciát, a (25) egyenlettel végzett iteráció is rendelkezik ezekkel a tulajdonságokkal. Konvergencia esetén a QHP = 0 és a PHQ = 0 feltételek teljesülnek, tehát H és P kommutál. Ezekb®l a konvergált P hermiticitása automatikusan következik [3].
5.2. Ljapunov-exponens számítás Ljapunov-exponens számításhoz a 4.2 fejezetben tárgyaltakhoz analóg meggondolásokat kell tennünk. A különbség abból adódik, hogy a P-iterációnak a hermiticitást is ki kell elégítenie, ezért a (25) egyenlet szerint kell az iterációt végezni. Írjuk fel a (25) egyenletet összegalakban:
P(n+1) = P(n) + ηQ(n) HP(n) + ηP(n) HQ(n) +
(26)
+ η 2 Q(n) HP(n) HQ(n) − η 2 P(n) HQ(n) HP(n) − − η 3 Q(n) HP(n) HQ(n) HP(n) Végezzük el a fenti egyenletben a P(n) = P+Π(n) , illetve a Q(n) = Q−Π(n) helyettesítéseket.
P az egzakt s¶r¶ségmátrix (az iteráció xpontja), Q pedig az ebb®l számított lyuk-mátrix.
27 Mivel P(n+1) = P + Π(n+1) , a (26) egyenlet a következ®képp alakul:
Π(n+1) = Π(n) + ηQHP − ηΠ(n) HP + ηQHΠ(n) + ηPHQ + ηΠ(n) HQ −
(27)
− ηPHΠ(n) + η 2 QHPHQ − η 2 Π(n) HPHQ + η 2 QHΠ(n) HQ − − η 2 QHPHΠ(n) − η 2 PHQHP − η 2 Π(n) HQHP + η 2 PHΠ(n) HP − − η 2 PHQHΠ(n) − η 3 QHPHQHP + η 3 Π(n) HPHQHP − η 3 QHΠ(n) HQHP + + η 3 QHPHΠ(n) HP − η 3 QHPHQHΠ(n) + O(Π2 ) A (21) egyenlet értelmében QHP = PHQ = 0, ezért a (27) egyenlet linearizálás után a
Π(n+1) = Π(n) − ηΠ(n) HP + ηQHΠ(n) + ηΠ(n) HQ − − ηPHΠ(n) + η 2 QHΠ(n) HQ + η 2 PHΠ(n) HP összefüggéssé egyszer¶södik. Végezzük el a Π(n) = Π(0) eλn helyettesítést. Mivel
Π(n+1) = Π(0) eλ(n+1) = Π(n) eλ , eλn -nel való osztás után a Π(0) eλ = Π(0) − ηΠ(0) HP + ηQHΠ(0) + ηΠ(0) HQ −
(28)
− ηPHΠ(0) + η 2 QHΠ(0) HQ + η 2 PHΠ(0) HP egyenletet kapjuk. Bevezetve az m := eλ jelölést a (12) egyenlettel analóg hipermátrix sajátérték problémát nyerünk:
AΠ(0) = mΠ(0) , ami ekvivalens a
X
(0)k
Aµν,λσ Πλσ = mk Π(0)k µν
λσ
egyenlettel, ha a mátrixszorzások indexeit beírjuk.
28 Részletesen, (28) alapján ³
AΠ(0)
´ µν
= Π(0) µν − η | {z }
|
X
{z
Qµλ Hλσ Π(0) σν + η {z
}
|
X
(0)
Πµλ Hλσ Qσν − η
λσ
3
X
{z
Pµλ Hλσ Π(0) σν +
λσ
}
4
5 (0) Pµλ Hλσ Πστ Hτ ξ Pξν ,
X
}
|
|
X
{z
2 Qµλ Hλσ Π(0) στ Hτ ξ Qξν + η
λστ ξ
|
(29)
}
2
λσ
+ η2
(0)
Πµλ Hλσ Pσν +
λσ
|
1
+ η
X
λστ ξ
{z
}
{z
6
}
7
ahol a tagokat a kés®bbi azonosítás céljából számoztuk be. Az A mátrix (µν, λσ) elemének (0)
meghatározásához indexcserék után a (29) egyenletb®l kiemelünk jobbról Πσλ -t. Az indexcseréket a következ®képp hajthatjuk végre: 1.
X
(0)
Πλσ δµλ δνσ =
X
λσ
(0)
δµσ δνλ Πσλ
σλ
2.
− η
X λσα
= −η
(0)
δµα Παλ Hλσ Pσν = −η X
X
(0)
δµσ Πσλ Hλ,α Pαν =
λσα (0) (HP)λν δµσ Πσλ
σλ
3. X
η
Qµλ Hλσ Π(0) σν = η
λσ
4.
η
X
X
(0)
δνλ (QH)µσ Πσλ
σλ
(0)
Πµλ Hλσ Qσν = η
λσ
X
(0)
δµσ (HQ)λν Πσλ
σλ
5.
− η
X λσ
= −η
Pµλ Hλσ Π(0) σν = −η X λσα
X
δνα Pµλ Hλσ Π(0) σα =
λσα (0) δνλ Pµα Hασ Πσλ = −η
X σλ
(0)
δνλ (PH)µσ Πσλ
29 6.
η2
X
2 Qµλ Hλσ Π(0) στ Hτ ξ Qξν = η
X
λστ ξ
= η
2
X
(0)
(QH)µσ Πσλ (HQ)λν =
σλ (0) (QH)µσ (HQ)λν Πσλ
σλ
7.
η2
X
2 Pµλ Hλσ Π(0) στ Hτ ξ Pξν = η
³
³
´
(0)
(PH)µσ (HP)λν Πσλ
λσ
λστ ξ
A AΠ(0)
X
(0)
µν
AΠ(0)
kifejezés Πσλ kiemelése után a
´ µν
=
X³
δµσ δνλ − ηδµσ (HP)λν + ηδµσ (QH)µσ + ηδµσ (HQ)λν −
σλ
´
(0)
− ηδνλ (PH)µσ + η 2 (QH)µσ (HQ)λν + η 2 (PH)µσ (HP)λν Πσλ alakban írható. Ebb®l A elemeire a következ® összefüggést kapjuk:
Aµν,λσ = δµλ δνσ − η(HP)σν δµλ + η(QH)µλ δνσ + η(HQ)σν δµλ −
(30)
− η(PH)µλ δνσ + η 2 (QH)µλ (HQ)σν + η 2 (PH)µλ (HP)σν A Ljapunov-exponenst most is az A hipermátrix m sajátértékeinek logaritmusa szolgáltatja.
5.3. Példák 5.3.1. Számolás ab initio Hartree-Fock szinten
A Hartree-Fock szint¶ P-iterációt 60 db víz molekula energiájának kiszámításával demonstráljuk (V. táblázat). Modellünkben a vízmolekulák lineáris láncot alkotnak, és az oxigén-oxigén távolság 3,09 Å. A számítást 6-31G bázison végeztük. A kapott energiaértékek nem nagyon változnak az η változtatásával, az ideális η értékr®l ezért az dönt, hogy mely esetben konvergált az iteráció a legkevesebb lépés után: ez alapján a legideálisabb az η = −0, 06 választás. Az ab initio Hartree-Fock szinten végzett P-iteráció el®nye, hogy bármilyen idempotens
P(0) mátrixból indítható, így az el®z® SCF ciklus során bekonvergáltatott P mátrixból is.
30 V. táblázat. 60 darab vízmolekula energiájának számítása Hartree-Fock szint¶ P-iterációval. A konvergenciához szükséges lépések száma NSCF .
η
E/Eh
NSCF
-0,05 -4557,781803
13
-0,06 -4557,781815
14
-0,08 -4557,781815
15
-0,15
divergál
Ez azért el®nyös, mert az el®z® ciklus P mátrixa egyre jobb közelítés, ahogy az SCF iteráció halad el®re, így a P-iteráció egyre kevesebb lépést igényel. Ezzel az el®nnyel nem minden s¶r¶ségmátrix-felépít® eljárás rendelkezik.
5.3.2. A Hückel-szint¶ számítás kiindulási paraméterei
Részletesebb elemzés céljából alkalmazzuk a P-iterációt a Hückel-közelítés bevezetésével. A közelítés konjugált kötésrendszer¶ molekulákra m¶ködik, lényege, hogy csak a π pályákat és az els® szomszéd kölcsönhatásokat vesszük gyelembe a Hamilton-mátrixban. Ennek megfelel®en
Hi,j = β , ha az i-edik és a j -edik atom között σ kötés van, Hi,j = 0, ha az i-edik és a j -edik atom között nincs σ kötés. A legegyszer¶bb esetben minden β mátrixelem egyforma, így alkalmas energiaegység választással a β = 1 egyszer¶sítéssel élhetünk. Így a H-mátrix például gy¶r¶k esetében
0 1 0
1 H= 0 . ..
0 1
0 ... 1 0 ... 0
1 0 1 ... .. . . . . . . . . . .
1 0 ... 0
0 .. .
1 0
alakú lesz. Nyílt láncú, konjugált kötésrendszer¶ molekulák esetén H1,N = HN,1 = 0, mivel a két széls® atom között nincs kötés.
31 Legyen az iteráció kiinduló P(0) mátrixa a molekula olyan Kekulé-mátrixa, amelynek diagonális elemei is 1-ek. A Kekulé-mátrix (i, j) eleme 1 ha az i-edik és a j -edik elem között kett®s kötés van, egyébként 0. A Kekulé-mátrixot a topologikus-mátrixból egy általunk írt program segítségével határoztuk meg.
1 1 0
1 0 (0) P = 0 .. .
1 0
0 ... 0 0 ... 0
0 1
1 ... 0
0 1 1 ... .. . . . . . . . . . .
0 0 ... 0
0 .. .
1 1
Az így kapott P(0) mátrix annak felel meg, mintha a konjugált rendszer egymástól független kett®s kötés¶ egységekb®l épülne fel. Ez természetesen nem jó közelítése a Hückel P mátrixnak, de matematikailag korrekt (azaz hermitikus és idempotens) kiindulóponot jelent az iteráció számára. A P-iteráció menetét egy nyílt láncú (butadién), és egy gy¶r¶s (benzol) konjugált rendszeren keresztül mutatjuk be. Végül szintén Hückel közelítésben modellszámítást végzünk egy nyílt láncú poliénre, illetve egy nanocs®re.
5.3.3. Butadién
A butadién, bár a legegyszer¶bb konjugált rendszer¶ molekula, lehet®séget ad arra, hogy az iteráció viselkedését az iterációs tér mellett a kon gurációs és a fázistérben is vizsgáljuk az η függvényében. Ha az energiát β egységekben adjuk meg, akkor a pozitívabb energiák felelnek meg a kedvez®bb állapotoknak, így a VI. táblázatban a pozitív η értékek vezetnek az alapállapot energiájához. Negatív η értékkel a legmagasaban gerjesztett állapothoz konvergál az iteráció. Ennek energiája csak el®jelben különbözik az alapállapot energiájától, ez az úgynevezett elektron-lyuk szimmetria következménye. A kapott eredmények hasonlók a Bloch-egyenlet példáinál tapasztaltakhoz: kis abszolútérték¶ η -k esetén az iterációnak egy xpontja van, majd növelve η értékét a xpont instabil lesz, újabb xpontok jelennek meg (5. ábra), majd
32 VI. táblázat. Energia számítás butadiénre a P-iteráció segítségével η
E/β -egység
jelenség
-0,100 -4,472135954 konvergál 664 lépésben 0,100 4,472135954 konvergál 408 lépésben 0,355 4,472135954 konvergál 2145 lépésben 0,356
-
oszcillációval konvergál
0,357
-
oszcillál
0,359
-
káosz
0,360
-
divergál
megjelenik a káosz (6. ábra), végül pedig a divergencia. 5. ábra. A butadién P-iterációja az iterációs térben ábrázolva η =0,357, illetve η =0,358 értékek esetén.
4.5 0,357
4.45
4.4
4.35
E
0,358 4.34 4.3 4.26 4.22 4.18 1200
4.3
4.25
1300
1400
4.2
4.15 0
500
1000
1500
2000
2500 N
3000
3500
4000
4500
5000
A kaotikus iterációk az iterációs térben ábrázolva struktúrálatlan képet nyújtanak (6. ábra), ezért célszer¶bb a 3. fejezetben bevezett kon gurációs-térbeli ábrát nézni. Ábrázoljuk ezért a s¶r¶ségmátrix P2,3 elemét a P1,2 függvényében az iteráció során (7., 8. ábra).
33 6. ábra. A butadién P-iterációja az iterációs térben ábrázolva η =0,359 paraméter érték esetén. 4.5 4.45 4.4
E
4.35 4.3 4.25 4.2 4.15 4.1 0
500
1000
1500
2000
2500 N
3000
3500
4000
4500
5000
Ez az ábrázolás a 3. fejezetben említett kon gurációs térben való ábrázolásnak felel meg. Vizsgáljuk meg alaposan a 7. ábrát! Az η =0,100 értékhez tartozó görbe a kon gurációs tér egy pontjába tart, a pontok ott torlódnak. Ez a görbe felel meg a konvergáló esetnek (vö: VI. táblázat). Ha megnézzük az η = 0,357 értékhez tartozó görbét, azt látjuk, hogy a pontok a kon gurációs tér két pontjába s¶r¶södnek, tehát ez oszcillációnak felel meg az iteráció terében. A másik két görbe kaotikus az iterációs térben. Figyeljük meg, hogy a két görbe hasonlít egymáshoz: mindkett® kis hurokból áll, melyek közül a kisebbik hurok metszi önmagát; η = 0, 3595 esetben pedig a két huroknak is van metszéspontja. A trajektóriák csak a kon gurációs térben metszhetik egymást, a fázistérben nem (9. ábra). A különbség tehát a két η értékhez tartozó görbe között lényegében annyi, hogy a nagyobb η -hoz tartozó görbén nagyobb hurkok vannak.
34 7. ábra. A butadién P-iterációja a kon gurációs térben ábrázolva η =0,100, η =0,357, η =0,359, illetve η =0,3595 értékek esetén.
0.9 0.8 0.7 0,357
0.6
P(2,3)
0.5 0,3595
0.4
0,359
0,100
0.3 0.2 0.1 0 -0.1 -0.2 -0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
P(1,2)
8. ábra. A butadién P-iterációja a kon gurációs térben ábrázolva η =0,358 érték esetén. 0.65
0.365 0.363
0.6
0.361 0.359
P(3,2)
0.55
0.5
0.646
0.357
0.642
0.355 1.101 1.103 1.105 1.107 1.109
0.638 0.634
0.45
0.63 0.37
0.39
0.41
0.43
0.45
0.4
0.35
0.4
0.5
0.6
0.7
0.8 P(1,2)
0.9
1
1.1
1.2
35 A 8. ábra az η =0,358 értékhez tartozó görbét mutatja, szintén a kon gurációs térben. A pontok nagyrésze a tér két pontján van, ez tehát egy oszcilláló iteráció (v.ö. 5. ábra). A kinagyított ábrákon látszik, hogy a xpontokba egy spirális pályán megy be a trajektória, ezért az ilyen típusú pontokat spirális attraktoroknak nevezik. Érdekesség, hogy a mechanikából jól ismert csillapított harmonikus oszcillátornak, melynek a valós térben ábrázolt görbéje exponenciálisan lecseng® szinusz-függvény, a fázistérben szintén spirális attraktora van [1]. Ha megnézzük az 5. ábrát, láthatjuk, hogy a P-iteráció a két xpontba exponenciálisan lecseng® szinusz függvény szerint konvergál. (n)
(n)
(n−1)
A ∆P1,2 = P1,2 − P1,2
(n)
vs P1,2 függvény megfelel a fázistérbeli ábrázolásnak (9. ábra).
Az η =0,355 görbe az iterációs térben oszcillálva konvergál, ennek megfelel®en a fázistérbeli trajektória egy egyenes, és látható, hogy a pontok a közepén s¶r¶södnek, a torlódási pont felel meg tehát a xpontnak a fázistérben. A másik két görbe ez esetben a kaotikus görbe. A különbség a kon gurációs térben ábrázolthoz képest, hogy a trajektóriák nem metszik egymást, a két hurok jól elkülönül egymástól. 9. ábra. A butadién P-iterációja a fázistérben ábrázolva η = 0, 358, 0, 359 és 0, 3595 értékek esetén a P1,2 mátrixelemre.
1.2 1 0.8
0,355
0,3595
dP(1,2)
0.6 0.4 0,359 0.2 0 -0.2 -0.4 -1.2
-1
-0.8
-0.6
-0.4
-0.2 P(1,2)
0
0.2
0.4
0.6
0.8
1
36 .
37 Vizsgáljuk meg az iteráció instabilitásának el®rejelezhet®ségét, a Ljapunov-exponenseket.
VII. táblázat. Az A mátrix sajátértékei és a Ljapunov-exponensek valós része a butadién molekulára
η
m1
λ1
m2
λ2
m3
λ3
0,100 0,35457 -1,03684 0,25908 -1,35062 0,16008 -1,83208 0,357 -0,72678 -0,31901 -0,87630 -0,13205 -0,97689 -0,02338 0,363 -0,74837 -0,28985 -0,90303 -0,10110 -1,00441 0,00440 0,386 -0,83296 -0,18277 -1,00391 0,00390 -1,10817 0,10271 0,431 -1,00104 0,00104 -1,19459 0,17780 -1,30316 0,26479
A VII. táblázatban használt m1 , m2 , m3 és λ1 , λ2 , λ3 rövidítések megegyeznek a Blochegyenletnél tárgyalt rövidítésekkel. A táblázatban két dolgot érdemes észrevenni: az egyik, hogy ebben az esetben, a héliumhoz hasonlóan, az m3 illetve a λ3 a leginkább megbízható, tehát ahhoz, hogy a legjobb jóslást adó Ljapunov-exponenst nyerjük, butadién esetén a teljes A mátrixot diagonalizálni kell Azonban ha a teljes mátrixot diagonalizáljuk, a spektrum degenerált lesz. Ha csak az alsó háromszöget, és a diagonálist vesszük gyelembe a diagonalizálásnál (m2 illetve λ2 ), akkor a spektrum nem degenerált, bár a predikció romlik. A másik felt¶n® változás a Bloch-egyenlet vizsgálatához képest, hogy az el®rejelezhet®ség romlik: míg a H2 esetében 0,001 pontossággal meg tudtuk mondani mely η értéknél romlik el a konvergencia, ebben az esetben a hiba az m3 esetében is 0,006. Ezzel együtt még mindig mondhatjuk, hogy az eljárás jó becslést ad az η kritikus értékére. 5.3.4. Benzol és más molekulák
A benzol P-iterációjának eredményeit vizsgálva (VIII. táblázat) a tapasztalat hasonló az el®z® példákban látottakhoz; az egyetlen szembet¶n® változás az, hogy semmilyen η értéknél sem válik kaotikussá az iteráció. Ez a benzol magas szimmetriájával magyarázható. Ismert, hogy a benzol molekula szimmetriája alapján a D6h pontcsoportba tartozik. Ez a magas szimmetria Hückel-szinten determinálja a megoldást. Az iterációra azért van csak szükség, hogy a P szimmetria-adaptált (6-os szimmetriájú) legyen, hiszen a Kekulé-mátrix nem az.
38 VIII. táblázat. Energia számítás benzolra a P-iteráció segítségével. η
E/β -egység
jelenség
0,100
8,000000
konvergál 3633 lépésben
0,250
8,000000
konvergál 162 lépésben
0,290
-
oszcilláció
0,291
-
divergál
0,292
-
divergál
A legjobb Ljapunov exponens az m3 értékb®l számítható, az így számolt spektrum azonban er®sen degenerált, ellentétben az m2 értékhez kapott spektrummal (IX. táblázat). Az iteráció xpontjának instabillá válásához tartozó η predikciójának hibája mindössze 0,003. IX. táblázat. Az A mátrix sajátértékei és a Ljapunov-exponensek valós része a benzol-molekulára η
m1
λ1
m2
λ2
m3
λ3
0,100 0,17410 -1,74813 0,11886 -2,12981 0,20000 -1,60944 0,290 -0,75363 -0,28285 -0,87390 0,13479 -0,98360 -0,01654 0,293 -0,76499 -0,26789 -0,88962 -0,11696 -1,0004 0,00040 0,315 -0,85295 -0,15905 -1,00345 0,00344 -1,12310 0,11609 0,350 -1,00480 0,00479 -1,17900 0,16466 -1,31000 0,27003
További számolási eredmények (X. táblázat) is meger®sítik azt a tapasztalatot, hogy a Piteráció esetén az m3 érték segítségével kapott kritikus η érték közelíti legjobban az iteráció instabillá válását. Ez alól egyedül a ciklooktatetraén kivétel, ez esetben az m2 érték jobb becslést ad. Összességében tehát azt mondhatjuk, hogy az ideális η becsléséhez a legjobb, ha a teljes
A mátrixot diagonalizáljuk, és az ekkor kapott m3 értékek alapján választunk η -t. Nagyobb molekulák esetén azonban (például nanocsövek) ez a számítás nagyon költséges lehet, mivel A négyindexes mátrix. Ezért jobb, ha az A mátrixnak csak az alsó háromszögét diagonalizáljuk (m1 ). Ha megelégszünk η egy-két század pontosságú becslésével, akkor ez a közelítés jó. Figyelembe véve, hogy a hiba pozitív irányú, ezt a tévedést korrekcióba tudjuk venni, ha a kapott η értéknél például 0,02-dal kisebb értéket választunk.
39 X. táblázat. Különböz® molekulák P-iterációjához tartozó tapasztalati (ηc ), illetve stabilitásvizsgálattal számolt kritikus η paraméterek. Az η alsó indexei a 16. oldalon ismertett módokra utal. molekula
ηc
η1
η2
η3
Butadién
0,357 0,431 0,386 0,363
Benzol
0,290 0,350 0,315 0,293
Ciklooktatetraén 0,310 0,330 0,314 0,293 Oktatetraén
0,310 0,331 0,324 0,312
5.3.5. Számítás lineáris poliénre
Láncok esetén lehet®ség van az ideális η érték empírikus meghatározására. Ha pár tíz szénatom hosszúságú láncra "trial and error " elven meghatározzuk az ideális η értéket, akkor tapasztalataink szerint ez az érték alkalmas lesz hosszabb láncok vizsgálatára is, ha a lánc egyéb paramétereit (pl. a kötéshosszt) nem változtatjuk. A számításokat egy 500 atomból álló lineáris poliénre végezzük, melyben minden C-C kötés 1,40 Å hosszúságú, tehát a lánc fémes. Megjegyezzük, hogy a legtöbb mások által korábban javasolt P mátrixot közvetlenül felépít® eljárás nem m¶ködik fémes láncok esetén. A XI. táblázatban gy¶jtöttük össze a számolás eredményeit. Az iterációt addig végeztük, amíg az energiaérték hetedik tizedesjegye már nem változott. Ezt az iterációs lépésszámot jelöltük N -nel. Az iteráció végén, a P mátrix hermiticitásának kis mérték¶ megsérülése miatt
szimmetrizálásra van szükség. A szimmetrizálás megsérti az idempotenciát, ezért szükség van a McWeeny által bevezetett ún. puri kációra [22]. A XI. táblázatban feltüntetett energiák a szimmetrizált és puri kált s¶r¶ségmátrixból számolt értékek. Az energiát kiszámoltuk a H mátrix közvetlen diagonalizálásával is. A kapott energia:
E = 635, 8940 . Ennek fényében tekintve a XI. táblázatra látható, hogy a kapott energia értékek még csak század pontosan sem egyeznek minden esetben. Ideális η értéknek ebben az esetben az η =0,15 tekinthet®, hiszen az csak 0,009-et téved az energiában, és körülbelül 100 lépés alatt konvergál.
40 XI. táblázat. Az 500 szénatomból álló fémes polién energiája az η paraméter különböz® értékeinél η
E/β -egység N
0,05
635,8933
755
0,10
635,8911
243
0,15
635,8849
106
0,20
635,8702
53
0,25
635,8917
1605
5.3.6. Számítás nanocs®re
A számítást egy úgynevezett (5,4) királis, 50 Å hosszúságú, 6,07 Å átmér®j¶, 368 atomból álló nanocs®re végeztük (10. ábra).
10. ábra. Az (5,4)-es, 50 Å hosszúságú nanocs®
A számolás eredményei analógok az eddig tapasztaltakkal: az iteráció kis η értékekre lassan konvergál, majd növelve a paraméter értékét egy optimum után divergál (XII. táblázat).
41 XII. táblázat. Az (5,4)-es királis, 50 Å hosszúságú nanocs® energiája P-iterációval számolva, az η paraméter különböz® értékeinél.
η
E/β -egység N
0,05 574,64611 878 0,10 547,64615 468 0,15 547,64616 308 0,20 574,64609 0,25
divergál
77 -
A H mátrix közvetlen diagonalizálásával is meghatároztuk a cs® energiáját:
E = 574, 64616 . Látható, hogy csövek esetén a P-iterációval kapott energia 0,001 pontossággal megegyezik a diagonalizálással kapott értékkel, ami nagyon jó egyezésnek mondható (ellentétben a fémes lineáris poliéneknél tapasztaltakkal (v.ö. XI. táblázat)). Nanocsövek esetén a P-iteráció jó közelítéssel megadja a Hamilton-mátrix diagonalizálásával kapott energiát. Ilyen nagy rendszerek esetén látszik a P-iteráció azon el®nye, hogy nem diagonalizációs eljárás: az iterációval gyorsabban lehetett energiaértékhez jutni.
42
6. ÖSSZEFOGLALÁS A dolgozatnak kett®s célja volt: egyrészt az iterációs káoszt tanulmányoztuk a Blochegyenlet és a P-iteráció példáján, mivel az iterációs paraméter (η ) bizonyos értékeinél ezek az iterációk kaotikusak. Másrészt azt vizsgáltuk, hogy ezek az iterációk milyen feltételek mellett alkalmazhatók energiaszámolásra. A káosz tanulmányozása során azt tapasztaltuk, hogy a kaotikus dinamikában ismertetett jelenségek analógjai az iteráció során is megjelennek. Az általunk vizsgált Blochtípusúegyenletek iteratív megoldásai kaotikus tulajdonságokat mutathatnak. Ez nem meglep®, hiszen az egyenletek alakja emlékeztet a kaotikus iterációk prototípusaként emlegetett logisztikus leképezére. Az iteráció eredetileg stabil xpontja az η paraméter változtatásával instabillá válik, megjelenik oszcilláció, a bifurkáció és a káosz. Az iterációs térben ábrázolt kaotikus iteráció rendezetlensége a kon gurációs- és a fázistérben elt¶nik, megjelennek a trajektóriák. A stabilitásvizsgálati analízis eredményeképp értelmezni tudtuk, hogy mely η értékeknél tapasztalunk konvergens iterációt. Ha ennél nagyobb abszolút érték¶ számot választunk η nak, akkor bifurkációkat, káoszt, vagy divergenciát tapasztalunk. Az iteráció jellegét indikáló Ljapunov-exponens egy négyindexes mennyiség (A) sajátértékéb®l származtatható. Az A felépítéséhez szükség van az iterált mátrix xpontbeli értékére. Az A konkrét alakját a Blochegyenlet és a P-iteráció kapcsán meghatároztuk. Az A elemei közötti redundanciákat kétféle módon küszöböltük ki, legmegbízhatóbbnak azonban azok az esetek bizonyultak, amelyekben a teljes A mátrixot diagonalizáljuk. A vizsgált példákban a Bloch-egyenlet esetén a xpont instabillá válását eredményez® kritikus η értéket 0,001 pontossággal, a P-iteráció esetén 0,005 pontossággal tudtuk megjósolni. A Bloch-egyenlet iterálása során azt tapasztaltuk, hogy ha nemcsak az alapállapot energiáját akarjuk számolni, hanem a gerjesztett állapotokét is, akkor az iteráció könnyen hamis xpontba jut, nem a megfelel® energiaértékeket kapjuk. Ezért a Bloch-egyenletnek ezt a formáját nem javasoljuk gerjesztett állapotok energiájának számolásra. Mivel a P-iteráció esetén elkerüljük a mátrixok diagonalizálását, olcsón tudunk energiát számolni viszonylag nagy molekulákra. Ezt sok száz atomot tartalmazó konjugált kett®kötés¶ molekulák energiaszámolásával mutattuk be.
43
[1] Tél Tamás és Gruiz Márton, Kaotikus Dinamika (Nemzeti Tankönyvkiadó, Budapest, 2002). [2] E. Ott, Chaos in Dynamical Systems (Cambridge University Press, Cambridge, 1993). [3] D. K®halmi and Á. Szabados and P. R. Surján, Phys. Rev. Letters 95, 013002 (2005). [4] I. Mayer, Simple Theorems, Proofs, and Derivations in Quantum Chemistry (Kluwer Acedenic/Plenum Publisher, New York, 2003). [5] Per-Olov Löwdin, Int J Quant Chem 72, 379 (1999). [6] V. Kvasni£ka, Chech J Phys B24, 605 (1974). [7] I. Lindgren and J. Morrison, Atomic Many-Body Theory (Springer, Berlin, 1986). [8] J. H. Wilkinson and C. Reinsch, Handbook for Automatic Computation (Springer, Berlin, 1971), Vol. 2. [9] R. A. Poirier and M. Peterson, Program MUNGAUSS, Dept.Chemistry, Memorial Univ.St.Johns, Canada (1989). [10] P. R. Surján, Program BP-MUNGAUSS, Dept.Theoretical Chemistry, Eötvös University, Budapest (2002). [11] Péter R. Surján, Second Quantized Approach to Quantum Chemistry (Springer-Verlag Heidelberg, Berlin, 1989). [12] I. A. Abrikosov, A. M. N. Niklasson, S. I. Simak, B. Johansson, A. V. Ruban, and H. L. Skriver, Phys. Rev. Letters 76, 4203 (1996). [13] R. Reidinger and M. Benard, J. Chem. Phys. 94, 1222 (1991). [14] R. McWeeny, Phys. Rev. 126, 1028 (1962). [15] X. P. Li, R. W. Nunes, and D. Vanderbilt, Phys. Rev. B 47, 10891 (1993). [16] M. Challacombe, J. Chem. Phys. 110, 2332 (1999). [17] A. H. R. Palser and D. E. Manolopoulos, Phys. Rev. B 58, 12704 (1998). [18] W. Yang, Phys. Rev. Letters 66, 1438 (1991). [19] W. Yang, Phys. Rev. A 44, 7823 (1991). [20] W. Yang and T.-S. Lee, J. Chem. Phys. 103, 5674 (1995). [21] K. Nemeth and G. E. Scuseria, J. Chem. Phys. 113, 6035 (2000).
44 [22] R.McWeeny, Rev.Mod.Phys 32, 335 (1960). [23] R. McWeeny, Phys. Rev. 126, 1028 (1961).
