From 5af3237ee6ba135299cc5ca3ee9e51003476a422 Mon Sep 17 00:00:00 2001 From: Philipp Rehner <69816385+prehner@users.noreply.github.com> Date: Fri, 14 Aug 2026 16:50:43 +0200 Subject: [PATCH 1/2] Fix calculation of bulk state when moles are specified (#382) --- crates/feos-dft/src/adsorption/pore.rs | 1 + crates/feos-dft/src/convolver/mod.rs | 61 +----- .../src/convolver/periodic_convolver.rs | 16 +- crates/feos-dft/src/geometry.rs | 14 +- crates/feos-dft/src/profile/mod.rs | 174 ++++++++++-------- crates/feos-dft/src/profile/properties.rs | 17 +- crates/feos-dft/src/solver.rs | 38 ++-- py-feos/src/dft/adsorption/pore.rs | 7 + py-feos/src/dft/profile.rs | 6 +- 9 files changed, 156 insertions(+), 178 deletions(-) diff --git a/crates/feos-dft/src/adsorption/pore.rs b/crates/feos-dft/src/adsorption/pore.rs index fc3cf93be..b5d37630c 100644 --- a/crates/feos-dft/src/adsorption/pore.rs +++ b/crates/feos-dft/src/adsorption/pore.rs @@ -162,6 +162,7 @@ where } pub fn update_bulk(mut self, bulk: &State) -> Self { + self.profile.specification = DFTSpecification::from_state(bulk); self.profile.bulk = bulk.clone(); self.grand_potential = None; self.interfacial_tension = None; diff --git a/crates/feos-dft/src/convolver/mod.rs b/crates/feos-dft/src/convolver/mod.rs index c047bcd0e..540ce4d97 100644 --- a/crates/feos-dft/src/convolver/mod.rs +++ b/crates/feos-dft/src/convolver/mod.rs @@ -1,10 +1,8 @@ use crate::geometry::{Axis, Geometry, Grid}; use crate::weight_functions::*; -use ndarray::linalg::Dot; use ndarray::prelude::*; use ndarray::{Axis as Axis_nd, RemoveAxis, Slice}; use num_dual::*; -use num_traits::Zero; use rustdct::DctNum; use std::ops::{AddAssign, MulAssign, SubAssign}; use std::sync::Arc; @@ -36,46 +34,6 @@ pub trait Convolver: Send + Sync { ) -> Array; } -pub(crate) struct BulkConvolver { - weight_constants: Vec>, -} - -impl + Copy + Send + Sync> BulkConvolver { - #[expect(clippy::new_ret_no_self)] - pub(crate) fn new(weight_functions: Vec>) -> Arc> { - let weight_constants = weight_functions - .into_iter() - .map(|w| w.weight_constants(Zero::zero(), 0)) - .collect(); - Arc::new(Self { weight_constants }) - } -} - -impl + Copy + Send + Sync> 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`. @@ -125,11 +83,10 @@ pub struct ConvolverFFT { cartesian_transforms: Vec>, } -impl ConvolverFFT +impl ConvolverFFT where T: DctNum + DualNum, D::Larger: Dimension, - D::Smaller: Dimension, ::Larger: Dimension, { /// Create the appropriate FFT convolver for the given grid. @@ -139,9 +96,11 @@ where lanczos: Option, ) -> Arc> { match grid { + Grid::Bulk => PeriodicConvolver::new_0d(weight_functions), Grid::Polar(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos), Grid::Spherical(r) => CurvilinearConvolver::new(r, &[], weight_functions, lanczos), Grid::Cartesian1(z) => Self::new(Some(z), &[], weight_functions, lanczos), + Grid::Periodical1(z) => PeriodicConvolver::new_1d(z, weight_functions, lanczos), Grid::Cylindrical { r, z } => { CurvilinearConvolver::new(r, &[z], weight_functions, lanczos) } @@ -553,11 +512,10 @@ struct CurvilinearConvolver { convolver_boundary: Arc>, } -impl CurvilinearConvolver +impl CurvilinearConvolver where T: DctNum + DualNum, D::Larger: Dimension, - D::Smaller: Dimension, ::Larger: Dimension, { #[expect(clippy::new_ret_no_self)] @@ -574,10 +532,9 @@ where } } -impl Convolver for CurvilinearConvolver +impl Convolver for CurvilinearConvolver where T: DctNum + DualNum, - D::Smaller: Dimension, D::Larger: Dimension, { fn convolve( @@ -587,9 +544,9 @@ where ) -> Array { // subtract boundary profile from full profile let profile_boundary = profile - .index_axis(Axis(0), profile.shape()[0] - 1) + .slice_axis(Axis(0), Slice::from(profile.shape()[0] - 1..)) .into_owned(); - for mut lane in profile.outer_iter_mut() { + for mut lane in profile.axis_chunks_iter_mut(Axis(0), 1) { lane.sub_assign(&profile_boundary); } @@ -597,14 +554,12 @@ where let mut result = self.convolver.convolve(profile, weight_function); // convolve boundary profile - let profile_boundary = profile_boundary.insert_axis(Axis(0)); let result_boundary = self .convolver_boundary .convolve(profile_boundary, weight_function); // Add boundary result back to full result - let result_boundary = result_boundary.index_axis(Axis(0), 0); - for mut lane in result.outer_iter_mut() { + for mut lane in result.axis_chunks_iter_mut(Axis(0), 1) { lane.add_assign(&result_boundary); } result diff --git a/crates/feos-dft/src/convolver/periodic_convolver.rs b/crates/feos-dft/src/convolver/periodic_convolver.rs index dfb12ff06..7678797a7 100644 --- a/crates/feos-dft/src/convolver/periodic_convolver.rs +++ b/crates/feos-dft/src/convolver/periodic_convolver.rs @@ -30,6 +30,18 @@ where D::Larger: Dimension, ::Larger: Dimension, { + pub fn new_0d(weight_functions: &[WeightFunctionInfo]) -> Arc> { + Self::new(&[], |_| {}, weight_functions, None) + } + + pub fn new_1d( + axis: &Axis, + weight_functions: &[WeightFunctionInfo], + lanczos: Option, + ) -> Arc> { + Self::new(&[axis], |_| {}, weight_functions, lanczos) + } + pub fn new_2d( axes: &[&Axis], angle: Angle, @@ -200,7 +212,7 @@ impl PeriodicConvolver { } fn forward_transform(&self, f: ArrayView) -> Array, D2> { - let offset = D2::NDIM.unwrap() - D::NDIM.unwrap(); + let offset = f.ndim() - self.k_abs.ndim(); let mut result = f.mapv(Complex::from); for (i, transform) in self.forward_transforms.iter().enumerate() { for r in result.lanes_mut(Axis_nd(i + offset)).into_iter() { @@ -211,7 +223,7 @@ impl PeriodicConvolver { } fn inverse_transform(&self, mut f: Array, D2>) -> Array { - let offset = D2::NDIM.unwrap() - D::NDIM.unwrap(); + let offset = f.ndim() - self.k_abs.ndim(); for (i, transform) in self.inverse_transforms.iter().enumerate() { for r in f.lanes_mut(Axis_nd(i + offset)).into_iter() { self.transform(transform, r); diff --git a/crates/feos-dft/src/geometry.rs b/crates/feos-dft/src/geometry.rs index e69409e9b..b4a14f5d1 100644 --- a/crates/feos-dft/src/geometry.rs +++ b/crates/feos-dft/src/geometry.rs @@ -6,7 +6,9 @@ use std::f64::consts::{FRAC_PI_3, PI}; /// Grids with up to three dimensions. #[derive(Clone)] pub enum Grid { + Bulk, Cartesian1(Axis), + Periodical1(Axis), Cartesian2(Axis, Axis), Periodical2(Axis, Axis, Angle), Cartesian3(Axis, Axis, Axis), @@ -27,7 +29,8 @@ impl Grid { pub fn axes(&self) -> Vec<&Axis> { match self { - Self::Cartesian1(x) => vec![x], + Self::Bulk => vec![], + Self::Cartesian1(x) | Self::Periodical1(x) => vec![x], Self::Cartesian2(x, y) | Self::Periodical2(x, y, _) => vec![x, y], Self::Cartesian3(x, y, z) | Self::Periodical3(x, y, z, _) => vec![x, y, z], Self::Spherical(r) | Self::Polar(r) => vec![r], @@ -37,7 +40,8 @@ impl Grid { pub fn axes_mut(&mut self) -> Vec<&mut Axis> { match self { - Self::Cartesian1(x) => vec![x], + Self::Bulk => vec![], + Self::Cartesian1(x) | Self::Periodical1(x) => vec![x], Self::Cartesian2(x, y) | Self::Periodical2(x, y, _) => vec![x, y], Self::Cartesian3(x, y, z) | Self::Periodical3(x, y, z, _) => vec![x, y, z], Self::Spherical(r) | Self::Polar(r) => vec![r], @@ -72,7 +76,11 @@ impl Grid { pub fn mesh(&self) -> Vec>> { match self { - Grid::Cartesian1(ax) | Grid::Spherical(ax) | Grid::Polar(ax) => { + Self::Bulk => vec![], + Grid::Cartesian1(ax) + | Self::Periodical1(ax) + | Grid::Spherical(ax) + | Grid::Polar(ax) => { vec![Length::from_reduced(ax.grid.clone()).into_dyn()] } Grid::Cartesian2(u, v) => mesh_2d(u, v, 90.0 * DEGREES), diff --git a/crates/feos-dft/src/profile/mod.rs b/crates/feos-dft/src/profile/mod.rs index 37f5e3aa9..6def8ded2 100644 --- a/crates/feos-dft/src/profile/mod.rs +++ b/crates/feos-dft/src/profile/mod.rs @@ -1,10 +1,10 @@ -use crate::convolver::{BulkConvolver, Convolver, ConvolverFFT}; +use crate::convolver::{Convolver, ConvolverFFT, PeriodicConvolver}; use crate::functional::HelmholtzEnergyFunctional; use crate::geometry::Grid; use crate::solver::{DFTSolver, DFTSolverLog}; use feos_core::{FeosError, FeosResult, ReferenceSystem, State}; use nalgebra::{DVector, Dyn, U1}; -use ndarray::{Array, Array1, ArrayBase, Axis as Axis_nd, Data, Dimension, RemoveAxis}; +use ndarray::{Array, Array1, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix0, arr1}; use num_dual::DualNum; use quantity::{_Volume, Density, Energy, Entropy, Length, Moles, Quantity, Temperature, Volume}; use std::ops::{Add, MulAssign}; @@ -22,7 +22,7 @@ const MAX_POTENTIAL: f64 = 50.0; #[derive(Clone)] pub enum DFTSpecification { /// DFT with specified chemical potential. - ChemicalPotential, + ChemicalPotential(Array1), /// DFT with specified number of particles. /// /// The solution is still a grand canonical density profile, but the chemical @@ -30,21 +30,42 @@ pub enum DFTSpecification { /// with the specified number of particles. Moles(Array1), /// DFT with specified total number of moles. - TotalMoles(f64), + TotalMoles(f64, Array1), } impl DFTSpecification { - fn calculate_bulk_density( - &self, - bulk_density: &Array1, - z: &Array1, - ) -> FeosResult> { + fn calculate_fugacity(&self, z: &Array1) -> FeosResult> { Ok(match self { - Self::ChemicalPotential => bulk_density.clone(), + Self::ChemicalPotential(fugacity) => fugacity.clone(), Self::Moles(moles) => moles / z, - Self::TotalMoles(total_moles) => bulk_density * *total_moles / (bulk_density * z).sum(), + Self::TotalMoles(total_moles, fugacity) => { + fugacity * *total_moles / (fugacity * z).sum() + } }) } + + pub fn from_state(state: &State) -> Self { + let component_index = state.eos.component_index().into_owned(); + let m = arr1(&state.eos.m()); + let partial_density = state.partial_density().into_reduced(); + let temperature = state.temperature.into_reduced(); + let bulk_density = component_index + .iter() + .map(|&i| partial_density[i]) + .collect(); + let bulk_convolver = + PeriodicConvolver::<_, Ix0>::new_0d(&state.eos.weight_functions(temperature)); + let (_, dfdrho_bulk) = state + .eos + .functional_derivative(temperature, &bulk_density, bulk_convolver.as_ref()) + .unwrap(); + let exp_dfdrho = (dfdrho_bulk / m).mapv(f64::exp); + let bonds = state + .eos + .bond_integrals(temperature, &exp_dfdrho, bulk_convolver.as_ref()); + let fugacity = bulk_density * exp_dfdrho * bonds; + Self::ChemicalPotential(fugacity) + } } #[derive(Clone)] @@ -80,10 +101,9 @@ impl DFTProfile { } } -impl DFTProfile +impl DFTProfile where D::Larger: Dimension, - D::Smaller: Dimension, ::Larger: Dimension, { /// Create a new density profile. @@ -147,7 +167,7 @@ where convolver, temperature: bulk.temperature, density, - specification: DFTSpecification::ChemicalPotential, + specification: DFTSpecification::from_state(bulk), external_potential, bulk: bulk.clone(), solver_log: None, @@ -167,7 +187,12 @@ where pub fn fix_total_moles(&mut self) { let rho = self.density.to_reduced(); let moles = self.integrate_reduced_comp(&rho).sum(); - self.specification = DFTSpecification::TotalMoles(moles); + let DFTSpecification::ChemicalPotential(fugacity) = + DFTSpecification::from_state(&self.bulk) + else { + unreachable!() + }; + self.specification = DFTSpecification::TotalMoles(moles, fugacity); } /// Return the external potential in SI units. @@ -293,37 +318,17 @@ where .weighted_densities(&self.density.to_reduced())) } - #[expect(clippy::type_complexity)] - pub fn residual(&self, log: bool) -> FeosResult<(Array, Array1, f64)> { - // Read from profile + pub fn residual(&self, log: bool) -> FeosResult<(Array, f64)> { let density = self.density.to_reduced(); - let partial_density = self.bulk.partial_density().into_reduced(); - let mut bulk_density = self - .bulk - .eos - .component_index() - .iter() - .map(|&i| partial_density[i]) - .collect(); - - let (res, res_bulk, res_norm, _, _) = - self.euler_lagrange_equation(&density, &mut bulk_density, log)?; - Ok((res, res_bulk, res_norm)) + let (res, res_norm, _, _) = self.euler_lagrange_equation(&density, log)?; + Ok((res, res_norm)) } #[expect(clippy::type_complexity)] - pub(crate) fn euler_lagrange_equation( + fn fugacity( &self, density: &Array, - bulk_density: &mut Array1, - log: bool, - ) -> FeosResult<( - Array, - Array1, - f64, - Array, - Array, - )> { + ) -> FeosResult<(Array, Array, Array1)> { // calculate reduced temperature let temperature = self.temperature.to_reduced(); @@ -336,21 +341,10 @@ where // calculate total functional derivative dfdrho += &self.external_potential; - // calculate bulk functional derivative - let bulk_convolver = BulkConvolver::new(self.bulk.eos.weight_functions(temperature)); - let (_, dfdrho_bulk) = self.bulk.eos.functional_derivative( - temperature, - bulk_density, - bulk_convolver.as_ref(), - )?; dfdrho .outer_iter_mut() - .zip(dfdrho_bulk) .zip(self.bulk.eos.m().iter()) - .for_each(|((mut df, df_b), &m)| { - df -= df_b; - df /= m - }); + .for_each(|(mut df, &m)| df /= m); // calculate bond integrals let exp_dfdrho = dfdrho.mapv(|x| (-x).exp()); @@ -358,22 +352,36 @@ where .bulk .eos .bond_integrals(temperature, &exp_dfdrho, self.convolver.as_ref()); - let mut rho_projected = &exp_dfdrho * bonds; + let z = &exp_dfdrho * bonds; - // calculate bulk density based on the given specification - let res_bulk = &*bulk_density - - self.specification.calculate_bulk_density( - bulk_density, - &self.integrate_reduced_comp(&rho_projected), - )?; - *bulk_density -= &res_bulk; + // calculate fugacity based on the given specification + let fugacity = self + .specification + .calculate_fugacity(&self.integrate_reduced_comp(&z))?; - // multiply bulk density + Ok((exp_dfdrho, z, fugacity)) + } + + #[expect(clippy::type_complexity)] + pub(crate) fn euler_lagrange_equation( + &self, + density: &Array, + log: bool, + ) -> FeosResult<( + Array, + f64, + Array, + Array, + )> { + // calculate functional derivatives and fugacity + let (exp_dfdrho, mut rho_projected, fugacity) = self.fugacity(density)?; + + // multiply fugacity rho_projected .outer_iter_mut() - .zip(bulk_density.iter()) - .for_each(|(mut x, &rho_b)| { - x *= rho_b; + .zip(fugacity.iter()) + .for_each(|(mut x, &f)| { + x *= f; }); // calculate residual @@ -394,7 +402,7 @@ where (density - &rho_projected).mapv(|x| x * x).sum().sqrt() / (res.len() as f64).sqrt(); if res_norm.is_finite() { - Ok((res, res_bulk, res_norm, exp_dfdrho, rho_projected)) + Ok((res, res_norm, exp_dfdrho, rho_projected)) } else { Err(FeosError::IterationFailed("Euler-Lagrange equation".into())) } @@ -405,25 +413,35 @@ where let solver = solver.cloned().unwrap_or_default(); // Read from profile - let component_index = self.bulk.eos.component_index().into_owned(); let mut density = self.density.to_reduced(); - let partial_density = self.bulk.partial_density().into_reduced(); - let mut bulk_density = component_index - .iter() - .map(|&i| partial_density[i]) - .collect(); // Call solver(s) - self.call_solver(&mut density, &mut bulk_density, &solver, debug)?; + self.call_solver(&mut density, &solver, debug)?; + + // Update bulk state + if !matches!(self.specification, DFTSpecification::ChemicalPotential(_)) { + // solve a bulk profile with the Newton solver + let mut bulk_profile = + DFTProfile::::new(Grid::Bulk, &self.bulk, None, None, None); + let (_, _, fugacity) = self.fugacity(&density)?; + bulk_profile.specification = DFTSpecification::ChemicalPotential(fugacity); + let solver = DFTSolver::new(None).newton(None, None, None, None); + bulk_profile.solve(Some(&solver), false)?; + + // create the state based on the results from the bulk profile + let component_index = self.bulk.eos.component_index(); + let mut partial_density = self.bulk.partial_density(); + bulk_profile + .density + .into_iter() + .enumerate() + .for_each(|(i, r)| partial_density.set(component_index[i], r)); + self.bulk = State::new_density(&self.bulk.eos, self.bulk.temperature, partial_density)?; + } // Update profile self.density = Density::from_reduced(density); - let mut partial_density = self.bulk.partial_density(); - bulk_density - .into_iter() - .enumerate() - .for_each(|(i, r)| partial_density.set(component_index[i], Density::from_reduced(r))); - self.bulk = State::new_density(&self.bulk.eos, self.bulk.temperature, partial_density)?; + Ok(()) } } diff --git a/crates/feos-dft/src/profile/properties.rs b/crates/feos-dft/src/profile/properties.rs index 1c0427b28..8e602b34b 100644 --- a/crates/feos-dft/src/profile/properties.rs +++ b/crates/feos-dft/src/profile/properties.rs @@ -1,11 +1,11 @@ #![allow(type_alias_bounds)] use super::DFTProfile; -use crate::convolver::{BulkConvolver, Convolver}; +use crate::convolver::{Convolver, PeriodicConvolver}; use crate::functional_contribution::FunctionalContribution; use crate::{ConvolverFFT, DFTSolverLog, HelmholtzEnergyFunctional, WeightFunctionInfo}; use feos_core::{Contributions, FeosResult, ReferenceSystem, Total, Verbosity}; use nalgebra::{DMatrix, DVector}; -use ndarray::{Array, Array1, Axis, Dimension, RemoveAxis}; +use ndarray::{Array, Array1, Axis, Dimension, Ix0, RemoveAxis}; use num_dual::{Dual64, DualNum}; use quantity::{ Density, Energy, Entropy, EntropyDensity, MolarEnergy, Moles, Pressure, Quantity, Temperature, @@ -298,17 +298,8 @@ where { fn density_derivative(&self, lhs: &Array) -> FeosResult> { let rho = self.density.to_reduced(); - let partial_density = self.bulk.partial_density().into_reduced(); - let mut rho_bulk = self - .bulk - .eos - .component_index() - .iter() - .map(|&i| partial_density[i]) - .collect(); - let second_partial_derivatives = self.second_partial_derivatives(&rho)?; - let (_, _, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, &mut rho_bulk, false)?; + let (_, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, false)?; let rhs = |x: &_| { let delta_functional_derivative = @@ -411,7 +402,7 @@ where .map(|&i| partial_density[i]) .collect(); let rho_bulk_dual = rho_bulk.mapv(Dual64::from); - let bulk_convolver = BulkConvolver::new(weight_functions); + let bulk_convolver = PeriodicConvolver::<_, Ix0>::new_0d(&weight_functions); let (_, dfdrho_bulk) = self.bulk .eos diff --git a/crates/feos-dft/src/solver.rs b/crates/feos-dft/src/solver.rs index 69d2e2757..2baba60fb 100644 --- a/crates/feos-dft/src/solver.rs +++ b/crates/feos-dft/src/solver.rs @@ -216,7 +216,6 @@ where pub(crate) fn call_solver( &mut self, rho: &mut Array, - rho_bulk: &mut Array1, solver: &DFTSolver, debug: bool, ) -> FeosResult<()> { @@ -225,13 +224,11 @@ where let mut log = DFTSolverLog::new(solver.verbosity); for algorithm in &solver.algorithms { let (conv, iter) = match algorithm { - DFTAlgorithm::PicardIteration(picard) => { - self.solve_picard(*picard, rho, rho_bulk, &mut log) - } + DFTAlgorithm::PicardIteration(picard) => self.solve_picard(*picard, rho, &mut log), DFTAlgorithm::AndersonMixing(anderson) => { - self.solve_anderson(*anderson, rho, rho_bulk, &mut log) + self.solve_anderson(*anderson, rho, &mut log) } - DFTAlgorithm::Newton(newton) => self.solve_newton(*newton, rho, rho_bulk, &mut log), + DFTAlgorithm::Newton(newton) => self.solve_newton(*newton, rho, &mut log), }?; converged = conv; iterations += iter; @@ -255,7 +252,6 @@ where &self, picard: PicardIteration, rho: &mut Array, - rho_bulk: &mut Array1, log: &mut DFTSolverLog, ) -> FeosResult<(bool, usize)> { let solver = if picard.log { @@ -266,8 +262,7 @@ where for k in 0..picard.max_iter { // calculate residual - let (res, _, res_norm, _, _) = - self.euler_lagrange_equation(&*rho, rho_bulk, picard.log)?; + let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, picard.log)?; log.add_residual(solver, k, res_norm); // check for convergence @@ -276,10 +271,9 @@ where } // apply line search or constant damping - let damping_coefficient = picard.damping_coefficient.map_or_else( - || self.line_search(rho, &res, rho_bulk, res_norm, picard.log), - Ok, - )?; + let damping_coefficient = picard + .damping_coefficient + .map_or_else(|| self.line_search(rho, &res, res_norm, picard.log), Ok)?; // update solution if picard.log { @@ -295,7 +289,6 @@ where &self, rho: &Array, delta_rho: &Array, - rho_bulk: &mut Array1, res0: f64, logarithm: bool, ) -> FeosResult { @@ -311,9 +304,7 @@ where } else { rho + alpha * delta_rho }; - let Ok((_, _, res2, _, _)) = - self.euler_lagrange_equation(&rho_new, rho_bulk, logarithm) - else { + let Ok((_, res2, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { continue; }; if res2 > res0 { @@ -326,9 +317,7 @@ where } else { rho + 0.5 * alpha * delta_rho }; - let Ok((_, _, res1, _, _)) = - self.euler_lagrange_equation(&rho_new, rho_bulk, logarithm) - else { + let Ok((_, res1, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { continue; }; @@ -359,7 +348,6 @@ where &self, anderson: AndersonMixing, rho: &mut Array, - rho_bulk: &mut Array1, log: &mut DFTSolverLog, ) -> FeosResult<(bool, usize)> { let solver = if anderson.log { @@ -382,8 +370,7 @@ where let m = resm.len() + 1; // calculate residual - let (res, _, res_norm, _, _) = - self.euler_lagrange_equation(&*rho, rho_bulk, anderson.log)?; + let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, anderson.log)?; log.add_residual(solver, k, res_norm); // check for convergence @@ -434,14 +421,13 @@ where &self, newton: Newton, rho: &mut Array, - rho_bulk: &mut Array1, log: &mut DFTSolverLog, ) -> FeosResult<(bool, usize)> { let solver = if newton.log { "Newton (log)" } else { "Newton" }; for k in 0..newton.max_iter { // calculate initial residual - let (res, _, res_norm, exp_dfdrho, rho_p) = - self.euler_lagrange_equation(rho, rho_bulk, newton.log)?; + let (res, res_norm, exp_dfdrho, rho_p) = + self.euler_lagrange_equation(rho, newton.log)?; log.add_residual(solver, k, res_norm); // check convergence diff --git a/py-feos/src/dft/adsorption/pore.rs b/py-feos/src/dft/adsorption/pore.rs index ab66684e2..d07c75b39 100644 --- a/py-feos/src/dft/adsorption/pore.rs +++ b/py-feos/src/dft/adsorption/pore.rs @@ -28,6 +28,13 @@ impl PyGrid { Self(Grid::Cartesian1(x)) } + /// Generate a 1D Cartesian grid with periodic boundary conditions on both sides. + #[staticmethod] + pub fn periodical_1d(n_points: usize, length: Length) -> Self { + let x = AxisDFT::new_cartesian(n_points, length, None); + Self(Grid::Periodical1(x)) + } + /// Generate a polar grid with radial axis. #[staticmethod] pub fn polar(n_points: usize, length: Length) -> Self { diff --git a/py-feos/src/dft/profile.rs b/py-feos/src/dft/profile.rs index 40104de89..e1acef783 100644 --- a/py-feos/src/dft/profile.rs +++ b/py-feos/src/dft/profile.rs @@ -18,9 +18,9 @@ macro_rules! impl_profile { &self, log: bool, py: Python<'py>, - ) -> PyResult<(Bound<'py, PyArrayDyn>, Bound<'py, PyArray1>, f64)> { - let (res_rho, res_mu, res_norm) = self.0.profile.residual(log).map_err(PyFeosError::from)?; - Ok((res_rho.view().into_dyn().to_pyarray(py), res_mu.view().to_pyarray(py), res_norm)) + ) -> PyResult<(Bound<'py, PyArrayDyn>, f64)> { + let (res_rho, res_norm) = self.0.profile.residual(log).map_err(PyFeosError::from)?; + Ok((res_rho.view().into_dyn().to_pyarray(py), res_norm)) } /// Solve the profile in-place. A non-default solver can be provided From 3a06cc61f4685f1cb586910b53de184b33a00769 Mon Sep 17 00:00:00 2001 From: Philipp Rehner Date: Sun, 26 Jul 2026 12:22:01 +0200 Subject: [PATCH 2/2] Use rayon to parallelize the evaluation of the functional --- crates/feos-dft/Cargo.toml | 2 +- crates/feos-dft/src/functional.rs | 21 +++++- .../feos-dft/src/functional_contribution.rs | 43 ++++++++++- crates/feos-dft/src/profile/mod.rs | 73 ++++++++++++------- crates/feos-dft/src/profile/properties.rs | 19 +++-- crates/feos-dft/src/solver.rs | 59 +++++++++++---- py-feos/src/dft/profile.rs | 2 +- py-feos/src/dft/solver.rs | 17 ++++- 8 files changed, 180 insertions(+), 56 deletions(-) diff --git a/crates/feos-dft/Cargo.toml b/crates/feos-dft/Cargo.toml index 400504e9e..67fd4c129 100644 --- a/crates/feos-dft/Cargo.toml +++ b/crates/feos-dft/Cargo.toml @@ -16,7 +16,7 @@ features = ["rayon"] [dependencies] quantity = { workspace = true, features = ["ndarray"] } -num-dual = { workspace = true } +num-dual = { workspace = true, features = ["ndarray"] } ndarray = { workspace = true } nalgebra = { workspace = true } rustdct = { workspace = true } diff --git a/crates/feos-dft/src/functional.rs b/crates/feos-dft/src/functional.rs index f584f2110..715af7511 100644 --- a/crates/feos-dft/src/functional.rs +++ b/crates/feos-dft/src/functional.rs @@ -11,6 +11,7 @@ use petgraph::graph::{Graph, UnGraph}; use petgraph::visit::EdgeRef; use std::borrow::Cow; use std::ops::{Deref, MulAssign}; +use std::time::{Duration, Instant}; impl HelmholtzEnergyFunctionalDyn for EquationOfState, F> @@ -112,17 +113,23 @@ pub trait HelmholtzEnergyFunctional: Residual { /// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{\beta F}}{\delta\rho_i(\mathbf{r})}$. #[expect(clippy::type_complexity)] - fn functional_derivative + Copy>( + fn functional_derivative>( &self, temperature: N, density: &Array, convolver: &dyn Convolver, - ) -> FeosResult<(Array, Array)> + ) -> FeosResult<(Array, Array, [Duration; 3])> where D: Dimension, D::Larger: Dimension, { + // calculate weighted densities + let start = Instant::now(); let weighted_densities = convolver.weighted_densities(density); + let t_wd = start.elapsed(); + + // calculate partial derivatives + let start = Instant::now(); let contributions = self.contributions(); let mut partial_derivatives = Vec::new(); let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0))); @@ -140,9 +147,17 @@ pub trait HelmholtzEnergyFunctional: Residual { partial_derivatives.push(pd); helmholtz_energy_density += φ } + let t_pd = start.elapsed(); + + // calculate functional derivative + let start = Instant::now(); + let functional_derivative = convolver.functional_derivative(&partial_derivatives); + let t_fd = start.elapsed(); + Ok(( helmholtz_energy_density, - convolver.functional_derivative(&partial_derivatives), + functional_derivative, + [t_wd, t_pd, t_fd], )) } diff --git a/crates/feos-dft/src/functional_contribution.rs b/crates/feos-dft/src/functional_contribution.rs index ff9e6f9a6..45879c181 100644 --- a/crates/feos-dft/src/functional_contribution.rs +++ b/crates/feos-dft/src/functional_contribution.rs @@ -1,10 +1,17 @@ use crate::weight_functions::WeightFunctionInfo; use feos_core::{FeosResult, StateHD}; use ndarray::RemoveAxis; +#[cfg(feature = "rayon")] +use ndarray::parallel::prelude::*; use ndarray::prelude::*; use num_dual::*; use num_traits::Zero; +// The chunk size is specifically chosen to be the default amount of grid points in a 1D density +// profile to switch to parallel evaluation only for 2D and 3D profiles. +#[cfg(feature = "rayon")] +const CHUNK_SIZE: usize = 2048; + /// Individual functional contribution that can be evaluated using generalized (hyper) dual numbers. pub trait FunctionalContribution: Sync + Send { /// Return the name of the contribution. @@ -31,6 +38,36 @@ pub trait FunctionalContribution: Sync + Send { weighted_densities: ArrayView2, ) -> FeosResult>; + fn helmholtz_energy_density_parallel>( + &self, + temperature: N, + weighted_densities: ArrayView2, + ) -> FeosResult> { + #[cfg(feature = "rayon")] + { + let n = weighted_densities.shape()[1]; + if n > CHUNK_SIZE { + let mut phi = Array::zeros(n); + let mut phi_view = phi.view_mut(); + let wd_iter = weighted_densities + .axis_chunks_iter(Axis(1), CHUNK_SIZE) + .into_par_iter(); + let phi_iter = phi_view + .axis_chunks_iter_mut(Axis(0), CHUNK_SIZE) + .into_par_iter(); + wd_iter.zip(phi_iter).try_for_each(|(wd, mut phi)| { + self.helmholtz_energy_density_parallel(temperature, wd.view()) + .map(|p| phi.assign(&p)) + })?; + Ok(phi) + } else { + self.helmholtz_energy_density(temperature, weighted_densities) + } + } + #[cfg(not(feature = "rayon"))] + self.helmholtz_energy_density(temperature, weighted_densities) + } + fn bulk_helmholtz_energy_density + Copy>( &self, state: &StateHD, @@ -52,7 +89,7 @@ pub trait FunctionalContribution: Sync + Send { .unwrap()[0] } - fn first_partial_derivatives + Copy>( + fn first_partial_derivatives>( &self, temperature: N, weighted_densities: Array2, @@ -66,7 +103,7 @@ pub trait FunctionalContribution: Sync + Send { for i in 0..wd.shape()[0] { wd.index_axis_mut(Axis(0), i) .map_inplace(|x| x.eps = N::one()); - phi = self.helmholtz_energy_density(t, wd.view())?; + phi = self.helmholtz_energy_density_parallel(t, wd.view())?; first_partial_derivative .index_axis_mut(Axis(0), i) .assign(&phi.mapv(|p| p.eps)); @@ -93,7 +130,7 @@ pub trait FunctionalContribution: Sync + Send { wd.index_axis_mut(Axis(0), i).map_inplace(|x| x.eps1 = 1.0); for j in 0..=i { wd.index_axis_mut(Axis(0), j).map_inplace(|x| x.eps2 = 1.0); - phi = self.helmholtz_energy_density(t, wd.view())?; + phi = self.helmholtz_energy_density_parallel(t, wd.view())?; let p = phi.mapv(|p| p.eps1eps2); second_partial_derivative .index_axis_mut(Axis(0), i) diff --git a/crates/feos-dft/src/profile/mod.rs b/crates/feos-dft/src/profile/mod.rs index 6def8ded2..c9085281f 100644 --- a/crates/feos-dft/src/profile/mod.rs +++ b/crates/feos-dft/src/profile/mod.rs @@ -9,6 +9,7 @@ use num_dual::DualNum; use quantity::{_Volume, Density, Energy, Entropy, Length, Moles, Quantity, Temperature, Volume}; use std::ops::{Add, MulAssign}; use std::sync::Arc; +use std::time::{Duration, Instant}; mod properties; @@ -34,14 +35,24 @@ pub enum DFTSpecification { } impl DFTSpecification { - fn calculate_fugacity(&self, z: &Array1) -> FeosResult> { - Ok(match self { + fn calculate_fugacity(&self, z: &Array1) -> Array1 { + match self { Self::ChemicalPotential(fugacity) => fugacity.clone(), Self::Moles(moles) => moles / z, Self::TotalMoles(total_moles, fugacity) => { fugacity * *total_moles / (fugacity * z).sum() } - }) + } + } + + pub(crate) fn delta_fugacity(&self, z: &Array1, delta_z: &Array1) -> Array1 { + match self { + Self::ChemicalPotential(fugacity) => Array1::zeros(fugacity.len()), + Self::Moles(_) => -delta_z / z, + Self::TotalMoles(_, fugacity) => { + -(fugacity * delta_z).sum() / (fugacity * z).sum() * Array1::ones(fugacity.len()) + } + } } pub fn from_state(state: &State) -> Self { @@ -55,7 +66,7 @@ impl DFTSpecification { .collect(); let bulk_convolver = PeriodicConvolver::<_, Ix0>::new_0d(&state.eos.weight_functions(temperature)); - let (_, dfdrho_bulk) = state + let (_, dfdrho_bulk, _) = state .eos .functional_derivative(temperature, &bulk_density, bulk_convolver.as_ref()) .unwrap(); @@ -216,7 +227,7 @@ where profile.sum() * functional_determinant } - fn integrate_reduced_comp, N: DualNum + Copy>( + pub(crate) fn integrate_reduced_comp, N: DualNum + Copy>( &self, profile: &ArrayBase, ) -> Array1 { @@ -306,7 +317,7 @@ where } } -impl DFTProfile +impl DFTProfile where D::Larger: Dimension, ::Larger: Dimension, @@ -320,20 +331,26 @@ where pub fn residual(&self, log: bool) -> FeosResult<(Array, f64)> { let density = self.density.to_reduced(); - let (res, res_norm, _, _) = self.euler_lagrange_equation(&density, log)?; + let (res, res_norm, _, _, _, _) = self.euler_lagrange_equation(&density, log)?; Ok((res, res_norm)) } #[expect(clippy::type_complexity)] - fn fugacity( + pub(crate) fn fugacity( &self, density: &Array, - ) -> FeosResult<(Array, Array, Array1)> { + ) -> FeosResult<( + Array, + Array1, + Array, + Array1, + [Duration; 4], + )> { // calculate reduced temperature let temperature = self.temperature.to_reduced(); // calculate intrinsic functional derivative - let (_, mut dfdrho) = + let (_, mut dfdrho, [t0, t1, t2]) = self.bulk .eos .functional_derivative(temperature, density, self.convolver.as_ref())?; @@ -341,6 +358,7 @@ where // calculate total functional derivative dfdrho += &self.external_potential; + let start = Instant::now(); dfdrho .outer_iter_mut() .zip(self.bulk.eos.m().iter()) @@ -352,14 +370,23 @@ where .bulk .eos .bond_integrals(temperature, &exp_dfdrho, self.convolver.as_ref()); - let z = &exp_dfdrho * bonds; + let mut rho_projected = &exp_dfdrho * bonds; + let z = self.integrate_reduced_comp(&rho_projected); // calculate fugacity based on the given specification - let fugacity = self - .specification - .calculate_fugacity(&self.integrate_reduced_comp(&z))?; + let fugacity = self.specification.calculate_fugacity(&z); - Ok((exp_dfdrho, z, fugacity)) + // multiply fugacity + rho_projected + .outer_iter_mut() + .zip(fugacity.iter()) + .for_each(|(mut x, &f)| { + x *= f; + }); + + let t3 = start.elapsed(); + + Ok((exp_dfdrho, z, rho_projected, fugacity, [t0, t1, t2, t3])) } #[expect(clippy::type_complexity)] @@ -371,18 +398,12 @@ where Array, f64, Array, + Array1, Array, + [Duration; 4], )> { // calculate functional derivatives and fugacity - let (exp_dfdrho, mut rho_projected, fugacity) = self.fugacity(density)?; - - // multiply fugacity - rho_projected - .outer_iter_mut() - .zip(fugacity.iter()) - .for_each(|(mut x, &f)| { - x *= f; - }); + let (exp_dfdrho, z, rho_projected, _, timings) = self.fugacity(density)?; // calculate residual let mut res = if log { @@ -402,7 +423,7 @@ where (density - &rho_projected).mapv(|x| x * x).sum().sqrt() / (res.len() as f64).sqrt(); if res_norm.is_finite() { - Ok((res, res_norm, exp_dfdrho, rho_projected)) + Ok((res, res_norm, exp_dfdrho, z, rho_projected, timings)) } else { Err(FeosError::IterationFailed("Euler-Lagrange equation".into())) } @@ -423,7 +444,7 @@ where // solve a bulk profile with the Newton solver let mut bulk_profile = DFTProfile::::new(Grid::Bulk, &self.bulk, None, None, None); - let (_, _, fugacity) = self.fugacity(&density)?; + let (_, _, _, fugacity, _) = self.fugacity(&density)?; bulk_profile.specification = DFTSpecification::ChemicalPotential(fugacity); let solver = DFTSolver::new(None).newton(None, None, None, None); bulk_profile.solve(Some(&solver), false)?; diff --git a/crates/feos-dft/src/profile/properties.rs b/crates/feos-dft/src/profile/properties.rs index 8e602b34b..a09dbe034 100644 --- a/crates/feos-dft/src/profile/properties.rs +++ b/crates/feos-dft/src/profile/properties.rs @@ -12,6 +12,7 @@ use quantity::{ }; use std::ops::{AddAssign, Div}; use std::sync::Arc; +use std::time::Duration; type DrhoDmu = ::Larger>> as Div>::Output; @@ -30,7 +31,7 @@ where // Calculate residual Helmholtz energy density and functional derivative let t = self.temperature.to_reduced(); let rho = self.density.to_reduced(); - let (mut f, dfdrho) = + let (mut f, dfdrho, _) = self.bulk .eos .functional_derivative(t, &rho, self.convolver.as_ref())?; @@ -59,13 +60,15 @@ where } /// Calculate the (residual) intrinsic functional derivative $\frac{\delta\mathcal{F}}{\delta\rho_i(\mathbf{r})}$. - pub fn functional_derivative(&self) -> FeosResult> { - let (_, dfdrho) = self.bulk.eos.functional_derivative( + #[expect(clippy::type_complexity)] + pub fn functional_derivative( + &self, + ) -> FeosResult<(Array, Array, [Duration; 3])> { + self.bulk.eos.functional_derivative( self.temperature.to_reduced(), &self.density.to_reduced(), self.convolver.as_ref(), - )?; - Ok(dfdrho) + ) } } @@ -299,7 +302,7 @@ where fn density_derivative(&self, lhs: &Array) -> FeosResult> { let rho = self.density.to_reduced(); let second_partial_derivatives = self.second_partial_derivatives(&rho)?; - let (_, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, false)?; + let (_, _, exp_dfdrho, _, _, _) = self.euler_lagrange_equation(&rho, false)?; let rhs = |x: &_| { let delta_functional_derivative = @@ -384,7 +387,7 @@ where .collect(); let convolver: Arc> = ConvolverFFT::plan(&self.grid, &weight_functions, self.lanczos); - let (_, mut dfdrho) = + let (_, mut dfdrho, _) = self.bulk .eos .functional_derivative(t_dual, &rho_dual, convolver.as_ref())?; @@ -403,7 +406,7 @@ where .collect(); let rho_bulk_dual = rho_bulk.mapv(Dual64::from); let bulk_convolver = PeriodicConvolver::<_, Ix0>::new_0d(&weight_functions); - let (_, dfdrho_bulk) = + let (_, dfdrho_bulk, _) = self.bulk .eos .functional_derivative(t_dual, &rho_bulk_dual, bulk_convolver.as_ref())?; diff --git a/crates/feos-dft/src/solver.rs b/crates/feos-dft/src/solver.rs index 2baba60fb..1269a367a 100644 --- a/crates/feos-dft/src/solver.rs +++ b/crates/feos-dft/src/solver.rs @@ -160,6 +160,10 @@ pub struct DFTSolverLog { residual: Vec, time: Vec, solver: Vec<&'static str>, + pub time_weighted_densities: Duration, + pub time_partial_derivatives: Duration, + pub time_functional_derivative: Duration, + time_euler_lagrange_equation: Duration, } impl DFTSolverLog { @@ -174,10 +178,20 @@ impl DFTSolverLog { residual: Vec::new(), time: Vec::new(), solver: Vec::new(), + time_weighted_densities: Duration::ZERO, + time_partial_derivatives: Duration::ZERO, + time_functional_derivative: Duration::ZERO, + time_euler_lagrange_equation: Duration::ZERO, } } - fn add_residual(&mut self, solver: &'static str, iteration: usize, residual: f64) { + fn add_residual( + &mut self, + solver: &'static str, + iteration: usize, + residual: f64, + [t_wd, t_pd, t_fd, t_el]: [Duration; 4], + ) { if iteration == 0 { log_iter!(self.verbosity, "{:-<59}", ""); } @@ -193,6 +207,10 @@ impl DFTSolverLog { time.as_secs_f64() * SECOND, residual, ); + self.time_weighted_densities += t_wd; + self.time_partial_derivatives += t_pd; + self.time_functional_derivative += t_fd; + self.time_euler_lagrange_equation += t_el; } pub fn residual(&self) -> ArrayView1<'_, f64> { @@ -208,7 +226,7 @@ impl DFTSolverLog { } } -impl DFTProfile +impl DFTProfile where D::Larger: Dimension, ::Larger: Dimension, @@ -262,8 +280,9 @@ where for k in 0..picard.max_iter { // calculate residual - let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, picard.log)?; - log.add_residual(solver, k, res_norm); + let (res, res_norm, _, _, _, timings) = + self.euler_lagrange_equation(&*rho, picard.log)?; + log.add_residual(solver, k, res_norm, timings); // check for convergence if res_norm < picard.tol { @@ -304,7 +323,8 @@ where } else { rho + alpha * delta_rho }; - let Ok((_, res2, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { + let Ok((_, res2, _, _, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) + else { continue; }; if res2 > res0 { @@ -317,7 +337,8 @@ where } else { rho + 0.5 * alpha * delta_rho }; - let Ok((_, res1, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) else { + let Ok((_, res1, _, _, _, _)) = self.euler_lagrange_equation(&rho_new, logarithm) + else { continue; }; @@ -370,8 +391,9 @@ where let m = resm.len() + 1; // calculate residual - let (res, res_norm, _, _) = self.euler_lagrange_equation(&*rho, anderson.log)?; - log.add_residual(solver, k, res_norm); + let (res, res_norm, _, _, _, timings) = + self.euler_lagrange_equation(&*rho, anderson.log)?; + log.add_residual(solver, k, res_norm, timings); // check for convergence if res_norm < anderson.tol { @@ -426,9 +448,9 @@ where let solver = if newton.log { "Newton (log)" } else { "Newton" }; for k in 0..newton.max_iter { // calculate initial residual - let (res, res_norm, exp_dfdrho, rho_p) = + let (res, res_norm, exp_dfdrho, z, rho_p, timings) = self.euler_lagrange_equation(rho, newton.log)?; - log.add_residual(solver, k, res_norm); + log.add_residual(solver, k, res_norm, timings); // check convergence if res_norm < newton.tol { @@ -447,8 +469,19 @@ where .zip(self.bulk.eos.m().iter()) .for_each(|(mut q, &m)| q /= m); let delta_i = self.delta_bond_integrals(&exp_dfdrho, &delta_functional_derivative); + let mut delta_exp_dfdrho = delta_functional_derivative - delta_i; + let delta_z = -self.integrate_reduced_comp(&(&delta_exp_dfdrho * &exp_dfdrho)); + + let delta_fugacity = self.specification.delta_fugacity(&z, &delta_z); + delta_exp_dfdrho + .outer_iter_mut() + .zip(delta_fugacity.iter()) + .for_each(|(mut z, &f)| { + z -= f; + }); + let rho = if newton.log { &*rho } else { &rho_p }; - delta_rho + (delta_functional_derivative - delta_i) * rho + delta_rho + delta_exp_dfdrho * rho }; // update solution @@ -479,7 +512,7 @@ where gamma[0] = (r0 * r0).sum().sqrt(); v.push(r0 / gamma[0]); - log.add_residual("GMRES", 0, gamma[0]); + log.add_residual("GMRES", 0, gamma[0], [Duration::ZERO; 4]); let mut iter = 0; for j in 0..max_iter { @@ -515,7 +548,7 @@ where gamma[j] *= c[j + 1]; // check for convergence - log.add_residual("GMRES", j + 1, gamma[j + 1].abs()); + log.add_residual("GMRES", j + 1, gamma[j + 1].abs(), [Duration::ZERO; 4]); if gamma[j + 1].abs() >= tol && j + 1 < max_iter { v.push(q / h[(j + 1, j)]); iter += 1; diff --git a/py-feos/src/dft/profile.rs b/py-feos/src/dft/profile.rs index e1acef783..0a91378fc 100644 --- a/py-feos/src/dft/profile.rs +++ b/py-feos/src/dft/profile.rs @@ -114,7 +114,7 @@ macro_rules! impl_profile { &self, py: Python<'py>, ) -> PyResult>> { - Ok(self.0.profile.functional_derivative().map_err(PyFeosError::from)?.view().into_dyn().to_pyarray(py)) + Ok(self.0.profile.functional_derivative().map_err(PyFeosError::from)?.1.view().into_dyn().to_pyarray(py)) } /// Calculate the entropy density of the inhomogeneous system. diff --git a/py-feos/src/dft/solver.rs b/py-feos/src/dft/solver.rs index 4ee5e7410..05990dff6 100644 --- a/py-feos/src/dft/solver.rs +++ b/py-feos/src/dft/solver.rs @@ -3,7 +3,7 @@ use feos_dft::{DFTSolver, DFTSolverLog}; use ndarray::Array1; use numpy::{PyArray1, ToPyArray}; use pyo3::prelude::*; -use quantity::Time; +use quantity::{SECOND, Time}; /// Settings for the DFT solver. /// @@ -180,4 +180,19 @@ impl PyDFTSolverLog { fn get_solver(&self) -> Vec<&'static str> { self.0.solver().to_vec() } + + #[getter] + fn get_time_weighted_densities(&self) -> Time { + self.0.time_weighted_densities.as_secs_f64() * SECOND + } + + #[getter] + fn get_time_partial_derivatives(&self) -> Time { + self.0.time_partial_derivatives.as_secs_f64() * SECOND + } + + #[getter] + fn get_time_functional_derivative(&self) -> Time { + self.0.time_functional_derivative.as_secs_f64() * SECOND + } }