feos_dft/adsorption/
external_potential.rs

1#[cfg(feature = "rayon")]
2use crate::adsorption::fea_potential::calculate_fea_potential;
3use crate::functional::HelmholtzEnergyFunctional;
4#[cfg(feature = "rayon")]
5use crate::geometry::Geometry;
6use libm::tgamma;
7use nalgebra::DVector;
8use ndarray::{Array1, Array2, Axis as Axis_nd};
9#[cfg(feature = "rayon")]
10use quantity::Length;
11use std::f64::consts::PI;
12use std::ops::Deref;
13
14const DELTA_STEELE: f64 = 3.35;
15
16/// A collection of external potentials.
17#[derive(Clone)]
18pub enum ExternalPotential {
19    /// Hard wall potential: $V_i^\mathrm{ext}(z)=\begin{cases}\infty&z\leq\sigma_{si}\\\\0&z>\sigma_{si}\end{cases},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$
20    HardWall { sigma_ss: f64 },
21    /// 9-3 Lennard-Jones potential: $V_i^\mathrm{ext}(z)=\frac{2\pi}{45} m_i\varepsilon_{si}\sigma_{si}^3\rho_s\left(2\left(\frac{\sigma_{si}}{z}\right)^9-15\left(\frac{\sigma_{si}}{z}\right)^3\right),~~~~\varepsilon_{si}=\sqrt{\varepsilon_{ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$
22    LJ93 {
23        sigma_ss: f64,
24        epsilon_k_ss: f64,
25        rho_s: f64,
26    },
27    /// Simple 9-3 Lennard-Jones potential: $V_i^\mathrm{ext}(z)=\varepsilon_{si}\left(\left(\frac{\sigma_{si}}{z}\right)^9-\left(\frac{\sigma_{si}}{z}\right)^3\right),~~~~\varepsilon_{si}=\sqrt{\varepsilon_{ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$
28    SimpleLJ93 { sigma_ss: f64, epsilon_k_ss: f64 },
29    /// Custom 9-3 Lennard-Jones potential: $V_i^\mathrm{ext}(z)=\varepsilon_{si}\left(\left(\frac{\sigma_{si}}{z}\right)^9-\left(\frac{\sigma_{si}}{z}\right)^3\right)$
30    CustomLJ93 {
31        sigma_sf: Array1<f64>,
32        epsilon_k_sf: Array1<f64>,
33    },
34    /// Steele potential: $V_i^\mathrm{ext}(z)=2\pi m_i\xi\varepsilon_{si}\sigma_{si}^2\Delta\rho_s\left(0.4\left(\frac{\sigma_{si}}{z}\right)^{10}-\left(\frac{\sigma_{si}}{z}\right)^4-\frac{\sigma_{si}^4}{3\Delta\left(z+0.61\Delta\right)^3}\right),~~~~\varepsilon_{si}=\sqrt{\varepsilon_{ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right),~~~~\Delta=3.35$
35    Steele {
36        sigma_ss: f64,
37        epsilon_k_ss: f64,
38        rho_s: f64,
39        xi: Option<f64>,
40    },
41    /// Steele potential with custom combining rules: $V_i^\mathrm{ext}(z)=2\pi m_i\xi\varepsilon_{si}\sigma_{si}^2\Delta\rho_s\left(0.4\left(\frac{\sigma_{si}}{z}\right)^{10}-\left(\frac{\sigma_{si}}{z}\right)^4-\frac{\sigma_{si}^4}{3\Delta\left(z+0.61\Delta\right)^3}\right),~~~~\Delta=3.35$
42    CustomSteele {
43        sigma_sf: Array1<f64>,
44        epsilon_k_sf: Array1<f64>,
45        rho_s: f64,
46        xi: Option<f64>,
47    },
48    /// Double well potential: $V_i^\mathrm{ext}(z)=\mathrm{min}\left(\frac{2\pi}{45} m_i\varepsilon_{2si}\sigma_{si}^3\rho_s\left(2\left(\frac{2\sigma_{si}}{z}\right)^9-15\left(\frac{2\sigma_{si}}{z}\right)^3\right),0\right)+\frac{2\pi}{45} m_i\varepsilon_{1si}\sigma_{si}^3\rho_s\left(2\left(\frac{\sigma_{si}}{z}\right)^9-15\left(\frac{\sigma_{si}}{z}\right)^3\right),~~~~\varepsilon_{1si}=\sqrt{\varepsilon_{1ss}\varepsilon_{ii}},~~~~\varepsilon_{2si}=\sqrt{\varepsilon_{2ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$
49    DoubleWell {
50        sigma_ss: f64,
51        epsilon1_k_ss: f64,
52        epsilon2_k_ss: f64,
53        rho_s: f64,
54    },
55    /// Free-energy averaged potential:
56    #[cfg(feature = "rayon")]
57    FreeEnergyAveraged {
58        coordinates: Length<Array2<f64>>,
59        sigma_ss: Array1<f64>,
60        epsilon_k_ss: Array1<f64>,
61        pore_center: [f64; 3],
62        system_size: [Length; 3],
63        n_grid: [usize; 2],
64        cutoff_radius: Option<f64>,
65    },
66
67    /// Custom potential
68    Custom(Array2<f64>),
69}
70
71/// Parameters of the fluid required to evaluate the external potential.
72pub trait FluidParameters {
73    fn epsilon_k_ff(&self) -> DVector<f64>;
74    fn sigma_ff(&self) -> DVector<f64>;
75}
76
77impl<C: Deref<Target = T>, T: FluidParameters> FluidParameters for C {
78    fn epsilon_k_ff(&self) -> DVector<f64> {
79        T::epsilon_k_ff(self)
80    }
81    fn sigma_ff(&self) -> DVector<f64> {
82        T::sigma_ff(self)
83    }
84}
85
86impl ExternalPotential {
87    // Evaluate the external potential in cartesian coordinates for a given grid and fluid parameters.
88    pub fn calculate_cartesian_potential<P: HelmholtzEnergyFunctional + FluidParameters>(
89        &self,
90        z_grid: &Array1<f64>,
91        fluid_parameters: &P,
92        #[cfg_attr(not(feature = "rayon"), expect(unused_variables))] temperature: f64,
93    ) -> Array2<f64> {
94        if let ExternalPotential::Custom(potential) = self {
95            return potential.clone();
96        }
97
98        // Allocate external potential
99        let m = fluid_parameters.m();
100        let mut ext_pot = Array2::zeros((m.len(), z_grid.len()));
101
102        for (i, &mi) in m.iter().enumerate() {
103            ext_pot.index_axis_mut(Axis_nd(0), i).assign(&match self {
104                Self::HardWall { sigma_ss } => {
105                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + *sigma_ss) * 0.5;
106                    z_grid.mapv(|z| if z < sigma_sf { f64::INFINITY } else { 0.0 })
107                }
108                Self::LJ93 {
109                    sigma_ss,
110                    epsilon_k_ss,
111                    rho_s,
112                } => {
113                    // combining rules
114                    let epsilon_k_sf =
115                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
116                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
117
118                    2.0 * PI * mi * epsilon_k_sf[i] * sigma_sf[i].powi(3) * rho_s / 45.0
119                        * (2.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
120                            - 15.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(3)))
121                }
122                Self::SimpleLJ93 {
123                    sigma_ss,
124                    epsilon_k_ss,
125                } => {
126                    // combining rules
127                    let epsilon_k_sf =
128                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
129                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
130
131                    epsilon_k_sf[i]
132                        * ((sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
133                            - (sigma_sf[i] / z_grid).mapv(|x| x.powi(3)))
134                }
135                Self::CustomLJ93 {
136                    sigma_sf,
137                    epsilon_k_sf,
138                } => {
139                    epsilon_k_sf[i]
140                        * ((sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
141                            - (sigma_sf[i] / z_grid).mapv(|x| x.powi(3)))
142                }
143                Self::Steele {
144                    sigma_ss,
145                    epsilon_k_ss,
146                    rho_s,
147                    xi,
148                } => {
149                    // combining rules
150                    let epsilon_k_sf =
151                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
152                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
153
154                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
155                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
156                        * (0.4 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(10))
157                            - (sigma_sf[i] / z_grid).mapv(|x| x.powi(4))
158                            - sigma_sf[i].powi(4)
159                                / ((3.0 * DELTA_STEELE)
160                                    * (z_grid + 0.61 * DELTA_STEELE).mapv(|x| x.powi(3))))
161                }
162                Self::CustomSteele {
163                    sigma_sf,
164                    epsilon_k_sf,
165                    rho_s,
166                    xi,
167                } => {
168                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
169                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
170                        * (0.4 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(10))
171                            - (sigma_sf[i] / z_grid).mapv(|x| x.powi(4))
172                            - sigma_sf[i].powi(4)
173                                / ((3.0 * DELTA_STEELE)
174                                    * (z_grid + 0.61 * DELTA_STEELE).mapv(|x| x.powi(3))))
175                }
176                Self::DoubleWell {
177                    sigma_ss,
178                    epsilon1_k_ss,
179                    epsilon2_k_ss,
180                    rho_s,
181                } => {
182                    // combining rules
183                    let epsilon1_k_sf =
184                        (fluid_parameters.epsilon_k_ff() * *epsilon1_k_ss).map(|e| e.sqrt());
185                    let epsilon2_k_sf =
186                        (fluid_parameters.epsilon_k_ff() * *epsilon2_k_ss).map(|e| e.sqrt());
187                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
188
189                    let bh_tail = (2.0 * PI * mi * epsilon2_k_sf[i] * sigma_sf[i].powi(3) * rho_s
190                        / 45.0
191                        * (2.0 * (2.0 * sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
192                            - 15.0 * (2.0 * sigma_sf[i] / z_grid).mapv(|x| x.powi(3))))
193                    .mapv(|x| x.min(0.0));
194
195                    bh_tail
196                        + &(2.0 * PI * mi * epsilon1_k_sf[i] * sigma_sf[i].powi(3) * rho_s / 45.0
197                            * (2.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
198                                - 15.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(3))))
199                }
200                #[cfg(feature = "rayon")]
201                Self::FreeEnergyAveraged {
202                    coordinates,
203                    sigma_ss,
204                    epsilon_k_ss,
205                    pore_center,
206                    system_size,
207                    n_grid,
208                    cutoff_radius,
209                } => {
210                    // combining rules
211                    let epsilon_k_sf =
212                        (fluid_parameters.epsilon_k_ff()[i] * epsilon_k_ss).map(|e| e.sqrt());
213                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + sigma_ss) * 0.5;
214
215                    calculate_fea_potential(
216                        z_grid,
217                        mi,
218                        coordinates,
219                        sigma_sf,
220                        epsilon_k_sf,
221                        pore_center,
222                        system_size,
223                        n_grid,
224                        temperature,
225                        Geometry::Cartesian,
226                        *cutoff_radius,
227                    )
228                }
229                _ => unreachable!(),
230            });
231        }
232        ext_pot
233    }
234
235    // Evaluate the external potential in cylindrical coordinates for a given grid and fluid parameters.
236    pub fn calculate_cylindrical_potential<P: HelmholtzEnergyFunctional + FluidParameters>(
237        &self,
238        r_grid: &Array1<f64>,
239        pore_size: f64,
240        fluid_parameters: &P,
241        #[cfg_attr(not(feature = "rayon"), expect(unused_variables))] temperature: f64,
242    ) -> Array2<f64> {
243        if let ExternalPotential::Custom(potential) = self {
244            return potential.clone();
245        }
246
247        // Allocate external potential
248        let m = fluid_parameters.m();
249        let mut ext_pot = Array2::zeros((m.len(), r_grid.len()));
250
251        for (i, &mi) in m.iter().enumerate() {
252            ext_pot.index_axis_mut(Axis_nd(0), i).assign(&match self {
253                Self::HardWall { sigma_ss } => {
254                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + *sigma_ss) * 0.5;
255                    r_grid.mapv(|r| {
256                        if r > pore_size - sigma_sf {
257                            f64::INFINITY
258                        } else {
259                            0.0
260                        }
261                    })
262                }
263                Self::LJ93 {
264                    sigma_ss,
265                    epsilon_k_ss,
266                    rho_s,
267                } => {
268                    // combining rules
269                    let epsilon_k_sf =
270                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
271                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
272
273                    (phi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size)
274                        - phi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size))
275                        * 2.0
276                        * PI
277                        * mi
278                        * epsilon_k_sf[i]
279                        * sigma_sf[i].powi(3)
280                        * *rho_s
281                }
282                Self::SimpleLJ93 {
283                    sigma_ss: _,
284                    epsilon_k_ss: _,
285                } => {
286                    unimplemented!()
287                }
288                Self::CustomLJ93 {
289                    sigma_sf: _,
290                    epsilon_k_sf: _,
291                } => {
292                    unimplemented!()
293                }
294                Self::Steele {
295                    sigma_ss,
296                    epsilon_k_ss,
297                    rho_s,
298                    xi,
299                } => {
300                    // combining rules
301                    let epsilon_k_sf =
302                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
303                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
304
305                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
306                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
307                        * (psi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size)
308                            - psi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size)
309                            - sigma_sf[i] / DELTA_STEELE
310                                * phi(
311                                    3,
312                                    &(r_grid / (pore_size + DELTA_STEELE * 0.61)),
313                                    sigma_sf[i] / (pore_size + DELTA_STEELE * 0.61),
314                                ))
315                }
316                Self::CustomSteele {
317                    sigma_sf,
318                    epsilon_k_sf,
319                    rho_s,
320                    xi,
321                } => {
322                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
323                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
324                        * (psi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size)
325                            - psi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size)
326                            - sigma_sf[i] / DELTA_STEELE
327                                * phi(
328                                    3,
329                                    &(r_grid / (pore_size + DELTA_STEELE * 0.61)),
330                                    sigma_sf[i] / (pore_size + DELTA_STEELE * 0.61),
331                                ))
332                }
333                Self::DoubleWell {
334                    sigma_ss,
335                    epsilon1_k_ss,
336                    epsilon2_k_ss,
337                    rho_s,
338                } => {
339                    // combining rules
340                    let epsilon1_k_sf =
341                        (fluid_parameters.epsilon_k_ff() * *epsilon1_k_ss).map(|e| e.sqrt());
342                    let epsilon2_k_sf =
343                        (fluid_parameters.epsilon_k_ff() * *epsilon2_k_ss).map(|e| e.sqrt());
344                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
345
346                    let bh_tail = ((phi(6, &(r_grid / pore_size), 2.0 * sigma_sf[i] / pore_size)
347                        - phi(3, &(r_grid / pore_size), 2.0 * sigma_sf[i] / pore_size))
348                        * 2.0
349                        * PI
350                        * mi
351                        * epsilon2_k_sf[i]
352                        * sigma_sf[i].powi(3)
353                        * *rho_s)
354                        .mapv(|x| x.min(0.0));
355
356                    bh_tail
357                        + &((phi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size)
358                            - phi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size))
359                            * 2.0
360                            * PI
361                            * mi
362                            * epsilon1_k_sf[i]
363                            * sigma_sf[i].powi(3)
364                            * *rho_s)
365                }
366                #[cfg(feature = "rayon")]
367                Self::FreeEnergyAveraged {
368                    coordinates,
369                    sigma_ss,
370                    epsilon_k_ss,
371                    pore_center,
372                    system_size,
373                    n_grid,
374                    cutoff_radius,
375                } => {
376                    // combining rules
377                    let epsilon_k_sf =
378                        (fluid_parameters.epsilon_k_ff()[i] * epsilon_k_ss).map(|e| e.sqrt());
379                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + sigma_ss) * 0.5;
380
381                    calculate_fea_potential(
382                        r_grid,
383                        mi,
384                        coordinates,
385                        sigma_sf,
386                        epsilon_k_sf,
387                        pore_center,
388                        system_size,
389                        n_grid,
390                        temperature,
391                        Geometry::Cylindrical,
392                        *cutoff_radius,
393                    )
394                }
395                _ => unreachable!(),
396            });
397        }
398        ext_pot
399    }
400
401    // Evaluate the external potential in spherical coordinates for a given grid and fluid parameters.
402    pub fn calculate_spherical_potential<P: HelmholtzEnergyFunctional + FluidParameters>(
403        &self,
404        r_grid: &Array1<f64>,
405        pore_size: f64,
406        fluid_parameters: &P,
407        #[cfg_attr(not(feature = "rayon"), expect(unused_variables))] temperature: f64,
408    ) -> Array2<f64> {
409        if let ExternalPotential::Custom(potential) = self {
410            return potential.clone();
411        }
412
413        // Allocate external potential
414        let m = fluid_parameters.m();
415        let mut ext_pot = Array2::zeros((m.len(), r_grid.len()));
416
417        for (i, &mi) in m.iter().enumerate() {
418            ext_pot.index_axis_mut(Axis_nd(0), i).assign(&match self {
419                Self::HardWall { sigma_ss } => {
420                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + *sigma_ss) * 0.5;
421                    r_grid.mapv(|r| {
422                        if r > pore_size - sigma_sf {
423                            f64::INFINITY
424                        } else {
425                            0.0
426                        }
427                    })
428                }
429                Self::LJ93 {
430                    sigma_ss,
431                    epsilon_k_ss,
432                    rho_s,
433                } => {
434                    // combining rules
435                    let epsilon_k_sf =
436                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
437                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
438
439                    PI * mi
440                        * epsilon_k_sf[i]
441                        * rho_s
442                        * (sigma_sf[i].powi(12) / 90.
443                            * ((r_grid - 9.0 * pore_size)
444                                / (r_grid - pore_size).mapv(|x| x.powi(9))
445                                - (r_grid + 9.0 * pore_size)
446                                    / (r_grid + pore_size).mapv(|x| x.powi(9)))
447                            - sigma_sf[i].powi(6) / 3.
448                                * ((r_grid - 3.0 * pore_size)
449                                    / (r_grid - pore_size).mapv(|x| x.powi(3))
450                                    - (r_grid + 3.0 * pore_size)
451                                        / (r_grid + pore_size).mapv(|x| x.powi(3))))
452                        / r_grid
453                }
454                Self::SimpleLJ93 {
455                    sigma_ss: _,
456                    epsilon_k_ss: _,
457                } => {
458                    unimplemented!()
459                }
460                Self::CustomLJ93 {
461                    sigma_sf: _,
462                    epsilon_k_sf: _,
463                } => {
464                    unimplemented!()
465                }
466                Self::Steele {
467                    sigma_ss,
468                    epsilon_k_ss,
469                    rho_s,
470                    xi,
471                } => {
472                    // combining rules
473                    let epsilon_k_sf =
474                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
475                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
476
477                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
478                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
479                        * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
480                            - sum_n(4, r_grid, sigma_sf[i], pore_size)
481                            - sigma_sf[i] / (3.0 * DELTA_STEELE)
482                                * (sigma_sf[i].powi(3)
483                                    / r_grid
484                                        .mapv(|r| (pore_size + 0.61 * DELTA_STEELE - r).powi(3))
485                                    + sigma_sf[i].powi(3)
486                                        / r_grid.mapv(|r| {
487                                            (pore_size + 0.61 * DELTA_STEELE + r).powi(3)
488                                        })
489                                    + 1.5
490                                        * sum_n(
491                                            3,
492                                            r_grid,
493                                            sigma_sf[i],
494                                            pore_size + 0.61 * DELTA_STEELE,
495                                        )))
496                }
497                Self::CustomSteele {
498                    sigma_sf,
499                    epsilon_k_sf,
500                    rho_s,
501                    xi,
502                } => {
503                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
504                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
505                        * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
506                            - sum_n(4, r_grid, sigma_sf[i], pore_size)
507                            - sigma_sf[i] / (3.0 * DELTA_STEELE)
508                                * (sigma_sf[i].powi(3)
509                                    / r_grid
510                                        .mapv(|r| (pore_size + 0.61 * DELTA_STEELE - r).powi(3))
511                                    + sigma_sf[i].powi(3)
512                                        / r_grid.mapv(|r| {
513                                            (pore_size + 0.61 * DELTA_STEELE + r).powi(3)
514                                        })
515                                    + 1.5
516                                        * sum_n(
517                                            3,
518                                            r_grid,
519                                            sigma_sf[i],
520                                            pore_size + 0.61 * DELTA_STEELE,
521                                        )))
522                }
523                Self::DoubleWell {
524                    sigma_ss,
525                    epsilon1_k_ss,
526                    epsilon2_k_ss,
527                    rho_s,
528                } => {
529                    // combining rules
530                    let epsilon1_k_sf =
531                        (fluid_parameters.epsilon_k_ff() * *epsilon1_k_ss).map(|e| e.sqrt());
532                    let epsilon2_k_sf =
533                        (fluid_parameters.epsilon_k_ff() * *epsilon2_k_ss).map(|e| e.sqrt());
534                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;
535
536                    let bh_tail = (2.0
537                        * PI
538                        * mi
539                        * epsilon2_k_sf[i]
540                        * sigma_sf[i].powi(2)
541                        * rho_s
542                        * (2.0 / 5.0 * sum_n(10, r_grid, 2.0 * sigma_sf[i], pore_size)
543                            - sum_n(4, r_grid, 2.0 * sigma_sf[i], pore_size)))
544                    .mapv(|x| x.min(0.0));
545
546                    bh_tail
547                        + &(2.0
548                            * PI
549                            * mi
550                            * epsilon1_k_sf[i]
551                            * sigma_sf[i].powi(2)
552                            * rho_s
553                            * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
554                                - sum_n(4, r_grid, sigma_sf[i], pore_size)))
555                }
556                #[cfg(feature = "rayon")]
557                Self::FreeEnergyAveraged {
558                    coordinates,
559                    sigma_ss,
560                    epsilon_k_ss,
561                    pore_center,
562                    system_size,
563                    n_grid,
564                    cutoff_radius,
565                } => {
566                    // combining rules
567                    let epsilon_k_sf =
568                        (fluid_parameters.epsilon_k_ff()[i] * epsilon_k_ss).map(|e| e.sqrt());
569                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + sigma_ss) * 0.5;
570
571                    calculate_fea_potential(
572                        r_grid,
573                        mi,
574                        coordinates,
575                        sigma_sf,
576                        epsilon_k_sf,
577                        pore_center,
578                        system_size,
579                        n_grid,
580                        temperature,
581                        Geometry::Spherical,
582                        *cutoff_radius,
583                    )
584                }
585                _ => unreachable!(),
586            });
587        }
588        ext_pot
589    }
590}
591
592fn phi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
593    let m3n2 = 3.0 - 2.0 * n as f64;
594    let n2m3 = 2.0 * n as f64 - 3.0;
595
596    (1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(m3n2)) * 4.0 * PI.sqrt() / n2m3
597        * sigma_r.powf(n2m3)
598        * tgamma(n as f64 - 0.5)
599        / tgamma(n as f64)
600        * taylor_2f1_phi(r_r, n)
601}
602
603fn psi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
604    (1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(2.0 - 2.0 * n as f64))
605        * 4.0
606        * PI.sqrt()
607        * tgamma(n as f64 - 0.5)
608        / tgamma(n as f64)
609        * sigma_r.powf(2.0 * n as f64 - 2.0)
610        * taylor_2f1_psi(r_r, n)
611}
612
613fn sum_n(n: i32, r_grid: &Array1<f64>, sigma: f64, pore_size: f64) -> Array1<f64> {
614    let mut result = Array1::zeros(r_grid.len());
615    for i in 0..n {
616        result = result
617            + sigma.powi(n) / (pore_size.powi(i) * r_grid.mapv(|r| (pore_size - r).powi(n - i)))
618            + sigma.powi(n) / (pore_size.powi(i) * r_grid.mapv(|r| (pore_size + r).powi(n - i)));
619    }
620    result
621}
622
623fn taylor_2f1_phi(x: &Array1<f64>, n: i32) -> Array1<f64> {
624    match n {
625        3 => {
626            1.0 + 3.0 / 4.0 * x.mapv(|x| x.powi(2)) + 3.0 / 64.0 * x.mapv(|x| x.powi(4))
627                - 1.0 / 256.0 * x.mapv(|x| x.powi(6))
628                - 15.0 / 16384.0 * x.mapv(|x| x.powi(8))
629        }
630        6 => {
631            1.0 + 63.0 / 4.0 * x.mapv(|x| x.powi(2))
632                + 2205.0 / 64.0 * x.mapv(|x| x.powi(4))
633                + 3675.0 / 256.0 * x.mapv(|x| x.powi(6))
634                + 11025.0 / 16384.0 * x.mapv(|x| x.powi(8))
635        }
636        _ => unreachable!(),
637    }
638}
639
640fn taylor_2f1_psi(x: &Array1<f64>, n: i32) -> Array1<f64> {
641    match n {
642        3 => {
643            1.0 + 9.0 / 4.0 * x.mapv(|x| x.powi(2))
644                + 9.0 / 64.0 * x.mapv(|x| x.powi(4))
645                + 1.0 / 256.0 * x.mapv(|x| x.powi(6))
646                + 9.0 / 16384.0 * x.mapv(|x| x.powi(8))
647        }
648        6 => {
649            1.0 + 81.0 / 4.0 * x.mapv(|x| x.powi(2))
650                + 3969.0 / 64.0 * x.mapv(|x| x.powi(4))
651                + 11025.0 / 256.0 * x.mapv(|x| x.powi(6))
652                + 99125.0 / 16384.0 * x.mapv(|x| x.powi(8))
653        }
654        _ => unreachable!(),
655    }
656}
feos_dft/adsorption/external_potential.rs

feos_dft/adsorption/
external_potential.rs
