Расчет термодинамических величин в методе молекулярной динамики

4.1. Общие соотношения

Когда после определенного числа шагов МД в системе установилось равновесие, можно “измерять” ее термодинамические параметры. Для этого значения соответствующих динамических величин усредняются во времени. Физические свойства динамических величин, как правило, являются функциями координат и скоростей. Допустим, мгновенное значение некоторой физической величины в момент времени t равно

. (4.1)

Тогда ее среднее значение равно

, (4.2)

где k – индекс, который пробегает шаги по времени от 1 до некоторого большого числа . Имеются два эквивалентных пути практической реализации такого усреднения.

1. рассчитывается на каждом шаге по времени в процессе “прогонки” программы МД. Сумма также обновляется после каждого шага. В конце работы программы МД среднее значение находится сразу делением суммы на число шагов. Этот способ предпочтителен, когда расчет величины прост или очень важен. К примеру, температура системы должна рассчитываться таким образом.

2. В процессе работы программы МД координаты (и/или скорости) периодически записываются в файл. Отдельная программа, запускаемая после МД, обрабатывает записанные данные, рассчитывая необходимые средние. Этот способ требует большого объема на диске, так как связан с сохранением величин (координат и проекций скоростей всех частиц) для каждого из большого числа шагов. Однако он удобен в другом отношении: программа МД прогоняется только один раз, а ее результаты могут быть много раз использованы для расчета различных термодинамических величин. Этот метод используется, если рассчитываемые величины сложны, связаны с использованием данных в разные моменты времени и т.д.

4.2. Потенциальная, кинетическая и полная энергия

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

. (4.3)

Хотя для интегрирования уравнений движения средняя потенциальная энергия не нужна, она нужна для расчета полной энергии и проверки ее сохранения. Последняя является важной проверкой любого МД-моделирования.

Расчет средней кинетической энергии очень прост:

. (4.4)

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

В действительности могут иметь место небольшие флуктуации полной энергии порядка 10^‑4 от среднего значения или меньше, связанные с погрешностями в интегрировании уравнений движения. Эти флуктуации могут быть уменьшены сокращением шага по времени.

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

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

4.3. Температура

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

. (4.5)

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

4.4. Калорическая кривая. Определение температуры плавления твердых тел

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

Одним из интересных и важных задач, решаемых с помощью МД, является исследование фазовых переходов, например, плавления кристалла, и определение температуры плавления . В принципе можно определять, увеличивая температуру кристалла до тех пор, пока не появится скачок в калорической кривой, указывая на поглощение латентной энергии. Однако температура, при которой произойдет этот скачок, неизменно выше температуры плавления. По определению, температура плавления – это температура, при которой твердая и жидкая фазы сосуществуют, имея одинаковую свободную энергию. Из-за отсутствия зародыша, на котором при нагревании росла бы жидкая фаза, всегда имеет место перегрев выше . В области перегрева система находится в термодинамически метастабильном состоянии, которая стабильна в течение времени моделирования. Когда достигается точка механической неустойчивости перегретого кристалла, он распадается. Эта точка может соответствовать исчезновению одного из модулей сдвига материала или подобным неустойчивостям, и, как правило, выше примерно на 20-30%.

Более точный метод определения температуры плавления состоит в следующем. Первоначально создается достаточно большая система, состоящая примерно на 50% из твердого тела и на 50% ‑ из жидкости. Такую систему можно создать, закрепив атомы одной части кристалла в фиксированных положениях, доведя до расплавления другую часть и сняв затем все ограничения. В дальнейшем с помощью МД достигается равновесие между жидкостью и твердым телом, и в равновесной системе измеряется температура указанным ранее способом.

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

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

4.5. Среднеквадратичные отклонения и диффузия

По определению, среднеквадратичное отклонение частиц в течение времени t равно

, (4.6)

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

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

. (4.7)

Коэффициент пропорциональности D называется коэффициентом диффузии. Таким образом, с помощью МД можно определить коэффициент диффузии:

. (4.8)

4.6. Расчет давления в молекулярной динамике

Давление в МД рассчитывается с помощью теоремы вириала. Для вывода этой теоремы рассмотрим величину

. (4.9)

где сумма берется по всем частицам системы, а ‑ полная сила, действующая на i- ю частицу. Посчитаем среднее по времени от этой величины:

, (4.10)

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

, (4.11)

где ‑ единичный вектор нормали к площадке . Таким образом, для второй части величины G имеем

, (4.12)

где была использована теорема Гаусса. В результате можно записать следующее равенство:

, (4.13)

или

, (4.14)

где теперь произведения в сумме берутся только с силами межчастичного взаимодействия.

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

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

. (4.15)

4.7. Расчет атомных напряжений в молекулярной динамике

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

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

. (4.16)

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

, (4.17)

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

, (4.18)

где по повторяющимся индексам подразумевается суммирование. Энергия деформированного твердого тела в единице объема определяется соотношением

. (4.19)

Отсюда следует, что тензор напряжений может быть определен как тензор с компонентами

. (4.20)

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

. (4.21)

Средние по системе напряжения тогда определяются как

. (4.22)

По полной аналогии с этим, локальное напряжение около атома i (атомное напряжение) определяется следующим соотношением:

, (4.23)

где вместо объема системы V теперь используется атомный объем .

В приближении парного взаимодействия энергия атома равна

. (4.24)

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

, (4.25)

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

. (4.26)

Отсюда легко определить атомное напряжение у атома i:

. (4.27)

При выводе последней части этой формулы было использовано соотношение

. (4.28)

При использовании МПА выражения для тензора напряжений несколько более сложные, но их расчет также не представляет больших трудностей.

Из выражения (4.27) видно, что для расчета локального напряжения в месте расположения данного атома необходимо знать объем, который он занимает. Когда атом расположен в бездефектной части материала, этот объем достаточно хорошо определен, и его можно определить, например, построением полиэдров Вороного. Однако вблизи дефектов, таких как трещины, поверхность, локальный объем плохо определен, и в этих местах невозможно определить напряжения. В отличие от напряжений, величины всегда однозначно определены, и, как правило, программы МД выводят значения именно этих величин.

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