Теория функционала плотности DFT: суть и приближение LDA

Теория функционала плотности (DFT, density functional theory) - метод квантовой химии, который заменяет волновую функцию системы из $N$ электронов гораздо более простым объектом: трёхмерной электронной плотностью $\rho(\mathbf{r})$. Вместо функции $3N$ переменных нужно знать всего одну функцию трёх координат, и по теореме Хоэнберга-Кона этого достаточно, чтобы в принципе найти любое свойство основного состояния - энергию, геометрию, спектры. Именно поэтому DFT стала рабочей лошадкой расчётной химии: она даёт приемлемую точность при вычислительных затратах на порядки ниже, чем у методов, явно строящих волновую функцию. Ниже разберём, откуда берётся такая экономия, что такое уравнения Кона-Шэма и обменно-корреляционный функционал, и как приближение локальной плотности (LDA) на модели однородного электронного газа даёт первые рабочие числа. Чтобы сразу увидеть, как энергия зависит от плотности, покрути калькулятор ниже - дальше разберём каждую формулу по отдельности.

Теоремы Хоэнберга-Кона: почему хватает электронной плотности

Первая теорема Хоэнберга-Кона утверждает: внешний потенциал $v(\mathbf{r})$, действующий на электроны (притяжение к ядрам), однозначно определяется плотностью основного состояния $\rho(\mathbf{r})$ - с точностью до аддитивной константы. Раз потенциал задаёт весь гамильтониан, а значит и волновую функцию, и все наблюдаемые, то плотность $\rho(\mathbf{r})$ фактически несёт в себе всю информацию о системе. Вторая теорема добавляет вариационный принцип: существует универсальный функционал энергии $E[\rho]$, который минимален именно на истинной плотности основного состояния. Формально энергию расписывают как

$E[\rho] = T_s[\rho] + E_{ne}[\rho] + J[\rho] + E_{xc}[\rho],$

где $T_s[\rho]$ - кинетическая энергия эталонной системы невзаимодействующих электронов с той же плотностью, $E_{ne}[\rho]$ - энергия притяжения к ядрам, $J[\rho]$ - классическое кулоновское отталкивание электронов, а $E_{xc}[\rho]$ - обменно-корреляционная энергия, куда спрятано всё остальное: квантовый обмен, корреляция движений электронов и поправка к кинетической энергии за счёт взаимодействия. Теоремы доказывают, что такой функционал существует, но не дают его явного вида - в этом и состоит главная трудность DFT.

Уравнения Кона-Шэма: от электронов к эффективным орбиталям

Прямая минимизация $E[\rho]$ по плотности неудобна, потому что кинетическую энергию взаимодействующей системы трудно выразить через $\rho(\mathbf{r})$ напрямую. Кон и Шэм предложили обходной манёвр: ввести вспомогательную систему невзаимодействующих электронов, которая создаёт ту же самую плотность $\rho(\mathbf{r})$, что и настоящая взаимодействующая система. Тогда $T_s[\rho]$ считается точно через одноэлектронные орбитали $\varphi_i$, а вся сложность взаимодействия остаётся в $E_{xc}[\rho]$. Минимизация энергии приводит к системе одноэлектронных уравнений Кона-Шэма:

$\left[-\frac{1}{2}\nabla^2 + v_{eff}(\mathbf{r})\right]\varphi_i(\mathbf{r}) = \varepsilon_i \varphi_i(\mathbf{r}),$

где эффективный потенциал $v_{eff}$ включает притяжение к ядрам, кулоновский потенциал электронной плотности и обменно-корреляционный потенциал $v_{xc} = \delta E_{xc}/\delta\rho$. Плотность восстанавливается суммой по занятым орбиталям, $\rho(\mathbf{r}) = \sum_i |\varphi_i(\mathbf{r})|^2$, а сами орбитали зависят от $\rho$ через $v_{eff}$ - отсюда и требование самосогласованности: уравнения решают итеративно, пока входная и выходная плотность не совпадут.

Обменно-корреляционный функционал и приближение LDA

Всё практическое различие между методами DFT сводится к тому, как именно приближают неизвестный функционал $E_{xc}[\rho]$. Простейшее и исторически первое приближение - LDA (local density approximation, локальная плотность): энергия обмена и корреляции в точке $\mathbf{r}$ берётся такой же, как в однородном электронном газе той же плотности $\rho(\mathbf{r})$ в этой точке, без учёта того, как плотность меняется в соседних точках.

Рассчитать для своей задачи

Для однородного газа обменная энергия на один электрон известна точно (обмен Дирака):

$\varepsilon_x(\rho) = -\frac{3}{4}\left(\frac{3}{\pi}\right)^{1/3}\rho^{1/3} = -C_x\,\rho^{1/3}, \qquad C_x \approx 0{,}7386.$

Это ключевая идея LDA: обменно-корреляционная энергия всей системы получается интегрированием локального выражения по всему пространству, $E_{xc}[\rho] \approx \int \varepsilon_{xc}(\rho(\mathbf{r}))\,\rho(\mathbf{r})\,d\mathbf{r}$ - функционал зависит только от значения плотности в каждой точке, а не от её градиентов (это уже уровень GGA, более точного, но более тяжёлого приближения).

Однородный электронный газ и формула Вигнера

Точного выражения для корреляционной энергии однородного газа нет - но есть надёжные интерполяции, подогнанные под предельные случаи и численные расчёты методом Монте-Карло. Одна из первых и самых простых - формула Вигнера:

$\varepsilon_c(r_s) = -\frac{0{,}44}{r_s + 7{,}8}\ \text{хартри},$

где $r_s$ - радиус Вигнера-Зейтца, средний радиус сферы, приходящейся на один электрон: $\rho = 3/(4\pi r_s^3)$. Малые $r_s$ отвечают плотному электронному газу (металлы с высокой плотностью проводимости), большие $r_s$ - разреженному газу.

Рассчитать для своей задачи

Полезно сразу дать кинетическую энергию невзаимодействующей эталонной системы в том же локальном духе - приближение Томаса-Ферми:

$t(\rho) = \frac{3}{10}(3\pi^2)^{2/3}\rho^{2/3} = C_F\,\rho^{2/3}, \qquad C_F \approx 2{,}871.$

Разберём числа: пусть $r_s = 4$ бор - типичный порядок для валентных электронов металла. Тогда плотность $\rho = 3/(4\pi\cdot 4^3) \approx 0{,}00373$ электрон/бор³, обменная энергия на электрон $\varepsilon_x \approx -3{,}12$ эВ, корреляционная $\varepsilon_c \approx -1{,}01$ эВ, а сумма $\varepsilon_{xc} \approx -4{,}13$ эВ. Кинетическая энергия по Томасу-Ферми на электрон здесь $t \approx 1{,}88$ эВ. Для модельной системы из $N = 8$ таких электронов полные вклады умножаются на $N$: кинетическая $\approx 15{,}0$ эВ, обменная $\approx -24{,}9$ эВ, корреляционная $\approx -8{,}1$ эВ, суммарная обменно-корреляционная $\approx -33{,}1$ эВ. Именно эти четыре числа при выбранных $r_s$ и $N$ показывает калькулятор выше.

Самосогласованный расчёт DFT: как ищут плотность на практике

На практике уравнения Кона-Шэма решают итерационно, потому что эффективный потенциал сам зависит от искомой плотности. Цикл выглядит так: берут стартовое приближение $\rho^{(0)}$ (например, суперпозицию атомных плотностей), строят по нему $v_{eff}$, решают одноэлектронные уравнения и получают новые орбитали и новую плотность $\rho^{(1)}$. Если $\rho^{(1)}$ заметно отличается от $\rho^{(0)}$, смешивают их (чтобы не расходиться) и повторяют цикл - до тех пор, пока плотность на входе и выходе не совпадут в пределах заданного порога. Такой самосогласованный расчёт (SCF, self-consistent field) - стандартная процедура почти всех методов квантовой химии, но в DFT она особенно дешёва: одноэлектронные уравнения Кона-Шэма решаются намного быстрее, чем уравнение для полной многоэлектронной волновой функции.

Где DFT ошибается и как это учитывают

LDA, при всей своей простоте, систематически завышает энергию связи молекул и занижает длины некоторых связей, потому что предполагает, что электронная плотность в молекуле локально ведёт себя как в однородном газе - а на самом деле она резко меняется вблизи ядер и на периферии электронного облака. Следующий уровень точности, GGA (обобщённое градиентное приближение), добавляет в функционал зависимость от градиента $\nabla\rho$, а гибридные функционалы (вроде B3LYP) подмешивают долю точного обмена Хартри-Фока. Все они наследуют одну и ту же архитектуру Кона-Шэма - меняется только конкретный вид $E_{xc}[\rho]$, и именно поэтому сравнение LDA-чисел с более точными функционалами - стандартное упражнение при знакомстве с DFT.

Частые ошибки

• Путать плотность $\rho(\mathbf{r})$ с волновой функцией. DFT работает с распределением заряда в пространстве, а не с многоэлектронной волновой функцией - теоремы Хоэнберга-Кона именно это и разрешают.
• Забывать про самосогласованность. Уравнения Кона-Шэма нельзя решить за один проход: эффективный потенциал зависит от плотности, которую сам же и определяет, поэтому нужен итеративный цикл.
• Путать LDA и точный обмен Хартри-Фока. Обмен Дирака в LDA - приближённая локальная формула для однородного газа, а не точный обменный интеграл Хартри-Фока; совпадают они только в пределе однородной плотности.
• Подставлять $r_s$ в ангстремах вместо боров. Радиус Вигнера-Зейтка традиционно задают в атомных единицах длины (борах), иначе константы 0,44 и 7,8 в формуле Вигнера не работают.
• Считать, что $E_{xc}$ - это только обмен. Обменно-корреляционный член объединяет квантовый обмен (антисимметрия) и корреляцию (совместное движение электронов) - оба слагаемых нужны для правильного порядка величины энергии.

FAQ

Чем плотность $\rho(\mathbf{r})$ лучше волновой функции для расчёта молекул? Волновая функция системы из $N$ электронов зависит от $3N$ координат (плюс спин), а плотность - всего от трёх. По теоремам Хоэнберга-Кона плотность основного состояния однозначно задаёт внешний потенциал и, значит, всю систему, поэтому расчёт можно вести через гораздо более компактный объект.

Что физически означает обменно-корреляционный функционал $E_{xc}[\rho]$? Это поправка к классической (кулоновской) энергии электрон-электронного отталкивания, которая учитывает квантовые эффекты: принцип Паули не даёт электронам с одинаковым спином находиться слишком близко (обмен), а кулоновское отталкивание заставляет электроны в среднем избегать друг друга сильнее, чем предсказывает классика (корреляция).

Почему LDA называют приближением однородного электронного газа? Потому что единственная система, для которой обменная энергия известна точно, а корреляционная - надёжно посчитана численно, это бесконечный газ электронов с постоянной плотностью (модель желе). LDA переносит эти локальные значения на реальную, неоднородную плотность молекулы точка за точкой.

Коротко

Теория функционала плотности заменяет волновую функцию $N$ электронов трёхмерной плотностью $\rho(\mathbf{r})$, которая по теоремам Хоэнберга-Кона однозначно определяет всю систему. Уравнения Кона-Шэма сводят задачу к самосогласованному набору одноэлектронных уравнений, а вся неизвестная физика взаимодействия спрятана в обменно-корреляционном функционале $E_{xc}[\rho]$. Простейшее приближение LDA берёт локальное значение энергии однородного электронного газа: обмен Дирака $\varepsilon_x = -C_x\rho^{1/3}$ и корреляцию по формуле Вигнера $\varepsilon_c = -0{,}44/(r_s+7{,}8)$ - с этих формул начинается количественная DFT.