Основы метода молекулярной динамики

3.1. Физическая основа метода молекулярной динамики

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

В классической МД решается система уравнений II закона Ньютона:

, (3.1)

где m_i ‑ масса i -й частицы, ‑ ее ускорение, и ‑ сила, действующая на нее. МД является детерминистическим методом: если известны начальные координаты и скорости частиц, в принципе может быть описана эволюция системы в любые моменты времени. В то же самое время, МД представляет собой метод статистической физики. Компьютер рассчитывает траекторию системы в 6 N -мерном фазовом пространстве (3 N координат и 3 N импульсов). Однако эта траектория сама по себе интереса не представляет. Главный смысл использования метода заключается в том, что генерируется совокупность конфигураций, распределенных по определенным законам распределения, то есть статистический ансамбль. Примерами таких ансамблей является микроканонический ансамбль, соответствующий распределению вероятностей в фазовом пространстве, обеспечивающем сохранение энергии, или канонический ансамбль, в котором постоянна температура. Траектория, полученная МД, дает такую совокупность конфигураций. Таким образом, измерение какой-либо физической величины путем моделирования сводится к усреднению различных мгновенных значений, принимаемых этой величиной в различные последовательные моменты времени.

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

3.2. Краткая история развития метода МД

Впервые метод МД был использован Алдером и Уэйнрайтом (Al-der, Weinwright) в 1957 г. В их работе была исследована фазовая диа-грамма системы твердых сфер, то есть переходы между жидким и твердым состояниями этих систем. Расчеты были проведены на компьютере IBM-704 (имеющем объем оперативной памяти 4 Кб).

В 1960 г. Гибсоном, Голандом, Милграмом и Виньярдом в Брукхавенской национальной лаборатории было проведено МД моделирование радиационного дефектообразования в меди. Это было первое использование МД с непрерывным потенциалом и конечно-разностным интегрированием уравнений движения. Расчеты для системы 500 атомов были проведены на компьютере IBM-704 и занимали около 1 мин. на один шаг по времени.

Рамэн (Rahman, Аргонская нац. лаб.) – один из известных пионеров МД. В 1964 г. он изучил с помощью потенциала ЛД свойства жидкого Ar. Рассматривалась система из 864 атомов на компьютере CDC-3600. Л. Верле (Verlet) в 1967 г. также изучал свойства жидкого Ar, в этих исследованиях он ввел понятие списка соседей Верле, а также метод интегрирования уравнений движения, называемый алгоритмом Верле, которые в настоящее время широко используются.

В 1981 г. Паринелло (Parinello) и Рамэн предложили метод МД при постоянном давлении, в 1984 г. Нозе (Nosé) предложил метод МД при постоянной температуре. Теперь эти методы являются классическими.

В увеличением возможностей компьютеров росли также размеры моделируемых систем и спектр задач, решаемых с ее помощью. В 1993 г. Ломдалем (Lomdahl) было реализовано моделирование роста трещины в системе, содержащей 250 млн. атомов, а уже в 2000 г. была смоделирована система из 5 млрд. атомов.

3.3. Обзор основных задач, решаемых с помощью МД

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

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

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

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

Изучение фазовых превращений, в том числе между агрегатными состояниями одного и того же вещества, построение фазовых диаграмм также составляют одно из направлений применения МД.

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

3.4. Ограничения классической МД

В МД квантовое уравнение Шредингера заменяется классическими уравнениями Ньютона. Эта замена справедлива, когда длина волны Де-Бройля l, соответствующая тепловому движению атома, много меньше кратчайшего межатомного расстояния 3´10^‑10 м. Имея для длины волны Де-Бройля и скорости теплового движения частиц формулы

, , (3.2)

где M ‑ масса частиц, T ‑ температура, из условия можно получить соотношение, определяющее применимость классического описания:

, (3.3)

где ‑ масса протона (равная примерно 1 а.е.м.). Подставляя значения физических постоянных в соотношение (3.3), можно показать, что классические уравнения справедливы, если . Это означает, что для всех атомов, за исключением наиболее легких, применение классических уравнений движения оправдано. Классическое МД неприменимо для таких легких элементов, как H₂, He, Ne.

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

Кроме этих, физически обоснованных условий, существуют также ограничения, связанные с возможностями интегрирования уравнений движения. Эти ограничения сводятся к ограничению размеров изучаемых систем и промежутков времени, в течение которых могут рассматриваться процессы. В типичных задачах на сегодняшний день рассматриваются системы, содержащие миллионы атомов, и процессы в течение времени порядка 1 нс. Для систем меньшего размера промежутки времени могут быть увеличены до сотен наносекунд, и, соответственно, для меньших промежутков времени размеры системы могут быть увеличены до нескольких миллиардов. Для того чтобы представить масштаб такой системы, отметим, что для никеля, атомный объем которого равен примерно 10^-2 нм³, система из 1 млрд. атомов занимает объем куба со сторонами примерно 300 нм.

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

3.5. Инициализация систем для моделирования в МД. Периодические граничные условия

Прежде чем начать моделирование эволюции состояния системы атомов во времени с помощью МД, необходимо эту систему инициализировать, то есть описать. Для инициализации необходимо выполнить следующие пункты: 1) описать силы (или потенциал) межатомного взаимодействия; 2) задать исходное состояние, то есть координаты и скорости частиц; 3) задать граничные условия. Первый пункт был рассмотрен выше,- результатом его является выбор, как правило, из имеющихся в наличии потенциалов, потенциала для данного вещества. Как это делается, мы рассмотрим при изучении конкретных пакетов, реализующих МД. Задание исходного состояния – важная часть моделирования, именно она определяется физическим подходом и постановкой задачи моделирования. Соответственно, эта часть определяется спецификой решаемой задачи. Как решается эта задача для конкретных приложений к нанотехнологии и наноматериалам, будет рассмотрено в соответствующей части в конце курса. Важным общим пунктом, который необходимо реализовать при моделировании любого объекта, является выбор граничных условий, поэтому рассмотрим это понятие подробно.

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

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

Однако чаще всего необходимо изучать части систем намного больших размеров, например, дефекты в макроскопических системах. В этом случае моделируемая система представляет собой только часть большой системы, содержащую интересующий нас дефект, а атомы на границе моделируемой системы в физической системе принадлежат объему кристалла. Если эти атомы считать свободными, то влияние поверхности на состояние системы будет значительным. Действительно, отношение числа атомов поверхности (пропорционального R ²) к числу атомов системы (пропорциональному R ³) уменьшается с уменьшением размеров системы как 1/ R. Для макросистем это отношение пренебрежимо мало, а при ограниченном числе атомов оно становится достаточно большим. Влияние поверхности для малых систем становится намного больше, чем было бы для макросистем.

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

Наиболее часто используемыми являются периодические граничные условия (ПГУ), в английской терминологии periodic boundary conditions (PBCs). Суть этих условий заключается в следующем.

Частицы заключаются в ячейке в форме параллелепипеда, называемой расчетной ячейкой. Эта ячейка повторяется бесконечное число раз путем трансляции во всех трех направлениях, заполняя все пространство (см. рис. 5). Иными словами, если одна частица расположена с координатами в расчетной ячейке, считают, что эта частица в действительности представляет бесконечную совокупность частиц, расположенных в точках , где ‑ целые числа и ‑ векторы, соответствующие ребрам расчетной ячейки. Все эти частицы-“изображения” движутся вместе, и только одна из них представляется в компьютерной программе. Ключевым является предположение, что каждая частица i в расчетной ячейке взаимодействует не только с остальными частицами в этой же ячейке, а также с их образами в соседних ячейках. То есть, взаимодействия могут осуществляться через границы расчетной ячейки. Это приводит к тому, что фактически (а) мы исключаем влияние поверхности на систему, и (б) положение границ ячейки относительно частиц не играет роли (то есть, трансляция ячейки по отношению к частицам не меняет сил). В процессе движения какая-либо частица системы может покинуть расчетную ячейку. В таком случае через противоположную границу ячейки в нее входит периодический образ этой частицы, так что число атомов в ячейке не изменяется.

Рис. 5. Схема, поясняющая периодические граничные условия

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

3.6. Правило ближайшей частицы

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

Чтобы это продемонстрировать, предположим, что i взаимодействует с двумя образами j1 и j2 частицы j. Два образа должны быть разделены вектором трансляции, переводящим одну ячейку в другую, длина которой, согласно нашему предположению, не меньше . Для того чтобы взаимодействовать с обеими частицами j1 и j2, частица i должна иметь расстояние от каждой из них, не превышающее , что невозможно.

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

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

3.7. Методы интегрирования уравнений движения

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

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

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

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

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

Имеются два популярных метода интегрирования, используемых в МД: алгоритм Верле (Verlet) и алгоритм предиктор-корректор.

3.8. Алгоритм Верле

Этот алгоритм, возможно, самый популярный в МД. Он вытекает из разложения радиуса-вектора частиц в два момента времени, и , в ряд Тейлора до третьей степени по . Учитывая, что первая производная по времени есть скорость , а вторая – ускорение , эти разложения можно записать в виде

, (3.4)

. (3.5)

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

. (3.6)

Это и есть алгоритм Верле. Поскольку мы интегрируем уравнения Ньютона, ускорения легко считаются через силы:

. (3.7)

Как видно из уравнения (3.7), ошибка усечения алгоритма Верле имеет четвертый порядок по шагу по времени. Кроме того, алгоритм очень прост и устойчив, поэтому он пользуется большой популярностью. Недостатком алгоритма является меньшая точность определения скоростей частиц. Хотя для интегрирования уравнений движения расчет скорости не обязателен, знание скоростей важно. Например, для контроля устойчивости решения уравнений движения необходимо проверять сохранение полной энергии , а для расчета кинетической энергии надо знать скорости. Выражение для скорости может быть получено путем вычитания уравнения (3.5) из (3.4):

. (3.8)

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

Для преодоления этого недостатка используются различные варианты алгоритма Верле. Эти алгоритмы обеспечивают эволюцию системы по той же траектории, но отличаются моментами времени, в которые в памяти сохраняются скорости и координаты. Наиболее простым из этих вариантов является алгоритм лягушачьих прыжков (leap-frog algorithm). В этом алгоритме вычисляются скорости в полуцелые шаги времени, которые используются для расчета нового положения частицы. Этот алгоритм получается из алгоритма Верле, если по формуле (3.8) мы посчитаем скорости в моменты времени и :

, . (3.9)

Из второго уравнения находим следующее выражение для расчета новых координат на основе старых координат и скоростей:

. (3.10)

Из алгоритма Верле находим следующее выражение для скорости:

. (3.11)

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

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

, (3.12)

. (3.13)

. (3.14)

. (3.15)

В этом алгоритме нужно сохранять 9 N величин – по 3 N координат, скоростей и ускорений, но только в один момент времени.

3.9. Списки соседей атомов

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

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

При большом количестве частиц в системе (N >1000 или более, в зависимости от радиуса обрезания потенциала), становится более предпочтительным другой метод. Кубическая расчетная ячейка (возможно распространение метода и на некубические ячейки) делится на субъячеек таких, что сторона каждой из них больше, чем . Списки атомов составляются для каждой из этих субъячеек, и при расчете взаимодействий любого атома перебираются только атомы, лежащие в той же и соседних субъячейках. Распределение атомов по субъячейкам требует расчетов, объем которых пропорционален N, поэтому может проделываться на каждом шаге без особого ущерба быстроте расчетов.

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