CC BY
9
ЧЕБЫШЕВСКИЙ СБОРНИК
Том 19. Выпуск 2
УДК 511.32 DOI 10.22405/2226-8383-2018-19-2-101-110
О повышении точности вычисления потенциала в системе взаимодействующих атомов
Заводинский Виктор Григорьевич — ведущий научный сотрудник Института материаловедения ХНЦ ДВО РАН, институт материаловедения ХНЦ ДВО РАН. e-mail: vzavod@mail.ru
Горкуша Ольга Александровна — старший научный сотрудник Хабаровского отделения института прикладной математики ДВО РАН, Хабаровское отделение института прикладной математики ДВО РАН. e-mail: 684-bmts@rambler.ru
Аннотация
Мы предлагаем высокоточный метод вычисления потенциала для многоатомной системы в прямом пространстве. Отличительная особенность метода состоит в разделении электронной плотности р и потенциала у та две части: р = р0 + р, у = + (р, где р0 — сумма сферических атомных плотностей, "р — результат взаимодействия атомов в многоатомной системе; потенциал f0 порождается плоти остью р0, потенциал (р, порожденный плотностью р, в нашей работе находится путем решения уравнения Пуассона.
Для нахождения граничных условий применяется мультипольное разложение потенциала. Для обеспечения высокой точности мы разделяем расчетное пространство на многогранники Вороного и применяем асимптотические оценки итераций при замене характеристической функции гладкими приближениями. Для численного решения уравнения Пуассона мы используем двух- сеточный метод и Фурье- преобразование на этапе начальной итерации.
Мы получили теоретические оценки точности метода О(ha-1), где h — шаг сетки, а — фиксированное число, большее 1.
Ключевые слова: уравнение Пуассона, электростатический потенциал, многогранникик Вороного, мультипольное разложение, двухсеточный метод.
Библиография: 19 названий. Для цитирования:
В. Г. Заводинский, О. А. Горкуша. О повышении точности вычисления потенциала в системе взаимодействующих атомов // Чебышевский сборник, 2018, т. 19, вып. 2, с. 101-110.
102
B. R 3aBO/piHCKHH, O. A. TopKynia
CHEBYSHEVSKII SBORNIK Vol. 19. No. 2
UDC 511.32 DOI 10.22405/2226-8383-2018-19-2-101-110
On the precision increasing in calculation of potential for the
systems of interactive atoms
Zavodinsky Viktor Grigorievich — leading researcher at the Institute of materials science of the Khabarovsk scientific center of the far eastern branch of the Russian academy of sciences, Institute of materials science of the Khabarovsk scientific center of the far Eastern Branch of the Russian academy of sciences e-mail: vzavod@mail.ru
Gorkusha Olga Aleksandrovna — senior researcher of the Khabarovsk branch of the institute of applied mathematics, far eastern branch of the Russian academy of aciences, Khabarovsk branch of the institute of applied mathematics of the far eastern branch of the Russian academy of sciences. e-mail: 684-bmts@rambler.ru
Abstract
We propose a high precision method of finding of potential for multi-atomic quantum-mechanical tasks in real space. The method is based on dividing of electron density and potential of a multi-atomic system into two parts. The first part of density is found as a sum of spherical atomic densities; the second part is a variation of density generated by interatomic interaction. The first part of potential is formed by the first part of density and may be calculated correctly using simple integrals. The second part of potential is found through a Poisson equation from the second part of density. To provide a high precision we divided a work space into Voronoy's polyhedrons and found the boundary conditions by means of a multi-pole distribution of potentials formed by local densities concentrated in these polyhedrons. Then we used double-grid approach, and fast Fourier transformations as initial functions for iterative solution of the Poisson's equation. We estimated accuracy of the offered method and carried out test calculations which showed that this method gives the accuracy several times better than accuracy of the fast Fourier transformation.
Keywords: Poisson equation, electrostatic potential, Voronoy polihedra, multipole expansion, double-grid method.
Bibliography: 19 titles. For citation:
V. G. Zavodinsky , O. A. Gorkusha, 2018, "On the precision increasing in calculation of potential for the systems of interactive atoms" , Chebyshevskii sbornik, vol. 19, no. 2, pp. 101-110.
1. Введение
Мы рассматриваем задачу вычисления электростатического потенциала <р(г), г € 3, формируемого электронной плотностью р(г) = р(г; К\,..., Кт) систем ы т взаимодействующих атомов, расположенных в точках К\,..., Кт (см. [1, 2, 3]):
*М = /я ^«
где О = {г € Н3|р(г) = 0} — открытая, связная, ограниченная область с гладкой границей.
Существуют разные подходы вычисления интеграла потенциала. Наиболее разработанными являются метод Фурье-преобразования — работы [4, 5, 6, 7], метод мультипольного разложения — работы [8, 9], и метод, опирающийся на решение уравнения Пуассона — работы [10, 11].
Каждый из перечисленных подходов имеет свои недостатки. Метод Фурье-преобразования обеспечивает хорошую точность лишь в обратном пространстве (пространстве волновых функций), метод мультипольного разложения пригоден лишь для случаев, когда потенциал находится в точках, расположенных за пределами существования электронной плотности. Использование третьего метода для многоатомных систем требует больших вычислительных затрат.
В нашей работе мы предлагаем комбинацию этих подходов, а также некоторые другие разработки последних лет.
2. Описание алгоритма
Традиционным способом нахождения электронной плостности в системе взаимодействующих атомов является итерационная процедура. На нулевой итерации электронная плотность атомной системы представляет собой сумму электронных плотностей невзаимодействующих атомов. Поскольку электронные плотности невзаимодействующих атомов сферичны, то нулевая плотность р0(г) есть сумма сферических плотностей центрированных на атомах, расположенных в точках ^ :
т
Р°(Г)= ^Ра8феге (1Г~Щ |). (2)
3 = 1
Плотности ра3феге полагаются известными, вычисленными с помощью подходящего квантово-механического метода, например, в рамках теории функционала плотности [2]. Поскольку эти плотности являются сферическими, то соответствующий электростатический потенциал вычисляется с высокой точностью с помощью одномерного интегрирования (см. (1)):
^рЬегеХ^) = J рарЬ,еге(г)Г ^Г + 2И J раРНеге(г)Г^.
Вследствии этого потенциал многоатомной системы, соответствующий начальной плотности (2), находится по формуле
т
<Р°(г)= Т.^РЬеге(1г - Я3 |). (3)
3 = 1
В процессе итераций плотность системы изменяется — обозначим это изменение через р>(г), и будем полагать, что для некоторого положительного числа А
тах \ р>(г \ < А. геп 1 1
Тогда р(г) = р°(г) + р(г), где р°(г) задается формулой (2). Потенциал, соответствующий р(г) приводится к виду
<р(г) = (р°(г) + <р(г), г е О, (5)
где ^°(г) определяется соотношением (3) .Следуя работе [12], для нахождения (р(г) мы решаем уравнение Пуассона
А(р(г) = -4жр(г), г е О (6)
с граничными условиями
двух- сеточным методом.
¡р(г) = / , dr', г е дÜ (7)
JП \г - г'\
3. Вычисление граничных значений потенциала
Поскольку атомы в изучаемой системе могут быть расположены совершенно произвольно, мы разобьем расчетную область Q на многогранники Вороного v(j) [14]:
т
П = U v{3\
3=1
v(j) = {г е ü\\r - Rj\ = min \r - Rk\}.
l^fc^m,
Таким образом, в каждом многограннике Вороного содержится только один атом. Затем при помощи характеристической функции области v(j)
Xj (r) = { J,
1, г е v(j);
0, в противном случае,
представим плотность 'р(г) в виде суммы плотностей f>j (г) = 'р(г) ■ \j (г). Тогда интеграл (7) запишется следующим образом:
т
ф(Г) = Y1 &(r),
3 = l
4>i (г) = i гЩ; dr'. (8)
Jv(i) \r - r'\
В таком представлении есть существенный недостаток: pj(г) имеет скачок на границе области у(i), чт0 отрицательно влияет на сходимость и устойчивость вычислительного процесса. Один из способов решения этой проблемы — заменить ступенчатую характеристическую функцию ее приближением, обладающим свойствами гладкости. В нашей работе мы использовали алгоритм построения такой функции, предложенный Becke в работе [12]. Здесь мы оценим точности этого подхода. Прежде всего, заметим, что
Xj(г)= П (1 - ®(v(r,i,j))),
l^i^rn i=3
где В — тета- функция Хэвисайда, rq(r,i,j) — гиперболическая координата точки г относительно фокусов в точках Ri и Rj. Затем представим B(t) как предел последовательности
{в к(¿)Ь>=1 возрастающих функций, производные которых по переменной £ образуют 5— последовательность. В работе [12] это представление выглядит так:
в к (I ) к( )
1 + а к ^) 2 , —1, —1;
Рк (г), |*| < 1, (Рк = Р оРо ...Р);
к
1,
Здесь Р(¿) — полином, удовлетворяющий условиям Р(±1) = ±1, Р'(t)|¿=±1 = 0 :
р 1 • *
Полагая
Ъ(г;к)= П (1 — вк(л(.г,г,3))),
г=3
получим гладкое приближение характеристической функции
Ъ (г;к) =
Уз (г; к) '3(1 ;к) = Е т=1 -Рг (г;к)
Ъ (г
хз(г) = ^ -Рз(г
Т>з (г
к), г €{П1,...,Пт};
к) + Щ-1 + £1(а), геди(3);
к) + £1(а), в остальных случаях,
где величина ^(г) равна числу областей, границы которых содержат точку г — это либо 2, либо 3, причем равенство ^(г) = 3 выполняется только в конечном числе точек, и
£1(а) = 0(а2 ), а € (0,1), к ^ то. (9)
Используя полученные асимптотические формулы, приведем интеграл (8) к виду:
& (г) = / р(Г У3'(Г']к) ¿г> + & (г) + е2(г), (10)
.¡иИ) К — I |
где
&(г) =1 • ^Г^(11)
3 .)ди а к — г'| а (г1) ' 1 ;
£2(г)= £1(а) •/ Г^Лг', (12)
и О^ Щ) — 5— окрестность точки Щ.
Обозначим через К(и(3^) радиус минимальной сферы с центром в точке Щ, покрывающей ( ). { € 3| | —
Щ| < Щ(и(3))} ПдО = 0. Тогда величину ■р—ту можно представить в виде сходящегося ряда по системе многочленов Лежандра Рг (сов в) [8, 9, 13, 15, 16]:
1 1 ^ г ^ / „ч 11"' — Щ |
=^ {совв)- *=ы,
с
с в = А(г — И^, — — И^). С учетом этого равенства, запишем интеграл, стоящий в правой части соотношения (10) в виде:
г тгэ(г>,к) (1г, = ^(г,к)+ез(г,.), (13)
„МЯ \г — г \
где
1 С ~ ^ (У - R•|V
&(г,к) = ^ р(г')Тд(г'; к) • £ р,
(14)
ез(г, Л = о(--А , •( П(и^ ,Ь (15)
,Л V\т — К3\—Е(и(^)) \\r-R,\) У 1 ;
В последнем выражении константа А задается соотношением (4). Интегралы (11), (12) оцениваются следующим образом:
, // м 8(дУ Щ ^ Л,
(г)\ ^ 1--Т}( (т • \р(г )\,
3 \г — Rj \ — R( и(з)) г' едиЬ)
£2(г) = £\(а),
где 5 — площадь границы области ь>(^. Из последних двух оценок и формул (10), (13), (15) следует, что вклад в асимптотику интеграла (р^ (г) с главным членом (14) вносит сумма
64(г, а, ]) = е3(г, ]) + £г(а), (16)
в которой оценки £\(а), £з(г, ]) вычисляются по формулам (9), (15). Подставляя полученное выражение для (р^ (г) в (8), приходим к равенству
(р(г) = (р'(г)+ £4(г, а, з), (17)
где индекс ] обозначает ту область граница которой содержит точку г и
&(г) = ^¡р(г, к) (18)
3 = 1
с функцией (р(г, к), заданной в (14).
4. Оценка численного решения уравнения Пуассона
Перейдем к решению уравнения Пуассона (6), (7). Учитывая (17) и (18), представим решение этой задачи в виде
Ср(г) =и(г) + и°(г), (19)
где и — решение задачи (6) с граничным условием
и| ап = (р', (20)
и°
ничным условием и°\эп = £4. Так как решение и° достигает максимума на границе, то согласно (16)
\и°(г)\ « £4(А) + £г(а),
тах ■ , х / х 1 1 1
геда \r — Rj \—R( и(з)) \\r — Rj \
£4(А) = о( тах-—„ А ти•(^),ь ^ж. (21)
Для численного решения задачи (6), (20) мы использовали двух-сеточный метод, описанный в работах [17, 13, 10].
Обозначим чере Qh сетку с шагом h на П, через Г^ множество узлов (ih,jh,kh), не принадлежащие Qh, с условием
дП С {(х,у, z) G Cubeh(ih,jh,kh)},
(ih,jh,kh)erh
где Cubeh(x,y, z) — куб с длиной стороны, равной h, и с центром в точке (х,у, z), при этом хотя бы одна из вершин куба содержится в Qh. Обозначим Qh = Qh U rh. Задаче (6), (20) поставим в соответствие разностную задачу
(Lh(uh))itjik = -4тгр( ih,jh, kh), (ih,jh,kh) G Qh,
uhjyk = 0(ih,jh,kh), (ih,jh, kh) G rh,
где
uhjk = uh( ih,jh,kh),
h/2 _ + Ш iThi hw u'2i-l,2j, 2k 2 u'2i,2j, 2k + u2i+l,2j, 2k . (L (u ))г,3,к = -(Щ2- +
h/2 _ 2~h/2 + h/2 . u'2i,2j-l,2k 2 u'2i,2j,2k + u'2i,2j+l,2k .
+ (h/2)2 +
h/2 _ 2~h/2 + h/2 u'2i,2j,2k-l 2 u'2i,2j,2k + u'2i,2j,2k+l /„„ч
+ -Щ2-,
/2 /2
и величины uii-k, — решения разностной задачи на предыдущей и текущей итерациях.
Известно, что разностная схема (22) сходится и погрешность относительно решения задачи (6), (20) имеет второй порядок аппроксимации. Учитывая (19), (21), получаем погрешность приближенного решения относительно задачи (6), (7):
0(г) = U\^k + 0(h2) + £l(h) + £a(Ä), г = (гh,jh, kh) G Üh.
Согласно (21), £^(А) ^ ^. Положив j = lnA/ ln h, и используя (5), получаем асимптотическую формулу для 0>(г) :
0(г) =р0(г) + ^к + 0(max{h2, h1-1}), г = (ih,jh,kh) G Üh, (23)
в которой потенциал задается в (2).
5. Численная схема решения уравнения Пуассона
Шаг итерационного процесса состоит из трех этапов.
1. Производим одну итерацию на Qh/2, вычисляя u2i'2j2k п0 Ф°РмУле (22). Новое приближенное решение к решению на Qh на этом этапе определяется соотношением
~h = ~h/2 ui,j,k = u2i,2j,2k.
2. Вычисляем параметр сходимости м итерационного процесса:
\\Lh(uh) - 4ттp|Uh
М
\\Lh(uh) - 4жр\у '
где оператор Ьн задается в (22), || ■ || — дискретная норма в П .
3. Если у близко к единице, то итерационный процесс завершен. Иначе мы интерполируем функцию и^ к с П на при помощи оператора интерполяции, точного для многочленов второй степени. И переходим к первому пункту итерационного процесса.
В итерационных вычислениях важную роль играет выбор начального приближения. В нашем случае в качестве начального приближения мы использовали потенциал, вычисленный с помощью БПФ, а именно:
и(г) ^ _с{тх,ту ,тг)_ А(тх^х+ту •у+тг •г) (9л\
и(,) 47Г А2 + (2тгЬ)2 ■ (ш1 +т2 + т2) , 1 ;
г = (х,у, г) е СиЬе(П) .
Здесь с(тх,ту,тг) — коэффициенты Фурье функции 'р(г), X — некоторая константа, значительно меньшая единицы [18], не влияющая на конечный результат.
6. Заключение
1. Данный подход был реализован в работах, описанных в [19]. В них показано, что разделение полной плотности р системы взаимодействующих атомов на начальную плотность ро, состоящую из суммы сферических плотностей, центрированных на отдельных атомах, и некоторую плотность р, которая образуется в результате взаимодействия атомов, обеспечивает более высокую точность по сравнению с вычислением потенциала напрямую от полной плотности.
2. Теоретическая оценка погрешности метода позволяет контролировать параметры численного расчета — амплитуда добавочной плотности должна быть сравнима с величиной h, с 7 е (1, 2).
СПИСОК ЦИТИРОВАННОЙ ЛИТЕРАТУРЫ
1. Ландау Л. Д., Лифшиц Е. М. Квантовая механика. Нерелятивистская теория. -М.:Физматгиз, 1963.- 768 с.
2. Kohn W., Sham J.L. Self-consistent equations including exchange and correlation effects // Phvs. Rev. 1965. Vol. 140, ЛМА. pp. A1133 Л113*.
3. Заводинский В. Г. Компьютерное моделирование наночастиц и наносистем. М.:Физматлит, 2013.- 137 с.
4. Skollermo G. A Fourier method for the Numerical Solution of Poisson Equation// Mathematics of Computation. 1975. Vol. 29, №131. pp. 697-711.
5. Chun- Min Chang, Yihan Shao, Jing Kong. Ewald mesh method for quantum mechanical calculations// J. Chem. Phvs. 2012. Vol. 136, №11 pp. 114112-114112-5.
6. Bobrov V. В., Trigger S.A. The problem of the universal density functional and the density matrix functional theory// Journal of Experimental and Theoretical Physics. 2013. Vol. 116, №4, pp. 635-640.
7. Chelikowskv J.R., Troullier N, Saad Y. Finite-difference-pseudopotential method: Electronic structure calculations without a basis//Phvs. Rev. Lett. 1994. Vol. 72, №8, pp. 1240-1243.
8. L. Kleinman L., Bvlander D. M. Efficacious Form for Model Pseudopotentials // Phvs. Rev. Lett. 1982. Vol. 48,'№20, pp. 1425-1428.
9. H. Cheng H., Greebgard L., Rokhlin V. A Fast Adaptive Multipole Algorithm in Three Dimensions //Journal of Computational Physics. 1999. Vol.155, №2, pp. 468-498.
10. Mortensen J.J, Hansen L.B., Jacobsen K.W. Real-space grid implementation of the projector augmented wave method //Phvs. Rev. В Condensed Matter. 2005. Vol.71, №3, pp. 035109-1035109-11.
11. Chelikowskv J.R., Wu K.,Troullier N., Saad Y. Higher-order finite-difference pseudopotential method: An application to diatomic molecules //Phvs. Rev. B. 1994. Vol.50, №16, pp. 1135511364.
12. Becke A. D. A multicenter numerical integration scheme for polyatomic molecules //J. Chem. Phvs. 1988. Vol.88, №4,pp. 2547-2553.
13. Kikuji Hirose, Tomova Ono, Yoshitaka Fujimoto, Shigeru Tsukamoto. First-Principles Calculations in Real-Space Formalism. - London: Imperial College Press, 2005. - 2253 c.
14. Atsuvuki Okabe, Barry Boots, Kokichi Sugihara, Sung Nok Chiu. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. - New York: Wiley, 2000. - 696 c.
15. Gonze X., Stumpf R., Scheffler M. Analysis of separable potentials //Phvs. Rev. B. 1991. Vol.44, №16, pp. 8503-8513.
16. Troullier N.,Martins J.L. Efficient pseudopotentials for plane-wave calculations //Phvs. Rev. 1991. Vol.43, №3, pp. 1993-2006.
17. Brandt A. Multi-level adaptive solutions to boundary-value problems //Mathematics of Computation. 1977. Vol. 31. №138. pp. 333-390.
18. Харрисон У. Электронная структура и свойства твердых тел.: Пер. с англ. - М.: Мир, 1983.
- Т.1. - 381 с.
19. Заводинский В. Г. Квантовое моделирование многоатомных систем без волновых функций.
- LAP LAMBERT Academic Publishing RU, 2017,- 56 с.
REFERENCES
1. Landau L.D., Lifshits E.M. (1963), Quantum mechanics. Nonrelativistic theory. [Kvantovava mekhanika. Nerelvativistskava teoriva.], Fizmatgiz, Moscow, 768 p.
2. Kohn W., Sham J.L. 1965, "Self-consistent equations including exchange and correlation effects", Phys. Rev., vol. 140, no.4A., pp. A1133-A1138.
3. Zavodinskv V. (2013) Computer modeling of nanoparticles and nanosvstems. [Komp'yuternove modelirovanive nanochastits i nanosistem.], Fizmatlit, Moscow, 137 p.
4. Skollermo G. 1975, "A Fourier method for the Numerical Solution of Poisson Equation", Mathematics of Computation, vol. 29, nol31. pp. 697-711.
5. Chun- Min Chang, Yihan Shao, Jing Kong. 2012, "Ewald mesh method for quantum mechanical calculations", J. Chem. Phys, vol. 136, noil, pp. 114112-114112-5.
6. Bobrov V. В., Trigger S.A. 2013, "The problem of the universal density functional and the density matrix functional theory", Journal of Experimental and Theoretical Physics, vol. 116, no4, pp. 635-640.
7. Chelikowskv J.R., Troullier N, Saad Y. 1994. "Finite-difference-pseudopotential method: Electronic structure calculations without a basis", Phys. Rev. Lett. vol. 72, no8, pp. 1240-1243.
8. L. Kleinman L., Bvlander D.M. 1982. "Efficacious Form for Model Pseudopotentials", Phys. Rev. Lett. vol. 48, no20, pp. 1425-1428.
9. H. Cheng H., Greebgard L., Rokhlin V. 1999. "A Fast Adaptive Multipole Algorithm in Three Dimensions", Journal of Computational Physics, vol.155, no2, pp. 468-498.
10. Mortensen J.J, Hansen L.B., Jacobsen K.W. 2005. "Real-space grid implementation of the projector augmented wave method ", Phys. Rev. В Condensed Matter, vol.71, no3, pp. 0351091-035109-11.
11. Chelikowskv J.R., Wu K.,Troullier N., Saad Y. 1994. "Higher-order finite-difference pseudopotential method: An application to diatomic molecules", Phys. Rev. B. vol.50, nol6, pp.11355-11364.
12. Becke A.D. 1988. "A multicenter numerical integration scheme for polyatomic molecules", J. Chem. Phys. vol.88, no4,pp. 2547-2553.
13. Kikuji Hirose, Tomova Ono, Yoshitaka Fujimoto, Shigeru Tsukamoto. (2005), First-Principles Calculations in Real-Space Formalism, Imperial College Press, London, 2253 p.
14. Atsuvuki Okabe, Barry Boots, Kokichi Sugihara, Sung Nok Chiu. (2000), Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, Wiley,New York, 696 p.
15. Gonze X., Stumpf R., Scheffler M. 1991. "Analysis of separable potentials", Phys. Rev. B. vol.44, nol6, pp. 8503-8513.
16. Troullier N.,Martins J.L. 1991. "Efficient pseudopotentials for plane-wave calculations", Phys. Rev. vol.43, №3, pp. 1993-2006.
17. Brandt A. 1977. "Multi-level adaptive solutions to boundary-value problems", Mathematics of Computation, vol. 31. nol38. pp. 333-390.
18. Walter A. Harrison. (1980), Electronic Structure and the Properties of Solids, W. H. Freeman and Company, San Francisco, 680 p.
19. Zavodinskv V. (2017) Quantum modeling of polyatomic systems without wave functions. [Kvantovoe modelirovanive mnogoatomnikh sistem bez volnovikh funktsii], LAP LAMBERT Academic Publishing RU, 56 p.
Получено 06.05.2018
Принято в печать 17.08.2018
