Skip to content

Commit 949ce5d

Browse files
prehnerprehner
committed
Partial derivatives of density profiles
1 parent 60403a1 commit 949ce5d

File tree

6 files changed

+259
-15
lines changed

6 files changed

+259
-15
lines changed

‎feos-dft/src/functional.rs‎

Lines changed: 37 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -441,6 +441,43 @@ impl<T: HelmholtzEnergyFunctional> DFT<T> {
441441
))
442442
}
443443

444+
#[allow(clippy::type_complexity)]
445+
pub(crate) fn functional_derivative_dual<D>(
446+
&self,
447+
temperature: f64,
448+
density: &Array<f64, D::Larger>,
449+
convolver: &Arc<dyn Convolver<Dual64, D>>,
450+
) -> EosResult<(Array<Dual64, D>, Array<Dual64, D::Larger>)>
451+
where
452+
D: Dimension,
453+
D::Larger: Dimension<Smaller = D>,
454+
{
455+
let temperature_dual = Dual64::from(temperature).derive();
456+
let density_dual = density.mapv(Dual64::from);
457+
let weighted_densities = convolver.weighted_densities(&density_dual);
458+
let contributions = self.contributions();
459+
let mut partial_derivatives = Vec::with_capacity(contributions.len());
460+
let mut helmholtz_energy_density = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
461+
for (c, wd) in contributions.iter().zip(weighted_densities) {
462+
let nwd = wd.shape()[0];
463+
let ngrid = wd.len() / nwd;
464+
let mut phi = Array::zeros(density.raw_dim().remove_axis(Axis(0)));
465+
let mut pd = Array::zeros(wd.raw_dim());
466+
c.first_partial_derivatives_dual(
467+
temperature_dual,
468+
wd.into_shape((nwd, ngrid)).unwrap(),
469+
phi.view_mut().into_shape(ngrid).unwrap(),
470+
pd.view_mut().into_shape((nwd, ngrid)).unwrap(),
471+
)?;
472+
partial_derivatives.push(pd);
473+
helmholtz_energy_density += &phi;
474+
}
475+
Ok((
476+
helmholtz_energy_density,
477+
convolver.functional_derivative(&partial_derivatives) * temperature_dual,
478+
))
479+
}
480+
444481
/// Calculate the bond integrals $I_{\alpha\alpha'}(\mathbf{r})$
445482
pub fn bond_integrals<D>(
446483
&self,

‎feos-dft/src/functional_contribution.rs‎

Lines changed: 28 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ use feos_core::{EosResult, HelmholtzEnergyDual, StateHD};
33
use ndarray::prelude::*;
44
use ndarray::RemoveAxis;
55
use num_dual::*;
6-
use num_traits::Zero;
6+
use num_traits::{One, Zero};
77
use std::fmt::Display;
88

99
macro_rules! impl_helmholtz_energy {
@@ -75,6 +75,7 @@ pub trait FunctionalContributionDual<N: DualNum<f64>>: Display {
7575
pub trait FunctionalContribution:
7676
FunctionalContributionDual<f64>
7777
+ FunctionalContributionDual<Dual64>
78+
+ FunctionalContributionDual<Dual<Dual64, f64>>
7879
+ FunctionalContributionDual<Dual<DualVec64<3>, f64>>
7980
+ FunctionalContributionDual<HyperDual64>
8081
+ FunctionalContributionDual<Dual2_64>
@@ -114,6 +115,31 @@ pub trait FunctionalContribution:
114115
Ok(())
115116
}
116117

118+
fn first_partial_derivatives_dual(
119+
&self,
120+
temperature: Dual64,
121+
weighted_densities: Array2<Dual64>,
122+
mut helmholtz_energy_density: ArrayViewMut1<Dual64>,
123+
mut first_partial_derivative: ArrayViewMut2<Dual64>,
124+
) -> EosResult<()> {
125+
let mut wd = weighted_densities.mapv(Dual::from_re);
126+
let t = Dual::from_re(temperature);
127+
let mut phi = Array::zeros(weighted_densities.raw_dim().remove_axis(Axis(0)));
128+
129+
for i in 0..wd.shape()[0] {
130+
wd.index_axis_mut(Axis(0), i)
131+
.map_inplace(|x| x.eps[0] = Dual64::one());
132+
phi = self.calculate_helmholtz_energy_density(t, wd.view())?;
133+
first_partial_derivative
134+
.index_axis_mut(Axis(0), i)
135+
.assign(&phi.mapv(|p| p.eps[0]));
136+
wd.index_axis_mut(Axis(0), i)
137+
.map_inplace(|x| x.eps[0] = Dual64::zero());
138+
}
139+
helmholtz_energy_density.assign(&phi.mapv(|p| p.re));
140+
Ok(())
141+
}
142+
117143
fn second_partial_derivatives(
118144
&self,
119145
temperature: f64,
@@ -161,6 +187,7 @@ pub trait FunctionalContribution:
161187
impl<T> FunctionalContribution for T where
162188
T: FunctionalContributionDual<f64>
163189
+ FunctionalContributionDual<Dual64>
190+
+ FunctionalContributionDual<Dual<Dual64, f64>>
164191
+ FunctionalContributionDual<Dual<DualVec64<3>, f64>>
165192
+ FunctionalContributionDual<HyperDual64>
166193
+ FunctionalContributionDual<Dual2_64>

‎feos-dft/src/geometry.rs‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -149,11 +149,9 @@ impl Axis {
149149
}
150150
let x0 = 0.5 * ((-alpha * points as f64).exp() + (-alpha * (points - 1) as f64).exp());
151151
let grid = (0..points)
152-
.into_iter()
153152
.map(|i| l * x0 * (alpha * i as f64).exp())
154153
.collect();
155154
let edges = (0..=points)
156-
.into_iter()
157155
.map(|i| {
158156
if i == 0 {
159157
0.0
@@ -166,7 +164,6 @@ impl Axis {
166164
let k0 = (2.0 * alpha).exp() * (2.0 * alpha.exp() + (2.0 * alpha).exp() - 1.0)
167165
/ ((1.0 + alpha.exp()).powi(2) * ((2.0 * alpha).exp() - 1.0));
168166
let integration_weights = (0..points)
169-
.into_iter()
170167
.map(|i| {
171168
(match i {
172169
0 => k0 * (2.0 * alpha).exp(),

‎feos-dft/src/profile.rs‎

Lines changed: 154 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3,14 +3,13 @@ use crate::functional::{HelmholtzEnergyFunctional, DFT};
33
use crate::geometry::Grid;
44
use crate::solver::{DFTSolver, DFTSolverLog};
55
use crate::weight_functions::WeightFunctionInfo;
6-
use feos_core::{Contributions, EosError, EosResult, EosUnit, EquationOfState, State};
6+
use feos_core::{Contributions, EosError, EosResult, EosUnit, EquationOfState, State, Verbosity};
77
use ndarray::{
88
Array, Array1, ArrayBase, Axis as Axis_nd, Data, Dimension, Ix1, Ix2, Ix3, RemoveAxis,
99
};
1010
use num_dual::Dual64;
11-
use quantity::si::{SIArray, SIArray1, SINumber, SIUnit};
11+
use quantity::si::{SIArray, SIArray1, SIArray2, SINumber, SIUnit};
1212
use quantity::Quantity;
13-
// use quantity::{Quantity, QuantityArray, SIArray1, SINumber};
1413
use std::ops::MulAssign;
1514
use std::sync::Arc;
1615

@@ -531,4 +530,156 @@ where
531530
)? * SIUnit::reference_pressure();
532531
Ok(self.integrate(&internal_energy_density))
533532
}
533+
534+
pub fn drho_dmu(&self) -> EosResult<SIArray<<D::Larger as Dimension>::Larger>> {
535+
let rho = self.density.to_reduced(SIUnit::reference_density())?;
536+
let second_partial_derivatives = self.second_partial_derivatives(&rho)?;
537+
538+
let shape = rho.shape();
539+
let shape: Vec<_> = std::iter::once(&shape[0]).chain(shape).copied().collect();
540+
let mut drho_dmu = Array::zeros(shape).into_dimensionality().unwrap();
541+
let rhs = |x: &_| {
542+
let delta_functional_derivative =
543+
self.delta_functional_derivative(x, &second_partial_derivatives);
544+
let mut xm = x.clone();
545+
xm.outer_iter_mut()
546+
.zip(self.dft.m().iter())
547+
.for_each(|(mut x, &m)| x *= m);
548+
// let delta_i = self.delta_bond_integrals(&exp_dfdrho, &delta_functional_derivative);
549+
// x + (delta_functional_derivative - delta_i) * rho
550+
xm + delta_functional_derivative * &rho
551+
};
552+
let mut log = DFTSolverLog::new(Verbosity::None);
553+
for (k, mut d) in drho_dmu.outer_iter_mut().enumerate() {
554+
let mut lhs = rho.clone();
555+
for (i, mut l) in lhs.outer_iter_mut().enumerate() {
556+
if i != k {
557+
l.fill(0.0);
558+
}
559+
}
560+
d.assign(&Self::gmres(rhs, &lhs, 200, 1e-13, &mut log)?);
561+
}
562+
Ok(drho_dmu
563+
* (SIUnit::reference_density() / SIUnit::reference_molar_entropy() / self.temperature))
564+
}
565+
566+
pub fn dn_dmu(&self) -> EosResult<SIArray2> {
567+
let drho_dmu = self.drho_dmu()?;
568+
let n = drho_dmu.shape()[0];
569+
let dn_dmu = SIArray2::from_shape_fn([n; 2], |(i, j)| {
570+
self.integrate(&drho_dmu.index_axis(Axis_nd(0), i).index_axis(Axis_nd(0), j))
571+
});
572+
Ok(dn_dmu)
573+
}
574+
575+
pub fn drho_dp(&self) -> EosResult<SIArray<D::Larger>> {
576+
let rho = self.density.to_reduced(SIUnit::reference_density())?;
577+
let partial_density = self
578+
.bulk
579+
.partial_density
580+
.to_reduced(SIUnit::reference_density())?;
581+
let rho_bulk = self.dft.component_index().mapv(|i| partial_density[i]);
582+
583+
let second_partial_derivatives = self.second_partial_derivatives(&rho)?;
584+
let (_, _, _, exp_dfdrho, _) = self.euler_lagrange_equation(&rho, &rho_bulk, false)?;
585+
586+
let rhs = |x: &_| {
587+
let delta_functional_derivative =
588+
self.delta_functional_derivative(x, &second_partial_derivatives);
589+
let mut xm = x.clone();
590+
xm.outer_iter_mut()
591+
.zip(self.dft.m().iter())
592+
.for_each(|(mut x, &m)| x *= m);
593+
let delta_i = self.delta_bond_integrals(&exp_dfdrho, &delta_functional_derivative);
594+
xm + (delta_functional_derivative - delta_i) * &rho
595+
};
596+
let mut lhs = rho.clone();
597+
let v = self
598+
.bulk
599+
.partial_molar_volume(Contributions::Total)
600+
.to_reduced(SIUnit::reference_volume() / SIUnit::reference_moles())?;
601+
for (mut l, &c) in lhs.outer_iter_mut().zip(self.dft.component_index().iter()) {
602+
l *= v[c];
603+
}
604+
let mut log = DFTSolverLog::new(Verbosity::None);
605+
Self::gmres(rhs, &lhs, 200, 1e-13, &mut log)
606+
.map(|x| x / (SIUnit::reference_molar_entropy() * self.temperature))
607+
}
608+
609+
pub fn dn_dp(&self) -> EosResult<SIArray1> {
610+
let dn_dp = self.integrate_comp(&self.drho_dp()?);
611+
let mut dn_dp_comp = Array1::zeros(self.dft.components()) * SIUnit::reference_moles()
612+
/ SIUnit::reference_pressure();
613+
for (i, &j) in self.dft.component_index().iter().enumerate() {
614+
dn_dp_comp.try_set(j, dn_dp.get(i)).unwrap();
615+
}
616+
Ok(dn_dp_comp)
617+
}
618+
619+
pub fn drho_dt(&self) -> EosResult<SIArray<D::Larger>> {
620+
let rho = self.density.to_reduced(SIUnit::reference_density())?;
621+
let t = self
622+
.temperature
623+
.to_reduced(SIUnit::reference_temperature())?;
624+
625+
let second_partial_derivatives = self.second_partial_derivatives(&rho)?;
626+
627+
// define rhs
628+
let rhs = |x: &_| {
629+
let delta_functional_derivative =
630+
self.delta_functional_derivative(x, &second_partial_derivatives);
631+
let mut xm = x.clone();
632+
xm.outer_iter_mut()
633+
.zip(self.dft.m().iter())
634+
.for_each(|(mut x, &m)| x *= m);
635+
// let delta_i = self.delta_bond_integrals(&exp_dfdrho, &delta_functional_derivative);
636+
// (xm + (delta_functional_derivative - delta_i) * &rho) * t
637+
(xm + delta_functional_derivative * &rho) * t
638+
};
639+
640+
// calculate temperature derivative of functional derivative
641+
let functional_contributions = self.dft.contributions();
642+
let weight_functions: Vec<WeightFunctionInfo<Dual64>> = functional_contributions
643+
.iter()
644+
.map(|c| c.weight_functions(Dual64::from(t).derive()))
645+
.collect();
646+
let convolver: Arc<dyn Convolver<_, D>> =
647+
ConvolverFFT::plan(&self.grid, &weight_functions, None);
648+
let (_, dfdrhodt) = self.dft.functional_derivative_dual(t, &rho, &convolver)?;
649+
650+
// calculate temperature derivative of bulk functional derivative
651+
let partial_density = self
652+
.bulk
653+
.partial_density
654+
.to_reduced(SIUnit::reference_density())?;
655+
let rho_bulk = self.dft.component_index().mapv(|i| partial_density[i]);
656+
let bulk_convolver = BulkConvolver::new(weight_functions);
657+
let (_, dfdrhodt_bulk) =
658+
self.dft
659+
.functional_derivative_dual(t, &rho_bulk, &bulk_convolver)?;
660+
661+
// solve for drho_dt
662+
let x = (self.bulk.dp_dni(Contributions::Total) * self.bulk.dp_dt(Contributions::Total)
663+
/ self.bulk.dp_dv(Contributions::Total))
664+
.to_reduced(SIUnit::reference_molar_entropy())?;
665+
let mut lhs = dfdrhodt.mapv(|d| d.eps[0]);
666+
lhs.outer_iter_mut()
667+
.zip(dfdrhodt_bulk.into_iter())
668+
.zip(x.into_iter())
669+
.for_each(|((mut lhs, d), x)| lhs -= d.eps[0] + x);
670+
lhs.outer_iter_mut()
671+
.zip(rho.outer_iter())
672+
.zip(rho_bulk.into_iter())
673+
.zip(self.dft.m().iter())
674+
.for_each(|(((mut lhs, rho), rho_b), &m)| lhs += &((&rho / rho_b).mapv(f64::ln) * m));
675+
676+
lhs *= &(-&rho);
677+
let mut log = DFTSolverLog::new(Verbosity::None);
678+
let drho_dt = Self::gmres(rhs, &lhs, 200, 1e-13, &mut log)?;
679+
Ok(drho_dt * (SIUnit::reference_density() / SIUnit::reference_temperature()))
680+
}
681+
682+
pub fn dn_dt(&self) -> EosResult<SIArray1> {
683+
Ok(self.integrate_comp(&self.drho_dt()?))
684+
}
534685
}

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

Lines changed: 34 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
#[macro_export]
22
macro_rules! impl_profile {
3-
($struct:ident, $arr:ident, $arr2:ident, $si_arr:ident, $si_arr2:ident, $py_arr2:ident, [$([$ind:expr, $ax:ident]),+]) => {
3+
($struct:ident, $arr:ident, $arr2:ident, $si_arr:ident, $si_arr2:ident, $py_arr2:ident, [$([$ind:expr, $ax:ident]),+]$(, $si_arr3:ident)?) => {
44
#[pymethods]
55
impl $struct {
66
/// Calculate the residual for the given profile.
@@ -178,6 +178,37 @@ macro_rules! impl_profile {
178178
self.0.profile.grand_potential_density()?,
179179
))
180180
}
181+
$(
182+
#[getter]
183+
fn get_drho_dmu(&self) -> PyResult<$si_arr3> {
184+
Ok(($si_arr3::from(self.0.profile.drho_dmu()?)))
185+
}
186+
)?
187+
188+
#[getter]
189+
fn get_dn_dmu(&self) -> PyResult<PySIArray2> {
190+
Ok((PySIArray2::from(self.0.profile.dn_dmu()?)))
191+
}
192+
193+
#[getter]
194+
fn get_drho_dp(&self) -> PyResult<$si_arr2> {
195+
Ok(($si_arr2::from(self.0.profile.drho_dp()?)))
196+
}
197+
198+
#[getter]
199+
fn get_dn_dp(&self) -> PyResult<PySIArray1> {
200+
Ok((PySIArray1::from(self.0.profile.dn_dp()?)))
201+
}
202+
203+
#[getter]
204+
fn get_drho_dt(&self) -> PyResult<$si_arr2> {
205+
Ok(($si_arr2::from(self.0.profile.drho_dt()?)))
206+
}
207+
208+
#[getter]
209+
fn get_dn_dt(&self) -> PyResult<PySIArray1> {
210+
Ok((PySIArray1::from(self.0.profile.dn_dt()?)))
211+
}
181212
}
182213
};
183214
}
@@ -192,7 +223,8 @@ macro_rules! impl_1d_profile {
192223
PySIArray1,
193224
PySIArray2,
194225
PyArray2,
195-
[$([0, $ax]),+]
226+
[$([0, $ax]),+],
227+
PySIArray3
196228
);
197229
};
198230
}

0 commit comments

Comments
 (0)