1.Введение

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

Исходная информация для расчёта равновесного состава раствора — концентрации компонентов раствора и термодинамические константы равновесий реакций взаимодействия между частицами в растворе. Приложение решает так называемую «прямую задачу» химических равновесий: по заданным начальным концентрациям компонентов смешанного электролитного раствора и заданным значениям термодинамических констант равновесий определяет полный набор равновесных концентраций всех частиц. Точность расчёта определяется только адекватностью уравнения, используемого для описания коэффициентов активности ионов. Постановка задачи, выбор базиса и типичные затруднения машинного счёта подробно разобраны в^[4], а проверка того, заслуживает ли доверия построенная модель равновесий, — в^[5].

2.Теоретическое обоснование

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

В настоящем приложении использован универсальный подход, основанный на формализме Бринкли — описании системы равновесий в матричной форме, что позволяет унифицировать схему расчёта независимо от вида и сложности задачи. Матричный подход впервые предложен для решения задачи химических равновесий в газовой фазе ^[1], рекомендован и практически использован для решения обратной задачи расчёта химических равновесий в растворах электролитов ^[2][3].

Базис системы

В равновесной системе, состоящей из $m$ частиц $\mathbf{A}$, выбирается минимальный набор из $q$ частиц $\mathbf{B}$, которые определяют базис системы. При известном механизме взаимодействия между частицами базис определяет общий состав всей системы. Как правило, в качестве компонентов базиса выбирают из всей совокупности частиц $\mathbf{A}$ простые частицы, из которых путём взаимодействия получаются более сложные.

Математически базис должен удовлетворять условиям полноты (все ионные и молекулярные формы сложной равновесной системы могут быть образованы по реакциям между компонентами базиса) и независимости (ни один из компонентов базиса не может быть образован из других). Химические реакции в таком случае представляются в виде линейной комбинации компонентов базиса:

где $\nu_{ij}$ ($m\times q$) — матрица стехиометрических коэффициентов реакции, $B_j$ — частицы базиса, $A_i$ — частицы системы.

Для описания системы равновесий выбор базиса не имеет значения и может быть любым, но для практической реализации расчётов, как правило, базис определяется набором частиц из компонентов раствора — так как расчёт основан на согласовании известных начальных концентраций частиц базиса и рассчитанных через них равновесных концентраций.

Каждая сложная форма $A_i$ представляется базисом и вектором, который определяет частицу в единицах базиса. Химические реакции рассматриваются как реакции образования частиц. Если в растворе не происходит образования новых частиц, то $q = m$.

Для единого описания всех частиц в растворе — в том числе частиц, входящих в сам базис ($\mathbf{B} \subset \mathbf{A}$) — реакции для них превращаются в тождества: если $i=j$ ($B_j \equiv A_i$), то $\nu_{ij}=1$, при $i\neq j$ — $\nu_{ij}=0$, то есть $\nu_{ij} = \delta_{ij}$, где $\delta_{ij}$ — символ Кронекера. Такие формальные реакции ($A_i \equiv B_j$) удобно располагать в верхней части стехиометрической матрицы — тогда верхняя квадратная часть её представлена единичной матрицей.

3.Матричный формализм

В качестве сквозного примера ниже приведено представление равновесий в смешанном растворе с тремя компонентами: двухосновной кислоты $H_2A$, соли этой кислоты $M_2A$ и слабой одноосновной кислоты $HX$. В растворе могут присутствовать частицы ($\mathbf{A}$): $H^+,\ A^{2-},\ M^+,\ X^-,\ HA^-,\ H_2A,\ HX,\ OH^-$. Система равновесий представляется стехиометрической матрицей при базисе ($\mathbf{B}$): $H^+,\ A^{2-},\ M^+,\ X^-$.

Первые четыре строки — формальные реакции базис↔базис (единичная подматрица). Отрицательный коэффициент у $OH^-$ отражает формальную реакцию $OH^- \equiv {H^+}^{-1}$ (см. ниже).

Закон действующих масс

Выражение закона действующих масс для любой $i$-й реакции имеет вид:

где $c_i$ — равновесная концентрация $i$-й частицы, $\gamma_i$ — её коэффициент активности, а $a_j$ — активность $j$-й частицы базиса, взятая в степени $\nu_{ij}$, соответствующей стехиометрическому соотношению в $i$-й реакции.

Постоянная активность: осадки и растворитель

Если концентрация (активность) некоторых частиц может быть принята постоянной (или её изменениями можно пренебречь), то такую активность включают в саму константу равновесия. Например, в гетерогенных реакциях активность твёрдой фазы считается постоянной (равной единице для чистого вещества в стандартном состоянии) и включается в произведение растворимости $K_{sp}$. Для малорастворимой соли $M_2A$ в равновесии с раствором:

Это приводит к формальной реакции образования иона $M^+$ из частицы базиса $A^{2-}$:

Это соответствует отдельной строке в стехиометрической матрице, определяющей равновесную концентрацию ионов $M^+$ малорастворимой соли: коэффициент при базисном анионе — $-\tfrac12$, константа равновесия в такой записи — $K_{sp}^{1/2}$.

Аналогичная ситуация — в случае расчёта концентрации гидроксильных ионов в водном растворе:

При постоянстве активности молекулярной формы воды (в случае её избытка) образование гидроксильных ионов представляется формальным уравнением: в строке стехиометрической матрицы в столбце иона водорода стоит элемент $-1$, а константа равновесия соответствует ионному произведению воды $K_w$.

Уравнения материального баланса

Зная равновесные концентрации всех частиц в растворе, можно составить для каждой из $q$ частиц базиса уравнение материального баланса. Для некоторой $l$-й частицы базиса с аналитической концентрацией $C_l$ в растворе уравнение материального баланса имеет вид:

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

Таким образом, всю систему уравнений можно записать в виде двух матричных уравнений в канонической форме:

где $\mathbf{c}$ — вектор равновесных концентраций всех $m$ частиц, $\mathbf{a}$ — вектор активностей $q$ частиц базиса, $\mathbf{C}$ — вектор их аналитических концентраций, $\mathbf{\nu}^{\mathsf T}$ — транспонированная стехиометрическая матрица.

Компонентная и концентрационная матрицы

Рис. 1. Шаг 2 мастера ввода — компонентная матрица. Два компонента (навески, которые действительно взвешивают) и базис из четырёх частиц; строка матрицы есть разложение навески по базису: Na₂HPO₄·2H₂O = 1·H⁺ + 2·Na⁺ + 1·PO₄³⁻, C₆H₈O₇·H₂O = 3·H⁺ + 1·Cit³⁻ (кристаллизационная вода в разложение не входит — она растворитель). Отсюда программа сама получает аналитические концентрации базиса для каждого раствора. Компоненты названы формулами, а не словами не для краткости: ровно эта строка стоит потом на оси абсцисс графика (рис. 5).

Рис. 2. Шаг 3 той же задачи — стехиометрическая матрица, вектор зарядов и вектор суммарных lg K. Первые четыре строки — частицы базиса: формально «базис через базис», поэтому там единичная подматрица и нулевые константы, и эти ячейки серые — их не правят. Ниже идут производные частицы, каждая — строка коэффициентов при частицах базиса: HPO₄²⁻ = H⁺+PO₄³⁻ (lg K 12,35), H₂PO₄⁻ = 2H⁺+PO₄³⁻ (19,55), H₃PO₄ (21,70), три ступени цитрата (6,40 / 11,16 / 14,29) и автопротолиз воды OH⁻ = −H⁺ (lg K −14,17). Жёлтые столбцы — заряд и константа; матрицу можно вставить готовым блоком из буфера обмена. Имена этих одиннадцати строк — это заголовки столбцов таблицы результатов (рис. 4) и имена кривых на графике (рис. 5).

Начальные (аналитические) концентрации частиц базиса рассчитываются на основе компонентной матрицы $\omega$ ($p\times q$, $p$ — количество компонентов раствора) и массива концентраций компонентов раствора $\mathbf{C}_{\text{комп}}$ ($n\times p$, $n$ — количество растворов серии):

Для приведённого выше примера ($H_2A$, $M_2A$, $HX$ — компоненты; $H^+,A^{2-},M^+,X^-$ — базис) компонентная матрица $\omega_{3\times4}$:

И, например, для массива из четырёх растворов со следующими концентрациями компонентов, моль/л ($\mathbf{C}_{\text{комп},\,4\times3}$):

матрица начальных (аналитических) концентраций частиц базиса $\mathbf{C}_{\text{базис}} = \mathbf{C}_{\text{комп}}\,\omega$ имеет вид:

Именно этот пересчёт (компонентная матрица → концентрации базиса) выполняется автоматически на шаге 4 мастера ввода данных — если не включён режим прямого ввода концентраций базиса, минуя компонентную матрицу.

4.Алгоритм решения

Коэффициенты активности

Поскольку параметрами равновесной системы выступают термодинамические константы равновесий реакций, необходимо использовать активности частиц базиса. Для незаряженных частиц с высокой степенью достоверности коэффициенты активности можно принять равными единице. Для описания коэффициентов активности ионов в достаточно широкой области концентраций (в водных растворах — вплоть до высоких концентраций, в отдельных случаях до 1 моль/л) используется полуэмпирическое уравнение на основе теории Дебая–Хюккеля (в модификации Дэвиса):

где $A_D$ и $B_D$ — теоретические параметры уравнения Дебая–Хюккеля, зависящие от температуры и диэлектрической проницаемости растворителя; $a^0$ — полуэмпирический параметр «максимального сближения ионов»; $b$ — эмпирический коэффициент (параметр Дэвиса); $I$ — ионная сила раствора; $z_i$ — заряд иона.

Сами параметры $A_D$ и $B_D$ вычисляются по температуре $T$, диэлектрической проницаемости $\varepsilon$ и плотности $\rho$ растворителя:

где $e$ — заряд электрона, $k_B$ — постоянная Больцмана, $N_A$ — число Авогадро (константы в системе СГСЭ). Эти значения рассчитываются программой автоматически и отображаются на шаге 1 мастера ввода данных сразу при изменении параметров растворителя.

Постановка прямой задачи

Целью прямой задачи химических равновесий является — для заданного вектора начальных (аналитических) концентраций частиц базиса $\mathbf{C} = (C_1,\dots,C_q)$ — определить вектор равновесных концентраций всех частиц $\mathbf{c} = (c_1,\dots,c_m)$, используя закон действующих масс и уравнения материального баланса при заданной стехиометрической матрице и векторе констант равновесия.

Метод Ньютона по логарифмам активностей базиса

Система из $q$ нелинейных уравнений материального баланса с $q$ неизвестными — в качестве которых выступают логарифмы активностей базисных частиц $\ln a_j$ — решается методом последовательных приближений (методом Ньютона). Вектор невязок материального баланса:

Якобиан этой системы (в силу $\partial c_i/\partial\ln a_k = \nu_{ik}\,c_i$):

На каждой итерации решается линейная система и уточняется вектор $\ln\mathbf a$:

где $N$ — номер итерации, $\Delta\ln\mathbf a$ — ньютоновская поправка, $\lambda$ — шаг (в текущей реализации применяется демпфирующий коэффициент $\lambda=0{.}5$ для устойчивости сходимости). Линейная система $\mathbf J\,\Delta\ln\mathbf a=-\mathbf g$ решается методом Гаусса с выбором ведущего элемента по столбцу.

Начальное приближение и сходимость

В качестве начального приближения принимается $\ln a_j^{(0)} = \ln C_j$ — то есть логарифм аналитической (начальной) концентрации соответствующей частицы базиса. Первые итерации проводятся при условии равенства единице всех коэффициентов активности заряженных частиц; далее, по мере уточнения равновесного состава, на каждом цикле пересчитывается ионная сила, и коэффициенты активности ионов всё точнее приближаются к истинным значениям.

Итерационный процесс завершается по достижении относительной погрешности уточнения логарифмов активности частиц базиса не более $10^{-7}$ либо по достижении предельного числа итераций. Метод Бринкли не имеет принципиальных ограничений на количество и стехиометрию равновесных реакций и может быть успешно использован для решения самых сложных задач химического равновесия.

5.Коэффициенты активности

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

5.1. Зачем нужна модель и когда её выбор безразличен

Ион-ионное взаимодействие входит в расчёт равновесия косвенно, но входит всегда. Закон действующих масс записан через активности $a_i = c_i\gamma_i$, ионная сила

зависит от найденного состава, а состав — от коэффициентов активности. Отсюда и берётся внешний итерационный цикл, описанный в разделе 4: сначала система решается при $\gamma_i \equiv 1$, затем по полученному составу считается $I$, по ней — коэффициенты активности, и всё повторяется до самосогласования.

Вопрос «насколько ответ обязан модели» рассуждением не решается. Влияние косвенное, а его величина зависит от конкретной задачи: от зарядов частиц, от того, есть ли ассоциация, и главное — от достигнутой ионной силы. Единственный способ узнать — посчитать одно и то же разными моделями. Именно для этого в программе есть пункт «Сравнение моделей» (см. 5.8), а в потенциометрии и кондуктометрии — вкладка сравнения по качеству подгонки (см. 5.9).

Порядок величин виден на простом примере: раствор $\mathrm{MgCl_2}$ в воде при 25 °C ($A = 0.5106$, $B = 0.3291\ \mathring{\mathrm{A}}^{-1}$), параметры расширенного уравнения $a_0 = 4.0$ Å, $b = 0$, радиусы для MSA $r(\mathrm{Mg^{2+}}) = 3.47$ Å, $r(\mathrm{Cl^-}) = 1.81$ Å. В таблице — коэффициент активности двухзарядного катиона, посчитанный всеми четырьмя моделями при двух ионных силах.

Читать эту таблицу нужно так. При $I \approx 0.01$ три «разумные» модели согласны в пределах трёх процентов, и выбор между ними безразличен: разница заведомо меньше погрешности констант равновесия. При $I = 0.30$ те же три расходятся уже на 38 %, а предельный закон отличается от них вчетверо — то есть перестаёт быть приближением вовсе. Именно поэтому модель стала выбором, а не была оставлена вшитой: задача комплекса — поднять порог применимости расчёта хотя бы до $I = 0.3$, а на этом уровне ответ уже заметно зависит от того, чем считать.

Чем модель отзывается в самом счёте

Внешний цикл устроен по-разному для разных моделей, и это не деталь реализации, а следствие их устройства. У всего семейства Дебая–Хюккеля коэффициент активности зависит от состава только через ионную силу: $\gamma_i = f(z_i, I)$. Поэтому сходимость по $I$ равносильна сходимости по $\gamma$, и критерий остановки — относительное изменение ионной силы.

У MSA это не так: в неё входят числовые плотности каждого сорта частиц по отдельности. Ионная сила может сойтись при ещё не сошедшихся коэффициентах активности — например, когда перераспределение между двумя частицами одинакового заряда не меняет $I$, но меняет состав твердосферной смеси. Поэтому для таких моделей к критерию добавляется норма изменения $\lg\gamma$. Программа определяет это по самой модели, а не по списку имён, так что новая модель с зависимостью от состава получит правильный критерий автоматически.

Практическое следствие: расчёт с MSA идёт заметно дольше, чем с уравнениями семейства Дебая–Хюккеля, — и это цена ион-специфичности, а не признак неполадки.

5.2. Лестница из четырёх моделей

Модель выбирается на шаге 1 мастера ввода данных. Модели расположены по возрастанию сложности; у первых двух подгоняемых параметров нет вовсе, у третьей их два, у четвёртой — один, и он же имеет физический смысл размера иона.

Границы указаны для водных растворов 1:1 и заметно сужаются с ростом зарядов; в неводной среде они смещаются ещё сильнее (см. 5.6). Двухбуквенные коды используются в заголовках колонок сводок сравнения — то же соглашение, что в разделе «Математическая обработка».

5.3. Термодинамическая согласованность

Набор индивидуальных коэффициентов активности не может быть произвольным. Он должен порождаться избыточной энергией Гиббса $G^E$ — одной функцией состава, из которой все $\gamma_i$ получаются дифференцированием:

Отсюда немедленно следует требование симметрии смешанных производных (условие Гиббса–Дюгема в дифференциальной форме):

Проверим на нашем семействе. Пусть $\ln\gamma_i = z_i^2\,g(I)$, то есть зависимость от сорта частицы входит только множителем $z_i^2$. Так как $\partial I/\partial c_j = z_j^2/2$, получаем

— выражение симметрично относительно перестановки $i \leftrightarrow j$, и всё в порядке. Более того, порождающая функция выписывается явно. Условие $\ln\gamma_i = z_i^2 g(I)$ не только достаточно, но и необходимо в классе моделей, зависящих от состава лишь через ионную силу: симметрия требует $f_i'(I)z_j^2 = f_j'(I)z_i^2$ для всех пар, то есть $f_i'(I)/z_i^2$ одинаково для всех сортов частиц.

Предельный закон и уравнение Дэвиса этому условию удовлетворяют: у обоих весь правый член, включая линейный, умножен на $z_i^2$. А у расширенного уравнения оно нарушено дважды:

• линейный член $b\,I$ не масштабирован зарядом — он одинаков у однозарядного и у трёхзарядного иона. Тогда $\partial\ln\gamma_i/\partial c_j \ni \ln\!10\cdot b\, z_j^2/2$, а симметричная производная даёт $\ln\!10\cdot b\, z_i^2/2$; они равны только при $z_i^2 = z_j^2$;
• разные $a_0$ у разных ионов делают $g$ зависящей от сорта частицы сверх множителя $z_i^2$, и та же проверка снова не проходит.

Всё это закодировано в программе. Осмотический коэффициент $\varphi$ и активность растворителя возвращаются как «не число», и в таблице сравнения на их месте стоит прочерк — по конвенции комплекса прочерк значит «величины не существует», а не «не посчитали». В колонке «Согласована» такая строка помечается «нет».

Замкнутые формулы $\varphi$ для согласованных моделей сверены с численным интегрированием по Гиббсу–Дюгему вдоль луча составов (замена $\lambda = u^2$ снимает корневую особенность подынтегральной функции в нуле). У MSA этот же тест проверяет согласованность сразу в четырёх местах — электростатику и твёрдые сферы, $\ln\gamma$ и избыточную энергию.

Почему на расчёт равновесия это не влияет

Рассуждение простое и стоит его привести, потому что иначе предыдущий абзац выглядит приговором расширенному уравнению. Пусть все коэффициенты активности сдвинуты произвольным множителем вида $\gamma_i \to \gamma_i\lambda^{z_i}$. В законе действующих масс для реакции образования $A_i = \sum_j \nu_{ij}B_j$ стоит комбинация

Но реакция сбалансирована по заряду: $z_i = \sum_j \nu_{ij} z_j$, поэтому показатель равен нулю и множитель уходит целиком. Равновесный состав к сдвигам этого класса нечувствителен — и именно они составляют произвол в определении индивидуальных коэффициентов активности.

Что затрагивается: всё, что выводится из растворителя (осмотический коэффициент, активность воды, понижение температуры замерзания), и всякий перенос модели в термодинамику — расчёт растворимости через $\Delta G$, согласование с калориметрией. А также — сравнение ваших параметров с чужими: несогласованный набор нельзя отождествлять с табличным.

5.4. Что отвергнуто и почему

Читатель, знакомый с литературой, немедленно спросит, почему в списке нет уравнений Питцера. Ответы на такие вопросы — не история разработки, а руководство к применению: каждое из перечисленных ниже уравнений хорошо в своей области, и понимание, почему оно не годится здесь, помогает понять, что именно делает комплекс.

Уравнения Питцера — нет

Модель Питцера ^[10] — общепринятый стандарт для концентрированных растворов, и в своей области ей нет равных. Но её вириальные коэффициенты $\beta^{(0)}, \beta^{(1)}, C^\varphi$ подобраны так, что уже вобрали в себя ионную ассоциацию: специфическое взаимодействие пары ионов описано эмпирическими членами вместо явного комплекса.

Комплекс же строит равновесие с явными комплексами — они вводятся как обычные частицы со своими $\lg K$ (в этом весь смысл формализма Бринкли). Соединить одно с другим значит учесть ассоциацию дважды: сначала явной частицей, потом вириальным членом. Результат будет выглядеть правдоподобно и будет неверен, а обнаружить это по невязкам почти невозможно — оба механизма дают качественно похожий вклад.

SIT — нет

Теория специфического ионного взаимодействия ^[11] проще Питцера и широко используется в ядерной химии. Её база параметров взаимодействия — водная (NEA-TDB и производные). Комплекс же рассчитан на любые среды: расчёт в ацетонитриле или в смешанном растворителе — штатная задача, а не экзотика.

В неводной среде табличный коэффициент взаимодействия взять неоткуда, и SIT вырождается в «попарный подгоночный $b$», то есть даёт ровно то же, что ион-специфичный Хюккель, уже имеющийся в списке, — но при этом создаёт иллюзию табличной обоснованности.

Робинсон–Стокс (гидратационная поправка) — нет

Модель гидратации ^[12] вводит число молекул растворителя, связанных ионом, и поправку на уменьшение количества свободного растворителя. Два возражения. Во-первых, это специфическая сольватация: число гидратации — параметр той же природы, что и подгоночный $b$, и для неводных сред не табулирован. Во-вторых и главное, модель выведена для среднего $\gamma_\pm$ одной соли и на многокомпонентную смесь не переносится без дополнительных допущений о том, как ионы делят между собой связанный растворитель.

Ла Мер–Гронуолл–Грайф — нет (факультативно)

Строгое решение нелинейного уравнения Пуассона–Больцмана в виде ряда ^[13]. Математически безупречно и исторически важно, но исправляет линеаризацию той же примитивной модели — точечных зарядов в сплошном диэлектрике. MSA делает то же самое лучше (учитывая ещё и объём ионов) и без медленно сходящихся рядов. Отвергнуто не как неверное, а как избыточное.

Константа Сеченова отдельным параметром — нет

Высаливание нейтральных частиц («эффект Сеченова», $\lg\gamma_0 = k_s I$) можно было бы ввести отдельным подгоночным коэффициентом $k_s$ на каждую нейтральную частицу. Это сознательно не сделано: такой коэффициент был бы чисто подгоночным и добавлял бы по параметру на частицу.

В MSA высаливание нейтральных получается даром, из члена твёрдых сфер: нейтральная частица конечного размера должна потесниться среди ионов, и работа этого вытеснения растёт с концентрацией. Никакого нового параметра при этом не появляется — только уже введённый радиус. См. 5.5.

5.5. MSA — что в ней происходит

Среднесферическое приближение решает уравнение Орнштейна–Цернике для смеси заряженных твёрдых сфер в электростатической среде. Коэффициент активности складывается из двух независимых вкладов:

Электростатическая часть

Определяется параметром экранирования $\Gamma$, который находится самосогласованно из уравнения, связывающего его с составом и размерами всех сортов частиц (уравнения (18)–(21) работы ^[9] — величины $\Delta$, $\Omega$, $P_n$, $N_i$). При равных диаметрах $\sigma$ уравнение решается в замкнутом виде,

где $\kappa$ — обратная дебаевская длина; при $\sigma \to 0$ отсюда $\Gamma \to \kappa/2$, и электростатическая часть переходит ровно в предельный закон. Это не совпадение, а обязательное свойство, и оно проверяется тестом.

Физический смысл $\Gamma$ стоит понимать: это обратная длина экранирования с поправкой на конечный размер ионов. В предельном законе роль экрана играет дебаевская длина $1/\kappa$, вычисленная для точечных зарядов; конечный размер мешает зарядам подойти вплотную, экран получается менее плотным, и $\Gamma < \kappa/2$. Именно поэтому MSA даёт коэффициенты активности выше предельного закона — что и видно в таблице 5.1 при обеих ионных силах.

Разница с расширенным уравнением Дебая–Хюккеля здесь принципиальная, хотя оба учитывают размер. В расширенном уравнении размер входит поправкой к готовой формуле: один параметр $a_0$ в знаменателе, отвечающий «расстоянию наибольшего сближения» — величине, которую нельзя измерить независимо. В MSA размеры входят в постановку задачи: это диаметры твёрдых сфер, из которых складывается уравнение состояния, и $\Gamma$ получается решением, а не подстановкой.

Часть твёрдых сфер

Берётся в форме Боублика–Мансури–Карнахана–Старлинга–Леланда (BMCSL) — стандартное уравнение состояния смеси твёрдых сфер ^[14]. Оно даёт химический потенциал внедрения частицы радиуса $r_i$ в среду сфер разных размеров и выражается через моменты распределения $\xi_n = \frac{\pi}{6}\sum_k \rho_k \sigma_k^n$.

Отсюда сразу два следствия, важных на практике:

• Высаливание нейтральных частиц. У незаряженной частицы электростатический вклад равен нулю, а твердосферный — нет: он и есть эффект Сеченова, полученный без единого нового параметра. У частицы с нулевым размером остаётся $-\ln(1-\xi_3)$ — работа внедрения точки в среду твёрдых сфер.
• Отказ при невозможной упаковке. Если заданные радиусы дают плотность упаковки $\xi_3 \ge 0.9$, программа явно отказывается считать и указывает на радиусы, а не выдаёт молча бессмысленное число.

Один радиус на $\gamma$ и на $\Lambda$

Радиус хранится один на всю задачу. Когда в «Кондуктометрии» показывается колонка «$r$, Å», собственная колонка MSA не добавляется — те же значения идут в модель равновесия. Это методическое приобретение, а не экономия полей: в совместном анализе параметров становится меньше, и появляется проверяемое утверждение — размер иона, найденный по электропроводности, обязан воспроизводить коэффициенты активности.

Концентрационно-зависимый диаметр (модель Симонена и подобные) сознательно не берётся: он вернул бы $\sigma$ в разряд подгоночных величин и убил бы то самое свойство, ради которого MSA и взята.

Уравнения для $\Gamma$ — те же, что в транспортной MSA раздела «Кондуктометрия»; копии самодостаточны (конвенция комплекса: каждый порт со своими вспомогательными функциями), а держит их вместе сквозной тест согласия. Остаточное расхождение по $\Gamma$ составляет $6.5\cdot 10^{-5}$ — это ровно разница в значении $e^2$ между $A$ и константами СГСЭ транспортной части.

Проверки MSA, каждая независима от реализации: $\sigma \to 0$ воспроизводит предельный закон; равные диаметры дают замкнутое $\Gamma$, выписанное выше; твердосферный член у одного компонента совпадает с химическим потенциалом Карнахана–Старлинга.

5.6. Неводные среды

Комплекс не привязан к воде, и это накладывает ограничения на выбор модели сильнее, чем принято думать. Диэлектрическая проницаемость входит в $A$ в степени $-3/2$: при переходе от воды ($\varepsilon = 78.3$) к среде с $\varepsilon = 20$ множитель $A$ вырастает примерно в 7.7 раза. Это значит, что при той же ионной силе отклонения от идеальности в 7.7 раза больше, и указанные в 5.2 границы применимости смещаются в сторону разбавления на порядки.

Отсюда три практических вывода:

• Справочные $a_0$ и $b$ не переносятся. Табличные значения расстояния наибольшего сближения подобраны по водным данным; в другой среде они становятся просто подгоночными числами. Заявлять «взято из справочника» в этом случае нельзя.
• Ассоциация становится главным эффектом, а не поправкой. Длина Бьеррума $l_B = e^2/(\varepsilon k_B T)$ растёт обратно пропорционально $\varepsilon$: в среде с $\varepsilon = 20$ она вчетверо больше водной, и ионные пары образуются при концентрациях, где в воде их ещё нет. Правильный ответ — ввести пару явной частицей со своей константой, а не пытаться спрятать её в коэффициент активности. Формализм Бринкли для этого и приспособлен.
• Модель, определённая через одни константы, безопаснее. Предельный закон и Дэвис не содержат ничего, кроме $A$, то есть кроме $\varepsilon$, $T$ и $\rho$ — величин, которые для неводной среды известны или измеримы. MSA добавляет только размеры ионов, а они от растворителя зависят слабее, чем эмпирические поправки.

5.7. Как выбор модели отражается в других приложениях

Выбор модели хранится в равновесной части проекта и потому доходит до всех приложений комплекса автоматически. Ниже — что он меняет в каждом из них.

Потенциометрия

Измеряемая ЭДС связана с активностью потенциалопределяющего иона: $E = E_0 + \sigma\,(S/n)\sum_i \nu_i \lg a_i$, $a_i = c_i\gamma_i$. Смена модели меняет $\lg\gamma$, то есть сдвигает расчётную ЭДС.

Осторожность, однако, нужна и в обратную сторону: если вместе с $E_0$ подгоняется ещё и $\lg K$, они вдвоём могут почти полностью скомпенсировать смену модели. Признак этого — высокая корреляция $E_0$ и $\lg K$ в матрице (см. 5.9); в типичной задаче с одной константой она нередко доходит до $|r| = 1.000$.

Осаждение (потенциометрическое титрование)

Условие насыщения записано через активности: $\prod_j a_j^{s_j} = K_{sp}$. Модель входит сюда напрямую, причём одинаково для точек до скачка и после — расчёт осадка идёт той же выбранной моделью, что и гомогенный расчёт.

Кондуктометрия

Здесь модель влияет двумя путями. Первый — обычный: через равновесный состав, то есть через концентрации свободных ионов, которые и переносят ток (принцип Фуосса–Джастиса: модель электропроводности вычисляется при свободных концентрациях, а не аналитических). Второй — специфический: под MSA радиус становится равновесно значимым.

Разница между этими случаями имеет прямое вычислительное следствие, и о ней стоит знать. В семействе Дебая–Хюккеля радиус на равновесие не влияет вовсе, поэтому «$\lambda^0$ и $r$» — дешёвые ручки: их варьирование не требует пересчёта равновесия, и подгонка идёт быстро. Под MSA радиус входит в состав, кэш равновесия на него реагирует, и та же подгонка становится на порядок дороже. Это не дефект, а цена ион-специфичности.

Набор подгоняемых параметров

В обоих приложениях список кандидатов на подгонку зависит от модели: $a_0$ и $b$ предлагаются, только если их читает выбранная модель. У предельного закона, Дэвиса и MSA этих параметров нет, и подгонка по ним двигала бы величину, ни на что не влияющую: матрица Гессе вышла бы вырожденной, а доверительные интервалы — недостоверными. Пользователь получил бы предупреждение вместо объяснения, поэтому строки просто не появляются.

5.8. Сравнение моделей по составу

Пункт «Сравнение моделей» на левом рельсе считает один и тот же проект несколькими моделями и раскладывает результаты рядом — тот же приём, что у прогонов сглаживания и аппроксимации в разделе «Математическая обработка». Предварительного «Рассчитать» он не требует: считает сам.

Проект при этом не изменяется — модель передаётся расчёту отдельным контекстом. Прерванное сравнение не оставит вам проект с чужой моделью.

Мера расхождения

Главная колонка сводки — «Δlg c к опорной»: наибольший по всем частицам и всем растворам сдвиг десятичного логарифма равновесной концентрации относительно выбранной опорной модели. Величина выбрана не случайно: сдвиг $\lg c$ — это ровно то, на сколько пришлось бы подвинуть $\lg K$, чтобы одна модель дала ответ другой. То есть она отвечает на вопрос «во что обходится выбор модели» в тех единицах, которыми пользователь и мыслит.

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

Отсечка ничтожных частиц

В сравнении не участвуют частицы, чья концентрация меньше $10^{-6}$ от наибольшей концентрации в этом растворе. Причина практическая: у частицы с концентрацией $10^{-28}$ отношение между моделями может быть каким угодно, и без отсечки сводка показывала бы численный мусор вместо содержательного расхождения.

Порог относительный, а не абсолютный: абсолютный значил бы разное в разбавленной и в концентрированной задаче.

Остальные колонки

• «Согласована» — существует ли порождающая $G^E$ (см. 5.3). «нет» означает, что осмотический коэффициент в этой строке будет прочерком.
• «Ионная сила» и «φ» — для выбранного в верхней полосе раствора.
• «Δlg γ к опорной» — расхождение самих коэффициентов активности. Оно всегда больше, чем Δlg c, и разница между ними как раз и показывает, сколько взаимно уничтожилось в произведении по реакции.
• «Замечание» — отказ модели с его причиной. Отказ одной модели не обрывает сравнение: остальные считаются, а отказавшая остаётся строкой таблицы.

Вкладки «Коэффициенты активности» и «Равновесный состав» показывают те же данные по частицам: колонка на модель, строка на частицу.

Разобранный пример

Ацетатный буфер: 0.05 M $\mathrm{CH_3COOH}$ + 0.05 M $\mathrm{CH_3COONa}$ на фоне 0.25 M $\mathrm{NaCl}$, $\lg K(\mathrm{HAc}) = 4.76$, $\lg K(\mathrm{OH^-}) = -14$, $a_0 = 4.0$ Å, $b = 0$. Ионная сила выходит 0.300 — та самая целевая величина, ради которой затевался выбор модели. Ниже — то, что печатает сводка сравнения (числа посчитаны самой программой, а не взяты из литературы).

Что здесь важно увидеть.

• Расхождение сосредоточено на одной частице. Концентрации $\mathrm{Na^+}$, $\mathrm{Cl^-}$, $\mathrm{Ac^-}$ и $\mathrm{HAc}$ у всех четырёх моделей совпадают в третьем знаке — они заданы балансом и почти не зависят от $\gamma$. Весь эффект пришёлся на $\mathrm{H^+}$, то есть на ту частицу, которая связана с остальными реакцией. Это общее правило: модель действует не на состав вообще, а на равновесия внутри него.
• Δlg c = 0.055 у Дэвиса против DH — это середина второго диапазона таблицы порогов: модель уже заметна и сопоставима с погрешностью хорошо определённой константы. Иначе говоря, $\lg K$, найденная в этой системе с одной моделью, отличалась бы от найденной с другой примерно на 0.06 — больше, чем обычно указывают в качестве погрешности.
• Предельный закон даёт 0.234, то есть ошибку в четверть логарифмической единицы — при $I = 0.3$ он неприменим, и это видно не рассуждением, а числом.
• В привычных единицах расхождение выглядит так: рассчитанный pH одного и того же буфера равен 4.59 (DH и MSA), 4.63 (Дэвис) и 4.48 (предельный закон). Разброс между тремя разумными моделями — 0.03 единицы pH, что близко к воспроизводимости хорошего измерения; отставание предельного закона — 0.12, что уже грубая ошибка.

Обратите внимание, что ионная сила у всех четырёх моделей вышла одинаковой (0.300): её задаёт фоновый электролит, а не равновесие. Это типично для буферных задач и объясняет, почему внешний цикл сходится здесь за считанные итерации.

5.9. Сравнение по качеству подгонки и корреляционная матрица

В «Равновесиях» опыта нет, поэтому сравнивать можно только состав. В «Потенциометрии» и «Кондуктометрии» есть измеренная величина — и появляется второй, независимый вопрос: какая модель лучше описывает эксперимент. На него отвечают уже данные, а не рассуждение. Обе вкладки живут в окне оптимизации этих приложений.

Вкладка «Сравнение моделей»: почему мерой служит AICc, а не сумма квадратов

Каждая отмеченная модель подгоняется заново, своим набором параметров. Наборы разные, и это не мелочь реализации, а суть лестницы: у предельного закона и Дэвиса подгоняемых параметров нет вовсе, у расширенного уравнения их два, у MSA — радиусы.

Поэтому основная мера — информационный критерий Акаике с поправкой на малую выборку ^[15]:

$n$ — число экспериментальных точек, $p$ — число подгонявшихся параметров. Постоянные слагаемые правдоподобия опущены (они одинаковы у сравниваемых моделей), поэтому сами числа сравнимы только между собой, а с чужой программой — лишь по разностям. Форма записи та же, что в разделе «Математическая обработка».

Поправка $2p(p+1)/(n-p-1)$ обязательна: в электрохимическом опыте отношение $n/p$ редко доходит до 40, а без поправки критерий в этой области систематически поощряет лишние параметры.

Пороги, по которым программа выносит вывод, — общепринятые границы Бёрнхема и Андерсона ^[16]:

Колонка «Вес Акаике» $w_i = e^{-\Delta_i/2}\big/\sum_j e^{-\Delta_j/2}$ говорит то, чего не говорит сам критерий: насколько лучшая модель лучше остальных, а не только что она лучшая.

Колонка RSS стоит рядом с AICc намеренно — чтобы расхождение было видно глазами. Если порядок моделей по сумме квадратов не совпадает с порядком по AICc, программа говорит об этом отдельной строкой предупреждения: это и есть наглядное доказательство того, ради чего критерий введён.

Колонка σ равна $\sqrt{\mathrm{RSS}/(n-p)}$. При единичных весах это разброс в единицах измеряемой величины (мВ, мСм/см), и его стоит сравнить с паспортной точностью прибора. При весах $1/\text{дисперсия}$ величина безразмерна и у правильной модели должна быть около единицы: заметно больше — модель или веса неверны, заметно меньше — веса занижены.

Взвешенные суммы квадратов сравнимы между моделями только потому, что веса приходят из опыта (повторные измерения, $1/\kappa^2$), а не из модели, и у всех строк таблицы одинаковы.

Вкладка «Корреляции»: чем она отличается от доверительных интервалов

Доверительный интервал говорит, насколько плохо определён параметр сам по себе. Корреляция отвечает на другой вопрос: не определена ли вместо двух параметров лишь их комбинация. Второго по интервалам не видно вовсе — оба могут выглядеть приемлемо, а поверхность суммы квадратов быть длинным узким оврагом, вдоль которого решение гуляет свободно.

Матрица получается из уже вычисляемого обратного гессиана даром:

множитель $\sigma^2$, входящий в дисперсии, здесь сокращается, поэтому вывод матрицы не стоит ни одного лишнего вызова целевой функции.

Ячейки с $|r| \ge 0.9$ подсвечиваются, с $|r| \ge 0.99$ — выделяются красным. Прочерк в ячейке означает, что величины не существует: дисперсия по этому параметру неположительна, то есть гессиан в найденной точке не положительно определён. Это само по себе диагноз — решение не является минимумом либо параметр неопределим.

Пары, типичные именно для этого комплекса:

• $\lg K \leftrightarrow E_0$ — самая частая. Константа сдвигает расчётные активности почти равномерно по серии, а $E_0$ сдвигает всю кривую ЭДС; данные, снятые в узком диапазоне ионной силы, различить их не могут. На примере простой слабой кислоты корреляция доходит до $|r| = 1.000$.
• $\lambda^0 \leftrightarrow \lg K$ в кондуктометрии: и увеличение подвижности, и уменьшение ассоциации повышают электропроводность. Это фундаментальная особенность задачи, а не дефект подгонки, и именно она делает необходимыми одиннадцать методов оптимизации: поверхность имеет искривлённые овраги.
• Параметры модели активности ($a_0$, $b$, радиусы) $\leftrightarrow$ и то, и другое. Это прямая причина, по которой переопределение $a_0$ и $b$ по отдельным ионам выключено по умолчанию.

5.10. Как выбрать модель: порядок действий

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

Типичные случаи

6.Работа с программой

Рис. 3. Главное окно приложения: слева рельс этапов (четыре шага ввода и результаты), справа — рабочая область. Всё происходит в одном окне: модальных диалогов у ввода данных нет, а переход по рельсу проверяет заполненный шаг и объясняет отказ словами.

Работа с программой представляет собой последовательное выполнение простых операций в виде 4-х шагов мастера ввода данных. Активировать мастер можно через меню «Равновесия → Ввод исходных данных» (F5), а также автоматически — сразу после открытия существующего файла задачи («Файл → Открыть проект», Ctrl+O): мастер в этом случае сразу появляется на шаге 1, а данные шагов 2–4 уже заполнены из файла и видны сразу, без необходимости вручную подгонять размерности таблиц — редактирование остаётся доступным при необходимости.

После проверки правильности ввода данных нажатие кнопки «Готово» запускает расчёт равновесий, и таблица результатов открывается автоматически — отдельно заходить в меню «Результаты → Таблица результатов» (F6) не требуется.

7.Вывод результатов

Рис. 4. Ответ на ту же задачу: 26 растворов справочной рецептуры, в каждом — равновесные концентрации всех одиннадцати частиц, затем lg активностей частиц базиса, их коэффициенты активности и ионная сила. Каждый столбец «[частица]» отвечает своей строке стехиометрической матрицы (рис. 2) — имена совпадают буквально, потому что берутся из одного места. Формат 0.000E00 — договорённость комплекса: большая точность в эксперименте нереальна. Ширина каждого столбца подобрана по самой длинной надписи, поэтому таблица шире окна и прокручивается — так числа не обрезаются.

Рис. 5. Тот же расчёт на графике. По оси абсцисс — аналитическая концентрация Na₂HPO₄·2H₂O, то есть первого компонента из рис. 1 (вдоль справочной рецептуры буфера меняется именно она); кривые — частицы из рис. 2, сняты только Na⁺ и служебные столбцы. Ось X выбирается списком (номер раствора, любая аналитическая концентрация или ионная сила), кривые — флажками. Видно то, ради чего расчёт и делается: при смещении рецептуры от кислого края к щелочному лимонная кислота последовательно теряет три протона, а фосфат их принимает, и преобладающая форма меняется трижды.

Соответствие сквозное: строка стехиометрической матрицы (рис. 2) = столбец таблицы (рис. 4) = кривая графика (рис. 5); строка компонентной матрицы (рис. 1) = ось абсцисс графика. Ни одно из этих имён не набрано дважды — программа берёт их из одного места, и разойтись им негде.

В таблице результатов представлены: равновесные концентрации всех частиц, логарифмы активностей частиц базиса, коэффициенты активности частиц базиса и ионная сила раствора — отдельная строка на каждый раствор серии. Таблица позволяет выделять и копировать строки и столбцы в буфер обмена (Ctrl+C / Ctrl+A) для переноса в другие программы.

Программа предоставляет возможность графического анализа зависимости равновесных характеристик от выбранного аргумента. Аргументом может служить номер раствора, любая аналитическая концентрация (компонента или, в режиме прямого ввода, частицы базиса) либо ионная сила раствора. Функции $Y$ выбираются пользователем: любая одна или несколько равновесных характеристик одновременно, с возможностью логарифмического масштаба по осям.

Полное сохранение задачи включает параметры растворителя, компонентную и стехиометрическую матрицы, концентрации и таблицу результатов — как в формате CSV (секциями, для переноса в любую программу), так и в формате XLSX (отдельный лист Excel на каждый раздел).

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

8.Управление файлами

Программа предусматривает сохранение исходных данных задачи в отдельный XML-файл (*.eqp.xml), а также последующую загрузку исходной информации из такого файла. Операции доступны через меню «Файл»:

Открытие проекта сразу переводит в мастер ввода данных на шаг 1 — заходить отдельно в меню «Равновесия → Ввод исходных данных» не требуется (см. раздел «Работа с программой»).

9.Список литературы

• [1]Brinkley S. R. Jr. Calculation of the equilibrium composition of systems of many constituents. — J. Chem. Phys., 1947, v. 15, N 2, p. 107–110. Работа, с которой начался матричный формализм: равновесный состав как решение системы «закон действующих масс + материальный баланс», записанной через базис. Именно эта схема и делает расчёт в приложении не зависящим от вида задачи.
• [2]Рубцов В. И., Александров В. В. Определение параметров равновесий в растворах слабых электролитов. Журн. физ. химии, 1977, т. 51, N 4, с. 972–974. Обратная задача — определение констант по измеренному составу — для растворов слабых электролитов; здесь формализм Бринкли впервые применён к растворам, а не к газовой фазе.
• [3]Круглов В. О., Бугаевский А. А. Развитие метода Бринкли для решения прямых и обратных задач равновесной химии. // Математика в химической термодинамике. Новосибирск: Наука, 1980, с. 36–47. Развитие метода Бринкли сразу на прямую и обратную задачи. Отсюда взято устройство итерации по ионной силе, которым замкнут внешний цикл расчёта (раздел 4).
• [4]Васильев В. П., Бородин В. А., Козловский Е. В. Применение ЭВМ в химических расчетах -М., Высш, шк., 1993. -112 с. Практическое пособие по машинному расчёту равновесий: постановка задачи, выбор базиса и типичные затруднения сходимости. Полезно как введение к разделам 3 и 4.
• [5]Холин Ю. В. Количественный физико-химический анализ комплексообразования в растворах и на поверхности химически модифицированных кремнезёмов - содержательные модели, математические методы и их приложения. — Харьков, Фолио, 2000. - 288 с. Монография о содержательных моделях комплексообразования и математических методах их проверки: чем «модель, описывающая данные» отличается от модели, которой можно верить. Тот же вопрос решают разделы 5.8–5.10.
• [6]Debye P., Hückel E. Zur Theorie der Elektrolyte. I. Gefrierpunktserniedrigung und verwandte Erscheinungen. — Phys. Z., 1923, Bd. 24, N 9, S. 185–206. Исходная работа теории сильных электролитов: линеаризованное уравнение Пуассона–Больцмана для точечных ионов и предельный закон $\lg\gamma = -Az^2\sqrt{I}$ — первая ступень лестницы моделей (раздел 5.2) и предел, к которому обязаны сходиться все остальные при разбавлении.
• [7]Davies C. W. Ion Association. — London: Butterworths, 1962. — 190 p. Уравнение Дэвиса — вторая ступень лестницы: эмпирический линейный член вместо индивидуального размера иона, зато без единого подгоняемого параметра. По этой причине оно и выбрано моделью иллюстративных расчётов комплекса.
• [8]Blum L. Mean spherical model for asymmetric electrolytes. I. Method of solution. — Mol. Phys., 1975, v. 30, N 5, p. 1529–1535. Аналитическое решение уравнения Орнштейна–Цернике для смеси заряженных твёрдых сфер в среднесферическом замыкании — теоретическая основа MSA (раздел 5.5).
• [9]Roger G. M., Durand-Vidal S., Bernard O., Turq P. Electrical conductivity of mixed electrolytes: modeling within the mean spherical approximation. — J. Phys. Chem. B, 2009, v. 113, N 25, p. 8670–8674. Рабочая форма MSA, по которой написан код: уравнения (18)–(21) для параметра экранирования $\Gamma$ и величин $\Delta$, $\Omega$, $P_n$, $N_i$. Те же уравнения использует транспортная MSA в «Кондуктометрии» — отсюда правило «один размер иона служит и $\gamma$, и $\Lambda$».
• [10]Pitzer K. S. Thermodynamics of electrolytes. I. Theoretical basis and general equations. — J. Phys. Chem., 1973, v. 77, N 2, p. 268–277. Вириальное разложение Питцера — общепринятый стандарт для концентрированных растворов. В комплексе сознательно не используется: его коэффициенты уже вобрали ионную ассоциацию, а здесь комплексы задаются явно, и учёт вышел бы двойным (раздел 5.4).
• [11]Ciavatta L. The specific interaction theory in evaluating ionic equilibria. — Ann. Chim. (Rome), 1980, v. 70, p. 551–567. Теория специфического ионного взаимодействия (SIT). Отвергнута по другой причине, чем Питцер: её база параметров водная, а комплекс рассчитан на любой растворитель (раздел 5.4).
• [12]Robinson R. A., Stokes R. H. Electrolyte Solutions. 2nd ed., revised. — London: Butterworths, 1965. — 571 p. Рус. пер.: Робинсон Р., Стокс Р. Растворы электролитов. — М.: Изд-во иностр. лит., 1963. — 647 с. Классическая монография по растворам электролитов; здесь — источник гидратационной поправки, которая тоже отвергнута: она выведена для среднего $\gamma_\pm$ одной соли и на смесь не переносится (раздел 5.4).
• [13]Gronwall T. H., La Mer V. K., Sandved K. Über den Einfluss der sogenannten höheren Glieder in der Debye–Hückelschen Theorie der Lösungen starker Elektrolyte. — Phys. Z., 1928, Bd. 29, S. 358–393. Там же по этому вопросу: Dutta M. On the calculation of the activity-coefficient. — Naturwissenschaften, 1953, Bd. 40, S. 481. Строгое решение нелинейного уравнения Пуассона–Больцмана рядом по заряду. Приведено как факультативная ступень: математически безупречно, но опирается на ту же примитивную модель среды, что и предельный закон, — MSA достигает большего и без рядов.
• [14]Mansoori G. A., Carnahan N. F., Starling K. E., Leland T. W. Equilibrium thermodynamic properties of the mixture of hard spheres. — J. Chem. Phys., 1971, v. 54, N 4, p. 1523–1525. Уравнение состояния смеси твёрдых сфер (BMCSL). Даёт твердосферную часть MSA, а вместе с ней — высаливание нейтральных частиц без константы Сеченова (раздел 5.5).
• [15]Hurvich C. M., Tsai C.-L. Regression and time series model selection in small samples. — Biometrika, 1989, v. 76, N 2, p. 297–307. Поправка Акаике на малую выборку (AICc). Обязательна здесь: в химическом опыте число точек редко превышает число параметров в сорок раз, а без поправки критерий систематически выбирает более богатую модель (раздел 5.9).
• [16]Burnham K. P., Anderson D. R. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd ed. — New York: Springer, 2002. — 488 p. Источник порогов, по которым программа выносит вывод о различии моделей (разности AICc в 2 и 10 единиц) и весов Акаике.