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

УДК 539.2, 539.19

РЕШЕНИЕ УРАВНЕНИЯ КОНА-ШЕМА ДЛЯ МОЛЕКУЛЫ ВОДОРОДА МЕТОДОМ ОПОРНОЙ ФУНКЦИИ А.Г. Шкловский

Белгородский государственный университет,

ул. Студенческая, 14, Белгород, 308007, Россия, e-mail: [email protected]

Аннотация. На основе метода опорной функции численно решается уравнение Кона-Шема для молекулы водорода. Предложенный автором общий метод позволяет за приемлемое время решать двумерные и трехмерные задачи, связанные с уравнением Кона-Шема. Изучена зависимость полной энергии и энергии связи молекулы водорода от расстояния между ядрами.

Ключевые слова: уравнение Кона-Шэма, молекула водорода, энергия связи, опорная функция.

1. Введение. Впервые данный метод предложен в работе [1]. В неё исследовалось уравнение Шредингера для иона молекулы водорода, так как наличие только одного электрона сводило потенциал Кона-Шема к известному потенциалу притяжения электрона к обоим ядрам.

Молекула водорода уже в течении восьмидесяти лет является эталоном для проверки всех новых методов проведения квантово-механических расчетов в молекулах. Это связано с тем, что хотя точного решения для молекулы водорода не существует, однако использование различных вариационных методов позволило получить решение практически с экспериментальной точностью. Тем не менее, все эти методы базировались на подборе удачной пробной волновой функции для системы из двух электронов с достаточно большим количеством варьируемых параметров. Конечно, различные варианты построения волновой функции, опробованные на молекуле водорода, нашли свое применение при систематическом приближенном вычислении электронных свойств молекул. Следует отметить, что прямых методов решения уравнения Кона-Шема используется очень мало, что объясняется быстрым нарастанием количества узлов в случае двух и особенно трехмерной сетки. Поэтому научная новизна настоящей работы связана с тем, что впервые применён метод опорной функции для прямого решения задачи Кона-Шема в молекуле водорода.

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

В работе [2] было выведено уравнение Кона-Шема для атома в приближении аппроксимированного локального потенциала (ПАЛП). Очевидно, что для молекулы водорода в ПАЛП это уравнение будет иметь аналогичный вид

[-Д/2 + Уто1(т)] ■ Ф1(г) = Е1Ф1(г), (1)

где

Vmol{r ) Vhm{r ) -\- Vxcp(.T )

ra П

т ^ ^е-1 п(гг) ^

укт(г = —----- (1гг, 3

N У | П — т |

к N'-1 Г п(п) (Г(п0)Г(г0) + ^(пь)Г(гь))

14ср(г) =----— J ------------рг—^-------------С1п . (4)

Здесь также использованы обозначения

га = г - Ет/2 , ть = т + Лт/2 , (5)

а Ят - расстояние между протонами в молекуле водорода. Функция Г(т), описывающая перераспределение плотности заряда за счёт процессов обмена и корреляции в ПАЛП, имеет вид

Г(х) = х 1п(ах) ■ ехр(—вх2) . (6)

В (4) функцию Г(та) отсчитываем от ядра , функцию Г (ть) отсчитываем от ядра Ь, параметр а в универсальной аппроксимирующей функции Г для молекулы водорода выбираем равным а = 1.085 — [(1.085 — 0.57)/2.2] ■ Кт. Закон сохранения полного заряда молекулы требует, чтобы функция Г(т) удовлетворяла условию

к

J т2Г(т)п(т)^т = 0 , (7)

0

что легко реализуется для заданного а подбором константы в при самосогласованном решении уравнения Кона-Шема (1).

Будем использовать в качестве опорной функции для молекулы водорода известную функцию Ф(т):

£3/2

= (Л (ехр + ехр • (8)

л/2П (1 + 5(р))

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

5(р) = (1 + р + р2/3) ■ ехр(—р), р = £Ят , (9)

Будем определять параметр £ из условия получения минимальной полной энергии молекулы

11

Ето1 = ‘2Е\ п(г!) (1//гт(г1) + Ухср{г\)) с1г\ + —— , (10)

2 ^ ^т

где Е1 - решение уравнения (1) с потенциалом (2) в основном состоянии.

Самосогласованную электронную плотность п(т) получаем из решения уравнения (1) по формуле

п(т) = 2 |Ф1(т )|2 . (11)

Множитель 2 появляется в (11) из-за спина.

Будем искать решение уравнения (1) итерационным способом. В нулевом приближении подставляем в формулу (11) вместо Ф^г) опорную функцию (8). В результате получается начальное приближение для электронной плотности

По (г)

£3

(exp(—2£га) + exp(—2£ть) + 2exp(-£ra - £ть)) . (12)

Подставим (12) в формулу (3) и получим Унт0(т):

^Ь,т,0(т )

1

— (1 - ехр(-2£г0)(1 +£г0)) +

Та

12

— (1 - ехр(-2£гь)(1 + £гь)) + -^-УніГа, п) Ть лт

(13)

Здесь введено обозначение для потенциала, созданного зарядом перекрывания, приведенное в работе [3],

Ун(Та,Тъ)

X

х \ 3{р) її 1 + ,5(-р)ЕІ (р(А + 1)) — 5(р)ЕІ (р(Л - 1)) - рехр(-рА) ) + (14)

+ 2 Г “ 3

Л£ (р) — exp(—рЛ)

Р л рЛ2

Ъ + Х + 12

В формуле (14) использованы обозначения Л = (га + гь)/Дт, р = (га — гь)/Дт, р = £■ К„ и Б(р) определены в (9), а интегральный логарифм Е1(я) определяется формулой

Бі(ж)

exp(—і)^і

Потенциал Ухср содержит кроме экспонент, которые были в электронной плотности, ещё и сложные функции ^(т). Это значит, что в аналитическом виде даже в нулевом приближении он не может быть рассчитан. Поэтому рассмотрим процедуру вычисления этого интеграла численно. Обменно-корреляционный потенциал (4) можно переписать в виде

1

Ухср(г) = -- {р {г а)ухср{г а) + Р(гЬ)пхср{гЬ))

где

[ п(гі)Г(гі0)

Ухср{га) = —агх.

(15)

(16)

9

1

Очевидно, что интеграл (16) является решением уравнения Пуассона. Переходя к цилиндрическим координатам (т, г), запишем

52^а(т,г) 52^а(т,г) ^(т,г)

Н-----------1---= -47гра(г, 2),

(17)

дт2 дг2 4т2

и0(г,с) = ьхср{га)\/г,

ра(г, Л) = '/?.0(г, 2)Р(га)ф-.

Согласно (5) та и Тъ запишутся в виде

= \/г‘2 + Кт/2)2 , гь = ^г2 + (л + Кт/2)2. (18)

Очевидно, что ра(т, г) описывает плотность заряда в молекуле водорода, перераспределённого за счёт процессов корреляции. Из (17) следует, что решив уравнение Пуассона и найдя г>а(т, г), мы легко вычислим потенциал VXCp(г) в (15). Приятным моментом является то, что на границе молекулы функции ^(та) и ^(ть) пренебрежимо малы. Следовательно, решать уравнение Пуассона (17) достаточно с нулевыми граничными условиями.

Будем решать уравнение Кона-Шема (1) методом опорной функции

Фі(г) = Ф(г) + у(г, г)!\[г-. (19)

Сначала полагаем у (г, г) = 0 и находим нулевое приближение для энергии

Ею = [ Ф(г)[—Д/2 + Ктог(г )]Ф(г )^г =

£2 1

— — — + Е00 + -Е^со

(20)

где введены обозначения £

00

1 + £ (р)

(ґ Л\ Лт ^ 1 I Л ( ^ 1 — (1 + Р) ехР (—2р)

С - 1 - 2 - 0(1 + Р) ехр(-р)---------------------------

р

(21)

£*с0 = По (г ) (^шо(г )+ ^хсР(г )) ^Г. (22)

Здесь потенциалы У^т0(г) и УХср(г) определены формулами (13) и (15). При выводе (20) было использовано уравнение

[—Д/2 + Уто1('Т)] ' Ф(т) =

2

£

— — + Уитіг) + Ухср(г)

Ф(г)+ УРт (г,г) , (23)

где

УРт (г, г)

£3/2

л/2тг (1 + 5(р))

------ - —) ехр (-£га) + ----- - — ) ехр (-£гЬ)

Га П/ \ ГЪ Га '

Уравнение (23) точное. Учитывая (19) и (23), преобразуем (1) к виду

[-А/2 + К го/(г)] • у (г, г)

£2

£-1 + — ^т(г) — Ухср(г)

Ф(г) - Урт (г, г)

(24)

(25)

(26)

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

к = гтах/ , г = к ■ п, г = Л, ■ т,

т = 1, 2,...,Ь1, п = 1, 2,...,Ь2 .

Здесь Ь1 - количество узлов по оси г, Ь2 - количество узлов по оси г, гтах - поперечный размер молекулы, для расчёта его выбираем равным 12 радиусам Бора.

В силу симметрии молекулы достаточно проводить расчет только для положительных г < гт = кЬ2, так как для г < 0 функция у (г, г) та же самая, что и для г > 0. Записываем (25) в виде разностного уравнения:

ут+1,п + ут-1,п + ут,п+1 + ут,п-1 4ут,п

к2

+

+

(27)

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

(рр1)т,п = 2\//Г-_т

Пусть

п

Ею + ^ (Унто)т,п (У*ср)т^п ) Фт,»г (Урт.)т.,п

4-

4т2

- 2к2 ^(Уто1)т,п - Ею)

-1

(28)

(29)

Воспользуемся следующей итерационной процедурой для нахождения неизвестной функции Ут,п:

ут,п ^т,п ут+1,п + ут-1,п + ут,п+1 + ут,п-1 + к (ЕЕ 1)г

(30)

2

1

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

значения ут,п. Процедура итерирования продолжается до достижения нужной точности. Так как это только нахождение первого приближения, то требуется примерно тысяча итераций. После получения первого приближения мы должны уточнить энергию Е10 и потенциалы Унт и Ухср, которые вычислялись с использованием плотности п0(г). Для этого вычисляется функция первого приближения Ф[1](г) = Ф(г) +г/1^(г, с) / у/г , в которой у(1)(г, г) вычислено по формулам (28)-(30). Далее по формуле (11) находим электронную плотность молекулы в первом приближении, проверяем её нормировку и перенормируем найденную итерационным образом функцию у(1)(г, г) так, чтобы полное количество электронов в молекуле водорода равнялось двум. При этом должно выполняться условие

\frdr (2Ф (г, z) + yW (г, л)) yW (г, z) dz = 0 .

(31)

Оно является критерием качества численного расчёта с использованием итерационной процедуры и перенормировки. Если этот критерий хорошо выполняется, можно уточнить энергию первого приближения Е11 по формуле

р ( \ _ Е00 - Еро (п) С2 1 Еи (п) г + ^ {п) - - + - (Ehc 1 - Ehc0)

(32)

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

Ipi (/?) = J \frdr J у^ (г, л) • Ф(г, z)dz , 0 0

(33)

Г

z

Г

z

m

Г max zm

Epo (/?) = J y/rdr J y{n) (r, z) ([Vhm(r, z) + Vxcp(r, z)] Ф(г, л) + Урт (г, л)) • Ф(г, z) dz . 00

(34)

В формулах (32-34) индекс n показывает число произведенных итераций с перенормировкой.

На рис. 1 представлена поправка к опорной функции, полученная после двадцатой итерации с перенормировкой. Видно, что отличные от нуля значения поправки сосредоточены в области примерно 6 б.р. При этом максимальное отклонение от опорной функции вблизи нуля около 0,06.

Рис. 1. Поправка к опорной функции у(20) (г, г), рассчитанная с учетом условия нормировки (31).

На рис. 2 изображена электронная плотность молекулы водорода, вычисленная с учетом поправки, представленной на рис. 1. Визуально рис. 2 практически не отличается от рисунка, где та же плотность была получена только с использованием опорной функции. Это как раз связано с малой величиной поправки у(20) (г, г). В то же время, несмотря на малую величину поправки, она важна при вычислении энергии связи молекулы водорода. Её учет приводит к существенному улучшению согласия рассчитанной энергии с экспериментальным значением.

Перейдем к рассмотрению полной энергии и энергии связи молекулы водорода. Запишем полную энергию молекулы водорода в виде

Е = 2Е1 + -^~] [ [ (1 - Р (гі а) Р (г2о) - Е (г1Ь) Е (г2Ь)) с1псВ2 . (35)

Лт 4 3 3 |гі - Г2І

Функция Е (г), аппроксимирующая обменно-корреляционную энергию в ПАЛП, даётся формулой (6), Е1 - энергия, полученная в результате решения уравнения Кона-Шема (1).

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

Рис. 3. Значения полной энергии Е(Кт) молекулы Щ в атомных единицах в зависимости от расстояния Кт [б.р.].

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

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

Для контроля правильности сделанного упрощения была дважды пересчитана обменно-корреляционная часть энергии. Обменно-корреляционная энергия считалась сначала для электронной плотности, построенной только из опорной функции, а второй раз с уточненной электронной плотностью. Тем не менее, уточнение электронной плотности привело к изменению энергии только в третьем знаке после запятой. Поэтому для полной энергии будем использовать формулу

1 1 /*

Е(Ит) = 2Е1 + —-------- п (г2) (г2 ) с1г2 , (36)

^т 4 о

УМЭ = УктО (г ) + Ухср (г ) • (37)

Здесь Унт0 (г) дается формулами (13) и (14), а вычисление потенциала Ухср (г ) подробно описывается формулами (15)-(17).

При решении уравнения Кона-Шема для фиксированного расстояния между ядрами параметры а и в задаются, и соответствующие потенциалы рассчитываются в виде двумерных массивов, которые могут быть использованы при численном вычислении поправки к полной энергии. Хотя при любом значении параметра £ точное нахождение поправки к опорной функции должно давать одну и ту же электронную плотность, а, следовательно, и одну и ту же энергию Кона-Шема, при реальном численном расчете проявляется зависимость от этого параметра вследствие того, что расчёт осуществляется всегда с конечным шагом. В настоящей работе выбран шаг к = 0, 03 б.р. Поэтому поправка не может быть получена с абсолютной точностью и электронная плотность будет несколько неточной. Согласно теореме Кона, энергия будет достигать минимума на правильной электронной плотности. В данной работе были проведены расчеты для разных параметров £, и при построении графика на рис. 3 выбиралось такое значение, при котором полученная энергия £1 была бы минимальной. В табл. 1 приведены расчетные данные для полной энергии £(Лт).

Таблица 1

Расчётные данные для полной энергии Е(Кт) молекулы И2 в атомных единицах

а = 1.085— [(1.085 —0.57)/2.2] -і?т £ Е(Ет)

1,5 0,733864 1,350 -1,166647

1,374 -1,166830

1,400 -1,166610

1,48 0,738545 1,350 -1,167000

1,375 -1,167292

1,400 -1,167183

1,46 0,743227 1,375 -1,167642

1,385 -1,167691

1,400 -1,167644

1,44 0,747910 1,385 -1,167973

1,395 -1,168000

1,405 -1,167970

1,42 0,752591 1,395 -1,168205

1,405 -1,168218

1,415 -1,168169

1,40 0,757273 1,400 -1,168310

1,410 -1,168340

1,420 -1,168307

1,39 0,759613 1,400 -1,168300

1,410 -1,168354

1,420 -1,168345

1,38 0,761954 1,410 -1,168324

1,420 -1,168340

1,430 -1,168294

1,37 0,764295 1,400 -1,168148

1,420 -1,168289

1,440 -1,168182

Ет - расстояние между протонами в молекуле Н2, а - параметр в обменно-корреляционном потенциале (6), £ - параметр в опорной функции (8). Жирным шрифтом выделены минимальные значения энергии, которые использовались для построения рис. 3.

Рис. 4. Функция Е(Кт) и её разложение в ряд Тейлора вокруг точки Ет = 1, 38961 б.р. в интервале от 1,38 до 1,40 б.р.

На рис. 4 в интервале от 1,38 до 1,40 б.р. показана функция Е(Ят)и её разложение в ряд Тейлора вокруг точки Ят = 1, 38961б.р. Видно, что совпадения кривых в области вблизи минимума полной энергии очень хорошее.

Для определения частоты нулевых колебаний на рис. 4 построена кривая Е(Ят) в интервале от 1,38 до 1,40 б.р. и её разложение вокруг точки Ят = 1, 38961 б.р. в ряд Тейлора. Эту кривую можно рассматривать как потенциальную энергию взаимодействия протонов между собой в зависимости от расстояния между ними. На рис. 4 вблизи минимума кривая потенциальной энергии разложена в ряд до второго порядка по малым отклонениям от положения равновесия. Протоны, находящиеся в потенциале Е(Ят), совершают нулевые колебания. Они полностью аналогичны нулевым колебаниям протонов в ионе молекулы водорода, которые были рассмотрены в работе [1].

Будем считать, что в хорошем приближении колебания протонов вокруг положения равновесия хорошо описываются в гармоническом приближении, как это следует из рис. 4. Тогда можно воспользоваться приближенным выражением

Е(Я,„) = -1,168354 + 0,135335 ■ (х - 1,3896)2

(38)

В этой задаче двух тел в качестве эффективной массы будет выступать половина массы протона тмэ = тр/2 =918,075 а.е. Энергия нулевых колебаний в атомной системе единиц будет иметь вид

к

1

тмд

(39)

где к1 - коэффициент жесткости молекулы водорода. Он определяется при разложении полной энергии (36) в ряд Тейлора вблизи положения равновесия Кт = Ят0 = 1, 3896

Ам = а2Е1Ет) • (40)

&гп — ^т0

Если воспользоваться приближенным выражением (38) то, очевидно, что к1 = 0,135335■

2 и частота нулевых колебаний ^0 рассчитывается по формуле (39):

0,0 = 27, 2117^0,9183075~ = °’ 467228 эВ ' (41)

Энергия нулевых колебаний равняется половине величины ^0, т.е. 0,234 эВ. Экспериментальное значение этой величины 0,27 эВ [3], что очень хорошо согласуется с полученным результатом.

Для получения энергии связи из полной энергии молекулы нужно отнять полную энергию двух изолированных атомов водорода. В атомной системе единиц сумма энергий двух изолированных атомов водорода равняется -1 Хартри. Поэтому энергия связи определяется формулой

Есв = Е(Кт0) + 1 + — . (42)

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

Есв = 27, 2117 • (-1,168354 + 1) + ^ = -4, 35 эВ . (43)

Экспериментальное значение энергии связи [3] равно -4,45 эВ, что на 0,1 эВ меньше полученного в (43) значения. Эта погрешность может быть связана как с недостаточно малым шагом решения к =0,03 б.р., так и с тем, что не учитывалось изменение электронной плотности для пересчёта потенциала Хартри. С другой стороны, полученный результат свидетельствует о высокой точности использованных численных методов, так как рассчитывалась не энергия связи, а полная энергия молекулы водорода, равная -31,9 эВ. Полученный результат (41) даёт величину -31,79 эВ, что отличается от правильного результата менее чем на 0,3%.

Автор благодарен Шкловской Н.А. за помощь в численных расчётах.

Литература

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

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

3. Слетер Дж. Электронная структура молекул / Дж. Слетер. - М.: Мир, 1965 -588 с.

THE KOHN-SHEM EQUATION SOLUTION OF HYDROGEN MOLECULE ON THE BASIC FUNCTION METHOD A.G. Shklovskiy

Belgorod State University,

Studencheskaya St., 14, Belgorod, 308007, Russia, e-mail: [email protected]

Abstract. This article describes Solution of the Kohn-Shem equation which is described the hydrogen molecule has been found using the basic function method. It provides to solve two- and three-dimensional problems numerically during the real time. It is investigated the dependence of the total energy and the hydrogen molecule bound energy on the distance between nuclei.

Key words: Kohn-Shem equation, hydrogen molecule, bound energy, basic function.