\documentclass[12pt,a4paper,oneside]{article}

\usepackage[slovak]{babel}
\usepackage[T1]{fontenc}
\usepackage[latin2]{inputenc}

\usepackage{color}
\usepackage[colorlinks,citecolor={blue},hyperindex]{hyperref}

\newcommand{\BBBB}{\begin{equation}}
\newcommand{\EEEE}{\end{equation}}
\newcommand{\BBAA}{\begin{eqnarray}}
\newcommand{\EEAA}{\end{eqnarray}}

\title{Návody na cvičenia z predmetu \\ \textbf{Počítačová fyzika \\ }
(letný semester 2006)}
\author{}
\date{\today}

\begin{document}

\maketitle

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\label{sec:uvod}Úvod}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pohybovou rovnicou sveta atómov, molekúl, fotónov atď. je 
\textit{Schr\"{o}dingerova rovnica} \cite{qm}
\BBBB
i \hbar \frac{\partial \Psi}{\partial t} = \hat{H} \Psi \ .
\label{SchrR}
\EEEE
Špeciálnym a veľmi dôležitým prípadom jej riešení sú stacionárne riešenia. 
Aby sme ich našli,
uvažujme rovnicu pre vlastné funkcie a vlastné hodnoty hamiltoniánu:
\BBBB
\hat{H} \Phi = \lambda \Phi \ .
\label{bezcasova}
\EEEE
$\Phi$ je vlastná funkcia (jedna z mnohých), ktorá nezávisí od času.
$\lambda$ je príslušná vlastná hodnota.
Zkonštruujme funkciu
\BBBB
\Psi(t) = e^{-\frac{i}{\hbar} \lambda t} \Phi \ .
\label{stacries}
\EEEE
Dosadením tejto funkcie do Schr\"{o}dingerovej rovnice zistíme, že táto funkcia
je jej riešením. Teda rovnica (\ref{bezcasova}), zatiaľ chápaná len ako 
nejaká matematická rovnica pre vlastnú sústavu (funkcií a hodnôt), 
má výnimočné postavenie v~tom, že jej vlastné funkcie po~vynásobení faktorom
$e^{-\frac{i}{\hbar} \lambda t}$ sú riešeniami Schr\"{o}dingerovej rovnice
(\ref{SchrR}). Dokonca tým dosadením zisťujeme, že pre~takýto typ riešení 
sa Schr\"{o}dingerova rovnica zjednodušuje na tvar rovnice (\ref{bezcasova}).
Z~tohto dôvodu a aj preto, že (\ref{bezcasova}) neobsahuje čas, nazývame
rovnicu (\ref{bezcasova}) \textit{bezčasová Schr\"{o}dingerova rovnica}.
Riešenia tvaru (\ref{stacries}) nazývame \textit{stacionárne stavy}. Na~ich
určenie teda stačí riešiť bezčasovú Schr\"{o}dingerovu rovnicu. To bude
náplňou týchto cvičení.

Ešte k vlastnej hodnote $\lambda$: jej fyzikálny význam spočíva v~tom, že
to je vlastná energia
hamiltoniánu prislúchajúca vlastnému stavu $\Phi$. Teda
\BBBB
\lambda = E = \int \Phi^* \hat{H} \Phi d\tau \ ,
\EEEE
kde integrovanie ide cez~celý priestor alebo tú jeho časť, kde sa popisované
častice môžu vyskytnúť. $d\tau$ je stručné označenie objemového elementu;
v~prípade jednej častice je $d\tau = dx dy dz = d^3 r$.


Bezčasová Schr\"{o}dingerova rovnica má nekonečne veľa riešení. My budeme
hľadať asi len to, ktoré prislúcha najnižšej energii - \textit{základný
stav}. Ten sa v~prírode aj v experimentoch vyskytuje najčastejšie, takže
je najdôležitejší (a ráta sa \textit{pomerne} jednoducho).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\label{sec:ries}Riešenie bezčasovej Schr\"{o}dingerovej rovnice}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Hľadanú neznámu funkciu, povedzme že závislú od priestorových premenných
$\vec{r}_1, \vec{r}_2, ..., \vec{r}_N$, stručne značených $\bar{r}$, zapíšeme
ako lineárnu kombináciu nejakých známych (vo všobecnosti komplexných) 
funkcií $f_1(\bar{r}), f_2(\bar{r}), ..., f_M(\bar{r})$:
\BBBB
\Phi(\bar{r}) = \sum_{j=1}^M c_j f_j(\bar{r}) \ ,
\label{linkomb}
\EEEE
kde $c_j$ sú neznáme komplexné koeficienty, ktoré treba určiť. 
Súbor funkcií $f_j(\bar{r})$ nazývame \textit{báza}.
Dosadením navrhnutého 
tvaru riešenia do (\ref{bezcasova}) postupne dostávame
$$
\sum_{j=1}^M c_j \hat{H} f_j(\bar{r}) = \sum_{j=1}^M c_j E f_j(\bar{r}) \ ;
$$
po vynásobení zľava funkciou $f_i^*(\bar{r})$ a po preintegrovaní cez všetky
premenné nám vyjde
\BBBB
\sum_{j=1}^M c_j \underbrace{\int f_i^*(\bar{r}) \hat{H} f_j(\bar{r}) 
d\tau}_{H_{ij}} = 
\sum_{j=1}^M E\, c_j \underbrace{\int f_i^*(\bar{r}) f_j(\bar{r}) 
d\tau}_{S_{ij}} \ .
\label{matice}
\EEEE
Zaviedli sme vhodné označenia pre integrály. Vidno, že $H_{ij}$ možno chápať
ako prvky nejakej $M \times M$ matice, podobne aj $S_{ij}$. Maticu $S$ nazývame
\textit{prekryvovou maticou}, kedže vyjadruje, ako sa bázové funkcie 
\glqq prekrývajú\grqq~ v~priestore. Teda bezčasovú Schr\"{o}dingerovu
rovnicu (ktorá je diferenciálnou rovnicou, lebo v hamiltoniáne vystupujú
derivácie podľa priestorových premenných) sme transformovali na sústavu
algebraických rovníc
\BBBB
\sum_{j=1}^M H_{ij} c_j = \sum_{j=1}^M E S_{ij} c_j \ \ ,\forall i=1,...,M \ ,
\EEEE
čiže
%\BBBB
$$
\left(
\begin{array}{cccc}
H_{11} & H_{12} & ... & H_{1M} \\
H_{21} & H_{22} & ... & H_{2M} \\
  .    &  .     & ... &  .     \\
  .    &  .     & ... &  .     \\
H_{M1} & H_{M2} & ... & H_{MM}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
...   \\
c_{M}
\end{array}
\right) =
E
\left(
\begin{array}{cccc}
S_{11} & S_{12} & ... & S_{1M} \\
S_{21} & S_{22} & ... & S_{2M} \\
  .    &  .     & ... &  .     \\
  .    &  .     & ... &  .     \\
S_{M1} & S_{M2} & ... & S_{MM}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
 .    \\
 .    \\
c_{M}
\end{array}
\right)
$$
\BBBB
\label{zovseob}
\EEEE
Je to tzv. zovšeobecnený problém vlastných vektorov ($c$) a hodnôt ($E$).
V~prípade, že by prekryvová matica $S$ bola jednotková, prešla by rovnica
(\ref{zovseob}) na~štandardnú úlohu hľadania vlastných vektorov a hodnôt.
V~každom prípade však matice $H$ a $S$ sú hermitovské, teda $H_{ji}^* =
H_{ij}$ a podobne pre~$S$.


Tvar Hamiltonovej matice $H$ aj prekryvovej matice $S$
samozrejme závisí (pre~daný hamiltonián)
len a len od~voľby bázy. Aké bázové
funkcie je vhodné voliť? V~každom prípade také, aby sa hľadaná vlastná 
funkcia dala v~tejto báze [viď (\ref{linkomb})] dostatočne presne
vyjadriť. Zrejme bude treba buď veľa funkcií v~báze (vysoké $M$) alebo 
veľmi šikovne navrhnutý tvar bázových funkcií (jedna z nich veľmi podobná
hľadanému riešeniu; potom už ďalších netreba veľa).
Ak chceme, aby problém bol výpočtovo nenáročný, šikovne si teda zvolíme
bázové funkcie a máme len $M$ neznámych konštánt $c_j$. Nesmierne zložitý
počiatočný  problém (\ref{bezcasova}), kde je neznámych nesmierne veľa
(sú to hodnoty vlastnej vlnovej funkcie $\Phi(\bar{r})$ vo všetkých bodoch
priestoru alebo aspoň vo~veľkom konečnom počte uzlov priestoru nasekaného na
veľmi malé kúsky), sme teda voľbou šikovnej bázy zredukovali na~problém, kde 
je neznámych len $M$, čo môže byť celkom malé číslo (povedzme jednotky až
desiatky) 
\footnote{Nie vždy však tvar bázových funkcií volíme tak, aby ich
bolo čo najmenej. To totiž prináša istý kompromis s presnosťou, aj keď často
nie zlý kompromis. Mnohokrát sa však hodí používať radšej \glqq menej 
šikovnú\grqq, teda väčšiu bázu, ktorá je na oplátku spoľahlivejšia. 
Takouto bázou
je (v~jednom rozmere) napr. súbor funkcií $1, \cos x , \sin x , \cos 2x,
\sin 2x, ... $.}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\label{sec:diag}Zovšeobecnený diagonalizačný problém}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ideme riešiť maticovú rovnicu (\ref{zovseob}). Na pomoc si zoberieme 
Numerické recepty \cite{numrec}. Hneď na začiatok poznamenajme, že 
\textbf{budeme
predpokladať takú bázu, v ktorej sú všetky tie maticové prvky $H_{ij}$ a
$S_{ij}$ reálne}, čiže tieto matice sú \textit{symetrické} (čo je špeciálny 
prípad hermitovosti). Reálne bázové funkcie v prípade atómu neznamenajú 
žiadne obmedzenie všeobecnosti alebo presnosti riešenia.


Riešenie má niekoľko krokov:
\begin{itemize}
\item[1.)] Choleského rozklad matice $S$ na~súčin dolnej a hornej 
trojuholníkovej matice
\item[2.)] transformácia pôvodného zovšeobecneného diagonalizačného problému 
na~zvyčajný diagonalizačný problém
\item[3.)] riešenie tohto problému niektorou zo~štandardných metód, napr. 
Jacobiho metódou
\end{itemize}
Takže poďme to rozobrať podrobnejšie.
%==============================================================================
\subsection{\label{ssec:cholesky}Choleského rozklad na súčin trojuholníkových 
matíc}
%==============================================================================
S hocijakou štvorcovou maticou $S$ sa dá spraviť tzv. LU rozklad 
\cite{numrec} na~súčin poddiagonálnej
matice $L$ a naddiagonálnej matice $U$, pričom na~ich diagonálach sú 
vo~šeobecnosti nenulové hodnoty, L (od~slova Lower) má všetky prvky 
nad~diagonálou nulové a U (od~Upper) zasa samé nuly pod~diagonálou.  Teda
$S = L U$. V~prípade, že matica $S$ je navyše symetrická a kladne 
definitná, matica $U$ sa dá nájsť ako $U = L^T$, čo dáva Choleského rozklad:
\BBBB
S = L L^T \ .
\label{cholesky}
\EEEE
Prvky matice $L$ sa rátajú nasledovne:
\BBAA
L_{ii} = \left(S_{ii} - \sum_{k=1}^{i-1} L_{ik}^2 \right)^{1/2} \ , \\
L_{ji} = \frac{1}{L_{ii}} \left(S_{ij} - \sum_{k=1}^{i-1} L_{ik} L_{jk}\right)
\ , \ \ j = i+1, i+2, ..., M \ .
\EEAA
(Ak vychádza horná hranica sumy nižšia ako dolná, t.j. pre~$i=1$, treba
sumu považovať za nulovú.)
Naddiagonálne prvky matice $L$ sú podľa definície nulové.
Nakoniec, aby sme elementy $L$ vedeli zrátať, musíme ich, ako 
sa dá zistiť, počítať v~správnom 
poradí; najprv $L_{11} = \sqrt{S_{11}}$, potom $L_{21}$, potom $L_{31}$, atď.
%==============================================================================
\subsection{Transformácia zovšeobecneného problému na~zvyčajný diagonalizačný
problém}
%==============================================================================
Zopakujeme, že máme riešiť problém (\ref{zovseob}), čiže
\BBBB
H c = E S c \ .
\EEEE
Veľmi jednoducho by sme ho mohli transformovať na~zvyčajný diagonalizačný 
problém vynásobením
zľava maticou inverznou k $S$:
$$
(S^{-1} H) c = E c \ .
$$
Teda by stačilo nájsť vlastnú sústavu matice $S^{-1} H$. Takto to však 
nebudeme robiť. Je to trochu náročnejšie, lebo matica $S^{-1} H$ už nie je
symetrická.
Avšak v prípade, že aspoň matice $H$ a $S$ sú symetrické a matica $S$ 
ešte navyše kladne definitná, čo v~našom prípade naozaj všetko platí, 
je možné úlohu transformovať na~problém
diagonalizácie symetrickej matice
s~využitím Choleského rozkladu $S = L L^T$.
Dostávame teda \cite{numrec}
$$
H c = E L L^T c \ .
$$
Vynásobením zľava maticou $L^{-1}$ máme
$$
L^{-1} H c = E L^T c \ .
$$
Medzi $H$ a $c$ vložíme jednotkovú maticu zapísanú ako $I = (L^T)^{-1} L^T$:
$$
L^{-1} H (L^T)^{-1} L^T c = E L^T c \ .
$$
Označíme
\BBBB
F \equiv L^{-1} H (L^T)^{-1}
\EEEE
a máme
\BBBB
F (L^T c) = E  (L^T c) \ .
\EEEE
Matica $F$ je symetrická. (Skúsiť dokázať.) 
\footnote{Pre každú nesingulárnu štvorcovú maticu platí $(A^{-1})^T = 
(A^T)^{-1}$. To sa dá ukázať dosť ľahko a potom aj to, že matica $F$ je
symetrická.}
Takže máme
už jednoduchší diagonalizačný problém pre~tie isté vlastné hodnoty $E$ 
ako bol ten pôvodný problém a pre~vlastné vektory $L^T . c$.
Ešte ako efektívne nájsť maticu $F$? V~dvoch krokoch:\\
a) ozn. $H (L^T)^{-1} \equiv Y \Longrightarrow L Y^T = H$, čo je 
sústava lineárnych algebraických rovníc pre~neznámu $Y^T$. Túto sústavu
vyriešime ľahko, lebo $L$ je trojuholníková matica.\\
b) Máme teda zatiaľ $F = L^{-1} Y \Longrightarrow L F = Y$, čo je opäť
ten istý typ problému, z ktorého $F$ už ľahko nájdeme.
%==============================================================================
\subsection{\label{ssec:jacobi}Jacobiho metóda diagonalizácie reálnej 
symetrickej matice}
%==============================================================================
Je v numerických receptoch \cite{numrec}, nebudeme to tu podrobnejšie 
rozpisovať.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\label{sec:vodik}Aplikácia na atóm vodíka}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Aj keď môže pripadať dosť trápne riešiť atóm vodíka numerickou metódou, ako 
prvé zoznámenie sa s praktickými metódami elektrónovej štruktúry v~ich 
najjednoduchšej podobe je to užitočné.


Hamiltonov operátor
\BBBB
\hat{H} = -\frac{\hbar^2}{2m_e} \vec{\nabla}^2 - \frac{1}{4\pi\varepsilon_0}
\frac{e^2}{r^2} \ .
\label{vodik-ham}
\EEEE
Úlohou je nájsť stacionárne riešenia bezčasovej Schr\"{o}dingerovej rovnice
(\ref{bezcasova}). Hneď aj uveďme, že presné analytické riešenia sú \cite{qm}
\BBBB
\Phi_{nlm}(\vec{r}) = R_{nl}(r,\vartheta) Y_{lm}(\vartheta,\varphi) 
\label{vodik-anries}
\EEEE
kde $R_{nl}(r)$ sú radiálne funkcie a $Y_{lm}(\vartheta,\varphi)$ sú guľové 
funkcie.
Celé čísla $n = 1, 2, 3, ...$, $l = 0, 1, ..., n-1$ a $|m| \leq l$.
Uvedieme len najnižšie funkcie.
\BBBB
R_{10}(r) = \left(\frac{1}{a_1}\right)^{3/2} e^{-r/a_1}\ ,
\label{radialna}
\EEEE
kde $a_1 = \hbar^2 / (m {e'}^2)$ je Bohrov polomer atómu vodíka.
\BBBB
Y_{00}(\vartheta,\varphi) = \frac{1}{\sqrt{4\pi}} \ .
\label{gulova}
\EEEE
Vlastné energie sú
\BBBB
E_n = - \frac{m_e {e'}^4}{2 \hbar^2 n^2}\ , \ \ \ n = 1, 2, 3, ...
\label{vodik-energie}
\EEEE
Značenie $e' = e^2 / (4 \pi \varepsilon_0)$.


V programe budeme používať \textit{atómové jednotky}, 
v ktorých jednotkou náboja je
náboj protónu (teda $e = 1$), jednotkou hmotnosti je hmotnosť elektrónu
(teda $m_e = 1$), jednotkou dĺžky je Bohrov polomer $a_1$, 
kladie sa $\hbar = 1$
a $\varepsilon_0 = 1/(4\pi)$. Potom vychádza, že základný stav atómu vodíka
má energiu $-1/2$. Atómová jednotka energie sa nazýva aj Hartree a je to 
energia $-m_e {e'}^4 / \hbar^2$.

%==============================================================================
\subsection{\label{ssec:baza}Báza}
%==============================================================================
Ako bázu použijeme funkcie Gaussovho typu.
Gaussova funkcia je
\BBBB
g(x) = e^{-x^2} \ .
\EEEE
Funkciami Gaussovho typu (skrátane gaussiánmi) nazveme funkcie 
tvaru \cite{carsky_urban}
(pre jednoduchosť ich uvádzame centrované v počiatku súradnicovej sústavy)
\BBBB
g(\vec{r}) = {\cal N} x^l y^m z^n e^{-\zeta r^2} \ ,
\label{gaussian}
\EEEE
kde ${\cal N}$ je nateraz nepodstatná normovacia konštanta.
$l, m, n$ sú nezáporné celé čísla a $\zeta$ je kladná 
konštanta. Dá sa ukázať, že hocijaká funkcia premenných $x, y, z$ sa dá 
vyjadriť ako
lineárna kombinácia gaussiánov. Teda gaussiány tvoria úplnú bázu; dokonca
\glqq nadúplnú\grqq~ (teda nejaká funkcia sa dá vyjadriť 
ako ich lineárna kombinácia 
nie jedným spôsobom ale mnohými), čo je trochu vidno aj z~toho, že nie sú 
ortogonálne.


Pre naše účely nateraz vystačíme len s maličkou podmnožinou všetkých 
gaussiánov. Pre jednoduchosť sa totiž najprv budeme zaujímat len o~základný
stav atómu vodíka. Ten má guľovo symetrickú vlnovú funkciu. Viď 
(\ref{vodik-anries}).
Preto vystačíme s~gaussiánmi závislými len od~veľkosti $\vec{r}$. Budeme ich
značiť
\BBBB
g_{n \alpha}(r) = {\cal N}_{n \alpha} r^{n-1} e^{-\zeta_{\alpha} r^2} \ .
\label{baza}
\EEEE
$n = 1, 2, 1, ..., n_{\rm max}$ nám má pripomínať hlavné kvantové číslo;
preto je tam nárokom mocnina $n-1$, aby číslovanie šlo od jednotky, tak ako 
hlavné kvantové číslo.
Koeficienty $\zeta_{\alpha}$ sú kladné čísla,
ktoré si musíme zvoliť nejak skusmo. Pre~poriadok ich indexujeme indexom 
$\alpha = 1, 2, ..., \alpha_{\rm max}$. Normovacie konštanty 
${\cal N}_{n \alpha}$ budeme vyberať
tak, aby 
\BBBB
\int g^2_{n \alpha}(r) d^3r = 1 \ ,
\label{normovanie}
\EEEE
aj keď \textbf{normovanie na~jednotku nie je nevyhnutné} - viď časť 
\ref{sec:ries}, 
z~ktorej vidno, že bázové funkcie síce musia byť normovateľné, ale 
normovať na~jednotku ich naozaj nemusíme a v~programe môžme pre~jednoduchosť
mať konštanty ${\cal N}$ rovné napr. 1.
Pre~konštrukciu (vyčíslenie) matice $S$ a matice $H$ 
musíme mať nielen bázové funkcie, ale
aj zadefinovať si ich poradie. To môže byť ľubovoľné, ale samozrejme je 
vhodné to poradie zadefinovať nejak prakticky. 
Najprv nám pôjde len o~jednoduché naprogramovanie, takže môžme bázové funkcie
voliť a očíslovať napr. takto:
\BBBB
g_{1 1}, g_{1 2}, \dots, g_{1 \alpha_{\rm max}}, g_{2 1}, \dots, 
g_{2 \alpha_{\rm max}}, \dots \dots, g_{n_{\rm max} \alpha_{\rm max}}.
\label{poradie}
\EEEE
Takáto báza má teda $n_{\rm max} \alpha_{\rm max} \equiv M$ funkcií, čiže
matice budú (štvorcové) o~rozmeroch $M \times M$.
Pre~stručnosť zápisu si zadefinujeme aj združený index 
\BBBB
i\equiv (m,\alpha)
\label{zdruzeny}
\EEEE
[alebo aj v znační $j \equiv (n,\beta)$] 
a potom bázové funkcie sú jednoducho $g_1$, $g_2$, \dots, $g_M$. 
Ešte raz, aby neboli pochybnosti: vzťah medzi $(m,\alpha)$ a indexom $i$ je 
taký, že $i$ vyjadruje poradie bázovej funkcie v rade (\ref{poradie}).
V~programe ten vzťah získame najjednoduchšie takto:
\begin{verbatim}
i = 0;
for (m=1; m<=nmax; m++) for (al=1; al<=alfamax; al++) {
        i++;
}
\end{verbatim}
A ak potrebujeme dvojitý cyklus (teda cez~$i$ aj $j$), tak
\begin{verbatim}
i = 0;
for (m=1; m<=nmax; m++) for (al=1; al<=alfamax; al++) {
        i++;
        j = 0;
        for (n=1; n<=nmax; n++) for (be=1; be<=alfamax; be++) {
                j++;
        }
}
\end{verbatim}
Združený index sa nám 
zíde najmä pri~zápise maticových prvkov, kde by sme inak museli písať napr.
$S_{m \alpha, n \beta}$, čo teraz stručne zapíšeme $S_{i j}$.
%==============================================================================
\subsection{\label{ssec:matS}Výpočet elementov prekryvovej matice}
%==============================================================================
Podľa (\ref{matice}) a (\ref{baza})
\BBBB
S_{m \alpha, n \beta} \equiv S_{i j} = \int g_{m \alpha}(r) g_{n \beta}(r) 
d^3 r \ .
\EEEE
Dosadením a transformovaním do~sférických súradníc vychádza
\begin{large}
\begin{center}
\framebox{$
S_{i j} = 4 \pi {\cal N}_{m \alpha} {\cal N}_{n \beta}
I_{m+n}(\zeta_{\alpha}+\zeta_{\beta}) \ ,
$}
\end{center}
\end{large}
\vskip-8mm
\BBBB
\label{matS}
\EEEE
\vskip+4mm
kde sme zaviedli označenie\\
\begin{large}
\begin{center}
\framebox{ $I_n(a) \equiv \int_0^\infty x^n e^{-a x^2} dx \ .$ }
\end{center}
\end{large}
\vskip-8mm
\BBBB
\label{integ_n}
\EEEE
\vskip+4mm
Tento integrál vieme po~použití per-partes rekurzívne zrátať, pričom vychádza
\begin{large}
\begin{center}
\framebox{ $I_n = \frac{n-1}{2a} I_{n-2} \ .$ }
\end{center}
\end{large}
\vskip-8mm
\BBBB
\label{integ_rekur}
\EEEE
\vskip+4mm
Potrebujeme k~tomu ešte prvé dva integrály. Tie sú
$$
I_0 = \frac{1}{2} \sqrt{\frac{\pi}{a}}
$$
a
$$
I_1 = \frac{1}{2a} \ .
$$
Potom
$$
I_{2k-1} = \frac{(2k-2)!!}{2^k a^k}
$$
a
$$
I_{2k} = \frac{\sqrt{\pi}(2k-1)!!}{2^{k+1} a^{k + 1/2}} \ ,
$$
pre $k = 1,2,3, ...$.

\noindent \underline{Poznámka} (praktická): V~programe nepoužívajte 
pre~výpočet integrálov $I_n(a)$ tieto explicitné vzťahy; omnoho efektívnejšie
je zo~známych hodnôt integrálov $I_0(a)$ a $I_1(a)$ určiť tie vyššie priamo
pomocou rekurentného vzťahu (\ref{integ_rekur}).


Ďalej môžme ešte určiť normovacie konštanty ${\cal N}_{n \alpha}$ 
napr. z~podmienky (\ref{normovanie}). Vychádzajú
\BBBB
{\cal N}_{n \alpha} = \frac{2^{2n-1/2} \zeta_{\alpha}^{n+1/2}}
{\pi^{3/2} (2n-1)!!} \ .
\EEEE
\noindent \underline{Poznámka} (praktická): V~programe radšej použijeme
inú definíciu ľubovoľných konštánt ${\cal N}_{n \alpha}$: nech sú všetky
rovné 1, aby sme mali menej práce. Je však možné, že použitím podmienky
(\ref{normovanie}) by sme na~oplátku dosiahli numericky stabilnejší kód.
%==============================================================================
\subsection{\label{ssec:matH}Výpočet maticových elementov hamiltoniánu}
%==============================================================================
Podľa (\ref{matice}) a (\ref{baza}) [viď aj značenie (\ref{zdruzeny})
združených indexov]
\BBBB
H_{m \alpha, n \beta} \equiv H_{i j} = \int g_{m \alpha}(r) \hat{H} 
g_{n \beta}(r) d^3 r \ ,
\EEEE
kde hamiltonián je daný vzťahom (\ref{vodik-ham}). 
Ešte si ho rozdelíme na operátory kinetickej a potenciálnej energie:
\BBBB
\hat{H} = \hat{T} + \hat{V}
\EEEE
a máme teda aj pre maticové elementy rozdelenie
\BBBB
H_{ij} = T_{ij} + V_{ij} \ .
\EEEE
Najprv vybavíme potenciálnu energiu. Tá sa ráta úplne podobne ako S.
Vychádza
\begin{large}
\begin{center}
\framebox{ $V_{ij} = -\frac{1}{4\pi\varepsilon_0} {\cal N}_{m \alpha} 
{\cal N}_{n \beta}4 \pi I_{m+n-1}(\zeta_{\alpha} + \zeta_{\beta}) \ .$ }
\end{center}
\end{large}
\vskip-10mm
\BBBB
\label{matV}
\EEEE
\vskip+6mm
\textbf{Pozor na konštantu:} v~programe budeme pracovať v~atómových 
jednotkách a v~nich je $1/(4\pi\varepsilon_0) = 1$.


Trochu ťažšie je to s maticou kinetickej energie. Ale ak si spomenieme 
na~vyjadrenie Laplaceovho operátora vo~sférických súradniciach \cite{qm}, 
ide aj toto
pomerne jednoducho. Laplaceov operátor je
\BBBB
\vec{\nabla}^2 = \frac{1}{r^2} \frac{\partial}{\partial r} 
\left(r^2 \frac{\partial}{\partial r}\right) + 
\frac{1}{r^2} \nabla^2_{\vartheta,\varphi} \ ,
\label{Laplace}
\EEEE
kde výlučne uhlovo závislý operátor
\BBBB
\nabla^2_{\vartheta,\varphi} = \frac{1}{\sin \vartheta} 
\frac{\partial}{\partial \vartheta} \left(\sin \vartheta 
\frac{\partial}{\partial \vartheta}\right) +
\frac{1}{\sin^2 \vartheta} \frac{\partial^2}
{\partial \vartheta^2} \ .
\EEEE
Kedže pracujme len s~radiálne závislými funkciami, druhý člen operátora
(\ref{Laplace}) po~derivovaní podľa uhlových premenných vyrobí nulu. Ostáva
teda len príspevok z~prvého člena. Ten vychádza
$$
T_{ij} = -\frac{\hbar^2}{2 m_e} {\cal N}_{m \alpha} {\cal N}_{n \beta}
\int r^{m-1} e^{-\zeta_\alpha r^2} \vec{\nabla}^2 \left(r^{n-1} 
e^{-\zeta_\beta r^2} \right) d^3 r = 
$$
$$
= -\frac{\hbar^2}{2 m_e} {\cal N}_{m \alpha} {\cal N}_{n \beta}
\int r^{m-1} e^{-\zeta_\alpha r^2}
\frac{1}{r^2} \frac{\partial}{\partial r}
\left[ r^2 \frac{\partial}{\partial r} 
\left(r^{n-1} e^{-\zeta_\beta r^2}\right) \right] d^3 r =
$$
(teraz per-partes)
$$
= \frac{\hbar^2}{2 m_e} {\cal N}_{m \alpha} {\cal N}_{n \beta} 4\pi
\int_0^\infty dr \frac{\partial}{\partial r}
\left(r^{m-1} e^{-\zeta_\alpha r^2}\right) r^2
\left(r^{n-1} e^{-\zeta_\beta r^2}\right) \ .
$$
Derivácie roznásobíme (trochu prácne ale nie zložité) a dostaneme 
členy s~už známymi
integrálmi typu (\ref{integ_n}). Výsledok bude
\begin{large}
\begin{center}
\framebox{$
T_{ij} = \frac{\hbar^2}{2 m_e} {\cal N}_{m \alpha} {\cal N}_{n \beta} 4 \pi {\cal I} \ ,
$}
\end{center}
\end{large}
\vskip-10mm
\BBBB
\EEEE
\vskip+6mm
kde
$$
{\cal I} = (m-1) (n-1) I_{m+n-2}(\zeta_\alpha+\zeta_\beta)
- 2 (n-1) \zeta_\alpha I_{m+n}(\zeta_\alpha+\zeta_\beta) +
$$
\BBBB
- 2 (m-1) \zeta_\beta  I_{m+n}(\zeta_\alpha+\zeta_\beta)
+ 4 \zeta_\alpha \zeta_\beta I_{m+n+2}(\zeta_\alpha+\zeta_\beta) \ .
\EEEE
(V~programe dáme $\hbar = 1$ a aj $m_e = 1$).
Funkcie ako napr. $I_{m+n}(\zeta_\alpha+\zeta_\beta)$ sú opäť len integrály
typu (\ref{integ_n}).
Integrály $I_{m+n}$ a $I_{m+n+2}$ vieme použitím (\ref{integ_rekur}) vyjadriť
pomocou $I_{m+n-2}$. Potom dostávame ešte jednoduchšie vyjadrenie
\begin{large}
\begin{center}
\begin{tabular}{|c|}
\hline
${\cal I} = \left\{(m-1)(n-1) - \left[(n-1)\zeta_\alpha + 
(m-1)\zeta_\beta\right] \frac{m+n-1}{\zeta_\alpha + \zeta_\beta}\right. + $ \\
$+ \left. (m+n+1)(m+n-1) \frac{\zeta_\alpha \zeta_\beta}{(\zeta_\alpha + 
\zeta_\beta)^2}\right\} I_{m+n-2}(\zeta_\alpha+\zeta_\beta) \ . $ \\
\hline
\end{tabular}
\end{center}
\end{large}
\vskip-3mm
\BBBB
\EEEE
\vskip+1mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\label{sec:molekula}Aplikácia na molekulu s jedným elektrónom}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Tu už ide o viaccentrový problém. 
Teória z častí \ref{sec:uvod}, \ref{sec:ries} a \ref{sec:diag} platí pre 
jednoelektrónový
problém všeobecne, teda aj pre (jeden) elektrón v hocijakom vonkajšom 
poli - napr.  aj v elektrostatickom poli viacerých atómových jadier.


Len stručne napíšeme výsledné vzťahy
pre maticové elementy z jednoduchých gaussiánov, ktoré
sme si odvodili na cvičení.
Gaussove funkcie teraz budeme značiť
\BBAA
g_A(\vec{r}) &= \exp\left[-\alpha (\vec{r}-\vec{A})^2\right]
\nonumber
\\
g_B(\vec{r}) &= \exp\left[-\beta (\vec{r}-\vec{B})^2\right] \ ,
\EEAA
kde $g_A(\vec{r})$ je bázová funkcia lokalizovaná okolo bodu $\vec{A}$
a obdobne $g_B(\vec{r})$ okolo bodu $\vec{B}$.
Platí, že súčin dvoch gaussiánov je opäť gaussián:
\begin{large}
\BBBB
\begin{array}{|c|}
\hline \\
\exp\left[-\alpha (\vec{r}-\vec{A})^2\right]
\exp\left[-\beta (\vec{r}-\vec{B})^2\right] = \\ \\
\exp\left[-\frac{\alpha \beta}{\alpha+\beta} (\vec{A}-\vec{B})^2\right]
\exp\left[-(\alpha+\beta) (\vec{r}-\vec{P})^2\right]
\\ \\
\hline
\end{array}
\label{sucin}
\EEEE
\end{large}
kde $\vec{P}$ je centrum súčinového gaussiánu,
\BBBB
\vec{P} = \frac{\alpha \vec{A} + \beta \vec{B}}{\alpha+\beta} \ .
\label{sucinove_centrum}
\EEEE
Prekryvová matica:
\begin{large}
\BBBB
\begin{array}{|c|}
\hline\\
S_{AB} =
\left(\frac{\pi}{\alpha+\beta}\right)^{3/2}
\exp\left[-\frac{\alpha \beta}{\alpha+\beta} (\vec{A}-\vec{B})^2\right] \ .
\\ \\
\hline
\end{array}
\label{Svysledok}
\EEEE
\end{large}
Matica kinetickej energie:
\begin{large}
\BBBB
\begin{array}{|c|}
\hline \\
T_{AB} = 
\frac{\alpha \beta}{\alpha + \beta} 
\left[3 - \frac{2 \alpha \beta}{\alpha + \beta} (\vec{A}-\vec{B})^2\right]
S_{AB}
\\ \\
\hline
\end{array}
\label{Tvysledok}
\EEEE
\end{large}
Matica Coulombovej interakcie pre jadro v bode $\vec{C}$ obsahujúce $Z_C$ 
protónov:
\begin{large}
\BBBB
\begin{array}{|c|}
\hline\\
V_{AB}^C = -Z_C \;
\frac{\textrm{erf}\left(\sqrt{\alpha+\beta}\ |\vec{P}-\vec{C}|\right)}
{|\vec{P}-\vec{C}|}\; S_{AB} \ ,
\\ \\
\hline
\end{array}
\EEEE
\end{large}
\underline{Poznámka}: Pre aplikáciu na $H_2^+$ vystačíme s~gaussiánmi 
centrovanými na~jadrách. Dve jadrá znamenajú dve centrá Coulombovej interakcie,
čiže budú dva rôzne vektory $\vec{C}$ značiace súradnice jadier. Centrá
$\vec{A}, \vec{B}$ gaussiánov ztotožníme s~centrami $\vec{C}$
Coulombovej interakcie. Ale je možné uvažovať aj o~bázových funkciách 
s~centrom mimo jadier, napr. Gaussovu väzbovú funkciu, čo je jednoduchý
gaussián typu $s$ centrovaný v~mieste povedzme uprostred chemickej väzby; 
v~prípade $H_2^+$ presne v~polovici medzi jadrami.
\footnote{Používajú sa aj tzv. plávajúce gaussiány \cite{carsky_urban}, ale
to je niečo trochu iné; pri~tom sa optimalizuje poloha centra gaussiánu
a aj šírka.}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{\label{sec:prakticka}Praktické poznámky k bázam\\
\small{(spracované podľa \cite{carsky_urban} a \cite{gaussian})}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%#############################################################################
\subsection{\label{ssec:zaklady}Základné pojmy}
%#############################################################################
\begin{itemize}
\item \textbf{minimálna báza}: jedna bázová funkcia pre každý obsadený atómový 
orbitál; napr. pre atóm kyslíka je minimálna báza tvorená bázovými funkciami
$1s, 2s, 2p_x, 2p_y$ a $2p_z$.
\item \textbf{double zeta báza}: dve bázové funkcie pre každý (zvyčajne
obsadený) atómový orbitál.
\item \textbf{rozšírená báza}: viac ako dve bázové funkcie pre každý 
atómový orbitál; zvyčajne obsahuje aj bázové funkcie pre neobsadené orbitály
(tzv. \textbf{polarizačné funkcie}); označenie napr. DZ+P. Často sa však
používajú aj pojmy \textbf{triple zeta}, \textbf{quadruple zeta} a
\textbf{quintuple zeta}.\\
\underline{Príklad}: DZ+P pre atóm vodíka je $(2s1p)$, čo znamená dve funkcie 
typu $s$, jednu 
funkciu typu p$_x$, jednu typu $p_y$ a jednu typu $p_z$.\\
\underline{Príklad}: DZ+P pre atóm kyslíka (ktorý má konfiguráciu
$1s^2, 2s^2, 2p^6$, čiže obsadené hladiny zahŕňajú dva rôzne orbitály 
typu $s$ a po jednom zo série
$(p_x, p_y, p_z)$ je $(4s2p1d)$, čo znamená štyri bázové funkcie typu $s$, 
po dve bázové funkcie typu $p_x$, $p_y$ a $p_z$ a jednu sadu bázových funkcií 
$d(xy, xz, yz, x^2-y^2, 3z^2-r^2)$.\\
\underline{Príklad}: DZ+P pre molekulu vody
je $(4s2p1d/2s1p)$, kde pred lomítkom je báza pre kyslík a za lomítkom
pre vodík.
\end{itemize}
%#############################################################################
\subsection{\label{ssec:kontrah}Kontrahované gaussovské bázy}
%#############################################################################
\underline{Príklad}: Molekula NH$_3$. Konfigurácia dusíka je
$1s^2, 2s^2, 2p^5$. Vhodnou bázou je $(7s3p/3s)$. Obsahuje 7 rôznych $s$ 
funkcií na atóme N, po 3 rôzne funkcie typov $p_x, p_y, p_z$ tiež na~N a 
3 rôzne $s$ funkcie na~každom atóme H. Teda celkovo 16 primitívnych 
gaussiánov len na dusíku a 3 na~každom vodíku.
Superpozíciu 7 gaussiánov na~N (bázu $7s$) nahradíme len dvomi zloženými
bázovými funkciami. Zaznačíme to
$(7s) \longrightarrow (12345) (67)$ alebo stručnejšie
$$
(7s) \longrightarrow [2s] \ .
$$
Podobne kontrahujeme aj každú trojicu rôznych $p_x$ na dusíku na~jednu
kontrahovanú funkciu typu $p_x$, čo zaznačíme $(3p) \longrightarrow [1p]$.
No a obdobne na~každom vodíku: $(3s) \longrightarrow [1s]$.
Výsledná kontrakcia pôvodnej bázy pre~NH$_3$ sa zapíše
$$
(7s3p/3s) \longrightarrow [2s1p/1s] \ .
$$
Občas sa používa aj iné značenie, vzťahujúce bázy k~jednotlivým atómom.
Napr. pre~dusík N$(7s3p/2s1p)$ a pre vodík H$(3s/1s)$, kde pred~lomítkom
je nekontrahovaná báza a za~lomítkom kontrahovaná. Ešte iné značenie je
N$(7s3p)/[2s1p]$ a H$(3s)/[1s]$.
%#############################################################################
\subsection{\label{ssec:STO}Aproximácia Slaterových orbitálov lineárnou
kombináciou gaussovských funkcií}
%#############################################################################
Nech sa pre každú funkciu STO použije rovnaký počet ($N$) gaussiánov
a orbitály $ns$ a $np$ nech majú pre~svoje gaussiány rovnaké rovnaké hodnoty
exponentov $\gamma_{n,k}$ v~nasledovnom zmysle:
\BBAA
\Phi_{1s}^{'}(1,r) &= \sum_k^N d_{1s,k} g_{1s}(\gamma_{1s,k}, r)
\nonumber\\
\Phi_{2s}^{'}(1,r) &= \sum_k^N d_{2s,k} g_{1s}(\gamma_{2s,k}, r)
\nonumber\\
\Phi_{2p}^{'}(1,r) &= \sum_k^N d_{2p,k} g_{2p}(\gamma_{2s,k}, r)
\nonumber
\EEAA
Takáto kontrahovaná báza sa značí ako 
$$
\textrm{STO-}N\textrm{G} \ .
$$
Je to pomerne jednoducho
zkonštruovaná báza (minimálna STO báza) a dáva dosť zlé výsledky.
Na~zlepšenie je potrebné zväčšiť flexibilitu bázy, teda spraviť 
dekontrakciu valenčnej šupky. Napr. pre~bázu STO-4G pre~atómy Li až F
zmiernime kontrakciu z~$(8s4p)/[2s/1p]$ na~$(8s4p)/[3s2p]$. 


Zo skupiny podobne rozštiepených báz sa používa napr. báza 
$$
\textrm{4-31G} \ .
$$
Význam jednotlivých znakov:
\begin{itemize}
\item[4] že každý orbitál vnútornej šupky je reprezentovaný jedinou 
kontrahovanou funkciou o~4 primitívnych gaussiánoch (tak isto ako u~STO-4G).
\item[-] oddeľuje vnútorné a valenčné šupky.
\item[31] že orbitály $ns$ a $np$ valenčnej šupky sú nahradené 
(dekontrahované) funkciami
$ns'$ a $np'$ o troch primitívnych gaussiánoch a funkciami 
$ns''$ a $np''$ (čo sú primitívne gaussiány).
\end{itemize}
Takáto kontrakcia sa značí aj $(4,3,1;\ 3,1)$.

Bázy ako 4-31G, teda s rozštiepenou valenčnou šupkou, však už majú málo 
spoločné so~Slaterovou bázou. 


Keď chceme vyššiu presnosť ako 4-31G, môžme
použiť napr. 6-31G. (Vnútorná šupka je reprezentovaná jednou kontrahovanou
funkciou o~6 primitívnych gaussiánoch). Táto báza má však flexibilitu
valenčnej šupky takú istú ako 4-31G, takže často neprináša žiadne významné
zlepšenie.
%#############################################################################
\subsection{\label{ssec:polar}Polarizačné funkcie}
%#############################################################################
(Viď aj odsek \ref{ssec:zaklady}). Napr. báza
$$
\textrm{6-31G}^* \ ,
$$
kde hviezdička znamená, že ku~báze 6-31G bol pridaný jeden súbor jednoduchých
gaussiánov typu $d$ (teda jedna kompletná šupka $d$) na~atómy Li až F.
Ak je aj druhá hviezdička, ako napr. v~báze
$$
\textrm{6-31G}^{**} \ ,
$$
tak bola pridaná polarizačná funkcia $p$ na atómy vodíka.
Namiesto dvojhviezdičkového značnia, ktoré môže byť mätúce, sa niekedy používa 
aj explicitnejšie: 6-31G($d$,$p$).
 
%#############################################################################
\subsection{\label{ssec:priklad}Priklad bázy}
%#############################################################################
Je to 6-31G$^*$ báza pre Cl a H \cite{gaussian}. SP znamená, že koeficienty
boli hľadané tak, že sa fitovali (metódou najmenších štvorcov) nie 
separátne na~$s$ a $p$, ale spolu. Prvý stĺpec sú exponenty, druhý 
rozvojové koeficienty. Báza bola optimalizovaná pre~teóriu hustotového
funkcionálu, konkrétne BLYP funkcionál.
\begin{small}
\begin{verbatim}
cl 0 
 S   6 1.00 
      0.2518010000D+05  0.1832959848D-02 
      0.3780350000D+04  0.1403419883D-01 
      0.8604740000D+03  0.6909739426D-01 
      0.2421450000D+03  0.2374519803D+00 
      0.7733490000D+02  0.4830339599D+00 
      0.2624700000D+02  0.3398559718D+00 
 SP   6 1.00 
      0.4917650000D+03 -0.2297391417D-02
      0.1169840000D+03 -0.3071371894D-01
      0.3741530000D+02 -0.1125280694D+00
      0.1378340000D+02  0.4501632776D-01
      0.5452150000D+01  0.5893533634D+00
      0.2225880000D+01  0.4652062868D+00
 SP   3 1.00  
      0.3186490000D+01 -0.2518280280D+00
      0.1144270000D+01  0.6158925141D-01
      0.4203770000D+00  0.1060184328D+01
 SP   1 1.00 
      0.1426570000D+00  0.1000000000D+01
 D   1 1.00  
      0.7500000000D+00 0.1000000000D+01 
**** 
h 0 
 S   3 1.00 
      0.1873113696D+02  0.3349460434D-01 
      0.2825394365D+01  0.2347269535D+00 
      0.6401216923D+00  0.8137573261D+00 
 S   1 1.00 
      0.1612777588D+00  0.1000000000D+01 
\end{verbatim}
\end{small}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibitem{qm}
J.~Pišút, L.~Gomolčák, V.~Černý: \textit{Úvod do~kvantovej mechaniky}
(Alfa, 2.~vydanie, Bratislava~1983).

\bibitem{numrec}
W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery:
\textit{Numerical Recipes in Fortran 77: The Art of Scientific Computing} 
(Cambridge University Press 1997). Dostupné cez 
\texttt{http://www.nr.com}.

\bibitem{carsky_urban}
P. Čársky, M. Urban: \textit{Ab initio výpočty v chemii} (SNTL Nakladatelství
technické literatury, Praha 1985).

\bibitem{gaussian}
\texttt{http://www.gaussian.com/g\_ur/k\_gen.htm}

\end{thebibliography}

\end{document}
