From fe8146d3359498de9c93cc1223d43cd51a2a0924 Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Fri, 2 Aug 2024 14:21:58 +0200 Subject: [PATCH] Add Sips isotherm MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Based on changes made by samuel.leweke@bayer.com Convert Sips to standard for Freundlich isotherms in CADET, which is ^(1/n) instead of ^n Improve stability by extending linearization range into negative numbers Update documentation Update units in documentation Extend tests Co-authored-by: Ronald Jäpel --- doc/interface/binding/index.rst | 1 + .../binding/multi_component_sips.rst | 85 ++++++ doc/literature.bib | 12 + doc/modelling/binding/index.rst | 5 + .../binding/multi_component_sips.rst | 32 +++ src/libcadet/BindingModelFactory.cpp | 2 + src/libcadet/CMakeLists.txt | 5 +- src/libcadet/model/binding/SipsBinding.cpp | 260 ++++++++++++++++++ test/BindingModels.cpp | 60 ++++ 9 files changed, 460 insertions(+), 2 deletions(-) create mode 100644 doc/interface/binding/multi_component_sips.rst create mode 100644 doc/modelling/binding/multi_component_sips.rst create mode 100644 src/libcadet/model/binding/SipsBinding.cpp diff --git a/doc/interface/binding/index.rst b/doc/interface/binding/index.rst index 3593970ba..186825533 100644 --- a/doc/interface/binding/index.rst +++ b/doc/interface/binding/index.rst @@ -69,6 +69,7 @@ This group also takes precedence over a possibly existing ``/input/model/unit_XX self_association freundlich_ldf multi_component_ldf_freundlich + multi_component_sips hic_water_on_hydrophobic_surfaces hic_constant_water_activity multi_component_colloidal diff --git a/doc/interface/binding/multi_component_sips.rst b/doc/interface/binding/multi_component_sips.rst new file mode 100644 index 000000000..f9da75dac --- /dev/null +++ b/doc/interface/binding/multi_component_sips.rst @@ -0,0 +1,85 @@ +.. _multi_component_sips_config: + +Multi Component Sips +~~~~~~~~~~~~~~~~~~~~ + +**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = MULTI_COMPONENT_SIPS** + +For information on model equations, refer to :ref:`multi_component_sips_model`. + +``IS_KINETIC`` + Selects kinetic or quasi-stationary adsorption mode: 1 = kinetic, 0 = + quasi-stationary. If a single value is given, the mode is set for all + bound states. Otherwise, the adsorption mode is set for each bound + state separately. + +=================== ========================= ========================================= +**Type:** int **Range:** {0,1} **Length:** 1/NTOTALBND +=================== ========================= ========================================= + +``SIPS_KA`` + Adsorption rate constants + +**Unit:** :math:`m_{MP}^3~mol^{-1}~s^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + +``SIPS_KD`` + Desorption rate constants + +**Unit:** :math:`s^{-1}` + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ================================== + +``SIPS_QMAX`` + Maximum adsorption capacities + +**Unit:** :math:`mol~m_{SP}^{-3}` + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\gt 0` **Length:** NCOMP +=================== ========================= ================================== + +``SIPS_EXP`` + Freundlich-type exponent + +=================== ========================= ================================== +**Type:** double **Range:** :math:`\gt 0` **Length:** NCOMP +=================== ========================= ================================== + +``SIPS_REFC0`` + Reference liquid phase concentration (optional, defaults to + :math:`1.0`) + + +**Unit:** :math:`mol~m_{MP}^{-3}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\gt 0` **Length:** 1 +=================== ========================= ========================================= + +``SIPS_REFQ`` + Reference solid phase concentration (optional, defaults to + :math:`1.0`) + + +**Unit:** :math:`mol~m_{SP}^{-3}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\gt 0` **Length:** 1 +=================== ========================= ========================================= + +``SIPS_LINEAR_THRESHOLD`` + Threshold for linearization. Originally, the liquid phase concentration + enters the adsorption rate via a power term. Below this threshold, a + quadratic approximation is used instead. This ensures numerical stability. + +**Unit:** :math:`mol~m_{MP}^{-3}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\gt 0` **Length:** 1 +=================== ========================= ========================================= diff --git a/doc/literature.bib b/doc/literature.bib index 7c7cdc16a..a4aad059a 100644 --- a/doc/literature.bib +++ b/doc/literature.bib @@ -634,3 +634,15 @@ @article{Xu2009 author = {Xuankuo Xu and Abraham M. Lenhoff}, keywords = {Protein adsorption isotherms, Competitive adsorption, Protein–surface interactions, Protein–protein interactions, Colloidal interactions, Ion-exchange chromatography}, } +@article{sips1948, + title = {On the {Structure} of a {Catalyst} {Surface}}, + author = {Sips, Robert}, + year = {1948}, + journal = {The Journal of Chemical Physics}, + pages = {490--495}, + volume = {16}, + issn = {0021-9606}, + doi = {https://doi.org/10.1063/1.1746922}, + number = {5}, + month = may +} \ No newline at end of file diff --git a/doc/modelling/binding/index.rst b/doc/modelling/binding/index.rst index d0866ece9..a7576699b 100644 --- a/doc/modelling/binding/index.rst +++ b/doc/modelling/binding/index.rst @@ -277,6 +277,11 @@ The models also differ in whether a mobile phase modifier (e.g., salt) is suppor - x - ✓ - x + * - :ref:`multi_component_sips_model` + - ✓ + - × + - ✓ + - × * - :ref:`multi_component_ldf_freundlich_model` - ✓ - × diff --git a/doc/modelling/binding/multi_component_sips.rst b/doc/modelling/binding/multi_component_sips.rst new file mode 100644 index 000000000..b213dc5db --- /dev/null +++ b/doc/modelling/binding/multi_component_sips.rst @@ -0,0 +1,32 @@ +.. _multi_component_sips_model: + +Multi Component Sips +~~~~~~~~~~~~~~~~~~~~~~~~ + +The Sips binding model is a combination of the :ref:`Freundlich` and the :ref:`Langmuir adsorption model`. + +.. math:: + + \begin{aligned} + \frac{\mathrm{d} q_i}{\mathrm{d} t} = k_{a,i}\: \left( \frac{c_{p,i}}{ c_{\text{ref}} }\right)^{1 / n_i}\: q_{\text{max},i} \left( 1 - \sum_{j=0}^{N_{\text{comp}} - 1} \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} \left( \frac{q_i}{q_{\text{ref}}} \right) && i = 0, \dots, N_{\text{comp}} - 1. + \end{aligned} + + +Here, :math:`c_{\text{ref}}` is a :ref:`reference concentration `, :math:`n_i` is the Freundlich exponent, :math:`k_{a,i}, k_{d,i}` are the adsorption and desorption rates, and :math:`q_{\text{max},j}` is the adsorption capacity. + +As for the :ref:`Freundlich` isotherm, the first order Jacobian :math:`\left(\frac{dq^*}{dc_p}\right)` tends to infinity as :math:`c_{p} \rightarrow 0` for :math:`n>1`. +Additionally, the isotherm is undefined for :math:`c_{p} < 0` if :math:`\frac{1}{n_i}` can be expressed as :math:`\frac{p}{q}` with :math:`p,q \in \mathbb{N}` where :math:`q` is an even number. +Negative concentrations can arise during simulations due to numerical fluctuations. +To address these issues an approximation of the isotherm is considered below a threshold concentration :math:`c_p < \varepsilon`. +This approximation matches the isotherm in such a way that :math:`q=0` at :math:`c_p=0` and also matches the value and the first derivative of the istotherm at :math:`c_p = \varepsilon`, where :math:`\varepsilon` is a very small number, for example :math:`1e-10`. +The form of approximation and its derivative is given below for :math:`c_p < \varepsilon`: + +.. math:: + + \begin{aligned} + c_{p,i,lin} &= \left(\frac{\varepsilon}{c_{p,\text{ref}}}\right)^{\frac{1}{n_i} - 2} \frac{c_{p,i}}{{c_{p,\text{ref}}}^2} \left( \left(2-\frac{1}{n_i}\right)\varepsilon + c_{p,i}\left(\frac{1}{n}-1\right) \right) \\ + \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i,lin} q_{\text{max},i} \left( 1 - \sum_j \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} \frac{q_i}{q_{i,\text{ref}}} + \end{aligned} + +For more information on model parameters required to define in CADET file format, see :ref:`multi_component_sips_config`. +For more information on the model and its origin, please refer to :cite:`sips1948`. diff --git a/src/libcadet/BindingModelFactory.cpp b/src/libcadet/BindingModelFactory.cpp index 3a0d92cd7..a2e9d0d90 100644 --- a/src/libcadet/BindingModelFactory.cpp +++ b/src/libcadet/BindingModelFactory.cpp @@ -45,6 +45,7 @@ namespace cadet void registerBiLangmuirLDFModel(std::unordered_map>& bindings); void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map>& bindings); void registerHICConstantWaterActivityModel(std::unordered_map>& bindings); + void registerSipsModel(std::unordered_map>& bindings); void registerMultiComponentLDFFreundlichModel(std::unordered_map>& bindings); } } @@ -76,6 +77,7 @@ namespace cadet model::binding::registerBiLangmuirLDFModel(_bindingModels); model::binding::registerHICWaterOnHydrophobicSurfacesModel(_bindingModels); model::binding::registerHICConstantWaterActivityModel(_bindingModels); + model::binding::registerSipsModel(_bindingModels); model::binding::registerMultiComponentLDFFreundlichModel(_bindingModels); registerModel(); } diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 3df020e5c..8765ebe40 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -100,6 +100,7 @@ set(LIBCADET_BINDINGMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/BiLangmuirLDFBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/SipsBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/MultiComponentLDFFreundlichBinding.cpp ) set(LIBCADET_EXCHANGEMODEL_SOURCES @@ -290,7 +291,7 @@ if (LAPACK_FOUND) add_library(libcadet_static STATIC $) set_target_properties(libcadet_static PROPERTIES OUTPUT_NAME cadet_static) target_link_libraries(libcadet_static PUBLIC CADET::CompileOptions CADET::LibOptions PRIVATE CADET::AD libcadet_nonlinalg_static SUNDIALS::sundials_idas ${SUNDIALS_NVEC_TARGET} ${TBB_TARGET} ${EIGEN_TARGET}) - + # --------------------------------------------------- # Build the shared library # --------------------------------------------------- @@ -298,7 +299,7 @@ if (LAPACK_FOUND) add_library(libcadet_shared SHARED $) set_target_properties(libcadet_shared PROPERTIES OUTPUT_NAME cadet) target_link_libraries (libcadet_shared PUBLIC CADET::CompileOptions CADET::LibOptions PRIVATE CADET::AD libcadet_nonlinalg_static SUNDIALS::sundials_idas ${SUNDIALS_NVEC_TARGET} ${TBB_TARGET} ${EIGEN_TARGET}) - + list(APPEND LIBCADET_TARGETS libcadet_nonlinalg_static libcadet_object libcadet_static libcadet_shared) unset(LIB_LAPACK_DEFINE) diff --git a/src/libcadet/model/binding/SipsBinding.cpp b/src/libcadet/model/binding/SipsBinding.cpp new file mode 100644 index 000000000..d4e340c9e --- /dev/null +++ b/src/libcadet/model/binding/SipsBinding.cpp @@ -0,0 +1,260 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-present: The CADET-Core Authors +// Please see the AUTHORS.md file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +#include "model/binding/BindingModelBase.hpp" +#include "model/ExternalFunctionSupport.hpp" +#include "model/ModelUtils.hpp" +#include "cadet/Exceptions.hpp" +#include "model/Parameters.hpp" +#include "LocalVector.hpp" +#include "SimulationTypes.hpp" +#include "MathUtil.hpp" + +#include +#include +#include +#include + +/* +{ + "name": "SipsParamHandler", + "externalName": "ExtSipsParamHandler", + "parameters": + [ + { "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "SIPS_KA"}, + { "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "SIPS_KD"}, + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "SIPS_QMAX"}, + { "type": "ScalarComponentDependentParameter", "varName": "n", "confName": "SIPS_EXP"} + ], + "constantParameters": + [ + { "type": "ReferenceConcentrationParameter", "varName": ["refC", "refQ"], "objName": "refConcentration", "confPrefix": "SIPS_"}, + { "type": "ScalarParameter", "varName": "linearThreshold", "confName": "SIPS_LINEAR_THRESHOLD"} + ] +} +*/ + +namespace cadet +{ + +namespace model +{ + +inline const char* SipsParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_SIPS"; } + +inline bool SipsParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("SIPS_KA, SIPS_KD, and SIPS_QMAX have to have the same size"); + + return true; +} + +inline const char* ExtSipsParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_SIPS"; } + +inline bool ExtSipsParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("EXT_SIPS_KA, EXT_SIPS_KD, and EXT_SIPS_QMAX have to have the same size"); + + return true; +} + + +/** + * @brief Defines the multi component Sips binding model + * @details Implements the Sips adsorption model: \f[ \begin{align} + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} \left(\frac{c_{p,i}}{c_{p,\text{ref}}} \right)^{\frac{1}{n_i}} q_{\text{max},i} \left( 1 - \sum_j \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} \frac{q_i}{q_{i,\text{ref}}} + * \end{align} \f] + * with a quadratic model approximation around the "linear" threshold for the Cp state: \f[ \begin{align} + * c_{p,i,lin} &= \left(\frac{\varepsilon}{c_{p,\text{ref}}}\right)^{\frac{1}{n_i} - 2} \frac{c_{p,i}}{{c_{p,\text{ref}}}^2} \left( \left(2-\frac{1}{n_i}\right)\varepsilon + c_{p,i}\left(\frac{1}{n}-1\right) \right) \\ + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i,lin} q_{\text{max},i} \left( 1 - \sum_j \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} \frac{q_i}{q_{i,\text{ref}}} + * \end{align} \f] + * Multiple bound states are not supported. + * Components without bound state (i.e., non-binding components) are supported. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ +template +class SipsBindingBase : public ParamHandlerBindingModelBase +{ +public: + + SipsBindingBase() { } + virtual ~SipsBindingBase() CADET_NOEXCEPT { } + + static const char* identifier() { return ParamHandler_t::identifier(); } + + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + + CADET_BINDINGMODELBASE_BOILERPLATE + +protected: + using ParamHandlerBindingModelBase::_paramHandler; + using ParamHandlerBindingModelBase::_reactionQuasistationarity; + using ParamHandlerBindingModelBase::_nComp; + using ParamHandlerBindingModelBase::_nBoundStates; + + template + int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, + CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + using CpStateParamType = typename DoubleActivePromoter::type; + + // Protein fluxes: -k_{a,i} * (c_{p,i} / c_{p,ref})^n_i * q_{max,i} * (1 - \sum_j q_j / q_{max,j}) + k_{d,i} * q_i / q_{ref} + ResidualType qSum = 1.0; + unsigned int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + qSum -= y[bndIdx] / static_cast(p->qMax[i]); + + // Next bound component + ++bndIdx; + } + + const ParamType refC = static_cast(p->refC); + const ParamType refQ = static_cast(p->refQ); + const double linearThreshold = static_cast(p->linearThreshold); + + bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + const ParamType n_i = static_cast(p->n[i]); + // Residual + if (yCp[i] <= linearThreshold) + { + // Linearize + const CpStateParamType cpLin = pow(linearThreshold / refC, 1.0 / n_i - 2.0) * yCp[i] / sqr(refC) * ((2.0 - 1.0 / n_i) * linearThreshold + yCp[i] * (1.0 / n_i - 1.0)); + res[bndIdx] = static_cast(p->kD[i]) * y[bndIdx] / refQ - static_cast(p->kA[i]) * cpLin * static_cast(p->qMax[i]) * qSum; + } + else + { + res[bndIdx] = static_cast(p->kD[i]) * y[bndIdx] / refQ - static_cast(p->kA[i]) * pow(yCp[i] / refC, 1.0 / n_i) * static_cast(p->qMax[i]) * qSum; + } + + // Next bound component + ++bndIdx; + } + + return 0; + } + + template + void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein fluxes: -k_{a,i} * (c_{p,i} / c_{p,ref})^n_i * q_{max,i} * (1 - \sum_j q_j / q_{max,j}) + k_{d,i} * q_i / q_{ref} + double qSum = 1.0; + int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + qSum -= y[bndIdx] / static_cast(p->qMax[i]); + + // Next bound component + ++bndIdx; + } + + const double refC = static_cast(p->refC); + const double refQ = static_cast(p->refQ); + const double linearThreshold = static_cast(p->linearThreshold); + + bndIdx = 0; + for (int i = 0; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + const double ka = static_cast(p->kA[i]); + const double kd = static_cast(p->kD[i]); + const double n = static_cast(p->n[i]); + + // dres_i / dc_{p,i} + if (yCp[i] <= linearThreshold) + { + // Linearized + jac[i - bndIdx - offsetCp] = -ka * static_cast(p->qMax[i]) * qSum * pow(linearThreshold / refC, 1.0 / n) / linearThreshold * (2.0 - 1.0 / n + 2.0 * (1.0 / n - 1.0) * yCp[i] / linearThreshold); + // Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}. + // This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}. + } + else + { + jac[i - bndIdx - offsetCp] = -ka * 1.0 / n * static_cast(p->qMax[i]) * qSum * pow(yCp[i] / refC, 1.0 / n - 1.0) / refC; + // Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}. + // This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}. + } + + double linearizedCp = 0; + + // Fill dres_i / dq_j + if (yCp[i] <= linearThreshold) + { + linearizedCp = pow(linearThreshold / refC, 1.0 / n - 2.0) * yCp[i] / sqr(refC) * ((2.0 - 1.0 / n) * linearThreshold + yCp[i] * (1.0 / n - 1.0)); + } + else + { + linearizedCp = pow(yCp[i] / refC, 1.0 / n); + } + + const double factor = ka * linearizedCp * static_cast(p->qMax[i]); + + int bndIdx2 = 0; + for (int j = 0; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + // dres_i / dq_j + jac[bndIdx2 - bndIdx] = factor / static_cast(p->qMax[j]); + // Getting to q_j: -bndIdx takes us to q_0, another +bndIdx2 to q_j. This means jac[bndIdx2 - bndIdx] corresponds to q_j. + + ++bndIdx2; + } + + // Add to dres_i / dq_i + jac[0] += kd / refQ; + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } +}; + +typedef SipsBindingBase SipsBinding; +typedef SipsBindingBase ExternalSipsBinding; + +namespace binding +{ + void registerSipsModel(std::unordered_map>& bindings) + { + bindings[SipsBinding::identifier()] = []() { return new SipsBinding(); }; + bindings[ExternalSipsBinding::identifier()] = []() { return new ExternalSipsBinding(); }; + } +} // namespace binding + +} // namespace model + +} // namespace cadet diff --git a/test/BindingModels.cpp b/test/BindingModels.cpp index c4c1eebfd..d0aae3c63 100644 --- a/test/BindingModels.cpp +++ b/test/BindingModels.cpp @@ -1024,6 +1024,66 @@ CADET_BINDINGTEST("MULTI_COMPONENT_LDF_FREUNDLICH", "EXT_MULTI_COMPONENT_LDF_FRE )json", \ 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_USED, CADET_DONT_COMPARE_BINDING_VS_NONBINDING) +CADET_BINDINGTEST("MULTI_COMPONENT_SIPS", "EXT_MULTI_COMPONENT_SIPS", (1, 1, 1, 1), (1, 0, 1, 1, 1), (1.0, 0.0, 2.0, -1e-4, 0.0, 0.0, 0.0, 0.0), (1.0, 3.0, 0.0, 2.0, -1e-4, 0.0, 0.0, 0.0, 0.0), \ + R"json( "SIPS_KA": [1.14, 2.0, 1.14, 2.0], + "SIPS_KD": [0.004, 0.008, 0.004, 0.008], + "SIPS_QMAX": [4.88, 3.5, 4.88, 3.5], + "SIPS_EXP": [1.2, 1.6, 1.2, 1.6], + "SIPS_REFC0": 2.0, + "SIPS_REFQ": 1.1, + "SIPS_LINEAR_THRESHOLD": 1e-10 + )json", \ + R"json( "SIPS_KA": [1.14, 1.0, 2.0, 1.14, 2.0], + "SIPS_KD": [0.004, 2.0, 0.008, 0.004, 0.008], + "SIPS_QMAX": [4.88, 3.0, 3.5, 4.88, 3.5], + "SIPS_EXP": [1.2, 1.5, 1.6, 1.2, 1.6], + "SIPS_REFC0": 2.0, + "SIPS_REFQ": 1.1, + "SIPS_LINEAR_THRESHOLD": 1e-10 + )json", \ + R"json( "EXT_SIPS_KA": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KA_T": [1.14, 2.0, 1.14, 2.0], + "EXT_SIPS_KA_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KA_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KD": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KD_T": [0.004, 0.008, 0.004, 0.008], + "EXT_SIPS_KD_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KD_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_QMAX": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_QMAX_T": [4.88, 3.5, 4.88, 3.5], + "EXT_SIPS_QMAX_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_QMAX_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_EXP": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_EXP_T": [1.2, 1.6, 1.2, 1.6], + "EXT_SIPS_EXP_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_EXP_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_REFC0": 2.0, + "EXT_SIPS_REFQ": 1.1, + "EXT_SIPS_LINEAR_THRESHOLD": 1e-10 + )json", \ + R"json( "EXT_SIPS_KA": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KA_T": [1.14, 1.0, 2.0, 1.14, 2.0], + "EXT_SIPS_KA_TT": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KA_TTT": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KD": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KD_T": [0.004, 2.0, 0.008, 0.004, 0.008], + "EXT_SIPS_KD_TT": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_KD_TTT": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_QMAX": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_QMAX_T": [4.88, 3.0, 3.5, 4.88, 3.5], + "EXT_SIPS_QMAX_TT": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_QMAX_TTT": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_EXP": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_EXP_T": [1.2, 1.5, 1.6, 1.2, 1.6], + "EXT_SIPS_EXP_TT": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_EXP_TTT": [0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_SIPS_REFC0": 2.0, + "EXT_SIPS_REFQ": 1.1, + "EXT_SIPS_LINEAR_THRESHOLD": 1e-10 + )json", \ + 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING) + + CADET_BINDINGTEST("MULTI_COMPONENT_SPREADING", "EXT_MULTI_COMPONENT_SPREADING", (2,2), (2,0,2), (1.2, 1.5, 0.1, 0.2, 0.3, 0.4), (1.2, 0.5, 1.5, 0.1, 0.2, 0.3, 0.4), \ R"json( "MCSPR_KA": [1.14, 2.0, 1.5, 1.9], "MCSPR_KD": [0.004, 0.008, 0.006, 0.002],