From 3e566d4323f05dfe36df74919ea22fb7c19631fe Mon Sep 17 00:00:00 2001 From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com> Date: Fri, 23 Oct 2026 20:15:09 +0800 Subject: [PATCH] Fix: missing last B-spline term in bsplinecoef scales electron number by (1-1/N!)^-3 The denominator loop in Structure_Factor::bsplinecoef missed the last integer-point sample M_{N+1}(N) = 1/N! of the Cardinal B-spline, so the G=0 deconvolution factor was 1/(1-1/N!) per direction instead of 1, and S(G=0) -- the integrated electron number -- was scaled by (1-1/N!)^-3 (8x at nbspline=2). The method is exact at G=0 at any order by partition of unity, so this was an off-by-one implementation error, not low-order interpolation error. Add a unit test asserting b(0)=1 for norder = 2,4,6,10,20 and a low-order OFDFT integration case 35_OF_LPS_bspline; regenerate references of the six existing nbspline=10 cases (energies shift by ~1e-4 eV, the intended behavior change). --- source/source_pw/module_pwdft/stru_fac.cpp | 6 ++-- .../test/structure_factor_test.cpp | 29 ++++++++++++++++++ tests/03_NAO_multik/scf_bspline/result.ref | 10 +++---- tests/07_OFDFT/25_OF_MD_LDA/result.ref | 10 +++---- tests/07_OFDFT/26_OF_MD_LibxcPBE/result.ref | 10 +++---- tests/07_OFDFT/30_TDOFDFT_Al/result.ref | 10 +++---- tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref | 10 +++---- tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref | 10 +++---- tests/07_OFDFT/35_OF_LPS_bspline/INPUT | 30 +++++++++++++++++++ tests/07_OFDFT/35_OF_LPS_bspline/KPT | 4 +++ tests/07_OFDFT/35_OF_LPS_bspline/README | 1 + tests/07_OFDFT/35_OF_LPS_bspline/STRU | 18 +++++++++++ tests/07_OFDFT/35_OF_LPS_bspline/result.ref | 8 +++++ tests/07_OFDFT/CASES_CPU.txt | 3 +- 14 files changed, 125 insertions(+), 34 deletions(-) create mode 100644 tests/07_OFDFT/35_OF_LPS_bspline/INPUT create mode 100644 tests/07_OFDFT/35_OF_LPS_bspline/KPT create mode 100644 tests/07_OFDFT/35_OF_LPS_bspline/README create mode 100644 tests/07_OFDFT/35_OF_LPS_bspline/STRU create mode 100644 tests/07_OFDFT/35_OF_LPS_bspline/result.ref diff --git a/source/source_pw/module_pwdft/stru_fac.cpp b/source/source_pw/module_pwdft/stru_fac.cpp index 20531268038..13526b7f4b6 100644 --- a/source/source_pw/module_pwdft/stru_fac.cpp +++ b/source/source_pw/module_pwdft/stru_fac.cpp @@ -311,7 +311,7 @@ void Structure_Factor::bsplinecoef(std::complex *b1, std::complex 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 *b1, std::complex 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 *b1, std::complex 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)); } diff --git a/source/source_pw/module_pwdft/test/structure_factor_test.cpp b/source/source_pw/module_pwdft/test/structure_factor_test.cpp index 129e398a77f..720a1bea2d3 100644 --- a/source/source_pw/module_pwdft/test/structure_factor_test.cpp +++ b/source/source_pw/module_pwdft/test/structure_factor_test.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #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> b1(nx); + std::vector> b2(ny); + std::vector> 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(); diff --git a/tests/03_NAO_multik/scf_bspline/result.ref b/tests/03_NAO_multik/scf_bspline/result.ref index 1d18434a3bb..df012bab01c 100644 --- a/tests/03_NAO_multik/scf_bspline/result.ref +++ b/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 diff --git a/tests/07_OFDFT/25_OF_MD_LDA/result.ref b/tests/07_OFDFT/25_OF_MD_LDA/result.ref index 0a65a37426c..b9ee156d7b5 100644 --- a/tests/07_OFDFT/25_OF_MD_LDA/result.ref +++ b/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 diff --git a/tests/07_OFDFT/26_OF_MD_LibxcPBE/result.ref b/tests/07_OFDFT/26_OF_MD_LibxcPBE/result.ref index 0869b1807d1..8367346d692 100644 --- a/tests/07_OFDFT/26_OF_MD_LibxcPBE/result.ref +++ b/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 diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/result.ref b/tests/07_OFDFT/30_TDOFDFT_Al/result.ref index b53afa1a10e..042a7e812aa 100644 --- a/tests/07_OFDFT/30_TDOFDFT_Al/result.ref +++ b/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 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref b/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref index 286a081395a..50f83a8774f 100644 --- a/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref +++ b/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 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref index a6e3cab7f3a..73db58e0e9e 100644 --- a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref +++ b/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 diff --git a/tests/07_OFDFT/35_OF_LPS_bspline/INPUT b/tests/07_OFDFT/35_OF_LPS_bspline/INPUT new file mode 100644 index 00000000000..e29e3dee920 --- /dev/null +++ b/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 diff --git a/tests/07_OFDFT/35_OF_LPS_bspline/KPT b/tests/07_OFDFT/35_OF_LPS_bspline/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/07_OFDFT/35_OF_LPS_bspline/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/07_OFDFT/35_OF_LPS_bspline/README b/tests/07_OFDFT/35_OF_LPS_bspline/README new file mode 100644 index 00000000000..d49a13fab27 --- /dev/null +++ b/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. diff --git a/tests/07_OFDFT/35_OF_LPS_bspline/STRU b/tests/07_OFDFT/35_OF_LPS_bspline/STRU new file mode 100644 index 00000000000..51c8570791a --- /dev/null +++ b/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 diff --git a/tests/07_OFDFT/35_OF_LPS_bspline/result.ref b/tests/07_OFDFT/35_OF_LPS_bspline/result.ref new file mode 100644 index 00000000000..0fb9395432f --- /dev/null +++ b/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 diff --git a/tests/07_OFDFT/CASES_CPU.txt b/tests/07_OFDFT/CASES_CPU.txt index d58f16da569..a3a7b5f7a4e 100644 --- a/tests/07_OFDFT/CASES_CPU.txt +++ b/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 \ No newline at end of file +34_OF_KE_extWT +35_OF_LPS_bspline