diff --git a/docs/theory/dft/derivatives.md b/docs/theory/dft/derivatives.md index 15b27c4cc..2ddffdf2f 100644 --- a/docs/theory/dft/derivatives.md +++ b/docs/theory/dft/derivatives.md @@ -1,31 +1,31 @@ -## Derivatives of density profiles +# Derivatives of density profiles For converged density properties equilibrium properties can be calculated as partial derivatives of thermodynamic potentials analogous to classical (bulk) thermodynamics. The difference is that the derivatives have to be along a path of valid density profiles (solutions of the [Euler-Lagrange equation](euler_lagrange_equation.md)). The density profiles are calculated implicitly from the Euler-Lagrange equation, which can be written simplified as -$$\Omega_{\rho_i}(T,\lbrace\mu_k\rbrace,[\lbrace\rho_k(\mathbf{r})\rbrace])=\mathcal{F}_{\rho_i}(T,[\lbrace\rho_k(\mathbf{r})\rbrace])-\mu_i+V^\mathrm{ext}(\mathbf{r})=0$$ (eqn:euler_lagrange) +$$\Omega_{\rho_i}(T,\lbrace\mu_k\rbrace,[\lbrace\rho_k(\mathbf{r})\rbrace])=F_{\rho_i}(T,[\lbrace\rho_k(\mathbf{r})\rbrace])-\mu_i+V_i^\mathrm{ext}(\mathbf{r})=0$$ (eqn:euler_lagrange) -Incorporating bond integrals can be done similar to the section on the [Newton solver](solver.md) but will not be discussed in this section. The derivatives of the density profiles can then be calculated from the total differential of eq. {eq}`eqn:euler_lagrange`. +Incorporating bond integrals can be done similar to the section on the [Newton solver](solver.md) but will not be discussed in this section. The derivatives of the density profiles can then be calculated from the total differential of eq. {eq}`eqn:euler_lagrange`, leading to $$\mathrm{d}\Omega_{\rho_i}(\mathbf{r})=\left(\frac{\partial\Omega_{\rho_i}(\mathbf{r})}{\partial T}\right)_{\mu_k,\rho_k}\mathrm{d}T+\sum_j\left(\frac{\partial\Omega_{\rho_i}(\mathbf{r})}{\partial\mu_j}\right)_{T,\mu_k,\rho_k}\mathrm{d}\mu_j+\int\sum_j\left(\frac{\delta\Omega_{\rho_i}(\mathbf{r})}{\delta\rho_j(\mathbf{r}')}\right)_{T,\mu_k,\rho_k}\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=0$$ -Using eq. {eq}`eqn:euler_lagrange` and the shortened notation for derivatives of functionals in their natural variables, e.g., $\mathcal{F}_T=\left(\frac{\partial\mathcal{F}}{\partial T}\right)_{\rho_k}$, the expression can be simplified to +Using eq. {eq}`eqn:euler_lagrange` and the shortened notation for derivatives of functionals in their natural variables, e.g., $F_T=\left(\frac{\partial F}{\partial T}\right)_{\rho_k}$, the expression can be simplified to -$$\mathcal{F}_{T\rho_i}(\mathbf{r})\mathrm{d}T-\mathrm{d}\mu_i+\int\sum_j\mathcal{F}_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=0$$ (eqn:gibbs_duhem) +$$F_{T\rho_i}(\mathbf{r})\mathrm{d}T-\mathrm{d}\mu_i+\int\sum_j F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=0$$ (eqn:gibbs_duhem) Similar to the Gibbs-Duhem relation for bulk phases, eq. {eq}`eqn:gibbs_duhem` shows how temperature, chemical potentials and the density profiles in an inhomogeneous system cannot be varied independently. The derivatives of the density profiles with respect to the intensive variables can be directly identified as -$$\int\sum_j\mathcal{F}_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial T}\right)_{\mu_k}\mathrm{d}\mathbf{r}'=-\mathcal{F}_{T\rho_i}(\mathbf{r})$$ +$$\int\sum_j F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial T}\right)_{\mu_k}\mathrm{d}\mathbf{r}'=-F_{T\rho_i}(\mathbf{r})$$ and -$$\int\sum_j\mathcal{F}_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial\mu_k}\right)_{T}\mathrm{d}\mathbf{r}'=\delta_{ik}$$ (eqn:drho_dmu) +$$\int\sum_j F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial\mu_k}\right)_{T}\mathrm{d}\mathbf{r}'=\delta_{ik}$$ (eqn:drho_dmu) Both of these expressions are implicit (linear) equations for the derivatives. They can be solved rapidly analogously to the implicit expression appearing in the [Newton solver](solver.md). In practice, it is useful to explicitly cancel out the (often unknown) thermal de Broglie wavelength $\Lambda_i$ from the expression where it has no influence. This is done by splitting the intrinsic Helmholtz energy into an ideal gas and a residual part. -$$\mathcal{F}=k_\mathrm{B}T\int\sum_im_i\rho_i(\mathbf{r})\left(\ln\left(\rho_i(\mathbf{r})\Lambda_i^3\right)-1\right)\mathrm{d}\mathbf{r}+\mathcal{\hat F}^\mathrm{res}$$ +$$F=k_\mathrm{B}T\int\sum_im_i\rho_i(\mathbf{r})\left(\ln\left(\rho_i(\mathbf{r})\Lambda_i^3\right)-1\right)\mathrm{d}\mathbf{r}+\mathcal{\hat F}^\mathrm{res}$$ -Then $\mathcal{F}_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')=m_i\frac{k_\mathrm{B}T}{\rho_i(\mathbf{r})}\delta_{ij}\delta(\mathbf{r}-\mathbf{r}')+\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')$ and eq. {eq}`eqn:drho_dmu` can be rewritten as +Then $F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')=m_i\frac{k_\mathrm{B}T}{\rho_i(\mathbf{r})}\delta_{ij}\delta(\mathbf{r}-\mathbf{r}')+\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')$ and eq. {eq}`eqn:drho_dmu` can be rewritten as $$m_i\frac{k_\mathrm{B}T}{\rho_i(\mathbf{r})}\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\mu_k}\right)_T+\int\sum_j\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial\mu_k}\right)_{T}\mathrm{d}\mathbf{r}'=\delta_{ik}$$ @@ -39,16 +39,16 @@ $$\mathrm{d}\mu_i=-s_i\mathrm{d}T+v_i\mathrm{d}p$$ which can be used in eq. {eq}`eqn:gibbs_duhem` to give -$$\left(\mathcal{F}_{T\rho_i}(\mathbf{r})+s_i\right)\mathrm{d}T+\int\sum_j\mathcal{F}_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=v_i\mathrm{d}p$$ +$$\left(F_{T\rho_i}(\mathbf{r})+s_i\right)\mathrm{d}T+\int\sum_j F_{\rho_i\rho_j}(\mathbf{r},\mathbf{r}')\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=v_i\mathrm{d}p$$ Even though $s_i$ is readily available in $\text{FeO}_\text{s}$ it is useful at this point to rewrite the partial molar entropy as -$$s_i=v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}-\mathcal{F}_{T\rho_i}^\mathrm{b}$$ +$$s_i=v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}-F_{T\rho_i}^\mathrm{b}$$ Then, the intrinsic Helmholtz energy can be split into an ideal gas and a residual part again, and the de Broglie wavelength cancels. $$\begin{align*} -&\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+\mathcal{F}_{T\rho_i}^\mathrm{res}(\mathbf{r})-\mathcal{F}_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right)\mathrm{d}T\\ +&\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+F_{T\rho_i}^\mathrm{res}(\mathbf{r})-F_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right)\mathrm{d}T\\ &~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~+m_i\frac{k_\mathrm{B}T}{\rho_i(\mathbf{r})}\delta\rho_i(\mathbf{r})+\int\sum_j\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')\delta\rho_j(\mathbf{r}')\mathrm{d}\mathbf{r}'=v_i\mathrm{d}p \end{align*}$$ @@ -60,7 +60,7 @@ and temperature $$\begin{align*} &m_i\left(\frac{\partial\rho_i(\mathbf{r})}{\partial T}\right)_{p,x_k}+\rho_i(\mathbf{r})\int\sum_j\beta\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')\left(\frac{\partial\rho_j(\mathbf{r}')}{\partial T}\right)_{p,x_k}\mathrm{d}\mathbf{r}'\\ -&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=-\frac{\rho_i(\mathbf{r})}{k_\mathrm{B}T}\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+\mathcal{F}_{T\rho_i}^\mathrm{res}(\mathbf{r})-\mathcal{F}_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right) +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=-\frac{\rho_i(\mathbf{r})}{k_\mathrm{B}T}\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+F_{T\rho_i}^\mathrm{res}(\mathbf{r})-F_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right) \end{align*}$$ follow. All derivatives $x_i$ shown here can be calculated from the same linear equation @@ -73,4 +73,4 @@ by just replacing the right hand side $y_i$. |-|-| |$\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\beta\mu_k}\right)_T$|$\rho_i(\mathbf{r})\delta_{ik}$| |$\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\beta p}\right)_{T,x_k}$|$\rho_i(\mathbf{r})v_i$| -|$\left(\frac{\partial\rho_i(\mathbf{r})}{\partial T}\right)_{p,x_k}$|$-\frac{\rho_i(\mathbf{r})}{k_\mathrm{B}T}\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+\mathcal{F}_{T\rho_i}^\mathrm{res}(\mathbf{r})-\mathcal{F}_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right)$| \ No newline at end of file +|$\left(\frac{\partial\rho_i(\mathbf{r})}{\partial T}\right)_{p,x_k}$|$-\frac{\rho_i(\mathbf{r})}{k_\mathrm{B}T}\left(m_ik_\mathrm{B}\ln\left(\frac{\rho_i(\mathbf{r})}{\rho_i^\mathrm{b}}\right)+F_{T\rho_i}^\mathrm{res}(\mathbf{r})-F_{T\rho_i}^\mathrm{b,res}+v_i\left(\frac{\partial p}{\partial T}\right)_{V,N_k}\right)$| \ No newline at end of file diff --git a/docs/theory/dft/enthalpy_of_adsorption.md b/docs/theory/dft/enthalpy_of_adsorption.md new file mode 100644 index 000000000..d111c0bd6 --- /dev/null +++ b/docs/theory/dft/enthalpy_of_adsorption.md @@ -0,0 +1,111 @@ +# Enthalpy of adsorption and the Clausius-Clapeyron relation + +## Enthalpy of adsorption +The energy balance of a simple adsorption process can be written as + +$$\frac{\mathrm{d}U}{\mathrm{d}t}=\sum_i\dot{n}_ih_i^\mathrm{b}+\dot{Q}$$ + +Here the balance is chosen to only include the fluid in the porous medium. The partial molar enthalpy $h_i^\mathrm{b}$ is the enthalpy of the fluid at the point where it enters or leaves the adsorber at which point it can be considered in a bulk state. The component balance is simply + +$$\frac{\mathrm{d}N_i}{\mathrm{d}t}=\dot{n}_i$$ + +and can be used in the energy balance to yield + +$$\dot{Q}=\frac{\mathrm{d}U}{\mathrm{d}t}-\sum_ih_i^\mathrm{b}\frac{\mathrm{d}N_i}{\mathrm{d}t}$$ + +or in differential form + +$$\delta Q=\mathrm{d}U-\sum_ih_i^\mathrm{b}\mathrm{d}N_i$$ + +The expression can be rewritten using the total differential of the internal energy (the volume of the adsorber is fixed and therefore not considered as a variable) + +$$\delta Q=\left(\frac{\partial U}{\partial T}\right)_{N_k}\mathrm{d}T-\sum_i\left(h_i^\mathrm{b}-\left(\frac{\partial U}{\partial N_i}\right)_T\right)\mathrm{d}N_i$$ + +The heat of adsorption $\delta Q$ can thus be split into a sensible part that depends on the change in temperature, and a latent part that depends on the change in loading. The expression can be simplified by using the definitions of the isochoric heat capacity $C_v=\left(\frac{\partial U}{\partial T}\right)_{N_k}$ and the **partial molar enthalpy of adsorption** + +$$\Delta h_i^\mathrm{ads}=h_i^\mathrm{b}-\left(\frac{\partial U}{\partial N_i}\right)_T$$ + +yielding + +$$\delta Q=C_v\mathrm{d}T-\sum_i\Delta h_i^\mathrm{ads}\mathrm{d}N_i$$ + +If the composition of the bulk phase is fixed, which can be a fair assumption for an adsorption process but is in general not the case for a desorption process, the heat of adsorption simplifies to + +$$\delta Q=C_v\mathrm{d}T-\Delta h^\mathrm{ads}\mathrm{d}N$$ + +with the **enthalpy of adsorption** + +$$\Delta h^\mathrm{ads}=\sum_ix_i\Delta h_i^\mathrm{ads}=h^\mathrm{b}-\sum_ix_i\left(\frac{\partial U}{\partial N_i}\right)_T$$ + +## Clausius-Clapeyron relation for porous media +The Clausius-Clapeyron relation relates the $p-T$ slope of a pure component phase transition line to the corresponding enthalpy of phase change. For a vapor-liquid phase transition, the exact relation is + +$$\frac{\mathrm{d}p^\mathrm{sat}}{\mathrm{d}T}=\frac{s^\mathrm{V}-s^\mathrm{L}}{v^\mathrm{V}-v^\mathrm{L}}=\frac{h^\mathrm{V}-h^\mathrm{L}}{T\left(v^\mathrm{V}-v^\mathrm{L}\right)}$$ (eqn:temp_dep_press) + +In this expression, the enthalpy of vaporization $\Delta h^\mathrm{vap}=h^\mathrm{V}-h^\mathrm{L}$ can be identified. The molar volumes $v$ of the two phases can be replaced by the compressibility factor $Z=\frac{pv}{RT}$. Then, eq. {eq}`eqn:temp_dep_press` simplifies to + +$$\frac{\mathrm{d}p^\mathrm{sat}}{\mathrm{d}T}=\frac{p}{R T^2}\frac{\Delta h^\mathrm{vap}}{Z^\mathrm{V}-Z^\mathrm{L}}$$ + +which can be compactly written as + +$$\frac{\mathrm{d}\ln p^\mathrm{sat}}{\mathrm{d}\frac{1}{RT}}=-\frac{\Delta h^\mathrm{vap}}{Z^\mathrm{V}-Z^\mathrm{L}}$$ (eqn:Clausius_Clapeyron_exact) + +Eq. {eq}`eqn:Clausius_Clapeyron_exact` is still an exact expression. In practice, the volume (and hence the compressibility factor) of the liquid phase can often be neglected compared to the volume of the gas phase. Additionally assuming an ideal gas phase ($Z^\mathrm{V}\approx1$), leads to the expression commonly referred to as Clausius-Clapeyron relation: + +$$\frac{\mathrm{d}\ln p^\mathrm{sat}}{\mathrm{d}\frac{1}{RT}}=-\Delta h^\mathrm{vap}$$ (eqn:Clausius_Clapeyron) + + +A similar relation can be derived for fluids adsorbed in a porous medium that is in equilibrium with a bulk phase. At this point it is important to clarify which variables describe the system +- The adsorbed fluid and the bulk phase are in equilibrium. Therefore, the temperature $T$ and chemical potentials $\mu_i$ are the same for both phases. +- The density profiles and hence the number of particles $N_i$ in the porous medium is determined by $T$ and $\mu_i$. The volume of the porous medium is not considered as a thermodynamic variable but rather as a (constant) property of the adsorbent. +- All intensive properties of the bulk phase are fully determined by $T$ and $\mu_i$. In practice it can be useful to relate these properties to measurable properties like the pressure $p$ and the composition $x_i$. + +To find an expression of the slope of an isostere (constant $N_i$), the pressure, which is only defined for the bulk phase, has to be related to properties of the adsorbed fluid. + +$$\frac{\mathrm{d}\ln p}{\mathrm{d}\frac{1}{RT}}=-\frac{RT^2}{p}\frac{\mathrm{d}p}{\mathrm{d}T}$$ + +First, the pressure can be replaced using the Gibbs-Duhem relation for the bulk phase (index $\mathrm{b}$) + +$$\frac{\mathrm{d}\ln p}{\mathrm{d}\frac{1}{RT}}=-\frac{RT^2}{pv^\mathrm{b}}\left(s^\mathrm{b}+\sum_ix_i\left(\frac{\partial\mu_i}{\partial T}\right)_{N_k}\right)$$ (eqn:clausius_clapeyron_intermediate) + +Here the directional derivative $\frac{\mathrm{d}\mu_i}{\mathrm{d}T}$ could be replaced with a partial derivative amongst the variables describing the adsorbed fluid. The partial derivative can then be replaced using a Maxwell relation based on the Helmholtz energy $F$ as follows + +$$\left(\frac{\partial\mu_i}{\partial T}\right)_{N_k}=\left(\frac{\partial^2 F}{\partial T\partial N_i}\right)=-\left(\frac{\partial S}{\partial N_i}\right)_T$$ + +Using the Maxwell relation together with the compressibility factor of the bulk phase $Z^\mathrm{b}=\frac{pv^\mathrm{b}}{RT}$ in eq. {eq}`eqn:clausius_clapeyron_intermediate` results in + +$$\frac{\mathrm{d}\ln p}{\mathrm{d}\frac{1}{RT}}=-\frac{T}{Z^\mathrm{b}}\left(s^\mathrm{b}-\sum_ix_i\left(\frac{\partial S}{\partial N_i}\right)_T\right)$$ + +Finally, using $h^\mathrm{b}=Ts^\mathrm{b}+\sum_ix_i\mu_i$ and $\mathrm{d}U=T\mathrm{d}S+\sum_i\mu_i\mathrm{d}N_i$ leads to + +$$\frac{\mathrm{d}\ln p}{\mathrm{d}\frac{1}{RT}}=-\frac{1}{Z^\mathrm{b}}\left(h^\mathrm{b}-\sum_ix_i\left(\frac{\partial U}{\partial N_i}\right)_T\right)=-\frac{\Delta h^\mathrm{ads}}{Z^\mathrm{b}}$$ (eqn:deriv_relation_hads) + +The relation is exact and valid for an arbitrary number of components in the fluid phase. + + +## Calculation of the enthalpy of adsorption from classical DFT +In a DFT context, the introduction of entropies and internal energies are just unnecessary complications. The most useful definition of the (partial molar) enthalpy of adsorption is + +$$\Delta h_i^\mathrm{ads}=T\left(s_i^\mathrm{b}+\left(\frac{\partial\mu_i}{\partial T}\right)_{N_k}\right)$$ + +The derivative at constant number of particles is problematic and has to be replaced. This is done starting from the total differential of the number of particles + +$$\mathrm{d}N_i=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\mathrm{d}\mu_j+\left(\frac{\partial N_i}{\partial T}\right)_{\mu_k}\mathrm{d}T$$ (eqn:dn) + +Calculating the derivative with respect to $T$ at constant $N_i$ leads to + +$$0=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\left(\frac{\partial\mu_j}{\partial T}\right)_{N_k}+\left(\frac{\partial N_i}{\partial T}\right)_{\mu_k}$$ (eqn:dndt_1) + +from which the unknown derivative $\left(\frac{\partial\mu_i}{\partial T}\right)_{N_k}$ can be calculated. In practice the expression has the disadvantage that $\left(\frac{\partial N_i}{\partial T}\right)_{\mu_k}$ depends on the (sometimes unknown) thermal de Broglie wavelength which cancels later with $s_i^\mathrm{b}$. This can be remedied by first calculating the derivative of eq. {eq}`eqn:dn` with respect to $T$ at constant (bulk) pressure and composition. + +$$\left(\frac{\partial N_i}{\partial T}\right)_{p,x_k}=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\left(\frac{\partial\mu_j}{\partial T}\right)_{p,x_k}+\left(\frac{\partial N_i}{\partial T}\right)_{\mu_k}$$ (eqn:dndt_2) + +From classical bulk thermodynamics we know $\left(\frac{\partial\mu_j}{\partial T}\right)_{p,x_k}=-s_j^\mathrm{b}$ and therefore, eq. {eq}`eqn:dndt_2` can be used in eq. {eq}`eqn:dndt_1` to give + +$$0=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\left(s_j^\mathrm{b}+\left(\frac{\partial\mu_j}{\partial T}\right)_{N_k}\right)+\left(\frac{\partial N_i}{\partial T}\right)_{p,x_k}$$ + +After multiplying with $T$, the following elegant expression remains + +$$0=\sum_j\left(\frac{\partial N_i}{\partial\mu_j}\right)_T\Delta h_j^\mathrm{ads}+T\left(\frac{\partial N_i}{\partial T}\right)_{p,x_k}$$ + +which is a symmetric linear system of equations due to $\left(\frac{\partial N_i}{\partial\mu_j}\right)_T=-\left(\frac{\partial^2\Omega}{\partial\mu_i\partial\mu_j}\right)_T$. The derivatives of the particle numbers are obtained by integrating over the respective derivatives of the density profiles which were discussed [previously](derivatives.md). \ No newline at end of file diff --git a/docs/theory/dft/euler_lagrange_equation.md b/docs/theory/dft/euler_lagrange_equation.md index 81291ec72..02902551c 100644 --- a/docs/theory/dft/euler_lagrange_equation.md +++ b/docs/theory/dft/euler_lagrange_equation.md @@ -1,17 +1,17 @@ # Euler-Lagrange equation -The fundamental expression in classical density functional theory is the relation between the grand potential $\Omega$ and the intrinsic Helmholtz energy $F$. +The fundamental expression in classical density functional theory is the relation between the grand potential functional $\Omega$ and the intrinsic Helmholtz energy functional $F$. $$\Omega(T,\mu,[\rho(r)])=F(T,[\rho(r)])-\sum_i\int\rho_i(r)\left(\mu_i-V_i^\mathrm{ext}(r)\right)\mathrm{d}r$$ -What makes this expression so appealing is that the intrinsic Helmholtz energy only depends on the temperature $T$ and the density profiles $\rho_i(r)$ of the system and not on the external potential $V_i^\mathrm{ext}(r)$. +What makes this expression so appealing is that the intrinsic Helmholtz energy functional only depends on the temperature $T$ and the density profiles $\rho_i(r)$ of the system and not on the external potentials $V_i^\mathrm{ext}(r)$. -For a given temperature $T$, chemical potentials $\mu$ and external potentials $V^\mathrm{ext}(r)$ the grand potential reaches a minimum at equilibrium. Mathematically this condition can be written as +For a given temperature $T$, chemical potentials $\mu_i$ and external potentials $V_i^\mathrm{ext}(r)$ the grand potential reaches a minimum at equilibrium. Mathematically this condition can be written as $$\left.\frac{\delta\Omega}{\delta\rho_i(r)}\right|_{T,\mu}=F_{\rho_i}(r)-\mu_i+V_i^{\mathrm{ext}}(r)=0$$ (eqn:euler_lagrange_mu) -where $F_{\rho_i}(r)=\left.\frac{\delta F}{\delta\rho_i(r)}\right|_T$ is short for the functional derivative of the intrinsic Helmholtz energy. In this context, eq. (1) is commonly referred to as the Euler-Lagrange equation, an implicit nonlinear integral equation which needs to be solved for the density profiles of the system. +where $F_{\rho_i}(r)=\left.\frac{\delta F}{\delta\rho_i(r)}\right|_T$ is short for the functional derivative of the intrinsic Helmholtz energy. In this context, eq. (1) is commonly referred to as the Euler-Lagrange equation, an implicit nonlinear integral equation which needs to be solved for the equilibrium density profiles of the system. -For a homogeneous (bulk) system, $V^\mathrm{ext}=0$ and we get +For a homogeneous (bulk) system, $V_i^\mathrm{ext}=0$ and we get $$F_{\rho_i}^\mathrm{b}-\mu_i=0$$ (eqn:euler_lagrange_bulk) @@ -50,8 +50,8 @@ For chain molecules that do not resolve individual segments (essentially the PC- $$\beta F^\mathrm{chain}=-\sum_i\int\rho_i(r)\left(m_i-1\right)\ln\left(\frac{y_{ii}\lambda_i(r)}{\rho_i(r)}\right)\mathrm{d}r$$ -Here, $m_i$ is the number of segments (i.e., the PC-SAFT chain length parameter), $y_{ii}$ the cavity correlation function at contact in the reference fluid, and $\lambda_i$ a weighted density. -The presence of $\rho(r)$ in the logarithm poses numerical problems. Therefore, it is convenient to rearrange the expression as +Where $m_i$ is the number of segments (i.e., the PC-SAFT chain length parameter), $y_{ii}$ is the cavity correlation function at contact in the reference fluid, and $\lambda_i$ is a weighted density. +The presence of $\rho_i(r)$ in the logarithm poses numerical problems. Therefore, it is convenient to rearrange the expression as $$\begin{align} \beta F^\mathrm{chain}=&\sum_i\int\rho_i(r)\left(m_i-1\right)\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\right)\mathrm{d}r\\ diff --git a/docs/theory/dft/functional_derivatives.md b/docs/theory/dft/functional_derivatives.md index 66fbda6f8..79817e4ab 100644 --- a/docs/theory/dft/functional_derivatives.md +++ b/docs/theory/dft/functional_derivatives.md @@ -25,7 +25,7 @@ F_{\rho_\alpha}(r)&=\int\sum_\gamma f_{n_\gamma}(r')\frac{\delta n_\gamma(r')}{\ &=\sum_\gamma\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r'-r)\mathrm{d}r \end{align}$$ -At this point the parity of the weight functions has to be taken into account. By construction scalar and spherically symmetric weight functions (the standard case) are even functions, i.e., $\omega(-r)=\omega(r)$. In contrast, vector valued weight functions, as they appear in fundamental measure theory, have odd parity, i.e., $\omega(-r)=-\omega(r)$. Therefore, the sum over the weight functions needs to be split into two contributions, as +At this point the parity of the weight functions has to be taken into account. By construction scalar and spherically symmetric weight functions (the standard case) are even functions, i.e., $\omega(-r)=\omega(r)$. In contrast, vector valued weight functions, as they appear in fundamental measure theory, are odd functions, i.e., $\omega(-r)=-\omega(r)$. Therefore, the sum over the weight functions needs to be split into two contributions, as $$F_{\rho_\alpha}(r)=\sum_\gamma^\mathrm{scal}\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r-r')\mathrm{d}r-\sum_\gamma^\mathrm{vec}\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r-r')\mathrm{d}r\tag{2}$$ diff --git a/docs/theory/dft/index.md b/docs/theory/dft/index.md index 1cd980279..3fb26446d 100644 --- a/docs/theory/dft/index.md +++ b/docs/theory/dft/index.md @@ -9,6 +9,7 @@ This section explains the implementation of the core expressions from classical functional_derivatives solver derivatives + enthalpy_of_adsorption pdgt ``` diff --git a/docs/theory/dft/solver.md b/docs/theory/dft/solver.md index e59450c12..45ee55ab9 100644 --- a/docs/theory/dft/solver.md +++ b/docs/theory/dft/solver.md @@ -52,7 +52,7 @@ The linear integral equation has to be solved for the step $\Delta\rho(r)$. Expl $$\int\sum_\beta\frac{\delta\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)}{\delta\rho_\beta(r')}\Delta\rho_\beta(r')\mathrm{d}r'\approx\frac{\mathcal{F}_\alpha\left(r;\left[\rho(r)+s\Delta\rho(r)\right]\right)-\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)}{s}$$ -However this approach requires the choice of an appropriate step size $s$ (something that we want to get away from in $\text{FeO}_\text{s}$) and also an evaluation of the full residual in every step of the linear solver. The solver can be sped up by doing parts of the functional derivative analytically beforehand. Using the definition of the residual in the rhs of eq. {eq}`eqn:newton` leads to +However this approach requires the choice of an appropriate step size $s$ (something that we want to avoid in $\text{FeO}_\text{s}$) and also an evaluation of the full residual in every step of the linear solver. The solver can be sped up by doing parts of the functional derivative analytically beforehand. Using the definition of the residual in the rhs of eq. {eq}`eqn:newton` leads to $$\begin{align*} q_\alpha(r)&\equiv-\int\sum_\beta\frac{\delta\mathcal{F}_\alpha\left(r;\left[\rho(r)\right]\right)}{\delta\rho_\beta(r')}\Delta\rho_\beta(r')\mathrm{d}r'\\ diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md index c09c28a3e..4fc58a151 100644 --- a/feos-dft/CHANGELOG.md +++ b/feos-dft/CHANGELOG.md @@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] ### Added - Added new methods `drho_dmu`, `drho_dp` and `drho_dt` that calculate partial derivatives of density profiles to every DFT profile. Also includes direct access to the integrated derivatives `dn_dmu`, `dn_dp` and `dn_dt`. [#134](https://github.com/feos-org/feos/pull/134) +- Added `enthalpy_of_adsoprtion` and `partial_molar_enthalpy_of_adsorption` getters to pores and adsorption isotherms. [#135](https://github.com/feos-org/feos/pull/135) ## [0.4.0] - 2023-01-27 ### Added diff --git a/feos-dft/src/adsorption/mod.rs b/feos-dft/src/adsorption/mod.rs index 921a30b2a..3f517a079 100644 --- a/feos-dft/src/adsorption/mod.rs +++ b/feos-dft/src/adsorption/mod.rs @@ -444,4 +444,39 @@ where } }) } + + pub fn partial_molar_enthalpy_of_adsorption(&self) -> SIArray2 { + let h_ads: Vec<_> = self + .profiles + .iter() + .map(|p| { + match p + .as_ref() + .ok() + .and_then(|p| p.partial_molar_enthalpy_of_adsorption().ok()) + { + Some(p) => p, + None => { + f64::NAN * Array1::ones(self.components) * SIUnit::reference_molar_energy() + } + } + }) + .collect(); + SIArray2::from_shape_fn((self.components, self.profiles.len()), |(j, i)| { + h_ads[i].get(j) + }) + } + + pub fn enthalpy_of_adsorption(&self) -> SIArray1 { + SIArray1::from_shape_fn(self.profiles.len(), |i| { + match self.profiles[i] + .as_ref() + .ok() + .and_then(|p| p.enthalpy_of_adsorption().ok()) + { + Some(p) => p, + None => f64::NAN * SIUnit::reference_molar_energy(), + } + }) + } } diff --git a/feos-dft/src/adsorption/pore.rs b/feos-dft/src/adsorption/pore.rs index ac25fd97f..5f3ac4501 100644 --- a/feos-dft/src/adsorption/pore.rs +++ b/feos-dft/src/adsorption/pore.rs @@ -9,7 +9,8 @@ use feos_core::{Contributions, EosResult, EosUnit, State, StateBuilder}; use ndarray::prelude::*; use ndarray::Axis as Axis_nd; use ndarray::RemoveAxis; -use quantity::si::{SIArray, SIArray2, SINumber, SIUnit}; +use num_dual::linalg::LU; +use quantity::si::{SIArray, SIArray1, SIArray2, SINumber, SIUnit}; use std::sync::Arc; const POTENTIAL_OFFSET: f64 = 2.0; @@ -128,6 +129,20 @@ where self.interfacial_tension = None; self } + + pub fn partial_molar_enthalpy_of_adsorption(&self) -> EosResult { + let a = self.profile.dn_dmu()?; + let a_unit = a.get((0, 0)); + let b = -self.profile.temperature * self.profile.dn_dt()?; + let b_unit = b.get(0); + + let h_ads = LU::new(a.to_reduced(a_unit)?)?.solve(&b.to_reduced(b_unit)?); + Ok(h_ads * b_unit / a_unit) + } + + pub fn enthalpy_of_adsorption(&self) -> EosResult { + Ok((self.partial_molar_enthalpy_of_adsorption()? * &self.profile.bulk.molefracs).sum()) + } } impl PoreSpecification for Pore1D { diff --git a/feos-dft/src/python/adsorption/mod.rs b/feos-dft/src/python/adsorption/mod.rs index 225acca10..99970cb55 100644 --- a/feos-dft/src/python/adsorption/mod.rs +++ b/feos-dft/src/python/adsorption/mod.rs @@ -248,6 +248,16 @@ macro_rules! impl_adsorption_isotherm { fn get_grand_potential(&mut self) -> PySIArray1 { self.0.grand_potential().into() } + + #[getter] + fn get_partial_molar_enthalpy_of_adsorption(&self) -> PySIArray2 { + self.0.partial_molar_enthalpy_of_adsorption().into() + } + + #[getter] + fn get_enthalpy_of_adsorption(&self) -> PySIArray1 { + self.0.enthalpy_of_adsorption().into() + } } }; } diff --git a/feos-dft/src/python/adsorption/pore.rs b/feos-dft/src/python/adsorption/pore.rs index bab576c6e..0bb8ef39f 100644 --- a/feos-dft/src/python/adsorption/pore.rs +++ b/feos-dft/src/python/adsorption/pore.rs @@ -122,6 +122,16 @@ macro_rules! impl_pore { fn get_interfacial_tension(&self) -> Option { self.0.interfacial_tension.map(PySINumber::from) } + + #[getter] + fn get_partial_molar_enthalpy_of_adsorption(&self) -> PyResult { + Ok(self.0.partial_molar_enthalpy_of_adsorption()?.into()) + } + + #[getter] + fn get_enthalpy_of_adsorption(&self) -> PyResult { + Ok(self.0.enthalpy_of_adsorption()?.into()) + } } /// Parameters required to specify a 3D pore. @@ -227,6 +237,16 @@ macro_rules! impl_pore { fn get_interfacial_tension(&self) -> Option { self.0.interfacial_tension.map(PySINumber::from) } + + #[getter] + fn get_partial_molar_enthalpy_of_adsorption(&self) -> PyResult { + Ok(self.0.partial_molar_enthalpy_of_adsorption()?.into()) + } + + #[getter] + fn get_enthalpy_of_adsorption(&self) -> PyResult { + Ok(self.0.enthalpy_of_adsorption()?.into()) + } } }; }