В окончательной редакции 19 августа 1998 г.) Стандартная схема ab initio расчета структуры кристаллов с использованием нелокальных псевдопотенциалов модифицирована для использования в криволинейных координатах. Найден метод решения уравнения Пуассона для кулоновского потенциала в искривленном пространстве в k-представлении. На примере расчетов кристаллов диэлектриков со структурой NaCl показано, что применение искривленного пространства позволяет очень значительно уменьшить необходимый размер базиса.
В последнее десятилетие в связи с развитием тех- В данной работе предложена модификация данного нологии и возрастанием интереса к описанию сложных подхода, не использующая подход молекулярной динамикро- и макро-структур появилась необходимость в ab мики, так как для небольшого числа базисных функций initio расчетах сложных структур, насчитывающих десят- данный подход не имеет преимуществ перед классичеки, сотни и тысячи атомов в ячейке. Включение динамики ским прямым методом диагонализации гамильтониана.
в такие расчеты [1Ц4], основанные на методе молекулярной динамики, предложенном в [5], а также нахождение отклика системы на различные возмущения (расчеты 1. Схема расчета фононных спектров и т. д.) потребовали наличия простых ab initio расчетных схем, легко модифицируемых приВ алгоритме расчета использовалась стандартная схеменительно к требуемым задачам. К сожалению, имеюма расчета по методу псевдопотенциала, изложенная щиеся методы являлись или достаточно сложными для в [8]. В качестве псевдопотенциалов использовались модификации (методы, основанные на решении уравнеnorm-conserving псевдопотенциалы, рассчитанные и прония Шредингера для МТ-сфер: ЛМТО (LMTO) метод, табулированные в [9]. Обменно-корреляционные эфFull Potential метод), или требовали больших машинных фекты учитывались в рамках формализма функционала ресурсов (в частности, метод псевдопотенциала для опиплотности (LDA-аппроксимация) с использованием апсания атомов с большими псевдопотенциалами требовал проксимации из [10,11].
использования базиса из более чем 1000 плоских волн Все физические величины (волновые функции, локальдаже для простых структур из нескольких атомов).
ные и нелокальные части потенциалов, гамильтонианы, В последние 5Ц7 лет появились новые эффективные методы расчета. Метод ultra-soft псевдопотенциала, пред- электронная плотность) рассчитывались с помощью ба ложенный в [6], в значительной мере уменьшил недо- зисных функций вида k(r) |k) = (r) eik(r), r статки классического метода псевдопотенциала ценой которые, как и обычные плоские волны, являются ортонекоторого усложнения схемы расчетов.
нормированными [1] Другой новый перспективный метод, предложенный в работе [7], является гибридом метода ЛМТО и метода 1 (r) -k)(r) псевдопотенциала и позволяет избавиться от многих (k|k ) = d3r ei(k недостатков обоих методов. r В то же время в [1,2] (см. также [3,4]) в рамках метода молекулярной динамики КараЦПаринелло был предложен новый подход к решению уравнений КонаЦШема в = d3(ei(k -k)) =k,k, рамках псевдопотенциальной схемы. Подход заключался во введении в качестве базисных функций вместо обыч ных плоских волн вида |k eikr Ч искривленных k (r) -1/rk здесь =gi j, где метрический тензор gi j=r j.
r i (r) плоских волн (ИПВ) вида |k) eik(r). Легко Переходя при вычислении матричных элементов гаr мильтониана, а также потенциалов и электронной плотпоказать, что такие волны образуют ортонормированный ности в данном базисе от пространства r к пространству базисный набор, и матричные элементы гамильтониана в этом базисе трансформируются в матричные элементы,, легко видеть, что выражения для них становятся рассчитываемые в базисе обычных плоских волн, но уже такими же, что и для случая плоских волн, но в пространв специальном искривленном пространстве, описывае- стве, за исключением форм лапласиана в операторе мом взаимно однозначным отображением = (r). кинетической энергии.
5 242 А.С. Федоров Уравнения КонаЦШема в обратном пространстве имеVlps(g) = d3r(Vlps((r))e-ig(r)). (5) ют вид H(g, g )(g ) =(g), Матричные элементы нелокальной части псевдопотенгде Ч собственный вектор (волновая функция), циала имеют вид n который для каждого раскладывается по базису k+g(r) (k|Vnl|k ) = k((r))|Ylm(r) Vlm(r) 1 (r) i,l,m ((r)) = k (g) ei(k+g)(r), (1) r g Ylm(r)|k ((r)) e-iRi (k-k ). (6) Ч собственное значение уравнения Шредингера с Электронная плотность выражается следующим волновым вектором k и спином в зоне n. В (1) образом:
= {k,, n}, k Ч волновой вектор в первой зоне Бриллюэна, g Ч вектор обратной решетки, n Ч номер зоны. ((r)) = ((r))((r)), (7) d Матричные элементы гамильтониана H = - +V 2m drв этом базисе имеют следующий вид:
(g) = k+g k+g+g, (8),g q|H|q = q|T +Vnl|q + d3(Vl((r))ei(q -q)), здесь и далее Ч функция заполнения уровня.
Вычисление кулоновского потенциала в базисе q|T|q = d3 [(qi - iAi)gi j(q j (r) eik(r) является нетривиальной задачей, 2m r так как если в ФплоскомФ пространстве лапласиан в k-представлении имеет вид диагональной матрицы + iAj)] ei(q -q), (2) (k + g)2g,g, то в искривленном пространстве (r) лапласиан имеет недиагональные компоненты (см.
где q|T |q описывает лапласиан в базисе правую часть (2)). Система уравнений Пуассона (r) записывается eik(r) (см. [1]) и где gi j обозначает r тензор, обратный к метрическому тензору gi j, а (k||k )(k |Vcoul|k ) =-4(k|((r))|k ) 1 r масштабный потенциал Ai = log возникает k при дифференцировании базисной функции. Полный -4(k - k ). (9) потенциал разбивается на нелокальную часть q|Vnl|q и на локальный потенциал Vl(q - q ) d3r(Vl((r)) Матрица лапласиана в любом полном базисе имеет минимальное собственное значение 0, равное 0, которое ei(q -q)(r)). Локальный потенциал в свою очередь соответствует решению-константе C((r)) = const из разбивается на кулоновский, обменно-корреляционный вырожденного подпространства собственных векторов в и локальную часть псевдопотенциала пространстве r. Выбор k k +g = 0 позволяет наинизVl((r)) Vcoul((r)) + Vxc(((r))) + Vlps((r)). шему собственному вектору лапласиана 0 описывать вектор C. Но если для диагональной матрицы исключеМатричные элементы обменно-корреляционного ние этого собственного значения является тривиальной потенциала и локальной части псевдопотенциала задачей, то для недиагональной матрицы это не так.
(k|V ((r))|k ) рассчитывались путем фурье-преобразоДля решения системы уравнений использовался факт вания, исходя из значений на равномерной сетке (r).
ортогональности нашего искомого решения Vcoul((r)) и Обменно-корреляционный потенциал Vxc в приближевектора из вырожденного подпространства нии локальной плотности и обменно-корреляционная энергия xc = d3r(xc()((r))), где xc() Ч (r) C((r)) = C(g) eig(r).
r плотность обменно-корреляционной энергии была взята g из [10]:
xc Это является следствием электронейтральности криVxc((r)) =, (3) сталла (r) Vxc(g) = d3r(Vxc((r))e-ig(r)), (4) d3r(Vcoul((r)) + Vpseudo((r))) = Физика твердого тела, 1999, том 41, вып. Использование криволинейных координат в ab initio расчетах диэлектриков... и выбора нормировки где каждый вклад описывался гауссовой функцией двух параметров: и d3r(Vcoul((r))) = 0, (10) r = - ( - 0)ie-i(-i )2.
i i откуда следует ортогональность величин Vcoul и constC. Здесь i характеризует радиус действия искривления Далее, выбирая в (9) k = 0 для получения минимально- пространства от соответствующего иона с координатой го собственного значения лапласиана и ортогонализуя 0, а i описывает амплитуду возмущения пространства i вектор правой части к вектору C(g) ( C), можем от данного иона.
использовать алгоритмы решения вырожденной системы Минимизация полной энергии E = E(i, i) проводиуравнений с правой частью, ортогональной вырожден- лась путем прямой минимизации с помощью квазиньюному собственному вектору (имеющей единственное тоновского метода минимизации.
решение), и перейти от решения системы (9) к (11).
2. Результаты (k|( - 0)|k )(k |Vcoul|k = 0) k Полученные результаты представлены на рис. 1Ц4 и в табл. 1 и 2.
= -4(k|((r))|k = 0) -4(k - 0). (11) На рис. 1 показана зависимость энергии связи E кристалла MgO в зависимости от объема элементарной В данных расчетах из-за неполноты базиса величина ячейки V и от числа базисных функций как для случая 0 3 10-4, что может служить косвенным фактором плоских волн (Xg {i, i} 0), так и для случая качества базиса, связанного с отображением (r).
ИПВ волн (Xg = 0). Видно, что введение искривленного Эффектом искривленного пространства, несмотря на пространства для базиса эквивалентно очень значительвыполнение соотношения (10), являются и следующие ному увеличению числа базисных функций в ФплоскомФ условия:
пространстве. Также видна очень медленная сходимость энергии связи с увеличением базиса в стандартном Vcoul(g = 0) = 0, Vloc(g = 0) = 0.
подходе. Здесь же показано экспериментальное значение равновесного объема элементарной ячейки V 0 и энергии Полная энергия связи кристалла имела вид связи из [17].
bZ На рис. 2 показано наилучшее отображение = (r) E/N = Emad/N + - Vcoul(g)(g) +Vcoul(g = 0) (для N = 339) для кристалла MgO в плоскости (001).
g Видно сгущение координатной сетки (а вместе с этим и (r) увеличение амплитудного множителя + Vps(g = 0) + d3r((xc() базисных r функций) в окрестности ионов O, где электронная плот ность максимальна, и разрежение координатной сетки в окрестности ионов Mg. Параметры кривизны для - xc())(r)) +, (12) N исследуемых кристаллов приведены в табл. 2.
где Emad Ч электростатическая энергия Маделунга, Z Ч Таблица 1. Равновесные параметры полный заряд ячейки, F Ч энергия Ферми, определяе1 bZ мая из условия Z =, Ч некулоновская N bZ 4ZeMgO 3.74 4.18 4.21 [13] 1.61 1.32 1.53 [15] = lim Ves(q) +. (13) q0 q2 BaO 4.03 5.19 5.54 [14] 8.46 0.45 0.74 [16] s NaCl 5.00 5.49 5.63 [15] 0.530 0.238 0.245 [15] PbS 5.76 5.77 5.92 [16] 0.96 0.763 0.62 [16] Отметим, что в гамильтониане не использовался фиктивный член, связанный с энергией деформации пространства, как в [1Ц4]. Таблица 2. Параметры искривления Отображение (r) выбиралось из предположения о (r) корреляции с электронной плотность и с целью Кристалл i i Размер базиса r описания отображения (r) меньшим числом парамеMgO 0.500/ - 0.503 0.991/0.688 тров. BaO 0.476/ - 0.191 0.332/0.937 Из этих соображений отображение выбиралось в виде NaCl 0.172/ - 0.352 0.110/0.177 PbS 0.101/ - 0.220 0.992/0.497 аддитивной суммы вклада от каждого атома в решетке, 5 Физика твердого тела, 1999, том 41, вып. 244 А.С. Федоров Рис. 1. Зависимость энергии связи E от объема элементарной ячейки V от числа плоских волн и волн ИПВ для MgO. На рис. 3 показано распределение электронной плот- сываются положительными значениями i (соответствуности (координата Z) в плоскости (001) (координаты ющими ФсгущениюФ пространства в местах с большой X, Y ) для MgO в пространстве (r) в пределах ближай- электронной плотностью), а катионы Ч отрицательными ших соседей. (соответствующими разрежению). На рис. 4 в аналогичных координатах X, Y показано (r) распределение якобиана перехода (координата r Z) в том же пространстве. Видна корреляция между электронной плотностью и якобианом, вызывающая значительное сглаживание в пространстве (r) электронной плотности, что приводит к уменьшению эффективного размера гамильтониана в пространстве (r) и позволяет резко уменьшить необходимый размер базиса в этом пространстве. В табл. 1 показаны равновесные параметры элементарной ячейки и величины модуля всестороннего сжатия в исследуемых кристаллах при использовании плоского и искривленного пространств и экспериментальные данные. Видны значительная недооценка параметра ячейки и погрешность в модуле сжатия для обычного подхода для всех исследуемых кристаллов и значительно лучшее согласие с экспериментом при использовании искривленных координатах. В табл. 2 показаны параметры кривизны i, i и соответствующий размер базиса для всех исследуемых Рис. 2. Наилучшее отображение = (r) (N = 339) для кристаллов. Видно, что в данных кристаллах анионы опи- MgO. Физика твердого тела, 1999, том 41, вып. Использование криволинейных координат в ab initio расчетах диэлектриков... Рис. 3. Распределение электронной плотности для MgO в пространстве (r). (r) Рис. 4. Распределение для MgO в пространстве (r). r Физика твердого тела, 1999, том 41, вып. 246 А.С. Федоров Таким образом, данный метод представляется перспективным в рамках псевдопотенциального подхода, так как позволяет значительно (в 10раз иболее) уменьшить размер базиса. При этом алгоритм расчета в пространстве меняется незначительно по сравнению со стандартным. Из полученных данных видно значительное понижение энергии в расчетах с искривленным базисом, что является следствием неточного описания псевдоволновых функций плоскими волнами вблизи атомов, где псевдопотенциал наиболее существен и где ФизмельчениеФ координатной сетки благодаря искривлению пространства позволяет более точно описать поведение псевдоволновых функций. Также видно существенно лучшее согласие с экспериментом при определении равновесных параметров ячейки. При этом желательны дополнительные исследования, направленные на более быстрое нахождение оптимальных характеристик искривленного пространства. Автор выражает глубокую благодарность В.И. Зиненко за плодотворные дискуссии. Работа поддержана Российским фондом фундаментальных исследований (проект 96-02-16542). Список литературы [1] F. Gygi. Europhys. Lett. 19, 617 (1992). [2] F. Gygi. Phys. Rev. B48, 11692 (1993). [3] D.R. Hamann. Phys. Rev. B51, 7337 (1995). [4] D.R. Hamann. Phys. Rev. B51, 9508 (1995). [5] R. Car, M. Parinello. Phys. Rev. Lett. 55, 2471 (1985). [6] D. Vanderbilt. Phys. Rev. B41, 7892 (1990). [7] P.E. Blchl. Phys. Rev. B50, 17953 (1994). [8] J. Ihm, A. Zunger, M.L. Cohen. J. Phys. C12, 4409 (1979). [9] G.B. Bachelet, D.R. Hamann, M. Schlter. Phys. Rev. B26, 4299 (1982). [10] J. Perdew, A. Zunger. Phys. Rev. B23, 5048 (1981). [11] D.M. Ceperley, V.J. Alder. Phys. Rev. B18, 3126 (1978). [12] Е.Г. Бровман, Ю.М. Каган. УФН 112, T3, 369 (1974). [13] G. Kalpana, B. Palanivel, M. Rajagopalan. Phys. Rev. B52, (1995). [14] S.A. Chang, C.W. Tompson, E. Gurnen, L.D. Muhlestein. J. Phys. Chem. Sol. 36, 769 (1975).