CC BY
34
1.3. ПРИЛОЖЕНИЕ БЕЗОРБИТАЛЬНОГО ПОДХОДА К МОДЕЛИРОВАНИЮ МНОГОАТОМНЫХ СИСТЕМ С РАЗЛИЧНЫМИ НАПРАВЛЕНИЯМИ МЕЖАТОМНЫХ СВЯЗЕЙ
Заводинский Виктор Григорьевич, доктор физ.-мат. наук, профессор, и. о. директора. Институт материаловедения Хабаровского научного центра Дальневосточного отделения Российской Академии Наук. E-mail: vzavod@mail.ru
Горкуша Ольга Александровна, канд. физ.-мат. наук, доцент, старший научный сотрудник. Институт прикладной математики ДВО РАН. E-mail: o_garok@rambler.ru
Аннотация: На примере трехатомных кластеров Al3, Si3 и С 3 показано, что безорбитальный вариант теории функционала плотности может быть использован для нахождения равновесных конфигураций многоатомных систем как с металлической, так и с ковалентной связью. Получены равновесные межатомные расстояния, энергии связи и углы между связями в хорошем согласии с известными данными.
Ключевые слова: Моделирование, функционал плотности, безорбитальный подход, тримеры, ковалент-ные связи.
1.3. APPLICATION OF ORBITAL-FREE APPROACH TO SIMULATION OF MULTI ATOMIC SYSTEMS WITH VARIOUS DIRECTIONS OF INTERATOMIC BONDS
Zavodinsky Victor Griroryevch, the Doctor of Science in Physics and Mathematics, professor. Institute for Material Studies, Khabarovsk Scientific Center, Far Eastern Branch of the Russian Academy of Sciences. E-mail: vzavod@mail.ru
Gorkusha Olga Alexandrovna, PhD in Physics and Mathematics, Institute of Applied Mathematics, Far Eastern Branch of the Russian Academy of Sciences. E-mail: o_garok@rambler.ru
Abstract: On the example of the three-atomic clusters Al3, Si3 and C3 it is shown that the orbital-free version of the density functional theory may be used for finding equilibrium configurations of multi atomic systems with both the metallic and covalent bonding. The equilibrium interatomic distances, interbonding angles and binding energies are found in good accordance with known data.
Index terms: Modeling, density functional, orbital free approach, trimers, covalent bonding.
1. ВВЕДЕНИЕ
В недавние годы появилась серия работ ряда авторов (например [1-7]), посвященных созданию так называемой безорбитальной (БО) версии теории функционала плотности (ТФП), отличающейся от версии Кона-Шэма (КШ) [8] тем, что в ней полностью отсутствуют волновые функции (орбитали). БО подход является последовательным развитием идеи Хоэнбер-га-Кона [9] о том, что основное состояние квантовой системы может быть полностью описано с помощью электронной плотности. В своих недавних статьях [10] мы описали основы нашего подхода к развитию БО метода, которые заключаются построении функционалов многоатомных систем с использованием функционалов плотности одиночных атомов,
известных из расчетов методом КШ. В этих работах, как и в работах других авторов (например [11]) достигнуты определенные успехи в моделировании атомных димеров. Однако существенным камнем преткновения на пути дальнейшего развития БО метода является тот факт, что электронная плотность одиночного (изолированного) атома является сферической, т.е. безорбитальный атом представляет собой шар, а шары, как известно, с неизбежностью формируют структуры с плотной упаковкой. Например, три одинаковых атома обязаны образовать равносторонний треугольник с углами в 60 градусов. В то же время, известно, что, скажем, в случае кремния три атома образуют равнобедренный треугольник с основным углом около 80 градусов [12], атомы
Заводинский В.
углерода выстраиваются в линеиную цепочку [13], а атомы алюминия действительно ведут себя как шары - формируют правильный равносторонний треугольник [14].
Настоящая работа представляет собой попытку развить методику, которая позволяла бы достаточно адекватно описывать геометрию межатомных связей в многоатомных системах в рамках БО подхода.
2. ОБЩЕЕ ОПИСАНИЕ БЕЗОРБИТАЛЬНОГО ПОДХОДА
Напомним, что согласно ТФП электронная энергия основного состояния любой квантовой системы (£е/) может быть найдена путем минимизации некоего функционала, зависящего только от полной электронной плотности данной системы р. Ограничиваясь для простоты случаем отсутствия спиновой поляризации, можно записать:
Е"[р] = Еш[р]+Еех[р] + Ес[р]+Ен[р] -\Уеа(г)р(г)йът, (1)
где \/ех1- внешний потенциал, Екп - кинетическая энергия, Еех - обменная энергия, Ес -корреляционная энергия, Ен - энергия Харт-ри, равная
Eh [р]= 11
If Р( r )P(r ) j 3
2J |r - r|
d rd r
1-1 " г-О
При этом величина полной энергии Е дается интегралом
Еш =|Ее[р(гV . (2)
Минимизация выражения (1) при условии |р(г= N означает решение уравнения
ЗЕе[Р1-, = 0. (3)
Ф
Здесь ц - параметр Лагранжа, имеющий в нашем случае смысл химического потенциала электрона. Его величина находится в процессе решения уравнения (3).
Вводя для удобства обозначение
F[p] = ^ -
6р
и, мы получаем уравнение
F[P] - ~KJr) + <(r) + цш(р) + цех(р) + цс(р) -¡и = 0, (4) где <r) = j^V, ^ = ^ ^ = ^
Ис(Р)
|r
SEC[ Р] 8р .
др
Основную проблему создает неизвестный член Цш(р), связанный с кинетической энергией, выражения для членов Цех(р) и цс(р) достаточно хорошо разработаны в различных приближениях, потенциал <р(r) может быть с хорошей точностью вычислен либо с помощью Фурье-преобразований, либо из уравнения Пуассона (в данной работе мы использовали решение уравнения Пуассона), а в качестве внешнего потенциала Vext(r), как правило, выступает суммарный потенциал атомных ядер или остовов, при псевдопотенциальном подходе. Как и ранее, мы будем вести рассмотрение именно в псевдопотенциальном приближении, используя псевдопотенциалы для различных s-, p- орбитальных состояний. (Детали описаны в работах [10].)
3. ОТ ДИМЕРОВ К ТРИМЕРАМ:
РОЛЬ КВАНТОВЫХ ОГРАНИЧЕНИЙ
Поскольку конструирование псевдопотенциалов сопровождается нахождением равновесных псевдоволновых функций, то всегда можно вычислить его полную равновесную электронную плотность p(r). Если бы мы знали вид функции Ццп(Р), мы могли бы вычислить энер-
I— ч/
гию E kin и найти полную энергию одиночного атома по формуле (2). Однако, во-первых, нет никаких оснований полагать, что существуют некие универсальные функции для кинетических функционалов и энергий, а во-вторых, проблема заключается не в том, чтобы найти энергию одиночного атома (это можно сделать более традиционными способами), а в том, чтобы вычислить энергию взаимодействия атомов друг с другом, а также определить соответствующую этому взаимодействию геометрию многоатомной системы.
3.1. Димеры. Для простоты рассмотрим димеры, состоящие из атомов одного типа. Простейшим приближением для электронной плотности такого димера pdim является сумма плотностей составляющих его атомов:
Pdmir) = Pat (r - RA ) + Pat (r - RB ) '
(5)
где [А и [в - координаты точек, в которых находятся атомы А и В с плотностями ра1. То-
гда энергия связи (на один атом) нашлась бы по формуле
Еь = \(Ет - 2Ва1), (6)
Eat= EГ =j Eel [pat (r )]d V,
где
Edim =J Ed [päm(r )]d 3Г + R-R
(7)
и - положительные заряды атомных остовов, равные по абсолютной величине заряду валентных электронов.
В качестве тестовых элементов возьмем А1, Б1, и С - из соображений, указанных выше: тримеры этих элементов имеют принципиально разные конфигурации. В качестве функции ¡лш(р), единой для б- и р-состояний
и универсальной для всех трех элементов примем функцию
^ш(р) = 09Р3 - 15p2
(8)
В качестве генератора псевдопотенциалов и псевдоволновых функций здесь и далее воспользуемся пакетом РИ!98рр [15], широко применяющимся в расчетах методом КШ. Обменные и корреляционные потенциалы будем вычислять в приближении локальной плотности [16, 17]. Размер кубической ячейки I, в которую помещались атомы, равнялся 52 а.е. (1 а.е. длины равна 0.529 А). Для интегрирования ячейка разбивалась на 128x18x128 элементарных субъячеек, что обеспечивало шаг интегрирования = 0.406 а.е. Результаты расчетов сопоставлялись с опубликованными данными.
Таблица 1.
Равновесные расстояния С и энергии связи Еь (абсолютные значения, на атом) для Б|2, А12 и С2 в сравнении с известными расчетными данными.
Димер Источник d, Ä Eb, eV
Si2 Наш метод 2.1 2.0
Другие данные 2.23a 2.21b 1.97a 1.599b
Alz Наш метод 2.4 1.2
Другие данные 2.46c 2.95e 1.0d 1.23e
С 2 Наш метод 1.3 4.8
Другие данные 1.247-1.367f 1.316g 4.7f 3.5g
Примечание: а [12], ь [19], с [20], " [14], е [21], ' [22], 8 [23]
В таблице 1 приведены значения межатомных расстояний и величин энергии связи для димеров А12, Б12, и С2 в сравнении с известными расчетными данными. Видно, что согласие достаточно хорошее. Особенно, если учесть,
что в литературе имеются данные, существенно отличающиеся друг от друга.
3.2. Тримеры. Для выработки методики выявления направленности сил межатомных взаимодействий обратимся к причинам возникновения таковой направленности в стандартном квантово-механическом рассмотрении с использованием волновых функций и электронных энергетических состояний. Например, в работе [18] указывается, что малые (содержащие несколько атомов) кластеры простых металлов имеют плотноупако-ванные структуры, а материалы с энергетической щелью Е8ар между занятыми и незанятыми состояниями (полупроводники типа кремния) характеризуются ковалентными направленными связями, углы между которыми определяются размером кластера и типом атомов. При этом указывается, что структура кластера Б13 с углом а около 80° определяется эффектом Яна-Теллера, обусловленного именно наличием энергетической щели Е8ар. Т.е. разница структур полупроводниковых и металлических малых кластеров связана с тем, что в первом случае имеются локализованные волновые функции, обеспечивающие направленность межатомных связей, а во втором - волновые функции рассредоточены почти однородно по пространству и направленность связей отсутствует.
В нашем случае волновых функций нет, электронных состояний тоже нет, а, следовательно, и об энергетической щели говорить не приходится. В безорбитальном подходе мы имеем дело только с электронной плотностью, которая единственно и определяет всю энергетику и структуру нашей многоатомной системы. И вот тут надо обратить внимание на некую тонкость описания кван-тово-механической системы с помощью волновых функций и электронных состояний. Помимо уравнений Шредингера (или Кона-Шэма), из которых эти функции и состояния выводятся, мы привычно оперируем принципом Паули, указывающим, что в одном квантовом состоянии может находиться лишь два электрона (без учета спина). В нашем случае это можно перефразировать так: ковалентная
Заводинский В.
связь формируется двумя электронами, общая волновая функция которых локализована в пространстве между двумя ближайшими атомами. При этом очевидно, что количество электронов, отвечающих за эту связь, не изменяется при изменении расстояния между атомами (если, конечно, связь вообще не разрывается и электронная структура не перестраивается полностью). В случае металлической связи состояния проводимости находятся близко друг от друга и при изменении атомной геометрии электроны могут «перетекать» из одного состояния в другое.
На языке электронной плотности вышесказанное можно переформулировать следующим образом: при изменении расстояния между атомами с ковалентной связью интеграл плотности между этими атомами, п^, должен сохраняться; в случае металлической связи величина интеграла может быть произвольной. Разумеется, сразу возникает вопрос: по какой области следует вести интегрирование? И как быть с промежуточными случаями, когда взаимодействуют атомы разных типов? Над этими вопросами следует думать, а пока мы попытаемся выяснить: имеется ли в нашем подходе принципиально разумное зерно, способное хотя бы качественно объяснить разницу в строении и энергетике ковалентных и металлических систем на примере гомогенных кластеров Al3, Si3 и С3?
Рис. 1. Схема расположения атомов в тримере. Пунктирными линиями показана область, по которой ведется интегрирование электронной плотности для определения числа электронов, вовлеченных в ковалентную связь; а - угол между связями одинаковой длины ^
Очевидно, что область интегрирования должна быть достаточно локальной и в то же время она должна давать информацию о количестве электронов, включенных в ковалентную связь. В настоящей работе мы использовали области пространства, имеющие
вид пластинки, находящейся между двумя соседними атомами (рис. 1) и ориентированной перпендикулярно к плоскости, в которой лежит треугольник тримера. Толщина пластинки составляла 2АЬ, где АЬ - расстояние между ближайшими точками, по которым ведется интегрирование. Поскольку при изменении конфигурации тримера количество точек интегрирования в пластинке может изменяться, то находилась величина интеграла, отнесенная к одной точке.
По аналогии с димерами найдем величины электронной плотности тримеров р^т:
Ргт (г) = ра1 (Г - Я.) + ра1 (Г - яв) + ра, (г - Яс) , (9 )
где ЯА, Яв и Яс - координаты точек, в которых находятся атомы А, В и С с плотностями р at. Энергия связи (на один атом) равна
ЕЬ з (ЕМш 3Еа1 V
где
Е,ПШ =| Ее1 [ргш (г Ж V +
Ял -Я„\ \Ял -Яс
%в^с
(10)
(11)
Примем за эталон межатомную связь в ди-мере. Тогда при изменении конфигурации кластера (и межатомных расстояний) величина Р = пы/п'Пт (которую условно можно назвать «силой» связи) не должна превышать единицы. На рисунке 2 приведены расчетные величины Р для кластеров А13, Si3 и С3 как функции угла между межатомными связями в случае, когда ограничения на эти величины отсутствуют. (Для каждого угла мы находили такие величины межатомных расстояний, которые соответствовали минимуму полной энергии кластера.) Видно, что во всех случаях величина Р равна примерно единице при а=180° и возрастает с приближением угла а к 60°. Максимальное значение (Р=1.40) наблюдается у углерода, межатомные расстояния у которого самые малые. Межатомные расстояния у алюминия и кремния близки друг к другу, поэтому неудивительно, что и «силы» связи для тримеров А13 и Si3 ведут себя примерно одинаково.
+
+
Яв - Яс
Interbonding angle a, degree
Рис. 2. «Сила» связи в тримерах Al 3, Si 3 и С3 как
функция угла между межатомными связями в случае, когда ограничения на межатомную плотность электронов отсутствует.
На рисунке 3 (кривые A) представлены результаты расчетов энергии связи для триме-ров Al3, Si3 и С3 без введения ограничения на величину «силы» связи. Из этих кривых видно, что во всех трех случаях максимальная (по абсолютной величине) энергия связи соответствует треугольным кластерам, «силы» связи в которых существенно превышают соответствующие величины, характерные для линейных цепочек. Для алюминия этот результат выглядит естественным (и согласуется с данными других расчетов и экспериментом), поскольку его состояния имеют металлический, не локализованный характер, но для кластеров Si3 и С3, обладающих ковалентными связями, необходимо ввести оговоренные выше ограничения на величины «сил» связи P. Найденные с учетом этого условия (P=1) зависимости энергии связи от угла между связями в кластерах Si3 и С3 приведены на этом же рисунке (кривые B). Мы видим, что атомы углерода, как им и подобает, стремятся образовывать линейные цепочки, в то время как для кремния ни линейная цепочка, ни равносторонний треугольник не являются энергетически выгодными: атомы кремния предпочитают образовывать равнобедренный треугольник с углом а, равным 80 градусов.
i—i—i—■—i—'—i—i—i—i—i—■—i— 60 00 100 120 140 160 ISO
Interbonding angle a, degree
Рис. 3. Расчетная зависимость энергии связи (на атом) для тримеров А13, Б13 и С 3. А - величины, полученные без ограничений на электронную плотность в межатомных связях, В - кривые, рассчитанные при условии постоянства интеграла плотности между соседними атомами.
В таблице 2 приведены равновесные значения межатомных расстояний, углов а и величин энергии связи для тримеров А13, Б13, и С3 (найденные при условии Р=1) в сравнении с известными данными. Видно, что согласие хорошее - в том числе и для углов между атомными связями. Таким образом, мы показали, что безорбитальный подход способен правильно описывать направленность межатомных связей в атомных кластерах, так же как величины межатомных расстояний и энергии связи.
Таблица 2.
Равновесные расстояния С, углы а и энергии связи Еь (абсолютные значения, на атом) для Б13, А13 и С3 в сравнении с известными расчетными данными.
Тример Источник Угол а, град d, Ä Eb, eV
Наш метод =80 2.1 3.7
Sis Другие данные 77.8a 78.10е 79.6е 2.26b 2.177c 2.51b 2.93d
Наш метод 60 2.3 2.4
Al з Другие дан- 60' 60h 2.51' 2.55h 1.74g
ные 1.96'
Наш метод 180 1.2 7.0
Сз Другие данные 180i 180k 1.29' 1.3k 1.316l 6.8i 5.01
Примечание: a [12], b [19], c [27], d [24], e [25], ' [21], g [14], h [26], [13], ' [13], k [22], 1 [23].
Заводинский В.
4. ЗАКЛЮЧЕНИЕ
Проведенный нами анализ поведения межатомной электронной плотности показал, что введение принципа ограничения на величины плотности (вытекающего из принципа Паули) позволяет правильно описывать угловые зависимости направления межатомных связей в многоатомных кластерах в рамках безорбитального варианта теории функционала плотности. В частности, удается показать, что для кластера Al3 выгодным является конфигурация в виде равностороннего треугольника, тример Si3 характеризуется равнобедренным треугольником с углами 80 и 50 градусов, а три атома углерода выстраиваются в прямолинейную цепочку. Вычисленные равновесные межатомные расстояний и величины энергии связи хорошо согласуются с известными расчетными данными.
Список литературы:
1. Wang Y.A., Carter E.A. Orbital-free kinetic-energy density functional theory. In: Progress in Theoretical Chemistry and Physics. Kluwer, Dordrecht. 2000, 117 p.
2. Huajie Chen, Aihui Zhou. Orbital-Free Density Functional Theory for Molecular Structure Calculations. Numerical Mathematics: Theory, Methods and Applications, 2008, 1, 1-28.
3. Baojing Zhou, Ligneres V.L., Carter E.A. Improving the orbital-free density functional theory description of covalent materials. Journal Chemical Physics, 2005, 122, 044103-044113.
4. Hung L., Carter E.A. Accurate Simulations of Metals at the Mesoscale: Explicit Treatment of 1 Million Atoms with Quantum Mechanics. Chemical Physics Letters, 2009, 475, 163-170.
5. Karasiev V.V., Trickey S.B. Issues and challenges in orbitalfree density functional calculations. Computational Physics Communications, 2012, 183, 2519-2527.
6. Karasiev V.V., Chakraborty D., Shukruto O.A., Trickey S.B. Nonempirical generalized gradient approximation free-energy functional for orbital-free simulations. Physical Review B, 88, 161108-161113(R).
7. Wesolowski T.A. Approximating the kinetic energy functional Ts[p]: lessons from four-electron systems. Molecular Physics, 2005, 103, 1165-1167.
8. Kohn W., Sham J.L. Self-Consistent Equations including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133-A1138.
9. Hohenbeg H., Kohn W. Inhomogeneous Electron Gas, Physical Review, 1964, 136, B864-B871.
10. В.Г. Заводинский, О.А. Горкуша. ФТТ, 56, 2253 (2014); V.G. Zavodinsky, O.A. Gorkusha. A practical way to develop the orbital-free density functional calculations. Physical Science International Journal, 2014, 4(6), 880-891; В.Г. Заводинский, О.А. Горкуша. На пути к моделированию больших наноси-стем на атомном уровне. Computational nanotechnology, 2014, 1, 11-16; V.G. Zavodinsky O.A. Gorkusha. A new OrbitalFree Approach for Density Functional Modeling of Large Mole-
cules and Nanoparticles. Modeling and Numerical Simulation of Material Science, 2015, 5, 39-47.
11. Junchao Xia, Chen Huang, Ilgyou Shin, Carter E.A. Can orbital-free density functional theory simulate molecules? The Journal of Chemical Physics, 2012, 136, 084102(13).
12. Raghavachari K., Logovinsky V. Structure and bonding in small silicon clusters. Phys. Rev. Lett. 1985, 55, 2853-2856.
13. Van Orden A., Saykally R.J. Small carbon clusters: spectroscopy, structure, and energetics. Chemical Review, 1998, 98, 2313-2357.
14. Feng-Chuan Chuang, Wang C.Z., Ho K.H. Structure of neutral aluminum clusters Aln (2<n<23): Genetic algorithm tight-binding calculations. Phys. Rev. B, 2006 ,73, 125431(7).
15. Fuchs M., Scheffler M. Ab initio pseudopotentials for electronic structure calculations of poly-atomic systems using density-functional theory, Computational Physics Communications, 1999), 119, 67-98.
16. Perdew J.P., Zunger A. Self-interaction correction to density functional approximation for many-electron systems, Physical Review B, 1981, 23, 5048-5079.
17. Ceperley D.M., Alder B.J. Ground state of the electron gas by a stochastic method, Physical Review Letters, 1980, 45. 566569.
18. Tomanek D., Schluter M.A. Structure and bonding of small semiconductor clusters. Phys. Rev. B, 1987 36, 1208-1217.
19. Mukhtarov A.P., Normurodov A.B., Sulaymonov N.T., Umarova F.T. Charge States of Bare Silicon Clusters up to Si8 by Non-Conventional Tight-Binding Method. Journal of nano- and electronic physics, 2015, 7, 01012(7).
20. Nayak S.K., Khanna S.N., Jena P.J. Evolution of bonding in AlnN clusters: A transition from nonmetallic to metallic character. Physical Review B, 1998, 57, 3787-3790.
21. Matrinez A., Vela A. Stability of charged aluminum clusters. Physical Review B, 1994, 49, 17464(4).
22. Karton A., Tarnopolsky A., Martin J.M.L. Atomization energies of the carbon clusters Cn (n=2-10) revisited by means of W4 theory as well as density functional, Gn, and CBS methods. International Journal of Interface between Chemistry and Physics, 2009, 107, 977-1003.
23. Mahdi Afshar, Mahboobeh Babaei, Amir Hossein Kordbacheh. First principles study on structural and magnetic properties of small and pure carbon clusters (Cn, n = 2-12) Journal of Theoretical and Applied Physics, 2014, 8, 103-108.
24. McCarthy M.C., Thaddeus P. Rotational spectrum and structure of Si3. Physical Review Letters, 2003, 90, 213003(4).
25. Liu B., Lu Z.Y., Pan B., Wang C.Z., Ho K. M., Shvartsburg A.A., Jarrold M.F. Ionization of medium-sized silicon clusters and the geometries of the cations. Journal of Chemical Physics, 1998, 109, 9401-9409.
26. Raghavachari K., Rohlfing C.M. Bonding and stabilities of small silicon clusters: A theoretical study of Si7-Si10. Journal of Chemical Physics, 1988, 89, 2219-2234.
