diff --git a/feos-dft/CHANGELOG.md b/feos-dft/CHANGELOG.md index 889275978..1fa46464d 100644 --- a/feos-dft/CHANGELOG.md +++ b/feos-dft/CHANGELOG.md @@ -17,6 +17,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Renamed the parameter `beta` of the Picard iteration and Anderson mixing solvers to `damping_coefficient`. [#75](https://github.com/feos-org/feos/pull/75) - Removed generics for units in all structs and traits in favor of static SI units. [#115](https://github.com/feos-org/feos/pull/115) +### Fixed +- Fixed the sign of vector weighted densities in Cartesian and cylindrical geometries. The wrong sign had no effects on the Helmholtz energy and thus on density profiles. [#120](https://github.com/feos-org/feos/pull/120) + ## [0.3.2] - 2022-10-13 ### Changed - 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) diff --git a/feos-dft/src/convolver/transform.rs b/feos-dft/src/convolver/transform.rs index de4a2379a..0b9d082cb 100644 --- a/feos-dft/src/convolver/transform.rs +++ b/feos-dft/src/convolver/transform.rs @@ -5,7 +5,7 @@ use num_dual::*; use rustdct::{DctNum, DctPlanner, TransformType2And3}; use rustfft::{num_complex::Complex, Fft, FftPlanner}; use std::f64::consts::PI; -use std::ops::{DivAssign, SubAssign}; +use std::ops::DivAssign; use std::sync::Arc; #[derive(Clone, Copy)] @@ -189,10 +189,10 @@ impl + DctNum + ScalarOperand> FourierTransform for Spherical if scalar { self.sine_transform(&f_r * &self.r_grid, f_k.view_mut(), false); } else { - self.cosine_transform(&f_r * &self.r_grid, f_k.view_mut(), false); let mut f_aux = Array::zeros(f_k.raw_dim()); - self.sine_transform(f_r, f_aux.view_mut(), false); - f_k.sub_assign(&(&f_aux / &self.k_grid)); + self.cosine_transform(&f_r * &self.r_grid, f_aux.view_mut(), false); + self.sine_transform(f_r, f_k.view_mut(), false); + f_k.assign(&(&f_k / &self.k_grid - &f_aux)); } f_k.assign(&(&f_k / &self.k_grid)); f_k[0] = T::zero(); @@ -202,10 +202,10 @@ impl + DctNum + ScalarOperand> FourierTransform for Spherical if scalar { self.sine_transform(&f_k * &self.k_grid, f_r.view_mut(), true); } else { - self.cosine_transform(&f_k * &self.k_grid, f_r.view_mut(), true); let mut f_aux = Array::zeros(f_r.raw_dim()); - self.sine_transform(f_k, f_aux.view_mut(), true); - f_r.sub_assign(&(&f_aux / &self.r_grid)); + self.cosine_transform(&f_k * &self.k_grid, f_aux.view_mut(), true); + self.sine_transform(f_k, f_r.view_mut(), true); + f_r.assign(&(&f_r / &self.r_grid - &f_aux)); } f_r.assign(&(&f_r / &self.r_grid)); } diff --git a/feos-dft/src/weight_functions.rs b/feos-dft/src/weight_functions.rs index 9184c8e01..871a0b59d 100644 --- a/feos-dft/src/weight_functions.rs +++ b/feos-dft/src/weight_functions.rs @@ -116,7 +116,7 @@ impl> WeightFunction { w_i.assign(&match self.shape { WeightFunctionShape::DeltaVec => { &rik.mapv(|rik| { - (rik.sph_j0() + rik.sph_j2()) * (-radius.powi(3) * 4.0 * FRAC_PI_3 * p) + (rik.sph_j0() + rik.sph_j2()) * (radius.powi(3) * 4.0 * FRAC_PI_3 * p) }) * &k_x } _ => unreachable!(),