diff --git a/doc/interface/binding/hic_constant_water_activity.rst b/doc/interface/binding/hic_constant_water_activity.rst index 8956c9a65..cbb10fc46 100644 --- a/doc/interface/binding/hic_constant_water_activity.rst +++ b/doc/interface/binding/hic_constant_water_activity.rst @@ -1,7 +1,7 @@ .. _hic_constant_water_activity_config: -Constant Water Activity -~~~~~~~~~~~~~~~~~~~~~~~ +HIC Constant Water Activity +~~~~~~~~~~~~~~~~~~~~~~~~~~~ **Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = HIC_CONSTANT_WATER_ACTIVITY** diff --git a/doc/interface/binding/hic_water_on_hydrophobic_surfaces.rst b/doc/interface/binding/hic_water_on_hydrophobic_surfaces.rst index 1044b87b7..c7730c853 100644 --- a/doc/interface/binding/hic_water_on_hydrophobic_surfaces.rst +++ b/doc/interface/binding/hic_water_on_hydrophobic_surfaces.rst @@ -1,7 +1,7 @@ .. _hic_water_on_hydrophobic_surfaces_config: -Water on Hydrophobic Surfaces -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +HIC Water on Hydrophobic Surfaces +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = HIC_WATER_ON_HYDROPHOBIC_SURFACES** diff --git a/doc/literature.bib b/doc/literature.bib index bea433e67..de1673281 100644 --- a/doc/literature.bib +++ b/doc/literature.bib @@ -414,7 +414,7 @@ @article{Benedikt2019 URL = { https://doi.org/10.1021/acs.est.8b05873}, eprint = { https://doi.org/10.1021/acs.est.8b05873} } -@article{WANG201671, +@article{Wang2016, title = {Water on hydrophobic surfaces: Mechanistic modeling of hydrophobic interaction chromatography}, journal = {Journal of Chromatography A}, volume = {1465}, @@ -425,7 +425,7 @@ @article{WANG201671 url = {https://www.sciencedirect.com/science/article/pii/S0021967316310263}, author = {Gang Wang and Tobias Hahn and Jürgen Hubbuch}, } -@article{JAPEL2022463408, +@article{Jaepel2022, title = {Bayesian optimization using multiple directional objective functions allows the rapid inverse fitting of parameters for chromatography simulations}, journal = {Journal of Chromatography A}, volume = {1679}, diff --git a/doc/modelling/binding/hic_constant_water_activity.rst b/doc/modelling/binding/hic_constant_water_activity.rst index f6ec34d3a..b1dc8fa73 100644 --- a/doc/modelling/binding/hic_constant_water_activity.rst +++ b/doc/modelling/binding/hic_constant_water_activity.rst @@ -1,12 +1,12 @@ .. _hic_constant_water_activity_model: -Constant Water Activity -~~~~~~~~~~~~~~~~~~~~~~~ -This model implemments the HIC Isotherm assuming a constant water activity as described by Jäpel and Buyel :cite:j 2022 +HIC Constant Water Activity +~~~~~~~~~~~~~~~~~~~~~~~~~~~ +This model implemments the HIC isotherm assuming a constant water activity as described by Jäpel and Buyel :cite:`Jaepel2022`. .. math:: \begin{align} - \beta&=\beta_0 e^{c_{p,0}\beta_1}\\ + \beta &= \beta_0 e^{c_{p,0}\beta_1} \\ \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - k_{d,i} q_i 0.1^{\nu_i \beta} \end{align} diff --git a/doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst b/doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst index 8cb3cbee9..f14deff60 100644 --- a/doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst +++ b/doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst @@ -1,21 +1,20 @@ .. _hic_water_on_hydrophobic_surfaces_model: -Water on Hydrophobic Surfaces -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +HIC Water on Hydrophobic Surfaces +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -This model implements a slightly modified version of the HIC Isotherm by Wang et al. based on their 2016 paper :cite:`WANG201671`. +This model implements a slightly modified version of the HIC isotherm by Wang et al. based on their 2016 paper :cite:`Wang2016`. A naive multicomponent version was added that reduces to the original formulation if only 1 binding species is present. .. math:: \begin{align} - \beta&=\beta_0 e^{c_{p,0}\beta_1}\\ - \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - - k_{d,i} q_i \left(\sum_j q_j \right)^{\nu_i \beta} + \beta &= \beta_0 e^{c_{p,0}\beta_1} \\ + \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - k_{d,i} q_i \left(\sum_j q_j \right)^{\nu_i \beta} \end{align} -Component :math:`c_0` is assumed to be salt without a bound state. -Multiple bound states are not supported. -Components without bound state (i.e., salt and non-binding components) are supported. +- Component :math:`c_0` is assumed to be salt without a bound state. +- Multiple bound states are not supported. +- Components without bound state (i.e., salt and non-binding components) are supported. For more information on model parameters required to define in CADET file format, see :ref:`hic_water_on_hydrophobic_surfaces_config`. diff --git a/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp b/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp index d152093f0..e1bf7f309 100644 --- a/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp +++ b/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp @@ -1,7 +1,7 @@ // ============================================================================= // CADET // -// Copyright © 2008-2021: The CADET Authors +// Copyright © 2008-2023: The CADET Authors // Please see the AUTHORS and CONTRIBUTORS file. // // All rights reserved. This program and the accompanying materials @@ -20,6 +20,7 @@ #include "LocalVector.hpp" #include "SimulationTypes.hpp" +#include #include #include #include @@ -37,6 +38,10 @@ { "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "HICCWA_KD"}, { "type": "ScalarComponentDependentParameter", "varName": "nu", "confName": "HICCWA_NU"}, { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "HICCWA_QMAX"} + ], + "constantParameters": + [ + { "type": "ScalarParameter", "varName": "waterActivity", "default": 0.1, "confName": "HICCWA_WATER_ACTIVITY"} ] } */ @@ -54,235 +59,228 @@ namespace cadet { - namespace model - { - - inline const char* HICCWAParamHandler::identifier() CADET_NOEXCEPT { return "HIC_CONSTANT_WATER_ACTIVITY"; } - - inline bool HICCWAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) - { - if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) - throw InvalidParameterException("HICCWA_KA, HICCWA_KD, HICCWA_NU, and HICCWA_QMAX have to have the same size"); - - return true; - } - - inline const char* ExtHICCWAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_HIC_CONSTANT_WATER_ACTIVITY"; } - - inline bool ExtHICCWAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) - { - if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) - throw InvalidParameterException("KA, KD, NU, and QMAX have to have the same size"); - - return true; - } +namespace model +{ +inline const char* HICCWAParamHandler::identifier() CADET_NOEXCEPT { return "HIC_CONSTANT_WATER_ACTIVITY"; } - /** - * @brief Defines the HIC Isotherm assuming a constant water activity as described by Jäpel and Buyel 2022 - * @details Implements the the HIC Isotherm assuming a constant water activity: \f[ \begin{align} - * \beta&=\beta_0 e^{c_{p,0}\beta_1}\\ - * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - k_{d,i} q_i 0.1^{\nu_i \beta} - * \end{align} \f] - * Component @c 0 is assumed to be salt without a bound state. Multiple bound states are not supported. - * Components without bound state (i.e., salt and non-binding components) are supported. - * - * See @cite Jäpel and Buyel 2022. - * @tparam ParamHandler_t Type that can add support for external function dependence - */ - template - class ConstantWaterActivityBindingBase : public ParamHandlerBindingModelBase - { - public: +inline bool HICCWAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("HICCWA_KA, HICCWA_KD, HICCWA_NU, and HICCWA_QMAX have to have the same size"); - ConstantWaterActivityBindingBase() { } - virtual ~ConstantWaterActivityBindingBase() CADET_NOEXCEPT { } + return true; +} - static const char* identifier() { return ParamHandler_t::identifier(); } +inline const char* ExtHICCWAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_HIC_CONSTANT_WATER_ACTIVITY"; } - virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) - { - const bool res = BindingModelBase::configureModelDiscretization(paramProvider, nComp, nBound, boundOffset); +inline bool ExtHICCWAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("KA, KD, NU, and QMAX have to have the same size"); - // Guarantee that salt has no bound state - if (nBound[0] != 0) - throw InvalidParameterException("HICCWA binding model requires exactly zero bound states for salt component"); + return true; +} - // First flux is salt, which is always quasi-stationary - _reactionQuasistationarity[0] = false; - return res; - } +/** + * @brief Defines the HIC Isotherm assuming a constant water activity as described by Jäpel and Buyel, 2022 + * @details Implements the the HIC Isotherm assuming a constant water activity: \f[ \begin{align} + * \beta &= \beta_0 e^{c_{p,0}\beta_1}\\ + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - k_{d,i} q_i 0.1^{\nu_i \beta} + * \end{align} \f] + * Component @c 0 is assumed to be salt without a bound state. Multiple bound states are not supported. + * Components without bound state (i.e., salt and non-binding components) are supported. + * + * See @cite Jaepel2022. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ +template +class ConstantWaterActivityBindingBase : public ParamHandlerBindingModelBase +{ +public: - virtual bool hasSalt() const CADET_NOEXCEPT { return true; } - virtual bool supportsMultistate() const CADET_NOEXCEPT { return false; } - virtual bool supportsNonBinding() const CADET_NOEXCEPT { return true; } - virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT { return false; } - virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + ConstantWaterActivityBindingBase() { } + virtual ~ConstantWaterActivityBindingBase() CADET_NOEXCEPT { } - virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const - { - typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); - return true; - } - - virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const - { - preConsistentInitialState(t, secIdx, colPos, y, yCp, workSpace); - } + static const char* identifier() { return ParamHandler_t::identifier(); } + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) + { + const bool res = BindingModelBase::configureModelDiscretization(paramProvider, nComp, nBound, boundOffset); - CADET_BINDINGMODELBASE_BOILERPLATE + // Guarantee that salt has no bound state + if (nBound[0] != 0) + throw InvalidParameterException("HICCWA binding model requires exactly zero bound states for salt component"); - protected: - using ParamHandlerBindingModelBase::_paramHandler; - using ParamHandlerBindingModelBase::_reactionQuasistationarity; - using ParamHandlerBindingModelBase::_nComp; - using ParamHandlerBindingModelBase::_nBoundStates; + // First flux is salt, which is always quasi-stationary + _reactionQuasistationarity[0] = false; - template - int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, - CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const - { - using CpStateParamType = typename DoubleActivePromoter::type; - using StateParamType = typename DoubleActivePromoter::type; + return res; + } - typename ParamHandler_t::ParamsHandle const _p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + virtual bool hasSalt() const CADET_NOEXCEPT { return true; } + virtual bool supportsMultistate() const CADET_NOEXCEPT { return false; } + virtual bool supportsNonBinding() const CADET_NOEXCEPT { return true; } + virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT { return false; } + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + return true; + } - const ParamType beta0 = static_cast(_p->beta0); - const ParamType beta1 = static_cast(_p->beta1); - const ParamType water_activity = static_cast(0.1); + virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + preConsistentInitialState(t, secIdx, colPos, y, yCp, workSpace); + } - const CpStateParamType beta = beta0 * exp(beta1 * yCp[0]); + CADET_BINDINGMODELBASE_BOILERPLATE - StateParamType 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; +protected: + using ParamHandlerBindingModelBase::_paramHandler; + using ParamHandlerBindingModelBase::_reactionQuasistationarity; + using ParamHandlerBindingModelBase::_nComp; + using ParamHandlerBindingModelBase::_nBoundStates; - qSum -= y[bndIdx] / static_cast(_p->qMax[i]); + template + int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, + CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const + { + using CpStateParamType = typename DoubleActivePromoter::type; + using StateParamType = typename DoubleActivePromoter::type; + using std::pow, std::exp; - // Next bound component - ++bndIdx; - } + typename ParamHandler_t::ParamsHandle const _p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + const ParamType beta0 = static_cast(_p->beta0); + const ParamType beta1 = static_cast(_p->beta1); - bndIdx = 0; + const CpStateParamType beta = beta0 * exp(beta1 * yCp[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; + StateParamType qSum = 1.0; - const ParamType kD = static_cast(_p->kD[i]); - const ParamType kA = static_cast(_p->kA[i]); - const ParamType nu = static_cast(_p->nu[i]); + 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; - CpStateParamType bulk_water = pow(water_activity, nu * beta); - StateParamType qSumPowNu = pow(qSum, nu); + qSum -= y[bndIdx] / static_cast(_p->qMax[i]); - res[bndIdx] = kD * y[bndIdx] * bulk_water - kA * qSumPowNu * yCp[i]; + // Next bound component + ++bndIdx; + } - // Next bound component - ++bndIdx; - } + 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; - return 0; - } + const ParamType kD = static_cast(_p->kD[i]); + const ParamType kA = static_cast(_p->kA[i]); + const ParamType nu = static_cast(_p->nu[i]); - 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); + const CpStateParamType bulkWater = pow(static_cast(_p->waterActivity), nu * beta); + const StateParamType qSumPowNu = pow(qSum, nu); - auto beta0 = static_cast(_p->beta0); - auto beta1 = static_cast(_p->beta1); - auto water_activity = static_cast(0.1); + res[bndIdx] = kD * y[bndIdx] * bulkWater - kA * qSumPowNu * yCp[i]; - auto beta = static_cast(beta0 * exp(beta1 * yCp[0])); + // Next bound component + ++bndIdx; + } - double 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; + return 0; + } - qSum -= y[bndIdx] / static_cast(_p->qMax[i]); + 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); - // Next bound component - ++bndIdx; - } + const double beta0 = static_cast(_p->beta0); + const double beta1 = static_cast(_p->beta1); + const double beta = beta0 * std::exp(beta1 * yCp[0]); - bndIdx = 0; - for (int i = 0; i < _nComp; ++i) - { + double qSum = 1.0; - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBoundStates[i] == 0) - continue; + 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; - const double kD = static_cast(_p->kD[i]); - const double kA = static_cast(_p->kA[i]); - const double nu = static_cast(_p->nu[i]); + qSum -= y[bndIdx] / static_cast(_p->qMax[i]); - const double bulk_water = pow(water_activity, nu * beta); + // Next bound component + ++bndIdx; + } - // dres_i / dc_{p,0} - jac[-bndIdx - offsetCp] = static_cast(std::log(water_activity)) * beta0 * beta1 * kD * y[bndIdx] * nu * exp(beta1 * yCp[0]) * bulk_water; + 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 kD = static_cast(_p->kD[i]); + const double kA = static_cast(_p->kA[i]); + const double nu = static_cast(_p->nu[i]); - // dres_i / dc_{p,i} - jac[i - bndIdx - offsetCp] = -kA * static_cast(pow(qSum, nu)); - // 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}. + const double bulkWater = std::pow(static_cast(_p->waterActivity), nu * beta); - // Fill dres_i / dq_j - // could start j at 1, because j=1&bndIdx2=0 is the first bound component anyways - 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 / dc_{p,0} + jac[-bndIdx - offsetCp] = std::log(static_cast(_p->waterActivity)) * beta0 * beta1 * kD * y[bndIdx] * nu * std::exp(beta1 * yCp[0]) * bulkWater; - // dres_i / dq_j - jac[bndIdx2 - bndIdx] = kA * yCp[i] * nu * static_cast(pow(qSum, nu - 1)) / (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. + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -kA * std::pow(qSum, nu); + // 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}. - ++bndIdx2; - } + // Fill dres_i / dq_j + // could start j at 1, because j=1&bndIdx2=0 is the first bound component anyways + 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; - // Add to dres_i / dq_i - jac[0] += kD * bulk_water; + // dres_i / dq_j + jac[bndIdx2 - bndIdx] = kA * yCp[i] * nu * std::pow(qSum, nu - 1) / (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. - // Advance to next flux and Jacobian row - ++bndIdx; - ++jac; - } + ++bndIdx2; } - }; + // Add to dres_i / dq_i + jac[0] += kD * bulkWater; + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } +}; - typedef ConstantWaterActivityBindingBase ConstantWaterActivityBinding; - typedef ConstantWaterActivityBindingBase ExternalConstantWaterActivityBinding; +typedef ConstantWaterActivityBindingBase ConstantWaterActivityBinding; +typedef ConstantWaterActivityBindingBase ExternalConstantWaterActivityBinding; - namespace binding - { - void registerHICConstantWaterActivityModel(std::unordered_map>& bindings) - { - bindings[ConstantWaterActivityBinding::identifier()] = []() { return new ConstantWaterActivityBinding(); }; - bindings[ExternalConstantWaterActivityBinding::identifier()] = []() { return new ExternalConstantWaterActivityBinding(); }; - } - } // namespace binding +namespace binding +{ + void registerHICConstantWaterActivityModel(std::unordered_map>& bindings) + { + bindings[ConstantWaterActivityBinding::identifier()] = []() { return new ConstantWaterActivityBinding(); }; + bindings[ExternalConstantWaterActivityBinding::identifier()] = []() { return new ExternalConstantWaterActivityBinding(); }; + } +} // namespace binding - } // namespace model +} // namespace model } // namespace cadet diff --git a/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp b/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp index 835a57844..33f1a75d7 100644 --- a/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp +++ b/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp @@ -1,7 +1,7 @@ // ============================================================================= // CADET // -// Copyright © 2008-2021: The CADET Authors +// Copyright © 2008-2023: The CADET Authors // Please see the AUTHORS and CONTRIBUTORS file. // // All rights reserved. This program and the accompanying materials @@ -20,6 +20,7 @@ #include "LocalVector.hpp" #include "SimulationTypes.hpp" +#include #include #include #include @@ -54,296 +55,286 @@ namespace cadet { - namespace model +namespace model +{ + +inline const char* HICWHSParamHandler::identifier() CADET_NOEXCEPT { return "HIC_WATER_ON_HYDROPHOBIC_SURFACES"; } + +inline bool HICWHSParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("HICWHS_KA, HICWHS_KD, HICWHS_NU, and HICWHS_QMAX have to have the same size"); + + return true; +} + +inline const char* ExtHICWHSParamHandler::identifier() CADET_NOEXCEPT { return "EXT_HIC_WATER_ON_HYDROPHOBIC_SURFACES"; } + +inline bool ExtHICWHSParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +{ + if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("HICWHS_KA, HICWHS_KD, HICWHS_NU, and HICWHS_QMAX have to have the same size"); + + return true; +} + + +/** + * @brief Defines the HIC Isotherm as described by Wang et al., 2016 + * @details Implements the "water on hydrophobic surfaces" model: \f[ \begin{align} + * \beta &= \beta_0 e^{c_{p,0}\beta_1}\\ + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} + * - k_{d,i} q_i \left(\sum_j q_j \right)^{\nu_i \beta} + * \end{align} \f] + * Component @c 0 is assumed to be salt without a bound state. Multiple bound states are not supported. + * Components without bound state (i.e., salt and non-binding components) are supported. + * + * See @cite Wang2016. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ +template +class HICWHSBase : public ParamHandlerBindingModelBase +{ +public: + + HICWHSBase() { } + virtual ~HICWHSBase() CADET_NOEXCEPT { } + + static const char* identifier() { return ParamHandler_t::identifier(); } + + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) { + const bool res = BindingModelBase::configureModelDiscretization(paramProvider, nComp, nBound, boundOffset); - inline const char* HICWHSParamHandler::identifier() CADET_NOEXCEPT { return "HIC_WATER_ON_HYDROPHOBIC_SURFACES"; } + // Guarantee that salt has no bound state + if (nBound[0] != 0) + throw InvalidParameterException("HICWHS model requires exactly zero bound states for salt component"); - inline bool HICWHSParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) - { - if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) - throw InvalidParameterException("HICWHS_KA, HICWHS_KD, HICWHS_NU, and HICWHS_QMAX have to have the same size"); + return res; + } - return true; - } + virtual bool hasSalt() const CADET_NOEXCEPT { return true; } + virtual bool supportsMultistate() const CADET_NOEXCEPT { return false; } + virtual bool supportsNonBinding() const CADET_NOEXCEPT { return true; } + virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT { return false; } + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } - inline const char* ExtHICWHSParamHandler::identifier() CADET_NOEXCEPT { return "EXT_HICWHS"; } + virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + return true; + } - inline bool ExtHICWHSParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) - { - if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) - throw InvalidParameterException("HICWHS_KA, HICWHS_KD, HICWHS_NU, and HICWHS_QMAX have to have the same size"); + virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + preConsistentInitialState(t, secIdx, colPos, y, yCp, workSpace); + } - 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 + { + using CpStateParamType = typename DoubleActivePromoter::type; + using StateParamType = typename DoubleActivePromoter::type; + typename ParamHandler_t::ParamsHandle const _p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); - /** - * @brief Defines the HIC Isotherm as described by Wang et al 2016 - * @details Implements the "water on hydrophobic surfaces" model: \f[ \begin{align} - * \beta&=\beta_0 e^{c_{p,0}\beta_1}\\ - * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - * - k_{d,i} q_i \left(\sum_j q_j \right)^{\nu_i \beta} - * \end{align} \f] - * Component @c 0 is assumed to be salt without a bound state. Multiple bound states are not supported. - * Components without bound state (i.e., salt and non-binding components) are supported. - * - * See @cite Wang2016. - * @tparam ParamHandler_t Type that can add support for external function dependence - */ - template - class HICWHSBase : public ParamHandlerBindingModelBase + const ParamType beta0 = static_cast(_p->beta0); + const ParamType beta1 = static_cast(_p->beta1); + + const CpStateParamType beta = beta0 * exp(beta1 * yCp[0]); + + StateParamType freeBindingSites = 1.0; + StateType qSum = 0.0; + + // qSumOverqMax could be calculated as 1-freeBindingSites, but is calculated explicitly for clarity + StateParamType qSumOverqMax = 0.0; + + unsigned int bndIdx = 0; + for (int i = 0; i < _nComp; ++i) { - public: + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; - HICWHSBase() { } - virtual ~HICWHSBase() CADET_NOEXCEPT { } + const StateParamType qOverQmax = y[bndIdx] / static_cast(_p->qMax[i]); - static const char* identifier() { return ParamHandler_t::identifier(); } + qSum += y[bndIdx]; + qSumOverqMax += qOverQmax; + freeBindingSites -= qOverQmax; - virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) - { - const bool res = BindingModelBase::configureModelDiscretization(paramProvider, nComp, nBound, boundOffset); + // Next bound component + ++bndIdx; + } - // Guarantee that salt has no bound state - if (nBound[0] != 0) - throw InvalidParameterException("HICWHS model requires exactly zero bound states for salt component"); + // Clip free binding sites to 0.0 + if (freeBindingSites < 0) + freeBindingSites = 0; - return res; - } + 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; - virtual bool hasSalt() const CADET_NOEXCEPT { return true; } - virtual bool supportsMultistate() const CADET_NOEXCEPT { return false; } - virtual bool supportsNonBinding() const CADET_NOEXCEPT { return true; } - virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT { return false; } - virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + const ParamType qMax = static_cast(_p->qMax[i]); + const ParamType kD = static_cast(_p->kD[i]); + const ParamType kA = static_cast(_p->kA[i]); + const ParamType nu = static_cast(_p->nu[i]); - virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + if (qSum <= 0.0) { - typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); - return true; + // Use Taylor series approximation + res[bndIdx] = kA * yCp[i] * nu * qSumOverqMax - kA * yCp[i]; } - - virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + else { - preConsistentInitialState(t, secIdx, colPos, y, yCp, workSpace); + res[bndIdx] = kD * y[bndIdx] * pow(qSum, nu * beta) - kA * pow(freeBindingSites, nu) * yCp[i]; } + // Next bound component + ++bndIdx; + } - CADET_BINDINGMODELBASE_BOILERPLATE + return 0; + } - protected: - using ParamHandlerBindingModelBase::_paramHandler; - using ParamHandlerBindingModelBase::_reactionQuasistationarity; - using ParamHandlerBindingModelBase::_nComp; - using ParamHandlerBindingModelBase::_nBoundStates; + 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); - template - int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, - CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const - { - using CpStateParamType = typename DoubleActivePromoter::type; - using StateParamType = typename DoubleActivePromoter::type; + auto beta0 = static_cast(_p->beta0); + auto beta1 = static_cast(_p->beta1); - typename ParamHandler_t::ParamsHandle const _p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + auto beta = static_cast(beta0 * exp(beta1 * yCp[0])); + double freeBindingSites = 1.0; + double qSum = 0.0; - const ParamType beta0 = static_cast(_p->beta0); - const ParamType beta1 = static_cast(_p->beta1); + // qSumOverqMax could be calculated as 1-freeBindingSites, but is calculated explicitly for clarity + double qSumOverqMax = 0.0; - const CpStateParamType beta = beta0 * exp(beta1 * yCp[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; - StateParamType freeBindingSites = 1.0; - StateParamType qSum = 0.0; - // qSumOverqMax could be calculated as 1-freeBindingSites, but is calculated explicitly for clarity - StateParamType qSumOverqMax = 0.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; + const double qOverQmax = y[bndIdx] / static_cast(_p->qMax[i]); - //auto q = static_cast(y[bndIdx]); + qSum += y[bndIdx]; + qSumOverqMax += qOverQmax; + freeBindingSites -= qOverQmax; - auto qMax = static_cast(_p->qMax[i]); - qSum += y[bndIdx]; - qSumOverqMax += y[bndIdx] / qMax; - freeBindingSites -= y[bndIdx] / qMax; + // Next bound component + ++bndIdx; + } - // Next bound component - ++bndIdx; - } + // Clip free binding sites to 0.0 + if (freeBindingSites < 0) + freeBindingSites = 0; - if (freeBindingSites < 0) - { - freeBindingSites = 0; - } + 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; - bndIdx = 0; + const double kA = static_cast(_p->kA[i]); + const double kD = static_cast(_p->kD[i]); + const double nu = static_cast(_p->nu[i]); + const double qMax = static_cast(_p->qMax[i]); - 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 qMax = static_cast(_p->qMax[i]); - const ParamType kD = static_cast(_p->kD[i]); - const ParamType kA = static_cast(_p->kA[i]); - const ParamType nu = static_cast(_p->nu[i]); - - // If qSum, i.e. the term raised to nu*beta is <= zero: use Taylor series approximation - if (qSum <= 0) - { - res[bndIdx] = kA * yCp[i] * nu * qSumOverqMax - kA * yCp[i]; - } - else - { - res[bndIdx] = kD * y[bndIdx] * pow(qSum, nu * beta) - - kA * pow(freeBindingSites, nu) * yCp[i]; - } - - // Next bound component - ++bndIdx; - } + if (qSum <= 0.0) + { + // Compute the Jacobian for the Taylor series defined above + + // dres_i / dc_{p,0} + jac[-bndIdx - offsetCp] = 0; - return 0; + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = kA * nu * qSumOverqMax - kA; } - - template - void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const + else { + // todo + // dres_i / dc_{p,0} + jac[-bndIdx - offsetCp] = kD * y[bndIdx] * beta * beta1 * nu * pow(qSum, nu * beta) * std::log(qSum); - typename ParamHandler_t::ParamsHandle const _p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); - - auto beta0 = static_cast(_p->beta0); - auto beta1 = static_cast(_p->beta1); + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -kA * static_cast(pow(freeBindingSites, nu)); + } - auto beta = static_cast(beta0 * exp(beta1 * yCp[0])); + // 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 freeBindingSites = 1.0; - double qSum = 0.0; - // qSumOverqMax could be calculated as 1-freeBindingSites, but is calculated explicitly for clarity - double qSumOverqMax = 0.0; + // Fill dres_i / dq_j + 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; - unsigned int bndIdx = 0; - for (int i = 0; i < _nComp; ++i) + if (qSum <= 0) { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBoundStates[i] == 0) - continue; - - auto q_i = static_cast(y[bndIdx]); - - qSum += q_i; - qSumOverqMax += q_i / static_cast(_p->qMax[i]); - freeBindingSites -= q_i / static_cast(_p->qMax[i]); - - // Next bound component - ++bndIdx; + jac[bndIdx2 - bndIdx] = yCp[i] * kA * nu / static_cast(_p->qMax[j]); } - - if (freeBindingSites < 0) + else { - freeBindingSites = 0; + // dres_i / dq_j + jac[bndIdx2 - bndIdx] = + -kA * yCp[i] * nu * pow(freeBindingSites, nu - 1) / (-static_cast(_p->qMax[j])) + + beta * kD * nu * y[bndIdx] * pow(qSum, nu * beta - 1); } - bndIdx = 0; - for (int i = 0; i < _nComp; ++i) - { + // Getting to q_j: -bndIdx takes us to q_0, another +bndIdx2 to q_j. This means jac[bndIdx2 - bndIdx] corresponds to q_j. - // 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 nu = static_cast(_p->nu[i]); - const double qMax = static_cast(_p->qMax[i]); - - // dres_i / dc_{p,0} - if (qSum <= 0) - { - // if q below zero the jacobian is computed for the Taylor series defined above - jac[-bndIdx - offsetCp] = 0; - } - else - { - // todo - jac[-bndIdx - offsetCp] = kD * y[bndIdx] * beta * beta1 * nu * pow(qSum, nu * beta) * std::log(qSum); - } - - // dres_i / dc_{p,i} - if (qSum <= 0) - { - // if qSum is below zero the jacobian is computed for the Taylor series defined above - jac[i - bndIdx - offsetCp] = kA * nu * qSumOverqMax - kA; - } - else - { - jac[i - bndIdx - offsetCp] = -kA * static_cast(pow(freeBindingSites, nu)); - } - - // 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}. - - - // Fill dres_i / dq_j - 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; - if (qSum <= 0) - { - jac[bndIdx2 - bndIdx] = yCp[i] * kA * nu / static_cast(_p->qMax[j]); - } - else - { - // dres_i / dq_j - jac[bndIdx2 - bndIdx] = - -kA * yCp[i] * nu * pow(freeBindingSites, nu - 1) / (-static_cast(_p->qMax[j])) - + beta * kD * nu * y[bndIdx] * pow(qSum, nu * beta - 1) - ; - } - - // 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 - if (qSum <= 0) - { - // if qSum is below zero the jacobian is computed for the Taylor series defined above - jac[0] = kA * yCp[i] * nu / qMax; - } - else - { - jac[0] += kD * pow(qSum, nu * beta); - } - - // Advance to next flux and Jacobian row - ++bndIdx; - ++jac; - } + ++bndIdx2; } - }; - typedef HICWHSBase HICWHS; - typedef HICWHSBase ExternalHICWHS; - - namespace binding - { - void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map>& bindings) + // Add to dres_i / dq_i + if (qSum <= 0) { - bindings[HICWHS::identifier()] = []() { return new HICWHS(); }; - bindings[ExternalHICWHS::identifier()] = []() { return new ExternalHICWHS(); }; + // if qSum is below zero the jacobian is computed for the Taylor series defined above + jac[0] = kA * yCp[i] * nu / qMax; } - } // namespace binding + else + { + jac[0] += kD * pow(qSum, nu * beta); + } + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } +}; + +typedef HICWHSBase HICWHS; +typedef HICWHSBase ExternalHICWHS; + +namespace binding +{ + void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map>& bindings) + { + bindings[HICWHS::identifier()] = []() { return new HICWHS(); }; + bindings[ExternalHICWHS::identifier()] = []() { return new ExternalHICWHS(); }; + } +} // namespace binding - } // namespace model +} // namespace model } // namespace cadet diff --git a/test/BindingModels.cpp b/test/BindingModels.cpp index 32b1c4165..a0b57a0cb 100644 --- a/test/BindingModels.cpp +++ b/test/BindingModels.cpp @@ -1385,7 +1385,7 @@ TEST_CASE("MULTI_COMPONENT_COLLOIDAL binding model analytic Jacobian vs AD witho } } -CADET_BINDINGTEST_SINGLE("HIC_WATER_ON_HYDROPHOBIC_SURFACES", (0,1,1), (0,1,1,0), (1.0, 2.0, 3.0, 0.0, 0.0), (1.0, 2.0, 3.0, 4.0, 0.0, 0.0), \ +CADET_BINDINGTEST("HIC_WATER_ON_HYDROPHOBIC_SURFACES", "EXT_HIC_WATER_ON_HYDROPHOBIC_SURFACES", (0,1,1), (0,1,1,0), (1.0, 2.0, 3.0, 0.0, 0.0), (1.0, 2.0, 3.0, 4.0, 0.0, 0.0), \ R"json( "HICWHS_KA": [0.0, 0.872767843365959, 1.74553568673192], "HICWHS_KD": [0.0, 44.9707701873943, 89.9415403747886], "HICWHS_BETA0": 0.0184390384521496, @@ -1400,9 +1400,58 @@ CADET_BINDINGTEST_SINGLE("HIC_WATER_ON_HYDROPHOBIC_SURFACES", (0,1,1), (0,1,1,0) "HICWHS_NU": [0.0, 10, 20, 30], "HICWHS_QMAX": [0.0, 1000, 2000, 3000] )json", \ + R"json( "EXT_HICWHS_KA": [0.0, 0.0, 0.0], + "EXT_HICWHS_KA_T": [0.0, 0.872767843365959, 1.74553568673192], + "EXT_HICWHS_KA_TT": [0.0, 0.0, 0.0], + "EXT_HICWHS_KA_TTT": [0.0, 0.0, 0.0], + "EXT_HICWHS_KD": [0.0, 0.0, 0.0], + "EXT_HICWHS_KD_T": [0.0, 44.9707701873943, 89.9415403747886], + "EXT_HICWHS_KD_TT": [0.0, 0.0, 0.0], + "EXT_HICWHS_KD_TTT": [0.0, 0.0, 0.0], + "EXT_HICWHS_BETA0": 0.0, + "EXT_HICWHS_BETA0_T": 0.0184390384521496, + "EXT_HICWHS_BETA0_TT": 0.0, + "EXT_HICWHS_BETA0_TTT": 0.0, + "EXT_HICWHS_BETA1": 0.0, + "EXT_HICWHS_BETA1_T": 0.000797098469630127, + "EXT_HICWHS_BETA1_TT": 0.0, + "EXT_HICWHS_BETA1_TTT": 0.0, + "EXT_HICWHS_NU": [0.0, 0.0, 0.0], + "EXT_HICWHS_NU_T": [0.0, 10, 20], + "EXT_HICWHS_NU_TT": [0.0, 0.0, 0.0], + "EXT_HICWHS_NU_TTT": [0.0, 0.0, 0.0], + "EXT_HICWHS_QMAX": [0.0, 0.0, 0.0], + "EXT_HICWHS_QMAX_T": [0.0, 1000, 2000], + "EXT_HICWHS_QMAX_TT": [0.0, 0.0, 0.0], + "EXT_HICWHS_QMAX_TTT": [0.0, 0.0, 0.0] + )json", \ + R"json( "EXT_HICWHS_KA": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_KA_T": [0.0, 0.872767843365959, 1.74553568673192, 2.61830353009788], + "EXT_HICWHS_KA_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_KA_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_KD": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_KD_T": [0.0, 44.9707701873943, 89.9415403747886, 134.912310562183], + "EXT_HICWHS_KD_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_KD_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_BETA0": 0.0, + "EXT_HICWHS_BETA0_T": 0.0184390384521496, + "EXT_HICWHS_BETA0_TT": 0.0, + "EXT_HICWHS_BETA0_TTT": 0.0, + "EXT_HICWHS_BETA1": 0.0, + "EXT_HICWHS_BETA1_T": 0.000797098469630127, + "EXT_HICWHS_BETA1_TT": 0.0, + "EXT_HICWHS_BETA1_TTT": 0.0, + "EXT_HICWHS_NU": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_NU_T": [0.0, 10, 20, 30], + "EXT_HICWHS_NU_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_NU_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_QMAX": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_QMAX_T": [0.0, 1000, 2000, 3000], + "EXT_HICWHS_QMAX_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_HICWHS_QMAX_TTT": [0.0, 0.0, 0.0, 0.0] + )json", \ 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_UNUSED, CADET_COMPARE_BINDING_VS_NONBINDING) - CADET_BINDINGTEST("HIC_CONSTANT_WATER_ACTIVITY", "EXT_HIC_CONSTANT_WATER_ACTIVITY", (0,1,1), (0,1,1,0), (1.0, 2.0, 3.0, 1.0, 1.0), (1.0, 2.0, 3.0, 4.0, 1.0, 1.0), \ R"json( "HICCWA_KA": [0.0, 0.474361388031419, 0.948722776062837], "HICCWA_KD": [0.0, 2045.88163309217, 4091.763266184343],