Математическое моделирование пластической деформации кристаллов (стр. 4 из 7)
Возможна и другая организация вычислений, которая будет удобна для организации параллельных вычислений. Именно для вычисления сил, действующих на данный атом, можно перейти от суммирования по близлежащим атомам, к суммированию по близлежащим субячейкам (см рис.1). Будем двигаться последовательно по субячейкам первого ряда. Дойдя до конца первого ряда, перейдем в начало второго ряда и т.д.
| 1 | 2 | 3 | 4 |
| 5 | 6 | 7 | 8 |
| 9 | 10 | 11 | 12 |
Рис.1 Схема поиска ближайших атомов.
Если в субячейке находится атом, то вычисляем силу, действующую на него, со стороны ближайших атомов, расположенных в близлежащих субячейках. Если же субячейка пуста, то переходим к следующей. Отметим при этом, что, например, для атома находящегося в субячейке 6 (см. рис.1) необходимо вычислить силу, действующую со стороны атомов расположенных в субячейках 1, 2, 3, 7. Силы, действующие со стороны атомов, расположенных в субячейках 5, 9, 10, 11 в силу третьего закона Ньютона, с точностью до знака уже известны. Они были вычислены, когда вычислялись силы, действующие на атомы, расположенные в этих субячейках. Таким образом, в данной организации вычислений, необходимо рассматривать лишь половину близлежайших субячеек. Далее, при переходе к смежной субячейке 7 нет необходимости исследовать все близлежащие субячейки для поиска находящихся в них близко расположенных атомов. Необходимо лишь исследовать ячейки 4 и 8. И к найденным в них атомам, добавить атомы, найденные для ячейки 6, за исключением атомов находящихся в субячейках 1 и 6. Таким образом, информация о ближайших атомах для данной субячейки не теряется, а используется при поиске ближайших атомов для смежной субячейки. Это естественно приводит к ускорению вычислений.
1.6. Периодичность
Число атомов, помещенных в ячейку моделирования, намного меньше числа атомов входящих в состав макроскопических систем. Чтобы результаты нашего моделирования можно было распространить на макроскопические тела, делают допущение, что макроскопические системы, состоят из бесконечного числа периодически повторяющихся ячеек моделирования. Такая периодичность может быть в одном, двух и трех направлениях в трехмерном случае и в одном и двух в двумерном случае (см. рис.2). В этой работе мы будем рассматривать только двумерные системы. Это связано как с повышенными требованиями к вычислительным ресурсам в случае трехмерных систем, так и с простотой визуализации результатов расчетов в двумерном случае. В двумерном случае ячейка моделирования представляет собой прямоугольник. В случае периодичности в одном направлении пара противолежащих сторон отождествляется, т.е. ячейку моделирования можно представить теперь как боковую поверхность цилиндра. В случае периодичности в двух направлениях отождествляются обе пары противоположных сторон и ячейку моделирования можно теперь представить как боковую поверхность тора. Если атом выходит за пределы ячейки моделирования, то вследствие периодичности он входит в ячейку с противоположной стороны.
1.7. Начальное состояние
В данной работе будут исследоваться с помощью МД кристаллы. Рассмотрим размещение совершенного кристалла в прямоугольной ячейке моделирования в случае периодичности в одном направлении. Периодическая структура самого кристалла накладывает ограничения на размер ячейки моделирования в направлении периодичности. Действительно, если в вершине, находящейся на одной из сторон ячейки моделирования находится атом кристалла, то эквивалентный атом кристалла должен быть в эквивалентной вершине, находящейся на другой из тождественных сторон. Это приводит к ограничениям на возможную длину ячейки моделирования в направлении периодичности и возможные ориентации кристаллографических осей кристалла относительно сторон ячейки моделирования. Возможны такие ориентации кристалла, при которых указанное выше требование выполнить точно невозможно.
Рис.2 Периодичность ячеек моделирования и размещение кристалла в ячейке моделирования.
Если же ориентация кристалла выбрана удачно, то длина ячейки моделирования может принимать значения кратные некоторой величине. Однако, на самом деле, эти ограничения не очень существенны. Для всех длин ячейки моделирования и ориентаций кристалла можно найти близкие к ним значения этих величин, для которых условие будет выполняться точно. Рецепт состоит в том, чтобы просто совместить указанную эквивалентную вершину с ближайшим эквивалентным атомом кристалла.
Если есть периодичность (см. рис. 2) и по второму направлению, то должно выполняться аналогичное требование и для второго направления. При этом необходимо заметить, что ориентация второй стороны для прямоугольной ячейки моделирования уже задана, поскольку она перпендикулярна первой стороне. Поэтому её длина будет кратна некоторой величине.
Если не принять специальных мер при подготовке начального состояния системы, то в ней возникают коллективные движения - колебания. Это связано с тем, что система может оказаться в сжатом или растянутом состоянии из-за несоответствия температуры системы с постоянной кристаллической решетки. Другими словами это тепловое расширение (сжатие) системы. Такие колебания имеют большой период и слабо затухают. Накладываясь на исследуемый процесс (например, деформирование системы) они смазывают картину этого исследуемого процесса. Следовательно, от этих колебаний необходимо избавиться. Это можно сделать несколькими способами. Во-первых, подождать пока колебания затухнут. Однако из-за большого периода и малого затухания это требует большого времени. Во-вторых, попытаться подогнать постоянную решетки кристалла к температуре. Опыт показывает, что, сделав несколько попыток, можно полностью исключить колебания. В-третьих, такую подгонку можно выполнить автоматически. О том, как это можно сделать, будет сказано ниже.
В МД моделировании часто возникает необходимость иметь систему в состоянии, характеризуемом определенной температурой. Однако, как мы можем получить систему с заданной температурой? Другими словами, как мы можем контролировать систему?
Для изменения температуры необходимо так изменить скорости частиц, чтобы получить желаемую температуру. В алгоритме Верле со скоростью, обсуждаемом выше, это может быть выполнено заменой уравнения
 | (9) |
на уравнение
 , | (10) |
где
желаемая температура, и 
текущая температура. Такая модификация означает, что мы больше не следуем уравнениям Ньютона и, что полная энергия больше не сохраняется.1.8. Начальное состояние для кристаллов с дефектами
С помощью МД можно исследовать деформирование, как совершенных кристаллов, так и кристаллов содержащих дефекты, например, кристаллов подвергнутых облучению. О том, как подготовить начальное состояние для совершенного кристалла, было сказано выше. Подготовка начального состояния для облученного кристалла намного более сложная задача. Однако, если известны доза и спектр первично выбитых атомов, МД позволяет выполнить моделирование каскада повреждений [9,10,11]и таким образом решить эту сложную задачу. При этом описанные выше потенциалы, необходимо дополнить, чтобы учесть отталкивание на малых расстояниях, например, гладко сшивая их с потенциалом Циглера-Бирсака-Литмарка [12]. Такой подход позволяет учесть многие явления, возникающие при облучении, но является достаточно сложным и лежит за рамками данной работы.
Можно также исследовать влияние определенных дефектов, возникающие при облучении ГПУ кристаллов на их пластические свойства. Например, можно исследовать влияние межузельных кластеров и дефектов Френкеля. Очевидно, что начальные состояния, содержащие такие дефекты, легко приготовить, стартуя с начального состояния для идеального кристалла. Для этого необходимо удалить (добавить, переместить) атомы кристалла так, чтобы получилась конфигурация кристалла с требуемыми дефектами. Кристалл при этом получается обычно в напряженном состоянии. Это справедливо особенно при добавлении атомов, так как для добавленных атомов расстояния до ближайших атомов кристалла обычно намного меньше, чем равновесные расстояния между атомами в кристалле. Из-за сильного роста потенциала межатомного взаимодействия на малых расстояниях такие атомы обладают большой потенциальной энергией. Если не принять специальных мер, это может вызвать разлет кристалла. Чтобы не допустить этого и обеспечить релаксацию напряжений можно использовать процедуру минимизации и последующий подогрев системы до нужной температуры.
1.9. Нагрузка
В данной роботе рассматривалось деформирование кристаллов путем одноосного растяжения. Поскольку вдоль направления растяжения наложены периодические граничные условия, то отсутствуют свободные границы, к которым можно было бы приложить нагрузку. Поэтому задается растяжение системы, и потом находится возникшее вследствие этого напряжение. МД и деформирование выполняются одновременно. После каждого шага по времени МД выполняется малое растяжение, обеспечивающее нужную скорость деформации (
на одном шаге). Растяжение выполнялось двумя способами. В первом, традиционно используемом [13], система растягивается равномерно по длине. При этом координаты атомов вдоль направления растяжения умножаются на масштабный множитель 
. На этот же множитель умножается длина ячейки моделирования. Согласно второму способу, предложенному в данной работе, растяжение концентрируется только возле торцов системы. Этот способ лучше соответствует экспериментальной ситуации, когда нагрузка прикладывается к торцам системы. При этом длина ячейки моделирования умножается на масштабный множитель, а координаты атомов не умножаются.