From bcfe1939a8de9a4d373b066dfa7bd9e199a82a31 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 6 Oct 2022 16:13:09 +0200 Subject: [PATCH 01/11] Rewrite and document core DFT functionality --- docs/theory.md | 63 +++++++++++ feos-dft/src/convolver/mod.rs | 41 +++++++ feos-dft/src/functional.rs | 7 +- feos-dft/src/interface/mod.rs | 12 +-- feos-dft/src/profile.rs | 195 +++++++++++++++------------------- src/gc_pcsaft/micelles.rs | 26 ++--- 6 files changed, 214 insertions(+), 130 deletions(-) create mode 100644 docs/theory.md diff --git a/docs/theory.md b/docs/theory.md new file mode 100644 index 000000000..6627c5a12 --- /dev/null +++ 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 index a558145b7..c0042b591 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>) -> Rc> { + let weight_constants = weight_functions + .into_iter() + .map(|w| w.weight_constants(Zero::zero(), 0)) + .collect(); + Rc::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..e84a76017 100644 --- a/feos-dft/src/interface/mod.rs +++ b/feos-dft/src/interface/mod.rs @@ -112,9 +112,9 @@ impl PlanarInterface { + 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 PlanarInterface { 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) } diff --git a/feos-dft/src/profile.rs b/feos-dft/src/profile.rs index 4b7f09970..cdcb12a21 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,43 @@ 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) + if rho_b.is_zero() { + Array::zeros(res.raw_dim()) + } else { + 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 +425,24 @@ 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 + println!("{bulk_density} {bulk_spec}"); + (bulk_density.mapv(f64::ln) - bulk_spec.mapv(f64::ln)).mapv(|r| { + if r.is_normal() { + r + } else { + 0.0 + } + }) } else { - chemical_potential.mapv(f64::exp) - mu_spec.mapv(f64::exp) + bulk_density - &bulk_spec }), ); @@ -478,40 +454,35 @@ where let solver = solver.cloned().unwrap_or_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 +501,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/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 } }) From 66598214ac936ba610b0a71c0bde98b12beb6532 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 6 Oct 2022 16:20:14 +0200 Subject: [PATCH 02/11] remove print and unnecessary ifs --- feos-dft/src/profile.rs | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/feos-dft/src/profile.rs b/feos-dft/src/profile.rs index cdcb12a21..6ffbe9b47 100644 --- a/feos-dft/src/profile.rs +++ b/feos-dft/src/profile.rs @@ -404,11 +404,7 @@ where .for_each(|(((mut res, df), &rho_b), rho)| { res.assign( &(if log { - if rho_b.is_zero() { - Array::zeros(res.raw_dim()) - } else { - rho.mapv(f64::ln) - rho_b.ln() - df.mapv(f64::ln) - } + rho.mapv(f64::ln) - rho_b.ln() - df.mapv(f64::ln) } else { &rho - rho_b * &df }), @@ -433,14 +429,7 @@ where res_bulk.assign( &(if log { - println!("{bulk_density} {bulk_spec}"); - (bulk_density.mapv(f64::ln) - bulk_spec.mapv(f64::ln)).mapv(|r| { - if r.is_normal() { - r - } else { - 0.0 - } - }) + bulk_density.mapv(f64::ln) - bulk_spec.mapv(f64::ln) } else { bulk_density - &bulk_spec }), From e5b453f01863375c97c5a3e79f6e42656f77a057 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 6 Oct 2022 17:01:50 +0200 Subject: [PATCH 03/11] make specification optional for planar interfaces --- feos-dft/src/interface/mod.rs | 23 +++++++++++++------ .../src/interface/surface_tension_diagram.rs | 5 +++- feos-dft/src/profile.rs | 8 ++++++- feos-dft/src/python/interface/mod.rs | 16 +++++++++---- .../interface/surface_tension_diagram.rs | 7 +++++- feos-dft/src/solver.rs | 14 +++++++++++ tests/gc_pcsaft/dft.rs | 5 ++-- tests/pcsaft/dft.rs | 14 ++++++----- 8 files changed, 70 insertions(+), 22 deletions(-) diff --git a/feos-dft/src/interface/mod.rs b/feos-dft/src/interface/mod.rs index e84a76017..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)?; @@ -112,14 +113,20 @@ impl PlanarInterface { + 0.5 * (rho_l + rho_v) }); - // // specify specification - // profile.profile.specification = - // DFTSpecifications::total_moles_from_profile(&profile.profile)?; + // specify specification + 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 { @@ -160,9 +167,11 @@ impl PlanarInterface { r, )?; - // // specify specification - // profile.profile.specification = - // DFTSpecifications::total_moles_from_profile(&profile.profile)?; + // specify specification + 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 6ffbe9b47..a44fb1b56 100644 --- a/feos-dft/src/profile.rs +++ b/feos-dft/src/profile.rs @@ -440,7 +440,13 @@ 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(); diff --git a/feos-dft/src/python/interface/mod.rs b/feos-dft/src/python/interface/mod.rs index f0787b4ab..9ae178bde 100644 --- a/feos-dft/src/python/interface/mod.rs +++ b/feos-dft/src/python/interface/mod.rs @@ -24,24 +24,29 @@ 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. /// /// 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 +59,18 @@ 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. /// /// 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..a6d6ee44b 100644 --- a/feos-dft/src/python/interface/surface_tension_diagram.rs +++ b/feos-dft/src/python/interface/surface_tension_diagram.rs @@ -19,6 +19,9 @@ 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. /// solver: DFTSolver, optional /// Custom solver options /// @@ -27,7 +30,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 +42,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 +52,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/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)?; From e19be620f491f5191694e4c68066262ba8db29c2 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 6 Oct 2022 22:48:43 +0200 Subject: [PATCH 04/11] include DFT guide in the documentation --- docs/index.md | 7 +++ docs/theory.md | 63 -------------------- docs/theory.rst | 154 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 161 insertions(+), 63 deletions(-) delete mode 100644 docs/theory.md create mode 100644 docs/theory.rst diff --git a/docs/index.md b/docs/index.md index f822f916b..89091e873 100644 --- a/docs/index.md +++ b/docs/index.md @@ -131,4 +131,11 @@ recipes/index rustguide/index rust_api +``` + +```{toctree} +:caption: Theory +:hidden: + +theory ``` \ No newline at end of file diff --git a/docs/theory.md b/docs/theory.md deleted file mode 100644 index 6627c5a12..000000000 --- a/docs/theory.md +++ /dev/null @@ -1,63 +0,0 @@ -# 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/docs/theory.rst b/docs/theory.rst new file mode 100644 index 000000000..a09d5c238 --- /dev/null +++ b/docs/theory.rst @@ -0,0 +1,154 @@ +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 :math:`\Omega` and the intrinsic Helmholtz energy :math:`F`. + +.. math:: + + \Omega(T,\mu,[\rho(r)])=F(T,[\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 :math:`T` and the density profiles :math:`\rho_i(r)` of the system and not on the external potential :math:`V_i^\mathrm{ext}(r)`. + +For a given temperature :math:`T`, chemical potentials :math:`\mu` and external potentials :math:`V^\mathrm{ext}(r)` the grand potential reaches a minimum at equilibrium. Mathematically this condition can be written as + +.. math:: + + \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 :math:`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, :math:`V^\mathrm{ext}(r)=0` and we get + +.. math:: + + F_{\rho_i}^\mathrm{b}-\mu_i=0 + +which can be inserted into (1) to give + +.. math:: + + 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: + +.. math:: + + \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 + +.. math:: + + \beta F_{\rho_i}=\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta F_{\rho_i}^\mathrm{res} + +.. math:: + + \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 + +.. math:: + + \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 + +.. math:: + + \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 + +.. math:: + + \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 + +.. math:: + + \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 + +.. math:: + + \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 + +.. math:: + + \beta F_{\rho_i}=m_i\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta\hat{F}_{\rho_i}^\mathrm{res} + +.. math:: + + \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 + +.. math:: + + \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 + +.. math:: + + \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 + +.. math:: + + 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 + +.. math:: + + \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: + +.. math:: + + \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 + +.. math:: + + \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) + +.. math:: + + 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: + +.. math:: + + \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$. From 31ae53c4a2379645ae0a98b8a6b25ab86c87c856 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 7 Oct 2022 14:43:06 +0200 Subject: [PATCH 05/11] include DFT theory in documentation --- docs/conf.py | 1 + docs/dft/convolution_integrals.md | 1 + docs/dft/euler_lagrange_equation.md | 98 ++++++++++++++++++ docs/dft/functional_derivatives.md | 1 + docs/dft/index.md | 12 +++ docs/index.md | 2 +- docs/theory.rst | 154 ---------------------------- 7 files changed, 114 insertions(+), 155 deletions(-) create mode 100644 docs/dft/convolution_integrals.md create mode 100644 docs/dft/euler_lagrange_equation.md create mode 100644 docs/dft/functional_derivatives.md create mode 100644 docs/dft/index.md delete mode 100644 docs/theory.rst 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/dft/convolution_integrals.md b/docs/dft/convolution_integrals.md new file mode 100644 index 000000000..6936d411c --- /dev/null +++ b/docs/dft/convolution_integrals.md @@ -0,0 +1 @@ +## Convolution integrals \ No newline at end of file diff --git a/docs/dft/euler_lagrange_equation.md b/docs/dft/euler_lagrange_equation.md new file mode 100644 index 000000000..e63af19d5 --- /dev/null +++ b/docs/dft/euler_lagrange_equation.md @@ -0,0 +1,98 @@ +## 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)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 $\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 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/docs/dft/functional_derivatives.md b/docs/dft/functional_derivatives.md new file mode 100644 index 000000000..93fd5753b --- /dev/null +++ b/docs/dft/functional_derivatives.md @@ -0,0 +1 @@ +## Functional derivatives \ No newline at end of file diff --git a/docs/dft/index.md b/docs/dft/index.md new file mode 100644 index 000000000..0bbc42c62 --- /dev/null +++ b/docs/dft/index.md @@ -0,0 +1,12 @@ +# Classical DFT +In this section, the implementation of the core expressions from classical density functional theory in $\text{FeO}_\text{s}$ is explained. + +```{eval-rst} +.. toctree:: + :maxdepth: 1 + :hidden: + + euler_lagrange_equation + functional_derivatives + convolution_integrals +``` \ No newline at end of file diff --git a/docs/index.md b/docs/index.md index 89091e873..0490849dc 100644 --- a/docs/index.md +++ b/docs/index.md @@ -137,5 +137,5 @@ rust_api :caption: Theory :hidden: -theory +dft/index ``` \ No newline at end of file diff --git a/docs/theory.rst b/docs/theory.rst deleted file mode 100644 index a09d5c238..000000000 --- a/docs/theory.rst +++ /dev/null @@ -1,154 +0,0 @@ -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 :math:`\Omega` and the intrinsic Helmholtz energy :math:`F`. - -.. math:: - - \Omega(T,\mu,[\rho(r)])=F(T,[\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 :math:`T` and the density profiles :math:`\rho_i(r)` of the system and not on the external potential :math:`V_i^\mathrm{ext}(r)`. - -For a given temperature :math:`T`, chemical potentials :math:`\mu` and external potentials :math:`V^\mathrm{ext}(r)` the grand potential reaches a minimum at equilibrium. Mathematically this condition can be written as - -.. math:: - - \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 :math:`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, :math:`V^\mathrm{ext}(r)=0` and we get - -.. math:: - - F_{\rho_i}^\mathrm{b}-\mu_i=0 - -which can be inserted into (1) to give - -.. math:: - - 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: - -.. math:: - - \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 - -.. math:: - - \beta F_{\rho_i}=\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta F_{\rho_i}^\mathrm{res} - -.. math:: - - \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 - -.. math:: - - \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 - -.. math:: - - \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 - -.. math:: - - \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 - -.. math:: - - \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 - -.. math:: - - \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 - -.. math:: - - \beta F_{\rho_i}=m_i\ln\left(\rho_i(r)\Lambda_i^3\right)+\beta\hat{F}_{\rho_i}^\mathrm{res} - -.. math:: - - \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 - -.. math:: - - \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 - -.. math:: - - \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 - -.. math:: - - 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 - -.. math:: - - \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: - -.. math:: - - \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 - -.. math:: - - \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) - -.. math:: - - 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: - -.. math:: - - \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$. From c232b56b23878172e5de5900b63f48eac731395d Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sat, 8 Oct 2022 12:35:21 +0200 Subject: [PATCH 06/11] add section on functional derivatives --- docs/dft/functional_derivatives.md | 40 +++++++++++++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/docs/dft/functional_derivatives.md b/docs/dft/functional_derivatives.md index 93fd5753b..93e17bc93 100644 --- a/docs/dft/functional_derivatives.md +++ b/docs/dft/functional_derivatives.md @@ -1 +1,39 @@ -## Functional derivatives \ No newline at end of file +## Functional derivatives + +In the last section the functional derivative + +$$\hat F_{\rho_\alpha}^\mathrm{res}(r)=\left(\frac{\delta 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 + +$$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$$ + +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. 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 From 68a373b02d784302df0f4d01a7f2a4b19bd027eb Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Tue, 18 Oct 2022 13:07:55 +0200 Subject: [PATCH 07/11] use Arc in BulkConvolver --- feos-dft/src/convolver/mod.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/feos-dft/src/convolver/mod.rs b/feos-dft/src/convolver/mod.rs index c0042b591..5cdaa56d2 100644 --- a/feos-dft/src/convolver/mod.rs +++ b/feos-dft/src/convolver/mod.rs @@ -41,12 +41,12 @@ pub(crate) struct BulkConvolver { } impl> BulkConvolver { - pub(crate) fn new(weight_functions: Vec>) -> Rc> { + pub(crate) fn new(weight_functions: Vec>) -> Arc> { let weight_constants = weight_functions .into_iter() .map(|w| w.weight_constants(Zero::zero(), 0)) .collect(); - Rc::new(Self { weight_constants }) + Arc::new(Self { weight_constants }) } } From 1be7bdd6840624be542d64298804d92f8f6e5180 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 10 Nov 2022 20:56:28 +0100 Subject: [PATCH 08/11] udpate theory guide toc --- docs/dft/index.md | 12 ------------ docs/index.md | 6 ++++-- docs/{ => theory}/dft/convolution_integrals.md | 0 docs/{ => theory}/dft/euler_lagrange_equation.md | 0 docs/{ => theory}/dft/functional_derivatives.md | 0 docs/theory/dft/index.md | 12 ++++++++++++ docs/theory/eos/index.md | 9 +++++++++ docs/theory/models/index.md | 9 +++++++++ 8 files changed, 34 insertions(+), 14 deletions(-) delete mode 100644 docs/dft/index.md rename docs/{ => theory}/dft/convolution_integrals.md (100%) rename docs/{ => theory}/dft/euler_lagrange_equation.md (100%) rename docs/{ => theory}/dft/functional_derivatives.md (100%) create mode 100644 docs/theory/dft/index.md create mode 100644 docs/theory/eos/index.md create mode 100644 docs/theory/models/index.md diff --git a/docs/dft/index.md b/docs/dft/index.md deleted file mode 100644 index 0bbc42c62..000000000 --- a/docs/dft/index.md +++ /dev/null @@ -1,12 +0,0 @@ -# Classical DFT -In this section, the implementation of the core expressions from classical density functional theory in $\text{FeO}_\text{s}$ is explained. - -```{eval-rst} -.. toctree:: - :maxdepth: 1 - :hidden: - - euler_lagrange_equation - functional_derivatives - convolution_integrals -``` \ No newline at end of file diff --git a/docs/index.md b/docs/index.md index 0490849dc..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} @@ -137,5 +137,7 @@ rust_api :caption: Theory :hidden: -dft/index +theory/eos/index +theory/dft/index +theory/models/index ``` \ No newline at end of file diff --git a/docs/dft/convolution_integrals.md b/docs/theory/dft/convolution_integrals.md similarity index 100% rename from docs/dft/convolution_integrals.md rename to docs/theory/dft/convolution_integrals.md diff --git a/docs/dft/euler_lagrange_equation.md b/docs/theory/dft/euler_lagrange_equation.md similarity index 100% rename from docs/dft/euler_lagrange_equation.md rename to docs/theory/dft/euler_lagrange_equation.md diff --git a/docs/dft/functional_derivatives.md b/docs/theory/dft/functional_derivatives.md similarity index 100% rename from docs/dft/functional_derivatives.md rename to docs/theory/dft/functional_derivatives.md 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 From 4433cbe2a28103698884f2e5e216b8fc2c00ba46 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 10 Nov 2022 20:57:16 +0100 Subject: [PATCH 09/11] remove stub --- docs/theory/dft/convolution_integrals.md | 1 - 1 file changed, 1 deletion(-) delete mode 100644 docs/theory/dft/convolution_integrals.md diff --git a/docs/theory/dft/convolution_integrals.md b/docs/theory/dft/convolution_integrals.md deleted file mode 100644 index 6936d411c..000000000 --- a/docs/theory/dft/convolution_integrals.md +++ /dev/null @@ -1 +0,0 @@ -## Convolution integrals \ No newline at end of file From 519405bac712ed953774d78a2981044857f3bb76 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Thu, 10 Nov 2022 21:00:52 +0100 Subject: [PATCH 10/11] Add changelog entry --- feos-dft/CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) 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 From a080fa487c03e0b1537ced556c855659f15fffa9 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Fri, 11 Nov 2022 16:09:29 +0100 Subject: [PATCH 11/11] incorporate comments from review --- docs/theory/dft/euler_lagrange_equation.md | 29 ++++++++++--------- docs/theory/dft/functional_derivatives.md | 15 ++++++---- feos-dft/src/python/interface/mod.rs | 2 ++ .../interface/surface_tension_diagram.rs | 1 + 4 files changed, 28 insertions(+), 19 deletions(-) diff --git a/docs/theory/dft/euler_lagrange_equation.md b/docs/theory/dft/euler_lagrange_equation.md index e63af19d5..9d9c5f051 100644 --- a/docs/theory/dft/euler_lagrange_equation.md +++ b/docs/theory/dft/euler_lagrange_equation.md @@ -2,7 +2,7 @@ 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)dr +\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}$. @@ -22,11 +22,11 @@ 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: +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)dr+\beta F^\mathrm{res}$$ +$$\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}$$ -The functional derivatives for an inhomogeneous and a bulk system follow as +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}$$ @@ -41,20 +41,23 @@ 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$$ +$$\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 chain length, $y_{ii}$ the cavity correlation function at contact in the reference fluid, and $\lambda_i$ a weighted density. +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 -$$\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}}$$ +$$\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)dr+\beta F^\mathrm{chain}+\beta F^\mathrm{res}$$ +$$\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)dr+\underbrace{\beta\hat{F}^\mathrm{chain}+\beta F^\mathrm{res}}_{\beta\hat{F}^\mathrm{res}}$$ +$$\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 @@ -67,20 +70,20 @@ 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 +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')dr$$ +$$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 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: +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)$$ @@ -88,7 +91,7 @@ At this point it can be numerically useful to redistribute the bulk contribution $$\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$$ +$$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: diff --git a/docs/theory/dft/functional_derivatives.md b/docs/theory/dft/functional_derivatives.md index 93e17bc93..ab1a145fc 100644 --- a/docs/theory/dft/functional_derivatives.md +++ b/docs/theory/dft/functional_derivatives.md @@ -2,13 +2,13 @@ In the last section the functional derivative -$$\hat F_{\rho_\alpha}^\mathrm{res}(r)=\left(\frac{\delta F^\mathrm{res}}{\delta\rho_\alpha(r)}\right)_{T,\rho_{\alpha'\neq\alpha}}$$ +$$\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$ +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}$$ @@ -18,11 +18,14 @@ For Helmholtz energy functionals that can be written in this form, the calculati $$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 +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 -$$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$$ +$$\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 +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}$$ @@ -34,6 +37,6 @@ With this distinction, the calculation of the functional derivative is split int 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. 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. +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/feos-dft/src/python/interface/mod.rs b/feos-dft/src/python/interface/mod.rs index 9ae178bde..1a3a305ec 100644 --- a/feos-dft/src/python/interface/mod.rs +++ b/feos-dft/src/python/interface/mod.rs @@ -27,6 +27,7 @@ macro_rules! impl_planar_interface { /// fix_equimolar_surface: bool, optional /// If True use additional constraints to fix the /// equimolar surface of the system. + /// Defaults to False. /// /// Returns /// ------- @@ -62,6 +63,7 @@ macro_rules! impl_planar_interface { /// fix_equimolar_surface: bool, optional /// If True use additional constraints to fix the /// equimolar surface of the system. + /// Defaults to False. /// /// Returns /// ------- diff --git a/feos-dft/src/python/interface/surface_tension_diagram.rs b/feos-dft/src/python/interface/surface_tension_diagram.rs index a6d6ee44b..7350b2f67 100644 --- a/feos-dft/src/python/interface/surface_tension_diagram.rs +++ b/feos-dft/src/python/interface/surface_tension_diagram.rs @@ -22,6 +22,7 @@ macro_rules! impl_surface_tension_diagram { /// 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 ///