Skip to content

Commit d10f83c

Browse files
prehnerprehner
authored
Add new 3d_dft feature to avoid libc dependency (feos-org#51)
* Add new `3d_dft` feature to avoid `libc` dependency * activate 3d_dft for python wheels * slightly rearrange feature dependencies * undo changes to main Cargo.toml * add 3d_dft feature to docs.rs metadata * update changelog
1 parent 32e0c3c commit d10f83c

File tree

11 files changed

+471
-468
lines changed

11 files changed

+471
-468
lines changed

‎feos-dft/CHANGELOG.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,8 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
55
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
66

77
## [Unreleased]
8+
### Changed
9+
- The 3D DFT functionalities (3D pores, solvation free energy, free-energy-averaged potentials) are hidden behind the new `3d_dft` feature. For the Python package, the feature is always enabled. [#51](https://github.com/feos-org/feos/pull/51)
810

911
## [0.3.1] - 2022-08-26
1012
### Changed

‎feos-dft/Cargo.toml‎

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -14,23 +14,24 @@ exclude = ["/.github/*", "*.ipynb"]
1414

1515
[package.metadata.docs.rs]
1616
rustdoc-args = [ "--html-in-header", "./docs-header.html" ]
17+
features = [ "3d_dft" ]
1718

1819
[dependencies]
1920
quantity = { version = "0.5", features = ["linalg"] }
2021
feos-core = { version = "0.3", path = "../feos-core" }
2122
num-dual = "0.5"
22-
ndarray = { version = "0.15", features = ["serde", "rayon"] }
23-
ndarray-stats = "0.5"
23+
ndarray = "0.15"
2424
rustdct = "0.7"
2525
rustfft = "6.0"
2626
ang = "0.6"
2727
num-traits = "0.2"
28-
libc = "0.2"
29-
gauss-quad = "0.1"
28+
libm = "0.2"
29+
gauss-quad = { version = "0.1", optional = true }
3030
petgraph = "0.6"
3131
numpy = { version = "0.16", optional = true }
3232
pyo3 = { version = "0.16", optional = true }
3333

3434
[features]
3535
default = []
36-
python = ["pyo3", "numpy", "feos-core/python"]
36+
3d_dft = ["gauss-quad", "ndarray/rayon"]
37+
python = ["pyo3", "numpy", "feos-core/python", "3d_dft"]

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

Lines changed: 21 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,14 @@
1+
#[cfg(feature = "3d_dft")]
12
use crate::adsorption::fea_potential::calculate_fea_potential;
23
use crate::functional::HelmholtzEnergyFunctional;
4+
#[cfg(feature = "3d_dft")]
35
use crate::geometry::Geometry;
46
use feos_core::EosUnit;
5-
use libc::c_double;
7+
use libm::tgamma;
68
use ndarray::{Array1, Array2, Axis as Axis_nd};
9+
#[cfg(feature = "3d_dft")]
710
use quantity::{QuantityArray2, QuantityScalar};
8-
use std::f64::consts::PI;
11+
use std::{f64::consts::PI, marker::PhantomData};
912

1013
const DELTA_STEELE: f64 = 3.35;
1114

@@ -49,6 +52,7 @@ pub enum ExternalPotential<U> {
4952
rho_s: f64,
5053
},
5154
/// Free-energy averaged potential:
55+
#[cfg(feature = "3d_dft")]
5256
FreeEnergyAveraged {
5357
coordinates: QuantityArray2<U>,
5458
sigma_ss: Array1<f64>,
@@ -61,6 +65,10 @@ pub enum ExternalPotential<U> {
6165

6266
/// Custom potential
6367
Custom(Array2<f64>),
68+
69+
/// Needed to keep `FreeEnergyAveraged` optional
70+
#[doc(hidden)]
71+
Phantom(PhantomData<U>),
6472
}
6573

6674
/// Parameters of the fluid required to evaluate the external potential.
@@ -69,6 +77,7 @@ pub trait FluidParameters: HelmholtzEnergyFunctional {
6977
fn sigma_ff(&self) -> &Array1<f64>;
7078
}
7179

80+
#[allow(unused_variables)]
7281
impl<U: EosUnit> ExternalPotential<U> {
7382
// Evaluate the external potential in cartesian coordinates for a given grid and fluid parameters.
7483
pub fn calculate_cartesian_potential<P: FluidParameters>(
@@ -183,6 +192,7 @@ impl<U: EosUnit> ExternalPotential<U> {
183192
* (2.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
184193
- 15.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(3))))
185194
}
195+
#[cfg(feature = "3d_dft")]
186196
Self::FreeEnergyAveraged {
187197
coordinates,
188198
sigma_ss,
@@ -211,7 +221,7 @@ impl<U: EosUnit> ExternalPotential<U> {
211221
*cutoff_radius,
212222
)
213223
}
214-
Self::Custom(_) => unreachable!(),
224+
_ => unreachable!(),
215225
});
216226
}
217227
ext_pot
@@ -348,6 +358,7 @@ impl<U: EosUnit> ExternalPotential<U> {
348358
* sigma_sf[i].powi(3)
349359
* *rho_s)
350360
}
361+
#[cfg(feature = "3d_dft")]
351362
Self::FreeEnergyAveraged {
352363
coordinates,
353364
sigma_ss,
@@ -376,7 +387,7 @@ impl<U: EosUnit> ExternalPotential<U> {
376387
*cutoff_radius,
377388
)
378389
}
379-
Self::Custom(_) => unreachable!(),
390+
_ => unreachable!(),
380391
});
381392
}
382393
ext_pot
@@ -537,6 +548,7 @@ impl<U: EosUnit> ExternalPotential<U> {
537548
* (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
538549
- sum_n(4, r_grid, sigma_sf[i], pore_size)))
539550
}
551+
#[cfg(feature = "3d_dft")]
540552
Self::FreeEnergyAveraged {
541553
coordinates,
542554
sigma_ss,
@@ -565,7 +577,7 @@ impl<U: EosUnit> ExternalPotential<U> {
565577
*cutoff_radius,
566578
)
567579
}
568-
Self::Custom(_) => unreachable!(),
580+
_ => unreachable!(),
569581
});
570582
}
571583
ext_pot
@@ -578,17 +590,17 @@ fn phi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
578590

579591
(1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(m3n2)) * 4.0 * PI.sqrt() / n2m3
580592
* sigma_r.powf(n2m3)
581-
* gamma(n as f64 - 0.5)
582-
/ gamma(n as f64)
593+
* tgamma(n as f64 - 0.5)
594+
/ tgamma(n as f64)
583595
* taylor_2f1_phi(r_r, n)
584596
}
585597

586598
fn psi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
587599
(1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(2.0 - 2.0 * n as f64))
588600
* 4.0
589601
* PI.sqrt()
590-
* gamma(n as f64 - 0.5)
591-
/ gamma(n as f64)
602+
* tgamma(n as f64 - 0.5)
603+
/ tgamma(n as f64)
592604
* sigma_r.powf(2.0 * n as f64 - 2.0)
593605
* taylor_2f1_psi(r_r, n)
594606
}
@@ -637,11 +649,3 @@ fn taylor_2f1_psi(x: &Array1<f64>, n: i32) -> Array1<f64> {
637649
_ => unreachable!(),
638650
}
639651
}
640-
641-
extern "C" {
642-
fn tgamma(x: c_double) -> c_double;
643-
}
644-
645-
fn gamma(x: f64) -> f64 {
646-
unsafe { tgamma(x) }
647-
}

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

Lines changed: 3 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
1+
use super::pore3d::{calculate_distance2, evaluate_lj_potential};
12
use crate::profile::{CUTOFF_RADIUS, MAX_POTENTIAL};
23
use crate::Geometry;
34
use feos_core::EosUnit;
45
use gauss_quad::GaussLegendre;
56
use ndarray::{Array1, Array2, Zip};
6-
use ndarray_stats::SummaryStatisticsExt;
77
use quantity::{QuantityArray2, QuantityScalar};
88
use std::f64::consts::PI;
99
use std::usize;
@@ -114,7 +114,7 @@ pub fn calculate_fea_potential<U: EosUnit>(
114114
let distance2 = calculate_distance2(point, &coordinates, system_size);
115115
let potential_sum: f64 = (0..sigma_sf.len())
116116
.map(|alpha| {
117-
mi * evaluate(
117+
mi * evaluate_lj_potential(
118118
distance2[alpha],
119119
sigma_sf[alpha],
120120
epsilon_k_sf[alpha],
@@ -125,44 +125,8 @@ pub fn calculate_fea_potential<U: EosUnit>(
125125
potential_2d[[i1, i2]] = (-potential_sum.min(MAX_POTENTIAL)).exp();
126126
}
127127
}
128-
*f = potential_2d.weighted_sum(&weights).unwrap();
128+
*f = (potential_2d * &weights).sum();
129129
});
130130

131131
-temperature * potential.map(|p| (p / weights_sum).ln())
132132
}
133-
134-
// -temperature * potential.map(|p| (p / weights_sum).ln())
135-
136-
/// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment
137-
fn evaluate(distance2: f64, sigma: f64, epsilon: f64, cutoff_radius2: f64) -> f64 {
138-
let sigma_r = sigma.powi(2) / distance2;
139-
140-
let potential: f64 = if distance2 > cutoff_radius2 {
141-
0.0
142-
} else if distance2 == 0.0 {
143-
f64::INFINITY
144-
} else {
145-
4.0 * epsilon * (sigma_r.powi(6) - sigma_r.powi(3))
146-
};
147-
148-
potential
149-
}
150-
151-
/// Evaluate the squared euclidian distance between a point and the coordinates of all solid atoms.
152-
fn calculate_distance2(
153-
point: [f64; 3],
154-
coordinates: &Array2<f64>,
155-
system_size: [f64; 3],
156-
) -> Array1<f64> {
157-
Array1::from_shape_fn(coordinates.ncols(), |i| {
158-
let mut rx = coordinates[[0, i]] - point[0];
159-
let mut ry = coordinates[[1, i]] - point[1];
160-
let mut rz = coordinates[[2, i]] - point[2];
161-
162-
rx = rx - system_size[0] * (rx / system_size[0]).round();
163-
ry = ry - system_size[1] * (ry / system_size[1]).round();
164-
rz = rz - system_size[2] * (rz / system_size[2]).round();
165-
166-
rx.powi(2) + ry.powi(2) + rz.powi(2)
167-
})
168-
}

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

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,16 @@ use std::iter;
1111
use std::rc::Rc;
1212

1313
mod external_potential;
14+
#[cfg(feature = "3d_dft")]
1415
mod fea_potential;
1516
mod pore;
1617
pub use external_potential::{ExternalPotential, FluidParameters};
17-
pub use pore::{Pore1D, Pore3D, PoreProfile, PoreProfile1D, PoreProfile3D, PoreSpecification};
18+
pub use pore::{Pore1D, PoreProfile, PoreProfile1D, PoreSpecification};
19+
20+
#[cfg(feature = "3d_dft")]
21+
mod pore3d;
22+
#[cfg(feature = "3d_dft")]
23+
pub use pore3d::{Pore3D, PoreProfile3D};
1824

1925
const MAX_ITER_ADSORPTION_EQUILIBRIUM: usize = 50;
2026
const TOL_ADSORPTION_EQUILIBRIUM: f64 = 1e-8;

0 commit comments

Comments
 (0)