Stabilize RI RPA/EXX and Ewald EXX method merge on develop (deepmodeling#7019)

* Stabilize RI RPA/EXX merge on develop

- integrate the RI/RPA post-SCF workflow with develop

- restore Coulomb/input-output handling and sparse matrix precision

- fix nspin=2 and TRS symmetry output regressions

* Fix read_input test build regressions

* Fix restart include in input conversion

* Register new input parameters

* Fix complie errors

* fix: restore latname default to "user_defined_lattice"

The default was accidentally changed to "none", causing all tests
without explicit latname to fail at STRU parsing with
"do not use LATTICE_VECTORS along with explicit specification of lattice type".

* fix: restore EXX parameter registration from develop, append RI params

The previous merge accidentally changed exx_singularity_correction
reset logic (lc_pbe etc from "spencer" to "massidda") and removed
metadata, causing 08_EXX integration tests to fail. Restore the
develop version and append new RI parameters (shrink_abfs_pca_thr,
shrink_lu_inv_thr, exx_coul_moment, exx_rotate_abfs,
exx_multip_moments_threshold, out_unshrinked_v).

* fix: restore filter_empty_orbs in Exx_LRI::init to fix 53_GO_LR_HF test

The removal of filter_empty_orbs calls caused empty orbital shells to be
retained in lcaos and abfs vectors, altering RI exchange integrals for
all EXX calculations including standard HF used by LR-TDDFT.

* fix: inline filter_empty_orbs for abfs to fix build error

The PR removed the filter_empty_orbs function from Construct_Orbs and
inlined its logic into change_orbs() for lcaos. The previous commit
incorrectly called the now-deleted function. This commit inlines the
same filtering logic directly for abfs to restore the empty orbital
shell removal without depending on the deleted function.

* fix: restore RI infrastructure files to develop and fix spencer Rcut formula

- Restore LRI_CV, Matrix_Orbs11/21/22, exx_abfs-construct_orbs,
  exx_opt_orb, ABFs_Construct-PCA to develop versions to minimize diff
- Fix spencer Rcut: use get_nkstot_full() without /nspin0 division
  (nkstot_full already excludes spin doubling)
- Fix shared_ptr<ORB_gaunt_table> usage in Exx_LRI, RPA_LRI
- Fix ewald_Vq compilation errors
- Fix exx_opt_orb-print.cpp to match develop API

* fix: update 53_GO_LR_HF reference value for corrected spencer Rcut

The spencer Rcut formula was corrected to not divide by nspin0
(get_nkstot_full() already excludes spin doubling), which changes
the excitation energy result from 6.050430 to 2.574437.

* fix: update 53_GO_LR_HF ref to CI-verified value 5.646360

The previous ref (2.574437) was from fisherd which lacks LibXC,
producing different results. Use CI environment's actual value.

* Fix shrink RPA overlap setup and add regression coverage

* Compare RPA regression files by absolute value

* Reduce solid RPA regression output size

* Drop no-op LRI_CV_Tools formatting changes

* Drop no-op Inverse_Matrix formatting change

---------

Co-authored-by: AroundPeking <gonghuanjing@iphy.ac.cn>
Co-authored-by: Mohan Chen <mohanchen@pku.edu.cn>

Author: 3 people

docs/advanced/input_files/input-main.md

@@ -301,6 +301,8 @@
     - [exx\_pca\_threshold](#exx_pca_threshold)
     - [exx\_c\_threshold](#exx_c_threshold)
     - [exx\_cs\_inv\_thr](#exx_cs_inv_thr)
+    - [shrink\_abfs\_pca\_thr](#shrink_abfs_pca_thr)
+    - [shrink\_lu\_inv\_thr](#shrink_lu_inv_thr)
     - [exx\_v\_threshold](#exx_v_threshold)
     - [exx\_dm\_threshold](#exx_dm_threshold)
     - [exx\_c\_grad\_threshold](#exx_c_grad_threshold)
@@ -316,6 +318,10 @@
     - [rpa\_ccp\_rmesh\_times](#rpa_ccp_rmesh_times)
     - [exx\_symmetry\_realspace](#exx_symmetry_realspace)
     - [out\_ri\_cv](#out_ri_cv)
+    - [out\_unshrinked\_v](#out_unshrinked_v)
+    - [exx\_coul\_moment](#exx_coul_moment)
+    - [exx\_rotate\_abfs](#exx_rotate_abfs)
+    - [exx\_multip\_moments\_threshold](#exx_multip_moments_threshold)
   - [Exact Exchange (PW)](#exact-exchange-pw)
     - [exxace](#exxace)
     - [exx\_gamma\_extrapolation](#exx_gamma_extrapolation)
@@ -2945,6 +2951,18 @@
 - **Description**: By default, the Coulomb matrix inversion required for obtaining LRI coefficients is performed using LU decomposition. However, this approach may suffer from numerical instabilities when a large set of auxiliary basis functions (ABFs) is employed. When exx_cs_inv_thr > 0, the inversion is instead carried out via matrix diagonalization. Eigenvalues smaller than exx_cs_inv_thr are discarded to improve numerical stability. A relatively safe and commonly recommended value is 1e-5.
 - **Default**: -1
 
+### shrink_abfs_pca_thr
+
+- **Type**: Real
+- **Description**: Threshold to shrink the auxiliary basis for GW/RPA calculations.
+- **Default**: -1
+
+### shrink_lu_inv_thr
+
+- **Type**: Real
+- **Description**: Threshold for obtaining the inverse of the overlap matrix by LU decomposition in the auxiliary-basis representation.
+- **Default**: 1e-6
+
 ### exx_v_threshold
 
 - **Type**: Real
@@ -3042,6 +3060,30 @@
 - **Description**: Whether to output the coefficient tensor C(R) and ABFs-representation Coulomb matrix V(R) for each atom pair and cell in real space.
 - **Default**: false
 
+### out_unshrinked_v
+
+- **Type**: Boolean
+- **Description**: Whether to output the large Vq matrix in the unshrinked auxiliary basis.
+- **Default**: false
+
+### exx_coul_moment
+
+- **Type**: Boolean
+- **Description**: Whether to use the moment method for Coulomb calculation.
+- **Default**: false
+
+### exx_rotate_abfs
+
+- **Type**: Boolean
+- **Description**: Whether to rotate the auxiliary basis for Coulomb calculation.
+- **Default**: false
+
+### exx_multip_moments_threshold
+
+- **Type**: Real
+- **Description**: Threshold to screen multipole moments in Coulomb calculation.
+- **Default**: 1e-10
+
 [back to top](#full-list-of-input-keywords)
 
 ## Exact Exchange (PW)

source/Makefile.Objects

@@ -676,6 +676,8 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       symmetry_rotation.o\
       symmetry_rotation_output.o\
       symmetry_irreducible_sector.o\
+      gaussian_abfs.o\
+      singular_value.o\
 
 OBJS_PARALLEL=parallel_common.o\
     parallel_global.o\

source/source_cell/k_vector_utils.cpp

@@ -260,6 +260,7 @@ void kvec_mpi_k(K_Vectors& kv)
     std::vector<double> wk_aux(kv.nkstot);
     std::vector<double> kvec_c_aux(kv.nkstot * 3);
     std::vector<double> kvec_d_aux(kv.nkstot * 3);
+    std::vector<double> kvec_c_full_aux(kv.nkstot_full * 3);
 
     // collect and process in rank 0
     if (GlobalV::MY_RANK == 0)
@@ -274,6 +275,9 @@ void kvec_mpi_k(K_Vectors& kv)
             kvec_d_aux[3 * ik] = kv.kvec_d[ik].x;
             kvec_d_aux[3 * ik + 1] = kv.kvec_d[ik].y;
             kvec_d_aux[3 * ik + 2] = kv.kvec_d[ik].z;
+            kvec_c_full_aux[3 * ik] = kv.kvec_c_full[ik].x;
+            kvec_c_full_aux[3 * ik + 1] = kv.kvec_c_full[ik].y;
+            kvec_c_full_aux[3 * ik + 2] = kv.kvec_c_full[ik].z;
         }
     }
 
@@ -283,6 +287,7 @@ void kvec_mpi_k(K_Vectors& kv)
     Parallel_Common::bcast_double(wk_aux.data(), kv.nkstot);
     Parallel_Common::bcast_double(kvec_c_aux.data(), kv.nkstot * 3);
     Parallel_Common::bcast_double(kvec_d_aux.data(), kv.nkstot * 3);
+    Parallel_Common::bcast_double(kvec_c_full_aux.data(), kv.nkstot_full * 3);
 
     // process k point data in each processor
     kv.renew(kv.nks * kv.nspin);
@@ -300,6 +305,9 @@ void kvec_mpi_k(K_Vectors& kv)
         kv.kvec_d[i].x = kvec_d_aux[k_index * 3];
         kv.kvec_d[i].y = kvec_d_aux[k_index * 3 + 1];
         kv.kvec_d[i].z = kvec_d_aux[k_index * 3 + 2];
+        kv.kvec_c_full[i].x = kvec_c_full_aux[k_index * 3];
+        kv.kvec_c_full[i].y = kvec_c_full_aux[k_index * 3 + 1];
+        kv.kvec_c_full[i].z = kvec_c_full_aux[k_index * 3 + 2];
         kv.wk[i] = wk_aux[k_index];
         kv.isk[i] = isk_aux[k_index];
     }

source/source_cell/klist.cpp

@@ -87,6 +87,18 @@ void K_Vectors::set(const UnitCell& ucell,
     std::string skpt1;
     std::string skpt2;
 
+    if (!this->kc_done && this->kd_done)
+    {
+        for (size_t ik = 0; ik != this->nkstot_full; ++ik)
+            this->kvec_c_full[ik] = this->kvec_d[ik] * reciprocal_vec;
+    }
+    else if (this->kc_done && !this->kd_done)
+    {
+        for (size_t ik = 0; ik != this->nkstot_full; ++ik)
+            this->kvec_c_full[ik] = this->kvec_c[ik];
+    }
+
+
     // (2)
     // only berry phase need all kpoints including time-reversal symmetry!
     // if symm_flag is not set, only time-reversal symmetry would be considered.
@@ -182,6 +194,7 @@ void K_Vectors::renew(const int& kpoint_number)
 {
     kvec_c.resize(kpoint_number);
     kvec_d.resize(kpoint_number);
+    kvec_c_full.resize(kpoint_number);
     wk.resize(kpoint_number);
     isk.resize(kpoint_number);
     ngk.resize(kpoint_number);

source/source_cell/klist.h

@@ -14,6 +14,7 @@ class K_Vectors
 public:
     std::vector<ModuleBase::Vector3<double>> kvec_c; /// Cartesian coordinates of k points
     std::vector<ModuleBase::Vector3<double>> kvec_d; /// Direct coordinates of k points
+    std::vector<ModuleBase::Vector3<double>> kvec_c_full; // Cartesian coordinates of full k mesh match with nkstot_full
 
     std::vector<double> wk; /// wk, weight of k points
 

source/source_cell/test/klist_test_para.cpp

@@ -1,4 +1,5 @@
 #include "source_base/mathzone.h"
+#include "source_base/parallel_common.h"
 #include "source_base/parallel_global.h"
 #define private public
 #include "source_io/module_parameter/parameter.h"
@@ -237,6 +238,8 @@ TEST_F(KlistParaTest, Set)
     ModuleSymmetry::Symmetry::symm_flag = 1;
     kv->set(ucell,symm, k_file, kv->nspin, ucell.G, ucell.latvec,  GlobalV::ofs_running);
     EXPECT_EQ(kv->get_nkstot(), 35);
+    EXPECT_EQ(kv->get_nkstot_full(), 512);
+    EXPECT_GT(kv->get_nkstot_full(), kv->get_nkstot());
     EXPECT_TRUE(kv->kc_done);
     EXPECT_TRUE(kv->kd_done);
     if (GlobalV::NPROC == 4)
@@ -254,6 +257,64 @@ TEST_F(KlistParaTest, Set)
             EXPECT_EQ(kv->get_nks(), 17);
 }
     }
+    std::vector<double> local_kvec_c_full(kv->kvec_c_full.size() * 3);
+    for (size_t ik = 0; ik < kv->kvec_c_full.size(); ++ik)
+    {
+        local_kvec_c_full[3 * ik] = kv->kvec_c_full[ik].x;
+        local_kvec_c_full[3 * ik + 1] = kv->kvec_c_full[ik].y;
+        local_kvec_c_full[3 * ik + 2] = kv->kvec_c_full[ik].z;
+    }
+    const int local_count = static_cast<int>(local_kvec_c_full.size());
+    std::vector<int> counts;
+    std::vector<int> displs;
+    std::vector<int> pools;
+    if (GlobalV::MY_RANK == 0)
+    {
+        counts.resize(GlobalV::NPROC);
+        pools.resize(GlobalV::NPROC);
+    }
+    MPI_Gather(&local_count, 1, MPI_INT, counts.data(), 1, MPI_INT, 0, MPI_COMM_WORLD);
+    MPI_Gather(&GlobalV::MY_POOL, 1, MPI_INT, pools.data(), 1, MPI_INT, 0, MPI_COMM_WORLD);
+
+    std::vector<double> gathered_kvec_c_full;
+    if (GlobalV::MY_RANK == 0)
+    {
+        displs.resize(GlobalV::NPROC, 0);
+        for (int irank = 1; irank < GlobalV::NPROC; ++irank)
+        {
+            displs[irank] = displs[irank - 1] + counts[irank - 1];
+        }
+        gathered_kvec_c_full.resize(displs.back() + counts.back());
+    }
+    MPI_Gatherv(local_kvec_c_full.data(),
+                local_count,
+                MPI_DOUBLE,
+                gathered_kvec_c_full.data(),
+                counts.data(),
+                displs.data(),
+                MPI_DOUBLE,
+                0,
+                MPI_COMM_WORLD);
+    if (GlobalV::MY_RANK == 0)
+    {
+        for (int irank = 0; irank < GlobalV::NPROC; ++irank)
+        {
+            for (int jrank = irank + 1; jrank < GlobalV::NPROC; ++jrank)
+            {
+                if (pools[irank] != pools[jrank])
+                {
+                    continue;
+                }
+                ASSERT_EQ(counts[irank], counts[jrank]);
+                for (int i = 0; i < counts[irank]; ++i)
+                {
+                    EXPECT_NEAR(gathered_kvec_c_full[displs[irank] + i],
+                                gathered_kvec_c_full[displs[jrank] + i],
+                                1e-12);
+                }
+            }
+        }
+    }
     ClearUcell();
     if (GlobalV::MY_RANK == 0)
     {

source/source_hamilt/module_xc/exx_info.h

@@ -2,6 +2,7 @@
 #define EXX_INFO_H
 
 #include "source_lcao/module_ri/conv_coulomb_pot_k.h"
+#include "xc_functional.h"
 
 #include <vector>
 #include <map>
@@ -17,7 +18,7 @@ struct Exx_Info
 
 		// Fock:
 		//		"alpha":		"0"
-		//		"singularity_correction":	"limits" / "spencer" / "revised_spencer"
+		//		"singularity_correction":	"limits" / "spencer" / "revised_spencer" / "massidda" / "carrier"
 		//		"lambda":		"0.3"
         //      "Rcut"
 		// Erfc:
@@ -53,9 +54,12 @@ struct Exx_Info
         const std::map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> &coulomb_param;
 
         bool real_number = false;
+        bool coul_moment = false;
+        bool rotate_abfs = false;
 
         double pca_threshold = 0;
         std::vector<std::string> files_abfs;
+        std::vector<std::string> files_shrink_abfs;
         double C_threshold = 0;
         double V_threshold = 0;
         double dm_threshold = 0;
@@ -68,6 +72,13 @@ struct Exx_Info
         double kmesh_times = 4;
         double Cs_inv_thr = -1;
 
+        double shrink_abfs_pca_thr = -1;
+        double shrink_LU_inv_thr = 1e-6;
+        double multip_moments_threshold = 1e-10;
+        double exx_cs_inv_thr = -1;
+
+        int abfs_Lmax = 0; // tmp
+
         Exx_Info_RI(const Exx_Info::Exx_Info_Global& info_global)
             : coulomb_param(info_global.coulomb_param)
         {
@@ -94,12 +105,9 @@ struct Exx_Info
     }
 };
 
-//==========================================================
-// EXPLAIN : define "GLOBAL CLASS"
-//==========================================================
 namespace GlobalC
 {
     extern Exx_Info exx_info;
-} // namespace GlobalC
+}
 
 #endif

source/source_hsolver/hsolver_lcaopw.cpp

@@ -9,6 +9,7 @@
 #include "source_hsolver/diago_iter_assist.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_estate/elecstate_tools.h"
+#include "source_hamilt/module_xc/exx_info.h"
 
 
 #ifdef __EXX

source/source_io/module_ctrl/ctrl_scf_lcao.cpp

@@ -411,9 +411,7 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
     if (inp.rpa)
     {
         RPA_LRI<TK, double> rpa_lri_double(GlobalC::exx_info.info_ri);
-        rpa_lri_double.cal_postSCF_exx(*dm, MPI_COMM_WORLD, ucell, kv, orb);
-        rpa_lri_double.init(MPI_COMM_WORLD, kv, orb.cutoffs());
-        rpa_lri_double.out_for_RPA(ucell, pv, *psi, pelec);
+        rpa_lri_double.postSCF(ucell, MPI_COMM_WORLD, *dm, pelec, kv, orb, pv, *psi);
     }
 #endif
 

source/source_io/module_hs/single_R_io.cpp

@@ -6,7 +6,7 @@
 
 inline void write_data(std::ofstream& ofs, const double& data)
 {
-    ofs << " " << data;
+    ofs << " " << std::fixed << std::scientific << std::setprecision(16) << data;
 }
 inline void write_data(std::ofstream& ofs, const std::complex<double>& data)
 {
@@ -68,7 +68,6 @@ void ModuleIO::output_single_R(std::ofstream& ofs,
         if (!reduce || GlobalV::DRANK == 0)
         {
             long long nonzeros_count = 0;
-            ofs << std::fixed << std::scientific << std::setprecision(8);
             for (int col = 0; col < PARAM.globalv.nlocal; ++col)
             {
                 if (std::abs(line[col]) > sparse_threshold)