Высокомолекулярные соединения (серия А), 2019, T. 61, № 5, стр. 466-480

ВВЕДЕНИЕ

Мембранное разделение – сравнительно новое направление промышленного разделения газовых смесей, быстро развивающееся благодаря ряду технологических преимуществ по сравнению с традиционными методами разделения газов (ректификацией, абсорбцией и адсорбцией) [1–4]. К достоинствам мембранного газоразделения относятся низкое энергопотребление, отсутствие фазовых переходов, малая материалоемкость, гибкость управления и модульность структуры технологических установок. Наибольшее распространение получили промышленные процессы мембранного газоразделения для выделения азота из воздуха, удаления углекислого газа из “кислых” природных газов и извлечение водорода из различных потоков химии и нефтепереработки.

Но в основе всего лежит квалифицированный выбор материала мембран и возможность приготовления мембран, которые были бы при этом и высокопроизводительными и селективными. В подавляющем большинстве процессов мембранного газоразделения используются полимеры. Синтез или выбор полимера и изучение его свойств и поведения в процессе газоразделения является предметом мембранного материаловедения. Важнейшая задача мембранного материаловедения – поиск новых высокоэффективных полимеров для решения разнообразных задач мембранной технологии. Для мембранного газоразделения решающую роль играют транспортные свойства полимерного материала – проницаемость ${{P}_{i}}$ и селективность ${{P}_{i}}{\text{/}}{{P}_{i}}$ для пары разделяемых газов $i$ и $j$. Возможность хотя бы приблизительно предсказать транспортные свойства гипотетического полимера еще до его синтеза существенно сокращает усилия и экономит время, позволяя отбросить заведомо неудачные варианты и сконцентрироваться на наиболее перспективных соединениях.

Традиционно важнейшим фактором, влияющим на транспортные свойства аморфного полимера, считается его химическое строение. Химическая природа мономерных звеньев и способы их объединения в полимерную цепь определяют взаимное расположение цепей и геометрию свободного объема внутри образца. Именно структура свободного объема полимерной матрицы определяет коэффициент проницаемости P_i молекул газа i через мембрану.

Как известно, коэффициент проницаемости P можно записать в виде произведения коэффициента растворимости S и коэффициента диффузии D:

Первый описывает движущую силу процесса переноса молекул газа, второй соответствует кинетической компоненте процесса. Во многих случаях именно коэффициент растворимости S, а также селективность растворимости ${{S}_{i}}{\text{/}}{{S}_{i}}$ по паре газов i и j определяют транспортные свойства и ценность полимера для мембранной технологии. Это обусловливает актуальность задачи предсказания S для вновь создаваемых полимерных материалов.

Согласно модели двойной сорбции [5], общее количество C газа, растворенного в объеме аморфного стеклообразного полимера в зависимости от давления p, описывается комбинацией закона Генри и сорбции ленгмюровского типа по формуле

где k_D – константа равновесия закона Генри, b – константа равновесия лэнгмюровской сорбции, $С_{H}^{'}$ – лэнгмюровская сорбционная емкость. Коэффициент растворимости при бесконечном разбавлении S вычисляется как отношение C/p при p → 0, т.е. из изотермы сорбции (2), S можно представить в виде суммы двух слагаемых: S = = k_D + $С_{H}^{'}b$.

В настоящей статье развивается подход [1] к предсказанию коэффициента растворимости S легких газов в стеклообразных аморфных полимерах, основанный на компьютерном моделировании конформаций короткой полимерной цепи. При той же точности и сравнимом числе объясняющих переменных усовершенствованный метод позволил устранить ряд недостатков подхода [1], а также предсказать константу закона Генри с относительной точностью 47%.

Как и в работе [1], для предсказания транспортных свойств генерируются несколько конформаций короткого отрезка цепи рассматриваемого полимера (обычно включающего несколько мономерных звеньев и от 200 до 600 атомов) путем релаксации в эмпирическом поле Дрейдинга из случайных начальных позиций атомов.

Виртуальные шаровые модели полученных конформеров “обкатываются” пенетрантом различного радиуса в диапазоне от 0 до 3 Å и строится зависимость площади доступной поверхности (площади контакта макромолекулы с пенетрантом) от радиуса обкатки. Параллельно вычисляются еще несколько геометрических химических индексов (см. детали в разделе 4). В отличие от [1], вместо эмпирического алгоритма [6], реализованного в Molecular Surface Area Plugin продукта ChemAxon, для “обкатки” используется (более точный и прозрачно реализованный) классический алгоритм Ли–Ричардса [7].

Полученные зависимости линеаризуются на диапазоне радиусов обкатки, который в отличие от работы [1] был адаптивно подобран для получения наилучших результатов. Линеаризованная зависимость значения i-го индекса от радиуса обкатки R имеет вид a_i + b_iR, и пары коэффициентов a_i, b_i используются в качестве основных объясняющих переменных для предсказания логарифма как растворимости S, так и константы равновесия закона Генри k_D.

Вес объясняющих переменных для предсказания lgS был подобран на основе того же, что и в работе [1] массива из 344 экспериментальных измерений (плюс еще 134 измерения были отобраны для тестовой выборки) 64 полимеров различных классов из Базы данных физико-химических свойств полимеров Лаборатории мембранного газоразделения ИНХС РАН. Регрессия для lgk_D настраивалась на основе 76 экспериментальных измерений (плюс 17 измерений, входящих в тестовую выборку) 40 полимеров.

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

ОБЗОР ЛИТЕРАТУРЫ

В литературе известно несколько подходов к предсказанию транспортных свойств полимеров. Наиболее хорошо разработанные и популярные основаны на атомистическом компьютерном моделировании (см. обзоры в [8, 9]). В частности, для предсказания коэффициентов проницаемости, диффузии и растворимости газов в полимерах использовались методы молекулярной динамики и большого канонического ансамбля [10, 11].

В методе молекулярной динамики [12] атомы перемещаются в соответствии с уравнениями движения Ньютона в эмпирическом силовом поле, моделирующим взаимодействия между молекулами газа и атомами полимерной цепи. Однако прямое применение метода молекулярной динамики к предсказанию диффузии невозможно, поскольку требует моделирования на колоссальном временном отрезке. Теория переходного состояния [13] позволила преодолеть эти трудности, что обусловило популярность данного подхода в моделировании процессов диффузии. В теории переходного состояния и ее многочисленных вариантах [14, 15] используется иерархия взаимосвязанных моделей процесса диффузии, причем каждая модель действует на своем временном и пространственном уровне [8], от поправочных членов к эмпирическому силовому полю до геометрического анализа границ “ям” потенциальной энергии (так называемых стабильных состояний) и марковской сетевой модели случайных переходов частицы газа между этими стабильными состояниями.

В отличие от метода молекулярной динамики широко используемого при моделировании диффузии, методы большого канонического ансамбля показали свою эффективность при моделировании растворимости. В тесте Видома [8] коэффициент растворимости вычисляется на основе потенциальной энергии молекулы газа, помещенной в случайное место полимерной фазы.

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

Стоит отметить, что в большинстве работ (за исключением, пожалуй, недавней статьи [16]), предполагается заранее известной плотность стеклообразной полимерной фазы, что делает проблематичным применение такого моделирования в рамках молекулярного дизайна, когда требуется моделировать гипотетические, ранее не синтезированные или не изученные полимеры.

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

Совершенно другой подход к предсказанию транспортных характеристик полимерных материалов, основанный на использовании топологических индексов и метода эмпирических групповых вкладов был предложен в [17, 18]. В методе эмпирических групповых вкладов каждому уникальному подграфу молекулярного графа (“группе”) присваивается такой статистический вес, чтобы взвешенная сумма числа вхождений в полимерное звено всевозможных молекулярных групп наилучшим образом коррелировала со значением интересующего физико-химического показателя на некоторой выборке экспериментальных данных. В частности, была показана эффективность метода для предсказания коэффициентов проницаемости и диффузии широкого круга газов в аморфных полимерах. К недостаткам стоит отнести относительно слабую предсказательную силу (способность к обобщению за пределы обучающей выборки), что является следствием сложности и нелинейности связи между топологией полимерного звена и структурой свободного объема полимерного образца.

Отметим еще несколько работ, в которых рассматривались родственные проблемы [19–24]. Монографию J. Bicerano [19], в которой декларировалось применение теории графов, но в действительности при вычислении конкретных параметров использовались подгоночные параметры, строго говоря, можно было бы исключить из рассмотрения. Важный вклад в вычисление многочисленных физико-химических параметров полимеров был внесен А. Аскадским [20, 21]. Однако среди большого числа предсказанных параметров (даже коэффициентов проницаемости для некоторых газов) коэффициенты растворимости газов отсутствуют. Метод, родственный используемому в настоящей работе, был применен И. Роновой с сотрудниками [22–24]. Он включает обкатку полимерных цепей сферическими зондами с размерами, моделирующими размеры молекул газовых пенетрантов. Но в этих работах рассматривались только корреляции экспериментальных значений Р и D с доступным свободным объемом у макромолекулярных цепей и никаких предсказаний не делалось. Такие корреляции были проанализированы в том числе и в работе [23], где они были построены и для коэффициентов растворимости, найденных как отношения экспериментальных параметров P/D.

В настоящей статье развивается относительно новый подход, впервые представленный в работе [1], который обладает следующими преимуществами: не требует экспериментальных данных о полимерном материале; применим к широкому кругу полимеров, в том числе еще не синтезированным; предъявляет весьма скромные требования к вычислительным мощностям.

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

Первый этап – моделирование короткого отрезка полимерной цепи и использование площади доступной поверхности идейно наиболее близок подходу [25]. Второй этап – использование статистических инструментов для обучения линейной регрессии ближе к подходу [17, 18]. В то же время есть и многочисленные отличия, обусловливающие научную новизну подхода настоящей статьи: использование большого числа (а именно, девяти) различных геометрических индексов, идея о том, что важными объясняющими переменными могут стать параметры зависимости площади доступной поверхности (и родственных ей индексов) от радиуса “обкатки”, уточненный (по сравнению с работой [26]) общий вид зависимости транспортных характеристик от параметров пенетранта и полимера, статистически корректная процедура отбора значимых переменных для линейной регрессии, позволяющая избежать эффекта переобучения.

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

МОЛЕКУЛЯРНО-ДИНАМИЧЕСКОЕ МОДЕЛИРОВАНИЕ

Подход к моделированию полимерной цепи в настоящей работе не претерпел изменений по сравнению с опубликованным в работе [1]. Как и ранее, для каждого полимера в среде молекулярного моделирования Instant JChem создается цепочка из нескольких мономерных звеньев, включающая от 200 до 600 атомов. Для этой цепочки средствами Conformer Plugin [27] пакета Instant JChem ChemAxon [28] генерируется 6 случайных конформаций различной геометрии, получаемой оптимизацией в эмпирическом силовом поле Дрейдинга из различных начальных позиций атомов. В следующем разделе обосновывается достаточность выбранной длины цепочки и количества генерируемых конформеров для устойчивости получаемых далее результатов.

Ниже приведено несколько полученных таким образом типичных конформаций для ряда “популярных” стеклообразных полимеров из Базы данных физико-химических свойств полимеров Лаборатории мембранного газоразделения ИНХС им. А.В. Топчиева [29].

ChemAxon Instant JChem считается довольно “всеядным” программным пакетом, не имеющим узкой специализации по классам моделируемых веществ. В то же время ряд построенных конформаций выглядел весьма нереалистично. Эти полимеры были пересчитаны в пакете PerkinElmer Chem3D версии 15.1.0.144 [30] по аналогичной методике; для задания случайного начального положения атомов к моделируемой полимерной цепи проводилось молекулярно-динамическое моделирование при температуре от 300 до 3000 K в течение 1000 итераций, после чего полученная структура оптимизировалась по свободной энергии в эмпирическом поле MM2.

Описанный подход позволяет достаточно точно моделировать локальную конфигурацию полимерной цепи, порождаемую взаимодействием соседних мономерных звеньев, но игнорирует влияние на форму молекулы далеко расположенных частей цепи и, что существенно, других полимерных цепей. Следовательно, от такого подхода не стоит ожидать точного предсказания ни свободного объема, ни распределения объема полостей по размерам, ни связности полостей – факторов, которые, как считается [31], сильно влияют на коэффициент проницаемости, но скорее через его кинетическую, диффузионную составляющую (см. формулу (1)).

ВЫЧИСЛЕНИЕ ГЕОМЕТРИЧЕСКИХ ИНДЕКСОВ

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

Площадь доступной поверхности (accessible surface area, ASA) молекулы – это площадь поверхности, описываемой центром шарообразного “зонда” заданного радиуса во всевозможных положениях его касания с ван-дер-ваальсовой поверхностью (van der Waals surface area VWDSA) этой молекулы. Иллюстрация к построению этой поверхности для простой молекулы из нескольких атомов приведена ниже.

ASA является геометрическим индексом [33], т.е. функцией, которая каждой конформации молекулы (пространственному расположению ее атомов) ставит в соответствие число, зависящее только от взаимного расположения атомов, но не от положения молекулы в пространстве. По сравнению с похожим индексом – площадью поверхности контакта (площадью Конолли) – ASA макромолекулы имеет преимущество: именно площадь доступной поверхности пропорциональна вероятности контакта между макромолекулой (полимера) и хаотически движущимся “зондом”, представляющим молекулу газа-пенетранта.

ASA для радиуса “зонда” 1.4 Å (соответствующего воде) широко используется в биохимии для изучения взаимодействия белковой макромолекулы с растворителями. В мембранной технологии газы-пенетранты отличаются своим эффективным радиусом, поэтому для целей мембранного газоразделения важно знать ASA макромолекул аморфных полимеров для различных радиусов “обкатки” R. Обкатка полимерных цепей сферическими зондами радиуса R, моделирующего размер молекулы газового пенетранта, применялась И. Роновой с сотрудниками [22–24] для вычисления доступного свободного объема (FAV) макромолекулярных цепей, однако, насколько нам известно, кривые зависимости ASA(R) полимерных материалов в задачах мембранной науки и технологии никогда не исследовались. В настоящей статье, как и в работе [1], именно параметры этой кривой и аналогичных ей позволяют предсказать коэффициент растворимости газа и константу равновесия закона Генри. В биохимии также широко используются несколько родственных ASA геометрических индексов (табл. 1) [33].

Перечисленные в табл. 1 индексы являются функциями радиуса “обкатки” R. Для построения зависимостей ASA(R), ASA⁺(R), и т.д. значение индекса вычислялось для каждого R от 0 до 3 Å с шагом 0.05 Å. Значение геометрического индекса зависит от конформации макромолекулы. Для получения стабильной характеристики полимерного материала “вообще”, не зависящей от конкретной конформации, значения каждого индекса для каждого радиуса “обкатки” усреднялись по шести различным конформациям макромолекулы. В работе [1] для десяти случайным образом выбранных полимеров из исследуемой выборки была построена зависимость средней удельной площади доступной поверхности ASA от числа усредняемых конформаций. Было показано (см. рис. 4 в работе [1]), что вариации значения ASA не превышают 1% при усреднении шести конформаций и более, на основании чего принимается, что шесть конформаций – достаточно большая выборка для обеспечения инвариантности результата усреднения.

Аналогичным образом в работе [1] проводилось исследование зависимости ASA от длины моделируемой “олигомерной” цепи. Было установлено (см. рис. 5 в работе [1]), что для 10 случайным образом выбранных полимеров значение ASA (а также и остальных вычисляемых индексов) варьируется в пределах нескольких процентов при длине олигомерной цепи более 200 атомов. Отдельным интересным направлением будущих исследований является изучение зависимости минимальной допустимой длины олигомерной цепи от длины сегмента Куна.

В первой версии нашего метода [1] для вычисления площади доступной поверхности использовался эмпирический алгоритм [6], реализованный в Molecular Surface Area Plugin продукта Chemaxon Instant JChem. Опыт его использования, однако, показал, что помимо преимуществ (высокой скорости расчета и корректного учета групповых пересечений) этот алгоритм имеет и ряд недостатков. В частности, оказалось, что радиусы атомов и стандартных остатков специальным образом откалиброваны, в первую очередь, для белковых молекул в водной среде, что делает вычисленные значения индекса менее предсказуемыми и объяснимыми.

В настоящей статье описывается модифицированный подход, в котором мы используем собственную реализацию (на языке Python) классического алгоритма Ли–Ричардса. Помимо большей точности (достигаемой, впрочем, за счет увеличения трудоемкости) такая реализация алгоритма позволяет находить вклад каждого атома в площадь доступной поверхности для вычисления индексов табл. 1, зависящих от частичных зарядов. На рис. 1 приведены несколько типичных зависимостей геометрических индексов от радиуса обкатки и их сравнение с индексами, рассчитанными с помощью эмпирического алгоритма [6].

Помимо отличия между алгоритмами, рисунок иллюстрирует разнообразие кривых. Имея общую тенденцию к возрастанию с увеличением R (что естественным образом следует из определения), зависимость площади доступной поверхности ASA от радиуса “обкатки” R может быть монотонной, иметь точку минимума, перегиб. Вид кривых для родственных индексов, зависящих от частичных зарядов, может быть и существенно более сложным, определяясь доступностью положительно/отрицательно заряженных участков молекулярной поверхности для “зондов” разного радиуса.

Совокупность этих восьми кривых дает довольно информативное описание поверхности макромолекулы с учетом частичных зарядов. Детальное исследование связи их вида с физико-химическими свойствами вещества является предметом будущих исследований. В настоящей работе мы будем пользоваться лишь их (достаточно грубой) линейной аппроксимацией. А именно, зададимся “эффективным” промежутком радиусов [R^–, R⁺], на котором построим методом наименьших квадратов лучшую аппроксимацию кривых для всех восьми индексов из табл. 1. Процесс выбора оптимального диапазона аппроксимации будет подробно описан в следующем разделе.

Линейная аппроксимация i-го индекса описывается парой числовых коэффициентов (рис. 2): константа c_i соответствует значению линеаризованного индекса при нулевом радиусе, тангенс угла наклона d_i описывает чувствительность (или селективность) индекса к радиусу пенетранта (имея, как позже будет показано, прямое отношение к селективности растворимости полимера). Именно эти числовые коэффициенты, описывающие характер взаимодействия макромолекулы полимера с пенетрантами различного радиуса, используются далее в качестве объясняющих переменных для предсказания транспортных характеристик полимера наряду с удельным (на единицу объема полимера) ван-дер-ваальсовым объемом VDWV (табл. 1). Последний описывает объемные характеристики макромолекулы и считается одним из важных детерминантов свободного объема полимерной матрицы [34].

ПРЕДСКАЗАНИЕ КОЭФФИЦИЕНТА РАСТВОРИМОСТИ

В работе [31] была обоснована линейная зависимость логарифма коэффициента растворимости (при бесконечном разбавлении) газа в полимере от площади доступной поверхности молекулы газа-пенетранта. В работе [1] было показано, что MaxPA – максимальная площадь проекции молекулы газа (где атомы представлены своими ван-дер-ваальсовыми сферами) – позволяет еще лучше описать линейную зависимость логарифма коэффициента растворимости газа g в полимере p для большинства популярных в мембранной технологии легких газов (за исключением окиси углерода, демонстрирующей девиантное поведение):

При этом коэффициенты a и b линейной зависимости (3) являются функциями полимера p, отличаясь для различных полимеров. Задача настоящей статьи состоит в выявлении зависимости коэффициентов a и b от введенных в предыдущем разделе характеристик полимера – коэффициентов c_i, d_i линейных аппроксимаций геометрических индексов макромолекулы полимера (см. табл. 1), а также от ван-дер-ваальсового объема VDWV макромолекулы.

Как и в работе [1], зависимость a и b от c_i, d_i (i = ASA, …, PNSA₃) предполагается линейной:

Числовые коэффициенты a₀, …, a₁₇, b₀, …, b₁₇ подбирали с помощью регрессионного анализа на основе экспериментальных данных из упомянутой Базы данных. Для обеспечения сопоставимости результатов использовали ту же выборку, что в предыдущей версии метода (см. работу [1]) с учетом удаления трех обнаруженных дублирующихся измерений. В результате осталось 478 уникальных измерений (344 в обучающей и 134 в тестовой выборке) по растворимости 18 различных легких газов (самым тяжелым из которых был бутан, CO в выборку не входил) в 64 полимерах различных классов (в основном полиимидов, полистиролов, полиакрилатов, виниловых полимеров). Коэффициент растворимости вычисляли по формуле (1) как отношение S = P/D экспериментально измеренных коэффициентов проницаемости и диффузии. Коэффициенты, полученные при разной температуре, приводили к единой температуре (308 K) по методике, описанной в работе [35]. Для расчета удельных (на единицу объема) значений геометрические индексы нормализовали на экспериментальные плотности (средние плотности, при наличии более одного источника) соответствующих образцов полимеров.

Подробное описание набора данных приведено в Дополнительных материалах, а критический анализ методики их подготовки описан в работе [1].

В связи с большим числом объясняющих переменных при небольшом объеме выборки для предотвращения переобучения использовали процедуру отбора переменных [36] на языке Python, реализующую алгоритм шаговой регрессии [37] на основе статистики Фишера F: переменная с наименьшей вероятностью p(F) включается в регрессию, если p(F) меньше заданного порога  f_in, переменная исключается, если для нее p(F) больше заданного порога  f_out.

Выбор оптимального диапазона [R^–, R⁺] линеаризации кривых ASA, …, PNSA₃ и оптимальных параметров шаговой регрессии f_in, f_out производился 5-блочной кросс-валидацией; а именно, обучающая выборка случайным образом разбивалась на пять частей и для каждой комбинации R^– < R⁺ в диапазоне от 0 до 3 Å, а также f_in < f_out в диапазоне от 0.01 до 0.1 строились пять шаговых регрессий на основе обучающей выборки с исключенным первым, вторым, …, пятым блоком и вычислялась средняя корреляция по исключенным блокам. В результате наилучшим оказался диапазон [2.2, 2.6] Å и параметры f_in = 0.06, f_out = 0.09, давшие в результате кросс-валидации максимальную среднюю корреляцию.

Данным значениям параметров соответствовала регрессия (настроенная уже на полной обучающей выборке):

задействующая девять объясняющих переменных, в том числе, значимыми оказались шесть различных параметров полимера. (Обучающая выборка: N = 344; R² = 0.85; скорректированный коэффициент детерминации q² = 0.84, RMSE = = 0.11, F = 207. Тестовая выборка: N = 134, R² = = 0.89, RMSE = 0.11). Значения lgS в обучающей выборке, а также остатки распределены нормально.

Таким образом,

Переменные в формулах (7) и (8) следуют в порядке убывания значимости, при этом наиболее значимыми оказались ${{d}_{{AS{{A}_{P}}}}}$ и ${{c}_{{DPS{{A}_{3}}}}}$, присутствующие, вдобавок, в обоих коэффициентах.

Коэффициент детерминации регрессии составил 0.85 (скорректированный R² = 0.84) на обучающей выборке (при корреляции 0.92) и 0.80 на тестовой выборке (при корреляции 0.89), что соответствует средней относительной ошибке предсказания S 56% на обучающей выборке и 77% на тестовой выборке. В результате модифицированный алгоритм дает в среднем ту же точность что и первая версия метода (70% на тестовой выборке), что считается неплохой точностью с учетом большой (до десятичного порядка) погрешности экспериментального измерения S.

На рис. 3 представлена диаграмма рассеяния для регрессии (6); в табл. 2 и 3 приведена средняя относительная ошибка предсказания S для различных классов полимеров и газов-пенетрантов соответственно.

ПРЕДСКАЗАНИЕ КОНСТАНТЫ ЗАКОНА ГЕНРИ

Как показано в предыдущем разделе, усовершенствованный подход [1] позволяет с удовлетворительной точностью предсказать коэффициенты растворимости около 20 газов в примерно 60 стеклообразных полимерах разной структуры. Межу тем, в этом подходе почти полностью игнорируется влияние свободного объема на коэффициенты растворимости, хотя известно, что многие полимеры, в том числе представляющие интерес в качестве мембранных материалов, отличаются повышенным свободным объемом и бóльшими коэффициентами растворимости, чем обычные стеклообразные полимеры [38–40].

Мерой свободного объема может служить ленгмюровская сорбционная емкость в модели двойной сорбции $C_{H}^{'}$ [32]. В рамках данной модели, по формуле (2) коэффициент растворимости при бесконечном разбавлении (низких давлениях сорбата) может быть представлен как S = k_D + + $C_{H}^{'}b$, где k_D – коэффициент растворимости “популяции Генри” молекул сорбата, т.е. параметр, характеризующий растворимость в более плотных областях полимерной матрицы, подчиняющийся линейной изотерме Генри, а параметр b является константой равновесия для молекул сорбатов в элементах свободного объема. Для высокопроницаемых полимеров выполняется неравенство $C_{H}^{'}b \gg {{k}_{D}}$. Таким образом, параметр k_D не должен зависеть от свободного объема в полимере. В связи с этим представляет интерес применение описанного в статье подхода для предсказания константы закона Генри k_D.

Был собран набор экспериментальных значений k_D из 94 измерений. Набор построен на основе анализа литературы и включает в себя данные по растворимости 13 газов в 40 полимерах. В табл. 2 и 3 приведена статистика по классам полимеров и газам-пенетрантам. Проблемой является то, что эксперименты проводились при разной температуре от 278 до 308 К. В большинстве работ отсутствовала информация о температурной зависимости k_D, и использованные в работе значения k_D являются усредненными в указанном температурном интервале.

К такому набору данных был применен описанный выше подход. Для каждого полимера из набора данных был построен короткий отрезок цепи размером от 200 до 700 атомов и рассчитано 6 конформаций средствами ChemAxon Instant JChem и PerkinElmer Chem3D версии 15.1.0.144. Для каждой конформации были вычислены значения геометрических индексов из табл. 1 для радиуса обкатки от 0 до 3 Å (за исключением VDWV, к которому понятие радиуса “обкатки” неприменимо). В связи с отсутствием экспериментальных величин плотности образцов значения индексов не нормализовались на плотность (что равнозначно предположению о средней плотности исследуемых стеклообразных полимеров 1 г/см³). В табл. 4 показано распределение числа образцов по классам полимеров, а в табл. 5 – по газам пенетрантам. Отметим отсутствие среди исследуемых газов CO (как и для растворимости), а также CO₂, который, предположительно, также демонстрирует индивидуальный характер.

Набор данных был разбит на обучающую (76 измерений) и тестовую (17 измерений) выборки. Предполагалось, что искомая регрессия имеет вид, аналогичный выражению (3):

где коэффициенты a и b зависят от полимера. Как и ранее, они подбирались в виде взвешенной суммы коэффициентов линейной аппроксимации кривых ASA(R), …, PNSA₃(R) на том же что для растворимости S отрезке [2.2 Å; 2.6 Å]. Также в качестве объясняющей переменной использовался коррелянт свободного объема, ван-дер-ваальсов объем VDWV макромолекулы (для проверки высказанного выше предположения о том, что для предсказания k_D показатели свободного объема менее значимы).

В результате применения шаговой регрессии (с теми же параметрами что и в предыдущем разделе) была получена формула

задействующая единственную переменную, характеризующую полимер – это чувствительность индекса PNSA₃ (площади отрицательно заряженной доступной поверхности с учетом величины частичного заряда) к радиусу “обкатки” (Обучающая выборка: N = 76, R² = 0,79; скорректированный коэффициент детерминации q² = 0.79, RMSE = 0.12, F = 155. Тестовая выборка: N = 17, R² = 0.81, RMSE = 0.09). Значения lgk_D в обучающей выборке, а также остатки распределены нормально.

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

Коэффициент детерминации регрессии составил 0.79 (скорректированный R² = 0.78) на обучающей выборке (при корреляции 0.89) и 0.81 на тестовой выборке (при корреляции 0.9), что соответствует средней относительной ошибке предсказания k_D 69% на обучающей выборке и 47% на тестовой выборке. Таким образом, несмотря на меньший объем выборки, отсутствие нормировки объясняющих переменных на экспериментальную плотность образца, а также на отсутствие нормализации значений k_D по температуре, удалось получить даже лучшую точность предсказания, чем для коэффициента растворимости S.

На рис. 4 приведена диаграмма рассеяния для регрессии (9), в табл. 4 и 5 – средняя относительная ошибка предсказания k_D для различных классов полимеров и газов-пенетрантов соответственно.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Изучение формулы (6) позволяет сделать несколько общих выводов относительно связи между геометрией полимерной цепи и коэффициентом растворимости S. Важен состав значимых характеристик полимера, а также знак и абсолютная величина коэффициентов при этих характеристиках, определяющие их влияние на процессы сорбции. Из формулы (3) легко видно, что поскольку для двух газов, g и g', ln[S(p, g)/S(p, g')] = b(p)[MaxPA(g) – MaxPA(g')], лишь коэффициент b(p) влияет на селективность растворимости.

Тогда из формулы (6) следует, что селективность растворимости определяется тремя значимыми параметрами – чувствительностью ${{d}_{{AS{{A}_{P}}}}}$ площади поляризованной доступной поверхности к радиусу обкатки, чувствительностью ${{d}_{{PNS{{A}_{3}}}}}$ площади отрицательно заряженной поверхности, а также площадью доступной поверхностью с учетом знака ${{c}_{{DPS{{A}_{3}}}}}$ (все это при линеаризации кривых зависимостей индекса от радиуса обкатки на диапазоне 2.2–2.6 Å). Интересно, что веса всех трех показателей устроены таким образом, что они всегда вносят отрицательный вклад в значение коэффициента b(p), который, тем не менее, всегда остается положительным (на рассматриваемом наборе данных он меняется от 0.1 до 0.3 lg[см³(н.у.)/см³(см.рт.ст.)]/Å). Так, полиимиды имеют относительно большую чувствительность ${{d}_{{AS{{A}_{P}}}}}$ поляризованной поверхности и соответственно меньшую селективность растворимости. Это довольно неожиданный эффект. Наоборот, ожидаемым является то, что переменные ${{d}_{{AS{{A}_{P}}}}}$ и ${{d}_{{PNS{{A}_{3}}}}}$, влияющие также и на коэффициент a(p), входят в формулу (7) с весами противоположного знака.

Стоит отметить, что ван-дер-ваальсов объем макромолекулы, коррелянт свободного объема, не оказался значимым фактором, влияющим на селективность растворимости. В то же время он оказался среди факторов, влияющих на коэффициент a(p) – “коэффициент растворимости при нулевом радиусе пенетранта”, влияющий на абсолютную величину коэффициента растворимости.

Интересно, что наилучший диапазон линеаризации [2.2; 2.6] Å оказался за пределами диапазона изменения эффективных (по Теплякову–Мирсу) радиусов пенетрантов, что может быть признаком того что с точки зрения процесса сорбции молекулы газа имеют больший эффективный радиус. Впрочем, потеря качества предсказания при использовании других диапазонов линеаризации невелика, и все они более или менее информативны.

Регрессия (10) задействует чувствительность индекса PNSA₃ (площади отрицательно заряженной доступной поверхности с учетом величины частичного заряда) к радиусу “обкатки”. При этом увеличение чувствительности индекса (в выборке присутствуют в основном положительные значения ${{d}_{{PPS{{A}_{3}}}}}$) приводит к уменьшению значения k_D.

Методика предсказания коэффициента растворимости и константы закона Генри на основе только компьютерных моделей конформаций макромолекул, описанная в настоящей статье, может быть использована при поиске и синтезе полимерных материалов с экстремальными транспортными характеристиками.

Для нового (в том числе еще не синтезированного) полимера в программной среде ChemAxon Instant JChem необходимо создать полимерную цепочку размером 200–600 атомов, рассчитать не менее шести конформаций средствами пакета Conformers Plugin среды Instant JChem или соответствующими средствами приложения ChemOffice Pro и вычислить геометрические индексы табл. 1. При этом достаточно ограничиться индексами, входящими в формулу (6) для радиусов обкатки от 2.2 до 2.6 Å и найти для них коэффициенты c_ASA, d_ASA, … наилучшей линейной аппроксимации по радиусу обкатки R.

После этого для вычисления предсказанного значения lgS достаточно подставить вычисленные значения индексов c_ASA, d_ASA, …, VDWV в формулу (6) и использовать полученное значение растворимости при отборе потенциально интересных для синтеза полимерных материалов с требуемыми транспортными характеристиками. Аналогично используется и формула для предсказания константы равновесия закона Генри k_D.

Исследование частично выполнено при финансовой поддержке Российского фонда фундаментальных исследований в рамках научного проекта 17-08-00164.

О.А. Милосердов выражает благодарность за финансовую поддержку Российскому фонду фундаментальных исследований в рамках научного проекта 18-37-00265.

Работа частично выполнена в рамках государственного задания ИНХС РАН.