CC BY
37
УДК 536.756, 538.953, 519.245
Вестник СПбГУ. Сер. 4. 2011. Вып. 4
И. А. Силантьева, П. Н. Воронцов-Вельяминов
ИССЛЕДОВАНИЕ РЕШЁТОЧНЫХ МОДЕЛEЙ ПОЛИМЕРНЫХ ЦЕПЕЙ И ЗВЁЗД МЕТОДОМ МОНТЕ-КАРЛО С ИСПОЛЬЗОВАНИЕМ АЛГОРИТМА ВАНГА—ЛАНДАУ
Введение. В последнее время полимерные звёзды являются предметом многочисленных теоретических и экспериментальных исследований [1—3]. Это связано с возможностью использования данных молекул для доставки лекарственных веществ и ДНК в клетки [4, 5]. Полимерные звёзды — простейший пример модели разветвлённых макромолекул. Поэтому знание их равновесных свойств и структуры может быть важно для исследования макромолекул более сложной архитектуры.
Целью представляемой работы является применение метода Монте-Карло (МК) с использованием алгоритма Ванга—Ландау (В-Л) для вычисления плотности состояний полимерных цепей и звёзд. Полученные распределения дают возможность рассмотреть в широком диапазоне температур равновесные свойства исследуемых макромолекул, такие как внутренняя энергия, энтропия, теплоёмкость. В поле зрения оказались цепи длиной N ^ 120 мономеров и звёзды с шестью лучами длиной Nл = 5,12, 20 сегментов каждый (общее число сегментов в звезде N = 30, 72,120). Аналогичная задача была успешно решена для полимерных цепей и колец в работах [6, 7]. Исследование полимерных цепей может быть полезно, несмотря на то, что они подробно исследованы ранее. Во-первых, сравнение результатов даёт уверенность в правильности работы алгоритма, во-вторых, возможность сопоставить результаты для цепей и для звёзд.
В атермическом случае рассчитана доля самонепересекающихся конформаций исследуемых молекул. В термическом случае отбирались самонепересекающиеся конфор-мации и вычислялись распределения по числу контактов между мономерами.
Метод. Для выполнения работы использовался метод энтропического моделирования [8, 9] (одна из модификаций метода Монте-Карло) по алгоритму В-Л, который был предложен в 2001 г. [10] и постепенно нашёл широкое применение для исследования полимерных и других систем. Алгоритм В-Л является реализацией метода энтропического моделирования. Он решает проблему подбора подходящих вероятностей перехода для получения требуемого при энтропическом моделировании равномерного посещения энергетических состояний и, соответственно, позволяет получить плотность состояний П(Е). В общем случае диапазон изменения энергии системы Ет\п ^ Е ^ Етах разбивается на конечное число (N5) равных отрезков («ящиков»). Заводится массив плотности состояний П, состоящий из (N5) элементов, каждый из которых соответствует отрезку разбиения энергии. Изначально все элементы П берутся равными единице. Процедуру можно сделать более удобной для машинного счёта, если перейти к энтропии состояний Б(Е) = 1пП(Е) (начальные состояния Б^ = 0). В процессе вычислительного эксперимента на каждом МК-шаге происходит изменение конформации системы. Пусть Е^ и Е^ — энергии системы до и после изменения. Каждая из них попадает в свой «ящик» — г-й и _7-й соответственно (номера г и ] могут совпадать). В таком случае изменения в системе принимаются с вероятностью
р(Е, ^ Ез) = шш(1, ) = тт(1,е^"^). В случае отказа система возвращается в исходное состояние.
© И.А.Силантьева, П.Н.Воронцов-Вельяминов, 2011
После принятия или непринятия новой конформации системы всё повторяется на новом МК-шаге. Каждый раз при посещении к-го «ящика» (в случае принятия изменений системы к = ], при отказе к = г) проводится изменение к-го элемента массива Б: к Б к добавляется слагаемое ДБ. Ряд таких элементарных шагов составляет серию. На протяжении серии МК-шагов величина ДБ остаётся неизменной. На каждой последующей серии ДБ уменьшается: ДБ ^ сДБ с инкрементом 0 < с < 1. Одновременно с массивом Б заводится массив посещений п, элементы которого изначально равны нулю. На каждом МК-шаге в ячейку пк, соответствующую посещению к-го «ящика», добавляется единица. Во время счёта необходимо следить, чтобы гистограмма посещений оставалась равномерной с достаточной степенью точности. В результате использования этого алгоритма происходит автоматическая настройка весов вероятности перехода, которые одновременно являются плотностями состояний. По окончании расчёта вычисляется массив 0.(Е) = ехр(Б(Е)) и нормируется на единицу. В качестве величины, по которой мы находим распределение, может выступать не обязательно энергия, а любая величина, зависящая от конформации системы, например число самопересечений.
Отметим, что за последние примерно пять лет резко возросло число работ, в которых данный алгоритм использовался для изучения полимеров [11—14]. Особенно эффективен он для исследования фазовых переходов, так как позволяет получить характеристики в широком диапазоне температур и при этом функция плотности состояний рассчитывается только один раз.
Модель. В рамках решёточной модели полимерная цепь представляется как случайное блуждание на 3-мерной решётке. Различают фантомное и полуфантомное блуждания [7]. В нашей работе использовалось полуфантомное блуждание. Преимущество полуфантомного блуждания состоит в том, что доля самонепересекающихся конфор-маций среди полуфантомных существенно больше, чем среди фантомных. Под цепью длиной N будем понимать цепь, имеющую N связей и N +1 мономер. Полимерная звезда из шести лучей представляет собой 6 полимерных цепей длиной Nл, закреплённых одним концом в центре. В нашем моделировании полимер описывается набором координат (хг, уг, хг), где г = [0,Ж]. В случае цепи её начальная точка и первый сегмент неподвижны, в случае звезды — центр и первый сегмент каждого луча звезды остаются неподвижными. Остальные сегменты генерируются с помощью случайного полуфантомного блуждания на решётке.
Чтобы изменить конформацию цепи, мы выбираем одну из точек от 1 до N — I (точку ко) и строим новый кусок полуфантомной цепи между точками ко и ко + I [6], где I — число сегментов в модифицируемом куске длиной [N/20]—[N/10]. Проверяется, не пересекается ли этот участок сам с собой и с участком цепи до точки ко. Оставшийся участок цепи подвергаем параллельному переносу, учитывая, что после сдвига цепь должна оставаться полуфантомной. Далее проверяем, нет ли самопересечений в модифицированной цепи. Если есть, то модифицируем заново тот же самый кусок, пока не получим конформацию без самопересечений. Таким образом, для отбора самонепересекающихся траекторий мы используем алгоритм безусловных случайных блужданий.
Для изменения конформации звезды сначала выбирался один из шести лучей, затем выбирался сегмент ко этого луча и производилась его модификация подобно тому, как это было проделано для цепи. В термическом случае для упрощения алгоритма луч модифицировался от сегмента ко + 1 до конца. После изменения конформации полимера вычисляем число контактов в полученной молекуле и производим процедуру согласно алгоритму В-Л. В большинстве случаев вычисления включают 30-60 серий (для цепи 30 серий) с начальными значениями Бг = 0, ДБг = 0,001 + 0,1 и инкрементом 0,9. Длина
серии зависит от длины цепи и составляет 1-10 млн шагов для N = 12 + 120. Время каждого расчёта составляло не более 15 ч на стандартном 4-ядерном процессоре.
Результаты.
Атермический случай. В атермическом случае взаимодействие сводится к запрету на самопересечение мономеров. Целью является вычисление доли самонепересекающихся конформаций в зависимости от числа N мономеров в молекуле. Для отбора самонепересекающихся конформаций мы использовали сначала два «ящика»: конформации без самопересечений и все остальные. Таким образом удалось получить долю самонепересекающихся конфромаций для цепей длиной N ^ 100 с погрешностью менее 1 %, а с использованием пяти «ящиков» (т. е. самопересекающиеся конформации складывались в несколько «ящиков») — до N ^ 300 с погрешностью 1 %. Погрешность вычислялась по отношению к результатам, рассчитанным c использованием формулы скейлинга. Известно, что число самонепересекающихся конфигураций для цепи длиной N на решётке определяется соотношением скейлинга [6]
CN = NY-1.
Для простой кубической решётки A = 1,2079, константа связи ц = 4,684043, восприимчивость у = 1,1568. Полное число конформаций для полуфантомной цепи длиной N есть z(z — 1)N-1, где z — число ближайших соседей на решётке. Для трёхмерной ПК решёти z = 6. Тогда доля самонепересекающихся конформаций для полуфантомной цепи
Q — -V4^1 (1)
Стоит отметить, что долю самонепересекающихся траекторий для свободной цепи можно вычислять помимо алгоритма В-Л, например методом ветвления траекторий [15] или используя Pivot-алгоритм [16].
Рассмотрим полимерную звезду. Если длина луча N < 10 звеньев (т. е. менее 60 звеньев в звезде), то долю самонепересекающихся конформаций Qo можно найти по алгоритму В-Л, используя только 2 «ящика»: для конформаций с пересечениями и без. Для N > 10 двух «ящиков» становится недостаточно. В этом случае используется следующий способ. С помощью алгоритма безусловного блуждания определяется граница нь диапазона [0, нь], в который попадает не менее 99 % событий от всего распределения по самопересечениям. Далее на отрезке [0,нь] проводится процедура В-Л (100 серий, по 0,5-1 млн шагов в серии, начальное значение AS = 0,1, c = 0,9) и, таким образом, определяется Qo. Конформации с числом самопересечений, превышающим нь, не учитываются (доля таких случаев < 3 %), полагается, что они не оказывают существенного влияния на конечный результат.
Данный алгоритм был протестирован в ситуации, когда не учитывалось взаимопересечение лучей. Результат сравнивался с формулой Q0 = й0л, где О0л — доля самонепересекающихся конформаций для цепи длиной Nл, вычисленная по формуле (1). Данные представлены в табл. 1.
Как видно из табл. 1, данные математического эксперимента хорошо согласуются с формулой скейлинга. Это даёт уверенность в том, что программу можно использовать для дальнейших исследований модели звезды. Затем было включено атермиче-ское взаимодействие между лучами (учёт взаимных пересечений) и рассчитана доля
0 г
-0,1 -0,2 -0,3 -0,4 -0,5 -0,6 -0,7
0
0,02
0,04 0,06
1/N
0,08
0,1
Рис. 1. Удельная избыточная энтропия (относительно энтропии фантомной цепи) в зависимости от обратного числа сегментов для звезды из шести лучей (закрашенные точки) — данные, полученные по алгоритму В-Л, аппроксимированные по методу наименьших квадратов (сплошная кривая); для цепи (длинный пунктир) — зависимость, рассчитанная по формуле (1); для колец [7] (короткий пунктир) — приведена для сравнения:
крестики — данные для одного луча звезды (алгоритм В-Л) при условии, что лучи не взаимодействуют между собой; Пол = (По)1/6; точки хорошо ложатся на кривую, построенную по формуле скейлинга; предельные значения при N ^ те, —0,2476 для свободной цепи, —0,2494 для звезды, —0,2495 для колец
самонепересекающихся конформаций по описанному алгоритму для звёзд с суммарной длиной сегментов N ^ 720 ^ 120). Зная долю самонепересекающихся конформаций, можно получить зависимость удельной избыточной энтропии ДS/N от 1/N (1). Полное число конформаций полуфантомной звезды с шестью лучами есть г\(г — 1)№-1)6 =
Таблица 1
Доля самонепересекающихся конформаций для полимерной
звезды при отсутствии взаимодействия между лучами
= z!(z - 1)
N-6
Тогда энтропия звезды в атермиче-
ском случае может быть сведена к формуле
По, В-Л пъ По, скейлинг
12 9,74Е-02 30 9,78Е-02
30 2,010Е-04 60 2,01Е-04
50 1,326Е-07 90 1,323Е-07
Б = - 1Г~ЬПП) = In z +(N- 1) In - + In
6
5 (z - 2) !
(z - 1)4
ln Пп
где N — полное число сегментов в звезде; По — доля самонепересекающихся конформаций среди полуфантомных. Первое слагаемое — энтропия фантомной цепи длиной N, второе слагаемое — энтропия полуфантомной цепи длиной N относительно фантомной цепи, третье слагаемое — добавка, связанная с тем, что первые 6 сегментов закреплены в центре. Удельная избыточная энтропия звезды относительно энтропии фантомной цепи
ДБ Л 1 \ , 5 1,0* — 2)! 1 — = ( 1 " - ) 1п - + - 1п
N
6 N (z - 1)4
+ TT In ПС
N
На рис. 1 представлена зависимость AS/N от 1/N для звезды, полученная в результате численного эксперимента, для цепей — построенная на основе формулы скейлинга
и колец — из работы [7]. Зависимость для звезды мы аппроксимировали по методу наименьших квадратов функцией / (х) = Ах 1п(х) + Вх + С, х = N-1. И получили для коэффициентов следующие значения: А = 0,7399, В = —1,4837, С = -0,2494. Предельное значение удельной избыточной энтропии при N ^ ж для звёзд С = —0,2494 с достаточной точностью совпадает с аналогичным значением для свободных цепей С8с = —0,2476 (скейлинг) и колец Ст-те = —0,2495. Таким образом, можно утверждать, что при N ^ ж удельная избыточная энтропия для звёзд, колец и цепей стремится к одному и тому же пределу.
Термический случай. В термическом случае полагается, что каждой паре несоседних по цепи мономеров, находящихся в контакте, отвечает некоторая энергия (е > 0 или е < 0). Конфигурационная энергия полимера есть Е = ет. Задача состоит в том, чтобы найти распределение самонепересекающихся конформаций полимера по числу контактов т. Число контактов полимерной молекулы меняется в пределах [0,ттах]. Максимальное число контактов ттах для цепи оценивалось по формуле, приведённой в работе [17]:
ттах = т8р;га1 < — + 1)Лс + ¿С],
Таблица 2 Максимальное число контактов
для самонепересекающихся конформаций полимерной цепи, полученное в процессе моделирования по алгоритму В-Л и рассчитанное по формуле (2)
(2)
N пгтах, В-Л пгтах, расчёт
12 9 10
20 18 20
30 32 33
50 59 61
100 133 137
Таблица 3 Максимальное число контактов
для самонепересекающихся конформаций полимерной цепи, полученное методом ветвления и рассчитанное по формуле (2)
N гпщах, ветвление пгтах, расчёт
50 59 61
100 136 137
150 216 217
200 297 300
300 460 468
где N — длина цепи; с! =3 — размерность пространства; ас = с! — 1, Дс = (! — 1)/!. В табл. 2 приведены данные для ттах, полученные в численном эксперименте и по вышеприведённой формуле. Данные моделирования не превышают значений расчёта, а значит, не противоречат ей.
Максимальное число контактов для более длинных цепей было оценено независимым методом ветвления траекторий [15]. Результаты приведены в табл. 3. Формула (2) даёт немного завышенное значение относительно численного эксперимента. Как отмечено в работе [17], максимальное число контактов важно как техническая стартовая точка в тестировании моделей самонепересекающихся блужданий и так как оно определяет внутреннюю энергию компактных полимеров (например, белков). Несоответствие между ттах и численной оценкой ттах может служить критерием, насколько хорошо смоделированы эти компактные конфигурации.
Максимальное число контактов для
звезды меньше либо равно максимальному числу контактов для цепей. Так, например, для цепи длиной N = 30 ттах = 32, в то время как для звезды с тем же суммарным количеством сегментов (^л = 5) ттах = 29. Разница связана с ограничениями, наложенными на сегменты в звезде, — пришивка лучей в общей точке. На рис. 2 приведено распределение по числу контактов для цепей длиной 12, 30, 50, 100 сегментов. Можно отметить, что вероятность генерации самонепересекающейся конформации со значениями т, близкими к ттах, достаточно мала
1,Е+00 1,Е-02 1,Е-04 1,Е-06 1,Е-08 1,Е-10 1,Е-12 1,Е-14 1,Е-16 1,Е-18 С? 1,Е-20 1,Е-22 1,Е-24 1,Е-26 1,Е-28 1,Е-30 1,Е-32 1,Е-34 1,Е-36 1,Е-38 1,Е-40
Рис. 2. Распределение самонепересекающихся конформаций по числу контактов для цепей длиной 12, 30, 50, 100 сегментов, полученное по алгоритму В-Л (закрашенные точки), и аналогичные данные из работы [6]: круги (М = 12), квадраты (М = 30), треугольники (М = 50), ромбы (М = 100)
0 10 20 30 40 50 60 70 80 90 100 110 120 130
т
Рис. 3. Диаграмма посещений, полученная при вычислении распределения по числу контактов для N = 30 (1), 50 (2), 100 (3):
для М = 12 диаграмма полностью плоская и здесь не приведена
0,04 0,035 0,03 0,025 - 0,02 0,015 0,01 0,005 0
10 30 50 70 90 110 130
т
и для N = 100 не удаётся уравновесить «ящики» посещений этих конформаций (в них попадает очень мало) (рис. 3). Так что на рис. 3 участок графика 2 для т > 128 N = 100) подлежит специальному рассмотрению. И даже для для N = 30 и 50 последняя точка диаграммы 3 сильно выбивается из равномерного посещения, т. е. сложно моделировать компактные конформации. Компактные конформации вносят существенный вклад в конфигурационную энергию при низких температурах и е < 0. В остальных случаях «хвост» распределения не будет оказывать значительного влияния на
1,E+01 1,E-01 1,E-03 1,E-05 1,E-07 1,E-09 1,E-11 1,E-13 1,E-15 1,E-17 1,E-19 cf 1'E-21 1,E-23 1,E-25 1,E-27 1,E-29 1,E-31 1,E-33 1,E-35 1,E-37 1,E-39 1,E-41 1,E-43 1,E-45
Рис. 4- Распределение по числу контактов для звёзд (чёрные точки) с длиной лучей 5, 12, 20 сегментов и цепей с таким же (как и в звезде) числом сегментов N = 30 (круги), 72 (квадраты), 120 (треугольники):
вертикальные пунктирные линии соответствуют максимальному числу контактов для цепей, полученному в процессе моделирования
результат для конфигурационной энергии. На рис. 2 для сравнения также представлены данные, полученные в работе [6] с помощью фантомного блуждания. Видно, что использование фантомных цепей не позволило получить информацию о распределении по числу контактов в маловероятных случаях. Особенно сильно результаты различаются для длинных цепей (N = 100).
На рис. 4 показано распределение по числу контактов для звёзд с длиной лучей N = 5,12, 30 и для цепей с таким же числом сегментов, как и у звёзд, т. е. N = = 30, 72,120. Зависимости для звёзд и цепей различаются, но достаточно близки друг к другу. Для звёзд, как и для цепей, распределение имеет один максимум.
Литература
1- Romiszowski P., Sikorski A. Temperature dependence of properties of star-branched polymers: a computer simulation study //J- Chem. Phys. 1998- Vol. 109- N 7- P. 2912-2920.
2. Cameron D. J. A., Shaver M. P. Aliphatic polyester polymer stars: synthesis, properties and applications in biomedicine and nanotechnology // Chem. Soc. Rev. 2011. Vol. 40. P. 1761-1776.
3. Бирштейн Т. М., Меркурьева А. А., Лирмэйкерс Ф. А. М., Рудь О. В. Конформации полимерных и полиэлектролитных звёзд // Высокомол. соединения. (А). 2008. Т. 50. № 8. C. 1-18.
L_l_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_и_I_I_I_I_I_I_I_I_I_I_I_I_I_1_L
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170
4. AryalS., Prabaharana M., PillaS., GongaS. Biodegradable and biocompatible multi-arm star amphiphilic block copolymer as a carrier for hydrophobic drug delivery // Int. J. Biol. Macro-molecules. 2009. Vol. 44. P. 346-352.
5. Georgiou T. K, VamvakakiM., Phylactou L. A., Patrickios C. S. Synthesis, characterization and evaluation as transfection reagents of double-hydrophilic star copolymers: effect of star architecture // Biomacromolecules. 2005. Vol. 6. N 6. P. 2290-2997.
6. Vorontsov-Velyaminov P. N., VolkovN. A., Yurchenko A. A. Entropic sampling of simple polymer models within Wang—Landau algorithm //J. Phys. (A). 2004. Vol. 37. P. 1573-1588.
7. Volkov N. A., Yurchenko A. A., Lyubartsev A. P., Vorontsov-Velyaminov P. N. Entropic sampling of free and ring polymer chains // Macromol. Theory and Simul. 2005. Vol. 14. P. 491-504.
8. BergB. A., Neuhaus T. Multicanonical ensemble: a new approach to simulate first-order phase transitions // Phys. Rev. Lett. 1992. Vol. 68. P. 9-12.
9. Lee J. New Monte Carlo algorithm: entropic sampling // Phys. Rev. Lett. 1993. Vol. 71. P. 211-214.
10. WangF., Landau D. P. Efficient, multiple-range random walk algorithm to calculate the density of states // Phys. Rev. Lett. 2001. Vol. 86. P. 2050-2053.
11. Parsons D. F., Williams D. R. M. An off-lattice Wang—Landau study of the coil-globule and melting transitions of a flexible homopolymer // J. Chem. Phys. 2006. Vol. 124. P. 221103.
12. Parsons D. F., Williams D. R. M. Globule transitions of a single homopolymer: A Wang—Landau Monte Carlo study // Phys. Rev. (E). 2006. Vol. 74. P. 041804.
13. SeatonD. T., Wust T., Landau D. P. A Wang—Landau study of the phase transitions in a flexible homopolymer // Computer Phys. Communications. 2009. Vol. 180. P. 587-589.
14. SeatonD. T., Wust T., Landau D. P. Collapse transitions in a flexible homopolymer chain: application of the Wang—Landau algorithm // Phys. Rev. (E). 2010. Vol. 81. P. 011802.
15. Melas V. B. Optimal Simulation Design by Branching Technique // Model-Oriented Data Analysis / ed. by W.G.Mueller, H. P. Wynn, A. A. Zhygliavsky. Heidelberg, 1993. P. 113-127.
16. ClisbyN. Efficient Implementation of the Pivot Algorithm for Self-avoiding Walks //J. Stat. Phys. 2010. Vol. 140. P. 349-392.
17. Douglas J., Ishinabe T. Self-avoiding-walk contacts and random-walk self-intersections in variable dimensionality // Phys. Rev. (E). 1997. Vol. 51. P. 1791-1817.
Статья поступила в редакцию 28 июня 2011 г.
