W poprzednim odcinku rozrysowywałem w Excelu jaki w zależności od temperatury, rozkładają się dwa fermiony, gdy każdy z nich ma do dyspozycji cztery możliwe stany. Podałem wówczas, że liczba stanów układu złożonego jest równa 6. W sytuacji kiedy mamy k fermionów i n stanów do obsadzenia przez każdy fermion, liczba stanów układu jest równa liczbie kombinacji bez powtórzeń:
Liczby te są elementami trójkąta Pascala:
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1 ← ta szóstka jest właśnie przypadkiem opisywanym w poprzedniej notce
Początkowo to niewielkie liczby, ale szybko robią się duże, bo przeliczenie sumy dla 10 cząstek na 20 możliwych poziomach da 184756 możliwych stanów układu złożonego. Nic mądrego w Excelu wyliczyć się chyba nie da, napisałem więc program w C++.
Model a rzeczywistość
Zadanie polega na znalezieniu wszystkich k-elementowych konfiguracji, gdy każda cząstka ma do dyspozycji n poziomów. Każda konfiguracja – inaczej mówiąc stan układu złożonego – ma swoją energię E. Trzeba wyliczyć wagi – unormowane exp{–E/kT} – dla każdego stanu, gdy będziemy chcieli liczyć średnie jakichś obserwabli – z notki „Stany mieszane kwantowo” dowiemy się, że owe średnie to właśnie średnie ważone.
W poprzedniej notce liczyłem średnią zajętość danego poziomu – n(i) – sumując iloczyn prawdopodobieństwa (wagi Boltznama) i zajętości (0 lub 1) dla każdej konfiguracji. Teraz zrobię to samo, tylko dla większych n i k. Jeśli nie zrozumiesz po co to całe liczenie, to pozostaje ci zajrzeć do poprzedniego bardziej „łopatologicznego” odcinka.
Podobnie jak w poprzedniej notce założę, że energie pojedynczej cząstki zmieniają się co stałą wartość: Ei=i∙E0. Czy taki wybór ma jakiś sens fizyczny? I tak i nie. W bardzo znanym modelu gdzie rozkład Fermiego-Diraca jest wykorzystywany czyli modelu elektronów w metalu, energia rośnie z grubsza kwadratowo wraz ze wzrostem numeru stanu (rolę liczby i pełni tam tzw. wektor falowy). Uwzględnienie specyfiki 3D – bo ciało stałe jest trójwymiarowe – powoduje, że ilość stanów do obsadzenia jest proporcjonalna do pierwiastka z energii poziomu. A jak dołożymy wpływ od sieci krystalicznej rozkład elektronów jeszcze bardziej się komplikuje. Tak więc model, o którym dziś piszę jest drastycznym uproszczeniem rzeczywistości.
Mimo to, można pokazać, że rozkład elektronów w ciele stałym jest w postaci iloczynu:
n(E) = N(E) f(E)
Funkcja N(E) to teoretyczny rozkład „wszystkie możliwe stany są pozajmowane”. Ponieważ żeby zająć wszystkie stany, trzeba by mieć nieskończenie wiele elektronów, więc ten prawdziwy, czyli n(E) będzie miał mniejsze wartości. Funkcja „obcinająca” f(E) jest właśnie funkcją rozkładu Fermiego-Diraca. Zależy ona od temperatury i jest dość nieczuła na to, jak poukładane są możliwe energie[1]. Ponieważ w moim modelu energie poukładane są równomiernie, gęstość N(E) jest stałą i obliczony rozkład n(E) będzie po prostu równy rozkładowi Fermiego-Diraca. Takiemu co to rzeczywiście ujawnia się w metalach.
Rozkład Fermiego-Diraca liczy się na kartce, za pomocą wzorów – można przeczytać o tym w dowolnym podręczniku dotyczącym własności ciał stałych[2]. Mimo to, chcę policzyć ów rozkład metodą „brute force” i zobaczyć czy to „co wyjdzie z obliczeń” jest podobne do wykresów jakie można spotkać w podręcznikach:
Wykres z „Fizyki ciała stałego” Ashcrofta i Mermina – „biblii” ciałostałowców.
Tak naprawdę, to wybrałem obserwable zajętości poziomów dlatego, że wyniki łatwo porównać z dowolnym podręcznikiem tego typu.
Pisarze do piór, amatorzy do klawiatur…
Zawołanie chyba nietrafione bo dziś pisarze też chyba częściej używają klawiatury niż pióra. Ale druga jego część jest jaka najbardziej w porządku. Jako amator zasiadłem więc przed klawiaturą i napisałem program tworzący wszystkie stany k-cząstkowe przez odpowiedni wybór który z n poziomów jest zajęty przez daną cząstkę. Czyli najpierw musiałem utworzyć wszystkie kombinacje bez powtórzeń. No cóż, chwilkę pozastanawiałem się jakiż to może być algorytm, w końcu znalazłem gdzieś w sieci jego opis.
Dla ustalenia uwagi pomyślmy, że mamy zbiór liczb 1,2,3, … N i musimy wygenerować wszystkie podzbiory K-elementowe. W programie powybierane liczby umieszczone będą w K-elementowej tablicy. Sztuka polega na tym, żeby wygenerować wszystkie możliwe wybory.
Pierwszą kombinacją będzie zestaw:
Potem będziemy generować kolejne kombinacje, tak żeby ostatnią była:
Po drodze program musi wytworzyć wszystkie możliwe zestawy. A oto sam przepis: Sprawdzamy od końca kolejne elementy tablicy (od ostatniego do pierwszego) i szukamy pierwszej liczby różniącej się od ustawienia końcowego. Jeśli na taką trafimy, to:
- zwiększamy ją o 1,
- w każdej następnej komórce ustawiamy liczbę o 1 większą od liczby z komórki poprzedzającej.
No i w ten sposób w tablicy mamy nową kombinację. Żeby dostać następną, powtarzamy powyższe czynności. Generowanie kolejnych zestawów kończy się, jeśli nie znajdziemy odstępstwa od spodziewanej ostatniej kombinacji.
Rzecz jasna przy pierwszych iteracjach zmieniany będzie tylko ostatni element tablicy, bo nie ma „następnych komórek”. Dość łatwo zauważyć, że algorytm zachowuje sortowanie elementów tablicy, pod warunkiem, że początkowa kombinacja jest posortowana. Fakt, że w ten sposób wygenerujemy wszystkie kombinacje nie jest już tak widoczny, ale działa. Sprawdziłem.
Wprowadziłem w życie powyższy algorytm, pamiętając, że tablice w C++ numerowane są od zera – pierwsza jest zerową, druga pierwszą itd. Poniżej cytuję kawałek tego programu, gdyby ktoś na przykład takiego potrzebował ☺:
const int N=6;
const int K=3;
struct Kombinacja
{
int tab[K];
void init() { for( int i=0; i<K; i++ ) tab[i]=i+1; }
Kombinacja() { init(); }
bool nast();
int const operator[]( int i ){ return tab[i]; }
};
bool Kombinacja::nast()
{
int buf;
for( int i=0; i<K; i++ )
{
if( tab[K-1-i] != N-i )
{
tab[K-1-i]++;
for( int j=K-i; j<K; j++ ) tab[j] = tab[j-1]+1;
return true;
}
}
return false;
}
Pozamykałem sobie całość kodu w klasę, bo zdaje się lepszy z tym porządek. Ustawienie początkowe zapewnia konstruktor (lub metoda init()), a dostęp do wygenerowanych liczb przeładowany operator indeksowania tablicy. Generowanie kolejnych kombinacji odbywa się sekwencyjnie: jak już nie ma co generować, bo doszliśmy do końca metoda nast() zwraca fałsz. A oto przykład zastosowania klasy – skutkuje wypisaniem każdej kombinacji w osobnej linijce:
Kombinacja komb;
do
{
for( int i=0; i<K; i++ ) cout << " " << komb[i];
cout << "\n";
} while( komb.nast() );
Jeszcze jedno: ponieważ wygenerowane liczby mają mówić, które poziomy w stanie są zajęte, to potrzebne mi będzie wybieranie liczb z zakresu od 0 do N-1, bo konfigurację stanu też zapiszę w tablicy, a tablice, jak już napisałem numeruje się od 0. Poprawki będą niewielkie. Trzeba zmienić linijkę przy generowaniu startowej kombinacji:
for( int i=0; i<K; i++ ) tab[i]=i;
oraz linijkę przy sprawdzaniu w którym miejscu bieżąca kombinacja różni się od końcowej:
if( tab[K-1-i] != N-1-i ) ...
Właściwy program
Wypisałem tyle szczegółów o kombinacjach bez powtórzeń, że tylko skrótowo opiszę jak wygląda zasadnicza część programu. Konfiguracje przechowuję w tablicy o rozmiarze N: Wartość równa 1 i-tej komórki oznaczała, że poziom o tym numerze jest zajęty. Wartość 0, że wolny. Mając kombinację opisaną powyżej tworzę konfigurację następująco:
for( int idx=0; idx<K; idx++ ) konf[komb[idx]]=1;
Tablica konf[] zawiera bieżącą konfigurację. Po uzyskaniu stanu złożonego, trzeba znaleźć jego energię:
energia=0;
for( int i=0; i<N; i++ ) energia += E[i]*konf[i];
jak łatwo się domyślić w tablicy E umieściłem wartości energii każdego ze stanów jednocząstkowych. Wagi znajdowałem jako:
pstwo = exp( (energia0-energia)/T );
Jak widać przeskalowałem liczone wagi (i tak na końcu trzeba je unormować). Wynikło to stąd, że owe wagi były bardzo małe, a w trakcie testowania programu chciałem je sobie „podglądać”. Ponieważ wyniki niekoniecznie mi coś mówiły, dla lepszej czytelności, przeskalowałem sobie wagi tak, by stan o najniższej energii miał wagę 1. Oczywiście stałą energia0 trzeba wcześniej obliczyć:
for( int i=0; i<K; i++ ) energia0 += E[i];
Zajętość i-tego poziomu przeliczana była na bieżąco. Podobnie zresztą jak suma Z, potrzeba do unormowania całości na koniec:
//dla każdej konfiguracji trzeba wykonać:
Z += pstwo;
for( int i=0; i<N; i++ ) rozklad[i] += pstwo*konf[i];
W ten sposób przeliczyłem rozkłady Fermiego-Diraca (trzeba było unormować wartości w tablicy rozklad[] dzieląc je przez Z) dla małych n i k. Wyniki widać na dwóch obrazkach:
Ten pierwszy (n=25, k=10) wymagał wygenerowania nieco ponad miliona kombinacji – trwało to chwilkę. Próba wykonania podobnego dla n=30 i k=15 okazało się średnio wykonalne. To znaczy program się uruchomił i 155 117 520 przeliczeń zajęło mu 1010 sekund (prawie 17 minut). Można się oczywiście zastanowić, gdzie można wprowadzić jakieś optymalizacje[3], by trwało to szybciej, ale składowe trójkąta Pascala rosną na tyle szybko, że nawet dwukrotne skrócenie czasu działania nie sprawdzi się dla większych n i k. Dlatego też w następnej notce spróbuję zastosować jakąś metodę przybliżoną.
[1] Słowa „dość nieczuła na to, jak poukładane są możliwe energie” oznaczają, że oprócz temperatury, drugim parametrem funkcji Fermiego-Diraca jest tzw. potencjał chemiczny, zapisywane literką μ. W tej notce i w następnych, gdy maluję wykres funkcji obliczona ze wzoru, podstawiam za potencjał chemiczny energię Fermiego, czyli największy zajęty poziom w konfiguracji odpowiadającej T=0.
[2] Całkiem uczciwe wyprowadzenie musi zahaczyć o tzw. wielki rozkład kanoniczny. A może i nie, skoro w cytowanej obrazkowo książce Ashcrofta i Mermina znajdziemy „nieortodoksyjne”, jak piszą autorzy, wyprowadzenie przy pomocy zespołu kanonicznego. Niemniej w tym ujęciu nie widać pełnego znaczenia potencjału chemicznego μ.
[3] I dało się nawet takie zrobić, dzięki czemu czas obliczeń (n=30 i k=15) skrócił się do 840 sekund. Program uruchamiałem na zwykłym PC-ie.
