feat(librpa): strict 2D/3D Coulomb kernels and 2D reader-v1 producer outputs

source/source_base/element_basis_index.cpp

@@ -5,6 +5,9 @@
 
 #include "element_basis_index.h"
 
+#include <stdexcept>
+#include <vector>
+
 namespace ModuleBase
 {
 
@@ -37,4 +40,63 @@ Element_Basis_Index::construct_index( const Range &range )
     return index;
 }
 
-}
+Element_Basis_Index::IndexLNM
+Element_Basis_Index::construct_index( const Range &range, const IndexPermutation &old_to_new )
+{
+	if( old_to_new.empty() )
+	{
+		return construct_index(range);
+	}
+
+	IndexLNM index;
+	index.resize( range.size() );
+	for( std::size_t T=0; T!=range.size(); ++T )
+	{
+		std::size_t count_size = 0;
+		for( std::size_t L=0; L!=range[T].size(); ++L )
+		{
+			count_size += range[T][L].N * range[T][L].M;
+		}
+
+		const bool use_permutation = (T < old_to_new.size() && !old_to_new[T].empty());
+		if( use_permutation )
+		{
+			if( old_to_new[T].size() != count_size )
+			{
+				throw std::invalid_argument("Element_Basis_Index::construct_index permutation size mismatch.");
+			}
+
+			std::vector<bool> seen(count_size, false);
+			for( const std::size_t new_index : old_to_new[T] )
+			{
+				if( new_index >= count_size || seen[new_index] )
+				{
+					throw std::invalid_argument("Element_Basis_Index::construct_index permutation is not bijective.");
+				}
+				seen[new_index] = true;
+			}
+		}
+
+		std::size_t count = 0;
+		index[T].resize( range[T].size() );
+		for( std::size_t L=0; L!=range[T].size(); ++L )
+		{
+			index[T][L].resize( range[T][L].N );
+			for( std::size_t N=0; N!=range[T][L].N; ++N )
+			{
+				index[T][L][N].resize( range[T][L].M );
+				for( std::size_t M=0; M!=range[T][L].M; ++M )
+				{
+					index[T][L][N][M] = use_permutation ? old_to_new[T][count] : count;
+					++count;
+				}
+			}
+			index[T][L].N = range[T][L].N;
+			index[T][L].M = range[T][L].M;
+		}
+		index[T].count_size = count_size;
+	}
+	return index;
+}
+
+}

source/source_base/element_basis_index.h

@@ -16,34 +16,36 @@ namespace Element_Basis_Index
 {
   //private:
 
-    struct NM
-    {
-        public:
-        std::size_t N;
-        std::size_t M;
-    };
-
-    class Index_TL: public std::vector<std::vector<std::size_t>>
-    {
-        public:
-        std::size_t N;
-        std::size_t M;
-    };
-
-    class Index_T: public std::vector<Index_TL>
-    {
-        public:
-        std::size_t count_size;
-    };
+	struct NM
+	{
+		public:
+		std::size_t N;
+		std::size_t M;
+	};
+
+	class Index_TL: public std::vector<std::vector<std::size_t>>
+	{
+		public:
+		std::size_t N;
+		std::size_t M;
+	};
+
+	class Index_T: public std::vector<Index_TL>
+	{
+		public:
+		std::size_t count_size;
+	};
 
   //public:
 
-    typedef std::vector<std::vector<NM>> Range; 						// range[T][L]
-    typedef std::vector<Index_T> IndexLNM;								// index[T][L][N][M]
+	typedef std::vector<std::vector<NM>> Range; 						// range[T][L]
+	typedef std::vector<Index_T> IndexLNM;								// index[T][L][N][M]
+	typedef std::vector<std::vector<std::size_t>> IndexPermutation;		// old_to_new[T][old_index]
 
-    extern IndexLNM construct_index( const Range &range );
+	extern IndexLNM construct_index( const Range &range );
+	extern IndexLNM construct_index( const Range &range, const IndexPermutation &old_to_new );
 }
 
 }
 
-#endif
+#endif

source/source_cell/module_symmetry/irred_sec.h

@@ -11,6 +11,8 @@
 
 namespace ModuleSymmetry
 {
+    class SymmetryRotationTestAccess;
+
     using Tap = std::pair<int, int>;
     using TC = std::array<int, 3>;
     using TapR = std::pair<Tap, TC>;
@@ -149,6 +151,7 @@ namespace ModuleSymmetry
 
         friend class Symmetry_rotation;
         friend class Symmetry_rotation_k;
+        friend class SymmetryRotationTestAccess;
     };
 }
 

source/source_cell/module_symmetry/symm_rotation_k.cpp

@@ -235,8 +235,9 @@ namespace ModuleSymmetry
     // calculate Wigner D matrix
     double Symmetry_rotation_k::wigner_d(const double beta, const int l, const int m1, const int m2) const
     {
-        auto factorial = [](int n) -> int {
-            int result = 1;
+        // The four-factorial product overflows int for l >= 5, so evaluate in double.
+        auto factorial = [](int n) -> double {
+            double result = 1;
             for (int i = 1;i <= n;++i) { result *= i;
 }
             return result;

source/source_hamilt/module_xc/exx_info.cpp

@@ -221,12 +221,24 @@ bool init_exx_info(Exx_Info& exx_info, const Input_para& inp)
         exx_info.info_ri.pca_threshold = inp.exx_pca_threshold;
         exx_info.info_ri.C_threshold = inp.exx_c_threshold;
         exx_info.info_ri.V_threshold = inp.exx_v_threshold;
+        exx_info.info_ri.V_threshold_long = inp.exx_v_threshold_long;
+        exx_info.info_ri.V_cd_threshold = inp.exx_vcd_threshold;
+        exx_info.info_ri.V_cd_stats_only = inp.exx_vcd_stats_only;
+        exx_info.info_ri.V_cd_short_only = inp.exx_vcd_short_only;
         exx_info.info_ri.dm_threshold = inp.exx_dm_threshold;
         exx_info.info_ri.C_grad_threshold = inp.exx_c_grad_threshold;
         exx_info.info_ri.V_grad_threshold = inp.exx_v_grad_threshold;
         exx_info.info_ri.C_grad_R_threshold = inp.exx_c_grad_r_threshold;
         exx_info.info_ri.V_grad_R_threshold = inp.exx_v_grad_r_threshold;
         exx_info.info_ri.ccp_rmesh_times = std::stod(inp.exx_ccp_rmesh_times);
+        exx_info.info_ri.ewald_lambda = inp.exx_ewald_lambda;
+        exx_info.info_ri.ewald_tail_check = inp.exx_ewald_tail_check;
+        exx_info.info_ri.ewald_tail_abs_tol = inp.exx_ewald_tail_abs_tol;
+        exx_info.info_ri.ewald_tail_rel_tol = inp.exx_ewald_tail_rel_tol;
+        exx_info.info_ri.ewald_tail_sum_rel_tol = inp.exx_ewald_tail_sum_rel_tol;
+        exx_info.info_ri.ewald_tail_guard_cells = inp.exx_ewald_tail_guard_cells;
+        exx_info.info_ri.ewald_tail_max_expansions = inp.exx_ewald_tail_max_expansions;
+        exx_info.info_ri.ewald_dimension = inp.exx_ewald_dimension;
         exx_info.info_ri.exx_symmetry_realspace = inp.exx_symmetry_realspace;
         exx_info.info_ri.Cs_inv_thr = inp.exx_cs_inv_thr;
         exx_info.info_ri.shrink_abfs_pca_thr = inp.shrink_abfs_pca_thr;

source/source_hamilt/module_xc/exx_info_ri.h

@@ -20,12 +20,24 @@ struct Exx_Info_RI
     std::vector<std::string> files_shrink_abfs;
     double C_threshold = 0;
     double V_threshold = 0;
+        double V_threshold_long = 0;
+        double V_cd_threshold = -1;
+        bool V_cd_stats_only = false;
+        bool V_cd_short_only = true;
     double dm_threshold = 0;
     double C_grad_threshold = 0;
     double V_grad_threshold = 0;
     double C_grad_R_threshold = 0;
     double V_grad_R_threshold = 0;
     double ccp_rmesh_times = 10;
+        double ewald_lambda = 1.0;
+        std::string ewald_tail_check = "warn";
+        double ewald_tail_abs_tol = 1e-12;
+        double ewald_tail_rel_tol = 1e-10;
+        double ewald_tail_sum_rel_tol = 1e-8;
+        int ewald_tail_guard_cells = 1;
+        int ewald_tail_max_expansions = 3;
+        int ewald_dimension = 3;
     bool exx_symmetry_realspace = true;
     double kmesh_times = 4;
     double Cs_inv_thr = -1;

source/source_hamilt/module_xc/xc_functional.cpp

@@ -2,8 +2,12 @@
 #include "source_base/global_function.h"
 #include "source_base/tool_title.h"
 #include "source_base/constants.h"
+#include "source_io/module_parameter/parameter.h"
 #include <cmath>
 
+#include <algorithm>
+#include <cctype>
+
 #ifdef __LIBXC
 #include "libxc_abacus.h"
 #endif
@@ -77,6 +81,23 @@ void XC_Functional::set_runtime_parameters(const XCFunctionalParameters& paramet
     runtime_parameters = parameters;
 }
 
+std::string XC_Functional::resolve_runtime_xc_type(const std::string& input_xc_func,
+                                                   const std::string& pseudo_xc_func)
+{
+    std::string input_xc = input_xc_func;
+    std::transform(input_xc.begin(), input_xc.end(), input_xc.begin(), (::toupper));
+    if (!input_xc.empty() && input_xc != "DEFAULT" && input_xc != "NONE")
+    {
+        return input_xc_func;
+    }
+    return pseudo_xc_func;
+}
+
+std::string XC_Functional::resolve_runtime_xc_type(const UnitCell& ucell)
+{
+    return resolve_runtime_xc_type(PARAM.inp.dft_functional, ucell.atoms[0].ncpp.xc_func);
+}
+
 void XC_Functional::set_xc_first_loop(const UnitCell& ucell)
 {
     ModuleBase::TITLE("XC_Functional", "set_xc_first_loop");

source/source_hamilt/module_xc/xc_functional.h

@@ -103,6 +103,9 @@ class XC_Functional
     };
 
     static void set_xc_type(const std::string xc_func_in);
+    static std::string resolve_runtime_xc_type(const std::string& input_xc_func,
+                                               const std::string& pseudo_xc_func);
+    static std::string resolve_runtime_xc_type(const UnitCell& ucell);
 
     // Temporary adapter for legacy Potential -> PotXC -> XC evaluation globals.
     // Restores the complete selection on scope exit, including exception unwinding.

source/source_io/module_parameter/input_parameter.h

@@ -630,33 +630,54 @@ struct Input_para
     bool exx_separate_loop = true;                          ///< if 1, a two-step method is employed, else it will start
                                                             ///< with a GGA-Loop, and then Hybrid-Loop
     std::string exx_singularity_correction = "default";     ///< set the scheme of Coulomb singularity correction
+    int exx_ewald_dimension = 3;                             ///< dimensionality used by Ewald Coulomb singularity correction
     int exx_hybrid_step = 100;                              ///< the maximal electronic iteration number in
                                                             ///< the evaluation of Fock exchange
     double exx_mixing_beta = 1.0;                           ///< mixing_beta for outer-loop when exx_separate_loop=1
     std::string exx_real_number = "default";                ///< exx calculated in real or complex
     double exx_pca_threshold = 0.0001;                      ///< threshold to screen on-site ABFs in exx
     double exx_c_threshold = 0.0001;                        ///< threshold to screen C matrix in exx
     double exx_v_threshold = 0.1;                           ///< threshold to screen C matrix in exx
+    double exx_v_threshold_long = 0.0;                      ///< threshold to screen long-range V matrix in exx
+    double exx_vcd_threshold = -1.0;                        ///< dynamic C/D-weighted threshold for merged-short EXX paths
+    bool exx_vcd_stats_only = false;                        ///< collect weighted short-screen stats without skipping paths
+    bool exx_vcd_short_only = true;                         ///< apply weighted short-screening only to the split short channel
     double exx_dm_threshold = 0.0001;                       ///< threshold to screen density matrix in exx
     double exx_c_grad_threshold = 0.0001;                   ///< threshold to screen nabla C matrix in exx
     double exx_v_grad_threshold = 0.1;                      ///< threshold to screen nabla V matrix in exx
     double exx_c_grad_r_threshold = 0.0001;                 ///< threshold to screen nabla C * R matrix in exx
     double exx_v_grad_r_threshold = 0.1;                    ///< threshold to screen nabla V * R matrix in exx
     std::string exx_ccp_rmesh_times = "default";            ///< how many times larger the radial mesh required for
                                                             ///< calculating Columb potential is to that of atomic orbitals
+    double exx_ewald_lambda = 1.0;                          ///< Gaussian decay coefficient for Ewald full Coulomb in RI-EXX
+    std::string exx_ewald_tail_check = "warn";              ///< off, warn, enlarge, or strict short-range tail handling
+    double exx_ewald_tail_abs_tol = 1e-12;                  ///< absolute Frobenius-norm tolerance for one guard-shell block
+    double exx_ewald_tail_rel_tol = 1e-10;                  ///< relative Frobenius-norm tolerance for one guard-shell block
+    double exx_ewald_tail_sum_rel_tol = 1e-8;               ///< relative sum-of-block-norm tolerance for the guard shell
+    int exx_ewald_tail_guard_cells = 1;                     ///< lattice-cell layers added in every periodic direction
+    int exx_ewald_tail_max_expansions = 3;                  ///< maximum adaptive real-space range expansions
     int exx_opt_orb_lmax = 0;                               ///< the maximum l of the spherical Bessel functions for opt ABFs
     double exx_opt_orb_ecut = 0.0;                          ///< the cut-off of plane wave expansion for opt ABFs
     double exx_opt_orb_tolerence = 0.0;                     ///< the threshold when solving for the zeros of spherical Bessel
                                                             ///< functions for opt ABFs
     bool exx_symmetry_realspace = true; ///< whether to reduce the real-space sector in when using symmetry=1 in EXX calculation
     double rpa_ccp_rmesh_times = 10.0;  ///< how many times larger the radial mesh required for
                                         ///< calculating Columb potential is to that of atomic orbitals
+    std::string rpa_abfs_preorth = "none";  ///< on-site auxiliary Coulomb preorthogonalization for RPA output
+    double rpa_abfs_preorth_threshold = 1.0e-2;  ///< residual Coulomb-norm threshold for RPA ABFS preorthogonalization
     double exx_cs_inv_thr = -1;         ///< threshold to inverse Vq in abfs for generating Cs
     double shrink_abfs_pca_thr = -1;    ///< threshold to shrink auxiliary basis for GW/RPA
     double shrink_LU_inv_thr = 1e-6; ///< threshold to get inverse of overlap matrix by LU decomposition in auxiliary basis representation
     bool out_ri_cv = false;          ///< Whether to output the coefficient tensor C and ABFs-representation Coulomb matrix V
     bool out_unshrinked_v = false;   ///< whether to output the large Vq matrix in unshrinked auxiliary basis
     int out_librpa_ver = 0; ///< 0: legacy ABACUS LibRPA text files; 1: LibRPA reader-v1 binary files
+    std::string out_librpa_2d_coulomb_method = "ewald";  ///< ewald or direct_mixed_fourier for reader-v1 full V
+    double out_librpa_2d_direct_ecut = 0.0;  ///< reciprocal cutoff in Ry for direct mixed-Fourier 2D Coulomb
+    int out_librpa_2d_direct_kz_order = 0;   ///< positive kz quadrature order for direct mixed-Fourier 2D Coulomb
+    int out_librpa_2d_direct_gamma_order = 8;  ///< positive even Gamma-plane quadrature order
+    std::string out_librpa_3d_coulomb_method = "ewald";  ///< ewald or direct_reciprocal for reader-v1 full V
+    double out_librpa_3d_direct_ecut = 0.0;  ///< reciprocal cutoff in Ry for direct 3D Coulomb
+    bool out_librpa_abf_overlap = false;  ///< output raw active-ABF overlap for LibRPA v1 PSD diagnostics
     bool exx_coul_moment = false;    ///< whether to use moment method for Coulomb calculation
     bool exx_rotate_abfs = false;    ///< whether to rotate auxiliary basis for Coulomb calculation
     double exx_multip_moments_threshold = 1e-10; ///< threshold to screen multipole moments in Coulomb calculation