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

feos_dft/profile/
mod.rs
