Диссертация (1102956), страница 6
Текст из файла (страница 6)
При коэффициенте kr существенно большем коэффициента kφ в выражении ( 4)(потенциал F3C) в координатах r − φ функция U intra будет «оврагом». Скорость сходимостиметода градиентов можно улучшить, используя динамический метод [66]. Динамическийметод реализуется следующим образом: выбираются 2 параметра - шаг по времени и коэффициент уменьшения скорости на каждом шаге β = 0.995. В данном случае нам не нужнонаходить истинную траекторию движения к минимуму, важна лишь скорость сходимости.Итерационный процесс представим в следующем виде:~ = v~n − 1vn+1iimi29∑ ~∇Uinj (dt)j<i(11)~ = β v~nvn+1ii(12)~ (dt)~ = x~n + vn+1xin+1ii(13)В методе молекулярной динамики рассчитываются классические траектории движенияатомов в силовом поле эмпирического атом - атомного потенциала, то есть моделируется детальная микроскопическая картина теплового движения.
Основу метода составляетчисленное решение уравнений движения:1~x¨ = ~fimi(14)x~0 = ~ri0(15)˙x~0 = v~0i(16)Начальные условия заданы в виде выражений ( 15) - ( 16), ~fi - равнодействующая всехсил, действующих на частицу:~fi = − 1mi∑ ~∇Uinj(17)j<i∑ j<i ~∇Ui j - потенциальная энергия системы.
Правая часть уравнения - сила, действующаяна данный атом со стороны всех остальных атомов. Взаимодействие между атомами являетсяпотенциальным, и поэтому сила записана как градиент потенциальной энергии системы.Потенциальную энергию системы можно представить в виде суммы вкладов от различныхтипов взаимодействий между атомами:U = Ur +Uφ +UKulon +UV dV(18)Ur - потенциальная энергия валентных связей, Uφ - валентных углов, UKulon - кулоновскихсил, UV dV - взаимодействий Ван-дер-Ваальса.
Энергия валентных взаимодействий и энергия колебаний валентных углов описывается параболическими потенциалами в моделинежёсткой воды. Графическое представление различных видов взаимодействия двух частицпредставлено на рис. 12.Для численного решения поставленной задачи (решения уравнений движения) возможноприменение различных численных алгоритмов:30Рис. 12: Зависимости Ui j от различных параметров• схема Эйлера2n~ = x~n + v~n dt + a~i (dt) + O(dt 3 )xin+1ii2(19)~ = v~n + a~n dt + O(dt 2 )vin+1ii(20)• алгоритм Верле.Разложим координаты i-ой частицы на (n + 1) и (n - 1) временных шагах в степеннойряд (Тейлора) по (dt) до 3-его порядка в точке 1,~vni ,~ani ~bni2, 6~ani (dt)2 ~bni (dt)3++ O(dt 4 )26~an (dt)2 ~bni (dt)3~xin−1 =~xin −~vni dt + i−+ O(dt 4 )26~xin+1 =~xin +~vni dt +(21)(22)⇓~ = 2~xn −~xn−1 +~an (dt)2 + O(dt 4 )xin+1iii(23)• метод Рунге-КуттыВ общем случае решается следующая задача:~y0 = f (x,~y)31(24)~y|x=x0 =~y0(25)~k1 = h f (xn ,~yn )(26)~~k2 = h f (xn + h ,~yn + k1 )22(27)~~k3 = h f (xn + h ,~yn + k2 )22(28)~~k4 = h f (xn + h ,~yn + k3 )22(29)⇓~yn+1 =~yn +~k1 + 2~k2 + 2~k3 +~k46(30)В нашем случае задачу можно переформулировать в виде системы дифференциальныхуравнений первого порядка:~x˙ = v(31)~v˙ = a(32)Каждое из уравнений решается предложенным методом независимо.2.2 Сравнение алгоритмовОсновным критерием оценки алгоритмов является время счёта при условии, что все вычисления проведены корректно и результаты при применении разных алгоритмов приведеныс одинаковой ошибкой.
Схема Эйлера ниже по порядку точности, следовательно, корректныевычисления требуют уменьшения порядка dt и увеличения времени счёта. Схемы Верле иРунге-Кутты одинаковы по порядку точности, но алгоритм Верле не учитывает скоростьявным образом (далее это различие будет подробно разобрано).322.3 Анализ стабильности получившейся структуры при увеличении температуры«Стабильность структуры» характеризуется следующими параметрами: сохранениемвсеми молекулами длин связей O-H и углов H-O-H, сохранением сетки водородных связей сзаданной точностью. Температура системы определяется следующим образом [64]:31kT N =∑ mi~v2i22 i=1,NT=∑i=1,N mi~v2i3kN(33)(34)где k - постоянная Больцмана, N - количество частиц.
В случае небольшого числа частицне так важно, в каком приближении это оценка корректна, температура является меройтеплового движения и свидетельствует о качественном поведении системы. В данном случаенас интересовало не только изменение координат частиц системы, но также изменение ихскоростей. Моделирование нагрева системы осуществлялось следующим образом: системав начальном положении (соответствующему минимуму энергии) релаксировала в течениенекоторого времени, затем скорости всех частиц одновременно увеличивались и снова системарелаксировала в течении некоторого длительного времени (аналог установления тепловогоравновесия).Разберём, в чём различие метода Рунге-Кутты и алгоритма Верле. При использованииалгоритма Верле рекуррентное выражение для координаты i-ой частицы выглядит так:~xin+1 = 2~xin −~xin−1 +~ani (dt)2 + O(dt 4 )(35)~vn+1=~vni +~ani dt + O(dt 2 )i(36)Пусть:В некоторый момент система получает некоторое количество энергии извне, что выражается следующим условием:~vn+1=~vn+1(1 + α),iiгде полгается, что α мало: α = 0.133(37)То есть, в данном примере частицы одновременно увеличивают все координаты скоростейна 10%.
Таким образом, в момент нагрева координата выражается следующим образом:~xin+1 =~xin+1 +~xin −~xin−1 +~ani dt 2+ O(dt 2 )1+α(38)Возможны и другие варианты моделирования нагревания системы: на каждом шагеувеличивать скорости (а, следовательно, и энергию), но незначительно, например, на 1%.В эти момент времени схема теряет 2 порядка точности, в то время как счёт по методуРунге-Кутты сохраняет 4-ый порядок точности. Схема Рунге-Кутты лучше подходит дляприменения в данном случае, поскольку она точнее.2.4 Описание кластеров, исследуемых в работеБыли рассмотрены следующие структуры молекул воды: T-кластер (рис.
13), L-L- кластер(рис. 14), фрагмент решётки льда (рис. 15), состоящий из 48 молекул, спираль 30/11 (рис. 16).В таблице 2 отображены их основные характеристики: количество молекул воды Nmolecules ,количество атомов Natoms , число водородных связей, NH−bonds , число дифференциальныхуравнений Neq , решаемых на каждом шаге.34Рис.
13: T-кластерРис. 14: L-L-кластер35Рис. 15: фрагмент решётки льдаРис. 16: Спираль 30/11Решётка льда, состоящая из 48 молекул, не является льдом в обычном представлении всвязи с малым количеством частиц. В этой структуре большая часть атомов находятся награнице, испытывая нескомпенсированное воздействие со стороны внутренних молекул. Этообщая особенность кластеров, содержащих малое число частиц.36Таблица 2: Характеристики параметрических структур и фрагмента решётки льдаСтруктураNmoleculesNatomsNH−bondsNeqT-кластер278140486L-L-кластер4012060720Лёд 48-кластер4814475864Спираль 30/1154162769722.5 Моделирование динамики водных кластеровРассмотрим динамику водных кластеров, начнём рассмотрение с T-кластера, так как онсильно отличается от других структур: данный кластер близок к сферической структурев отличие от других рассматриваемых кластеров.
Приведём некоторые количественныехарактеристики исследуемых водных кластеров после процесса минимизации (то есть система достигла локального минимума с заданной точностью) в табл. 3: общее количествомолекул воды в кластере Nmolecules , общее количество водородных связей NH−bonds , значениепотенциальной энергии UKul +UV dV , удельное значение потенциальной энергииUKul +UV dVNH−bonds:Таблица 3: Количественные характеристики водных кластеров после процесса минимизацииUKul +UV dVNH−bondsСтруктураNmoleculesNH−bondsUKul +UV dV , ккал/моль, ккал/мольT-кластер2739-257-5.8L-L-кластер4060-477-7.95Лёд 48-кластер4875-525-7Спираль 30/115476-566-7.45T-кластер наиболее приближен к шаровой форме, которая является фигурой с наименьшей площадью поверхности при одинаковом объёме, следовательно, доля граничныхатомов, на которые оказывается нескомпенсированное воздействие, в T-кластере минимально.
Зададим некоторый начальный уровень энергии, понаблюдаем за системой безнагрева, проверим консистентность модели. На рис. 17 представлен график зависимостиполной энергии T-кластера от времени для процесса свободной динамики при постояннойтемпературе. Полная энергия представляет собой постоянную функцию, что говорит онепротиворечивости модели.Было важно изучить поведение молекулярной системы при увеличении температуры.При нагревании амплитуды колебания частиц увеличиваются, и система отклоняется от37Рис. 17: T-кластер, график зависимости энергии при постоянной температуреравновесного положения.
Алгоритм нагрева был реализован следующим образом: черезкаждые 100 фс, за которые, как мы полагаем, система успевает релаксировать, скоростивсех частиц увеличивали на 1%. Наблюдали, как при нагревании перестраивается сетьводородных связей, а затем отдельные молекулы удаляются на бесконечность. Сделатьутверждение о значении температуры плавления кластера нельзя в виду малого количествачастиц и невозможности корректной оценки реальных температур.Введём понятие ”максимальный диаметр кластера” - это максимальное расстояниемежду атомами кислородов, принадлежащим водном кластеру, т.е.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
