cleanup and documentation

feos-org · prehner · Feb 24, 2023 · Feb 14, 2023 · Feb 15, 2023 · Feb 15, 2023
commit a5e42f8fcda33f6b28f29840d3bcd4d6a4f73863
diff --git a/docs/theory/dft/derivatives.md 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)$|
diff --git 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
@@ -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`.
 
 ## [0.4.0] - 2023-01-27
 ### Added

diff --git a/feos-dft/src/profile.rs b/feos-dft/src/profile.rs
@@ -268,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<S: Data<Elem = f64>>(
+        &self,
+        profile: &Quantity<ArrayBase<S, D::Larger>, 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,38 +532,50 @@ where
         Ok(self.integrate(&internal_energy_density))
     }
 
-    pub fn drho_dmu(&self) -> EosResult<SIArray<<D::Larger as Dimension>::Larger>> {
+    fn density_derivative(&self, lhs: &Array<f64, D::Larger>) -> EosResult<Array<f64, D::Larger>> {
         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 shape = rho.shape();
-        let shape: Vec<_> = std::iter::once(&shape[0]).chain(shape).copied().collect();
-        let mut drho_dmu = Array::zeros(shape).into_dimensionality().unwrap();
         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);
-            // x + (delta_functional_derivative - delta_i) * rho
-            xm + delta_functional_derivative * &rho
+            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<SIArray<<D::Larger as Dimension>::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 = rho.clone();
+            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::gmres(rhs, &lhs, 200, 1e-13, &mut log)?);
+            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<SIArray2> {
         let drho_dmu = self.drho_dmu()?;
         let n = drho_dmu.shape()[0];
@@ -572,71 +585,34 @@ where
         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<SIArray<D::Larger>> {
-        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 lhs = rho.clone();
+        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];
         }
-        let mut log = DFTSolverLog::new(Verbosity::None);
-        Self::gmres(rhs, &lhs, 200, 1e-13, &mut log)
+        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<SIArray1> {
-        let dn_dp = self.integrate_comp(&self.drho_dp()?);
-        let mut dn_dp_comp = Array1::zeros(self.dft.components()) * SIUnit::reference_moles()
-            / SIUnit::reference_pressure();
-        for (i, &j) in self.dft.component_index().iter().enumerate() {
-            dn_dp_comp.try_set(j, dn_dp.get(i)).unwrap();
-        }
-        Ok(dn_dp_comp)
+        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<SIArray<D::Larger>> {
         let rho = self.density.to_reduced(SIUnit::reference_density())?;
         let t = self
             .temperature
             .to_reduced(SIUnit::reference_temperature())?;
 
-        let second_partial_derivatives = self.second_partial_derivatives(&rho)?;
-
-        // define rhs
-        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) * t
-            (xm + delta_functional_derivative * &rho) * t
-        };
-
         // calculate temperature derivative of functional derivative
         let functional_contributions = self.dft.contributions();
         let weight_functions: Vec<WeightFunctionInfo<Dual64>> = functional_contributions
@@ -659,27 +635,29 @@ where
                 .functional_derivative_dual(t, &rho_bulk, &bulk_convolver)?;
 
         // solve for drho_dt
-        let x = (self.bulk.dp_dni(Contributions::Total) * self.bulk.dp_dt(Contributions::Total)
-            / self.bulk.dp_dv(Contributions::Total))
+        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);
+            .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);
-        let mut log = DFTSolverLog::new(Verbosity::None);
-        let drho_dt = Self::gmres(rhs, &lhs, 200, 1e-13, &mut log)?;
-        Ok(drho_dt * (SIUnit::reference_density() / SIUnit::reference_temperature()))
+        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<SIArray1> {
-        Ok(self.integrate_comp(&self.drho_dt()?))
+        Ok(self.integrate_segments(&self.drho_dt()?))
     }
 }
diff --git 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,
-Original file line number
+Diff line change
@@ Expand Up @@
        euler_lagrange_equation
        functional_derivatives
        solver
+       derivatives
        pdgt
     ```
@@ Expand Down @@