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

Схема алгоритма

Первая часть алгоритма моделирования по методу молекулярной динамики такая же, как в методе Монте-Карло. Строится исходная конфигурация атомов в кубе, размер которого определяется плотностью моделируемого материала. Координаты атомов у исходной конфигурации получаются с помощью генератора случайных чисел. Задаётся парный потенциал, описывающий взаимодействие атомов. Обычно используются такие же потенциалы, как в методе Монте-Карло (2.3, 2.4). Основная часть алгоритма – релаксация в методах Монте-Карло и молекулярной динамики существенно различны.

В методе молекулярной динамики релаксация модели производится путём решения уравнений классической механики Ньютона для каждого из атомов. При этом атомы рассматриваются как взаимодействующие точки. Каждый атом движется в поле сложной формы, создаваемом всеми остальными атомами [4,6,19-21].

Столь сложная задача может быть решена только приближённо. Шкала времени разбивается на равные интервалы , которые настолько малы, что на протяжении одного интервала движение каждого атома может рассматриваться как равноускоренное. Изменения координат атома с номером “” на протяжении интервала с номером “” описываются следующими выражениями:

(2.9)

Здесь – векторы положения, скорости и ускорения атома номером “”. Скорость на “” интервале принимает значение согласно:

(2.10)

Ускорение вычисляется путём деления на массу “”­го атома суммарной силы, действующей со стороны всех остальных атомов

(2.11)

Вышеприведённый алгоритм был предложен Верле, по имени которого он и получил название. При его реализации используются такие же граничные условия, как и в описанном выше методе Монте-Карло.

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

По значению средней кинетической энергии определяется температура модели.

(2.12)

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

Одно из принципиальных отличий метода молекулярной динамики от метода Монте-Карло состоит в том, что при выполнении релаксации по методу молекулярной динамики имитируется реальное движение атомов в реальном времени. В методе Монте-Карло при использовании алгоритма Метрополиса, строго говоря, время отсутствует. Хотя в ряде случаев имеет смысл рассматривать последовательность пробных сдвигов как последовательность временных интервалов. Такая интерпретация «напрашивается», например, при рассмотрении результатов по релаксации стекла (стр. В этом случае используется термин «марковское время», так как последовательность структур, получаемая при пробных сдвигах, образует марковскую цепь.

Рис 2.8 Индивидуальные ФРР для кварцевого стекла. Сравнение данных для модели, построенной методом молекулярной динамики (сплошная линия), с результатами анализа рассеяния рентгеновских лучей (пунктирная линия). число атомов кислорода, расположенных внутри сферы радиуса вокруг атома кремния [19]. На этом рисунке Si ­ O, Si ­ Si и O ­ O то же, что , и согласно соотношению (2.1).

Реальное время, которое затрачивается на релаксацию по методу молекулярной динамики, определяется длительностью элементарного временного интервала. Обычно, считают, что она должна быть меньше периода тепловых колебаний. В одной из первых работ по моделированию методом молекулярной динамики структуры стёкол для кварцевого стекла полагалось, что =10^-14 с.[19]. Общее число таких элементарных шагов, затраченное на весь процесс релаксации, составляло 10⁵. При этом температура модели уменьшилась на насколько тысяч градусов. Таким образом, скорость охлаждения модели превышала 10¹¹­10¹² град./сек. Эта величина на много порядков превышает достигнутые в настоящее время самые большие скорости охлаждения при закалке. Следовательно, структура модели должна сопоставляться с структурой очень сильно закалённых стёкол.

Структурные характеристики моделей

Сравнение ФРР, рассчитанных для модели, с результатами анализа дифракционных экспериментов показало, что модель вполне удовлетворительно описывает эксперимент. Особенно выразительно выглядит плоский участок зависимости от обобщённого координационного числа , соответствующий значению 4. Это служит убедительным доказательством автоматического образования тетраэдрической координации кремния в чисто ионной модели.

Большое достоинство ионных моделей, при использовании методов Монте-Карло и молекулярной динамики, состоит в том, что не возникает каких – либо проблем при переходе от рассмотрения стёкол простейшего состава, содержащих один стеклообразующий окисел, к стёклам сколь угодно сложного состава, содержащих несколько окислов.

Расчёт колебательного спектра

Так как в моделировании методом молекулярной динамики имитируются реальные движения атомов, он позволяет рассчитывать динамические характеристики материала.

Примером одной из самых важных таких характеристик является колебательный спектр . Для его получения вычисляется автокорреляционная функция скоростей согласно следующему выражению [20]:

(2.13)

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

Колебательный спектр определялся путём вычисления интеграла

(2.14)

Результаты такого расчёта для кварцевого стекла представлены на рис. 2.9. Величина временного шага в этом расчёте составляла 10^-15с. Отметим, что под колебательным спектром здесь понимается функция плотности колебательных состояний. При сравнении ИК спектр поглощения света при возбуждении этих состояний

Рис 2.9 ФРР и колебательный спектр для аморфного кремния [20]

Коэффициент самодиффузии.

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

(2.15)

Здесь – временной интервал, на котором определяется диффузия; число атомов определённого типа (в случае кварца – кремния или кислорода. положения атома с номером в начале и в конце численного эксперимента. Усреднение проводится по множеству различных начальных моментов отсчёта времени.

Эхо гашения

Компьютерное моделирование позволяет обнаруживать новые нетривиальные явления, которые невозможно наблюдать в реальных экспериментах. Примером такого явления, является «эхо гашения», которое описывается в настоящем разделе [22].

Явление было неожиданно обнаружено при моделировании методом молекулярной динамики колебательных движений в расплавах и стёклах очень простого модельного состава. Модель состояла из 500 частиц, взаимодействующих согласно потенциалу Ленарда-Джонса

(2.16)

Здесь – глубина потенциала, – радиус частицы, а – расстояние между частицами. Алгоритм интегрирования использовал шаг по времени 10^-14 с. Мгновенная температура определялась по значению средней кинетической энергии частицы

(2.17)

В последнем выражении среднее берётся по ансамблю частиц.

Рис. 2.10 Зависимость от времени температуры модели из 500 частиц, дважды подвергавшейся гашению в моменты и . Сигнал эха наблюдается в момент [22].

Эхо гашения возникает при вполне определённом режиме охлаждения модели.

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

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

(в) Затем в момент времени скорости всех частиц зануляются вторично. Казалось бы, что после второго гашения должно всё происходить точно так же, как после первого. Сначала всё так и происходит.

(г) Однако, в момент второе гашение вызывает совершенно неожиданный результат: температура понижается весьма заметно сама по себе без каких-либо внешних воздействий. Такое спонтанное понижение температуры получило название «эха гашения».

Этот любопытный эффект получил объяснение на основе очень простой модели (рис 2.11). Модель рассматривается как ансамбль осцилляторов, различающихся по частоте колебаний. Каждый осциллятор движется в своей потенциальной яме, имеющей форму параболы. Кривизна парабол и перед колебаний у осцилляторов варьируют в широких пределах.

2.11 Схема, иллюстрирующая происхождение эха гашения

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

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

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

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

В частности, в результате гашения «выживают» моды, для которых период задаётся режимом охлаждения. Если процедуру повторять многократно с одним и тем же интервалом , то сохранится колебание только одной частоты (или её гармоники). Меняя интервал , можно получит весь спектр нормальных колебаний.