Add Sips isotherm

Based on changes made by [email protected]
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 <[email protected]>

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

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
+===================  =========================  =========================================

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
+}

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`
      - ✓
      - ×

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<freundlich_ldf_model>` and the :ref:`Langmuir adsorption model<multi_component_langmuir_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 <reference_concentrations>`, :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<freundlich_ldf_model>` 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`.

src/libcadet/BindingModelFactory.cpp

@@ -45,6 +45,7 @@ namespace cadet
 			void registerBiLangmuirLDFModel(std::unordered_map<std::string, std::function<model::IBindingModel* ()>>& bindings);
 			void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings);
 			void registerHICConstantWaterActivityModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings);
+			void registerSipsModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings);
 			void registerMultiComponentLDFFreundlichModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& 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<model::SimplifiedMultiStateStericMassActionBinding>();
 	}

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,15 +291,15 @@ if (LAPACK_FOUND)
 	add_library(libcadet_static STATIC $<TARGET_OBJECTS:libcadet_object>)
 	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
 	# ---------------------------------------------------
 
 	add_library(libcadet_shared SHARED $<TARGET_OBJECTS:libcadet_object>)
 	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)