CC BY
45
Вестник Челябинского государственного университета. 2013. № 25 (316). Физика. Вып. 18. С. 75-79.
МЕТОДИЧЕСКИЕ ЗАМЕТКИ
А. В. Береговой, А. Г. Шкловский
РЕШЕНИЕ УРАВНЕНИЯ КОНА — ШЕМА ДЛЯ ЛИНЕЙНЫХ СИММЕТРИЧНЫХ МОЛЕКУЛ
Предложен вариант модифицированного локального потенциала, пригодный для вычисления электронного спектра, полной энергии и энергии связи многоэлектронных симметричных молекул. Качество предложенных алгоритмов проверено при вычислении полной энергии и энергии связи молекулы углерода.
Ключевые слова: уравнение Кона — Шема, метод опорной функции, модифицированный локальный обменно-корреляционный потенциал.
1. Введение
Будем пользоваться атомной системой единиц И = є2 = те = 1. Расположим ядра атомов вдоль оси 2. Начало координат выберем в центре молекулы Расстояния от ядер до произвольной точки г не зависят от угла ф . Будем предполагать, что молекулярный потенциал ¥(г) тоже не зависит от угла ф:
¥(г) = Г(р,г) . (1)
Это позволяет свести решение трехмерного уравнения Кона — Шема
(2)
- - + V (г) |Ч я (г) = Еп Ч п (г)
к численному решению двумерного уравнения . В цилиндрической системе координат уравнение
(2) с учетом условия (1) можно переписать в следующем виде [1]:
д2 д2 ш2 - 0,25
др2 дг2
р
= 2Е и
п,т п,т
+ 2Г (р, г) (Р.2 )•
ип,ш (р, г) =
(3)
При этом решение уравнения (3) для функции ¥ (г) имеет вид
л/р
А (фХ
(4)
где Ат(ф) — система ортогональных функций:
Ат(ф>) = шз(т • ф), Ая(ф) = sin(m • ф); (5)
т — целое число .
Из выражения (4) видно, что при р = 0, ип т(0,£) = 0 при любом г. Очевидно, что на машине мы не можем реализовать бесконечные отрезки. Поэтому вводятся числа рт и гт, за пределами которых можно считать решения уравнения (3) равными нулю . Следовательно, будем решать уравнение (3) с нулевыми граничными условиями:
ии,т(Р,гт) = 0 и„,т(Рт,г) = 0 и (р,-г ) = 0, и (0,г) = 0 . (6)
7 т' 7 п.тк 7 ' 4 '
Если сразу решать уравнение (3) численно, то возникают две трудности. Во-первых, уравнение
(3) является уравнением на собственные функции и значения . Качественных двумерных численных алгоритмов для решения таких задач до последнего времени не было . Во-вторых, замена второй производной на ее дискретный аналог может быть произведена с не очень высокой точностью . Это связано с резким увеличением количества узлов сетки при уменьшении шага, или с уменьшением порядка аппроксимации при использовании неравномерной сетки. Обе эти трудности и приводили ранее к тому, что вместо прямых численных алгоритмов использовались проекционные методы
Кроме этого, для молекул с одинаковыми атомами есть дополнительная сложность, связанная с тем, что гамильтониан молекулы коммутирует с оператором отражения в плоскости ХУ. Поэтому функция и (рг) должна быть собственной функцией этого оператора. Следовательно, либо функция и (рг) симметрична: и (рг) = и (р,-г), либо
п,тК1 7 ' А з,п,тУ1 7 у з,п,т41 7 '7
и (рг) антисимметрична: и (р,г) = и (р,-г) .
п,т А ап,т а,п,т41 7 '
Эти трудности позволяет преодолеть метод опорной функции [2]. Поэтому далее именно этот метод использован для численного решения уравнения (3)
2. Вычисление молекулярного потенциала
Для придания дальнейшим вычислениям конкретности рассмотрим молекулу углерода В качестве базисных функций УДрХ), входящих в симметричные или антисимметричные линейные комбинации [2], будем использовать атомные функции
(у1р2+22).РГ Ґ \ z
/ 2 . 2
1л/Р +2 J
л/р- (7)
Здесь Р^х) — присоединенный полином Лежандра, а Яп (г) решение радиального уравнения Кона — Шема для атома углерода:
1 д_ 2r2 дг
дКг (г)
дг
Rn,l (Г) = 6n,l ' Rn,l (r)•
/ • (/ +1) 2r2
(8)
В работе [3] был введен модифицированный локальный потенциал (МЛП), с помощью которого удалось очень точно вычислить как полную энергию, так и энергию ионизации релятивистских атомов В рамках теории функционала электронной плотности в работах [4; 5] было получено выражение для полной нерелятивистской энергии атома:
Ыеа
Еа = Е ^ “I ¥а (Г)Па (ГМГ + ЕЬаг [Па ] +
( =1
+ЕХС [Па ] + _[ ¥е (г)Па (Г)*. (9)
Здесь энергия Хартри Екаг[па] дается обычным выражением; i нумерует решения уравнений Кона — Шема (8); Уа(г) — сферически симметричный потенциал атома. Для атома углерода, как и в работе [2], считаем, что уравнение Кона — Шема (8) описывает систему невзаимодействующих квазичастиц с энергией еп находящихся во внешнем потенциале К (г) .
В приближении МЛП [3] функционал обменнокорреляционной энергии Е [п ] имеет вид
Exc К ] = ■
1
2 N.
3Р„ (Nea -1)
х f dr1 • na(r1 )3 + f dr1dr2 Па(fl)Wa (Г2)
I r1 - r2 |
+ k f dr1dr2 F (r1)na (,1) F W' J if --- f l
(10)
где к = (3 + ^Жа-1);
^(х) — аппроксимирующая функция, введенная в работе [2]:
^(х) = хРа(х)[а3 + а4Ра(х) + (1 - а3) Ра2(х)] х
х ехр(-а2х), (11)
где
(ха1 -1)
Pa(x) = -
(12)
(x + 0,5)
Отметим, что в этом приближении электрон-
ную плотность атома
Па (Г) = 2Z Rn,l (r )2
n,l
(13)
мы получаем из обычной формулы Кона — Шема усреднением по углам 0 и ф . Процедура усреднения очень важна для приближения МЛП, т к универсальная функция ^ подбирается для сферически симметричной плотности атома В работе
[5] приведены веские соображения в пользу того, что хорошее качество локального приближения в теории функционала электронной плотности для атома связано как раз с использованием для нее сферически симметричного приближения
В формуле (11) коэффициент а2 определяет скорость убывания функции ^(х) вдали от ядра атома или иона, а коэффициент а вычисляется исходя из условия
| F (r )n(r )r2 dr = 0,
(14)
которое обеспечивает выполнение закона сохранения перераспределенного заряда обменнокорреляционной дырки, усредненной по сфере Здесь — радиус сферы, вне которой электронной плотностью атома или иона можно пренебречь . Обычно при численных расчетах в атоме мы выбираем Яс = 12 радиусов Бора . В молекулярных расчетах для опорных функций достаточно выбрать Яс = 8 .
Вернемся к вопросу о нахождении потенциала У(р^) для молекулы . При переходе от атома к молекуле приходится несколько изменить выражение (9) для полной энергии:
E=£ E
)=1
+E,,
-j(V (г) - Ve (г) )и(г )d
г +
r [n] + Exc [n] + :
(15)
где Net — число электронов;
E — энергии молекулярного спектра, полученные в результате решения уравнения Кона — Шема (2);
V(r) — молекулярный потенциал, входящий в уравнение (2);
n(r) — электронная плотность в молекуле;
|а| — расстояние между ядрами .
Энергия Хартри в (15) дается обычным выражением, которое отличается от E [n] в формуле (9) только заменой электронных плотностей в атоме на электронные плотности в молекуле . Постоянная кулоновская энергия отталкивания ядер друг от друга вводится для того, чтобы молекула была электронейтральна
При приближении к реальному межъядерному расстоянию в молекуле электронная плотность n(r) значительно отличается от суммы атомных плотностей . Чтобы это учесть, запишем ее в виде
n(r) = 2^\W и (r)|2 =
= n„
f a f a
r — 1 + Па 2 r + —
V 2 V 2
+4d (r). (16)
0
a
n
Здесь ¥п(г) является решением уравнения (2), а п/г) показывает отличие реальной плотности молекулы от суммы сферических атомных электронных плотностей п п2 и находится самосогласованным образом. Приближение для обменно-корреляционной энергии Е [п] тоже изменяется по сравнению с (17) . Дело в том, что плотность п(г) в молекуле будет обладать уже не сферической, а цилиндрической симметрией. Эта плотность, как видно из (16), имеет два центра сгущения. Поэтому введем 2 псевдоатома в молекулу При этом для псевдоатомов с ядрами в точках г = ±|а|/2 имеем электронные плотности ПЬ(р,г) и па(р,г), определяемые выражениями
( I 1\2^
2 I N I
_ а • р + 2 л—
I 2)
Л b (P, z) = n(P, z) • exP
Ла (P, Z) = n(P, Z) - ЄХР
-а<
P2 +
z-
2\ і
(17)
г a і 2 2 г a і
2 + Р + z — — , Г =4 р + z +—
V 1 2 J 1 2 J
(21)
Функция Щ(г) берется в виде Щг) = г1 ■ Ра(г1)[а3 + а4Ра(г1) + (1 - а3)Ра2(г1)] х X ехр(-а2г1 - а/2) + Г2 ■ Ра(г2)[а3 + а4Ра(г2) +
+ (1 - а3)Ра2(г2)] ■ ехр^^ - а/2), (22)
т. к. для молекулы с одинаковыми атомами все значения параметров совпадают.
Напомним, что в приближении МЛП мы сначала производим усреднение электронной плотности по сфере, и лишь после этого вводим аппроксимирующую функцию Щ (х) и потенциал .
Введем усредненные по сферам с центрами в точках р = 0, г = ±|а|/2 плотности псевдоатомов п(г) и Пъ(т):
я
Па 00 = 2я\Па (Р. 2)вІП ©!• (23)
Очевидно, что функция пЬ(Р,г) описывает электронную плотность псевдоатома, центрированного в точке р = 0, г = — я/2 и получается из молекулярной электронной плотности п(р, г) с учетом вероятности того, что электрон в точке (р,г) принадлежит к псевдоатому Ь. Максимум этой вероятности равен 1 и достигается в точке (0,—а/2) . Как и раньше введем граничные значения рт и г за пределами которых можно считать, что электронная плотность псевдоатома обратилась в ноль
Эти плотности нормируются на N — число электронов в псевдоатомах а и Ь:
Здесь
2П Jpdp J Ца (p, z)dz = Nei
(18)
Здесь рт — радиус, гт — высота цилиндра, содержащего псевдоатомы в молекуле Для углерода Иеа = 6 . Параметр а5 однозначно определяется из условий (18) . Для молекулы теперь можно ввести обменно-корреляционную энергию в приближении МЛП:
Ехс [п] = | п(г1) • ^хс [п(Г1 )] • , (19)
где ц [п] дается выражением
z)
1
n(p, z)3 н----------x
<f dr*?!)- н k. F (p, z) fF (r'-)n(r'-)dv1
J It* -- f I J I I* -- I* I
. (20)
Здесь k = (3 + N 1/3)(N -1)/N , а P(pz) имеет вид
P(P,z)=
в • (Nea -1)
N„„
[exP(-a5 • Г?) + exp(-a5 • r22)] ,
Z -~^- = rj • cos 0J, P = r • sin 0J,
z +2= r2'cos ®2 ’ P = r2' sin ®2 •
Законы сохранения усредненных по сфере перераспределенных зарядов обменно-корреляционных дырок
J Fa (r )na (r )r2 dr = 0.
(24)
Из этих условий и определяются новые коэффициенты а1 в формуле (12) .
Физический смысл ц [п] очевиден — это обменно-корреляционная энергия в приближении МЛП, приходящаяся на один электрон в молекуле . С учетом выражений (19)-(24) функционал полной энергии молекулы (15) известен с точностью до аппроксимирующих констант а и в .
Согласно теореме Хоэнберга и Кона [7], функционал полной энергии (15) для точной электронной плотности (16) имеет минимум . Так как минимум функционала (15) определяется из равенства нулю его первой вариации, то при очень малом изменении электронной плотности п(г), например, за счет малого изменения потенциала К(г), поправки к полной энергии Е будут второго порядка малости Поэтому вблизи от правильной электронной плотности модифицированный функционал можно считать точным вплоть до второго порядка малости Следовательно, можно брать его вариацию по электронной плотности, считая параметры
о
О
о
аппроксимации постоянными. Вычисляя вариацию функционала (15) по электронной плотности п(г) и приравнивая ее к нулю, получим выражение для неизвестного молекулярного потенциала Е(г): V (г) = ¥е (г) + ¥хса (X!) + ¥хсЬ (л2) +
+ Паг (г) • (Мет - !)/Нет + Vа (г). (25)
Здесь Е(г) — потенциал притяжения электронов к покоящимся ядрам в молекуле, выражение для потенциала ¥Ыг(г) с учетом (16) имеет обычный вид
^ (г) = .[^Ц <1г2Меа х г - г,
Ґ ( \\ Ґ ( \\
а і а 11
1 ~ Чкаг Г — 9 1 а 1 Г + — 9 )
V V 2 )) V V 2 ))
—+ -
а а
Г — Г + —
2 2
(26)
Здесь сферический заряд дка(г) дается формулой
Яьаг ) = 4П -) | Х1па (Х) (\ - -1 ^ (27)
а дополнительный потенциал V (г) определяется формулой
(Г2).
V,
:йг2.
(28)
Потенциал Слэтера ^;(р^) = Р(р,г)^и(р,г)1/3, а
Кса(Хі) и КсЬ(Х2) имеют вид
^(X!) = к ■ Еа(X!) І Е“{Гі)йг.
а
2
4. Вычислительный эксперимент
С учетом (19)—(21) и (25)-(29) полная энергия молекулы углерода Е (|а|) может быть вычислена по формуле (15)
Для получения энергии связи молекулы нужно осуществить вычисление полной энергии молекулы углерода Есс(|а|) по формуле (15) для различных расстояний между ядрами |а| и из полной энергии молекулы вычесть полную энергию двух изолированных атомов углерода 2Ес. При вычислении полной энергии по формуле (15) были использованы следующие значения параметров:
а3 = 0,8; а4 = -0,2669;
а =
0,46519262673• а 5 0 4 0 0 5 ,0 0, +
а -1,1996914627
в = -0,39292.
Параметры а1 и а5 рассчитывались по формулам (26) и (18) соответственно . На рис . 1 приведен график зависимости Е (|а|) - 2Ес в интервале от 2,15 до 2,55 радиуса Бора.
На графике зависимости Е (|а|) - 2Ес есть минимум . Вблизи от этого минимума график можно аппроксимировать параболой
Есс(|а|) - 2Ес - -0,23558 + 0,39888 • (|а| - 2,341)2 .
Хотя мы и описывали ядра атомов углерода в молекуле в адиабатическом приближении как покоящиеся классические частицы, на самом деле они участвуют в так называемых нулевых колебаниях. В этом случае должна решаться квантовая задача о колебании двух частиц вокруг положения равновесия в потенциале Е (|а|). В качестве (29) эффективной массы выступает М = 0,5 М, где Мс = 22214,8 — масса атома углерода . Энергия нулевых колебаний будет иметь вид
А
М„
• 27,211396 = 0,2306 эВ,
(30)
где k1 = 2 • 0,39888 — коэффициент жесткости молекулы углерода
График зависимости Есс(|а|) - 2Ес и аппроксимирующая парабола
Поэтому энергия связи ECB определяется формулой
ЕСВ =ЕсМ ) + — ~2Ес =
СВ cc\J Imrn / 2
0,2306
-0,23558-27,211396 +
(31)
= -6,2951 эВ, где |а|ш1п = 2,341 радиуса Бора.
Этот результат, как и значение энергий нулевых колебаний, находится в хорошем согласии с экспериментальными результатами, взятыми из работы [8] . Значение энергии связи (31) совпадает с экспериментальным значением с точностью до сотых долей эВ
Список литература
1. Береговой, А . В . Решение уравнения Кона — Шема для цилиндрических атомов методом опорной функции / А. В . Береговой, А . Г. Шкловский // Науч . ведомости Белгород . гос . ун-та. Сер . Математика. Физика. 2013. № 5 (138) . Вып. 30 . С . 121-134 .
2 . Береговой, А. В . Метод опорной функции для уравнений Шрёдингера и Кона — Шема / А . В . Береговой, А . Г Шкловский // Вести. Челяб . гос . ун-та. 2012 . № 9 (300) . Физика. Вып. 16 . С. 60-69.
3 . Береговой, А. В . Приближение локального функционала плотности с обменно-корреляционной энергией для релятивистских атомов / А. В . Береговой, А. А . Плесканев, А. Г Шкловский // Науч. ведомости Белгород, гос . ун-та. Сер . Математика. Физика. 2012 . № 23 (142), вып. 29 . С . 17-42 .
4 . Старовойтов, А . С . Модифицированный локальный потенциал для вычисления энергии ионизации атомов / А. С . Старовойтов, А. Г. Шкловский // Науч ведомости Белгород гос ун-та Сер Математика. Физика. 2010 . № 11 (82), вып. 19 . C . 126-134.
5 . Теория неоднородного электронного газа / под ред. С . Лундквиста, Н. Марча. М . : Мир, 1987.
6 Бабичев, А П Физические величины : справочник / А. П. Бабичев, И. А. Бабушкина, А . М . Братков-ский и др . М. : Энергоатомиздат, 1991.
7. Hohenberg, P. Inhomogeneous Electron Gas / P. Hohenberg, W. Kohn // Phys . Rev. B . 1964. Vol . 136 . P. 864-871.
8 . Becke, A. D . Completely numerical calculations on diatomic molecules in the local-density approximation // Phys . Rev. A — Gen . Phys . 1986. Vol . 33, № 4 . P. 2786-2788 .
