FULL-ELECTRON ORBITAL-FREE MODELING METHOD FOR ATOMIC SYSTEMS: THE FIRST STEP

Abstract

We studied an opportunity to develop a full-potential orbital-free method for modeling of multi-atomic systems using results of Kohn-Sham calculations for single atoms. We have obtained equilibrium bond lengths and binding energies for dimers Li2, Be2, B2, C2, N2, O2, F2, Na2, Mg2, Al2, Si2, P2, S2 & Cl2, as well as for C3, C24 and C60 systems in good accordance to other theoretical and experimental data.

Full Text

Введение Современные технологии - нанотехнологии, биотехно- логии, производство лекарств и так далее, нуждаются в мощ- ных инструментах для прогнозирования свойств систем, со- εex-c и εkin - обменно-корреляционная и кинетическая энер- гии (на электрон). Минимизация уравнения (1) означает ре- шение следующего уравнения:  держащих сотни тысяч и миллионы атомов. Традиционные квантово-механические подходы, такие как теория Хартри- image F     V r   r   ex-c  kin   0, (2) Фока, методы квантовой химии и теория функционала плот- ности (ТФП), не дают возможности работать с большим ко- личеством атомов; их пределы не превышают нескольких тысяч атомов даже с использованием псевдопотенциалов при условии ∫ρ (r)dr = N, где N - количество электронов в си- стеме, а величины    ex-c     kin  image и суперкомпьютеров. Методы эмпирических потенциалов и Монте-Карло позволяют работать с большими системами, ex-c image  и kin   но они не обеспечивают нужную достоверность результатов. С другой стороны, возможности для увеличения компьютер- ной скорости близки к физическому пределу, и, следователь- но, бесполезно надеяться на решение проблемы этим путем. Таким образом, назрела настоятельная потребность в новом методе моделирования, который сочетал бы квантово-ме- ханическую точность с возможностью работать с огромным количеством атомов. Идея такого метода возникла еще в 1964 г., когда Хоэн- берг и Кон сформулировали теорему [1] о том, что энергия основного состояния любой квантовой системы полностью определяется ее электронной плотностью. В той же рабо- те они заявили, что существует некоторый универсальный являются так называемыми кинетическим и обменно-корре- ляционным потенциалами. Потенциал Хартри может быть вычислен с использовани- ем преобразований Фурье или уравнений Пуассона, внеш- ний потенциал обычно состоит из атомных потенциалов или псевдопотенциалов. Существуют также некоторые реа- листичные аппроксимации для обменно-корреляционного потенциала (например, [2-4]). Единственная реальная про- блема заключается в кинетическом потенциале (или кинети- ческой энергии). Были попытки использовать аппроксима- цию Томаса и Ферми (TF) [5; 6] на основе подхода свободных электронов: функционал E [ρ], минимизация которого должна приводить TF 3 image 2 2 3 5 3 TF image 1 2 2 3 1 3 к равновесной электронной плотности ρ и полной энергии Etot. Функционал E [ρ] может быть записан в следующем виде: image kin   3  10 image  ; kin   3   . 2 image 2 E   dr V r r dr  1 r r dr  ex-c dr kin dr , (1) Этот функционал оказался абсолютно неадекватным (все молекулы оказывались нестабильными), а поправка Вайцзе- кера (Weizsacker (W) [7]) где ε (ρ) - плотность энергии; V (r) - внешний потенциал; φ - электростатический электронный потенциал Хартри, kin  W   1 8 image r  2 r  dr r image    r  r dr; не решила проблему (величины энергии связи получались неверные). 80 Других серьезных нововведений в этой области не было, поэтому Кон и Шэм предложили компромиссный подход [8]. Они стали вычислять кинетическую энергию Ekin путем реше- ния некоторого одноэлектронного уравнения, гамильтониан которого зависел только от электронной плотности: полуэмпирических. В настоящей работе мы описываем по- пытку развития безорбитального полноэлектронного подхо- да, работающего без псевдопотенциалов. Основные положения Прежде всего, рассмотрим одиночный атом, равновес- image  1  r   V 2 i eff r i r   i r i r ; ная полная электронная плотность которого может быть лег- ко вычислена методом КШ. В соответствии с (2) мы можем E 1 image kin  Veff r   V r   r   ex-c ; i r i r dr  i i r Veff i r dr , kin написать уравнение для нахождения одноатомного кинети- ческого потенциала μ(1) : Z 2 i i i где image image 2 kin r    r   ex-c r  , image r image kin kin где Z - полный ядерный заряд атома. Если мы знаем μ(1) (r) и ρ (r), то мы легко можем найти μ(1) (ρ). Исходя из резуль- r    i r  ; i ψj (r) - волновые функции, или орбитали; εi - энерг(и3я) i-го состояния. Уравнение (3) было названо уравнением Кона-Шэма (КШ). Уравнение КШ стало широко известным; на его осно- ве создано много эффективных вычислительных программ и было решено много задач моделирования многоатомных систем; однако, как было сказано выше, его возможности практически исчерпаны. Безорбитальный (БО) подход как версия теории функционала плотности мог бы составить аль- тернативу методу КШ. Он является последовательным разви- тием идеи Хоэнберга-Кона [1] о том, что основное состояние квантовой системы можно полностью описать с помощью электронной плотности. Преимущество этого подхода оче- видно: работа только с электронной плотностью вместо мно- татов полноэлектронных КШ расчетов (с помощью пакета FHI98pp [22]), мы сконструировали кинетические потенциа- лы для атомов B, C, N, O. Они показаны на рис. 1. image Кинетический потенциал, усл. ед. 10 000 O 8000 N 6000 C 4000 B 2000 0 0 50 100 150 200 250 300 image N O B C 14 image image image image 12 10 8 6 4 2 гочисленных волновых функций позволяет резко увеличить скорость вычислений и включить в расчет огромное количе- ство атомов. 0 0,00 0,25 0,50 0,75 1,00 1,25 1,50 1,75 2,00 Плотность, усл. ед. Первые попытки разработки БО-метода моделирования начались около 20 лет назад. Это были имитации жидких металлов в приближения к желе [9; 10]. Затем последова- ли работы других исследователей (см., например, обзоры и оригинальные статьи [11-14]), которые были применены к некоторым простым молекулам и твердому веществу. Все они были основаны на использовании специальных псев- допотенциалов и большинство из них пытались исполь- зовать аппроксимации TF и W для кинетической энергии в некоторых комбинациях. Однако эти попытки не имели большого успеха. Нам кажется, что главная причина их не- достаточной неэффективности в том, что они пытались (и пытаются) использовать некий универсальный функ- ционал кинетической энергии для всех систем. Однако не- давно было показано [15, 16], что идея Хоэнберга-Кона о су- ществовании универсального функционала плотности, веду- щего к энергетическому минимуму, не была строго доказана. В наших недавних работах [17-21] мы описали безорби- тальный псевдопотенциальный подход для моделирования наносистем, содержащих атомы с s-, p- и d-электронами. Ключевым моментом подхода было нахождение кинетиче- ской энергии с использованием некоторых функций, специ- ально подобранных для каждого типа атомов. Этот подход был испытан на кластерах, содержащих атомы C, Al, Si, O, Ti и Cu, и продемонстрировал хорошее согласие с методом Кон-Шэма и экспериментальными данными. Однако по- строение псевдопотенциалов является довольно неодно- значной операцией и приближает этот подход к категории Рис. 1. Кинетические потенциалы, рассчитанные для равновесных одноатомных полноэлектронных плотностей. kin На верхней панели показан общий вид μ(1) (ρ). kin Нижняя панель показывает μ(1) (ρ) для малых плотностей Можно видеть, что построенные кривые очень сильно отличаются друг от друга, подтверждая вывод работ [15; 16], что не существует универсальный кинетический потенциал для различных квантовых систем. Однако возникает вопрос: если у нас есть одноатомные кинетические потенциалы для каждого вида атомов, почему мы не можем использовать их для многоатомных систем? В чем разница между одиноч- ным атомом и, например, димером? Рассмотрим димер, состоящий из двух атомов бора с использованием полноэлектронного КШ-кода Elk [23]. Равновесная плотность электронов этого димера показана на рис. 2. image Плотность, усл. ед. 14 12 10 8 6 4 I II I 2 0 11 12 13 14 15 16 17 18 Межатомное расстояние, ат. ед. Рис. 2. Равновесная электронная плотность димера бора. Межатомное расстояние равно 3 ат. ед. (1,59 Ǻ) kin Исходя из найденной равновесной плотности, мы можем найти двухатомный кинетический потенциал μ(1) (r): центрированный на атомном ядре. Такие пики соответствуют локализованным остовным состояниям, которые не участву- ют в межатомных взаимодействиях и обычно рассматрива- 2 r   ZB ZB       ются как «замороженные» (так делается, например, в пакете image kin image r  R r  R  r ex-c r , Elk [23]). Такое разделение плотности электронов помогает из- 1 2 где ZB = 5; R1, R2 - координаты первого и второго атомов ди- мера. Из рис. 2 ясно, что в случае димера существуют два раз- личных типа областей. В первом из них (I) плотность растет от нуля на большом расстоянии от атома до максимально- го значения при атомном ядре. Вторым типом области (II) является пространство между атомами. Наши расчеты по- казали (рис. 3), что в области I двухатомный кинетический потенциал ведет себя идентично одноатомному потенциалу. Это естественно, так как электрон в этой области находит- ся далеко от атома 2 и взаимодействует только с атомом 1. В области II электрон взаимодействует с обоими атомами, кинетический потенциал уменьшается, и уровень уменьше- ния зависит от межатомного расстояния: чем меньше рассто- яние, тем меньше кинетический потенциал. image 5 Кинетический потенциал, усл. ед. 4 бежать операций с интенсивными резкими пиками и постро- ить реалистичные компьютерные коды. В нашей работе мы также следуем этой технике и разделяем атомную плотность на остовную и валентную части. В частности, мы считаем, что атом бора имеет два остовных электрона и три валентных, и распределение их плотностей показано на рис. 4. image 0,5 Плотность, усл. ед. 0,4 0,3 0,2 0,1 3 2 1 0 0,0 I II 0,1 0,2 0,3 0,0 Плотность, усл. ед. 0,0 0,0 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 Межатомное расстояние, ат. ед. Рис. 4. Плотность электронов одного атома бора: сплошная линия - полная атомная плотность; штриховая линия - плотность электронов 2s22p1, пунктирная линия - валентная плотность, Рис. 3. Поведение кинетического потенциала в различных областях димера бора Очевидно, что прямой способ получения правильных результатов для безорбитальных расчетов заключается в вы- используемая в нашем БО-подходе Функционал F12 произвольного гомоатомического диме- ра может быть записан следующим образом: полнении КШ-расчетов для различных атомных позиций, по- image image F  Z  Z   r       (2)  , (4) иске кинетического потенциала в каждой точке пространства 12 r  R r  R 12 ex-c 12 kon 12 и использовании его для БО-расчетов. Однако такой способ очень дорог и не имеет смысла. Нам кажется, оптимальный вариант - найти некоторые закономерности из КШ-вычисле- 1 2 где Z - атомный ядерный заряд; core core val  r ний на уровне димера и распространить их на многоатомные системы в БО-подходе. Простейший способ состоит в том, 12  1  2  21 ; 12 r    image 12 r  r dr; чтобы ввести некоторый средний кинетический потенциал, величина μ (ρ ) рассчитывается с использованием хорошо действующий в многоатомной системе и дающий правиль- ную связывающую энергию и атомную конфигурацию. Для димеров мы сконструировали следующее выражение для двухатомного кинетического потенциала: ex-c 12 известных подходов (в нашем случае это приближение ло- кальной плотности). Поскольку мы начинали с равновесных состояний для единичных атомов 1 и 2, мы можем записать image    (2) r   (1) r  1,0   expN2  , Z a r   (1) a    a ; (5) d kin kin  val  r  R1 image 1 kin 1 ex-c 1 где d - длина димера; α и β - константы подгонки; Nval - число 2 Z a r   (1) a    a , (6) валентных электронов в атоме. Параметры α и β контроли- руют значения энергии связывания и длины связи. Если они приспособлены для одной системы (например, для димера где ρa image r  R21 image a kin 2 ex-c 2 1 и ρ2 - равновесные электронные плотности электронов атомов 1 и 2; ϕa и ϕa - электростатические потенциалы, обра- B2) и удовлетворительно работают для других димеров, этот 1 2 этап нашего подхода можно считать успешным, и мы можем перейти к моделированию более сложных систем. зованные этими плотностями. Подставляя (5) и (6) в (4), получаем F   r   a r   a r       1 a   Методология 12 12 1 2 ex-c 12 kin 1  (1) a   2     a    a . (7) Полноэлектронное рассмотрение атомов обладает неко- kin 2 kon 12 ex-c 1 ex-c 2 торым особенностям, которые делают его довольно сложным. Одна из них - резкий интенсивный пик плотности электронов, Наша цель - найти такую плотность ρ12, которая превра- щала бы выражение (7) в ноль. 12 С учетом того, что остовные плотности не изменяются при межатомных взаимодействиях, мы можем записать ите- рационное уравнение для валентной плотности ρval: системы: Li-F и Na-Cl. Результаты для энергий связи (на атом) показаны на рис. 5, длины димеров собраны в табл. 1. image 10 val i   val i    K F i   val i   , 9 1 12 12 1 iter 12 1 12 1 2 Энергия связи, эВ/атом 8 3 где Kiter - итерационный параметр, управляющий сходимо- 7 4 стью процедуры. Начальный шаг (i = 0) означает, что плот- 6 12 ность ρval (0) представляет собой сумму равновесных атом- 5 ных валентных плотностей: 4 val (0) = (ρval)a + (ρval)a. 3 ρ12 1 2 2 Отметим, что полная энергия димера (Etot)12 является сум- мой энергии отталкивания ядер (Erep)12, кулоновской энергии (EC)12, энергии Хартри (EH)12, обменно-корреляционная энер- гии (Eex-c)12 и кинетической энергии (Ekin)12: image image E   Z  Z ; rep 12 R  R 1 2 1 0 -1 Li Be B C N O F Na Mg Al Si P S Cl Типы атомов Рис. 5. Энергия связи для исследованных димеров: 1 - эксперимент [24]; 2 - полноэлектронные КШ-вычислений с использованием кода Elk [23]; 3 - результаты КШ-вычислений с использованием псевдопотенциального кода FHI96md [25]; E    Z Z  image image image image    r dr; 4 - наши БО-расчеты C 12   r  R  12   r  R  1 2 Из рис. 5 ясно, что наши энергетические результаты хо- 1 рошо соответствуют экспериментальным данным, во многих image EH 12  2 12 12 r dr; Eex-c 12  ex-c 12 dr; Ekin 12  kin 12 dr. Обсуждение результатов Напомним, что мы выбрали димер бора (B2) в каче- стве тестового объекта. Мы нашли для него энергию связи Eb = 1,8 эВ и равновесную длину d = 1,59 Ǻ, взяв параметры α = 1,08 и β = 0,13. Затем мы использовали эти значения параметров для всех димеров двух рядов Периодической случаях они демонстрируют даже лучшую точность по срав- нению с КШ-вычислениями. Длины димеров, найденные по БО-методу (см. табл. 1) также удовлетворительно согласу- ются с экспериментальными величинами, хотя это согласие несколько хуже, чем в случае КШ-расчетов. В общем, мы можем сказать, что с учетом того факта, что подгоночные параметры α и β были подобраны только для одного типа димеров (B2), наш БО-подход демонстрирует довольно хорошую способность описывать взаимодействие атомов. Однако мы должны были проверить подход на при- мерах более крупных систем. С этой целью мы провели рас- четы для углеродных систем: C3 (равносторонний треуголь- ник), C24 (фрагмент нанотрубки) и C60 (фуллерен). Результаты собраны в табл. 2. Равновесные длины связи d(Ǻ) для изученных димеров Таблица 1 Метод Li2 Be2 B2 C2 N2 O2 F2 Na2 Mg2 Al2 Si2 P2 S2 Cl2 БО 3,02 2,43 1,59 1,16 1,00 1,16 1,59 3,02 4,02 2,33 2,16 1,91 1,91 1,85 КШ-полный 3,07 2,46 1,63 1,23 1,11 1,18 1,37 3,07 3,89 2,54 2,22 1,83 1,90 1,96 КШ-псевдо 2,92 2,48 1,64 1,27 1,11 1,20 1,40 3,04 3,91 2,62 2,22 1,80 1,91 1,98 Экспримент [24] 3,09 2,47 1,59 1,24 1,10 1,15 1,417 3,07 3,891 2,56 2,32 1,90 1,88 1,99 О б о з н а ч е н и я: БО - наши БО-расчеты; КШ-полный - полноэлектронные вычисления с использованием сода Elk [23], КШ-псевдо-результа- ты КШ-вычислений с использованием псевдопотенциального кода FHI96md [25]. Величины энергии связи Eb C3 C24 C60 Eb, еВ Наши БО-расчеты 6,7 8,1 11,0 КШ-полный (Elk [23]) 6,8 7,8 10,5 для углеродных систем Таблица 2 Можно видеть, что в этом случае наш БО-метод также демонстрирует хорошее согласие с результатами, получае- мыми в рамках подхода Кона-Шэма. Выводы В этой работе мы продемонстрировали, что возможно создать полноэлектронный безорбитальный метод для мо- делирования многоатомных наносистем с использовани- ем одноатомных кинетических потенциалов, полученных из расчетов Кона-Шэма. Мы предложили способ построения двухатомного кинетического потенциала для димера бора (B2) и успешно применили этот подход для гомоатомных димеров, состоящих из атомов с s- и p-электронами. Затем на примерах кластеров C3, C24 и C60 мы продемонстрировали принципиальную способность нашего БО-подхода описы- вать многоатомные системы.