CC BY
44
МАТЕМАТИЧЕСКАЯ ФИЗИКА, МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
УДК 519.673, 539.182
ПРИБЛИЖЕНИЕ ЛОКАЛЬНОГО ФУНКЦИОНАЛА ПЛОТНОСТИ С ОБМЕННО-КОРРЕЛЯЦИОННОЙ ЭНЕРГИИ ДЛЯ РЕЛЯТИВИСТСКИХ АТОМОВ
А. В. Береговой, А.А. Плесканев, А.Г. Шкловский
Белгородский государственный университет, ул. Победы, 85, Белгород, 308015, Россия, e-mail: Shklovsky@bsu.edu.ru
Аннотация. Описан метод расчета полной энергии сферически симметричного релятивистского атома в приближении локального функционала электронной плотности с аппроксимированной обменно-корреляционной энергией. Проведено сравнение с экспериментом рассчитанной полной энергии для атомов от гелия до криптона.
Ключевые слова: функционал электронной плотности, обменно-корреляционный потенциал, уравнение Кона-Шема-Дирака, релятивистский атом, метод опорной функции.
Введение. В настоящее время для расчета свойств атомов, молекул и твердых тел широко применяется теория функционала электронной плотности. В данной работе дается обзор общей теории в основном с использованием приближения локального функционала для однокомпонентных не спин-поляризованных электронных систем в основном состоянии. Предпринимается попытка улучшить приближение локального потенциала путем объединения двух методов: метода обменного потенциала Слетера [1] и приближения аппроксимированного локального потенциала [2]. Преимущества этого подхода иллюстрируются на примере расчета полной энергии, энергии ионизации и рентгеновских спектров релятивистских атомов. Релятивистская теория функционала электронной плотности излагается в широко известном обзоре Раджагопала [3]. При работе с атомами, молекулами и твердыми телами эта теория может быть упрощена в связи с тем, что в адиабатическом приближении ядра атомов считаются покоящимися. Поэтому наше приближение для релятивистских атомов сводится к тому, что вместо уравнения Кона-Шема [4], используемого в нерелятивистской теории, мы будем применять уравнение Кона-Шема-Дирака (КШД):
Hat(f )ф(г ) = £jф](г )
(1)
где
Hat(f)
Vat(r ) + С2 -2c(sp ) -2c(sp ) Vat (Г) - c2
Здесь и далее используется атомная система единиц К = е2 скорость света, компоненты векторного оператора спина 8-
Гх
1 0 0 1
Гу
0 -i i 0
(z
= rne = Г/2:
10 01
(2)
1, c = 137,03598
(3)
Четырех-компонентный биспинор фг(г) представляется в виде
Щ(г)
ф (Г) = (^), (4)
где Uj (г) и ю? (г) - обычные двух-компонентные спиноры.
Уравнение КШД описывает систему невзаимодействующих друг с другом квазичастиц с энергией ej, которую будем отсчитывать от энергии покоящегося электрона ej = (£j — с2) , находящихся во внешнем потенциале Уаг(т ). Потенциал Уаг(т ) выбирается так, чтобы электронная плотность квазичастиц п(г) :
п(г) = X] ф+(г )ф(г) (5)
j=1
была та же, что и у электронов в атоме. Поэтому для нее имеет место условие нормировки
J п(г = N (6)
и функционал кинетической энергии для этих квазичастиц Те [п] имеет вид:
N
п] = У~] ej — I Усл(т )п(г )<1г. (7)
Те[п] = X е? — ^(г )п(г j=l ^
Рассмотрим систему релятивистских электронов во внешнем поле и запишем их гамильтониан в виде:
Н = Т + Н', (8)
где Т оператор кинетической энергии, а Н' содержит взаимодействие с полем ядра и кулоновское взаимодействие электронов между собой. Рассмотрим эту систему для случая, когда квадрат заряда электрона сначала равен нулю, а затем, адиабатически возрастает до нормальной величины (0 < Л < 1):
Й'х = Ч ^Щг
2 у — г'|
п(г') — 5(г — г') + ¿гу\(г )п(г ). (9)
Здесь оператор плотности п(г ) задается суммой по спиновым индексам а произведений операторов рождения и уничтожения релятивистских электронов:
3
(г ) ф+(г )фа(г ). (10)
п(Г ) = , )Га
а=0
Для нерелятивистских электронов эту задачу рассматривали Лангрет и Педью [5]. При этом наблюдаемая электронная плотность была функцией трех переменных и определялась выражением
п(г) = (дп^|п(г )|дпд), (11)
где Ignd) - вектор основного состояния, а скалярное произведение этих векторов, как обычно, означает усреднение по этому состоянию. Как показано Раджагопалом [3] и для релятивистских электронов основное состояние является однозначным функционалом от внешнего потенциала v(r), а следовательно и наблюдаемая электронная плотность также является однозначным функционалом от v(r). Поэтому, как и в нерелятивистском случае справедливы обе теоремы Кона-Хоэнберга [6] и можно утверждать, что минимальное значение функционала энергии E[n] есть точная энергия основного состояния для внешнего поля v(r) и числа частиц Ne.
Для атома
v(r) = —
Z
(12)
где Z - заряд ядра. Тогда, возвращаясь к задаче (9) будем считать, что внешний потенциал v(r) заменен на потенциал v\(r) так, чтобы наблюдаемая плотность в основном состоянии n(r) не зависела от А и поэтому v\(r) = Vat(r) при А = 0, а при А = 1 v\(r) совпадает с v(r). Энергию основного состояния, вычисленную для любого промежуточного значения А, обозначим E\. В соответствии с теоремой Файнмана можно записать
dEx d\
gnd
дИ[
дА
gnd
drdr' |r — r'1
gnd
П(Г )[П(г') — ö(r — r')] д Г
gnd) +
+ ~д\ I (13)
После интегрирования по А получаем:
E = {gndH Ignd) = Ei = Eq +
dE> ~dX
d\.
(14)
Значение энергии основного состояния Е0 при А = 0 представляет собой энергию не взаимодействующих квазичастиц, находящихся во внешнем поле Уаг(?), и легко определяется из решения уравнения КШД (1):
Ne
e0 = X] ej • j=1
(15)
При этом кинетическая энергия этих квазичастиц Те[п] дается выражением (7). Эта кинетическая энергия возникает при вычислении интеграла в (14) от последнего слагаемого из (13):
Е = Е0+ [ (v(r) — Vatif)) n(r)dr + 7- [ drdr'+ Ехс[п].
J * J | r — r' I
r
X
1
В (16) введено обозначение Exc[n] для обменно-корреляционной энергии:
1
1
(17)
Как и в работе Лангрета и Педью [5], неизвестная функция к(т,г') может быть выражена через парную корреляционную функцию, вычисление которой может быть осуществлено только для простейших систем.
Очевидно, что сама попытка применения схемы КШД, как инструмента для решения многочастичного уравнения, имеет смысл, только если возможны явные приближения для функционала обменно-корреляционной энергии Ехс. Качество метода функционалов плотности зависит исключительно от точности выбранного приближения к Ехс. Поэтому, поиск все лучших и лучших функционалов составляет самую суть теории функционалов плотности.
Еще более десяти лет назад в обзоре [7] отмечалось: «Не только физика, лежащая в основе успеха существующих функционалов, далека от ясности, но мы просто не имеем ни малейшей идеи о том, как искать приближения, более близкие к точному функционалу. Все поиски лучших функционалов опираются на физическую или математическую интуицию и в значительной степени представляют из себя метод проб и ошибок».
Отметим, что, в отличие от традиционных методов основанных на волновой функции, в теории функционалов плотности нет систематического подхода к получению улучшенных приближенных функционалов, что фактически является одним из главных недостатков данной теории.
Существуют, однако, определенные физические ограничения, которым правдоподобный функционал должен соответствовать. Среди них заметное место принадлежит правилам сумм, верным для точных обменно-корреляционных дырок. Конечно, чем ближе модель дырки, возникающая из приближенного обменно-корреляционного функционала, отражает реальную дырку, тем лучше этот функционал будет способен учитывать неклассические эффекты. Другие свойства включают в себя условие отсутствия корреляции электронной дырки при больших расстояниях между двумя электронами, определенные условия масштабирования обменно-корреляционной энергии и асимптотические свойства соответствующих обменно-корреляционных потенциалов. Однако, мы не должны ожидать большой помощи от таких формальных граничных условий, т.к. одна из странных особенностей теории приближенного функционала плотности в том, что функционалы, жестко удовлетворяющие этим условиям, не обязательно лучше других. Фактически, некоторые из самых удачных приближенных функционалов нарушают некоторые из этих условий. Поэтому очень важно внимательно изучить поведение конкретного функционала относительно подходящих контрольных данных. Конечно, наиболее убедительные тесты, доступные сейчас для новых функционалов, полностью эмпирические и включают в себя сравнение с верными контрольными данными, таки-
ми как энергии диссоциации, ионизации и рекомбинации, структурные данные и тому подобное. Наиболее часто используемый набор контрольных данных энергий - это, возможно, так называемая G2 термодинамическая база данных, которая содержит более 50 экспериментально хорошо установленных энергий диссоциации маленьких молекул, содержащих элементы главной группы, первоначально собранная Кертисом и др. в 1991 [8]. Способность воспроизводить значения энергий этой базы или ее расширений де-факто стало стандартом для измерения точности нового численного метода. Желаемая точность - это так называемая химическая точность, которая соответствует средней абсолютной погрешности около 0.1 эВ или 2 ккал/моль.
Однако это очень амбициозная цель и до сих пор только очень немногие и очень ресурсоемкие традиционные квантовые химические стратегии были способны достичь такой точности. Поведение имевшихся до 2000 г. приближенных функционалов относительно G2 и соответствующие контрольные наборы детально рассмотрены в гл. 9 обзора [7]. В настоящей работе мы предлагаем другой подход к поиску приближенного выражения для обменно корреляционной энергии.
1. Аппроксимация обменно-корреляционной энергии. Рассмотрим текущее состояние дел относительно приближенных функционалов для Exc. Наиболее широко используемым приближением для Exc является приближение локальной плотности (ПЛП) [9]
Exc[n] = j £xc(n(r ))n(r )dr, (18)
где £xc(n) - обменно-корреляционная энергия, приходящаяся на одну частицу. Были построены различные аналитические аппроксимации для £xc(n), которые обеспечивали точность порядка 1 — 2%. Так Гуннарсон и Лундквист [10] предложили выражение
0 458 / г \
ехс(п) = - 0,0666G Ш , (19)
где гs - радиус Вигнера-Зейтца, определяемый формулой
4П
Yrs = n~1(r), (20)
а функция G(x) имеет вид
ад = \
(1 + 1п(1 + аГ1) - + | - ^
2 3
(21)
При этом обменно корреляционный потенциал vxc(n) имеет вид
в
Ухс(п) = —[п£хс(п)] . (22)
Вычисление полной энергии положительного иона натрия по формуле (17) дает в этом приближении -4398,13 эВ, а экспериментально полученное значение -4414, 81 эВ. Абсолютная погрешность составляет 16, 68 эВ, а относительная погрешность всего 0, 38%.
Такой подход привлекает своей простотой и универсальностью. Если не требуется более высокой точности при вычислении полной энергии этот подход является вполне удовлетворительным.
Следует отметить,что для элементов с большим зарядом ядра полная энергия может составлять десятки и даже сотни тысяч электрон вольт, поэтому погрешность даже менее 1% будет давать ошибку в десятки и даже сотни эВ. Для вычисления энергии ионизации необходимо из вычисленной полной энергии атома отнять вычисленную энергию иона. Поскольку при вычислении этих энергий, они имеют погрешность в десятки электрон вольт, то полученное значение энергии ионизации вызывает большие сомнения, хотя иногда и оказывается близким к экспериментальному значению. Поэтому считается, что ранние варианты ПЛП [10-13] недостаточно точны для серьезных расчетов. Для улучшения получаемых результатов обычно пытаются учесть градиентные поправки [7,14].
В настоящей работе мы рассмотрим другой подход, направленный на улучшение ПЛП для релятивистских атомов. Он связан с непосредственным использованием представления обменно-корреляционной энергии в виде (17).
Напомним, что, согласно второй теореме Хоэнберга-Кона [6], для правильной электронной плотности функционал полной энергии (16) имеет минимум. Таким образом, вычисляя вариацию функционала (16) по электронной плотности и приравнивая ее к нулю, получим выражение для неизвестного потенциала Vat:
Vat(r ) = v(r ) + Vhar [n] + VXc[n], (23)
где Vhar - потенциал Хартри
Vhar N = [ XdP, (24)
J |r — r' I
а Vxc - обменно-корреляционный потенциал
VWn] = ■ (25)
Так как минимум функционала определяется из равенства нулю его первой вариации, то при очень малом изменении электронной плотности, например, за счет малого изменения потенциала Vat, поправки к полной энергии E будут второго порядка малости. Поэтому возникает идея аппроксимировать функцию h в (17) с помощью некой универсальной функции, содержащей несколько констант аппроксимации. Параметры аппроксимации подбираются из условия равенства полной экспериментальной энергии для атома, значению функционала полной энергии (16) в соответствующем приближении, полученному для электронной плотности, которая определяется самосогласованным образом. Поэтому, вблизи от правильной электронной плотности, приближенный функционал можно считать точным вплоть до второго порядка малости. Следовательно, можно брать его вариацию по электронной плотности, считая параметры аппроксимации функции h постоянными.
Будем искать приближенное выражение для функции к, исходя из следующих соображений: она должна, во-первых, быть симметричной
к (г, г') = к(г',г), (26)
во-вторых, на достаточно больших расстояниях от ядра обращаться в ноль, и в-третьих, удовлетворять правилу сумм:
У п(Г )к(г, г' )<г = —1. (27)
Как показано в работе [15], при рассмотрении ограниченных систем (атомов или молекул), из функции к можно выделить постоянный член, описывающий компенсацию усредненного самодействия, которое входит в потенциал Хартри (24). Затем в работе [2] было показано, что в функции к можно выделить член не зависящий от электронной плотности и удовлетворяющий всем указанным выше условиям:
Ь.(г, ?) = _1 + №-1)((3 + №У(г-)Пг') + ^.г'))е(д _ |г- |)е{д _ |Р|) (28)
1Уе
"в
Здесь Я - радиус атома, в - функция Хэвисайда
е^ЦХ»0. <29>
Условие (27) накладывает ограничения на функции к1 (г) и п(г,г'):
У п(г (г )<г = 0, (30)
У па(г,г' )<г' = 0, (31)
где п^(г, г') = п(г')п(г, г') - плотность остаточной обменно-корреляционной дырки.
Явный вид аппроксимирующей функции к (г) будет обсуждаться ниже. Простейшее выражение для нее, содержащие только один независимый параметр, было предложено в [16].
Перепишем формулу (14) с учетом (15), (16) и (27)
] Мв I 2 У |г — г' |
3=1
I + ^ _ (32)
= -\( '^0-Лг'п(г)с1г, (33)
2 7 |г — г' |
Если считать в выражении для функционала полной энергии (32) функцию Г (г) известной, то неопределенным остается только функционал остаточной обменно-корре-ляционной энергии Ехсз (33). Формулу (33) можно переписать в виде
Ex
n(r )£xcs(r )dr,
(34)
где
s(r-) = -- [ Ud^rXir>
2 J |r — r'|
(35)
- энергия остаточной обменно-корреляционной дырки в точке r.
Рассмотрим простую модель, позволяющую определить функционал (33) с точностью до константы. Функция F(r ) дает усредненное описание обменно-корреляционной дырки для всей плотности n(r'). Поэтому nd(r, r') описывает отличие от этого среднего именно в малой окрестности точки r. Для выполнения условия (31) функция r') вне малой окрестности точки r должна быть быстро осциллирующей и малой по модулю. Следовательно, в качестве простейшей модели может быть принята функция
nd(r, ?) = C(r)n(r)@(Rd(r) - [f> - f I) = - \? - г I). (36)
4nRd(r)
Здесь qd(r) - заряд остаточной дырки, сосредоточенной внутри сферы малого радиуса Rd(r). Остальной частью функции n(r,r') пренебрегаем, т.к. она вносит малый вклад в интеграл (35). Из (36) следует, что
Rd(r) =
" 3 Qd(r) '
4п С(г) \
1/3
n
-1/3
(r).
(37)
Знак модуля в (37) учитывает, тот факт, что хотя в силу (31) функции qd(r) и С (г) могут принимать как положительные, так и отрицательные значения, их отношение должно быть всегда положительно.
Интегрирование в (35) с учетом (36) сводится к интегрированию по сфере малого радиуса ) с центром в точке г, в результате которого получаем
£xcs = —nC(r )n(r )Rd(r) = — C(r ) Подставляя (38) в (34), получаем
2
qd(r)
9п (qd(r)
16 V C(r)
n 3 (r )
(38)
Ex
9тг\ 3
ley
C (r)
C (r)
4 3/4
n^{f)dr— ——/3 I nz{r)dr.
(39)
Так как электронная плотность в атоме очень быстро возрастает по мере приближения к ядру, то наибольший вклад в интеграл (39) вносит именно эта область. В этой
£
2
области функция C{f)zqd(f)z является медленно меняющейся, фактически, оставаясь постоянной. Это и учтено в выражении (39). Оценим величину константы в:
*=(Тj1—гт——■ <40>
\ / I пз(г)с1г
Из (36) следует, что функции C(r ) и qd(r ) много меньше единицы. Возьмем для грубой оценки сверху C(r ) = 1/3 и qd(r ) = 1/3, тогда получаем, что в < 0.54. Обычно значение константы |в| лежит в интервале 0.1 — 0.4.
Теперь все члены в выражении для полной энергии (32) определены. 3. Решение уравнения Кона-Шема-Дирака. Вычислим первую вариацию функционала (32) и приравняем ее к нулю. Из этого уравнения можно определить неизвестный потенциал Vat, входящий в гамильтониан КШД (2):
Vat(r ) = -- + \ [ -^X-d? - (3 + N¡)F(r) [ Hr>ír'\l? - ¡3nh {r) } . (41)
r Ne \J |r - r'| J |r - r'|
Чтобы решить уравнение КШД (1) сделаем два дополнительных приближения. Во-первых, будем считать, что электронная плотность, входящая в Vat, сферически симметрична. Для этого электронную плотность, вычисленную по формуле (5), следует усреднить по углам в и Во-вторых, функцию F(r) тоже усредняем по углам и выбираем ее в виде
F(r) = rPa(r) [аз + a4Pa(r) + (1 — a3)Pa2(r)] exp(-a2r), (42)
где
Mr) = • (43)
r + 0.5
Здесь константа ai определяет значение r, для которого Pa(r) = 0. Следовательно для этого r в нуль обращается и функция F(r). В зависимости от коэффициентов a3 и a4, в полиноме третьего порядка от Pa(r), таких внутренних точек, в которых F(r) обращается в нуль, может быть одна, две или три. Кроме них, F(r) = 0 при r = 0 и r — то. Коэффициент a2, определяющий скорость убывания F(r) на бесконечности, не является свободным и находится из условия (30).
В силу введенных выше предположений, потенциал Vhar (r) в (24) легко вычисляется:
r i
ЛГ ( \ О Г n(r2)r^dr2 sin e2de2 f . л 2 , f dx
Vhar(r) = 2тг / === = 2тг / n(r2)r2dr2 / = , (44)
j \jr2 + r2 — 2rr2 cos в2 j j \jr2 + r2 — 2rr2x
2 0 -1 v 2
В (44) сделана замена переменной интегрирования cos в2 = x . Далее, используем формулу
1
f dx f 2/r, r2 < r ,
,_- Л 0 , (45)
1 л/г2 + r\ - 2rr2x l 2/r'2, Г <r2.
Подставляя(45) в (44), окончательно получаем:
N • (1 - днат (г))
Укат (г) = 4п
1 Г 2 Г п(г2)г^г2
~ / п(Г2)Г2аГ2 + / -
г У J Г 2
. 0 т
(46)
В формуле (46) введено обозначение дкат для среднего заряда Хартри, приходящегося на один электрон в атоме:
я
4п0(Я - г) Г 2 . . / г
фтг(г) =-—-j X п(х) (^1 --)йх. (47)
т
Аналогично может быть вычислен усредненный сферически симметричный потенциал обменно-корреляционной дырки. При этом формула (41) принимает вид:
К*(г) = -- + ^ • (1 - дНаЛт) + дхср(г) + дхс8(г)) = . (48)
г г г
Здесь введено обозначение Zeff (г) для усредненного сферически симметричного эффективного заряда атома:
Zeff (г) = Z - (Ме - 1) • (1 - днат(г) + дхср(г) + дХС8(г)). (49)
Усредненный обменно-корреляционный заряд дхср , приходящийся на один электрон в сферически симметричном атоме, задается формулой:
1 /3 Яа£
Ягср(г) = 47Г'(3 + ^ I .гПх)Щх) (1 - £) ¿г. (50)
Остаточный обменно-корреляционный заряд Слетера дхсз(г), приходящийся на одну частицу, имеет вид:
яиг) = ~ ■ (»(г))1/3- (51)
Формула (49) для Zeff (г) дает правильную асимптотику в нуле и на бесконечности, а для водородоподобных атомов (Ме = 1) Zeff = Z.
Таким образом, уравнение КШД (1) свелось к решению задачи о релятивистской частице, помещенной в сферически симметричное поле Уаг(г). В стандартном представлении (2) уравнение КШД расщепляется на пару двухкомпонетных уравнений. Гамильтониан КШД Я„, определенный в (2), коммутирует с операторами ,72 = (£ +§) и
= ^+ . Пусть двухкомпонентные спиноры и^ и и^ ^ являются собственными спинорами операторов 12 и 1 [17]:
/+1/2+™^ уШ^-1/2 1 1
при ]=1 + 1/2 о« , = ¡Кг) I I • «
1+1/2-т} —1/2 1 1
■ч»' =1 - ^ «К, = ^ ) ■ <53)
Здесь 1 - собственное значение оператора Ь2, т^ - собственное значение оператора , ут] 1/2 и +1/2 - сферические функции от углов 9 и р. Произвольные функции д(г) и /(г) должны быть нормированы.
Тогда собственными векторами операторов квадрата момента и проекции момента являются являются также биспиноры ф , определенные соотношениями
ф = (£!) и ф = (£5) (54)
Подставляя эти биспиноры в уравнение (1), получим две системы радиальных уравнений КШД:
при,) = / + 1/2 < Г # (55)
/--/ - г-д = 0,
г с
, , 3 + 3/2 .е - К; +
2
, +-9 ~ 1-:-/ =
при) = 1-1/2 < Г (56)
/--/ - г-£ = 0.
г с
Если предположить, что потенциал Уаг(г) = —Zeff (г)/г - заданная функция, то системы уравнений (55)-(56) представляют собой системы радиальных уравнений Дирака, для которых каждому собственному значению е отвечает пара собственных функций /(г) и д(г). Существование и единственность собственных функций и собственных значений для такого потенциала гарантировано известными теоремами.
При решении задачи о релятивистском атоме нас интересует нахождение тех значений дискретного спектра, которые относятся к медленно меняющимся функциям /(г) и д(г). Вводя обозначения д = С, / = гФ и е = с2 + е, для системы (56) получим:
(57)
Будем искать решение системы уравнений (57) методом последовательных приближений. Введем обозначения:
г
с
г
с
где е(0) - заданное число, близкое к собственному значению е. Из (57) выразим приближенно функцию Ф через функцию О и ее производную С:
Ф(г) = (а + ¿-^С. ) . (59)
сьцт) \ т
Найдем производную Ф ':
™ = ^ - ^ ■ (»)
Отметим, что чем ближе е(0) к е, тем точнее формулы (59) и (60).
Подставив (59) и (60) во второе уравнение системы (57), вводя обозначение О (г) = у(т)/т и учитывая, что ] = I — 1/2, получим
Ь(-)у(т) = е(-)у(г), (61)
= + + («о
Здесь введены обозначения е(-) для собственного значения линейного оператора Ь(-), наиболее близкого к е(0) и Р1{т)
Уд^т) 2с2Ь1{г) '
РЛг) = . (63)
Вводя обозначение Ф = г(т)/т, из уравнения (59) имеем
2{г) = ^ ^у'(г) + 1-у{г)^ . (64)
Таким образом, решая уравнение (61), находим у(т), а затем по формуле (64) определяем г(т). Из (64) видно, что г(т) - малая компонента релятивистской волновой функции, а у(т) - большая.
Придерживаясь такого же соглашения для обозначения малой и большой компонент релятивистской волновой функции, для случая ] = I + 1/2 из системы уравнений (55) получаем
Ь(+)у(т) = е(+)у(т), (65)
где
= - Л(г,1 + + ,,, + ^±11 _ М). (66)
z(r) = ^^у (у1 (г) - ЦгУ(г)^ • (67)
Отметим, что хотя уравнения (61) и (65) получены из систем уравнений Дирака (55)-(56) для радиальных функций, они отличаются от них тем, что относятся только к одному собственному значению. Для нахождения другого собственного значения и собственной функции нужно заново решать систему с другим значением е0. Для нахождения собственного значения применяется итерационная процедура, на входе которой задается число е0, а на первом шаге итерации получаем ближайшее к нему собственное значение е(-) или е(+) в зависимости от решаемой системы.
На следующем шаге, если вместо е0 просто подставить полученное собственное значение, то применение такой итерационной процедуры оказывается не эффективным. Это связано с тем, что при прямой итерации можно выйти за пределы области сходимости к данному собственному значению. В работе [18] описан простой алгоритм приближенного решения уравнения Дирака с заданным потенциалом. При реализации итерационной процедуры в [18] применен алгоритм Вегстейна [19], который позволяет при попадании е0 в область сходимости, получить собственное значение с точностью лучше чем 10-5 за 5-6 итераций. В системах компьютерной математики, например в МАТЬАБ, предусмотрены стандартные функции для нахождении собственных векторов и собственных значений, ближайших к заданному числу. Эти функции дают правильный результат, если матрицы, для которых ищутся собственные вектора и собственные значения, являются практически симметричными, т.е. их соответствующие верхние и нижние диагонали незначительно отличаются друг от друга.
При численной реализации второй производной на равномерной сетке получается симметричная трех-диагональная матрица. При численной реализации первой производной на равномерной сетке получается антисимметричная трех диагональная матрица. Заметим, что в (62) и (66) первая производная входит с малым коэффициентом, поскольку Р\ содержит с2 в знаменателе. При использовании неравномерной сетки вторая производная дает сильно несимметричную матрицу и стандартные алгоритмы перестают работать. Следовательно, недостатком примененного в [18] алгоритма, является необходимость использования равномерной сетки. Это ограничивает точность решения уравнений (61), (63) примерно четырьмя, пятью значащими цифрами. Для релятивистских атомов значения энергий 1з состояний может достигать десятков тысяч эВ. Чтобы получить при этом точность решения в одну сотую эВ нужно использовать более точный метод.
В приведенных выше рассуждениях потенциал Уаг(т) предполагался известным. Заметим, что в системе уравнений КШД потенциал Уаг(т) находится самосогласованным образом. Следовательно, необходима еще одна итерационная процедура, связанная с вычислением самого потенциала (41). В работах [20] и [21] был предложен метод опорной функции, который позволил увеличить точность расчета энергий молекул. Этот метод может быть применен и для уточнения решений уравнений КШД (61) и (64), полученных методом, рассмотренным в работе [18].
Задача ставится следующим образом. Пусть нам известен некий потенциал
7(°) (т)
= (68)
Функцию (г) зададим в виде сеточной функции у.(0) на системе узлов Xк, где к = 1, 2, . . . , К. Размерность этого массива К можно взять примерно 200 - 300. Заметим, что у(0) = Z, при Х1 = 0, где Z - заряд ядра. С помощью сплайн-интерполяции эта сеточная функция может быть пересчитана на любую систему равномерных узлов, число которых N1 может быть велико.
Производим достаточно точное решение систем уравнений (61), (64) и (65), (67) с потенциалом (68), применяя алгоритм работы [18]. В результате, мы имеем набор приближенных собственных значений е и соответствующие им массивы приближенных собственных векторов уП1, п1 = 1, 2,..., N1.
Пусть главное квантовое число п указывает номер оболочки, для которой получено указанное приближенное решение. Проведем классификацию решений по значению углового момента. Для в состояний (I = 0) мы имеем только систему уравнений (65)-(67)
и при этом ] = —, а т^ = =(=—. Следовательно, как и в нерелятивистском случае, в в
состоянии могут быть размещены два электрона.
В р состояниях (/ = 1) следует решать обе системы уравнений с j = — и т^ = ;
• 3 13 22
и с ;] = — и = Т^Т^- Поэтому шесть нерелятивистских электронов разделяются
на две подгруппы из двух (3 = 1/2) и четырех электронов (3 = 3/2). Эти подгруппы будем обозначать р1/2 и р3/2. Аналогичные обозначения ¿3/2 и ¿5/2 относятся к двум
подгруппам с / = 2 и з = — , з = — соответственно. Эта система обозначений позволяет классифицировать энергии связи электронов из внутренних оболочек основного состояния релятивистских атомов. Таким образом, запись 2р3/2 обозначает состояние электрона с п = 2,1 = 1,3 = 3/2, т^.
Пронумеруем все занятые состояния (п,1,з,т^) индексом который пробегает значения от 1 до Ne. Тогда, для электронной плотности (5) можем записать:
4пг2п(г) = ^(у|(г) + (г)). (69)
?=1
Здесь п(г) - радиальная электронная плотность релятивистских электронов в атоме. Подставляя ее в (49) можно вычислить следующее приближение для функции Zeff (г) и создать новый массив у.(1). Повторяя эту процедуру несколько раз, можно построить достаточно хорошие приближения для Zeff (г) и для функций У (г) = у% (г) и (г).
Полученная точность порядка 10-4 позволяет нам пренебречь разницей между приближенной энергией и точной в формулах (58),(63) т.к. в этих формулах погрешность будет делиться на величину с2 и следовательно будет менее 10-8. Для нас это очень мало, так как мы рассматриваем алгоритм, который должен дать решение с точностью до 10-7. Те же самые рассуждения касаются и разницы между функциями Zeff (г) и (г), входящими в формулы (58),(63). Поэтому главным отличием приближенного оператора Ь(0) от оператора Ь(-) является потенциал Уд(г):
Ь(-) = Ь(0) + Уд(г), (70)
где
(т) — (т))
У(г) = - т '-. (71)
т
Итак, будем считать, что у нас есть точное решение задачи:
Ь(0)У(—) (т) = е(0)У(—) (т). (72)
Требуется найти функцию у(т), являющуюся приближенной с точностью до 10-7 собственной функцией оператора (62). При этом собственное значение е(-) должно быть получено с такой же точностью. Будем искать у(т) в виде:
у(т) = У(-)(т)+ ур(т). (73)
Здесь опорная функция У(-)(т) является точной собственной функцией оператора 1(0 с точным собственным значением е(0). Малая функция ур(т) является решением уравнения:
£(-)уР(т) = е(-)ур(т) + РР (-)(т), (74)
где
рр (-)(т) = [(е(-) — е(0)) — К (т)]У (-)(т). (75)
Прежде чем обсуждать алгоритм решения уравнения (74) необходимо разобраться с рядом математических вопросов. Отметим, что само уравнение (74), получено из дифференциального уравнения (61) с помощью простых алгебраических преобразований. Поэтому, вопрос о существовании решения уравнения (74) не стоит, однако не очевидно, что любое решение уравнения (74) будет также решением уравнения (61). Может также оказаться, что одному решению уравнения (61) будет соответствовать несколько решений уравнения (74). Пусть у нас существуют две функции д1 и д2 , которые являются различными решениями уравнения (74), относящимися к одной и той же энергии е(-) . Для простоты рассмотрим случай, когда в исходном уравнении (61) эта энергия является невырожденной. Вычитая из уравнения (74) для функции д1 то же самое уравнение для функции д2 и приводя подобные члены, получим линейное однородное уравнение:
Ь(-)уа(т) = е(-)уа(т), (76)
где
у^(т) = д1(т) — д2(т)- (77)
Так как для уравнения (76) е(-) является невырожденным собственным значением, то уа(т) либо равна нулю, либо пропорциональна у(т):
у^т) = Лу(т). (78)
Следовательно,
д1(т) = д2(т) + Лу(т)
(79)
Так как нас интересует нормированное решение уравнения (61), то подставим (79) в
(73):
|у (г) | = |У (г) + 01 (г) | = |У (г) + 02 (г) + Лу(г)| = |(1 + А)у(г) | (80)
и, далее, подставим (80) в условие нормировки для функции у(г) и получим уравнение:
(1 + Л)2 = 1,Л1 = 0,Л2 = -2 . (81)
Таким образом, мы видим, что, в отличие от уравнения (61), уравнение (74) имеет два различных решения. Первое решение д1 действительно является малой добавкой к опорной функции У (г). Второе решение д2, соответствующее Л2 = —2 , приводит согласно (80) к функции -у(г). Эта функция в уравнении (61) не давала нового решения, так как решение линейного однородного уравнения находится с точностью до произвольной постоянной. Уравнение (74) не однородно и, очевидно, что функции д1 и д2 -различные решения. Поэтому при построении алгоритма численного расчета следует позаботиться о том, чтобы решение ур(г) было мало по сравнению с опорной функцией. В этом случае второе решение окажется исключенным. Процедура нормировки позволяет уменьшить ошибки, возникающие при организации итерационного процесса и не допустить его срыва.
Перейдем к построению итерационной процедуры решения уравнений (74). Эти уравнения для операторов Ь- и Ь+ одного типа. Метод решения для них один и тот же. Нужные решения обоих уравнений удовлетворяют нулевым граничным условиям при г = 0 и при г = то. Рассмотрим сначала простейший случай равномерной сетки. При этом вторая производная заменяется на разностную формулу:
#Ур(г) УР((п 1 + ОД + Ур((п! - 1)к) - 2ур(щк) , ^^ . .
-Р-+ 0(/°" (82)
Аналогично первая производная с той же точностью заменяется на разностную формулу:
¿г 2Н
Подстановка (82) и (83) в (74) позволяет получить итерационную процедуру:
, м АчУрЦп 1 + ВД + В^иДщ - 1)к) + /^^(щ/р Ур{п\Ь) =-- . (84)
Шп1
Здесь введены обозначения
Ли = 1 + , (85)
В, = . (86)
д1_, = 1 + Щ + 1 - + „2 Л,1(1!1/1) _ е,_, _ М
2п1 2
Р1{ЩЬ) -4с^щк^щк)- • (88)
Аналогичные формулы могут быть выведены и для уравнения
Ь(+)ур(т) = е(+)ур(т) + РР(+)(т) , (89)
где
РР(+) (т) = [(е(+) — е(0)) — (т)]У(+) (т). (90)
Подстановка (82) и (83) в (89) позволяет получить итерационную процедуру:
, м АП1ур{{гы + 1)к) + ВП1ур((П1 - 1)к) + /^РР^СЩк) Ур{П\Ь) =-—^- . (91)
Здесь введено обозначение
^ = ь + (1 + 1)(1 + 2(п,ЦР,(п,Ц) + н2 _ е<+) _ ЩЩ (92)
2п^ \ 2 /
При применении итерационных процедур равномерной сеткой с шагом к = 0, 000005 пользуемся только вблизи ядра п1 < 400. Дальше следует использовать сетку с возрастающим шагом:
к1(П2 + 1) = Мп2) ■ Ън . (93)
Для достижения нужной точности положим к1(1) = 0, 000005 и Ън = 1, 00344. При этом узлы сетки связаны простым соотношением:
т(п2 + 1) = т(п2) + к1(П2) . (94)
Здесь т(1) = 0,002 и вместо формулы (82) для п2 > 2 вторая производная дается выражением:
й2ур(т) _ 2 ■ (ур(т(п2 + 1)) + Ън ■ ур(т(п2 — 1)) — (1 + Ън) ■ уР(т(п2)))
(95)
йт2 к1(п2) ■ к1(п2 — 1) ■ (1 + Ън)
Аналогично первая производная с той же точностью заменяется на разностную формулу: 2 2
<1Ур{г) _ Ур(г(п2 + !))-&£• ур(г(п2 - 1)) + {Ь{ - 1) • ур(г(п2)) ^
йт к1(п2) ■ (1 + Ън)
При Ън = 1 формулы (95) и (96) переходят в формулы (82) и (83). Подстановка (95) и (96) в (74) позволяет получить итерационную процедуру:
АП2 ■ Ур(г(п2 + 1)) + ВП2 ■ Ур(г(п2 - 1)) + - I))2 • РР("}(г(/г2))
Р>1
У (г(п2)) = ~ У"-Н"2 - }-)) ± ± у/ \П:2)) ^^
'П2
Здесь введены обозначения:
Ап
Вп
= р - (Ып2 - I))2
ЬН
1 + Рг(г(п2)) • кгЫ
1 ~ РМщ)) • Ы(п2) • Ън % + 1)
(98)
(99)
¿(¿ + 1 " 2г(п2) ■ Рг(г(п2)))
РМп2)) • (Ън - 1) • Ыщ - 1) _ Ън
2 ■ г(п2)2
Для уравнения (89) рекуррентная формула для неравномерной сетки имеет вид: Лг2 • Ур(г(п2 + 1)) + ВП2 ■ ур(г(п2_- 1)) + (/н('/г2 - I))2 •
П
у. (г(п2)) = ' ' ~ У"Н"2 ~ '-)) ' 1 1 у У""2))
п2
где
\2
= т1 - Ып, -1))2 • к(гЫ) - е<+>) -
Ьн
(1 + 1)(1 + 2г(п2)-Р1(г(п2)))
РМп2)) • - 1) • - 1) _ Ьн
2 ■ г(п2)2
При п2 = 1 используются формулы для равномерной сетки, т.к. к1(1) = к. В правую часть уравнений (84), (91), (97) и (101) подставляем массивы из предыдущей итерации, а в левой части возникает новый массив ур(п1к). Чтобы массив ур(п1к) состоял из малых чисел, итерационную процедуру начинаем с нулевого массива.
В итерационные формулы (84) и (97) входит энергия е(-). Для ее вычисления умножим уравнение (74) слева на У-(г) и проинтегрируем по г. Дважды интегрируя по частям, под знаком интеграла можно перенести вторую производную с ур(г) на У-(г) и воспользоваться уравнением (72). В результате, получим следующее выражение для е(-):
оо оо
/ у(г)Уд(г)У(-)(г)^г + / Р1(г)(ур(г)У'(-)(г) — У(-)(г)уР(г))^г
е<"> = е<°> + ^-1- . (ЮЗ)
/ у(г)У(-\г)йг 0
Для энергии е+, входящей в (91) и (101), используя (89), аналогично получим выражение:
оо оо
/ у(г)Уд(г)У(+)(г)^г + / Р1(г)(ур(г)У'(+)(г) — У(+)(г)у'р(г))йг
= е<°> + ^-^- . (104)
/ у(г)У(+\г)в>г 0
Для численного расчета интегралов в (103) и (104) использовалась формула
г(П2+2)
f (т)в,т
(1 + ЪН)К
■П2
6
г(п2)
(1 + Ън)
2
/(гЫ)(2 - Ън) + 4 /(г(»2 + 1))
Ън
2Ъ _1
+ -/(г(П2 + 2))
Ън
(105)
которая обеспечивает нужную точность вычисления и в случае равномерной сетки переходит в формулу Симпсона. Формулы (103) и (104) можно использовать на любом этапе итерации т.к. даже при ур(г) = 0 они имеют смысл. При этом они переходят в формулы первого порядка теории возмущений. При ур(г) < 1 вторым интегралом в числителе можно пренебречь, т.к. функция Р\(т) почти везде также пренебрежимо мала.
Таблица 1
Результаты расчета полной энергии и энергий ионизации легких атомов и соответствующие им экспериментальные энергии, эВ
7 Ат. -К -^эксп -Е -ЁТэксп ~Е1 о1 оЗ о4
2 Не 79.0056 79.0052 24.5876 25.6368 0.0778 1.312 -0.62 0.158
3 и 203.4828 203.4882 5.3918 5.3922 -0.15327 0.64 0.94 -0.043
4 Ве 399.036 398.9891 9.32 9.3226 0.0879 1.011 1.3 -0.357
5 В 670.9941 670.99 8.2981 8.2982 0.1076 0.9531 1.3 -0.3342
6 С 1030.11 1030.12 11.2643 11.2643 0.2997 1.5092 2.0 -0.7431
7 N 1486.0671 1486.07 14.5341 14.5339 0.4003 1.1201 0.9 0.042
8 О 2043.866 2043.87 13.618 13.6181 0.1201 1.8340 -0.03 -0.64
9 Е 2715.878 2715.87 17.423 17.4235 0.2417 1.7757 -0.01 -0.34
10 N6 3511.598 3511.60 21.565 21.5650 0.3534 2.121 -0.026 -0.34
И Ма 4419.9461 4419.96 5.1391 5.1392 0.1913 1.2976 0.34 -0.017
12 М§ 5450.9453 5450.92 7.6463 7.6463 0.1854 1.4692 0.12 -0.09
13 А1 6613.3028 6613.30 5.9858 5.9816 0.2608 1.5072 0.3 -0.07
14 7888.4027 7888.41 8.1517 8.1517 0.3378 1.7742 0.24 -0.1
15 Р 9305.82 9305.81 10.4868 10.4858 0.4346 2.0574 0.24 -0.09
16 Б 10858.283 10858.30 10.3600 10.3598 0.2551 1.3219 0.4 0.149
17 С1 12555.362 12555.36 12.9680 12.9675 0.3363 1.5458 0.66 0.04
18 Аг 14397.801 14397.81 15.7600 15.7597 0.4294 1.4529 0.66 0.24
4. Вычислительный эксперимент. Проведем апробацию выведенных формул, на примере расчетов полных энергий релятивистских атомов от Не до Кг. Для нахождения усредненной по сфере электронной плотности п(т) и занятых уровней энергии е% самосогласованные уравнения Дирака (61) и (65) с потенциалом (48) решается итерационным методом опорной функции, изложенным выше. Для вычисления полной энергии будем пользоваться формулой (32), в которую подставим выражение для Уаг(т) из (48):
N Я
Е = ^ - 2п(Же - 1) / тп(т)(1 - ^(т) + 0.5дхсз(т) + дхср(т)^т . (106) ?=1 0
В табл. 1 и 2 приведены результаты расчета полной энергии Е атомов, а также соответствующие им экспериментальные значения энергий Еэксп. За экспериментальную полную энергию мы приняли взятую с обратным знаком сумму потенциалов ионизации атомов и атомных ионов, приведенных в табл. 19.1 и 19.2 справочника [22].
За экспериментальное значение энергия ионизации Е1эксп мы приняли взятое с обратным знаком значение потенциала ионизации из табл. 19.1 справочника [22]. Эти значения приведены с точностью от 10-2 до 10-4 эВ. За вычисленные значения энергии ионизации Ех мы приняли энергию последнего занятого уровня в атоме.
В табл. 1 и 2 приведены значения безразмерных параметров —в и а1, определяющих вид обменно-корреляционного потенциала. Они подбирались так, чтобы вычисленное значение полной энергии Е совпадало с Еэксп с точностью лучше 0.05 эВ, а значение Ех совпадало с Е\эксп с точностью лучше 10-3 эВ. Здесь же приведены значения параметров а3 и а4, которые подбирались так, чтобы значения энергий, полученных при решении уравнений КШД, были близки к соответствующим экспериментальным значениям энергий внутренних оболочек атомов.
В табл. 3 приведены результаты расчетов, позволяющие провести сравнение расчетных и экспериментальных энергий Еэксп внутренних оболочек атомов элементов второго и третьего периодов таблицы Менделеева. Для сравнения, наряду с нашими результатами Ерасч, приведены значения Ехф, полученные из решений уравнений Хартри-Фока, взятые из монографии [23].
Таблица 2
Результаты расчета полной энергии и энергии ионизации переходных элементов, и соответствующие им экспериментальные энергии, эВ
г Ат. -К -^эксп -Е -ЁТэксп ~Е1 о1 оЗ о4
19 К 16379.71 16379.72 4.3407 4.3407 0.0002 1.38914 0.34 -0.0133
20 Са 18508.08 18508.07 6.1132 6.1129 0.1159 1.7157 0.12 -0.09
21 Бс 20786.15 20786.15 6.5615 6.5614 0.33204 2.34974 0.2 -0.31
22 Т4 23221.41 23221.404 6.8200 6.8199 0.4734 2.9272 -0.4 -0.18
23 V 25820.80 25820.801 6.7400 6.7401 0.40056 1.91496 -0.14 0.26
24 Сг 28586.01 28586.01 6.7660 6.7661 0.33846 1.97585 -0.19 0.17
25 Мп 31536.74 31536.74 7.4340 7.4337 0.4081 1.3061 0.04 0.57
26 Ее 34651.58 34651.59 7.9024 7.9014 0.3444 2.0767 -0.14 0.112
27 Со 37896.04 37896.04 7.8600 7.8597 0.2886 2.0593 -0.14 0.112
28 N1 41381.31 41381.30 7.6370 7.63698 0.25036 1.8948 -0.14 0.19
29 Си 44956.26 44956.26 7.7264 7.7260 0.2016 2.0194 0.06 0.04
30 Ъ\\ 48785.65 48785.64 9.3943 9.3934 0.1387 1.9074 -0.14 0.16
31 ва 52815.44 52815.45 5.9993 5.9993 0.16562 2.10077 -0.046 0.03
32 ве 57009.25 57009.26 7.8995 7.8995 0.27416 2.33704 -0.046 0.03
33 Ав 61402.48 61402.48 9.7890 9.7889 0.3971 2.9085 0.18 -0.18
34 Бе 65982.01 65982.00 9.7520 9.7519 0.2140 2.4547 0.18 -0.16
35 Вг 70758.92 70758.90 11.8140 11.81395 0.29701 2.8919 -0.004 -0.22
36 Кг 75724.15 75724.15 13.9997 13.9999 0.4243 2.0277 0.66 0.24
Таблица 3
Сравнение расчетных и экспериментальных энергий внутренних оболочек атомов элементов второго и третьего периодов таблицы Менделеева, эВ
4 Ве 5 В 6 С
-Е 1 ^|>1 (г г —Ехф -Е -^эксп -Е -^ра.сч —Ехф -Е -^эксп -Е 1 ^|>1 (г г —Ехф -Е -^эксп
ь 126.11 128.78 123.60 204.74 209.40 201 296.76 308.20 296
2в 9.3226 8.4165 9.3227 12.97 13.46 12.93 18.82 19.20 16.59
Е 398.989 396.552 399.036 670.997 667.470 670.994 1030.12 1025.56 1030.12
7 N 8 О 9 Е
-Е 1 ^|>1 (г г —Ехф -Е -^эксп -Е -^ра.сч —Ехф -Е -^эксп -Е 1 ^|>1 (г г —Ехф -Е -^эксп
Ь 401.835 425.287 403 564.65 562.43 538 708.93 717.92 694
2в 24.799 25.723 20.330 29.15 33.86 28.48 36.26 42.79 37.86
2Р1/2 14.534 15.445 14.534 13.618 17.195 13.618 36.26 42.79 37.86
Е 1486.07 1480.33 1486.07 2043.87 2035.67 2043.87 2715.87 2705.07 2715.88
10 N6 И Ма 12
-Е 1 ^|>1 (г г —Ехф -Е -^эксп -Е -^ра.сч —Ехф -Е -^эксп -Е 1 ^|>1 (г г —Ехф -Е -^эксп
Ь 868.45 891.77 870.27 1079.61 1101.49 1079 1327.63 1334.23 1311.30
2в 44.64 52.53 48.47 59.93 76.11 70.90 83.57 102.52 96.50
2Р1/2 21.681 23.141 - 35.35 41.31 38.46 54.00 62.10 57.60
2^3/2 21.565 23.141 21.565 35.11 41.31 38 53.64 62.10 57.60
Зв - - - 5.1392 4.9552 5.1391 7.646 6.884 7.646
Е 3511.601 3497.948 3511.598 4419.96 4404.41 4419.95 5450.92 5431.79 5450.95
13 А1 14 15 Р
-Е 1 ^|>1 (г г —Ехф -Е -^эксп -Е -^ра.сч —Ехф -Е -^эксп -Е 1 ^|>1 (г г —Ехф -Е -^эксп
Ь 1569.10 1591.89 1567 1842.43 1872.47 1844 2132.01 2176.09 2148
2в 113.50 133.63 126 145.31 167.53 154 181.30 204.38 191
2^1/2 80.29 87.57 81 106.11 115.81 104 136.02 146.97 135
2рз/2 79.75 87.57 80 105.36 115.81 104 135.00 146.97 134
Зв 10.96 10.70 10.62 14.92 14.69 13.46 19.12 18.95 16.15
ЗР1/2 5.982 5.714 5.986 8.152 8.082 8.152 10.49 10.66 10.49
16 Б 17 С1 18 Аг
-Е 1 ^|>1 (г г —Ехф -Е -^эксп -Е -^ра.сч —Ехф -Е -^эксп -Е 1 ^|>1 (г г —Ехф -Е -^эксп
ь 2445.89 2503.56 2476 2773.71 2853.93 2829 3128.88 3227.54 3206
2в 224.43 245.02 232 268.82 288.66 277 315.11 335.30 326.37
2^1/2 179.27 181.84 170 218.10 219.66 208 256.80 260.45 250.60
2рз/2 177.81 181.84 168 216.20 219.66 206 254.41 260.45 248.60
Зв 20.12 23.94 20.20 24.57 29.20 24.59 29.48 34.76 29.24
3^1/2 10.45 11.90 10.36 13.100 13.783 12.968 15.95 16.08 15.76
3^3/2 10.36 11.90 10.36 12.967 13.783 12.968 15.76 16.08 15.76
В табл. 4-6 приведены аналогичные данные для элементов четвертого периода таблицы Менделеева. Результаты расчетов Ехф также взяты из монографии [23]. Как видно из этих таблиц, полученные нами значения в большинстве случаев ближе к экспериментальным чем аналогичные результаты расчетов, полученные по методу Хартри-Фока.
Таблица 4
Сравнение расчетных и экспериментальных рентгеновских спектров элементов четвертого периода таблицы Менделеева, эВ
19 К ь 2в 2^1/2 2^3/2 Зв 3^1/2 Зрз/2 4в
— Е 3585.99 374.67 324.46 320.98 32.17 17.43 17.14 4.34
—Ехф 3633.54 394.29 313.45 313.45 47.58 25.97 25.97 4.01
-Е -^эксп 3614 384 303.3 300.7 37 24.82 24.49 4.34
20 Са Ь 2в 2Р1/2 2рз/2 Зв 3^1/2 Зрз/2 4в
— Е 4051.66 423.69 367.11 362.69 45.75 28.84 28.40 6.11
—Ехф 4064.29 457.78 370.86 370.86 61.1 36.48 36.48 5.32
-Е -^эксп 4048 447 360 356 48 34.7 34.3 6.11
21 5с Ь 2в 2^1/2 2рз/2 Зв 3^1/2 ЗРЗ/2 4в
— Е 4489.12 482.55 413.04 407.83 58.01 38.45 37.88 7.75
—Ехф 4514.37 519.22 426.35 426.35 69.86 42.85 42.85 5.72
-Е -^эксп 4494 503 408 404 56 33 33 -
22 Т% Ь 2в 2^1/2 2^3/2 Зв 3^1/2 Зрз/2 4в
— Е 4959.01 539.40 455.41 449.27 63.98 41.62 40.91 7.94
—Ехф 4987.03 582.94 484.12 484.12 78.19 48.85 48.85 6.0
-Е -^эксп 4970 567 465 459 64 39 38 -
23 V Ь 2в 2^1/2 2рз/2 Зв 3^1/2 Зрз/2 4в
— Е -^Орэ.сч 5448.31 606.39 526.01 518.22 70.99 47.10 46.19 8.19
—Ехф 5483.09 649.67 544.82 544.82 86.62 54.95 54.95 6.27
-Е -^эксп 5470 633 525 518 72 44 43 -
24 Сг Ь 2в 2^1/2 2^3/2 Зв 3^1/2 Зрз/2 4в
— Е -^Орэ.сч 5993.54 676.66 595.43 585.77 78.39 52.59 51.44 8.37
—Ехф 5997.12 713.18 602.46 602.46 89.39 55.80 55.80 6.04
-Е -^эксп 5995 702 589 580 80 49 48 -
В табл. 7 сравнивались результаты наших нерелятивистских расчетов с результатами расчетов, полученных по методу Хартри-Фока для атомов от Ga до Кг. При этом для каждого из них указано новое значение параметра а\, подобранное так, чтобы энергия последнего занятого уровня в соответствующем атоме совпадала с £-1ЭКСП с точностью лучше 10-3 эВ. Параметры в, а3 и а4 оставлены такими же, как и при релятивистских расчетах. Полученная при нерелятивистских расчетах полная энергия, естественно не совпадает с экспериментальной. Эти значения, наряду с соответствующими значениями полных энергий Хартри-Фока, приведены в последней строке табл. 7. Видно, что рассчитанные нами значения полных энергий, ближе к экспериментальным, чем энергии Хартри-Фока. В большинстве случаев рассчитанные нами значения энергий внутренних оболочек также оказываются ближе к экспериментальным, чем соответствующие значения, рассчитанные из уравнения Хартри-Фока.
Таблица 5
Сравнение расчетных и экспериментальных рентгеновских спектров элементов четвертого периода таблицы Менделеева, эВ
25 Mn 26 Fe 27 Co
-e —ex Ф -e ^эксп -e -^ра.сч —ехф -e ^эксп -e -^ра.сч —ехф -e -^эксп
1 6478.34 6545.16 6544 7117.11 7112.24 7117 7712.72 7702.73 7715
2s 758.48 792.12 755 831.17 869.02 851 913.27 948.83 931
2Р1/2 670.82 675.20 663 743.02 745.97 726 823.12 819.61 800
2^3/2 659.54 675.20 652 729.11 745.97 713 806.58 819.61 785
3s 86.77 103.85 89 94.72 113.46 98 102.75 123.12 107
3pi/2 58.12 67.47 55 64.36 74.62 61 70.09 81.81 68
Зрз/2 56.79 67.47 53 62.68 74.62 59 68.09 81.81 66
4s 8.37 6.74 - 9.04 7.03 - 9.10 7.28 -
3ci3/2 7.434 17.38 7.434 8.03 17.60 - 8.03 18.38 -
3ci5/2 - - - 7.90 17.60 7.90 7.86 18.38 7.86
28 Ni 29 Cu 30 Zn
-e —ex Ф -e ^эксп -e -^ра.сч —ехф -e ^эксп -e -^ра.сч —ехф -e -^эксп
ls 8337.39 8316.35 8338 8951.42 8946.83 8986 9656.94 9613.78 9663
2s 1005.20 1031.80 1015 1098.94 1110.71 1103 1197.83 1207.15 1198
2^1/2 915.42 896.40 877 1003.00 969.19 958 1104.81 1059.20 1052
2^3/2 895.78 896.40 860 980.39 969.19 938 1078.02 1059.20 1029
3s 111.50 133.00 117 119.10 136.36 127 132.70 153.41 141
3^1/2 76.41 89.19 75 81.34 90.46 82 92.49 1.0447 98.70
3i?3/2 74.04 89.19 73 78.57 90.46 80 89.24 1.0447 96.10
4s 8.99 7.52 - 9.17 6.48 - 9.393 7.96 9.394
3^3/2 7.857 19.23 - 8.006 13.35 - 10.74 21.30 -
3ci5/2 7.637 19.24 7.637 7.726 13.35 7.726 10.37 21.30 -
Заключение. Рассмотренное в настоящей статье применение метода опорной функции, позволило на порядок улучшить точность решения уравнения КШД. Предложенная аппроксимации обменно-корреляционного потенциала вывела локальное приближение в один ряд с наиболее точными из существующих приближений метода функционала электронной плотности.
Отметим, что потенциалы ионизации ионов приведены в [22] с точностью от тысячных долей эВ до 0.5 эВ. Например, для криптона 19 из 36 значений потенциалов ионизации ионов приведены с точностью до 0.5 эВ. Т.о. в этом случае погрешность измерения полной энергии составляет порядка 9 эв. До настоящего времени это не играло большого значения, т.к. погрешность существовавших ранее методов составляла десятки электрон вольт. Для метода опорной функции, приведенные в табл. 1 и 2 значения экспериментальных полных энергий являются условными числами, по отношению к которым точность рассчитанных нами энергий составляет менее 0.05 эВ. При появлении более точных экспериментальных данных, константы а\, ß, а3 и а4 могут быть пересчитаны без изменения описанного метода.
Авторы благодарят А.С. Старовойтова за обсуждения и ценные замечания.
Таблица 6
Сравнение расчетных и экспериментальных рентгеновских спектров элементов четвертого периода таблицы Менделеева, эВ
31 Ga 32 Ge 33 As
-e —Ex ф -e ^эксп -e —Ex ф -e ^эксп -e —Ехф -e ^эксп
Is 10354.60 10290.80 10371 11069.31 11027.14 11107 11790.81 11771.38 11871
2s 1303.60 1310.72 1302 1401.74 1419.07 1413 1508.63 1532.27 1531
2P1/2 1203.87 1156.35 1146 1287.92 1258.15 1251 1377.24 1364.76 1362
2^3/2 1173.16 1156.35 1119 1253.28 1258.15 1220 1338.64 1364.76 1327
3s 149.70 174.01 162 170.57 195.67 184 193.70 218.50 208
3^1/2 106.61 121.97 111 124.14 140.45 130 143.50 160.02 151
3i?3/2 102.78 121.97 107 119.65 140.45 125 138.30 160.02 145
4s 12.054 11.551 11.870 15.61 15.06 14.28 19.53 18.66 18.96
3ci3/2 18.39 32.47 21 29.92 44.49 33 44.14 57.49 46
3ci5/2 17.91 32.47 20 29.32 44.49 32 43.39 57.49 45
4^1/2 5.999 - 5.999 7.899 7.818 7.899 10.12 10.05 -
4^3/2 - - - - - - 9.789 10.05 9.789
34 Se 35 Br 36 Ivr
-e —Ex ф -e ^эксп -e —Ex ф -e ^эксп -e —Ехф -e ^эксп
Is 12592.21 12540.91 12662 13432.59 13335.21 13481 14091.93 14154.55 14327
2s 1653.84 1650.89 1656 1761.45 1774.18 1787 1896.13 1902.16 1923
2pi/2 1532.37 1476.73 1479 1629.57 1593.34 1602 1730.04 1714.59 1731
2^3/2 1487.34 1476.73 1439 1578.53 1593.34 1556 1676.06 1714.59 1678
3s 212.40 243.06 234 239.03 268.63 262 266.46 295.22 293
3pi/2 159.86 181.28 173 183.28 203.49 197 205.21 226.71 222
Зрз/2 1.5380 181.28 166 176.20 203.49 189 197.13 226.71 214.60
4s 20.19 22.79 22.19 24.04 27.01 23.80 27.92 31.37 27.51
3ci3/2 51.16 72.10 61 66.40 87.63 77 88.72 104.09 95
3^5/2 50.20 72.10 60 65.29 87.63 76 87.28 104.09 93.80
4^1/2 10.12 10.97 - 12.34 12.43 - 14.66 14.26 -
4p3/2 9.752 10.97 9.752 11.814 12.43 11.814 13.999 14.26 13.999
Литература
1. Слэтер Дж. Методы самосогласованного поля для молекул и твердых тел / М.: Мир, 1978. - 658 с.
2. Шкловский А.Г. Аппроксимация обменно-корреляционного потенциала в методе функционала электронной плотности // Научные ведомости БелГУ. Серия Математика. Физика. - 2007. - №6(37);13. - С.150-155.
3. Rajagopal A. K. Advances in Chemical Physics / ed. I. Prigogine, S.A. Rice / 41, New York: Wiley, 1980 - P.59.
4. Теория неоднородного электронного газа / под ред. С. Лундквиста, Н. Марча /-М.: Мир, 1987. - 400 с.
5. Langreth D.C., Perdew J.P. // Phys. Rev. - 1977. - B 15. - P.2884.
6. Hohenberg P., Kohn W. Inhomogeneous Electron Gas // Phys. Rev. B. -1964. - 136. - P.864-871.
7. Koch W., Holthausen M.C. A Chemists Guide to Density Functional Theory / Wiley-VCH Verlag GmbH, 2001. - 293 с.
8. Curtiss L.A., Raghavachari K., Trucks G.W., Pople J.A. Gaussian-2 Theory for Molecular Energies of First-and Second-Row Compounds //J. Chem. Phys. - 1991 - 94. - P.7221.
9. Kohn W. Sham L.J. Self-Consistent Equations Including Exchange and Correlation Effects // Phys. Rev. - 1965. - A140. - P.1133-1138.
10. Gunnarsson O., Lundqvist B.I. Exchange and Correlation in Atoms, Molecules, and Solids by the Spin-Density-Functional Formalism // Phys. Rev. B. - 1976. - 13. - P.4274.
11. von Barth U. Local-Density Theory of Multiplet Structure // Phys. Rev. - 1979. - A20. -P.1693.
12. Vosko S.J., Wilk L., Nusair M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis // Can. J. Phys. - 1980. - 58. -P.1200.
13. von Barth, U., Hedin, L. // J. Phys. - 1972. - 5. - P.1629.
14. von Barth U. Basic Density-Functional Theory-Overview // Physica Scripta. - 2004. - 109. -P.9-39.
15. Старовойтов А.С., Шкловский А.Г. Обменно корреляционная энергия в методе локального функционала электронной плотности // Научные ведомости БелГУ. Серия Математика. Физика. - 2006. - №6(26);12. - C.114-121.
16. Старовойтов А.С., Шкловский А.Г. Модифицированный локальный потенциал для вычисления энергии ионизации атомов // Научные ведомости БелГУ. Серия Математика. Физика, 2010. - №11(82);19. - С.126-134.
17. Бьеркен Дж. Д., Дрелл С.Д. Релятивистская квантовая теория, т. 1 / М.: Наука, 1978. -296 с.
18. Сидельников Г.Л., Шкловский А.Г. Алгоритм численного решения уравнения Дирака в центральносимметричном поле // Электромагнитные волны и электронные системы. -2004. - 9;12. - С.4-9.
19. Вержбицкий В.М. Численные методы (линейная алгебра и нелинейные уравнения) / Учеб. пособие для вузов. 2-е изд. / М.: ООО Издательский дом «ОНИКС 21 век», 2005. -432 с.
20. Шкловский А.Г., Шкловская М.А. Решение уравнения Кона-Шема для молекулы методом опорной функции // Научные ведомости БелГУ. Серия Физика. Математика. -2008. - №(49);14. - C.123-136.
21. Шкловский А.Г. Решение уравнения Кона-Шема для молекулы водорода методом опорной функции // Научные ведомости БелГУ. Серия Математика. Физика. - 2011. - №11(106);23. -С.52-64.
22. Бабичев А. П., Бабушкина И.А., Братковский А.М., и др. Физические величины / Справочник / М.: Энергоатомиздат, 1991. - 1232 с.
23. Фудзинага С. Метод молекулярных орбиталей / М.: Мир, 1983. - 462 с.
APPROXIMATION OF LOCAL DENSITY FUNCTIONAL
WITH EXCHANGE-CORRELATION ENERGY OF RELATIVISTIC ATOM
A.V. Beregovoj, A.A. Pleskanev, A.G. Shklovskij
Belgorod State University, Pobedy St., 85, Belgorod, 308015, Russia, e-mail: Shklovsky@bsu.edu.ru
Abstract. The method of total energy calculation of spherical symmetric relativistic atom by local electronic density functional with approximated exchange-correlation energy is described. The comparison of computed values of full energy of atoms from He to Kr with the experimental ones is done.
Key words: density functional theory, exchange-correlation potential, Kohn-Sham-Dirac equation, relativistic atom, support function method.
Таблица 7
Сравнение нерелятивистских расчетных и экспериментальных рентгеновских спектров элементов четвертого периода таблицы Менделеева, эВ
31 С-й, а1 = 1.9877 32 ве, в! = 2.1309 33 Ав,01 = 2.2338
-Е -^ра.сч —Ехф -Е -^эксп -Е 1 ^|>1 (г г —Ехф -Е -^эксп -Е -^ра.сч —Ехф -Е -^эксп
ь 10278.29 10290.80 10371 10987.12 11027.14 11107 11757.63 11771.38 11871
2в 1302.65 1310.72 1302 1403.06 1419.07 1413 1546.60 1532.27 1531
2р 1212.80 1156.35 1119 1303.24 1.258.15 1220 1446.99 1364.76 1327
Зв 152.51 174.01 162 173.57 195.67 184 196.75 218.50 208
3р 113.74 121.97 107 131.20 140.45 125 151.32 160.02 145
4в 11.97 11.55 11.87 15.43 15.06 14.28 18.59 18.66 18.96
М 26.97 32.47 20 38.86 44.49 32 52.88 57.49 45
4р 5.997 - 5.999 7.899 7.818 7.899 9.79 10.05 9.79
Е 52410.45 52334.88 52815.44 56589.64 56473.43 57009.25 61321.91 60796.72 61402.48
34 Бе, в! = 2.163 35 Вг, в! = 2.3921 36 Кг,си = 1.7391
-Е -^ра.сч —Ехф -Е -^эксп -Е 1 ^|>1 (г г —Ехф -Е -^эксп -Е -^ра.сч —Ехф -Е -^эксп
Ь 12480.45 12540.91 12662 13348.22 13335.21 13481 13877.70 14154.55 14327
2в 1665.20 1650.89 1656 1783.82 1774.18 1787 1878.71 1902.16 1923
2р 1563.02 1476.73 1439 1683.55 1593.34 1556 1722.12 1714.59 1678
Зв 215.35 243.06 234 243.08 268.63 262 269.71 295.22 293
3р 167.29 181.28 166 192.12 203.49 189 212.50 226.71 214.60
4в 19.36 22.79 22.19 22.91 27.01 23.80 26.70 31.37 27.51
■ы 63.17 72.10 60 78.37 87.63 76 106.22 104.09 93.80
Ар 9.74 10.97 9.75 11.813 12.436 11.814 14.00 14.26 14.00
Е 65561.02 65303.69 65982.01 70424.04 69999.75 70758.92 74958.36 74887.21 75724.15
