5. Energia stanu podstawowego atomu wodoru według metody Monte Carlo. Połączone wprawki II-VII - Aktualizacja 01.2020.

Poprzednia wprawka pozwoliła mi zrozumieć istotę metody Monte Carlo. Teraz postanowiłem zastosować ją do szkolnego, ale nieco bardziej zaawansowanego zagadnienia. Stany energetyczne atomu wodoru można obliczyć z wykorzystaniem wzorów Bohra lub równania Schroedingera. Tu zastosowałem metoda Monte Carlo, która jest popularnym w fizyce i chemii sposobem obliczania stanów rozmaitych układów kwantowych. W celu obliczenia stanu podstawowego atomu wodoru wykorzystałem algorytm opublikowany przez prof. A. Dawida (obecnie Wyższa Szkoła Biznesu w Dąbrowie Górniczej). Sporo czasu zajęło mi rozpracowanie oryginalnego algorytmu i zrozumienie działania jego składowych. Po kosmetycznych zmianach kodu zapisałem go w języku VBA dla Excela. Poniżej zapisany program oblicza zależność energii stanu podstawowego atomu wodoru od wartości zmiennej wariacyjnej. Minimalna wartość energii dla obliczenia ścisłego wynosi -0.5 hartree.

Public c, pi As Double
Sub atomH() 'procedura oblicza energię stanu podstawowego atomu wodoru
Dim x, y, z, h, h2 As Double
pi = 3.14159265
llos = 1000000            'liczba losowań
h = 0.01                 'krok różniczkowania
h2 = h * h               'kwadrat kroku różniczkowania
j = 1                    'licznik rzędu w tablicy Excela
For c = 0.2 To 1.9 Step 0.1 'obl. energii w zal. od zm. wariacyjnej
  En = 0                   'zerowanie energii do sumowania
  x = 0: y = 1.2: z = 0    'punkt początkowy
    For i = 1 To llos      'pętla obliczeniowa
    xzapas = x             'zabezpieczenie starych współrzędnych
    yzapas = y
    zzapas = z

    psi = fufal(x, y, z)    'wartość psi w punkcie (x,y,z)
    psikw = psi * psi       'kwadrat psi

           'Losowanie nowego punktu  (xn,yn,zn)
    xn = 2 * Rnd - 1        'zakres (-1, 1)
    yn = 2 * Rnd - 1
    zn = 2 * Rnd - 1

    pr = Sqr(xn * xn + yn * yn + zn * zn)  'promień punktu (xn,yn,zn)
    skok = 0.6 * GRNF        'skok przypadkowy wg rozkładu Gaussa

    x = x + xn * skok / pr  'poprawione współrzędne nowego punktu
    y = y + yn * skok / pr
    z = z + zn * skok / pr
 
    psinowa = fufal(x, y, z)     'wartość psi w nowym punkcie
    psinowakw = psinowa * psinowa

'obliczanie warunku przejścia łańcucha do nowego punktu

   rel = psinowakw / psikw   'stosunek nowej i starej psi
    If rel < 1 Then          'przejście bezwarunkowe gdy rel >= 1
      If Rnd > rel Then      'przejście odrzucone
        x = xzapas
        y = yzapas
        z = zzapas
        psinowa = fufal(x, y, z)
      End If
   End If

 'obliczanie sumy energii lokalnych
 'obliczanie drugich pochodnych cząstkowych w punkcie (x,y,z)
  fufal2 = 2 * fufal(x, y, z)
  d1 = (fufal(x - h, y, z) - fufal2 + fufal(x + h, y, z)) / h2
  d2 = (fufal(x, y - h, z) - fufal2 + fufal(x, y + h, z)) / h2
  d3 = (fufal(x, y, z - h) - fufal2 + fufal(x, y, z + h)) / h2
  Lappsi = -(d1 + d2 + d3) / 2 'druga poch. psi w punkcie (x,y,z)
  V = -1 / (Sqr(x * x + y * y + z * z))   'potencjal, tu  psi/psi=1
  En = En + Lappsi / psinowa + V  'sumowanie energii lokalnych
 Next i
  Cells(j, 1) = c
  Cells(j, 2) = En / llos 'zapis energii średniej
  j = j + 1
Next c
End Sub

Function fufal(x, y, z) 'oblicza wartość funkcji psi w punkcie x,y,z
  r = Sqr(x * x + y * y + z * z)
  fufal = Exp(-c * r)
End Function

Function GRNF() 'generuje liczbę pseudolosową z rozkładu Gaussa
  r1 = Sqr(-Log(1 - Rnd))
  GRNF = r1 * Sin(2 * pi * Rnd)
End Function

Wykres przedstawia wyniki obliczeń dla miliona losowań. Minimum znajduje się przy c = 1 i wynosi prawie dokładnie -0.5 hartree.

W poprzednim wpisie przedstawiłem algorytm obliczania energii atomu wodoru z wykorzystaniem metody Monte Carlo. Energie lokalne liczone były za pomocą numerycznego różniczkowania funkcji falowej w punkcie (x,y,z). Celem tego różniczkowania było obliczenie wartości laplasjanu w tym punkcie.
Ponieważ różniczkowanie takie można wykonać symbolicznie, procedurę można uprościć wpisując do programu symboliczną postać wyniku działania hamiltonianu na wybraną funkcję psi. Takie podejście jest mniej ogólne od numerycznego wyznaczania wartości laplasjanu, ale może ułatwić zrozumienie sposobu działania algorytmu Metropolisa.
Hamiltonian działający na funkcję psi ma postać (jednostki atomowe):

Dla funkcji psi postaci

Wynik działania hamiltonianu na tę funkcję wyraża się wzorem:

Wobec powyższego wzór na energię lokalną, używaną w algorytmie Monte Carlo, ma postać:

Fragment programu umieszczony w poprzednim wpisie, a odnoszący się do modułu całkowania (właściwie sumowania) upraszcza się do postaci:

...
'obliczanie cakłi z energii lokalnych
  r1 = SQR(x * x + y * y + z * z)
  En = En + (-c * c / 2 + (c - 1) / r1)
...

Wartość energii -0.5, jako minimum dla wartości c=1 otrzymuje się po dokonaniu sumowania 20000 konfiguracji.

W poprzednim poście pokazałem zapis modułu całkowania energii lokalnych z wykorzystaniem analitycznej postaci wyniku działania hamiltonianu dla atomu wodoru na funkcję exp(-cr). W niniejszym poście wykorzystałem to samo podejście i wyznaczyłem postać wyniku działania hamiltonianu dla atomu wodoru na funkcję gaussowską exp(-cr²):

Wzór na energię lokalną ma postać, jak niżej:

Funkcja programu QBASIC definiująca funkcję Gaussa w programie ma postać:

FUNCTION fufal (x, y, z)
  r = SQR(x * x + y * y + z * z)
  fufal = EXP(-c * r * r)
END FUNCTION

Odpowiedni moduł całkowania numerycznego sprowadza się do następujących instrukcji:

...
'obliczanie calki z energii lokalnych
  r1 = SQR(x * x + y * y + z * z)
  El = -2 * c * c * r1 * r1 - 1 / r1 + 3 * c
  En = En + El
...

Dla 500000 konfiguracji minimum dla c=0.3 wynosi -0.421 j.a., co w przeliczeniu na elektronowolty daje wartość -11.45 eV. Błąd wyliczenia w stosunku do dokładnej wartości energii wawynosi +15.8%.

W celu dokładnego zrozumienia sposobu obliczania energii atomu wodoru z wykorzystaniem metody Monte Carlo i algorytmu Metropolisa, napisałem program w języku QBASIC. Program wykorzystuje znany generator liczb pseudolosowych o rozkładzie normalnym. Aby przekonać się, że tak się dzieje w rzeczywistości wyniki losowań zapisałem w pliku, z którego wykonałem wykres za pomocą programu Excel. W programie umieściłem kod, który pozwala zliczać losowania w wektorze x:

'program GRNF generator lipselo o rozkladzie Gaussa
'zakres GRNF od -3.4 do 3.11 dla miliona losowan

pi = 3.14159265#

DIM x(-350 TO 350)

FOR i = 1 TO 100000
  r1 = SQR(-LOG(1 - RND))       'generowanie lipselo
  GRNF = r1 * SIN(2 * pi * RND)
  liczba = INT(GRNF * 100)      'zliczanie
  x(liczba) = x(liczba) + 1
NEXT i

OPEN "gauss.dat" FOR APPEND AS #1

FOR i = -350 TO 350
  PRINT #1, i, x(i) 'wpis wyniku do pliku
NEXT i

CLOSE

Zamieszczony poniżej wykres pokazuje, że generowane liczby pseudolosowe mieszczą się w zakresie od około -3.4 do około 3.2. Wartość średnia losowań jest bliska zera a maksymalna liczebność losowań średniej wynosi około 0.6% całkowitej liczby losowań:

Przedstawiony poprzednio wykres zależności częstości występowania promienia od jego wartości przeanalizowałem z wykorzystaniem postaci analitycznych funkcji. Kilka prób wykonanych za pomocą programu Derive pokazało, że dobrym przybliżeniem jest iloczyn funkcji Gaussa i kwadratu zmiennej niezależnej z odpowiednimi współczynnikami. Uzyskany wynik pokazuje poniższy wykres:

Uwaga. Ta aproksymacja nie jest zbyt dobra. Lepsza została pokazana w następnym wpisie.

Metoda Monte Carlo obliczania energii stanów atomu wodoru wykorzystuje algorytm Metropolisa wyboru punktów w przestrzeni wokółjądrowej. Dzięki temu wyborowi suma energii lokalnych aproksymuje energię całkowitą badanego stanu. Intrygowało mnie, które punkty są wybierane przez algorytm Metropolisa. Ze względu na fakt, że rozważany atom wodoru w stanie 1s jest kulistosymetryczny, wystarczyło zapytać, które wartości promienia wodzącego są wybierane i z jaką częstością.
W tym celu zmodyfikowałem program do obliczenia energii, umieszczony we Wprawce II. Usunąłem fragmenty obliczające energię lokalną i pozostawiłem tylko algorytm wyboru promienia. Dodałem procedurę zliczającą promienie. Zakres wektora liczności (6.2*100) ustaliłem eksperymentalnie dla 500000 losowań. Kod tak zmodyfikowanego programu ma następującą postać:

DECLARE FUNCTION fufal! (x!, y!, z!)
DECLARE FUNCTION GRNF! ()  'gaussowski generator pseudolosowy
                           'jednostki atomowe
CLS
SCREEN 12
llos = 500000           'liczba losowan
DIM SHARED c             'zmienna wariacyjna globalna
DIM x(620)               'wektor licznosci promieni
c = 1.5                  'wybrana wartosc zm. wariacyjnej
OPEN "dane.DAT" FOR APPEND AS #1

x = 0: y = 1.2: z = 0    'poczatkowy punkt
FOR i = 1 TO llos
  xzapas = x             'zabezpieczenie starych wspolrzednych
  yzapas = y
  zzapas = z

  psi = fufal(x, y, z)    'wartosc psi w punkcie (x,y,z)
  psikw = psi * psi       'kwadrat psi

           'Losowanie nowego punktu  (xn,yn,zn)
  xn = 2 * RND - 1        'zakres (-1, 1)
  yn = 2 * RND - 1
  zn = 2 * RND - 1

  pr = SQR(xn * xn + yn * yn + zn * zn)  'promien punktu (xn,yn,zn)
  skok = .6 * GRNF         'skok przypadkowy wg rozkladu Gaussa

  x = x + xn * skok / pr  'poprawione wspolrzedne nowego punktu
  y = y + yn * skok / pr
  z = z + zn * skok / pr
 
  psinowa = fufal(x, y, z)     'wartsc psi w nowym punkcie
  psinowakw = psinowa * psinowa

'obliczanie warunku przejscia lancucha do nowego punktu

 rel = psinowakw / psikw   'stosunek nowej i starej psi
  IF rel < 1 THEN          'przejscie bezwarunkowe gdy rel >= 1
    IF RND > rel THEN      'przejscie odrzucone
      x = xzapas           'przywrocenie starych wartosci
      y = yzapas
      z = zzapas
      psinowa = fufal(x, y, z)
    END IF
 END IF

  rn = SQR(x * x + y * y + z * z) 'promien do zliczania
  liczba = INT(rn * 100)     ' definicja adresu wektora promieni
  x(liczba) = x(liczba) + 1  ' zliczanie promieni
 NEXT i

FOR i = 1 TO 620
  PRINT #1, i, x(i)
NEXT i
CLOSE
END

FUNCTION fufal (x, y, z)
  r = SQR(x * x + y * y + z * z)
  fufal = EXP(-c * r)
END FUNCTION

FUNCTION GRNF
  r1 = SQR(-LOG(1 - RND))
  GRNF = r1 * SIN(2 * 3.14159265# * RND)
END FUNCTION

Program generuje zbiór danych, które można przedstawić graficznie. Poniżej zamieściłem wykres uzyskany za pomocą programu Excel:

Jak widać na wykresie, algorytm Metropolisa traktuje promienie w sposób wybiórczy. Uzyskana krzywa przypomina funkcję dystrybucji radialnej (radialnej gęstości elektronowej) stanu podstawowego atomu wodoru z maksimum poniżej jedności (w obliczeniach użyłem jednostki atomowe).