Fix 007_PW_UPF201_USPP_Fe in 01_PW: G=0 indexing in nspin=2 mixing inner products

source/source_estate/module_charge/mix_resid.cpp

@@ -171,14 +171,23 @@ double Charge_Mixing::inner_product_recip_rho(std::complex<double>* rho1, std::c
         // including |G|=0 term.
         double sum2 = 0.0;
 
-        sum2 += fac2 * (conj(rhog1[0][0] - rhog1[1][0]) * (rhog2[0][0] - rhog2[1][0])).real();
+        // The G=0 component is the ig_gge0-th element of the local G-list on the
+        // rank that owns it, not necessarily element 0: the local G-list is built
+        // by scanning (x,y) sticks in grid order, so element 0 is the first plane
+        // wave of the first owned stick. Using a hardcoded index 0 made the inner
+        // product partition-dependent for pools with more than one rank.
+        if (ig0 >= 0)
+        {
+            sum2 += fac2 * (conj(rhog1[0][ig0] - rhog1[1][ig0]) * (rhog2[0][ig0] - rhog2[1][ig0])).real();
+        }
 
         double mag = 0.0;
 #ifdef _OPENMP
 #pragma omp parallel for reduction(+ : mag)
 #endif
         for (int ig = 0; ig < this->rhopw->npw; ig++)
         {
+            if (ig == ig0) { continue; }
             mag += (conj(rhog1[0][ig] - rhog1[1][ig]) * (rhog2[0][ig] - rhog2[1][ig])).real();
         }
         mag *= fac2;
@@ -342,14 +351,21 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex<double>* rhog1, s
         // including |G|=0 term.
         double sum2 = 0.0;
 
-        sum2 += fac2 * (conj(rhog1[0 + this->rhopw->npw]) * rhog2[0 + this->rhopw->npw]).real();
+        // Same G=0 indexing remark as in inner_product_recip_rho: use ig_gge0
+        // instead of a hardcoded index 0, otherwise the inner product (and hence
+        // the DIIS mixing coefficients) depends on how the pool is divided.
+        if (ig0 >= 0)
+        {
+            sum2 += fac2 * (conj(rhog1[ig0 + this->rhopw->npw]) * rhog2[ig0 + this->rhopw->npw]).real();
+        }
 
         double mag = 0.0;
 #ifdef _OPENMP
 #pragma omp parallel for reduction(+ : mag)
 #endif
         for (int ig = 0; ig < this->rhopw->npw; ig++)
         {
+            if (ig == ig0) { continue; }
             mag += (conj(rhog1[ig + this->rhopw->npw]) * rhog2[ig + this->rhopw->npw]).real();
         }
         mag *= fac2;

tests/01_PW/007_PW_UPF201_USPP_Fe/result.ref

@@ -1,5 +1,5 @@
-etotref -673.8349306935572258
-etotperatomref -673.8349306936
+etotref -673.8349346996149052
+etotperatomref -673.8349346996
 totalforceref 0.000000
-totalstressref 66620.323844
-totaltimeref 1.34
+totalstressref 66620.327426
+totaltimeref 1.18

tests/01_PW/096_PW_PBE0_AFM/threshold

@@ -0,0 +1 @@
+threshold 0.001

tests/01_PW/CASES_GPU.txt

@@ -8,7 +8,7 @@ scf_out_elf
 004_PW_UPF201_Si
 005_PW_UPF201_UPF100
 006_PW_UPF201_Eu
-#007_PW_UPF201_USPP_Fe
+007_PW_UPF201_USPP_Fe
 008_PW_UPF201_USPP_NaCl
 009_PW_UPF201_USPP
 010_PW_0TYPE