CC BY
2Рос. хим. ж. (Ж. Рос. хим. об-ва им. Д.И. Менделеева), 2001, т. XLV, № 3
УДК (546.212 + 548.562): 541.68 + 54-171
Модули упругости и границы стабильности льдов и клатратных гидратов кубической структуры I
В. Р. Белослудов, Т. М. Инербаев, В. П. Шпаков, Д. С. Це, Р. В. Белослудов, Е. Кавазое
ВЛАДИМИР РОМАНОВИЧ БЕЛОСЛУДОВ — доктор физико-математических наук, старший научный сотрудник, заведующий лабораторией Института неорганической химии СО РАН. Область научных интересов: статистическая физика, исследование материалов, фазовые переходы и механические и термодинамические свойства кристаллических и аморфных веществ.
ТАЛГАТ МУРАТОВИЧ ИНЕРБАЕВ — младший научный сотрудник Института неорганической химии СО РАН. Область научных интересов: построение фазовых диаграмм льдов и клатратных гидратов.
ВЛАДИМИР ПЕТРОВИЧ ШПАКОВ — кандидат физико-математических наук Института неорганической химии СО РАН. Область научных интересов: статистическая физика, исследование материалов, фазовые
переходы и механические и термодинамические свойства кристаллических и аморфных веществ.
ДЖОН С. ЦЕ — профессор, заведующий лабораторией Стишевского института молекулярных исследований (Государственный совет Канады по научным исследованиям). Область научных интересов: квантовая механика и статистическая физика, исследование материалов, фазовые переходы и механические и термодинамические свойства кристаллических и аморфных веществ.
РОДИОН ВЛАДИМИРОВИЧ БЕЛОСЛУДОВ — кандидат физико-математических наук, младший научный сотрудник Института неорганической химии СО РАН. Область научных интересов: квантовая механика и статистическая физика, исследование клатратных соединений, дизайн новых материалов.
ЕОШИЮКИ КАВАЗОЕ — профессор, заведующий лабораторией Института исследования материалов (Тохоку университет). Область научных интересов: квантовая механика и статистическая физика, исследование материалов, дизайн новых материалов.
630090 Новосибирск, просп. акад. Лаврентьева 3, ИНХ СО РАН, тел. (8-3832)34-30-57, факс (8-3832)34-44-89, E-mail [email protected]
100 Sussex Drive, Ottawa, Ontario Canada K1A 0R6, тел. +1 613-991-1237, факс +1 613-947-2838, E-mail John. Tse@nrc. ca
Institute for Materials Research, Tohoku University, Sendai 980-77, Japan, тел. 09023621516, факс +81-22-215-2052, [email protected]
Экспериментальные исследования обнаруживают сложный характер фазовых переходов во льдах [1, 2] и гидратах [3]. При низких температурах переход «кристалл—аморфное тело», вызываемый давлением, не может быть описан как обычный процесс термодинамического плавления при повышении давления, поскольку не происходит смены фаз при переходе линии равенства термодинамических потенциалов кристаллической и аморфной фаз. Исследования кривой плавления кристаллической модификации льда Iм [4] показывают, что термодинамическое плавление должно начинаться при гораздо более низких давлениях, чем наблюдается в эксперименте.
В ряде работ [5—8] высказано предположение, что механизм аморфизации кристаллической фазы под давлением связан с динамической неустойчивостью решетки кристалла. Эта неустойчивость проявляется в наличии в спектре упругих колебаний кристаллической решетки мнимых частот, когда давление или температу-
ра превышает критическое значение [6—8]. Другой возможный механизм фазового перехода может быть связан с механической неустойчивостью твердой фазы по отношению к малой однородной деформации решетки кристалла при постоянной температуре или постоянной энтропии [5]. В этом случае нарушается положительная определенность тензора изотермических или адиабатических модулей упругости [9, 10].
Для исследования стабильности кристаллов разработан метод решеточной динамики [11]. Этот метод изначально создавался применительно к ионным кристаллам, и полученные выражения для модулей упругости ионных кристаллов не подходят для молекулярных кристаллов, поскольку в рамках этого метода модель кристалла рассматривается как точечный объект и не учитываются вращательные степени свободы молекулы.
Основной целью представленных в статье исследований является развитие метода решеточной динамики
для модельного описания на атомном уровне структурных, динамических, термодинамических и механических свойств с целью изучения стабильности льдов и клат-ратных гидратов в широком интервале давлений и установления природы фазового перехода «кристаллическая фаза—аморфная фаза» при изменении давления. На первом этапе ставится задача построения границ абсолютной стабильности льдов и гидратов.
Сравнение границ механической стабильности льдов и гидратов позволяет выявлять различие в свойствах обычных кристаллов и клатратов и исследовать роль молекул-гостей в устойчивости решетки хозяина. Построение границ стабильности необходимо для проведения экспериментальных исследований фазовых превращений при низких температурах и высоких давлениях. При этих условиях времена релаксации системы больше периода наблюдения и существует возможность приблизиться к границе абсолютной устойчивости кристалла.
Подходы к исследованию стабильности молекулярных кристаллов. Метод решеточной динамики
Модули упругости и критерий Борна. Механическая стабильность кристаллической решетки, означающая устойчивость ее по отношению к малым искажениям, определяется из обобщенного критерия стабильности Борна [11]. Согласно этому критерию, при данной температуре (энтропии) все детерминанты йк главных миноров матрицы Су (/, у = 1,...,к) изотермических (адиабатических) модулей упругости кристалла С у (/, У = 1,...,6) в обозначениях Фогта должны быть положительны.
Расчет упругих констант твердого тела для различных термодинамических условий может быть проведен разными способами. Адиабатические модули, возникающие под давлением, можно найти методом молекулярной динамики [12]. Этот метод был успешно использован для расчета модулей упругости твердого аргона [13] и льда !(! [4]. Серьезным недостатком этого подхода является большая статистическая ошибка и длительность расчета для получения результатов с высокой точностью [4, 13,]. Другим методом расчета модулей упругости молекулярных кристаллов может быть метод решеточной динамики в квазигармоническом приближении [9]. Свободная энергия кристалла Рць в этом подходе рассчитывается по формуле:
РФ = и + ГтЬ (1)
где и — потенциальная энергия; — колебательный вклад в свободную энергию. В отличие от ионных кристаллов для молекулярных кристаллов колебательный вклад в свободную энергию значителен и его необходимо учитывать. Если для ионных кристаллов он составляет порядка процента, то для льдов он уже равен порядка несколько десятков процентов.
Выражения для модулей упругости находятся как вторые производные от свободной энергии по компонентам тензора деформации кристалла [14]. Для молекулярных кристаллов в приближении жестких молекул выражение для модулей упругости получено в [9, 15].
Для определения упругих констант необходимо найти производные частот колебаний кристалла по деформациям. Эти производные можно численно рассчитать с помощью алгоритмов [16], которые дают значения час-
тот с гарантированной точностью для любой конфигурации кристаллической решетки.
Оптимизация геометрии молекулярного кристалла. В рамках метода решеточной динамики схема оптимизации геометрии хорошо разработана для ионного кристалла [17, 18]. Для молекулярного кристалла процедура оптимизации геометрии (внешних координат) существенно усложняется, так как необходимо учитывать вращательные степени свободы кристалла. В этом случае динамическая матрица кристалла включает трансляционные, ротационные и смешанные члены.
Особенность предлагаемого нами подхода к описанию молекулярных кристаллов заключается в том, что все структурные параметры оптимизируются одновременно с расчетом динамической матрицы и модулей упругости [19]. Этот подход позволяет находить структурные изменения и исследовать стабильность сложных соединений с наименьшими вычислительными усилиями. Для определения параметров элементарной ячейки при заданном давлении и температуре в нашем подходе используется процедура последовательного расчета параметров ячейки, при которых тензор напряжений имеет диагональный вид. При изменении давления Ар (однородного всестороннего сжатия) или температуры АТ находятся равновесные деформаций элементарной ячейки. Относительно найденных новых параметров элементарной ячейки определяется тензор напряжений. Если он окажется не диагональным, то для его диагона-лизации проводится последовательная корректировка равновесных деформаций элементарной ячейки, пока недиагональные элементы тензора напряжений станут меньше некоторой малой величины, определяющей точность расчета.
Для оптимизации внутренних координат положений центра масс и позиций молекул в элементарной ячейке используется метод Ньютона—Рапсона [20]. В случае кристаллов с большим числом молекул в элементарной ячейке применяется приближение внутренней статической деформации [21]. В рамках этого приближения при нахождении новых положений равновесия молекул при изменении температуры и объема вкладом в свободную энергию от колебаний кристаллической решетки пренебрегают. Задача сводится к минимизации потенциальной энергии молекулярного кристалла по отношению к внутренним координатам. Детали метода представлены в работе [19].
Критерий Линдемана для молекулярного кристалла. В качестве критерия плавления молекулярного кристалла мы выбрали критерий Линдемана [22], в котором вместо амплитуд колебаний атомов используются амплитуды колебаний центров масс молекул [23]. Согласно этому критерию в области стабильности кристалла вплоть до температуры плавления То должно выполняться неравенство
П(Т) = ((и2)/а2)1/2 < п(То) = По (2)
где а — ближайшее расстояние между центрами масс молекул; (и2) — средний квадрат амплитуд колебаний центров масс молекул.
В случае одноатомных кристаллов собственные вектора колебаний атомов в силу условий их ортогональности [24] не входят в выражение для среднего квадрата смещений атомов и здесь для нахождения (и2) необходимо рассчитать только собственные частоты
колебаний кристалла. Для молекулярного кристалла величина {u2) включает только собственные вектора колебаний центра масс молекул и не содержит собственных векторов, описывающих вращательные колебания молекул. В этом случае нельзя воспользоваться условием ортогональности полного набора собственных векторов и для расчета {u2) необходимо наряду с частотами рассчитывать и собственные вектора колебаний кристалла [23].
Термодинамическая и механическая стабильность льдов и VIII
В рамках предлагаемого метода самосогласованно находятся структура кристалла, его термодинамические функции и модули упругости. Метод был использован нами для численного построения кривых плавления и механической стабильности льдов ^ и VIII [25, 26].
Лед !|1. Для описания взаимодействия между молекулами воды применялся модифицированный Т1Р4Р потенциал [19], который получали из стандартного Т1Р4Р потенциала [28] путем изменения параметров взаимодействия с использованием перенормировочной константы K = 1,0066. При этом эффективные заряды умножали на K~2, параметр о, описывающий короткодействующее взаимодействие между атомами кислорода — на К-1, энергетический параметр е — на К3, а все расстояния между центрами взаимодействия на молекуле воды — на К-1. В случае использования модифицированного потенциала результаты расчета постоянной решетки льда I гораздо лучше согласуются с экспериментом (без ухудшения точности расчета других характеристик).
Расчеты были выполнены для суперячейки, содержащей 64 молекулы воды. Выбор такой большой ячейки позволяет моделировать протонное разупорядочение в структуре льда 1и. Расположение протонов на водородных связях принято в соответствии с правилом Берна-ла—Фовлера [27] так, чтобы дипольный момент элементарной ячейки был равен нулю. Влияние дально-действующего электростатического взаимодействия рассчитывали методом Эвальда. Свободную энергию и ее производные определяли численным методом. Вибрационная часть свободной энергии рассчитывалась по 2x2x2 точкам зоны Бриллюэна.
Ранее в работе [26] с целью сравнения точности метода решеточной динамики с методом молекулярной динамики и экспериментом были рассчитаны модули упругости льда ^ при 80 К и давлении до 3 кбар с использованием стандартного Т1Р4Р потенциала. Аналогичные вычисления были проведены с помощью модифицированного Т1Р4Р потенциала [29]. Как показывает рис. 1, относительное расположение расчетных модулей и наклоны кривых достаточно хорошо воспроизводят экспериментальные данные. Абсолютные значения рассчитанных и экспериментально найденных модулей упругости совпадают в пределах 20%.
Борновский критерий механической стабильности
гексагонального кристалла, каковым является структура льда 1и, включает условия: С44 > 0; Сц > |С12|; (Сц + С12) С33 > 2С132. Расчеты показали, что условие стабильности С11 > |С12| нарушается при повышении давления, так что это условие и определяет границу механической устойчивости льда 1и. Полученные результаты хорошо согласуются с экспериментальными данными по аморфизации льда при низких температурах и высоких давлениях. Среднеквадратичное смещение
160 -
&
MS
s" н о о
и
^
&
120 -
80 -
^ 40 "
Рис. 1. Зависимость адиабатических модулей упругости Оц льда !и от давления при 238 К.
Зачерненные значки — теоретические значения, светлые значки — эксперимент
центра масс молекулы воды в точке плавления Т = 273 К и р = 1 бар льда I согласно расчету по формуле (2) составляет П0 = 0,16. Построенная кривая термодинамического плавления, описывающая плавление льда I при низких давлениях и высоких температурах, хорошо согласуется с экспериментальными данными вплоть до давлений 4 кбар [25].
Лед VIII. Для кристаллической модификации льда VIII в широком интервале температур и давлений была построена кривая плавления, рассчитаны модули упругости как функции давления и определена граница механической стабильности. Элементарная ячейка тетрагонального протонно-упорядоченного кристалла льда VIII содержит 8 молекул. В расчетах использовались экспериментально полученные (при р = 24 кбар и Т = 15 К) параметры элементарной ячейки: а = 4,656 Е, с = 6,775 Е как начальные и затем проводилась оптимизация толь ко внутренних координат молекул и не оптимизировалась форма ячейки.
На рис. 2 представлены зависимости модулей упругости от давления при 160 К. Как видно, модули упругости С11, С13, С33, С44 и С66 уменьшаются при понижении давления, оставаясь положительными во всем интерва-
500
400
&
s
о 300 о
и
^
&
^ 200
100
Сп
С33
__С13
1С44
С12 -4
—«— С66
10
14 р, кбар
18
22
Рис. 2. Зависимость расчетных изотермических модулей упругости C/j льда VIII от давления при 160 K
0
0
6
400 300 ^ 200
100 -
600
500 -
«
400
300
10
15 20 р, кбар
25
30
35
10
20 30 40 50 р, кбар
60
70
Рис. 3. Линии плавления кристаллических модификаций льда VIII и VII.
Для льда VIII: экспериментальная граница механической устойчивости по критерию Борна ( А ); расчет по критерию Линдемана (о); расчет по уравнению Клапейрона—Клаузиуса (■). Для льда VII: экспериментальные линии (+) [31] и (х) [33]
ле давлений. Для широких интервалов температур (80— 300 К) и давлений расчеты изотермических модулей упругости для льда V!!! даны в [26].
Условия стабильности Борна для льда V!!! проанализированы для нескольких температур [26]. Критерием стабильности для тетрагонального кристалла являются условия [31]:
С11 > |С-|2|, (С11 + С12С33 > 2С 213, С44 > 0, Сбб > 0.
Расчеты стабильности льда V!!! выявляют нарушение условия С11 > |С-|2| при понижении давления. Следовательно, условие С11 = |С-|2| определяет границу механической стабильности этой модификации льда. Найденная граница стабильности льда V!!! по Борну приведена на рис. 3, здесь же дана линия плавления льда V!!!, рассчитанная по критерию Линдемана, а также гипотетическая линия плавления льда V!!!, построенная с помощью уравнения Клапейрона—Клаузиуса [12]. Значение П0 = 0,1468 в тройной точке равновесия жидкость—лед V!—лед V!!! при Т = 358 К и р = 24 кбар определено экстраполяцией линий фазовых равновесий жидкость—лед V! [31] и лед V!—лед V!!! [32]. На этих же рисунках приводятся экспериментально найденные линии плавления льда V!! [31, 33] (две ветви точек «+» на правом рисунке обусловлены разной точностью проведения эксперимента).
Рассчитанная линия плавления льда V!!! близка к границе механической стабильности и хорошо совпадает с экспериментальной линией плавления льда V!! при высоких давлениях.
Механическая стабильность клатратных гидратов
лостей. За исходную конфигурацию была принята элементарная ячейка гидрата структуры !, содержащая 46 молекул воды и 8 молекул гостей с параметром решетки а = 12Е при Т=245 К. Позиции атомов кислорода и координаты молекул-гостей были определены из результатов структурного анализа [34]. Расположение протонов выбрано таким, чтобы, как и в случае льда, полный дипольный момент элементарной ячейки равнялся нулю. Для описания взаимодействия между молекулами воды также использовали модифицированный Т!Р4Р потенциал. Молекулы-гости рассматриваются как сферически симметричные леннард-джонсоновские частицы.
Результаты расчета границы механической стабильности пустой гидратной решетки структуры ! даны на рис. 4. Показанные для сравнения на рис. 4 спинодали льда !(! при отрицательных и положительных давлениях близки к соответствующим границам пустой гидратной
Т, К -250 II
-150 ■50
_i_I_
0
0
5
Как и в структурах льда, в клатратных гидратах каждый атом кислорода связан водородными связями с четырьмя ближайшими атомами кислорода, образуя льдоподобные кристаллические структуры (решетка хозяина). В отличие от кристаллических структур льдов в этих структурах существуют полости, размеры которых позволяют включать только отдельные атомы или молекулы (гости).
В рамках предложенного нами подхода проведен расчет упругих констант и исследована механическая стабильность пустой гидратной кубической решетки и с метановым и ксеноновым заполнением гидратных по-
-16 -12 -8 -4 0 4 8 12 16 р, кбар
Рис. 4. Границы механической устойчивости льда I (■) и пустой гидратной решетки кубической структуры I (□) при положительных и отрицательных давлениях.
♦ _— рассчитанная линия фазового перехода первого рода между этими фазами, ! — область термодинамической стабильности пустой решетки гидрата, !! — область стабильности льда 11-,. Для сравнения приведены данные для спинодалей льда 1(, ( А ), рассчитанные методами молекулярной динамики [34]
решетки. Это говорит о том, что пустая гидратная решетка может существовать как механически устойчивая, но метастабильная фаза. Результаты расчетов, полученные методом молекулярной динамики [35], приведены для сравнения.
При повышении температуры характер изменения механической устойчивости кристалла в области положительных давлений имеет такую особенность — сначала устойчивость незначительно растет, а затем слабо уменьшается. Это может быть связано с температурными зависимостями параметров элементарной ячейки и модулей упругости. При нагревании кристалла решетка расширяется и для того, чтобы приблизиться к границе механической устойчивости кристалла, необходимо прикладывать большие давления. При повышении температуры меняется также соотношение между модулями упругости, что вызывает снижение устойчивости. Эти два эффекта и определяют поведение границы абсолютной устойчивости кристалла. В области отрицательных давлений расширение кристалла при повышении температуры приводит к понижению устойчивости кристалла.
Область существования решетки структуры I лежит внутри области стабильности льда Iм. Расчеты показывают, что пустая решетка хозяина может существовать при положительных давлениях, но как метастабильная фаза по отношению ко льду, и только при отрицательных давлениях становится стабильной. Линия фазового перехода «лед I м—пустая решетка гидрата структуры I» находится в области отрицательных давлений (см. рис. 4), что означает возможность фазового перехода в другую, более рыхлую кристаллическую структуру с полостями атомного размера при равномерном растяжении кристаллической решетки льда. Несмотря на то, что такая структура термодинамически неустойчива при положительных давлениях, она становится стабильной при заполнении полостей атомами подходящего размера. В работе [36] при использовании обобщенного уравнения Ван-дер-Ваальса построена линия структурного фазового перехода первого рода «лед ^—пустая решетка гидрата структуры II», что указывает на возможность существования при отрицательных давлениях и пустой решетки гидратов структуры II.
Результаты расчета модулей упругости Сц, С12, С44 для метанового и ксенонового гидратов при 80 К и давлении до 3,5 кбар представлены на рис. 5. С увеличени-
160 -
&
120
н о о
и
Л
ч
80
40 -
0,5 1,0
1,5 2,0 2,5 р, кбар
3,0 3,5
ем давления модуль Сц для обоих гидратов возрастает в равной степени (наклоны кривых равны с(Сц/ф = 6,3), модуль С12 возрастает несколько сильнее, причем в разной степени (с(С12/ф ° 7,45 для ксенонового гидрата и сС|2/ф ° 8,09 для метанового гидрата), модуль С44 слабо возрастает для ксенонового гидрата (dC44/dр ~ 0,245) и слабо снижается для метанового гидрата (dC 44/dр = 0,328).
Для обоих гидратов построены границы механической стабильности и спинодали (рис. 6) . Для ксеноно-вого гидрата граница механической стабильности и спинодали совпадают при отрицательных давлениях, в области положительных давлений граница стабильности начинается с давления 25 кбар, а спинодаль — при более высоких давлениях (63 кбар). Для кристаллов физический смысл имеет только граница механической стабильности. В большинстве случаев, как показали расчеты, спинодаль совпадает с границей механической стабильности, однако есть случаи, когда эти границы не совпадают. Поэтому для расчета границы абсолютной устойчивости кристалла необходимо использовать только условия Борна, а не условие dр/dV = 0. Как видно из рис. 6, граница стабильности для метанового гидрата при положительных давлениях лежит между границами стабильности пустой решетки и ксенонового гидрата. Это говорит о том, что введение в полость водного каркаса молекул-гостей приводит к росту механической устойчивости решетки. В то же время, поскольку ван-дер-ваальсов радиус молекулы метана заметно меньше, чем аналогичная величина для ксенона, можно предположить, что при увеличении размера молекулы-гостя область устойчивости гидрата расширяется. Вычисленные давления начала механической неустойчивости метанового и ксенонового гидратов соответствуют значениям, наблюдаемым в отношении других гидратов [3].
Исследование проведено для кристаллов с водородными связями, но можно предположить, что аналогичные фазовые переходы реализуются и для других кристаллов. Проведение экспериментов при отрицательных давлениях (равномерного растяжения, например, кристалла с достаточно большой концентрацией атомов
Рис. 5. Зависимость адиабатических модулей упругости Оц для ксенонового (зачерненные символы) и метанового (светлые символы) гидратов от давления при 80 К
0 5 10 15 20 25 30 р, кбар
Рис. 6. Границы механической устойчивости пустой решетки гидрата кубической структуры I ( ), ксенонового (■) и метанового (♦) гидратов при положительных и отрицательных давлениях
0
0
другого сорта) позволит выявить существование новых кристаллов клатратного типа, область устойчивости
которых значительно больше, чем исходного кристалла.
* * *
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 96-032-928) и Президиума Сибирского отделения Российской академии наук (грант № 76 «Газовые гидраты в природных экосистемах»).
ЛИТЕРАТУРА
1. Mishima O., Calvert L.D., Whalley E. Nature, 1984, v. 310, p. 393.
2. Mishima O. Ibid., 1996, v. 384, p. 546.
3. Handa Y.P., Tse I.S., Klug D.D., Whalley E. J. Chem. Phys., 1991, v. 94, p. 623.
4. Whalley E., Klug D.D., Handa Y.P. Nature, 1989, v. 342, p. 782.
5. Tse J.S. J. Chem. Phys., 1992, v. 96, p. 5482.
6. Binggeli N., Chelikowsky J.R. Phys. Rev. Lett., 1992, v. 69, p. 2220.
7. Binggeli N., Keskar N.R., Chelikowsky J.R. Phys.Rev.B, 1994, v.
49, p. 3075.
8. Keskar N.R., Chelikowsky J.R., Wentzcovitch R.M. Ibid., 1994, v.
50, p. 9072.
9. Shpakov V.P., Tse J.S., Belosludov V.R., Belosludov R.V. J. Phys.: Cond. Mat., 1997, v. 9, p. 5853.
10. Zhou Z., Joos B. Phys. Rev., 1996, v. B54, p. 3841.
11. Born M., Huang H. Dynamical Theory of Crystal Lattices. Oxford: Clarendon Press, 1954.
12. Parrinello M., Rahman A. J. Chem. Phys., 1982, v. 76, p. 2662.
13. Sprik M., Imprey R.W., Klein M.L. Phys. Rev. B, 1983, v. 29, p. 4368.
14. Leibfried G., Ludwig W. Theory of Anharmonic Effects in Crystals. New York—London: Academic Press Inc., 1961.
15. Белослудов В.Р., Шпаков В.П., Це Д.С. и др. Химия в интересах устойчивого развития, 1998, т. 6, с. 75.
16. Годунов С.К. Решение системы линейных уравнений, Новосибирск: Наука, 1980.
17. Kantorovich L.N. Phys. Rev. B, 1995, v. 51, p. 3520.
18. Taylor M.B., Barrera G.D., Allan N.L., Barron T.H.K. Ibid., 1997, v. 56, p. 14380.
19. Belosludov V.R., Shpakov V.P., Tse J.S., Belosludov R.V., Kawazoe Y. Ann. N.Y. Acad. Sci., 2000, v. 912, p. 993.
20. Gill P.E., Mirray W,, Wright M.H. Practical Optimization. London : Academic Press Inc., 1981.
21. Watson G.M., Tschaufeser P., Wall A., Jackson R.A., Parker S.C. In: Computer Modeling in Inorganic Crystallography. Ed. C.R.A. Catlow. San Diego: Academic, 1997, p. 55.
22. Lindemann. F.A. Z. Phys., 1910, Bd. 11, S. 609.
23. Belosludov R.V., Kawazoe Y., Grachev E.V., Dyadin Yu.A., Belosludov V.R. Solid State Comm., 1999, v. 109, p. 157.
24. Беттер Х. Принципы динамической теории решетки. М.: Мир, 1986, c. 396.
25. Tse J.S., Klug D.D., Tulk C.A. e. a. Nature, 1999, v. 400, p. 647.
26. Tse J.S., Shpakov VP,, Belosludov V.R. J. Chem. Phys., 1999, v. 111, p. 11111.
27. Bernal J.D., FowlerR.H. Ibid., 1933, v. 1, p. 515.
28. Jorgensen W.L., Chandrasekhar J., Madura J.D. e. a. Ibid., 1983, v. 79, p. 926.
29. Gagnon R.E., Kiefte H., Clouter M.J., Whalley E. Ibid., 1988, v. 89, p. 4522.
30. Fedorov F.I. Theory of Elustic Waves in Crystal. New York: Plenum Press, 1968.
31. Tkachev S.N., Nasimov R.M., Kalinin V.A. J. Chem. Phys., 1996, v. 105, p. 3722.
32. Brown A.J., Whalley E. Ibid., 1966, v. 45, p. 4360.
33. Fein Y., Mao H., Hemley R.J. Ibid., 1993, v. 99, p. 5369.
34. McMullan R.K., Jeffrey G.A. Ibid., 1965, v. 42, p. 2725.
35. Sciortino F., Essmann U., Stanley H.E. e. a. Phys. Rev. B, 1995, v. 52, p. 6485.
36.Косяков В.И., Шестаков В.А. Докл. АН, 2000, т. 376, № 6, с. 1—3.
