Методы Молекулярной статики

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

Математическая задача минимизация функции многих переменных может быть сформулирована чисто формальным образом:

Задана функция U, которая однозначно зависит от некоторого числа N независимых переменных x₁, x₂,..., x_N. Конкретно для атомных систем N равно утроенному числу атомов в системе.

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

, (14.6)

а все собственные значения матрицы вторых производных (или матрицы Гессе),

, (14.7)

были положительны.

Строго говоря, в приложениях задача минимизации часто рассматривается в смягченном варианте и требуется найти значения переменных, обеспечивающих локальный минимум потенциальной энергии, наиболее близкий к некоторой точке фазового пространства (x₁*, x₂*,..., x_N*), которая используется в качестве первоначального предположения о положении минимума. Хотя, используя различные начальные предположения, можно, в принципе, просканировать всю поверхность потенциальной энергии, найти все минимумы и выбрать среди них самый глубокий, на практике такая тактика используется достаточно редко и решение ограничивается рассмотрением одного или нескольких минимумов для исходных конфигураций, выбираемых из физических соображений (например, из анализа симметрии кристалла).

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

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

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

• Третья группа методов получила название методов Ньютона. Это наиболее быстро сходящиеся численные алгоритмы, но для их использования необходимо знание самой функции, ее первых и вторых производных. Платой за скорость является объем памяти, необходимый для хранения матрицы Гессе. Он пропорционален N², что может быть непосильным для моделирования больших кристаллов. Поэтому была разработана группа алгоритмов – производных метода Ньютона, которые называются квази-ньютоновские методы. Основной идеей этих алгоритмов является использование не фактической матрицы Гессе, но ее приближенных значений. Типичными представителями этой группы являются алгоритмы Давидона-Флетчера-Пауэлла (DFP) или, например, Бройдена-Флетчера-Гольдфарба-Шанно (BFGS).

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

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

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

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

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

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

Минимизация методом наискорейшего спуска. Метод основан на аналогии с механикой, где движение системы материальных точек происходит в каждый момент времени в направлении действия приложенных к каждой точке силы. В методе наискорейшего спуска непрерывное движение разбивается на последовательность дискретных шагов, причем на каждом шаге атомы смещаются вдоль одного и того же направления в фазовом пространстве, совпадающего с направлением обобщенной «силы», действующих на атомы в начальной точке шага (направление наискорейшего спуска). В 3N мерном пространстве координат атомов, x = (x ₁, …, x _{3 N}), обобщенная сила определяется как взятый с обратным знаком градиент потенциальной энергии,

(14.8)

Итерационный шаг, для некоторой конфигурации системы x ⁿ (где индекс n обозначает номер шага процедуры минимизации), начинается с расчета градиента энергии в этой точке, g ⁿ. После чего атомы смещаются вдоль направления в 3N -мерном фазовом пространстве, где bⁿ - варьируемый параметр, который последовательно увеличивается до тех пор, пока не будет достигнута точка, где потенциальная энергия (как функция единственного свободного параметра bⁿ) достигает минимума. Значение , соответствующее локальному минимуму энергии, может быть найдено по любой стандартной методике определения минимума одномерной функции и обычно не требует повторной оценки градиента энергии, ограничиваясь использованием алгоритмов поиска. Полученная новая конфигурации системы, , используется в качестве начальной позиции для следующего шага итерации.

Метод наискорейшего спуска всегда приводит к минимуму энергии. Он прост в алгоритмической реализации и требует хранения только данных о градиенте. Тем не менее, эта схема минимизации далеко не всегда оптимальна. Для понимания причин этого воспользуемся простой житейской аналогией. Если спускаться по пересеченной местности, то направление наибольшего спуска приводит обычно на дно ближайшего оврага, вдоль которого и происходит дальнейший спуск. Аналогичная ситуация наблюдается и при минимизации функций: уже первые несколько шагов (или даже просто первый шаг) приводят в долину на потенциальной поверхности, а дальнейшие шаги совершаются более-менее в общем направлении дна долины. При этом, на каждом шаге при движении по прямой дно долины пересекается. Иными словами, спуск вдоль дна долины происходит по зигзагообразной траектории, осциллирующей вблизи дна долины. При достижении локального минимума энергии происходит поворот строго перпендикулярно предыдущему направлению. Это следует из того, что в методе наискорейшего спуска направление движения на шаге n +1, h ^{n +1}, совпадает с направлением градиента энергии g ^{n +1} в точке , а в этой точке проекция g ^{n +1} на направление минимизации равна нулю, так что

(14.9)

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

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

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

(1) h ⁿ будет выражено, согласно алгоритму Флетчера-Ривса, в виде простого рекуррентного соотношения,

, (4.10)

где gⁿ - скалярный фактор, и

(2) будет предположено, что вектора градиента в начале и в конце каждого шага - ортогональны,

. (4.11)

Заметим, что, хотя рассмотренный метод минимизации называется методом " сопряженных градиентов", минимизация идет вдоль сопряженных последовательных направлений смещения системы в фазовом пространстве, а не градиентов.

Условия ортогональности и сопряженности позволяют выразить скалярные константы и gⁿ как

(4.12)

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

В дальнейшем для улучшения сходимости схемы сопряженных градиентов, Полак и Рибьер предложили выбирать скалярные константы gⁿ в виде:

. (4.13)

Вблизи от точки минимума разницы между двумя определениями gⁿ практически нет, но вдали от нее метод Полака-Рибьера работает, как правило, лучше, чем оригинальный алгоритм Флетчера-Ривса. Когда потенциальная поверхность в непосредственной близости от точки минимума описывается чисто квадратичной функцией, эти константы совпадают с константами вида (4.12), так как все градиенты будут ортогональны. Однако на самом деле функция энергии только примерно квадратична и поэтому метод Полака-Рибьера работает лучше, чем оригинальный алгоритм Флетчера-Ривса. Метод сопряженных градиентов является одним из самых популярных методов минимизации функционала энергии Кона-Шэма в квантово-механических расчетах, так как обеспечивают сравнительно быструю сходимость и использование только первых производных энергии системы.

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

Идею метода можно проиллюстрировать, если рассмотреть одномерную функцию потенциальной энергии U(х), которая известна в некоторой точке xⁿ вместе со своими первыми и вторыми производными. Для малого отклонения расстояния dx от xⁿ потенциальная энергия и ее градиент могут быть представлены в виде ряда Тейлора:

, (14.14)

, (14.15)

где H (x) = d ² U / dx ². Когда точка xⁿ близка к точке x ⁰, в которой g обращается в нуль, нелинейными членами разложения можно пренебречь, и расстояние dx=x ⁰- xⁿ может быть напрямую определено из условия g (xⁿ + dx) = 0, что приводит к выражению

. (14.16).

Если, однако, точность нахождения нуля функции с помощью линейной аппроксимации градиента энергии недостаточна, мы можем пересчитать первую и вторую производные энергии в точке x^{n +1}= xⁿ + dx и сделать еще один шаг вдоль касательной к функции g. Учитывая, что n -номер шага в итерационной процедуре, получаем следующую рекуррентную формулу для нахождения минимума энергии:

. (14.17)

Те же рассуждения могут быть легко обобщены для многомерного случая, для которого выражение (4.17). принимает вид

, (14.18).

Где и - многомерный градиент и обратная матрица Гессе, соответственно.

Длина шага в методе Ньютона-Рафсона уменьшается с каждой итерацией

. (14.19).

Алгоритм Ньютона-Рафсона сходится квадратично, увеличивая количество значащих цифр в x ⁰ примерно в два раза после каждого шага. Эта особенность делает этот метод хорошо применимым для нахождения минимума функции потенциальной энергии в случае, если ее первые и вторые производные могут быть легко определены.

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

Существует набор вариаций метода Ньютона-Рафсона, целью разработки которых являлась попытка обойти необходимость расчета полной матрицы вторых производных. Например, семейство методов, называемых квази-ньютоновскими методами, требует расчета только первых производных и пошагового построения обратной матрицы Гессе в процессе процедуры минимизации. Приставка «квази» добавляется потому, что в расчетах используется не реальная матрица Гессе, а ее приближение для данного текущего шага.

Один из простых способов, которые позволяют ускорить процедуру минимизации методом Ньютона-Рафсона, заключается в использовании одной и той же матрицы Гессе для нескольких последовательных шагов алгоритма, пересчитывая на каждой итерации только градиенты. Одним из наиболее широко используемых алгоритмов является блочно-диагональный метод Ньютона-Рафсона, в котором на каждой итерации смещается только один атом. Это позволяет оставить только те члены, которые включают в себя координаты перемещаемого атома, и задача сводится к тривиальному обращению матрицы 3x3. Тем не менее, блочно-диагональной подход может быть не эффективным, если движения некоторых атомов системы строго коррелированы.

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

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