Non orthogonal dft (#176)

feos-core/Cargo.toml

@@ -35,6 +35,7 @@ typenum = "1.16"
 approx = "0.4"
 regex = "1.9"
 once_cell = "1.18"
+ang = "0.6"
 
 [features]
 default = []

feos-core/src/si/mod.rs

@@ -1,3 +1,4 @@
+use ang::{Angle, Degrees, Radians};
 use num_traits::Zero;
 use std::marker::PhantomData;
 use std::ops::{Div, Mul};
@@ -140,6 +141,11 @@ pub const HOUR: Time<f64> = Quantity(3600.0, PhantomData);
 pub const LITER: Volume<f64> = Quantity(1e-3, PhantomData);
 pub const MINUTE: Time<f64> = Quantity(60.0, PhantomData);
 
+/// Angle unit radian $\\left(\text{rad}\\right)$
+pub const RADIANS: Angle = Radians(1.0);
+/// Angle unit degree $\\left(1^\\circ=\frac{\pi}{180}\\,\text{rad}\\approx 0.0174532925\\,\text{rad}\\right)$
+pub const DEGREES: Angle = Degrees(1.0);
+
 /// Boltzmann constant $\\left(k_\text{B}=1.380649\times 10^{-23}\\,\\frac{\text{J}}{\text{K}}\\right)$
 pub const KB: Entropy<f64> = Quantity(1.380649e-23, PhantomData);
 /// Avogadro constant $\\left(N_\text{A}=6.02214076\times 10^{23}\\,\text{mol}^{-1}\\right)$

feos-dft/CHANGELOG.md

@@ -6,7 +6,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 ### Added
-- Implemented `HelmholtzEnergyFunctional` for `EquationOfState` to be able to use functionals as equations of state. [#158](https://github.com/feos-org/feos/pull/158) 
+- Implemented `HelmholtzEnergyFunctional` for `EquationOfState` to be able to use functionals as equations of state. [#158](https://github.com/feos-org/feos/pull/158)
+- Added the possibility to specify the angles of non-orthorombic unit cells. Currently, the external potential must be specified if non-orthorombic unit cells are calculated. [#176](https://github.com/feos-org/feos/pull/176)
 
 ### Changed
 - `HelmholtzEnergyFunctional`: added `Components` trait as trait bound and removed `ideal_gas` method. [#158](https://github.com/feos-org/feos/pull/158)
@@ -18,6 +19,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Removed `DefaultIdealGasContribution` [#158](https://github.com/feos-org/feos/pull/158)
 - Removed getters for `chemical_potential` (for profiles) and `molar_gibbs_energy` (for `Adsorption1D` and `Adsorption3D`) from Python interface. [#158](https://github.com/feos-org/feos/pull/158)
 
+### Fixed
+- Added an error, if the unit cells in 3D DFT are too small and violate the minimum image convention. [#176](https://github.com/feos-org/feos/pull/176)
+
 ### Packaging
 - Updated `num-dual` dependency to 0.7. [#137](https://github.com/feos-org/feos/pull/137)
 

feos-dft/src/adsorption/mod.rs

@@ -14,8 +14,10 @@ mod external_potential;
 #[cfg(feature = "rayon")]
 mod fea_potential;
 mod pore;
+mod pore2d;
 pub use external_potential::{ExternalPotential, FluidParameters};
 pub use pore::{Pore1D, PoreProfile, PoreProfile1D, PoreSpecification};
+pub use pore2d::{Pore2D, PoreProfile2D};
 
 #[cfg(feature = "rayon")]
 mod pore3d;
@@ -212,7 +214,7 @@ where
                 _ => unreachable!(),
             };
         }
-        let profile = pore.initialize(&bulk, None, None).solve(solver)?.profile;
+        let profile = pore.initialize(&bulk, None, None)?.solve(solver)?.profile;
         let external_potential = Some(&profile.external_potential);
         let mut old_density = Some(&profile.density);
 
@@ -229,8 +231,8 @@ where
                     .clone();
             }
 
-            let p = pore.initialize(&bulk, old_density, external_potential);
-            let p2 = pore.initialize(&bulk, None, external_potential);
+            let p = pore.initialize(&bulk, old_density, external_potential)?;
+            let p2 = pore.initialize(&bulk, None, external_potential)?;
             profiles.push(p.solve(solver).or_else(|_| p2.solve(solver)));
 
             old_density = if let Some(Ok(l)) = profiles.last() {
@@ -277,8 +279,8 @@ where
             .vapor()
             .build()?;
 
-        let mut vapor = pore.initialize(&vapor_bulk, None, None).solve(solver)?;
-        let mut liquid = pore.initialize(&bulk_init, None, None).solve(solver)?;
+        let mut vapor = pore.initialize(&vapor_bulk, None, None)?.solve(solver)?;
+        let mut liquid = pore.initialize(&bulk_init, None, None)?.solve(solver)?;
 
         // calculate initial value for bulk density
         let n_dp_drho_v = (vapor.profile.moles() * vapor_bulk.dp_drho(Contributions::Total)).sum();

feos-dft/src/adsorption/pore.rs

@@ -53,7 +53,7 @@ pub trait PoreSpecification<D: Dimension> {
         bulk: &State<DFT<F>>,
         density: Option<&Density<Array<f64, D::Larger>>>,
         external_potential: Option<&Array<f64, D::Larger>>,
-    ) -> PoreProfile<D, F>;
+    ) -> EosResult<PoreProfile<D, F>>;
 
     /// Return the pore volume using Helium at 298 K as reference.
     fn pore_volume(&self) -> EosResult<Volume<f64>>
@@ -64,7 +64,7 @@ pub trait PoreSpecification<D: Dimension> {
             .temperature(298.0 * KELVIN)
             .density(Density::from_reduced(1.0))
             .build()?;
-        let pore = self.initialize(&bulk, None, None);
+        let pore = self.initialize(&bulk, None, None)?;
         let pot = Dimensionless::from_reduced(
             pore.profile
                 .external_potential
@@ -151,7 +151,7 @@ impl PoreSpecification<Ix1> for Pore1D {
         bulk: &State<DFT<F>>,
         density: Option<&Density<Array2<f64>>>,
         external_potential: Option<&Array2<f64>>,
-    ) -> PoreProfile1D<F> {
+    ) -> EosResult<PoreProfile1D<F>> {
         let dft: &F = &bulk.eos;
         let n_grid = self.n_grid.unwrap_or(DEFAULT_GRID_POINTS);
 
@@ -191,11 +191,11 @@ impl PoreSpecification<Ix1> for Pore1D {
         let weight_functions = dft.weight_functions(t);
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
 
-        PoreProfile {
+        Ok(PoreProfile {
             profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), density),
             grand_potential: None,
             interfacial_tension: None,
-        }
+        })
     }
 }
 

feos-dft/src/adsorption/pore2d.rs

@@ -0,0 +1,53 @@
+use super::{FluidParameters, PoreProfile, PoreSpecification};
+use crate::{Axis, ConvolverFFT, DFTProfile, Grid, HelmholtzEnergyFunctional, DFT};
+use ang::Angle;
+use feos_core::si::{Density, Length};
+use feos_core::{EosResult, State};
+use ndarray::{Array3, Ix2};
+
+pub struct Pore2D {
+    system_size: [Length<f64>; 2],
+    angle: Angle,
+    n_grid: [usize; 2],
+}
+
+pub type PoreProfile2D<F> = PoreProfile<Ix2, F>;
+
+impl Pore2D {
+    pub fn new(system_size: [Length<f64>; 2], angle: Angle, n_grid: [usize; 2]) -> Self {
+        Self {
+            system_size,
+            angle,
+            n_grid,
+        }
+    }
+}
+
+impl PoreSpecification<Ix2> for Pore2D {
+    fn initialize<F: HelmholtzEnergyFunctional + FluidParameters>(
+        &self,
+        bulk: &State<DFT<F>>,
+        density: Option<&Density<Array3<f64>>>,
+        external_potential: Option<&Array3<f64>>,
+    ) -> EosResult<PoreProfile<Ix2, F>> {
+        let dft: &F = &bulk.eos;
+
+        // generate grid
+        let x = Axis::new_cartesian(self.n_grid[0], self.system_size[0], None);
+        let y = Axis::new_cartesian(self.n_grid[1], self.system_size[1], None);
+
+        // temperature
+        let t = bulk.temperature.to_reduced();
+
+        // initialize convolver
+        let grid = Grid::Periodical2(x, y, self.angle);
+        let weight_functions = dft.weight_functions(t);
+        let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
+
+        Ok(PoreProfile {
+            profile: DFTProfile::new(grid, convolver, bulk, external_potential.cloned(), density),
+            grand_potential: None,
+            interfacial_tension: None,
+        })
+    }
+}

feos-dft/src/adsorption/pore3d.rs

@@ -4,14 +4,16 @@ use crate::convolver::ConvolverFFT;
 use crate::functional::{HelmholtzEnergyFunctional, DFT};
 use crate::geometry::{Axis, Grid};
 use crate::profile::{DFTProfile, CUTOFF_RADIUS, MAX_POTENTIAL};
-use feos_core::si::{Density, Length};
-use feos_core::State;
+use ang::Angle;
+use feos_core::si::{Density, Length, DEGREES};
+use feos_core::{EosError, EosResult, State};
 use ndarray::prelude::*;
 use ndarray::Zip;
 
 /// Parameters required to specify a 3D pore.
 pub struct Pore3D {
     system_size: [Length<f64>; 3],
+    angles: Option<[Angle; 3]>,
     n_grid: [usize; 3],
     coordinates: Length<Array2<f64>>,
     sigma_ss: Array1<f64>,
@@ -23,6 +25,7 @@ pub struct Pore3D {
 impl Pore3D {
     pub fn new(
         system_size: [Length<f64>; 3],
+        angles: Option<[Angle; 3]>,
         n_grid: [usize; 3],
         coordinates: Length<Array2<f64>>,
         sigma_ss: Array1<f64>,
@@ -32,6 +35,7 @@ impl Pore3D {
     ) -> Self {
         Self {
             system_size,
+            angles,
             n_grid,
             coordinates,
             sigma_ss,
@@ -51,22 +55,27 @@ impl PoreSpecification<Ix3> for Pore3D {
         bulk: &State<DFT<F>>,
         density: Option<&Density<Array4<f64>>>,
         external_potential: Option<&Array4<f64>>,
-    ) -> PoreProfile3D<F> {
+    ) -> EosResult<PoreProfile3D<F>> {
         let dft: &F = &bulk.eos;
 
         // generate grid
         let x = Axis::new_cartesian(self.n_grid[0], self.system_size[0], None);
         let y = Axis::new_cartesian(self.n_grid[1], self.system_size[1], None);
         let z = Axis::new_cartesian(self.n_grid[2], self.system_size[2], None);
 
-        // move center of geometry of solute to box center
-        let coordinates = Array2::from_shape_fn(self.coordinates.raw_dim(), |(i, j)| {
-            (self.coordinates.get((i, j))).to_reduced()
-        });
+        let coordinates = self.coordinates.to_reduced();
 
         // temperature
         let t = bulk.temperature.to_reduced();
 
+        // For non-orthorombic unit cells, the external potential has to be
+        // provided at the moment
+        if let (Some(_), None) = (self.angles, external_potential) {
+            return Err(EosError::UndeterminedState(
+                "For non-orthorombic unit cells, the external potential has to provided!".into(),
+            ));
+        }
+
         // calculate external potential
         let external_potential = external_potential.map_or_else(
             || {
@@ -82,19 +91,19 @@ impl PoreSpecification<Ix3> for Pore3D {
                     t,
                 )
             },
-            |e| e.clone(),
-        );
+            |e| Ok(e.clone()),
+        )?;
 
         // initialize convolver
-        let grid = Grid::Periodical3(x, y, z);
+        let grid = Grid::Periodical3(x, y, z, self.angles.unwrap_or([90.0 * DEGREES; 3]));
         let weight_functions = dft.weight_functions(t);
         let convolver = ConvolverFFT::plan(&grid, &weight_functions, Some(1));
 
-        PoreProfile {
+        Ok(PoreProfile {
             profile: DFTProfile::new(grid, convolver, bulk, Some(external_potential), density),
             grand_potential: None,
             interfacial_tension: None,
-        }
+        })
     }
 }
 
@@ -108,7 +117,7 @@ pub fn external_potential_3d<F: FluidParameters>(
     cutoff_radius: Option<Length<f64>>,
     potential_cutoff: Option<f64>,
     reduced_temperature: f64,
-) -> Array4<f64> {
+) -> EosResult<Array4<f64>> {
     // allocate external potential
     let m = functional.m();
     let mut external_potential = Array4::zeros((
@@ -128,6 +137,12 @@ pub fn external_potential_3d<F: FluidParameters>(
         .unwrap_or(Length::from_reduced(CUTOFF_RADIUS))
         .to_reduced();
 
+    if system_size.iter().any(|&s| s < 2.0 * cutoff_radius) {
+        return Err(EosError::UndeterminedState(
+            "The unit cell is smaller than 2*cutoff".into(),
+        ));
+    }
+
     // square cut-off radius
     let cutoff_radius2 = cutoff_radius.powi(2);
 
@@ -163,7 +178,7 @@ pub fn external_potential_3d<F: FluidParameters>(
         }
     });
 
-    external_potential
+    Ok(external_potential)
 }
 
 /// Evaluate LJ12-6 potential between solid site "alpha" and fluid segment

feos-dft/src/convolver/mod.rs

@@ -146,10 +146,12 @@ where
                 CurvilinearConvolver::new(r, &[z], weight_functions, lanczos)
             }
             Grid::Cartesian2(x, y) => Self::new(Some(x), &[y], weight_functions, lanczos),
-            Grid::Periodical2(x, y) => PeriodicConvolver::new(&[x, y], weight_functions, lanczos),
+            Grid::Periodical2(x, y, alpha) => {
+                PeriodicConvolver::new_2d(&[x, y], *alpha, weight_functions, lanczos)
+            }
             Grid::Cartesian3(x, y, z) => Self::new(Some(x), &[y, z], weight_functions, lanczos),
-            Grid::Periodical3(x, y, z) => {
-                PeriodicConvolver::new(&[x, y, z], weight_functions, lanczos)
+            Grid::Periodical3(x, y, z, angles) => {
+                PeriodicConvolver::new_3d(&[x, y, z], *angles, weight_functions, lanczos)
             }
         }
     }

feos-dft/src/convolver/periodic_convolver.rs

@@ -1,6 +1,7 @@
 use super::{Convolver, FFTWeightFunctions};
 use crate::geometry::Axis;
 use crate::weight_functions::{WeightFunction, WeightFunctionInfo};
+use ang::Angle;
 use ndarray::Axis as Axis_nd;
 use ndarray::*;
 use num_dual::DualNum;
@@ -30,8 +31,45 @@ where
     D::Larger: Dimension<Smaller = D>,
     <D::Larger as Dimension>::Larger: Dimension<Smaller = D::Larger>,
 {
-    pub fn new(
+    pub fn new_2d(
         axes: &[&Axis],
+        angle: Angle,
+        weight_functions: &[WeightFunctionInfo<T>],
+        lanczos: Option<i32>,
+    ) -> Arc<dyn Convolver<T, D>> {
+        let f = |k: &mut Array<f64, D::Larger>| {
+            let k_y =
+                (&k.index_axis(Axis(0), 1) - &k.index_axis(Axis(0), 0) * angle.cos()) / angle.sin();
+            k.index_axis_mut(Axis(0), 1).assign(&k_y);
+        };
+        Self::new(axes, f, weight_functions, lanczos)
+    }
+
+    pub fn new_3d(
+        axes: &[&Axis],
+        angles: [Angle; 3],
+        weight_functions: &[WeightFunctionInfo<T>],
+        lanczos: Option<i32>,
+    ) -> Arc<dyn Convolver<T, D>> {
+        let f = |k: &mut Array<f64, D::Larger>| {
+            let [alpha, beta, gamma] = angles;
+            let [k_u, k_v, k_w] = [0, 1, 2].map(|i| k.index_axis(Axis(0), i));
+            let k_y = (&k_v - &k_u * gamma.cos()) / gamma.sin();
+            let xi = (alpha.cos() - gamma.cos() * beta.cos()) / gamma.sin();
+            let zeta = (1.0 - beta.cos().powi(2) - xi * xi).sqrt();
+            let k_z = ((gamma.cos() * xi / gamma.sin() - beta.cos()) * &k_u
+                - xi / gamma.sin() * &k_v
+                + &k_w)
+                / zeta;
+            k.index_axis_mut(Axis(0), 1).assign(&k_y);
+            k.index_axis_mut(Axis(0), 2).assign(&k_z);
+        };
+        Self::new(axes, f, weight_functions, lanczos)
+    }
+
+    pub fn new<F: Fn(&mut Array<f64, D::Larger>)>(
+        axes: &[&Axis],
+        non_orthogonal_correction: F,
         weight_functions: &[WeightFunctionInfo<T>],
         lanczos: Option<i32>,
     ) -> Arc<dyn Convolver<T, D>> {
@@ -59,11 +97,18 @@ where
         let mut dim = vec![k_vec.len()];
         k_vec.iter().for_each(|k_x| dim.push(k_x.len()));
         let mut k: Array<_, D::Larger> = Array::zeros(dim).into_dimensionality().unwrap();
-        let mut k_abs = Array::zeros(k.raw_dim().remove_axis(Axis_nd(0)));
         for (i, (mut k_i, k_x)) in k.outer_iter_mut().zip(k_vec.iter()).enumerate() {
             k_i.lanes_mut(Axis_nd(i))
                 .into_iter()
                 .for_each(|mut l| l.assign(k_x));
+        }
+
+        // Correction for non-orthogonal coordinate systems
+        non_orthogonal_correction(&mut k);
+
+        // Calculate the absolute value of the k vector
+        let mut k_abs = Array::zeros(k.raw_dim().remove_axis(Axis_nd(0)));
+        for k_i in k.outer_iter() {
             k_abs.add_assign(&k_i.mapv(|k| k.powi(2)));
         }
         k_abs.map_inplace(|k| *k = k.sqrt());