Scientific journal
Modern problems of science and education
ISSN 2070-7428
"Перечень" ВАК
ИФ РИНЦ = 1,006

THE METHOD FOR CALCULATING THE EFFECTIVE CHARGE OF SPHERICAL COLLOIDAL PARTICLE

Kukhtetskiy S.V. 1 Matsulev A.N. 1
1 Institute of Chemistry and Chemical Technology SB RAS
The method for calculating the effective charge of a spherical colloidal particles using counterions distributions obtained from Langevin dynamics simulation was proposed. Algorithm is based on effective charge and Debye screening length fitting in order to minimize the root-mean-square deviation of calculated counterion´s potential energy with respect to screened Coulomb potential, obtained from the linear Debye – Huckel theory. Fitting was done for some part of counterions cloud are located far from the centre of cell, if the number of such counterions is equal to effective charge calculated by fitting. The comparison of the proposed algorithm with other algorithms for effective charge calculating has been done. A good correlation between calculated values and experimental data from literature has been shown.
counterion condensation
affective charge
primitive model
molecular dynamics
Langevin dynamics

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

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

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

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

Классическая теория Дерягина – Ландау – Вервея – Овербека (ДЛВО) [1,23], разработанная еще в 40-х годах прошлого века, предлагает простейший вариант подходящего потенциала макроиона, электростатическая часть которого имеет вид потенциала Юкавы с коррекцией на конечный радиус макроиона:

, (1)

где kB – постоянная Больцмана, T – температура среды, Z – зарядовое число макроиона, a – радиус макроиона, D – дебаевский радиус экранирования и λB – длина Бьеррума:

, (2)

e – элементарный заряд, ε0 – диэлектрическая постоянная, ε – относительная диэлектрическая проницаемость среды.

Теория ДЛВО позволяет свести анализ поведения коллоидной системы к задаче многих тел с парным потенциалом достаточно простого вида, но, к сожалению, эта теория разработана при определенных, достаточно сильных допущениях, которые не всегда легко выполнить. В частности, это условие малости кулоновского взаимодействия по сравнению с тепловой энергией, которое позволяет линеаризовать задачу, свести ее к линейному уравнению Дебая – Хюккеля и получить аналитические решения типа (1).

Тем не менее уже в первых работах [5], связанных с более общими моделями, чем ДЛВО, было обнаружено, что при определенных условиях существенная часть противоионов концентрируется в непосредственной близости от макроиона (в данном случаей линейной молекулы полииона). Такие ионы затем стали называться «сконденсированными», а сама концепция конденсации противоионов для линейных полиэлектролитов получила систематическое развитие в работах Мэннинга [17] (см. обзор [2] и ссылки в нем). Для сферических частиц, в отличие от линейных полиэлектролитов, конденсация противоионов выглядит несколько сложнее. В работе [24] было показано, что при бесконечном разбавлении в обессоленном растворителе конденсации не происходит – по мере разбавления коллоида количество связанных противоионов стремится к нулю. Но при конечных разбавлениях конденсация противоионов на сферических макроионах тоже имеет место [4,11]. Более того, не так давно в работе [18] было показано, что конденсация на сферах наблюдается и при бесконечном разбавлении, но при ненулевой концентрации соли в растворе.

Привлекательной стороной концепции конденсации противоионов является то, что макроион и сконденсироваванные на нем противоионы можно рассматривать как единую частицу, обладающую неким «эффективным» зарядом Zeff, меньшим, чем собственный заряд макроиона Z. Причем, при корректном выборе эффективного заряда появляется возможность использовать потенциал парного взаимодействия вида (1) даже вне рамок применимости линейной теории ДЛВО. Т.е. решения, получаемые по более общим моделям, вполне можно аппроксимировать решением (1) линейной теории. Нужно только в этой формуле заменить собственный заряд Z и дебаевский радиус экранирования D на их эффективные значения Zeff и D.

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

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

Другая группа методов вычисления Zeff основана на совсем другом принципе. Идея заключается в разбиении множества противоионов на два класса: сконденсированные (связанные) и несконденсированные (свободные). В этом случае эффективный заряд макроиона равен суммарному заряду свободных противоионов. Само же разбиение ведется по какому-нибудь критерию, связанному либо с расположением противоионов в пространстве, либо с их энергией. Например, в работах [7,10] используется критерий, основанный на функции распределения противоионов N(r), представляющей собой количество противоионов, находящихся внутри сферы радиуса r, центр которой совпадает с центром макроиона. Радиус, отделяющий связанные противоионы от свободных, определяется по точке перегиба функции N(r) или N(1/r). В работе [20] тоже используется критерий, связанный с концентрацией противоионов, только граница между связанными и свободными противоионами находится на радиусе, где локальная концентрация макроионов равна их средней концентрации в ячейке. В работах [8,9,21] используются «энергетические» критерии деления противоионов на свободные и связанные. Противоион считается свободным, если его потенциальная энергия в самосогласованном поле макроиона выше некоторого уровня, значение которого пропорционально кинетической энергии теплового движения противоиона. Т.е. для свободного противоиона U(r) > - χkBT. Эффективный заряд макроиона равен суммарному заряду всех свободных противоионов в ячейке Вигнера –Зейтца. В работе [8] для моделирования распределения противоионов использовался метод Монте-Карло и было показано, что при χ=2 наблюдается хорошее согласие с асимптотическим методом, предложенным в [4]. В то же время в работе [21] для моделирования использовался метод броуновской динамики, и в потенциальной энергии каждого противоиона уже учтено влияние других противоионов. Было установлено, что в этом случае корректным критерием является χ=1.

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

Модель

При моделировании коллоидных систем широкое распространение получила так называемая «ячеистая» модель, в которой с каждым макроионом ассоциирован некий объем, состоящий из элементов среды, находящихся к данному макроиону ближе, чем к остальным. Эти объемы представляют собой многогранники, называющиеся, в зависимости от предметной области, многогранниками Вороного, ячейками Дирихле или ячейками Вигнера – Зейтца. В коллоидной тематике обычно используется последнее название. Если макроионы распределены в объеме достаточно равномерно, то ячейки Вигнера – Зейтца хорошо аппроксимируются сферами.

Пусть R – радиус ячейки, a – радиус макроиона и Z – его зарядовое число. Если ограничиться сильно разбавленными коллоидами (a/R << 1), то за счет эффективного экранирования периферийные противоионы можно считать относительно свободными. Это дает возможность рассматривать каждую ячейку независимо и свести задачу к моделированию одной сферической ячейки Вигнера – Зейтца, содержащей один макроион в центре ячейки и соответствующее количество противоионов в ее объеме. Будем рассматривать однозарядные противоионы, поэтому их количество в ячейке будет равно Z. Поскольку масса макроиона существенно больше массы противоиона, то смещением макроиона за время установления равновесия в облаке противоионов можно пренебречь. Т.е. макроион, находящийся в центре ячейки, считается неподвижным. Ограничимся достаточно крупными макроионами, такими, чтобы радиусом (и, естественно, объемом) противоиона можно было бы пренебречь по сравнению с радиусом макроиона. Т.е. противоионы считаются точечными. Поверхность макроиона и границы ячейки считаются непроницаемыми для противоионов. Поэтому вблизи их границ для противоионов вводится дополнительный короткодействующий потенциал, аналогичный отталкивательной части потенциала Леннарда – Джонса для стенки и точечной частицы (1/r9).

При обезразмеривании в качестве единицы длины выбрана длина Бьеррума для чистой воды при комнатной температуре, за единицу массы взята масса противоиона, в качестве единицы энергии – kBT0, где T0 = 300K. Таким образом, основные единицы в данной задаче следующие: длина x0 = 6.88·10-10 м, время t0 = 1.9·10-12 сек, масса m0 = 3.16·10-26 кг и T0 = 300 K.

Движение противоионов в ячейке моделировалось в рамках динамики Ланжевена. В качестве расчетной схемы была взята схема с аббревиатурой «BAOAB», предложенная и исследованная в работах [15,16]. В этой схеме каждый временной шаг расщепляется на 5 шагов, имеющих следующий вид (для каждого противоиона):

(3)

где n – номер шага по времени, dt – его величина, x, ν, a – координаты, компоненты скорости и ускорения противоиона, γ – коэффициент трения, T – температура, Rn – вектор случайных чисел со стандартным нормальным распределением. После выполнения четвертого шага вычислялись силы и ускорения каждого противоиона с учетом всех остальных противоионов, макроиона, границ ячейки и макроиона. При нулевых γ и Rn эта схема эквивалентна широко известной скоростной форме алгоритма Верле.

Практическая реализация этой схемы была выполнена с использованием технологии CUDA для графических процессоров фирмы NVIDIA [3] с распараллеливанием по частицам (противоионам). Для генерации векторов Rn использовался параллельный генератор случайных чисел со стандартным нормальным распределением из библиотеки cuRAND [19]. Визуализация частиц производилась при помощи приложения AtomEye 3 [6].

Алгоритм расчета эффективного заряда

В процессе моделирования ансамбля противоионов известны положение и потенциальная энергия (в том числе и электростатическая составляющая) каждого из них. При комнатной температуре облако противоионов распределено (в общем случае неоднородно) по всему объему ячейки. Поэтому, интерпретируя противоионы как «пробные частицы», мы получаем значения самосогласованного экранированного потенциала макроиона U(r) в дискретном множестве точек пространства ячейки ri, где ri – радиус-вектор i-го противоиона (предполагается, что начало координат находится в центре макроиона). Далее, учитывая сферическую симметрию задачи, мы можем построить и одномерную зависимость потенциала макроиона Ui = U(ri) на дискретном нерегулярном множестве точек ri, являющихся расстоянием от центра ячейки до i-го противоиона. Теперь можно провести сортировку множества пар (ri,Ui) в порядке убывания ri и решить следующую задачу нелинейной оптимизации:

, (4)

где U(r,Z,D) – потенциал Юкавы, задаваемый выражением (1), Zv и Dv – варьируемые вещественные параметры, Nfit – целочисленный параметр, задающий количество периферийных точек (противоионов), участвующих в фиттинге. При уменьшении Nfit, начиная с Nfit = Z, функция

(5)

в общем случае несколько раз пересекает ноль. Для получения оптимального значения эффективного заряда Zeff берется самый первый корень, т.к. при таком выборе количество противоионов, участвующих в фиттинге, максимально. Точнее, в качестве оценки Zeff выбирается целочисленное значение Nfit, ближайшее к точке первого пересечения нуля функцией δ(Nfit). После того, как Zeff определен, производится уточнение эффективного дебаевского радиуса экранирования Deff путем минимизации (6) по Deff уже при фиксированном Zeff:

. (6)

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

Во всех процедурах минимизации использовались функции нелинейной оптимизации из библиотеки GNU Scientific Library (GSL) [12] (алгоритм Левенберга – Марквардта).

Результаты

В описанной выше постановке задачи есть четыре независимых параметра: температура среды T, собственный заряд частицы Z, радиус макроиона a и радиус ячейки R. Температура среды во всех расчетах была зафиксирована и равна комнатной. Можно зафиксировать еще один параметр, если учесть, что при сильных разбавлениях коллоида (φ = (a/R)3 <<1) эффективный заряд Zeff довольно слабо зависит от радиуса ячейки R. Это легко видеть на рис.1, где для иллюстрации приведена одна из таких зависимостей. Поэтому для целей данной работы вполне можно ограничиться каким-нибудь одним, достаточно большим значением радиуса ячейки (таким, чтобы выполнялось условие φ ≤ 0.01).

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

Рис.1. Зависимость эффективного заряда от радиуса ячейки. Z = 2000, a = 20

Первый тест заключается в обнаружении самого факта начала конденсации противоионов: при уменьшении радиуса макроиона a меньше некоторого критического значения (при этом собственный заряд макроиона Z зафиксирован), Zeff становится меньше Z, т.е. начинается конденсация противоионов. Поскольку молекулярная динамика дает максимально полную (в рамках выбранной модели) информацию о системе, то можно рассчитать значения эффективного заряда Zeff по нескольким алгоритмам на одних и тех же распределениях противоионов и сравнить результаты в рамках этого же теста. Они представлены на рис.2.

Рис. 2. Зависимости эффективного заряда от радиуса макроиона, рассчитанные по разным алгоритмам. Кривая 1 – данная работа, 2 – критерий U>-2kBT [8], 3 – по точке перегиба функции распределения противоионов [7,10], 4 – критерий U>-kBT [21].

Z = 500, R = 500

Из рисунка видно, что при относительно малых радиусах (при заданных данных Z и R), когда конденсация противоионов происходит достаточно эффективно, различие между значениями эффективного заряда, рассчитанными по различным критериям, невелика. Различие становится заметным при больших радиусах, когда граница между сконденсированными и свободными противоионами становится не совсем очевидной. При радиусах макроиона a>80 по алгоритму, предлагаемому в данной работе, конденсация противоионов уже отсутствует (т.е. все макроионы свободны: Zeff = Z). Из графиков, рассчитанных по другим критериям, видно, что от 10 до 20 % противоионов еще остаются сконденсированными. Тем не менее, если посмотреть на реальное распределения противоионов в ячейке, представленное на стереопаре (рис. 3), то видно, что никакого заметного концентрирования противоионов у поверхности макроина (в данном случае при a=100) не наблюдается. Всего лишь несколько противоионов, находящихся непосредственно на поверхности макроиона, выделены цветом.

Таким образом, предлагаемый в данной работе критерий, с точки зрения деления противоионов на свободные и сконденсированные, интуитивно представляется более корректным. Для просмотра стереопары глаза должны быть перекрещены.

Рис. 3. Пример реального распределения противоионов в пространстве ячейки

Итак, предложенный в данной работе алгоритм вполне корректно «фиксирует» сам факт конденсации противоионов. Следующий тест связан с так называемым явлением «насыщения» эффективного заряда. Суть явления в том, что при увеличении собственного заряда макроиона Z (при фиксированном радиусе макроиона a) эффективный заряд Zeff сначала растет линейно. Затем (при дальнейшем увеличении Z) начинается конденсация противоионов и рост Zeff замедляется. В конце концов рост останавливается на некотором уровне Zsat и далее уже не меняется. Т.е., образно говоря, в режиме насыщения все «лишние» противоионы конденсируются на макроионе, а количество свободных противоионов уже не зависит от Z. На рис. 4 слева представлены диаграммы, иллюстрирующие этот эффект.

Рис. 4. Явление насыщения эффективного заряда для нескольких значений радиуса макроиона (цифры у кривых на диаграммах слева). Справа – зависимость насыщенного значения эффективного заряда от радиуса макроиона

Из правого графика на рис. 4 видно, что насыщенное значение эффективного заряда Zsat зависит от радиуса частицы линейно. Провести непосредственное количественное сопоставление полученных значений Zsat c данными, имеющимися в литературе [4,22], затруднительно из-за сильно различающихся концентраций коллоида, но качественное поведение кривых совпадает.

Последний тест заключается в сравнении рассчитанных значений эффективного заряда с экспериментальными данными, имеющимися в литературе. В таблице 3 работы [20] приведена подборка экспериментально измеренных значений собственного и эффективного зарядов сферических частиц из нескольких литературных источников. Большая часть приведенных результаты относится к супензиям латексов. Эти экспериментальные данные были преобразованы к используемым здесь безразмерным единицам и представлены на рис.5 вместе со значениями Zeff, рассчитанными в данной работе. Эффективный заряд вычислялся по распределениям противоионов, полученных в диапазоне изменения параметров задачи: Z – от 100 до 2000 и a – от 10 до 100. Радиус ячейки во всех вариантах был зафиксирован R=500. На графиках изображена зависимость относительного значения эффективного заряда Zeff/Z от параметра Z/a. Видна неплохая корреляция рассчитанных и экспериментальных данных.

Рис. 5. Сравнение рассчитанных и экспериментальных значений эффективного заряда. Точки – данная работа (расчет), крестики – таблица 3 из [20] (эксперимент)

В заключение отметим, что рассчитанные в данной работе значения эффективного заряда в диапазоне изменения параметров Z – 100...2000 и a – 100...100 хорошо описываются формулой:

. (7)

Соответствующая кривая изображена на рис. 5 пунктиром.

Обсуждение и выводы

В данной работе предложен и опробован алгоритм расчета эффективного заряда коллоидной частицы Zeff и эффективного дебаевского радиуса экранирования Deff по распределениям противоионов, получаемым в результате моделирования методом динамики Ланжевена. Значения Zeff и Deff получаются в результате фитинга, в котором участвует только часть противоионов, наиболее удаленных от макроиона. На количество противоионов, участвующих в фиттинге, накладывается условие его равенства эффективному заряду макроиона, которое получается в результате фиттинга. Это условие представляется вполне естественным, т.к. именно эти противоионы считаются «свободными» и экранируют поле макроиона приблизительно так, как это происходило бы в рамках линейной задачи.

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

В заключение можно сформулировать следующие выводы:

  1. Предложен однозначный алгоритм расчета эффективного заряда и дебаевского радиуса экранирования по распределениям противоионов в ячейке, получаемым методами динамики Ланжевена.
  2. Показана работоспособность алгоритма в широком диапазоне изменения параметров задачи: Z, R и a.
  3. Продемонстрирована хорошая корреляция рассчитанных значений эффективного заряда с существующими в литературе экспериментальными результатами.

Рецензенты:

Тарабанько В.Е., д.х.н., профессор, заведующий лабораторией Комплексной переработки биомассы Института химии и химической технологии СО РАН, г. Красноярск.

Белобров П.И., д.ф.-м.н., с.н.с., профессор кафедры биофизики СФУ, в.н.с. Института биофизики СО РАН, г. Красноярск.