diff --git a/docs/theory/dft/derivatives.md b/docs/theory/dft/derivatives.md new file mode 100644 index 000000000..15b27c4cc --- /dev/null +++ b/docs/theory/dft/derivatives.md @@ -0,0 +1,76 @@ +## 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) + +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`. + +$$\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 + +$$\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) + +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})$$ + +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) + +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}$$ + +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 + +$$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}$$ + +In practice, the division by the density should be avoided for numerical reasons and the energetic properties are reduced with the factor $\beta=\frac{1}{k_\mathrm{B}T}$. The final expression is + +$$m_i\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\beta\mu_k}\right)_T+\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\beta\mu_k}\right)_{T}\mathrm{d}\mathbf{r}'=\rho_i(\mathbf{r})\delta_{ik}$$ + +For the temperature derivative, it is more convenient to express eq. {eq}`eqn:gibbs_duhem` in terms of the pressure of a bulk phase that is in equilibrium with the inhomogeneous system. In the following, only paths along **constant bulk composition** are considered. With this constraint, the total differential of the chemical potential simplifies to + +$$\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$$ + +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}$$ + +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\\ +&~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~+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*}$$ + +Finally the expressions for the derivatives with respect to pressure + +$$m_i\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\beta p}\right)_{T,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\beta p}\right)_{T,x_k}\mathrm{d}\mathbf{r}'=v_i\rho_i(\mathbf{r})$$ + +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) +\end{align*}$$ + +follow. All derivatives $x_i$ shown here can be calculated from the same linear equation + +$$m_ix_i+\rho_i(\mathbf{r})\int\sum_j\beta\mathcal{\hat F}_{\rho_i\rho_j}^\mathrm{res}(\mathbf{r},\mathbf{r}')x_i\mathrm{d}\mathbf{r}'=y_i$$ + +by just replacing the right hand side $y_i$. + +|derivative|right hand side| +|-|-| +|$\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 diff --git a/docs/theory/dft/index.md b/docs/theory/dft/index.md index ab255620e..1cd980279 100644 --- a/docs/theory/dft/index.md +++ b/docs/theory/dft/index.md @@ -8,6 +8,7 @@ This section explains the implementation of the core expressions from classical euler_lagrange_equation functional_derivatives solver + derivatives pdgt ``` diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md index e331e2b23..c09c28a3e 100644 --- a/feos-dft/CHANGELOG.md +++ b/feos-dft/CHANGELOG.md @@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [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) ## [0.4.0] - 2023-01-27 ### Added diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs index 0a19d8a19..2c037060a 100644 --- a/feos-dft/src/functional.rs +++ b/feos-dft/src/functional.rs @@ -441,6 +441,43 @@ impl DFT { )) } + #[allow(clippy::type_complexity)] + pub(crate) fn functional_derivative_dual( + &self, + temperature: f64, + density: &Array, + convolver: &Arc>, + ) -> EosResult<(Array, Array)> + where + D: Dimension, + D::Larger: Dimension, + { + let temperature_dual = Dual64::from(temperature).derive(); + let density_dual = density.mapv(Dual64::from); + let weighted_densities = convolver.weighted_densities(&density_dual); + let contributions = self.contributions(); + let mut partial_derivatives = Vec::with_capacity(contributions.len()); + let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0))); + for (c, wd) in contributions.iter().zip(weighted_densities) { + let nwd = wd.shape()[0]; + let ngrid = wd.len() / nwd; + let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0))); + let mut pd = Array::zeros(wd.raw_dim()); + c.first_partial_derivatives_dual( + temperature_dual, + wd.into_shape((nwd, ngrid)).unwrap(), + phi.view_mut().into_shape(ngrid).unwrap(), + pd.view_mut().into_shape((nwd, ngrid)).unwrap(), + )?; + partial_derivatives.push(pd); + helmholtz_energy_density += φ + } + Ok(( + helmholtz_energy_density, + convolver.functional_derivative(&partial_derivatives) * temperature_dual, + )) + } + /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$ pub fn bond_integrals( &self, diff --git a/feos-dft/src/functional_contribution.rs b/feos-dft/src/functional_contribution.rs index b458d386b..67541246d 100644 --- a/feos-dft/src/functional_contribution.rs +++ b/feos-dft/src/functional_contribution.rs @@ -3,7 +3,7 @@ use feos_core::{EosResult, HelmholtzEnergyDual, StateHD}; use ndarray::prelude::*; use ndarray::RemoveAxis; use num_dual::*; -use num_traits::Zero; +use num_traits::{One, Zero}; use std::fmt::Display; macro_rules! impl_helmholtz_energy { @@ -75,6 +75,7 @@ pub trait FunctionalContributionDual>: Display { pub trait FunctionalContribution: FunctionalContributionDual + FunctionalContributionDual + + FunctionalContributionDual> + FunctionalContributionDual, f64>> + FunctionalContributionDual + FunctionalContributionDual @@ -114,6 +115,31 @@ pub trait FunctionalContribution: Ok(()) } + fn first_partial_derivatives_dual( + &self, + temperature: Dual64, + weighted_densities: Array2, + mut helmholtz_energy_density: ArrayViewMut1, + mut first_partial_derivative: ArrayViewMut2, + ) -> EosResult<()> { + let mut wd = weighted_densities.mapv(Dual::from_re); + let t = Dual::from_re(temperature); + let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0))); + + for i in 0..wd.shape()[0] { + wd.index_axis_mut(Axis(0), i) + .map_inplace(|x| x.eps[0] = Dual64::one()); + phi = self.calculate_helmholtz_energy_density(t, wd.view())?; + first_partial_derivative + .index_axis_mut(Axis(0), i) + .assign(&phi.mapv(|p| p.eps[0])); + wd.index_axis_mut(Axis(0), i) + .map_inplace(|x| x.eps[0] = Dual64::zero()); + } + helmholtz_energy_density.assign(&phi.mapv(|p| p.re)); + Ok(()) + } + fn second_partial_derivatives( &self, temperature: f64, @@ -161,6 +187,7 @@ pub trait FunctionalContribution: impl FunctionalContribution for T where T: FunctionalContributionDual + FunctionalContributionDual + + FunctionalContributionDual> + FunctionalContributionDual, f64>> + FunctionalContributionDual + FunctionalContributionDual diff --git a/feos-dft/src/geometry.rs b/feos-dft/src/geometry.rs index e28c292ef..a596b79c1 100644 --- a/feos-dft/src/geometry.rs +++ b/feos-dft/src/geometry.rs @@ -149,11 +149,9 @@ impl Axis { } let x0 = 0.5 * ((-alpha * points as f64).exp() + (-alpha * (points - 1) as f64).exp()); let grid = (0..points) - .into_iter() .map(|i| l * x0 * (alpha * i as f64).exp()) .collect(); let edges = (0..=points) - .into_iter() .map(|i| { if i == 0 { 0.0 @@ -166,7 +164,6 @@ impl Axis { let k0 = (2.0 * alpha).exp() * (2.0 * alpha.exp() + (2.0 * alpha).exp() - 1.0) / ((1.0 + alpha.exp()).powi(2) * ((2.0 * alpha).exp() - 1.0)); let integration_weights = (0..points) - .into_iter() .map(|i| { (match i { 0 => k0 * (2.0 * alpha).exp(), diff --git a/feos-dft/src/profile.rs b/feos-dft/src/profile.rs index fcf4b3b33..86ab4479e 100644 --- a/feos-dft/src/profile.rs +++ b/feos-dft/src/profile.rs @@ -3,14 +3,13 @@ use crate::functional::{HelmholtzEnergyFunctional, DFT}; use crate::geometry::Grid; use crate::solver::{DFTSolver, DFTSolverLog}; use crate::weight_functions::WeightFunctionInfo; -use feos_core::{Contributions, EosError, EosResult, EosUnit, EquationOfState, State}; +use feos_core::{Contributions, EosError, EosResult, EosUnit, EquationOfState, State, Verbosity}; use ndarray::{ Array, Array1, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix1, Ix2, Ix3, RemoveAxis, }; use num_dual::Dual64; -use quantity::si::{SIArray, SIArray1, SINumber, SIUnit}; +use quantity::si::{SIArray, SIArray1, SIArray2, SINumber, SIUnit}; use quantity::Quantity; -// use quantity::{Quantity, QuantityArray, SIArray1, SINumber}; use std::ops::MulAssign; use std::sync::Arc; @@ -269,21 +268,22 @@ where }) } - /// Return the number of moles of each component in the system. - pub fn moles(&self) -> SIArray1 { - let rho = self - .density - .to_reduced(SIUnit::reference_density()) - .unwrap(); - let mut d = rho.raw_dim(); - d[0] = self.dft.components(); - let mut density_comps = Array::zeros(d); + /// Integrate each segment individually and aggregate to components. + pub fn integrate_segments>( + &self, + profile: &Quantity, SIUnit>, + ) -> SIArray1 { + let integral = self.integrate_comp(profile); + let mut integral_comp = Array1::zeros(self.dft.components()) * integral.get(0); for (i, &j) in self.dft.component_index().iter().enumerate() { - density_comps - .index_axis_mut(Axis_nd(0), j) - .assign(&rho.index_axis(Axis_nd(0), i)); + integral_comp.try_set(j, integral.get(i)).unwrap(); } - self.integrate_comp(&(density_comps * SIUnit::reference_density())) + integral_comp + } + + /// Return the number of moles of each component in the system. + pub fn moles(&self) -> SIArray1 { + self.integrate_segments(&self.density) } /// Return the total number of moles in the system. @@ -531,4 +531,133 @@ where )? * SIUnit::reference_pressure(); Ok(self.integrate(&internal_energy_density)) } + + fn density_derivative(&self, lhs: &Array) -> EosResult> { + let rho = self.density.to_reduced(SIUnit::reference_density())?; + let partial_density = self + .bulk + .partial_density + .to_reduced(SIUnit::reference_density())?; + let rho_bulk = self.dft.component_index().mapv(|i| partial_density[i]); + + let second_partial_derivatives = self.second_partial_derivatives(&rho)?; + let (_, _, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, &rho_bulk, false)?; + + let rhs = |x: &_| { + let delta_functional_derivative = + self.delta_functional_derivative(x, &second_partial_derivatives); + let mut xm = x.clone(); + xm.outer_iter_mut() + .zip(self.dft.m().iter()) + .for_each(|(mut x, &m)| x *= m); + let delta_i = self.delta_bond_integrals(&exp_dfdrho, &delta_functional_derivative); + xm + (delta_functional_derivative - delta_i) * &rho + }; + let mut log = DFTSolverLog::new(Verbosity::None); + Self::gmres(rhs, lhs, 200, 1e-13, &mut log) + } + + /// Return the partial derivatives of the density profiles w.r.t. the chemical potentials $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial\mu_k}\right)_T$ + pub fn drho_dmu(&self) -> EosResult::Larger>> { + let shape = self.density.shape(); + let shape: Vec<_> = std::iter::once(&shape[0]).chain(shape).copied().collect(); + let mut drho_dmu = Array::zeros(shape).into_dimensionality().unwrap(); + for (k, mut d) in drho_dmu.outer_iter_mut().enumerate() { + let mut lhs = self.density.to_reduced(SIUnit::reference_density())?; + for (i, mut l) in lhs.outer_iter_mut().enumerate() { + if i != k { + l.fill(0.0); + } + } + d.assign(&self.density_derivative(&lhs)?); + } + Ok(drho_dmu + * (SIUnit::reference_density() / SIUnit::reference_molar_entropy() / self.temperature)) + } + + /// Return the partial derivatives of the number of moles w.r.t. the chemical potentials $\left(\frac{\partial N_i}{\partial\mu_k}\right)_T$ + pub fn dn_dmu(&self) -> EosResult { + let drho_dmu = self.drho_dmu()?; + let n = drho_dmu.shape()[0]; + let dn_dmu = SIArray2::from_shape_fn([n; 2], |(i, j)| { + self.integrate(&drho_dmu.index_axis(Axis_nd(0), i).index_axis(Axis_nd(0), j)) + }); + Ok(dn_dmu) + } + + /// Return the partial derivatives of the density profiles w.r.t. the bulk pressure at constant temperature and bulk composition $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial p}\right)_{T,\mathbf{x}}$ + pub fn drho_dp(&self) -> EosResult> { + let mut lhs = self.density.to_reduced(SIUnit::reference_density())?; + let v = self + .bulk + .partial_molar_volume(Contributions::Total) + .to_reduced(SIUnit::reference_volume() / SIUnit::reference_moles())?; + for (mut l, &c) in lhs.outer_iter_mut().zip(self.dft.component_index().iter()) { + l *= v[c]; + } + self.density_derivative(&lhs) + .map(|x| x / (SIUnit::reference_molar_entropy() * self.temperature)) + } + + /// Return the partial derivatives of the number of moles w.r.t. the bulk pressure at constant temperature and bulk composition $\left(\frac{\partial N_i}{\partial p}\right)_{T,\mathbf{x}}$ + pub fn dn_dp(&self) -> EosResult { + Ok(self.integrate_segments(&self.drho_dp()?)) + } + + /// Return the partial derivatives of the density profiles w.r.t. the temperature at constant bulk pressure and composition $\left(\frac{\partial\rho_i(\mathbf{r})}{\partial T}\right)_{p,\mathbf{x}}$ + /// + /// Not compatible with heterosegmented DFT. + pub fn drho_dt(&self) -> EosResult> { + let rho = self.density.to_reduced(SIUnit::reference_density())?; + let t = self + .temperature + .to_reduced(SIUnit::reference_temperature())?; + + // calculate temperature derivative of functional derivative + let functional_contributions = self.dft.contributions(); + let weight_functions: Vec> = functional_contributions + .iter() + .map(|c| c.weight_functions(Dual64::from(t).derive())) + .collect(); + let convolver: Arc> = + ConvolverFFT::plan(&self.grid, &weight_functions, None); + let (_, dfdrhodt) = self.dft.functional_derivative_dual(t, &rho, &convolver)?; + + // calculate temperature derivative of bulk functional derivative + let partial_density = self + .bulk + .partial_density + .to_reduced(SIUnit::reference_density())?; + let rho_bulk = self.dft.component_index().mapv(|i| partial_density[i]); + let bulk_convolver = BulkConvolver::new(weight_functions); + let (_, dfdrhodt_bulk) = + self.dft + .functional_derivative_dual(t, &rho_bulk, &bulk_convolver)?; + + // solve for drho_dt + let x = (self.bulk.partial_molar_volume(Contributions::Total) + * self.bulk.dp_dt(Contributions::Total)) + .to_reduced(SIUnit::reference_molar_entropy())?; + let mut lhs = dfdrhodt.mapv(|d| d.eps[0]); + lhs.outer_iter_mut() + .zip(dfdrhodt_bulk.into_iter()) + .zip(x.into_iter()) + .for_each(|((mut lhs, d), x)| lhs -= d.eps[0] - x); + lhs.outer_iter_mut() + .zip(rho.outer_iter()) + .zip(rho_bulk.into_iter()) + .zip(self.dft.m().iter()) + .for_each(|(((mut lhs, rho), rho_b), &m)| lhs += &((&rho / rho_b).mapv(f64::ln) * m)); + + lhs *= &(-&rho / t); + self.density_derivative(&lhs) + .map(|x| x * (SIUnit::reference_density() / SIUnit::reference_temperature())) + } + + /// Return the partial derivatives of the number of moles w.r.t. the temperature at constant bulk pressure and composition $\left(\frac{\partial N_i}{\partial T}\right)_{p,\mathbf{x}}$ + /// + /// Not compatible with heterosegmented DFT. + pub fn dn_dt(&self) -> EosResult { + Ok(self.integrate_segments(&self.drho_dt()?)) + } } diff --git a/feos-dft/src/python/profile.rs b/feos-dft/src/python/profile.rs index 36cd8a00d..4223f4dab 100644 --- a/feos-dft/src/python/profile.rs +++ b/feos-dft/src/python/profile.rs @@ -1,6 +1,6 @@ #[macro_export] macro_rules! impl_profile { - ($struct:ident, $arr:ident, $arr2:ident, $si_arr:ident, $si_arr2:ident, $py_arr2:ident, [$([$ind:expr, $ax:ident]),+]) => { + ($struct:ident, $arr:ident, $arr2:ident, $si_arr:ident, $si_arr2:ident, $py_arr2:ident, [$([$ind:expr, $ax:ident]),+]$(, $si_arr3:ident)?) => { #[pymethods] impl $struct { /// Calculate the residual for the given profile. @@ -178,6 +178,37 @@ macro_rules! impl_profile { self.0.profile.grand_potential_density()?, )) } + $( + #[getter] + fn get_drho_dmu(&self) -> PyResult<$si_arr3> { + Ok(($si_arr3::from(self.0.profile.drho_dmu()?))) + } + )? + + #[getter] + fn get_dn_dmu(&self) -> PyResult { + Ok((PySIArray2::from(self.0.profile.dn_dmu()?))) + } + + #[getter] + fn get_drho_dp(&self) -> PyResult<$si_arr2> { + Ok(($si_arr2::from(self.0.profile.drho_dp()?))) + } + + #[getter] + fn get_dn_dp(&self) -> PyResult { + Ok((PySIArray1::from(self.0.profile.dn_dp()?))) + } + + #[getter] + fn get_drho_dt(&self) -> PyResult<$si_arr2> { + Ok(($si_arr2::from(self.0.profile.drho_dt()?))) + } + + #[getter] + fn get_dn_dt(&self) -> PyResult { + Ok((PySIArray1::from(self.0.profile.dn_dt()?))) + } } }; } @@ -192,7 +223,8 @@ macro_rules! impl_1d_profile { PySIArray1, PySIArray2, PyArray2, - [$([0, $ax]),+] + [$([0, $ax]),+], + PySIArray3 ); }; } diff --git a/feos-dft/src/solver.rs b/feos-dft/src/solver.rs index 0fbce2cce..3da93b1f8 100644 --- a/feos-dft/src/solver.rs +++ b/feos-dft/src/solver.rs @@ -150,6 +150,7 @@ impl DFTSolver { } } +/// A log that stores the residuals and execution time of DFT solvers. #[derive(Clone)] pub struct DFTSolverLog { verbosity: Verbosity, @@ -160,7 +161,7 @@ pub struct DFTSolverLog { } impl DFTSolverLog { - fn new(verbosity: Verbosity) -> Self { + pub(crate) fn new(verbosity: Verbosity) -> Self { log_iter!( verbosity, "solver | iter | time | residual " @@ -474,7 +475,7 @@ where Ok((false, newton.max_iter)) } - fn gmres( + pub(crate) fn gmres( rhs: R, r0: &Array, max_iter: usize, @@ -547,7 +548,7 @@ where Ok(x) } - fn second_partial_derivatives( + pub(crate) fn second_partial_derivatives( &self, density: &Array, ) -> EosResult::Larger>>> { @@ -577,7 +578,7 @@ where Ok(second_partial_derivatives) } - fn delta_functional_derivative( + pub(crate) fn delta_functional_derivative( &self, delta_density: &Array, second_partial_derivatives: &[Array::Larger>], @@ -607,7 +608,7 @@ where .functional_derivative(&delta_partial_derivatives) } - fn delta_bond_integrals( + pub(crate) fn delta_bond_integrals( &self, exponential: &Array, delta_functional_derivative: &Array, @@ -770,7 +771,7 @@ impl DFTSolver { newton.tol, ), }; - res += &format!("\n|{}|{}|{:e}|", solver, max_iter, tol); + res += &format!("\n|{solver}|{max_iter}|{tol:e}|"); } res }