Реферат: Математическое моделирование пластической деформации кристаллов

--PAGE_BREAK--Конечно, эти схемы приближенными, и, поэтому, существуют ошибки, связанные с ними. Они классифицируются так:
Ошибки обрывания, связанные с точностью метода конечных разностей по отношению к истинному решению. Метод конечных разностей обычно базируется на ряде Тейлора, оборванном на некотором члене. Эти ошибки не зависят от программной реализации метода, они присущи самому алгоритму.
Ошибки округления, связаны с ошибками, возникающими при программной реализации алгоритма. Например, такие ошибки возникают из-за конечного числа цифр, используемых в компьютерной арифметике.
Оба типа ошибок можно уменьшить, уменьшая <shape id="_x0000_i1058" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image061.wmz» o:><img width=«20» height=«19» src=«dopb219315.zip» v:shapes="_x0000_i1058">. Для больших <shape id="_x0000_i1059" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image063.wmz» o:><img width=«20» height=«19» src=«dopb219315.zip» v:shapes="_x0000_i1059"> ошибки обрывания доминируют, но они быстро уменьшаются, когда <shape id="_x0000_i1060" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image063.wmz» o:><img width=«20» height=«19» src=«dopb219315.zip» v:shapes="_x0000_i1060"> уменьшается. Например, алгоритм Верле имеет ошибки обрывания пропорциональные <shape id="_x0000_i1061" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image064.wmz» o:><img width=«27» height=«21» src=«dopb219316.zip» v:shapes="_x0000_i1061"> для каждого временного шага интегрирования. Ошибки округления падают более медленно с уменьшением <shape id="_x0000_i1062" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image063.wmz» o:><img width=«20» height=«19» src=«dopb219315.zip» v:shapes="_x0000_i1062"> и доминируют в пределе малых <shape id="_x0000_i1063" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image063.wmz» o:><img width=«20» height=«19» src=«dopb219315.zip» v:shapes="_x0000_i1063">. Использование 64-битной точности (соответствующую “двойной точности” в Fortrane) помогает сохранить ошибки округления минимальными.
В молекулярной динамике наиболее часто используемым алгоритмом интегрирования по времени является, вероятно, так называемый алгоритм Верле [5]. Основная идея состоит в том, чтобы записать разложение Тейлора до третьего порядка вперед и назад по времени. Пусть <shape id="_x0000_i1064" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image066.wmz» o:><img width=«13» height=«19» src=«dopb219317.zip» v:shapes="_x0000_i1064"> обозначает скорость, <shape id="_x0000_i1065" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image068.wmz» o:><img width=«13» height=«19» src=«dopb219318.zip» v:shapes="_x0000_i1065"> - ускорение и <shape id="_x0000_i1066" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image070.wmz» o:><img width=«15» height=«23» src=«dopb219319.zip» v:shapes="_x0000_i1066"> - третью производную от <shape id="_x0000_i1067" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image072.wmz» o:><img width=«13» height=«17» src=«dopb219320.zip» v:shapes="_x0000_i1067"> по <shape id="_x0000_i1068" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image057.wmz» o:><img width=«11» height=«16» src=«dopb219313.zip» v:shapes="_x0000_i1068">. Тогда имеем:
Складывая эти 2 выражения получаем
Это основная формула алгоритма Верле. Так как мы интегрируем уравнения Ньютона, то <shape id="_x0000_i1072" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image080.wmz» o:><img width=«28» height=«23» src=«dopb219324.zip» v:shapes="_x0000_i1072"> есть просто сила, деленная на массу, и сила в свою очередь есть функция положения <shape id="_x0000_i1073" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image082.wmz» o:><img width=«28» height=«23» src=«dopb219325.zip» v:shapes="_x0000_i1073">:
Видно, что ошибка обрывания алгоритма, когда система эволюционирует в течении времени <shape id="_x0000_i1075" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image086.wmz» o:><img width=«20» height=«19» src=«dopb219315.zip» v:shapes="_x0000_i1075">, есть величина порядка <shape id="_x0000_i1076" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image087.wmz» o:><img width=«27» height=«21» src=«dopb219316.zip» v:shapes="_x0000_i1076">, даже если третья производная не появляется в явном виде. Этот алгоритм в тоже время является простым для программной реализации, точным и стабильным, что объясняет его большую популярность при МД моделировании.
Проблема с этой версией алгоритма Верле состоит в том, что скорости явно не вычисляются. Хотя они не нужны для временной эволюции, но их знание иногда необходимо. Кроме того, они нужны для вычисления кинетической энергии <shape id="_x0000_i1077" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image088.wmz» o:><img width=«17» height=«17» src=«dopb219327.zip» v:shapes="_x0000_i1077">, чья оценка необходима чтобы проверить сохранение полной энергии <shape id="_x0000_i1078" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image090.wmz» o:><img width=«73» height=«19» src=«dopb219328.zip» v:shapes="_x0000_i1078">. Это один из наиболее важных тестов, указывающих, что МД моделирование выполняется корректно. Можно вычислить скорости из положений использую формулу
Однако, ошибки, которые дает это выражение, порядка <shape id="_x0000_i1080" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image094.wmz» o:><img width=«27» height=«21» src=«dopb219330.zip» v:shapes="_x0000_i1080"> а не <shape id="_x0000_i1081" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image096.wmz» o:><img width=«27» height=«21» src=«dopb219316.zip» v:shapes="_x0000_i1081">.
Чтобы преодолеть эту трудность, были развиты варианты алгоритма Верле. Они дают точно ту же траекторию и отличаются переменными, которые хранятся в памяти. Leap-frog алгоритм есть один из таких вариантов.
Лучший реализацией того же основного алгоритма есть так называемый алгоритм Верле со скоростью, когда положение скорости и ускорения в момент времени <shape id="_x0000_i1082" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image097.wmz» o:><img width=«41» height=«19» src=«dopb219314.zip» v:shapes="_x0000_i1082"> получается из тех же величин в момент времени <shape id="_x0000_i1083" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image098.wmz» o:><img width=«11» height=«16» src=«dopb219313.zip» v:shapes="_x0000_i1083">следующим образом:
Заметим, что необходимо 9N ячеек памяти, чтобы сохранить 3N положений, скоростей и ускорений, но нам не нужно одновременно хранить значения любой из этих величин для двух различных времен.

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

1.5. Вычисление сил Наибольших вычислительных усилий требует вычисление сил, действующих между атомами. Поэтому оптимизации алгоритма вычисления сил необходимо уделить особое внимание. Один из шагов в этом направлении состоит в замене сложных для вычисления выражений для сил (например, содержащих экспоненту) на легко вычисляемые выражения (например, сплайны третьего порядка). Второй шаг состоит в использовании потенциалов с ограниченным радиусом действия, или, как указывалось выше, в обрезании несущественной области потенциала, если радиус действия потенциала бесконечен. При этом необходимо вычислить только силы, действующие со стороны ближайших атомов, т.е. находящихся внутри сферы (окружности в двумерном случае) с радиусом равным радиусу обрезания <shape id="_x0000_i1088" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image107.wmz» o:><img width=«20» height=«24» src=«dopb219309.zip» v:shapes="_x0000_i1088">.
Третий шаг состоит в оптимизации алгоритма поиска атомов, ближайших к данному атому. Дело в том, что прямолинейный перебор всех атомов, вычисление расстояний до них и отбрасывание тех атомов, расстояние до которых превышает радиус обрезания <shape id="_x0000_i1089" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image108.wmz» o:><img width=«20» height=«24» src=«dopb219309.zip» v:shapes="_x0000_i1089">, требует количества операций пропорционального <shape id="_x0000_i1090" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image109.wmz» o:><img width=«24» height=«21» src=«dopb219335.zip» v:shapes="_x0000_i1090">, где <shape id="_x0000_i1091" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image111.wmz» o:><img width=«19» height=«19» src=«dopb219336.zip» v:shapes="_x0000_i1091"> - число атомов в системе. Следовательно, с  ростом <shape id="_x0000_i1092" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image113.wmz» o:><img width=«19» height=«19» src=«dopb219336.zip» v:shapes="_x0000_i1092"> число требуемых операций быстро возрастает, и поэтому выполнение вычислений сильно замедляется, а, для больших <shape id="_x0000_i1093" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image114.wmz» o:><img width=«19» height=«19» src=«dopb219336.zip» v:shapes="_x0000_i1093">, делается практически невыполнимым. Таким образом, чтобы избежать этого замедления нужен алгоритм, для которого число требуемых операций росло бы с <shape id="_x0000_i1094" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image115.wmz» o:><img width=«19» height=«19» src=«dopb219336.zip» v:shapes="_x0000_i1094"> линейно, а не квадратично. В принципе такой алгоритм прост – нужно перебирать не все атомы, а только достаточно близкорасположенные. Такое утверждение представляет собой тавтологию, пока не конкретизировано понятие близкорасположенных атомов. Чтобы сделать это, разобьем ячейку моделирования на более мелкие субячейки. Тогда близкорасположенные к данному атому будут атомы, которые расположены в субячейках, соседних с субячейкой, содержащей данный атом или в субячейках соседних с соседними.
Удобно разбить ячейку моделирования на субячейки – параллелепипеды (прямоугольники в двумерном случае). Вследствие сильного отталкивания на малых расстояниях, атомы не могут подходить близко друг к другу. Поэтому можно выбрать такие размеры субячеек, что в каждой из них будет находится не более одного атома.
Таким образом, алгоритм поиска атомов, удаленных от данного атома на расстояние не больше радиуса обрезания <shape id="_x0000_i1095" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image116.wmz» o:><img width=«20» height=«24» src=«dopb219309.zip» v:shapes="_x0000_i1095">, выглядит следующим образом. По номеру атома находим координаты атома и по ним субячейку, в которой находится атом. Затем находим субячейки, удаленные от нее на расстояние не более чем <shape id="_x0000_i1096" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image117.wmz» o:><img width=«20» height=«24» src=«dopb219309.zip» v:shapes="_x0000_i1096">. Атомы, расположенные в этих субячейках, и будут искомыми (см. рис.1). Чтобы найти номер атома, хранящегося в заданной субячейке, удобно ввести массив, каждый элемент которого соответствует определенной субячейке. В этом элементе массива будет хранится номер атома, расположенного в этой субячейке, или нуль, если субячейка пуста. Элементы этого массива обновляются на каждом шаге по времени МД. Ясно, что изложенный алгоритм обеспечивает линейный рост числа операций с ростом числа атомов <shape id="_x0000_i1097" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image118.wmz» o:><img width=«19» height=«19» src=«dopb219336.zip» v:shapes="_x0000_i1097"> в системе. Вариации этого алгоритма используются в программах МД “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. Таким образом, информация о ближайших атомах для данной субячейки не теряется, а используется при поиске ближайших атомов для смежной субячейки. Это естественно приводит к ускорению вычислений.


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

1.7. Начальное состояние В данной работе будут исследоваться с помощью МД кристаллы. Рассмотрим размещение совершенного кристалла в прямоугольной ячейке моделирования в случае периодичности в одном направлении. Периодическая структура самого кристалла накладывает ограничения на размер ячейки моделирования в направлении периодичности. Действительно, если в вершине, находящейся на одной из сторон ячейки моделирования находится атом кристалла, то эквивалентный атом кристалла должен быть в эквивалентной вершине, находящейся на другой из тождественных сторон. Это приводит к ограничениям на возможную длину ячейки моделирования в направлении периодичности и возможные ориентации кристаллографических осей кристалла относительно сторон ячейки моделирования. Возможны такие ориентации кристалла, при которых указанное выше требование выполнить точно невозможно.
<imagedata src=«45520.files/image119.wmz» o:><img width=«449» height=«357» src=«dopb219337.zip» v:shapes="_x0000_i1098">
Рис.2 Периодичность ячеек моделирования и размещение кристалла в ячейке моделирования.
Если же ориентация кристалла выбрана удачно, то длина ячейки моделирования может принимать значения кратные некоторой величине. Однако, на самом деле, эти ограничения не очень существенны. Для всех длин ячейки моделирования и ориентаций кристалла можно найти близкие к ним значения этих величин, для которых условие будет выполняться точно. Рецепт состоит в том, чтобы просто совместить указанную эквивалентную вершину с ближайшим эквивалентным атомом кристалла.
Если есть периодичность (см. рис. 2) и по второму направлению, то должно выполняться аналогичное требование и для второго направления. При этом необходимо заметить, что ориентация второй стороны  для прямоугольной ячейки моделирования уже задана, поскольку она перпендикулярна первой стороне. Поэтому её длина будет кратна некоторой величине.
Если не принять специальных мер при подготовке начального состояния системы, то в ней возникают коллективные движения — колебания. Это связано с тем, что система может оказаться в сжатом или растянутом состоянии из-за несоответствия температуры системы с постоянной кристаллической решетки. Другими словами это тепловое расширение (сжатие) системы. Такие колебания имеют большой период и слабо затухают. Накладываясь на исследуемый процесс (например, деформирование системы) они смазывают картину этого исследуемого процесса. Следовательно, от этих колебаний необходимо избавиться. Это можно сделать несколькими способами. Во-первых, подождать пока колебания затухнут. Однако из-за большого периода и малого затухания это требует большого времени. Во-вторых, попытаться подогнать постоянную решетки кристалла к температуре. Опыт показывает, что, сделав несколько попыток, можно полностью исключить колебания. В-третьих, такую подгонку можно выполнить автоматически. О том, как это можно сделать, будет сказано ниже.
В МД моделировании часто возникает необходимость иметь систему в состоянии, характеризуемом определенной температурой. Однако, как мы можем получить систему с заданной температурой? Другими словами, как мы можем контролировать систему?
Для изменения температуры необходимо так изменить скорости частиц, чтобы получить желаемую температуру. В алгоритме Верле со скоростью, обсуждаемом выше, это может быть выполнено заменой уравнения
на уравнение
где <shape id="_x0000_i1101" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image124.wmz» o:><img width=«17» height=«24» src=«dopb219339.zip» v:shapes="_x0000_i1101"> желаемая температура, и <shape id="_x0000_i1102" type="#_x0000_t75" o:ole="" fillcolor=«window»><imagedata src=«45520.files/image126.wmz» o:><img width=«29» height=«23» src=«dopb219340.zip» v:shapes="_x0000_i1102"> текущая температура. Такая модификация означает, что мы больше не следуем уравнениям Ньютона и, что полная энергия больше не сохраняется.

1.8. Начальное состояние для кристаллов с дефектами С помощью МД можно исследовать деформирование, как совершенных кристаллов, так и кристаллов содержащих дефекты, например, кристаллов подвергнутых облучению. О том, как подготовить начальное состояние для совершенного кристалла, было сказано выше. Подготовка начального состояния для облученного кристалла намного более сложная задача. Однако, если известны доза и спектр первично выбитых атомов, МД позволяет выполнить моделирование каскада повреждений [9,10,11]и таким образом решить эту сложную задачу. При этом описанные выше потенциалы, необходимо дополнить, чтобы учесть отталкивание на малых расстояниях, например, гладко сшивая их с потенциалом Циглера-Бирсака-Литмарка [12]. Такой подход позволяет учесть многие явления, возникающие при облучении, но является достаточно сложным и лежит за рамками данной работы.
    продолжение
--PAGE_BREAK--