33. DFT ab initio

Tytuł jest nieco przewrotny i ma sugerować, że rozważam teorię DFT od początku.

    Zastanawiam się, skąd się wzięła idea DFT. Podejrzewam, że genezą tej teorii jest dość prosta konstatacja, że podstawą obliczeń kwantowych powinny być wielkości obserwowalne. Równanie Schroedingera opiera się na abstrakcyjnej funkcji falowej, która w ogólności jest funkcją zespoloną. Wydaje się, że funkcje zespolone nie są tworami obserwowanymi w przyrodzie (choć istnieją próby dowodu odwrotnego). Dopiero utożsamienie przez Borna kwadratu modułu funkcji falowej z gęstością prawdopodobieństwa znalezienia cząstki w określonym obszarze pozwoliło (choć z oporami) "udomowić" owo równanie. Kwadrat modułu jest wielkością nieujemną obrazuje zatem istotną dla chemii gęstość elektronową. 

    Parr i Yang założyli, że podstawią obliczeń powinna być wielkość obserwowalna, jaką jest gęstość elektronowa [R.G. Parr, W. Yang, Density-Functional Theory of Atoms and Molecules, in The International Series of Monographs on Chemistry, Oxford University Press, Oxford, New York 1989]. Obowiązuje nadal schemat obliczeniowy równania Schroedingera, ale tak przekształcony, żeby w miejscu funkcji falowej pojawiła się gęstość i jej wielkości pochodne. Dla odróżnienia od formalizmu Schroedingera teorię tę nazwano teorią funkcjonałów, choć chodzi oczywiście o aparat matematyczny związany z obliczaniem i własnościami funkcji, której argumentem jest inna funkcja. Można postawić pytanie, co jest źródłem informacji o gęstości elektronowej w omawianej teorii. Otóż źródłem tym jest kwadrat modułu arbitralnie dobranej funkcji falowej. Do tej pory nie znaleziono innego źródła. Wydaje mi się jednak, że powinna istnieć inna procedura matematyczna, która mogłaby zastąpić funkcję falową. Procedura ta musi dawać w wyniku wielkość, którą można utożsamić z gęstością elektronową....

32. Automaty komórkowe

Tekst napisany w roku 2004

    Automaty komórkowe zainteresowały mnie na chwilę, przy okazji lektury książki Petera Coveneya i Rogera Higfielda Granice złożoności, Prószyński i S-ka, Warszawa 1997. Napisałem kilka programów w QBASICu, które odtwarzały zestaw reguł rządzących automatem. Poniżej przytaczam program, który losuje położenia punktów w "świecie automatu" i następnie realizuje schemat ewolucyjny według zadanych reguł.

'automat komórkowy. odtwarza dane książkowe. Zbudowany 2 marca 2001.
bok = 50        ' bok świata
lewol = 1000    ' liczba ewolucji
screen 10
 dim a(1 to bok, 1 to bok)
 dim b(1 to bok, 1 to bok)
'przypadkowy układ punktów
for j = 1 to 723
  x = int(rnd(1) * bok) + 1
  y = int(rnd(1) * bok) + 1
  a(x, y) = 1: pset (x * 2, y * 2)
next j
'pętla główna
for k = 1 to lewol
  for y = 2 to bok - 1
    for x = 2 to bok - 1
      gosub sasiedzi:
      if a(x, y) = 1 then
        if (suma = 2) or (suma = 3) then b(x, y) = 1
      else
        b(x, y) = 0
      end if
        if suma = 3 then b(x, y) = 1
    next x
  next y

'nowa tablica a i malowanie tablicy b na ekranie, zerowanie tablicy b
 for y = 1 to bok
   for x = 1 to bok
      a(x, y) = 0                   'zerowanie tablicy a
      a(x, y) = b(x, y)             'nadawanie nowych wartoœci tablicy a
      pset (x * 2, y * 2), 0        'mazanie starego ekranu
     if b(x, y) = 1 then pset (x * 2, y * 2), 1    'malowanie nowego ekranu
     b(x, y) = 0                    'zerowanie tablicy b
   next x
 next y
next k
end

sasiedzi:
'oblicza liczbę "pionków" wokół badanego punktu
  suma = 0
      for i = x - 1 to x + 1
        for j = y - 1 to y + 1
          if (i = x) and (j = y) then goto tu
          if a(i, j) = 1 then suma = suma + 1
tu:      next j
      next i
return

  Obserwacja zmian na ekranie jest ciekawa. Oczekuje się rychłego ustabilizowania się układu ewoluujących punktów, a tymczasem całość długo zmienia swoją postać. 

    Poniżej pokazuję ekran monitora po pewnej liczbie iteracji.

  Pewne układy początkowe stabilizują się po pewnej liczbie iteracji. Ciekawe zjawisko, zwłaszcza, że modyfikuje się pole wokół obiektu. Czuje się tu jakiś związek z fizyką cząstek. Tylko jaki?

31. Przygoda z chaosem

    

Moja przygoda z chaosem rozpoczęła się chyba podobnie, jak analiza zachowania się populacji zwierzaczków, robiona w 1975 roku przez Mitchella Feigenbauma, od zabawy z kalkulatorem. Dawno temu, może w czwartej klasie szkoły średniej, dostałem od Rodziców taką maszynkę. Ówcześnie to był czad! Miałem w rękach czarny, zgrabnie zaprojektowany i wykonany Texas Instruments z czerwonym wyświetlaczem. Byłem zdumiony, gdy zauważyłem, że na klawiaturze brakuje klawisza ze znakiem "równa się". Wtedy zetknąłem się po raz pierwszy z Odwrotną Notacją Polską. Dowiedziałem się o niej za pomocą amerykańskiej technologii, a nie przy pomocy polskich nauczycieli. Dziwna jakaś była ta edukacja w naszym Kraju.

    W ramach zabawy z nowym urządzeniem robiłem próby z iteracyjnym obliczaniem rozmaitych funkcji, ale bez żadnego mądrego podtekstu, chyba w ramach odpoczynku pomiędzy kolejnymi zadaniami z technologii chemicznej. Wprowadzałem pomyślaną liczbę i naciskałem klawisz jakiejś funkcji (sinusa, pierwiastka, itp.) tyle razy aż wynik się ustalił, albo wyświetlił się znak błędu. Jednak zachowanie się jednej z tych funkcji było odmienne od pozostałych.

    Gdy logarytmowałem jakiś ułamek, otrzymywałem wynik, który był oczywiście liczbą ujemną. Wydawało się, że to koniec zabawy, gdyż próba logarytmowania liczy ujemnej na kalkulatorze daje sygnał błędu. Zrobiłem jednak jeden krok dalej. Po prostu usunąłem znak minus sprzed wyniku i logarytmowanie wykonałem ponownie. Powtarzałem tę zabawę, ale okazało się, że nie umiałem zakończyć tak powtarzanej operacji. Ciągle pojawiały się jakieś liczby, mniejsze, albo większe od jedności. Operację mogłem powtarzać godzinami i nie udawało mi się dojść do jakiejś granicy... 

    Wszystko się jednak nudzi po jakimś czasie. Zapomniałem o tej zabawie na rok lub dwa. W tym czasie pojawiły się komputery ZX81 a potem Spectrum.

    Na studiach chemicznych nauczyłem się konstruować schematy blokowe procedur obliczeniowych i tworzyć odpowiednie programy komputerowe z wykorzystaniem języka BASIC. Lubię ten język... Wtedy przypomniałem sobie o niesfornym logarytmie. Schemat iteracji zapisałem następująco:

x_n+1 = abs(log(x_n))

    Za x₁ wybierałem liczbę większą od zera i "puszczałem" nieskończoną pętlę. Na ekranie monitora pojawiały się ciągi jakichś liczb, które nie chciały dążyć do żadnej konkretnej liczby. Program komputerowy zapisany w QBASICu* wyglądał tak:

'Program 00

x = .2 'na przykład 0.2 oczywiście

10 x = ABS(LOG(x))

PRINT x

GOTO 10

    Liczby migały na ekranie i, tak naprawdę, niewiele widziałem. Zmodyfikowałem nieco program, aby pokazywał generowane liczby w postaci punktów na osi x:

'Program 01 

SCREEN 12

x = .2

10 x = ABS(LOG(x))

PSET (x * 70, 100)

GOTO 10

    Powoli (na tamtych komputerach to było powoli) wyłaniała się kropkowana prosta, która miała dużo punktów przy początku osi. Punkty stawały się coraz rzadsze w miarę wzrostu wartości x. Po pewnym czasie obraz ustalał się:

Rysunek kreseczki, która pojawiła się po wykonaniu programu 01

    Tu przyszła mi do głowy pewna, całkiem kartezjańska myśl. Każdy punkt x_n+1 może być potraktowany jako punkt na osi prostopadłej do osi punktu x_n. To oznacza tyle, że wolno mi przyjąć, iż y = x_n+1, a każde y może być następnym x_n. Mogę po prostu narysować pojawiające się punkty w dwuwymiarowym układzie współrzędnych tak, jak rysuje się funkcje. Myśl tę zobrazowałem za pomocą prostego programu:

'Program 02

SCREEN 12

y = .2              'na przykład 0.2 oczywiście

skala = 30

10 x = y

y = ABS(LOG(x))

PSET (x * skala, y * skala)

GOTO 10

    Otrzymałem dwie krzywe, schodzące się w jednym punkcie. Krzywe te podobnie, jak poprzednio były gęsto usiane punktami przy początku układu współrzędnych. Dalej punkty stawały się coraz rzadsze.

Rysunek krzywych, powstałych w wyniku działania programu 02

    Do programu dodałem drobną modyfikację, która usuwa znak minus sprzed zmiennej y podczas wydruku punktu oraz usunąłem brzydką komendę 'goto', zastępując ją elegancką pętlą:

'Program 03

SCREEN 12

LINE (50, 300)-(200, 300), 5

LINE (100, 250)-(100, 350), 5 ' rysuje początkowy fragment układu

y = .2

skala = 10

FOR i = 1 TO 5000

  x = ABS(y)

  y = LOG(x)

  PSET (x * skala + 100, 300 - y * skala)

NEXT i

    Program odtwarza funkcję log(x). Dla ułatwienia interpretacji rysunku dopisałem do programu komendy rysujące fragment osi układu współrzędnych:

Rysunek funkcji z programu 03

    Odtworzyłem w jakiś sposób funkcję logarytmiczną, ale nie całą, jedynie jej część dość bliską początkowi układu współrzędnych. Kolejne się iteracje nie wpływały już na rysunek. Układ osiągnął coś w rodzaju stanu stacjonarnego. Punkty są gęsto usiane w pobliżu początku układu współrzędnych, dalej są rozmieszczone coraz rzadziej.

    Zrobiłem jeszcze jedną modyfikację. Dodałem zmienną c, przez którą mnożyłem x w każdej pętli FOR...NEXT:

'Program 04

SCREEN 12

  y = .2

  skala = 30

FOR c = .01 TO 1 STEP .01

  FOR i = 1 TO 1000

    x = ABS(y)

    y = LOG(c * x)

    PSET (x * skala + 100, 200 - y * skala)

  NEXT i

NEXT c

    Okazało się, że zmienna c wpływa na postać wykresu. Gdy c jest w pobliżu 0,01 ukazuje się tylko kilka rozrzuconych punktów. Potem punkty te skupiają się, a od wartości c około 0,35 skromny układ punktów rozwija się w znane już nam krzywe logarytmiczne. Zastanawiałem się, czy nie mam tu przypadkiem problemu z dokładnością obliczeń, i że niektóre punkty być może są wynikiem kumulacji błędów. Podejrzewałem jednak, że w ogólności otrzymany obraz, chyba pokazuje jakąœ prawidłowość, ale jaką?

Rodzina krzywych otrzymanych z programu 04 przy różnych wartościach zmiennej c.

    Nieco później dowiedziałem się o teorii chaosu. Klasyczny sposób popularnego wprowadzenia do tego zagadnienia, polega na pokazaniu pojawiania się chaosu w funkcji logistycznej:

x_n+1 = cx_n(1-x_n)

    W zależności od (ustalonej wartości) zmiennej c, iteracja prowadzi albo do pojedynczego rozwiązania, zwanego atraktorem, albo rozwiązanie rozdwaja się w punkcje zwanym bifurkacją (rysowane są dwa punkty), albo wreszcie rozwiązania stają się chaotyczne i trudno przewidzieć, jaka wartość będzie wynikiem następnej iteracji. Poniższy program pozwala na zaobserwowanie opisanych zależności:

'Program 05

'Atraktory, bifurkacje i chaos w r. logistycznym

DIM x(30)

SCREEN 12

k = 1000

FOR c = .9 TO 4 STEP .01

  x = .1

  FOR n = 1 TO k

    x = c * x * (1 - x)                   'równanie logistyczne

    IF n >= (k - 30) THEN x(k - n) = x

  NEXT n

  FOR i = 1 TO 30

    PSET (100 * c, 400 - x(i) * 300)

  NEXT i

NEXT c

Atraktory, bifurkacje i chaos dla funkcji logistycznej

    Poniżej przedstawiam znacznie późniejszą wersję tego programu (z roku 2004), która pokazuje o wiele więcej punktów chaosu (ale czy to jest zaleta, nie wiem) i jest prostsza od poprzedniej.

'Program 06

'atraktory, bifurkacje i chaos funkcji logistycznej

SCREEN 12

FOR c = .9 TO 4 STEP .01

  x = .2

  min = 60                            'liczba iteracji początkowych

  FOR i = 1 TO 1000

    x = c * x * (1 - x)              'równanie logistyczne

    IF i > min THEN PSET (c * 100, 400 - x * 300)

  NEXT i

NEXT c

    W tym miejscu zdałem sobie sprawę, że mogę przedstawić równanie logistyczne w prostokątnym układzie współrzędnych opartym o osie x oraz x_n+1. Kiedyś już to robiłem dla logarytmu. Ciekawiło mnie, jaka krzywa powstanie. Program, którym to policzyłem, różnił się niewiele od przedstawionego powyżej.

'Program 07

'atraktory, bifurkacje i chaos

SCREEN 12

FOR c = .9 TO 4 STEP .01

  y = .2

  min = 100                       'liczba iteracji początkowych

  FOR i = 1 TO 2000

    x = y

    y = c * x * (1 - x)           'r. logistyczne

    IF i > min THEN PSET (x * 200 + 100, 350 - y * 300)

  NEXT i

NEXT c

    Wynik działania tego programu przedstawiłem poniżej. Pojawiają się nowe informacje o zachowaniu się atraktorów i rozwiązań chaotycznych.

    Linia ukośna wewnątrz łuku paraboli pokazuje sytuację, w której rozwiązaniem jest atraktor przy wzrastających wartościach c. Potem pojawia się rozgałęzienie na dwa ramiona, co odpowiada pierwszej bifurkacji. Ramiona są kawałkiem jakiejś krzywej, co jest wyraźnie widoczne. Kolejne rozgałęzienie jest mniejsze. Następne rozgałęzienia stają się jeszcze mniejsze aż układ "rozwija się" nagle w parabolę! Kolejne wartości c zmieniają wysokość tej paraboli.

    Za pomocą następnego programu można prześledzić powstawanie jednej z takich parabol przy ustalonej wartości c, oraz po różnej liczbie iteracji, większej od minimalnej liczby iteracji, która zapewnia zbieżność do lokalnego atraktora.

'Program 08

'atraktory, bifurkacje i chaos

SCREEN 12

c = 4

  y = .2

min = 100                              'liczba iteracji początkowych

FOR i = 1 TO 120                       ' potem 200, 400 oraz 2000

  x = y

  y = c * x * (1 - x)                        'r. logistyczne

  IF i > min THEN PSET (x * 200 + 100, 350 - y * 300)

NEXT i

Rysunki paraboli po 120, 200 i 4000 iteracjach

    Nasuwają się ciekawe wnioski. Rozwiązania chaotyczne niosą najwięcej informacji o przestrzeni, w której rozgrywa się dramat iteracji równania nieliniowego. Oczywiście dzieje się to pod warunkiem, że przeprowadzimy dużą liczbę iteracji. Widzimy wyraźnie, że odtwarza się parabola (albo wcześniej funkcja logarytmiczna). Rozwiązania chaotyczne odtwarzają tylko część krzywej matematycznej w pobliżu początku układu współrzędnych. Dokładniej, część leżącą w zakresie nieujemnych wartości funkcji logistycznej. Chaos jest jednak informatywny w odróżnieniu od rozwiązań zawierających atraktory. Zdaje się również, że ustala się pewien rozkład punktów w odtwarzanej przestrzeni, co oznacza, że w ramach przedstawionego algorytmu chaos nie jest idealny (prawdopodobnie rozmieszczenie punktów nie jest równomierne). Powyżej pewnej liczby iteracji następne iteracje nie zmieniają rysunku, co może oznaczać, że został osiągnięty pewien "stan stacjonarny". To ostatnie spostrzeżenie jest intrygujące. Załóżmy, że istnieje pewien proces fizyczny, opisany jakąś funkcją, której przebieg można odtworzyć za pomocą iteracji chaotycznej. Wolno się spodziewać, że taki układ fizyczny osiągnie wspomniany stan stacjonarny. Dodatkowo, jeśli nawet dziedzina funkcji rozciągnięta jest na całą przestrzeń, to rozwiązania chaotyczne są lokalne i skupione w pobliżu początku układu współrzędnych. 

  Taka lokalność jest dla mnie szczególnie intrygująca i prowadzi do pytania, czy procesy dotychczas opisywane funkcjami, których dziedzina jest nieskończona, muszą koniecznie "włączać" ową nieskończoność do tego opisu?  Mogę zapytać inaczej, czy procesy kwantowe opisane ciągłymi funkcjami falowymi mogą posiadać chaotyczne rozwiązania, które znacznie lepiej odpowiadają idei kwantów niż owa ciągłość. Być może nie jest konieczne angażowanie całej przestrzeni do obliczenia właściwości układów kwantowych? Być może włączenie chaosu do takich rozważań zamieni niepokojącą nielokalność w lokalność, która jest bliższa naszemu ograniczonemu pojmowaniu świata. Ciekawi mnie to.  

* Współcześnie można uruchomić te paragramy przez emulację środowiska DOSu. Potrzebne informacje są na stronie QBASIC

30. Obliczanie energii stanu podstawowego atomu wodoru - równanie Schrödingera

Równanie Schrödingera dla atomu wodoru

    W celu obliczenia energii atomu wodoru weźmy równanie Schrödingera w przybliżeniu nieruchomego i nieskończenie ciężkiego jądra. Nieruchome jądro pozwala zaniedbać energię kinetyczną protonu. Nieskończenie ciężkie jądro umożliwia pewne uproszczenie obliczeń przez zastosowanie jedynie  masy elektronu. Tak skonstruowane równanie ma następującą postać:

Gdzie:

Dodatkowo, operator promienia wodzącego elektronu dany jest wzorem:

    W tym: h - stała Plancka (tu podzielona przez 2pi); m - masa elektronu; e - ładunek elementarny; epsilon_zero - stała dielektryczna próżni; r - promień wodzący elektronu.

Dobór funkcji falowej psi

    W celu rozwiązania tego równania należy dobrać funkcję próbną. niech to będzie funkcja eksponencjalna, która posiada element elastyczny, czyli zmienną wariacyjną c:

Wykonanie działania hamiltonianu na funkcję psi

     Działanie to sprowadza się do obliczenia postaci drugich pochodnych funkcji psi po zmiennych x,y,z. Działanie operatora promienia sprowadza się do pomnożenia operatora przez stałe. Do obliczenia drugich pochodnych proponuję zastosować jeden z programów do obliczeń symbolicznych np: Mathcad, Maxima, Mathematica, itp.

    Obliczenia prowadzą do następujących wzorów:

Pamiętajmy, że pierwiastek z sumy kwadratów współrzędnych jest promieniem wodzącym elektronu r. Możemy dokonać tej podmiany i zsumować drugie pochodne. W wyniku tego działania i dodania członu operatora energii potencjalnej otrzymuje się następujące wyrażenie:

 (patrz dopisek na końcu)

   Wyciągnięcie exp(-cr) przed nawias prowadzi do bardziej zwartego wzoru:

Mnożenie lewostronne

   W tym miejscu uzyskane równanie należy pomnożyć lewostronnie przez funkcję sprzężoną do psi. tutaj użyjemy funkcji rzeczywistych, więc funkcja sprzężona jest równa danej funkcji. Proszę zwrócić uwagę, że fakt, czy mnożymy lewo, czy prawostronnie nie ma już znaczenia. Tym niemniej dla zachowania porządku poniższy zapis uwzględnia lewostronność mnożenia:

    Po oczywistym wymnożeniu zapis można przedstawić następująco:

Całkowanie równania

    Tak uzyskane równanie należy scałkować po całej dziedzinie zmiennych x,x,z.

    Tego typu całkę we współrzędnych kartezjańskich trudno jest rozwiązać analitycznie. Wygodnej jest przejść do współrzędnych sferycznych. W tym celu zmienne kartezjańskie należy zastąpić zmiennymi sferycznymi, a równanie pomnożyć obustronnie przez jakobian przekształcenia. Dokonanie tych zmian prowadzi do następującego równania:

    Po podstawieniu jakobianu i uporządkowaniu znaków równanie przyjmuje postać:

  Całkę potrójną należy przedstawić w postaci iloczynu trzech całek pojedynczych:

    

    Całkowanie po kątach fi i teta daje współczynnik 4π. Wolno zatem ten współczynnik skrócić, gdyż występuje po obu stronach powyższego równania (oznacza to, że rozwiązania są niezależne od kątów, czyli kulistosymetryczne). Pozostaje zatem całkowanie po promieniu wodzącym elektronu. Promień ten rozciąga się od zera do nieskończoności:

    Wzór na energię przybiera następującą postać:

Po uporządkowaniu stałych i wyciągnięciu zmiennej r przed nawias, oraz po dokonaniu całkowania otrzymujemy (drugi wzór):

   

Uproszczenie wyrażenia na energię prowadzi do prostego wzoru, który jest uzależniony od wartości zmiennej c oraz wartości stałych a  oraz b oczywiście:

Szukanie minimalnej wartości energii

   W celu policzenia energii minimalnej względem c należy obie strony równania zróżniczkować po c:

    Stąd już łatwo jest obliczyć wyrażenie na energię minimalną atomu wodoru przez podstawienie wyrażenia na c do równania na energię:

   Można się przekonać po podstawieniu wartości a  oraz b, że uzyskany wynik odpowiada energii stanu podstawowego atomu wodoru:

     Warto wykonać podobne obliczenia dla elastycznej funkcji Gaussa typu psi = exp(-cr²).

   
    Na marginesie dodam, że jeśli obliczenia wykonać w jednostkach atomowych, to wówczas stała a przybierze wartość 1/2 oraz b = 1, wówczas c = 1, oraz E = -1/2. Energia podana jest tu w jednostkach hartree.

Dopisek. Należy sobie zadać w tym miejscu pytanie, dlaczego nie skraca się obu stron równania przez czynnik exp(-cr). Otóż robi się tak w pewnym sposobie obliczania energii układów kwantowych, który nazywa się metodą energii lokalnych E = Hpsi/psi. Wykorzystuje się tu najczęściej procedurę Monte Carlo z algorytmem Metropolisa. Omówiłem to zagadnienie w innym miejscu.