Energetics and Elastic Properties of Large Nano-objects: Orbital-free Approach on the Basis of the Density Functional Theory

Full Text

1. Введение Современные нанотехнологии нуждаются в мощных инструментах, способных предсказывать свойства систем, содержащих сотни тысяч и даже миллионы атомов. Традиционные квантово-механические подходы, такие как версия Кона-Шэма (КШ) теории функционала плотности (ТФП), не дают возможности оперировать большими количествами атомов, их лимит не превышает тысячу атомов даже при использовании псевдопотенциалов. Методы эмпирических потенциалов позволяют работать с большими системами, однако они не обеспечивают достоверности и надежности результатов. С другой стороны, возможности повышения быстродействия компьютеров приближаются к своему физическому пределу, поэтому трудно ожидать решения проблемы этим способом. Необходим новый метод моделирования, который соединял бы квантово-механическую точность с возможностью исследования огромного количества атомов. Идея такого метода появилась еще в 1964 г., когда Хоэнберг и Кон показали [1], что энергия основного состояния любой квантовой системы полностью определяется ее электронной плотностью. В этой же работе они объявили, что существует некий универсальный функционал Ee[ρ], минимизация которого приводит к равновесной электронной плотности ρ и к равновесной полной электронной энергии Ee. Функционал Ee[ρ] может быть записан в следующем виде: (1) где ε(ρ) - плотность полной электронной энергии; V(r) - внешний потенциал; φ(r) - электростатический потенциал Хартри; εex-c(ρ) и εkin(ρ) - плотности обменно-корреляционной и кинетической энергий. Минимизация выражения (1) при условии сохранения полного числа электронов Nel(∫ρ(r)drel) означает решение следующего уравнения: (2) где μ(r) - электронный химический потенциал, который в конечной системе зависит от координат. Из уравнения (2) следует: (3) Чтобы его решить, нам необходимо знать все члены, входящие в ε(ρ) (см. (1)), а также потенциал μ(r). Рассмотрим эти члены. Внешний потенциал V(r) представляет собой сумму полноэлектронных атомных потенциалов: где Zi - заряд ядра атома с номером i. Потенциал Хартри может быть вычислен с помощью преобразований Фурье или из уравнения Пуассона. Существуют также достаточно реалистичные аппроксимации для обменно-корреляционной энергии εex-c(ρ) (например приближение локальной плотности LDA [2; 3], которое мы используем в данной работе). Единственная серьезная проблема - это кинетическая энергия εkin(ρ). В свое время были сделана попытка применить для вычисления кинетической энергии приближение Томаса-Ферми [4; 5], основанное на теории свободных электронов. Эта попытка оказалась неудачной, все молекулы оказывались нестабильными. Поправка Вейцзекера [6], учитывающая неоднородность электронной плотности тоже не исправила ситуацию коренным образом. Кон и Шэм представили компромиссный вариант [7]. Они предложили находить кинетическую энергию Ekin из решения некоего одноэлектронного уравнения, известного нынче как уравнение Кона-Шэма. Этот подход сделался очень популярным, на его основе разработано множество эффективных вычислительных программ и решено много важных и интересных задач, однако, как было замечено выше, его возможности уже практически исчерпаны. Безорбитальный (БО) подход является развитием теории функционала плотности и представляет собой альтернативу методу КШ. Достоинство этого подхода очевидно: работая только с электронной плотностью вместо многочисленных волновых функций, он позволяет резко увеличить скорость вычислений и включить в рассмотрение огромное число атомов. Первые попытки развить БО метод начались более двадцати лет назад и касались моделирования жидких металлов в приближении «желе» [8]. Затем появились работы других исследователей (см. например статьи и обзоры [9-12]) в приложении к простым молекулам и твердым телам. Все эти работы использовали псевдопотенциалы и в большинстве своем пытались применять приближения Томаса-Ферми и Вейцзекера в различных комбинациях. Однако эти попытки не имели большого успеха и не получили широкого распространения. Нам кажется, что главная причина их неудач заключается в том, что разработчики не отходили от идеи существования универсального функционала энергии (в том числе кинетической) для любых систем. Недавно было показано [13; 14], что утверждение Хоэнберга и Кона о существовании универсального функционала, приводящего к минимуму энергии, не было строго доказано, и более того - неверно. Поэтому является актуальным поиск и использование специальных видов кинетических потенциалов для различных атомных систем. В недавних работах [15; 15а] мы описали наш безорбитальный полноэлектронный (БО-ПЭ) подход и продемонстрировали, что он адекватно описывает гомоатомные и гетероатомные димеры - их энергии и длины связей. В данной работе мы распространяем этот подход на описание больших наносистем, содержащих тысячи и десятки тысяч атомов. 2. Основные особенности подхода Рассмотрим для простоты некую систему, состоящую из N одинаковых атомов, и введем для этой системы функционал F(r): (4) где Наша задача найти такую плотность ρ(r), которая удовлетворяла бы уравнению F(ρ) - μ(ρ) = 0. (5) Поиск такой плотности мы будем осуществлять с помощью итерационной процедуры: (6) где K - численный параметр, контролирующий сходимость процедуры. В качестве начальной плотности возьмем сумму атомных плотностей: Здесь имеется две проблемы. Во-первых, мы не знаем кинетический потенциал μkin(ρ), входящий в F(ρ). Во-вторых, нам не известен химический потенциал μ(r). Вторая проблема решается просто: мы находим равновесное значение F(ρ) и получаем μ(r) из уравнения F(ρ) - μ(r) = 0. Однако для этого нам надо знать μkin(ρ), и эта проблема является принципиальной. В наших работах [15, 16] было предложено искать μkin(ρ) для двухатомных систем, используя кинетические потенциалы одиночных атомов μakin(ρ), в следующем виде: μkin(ρ) = μakin(ρ)f, (7) где f - некая функция межатомного расстояния d. Для многоатомных систем мы предлагаем использовать подобное приближение, однако в этом случае функция f должна зависеть от среднего расстояния между ближайшими соседями dav. Одноатомный потенциал μakin(ρ), мы можем легко найти с помощью стандартных вычислений методом КШ следующим образом. Для равновесного состояния одиночного атома мы можем для простоты положить μa = 0 и записать (8) откуда следует (9) Далее, поскольку мы знаем равновесную величину ρa(r), мы легко делаем замену переменных и находим μakin(ρ). Решив уравнение (5) и найдя тем самым плотность ρ(r) для многоатомной системы, мы можем вычислить полную энергию Etot: Etot = EH + EC + Eex-c + Ek + Er, (10) где EH = 1/2∫φ(r)ρ(r)dr - энергия Хартри; EC = ∫V(r)ρ(r)dr - кулоновская энергия; Eex-c = ∫εex-c(ρ)dr - обменно-корреляционная энергия; Ekin = ∫εkin(ρ)dr - кинетическая энергия; Er - энергия отталкивания атомных ядер, Для нахождения равновесной атомной геометрии системы мы вычисляли силы на атомах через производные от полной энергии по координатам атомов и cдвигали атомы до тех пор, пока энергия не достигала минимума. После этого мы находили равновесную плотность, следуя итерационной процедуре (6) и вычисляли окончательную равновесную полную энергию по формуле (10). 3. Энергия когезии Мы выполнили вычисления энергии когезии для углерода, алюминия, кремния и титана. Углеродная и кремниевая системы были взяты как кубические фрагменты кристалла со структурой алмаза, алюминий - как фрагмент кристалла с гранецентрированной кубической ячейкой, а титан был для простоты исследован в кубической объемно-центрированной фазе. Для нахождения равновесных плотностей одиночных атомов мы использовали пакет FHI98pp [17], который обычно применяется для конструирования псевдопотенциалов, однако в процессе нахождения псевдопотенциала он вычисляет полноэлектронную структуру атома, из которой легко строится полная атомная плотность. Опуская технические детали вычислений, остановимся на фундаментальном моменте: на нахождении кинетических потенциалов. Как говорилось выше, наш подход требует для этого задания неких функций f(dav). В работе [18] указывалось, что энергия системы зависит не только от среднего расстояния между ближайшими соседями dav, но и от среднего числа этих соседей Nn. Мы приняли это указание во внимание и использовали следующие функции для исследованных нами материалов, приводящие не только к достоверным значениям энергии когезии, но и к разумным величинам модуля упругости: Среднее число ближайших соседей Nn зависит от структуры системы и, как правило, автоматически возрастает с ростом системы. Численные параметры, входящие в функции, не зависят от количества атомов, то есть, приведенные формулы справедливы для систем любой величины. Все системы были изучены в кубических ячейках таких размеров, которые сохраняли достаточное пустое пространство вокруг атомной системы; для расчетов использовался обыкновенный персональный компьютер. Скорость вычислений зависела от плотности трехмерной сетки, на которую разделялась рабочая ячейка для проведения процедур интегрирования. Для каждого материала мы использовали такую сетку, которая обеспечивала бы достаточную точность вычислений в сочетании с разумными затратами времени. Для углерода была взята сетка 130 × 130 × 130, для кремния - 100 × 100 × 100, для алюминия 50 × 50 × 50, для титана - 150 × 150 × 150. (При вычислении модулей упругости плотность сеток была увеличена для обеспечения нужной точности.) Мы исследовали зависимость величины энергии когезии от размера системы, то есть от числа атомов. Результаты вычислений представлены на рис. 1. Легко видеть, что во всех случаях вычисленные величины энергии возрастают с ростом системы в соответствии с известными данными для малых кластеров углерода [19; 20], алюминия [21; 22], титана [18; 23] и кремния [24; 25]. Рис. 1. Зависимость величины энергии когезии от размера системы. Белые кружки демонстрируют результаты, полученные без оптимизации плотности, черные - с оптимизированной плотностью. Прерывистые линии показывают экспериментальные значения величин энергии для кристаллов Fig. 1. Dependence of cohesive energy on the system size. Open circles demonstrate results obtained without density optimization; solid circles correspond to results with density optimization. Dashed lines show experimental energy values for bulk crystals В случае углерода и кремния энергия когезии приближается к экспериментальной величине для кристалла при числе атомов около 3000. Титан показывает хорошую сходимость уже при 1300 атомах - после оптимизации плотности. Алюминий демонстрирует самую медленную сходимость: лишь при n = 15 000-20 000 энергия когезии становится близкой к экспериментальному значению в кристалле. Отметим также, что энергия титановой системы существенно изменяется в процессе оптимизации электронной плотности, что косвенно свидетельствует о серьезной перестройке электронной структуры, характерной для взаимодействия атомов с d-электронами. 4. Объемный модуль упругости Главный недостаток моделирования многоатомных систем с помощью безорбитального подхода заключается в том, что с его помощью невозможно получить информацию об электронной структуре системы, о ее электрических и оптических свойствах. Однако он может быть эффективно использован для изучения механических свойств. Для этого мы должны убедиться, что с его помощью возможно не только нахождение правильных значений энергии, но и корректное описание сил межатомных взаимодействий, а последнее наиболее отчетливо проявляется в упругих свойствах. Обычно такое тестирование заключается в расчетах объемного модуля упругости B. Для нахождения модуля упругости мы вычисляли полную энергию кубических частиц Cn, Sin, Aln и Tin, представляющих собой фрагменты соответствующих кристаллов, изменяя величину параметра решетки кристалла и фиксируя положение всех граничных атомов. Найдя таким образом зависимость энергии от объема частицы, мы аппроксимировали эту зависимость вблизи минимума параболой и находили величину В по формуле где Ea и va - энергия и объем, приходящиеся на один атом. На рис. 2 представлены зависимости величин объемного модуля упругости для исследованных материалов от размера частиц. Рис. 2. Зависимость величины объемного модуля упругости В для частиц Cn, Sin, Aln и Tin, от числа атомов п Fig. 2. Dependence of bulk modulus B for Cn, Sin, Aln and Tin particles on number of atoms п Как видно из рис. 2, для всех исследованных материалов модуль упругости растет с увеличением размера частицы, приближаясь к значениям Bсr, характерным для массивных кристаллов. Для углерода Bсr = 445 ГПа, для кремния - 103 ГПа, для алюминия - 72 ГПа, для титана - 120 ГПа. Существует ряд публикаций (например [26-30]), в которых приводятся данные, свидетельствующие о том, что величины модулей упругости наночастиц значительно превышают соответствующие значения тех же модулей у массивных кристаллов. Однако в этом нет никаких противоречий с нашими результатами. Дело в том, что структура наночастиц, о которых идет речь в таких работах, существенно отличается от структуры массивного кристалла, атомы связаны между собой энергетически более выгодно, что и ведет к увеличению упругих сил. В работе [31], где изучалась зависимость упругих свойств углерода, кремния и германия, взятых в виде фрагментов алмазных решеток, получены результаты, аналогичные нашим: модуль упругости растет с увеличением размера фрагмента. Упомянем еще, что один из авторов данной работы (В.Г. Заводинский) проводил ранее расчеты упругих свойств (модуля Юнга) частиц кремния. В том случае, когда исследовались частицы, у которых атомная структура оптимизировалась [26] и окружение атомов не имело ничего общего с их окружением в кристалле, модуль упругости значительно превышал его значение в кристалле. Когда же в качестве объекта исследования был взят фрагмент кристалла кремния [33], модуль упругости оказался существенно ниже. Таким образом, мы можем утверждать, что наш метод вполне адекватно описывает упругие свойства больших нано-объектов и пригоден для исследования их механических характеристик. 5. Быстродействие и сходимость Важными характеристиками любого метода моделирования являются скорость вычислений и максимально возможный размер изучаемых систем. Как правило, эти величины связаны между собой, и в идеале исследователь хотел бы иметь метод, позволяющий максимально быстро производить расчеты для огромных систем. Все это в особенности касается безорбитального подхода, который предназначен для работы именно с большими системами. На рис. 3 мы приводим в качестве примера зависимость вычислительного времени, затрачиваемого на одну итерацию в процессе оптимизации структуры алюминия. Во-первых, мы видим, что эта зависимость линейна. Это очень важный факт, поскольку у всех существующих методов моделирования аналогичная зависимость носит нелинейный характер: квадратичный, или выше. Во-вторых, скорость вычислений огромная: на одну итерацию для частиц из 20 тысяч атомов алюминия затрачивается лишь 20 минут. Рис. 3. Зависимость времени выполнения одной итерации при оптимизации структуры наносистем Aln от числа атомов п Fig. 3. Dependence of the one-iteration structure optimization time for the Aln от nanosystem on number of atoms n Выше речь шла о скорости выполнения одной итерации. Возникает вопрос: а сколько же итераций требуется для полной оптимизации системы? На рис. 4 представлена зависимость величины энергии когезии частицы углерода (1000 атомов) как функция числа итераций. Как мы видим, процесс оптимизации сходится достаточно быстро, за двадцать итераций. Рис. 4. Зависимость энергии когезии частицы углерода (1000 атомов) от числа итераций Fig. 4. Dependence of cohesive energy of the carbon particles (1000 atoms) on number of iterations Как уже говорилось, скорость вычислений зависит от плотности сетки разбиения рабочей ячейки: чем ниже плотность, тем выше скорость. Однако уменьшение плотности сетки обычно приводит к потерям в точности вычислений. Этот эффект мы продемонстрировали на примере титана на рис. 5, из которого можно видеть, что использование сеток 125 × 125 × 125 и 150 × 150 × 150 дает примерно одинаковые результаты, а сетки 100 × 100 × 100 и тем более 50 × 50 × 50 непригодны для расчетов. Рис. 5. Зависимость величины энергии когезии частиц Tin от плотности сетки разбиения рабочей ячейки Fig. 5. Dependence of cohesive energy of Tin particles on the used mesh 6. Заключение На примерах углерода, кремния, алюминия и титана мы продемонстрировали, что наш безорбитальный полноэлектронный подход позволяет описать энергетику и упругие свойства больших нано-объектов, содержащих тысячи и десятки тысяч атомов. Наш метод обладает высоким быстродействием, а вычислительное время зависит от размера системы линейно. В данной работе мы представили вычисления только для гомоатомных систем, однако метод пригоден и для систем, содержащих атомы различных типов. В наших работах [15; 16] описано, как распространить данный подход на гетероатомные системы, используя весовые функции, которые конструируют кинетические потенциалы в пространстве между атомами разных типов.

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.), Profes-sor; leader-researcher at the Khabarovsk Department 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 re-searcher at the Khabarovsk Department Khabarovsk, Russian Federation

References

Supplementary files