Система матрично-векторных уравнений в многоконфигурационном методе Хартри-Фока Текст научной статьи по специальности «Математика»

УДК 539.183

СИСТЕМА МАТРИЧНО-ВЕКТОРНЫХ УРАВНЕНИЙ В МНОГОКОНФИГУРАЦИОННОМ МЕТОДЕ

ХАРТРИ-ФОКА

М. С. Лицарсв, О. В. Иванов

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

Ключевые слова: многоконфигурационный метод Хартри Фока матрично-векторного уравнения.

Многоконфигурационный метод Хартри Фока (МКХФ-метод) применяется во многих областях физики конденсированного состояния вещества, квантовой химии, атомной спектроскопии, как правило в тех случаях, когда необходимо достичь высокой точности расчетов электронной структуры атомов или ионов.

Являясь вариационным. МКХФ-метод требует решения системы интегро-дифференциальньтх уравнений. Применение конечно-разностных схем [1. 2] не может гарантировать в общем случае сходимости решения на отдельном тттаге итерационного МКХФ-процесса. а получаемые таким способом волновые функции не обладают заданной степенью гладкости, что необходимо для ряда твердотельных приложений (например, при построении псевдопотенциалов [4. 5]. Существующие методы решения уравнений Хартри Фока, основанные на разложении по наборам базисных функций [6] (в которых, как правило, используются слэтеровские или гауссовы орбитали) не всегда позволяют решить уравнения Хартри Фока Рутана [7] с З&ДШШОИ ТОЧНОСТЬЮ.

Получение МКХФ-уравнений при расширении многоэлектронного базиса и5соответ-ственно. при увеличении числа одноэлектронных состояний, которые необходимо определять. является отдельной, и при том весьма трудоемкой задачей. В случае метода Хартри Фока, вывод уравнений осуществляется аналитически [1 3]. МКХФ-уравнения

строятся на основе уравнений Хартри Фока по усложненным правилам [8], которые с алгоритмической точки зрения крайне сложно формализовать и обобщить для случая произвольного многоэлектронного базиса.

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

В данной работе для определения одноэлектронньтх состояний, входящих в много-электронньтй МКХФ-базис. получена система матрично-векторньтх уравнений и показано. что она эквивалентна уравнениям МКХФ-метода. Для этой системы разработаны алгоритмы построения матрично-векторньтх уравнений, найден итерационный быстро сходящийся метод решения для случая произвольного многоэлектронного базиса. Приводятся примеры расчетов для различных атомов.

Постановка задачи. Рассмотрим нерелятивистское уравнение Шредингера для атома с зарядом ^ и числом электронов N которое в атомных единицах имеет вид

1 N N 1 N 1

Z ^ ^ r + ^ ^ r. .

I г

2

i=1

i=1

г<j

ij

Ф (q) = EФ (q) .

(1)

Здесь гг - радиус-вектор г-ro электрона, гг = |гг|, Г] = |гг — г]|, q = {q1,q2, ■ ■ ■ ,qN}■ Через qi обозначена совокупность радиус-вектора гг и спиновой переменной аг, которую, несмотря на то, что уравнение (1) явным образом не зависит от спина, приходится вводить вследствие принципа Паули [9].

Обычно в МКХФ-методе решение уравнения (1) ищется в виде конечного разложения по многоэлектронному CSF-базису (configuration state functions) [1]

jmax

*(q) = £ CjФ] (q). (2)

j=1

Каждый его элемент Ф] (q) соответствует определенной электронной конфигурации

[(Ulll)W 1 Ы2Г ■ ■ ■ (nv lv )WV ]j

с числом оболочек v и числом электронов wi на г-ой оболочке, Y^i=1 ставляет собой линейную комбинацию слетеровских детерминантов

ф] (q) = Y. ldeta!j ...aN ). k=1

(3)

N, и пред-

(5)

±1, если

(6)

|а) = \nlmims). (7)

Здесь п - главное квантовое число, I - орбитальный угловой момент (азимутальное квантовое число), т1 - магнитное квантовое число, т3 - проекция спина на выделенную ось г. В выражении (6) переменные (г,9,ф) обозначают сферические координаты, У1т1 (9, у>) - сферическая функция [10], \та (&) ~ спиновая часть одноэлектронного соСТОЯНИЯ.

В данной работе предполагаются известными [11] число ]тах, входящее в формулу (2), коэффициенты А— и наборы квантовых чисел а1 ,..., а^- одноэлектронных состояний, формирующих слэтеровские детерминанты в выражении (4); У] € [1, ]тах], У к € [1,^ ]• Кроме этого, предполагается, что неизвестные коэффициенты С) разложения (2) определяются в рамках самосогласованной МКХФ-процедурьт стандартным способом [1].

Таким образом, необходимо определить все неизвестные радиальные функции Рп1(т) с различными п1, входящие в базис (4). При этом, так как волновые функции (6) ор-тонормированы, на радиальные функции Рп1(т) налагаются дополнительные условия ортонормировки

/ Рп^РпП^йт = 5пп', I = 0, 1 ...1тах- (8)

./о

Одночастичный базис. Будем искать одночастичные радиальные функции Рп1(т) в виде подкласса функций, представимьтх в виде конечного разложения

К1тах

Рщ(т)=^ъп1 Як1 (т) (9)

к=0

Детерминант Слэтера здесь и далее обозначен через

1 "

т 1=1

суммирование ведется по всем возможным перестановкам т индексов, ет = перестановка т четная или нечетная соответственно. Одноэлектронное состояние

Фа(Я) = 1 Рп1(т)У1т{ (9,^)Хше И

в обозначениях Дирака записывается в виде

по ортонормированному базису

Q»(r) = Jï+bw- г'+'е-2 L'r (г) ■ (10)

Здесь La(x) - многочлены Чебышева-Лагерра [10, 12], l G [0,lmax]. Базис (9) применяется при решении различных задач математической физики [12].

Физический СМЫСЛ раДИаЛЬНЫХ функций Pni (r) состоит в том, что они соответствуют распределению электронной плотности в атоме, которая экспоненциально спадает при увеличении радиуса r [9]. Поэтому при проведении вычислений там, где функции Pnl(r) необходимо использовать в явном виде, область r G [0, ж) заменяется областью r G

[0, rmax] •

Вариационные соотношения. В выражение для полной энергии E = {^\И\Ф) уравнения (1) в соответствие с общими правилами вычисления матричных элементов многоэлектронного нерелятивистского гамильтониана между слэтеровскими детерминантами [3, 6] входят одноэлектронньте

U = [ pnai(r)Dipnbi(r)dr, Di = --d-2 + - r> (H)

где la = lb = l за счет произведения сферических гармоник, a = {nl}, и двухэлектронные

ИНТ6ГрШ1Ы

Rkbcd = ff Pa(r1)Pb(r2)Pc(r1)Pd(r2) dr2dn. (12)

Последние записываются как

рж i рж i

Rkbcd = Pa (r1)Pc(n) - Ybkd(r1)dr1 = Pb(r2)Pd(r2) - Ykc(r2)dr2, (13) Jo rl Jo r2

с помощью обозначения

Ykb(ri) = ri I Pnaia (r2)Pnbib (r2)dr2, (14)

k

./ т>

которое после элементарных упрощений принимает вид

У^П) = £ ^ Рпага(т2)РПъгъ(г2)¿г2 + Ц + ^ЫРщъ(15)

Подставляя разложение (9) в подынтегральные выражения (11) и (12) и производя варьирование по коэффициентам 1°, получаем

ЯТ Г-Т hni\ 1 Km ax

U±ab \_{hk }\ = во ti (Xni0 hnbi + £nÎ0 hnai )

nio10 = °i 'kko (0na hk + 0nb hk ),

,„ „0 i / j SkkoVUna uk ' Unb uk

),

dhk0 k=i

-дк Г^ШЦ Ктах

д-паЬе4 1_{ък }_| = Н0 то \ л гг\а1с,к,Ъйъпс1с + Но ¿по \ л п1ъ1а,к,асъпл1л + /17\

Й,по1о = °1а°па ''кокг ък1 + °1ъ °пь Чкокг ъкг +

дЪпо1о 1а па

дъко к1 = 1 к1 = 1

К

1с 1а,к,Ь4ъпа1а "пс / у Чкок1 ък1 • ~1Л~ пл / ; Чкок кг = 1 к! = 1

+81£япоТ, <1[кМъпк:1а + ¿1о8% £ ъй^ъпь4.

Здесь

ск1к2 = Яы (т)Ьг Як21 (т)йт = Ск2к1, (18)

•Уо

П—- = [ Як1н (т)Як212 (т)1 У-ь(т^т, (19)

при ЭТОМ

1г12,к,аЬ 1211,к,аЬ ^ 1\12,к,аЬ / 1г12,к,аЬ /пп\

ПккУ = Пк2к\ , но, вообще говоря, Пкг к'2 = Пк2к\ . (2°)

кгк2 ~ Чк2к1 ' пи> Чк1к2 ~Г Чк2к1

Входящее в функцию Лагранжа условие ортонормировки (8)

I N1 К1

1тах шах шах

С №1}} = ££ ^ (£ ъ% ъп11 - 8г]). (21)

1=0 г=1 кг = 1

1=

при варьировании по ъкгоо 1о дает следующий вклад в уравнение на коэффициенты Ь.

=£ (1 + г* К:1о. (22)

дъко° 3=1

Атом гелия в состоянии 1Б. Рассмотрим, в качестве примера, систему матрично-векторньтх уравнений для атома гелия в основном состоянии в базисе, состоящем из четырех конфигураций {1з2,2в2, 2р2}. В этом случае существует четыре оазисных состояния (4)

Ф1(Я1,Я2) = № 1001,1001), (23)

ЫЯ1,Я2) = -—=\Ы 1001, 2001) + —\Ы 1001, 2001), (24)

22

ЫЯ1,Я2) = № 2001, 2001), (25)

Ф4(Я1,Я2) = —^21-11, 2111) + — № 21 -11, 2111) + — \dett 210|, 2101). (26) у3 у3 у3

оо

Полная энергия системы имеет вид

Е = (2С1 + С|)/ю,ю + (2Сз + С|)12о,2о - 2л/2(С\С2 + С2Сз)110,20 + 2С2121,21+

+ ^10,10,10,10 + С3Я00,20,20,20 - 2^/2С1С2ЯЮ,10,10,20 - 2л/2С2С3ЯЮ,20,20,20 +

20 2 0 20 2 22 + С2 ЯЮ,20,10,20 + (С2 + 2 С1С3)ЯЮ, 10,20,20 + С4 Яц, ццц + 5 С4 Я21,21,21,21 + (27)

2 1 2\[2 1 2 1

+ С1С4 Я10,10,21,21 С2С4Я10,20,21,21 + С3С4 Я20,20,21,21 •

Ей соответствует система матрично-векторньтх уравнений

АПЪ10 + АцЪ20 + А^Ъ21 = -2А01Ь10 - Л?2Ъ20, (28)

АцЪ10 + АцЪ20 + А2зЪ21 = -А^Ъ10 - 2А12Ъ20, (29)

Аз1Ъ10 + А32 Ъ20 + АззЪ21 = -2А11Ъ21, (30)

где матрицы, стоящие перед векторами Ъга1, равны

(Ап)гз = 2(2С2 + С2)^ + 4Cln00,O,1O1O - ^ССгЦ0,0,1020 + 2^^и0,0,2020} (31)

(Аи)гз = -2\/2(С1С2 + С2Сз)& - 2V2ClC2V00,O,1O1O- (32)

112= 2(С1С2 + С2С3)Я,г] - 2 V 2С1С2П^

-2^2С2СзП00,0,2020 + 2(С2 + 2С1Сз)п00,0,1020

( Л \ 4 Г< Г< 01,1,1021 2\12 01,1,2021 /00ч

1 = С1С4Па С2С4Пц , (33)

( Л г,/пп2 I , ЛГ12 00,0,2020 л ГКг< п 00,0,1020 , о^<2 00,0,1010 Юл\

(А22)Ц = 2(2С3 + С2 Кгз + 4С3 Пц - 4У2С2СзП^ + 2С2 Пц , (34)

( Л \ 2^2 ^ 01,1,1021 . 4 П П 01,1,2021 /осЛ

(А23)Ц = С2С4'Цц + ^ С3 С4 , (35)

(А \ лп2<-1 , лп2 11,0,212^1 ^^»2 11,2,2121 /ог\

(А33)Ц = 4С4 Сгз + 4С4 Пц + 5 С4 Пц • (Щ

При этом А12 = А21, (А13)т = А31 и (А23)т = А32. Кроме этого необходимо понимать, что размерность Ктах векторов Ъ10 и Ъ20, вообще говоря, отличается от размерности К^ах вектор а Ъ21, то есть матр ицы А13 ж А23 - прямоугольные.

Решение системы, матрично-векторных уравнений. Рассмотрим систему матрично-векторньтх уравнений, записанную в общем виде

kmax (i) j 1 k=kmin(i)

Aij Xj = \lv{i)k Xk, kmin(i) < v(i) < kmax(i), i = 1, 2 ...m. (37)

Векторы Xj, число которых предполагаете я равным m, представляют собой коэффпцн-енты разложения атомных волновых функций (9): Xj ^ bnl (в примере данной работы X! ^ b10, x2 ^ b20, x3 ^ Ъ21) и считаются упорядоченными по l - орбитальному квантовому числу, так что (Xi; Xj) = öij и li = lj = const, Vi,j G [kmin(i), kmax(i)].

Система (37) представляет собой задачу на собственные значения в некотором обоб-

l

Xj

числа Xj, чтобы вы-Xj

приблизительно 100. Это означает, что необходимо решать систему уравнений с числом переменных от 100 до 1000, что осуществимо только численно, с использованием ЭВМ.

Будем искать решение с помощью итерационного алгоритма, который построим следующим образом. Введем величины

n

qi = AijXj, (38)

j=i

kmax (i)

dqi = qi - (qi, Xk)Xk. (39)

k=kmin (i)

Геометрический смысл невязки dqi заключается в том, что она представляет собой вектор, ортогональный подпространству Span(Xj), определяемому базисом, состоящим из векторов Xj j G [kmin(i), kmax(i)]. Тогда, если набор {Xj}m=1 - решение, то каждый вектор qi лежит только в пространстве Span(Xj), как это следует из (37), и, следовательно, dqi = 0. Таким образом, нужно минимизировать ортогональную невязку (39)

dqi ^ min , i = 1,2,...,m. (40)

Или, что то же самое, найти итерационное решение уравнения

dqi ({Xj}m=i) =0, i = 1, 2,... ,m. (41)

10000 100 1

0.01 о

55 о.ооо1

1е-006

3

й ■8 Й О

0.01

■ ! , , Г 1 , ■ 1 '

• | . 1 \ , 1_

0.1 1 ЯасЦш, а.и.

10

Рис. 1: Электронная плотность атома Ыд.

Задачи типа (41) решаются следующим образом [13, 14].

Пусть имеется некоторое начальное приближение {х0. Будем вычислять последующее приближение по правилу

zkj+1 = хк — ^ , хк+1 = ък+1 ¡\ък+1 ] =1, 2)...,т) (42)

пока процесс не сойдется. Здесь а^ - некоторая малая постоянная матрица с доминирующими диагональными элементами 10"2 ^ 10"5), которая для одних и тех же I одинакова,

К

|х\ = [ £ х2

(43)

г=1

В результате, получающиеся по схеме (42) решения х^ попарно для каждого I ортонор-мированы, а отвечающие им числа как это следует из вариационного принципа, -симметричны по г,].

В качестве начальных значений {х0}т=1 следует выбирать (п^ — I^ — 1)-ые собственные векторы диагональных матриц А^^ (считаем, что нижнее собственное значение индексируется с нуля).

Увеличение скорости сходимости решения следует проводить в соответствие с указаниями работы [15], в которой показано, что в качестве матрицы а^ в уравнениях (42) следует брать диагональную матрицу, элементы которой равны обратным матричным элементам оператора кинетической энергии, вычисленным по базисным функциям (9). Это позволяет подавить неравномерно растущие компоненты итерируемых векторов, и таким образом ускорить сходимость процесса приблизительно в 10 раз.

optimized + natural х

0.1

0.01

^ 0.001 g

w 0.0001

1е-005

1е-006

1е-007

10

100

1000

iter number

Рис. 2: Зависимость величины невязки |dq| от номера итерации для атома Mg.

Кроме этого, иримеиие метода прямого обращения итерированного подпространства (direct inversion of iterative subspace - DUS), изложенного в работе [16], ускоряет сходимость итерационного процесса приблизительно в yU раз, где n - число итераций до ускорения. На рис. 1 и рис. 2 представлены примеры расчетов для атома Mg.

Работа выполнена при частичной финансовой поддержке научных программ Президиума РАН, ОФН РАН и РФФИ (грант N08-02-00757). Авторы выражают искреннюю благодарность Е. Г. Максимову за ценные советы и замечания в ходе подготовки статьи.

[1] С. Froese Fisher, Т. Brage, and P. Jonsson, Computational atomic structure; An MCHF approach. (Institute of phisics publishing, Bristol and Philadelphia 2003).

[2] Д. Хартри Расчеты атомных структур. (Изд. иностранной литературы, Москва, 1960).

[3] А. Бете, Квантовая механика. (Мир, Москва, 1965).

[4] У. Харрисон, Теория твердого тела. (Мир, Москва, 1972).

[5] С. Hartwigsen, S.Goedecker, and J.Hutter, Phys. Rev. В 58, 3641 (1998).

[6] С. Фудзинага, Метод молекулярных орбиталей. (Мир, Москва, 1983).

[7] С. С. Roothan, Rev. Mod. Phys., 23, 69 (1951).

[8] G. Gaigalas and C. Froese Fisher, Сотр. Phys. Comm. 98, 255 (1996).

[9] Л.Д.Ландау, E.M. Лифшиц, Квантовая механика, том III. (Физматлит, Москва,

ЛИТЕРАТУРА

2001).

[10] А.Ф.Никифоров. В.Б.Уваров. Основы, теории специальных функций. (Наука. Москва, 1974).

[11] М. С. Лицарев, О.В.Иванов, ЖЭТФ 138, 28 (2010).

[12] В.Я.Арсенин, Математическая, физика. Основные уравнения, и специальные функции. (Наука, Москва, 1966).

[13] Дж.Трауб, Итерационные методы, решения, уравнений. (Мир, Москва, 1985).

[14] В. М. Вержбицкий, Численные методы,. Линейная, алгебра и нел/иненые уравнения. (Высшая ттткола, Москва, 2000).

[15] М. С. Payne et al., Rev. of Mod. Phys. 64, 1045 (1992).

[16] P.Pulay, Chem. Phys. Lett. 73, 393 (1980).

Поступила в редакцию 19 июля 2010 г.