(*Napisane przez Pawła Czachorowskiego. Wykorzystanie dowolne.*)
(*Wszelkie uwagi mile widziane.*)

(*Kod ilustrujący rozwiązanie zadania 2 z pliku Ritz.pdf (baza 2 Gaussów z zadanymi wykładnikami. *)

(*Rozmiar bazy: 2. Możecie spróbować zwiększyć,tylko musicie zdefiniować kolejne wykładniki: a[3], a[4] itp.*)

SIZ=2

(*aprev - parametr funkcji Gaussa wyznaczony w poprzednim zadaniu - weźmy go sobie do pierwszej funkcji*)

a[1]=0.282942;

(*Bierzemy nową wartość wykładnika "z kapelusza". Odpowiada on dość "stromemu" Gaussowi - liczymy, że dobrze odda "szpic" funkcji*)

a[2]=2.0;
Plot[{Exp[-a[1]r^2],Exp[-a[2]r^2],Exp[-r]},{r,0,5},PlotLegends->"Expressions"]

(*Definicje całek*)

S[e1_,e2_]:=1.*(Pi/(a[e1]+a[e2]))^(3/2)

T[e1_,e2_]:=3.*Pi^(3/2)*a[e1]*a[e2]/(a[e1]+a[e2])^(5/2)

V[e1_,e2_]:=-2.*Pi/(a[e1]+a[e2])

H[e1_,e2_]:=T[e1,e2]+V[e1,e2]

(*Macierz hamiltonianu, w bazie nieznormalizowanej.*)

HMatNNorm=Table[H[i,j],{i,1,SIZ},{j,1,SIZ}];

MatrixForm[HMatNNorm]

(*Macierz całek nakładania, w bazie nieznormalizowanej.*)

SMatNNorm=Table[S[i,j],{i,1,SIZ},{j,SIZ}];

MatrixForm[SMatNNorm]

(*Znormalizujmy bazę. Tak naprawdę nie musimy tego robić! Ale skoro na ćwiczeniach to zrobiliśmy, to zróbmy i tu.*)

SMat=Table[S[i,j]/Sqrt[S[i,i]]/Sqrt[S[j,j]],{i,1,SIZ},{j,SIZ}];

MatrixForm[SMat]

(*Diagonalizacja Jacobiego macierzy symetrycznej M*)
(*Znajdujemy macierz diagonalizującą U *)
DiagM[M_]:={{Cos[alfa[M]],Sin[alfa[M]]},{-Sin[alfa[M]],Cos[alfa[M]]}}
alfa[M_]:=1/2ArcCot[(M[[2]][[2]]-M[[1]][[1]])/(2*M[[1]][[2]])]
(*Dla takiej prostej macierzy jak macierz całek nakładania S w bazie znormalizowanej (1 na diagonali) jest ona sama wyjątkowo prosta.*)
alfa[SMat]/Degree

U=DiagM[SMat]

MatrixForm[U]

SDiag=Inverse[U].SMat.U

(*Polecenie Chop usuwa śmieci numeryczne, które nie są zerami ze względu na skończoną precyzję komputera.*)

SDiag=Chop[SDiag]

(*"Pierwiastkujemy" macierz:*)

SDiag12=Sqrt[SDiag]

(*Odwracamy - po prostu 1/(dany element macierzy)*)
SDiagMinus12=Inverse[SDiag12]

(*"Oddiagonalizowujemy" w celu dostania S do minus 1/2*)

Sm12=U.SDiagMinus12.Inverse[U]

(*Teraz kolej na macierz hamiltonianu. Skoro postanowiliśmy uwzględnić normalizację funkcji w S, konsekwentnie musimy i tutaj.*)

HMat=Table[H[i,j]/Sqrt[S[i,i]]/Sqrt[S[j,j]],{i,1,SIZ},{j,1,SIZ}];

MatrixForm[HMat]

(*Obkładamy macierz H macierzami S do minus 1/2*)

HNew=Sm12.HMat.Sm12

(*...a teraz ją diagonalizujemy:*)

(*Znajdźmy U... *)

U=DiagM[HNew]

(*...i użyjmy do zdiagonalizowania macierzy HNew.*)

HDiag=Inverse[U].HNew.U

HDiag=Chop[HDiag]

MatrixForm[HDiag]

(*Dostaliśmy dwie energie - niższa jest tym, czego szukamy - przybliżeniem do stanu podstawowego H. Wyższa jest zupełnie beznadziejnym przybliżeniem do stanu wzbudzonego o tej samej symetrii. Ma wartość dodatnią - co jest absurdalne! Nam zależy jednak na stanie podstawowym, więc nas to nie martwi. :-P *)

(*Jak pamiętamy (mam nadzieję :-) ), kolumny macierzy U, którą zdiagonalizowaliśmy macierz hamiltonianu, to wektory własne H, czyli współczynniki rozwinięcia funkcji w naszej bazie.*)

cTylda=Transpose[U][[1]]

(*UWAGA! Aby uzyskać współczynniki rozwinięcia w bazie wyjściowej, musimy ten wektor jeszcze potraktować S do minus 1/2 (patrz Ritz.pdf - przecież przekształciliśmy uogólnione równanie własne).*)

c=Sm12.cTylda

(*Zatem nasza funkcja to (pamiętajmy, że używaliśmy w obliczeniach funkcji znormalizowanych):*)

FUN=c[[1]]/Sqrt[S[1,1]]Exp[-a[1]r^2]+c[[2]]/Sqrt[S[2,2]]Exp[-a[2]r^2]

(*Sprawdźmy, czy jest ok, wstawiając do równania wiekowego epsilon i współczynniki c:*)

Chop[(HMat-HDiag[[1]][[1]]*SMat).c]

(*Wyszedł wektor zerowy,tak jak powinien po prawej stronie równania wiekowego. *)

(*Możemy sobie jeszcze dla świętego spokoju porównać z wynikiem uzyskanym od początku do końca przez Mathematikę*)

WYNIK=Eigensystem[{HMatNNorm,SMatNNorm}]

MatrixForm[WYNIK]

(*W górnym wierszu są energie. Pod nimi wektory własne (współczynniki). Wydaje się, że są inne, ale pamiętajmy, że rozwiązanie jest określone z dokładnością do stałej i tak naprawdę liczy się stosunek tych współczynników. Sprawdźmy:*)
WYNIK[[2]][[2]][[1]]/WYNIK[[2]][[2]][[2]]

(*A dla naszego rozwiązania:*)
(c[[1]]/Sqrt[S[1,1]])/(c[[2]]/Sqrt[S[2,2]])

(*Czyli tyle samo. Yay! *)

(*Sprawdźmy, jak wygląda nasza funkcja w porównaniu z dokładną oraz tą z poprzedniego zadania:*)

Plot[{FUN, 1/Sqrt[S[1,1]]*Exp[-a[1]r^2],Sqrt[1/Pi]Exp[-r]},{r,0,5},PlotLegends->"Expressions"]

(*Możemy sobie jeszcze wyrysować (bo czemu nie?) jak "wygląda" atom w opisywanym przez nas stanie.
Kwadrat modułu funkcji falowej to gęstość prawdopodobieństwa znalezienia elektronu w danym punkcie. Kolorowanie wykresu dość arbitralne (w końcu tak naprawdę funkcja nie urywa się gdzieś koło 2a.u., tylko rozciąga w nieskończoność)*)
box=2.5;
DensityPlot3D[Abs[FUN]^2/.r->Sqrt[x^2+y^2+z^2],{x,-box,box},{y,-box,box},{z,-box,box},OpacityFunction->Function[f,(f)^1],PlotPoints->30,BoxRatios->Automatic,PlotRange->Full]

(*Mieliśmy przy okazji zobaczyć w praniu ortogonalizację Lowdina. Tak naprawdę przekształcając uogólnione równanie własne do równania własnego (z jednostkową macierzą całek nakładania) w sposób niejawny wykonaliśmy ten właśnie rodzaj ortogonalizacji. Informacja o bazie ortogonalnej siedzi w obliczonym przez nas Sm12*)
MatrixForm[Sm12]

(*Widać teraz, dlaczego jest to ortogonalizacja *symetryczna*! Nowe funkcje powstały ze starych poprzez jednakową (symetryczną) kombinację starych. *)

(*Możemy sprawdzić, czy nasza baza rzeczywiście jest ortogonalna. Znów musimy pamiętać, że znormalizowaliśmy bazę. *)

forto1=Sm12[[1]][[1]]/Sqrt[S[1,1]]*Exp[-a[1]r^2]+Sm12[[1]][[2]]/Sqrt[S[2,2]]*Exp[-a[2]r^2]

forto2=Sm12[[2]][[1]]/Sqrt[S[1,1]]*Exp[-a[1]r^2]+Sm12[[2]][[2]]/Sqrt[S[2,2]]*Exp[-a[2]r^2]

(*Całka nakładania dwóch różnych funkcji wynosi zero!*)

Chop[4*Pi*Integrate[forto1*forto2*r^2,{r,0,Infinity}]]

(*Sprawdźmy też normalizację funkcji - ortogonalizacja Lowdina powinna generować funkcje nie tylko ortogonalne, ale nawet ortonormalne:*)

Chop[4*Pi*Integrate[forto1*forto1*r^2,{r,0,Infinity}]]

Chop[4*Pi*Integrate[forto2*forto2*r^2,{r,0,Infinity}]]