Исследование неидеальной электрон-ионной плазмы методом динамики волновых пакетов

Во Введении обоснована актуальность диссертационной работы, сфор­ мулирована цель и аргументирована научная новизна исследований, показана практическая значимость полученных результатов, представлены выносимые на защиту научные положения.
В первой главе проводится обзор литературы по имеющимся экспери­ ментальным методам исследования неидеальной плазмы. Помимо этого, раздел включает описание методов атомистического моделирования, применяемых для теоретического изучения разогретого плотного вещества, таких как классиче­ ская молекулярная динамика, динамика волновых пакетов, теория функциона­ ла плотности и метод Монте-Карло в терминах интегралов по траекториям.
Метод классической молекулярной динамики (МД) основывается на ре­ шении уравнений движения Ньютона для точечных частиц, представляющих собой молекулы, атомы, электроны и ионы, взаимодействие между которыми определяется потенциалом взаимодействия. В случае неидеальной плазмы в ка­ честве потенциалов используют функции, совпадающие с потенциалом Кулона на больших расстояниях и имеющих финитные значения энергии взаимодей­ ствия вблизи нуля. Такой подход в классической МД успешно применяется для исследования свойств неидеальной плазмы в области низких плотностей и боль­ ших температур. Недостатками классической МД является низкая точность учета квантово-механических эффектов, ограниченная предположениями, сде­
ланными при конструировании потенциалов взаимодействия.
Дальнейшим развитием классической МД в сторону более высокой точ­ ности моделирования межчастичного взаимодействия в плазме является метод молекулярной динамики волновых пакетов (МДВП). В этом методе рассматри­ вается система точечных классических ионов и квантовых электронов, волно­ вая функция которых параметризована модельной функцией. Наиболее частым вариантом параметризации является изотропный Гауссовый волновой пакет пе­ ременной ширины, характеризуемый положением центра пакета ⃗r, его шириной
s и сопряженными с ними значениями импульсов p⃗ и ps соответственно 3 3/4 3ips 2i
φ(⃗x)= 2 exp − 2πs
4s
2 − (x−r) + p·(x−r) . (1) 2 s
Эволюция параметров волновых пакетов следует из вариационного прин­ ципа и в некоторых случаях она может быть описана уравнениями гамильтоно­ вой динамики для параметров волновых пакетов. Так, в случае конструирова­ ния многоэлектронной волновой функции как произведения волновых пакетов в приближении Хартри, уравнения, описывающие эволюцию во времени поло­ жения и ширины волнового пакета, а также их сопряженных импульсов, имеют вид
r ̇ =∂H, p ̇ =−∂H, s ̇ =∂H, p ̇ =−∂H. (2) k ∂p k ∂r k ∂p sk ∂s
k k sk k
Данный подход позволяет избавиться от необходимости задавать потен­ циал взаимодействия между частицами и повышает точность описания элек­ тронной подсистемы, но такая модель все еще имеет проблемы с описанием связанных состояний атомов и молекул с несколькими электронами, а также состояний с высокой плотностью электронов, когда эффектами вырождения электронного газа нельзя пренебречь. Более сложные приближения (антисим­ метризованные волновые пакеты — АМДВП, Electron Force Field — eFF) повы­ шают точность описания квантовых явлений в веществе, но требуют большого объема вычислений для получения статистически значимых результатов.
Несмотря на большое разнообразие вариаций молекулярной динамики вол­ новых пакетов, получение численного метода с минимальным количеством при­ ближений и хорошей производительностью все еще остается актуальной зада­ чей, так как существующие варианты являются либо вычислительно сложными (АМДВП), либо требуют ввода эмпирических функционалов (eFF). Более эф­ фективным может быть подход, основанный на объединении метода МДВП в приближении Хартри и теории функционала плотности, разработка которого являлась одной из целей настоящей работы.
В первой главе приведены также описания первопринципных подходов, основанных на теории функционала плотности (DFT) и методе Монте-Карло в терминах интегралов по траекториям (PIMC). Эти подходы обладают высо­ кой точностью описания квантово-механических явлений, но требуют большого

⟨ne(r)⟩V/8,t = V8
объема вычислений и не способны моделировать неадиабатическую динамику электронов и ионов.
Во второй главе представлен новый поход к описанию системы электро­ нов в плазме, основанный на динамике волновых пакетов и теории функционала плотности — МДВП-ФП. В работе приведено теоретическое обоснование это­ го подхода, причем особое внимание уделено граничным условиям, к которым семейство методов динамики волновых пакетов очень чувствительно. В допол­ нение предложены способы применения метода МДВП-ФП для неравновесного расчета изоэнтропы сжатия и ударной адиабаты вещества.
Использование МДВП в приближении Хартри дает легко решаемые урав­ нения динамики волновых пакетов, а учет обменных и корреляционных эффек­ тов проводится в рамках теории функционала плотности путем добавления в гамильтониан моделируемой системы H дополнительного члена, выражающего обменно-корреляционную энергию через электронную плотность
Hwpmd−dft = H + (Ts[n] − Ts[ni]) + (Exc[n] − Exc[ni]), (3) ii
где n — суммарная локальной плотность электронов, ni = φi(⃗r)φ∗i (⃗r) — вклад в электронную плотность от i-го волнового пакета, Ts[n] — функционал, опи­ сывающий кинетическую энергию, Exc[n] — обменно-корреляционный функци­ онал, аналогичный используемым в DFT. В качестве обменно-корреляционного функционала в перелагаемом методе предлагается использовать приближение локальной электронной плотности (LDA) как дающее оптимальный баланс точ­ ности и сложности вычислений.
В качестве граничных условий для метода предложены отражающие стен­ ки (рис. 1). В отличие от периодических граничных условий, типично применя­ емых для классической МД, отражающие граничные условия имеют ряд недо­ статков, однако позволяют избежать неограниченного уширения волновых па­ кетов для слабовзаимодействующих электронов [21].
Использование отражающего потенциала в качестве граничных условий приводит к возникновению внешнего воздействия на моделируемую систему, которое искажает профили электронной плотности и вводит систематическую погрешность в термодинамические параметры изучаемой системы. Минимизи­ ровать возникающие эффекты можно путём увеличения размеров моделируе­ мого вещества или же расчётом интересующих параметров внутри ячейки мо­ делирования, вдали от приграничного слоя. Дополнительно из-за искажения профиля электронной плотности предлагается в процессе расчёта вычислять среднюю электронную плотность только во внутренней области моделируемой ячейки
Lx/4 Ly/4 Lz/4 Ne
φk(r)φ∗k(r) dr. (4) −Lx/4 −Ly/4 −Lz/4 k=1
Рис. 1. Ячейка моделирования с ограничивающим потенциалом (темно-красные кривые). Средняя электронная плотность для равномерно распределенных электронов показана синей пунктирной линией, зеленая пунктирная линия — пример профиля электронной плотности, полученного из моделирования. На рисунке также показаны соотношения средней ширины волновых пакетов в центе и на границе ячейки.
Выбранный тип граничных условий позволяет выполнять расчёт некото­ рых термодинамических параметров, таких как давление, вычисление которо­ го с использованием вириальной теоремы не вполне обосновано из-за наличия внешнего воздействия на систему от граничных условий.
В третьей главе описано построение численного алгоритма на базе раз­ работанного метода МДВП-ФП. Представлено подробное описание алгоритмов расчета обменно-корреляционной энергии с использованием регулярных и адап­ тивных пространственных сеток. Приведен алгоритм построения адаптивной сетки с изменяемым шагом, опирающийся на параметры волновых пакетов. Описаны основные оптимизационные подходы, применяемые для увеличения быстродействия. В заключение представлены результаты тестов производитель­ ности и масштабирования метода, выполнено сравнение характерных времён работы алгоритмов МДВП-ФП и МКВП-ФП (аналог МДВП-ФП с использова­ нием алгоримта Монте-Карло вместо динамики пакетов) с аналогичными рас­ четами, проведенными с помощью квантовой молекулярной динамики.
Расчет обменно-корреляционного функционала и его производных в мето­ де МДВП-ФП выполняется путем численного интегрирования по трёхмерной пространственной сетке, содержащей значения локальной плотности электро­ нов, полученные из координат и ширин волновых пакетов. Пространственные сетки могут быть с регулярным шагом или же подстраивающимся к гради­ енту электронной плотности — адаптивные сетки. Последние являются более эффективными, так как минимизируют количество избыточных вычислений. Построение адаптивной сетки выполняется, основываясь на параметрах волно­ вых пакетов и не требует традиционной для данного типа сеток процедуры

12 разбиения и объединения ячеек регулярной сетки.
Рис. 2. Время выполнения одного шага алгоритмов МКВП и МКВП-ФП в зависимости от числа частиц (a) и ускорение от количества вычислительных ядер при Ne = 1024 (b). Крас­ ные кружки соответствуют оригинальному методу МКВП, синие квадраты и треугольники — методу МКВП-ФП с регулярной сеткой, выполняющемуся на ЦП (квадраты) и ГУ (тре­ угольники). Голубые квадраты и треугольники соответствуют методу МКВП-ФП с адаптив­ ной пространственной сектой с использованием ЦП и ГУ для вычислений соответственно. Верхняя и нижняя штриховые линии на рисунке (а) отображают линейную и квадратичную зависимости от числа частиц. Штриховая линия на рисунке (б) — теоретическая зависи­ мость ускорения R = (p2(N − 1))/(p(2N − 1) − N) от числа ядер p, штрихпунктирная линия — идеальная кривая эффективности распараллеливания R = p.
Использование графических ускорителей (ГУ) и суперкомпьютеров на их базе позволяет существенно увеличить производительность расчетов для мно­ гих алгоритмов, основанных на решении уравнений движения системы частиц и обеспечить ускорение вычислений на несколько порядков в сравнении с цен­ тральным процессором (ЦП). В случае МКВП-ФП, как показано на рис 2а, сово­ купное применение адаптивных сеток и ГУ позволило ускорить вычисления на 3−4 порядка. Также рисунок демонстрирует, что время расчета для алгоритма МКВП-ФП растет с увеличением числа частиц не быстрее, чем O(Ne2).
Методы МДВП-ФП и МКВП-ФП были адаптированы для вычислений как на ГУ, так и на многопроцессорных вычислительных кластерах. При создании параллельной версии в качестве основы использовался пакет МД моделирова­ ния LAMMPS. Как видно из рисунка 2б, при использовании 128 вычислитель­ ных узлов расчеты ускоряются в 40 раз по сравнению с выполнением програм­ мы на одном узле. Ограничением эффективности масштабирования является часть алгоритма, ответственная за расчет кулоновского взаимодействия, кото­ рая требует пересылки большого количества данных между процессами.
В четвертой главе представлены результаты моделирования неидеаль­ ной водородной и дейтериевой плазмы методом МДВП-ФП и сравнения с ме­ тодами классической МД, МДВП в приближении Хартри, а также PIMC. На
основании этих расчётов определены границы применимости методов МДВП и МДВП-ФП. Помимо этого, в главе дано описание расчетов, посвященных по­ строению ударной адиабаты дейтерия методом МДВП-ФП и их сравнение с экспериментальными данными и другими теоретическими моделями.
Результатом добавления обменно-корреляционной энергии в метод МДВП является улучшение точности описания связанных молекулярных состояний. Так, из-за введения дополнительной энергии метод МДВП-ФП наиболее точно из всего семейства методов МДВП воспроизводит энергию диссоциации моле­ кулы водорода, расчетное значение которой составляет 0.155 Eh по сравнению с точным значением 0.166 Eh (рис. 3).
Рис. 3. Энергия связи молекулы H2 в зависимости от межатомного расстояния, рассчитан­ ная с использованием исходного алгоритма МДВП (пунктирная красная линия), МДВП c коррекцией на энергию изолированного состояния двух атомов H (сплошная красная линия), eFF с коррекцией (пунктирная зеленая линия) и МДВП-ФП с коррекцией (сплошная синяя линия). Точное значение показано черной сплошной линией.
В данной работе метод МДВП-ФП использовался для расчета энергии во­ дородной плазмы в диапазоне температур T = 104 − 5 · 104 К и концентра­ ций электронов ne = 1020 − 5 · 1024 см−3. На рис. 4а результаты расчетов для T = 3 · 104 К сопоставлены с расчетами методами классической МД, МДВП без антисимметризации, eFF и PIMC. Как видно из рисунка, учет обменно­ корреляционной энергии в методе МДВП-ФП позволил получить хорошее ка­ чественное согласие с первопринципными расчетами методом PIMC в области высоких плотностей ne = 1022 −1023 см−3, где метод МДВП, как и классическая МД, оказался неприменимым. Разница абсолютных значений энергий в методах МДВП-ФП и PIMC может быть связана с различным способом учета энергии связанных состояний электронов и ионов.
Рис. 4. Внутренняя энергия водородной плазмы на атом (a) и адиабата ударного сжатия (адиабата Гюгонио) дейтерия (b). На рисунке (а) синяя кривая соответствует значениям энергии, полученным в МКВП-ФП, коричневая — МКВП в приближении Хартри, зеленые треугольники — классическая МД с потенциалом Кельба [22], фиолетовые крестики — eFF, черная кривая — PIMC. Вертикальная пунктирная линия соответствует границе вырожде­ ния neλ3B = 1. Красная штриховая линия показывает область некорректных результатов для метода МКВП. На рисунке (b) результаты моделирования методом МДВП-ФП представлены красными кружками, расчеты методом АМДВП [23] — черные крестики, расчеты методами квантового Моте-Карло [24] и Монте-Карло в терминах интегралов по траекториям [25] пред­ ставлены треугольниками и черными звездами соответственно. Экспериментальные данные Сандийской национальной лаборатории [5,26] отражены в виде ромбов зеленого и бирюзового цвета. Эксперименты по лазерному сжатию вещества [6] представлены синими квадратами.
Ударная адиабата дейтерия, рассчитанная методом МДВП-ФП, в отличие от результатов МДВП с антисимметризацией с хорошей точностью воспроизво­ дит предельную степень сжатия, наблюдаемую в эксперименте, а также описы­ вает переход от молекулярного к атомарному состоянию дейтерия. Причиной завышения давления в области низких температур может быть недостаточная точность параметризации волновой функции электронов и используемого об­ менно-корреляционного потенциала, что ограничивает применимость МДВП­ ФП для низких температур (T < 5000 К). В пятой главе продемонстрировано применение метода МДВП-ФП для исследования динамических и релаксационных процессов: изоэнтропического сжатия и электрон-ионной релаксации в неизотермической плазме. Для построения изоэнтропы дейтерия был применен оригинальный под­ ход, заключающийся в ступенчатом сжатии вещества с контролируемой скоро­ стью сжатия. Выбор последней основывался на необходимости обеспечения изо­ энтропийности процесса и проводился путем исследования сходимости кривой P(ρ) при уменьшении скорости сжатия. В дополнение к динамическому под­ ходу, также был применен традиционный метод Ферми-Зельдовича. Кривые, полученные методом МДВП-ФП (рис. 5) с хорошей точностью согласуются с данными, полученными в экспериментах по ударному сжатию, хотя значения температур, предсказанных методом МДВП-ФП, являются несколько занижен­ ными в сравнении с другими теоретическими моделями. Рис. 5. Изоэнтропы ударно сжатого дейтерия в плоскости P − ρ, рассчитанные различными численными методами. Синие квадраты — статические расчеты МДВП-ФП для температур, рассчитанных методом QMD [27], зеленая и красная линии — динамические расчеты МДВП­ ФП для различных начальных точек. Коричневая штрихпунктирная кривая — численное решение уравнения Ферми-Зельдовича на сетке термодинамических параметров, определен­ ных с помощью МДВП-ФП (ρ0 = 1.09 г/см3, T0 = 6900 K). Черная штриховая линия — расчеты методом QMD [27]. Черные ромбы и кружки соответствуют экспериментальным данным [28,29,29–32]. Верхняя и нижняя пунктирные линии показывают давление идеаль­ ного газа Больцмана и Ферми соответственно. Описанный динамический подход обеспечивает расчет изоэнтропы, исполь­ зуя лишь одну траекторию системы, поэтому скорость его работы значительно превосходит метод, основанный на расчете термодинамических параметров на сетке ρ−T с последующим численным решением уравнения Ферми-Зельдовича. Постановка численного эксперимента для моделирования электрон-ион­ ной релаксации соответствовала работе [33]. Характерное время температур­ ной релаксации, полученное методом МДВП-ФП при ne = 1022 см−3, составило 58 ± 8 фс для плазмы с модельным соотношением масс ионов и электронов mi/me = 20, и 250 ± 20 фс для mi/me = 100 (рис. 6). В данном случае методы МДВП и МДВП-ФП показали близкие результаты, что подтверждает возмож­ ность применения МДВП-ФП для исследования неравновесных систем и про­ цессов. Напомним, что при более высоких плотностях плазмы метод МДВП ока­ зывается неприменимым, а методы квантовой молекулярной динамики (QMD) и PIMC не имеют возможности моделирования неадиабатической динамики электронов и ионов. Таким образом, метод МДВП-ФП остается практически единственным подходом для моделирования релаксационных и динамических процессов при ne > 1022 см−3.
3.2 3 2.8 2.6 2.4 2.2 2 1.8 1.6
1.4
0 100
200 300
T, eV
(a)
t/τpl
∆T, eV
(b)
3
t/τpl
2.5 2
1.5
0.5
0 10 20 30 40 50 60 7075
Рис. 6. (a) Зависимость температур электронов (1) и ионов (2) от времени, полученная методом МДВП в приближении Хартри, для концентрации электронов ne = 1022 см−3 и модельного соотношения масс ионов и электронов mi/me = 20. (b) Эволюция разности тем­ ператур электронов и ионов от времени: МДВП-ФП для mi/me = 20 (1) и mi/me = 100 (2), МДВП в приближении Хартри для mi/me = 20 (3). Сплошные линии — экспоненциальная аппроксимация.
ОСНОВНЫЕ РЕЗУЛЬТАТЫ И ВЫВОДЫ
1. Предложена модель неидеальной электрон-ионной плазмы МДВП-ФП, ос­ нованная на представлении электронов в форме гауссовских волновых пакетов и учете обменно-корреляционного взаимодействия, выраженного через функционал электронной плотности. Показано, что использование отражающих граничных условий позволяет решить проблему расплыва­ ния волновых пакетов, а влияние этого типа граничных условий на изме­ ряемые характеристики плазмы можно учесть за счет анализа профиля электронной плотности внутри расчетной ячейки.
2. На базе модели МДВП-ФП разработан метод компьютерного моделирова­ ния неидеальной плазмы и плотного разогретого вещества. Показано, что расчет обменно-корреляционного взаимодействия может быть ускорен за счета использования адаптивной пространственной сетки, а также приме­ нения графических ускорителей (ГУ). Совокупное применение адаптив­ ных сеток и ГУ позволяет сократить время вычислений на 3 − 4 порядка величины, а интеграция алгоритма МДВП-ФП в пакет LAMMPS позволя­ ет выполнять расчеты на суперкомпьютерных кластерах с ускорением до 40 раз на 128 вычислительных узлах. Асимптотика зависимости времени расчета от числа частиц ограничена квадратичной функцией.
3. С помощью метода МДВП-ФП исследованы равновесные свойства неиде­ альной плазмы водорода и дейтерия. На примере молекулы водорода по­ казано, что учет обменно-корреляционного взаимодействия улучшает точ­ ность моделирования многоэлектронных систем. Расчет уравнения состо­ яния неидеальной плазмы водорода показывает преимущества метода в области концентраций электронов 1022 − 1024 см−3 и дает качественное совпадение с расчетами более высокого уровня. Расчет ударной адиабаты дейтерия демонстрирует хорошее согласие с имеющимися эксперименталь­ ным данными.
4. На примере расчетов изоэнтропы сжатия плазмы дейтерия и электрон­ ионной релаксации в неизотермической неидеальной плазме водорода про­ демонстрирована способность метода МДВП-ФП моделировать динамиче­ ские и неравновесные процессы вне адиабатического приближения. Дан­ ная способность, в совокупности с удовлетворительным согласием с имею­ щимися экспериментальными и теоретическими данными по равновесным свойствам неидеальной плазмы, делает МДВП-ФП уникальным расчет­ ным методом в области высоких плотностей, где подходы классической молекулярной динамики и метода волновых пакетов без учета обменно­ корреляционных эффектов оказываются неприменимыми.