CC BY
52
УДК 519.673, 539.182
РЕШЕНИЕ УРАВНЕНИЯ КОНА-ШЕМА ДЛЯ МОЛЕКУЛЫ МЕТОДОМ ОПОРНОЙ ФУНКЦИИ
А.Г. Шкловский1^ М.А. Шкловская2)
^ Белгородский государственный университет, 308007, г. Белгород, ул. Студенческая, 14 2) Московский Государственный Технический университет им. Н.Э. Баумана,
105005, г. Москва, ул. Вторая Бауманская, 5, e-mail: [email protected]
В данной работе предлагается обобщение предложенной ранее аппроксимации обменно-корреляционной энергии в теории функционала электронной плотности. Это обобщение используется для нахождения обменно-корреляционного потенциала двухатомной молекулы. Предложен новый алгоритм решения уравнения Кона-Шема в молекуле - метод опорной функции. Для апробации метода проведен расчет энергии диссоциации иона молекулы водорода.
Ключевые слова: функционал электронной плотности, обменно-корреляционный потенциал, уравнение Кона-Шема, метод опорной функции.
Введение
В работе [1] описывается метод, разработанный для того, чтобы оптимальным образом аппроксимировать обменно-корреляционную энергию в атоме. Этот метод будет использован в настоящей работе для аппроксимации обменно-корреляционной энергии и обменно-корреляционного потенциала в молекуле.
Предлагается новый алгоритм решения уравнения Кона-Шема в молекуле -метод
опорной функции. Явный вид опорных функций // (г) выбирается в зависимости от конкретной молекулы так, чтобы она была известна с высокой точностью, нормирована и мало отличалась от молекулярной волновой функции. Обычно удобно брать в качестве опорной функции собственные функции соответствующих сжатых сферических атомов в молекуле.
Для проверки точности метода опорной функции была рассмотрена тестовая задача, имеющая точное решение. Энергия диссоциации иона молекулы водорода известна и есть независимые методы решения уравнения Шредингера для электрона в
поле двух протонов, находящихся на заданном расстоянии Rm. Поэтому предложенный новый алгоритм применяется к вычислению энергии диссоциации иона молекулы водорода.
Далее будет использоваться атомная система единиц h = е2 = me = 1.
1. Обменно-корреляционный потенциал для молекулы
В работе [2] приведена формула для вычисления полной энергии в рамках теории функционала электронной плотности:
Е = ТеЫ + И + EJn]+JVe(r)• n(r)dr, (1.1)
где Te [n ] - функционал кинетической энергии, Е har [ n ] - функционал Хартри:
г [ ] 1 f n(Tj)n(Я,)
har [ n ] = — J —F-------ry~ dr1 dr2 (12)
2 J |Г1 - Г2 I
Ve (r) - кулоновский потенциал взаимодействия электрона с ядрами, а обменно-
корреляционная энергия Е хс [ п ] дается выражением:
1 С п (к (^ г2)п (г2)
Г1 - Г2
<г1 <г2
Здесь
I
к (^ Г2) = | <Л( gx( ^ Г2) -1),
(1.3)
(1.4)
где gл( Г1, г2 ) - парная корреляционная функция для системы с электрон - электронным взаимодействием, измененным в X раз, где X пробегает значения от 0 (для не взаимодействующей системы электронов) до 1 (электрон - электронное взаимодействие включено полностью).
Формулы (1.1 -1.4) точные, но функция к (г1, г2) не известна. В приближении локального потенциала (ПЛП) вводится обменно-корреляционная энергия, приходящаяся на одну частицу £ хс (п ( г )) , и полная обменно-корреляционная энергия (1.3) записывается в виде:
Ехс [п] = \£хс (п(Г X Г) •п(г )dr . (1.5)
В работе [1] предложена аппроксимация функции к(^,г2) в выражении (1.3) и приведен вывод уравнения Кона - Шема для атома в приближении аппроксимированного локального потенциала (ПАЛП). Для молекулы минимум функционала (1.1) ищется при условии
г
] П(Г )<Г = Ыет , П(Г ) = 2 фа (Г)
а=1
(16)
Здесь п(г ) - плотность электронов в молекуле, Ыет - число электронов в молекуле, а функция Фа(г) является решением уравнения Кона - Шема:
[ - (1/2) А + Умоя (г)] •Фа (г) = Е ,Ф. (г). (1.7)
Как обычно, считается, что уравнение Кона - Шема описывает систему невзаимодействующих квазичастиц с энергией Еа, находящихся во внешнем потенциале:
V (г) = У (г) + У . (г)
мод\ у тхс V / ткаг V /
Я 2
г + Я 2
(18)
Здесь для потенциала Хартри используется точное выражение:
п(г)
V
ткаг
(Г) = | <2
г - г
(19)
Электронная плотность квазичастиц совпадает с молекулярной плотностью п ( г ) и может быть представлена в виде:
п(Г) = п1 (|Г - Ят /2|) + п2 (
г + Я /2
) + £п(г ).
(110)
0
2
Здесь Ят - вектор, идущий от ядра первого атома до ядра второго атома, а 8п ( г ) показывает отличие молекулярной плотности п (г ) от суммы плотностей
сферических сжатых атомов в молекуле п1 (|г — Ят / 2 ) и п2 (|г + Ят /2 ) и находится
самосогласованным образом.
Чтобы найти электронную плотность молекулы по формуле (1.6) необходимо вычислить волновую функцию Фа (г ) :
Ф а ( г ) = ¥ а ( г ) + а ( г ). (1.11)
Здесь ¥а (г) описывает отличие молекулярной волновой функции Фа(г ) от опорной функции (г) . Эти опорные функции используются для построения
плотностей пД г ) и п2 (г ). Также как и для электронной плотности, поправку к
опорной функции необходимо находить самосогласованным образом.
Подставляя выражение (1.10) в формулу (1.9), получим потенциал Хартри в молекуле в виде:
дп(^2)
Vmhar <7 ) = J dr2 + V1har (|Г - Rm / 2|) + V2har (|Г + Rm / 2|) • (1.12)
Здесь введены обозначения для потенциалов Хартри сферических атомов 1 и 2:
V1har (r ) = — (1 - Ui(r)), U 1(г ) = -^f x2ni(X)i1 - -Ът. (1.13)
r Nel \ V X J
N 4 Я 2 z' \
V2har (r) = —i- (1 - U 2 (r ))’ U 2(r) = -Г—• f X 2 n2(x)[1 - ^ Ц (114)
r Ne2 i V x J
Обменно-корреляционный потенциал молекулы Vmxc (r ) точно не известен. Для
молекулы, также как и для атома берётся выражение:
(-) f л- й(Г) • h (r, r)
Vmxc (r ) = J dr2---2 - 2 •
J r - r2
r 2 (115)
Для электронной плотности n(r ) используется формула (1.10).
При аппроксимации функции h (r1, r2) для молекулы учтем, что она должна быть симметрична
Щ,г2) = h(rL rd, (1.16)
на больших расстояниях от ядер должна обращаться в ноль
1-m h ( -2) = 0, (1.17)
|Г1|,|г2 œ ’ v 7
и для r , лежащих внутри молекулы, удовлетворять правилу сумм [2] :
J ¿А п(Я,) • h (Г, Я,) = -1. (118)
В методе опорной функции для нахождения электронной плотности ôn(r ) в мо-
лекуле удобно пользоваться цилиндрической системой координат. Рассмотрим снача-
ла опорную функцию ¥а (г ), взятую из решения уравнения Кона - Шема для сфери-чески-симметричного первого атома. Для краткости будем называть их решениями первого типа. В случае решения первого типа вместо условия (1.18) будем использовать:
\ <2 п1(г2) • к (Г, г) = -1,
J'
J<
(119)
<2 п2(г2) • к (г, г2) = 0 <т2бп(г2) • к (г, г) = 0.
Воспользуемся формулой (1.10) для электронной плотности п(г) и перепишем выражение (1.15) в виде суммы вкладов от отдельных интегралов:
(Г ) = | <г2
п 2 ( г2) • к (г > г2) , Г 8п (г2) • к ( г , г2) • (120)
n1 (r2) • h (r, r2)
(r ) = I dr 2 ГЦ—— +
mxc
r2
f dr2 П 2( ^ ) ^ h (r , r2 ) + f df2
J 2 IV — Vi i 2
r - r2 r - r2
В первом интеграле, который берётся по шару с центром в ядре первого атома, аппроксимация функции h (r1, r2) очевидна:
- f 1 ^
h (^ r2) = - — +П( ^ r2) •© (R1 - ИЬ •© (R1 - \r2\)- (121)
V Ne1 J
Здесь Ne1 - число электронов в первом атоме, Ri - радиус этого атома, Í1, x > 0,
0(х) = < - функция, обеспечивающая наличие обмена и корреляции только
[0, х < 0
внутри атома, а n(r1,r2) - функция, описывающая перераспределение плотности заряда
за счёт процессов обмена и корреляции. Закон сохранения перераспределенного заряда (1.19) имеет вид:
J dr2n1(r)n(r1, r2) = 0- (1.22)
В [1] для сферически симметричного атома n1(r) = n1(r) и функция п1(Г1, r2) выбирается в виде:
П(Г r2) = -Ne1 -1) • F1(r1) • F1(ri), (1.23)
где аппроксимирующая функция Fi(r) зависит только от модуля r и не зависит от углов. Универсальная аппроксимирующая функция F(r), описывающая перераспределение плотности заряда за счёт процессов обмена и корреляции в ПАЛП, имеет вид:
F(x) = х • ln(jx) • exp(-Дх2). (1.24)
Здесь введены две константы. Константа у определяет точку, в которой натуральный логарифм, а, следовательно, и функция F(r), меняет знак. Равенство (1.23) гарантирует выполнение свойства (1.16). При любом положительном в выполняется свойство (1.17). Множитель (Ne1 -1) в (1.23) связан с тем, что учитывается средний вклад в перераспределение зарядовой плотности от каждого электрона, участвующего
в процессах обмена и корреляции. Закон сохранения перераспределенного заряда (1.22) требует, чтобы функция ^(г) удовлетворяла условию
Я1
| г2 F1(r)n1(r ")Сг = 0, (1.25)
0
что легко реализуется для заданного у подбором константы в при самосогласованном решении уравнения Кона - Шема для первого атома.
Во втором интеграле, который берётся по шару с центром в ядре второго атома, аппроксимация функции к (г1, г2) также очевидна. Из второго интеграла в (1.19) следует,
что к(г1,г2) = -Nе2 • F2(r1)• F2(r2), где F2(г) берётся для для второго атома. Поправка 5п ( г ) в формулах (1.20), (1.19) достаточно мала в тех областях интегрирования, где величины пС (|г - гс |) и пО (|г - г0 |) являются основными. Так как заряд, связанный с величиной 5п (г ), равен нулю и третий интеграл в (1.19) тоже равен нулю, то в третьем интеграле в (1.20) заменим функцию к(гх,г2) на её среднее значение по молекуле,
равное 1/№ет. Это соответствует исключению из потенциала Хартри усреднённого са-модействия.
Решения уравнения Кона - Шема с опорной функцией, взятой из решения для второго атома, будем называть решениями второго типа. В случае решений второго типа вместо условий (1.19) будем полагать:
Ог2 п1(г2) • к (г,г) =
/■
I
I
СЩП2 (г2 ) • к (г, г2) = -1, (1>26)
<іг2дп(г2) • к (г, г2) = 0.
В этом случае в первом интеграле в (1.20), который берётся по шару с центром вокруг первого атома, аппроксимация функции к (гх, г2) с учётом (1.26) очевидна. Из первого интеграла в (1.26) следует, что к (^, г2) = -Ые1 • Fl(гl) • ^(г2), где ^(г) берётся для первого атома. Во втором интеграле, который берётся по шару с центром вокруг
второго атома, аппроксимация функции к (г, г2 ) также очевидна. Она даётся формулами (1.21), (1.23) и (1.24) с константами у, в и N е 2 для второго атома. Поправка 5п (г) в формулах (1.20), (1.26) достаточно мала в тех областях интегрирования, где
величины пх ( г - Ят /2 ) и п2 (|г + Ят / 2|) являются основными. Так как заряд, связанный с величиной 8п (г ), равен нулю и третий интеграл в (1.26) тоже равен нулю, то в
третьем интеграле в (1.20) заменим функцию к (^, г2) на её среднее значение по молекуле равное 1/№ет. Это опять соответствует исключению из потенциала Хартри усреднённого самодействия.
Теперь можно переписать приближенный молекулярный потенциал для уравнения Кона - Шема первого типа в виде:
V,», (Г ) = У„,1 (|Г - ‘т / 2|) + Уа12 (|г + Ят / 2|) + | ІГ (1.27)
N і г - г
1Ует V 2
Здесь введены обозначения:
(^ -1) •(! - Щг) + к,с,1(г))- 21
^1(г )
V
хср1
Я
(г) = 4^F1( г) I х2 F1( х) • п1( х) -| 1 - —
г ^
Ne2 •(! - и2 (г) + Ухср2 (г))- 2
Уаі2 (г) =
х у
йх
Я2
Кср2 (г) = 4пЪ2(г) I х2Ъ2(х) • п2(х)
'1 -
V
х
йх
(1.28)
(129)
(130) (1.31)
Аналогично можно переписать аппроксимированный молекулярный потенциал для уравнения Кона - Шема второго типа в виде (1.27), в котором введены обозначения:
N1 •(! - и,(г) + УХСрЛ(г))- 7,
Г
V 2 (г) =
(Ne2 - 1)^(1 - ^2 (г) + Ухср2 (г))- 7
Г
(132)
(1.33)
2. Решение уравнения Кона-Шема для молекулы методом опорной функции
Используя метод опорной функции, следует переписать уравнение Кона-Шема
для 8ц/а{г ), выделив /]а (г )
первого типа:
[ - (1/2) А + V,», ( г )] •^Ра (г ) = Еа-б^а (? ) +
(Еа - Є) - Уаі2 (
г + Я /2
)-^к(г)) •/а(г)
Здесь введено обозначение:
ёК(г) = (^ 1) Iіг2 бп(г)
N „
г - г
(2.1)
(2.2)
а е1 - собственное значение, соответствующее опорной функции / (г) первого типа,
полученной при оптимальном решении задачи о первом сферическом сжатом атоме в молекуле.
Аналогично получим уравнение Кона-Шема для 5/а(г) , выделив опорную
функцию /Ох (г ) второго типа:
[ - (1/2) А + V,», ( г )] • бРа ( г ) = Еа-ёР а ( ? ) +
(Еа - Є) - Уаі 1(
г - Я /2
)-бк (г)) •/о (г)
(2.3)
где £;2 - собственное значение, соответствующее опорной функции // (г) второго типа, полученной при оптимальном решении задачи о втором сферическом сжатом атоме в молекуле.
г
г
г
Если получено решение уравнений (2.1) и (2.3) для всех а, то можно вычислить поправку к молекулярной электронной плотности 5п (г ) :
5п (Г ) = Х Па-^а ( г ) ■ (Жа ( Г ) + 2-^ ( Г )) . (2.4)
а
Как было указано выше, поправка (2.4) находится самосогласованным образом. Сначала полагаем ее равной нулю. В этом случае из формулы (2.2) следует, что (5К(г) тоже равно нулю. Решая приближенные уравнения (2.1) и (2.3) для всех а, вычисляем
поправку к молекулярной электронной плотности 5п (г)
. Подставляя полученную
поправку в формулу (2.2) получаем новое приближение для 5п (г) . Процедуру продолжаем до тех пор, пока начальное и конечное 5п ( г ) не совпадут с заданной точностью.
В уравнениях (2.1) и (2.3) удобно перейти в цилиндрическую систему координат с осью 2 направленной вдоль линии, соединяющей ядра атомов в молекуле. Тогда потенциалы (1.27)—(1.33) не зависят от угла ф. Для уравнения (2.1) удобно ввести обозначение для функции 5/а (г ), не зависящей от угла ф:
/(г) = Уа(г, 2)/^г , (2.5)
и переписать уравнение (2.1) в виде:
d y (r z) д y (r z) / /9
Л. V , ) + Усе ) + (0.25/r
(2.6)
dr2 dz2
-2• [V (r)-Е ])• y (r,z) = -F1 (r,z)
L мол V / a J / s а \5/ а V?/
Здесь введено обозначение:
F1 а (r, z) = 2 ■ ^(E а - SJ - V2 (|r + Rm 1 2^ - SV(r, z)) i (r > z) • (2.7)
Уравнение (2.3) также можно записать аналогично (2.6), если ввести вместо F1 а(r, z) обозначение:
F2 а (r, z) = 2 ■ ^(E а - Sj - Vat1 (|Г - Rm 1 2^ - SV(r, z)) a (r, z)- (2.8)
Для опорных функций, являющихся решениями задачи для сферического атома с сферическими гармониками, зависящими от угла ф, следует вместо (2.5) ввести 5ща (r) с
такой же зависимостью от угла ф. Например, для p состояния можно ввести 5ща (r) в виде:
0¥а (Г) = Уа (r> z) ‘ C0S (P/Jr , (2.9)
или
0¥а (Г) = Уа (r, z) ‘ Sln P/Jr . (210)
Подставляя (2.9) или (2.10) в уравнение (2.3), вместо (2.6) получим уравнение:
~ ,_/--((m2 -0.25)/r2 +
. (2.11)
d 2 У а ( r, z ) . d 2 -У а ( r, z ) dr2 dz2
+2 ■ [Кол (r ) - Е а ])‘ У а (r, z) = -F1 а (r, z)
Под величиной ¥ За (г, г) в (2.7) в этом случае следует понимать часть атомной
функции, не зависящую от угла ф.
Очевидно, что подстановка (2.10) или (2.9) в уравнение второго типа снова приведет к уравнению (2.11) с функцией Е2 а(г,г), где под величиной ¥]а(г,г) в (2.8) опять следует понимать часть атомной функции, не зависящую от угла ф. Следовательно, для опорной функции р типа решение у а (г, г) соответствует дважды вырожденной энергии. Аналогично, для числа т, не равного нулю, возникнет двукратное вырождение из-за зависимости от угла ф для каждой опорной фукнции.
Рассмотрим решение уравнений (2.6) и (2.11). Это уравнения одного типа. Все они удовлетворяют нулевым граничным условиям:
У а (г> 2шах ) = У а (Г>-2шах ) = У а (гшах >2) = У а (0,2) = 0 . (2.12)
Поэтому метод решения для них один и тот же.
Будем решать методом последовательных приближений. Перейдем от дифференциальных уравнений к разностным. Введем одинаковый шаг по осям г и z:
Г
к = ^^, г = -г + к ■ п2, г = к ■ п1,
Ы шах .
п1 = 1,2,..., Д п2 = 1,2,..., Ь2
Тогда, с точностью до второго порядка по h уравнение (2.6) примет вид:
Уа + Уа + Уа + Уа - 4 ■ уа +(0.25/п12 -
^ п 1+1 ,п 2 ^ п1-1,п 2 ^ п1,п 2+1 ^ п1,п 2-1 ^ п1,п 2 V '
(2.13)
-2 • [(Умол ) n1,n 2 - Е а ] • h 2 )• У П1,п 2 _ - ( F1 а ) n1,n 2 ' h 2
(2.14)
Так как сама величина yn\,n2 мала по сравнению с Va (r ) , то такая точность приемлема. Введем обозначение:
DnU2 _ 4-(0.25/n12 -2 Wmo,)n1,n2 -Еа] • И) . (2.15)
Теперь уравнение (3.41) можно переписать в виде:
а _ ( Уп 1, n 2 +1 + yn1, n 2-1 + У n1+1, n 2 + У nl-l, n 2) + h ' ( F 1 a ) n1, n 2 У n\,n 2 d " (216)
n 1, n 2
Вид (2.16) удобен для применения итерационного метода решения. В (2.15) входит неизвестная энергия Еа. Так как используется метод последовательных приближений, то при нахождении первоначального приближения полагаем SV(г ) _ 0 и
Еа -£j _ jvj ( Г ) Va, 2 (|Г + Rm |)vj (Г № . (2.17)
Это соответствует применению первого порядка теории возмущений к потенциалу vat2 ( r + Rm ) для уравнения (2.6). При решении уравнения второго типа среднее значение по невозмущенным функциям берется от потенциала vat1 ( r - Rm ). Остальные изменения в (2.16) при переходе от (2.6) к (2.11) или к аналогичному уравнению второго типа очевидны.
Рассмотрим формулу (2.4) для добавочной молекулярной плотности Sn(r, z) . Перепишем её с учетом обозначений (2.5):
»(г, г) = ,
Г
М(г,г) = Х па • У а (г^г) ■ (У а (г>2) + Ха (г> 2Яг ■ Аа (г>г)У (2'18)
Здесь использованы обозначения Аа(г,г) для не зависящей от ф части сферической функции, деленной на л/г2 + г2 , и ха (г, г) для нормированной радиальной атомной волновой функции. Как уже отмечалось выше при вычислении электронной плотности молекулы зависимость от ф, даже если она есть в волновой функции, исчезает. Функция ¡Жт, г) может быть табулирована в массив, если будут вычислены все функции У а (гг).
Вычисление функций Уа(г,г) проводится самосогласованным образом. На первом этапе полагается дУ (г) = 0 . После вычисления функций Уа (г, г) в первом приближении, по формуле (2.18) вычисляется добавочная молекулярная электронная плотность 8п (г, г) в первом приближении. После этого по формуле (2.2) вычисляется дУ (г) тоже в первом приближении.
Когда задан массив дУ (г, г) , то описанным выше методом могут быть найдены массивы Уа(г,г) в следующем приближении. Самосогласование будет производиться до тех пор, пока массивы Уа(г,г) не совпадут с заданной точностью. Трехмер-
ный интеграл (2.2), определяющий массив дУ (г, г), является несобственным интегралом, поэтому обычные алгоритмы для его вычисления нельзя использовать. Поэтому будем вычислять несобственный интеграл (2.2) только для границ: г=0, г=гшах, z=zшax , z= ^шах. Так как на каждой границе число узлов не велико, то общее вычисление значений интеграла во всех узлах не займет много времени. Когда граничные значения вычислены, можно легко определить значения потенциала внутри области -гш< z <гш, 0<г<гш. Для этого достаточно в цилиндрических координатах решить двумерное уравнение Пуассона:
д2 Vd (г, г) д2 Vd (г, г) 0.25 ¡-
----—2-----+-------—2-----+ —Т~ ■ vd (г, г) = -4п ■ Ж(г, г)^г . (2.19)
дг дг г
Здесь введено обозначение:
vd (г, г) = дУ (г, г) ■л/Т . (2.20)
Алгоритм решения уравнения (2.19) аналогичен тому, который был использован при решении уравнения (2.6).
3. Вычисление энергии связи молекулы
В первом разделе данной статьи обсуждалась аппроксимация обменнокорреляционной энергии и потенциала в молекуле. Будем использовать следующее оп-
ределение для энергии связи:
Есв = Е - Е2 - ^ (31)
где Е1 и Е2 - полные энергии атомов, составляющих молекулу.
В рамках метода функционала электронной плотности для полной энергии молекулы Е вместо (1.1) имеем следующее выражение:
Е = Е Па ' Е а - | Уол (Т') 'П(Г№ + Е!аг [п] +
а
г z • Z
Ес [”]+| V(г) •п(г )dr +"JRJ'
(3.2)
В первом разделе данной статьи была введена аппроксимация обменнокорреляционного потенциала в молекуле. В рамках этой аппроксимации можно записать сумму энергии Хартри (1.2) и обменно-корреляционной энергии (1.3) в виде:
корреляционного потенциала (1.20) и потенциала Хартри в молекуле (1.12). Как и ра-
зоваться явным видом потенциалов Умол (г) и Уе (г) и привести подобные члены, то формулу (3.2) можно преобразовать к виду:
В результате подстановки электронной плотности (1.10) в формулу (3.4) возникают три интеграла. Первый интеграл берется по объёму первого сферического атома, второй интеграл берется по объему второго сферического атома, а третий - от малой добавки электронной плотности, по всему объёму молекулы. Поэтому в первом интеграле для обменно-корреляционного потенциала (1.20) и потенциала Хартри в молекуле (1.12) используем ту аппроксимацию, которая была применена при решении уравнения Кона - Шема для волновых функций первого типа. Во втором интеграле для обменно-корреляционного потенциала (1.20) и потенциала Хартри в молекуле (1.12) используем ту аппроксимацию, которая была применена при решении уравнения Кона -Шема для волновых функций второго типа. В третьем интеграле для обменнокорреляционного потенциала используем простейшую аппроксимацию для компенсации усредненного самодействия в потенциале Хартри.
В результате подстановки всех необходимых выражений и приведения подобных членов в формуле (3.4), имеем:
Здесь использованы обозначения утхс (г / и утиаг(г) для обменно-
нее, для электронной плотности п(г2) используется выражение (1.10). Если восполь-
(3.4)
(3.5)
]
Здесь введены обозначения:
Е..
и
-1)
N.
1 Г
I ( I----------------------\
| йЪ20Ф2 и(г2,г2)-\ут<г(^г2 + (22 +К/2)2)+
\2
(3.6)
+У21,Г (>/ Г2 + (22 - К/2)2 ) + V (Г2> ))
Е1 И Е2 - полные энергии сферических сжатых атомов в молекуле с опорными функциями у/а (г) первого и, соответственно, второго типа, полученными для оптимального решения задачи Кона-Шема, и Епер - энергия перекрытия, которая задается
выражением:
Е = -
пер 2
V
К
П2 ( Г1)-(и 1 ( гі + К ) - У1хср (|Г + Кт )),
2 ( бТ ( гі - гі + Кт Кт )- ) ^2 хср ( Г - Кт 1 ))!
( Г - Кт ) і
(3.7)
+^ 21 (Г1
Подставляя Е из (3.5) в формулу (3.1), получим для энергии связи окончательное выражение:
а
К
К
У пІ ■ + (Е- + (Е2 + Епер - Еи
. (3.8)
Здесь введены обозначения — Е2 — Е2 и ,Е^ — £1 £1 для разности
полных энергий сферических сжатых атомов в молекуле и соответствующих свободных атомов.
Отдельные компоненты, входящие в (3.8) зависят от выбранного радиуса сжатого атома. Для того, чтобы определить, какие радиусы R1 и R2 являются оптимальными для численных расчетов энергии связи следует провести расчеты энергии связи по формуле (3.8) для разных R1 и R2 и посмотреть как от этого будет зависеть результат.
4. Энергия связи иона молекулы водорода
Для проверки точности метода опорной функции была рассмотрена тестовая задача, имеющая точное решение. Энергия диссоциации иона молекулы водорода известна и есть независимые методы решения уравнения Шредингера для электрона в
поле двух протонов, находящихся на заданном расстоянии Ят. Поэтому можно применить описанную выше теорию к вычислению этой энергии. Для энергии диссоциации Е, воспользуемся формулой:
-1
т
Еа — Есв + 412 (4.1)
'а
Здесь 0)к - частота колебаний протонов в ионе молекулы водорода, которая вычисляется в приближении Борна-Оппенгеймера. Так как в ионе молекулы водорода только один электрон, то гамильтониан молекулы Нмол имеет вид:
— —(1/2)А + У,„ (г) (4.2)
и уравнение Шредингера совпадает с уравнением Кона-Шема:
[—(1/2) А + Умо, (г)]Ф,(г) — Е1Ф 1(г), (4.3)
где
1 1
У (г) —
мол \ /
г — Я /2
г + Я /2
(4.4)
Для полной энергии имеем простое выражение:
Е — ТТ + Е1- (45)
Ят
В работе [3] предложено удобное приближение для волновой функции иона молекулы водорода:
3
4
^1(Г) — ( 6ХР ^ Га ) + 6ХР ^ Г )) .
(4.6)
Здесь использованы обозначения:
' ' ' Р — 4Ят,
га — г — Ят /2, гь — г + Ят /2
а т Р Ь т
(4.7)
5 — (1+ Р + Р2/3) ехр (—р).
Будем ее использовать в качестве опорной функции ^ (г ).
Для начала итерационного процесса необходимо задать хорошее приближение для величины Е1. Для этого воспользуемся диагональным матричным элементом гамильтониана по опорной функции:
ех —1^1 (г) н ^1 (г )аг ——4 /2+
+(4(4—1)—4<+4(4 — 2 )т)/(1+5)' (48)
где
С — 1 Р(1 — (1 + Р) ехр(—2р)),
(4 9)
г —(1 + р) ехр(—р).
Как и ранее во втором разделе данной статьи будем искать решение уравнения (4.3) с потенциалом (4.4) в виде:
Ф1(г) — ^(г) + у1(г,г)Д/г . (4.10)
Была составлена программа для нахождения функции ,У1(г, %) с учетом вида опорной функции (4.6) и действия на нее гамильтониана Нмол (4.2). Результаты расчетов приведены на рис. 1.
Расстояние между протонами измеряется в радиусах Бора. Видно, что минимум достигается при Ят = 1.9975 □ 0.1057нм, что хорошо согласуется с экспериментальным значением [4].
На рис. 1 по рассчитанным точкам построена парабола:
Е1 =-0.600137 + 0.050586(Кт -1.9975)2
Рис.1. График зависимости полной энергии иона молекулы водорода Е от расстояния между протонами я ■
Если считать, что потенциальная энергия колеблющихся протонов в ионе молекулы водорода описывается этой параболой, то нулевая энергия колебаний протонов
сок!2 —13.6058 О 05050765 2 — 0 143еВ . С учетом этой поправки энергия связи иона молекулы водорода получается равной -2.582 еВ. Экспериментальное значение энергии связи равно -2.65 еВ. Разница составляет около -0.068 еВ. Эта разница связана с тем, что при решении использовался шаг равный 0.03 радиуса Бора.
Заключение
В настоящей работе предложено обобщение приближения аппроксимированного локального потенциала на случай молекул. Для решения уравнения Кона-Шема разработан новый алгоритм - метод опорной функции. Этот метод опробован на тестовой задаче: вычисление энергии диссоциации иона молекулы водорода.
Для шага 0.03 радиуса Бора точность вычисления энергии диссоциации составила
0.07 еВ. Если учесть, что на самом деле вычисляется не энергия диссоциации, а полная энергия молекулы, которая составляет - 16.37 еВ, то относительная точность расчета оказывается лучше 0.4%. Возможно, что уменьшение шага и увеличение числа итераций позволило бы еще лучше согласовать теоретическое и экспериментальное значения энергии связи. Все же следует признать, что и для такого грубого шага метод опорной функции оказывается достаточно точным.
Литература
1. Шкловский А.Г. Аппроксимация обменно-корреляционного потенциала в методе функционала электронной плотности// Научные ведомости БелГУ. Серия Физико-математические науки, 2007. -№6(37). - Вып. 13. - С. 150-155.
2. Теория неоднородного электронного газа / Под ред. С. Лундквиста, Н. Марча, - М.: Мир, 1987. -400 с.
3. Флюгге З. Задачи по квантовой механике. Т. 1, - М.: Мир, 1974. - 341 c.
4. Madsen M.M., Peak I.M. Eigenparameters for the Lowest Twenty Electronics States of the Hydrogen Molecule Ion. - «Atomic Data», 1971, v.2(3), p.171.
SOLVING KOHN-SHAM EQUATION FOR MOLECULE WITH METHOD OF SUPPORT FUNCTION
A.G. Shklovskij1), M.A. Shklovskaya2)
1 Belgorod State University, Studencheskaja street, 14, Belgorod, 308007, Russia 2) Moscow State Technical University named after N.A. Bauman,
Vtoraja Baumanskaja street, 5, Moscow, 105005, e-mail: [email protected]
The generalization of approximation of exchange-correlation energy in the density functional theory was described in this article. This generalization was used for finding exchange-correlation potential of diatomic molecule. New algorithm of solving Kohn-Sham equation for molecule, method of support function, was given. Calculation of dissociation energy of hydrogen ion in molecule was applied for testing this method.
Key words: density functional theory, exchange-correlation potential, Kohn-Sham equation, method of support function.
