Квантовая химия. Теория самосогласованного поля.#

Автор(ы):

Описание лекции#

Из этой лекции узнаем:

  • как с помощью квантовой химии предсказать спектр атома водорода “из первых принципов”;

  • как посчитать энергию атома с помощью Python;

  • какие бывают волновые функции электронов и как их вычислять методом Self-Consistent Field;

  • как посчитать энергию спирта.

Введение#

Боровская теория оказалось недостаточной для описания таких сложных систем, как молекулы или не-водородоподобные атомы. В этой лекции разберемся с теорией самосогласованного поля – это один из первых вычислительных методов квантовой химии, пригодный на практике для моделирования разных систем (атомов, молекул, кристаллов).

Несмотря на то, что основа метода была разработана чуть ли не сто лет назад, он активно используется и сегодня, особенно с различными модификациями и дополнениями.

Что ищем?#

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

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

Долгое и правильное решение, часть 1#

Начнем с атома водорода и определения его уровней энергии – их можно сравнить с экспериментальными данными, а также результатами Боровской теории. Квантовая механика позволяет решить эту задачу, исходя только из “первых принципов”, то есть аксиоматики квантовой механики. При этом после некоторых усилий решения успешно обобщаются с атома водорода на любой другой.

Давайте опишем атом водорода – систему из электрона и протона – на языке квантовой механики, то есть в виде уравнения Шредингера, и посмотрим, что получится.

Электрон находится в потенциале протона, и его волновая функция должна зависеть от расстояния до протона.

Уравнение Шредингера для электрона будет следующим:

\[ \dot{\imath} \hbar \frac{d}{dt} {\ket{\Psi(t, \vec{r})}} = (\frac{ \hat{p}^{2} }{2 m} + \hat{V}(\vec{r})) {\ket{\Psi (t, \vec{r})}} \]

Мы ищем решение, в котором электрон остается в атоме, то есть решаем стационарное уравнение, в котором \(\ket{\Psi}\) не зависит от времени. Тогда оператор эволюции (в левой части уравнения Шредингера) при применении к \(\ket{\Psi}\) должен вернуть нам тот же вектор \(\ket{\Psi}\), умноженный на \(E\) (энергия частицы).

С точки зрения математики, искомая волновая функция является собственным вектором оператора эволюции, а энергия – собственным значением (подробнее в разделе Матрицы). Если объяснять “на пальцах”, то оператор эволюции при применении к волновой функции должен вернуть нам новую (эволюционировавшую или изменившуюся во времени) волновую функцию. Если ищем стационарную – не меняющуюся во времени – волновую функцию, то при применении к ней оператора эволюции она не должна изменяться, иначе будет уже не стационарной.

\[ \dot{\imath} \hbar \frac{d}{dt} {\ket{\Psi(t, \vec{r})}} \equiv \dot{\imath} \hbar \frac{d}{dt} {\ket{\Psi(\vec{r})}} \equiv E { \ket{\Psi(\vec{r})} } \]

Оператор импульса раскрывается через градиент так:

\[ \hat{p} = \dot{\imath} \hbar \nabla \]

Потенциал взаимодействия двух частиц с противоположными единичными зарядами раскрывается так:

\[ \hat{V}(\vec{r}) = -\frac{e^2}{r} \]

Итого имеем:

(20)#\[ (-\frac{\hbar^2}{2m}\nabla^2 - \frac{e^2}{r} ){\ket{\Psi(\vec{r})}} = E {\ket{\Psi(\vec{r})}} \]

Прежде чем погружаться глубже в решение этого уравнения, нам надо разобраться с квантовыми числами.

Про квантовые числа#

Теория Бора базируется на идее квантования импульса, и ее следствием является возникшее в формуле энергии (19) число \(n\). Каждое число \(n\) соответствует определенному состоянию, в котором может находиться электрон, и эти состояния отличаются энергией.

Идея, что электрон может находиться в одном состоянии из некоторого конечного набора вариантов, характеризуемого дискретными (квантовыми) числами, оказалась правильной. В дальнейшем она получила развитие в квантовой механике.

Число \(n\) называется главным квантовым числом, оно определяет энергетический уровень электрона.

Из анализа свойств операторов гамильтониана и момента импульса в квантовой механике выводятся еще три квантовых числа:

  • \(\ell\) – орбитальное квантовое число, \(0 \leq \ell \leq n-1\);

  • \(m\) – магнитное квантовое число, \(-\ell \leq m \leq \ell\);

  • \(s\) – спиновое квантовое число, \(s = 0|1\) (для атома водорода оно не играет роли, так как электрон только один).

В химии они используются в другой форме записи для описания орбиталей электронов в атоме.

Эти числа нам скоро понадобятся, чтобы описывать электрон в атоме.

Любая функция \(\ket{\Psi}\), для которой уравнение (20) верно, описывает электрон, который стабильно находится где-то около ядра и не покидает его, то есть образует с ним атом. При этом функций-решений у уравнения на самом деле много, что физически соответствует тому, что электрон может находиться на разных орбиталях. Каждая орбиталь характеризуется квантовыми числами – \(n\), \(\ell\), \(m\) и \(s\), и обозначается как \(\ket{\Psi_{n \ell m}} \). \(\ket{\Psi}\) (орбиталь) с минимальной энергией \(E\) соответствует основному состоянию (ground state) – она описывает невозбужденный электрон. У водорода только один электрон, поэтому единственная \(\ket{\Psi}\) с минимальной энергией соответствует невозбужденному атому водорода.

В целом, основная задача квантовой химии – найти ground state произвольной системы частиц. В реальности материя редко находятся непосредственно в ground state состоянии, т.к. при ненулевой температуре какие-то флуктуации неизбежно возникают. Тем не менее, ground state обычно в наибольшей степени определяет поведение системы, а другие состояния являются “поправками” к нему.

Частицы могут переходить из одного состояния в другое: если электрону в основном состоянии придать энергии, например, попасть в него другой частицей, то он может перейти на другую орбиталь, при этом оставшись частью атома. Со временем возбужденные атомы рано или поздно переходят в основное состояния, выбрасывая лишнюю энергию в виде фотона, что и дает спектр излучения атома.

Итого: \(\ket{\Psi_{n \ell m}(\vec{r})}\) – это функция, соответствующая какой-то стабильной “траектории” (распределению плотности вероятности) электрона вокруг ядра, она является решением уравнения Шредингера, то есть собственной функцией гамильтониана. Соответствующее ей собственное число – энергия электрона на этой орбитали. Для всех волновых функций с одним \(n\) энергия одинакова, но одной энергии может соответствовать несколько волновых функций, каждая - со своим уникальным набором волновых чисел.

Если электрон описывается такой волновой функцией, то он часть атома, и если какой-либо электрон – часть атома, то он описывается такой волновой функцией. Электрон может переходить между этими волновыми функциями, получая и отдавая энергию, оставаясь при этом частью атома.

Волновая функция \(\ket{\Psi_{100}}\) (т.е. с числами \(n=1, \ell=0, m=0\)) соответствует минимальной энергии для единственного электрона в атоме водорода. Т.к. в атоме водорода есть лишь один электрон, и других вкладов в энергию нет – весь атом водорода в такой конфигурации имеет минимальную энергию.

Долгое и правильное решение, часть 2#

А зачем нам вообще сдался спектр и энергии?

Вслед за физиками двадцатого века нам приходится продираться через дебри уравнений, абстракций и формул. Дело это непростое и может возникнуть резонный вопрос – а зачем оно все надо?

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

Если перейти в уравнении (20) в сферические координаты со следующей параметризацией \(\ket{\Psi(\vec{r})} = \frac{\chi(r)}{r} Y(\theta, \phi)\) и воспользоваться несколькими волшебными выводами квантмеха [89], то получится:

\[ -\frac{{\hbar}^2}{2 m} \frac{d^2 \chi}{d r^2} + (-\frac{e^2}{r} + \frac{\hbar^2 \ell (\ell + 1)}{2mr^2} -E) \chi(\vec{r}) = 0 \]

Все константы (массы, заряда, импульса, энергии) можно убрать, если перейти в кулоновские единицы измерений, где они приняты за единицу (то есть 1 единица заряда = заряд электрона, 1 единица массы = масса электрона).

\[ -\frac{\chi''(r)}{2} + (\frac{\ell(\ell+1)}{2r^2} -\frac{1}{r})\chi(r) = 0 \]

Опустим несколько страниц выкладок [15a], учтем граничные условия и получим следующее решение:

\[\begin{split} R_{n\ell}(r) = r^\ell \cdot e^{-r/n} \cdot \sum_{k=0}^{n -\ell -1} \frac{(-2r/n)^k}{ (2\ell + 2 +k)! (n-\ell-k-1)! k!} \cdot C \\ \ket{\Psi(r, \theta, \phi)} = R_{n\ell}(r) Y_{\ell m} (\theta, \phi), \end{split}\]

где \(Y_{\ell m}(\theta, \phi)\) – сферические функции.

Если подставить это решение в уравнение Шредингера и найти энергию, то получим:

\[ E_{n} = -\frac{1}{n^2} \frac{me^4}{2\hbar^2} \]

То есть получим ту же формулу, что и в теории Бора, и тот же численный результат – 13.6 eV.

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

От теории к практике#

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

Сегодня реальные научные задачи решаются численно – и для этого написано множество высокопроизводительных пакетов, позволяющих на основе входных данных и ограничений вычислить определенным методом желаемые характеристики.

Как и в задачах машинного обучения, в первом приближении достаточно подготовить данные, скормить их волшебному программному пакету-вычислителю и забрать ответ, но без понимания происходящего под капотом есть шансы получить что-то неправильное.

В этой лекции воспользуемся Python-пакетом pyscf, который реализует алгоритмы квантовой химии с верхнеуровневым API на Python.

Note

Другим известным Python-пакетом с численными алгоритмами квантовой химии является psi4. Однако он имеет ряд сложностей с установкой, поэтому интерактивные ячейки в этой лекции будут приведены на pyscf, а код на psi4 будет приведен как альтернатива. Библиотека psi4 не устанавливается вместе с данной книгой и, если хотите попробовать ее, то рекомендуется создать новое виртуальное окружение (подробнее о виртуальных окружениях было в вводных главах про Python).

Посчитаем c помощью pyscf энергию основного состояния атома водорода. Некоторые параметры сейчас придется использовать, “поверив на слово”. Их смысл будет объяснен в дальнейшем.

from pyscf import gto, scf

h_atom = gto.Mole()
h_atom.build(
    atom="H 0 0 0",
    basis="STO-3G",
    spin=1,
)

<pyscf.gto.mole.Mole at 0x7f7e54a56580>

Задали атом водорода, по умолчанию атом помещается в начало координат. Можно заметить параметр basis="STO-3G" в задании атома. Тут уже поинтереснее – объяснение этих параметров пока отложим и вернемся после объяснения теории. В целом они определяют, каким именно методом и в каком базисе нужно численно решить уравнение Шредингера.

from scipy.constants import physical_constants

mf = scf.ROHF(h_atom)
e_in_ht = mf.kernel()

h2ev = physical_constants["hartree-electron volt relationship"]

def e_in_ev(energy_in_ht: float) -> float:
    return energy_in_ht * h2ev[0]


print(f"Hydrogen ground state energy: {e_in_ev(e_in_ht)} eV")

converged SCF energy = -0.466581849557275

Hydrogen ground state energy: -12.696338923670483 eV

Здесь посчитали энергию в единицах Хартри – специальной физической системе единиц, где истинная энергия атома водорода равна 1/2, и перевели ее в электрон-вольты. Результат не очень точный (правильный, как помним, равен 13.6 eV), и мы его улучшим после того как разберемся, что и как только что посчитали. Это будет удобнее сделать на примере атома гелия, потому что в атоме водорода есть только один электрон, а в любой реальной системе – больше одного.

Решение на psi4#

Такой же результат можно получить и в библиотеке psi4. Правда в данном случае Python API выглядит немного иначе.

import psi4
psi4.core.be_quiet() # отключаем логирование в stdout

h_atom = psi4.geometry("H")

psi4.set_options({
  "basis": "STO-3G",
  "reference": "rohf",
})

from scipy.constants import physical_constants

h2ev = physical_constants["hartree-electron volt relationship"]

def e_in_ev(energy_in_ht):
    return energy_in_ht * h2ev[0]

# энергия в единицах Hartree
e_in_ht = psi4.energy("scf")

print(f"Hydrogen ground state energy: {e_in_ev(e_in_ht)} eV")

Hydrogen ground state energy: -12.696338923670487 eV

Теория самосогласованного поля#

Атом гелия#

Следующим “по простоте” после атома водорода идет атом гелия – как говорит нам школьная химия, это атом из двух протонов, двух нейтронов и двух электронов. Протоны и нейтроны находятся близко друг к другу в ядре и имеют почти одинаковую массу, так что можно просто считать, что есть ядро с зарядом +2 и массой 4. А вот с электронами все сложнее: с одной стороны, это независимые частицы, а с другой – они взаимодействуют друг с другом по закону Кулона, так как оба имеют отрицательный заряд.

Note

Для гравитационного взаимодействия “проблема трех тел” не имеет известного аналитического решения. Это означает, что если мы знаем, что где-то в глубоком космосе вдалеке от остального мира есть три объекта с известными массами, импульсами и координатами, то, увы, в общем случае не сможем предсказать их движение аналитически (хотя сможем предсказать численно, либо найти приближенное аналитическое решение, если масса одного объекта много больше других, например).

Для трех классических тел с кулоновским потенциалом все тоже сложно – можно посмотреть тут, как поведет себя система трех тел с различными зарядами.

Запишем уравнение Шредингера для системы из ядра и двух электронов:

\[ \dot{\imath} \hbar \frac{d}{dt} {\ket{\Psi(t, \vec{r_A}, \vec{r_B})}} = ( \frac{ \hat{p_A}^{2} }{2m} - \frac{2e \cdot e}{r_A} + \frac{ \hat{p_B}^{2} }{2m} - \frac{2e \cdot e}{r_B} + \frac{e \cdot e}{|r_A - r_B|} ) {\ket{\Psi(t, \vec{r_A}, \vec{r_B})}} \]

Оно состоит из кинетической и потенциальной энергии двух электронов в поле ядра и (последнее слагаемое) потенциала взаимодействия двух электронов. Учитывая, что уравнение примерно такое же, как у водорода, только еще более сложное, решение стоит искать численно и с использованием приближений, а не аналитически.

Если бы не последний член гамильтониана, то можно было бы разбить все выражение на две независимых части – с переменными электрона \(А\) и с переменными электрона \(B\). Так как это дифференциальное уравнение, можно было бы воспользоваться разделением переменных и найти отдельные решения для двух электронов – задача аналогична атому водорода, а ее мы уже решили.

Но из-за потенциала взаимодействия решение существенно усложняется, поскольку электроны влияют друг на друга. Придется прибегнуть к упрощениям – и одним из наиболее популярных подходов является теория самосогласованного поля (Self-Consistent Field).

Теория самосогласованного поля#

Теория самосогласованного поля (self-consistent field theory) – это подход итеративного решения уравнения Шредингера для многочастичной системы, на основе которого построено много квантово-химических методов, наиболее известный из которых – метод Хартри-Фока. В сниппете выше строчка scf.ROHF(h_atom) (psi4.energy('scf')) означает, что энергия посчитана этим методом.

Основная идея теории заключается в следующем.

Есть несколько частиц, которые взаимодействуют между собой, и найти цельное решение уравнения Шредингера для всех сразу не получается. Тогда вместо этого будем рассматривать частицы по очереди и считать, что все остальные действуют “в среднем” на выбранную частицу. То есть будем считать усредненный по пространству потенциал вместо точного.

Для электрона \(А\) в атоме гелия нам нужно учесть усредненное влияние электрона \(B\). Можем для электрона \(B\) взять волновую функцию от атома водорода, посчитать на ее основе усредненное влияние на электрон \(А\):

\[ (h_A + \hat{V}_{eff}) \ket{\psi_A} = (h_A + \bra{\psi_B} \frac{e^2}{r_{AB}} \ket{\psi_B}) \ket{\psi_A} = \epsilon_A \ket{\psi_A}, \]

где \(h_A\) – это кинетическая энергия и потенциал ядра для электрона \(A\), \(V_{eff}\) – влияние электрона \(B\) на электрон \(А\) и вычисляется как влияние усредненной электронной плотности, распределенной в соответствии с \(\psi_B\).

На этом шаге считаем, что \(\psi_B\) – фиксированная волновая функция и находим “переменную” \(\psi_A\). Однако \(\psi_A\) тоже будет влиять на \(\psi_B\) и, записав аналогичное уравнение для частицы \(B\), следующим шагом найдем новую \(\psi_B\).

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

Все вместе это создает итеративную процедуру:

  1. на основе имеющихся волновых функций посчитать среднее поле, которое создают частицы (например, электрон \(B\) для электрона \(А\));

  2. решить уравнение Шредингера с потенциалом, учитывающим среднее поле, то есть вычислить энергии и новые волновые функции;

  3. вернуться к шагу 1.

Для инициализации можно взять какие-то приблизительные волновые функции для всех частиц (двух электронов гелия), например, решения из уравнения водорода, то есть водородоподобные волновые функции. Далее итерацию повторяют, пока волновые функции и их энергии не перестанут изменяться, то есть самосогласуются.

Система таких уравнений, записанных для каждого электрона, называется уравнениями Хартри. На основе таких волновых функций для отдельных электронов можно собрать общую волновую функцию \(\Psi(\vec{r_A}, \vec{r_B})\), самый простой вариант – это \(\Psi(\vec{r_A}, \vec{r_B}) = \psi_A(\vec{r_A}) \psi_B(\vec{r_B})\), он и был предложен первоначально. Однако есть проблема: такая \(\Psi\) получается не антисимметричной, а только такие волновые функции для системы электронов являются “физичными”, то есть могут существовать в реальности.

Note

Один из основополагающих принципов квантовой механики – тождественность частиц и их неразличимость. Принцип означает, что все частицы одного типа (например, электроны) одинаковы и характеризуются только своим состоянием. Например, если мы “переставим местами” два электрона в атоме (не только пространственно, но и в смысле их состояния и энергий), то получившийся атом будет неотличим от исходного.

Из принципа тождественности следует, что при перестановке не должна меняться плотность вероятности:

\[ | {\Psi}(X_A, X_B) |^2 = | {\Psi}(X_B, X_A) |^2 \]

Кроме того, представим что мы обменяли частицы дважды: \(A\) с \(B\) и обратно. Никакие физические свойства системы из-за этого измениться не должны. В общем случае, из этого не следует что волновая функция не изменилась – например, так как глобальная фаза волновой функции неизмерима, то могли получить \(\Psi' = e^{2 i \pi \theta} \Psi\).

Но, если работаем больше чем в двух пространственных измерениях, то такой двойной обмен эквивалентен отсутсвию обмена. Для начала представим, что обмениваем частицы медленно, описывая полукруг частицей \(A\) вокруг частицы \(B\) и потом сдвигая обе частицы. Тогда двойной обмен значит, что частица \(A\) описывает полный круг вокруг частицы \(B\). Также, если непрерывно изменим ее маршрут не приближая ее к частице \(B\), то ожидаем получить тот же результат.

Итак, какие замкнутые маршруты можем получить непрерывно деформируя маршрут в евклидовом пространстве без точки (\(B\))? Ответ для любого измерения больше 2 – любые. Следовательно двойной обмен должен давать такой же результат как и если бы просто оставили частицу \(А\) на месте [NSS+08]. Это свойство задает ограничение на то, какими могут быть волновые функции.

У этого ограничения есть два решения: либо \({\Psi}(X_A, X_B) = {\Psi}(X_B, X_A)\) (симметричность), либо \(\Psi(X_A, X_B) = - \Psi(X_B, X_A)\) (антисимметричность). У антисимметричных функций есть интересное свойство: если функция \(f(x_1, x_2)\) антисимметрична, то \(f(x_1=X, x_2=X) = 0\), то есть антисимметричная функция равна нулю, если ее аргументы одинаковы. В этом легко убедиться на примере \(f = x_1 - x_2\).

Для антисимметричной волновой функции это означает, что две частицы не могут иметь полностью одинаковое состояние – волновая функциия (и вероятность) такой конфигурации равна нулю.

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

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

Чтобы сделать волновую функцию системы из двух электронов антисимметричной, используется такой прием:

\[ \Psi(X_A, X_B) = \psi_A(X_A) \psi_B(X_B) - \psi_B(X_A) \psi_A(X_B) \]

Во-первых, легко проверить, что если поменять местами \(X_A\) и \(X_B\), то вся функция просто изменит знак. Во-вторых, можно заметить, что формулу можно записать как определитель матрицы:

\[\begin{split} \Psi(X_A, X_B) = \begin{vmatrix} \psi_A(X_A) & \psi_B(X_A) \\ \psi_A(X_B) & \psi_B(X_B) \end{vmatrix} \end{split}\]

Из курса линейной алгебры можно вспомнить, что определитель меняет знак при перестановке двух столбцов или двух строк – это свойство позволяет делать антисимметричные волновые функции систем из \(N\) волновых функций отдельных электронов, если использовать метод Хартри не для атома гелия, а для системы с большим числом электронов. Для этого составляется определитель \(N \times N\) по аналогии с формулой выше: элемент в строке \(i\), столбце \(j\) – это \(i\)-я волновая функция с параметрами \(j\)-го электрона в качестве аргумента. Волновая функция системы частиц вычисляется как нормированный определитель матрицы, а определитель всегда антисимметричен. Такой вариант сборки волновой функции системы частиц называется “определитель Слэтера”.

Так как весь подход является аппроксимацией, не любая система может быть точно представлена таким детерминантом, но он является очень распространенным методом “сборки” волновой функции системы электронов в квантовой химии. Его использование также немного меняет вид одноэлектронных уравнений: чтобы корректно учесть антисимметрию, в эффективный потенциал добавляется так называемое “обменное взаимодействие”.

Все вместе составляет метод Хартри-Фока:

  • итеративная процедура самосогласованного поля;

  • усредненное действие электронов друг на друга, учет обменного взаимодействия;

  • детерминант Слетера.

SCF в python#

Теперь можно вернуться к коду и взглянуть на него чуть более осмысленно. При вычислении энергии явно передаем, что хотим посчитать ее методом Self-Consistent Field:

h_atom = gto.Mole()
h_atom.build(
    atom="H 0 0 0",
    basis="STO-3G",
    spin=1,
)

mf = scf.ROHF(h_atom)
e_in_ht = mf.kernel()

converged SCF energy = -0.466581849557275

Но что тут происходит, пока по-прежнему неясно.

Начнем с параметра basis. STO-3G – не стандарт связи, а Slater Type Orbital с 3 Гауссианами в базисном наборе, то есть базис на основе детерминанта Слетера. В описании метода SCF собирались начинать итерации с водородоподобных волновых функций, но так как весь метод является аппроксимирующим, никто не мешает выбрать другие волновые функции, если результаты лучше согласуются с экспериментом. Выбор базиса может существенно влиять на результат вычислений и современные базисы сложнее, чем Слетеровский детерминант – он просто один из первых и наиболее популярных.

Параметр spin=1 задает общий спиновый момент атома. Это очень важно, так как любые ограничения улучшают сходимость метода SCF. Ну а когда пишем scf.ROHF, то (в это случае) используем Restricted Open Shell Hartree-Fock, так как у атома водорода только один электрон и его оболочка не заполнена (на уровне энергии \(n = 1\) для этого нужно 2 электрона).

Повторим вычисления, например, с более “современным” базисом.

h_atom = gto.Mole()
h_atom.build(
    atom="H 0 0 0",
    basis="augccpv5z",
    spin=1,
)

mf = scf.ROHF(h_atom)
e_in_ht = mf.kernel()
print(f"Better hydrogen ground state energy: {e_in_ev(e_in_ht)} eV")

converged SCF energy = -0.49999478458371

Better hydrogen ground state energy: -13.605551204286888 eV

Та-дам! Используя более прокачанные базисы, получили правильный ответ.

Посмотрим, что еще можно сделать с помощью self-consistent field.

SCF с использованием psi4#

Используем более “модный” базис в psi4 и также получим более точное решение для атома водорода:

psi4.core.clean()

h_atom = psi4.geometry("H")

psi4.set_options({
  "basis": "d-aug-cc-pv5z", # разбор этого базиса выходит за рамки этого интро
  "scf_type": "pk",
  "reference": "rohf",
})

e_in_ht = psi4.energy("scf")

print(f"Better hydrogen ground state energy: {e_in_ev(e_in_ht)} eV")

Better hydrogen ground state energy: -13.605551648965216 eV

Атом Гелия (численно)#

Раз разобрали SCF на примере атома гелия, то наверняка можно посчитать его энергию в pyscf.

he_atom = gto.Mole()
he_atom.build(
    atom="He 0 0 0",
    basis="STO-3G",
    spin=0,
)

mf = scf.ROHF(he_atom)
e_in_ht = mf.kernel()
print(f"Helium ground state energy: {e_in_ev(e_in_ht)} eV")

converged SCF energy = -2.80778395753997

Helium ground state energy: -76.403693763909 eV

Атом гелия в psi4#

Ну и то же самое в psi4:

psi4.core.clean()

he_atom = psi4.geometry("He")

psi4.set_options({
  "basis": "STO-3G",
  "reference": "rohf",
})

e_in_ht = psi4.energy("scf")

print(f"Helium ground state energy: {e_in_ev(e_in_ht)} eV")

Helium ground state energy: -76.403693763909 eV

Экспериментальное значение энергии атома гелия равно -79.0 eV.

Молекула водорода#

Пока рассматривали только атомы, но SCF можно использовать и для молекул – потенциалы становятся сложнее, электронов больше, но общая логика не меняется.

h_mol = gto.Mole()
h_mol.build(
    atom="H 0 0 0; H 0 0 0.74",
    basis="STO-3G",
    spin=0,
)

mf = scf.RHF(h_mol)
e_in_ht = mf.kernel()
print(f"Hydrogen ground state energy: {e_in_ev(e_in_ht)} eV")

converged SCF energy = -1.11675930739643

Hydrogen ground state energy: -30.388568857366177 eV

Здесь задали явно координаты обоих атомов водорода в молекуле и энергия электронов была высчитана в предположении, что ядра водородов неподвижны. Здесь расстояние в 0.74 Ангстрема взято из экспериментальных данных. Если бы задали неправильные координаты, то рассчитанная энергия окажется неверной. Точнее, она соответствовала бы нефизичной ситуации, когда неведомая сила “удерживает” ядра водорода на месте.

Молекула водорода в psi4#

Несложно сделать тоже самое и в psi4:

psi4.core.clean()

# задаем 2 атома водорода с явными координатами
h_mol = psi4.geometry("""
  H 0 0 0
  H 0 0 0.74
""")

psi4.set_options({
  "basis": "STO-3G",
  "reference": "rohf",
})

e_in_ht_h = psi4.energy("scf", molecule=h_mol)

print(f"Hydrogen ground state energy: {e_in_ev(e_in_ht_h)} eV")

Hydrogen ground state energy: -30.388568856869096 eV

Оптимизация геометрии молекулы#

В pyscf (а также и в psi4) есть встроенная возможность оптимизации молекулы: pyscf.geomopt.berny_solver.optimize (psi4.optimize). В этом случае оптимизатор не фиксирует положение ядер и возвращает минимальную возможную энергию с учетом вариации положения атомов. Вычисления с оптимизацией геометрии занимают значительно больше времени.

Note

В pyscf оптимизация выполняется при помощи сторонней библиотеки Pyberny, которая не является зависимостью и должна быть установлена отдельно: pip install -U pyberny.

from pyscf.geomopt.berny_solver import optimize

h_mol_bad = gto.Mole()
h_mol_bad.build(
    atom="H 0 0 0; H 0 0 1.5",
    basis="STO-3G",
    spin=0,
)

mf = scf.RHF(h_mol_bad)
e_in_ht_h_bad = mf.kernel()

mol_eq = optimize(mf)
mf_right = scf.RHF(mol_eq)
e_in_ht_h_optimized = mf_right.kernel()

print(f"Hydrogen molecule, incorrect ground state energy: {e_in_ev(e_in_ht_h_bad)} eV")
print(f"Hydrogen molecule, optimized ground state energy: {e_in_ev(e_in_ht_h_optimized)} eV")

converged SCF energy = -0.910873554594387

Geometry optimization cycle 1

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.000000    0.000000  0.000000  0.000000
   H   0.000000   0.000000   1.500000    0.000000  0.000000  0.000000
converged SCF energy = -0.910873554594387

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000    -0.1583976826
1 H     0.0000000000     0.0000000000     0.1583976826
----------------------------------------------

cycle 1: E = -0.910873554594  dE = -0.910874  norm(grad) = 0.224008

Geometry optimization cycle 2

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.079377    0.000000  0.000000  0.079377
   H   0.000000   0.000000   1.420623    0.000000  0.000000 -0.079377

WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.

converged SCF energy = -0.959976948556178

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000    -0.1678535564
1 H     0.0000000000     0.0000000000     0.1678535564
----------------------------------------------

cycle 2: E = -0.959976948556  dE = -0.0491034  norm(grad) = 0.237381

Geometry optimization cycle 3

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.238130    0.000000  0.000000  0.158753
   H   0.000000   0.000000   1.261870    0.000000  0.000000 -0.158753

WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.

converged SCF energy = -1.05934476925058

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000    -0.1528552490
1 H     0.0000000000     0.0000000000     0.1528552490
----------------------------------------------

cycle 3: E = -1.05934476925  dE = -0.0993678  norm(grad) = 0.21617

Geometry optimization cycle 4

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.481967    0.000000  0.000000  0.243837
   H   0.000000   0.000000   1.018033    0.000000  0.000000 -0.243837

WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.

converged SCF energy = -1.07060402057946

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000     0.3424691642
1 H     0.0000000000     0.0000000000    -0.3424691642
----------------------------------------------

cycle 4: E = -1.07060402058  dE = -0.0112593  norm(grad) = 0.484325

Geometry optimization cycle 5

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.331032    0.000000  0.000000 -0.150935
   H   0.000000   0.000000   1.168968    0.000000  0.000000  0.150935

WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.

converged SCF energy = -1.104755003306

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000    -0.0954526877
1 H     0.0000000000     0.0000000000     0.0954526877
----------------------------------------------

cycle 5: E = -1.10475500331  dE = -0.034151  norm(grad) = 0.13499

Geometry optimization cycle 6

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.364363    0.000000  0.000000  0.033331
   H   0.000000   0.000000   1.135637    0.000000  0.000000 -0.033331

WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.

converged SCF energy = -1.11432815053439

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000    -0.0537746620
1 H     0.0000000000     0.0000000000     0.0537746620
----------------------------------------------

cycle 6: E = -1.11432815053  dE = -0.00957315  norm(grad) = 0.0760489

Geometry optimization cycle 7

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.407339    0.000000  0.000000  0.042976
   H   0.000000   0.000000   1.092661    0.000000  0.000000 -0.042976

WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.

converged SCF energy = -1.11672330760192

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000     0.0316379732
1 H     0.0000000000     0.0000000000    -0.0316379732
----------------------------------------------

cycle 7: E = -1.1167233076  dE = -0.00239516  norm(grad) = 0.0447429

Geometry optimization cycle 8

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.391420    0.000000  0.000000 -0.015919
   H   0.000000   0.000000   1.108580    0.000000  0.000000  0.015919

WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.

/home/runner/work/qmlcourse/qmlcourse/.venv/lib/python3.8/site-packages/pyscf/dft/libxc.py:771: UserWarning: Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, corresponding to the original definition by Stephens et al. (issue 1480) and the same as the B3LYP functional in Gaussian. To restore the VWN5 definition, you can put the setting "B3LYP_WITH_VWN5 = True" in pyscf_conf.py
  warnings.warn('Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, '

converged SCF energy = -1.11748127474604

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000    -0.0052580454
1 H     0.0000000000     0.0000000000     0.0052580454
----------------------------------------------

cycle 8: E = -1.11748127475  dE = -0.000757967  norm(grad) = 0.007436

Geometry optimization cycle 9

Cartesian coordinates (Angstrom)
 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.393689    0.000000  0.000000  0.002269
   H   0.000000   0.000000   1.106311    0.000000  0.000000 -0.002269

WARN: Large deviations found between the input molecule and the molecule from chkfile
Initial guess density matrix may have large error.

converged SCF energy = -1.11750572804771

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000    -0.0004239054
1 H     0.0000000000     0.0000000000     0.0004239054
----------------------------------------------

cycle 9: E = -1.11750572805  dE = -2.44533e-05  norm(grad) = 0.000599493

Geometry optimization cycle 10

Cartesian coordinates (Angstrom)

 Atom        New coordinates             dX        dY        dZ
   H   0.000000   0.000000   0.393888    0.000000  0.000000  0.000199
   H   0.000000   0.000000   1.106112    0.000000  0.000000 -0.000199

converged SCF energy = -1.11750588508088

--------------- RHF_Scanner gradients ---------------

         x                y                z
0 H     0.0000000000     0.0000000000     0.0000063561
1 H     0.0000000000     0.0000000000    -0.0000063561
----------------------------------------------

cycle 10: E = -1.11750588508  dE = -1.57033e-07  norm(grad) = 8.9889e-06

converged SCF energy = -1.11750588508088

Hydrogen molecule, incorrect ground state energy: -24.78613211532389 eV

Hydrogen molecule, optimized ground state energy: -30.40888427110039 eV

Для некорректной геометрии получилась завышенная энергия, а после оптимизации – почти что такая же энергия, как при вычислении с фиксированным расстоянием 0.74. В оптимальном состоянии энергия системы должна быть минимальна, так что результаты вполне разумны.

Подобный метод можно использовать и для поиска геометрии куда более сложных молекул.

Оптимизация в psi4#

psi4.core.clean()

h_mol_bad = psi4.geometry("""
  H 0 0 0
  H 0 0 1.5
""")
# 1.5 - неверное расстояние в ангстремах, верное - 0.74

psi4.set_options({
  "basis": "STO-3G",
  "reference": "rohf",
})

# рассчитываем энергию "в точке" с неправильной геометрией
e_in_ht_h_bad = psi4.energy("scf", molecule=h_mol_bad)

# рассчитываем энергию, оптимизируя по ходу геометрию
e_in_ht_h_optimized = psi4.optimize("scf", molecule=h_mol_bad)

print(f"Hydrogen molecule, incorrect ground state energy: {e_in_ev(e_in_ht_h_bad)} eV")
print(f"Hydrogen molecule, optimized ground state energy: {e_in_ev(e_in_ht_h_optimized)} eV")

Optimizer: Optimization complete!
Hydrogen molecule, incorrect ground state energy: -24.786132109551886 eV
Hydrogen molecule, optimized ground state energy: -30.408884269506693 eV

Молекула этилового спирта#

Молекула водорода – это все еще почти игрушечный пример. Давайте попробуем обсчитать молекулу этанола.

Задавать руками геометрию молекулы \(C_2H_5OH\) можно, но будет явно сложнее, чем для молекулы водорода. К счастью, это необязательно: psi4 умеет скачивать геометрию из базы данных PubChem по номенклатурному имени либо уникальному ChemId.

Note

Возможность поиска геометрии в PubChem – это уникальная “фишка” psi4. В pysc, в данном случае, пришлось бы искать молекулу в базе вручную, а потом самому забивать в код все координаты.

psi4.core.clean()

eth = psi4.geometry("pubchem:ethanol")

psi4.set_options({
  "basis": "STO-3G",
  "reference": "rohf",
})

e_in_ht_eth = psi4.energy("scf", molecule=eth)

print(f"Ethanol ground state energy: {e_in_ev(e_in_ht_eth)} eV")

	Searching PubChem database for ethanol (single best match returned)
	Found 1 result(s)
Ethanol ground state energy: -4139.64590826589 eV

Что узнали?#

Разобрались с базовой теорией квантовой химии:

  • как записать уравнение Шредингера для атома;

  • какое получается аналитическое решение для атома водорода;

  • как устроен метод Self-Consistent Field для вычисления волновых функций и энергии для задачи многих тел;

  • как использовать SCF в python пакетах pyscf и psi4.

В примерах везде вычисляли ground state энергию, но, конечно, зная волновые функции, можно посчитать много чего еще. Например, можно вычислить спектр поглощения и энергию ионизации (энергии возбужденных состояний), моделировать взаимодействие молекул (найти равновесное состояние для двух систем), с помощью плагинов можно смоделировать рассеяние рентгена на молекуле… Квантовая механика постулирует, что любая измеримая величина является усреднением определенного оператора по волновой функции, поэтому в теории любую физическую величину можно посчитать, зная волновую функцию. Возможности ограничены в основном вычислительной сложностью, а не теорией.

Для более глубокого погружения в практику квантовой химии можно пройти лабораторные работы psi4: раз, два.