CC BY
26
Силовые поля для молекулярно-динамического моделирования процесса напыления пленок диоксида кремния
Ф. В. Григорьев1,2
1 Московский государственный университет имени М.В. Ломоносова, Научно-исследовательский вычислительный центр.
Россия, 119991, Москва, Ленинские горы, д. 1, стр. 4.
2 Национальный исследовательский ядерный университет «МИФИ».
Россия, 115409, Москва, Каширское ш., д. 31.
E-mail: afedor.grigoriev@gmail.com
Статья поступила 22.06.2015, подписана в печать 11.08.2015.
Проведено сравнение двух силовых полей, предназначенных для молекулярно-динамического моделирования процесса напыления тонких пленок диоксида кремния. Проведен анализ структурных характеристик (плотность, радиальная функция распределения) кластера стеклообразного диоксида кремния, используемого в качестве подложки, и напыленной пленки. Показано, что для моделирования процесса напыления наиболее подходит силовое поле DESIL, в котором взаимодействие Ван-дер-Ваальса описывается потенциалом Леннарда-Джонса.
Ключевые слова: суперкомпьютерное моделирование, молекулярная динамика, тонкие пленки, диоксид кремния, напыление, силовое поле.
УДК: 539.231. PACS: 81.15.Aa.
Введение
Многослойные оптические покрытия широко используются в современных оптических приборах и устройствах. Как правило, такие покрытия состоят из чередующихся слоев с высоким и низким показателями преломления [1]. Одним из наиболее распространенных материалов для слоев с низким показателем преломления является диоксид кремния.
Структурные и оптические свойства пленок существенным образом зависят от технологических параметров процесса напыления — энергии и углового распределения скоростей атомов, напыляемых на подложку, ее температуры, давления и состава газа в вакуумной камере и т. д. В настоящее время одним из наиболее перспективных методов, позволяющих получить однородные и плотные пленки, является метод высокоэнергетического напыления (ion beam sputtering, IBS), в котором энергия атомов кремния, падающих на подложку и ранее напыленные слои пленки, достигает нескольких десятков элек-тронвольт [2]. Актуальной задачей с точки зрения совершенствования технологии высокоэнергетического напыления является исследование зависимости структурных и оптических свойств пленки от величин технологических параметров процесса напыления. Экспериментальное исследование такой зависимости затруднено вследствие малой толщины и неупорядоченной структуры пленки, поэтому использование методов атомистического моделирования в той области представляется целесообразным.
При выборе метода моделирования следует принимать во внимание, что структурные неоднородности пленки, влияющие на ее оптические свойства, имеют характерные размеры порядка нескольких
десятков нанометров. Кубический кластер диоксида кремния с характерным размером порядка 100 нм сдержит около 108 атомов. В настоящий момент моделирование кластеров с таким числом атомов возможно только в рамках классической молекулярной динамики (МД) или метода Монте-Карло (МС) с использованием технологий параллельных вычислений, реализованных на суперкомпьютерных комплексах.
Моделирование роста пленок методами атомистического уровня проводится с 1970-х гг. (см. работу [3] и цит. лит. 3-17). В последние годы в связи с развитием вычислительных средств — прежде всего многопроцессорных кластеров — появилась возможность моделирования кластеров относительно больших размеров и использования трехчастич-ных силовых полей, позволяющих описывать химические реакции на поверхности пленки. В работе [4] с использованием таких силовых полей моделировался рост пленки диоксида кремния. Был продемонстрирован эффект увеличения плотности пленки с ростом кинетической энергии осаждаемых атомов кремния. В работе [5] методом МД моделировались поверхностные свойства пленок диоксида титана с кластерами, включающими до 2.4 • 104 атомов. В работе [6] структурные характеристики пленок (плотность, шероховатость, концентрация пор) исследовались методами МД и МС на масштабах порядка десятков нанометров с числом атомов в кластерах около 106 .
Адекватность МД- или МС-моделирования, так же как и его численная эффективность, в значительной степени определяются силовым полем, выбранным для расчета потенциальной энергии взаимодействия атомов пленки и подложки. В работе [7]
дан достаточно подробный обзор различных силовых полей, созданных для диоксида кремния, и сделан вывод, что для моделирования напыления наиболее подходящим является двухчастичное силовое поле (каждое слагаемое в выражении для потенциальной энергии зависит от расстояния между двумя атомами) с относительно простой функциональной формой. В [7] описано оригинальное силовое поле ОЕ81Ц в рамках которого потенциальная энергия пары атомов представляется суммой кулоновского слагаемого, описывающего взаимодействие точечных зарядов, центрированных на атомах, и слагаемого с потенциалом Леннарда-Джонса, описывающим неэлектростатическую часть потенциальной энергии взаимодействия. ЭЕ81Ь с хорошей точностью воспроизводит ключевые структурные характеристики стеклообразного диоксида кремния (плотность при нормальных условиях, тетраэдрический мотив в структуре непрерывной неупорядоченной сетки, длину связи Б1—0, валентные углы 81-0-81 и 0-81-0, положения первого и второго пиков радиальной функции распределения). В то же время в потенциале Леннарда-Джонса энергия обменного взаимодействия на малых расстояниях описывается слагаемым пропорциональным 1/г12 (г — расстояние между атомами) вместо более точной экспоненциальной зависимости от расстояния. Такая неточность может сказываться на описании кинематики взаимодействия атомов вследствие большой начальной энергии атомов (несколько десятков элек-тронвольт), сближающихся на малые расстояния.
В настоящей работе представлено силовое поле ЭЕ81Ь_Б, в рамках которого для описания неэлектростатической части потенциальной энергии взаимодействия атомов используется потенциал Букин-гема [8], содержащий экспоненциальное слагаемое для расчета обменного взаимодействия. По результатам моделирования кластеров различных размеров получены структурные характеристики стеклообразного диоксида кремния в рамках ЭЕ81Ь_Б. Обсуж-
даются преимущества и недостатки силовых полей ЭЕ81Ь и ЭЕ81Ь_Б.
1. Метод моделирования
Моделирование проводилось методом классической МД. Процедура напыления (инжекция атомов кремния и кислорода, распределение их начальных координат и скоростей, граничные условия, ограничения на движение атомов подложки), реализованная в оригинальной программе КиУАЬЭА, подробно описана в [7]. Для молекулярно-динамической части процедуры КиУАЬЭА использует программу 0ДОМАС8 [9] в качестве внешней.
Моделирование проводилось на суперкомпьютере «Ломоносов» МГУ имени М.В. Ломоносова [10].
2. Результаты и обсуждение
Потенциальная энергия парного взаимодействия с использованием потенциала Букингема записывается в виде
ив = Ш + Ац еМ-ЬцГц) - Ц.
(1)
Силовое поле для диоксида кремния с функциональной формой (1) предложено в [11]. Параметры определялись по результатам квантово-химического моделирования небольших молекулярных кластеров и позволили воспроизвести их геометрические характеристики. Однако для моделирования процесса напыления силовое поле должно воспроизводить и объемные структурные характеристики, такие как плотность и радиальная функция распределения (РФР). Параметры [11] были взяты в качестве начальных и затем изменялись с тем, чтобы удовлетворить этим требованиям. Процедура подгонки параметров описана в [3].
Зависимость потенциальной энергии взаимодействия атомов кремния и кислорода от расстояния между ними показана на рис. 1. Видно (рис. 1, а),
С/104, кДж/моль
0.05 0.07 0.09 0.11 Я, нм -4
0.18 Я та
-4
Рис. 1. Потенциальная энергия и взаимодействия атомов кремния и кислорода в зависимости от расстояния к между ними. Сплошная линия — силовое поле ЭЕ81Ь с потенциалом Леннарда-Джонса,
штриховая — силовое поле ЭЕ81Ь_Б с потенциалом (1)
что при малых расстояниях потенциал Букингема имеет нефизическое поведение, что связано с резким возрастанием слагаемого, обратно пропорционального шестой степени расстояния. Применительно к задаче напыления тонких пленок такое нефизическое поведение создает трудности, поскольку при энергии налетающих на пленку атомов около 10-100 эВ атомы могут сближаться на малое расстояние. Барьер в потенциале Букингема должен быть достаточно большим, чтобы при взаимодействии высокоэнергетических атомов кремния с подложкой и ранее напыленными слоями пленки межатомное расстояние не уменьшилось до величин, соответствующих нефизическому поведению потенциала. Значение энергии 100 эВ соответствует 9.65 • 103 кДж/моль, что составляет менее трети от величины барьера — порядка 3.5 • 104 кДж/моль (рис. 1, а). Как показало последующее моделирование, этого достаточно для предотвращения попадания в нефизическую область потенциала Букингема.
Полученные нами (силовое поле ЭЕ81Ь_В) и известные ранее из [11] параметры потенциала (1) приведены в табл. 1. Заряды в рамках ЭЕ81Ь_В существенно меньше, чем в [11], в то же время величина коэффициента сц для пары Б1—0 заметно выше. В рамках ЭЕ81Ь_В потенциал Букингема убывает быстрее, чем в [11], из-за большей величины параметра Ьц . Параметры силового поля ЭЕ81Ь с потенциалом Леннарда-Джонса приведены в табл. 2.
При локальной оптимизации с периодическими граничными условиями в рамках силового поля ЭЕ81Ь_В от начальной структуры, соответствующей кристаллу а -кварца, геометрия кластера в основном сохраняется (рис. 2), некоторые искажения связаны с тем, что в рамках ЭЕ81Ь_В длина связи Б1-0, соответствующая минимуму потенциальной энергии
Рис. 2. Атомистическая структура кристаллического диоксида кремния в рамках силового поля ЭЕ81Ь_В
кластера, равна 0.168 нм, что на 0.004 нм больше экспериментальной величины. Плотность при локальной оптимизации не меняется и соответствует значению 2.65 г/см3 кристаллической модификации а -кварца.
Структура подложки стеклообразного диоксида кремния была получена с использованием силового поля ЭЕ81Ь_В и МД-процедуры, описанной в [7, 12]. Кластер подложки толщиной 2 нм содержал 9 • 104 атомов, размеры кластера по горизонтали 26 х 23 нм.
Радиальная функция распределения атомов подложки показана на рис. 3. При использовании потенциала Букингема (силовое поле ЭЕ81Ь_В) высота первого пика почти вдвое меньше, чем у РФР для потенциала Леннарда-Джонса (силовое поле ЭЕ81Ь), а его положение смещено вправо примерно на 0.004 нм. В то же время высота второго пика
Таблица 1
Параметры силовых полей DESIL_B и поля [11]
Для описания взаимодействия Ван-дер-Ваальса используется потенциал Букингема (1)
1-1 А', кДж/моль Ьц , 1/нм 011, кДж-нм6/моль Атомные заряды, е
ЭЕ81Ь_В [10] ЭЕ81Ь_В [10] ЭЕ81Ь_В [10] ЭЕ81Ь_В [10]
0-0 5 • 105 1.34 • 105 37.0 27.6 0 1.69 • 10~2 1.6 (81) -0.8 (0) 2.4 (81) -1.2 (0)
81-81 6 • 105 — 38.0 — 0 —
81-0 2.1 • 106 1.737 • 106 60.8 48.73 2.0 • 10~3 1.29 • 10~2
Таблица 2
Параметры силового поля DESIL [7]
Для описания взаимодействия Ван-дер-Ваальса используется потенциал
Леннарда-Джонса
81-0 81-81, 0-0 Атомные заряды, е
с12, кДж-нм12/моль 4.6 • 10~8 1.5 • 10~6 ?81 = 1-3 ц0 = -0.65
с6, кДж-нм6/моль 4.2 • 10~3 5.0 • 10~5
16
12
T(R), отн. ед. Si-O
О-О
Si-Si
0.16 0.20 0.24 0.28
0.32 0.36 R, нм
3.0
2.0
1.0
d, г/см3
1.8
3.6
5.4
7.2 H, нм
Рис. 3. Радиальная функция распределения T (отн. ед.) для стеклообразного диоксида кремния в зависимости от расстояния R между атомами. Сплошная линия — силовое поле ЭЕ81Ь, штриховая — силовое поле ЭЕ81Ь_Б
в рамках ЭЕ81Ь_Б существенно ниже, их положения примерно совпадают. Как первый, так и второй пик РФР в случае потенциала Букингема более размыты, что означает большую флуктуацию длин связей Б1—0 (первый пик) и расстояний между ближайшими атомами кислорода (второй пик). Площади под первым пиком РФР совпадают для обоих потенциалов, что означает одинаковую концентрацию четырехкоординированных атомов кремния, как и должно быть. Положения вторых пиков совпадают с экспериментальными [13], положение первого пика РФР с потенциалом Леннарда-Джонса ближе к экспериментальному.
В соответствии с методикой, описанной в [7], было проведено моделирование процесса напыления тонкой пленки. Группы атомов кремния и кислорода инжектировались в область моделирования каждые
6 пикосекунд, температура подложки поддерживалась равной 500 К, энергия осаждаемых атомов кремния 10 эВ, энергия атомов кислорода соответствует температуре 1000 К.
Зависимости плотности от суммарной толщины подложки и пленки показаны на рис. 4. Для силового поля ЭЕ81Ь зависимость близка к той, что получена ранее в [7]: плотность пленки не превышает 2.5 г/см3, толщина переходных слоев подложка-пленка и пленка-вакуум около 1.5 нм. Для зависимости, полученной с использованием силового поля ЭЕ81Ь_Б, характерно существенное — на 1 г/см 3 превышение плотности пленки над плотностью подложки.
Анализ структуры показал, что это превышение обусловлено появлением большого количества — до 40% от общего числа — пятикоординированных атомов кремния, при этом среднее расстояние связи 81-0 повышается от 0.168 до 0.178 нм. Наличие атомов кремния с такими характеристиками связей 81-0 не имеет экспериментального подтверждения и, по всей видимости, является артефактом силового
Рис. 4. Зависимость плотности d напыленной пленки с учетом подложки (начальная толщина 2 нм) от расстояния H от дна подложки. Сплошная линия — силовое поле ЭЕ81Ь, штриховая — силовое поле ЭЕ81Ь_Б
поля. Возможно, образование конфигураций с пяти-координированными атомами 81 в большой концентрации связано с пологой формой потенциала Бу-кингема вблизи минимума энергии взаимодействия атомов кислорода и кремния (рис. 1, б), что приводит к большему разнообразию геометрии ближнего порядка в напыленной пленке в сравнении со стеклообразным состоянием. Отметим, что при напылении с использованием потенциала Леннарда-Джонса концентрация пятикоординированных атомов кремния не превышает четверти процента. Концентрация однокоординированных атомов кислорода и трехко-ординированных атомов кремния при использовании потенциала Леннарда-Джонса в несколько раз ниже, чем при использовании потенциала Букингема.
Таким образом, для моделирования процесса напыления тонких пленок 8Ю2 более предпочтительным выглядит силовое поле ОЕ81Ц в то время как для моделирования объемных образцов стеклообразного 8Ю2 могут быть использованы оба силовых поля.
Было проведено сравнение численной эффективности силового поля ЭЕ81Ь с потенциалом Леннар-да-Джонса и силового поля ЭЕ81Ь_Б с потенциалом Букингема. Результаты представлены на рис. 5.
N
Рис. 5. Зависимость времени моделирования t, подложки от числа ядер N. Сплошная линия — силовое поле ЭЕ8!Ь, штриховая — силовое поле ЭЕ8!Ь_Б
0
Длина траектории МД-моделирования кластера подложки, содержащего 9 • 104 атомов, составляла 500 пс с шагом 1 фс. Моделирование проводилось с использованием программы 0Я0МЛС8 в ЫРТ ансамбле (постоянные число частиц, давление и температура) с периодическими граничными условиями.
Из зависимостей на рис. 5 видно, что при прочих равных условиях время расчета с потенциалом Букингема примерно вдвое выше времени расчета с потенциалом Леннарда-Джонса. Такой результат ожидаем, поскольку вычисление значения экспоненциальной функции требует большего времени, чем вычисление степенной. Численная эффективность, рассчитанная как а(Д) = Тл^а6 для потенциала Букингема падает от а(16) = 1 до 0.7 при N = 32 и далее остается примерно на том же уровне. Для потенциала Леннарда-Джонса а(32) = 0.8 и далее снижается в пределах 0.05 с ростом числа ядер до 128.
Заключение
В работе проводится сравнение двух силовых полей для атомистического моделирования процесса напыления тонких пленок диоксида кремния. В рамках разработанного ранее [7] силового поля ЭЕ81Ь для описания взаимодействия Ван-дер-Ва-альса используется потенциал Леннарда-Джонса, в рамках представленного в настоящей работе поля ЭЕ81Ь_Б используется потенциал Букингема [8], содержащий экспоненциально зависящее от расстояния между атомами слагаемое. Показано, что оба силовых поля хорошо воспроизводят структуру (длина связи, тетраэдры 8Ю4, плотность, положение первых пиков радиальной функции распределения) стеклообразного диоксида кремния. При моделировании напыления тонкой пленки в рамках поля ЭЕ81Ь_Б было найдено, что ее плотность достигает 3.2 г/см3, что превышает экспериментальные
значения 2.1-2.4 г/см3. Анализ структуры выявил, что увеличение плотности обусловлено большой — до 40 % — долей пятикоординированых атомов кремния.
Таким образом, несмотря на то, что поле DESIL_B хорошо воспроизводит структуру стеклообразного диоксида кремния, для моделирования процесса напыления предпочтительнее использовать силовое поле DESIL с потенциалом Леннар-да-Джонса.
Работа выполнена при финансовой поддержке Российского научного фонда (грант 14-11-00409).
Список литературы
1. Piegari A., Flory F. Optical Thin Films and Coatings. Cambridge, 2013.
2. Tabata A., Matsuno N., Suzuoki Y. // Thin Solid Films. 199б. 289, N 1-2. P. 84.
3. Smith R.W., Srolovitz D.J. // Appl. Phys. 199б. 79, N 3. P. 1448.
4. Taguchi M., Hamaguchi S. // Thin Solid Films. 2007. 515, N 12. P. 4879.
5. Bahramian A. // Surf. Interface Anal. 2013. 45, N 11-12. P. 1727.
6. Turowski M., Jupea M., Melzig T. et al. // Thin Solid Films. 2015. 592. P. 240.
7. Grigoriev F.V., Sulimov A.V., Kochikov I.V. et al. // Int. J. High Perf. Comp. Appl. 2015. 29, N 2. P. 184.
8. Buckingham R.A. // Proc. Royal Soc. London. Ser. A. Math. and Phys. Sci. 1938. 168, N 933. P. 2б4.
9. http://www.gromacs.org
10. http://parallel.ru/cluster
11. Best B.W.H., Kramer G.R., Santen R.A. // Phys. Rev. Lett. 1990. 64, N 1б. P. 1955.
12. Григорьев Ф.В., Сулимое В.Б., Кондакова O.A. и др. // Вестн. Моск. ун-та. Физ. Астрон. 2013. № 3, С. 80. (Grigoriev F.V., Sulimov V.B., Kondakova O.A. et al. // Moscow University Phys. Bull. 2013. 68, N 3. P. 259.)
13. Johnson P.A.V., Wright A.C., Sinclair R.N. // J. Non-Cryst. Solids. 1983. 58, N 1. P. 109.
Force fields for molecular dynamics simulation of the deposition of a silicon dioxide film F. V. Grigoriev1,2
1 Research Computing Center, Lomonosov Moscow State University, Moscow 119991, Russia.
2 National Engineering Physics Institute «MEPhI», Moscow 115409, Russia. E-mail: a fedor.grigoriev@gmail.com.
In this work we compare two force fields that are intended for the molecular dynamics simulation of the process of the deposition of silicon dioxide thin films. Analysis of the structural characteristics (the density and radial distribution function) of a glassy silicon dioxide cluster that was used as a substrate and a deposited film is carried out. It is shown that the DESIL force field in which the Van der Waals interaction is described by the Lennard-Jones potential turns out to be more suitable for modeling the process of deposition.
Keywords: high-performance simulation, molecular dynamics, thin film growth, deposition process, silicon dioxide force fields. PACS: 81.15.Aa. Received 22 June 2015.
English version: Moscow University Physics Bulletin 6(2015). Сведения об авторе
Григорьев Федор Васильевич — канд. хим. наук, вед. науч. сотрудник; тел.: (495) 939-32-53, e-mail: fedor.grigoriev@gmail.com.
