Атомистическое компьютерное моделирование кристаллической структуры и свойств минералов, их дефектов и твердых растворов

Глава 2. Разработка и усовершенствование методов атомистического компьютерного моделирования.

2.1. Метод минимизации энергии атомизации, основанный на поиске оптимальных зарядов на атомах.

Проблема учета энергии переноса заряда при ионизации катионов, содержащих в валентной электронной оболочке электроны разных типов, потребовала разработки определенной физически обоснованной процедуры смешения энергий отдельных валентных уровней. Такая процедура учета энергии переноса нецелочисленного заряда катиона, находящегося в гибридном валентном состоянии, при любом значении f разрабатывались при непосредственном участии автора на протяжении ряда лет [Urusov V.S. and Eremin N.N. (1995), (1996), (1997), Еремин Н.Н (1996), Урусов В.С., Еремин Н.Н., Оганов А.Р. (1997), Урусов В.С., Оганов А.Р., Еремин Н.Н. (1999)]. Основным принципом процедуры является метод выравнивания орбитальных электроотрицательностей (EEM) отдельных электронных уровней одного ионизируемого атома. Под орбитальной электроотрицательностью понимается первая производная энергии орбитали по ее заселенности: χ_i = dE_i / dq_i, где q_i - число электронов. Принцип выравнивания сводится к требованию: χ_s = χ_p = χ_d = χ_f, что является основанием для определения заселенности s, p, d, f-орбиталей при заданной величине , т.е. для восстановления хода ионизации. На рис. 2. в качестве примера приведена синтезированная кривая ионизации атома Cr в spd-состоянии.

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

2.2. Разработка методик структурного моделирования твердых растворов минералов.

Теоретические подходы к моделированию структур и свойств концентрированных твердых растворов основываются на феноменологических моделях, либо на атомистических расчетах или расчетах . Для учета всех локальных конфигураций, которые возникают в изоморфной смеси требуется применение моделей, принципиально отличных от методики вложенных сфер, используемой для расчетов изолированных точечных дефектов [Mott N. F., Littleton M.J. (1938)]. На рис. 4 приведена общая схема оригинальной методики моделирования концентрированных твердых растворов замещения во всем диапазоне составов твердого раствора.

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

На втором этапе выбирается сверхструктура, в пределах которой задаются конфигурации различного состава, максимально приближенные к неупорядоченному твердому раствору. Описаны проблемы, связанные с выбором неупорядоченной конфигурации, и предлагается новый способ их решения, подробно описанный в недавней работе автора [Еремин Н.Н., Деянов Р.З., Урусов В.С. (2008)]. В качестве критерия степени неупорядоченности конфигурации использовалась величина квадратов отклонений числа разнородных пар атомов во второй координационной сфере для случайной конфигурации от статистической теоретической гистограммы (критерий согласия Пирсона χ²). Множество случайных конфигураций анализировалось по величине отклонения от идеальной статистической гистограммы частоты встречаемости разнородных вторых соседей для каждого состава (рис. 5). Результирующая гистограмма распределения (серые столбцы на рис. 5) получается как суперпозиция двух и более вкладов от атомов различных сортов, изоморфно замещающих друг друга в одной атомной позиции. К достоинствам предложенного способа нахождения оптимальной конфигурации атомов, реализованного в программе BINAR [Еремин Н.Н., Деянов Р.З, Урусов В.С. (2008)] можно отнести: наличие количественной оценки качества конфигурации, возможность в рамках сравнительно небольшой сверхячейки максимально приблизиться к неупорядоченному статистическому распределению в кристалле, возможность в дальнейших расчетах обойтись одной конфигурацией, что немаловажно для минимизации расчетного времени.

Определение числа k+1 интервалов разбиения гистограммы (задание рационального радиуса второй координационной сферы для позиции замещения) иногда становится достаточно сложным и неоднозначным. Если для ряда простых структур, построенных по законам плотнейших шаровых упаковок число вторых соседей равно в точности 12, координационный многогранник - правильный кубоктаэдр, то для менее симметричных структур могут возникнуть неоднозначности. При выборе второго координационного числа предлагается руководствоваться следующими правилами:
а) вторые соседи должны окружать центральный атом как можно более равномерно со всех сторон (образовывать координационный полиэдр таким образом, чтобы центральный атом находился внутри него и вблизи от его <центра тяжести> с минимальной несферичностью);
б) число вторых соседей должно находиться, как правило, в интервале 10-14 (12 - для ГЦК и ГПУ и 14 - для ОЦК-упаковок). Отметим, что число 10 является минимально рекомендованной величиной для критерия Пирсона χ².

Так, для структурного типа рутила TiO₂ при построении замкнутого полиэдра вокруг центрального атома необходимо учитывать, помимо ближайших, 8 следующих вторых соседей. Следовательно, общее число вторых соседей катиона равно 10 (2+8) (рис. 6а). В структурном типе корунда ситуация более сложная, и, как показал проведенный геометрический анализ, наиболее обоснованное число вторых соседей катиона равно 13 (рис. 6б).

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

Для практической реализации предложенного метода моделирования твердых растворов замещения автором совместно с Деяновым Р.З. и Урусовым В.С. был написан пакет компьютерных программ ODSS, состоящий, помимо программы BINAR, еще из трех программ: RELAX, GISTOGRAMMA и VOLUME, используемых на пятом этапе для анализа локальной структуры твердого раствора. Так как сверхъячейка в рамках вышеописанного подхода описывается симметрией P1, расположение атомов не закреплено элементами симметрии и при оптимизации они смещаются со своих первоначальных позиций. Смещения атомов характеризуют релаксацию структуры твердого раствора. Автором предложена величина <сдвигаемость атома> ω_i, определяемая как усредненная по всем парам i-j среднеквадратичная разность их начального и конечного расстояний:

(2)

Здесь n - число частиц в сфере, внутри которой анализируются атомные смещения, R_ij и R'_ij начальное и конечное расстояние от атома i до атома j, соответственно. Помимо расчета сдвигаемостей атомов, в программе RELAX можно оценить податливость каждой катионной позиции C_s, определяемая по аналогии с определением Долласа [Dollase W.A. (1980)] как реальная доля увеличения (или уменьшения) конкретных длин связей в полиэдре относительно разности длин связи чистых компонентов:

(2)

Здесь - среднее значения межатомного расстояния для конкретного катиона со своим ближайшим окружением в твердом растворе, R_addit - аддитивное значение межатомного расстояния по правилу Вегарда (приближение виртуального кристалла), а R - расстояние в чистом беспримесном кристалле. Отметим, что понятия <сдвигаемость атома> и <податливость позиции> описывают несколько различные типы искажений в кристаллической структуре. Проиллюстрируем это на примере искажения октаэдра MO6 (рис. 7).

Так, при смещении центрального катиона из геометрического центра октаэдра может возникнуть вариант, когда три связи катион-анион удлинятся, а три укоротятся практически на одну и ту же величину ΔR. Такая картина может наблюдаться, например, в корундовых системах при смещении катиона Al (Cr, Fe и т.д.) вдоль оси c элементарной ячейки. В этом случае величина податливости позиции C_s будет либо равна нулю, либо крайне мала, что может быть ошибочно интерпретировано как отсутствие атомных смещений вообще. Напротив, значения ω в этом случае будут отличны от нуля, ω центрального катиона будет существенно превышать величины сдвигаемостей анионов, что позволит сделать вывод о том, что из стартовой позиции сдвигается именно катион.

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

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

(4)

Здесь - усредненное по всем катионам одного сорта в сверхячейке значение объема координационного полиэдра, V_addit - аддитивное значение объема полиэдра в соответствии с правилом Ретгерса, а V - объем полиэдра в чистом беспримесном кристалле. Для оценки объемов тетраэдров и октаэдров на основании алгоритма Робинсона [Robinson K. et al (1971)] автором была написана программа VOLUME, проводящая геометрический анализ (включая расчет объемов) полиэдров. Подробное описание всех программ доступно по адресу http://cryst.geol.msu.ru/odss/.

Разработанный пакет программ был использован при компьютерном моделировании ряда изоморфных систем минералов [Урусов В.С. и др. (2007), Eremin N.N., Deyanov R.Z., Urusov V.S. (2008), Урусов В.С., Петрова Т.Г., Еремин Н.Н. (2008), Еремин Н.Н., Талис Р.А., Урусов В.С. (2008)]. Конкретные результаты расчетов свойств смешения и локальной структуры твердых растворов представлены в главах 5 и 6.