Fix missing last term in B-spline factor for structure factor

source/source_pw/module_pwdft/stru_fac.cpp

@@ -311,7 +311,7 @@ void Structure_Factor::bsplinecoef(std::complex<double> *b1, std::complex<double
     for(int ix = 0 ; ix < nx ; ++ix)
     {
         std::complex<double> fracx=0;
-        for(int io = 0 ; io < norder - 1 ; ++io)
+        for(int io = 0 ; io < norder ; ++io)
         {
             fracx += bsp.bezier_ele(io)*ModuleBase::libm::exp(ci_tpi*double(ix)/double(nx)*double(io));
         }
@@ -323,7 +323,7 @@ void Structure_Factor::bsplinecoef(std::complex<double> *b1, std::complex<double
     for(int iy = 0 ; iy < ny ; ++iy)
     {
         std::complex<double> fracy=0;
-        for(int io = 0 ; io < norder - 1 ; ++io)
+        for(int io = 0 ; io < norder ; ++io)
         {
             fracy += bsp.bezier_ele(io)*ModuleBase::libm::exp(ci_tpi*double(iy)/double(ny)*double(io));
         }
@@ -335,7 +335,7 @@ void Structure_Factor::bsplinecoef(std::complex<double> *b1, std::complex<double
     for(int iz = 0 ; iz < nz ; ++iz)
     {
         std::complex<double> fracz=0;
-        for(int io = 0 ; io < norder - 1 ; ++io)
+        for(int io = 0 ; io < norder ; ++io)
         {
             fracz += bsp.bezier_ele(io)*ModuleBase::libm::exp(ci_tpi*double(iz)/double(nz)*double(io));
         }

source/source_pw/module_pwdft/test/structure_factor_test.cpp

@@ -3,6 +3,7 @@
 #include <string>
 #include <cmath>
 #include <complex>
+#include <vector>
 #include "source_cell/unitcell.h"
 #include "source_estate/module_dm/unittests/prepare_unitcell.h"
 #include "source_io/module_parameter/parameter.h"
@@ -114,6 +115,34 @@ TEST_F(StructureFactorTest, setup_structure_factor_float)
     }
 }
 
+TEST_F(StructureFactorTest, bsplinecoef_zero_g_component)
+{
+    // The deconvolution factor at G = 0 must be exactly 1 for any spline
+    // order: its denominator sums the Cardinal B-spline values at integer
+    // points, which is a partition of unity.
+    // Regression test for a missing last term M_{n+1}(n) = 1/n! in that
+    // sum, which scaled S(G = 0) -- and hence the integrated electron
+    // number -- by (1 - 1/n!)^-3.
+    const int nx = 8;
+    const int ny = 12;
+    const int nz = 16;
+    std::vector<std::complex<double>> b1(nx);
+    std::vector<std::complex<double>> b2(ny);
+    std::vector<std::complex<double>> b3(nz);
+    // norder must be even: Bspline::init asserts norder % 2 == 0, and the
+    // production call site rounds odd nbspline values up to even.
+    for (const int norder : {2, 4, 6, 10, 20})
+    {
+        SF.bsplinecoef(b1.data(), b2.data(), b3.data(), nx, ny, nz, norder);
+        EXPECT_NEAR(b1[0].real(), 1.0, 1e-14);
+        EXPECT_NEAR(b1[0].imag(), 0.0, 1e-14);
+        EXPECT_NEAR(b2[0].real(), 1.0, 1e-14);
+        EXPECT_NEAR(b2[0].imag(), 0.0, 1e-14);
+        EXPECT_NEAR(b3[0].real(), 1.0, 1e-14);
+        EXPECT_NEAR(b3[0].imag(), 0.0, 1e-14);
+    }
+}
+
 int main()
 {
     testing::InitGoogleTest();

tests/03_NAO_multik/scf_bspline/result.ref

@@ -1,8 +1,8 @@
-etotref -4883.8847047255100000
-etotperatomref -2441.9423523628
-totalforceref 771.961380
-totalstressref 60277.376778
+etotref -4883.8848761152030420
+etotperatomref -2441.9424380576
+totalforceref 771.960378
+totalstressref 60277.381032
 pointgroupref C_3v
 spacegroupref C_3v
 nksibzref 4
-totaltimeref +0.65073
+totaltimeref 0.82

tests/07_OFDFT/25_OF_MD_LDA/result.ref

@@ -1,5 +1,5 @@
-etotref -231.6926210634764800
-etotperatomref -57.9231552659
-totalforceref 0.787244
-totalstressref 107.895762
-totaltimeref +1.49612
+etotref -231.6927874142393
+etotperatomref -57.9231968536
+totalforceref 0.787242
+totalstressref 107.899521
+totaltimeref 0.75

tests/07_OFDFT/26_OF_MD_LibxcPBE/result.ref

@@ -1,5 +1,5 @@
-etotref -228.7238739795837716
-etotperatomref -57.1809684949
-totalforceref 0.924446
-totalstressref 28.586116
-totaltimeref +3.98101
+etotref -228.7240427619324
+etotperatomref -57.1810106905
+totalforceref 0.924444
+totalstressref 28.582218
+totaltimeref 1.83

tests/07_OFDFT/30_TDOFDFT_Al/result.ref

@@ -1,5 +1,5 @@
-etotref 226.9979568201541724
-etotperatomref 113.4989784101
-totalforceref 562.524006
-totalstressref 194931.772028
-totaltimeref 0.45
+etotref 226.9983065043999
+etotperatomref 113.4991532522
+totalforceref 562.524754
+totalstressref 194931.658778
+totaltimeref 0.06

tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref

@@ -1,5 +1,5 @@
-etotref 225.6175831066692297
-etotperatomref 112.8087915533
-totalforceref 488.530072
-totalstressref 185999.479726
-totaltimeref 0.09
+etotref 225.6177626465156
+etotperatomref 112.8088813233
+totalforceref 488.530464
+totalstressref 185999.439893
+totaltimeref 0.06

tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref

@@ -1,5 +1,5 @@
-etotref 223.1193797466232809
-etotperatomref 111.5596898733
-totalforceref 386.437466
-totalstressref 181323.311310
-totaltimeref 0.08
+etotref 223.1195553368393
+etotperatomref 111.5597776684
+totalforceref 386.436746
+totalstressref 181323.302533
+totaltimeref 0.06

tests/07_OFDFT/35_OF_LPS_bspline/INPUT

@@ -0,0 +1,30 @@
+INPUT_PARAMETERS 
+#Parameters (1.General)
+suffix            autotest
+calculation       scf
+esolver_type      ofdft
+
+symmetry          1
+pseudo_dir        ../../PP_ORB
+pseudo_rcut       10
+nspin             1
+cal_force         1
+test_force        1
+cal_stress        1
+test_stress       1
+
+#Parameters (2.Iteration)
+ecutwfc           20
+scf_nmax          50
+
+#OFDFT
+of_kinetic        wt
+of_method         tn
+of_conv           energy
+of_tole           2e-10
+of_full_pw        1
+of_full_pw_dim    1
+nbspline          2
+
+#Parameters (3.Basis)
+basis_type        pw

tests/07_OFDFT/35_OF_LPS_bspline/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0

tests/07_OFDFT/35_OF_LPS_bspline/README

@@ -0,0 +1 @@
+Test the energy, force, and stress of OFDFT with Cardinal B-spline interpolation for the structure factor (nbspline = 2), using the Wang-Teter (WT) kinetic energy functional with LPS, symmetry=on; of_full_pw = 1 and of_full_pw_dim = 1 are set to keep the FFT dimensions odd, as required by the interpolation.

tests/07_OFDFT/35_OF_LPS_bspline/STRU

@@ -0,0 +1,18 @@
+ATOMIC_SPECIES
+Al 26.98 al_lps.cpi blps
+
+LATTICE_CONSTANT
+7.50241114482312  // add lattice constant
+
+LATTICE_VECTORS
+0.000000000000    0.500000000000    0.500000000000    
+0.500000000000    0.000000000000    0.500000000000    
+0.500000000000    0.500000000000    0.000000000000    
+
+ATOMIC_POSITIONS
+Direct
+
+Al
+0
+1
+    0.000000000000    0.000000000000    0.000000000000 1 1 1 

tests/07_OFDFT/35_OF_LPS_bspline/result.ref

@@ -0,0 +1,8 @@
+etotref -55.9081453854524568
+etotperatomref -55.9081453855
+totalforceref 0.000000
+totalstressref 282.762177
+pointgroupref O_h
+spacegroupref O_h
+nksibzref 1
+totaltimeref 0.40

tests/07_OFDFT/CASES_CPU.txt

@@ -31,4 +31,5 @@
 31_TDOFDFT_Al_CD
 32_TDOFDFT_Al_mCD
 33_OF_out_chg
-34_OF_KE_extWT
+34_OF_KE_extWT
+35_OF_LPS_bspline