Размещено на http://www.allbest.ru/

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ

«ВЫСШАЯ ШКОЛА ЭКОНОМИКИ»

Московский институт электроники и математики

Департамент прикладной математики

Курсовая работа

Анализ результатов молекулярно-динамического моделирования электронных колебаний в кластерной наноплазме

Рецензент

В.В. Стегайлов

Научный руководитель

И.В. Морозов



Оглавление

Введение

1. Метод молекулярной динамики

1.1 Решение уравнений движения

1.2 Потенциалы взаимодействия

2. Оптимизация моделирования систем заряженных частиц

2.1 Приближенные методы моделирования систем заряженных частиц

2.2 Расчет взаимодействия частиц методом PPPM

2.3 Применение графических ускорителей

3. Моделирование неидеальной плазмы, образованной в результате ионизации наноразмерных кластеров

3.1 Постановка задачи

3.2 Численная модель моделирования кластерной наноплазмы

3.3 Спектры электронных колебаний в ионизованном нанокластере

4. Анализ колебаний электронов в кластерной наноплазме

4.1 Теоретическая оценка частоты электронных колебаний

4.2 Параметризация профиля электронной плотности

Выводы

Список цитированной литературы

Аннотация



Введение

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

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

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

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

Также в работе оценивается применимость метода Particle-Particle-Particle-Mesh для моделирования указанного типа систем заряженных частиц.

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



1. Метод молекулярной динамики

1.1 Решение уравнений движения

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

,

,

где - координаты частиц, - скорости частиц, - силы, действующие на частицы, - массы частиц.

Уравнения требуют задания начальных условий: координат частиц и их скоростей:

, , i = 1,…,N.

Затем уравнения движения интегрируются с заданным шагом по времени. Каждый шаг требуется вычисления сил, действующих на частицы, которые выражаются через заданный потенциал взаимодействия, определяющий свойства системы [1-3]:

,

Метод классической молекулярной динамики позволяет рассматривать системы, состоящие из более чем триллионов атомов на временах до нескольких пикосекунд. Однако, максимальный размер системы сильно зависит от сложности потенциала взаимодействия и применения различных методов оптимизации, которые будут рассмотрены в дальнейшем. Системы больших масштабов и на больших временах могут быть рассмотрены в рамках так называемой многомасштабной модели [12].

1.2 Потенциалы взаимодействия

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

В общем случае с ростом числа частиц сложность вычислений растёт как O(N2). Однако, для короткодействующих потенциалов существует много методов оптимизации вычислений, например, вычисление списков ближайших соседей. Суть метода в том, что для каждой частицы создаётся список, в который включаются частицы, находящиеся в определённом расстоянии вокруг неё. Считается, что с другими частицами она не взаимодействует. На каждом шаге проводятся расчёты взаимодействий только внутри списка, а сами списки пересчитываются через определённые интервалы по времени. Применение такого метода (часто в совокупности с методом связанных списков) позволяет добиться роста сложности как О(N).

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

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



2. Оптимизация моделирования систем заряженных частиц

2.1 Приближенные методы моделирования систем заряженных частиц

Для ускорения моделирования заряженных частиц, взаимодействующих посредством дальнодействующего куловноского потенциала, в больших системах существуют различные вычислительные методы. Например, метод сумм Эвальда [1, 13] очень популярен в современной молекулярной динамике, однако он применим только к периодическим системам. Этот метод является частным случаем формулы суммирования Пуассона, в которой суммирование энергий взаимодействия в реальном пространстве заменяется на эквивалентное суммирование в пространстве Фурье.

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

Метод PPPM (Particle-Particle, Particle-Mesh) содержит в себе плюсы точных молекулярно-динамических расчётов и вычислений на сетках [1]. Сила взаимодействия разделяется на два слагаемых: Fij = FijPP + FijPM , где РР (Particle-Particle) - непосредственное межчастичное взаимодействие, вычисляемая до определённого радиуса отсечки, а РМ (Particle-Mesh) считается достаточно гладкой функцией для вычисления на сетке. Вычисления на сетке производятся путем численного решения уравнения Пуассона в Фурье-пространстве.

Для проверки применимости метода PPPM к задаче моделирования кластерной наноплазмы был использован пакет программ LAMMPS (англ. Large-scale Atomic/Molecular Massively Parallel Simulator) [14]. Это свободно распространяемый пакет программ для классической молекулярной динамики, разработанный в Сандийской национальной лаборатории США. LAMMPS может применяться как для небольших расчётов, так и для изучения взаимодействия миллионов частиц. В этом пакете реализовано большое количество потенциалов и методов, применяемых для оптимизации вычислений, например, пересчёт списков ближайших соседей. Так же LAMMPS поддерживает параллельные вычисления. В частности, в пакете LAMMPS реализован метод PPPM. В данной работе был проведён тестовый расчет на LAMMPS с целью выяснения, подходят ли использованные в нем приближённые вычисления для задачи моделирования наноплазмы.

Также следует отметить метод TreeMD [1], использованный, в частности, и для задач моделирования кластерной наноплазмы15]. Тем не менее, результаты, полученные этим методом имеют явные отличия, от результатов, полученных точным МД моделированием [8,11]

2.2 Расчет взаимодействия частиц методом PPPM

При помощи LAMMPS было смоделировано движение двух отрицательно заряженных частиц с кулоновским потенциалом взаимодействия. За образец взят расчёт без применения PPPM. По траектории движения частиц построена зависимость потенциальной энергии от расстояния.

Для вычислений был использован потенциал coul/cut с отсечкой на расстоянии 75 атомных единиц. Во всех расчётах сделано 2000 шагов по 0.007фс.

При использовании PPPM на результат влияют несколько параметров. Одним из них является размер ячейки, в котором проводится расчёт. При уменьшении ячейки возрастает погрешность относительно исходного потенциала, так как меньшее количество частиц попадают в зону непосредственного расчёта взаимодействия. Для проверки взяты ячейки размером от 100 до 1000, измерялась средняя ошибка. Результаты представлены на рис. 2.2.

Рис Зависимость средней погрешности от размера ячейки

Ещё одним параметром, влияющим на результат, является точность Фурье-преобразования, указываемая в параметре kspace. Зависимость погрешности определения потенциала взаимодействия от этого параметра приведена на рис. 2.3. На рисунке видно, что в точке 0.0001 погрешность увеличивается. Это связано с автоматическим выбором размера сетки. При 0.0001 программа автоматически выбирает размер сетки равным 4 а.е., а при других параметрах больше. Это связано с тем, что параметр размера сетки должен быть кратным степеням чисел 2, 3, 5. То есть, если на основе остальных параметров программой автоматически был выбран неподходящий размер, резко возрастает погрешность.



Рис. 2.3 Зависимость средней погрешности от точности kspace.

Расчёты на PPPM занимали достаточно много времени: взаимодействие двух частиц с точностью 0. 001 занимало вплоть до 34.3807с. Это связано с тем, что изучаемая система не была периодической, частицы могут вылетать достаточно далеко и не должны возвращаться с другой стороны расчётной области, поэтому приходится брать размер области очень большим. При этом укрупнение сетки приводило к большой погрешности.

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

2.3 Применение графических ускорителей

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

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

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

Гибридные кластеры на основе GPU зачастую показывают быстродействие, во много раз превышающее быстродействие систем на универсальных центральных процессорах.

Несмотря на преимущества вычислений на GPU, нужно помнить, что создание программ, эффективно использующих их вычислительные возможности - это сложный процесс, требующий учета особенностей архитектуры конкретного устройства применительно к конкретной задаче, без которого производительность может оказаться недостаточно высокой. На рисунке 3.1 можно видеть, что наибольшая эффективность в сравнении с вычислениями на СPU при расчётах МД достигается на системах с большим количеством частиц, когда загружено большое количество ядер GPU, а так же в случае применения оптимизированного алгоритма [5].

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

Максимальное ускорение для алгоритма моделирования кластерной наноплазмы, достигнутое на GPU Nvidia Tesla M2050 по сравнению с одним ядром центрального процессора (CPU) Intel Xeon E5520, составило 100 раз (рис. 3.2). Как видно из рисунка, эффективность применения GPU существенным образом снижается при переходе к малому числу частиц N < 104. Это объясняется малой загрузкой GPU при использовании одного потока команд (thread) на частицу.

Таким образом, применение GPU оказывается эффективным в случае моделирования систем заряженных частиц. Это избавляет от необходимости применять приближенные методы типа PPPM, которые, как показано в предыдущем параграфе, плохо применимы для расчета кластерной наноплазмы.

Малое суммарное число потоков не позволяет задействовать внутренний механизм GPU, скрывающий задержки доступа к памяти за счет быстрого аппаратного переключения между ждущими и готовыми к выполнению потоками. Для оптимизации полученного алгоритма применяет метод расчета, увеличивающий число потоков на частицу для малых N.



Отношение быстродействия GPU Nvidia Tesla M2050 к одному ядру CPU Intel Xeon E5520 при моделировании неидеальной плазмы в зависимости от количества частиц в системе. Крестики - алгоритм с одним потоком на частицу, треугольники - алгоритм с несколькими потоками на частицу.



3. Моделирование неидеальной плазмы, образованной в результате ионизации наноразмерных кластеров

3.1 Постановка задачи

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

При помощи фемтосекундных лазерных установок можно изучать колебания электронов в наноплазме, образующейся в результате быстрой ионизации кластеров короткими лазерными импульсами [6].

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

Для построения теоретической модели указанных выше процессов была написана программа молекулярно-динамического моделирования неидеальной электрон-ионной плазмы [7,11]. В ней были смоделированы колебания электронов в неидеальной наноплазме при неподвижных ионах.



3.2 Численная модель моделирования кластерной наноплазмы

Классический кулоновский потенциал электрон-ионного взаимодействия на малых расстояниях быстро убывает, из-за этого на короткой дистанции при расчётах получаются очень большие по модулю отрицательные значения энергии. Это противоречит квантово-механическим представлениям о наличии минималного уровня энергии (основного состояния) в системе электрона и иона. Для приближенного учета этого явления куловновский потенциал взаимодействия принято заменять псевдопотенциалом (сравнение на рис. 4.2), который при нулевом расстоянии между атомами имеет определённое значение. Использованный в [9,10] потенциал имеет вид:

,

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

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

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



Рис. 4.2. Потенциалы электрон-ионного взаимодействия на коротких расстояниях:

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

3.3 Спектры электронных колебаний в ионизованном нанокластере

Анализ колебаний электронов в электронной подсистеме проводился путем вычисления спектра автокорреляционной функции тока [8]:

,

,

где - полный электронный ток в кластере, - нормировочный объем, - диэлектрическая проницаемость.

Рис. 4.3. Спектры автокорреляционной функции тока для различного числа ионов в кластере Ni (размера кластера R). Пунктирные линии показывают теоретические значения частоты Ми и плазменной частоты .

Теоретическое значение частоты Ми имеет вид . Это значение получено на основе модели, в которой электроны представляют собой сферическое облако, которое колеблется около примерно такого же по размеру облака ионов. Однако в результате удара лазером некоторое количество электронов вылетело, следовательно, кластер имеет положительный заряд, поэтому радиус электронного облака меньше. Это объясняет сдвиг резонанса Ми для маленьких кластеров [5].

Вставка на рисунке 4.3 показывает, что амплитуда колебаний электронов, соответствующих второму максимуму, в разных частях кластера имеет разный знак. Далее будем называть этот пик плазменным резонансом. Теоретическое значение его частоты имеет вид . С увеличением размера кластера амплитуда этого резонанса резко возрастает, и при он становится доминирующим, в то время как резонанс Ми затухает. Как видно из рисунка, частоты первого и второго резонансов стремятся к частоте Ми при увеличении размера кластера, а для частоты третьего резонанса предельным значением является плазменная частота.

Положения основных резонансов в зависимости от размера кластера показаны на рис. 4.4. Определение плазменной частоты становится неточным для кластеров с .

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

молекулярный ионизация наноразмерный кластер



4. Анализ колебаний электронов в кластерной наноплазме

4.1 Теоретическая оценка частоты электронных колебаний

Результаты молекулярно-динамического моделирования неидеальной плазмы потребовали дополнительного теоретического обоснования. Был проведён анализ и проверка вывода теоретических формул, а так же анализ результатов молекулярно-динамического моделирования и дополнение теоретической базы по проведённому эксперименту.

Выведем формулу коллективных колебаний электронов в кластере. Пусть облако электронов имеет фиксированную форму. Допустим, что плотность электронов n(r) и потенциальная энергия электрона в кластере V(r)=-eц(r) сферически симметричны и их центры находятся в равновесии. Предположим сдвиг z относительно минимума потенциала. Необходимая частота в гармоническом приближении:

,

где U - энергия всей системы, а m=Neme - это масса Ne электронов в облаке. В то же время ДU может быть выражено через V(r) и n(r).

,

,

Если потенциал непрерывный, можно продифференцировать его перед интегрированием



,

,

Полагая z=0, с=r, можно получить:

,

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

В работе [8] сделано предположение о том, что распределения электронов и ионов представляют собой однородные сферы разного радиуса (за счет вылета части электронов). На основе этого предположения выведена формула:

,

При помощи программы Wolfram была проверена корректность вывода этой формулы. График частоты, вычисленной по этой формуле, изображён на рис. 3.4 сплошной линией. На нём видно, что на эту кривую не попадают точки, обозначающие результаты молекулярной динамики, следовательно, требуется ввести в неё дополнительные параметры, улучшающие точность.

Для того чтобы аналитические расчёты давали результат, лучше согласованный с результатами молекулярной динамики, было предложено учитывать реальный профиль электронной плотности. Это означает, что электронное облако считается шаром с плотностью, зависящей от радиуса, в то время как ионы остаются замороженными, и их распределение по-прежнему считается однородной сферой. Необходимый для этого профиль электронной плотности брался из результатов МД моделирования путем усреднения средней плотности электронов в различных точках на протяжении достаточно длинной МД траектории, пример его сравнения с постоянной плотностью на рис. 4.1.

Рис. 4.1. Пример сравнения реального профиля электронной плотности с постоянной электронной плотностью. Размер кластера - 3150 частиц.

Расчёты проводились в программе Wolfram Mathematica. Для подстановки в формулу колебаний (5.1) был взят коллективный потенциал кластера: , где б = e, i - тип частиц. После подстановки erf-like потенциала 3.1 получаем (здесь и далее код из программы Mathematica):

,



Для того чтобы программа могла корректно аналитически проинтегрировать потенциал, он был разбит на две составляющие: первая для случая r > r1, вторая для r1 > r:

Затем аналитически было получено выражение fV(r) = и взят его интеграл.

В результате этого получилась следующая функция:

Следовательно, выражение для частоты приняло следующий вид:

,

где n(r) - электронная плотность.

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

Рис. 4.2 Зависимость резонансных частот от числа электронов в кластере. Ромбы - результаты численного интегрирования (4.5) с учётом реальной электронной плотности, полученной из МД моделирования.

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

4.2 Параметризация профиля электронной плотности

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

При изучении графиков можно увидеть закономерность, с которой меняется форма профиля при увеличении числа ионов в кластере. Было выдвинуто предположение, что возможно вывести аналитическую формулу, описывающую зависимость электронной плотности в кластере от радиуса и температуры. Вывод осуществлялся при помощи программы Wolfram Mathematica, которая позволяет наглядно отслеживать зависимость искомой функции от параметров и минимизировать отклонения от референтных данных МД моделирования.

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

,

где fFermi(x,x0,s) - это функции Ферми вида . Здесь и далее расстояния измеряются в нм, концентрация электронов в нм-3, частота в фс-1.

Рис. 4.3 Сравнение выбранной функции с исходными данными



Подбор параметров при помощи опции Manipulate в Wolfram Mathematica. Зеленая линяя - функция, синяя - данные из массива результатов моделирования. N = 55.

Требовалось сократить количество параметров для того, чтобы можно было определять профиль плотности, зная только число частиц. Для этого были подобраны оптимальные наборы параметров такие, что отклонение функции от исходного набора данных было минимальным. Подбор проводился для результатов расчётов в кластерах размером 55 ионов (рис. 4.4), 147 (рис. 4.5), 309 (рис. 4.6), 1000 (рис. 4.7), 3150 (рис. 4.8), 7100 (рис. 4.9).

Рис. 4.5 Сравнение выбранной функции с исходными данными. N = 309.

Рис. 4.6 Сравнение выбранной функции с исходными данными. N = 1000.



Рис. 4.7 Сравнение выбранной функции с исходными данными. N = 3150.

Рис. 4.8 Сравнение выбранной функции с исходными данными. N = 7100.

После подбора оптимальных параметров была установлена их зависимость друг от друга или от размера кластера и температуры. Параметр c1 отвечал за сдвиг вершины последнего пика и являлся, по сути, главным для повторения формы реального профиля, зависимость получилась квадратичной (рис. 4.9). Точка переключения между функциями Ферми отвечала за переход к асимптоте на больших расстояниях, рис. 4.10. Гладкость переключения и частота возмущения для всех размеров кластеров совпали со значениями 0.08 и 8 соответственно. Параметр амплитуды пика отвечал за интенсивность уплотнения на краях кластера (рис. 4.11), а степень сглаживания определяла плавность перехода от относительно равномерного распределения плотности в центре кластера к неравномерной плотности по краям (рис. 4.12).

Затем была построена итоговая функция:

,

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

,

,

,

,

,

,

Эта функция была подставлена в уравнение (4.5), в результате чего получилось окончательное выражение для численного интегрирования:

,

,



Рис. 4.13 Зависимость частоты резонанса от числа ионов.

В результате интегрирования получилась кривая значений частоты, нанесённая на график вместе с предыдущими результатами на рис. 4.14. Эта кривая гораздо лучше описывает результаты МД моделирования.

Рис. 4.14 Зависимость резонансных частот от числа электронов в кластере. Оранжевая линия - результат.



Выводы

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

Проведено сравнение двух методик ускорения МД моделирования систем заряженных частиц: приближение PPPM и полноатомная модель с расчетом сил и решением уравнений движения на графических ускорителях. Показаны явные преимущества использования полноатомной модели.



Список цитированной литературы

1. Gibbon P., Sutmann G. Long-Range Interactions in Many-Particle Simulation. In: Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms eds. J. Grotendorst, et al, Julich: NIC. 2002. V. 10. P. 467-506.

2. Allen M.P., Tildesley D.J. Computer Simulation of Liquids. Oxford : Clarendon Press, 1989.

3. Frenkel D., Smit B. Understanding Molecular Simulation: From Algorithms to Applications. San Diego: Academic Press, 2002.

4. I.V. Morozov, A.M. Kazennov, R.G. Bystryi, G.E. Norman, V.V. Pisarev, V.V. Stegailov. Comput. Phys. Commun. 2011. V. 182. P. 1974.

5. T. Raitza, H. Reinholz, G. Rцpke, I. Morozov, E. Suraud. Laser excited expanding small clusters: Single time distribution functions // Contributions to Plasma Physics. 2009. V. 49. P. 496-506.

6. T. Doeppner, T. Fennel, P. Radcliffe, J. Tiggesbaumker, and K.-H. Meiwes-Broer, Phys. Rev. A 73, 031202 (2006).

7. H. Reinholz, T. Raitza, G. Roepke, and I. V. Morozov, Int. J. Mod. Phys. B 22, 4627 (2008).

8. T. Raitza, G. Roepke, G. Reinholz, and I. V. Morozov, Phys. Rev. E 84, 036406 (2011).

9. G. Zwicknagel and T. Pschiwul, Contrib. Plasma Phys. 43, 393 (2003).

A. Filinov, M. Bonitz, and W. Ebeling, J. Phys. A 36, 5957 (2003).

10. R.G. Bystryi, I.V. Morozov. Electronic oscillations in ionized sodium nanoclusters // J. Phys. B: At. Mol. Opt. Phys. 2015. V. 48. P. 015401.

11. Тимофеев А.В., Морозов И.В., Стегайлов В.В., Норман Г.Э., Куксин А.Ю., Ланкин А.В., Орехов Н.Д., Писарев В.В., Смирнов Г.С., Стариков С.В. Зачем и какие нужны суперкомпьютеры эксафлопсного класса? Предсказательное моделирование свойств и многомасштабных процессов в материаловедении // Программные системы: теория и приложения. 2014. Т. 5. № 1. С. 191-244.»

12. Ewald, P (1921). "Die Berechnung optischer und elektrostatischer Gitterpotentiale". Ann. Phys. 369 (3): 253-287.

13. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J Comp Phys, 117, 1-19 (1995).

14. Winkel M., Gibbon P. Spatially Resolved Electronic Correlations in Nanoclusters // Contrib. Plasma Phys. 2013. V. 53. P. 254-262.



Аннотация

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

В работе предложена теоретическая модель, позволяющая объяснить смещение частоты колебаний Ми в зависимости от размера ионизованного кластера, полученное ранее на основе молекулярно-динамического моделирования и согласующееся с результатами экспериментов. Модель основана на анализе неоднородности профиля электронной плотности. Получена интерполяционная формула для профиля электронной плотности в кластере произвольного размера. Результаты применения этой формулы показывают хорошее совпадение с данными МД моделирования по частоте колебаний Ми. Исследование проводилось при помощи программы Wolfram Mathematica, использованной для численного интегрирования и оптимизации параметров интерполяционной функции.

Объём работы 38 страниц с 23 иллюстрациями, 1 таблицей и 15 источниками. Ключевые слова: наноразмерные кластеры, неидеальная плазма, молекулярная динамика, электронные колебания, колебания Ми.



Irradiation of nanosized metallic clusters by femtosecond laser pulses of moderate intensities (1013 - 1016 W/cm2) leads to creation of the so-called electron-ion nanoplasma with size-dependent properties. Electron oscillations in such plasma are of interest in view of resonance absorption or probing of plasma dynamics in experiments. Results of molecular dynamics simulations are analyzed. The electron oscillation spectra dependence on the cluster size is studied for clusters from 100 to 100000 atoms.

A theoretical model explaining Mie frequency shift according to ionized cluster size is proposed. The model is based on the analysis of electron density profile heterogeneity based on the molecular dynamics simulations results. An interpolation formula for an electron density profile in a cluster of arbitrary size is derived. Results of this formula as applied to frequency calculations are in a good agreement with MD simulation results for the Mie frequency. Wolfram Mathematica package was used for numeric integrations and optimization of interpolation function parameters.

The article consists of 38 pages with 23 illustrations, 1 table and 15 sources.

Keywords: molecular dynamics, nonideal plasma, metal cluster, nanocluster, electron oscillations, Mie resonance.

Размещено на Allbest.ru