CC BY
63
ФИЗИКА ТВЁРДОГО ТЕЛА
УДК 548.4 + 519.24
СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РОСТА КРИСТАЛЛА С УЧЁТОМ ЕГО АТОМНОЙ СТРУКТУРЫ
© 2009 г. С.С. Гуськов, М.А. Фаддеев, Е.В. Чупрунов
Нижегородский госуниверситет им. Н.И. Лобачевского guskovss@mail.ru
Поступила в редакцию 28.07.2009
Предложен оригинальный алгоритм численного моделирования роста грани кристалла с произвольной атомной структурой. Кристалл может включать точечные дефекты и дислокации. Алгоритм допускает введение в ростовую систему примеси, замещающей атомы определенных сортов. Предусмотрена возможность использования модели для исследования механизмов образования зонарной микропериоди-ческой структуры монокристаллов, связанных с неоднородным распределением примеси.
Ключевые слова: рост монокристаллов, атомная структура, потенциал межатомного взаимодействия, дефекты, зонарная микроструктура, примесь, математическое моделирование.
Введение
Разнообразие процессов, происходящих в ходе роста кристаллов, приводит к затруднениям при теоретическом описании кинетики роста. Как следствие, часто используется исследование роста кристалла путем прямого статистического моделирования методами Монте-Карло [1-3]. Это позволяет учитывать вероятностный характер присоединения и отрыва атомов в разных положениях на растущей грани кристалла.
Обычно модели роста кристаллов строятся для конкретного типа кристаллической структуры [4]. Первые опыты проводились с простейшей кубической структурой (кристалл Кос-селя) [5]. Известны модели, рассматривающие простейшие грани ОЦК и ГЦК одноатомных решеток [6] и практически важные грани некоторых полупроводников (как правило, с алмазоподобной структурой) [7, 8]. Нами применен другой подход, пригодный для кристаллов не только кубической, но и более низких синго-ний. В данной работе предлагается оригинальный алгоритм статистического моделирования роста произвольной грани кристалла, для которого известна атомная структура, а также вид и численные параметры функций, аппроксимирующих потенциалы межатомного взаимодействия. Исследуемый кристалл может содержать дефекты структуры, такие как атомы примеси или дислокации с заданным направляющим
вектором и вектором Бюргерса. Наша модель может быть использована для исследований особенностей роста кристалла в присутствии примесей, которые играют ключевую роль в процессе образования зонарной микроструктуры монокристаллов [9, 10].
Геометрия области расчета
Параметры элементарных ячеек и координаты атомов для многих кристаллических структур определены с высокой точностью с помощью методов рентгеноструктурного анализа. Описания структур десятков тысяч неорганических и сотен тысяч органических соединений доступны в соответствующих банках данных [11, 12]. Полагаем известными: а) компоненты базисных векторов а, Ъ, с элементарной ячейки в кристаллофизической (т.е. декартовой с ортами е\, е2, е3) системе координат (рис. 1); б) количество атомов в элементарной ячейке N в) безразмерные кристаллографические координаты центров атомов р (г = 1, 2, ..., Ж). Рассматриваемая грань кристаллической структуры определяется индексами Миллера (Ик/).
Будем рассматривать область пространства, ограниченную цилиндрической поверхностью и двумя параллельными плоскостями - основаниями цилиндра (рис. 2). Ось цилиндра параллельна нормали к семейству плоскостей (ИИ) и проходит через начало декартовой системы ко-
ряют условиям (4) и (5), получим набор координат всех атомных позиций, входящих в рассматриваемый цилиндрический фрагмент кристаллической структуры.
Энергии межатомного взаимодействия
Взаимодействие между двумя атомами будем описывать функцией потенциальной энергии Цт), где г - межатомное расстояние. Если между атомами существует устойчивая химическая связь, то кривая и(г) имеет минимум г = г0, соответствующий равновесному межатомному расстоянию [13]. Величина г0 может быть интерпретирована как длина химической связи.
Для проведения расчетов в рамках нашей модели необходимо знание энергий связи еАв всех возможных пар атомов АВ. В случае ионного характера химической связи для энергии взаимодействия между двумя атомами А и В используется выражение следующего вида [13]:
(
% А % В е
ГАВ
(
+ Ь ехр
ГАВ ' ГЬ )
(6)
ординат. Основания цилиндра параллельны семейству плоскостей (Ик/).
Уравнение оси цилиндра запишем в виде
[пшт] = 0. (1)
Здесь т - радиус-вектор точки пространства, а пш - единичный вектор нормали к семейству плоскостей (Ик/).
Уравнения двух параллельных плоскостей -оснований цилиндра имеют вид
пшт = 1ъ пшт = d2, (2)
где 11 и 12 - расстояния от начала координат до нижнего и верхнего оснований цилиндра соответственно. Для высоты цилиндра В справедливо соотношение
В = (12 - (3)
Условия принадлежности точки М с радиус-вектором тм рассматриваемой области могут быть записаны в следующем виде:
| [пштм] | < Я, (4)
1 < пштм < 12, (5)
где Я - радиус цилиндра.
Размножив элементарную ячейку кристаллической структуры с помощью векторов трансляций Тип, = иа + уЪ + wc (и, у, w - независимые целые числа) и выбрав из полученного множества точек те точки, которые удовлетво-
Здесь ZAe = ца и Zвe = цв (е - элементарный заряд) представляют собой электрические заряды ионов, находящихся на расстоянии гАв друг от друга; К, Ь и гь - постоянные. Первое слагаемое (6) соответствует потенциалу кулоновского взаимодействия, а второе - потенциалу сил отталкивания, возникающих при сближении атомов. Исследования показали, что среднее значение гь = 0.345 А, а отклонения от среднего значения для большинства ионных кристаллов не превышают 6% [14].
Введем параметр р, характеризующий максимальное расстояние межатомного взаимодействия. Если расстояние между двумя атомами г превышает параметр р, то энергией взаимодействия данной пары атомов в рамках нашей модели будем пренебрегать.
Неполярные межатомные связи будем рассматривать в приближении ближайших соседей. При этом параметр р равен расстоянию между ближайшими соседями, а соответствующая энергия межатомного взаимодействия выражается определенной величиной еАв. Таким образом, ионное взаимодействие полагается дально-действующим, а другие типы взаимодействий рассматриваются как короткодействующие.
Полная энергия взаимодействия определенного атома А с остальными атомами вычисляется следующим образом. Находятся все атомы В,-, лежащие от данного на расстоянии, меньшем либо равном принятой дальности действия потенциала р. Для каждого атома В, вычисляется значение энергии вАв,- взаимодействия с атомом
2
А. Заметим, что, вообще говоря, значения вдвь могут быть как положительными, так и отрицательными, причем отрицательные значения вАв,-соответствуют притяжению атомов А и В,, а положительные - отталкиванию. Полная энергия данного атома вА представляет собой сумму вкладов всех атомов В,-:
8 А = X 8 АЫ . (7)
г
Рабочая область и элементарные события
Первый этап моделирования состоит в формировании исходной конфигурации поверхности кристалла, т.е. в заполнении части позиций атомами. Простейшая (плоская) исходная конфигурация соответствует заполнению атомами всех позиций, удовлетворяющих условию 7 < 2р, где 7 - расстояние от данной точки до нижнего основания цилиндрической области моделирования, а р - введенный ранее характерный радиус действия потенциала межатомного взаимодействия.
Фронтом роста будем называть поверхность, отделяющую занятые атомами позиции от незанятых позиций в области моделирования. Вообще говоря, фронт роста не является плоскостью.
Рабочей областью будем называть часть области моделирования, ограниченную боковой цилиндрической поверхностью, а также двумя поверхностями, полученными смещением фронта роста вдоль оси области на величины ±р. Для простоты расчетов в качестве нижней и верхней границ рабочей области можно рассматривать плоскости, перпендикулярные оси области расчета и расположенные от нижнего основания на расстояниях 1тп - р и 1тах + р. Величины 1тп и 1тах показаны на рис. 3.
Процесс роста кристалла представим как последовательность элементарных событий. При
определенном положении фронта роста в каждой позиции рабочей области в данный момент времени может произойти то или иное элементарное событие.
Для позиции рабочей области, занятой атомом, введем элементарные события двух типов - отрыв атома и пустое событие. Для незанятой позиции возможны следующие элементарные события - присоединение атома какого-либо химического элемента, находящегося в маточной среде, и пустое событие.
Вероятности элементарных событий определяются физико-химическими процессами на границе фаз. В [3] показано, что вероятность Ж12 перехода системы из состояния 1 в состояние 2 может выражаться через потенциальные энергии системы Е1 и Е2 в состояниях 1 и 2 соответственно следующим образом:
Щ2--------------1----------, (8)
1 + ехр[(Е2 - Е)/ кТ]
где к - постоянная Больцмана; Т - абсолютная температура.
При некоторых разумных допущениях [15] вероятность присоединения атома 7-го типа в данную позицию Щ+ и вероятность отрыва атома 7-го типа из данной позиции Wj. могут быть записаны в виде
WJ + =
Wj- =
1 + /(Дц) ехр[- Щ / кТ] ’
1
1 + (/(Ац))-1 ехр[ £ /кТ]
(9)
где вj - энергия атома 7-го типа, вычисляемая с помощью (7); С - относительная концентрация данного типа ростовых единиц в слое среды у поверхности кристалла. Функция ХЛц) определяется следующим выражением:
/(Дц) = ехр[(Б1/2у - Дц)/кТ]. (10)
2
Рис. 3. Фронт роста и рабочая область. Цифрами обозначены: 1 - позиции, заполненные атомами; 2 - незаполненные позиции; 3 - фронт роста; 4 - точные границы рабочей области; 5 - приближенные границы рабочей области. - наименьшее из расстояний от нижнего основания области расчета до не заполненных атомами позиций, /ітах - наибольшее из расстояний от нижнего основания области расчета до позиций, заполненных атомами
Здесь вуу - энергия атома 7-го типа в положении полукристалла. Эту величину можно вычислить, разделив на два полную энергию атома 7-го типа в объеме кристалла. Такое вычисление для каждой позиции и для всех возможных типов атомов проводится один раз перед началом моделирования роста. Величина Лц представляет собой разность химических потенциалов маточной среды и кристалла. Этой величиной принято характеризовать отклонение физической системы «кристалл - среда» от равновесия [16]. Значение Лц и концентрации атомов всех типов в среде задаются перед началом расчета.
Отклонение от равновесия определяется в случае роста из паровой фазы давлением пара, в случае роста из раствора - концентрацией кристаллизующегося вещества в растворе, при росте из расплава - температурой расплава [16]. Например, для пара - идеального газа и идеального (разбавленного) раствора выполняется соотношение
Дц = кТ 1п(ст +1) . (11)
В зависимости от типа среды параметр о представляет собой относительное пересыщение пара (Р - Р0)/Р0 или относительное пересыщение раствора (С - С0)/С0 , где Р0 и Со - равновесные значения давления пара и концентрации раствора соответственно. Для переохлажденной жидкости (расплава)
Дц = ДЯст, (12)
где ЛН - удельная теплота плавления в расчете на одну частицу, параметр о = (Т0 - Т)/Т0 - относительное переохлаждение расплава, То -температура плавления.
Алгоритм моделирования роста грани
Перед началом моделирования необходимо сделать выбор структуры, задать список химических элементов (ростовых единиц), входящих в состав данной структуры, функции потенциалов взаимодействия атомов, индексы исследуемой грани, размеры области моделирования и количество шагов алгоритма Монте-Карло (МКШ), в течение которых будет выполняться расчет.
После этого осуществляется инициализация модели, в ходе которой вычисляются координаты атомных позиций в области моделирования и создается начальная конфигурация грани кристалла. Затем происходит выполнение заданного количества МКШ, каждый из которых представляет собой последовательность следующих действий:
1. Составляется список атомных позиций в рабочей области (рабочих точек).
2. Для каждой рабочей точки (в зависимости от того, занята она атомом или не занята) опре-
деляется набор возможных элементарных событий и вычисляются их вероятности с помощью выражений (9).
3. В зависимости от значения сгенерированного псевдослучайного числа выполняется то или иное элементарное событие.
МКШ завершен. На следующем шаге повторяется та же последовательность действий.
После завершения моделирования строится изображение рассматриваемого участка структуры, оцениваются некоторые численные параметры структуры, такие как количество присоединений, отрывов, средняя скорость роста.
Моделирование роста некоторых кристаллических образцов
Разработанный алгоритм предназначен для моделирования развития граней кристаллического образца в процессе его роста. Для иллюстрации работы алгоритма рассмотрим кинетику роста граней различной ориентации у кристалла серебра. Данная кристаллическая структура состоит из атомов одного сорта, расположенных в узлах гра-нецентрированной кубической решетки. Будем учитывать только взаимодействие ближайших соседей. Величину энергии межатомного взаимодействия примем равной -0.2 эВ, что характерно для кристаллизации металлов [16].
Зависимости нормальной скорости роста у от отклонения от равновесия Лц/кТ для граней (100), (110) и (111) данной кристаллической структуры, полученные с помощью нашего алгоритма, приведены на рис. 4. Известно, что в огранке кристаллов доминируют грани с наименьшей скоростью роста [16]. Из зависимостей на рис. 4 видно, что появление на кристаллическом образце граней, симметрично эквивалентных (111), более вероятно, чем граней, симметрично эквивалентных (110) и (100). Вид зависимостей у(Лц/кТ) при близких к нулю значениях аргумента свидетельствует о том, что грани (110) и (100) являются шероховатыми и их рост осуществляется по нормальному механизму, в то время как грань (111) является гладкой и растет послойно. Такой вывод подтверждают и изображения атомных конфигураций на поверхностях граней (111), (110) и (100) (рис. 5).
При кристаллизации дигидрофосфата калия (КЭР) в качестве ростовых единиц могут рассматриваться ионы калия с зарядом +е и группы Р04 с суммарным зарядом -е, образующие связи ионного характера. В качестве тестирования предлагаемого алгоритма нами было проведено моделирование движения ступени по поверхности грани призмы кристалла КЭР в указанном выше предположении о структуре ростовых единиц.
Рис. 4. Зависимости нормальной скорости роста V от отклонения от равновесия А^/кТ, построенные по результатам моделирования роста некоторых граней одноатомного неполярного кристалла с ГЦК-структурой
Рис. 6. Полученная в результате моделирования зависимость скорости движения ступени вдоль направления [001] от отклонения от равновесия А^/кТ для грани (100) кристалла КЭР (черные кружки, скорость выражена в относительных единицах) в сравнении с зависимостью скорости ступени от пересыщения, полученной экспериментально [17] (белые кружки, скорость выражена в см/с)
Полученные зависимости скорости движения ступени вдоль направления [001] от отклонения от равновесия Дц/kT качественно совпадают с экспериментальными кривыми скорости движения ступени в зависимости от пересыщения раствора (рис. 6). Значительный разброс результатов моделирования объясняется вероятностным характером используемого метода Монте-Карло.
Заключение
Предложенная модель может найти применение для исследования роста различных граней кристаллов, в том числе и с дефектами структуры. При наличии подробной информации об энергиях взаимодействия ростовых единиц разработанный алгоритм может использоваться для определения конфигурации вици-нальных и дислокационных холмиков, оценки плотности атомов в ростовой ступени и на ее торце, для исследования взаимодействия поверхности и примеси (отравление изломов, адсорбция и т. п.), а также для изучения формирования макроступеней.
Список литературы
1. Landau D.P., Binder K.A Guide to Monte Carlo simulations in statistical physics. Cambridge, 2005. 435 p.
2. Fishman G. Monte Carlo: concepts, algorithms and applications. Springer-Verlag, 1996. 698 p.
3. Методы Монте-Карло в статистической физике / Ред. К. Биндер. М.: Мир, 1982. 400 с.
4. Piana S., Gale J.D. // Journal of Crystal Growth. 200б. V. 294. P. 4б-52.
5. Черепанова Т.А. // Рост кристаллов. 1980. Т. 13. С. 143-153.
6. Gilmer G.H. // Journal of Crystal Growth. 1977. V. 42. P. 3-10.
7. Александров Л. Н. Моделирование роста и легирования полупроводниковых пленок методом Монте-Карло. Новосибирск: Наука, 1991. 125 с.
8. Зверев А.В., Неизвестный И.Г., Шварц Н.Л. и др. // Физика и техника полупроводников. 2001. Т. 35. Вып. 9. С. 10б7-1074.
9. Smolsky I.L., Voloshin A.E., Zaitseva N.P. et al. // Philosophical Transactions of the Royal Society of London. 1999. V. 357. P. 2б31-2б49.
10. Tadasuke K., Tetsuo I., ShinichiroY. et al. // Journal Cryst. Growth. 2002. V. 243. N 3-4. P. 517-521.
11. Кристаллографический банк данных неорганических структур ICSD (Карлсруэ, Германия). URL: http://icsd.fiz-Karlsruhe.de/ (дата обращения 20.07.2009).
12. Кембриджский банк структурных данных CCDC (Кембридж, Великобритания). URL: http:// www.ccdc.cam.ac.uk/ (дата обращения 20.07.2009).
13. Чупрунов Е.В., Хохлов А.Ф., Фаддеев М.А. Кристаллография. М.: Физматлит, 2000. 49б с.
14. Павлов П. В., Хохлов А. Ф. Физика твердого тела. Н. Новгород: ННГУ, 1993. 490 с.
15. Гуськов С.С., Фаддеев М.А., Чупрунов Е.В. // Вестник Нижегородского университета им. Н. И. Лобачевского. 2008. № б. С. 39-45.
16. Современная кристаллография. Т. 3. Образование кристаллов. / Ред. А. А. Чернов. М.: Наука, 1980. 408 с.
17. Rashkovich L.N., Kronsky N.V. // Journal of Crystal Growth. 1997. V. 182. P. 434-441.
STATISTICAL MODELLING OF CRYSTAL GROWTH WITH THE ACCOUNT OF ITS ATOMIC STRUCTURE
S.S. Guskov, M.A. Faddeev, E. V. Chuprunov
An original numerical algorithm to simulate the face growth of a crystal with an arbitrary atomic structure has been proposed. The crystal may have dot defects and dislocations. The algorithm allows for introduction in the growth system of an impurity replacing certain sorts of atoms. The model can be used to study the formation mechanisms of the monocrystal oscillatory zoning pattern which are associated with the impurity non-uniform distribution.
Keywords: crystal growth, atomic structure, interatomic potential, defects, oscillatory zoning, impurity, mathematical modelling.
