Изучение равновесных свойств гибкого решеточного полиэлектролита методом Монте-Карло с помощью алгоритма Ванга-Ландау Текст научной статьи по специальности «Физика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 4 2007 Вып. 3

УДК 536.756, 538.953, 519.245

H.A. Волков, П.Н. Воронцов-Вельяминов, А.П. Любарцев

ИЗУЧЕНИЕ РАВНОВЕСНЫХ СВОЙСТВ ГИБКОГО РЕШЕТОЧНОГО ПОЛИЭЛЕКТРОЛИТА МЕТОДОМ МОНТЕ-КАРЛО С ПОМОЩЬЮ АЛГОРИТМА ВАНГА-ЛАНДАУ

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

Теоретическое описание гибких полиэлектролитов представляет собой сложную задачу. Существуют теории, описывающие незаряженные полимерные цепи, но они, как правило, неприемлемы для описания полиэлектролитов из-за присутствия дальнодействующих электростатических взаимодействий. В лучшем случае, взаимодействие между мономерами полииона описывается при помощи приближения Де-бая—Хюккеля без явного учета контрионов. Для сильно заряженных полиэлектролитов даже фундаментальные скейлинговые свойства остаются малоизученными. Таким образом, компьютерное моделирование становится важным инструментом для исследования подобных систем [1].

Число работ по компьютерному моделированию гибких полиэлектролитов при явном учете контрионов все еще достаточно ограничено (см., например [2-6]). В последние годы появились некоторые новые работы в этой области [7-9], в то же время остается множество вопросов, требующих дальнейшего рассмотрения. Существуют определенные трудности при получении низкотемпературных свойств полиэлектролитов, таких, например, как средний размер полииона. Это связано с тем, что эффективность стандартных методов моделирования ухудшается при усилении электростатических взаимодействий или при уменьшении эффективной температуры. Другое важное свойство—это свободная энергия, которая тесно связана со свойствами систем, наблюдаемыми экспериментально. Основные трудности ее вычисления определяются тем, что свободная энергия не может быть получена простым усреднением по конформациям системы, т. к. она является свойством всего статистического ансамбля. Поэтому используют специальные методы для получения свободной энергии.

Настоящая работа посвящена энтропическому моделированию (ЭМ) [10, 11] гибкого полиэлектролита с явным учетом мобильных ионов. Для реализации ЭМ-метода используется процедура самонастройки статвесов—алгоритм Ванга-Ландау (BJI) [12, 13]. В работе изучаются некоторые важные характеристики полиэлектролитов: зависимости от температуры и длины полииона внутренней и свободной энергии, а также среднего размера полииона. Хорошо известно, что в теории незаряженных полимеров такие фундаментальные свойства полимерных растворов слабо зависят от деталей модели. В нашей работе мы используем решеточную модель.

© H.A. Волков, П.Н. Воронцов-Вельяминов, А. П.Любарцев, 2007

Модель и метод моделирования. Рассматриваемая модель, предложенная ранее в работе [2], состоит из гибкого полииона, расположенного на ПК-решетке. Введены периодические граничные условия.

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

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

Объем основной ячейки для каждой решетки У=ЬЪ, где Ь—длина ребра ячейки. Рассматривались полиионы длины 10, 30, 50 и 80 звеньев (/V =11, 31, 51, 81) в ячейке с Ь = 30, 42, 50, 58 соответственно. Величина Ь выбирается таким образом, чтобы объемные доли контрионов и полииона с. = с V =4.07-10^ оставались постоянными.

. Взаимодействие мономеров полииона и мобильных ионов описывается кулонов-ским потенциалом У(г) = qfj.fr, где заряды = д > 0. Естественные единицы для нашей модели—это постоянная решетки а и заряд д. Таким образом, энергия измеряется в единицах с? I а, а температура—в ф / (акв). Помня определение эффективной плотности заряда полииона 4 как отношение Бьеррумовской длины (расстояния, на котором взаимодействие двух единичных зарядов равно их тепловой энергии) к расстоянию между зарядами, расположенными на полиионе, можно заключить, что 1 / Г для выбранной системы единиц [2].

Ключевая цель расчетов—получить для исследуемой системы плотность энергетических состояний которая позволяет вычислить путем простого интегрирования (суммирования) канонические средние в широком диапазоне температур. Мы используем тот же метод, что и наших предыдущих работах [14, 15], а именно, энт-ропическое моделирование [10, 11] с использованием алгоритма Ванга-Ландау [12, 13] с некоторыми модификациями.

Правильный учет дальнодействующих электростатических взаимодействий исключительно важен при моделировании полиэлектролитов. Общепринято, что метод Эвальда [16, 17]—это наиболее адекватный путь их учета. Мы сравнивали результаты, которые дает метод Эвальда, с данными, полученными при помощи метода минимального образа. Последний требует меньших затрат машинного времени и, в нашем случае, дает практически такие же результаты для плотности энергетических состояний при N = 11, 31 и рассмотренной объемной концентрации, как и метод Эвальда. Для более длинных цепей, моделируемых с большим количеством мобильных ионов, становится сильнее эффект экранирования, и предполагается, что отклонение этих распределений не должно быть больше. Поэтому в основной части наших расчетов при учете взаимодействий мы использовали метод минимального образа.

Так как 0.(Е), являясь целью моделирования, изначально неизвестна, для ее оценки используется итерационная процедура—ВЛ-алгоритм, который представляет собой самонастраивающуюся процедуру для получения 0.(Е) в ЭМ-методе. Вводится два набора ящиков: один из них соответствует конфигурациям системы

без самоналожений, а во втором (вспомогательном) учитываются конформации с самоналожениями. Для конфигураций без самоналожений интервал энергий Ёт.п<Е < Ётах поделен на пЕ равных интервалов (ящиков). В наших вычислениях число пЕ варьировалось от 100 до 1000. в зависимости от Ыр. ЁтЫ выбиралось для каждого N с таким расчетом, чтобы в процессе моделирования оставалось незадействованным некоторое количество ящиков с энергиями Е>ЁтЫ, близкими к Ё . Крайний «левый» ящик по шкале энергий, посещаемый системой при моделировании, соответствует минимальной энергии Е . > Ё ,, и значение Е . может рассматриваться как основное

г тт тт? т\п г г

состояние для данной модели.

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

Поскольку нам важно получить долю конфигураций без самоналожений Г20 (эффект исключенного объема), мы можем, в принципе, использовать всего один ящик для всех самопересекающихся конфигураций и получить в процессе вычислений его нормированный вес (1-П0). Однако, как это было показано в нашей предыдущей работе [14], такая процедура оказывается нестабильной. Поэтому необходимо ввести набор ящиков для конфигураций с самопересечениями. Каждый ящик из этого набора соответствует конфигурациям с определенным числом самопересечений полииона и наложений мобильных ионов. Число ящиков в данном наборе зависит от длины полииона и выбирается в диапазоне пт- 10-н 100. Отдельный ящик соответствует каждому возможному числу наложений меньше граничного значения (и < и у), и один дополнительный ящик (« = п^ используется для числа наложений, превосходящих его.

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

При реализации ВЛ-процедуры [12, 13] удобно перейти от числа состояний Ц к распределению по энтропии (как делалось в работах [14, 15]): Б. = 1пЦ (изначально все £ = 0). Вводятся два набора счетчиков: в одном накапливаются значения Б. для каждого ящика (счетчики энтропии), в другом—числа попаданий в каждый ящик п., которые в конце каждой серии вычислений дают нормированные вероятности посещений.

В процессе моделирования выполняются два типа шагов: изменяется конформация цепи или сдвигается мобильный ион. Вероятность каждого типа шага принималась равной 50%. Чтобы изменить положение выбранного случайным образом иона, мы делали его сдвиг на ±Д вдоль одной из осей ОХ, О У или ОХ, где Д= 1, 2, 3.

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

конформациями цепи велика, вероятность перехода обычно становится малой, что негативно сказывается на работе метода. Поэтому в большинстве наших расчетов I принималось равным 1 (как и в работе [2]), чтобы изменение конформации цепи на каждом шаге не было бы слишком большим. 2) Если два соседних звена цепи перпендикулярны друг другу, можно совершить так называемый шаг Г-типа. Направления этих двух звеньев меняются местами, т. е. новое направление первого звена принимается равным старому направлению второго, и наоборот. 3) Также совершается рептация цепи: отсекается звено с одного конца (который выбирается случайным образом) и пристраивается звено с другого конца цепи в случайно выбранном направлении. Шаги типа 1), 2) и 3) производились с вероятностями 50%, 40% и 10% соответственно.

Пробные шаги принимаются с вероятностью перехода [12, 13]:

Если это условие выполнено, пробное состояние V считается принятым. В противном случае система остается в старом состоянии /. Энтропия принятого состояния (5 или $г) увеличивается на величину А5, и соответствующий счетчик числа посещений (п. или увеличивается на 1.

Выбор вероятностей перехода согласно ВЛ-алгоритму (1) реализует метод эн-тропического моделирования [10, 11]. Метод состоит в следующем: совершая случайное блуждание в пространстве энергий с вероятностями перехода, пропорциональными обратной плотности состояний 1 / С1(Е), система равномерно посещает все энергетические ящики. Иными словами используя процедуру самонастройки ВЛ-алгоритма можно подобрать вероятности перехода такими, чтобы гистограмма посещений стала равномерной (плоской), и, следовательно, получить изначально неизвестную плотность состояний

После каждой серии шагов, согласно ВЛ-алгоритму [12, 13], значение Д£ уменьшают: Д5 —у сД£, где 0<с<1 (в работе [12] с = 0,5 и начальное значение Д5= 1). Количество шагов в серии должно выбираться таким, чтобы гистограмма посещений на каждой серии оставалась равномерной. После нескольких серий распределение 5. сформировано и настроено во всем интервале энергий.

В нашей работе множитель с = 0,9; вследствие этого, требовалось около 120 серий для достижения значений Д5, достаточно малых для точной настройки распределения П(Е). Число шагов в первой серии выбиралось равным 10б. В каждой последующей серии число шагов увеличивалось в 1,03 раза по сравнению с предыдущей. Такое число шагов оказывалось достаточным для достижения равномерности гистограммы посещений в конце каждой серии.

При помощи ВЛ-процедуры, описанной выше, мы получаем долю конфигураций без наложений и нормированную плотность энергетических состояний &(Е) = ехр(Я/Г.)). Распределение С1(Е) = Ц используется для вычисления канонических средних в широком диапазоне температур согласно стандартным соотношениям (2)—(6).

В каноническом ансамбле для внутренней энергии мы имеем:

->/") = тт[1,ехр(5;-5г)].

(1)

<Е> (Т)

(2)

Аналогично вычисляется <Е2>, давая в результате зависимость теплоемкости от температуры:

С(Т)=Т-2(<Е2>(Т)-(<Е>(Т))2). (3)

В выражении для энергии (2) нормировка Ц не важна, поскольку £2 стоит как в числителе, так и в знаменателе, и постоянный множитель сокращается. Для вычисления свободной энергии мы должны использовать полное число состояний, т. е. значение \У р.£1р равное числу конформаций без самоналожений, имеющих энергию Е.. 1¥рН—это полное число конформаций системы, Ц,-—доля конформаций без самоналожений относительно Ц—плотность энергетических состояний, нормированная на 1. Таким образом, для свободной энергии мы получаем:

ЕехР

V '

к Ту

(4)

где Ррк = -Пп Швк—свободная энергия фантомной системы; А/^= -Пп Ц,—избыточная

■Ри

свободная энергия в атермическом случае;

Д7Г = -Г1п

ЕехР

V I

а

(5)

избыточная каноническая свободная энергия (исчезает при отсутствии взаимодействий).

Аналогично выражается и каноническая энтропия 5(7) = (Е(Т) - Р{Т)) / Т.

Величина Ш ь для рассматриваемой системы вычисляется аналитически. Значение Ц,, доля конформаций без самоналожений для данной системы, вычисляется в процессе моделирования (также см. [14]). С другой стороны, его можно получить при помощи соотношения скейлинга [18, 19] для числа блужданий на решетке без самопересечений. Согласно нашим расчетам, Ц, = 0,146, 0,00123, 9,19-Ю"6, 6,1 МО"9 для N = 11, 31, 51, 81.

В процессе моделирования мы вычисляем среднии квадрат расстояния между концами полииона </?2>( для каждого энергетического интервала и получаем канонические средние стандартным образом:

>, ехр Ц

<Я2>(Т) = -

■=« Я

2^ехр

(6)

£2,

Результаты. Основной результат наших вычислений—плотность состояний 0.(Е) для N = 11, 31, 51, 81—представлен на рис. 1. Видно, что это гладкие кривые, которые опускаются до очень малых значений в области низких энергий. Для N =81 вероятность того, что система попадет в ящик с минимальной энергией, составляет приблизительно 10~285. Таким образом, ВЛ-алгоритм позволяет вычислять вероятности крайне редких состояний.

Поскольку мы используем шкалу удельных энергий Е / N, на рис. 1 становится видно, что значения минимальных удельных энергий для разных длин полииона близки по величине: Е . Ш =-1,754, -1,840, -1,809, -1,695 для # = 11, 31, 51, 81 соот-

т}П р 7 7 7 7 7 7 7 ^ р 7 7 7

ветственно. Эти значения, по-видимому, близки к наименьшим значениям, возможным в принципе для рассматриваемой модели (основным состояниям).

E/N

p

Рис. 1. Нормированные распределения Q(E) в зависимости от удельной энергии EIN для полиэлектролита с Л^ =11 (ромбы), 31 (треугольники), 51 (квадраты), 81 (круги).

На рис. 2 представлено сравнение зависимостей AF(T) и Е(Т), полученных согласно (2) и (5). Можно заметить, что Е{Т) и AF(T) имеют одинаковые горизонтальные асимптоты при Г-^О и при Т —*■ оо для каждого N. Это объясняется поведением внутренней (2) и свободной (5) энергий в этих предельных случаях.

Интересно отметить, что кривые AF, соответствующие разным N, пересекаются в одной точке. Это также справедливо и для зависимостей Е. Для АF эта точка, примерно, соответствует AF ~ О, Т~ 0,2. Иными словами избыточная каноническая свободная энергия при данной температуре практически не зависит от длины полии-она. Заметим также, что высокие пики на графиках теплоемкости (рис. 3) возникают, приблизительно, при той же температуре.

Мы также используем метод Метрополиса для вычисления Е(Т) для // = 11 и N=81 для нескольких температур и приводим сравнение с соответствующими BJI-результатами. Хорошее соответствие результатов, полученных двумя разными методами, наблюдается в области температур Т- 0,05 1000 для // = 11 (не показано на рисунке) и в области Т= 0,2 1000 для N = 81.

На рис. 3 представлены температурные зависимости для удельной теплоемкости. Для коротких цепей iV = 11, 31, кривые С(Т) имеют множество небольших максимумов, а для более длинных, N =51, 81, появляется высокий пик. Этот узкий пик может быть рассмотрен как свидетельство перехода полииона в глобулярное состояние. Зависимости С(Т) стремятся к нулю при высоких и низких температурах. Сложный рельеф кривых, имеющих несколько максимумов, по-видимому, вызван перестройкой конформаций цепи и мобильных ионов в соответствии с разными температурными условиями. Наши С(7)-кривые для N = 51, 81 (рис.3) качественно имеют тот же характер, что и кривые, полученные в работе [20] для решеточной модели незаряженного полимера. Конечномерный скейлинг для максимумов С{Т) зависимостей,

Рис. 2. Внутренняя энергия Е (пунктир) и избыточная каноническая свободная энергия Д/7 (сплошные линии) как функции от Т. Квадраты соответствуют данным, полученным методом

Метрополиса для N = 81.

С

1000-800 600

- 1 N =81 • р

/ ........

0.1 1 т 1 1 • 1 10

Рис. 3. Удельные теплоемкости как функции от Г для Л^ = 11 (точки), 31 (пунктир), 51 (сплошная линия), 81 (штрих-пунктир).

Рис. 4. Зависимость среднего расстояния между концами цепи \< Я2 > от температуры

для N = 11 (точки), 31 (пунктир), 51 (сплошная линия), 81 (штрих-пунктир). Крестики соответствуют данным, полученным методом Метрополиса для N =11,81.

проведенный в работе [20], указывает на две стадии перехода для цепей конечной длины: переход «клубок—глобула» и кристаллизация. К сожалению, наших данных недостаточно для проведения анализа, проделанного в работе [20], и необходимы данные для более длинных цепей.

Зависимость среднего расстояния между концами цепи > от температуры представлена на рис. 4. Для каждого N кривые имеют максимумы в окрестности 1. Этот факт согласуется с результатами, полученными методом Метрополиса в более ранней работе [2], а также с результатами, представленными другими авторами [7].

Видно, что максимум кривой

>(7) смещается в область более высоких температур при увеличении длины полииона. Для А^ =11, 31, 51, 81 соответствующие температуры Т0 = 0,92, 1,18, 1,28, 1,36, и максимальные значения \1<Я2>(Т0) = 5,60, 15,94, 26,75, 40,38. Таким образом, при Т~ 1 полиион растянут больше, чем при других

температурах. При высоких температурах у1<Л2> стремится к значению расстояния между концами незаряженного полимера без самопересечений =а№, где N=^-1 (согласно нашим вычислениям [14], а= 1,025). На рис. 4 видно, что полученные кривые асимптотически стремятся к этим уровням.

При Т< 1 конденсация контрионов становится очень сильной, что приводит к резкому уменьшению размера цепи (переход в глобулу). Для N =51,81 этот переход

имеет вид резкого скачка при температуре Т ~ 0,1 ^ 0,2. Наряду с результатами для удельной теплоемкости (рис. 3), это может служить свидетельством фазового перехода. Сравнение с данными, полученными стандартным МК-методом (метод Метрополиса), показывает полное согласие обоих подходов для случая N = 11 во всем диапазоне рассмотренных температур. Для наиболее длинного полииона, А^ — 81, хорошее согласие данных наблюдается в широком диапазоне температур. Только в области коллапса, 0,2 < Г <0,5, имеется заметное отличие результатов, полученных разными методами, хотя качественно результаты согласуются.

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

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

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

Данная работа была проведена при поддержке Российского Фонда Фундаментальных Исследований (грант РФФИ 05-02-17428) и Шведской Королевской Академии Наук.

Summary

VolkovN.A., Vorontsov-Velyaminov P.N., Lyubartsev А. P. Study of equilibrium properties of the flexible lattice polyelectrolyte by Monte Carlo method within Wang-Landau algorithm.

We consider a model of the flexible polyelectrolyte model on a lattice. The Coulomb potential and the excluded volume condition between different ions/beads are taken into account. An efficient Monte Carlo method within the numeric algorithm of Wang and Landau is used. The obtained energy distributions provide further calculation of canonical properties such as conformational energy, heat capacity, entropy, free energy in a wide temperature range by simple integration. Wang-Landau algorithm allows us to obtain statistics for the states with extremely low probabilities. We also test our results comparing them with the data obtained by standard Monte Carlo method.

Литература

1. Lyubartsev A. P., Nordenskiold L. Computer Simulation of Polyelectrolytes//Handbook of Polyelectrolytes and Their Applications. 2002. Ch. 11. P. 309-325. 2. Любарцев А.П., Воронцов-Вельяминов П.H. //Высокомол. соедин. 1990. Vol. 32. Р. 721. 3. Severin М. Hi. Chem. Phys. 1993. Vol. 99. P. 628. 4. Stevens M. J., KremerK. Hi. Chem. Phys. 1995. Vol. 103. P. 1669. 5. Stevens M. J., Plimpton S. J. // Eur. Phys. J. 1998. Vol. 2. P. 341. 6. Hsiao P-Y. II J. Chem. Phys. 2006. Vol. 124. P. 044904. 7. KlosJ., and Pakula Т. И J. Chem. Phys. 2004. Vol. 120. P. 2496. 8. J. Klos, and T. Pakulall J. Chem.

Phys. 2004. Vol. 120. P. 2502. 9. KlosJ., PakulaT.IIJ. Chem. Phys. 2005. Vol. 122. P. 134908. 10. BergB.A., Neuhaus T. //Phys. Rev. Lett. 1992. Vol. 68. P. 9. 11. Lee J. //Phys. Rev. Lett. 1993. Vol. 71. P. 211. 12. Wang F., Landau D.P.H Phys. Rev. Lett. 2001. Vol.86. P. 2050. 13. Wang F., Landau D. P. II Phys. Rev. 2001. Vol. 64. P. 056101. 14. Vorontsov-Velyaminov P.N., VolkovN.A., Yurchenko A.A. //Journ. Phys. A: Math. Gen. 2004. Vol. 37. P. 1573. 15. VolkovN.A., Yurchen-koA.A., LyubartsevA.P, Vorontsov-Velyaminov P.N. //Macromol. Theory Simul. 2005. Vol. 14. P. 491. 16. Ewald P. II Asm. Phys. 1921. Vol. 64. P. 253. 17. Allen M. P., Tildesley D.J. //Computer Simulations of Liquids, 2nd ed. Clarendon. Oxford, 1987. 18. Douglas J., Guttman C.M., MahA., Ishina-be T. IIPhys. Rev. 1997. Vol. 55. P. 738. 19. GennesP.G. dell Scaling Concepts in Polymer Science, Cornell University Press. London, 1979. 20. RampfF., Paul W., Binder K. II Europhys. Lett. 2005. Vol. 70. P. 628.

Статья принята к печати 20 февраля 2007 г.