Задание №4

Решение уравнения Шрёдингера для одноэлектронного атома в сферических координатах выглядит следующим образом: \begin{equation*} \psi_{n,l,m}(r, \theta, \phi) = \sqrt{\left( \frac{2}{na_0} \right)^3 \frac{(n - l - 1)!}{2n(n+l)!}} e^{-r/na_0} \left(\frac{2r}{na_0}\right)^l L_{n-l-1}^{2l+1}\left(\frac{2r}{na_0}\right) Y_{l,m}(\theta, \phi) \end{equation*}

Здесь $n, l$ и $m$ - квантовые числа,
$a_0$ - боровский радиус,
$L$ - обобщенные многочлены Лагерра,
$Y_{l,m}(\theta, \phi)$ - сферические гармоники.

In [1]:
import numpy
import npy2cube
import scipy.special
import scipy.misc
import npy2cube
import ipyvolume as ipv
from IPython.display import display

Процедура для вычисления значения волновой функции одноэлектронного атома в узлах некоторой сетки представлена ниже:

In [2]:
def wave_function(n,l,m,d):
    x,y,z = numpy.mgrid[-d:d:128j,-d:d:128j,-d:d:128j] 

    r = lambda x,y,z: numpy.sqrt(x**2+y**2+z**2) #Euclidean distance
    theta = lambda x,y,z: numpy.arccos(z/r(x,y,z)) # Compute polar angle
    phi = lambda x,y,z: numpy.arctan(y/x) #Compute azimuth

    a0 = 0.53

    R = lambda r,n,l: (((2. / n / a0) ** 3) * scipy.special.factorial(n - l - 1) / 2. / n / scipy.special.factorial(n + l)) ** 0.5 * numpy.exp(-r/n/a0) * ((2*r/n/a0) ** l) * scipy.special.genlaguerre(n-l-1,2*l+1)(2*r/n/a0) #Radial part of wave function
    WF = lambda r,theta,phi,n,l,m: R(r,n,l) * scipy.special.sph_harm(m,l,phi,theta) # wave function
    absWF = lambda r,theta,phi,n,l,m: numpy.absolute(WF(r,theta,phi,n,l,m))**2 # Square of wave function

    return absWF(r(x,y,z),theta(x,y,z),phi(x,y,z),n,l,m)

Нарисуем электронную плотностью для первых нескольких квантовых чисел с помощью пакета ipyvolume:

In [4]:
d = 15
for n in range(1,4):
    for l in range(0,n):
        for m in range(0,l+1):
            print 'n = ' + str(n) + ', l = ' + str(l) + ', m = ' + str(m)
            ipv.pylab.figure()
            ipv.pylab.volshow(wave_function(n, l, m, d))
            ipv.pylab.show()

n = 1, l = 0, m = 0

n = 2, l = 0, m = 0

n = 2, l = 1, m = 0

n = 2, l = 1, m = 1

n = 3, l = 0, m = 0

n = 3, l = 1, m = 0

n = 3, l = 1, m = 1

n = 3, l = 2, m = 0

n = 3, l = 2, m = 1

n = 3, l = 2, m = 2

Попробуем с помощью ORCA получить изображение тех же электронных плотностей. Мы не можем использовать для этих целей атом водорода, поскольку у него в невозбужденном состоянии занята лишь нижняя орбиталь и ORCA, соответственно откажется рисовать более высокие орбитали. Используем для наших целей атом натрия, у которого занято достаточное количество орбиталей для того, чтоб получить изображения электронных плотностей для первых нескольких квантовых чисел.

В качестве входного файла для ORCA используем следующее:

In [ ]:
! UHF SVP XYZFile
%plots Format Cube
 dim1 100
 dim2 100
 dim3 100
 MO("H-0.cube",0,0);
 MO("H-1.cube",1,0);
 MO("H-2.cube",2,0);
 MO("H-3.cube",3,0);
 MO("H-4.cube",4,0);
 MO("H-5.cube",5,0);
 MO("H-6.cube",6,0);
 MO("H-7.cube",7,0);
 MO("H-8.cube",8,0);
 MO("H-9.cube",9,0);
end
* xyz 0 2
 Na 0 0 0
*

Визуализируем полученные cube-файлы с помощью ipyvolume:

In [5]:
for i in xrange(10):
    with open('orca_densities/H-'+ str(i) +'.cube') as f:
        lines = f.readlines()[8:]
        grid = numpy.array(map(lambda x: abs(float(x)), ' '.join(lines).split())).reshape((100, 100, 100))
        ipv.pylab.figure()
        ipv.pylab.volshow(grid)
        ipv.pylab.show()

Видно, что полученные изображения практически в точности совпадают с изображениями, полученными непосредственной визуализацией решений уравнения Шрёдингера.