Решение уравнениякона Шема для цилиндрических атомов методом опорной функции Текст научной статьи по специальности «Физика»

МАТЕМАТИЧЕСКАЯ ФИЗИКА, МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

MSC 81V45

РЕШЕНИЕ УРАВНЕНИЯ КОНА ШЕМА ДЛЯ ЦИЛИНДРИЧЕСКИХ АТОМОВ МЕТОДОМ ОПОРНОЙ ФУНКЦИИ

А.В. Береговой, А.Г. Шкловский

Белгородский государственный университет, ул. Победы, 85, Белгород, 308015, Россия, e-mail: Shklovsky@bsu.edu.ru

Аннотация. Рассматривается решение уравнения Кона-Шема для цилиндрических атомов методом опорной функции на примере атомов углерода и кислорода.

Ключевые слова: уравнение Кона-Шема, цилиндрический атом, метод опорной функции, локальный обменно-корреляционный потенциал.

1. Введение. В настоящей работе рассматривается решение уравнения Кона-Шема для цилиндрических атомов методом опорной функции на примере атомов углерода и кислорода. Обычно для атомов используется предложенное Слэтером [1] приближение сферического атома. Это приближение предполагает, что вместо реальной электронной плотности следующей из формулы Кона-Шема [2], следует брать электронную плотность усреднённую по углам. Конечно для многих атомов, на пример атомов щелочных металлов или благородных газов это приближение является точным. Более того, если учесть, что в методе Кона-Шема не требуется чтобы числа заполнения являлись целыми, то всегда можно подобрать числа заполнения так, что бы получившаяся электронная плотность была сферически симметричной. Например для углерода можно взять число заполнения г орбитали равное 2/3 и такие же числа заполнения взять для х и у орбиталей. При этом электронная плотность окажется сферически симметричной и полное число заполнения 2р орбитали будет равно 2.

В случае использования для заполнения только целых чисел оба электрона в углероде должны занять тот уровень, энергия которого окажется ниже. В этом случае атом углерода будет иметь уже не сферическую электронную плотность, а цилиндрическую, так как соответствующую орбиталь всегда можно назвать г орбиталью.

Конечно, решение цилиндрической задачи требует гораздо больших вычислительных усилий и может быть оправдано только в том случае, если для обменнокорреляционного потенциала Кона-Шема есть достаточно точное приближение. В работе [3] предложено новое приближение для локального обменно-корреляционного потенциала релятивистских атомов. Естественно при этом использовалось сферическое приближение для электронной плотности. Для атомов углерода и кислорода релятивистские поправки малы, поэтому будем решать уравнение Кона-Шема в нерелятивистском приближении и подбирать константы в аппроксимирующем потенциале так,

чтобы энергии ионизации и энергии внутренних оболочек, рассчитанные с использованием приближённого потенциала, были по возможности близки к экспериментальным энергиям.

Далее будет использоваться атомная система единиц К = е2 = те = 1.

Расположим ядро атома углерода или кислорода в начале координат. Будем предполагать, что атомный потенциал V(г) не зависит от угла

V(r) = V(r, z).

Это позволяет свести трёхмерное уравнение Кона-Шема:

-у + V(r) ) Фп(г) = ЕпЧ>п(г)

(1)

(2)

к численному решению двумерного уравнения.

В цилиндрической системе координат уравнение (2), с учётом условия (1), можно переписать в виде:

1/1 д ( д\ 1 д2 д2 . ч . т ^ ^ т ^

“2 Ь?) + 7'W + а? 1 + v(г-1 ф”(г1 =

(З)

Будем искать решение уравнения (3) для функции Фп(г) методом разделения переменных:

фп(Г) = Хп,т(г, г)Лт(ф). (4)

Для Ат(^) легко получается уравнение:

д2Ат(ф 2 л г \

д(р2 = ~т Агп{ф).

Здесь Am(^) удовлетворяет очевидному граничному условию:

Am(^ + 2п) Am(^)-

(5)

(б)

Решение уравнения (5) с граничными условиями (6) очевидно и приводит к системе ортогональных функций:

Am(<^) = cos (m ■ ф), Am(^) = sin (m ■ ф).

(Т)

Подставляя (7) в (6), получаем, что для выполнения этих граничных условий m должно быть целым числом.

Подставим в уравнение (3) выражение (4) с учётом (7). После сокращения sin или cos уравнение для функции Xn,m{r, z) принимает вид:

4 (І¥г(гїії)-ф + їЬ} + УІГ-:)

Xn,m(r, z) — En,mXn,m(r, z) .

(8)

Оставшееся уравнение содержит только две переменные г и г. Для того чтобы избавиться в (8) от первых производных по г запишем:

( ^ т(г, л) ( ^

Хп,т(Г, z) = ------=--- . (9)

Из соотношения (9) видно, что при г = 0, ип,т(0, г) = 0 при любом г. Подставляя (9) в (8), для функции ип,т(г,г) получаем уравнение:

д2 д2 т2 — 0,25 \

4~ 2У(г, л) I и>п,т{г 1 ~0 ^Еп^тип^т(т, л). (Ю)

дг2 дг2 г2

Очевидно, что на машине мы не можем реализовать бесконечные отрезки. Поэтому вводятся числа гс и гс за пределами которых можно считать решения уравнения (10) равными нулю. Следовательно, будем решать уравнение (10) с нулевыми граничными условиями:

ип,т(г, гс) = 0, ип,т(гс, г) = 0, ип,т(г, -гс) = 0, ип,т(0, г) = 0. (11)

Если сразу решать уравнение (10) численно, то возникают две трудности. Во-первых, уравнение (10) является уравнением на собственные функции и собственные значения. Качественных двумерных численных алгоритмов для решения таких задач до последнего времени не было. Во-вторых, замена второй производной на её дискретный аналог может быть произведена с не очень высокой точностью. Это связано с резким увеличением количества узлов сетки при уменьшении шага, или с уменьшением порядка аппроксимации при использовании неравномерной сетки. Обе эти трудности и приводили ранее к тому, что вместо прямых численных алгоритмов использовались проекционные методы. Эти трудности позволяет преодолеть метод опорной функции. Поэтому далее именно этот метод будет использован для численного решения уравнения (10) .

2. Метод опорной функции. Рассмотрим вариант метода опорной функции, применяемый для решения уравнения Кона-Шема в цилиндрических атомах. Для этого представим функцию ип,т(г, г) в виде:

ип,т(г,г)= Уп(г,г)+ Уп(г,г). (12)

Здесь функция уп(г, г) описывает отличие более точной волновой функции ип,т(г, г) от опорной функции Уп(г, г). Явный вид опорных функций Уп(г, г) выбирается из следующих соображений. Во-первых, эта функция должна быть или аналитической, или задаваться эффективным численным алгоритмом, гарантирующим очень точное задание этой функции в нужных нам узлах решётки. Во-вторых, эта функция должна достаточно хорошо аппроксимировать функцию ип(г, г) чтобы можно было считать функцию уп(г, г) малой добавкой.

С учётом условия (7) для того случая, когда полная функция не зависит от угла ф, для уравнения (10) в качестве опорной функции будем брать линейную комбинацию:

Фп,о(г,г)= ^ < ■ Ук(г, г), (13)

1< к< 3

где для к = 1, 2

Для k = З имеем

п{г,;) = *(£Ь£ (14)

x

Уз{г,г) =-----------5-----• (15)

x2

Здесь введено обозначение х для расстояния от ядра до электрона:

X = л/г2 + Z2 . (16)

При этом ап - нормировочный коэффициент, а Кк(г) - решение радиального уравнения Кона-Шема для сжатого атома углерода или кислорода:

d2 / , ч l(l + 1)\

~d?+'2' (V^W + -2TP_)

Rk (x) = 2 ■ ЄkRk (x). (1Т)

Здесь k =1 соответствует ls-состоянию, k = 2 - 2s-состоянию, а k = 3 - 2р-состоянию сжатого атома углерода или кислорода. Потенциал сферически симметричного атома Vsph(x) был описан в работе [3].

Для случая, когда полная функция пропорциональна cos p или sin p имеем более

простую линейную комбинацию сжатых 2р атомных функций углерода или кислорода:

*t(r.=)=R,lx)\r''/f. (18)

x2

Видно, что решение для опорной функции (18) дважды вырождено, так как одна и та же опорная функция, а, следовательно, и одно и то же решение уравнения (10) получаются для различных полных функций. Первая полная функция будет отличаться на cos p, а вторая на sin p.

Перепишем уравнение Кона-Шема (10) для функции yn,m(r, z) в случае m = 0: д2 д2 О, 25

+ 2 • (V(г, z) - Еп<о) • у„,0(г, z) = Fn0(r, z). (19)

дт-2 дz2 r2

Здесь введены обозначения:

Fno(r, z) = 2 ^2 at ■ (En,o _ Єk _ Q(r, z)) ■ Yk(r, z). (20)

1<k<3

Q(r,z) = V(r,z) _ VsPh(x). (21)

Аналогично получим уравнение Кона-Шема для функции yn,1(r, z):

д2 д2 0 75 \

~~0^2 ~ ~qZ2 + 2 ' (^(г’ ~) “ Еп-1) ) ' 2/™д(г’ ~) = Fnl(r, z). (22)

Здесь введено обозначение:

Р*п1 (г, г) = 2 (Еп,1 - £3 - я(г, г)) ■ Ф4(г, г). (23)

В нулевом приближении можно положить уп(р, г) равным нулю, тогда описанные выше опорные функции соответствуют известному методу сильной связи. Интегралы перекрытия Ьтп и <!тп вычисляются в нулевом приближении обычным способом. Это позволяет найти коэффициенты ап в этом приближении.

Прежде чем обсуждать алгоритм решения уравнений (19) или (22) необходимо разобраться с рядом математических вопросов. Отметим, что само уравнение (19), получено из дифференциального уравнения (10) с помощью простых алгебраических преобразований. Поэтому, вопрос о существовании решения уравнения (19) (или (22)) не стоит, однако не очевидно, что любое решение уравнения (19) будет также решением уравнения (10). Может также оказаться, что одному решению уравнения (10) будет соответствовать несколько решений уравнения (19) или (22).

Пусть у нас существуют две функции д1 и д2, которые являются различными решениями уравнения (19), относящимися к одной и той же энергии Еп,т. Для простоты рассмотрим случай, когда в исходном уравнении (10) эта энергия является не вырожденной. Вычитая из уравнения (19) для функции д1 то же самое уравнение для функции д2 и приводя подобные члены, получим линейное уравнение:

д2 д2 ^т2 0 25 \

2У(г, л) I <7га,т(?”) ~0 2Еп тдп т(г, л). (-^4)

дг2 дг2 г2

Здесь дпт(г, г) = д1>п,т(г, г) _ д2>п,т(г, г). Так как для уравнения (24) Еп,т является не вырожденным собственным значением, то дп,т(г, г) либо равна нулю, либо пропорциональна пп>т(г, г):

дп,т(г, г) = Л ■ Пп,т(г, г),

д1,п,т(г, г) = Л ■ Пп,т(г, г) + д2,п,т(г, г). (25)

Так как нас интересует нормированное решение уравнения (10), то подставим (25) в (12):

|пп,т(г, г) | = |Уп(г, г) + Л ■ п,п,т(г, г) + дч,п,т(г, г) | = | (1 + Л) ■ п,п,т(г, г) | . (26)

Отсюда, используя условие нормировки для функции пп,т(г, г), получим уравнение:

(1 + Л)2 = 1, Л1 = 0, Л2 = _2 . (27)

Таким образом, мы видим, что в отличие от уравнения (10), уравнение (19) имеет два различных решения. Первое решение д1 действительно является малой добавкой к опорной функции. Второе решение д2, соответствующее Л2 = _2, приводит согласно (26) к функции —пп(г, г). Эта функция в уравнении (10) не давала нового решения, так как решение линейного однородного уравнения находится с точностью до произвольного коэффициента. Уравнение (19) не однородно и, очевидно, что функции д1 и д2 - это

разные решения. Поэтому при построении алгоритма следует позаботиться о том, чтобы решения уп,т(г, г) по сравнению с опорными функциями были малы. В этом случае второе решение окажется исключенным.

Перейдем к построению итерационной процедуры решения уравнений (19) и (22). Эти уравнения для различных т и п одного типа. Метод решения для них один и тот же. Все они удовлетворяют нулевым граничным условиям:

Уп,т(г, г с) = Уп,т(г, _гс) = Уп,т(г с, г) = Уп,т(0, г) = 0 . (28)

Будем решать методом последовательных приближений. Перейдем от дифференциальных уравнений (19) к разностным. Введем неравномерные сетки по осям г и г:

П1 = 1, 2,N1, П2 = 1, 2,..., N2, хг(пл) = гп1, гс(п) = гт,

хг(п1 + 1) = хг(п1) + Нг(п1), гс(п2 + 1) = гс(п2) + Нг(п2). (29)

Здесь массивы узлов гп1 с различными номерами связаны друг с другом соотношением (29). При этом хг(1) = 0 и хг(^) = гс Аналогично для массивов гп2 имеем гс(1) = _гс, гс(^) = гс. Так как мы собираемся находить функцию уп,т(г, г) приближенно с помощью итерационной процедуры и следить за тем чтобы она была мала по сравнению с пп(г, г) , то в качестве нулевого приближения следует положить уп,т(г, г) = 0 . Введем два массива ат(п1, п2) и Ьт(п1, п2) . Первый массив будет задавать сеточную функцию приближенно описывающую функцию уп,т(г, г) на нечетных итерациях, а второй - на четных. Тогда вместо уравнения (19) получим разностное уравнение:

Ьт(п1, п2) = ^Аг(п1, п2) ■ ат(п1 + 1, п2) + Вг(п1, п2) ■ ат(п1 _ 1, п2) +

Аг(п1, п2) ■ ат(п1, п2 + 1) + Вг(п1, п2) ■ ат(п1, п2 _ 1) + Ет(п1, п2)^ ■ Бт(п1, п2) (30) Здесь введены обозначения:

2 ■ кг(п2)

(Нг(п-]) + Нг(п1 _ 1))

, , hrz(nl,n2)

г{пъ п2) - + Нг(щ _ . Нг(щ _ .

Нгг(п1, п2) = 2 ■ кг(п2) ■ кг(п1). (31)

.4фь»,)- 2'ЛГ(Л1)

(Нг(п2) + кг(п2 _ 1)) :

\- hrz(n1,n2) . л

^/?1,/?2) (ЬЫ + М«2 - 1)) • М«2 - 1) ‘

г 0.125 1-1

Dm.ini, п2) = БгОпг, п2) + Dz(nl,n2) + (V(п1, п2) - Еп<0-------------т——гт) • Ьг^{;п1,п2)

I хг(п1 + 1)

Гт(п1,щ) = ЕП}о(г,П1 ,гт) ■ Нгг(п1,щ). (33)

В (33) введены обозначения:

Г)г(,н • "2) = Ипї’-ч ■ 0:(П1- "2) = Ьфі-І) ■ (34)

Система переменных узлов может быть выбрана так, чтобы обеспечить для рекуррентной формулы (30) относительную погрешность 6 < 10-3 . Так как мы строим алгоритм, для которого сама величина уП^1П2 мала по сравнению с иа(т,г) , то такая точность приемлема.

Уравнение (22) может быть записано в виде разностного аналогично тому, как разностное уравнение (30) описывало уравнение (19):

Ьтт(п1,п2) = (Ат(и\, п2) ■ атт(п1 + 1, п2) + Бт(иі, п2) ■ атт(п1 — 1, п2) +

+ Аг(п1, п2) ■ атт(п1, п2 + 1) + Бг(п1, п2) ■ атт(п1, п2 — 1) + Етт(п1, п2)) ■ Бтт(п1, п2).

(35)

Здесь введены обозначения:

г 0.375 і -і

Бішіпі, п2) = ЇМ/?і, п2) -\-Dzini, п2) + (V(щ, п2) -ЕпЛ Н-----—-----—) • hrz(nl,n2) ,

I хт(п1 + 1)2 1

Етт(пі, щ) = ЕП1 (тП1 ,гП2) ■ Нтг^щ) . (36)

Очевидно, что представление разностных уравнений в виде (30) и (35) удобно для применения итерационного метода решения. В (33) и (36) входит неизвестная энергия атома Еп,т = Еа . Для нахождения энергии Еп,т воспользуемся уравнением (19). Для этого умножим его слева на ип,т(т, г) и проинтегрируем по т и г . Учитывая, что уравнение (10) является самосопряжённым, и сокращая подобные члены, получим:

£ а“[(Е„ — е,) ■ Щ — д“] =0. (37)

3

Здесь введены обозначения:

%гп г т

,=

Щ = йг па(т,г)У,(т,г)йт, (38)

— Zm 0

О* = j dz j па(г,г) • О(г, г) • У^(г,г)йг. (39)

-%гп О

Система уравнений (37) является однородной и имеет ненулевое решение, если ее определитель равен нулю. Очевидно, что эту систему следует решать самосогласованным способом. В нулевом приближении па(г,г) совпадает с Ук(г,г) , таким, что энергия Еа наиболее близка к ек . При этом формула (38) дает интеграл перекрытия из приближения сильной связи. Если в формулу (39) подставить О (г, г) нулевого приближения, то

т

т

получится система уравнений сильной связи в первом приближении. Приравнивая определитель системы (37) с коэффициентами Щ и О°- нулевого приближения к нулю, можно получить энергии первого приближения Еа и коэффициенты аа нулевого приближения. Используя эти коэффициенты в формуле (13) можно начать процесс итерации, полагая в нулевом приближении уф П2 =0 Формула (37) позволяет начать итерационный процесс. Какое-то количество итераций по этой формуле можно сделать не вычисляя заново коэффициенты аа . На следующем этапе необходимо уточнить этот набор коэффициентов. Для того, чтобы это сделать следует заново вычислить коэффициенты (38) и (39), используя для па(г,г) формулы (12), (13)-(15) и (18).

Интегралы численно вычисляются гораздо точнее, чем производные. В оба интеграла входят функции ппт(г, г), которые с каждой итерацией находятся все точнее. Поэтому и вычисление интегралов в (38) и (39) тоже можно осуществить все более точно. Следовательно, осуществляя через какое-то количество итераций вычисление интегралов в (38) и (39), можно методом последовательных приближений с высокой точностью получить энергетический спектр цилиндрического атома.

3. Вычисление атомного потенциала. Вычисление атомного потенциала. В работе [4] был введен локальный модифицированный потенциал Слетера (МЛПС), с помощью которого удалось очень точно вычислить как полную энергию, так и энергию ионизации нерелятивистских атомов. В рамках теории функционала электронной плотности [5] в работе [4] было рассмотрено выражение для вычисления полной нерелятивистской энергии атома:

Е = Те[п] + ЕНаг[п] + Ехс[и] + J Уе(т) • n(r)dr. (40)

Для атомов углерода и кислорода, как и в работе [4], считаем, что уравнение Кона-Шема (2) описывает систему невзаимодействующих квазичастиц с энергией Еп, находящихся во внешнем потенциале У (г). Электронная плотность квазичастиц п(г) та же. что и у нерелятивистских электронов в атомах углерода или кислорода, поэтому для функционала Те[п] из (40) имеем выражение:

N „

Те[п] = ^ Еп — У (г)n(f)df. (41)

п=1 ^

Для удобства используем формулу (41) и перепишем выражение (40) для полной энергии Е в виде:

Ме I I

Е = ^2 Еп — У(г)n(r)dr + ЕНаг[п] + Ехс[п] + Уе(т)п(^т. (42)

п=1

Здесь энергия Хартри Е^аг [пс] дается обычным выражением:

= (43)

2] |г1 — г2|

В приближении МЛПС [4] функционал обменно-корреляционной энергии Ехс[пс] имеет вид:

где к = (3 + М(?)(Ме — 1) , Р(х) - аппроксимирующая функция, введённая в работе [3]:

Здесь коэффициент а2 определяет скорость убывания функции Г (х) вдали от ядра углерода или кислорода, и вычисляется, исходя из условия

которое обеспечивает выполнение закона сохранения перераспределённого заряда усреднённой обменно-корреляционной дырки.

Параметр а\ определяет значение х для которого Ра(х) = 0 , а константы а3 и а4 определяют количество дополнительных нулей в полиноме Ра(х) и подбираются так, чтобы численные значения полной энергии и энергии ионизации атомов углерода или кислорода совпали с экспериментальными данными. При этом остальные вычисленные энергии внутренних оболочек, должны быть как можно ближе к экспериментальным.

Теперь можно вернуться к вопросу о нахождении потенциала V(р, г) для цилиндрических атомов углерода и кислорода. Электронная плотность цилиндрического атома может быть представлена в виде:

Здесь Фп(г) является решением уравнения (2), а щ(?) показывает отличие реальной электронной плотности цилиндрического атома от электронной плотности сферически симметричного атома пс(|г|) и находится самосогласованным образом. При этом приближение для обменно-корреляционной энергии Ехс[п] не изменяется по сравнению с (44). Дело в том, что в эту формулу входит самосогласованная плотность п(г) обладающая уже не сферической, а цилиндрической симметрией.

Согласно теореме Хоэнберга и Кона [6] функционал полной энергии (42) для точной электронной плотности (48) имеет минимум. Так как минимум функционала (42)

/

4

пз (п)с1п ■ +

/

Г(х) = х ■ Ра(х)[а3 + а4Ра(х) + (1 — а3)Ра2(х)] ■ ехр(—а2х),

(45)

где

(46)

(47)

(48)

п

определяется из равенства нулю его первой вариации, то при очень малом изменении электронной плотности п(г) , например, за счёт малого изменения потенциала V(г) , поправки к полной энергии Е будут второго порядка малости. Поэтому вблизи от правильной электронной плотности модифицированный функционал можно считать точным вплоть до второго порядка малости. Следовательно, можно брать его вариацию по электронной плотности, считая параметры аппроксимации постоянными. Вычисляя вариацию функционала (46) по электронной плотности п(г) и приравнивая её к нулю, получим выражение для неизвестного потенциала цилиндрического атома V(г) :

V(г) = Ve(r) + Vhar(г) + ^(г). (49)

Здесь Ve(т) - потенциал притяжения электронов к ядру, выражение для потенциала Vhar (т) с учетом (48) имеет обычный вид:

УнаГ{г) = [ ^Г21. (1Г*2 = УНа,{\Л) + УМ{г) . (50)

3 \Г - Г'2 \

Сферическая часть потенциала Хартри У}1аг (г) даётся формулой:

г Ко

(\ Л (1 [ ( \ 2А , [ Пс(г2)г22(1г2\ Ме • (1 - дЛог(г)) . ,

”наг{г) = 4тг - пс{г2)г2аг2 + ------------- =--------------------, (51)

\г / г2 / г

Я

4п ■ 0(Яс — г) [ ,_л(л Г

Фшг(г) = ------—------- J X Пс(х){1 ~ -}(1х. (52)

г

Очевидно, что добавочный потенциал Хартри имеет вид:

Уы{г) = [ 1]^Г21 йг2 ■ (53)

3 \Г — Г21

Обменно-корреляционный потенциал псевдоатома Ухс{т) состоит из трёх частей:

+ [ + к • Щг1) ■ / . (54)

3 \Г — Гі\ ] \г — Гі\

Первый член в квадратных скобках - это потенциал Слетера, второй член компенсирует усредненное самодействие в потенциале Хартри, а третий член аппроксимирует остаточную часть обменно-корреляционной дырки.

В разделе 2 описано, как решение уравнения Кона-Шема методом опорной функции осуществляется с использованием итерационной процедуры. В начале итерационной процедуры нам не известна точная электронная плотность цилиндрического атома (48). Поэтому вместо потенциала V(г), вычисляемого по формуле (49), берется потенциал

нулевого приближения Уо(г), который находится следующим образом. В формулах (51)-(54) вместо электронной плотности п(т) используется электронная плотность нулевого приближения По (г), выраженная только через опорные функции:

п0(г,г) = 2^ \Уп(т,г)\2. (55)

п

Вычисление потенциала нулевого приближения V0(т) позволяет начать итерационный процесс, описанный в разделе 2. При нахождении всех фп(г) в первом приближении и вычислении в этом приближении поправки щ(г) можно одновременно найти поправку Vhd(r) к потенциалу Хартри (53). Заметим, что этот интеграл есть решение уравнения Пуассона:

△Ц*(г) = -4пщ(т). (56)

Так как щ(Г) не зависит от угла то будем численно решать уравнение Пуассона в цилиндрической системе координат. Перепишем уравнение (56) в виде:

+ = (57>

Введем новую функцию:

[/(г, л) = 1/ы(г, л) • у/г . (58)

Для неё получаем следующее уравнение:

д2 и д2и и

-^ + ^ + ^-2 = г). (59)

Для нахождения граничных условий в точках гс и гс вычислим Vhd в произвольной точке т за пределами атома углерода. В этой области щ(т) = 0. В формуле (53) разложим \т — т2\-1 в ряд по полиномам Лежандра для |г2| < |т|:

1 (60)

Здесь y - угол между r2 и г . Так как для основной части плотности щ{г2) имеем |r^21 ^ И ограничимся в этом разложении тремя первыми членами. Учтем также, что:

гт2 ХХ2 + УУ2 + ZZ2

cos 7=^— =----------------------. 61

ИИ | ^|| Г 21

Отсюда следует, что с точностью до членов |г2|2/ |г|3

1 1 , \Г2\ , |rl|2 3COS27 — 1

БТ + 7^2 cos 7 + Т^Гз"------о-----• (62)

|г -г2| |f| |f!2 ' |f!3 2

Подставим (60) в (53):

Vhd{f) = Щ J ш{г2)сЩ + Щз J nd(f?2)(xx2 + уу2 + zz2)df2+

+ |rl|2(3cos27-l)drl. (63)

Первый интеграл в (63) равен нулю, так как он численно равен заряду, создаваемому плотностью nd(r) • Рассмотрим второй интеграл:

J Щ(Т2)(ХХ2 + УУ2 + ZZ2)df*2 = xj Пй(г"2)ж2^Г"2 +

Гс Zc

+ lyJ шЫшЛъ + zj = 2nzj г**2! Z2) ■ ^ - о. («)

0 -Zc

Этот интеграл равен нулю в силу симметрии nd(r2,z2) = Пл(г2, -z2) . Таким образом получим:

Vhd(r) = -7j4^(3 cos2 7 - 1) • iid{f2)df2 ■ (65)

21 г |3

Подстановка (61) в (65) и переход в интеграле к цилиндрической системе координат, позволяет переписать (65) в виде:

Т/ ( ^ 7Г(Г’2 “ 2z^ Т ГмЛ

ум(г, z) = —---------у- • Ihd , (66)

2(Н + z2) 2

где

Г -

Vd(r2, z2){rl - 2z%)r2dr2dz2 (r2 + ,2)

J 4d\l 2, ~2 \l 2 - *~2 '2Ш 2U~2 (^

!hd= I I ------------------,_2 I ^---------------• (67)

-Zc 0

По формуле (67) мы можем найти потенциал вне атома углерода, а с ним и величину

и {гс, л) :

г ту ^ к(r■c-‘2z2)vftгc

и{гс^)= -1Ы. (68)

2(г2 + z2) 2

Аналогично вычисляется и величина и (г, гс). На линии г = 0 имеем условие и (0, г) = 0.

ди(г, г)

На линии л = 0 в силу симметрии имеем —-—к=0 = 0. Таким образом решение урав-

дг

нения (59) с указанными граничными условиями дает нам возможность вычислить величину и (г, г), а следовательно, и Ум(р, г) во всех внутренних точках молекулы. Алгоритм решения уравнения (59) аналогичен тому, который был использован при решении уравнения (19).

4. Вычисление полной энергии и энергий внутренних оболочек цилиндрического атома. С учётом (49)-(54) полная энергия цилиндрического атома может

быть записана в виде:

Подставим электронную плотность в виде (48) в выражение для энергии Хартри Енаг :

В (70) введено обозначение Е^рн для энергии Хартри сферического атома, приходящейся на один электрон и вычисляемой по формуле:

С другой стороны, используя (51) и (53), можно записать энергию Хартри в виде:

При решении уравнения Кона-Шема потенциалы Уарг (г~1) и Уы(г~1) вычисляются по формулам (74) и (53). К моменту окончания итерационной процедуры соответствующие двумерные массивы уже существуют в рабочей области программы и могут быть использованы для вычисления соответствующих двумерных интегралов в формуле (73).

5. Вычислительный эксперимент. Результаты апробации выведенных нами формул мы дадим на примере расчётов полных нерелятивистских энергий цилиндрических атомов углерода и кислорода. Дело в том, что для цилиндрического атома углерода мы должны разместить два р электрона либо на орбитали рх (с противоположными спинами), либо на рх и ру орбиталях. Если оба электрона расположены на рх орбитали, то рх и ру орбитали полностью свободны. Проведённый расчёт показал, что при этом энергия свободных орбиталей оказывается меньше чем энергия занятой орбитали. Если теперь переместить один электрон на рх а другой на ру орбиталь, то при таком заполнении электронной плотности после самосогласования энергия рх орбитали окажется лежащей ниже. Так как заполняться по теореме Кона-Шема должны орбитали с самой низкой энергией, то мы приходим к противоречию. Единственным способом разрешить это противоречие является отказ от целых чисел заполнения орбитали. Если положить что по две трети электрона заполняют рх,ру, и рх орбитали, то цилиндрический атом с таким заполнением сводится к сферическому атому, в котором все энергии р орбиталей совпадают.

Уи,ш(г, г с) = Уп,ш(г, —гс) = Уп,т(гс, г) = уп,т(0, г) = 0

(70)

Упт(г, гс) = Уп,ш(г, —гс) = Уп,ш(г'с, г) = Уп,ш(0, г) = 0

(71)

Уп,ш(г, гс) = Уп,ш(г, —гс) = Уп,ш(г'с, г) = Уп,ш(0, г) = 0 Подставим (72) в (69) и воспользуемся (44):

Уп,ш(г, гс) = Уп,ш(г, —гс) = Уп,ш(гс, г) = Уп,ш(0, г) = 0 . Здесь введено обозначение:

(72)

(73)

Уп,ш(г, г с) = Упш(г, —гс) = Упш^с, г) = Уп,ш(0, г) = 0 .

(74)

Совсем другая ситуация возникает для цилиндрического атома кислорода. В этом случае мы можем разместить на pz орбитали только два электрона с противоположно направленными спинами. Ещё по одному электрону можно разместить на px и py орбиталях, у которых энергии одинаковы и лежат выше энергии pz орбитали.

Литература

1. Шкловский А.Г., Шкловская М.А. Решение уравнения Кона-Шема для молекулы методом опорной функции // Научные ведомости БелГУ. Серия Физика. Математика. -2008. - №9(49);14. - C.123-136.

2. Старовойтов А.С., Шкловский А.Г. Модифицированный локальный потенциал для вычисления энергии ионизации атомов // Научные ведомости БелГУ. Серия Математи-ка.Физика. - 2010. - №11(82);19. - C.126-134.

3. Теория неоднородного электронного газа / Под ред. С. Лундквиста, Н. Марча, - М.: Мир, 1987. - 400 с

4. Шкловский А.Г. Аппроксимация обменно-корреляционного потенциала в методе функционала электронной плотности // Научные ведомости БелГУ. Серия Физикоматематические науки. - 2007. - №6(37);13. - С.150-155.

5. Бабичев А.П., Бабушкина И.А., Братковский А.М. и др. Физические величины / Справочник / М.: Энергоатомиздат, 1991. - 1232 с.

6. Rajagopal A.K. in Advances in Chemical Physics / ed. I. Prigogine, S.A. Rice, V.41 / New York: Wiley, 1980. - P.59.

7. Шкловский А.Г. Алгоритм численного решения уравнения Дирака в центрально симметричном поле // Электромагнитные волны и электронные системы. - 2004. - 9;12. -С.4-9.

8. Hohenberg P., Kohn W. Inhomogeneous Electron Gas // Phys. Rev. B. - 1964. - 136. -P.864-871.

9. Langreth D.C., Perdew J.P. // Phys. Rev. - 1977. - B15. - P.2884.

THE SOLUTION OF THE KOHN-SHEM EQUATION FOR CYLINDRICAL ATOMS BASED ON THE METHOD OF SUPPORT FUNCTION

A.V. Beregovoy, A.G. Shklovsky

Belgorod State University,

Pobedy St., 85, Belgorod, 308015, Russia, e-mail: Shklovsky@bsu.edu.ru

Abstract. It is proposed the method of numerical solution of the Kohn-Shem equation for cylindrical atoms. It is based on the method of support function and applied for the carbon and oxygen atoms.

Key words: Kohn-Shem’s equation, cylindrical atom, method of support function, local exchange-correlation potential.