Потенциал

1.2. Потенциал

Для моделирования материала необходимо задать потенциал взаимодействия составляющих его атомов. Наиболее простым является парный потенциал Леннарда-Джонса

Здесь,  - расстояние между атомами,  - глубина потенциальной ямы и  связано с положением минимума потенциала . Потенциал Леннарда-Джонса качественно правильно описывает взаимодействие между атомами – сильное отталкивание на малых расстояниях, обусловленное первым слагаемым в скобках, и притяжение на больших расстояниях, за которое отвечает второе слагаемое в скобках. Он хорошо описывает ван-дер-ваальсовское взаимодействие между атомами кристаллов благородных газов, но, вследствие своей простоты, часто используется для качественного описания взаимодействия других атомов. С потенциалом Леннарда-Джонса проведено большое количество вычислений. Он является стандартным в вычислениях методом МД.

Основными материалами реакторостроения являются металлы – сталь, цирконий и т.д. В металлах природа сил взаимодействия между атомами не двухчастичная (парная) а многочастичная. Effective Medium Theory (EMT) дает реалистическое описание металлической связи с учетом её многочастичной природы [3,4]. EMT- потенциал, с вычислительной точки зрения, не намного сложнее парного потенциала, но дает намного более реалистическое описание свойств материалов. Поскольку в данной работе не ставится задача изучения пластических свойств конкретного материала мы будем использовать модельный потенциал Леннарда-Джонса.

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

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

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

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

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

Конечно, эти схемы приближенными, и, поэтому, существуют ошибки, связанные с ними. Они классифицируются так:

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

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

Оба типа ошибок можно уменьшить, уменьшая . Для больших  ошибки обрывания доминируют, но они быстро уменьшаются, когда  уменьшается. Например, алгоритм Верле имеет ошибки обрывания пропорциональные  для каждого временного шага интегрирования. Ошибки округления падают более медленно с уменьшением  и доминируют в пределе малых . Использование 64-битной точности (соответствующую “двойной точности” в Fortrane) помогает сохранить ошибки округления минимальными.

В молекулярной динамике наиболее часто используемым алгоритмом интегрирования по времени является, вероятно, так называемый алгоритм Верле [5]. Основная идея состоит в том, чтобы записать разложение Тейлора до третьего порядка вперед и назад по времени. Пусть  обозначает скорость,  - ускорение и  - третью производную от  по . Тогда имеем:

Складывая эти 2 выражения получаем

Это основная формула алгоритма Верле. Так как мы интегрируем уравнения Ньютона, то  есть просто сила, деленная на массу, и сила в свою очередь есть функция положения :

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

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

Однако, ошибки, которые дает это выражение, порядка  а не .

Чтобы преодолеть эту трудность, были развиты варианты алгоритма Верле. Они дают точно ту же траекторию и отличаются переменными, которые хранятся в памяти. Leap-frog алгоритм есть один из таких вариантов.

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

Заметим, что необходимо 9N ячеек памяти, чтобы сохранить 3N положений, скоростей и ускорений, но нам не нужно одновременно хранить значения любой из этих величин для двух различных времен.

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

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

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

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

Таким образом, алгоритм поиска атомов, удаленных от данного атома на расстояние не больше радиуса обрезания , выглядит следующим образом. По номеру атома находим координаты атома и по ним субячейку, в которой находится атом. Затем находим субячейки, удаленные от нее на расстояние не более чем . Атомы, расположенные в этих субячейках, и будут искомыми (см. рис.1). Чтобы найти номер атома, хранящегося в заданной субячейке, удобно ввести массив, каждый элемент которого соответствует определенной субячейке. В этом элементе массива будет хранится номер атома, расположенного в этой субячейке, или нуль, если субячейка пуста. Элементы этого массива обновляются на каждом шаге по времени МД. Ясно, что изложенный алгоритм обеспечивает линейный рост числа операций с ростом числа атомов  в системе. Вариации этого алгоритма используются в программах МД “Gromex”[6], “MOLDY”[7], “DL_POLY”[8] и др.

Возможна и другая организация вычислений, которая будет удобна для организации параллельных вычислений. Именно для вычисления сил, действующих на данный атом, можно перейти от суммирования по близлежащим атомам, к суммированию по близлежащим субячейкам (см рис.1). Будем двигаться последовательно по субячейкам первого ряда. Дойдя до конца первого ряда, перейдем в начало второго ряда и т.д.

Рис.1 Схема поиска ближайших атомов.

Если в субячейке находится атом, то вычисляем силу, действующую на него, со стороны ближайших атомов, расположенных в близлежащих субячейках. Если же субячейка пуста, то переходим к следующей. Отметим при этом, что, например, для атома находящегося в субячейке 6 (см. рис.1) необходимо вычислить силу, действующую со стороны атомов расположенных в субячейках 1, 2, 3, 7. Силы, действующие со стороны атомов, расположенных в субячейках 5, 9, 10, 11 в силу третьего закона Ньютона, с точностью до знака уже известны. Они были вычислены, когда вычислялись силы, действующие на атомы, расположенные в этих субячейках. Таким образом, в данной организации вычислений, необходимо рассматривать лишь половину близлежайших субячеек. Далее, при переходе к смежной субячейке 7 нет необходимости исследовать все близлежащие субячейки для поиска находящихся в них близко расположенных атомов. Необходимо лишь исследовать ячейки 4 и 8. И к найденным в них атомам, добавить атомы, найденные для ячейки 6, за исключением атомов находящихся в субячейках 1 и 6. Таким образом, информация о ближайших атомах для данной субячейки не теряется, а используется при поиске ближайших атомов для смежной субячейки. Это естественно приводит к ускорению вычислений.

Похожие работы

Моделирование как метод физической мезомеханики

Скачать

40807

0

0

... базе Института физики прочности и материаловедения СО РАН издается на русском и английском языках международный журнал «Физическая мезомеханика». 4. Моделирование как средство экспериментального исследования Моделирование всегда используется в комплексе с другими общенаучными и специальными методами. Теснее всего моделирование связано с экспериментом. Попробуем разобраться, в чем отличие ...

Основы теории и технологии контактной точечной сварки

Скачать

318063

13

95

... , при обработке металлов давлением. Экспериментальные исследования процессов пластической деформации металла в зоне формирования соединения при контактной точечной сварке по этой методике проводятся на натурных образцах с предварительно нанесенной координатной сеткой, технология изготовления которых предложена и описана в работе [128]. При исследованиях пластических деформаций в плоскостях ...

Снижение энергозатрат толстолистовых станов

Скачать

59925

3

5

... повысить пластичность поверхностных слоев на 40%. Это позволяет на 10-15% снизить энергозатраты на прокатку.   Применение технологических смазок   Применение традиционных смазочных материалов (минерального, растительного масла, синтетических смазок с содержанием поверхностно-активных веществ и др.) на толстолистовых станах трудно реализуемо из-за дефицитности и недостаточной эффективности, ...

Основы проектирования и конструирования

Скачать

460103

24

39

... ребрами) изображают конструктивные и потоковые функциональные структуры [14]. Принципы построения функциональных структур технических объектов рассматриваются в последующих главах курса "Основы проектирования им конструирования" не включенных в настоящее пособие. Для систем управления существуют характеристики, которые можно использовать в качестве критериев для оценки структур. Одна из них - ...