Моделирование двумерного полимерного расплава методом молекулярной динамики Текст научной статьи по специальности «Физика»

ВЫСОКОМОЛЕКУЛЯРНЫЕ СОЕДИНЕНИЯ, Серия А, 2002, том 44, № 7, с. 1228-1239

ТЕОРИЯ

УДК 541.64:539.199

МОДЕЛИРОВАНИЕ ДВУМЕРНОГО ПОЛИМЕРНОГО РАСПЛАВА МЕТОДОМ МОЛЕКУЛЯРНОЙ ДИНАМИКИ1

© 2002 г. Н. К. Балабаев*, А. А. Даринский**, И. М. Неелов**, Н. В. Лукашева**, I. Emri***

* Институт математических проблем биологии Российской академии наук 142292 Пущино Московской обл., пр. Науки, 4

**Институт высокомолекулярных соединений Российской академии наук 199004 Санкт-Петербург, Большой пр., 31

***Center of Experimental Mechanics, University of Ljubljana Costa na Brdo 49, SI-1125 Ljubljana, Slovenia

Поступила в редакцию 06.08.2001 г.

Принята в печать 28.01.2002 г.

Методом молекулярной динамики проведено моделирование расплава двумерных полимерных цепей. В качестве модели используется цепь, состоящая из частиц, соединенных упругими пружинами-связями. Взаимодействие между частицами описывается модифицированным потенциалом Леннар-да-Джонса, учитывающим только отталкивание частиц. Моделировали расплавы, состоящие из цепей с N = 20-60 звеньев. Цепи в расплаве коллапсированы и сегрегированы. Среднеквадратичные размеры цепей пропорциональны N, а функция распределения цепей по расстояниям между концами близка к гауссовой. В области смещений, меньших радиуса инерции, временные зависимости среднеквадратичных смещений мономеров характеризуются наклоном, несколько превышающим величину 1/2, характерную для модели Рауза. Время релаксации к-й нормальной моды цепей пропорционально (N + 1 )2/к2, как и для модели Рауза. Однако диффузия цепи как целого происходит несколько быстрее, чем предсказывает модель Рауза. Рассчитаны также характеристики локальной и глобальной ориентационной подвижности цепей.

ВВЕДЕНИЕ

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

1 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 99-03-33314), фонда ЮТАБ (грант 99-1114), Федеральной целевой программы "Интеграция" (проект 326.38) и Программы ЕЭР "5иРЕ1иЧЕТ'.

E-mail: adar@imc.macro.ru (Даринский Анатолий Анатольевич).

пересечений (цепь в хорошем растворителе) среднеквадратичные размеры цепи растут с увеличением числа звеньев как (Я2) ~ /У3/2, в отличие от закона (/?2) ~ Л/6/5 для трехмерного случая. Выводы теории были подтверждены и компьютерным моделированием. Так, в работе [4] методом Монте-Карло проводили моделирование структуры и динамики разбавленного раствора полимерных цепей с учетом объемных взаимодействий в узкой поре. Было показано, что с уменьшением ширины щели зависимость размеров цепи от ее длины N изменяется от (Я2) ~ М6/5 для трехмерных цепей до (К2) ~ Ы312 для модели двумерных цепей с запретом самопересечений. Максимальное время релаксации для такой цепи возрастает с увеличением N согласно скелинговой зависимости для протекаемой цепи ттах ~ М1'^*1 = А/2-5. Поскольку для цепи на поверхности гидродинамиче-

ские взаимодействия практически полностью экранируются, использование протекаемой модели вполне оправдано. Недавно было проведено экспериментальное исследование разбавленных и полуразбавленных растворов ДНК на поверхности жидких липидных слоев [5]. Такие системы представляют собой пример двумерного полимерного раствора. Полученные результаты хорошо согласуются с предсказаниями теории и результатами моделирования.

Различия между трехмерной и двумерной цепями проявляется и в поведении цепей в 8-раство-рителе (когда второй вириальный коэффициент равен нулю). Для двумерной цепи теория и компьютерное моделирование дают зависимость {К2) ~ М8/7 [6-9], в отличие от трехмерного случая </?2) ~ N. При переходе к концентрированному раствору и расплаву роль объемных взаимодействий ослабевает. Уже довольно давно установлено, что в трехмерном полимерном расплаве объемные взаимодействия экранируются [10], и конформация цепи в расплаве соответствует гауссовой конфор-мации цепи в в-растворителе. В двумерном расплаве из-за топологических ограничений цепи остаются сегрегированными и не проникают друг в друга. Поэтому звенья выделенной цепи взаимодействуют в основном со звеньями той же цепи. В работе [11] было проведено моделирование структуры и динамики двумерных концентрированных растворов полимеров вплоть до расплава на решеточной модели цепи. Использовали вариант метода Монте-Карло (метод флуктуирующих связей). Моделировали системы, состоящие из цепей с N = 20-100. Было показано, что, несмотря на сегрегирован-ность цепей в расплаве, зависимость их размеров от N также соответствует идеальной гауссовой модели </?2) ~ N. Таким образом, конформация цепи в расплаве скорее близка к конформации двумерной глобулы, где также (Я2) ~ N.

Сходство с идеальной гауссовой цепью наблюдали и для динамических свойств цепи. Было показано, что на временной зависимости среднеквадратичного смещения (Аг*(0) мономера, характеризующей его поступательную подвижность, наблюдается область (Аг2) ~ г1/2, соответствующая модели субцепей Рауза. Нормальные моды для модели субцепей Рауза оказались нормальными модами и для двумерной цепи в расплаве. Более того, зависимость времен релаксации нормальных мод хк от номера моды к и от N оказывается аналогичной соответствующей зависимости для крупномасштабных мод модели Рауза хк ~ (к№)2. Следует отметить, однако, что динамика цепи, реализуемая с помощью метода флуктуирующих связей, обладает

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

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

МОДЕЛЬ И МЕТОД

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

ад = 4г0[^)^)6] + е0 при г <2%

U(r) = 0 при г>2тс

Здесь г - расстояние между частицами. Такой потенциал в литературе иногда называют потенциалом Уикса-Чандлера-Андерсена [14]. Для задания связей между соседними по цепи частицами вводится дополнительный потенциал, соответствующий френкелевской пружине Ub = 0.5Kb(r-a)2. Расчеты проводили для значения силовой константы Кь = 50.

Для моделирования системы при заданном внешнем давлении Pref (т.е. при переменной площади) и температуре Tref использовали форму

уравнений движения, предложенной Берендсе-ном и др. [15] для NPN-ансамбля

^jf = v[.„ + pPii>(Ptl„-P„/K,,„ (1) - ,)„,,„. (2)

dt

(a)

; < г > J. з

» iw i r*^ v ""w • W» w* * i, / / t /

L/ v } \ j

i f S f 1

t r i t i 1

• » • » - » -.y». ». «■ - %

(6)

(B)

Рис. 1. Конфигурации системы на последовательных стадиях приготовления. Число шагов 1200 (а), 5600 (б), 100000 (в). Рисунки (а)-(в) отвечают разным размерам ячейки.

где г = 1,и, а = Х или У; Г - мгновенная температура, определяемая из отношения полной кинетической энергии всей системы ЕШп к числу степе-

1 £ ней свободы - къТ = ^^'

Компоненты Раа мгновенного тензора давления рассчитывали с помощью теоремы вириала с учетом периодических граничных условий

(3)

Здесь 5 - площадь расчетной ячейки. Суммирование в выражении (3) проводится по всем парам частиц. Для каждой частицы в качестве пары берется ближайший имидж другой частицы. При моделировании системы при постоянном объеме (ИУТ-ансамбль) величины и в уравнениях (1), (2) полагали равным нулю.

Численное интегрирование уравнений движения проводили с помощью алгоритма Верле [16], который дает точность порядка 0(к3) для скоростей частиц на каждом шаге. При выводе формул, используемых при расчетах, предлагали, что температура и давление в системе меняются медленнее, чем динамические переменные (координаты и скорости частиц). Поэтому вместо мгновенных температуры и давления, определяемых уравнениями (1) и (2), использовали значения этих величин, усредненные по предшествующим двадцати шагам, для расчета сил парных взаимодействий мы применяли метод связанных листов [16, с. 149]. Необходимые для реализации этого алгоритма число операций и память растут линейно с увеличением числа частиц в системе. В качестве единиц длины и энергии были выбраны, как обычно, параметры потенциала Леннарда-Джонса а и е^ соответственно, а в качестве единицы времени - характерное

время потенциала Леннарда-Джонса т0

- с

V£o'

Температура Т задается в единицах Т0 = е$/кв (кв -константа Больцмана). Все расчеты проводили при Т= 0.5.

ПРИГОТОВЛЕНИЕ СИСТЕМЫ

Систему готовили по методике, сходной с использованной в работах [12,13]. Этапы приготовления системы проиллюстрированы на рис. 1. Начальная конфигурация системы соответствует регулярной сетке, узлы которой располагаются на квадратной решетке и цепи между узлами вытянуты вдоль ее двух основных направлений. На первом этапе разрываются связи между частицами

в узлах и цепями. В результате мы имеем пс цепей и ns отдельных частиц. Наличие в системе частиц, не включенных в цепи, позволяет изучать диффузию этих частиц в полимерной матрице. Полное число частиц в системе п = ncN + ns. Число звеньев N в цепях варьировалась от 20 до 60. Расчеты выполняли для систем, состоящих из пс = 50 цепей и пс = 25 отдельных частиц. На следующем этапе система сжималась при постоянном внешнем давлении до заданной (численной) плотности р = n/V. При выборе плотности изучаемой системы стремились к тому, чтобы наша система была близка по плотности к модели расплава, рассмотренной в работе [11], где использовали модель цепи на квадратной решетке. Самая высокая степень замещения узлов решетки, при которой проводили моделирование, составляла 0.8 от максимальной. Если частицу в нашей модели рассматривать как

сферу с объемом ^ су3, то максимально плотной

упаковке таких сфер на плоскости соответствует доля занятого объема (поверхности), равная 0.78. Это отвечает р = 1. Моделирование проводили для систем с р = 0.78.

СТАТИЧЕСКИЕ СВОЙСТВА

На рис. 1 представлена одна из конфигураций системы при р = 0.78, снятая с экрана дисплея. Видно, что цепи сегрегированы и кол л апсированы, что и следовало ожидать в силу невозможности взаимопересечений на плоскости. Для характеризации статических свойств цепей в системе были рассчитаны среднеквадратичные проекции расстояния между

концами цепи а (Л/)) = ((г, а - гк „)2), где а = X, У и компоненты среднеквадратичного радиуса инер-

N

2 1 _ ции (Я^ а (Л0> = « - ГЛГ, а)2}- где гк а, г^ а -

* = 1

координаты к-то мономера и центра масс соответственно. Результаты расчетов приведены в табл. 1.

Наклон рассчитанных зависимостей ^ (/га) и

^ (/?й а) от 1ЕЛГ (рис. 2) близок к единице в согласии с предсказаниями теории и результатами моделирования. Известно, что для гауссовых цепей

отношение (/г2)/( ) = 6. Рисунок 2 показывает, что это отношение близко к шести и для рассмотренной модели двумерной цепи в расплаве. Заметим, что в работе [11] величина отношения была несколько меньше. Более подробную инфор-

Таблица 1. Характеристика размеров цепей

N (h2) Сh2)/(R¡>

20 43 7.69 5.59

40 89.2 16.7 5.34

60 12.9 23.08 5.59

мацию о конформационных свойствах цепи дает функция распределения. На рис. 3 приведены рассчитанные функции распределения Р(1гх) для х-проекций вектора расстояния между концами. Видно, что для всех длин цепей наблюдается

lg«*,2» lg«*2» 2.0

6 -

(б)

1

я-_ --■

i ----- i 2

20

40

60

N

Рис. 2. Среднеквадратичные размеры цепей как функции их длины, а - зависимости ^-проекций среднего квадрата расстояния между концами

цепей ((Лх)) (1) и среднего квадрата радиуса

инерции цепи ((Ях)) (2) от N в логарифмическом

2 2

масштабе; б - отношение (Ах )/(/?х) как функция

N для моделируемой системы (Ли гауссового клубка (2).

Р(Ь)

М{/!2))1/2

Рис. 3. Функции распределения х-проекций расстояний между концами цепей с N = 20 (1), 40 (2) и 60 (3), представленные в зависимости от

к'/({кгх ))|/2. Штриховыми линиями показаны гауссовы функции с рассчитанными значениями (и] > 21.5 = 20), 44.6 (ЛГ = 40) и 74.3 (ЛГ = 60).

((х({) - л'(0))2)

5 10

?х 10"3

Рис. 4. Среднеквадратичные смещения свободных, не включенных в цепь частиц для N=20 (1), 40 (2) и 60 (3). Штриховыми линиями показаны линейные аппроксимации. Для N = 20 и 40 максимальное время 2000.

очень хорошее согласие с гауссовой функций

1/2 ( ,24

ехр

К

, рассчитанной при

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

Таблица 2. Характеристики поступательной подвижности отдельных частиц, мономеров и цепей

N ост А)

20 0.013 0.0035 0.07 0.043

40 0.01 0.00175 0.07 0.049

60 0.01 0.00175 0.105 0.035

ДИНАМИЧЕСКИЕ СВОЙСТВА

Локальная поступательная подвижность

Для характеристики локальной поступательной подвижности проанализировали среднеквадратичные смещения вдоль координатных осей (Дл:2(0) как для отдельных частиц, не включенных в цепь, так и для мономеров в цепи. Для остальных частиц, не включенных в цепи, зависимости (Ах2(0) (рис. 4) при малых г < 50-100 обнаруживают некоторую кривизну. Это указывает на то, что в данной области времен инерционные эффекты играют заметную роль. При больших г зависимости (Дх2(0) выходят на линейный режим, что указывает на чисто диффузионный характер движения. Значения коэффициентов диффузии /)5 для отдельных частиц в трех рассмотренных системах, определенные из наклонов этих зависимостей, представлены в табл. 2. Диффузия ускоряется с уменьшением N. Для наших систем с постоянной поверхностной плотностью это ускорение можно связать с увеличением удельного свободного объема (поверхности), обусловленно-

го ростом числа свободных концов, вблизи которых свободный объем выше. Однако временные зависимости (Ajtfy)) для мономеров, расположенных как в середине, так и на концах цепи, криволинейны и в этой области времен, что отражает эффект их включения в цель. Поступательная подвижность концевых мономеров, естественно, выше, чем подвижность центральных мономеров, но начальные наклоны для обеих зависимостей близки (на рис. 5 показаны эти зависимости для N = 40). Для обеих зависимостей (как для центральных, так и для концевых моно-

меров) в координатах lg(Ajc (/)) , lg t (рис. 6) наблюдаются две временные области. В области смещений, меньших радиуса инерции цепи, обе зависимости характеризуется наклоном, несколько превышающим величину 1/2, характерную для модели гауссовых субцепей. Отметим, что в работе [11] в этой области смещений наблюдали наклон, равный 1/2. Возможны две причины такого расхождения: либо изучаемые нами цепи недостаточно длинные (в работе [11] рассматривали цепи с N до 100), либо проявляется эффект корреляций между движениями соседних мономеров, который, как уже говорилось, не учитывается в рамках алгоритма Монте-Карло. При больших

временах и смещениях зависимости lg(Ax (0) от lgí характеризуются наклоном, близким к единице. В этой области проявляется смещение цепей как целого.

Локальная ориентационная подвижность

Для характеристики локальной ориентацион-ной подвижности звеньев цепей были рассчитаны ориентационные автокорреляционные функции

Pl(t)nP2(t)

/>,(') = <cos9(í)>

P2{t) = 2(cos20(í) - 1/2) Здесь 0(0 - угол поворота звена за время t.

Соответствующие временные зависимости для концевых и центральных звеньев приведены на рис. 7. Заметим, что для отдельной жесткой гантели, состоящей из двух центров вязкого сопротивления с коэффициентами трения ¿¡о, соединенных жесткой связью длины /, совершающей броуновское движение в вязкой среде на плоскости, Pj(í) и P2(t) описываются простыми экспоненциальными функциями Ру 2 - exp(-í/t12), где

Рис. 5. Среднеквадратичные смещения свободных, не включенных в цепи частиц (У), концевых (2) и центральных мономеров цепи (5) для М = 40.

1пг

Рис. 6. Зависимости логарифмов среднеквадратичных смещений 1п(((*(/) - л(0))2)) от 1п(/) для средних мономеров цепи. N=20 (Г), 30 (2), 40 (3) и 60 (4). Штриховыми линиями показаны линии с наклонами 1/2, характерным для гауссовых субцепей, и 1, характерным для смещения цепей как целого.

Для звеньев, включенных в цепь, экспоненциальный характер релаксации нарушается. Для рассмотренной модели зависимости Р^г) и Р20) неплохо представляются так называемыми дроб-

Р,

ными экспонентами ехр(-г/т,) , где / = 1,2. Времена релаксации Т! и т2 и параметры 3, и (32 для соответствующих корреляционных функций Ру{1) и Р2(0 (табл. 3) практически не зависят от длины цепи, что отражает локальный характер ориента-ционного движения. Времена релаксации т1 для

Р1« 1.0-

0.6

0.2

(а)

□ А 7 Б ж В

Р-ф 1.01-

0.6

1 I | 0.2

' ■ * • я

(б)

* 1

п »1

2 "1

200

400

*

Я

*9*Ш*ЧШШШШШШЮ

20 40 60 г

Рис. 7. Временные зависимости ориентационных корреляционных функций /^(г) (а) и Р2(г) (б) для центральных. (1) и концевых (2) звеньев. N = 20 (Л), 40 (Б) и 60 (В).

концевых звеньев заметно (в 3-3.5 раза) меньше, чем для центральных звеньев. Что же касается времен т2, то их значения для центральных и концевых звеньев близки. Значения параметров р[ и р2 для концевых звеньев выше, чем для централь-

Таблица 3. Времена релаксации и параметры Р для корреляционных функций Р[(г) и Р2(0

N Ч Р1 02 1,/Хг

Конец

20 22 3.60 0.50 0.58 6.2

40 21 3.06 0.48 0.52 7.0

60 23 2.86 0.51 0.50 8.0

Середина

20 72 3.05 0.32 0.45 24.0

40 72 2.71 0.30 0.42 26.8

60 68 2.09 0.31 0.39 34.0

ных, что отражает более широкий спектр для центральных мономеров. В табл. 3 представлены также рассчитанные отношения т^. Для концевых звеньев эти отношения лишь немного выше теоретического значения, равного 4, для отдельной жесткой гантели. Для центральных звеньев, однако, эта величина почти на порядок выше. Следует иметь в виду, что соотношение (4) между временами т, и т2 получено для модели гантели в вязкой среде, когда инерционными эффектами можно пренебречь. Такое условие наверняка не выполняется в области наблюдаемых времен (т2 ~ ~ 2-3) (см. выше обсуждение поступательной диффузии свободных частиц).

Нормальные моды

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

1п(сог?) Ок-

Рис. 8. Корреляционные функции нормальных мод с к = 1-7 для цепей с N=60. Штриховыми линиями показаны линейные аппроксимации.

-1

-2

-3

-4

' *

' ж

/ V

/ ж

v

□ I

V 2

ж 5

-1.8

-1.2

-0.6 0

Рис. 9. Логарифмы времен релаксации нормальных мод как функции +1)) для первых 18 нормальных мод для цепей с N = 20 (1), 40 (2) и 60 0). Наклон штриховой линии равен двум.

ствуют линейные комбинации проекций звеньев цепи ир [17]

Як

р= 1

ккр N +1'

(5)

С(Як> Ч = -—2-

\Як)

т.. -

и

равен 12), где С,т - мономерный коэффициент трения, I2 - среднеквадратичная длина субцепи.

гт - к

Для малых значении м 1 получаем

которые являются нормальными модами с номерами к. Это означает, что корреляционные функции

Ъ =

N+1

(н+1)2Ы2

Ч ~ 2

2 пкиТк

(7)

(6)

релаксируют по простому экспоненциальному закону с временами релаксации хк

С(Яь 0 = ехр(-//т*)

Для двумерного варианта модели Рауза времена релаксации даются выражением

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

Рассчитанные временные зависимости логарифмов корреляционных функций (6) нашей системы показаны на рис. 8. Эти зависимости близки к линейным, т.е. можно утверждать, что линейные комбинации (5) в хорошем приближении являются нормальными модами и для двумерной цепи в расплаве. Из наклонов временных зависимостей 1п г) определены времена релаксации хк. На

рис. 9 построены зависимости и ^д^ первых десяти мод для N = 20 и первых 18 мод для N = 40 и 60. Все зависимости ложатся практически на одну прямую с наклоном, близким к двум, т.е. значению, соответствующему модели Рауза. Этот результат хорошо согласуется с данными моделирования методом Монте-Карло [11]. Таким образом, дисперсионная зависимость нор-

<(хст(0-хст(0))2> 18

Рис. 10. Зависимости среднеквадратичных смещений центров масс (((хст(1) - хст(0))2)) отдельных цепей от г. N=20 (7), 40 (2) и 60 (3). Штриховыми линиями показаны линейные аппроксимации.

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

Согласно выражению (7), из наклона зависи-

к ^2

МОСТИ Ть от

N

можно оценить коэффици-

ент диффузии мономера й„

к Т

= -у—. Величина I2

'э т

оценивалась из соотношения (И2) = Ш2 при использовании рассчитанных величин (Л2) (табл. 1). Значения От приведены в табл. 2. Интересно, что От оказывается существенно больше, чем коэффициент диффузии для отдельных частиц, не включенных в цепь. На первый взгляд это кажется парадоксальным. Отметим, однако, что начальные наклоны для среднеквадратичных смещений отдельных частиц и частиц, включенных в цепь, близки (рис. 5). Это значит, что при очень малых смещениях частицы не чувствуют связи в цепи. При больших смещениях начинают играть роль жесткие связи между частицами цепи. Величина йт рассчитана на одну частицу цепи. Боль-

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

Диффузия цепей как целого

Для характеристики диффузии цепей как целого

были рассчитаны {Ахсет (г)) для центров масс отдельных цепей (рис. 10). Как и для отдельных частиц, не включенных в цепь, эти зависимости при больших г линейны. Коэффициент самодиффузии цепей £>ст определяли из наклона зависимостей при

2 2

больших смещениях (Ахсет (г)) > ). Для наиболее длинных цепей с N = 60 удается получить сме-щения лишь до 0.7 ). Поэтому Ост для длинных цепей недостаточно достоверны. Рассчитанные значения Ост приведены в табл. 2.

Известно, что для трехмерных расплавов гибких цепей в дорептационной области молекулярных масс Ост ~ М-1 [18]. Для нашей двумерной модели значения Ост для N = 20 и 40 соответствуют зависимости Для N = 60 получены, по-видимому, завышенные значения Д.т. Отметим, что и в работе [211] четкая зависимость йст ~ Ы~1 не наблюдалась. Авторы также связывали этот результат с недостаточной статистикой при больших N. Как и авторы [11], мы полагаем, однако, что такая зависимость будет выполняться, если проводить моделирование достаточно долго. Из соотношения Д, = Г)ст№ можно оценить эффективный коэффициент диффузии одного мономера £>0 (табл. 2). Интересно, что определенный таким образом эффективный коэффициент диффузии несколько превышает мономерный коэффициент диффузии йст, полученный из данных для времен релаксации нормальных мод. Таким образом, поступательная диффузия цепи как целого происходит быстрее, чем подсказывает модель субцепей.

Аналогичное поведение было обнаружено в работе [11], где величины эффективного и мономерного коэффициентов диффузии различались более, чем на порядок. Там же была предложена интерпретация такого различия: в силу невозможности взаимопроникновения цепей в двумерной модели диффузии цепи как целого происходит путем флуктуаций формы цепи. Этот механизм для двумерной диффузии оказывается более

1п(согА) О

cor(Rl), cor(/z2)

1.0b-

0.8

0.6

0.4

0 „

0 0 Я n

® 5 8 □ g

ё 8 □ * D

z □ д д 5 □ n

s п _ д * 9 □

ö

д °

Д □ „ д □

' 4 1

2S°g

□ а а 5

Рис. 11. Временная зависимость логарифма корреляционной функции для вектора расстояния между концами цепи для N = 20(1, Г), 40 (2,2') и 60 (3, 3'). Сплошными линиями представлены результаты моделирования, штриховыми - теоретические кривые, рассчитанные по временам релаксации нормальных мод.

I

2 4

txlOr2

Рис. 12. Временные зависимости крупномасштаб-

ных корреляционных функций: для квадрата рас-

стояния между концами (а) и квадрата радиуса

инерции цепей (б). N=20 (7), 40 (2) и 60 (3).

быстрым, чем механизм внутрицепных перестроек, определяющих релаксацию раузовских мод. Количественные различия зависят, по-видимому, от используемой модели. В работе [11] рассматривалась модель цепи на квадратной решетке, и движение в цепи осуществлялось вследствие перестроек кинетических единиц заданной структуры. Мы анализируем континуальную модель, в которой механизм движений не постулируется заранее.

Глобальная ориентационная подвижность

Для характеристики ориентационной подвижности цепей как целого рассчитали корреляционные функции для вектора

_ (h(O)h(t)) Chit) ~

и квадрата вектора расстояния между концами

(h2(0)h\t))-{h2)2 (h4)-{h2)2 '

а также корреляционную функцию для квадрата радиуса инерции

С „ч _ (Д»(0Я»(°)> - <*,У

r( '--1 Тг-

<*t>-<*v

Временные зависимости логарифмов данных корреляционных функций показаны на рис. 11 и 12. Эти функции хорошо представляются суммой двух экспонент с малым т' и большим т" временами релаксации

C(t) = Aexp(-th') + (1 - А)ехр(-г/т")

В табл. 4 приведены значения А, т' и т". Времена т' и т" различаются более, чем на порядок. Короткие времена т' меняются мало при изменении N. Главный вклад в релаксацию ориентационных корреляционных функций вносят большие времена т". Зависимости Igt" от IgN (рис. 13) показывают, что ~ N1-25, т"2 ~ N135, х'\ ~ /V1-45. От-

h R

метим, что в работе [11] была получена более сильная зависимость tr2 ~ N2. Времена для квадратов расстояния между концами т"> и радиуса

Таблица 4. Времена релаксации для крупномасштабных корреляционных функций

N Ah

20 0.09 63 1240 0.2 47 440 0.2 46 480 2.82

40 0.05 75 2980 0.1 54 1130 0.1 59 1280 2.63

60 0.03 69 4690 0.07 67 1890 0.09 97 2110 2.50

инерции близки. Отношение времен хЦ /х'\ =

= 2.6-2.8 близко к теоретическому значению, равному 3, для модели Рауза. Теория для гауссовых субцепей [17] предсказывает связь между корреляционной функцией для вектора расстояния между концами цепи и временами релаксации нормальных мод

N

С"(0 = N(N+1) X ^ехр(-г/т,), (8)

к = 1

где суммирование проводится лишь по нечетным значениям номера моды к. На рис. 11 результаты моделирования сопоставлены с теоретическими

Ч^Л&хЪЛш^)

№

Рис. 13. Зависимости времен релаксации размеров цепей: расстояния между концами цепи хк (/), квадрата расстояния между концами цепи т 2 (2) и

квадрата радиуса инерции т 2 0) от Ы, представ-

Л

ленные в двойном логарифмическом масштабе. Наклон штриховой линии равен двум.

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

СПИСОК ЛИТЕРАТУРЫ

1. Edwards S.F. // Proc. Phys. Soc. London. 1965. V. 86. P. 613.

2. Pereira GJ. // Phys. A. 1997. V. 30. P. 467.

3. Cardy J.L., Hamber H.W. // Phys. Rev. Lett. 1980. V. 45. P. 499.

4. Milchev A., Binder К. // J. Phys. France. 1996. V. 6. P. 21.

5. Maier В., Radler J.O. // Macromolecules. 2000. V. 33. P. 7185.

6. Jan N„ Coniglio A., Majid /., Stanley H.E. Growth and Form-fractal and Non-fractal Patterns in Physics. Dordrecht: Springer, 1986. P. 263.

7. Weintrib A., Trugman SA. // Phys. Rev. B. 1985. V. 31. № 5. P. 2993.

8. Бириипейн T.M., Булдырев C.B., Ельяшевич AM. // Высокомолек. соед. А. 1986. Т. 28. № 3. С. 634.

9. Birshtein T.M.,Buldyrev S.V., Eliyashevich AM. // Polymer. 1985. V. 26. № 11. P. 1814.

10. Flory PJ. Principles of Polymer Chemistry. New York: Cornel Univ. Press, 1953.

11. Carmesin /., Kremer K. // J. Phys. France. 1990. V. 51. P. 915.

12. Даринский A.A., Неелов И.М. Балабаев Н.К., Sundholm F. // Высокомолек. соед. А. 1998. Т. 40. №7. С. 1110.

13. Неелов U.M., Даринский A.A., Балабаев Н.К., Sundholm F. // Высокомолек. соед. А. 1998. Т. 40. № 12. С. 1963.

14. Weeks J.D., Chandler D., Andersen H.C. // J. Chem. Phys. 1971. V. 54. P. 5237.

15. Berendsen HJ.C., PostmaJ.P.M., van Gunsteren W.F., DiNolaA.,HaakJ.R. //J. Chem. Phys. 1984. V. 81. № 8. P. 3684.

16. Allen M.P., Tidsley DJ. Computer Simulations of Liquids. Oxford: Clarendon, 1987.

17. Готлиб Ю.Я., Даринский A.A., Светлов Ю.Е. Физическая кинетика макромолекул. Л.: Химия, 1986.

18. Дои М„ Эдварде С. Динамическая теория полимеров. М.: Мир, 1998.

Molecular Dynamics Simulation of a Two-Dimensional Polymer Melt

N. K. Balabaev*, A. A. Darinskii**, I. M. Neelov**, N. V. Lukasheva**, and I. Emri***

*Institute of Mathematical Problems of Biology, Russian Academy of Sciences, pr. Nauki 4, Pushchino, Moscow oblast, 142292 Russia

**Institute of Macromolecular Compounds, Russian Academy of Sciences, Bol'shoipr. 31, St. Petersburg, 199004 Russia

***Center of Experimental Mechanics, University of Ljubljana, Costa na Brdo 49, SI-1125 Ljubljana, Slovenia

Abstract—A melt of two-dimensional polymer chains was simulated by molecular dynamics. A chain consisting of particles linked by elastic strings (bonds) was used as a model. The interaction of the particles was described by the modified Lennard-Jones potential that accounts only for repulsion of particles. Melts of chains consisting of N = 20-60 units were modeled. The chains in a melt are collapsed and segregated. The average square sizes of chains are proportional to N and the distribution function of the end-to-end distances of chains is close to Gaussian. In the range of displacements smaller than the gyration radius, the time dependences of the average square displacements of monomers are characterized by a slope slightly larger than 1/2 typical of the Rouse model. As in the case of the Rouse model, the relaxation time of the &th normal mode of chains is proportional to (N + l )2/k2. However, the diffusion of a chain as a whole is slightly faster than predicted by the Rouse model. The characteristics of the local and global orientational mobilities of chains were also calculated.