1. На главную
  2. Электронные версии
  3. Доклады Российской академии наук. Математика, информатика, процессы управления
  4. T. 509, № 1, 2023

Доклады Российской академии наук. Математика, информатика, процессы управления, 2023, T. 509, № 1, стр. 23-27

Способ преобразования спектра оператора в уравнениях Хартри–Фока и Кона–Шэма

А. А. Даньшин 1, *, член-корреспондент РАН А. А. Ковалишин 1, **

1 Национальный исследовательский центр “Курчатовский институт”
Москва, Россия

* E-mail: danshin_aa@nrcki.ru
** E-mail: kovalishin_aa@nrcki.ru

Поступила в редакцию 23.09.2022
После доработки 20.10.2022
Принята к публикации 20.12.2022

Полный текст (PDF)

Аннотация

В работе предложен способ предварительного преобразования спектра оператора уравнений метода Хартри–Фока и теории функционала плотности, позволяющий перейти от решения полной проблемы собственных значений к частичной, причем собственные функции оказываются упорядоченными удобным для расчета способом. Эти преобразования позволяют сделать старую идею сеточной аппроксимации решения конкурентоспособной с точки зрения скорости вычислений по сравнению с широко используемыми подходами на основе базисных наборов.

Ключевые слова: метод Хартри–Фока, теория функционала плотности, проблема собственных значений

1. ВВЕДЕНИЕ

Теория функционала плотности [1–3] и метод Хартри–Фока [4, 5] являются основными инструментами в теории электронной структуры в химии, материаловедении и в других смежных областях. Традиционные подходы к решению уравнений Кона–Шэма и Хартри–Фока включают методы с использованием явных базисных наборов (например, базисы плоских волн [6, 7], локализованные базисные наборы [8, 9]) и методы без их использования, т.е. методы реального пространства [10, 11]. Достоинства методов реального пространства понятны: такие подходы не приводят к возникновению дополнительных ошибок, связанных с центрированием базисных наборов; концептуально просты для реализации в высокопараллельных средах и позволяют избежать глобальных взаимодействий; не навязывают искусственной периодичности и др.

Идея проводить расчеты без базисных наборов, а опираясь только на конечные разности в реальном пространстве, исследуется довольно давно. Гамильтоновы матрицы, получаемые в результате дискретизации уравнений в реальном пространстве, есть матрицы гораздо большей размерности, чем матрицы, полученные в результате разложения по базисным наборам. Наиболее затратной частью вычислений является решение задачи на собственные значения для каждого электрона на каждом шаге цикла по самосогласованию (SCF). При этом те методы диагонализации, которые применяются для решения задачи на собственные значения, требуют значительных вычислительных ресурсов для систем с большим числом электронов. Применение различных приемов для ускорения, таких как введение в рассмотрение только валентных электронов [12, 13], или использование чебышевских методов для ускорения поиска собственных значений [14, 15] не привело к тому, что методы реального пространства могли бы конкурировать с методами базисных наборов в плане вычислительных затрат.

Целью настоящей работы является создание асимптотически точного в смысле шага расчетной сетки алгоритма решения уравнений метода Хартри–Фока и теории функционала плотности методами реального пространства, вычислительная сложность которого будет сопоставима со сложностью алгоритмов, реализующих метод базисных наборов. Здесь мы так же, как и в большинстве работ, посвященных этому направлению, опираемся на конечно-разностные уравнения, но предварительно тождественно преобразовываем спектр конечно-разностного оператора, и, благодаря этому, поиск собственных функций осуществляется последовательно, начиная с основного состояния. Такой подход позволяет значительно ускорить вычисления за счет отказа от методов диагонализации. Используемые при этом математические приемы хорошо известны из задач переноса нейтронов в физике ядерных реакторов. Представленный алгоритм верифицирован на полном наборе элементов таблицы Менделеева.

2. ПРОБЛЕМА СОБСТВЕННЫХ ЗНАЧЕНИЙ

Уравнения Кона–Шэма и Хартри–Фока есть уравнения вида:

(1)
$\begin{gathered} \Delta {{\psi }_{k}}({\mathbf{r}}) - J({\mathbf{r}}){{\psi }_{k}}({\mathbf{r}}) + \\ + \;U({\mathbf{r}}){{\psi }_{k}}({\mathbf{r}}) + K({\mathbf{r}},{{\psi }_{k}}) = {{E}_{k}}{{\psi }_{k}}({\mathbf{r}}), \\ k = \overline {1,N} , \\ \end{gathered} $
где:

${{\psi }_{k}}({\mathbf{r}})$ – одноэлектронная орбитальная функция координат в трехмерном пространстве с номером $k$, собственная функция задачи (1), $N$ – число электронов в системе;

${{E}_{k}}$ – энергия $k$-го электрона в Ридбергах, взятая с обратным знаком, собственное значение задачи (1);

$J({\mathbf{r}}) = 2\int \frac{{n({\mathbf{r}}{\kern 1pt} ')}}{{{\text{|}}{\mathbf{r}} - {\mathbf{r}}{\kern 1pt} '{\text{|}}}}d{\mathbf{r}}{\kern 1pt} '$ – кулоновский потенциал системы электронов, распределенных с плотностью $n$;

$U({\mathbf{r}}) = \sum\nolimits_{j = 1}^M {\frac{{2{{Z}_{j}}}}{{{\text{|}}{\mathbf{r}} - {{{\mathbf{a}}}_{j}}{\text{|}}}}} $ – потенциал ядер, $M$ – число ядер в системе, ${{{\mathbf{a}}}_{j}}$ – координаты ядер, ${{Z}_{j}}$ – заряд ядер;

$K({\mathbf{r}},{{\psi }_{k}})$ – обменно-корреляционный член, зависящий от используемого приближения.

Наиболее затратной частью вычислений является повторное решение задачи на собственные значения (1) для каждого $k$-го электрона на каждом $m$-м шаге итераций по самосогласованию:

(2)
$H\psi _{k}^{{(m)}} = {{E}_{k}}\psi _{k}^{{(m)}},$
где $H \in {{\mathbb{R}}^{{L \times L}}}$, ${{\psi }_{k}} \in {{\mathbb{R}}^{L}}$, $L$ – размерность задачи (число расчетных узлов после дискретизации уравнения (1)). Вычисление собственных пар напрямую из $H$ стоит $O\left( {{{L}^{3}}} \right)$, что очень дорого с вычислительной точки зрения для больших $L$. Собственные решатели для разреженных матриц имеют сложность $O\left( {{{L}^{2}}N} \right)$, что по-прежнему дорого, а задача (2) решается много раз на одной итерации SCF. Также для данной задачи развит подход спектральных фильтров на основе полиномов Чебышева, однако эффективность построенного фильтра зависит от выбора границ фильтрации [15]. Использование простого сдвига спектра не дает существенного выигрыша из-за плохой обусловленности задачи.

Качественно структура собственных значений уравнений (1) известна – это знакопеременный спектр, ограниченный снизу минимальным значением (отрицательная энергия основного состояния) и неограниченный сверху (положительная энергия свободных состояний). При дискретизации уравнений неограниченный спектр заменяется ограниченным: за счет конечного числа узлов сетки качество аппроксимации собственных функций верхней части спектра падает. При этом качество аппроксимации собственных функций, соответствующих связанным состояниям, остается достаточно высоким благодаря их относительной гладкости.

Теперь необходимо сделать несколько пояснений относительно дальнейшего изложения, поскольку это место является ключевым в данной работе. Уравнение Шредингера и соответствующие ему системы уравнений Хартри–Фока и теории функционала плотности с точки зрения структуры аналогичны уравнению диффузии нейтральных частиц в среде с поглощением и размножением. В случае нейтронной задачи с целью расчета ядерных реакторов были разработаны алгоритмы и приемы преобразования спектра для обеспечения устойчивости и лучшей сходимости итерационных методов. Мы используем эти подходы, хорошо изложенные, например, в [16, 17], для решения уравнений вида (1). Проводя аналогию с уравнением диффузии нейтронов, члены $J$ и ${{E}_{k}}$ в (1) отвечают за поглощение (они неотрицательны и “вычитаются” из лапласиана), а члены $K$ и $U$ отвечают за размножение (они неотрицательны и “складываются” с лапласианом). Далее, следуя аналогии, введем вспомогательное собственное значение ${{k}_{{{\text{eff}}}}}$ следующим образом:

(3)
$\begin{gathered} - \Delta {{\psi }_{k}}({\mathbf{r}}) + J({\mathbf{r}}){{\psi }_{k}}({\mathbf{r}}) + {{E}_{k}}{{\psi }_{k}}({\mathbf{r}}) = \\ = \frac{1}{{{{k}_{{{\text{eff}}}}}}}\left( {U({\mathbf{r}}){{\psi }_{k}}({\mathbf{r}}) + K({\mathbf{r}},{{\psi }_{k}})} \right). \\ \end{gathered} $

При фиксированном ${{E}_{k}}$ собственные значения ${{k}_{{{\text{eff}}}}}$ в данном случае будут упорядочены удобным для нас способом с сохранением свойств (1): теперь спектр устроен так, что поиск собственных функций осуществляется последовательно, начиная с основного состояния. Это ключевое положение асимптотической теории диффузии нейтронов, а ${{k}_{{{\text{eff}}}}}$ носит название эффективного коэффициента размножения. Очевидно, что уравнение (3) будет совпадать с уравнением (1) при ${{k}_{{\operatorname{eff} }}} = 1$. Построение эффективных устойчивых методов решения уравнений вида (3) хорошо известно и описано, например, в [18]. В основе алгоритма лежит итерационный процесс, когда после решения уравнения (3) мы “поправляем” значение ${{E}_{k}}$ так, чтобы получить ${{k}_{{\operatorname{eff} }}} = 1$.

Запишем (3) в виде:

(4)
$\widehat T{{\psi }_{k}} = \frac{1}{{{{k}_{{\operatorname{eff} }}}}}\widehat Q{{\psi }_{k}}.$

Пусть задача (4) решена для заданного ${{\tilde {E}}_{k}} = {{E}_{k}} + \delta {{E}_{k}}$, т.е. найдено решение $\frac{1}{{{{{\widetilde k}}_{{\operatorname{eff} }}}}}$ = $\frac{1}{{{{k}_{{\operatorname{eff} }}}}}$ + + $\delta \frac{1}{{{{k}_{{\operatorname{eff} }}}}}$ и ${{\tilde {\psi }}_{k}} = {{\psi }_{k}} + \delta {{\psi }_{k}}$:

(5)
$\left( {\widehat T + \delta \widehat T} \right){{\tilde {\psi }}_{k}} = \frac{1}{{{{{\widetilde k}}_{{eff}}}}}\widehat Q{{\tilde {\psi }}_{k}}.$

Умножим скалярно (4) на ${{\tilde {\psi }}_{k}}$, (5) – на ${{\psi }_{k}}$, и вычтем одно из другого. Тогда с учетом $\delta \widehat T \equiv \delta {{E}_{k}}$ в первом порядке теории возмущений:

(6)
$\left\langle {{{\psi }_{k}}} \right|\delta {{E}_{k}}\left| {{{\psi }_{k}}} \right\rangle = \delta \frac{1}{{{{k}_{{\operatorname{eff} }}}}}\left\langle {{{\psi }_{k}}} \right|\widehat Q\left| {{{\psi }_{k}}} \right\rangle .$

Окончательно (с учетом ${{k}_{{\operatorname{eff} }}} \equiv 1$ для точного решения (1)), в соответствии с теорией малых возмущений вычисляется новое значение:

(7)
${{E}_{k}} = {{\tilde {E}}_{k}} + \left( {1 - \frac{1}{{{{{\widetilde k}}_{{\operatorname{eff} }}}}}} \right)\frac{{\left\langle {{{\psi }_{k}}} \right|\widehat Q\left| {{{\psi }_{k}}} \right\rangle }}{{\left\langle {{{\psi }_{k}}|{{\psi }_{k}}} \right\rangle }}.$

Затем процедура решения (3) повторяется. Организованный таким образом итерационный процесс приведет к тому, что ${{k}_{{\operatorname{eff} }}} \to 1$, а собственные значения и собственные функции сойдутся к искомым для уравнений (1). В этом случае аппроксимацию операторов можно выполнить с помощью конечных разностей. Тогда мы получаем разреженную матрицу, которую удобно использовать в итерационных процедурах благодаря экономичности операций произведения разреженной матрицы на вектор. Дополнительное преимущество в скорости можно получить, если использовать конечно-разностные схемы высокого порядка аппроксимации.

3. АЛГОРИТМ

Алгоритм может быть представлен в следующем виде:

$\left[ \begin{gathered} {\text{Цикл по самосогласованию}} \hfill \\ \left[ \begin{gathered} {\text{Цикл по орбитальным функциям}} \hfill \\ \left[ \begin{gathered} {\text{Расчет кулоновского потенциала}}\;J \hfill \\ \left[ \begin{gathered} {\text{Цикл итераций по теории малых возмущений}} \hfill \\ \left[ {{\text{Решение задачи (3) методом итераций по источнику}}} \right. \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered} \right.$

Выход из цикла по самосогласованию осуществляется, когда достигнута сходимость орбитальных функций в соседних итерациях по отношению к заданному значению параметра сходимости. В цикле по орбитальным функциям решается задача (3) для заданного числа $N$ электронов. Цикл итераций по теории малых возмущений организован в соответствии с (7).

Теперь рассмотрим решение задачи (3) методом итераций по источнику. Введем вектор $\Psi $ размера $L$, каждый элемент которого есть значение искомой функции в соответствующем узле сетки, и запишем уравнения (3) в матричной форме:

(8)
${\mathbf{T}}\Psi = \frac{1}{{{{k}_{{\operatorname{eff} }}}}}{\mathbf{Q}}\Psi .$

Пусть первые $i$ орбитальных функций ${{\Psi }_{k}}$, $k = \overline {1,i} $ найдены (при этом для остальных орбитальных функций берутся значения, вычисленные на предыдущей итерации цикла по самосогласованию), а $n - 1$ и $n$ – номера итераций по источнику. Тогда следующий алгоритм используется для поиска cобственного значения и собственной функции $(i + 1)$-го электрона:

1. Нормировка: $\left\langle {{{\Psi }^{{(n - 1)}}}\,{\text{|}}\,{{\Psi }^{{(n - 1)}}}} \right\rangle $ = 1.

2. Формирование источника: ${{{\mathbf{S}}}^{{(n - 1)}}} = {\mathbf{Q}}{{\Psi }^{{(n - 1)}}}$.

3. Решение неоднородной задачи: ${\mathbf{T}}{{\widetilde \Psi }^{{(n)}}} = {{{\mathbf{S}}}^{{(n - 1)}}}$, где ${{\widetilde \Psi }^{{(n)}}}$ – искомая функция.

4. Ортогонализация к первым $\widetilde i$ орбитальным функциям электронов c тем же значением проекции спина: ${{\Psi }^{{(n)}}} = {{\widetilde \Psi }^{{(n)}}}$ – $\sum\limits_k^{\tilde {i}} \left\langle {{{{\widetilde \Psi }}^{{(n)}}}\,{\text{|}}\,{{\Psi }_{k}}} \right\rangle {{\Psi }_{k}}$.

5. Оценка ${{k}_{{\operatorname{eff} }}}$: ${{k}_{{\operatorname{eff} }}} = \left\langle {{{\Psi }^{{(n)}}}\,{\text{|}}\,{{\Psi }^{{(n - 1)}}}} \right\rangle $.

6. Оценка сходимости по ${{k}_{{\operatorname{eff} }}}$ и собственной функции. Если заданная точность по разнице значений ${{k}_{{\operatorname{eff} }}}$ и собственной функции в соседних итерациях не достигнута, повтор, начиная с пункта 1.

При выполнении пункта 3 необходимо решать неоднородную систему алгебраических уравнений с разреженной матрицей. Здесь можно использовать стандартные итерационные процедуры для решения таких систем, например, метод верхней релаксации [19] или чебышевские итерационные методы [20]. Значения $J$ и $K$ вычисляются стандартным способом через решения уравнений Пуассона с соответствующей правой частью [10–13].

4. ВЕРИФИКАЦИЯ

Алгоритм реализован в виде компьютерной программы. Мы использовали простейшую схему второго порядка для аппроксимации оператора Лапласа. Для верификации использовались доступные значения полных энергий основного состояния атомов из работ [21–23], в которых фактически был достигнут хартри-фоковский предел в результате вычислений на основе метода базисных наборов. Получено, что с уменьшением шага сетки решение сходится к значениям работ [21–23]. Так, например, для шага $h = {{10}^{{ - 3}}}{{a}_{0}}$ (${{a}_{0}}$ – боровский радиус) разница в решении для элементов от He до U лежит в диапазоне от ${{10}^{{ - 5}}}$ до 210 эВ, а для шага $h = {{10}^{{ - 4}}}{{a}_{0}}$ – от ${{10}^{{ - 7}}}$ до 2.1 эВ. Результаты отклонений рассчитанных значений полных энергий основного состояния от значений работ [21–23] представлены на рис. 1. При дальнейшем уменьшении шага сетки решение асимптотически сходится к точному для рассматриваемых уравнений.

Рис. 1.

Разница в энергии основного состояния между рассчитанным значением и значением, полученным в работах [21–23], в зависимости от номера элемента N.

5. ВЫВОДЫ

Авторами создан и верифицирован асимптотически точный в смысле шага расчетной сетки алгоритм решения уравнений метода Хартри–Фока и теории функционала плотности методами реального пространства. Предварительное тождественное преобразование спектра конечно-разностного оператора позволило сразу находить решение для основного состояния системы, что привело к значительному ускорению за счет отказа от методов диагонализации. Данный подход делает старую идею сеточной аппроксимации решения конкурентоспособной с точки зрения скорости вычислений по сравнению с широко используемыми подходами на основе метода базисных наборов. Отметим, что реализация схем более высокого порядка позволит дополнительно ускорить вычисления.

Список литературы

  1. Hohenberg P., Kohn W. Inhomogeneous Electron Gas // Phys. Rev. 1973. V. 136. P. B864–B871.

  2. Kohn W., Sham L.J. Self-Consistent Equations Including Exchange and Correlation Effects // Phys. Rev. 1965. V. 140. P. A1133–A1138.

  3. Kohn W. Nobel Lecture: Electronic structure of matter–wave functions and density functionals // Rev. Mod. Phys. 1999. V. 71. P. 1253–1266.

  4. Фок В.А. Приближенный способ решения квантовой задачи многих тел // Успехи физических наук. 1967. Т. 93. С. 342–363.

  5. Pople J.A. Nobel Lecture: Quantum chemical models // Rev. Mod. Phys. 1999. V. 71. P. 1267–1274.

  6. Payne M.C., Teter M.P., Allan D.C., Arias T.A., Joannopoulos J.D. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients // Rev. Mod. Phys. 1992. V. 64. P. 1045–1097.

  7. Kresse G., Furthmuller J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set // Phys. Rev. B. 1996. V. 54. P. 11169–11186.

  8. Koch W., Holthausen M.C. A Chemist’s Guide to Density Functional Theory, 2nd Edition. Wiley-VCH, Weinheim, 2001. 313 p.

  9. Martin R.M. Electronic Structure: Basic Theory and Practical Methods. Cambridge University Press, 2004. 624 p.

  10. Beck T.L. Real-space mesh techniques in density-functional theory // Rev. Mod. Phys. 2000. V. 72. P. 1041–1080.

  11. Torsti T. et al. Three real-space discretization techniques in electronic structure calculations // physica status solidi (b). 2006. V. 243. P. 1016–1053.

  12. Marques M., Castro A., Bertsch G., Rubio A. octopus: a first-principles tool for excited electron-ion dynamics // Computer Physics Communications. 2003. V. 151. P. 60–78.

  13. Kronik L. et al. PARSEC – the pseudopotential algorithm for real-space electronic structure calculations: recent advances and novel applications to nano-structures // physica status solidi (b). 2006. V. 243. P. 1063–1079.

  14. Zhou Y., Saad Y., Tiago M., Chelikowsky J. Parallel self-consistent-field calculations via Chebyshev-filtered subspace acceleration // Phys. Rev. E. 2006. V. 74. P. 066704.

  15. Zhou Y., Chelikowsky J., Saad Y. Chebyshev-filtered subspace iteration method free of sparse diagonalization for solving the Kohn-Sham equation // Journal of Computational Physics. 2014. V. 274. P. 770–782.

  16. Шихов С.Б. Вопросы математической теории реакторов. Линейный анализ. М.: Атомиздат, 1973. 376 с.

  17. Крянев А.В., Шихов С.Б. Вопросы математической теории реакторов: Нелинейный анализ. М.: Энергоатомиздат, 1983. 280 с.

  18. Марчук Г.И., Лебедев В.И. Численные методы в теории переноса нейтронов. М.: Атомиздат, 1981. 456 с.

  19. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 432 с.

  20. Лебедев В.И. Функциональный анализ и вычислительная математика. М.: Физматлит, 2005. 296 с.

  21. Bunge C.F., Barrientos J.A., Bunge A.V., Cogordan J.A. Hartree-Fock and Roothaan-Hartree-Fock energies for the ground states of He through Xe // Phys. Rev. A. 1992. V. 46. P. 3691–3696.

  22. Koga T., Tatewaki H., Thakkar A. Roothaan-Hartree-Fock wave functions for atoms with Z$ \leqslant $54 // Phys. Rev. A. 1993. V. 47. P. 4510–4512.

  23. Koga T., Thakkar A. Roothaan-Hartree-Fock wave functions for atoms from Cs through U // Phys. Rev. A. 1993. V. 48. P. 4775–4777.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Доклады Российской академии наук. Математика, информатика, процессы управления

Архивы выпусков Информация о журнале Отправить рукопись в журнал