From 6aff14ed6baf76d0abfdc56335e42c8f434a9b50 Mon Sep 17 00:00:00 2001 From: "r.jaepel" Date: Thu, 25 Apr 2024 15:38:59 +0200 Subject: [PATCH] Add MMC isotherm by Nfor 2010 --- doc/interface/binding/mmc_nfor.rst | 131 ++++++ .../hic_water_on_hydrophobic_surfaces.rst | 2 +- doc/modelling/binding/mmc_nfor.rst | 39 ++ src/libcadet/BindingModelFactory.cpp | 2 + src/libcadet/CMakeLists.txt | 1 + src/libcadet/model/binding/MMCNfor.cpp | 389 ++++++++++++++++++ test/BindingModels.cpp | 73 ++++ 7 files changed, 636 insertions(+), 1 deletion(-) create mode 100644 doc/interface/binding/mmc_nfor.rst create mode 100644 doc/modelling/binding/mmc_nfor.rst create mode 100644 src/libcadet/model/binding/MMCNfor.cpp diff --git a/doc/interface/binding/mmc_nfor.rst b/doc/interface/binding/mmc_nfor.rst new file mode 100644 index 000000000..451e09b96 --- /dev/null +++ b/doc/interface/binding/mmc_nfor.rst @@ -0,0 +1,131 @@ +.. _mmc_nfor_config: + +MMC Nfor 2010 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = MMC_NFOR** + +For information on model equations, refer to :ref:`_mmc_nfor_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 +=================== ========================= ======================= + +``MMCNFOR_KA`` + Adsorption rate constant + +**Unit:** :math:`m_{MP}^{3}~m_{SP}^{-3}~s^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + +``MMCNFOR_KD`` + Desorption rate constant + +**Unit:** :math:`s^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + +``MMCNFOR_KP`` + Protein-salt interaction factor (dependence of activity of solved protein from salt concentration) + +**Unit:** :math:`s^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + +``MMCNFOR_KS`` + Protein-protein interaction factor (dependence of activity of solved protein from protein concentration) + +**Unit:** :math:`s^{-1}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + +``MMCNFOR_NU`` + Characteristic charges of the protein; The number of sites + :math:`\nu` that the protein interacts with on the resin surface + +**Unit: [-]** + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + + +``MMCNFOR_N`` + Number of hydrophobic ligands occupied during binding + +**Unit: [-]** + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + + +``MMCNFOR_SIGMA`` + Steric factors for ionic binding of the protein; The number of sites :math:`\sigma` on + the surface that are shielded by the protein and prevented from + exchange with the salt counterions in solution + +**Unit: [-]** + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + +``MMCNFOR_S`` + Steric factor for hydrophobic binding of the protein;; The number of sites :math:`\sigma` on + the surface that are shielded by the protein and prevented from + exchange with the salt counterions in solution + +**Unit: [-]** + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** NCOMP +=================== ========================= ========================================= + + +``SMA_LAMBDA`` + Stationary phase capacity (monovalent salt counterions); The total + number of binding sites available on the resin surface + +**Unit:** :math:`mol~m_{SP}^{-3}` + +=================== ========================= ========================================= +**Type:** double **Range:** :math:`\ge 0` **Length:** 1 +=================== ========================= ========================================= + +``SMA_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 +=================== ========================= ========================================= + +``SMA_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 +=================== ========================= ========================================= diff --git a/doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst b/doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst index f14deff60..ba08f0430 100644 --- a/doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst +++ b/doc/modelling/binding/hic_water_on_hydrophobic_surfaces.rst @@ -1,4 +1,4 @@ -.. _hic_water_on_hydrophobic_surfaces_model: +.. _mmc_nfor_model: HIC Water on Hydrophobic Surfaces ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/doc/modelling/binding/mmc_nfor.rst b/doc/modelling/binding/mmc_nfor.rst new file mode 100644 index 000000000..d7bb2e9ff --- /dev/null +++ b/doc/modelling/binding/mmc_nfor.rst @@ -0,0 +1,39 @@ +.. _hic_water_on_hydrophobic_surfaces_model: + +MMC Nfor 2010 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This model implements the mixed mode chromatography isotherm described by Nfor et al 2010 :cite:`Nfor2010`. + +.. math:: + + \begin{aligned} + \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \widetilde{\gamma_i} \left[ \Lambda - \sum_j\left( \nu_j + \sigma_j \right) q_j \right]^{\nu_i}\left[ \Lambda - \sum_j\left( n_j + s_j \right) q_j \right]^{n_i} - k_{d,i} q_i c_{s}^{\nu_i} \\ + c_s &= c_{p,0}\\ + \widetilde{\gamma_i} &= e^{K_{p,i}c_{p,i} + K_{s,i}c_{s}} + \end{aligned} + +where :math:`c_{p,0}` and :math:`q_0` denote the salt concentrations in the liquid and solid phase of the beads, respectively. +The number of free binding sites + +.. math:: + + \begin{aligned} + \bar{q}_0 = \Lambda - \sum_{j=1}^{N_{\text{comp}} - 1} \left( \nu_j + \sigma_j \right) q_j = q_0 - \sum_{j=1}^{N_{\text{comp}} - 1} \sigma_j q_j + \end{aligned} + +is calculated from the number of bound counter ions :math:`q_0` by taking steric shielding into account. +In turn, the number of bound counter ions :math:`q_0` (electro-neutrality condition) is given by + +.. math:: + + \begin{aligned} + q_0 = \Lambda - \sum_{j=1}^{N_{\text{comp}} - 1} \nu_j q_j, + \end{aligned} + +which also compensates for the missing equation for :math:`\frac{\mathrm{d} q_0}{\mathrm{d}t}`. + +The concept of reference concentrations (:math:`c_{\text{ref}}` and :math:`q_{\text{ref}}`) is explained in the respective paragraph in Section :ref:`reference_concentrations`. + + +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/BindingModelFactory.cpp b/src/libcadet/BindingModelFactory.cpp index b937462ea..ca2a49403 100644 --- a/src/libcadet/BindingModelFactory.cpp +++ b/src/libcadet/BindingModelFactory.cpp @@ -44,6 +44,7 @@ namespace cadet void registerBiLangmuirLDFModel(std::unordered_map>& bindings); void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map>& bindings); void registerHICConstantWaterActivityModel(std::unordered_map>& bindings); + void registerMMCNforModel(std::unordered_map>& bindings); } } @@ -73,6 +74,7 @@ namespace cadet model::binding::registerBiLangmuirLDFModel(_bindingModels); model::binding::registerHICWaterOnHydrophobicSurfacesModel(_bindingModels); model::binding::registerHICConstantWaterActivityModel(_bindingModels); + model::binding::registerMMCNforModel(_bindingModels); registerModel(); } diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index c1c550a0b..aa98c233e 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/MMCNfor.cpp ) # LIBCADET_REACTIONMODEL_SOURCES holds all source files of reaction models diff --git a/src/libcadet/model/binding/MMCNfor.cpp b/src/libcadet/model/binding/MMCNfor.cpp new file mode 100644 index 000000000..2bad58916 --- /dev/null +++ b/src/libcadet/model/binding/MMCNfor.cpp @@ -0,0 +1,389 @@ +/*********************************************************** + * DO NOT EDIT. This file was generated by templateCodeGen from + * C:/Users/ronal/Documents/CADET-newbuild/src/libcadet/model/binding/MMCNfor.cpp + * using template + * C:/Users/ronal/Documents/CADET-newbuild/src/libcadet/model/binding/ExternalFunctionTemplate.cpp + * Please edit either the data or the applied template. + **********************************************************/ + + // ============================================================================= + // CADET + // + // Copyright © 2008-2024: The CADET Authors + // Please see the AUTHORS and CONTRIBUTORS 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/binding/BindingModelMacros.hpp" +#include "model/binding/RefConcentrationSupport.hpp" +#include "model/ModelUtils.hpp" +#include "cadet/Exceptions.hpp" +#include "model/Parameters.hpp" +#include "LocalVector.hpp" +#include "SimulationTypes.hpp" + +#include +#include +#include +#include + +/* + { + "name": "MMCNFORParamHandler", + "externalName": "ExtMMCNFORParamHandler", + "parameters": + [ + { "type": "ScalarParameter", "varName": "lambda", "confName": "MMCNFOR_LAMBDA"}, + { "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "MMCNFOR_KA"}, + { "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "MMCNFOR_KD"}, + { "type": "ScalarComponentDependentParameter", "varName": "kS", "confName": "MMCNFOR_KS"}, + { "type": "ScalarComponentDependentParameter", "varName": "kP", "confName": "MMCNFOR_KP"}, + { "type": "ScalarComponentDependentParameter", "varName": "nu", "confName": "MMCNFOR_NU"}, + { "type": "ScalarComponentDependentParameter", "varName": "n", "confName": "MMCNFOR_N"}, + { "type": "ScalarComponentDependentParameter", "varName": "sigma", "confName": "MMCNFOR_SIGMA"}, + { "type": "ScalarComponentDependentParameter", "varName": "s", "confName": "MMCNFOR_S"} + ], + "constantParameters": + [ + { "type": "ReferenceConcentrationParameter", "varName": ["refC0", "refQ"], "objName": "refConcentration", "confPrefix": "MMCNFOR_"} + ] + } +*/ + +/* Parameter description + ------------------------ + lambda = Capacity of the resin + kA = Adsorption rate + kD = Desorption rate + kS = Protein-salt interaction factor (dependence of activity of solved protein from salt concentration) + kP = Protein-protein interaction factor (dependence of activity of solved protein from protein concentration) + nu = Characteristic charge / i.e. Number of ionic ligands occupied during binding + n = Number of hydrophobic ligands occupied during binding + sigma = Steric factor for ionic binding + s = Steric factor for hydrophobic binding + refC0, refQ = Reference concentrations +*/ + +namespace cadet +{ + + namespace model + { + + inline const char* MMCNFORParamHandler::identifier() CADET_NOEXCEPT { return "MMC_NFOR"; } + + inline bool MMCNFORParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ( + (_kA.size() != _kD.size()) || + (_kA.size() != _kP.size()) || + (_kA.size() != _kS.size()) || + (_kA.size() != _nu.size()) || + (_kA.size() != _n.size()) || + (_kA.size() != _sigma.size()) || + (_kA.size() != _s.size()) || + (_kA.size() < nComp) + ) + throw InvalidParameterException("MMCNFOR_KA, MMCNFOR_KD, MMCNFOR_KP, MMCNFOR_KS, MMCNFOR_NU, MMCNFOR_N, MMCNFOR_SIGMA, and MMCNFOR_S have to have the same size"); + // Assume monovalent salt ions by default + if (_nu.get()[0] <= 0.0) + _nu.get()[0] = 1.0; + + return true; + } + + inline const char* ExtMMCNFORParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MMC_NFOR"; } + + inline bool ExtMMCNFORParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ( + (_kA.size() != _kD.size()) || + (_kA.size() != _kP.size()) || + (_kA.size() != _kS.size()) || + (_kA.size() != _nu.size()) || + (_kA.size() != _n.size()) || + (_kA.size() != _sigma.size()) || + (_kA.size() != _s.size()) || + (_kA.size() < nComp) + ) + throw InvalidParameterException("KA, KD, KP, KS, NU, N, SIGMA and S have to have the same size"); + + // Assume monovalent salt ions by default + if (_nu.base()[0] <= 0.0) + _nu.base()[0] = 1.0; + + return true; + } + + + /** + * @brief Defines the Mixed Mode Chromatrography (MMC) Isotherm as described by Nfor 2010 + * @details Implements the the MMC Nfor 2010 Isotherm: \f[ \begin{align} + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \widetilde{\gamma_i} \left[ \Lambda - \sum_j\left( \nu_j + \sigma_j \right) q_j \right]^{\nu_i}\left[ \Lambda - \sum_j\left( n_j + s_j \right) q_j \right]^{n_i} - k_{d,i} q_i c_{s}^{\nu_i} \\ + * c_s &= c_{p,0}\\ + * \widetilde{\gamma_i} &= e^{K_{p,i}c_{p,i} + K_{s,i}c_{s}} + * \end{align} \f] + * Component @c 0 is assumed to be salt. Multiple bound states are not supported. + * Components without bound state (i.e., non-binding components) are supported. + * + * See @cite Nfor2010. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ + template + class MMCNforBindingBase : public ParamHandlerBindingModelBase + { + public: + + MMCNforBindingBase() { } + virtual ~MMCNforBindingBase() 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); + + // Guarantee that salt has exactly one bound state + if (nBound[0] != 1) + throw InvalidParameterException("MMC Nfor binding model requires exactly one bound state for salt component"); + + // First flux is salt, which is always quasi-stationary + _reactionQuasistationarity[0] = true; + + return res; + } + + 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 true; } + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return false; } + + 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); + + // Compute salt component from given bound states q_j + + // Salt equation: nu_0 * q_0 - Lambda + Sum[nu_j * q_j, j] == 0 + // <=> q_0 == (Lambda - Sum[nu_j * q_j, j]) / nu_0 + y[0] = static_cast(p->lambda); + + unsigned int bndIdx = 1; + for (int j = 1; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + y[0] -= static_cast(p->nu[j]) * y[bndIdx]; + + // Next bound component + ++bndIdx; + } + + y[0] /= static_cast(p->nu[0]); + + 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); + } + + + 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; + using std::pow, std::exp; + + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Salt flux: nu_0 * q_0 - Lambda + Sum[nu_j * q_j, j] == 0 + // <=> nu_0 * q_0 == Lambda - Sum[nu_j * q_j, j] + // Also compute \bar{q}_0 = nu_0 * q_0 - Sum[sigma_j * q_j, j] + // Also compute \free_hic_binding_sites = Lambda - Sum[(s_j + n_j) * q_j, j] + res[0] = static_cast(p->nu[0]) * y[0] - static_cast(p->lambda); + StateParamType q0_bar = static_cast(p->nu[0]) * y[0]; + StateParamType free_hic_binding_sites = static_cast(p->lambda); + + unsigned int bndIdx = 1; + for (int j = 1; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + res[0] += static_cast(p->nu[j]) * y[bndIdx]; + q0_bar -= static_cast(p->sigma[j]) * y[bndIdx]; + free_hic_binding_sites -= (static_cast(p->s[j]) + static_cast(p->n[j])) * y[bndIdx]; + + // Next bound component + ++bndIdx; + } + + const ParamType refC0 = static_cast(p->refC0); + const ParamType refQ = static_cast(p->refQ); + const CpStateParamType yCp0_divRef = yCp[0] / refC0; + const StateParamType q0_bar_divRef = q0_bar / refQ; + const StateParamType free_hic_binding_sites_divRef = free_hic_binding_sites / refQ; + + // Protein fluxes: -k_{a,i} * c_{p,i} * \bar{q}_0^{nu_i / nu_0} * \free_hic_binding_sites^{n_i} + k_{d,i} * q_i * c_{p,0}^{nu_i / nu_0} + bndIdx = 1; + for (int i = 1; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + const ParamType nu_over_nu0 = static_cast(p->nu[i]) / static_cast(p->nu[0]); + const CpStateParamType c0_pow_nu = pow(yCp0_divRef, nu_over_nu0); + const StateParamType q0_bar_pow_nu = pow(q0_bar_divRef, nu_over_nu0); + const StateParamType free_hic_binding_sites_pow_n = pow(free_hic_binding_sites_divRef, static_cast(p->n[i])); + const CpStateParamType gamma = exp(static_cast(p->kP[i]) * yCp[i] + static_cast(p->kS[i]) * yCp[0]); + + // Residual + res[bndIdx] = static_cast(p->kD[i]) * y[bndIdx] * c0_pow_nu + - static_cast(p->kA[i]) * yCp[i] * q0_bar_pow_nu * free_hic_binding_sites_pow_n * gamma; + + // 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); + + double q0_bar = static_cast(p->nu[0]) * y[0]; + double free_hic_binding_sites = static_cast(p->lambda); + + //Salt flux: nu_0 * q_0 - Lambda + Sum[nu_j * q_j, j] == 0 + // <=> nu_0 * q_0 == Lambda - Sum[nu_j * q_j, j] + //Also compute \bar{q}_0 = nu_0 * q_0 - Sum[sigma_j * q_j, j] + jac[0] = static_cast(p->nu[0]); + int bndIdx = 1; + for (int j = 1; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + jac[bndIdx] = static_cast(p->nu[j]); + + // Calculate \bar{q}_0 = nu_0 * q_0 - Sum[sigma_j * q_j, j] + q0_bar -= static_cast(p->sigma[j]) * y[bndIdx]; + free_hic_binding_sites -= (static_cast(p->s[j]) + static_cast(p->n[j])) * y[bndIdx]; + + // Next bound component + ++bndIdx; + } + + // Advance to protein fluxes + ++jac; + + const double refC0 = static_cast(p->refC0); + const double refQ = static_cast(p->refQ); + const double yCp0_divRef = yCp[0] / refC0; + const double q0_bar_divRef = q0_bar / refQ; + const double free_hic_binding_sites_divRef = free_hic_binding_sites / refQ; + + // Protein fluxes: -k_{a,i} * c_{p,i} * \bar{q}_0^{nu_i/ nu_0} * free_hic_binding_sites^{n_i} + k_{d,i} * q_i * c_{p,0}^{nu_i/ nu_0} + // For the derivative calculator: k_{ai} c_{pi} e^{K_{pi}c_{pi} + K_{si}c_{s}}\left[ \nu_s*q_{s} - \sigma_j q_j \right]^{\nu_i}\left[ \lambda - \left( n_j + s_j \right) q_j \right]^{n_i} - k_{d} q_i c_{s}^{\nu_i} + // We have already computed \bar{q}_0 and free_hic_binding_sites in the loop above + bndIdx = 1; + for (int i = 1; i < _nComp; ++i) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[i] == 0) + continue; + + // Getting to c_{p,0}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0}. This means jac[-bndIdx - offsetCp] corresponds to c_{p,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}. + + const double kA = static_cast(p->kA[i]); + const double kD = static_cast(p->kD[i]); + const double kS = static_cast(p->kP[i]); + const double kP = static_cast(p->kS[i]); + const double nu = static_cast(p->nu[i]) / static_cast(p->nu[0]); + const double n = static_cast(p->n[i]); + const double gamma = std::exp(kP * yCp[i] + kS * yCp[0]); + + const double c0_pow_nu = pow(yCp0_divRef, nu); + const double q0_bar_pow_nu = pow(q0_bar_divRef, nu); + const double free_hic_binding_sites_pow_n = pow(free_hic_binding_sites_divRef, n); + const double c0_pow_nu_m1_divRef = pow(yCp0_divRef, nu - 1.0) / refC0; + const double q0_bar_pow_nu_m1_divRef = nu * pow(q0_bar_divRef, nu - 1.0) / refQ; + const double free_hic_binding_sites_pow_n_m1_divRef = n * pow(free_hic_binding_sites_divRef, n - 1.0) / refQ; + + // dres_i / dc_{p,0} + jac[-bndIdx - offsetCp] = -kS * yCp[i] * kA * q0_bar_pow_nu * free_hic_binding_sites_pow_n * gamma + + kD * y[bndIdx] * nu * c0_pow_nu_m1_divRef; + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -kA * q0_bar_pow_nu * free_hic_binding_sites_pow_n * gamma * (kP * yCp[i] + 1); + // dres_i / dq_0 + jac[-bndIdx] = -kA * yCp[i] * q0_bar_pow_nu_m1_divRef * static_cast(p->nu[0]) * gamma * free_hic_binding_sites_pow_n; + + // Fill dres_i / dq_j + int bndIdx2 = 1; + for (int j = 1; 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] = -kA * yCp[i] * gamma * q0_bar_pow_nu_m1_divRef * (-static_cast(p->sigma[j])) * free_hic_binding_sites_pow_n + - kA * yCp[i] * gamma * free_hic_binding_sites_pow_n_m1_divRef * q0_bar_pow_nu * (-static_cast(p->s[j]) - static_cast(p->n[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 * c0_pow_nu; + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } + }; + + typedef MMCNforBindingBase MMCNforBinding; + typedef MMCNforBindingBase ExternalMMCNforBinding; + + namespace binding + { + void registerMMCNforModel(std::unordered_map>& bindings) + { + bindings[MMCNforBinding::identifier()] = []() { return new MMCNforBinding(); }; + bindings[ExternalMMCNforBinding::identifier()] = []() { return new ExternalMMCNforBinding(); }; + } + } // namespace binding + + } // namespace model + +} // namespace cadet diff --git a/test/BindingModels.cpp b/test/BindingModels.cpp index 3f780d5c6..6e5feafcc 100644 --- a/test/BindingModels.cpp +++ b/test/BindingModels.cpp @@ -1520,3 +1520,76 @@ CADET_BINDINGTEST("HIC_CONSTANT_WATER_ACTIVITY", "EXT_HIC_CONSTANT_WATER_ACTIVIT )json", \ 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_USED, CADET_COMPARE_BINDING_VS_NONBINDING) + CADET_BINDINGTEST("MMC_NFOR", "EXT_HIC_MMC_NFOR", (1, 1, 1), (1, 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( "MMCNFOR_KA": [0.0, 1, 2], + "MMCNFOR_KD": [0.0, 1, 1], + "MMCNFOR_KP": [0.0, 0.01, 0.02], + "MMCNFOR_KS": [0.0, 1, 2], + "MMCNFOR_NU": [0.0, 2, 5], + "MMCNFOR_N": [0.0, 2, 5], + "MMCNFOR_QMAX": [0.0, 10, 20] + "MMCNFOR_SIGMA": [0.0, 11.83, 10.6], + "MMCNFOR_S": [0.0, 11.83, 10.6], + "MMCNFOR_LAMBDA": 100.0, + "MMCNFOR_REFC0": 1.0, + "MMCNFOR_REFQ": 1.1 +)json", \ +R"json( "MMCNFOR_KA": [0.0, 0.474361388031419, 0.948722776062837, 1.42308416409426], + "MMCNFOR_KD": [0.0, 2045.88163309217, 4091.763266184343, 6137.64489927651], + "MMCNFOR_BETA0": 0.326934960685602, + "MMCNFOR_BETA1": 0.000221461823367185, + "MMCNFOR_NU": [0.0, 10, 20, 30], + "MMCNFOR_QMAX": [0.0, 10, 20, 30] +)json", \ +R"json( "EXT_MMCNFOR_KA": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_KA_T": [0.0, 0.474361388031419, 0.948722776062837], + "EXT_MMCNFOR_KA_TT": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_KA_TTT": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_KD": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_KD_T": [0.0, 2045.88163309217, 4091.7632661843], + "EXT_MMCNFOR_KD_TT": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_KD_TTT": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_BETA0": 0.0, + "EXT_MMCNFOR_BETA0_T": 0.326934960685602, + "EXT_MMCNFOR_BETA0_TT": 0.0, + "EXT_MMCNFOR_BETA0_TTT": 0.0, + "EXT_MMCNFOR_BETA1": 0.0, + "EXT_MMCNFOR_BETA1_T": 0.000221461823367185, + "EXT_MMCNFOR_BETA1_TT": 0.0, + "EXT_MMCNFOR_BETA1_TTT": 0.0, + "EXT_MMCNFOR_NU": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_NU_T": [0.0, 10, 20], + "EXT_MMCNFOR_NU_TT": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_NU_TTT": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_QMAX": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_QMAX_T": [0.0, 10, 20], + "EXT_MMCNFOR_QMAX_TT": [0.0, 0.0, 0.0], + "EXT_MMCNFOR_QMAX_TTT": [0.0, 0.0, 0.0] +)json", \ +R"json( "EXT_MMCNFOR_KA": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_KA_T": [0.0, 0.474361388031419, 0.948722776062837, 1.42308416409426], + "EXT_MMCNFOR_KA_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_KA_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_KD": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_KD_T": [0.0, 2045.88163309217, 4091.7632661843, 6137.64489927651], + "EXT_MMCNFOR_KD_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_KD_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_BETA0": 0.0, + "EXT_MMCNFOR_BETA0_T": 0.326934960685602, + "EXT_MMCNFOR_BETA0_TT": 0.0, + "EXT_MMCNFOR_BETA0_TTT": 0.0, + "EXT_MMCNFOR_BETA1": 0.0, + "EXT_MMCNFOR_BETA1_T": 0.000221461823367185, + "EXT_MMCNFOR_BETA1_TT": 0.0, + "EXT_MMCNFOR_BETA1_TTT": 0.0, + "EXT_MMCNFOR_NU": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_NU_T": [0.0, 10, 20, 30], + "EXT_MMCNFOR_NU_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_NU_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_QMAX": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_QMAX_T": [0.0, 10, 20, 30], + "EXT_MMCNFOR_QMAX_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_MMCNFOR_QMAX_TTT": [0.0, 0.0, 0.0, 0.0] + )json", \ + 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_USED, CADET_COMPARE_BINDING_VS_NONBINDING) +