К теории электронных возбуждений в молекулах на основе уравнения Петерсилки-Госсманна-Гросса

УДК 535.343.2

ББК 22.343

К 56

П.Г. Коваль, Н.Е. Коваль, Д. Фёрстер К теории электронных возбуждений в молекулах на основе уравнения Петерсилки-Госсманна-Гросса

Аннотация

Представлен простой и эффективный метод расчета спектра электронных возбуждений молекул в теории функционала плотности на основе уравнения Петерсилки-Госсманна-Гросса.

Ключевые слова: уравнение Петерсилки-Госсманна-Гросса, функционал электронной плотности.

P.G. Koval, N.E. Koval, D. Foerster On the theory of electronic excitation in molecules on the basis of Petersilka-Gossmann-Gross equation

Abstract

We present a simple and effective method for computation of electronic excitation spectra of molecules. The method is formulated within the density-functional theory, making use of Petersilka-Gossmann-Gross equation for linear density response.

Key words: Petersilka-Gossmann-Gross equation, noninteracting response function.

1. Введение

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

Известно, что в качестве фундаментального объекта можно использовать не только волновую функцию, но и электронную плотность. Основанием для этого (в случае основного состояния) является теорема Хоэнберга-Кона [3, 4], которая утверждает, что внешний потенциал, действующий в шредингеровской системе есть функционал электронной плотности с точностью до несущественной аддитивной постоянной. Так как внешний потенциал определяет гамильтониан системы, то в силу уравнения Шредин-гера, не только внешний потенциал, но и сама волновая функция и все наблюдаемые величины будут функционалами электронной плотности1.

1Теорема Хоэнберга--Кона неоднократно публиковалась как на английском [3-5], так и на русском языках [2, 6].

- 3S -

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

Очевидно, что электронная плотность значительно менее “дорогой” объект по сравнению с волновой функцией, ведь электронная плотность всегда есть функция в трехмерном пространстве времени, а не функция в многомерном гильбертовом пространстве. Более того, плотность, в принципе, измеримая величина, в противоположность волновой функции.

Современная теория функционалов плотности (ТФП) имеет и существенные недостатки. Наиболее серьезный недостаток состоит в том. что неизвестны функционалы плотности измеримых величин. Теорема Хоэнберга-Кона. к сожалению, не дает инструментов для конструирования даже функционала плотности внешнего потенциала -она только доказывает его существование. Для формулирования конструктивной ТФП необходимо знать функционалы энергий: кинетической энергии взаимодействующих электронов, потенциальной энергии электронов во внешнем поле и энергии взаимодействия электронов внутри системы. Если такие функционалы найдены, то энергию и электронную плотность основного состояния можно найти исходя из принципа минимума полной энергии системы.

В настоящее время наиболее популярной (но конечно приближенной) схемой реализации ТФП является схема Кона и Шэма [2. 3. 4]. Схема Кона-Шэма основана на идее эффективных потенциалов, которые действуют в искусственной многоэлектронной системе без взаимодействий и воспроизводят электронную плотность многоэлектронной системы со взаимодействиями. Напомним, что теорема Хоэнберга-Кона утверждает существование однозначного функционала внешнего потенциала от электронной плотности в любой шредингеровской системе. Поэтому можно сделать вывод, что эффективный потенциал, воспроизводящий данную электронную плотность в искусственной системе существует и однозначен.

Таким образом, математической основой формализма Кона-Шэма является уравнение Шредингера из которого изъят потенциал взаимодействия между электронами, но добавлен эффективный потенциал у8(г), призванный воспроизвести плотность системы со взаимодействиями. В этом случае уравнение Шредингера для многоэлектронной системы без взаимодействий распадается на одноэлектронные уравнения. Система этих одноэлектронных уравнений вместе с принятым функционалом эффективного потенциала формируют уравнения Кона-Шэма2

= Егфг

осс

п

г

ув

Эффективный потенциал Кона-Шэма у8(г) разлагается на внешний потенциал (г),

2 В этом уравнении и далее в статье используются атомные единицы, если не указано другое.

потенциал кулоновского отталкивания электронов (потенциал Хартри) г>н(г), потенциал обменного взаимодействия электронов ух(г) и потенциал корреляционного взаимодействия ус(г). Внешний (для электронов) потенциал эквивалентен внешнему потенциалу в формализме Шредингера. Три оставшихся потенциала являются функционалами электронной плотности п(г). Их записывают как вариационные производные соответствующих энергий по плотности

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

Полная энергия моделируемой многоэлектронной системы в схеме Кона-Шэма представляется суммой

а Ееу± - потенциальная энергия электрического заряда с плотностью п(г) во внешнем поле уех1(г). Энергия корреляционного взаимодействия формально определяется в виде

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

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

Статья организована следующим образом: во втором параграфе выводится уравнение Петерсилки-Госсманна-Гросса; в параграфе 1 описывается стандартный расчет спектров в приближении линейного отклика методом Касиды; в параграфе 1 формулируется наш метод расчета спектров, в параграфе 1 приводится сравнение предложенного метода с методом Касиды. В параграфе 1 делаются обобщающие выводы.

Ætot — Т + Eexi + Ец + Ех + Ес,

где Т - кинетическая энергия невзаимодействующих электронов

Ес — Ætot — Т — Eexi — Ец — Ех

[8].

2. Уравнение Петерсилки-Госсманна-Гросса

При малых возмущениях изменение электронной плотности в некоторой точке пространства г и времени t пропорционально возмущающему потенциалу

6п

ÖTtiv t)

где функция = ----—^—- называется функцией отклика электронной си-

8vext {r',t')

стемы на изменение внешнего потенциала 8vext(r,t). Определив эту функцию, мы получаем возможность, например, расчета дипольпого момента d(t) = f dr rSn(r,t) и, следовательно, спектра излучаемых (или поглощаемых) фотонов.

Функция отклика плотности на изменение кон-шэмовского потенциала vs

8vB

8n

где мы обозначили г>нХс = + vx + vc. Определяя функции отклика х = т--т—Xs =

8vext(r,t)'

8n(r', t') SvRxc(r,t) ,

——;----r- и обозначая ——---------------— = tnxc{r,t,r ,t ). мы можем сразу записать уравнение

övs

Петерсилки-Госсманна-Г росса

X-1 = Xs“1 - /Нхс-

Чтобы уяснить смысл оператора инверсии в этом уравнении, подставим, например, вариацию внешнего потенциала Svext(r,t) = J dr'dt' x_1(r, i, г', i')5n(r', i') в уравнение (9). тогда получим

Таким образом, инверсию в уравнении (11) следует понимать в матричном смысле dr,dt,x(r,t,r,,t,)x~1{r,,t,,r”,t") = 8(г — г")8(г - ¿"), а само уравнение (11) нетрудно переписать в виде уравнения Дайсона [9]

Х(г, г', ¡') = х,(г,»,г', *)+]] с1г"си"Лг"'си"'х,(г,«, г", «")/н*с(г", Г, г'", Г)Х(г"', Г, г', (

Собственную энергию /нхс можно приближенно рассчитать, если задана модель обменного и корреляционного потенциалов. В этой статье используется простая модель, для

которой собственные энергии обменного и корреляционного взаимодействий /хс известны как явные функции плотности fxc(r, t, г', t') = S(t — t')S(r — r')/xc(n(r)). Зная определение потенциала кулоповского взаимодействия, можно рассчитать и его собственную энергию fu{r, Í, г', t') = S(t — t')\r — r'|_1.

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

Х{г,г',ш) = Хв{г,г',ш) + JJ dr"dr/"xs

Для изложения метода Касиды (Mark Casida) и нашего метода расчета спектров необходимо перейти к матричному представлению в последнем уравнении (13).

3. Расчет спектров методом Касиды

Функция отклика стационарной кон-шэмовской системы невзаимодействующих электронов на изменение эффективного потенциала в пространстве частот имеет явный вид [4, 5. 10]

Х,(гУ,и,)= £ фЕ(г)МгШг')Мг') уи1_{Е_Р) + к-и; + {Е_Р) + к

Е<0,Р>0 v ’

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

Объявим произведения '0£(г)'гДр(г) новыми базисными функциями и перейдем в уравнении (13) к базису произведений собственных состояний (или базису электроннодырочных возбуждений). Тогда, подставляя кон-шэмовскую функцию отклика (14) и формальное разложение функции отклика х =

ХЕР,Е'Р'Фе{г)Фр(г)Фе'{г)Фе1(г) в уравнение (13). мы получим уравнение (11) и точные выражения для матриц х* и ^нхе

xW'M - Wífí-|^_

КяжЬ'Р'

Производя матричные операции [10], можно показать, что спектр электронных возбуждений представляет собой собственные значения матричного уравнения где матрица задана в пространстве произведений собственных функций уравнения Кона-Шэма. Явный вид ее можно найти во многих монографиях и статьях [4. 5. 10].

/'“ч Е1^ Е1

Для нас важно только запомнить, что матрица содержит ядра в явном ви-

де. Определяя собственные значения и собственные вектора .?/. можно рассчитать

осцилляторные константы переходов //. а затем и спектр динамической поляризуемости Р

Если при расчете ядер (16) мы используем метод вспомогательных функций (auxilary functions) [10] или аналогичный метод приближения плотности (density fitting) [5, 11], то общая сложность расчета спектра методом Касиды будет пропорциональна кубу числа атомов в молекуле (iV3).

4. Расчет спектров методом последовательных приближений

Метод Касиды имеет много достоинств. В частности, он позволяет классифицировать переходы исходя из симметрии состояний, позволяет провести простой расчет осцил-ляторных констант и энергий запрещенных переходов. Благодаря этому метод Касиды стал фактически стандартом в расчетах свойств возбужденных состояний для молекулярных систем - практически все программы для квантово-химических расчетов реализуют метод Касиды [12, 13].

В данной работе мы представляем альтернативный метод расчета динамической поляризуемости Р(ш). Достоинство этого метода заключается в его относительной простоте реализации при такой же масштабируемости вычислительной сложности, как и у метода Касиды. Более того, предлагаемый итеративный расчет поляризуемости был успешно применен в методе с еще более выгодной масштабируемостью вычислительной сложности [14].

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

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

где функция отклика %(г, г', ш) определяется уравнением Петерсилки-Госсманна-Гросса (11). Точные выражения для матриц Xя (ш) и ^Нхс задаются формулами (15, 16). Подставляя разложение функции отклика х(г,г',а>) в базисе произведений собственных функции в уравнение (18), мы приходим к выражению

EF,E'F'

где матричные элементы образуют вектор в пространстве электронно-дырочных возбуждений Е < 0, .Г > 0

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

4.2. Метод подпространств Крылова для поляризуемости

Так как решение линейного уравнения требует 0(М3) операций, где N есть размер матрицы, то очевидно, что прямое решение системы линейных уравнений для вектора Хк возможно только в случае малых молекул. Однако можно заметить, что поляризуемость Ргук(ш) есть скалярное произведение инвертированной матрицы А-1 = (1 — Хо/нхс)-1 на векторы (¿г и Хо^к — т0 есть по сути является усредненной величиной, которую можно легко рассчитать методами подпространств Крылова [14].

Идея метода Крылова состоит в том. чтобы построить простое приближение для матрицы А = ^2^ используя набор векторов ур = Ару, где степень р изменяется

в малом интервале р = 0.1,2 ...Р. Р <С N. а начальный вектор V выбран случайным образом. Пространство, построенное таким образом, позволяет тем более простой расчет обратной матрицы, чем короче наименьший характеристический многочлен <?(£) матрицы А.

Основное свойство характеристического многочлена есть д(А) = ^2р=о арАр = 0. Следовательно,

причем можно показать [17], что, как правило, Р <С N. Отсюда становится понятным, что инвертирование матрицы можно провести с помощью натуральных степеней матрицы и, следовательно, метод построения подпространств Крылова представляется естественным.

В нашем случае мы определяем поляризуемость (20), что проще решения системы линейных уравнений, так как результатом является число Р^(ш), а не вектор Хк- Метод подпространств Крылова в виде предложенном самим Крыловым часто неэффективен, так как базисные векторы ур = Ару быстро становятся линейно-зависимыми и, следовательно, бесполезными. Арнольди предложил модифицированный метод Крылова [18], в котором каждый последующий вектор |'«р+1) претерпевает дополнительную явную ортогонализацию к предыдущим векторам \ур). например по методу

Грама-Шмидта. Если матрица А является самосопряженной (А = Ат*), то процесс ортогонализации нужно проделать только для двух предыдущих векторов \ур). 1^-1), после чего вектор |'Ур+1) становится ортогональным ко всем предыдущим векторам автоматически. Проиллюстрируем этот факт на простейшем примере с использованием ортогонализации Грама-Шмидта. Построим подпространство Крылова |г>п) по правилу

Здесь вектор («п+х) выбран ортогональным только к двум предыдущим векторам |г>п) и ^п-х), однако будет ортогонален и к |г>п_2). Проверяя это утверждение, вычислим (г)та_2|«та+1). Согласно рекуррентному правилу вектор ортогонален к векторам

|г>п_1) и |г>п), поэтому

Р-1

К+1> = Ии„) - К)<г1„|.4и„) - |1>о) = 0.

Выразим \Avn_2) из рекуррентного соотношения (21) для |fn-i), тогда получим

(vn_2\vn+l) = (vn\vn-i) + {vn\vn_2)(vn_2\A\vn_2) + (vn\vn_3)(vn_3\A\vn_2).

То есть, мы свели скалярное произведение {vn_2\vn+i) к (vn-3\vn)- Повторяя это рассуждение мы приходим к выводу, что {vn_2\vn+i) пропорционально (г)о|«з) и. следовательно, равно нулю.

В нашем случае матрица А не самосопряженная, поэтому мы применили более общий биортогональный метод Ланцоша (Cornelius Lanczos) [18]. Матрица А = (1~хо^) будет представлена в виде

ч

где начальные векторы и (г^х| выбраны в направлении |хо^} и (^| и нормированы согласно (^[гУх) = 1 :

|-0- }Х°*> М

л/Щхой' у/Щхо\Ф

Последующие векторы |г>п) и (го„| получены при помощи рекуррентного соотношения [18]. которое ортогопалнзирует правые и левые векторы = 8ц. Таким образом,

выражение для поляризации будет иметь вид

рЫ =

ч

В силу биортогональности базисных векторов, сумма по ним сводится к единственному члену г = .7 = 1, а поляризуемость Р(ш) будет пропорциональна поляризуемости Ро (ш) = (в,\xo\d)-, вычисленной при помощи функции отклика Хо для невзаимодействующих электронов: Р(ш) = Ро(со)^~1)ц. Эффективность метода целиком зависит от размера матрицы = (т^А^^ по сравнению с матрицей А. В наших примерах мы достигаем приемлемой точности для Р(ш) уже с маленькими подпространствами Крылова А^к ~ 10 для всех рассмотренных молекул от метана СН4 (размер матрицы А, Ид = 120) до нафталина СюН8 (размер матрицы А, Ид = 3504).

4.3. Расчет матриц взаимодействия Кн и Кхс

Определим собственную энергию /нхе в базисе произведений собственных состояний Фе(г)Фе(г). Для этого мы должны рассчитать интегралы

КГ-Е'р' = [[ ;

KxcF,E'F' = Jj d3rd3rf фЕ

Собственные состояния Фе(г) являются решениями уравнений Кона-Шэма для определенной энергии Е. При этом для расчета собственных состояний мы используем широко известный метод линейной комбинации локализованных атомных орбиталей fa{r)-В результате диагонализации гамильтониана Кона-Шэма собственные состояния получаются в виде

Фе

а

- 45 -

где Хр - собственные вектора, а /а(г) - атомные орбитали — элементарные функции. обладающие сферической симметрией /а(г) = /га(у>)5'гата(г) — произведение радиальной функции и сферической гармоники3. Подставляя собственные состояния (25) в уравнение (23). мы можем выразить ядро кулоновского взаимодействия К-ц через так называемые четырехцентровые интегралы

(оЬ\с<1) = ЦЛ?'гЛ3г' /„(гЩг) 1 /с(г')1а(г'),

КЕР,Е'Р' = '£x*X[X?Xf

abed,

Прямой расчет собственной энергии Хартри по формуле (26) чрезвычайно неэффективен. Для эффективного расчета матриц К-ц и Кхс используют либо метод вспомогательных функций (auxiliary functions) [10], либо метод аппроксимации плотностей (density fitting) [5, 11]. Мы используем аналогичный метод доминирующих произведений [16]. Метод доминирующих произведений позволяет разложить произведение локализованных орбиталей fa{r)fb{r) в базисе ортогональных локализованных функций F^{r) — доминирующих произведений

/,МЛ

Учитывая это разложение, мы получим выражение для ядра Хартри

КЕР,Е'Р' = Y^

abedfiv

fs = Jf*rtуМр.

Обозначая WffF = , мы можем переписать уравнение (28) в матричных

ЕF Е1 Ff т1

обозначениях Кн ’ = Wf-^W . Для собственной энергии обменно-корреляционного

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

/£- = J d3rF^

В завершение следует рассчитать интегралы и

Доминирующие произведения конструируются отдельно для каждой пары атомов. Поэтому возможны два вида произведений: локальные (когда атомы в паре совпадают) и билокальные (когда атомы в паре различаются). В случае локальных произведений, доминирующие функции F^{r) обладают сферической симметрией F^(r) = F^^Sj m (г), так как являются линейной комбинацией исходных сферически-симметричных орбиталей. В случае билокальных произведений, доминирующие функции F^{r) сохраняют только аксиальную симметрию. То есть выбирая систему координат с осью z

3Для лучшего быстродействия программы мы используем действительные сферические гармоники.

вдоль линии соединяющей атомы, мы можем сохранить магнитное число Мд для данной доминирующей функции ^(г) = ^^(г)^^ (Иг). Здесь матрица поворота И

преобразует координаты из глобальной системы координат в систему координат, связанную с данной парой атомов. Используя выражения для локальных и билокальных доминирующих функций, мы можем свести кулоновский интеграл (29) к суммам по элементарным кулоновским интегралам [15. 19]

где Fj(r) - радиальные функции, a Yjm(r) (комплексные) сферические гармоники. В нашем случае, радиальные функции доминирующих произведений F^(r) и Fl,J'(r) заданы численно на логарифмической сетке, поэтому для вычисления кулоновского интеграла (31) мы будем следовать методу Талмана (James Taiman) [15]. Суть метода состоит в переходе от координатного к импульсному представлению в этом интеграле. Записывая Fjm(r) = Fj(r)Yjm{r) как обратное преобразование Фурье от соответствующей функции в пространстве импульсов Fjm(p) и используя разложения кулоновского взаимодействия и 5-функции Дирака в ряд Фурье, получим

Для преобразования радиальных функций из пространства координат в пространство импульсов Fj(r) —> Fj(p) мы применяем быстрое преобразование Бесселя [15]. связанное с быстрым преобразованием Фурье. В результате, расчет матрицы может быть проведен за 0(./V2) математических операций и занимает, например, не более 2 минут для молекулы бензола на современном персональном компьютере.

Расчет собственной энергии обменно-корреляционного взаимодействия лучше провести в координатном пространстве, так как интеграл (30) локальный. Доминирующие функции ^(г) отличны от нуля только в области пространства в виде линзы, а произведение двух функций ^(г)^(г) будет отлично от нуля в области, которая является перекрытием двух линз. Несмотря на это. мы нашли достаточным использование обычных сферических координат с центром в середине линии, соединяющей центры ^(г) и ^(г). Используя для интегрирования по телесному углу метод Лебедева [20]. а для интегрирования по радиальной координате метод Гаусса-Лежандра, мы можем проводить расчет интеграла только в порядке тысячи точек, сохраняя приемлемую точность результата. Погрешность интегрирования оценивалась по интегральному критерию - сходимости поляризуемости Р(ш). Сложность вычисления собственной энергии обменно-корреляционного взаимодействия растет как первая степень числа атомов N в силу его локальности в координатном пространстве. При расчете молекулы бензола наша программа требует порядка 3 минут.

4.4. Вычислительная сложность алгоритма

Оценим количество операций для расчета поляризуемости Р(ш) итеративным способом. Как мы показали выше, расчет матриц кулоновского взаимодействия /н в базисе доминирующих произведений может быть проведен за 0(М2) операций, а расчет матриц /хс в этом же базисе - за О(М) операций. В то же время функция отклика для невзаимодействующих электронов ХоР’Е Р требует 0(М2) операций, так как является диагональной матрицей и пространство электронных возбуждений Е < 0, Е > 0

(.jmj’m') =

drdr'Fj(r)Yjm(r)\r - г'\ 1Е0,(г')¥0,т,(г'),

имеет размер Л^осс-^у1й ~ ^2. Для построения подпространства Крылова мы будем последовательно умножать матрицы ^нхс и Xя на векторы |^¿), |г^). Умножение вектора на диагональную матрицу I) происходит за 0(N0) математических операций, то есть 0(М2) в нашем случае. Используя представление собственной энергии ^нхс в виде К-цхс = №/лхсЦгТ, мы потратим 0(М3) операций, так как IV является прямоугольной матрицей (]У2 х Л/”). Расчет самой матрицы IV можно провести так же за 0(М3) операций \¥рр = ХРХРУ£ь, так как вертекс У£ь является “пустым” массивом, для которого только несколько элементов отличны от нуля для данного индекса ¡1. Следовательно расчет поляризуемости предложенным методом можно провести за NzNкNU} операций, где N^3 - это количество частот, для которых нужна поляризуемость, а Nк - размерность подпространства Крылова. Проведенное нами компьютерное моделирование показывает. что Нк может быть выбрано значительно меньшим по сравнению с размером исходных матриц и практически независящим от количества атомов N. Таким образом, сложность алгоритма будет пропорциональна Л^А^. Практически наша программа затрачивает порядка 7 минут для спектра бензола, состоящего из 256 частот и порядка 30 минут для спектра нафталина такого-же разрешения.

Предложенный итерационный метод расчета поляризуемости может быть с успехом применен в комбинации с другими базисами, нежели предлагаемый в этой статье -базис электронно-дырочных возбуждений. Например, использование базиса доминирующих произведений и для расчета кон-шэмовской функции отклика и и для расчета собственных энергий дает возможность построить алгоритм еще меньшей сложности 0(7У2ЛУ [14].

5. Примеры спектров электронных возбуждений

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

5.1. Бензол

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

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

э

с

н

со, еУ

Рис. 1. Поляризуемость молекулы бензола. Экспериментальные данные взяты из [23]. В расчетах использован базис с двойными поляризованными орбиталями Б2УР доступный по умолчанию в программном пакете с1еМоп2к [12]. В качестве потенциала обменно-корреляционного взаимодействия выбрана параметризация Пердева-Цунгера [7] данных Монте-Карло симуляций Цеперли и Алдера [8]. Теоретические расчеты проведены стандартным методом Касиды. и предлагаемым в данной статье методом. Эти результаты очень близки друг к другу.

следствием электронных переходов между связывающими и антисвязывающими 7г орбиталями (7г —> 7г*) углеродного кольца [23]. в то время как непрерывный спектр в районе 12 — 30 эВ есть следствие переходов между связывающими и антисвязывающими а орбиталями. Теория предсказывает несколько большую энергию перехода для коллективного пика (7 эВ). что представляется хорошим совпадением. Непрерывный спектр переходов в высокочастотной области 12 — 30 эВ невозможно воспроизвести с базисом, состоящим из орбиталей дискретного спектра, поэтому мы видим набор резонансов. качественно повторяющий экспериментальную поляризуемость.

5.2. Пиррол

На рисунке 2 представлены результаты двух теоретических расчетов и экспериментальные данные для поляризуемости молекулы пиррола. Оба теоретических расчета выполнены нашим методом с использованием базиса DZVP (95 орбиталей) и Т7УР (130 орбиталей). Экспериментальные данные взяты из [24]. подробный анализ природы переходов проведен там же. Наши расчеты качественно свидетельствуют о наличии низкочастотной и высокочастотной зон спектра. Относительная интенсивность низкочастотной и высокочастотной зон также воспроизводится удовлетворительно. Вместе с тем в наших расчетах отсутствует среднечастотная зона, где наблюдаются хорошо выраженные электронные переходы.

Можно предположить, что для получения хороших результатов в среднечастотной области необходимо использовать еще более расширенный базис электронных орбита-лей. как это сделано в [24]. В подтверждение этой гипотезы, мы сравниваем на рисунке 3 поляризуемость, построенную по теоретически полученным данным из [24] с нашим расчетом.

<2 Ь5 е

Рис. 2. Поляризуемость молекулы пиррола. Экспериментальные данные взяты из [24]. В расчетах использован базис с двойными (DZVP) и тройными (TZVP) поляризованными орбиталями, доступными в программном пакете с1еМоп2к [12. 22]. В качестве потенциала обменно-корреляционного взаимодействия выбрана параметризация Пердева-Цунгера [7] данных Монте-Карло симуляций Цеперли и Алдера [8]. Использование расширенного базиса TZVP не позволило описать спектр в среднечастотной зоне, но скорректировало положение

со. еУ

первого резонанса в правильном направлении.

Our method TZVP ---------

SAC-CI method .......

Experimental -------

5 5.5 6 6.5 7 7.5 8 8.5

0), eV

Рис. 3. Поляризуемость молекулы пиррола. Вторая теоретическая кривая построена с использованием данных из [24]. Видно, что базис, использованный в [24] в сочетании с методом адаптированного по симметрии кластерного конфигурационного взаимодействия (symmetry adapted cluster-configuration interaction method), воспроизводит спектр электронных возбуждений более качественно.

6. Выводы

В настоящей статье представлен простой метод расчета поляризуемости, применимый к молекулярным системам. Основой метода является уравнение Петерсильки-Госсманна-Гросса для функции отклика. Базис доминантных произведений и итеративный метод расчета поляризуемости позволяют проводить расчет поляризуемости, затрачивая 0(N^N^) операций. Предложенный итеративный метод расчета поляризуемости был использован в конструировании альтернативного метода расчета поляризуемости с еще более выгодным ростом сложности 0(М2Мш).

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

Благодарности

Работа была выполнена в рамках проекта ANR NOSSI в университете Бордо.

Мы выражаем благодарность Джеймсу Талману (Лондон. Канада) и Марку Касиде (Гренобль. Франция) за плодотворные дискуссии, необходимые алгоритмы и программы; Андрею Постникову (Метц. Франция) за советы по использованию программного пакета Siesta; Михаилу Маматкулову (Марсель. Франция). Ланжу Шекояну (Москва. Россия) за критические замечания к препринту данной работы.

Приложение: функция отклика для одного электрона

Чтобы обосновать выражение (14) для функции отклика системы невзаимодействующих электронов, приведем элементарный вывод функции отклика одно-электронной системы. Функция отклика есть вариационная производная от электронной плотности по потенциалу xo(r,t,r',t') = Sn(r,t)/SV(r',t'). Поэтому мы рассчитаем волновую функцию ф(гЛ) для возмущения SV(r.t). а затем и изменение электронной плотности 8n(r,t) = 8(ф*(г,*)ф(г,г)).

Пусть возмущенная функция ф(г^) удовлетворяет уравнению Шредингера

Будем искать эту функцию в виде ф(г^) = фо(г^) + 5ф(г^), где фо(г^) - волновая

функция певозмущеппого состояния удовлетворяет уравнению

Подставляя ф(г^) = фо(г^) +5ф(г^) в уравнение (33). мы получим для поправки 5ф(г^) следующее неоднородное уравнение

5ф

Решить последнее уравнение можно стандартным методом функции Грина. Невозму-щеппая функция Грина (зависящая от времени) удовлетворяет уравнению

- H0^j Go

Функция Грина есть ядро интегрального оператора обратного — //о) и может быть использована для решения уравнения

8ф{г,£) = с13г'сИ' Сго

Записывая далее изменение электронной плотности, мы получим

8п

Поэтому функция отклика х(г,£,г',£') будет

Х{т,г,г',1!) = ^ =фп

Выражение для функции Грина через сумму по всем собственным состояниям фоп(г, невозмугценного уравнения Шредингера необходимо в случае не зависящего от времени гамильтониана Но- В этом случае уравнение (36) можно регуляризовать путем добав-

(.д \

ления малой мнимой константы к оператору 1—— Но . Прямой подстановкой

дЬ

можно убедится. ЧТО Со(г, г', ¿') имеет вид [26]

= -1^2ф0п

,, . . 0,ж>0 ( л Г1, ж > 0

где 0+(ж) = < 1 „ , 0 (ж) = < „ „ .

4 ' 1 -1,х < 0 ’ к ' 1 0,х < 0

Тогда, подставляя явное выражение для функции Грина (40) в формулу (39). полу-

чим

х(г,г,г',1?) = -[ф*{г,г)^2ф0п(г,ф*п{г',г')е±^ *,)0±(г-г/)^о(»Л*/)+

п

+ 1ф0

п

Из этого выражения уже явно видно, что функция отклика х(г, г', £') зависит только от разницы 1 — ^. а не от обоих аргументов по отдельности. Поэтому, мы можем найти преобразование Фурье для одного аргумента частоты. Выполняя преобразование Фурье для функции отклика (41). получим

Х(г г' ш) = ±У(_ Мг)Фп(г)Мг'Ш(г')\ ,42ч

^ \ ш - (Е0~ Еп)±1е ш + (Е0 - Еп) ±1в ) '

Этот результат можно сравнить с выражением для функции отклика (14). Кроме этого. приведенный вывод дает возможность критически оценить выражение для функции отклика (14). Именно сумма по полному спектру ^п в полученной формуле (42) заменяется в выражении (14) на сумму по дискретному спектру, полученному путем прямой диагопализации кон-шэмовского гамильтониана. Состояния непрерывного спектра моделируются при этом с помощью связанных состояний, что ведет к плохому совпадению с экспериментальными спектрами. Лучшее описание непрерывного спектра может быть достигнуто, например, в рамках прямого решения, зависящего от времени уравнения Кона-Шэма или при использовании методов, основанных на решении уравнения определяющего функцию Грина (36).

Примечания:

1. Ландау Л.Д., Лифшиц ЕМ., Квантовая механика (нерелятивистская теория). М.: Физматлит, 2004. 800 с.

2. Кон В. Электронная структура вещества - волновые функции и функционалы плотности // УФН. 2002. Том 172, № 3. С. 336-348.

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

4. A Primer in Density Functional Theory / ed. by C. Fiolhais, F. Nogueira, M.A.L. Marques. Berlin: Springer, 2003. Vol. 620.

5. Time-Dependent Density Functional Theory / ed. by M.A.L. Marques, C.A. Ullrich, F. Nogueira, A. Rubio, K. Burke, E.K.U. Gross. Berlin: Springer, 2008. Vol. 706.

6. , -

ры ряда углеродных и неуглеродных нанокластеров и их взаимодействия с легкими элемен-/ . . , . . , . . -мов, СТ. Овчинников. Новосибирск: СО РАН, 2006. URL: http://www.kirensky.ru/master/arti-cles/monogr/Book/

7. Perdew J.P., Zunger A. Self-interaction correction to density-functional approximations for many-electron systems // Phys. Rev. 1981. Vol. B 23. P. 5048-5079.

8. Ceperley D.M., Alder B.J. Ground State of the Electron Gas by a Stochastic Method // Phys. Rev. Lett. 1980. Vol. 45. P. 566-569.

9. Petersilka M., Gossmann U.J., Gross E.K.U. Excitation Energies from Time-Dependent Density-Functional Theory // Phys. Rev. Lett. 1996. Vol. 76. P. 1212-1215.

10. Casida M.E. Time-Dependent Density Functional Response Theory for Molecules // Recent Advances in Density Functional Theory / ed. by D.P. Chong. Singapore: World Scientific, 1995. P. 155-193.

11. Chemistry with ADF / G.Te Velde, F.M. Bickelhaupt, E.J. Baerends, C. Fonseca Guerra, S.J.A. Van Gisbergen, J.G. Snijders, T. Ziegler // J. of Comput. Chem. 2001. Vol. 22, Iss. 9. P. 931-967.

12. The Grenoble development version of deMon2k

which is based upon deMon2k, version 2.2 / A.M. Koster, P. Calaminici, M.E. Casida, R. FloresMoreno, G. Geudtner, A. Goursot, T. Heine, A. Ipatov, F. Janetzko, J.M. del Campo, S. Patchovskii, J.U. Reveles, D.R. Salahub, A. Vela. Mexico: Cin-vestav, 2006. URL: http://dcm.ujf-

greno-

ble.fr/PERSONNEL/CT /casida/deMonaGrenoble

13. Scientific Computing & Modelling (SCM) NV

corporation - molecular modeling company -http://www.scm.com; Time-dependent density functional theory (TDDFT) -

http://www.tddft.org}; The Official Gaussian Website - http://www.gaussian.com

References:

1. Landau L.D., Lifshits E.M. The quantum mechanics (nonrelativistic theory). M.: Fizmatlit, 2004. 800 pp.

2. Kon V. Electronic substance structure - wave functions and density functionals // UFN. 2002. Vol. 172, No. 3. P. 336-348.

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

4. A Primer in Density Functional Theory / ed. by C. Fiolhais, F. Nogueira, M.A.L. Marques. Berlin: Springer, 2003. Vol. 620.

5. Time-Dependent Density Functional Theory / ed. by M.A.L. Marques, C.A. Ullrich, F. Nogueira, A. Rubio, K. Burke, E.K.U. Gross. Berlin: Springer, 2008. Vol. 706.

6. Modelling of properties, electronic structure of some carbon and noncarbon nanoclasters and their interaction with easy elements / A. S. Feodorov, P.B.Sorokin, P.V.Avramov, S.G. Ovchinnikov. Novosibirsk: the Siberian Branch of the Russian Academy of Science, 2006. URL: http://www.kirensky.ru/master/articles/monogr/Book/

7. Perdew J.P., Zunger A. Self-interaction correction to density-functional approximations for many-electron systems // Phys. Rev. 1981. Vol. B 23. P. 5048-5079.

8. Ceperley D.M., Alder B.J. Ground State of the Electron Gas by a Stochastic Method // Phys. Rev. Lett. 1980. Vol. 45. P. 566-569.

9. Petersilka M., Gossmann U.J., Gross E.K.U. Excitation Energies from Time-Dependent Density-Functional Theory // Phys. Rev. Lett. 1996. Vol. 76. P. 1212-1215.

10. Casida M.E. Time-Dependent Density Functional Response Theory for Molecules // Recent Advances in Density Functional Theory / ed. by D.P. Chong. Singapore: World Scientific. 1995. P. 155-193.

11. Chemistry with ADF / G.Te Velde, F.M. Bickelhaupt, E.J. Baerends, C. Fonseca Guerra, S.J.A. Van Gisbergen, J.G. Snijders, T. Ziegler // J. of Comput. Chem. 2001. Vol. 22, Iss. 9. P. 931-967.

12. The Grenoble development version of deMon2k

which is based upon deMon2k, version 2.2 / A.M. Koster, P. Calaminici, M.E. Casida, R. FloresMoreno, G. Geudtner, A. Goursot, T. Heine, A. Ipatov, F. Janetzko, J.M. del Campo, S. Patchovskii, J.U. Reveles, D.R. Salahub, A. Vela. Mexico: Cin-vestav, 2006. URL: http://dcm.ujf-

greno-

ble.fr/PERSONNEL/CT/casida/deMonaGrenoble

13. Scientific Computing & Modelling (SCM) NV

corporation - molecular modeling company -http://www.scm.com; Time-dependent density functional theory (TDDFT) -

http://www.tddft.org}; The Official Gaussian Website - http://www.gaussian.com

14. On the Kohn-Sham density response function / D. Foerster, P. Koval // J. Chem. Phys. 2009. Vol. 131. 044103.

15. Numerical calculation of four-center Coulomb integrals // J. Chem. Phys. 1984. Vol. 80. Iss. 5. 2000; Talman J.D. Numerical Fourier and Bessel transforms in logarithmic variables // Journal of Computational Physics. 1978. Vol. 29. Iss. 1. P. 35-48; Talman J.D. LSFBTR: A subroutine for calculating spherical bessel transforms // Computer Physics Communications. 1983. Vol. 30. Iss. 1. P. 93-99; Talman J.D. NumSBT: A subroutine for calculating spherical Bessel transforms numerically // Computer Physics Communications. 2009. Vol. 180. Iss. 2. P. 332-338.

16. Foerster D. On an «interaction by moments» property of four center integrals // J. Chem. Phys. 2008. Vol. 128; Talman J.D. The implementation of the method proposed in this paper has been improved thanks to new algorithm // J. Quant. Chem. 2007. Vol. 107. P. 1578.

17. Ipsen I.C.F., Meyer C.D. The Idea behind Krylov Methods // The American Mathematical Monthly. 1998. Vol. 105. P. 889.

18. Saad Y. Iterative Methods for Sparse Linear Systems. Philadelphia: Siam. 2003. 552 p.; Баландин М.Ю., Шурина ЭЛ. Методы решения СЛАУ большой размерности. Новосибирск:

- , 2000. 70 .

19. Blanco M.A., Florez M., Bermejo M. // J. Mol. Struct. (THEOCHEM). 1997. Vol. 419. P. 19-27.

20. Lebedev V.I. et al. // Rus. Acad. Sci. Dokl. Math.

1995. URL: http://www.ccl.net/cca/software/

SOURCES/FORTRAN/Lebedev-Laikov-Grids/

21. Godbout N., Salahub D.R., Andzelm J., Wimmer E. Optimization of Gaussian-type basis sets for local spin density functional calculations. Part I. Boron through neon, optimization technique and validation // Canadian Journal of Chemistry. 1992. Vol. 70, No 2. P. 560-571.

22. URL:http://www.emsl.pnl.gov/forms/basisform.html

23. Yabana K., Bertsch G.F. Time-dependent local-density approximation in real time: Application to conjugated molecules // International Journal of Quantum Chemistry. 1999. Vol. 75. Iss. 1. P. 55-66. URL: http://arxiv.org/abs/physics/9808015

24. Electronic excitation spectra of furan and pyrrole: Revisited by symmetry adapted cluster-configuration method / J. Wan, J. Meller, M. Hada, M. Ehara, H. Nakatsuji // The Journal of Chemical Physics. 2000. Vol. 113. Iss. 18. P. 7853.

25. - . URL:

http://en. citizendium.org/wiki/Rayleigh-Ritz_method

26. Harrison W.A. Solid State Theory. N.Y.: McGraw-Hill. 1970; имеется перевод: Харрисон У. Теория твердого тела. М.: Мир, 1972. 616 с.

14. On the Kohn-Sham density response function / D. Foerster, P. Koval // J. Chem. Phys. 2009. Vol. 131. 044103.

15. Numerical calculation of four-center Coulomb integrals // J. Chem. Phys. 1984. Vol. 80. Iss. 5. 2000; Talman J.D. Numerical Fourier and Bessel transforms in logarithmic variables // Journal of Computational Physics. 1978. Vol. 29. Iss. 1. P. 35-48; Talman J.D. LSFBTR: A subroutine for calculating spherical bessel transforms // Computer Physics Communications. 1983. Vol. 30. Iss. 1. P. 93-99; Talman J.D. NumSBT: A subroutine for calculating spherical Bessel transforms numerically // Computer Physics Communications. 2009. Vol. 180. Iss. 2. P. 332-338.

16. Foerster D. On an «interaction by moments» property of four center integrals // J. Chem. Phys. 2008. Vol. 128; Talman J.D. The implementation of the method proposed in this paper has been improved thanks to new algorithm // J. Quant. Chem. 2007. Vol. 107. P. 1578.

17. Ipsen I.C.F., Meyer C.D. The Idea behind Krylov Methods // The American Mathematical Monthly. 1998. Vol. 105. P. 889.

18. Saad Y. Iterative Methods for Sparse Linear Systems. Philadelphia: Siam. 2003. 552 p.; Баландин М.Ю., Шурина ЭЛ. Методы решения СЛАУ большой размерности. Новосибирск:

- , 2000. 70 .

19. Blanco M.A., Florez M., Bermejo M. // J. Mol. Struct. (THEOCHEM). 1997. Vol. 419. P. 19-27.

20. Lebedev V.I. et al. // Rus. Acad. Sci. Dokl. Math.

1995. URL: http://www.ccl.net/cca/software/

SOURCES/FORTRAN/Lebedev-Laikov-Grids/

21. Godbout N., Salahub D.R., Andzelm J., Wimmer E. Optimization of Gaussian-type basis sets for local spin density functional calculations. Part I. Boron through neon, optimization technique and validation // Canadian Journal of Chemistry. 1992. Vol. 70, No 2. P. 560-571.

22. URL:http://www.emsl.pnl.gov/forms/basisform.html

23. Yabana K., Bertsch G.F. Time-dependent local-density approximation in real time: Application to conjugated molecules // International Journal of Quantum Chemistry. 1999. Vol. 75. Iss. 1. P. 55-66. URL: http://arxiv.org/abs/physics/9808015

24. Electronic excitation spectra of furan and pyrrole: Revisited by symmetry adapted cluster-configuration method / J. Wan, J. Meller, M. Hada, M. Ehara, H. Nakatsuji // The Journal of Chemical Physics. 2000. Vol. 113. Iss. 18. P. 7853.

25. Rayleigh-Ritz method. URL:

http://en.citizendium.org/wiki/Rayleigh-Ritz_method

26. Harrison W.A. Solid State Theory. N.Y.: McGraw-Hill. 1970; trans.: Harrison W.A. Solid State Theory. М.: Mir, 1972. 616 pp.