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

Глава 3. Уточнение параметров межатомных потенциалов, структурное моделирование, расчет и предсказание свойств кристаллов различных минералогических классов.

В главе приведены конкретные результаты расчетов в различных приближениях кристаллических структур и свойств минералов и неорганических соединений. В 3-1 возможности и ограничения чисто ионного приближения демонстрируются двумя примерами - структурным моделированием купратов щелочноземельных металлов, обладающих сверхпроводящими свойствами, и структурной оптимизацией известных минералов - кианита, андалузита и силлиманита, являющимися полиморфными модификациями Al₂SiO₅.

3-1-1. Моделирование купратов щелочноземельных металлов.

Благодаря исследованиям высокотемпературной сверхпроводимости, открытой в 1986 году, для класса купратов в последние два десятилетия был накоплен значительный объем кристаллохимических данных [Leonyuk L.I. et al (1998)]. Однако работы по структурному моделированию этих соединений были крайне малочисленны, что связано с тем, что CuO₄ квадратные группировки практически не поддаются корректному описанию с использованием традиционных приемов структурного моделирования. В работе в рамках модели формальных зарядов был разработан набор потенциалов, способный корректно описать кристаллические структуры купратов Ca, Sr и Ba. Проведенное структурное моделирование Ca₂CuO₃, CaCuO₂, SrCuO₃, (Sr_1.19Ca_0.73)Cu₂O₄ и BaCuO₂ позволило предсказать ряд их неизвестных энергетических и упругих свойств. Разработанный набор потенциалов оказался хорошо переносимым в схожие структуры, несмотря на различия в мотивах соединения группировок Cu^-O₄. Впоследствии он был рекомендован авторами программы GULP для структурного моделирования купратов, содержащих CuO₄ квадратные группировки [Gale J.D., Rohl A.L. (2003)].

3-1-2. Система кианит-андалузит-силлиманит.

На примере полиморфных модификаций Al₂SiO₅ (кианита, андалузита, силлиманита и двух возможных фаз высокого давления) было показано, как с помощью ионной модели можно относительно успешно предсказывать структуру, упругие свойства минералов, уравнение состояния, фононные спектры и термодинамических свойства кристаллических веществ. Однако, показан целый ряд систематических ошибок ионной модели, связанный с завышением вклада электростатического взаимодействия. В частности, показаны заметные искажения структуры силлиманита и гипотетической фазы высокого давления со структурой V₃O₅. Для упругих свойств кристаллов получены хорошие (в целом несколько завышенные) полуколичественные оценки, а для диэлектрических свойств - качественные. С удовлетворительной точностью были рассчитаны фононные спектры, в частности, фононный спектр кианита. Вычисленные термодинамические функции неплохо согласуются с экспериментом. Вместе с тем, попытка воспроизведения фазовой диаграммы Al₂SiO₅ в рамках ионной модели оказалась неудачной. Для исправления этой и других ошибок следовало учесть частичную ковалентность химических связей, что было сделано в 3.3.2 для кианита.

В 3-2 возможности ковалентного приближения демонстрируются на разнообразных кристаллохимических объектах, которые невозможно корректно описать ионной моделью.

3-2-1. Структурное моделирование фосфатов структурного типа насикона.

Структурный тип NaZr₂(PO₄)₃ (NZP, насикон) широко распространен среди неорганических солей с тетраэдрическими анионами PO₄^3-. Благодаря своим свойствам (термической, химической и радиационной устойчивости, малому тепловому расширению, низкой теплопроводности, хорошей ионной проводимости и т.д.) они являются предметом интенсивных теоретических и экспериментальных исследований. Однако исследования по их структурному моделированию до последнего времени не проводились. В работе была разработана модель межатомных потенциалов и проведено компьютерное моделирование кристаллической структуры щелочно-циркониевых фосфатов МZr₂(PO₄)₃, где М = Li, Na, K, Rb, Cs и Fr, с использованием частично-ковалентного приближения. Был предсказан ряд их неизвестных энергетических и упругих свойств.

3-2-2. Разработка набора межатомных потенциалов для структурного моделирования дигидрофосфата калия KDP.

Дигидрофосфат калия KH₂PO₄ (KDP) является одним из важнейших сегнетоэлектриков. Простота получения из низкотемпературных водных растворов кристаллов большого размера, высокое оптическое качество и стойкость к лазерному излучению ставит KDP на одно из первых мест при использовании в квантовой электронике и лазерной технике. Для моделирования структуры KDP было выбрано частично ковалентное приближение, как соответствующее реальному характеру химической связи в соединении. Параметры парных потенциалов межатомного взаимодействия в форме Букингема и Морзе оптимизировались с использованием экспериментальной информации о кристаллической структуре и упругих свойствах соединения. Разработанная модель позволила воспроизвести параметры элементарной ячейки KDP с ошибкой не более 1.8%, а атомные координаты - не более 2.5%. Добавочным параметром оценки качества проведенного моделирования и предсказательной возможности оптимизированного набора потенциалов явилось согласие рассчитанной и экспериментальной величины энтальпии образования соединения. Это позволило применить разработанный набор потенциалов в расчетах дефектных областей структуры KDP (глава 4).

3-2-3. Особенности структурного моделирования соединений, содержащих органические комплексы на примере бифталата калия C₈H₅O₄K.

Один из известнейших нелинейно-оптических кристаллов бифталат калия C₈H₅O₄K (БФК) (рис. 8) содержит в своей кристаллической структуре как органический комплекс C₈H₅O₄, так и неорганический ион K. В связи с этим, роль кулоновского взаимодействия достаточно велика по сравнению с чисто органическими соединениями. Следовательно, использование потенциалов, разработанных для органических кристаллов, не приводит (и не может привести) к удовлетворительному воспроизведению кристаллической структуры. С другой стороны, очевидно, что атом калия в структуре БФК не несет на себе полный положительный заряд +1, а находится в некотором промежуточном состоянии, и использование известных из моделирования структур неорганических соединений потенциалов парного взаимодействия K-O также не может привести к положительному результату. В связи с этим структурное моделирование БФК осуществлялось в рамках частично ковалентного приближения. Помимо кулоновского взаимодействия, использовались потенциалы Морзе и Букингема, а также угловые трехчастичные гармонические потенциалы и потенциалы кручения. Использование разработанного оптимизированного набора потенциалов позволило превосходно воспроизвести структуру БФК: максимальные расхождения в параметрах элементарной ячейки составили всего лишь 0.23%, а атомные координаты отличались от экспериментальных значений не более чем на 0.012 Å (в среднем - 0.001 Å).

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

3.3.1. Минимизация энергии атомизации кварца, корунда и стишовита

В качестве объектов исследования для отработки методики минимизации энергии атомизации были выбраны следующие соединения кремния и алюминия: α-кварц и стишовит SiO₂ и корунд Al₂O₃. Для их структурного моделирования была выбрана модель парных потенциалов с введением квадратичных весовых функций для потенциалов Букингема и Морзе. Оказалось, что структуру кварца можно удовлетворительно воспроизвести лишь в очень узком интервале степеней ионности f=50(6)%, в то время как структура стишовита воспроизводима при степенях ионности приблизительно от 45 до 70%. Структура корунда может быть с успехом рассчитана в пределах f=60(±15)%. Оптимальные величины атомных зарядов, соответствующие минимумам на кривых зависимости энергии атомизации от степени ионности связи равны: Z(Si^IV) = +1.92e, Z(O) = -0.96e для кварца; Z(Si^VI) = +2.04e, Z(O) = -1.02e для стишовита и Z(Al)= +1.86e, Z(O) = -1.24e для корунда. Особенно убедительно преимущества модели с варьируемыми зарядами проявляются в предсказании энергии атомизации кристаллов, которая при оптимальных зарядах атомов оказывается весьма близкой к эксперименту. Неплохо предсказывается и энтальпия перехода кварц-стишовит (расчет - 0.34 эВ, эксперимент - 0.51 эВ).

3.3.2. Минимизация энергии атомизации кианита Al₂SiO₅

Одним из критериев разумности и результативности предложенной модели с варьируемыми атомными зарядами является ее предсказательная способность. Для проверки было осуществлено структурное моделирование кианита, в котором варьирование зарядов осуществлялось одновременно на атомах Si и Al. Структура и свойства кианита были рассчитаны с использованием потенциалов, разработанных для кварца и корунда. Полученные для параметров ячейки, координат атомов и энергии сцепления результаты показали превосходное согласие с экспериментом. Более того, на момент проведения расчетов экспериментально упругие свойства кианита не изучались. Однако в 1997 году в работе [Yang et al. (1997)] сжимаемость кианита была экспериментально изучена и экспериментальная величина модуля всестороннего сжатия Kэксп, равная 193(1) гПа, оказалась в великолепном согласии с предсказанной по оптимальной модели Kрасч = 191 гПа.

На рис. 9 показан ландшафт поверхности энергии атомизации кианита Al2SiO5 в координатах степеней ионности связей Al-O и Si-O. Обратим внимание на то, что оптимальные значения f(Al-O) и f(Si-O) близко совпадают с соответствующими значениями для корунда и кварца, а значение эффективного заряда кислорода (-1.15e) является промежуточным между его зарядами в корунде и кварце. Величина эффективного заряда Al (+1.83e) крайне незначительно отличается от заряда Al в корунде (+1.86e). Кроме того, нужно отметить чрезвычайно близкое совпадение рассчитанного (-50.96 эВ) и экспериментального (-51.27 эВ) значений энергии атомизации кианита.

3.3.3. Минимизация энергии атомизации малаяита CaSnOSiO₄.

В отличие от изоструктурного титанита CaTiOSiO₄, в структуре малаяита атомы олова уже при комнатной температуре находятся в центре SnO₆-октаэдров, что обуславливает кристаллизацию минерала в пространственной группе A2/a. Моделирование малаяита осуществлялось в частично ковалентном приближении методом минимизации энергии атомизации. Такое моделирование представляло самостоятельный интерес, так как упругие и диэлектрические свойства минерала экспериментально не исследовались. Проведенные расчеты позволили успешно воспроизвести структурные особенности малаяита, предсказать значения упругих, диэлектрических и энергетических характеристик соединения, которые согласуются с имеющейся ограниченной экспериментальной информацией о его свойствах. Поверхность E(Z_Sn,Z_Si) имеет четко выраженный минимум при зарядах Si=+2.40(4)e и Sn=+3.20(4)e, что согласуется с мессбауэровскими данными, согласно которым заряд на ядре олова в малаяите выше, чем в SnO₂ (для последнего расчет дает Z_Sn=+2.96e).

В качестве модельных соединений для тестирования модифицированной версии программы VALUES, позволяющей рассчитывать кривые энергии переноса заряда для атомов, имеющих во внешней валентной оболочке электроны трех различных типов, были выбраны оксиды хрома и марганца: эсколаит Cr₂O₃, пиролюзит MnO₂, курнакит Mn₂O₃ и манганозит MnO. Во всех случаях в точке зарядового минимума наблюдалось близкое согласие рассчитанных упругих свойств кристаллов с экспериментально наблюдаемыми значениями модулей всестороннего сжатия, что указывало на достаточную корректность синтезированных кривых переноса заряда для атомов Cr и Mn.

В 3.4 представлены результаты разработки согласованных потенциальных наборов для ряда оксидных и галоидных изоструктурных соединений с целью последующего изучения их твердых растворов.

3.4.1. Оксиды щелочноземельных металлов со структурой галита.

В структурном типе NaCl кристаллизуются несколько десятков природных и синтетических оксидов с формулой MO. Для них феноменологическая теория твердых растворов замещения разработана наиболее полно (Урусов, 1977). Естественно, что эти соединения являются одними из наиболее изученных объектов методом межатомных потенциалов, тем более что для них наиболее простая ионная модель может быть применена совершенно обоснованно. Однако, парные потенциалы с учетом неполной ионности химической связи гораздо лучше описывают упругие свойства кристаллов, поэтому в нашей работе предпочтение было отдано модели потенциалов со степенью ионности f=0.85.

3.4.2. Модель межионных потенциалов для галита и сильвина.

Так как связь в галогенидах носит еще более ионный характер, чем в изоструктурных им оксидах, было принято решение не вводить в межатомные потенциалы для этой системы ковалентную составляющую. Пробные расчеты показали, что отталкивание во второй координационной сфере Me-Me и Cl-Cl может быть описано исключительно электростатическими силами без использования короткодействующих потенциалов, поэтому в модели короткодействующих потенциалов были использованы только потенциалы Букингема Mе-Cl. Разработанная модель обеспечивает точное описание структуры и упругих свойств NaCl и KCl. Температурные зависимости теплоемкости хорошо согласуются с экспериментом, однако зависимости энтропии от температуры испытывают систематические отклонения от экспериментальных величин, что связано, вероятно, с чрезмерной простотой набора парных потенциалов. Тем не менее, разность SKCl-SNaCl почти в точности соответствует экспериментально наблюдаемой, что позволяет использовать предлагаемую модель ионных потенциалов для дальнейшего моделирования твердых растворов в системе галит-сильвин.

3.4.3. Разработка согласованной модели потенциалов для полуторных оксидов M₂O₃ со структурой корунда.

Несмотря на то, что корунд неоднократно подвергался структурному моделированию ранее, критический анализ предлагаемых моделей показал необходимость разработать оригинальную согласованную модель для последующего моделирования твердых растворов состава M₂O₃. Заряд металла в разработанной модели был принят равным величине +1.89e, а совместное использование потенциала Букингема и Морзе позволило существенно улучшить описание взаимодействия катион-анион. Некоторые рассчитанные структурные, упругие и термодинамические свойства корунда, гематита и эсколаита приведены в таблице 1 в сравнении с экспериментальными данными и расчетными данными из работы [Catlow C.R.A. et al (1993)].

Таблица 1. Рассчитанные свойства корунда, гематита и эсколаита в сравнении с экспериментальными величинами.
Соединение	а, Å	c, Å	V, Å3	K, ГПа
Al2O3 (эксперимент) [Gale J.D. et al (1992)]	4.760	12.993	254.977	253.60 [Goto T. et al (1989)]
Al2O3 (расчет)	4.753	13.016	254.613	296.58
Al2O3 (расчет) [Catlow C.R.A. et al (1993)]	4.7688	12.9869	255.77	333.22
Cr2O3 (эксперимент) [Ishimata N. et al (1980)]	4.954	13.566	288.356	238 [Finger L., Hazen R.M. (1980)]
Cr2O3 (расчет)	4.965	13.538	288.993	230.26
Fe2O3 (эксперимент) [Antipin M.Iu. et al (1985)])	5.050	13.720	303.017	225 [Finger L., Hazen R.M. (1980)]
Fe2O3 (расчет)	5.053	13.730	303.537	213.82

Как видно из таблицы, наша модель существенно успешнее описывает упругие свойства корунда, чем модель работы Кэтлоу, которая, более того, оказалась совершенно непригодной для моделирования термодинамических свойств соединения. С нашей точки зрения, это вызвано тем обстоятельством, что авторы вышеназванной работы выбрали несколько завышенные значения эффективных атомных зарядов +2.08 (Al) и -1.39 (O) (взятых из квантовомеханического расчета небольшого кластера [Gale J.D. et al (1992)]).

Как видно из рис. 10 модель с зарядами +1.89 (Me) и -1.26 (O) успешно описывает энтропию корунда, во всем интересующем нас диапазоне температур. Таким образом, в результате проведенных расчетов был получен согласованный набор парных потенциалов, удобный для последующего моделирования дефектов и твердых растворов замещения в системе Al₂O₃-Fe₂O₃-Сr₂O₃.

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

Методом атомистических потенциалов соединения бериллия практически не изучались, возможно, из-за отсутствия до сих пор рекомендованного потенциала межатомного взаимодействия Be-O в известных банках межатомных потенциалов. Оптимизация параметров короткодействующих потенциалов Be-O и O-Be-O в рамках частично-ковалентного приближения проводилась по экспериментальным значениям параметров ячейки, координат атомов и упругих констант бромеллита BeO, что позволило точнее воспроизвести кристаллическую структуру BeO, чем расчет ab-initio [Шеин И.Р. и др. (2007)]. Рассчитанный модуль всестороннего сжатия отличается от экспериментальной величины всего на 5% (расчет ab-initio - 21%). С целью проверки трансферабельности разработанного набора параметров потенциалов Be-O он был использован при моделировании структур и свойств хризоберилла BeAl₂O₄ и его изоструктурных аналогов: BeСr₂O₄ и неизученного до сих пор рентгеновскими методами BeFe₂O₄. Рассчитанные параметры элементарной ячейки хризобериллов отличаются от экспериментальных величин не более чем на 1%, а координаты атомов не более чем на 1.75%. В целом, хорошо воспроизвелся и тензор упругих констант BeAl₂O₄. Как видно из рис. 11 энтропия соединений согласуется с экспериментальными значениями в пределах 3% в интервале температур от 298 до 1800 К.

Рассчитанные структурные и упругие свойства минералов также хорошо согласуются с имеющейся экспериментальной информацией, что позволяет считать ряд рассчитанных свойств этих соединений и модельную кристаллическую структуру Fe-хризоберилла (табл. 2) достоверным предсказанием. Это дает основание, во-первых, использовать разработанный набор потенциалов для моделирования свойств смешения и локальной структуры твердых растворов в хризобериллах (главы 5 и 6), а во-вторых, использовать его для моделирования других соединений, содержащих BeO₄ тетраэдрические группировки, например фенакита Be₂SiO₄ и берилла Be₃Al₂[Si₆O₁₈]. Кристаллические структуры этих силикатов хорошо известны и могут быть описаны как смешанные каркасы из тетраэдров BеO₄ и SiO₄. Для их моделирования с использованием разработанного набора потенциалов потребовалось оптимизировать потенциал взаимодействия Si-O для заряда на атоме Si +2.52e (для соблюдения электронейтральности элементарных ячеек). Потенциалы Si-O в форме Букингема и Si-O-O в трехчастичной гармонической форме оптимизировались на кристаллической структуре α-кварца.

Отметим, что предложенная зарядовая модель не противоречит значениям эффективных зарядов в этих соединениях, полученных из ab-initio расчетов и из различных моделей обработки прецизионных рентгеновских данных. Как показали проведенные расчеты, вышеописанный набор потенциалов отлично воспроизвел кристаллические структуру берилла и фенакита (отклонения параметров элементарных ячеек менее 1%) при разумном описании его упругих и термодинамических свойств. Для фенакита модуль всестороннего сжатия расходится с экспериментом не более 16%, а энтропия - не более 3% (T= 300-1800 К). Для берилла модуль всестороннего сжатия отличается от экспериментальной величины на 15%, а отклонения энтропии от экспериментальных значений не превышают 5% при температурах выше 900 К, и лишь при более низких температурах значения энтропии занижены по сравнению с экспериментальными значениями.

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

Рутил и касситерит неоднократно и с переменным успехом подвергались структурному моделированию методом атомных потенциалов. Проведенный в работе анализ предыдущих исследований продемонстрировал необходимость разработки самостоятельной согласованной модели потенциалов. Наилучшее согласие экспериментальных и расчетных величин структурных, упругих и термодинамических характеристик этих оксидов удалось достичь при зарядах Ti и Sn равным +2.8, а кислорода - -1.4. Как видно из таблицы 3, разработанная для рутила и касситерита согласованная модель межатомных потенциалов неплохо описывает их структурные, упругие (а также термодинамические) свойства и может быть использована для проведения дальнейших расчетов в изоморфной системе (Sn,Ti)O₂.

Таблица 3. Сравнение вычисленных и экспериментальных [Park et al. (1985), Ahrens (1995)] структурных и упругих свойств TiO2 и SnO2.
Параметр	TiO2		SnO2
	эксперимент	расчет	эксперимент	расчет
a, с (Å)	4.5943, 2.9586	4.5993, 2.9857	4.737, 3.187	4.7623, 3.1591
х (O)	0.305	0.3084	0.306	0.3059
V, Å3	62.449	63.158	71.559	71.645
C11, C12 (ГПа)	26.9, 17.7	32.0, 10.8	26.3, 17.8	31.7, 16.6
C13, C33 (ГПа)	14.8, 48.2	18.6, 56.4	15.6, 45.0	14.4, 49.2
C44, C66 (ГПа)	12.4, 19.3	13.7, 15.9	10.3, 20.8	13.5, 17.6
K, ГПа	218.6	226.6	217.7	223.2