CC BY
301
УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА
Том 151, кн. 3
Физико-математические пауки
2009
УДК 539.21^544.255
АВ INITIO РАСЧЕТЫ СТРУКТУРНЫХ
И ЭЛЕКТРОННЫХ СВОЙСТВ КРИСТАЛЛИЧЕСКИХ ТВЕРДЫХ ТЕЛ
В ПРИБЛИЖЕНИИ ФУНКЦИОНАЛА ПЛОТНОСТИ И ПСЕВДОПОТЕНЦИАЛА В ИМПУЛЬСНОМ ПРОСТРАНСТВЕ: ДЕТАЛИ И ПРИМЕРЫ
В. В. Клековкипа, P.M. Амипова
Аннотация
В статье обсуждаются методы теории функционала плотности, теории псевдопотеп-цпала и базиса плоских воли, которые составляют основу ab initio изучения строения и физико-химических свойств кристаллов. Приводятся результаты расчетов структурных и электронных свойств для кристаллов углерода, кремния, германия и кристалла никеля, которые сравниваются с экспериментальными данными.
Ключевые слова: теория функционала плотности, приближение локальной плотности. метод псевдопотепциала, базис плоских волп.
Введение
Неэмпирические методы изучения строения и разнообразных свойств вещества начали развиваться с конца 30-х годов двадцатого столетия, сразу после появления нерелятивистского волнового уравнения Шредингера НФ = ЕФ. Развитие вычислительных методов шло по двум направлениям: в направлении разработки методов Хартри Фока Рутана. с одной стороны, и методов, основанных на теории функционала плотности, с другой.
Методы Хартри Фока Рутана заключаются в решении соответствующих дифференциальных уравнений методом самосогласованного поля в приближении независимого электрона. При этом каждый отдельный электрон рассматривается в некотором усредненном эффективном потенциальном поле, созданном ядрами и остальными электронами. Однако возрастающие требования к точности расчетов стимулировали поиск более корректного решения многоэлектронной задачи с учетом эффектов электронной корреляции. Рост компьютерных возможностей и разработка более совершенных математических алгоритмов способствовали тому, что методы Хартри Фока Рутана в сочетании с методами учета электронной корреляции достигли уровня предсказательной силы при изучении физико-химических свойств молекул в газовой фазе и растворах.
По сравнению с методами Хартри Фока Рутана в рамках теории функционала плотности проблема учета электронных корреляций решается более просто. Варьируемой величиной в этой теории является электронная плотность функция. зависящая от трех пространственных координат (без учета спина электрона). Методы теории функционала плотности получили широкое развитие благодаря возможности простой и наглядной интерпретации результатов вычислений, в отличие от методов, основанных на формализме многоэлектронных волновых функций. Кроме того, немаловажным, на наш взгляд, оказывается тот факт, что проведение расчетов в рамках методов теории функционала плотности проще, чем в рамках
других методов учета электронных корреляций, и точность расчетов подобного рода сравнима с экспериментальной точностью определения соответствующих величии.
Развитие вычислительных методов применительно к кристаллическим материалам произошло несколько позднее [1]. чем для молекулярных систем. Основополагающую роль при этом сыграло появление теории псевдопотенциала и поиск эффективных форм для него, применение которых допускает проведение расчетов с использованием базисного набора, состоящего из сравнительно небольшого числа плоских волн.
Широкому распространению вычислительных методов в квантовой теории твердого тела, безусловно, способствовал недавний значительный прогресс в увеличении доступной компьютерной мощности, что позволило многим исследователям проводить расчеты для сложных материалов, изучение которых теоретическими методами еще 10 20 лет назад представлялось совершенно невозможным. Следует подчеркнуть, что фундаментом для решения в настоящее время широким кругом исследователей целого спектра различных задач стали результаты упорного труда многих известных теоретиков, разработавших эффективные методы и алгоритмы расчетов «из первых принципов» для изучения свойств материи (атомов, молекул и твердых тел) и усилий программистов, оформивших эти методики в виде доступных компьютерных программ.
Точность расчетов физико-химических свойств кристаллов достигла такого уровня, что сегодня является общепризнанным тот факт, что методы компьютерного моделирования могут быть использованы и используются для дизайна и разработки материалов со специфическими свойствами и для предсказания поведения твердых тел в экстремальных условиях при высоких давлении и температуре (некоторые недавние работы [2 6]).
В настоящей статье нами детально обсуждается применение методов теории функционала плотности для описания электрон-электронных взаимодействий и теории псевдопотенциала для описания электрон-ионных взаимодействий в твердом теле, что в силу естественных причин, как правило, опускается в оригинальных работах, представляющих новые научные результаты. Заметим, что ясное понимание этих ключевых концепций необходимо при проведении расчетов параметров кристаллов с использованием компьютерных программ. Дальнейшие подробности можно почерпнуть из имеющегося множества обзоров и монографий, приведенных в списке литературы (например. [7 11]).
В разделе 1 нами рассмотрены основы метода теории функционала плотности, описываются оригинальный вариационный принцип Хоэнберга Кона, самосогласованные уравнения Кона Шэма. приближение локальной плотности.
В разделе 2 обсуждается формализм импульсного пространства для периодических систем, приводится вид уравнений Кона Шэма при использовании базиса плоских волн.
Метод псевдопотенциала и его современные формы подробно обсуждаются разделе 3.
Результаты расчетов для простых кристаллических систем (углерода, кремния, германия и никеля), проведенных нами в рамках изложенных в предыдущих разделах подходов, описаны в разделе 5.
Итоги и выводы приведены в заключении.
1. Теория функционала плотности
В этом разделе излагаются основные положения теории функционала плотности (Density Functional Theory, DFT), которая широко применяется при изучении
электронной структуры различных материалов и приобрела огромную популярность в последнее время, поскольку ИП предоставляет достаточно надежный метод для описания эффектов обмена и корреляции в электронном газе.
Освещению вопросов, касающихся применения БЕТ, посвящено много обзоров и книг (например. [13 21]).
1.1. Теорема и вариационный принцип Хоэнберга —Кона. В основе методов БЕТ лежит фундаментальная теорема, доказанная П. Хоэнбергом и В. Коном [22] в 1964 г.1 Согласно этой теореме плотность основного состояния2 п(г) связанной системы взаимодействующих электронов в некотором внешнем потенциале V(г) однозначно (вернее, с точностью до не представляющей интереса аддитивной константы) определяет этот потенциал. В случае вырожденного основного состояния это утверждение справедливо для плотности п(г) любого основного состояния. Здесь мы не приводим доказательство этой теоремы: для невырожденного основного состояния оно довольно простое и его можно найти в любом обзоре, посвященном БРТ. Для вырожденного случая доказательство приведено в книге [23, р. 4].
Поскольку плотность п(г), очевидно, определяет число частиц N = J п(г) ¿г
(здесь для простоты мы рассматриваем систему из конечного числа электронов) и согласно теореме Хоэнберга-Кона потенциал V (г), то данной плотности, таким образом, можно сопоставить полный гамильтониан Н электронной системы, нерелятивистская форма которого в рамках адиабатического приближения имеет
N
Н = Т + и + £ V (гО, (1)
¿=1
1 N N 1
где Т = - - ^^ V2 - оператор кинетической энергии электронов, и = --- -
2 ¿=1 ¿<з |г¿ 1
энергия электростатического отталкивания электронов, V(г) - внешний, в рассматриваемом сейчас случае, скалярный потенциал (включая потенциал взаимодействия с ядрами), ^ - координата г-го электрона, N - полное число электронов в системе. В настоящей работе при записи формул нами используются атомные единицы Хартри (К =1, те = 1, е =1). п( )
Н
уравнения Шредннгера (даже при наличии дополнительных возмущений, таких, как электромагнитное поле). Таким образом, среднее значение любой физической величины может быть выражено в виде функционала электронной плотности.
В терминах электронной плотности среднее значение гамильтониана (1) в основном состоянии Фо имеет следующий вид:
N
(Фо|Н |Фо> = (Фо|Т + и |Фо> + (*о^ (^)|*о> =
¿=1
= Е[по(г)] = ^[по(г)] + У по(т^(г) <к, (2)
1Вклад В. Кона в «развитие теории функционала плотности» был отмечен Нобелевской премией по химии в 1998 г.
2 Вслед за Хоэнбергом и Коном мы используем обозначение га(г) для электронной плотности, что является общепринятым. Заметим, что в научной литературе довольно часто также встречается обозначение р(г). Ради краткости мы будем говорить просто плотность, подразумевая полную зарядовую электронную плотность, если но оговорено особо.
где no(r) = N j |Ф0(г, г2,... , гN)|2 dr2 ... drN - электронная одночастичная плотность основного состояния (через j dr здесь обозначено интегрирование по пространственным и суммирование по спиновым переменным). В уравнении (2) F[«(г)] обозначает универсальный функциональный оператор плотности, который будет определен ниже. Форма этого функционала одинакова для всех систем, однако точный вид этого функционала неизвестен, поскольку невозможно найти точное решение многоэлектронной задачи.
В работе [22] был впервые сформулирован принцип минимума полной энергии но через пробные волновые функции (как в вариационном принципе Релоя Ритца), а через пробные плотности «(г). Более простой и общий вывод вариационного принципа был предложен М. Лови [24] и Е. Либом [25. р. 111]. Строго говоря, класс физически реализуемых пробных плотностей «(г) ограничен. Вариационный принцип Хоэнборга Кона был сформулирован для так называемых V-иредставимых плотностей, то есть таких плотностей, которые могут быть ре-
N
продставимыо плотности, то есть плотности, которые соответствуют некоторым антисимметричным волновым функциям. Пространство V-представимых плотно-
N
метим только, что изучение топологии областей V-представимости в абстрактном пространстве всех возможных функций плотности «(г) продолжается и в настоящее время, однако нерешенность этой задачи, как будет показано ниже, не является препятствием для практических применений DFT.
V
представимой пробной плотности «(г) справедливо следующее неравенство:
E[«o(r)] < E[«(г)]. (3)
Формально этот принцип сводит задачу минимизации энергии по многочастичной волновой функции (как в вариационном принципе Релея Ритца) к кажущейся тривиальной задаче минимизации по плотности (функции трех переменных), которая должна решаться методом неопределенных множителей Лагранжа с дополнительным условием j «(г) dr = N. Однако остается открытым вопрос о выборе V
Эта проблема, как отмечено выше, была решена Лови [24] и Либом [25], которые предложили так называемый метод условного поиска. В этом методе минимизация электронной энергии проводится в два этапа. Во-первых, при фиксированной пробной плотности «(г) проводится поиск условного минимума энергии на классе пробных функций ^(г)' поРожДа1°Щих данную плотность. На втором этапе проводится минимизация по плотности:
E = min
п (г)
min (Vj?(P)|T + U|*?(Г)> + / «(r)V(r) dr)
F[«(r)] + I «(r)V(r) dr
min
n(r)
(4)
где функционал ^[п(г)], который фигурировал в уравнении (2), имеет вид
F[«(r)] = min <^(г)1 T + U|Ф0?^)>.
При таком подходе требование V-представимости заменяется на достаточно мягкое требование Ж-представимости, однако поиск Г[п(г)] приводит снова к трудной задаче минимизации по 3Ж-мерным пробным функциям.
Приведенные выражения формально являются совершенно строгими, хотя явный вид функционала Г[п(г)] неизвестен.
Отметим, что частным случаем БЕТ является теория Томаса Ферми, в которой для кинетической энергии Т используется следующее приближение:
Т = | п(г)10к2(п(г)) ^ (5)
где к? (п) = [3п2п]1/3 фермиевский волновой вектор для однородного электронного газа с плотностью п, (3/10)к|(п) - средняя кинетическая энергия, приходящаяся па один электрон в таком газе. Энергия взаимодействия электронов V записывается в классическом приближении:
1 ¡' п(г)п(г') , , ,
7 =2/ Т^ГЛтЛт• (6)
Обменная поправка была введена в теорию Томаса Ферми Дираком:
К =3( 3)1/3/п 3 (г)йг. (7)
Известно, что погрешность теории Томаса Ферми обусловлена в основном
Т
тельной степени был исправлен в самосогласованных одночастичных уравнениях Кона Шэма.
1.2. Самосогласованные уравнения Кона —Шэма. В. Кон и Л. Шэм в
работе [29] вывели одночастичные уравнения, путем решения которых можно определить электронную плотность, доставляющую минимум функционалу энергии. Кратко опишем вывод этих уравнений, который, как нам представляется, является весьма оригинальным. Фактически. Кон и Шэм получили уравнения Хартри из вариационного принципа Хоэнберга Кона для энергии (уравнение (3)). который формально является точным и потому содержит в себе уравнения Хартри вместе с уточнениями этих уравнений.
Кон и Шэм рассмотрели систему невзаимодействующих электронов, записав для нее вариационный принцип Хоэнберга Кона и соответствующее уравнение Эйлера Лагранжа. Для системы невзаимодействующих электронов одночастичные уравнения, из которых может быть определена энергия основного состояния, хорошо известны. Поэтому, эти авторы явным образом выделили кинетическую энергию невзаимодействующих электронов и классическую часть кулоновской энергии взаимодействия электронов (ур. (5), (6)), записав функционал Г[п(г)] в следующем виде:
Г[п(г)] = Т Ж] + 7[п(г)] + Ее Ж], (8)
Ее Ж] = Т[п(г)] - Т Ж] + V [п(г)] - 7[п(г)]. (9)
Функционал Ее [п(г)] называется обменно-корреляционным функционалом, и выражение (9) является для него определением. Ее [п(г)] содержит вклады в энергию, которые трудны для расчета, а именно: вклад, обусловленный фермиевской
и кулоновской корреляциями, поправку на самодействие, вносимое классическим кулоновским потенциалом.
Затем Кон и Шэм записали вариационный принцип Хоэнберга Кона и уравнение Эйлера Лагранжа уже для системы взаимодействующих электронов с учетом представления (8). Оказалось, что уравнение Эйлера Лагранжа для взаимодействующих электронов по своему виду идентично уравнению Эйлера Лагранжа для невзаимодействующих электронов, которые находятся во внешнем эффективном потенциале V®:
Ъ "> = V « + / ^ *' + ё-Ш- <10>
Поэтому движение электронов во внешнем потенциале V(г) с учетом их взаимодействия должно описываться теми же одночастичными уравнениями, которые справедливы в модели невзаимодействующих электронов, находящихся в потенциале V®:
1 v2 +Veg(г)) ^(г) = £^(г), г = 1,...,Ж. (И)
Здесь фг(т) обозначает одноэлектронные волновые функцииКона-Шэма, ег - собственные значения3 уравнения (11). Эффективный потенциал (г) определяется распределением зарядовой плотности согласно (10) и является суммой внешнего потенциала (с учетом потенциала ядер), кулоновской (или. точнее, хартриевской) части электрон-электронного взаимодействия и обменно-корреляционной части, которая равна функциональной производной обменно-корреляционного функционала Ехс[п] по электронной плотности п(г).
Электронная плотность вычисляется непосредственно из орбиталей Кона Шэма:
осе
п(г)=£ Шг)|2, (12)
г
индекс г пробегает только заполненные состояния. Уравнения (11), (12) образуют систему уравнений Кона Шэма. которая решается путем самосогласования.
Заметим, что эти уравнения по форме аналогичны уравнениям метода Хартри, но они много проще. Их решение также оказывается не намного сложнее уравнений Хартри Фока, в которых фигурирует нелокальный потенциал. Кроме того, уравнения Кона Шэма являются точными и при известном функционале обменно-корреляционной энергии позволяют получить точную энергию системы, тогда как уравнения Хартри Фока являются приближенными.
Полная энергия основного состояния вычисляется по следующей формуле:
Е[п(г)] = ±ег - 1 Ц **' - / п(0 * + Ее[п(г)]. (13)
До сих пор речь шла о системах, для которых п-|-(г) = п^(г) = п(г)/2 (п^ и п - плотности электронов с проекциями спина +1/2 и -1/2 соответственно). Очевидно, что в общем случае функции распределения плотности электронов различны, поскольку состояниям электронов с разными значениями проекции спина соответствуют разные наборы пространственных орбиталей и функционалы будут зависеть от обеих функций распределения электронных плотностей. Поэтому теория ИП легко обобщается на случаи магнитных систем или систем, находящихся в магнитном поле [27].
З3аметим, что £1 не являются собственными значениями энергии, а, как оказалось, несут смысл производных энергии но числам заполнения [26].
Следует отметить, что при вычислении энергии электронной системы можно использовать как уравнения Кона Шэма. так и вариационный прицип Хоэнберга Кона. Однако практическое использование уравнений Кона Шэма оказывается много проще.
При использовании изложенной теории существенное значение играет конкретный вид обменно-корреляционного вклада Ec [n], для которого был предложен целый ряд приближений [28]. Одним из первых приближений для Ec [n] является приближение локальной плотности. Забегая вперед, отметим, что это приближение часто используется при проведении расчетов полной энергии кристаллов в рамках метода псевдопотенциала.
1.3. Приближение локальной плотности. Приближения для функцио-
E [n]
физической природы электронной структуры конкретных атомных систем. Наиболее важные приближения имеют квазилокальный вид:
Ec[n(r)] = J ехс(г; [n(r)])n(r) dr,
где ехс(г; [n(r)]) - обменно-корреляционная энергия, приходящаяся па одну частицу в точке г и представляющая собой функционал от распределения плотности n(r) в точках г, близких4 к точке г.
В приближении локальной плотности (Local Density Approximation, LDA) [29] предполагается, что обменно-корреляционная энергия ехс, приходящаяся на один электрон в точке г в электронном газе, равна обменно-корреляционной энергии, приходящейся на один электрон в однородном электронном газе > который имеет такую же плотность, что и рассматриваемый электронный газ в данной точке г. Таким образом, в этом приближении ехс является функцией плотности.
Обменно-корреляционная энергия в приближении LDA изначально была разложена на обменное и корреляционное слагаемые5:
ехс(r; n(r)) = ех(r; n(r)) + ес(r; n(r)).
Плотность обменной энергии ех (r; n(r)) в приближении LDA находится с помощью модели идеального вырожденного электронного газа в потенциальном ящике и была получена, как отмечалось выше, Дираком в качестве поправки к теории Томаса Ферми (уравнение (7)):
3 / 3 \ 1/3 ех(r; n(r)) = - Д-J n1/3(r).
Точная форма ес(r; n(г)) неизвестна. Корреляционный вклад в энергию однородного газа вычислялся многими авторами (см., например, [29, 30]), их результаты согласуются между собой. Среди часто используемых приближений можно назвать формулу Вигнера:
( < W О.44 £с(r; n(r)) = - ^ГТв,
r
цем и определяемый условием (4n/3)r3 = 1/n. Наиболее точные значения функции
4Близость точек г и г означает, что |г — г | - величина порядка локальной фермиевской длины волны.
5 Это разделение довольно произвольно, однако оно оказывается удобным в этом приближении, поскольку известно точное выражение для ех(r; n(r)).
ес(n(r)) были получены Сиперли [31] путем численного моделирования методом Монте-Карло.
Пренебрегая в рамках приближения LDA корреляционной составляющей энергии ес, можно получить точные уравнения для однородного электронного газа в предположении отсутствия кулоновской корреляции. Отметим, что схожее построение было сформулировано Слэтером [32] еще до появления теории DFT: оно известно как метод Xa, или метод Хартри-Фока-Слэтера.
Эффективность LDA была доказана posteriori после удовлетворительных исследований многих материалов и систем6. Оказалось, что полная энергия может быть вычислена в рамках LDA с точностью до нескольких процентов, геометрнческне параметры с точностью до 0.1 А. энергия связи, однако, может иметь ошибку до 10%. Приближения, которые, например, учитывают поправку к обменно-корреляцпонной энергии, обусловленную неоднородностью электронной плотности (обобщенное градиентное приближение Generalized Gradient Approximation). и другие приближения, как правило, несущественно улучшают результаты, полученные в рамках простого приближения LDA. Однако для многих систем учет этих поправок оказывается крайне необходим. В связи с этим поиск эффективных форм функционалов является актуальной задачей и в настоящее время [28].
Таким образом, в теории функционала плотности осуществляется, в принципе. точное отображение задачи об электронном газе с сильными взаимодействиями между электронами на задачу о движении одной частицы в эффективном потенциале. Хотя обмеиио-корреляциоииый потенциал точно неизвестен, для него имеются эффективные приближения, которые позволяют воспроизводить свойства системы в основном электронном состоянии с точностью до нескольких процентов по сравнению с экспериментом.
2. Формализм импульсного пространства для периодических систем
Хорошо известно, что при рассмотрении физико-химических свойств периодических систем немаловажное значение имеет учет далыгодействующих взаимодействий. При рассмотрении возможности расчета энергии периодических систем (твердых тел) в рамках кваитовохимических методов расчета энергии систем. состоящих из конечного числа атомов, возникают две основные трудности. Во-первых, волновая функция должна быть вычислена для огромного числа электронов в системе. Во-вторых, поскольку электронные волновые функции отличны от нуля во всем твердом теле, для их разложения требуется, вообще говоря, базисный набор, включающий атомные функции, центрированные на всех атомах, составляющих твердое тело. Эти проблемы, кажущиеся непреодолимыми на первый взгляд, решаются в рамках следующих подходов. В рамках кластерного подхода твердое тело моделируется некоторым кластером, состоящим из конечного числа атомов [35 37]. Однако такой подход оказывается эффективным, скорее, при изучении апериодических систем, в частности дефектов и примесей, нарушающих трансляционную симметрию кристаллов. В настоящее время при изучении кристаллов. как правило, применяется стандартный подход в теории твердого тела расчеты выполняются с учетом периодических условий и с применением теоремы Блоха для электронных волновых функций.
6Этот факт объясняется случайным совпадением, заключающимся в том, что LDA удовлетворяет правилу сумм для обмешю-корреляциошюй дырки, которая характеризует перераспределение плотности заряда вблизи электрона (перенормировку взаимодействия) [33, 34].
2.1. Базис плоских волн. Согласно теореме Блоха [38] волновая функция электрона в кристалле может быть записана в виде:
^nk = ^ cnk(G) exp[i( k + G)r], (14)
G
где G - вектор обратной решетки, к - волновой вектор, n - индекс зоны, суммирование проводится по всем векторам обратной решетки. Таким образом, электронная волновая функция может быть разложена по дискретному бесконечному базисному набору, составленному из плоских волн.
При проведении расчетов базис плоских волн всегда ограничивают. В базисный набор включаются только такие плоские волны, которые описывают состояние электрона, кинетическая энергия 1/2|k + G |2 которого меньше некоторого значения - Ecutoff (cutoff energy). Хотя такой базис, очевидно, не является полным, заметную величину, как правило, имеют только коэффициенты Cn(k + G) для плоских воли, соответствующих небольшим значениям кинетической энергии. Ошибка, вносимая в расчет полной энергии из-за ограничения базисного набора, всегда может быть уменьшена путем расширения базисного набора за счет добавления в него функций, соответствующих более коротким длинам волн. Одной из сложностей, связанных с использованием базиса плоских волн, является то, что число базисных функций не меняется непрерывно с изменением значения Ecutoff и требуется разное число плоских волн для разных точек к-пространства.
На практике при проведении расчетов приходится выбирать оптимальное значение Ecutoff для каждого вещества, при этом критерием, как правило, служит компромисс между сходимостью результатов для полной энергии и сложностью вычислений.
Заметим, что при изучении твердых тел кроме базиса, включающего только плоские волны, широко используются так называемые смешанные базисные наборы, включающие также орбитали, локализованные на атомах, что оказывается весьма эффективным при рассмотрении, например, систем, содержащих d-и /-элементы [39, 40]. Довольно часто изучение свойств кристаллов проводится с использованием полностью локализованного базисного набора, состоящего из атомных орбиталей [41, 42].
2.2. Специальные точки для интегрирования по зоне Бриллюэна. Решение ряда задач в теории твердого тела (самосогласованные расчеты энергетического спектра и электронной плотности в кристалле, определение полной энергии и др.) связано с суммированием по состояниям с различными значениями волнового вектора к, изменяющегося в зоне Бриллюэна (Brillouin Zone, BZ). В практических расчетах при макроскопических размерах рассматриваемой циклической системы (при введении основной области кристалла) суммирование заменяется интегрированием по BZ:
(k)=(2^/P 00
к BZ
Здесь П - объем элементарной ячейки, P(к) - некоторая, как правило, достаточно гладкая и периодическая функция в к-пространстве с периодом обратной решетки, инвариантная относительно операций точечной группы симметрии кристалла G (например, одноэлектронные блоховскне функции кристалла при определенном выборе фазовых множителей, квадраты их модулей, матричные элементы различных операторов па блоховских функциях, то есть многие величины, которые
фигурируют в теории твердого тела):
р(к) = р(к + ь<) = р(дк), д е С.
Через Ь^ обозначены векторы основных трансляций обратной решетки, группа С
является группой симметрии ВЕ (совпадает с точечной группой симметрии прямой
д
Хотя условия периодичности и инвариантности относительно групповых преоб-
р( )
чения в точках на поверхности ВИ, отстоящих друг от друга на векторы обратной решетки, и для векторов соответствующей звезды, тем не менее для вычисления таких интегралов, в принципе, требуется проводить расчет волновых функций для достаточно плотной сетки к-точек и экстраполировать их значения для промежуточных волновых векторов к. Это, разумеется, требует значительных вычислительных затрат, поскольку расчет значений волновых функций при каждом значении к является трудоемкой задачей.
Однако существуют методы приближенного вычисления интегралов по ВЕ, которые требуют расчета значений этой функции для небольшого набора специальных к-точек. Формула приближенного численного интегрирования имеет следующий вид:
р(к> л = £ (к.),
где - вес, с которым специальная точка кд входит в данный набор
^^^ Шк. = , индекс в нумерует специальные точки.
Методы построения специальных точек описаны, например, в книге [43]. Эффективные наборы специальных к-точек приведены в работах [44 46].
Используя набор из небольшого числа специальных к-точек, можно получить достаточно точное приближение для полной энергии диэлектриков или полупроводников. Для металлов требуется более плотный набор к-точек для точного определения поверхности Ферми. В этом случае необходимо использовать однородную сетку к-точек.
Величина ошибки, вносимой в полную энергию за счет использования редкой сетки к-точек, всегда может быть уменьшена за счет использования более плотного набора к-точек, что, разумеется, увеличивает объем расчетов.
2.3. Уравнения Кона —Шэма в к-пространстве. При использовании базиса плоских воли для разложения электронной волновой функции (см. уравнение (14)), уравнения Кона Шэма (11) принимают вид уравнений для Фурье-компонент Сп(к + в) волновой функции:
Е
(к + С)2£сс ' + Уе,г(к + в , к + С')
Спк (С') = £п(к)Спк (С),
(15)
где Еп(к) — энергия п-й зоны для приведенного к первой ВЕ волнового вектора к, ¿оо' _ символ Кронекера. Через У®(к + в, к + в') обозначены Фурье-компоненты эффективного потенциала
Уе№(к + в, к + в') =
1
N ^
ехр[—¿(к + С)г]Уей- (г) ехр[—¿(к + в >] а г, (16)
где N - число элементарных ячеек в кристалле.
Рис. 1. Блок-схема, описывающая вычислительную процедуру расчета полной энергии
Выражение (12) для плотности в Фурье-представлении имеет следующий вид:
п(г)=(2Л)з/ П(г) (17)
BZ
occ
nk(г) = £ ^k|2. (18)
n
Решение выписанных уравнений проводится как обычно путем процедуры самосогласования, примерная блок-схема которой показана на рис. 1.
В этой процедуре выбирается некоторая пробная плотность n(0), использую которую, рассчитываются потенциал Хартри и обменно-корреляционный потенциал Vic0) (второе и третье слагаемые в уравнении (10)) на нулевом шаге. Затем решается основное уравнение (15), например, путем диагонализации матрицы Гамильтона, матричные элементы которой ffk+G ,k+G' определяются выражением в квадратных скобках в уравнение (15). На основании найденных коэффициентов Cnk(G) вычисляется новая плотность n(1). Далее, уже исходя из найденной новой плотности n(1), вычисляются V^^ и V(c) на первом шаге самосогласования и снова решается уравнение (15). Эта процедура повторяется до тех пор, пока потенциал экранирования, полученный на некотором шаге, не окажется равным с некоторой конечной точностью потенциалу, вычисленному на предыдущем шаге.
Интересно заметить, что размер матрицы Hk+qk+q' зависит от выбора значения Ecutoff и будет настолько большим для системы, которая содержит как валентные электроны, так и электроны остова (core), что осуществление диагонализации этой матрицы станет практически невозможным7. Однако при использовании ме-
7 Для точного описания сильно локализованных волновых функций электронов остова и быстрых осцилляции волновых функции валентных электронов вблизи ядер требуется колоссально
тода псевдопотенциала, описанного ниже, оказывается возможным свести эту задачу к решению сокулярного уравнения приемлемой размерности.
Отметим, что время расчета при использовании обычных методов диагонали-зации матриц возрастает как третья степень размера базиса, поэтому на практике применяются более совершенные методы минимизации функционала энергии, которые требуют меньше вычислительного времени и оперативной памяти. В частности. одним из наиболее широко применяемых в настоящее время методов прямой минимизации функционала Кона Шэма является метод сопряженных градиентов (сощидаЛе-дпмЫетйя). который подробно описан, например, в [47 51].
2.4. Симметризованная зарядовая плотность. Учет не только трансляционной симметрии, но также и точечной группы симметрии кристалла позволяет значительно упростить проведение расчетов. Поскольку одной из основных величин в формализме функционала плотности является электронная плотность п(г), приведем явное выражение для плотности, которое применяется в расчетах при использовании специальных к-точек с учетом пространственной симметрии кристалла. Использование специальных к-точек, как отмечалось выше, позволяет заменить интегрирование по ВЕ в уравнении (17) суммой по нескольким специальным точкам к8 с соответствующими весовыми множителями и>кя:
п(г) = Щ ^ (пк. (г))с . к3
Здесь (пкя (г))с - симметризованная зарядовая плотность в точке к8:
(пК ^У)а = пзк-(г)' д
С .......... точечная группа симметрии кристалла, — ^^ _ операция усреднения по груп-
д
пе. Фурье-компоненты зарядовой плотности вычисляются по следующей формуле:
осе
п(С) = £ - ]Т ехр[г^ • ^ Е /пк Е (С' + ^ )СПк. (С'),
к3 д п в'
где Ь(д) - трансляция, связанная с пространственной групповой операцией симметрии д, /пк - соответствующие числа заполнения, звездочкой обозначено комплексное сопряжение. Соотношение, эквивалентное этому, было впервые дано в работе [52], которая, как отмечено в [11], содержит некоторые опечатки.
2.5. Вычисление полной энергии. В представлении плоских волн кинетическая энергия имеет диагональный вид и вычисляется по формуле
осе
т = ]Т ^]Т/пк]Г |с + к |2Сп(с + к)|2.
п
В уравнении (10), как отмечалось, V(г) обозначает внешний потенциал, в котором находится рассматриваемая система электронов. Будем теперь обозначать
большое число плоских воли. Легко убедиться в том, что, например, для описания волновых функций электронов атома алюминия потребуется ~ 106 плоских волн!
его через Цоп (г), указывая тем самым, что этот потенциал создан ядрами (ионными остовами в рассматриваемом ниже методе псевдопотенциала). Энергия взаимодействия электронов с потенциалом (г — И,/), где II/ - радиусы-векторы,
I
характеризующие положения ядер атомов, вычисляется по формуле
= Н Е п(—с (с)я (с),
н с
где Цоп(в) - Фурье-образ потного потенциала, Я(в) = У~^ехр[—¡С] - струк-
I
турный фактор. В методе псевдопотенциала Цоп (г) является, вообще говоря, нелокальным потенциалом.
Энергия Хартри в к-пространстве вычисляется как
^ = н ь ~сГ ■
Обмеиио-корреляциоииая энергия вычисляется по формуле
5
Кулоиовская энергия взаимодействия ионов (энергия Маделунга) вычисляется в рамках обычного метода Эвальда.
3. Метод псевдопотенциала
Как указывалось выше, проведение расчетов с учетом всех электронов в базисе плоских воли оказывается невозможным, поскольку это потребовало бы решения секулярного уравнения огромной размерности. Эта проблема успешно решается в развитой теории псевдопотенциала, которая показала себя настолько эффективной, что уже в книге [7] отмечалась возможность получения результатов расчетов электронной структуры и, соответственно, химических, электрических и оптических характеристик ряда соединений, обладающих предсказательной силой.
3.1. Общая концепция метода. Хорошо известно, что большинство физико-химических свойств твердых тел определяется главным образом валентными электронами и в меньшей степени электронами остова. Поэтому представляет интерес рассмотрение только валентных электронов.
Основная задача теории псевдопотенциала (pseudopotential theory) заключается в подборе такого эффективного потенциала псевдопотенциала который позволил бы описать все основные особенности движения валентного электрона в твердом теле, в том числе релятивистские эффекты [7, 53, 54].
Идея появилась довольно давно: она отмечалась уже в работах [55, 56], фактически построение псевдопотенцналов проводится уже в методах присоединенных н ортогонализованных плоских волн [57, 58]. Практически идентичное построение применяется для тяжелых атомов в молекулярных расчетах (оно носит название эффективного потенциала остова effective соте potential) [59]. Однако только после работ Дж. Филлипса и Л. Клейнмана [60, 61] начался поиск различных форм псевдопотенциала, что представляет важную задачу и составляет предмет современной теории.
^рвешЬ у' / \/
гс г
/' ' ^рвешк)
/2/г
Рис. 2. Схематическое гоображепие потенциалов и соответствующих волновых функций в расчетах с учетом всех электронов (сплошные лилии) и в расчетах по методу псевдопо-тепциала (пунктирные лилии)
Псевдопотенциал подбирается таким образом, чтобы рассеяние валентных электронов на псевдопотенциале происходило так же. как и рассеяние валентных электронов на реальном потенциале, созданном ядрами и внутренними электронными оболочками.
Как показано на рис. 2. псевдоволновые функции не имеют узлов внутри области электронного остова (сферы радиуса гс), поскольку псевдопотенциал оказывается слабым в отличие от правильной волновой функции валентного электрона, которая быстро осциллирует в этой области.
Интересно отметить, что таким образом в методе псевдопотенциала находится разумное объяснение того факта, что зонная структура многих твердых тел хорошо описывается в рамках более раннего приближения почти свободного электрона [38]. в котором необоснованно делается предположение о слабости потенциала, действующего на электрон.
В настоящее время подбор параметров псевдопотенциалов осуществляется на основе результатов компьютерных расчетов, в которых учитываются все электроны. Параметры псевдопотенциала подбираются таким образом, чтобы псевдоволновая функция, соответствующая электрону, движущемуся под действием псевдопотенциала. совпадала с волновой функцией электрона, движущегося в реальном потенциале (созданном ядром и электронами внутренних заполненных оболочек) вне сферы некоторого радиуса гс . Блок-схема процедуры подбора параметров псевдопотенциала представлена на рис. 3.
3.2. Построение Филлипса — Клейнмана. Как уже было сказано, метод псевдопотенциала возник изначально как обобщение метода ортогонализованных плоских воли [58]. в котором волновые функции валентных электронов раскладываются по базису ортогонализованных плоских волн:
фОРШ(к + с) = ехр[г(к + с)г] - ^С | ехр[г(к + С)г]М>, (19)
с
где - волновая функция электрона остова, суммирование проводится по всем состояниям электронов остова. Волновая функция валентного электрона в базисе
<^Ввод атомного номера]
Выборг [л]
Решение волнового уравнения исходной задачи_
Равны собственные значения для валентных электронов в задаче с псевдопотенциалом [в исходной задаче?
Да
Равны псевдоволновые функции и собственные волновые^ Нет функции валентных электронов дне сферы радиуса г?
Да
Псевдопотенциал найден
Выбор формы параметризации для псевдопотенциала ^р8еш1о(г)
Выбор значений параметров
Решение волнового уравнения с псевдопотенциалом_
Рис. 3. Блок-схема, описывающая процедуру подбора эффективных параметров псевдо-потепциала
ортогонализованных плоских волн имеет вид
К = "(к
,OPW(b + G) = ^pseudo IT"^)' >,
G
где = С ехр[г(к + в)г] - так называемая псевдоволновая функция.
Эта функция удовлетворяет следующему уравнению [60]:
-1v2 + VH + Vxc +
т т \ / pseudo v / pseudo
Vpseudo = £kVk >
где
pseudo
Von
с
(4 - 4)' Ж
(20) (21)
представляет собой псевдопотенциал.
Снаружи электронного остова, где волновые функции внутренних оболочек обращаются в нуль, потенциал ^8еиао равен Цоп. Псевдопотенциал Ц>8еис1о оказывается намного слабее исходного потенциала Цоп, поскольку слагаемые в уравнении (21) имеют разный знак, поэтому псевдоволновая функция, являющаяся довольно плавной, может быть разложена по базису, составленному из относительно небольшого числа плоских волн.
В настоящее время применяются другие, более эффективные, формы псевдопотенциалов, а именно так называемые поггп-сопяегтпд (N0) и иШшоД ( £/£>) псевдопотенциалы. Эффективный псевдопотеициал должен соответствовать следующим основным требованиям: псевдопотенциал должен быть достаточно слабым, чтобы для корректного разложения псевдоволновой функции требовалось как можно меньше плоских волн: псевдопотеициал, созданный для отдельного атома, должен «работать» при рассмотрении твердых тел: зарядовая плотность, вычисленная на основе псевдоволновых функций, должна как можно точнее воспроизводить реальную плотность заряда.
3.3. Формы псевдопотенциала. Метод построения NC псевдопотенциалов был предложен в работе [62] и усовершенствован в работе [63]. в которой затабу-лпрованы псевдопотенциалы для всех элементов периодической таблицы. В работе [64] представлено альтернативное приближение, основанное на простом аналитическом представлении псевдоволновых функций внутри сферы радиуса rc .
Так. при использовании NC псевдопотенциала псевдоволновая функция (и потенциал) равна истинной волновой функции валентного электрона вне сферы некоторого радиуса rc. Внутри сферы радиуса rc псевдоволновая функция отличается от истинной волновой функции, но норма остается той же самой:
У ^pseudo(r)^pseudo(r) dr = J (г) dr.
r<rc r<rc
При таком подходе псевдопотенциал имеет квазилокальную форму:
^seudо "iseudoPi,
i
где / - угловой момент, Pi - оператор проектирования. Оказывается, что эта сумма быстро сходится, и на практике при построении псевдопотенциала используется приближенная формула:
i
^seudо — Мое + "pseudoPb i=0
где Мое - локальный потенциал, в качестве /тах обычно достаточно выбирать 1 или 2.
Матричные элементы при использовании базиса плоских волн имеют вид [10]:
e-i(k+G)r"¿eudo(r)ei(k+G> dr =
2/ + 1 С
= ^ 4nPi(cos y) J r2M^ud0(r)ji(|k + G|r)ji (|k + G'|r) dr, (22)
где Pi и j - полиномы Лежапдра и сферические функции Бесселя, y - угол между векторами (к + G) и (k + G').
С вычислительной точки зрения при использовании NC псевдопотенциалов имеется несколько довольно неприятных моментов. Во-первых, поскольку псев-допотеициал нелокален, матричные элементы, которые приходится вычислять, зависят не от разности (k + G) — (k+G') = G — G', а отдельно от (k + G) и (k + G'). Это означает, что если базис состоит из m плоских волн, то в общем случае должно быть рассчитано m(m + 1)/2 независимых матричных элементов вместо 8m (как было бы в случае локального потенциала). Во-вторых, такой вид псевдопотенциала довольно неэффективен для расчетов, поскольку в подынтегральном выражении имеется сложная зависимость от волновых векторов. И, наконец, матричные элементы для псевдопотенциала в этом случае зависят от волнового вектора к и, следовательно, должны рассчитываться отдельно для каждой точки зоны Брил-люэна.
Однако, несмотря на эти трудности, полулокальные потенциалы такого типа используются очень широко, что позволяет получить очень хорошие результаты.
В работах [65 67] предложена другая концепция построения псевдопотенциалов. При написании US псевдопотенциалов требуется равенство псевдоволновой
функции и истинной волновой функции вне сферы электронного остова, как и для N0 псевдопотенциалов. Но внутри сферы радиуса гс требуется, чтобы псевдоволновая функция была настолько плавной, насколько это вообще возможно. При этом, разумеется, требование сохранения нормы внутри сферы радиуса гс не может быть удовлетворено. Это обстоятельство вносит некоторые усложнения, однако такой ценой может быть значительно уменьшен размер необходимого базисного набора плоских волн.
Поскольку псевдоволновая функция, очевидно, вообще не является нормированной, то в секулярном уравнении появляются нетривиальные интегралы перекрывания. Кроме того, псевдозарядовая плотность уже не равна просто | ^рвеис!о |2 ) псевдопотенциалы такого типа, подобранные для свободного атома, не всегда дают хорошие результаты при использовании их для атомов внутри твердого тела. Не вдаваясь в пространные объяснения, что выходит за рамки настоящей работы, приведем основные соотношения, которые используются для соответствующих вычислений с применением этой формы псевдопотенцнала. Секулярное уравнение (20) при таком подходе записывается в виде:
(-2V2 + Ун + Ме + ^рвеисю) = (23)
где псевдопотенциал вычисляется по формуле
Мрвеис!о = Мое(г) Опш1ви){вш1-
тп
Оператор перекрывания есть
# = 1 + ^ ?пт|вп)(вт|-тп
Псевдозарядовая плотность определяется выражением:
n(r) = Y1
осе
|"0pseudo(r) | + ^ ] Qnm (г) (^pseudо \Pn) (Pm \^pseudо)
Функции вп = ^2 ^(гУ^Ьт(д, ф) И величины Япт,Япт, Опт определяются в ходе
1т
построения псевдопотенциала.
Таким образом, применение современных формулировок метода псевдопотенциала позволяет проводить самосогласованные расчеты волновых функций валентных электронов внутри кристаллов с использованием базисного набора из относительно небольшого числа плоских воли.
Имеются опубликованные таблицы псевдопотенциалов [63, 68, 69] и открытые базы данных в Интернете [70].
4. Результаты расчетов для простых кристаллических систем
В данном разделе приведены результаты проведенных нами расчетов в рамках описанных выше приближений структурных, электронных, упругих и оптических свойств некоторых простых кристаллов, свойства которых хорошо изучены теоретически и экспериментально, а именно типичных полупроводников углерода, кремния, гермаиия и кристалла никеля. Расчеты были выполнены с использованием программного пакета Quantum ESPRESSO [71, 72]. Псевдопотенциалы для атомов взяты из открытой базы данных [70].
4.1. Зонная структура углерода, кремния и германия. Нами проведены расчеты равновесных параметров решетки, объемных модулей упругости и зонной структуры простых кристаллов типичных полупроводников углерода, кремния н германия.
Эти вещества имеют структуру алмаза, то есть кубическую гранецентрирован-ную кристаллическую решетку с примитивным базисом, состоящим из двух атомов с координатами (0, 0. 0) и (0.25. 0.25. 0.25).
Псевдопотенциалы для атомов углерода С.рг-уЬс.иРГ, кремния Si.pz-vbc.UPF и германия Ge.pz-bhs.UPF [70]. которые использовались в расчетах, построены в приближении локальной плотности с корреляционной энергией в форме Педыо Зунгера [30] в нерелятнвнстском приближении для атомов С и и скалярном релятивистском приближении для атома Се. Псевдопотенцналы имеют ЖС-форму и построены с учетом нелинейной поправки. В расчет включены четыре валентных электрона от каждого атома (пв2пр2, п = 2 , 3 ,4 для С, и Се соответственно).
Значение Еси^ выбрано нами равным 40 Т1у (это соответствует базису, состоящему из ~ 1200 плоских волн) на основании проведенных нами предварительных расчетов, показавших, что такое ограничение базисного набора обеспечивает точность вычисления полной энергии, которая составляет 10-3 Т1у.
Для вычисления интегралов в к-пространстве использовался следующий набор из десяти специальных точек [44] (указаны координаты к-точек и весовые множители):
кх = кз = кб = кг = кд =
7 3 1 8' 8' 8 5 5 1 " 8' 8 3 1 8' 8 3 3 8' 8
3 1 1 8' 8' 8
=
адз =
=
=
=
3 16' 3 ;
32' 3 16' 1 ' 32' 3 ' 32'
2
к4
ке кв кю
7 1 1 8' 8' 8 5 3 3 8' 8' 8 5 1 1 8' 8' 8 3 3 1 8' 8' 8 1 1 1 8' 8' 8
^2 =
и>4 =
аде =
=
^10 =
3 '
32' 3 ' 32' 3 ' 32' 3 ' 32' 1
32'
При расчетах учитывалось фиксированное заполнение нижних зон. Расчеты проводились в спии-иеполяризоваииом варианте метода БИ?.
Нами была рассчитана зонная структура рассматриваемых полупроводников путем самосогласованного решения основного уравнения для точек к-пространства. Зонная структура углерода, кремния и германия показана на рис. 4, 5 и 6 соответственно.
Для кристалла углерода как верхний край валентной зоны, так и нижний край зоны проводимости находятся в точке Г. Ширина запрещенной зоны Ед согласно нашим расчетам составляет 5.629 эВ (экспериментальное значение [75] — 5.40 эВ).
Верхние края валентной зоны кристаллов кремния и германия согласно представленных расчетов находятся в центре зоны Бриллюэна, однако нижние края зоны проводимости лежат в точке Ь зоны Бриллюэна, что соответствует, как отмечено в [75], экспериментальным данным по циклотронному резонансу и оптическому поглощению. Теоретически рассчитанная ширина запрещенной зоны равна 1 . 549 0. 445
Ед 1 . 17 0. 74
W ь г X к г
Рис. 4. Зоппая структура углерода
W L Г X К Г
Рис. 5. Зоппая структура кремния
Рис. 6. Зоппая структура германия
Отмстим, что определение ширины запрещенной зоны проводилось нами в приближении «замороженных» электронный орбиталой, при этом но учитывается около 10% энергии, которая необходима для переброса электрона в зону проводимости, а именно эта величина определяется на эксперименте.
Поиск равновесных параметров решетки проводился методом сопряженных градиентов [51] для ионов для случая варьируемой ячейки. Согласно проведенным расчетам равновесное значение параметра решетки составляет 6.679 (6.74), 10.209 (10.26), 10.623 (10.68) ат. ед. для кристаллов углерода, кремния и германия соответственно (в скобках указаны экспериментальные данные [75]).
Объемные модули упругости определялись по формуле В = — V(ЗР/ЗУ), где V .......... объем кристалла, Р — внутреннее электронное давление. Рассчитанные значения объемных модулей упругости В составляют 4.498 (4.43), 0.879 (0.99), 0.649 0. 77
ках указаны экспериментальные данные [75]).
4.2. Фононная дисперсия кремния. Определение фононной дисперсии предполагает помимо необходимого предварительного расчета энергии кристалла также вычисление динамической матрицы. Колебательные частоты определяются путем нахождения собственных значений матрицы Гессе (вторые производные энергии по смещениям атомов в масс-взвошонных координатах):
det
1 t(H.) 2 , 0
— w =0.
ч-/м7м7 Ж/ жз
Здесь Е^ = Еее + Еем + - полная энергия Борна-Оппенгеймера.
Рис. 7. Фоноииая дисперсия кремния: — по расчетам настоящей работы, • • • по данным работ [83, 84]. в которых проведена некоторая подгонка под экспериментальные данные [85], отмеченные символом о для поперечной моды и символом Д для продольной моды
Гессиан получается дифференцированием сил Гельмана Фейнмана по координатам ядер:
d2Et0t(H-) = Г dn(r;R) dEeN (г; R) d + Г ( d2EeN(r;R) d + d2ENN (R) dRi BRJ J an J dR i +J OR! OR j + SRi OR J '
Таким образом, определение гессиана в рамках метода функционала плотно-
n( ; )
же и линейного отклика на смещения ядер Зп(т; R)/dRi, вычисление которого проводится в рамках теории возмущений для функционала плотности Density-Functional Perturbation Theory (подробно она описана в работах [76 81]).
Литературные данные по результатам a,b initio расчетов фононной дисперсии в кристалле кремния можно найти, например, в работе [82]. Расчеты оказываются настолько точными, что экспериментальные точки практически полностью «ложатся» на рассчитанные кривые. Поэтому в настоящей работе мы сочли необходимым с целыо большей наглядности для демонстрации точности современных расчетов и значительного успеха в развитии вычислительных методов привести результаты не только наших расчетов, но также и ранних расчетов фононной дисперсии кремния.
Так, на рис. 7 приведены результаты ранних расчетов фононной дисперсии для кремния [83, 84], для которых проведена некоторая подгонка к экспериментальным данным [85], и результаты проведенных памп расчетов. Этот рисунок наглядно показывает, насколько точны современные расчеты, которые были проведены «из первых принципов», без каких-либо подгонок. В частности, проведенные памп расчеты гораздо лучше соответствуют экспериментальным данным для акустических мод.
4.3. Магнитные свойства кристалла никеля. Интересным представляется применение описанного подхода для изучения магнитных свойств материалов. Расчеты, результаты которых приведены в этом параграфе, выполнялись нами с целыо определения равновесной намагниченности кристалла никеля.
Свободный атом Ni [3d84s2] обладает магнитным моментом, обусловленным характером заполнения внутренней 3й-оболочки. Хорошо известно, что по мере того,
Рис. 8. Энергетические зоны ферромагнитного никеля (сплошная кривая соответствует электронам со спппом «вверх», пунктирная со спилом «вниз»)
2-
2-
Рис. 9. Плотность состояний для электронов со спилом «вверх» (сплошная лилия) и со спилом «вниз» (пунктирная лилия) в кристалле ферромагнитного никеля
как атомы объединяются в кристалл, происходит уширение частично заполненных зон н возрастает эффект, вызывающий увеличение энергии с ростом намагниченности. Ферромагнетизм кристаллов в равновесных условиях соответствует ситуации, когда частично заполненные зоны, приводящие к намагниченности, достаточно узки. Поэтому намагниченность кристалла является довольно чувствительным параметром по отношению к точности описания электронной структуры и значению параметра решетки, как отмочено в книге [74].
Расчеты проводились в приближении локальной спиновой плотности метода В¥Т. Интегрирование в к-пространство проводилось с использованием однородной сетки к-точек (144 точки).
Никель имеет гранецентрированную кубическую кристаллическую решетку. Псевдопотонциал для атома никеля NiUS.RRKJ3.UPF [70], который использовался нами, построен на основе скалярных релятивистских расчетов в локальном приближении метода функционала плотности с корреляцией в форме Педыо Зунгера [30]. Псевдопотонциал имеет ив форму и построен с учетом нелинейной поправки. Значение параметра выбрано равным 24 11у (~140 плоских волн
в базисе).
Зонная структура, рассчитанная для решетки с параметром а = 6.48 ат. ед., показана на рис. 8. Уровень Ферми пересекает частично заполненные зоны, причем заполнение зон, в которых находятся электроны со спином «вверх» и спином «вниз», оказывается различным, и кристалл никеля обладает отличным от нуля магнитным моментом.
Для определения равновесной намагниченности кристалла нами стандартным методом [75] была вычислена плотность электронных состояний (рис. 9). Согласно проведенным вычислениям магнитный момент, приходящийся па элементарную ячейку ферромагнитного никеля, составляет 0.54мв (мв - магнетон Бора), что находится в хорошем количественном согласии с экспериментальной величиной
0.6. [74].
Заключение
Подводя итоги, нам хотелось бы отметить, что описанный в настоящей работе подход для проведения расчетов «из первых принципов» дает возможность изучения структурных н электронных свойств кристаллов без привлечения каких-либо дополнительных сведений, кроме их атомной структуры. Результаты расчетов имеют хорошее количественное согласие с экспериментальными данными. Это открывает широкие возможности для теоретической интерпретации многих «твердотельных» экспериментов: ЯМР- и ЭПР-спектроскопнн, рентгеноструктурного анализа н оптической спектроскопии.
Авторы благодарят Д.В. Чачкова, который взял на себя труд, связанный с компиляцией программы.
Работа выполнена при финансовой поддержке специального стипендиального фонда Республики Татарстан (указ Президента Республики Татарстан УП-17 от 19 января 2009 г.).
Summary
V. V. Klekovkina, R.M. Aminova. Ab initio Calculations of Structural and Electronic Properties of Crystal Solids in Density Functional and Pseudopotential Approach in Momentum Space: Details and Examples.
The paper discusses the methods of density functional theory, pseudopotential theory and plane-wave basis set approach. These theories serve as a basis for quantum-chemical study of structural and physical-chemical properties of crystals. Structural and electronic properties of carbon, silicon, germanium and nickel crystals have been calculated. The results of calculations are compared with experimental data.
Key words: density functional theory, local density approximation, pseudopotential method, plane-wave basis set.
Литература
1. Catlow C.R.A., Price G.D. Computer modeling of solid-state inorganic materials // Nature. 1990. V. 347. P. 243 248.
2. Zhang L., Niu Y., Li Q. et al. Ab initio prediction of superconductivity in molecular metallic hydrogen under high pressure // Solid State Commuu. 2007. V. 141. P. 610 614.
3. Mattioli G., Filippone F., Giannozzi P. et al. Theoretical design of coupled organic-inorganic systems // Pliys. Rev. Lett. 2008. V. 101. Art. 126805.
4. Feng H.-J,, Liu F.-M. First-principles prediction of coexistence of magnetism and ferro-electricity in rhombohedral Bi2FeTiC>6 // Phys. Lett. A. - 2008. - V. 372. - P. 1904-1909.
5. Zhang L., Wang Y., Cui T. et al. First- principles study of the pressure-induced rutile-
22
6. An J., Subetli A., Singh D.J. Ab initio plionon dispersions for PbTe // Solid State Commun. 2008. V. 148. P. 417 419.
7. Хейне В., Ko:m М., У:тр Д. Теория исевдопотепциала. М.: Мир. 1973. 557 с.
8. Займам Дон:. Вычисление блоховских функций. М.: Мир. 1973. 159 с.
9. Singh D.J., Nordstriim L. Planewaves, Pseudopotentials and LAPW Method. N. Y.: Springer. 2006. 134 p.
10. Ihm J., Zunger A., Cohen M.L. Momentum-space formalism for the total energy of solids // J. Phys. C.: Solid State Phys. 1979. V. 12. P. 4409 4422.
11. Denteneer P.J.H., van H tiering en W. The pseudopotential-density-functional method in momentum space: details and test cases // J. Phys. C.: Solid State Phys. 1985. V. 18. P. 4127 4142.
12. Peressi M., Baltleresehi A., Baroni S. Ab initio studies of structural and electronic properties // Characterization of Semiconductor Het.erost.ruct.ures and Nanostructures / Ed. C. Lamberti. Elsevier, 2008. P. 17 54.
13. Jones R.O., Gunnarsson O.P. The density functional formalism, its applications and prospects // Rev. Mod. Phys. 1989. V. 61, No 3. P. 689 740.
14. Kohn W., Beeke A.D., Parr R.G. Density Functional Theory of Electronic Structure // J. Chem. Phys. 1996. V. 100. P. 12974 12980.
15. Nagy A. Density functional. Theory and application to atoms and molecules // Phys. Rep. 1998. V. 298. P. 1 79.
16. Kohn W. Nobel Lecture: Electronic structure of matter wave functions and density functional // Rev. Mod. Phys. 1999. V. 71, No 5. P. 1253 1266. = Кои В. Электронная структура вещества волновые функции и функционалы плотности // Усп. физ. паук. 2002. Т. 172, Л» 3. С. 336 348.
17. Lecture Notes in Physics: Density Functional Theory. N. Y.: Springer, 1983. 301 p.
18. Теория неоднородного электронного газа. М.: Мир, 1987. 400 с.
19. Topics in Current Chemistry: Density Functional Theory I. Functionals and Effective Pot.ent.als. Berlin Heidelberg: Springer, 1996. 224 p.
20. Topics in Current Chemistry: Density Functional Theory II. Relat.ivist.ic and Time Dependent Extensions. Berlin Heidelberg: Springer, 1996. 209 p.
21. Lecture Notes in Physics. Density Functionals: Theory and Applications. N. Y.: Springer, 1998. 194 p.
22. Hohenberg P., Kohn W. Inhomogeneous Electron Gas // Phys. Rev. 1964. V. 136, No 3B. P. 864 876.
23. Highlight of Condensed-Matter Theory / Ed. by F. Bassani, F. Fumi, M.P. Tosi. Elsevier Science Ltd, 1985. 940 p.
24. Levy M. Electron densities in search of Hamilt.oiiians // Phys. Rev. A. 1982. V. 24, No 2. P. 1200 1208.
25. Physics as Natural Philosophy. Essays in Honor of Laszlo Tisza on his 75t.li Birthday / Ed by A. Sliimoiiy, H. Feslibacli. Cambridge, Mass.: MIT Press, 1982. 448 p.
26. Janak J.F. Proof that dE/dni = et in density functional theory // Phys. Rev. B. -
1978. V. 18, No 12. P. 7165 7168.
27. Von Barth U., Hedin L. Local-density theory of mult.plet. structure // Phys. Rev. A.
1979. V. 20, No 4. P. 1693 1703.
28. Soma S.F., Fernandas P.A., Ramos M.J. General Performance of Density Functional // J.Phys. Chem. A. 2007. V. 111. P. 10439 10452.
29. Kohn W., Sham L. Self-consistent, equations including exchange and correlation effects // Phys. Rev. 1965. V. 140, No 4A. P. 1133 1138.
30. Perdew J.P., Zunger A. Self-interaction correction to density-functional approximations for many-electron systems // Phys. Rev. B. 1981. - V. 23, No 10. P. 5048 5079.
31. Ceperley D.M. Ground state of the fermion one-component, plasma: A Monte-Carlo study in two and three dimension // Phys. Rev. B. 1978. V. 18, No 7. P. 3126 3138.
32. Slater J. C. Simplification of the Hart.ree Fock Method//Phys. Rev. B. 1951. V. 81, No 3. P. 385 390.
33. Gunnarsson O., Lundqvist B.I. Excange and correlation in atoms, molecules, and solids by the spin-densit.y-funct.ional formalism // Phys. Rev. B. 1976. V. 13, No 10. P. 4274 4298.
34. Langreth B.C., Perdew J.P. Excange-correlat.ion energy of a metallic surface: Wave-vector analysis // Phys. Rev. B. 1977. V. 15, No 6. P. 2884 2901.
35. Эва/рестов P.А. Кваптовохимические методы в теории твердого тела. Л.: Изд-во Лепипгр. ун-та, 1982. 280 с.
36. Немогикаленко В.В., Кучеренко Ю.Н. Методы вычислительной физики в теории твердого тела. Электронные состояния в пеидеальпых кристаллах. Киев: Наукова думка, 1986. 296 с.
37. Fundamental Materials Research: Electronic Properties of Solids Using Cluster Methods. N. Y.: Kluwer Academic, 2002. 202 c.
38. Киттель Ч. Квантовая теория твердых тел. М.: Наука, 1967. 491 с.
39. Louie S.G., Но К. М., Cohen M.L. Self-consistent, mixed-basis approach to the electronic structure in solids // Phys. Rev. B. 1979. V. 19, No 4. P. 1774 1782.
40. Singh D. Ground-state properties of lanthanum: Treatment, of extended-core states // Phys. Rev. B. 1991. V. 43, No 8. P. 6388 6392.
41. Evarestov R.A. Quantum Chemistry of Solids. The LCAO First-Principles Treatment, of Crystals. Berlin Heidelberg: Springer, 2007. 571 p.
42. Бандура А.В., Эоарестоо P.А. Неэмпирические расчеты кристаллов в атомном базисе (с использованием Интернет-сайтов и параллельных вычислений). СПб.: Изд-во С.-Петерб. ун-та, 2004. 232 с.
43. Эоарестоо Р.А., Смирнов П.М. Методы теории групп в квантовой химии твердого тела. Л.: Изд-во Лепипгр. ун-та, 1987. 375 с.
44. Chadi D.J., Cohen M.L. Special Points in the Brillion Zone // Phys. Rev. B. 1973. V. 8, No 12. P. 5747 5753.
45. Balderesehi A. Mean-Value Point, in the Brillion Zone // Phys. Rev. B. 1973. V. 7, No 12. P. 5212 5215.
46. Monkhorst H.J., Pack J.D. Special points for Brillion-zone integrations // Phys. Rev. B. 1976. V. 13, No 12. P. 5188 5192.
47. Gill P.E., Murray W., Wright M.H. Practical Optimization. London: N. Y.: Acad. Press, 1981. 417 p.
48. Gillan M.J. Calculation of the vacancy formation energy in aluminium // J. Pliys.: Condens. Matter. 1989. V. 1. P. 689 711.
49. Stick I., Car R., Parrinellu M., Barirni S. Conjugate gradient minimization of the energy functional: A new method for electronic structure calculation // Pliys. Rev. B. 1989. V. 39, No 8. P. 4997 5004.
50. Teter M.P., Payne M.C., Allan B.C. Solution of Sclirodinger's equation for large systems // Pliys. Rev. B. 1989. V. 40, No 18. P. 12255 12263.
51. Payne M.G., Tatar M.P., Allan D.C. et al. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients // Rev. Mod. Pliys. 1992. V. 64, No 4. P. 1046 1096.
52. Holzschuh E. Convergence of momentum space, pseudopotential calculations for Si // Pliys. Rev. B. 1983. V. 28, No 12. P. 7346 7348.
53. Pickett W.E. Pseudopotential methods in condenced matter applications // Computer Pliys. Rep. 1989. V. 9. P. 115 197.
54. Austin B.J., Heine V., Sham L.J. General Theory of Pseudopotentials // Pliys. Rev. 1962. V. 127, No 1. P. 276 282.
55. Helhnann H. A new approximation method in the problem of many electrons // J. Cliem. Pliys. 1935. V. 3. P. 61 63.
56. Gombás P. Uber die metallisclie Biiidung // Z. Pliys. 1935. V. 94. P. 473 475.
57. Slater J.C. Wave functions in a Periodic Potential // Pliys. Rev. 1937. V. 51. P. 846 851.
58. Herring C. A New Method for Calculating Wave Functions in Crystals // Pliys. Rev. 1940. V. 57. P. 1169 1177.
59. Preuss H. Untersuchungen zum kombiniert.en Nalieruiigsverfahreii // Z. Nat.urf. 1955. V. 10a. P. 365 367.
60. Phillips J.C. Energy-Band Extrapolation Scheme Based on Pseudopotentials // Pliys. Rev. 1958. V. 112, No 3. P. 685 695.
61. Phillips J.C., Kleinman L. New Method for Calculating Wave Functions in Crystals and Molecules // Pliys. Rev. 1959. V. 116, No 2. P. 287 294.
62. Hamann D.R., Shlüter M., Chiang C. Norm- Conserving Pseudopotentials // Pliys. Rev. Lett. 1979. V. 43, No 20. P. 1494 1497.
63. Bachelet G.B., Hamann D.R., Shlüter M. Pseudopotential that work: From H to Pu // Pliys. Rev. B. 1982. V. 26, No 8. P. 4199 4228.
64. Kerker G.P. Non-singular atomic pseudopotentials for solid state applications // J. Pliys. C.: Solid St. Pliys. 1980. V. 13. P. L189 L194.
65. Vanderbilt D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism // Pliys. Rev. B. 1990. V. 41, No 11. P. 7892 7895.
66. Laasunen K., Car R., Lee C., Vanderbilt D. Implementation of ultrasoft pseudopotentials in ab initio molecular dynamics // Pliys. Rev. B. 1991. V. 43, No 8. P. 6796 6799.
67. Laasunen K., Pasquarellu, Car R., Lee Ch., Vanderbilt D. Car-Parrinello molecular dynamics with ultrasoft pseudopotentials // Pliys. Rev. B. 1993. V. 47, No 16. P. 10142 10153.
68. Gome X, Stumpf R., Scheffier M. Analysis of separable potentials // Pliys. Rev. B. 1991. V. 44, No 16. P. 8503 8513.
69. Goedecker S., Tatar M., Huttar J. Separable dual-space Gaussian pseudopotentials // Pliys. Rev. B. 1996. V. 54, No 3. P. 1703 1710.
70. Pseudopotentials Database. URL: http://www.pwscf.org/pseudo.lit.iiil, свободный.
71. Scantlolo S., Giannozzi P., C. Cavazzoni C., de Gironcoly S., et al. First-principles codes for Computational Crystallography in the Quantum-ESPRESSO package // Z. Krist.al-logr. 2005. V. 220. P. 574 579.
72. Giannozzi P. et al. Quantum ESPRESSO, ver. 4.0.4. 2008. URL: lit.tp://www.quantum-espresso.org/, свободный.
73. Вычислительные методы в теории твердого тела. М.: Мир, 1975. 400 с.
74. Слэтер Дон:. Методы самосогласованного поля для молекул и твердых тел. М.: Мир, 1978. 662 с.
75. Киттель Ч. Введение в физику твердого тела. М.: Наука, 1978. 791 с.
76. Gome X. First-principles responses of solids to atomic displacements and homogeneous electric fields: Implementation of a conjugate-gradient, algorithm // Pliys. Rev. B. 1997. V. 55, No 16. P. 10337 10354.
77. Gonze X., Lee C. Dynamical matrices. Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory // Pliys. Rev. B. 1997. V. 55, No 16. P. 10355 10368.
78. Dal Corso A., Pasquarello A., Baldereschi A. Density-functional perturbation theory for lattice dynamics with ult.rasoft. pseudopotentials // Pliys. Rev. B. 1997. V. 56, No 18. P. 369 372.
79. Dal Corso A. Densit.y-funct.ional perturbation theory with ult.rasoft. pseudopotentials // Pliys. Rev. B. 2001. V. 64. Art.. 235118.
80. Baroni S., de Gironcoli S., Dal Corso A., Giannozzi P. Plionons and related crystal properties from densit.y-funct.ional perturbation theory // Rev. Mod. Pliys. 2001. V. 73. P. 516 562.
81. Giannozzi P., Baroni S. Density-Functional Perturbation Theory // Handbook of Materials Modeling. N. Y.: Springer, 2005. P. 195 214.
82. Giannozzi P., de Gironcoli S., Pavone P., Baroni S. Ab initio calculation of plioiion dispersions in semiconductors // Pliys. Rev. B. 1991. V. 43, No 9. P. 7231 7242.
83. Martin R.M. Lattice vibrations in silicon: microscopic dielectric model // Pliys. Rev. Lett. 1968. V. 21, No 8. P. 536 539.
84. Martin R.M. Dielectric Screening Model for Lattice Vibrations of Diamond-Structure Crystals // Pliys. Rev. 1969. V. 186., No 3. P. 871 884.
85. Palevsky H., Hughes D.J., Kley W. et al. Lattice vibrations in silicon by scattering of cold neutrons // Pliys. Rev. Lett.. 1959. V. 2, No 6. P. 258 259.
Поступила в редакцию 16.03.09
Клековкина Вера Вадимовна магистрант направления «Теоретическая и математическая физика» Казанского государственного университета. Е-шаП: ьега_ klekovkina0mail.ru
Аминова Роза Мухаметовна доктор химических паук, профессор кафедры химической физики Казанского государственного университета. Е-шаП: Roza.Aminovaeksu.ru
