Диссертация (1097582), страница 14
Текст из файла (страница 14)
Так вот, этоусловие выполняется в алгоритме Ванга-Ландау только в предельном случае, когдамножитель f стремится к 1. По этой причине конечное значение ffinal для множителя fвыбирается достаточно близким к 1 (например, с точностью до седьмого знака послезапятой, то есть конечное значение множителя f равно 1 в пределах точности обычногодействительного числа на компьютере).
Таким образом, производить набор статистикидля вычисляемых в ходе компьютерного эксперимента физических величин имеет смыслтолько после окончания итерационной процедуры для расчета функции плотностисостояний, то есть когда моделирование проводится уже с неизменяемой функцией g (втом числе, в пункте 2). Для алгоритма ВЛ по внешнему параметру условие детальногобаланса должно фактически выполняться только при одном значении внешнего параметрарасширения (при том, который соответствует реальной физической системе), а это должен53обеспечить «основной» алгоритм МК для данного расширенного ансамбля (лежащий в егооснове еще до процедуры расширения).Второе замечание – чисто техническое.
При расчете функции g необходимопроизводить все операции с действительными числами с двойной точностью, а такженакапливать не саму весовую функцию g, а ее логарифм, ибо иначе слишком быстроможно достигнуть максимально допустимого (на данном компьютере и для данногокомпилятора) значения действительного числа. Кстати, именно поэтому начальноезначение множителя f для весовой функции выбирается равным f0 = e (на первом шагеитерационной процедуры к логарифму g прибавляется 1 при каждом посещении данногозначения величины A).Третье замечание касается разбиения полного интервала исследуемых значенийпараметра на «окна», или подынтервалы.
Если интервал доступных значений энергиислишком широк, приходится считать функцию плотности состояний g(E) отдельно внескольких интервалах во всем диапазоне возможных значений полной энергии. Послеэтого производится «сшивка» функции плотности состояний и «сшивка» гистограммэнергии и других наблюдаемых физических величин (радиуса инерции, параметрапорядка и т.д.) для всех интервалов. Далее пересчет температурных зависимостей игистограмм физических величин для канонического ансамбля производится стандартнымспособом.Четвертое замечание: могут возникнуть проблемы со сходимостью алгоритма ВЛ,особенно, в сложных полимерных системах.Пятое замечание: в литературе есть примеры использования двумерного алгоритмаВЛ, когда строится функция плотности состояний по двум параметрам [27], однако такаяпроцедура требует намного больше компьютерного времени.2.4.
Методы расчета давления в решеточных моделях Монте-Карло для расчетауравнения состояния полимерного раствораВ этом разделе изложены методы измерения давления в компьютерноммоделировании с помощью решеточных алгоритмов Монте-Карло, часть из которых былавпервые разработана в рамках данной диссертации. Было также детально исследовановлияния конечных размеров системы на величину рассчитанного давления (результатыподробно изложены в работах [193] и [233].Проблема точного вычисления давления в компьютерном моделированииполимерных растворов и расплавов является очень важной при изучении фазовых54переходовпервогорода(например,нематическогоупорядоченияврастворежесткоцепных макромолекул).
Температура перехода и соответствующие плотности фаз(изотропной и нематической) в точке перехода могут быть определены из условияравенства химических потенциалов и давлений в разных фазах. Так, например, в работах[188, 190] изучался переход изотроп-нематик в растворах жесткоцепных макромолекул спотенциалом притяжения между мономерными звеньями в рамках решеточной модели ибольшого канонического ансамбля.
Важный методический аспект этих работ связан сметодом определения перехода и фазовой диаграммы. Из трех условий равновесия [191]между фазами 1 и 2: T1 = T2 , µ1 = µ2 , P1 = P2 (где T – температура, µ - химическийпотенциал, P – давление), первые два выполняются автоматически в большомканоническом ансамбле. Задача, таким образом, сводится к нахождению значенияхимического потенциала µ ∗ , µmin < µ ∗ < µmax , при котором значения давления в обеих( )( )фазах равны, т.е.
P1 µ ∗ = P2 µ ∗ . Таким образом, возникает необходимость измерениядавления в полимерных растворах различной плотности, т.е. уравнения состояниярастворов макромолекул.Рассчитать давление при моделировании методами МД и МК с использованиемконтинуальных моделей (и вычислить уравнение состояния, т.е. зависимость давления отплотности) можно несколькими способами: 1) при помощи использования теоремывириала, 2) путем моделирования в изотермо-изобарическом ансамбле, 3) методомвстраивания пробной частицы (цепи), 4) на основе межмолекулярных корреляционныхфункций (детальное описание этих методов можно найти в [35,192]).
Например, согласнотеореме вириала [192], давление в каноническом ансамбле может быть получено какP = ρk BT +1dV∑ f (r ) ⋅ rijij,(21)i< jгде ρ=N/V – плотность, N – число частиц, V – объем, d – размерность пространства,f ( rij ) - сила, действующая между частицами i и j, находящимися на расстоянии rij , суммаберется по всем парам частиц, а усреднение проводится по полному временимоделирования после достижения системой равновесия. Этот метод неприменим врешеточных моделях в силу их дискретной природы.Детальная разработка методов расчета давления в полимерных растворах длярешеточных моделей МК (включая сравнение различных методик) была выполнена вработах [193, 205].
Эти исследования позволили изучить зависимость уравнениясостояния раствора жесткоцепных макромолекул от жесткости цепей, а также изучить55поверхности раздела изотропной и нематической фаз и эффекты смачивания плоскойповерхности нематическим раствором полимера [205].В случае решеточных моделей наиболее часто используются метод встраиванияпробной цепи (МВПЦ, test chain insertion method, TCIM) [194, 195] и методотталкивающей поверхности (МОП, repulsive wall method, RWM) [196-200].
Используя этиметоды для расчета давления при разных плотностях системы, мы можем получитьуравнение состояния раствора макромолекул. В конце 80-х годов Р.Дикманом ссоавторами был опубликован ряд работ [194-204], посвященных вычислению давления ватермических полимерных растворах (где присутствуют взаимодействия, обусловленныетолько исключенным объемом). Были представлены теоретические расчеты и результатыкомпьютерного моделирования методом МК в каноническом ансамбле для решеточной иконтинуальной моделей в двух- и трехмерном пространстве с помощью обоих этихметодов.2.4.1. Метод термодинамического интегрирования.Метод вставки пробной цепи (рис.6) является методом термодинамическогоинтегрирования, при котором требуется определить вероятность вставки пробной цепи враствор с объемной долей полимера φ . Осмотическое давление π (φ ) ≡Pдля системыk BTn-меров рассчитывается затем по формулеφφ1π (φ ) = [1 − ln p(φ , n)] + ∫ ln p(φ ", n)dφ " ,nn0(22)где p(φ , n) есть вероятность вставить n-мер в раствор n-меров с объемной долей φ .Следует отметить, что метод вставки пробной цепи трудоемок и применим лишь дляразбавленных и полуразбавленных растворов и не слишком длинных цепей, т.к.
прибольших длинах цепей вероятность вставки n-мера становится ничтожно малой. Этоограничение можно обойти, если использовать алгоритм МК с конформационнымсмещением выборки (CBMC, [28-34]).В работе [205] описан метод термодинамического интегрирования, который такжеможет применяться для вычисления уравнения состояния полимерного раствора вкомпьютерном эксперименте. Идея метода состоит в том, что моделирование проводитсяв большом каноническом ансамбле при различных значениях химического потенциала, тоесть пробная цепь действительно добавляется или, наоборот, удаляется из системы при56выполнении определенных условий на принятие пробного шага добавления или удаленияцепи.2.4.2. Метод отталкивающей поверхностиВ то же время метод отталкивающей поверхности (МОП) остается эффективным ипри больших плотностях растворов, а также в расплавах.
В ячейке моделированияразмером L × L × Hвводится отталкивающий потенциал λ = exp{−ε1 / kBT}, 0 ≤ λ ≤ 1,действующий на частицы, находящиеся в пристеночном слое z = H . Тогда осмотическоедавление может быть получено по формуле1π (φ ) = ∫0dλλφ H (λ ) ,(23)где φ H (λ ) - доля занятых ячеек в слое z = H при данном потенциале отталкиванияот стенки. Таким образом, этот метод предполагает серию запусков программы дляразличных значений отталкивающего потенциала для последующего интегрирования ирасчета давления раствора макромолекул с данным значением плотности.Отметим, что при использовании любого из вышеописанных методов для того,чтобы получить уравнение состояния раствора макромолекул (т.е.
зависимость давленияот плотности) необходимо проводить серии запусков программы для ряда значенийплотности системы. Это, в свою очередь, приводит к большим затратам по времени. Болеетого, в работе [193] было показано, что методу отталкивающей стенки присущизначительные эффекты конечного размера (в особенности при больших значенияхплотности), которых можно избежать, лишь проводя моделирование в большомканоническом ансамбле (более подробно этот материал изложен в разделе 4.1.4 ниже).2.4.3.