Межатомные взаимодействия в конденсированных средах

1.1.Квантовомеханическое обоснование понятия межатомного потенциала

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

. (1.1)

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

, (1.2)

где ‑ полная волновая функция системы, являющаяся функцией координат ядер и электронов (здесь и далее индексы, обозначенные латинскими буквами, пробегают номера всех ядер, а греческие индексы – номера всех электронов). Оператор Гамильтона системы H равен

, (1.3)

где ‑ оператор кинетической энергии ядер ( ‑ масса ядер), ‑ оператор кинетической энергии электронов ( ‑ масса электрона), а

(1.4)

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

, (1.5)

. (1.6)

Уравнение Шредингера с невозмущенным гамильтонианом ,

(1.7)

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

Рассматривая теперь уравнение Шредингера для возмущенной системы,

, (1.2а)

в соответствии с теорией возмущений решение следует искать в виде разложения по собственным функциям оператора :

. (1.8)

Подставляя (1.8) в (1.2а), умножая на комплексно-сопряженную волновую функцию и интегрируя по координатам электронов, получим систему уравнений

, (1.9)

где оператор имеет вид

, (1.10)

и ‑ элементарный объем в пространстве координат всех электронов.

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

(1.11)

для каждого состояния движения электронной подсистемы, определяемого квантовыми числами m.

Из уравнений (1.11) видно, что движение ядер характеризуется потенциальной энергией , которая соответствует энергии электронов в уравнении (1.7) при фиксированных координатах ядер. Эта энергия при переходе к классическому описанию движения ядер и соответствует потенциальной энергии межатомного взаимодействия . Таким образом, понятие межатомного потенциала находит квантовомеханическое обоснование в адиабатическом приближении.

1.2. Парные потенциалы

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

, (1.12)

где ‑ расстояние между атомами пары i-j.

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

Наиболее распространенными парными потенциалами являются потенциалы Леннарда-Джонса и Морзе.

Потенциал Леннарда-Джонса (ЛД) имеет функциональный вид

, (1.13)

где s и e ‑ постоянные коэффициенты, имеющие размерности длины и энергии, соответственно.

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

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

Потенциал Морзе имеет форму

. (1.14)

Наличие трех параметров e, a и r ₀ позволяет подгонять потенциал, вдобавок к энергии связи и параметру решетки, к модулю всестороннего сжатия B кристалла, связанному с крутизной потенциала около минимума. Этот потенциал часто использовался для моделирования металлов с г.ц.к. и г.п.у. решетками.

Более сложные парные потенциалы, имеющие больше параметров и, соответственно, подогнанные к большему количеству свойств веществ, были получены путем “сшивки” нескольких полиномиальных функций (так называемые сплайн-потенциалы). С помощью парных потенциалов в 1970-80-х гг. были получены многие фундаментальные результаты о свойствах точечных, линейных и плоских дефектов в кристаллах.

Для составления и решения уравнений движения атомов необходимо рассчитывать силы, действующие на них. Это делается на основании уравнения (1.1). Для случая парных потенциалов силы рассчитываются следующим образом.

(1.15)

Для потенциала ЛД

. (1.16)

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

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

1. Число пар взаимодействующих атомов растет с числом атомов N как . Для системы из N =100 атомов имеется 5000 пар, а при N =3000 ‑ 4500000. При еще большем увеличении числа атомов расчет взаимодействия всех атомов становится невыполнимой задачей. В этом нет и необходимости, так как силы взаимодействия атомов достаточно быстро убывают (за исключением ионов в ионных кристаллах). Использование радиуса обрезания 8‑10 Å сокращает число соседей, взаимодействующих с каждым атомом, до примерно 50, поэтому при N =3000 силы взаимодействия рассчитываются только примерно 50 N =150000 раз.

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

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

. (1.17)

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

Например, в работе Стоддарда и Форда (Phys. Rev.A, 1973, v.8, p. 1504) использован обрезанный потенциал Леннарда-Джонса следующего вида:

. (1.18)

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

Парные потенциалы имеют существенные недостатки.

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

Во-вторых, парные потенциалы не учитывают зависимость силы связи от направления, которая ярко выражена в ковалентных кристаллах, а также в переходных металлах, в которых в связь дают вклад d -электроны. Парные потенциалы работают лучше для металлов, в которых связь осуществляется s - и p -электронами.

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

В-четвертых, парные потенциалы неправильно описывают упругие свойства кристаллов. При их использовании между упругими константами кубических кристаллов существует так называемое соотношение Коши , которое для металлов никогда не соблюдается (например, 1,5 для Cu; 1,9 для Ag и 3,7 для Au). Рассмотрим этот вопрос более подробно.

Плотность энергии парного взаимодействия системы N атомов равна

, (1.19)

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

, (1.20)

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

, (1.21)

где

‑ (1.22)

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

Тензор упругих постоянных определяется выражением

(1.23)

Для парных потенциалов

, (1.24)

(1.25)

. (1.26)

В силу того, что в равновесии , второе слагаемое в фигурных скобках равно нулю. Оставшийся член полностью симметричен по отношению к перестановкам индексов a, b, g и d. Отсюда следует, что тензор упругих постоянных симметричен по отношению к перестановке всех индексов. С использованием обозначений Фогта (, , , и т.д., получаем соотношение Коши или для кубических кристаллов, где и .

Современные потенциалы, учитывающие многочастичный характер взаимодействия атомов, строятся из квантовой теории электронной структуры кристаллов, или, как говорят, из первых принципов (ab initio).Эта теория основана на теории функционала плотности многоэлектронных систем, основные элементы которой мы сейчас рассмотрим.

1.3. Основы теории функционала плотности многоэлектронных систем

Как мы видели из рассмотрения адиабатического приближения, состояние системы электронов в твердом теле определяется волновой функцией , параметрически зависящей от координат ядер и зависящей от координат всех электронов . Определение этой ВФ для кристалла, содержащего 10²³ электронов в 1 см³, является безнадежной задачей. Путь описания состояний многоэлектронных систем был найден Хоэнбергом и Коном в 1964 г. (за это открытие Кон был удостоен Нобелевской премии по химии 1998 г.), которые разработали основы теории функционала плотности (Density functional theory, DFT).

Хоэнберг и Кон доказали теорему, утверждающую, что:

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

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

Этот функционал имеет вид

. (1.27)

Здесь ‑ функционал кинетической, обменной и корреляционной энергии электронов.

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

, (1.28)

где суммирование по i и j производится по ядрам твердого тела, и ‑ заряд и радиус-вектор i -го ядра, интегралы берутся по координатам электронов, ‑ суммарная энергия изолированных атомов.

1.4. Межатомные потенциалы, основанные на теории функционала плотности. Метод внедренного атома

Начиная с 1980-х гг. был предложен ряд методов описания межатомных взаимодействий, в основе которых, по сути, лежит теория функционала плотности. Изначально эти методы были обоснованы исходя из различных физических соображений, но все они приводят к похожим выражениям для энергии твердого тела. Эти методы следующие: метод погруженного (или внедренного) атома МПА (Embedded atom method, EAM), метод эффективной среды (Effective medium theory, EMT), потенциал Финниса-Синклера, клеевая модель Эрколесси и т.д. Во всех этих моделях полная потенциальная энергия твердого тела, состоящего из N атомов, представляется в виде

, (1.29)

где

(1.30)

Здесь ‑ определенные функции, конкретная интерпретация которых зависит от метода. Например, в методе эффективной среды и МПА плотность электронов в узле с номером i, энергия внедрения атома i в эту электронную плотность, энергия парного взаимодействия между атомами i и j, причем плотность электронов в каждом узле является суперпозицией электронных плотностей отдельных атомов, то есть представляет собой электронную плотность, создаваемую одним атомом. Энергия парного взаимодействия представляет собой положительную энергию отталкивания ионных остовов атомов, которые заряжены положительно. Составляющие потенциальной энергии в МПА схематически изображены на рис.1.

Рис.1.Схематическое изображение составляющих потенциала в методе погруженного атома

Конкретная функциональная форма и параметры потенциалов вида (1.29)-(1.30) могут быть определены двумя способами. Первый путь – это расчеты из первых принципов, то есть из теории функционала плотности. Однако расчеты такого типа пока не позволяют строить достаточно точные потенциалы. Второй путь – это подбор функциональной формы и параметров таким образом, чтобы получить некоторый набор свойств кристаллов, совпадающий с экспериментальными значениями. Потенциалы МПА подгоняются под следующие свойства: параметр(ы) решетки, энергию связи кристалла , упругие константы (для г.ц.к. металлов , и ) и энергию образования вакансий .

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

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

В качестве примера рассмотрим аналитический потенциал МПА, разработанный Джонсоном (Phys. Rev., 1988). Джонсоном был предложен следующий простой функциональный вид плотности электронов, парного потенциала и энергии внедрения:

, (1.31)

, (1.32)

. (1.33)

В этих выражениях , где B ‑ модуль всестороннего сжатия, ‑ атомный объем. Кроме того, и . Имеются параметры , которые определяются подгонкой к характеристикам материала ‑ параметру решетки или атомному объему , энергии связи , энергии образования вакансии , модулю всестороннего сжатия B и модулю сдвига G.

Ниже приводятся графики трех функций, подогнанных Джонсоном для никеля.

Как видно из рисунков, электронная плотность представляет собой положительную убывающую функцию, а энергия внедрения – функцию с положительной кривизной. Поэтому взаимодействие становится более отталкивательным при увеличении плотности. Это приводит, в частности к следующему эффекту. В кристалле Ni равновесное расстояние между атомами равно =2,49 Å, тогда как в димере Ni-Ni оно значительно меньше: 2,07 Å (рис. 3).

Рис. 2. Графики функций, характеризующих потенциал МПА для Ni

а б

Рис. 3. Зависимости энергии связи на один атом от параметра решетки в кристалле (а) и расстояния между атомами в димере (б) никеля

1.5. Межатомные потенциалы для сплавов

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

, (1.34)

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

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

Например, в реализации МПА, предложенном авторами метода До, Басксом и Фойлзем (Daw, Baskes, Foiles), эта часть рассчитывается следующим образом:

, (1.35)

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

, (1.36)

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

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

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

1.6. Потенциалы для ковалентных кристаллов

Элементы четвертой группы (C, Ge, Si) образуют решетку алмаза с четырьмя ковалентными связями на каждый атом, которые направлены к вершинам тетраэдра и образуют угол 109,5° друг с другом. Главным требованием к потенциалам для этих элементов является описание устойчивости этой решетки.

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

, (1.37)

. (1.38)

Первый множитель в выражении для парной части потенциала (1.38) имеет форму потенциала ЛД, а второй, экспоненциальный, представляет собой “обрезающий” член, который плавно сводит потенциал к нулю на расстоянии . Трехчастичный член потенциала имеет вид

, (1.39)

где ‑ угол с вершиной на атоме i и h – следующая функция с тремя параметрами :

. (1.40)

При зависящий от угла член (1.39) проводит к энергетически выгодной кристаллической структуре алмаза с тетраэдрическими углами для Si.

Потенциал Стиллингера-Вебера завоевал большую популярность при моделировании кремния и был использован во многих исследованиях. Этому способствовали относительная простота потенциала и довольно реалистическое описание кристаллического кремния. Однако этот потенциал имеет несколько серьезных недостатков.

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

Поэтому для элементов IV группы были разработаны так называемые потенциалы кратной связи (bond order potentials), которые основаны на учете зависимости прочности связи от локального окружения (потенциалы Терсоффа-Абеля для Si и Ge, потенциал Бреннера для C). Эти потенциалы являются в настоящее время наиболее употребительными при моделировании ковалентных кристаллов.

Дата добавления: 2015-10-02