feos-org · prehner · Jun 10, 2022 · May 24, 2022 · May 24, 2022 · Jun 10, 2022
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 - Made FMT functional more flexible w.r.t. the shape of the weight functions. [#31](https://github.com/feos-org/feos-dft/pull/31)
+- Changed interface of `PairCorrelationFunction` to facilitate the calculation of pair correlation functions in mixtures. [#29](https://github.com/feos-org/feos-dft/pull/29)
 
 ## [0.2.0] - 2022-04-12
 ### Added

diff --git a/build_wheel/src/lib.rs b/build_wheel/src/lib.rs
@@ -52,6 +52,7 @@ impl_pair_correlation!(FMTFunctional);
 #[pymodule]
 pub fn feos_dft(py: Python<'_>, m: &PyModule) -> PyResult<()> {
     m.add_class::<Verbosity>()?;
+    m.add_class::<Contributions>()?;
 
     m.add_class::<PyExternalPotential>()?;
     m.add_class::<Geometry>()?;

diff --git a/examples/FundamentalMeasureTheory.ipynb b/examples/FundamentalMeasureTheory.ipynb
diff --git a/src/fundamental_measure_theory.rs b/src/fundamental_measure_theory.rs
@@ -352,9 +352,10 @@ impl HelmholtzEnergyFunctional for FMTFunctional {
 }
 
 impl PairPotential for FMTFunctional {
-    fn pair_potential(&self, r: &Array1<f64>) -> Array2<f64> {
-        Array::from_shape_fn((self.properties.sigma.len(), r.len()), |(i, j)| {
-            if r[j] > self.properties.sigma[i] {
+    fn pair_potential(&self, i: usize, r: &Array1<f64>, _: f64) -> Array2<f64> {
+        let s = &self.properties.sigma;
+        Array::from_shape_fn((s.len(), r.len()), |(j, k)| {
+            if r[k] > 0.5 * (s[i] + s[j]) {
                 0.0
             } else {
                 f64::INFINITY

diff --git a/src/python/solvation.rs b/src/python/solvation.rs
@@ -81,6 +81,8 @@ macro_rules! impl_pair_correlation {
         /// ----------
         /// bulk : State
         ///     The bulk state in equilibrium with the profile.
+        /// test_particle : int
+        ///     The index of the test particle.
         /// n_grid : int
         ///     The number of grid points.
         /// width: SINumber
@@ -91,16 +93,21 @@ macro_rules! impl_pair_correlation {
         /// PairCorrelation
         ///
         #[pyclass(name = "PairCorrelation", unsendable)]
-        #[pyo3(text_signature = "(bulk, n_grid, width)")]
+        #[pyo3(text_signature = "(bulk, test_particle, n_grid, width)")]
         pub struct PyPairCorrelation(PairCorrelation<SIUnit, $func>);
 
         impl_1d_profile!(PyPairCorrelation, [get_r]);
 
         #[pymethods]
         impl PyPairCorrelation {
             #[new]
-            fn new(bulk: PyState, n_grid: usize, width: PySINumber) -> PyResult<Self> {
-                let profile = PairCorrelation::new(&bulk.0, n_grid, width.into())?;
+            fn new(
+                bulk: PyState,
+                test_particle: usize,
+                n_grid: usize,
+                width: PySINumber,
+            ) -> PyResult<Self> {
+                let profile = PairCorrelation::new(&bulk.0, test_particle, n_grid, width.into())?;
                 Ok(PyPairCorrelation(profile))
             }
 

diff --git a/src/solvation/pair_correlation.rs b/src/solvation/pair_correlation.rs
@@ -11,7 +11,8 @@ use quantity::QuantityScalar;
 /// The underlying pair potential, that the Helmholtz energy functional
 /// models.
 pub trait PairPotential {
-    fn pair_potential(&self, r: &Array1<f64>) -> Array2<f64>;
+    /// Return the pair potential of particle i with all other particles.
+    fn pair_potential(&self, i: usize, r: &Array1<f64>, temperature: f64) -> Array2<f64>;
 }
 
 /// Density profile and properties of a test particle system.
@@ -36,6 +37,7 @@ impl<U: Clone, F> Clone for PairCorrelation<U, F> {
 impl<U: EosUnit, F: HelmholtzEnergyFunctional + PairPotential> PairCorrelation<U, F> {
     pub fn new(
         bulk: &State<U, DFT<F>>,
+        test_particle: usize,
         n_grid: usize,
         width: QuantityScalar<U>,
     ) -> EosResult<Self> {
@@ -46,7 +48,7 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + PairPotential> PairCorrelation<U
 
         // calculate external potential
         let t = bulk.temperature.to_reduced(U::reference_temperature())?;
-        let mut external_potential = dft.pair_potential(&axis.grid) / t;
+        let mut external_potential = dft.pair_potential(test_particle, &axis.grid, t) / t;
         external_potential.map_inplace(|x| {
             if *x > MAX_POTENTIAL {
                 *x = MAX_POTENTIAL
@@ -71,8 +73,16 @@ impl<U: EosUnit, F: HelmholtzEnergyFunctional + PairPotential> PairCorrelation<U
         self.profile.solve(solver, debug)?;
 
         // calculate pair correlation function
-        self.pair_correlation_function =
-            Some(self.profile.density.to_reduced(self.profile.bulk.density)?);
+        self.pair_correlation_function = Some(Array::from_shape_fn(
+            self.profile.density.raw_dim(),
+            |(i, j)| {
+                self.profile
+                    .density
+                    .get((i, j))
+                    .to_reduced(self.profile.bulk.partial_density.get(i))
+                    .unwrap()
+            },
+        ));
 
         // calculate self solvation free energy
         self.self_solvation_free_energy = Some(self.profile.integrate(

diff --git a/src/solver.rs b/src/solver.rs
@@ -311,7 +311,7 @@ impl fmt::Display for SolverParameter {
 impl fmt::Display for DFTSolver {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         for algorithm in &self.parameters {
-            write!(f, "{}\n", algorithm)?;
+            writeln!(f, "{algorithm}")?;
         }
         Ok(())
     }