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_base/parallel_common.cpp

@@ -108,4 +108,35 @@ void bcast_char(char* object, const int n)
 #endif
 }
 
+#ifdef __MPI
+void bcast_int(int* object, const int n, const MPI_Comm mpi_comm, const int root)
+{
+    MPI_Bcast(object, n, MPI_INT, root, mpi_comm);
+}
+
+void bcast_double(double* object, const int n, const MPI_Comm mpi_comm, const int root)
+{
+    MPI_Bcast(object, n, MPI_DOUBLE, root, mpi_comm);
+}
+
+int my_rank(const MPI_Comm mpi_comm)
+{
+    int rank = 0;
+    MPI_Comm_rank(mpi_comm, &rank);
+    return rank;
+}
+
+int nproc(const MPI_Comm mpi_comm)
+{
+    int size = 1;
+    MPI_Comm_size(mpi_comm, &size);
+    return size;
+}
+
+void barrier(const MPI_Comm mpi_comm)
+{
+    MPI_Barrier(mpi_comm);
+}
+#endif
+
 } // namespace Parallel_Common

source/source_base/parallel_common.h

@@ -25,6 +25,19 @@ void bcast_double(double& object);
 void bcast_int(int& object);
 void bcast_bool(bool& object);
 
+#ifdef __MPI
+//(3) communicator- and root-aware counterparts for the sub-communicator
+//    stages (RI/RPA) that run outside MPI_COMM_WORLD and outside root 0
+void bcast_int(int* object, const int n, const MPI_Comm mpi_comm, const int root);
+void bcast_double(double* object, const int n, const MPI_Comm mpi_comm, const int root);
+/// @brief Rank of the caller on an arbitrary communicator.
+int my_rank(const MPI_Comm mpi_comm);
+/// @brief Number of ranks on an arbitrary communicator.
+int nproc(const MPI_Comm mpi_comm);
+/// @brief Barrier on an arbitrary communicator.
+void barrier(const MPI_Comm mpi_comm);
+#endif
+
 } // namespace Parallel_Common
 
 #endif

source/source_base/parallel_reduce.cpp

@@ -11,6 +11,7 @@ const MPI_Datatype Parallel_Reduce::MPI_Type<float>::value = MPI_FLOAT;
 const MPI_Datatype Parallel_Reduce::MPI_Type<std::complex<double>>::value = MPI_DOUBLE_COMPLEX;
 const MPI_Datatype Parallel_Reduce::MPI_Type<std::complex<float>>::value = MPI_C_FLOAT_COMPLEX;
 const MPI_Datatype Parallel_Reduce::MPI_Type<long long>::value = MPI_LONG_LONG;
+const MPI_Datatype Parallel_Reduce::MPI_Type<unsigned long long>::value = MPI_UNSIGNED_LONG_LONG;
 #endif
 
 template <typename T>
@@ -126,6 +127,67 @@ template void Parallel_Reduce::reduce_max<double>(double&);
 template void Parallel_Reduce::reduce_min_pool<double>(const int&, double&);
 template void Parallel_Reduce::reduce_max_pool<double>(const int&, double&);
 
+#ifdef __MPI
+template <typename T>
+void Parallel_Reduce::reduce_all(const MPI_Comm mpi_comm, T& object)
+{
+    MPI_Allreduce(MPI_IN_PLACE, &object, 1, MPI_Type<T>::value, MPI_SUM, mpi_comm);
+}
+
+template <typename T>
+void Parallel_Reduce::reduce_all(const MPI_Comm mpi_comm, T* object, const int n)
+{
+    MPI_Allreduce(MPI_IN_PLACE, object, n, MPI_Type<T>::value, MPI_SUM, mpi_comm);
+}
+
+template <typename T>
+void Parallel_Reduce::reduce_max(const MPI_Comm mpi_comm, T& object)
+{
+    MPI_Allreduce(MPI_IN_PLACE, &object, 1, MPI_Type<T>::value, MPI_MAX, mpi_comm);
+}
+
+template <typename T>
+void Parallel_Reduce::reduce_min(const MPI_Comm mpi_comm, T& object)
+{
+    MPI_Allreduce(MPI_IN_PLACE, &object, 1, MPI_Type<T>::value, MPI_MIN, mpi_comm);
+}
+
+void Parallel_Reduce::reduce_max_loc(const MPI_Comm mpi_comm, double& value, int& rank)
+{
+    struct MaxLoc
+    {
+        double value;
+        int rank;
+    } local{value, rank}, global{};
+    MPI_Allreduce(&local, &global, 1, MPI_DOUBLE_INT, MPI_MAXLOC, mpi_comm);
+    value = global.value;
+    rank = global.rank;
+}
+
+void Parallel_Reduce::gather_int_all(const MPI_Comm mpi_comm, int& v, int* all)
+{
+    assert(all != nullptr);
+    MPI_Allgather(&v, 1, MPI_INT, all, 1, MPI_INT, mpi_comm);
+}
+
+template void Parallel_Reduce::reduce_all<int>(const MPI_Comm, int&);
+template void Parallel_Reduce::reduce_all<double>(const MPI_Comm, double&);
+template void Parallel_Reduce::reduce_all<unsigned long long>(const MPI_Comm, unsigned long long&);
+
+template void Parallel_Reduce::reduce_all<int>(const MPI_Comm, int*, const int);
+template void Parallel_Reduce::reduce_all<double>(const MPI_Comm, double*, const int);
+template void Parallel_Reduce::reduce_all<unsigned long long>(const MPI_Comm, unsigned long long*, const int);
+
+template void Parallel_Reduce::reduce_max<int>(const MPI_Comm, int&);
+template void Parallel_Reduce::reduce_max<double>(const MPI_Comm, double&);
+template void Parallel_Reduce::reduce_max<unsigned long long>(const MPI_Comm, unsigned long long&);
+
+template void Parallel_Reduce::reduce_min<int>(const MPI_Comm, int&);
+template void Parallel_Reduce::reduce_min<double>(const MPI_Comm, double&);
+template void Parallel_Reduce::reduce_min<unsigned long long>(const MPI_Comm, unsigned long long&);
+#endif
+
+
 void Parallel_Reduce::reduce_double_grid(double* object, const int n)
 {
 #ifdef __MPI

source/source_base/parallel_reduce.h

@@ -50,6 +50,12 @@ struct MPI_Type<long long>
 {
     static const MPI_Datatype value;
 };
+
+template <>
+struct MPI_Type<unsigned long long>
+{
+    static const MPI_Datatype value;
+};
 #endif
 
 /// reduce in all process
@@ -82,6 +88,25 @@ void reduce_double_allpool(const int& npool, const int& nproc_in_pool, double* o
 
 void gather_int_all(int& v, int* all);
 
+#ifdef __MPI
+/// Communicator-aware counterparts of the MPI_COMM_WORLD-bound helpers above.
+/// The RI and RPA stages run on caller-provided sub-communicators, which those
+/// helpers cannot express; route every reduction on such a communicator
+/// through these overloads so the raw MPI calls stay centralized here.
+template <typename T>
+void reduce_all(const MPI_Comm mpi_comm, T& object);
+template <typename T>
+void reduce_all(const MPI_Comm mpi_comm, T* object, const int n);
+template <typename T>
+void reduce_max(const MPI_Comm mpi_comm, T& object);
+template <typename T>
+void reduce_min(const MPI_Comm mpi_comm, T& object);
+/// @brief Largest value together with the rank that supplied it (MPI_MAXLOC).
+void reduce_max_loc(const MPI_Comm mpi_comm, double& value, int& rank);
+/// @brief Allgather one int per rank on an arbitrary communicator.
+void gather_int_all(const MPI_Comm mpi_comm, int& v, int* all);
+#endif
+
 } // namespace Parallel_Reduce
 
 #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_esolver/esolver_lr_init.cpp

@@ -227,7 +227,7 @@ void ModuleESolver::ESolver_LR<T, TR>::initialize_from_ks_(UnitCell& ucell, cons
             exx_info.info_ri.files_abfs = ucell.abfs_orbital_files;
             exx_info.info_opt_abfs.files_abfs = ucell.abfs_orbital_files;
             exx_info.info_opt_abfs.files_jles = ucell.jle_orbital_files;
-            this->exx_lri = std::make_shared<Exx_LRI<T>>(exx_info.info_ri);
+            this->exx_lri = std::make_shared<Exx_LRI<T>>(exx_info.info_ri, make_exx_lri_runtime());
             this->exx_lri->init(MPI_COMM_WORLD, ucell, this->kv, ks_sol.orb_);
             this->exx_owned_ = true;   // the position-dependent `cal_exx_ions` is left to `refresh_from_ks_`
         }
@@ -487,7 +487,7 @@ void ModuleESolver::ESolver_LR<T, TR>::initialize_from_unitcell_(UnitCell& ucell
         exx_info.info_ri.files_abfs = ucell.abfs_orbital_files;
         exx_info.info_opt_abfs.files_abfs = ucell.abfs_orbital_files;
         exx_info.info_opt_abfs.files_jles = ucell.jle_orbital_files;
-        this->exx_lri = std::make_shared<Exx_LRI<T>>(exx_info.info_ri);
+        this->exx_lri = std::make_shared<Exx_LRI<T>>(exx_info.info_ri, make_exx_lri_runtime());
         this->exx_lri->init(MPI_COMM_WORLD, ucell,this->kv, orb);
         this->exx_lri->cal_exx_ions(ucell,this->inp_->out_ri_cv);
     }

source/source_hamilt/module_xc/exx_info.cpp

@@ -222,12 +222,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;