Rewrite and document core DFT functionality

feos-org · prehner · Nov 11, 2022 · Oct 6, 2022 · Oct 6, 2022 · Oct 6, 2022
commit bcfe1939a8de9a4d373b066dfa7bd9e199a82a31
diff --git a/docs/theory.md b/docs/theory.md
@@ -0,0 +1,63 @@
+# Classical Density Functional Theory
+In this section, the implementation of the core expressions in classical density functional theory is explained.
+
+## 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$.
+$$\Omega(T,\bm\mu,[\bm\rho(r)])=F(T,[\bm\rho(r)])-\sum_i\int\rho_i(r)\left(\mu_i-V_i^\mathrm{ext}(r)\right)dr$$
+What makes this expression so appealing is that the intrinsic Helmholtz energy does only depend on the temperature $T$ and the density profiles $\rho_i(r)$ of the system and not on the external potential $V_i^\mathrm{ext}$.
+
+For a given temperature $T$, chemical potentials $\bm\mu$ and external potentials $\bm V^\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\tag{1}$$
+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.
+
+For a homogeneous (bulk) system, $\bm V^\mathrm{ext}=0$ and we get
+$$F_{\rho_i}^\mathrm{b}-\mu_i=0$$
+which can be inserted into (1) to give
+$$F_{\rho_i}(r)=F_{\rho_i}^\mathrm{b}-V_i^\mathrm{ext}(r)\tag{2}$$
+
+### Spherical molecules
+In the simplest case, the molecules under considerations can be described as spherical. Then the Helmholtz energy can be split in to an ideal and a residual part:
+$$\beta F=\sum_i\int\rho_i(r)\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\right)dr+\beta F^\mathrm{res}$$
+The functional derivatives for an inhomogeneous and a bulk system follow as
+$$\beta F_{\rho_i}=\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta F_{\rho_i}^\mathrm{res}$$
+$$\beta F_{\rho_i}^\mathrm{b}=\ln\left(\rho_i^\mathrm{b}\Lambda_i^3\right)+\beta F_{\rho_i}^\mathrm{b,res}$$
+Using these expressions in eq. (2) and solving for the density results in
+$$\rho_i(r)=\rho_i^\mathrm{b}e^{\beta\left(F_{\rho_i}^\mathrm{b,res}-F_{\rho_i}^\mathrm{res}(r)-V_i^\mathrm{ext}(r)\right)}$$
+which is the common form of the Euler-Lagrange equation for spherical molecules.
+
+### Homosegmented chains
+For chain molecules that do not resolve individual segments (essentially the PC-SAFT Helmholtz energy functional) a chain contribution is introduced as
+$$\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)dr$$
+Here, $m_i$ is the chain length, $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
+$$\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)dr\underbrace{-\sum_i\int\rho_i(r)\left(m_i-1\right)\left(\ln\left(y_{ii}\lambda_i(r)\Lambda_i^3\right)-1\right)dr}_{\beta\hat{F}^\mathrm{chain}}$$
+Then the total Helmholtz energy
+$$\beta F=\sum_i\int\rho_i(r)\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\right)dr+\beta F^\mathrm{chain}+\beta F^\mathrm{res}$$
+can be rearranged to
+$$\beta F=\sum_i\int\rho_i(r)m_i\left(\ln\left(\rho_i(r)\Lambda_i^3\right)-1\right)dr+\underbrace{\beta\hat{F}^\mathrm{chain}+\beta F^\mathrm{res}}_{\beta\hat{F}^\mathrm{res}}$$
+The functional derivatives are then similar to the spherical case
+$$\beta F_{\rho_i}=m_i\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta\hat{F}_{\rho_i}^\mathrm{res}$$
+$$\beta F_{\rho_i}^\mathrm{b}=m_i\ln\left(\rho_i^\mathrm{b}\Lambda_i^3\right)+\beta\hat{F}_{\rho_i}^\mathrm{b,res}$$
+and lead to a slightly modified Euler-Lagrange equation
+$$\rho_i(r)=\rho_i^\mathrm{b}e^{\frac{\beta}{m_i}\left(\hat F_{\rho_i}^\mathrm{b,res}-\hat F_{\rho_i}^\mathrm{res}(r)-V_i^\mathrm{ext}(r)\right)}$$
+
+### Heterosegmented chains
+Thex expressions are more complex for models in which density profiles of individual segments are considered. A derivation is given in the appendix of ***. The resulting Euler-Lagrange equation is given as
+$$\rho_\alpha(r)=\Lambda_i^{-3}e^{\beta\left(\mu_i-\hat F_{\rho_\alpha}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)$$
+with
+$$I_{\alpha\alpha'}(r)=\int e^{-\beta\left(F_{\rho_{\alpha'}}(r')+V_{\alpha'}^\mathrm{ext}(r')\right)}\left(\prod_{\alpha''\neq\alpha}I_{\alpha'\alpha''}(r')
+\right)\omega_\mathrm{chain}^{\alpha\alpha'}(r-r')dr$$
+The index $\alpha$ is used for every segment on component $i$, $\alpha'$ refers to all segments bonded to segment $\alpha$ and $\alpha''$ to all segments bonded to $\alpha'$. 
+For bulk systems the expressions simplify to
+$$\rho_\alpha^\mathrm{b}=\Lambda_i^{-3}e^{\beta\left(\mu_i-\sum_\gamma\hat F_{\rho_\gamma}^\mathrm{b,res}\right)}$$
+which shows that, by construction, the density of every segment on a molecule is identical in a bulk system. The index $\gamma$ refers to all segments on moecule $i$. The expressions can be combined in a similar way as for the molecular DFT:
+$$\rho_\alpha(r)=\rho_\alpha^\mathrm{b}e^{\beta\left(\sum_\gamma\hat F_{\rho_\gamma}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)$$
+At this point it can be numerically useful to redistribute the bulk contributions back into the bond integrals
+$$\rho_\alpha(r)=\rho_\alpha^\mathrm{b}e^{\beta\left(\hat F_{\rho_\alpha}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)$$
+$$I_{\alpha\alpha'}(r)=\int e^{\beta\left(\hat F_{\rho_{\alpha'}}^\mathrm{b,res}-\hat F_{\rho_{\alpha'}}^\mathrm{res}(r')-V_{\alpha'}^\mathrm{ext}(r')\right)}\left(\prod_{\alpha''\neq\alpha}I_{\alpha'\alpha''}(r')
+\right)\omega_\mathrm{chain}^{\alpha\alpha'}(r-r')dr$$
+
+### Combined expression
+To avoid having multiple implementations of the central part of the DFT code, the different descriptions of molecules can be combined in a single version of the Euler-Lagrange equation:
+$$\rho_\alpha(r)=\rho_\alpha^\mathrm{b}e^{\frac{\beta}{m_\alpha}\left(\hat F_{\rho_\alpha}^\mathrm{b,res}-\hat F_{\rho_\alpha}^\mathrm{res}(r)-V_\alpha^\mathrm{ext}(r)\right)}\prod_{\alpha'}I_{\alpha\alpha'}(r)$$
+If molecules consist of single (possibly non-spherical) segments, the Euler-Lagrange equation simplifies to that of the homosegmented chains shown above. For heterosegmented chains, the correct expression is obtained by setting $m_\alpha=1$.
diff --git a/feos-dft/src/convolver/mod.rs b/feos-dft/src/convolver/mod.rs
@@ -1,8 +1,10 @@
 use crate::geometry::{Axis, Geometry, Grid};
 use crate::weight_functions::*;
+use ndarray::linalg::Dot;
 use ndarray::prelude::*;
 use ndarray::{Axis as Axis_nd, RemoveAxis, ScalarOperand, Slice};
 use num_dual::*;
+use num_traits::Zero;
 use rustdct::DctNum;
 use std::ops::{AddAssign, MulAssign, SubAssign};
 use std::sync::Arc;
@@ -34,6 +36,45 @@ pub trait Convolver<T, D: Dimension>: Send + Sync {
     ) -> Array<T, D::Larger>;
 }
 
+pub(crate) struct BulkConvolver<T> {
+    weight_constants: Vec<Array2<T>>,
+}
+
+impl<T: DualNum<f64>> BulkConvolver<T> {
+    pub(crate) fn new(weight_functions: Vec<WeightFunctionInfo<T>>) -> Rc<dyn Convolver<T, Ix0>> {
+        let weight_constants = weight_functions
+            .into_iter()
+            .map(|w| w.weight_constants(Zero::zero(), 0))
+            .collect();
+        Rc::new(Self { weight_constants })
+    }
+}
+
+impl<T: DualNum<f64>> Convolver<T, Ix0> for BulkConvolver<T>
+where
+    Array2<T>: Dot<Array1<T>, Output = Array1<T>>,
+{
+    fn convolve(&self, _: Array0<T>, _: &WeightFunction<T>) -> Array0<T> {
+        unreachable!()
+    }
+
+    fn weighted_densities(&self, density: &Array1<T>) -> Vec<Array1<T>> {
+        self.weight_constants
+            .iter()
+            .map(|w| w.dot(density))
+            .collect()
+    }
+
+    fn functional_derivative(&self, partial_derivatives: &[Array1<T>]) -> Array1<T> {
+        self.weight_constants
+            .iter()
+            .zip(partial_derivatives.iter())
+            .map(|(w, pd)| pd.dot(w))
+            .reduce(|a, b| a + b)
+            .unwrap()
+    }
+}
+
 /// Base structure to hold either information about the weight function through
 /// `WeightFunctionInfo` or the weight functions themselves via
 /// `FFTWeightFunctions`.

diff --git a/feos-dft/src/functional.rs b/feos-dft/src/functional.rs
@@ -445,7 +445,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
     pub fn bond_integrals<D>(
         &self,
         temperature: f64,
-        functional_derivative: &Array<f64, D::Larger>,
+        exponential: &Array<f64, D::Larger>,
         convolver: &Arc<dyn Convolver<f64, D>>,
     ) -> Array<f64, D::Larger>
     where
@@ -468,7 +468,6 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
             }
         }
 
-        let expdfdrho = functional_derivative.mapv(|x| (-x).exp());
         let mut i_graph: Graph<_, Option<Array<f64, D>>, Directed> =
             bond_weight_functions.map(|_, _| (), |_, _| None);
 
@@ -495,7 +494,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
                     if edges.clone().all(|e| e.weight().is_some()) {
                         edge_id = Some(edge.id());
                         let i0 = edges.fold(
-                            expdfdrho
+                            exponential
                                 .index_axis(Axis(0), edge.target().index())
                                 .to_owned(),
                             |acc: Array<f64, D>, e| acc * e.weight().as_ref().unwrap(),
@@ -514,7 +513,7 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
             }
         }
 
-        let mut i = Array::ones(functional_derivative.raw_dim());
+        let mut i = Array::ones(exponential.raw_dim());
         for node in i_graph.node_indices() {
             for edge in i_graph.edges(node) {
                 i.index_axis_mut(Axis(0), node.index())

diff --git a/feos-dft/src/interface/mod.rs b/feos-dft/src/interface/mod.rs
@@ -112,9 +112,9 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
                     + 0.5 * (rho_l + rho_v)
             });
 
-        // specify specification
-        profile.profile.specification =
-            DFTSpecifications::total_moles_from_profile(&profile.profile)?;
+        // // specify specification
+        // profile.profile.specification =
+        //     DFTSpecifications::total_moles_from_profile(&profile.profile)?;
 
         Ok(profile)
     }
@@ -160,9 +160,9 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional> PlanarInterface<U, F> {
             r,
         )?;
 
-        // specify specification
-        profile.profile.specification =
-            DFTSpecifications::total_moles_from_profile(&profile.profile)?;
+        // // specify specification
+        // profile.profile.specification =
+        //     DFTSpecifications::total_moles_from_profile(&profile.profile)?;
 
         Ok(profile)
     }