Использование метода мультимасштабного моделирования для прогнозирования свойств: сверхсшитый полистирол*

Abstract

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

-

ной макромолекуле в вакууме (полистирол [ 21 ] и сополимер полистирол-бензоциклобутен [1920]).

Данная работа посвящена мультимасштабному компьютерному моделированию синтеза сверхсшитого полистирола из линейных молекул полистирола, растворенных в дихлорэтане с использованием монохлордиметилового эфира в качестве сшивателя. Сначала проводилось атомистическое моделирование раствора полистирола в дихлорэтане. Его результаты использовались для создания огрубленной мезоскопической модели, которая применялась при моделировании реакции сшивки. Результаты моделирования в огрубленной модели превращались с помощью обратного картирования в атомистическую структуру сверхсшитой сетки и затем исследовались. 2. Модель 2.1 Атомистическое моделирование

Атомистическое моделирование проводилось с помощью пакета программ GROMACS с использованием набора параметров 53A6 силового поля GROMOS [ 22 ]. При подготовке потенциалов связей, валентных и торсионных углов использовался ATB [ 23 ]. Частичные заряды мономерных звеньев полистирола брались из результатов ab initio вычислений для 1,3,5-трифенил гексана, олигомера полистирола, состоящего из трех звеньев. Конфигурация для ab initio вычислений определялась методом сопряженных градиентов[ 24 ] для двух различных полуэмпирических потенциалов AM1[ 25 ] и PM3[26,27]. В ab initio вычислениях использовался базовый набор функций 3-21G[ 28 ] для полистирола и 6-31G* для 1,2-дихлорэтана.

Ячейка моделирования содержала 16 цепей полистирола, каждая из которых состоит из 20 мономерных звеньев, и молекул дихлорэтана, число молекул которого варьировалось от 10 до 1900, так что концентрация полимера изменялась от 103 до 17%w/v. Моделирование проводилось в пакете для молекулярного моделирования GROMACS[ 29-33 ] с шагом интегрирования 0.001 пс. Радиус обрезки Ван-дер-Ваальсовых взаимодействий был равен 1.6 нм. Электростатические взаимодействия обсчитывались алгоритмом PME[ 34 ] с радиусом обрезки 1.6 нм. Моделирование проводилось в каноническом ансамбле. Постоянная температура поддерживалась методом перенормировки скоростей [ 35 ] с временной константой 0.1 пс. Постоянное давление поддерживалось методом Берендсена [ 36 ] с временной константой 0.5 пс.

Начальные конфигурации приготовлялись следующим образом: молекулы полистирола сначала уравновешивались в вакууме, затем помещались в большую ячейку, затем в нее добавлялись молекулы 1,2-дихлорэтана. Моделирование начиналось при начальной плотности около 700 кг/м3, при температуре 350 K и нормальном давлении. На 6 нс включается высокое давление в 100 атм, а затем моделирование проводится в течение 7.5 нс при нормальном давлении. Затем для ускорения релаксации и перемешивания проводится процедура отжига[ 37,38 ]. При этом чередуется моделирование при повышенной (T=500 K) и заданной (T=350 K) температурах, в каждом случае в течение 2.5 нс. Процедура повторяется три раза. После уравновешивания производится моделирование в течении 5 нс. При этом конфигурация системы записывается каждую пикосекунду и полученная траектория используется для получения средних значений и огрубленного потенциала для мезоскопического моделирования.

Получена зависимость плотности системы от концентрации полистирола. Значения плотности находятся между плотностью чистого полистирола (1040 кг/м3) и чистого дихлорэтана (1180 кг/м3) при T=350 K, что подтверждает соответствие результатов моделирования экспериментальным данным. Траектория системы с раствором концентрации 17%w/v и 28%w/v использовались для получения параметров силового поля огрубленной модели. 2.2 Картирование в огрубленную модель

При картировании использовалась модель, предложенная Кремером и Хармандарисом сначала для атактического полистирола и расширенная впоследствии на полистирола различной регулярности [ 12,13,18,39,40 ]. Полимер картируется в линейную цепь бусинок. Каждое мономерное звено состоит из двух бусинок разных типов. Бусинка типа A содержит атомы углерода, соединяющие два последовательных углеродных кольца и атомы водорода, соединенные с ними. CH группа остова, к которой прикреплено фенильное кольцо, принадлежит обеим соседним бусинам типа A. Центр масс бусины A содержит центр масс CH2 группы и две CH группы, взятые с половинным весом. Бусина Б содержит атомы фенильной группы. Цепь представляет собой альтернирующий АБ-сополимер.

Огрубленное силовое поле было получено из атомистической траектории с помощью программного пакета VOTCA [ 41 ]. Огрубленное представление содержит только молекулы полистирола, растворитель учитывается косвенно. Потенциалы на длину связей, валентные углы и углы внутреннего вращения были получены с помощью больцмановского обращения [ 41 ]. Потенциалы парного взаимодействия огрубленной модели были получены с помощью итеративного алгоритма больцмановского обращения [ 41 ]. В качестве теста огрубленного силового поля было проведено сравнения распределений длин связей, валентный углов и углов внутреннего вращения при огрубленном моделировании раствора с полученным силовым полем с аналогичными распределениями, полученными из атомистической траектории. Распределения находятся в хорошем соответствии. Аналогичное сравнении радиальных функций распределения показывает отличное соответствие. 2.3 Синтез полимерного геля

Сшивка молекул полистирола проводилась в рамках огрубленной модели. Сшиватель моделировался сферической бусинкой типа В. Бусинки В свободно двигаются в ячейке моделирования и могут химически связываться с бусинками Б, если расстояние между ними становится меньше критического расстояния rlink = 0.288 нм. Бусинка В может связаться в двумя различными бусинками Б, образуя сшивку между двумя молекулами полистирола.

Моделирование проводилось с помощью пакета программ LAMMPS[42]. Получение силовое поле для молекул полистирола описано выше. Молекулы сшивающего агента В имели размер 0.14 нм и взаимодействовали друг с другом и бусинками А и Б посредством потенциала исключенного объема. После реакции сшивания несвязанные взаимодействия с возникшими 12 и 1-3 соседями заменялись связанными взаимодействиями на длину связи Б-В и угол Б-В-Б, равновесные значения которых были посчитаны в атомистической модели. Потенциалы для углов А-Б-В и В-Б-В устанавливались таким образом, чтобы взаимодействие между соответствующими бусинками не изменилось после замены несвязанного взаимодействия на связанное. Проверка расстояний между бусинками Б и В проводилась каждые 1000 шагов интегрирования. При наличии удовлетворяющих критериям пар интегрирование останавливалось и в структуру системы вносились вышеописанные изменения.

Моделирование проводилось в ячейки с n = 125 молекулами полистирола длиной в N = 20 звеньев. Число молекул сшивателя Nc равнялось Nc = xnN, где x - отношение количества сшивателя к количеству мономерных звеньев. Моделирование проводилось в кубической ячейке с периодическими граничными условиями в каноническом ансамбле. Размер ячейки подбирался таким образом, чтобы концентрация полимера совпадала с его концентрацией в атомистическом моделировании. Процесс сшивки завершался в случае полной конверсии, либо если в течении 105 шагов интегрирования не было образовано ни одной новой сшивки. Затем гель уравновешивался в течение 2∙106 шагов. Над полученными конфигурациями проводилась процедура обратного картирования. 2.4 Обратное картирование

На данном этапе проводилось преобразование сверхсшитого геля в огрубленной модели в атомистическое представление. Сначала восстанавливался углеродный остов из позиций центров масс звеньев типа А. Затем ароматические кольцо прикрепляются к нем в позициях Б звеньев и мостики - атомы углерода вставляются в позициях бусинок типа В. В конце к углеродному остову прикрепляются атомы водорода. Потом, чтобы убрать локальные напряжения и близкие контакты между атомами, проводилась минимизация энергии методом сопряженных градиентов [ 24 ]. Полученная сверхсшитая сетка исследовалась с помощью пакета программ GROMACS с силовым полем, описанным в п. 2.1, при температуре T = 350 K и давлением P = 1 атм.

Все описанные компьютерные эксперименты проводились на суперкомпьютером комплексе МГУ[ 43 ]. 3. Результаты

Была проведена реакция сшивания полистирола при различной концентрации сшивателя x = 0.25, 0.5, 1.0. При концентрации сшивателя x = 0.25 образовалось несколько разделенных кластеров. При оставшихся концентрациях сшивателя образовался единый кластер с перколяцией по всем трем осям, т.е. можно говорить об образовании макроскопического геля. При этом конверсия в обоих случаях была выше 99%, т.е. можно говорить об образовании сверхсшитого геля. Обнаружено, что в процессе сшивания сшивки распределяются неравномерно в объеме системы. Значительное число фенильных колец имеет больше сшивок, чем среднее значение, либо не имеет сшивок вообще. Пространственное распределение сшивок также неоднородно, как и распределение полимера в целом. Можно говорить об эффективном микрофазном расслоении полимера.

Атомистическое моделирование сшитой сетки с растворителем в изобарических условиях не приводит к значительному изменению размера системы, перераспределение полимера также не наблюдается, что служит косвенным подтверждением корректности процедуры картирования/обратного картирования. Сухой гель меняет структуру, теряя ажурность, и значительно уменьшается в размерах. Сравнение распределений валентных углов показывает, что в сухом геле среднее значение находится дальше от равновесного. Можно сказать, что отсутствие растворителя вызывает внутреннее напряжение в сверхсшитом геле, что позволяет объяснить абсорбцию даже термодинамически невыгодных растворителей.

Были вычислены плотность, удельная поверхность и модуль упругости образцов сверхсшитого геля. Удельная поверхность вычислялась как площадь поверхности, доступной растворителю (SASA) [ 44 ] с радиусом частицы 0.14 нм. Модуль упругости определялся по диаграмме деформирования, полученной последовательным одноосным растяжением системы на 1% за шаг с последующей релаксацией, что эквивалентно постоянной скорости истинной деформации в 0.02 нс-1. Для образца с x=0.5 плотность было получено значение плотности 1.19 г/см3 , тогда как в литературе для подобных гелей упоминаются экспериментальные значения от 0.8 г/см3 [ 6 ] до 1.37 г/см3 [ 45 ]. Удельная поверхность этого образца составила 720 м2/г. Для подобных сухих гелей экспериментальные значения варьируются от 480 м2/г до 740 м2/г. Значение модуля упругости составило около 50 МПа. Экспериментально модуль упругости определялся контактной деформацией двух сферических гелей, набухших в толуоле. Полученные значения составили от 50 до 900 МПа [ 46 ]. Экспериментальные данные свидетельствуют о повышении модуля упругости при увеличении концентрации сшивателя x. Данный эффект качественно подтверждается моделированием.

Описанные результаты опубликованы в ряде статей [ 47,48 ]. В настоящее время идет изучение влияние скорости реакции сшивки на свойства получаемых гелей. Литература 8. Y. He, T. R. Lutz, M. D. Ediger, C. Ayyagari, D. Bedrov, G. D. Smith. Macromolecules, 37, 5032 (2004). 11. R. Faller. Phys. Chem. Chem. Phys., 4, 2269 (2002). 14. Q. Sun, R. Faller. Macromolecules, 39, 812 (2006). 26. J. J. P. Stewart. J. Comput. Chem., 10, 209 (1989). 27. J. J. P. Stewart. J. Comput. Chem., 10, 221 (1989). 42. S. J. Plimpton. J. Comp. Phys., 117, 1 (1995). Multiscale modelling approach to property prediction: hypercrosslinked polystyrene.

Alexei Lazutin, Anna Glagoleva, Mikhail Glagolev and Valentina Vasilevskaya Keywords: multiscale modelling, hypercrosslinked polystyrene, mapping/reverse mapping A multiscale modelling algorithm is suggested to simulate the crosslinking of the polystyrene dissolved in dichloroethane by monochlorodimethyl ether. The algorithm comprises consecutive stages: molecular dynamics atomistic simulation of a polystyrene solution, the mapping of atomistic structure onto coarse-grained model, the crosslink formation, the reverse mapping, and finally relaxation of the structure. The calculated values of the specific surface and the elastic modulus are in reasonable quantitative correspondence with experimental data.

References

1. K. Duzek , D. Patterson . J. Polym. Sci. A2 , 6 , 1209 ( 1968 ).

2. A. R. Khokhlov . Polymer, 21 , 376 ( 1980 ).

3. A. Yu. Grosberg , A. R. Khokhlov . Statistical Physics of Macromolecules . AIP Press, New York ( 1996 )

4. A. R. Khokhlov , S. G. Starodubtsev , V. V. Vasilevskaya . Conformational transitions in polymer gels: theory and experiment . In Responsive Gels: Volume Transitions I, ed. by K. Dusek, Springer-Verlag, Berlin., 109 , 123 ( 1993 ).

5. S. V. Rogozhin , V. A. Davankov , M. P. Tsyurupa , Patent USSR 299165 ( 1969 ).

6. M. P. Tsyurupa , V. A. Davankov . Reactive & Functional Polymers , 66 , 768 ( 2006 ).

7. J. Germain , J. M. J. Frechet , F. Svec . SMALL, 5 , 1098 ( 2009 ).

9. S. Rao , X. Li , H. Liang . Macromolecular Research , 15 , 610 ( 2007 ).

10. B. Bayramoglu , R. Faller . Macromolecules, 45 , 9209 ( 2012 ).

12. V. A. Harmandaris , K. Kremer . Macromolecules, 42 , 791 ( 2009 ).

13. D. Fritz , V. A. Harmandaris , K. Kremer , N. F. A. van der Vegt . Macromolecules, 42 , 7579 ( 2009 ).

15. A.V. Lyulin , M.A.J. Michels . Macromolecules, 35 , 1463 ( 2002 ).

16. A.V. Lyulin , N.K. Balabayev , M.A.J. Michels . Macromolecules, 35 , 9595 ( 2002 ).

17. A. V. Lyulin , B. Vorselaars , M. A. Mazo , N. K. Balabaev and M. A. J. Michels . Europhys. Lett., 71 , 618 ( 2005 ).

18. T. Mulder , V. A. Harmandaris , A. V. Lyulin , N. F. A. van der Vegt , K. Kremer , M. A. J. Michels . Macromolecules, 42 , 384 ( 2009 ).

19. J. W. Liu, M. E. Mackay , P. M. Duxbury . Europhys. Lett., 84 , 46001 ( 2008 ).

20. J. W. Liu, M. E. Mackay , P. M. Duxbury . Macromolecules, 42 , 8534 ( 2009 ).

21. F. Ferrante , F. L. Celso , D. Duca . Colloid Polym. Sci., 290 , 1443 ( 2012 ).

22. C. Oostenbrink , A. Villa , A. E. Mark , W. F. van Gunsteren. J. Comput . Chem., 25 , 1656 ( 2004 ).

23. A. K. Malde , L. Zuo , M. Breeze , M. Stroet , D. Poger , P. C. Nair , C. Oostenbrink , A. E. Mark . Journal of Chemical Theory and Computation , 7 , 4026 ( 2011 ).

24. R. Fletcher , C. M. Reeves . The Computer Journal , 7 , 149 ( 1964 ).

25. M. J. S. Dewar , E. G. Zoebisch , E. F. Healy , J. J. P. Stewart . Journal of the American Chemical Society , 107 , 3902 ( 1985 ).

28. E. R. Davidson , D. Feller . Chem. Rev., 86 , 661 ( 1988 ).

29. H. J. C. Berendsen , D. van der Spoel , R. van Drunen. Comp. Phys. Comm. , 91 , 43 ( 1995 ).

30. E. Lindahl , B. Hess , D. van der Spoel. J. Molecular Mod ., 7 , 306 ( 2001 ).

31. D. van der Spoel , E. Lindahl, B. Hess , G. Groenhof , A. E. Mark , H. J. Berendsen , H. J. J. Comp . Chem., 26 , 1701 ( 2005 ).

32. B. Hess , C. Kutzner , D. van der Spoel, E. Lindahl. J. Chem . Theory Comput. , 4 , 435 ( 2008 ).

33. S. Pronk , S. Páll , R. Schulz , P. Larsson , P. Bjelkmar , R. Apostolov , M. R. Shirts , J. C. Smith , P. M. Kasson , D. van der Spoel, B. Hess , E. Lindahl . Bioinformatics, 29 , 845 ( 2013 ).

34. U. Essmann , L. Perera , M. L. Berkowitz , T. Darden , H. Lee , L. G. Pedersen. J. Chem . Phys., 103 , 8577 ( 1995 ).

35. G. Bussi, D. Donadio , M. Parrinello . J. Chem. Phys. , 126 , 014101 ( 2007 ).

36. H. J. C. Berendsen , J. P. M. Postma , W. F. van Gunsteren , A. DiNola , J. R. Haak . Journal of Chemical Physics , 81 , 3684 ( 1984 ).

37. S. Kirkpatrick , C. D. Gelatt , M. P. Vecchi . Science, 220 , 671 ( 1983 ).

38. D. Frenkel , B. Smit . Understanding Molecular Simulation: From Algorithms to Applications. 2nd Edition , Academic Press: San Diego ( 2002 ).

39. V. A. Harmandaris , N. P. Adhikari , N. F. A. van der Vegt , and K. Kremer . Macromolecules, 39 , 6708 ( 2006 ).

40. V. A. Harmandaris , D. Reith , N. F. A. van der Vegt , K. Kremer . Macromol. Chem. Phys. , 208 , 2109 ( 2007 ).

41. V. Rühle , C. Junghans , A. Lukyanov , K. Kremer , D. Andrienko . J. Chem. Theory Comput. , 5 , 3211 ( 2009 ).

43. V. Sadovnichy , A. Tikhonravov , Vl. Voevodin, V. Opanasenko . "Lomonosov" : Supercomputing at Moscow State University. In Contemporary High Performance Computing: From Petascale toward Exascale , 283 , Boca Raton , USA, CRC Press ( 2013 ).

44. F. Eisenhaber , P. Lijnzaad , P. Argos , C. Sander , M. Scharf . J. Comput. Chem ., 16 , 273 ( 1995 ).

45. M. Xu , Z. Shi , B. He. Chinese J. React . Polym., 2 , 119 ( 1993 ).

46. A. V. Pastukhov . Physicochemical properties and structural mobility of hypercrosslinked polystyrenes . Doctoral thesis , A. N. Nesmeyanov Institute of Elementoorganic Compounds ( 2008 ).

47. A.A. Lazutin , M.K. Glagolev , V.V. Vasilevskaya , Khokhlov A.R. Hypercrosslinked polystyrene networks: An atomistic molecular dynamics simulation combined with a mapping/reverse mapping procedure // J. Chem. Phys . - 2014 . - V. 140 - P. 134903

48. M.K. Glagolev , А.А. Lazutin , V.V. Vasilevskaya Macroscopic properties of hypercrosslinked polystyrene networks: an atomistic and coarse-grained molecular dynamics simulation // Macromolecular Symposia - 2015 . - V. 348, Is . 1 - P. 14 - 24