prehner · Oct 26, 2026 · Oct 15, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -10,6 +10,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added the `Grid` and `PoreProfile` classes as an entry points to DFT in various coordinate systems in Python. [#376](https://github.com/feos-org/feos/pull/376)
 - `PoreProfile.solver_log` now stores a timing profile of the different solver steps. [#385](https://github.com/feos-org/feos/pull/385)
 - Added `Interface::curved` to calculate cylindrical and spherical interfaces from a planar interface. [#387](https://github.com/feos-org/feos/pull/387)
+- Added `drho_dn` and `drho_dt_n` to DFT profiles to calculate derivatives of density profiles at constant number of particles. [#403](https://github.com/feos-org/feos/pull/403)
+- Added the association functionals of Yu and Wu, WDA, and iSAFT with the fractions of non-bonded sites as additional variables of the density profile. The version is selected with `AssociationVersion` (default iSAFT). [#403](https://github.com/feos-org/feos/pull/403)
+- Added `non_bonded_fractions` to DFT profiles to access the fractions of non-bonded association sites. [#403](https://github.com/feos-org/feos/pull/403)
 
 ### Changed
 - Changed data type of initial temperatures or pressure for phase equilibrium calculations (`TemperatureOrPressure::Other`) from `D` to `f64`. [#369](https://github.com/feos-org/feos/pull/369)
@@ -20,11 +23,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Renamed `PlanarInterface` to `Interface` and its constructors `from_tanh` and `from_pdgt` to `planar_from_tanh` and `planar_from_pdgt`. [#387](https://github.com/feos-org/feos/pull/387)
 - Replaced `PlanarInterface.from_density_profile` in Python with the `Interface(grid, vle, density)` constructor in analogy to `PoreProfile`. [#387](https://github.com/feos-org/feos/pull/387)
 - `drho_dp`, `dn_dp`, `drho_dt`, and `dn_dt` of DFT profiles now take the bulk state as an argument (methods instead of properties in Python). [#387](https://github.com/feos-org/feos/pull/387)
+- Renamed `drho_dt` and `dn_dt` of DFT profiles to `drho_dt_p` and `dn_dt_p` (derivatives at constant pressure and bulk composition). [#403](https://github.com/feos-org/feos/pull/403)
+- The Newton solver is now the default solver for DFT (previously Anderson mixing). It always iterates the logarithms of the density profiles (and of the fractions of non-bonded sites), so the `log` argument of `DFTSolver.newton` is removed. Jacobian-vector products are calculated exactly from the linearized residual, and the step size is controlled by a line search. [#403](https://github.com/feos-org/feos/pull/403)
 
 ### Removed
 - Removed the `DFTSpecification` trait in favor of only using the `DFTSpecification` enum (renamed from `DFTSpecifications`). [#371](https://github.com/feos-org/feos/pull/371)
 - Removed the `Pore2D` and `Pore3D` interfaces including the (limited) calculation of external potentials for complex pore grometries. [#376](https://github.com/feos-org/feos/pull/376)
 - Removed the free energy-averaged external (FEA) potential. [#376](https://github.com/feos-org/feos/pull/376)
+- Removed predictive density gradient theory (pDGT) including `Interface::planar_from_pdgt`. Interfaces are initialized with `Interface::planar_from_tanh`. [#403](https://github.com/feos-org/feos/pull/403)
+- Removed the association functional of Yu and Wu as functional contribution (`YuWuAssociationFunctional`). Association in DFT is now handled by `SaftAssociationFunctional` with the fractions of non-bonded sites as additional variables. [#403](https://github.com/feos-org/feos/pull/403)
 
 ### Fixed
 - The Newton solver for cDFT now correctly handles specifications other than the chemical potential. [#384](https://github.com/feos-org/feos/pull/384)

diff --git a/crates/feos-derive/src/association.rs b/crates/feos-derive/src/association.rs
@@ -0,0 +1,75 @@
+use quote::quote;
+use syn::{Ident, Type, Visibility};
+
+/// Enum `AssociationFunctionalVariant` of the association functionals of the given variants
+/// (name and type of the Helmholtz energy functional) that dispatches the methods of
+/// `AssociationFunctional`, or `NoAssociation` if there are no such variants.
+pub(crate) fn impl_association_variant(
+    vis: &Visibility,
+    variants: &[(&Ident, &Type)],
+) -> proc_macro2::TokenStream {
+    if variants.is_empty() {
+        return quote! {
+            #vis type AssociationFunctionalVariant<'a> = feos_dft::NoAssociation;
+        };
+    }
+
+    let enum_variants = variants.iter().map(|(name, ty)| {
+        quote! {
+            #name(<#ty as HelmholtzEnergyFunctionalDyn>::Association<'a>)
+        }
+    });
+    let arms = |method: proc_macro2::TokenStream| {
+        variants
+            .iter()
+            .map(|(name, _)| {
+                quote! {
+                    Self::#name(ref association) => association.#method
+                }
+            })
+            .collect::<Vec<_>>()
+    };
+    let sites = arms(quote!(sites()));
+    let weight_functions = arms(quote!(weight_functions(temperature)));
+    let contact_distances = arms(quote!(contact_distances(temperature)));
+    let contact_values = arms(quote!(contact_values(temperature, weighted_densities)));
+    let association_strength = arms(quote!(association_strength(temperature)));
+
+    quote! {
+        #vis enum AssociationFunctionalVariant<'a> {
+            #(#enum_variants,)*
+        }
+
+        impl feos_dft::AssociationFunctional for AssociationFunctionalVariant<'_> {
+            fn sites(&self) -> feos_dft::AssociationSites<'_> {
+                match *self {
+                    #(#sites,)*
+                }
+            }
+            fn weight_functions<N: DualNum<Primitive = f64> + Copy>(&self, temperature: N) -> feos_dft::WeightFunctionInfo<N> {
+                match *self {
+                    #(#weight_functions,)*
+                }
+            }
+            fn contact_distances<N: DualNum<Primitive = f64> + Copy>(&self, temperature: N) -> Option<ndarray::Array1<N>> {
+                match *self {
+                    #(#contact_distances,)*
+                }
+            }
+            fn contact_values<N: DualNum<Primitive = f64> + Copy>(
+                &self,
+                temperature: N,
+                weighted_densities: ndarray::ArrayView2<N>,
+            ) -> (ndarray::Array2<N>, ndarray::Array2<N>) {
+                match *self {
+                    #(#contact_values,)*
+                }
+            }
+            fn association_strength<N: DualNum<Primitive = f64> + Copy>(&self, temperature: N) -> ndarray::Array1<N> {
+                match *self {
+                    #(#association_strength,)*
+                }
+            }
+        }
+    }
+}
diff --git a/crates/feos-derive/src/dft.rs b/crates/feos-derive/src/dft.rs
@@ -1,3 +1,4 @@
+use crate::association::impl_association_variant;
 use crate::{OPT_IMPLS, implement};
 use quote::quote;
 use syn::DeriveInput;
@@ -10,15 +11,34 @@ pub(crate) fn expand_helmholtz_energy_functional(
     };
 
     let functional = impl_helmholtz_energy_functional(&input.ident, variants)?;
+    let association = impl_association(&input.vis, variants)?;
     let fluid_parameters = impl_fluid_parameters(&input.ident, variants)?;
     let pair_potential = impl_pair_potential(&input.ident, variants)?;
     Ok(quote! {
         #functional
+        #association
         #fluid_parameters
         #pair_potential
     })
 }
 
+/// The enum of the association functionals of all variants that implement `association`.
+fn impl_association(
+    vis: &syn::Visibility,
+    variants: &syn::punctuated::Punctuated<syn::Variant, syn::token::Comma>,
+) -> syn::Result<proc_macro2::TokenStream> {
+    let mut association_variants = Vec::new();
+    for v in variants.iter() {
+        if implement("association", v, &OPT_IMPLS)? {
+            let syn::Fields::Unnamed(syn::FieldsUnnamed { unnamed, .. }) = &v.fields else {
+                panic!("All variants must be tuple structs!")
+            };
+            association_variants.push((&v.ident, &unnamed.first().unwrap().ty));
+        }
+    }
+    Ok(impl_association_variant(vis, &association_variants))
+}
+
 pub(crate) fn impl_helmholtz_energy_functional(
     ident: &syn::Ident,
     variants: &syn::punctuated::Punctuated<syn::Variant, syn::token::Comma>,
@@ -54,9 +74,20 @@ pub(crate) fn impl_helmholtz_energy_functional(
         }
     }
 
+    let mut association = Vec::new();
+    for v in variants.iter() {
+        if implement("association", v, &OPT_IMPLS)? {
+            let name = &v.ident;
+            association.push(quote! {
+                Self::#name(functional) => functional.association().map(AssociationFunctionalVariant::#name)
+            });
+        }
+    }
+
     Ok(quote! {
         impl HelmholtzEnergyFunctionalDyn for #ident {
             type Contribution<'a> = FunctionalContributionVariant<'a>;
+            type Association<'a> = AssociationFunctionalVariant<'a>;
             fn molecule_shape(&self) -> feos_dft::MoleculeShape<'_> {
                 match self {
                     #(#molecule_shape,)*
@@ -73,6 +104,12 @@ pub(crate) fn impl_helmholtz_energy_functional(
                     _ => petgraph::Graph::with_capacity(0, 0),
                 }
             }
+            fn association(&self) -> Option<Self::Association<'_>> {
+                match self {
+                    #(#association,)*
+                    _ => None,
+                }
+            }
         }
     })
 }

diff --git a/crates/feos-derive/src/functional_contribution.rs b/crates/feos-derive/src/functional_contribution.rs
@@ -36,12 +36,6 @@ fn impl_functional_contribution(
             Self::#name(functional_contribution) => functional_contribution.weight_functions(temperature)
         }
     });
-    let weight_functions_pdgt = variants.iter().map(|v| {
-        let name = &v.ident;
-        quote! {
-            Self::#name(functional_contribution) => functional_contribution.weight_functions_pdgt(temperature)
-        }
-    });
     let helmholtz_energy_density = variants.iter().map(|v| {
         let name = &v.ident;
         quote! {
@@ -61,11 +55,6 @@ fn impl_functional_contribution(
                     #(#weight_functions,)*
                 }
             }
-            fn weight_functions_pdgt<N: DualNum<Primitive = f64> + Copy>(&self, temperature: N) -> feos_dft::WeightFunctionInfo<N> {
-                match self {
-                    #(#weight_functions_pdgt,)*
-                }
-            }
             fn helmholtz_energy_density<N: DualNum<Primitive = f64> + Copy>(
                 &self,
                 temperature: N,

diff --git a/crates/feos-derive/src/lib.rs b/crates/feos-derive/src/lib.rs
@@ -8,21 +8,23 @@ use ideal_gas::expand_ideal_gas;
 use proc_macro::TokenStream;
 use residual::expand_residual;
 use subset::expand_subset;
-use syn::{parse_macro_input, DeriveInput};
+use syn::{DeriveInput, parse_macro_input};
 
+mod association;
 mod dft;
 mod functional_contribution;
 mod ideal_gas;
 mod residual;
 mod subset;
 
 // possible additional traits to implement
-const OPT_IMPLS: [&str; 7] = [
+const OPT_IMPLS: [&str; 8] = [
     "molar_weight",
     "parameter_info",
     "entropy_scaling",
     "functional",
     "bond_lengths",
+    "association",
     "fluid_parameters",
     "pair_potential",
 ];