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

Вестник Челябинского государственного университета. 2013. № 9 (300). Физика. Вып. 16. С. 60-69.

МЕТОДИЧЕСКИЕ ЗАМЕТКИ

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

МЕТОД ОПОРНОЙ ФУНКЦИИ ДЛЯ РЕШЕНИЯ УРАВНЕНИЙ ШРЁДИНГЕРА И КОНА — ШЕМА

Рассматривается численное решение одномерных и двумерных уравнений Шрёдингера и Кона — Шема с нулевыми граничными условиями на бесконечности с использованием нового метода опорной функции. При решении уравнения Кона — Шема применён новый вариант модифицированного локального обменно-корреляционного потенциала. Дано сравнение с экспериментом.

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

1. Введение. В квантовой механике большое количество стационарных задач по нахождению связанных состояний решается с использованием теории возмущений Шрёдингера или Гейзенберга [1]. Типичной постановкой задачи в этом случае является наличие некоего гамильтониана Н0, для которого имеется полный набор собственных функций и собственных значений:

Н0¥ я (г ) = Еп Т п (г ). С1)

Необходимо получить некоторые собственные функции и собственные значения оператора Н = Н 0 + V (Г):

(Н0 + V (г))Ф п (г ) = Еп Фп (г ). (2)

Если достаточно рассмотреть только поправку первого порядка к энергии или к волновой функции, то применение теории возмущений не вызывает затруднений.

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

В настоящей работе рассматривается решение уравнений Шрёдингера и Кона — Шема методом опорной функции. Этот метод позволяет воспользоваться вместо теории возмущений набором простых алгоритмов и получить достаточно быстрое решение задачи (2) для наиболее низких энергетических уровней. Рассмотрим в начале случай, когда уравнение (2) является обычным нерелятивистским уравнением Шрёдингера.

2. Метод опорной функции для нерелятивистского уравнения Шрёдингера. Рассмот-

рим вариант метода опорной функции, применяемый для решения уравнения Шрёдингера в заданном потенциале. Для этого представим функцию Ф п (г ) в виде

Ф „(г ) = Т „(г ) + у„(Г ). (3)

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

уп (г) малой добавкой.

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

Тп(г) = I а, ■ У,(г), (4)

1<; <;т

где У, (г) — выбранный для решения задачи (2) набор базисных функций:

(Н о + Гок (г )У (г ) = екУк (г ). (5)

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

Перепишем уравнение Шрёдингера (2) для

функции Уп (г ) :

(Но + (V (?) - Еп ))уп (г) = Оп (г). (6)

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

°п(г) = Iа; -(Еп - ек - а(г))х

1<;<;т

хУ; (гХ Qk (г) = V (г) - К; (г). (7)

Для нахождения энергии Еп воспользуемся уравнением (6). Для этого умножим его слева на Фп (г ) и проинтегрируем по г . Учитывая, что

уравнение (2) является самосопряжённым, то, сокращая подобные члены, получим

I»' [(£,-г1УЬ--Q■J-] = 0. (8)

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

Ьп} = |Фп(г)Уу {г)ёг, Q“ =

= /Фп (г О/ (г У (г ^. (9)

Система уравнений (8) является однородной и имеет ненулевое решение, если её определитель равен нулю. Очевидно, что эту систему следует решать самосогласованным способом. В нулевом приближении Ф (г) совпадает с Ук (г ), таким, что энергия Еп наиболее близка к г,. При этом формула (9) даёт интеграл перекрытия из приближения сильной связи. Если в формулу (9) подставить Qj (г) нулевого приближения, то получится система уравнений сильной связи в первом приближении. Приравнивая определитель системы (8) с коэффициентами Ь^ и Qrп нулевого приближения к нулю, можно получить энергии первого приближения Еп и коэффициенты агп нулевого приближения. Используя эти коэффициенты в формуле (4), можно начать процесс итерации, полагая в нулевом приближении уа = 0. Формула (8) позволяет начать итерационный процесс. Какое-то количество итераций по этой формуле можно сделать не вычисляя заново коэффициенты а"п . На следующем этапе

необходимо уточнить этот набор коэффициентов.

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

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

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

(Н0 + (V(г) - Еп))gn(г) = 0. (10)

Здесь gnnп) = п,ипг)-ё2,„пп) . Так как для

уравнения (10) Еп является невырожденным собственным значением, то gn (г) либо равна нулю,

либо пропорциональна Ф п (г):

gn (г ) = ^ ■ Фп (гX gl,п (г ) = ^ ■ Фп (г ) + g2,п (г ). (11)

Так как нас интересует нормированное решение уравнения (2), то подставим (11) в (3):

|Фп (г^ = |Тп (г) + ^ ■ Фп (г) + g2,п (г)| =

= |(1 + Ф п (г )|. (12)

Подставим (12) в условие нормировки для функции Ф п (г ) и получим уравнение

(1 + Я)2 = 1, Х1 = 0, X 2 = -2. (13)

Таким образом, мы видим, что в отличие от уравнения (2) уравнение (6) имеет два различных решения. Первое решение g1 действительно является малой добавкой к опорной функции. Второе решение g2 , соответствующее Х2 = -2,

приводит согласно (3) к функции -Фп (г). Эта

функция в уравнении (2) не давала нового решения, так как решение линейного однородного уравнения находится с точностью до произвольного коэффициента. Уравнение (6) не однородно, и очевидно, что функции g1 и g2 — это разные решения. Поэтому при построении алгоритма следует позаботиться о том, чтобы решения уп (г) по сравнению с опорными функциями

были малы. В этом случае второе решение окажется исключённым.

Перейдём к построению итерационной процедуры решения уравнения (6). Будем считать, что уравнение (6) удовлетворяет нулевым граничным условиям: yn(rc) = 0.

В качестве примера такого уравнения Шрёдингера рассмотрим задачу об ионе молекулы водорода. Это тестовая задача, имеющая точное решение. Энергия диссоциации иона молекулы водорода известна, и есть независимые методы решения уравнения Шрёдингера для электрона в поле двух протонов, находящихся на заданном расстоянии Rm, поэтому в работе [2] описанная выше теория была применена к вычислению этой энергии. Для энергии диссоциации Ed использовалась формула

(14)

Ed = E

d св

-®k/2.

Здесь ю, — частота колебаний протонов в ионе молекулы водорода, которая вычисляется в приближении Борна — Оппенгеймера. Так как в ионе молекулы водорода только один электрон, то гамильтониан молекулы Нмол имеет вид

Нюл = -(1 / 2) А + (г )

(15)

и уравнение Шрёдингера для основного состояния имеет вид

[-{1/2 )А + ^ (r )] x x^{f) = Ei^(f X

где

Кол (r ) = -1/

r -R /2

-1/

r + R /2

(1б)

. (і?)

Для полной энергии имеем простое выражение

E = УR + Єі.

(IS)

В работе [2] в качестве опорной функции Т1(г) использовалось удобное приближение

для волновой функции иона молекулы водорода, описанное, например, в [3]:

Ti(f) = \ 7^2п (1 + S) x x(exP (-^ r»i ) + exp (-^- Га2 )). {19)

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

г»1 = \Г - К /2, га 2 = \Г + К /2,

Р = ^ , 5 =(1 + Р + Р 73)еХР (-Р). (20)

Для начала итерационного процесса необходимо задать хорошее приближение для величины Е1. Для этого воспользуемся диагональным матричным элементом гамильтониана по опорной функции:

81 = Ч2/2 +

+ ( (;- 1)-^ + £ (£-2) т)(1 + 5), (21)

где

Z = V p (V - (1 + p )exP{-2p) ), т = ( + p )exp{-p).

(22)

Как и ранее в (3), будем искать решение уравнения (16) с потенциалом (17) в виде

Фі(г ) = ^(г ) + у (г, 2)! 4т . (23)

Перепишем уравнение (16) в виде

д2У(т, 2) + д2y(т, 2) + дг2 д22

+ (0,25/г2 - 2[ГМоЛ (Г) - ^])х

X у(г, 2) = -Ц(г, 2). (24)

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

GV{r,z) = 2-47- 42Д/2п(1 + S) x x I {Ei -Si - Qk{Г)) - exP ("4- r»k ) (25)

Qk {r ) = Кол {Г ) + V rak .

(2б)

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

п, = 1, 2,...,N1, п2 = 1, 2,...,

N2, хг(П1) = гп1, 2С (п2) = 2п2,

хг (п1 +1) = хг (п1) + Нг (п1),

2С( п2 + 1) = 2С( п2) + Н2( п2). (27)

Здесь массивы узлов гп^ с различными номерами связаны друг с другом соотношением (27). При этом хг(1) = 0 и хг (N1 ) = гС. Аналогично для

массивов 2 имеем 2с(1) = -2С, 2c(N2) = 2С. Так как мы собираемся находить функцию У„,я (г>2)

приближённо с помощью итерационной процедуры и следить за тем, чтобы она была мала по сравнению с Ф1 (г ) , то в качестве нулевого приближения следует положить упт (г, г) = 0. Введём два массива ат(п1,п2) и Ьт(п1,п2).

Первый массив будет задавать сеточную функцию, приближённо описывающую функцию уп т (г, г) на нечётных итерациях, а второй — на

чётных. Тогда вместо уравнения (24) получим разностное уравнение:

Ьт(п1, п2) = (Аг (п1, п2)ат(п1 +1, п2) +

+Вг (п1, п2)ат(п1 -1, п2) +

+Аг(п1, п2 )ат(п1, п2 +1) +

+Вг(п1, п2 )ат(п1, п2 -1) +

+Gm (п1, п2 ))^т (п1, п2).

Здесь введены обозначения: 2Н2(п2)

Аг(nl, п2) = ВГ(п^ п2) =

(Нг (п1) + Нг (п1 -1)) Нг2 (п1, п2)

(Нг (п1) + Нг (п1 - 1))Нг (п1 -1) Нг2(п1, п2) = 2Н2(п2)Нг(п1).

2Нг (п1)

А2(п1, п2) = В2(п1, п2) =

(Н2(п2) + Н2(п2 - 1)) Нг2(п1, п2)

(Н2(п2) + Н2(п2 - 1))Н2(п2 - 1)

Бш(п1, п2) = У [ Вг (п1, п2) + В (п1, п2) + 0,125

+(К (п1, пг)-Ех —

г) Нг2(щ, п2)],

хг (п1 +1)2' бмЦ, п2) = бЮ;, 2^™ (п1, п^).

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

2Н2(п2)

(31)

Вг п п2) =

Щ^, п2) =

Нг (п1 -1) ’ 2Нг (п1) Н2 (п2 - 1)

(32)

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

(28) Ф1(г, г), то такая точность приемлема.

Результаты расчётов приведены на рис. 1.

Расстояние между протонами измеряется в радиусах Бора. Видно, что минимум достигается при Ят = 1,9975 = 0,1057 нм, что хорошо согласуется с экспериментальным значением.

На рис. 1 по рассчитанным точкам построена

(29) парабола:

Е1 = -0,600137 + 0,050586( -1,9975)2. (33)

Если считать, что потенциальная энергия ко-

(30) леблющихся протонов в ионе молекулы водорода описывается этой параболой, то нулевая энергия колебаний протонов

Рис. 1. График зависимости полной энергии иона молекулы водорода Е от расстояния между протонами Я.т

&к/2 = 13,6058,

0,050586 • 2

918,075

= 0,143 еВ.

С учётом этой поправки энергия связи иона молекулы водорода получается равной -2,582 еВ. Экспериментальное значение энергии связи равно -2,65 еВ. Разница составляет около -0,068 еВ.

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

Ехс [ п ] = -

каг

N.

3в(Ые -1) 4 N.

| Сгх п ( г)

(34)

где к = (Ь+Ы/ъ)(Ые -1)/Ые, Ше — число электронов в атоме; ^(х) — аппроксимирующая функция, введённая в работе [4]:

Е (х) = х • Ра( х)[а3 + а4 Ра( х) +

+(1 - а3)Ра 2( х)] • ехр(-а2 х), (35)

где

Ра( х) = (х • а1 -1)/( х + 0,5). (36)

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

*с

| ^ (г )п(г )г 2йг = 0, (37)

0

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

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

атомов или ионов совпали с экспериментальными данными. При этом остальные вычисленные энергии внутренних оболочек должны быть как можно ближе к экспериментальным.

В рамках нерелятивистской теории функционала электронной плотности [5] в работе [6] было рассмотрено выражение для вычисления полной нерелятивистской энергии атома:

Е = Те[п] + Екаг[п] + ЕХс[п] +1 Уе(г) • п(г)аг. (38)

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

после усреднения по углам имеет вид

<г) = XIК і(г )|2 • ч„, і .

(39)

Здесь чп,і — числа заполнения электронных оболочек атома или иона с учётом суммирования по магнитным квантовым числам т и спину. Для функционала Те[п], входящего в (38), имеем выражение

Ые

те [п] = X Еп - Г V(г)п(Г)СГ.

п=1

(40)

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

Ые

Е = X Еп -| ¥(г)п(г)СГ +

п=1

+Екаг [п] + Ехс [п] + / ¥е (г)п(г)СГ. (41)

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

Е,„, [ п] = 2 Г

1 Гп(г1)п(г2)

г1 - г2

сгхсг2.

(42)

Согласно теореме Хоэнберга и Кона [7] функционал полной энергии (41) для электронной плотности (39) имеет минимум. Так как минимум функционала (41) определяется из равенства нулю его первой вариации, то при очень малом изменении электронной плотности п(г), например, за счёт малого изменения потенциала У(г), поправки к полной энергии Е будут второго порядка малости. Поэтому вблизи от правиль-

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

V (г) = V; (г) + Гкаг (г) + Гхс (г).

(43)

Здесь Ve(г) — потенциал притяжения электронов к ядру, выражение для потенциала Vhaг(r) с учётом (39) имеет обычный вид:

Каг (Г) = |п(г2 )/\Г - ^1 ^2 =

Ме/г •(1 - Чмг(г)), (44)

Чъаг (г) = 4П/ Ме - г) Х

К

х | х2п(х) (1 - г/х^ёх. (45)

г

Кс (г) = - [(1 - Чъаг (г) - к • Чхс (г))/г +

+ р/К (К-1) п3(г) , (46)

Чхс (г) = 4П Р(г)• ©(Кс - г) Х

К

х | х2п(х) (1 - г/х) Р(х)ёх. (47)

г

Получившийся потенциал V(г) используется для решения уравнения Кона — Шема в атоме [8] и отличается от потенциала, использованного в работе [6], только выражением (35) для аппроксимирующей функции Р(г). При этом, как и в релятивистском случае, для большинства многоэлектронных атомов удаётся получить лучшие, чем в приближении Хартри — Фока значения энергий внутренних оболочек атомов.

В случае малого числа электронов и особенно для систем с двумя электронами результат не столь хорош. Рассмотрим этот случай более подробно. Изменяя параметры в обменно-корре-

ляционном потенциале для многоэлектронных атомов, мы изменяем значения энергий внутренних оболочек не зависимо от значения энергии ионизации. Это позволяет точно подогнать две величины: полную энергию и энергию ионизации. Для двух электронных систем это невозможно, так как имеется только одна 5-оболочка. Поэтому мы можем подогнать под эксперимент энергию 15-состояния либо полную энергию. В таблице приведены полные энергии и энергии 15-состояния двухэлектронных систем: атома гелия, отрицательного иона водорода и положительного иона лития.

В таблице использованы обозначения: Еехр— экспериментальная полная энергия (эВ) [11], Е — полная энергия атома по формуле (41), Е^ — полная энергия в приближении Хартри Фока [12], Е1ехр — экспериментальная энергия ионизации, Е51 — энергия 15-уровня для уравнения Кона — Шема, Е51^ — энергия 15-уровня для уравнения Хартри — Фока.

При решении уравнения Кона — Шема в двухэлектронных атомах и ионах были использованы одинаковые значения параметров а3 = 0,220 и а4 = 0,258. При заданном в таблице параметре а2 параметр а1 определялся из формулы (37).

Интересной системой с двумя электронами является молекула водорода. В работе [9] было получено решение уравнения Кона — Шема для этой молекулы с использованием обменно-корреляционного потенциала, предложенного там же. При этом была достигнута точность решения около 0,1 эВ. В настоящей работе решим уравнение Кона — Шема для молекулы водорода с обсуждаемой выше модификацией обменно-корреляционного потенциала.

3. Вычисление молекулярного потенциала. Теперь можно вернуться к вопросу о нахождении потенциала У(г, 2) для двухатомной молекулы. Чтобы не усложнять изложение переходом к релятивистскому уравнению для молекулы, будем использовать нерелятивистские уравнения (15) и (16), но для потенциала молекулы V™ (г) воспользуемся выражением

Полные энергии и энергии 1*-состояния двухэлектронных систем

Ион Еехр эВ Е, эВ ЕЪр эВ ^ер эВ Е51, эВ £51м , эВ а2 в

Н+ -14,3599 -14,3562 -14,0262 -0,7542 -1,6378 0,33 0,694 0,212

Не -79,006 -79,001 -77,9 -24,587 -25,6198 -24,98 0,694 0,0749

ы- -198,091 -198,083 -196,791 -75,641 -76,736 - 0,94 0,0436

Кол (Г ) = К (Г ) + Каг (г ) + Кс (г ). (48)

Здесь Ге (г ) — потенциал притяжения электронов к покоящимся ядрам в молекуле, выражение для потенциала Vhaг (г ) имеет обычный

вид:

Каг (г ) = |П (А)/\г - А\ Сг2 ,

(49)

а обменно-корреляционный потенциал определяется выражением

VXc (1) = 5ЕХс [ п]/5п(1). (50)

При приближении к реальному межъядерно-му расстоянию в молекуле электронная плотность значительно отличается от суммы атомных плотностей. Чтобы это учесть, запишем её в виде

п(г) = 2£|Фп (г)|2 =

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

При переходе от атома к молекуле приходится несколько изменить выражение (3) для полной энергии:

Е = Те [п] + ЕЫг [п] + Ехс [п] +

+| Ve (г) •п(г)Сг-

7 7

^1 2

(54)

Здесь а — расстояние между ядрами атомов

в молекуле. Постоянная кулоновская энергия отталкивания ядер друг от друга вводится для того, чтобы молекула была электронейтральна.

Напомним, что согласно теории функционала электронной плотности [5] кинетическая энергия молекулы задаётся формулой

Ме

Те [п] = 2 • X Еп -1 Умол (г )п(г )Сг. (55)

п=1

па1(

1 а

г----

2

) + па 2 (

) + Пс (г). (51)

Здесь Фп (г) является решением уравнения

(16) с энергией Еп, а Пс (г) показывает отличие реальной плотности молекулы от суммы атомных электронных плотностей и находится самосогласованным образом. Приближение для обменно-корреляционной энергии Ехс [п] тоже изменяется по сравнению с (34). Дело в том, что плотность п(г ) в молекуле будет обладать уже

не сферической, а цилиндрической симметрией. Эта плотность, как видно из (51), имеет два центра сгущения. Поэтому вместо (34) имеем

4

и Л/ _ I »

Ехс [п] = -

Еъаг , 3(Ме - 1)

| й/1Р(/1)п(/1)3 +

А |Р1(г1)п(г1)р1 (г2)п(г2) _

2 12 I 1 - г2 I

(г )п(г ) р2(12)п(12)

г - г

1 '2

(52)

Подставим (51) в (49):

Vhar (г ) = ^ Ъаг (

1 а

г----

2

) +

+^ Ъаг (

) +

(53)

Здесь Еп — энергии молекулярного спектра, полученные в результате решения уравнения Кона — Шема (2), а VмOIl (г) — молекулярный

потенциал, входящий в уравнение (2). Энергия Хартри, как обычно, даётся выражением вида (42), которое отличается только заменой электронных плотностей в атоме на электронные плотности в молекуле. Функционал полной энергии молекулы (54) с учётом выражений (48), (49) и (55) с точностью до аппроксимирующих констант известен.

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

Vxc (г) = -[^г (г)/Мет + (Мет - 1)/

Метв(г )п(г )1 + кр (г - 2) Х

х\п(г1\

г - г

•Р1(г -

2сг1

г - г2

+£2 Р2(г + 2) |

ач гп(г1).

г - г

1Р2(г1

21

. (56)

Здесь Ыет — число электронов в двухатомной молекуле.

Если потенциал ¥мол (г) известен, то уравнение Кона — Шема не отличается от уравнения Шрёдингера (6). Поэтому для его решения можно применить те же формулы (3)-(9). Вернёмся от общей теории двухатомных молекул к молекуле водорода. Учитывая цилиндрическую симметрию этой молекулы, перепишем уравнение (6) для Сигма электронов в виде

д2у(г, х) д2у(г, х) + дг2 д22

(0,25/г2 - 2[КМоЛ (г, х) - £!])у(г, х) = -в^г, 2). (57)

Так как потенциал Умол (г, х) не зависит от угла ф, то и волновая функция основного состояния тоже не зависит от угла ф.

При нахождении Фх(г, х) и вычислении в первом приближении поправки ц^г, 2) можно одновременно найти поправку Укс(г, 2) к потенциалу Хартри (53). Заметим, что этот интеграл есть решение уравнения Пуассона:

Л^ы(г) = -4пП№ (г). (58)

Так как цс1 (г ) не зависит от угла ф, то будем

численно решать уравнение Пуассона в цилиндрической системе координат. Перепишем уравнение (58) в виде

д2и д2и и , Л

- + —- +—- = -4 (г, 2). (59)

Здесь

и (г, х) = ¥м (г, х) • 4г.

(60)

В приложении показывается, как находятся граничные условия на линиях г = гт и 2 = 2т. На линии г = 0 имеем условие и(0, 2) = 0. На линии 2 = 0 в силу симметрии имеем

ди (г, 2),

дх

= 0 . Таким образом, решение урав-

нения (59) с указанными граничными условиями даёт нам возможность вычислить величину и(г, 2), а следовательно, и поправку Vhd(г, 2) к потенциалу Хартри во всех внутренних точках молекулы. Алгоритм решения уравнения (59) аналогичен алгоритму (27)-(32), который был применён при решении уравнения (24).

Для получения энергии связи молекулы нужно осуществить вычисление полной энергии молекулы водорода для различных расстояний между ядрами по формуле

Е = 2Е! - Еш [п] + Ехс [п] -

-{ Кс <7 )•п(г № + V К.

(61)

Для вычисления полной энергии по формуле (61) в обменно-корреляционном потенциале (56) использовались значения параметров:

Р = -0,2226, к1 =

= к 2 = (3 + 213 )Д, а 3 = 0,9, а4 = 0,042, а2 = 0,28728 +

+0,96113(Ят -1,397);

х=0

Рис. 2. График зависимости полной энергии молекулы водорода Е от расстояния между протонами Ят и график параболы Е = -1,17031 + 0,11013(Кт -1,397)2

Отметим, что при изменении расстояния Ят мы меняли только параметр а2 и, естественно, пересчитывали при этом параметр ах. Точность расчёта энергии при решении уравнения Кона — Шема составляла 0,01 эВ. Интересно, что оказалась достаточно линейной интерполяция параметра а2 вблизи минимума полной энергии. Вблизи от этого минимума график можно аппроксимировать параболой

Е = -1,17031 + 0,П013(Лт — 1,397)2.

Хотя мы и описывали ядра атомов водорода в молекуле в адиабатическом приближении как покоящиеся классические частицы, на самом деле они участвуют в так называемых нулевых колебаниях. В этом случае должна решаться квантовая задача о колебании двух частиц вокруг положения равновесия в потенциале Е(Ят). В качестве эффективной массы выступает Ме^= 0,5М, где Мс — масса атома водорода. Энергия нулевых колебаний в атомной системе единиц ю0 будет иметь вид

M

ef

0,22026

27,21139б(еВ) =

'918,075

= 0,4215(еВ).

(б2)

Для получения энергии связи нужно из полной энергии молекулы вычесть полную энергию двух изолированных атомов водорода, поэтому энергия связи ЕСВ определяется формулой

Есв = Е (Яттт) + ^ - 2 Ен =

= -4,4237(еВ). (63)

Экспериментальная энергия связи Еэксп = -4,45 (эВ). Видно, что полученный результат совпадает примерно с точностью, использованной при решении уравнения Кона — Шема.

Список литературы

1. Ландау, Л. Д. Квантовая механика. Нерелятивистская теория / Л. Д. Ландау, Е. М. Лифшиц. М., 1963. 702 с.

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

3. Флюгге, З. Задачи по квантовой механике. Т. 1. М. : Мир, 1974. 341 с.

4. Береговой, А. В. Локальное приближение функционала обменно-корреляционной энергии для релятивистских атомов / А. В. Береговой, А. А. Плесканев, А. Г. Шкловский // Науч. ведомости БелГУ. Сер. Математика. Физика. 2012. № 23 (142). Вып. 29.

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

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

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

8. Kohn, W. Self-Consistent Equations Including Exchange and Correlation Effects / W. Kohn, L. J. Sham // Phys. Rev. A140. 1965. Р 1133-1138.

9. Шкловский, А. Г. Решение уравнения Кона — Шема для молекулы водорода методом опорной функции // Науч. ведомости БелГУ. Сер. Математика. Физика. 2011. № 11 (106). Вып. 23. С. 52-63.

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

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

12. Фудзинага, С. Метод молекулярных орбиталей. М. : Мир, 1983. 461с.

Приложение

Вычисление граничных условий для потенциала Ука

Для нахождения граничных условий в точках рт и zm вычислим Уы в произвольной точке г за

1

пределами молекулы водорода. В этой области % (г ) = 0. В формуле (53) разложим полиномам Лежандра для |Я,| < |г|:

в ряд по

1 w |Г - Ч = ' I ГІ r п 1 2 1

r П"0 r VI 1J

Рп (cos у). (П1)

Здесь y — угол между r2 и r . Так как для основной части плотности nd (r2) имеем \r21 ^ |r , ограничимся в этом разложении тремя первыми членами. Учтём также, что

rr xx + yy. + zzn

cos у = = 2 /Г2-------2. (П2)

\r\%\ \r\\r2\

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

1 1 kl |r2|2 3cos2 y-1

-cos Y + J—3--------------. (П3)

Подставим (П3) в поправку к потенциалу Хартри (53):

¥М (Г) = 1 І Па (Г2 )аГ2 + 7“Г ІПа (Г2 )(ХХ2 + УУ2 + 222 )аГ2 + 11 ^2 |2(3 С0^ У - 1)аГ2 • (П4)

Г |г| 2 |г|

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

І Па (Г2 )(хх2 + УУ2 + 222 Щ =

= х І Па ( Г2 ) х2 аГ2 + УІ Па ( Г2 ) У2 аГ2 + 2 І Па ( Г2 ) 22 аГ2 = 21 Па (Г2 ) 22 аГ2 =

рт 2т

= 2п2 І Р2аг2 І Па (Р2,22) ■ 22а22 = 0. (П5)

0 -2т

Этот интеграл равен нулю в силу симметрии па(р2, 22) = Па(Р2, -22), поскольку атомы углерода, образующие молекулу, одинаковы. Таким образом получим

¥ка(Г) = - лз 1Г|2(3С0э2У -1■Па(г2)а2• (6)

21 г\

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

¥ы(г,2) = (П<)■ (Г2 -222)/(г2 + 22)2 ■ Ім, (П7)

где

2 ~ 2

“m m ■

I _ гг Пё (р2 , 22 )(р2 2 22)р2ё р2 ё22 (8)

“ .1 г ,2,2ч . ' '

-?ш 0 (р2 + 22 )

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

и(гт, 2) _ (у2)(гт- ъ'1)4гтт1 от+^2)2 • ^. (П9)

Аналогично вычисляется и величина и(гт, 2):

и ^, 2т ) _ (П? )(г 2 - 2 (Г 2 + 4)2 • 1Иё . (10)