НАУКА и ОБРАЗОВАНИЕ

УДК.539.192   

МГТУ им. Н.Э. Баумана

vladimirglushkov@yandex.ru

Введение

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

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

В основе метода многочастичных функционалов плотности лежит обобщенная теорема Хоэнберга-Кона, утверждающая, что полная энергия основного состояния многочастичной нерелятивистской  квантовой системы с фиксированным числом частиц N представляет собой однозначный функционал многочастичной функции плотности n_m, m<N, минимум которого реализуется на функции n_m0, соответствующей распределению  частиц в основном состоянии системы [1, 2].

Постановка задачи

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

При исследовании систем, имеющих прикладное значение и состоящих из двух разнородных подсистем, например, молекул или наноразмерных кластеров (в которых роль «подсистемы 1» играют электроны, а роль «подсистемы 2» - атомные ядра или ионные остовы) на практике, как правило, движение частиц одной из подсистем игнорируется (приближение Борна - Оппенгеймера). Иными словами, при изучении подобных систем постановка задачи сводится к исследованию движения частиц одной из подсистем (допустим, «подсистемы 1») в поле, созданном неподвижными частицами «подсистемы 2». При этом исследование проводится исходя из вариационного принципа, сформулированного в теореме Хоэнберга-Кона [1]. При рассмотрении многокомпонентных систем такое приближение может привести к ошибке в определении основных параметров системы, в этом случае используют другие приближения (non-Born-Oppenheimer) [3].

Гамильтониан системы состоящей из двух типов нетождественных ферми-частиц имеет вид:

.                                                  (1)

Здесь

,

,

где  - оператор кинетической энергии i-ой частицы первого (второго) типа,  и  – потенциалы взаимодействия i-ых частиц первого и второго типов с внешним полем,  – взаимодействие i-ой и j-ой частиц первого типа между собой,   – взаимодействие i-ой и j-ой частиц второго типа между собой,  – взаимодействие i-ой частицы первого типа и j-ой частиц второго типа.

Функция плотности двухчастичной системы определяется выражением:

,

где R – вся совокупность координат,  - биномиальные коэффициенты [4].

Для двухкомпонентных систем с гамильтонианом вида (1) справедлива обобщенная теорема Хоэнберга-Кона:

1) ожидаемое значение любой  наблюдаемой  величины, описываемой оператором F, является однозначным функционалом m₁m₂ - частичной функции плотности основного состояния:

;

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

.

В качестве примера, с помощью которого можно оценить разницу в определении основной энергии системы, используя различные приближения, можно взять двухатомные молекулы. Данный выбор связан с тем, что двухатомные молекулы имеют относительно простую структуру: «подсистема 1» - 2 частицы (электроны) и «подсистема 2» - 2 частицы (ядра). Число частиц в данной системе невелико. Также для таких структур имеется большое число экспериментальных данных. Для выполнения поставленной задачи необходимо рассмотреть данную структуру в рамках теории функционала плотности, используя  приближение Борна-Оппенгеймера и приближение локальной плотности.

Техника расчетов

В общем случае изучение любой квантовомеханической задачи сводится к нахождению решения уравнения Шредингера. Для ее решения используют ряд приближений.

Первое из них – приближение Борна-Оппенгеймера. Оператор Гамильтона молекулы с N ядрами и n электронами содержит члены кинетической энергии электронов, потенциальной энергии притяжения электронов к ядрам, а также члены, обуславливающие межэлектронное отталкивание. Кроме того, по сравнению с гамильтонианом атома появляются слагаемые, описывающие электростатическое отталкивание ядер и их кинетические энергии:

.       (2)

Здесь индексы α и β принадлежат атомным ядрам, а индексы i и j относятся к электронам. Введены обозначения ,  и .

Так как гамильтониан молекулы (2) зависит не только от координат электронов, но и от ядерных координат, полная волновая функция системы должна содержать как электронные (r), так и ядерные (R) координаты. Это значительно усложняет задачу математического поиска волновой функции. Поэтому в конкретных расчетах молекулярных свойств стремятся обычно к раздельному рассмотрению движения ядер и электронов.

Вид гамильтониана (2) существенно усложнен по сравнению с гамильтонианом многоэлектронного атома главным образом из-за наличия слагаемого, описывающего кинетическую энергию ядер. Однако масса ядра значительно превышает массу электрона  (масса самого легчайшего ядра протона в 1836 раз больше массы электрона) [5]. Соответственно скорость движения ядер мала по сравнению со скоростью движения электронов. В результате медленно движущиеся ядра образуют электростатическое поле, в котором с намного большей скоростью движутся электроны, успевающие почти мгновенно подстроится к любому изменению координат ядер. Поэтому в первом приближении можно считать ядра атомов фиксированными и рассматривать только движение электронов. На языке квантовой механики такое приближение эквивалентно допущению, что полная волновая функция молекулы Ψ(r,R) может быть выражена в виде произведения электронной Ψ_э(r,R) и ядерной Ψ_я(R) функций:

                                                      (3)

Координаты ядер R входят Ψ_э(r,R) в качестве параметров.

Рассмотрим условия, при которых справедливо допущение (3). Запишем уравнение Шредингера для молекулы с гамильтонианом (2) и волновой функцией (3):

                (4)

где  - энергия отталкивания ядер

 - энергия притяжения электронов к ядрам

 - энергия отталкивания электронов.

Обозначим

                                           (5)

.

Электронная функция Ψ_э(r,R) определяется как собственная функция оператора H_э, т.е.

                                                      (6)

где E_э - электронная энергия, обусловленная движением n электронов в поле N ядер молекулы, плюс энергия взаимодействия между ядрами V_яя. Величину E_э называют адиабатическим электронным термом молекулы или адиабатическим потенциалом.

Учитывая, что

,

,

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

Пренебрегая выражением в первых круглых скобках

                           (7)

и используя соотношения (5) и (6), получим

                                          (8)

Разделив все члены уравнения (8) на Y_э и принимая во внимание (6), получим уравнение для определения Y_я

.                                                     (9)

Условие (7) означает, что электронная волновая функция Ψ_э должна быть настолько медленно меняющейся функцией ядерных координат R, что можно пренебречь ее первой и второй производными по этим координатам. М. Борн и Р. Оппенгеймер впервые показали, что электронные волновые функции обычно подчиняются этому условию с требуемой степенью точности [6].

Приближение (3) является весьма существенным для квантовой химии [7], его называют приближением Борна-Оппенгеймера или простым адиабатическим приближением. В этом приближении полная энергия молекулы представляет собой сумму электронной энергии, вычисленной при фиксированной конфигурации ядер, и колебательно-вращательной энергии ядер:

.                                                                    (10)

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

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

Тогда для того чтобы решить уравнение

                                                                 (11)

с гамильтонианом H (2), представим Y(r,R) в виде ряда

                                                  (12)

где  и  - ортонормированные волновые функции ядер и электронов соответственно в состоянии l для данной конфигурации ядер R. Это означает, что выполняется условие

,

где интегрирование производится по координатам всех электронов.

Подставим функцию (12) в уравнение (11). Учитывая правила дифференцирования, получим уравнение

               (13)

Поскольку оператор H_э не включает операции дифференцирования по R, после умножения предыдущего выражения слева на функцию , интегрирования по электронным координатам r и учета уравнения (13) получим систему уравнений

                              (14)

где оператор  определен как

Оператор  зависит от ядерных координат и определяет так называемое вибронное взаимодействие ядер и электронов.

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

.                                 (15)

Система уравнений (15) отличается от (9) на величину , и приближение, основанное на решении системы уравнений (15) называют адиабатическим в отличие от приближения Борна-Оппенгеймера. С точки зрения теории возмущения, приближение Борна-Оппенгеймера является нулевым приближением к точному решению уравнения (14), учет адиабатической поправки  отвечает учету возмущения первого порядка, а учет недиагональных элементов , отвечающих взаимодействию различных электронных состояний, характеризует поправку более высокого порядка. Применение адиабатического приближения связано с малостью возмущающего вклада . Для случая устойчивых многоатомных молекул существует простой критерий применимости адиабатического приближения [9]

,                                                                    (16)

где ν – наибольшая из частот малых колебаний ядер вблизи точки равновесия,  и  - энергии двух соседних электронных состояний. Критерий (16) выполняется для большинства молекул, однако при сближении поверхностей потенциальной энергии различных электронных состояний молекулярной системы, когда разность между ними становится сравнима с колебательным квантом, приближение Борна-Оппенгеймера становится неудовлетворительным. В области сближения, касания или пересечения поверхности потенциальной энергии происходит смешивание электронных состояний вследствие сильного взаимодействия электронного и ядерного движений. С точки зрения классической механики, в этой области сближения поверхностей потенциальной энергии скорость движения ядер приближается к скорости движения электронов. Квантово-механически это означает, что в областях пересечения или сближения поверхностей потенциальной энергии нельзя пренебрегать оператором кинетической энергии ядер и необходимо решать общее электронно-ядерное уравнение (14), где по крайней мере некоторые из диагональных элементов  отличны от нуля [6, 7].

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

Следующее важное приближение вытекает из того, что реально «свободными» являются только электроны верхнего уровня. Остальные прочно связаны с ядрами своих атомов, что позволяет исключить их из рассмотрения, заменив ядра ионами, заряд которых равен элементарному заряду, помноженному на w– число валентных электронов в атоме.

Введя волновую функцию системы Ψ в виде определителя матрицы  ортонормированных одночастичных волновых функций φ_i и решая задачу об условном экстремуме, можно перейти к системе интегро-дифференциальных уравнений Хартри-Фока [11]:

          (17)

Здесь  - не зависящая от спина одночастичная матрица плотности,  - ее диагональный элемент, иначе называемый функцией плотности (функцией электронной плотности), V_ext(r) – потенциал, создаваемый ионами кристаллической решетки.

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

Основному состоянию соответствует набор волновых функций, доставляющих минимум энергии системы.

С другой стороны, согласно теореме Хоэнберга-Кона, энергия системы полностью определяется функцией электронной плотности . Тогда реальная система может быть заменена совокупностью электронов, представляемых набором пробных волновых функций ψ_i(r), дающих ту же общую плотность, что и исходная система. Такое рассмотрение приводит к уравнениям Кона-Шэма [11]:

                            (18)

Здесь обменно-корреляционный вклад V_XC(r) может быть рассмотрен как возмущение или, если известен вид функционала обменно-корреляционной энергии E_XC[ρ], как

                                            (20)

Эта система уравнений является нелинейной и решается самосогласованно. Расчет существенно упрощается, если перейти к рассмотрению функционалов только электронной плотности n(r). Тогда решение задачи может быть сведено к исследованию функционала полной энергии системы

                  (21)

Здесь

 - энергия электрон-ионного взаимодействия,

 - энергия ион-ионного взаимодействия;

 - энергия электрон-электронного взаимодействия; V(r) – потенциал, создаваемый положительными зарядами с плотностью ρ₊(r).

Если рассматривать электронный газ как однородный, то можно ограничиться  приближением Томаса-Ферми [11]. Функционал кинетической энергии в этом приближении имеет вид

  .                                        (22)

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

Функционал обменной энергии в приближении локальной плотности равен

.                                                           (23)

Он может быть представлен в виде суммы дираковского обменного вклада [4]

  .                   (24)

и корреляционного вклада, который представляют в форме Вигнера [11]

 .                                       (25)

Результаты расчетов

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

,                            (26)

где  μ, ν - эллипсоидальные координаты, γ - вариационный параметр.

Общий вид функционала полной энергии можно построить в виде

 .                                            (27)

Для функционала T[n₂] применяется метод обобщенного градиентного приближения, хорошо зарекомендовавшее себя в одночастичном методе функционалов плотности [13].

Функционал взаимодействия электронов с ядрами V[n₂]:

,           (28)

где .

Функционал, описывающий взаимодействие электронов между собой W[n₂], в отличие от одночастичных подходов, может быть точно представлен в виде функционала:

            .                     (29)

Для представления данного функционала в эллиптических координатах использовалось разложение Неймана [14].

В рамках предложенного подхода можно получить выражение для E[n₂], минимизация которого позволит найти верхнюю оценку для энергии основного состояния системы, а также наилучшее приближение к функции , соответствующей пространственному распределению частиц в основном состоянии системы.

Равновесное значение расстояния между ядрами R₀=1,50 (данные представлены в атомной системе единиц). Рассчитав равновесное состояние, то есть, получив значение вариационного параметра, можно усовершенствовать двухчастичную функцию плотности

,             (30)

где α₁, α₂, α₃ - вариационные параметры.

Функция (30) также неотрицательна и обладает симметрией относительно перестановок .

Параметры функции плотности α₁, α₂, α₃ и связанные с ним характеристики системы могут быть определены вариационными методами по аналогии с предыдущим вычислением.

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

Таблица равновесных значений параметров.

Распределение электронной плотности вдоль прямой соединяющей ядра атомов представлено на рис. 1.

Для данной структуры имеем следущие результаты: равновесное расстояние между атомами , энергия ионизации  (результаты расчетов представлены в атомной системе единиц).

Рис. 1 Распределение электронной плотности вдоль прямой соединяющей ядра атомов

Зависимость полной энергии системы от межъядерного расстояния представлена на рис. 2

Рис. 2. Зависимость полной энергии от межъядерного расстояния

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

.                                       (31)

Если рассматривать электронный газ как однородный, то можно ограничиться  приближением Томаса-Ферми [11]. Функционал кинетической энергии в этом приближении имеет вид:

 .                        (32)

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

Проверка полученных результатов, проводилась методом функционала плотности с использованием псевдопотенциалов и разложения волновых функций по плоским волнам, которые реализованы в программе ABINIT[10, 15]. Обменно-корреляционное взаимодействие описывалось в приближении локальной плотности (LDA). В качестве псевдопотенциалов использовались оптимизированные сепарабельные нелокальные псевдопотенциалы [8].

Параметры системы с гамильтонианом вида (1), полученные в результате численного расчета, с небольшой погрешностью (порядка 0,05 % относительно энергии ионизации) совпадают с численно-аналитическими результатами, что позволяет сделать вывод о применимости данного подхода в описании двухкомпонентных систем.

Заключение

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

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

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

Литература

1. HohenbergP.,  Kohn  W.  Inhomogeneous  electron  gas. Phys.Rev. - 1964. - V.136: B864-B871.

2.                             Levy M., Freed K. Direct first principles algorithm for the electron density functional. J. Chem. Phys. – 1982. – 77: 396-398.

3. Nalewajski R.F., Capitani J.F. Density functional theory: non Born-Oppenheimer Legendre transformations and Maxwell relations, equilibrium and stability conditions. J. Chem. Phys. – 1982. – 77: 2514-2526.
4. Еркович О.С., Метод многочастичных функционалов плотности в нерелятивистской квантовой механике. – М.: МГТУ им. Баумана, 2001. –100 с.
5. Физическая энциклопедия, т. 5. Стробоскопические приборы - Яркость / Гл. ред. А.М. Прохоров. Ред. кол.: Д.М. Алексеев, А.М. Балдин, А.М. Бонч-Бруевич и др. - М.: Большая Российская энциклопедия, 1998. - 760 c.: ил.
6. Минкин В.И., Симкин Б.Я., Миняев Р.М. Теория строения молекул. - Ростов-на-Дону: Изд-во «Феникс», 1997 г. - 560 с.
7. Спектры и строение двухатомных молекул. / Под ред. В.Н. Кондратьева. - М.: Издательство иностранной литературы, 1949. - 420 c.
8. J.P. Perdew, A. Zunger. Phys. Rev. B 23, 5048 (1981).
9. Киржниц Д.А. Полевые методы теории многих частиц. - М.: Госатомиздат, 1963. - 344 с.Давыдов А.С. Квантовая механика.- М.: Наука. - 1973.- 704 с.
10. Richard M. Martin. Electronic Structure. Basic Theory and Practical Methods. Cambridge University Press. 2004. 624 р.
11. Теория неоднородного электронного газа. / Под ред. С. Лундквиста и Н. Марча. - М.: Мир, 1987. - 400 c.
12. Parr R.G., Yang W. Density-functional Theory of Atoms and Molecules. New York: Oxford University Press. 1989. 333 р.

13.  Dreizer R.M., Gross E.K.U. Density Functional Theory. Berlin: Springer-Verlag. 1990. 303 p.

14.  Соколов А.А., Лоскутов Ю.М., Тернов И.М. Квантовая механика. – М.: Просвещение, 1965. – 640 с.

15. X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. Ghosez, J.-Y. Raty, D.C. Allan. Comput. Mater. Sci. 25, 478 (2002).