A study of carbon nanotubes energetics using orbital free method in the frame-work of the density functional theory

Full Text

Введение Углеродные нанотрубки (УНТ), открытые в конце прошлого столетия, находятся в центре внимания многих исследователей, которые интересуются их атомной и электронной структурой, механическими свойствами, их взаимодействием друг с другом и с другими материалами и т.п. Однако, несмотря на огромное число работ, посвященных изучению УНТ, до сих пор остаются непонятными их некоторые важные особенности. В частности, без ответа остается вопрос: чем определяется диаметр нанотрубок? Наиболее часто в литературе просто констатируется, что диаметр УНТ лежит в пределах 0,7-1,6 нм. В то же время существуют данные, что можно получить УНТ с минимальным диаметром 0,3 нм [1] и с максимальным диаметром 12 нм [2]. Экспериментальное распределение УНТ по величине диаметра от 1,0 до 1,5 нм дано в работах [3, 4]. В нашей недавней работе [5] мы исследовали энергетику УНТ методом Кона-Шэма [6] и нашим безорбитальным (БО) методом с использованием псевдопотенциалов. В данной работе мы проводим такое же исследование, с целью уточнения результатов, используя более однозначный полноэлектронный подход и системы с большим числом атомов. 1. Основные положения безорбитального подхода Теория функционала плотности (ТФП) [7] основывается на представлении о том, что основное (невозбужденное) состояние квантовой системы может описано без использования волновых функций (орбиталей), с помощью единственной функции - электронной плотности ρ, которая находится путем минимизации функционала энергии E [ρ]: где ε(ρ) - плотность энергии; V (r) - внешний потенциал; φ - электростатический электронный потенциал Хартри, εex-c и εkin - обменно-корреляционная и кинетическая энергии (на электрон). Минимизация уравнения (1) при условии ∫ρ(r)dr = N, где N - количество электронов в системе, означает решение следующего уравнения: (2) где μ - электронный химический потенциал; и ρex-c и ρkin - обменно-корреляционный и кинетический потенциалы, соответственно, Внешний потенциал V (r) в нашем полноэлектронном случае есть сумма атомных потенциалов. Для вычисления обменно-корреляционного потенциала существуют достаточно реалистические приближения (например, [8; 9; 10]). Потенциал Хартри может быть вычислен с помощью Фурье-преобразований или из уравнения Пуассона. Что касается кинетического потенциала (и соответственно кинетической энергии), то здесь теория функционала плотности столкнулась с проблемой, которая по сути не решена до сих пор. Кон и Шэм [6] предложили компромиссный способ вычисления кинетической энергии с использованием одночастичных волновых функций путем решения уравнения, которое получило название уравнение Кона-Шэма (КШ): (3) где ψi (r) - волновые функции, или орбитали; εi - энергия i-го состояния. Подход КШ оказался весьма эффективным и до сих пор широко применяется для моделирования многоатомных систем. Однако, в силу использования волновых функций, число которых растет с ростом числа атомов, его вычислительные возможности весьма ограничены и к настоящему времени практически достигли своего предела - несколько сотен атомов (при использовании псевдопотенциалов) и несколько десятков атомов (в полноэлектронном варианте). В то же время, современные нанотехнологии заинтересованы в моделирования систем, состоящих из десятков и сотен тысяч атомов. Рассмотрим одиночный атом, равновесная полная электронная плотность которого может быть легко вычислена методом КШ. В соответствии с уравнением (2) мы можем написать уравнение для нахождения одноатомного кинетического потенциала μatkin: (4) где Z - полный заряд ядра атома. Электронный химический потенциал одиночного атома мы полагаем равным нулю; то есть, электронный химический потенциал системы взаимодействующих атомов будем вычислять относительно системы изолированных атомов. Если мы знаем μatkin(r) и ρ(r) для одиночного атома (а их мы легко находим методом КШ с помощью подходящей полноэлектронной программы), то мы можем найти для одиночного атома потенциал μatkin(ρ). В работах [11; 12] мы, следуя этой методике, нашли кинетические потенциалы для атомов B, C, N, O. В этих же работах был рассмотрен вопрос о том, как, зная кинетический потенциал одиночного атома, найти кинетический потенциал системы атомов. Мы убедились, что в случае системы одинаковых атомов хорошо работает простая формула: (5) где d - среднее расстояние между ближайшими атомами; γ - константа, характерная для данного типа атомов; Nval - число валентных электронов в одиночном атоме. Для углерода мы использовали γ = 12. Выше было сказано, что существуют реалистичные приближения для вычисления обменно-корреляционных потенциалов. Для систем, содержащих атомы малых размеров (а углерод относится к их числу), наиболее реалистичным является приближение обобщенных градиентов (Generalized Gradient Approximation - GGA) [10]. Однако применение GGA для полноэлектронных задач сталкивается с определенным трудностями, связанными с необходимостью выполнять численное дифференцирование функций сильно локализованных в пространстве. Такие функции соответствуют остовным состояниям, которые не участвуют в межатомных взаимодействиях и обычно рассматриваются как «замороженные». Такое разделение плотности электронов помогает избежать операций с интенсивными резкими пиками и построить реалистичные компьютерные коды (например, пакет Elk [13], который в данной работе мы используем для нахождения полноэлектронного кинетического потенциала одиночного атома углерода). Следуя этой технике, мы разделяем атомную плотность на остовную и валентную части. Однако это не позволяет избежать всех проблем при использовании GGA, поскольку валентная плотность атома углерода все-таки локализована достаточно сильно. Для преодоления этих проблем мы разработали специальную методику дискретной минимизации функционала энергии, которая описана в приложении. 2. Вычисления и обсуждение результатов Минимум энергии атомной системы мы находим путем нахождения плотности ρ, которая обращает в ноль выражение (3). С учетом того, что остовные плотности не изменяются при межатомных взаимодействиях, мы можем записать итерационное уравнение для валентной плотности ρval: ρval (r) = ρval (i - 1) + KF (i - 1) ρval (i - 1), где K - итерационный параметр, управляющий сходимостью процедуры. На начальном шаге (i = 0) плотность ρval(0) представляет собой сумму равновесных атомных валентных плотностей: Отметим, что полная энергия системы атомов Etot является суммой энергии отталкивания ядер Erep, кулоновской энергии EC, энергии Хартри EH, обменно-корреляционная энергии Eex-c и кинетической энергии Ekin: Мы рассматривали нанотрубки конечной длины с открытыми концами типа «кресло» (armchair). Диаметр трубок D варьировался от 0,68 до 1,50 нм, а количество входящих в них атомов изменялось от 80 до 320. Были рассмотрены три набора трубок: длина трубок в первом наборе составляла 0,87 нм, во втором - 1,36 нм, в третьем - 1,86 нм. Моделирование проводилось в кубическом объеме со стороной 5 нм, который разбивался при интегрировании на ячейки с помощью трехмерной сетки 200 × 200 × 200. Геометрические параметры исследованных трубок приведены в табл. 1. Таблица 1 Геометрические параметры исследованных УНТ [Geometry parameters of studied CNTs] Диаметр, нм [Diameter, nm] Длина, нм [Length, nm] 0,87 Набор № 1 [Set 1] 1,36 Набор № 2 [Set 2] 1,86 Набор № 3 [Set 3] Количество атомов [Number of atoms] 0,68 80 120 160 0,82 96 144 192 0,95 112 168 224 1,06 128 192 256 1,22 144 216 288 1,36 160 240 320 Результаты расчетов энергии связи представлены на рис. 1. Рис. 1. Зависимость энергии связи (на атом) от диаметра нанотрубок. Цифры 1, 2, 3 соответствуют номерам исследованных наборов УНТ Fig. 1. Dependence of the binding energy (per atom) on the nanotubes diameter. Marks 1, 2, 3 correspond to the set numbers Из рис. 1 видно, что во всех случаях наблюдается минимум энергии связи при величине диаметра нанотрубки около одного нанометра. При этом величина энергии связи растет с увеличением длины трубок, оставаясь вблизи 8 эВ. В работе [5] мы провели методом Кона-Шэма расчет энергии УНТ бесконечной длины и нашли величину энергии связи около 10 эВ. На рис. 1 обращает на себя внимание тот факт, что энергия нанотрубок слабо изменяется при изменении диаметра. Это объясняет большой разброс экспериментальных значений диаметров УНТ, о котором говорилось во Введении. Следует также отметить, что полученные нами ранее результаты с использованием псевдопотенциального приближения [5] также свидетельствуют о том, что диаметр D ≈ 1 нм является наиболее выгодным для углеродных нанотрубок. 4. Выводы Таким образом, нами показано, что наш безорбитальный метод позволяет описать энергетику углеродных нанотрубок в согласии с экспериментом. А именно: мы нашли, что энергия связи, приходящаяся на один атом, имеет минимум при диаметре D ≈ 1 нм, а ее абсолютная величина равна 8 эВ. Приложение Дискретный подход к решению вариационной задачи теории функционала плотности в реальном пространстве Теория функционала плотности (ТФП) является одним из подходов к нахождению решения задач квантовой механики для многоатомных систем. Напомним, что основным положением ТФП является утверждение о том, что стационарное состояние квантовой системы определяется ее равновесной электронной плотностью ρ(r) (r ∈ R3), которая минимизирует функционал электронной энергии данной системы Eel: (П.1) при условии постоянства количества электронов Nel = ∫ρ(r)dr в системе, то есть удовлетворяет вариационному уравнению (П.2) из которого следует (П.3) Здесь ε(ρ) - плотность полной электронной энергии; μ(r) - химический потенциал электрона, т.е. потенциал одного электрона в системе. Для бесконечных квантовых систем μ - константа, для конечных же систем (электронная плотность которых стремится на бесконечности к нулю) μ = μ(r) есть функция координат. Величина ε(ρ) складывается из нескольких частей: (П.4) где V (r) - внешний потенциал; φ(r) - электростатический потенциал Хартри, εkin(ρ) и εex-c(ρ) - плотности кинетической и обменно-корреляционной энергий. Используя (П.3) и (П.4), получаем V (r) + φ(r) + μkin(r) + μex-c(r) - μ(r) = 0, (П.5) где В зависимости от типа атомов в системе для вычисления обменно-корелляционной энергии могут быть использованы различные приближения. Если электронная плотность в системе слабо зависит от координат, можно ограничиться приближением локальной плотности (Local Density Approximation, LDA). В этом случае (П.6) где γ = 0,1423; β1 = 1,0529; β2 = 0,3334; A = 0,0311; B = -0,048; C = 0,002; D = -0,0116. Однако это приближение плохо работает для атомов с сильно локализованными электронами (типа углерода) и еще хуже для атомов с d-электронами (типа переходных металлов). В этих случаях необходимо использование приближения, содержащего производные плотности (градиенты) (Generalized Gradient Approximation, GGA). Возникающая при этом необходимость нахождения производных плотности приводит к большим вычислительным проблемам, которые усугубляются тем, что при решении вариационной задачи мы имеем дело не только с обменно-корелляционной энергией εex-c(ρ), но и с ее производной по плотности, то есть с обменно-корреляционным потенциалом, содержащим вторые производные плотности: εex-c(ρ) = εex(ρ) + εc(ρ); εex(ρ) = ρ(r)εexunifFex; (П.7) εc(ρ) = ρ(r)(εcunif + H), где Для сравнения мы провели расчеты обменно-корреляционного потенциала c использованием приближений как в GGA, так и в LDA для димера C2. В случае LDA вычисление потенциала не вызывает особых сложностей, и может быть проведено в декартовых координатах с необходимой точностью. Вычисление же GGA-потенциала, содержащего градиенты плотности, сталкивается с известной проблемой невысокой точности численного дифференцирования, для преодоления которой необходимо уменьшать шаг дифференцирования: мы используем шаг 0,1 атомной единицы (1 а.е. = 0,0529 нм) и дальнейшее его уменьшение технически невозможно. На рис. П.1 показано поведение обменно-корреляционного потенциала в LDA и GGA приближениях для димера C2. Из этого рисунка видно, что в случае LDA обменно-корреляционный потенциал ведет себя гладко, а в случае GGA на соответствующей кривой имеются острые пики, обусловленные погрешностями численного дифференцирования. А так как численное решение вариационной задачи (П.3) производится итерационным способом, то погрешности численного дифференцирования накапливаются в процессе итераций, что приводит к расходимости итерационного процесса. Рис. П.1. Сравнение поведения обменно-корреляционного потенциала в приближениях LDA и GGA для димера C2 Fig. A.1. Comparing of the LDA and GGA exchange-correlation potentials for the C2 dimer Поэтому для вычисления GGA-потенциала и его использования в итерационной процедуре мы разработали специальную (дискретную) методику. Обозначим через F (ρ(r)) полный потенциал атомной системы, состоящей из Nat атомов, в которой n-атом расположен в точке Rn ∈ Ω ⊂ R3 и относится к типу Tn (тип - это химический элемент): (П.8) где n - номер атома в системе; Z (Tn) - полный ядерный заряд атома с номером n. Перепишем вариационную задачу (П.3) в терминах F (ρ(r)): F (ρ(r)) - μ(r) = 0. (П.9) Наша цель - найти плотность ρ которая удовлетворяет этому соотношению. Для нахождения ρ(r) будем использовать метод последовательных приближений (метод простых итераций): где ρ(r; i) - плотность на i-й итерации; K - параметр, контролирующий процедуру сходимости, а в качестве химического электронного потенциала μ(r; i) берется среднее значение полного потенциала, вычисленного на итерациях i и i - 1: (П.11) На нулевой итерации плотность системы задается как сумма сферических плотностей одиночных невзаимодействующих атомов (П.12) где ρsphere (n; radius) - равновесная плотность атома с номером n, заданная в сферической системе координат. Найти атомные плотности можно при помощи любой программы, обеспечивающей нахождение равновесного состояния одиночного атома в рамках полноэлектронного метода КШ. Теперь вычислим полный потенциал атомной системы F (r; i) при i ≥ 0. На нулевой итерации воспользуемся аддитивностью электростатического потенциала. Тогда, согласно (П.8) (П.13) Для нахождения кинетического потенциала будем использовать подход, описанный в работах [11; 12], в основе которого лежат расчеты одиночных атомов методом КШ и использование для них условия F (ρ(r)) = 0. Для нахождения μex-c(ρ(r; 0)), куда входит суммарная плотность, запишем равенство (П.12) в виде суммы сферической плотности одного произвольного атома и вклада всех остальных: (П.14) где n (r) - номер ближайшего к точке r атома. Тогда Представляя μex-c(ρ(r; 0)) в виде ряда Тейлора по функции ρ(r; 0), в каждой точке r ∈ Ω, получаем (П.15) В этой формуле величины вычисляются в сферической системе координат, а μex-c(ρsphere) вычисляется по формулам (П.5) и (П.7). На последующих итерациях i ≥ 1, согласно (П.8) и (П.10), (П.16) где (П.17) Запись [A] означает значение логического выражения A. Рассмотрим вычисление μex-c(ρ(r; i)) в GGA-приближении. При i ≥ 1 плотность ρ(r; i) уже не является суммой сферических плотностей. Разложим GGA-потенциал в ряд Тейлора по функции ρ используя (П.10) и (П.17) и вычисляя dμex-c(ρ)/dρ в декартовой системе координат: Поставим в соответствие итерационным соотношениям (П.10)-(П.18) дискретный аналог. Для этого зададим в расчетной области Ω ∈ R3 две сетки Sphs (Ω), Carth (Ω): где последовательность радиусов (rir > 0) - геометрическая прогрессия с шагом s > 1; последовательности углов получены разбиением сферы с помощью вписанного икосаэдра, а затем дроблением граней икосаэдра на более мелкие грани. Построенные таким образом сетки Sphs (Ω; n) обладают однородностью - все вершины имеют одинаковое число соседей. Введем обозначение: In(f) - интерполяционный полином для функции f совпадающий с f, в n точках из множества I. В каждом узле (n, ir, j) сетки Sphs (Ω) разностное уравнение, соответствующее (П.15), представим в виде (П.19) где, согласно (П.10) и (П.14), (П.20) В узлах (k, l, m) декартовой сетки Carth (Ω) величины GGA-потенциал и Δρ вычислим следующим образом. Положим где куб с длиной стороны, равной h, и с центром в точке (kh, lh, mh). Определим также множество ℜ (k, l, m) - множество ближайших к вершинам куба узлов из сетки Sphs (Ω). Тогда на нулевой итерации возьмем значения GGA-потенциала и Δρ в узлах декартовой сетки как средние значения в соответствующих областях: если N (k, l, m) > 0, то в противном случае (П.22) Верхние индексы в величинах в правых частях равенств указывают на номер итерации). Плотность ρ на нулевой итерации в узлах декартовой сетки находится по формуле (П.12): (П.23) Как видим, на нулевой итерации в разностных формулах (19)-(22) обменно-корреляционный GGA-потенциал определяется через GGA-потенциал отдельного атома μex-c(ρsphere). Он находится с высокой точностью в сферической системе координат, поскольку сама плотность одного атома зависит только от радиуса сферы с центром в точке расположения атома, и в областях с большим градиентом плотности (вблизи расположения атома) расчетный шаг изменения радиуса очень мал. Для вычисления GGA-потенциала на итерации с номером i ≥ 1 обратимся к соотношениям (П.12), (П.13), (П.16)-(П.18). Положим Тогда разностные уравнения, соответствующие соотношениям (П.11), (П.13), (П.16), (П.17) и (П.18), запишем в виде (П.24) (П.25) (П.26) (П.27) При вычислении [F]0k, l, m используется φsphere(n; radius) - потенциал Хартри атома с номером n, заданный в сферической системе координат. Найти его можно так же, как равновесную плотность одного атома при помощи любой программы, обеспечивающей нахождение равновесного состояния одиночного атома в рамках полноэлектронного метода КШ. Таким образом, задача сводится к нахождению рекуррентных последовательностей в каждой точке r декартовой сетки Carth (Ω) в следующем порядке: на 0 - итерации по формулам (П.21), (П.22), (П.23), (26) вычисляются [ρ]0r; [Δρ]0r; [μex-c (ρ)]0r; [F]0r; на i - итерации при i > 0 по формулам (П.25) вычисляются [Δρ]ir и [ρ]ir. Далее вычисляются [μex-c(ρ)]ir, [F]ir, [μ]ir, - соотношения (П.24), (П.26), (П.27). В качестве примера на рис. П.2 показано поведение величины F (ρ (r)) - μ (r) для димера углерода в зависимости от числа итераций с использованием дискретного вычисления GGA-потенциала. Рис. П.2. Поведение F (ρ(r)) - μ(r) для GGA-приближения в процессе итераций Fig. A.2. Behavior of F (ρ(r)) - μ(r) for the GGA approach in the iteration process Рис. П.3. Зависимость энергии связи атом в димере С2 от числа итераций Fig. A.3. Dependence of the binding energy on numbers of iterations for the C2 dimer Из рис. П.2 видно, что итерационная процедура сходится быстро - уже на третьей итерации величина F (ρ(r)) - μ(r) близка к нулю. Данный подход был опробован на моделировании димеров B2, C2, N2, O2 [11; 12]. Полученные результаты хорошо согласовались с экспериментальными данными и расчетами по стандартному методу КШ. Одно из главных преимуществ БО подхода заключается в высокой скорости вычислений, которая проявляется в быстрой сходимости решения вариационной задачи. На рис. П.3 приведен график величины энергии связи в димере С2 в зависимости от числа итераций. Из этого рисунка видно, что равновесная энергия связи достигается уже на четвертой итерации. Аналогичные расчеты, проведенные с использованием пакета ELK [13] методом КШ, демонстрируют сходимость энергии только на сороковой итерации. То есть, можно говорить, что наш подход увеличивает скорость расчетов на порядок.

About the authors

Victor G. Zavodinsky

Institute of Applied Mathematics of the Russian Academy of Sciences

Email: vzavod@mail.ru
Ph.D, Dr. Sci. (Phys.-Math.), Professor; leader-researcher; Editorial Board member for the Computational nanotechnology Khabarovsk, Russian Federation

Olga A. Gorkusha

Institute of Applied Mathematics of the Russian Academy of Sciences

Email: o_garok@rambler.ru
Cand. Sci. (Phys.-Math.); senior researcher at the Khabarovsk Department of Institute of Applied Mathematics Khabarovsk, Russian Federation

References

Supplementary files