Fix the sign of vector weighted densities in cartesian and cylindrical geometry (feos-org#120)

prehner · web-flow · commit 0cfccdbcc39a · 2023-01-23T17:08:23.000+01:00
* Fix the sign of vector weighted densities in cartesian and cylindrical geometry

* add changelog entry
diff --git 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
@@ -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<T: DualNum<f64> + DctNum + ScalarOperand> FourierTransform<T> 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<T: DualNum<f64> + DctNum + ScalarOperand> FourierTransform<T> 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
@@ -116,7 +116,7 @@ impl<T: DualNum<f64>> WeightFunction<T> {
                 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!(),

Original file line number	Diff line number	Diff line change
`@@ -116,7 +116,7 @@ impl<T: DualNum<f64>> WeightFunction<T> {`
`116`	`116`	`w_i.assign(&match self.shape {`
`117`	`117`	`WeightFunctionShape::DeltaVec => {`
`118`	`118`	`&rik.mapv(\|rik\| {`
`119`		`- (rik.sph_j0() + rik.sph_j2()) * (-radius.powi(3) * 4.0 * FRAC_PI_3 * p)`
	`119`	`+ (rik.sph_j0() + rik.sph_j2()) * (radius.powi(3) * 4.0 * FRAC_PI_3 * p)`
`120`	`120`	`}) * &k_x`
`121`	`121`	`}`
`122`	`122`	`_ => unreachable!(),`