diff --git a/docs/conf.py b/docs/conf.py index cf31484c8..a95c1edbe 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -30,6 +30,7 @@ myst_enable_extensions = [ "colon_fence", "deflist", + "dollarmath" ] myst_heading_anchors = 3 diff --git a/docs/index.md b/docs/index.md index f822f916b..717e90807 100644 --- a/docs/index.md +++ b/docs/index.md @@ -120,9 +120,9 @@ help_and_feedback :caption: Python :hidden: -api/index tutorials/index recipes/index +api/index ``` ```{toctree} @@ -131,4 +131,13 @@ recipes/index rustguide/index rust_api +``` + +```{toctree} +:caption: Theory +:hidden: + +theory/eos/index +theory/dft/index +theory/models/index ``` \ No newline at end of file diff --git a/docs/theory/dft/euler_lagrange_equation.md b/docs/theory/dft/euler_lagrange_equation.md new file mode 100644 index 000000000..9d9c5f051 --- /dev/null +++ b/docs/theory/dft/euler_lagrange_equation.md @@ -0,0 +1,101 @@ +## 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,\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 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 $\mu$ and external potentials $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, $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 consideration can be described as spherical. Then the Helmholtz energy can be split into 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)\mathrm{d}r+\beta F^\mathrm{res}$$ + +with the de Broglie wavelength $\Lambda_i$. 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)\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 + +$$\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\\ +&\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)\mathrm{d}r}_{\beta\hat{F}^\mathrm{chain}} +\end{align}$$ + +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)\mathrm{d}r+\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)\mathrm{d}r+\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 +The expressions are more complex for models in which density profiles of individual segments are considered. A derivation is given in the appendix of [Rehner et al. (2022)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.105.034110). 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')\mathrm{d}r$$ + +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 in 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')\mathrm{d}r$$ + +### 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/docs/theory/dft/functional_derivatives.md b/docs/theory/dft/functional_derivatives.md new file mode 100644 index 000000000..ab1a145fc --- /dev/null +++ b/docs/theory/dft/functional_derivatives.md @@ -0,0 +1,42 @@ +## Functional derivatives + +In the last section the functional derivative + +$$\hat F_{\rho_\alpha}^\mathrm{res}(r)=\left(\frac{\delta\hat F^\mathrm{res}}{\delta\rho_\alpha(r)}\right)_{T,\rho_{\alpha'\neq\alpha}}$$ + +was introduced as part of the Euler-Lagrange equation. The implementation of these functional derivatives can be a major difficulty during the development of a new Helmholtz energy model. In $\text{FeO}_\text{s}$ it is fully automated. The core assumption is that the residual Helmholtz energy functional $\hat F^\mathrm{res}$ can be written as a sum of contributions that each can be written in the following way: + +$$F=\int f[\rho(r)]dr=\int f(\lbrace n_\gamma\rbrace)dr$$ + +The Helmholtz energy density $f$ which would in general be a functional of the density itself can be expressed as a *function* of weighted densities $n_\gamma$ which are obtained by convolving the density profiles with weight functions $\omega_\gamma^\alpha$ + +$$n_\gamma(r)=\sum_\alpha\int\rho_\alpha(r')\omega_\gamma^\alpha(r-r')dr'\tag{1}$$ + +In practice the weight functions tend to have simple shapes like step functions (i.e. the weighted density is an average over a sphere) or Dirac distributions (i.e. the weighted density is an average over the surface of a sphere). + +For Helmholtz energy functionals that can be written in this form, the calculation of the functional derivative can be automated. In general the functional derivative can be written as + +$$F_{\rho_\alpha}(r)=\int\frac{\delta f(r')}{\delta\rho_\alpha(r)}dr'=\int\sum_\gamma f_{n_\gamma}(r')\frac{\delta n_\gamma(r')}{\delta\rho_\alpha(r)}dr'$$ + +with $f_{n_\gamma}$ as abbreviation for the *partial* derivative $\frac{\partial f}{\partial n_\gamma}$. Using the definition of the weighted densities (1), the expression can be rewritten as + +$$\begin{align} +F_{\rho_\alpha}(r)&=\int\sum_\gamma f_{n_\gamma}(r')\frac{\delta n_\gamma(r')}{\delta\rho_\alpha(r)}dr'=\int\sum_\gamma f_{n_\gamma}(r')\sum_{\alpha'}\int\underbrace{\frac{\delta\rho_{\alpha'}(r'')}{\delta\rho_\alpha(r)}}_{\delta(r-r'')\delta_{\alpha\alpha'}}\omega_\gamma^{\alpha'}(r'-r'')dr''dr'\\ +&=\sum_\gamma\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r'-r)dr +\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 + +$$F_{\rho_\alpha}(r)=\sum_\gamma^\mathrm{scal}\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r-r')dr-\sum_\gamma^\mathrm{vec}\int f_{n_\gamma}(r')\omega_\gamma^\alpha(r-r')dr\tag{2}$$ + +With this distinction, the calculation of the functional derivative is split into three steps + +1. Calculate the weighted densities $n_\gamma$ from eq. (1) +2. Evaluate the partial derivatives $f_{n_\gamma}$ +3. Use eq. (2) to obtain the functional derivative $F_{\rho_\alpha}$ + +A fast method to calculate the convolution integrals required in steps 1 and 3 is shown in the next section. + +The implementation of partial derivatives by hand can be cumbersome and error prone. $\text{FeO}_\text{s}$ uses automatic differentiation with dual numbers to facilitate this step. For details about dual numbers and their generalization, we refer to [our publication](https://www.frontiersin.org/articles/10.3389/fceng.2021.758090/full). The essential relation is that the Helmholtz energy density evaluated with a dual number as input for one of the weighted densities evaluates to a dual number with the function value as real part and the partial derivative as dual part. + +$$f(n_\gamma+\varepsilon,\lbrace n_{\gamma'\neq\gamma}\rbrace)=f+f_{n_\gamma}\varepsilon$$ \ No newline at end of file diff --git a/docs/theory/dft/index.md b/docs/theory/dft/index.md new file mode 100644 index 000000000..67b79952e --- /dev/null +++ b/docs/theory/dft/index.md @@ -0,0 +1,12 @@ +# Classical density functional theory +This section explains the implementation of the core expressions from classical density functional theory in $\text{FeO}_\text{s}$. + +```{eval-rst} +.. toctree:: + :maxdepth: 1 + + euler_lagrange_equation + functional_derivatives +``` + +It is currently still under construction. You can help by [contributing](https://github.com/feos-org/feos/issues/70). \ No newline at end of file diff --git a/docs/theory/eos/index.md b/docs/theory/eos/index.md new file mode 100644 index 000000000..8112bcdd7 --- /dev/null +++ b/docs/theory/eos/index.md @@ -0,0 +1,9 @@ +# Equations of state +This section explains the thermodynamic principles and algorithms used for equation of state modeling in $\text{FeO}_\text{s}$. + +It is currently still under construction. You can help by [contributing](https://github.com/feos-org/feos/issues/70). + +```{eval-rst} +.. toctree:: + :maxdepth: 1 +``` \ No newline at end of file diff --git a/docs/theory/models/index.md b/docs/theory/models/index.md new file mode 100644 index 000000000..148f360b4 --- /dev/null +++ b/docs/theory/models/index.md @@ -0,0 +1,9 @@ +# Models +This section describes the thermodynamic models implemented in $\text{FeO}_\text{s}$. + +It is currently still under construction. You can help by [contributing](https://github.com/feos-org/feos/issues/70). + +```{eval-rst} +.. toctree:: + :maxdepth: 1 +``` \ No newline at end of file diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md index 8cf352dc0..793dc115c 100644 --- a/feos-dft/CHANGELOG.md +++ b/feos-dft/CHANGELOG.md @@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - Added `Send` and `Sync` as supertraits to `HelmholtzEnergyFunctional` and all related traits. [#57](https://github.com/feos-org/feos/pull/57) - Renamed the `3d_dft` feature to `rayon` in accordance to the other feos crates. [#62](https://github.com/feos-org/feos/pull/62) +- internally rewrote the implementation of the Euler-Lagrange equation to use a bulk density instead of the chemical potential as variable. [#60](https://github.com/feos-org/feos/pull/60) ## [0.3.2] - 2022-10-13 ### Changed diff --git a/feos-dft/src/convolver/mod.rs b/feos-dft/src/convolver/mod.rs index a558145b7..5cdaa56d2 100644 --- 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: Send + Sync { ) -> Array; } +pub(crate) struct BulkConvolver { + weight_constants: Vec>, +} + +impl> BulkConvolver { + pub(crate) fn new(weight_functions: Vec>) -> Arc> { + let weight_constants = weight_functions + .into_iter() + .map(|w| w.weight_constants(Zero::zero(), 0)) + .collect(); + Arc::new(Self { weight_constants }) + } +} + +impl> Convolver for BulkConvolver +where + Array2: Dot, Output = Array1>, +{ + fn convolve(&self, _: Array0, _: &WeightFunction) -> Array0 { + unreachable!() + } + + fn weighted_densities(&self, density: &Array1) -> Vec> { + self.weight_constants + .iter() + .map(|w| w.dot(density)) + .collect() + } + + fn functional_derivative(&self, partial_derivatives: &[Array1]) -> Array1 { + 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 index ea07db2da..58f56d246 100644 --- a/feos-dft/src/functional.rs +++ b/feos-dft/src/functional.rs @@ -445,7 +445,7 @@ impl DFT { pub fn bond_integrals( &self, temperature: f64, - functional_derivative: &Array, + exponential: &Array, convolver: &Arc>, ) -> Array where @@ -468,7 +468,6 @@ impl DFT { } } - let expdfdrho = functional_derivative.mapv(|x| (-x).exp()); let mut i_graph: Graph<_, Option>, Directed> = bond_weight_functions.map(|_, _| (), |_, _| None); @@ -495,7 +494,7 @@ impl DFT { 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, e| acc * e.weight().as_ref().unwrap(), @@ -514,7 +513,7 @@ impl DFT { } } - 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 index 0c9841790..4aa8b90ab 100644 --- a/feos-dft/src/interface/mod.rs +++ b/feos-dft/src/interface/mod.rs @@ -91,6 +91,7 @@ impl PlanarInterface { n_grid: usize, l_grid: QuantityScalar, critical_temperature: QuantityScalar, + fix_equimolar_surface: bool, ) -> EosResult { let mut profile = Self::new(vle, n_grid, l_grid)?; @@ -113,13 +114,19 @@ impl PlanarInterface { }); // specify specification - profile.profile.specification = - DFTSpecifications::total_moles_from_profile(&profile.profile)?; + if fix_equimolar_surface { + profile.profile.specification = + DFTSpecifications::total_moles_from_profile(&profile.profile)?; + } Ok(profile) } - pub fn from_pdgt(vle: &PhaseEquilibrium, 2>, n_grid: usize) -> EosResult { + pub fn from_pdgt( + vle: &PhaseEquilibrium, 2>, + n_grid: usize, + fix_equimolar_surface: bool, + ) -> EosResult { let dft = &vle.vapor().eos; if dft.component_index().len() != 1 { @@ -161,8 +168,10 @@ impl PlanarInterface { )?; // specify specification - profile.profile.specification = - DFTSpecifications::total_moles_from_profile(&profile.profile)?; + if fix_equimolar_surface { + profile.profile.specification = + DFTSpecifications::total_moles_from_profile(&profile.profile)?; + } Ok(profile) } diff --git a/feos-dft/src/interface/surface_tension_diagram.rs b/feos-dft/src/interface/surface_tension_diagram.rs index 0122432ee..4db113083 100644 --- a/feos-dft/src/interface/surface_tension_diagram.rs +++ b/feos-dft/src/interface/surface_tension_diagram.rs @@ -19,6 +19,7 @@ impl SurfaceTensionDiagram { n_grid: Option, l_grid: Option>, critical_temperature: Option>, + fix_equimolar_surface: Option, solver: Option<&DFTSolver>, ) -> Self { let n_grid = n_grid.unwrap_or(DEFAULT_GRID_POINTS); @@ -31,17 +32,19 @@ impl SurfaceTensionDiagram { 10, 100.0 * U::reference_length(), 500.0 * U::reference_temperature(), + fix_equimolar_surface.unwrap_or(false), ) } else { // initialize with pDGT for single segments and tanh for mixtures and segment DFT if vle.vapor().eos.component_index().len() == 1 { - PlanarInterface::from_pdgt(vle, n_grid) + PlanarInterface::from_pdgt(vle, n_grid, false) } else { PlanarInterface::from_tanh( vle, n_grid, l_grid.unwrap_or(100.0 * U::reference_length()), critical_temperature.unwrap_or(500.0 * U::reference_temperature()), + fix_equimolar_surface.unwrap_or(false), ) } .map(|mut profile| { diff --git a/feos-dft/src/profile.rs b/feos-dft/src/profile.rs index 4b7f09970..a44fb1b56 100644 --- a/feos-dft/src/profile.rs +++ b/feos-dft/src/profile.rs @@ -1,4 +1,4 @@ -use crate::convolver::{Convolver, ConvolverFFT}; +use crate::convolver::{BulkConvolver, Convolver, ConvolverFFT}; use crate::functional::{HelmholtzEnergyFunctional, DFT}; use crate::geometry::Grid; use crate::solver::DFTSolver; @@ -11,6 +11,7 @@ use ndarray::{ Ix1, Ix2, Ix3, RemoveAxis, }; use num_dual::Dual64; +use num_traits::Zero; use quantity::{Quantity, QuantityArray, QuantityArray1, QuantityScalar}; use std::ops::MulAssign; use std::sync::Arc; @@ -25,12 +26,11 @@ pub(crate) const CUTOFF_RADIUS: f64 = 14.0; /// for more general systems, this trait provides the possibility to declare additional /// equations for the calculation of the chemical potential during the iteration. pub trait DFTSpecification: Send + Sync { - fn calculate_chemical_potential( + fn calculate_bulk_density( &self, profile: &DFTProfile, - chemical_potential: &Array1, + bulk_density: &Array1, z: &Array1, - bulk: &State>, ) -> EosResult>; } @@ -44,11 +44,8 @@ pub enum DFTSpecifications { /// potentials are iterated together with the density profile to obtain a result /// with the specified number of particles. Moles { moles: Array1 }, - /// DFT with specified total number of moles and chemical potential differences. - TotalMoles { - total_moles: f64, - chemical_potential: Array1, - }, + /// DFT with specified total number of moles. + TotalMoles { total_moles: f64 }, } impl DFTSpecifications { @@ -78,37 +75,26 @@ impl DFTSpecifications { where ::Larger: Dimension, { - let rho = profile - .density - .to_reduced(U::reference_density())? - .sum_axis(Axis_nd(0)); - Ok(Arc::new(Self::TotalMoles { - total_moles: profile.integrate_reduced(rho), - chemical_potential: profile.reduced_chemical_potential()?, - })) + let rho = profile.density.to_reduced(U::reference_density())?; + let moles = profile.integrate_reduced_comp(&rho).sum(); + Ok(Arc::new(Self::TotalMoles { total_moles: moles })) } } impl DFTSpecification for DFTSpecifications { - fn calculate_chemical_potential( + fn calculate_bulk_density( &self, - profile: &DFTProfile, - chemical_potential: &Array1, + _profile: &DFTProfile, + bulk_density: &Array1, z: &Array1, - _: &State>, ) -> EosResult> { - let m: &Array1 = &profile.dft.m(); Ok(match self { - Self::ChemicalPotential => chemical_potential.clone(), - Self::Moles { moles } => (moles / z).mapv(f64::ln) * m, - Self::TotalMoles { - total_moles, - chemical_potential, - } => { - let exp_mu = (chemical_potential / m).mapv(f64::exp); - ((&exp_mu * *total_moles) / (z * &exp_mu).sum()).mapv(f64::ln) * m + Self::ChemicalPotential => bulk_density.clone(), + Self::Moles { moles } => moles / z, + Self::TotalMoles { total_moles } => { + bulk_density * *total_moles / (bulk_density * z).sum() } }) } @@ -305,23 +291,6 @@ where pub fn chemical_potential(&self) -> QuantityArray1 { self.bulk.chemical_potential(Contributions::Total) } - - fn reduced_chemical_potential(&self) -> EosResult> { - let temperature = self - .bulk - .temperature - .to_reduced(U::reference_temperature())?; - let lambda_de_broglie = self - .dft - .ideal_gas() - .de_broglie_wavelength(temperature, self.bulk.eos.components()); - let mu_comp = self - .chemical_potential() - .to_reduced(U::reference_molar_energy())? - / temperature - - lambda_de_broglie; - Ok(self.dft.component_index().mapv(|c| mu_comp[c])) - } } impl Clone for DFTProfile { @@ -366,49 +335,37 @@ where // Read from profile let temperature = self.temperature.to_reduced(U::reference_temperature())?; let density = self.density.to_reduced(U::reference_density())?; - let chemical_potential = self.reduced_chemical_potential()?; - let mut bulk = self.bulk.clone(); + let partial_density = self + .bulk + .partial_density + .to_reduced(U::reference_density())?; + let bulk_density = self.dft.component_index().mapv(|i| partial_density[i]); // Allocate residual vectors let mut res_rho = Array::zeros(density.raw_dim()); - let mut res_mu = Array1::zeros(chemical_potential.len()); + let mut res_bulk = Array1::zeros(bulk_density.len()); self.calculate_residual( temperature, &density, - &chemical_potential, - &mut bulk, + &bulk_density, res_rho.view_mut(), - res_mu.view_mut(), + res_bulk.view_mut(), log, )?; - Ok((res_rho, res_mu)) + Ok((res_rho, res_bulk)) } fn calculate_residual( &self, temperature: f64, density: &Array, - chemical_potential: &Array1, - bulk: &mut State>, + bulk_density: &Array1, mut res_rho: ArrayViewMut, - mut res_mu: ArrayViewMut1, + mut res_bulk: ArrayViewMut1, log: bool, ) -> EosResult<()> { - // Update bulk state - let lambda_de_broglie = self - .dft - .ideal_gas() - .de_broglie_wavelength(temperature, bulk.eos.components()); - let mut mu_comp = Array::zeros(bulk.eos.components()); - for (s, &c) in self.dft.component_index().iter().enumerate() { - mu_comp[c] = chemical_potential[s]; - } - bulk.update_chemical_potential( - &((mu_comp + lambda_de_broglie) * temperature * U::reference_molar_energy()), - )?; - // calculate intrinsic functional derivative let (_, mut dfdrho) = self.dft @@ -417,26 +374,39 @@ where // calculate total functional derivative dfdrho += &self.external_potential; + // calculate bulk functional derivative + let bulk_convolver = BulkConvolver::new(self.dft.weight_functions(temperature)); + let (_, dfdrho_bulk) = + self.dft + .functional_derivative(temperature, bulk_density, &bulk_convolver)?; + dfdrho + .outer_iter_mut() + .zip(dfdrho_bulk.into_iter()) + .zip(self.dft.m().iter()) + .for_each(|((mut df, df_b), &m)| { + df -= df_b; + df /= m + }); + // calculate bond integrals + let mut exp_dfdrho = dfdrho.mapv(|x| (-x).exp()); let bonds = self .dft - .bond_integrals(temperature, &dfdrho, &self.convolver); + .bond_integrals(temperature, &exp_dfdrho, &self.convolver); + exp_dfdrho *= &bonds; // Euler-Lagrange equation - let m = &self.dft.m(); res_rho .outer_iter_mut() - .zip(dfdrho.outer_iter()) - .zip(chemical_potential.iter()) - .zip(m.iter()) + .zip(exp_dfdrho.outer_iter()) + .zip(bulk_density.iter()) .zip(density.outer_iter()) - .zip(bonds.outer_iter()) - .for_each(|(((((mut res, df), &mu), &m), rho), is)| { + .for_each(|(((mut res, df), &rho_b), rho)| { res.assign( &(if log { - rho.mapv(f64::ln) - (mu - &df) / m - is.mapv(f64::ln) + rho.mapv(f64::ln) - rho_b.ln() - df.mapv(f64::ln) } else { - &rho - &(((mu - &df) / m).mapv(f64::exp) * is) + &rho - rho_b * &df }), ); }); @@ -451,22 +421,17 @@ where } }); - // Additional residuals for the calculation of the chemical potential - let z: Array1<_> = dfdrho - .outer_iter() - .zip(m.iter()) - .zip(bonds.outer_iter()) - .map(|((df, &m), is)| self.integrate_reduced((-&df / m).mapv(f64::exp) * is)) - .collect(); - let mu_spec = - self.specification - .calculate_chemical_potential(self, chemical_potential, &z, bulk)?; + // Additional residuals for the calculation of the bulk densitiess + let z = self.integrate_reduced_comp(&exp_dfdrho); + let bulk_spec = self + .specification + .calculate_bulk_density(self, bulk_density, &z)?; - res_mu.assign( + res_bulk.assign( &(if log { - chemical_potential - &mu_spec + bulk_density.mapv(f64::ln) - bulk_spec.mapv(f64::ln) } else { - chemical_potential.mapv(f64::exp) - mu_spec.mapv(f64::exp) + bulk_density - &bulk_spec }), ); @@ -475,43 +440,44 @@ where pub fn solve(&mut self, solver: Option<&DFTSolver>, debug: bool) -> EosResult<()> { // unwrap solver - let solver = solver.cloned().unwrap_or_default(); + let solver = solver.cloned().unwrap_or_else(|| { + if self.bulk.molefracs.iter().any(|x| x.is_zero()) { + DFTSolver::default_no_log() + } else { + DFTSolver::default() + } + }); // Read from profile + let component_index = self.dft.component_index(); let temperature = self.temperature.to_reduced(U::reference_temperature())?; let mut density = self.density.to_reduced(U::reference_density())?; - let mut chemical_potential = self.reduced_chemical_potential()?; - let mut bulk = self.bulk.clone(); + let partial_density = self + .bulk + .partial_density + .to_reduced(U::reference_density())?; + let mut bulk_density = component_index.mapv(|i| partial_density[i]); // initialize x-vector let n_rho = density.len(); - let mut x = Array1::zeros(n_rho + density.shape()[0]); + let mut x = Array1::zeros(n_rho + bulk_density.len()); x.slice_mut(s![..n_rho]) .assign(&density.view().into_shape(n_rho).unwrap()); - x.slice_mut(s![n_rho..]) - .assign(&chemical_potential.mapv(f64::exp)); + x.slice_mut(s![n_rho..]).assign(&bulk_density); // Residual function let mut residual = |x: &Array1, mut res: ArrayViewMut1, log: bool| -> EosResult<()> { // Read density and chemical potential from solution vector density.assign(&x.slice(s![..n_rho]).into_shape(density.shape()).unwrap()); - chemical_potential.assign(&x.slice(s![n_rho..]).mapv(f64::ln)); + bulk_density.assign(&x.slice(s![n_rho..])); // Create views for different residuals let (res_rho, res_mu) = res.multi_slice_mut((s![..n_rho], s![n_rho..])); let res_rho = res_rho.into_shape(density.raw_dim()).unwrap(); // Calculate residual - self.calculate_residual( - temperature, - &density, - &chemical_potential, - &mut bulk, - res_rho, - res_mu, - log, - ) + self.calculate_residual(temperature, &density, &bulk_density, res_rho, res_mu, log) }; // Call solver(s) @@ -530,7 +496,15 @@ where // Update profile self.density = density * U::reference_density(); - self.bulk = bulk; + let volume = U::reference_volume(); + let mut moles = self.bulk.moles.clone(); + bulk_density + .into_iter() + .enumerate() + .try_for_each(|(i, r)| { + moles.try_set(component_index[i], r * U::reference_density() * volume) + })?; + self.bulk = State::new_nvt(&self.bulk.eos, self.bulk.temperature, volume, &moles)?; Ok(()) } diff --git a/feos-dft/src/python/interface/mod.rs b/feos-dft/src/python/interface/mod.rs index f0787b4ab..1a3a305ec 100644 --- a/feos-dft/src/python/interface/mod.rs +++ b/feos-dft/src/python/interface/mod.rs @@ -24,24 +24,30 @@ macro_rules! impl_planar_interface { /// critical_temperature: SINumber /// An estimate for the critical temperature of the system. /// Used to guess the width of the interface. + /// fix_equimolar_surface: bool, optional + /// If True use additional constraints to fix the + /// equimolar surface of the system. + /// Defaults to False. /// /// Returns /// ------- /// PlanarInterface /// #[staticmethod] - #[pyo3(text_signature = "(vle, n_grid, l_grid, critical_temperature)")] + #[pyo3(text_signature = "(vle, n_grid, l_grid, critical_temperature, fix_equimolar_surface=None)")] fn from_tanh( vle: &PyPhaseEquilibrium, n_grid: usize, l_grid: PySINumber, critical_temperature: PySINumber, + fix_equimolar_surface: Option, ) -> PyResult { let profile = PlanarInterface::from_tanh( &vle.0, n_grid, l_grid.into(), critical_temperature.into(), + fix_equimolar_surface.unwrap_or(false), )?; Ok(PyPlanarInterface(profile)) } @@ -54,15 +60,19 @@ macro_rules! impl_planar_interface { /// The bulk phase equilibrium. /// n_grid : int /// The number of grid points. + /// fix_equimolar_surface: bool, optional + /// If True use additional constraints to fix the + /// equimolar surface of the system. + /// Defaults to False. /// /// Returns /// ------- /// PlanarInterface /// #[staticmethod] - #[pyo3(text_signature = "(vle, n_grid)")] - fn from_pdgt(vle: &PyPhaseEquilibrium, n_grid: usize) -> PyResult { - let profile = PlanarInterface::from_pdgt(&vle.0, n_grid)?; + #[pyo3(text_signature = "(vle, n_grid, fix_equimolar_surface=None)")] + fn from_pdgt(vle: &PyPhaseEquilibrium, n_grid: usize, fix_equimolar_surface: Option) -> PyResult { + let profile = PlanarInterface::from_pdgt(&vle.0, n_grid, fix_equimolar_surface.unwrap_or(false))?; Ok(PyPlanarInterface(profile)) } diff --git a/feos-dft/src/python/interface/surface_tension_diagram.rs b/feos-dft/src/python/interface/surface_tension_diagram.rs index 0f2687b8a..7350b2f67 100644 --- a/feos-dft/src/python/interface/surface_tension_diagram.rs +++ b/feos-dft/src/python/interface/surface_tension_diagram.rs @@ -19,6 +19,10 @@ macro_rules! impl_surface_tension_diagram { /// critical_temperature: SINumber, optional /// An estimate for the critical temperature, used to initialize /// density profile (default: 500 K) + /// fix_equimolar_surface: bool, optional + /// If True use additional constraints to fix the + /// equimolar surface of the system. + /// Defaults to False. /// solver: DFTSolver, optional /// Custom solver options /// @@ -27,7 +31,7 @@ macro_rules! impl_surface_tension_diagram { /// SurfaceTensionDiagram /// #[pyclass(name = "SurfaceTensionDiagram")] - #[pyo3(text_signature = "(dia, init_densities=None, n_grid=None, l_grid=None, critical_temperature=None, solver=None)")] + #[pyo3(text_signature = "(dia, init_densities=None, n_grid=None, l_grid=None, critical_temperature=None, fix_equimolar_surface=None, solver=None)")] pub struct PySurfaceTensionDiagram(SurfaceTensionDiagram); #[pymethods] @@ -39,6 +43,7 @@ macro_rules! impl_surface_tension_diagram { n_grid: Option, l_grid: Option, critical_temperature: Option, + fix_equimolar_surface: Option, solver: Option, ) -> Self { let x = dia.into_iter().map(|vle| vle.0).collect(); @@ -48,6 +53,7 @@ macro_rules! impl_surface_tension_diagram { n_grid, l_grid.map(|l| l.into()), critical_temperature.map(|c| c.into()), + fix_equimolar_surface, solver.map(|s| s.0).as_ref(), )) } diff --git a/feos-dft/src/solver.rs b/feos-dft/src/solver.rs index c7a0ee5d1..cab7e6d65 100644 --- a/feos-dft/src/solver.rs +++ b/feos-dft/src/solver.rs @@ -11,6 +11,13 @@ const DEFAULT_PARAMS_PICARD: SolverParameter = SolverParameter { tol: 1e-11, beta: 0.15, }; +const DEFAULT_PARAMS_PICARD_INIT: SolverParameter = SolverParameter { + solver: DFTAlgorithm::PicardIteration(1.0), + log: false, + max_iter: 50, + tol: 1e-5, + beta: 0.15, +}; const DEFAULT_PARAMS_ANDERSON_LOG: SolverParameter = SolverParameter { solver: DFTAlgorithm::AndersonMixing(100), log: true, @@ -66,6 +73,13 @@ impl DFTSolver { } } + pub fn default_no_log() -> Self { + Self { + parameters: vec![DEFAULT_PARAMS_PICARD_INIT, DEFAULT_PARAMS_ANDERSON], + verbosity: Verbosity::None, + } + } + /// Add a Picard iteration to the solver. pub fn picard_iteration(mut self, max_rel: Option) -> Self { let mut algorithm = DEFAULT_PARAMS_PICARD; diff --git a/src/gc_pcsaft/micelles.rs b/src/gc_pcsaft/micelles.rs index f68af8201..3c16a1690 100644 --- a/src/gc_pcsaft/micelles.rs +++ b/src/gc_pcsaft/micelles.rs @@ -32,29 +32,31 @@ pub enum MicelleSpecification { impl DFTSpecification for MicelleSpecification { - fn calculate_chemical_potential( + fn calculate_bulk_density( &self, profile: &DFTProfile, - chemical_potential: &Array1, + bulk_density: &Array1, z: &Array1, - bulk: &State>, ) -> EosResult> { Ok(match self { - Self::ChemicalPotential => chemical_potential.clone(), + Self::ChemicalPotential => bulk_density.clone(), Self::Size { delta_n_surfactant, pressure, } => { - let m: &Array1<_> = &profile.dft.m(); - let rho_s_bulk = bulk.partial_density.get(1); - let rho_w_bulk = bulk.partial_density.get(0); + let rho_s_bulk = bulk_density[1]; + let rho_w_bulk = bulk_density[0]; + let volume = U::reference_volume(); + let moles = arr1(&[rho_w_bulk, rho_s_bulk]) * U::reference_density() * volume; + let bulk = State::new_nvt(&profile.dft, profile.temperature, volume, &moles)?; let f_bulk = bulk.helmholtz_energy(Contributions::Total) / bulk.volume; - let mu_s_bulk = bulk.chemical_potential(Contributions::Total).get(1); + let mu_bulk = bulk.chemical_potential(Contributions::Total); + let mu_s_bulk = mu_bulk.get(1); + let mu_w_bulk = mu_bulk.get(0); let n_s_bulk = (rho_s_bulk * profile.volume()).to_reduced(U::reference_moles())?; - let mut spec = ((delta_n_surfactant + n_s_bulk) / z).mapv(f64::ln) * m; - spec[0] = ((pressure + f_bulk - rho_s_bulk * mu_s_bulk) / rho_w_bulk) - .to_reduced(U::reference_molar_energy())? - / bulk.temperature.to_reduced(U::reference_temperature())?; + let mut spec = (delta_n_surfactant + n_s_bulk) / z; + spec[0] = ((pressure + f_bulk - rho_s_bulk * mu_s_bulk) / mu_w_bulk) + .to_reduced(U::reference_density())?; spec } }) diff --git a/tests/gc_pcsaft/dft.rs b/tests/gc_pcsaft/dft.rs index 6cee26de1..ed6d11c97 100644 --- a/tests/gc_pcsaft/dft.rs +++ b/tests/gc_pcsaft/dft.rs @@ -193,7 +193,7 @@ fn test_dft() -> Result<(), Box> { let points = 2048; let tc = State::critical_point(&func, None, None, Default::default())?.temperature; let vle = PhaseEquilibrium::pure(&func, t, None, Default::default())?; - let profile = PlanarInterface::from_tanh(&vle, points, w, tc)?.solve(None)?; + let profile = PlanarInterface::from_tanh(&vle, points, w, tc, false)?.solve(None)?; println!( "hetero {} {} {}", profile.surface_tension.unwrap(), @@ -238,7 +238,8 @@ fn test_dft_assoc() -> Result<(), Box> { let w = 100.0 * ANGSTROM; let points = 4096; let vle = PhaseEquilibrium::pure(&func, t, None, Default::default())?; - let profile = PlanarInterface::from_tanh(&vle, points, w, 600.0 * KELVIN)?.solve(None)?; + let profile = + PlanarInterface::from_tanh(&vle, points, w, 600.0 * KELVIN, false)?.solve(None)?; println!( "hetero {} {} {}", profile.surface_tension.unwrap(), diff --git a/tests/pcsaft/dft.rs b/tests/pcsaft/dft.rs index 45ce81c55..56f5e458a 100644 --- a/tests/pcsaft/dft.rs +++ b/tests/pcsaft/dft.rs @@ -112,9 +112,10 @@ fn test_dft_propane() -> Result<(), Box> { let vle_pure = PhaseEquilibrium::pure(&func_pure, t, None, Default::default())?; let vle_full = PhaseEquilibrium::pure(&func_full, t, None, Default::default())?; let vle_full_vec = PhaseEquilibrium::pure(&func_full_vec, t, None, Default::default())?; - let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc)?.solve(None)?; - let profile_full = PlanarInterface::from_tanh(&vle_full, points, w, tc)?.solve(None)?; - let profile_full_vec = PlanarInterface::from_tanh(&vle_full_vec, points, w, tc)?.solve(None)?; + let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc, false)?.solve(None)?; + let profile_full = PlanarInterface::from_tanh(&vle_full, points, w, tc, false)?.solve(None)?; + let profile_full_vec = + PlanarInterface::from_tanh(&vle_full_vec, points, w, tc, false)?.solve(None)?; let _ = func_pure.solve_pdgt(&vle_pure, 198, 0, None)?; println!( "pure {} {} {} {}", @@ -231,8 +232,9 @@ fn test_dft_water() -> Result<(), Box> { let tc = State::critical_point(&func_pure, None, None, Default::default())?.temperature; let vle_pure = PhaseEquilibrium::pure(&func_pure, t, None, Default::default())?; let vle_full_vec = PhaseEquilibrium::pure(&func_full_vec, t, None, Default::default())?; - let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc)?.solve(None)?; - let profile_full_vec = PlanarInterface::from_tanh(&vle_full_vec, points, w, tc)?.solve(None)?; + let profile_pure = PlanarInterface::from_tanh(&vle_pure, points, w, tc, false)?.solve(None)?; + let profile_full_vec = + PlanarInterface::from_tanh(&vle_full_vec, points, w, tc, false)?.solve(None)?; println!( "pure {} {} {}", profile_pure.surface_tension.unwrap(), @@ -306,7 +308,7 @@ fn test_entropy_bulk_values() -> Result<(), Box> { )?; let func = Arc::new(PcSaftFunctional::new(Arc::new(params))); let vle = PhaseEquilibrium::pure(&func, 350.0 * KELVIN, None, Default::default())?; - let profile = PlanarInterface::from_pdgt(&vle, 2048)?.solve(None)?; + let profile = PlanarInterface::from_pdgt(&vle, 2048, false)?.solve(None)?; let s_res = profile .profile .entropy_density(Contributions::ResidualNvt)?;