feos_dft/profile/
mod.rs

1use crate::convolver::{BulkConvolver, Convolver, ConvolverFFT};
2use crate::functional::HelmholtzEnergyFunctional;
3use crate::geometry::Grid;
4use crate::solver::{DFTSolver, DFTSolverLog};
5use feos_core::{FeosError, FeosResult, ReferenceSystem, State};
6use nalgebra::{DVector, Dyn, U1};
7use ndarray::{
8    Array, Array1, Array2, Array3, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix1, Ix2, Ix3,
9    RemoveAxis,
10};
11use num_dual::DualNum;
12use quantity::{_Volume, DEGREES, Density, Length, Moles, Quantity, Temperature, Volume};
13use std::ops::{Add, MulAssign};
14use std::sync::Arc;
15
16mod properties;
17
18pub(crate) const MAX_POTENTIAL: f64 = 50.0;
19#[cfg(feature = "rayon")]
20pub(crate) const CUTOFF_RADIUS: f64 = 14.0;
21
22/// General specifications for the chemical potential in a DFT calculation.
23///
24/// In the most basic case, the chemical potential is specified in a DFT calculation,
25/// for more general systems, this trait provides the possibility to declare additional
26/// equations for the calculation of the chemical potential during the iteration.
27pub trait DFTSpecification<D: Dimension, F>: Send + Sync {
28    fn calculate_bulk_density(
29        &self,
30        profile: &DFTProfile<D, F>,
31        bulk_density: &Array1<f64>,
32        z: &Array1<f64>,
33    ) -> FeosResult<Array1<f64>>;
34}
35
36/// Common specifications for the grand potentials in a DFT calculation.
37pub enum DFTSpecifications {
38    /// DFT with specified chemical potential.
39    ChemicalPotential,
40    /// DFT with specified number of particles.
41    ///
42    /// The solution is still a grand canonical density profile, but the chemical
43    /// potentials are iterated together with the density profile to obtain a result
44    /// with the specified number of particles.
45    Moles { moles: Array1<f64> },
46    /// DFT with specified total number of moles.
47    TotalMoles { total_moles: f64 },
48}
49
50impl DFTSpecifications {
51    /// Calculate the number of particles from the profile.
52    ///
53    /// Call this after initializing the density profile to keep the number of
54    /// particles constant in systems, where the number itself is difficult to obtain.
55    pub fn moles_from_profile<D: Dimension, F: HelmholtzEnergyFunctional>(
56        profile: &DFTProfile<D, F>,
57    ) -> Self
58    where
59        D::Larger: Dimension<Smaller = D>,
60    {
61        let rho = profile.density.to_reduced();
62        Self::Moles {
63            moles: profile.integrate_reduced_comp(&rho),
64        }
65    }
66
67    /// Calculate the number of particles from the profile.
68    ///
69    /// Call this after initializing the density profile to keep the total number of
70    /// particles constant in systems, e.g. to fix the equimolar dividing surface.
71    pub fn total_moles_from_profile<D: Dimension, F: HelmholtzEnergyFunctional>(
72        profile: &DFTProfile<D, F>,
73    ) -> Self
74    where
75        D::Larger: Dimension<Smaller = D>,
76    {
77        let rho = profile.density.to_reduced();
78        let moles = profile.integrate_reduced_comp(&rho).sum();
79        Self::TotalMoles { total_moles: moles }
80    }
81}
82
83impl<D: Dimension, F: HelmholtzEnergyFunctional> DFTSpecification<D, F> for DFTSpecifications {
84    fn calculate_bulk_density(
85        &self,
86        _profile: &DFTProfile<D, F>,
87        bulk_density: &Array1<f64>,
88        z: &Array1<f64>,
89    ) -> FeosResult<Array1<f64>> {
90        Ok(match self {
91            Self::ChemicalPotential => bulk_density.clone(),
92            Self::Moles { moles } => moles / z,
93            Self::TotalMoles { total_moles } => {
94                bulk_density * *total_moles / (bulk_density * z).sum()
95            }
96        })
97    }
98}
99
100#[derive(Clone)]
101/// A one-, two-, or three-dimensional density profile.
102pub struct DFTProfile<D: Dimension, F> {
103    pub grid: Grid,
104    pub convolver: Arc<dyn Convolver<f64, D>>,
105    pub temperature: Temperature,
106    pub density: Density<Array<f64, D::Larger>>,
107    pub specification: Arc<dyn DFTSpecification<D, F>>,
108    pub external_potential: Array<f64, D::Larger>,
109    pub bulk: State<F>,
110    pub solver_log: Option<DFTSolverLog>,
111    pub lanczos: Option<i32>,
112}
113
114impl<F> DFTProfile<Ix1, F> {
115    pub fn r(&self) -> Length<Array1<f64>> {
116        Length::from_reduced(self.grid.grids()[0].to_owned())
117    }
118
119    pub fn z(&self) -> Length<Array1<f64>> {
120        Length::from_reduced(self.grid.grids()[0].to_owned())
121    }
122}
123
124impl<F> DFTProfile<Ix2, F> {
125    pub fn edges(&self) -> [Length<Array1<f64>>; 2] {
126        [
127            Length::from_reduced(self.grid.axes()[0].edges.to_owned()),
128            Length::from_reduced(self.grid.axes()[1].edges.to_owned()),
129        ]
130    }
131
132    pub fn meshgrid(&self) -> [Length<Array2<f64>>; 2] {
133        let (u, v, alpha) = match &self.grid {
134            Grid::Cartesian2(u, v) => (u, v, 90.0 * DEGREES),
135            Grid::Periodical2(u, v, alpha) => (u, v, *alpha),
136            _ => unreachable!(),
137        };
138        let u_grid = Array::from_shape_fn([u.grid.len(), v.grid.len()], |(i, _)| u.grid[i]);
139        let v_grid = Array::from_shape_fn([u.grid.len(), v.grid.len()], |(_, j)| v.grid[j]);
140        let x = Length::from_reduced(u_grid + &v_grid * alpha.cos());
141        let y = Length::from_reduced(v_grid * alpha.sin());
142        [x, y]
143    }
144
145    pub fn r(&self) -> Length<Array1<f64>> {
146        Length::from_reduced(self.grid.grids()[0].to_owned())
147    }
148
149    pub fn z(&self) -> Length<Array1<f64>> {
150        Length::from_reduced(self.grid.grids()[1].to_owned())
151    }
152}
153
154impl<F> DFTProfile<Ix3, F> {
155    pub fn edges(&self) -> [Length<Array1<f64>>; 3] {
156        [
157            Length::from_reduced(self.grid.axes()[0].edges.to_owned()),
158            Length::from_reduced(self.grid.axes()[1].edges.to_owned()),
159            Length::from_reduced(self.grid.axes()[2].edges.to_owned()),
160        ]
161    }
162
163    pub fn meshgrid(&self) -> [Length<Array3<f64>>; 3] {
164        let (u, v, w, [alpha, beta, gamma]) = match &self.grid {
165            Grid::Cartesian3(u, v, w) => (u, v, w, [90.0 * DEGREES; 3]),
166            Grid::Periodical3(u, v, w, angles) => (u, v, w, *angles),
167            _ => unreachable!(),
168        };
169        let shape = [u.grid.len(), v.grid.len(), w.grid.len()];
170        let u_grid = Array::from_shape_fn(shape, |(i, _, _)| u.grid[i]);
171        let v_grid = Array::from_shape_fn(shape, |(_, j, _)| v.grid[j]);
172        let w_grid = Array::from_shape_fn(shape, |(_, _, k)| w.grid[k]);
173        let xi = (alpha.cos() - gamma.cos() * beta.cos()) / gamma.sin();
174        let zeta = (1.0_f64 - beta.cos().powi(2) - xi * xi).sqrt();
175        let x = Length::from_reduced(u_grid + &v_grid * gamma.cos() + &w_grid * beta.cos());
176        let y = Length::from_reduced(v_grid * gamma.sin() + &w_grid * xi);
177        let z = Length::from_reduced(w_grid * zeta);
178        [x, y, z]
179    }
180
181    pub fn x(&self) -> Length<Array1<f64>> {
182        Length::from_reduced(self.grid.grids()[0].to_owned())
183    }
184
185    pub fn y(&self) -> Length<Array1<f64>> {
186        Length::from_reduced(self.grid.grids()[1].to_owned())
187    }
188
189    pub fn z(&self) -> Length<Array1<f64>> {
190        Length::from_reduced(self.grid.grids()[2].to_owned())
191    }
192}
193
194impl<D: Dimension + RemoveAxis + 'static, F: HelmholtzEnergyFunctional> DFTProfile<D, F>
195where
196    D::Larger: Dimension<Smaller = D>,
197    D::Smaller: Dimension<Larger = D>,
198    <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
199{
200    /// Create a new density profile.
201    ///
202    /// If no external potential is specified, it is set to 0. The density is
203    /// initialized based on the bulk state and the external potential. The
204    /// specification is set to `ChemicalPotential` and needs to be overriden
205    /// after this call if something else is required.
206    pub fn new(
207        grid: Grid,
208        bulk: &State<F>,
209        external_potential: Option<Array<f64, D::Larger>>,
210        density: Option<&Density<Array<f64, D::Larger>>>,
211        lanczos: Option<i32>,
212    ) -> Self {
213        // initialize convolver
214        let t = bulk.temperature.to_reduced();
215        let weight_functions = bulk.eos.weight_functions(t);
216        let convolver = ConvolverFFT::plan(&grid, &weight_functions, lanczos);
217
218        // initialize external potential
219        let external_potential = external_potential.unwrap_or_else(|| {
220            let mut n_grid = vec![bulk.eos.component_index().len()];
221            grid.axes()
222                .iter()
223                .for_each(|&ax| n_grid.push(ax.grid.len()));
224            Array::zeros(n_grid).into_dimensionality().unwrap()
225        });
226
227        // initialize density
228        let density = if let Some(density) = density {
229            density.to_owned()
230        } else {
231            let exp_dfdrho = (-&external_potential).mapv(f64::exp);
232            let mut bonds = bulk.eos.bond_integrals(t, &exp_dfdrho, convolver.as_ref());
233            bonds *= &exp_dfdrho;
234            let mut density = Array::zeros(external_potential.raw_dim());
235            let bulk_density = bulk.partial_density().into_reduced();
236            for (s, &c) in bulk.eos.component_index().iter().enumerate() {
237                density.index_axis_mut(Axis_nd(0), s).assign(
238                    &(bonds.index_axis(Axis_nd(0), s).map(|is| is.min(1.0)) * bulk_density[c]),
239                );
240            }
241            Density::from_reduced(density)
242        };
243
244        Self {
245            grid,
246            convolver,
247            temperature: bulk.temperature,
248            density,
249            specification: Arc::new(DFTSpecifications::ChemicalPotential),
250            external_potential,
251            bulk: bulk.clone(),
252            solver_log: None,
253            lanczos,
254        }
255    }
256}
257
258impl<D: Dimension, F: HelmholtzEnergyFunctional> DFTProfile<D, F>
259where
260    D::Larger: Dimension<Smaller = D>,
261{
262    fn integrate_reduced<N: DualNum<f64> + Copy>(&self, mut profile: Array<N, D>) -> N {
263        let (integration_weights, functional_determinant) = self.grid.integration_weights();
264
265        for (i, w) in integration_weights.into_iter().enumerate() {
266            for mut l in profile.lanes_mut(Axis_nd(i)) {
267                l.mul_assign(&w.mapv(N::from));
268            }
269        }
270        profile.sum() * functional_determinant
271    }
272
273    fn integrate_reduced_comp<S: Data<Elem = N>, N: DualNum<f64> + Copy>(
274        &self,
275        profile: &ArrayBase<S, D::Larger>,
276    ) -> Array1<N> {
277        Array1::from_shape_fn(profile.shape()[0], |i| {
278            self.integrate_reduced(profile.index_axis(Axis_nd(0), i).to_owned())
279        })
280    }
281
282    pub(crate) fn integrate_reduced_segments<S: Data<Elem = N>, N: DualNum<f64> + Copy>(
283        &self,
284        profile: &ArrayBase<S, D::Larger>,
285    ) -> DVector<N> {
286        let integral = self.integrate_reduced_comp(profile);
287        let mut integral_comp = DVector::zeros(self.bulk.eos.components());
288        for (i, &j) in self.bulk.eos.component_index().iter().enumerate() {
289            integral_comp[j] = integral[i];
290        }
291        integral_comp
292    }
293
294    /// Return the volume of the profile.
295    ///
296    /// In periodic directions, the length is assumed to be 1 Å.
297    pub fn volume(&self) -> Volume {
298        let volume: f64 = self.grid.axes().iter().map(|ax| ax.volume()).product();
299        Volume::from_reduced(volume * self.grid.functional_determinant())
300    }
301
302    /// Integrate a given profile over the iteration domain.
303    pub fn integrate<S: Data<Elem = f64>, U>(
304        &self,
305        profile: &Quantity<ArrayBase<S, D>, U>,
306    ) -> Quantity<f64, <_Volume as Add<U>>::Output>
307    where
308        _Volume: Add<U>,
309    {
310        let (integration_weights, functional_determinant) = self.grid.integration_weights();
311        let mut value = profile.to_owned();
312        for (i, &w) in integration_weights.iter().enumerate() {
313            for mut l in value.lanes_mut(Axis_nd(i)) {
314                l.mul_assign(w);
315            }
316        }
317        Volume::from_reduced(functional_determinant) * value.sum()
318    }
319
320    /// Integrate each component individually.
321    pub fn integrate_comp<S: Data<Elem = f64>, U>(
322        &self,
323        profile: &Quantity<ArrayBase<S, D::Larger>, U>,
324    ) -> Quantity<DVector<f64>, <_Volume as Add<U>>::Output>
325    where
326        _Volume: Add<U>,
327    {
328        Quantity::from_fn_generic(Dyn(profile.shape()[0]), U1, |i, _| {
329            self.integrate(&profile.index_axis(Axis_nd(0), i))
330        })
331    }
332
333    /// Integrate each segment individually and aggregate to components.
334    pub fn integrate_segments<S: Data<Elem = f64>, U>(
335        &self,
336        profile: &Quantity<ArrayBase<S, D::Larger>, U>,
337    ) -> Quantity<DVector<f64>, <_Volume as Add<U>>::Output>
338    where
339        _Volume: Add<U>,
340    {
341        let integral = self.integrate_comp(profile);
342        let mut integral_comp = Quantity::new(DVector::zeros(self.bulk.eos.components()));
343        for (i, &j) in self.bulk.eos.component_index().iter().enumerate() {
344            integral_comp.set(j, integral.get(i));
345        }
346        integral_comp
347    }
348
349    /// Return the number of moles of each component in the system.
350    pub fn moles(&self) -> Moles<DVector<f64>> {
351        self.integrate_segments(&self.density)
352    }
353
354    /// Return the total number of moles in the system.
355    pub fn total_moles(&self) -> Moles {
356        self.moles().sum()
357    }
358}
359
360impl<D: Dimension, F> DFTProfile<D, F>
361where
362    D::Larger: Dimension<Smaller = D>,
363    <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
364    F: HelmholtzEnergyFunctional,
365{
366    pub fn weighted_densities(&self) -> FeosResult<Vec<Array<f64, D::Larger>>> {
367        Ok(self
368            .convolver
369            .weighted_densities(&self.density.to_reduced()))
370    }
371
372    #[expect(clippy::type_complexity)]
373    pub fn residual(&self, log: bool) -> FeosResult<(Array<f64, D::Larger>, Array1<f64>, f64)> {
374        // Read from profile
375        let density = self.density.to_reduced();
376        let partial_density = self.bulk.partial_density().into_reduced();
377        let bulk_density = self
378            .bulk
379            .eos
380            .component_index()
381            .iter()
382            .map(|&i| partial_density[i])
383            .collect();
384
385        let (res, res_bulk, res_norm, _, _) =
386            self.euler_lagrange_equation(&density, &bulk_density, log)?;
387        Ok((res, res_bulk, res_norm))
388    }
389
390    #[expect(clippy::type_complexity)]
391    pub(crate) fn euler_lagrange_equation(
392        &self,
393        density: &Array<f64, D::Larger>,
394        bulk_density: &Array1<f64>,
395        log: bool,
396    ) -> FeosResult<(
397        Array<f64, D::Larger>,
398        Array1<f64>,
399        f64,
400        Array<f64, D::Larger>,
401        Array<f64, D::Larger>,
402    )> {
403        // calculate reduced temperature
404        let temperature = self.temperature.to_reduced();
405
406        // calculate intrinsic functional derivative
407        let (_, mut dfdrho) =
408            self.bulk
409                .eos
410                .functional_derivative(temperature, density, self.convolver.as_ref())?;
411
412        // calculate total functional derivative
413        dfdrho += &self.external_potential;
414
415        // calculate bulk functional derivative
416        let bulk_convolver = BulkConvolver::new(self.bulk.eos.weight_functions(temperature));
417        let (_, dfdrho_bulk) = self.bulk.eos.functional_derivative(
418            temperature,
419            bulk_density,
420            bulk_convolver.as_ref(),
421        )?;
422        dfdrho
423            .outer_iter_mut()
424            .zip(dfdrho_bulk)
425            .zip(self.bulk.eos.m().iter())
426            .for_each(|((mut df, df_b), &m)| {
427                df -= df_b;
428                df /= m
429            });
430
431        // calculate bond integrals
432        let exp_dfdrho = dfdrho.mapv(|x| (-x).exp());
433        let bonds = self
434            .bulk
435            .eos
436            .bond_integrals(temperature, &exp_dfdrho, self.convolver.as_ref());
437        let mut rho_projected = &exp_dfdrho * bonds;
438
439        // multiply bulk density
440        rho_projected
441            .outer_iter_mut()
442            .zip(bulk_density.iter())
443            .for_each(|(mut x, &rho_b)| {
444                x *= rho_b;
445            });
446
447        // calculate residual
448        let mut res = if log {
449            rho_projected.mapv(f64::ln) - density.mapv(f64::ln)
450        } else {
451            &rho_projected - density
452        };
453
454        // set residual to 0 where external potentials are overwhelming
455        res.iter_mut()
456            .zip(self.external_potential.iter())
457            .filter(|&(_, &p)| p + f64::EPSILON >= MAX_POTENTIAL)
458            .for_each(|(r, _)| *r = 0.0);
459
460        // additional residuals for the calculation of the bulk densities
461        let z = self.integrate_reduced_comp(&rho_projected);
462        let res_bulk = bulk_density
463            - self
464                .specification
465                .calculate_bulk_density(self, bulk_density, &z)?;
466
467        // calculate the norm of the residual
468        let res_norm = ((density - &rho_projected).mapv(|x| x * x).sum()
469            + res_bulk.mapv(|x| x * x).sum())
470        .sqrt()
471            / ((res.len() + res_bulk.len()) as f64).sqrt();
472
473        if res_norm.is_finite() {
474            Ok((res, res_bulk, res_norm, exp_dfdrho, rho_projected))
475        } else {
476            Err(FeosError::IterationFailed("Euler-Lagrange equation".into()))
477        }
478    }
479
480    pub fn solve(&mut self, solver: Option<&DFTSolver>, debug: bool) -> FeosResult<()> {
481        // unwrap solver
482        let solver = solver.cloned().unwrap_or_default();
483
484        // Read from profile
485        let component_index = self.bulk.eos.component_index().into_owned();
486        let mut density = self.density.to_reduced();
487        let partial_density = self.bulk.partial_density().into_reduced();
488        let mut bulk_density = component_index
489            .iter()
490            .map(|&i| partial_density[i])
491            .collect();
492
493        // Call solver(s)
494        self.call_solver(&mut density, &mut bulk_density, &solver, debug)?;
495
496        // Update profile
497        self.density = Density::from_reduced(density);
498        let mut partial_density = self.bulk.partial_density();
499        bulk_density
500            .into_iter()
501            .enumerate()
502            .for_each(|(i, r)| partial_density.set(component_index[i], Density::from_reduced(r)));
503        self.bulk = State::new_density(&self.bulk.eos, self.bulk.temperature, partial_density)?;
504
505        Ok(())
506    }
507}
feos_dft/profile/mod.rs

feos_dft/profile/
mod.rs
