Программная реализация метода функционала плотности для взаимодействующего электронного газа на однородном компенсирующем фоне

Алексей Кожарин

Программная реализация метода функционала плотности для взаимодействующего электронного газа на однородном компенсирующем фоне

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

Физика

Дипломы

Вуз: Московский физико-технический институт (государственный университет) (МФТИ)

ID: 5f2abe00cd3d3e0001b3bd94

UUID: 70eb8140-b953-0138-11d7-0242ac180006

Язык: Русский

Опубликовано: больше 4 лет назад

Просмотры: 13

10.26

Алексей Кожарин

Московский физико-технический институт (государственный университет) (МФТИ)

Комментировать 12

Рецензировать 0

Скачать - 708,8 КБ

Поделиться работой

Московский физико-технический институт Работа допущена к защите зав. кафедрой В.Е. Фортов « » 2020 г. ВЫПУСКНАЯ РАБОТА БАКАЛАВРА Тема: Программная реализация метода функционала плотности для взаимодействующего электронного газа на однородном компенсирующем фоне Направление: 03.03.01 – Прикладная математика и физика Выполнил студент гр. 642б Кожарин Алексей Сергеевич Научный руководитель, Левашов Павел Ремирович к. ф.-м. н. Москва – 2020

2 Оглавление Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Глава 1. Идеальный ферми–газ . . . . . . . . . . . . . . . . . . . . 6 1.1. Химический потенциал и его производные . . . . . . . . . . . 6 1.2. Термодинамические функции ИФГ . . . . . . . . . . . . . . . 8 1.3. Асимптотики при низких температурах . . . . . . . . . . . . . 9 1.4. Асимптотики при высоких температурах . . . . . . . . . . . . 11 1.5. Программная реализация модели ИФГ . . . . . . . . . . . . . 13 1.6. Результаты расчетов . . . . . . . . . . . . . . . . . . . . . . . 14 Глава 2. Модель взаимодействующего электронного газа на одно родном компенсирующем фоне . . . . . . . . . . . . . . . . . . . . 16 2.1. Классическая модель однокомпонентной плазмы . . . . . . . 16 2.2. Суммирование Эвальда . . . . . . . . . . . . . . . . . . . . . . 18 2.3. Гамильтониан модели ВЭГ . . . . . . . . . . . . . . . . . . . . 21 2.4. Усреднение по направлениям и изотропная модель ВЭГ . . . 26 Глава 3. Программная реализация метода функционала плотности для модели ВЭГ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.1. Метод функционала электронной плотности . . . . . . . . . . 30 3.2. Метод функционала плотности для модели ВЭГ . . . . . . . . 33 3.3. Алгоритм моделирования ВЭГ методом функционала плот ности при конечной температуре . . . . . . . . . . . . . . . . 35 Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Словарь терминов . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3 Приложение А. Пример расчета термодинамических свойств ИФГ . . . . . . . . . . . . . . . . . . 42

4 Введение Свойства электронного газа необходимо знать для решения широкого круга задач астрофизики, физики твердого тела, физики плазмы, химии и мо лекулярной биологии; кроме того, модель электронного газа играет важную роль в статистической физике. Учет квантовых свойств электронов позво лил успешно объяснить такие явления, как эволюция и коллапс звезд, взрыв новых и сверхновых, проводимость полупроводников и металлов, сверхпро водимость и многие другие. Широкую известность получила модель невза имодействующего (идеального) ферми-газа в силу своей простоты. Между тем, учет взаимодействия для системы квантовых частиц представляет со бой чрезвычайно сложную задачу статистической физики, не решенную до настоящего времени. Особые сложности представляет кулоновское взаимо действие, так как в этом случае необходимо учитывать дальнодействующий характер потенциала. Поэтому особое значение приобретают модельные си стемы, на примере которых можно предсказывать свойства реальных веществ. Наиболее популярной моделью такого рода является модель «желе» — элек тронов на однородном несжимаемом компенсирующем положительном фоне. Эта модель позволяет в чистом виде вычислить обменно–корреляционную энергию, которая необходима для построения так называемых обменно–кор реляционных функционалов для метода функционала плотности (МФП). В настоящее время МФП очень активно применяется для моделирования раз личных свойств реальных веществ и играет большое значение в квантовой химии, физике твердого тела и физике неидеальной плазмы. Модель «желе» хорошо изучена при нулевой температуре, однако исследования этой модели при конечных температурах активно продолжаются и в наше время, и на се годняшний день далеки от завершения. Результатом изучения модели «желе» при конечной температуре должно стать создание надежных обменно–кор реляционных функционалов с явной зависимостью от температуры. Большой

5 интерес представляют также и другие свойства модели «желе», в частно сти, диэлектрическая функция. Поэтому весьма актуальной является про граммная реализация метода функционала плотности для взаимодействую щего электронного газа, позволяющая, помимо термодинамических свойств, рассчитать также транспортные и оптические свойства. Цель работы — получение аналитических выражений для термодинами ческих функций и термодинамических коэффициентов модели идеального ферми–газа, а также разработка алгоритма программной реализации мето да функционала плотности для взаимодействующего электронного газа на однородном компенсирующем фоне положительного заряда.

6 Глава 1 Идеальный ферми–газ Идеальный Ферми-газ (ИФГ) – это система невзаимодействующих фер мионов (частиц, имеющих полуцелый спин). Согласно принципу запрета Па ули [1; 2], в заданном квантовом состоянии может находиться только один фермион. Это приводит приводит к нескольким следствиям. Во-первых, пол ная энергия ИФГ при нулевой температуре отлична от нуля. Во-вторых, дав ление ИФГ не равно нулю при нулевой температуре, в отличие от идеального больцмановского газа. В дальнейшей части работы, если не оговорено иное, используется атом ная система единиц, в которой редуцированная постоянная Планка, масса и заряд фермиона равны единице. Помимо этого, постоянная Больцмана тоже положена равной единице. 1.1. Химический потенциал и его производные Для системы фермионов распределение частиц по энергиям задается статистикой Ферми–Дирака [3]: {eq:fd} 𝑛𝑘 = 1 . exp{[(𝜀𝑘 − 𝜇)/𝑇 ]} + 1 (1.1) Здесь 𝜇 есть химический потенциал, 𝑇 – температура системы, а 𝜀𝑘 – энергия 𝑘-ого квантового состояния. Полное число частиц получается суммировани ем (1.1) по всем возможным значениям энергии. Полагая 𝜀𝑘 = 𝑝2𝑘 /(2𝑚) и переходя к квазинепрерывному спектру по 𝑘, после сведения к интегрирова нию по энергии получим: {eq:conc} 𝑁 𝑔 =√ 𝑉 2𝜋 2 ∫︁∞ 0 √ 𝜀𝑑𝜀 . 𝑒(𝜀−𝜇)/𝑇 + 1 (1.2)

7 Здесь 𝑉 — объем системы, 𝑔 — фактор вырождения по спину (для электрона 𝑔 = 2). Вводя безразмерный химический потенциал 𝑦 = 𝜇/𝑇 и объем, при ходящийся на одну частицу 𝑣 = 𝑉 /𝑁 , можно переписать (1.2) в следующем виде: mu_equation} _definition} 1 𝑔𝑇 3/2 =√ 𝐼1/2 (𝑦). 𝑣 2𝜋 2 Здесь введен так называемый интеграл Ферми-Дирака: ∫︁∞ 𝑡𝑗 𝐼𝑗 (𝑦) = 𝑑𝑡, 𝑒𝑡−𝑦 + 1 (1.3) (1.4) 0 обладающий следующим важным свойством при дифференцировании по па раметру: g_diff_rule} 𝑑 𝐼𝑘 (𝑥) = 𝑘𝐼𝑘−1 (𝑥) (1.5) 𝑑𝑥 Уравнение (1.3) задает химический потенциал как функцию от темпе ратуры и удельного объема: 𝜇 = 𝜇(𝑣, 𝑇 ). В разделах 1.3, 1.4 приводится его решение в пределе низких и высоких температур. В дальнейшем нам понадобятся первые и вторые производные 𝜇(𝑣, 𝑇 ). Для их получения продифференцируем (1.3) по соответствующим перемен ным с учетом свойства (1.5) и разрешим получающиеся уравнения относи тельно производных. Это даст следующие результаты: {eq:mu_T} 𝑦𝑇′ = − 3𝐼1/2 (𝑦) 𝑇 𝐼−1/2 (𝑦) (1.6) {eq:mu_v} 𝑦𝑣′ = − 2𝐼1/2 (𝑦) 𝑣𝐼−1/2 (𝑦) (1.7) ′′ 𝑦𝑣𝑣 {eq:mu_vv} {eq:mu_TT} {eq:mu_vT} 𝑦𝑇′′ 𝑇 ′′ 𝑦𝑣𝑇 = 2 2𝐼−3/2 (𝑦)𝐼1/2 (𝑦) 4𝐼1/2 (𝑦) = 2 + 3 𝑣 𝐼−1/2 (𝑦) 𝑣 2 𝐼−1/2 (𝑦) (1.8) 2 9𝐼−3/2 (𝑦)𝐼1/2 (𝑦) 15𝐼1/2 (𝑦) = + 3 2𝑇 2 𝐼−1/2 (𝑦) 2𝑇 2 𝐼−1/2 (𝑦) (1.9) 𝑦𝑇′′ 𝑣 2 3𝐼−3/2 (𝑦)𝐼1/2 (𝑦) 3𝐼1/2 (𝑦) = + 3 (𝑦) 𝑇 𝑣𝐼−1/2 (𝑦) 𝑇 𝑣𝐼−1/2 (1.10)

8 1.2. Термодинамические функции ИФГ Согласно [4], выражение для свободной энергии Гельмгольца ИФГ имеет вид: z_potential} (︂ )︂ [︂ ]︂ 𝑔 2 2 𝐼3/2 (𝑦) 5/2 𝐹 =√ 𝑇 𝑣 𝑦 𝐼1/2 (𝑦) − 𝐼3/2 (𝑦) = 𝑇 𝑦 − . 3 3 𝐼1/2 (𝑦) 2𝜋 2 (1.11) Последнее равенство получается подстановкой 𝑣 из (1.3). Из выражения для 𝐹 можно получить остальные термодинамические величины: 𝑃 = −𝐹𝑣′ ; {eq:PES} 𝐶𝑉 = −𝑇 𝐹𝑇′′𝑇 ; {eq:CVCP} {eq:CTCS} 𝐸 = 𝐹 − 𝑇 𝐹𝑇′ ; 𝐶𝑇2 =𝑣 2 ′′ 𝐹𝑣𝑣 ; 𝐶𝑆2 =𝑣 2 ′′ 𝐹𝑣𝑣 𝑆 = −𝐹𝑇′ ; 𝐶𝑃 = −𝑇 𝐹𝑇′′𝑇 + ′′ 2 ) 𝑇 (𝐹𝑣𝑇 ; ′′ 𝐹𝑣𝑣 ′′ 𝑉 𝐹𝑣𝑇 𝛾= ; ... 𝑇 𝐹𝑇′′𝑇 ′′ 2 𝑣 2 (𝐹𝑣𝑇 ) ; − 𝐹𝑇′′𝑇 (1.12) (1.13) (1.14) Здесь 𝑃 – давление, 𝐶𝑉 и 𝐶𝑃 – теплоемкости при постоянном объеме и дав лении, соответственно, 𝐶𝑇 и 𝐶𝑆 – изотермическая и адиабатическая скорости звука, 𝛾 – параметр Грюнайзена. Используя формулы (1.6) – (1.10), получим для термодинамических функций ИФГ: eq:pressure} {eq:energy} {eq:entropy} (1.15) 𝐼3/2 (𝑦) 𝑔𝑣 5/2 𝐸=√ 𝑇 𝐼3/2 (𝑦) = 𝑇, 𝐼1/2 (𝑦) 2𝜋 2 (1.16) √ (︀ )︀ 𝑔 2𝑇 3/2 𝑣 3𝐼1/2 (𝑦)𝑦 − 5𝐼3/2 (𝑦) 3𝐼1/2 (𝑦) − 5𝐼3/2 (𝑦) 𝑆=− = − , 6𝜋 2 3𝐼1/2 (𝑦) √ 𝑔 2𝑇 𝐶𝑣 = {eq:C_v} √ 𝑔 2𝑇 5/2 𝐼3/2 (𝑦) 2𝑇 𝐼3/2 (𝑦) 𝑃 = = , 3𝜋 2 3𝑣 𝐼1/2 (𝑦) 3/2 (︁ 𝑣 5𝐼−1/2 (𝑦)𝐼3/2 (𝑦) − 4𝜋 2 𝐼−1/2 (𝑦) )︁ 2 9𝐼1/2 (𝑦) = (1.17) 5 𝐼3/2 (𝑦) 9 𝐼1/2 (𝑦) − , 2 𝐼1/2 (𝑦) 2 𝐼−1/2 (𝑦) (1.18)

9 √ 5𝑔 2𝑇 𝐶𝑃 = 3/2 (︁ 𝑣 5𝐼−1/2 (𝑦)𝐼3/2 (𝑦) − 2 (𝑦) 36𝜋 2 𝐼1/2 {eq:C_S} 𝐼3/2 (𝑦) = 2 25 𝐼3/2 (𝑦)𝐼−1/2 (𝑦) 5 𝐼3/2 (𝑦) − , (1.19) 3 (𝑦) 18 𝐼1/2 2 𝐼1/2 (𝑦) {eq:C_P} {eq:C_T} )︁ 2 9𝐼1/2 (𝑦) 𝐶𝑇2 = √ 2 2𝑔 𝑇 5/2 𝑣𝐼1/2 (𝑦) 𝜋2 √ 5 2𝑔 𝐼−1/2 (𝑦) = 2𝑇 𝐼1/2 (𝑦) , 𝐼−1/2 (𝑦) (1.20) 10𝑇 𝐼3/2 (𝑦) . (1.21) 9𝜋 2 9𝐼1/2 (𝑦) В разделах 1.3 и 1.4 обсуждается асимптотическое поведение данных 𝐶𝑆2 = 𝑇 5/2 𝑣𝐼3/2 (𝑦) = величин в пределе низких и высоких температур. 1.3. Асимптотики при низких температурах c:asymp_low} При низких температурах выполняется условие 𝑦 = 𝜇/𝑇 ≫ 1. В этом пределе справедливо разложение [5]: [︂ ]︂ 𝜋2 7𝜋 4 2𝑦 3/2 {eq:Ihalf} 1+ 2 + + ... , 𝐼1/2 (𝑦) ≈ 3 8𝑦 640𝑦 4 (1.22) [︂ ]︂ 2𝑦 5/2 5𝜋 2 7𝜋 4 {eq:I3half} 𝐼3/2 (𝑦) ≈ 1+ 2 − + ... . (1.23) 5 8𝑦 384𝑦 4 Если оставить только первый член в (1.22), то выражение (1.3) станет незави q:mulowtemp} симым от температуры. Это соответствует случаю 𝑇 = 0, в котором химиче ский потенциал равен энергии Ферми: 𝜇|𝑇 =0 = 𝜀𝐹 = (︂ 3𝜋 2 √ 2𝑔𝑣 )︂2/3 . (1.24) Если же оставить первые два члена в (1.22) и подставить 𝑦 = 𝜀𝐹 /𝑇 во второй из них, то можно получить первую поправку для химического потенциала: [︃ (︂ )︂2 ]︃ 2 𝜋 𝑇 𝜇 ≈ 𝜀𝐹 1 − . (1.25) 12 𝜀𝐹

10 Подставляя (1.22), (1.23) в (1.11) и оставляя первые два члена, получим: [︂ [︂ ]︂ ]︂ 2 5/2 2 𝑔 𝜋 3𝑇 2 𝜋 2 𝐹 ≈√ eq:Flowtemp} 𝑇 5/2 𝑣 𝑦 5/2 − 𝑦 1/2 = 3/2 𝑦 5/2 − 𝑦 1/2 . (1.26) 2 5 12 5 12 2𝜋 2𝜀𝐹 Наконец, подстановка (1.25) в (1.26) даст: [︃ (︂ )︂2 ]︃ 2 3 5𝜋 𝑇 𝛽 {eq:F_lt} 𝐹 ≈ 𝜀𝐹 1 − = 𝐴𝑣 −2/3 − 𝑇 2 𝑣 2/3 . 5 12 𝜀𝐹 2 (1.27) Здесь 𝛽 = (𝑔𝜋/6)3/2 — так называемый коэффициент электронной теплоем кости [3], {eq:A_def} )︂2/3 3𝜋 2 √ . (1.28) 2𝑔 Используя формулы (1.12), (1.13), (1.14), получим в пределе низких тем 3 𝐴= 5 (︂ ператур: {eq:P_lt} 2 𝛽 𝑃 = 𝐴𝑣 −5/3 + 𝑇 2 𝑣 −1/3 , 3 3 (1.29) {eq:E_lt} 3 𝛽 𝛽 𝐸 = 𝐴𝑣 −2/3 + 𝑇 2 𝑣 2/3 = 𝜀𝐹 + 𝑇 2 𝑣 2/3 , 2 5 2 (1.30) {eq:S_lt} 𝑆 = 𝛽𝑇 𝑣 2/3 , (1.31) {eq:cv_lt} 𝐶𝑣 = 𝛽𝑇 𝑣 2/3 , (1.32) {eq:cp_lt} 𝐶𝑃 = 𝛽𝑇 𝑣 2/3 , (1.33) {eq:ct2_lt} 𝐶𝑇2 = 10𝐴 −2/3 𝛽 2 2/3 𝑣 + 𝑇 𝑣 , 9 9 (1.34) {eq:cs2_lt} 𝐶𝑆2 = 10𝐴 −2/3 5𝛽 2 2/3 𝑣 + 𝑇 𝑣 . 9 9 (1.35)

11 1.4. Асимптотики при высоких температурах :asymp_high} При высоких температурах реализуется случай 𝑦 = 𝜇/𝑇 ≪ −1. По скольку на высоких температурах ИФГ переходит в идеальный больцманов ский газ (ИБГ), обсуждаемые свойства ИФГ должны стремиться к классиче скому результату при 𝑇 → ∞. Чтобы это показать, необходимо сделать два замечания. Во-первых, при 𝑦 ≪ −1 можно пренебречь единицей в знаменателе (1.4) и получить: 𝐼𝑘 (𝑦) ≈ Γ(𝑘 + 1)𝑒𝑦 , fg_hightemp} mu_boltzman} (1.36) Здесь Γ(𝑥) – гамма-функция. В частности, 𝐼−1/2 (𝑦) ≈ √ 𝜋𝑒𝑦 ; √ 3 𝜋 𝑦 𝐼3/2 (𝑦) ≈ 𝑒 . 4 √ 𝜋 𝑦 𝐼1/2 (𝑦) ≈ 𝑒 ; 2 Во-вторых, правая часть (1.3) должна оставаться конечной при 𝑇 → ∞. Поэтому 𝐼1/2 (𝑦) должен стремиться к нулю или, что равносильно, 𝜇/𝑇 → −∞. Это означает, что 𝜇 стремится к −∞, причем быстрее, чем 𝑇 в первой степени. Подстановка (1.36) в (1.3) приводит к классическому результату для идеального больцмановского газа: [︃ 1 𝜇 = 𝜇𝐵 = 𝑇 ln 𝑔𝑣 (︂ 2𝜋 𝑇 )︂3/2 ]︃ (1.37) что, очевидно, удовлетворяет вышеупомянутым требованиям. Подстановка разложения (1.36) в выражения для энергии Гельмгольца (1.11), давления (1.15), энергии (1.16), энтропии (1.17), теплоемкостей (1.18), (1.19) и скоростей звука (1.20), (1.21) дает: [︃ 𝐹 ≈ (𝑦 − 1)𝑇 = 𝜇 − 𝑇 = −𝑇 ln 𝑔𝑒𝑣 𝑃 ≈ 𝑇 , 𝑣 𝐸≈ 3 𝑇, 2 𝑆≈ (︂ 𝑇 2𝜋 5 − 𝑦, 2 )︂3/2 ]︃ , (1.38) (1.39)

12 3 𝐶𝑉 ≈ , 2 5 𝐶𝑃 ≈ , 2 𝐶𝑇2 ≈ 𝑇, (1.40) 5 𝑇. 3 𝐶𝑆2 ≈ (1.41) Полученные результаты соответствуют идеальному больцмановскому газу. Чтобы получить первую поправку к ним, необходимо более точное раз ложение (1.4): 𝑒2𝑦 𝐼𝑗 (𝑦) ≈ 𝑒 Γ(𝑗 + 1) − 𝑗+1 Γ(𝑗 + 1), 2 𝑦 emp_precise} (1.42) В частности, 𝐼1/2 (𝑦) ≈ √ 𝜋𝑒𝑦 2 (︂ )︂ )︂ √ (︂ 1 𝑦 1 𝑦 3 𝜋 1 − 3/2 𝑒 , 𝐼3/2 (𝑦) ≈ 1 − 5/2 𝑒 . 4 2 2 Подставляя (1.42) в уравнение для 𝜇 (1.3) и затем заменяя 𝑦 на 𝜀𝐹 /𝑇 во втором члене, получим: 1 𝑔𝑣 tion_high_T} (︂ 2𝜋 𝑇 )︂3/2 = exp (︁ 𝜇 )︁ 𝐵 𝑇 𝑦 =𝑒 (︂ 𝜋 3/2 1− 𝑔𝑣𝑇 3/2 )︂ . (1.43) Второй член в скобках мал по сравнению с единицей, поэтому можно записать: (︂ )︂ 3/2 𝜋 𝑒𝑦 ≈ 𝑒𝜇𝐵 /𝑇 1 + , 𝑔𝑣𝑇 3/2 {eq:muhighT} С той же точностью: (1.44) 𝜋 3/2 𝜇𝐵 + . 𝑦= 𝑇 𝑔𝑣𝑇 3/2 Подставим теперь (1.42) в общую формулу для давления (1.15), затем заменим 𝑒𝑦 в соответствии с (1.44) и оставим только члены 𝑇 и 𝑇 −1/2 Это даст: 𝑇 𝜋 3/2 𝑃 ≈ + √ {eq:p_ht} . (1.45) 𝑣 2 𝑇 𝑔𝑣 2 Заметим, что знак «+» перед вторым членом соответствует отталкиванию фермионов.

13 Следуя этому алгоритму, можно получить асимптотическое выражение для удельной энергии Гельмгольца (опять же, оставляя только члены с 𝑇 в степени, большей или равной −1/2): 𝜋 3/2 𝐹 = 𝜇 − 𝑃 𝑣 = (𝜇𝐵 − 𝑇 ) + . 2𝑔𝑣𝑇 1/2 {eq:FhighT} {eq:EShighT} q:CVCPhighT} q:CTCShighT} (1.46) Остальные асимптотические выражения получаются дифференцирова нием (1.46) в соответствии с общими выражениями (1.16), (1.17), (1.18), (1.19), (1.20), (1.21): 3/2 3𝑇 3𝜋 𝐸= + √ ; 2 4 𝑇 𝑔𝑣 [︃ 5 1 𝑆 = − ln 2 𝑔𝑣 3 3𝜋 3/2 𝐶𝑉 = − ; 2 8𝑇 3/2 𝑔𝑣 (︂ 2𝜋 𝑇 )︂3/2 ]︃ 𝜋 3/2 ; + 4𝑇 3/2 𝑔𝑣 15𝜋 3/2 5 𝐶𝑃 = − ; 2 8𝑇 3/2 𝑔𝑣 (1.47) (1.48) 5𝑇 𝜋 3/2 5𝜋 3/2 2 ; 𝐶𝑆 = . (1.49) =𝑇+√ + √ 3 𝑇 𝑔𝑣 6 𝑇 𝑔𝑣 Из полученных выражений видно, что эффекты вырождения приводят к 𝐶𝑇2 уменьшению теплоемкости и увеличению скорости звука. 1.5. Программная реализация модели ИФГ Формулы (1.15), (1.16), (1.17), (1.18), (1.19), (1.20), (1.21) дают возмож ность напрямую вычислять характеристики ИФГ через интеграл Ферми-Дира ка (1.4). Для языка Python существует модуль fdint, осуществляющий эф фективный подсчет упомянутого интеграла [6]. Формулы (1.17), (1.18), (1.19) содержат операцию вычитания и поэтому их использование в пределах высоких или низких температур будет приводить к большим ошибкам в расчете. Для решения этой проблемы в упомянутом диапазоне температур применяются асимптотические формулы для расчета соответствующих величин – (1.47), (1.31) для энтропии, (1.48), (1.32), (1.33) для теплоемкостей.

14 Расчет свойств был реализован в виде модуля на Python и выложен на PyPi под названием ifg [7]. Модуль поддерживает расчет химического по тенциала 𝜇, энергии 𝐸, энтропии 𝑆, теплоемостей 𝐶𝑣 , 𝐶𝑃 , а также скоростей звука 𝐶𝑇 , 𝐶𝑆 . Возможно использование как атомной системы единиц, так и СИ. Расчеты можно выполнять в широком диапазоне температур и объемов. Из-за упомянутых численных проблемы диапазон допустимого ввода ограничен: 𝑇 ∈ [10−49 , 1049 ] для температур и 𝑣 ∈ [10−30 , 1020 ] для удельных объемов. В этих диапазонах гарантируется точность 𝛿 = 10−7 (это было протестировано с использованием фреймворка Hypothesis [8]). 1.6. Результаты расчетов На графике 1.1 проиллюстрированы химический потенциал (а) и дав ление (б) на трех изохорах, полученные численно применением полученных выше формул с параметрами 𝑔 = 2, 𝑚𝑟 = 1. Видно, что химический потен циал стремится к 𝜀𝐹 при 𝑇 → 0 и быстро убывает в сторону отрицательных значений с повышением температуры. 102 104 103 0 µ 10 P , atomic units 101 l_potential} 10−1 10−2 10−3 −2 10 v = 0.1 v=1 v = 10 low-T asymp. high-T asymp. 10−1 101 v = 0.1 v=1 v = 10 low-T asymp. high-T asymp. 100 10−1 100 T , Ha (а) 102 101 102 10−2 −2 10 10−1 100 101 102 103 104 T , Ha (б) Рис. 1.1. Химический потенциал (a) и давление (б) ИФГ на трех изохорах. Теплоемкости при постоянном объеме и давлении показаны на графике 1.2. При низких температурах оба 𝐶𝑉 и 𝐶𝑃 линейно зависят от температуры

15 — их различия проявляются только в членах, пропорциональных 𝑇 3 [3]. На высоких температурах 𝐶𝑉 → 3/2 и 𝐶𝑃 → 5/2. 101 CV , CP , atomic units CV = βv 2/3 T CP = 5/2 CV = 3/2 CV : CV : CV : CP : CP : CP : 100 v = 0.1 v=1 v = 10 v = 0.1 v=1 v = 10 CV : low-T asymp. CV : high-T asymp. CP : high-T asymp. 10−1 10−2 10−1 100 101 102 103 T , Ha Рис. 1.2. Теплоемкости при постоянном объеме и давлении ИФГ на трех изохорах. at_capacity} Отношение 𝐶𝑃 /𝐶𝑉 и адиабатическую скорость звука можно найти на графике 1.3 в виде тех же трех изохор 𝑣 = 0.1, 1 и 10. Заметим, что 𝐶𝑃 /𝐶𝑉 не является константой для ИФГ и стремится к зависящему от 𝑣 значению при 𝑇 → 0. Адиабатическая скорость звука проявляет схожее проведение; при высоких 𝑇 она стремится к независящему от объема значению 5/3𝑇 . 101 2.00 CS = 1.75 CP /CV = 5/3 CS , atomic units 1.50 CP /CV 1.25 nd_velocity} 1.00 0.75 0.50 0.25 0.00 10−2 10−1 100 101 T , Ha (а) 102 103 104 10−2 (5/3)T v = 0.1 v=1 v = 10 low-T asymp. high-T asymp. 100 v = 0.1 v=1 v = 10 High-T asymp. p 10−1 100 101 102 T , Ha (б) Рис. 1.3. Отношение теплоемкостей (a) и адиабатическая скорость звука (б) ИФГ на трех изохорах.

16 Глава 2 Модель взаимодействующего электронного газа на однородном компенсирующем фоне 2.1. Классическая модель однокомпонентной плазмы Классическая модель однокомпонентной плазмы (ОКП) — это модель точечных ионов, помещенных для обеспечения электронейтральности в рав номерно распределенную среду заряда противоположного знака. Такая мо дель является хорошим приближением для плазмы сверхвысоких давлений, реализующихся, например, в центре белых карликов и тяжелых планет ти па Юпитера. В этих случаях под действием давления вещество полностью ионизовано, а вырожденные электроны обладают достаточной кинетической энергией 𝜀𝑘 ≈ 𝜀𝐹 , чтобы образовать почти однородное фоновое распределе ние зарядовой плотности. Модель ОКП является простейшей нетривиальной моделью плазмы, так как вид потенциала взаимодействия здесь не вызывает сомнений, а отсутствие квантовых эффектов позволяет исключить из рас смотрения образование связанных состояний (молекул, атомов, ионов) и вли яние вырождения и интерференции. Важным свойством модели ОКП являет ся зависимость ее свойств только от одного параметра 𝛾 = 𝑞𝑖2 /(𝑘𝑇 𝑟¯) в силу однородности кулоновского потенциала, где 𝑞𝑖 — заряд иона, 𝑟¯ — среднее рас стояние между ионами, обычно определяемое соотношением 4𝜋¯ 𝑟3 /3 = 1/𝑛𝑖 , 𝑛𝑖 — концентрация ионов. Более подробно с моделью и ее основными резуль татами можно ознакомиться в [9]. Большое количество результатов по ОКП получено методом Монте– Карло. Основным параметром в этих результатах является бинарная корре ляционная функция 𝑔(𝑟). На ее основании возможно вычислить внутреннюю

17 энергию: al_energy_g} 𝑢 = 3𝑛i 𝑘𝑇 /2 + 𝑢кор , (2.1) eq:u_corr_g} ∫︁ (2.2) corr_approx} eq:p_approx} -gas-approx} 𝑢кор / (𝑛i 𝑘𝑇 ) = (𝑛i /2𝑘𝑇 ) (︀ )︀ 𝑑r 𝑍 2 𝑒2 /𝑟 [𝑔(𝑟) − 1]. Результаты вычислений методом Монте-Карло для 1 ≤ 𝛾 ≤ 160 были аппроксимированы в [10] с погрешностью 3 · 10−5 следующим выражением: 𝑢кор / (𝑛i 𝑘𝑇 ) = 𝑎𝛾 + 𝑏𝛾 1/4 + 𝑐𝛾 −1/4 + 𝑑 (2.3) где 𝑎 = −0,89752, 𝑏 = 0,94544, 𝑐 = 0,17954, 𝑑 = −0,80049. Давление можно получить по теореме вириала 𝑝𝑉 = 𝑢/3. )︂ 𝑛i 𝑘𝑇 (︂ 3 𝑝= = + 𝑎𝛾 + 𝑏𝛾 1/4 + 𝑐𝛾 −1/4 + 𝑑 3𝑉 3𝑉 2 𝑢 (2.4) Интегрируя (2.3) по 𝛾, можно получить свободную энергию ОКП: 𝑓 = 𝐹/ (𝑛i 𝑘𝑇 ) = (︁ )︁ 1/4 −1/4 = 𝑎𝛾 + 4 𝑏𝛾 − 𝑐𝛾 + (𝑑 + 3) ln 𝛾 − (𝑎 + 4𝑏 − 4𝑐 + 1,135) (2.5) Заметим, что в соответствии с (2.4) высокое значение 𝛾 приводит к отрицательному давлению компонента (поскольку 𝑎 < 0). Это не приводит к неустойчивости ОКП: полное давление оказывается положительным за счет большого давления электронного газа. Для ОКП возможно явление кристаллизации, впервые рассмотренное Вигнером в [11]. Было показано, что классический электронный газ с концен трацией 𝑛𝑒 на фоне компенсирующего заряда должен при достаточно низких концентрациях образовывать упорядоченную структуру. Это вызвано тем, 1/3 что стабилизирующая решетку кулоновская энергия 𝑉C ∼ 𝑒2 𝑛e ∼ 𝑒2 /¯ 𝑟 при расширении плазмы уменьшается медленнее, чем разрушающая решетку ки 2/3 нетическая энергия 𝜀к ∼ 𝜀𝐹 ∼ ~2 𝑛e /2𝑚. Поэтому при малых плотностях 2/3 1/3 кинетическая энергия 𝜀к ∼ 𝑛e становится меньше потенциальной 𝑉C ∼ 𝑛e

18 и не способна разрушить упорядоченную структуру электронов, возникшую из-за отталкивания. В работе [10] были получены результаты вычислений для кристалли ческой ОКП с объемно-центрированной решеткой. Избыточная энергия при 160 . 𝛾 . 300 в [10] аппроксимирована выражением: 𝑢кор /𝑛i 𝑘𝑇 = 𝑎BCC 𝛾 + olid_approx} 3 −2 𝑏𝛾 = 0,895929𝛾 + 1,5 + 2980/𝛾 2 2 (2.6) Из выражения (2.6) следует выражение для плотности свободной энергии решетки: 𝑓 (𝛾) = −0,895929 + 9𝛾/2 − 1,8856 − 1490/𝛾 2 . olid-approx} (2.7) Зависимости свободных энергий газовой (2.5) и твердой (2.7) фаз пересе каются при 𝛾𝑚 = 165 [10]. При этом значении параметра неидеальности происходит вигнеровская кристаллизация ОКП. 2.2. Суммирование Эвальда {sec:ewald} tention_def} Для расчета энергии в системах с дальнодействующими потенциалами широко используется метод суммирования по Эвальду [12]. В данном методе потенциал взаимодействия разбивается на два слагаемых: короткодействую щее и дальнодействующее. Первое слагаемое вычисляется в действительном пространстве, в то время как для расчета второго используется преобразова ние Фурье. По сравнению с прямым подсчетом данный метод дает быструю сходимость по энергии, что обеспечивает высокую точность и скорость рас четов. Более подробно, в методе суммирования по Эвальду предлагается пере писать потенциал взаимодействия 𝜙(r) в виде: def 𝜙(r) = 𝜙𝑠𝑟 (r) + 𝜙𝑙𝑟 (r), (2.8) где 𝜙𝑠𝑟 (r) отвечает за короткодействующую часть, энергия которой быстро

d:energy_lr} 19 сходится, а 𝜙𝑙𝑟 (r) соответствует дальнодействующей части, энергия которой суммируется в Фурье-пространстве. Предполагается, что энергия короткодействующей части потенциала суммируется непосредственно и основная проблема состоит в суммировании дальнодействующего вклада в потенциал. Также из-за использования Фурье– преобразования метод неявно подразумевает, что рассматриваемая система является бесконечно периодической, что является разумным предположени ем для кристаллических структур. Будем называть примитивной ячейкой минимальную повторяющуюся часть периодической системы. Также для определенности выберем одну ячей ку и назовем ее центральной, а остальные будем называть образами. Энергия дальнодействующего взаимодействия есть сумма энергий вза имодействия между зарядами центральной ячейки и всех остальных зарядов решетки. Ее можно записать в виде: ∫︁ ∫︁ 𝐸𝑙𝑟 = 𝑑r𝑑r′ 𝜌𝑡𝑜𝑡 (r)𝜌𝑢𝑐 (r′ ) 𝜙𝑙𝑟 (r − r′ ) , где введена плотность заряда примитивной ячейки 𝜌𝑢𝑐 (r), представляющая из себя сумму по зарядам внутри одной ячейки: ∑︁ def :rho_uc_def} 𝜌𝑢𝑐 (r) = 𝑞𝑘 𝛿 (r − r𝑘 ) (2.10) charges 𝑘 и суммарная плотность заряда 𝜌𝑡𝑜𝑡 (r), представляющая из себя ту же сумму, но с добавлением зарядов от образов: ∑︁ ∑︁ def 𝜌𝑡𝑜𝑡 (r) = rho_tot_def} 𝑞𝑘 𝛿 (r − r𝑘 − 𝑛1 a1 − 𝑛2 a2 − 𝑛3 a3 ) . unction_def} (2.9) (2.11) 𝑛1 ,𝑛2 ,𝑛3 charges 𝑘 Здесь 𝛿(x) — дельта-функция Дирака, a1 , a2 , a3 — вектора решетки и 𝑛1 , 𝑛2 , 𝑛3 пробегают по всем целым числам. На (2.11) можно смотреть как на свертку (2.10) с функцией ячейки 𝐿(r): def 𝐿(r) = ∑︁ 𝑛1 ,𝑛2 ,𝑛3 𝛿 (r − 𝑛1 a1 − 𝑛2 a2 − 𝑛3 a3 ) . (2.12)

20 Для полученной свертки имеем: ˜ 𝜌𝑢𝑐 (k), 𝜌˜𝑡𝑜𝑡 (k) = 𝐿(k)˜ t_conv_prod} где тильдами обозначены Фурье-образы. Для функции ячейки 𝐿(r) Фурье-образ есть: (2𝜋)3 ∑︁ ˜ 𝐿(k) = 𝛿 (k − 𝑚1 b1 − 𝑚2 b2 − 𝑚3 b3 ) , Ω 𝑚 ,𝑚 ,𝑚 ction_image} ewald:v_def} 1 где b1 = 2 (2.14) 3 a2 × a3 , b2 и b3 получаются циклическими перестановками, ин Ω дексы 𝑚1 , 𝑚2 , 𝑚3 по-прежнему пробегают все целые числа и Ω есть объем ˜ центральной ячейки: Ω = a1 · (a2 × a3 ). Стоит заметить, что и 𝐿(r), и 𝐿(r) являются действительными четными функциями. Вводя эффективный одночастичный потенциал ∫︁ def 𝑣(r) = 𝑑r′ 𝜌𝑢𝑐 (r′ ) 𝜙𝑙𝑟 (r − r′ ) , (2.15) можно переписать (2.9) в виде: 𝐸𝑙𝑟 = y_v_rewrite} ∫︁ 𝑑r𝜌𝑡𝑜𝑡 (r)𝑣(r) (2.16) Заметим, что (2.15) тоже является сверткой, что позволяет записать: def 𝑉˜ (k) = 𝜌˜𝑢𝑐 (k)Φ̃(k), r_conv_prod} d:v_fourier} (2.13) (2.17) где, как и в (2.13), тильдами обозначены Фурье-образы и введен Фурье-образ (2.15): 𝑉˜ (k) = ∫︁ 𝑑r𝑣(r)𝑒−𝑖k·r . (2.18) Согласно теореме Планшереля, суммирование для получения 𝐸𝑙𝑟 можно провести в Фурье-пространстве, что дает: ∫︁ ∫︁ 𝑑k * 𝑑k ˜ * ˜ 𝐸𝑙𝑟 = 𝜌 ˜ (k) 𝑉 (k) = 𝐿 (k) |˜ 𝜌𝑢𝑐 (k)|2 Φ̃(k) = 𝑡𝑜𝑡 3 3 (2𝜋) (2𝜋) 1 ∑︁ = |˜ 𝜌𝑢𝑐 (k)|2 Φ̃(k), Ω 𝑚 ,𝑚 ,𝑚 1 2 3

21 1 ∑︁ 𝐸𝑙𝑟 = |˜ 𝜌𝑢𝑐 (k)|2 Φ̃(k), Ω 𝑚 ,𝑚 ,𝑚 nergy_final} 1 2 (2.19) 3 где в суммировании k = 𝑚1 b1 + 𝑚2 b2 + 𝑚3 b3 . Формула (2.19) есть основной результат. Фурье-образ 𝜌˜𝑢𝑐 (k) можно под считать, к примеру, алгоритмами БПФ (англ. FFT). Если найден фурье-образ 𝜌˜𝑢𝑐 (k), то суммирование (или интегрирование) по k становится простым и приобретает быструю сходимость. Стоит отметить, что используемая прими тивная ячейка должна быть электронейтральной во избежание бесконечных сумм. На практике суммирование (2.19) ведется в конечном диапазоне зна чений 𝑚1 , 𝑚2 , 𝑚3 . Благодаря особенностям метода отбрасывание остальных членов приводит к небольшой потере точности, однако значительно уменьша ет вычислительное время. Оценку возникающих из-за отбрасывания ошибок можно найти, например, в [13]. 2.3. Гамильтониан модели ВЭГ Взаимодействующий электронный газ (ВЭГ) — это физическая модель электронного газа, в которой учитывается электростатическое взаимодей ствие между всеми зарядами в системе. В данной работе рассматривается важный частный случай этой модели — так называемая модель «желе», в которой электронейтральность обеспечивается компенсирующим фоном по ложительного заряда с постоянной полностью (так, чтобы суммарный заряд системы был равен нулю). Фон предполагается однородным и несжимаемым. При изложении данной главы будем использовать систему единиц СГС. Для модели «желе» из 𝑁 электронов, заключенных в объеме 𝑉 , имеющих ∑︀𝑁 концентрацию 𝜌(r) = 𝑖=1 𝛿(r − ri ) при концентрации фоновых зарядов

22 𝑛(R) = 𝑁/𝑉 , гамильтониан имеет вид [14]: ˆ =𝐻 ˆ 𝑒𝑙 + 𝐻 ˆ 𝑏𝑎𝑐𝑘 + 𝐻 ˆ 𝑒𝑙−𝑏𝑎𝑐𝑘 . 𝐻 am_four_sum} (2.20) ˆ 𝑒𝑙 — это электронный гамильтониан, состоящий из кинетической Здесь 𝐻 энергии и электрон-электронного отталкивания: ˆ 𝑒𝑙 = 𝐻 :ham_el_def} 𝑁 𝑁 ∑︁ ∑︁ p̂2𝑖 𝑒2 + . 2𝑚 |r̂ − r̂ | 𝑖 𝑗 𝑖=1 𝑖<𝑗 (2.21) ˆ 𝑏𝑎𝑐𝑘 есть гамильтониан положительного фонового заряда, описывающий его 𝐻 ˆ 𝑏𝑎𝑐𝑘 есть: взаимодействие с самим собой. Усредненное значение 𝐻 ∫︁ ⟨ ⟩ 𝑒2 ∫︁ ′ 𝑁2 ′ 𝑛(R)𝑛 (R ) ˆ 𝐻𝑏𝑎𝑐𝑘 = 𝑑R 𝑑R back_expect} = lim 𝑣𝑞 , 2 |R − R′ | 2𝑉 q→0 𝑉 (2.22) 𝑉 где введена фурье-компонента кулоновского потенциала 𝑣𝑞 = 4𝜋𝑒2 𝑞2 . Вклад электростатического взаимодействия электрона с фоном описывается членом 𝐻𝑒𝑙−𝑏𝑎𝑐𝑘 : ⟨ ⟩ ˆ 𝑒𝑙−𝑏𝑎𝑐𝑘 = − 𝐻 ∫︁ 𝑉 back_expect} 𝑑r ∫︁ 𝑉 𝑒2 𝜌(r)𝑛(R) 𝑑R = |r − R| ⟩ ⟨ 𝑁2 ˆ =− lim 𝑣𝑞 = −2 𝐻𝑏𝑎𝑐𝑘 . (2.23) 𝑉 q→0 Для конечных систем результаты (2.22) и (2.23) получаются конечными, поскольку 𝑞 достигает отличное от нуля минимальное значение, зависящее от объема системы. Для систем с периодическими граничными условиями суммирование кулоновского взаимодействия можно провести упомянутым выше методом Эвальда (см. раздел 2.2). Учтем также возможное размытие заряда относи тельно его равновесного положения в форме функции Гаусса [15]. Введем потенциал в точке r, создаваемый всей решеткой без одного

23 лишь заряда 𝑗, но с учетом его образов: 𝑁 1 ∑︁ ∑︁ ′ 𝑞𝑖 𝜙−𝑗 (r) = , 𝜀 n 𝑖=1 |r − r𝑖 + n𝐿| el:phi_no_j} где штрих у суммы обозначает, что 𝑖 ̸= 𝑗 при n = 0. Здесь 𝑞𝑖 указан для того, чтобы получающиеся результаты не были привязаны к типам и значениям зарядов. Используя выражения для полной энергии: 𝑁 oumb_energy} и плотности заряда: (2.25) 𝜌𝑖 (r) = 𝑞𝑖 · 𝛿 (r − ri ) , (2.26) 𝑁 𝑁 1 ∑︁ ∑︁ ∑︁ ′ 𝐸= 𝜀 n 𝑖=1 𝑗=1 ∫︁ ∫︁ 𝜌𝑖 (r) · 𝜌𝑗 (r′ ) 3 3 ′ 𝑑 r𝑑 r . |r − r′ + n𝐿| (2.27) Разбиваем заряды на части: 𝜌𝑖 (r) = 𝜌𝑆𝑖 (r) + 𝜌𝐿𝑖 (r), 𝜌𝑆𝑖 (r) = 𝑞𝑖 𝛿(r − ri ) − 𝑞𝑖 𝐺𝜎 (r − ri ), _separation} (2.28) 𝜌𝐿𝑖 (r) = 𝑞𝑖 𝐺𝜎 (r − ri ), где введено гауссово размытие заряда 𝐺𝜎 (r) = :gauss_blur} phi_poisson} 1 ∑︁ 𝑞𝑗 𝜙−𝑗 (r𝑗 ) 𝐸= 2 𝑗=1 получим полную энергию в виде: loumb_total} _separation} (2.24) ⎡ |r|2 ⎤ 1 exp ⎣− 2 ⎦ . 2𝜎 (2𝜋𝜎 2 )3/2 (2.29) Дисперсия здесь будет играть роль некоторого параметра. Выражение для энергии (2.25) с учетом (2.28) перепишется в виде: 𝑁 𝑁 𝑁 1 ∑︁ 1 ∑︁ 1 ∑︁ 𝑆 𝐿 𝑞𝑖 𝜙−𝑖 (ri )+ 𝑞𝑖 𝜙 (ri )− 𝑞𝑖 𝜙𝐿𝑖 (ri ) = 𝐸 𝑆 +𝐸 𝐿 −𝐸 𝑠𝑒𝑙𝑓 . (2.30) 𝐸= 2 𝑖=1 2 𝑖=1 2 𝑖=1 Уравнение Пуассона для случая, когда заряд распределен по функции Гаусса, может быть решено аналитически. Само уравнение имеет вид: ∆𝜙𝑖 (r) = − 4𝜋 𝐺𝜎 (r), 𝜀 (2.31)

24 а его решение выражается через функцию ошибок erf(𝑥) = дующим образом: 1 erf 𝜙𝜎 (𝑟) = 𝜀𝑟 hi_solution} (︂ 𝑟 √ 2𝜎 )︂ √2 𝜋 ∫︀ 𝑥 0 . 2 𝑒−𝑡 𝑑𝑡 сле (2.32) Таким образом, ]︂ [︂ 1 𝑞𝑖 |r − ri | = , erfc √ 𝜀 |r − ri | 2𝜎 [︂ ]︂ 1 𝑞 |r − r | 𝑖 i 𝜙𝐿𝑖 (r) = erf √ , 𝜀 |r − ri | 2𝜎 𝜙𝑆𝑖 (r) _S_solution} (2.33) где erfc(𝑥) = 1 − erf(𝑥). Теперь, с учетом размытия, можно переписать 𝜙𝑆−𝑖 (r) в виде: )︂ (︂ 𝑁 ∑︁ ∑︁ 1 𝑞 |r − r + n𝐿 𝑗 ′ √j , erfc 𝜙𝑆−𝑖 (r) = 𝜀 n 𝑗=1 |r − rj + n𝐿| 2𝜎 а энергию: (︂ )︂ 𝑁 𝑁 𝑁 ∑︁ ∑︁ ∑︁ ∑︁ 1 𝑞 𝑞 1 |r − r + n𝐿| 𝑖 𝑗 i j ′ √ 𝐸𝑆 = 𝑞𝑖 𝜙𝑆−𝑖 (ri ) = erfc . 2 𝑖=1 2𝜀 n 𝑖=1 𝑗=1 |ri − rj + n𝐿| 2𝜎 Это выражение аналогично исходному лишь с той разницей, что сумма обре зается функцией erfc, что упрощает расчет. Используя erfc(𝑥) −−→ 𝑥→0 √2 𝑥, 𝜋 можно получить для 𝐸 𝑠𝑒𝑙𝑓 : 𝑁 𝐸 _self_final} 𝑠𝑒𝑙𝑓 1 1 ∑︁ 2 = √ 𝑞 . 𝜀 2𝜎 𝑖=1 𝑖 (2.34) Таким образом, осталось посчитать третью, дальнодействующую, часть. Как было отмечено в разделе 2.2, дальнодействующая часть считается в об ратном пространстве. Записывая плотность заряда в виде суммы по всем зарядам системы, периодически повторяющимся в пространстве: 𝐿 𝜌 (r) = 𝑁 ∑︁ ∑︁ n 𝑗=1 𝑞𝑗 𝐺𝜎 (r − rj + n𝐿),

25 имеем для Фурье-образа: 𝐿 𝜌˜ (k) = ∫︁ ∑︁ 𝑁 𝑗=1 𝑉 𝑞𝑗 𝐺𝜎 (r − rj + n𝐿)𝑒 l:E_L_final} 𝑞𝑗 𝑒−𝑖krj 𝑒−𝜎 2 2 𝑘 /2 𝑑r. 𝑗=1 exp(−𝑖kn𝐿) = 1. Тогда для Фурье-образа потенциала имеем: 𝑁 4𝜋 ∑︁ −𝑖krj 𝑒−𝜎 𝑘 𝜙˜ (k) = 𝑞𝑗 𝑒 𝜀 𝑗=1 𝑘2 2 2 𝐿 /2 , (2.35) а сам потенциал (через обратное преобразование): 𝑁 4𝜋 ∑︁ ∑︁ 𝑞𝑗 𝑖k(r−rj ) −𝜎2 𝑘2 /2 𝜙 (r) = 𝑒 𝑒 . 2 𝑉𝜀 𝑘 𝑗=1 𝐿 (2.36) k̸=0 Вклад члена с k = 0 нулевой, так как полный заряд примитивной ячейки предполагается нулевым. Для энергии 𝐸 𝐿 получаем: 𝑁 𝑁 4𝜋 ∑︁ ∑︁ ∑︁ 𝑞𝑖 𝑞𝑗 𝑖k(ri −rj ) −𝜎2 𝑘2 /2 𝑒 𝑒 . 𝐸 = 2 2𝑉 𝜀 𝑘 𝑖=1 𝑗=1 𝐿 (2.37) k̸=0 Определяя структурный фактор 𝑆(k) = ∑︀𝑁 𝑖kri , 𝑖=1 𝑞𝑖 𝑒 дальнодействую щую часть энергии возможно переписать в виде: 4𝜋 ∑︁ 𝑒−𝜎 𝑘 𝐿 𝐸 = 2𝑉 𝜀 𝑘2 2 2 :E_L_from_S} nergy_final} 𝑑r = 𝑁 ∑︁ Здесь использован тот факт, что 𝑘 есть вектор обратной решетки и потому i_L_fourier} i_L_spatial} −𝑖kr k̸=0 /2 |𝑆(k)|2 (2.38) и получить окончательное выражение для полной энергии: 𝐸 = 𝐸 𝑆 + 𝐸 𝐿 − 𝐸 𝑠𝑒𝑙𝑓 = (︂ )︂ 𝑁 𝑁 1 ∑︁ ∑︁ ∑︁ ′ 𝑞𝑖 𝑞𝑗 |ri − rj + n𝐿| √ = erfc + 2𝜀 n 𝑖=1 𝑗=1 |ri − rj + n𝐿| 2𝜎 4𝜋 ∑︁ 𝑒−𝜎 𝑘 + 2𝑉 𝜀 𝑘2 2 2 k̸=0 /2 𝑁 1 1 ∑︁ 2 |𝑆(k)| − √ 𝑞 . 𝜀 2𝜋𝜎 𝑖=1 𝑖 2 (2.39)

26 Собирая результаты воедино, получим гамильтониан модели ВЭГ: ˆ = 𝐻 )︂ (︂ 𝑁 𝑁 𝑞𝑖 𝑞 𝑗 1 ∑︁ ∑︁ ∑︁ ′ |r̂𝑖 − r̂𝑗 + n̂𝐿| √ + + erfc 2𝑚 2𝜀 n 𝑖=1 𝑗=1 |r̂𝑖 − r̂𝑗 + n̂𝐿| 2𝜎 𝑁 p̂2 ∑︁ 𝑖 𝑖=1 4𝜋 ∑︁ 𝑒−𝜎 𝑘 + 2𝑉 𝜀 𝑘2 2 2 l:ham_final} k̸=0 /2 𝑁 1 1 ∑︁ 2 2 ˆ |𝑆(k)| − √ 𝑞 . (2.40) 𝜀 2𝜋𝜎 𝑖=1 𝑖 𝑁 ∑︀ ˆ Здесь 𝑆(k) = 𝑞𝑖 𝑒𝑖k^r𝑖 , а оператор n̂𝐿 соответствует координатному сдвигу 𝑖=1 на целое число ячеек. 2.4. Усреднение по направлениям и изотропная модель ВЭГ Непосредственное суммирование по Эвальду в том виде, который при веден в секции 2.2, имеет два практических недостатка: 1. Вычислительная сложность (нагрузка на процессор компьютера) быстро растет с ростом числа частиц в главной ячейке (как 𝑁 2 ). 2. Комбинирование дальнодействующего кулоновского взаимодействия с периодическими граничными условиями может привести к появлению неизотропического электрического поля, имеющего кубическую сим метрию в кристаллической решетке, составленной из главных ячеек. Поэтому любая процедура суммирования кулоновских сил в неупоря доченной системе может считаться подходящей лишь в том случае, ес ли погрешность результата, вызванная этой искуственной симметрией, пренебрежимо мала. В [16] предлагается новый метод суммирования, основанный на предвари тельном усреднении вдоль всех направлений главной ячейки. Предполагает ся, что ячейка имеет вид куба с объемом 𝑉 = 𝐿3 и на систему наложены периодические граничные условия вместе с условиями электронейтрально сти.

27 Для начала заметим, что подбором параметра размытия 𝜎 можно до биться возможности отбросить значительную часть членов суммирования в (2.39) [17]. В дальнейшем будем предполагать, что используемый параметр размытия позволяет осуществить такое упрощение. Приведем (2.39) к виду, полученному в оригинальной работе Эвальда [12]. Для этого необходимо сделать замену 𝛿 = 2𝜋 2 𝜎 , а |𝑆(k)| явно раскрыть с (︀ )︀ выделением тригонометрических функций. Члены порядка exp −𝛿 2 𝐿2 будут давать малый вклад в энергию, если выбрать типичное значение 𝛿 = 5/𝐿 [17]. Рассматривая энергию кулоновского взаимодействия 𝑁 ионов в главной ячейке, получим: 𝐸𝑁 = b_main-cell} :phi_as_sum} √ 𝑁 ∑︁ 𝑞𝑖 𝜙(r𝑖 ), (2.41) 𝑖=1 где 𝜙(r𝑖 ) — электростатический потенциал 𝑖-го иона в точке r𝑖 . Согласно [12], 𝜙 имеет вид: 𝑁 1 ∑︁ 𝜙(r𝑖 ) = 𝜙1 (r𝑖 ) + 𝜙2 (r𝑖 , r𝑗 ), 2 (2.42) 𝑗̸=𝑖 где в отстутствии внешнего поля 𝜙1 есть константа: ⎛ ⎞ 𝑞𝑖 𝛿 1 ∑︁ 1 −𝜋2 𝑛2 /𝛿2 ⎠, √ 𝜙1 = ⎝ 𝑒 − :mean:phi_1} 𝐿 2𝜋 𝑛>0 𝑛2 𝜋 (2.43) а парный потенциал получается следующим: ⎡ ⎛ ⎞ ⎤ (︂ )︂ 𝑟𝑖𝑗 1 1 ∑︁ 1 −𝜋2 𝑛2 /𝛿2 2𝜋 :mean:phi_2} 𝜙2 = 𝑞𝑗 ⎣ erfc ⎝𝛿 ⎠ + 𝑒 cos n · r𝑖𝑗 ⎦ . (2.44) 2 𝑟𝑖𝑗 𝐿 2𝜋𝐿 n>0 𝑛 𝐿 Здесь n/𝐿 есть трехмерный взаимный вектор решетки (𝑛 = |n|), r𝑖𝑗 = r𝑖 − r𝑗 , 𝑟𝑖𝑗 = |r𝑖𝑗 |. Параметр 𝛿/𝐿 обычно называют параметром Эвальда. Учитывая, что все направления главной решетки в изотропической жид кости должны быть равноправны, обе части (2.44) можно усреднить вдоль

28 всех всех направлений вектора n. Обозначая 𝜙2 (𝑟𝑖𝑗 ) ≡ ⟨𝜙2 (r𝑖𝑗 )⟩, получим: ⎛ ⎞ ⎡ ⎤ (︂ )︂ 𝑟𝑖𝑗 𝑞𝑗 ⎣ 1 ∑︁ 1 −𝜋2 𝑛2 /𝛿2 2𝜋 ⎝ ⎠ 𝜙2 (𝑟𝑖𝑗 ) = + 2 _2_averaged} erfc 𝛿 𝑒 sin 𝑛𝑟𝑖𝑗 ⎦ . (2.45) 𝑟𝑖𝑗 𝐿 2𝜋 n>0 𝑛3 𝐿 Поскольку erfc(𝑥) − 1 и sin(𝑥) являются нечетными функциями, возможно следующее разложение в ряд: 𝜙2 (𝑟𝑖𝑗 ) = hi_2_series} an:C_coeffs} 𝑞𝑗 𝑟𝑖𝑗 (︃ 1+ ∑︁ )︃ 𝐶𝑘 𝑟𝑖𝑗2𝑘+1 , 𝑘≥0 (2.46) где коэффициенты 𝐶𝑘 можно найти прямым разложением в ряд Маклорена (2.45) с использованием формулы Эйлера-Маклорена, обобщенной на случай суммирования по трехмерному пространству. Это в результате даст следую щие выражения для членов ряда [16]: 1 ∑︁ 1 −𝜋2 𝑛2 /𝛿2 2𝛿 √ 𝐶0 = , 𝑒 − 𝜋 n>0 𝑛2 𝜋 𝐶1 = 2𝜋 , 3𝐿3 𝐶𝑘 = 0, (2.47) 𝑘 > 1. Если теперь учесть условие электронейтральности, то окажется, что член в (2.46), не зависящий от расстояния (пропорциональный 𝐶0 ), уничтожает вклад от 𝜙1 в (2.42). Это означает, что суммарная энергия кулоновского взаи ∑︀ модействия в главной ячейке может быть описана суммой 𝐸𝑁 = 𝑁 𝑖=1 𝜙(𝑟𝑖𝑗 ), где 𝑞𝑖 𝑞𝑗 𝜙(𝑟𝑖𝑗 ) = 𝑟𝑖𝑗 n:phi_final} [︃ 1 1+ 2 (︂ 𝑟𝑖𝑗 𝑟𝑚 )︂3 ]︃ (2.48) и 𝑟𝑚 = (3/4𝜋)1/3 𝐿 — радиус эквивалентной по объему сферы главной ячейки: 4 3 3 𝜋𝑟𝑚 = 𝐿3 . Результат (2.48) и есть усредненный по направлениям потенциал. Он обладает следующими свойствами: 1. Он стремится к чистому кулоновскому парному потенциалу при малых межионных расстояниях;

29 2. Его первая производная равна нулю в 𝑟𝑚 ; 3. Его значение в минимуме 𝑟 = 𝑟𝑚 отлично от нуля и равно 𝜙(𝑟𝑚 ) = 3𝑞𝑖 𝑞𝑗 /2𝑟𝑚 . Значение потенциала в точке 𝑟𝑚 можно обнулить, прибавив к нему посто янную добавку −𝜙(𝑟𝑚 ). Воспользовавшись условием электронейтральности, получаем окончательное выражение суммарной кулоновской энергии главной ячейки в следующем виде: 𝑁 𝑁 𝑁 ∑︁ 3𝑞𝑖2 1 ∑︁ ∑︁ 𝐸𝑁 = − 𝜙(𝑟 ˜ 𝑖𝑗 ), + 4𝜋𝑟 2 𝑚 𝑖=1 𝑖=1 -cell_final} n:phi_tilde} n:ham_final} (2.49) 𝑗=1,𝑗̸=𝑖 где ⎧ ⎪ ⎪ 𝑞 𝑖 𝑞𝑗 ⎪ ⎪ ⎨ 𝑟 𝜙(𝑟) ˜ = ⎪ ⎪ ⎪ ⎪ ⎩0, ⎧ ⎪ ⎨ ⎤⎫ ⎞ ⎡⎛ ⎞2 ⎪ ⎬ 1 𝑟 𝑟 ⎢ ⎥ 1 + ⎝ ⎠ ⎣⎝ ⎠ − 3⎦ , ⎪ ⎪ 2 𝑟𝑚 𝑟𝑚 ⎩ ⎭ ⎛ 𝑟 < 𝑟𝑚 (2.50) 𝑟 ≥ 𝑟𝑚 есть эффективный короткодействующий потенциал межчастичного взаимо действия, который равен нулю вместе со своими производными в точке 𝑟𝑚 и остается таковым при 𝑟 > 𝑟𝑚 . Объединяя полученное выражение с кинетической энергией, получим гамильтониан изотропной модели ВЭГ в координатном представлении: ˆ = 𝐻 𝑁 𝑁 𝑁 𝑁 ∑︁ ∑︁ p̂2𝑖 1 ∑︁ ∑︁ 3𝑞𝑖2 + − 𝜙(ˆ ˜ 𝑟𝑖𝑗 ). 2𝑚 4𝜋ˆ 𝑟 2 𝑚 𝑖=1 𝑖=1 𝑖=1 𝑗=1,𝑗̸=𝑖 (2.51)

30 Глава 3 Программная реализация метода функционала плотности для модели ВЭГ 3.1. Метод функционала электронной плотности Метод функционала плотности (МФП) представляет собой точную квантовомеханическую теорию для системы взаимодействующих квантовых частиц во внешнем потенциале 𝑉𝑒𝑥𝑡 (r). Сам метод основан на двух строго доказанных теоремах [18]: 1. Для любой системы взаимодействующих частиц во внешнем потенциа ле 𝑉𝑒𝑥𝑡 (r), потенциал 𝑉𝑒𝑥𝑡 (r) с точностью до произвольной постоянной определяется электронной плотностью основного состояния 𝑛0 (r). 2. Энергия невырожденного основного состояния системы для любого внешнего потенциала 𝑉𝑒𝑥𝑡 (r) является функционалом электронной плот ности 𝐸[𝑛(r)]. Основное состояние системы является минимумом дан ного функционала, который достигается при плотности, соответствую щей основному состоянию системы 𝑛0 (r). Из первой теоремы следует, что гамильтониан системы с точностью до аддитивной константы определяется электронной плотностью основного со стояния 𝑛0 (r). Следовательно, определены многочастичные волновые функ ции для всех состояний (основного и возбужденных). Таким образом, все свойства системы полностью определены, если известна электронная плот ность основного состояния 𝑛0 (r). Из второй теоремы следует, что при извест ном функционале 𝐸[𝑛] можно определить плотность и энергию основного состояния.

31 Функционал энергии в формулировке [18] можно записать следующим образом: 𝐸𝐻𝐾 [𝑛] = 𝑇 [𝑛] + 𝐸𝑖𝑛𝑡 [𝑛] + ∫︁ 𝑑3 𝑟𝑉𝑒𝑥𝑡 (r)𝑛(r) + 𝐸𝐼𝐼 ≡ 𝐹𝐻𝐾 [𝑛] + ft:E_HK-def} t:Omega-min} t:rho_0-bka} 𝑑3 𝑟𝑉𝑒𝑥𝑡 (r)𝑛(r) + 𝐸𝐼𝐼 , (3.1) где 𝐸𝐼𝐼 — энергия взаимодействия ядер. Функционал 𝐹𝐻𝐾 [𝑛] включает в себя кинетическую и потенциальную энергию системы взаимодействующих электронов, 𝐹𝐻𝐾 [𝑛] = 𝑇 [𝑛] + 𝐸𝑖𝑛𝑡 [𝑛]. ft:F_HK-def} t:Omega_rho} ∫︁ (3.2) Минимум функционала (3.1) опрелеляет основное состояние энергии систе мы и ее электронную плотность. Теория функционала плотности допускает обобщение на случай конеч ных температур [19]. Функционал энергии от электронной плотности заме няется на функционал большого термодинамического потенциала Ω от опе ратора плотности 𝜌ˆ: [︂ ]︂ 1 ˆ − 𝜇𝑁 ˆ ) + ln 𝜌ˆ . Ω[ˆ 𝜌] = Sp 𝜌ˆ(𝐻 𝛽 (3.3) Минимум этого функционала совпадает с термодинамически равновесным выражением для большого термодинамического потенциала: ^ ^ Ω = Ω[𝜌ˆ0 ] = − ln Sp 𝑒−𝛽(𝐻−𝜇𝑁 ) , (3.4) где 𝜌ˆ0 — оператор плотности в большом каноническом ансамбле (БКА), ^ 𝜌ˆ0 = ^ 𝑒−𝛽(𝐻−𝜇𝑁 ) ^ ^) 𝑁 Sp 𝑒−𝛽(𝐻−𝜇 . (3.5) Теормера Мермина [19] утверджает, что не только энергия, но также и все термодинамические функции (энтропия, теплоемкость) являются функ ционалами равновесной плотности.

ft:E_KS-def} ensity-kohn} dft:T_s-def} 32 Одна из основных проблем МФП — это отсутствие в общем случае пря мой связи между кинетической энергией и функцией электронной плотности. Для решения этой проблемы в [20] высказано предположение, что основное состояние взаимодействующей системы частиц совпадает с основным состо янием эквивалентной системы невзаимодействующих частиц, а само взаи модействие учитывается с помощью так называемого «обменно-корреляци онного» функционала, зависящего от электронной плотности. Такой прием показал хорошие результаты и на данный момент все существующие методы расчета на основе МФП используют это приближение. Функционал энергии 𝐸[𝑛] в описываемом приближении имеет вид: ∫︁ (3.6) 𝐸𝐾𝑆 = 𝑇𝑠 [𝑛] + 𝑑r𝑉𝑒𝑥𝑡 (r)𝑛(r) + 𝐸𝐻𝑎𝑟𝑡𝑟𝑒𝑒 [𝑛] + 𝐸𝐼𝐼 + 𝐸𝑥𝑐 [𝑛], где электронная плотность определяется выражением 𝑛(r) = ∑︁ 𝜎 𝑛(r, 𝜎) = 𝜎 𝑁 ∑︁ ∑︁ 𝜎 𝑖=1 |𝜓𝑖𝜎 (r)|2 , кинетическая энергия в приближении невзаимодействующих частиц: 𝜎 𝜎 𝑁 𝑁 1 ∑︁ ∑︁ 𝜎 2 𝜎 1 ∑︁ ∑︁ 𝑇𝑠 = − ⟨𝜓𝑖 | ∇ |𝜓𝑖 ⟩ = |∇𝜓𝑖𝜎 |2 , 2 𝜎 𝑖=1 2 𝜎 𝑖=1 а энергия кулоновского взаимодействия электронов ∫︁ ′ 1 3 3 ′ 𝑛(r)𝑛(r ) Hartree-def} 𝐸𝐻𝑎𝑟𝑡𝑟𝑒𝑒 [𝑛] = 𝑑 𝑟𝑑 𝑟 . 2 |r − r′ | ft:E_xc-def} (3.7) (3.8) (3.9) Здесь 𝜎 — суммарная проекция спина системы, 𝑁 𝜎 — число состояний при заданной проекции спина 𝜎, характеризующихся волновыми функциями 𝜙𝜎𝑖 (r) и собственными значениями энергии 𝜀𝜎𝑖 . Обменно-корреляционный функционал содержит многочастичные об менные и корреляционные эффекты, а также часть кинетической энергии, связанной со взаимодействием: 𝐸𝑥𝑐 [𝑛] = ⟨𝑇ˆ⟩ − 𝑇𝑠 [𝑛] + ⟨𝑉ˆ𝑖𝑛𝑡 ⟩ − 𝐸𝐻𝑎𝑟𝑡𝑟𝑒𝑒 [𝑛]. (3.10)

equations-1} 33 Уравнения для определения волновых функций 𝜙𝜎𝑖 (r) и собственных значений энергии 𝜀𝜎𝑖 могут быть получены с помощью минимизации 𝐸𝐾𝑆 [𝑛] ⟨︀ ⃒ ′ ⟩︀ при условии 𝜓𝑖𝜎 ⃒𝜓𝑗𝜎 = 𝛿𝑖𝑗 𝛿𝜎𝜎′ [20]: equations-2} 𝜎 𝑉𝐾𝑆 𝜎 (𝐻𝐾𝑆 − 𝜀𝜎𝑖 )𝜓𝑖𝜎 (r) = 0, (3.11) 1 𝜎 𝜎 𝐻𝐾𝑆 (r) = − ∇2 + 𝑉𝐾𝑆 (r) (3.12) 2 𝛿𝐸𝐻𝑎𝑟𝑡𝑟𝑒𝑒 𝛿𝐸𝑥𝑐 = 𝑉𝑒𝑥𝑡 (r) + + = 𝑉𝑒𝑥𝑡 (r) + 𝑉𝐻𝑎𝑟𝑡𝑟𝑒𝑒 (r) + 𝑉𝑥𝑐𝜎 (r). 𝛿𝑛(r, 𝜎) 𝛿𝑛(r, 𝜎) (3.13) equations-3} ft:E_xc-LDA} Уравнения (3.11)-(3.13) называются уравнениями Кона-Шэма. При извест ном точном выражении для обменно-корреляционного функционала 𝐸𝑥𝑐 [𝑛] решение системы (3.11)-(3.13) дает точные значения энергии и электронной плотности основного состояния. Обменно-корреляционный функционал 𝐸𝑥𝑐 [𝑛] может быть с хорошей точностью аппроксимирован локальным или почти локальным функционалом плотности (LDA): 𝐸𝑥𝑐 [𝑛] = ∫︁ 𝑑r𝑛(r)𝜀𝑥𝑐 ([𝑛], r), (3.14) где 𝜀𝑥𝑐 ([𝑛], r) — энергия на один электрон в точке r, зависящая только от плотности 𝑛(r, 𝜎) в некоторой окрестности точки r. 3.2. Метод функционала плотности для модели ВЭГ Традиционные реализации МФП обычно подразумевают в качестве внеш него потенциала потенциал точечных ионов. Для ВЭГ вместо этого использу ется однородный несжимаемый фон положительных зарядов. В модели ВЭГ электроны не могут образовывать связанных состояний, поэтому для описа ния волновых функций электронов логично использовать базис плоских волн [21]. Рассмотрим кубическую ячейку c объёмом 𝑉 , в которой содержится 𝑁 электронов на однородном несжимаемом компенсирующем фоне положи

matrix_form} 34 тельного заряда; на ячейку наложены периодические граничные условия. Га мильтониан такой системы определяется формулой (2.40) или, с учетом изо тропности, (2.51). Волновые функции в уравнениях Кона–Шэма (3.11)–(3.13) записываются в виде разложения (индекс 𝜎 далее опущен): 𝜓𝑖 (r) = ∑︁ q ∑︁ 1 𝑐𝑖,q × √ exp(iq · r) ≡ 𝑐𝑖,q × |q⟩, 𝑉 q (3.15) где 𝑐𝑖,q — коэффициенты разложения волновых функций в базисе ортонор мированных плоских волн: ∫︁ 1 ′ 𝑑r exp(−iq′ · r) exp(iq · r) = 𝛿q′ ,q . ⟨q |q⟩ ≡ 𝑉 (3.16) 𝑉 Так как ячейка периодически повторяется в пространстве, потенциал 𝑉𝐾𝑆 также является периодическим и может быть разложен в ряд Фурье: 𝑉𝐾𝑆 (r) = ∑︁ 𝑚 𝑉𝐾𝑆 (G𝑚 ) exp(iG𝑚 · r), где G𝑚 — векторы обратной решетки, ∫︁ 1 𝑉𝐾𝑆 (G) = 𝑉𝐾𝑆 (r) exp(−iG × r)𝑑r, 𝑉 (3.17) (3.18) 𝑉 а матричные элементы потенциала принимают вид: ′ ⟨q |𝑉𝐾𝑆 |q⟩ = ∑︁ 𝑚 𝑉𝐾𝑆 (G𝑚 )𝛿q′ −q,G𝑚 , (3.19) и не равны нулю только если q и q′ отличаются на некоторый вектор обратной решётки G𝑚 . Если теперь определить векторы q и q′ следующим образом: q = k+G𝑚 и q′ = k + G𝑚′ (G𝑚 и G𝑚′ отличаются на вектор обратной решетки), то уравнения Кона–Шэма для любого заданного k преобразуются к матричному уравнению на коэффициенты базисных функций: ∑︁ 𝑚′ 𝐻𝑚,𝑚′ (k)𝑐𝑖,𝑚′ (k) = 𝜀𝑖 (k)𝑐𝑖,𝑚 (k), (3.20)

35 где 1 𝐻𝑚,𝑚′ (k) = ⟨k + G𝑚 | − ∇2 + 𝑉𝐾𝑆 |k + G𝑚′ ⟩. 2 (3.21) Матричное уравнение (3.20) решается итерационными методами; в ре зультате для каждого k находятся собственные значения энергии 𝜀𝑖 (k) и ко эффициенты разложения одночастичных волновых функций в базисе плоских волн 𝑐𝑖,𝑚 . Далее можно рассчитать распределение электронной плотности 𝑛(r) и некоторые термодинамические функции. При нулевой температуре обыч но рассчитывается полная энергия, при ненулевой температуре — свободная энергия: из полной энергии вычитается 𝑇 𝑆, где [︃ ]︃ ∑︁ ∑︁ 𝑆=− 𝑓𝑖 ln 𝑓𝑖 + (1 − 𝑓𝑖 ) ln(1 − 𝑓𝑖 ) , 𝑖 (3.22) 𝑖 𝑓𝑖 — числа заполнения, определяемые функцией Ферми–Дирака для уровня с энергией 𝜀𝑖 . 3.3. Алгоритм моделирования ВЭГ методом функционала плотности при конечной температуре Для моделирования ВЭГ необходимо найти способ генерации конфи гураций электронов, по которым затем будет производиться усреднение для расчета термодинамических функций. Для классических частиц (ионов) для этой цели можно воспользоваться классическими уравнениями Ньютона: по теореме Фейнмана–Гельмана по заданному распределению электронной плот ности можно рассчитать силы, действующие на ионы, а затем проинтегриро вать уравнения движения с заданным временным шагом [21]. Для квантовых частиц (электронов) такой способ не подходит. Предлагается следующий ал горитм. Перед началом моделирования электроны можно разместить в узлах гранецентрированной кубической ячейки; такая конфигурация обладает до

36 статочно малой энергией, кроме того, соседние электроны при таком располо жении не смогут находиться на малом расстоянии друг от друга. Далее можно использовать алгоритм Метрополиса [22]: выбирается случайный электрон и сдвигается на случайный вектор в пространстве. Новая конфигурация при нимается с вероятностью {︂ (︂ )︂}︂ ∆𝐸 𝑊 = min 1, exp − , 𝑇 (3.23) где ∆𝐸 = 𝐸new − 𝐸old , 𝐸old — энергия старой конфигурации, 𝐸new — энер гия новой конфигурации, 𝑇 — температура моделирования. Для вычисления энергии необходимо решить уравнения Кона–Шэма. Алгоритм решения уравнений Кона–Шэма (3.11)–(3.13) для ВЭГ выгля дит следующим образом: 1. Выбирается начальное приближение для электронной плотности 𝑛(r) текущей конфигурации. 2. Вычисляется потенциал 𝑉𝐾𝑆 (r) = 𝑉𝑒𝑥𝑡 (r) + 𝑉𝐻𝑎𝑟𝑡𝑟𝑒𝑒 [𝑛] + 𝑉𝑥𝑐 [𝑛]. 3. Решается задача на собственные значения (3.20), определяются соб ственные уровни энергии и коэффициенты разложения по базису для каждого значения k. 4. Вычисляется новая электронная плотность 𝑛1 (r) = ∑︀ 2 𝑖 𝑓𝑖 |𝜓𝑖 (r)| . 5. Если ||𝑛1 (r)−𝑛(r)|| > 𝜀𝑐 , то 𝑛(r) = 𝛼𝑛1 (r)+(1−𝛼)𝑛(r) и возвращаемся к пункту 2. 6. Иначе вычисляем требуемые свойства системы, в том числе полную энергию. Здесь 𝜀𝑐 — заданная невязка для электронной плотности, 𝛼 — так назы ваемый коэффициент линейного смешивания. Часто используют более общий метод смешивания, подразумевая, что 𝛼 не является константой, а зависит от значений электронной плотности на предыдущих итерациях.

37 Заключение В работе получены следующие основные результаты: 1. Получены аналитические выражения для вторых производных термоди намического потенциала и термодинамических коэффициентов ИФГ. 2. Получены асимптотические выражения для всех термодинамических функций ИФГ в пределе низких и высоких температур. 3. Разработана общедоступная программная реализация для модели ИФГ на языке Python. 4. Получен гамильтониан изотропной модели ВЭГ. 5. Разработан алгоритм программной реализации МФП для модели ВЭГ.

38 Список литературы 1. Pauli W. On the connexion between the completion of electron groups in an atom with the complex structure of spectra // Zeitschrift für Physik. — 1925. — т. 31. — с. 765. 2. Pauli W. The connection between spin and statistics // Physical Review. — 1940. — т. 58, № 8. — с. 716. 3. Ландау Л. Д., Лифшиц Е. М. Теоретическая физика. Статистическая физика. т. 5. ч. 1. — Москва : Физматлит, 2018. — с. 616. 4. Киржниц Д. А., Лозовик Ю. Е., Шпатаковская Г. В. Статистическая модель вещества // Усп. физ. наук. — 1975. — т. 117, № 9. — с. 3—47. 5. Thermal contribution to thermodynamic functions in the Thomas–Fermi model / O. Shemyakin [и др.] // Journal of Physics A: Mathematical and Theoretical. — 2010. — т. 43, № 33. — с. 335003. 6. fdint library. — URL: https://pypi.org/project/fdint. 7. IFG module. — URL: https://pypi.org/project/ifg. 8. MacIver D., Hatfield-Dodds Z., Contributors M. Hypothesis: A new approach to property-based testing // Journal of Open Source Software. — 2019. — 21 нояб. — т. 4, № 43. — с. 1891. — ISSN 2475-9066. 9. Фортов В. Е., Храпак А. Г., Якубов И. Т. Физика неидеальной плазмы. Учебное пособие. — Москва : Физматлит, 2004. — с. 528. — ISBN 5-9221-0173-0. 10. Slattery W. L., Doolen G. D., DeWitt H. E. Improved equation of state for the classical one-component plasma // Phys. Rev. A. — 1980. — июнь. — т. 21, вып. 6. — с. 2087—2095. 11. Wigner E. On the Interaction of Electrons in Metals // Phys. Rev. — 1934. — дек. — т. 46, вып. 11. — с. 1002—1011.

39 12. Ewald P. P. Die Berechnung optischer und elektrostatischer Gitterpotentiale. — 1921. — янв. 13. Kolafa J., Perram J. W. Cutoff Errors in the Ewald Summation Formulae for Point Charge Systems // Molecular Simulation. — 1992. — т. 9, № 5. — с. 351—368. 14. Mahan G. D. Many Particle Physics, Third Edition. — New York : Plenum, 2000. 15. Блинов В. Н. Дальнодействующие взаимодействия в компьютерном мо делировании систем в конденсированном состоянии // Наноструктуры. Математическая физика и моделирование. — 2014. — т. 10, № 1. — с. 5—27. 16. Yakub E., Ronchi C. An efficient method for computation of long-ranged Coulomb forces in computer simulation of ionic fluids // The Journal of Chemical Physics. — 2003. — т. 119, № 22. — с. 11556—11560. 17. Rapaport D. C. The Art of Molecular Dynamics Simulation. — 2-е изд. — Cambridge University Press, 2004. 18. Hohenberg P., Kohn W. Inhomogeneous Electron Gas // Phys. Rev. — 1964. — нояб. — т. 136, 3B. — B864—B871. — URL: https://link.aps.org/ doi/10.1103/PhysRev.136.B864. 19. Mermin N. D. Thermal Properties of the Inhomogeneous Electron Gas // Phys. Rev. — 1965. — март. — т. 137, 5A. — A1441—A1443. — URL: https://link.aps.org/doi/10.1103/PhysRev.137.A1441. 20. Kohn W., Sham L. J. Self-Consistent Equations Including Exchange and Correlation Effects // Phys. Rev. — 1965. — нояб. — т. 140, 4A. — A1133— A1138. — URL: https://link.aps.org/doi/10.1103/PhysRev.140. A1133.

40 21. Martin R. M., Martin R. M. Electronic structure: basic theory and practical methods. — Cambridge university press, 2004. 22. Замалин В., Норман Г., Филинов В. Метод Монте-Карло в статистиче ской термодинамике. — Москва : Наука, 1977.

41 Словарь терминов FFT Fast Fourier transform. 21 LDA Local density approximation. 33 PyPi Хранилище пакетов для Python. 14, 42 Python Скриптовый язык программирования общего назначения. 13, 14, 37 БКА Большой канонический ансамбль. 31 БПФ Быстрое преобразование Фурье. 21 ВЭГ Взаимодействующий электронный газ. 2, 21, 26, 29, 30, 33, 35–37 ИБГ Идеальный больцмановский газ. 11 ИФГ Идеальный Ферми-газ. 2, 3, 6, 8, 11, 13–15, 37, 42, 43 МФП Метод функционала плотности. 4, 30, 32, 33, 37 ОКП Однокомпонентная плазма. 16–18 интеграл Ферми-Дирака Функция, определяемая выражением (1.4). 7, 13

42 Приложение А Пример расчета термодинамических свойств ИФГ Модуль легко ставится через PyPi: $ pip install numpy $ pip install ifg Основной функционал содержится в классе IfgCalculator. from ifg import IfgCalculator import numpy as np specific_volumes = [0.33, 1., 10.,] temperatures = np.linspace(1e-4, 1e4, 10**6) calculator = IfgCalculator( temperatures=temperatures, specific_volumes=specific_volumes, input_in_si=False, # False даст в атомной системе output_in_si=False ) После этого все свойства будут доступны как поля экземпляра calculator: >>> calculator.p array([[1.21465823e+01, 1.91415600e+00, 4.12392425e-02], [1.21466330e+01, 1.91419106e+00, 4.12555143e-02], [1.21467832e+01, 1.91429487e+00, 4.13036646e-02], ..., [3.03030975e+04, 9.99999392e+03, 9.99998139e+02], [3.03031278e+04, 1.00000039e+04, 9.99999139e+02],

43 [3.03031581e+04, 1.00000139e+04, 1.00000014e+03]]) >>> calculator.C_P array([[4.92449475e-05, 1.03122265e-04, 4.78651124e-04], [4.97374654e-03, 1.04153735e-02, 4.83453163e-02], [9.89825789e-03, 2.07277142e-02, 9.62203561e-02], ..., [2.49998418e+00, 2.49999478e+00, 2.49999948e+00], [2.49998418e+00, 2.49999478e+00, 2.49999948e+00], [2.49998418e+00, 2.49999478e+00, 2.49999948e+00]]) Во втором примере видим, как теплоемкость 𝐶𝑃 стремится к 5/2 при высоких температурах — ожидаемое поведение для ИФГ. Можно также проверить, что давление стремится к конечному пределу при низких температурах. Изменим входные данные: temperatures = np.linspace(1e-10, 1e-8, 10**4) calculator = IfgCalculator( temperatures=temperatures, specific_volumes=specific_volumes, input_in_si=False, output_in_si=False ) А затем запустим расчет давления: >>> calculator.p array([[12.14658227, 1.914156 , 0.04123924], [12.14658227, 1.914156 , 0.04123924], [12.14658227, 1.914156 , 0.04123924], 1.914156 , 0.04123924], ..., [12.14658227,

44 [12.14658227, 1.914156 , 0.04123924], [12.14658227, 1.914156 , 0.04123924]]) Получены конечные пределы. Больше примеров можно найти по ссылке https://ifg-py.readthedocs. io/en/latest/class_desc.html. Код графиков 1.1, 1.2, 1.3 доступен в репо зитории: https://github.com/alekseik1/ifg-py/tree/master/examples/ plots.