Атомистическое моделирование окисления углеродных наноструктур тема автореферата и диссертации по физике, 01.04.17 ВАК РФ



На правах рукописи УДК 539.194

Валуев Илья Александрович

Атомистическое моделирование окисления углеродных

наноструктур

Специальность 01.04.17 — химическая физика, горение и взрыв, физика экстремальных состояний веществ

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

О с с^З 2011

Москва - 2010

4853766

Работа выполнена в Государственном образовательном учреждении высшего профессионального образования "Московский физико-технический институт (национальный исследовательский университет)" на кафедре химической физики (базовый институт — Институт химической физики им. H.H. Семенова РАН) и в Объединенном институте высоких температур РАН.

Научный руководитель: доктор физико-математических наук,

профессор Генри Эдгарович Норман

Официальные оппоненты: доктор химических наук,

профессор Александр Владимирович Немухин

доктор физико-математических наук, профессор Станислав Яковлевич Уманский

Ведущая организация: Институт общей и неорганической химии

им. Н. С. Курнакова РАН

Защита состоится 16 февраля 2011 г. в 14 часов на заседании диссертационного совета Д.002.012.02 при Учреждении Российской академии наук Институте химической физики им. H.H. Семенова РАН по адресу: 119991 г. Москва, ул. Косыгина, д. 4.

С диссертацией можно ознакомиться в библиотеке ИХФ РАН.

Автореферат разослан 15 января 2011 г.

Ученый секретарь

диссертационного совета

кандидат физико-математических наук

М. Г. Голубков

1 Общая характеристика работы

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

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

Отдельные механизмы реакций окисления углерода успешно изучались с помощью метода функционала плотности [3-5]. Моделирование окисления требует расчетов с большим числом частиц, недоступным для первонринцип-ных методов. Атомистическое описание систем С+О на основе полуэмпирических потенциалов взаимодействия затрудняется тем, что реакции с участием кислорода происходят с большим переносом заряда, и требуют учета спина. Создание вычислительно эффективной модели потенциальной энергии системы углерод+кислород, учитывающей разрыв и образование химических связей, является важной задачей.

Развитие вычислительной техники позволяет увеличивать доступные размеры систем для первопринципных численных атомистических моделей, однако, потребность в быстрых/гибридных потенциальных моделях не уменьшается [6]. В работе на примере ограниченного класса реакций демонстрируется построение быстрой потенциальной модели. Развиваемый формальный подход может быть использован как для расширения модели реакций окисления углерода, так и для построения потенциалов для других систем.

Цели работы состоят в: 1) разработке программной среды атомистического моделирования, позволяющей создавать и оптимизировать новые модели взаимодействия; 2) создании универсального подхода для конструирования потенциалов взаимодействия, описывающих разрыв и образование химических связей; 3) систематизации представлений о механизмах окисления углеродного каркаса; 4) разработке быстрого и точного полуэмпирического потенциала, способного описывать реакции окисления.

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

Положения, выносимые на защиту.

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

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

3. Разработан алгоритм кластерной диапазонной интерполяции для построения потенциалов взаимодействия, описывающих разрыв и образование химических связей.

4. Продемонстрирована возможность построения "быстрого" полуэмнири-ческого потенциала взаимодействия атомов углерода и кислорода, воспроизводящего результаты, полученные ab initio.

5. Создан пакет программ атомистического моделирования, предназначенных для постановки распределенных численных экспериментов.

Апробация работы. Основные результаты диссертационной работы были доложены па конференциях "Современная химическая физика" (Туапсе, 2001, 2008, 2009), "Современные проблемы фундаментальных и прикладных наук (МФТИ, 1996-1999), "Уравнения состояния вещества" (и. Эльбрус, 2002, 200С, 2008), "Воздействие интенсивных потоков энергии па вещество" (п. Эльбрус, 2003), "43rd International Field Emission Symposium" (Москва, 1996), "International Conference on Algorithms and Architectures for Parallel Processing" (Австралия, Мельбурн 2005), "International Conference on Computational Science" (КНР, Пекин, 2007), "International Confcrcnce on Computational Science and its Applications" (Малайзия, Куала-Лумнур, 2007), "Conference on Computational Physics" (Италия, Генуя 2004; Бельгия, Брюссель 2007; Тайвань, Гаосюн 2009) и на других российских и международных конференциях.

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

Структура и объем диссертации. Диссертация состоит из введения, семи глав и заключения, изложена на 122 страницах, включая 25 рисунков, 3 таблицы и 109 наименований цитируемой литературы.

2 Содержание работы

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

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

Во второй главе описывается исследование путей реакции взаимодействия углеродных бездефектных наноструктур с различной кривизной поверхности (графен, тубулен 4x4, фуллерен Cqo) с атомарным и молекулярным кислородом с помощью метода функционала плотности. Из результатов расчетов следует, что энергетические характеристики реакций окисления бездефектных углеродных с структур атомарным и молекулярным кислородом зависят от кривизны поверхности структуры. Чем больше кривизна поверхности, тем более экзотермической является реакция и тем меньше ее барьер. Этот вывод согласуется с экспериментальными данными и другими расчетами. Реакция с молекулярным кислородом является экзотермической только для структур с большой кривизной поверхности. Результаты расчетов использовались в дальнейшем при построении полуэмпирического потенциала.

Для расчетов применялась лицензионная версия программы VASP (4.6.31 MPI), базис плоских волн с обрезанием 400 эВ (средняя точность), приближение обобщенных градиентов с функционалом PW91 и псевдопотенциалами PAW. При организации распределенных вычислений программа VASP вызывалась из приложения GridMD, реализующего сценарий сканирования параметров. Атомы и малые молекулы могут присоединяться к углеродному каркасу в разных положениях. Рассматривались наиболее выгодные но энергии конфигурации присоединения, когда атом адсорбируемой системы находится непосредственно над одним из атомов углерода (позиция "top") или над одной из связей углеродного каркаса (позиция "bond"). Для двухатомной молекулы рассматривается также позиция "bridge" (мостик), когда оба атома кислорода находятся над связанными атомами углерода. Для исследуемого места присоединения и заданного спинового состояния строилась зависимость энергии от расстояния адсорбированной молекулы до поверхности углеродного каркаса при фиксированных положениях атомов каркаса и параллельном смещении адсорбируемого атома (молекулы) перпендикулярно поверхности. Смещение производилось с шагом 0,1 Â в диапазоне расстояний до поверхности 0,9-4,5 À. Для каждого из смещений затем проводилась полная оптимизация энергии системы относительно положений атомов (оптимизация геометрии) (Рис. 1).

При наличии барьера, его значение вычислялось методом "подталкивае-

'О, мостик над 6,6 'О, лерп, связи 6,5 'О, мостик над 6,6 поиск барьера в 'О, вдоль пути О

Юг мостик над 6.6 оптимизация "Т-

з

г, А

° 492 -К

5

C,>sO,(GGA)

нанотр.+'О (GGA)

(НН)'0 к атому

.....!02 мосмна npen. связью

'О, поиск барьера

H3HOTp.+J02(GGA)

мостик на лерп. связью - оптимизация —,-1-1-1-1-.

г, А

Рис. 1: Слева: Энергия присоединении молекулярного кислорода к фуллерену Свд и зависимости от (минимального) расстояния от молекулы до поверхности; "мостик" - молекула кислорода ориентирована параллельно связи 0,6; "нерп, связи 6,6:' — молекула кислорода ориентирована перпендикулярно (короткой) связи между шестичлспными циклами. Стрелками отмечены конфигурации локальных минимумов энергии, включая основное состояние реагентов: триплстный кислород) синглетный Соо- Расчет энергий триплетного состояния производился вдоль оптимального пути реакции, найденного для синглетного состояния 'Ог.

Справа: Энергия присоединения атомарного (над атомом Сив центре связи С-С, перпендикулярной оси иапотрубки) и молекулярного кислорода (вдоль связи С-С, перпендикулярной оси напотрубки) к напотрубке 4x4, в зависимости ог (минимального) расстояния до поверхности. Стрелками отмочены конфигурации локальных минимумов энергии, включая основное состояние реагентов: триплетный кнелород+синглетная нанотрубка. Поиск оптимального пути реакции производился для синглетного состояния.

мой упругой ленты" (nudged elastic band), встроенным в программу VASP. При этом оптимизированные состояния реагентов и продуктов соединяются "путем" из 5-10 промежуточных конфигураций. Оптимизация "пути" с дополнительными ограничениями (упругая лента) приводит к нахождению седловой точки с минимальной энергией, которая используется для оценки барьера реакции.

Основное состояние исследуемых углеродных структур с замкнутой оболочкой является синглетным. Поскольку как для атомарного, так и для молекулярного кислорода основное состояние — триплет, в обоих случаях отдельно рассчитывались синглетное и триплетное спиновое состояние всей системы в целом. Для исследуемых систем минимальное по энергии состояние реагентов всегда триплетное, а связанное состояние продуктов реакции — синглетное. Таким образом, реакция с минимальным барьером предполагает изменение спина системы. Такое изменение запрещено для изолированной

Система Выделяемая энергия синглет (триплет), эВ Барьер присоединения сииглет (триплет), эВ Оптимальное присоединение кислорода к углеродному каркасу

Сбо+О 3,9 (2,3) 0(0) над связью 6,6

панотрубка 4x4+0 3,2 (1,6) 0(0) над связью, перпендикулярной оси папотрубки

графен+О 3,0 (1,4) 0(0) над связью С-С

С'боЮг 1,0 (0,06) 0,5±0,1 (1,0) вдоль связи 6,6

панотрубка 4х4+02 0,26 (-0,8) 0,7±0,1 (1,5) вдоль связи, перпендикулярной оси нанотрубки

графси+Ог < -2 > 2,0 вдоль связи С-С

Таблица 1: Энергетические характеристики реакций исследуемых структур с атомарным и молекулярным кислородом.

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

Данные расчетов энергетических параметров реакций исследуемых структур углерода для конфигураций, соответствующих глобальному минимуму, приведены в Таблице 1. Оптимальное присоединение атомарного кислорода к фуллерену достигается при адсорбции кислорода вдоль двойной связи С—С (6,6). Этот вывод согласуется с результатами экспериментов и других расчетов, согласно которым оксид фуллерена имеет "эпоксидную" структуру над короткой связью 6,6. При этом присоединение кислорода приводит к увеличению длины связи на 0,09 А. Состояние глобального минимума энергии — синглетное, энергетический выход относительно триплстного (основного) состояния реагентов составляет 2,3 эВ, относительно возбужденного (синглет-ного) состояния — 3,9 эВ. Присоединение как синглетиого, так и триплстного атома кислорода происходит без барьера. Согласно экспериментальным данным, энергия присоединения атома кислорода к молекуле Соо составляет 3,9 эВ, что согласуется с нашим расчетом для синглетиого состояния. Согласно экспериментальным данным [7], реакция фуллерена с молекулярным кислородом происходит при температурах начиная с 500 К. Наша оценка барьера составляет 0,5 эВ при переходе из синглетиого состояния и 1 эВ при переходе из триплетного состояния.

В третьей главе вводится метод сильной связи и описывается его оригинальная численная реализации. Этот метод в дальнейшем применен в работе для построения быстрого потенциала взаимодействия углерод-кислород. Ме-

тод сильной связи является простым полуэмпирическим квантовым методом для определения электронной структуры и полной энергии связи системы электронов и атомных ядер. Исторически его основы заложены в приближении линейной комбинации атомных орбиталей и работах Слэтера, Костера, Чади, Левина и Харриеопа [8,9]. Простейшей формой метода сильной связи для определения электронной структуры можно считать расширенный метод Хюккеля. В современной формулировке основные определения метода сильной связи выглядят следующим образом.

При фиксированных положениях атомных центров (ионов) одночастич-ный гамильтониан Н системы валентных электронов записывается в базисе "квазиатомных" орбиталей ■фыа, локализованных на ионах: <р = ^^таФьпа-Собственные значения еп получаются как решения сскулярного уравнения, а электронная часть энергии основного состояния системы Е^ — как сумма энергий заполненных электронных уровней:

(Я-еп5)(7<'г> = 0 (1)

¿¡М = Т)е„. (2)

Здесь С^ — вектор неизвестных коэффициентов, Н и 5 — матрицы гамильтониана и перекрытия орбиталей соответственно, дп{ц, Т) — кратности заполнения одноэлектронных уровней, зависящие в общем случае от химического потенциала системы /х и температуры Т. Орбитали фша, как и орбитали изолированных атомов, обладают симметрией, позволяющей классифицировать их по собственным состояниям орбитального момента (индексам I и га). Здесь и в дальнейшем греческие индексы (например, а) нумеруют ионы. Постулируется предположение, что Н может с хорошей точностью описываться небольшим набором орбиталей.

В отличие от матрицы Фока, используемой в методах МО ЛКАО для решения секулярного уравнения, гамильтониан Н не содержит многоцентровых интегралов кулоновского отталкивания между электронами вида

(аЬ\с<1) = У У ^(1)^(1)^(2)^(2)^^2, (3)

поскольку предполагается, что их усредненное значение компенсируется притяжением ядер: и ^ ZaZ|>R~^■

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

Базисные орбитали называются "квазиатомными", поскольку они отличаются от орбиталей изолированных атомов. В идеологии метода сильной связи их конкретная пространственная форма (кроме угловой симметрии) не выясняется, а элементы гамильтониана и матрицы перекрытия 5 записываются параметрически как функции расстояния гаа' = |г<* — га<| между двумя ионами: ií/ma/-mv = <Sw'mV - StmUm>(raa>). Здесь сделано до-

полнительное предположение о парном характере зависимости матричных элементов одночастичного гамильтониана от положений ионов.

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

П2

Hlmal'm'a' = ~^Ulml2rISK{^aa')} {Гаа') ^ Himnima = -e¡a S = /, (4)

VsstT = -1,4 Vspa = 1,84 Vvpa = 3,24 = -0,81 /(r) = i. (5)

Здесь ?75/í(naa') — универсальные угловые коэффициенты перекрытия Слэтера-Костера, зависящие от ориентации вектора, соединяющего центры атомов naa' — raa'/raa'i относительно базиса, в котором измеряется угловой момент орбиталей; e¡a — энергии одноэлектронных уровней изолированных атомов в приближении Хартри-Фока, Уц>тп — универсальные коэффициенты Харрисона. Как мы видим, парные матричные элементы определяются универсальными зависимостями, а от конкретных химических элементов зависят только диагональные члены гамильтониана. В случае, когда требуется знание не только электронной структуры, но и полной энергии связи Етв, к электронной энергии добавляется отталкивательный член, параметризующий одновременно кулоновское взаимодействие между ионами и часть обменно-корреляционной энергии не учтенной (учтенной дважды) в электронном члене:

Етв = £o,ei + l/«p(R). (6)

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

Urep{R) = UaQ>(r), Uaa'(r) = АаЫ ехр(^- + Саа'Г), (7)

ct<a'

Рис. 2: Экспериментальная потенциальная кривая молекулы Сг, аппрокпмированная потенциалом Морзе (Слева) и полученная методом сильной связи в параметризации Харрисопа с отталкиванием в форме (7) (Сгцша)

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

Разработана универсальная численная реализация метода сильной связи, которая допускает ортогональную и неортогопальную параметризацию. Рассмотрены вычислительные процедуры, необходимые для определения собственных значений и векторов. В параметризации Харрисопа с отталкиванием в форме (7) автором диссертации произведен ряд расчетов свойств углеродных наноструктур, краткое описание которых приводится в третьей главе. Эти расчеты можно рассматривать как апробацию программной реализации метода сильной связи. Проанализирована структура электронных оболочек фуллерена, распадающаяся на р— и 7Г— подсистемы. По энергии связи и величине щели между верхним заполненным и нижним вакантным состояниями и величине полной энергии предсказаны ряды стабильности фуллеренов Ою, Сто, С78 и их изомеров. Рассчитаны энергия связи и геометрии интеркалиро-вапных фуллеренов и нанотрубок.

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

ниане Харрисона ННагг), предложенной Меноном и Саббасвами [10]:

Н/тЫ'т'а' = + К 1 — Б^) Н1та[та = —е;а (8)

<5 = Зсю>Ц1ТП12Г18к{Паа')/('Гаа')> а'11!т,п = ~ГГ,-Г7—-Г (9)

К{Гаа>)(е1а +

52 = (Бж - - 35ии)/4 (10)

К(г) = К0 ехр(С(г - с?0)2) /(г) = /о ехр(-Л(г - (11)

= Щ ехр(~В(гаа, - ¿о)) В = 44. (12)

Здесь универсальным способом учтено перекрытие орбиталей. Основными изменяемыми параметрами модели служат равновесные расстояния ¿о и константа отталкивания Щ. Постулируется, что остальные параметры модели (Ко, А, С, /о) можно считать универсальными (не зависящими от конкретных химических элементов). Данная параметризация была успешно использована группой М.Мепона при моделировании кремниевых кластеров и оксидов фуллерепов и напотрубок.

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

Отправной точкой для КДИ является состав системы (множество атомов и множество связей, которые могут быть между ними сформированы). Мы будем характеризовать "прочность" каждой связи некоторым значением р в интервале от 0 до 1. Это значение может трактоваться как вероятность формирования данной связи (связывающее состояние). Вероятность того, что данная связь не образовалась (несвязывающее состояние) выражается как 1 — р. Функция р(г) зависит от расстояния между атомами, образующими связь (и в общем случае от наличия других связей). Характер межатомного взаимодействия может быть включен в эту зависимость. Например, можно строить р как функцию характерного матричного элемента гамильтониана между двумя атомами или на основе порядка связи.

Относительное растяжение с: Ц]аг МД

Рис. 3: Тестовые расчеты КДП с потенциалом ил,(11) = — 0.5Л/(Л/ — 1) Сякьа: Зависимость энергии и весовых коэффициентов P¡ от относительного изменения длины при равномерном расширении линейной молекулы: штриховая линия - учет 1 конфигурации, пунктирная линия - учет 3 конфигураций, сплошная лшшя - все конфигурации.

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

Рассмотрим систему и:? N частиц, для которой известны гладкие потенциалы Í7'V(R), описывающие взаимодействие внутри связанной подсистемы (кластера) из А/ < N частиц, а также задано взаимодействие па больших расстояниях U'NB(r) при отсутствии химических связей (например, силы Ван-дер-Вальса). Под кластером мы будем понимать множество частиц, в котором для каждой частицы существует другая из того же множества, находящаяся от первой на расстоянии не более заданного характерного расстояния образования химической связи rr¡. Энергия U(R) для этой системы, являющаяся гладкой интерполяцией между UM{R) с разными А/ и (7хв(г), согласно КДИ записывается следующим образом:

U( R) = £ Ps(R)Us(R), £ Ps( R) = 1, (13)

5 S

E и™ш, (14)

где S индексирует кластерные конфигурации по типу связности, Cs индекси-

руст все кластеры в составе данной конфигурации 5, Лгр5 — это число атомов в кластере Св, IIм(В) — потенциалы химической связи для системы из М атомов, и™(Япъ) — парно-аддитивные несвязывающие потенциалы для атомов а и Ь, II — множество координат всех атомов, Наь — расстояние между атомами а и 6, /^(Л) — гладкие весовые функции. Переключающий множитель Рз(К) (вероятность данной кластерной конфигурации) определяется согласно КДИ как произведение элементарных множителей для отдельных связей:

где = (П.* — — расстояние между атомами г и j. т%(Я) = р{В) если связь между атомами г и у существует в данной конфигурации и и>у(Д) = 1 — р(Н), если связи не т. Функция р(Я) (элементарный переключатель) — гладкая монотонно убывающая от 0 до 1. Для построения смеси кластерных конфигураций, используемой о (13) для вычисления энергии, необходимо сгенерировать набор конфигураций связей, в которых некоторые связи образовались, а некоторые нет, при этом каждая конфигурация должна иметь свой вычисляемый вес в ансамбле. Прямой путь — это вариация по всем 2В конфигурациям, которые могут быть сгенерированы путем перебора всех свя-зывающих/несвязываюгцих состояний В связей. Генерация всех возможных конфигураций может быть осуществлена на практике только если число связей в системе невелико. Построены алгоритмы, позволяющие учесть только значимые конфигурации с контролируемой точностью. На Рис. 3 иллюстрируется применение КДИ для простейшей модели связи М частиц, в которой энергия зависит только от числа частиц в кластере,

В пятой главе с целью демонстрации возможностей метода кратко описана апробация метода КДИ, который применялся для молекулярцо-динамического (МД) моделирования водородной плазмы. В такой модели электроны являются одним из компонентов в "химической реакции", например, динамически моделируются процессы вида.:

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

(15)

Н ^ Н+ + е

р(£й)

н,

н

(а)

I i

' I

I

i

Рис. 4: Энергетический спектр р(Е\) (а) и парные корреляционные функции (Ь) - (<1) для равновесной водородной плазмы при Г, = 1 и различных потенциалах взаимодействия, Т — 10000 К. Сплошная линия - потенциал КДИ, штриховая линия — кулоновский потенциал с обрезанием — 13,5 эВ, штрих пунктирная линия — кулоновский потенциал с обрезанием —ХГ

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

Для описания взаимодействий на малых расстояниях использовались следующие элементарные модели. Взаимодействие между двумя протонами, окруженными одним, двумя или тремя электронами (водящими с этими электронами в один кластер) моделировалось потенциалами Морзе, подобранными для описания равновесных квантовых энергий связи и расстояний в молекулах (молекулярных ионах) Нг, and Hj соответственно. Для целей данной модели энергия взаимодействия в кластерах, содержащих более двух протонов или более трех электронов полагалась настолько большой (+10 эВ на

частицу), чтобы запретить образование подобных связанных состояний в системе (область металлоподобиого водорода исключалась из рассмотрения). Потенциальная поверхность для 1-3 электронов, связанных вблизи 1-2 протонов полагалась плоской на уровне квантовой энергии связи кластера (взаимодействие между электронами "внутри" атомов или молекул полагалось равным нулю). Анализ спектров энергии и парных корреляционных функций, полученных усреднением по равновесной МД-траектории при заданной температуре Т позволяет сделать вывод о наличии в системе связанных состояний, соответствующих атомам и молекулам, которые не воспроизводятся простыми псевдопотенциальнымн моделями (Рис. 4).

В шестой главе описывается построение полуэмпирического потенциала взаимодействия атомарного и молекулярного кислорода с бездефектным углеродным каркасом на основе метода сильной связи. Обсуждаются ограничения метода сильной связи в формулировке (1-2), который формально не подходит для описания атомных подсистем, содержащих разные химические элемента и не взаимодействующих друг с другом. Модель сильной связи интересна как основа полуэмпирпческого универсального потенциала для систем углерод+кислород, однако требует некоторых доработок и поправок. Требуется устранение или компенсация следующих систематических или локальных (зависящих от конкретного расположения атомов) ошибок метода: 1) проблема перетекания заряда и зависимости энергии связи одной удаленной подсистемы от наличия другой (систематическая ошибка); 2) ошибка заполнения несвязанных орбиталей в связанном состоянии подсистемы (локальная ошибка); 3) выбор нуля взаимодействия для сравнения конфигураций с разной координацией атомов (систематическая ошибка); 4) неверный учет взаимного влияния электронов друг на друга (локальная ошибка);5) отсутствие описания для мультиплетиых состояний (систематическая ошибка). При этом требуется сохранить быстродействие, метода и возможность аналитически вычислять производные энергии по положениям ядер атомов.

Отметим, что некоторые из перечисленных ошибок могут быть устранены явным введением самосогласования по кулоновскому и обменному взаимодействию в схему сильной связи [11]. Это приводит к существенному снижению быстродействия метода и не рассматривалось в рамках данной работы. В ней за основу борется несамосогласованный вариант метода сильной связи и демонстрируется возможность компенсации ошибок с помощью локальных поправок и КДИ.

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

модифицированной электронной части энергии сильной связи произвольной системы атомов С и О следующего вида:

Сумма в выражении (17) берегся по всем значимым конфигурациям определенной связности. Переключающий множитель рс>ши.г(К-) (вероятность данной кластерной конфигурации) имеет форму (15). Здесь Етв(Пс), Ет в (&<.)) - энергии сильной связи углеродной и кислородной подсистем, рассчитанные без учета взаимного влияния, [¿?гд(К-с)]ксС=эо и [£тв(1*о)]|11оо==о ~ пУле" вьте энергии отдельных атомов кислородной и углеродной подсистем соответственно. В силу одинаковости одноэлектронных атомных уровней внутри атомов подсистемы, в каждой из подсистем проблема перетекания заряда не возникает. Величина Еша-ась^с, Ко) вычисляется с помощью (18) как разность энергии сильной связи, рассчитанной обычном способом и энергии, полученной в предположении, что все парные расстояния углерод-кислород равны бесконечности.

В выражении (18) допускается перетекание заряда, но при разведении на бесконечность углеродной и кислородной подсистем, величина Дпитась ста~ нет равной нулю и неправильное распределение электронов не повлияет па энергии углеродной и кислородной подсистем. Таким образом, мы выбрали специальную процедуру определения нуля энергии для системы со связанными подсистемами из атомов кислорода и углерода. Чтобы определить конфигурацию атомов с нулевой энергией из произвольной связанной конфигурации, нужно сначала развести на бесконечность подсистемы С и О, вычислить энергию Е-ти,гпЛ, а затем стандартным образом развести все атомы на бесконечность в подсистемах С и О. Как и в ^модифицированном методе сильной связи, возможно быстрое аналитическое дифференцирование Е$ер(И) по К с применением теоремы Гельмана- Фейнмана для квантовых составляющих Етв• Получающиеся силы являются непрерывными функциями координат и могут использоваться в методе молекулярной динамики.

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

Е,ер(К) = £ЫКс)-[ёЫКс)К=оо +

+ Етв{Шо) - [£'Гл(К-о)]коо=оо +

(1С)

(17)

(18)

ЕггйегшлО^-) = Етв(&) ~ [■Етв(Нс) + -Етв(Ко)] кс<

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

CWR) = Al{l-Pco^)){l-Poo{Roo)) +

+Poo{Roo)J2f(R"o,Rco)/2- (19)

с

Здесь рсо(х), Роо(х) — стандартные функции-переключатели КДИ, в частности, роо соответствует вероятности связи между двумя атомами кислорода. Первый член сумме (19) соответствует конфигурациям с одним атомом кислорода и компенсирует недооценку энергии для двойной связи С—О. Здесь для всех расстояний Rco Для данного атома кислорода выбираются два наименьших R\ и i?2- Если различие между ними Д = — i?i)/T?2 больше 10%, и атом кислорода не является атомом молекулы кислорода, то применяется дополнительная поправка величиной в А\ = 1 эВ. Кроме того, для атомов кислорода, связанных в молекулу, применяется поправка, парно-аддитивная по связям С-0 (второй член в сумме (19). Функция f(x,y) представляется в форме двумерного сплайна, коэффициенты которого находятся минимизацией невязки полуэмпирической модели относительно первопринципной ППЭ для системы Соо+Ог. Такая ППЭ была представлена в двумерном виде как функция расстояния от молекулы фуллерепа до центра связи в молекуле кислорода и длины связи в молекуле кислорода. Результирующая максимальная невязка составляет до 0.05 эВ на атом кислорода. Эта же двумерная функция используется для описания различия между триплетным и синглетным листом ППЭ (для разных листов получены разные сплайновые коэффициенты). Проверка построенного потенциала на исследованных углеродных наноструктурах и при вычислении барьеров показала, что он с хорошей точностью (до 0.1 эВ) воспроизводит результаты, приведенные в Таблице 1.

В седьмой главе кратко описываются принципы построения оригинальной программной среды GridMD, с помощью которой в рамках диссертационной работы проводилось численное моделирование. GridMD — это системно-независимая компонентная библиотека языка C-I-+, предназначенная для быстрого конструирования распределенных приложений для атомистических расчетов. Аббревиатура GridMD сформирована из английских слов Grid-enabled Molecular Dynamics (молекулярная динамика для распределенных вычислительных сетей).

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

Рис. 5: Базовые элементы графа исполнения: а - узлы со связью по вычислениям, 6 - цепочка узлов с генерацией промежуточных данных и контрольной точкой. Жесткие логические связи показаны тонкими сплошными линиями, связи по данным - штриховыми линиями, связи по вычислениям - сплошными жирными линиями; различные типы сценариев и их графы исполнения: с - ветвление; й - ветвление с контрольной точкой; е - частично определенный граф общего вида.

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

Новизна подхода Сп<1МО заключается в предоставлении пользовательких функций второго уровня. В отличие от традиционного подхода, где сценарий численного эксперимента реализуется вручную либо с помощью управляющего системно-зависимого скрипта, сценарий СпсШБ организуется системно-независимым способом средствами языка С++. Сценарий организован в виде графа исполнения, состоящего из узлов (элементов программы), допускающих связи по управлению либо связи по данным (Рис. 5). Связи по данным являются источником параллелизма, так как допускают производство данных на физически разных вычислительных системах с последующей их передачей по сети. Постановка задач на удаленной вычислительной системе и передача данных производятся стандартным способом заранее разработанными специализированными компонентами. Детали работы сервисных компонент скрыты от пользователя библиотеки - специалиста-физика, программирующего численный эксперимент. В диссертации приводится описание алгоритма

int griduid_main(int argc.char* argv[]) {

gmExpeiimeiit.init(argc,argv); // initialization begin_distributed();

// Define the skeleton and internal data link types

gmFork<void,val_t,void> forkl("loop");

foikl.begin.hereO; 11 marks the loop 'begin' node

int nterms =3; // number of terms in the series

val.t sum =0.; // 'sua' accumulates the result

for(int i=0; i<nterms; i++){

// Create a new 'split' node and define its output

if (f orkl. split O)

forkl.vsplit_out() = pow(i, (vai_t)i) / i;

// Define the action of the 'merge' node if(forkl.merge()) sum += forkl.vmerge_in(); // accumulation of the terms

>

if(end_distributed()) // loop 'end' node is optional

printf("The result is V,g\n", sum); return 0;

loop.split(0)

]

kK>p.split(l) ]

loop.splii(2)

1 l_0.oul/2_0.in '3_0.oul/4_0.in 15_0,oul/6_0.in

J_

jL

loop.mcrgeO)

lnop.mcrgc(2)

loop.mcrgeO)

Рис. 6: Пример использования распределенного алгоритмического шаблона ветвления gmFork в Сп<ШБ. Слева: Программа с вызовами функций СпсЛГО, вычисляющая ряд 1п(г). Справа: Граф исполнения задачи, автоматически сгенерированный для исходной программы. Жесткие логические связи показаны тонкими сплошными линиями, связи по данным - штриховыми линиями. В данном примере СпсЗЛГО сформирует и выполнит 3 независимых задания из узлов (0,1), (0,3), (0,5) и завершающее задание из узлов (2,4,6,7)

разбора графа исполнения для определения независимых стадий сценария. GridMD является свободно распространяемой программой и предоставлена в открытый доступ на сайте http : //gridmd. sourcef orge. net/.

Большинство масштабных молекулярно-динамических расчетов обладают характеристиками, позволяющими разбить их на ряд независимых стадий для выполнения на распределенных вычислительных ресурсах. Так, при определении параметров полуэмпирического потенциала взаимодействия по данным первопринципных расчетов потребовалось вычисление порядка 103 независимых точек на первопринципной поверхности потенциальной энергии. Для стандартных распределенных сценариев в GridMD разработан ряд скелетных конструкций, например ветвление (Рис. 6). Также с помощью GridMD был реализован распределенный вариант метода "подталкиваемой упругой ленты" для нахождения потенциальных барьеров реакций присоединения молекулярного кислорода к углеродному каркасу.

3 Достоверность результатов

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

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

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

Достоверность результатов вычисления поверхности потенциальной энергии системы кислород [ углерод методом функционала плотности проверялась сравнением с экспериментальными данными по энергиям окисления и геометрической конфигурации известных оксидов фуллерена и нанотру-бок. Установлено совпадение энергий присоединения атомарного кислорода к фуллеренам и нанотрубкам (2-4 эВ) с экспериментальными оценками и совпадение геометрии оксида фуллерспа, имеющего эпоксидную структуру [12,13]. Результат о безбарьерном присоединении как триплетного, так и синглетиого атомарного кислорода к бездефектному углеродному каркасу совпадает с экспериментальными фактами. Результат о наличии барьера не менее 1 эВ при присоединении молекулярного триплетного кислорода к поверхности фуллерена качественно совпадает с экспериментальной оценкой [7]. Уменьшение барьера до 0,5 эВ для синглстпого кислорода согласуется с экспериментальным фактом об уменьшении энергии активации реакции с углеродным каркасом при возбуждении кислорода излучением. Наличие барьера не менее 3 эВ для присоединения триплетного кислорода к базовой поверхности графита согласуется с экспериментальным фактом о пассивности базовой поверхности и о роли дефектов и краев в процессе окисления [14].

Проводило«) сравнение полученных результатов с другими расчетами методом функционала плотности [3,5,11]. Результаты данной работы на уровне качественных выводов совпадают с результатами расчетов другими авторами. Различия до 0,5 эВ в оценке энергетических параметров реакций обусловлено применением разных вариантов метода функционала плотности.

Достоверность результатов представления поверхности потенциальной энергии системы кислород-Ьуглерод с помощью построенного полуэмпирического потенциала проверялась сравнением энергии различных систем вычисленных с помощью построенного потенциала и напрямую методом функционала плотности. Сравнение проводилось для конфигураций, не входящих в набор для подгонки параметров. Сравнение показывает хорошую точность представления поверхности потенциальной энергии (до 0,05 эВ на атом кислорода вблизи равновесия и до 0,1 эВ на атом кислорода при отклонении от равновесия до 5 эВ).

Основные результаты и выводы работы

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

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

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

3. Для исследуемых систем (графена, нанотрубки 4x4 и фуллерена Cßo) минимальное по энергии состояние реагентов всегда триплетное, а связанное состояние продуктов реакции —■ синглетное.

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

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

6. Продемонстрирована возможность построения "быстрого" полуэмпирического потенциала взаимодействия атомов углерода и кислорода на основе МСС в параметризации Менона-Саббасвами и КДИ.

7. Создана библиотека GridMD для постановки распределенных численных экспериментов (http ://gridmd. sourceforge.net/).

Публикации автора по теме диссертации

В реферируемых журналах

1. Валуев И. А., Норман Г. Э., Шуб Б. Р. Механизмы окисления бездефектных поверхностей углеродных наноструктур: влияние кривизны поверхности // Химическая физика,— 2011. — Т. 30, № 1.-С. 82.

2. Morozov I., Valuev I. Distributed applications from scratch: Using gridmd workflow patterns // Lecture Notes in Computer Science (LNCS). - 2007. - Vol. 4489. P. 190.

3. Valuev I. Simulation of hydrogen plasma with cluster multi-range interpolation // J. Phys. A: Math. Gen. - 2006,- Vol. 39.- P. 4465.

4. Valuev I. A. Reactive potentials for molecular dynamics with cluster multi-range interpolation // Computer Physics Communications. — 2005. — Vol. 1G9, no. 13. — P. CO.

5. Valuev I. A. Gridmd: Program architecture for distributed molecular simulation // Lecture Notes in Computer Science (LNCS). - 2005. - Vol. 3719. - P. 309.

ß. Kuksin A. K, Morozov I. V., Norman G. E., Stegailov V. V., Valuev 1. A. Standards for molecular dynamics modelling and simulation of relaxation // Molecular Simulation. — 2005. — Vol. 31,- P. 1005.

7. Валуев И. А., П.H.Дьячков, Каклюгип А. С., Норман Г. Э. Квантовохимическое моделирование стросция, электронной структуры и энергии образования интеркалированного литием и чистого фуллерена // Ж. неорг. хим. — 2000. — Т. 44, № 3. — С. 472.

8. Валуев И. А., Каклюгип А. С., Норман Г. Э. Численное моделирование распределения электронной плотности поверхностных электронных состояний // Поверхность. Рентгеновские, сипхротропные и нейтронные исслсд. — 1999. — № 7. — С. 56.

9. Валуев И. А., Дьячков П. Н., Каклюгин А. С., Норман Г. Э. Квантово-химичеекая трактовка капиллярных явлений при кптеркалироваиии нанотрубок // Известия РАН. Серия физическая. - 1998. - Т. 62, № 6. - С. 116.

10. Валуев II. А., Каклюгип А. С., Норман Г. Э., Погудим А. В. Электронная структура и энергии связи фуллеренов Сбо, С70, С76 и С78 // Химическая физика. — 1997,— Т. 16, .Ys 6.-С. 88.

11. Валуев И. А., Каклюгин А. С., Норман Г. Э. Структура одпоэлскгрошшх уровней л оболочек фуллерена Сел // Химическая физика. — 1996. — Т. 15, Л4 5. — С. 26.

12. Валуев И. Л., Каклюгин А. G., Норман Г. Э. Электронная структура и химические свойства фуллерена Сео // Известия РАН. Серия физическая. — 1996. — Т. 60, JY« 5. — С. 24.

13. Валуев И. А., Каклюгин А. С., Норман Г. Э. Энергии и электронные волновые функции больших молекул с ковалентными связями // Химическая физика. — 1996. — Т. 15, № 95. — С. 53.

В сборниках

1. Валуев И. А., Норман Р. Э., Шуб Б. Р. Полуэмпирическая модель взаимодействия О и О2 споверхностью углеродных наноструктур // Тезисы XXI симпозиума Современная химическая физика. — 2009. - С. 7.

2. Valxiev I. A. Combined Multi-Range Interpolated Potential for MD Simulations // Computational Physics / Ed. by X. G. Zhao, S. Jiang, X. J. Yu.-- Rinton Press Inc., USA, 2005,-P. 154.

3. Морозов И. В., Валуев И. А. Пакет программ молекулярно-динамического моделирования с интерфейсом для коллективного доступа // Научный сервис в сети Интернет: технологии распределенных вычислений. — Москва: МГУ, 2005. — С. 134.

4. Богданов А. В., Валуев И. А., Городничев М. А., Евлампиев А. А., Корхов В. В., Лузин 77. И. Опыт создания экспериментального сегмента Grid для тяжелых инженерных приложений // Интернет: технологии распределенных вычислений. — Москва: МГУ, 2005. — С. 33.

Список литературы

[1] Haynes В. S. A turnover model for carbon reactivity i. development // Combustion and Flame. — 2001. - Vol. 12G. - Pp. 1421-1432.

[2] Дьячков П. Я. Углеродные нанотрубки. — М.: Бином Лаборатория знаний, 2006.

[3| Хи Y.-J., Ы J.-Q. The interaction of molecular oxygen with active sites on graphite: a theoretical study // Chem. Phys. Lett. - 2004. - Vol. 400. - Pp. 406-412.

[4] Montoya A., Mondragon F., Truong T. N. First-principles kinetics of CO desorption from oxygen species on carbonaceous surface // J. Phys. Chem. B. — 2002. — Vol. 106, no. 16. — Pp. 4236-4239.

[5] Proudakis G. E., Schnell M., Mithlhiiuser M., Peyerimhoff S. D., Andriotis A. N., Menon M., Sheetz R. M. Pathways for oxygen adsorption on single-wall carbon nanotubes // Phys. Rev. B. — 2003. — Vol. 68. — Pp. 115435-1-5.

[6] Nemukhin A. V., Grigorenko B. L., Topol 1. A., Burt S. K. Flexible effective fragment QM/MM method: Validation through the challenging tests //J. Comput. Chem. — 2003. — Vol. 24, no. 12. - Pp. 1410-1420.

[7] Елецкий А. В., Смирнов Б. M. Фуллерены и структуры углерода // УФН. — 1995. — Т. 165, № 9. - С. 977-1009.

[8] Левин А. А. Введение в квантовую химию твердого тела. — М.: Химия, 1974.

[9] Harrison W. Electronic Strucure and Properties of Solids. — W.H. Freeman and Co., San Francisco, 1980. - Vol. 1.

[10] Menon M., Subbaswamy K. R. Optimized structures of СбоО and C60O2 calculated by a damped molecular dynamics optimization scheme // Chem. Phys. Lett. — 1993.— Vol. 201, no. 13.— Pp. 321-325.

[11] Zhu X. Y., Lee S. M., Lee Y. H., Frauenheim T. Adsorption and desorption of an O2 molecule on carbon nanotubes // Phys. Rev. Lett. - 2000. — Vol. 85, no. 13. - Pp. 1003-1008.

[12] Heymann D., Baehilo S. M., Weisman R. В., Cataldo F., Fokkens R. H., et al. СеоОз, afullerene ozonide: synthesis and dissociation to СбоО and O2 // J. Am. Chem. Soc. — 2000. — Vol. 24. — Pp. 11473—11479.

[13] К. M. Creegan, J. L. Robbins, W. K. Robbins, J. M. Millar, R. D. Sherwood, P. J. Tindall, D. M. Cox, et al // J. Am. Chem. Soc. - 1992. - Vol. 114. - P. 1103.

[14] Lee S. M., Lee Y. H., Hwang Y. G., Hahn J. R., Kang H. Defect-induced oxidation of graphite // Phys. Rev. Lett. - 1998. - Vol. 82, no. 1. - Pp. 217 -220.

Подписано в печать 08 ноября 2010 г. Объем 1,2 п.л. Тираж 100 экз. Заказ № 991 Отпечатано в Центре оперативной полиграфии ООО «Ол Би Принт» Москва, Ленинский пр-т, д.37

Введение

1 Обзор литературы

1.1 Современные методы атомистического моделирования

1.2 Исследование окисления углерода.

1.3 Методы программирования и организации сложных вычислительных экспериментов

2 Построение ППЭ для систем углерод-кислород методом функционала плотности

2.1 Введение.

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

2.3 Модель

2.4 Поверхность потенциальной энергии.

2.5 Обсуждение.

2.6 Выводы.

3 Метод сильной связи

3.1 Введение.

3.2 Матричные элементы Левина-Харрисона.

3.3 Моделирование наноструктур углерода методом сильной связи

3.4 Выводы и обсуждение.

4 Интерполяционный потенциал на основе топологии связей

4.1 Связанный и несвязанный предельные случаи.

4.2 Кластеризация и смешивание состояний.

4.3 Построение весовых функций.

4.4 Тестовые примеры.

5 Применение схемы КДИ для описания связанных состояний в неидеальной плазме

5.1 Введение.

5.2 Интерполяционная схема.

5.3 Модель водородной плазмы с учетом ионизации и связанных состояний.

5.4 Молекулярно-динамическое моделирование.

5.5 Выводы.

6 Построение полуэмпирического потенциала для описания реакций наноструктур углерода с кислородом

6.1 Обоснование выбора модели сильной связи в качестве основы полуэмпирического потенциала.

6.2 Ограничения модели сильной связи.

6.3 Калибровка метода сильной связи по полной энергии

6.4 Решение проблемы неверного переноса заряда.

6.5 Дополнительные поправки к электронной энергии.

6.6 Тестирование потенциала.

7 Разработка программной среды для атомистического моделирования

7.1 Введение.

7.2 Концепция библиотеки СпсИУГО и используемые технологии

7.3 Краткое описание архитектуры.

7.4 Пример программы, использующей Спс1МО

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

Углерод характеризуется большим разнообразием форм, в которых он представлен в природе. Помимо кристаллических решеток графита (наиболее энергетически стабильная структура) и алмаза, существует громадное число структур различного строения с близкими к графиту по энергии связи в расчете на один атом. Эти структуры, как правило, обладают поверхностью, в которую выстраиваются пяти- или шестичленные циклы из связанных атомов углерода. К таким структурам относятся фуллерены, нанотрубки, графеновые листы. Пространственные масштабы, характеризующие такие образования - от нескольких ангстрем до сотен нанометров, поэтому мы пользуемся общепринятым термином «наноструктура» для их обобщенного названия.

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

Отдельные механизмы реакций окисления углерода успешно изучались с помощью метода функционала плотности [3-5]. Моделирование окисления требует расчетов с большим числом частиц, недоступным для первопринципных методов. Атомистическое описание систем С+О на основе полуэмпирических потенциалов взаимодействия затрудняется тем, что реакции с участием кислорода происходят с большим переносом заряда, и требуют учета спина. Создание вычислительно эффективной модели потенциальной энергии системы углерод+кислород, учитывающей разрыв и образование химических связей, является важной задачей.

Развитие вычислительной техники позволяет увеличивать доступные размеры систем для первопринципных численных атомистических моделей, однако, потребность в быстрых/гибридных потенциальных моделях не уменьшается [6]. В работе на примере ограниченного класса реакций демонстрируется построение быстрой потенциальной модели. Развиваемый формальный подход может быть использован как для расширения модели реакций окисления углерода, так и для построения потенциалов для других систем.

Цели работы состоят в 1) разработке программной среды атомистического моделирования, позволяющей создавать и оптимизировать новые модели взаимодействия; 2) создании универсального подхода для конструирования потенциалов взаимодействия, описывающих разрыв и образование химических связей; 3) систематизации представлений о механизмах окисления углеродного каркаса; 4) разработке быстрого и точного полуэмпирического потенциала, способного описывать реакции окисления.

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

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

Разработан и реализован программно оригинальный метод кластерной диапазонной интерполяции (КДИ) для конструирования потенциалов взаимодействия.

На основе метода сильной связи в параметризации Менона-Саббасвами построена полуэмлирическая модель поверхности потенциальной энергии для системы углерод+кислород.

Разработан оригинальный молекулярно-динамический пакет ОпсШБ, использующий новые архитектурные принципы, в частности, автоматическую генерацию распределенных сценариев исполнения.

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

Положения, выносимые на защиту.

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

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

3. Разработан алгоритм кластерной диапазонной интерполяции для построения потенциалов взаимодействия, описывающих разрыв и образование химических связей.

4. Продемонстрирована возможность построения "быстрого" полуэмпирического потенциала взаимодействия атомов углерода и кислорода, воспроизводящего результаты, полученные ab initio.

5. Создан пакет программ атомистического моделирования, предназначенных для постановки распределенных численных экспериментов.

Апробация работы. Основные результаты диссертационной работы были доложены на конференциях "Современная химическая физика" (Туапсе, 2001, 2008, 2009), "Современные проблемы фундаментальных и прикладных наук (МФТИ, 1996-1999), 'Уравнения состояния вещества" (п. Эльбрус, 2002, 2006, 2008), "Воздействие интенсивных потоков энергии на вещество" (п. Эльбрус, 2003), "43rd International Field Emission Symposium" (Москва, 1996), "International Conference on Algorithms and Architectures for Parallel Processing" (Австралия, Мельбурн 2005), "International Conference on Computational Science" (КНР, Пекин, 2007), "Conference on Computational Physics" (Италия, Генуя 2004; Бельгия, Брюссель 2007; Тайвань, Гаосюн 2009) и на других российских и международных конференциях.

Публикации. По материалам диссертации опубликовано 13 работ в реферируемых научных изданиях [7-19], 4 работы в сборниках [20-23] и тезисы российских и международных конференций.

Объем и структура диссертации. Диссертация состоит из введения, семи глав и заключения, изложена на 122 страницах, включая 25 рисунков, 3 таблицы и 109 наименований цитируемой литературы.

8.1. Основные результаты и выводы работы

Создан оригинальный пакет программ МД моделирования, предназначенный для постановки крупномасштабных численных экспериментов. Пакет использован практически для различных задач моделирования [7, 9, 10, И, 14, 16, 17, 19-22].

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

Разработан алгоритм кластерной диапазонной интерполяции для построения потенциалов взаимодействия, описывающих разрыв и образование химических связей. Продемонстрирована работоспособность алгоритма КДИ, в частности, показано, что сохранение полной энергии в модели обеспечивается с контролируемой точностью [12,13,15]. Модель КДИ применена для тестовой задачи (эмпирические связанные состояния в водородной плазме). Построены равновесные парные корреляционные функции частиц и распределения частиц по энергии. Изучены времена жизни связанных состояний [18].

На основе метода функционала плотности построена и исследована модель потенциальной энергии взаимодействия графена, нанотрубки 4x4 и фуллерена Сбо с атомарным и молекулярным кислородом в синглетном и триплетном состоянии [7].

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

2. Оптимизация геометрии системы показывает, что присоединение атомарного кислорода наиболее энергетически выгодно над самой короткой связью углеродного каркаса (эпоксидная структура), а присоединение молекулы кислорода - параллельно самой короткой связи С-С.

3. Наличие пятичленных циклов («двойных связей») приводит к большей активности поверхности относительно присоединения атомарного или молекулярного кислорода.

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

5. Для исследуемых систем минимальное по энергии состояние реагентов всегда трип летное, а связанное состояние продуктов реакции — сингл етное.

Продемонстрирована возможность построения «быстрого» полуэмпирического потенциала взаимодействия атомов углерода и кислорода, воспроизводящего результаты, полученные ab initio. Такой потенциал создается на основе метода сильной связи и схемы кластерной диапазонной интерполяции.

8.2. Достоверность результатов

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

Достоверность результатов вычисления поверхности потенциальной энергии системы кислород+углерод методом функционала плотности проверялась сравнением с экспериментальными данными по энергиям окисления и геометрической конфигурации известных оксидов фуллерена и нанотрубок. Установлено совпадение энергий присоединения атомарного кислорода к фуллеренам и наиотрубкам (2-4 эВ) с экспериментальными оценками и совпадение геометрии оксида фуллерена, имеющего эпоксидную структуру. Результат о безбарьерном присоединении как триплетно-го, так и синглетного атомарного кислорода к бездефектному углеродному каркасу совпадает с экспериментальными фактами. Результат о наличии барьера не менее 1 эВ при присоединении молекулярного триплетного кислорода к поверхности фуллерена качественно совпадает с экспериментальной оценкой. Уменьшение барьера до 0,5 эВ для сингл етного кислорода согласуется с экспериментальным фактом об уменьшении энергии активации реакции с углеродным каркасом при возбуждении кислорода излучением. Наличие барьера не менее 3 эВ для присоединения триплетного кислорода к базовой поверхности графита согласуется с экспериментальным фактом о пассивности базовой поверхности и о роли дефектов и краев в процессе окисления.

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

Достоверность результатов представления поверхности потенциальной энергии системы кислород+углерод с помощью построенного полуэмпирического потенциала проверялась сравнением энергии различных систем вычисленных с помощью построенного потенциала и напрямую методом функционала плотности. Сравнение проводилось для конфигураций, не входящих в набор для подгонки параметров. Сравнение показывает хорошую точность представления поверхности потенциальной энергии (до 0,05 эВ на атом кислорода вблизи равновесия и до 0,1 эВ на атом кислорода при отклонении от равновесия до 5 эВ).

8. Заключение

1. Haynes В. S. A turnover model for carbon reactivity i. development // Combustion and Flame. — 2001.-Vol. 126.-Pp. 1421-1432.

2. Дьячков П. H. Углеродные нанотрубки. — М.: Бином Лаборатория знаний, 2006.

3. Хи Y.-J., Li J.-Q. The interaction of molecular oxygen with active sites on graphite: a theoretical study // Chem. Phys. Lett. — 2004. — Vol. 400. — Pp. 406-412.

4. Montoya A., Mondragon F., Truomj T. N. First-principles kinetics of CO desorption from oxygen species on carbonaceous surface //J. Phys. Chem. B. — 2002. — Vol. 106, no. 16. — Pp. 4236-4239.

5. Froudakis G. E., Schnell M., Miihlhauser M., Peyerimhoff S. D., Andriotis A. N., Menon M., Sheetz R. M. Pathways for oxygen adsorption on single-wall carbon nanotubes // Phys. Rev. B. 2003. - Vol. 68. - Pp. 115435-1-5.

6. Nemukhin A. V., Grigorenko B. L., Topol I. A., Burt S. K. Flexible effective fragment QM/MM method: Validation through the challenging tests //J. Comput. Chem. — 2003.— Vol. 24, no. 12.- Pp. 1410-1420.

7. Валуев И. А., Норман Г. Э., Шуб Б. Р. Механизмы окисления бездефектных поверхностей углеродных наноструктур: влияние кривизны поверхности // Химическая физика. — 2011. Т. 30, № 1. - С. 82-88.

8. Morozov I., Valuev I. Distributed applications from scratch: Using gridmd workflow patterns // Lecture Notes in Computer Science (LNCS). — 2007. — Vol. 4489. — Pp. 199-203.

9. Valuev I. Simulation of hydrogen plasma with cluster multi-range interpolation //J. Phys. A: Math. Gen. 2006. — Vol. 39. — Pp. 4465-4468.

10. Valuev I. A. Reactive potentials for molecular dynamics with cluster multi-range interpolation // Computer Physics Communications. — 2005. — Vol. 169, no. 1-3. — P. 60.

11. Valuev I. A. Gridmd: Program architecture for distributed molecular simulation // Lecture Notes in Computer Science (LNCS). — 2005. — Vol. 3719. — P. 309.

12. Kuksin A. Y., Morozov I. V., Norman G. E., Stegailov V. V., Valuev I. A. Standards for molecular dynamics modelling and simulation of relaxation // Molecular Simulation. — 2005. — Vol. 31.-Pp. 1005-1017.

13. Валуев И. А., П.Н.Дьячков, Каклюгин А. С., Норман Г. Э. Квантовохимическое моделирование строения, электронной структуры и энергии образования интеркалированного литием и чистого фуллерена // Ж. неорг. хим. — 2000. — Т. 44, № 3. — С. 472-477.

14. Валуев И. А., Каклюгин А. С., Норман Г. Э. Численное моделирование распределения электронной плотности поверхностных электронных состояний // Поверхность. Рентгеновские, синхротронные и нейтронные исслед.— 1999. — № 7.— С. 56.

15. Валуев И. А., Каклюгин А. С., Норман Г. Э. Структура одноэлектронных уровней и оболочек фуллерена Cgo // Химическая физика.— 1996.— Т. 15, № 5.— С. 26-34.

16. Валуев И. А., Каклюгин А. С., Норман Г. Э. Электронная структура и химические свойства фуллерена Cgo // Известия РАН. Серия физическая. — 1996. — Т. 60, № 5. — С. 2432.

17. Валуев И. А., Каклюгин А. С., Норман Г. Э. Энергии и электронные волновые функции больших молекул с ковалентными связями // Химическая физика.— 1996.— Т. 15, № 95.- С. 53.

18. Валуев И. А., Норман Г. Э., Шуб Б. Р. Полуэмпирическая модель взаимодействия О и 02 споверхностью углеродных наноструктур // Тезисы XXI симпозиума Современная химическая физика. — 2009. — С. 7.

19. Valuev I. A. Combined Multi-Range Interpolated Potential for MD Simulations // Computational Physics / Ed. by X. G. Zhao, S. Jiang, X. J. Yu.~ Rinton Press Inc., USA, 2005.— Pp. 154-157.

20. Морозов И. В., Валуев И. А. Пакет программ молекулярно-динамического моделирования с интерфейсом для коллективного доступа // Научный сервис в сети Интернет: технологии распределенных вычислений. — Москва: МГУ, 2005. — С. 134-135.

21. Busnengo H. F., Salin A., Dong W. Representation of the 6d potential energy surface for a diatomic molecule near a solid surface //J. C%em. Phys. — 2000. — Vol. 112. — P. 7641.

22. Stuart S. J., Tutein A. B., Harrison J. A. A reactive potential for hydrocarbons with intermolecular interactions // J. Chem. Phys. — 2000. — Vol. 112. — P. G472.

23. Sorbie K., Murrell J. Analytical potentials for triatomic molecules from spectroscopic data // Mol. Phys. 1975. - Vol. 29, no. 5. - P. 1387.

24. Papaconstantopoulos D. A., Mehl M. J. The Slater-Koster tight-binding method: a computationally efficient and accurate approach // J. Phys.: Condens. Matter. — 2003.—- Vol. 15.— P. R413.

25. Dewar M. J. S., Thiel W. Ground states of molecules. 38.the mndo method, approximations and parameters // J. Am. Chem. Soc.— 1977. — Vol. 99, no. 15. —P. 4899.

26. Stewart J. Optimization of parameters for semiempirical methods i. method //J. Comput. Chem. 1989. - Vol. 10, no. 2. - P. 209.

27. Kuger T., Elstner M., Schiffels P., Frauenheim T. Validation of the density-functional based tight-binding approximation method for the calculation of reaction energies and other data // JCP.~ 2005. -Vol. 122.-Pp. 114110-1.

28. Lee S. M., Lee Y. H., Hwang Y. G., Hahn J. R., Kang H. Defect-induced oxidation of graphite // Phys. Rev. Lett. — 1999. — Vol. 82, no. 1.— Pp. 217-220.

29. Hurt R. H., Haynes B. S. On the origin of power-law kinetics in carbon oxidation // Proceedings of the Combustion Institute. — 2005. — Vol. 30. — Pp. 2161-2168.

30. Jelea A., Marinelli F., Ferro Y., Allouche A., Brosset C. Quantum study of hydrogen-oxygen-graphite interactions // J. Comput. Chem. — 2004. — Vol. 42. — Pp. 3189-3198.

31. Sendt K., Haynes B. S. Density functional study of the chemisorption of O2 on the armchair surface of graphite // Proceedings of the Combustion Institute. — 2005. — Vol. 30. — Pp. 21412149.

32. Radovic L. R. The mechanism of CO2 chemisorption on zigzag carbon active sites: A computational chemistry study // Carbon. — 2005. — Vol. 43. — Pp. 907-915.38.