CC BY
36
УДК 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
= Егфг(г), (1)
осс
п(г) = ^Шг)\2; (2)
г
ув(г) = Уехь(г) + ин(У) + Ух(г) + Ус(г). (3)
Эффективный потенциал Кона-Шэма у8(г) разлагается на внешний потенциал (г),
2 В этом уравнении и далее в статье используются атомные единицы, если не указано другое.
потенциал кулоновского отталкивания электронов (потенциал Хартри) г>н(г), потенциал обменного взаимодействия электронов ух(г) и потенциал корреляционного взаимодействия ус(г). Внешний (для электронов) потенциал эквивалентен внешнему потенциалу в формализме Шредингера. Три оставшихся потенциала являются функционалами электронной плотности п(г). Их записывают как вариационные производные соответствующих энергий по плотности
Обменный и корреляционный потенциалы вычисляются в некотором приближении. Например. в данной статье используются потенциалы, представляющие собой явные функции от локального значения плотности, а функционалы энергии будут заданы через квадратуры функций плотности. Разделение обменного и корреляционного потенциалов на два потенциала пояснено ниже.
Полная энергия моделируемой многоэлектронной системы в схеме Кона-Шэма представляется суммой
а Ееу± - потенциальная энергия электрического заряда с плотностью п(г) во внешнем поле уех1(г). Энергия корреляционного взаимодействия формально определяется в виде
Она абсорбирует в себя разницу кинетической энергии системы взаимодействующих электронов и кинетической энергии системы невзаимодействующих электронов, а так же неточность в функционале обменного взаимодействия. На практике используются приближенные конструкции. Например, в данной статье мы рассчитываем корреляционный потенциал, используя параметризацию Пердева и Цунгера (John Perdew. Alex Zunger) [7] для Монте Карло расчетов Цеперли и Алдера (David Ceperley. Berni Alder)
Решающим элементом схемы Кона-Шэма является приближенная замена кинетической энергии взаимодействующих электронов на кинетическую энергию невзаимодействующих электронов. Уравнения Кона-Шэма с успехом применяются для многих электронных систем, начиная от атомных и заканчивая твердыми телами. В настоящей статье мы применим формализм Кона-Шэма для молекул и рассмотрим простой метод расчета спектров электронных возбуждений в приближении линейного отклика.
Статья организована следующим образом: во втором параграфе выводится уравнение Петерсилки-Госсманна-Гросса; в параграфе 1 описывается стандартный расчет спектров в приближении линейного отклика методом Касиды; в параграфе 1 формулируется наш метод расчета спектров, в параграфе 1 приводится сравнение предложенного метода с методом Касиды. В параграфе 1 делаются обобщающие выводы.
(5)
Ætot — Т + Eexi + Ец + Ех + Ес,
(6)
где Т - кинетическая энергия невзаимодействующих электронов
(7)
Ес — Ætot — Т — Eexi — Ец — Ех
(8)
[8].
2. Уравнение Петерсилки-Госсманна-Гросса
При малых возмущениях изменение электронной плотности в некоторой точке пространства г и времени t пропорционально возмущающему потенциалу
6п(г,г)= / е/г'сй' (9)
ÖTtiv t)
где функция = ----—^—- называется функцией отклика электронной си-
8vext {r',t')
стемы на изменение внешнего потенциала 8vext(r,t). Определив эту функцию, мы получаем возможность, например, расчета дипольпого момента d(t) = f dr rSn(r,t) и, следовательно, спектра излучаемых (или поглощаемых) фотонов.
Функция отклика плотности на изменение кон-шэмовского потенциала vs(r,t) для КОН-ШЭМОВСКОЙ многоэлектронной системы невзаимодействующих электронов Xs{ri r' 1 ш) имеет простой вид [5] (см. также ниже формулу (14) и приложение). Однако, функция отклика Хз не эквивалентна отклику на изменение внешнего потенциала %. а уравнение Петерсильки-Госсманна-Гросса (Martin Petersilka. Ullrich Gossmarm. Eberhard Gross) связывает функции отклика х и Xs- Чтобы вывести это уравнение, запишем вариационную производную от эффективного потенциала (3) по плотности n(r,t)
8vB(r,t) _8vext(r,t) 8vKxc(r,t)
8n(r',t') 8n(r',t') 8n(r',t')
где мы обозначили г>нХс = + vx + vc. Определяя функции отклика х = т--т—Xs =
8vext(r,t)'
8n(r', t') SvRxc(r,t) ,
——;----r- и обозначая ——---------------— = tnxc{r,t,r ,t ). мы можем сразу записать уравнение
övs(r,t) 8n(r',t')
Петерсилки-Госсманна-Г росса
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г"си"Лг"'си"'х,(г,«, г", «")/н*с(г", Г, г'", Г)Х(г"', Г, г', (
(12)
Собственную энергию /нхс можно приближенно рассчитать, если задана модель обменного и корреляционного потенциалов. В этой статье используется простая модель, для
которой собственные энергии обменного и корреляционного взаимодействий /хс известны как явные функции плотности 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(r,r",D)fIlxc(r",r"l)x(r,",r,,u). (13)
Для изложения метода Касиды (Mark Casida) и нашего метода расчета спектров необходимо перейти к матричному представлению в последнем уравнении (13).
3. Расчет спектров методом Касиды
Функция отклика стационарной кон-шэмовской системы невзаимодействующих электронов на изменение эффективного потенциала в пространстве частот имеет явный вид [4, 5. 10]
Х,(гУ,и,)= £ фЕ(г)МгШг')Мг') уи1_{Е_Р) + к-и; + {Е_Р) + к
Е<0,Р>0 v ’
(14)
где функции Фе{^) - собственные функции уравнения Кона-Шэма. Уровни энергии Е, Е сдвинуты таким образом, что состояния с отрицательными энергиями заполнены. а с положительной энергией не заполнены (виртуальные состояния). Малая положительная константа е предотвращает пулевые значения знаменателя и приводит к лоренцевской форме резонансов.
Объявим произведения '0£(г)'гДр(г) новыми базисными функциями и перейдем в уравнении (13) к базису произведений собственных состояний (или базису электроннодырочных возбуждений). Тогда, подставляя кон-шэмовскую функцию отклика (14) и формальное разложение функции отклика х =
ХЕР,Е'Р'Фе{г)Фр(г)Фе'{г)Фе1(г) в уравнение (13). мы получим уравнение (11) и точные выражения для матриц х* и ^нхе
xW'M - Wífí-|^_(b_í,) + í£. а, + (я_^) + к)' <15)
КяжЬ'Р'(и) = JJ (¡Ыг'фЕ{г),фР(г)КВхс(г,г')'фЕ'{г),фр'(г). (16)
Производя матричные операции [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), мы приходим к выражению
(17)
(18)
(19)
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) и |г>п), поэтому
(20)
Р-1
К+1> = Ии„) - К)<г1„|.4и„) - |1>о) = 0.
(21)
Выразим \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ц. Таким образом,
выражение для поляризации будет иметь вид
рЫ = (22)
ч
В силу биортогональности базисных векторов, сумма по ним сводится к единственному члену г = .7 = 1, а поляризуемость Р(ш) будет пропорциональна поляризуемости Ро (ш) = (в,\xo\d)-, вычисленной при помощи функции отклика Хо для невзаимодействующих электронов: Р(ш) = Ро(со)^~1)ц. Эффективность метода целиком зависит от размера матрицы = (т^А^^ по сравнению с матрицей А. В наших примерах мы достигаем приемлемой точности для Р(ш) уже с маленькими подпространствами Крылова А^к ~ 10 для всех рассмотренных молекул от метана СН4 (размер матрицы А, Ид = 120) до нафталина СюН8 (размер матрицы А, Ид = 3504).
4.3. Расчет матриц взаимодействия Кн и Кхс
Определим собственную энергию /нхе в базисе произведений собственных состояний Фе(г)Фе(г). Для этого мы должны рассчитать интегралы
КГ-Е'р' = [[ ; (23)
KxcF,E'F' = Jj d3rd3rf фЕ(г)фР(г)/хс(г)6(г - г')фЕ'(г')фРі{г'). (24)
Собственные состояния Фе(г) являются решениями уравнений Кона-Шэма для определенной энергии Е. При этом для расчета собственных состояний мы используем широко известный метод линейной комбинации локализованных атомных орбиталей fa{r)-В результате диагонализации гамильтониана Кона-Шэма собственные состояния получаются в виде
Фе(г) = J2xpfa(r); (25)
а
- 45 -
где Хр - собственные вектора, а /а(г) - атомные орбитали — элементарные функции. обладающие сферической симметрией /а(г) = /га(у>)5'гата(г) — произведение радиальной функции и сферической гармоники3. Подставляя собственные состояния (25) в уравнение (23). мы можем выразить ядро кулоновского взаимодействия К-ц через так называемые четырехцентровые интегралы
(оЬ\с<1) = ЦЛ?'гЛ3г' /„(гЩг) 1 /с(г')1а(г'),
КЕР,Е'Р' = '£x*X[X?Xf(ab\cd). (26)
abed,
Прямой расчет собственной энергии Хартри по формуле (26) чрезвычайно неэффективен. Для эффективного расчета матриц К-ц и Кхс используют либо метод вспомогательных функций (auxiliary functions) [10], либо метод аппроксимации плотностей (density fitting) [5, 11]. Мы используем аналогичный метод доминирующих произведений [16]. Метод доминирующих произведений позволяет разложить произведение локализованных орбиталей fa{r)fb{r) в базисе ортогональных локализованных функций F^{r) — доминирующих произведений
/,МЛ(г) = £^-рт(>-)- (27)
Учитывая это разложение, мы получим выражение для ядра Хартри (23)
КЕР,Е'Р' = Y^(x^\fx[)(xfv;dx^)sSJ.. (28)
abedfiv
fs = Jf*rtуМр. (29)
Обозначая WffF = , мы можем переписать уравнение (28) в матричных
ЕF Е1 Ff т1
обозначениях Кн ’ = Wf-^W . Для собственной энергии обменно-корреляционного
взаимодействия выражения очевидны - они получаются заменой /н —> /хс в последней формуле. При этом, если мы ограничиваемся функционалами vxc(r), зависящими от локального значения плотности п(г). то матрица обменно-корреляционного взаимодействия примет вид
/£- = J d3rF^(r)fxc(r)Ft'(r). (30)
В завершение следует рассчитать интегралы и
Доминирующие произведения конструируются отдельно для каждой пары атомов. Поэтому возможны два вида произведений: локальные (когда атомы в паре совпадают) и билокальные (когда атомы в паре различаются). В случае локальных произведений, доминирующие функции 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,т,(г'),
(31)
(32)
имеет размер Л^осс-^у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(ф*(г,*)ф(г,г)).
Пусть возмущенная функция ф(г^) удовлетворяет уравнению Шредингера
(i^-H0-6V^(r,t)=0. (33)
Будем искать эту функцию в виде ф(г^) = фо(г^) + 5ф(г^), где фо(г^) - волновая
функция певозмущеппого состояния удовлетворяет уравнению
(i^j. - Но^ Фо(г, t) = 0. (34)
Подставляя ф(г^) = фо(г^) +5ф(г^) в уравнение (33). мы получим для поправки 5ф(г^) следующее неоднородное уравнение
5ф(г,г) = 6Уф0(г,г). (35)
Решить последнее уравнение можно стандартным методом функции Грина. Невозму-щеппая функция Грина (зависящая от времени) удовлетворяет уравнению
- H0^j Go(r, t, r', t') = 5(г - r')5(t - t'). (36)
Функция Грина есть ядро интегрального оператора обратного — //о) и может быть использована для решения уравнения (35)
8ф{г,£) = с13г'сИ' Сго(г,^г/,^)5У(г/,^)'0о(г',^). (37)
Записывая далее изменение электронной плотности, мы получим
8п(г,£) = ф*(г^)5ф(г,г)±ф0(г,г)5ф*(г11). (38)
Поэтому функция отклика х(г,£,г',£') будет
Х{т,г,г',1!) = ^ =фп(г^)О0(г,Ь,г'^')ф0(г'^') + с.с. (39)
Выражение для функции Грина через сумму по всем собственным состояниям фоп(г, невозмугценного уравнения Шредингера необходимо в случае не зависящего от времени гамильтониана Но- В этом случае уравнение (36) можно регуляризовать путем добав-
(.д \
ления малой мнимой константы к оператору 1—— Но . Прямой подстановкой
дЬ
можно убедится. ЧТО Со(г, г', ¿') имеет вид [26]
= -1^2ф0п(г^)ф*п(г'^')е±е(-1 (40)
,, . . 0,ж>0 ( л Г1, ж > 0
где 0+(ж) = < 1 „ , 0 (ж) = < „ „ .
4 ' 1 -1,х < 0 ’ к ' 1 0,х < 0
Тогда, подставляя явное выражение для функции Грина (40) в формулу (39). полу-
чим
х(г,г,г',1?) = -[ф*{г,г)^2ф0п(г,ф*п{г',г')е±^ *,)0±(г-г/)^о(»Л*/)+
п
+ 1ф0(г,г)^2ф*п(г^)ф0п{г,^,)е±е^')в±{1-^)ф*(г1^1). (41)
п
Из этого выражения уже явно видно, что функция отклика х(г, г', £') зависит только от разницы 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.
