feos-org · prehner · Feb 24, 2023 · Feb 14, 2023 · Feb 15, 2023 · Feb 15, 2023
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`. [#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
@@ -441,6 +441,43 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
         ))
     }
 
+    #[allow(clippy::type_complexity)]
+    pub(crate) fn functional_derivative_dual<D>(
+        &self,
+        temperature: f64,
+        density: &Array<f64, D::Larger>,
+        convolver: &Arc<dyn Convolver<Dual64, D>>,
+    ) -> EosResult<(Array<Dual64, D>, Array<Dual64, D::Larger>)>
+    where
+        D: Dimension,
+        D::Larger: Dimension<Smaller = D>,
+    {
+        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 += &phi;
+        }
+        Ok((
+            helmholtz_energy_density,
+            convolver.functional_derivative(&partial_derivatives) * temperature_dual,
+        ))
+    }
+
     /// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$
     pub fn bond_integrals<D>(
         &self,

diff --git 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<N: DualNum<f64>>: Display {
 pub trait FunctionalContribution:
     FunctionalContributionDual<f64>
     + FunctionalContributionDual<Dual64>
+    + FunctionalContributionDual<Dual<Dual64, f64>>
     + FunctionalContributionDual<Dual<DualVec64<3>, f64>>
     + FunctionalContributionDual<HyperDual64>
     + FunctionalContributionDual<Dual2_64>
@@ -114,6 +115,31 @@ pub trait FunctionalContribution:
         Ok(())
     }
 
+    fn first_partial_derivatives_dual(
+        &self,
+        temperature: Dual64,
+        weighted_densities: Array2<Dual64>,
+        mut helmholtz_energy_density: ArrayViewMut1<Dual64>,
+        mut first_partial_derivative: ArrayViewMut2<Dual64>,
+    ) -> 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<T> FunctionalContribution for T where
     T: FunctionalContributionDual<f64>
         + FunctionalContributionDual<Dual64>
+        + FunctionalContributionDual<Dual<Dual64, f64>>
         + FunctionalContributionDual<Dual<DualVec64<3>, f64>>
         + FunctionalContributionDual<HyperDual64>
         + FunctionalContributionDual<Dual2_64>

diff --git 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(),
-Original file line number
+Diff line change
@@ Expand Up @@
        euler_lagrange_equation
        functional_derivatives
        solver
+       derivatives
        pdgt
     ```
@@ Expand Down @@