Skip to content

Commit 6d34c8d

Browse files
prehnerprehner
authored
Always use the same lanczos sigma factor within one DFT profile (#272)
1 parent 1ced199 commit 6d34c8d

File tree

9 files changed

+38
-57
lines changed

9 files changed

+38
-57
lines changed

‎feos-dft/src/adsorption/pore.rs

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -167,7 +167,8 @@ where
167167
let weight_functions: Vec<WeightFunctionInfo<N>> = functional_contributions
168168
.map(|c| c.weight_functions(temperature))
169169
.collect();
170-
let convolver = ConvolverFFT::<_, D>::plan(&self.profile.grid, &weight_functions, None);
170+
let convolver =
171+
ConvolverFFT::<_, D>::plan(&self.profile.grid, &weight_functions, self.profile.lanczos);
171172
let bonds = self
172173
.profile
173174
.dft
@@ -230,14 +231,11 @@ impl PoreSpecification<Ix1> for Pore1D {
230231
|e| e.clone(),
231232
);
232233

233-
// initialize convolver
234+
// initialize grid
234235
let grid = Grid::new_1d(axis);
235-
let t = bulk.temperature.to_reduced();
236-
let weight_functions = dft.weight_functions(t);
237-
let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
238236

239237
Ok(PoreProfile {
240-
profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), density),
238+
profile: DFTProfile::new(grid, bulk, Some(external_potential), density, Some(1)),
241239
grand_potential: None,
242240
interfacial_tension: None,
243241
})

‎feos-dft/src/adsorption/pore2d.rs

Lines changed: 3 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
use super::{FluidParameters, PoreProfile, PoreSpecification};
2-
use crate::{Axis, ConvolverFFT, DFTProfile, Grid, HelmholtzEnergyFunctional, DFT};
3-
use feos_core::{EosResult, ReferenceSystem, State};
2+
use crate::{Axis, DFTProfile, Grid, HelmholtzEnergyFunctional, DFT};
3+
use feos_core::{EosResult, State};
44
use ndarray::{Array3, Ix2};
55
use quantity::{Angle, Density, Length};
66

@@ -29,22 +29,13 @@ impl PoreSpecification<Ix2> for Pore2D {
2929
density: Option<&Density<Array3<f64>>>,
3030
external_potential: Option<&Array3<f64>>,
3131
) -> EosResult<PoreProfile<Ix2, F>> {
32-
let dft: &F = &bulk.eos;
33-
3432
// generate grid
3533
let x = Axis::new_cartesian(self.n_grid[0], self.system_size[0], None);
3634
let y = Axis::new_cartesian(self.n_grid[1], self.system_size[1], None);
37-
38-
// temperature
39-
let t = bulk.temperature.to_reduced();
40-
41-
// initialize convolver
4235
let grid = Grid::Periodical2(x, y, self.angle);
43-
let weight_functions = dft.weight_functions(t);
44-
let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
4536

4637
Ok(PoreProfile {
47-
profile: DFTProfile::new(grid, convolver, bulk, external_potential.cloned(), density),
38+
profile: DFTProfile::new(grid, bulk, external_potential.cloned(), density, Some(1)),
4839
grand_potential: None,
4940
interfacial_tension: None,
5041
})

‎feos-dft/src/adsorption/pore3d.rs

Lines changed: 1 addition & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,5 @@
11
use super::pore::{PoreProfile, PoreSpecification};
22
use crate::adsorption::FluidParameters;
3-
use crate::convolver::ConvolverFFT;
43
use crate::functional::{HelmholtzEnergyFunctional, DFT};
54
use crate::geometry::{Axis, Grid};
65
use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
@@ -93,14 +92,10 @@ impl PoreSpecification<Ix3> for Pore3D {
9392
},
9493
|e| Ok(e.clone()),
9594
)?;
96-
97-
// initialize convolver
9895
let grid = Grid::Periodical3(x, y, z, self.angles.unwrap_or([90.0 * DEGREES; 3]));
99-
let weight_functions = dft.weight_functions(t);
100-
let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
10196

10297
Ok(PoreProfile {
103-
profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), density),
98+
profile: DFTProfile::new(grid, bulk, Some(external_potential), density, Some(1)),
10499
grand_potential: None,
105100
interfacial_tension: None,
106101
})

‎feos-dft/src/convolver/mod.rs

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -293,7 +293,7 @@ where
293293
.zip(result.lanes_mut(Axis_nd(0)))
294294
{
295295
self.transform
296-
.forward_transform(f, r, vector_index.map_or(true, |ind| ind != 0));
296+
.forward_transform(f, r, vector_index != Some(0));
297297
}
298298
for (i, transform) in self.cartesian_transforms.iter().enumerate() {
299299
dim[i + 1] = self.k_abs.shape()[i + 1];
@@ -303,7 +303,7 @@ where
303303
.into_iter()
304304
.zip(res.lanes_mut(Axis_nd(i + 1)))
305305
{
306-
transform.forward_transform(f, r, vector_index.map_or(true, |ind| ind != i + 1));
306+
transform.forward_transform(f, r, vector_index.is_none_or(|ind| ind != i + 1));
307307
}
308308
result = res;
309309
}
@@ -339,8 +339,7 @@ where
339339
.into_iter()
340340
.zip(res.lanes_mut(Axis_nd(0)))
341341
{
342-
self.transform
343-
.back_transform(f, r, vector_index.map_or(true, |ind| ind != 0));
342+
self.transform.back_transform(f, r, vector_index != Some(0));
344343
}
345344
for (i, transform) in self.cartesian_transforms.iter().enumerate() {
346345
dim[i + 1] = result.shape()[i + 1];
@@ -350,7 +349,7 @@ where
350349
.into_iter()
351350
.zip(res2.lanes_mut(Axis_nd(i + 1)))
352351
{
353-
transform.back_transform(f, r, vector_index.map_or(true, |ind| ind != i + 1));
352+
transform.back_transform(f, r, vector_index.is_none_or(|ind| ind != i + 1));
354353
}
355354
res = res2;
356355
}

‎feos-dft/src/interface/mod.rs

Lines changed: 1 addition & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
11
//! Density profiles at planar interfaces and interfacial tensions.
2-
use crate::convolver::ConvolverFFT;
32
use crate::functional::{HelmholtzEnergyFunctional, DFT};
43
use crate::geometry::{Axis, Grid};
54
use crate::profile::{DFTProfile, DFTSpecifications};
@@ -64,18 +63,11 @@ impl<F: HelmholtzEnergyFunctional> PlanarInterface<F> {
6463

6564
impl<F: HelmholtzEnergyFunctional> PlanarInterface<F> {
6665
pub fn new(vle: &PhaseEquilibrium<DFT<F>, 2>, n_grid: usize, l_grid: Length) -> Self {
67-
let dft = &vle.vapor().eos;
68-
6966
// generate grid
7067
let grid = Grid::Cartesian1(Axis::new_cartesian(n_grid, l_grid, None));
7168

72-
// initialize convolver
73-
let t = vle.vapor().temperature.to_reduced();
74-
let weight_functions = dft.weight_functions(t);
75-
let convolver = ConvolverFFT::plan(&grid, &weight_functions, None);
76-
7769
Self {
78-
profile: DFTProfile::new(grid, convolver, vle.vapor(), None, None),
70+
profile: DFTProfile::new(grid, vle.vapor(), None, None, None),
7971
vle: vle.clone(),
8072
surface_tension: None,
8173
equimolar_radius: None,

‎feos-dft/src/profile/mod.rs

Lines changed: 19 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,11 @@
1-
use crate::convolver::{BulkConvolver, Convolver};
1+
use crate::convolver::{BulkConvolver, Convolver, ConvolverFFT};
22
use crate::functional::{HelmholtzEnergyFunctional, DFT};
33
use crate::geometry::Grid;
44
use crate::solver::{DFTSolver, DFTSolverLog};
55
use feos_core::{Components, EosError, EosResult, ReferenceSystem, State};
66
use ndarray::{
77
Array, Array1, Array2, Array3, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix1, Ix2, Ix3,
8+
RemoveAxis,
89
};
910
use num_dual::DualNum;
1011
use quantity::{Density, Length, Moles, Quantity, Temperature, Volume, _Volume, DEGREES};
@@ -107,6 +108,7 @@ pub struct DFTProfile<D: Dimension, F> {
107108
pub external_potential: Array<f64, D::Larger>,
108109
pub bulk: State<DFT<F>>,
109110
pub solver_log: Option<DFTSolverLog>,
111+
pub lanczos: Option<i32>,
110112
}
111113

112114
impl<F> DFTProfile<Ix1, F> {
@@ -189,9 +191,11 @@ impl<F> DFTProfile<Ix3, F> {
189191
}
190192
}
191193

192-
impl<D: Dimension, F: HelmholtzEnergyFunctional> DFTProfile<D, F>
194+
impl<D: Dimension + RemoveAxis + 'static, F: HelmholtzEnergyFunctional> DFTProfile<D, F>
193195
where
194196
D::Larger: Dimension<Smaller = D>,
197+
D::Smaller: Dimension<Larger = D>,
198+
<D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
195199
{
196200
/// Create a new density profile.
197201
///
@@ -201,13 +205,18 @@ where
201205
/// after this call if something else is required.
202206
pub fn new(
203207
grid: Grid,
204-
convolver: Arc<dyn Convolver<f64, D>>,
205208
bulk: &State<DFT<F>>,
206209
external_potential: Option<Array<f64, D::Larger>>,
207210
density: Option<&Density<Array<f64, D::Larger>>>,
211+
lanczos: Option<i32>,
208212
) -> Self {
209213
let dft = bulk.eos.clone();
210214

215+
// initialize convolver
216+
let t = bulk.temperature.to_reduced();
217+
let weight_functions = dft.weight_functions(t);
218+
let convolver = ConvolverFFT::plan(&grid, &weight_functions, lanczos);
219+
211220
// initialize external potential
212221
let external_potential = external_potential.unwrap_or_else(|| {
213222
let mut n_grid = vec![dft.component_index().len()];
@@ -221,7 +230,6 @@ where
221230
let density = if let Some(density) = density {
222231
density.to_owned()
223232
} else {
224-
let t = bulk.temperature.to_reduced();
225233
let exp_dfdrho = (-&external_potential).mapv(f64::exp);
226234
let mut bonds = dft.bond_integrals(t, &exp_dfdrho, &convolver);
227235
bonds *= &exp_dfdrho;
@@ -245,9 +253,15 @@ where
245253
external_potential,
246254
bulk: bulk.clone(),
247255
solver_log: None,
256+
lanczos,
248257
}
249258
}
259+
}
250260

261+
impl<D: Dimension, F: HelmholtzEnergyFunctional> DFTProfile<D, F>
262+
where
263+
D::Larger: Dimension<Smaller = D>,
264+
{
251265
fn integrate_reduced<N: DualNum<f64> + Copy>(&self, mut profile: Array<N, D>) -> N {
252266
let (integration_weights, functional_determinant) = self.grid.integration_weights();
253267

@@ -358,6 +372,7 @@ impl<D: Dimension + Clone, F> Clone for DFTProfile<D, F> {
358372
external_potential: self.external_potential.clone(),
359373
bulk: self.bulk.clone(),
360374
solver_log: self.solver_log.clone(),
375+
lanczos: self.lanczos,
361376
}
362377
}
363378
}

‎feos-dft/src/profile/properties.rs

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -113,7 +113,7 @@ where
113113
let weight_functions: Vec<WeightFunctionInfo<Dual64>> = functional_contributions
114114
.map(|c| c.weight_functions(temperature_dual))
115115
.collect();
116-
let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, None);
116+
let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, self.lanczos);
117117

118118
let density = self.density.to_reduced();
119119

@@ -196,7 +196,7 @@ where
196196
let weight_functions: Vec<WeightFunctionInfo<Dual64>> = functional_contributions
197197
.map(|c| c.weight_functions(temperature_dual))
198198
.collect();
199-
let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, None);
199+
let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, self.lanczos);
200200

201201
let density = self.density.to_reduced();
202202

@@ -239,7 +239,7 @@ where
239239
let weight_functions: Vec<WeightFunctionInfo<Dual64>> = functional_contributions
240240
.map(|c| c.weight_functions(temperature_dual))
241241
.collect();
242-
let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, None);
242+
let convolver = ConvolverFFT::plan(&self.grid, &weight_functions, self.lanczos);
243243

244244
let density = self.density.to_reduced();
245245

@@ -358,7 +358,7 @@ where
358358
.map(|c| c.weight_functions(t_dual))
359359
.collect();
360360
let convolver: Arc<dyn Convolver<_, D>> =
361-
ConvolverFFT::plan(&self.grid, &weight_functions, None);
361+
ConvolverFFT::plan(&self.grid, &weight_functions, self.lanczos);
362362
let (_, mut dfdrho) = self
363363
.dft
364364
.functional_derivative(t_dual, &rho_dual, &convolver)?;

‎feos-dft/src/solvation/pair_correlation.rs

Lines changed: 1 addition & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
11
//! Functionalities for the calculation of pair correlation functions.
2-
use crate::convolver::ConvolverFFT;
32
use crate::functional::{HelmholtzEnergyFunctional, DFT};
43
use crate::profile::MAX_POTENTIAL;
54
use crate::solver::DFTSolver;
@@ -49,14 +48,10 @@ impl<F: HelmholtzEnergyFunctional + PairPotential> PairCorrelation<F> {
4948
*x = MAX_POTENTIAL
5049
}
5150
});
52-
53-
// initialize convolver
5451
let grid = Grid::Spherical(axis);
55-
let weight_functions = dft.weight_functions(t);
56-
let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
5752

5853
Self {
59-
profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), None),
54+
profile: DFTProfile::new(grid, bulk, Some(external_potential), None, Some(1)),
6055
pair_correlation_function: None,
6156
self_solvation_free_energy: None,
6257
structure_factor: None,

‎feos-dft/src/solvation/solvation_profile.rs

Lines changed: 1 addition & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,4 @@
11
use crate::adsorption::FluidParameters;
2-
use crate::convolver::ConvolverFFT;
32
use crate::functional::{HelmholtzEnergyFunctional, DFT};
43
use crate::geometry::{Axis, Grid};
54
use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
@@ -103,13 +102,10 @@ impl<F: HelmholtzEnergyFunctional + FluidParameters> SolvationProfile<F> {
103102
t,
104103
)?;
105104

106-
// initialize convolver
107105
let grid = Grid::Cartesian3(x, y, z);
108-
let weight_functions = dft.weight_functions(t);
109-
let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
110106

111107
Ok(Self {
112-
profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), None),
108+
profile: DFTProfile::new(grid, bulk, Some(external_potential), None, Some(1)),
113109
grand_potential: None,
114110
solvation_free_energy: None,
115111
})

0 commit comments

Comments
 (0)