Так как размеры наноструктур сопоставимы с межатомными расстояниями, то расположение атомов в них подвержено значительному влиянию поверхности, границ раздела и внутренних точечных дефектов. Экспериментальное определение связанной с этим релаксации кристаллической решетки крайне затруднено, поскольку требует использования методов анализа с атомным разрешением. Поэтому эффективным методом исследования и описания атомных конфигураций наноразмерных структур является их теоретическое моделирование с использованием современной вычислительной техники.
Поиск оптимального расположения атомов в наноразмерной структуре с учетом их взаимодействия друг с другом и влияния на них внешних сил осуществляется путем нахождения такого взаимного расположения, которое соответствует минимальной энергии всей структуры. Для решения этой задачи могут быть использованы два принципиально различных подхода. В их основе лежат представления о твердом теле соответственно как о механической и о квантово-химической системах. В первом случае используются методы молекулярной динамики (molecular dynamics) и молекулярной механики (molecular mechanics), во втором – методы квантовой химии. В данном разделе мы ограничимся рассмотрением теоретических основ лишь первой группы методов. Квантово-химические подходы будут описаны в последующих разделах, где мы рассмотрим фундаментальные электронные свойства наноразмерных структур. Именно для моделирования электронных свойств твердых тел они главным образом и разрабатывались.
В методе молекулярной динамики поиск оптимальной атомной конфигурации осуществляется путем пошагового движения атомов в положения, обеспечивающие минимум потенциальной энергии анализируемой структуры. Рассматривается конечная область кристалла – кристаллит. Движение атомов описывается уравнениями классической механики. Взаимодействие атомов с поверхностью или границей раздела имитируется с помощью направленных внутрь рассматриваемого объема постоянных сил, действующих только на атомы в этой граничной области. Расчет атомной конфигурации начинают с задания предполагаемых координат атомов. Затем, меняя положение атомов, находят конфигурацию, соответствующую минимуму потенциальной энергии. При этом решают систему связанных уравнений движения для каждого атома,
|
|
(2.2) |
|
|
|
(2.3) |
|
;
где ri(t) – координата i-го атома; vi(t) – его скорость; N – общее количество атомов; M – масса атома.
Сила Fi, действующая на i-й атом, является результатом взаимодействия этого атома с его соседями. Она определяется потенциалом многочастичного взаимодействия. Однако если последний неизвестен, используют потенциалы парного взаимодействия φij между i-м атомом и всеми другими j-ми атомами:
|
|
(2.4) |
,где сила, действующая на, i-й атом со стороны j-го атома, равна:
|
|
(2.5) |
.Здесь rij = ri – rj. Потенциал парного взаимодействия (φij) может быть задан в форме Ленарда-Джонса, Борна-Мейера, Томаса-Ферми, экранированного кулоновского потенциала или любой другой форме, адекватно описывающей особенности межатомного
взаимодействия в моделируемой системе. Следует иметь в виду, что полный потенциал межатомного взаимодействия очень редко сводится к парным взаимодействиям, по крайней мере, в твердых телах. Поэтому испол
ьзование потенциалов многочастичного взаимодействия дает более точный результат.
Решение системы уравнений (2. 2), (2. 3) с учетом уравнений (2.4) и (2.5) возможно только путем численного интегрирования. При этом следует учитывать следующее:
1) метод интегрирования нужно выбирать с учетом размера моделируемой области, желательной точности расчетов и вычислительной мощности используемого компьютера;
2) для получения объективного результата важно выбрать корректное описание потенциала межатомного взаимодействия в конкретной моделируемой структуре;
3) учет тепловых колебаний атомов позволяет оценить динамическую стабильность найденных равновесных атомных конфигураций;
4) «сшивка» границ выделенного кристаллита с остальным объемом должна быть корректной. Для описания гетерофазных границ необходимо использовать периодические граничные условия в соответствующих направлениях;
5) в процессе моделирования необходимо уменьшать скорости движения всех атомов на каждой итерации. В противном случае атомы будут просто колебаться, но не релаксировать к положениям равновесия.
Процедура динамической релаксации моделируемого ансамбля атомов контролируется межатомными потенциалами и законами ньютоновской механики, что дает возможность отыскать такое расположение атомов, которое соответствует состоянию с наименьшей потенциальной энергией.
B основе метода молекулярной механики лежит процедура статической релаксации. Анализируется полная потенциальная энергия системы, которая рассчитывается для различных атомных конфигураций. Конфигурация с наименьшей потенциальной энергией считается равновесной.
Если вычислительные возможности используемого компьютера не позволяют моделировать поведение структуры в целом, то для анализа выделяется часть исследуемой структуры – кристаллит. Атомы внутри кристаллита рассматриваются как взаимодействующие друг с другом частицы. Между ними задается центральное парное взаимодействие, характеризуемое потенциалом φij. При этом учитываются соседи из ближайших координационных сфер. Как отмечено ранее, описание атомной структуры твердых тел с использованием парного потенциала является приближением к реальному многочастичному взаимодействию. Кроме этого, соседи только из ближайших координационных сфер учитываются «как правило», но не всегда.
Точность предсказания равновесной атомной конфигурации методом молекулярной механики во многом определяется интуицией и опытом исследователя, поскольку только выбор оптимального алгоритма изменения координат атомов позволяет за разумное время найти их расположение, отвечающее глобальному (а не локальному) минимуму потенциальной энергии. Расчет потенциальной энергии может быть осуществлен как в рамках классической механики, так и с использованием квантово-механических методов. При этом применение процедуры Монте-Карло делает результат статистически более надежным.
Одной из важнейших проблем практического применения методов молекулярной динамики и механики является ограничение на количество атомов в анализируемой структуре. Оно связано с ограниченной вычислительной мощностью используемого компьютера. Анализ взаимодействия N атомов со своими N – 1 соседями требует нахождения N(N – 1)/2 парных взаимодействий. Если не используется никакой специальный алгоритм, то время, необходимое для расчетов, пропорционально N2. Это не яв
ляется большой проблемой для атомных кластеров (нульмерных структур), включающих до нескольких тысяч атомов. Но в случае квантовых шнуров и квантовых пленок мощность используемого компьютера будет ограничивать максимальный размер структуры, для которой проводится моделирование.
Существуют приемы, позволяющие преодолевать указанные ограничения на размер моделируемых структур. Один из них заключается в использовании периодических граничных условий. При этом из структуры «вырезается» представительная область, состояния на границах которой описываются периодическими граничными условиями. Это эквивалентно рассмотрению бесконечного пространства, заполненного идентичными копиями моделируемой области. Такая периодичность имеет два важных следствия. Во-первых, любой атом, покидающий моделируемую область через границу, неизбежно входит в нее же с противоположной стороны. Во-вторых, атомы, расположенные на границе, взаимодействуют с атомами в соседней копии моделируемой области, что эквив
алентно взаимодействию между атомами, находящимися у противоположных сторон самой моделируемой области – так называемый циклический эффект. Этот эффект необходимо учитывать при расчете межатомных взаимодействий и численном интегрировании уравнений движения.
Другой подход предполагает замену внешней (пассивной) области моделируемой структуры упругой средой с соответствующими интегральными свойствами. Деформация упругой среды и находящихся в ней атомов описывается сферически симметричным векторным полем:
|
|
(2.6) |
,где C – параметр, характеризующий среду.
В таком приближении энергия кристаллита (E) равна:
|
|
(2.7) |
где член aC описывает работу против сил, удерживающих совершенную решетку в равновесии, а слагаемое bC2 представляет собой энергию, запасённую упругим полем. Параметры С, a и b определяются из условия равенства нулю результирующей силы (FΣ), действующей со стороны упругой среды на область внутри нее:
|
|
(2.8) |
Здесь индексы i и k относятся, соответственно, к атомам вне и внутри среды.
После введения и описания внешней упругой среды к внутренней выделенной области можно применить процедуру динамической или статической релаксации. Рассмотренные приемы следует использовать с известной долей осторожности, поскольку они не полностью исключают эффекты, связанные с конечными размерами моделируемой структуры.
Сравнивая процедуры динамической и статической релаксации, следует подчеркнуть, что динамическая модель дает возможность атомам практически сразу занимать положения, отвечающие наименьшей потенциальной энергии. Процедура же статической релаксации позволяет определить, какая атомная конфигурация имеет наименьшую энергию (а следовательно, является наиболее устойчивой), только после вычисления энергии различных возможных метастабильных конфигураций. Прогресс в использовании обоих подходов неразрывно связан с повышением вычислительной мощности (быстродействия и объема оперативной памяти) используемых для этих целей компьютеров.
