diff --git a/source/module_basis/module_nao/real_gaunt_table.cpp b/source/module_basis/module_nao/real_gaunt_table.cpp index 42c28319d89..a2567acb2e3 100644 --- a/source/module_basis/module_nao/real_gaunt_table.cpp +++ b/source/module_basis/module_nao/real_gaunt_table.cpp @@ -1,8 +1,8 @@ #include "module_basis/module_nao/real_gaunt_table.h" +#include #include #include -#include #include "module_base/constants.h" @@ -16,6 +16,10 @@ void RealGauntTable::build(const int lmax) return; } + // TODO + // If the table already exists and lmax is larger than the current lmax_, + // we should extend the table instead of rebuilding it from scratch. + // build the standard Gaunt table (with symmetry & selection rule considered) for (int l1 = 0; l1 <= 2 * lmax; ++l1) { @@ -67,7 +71,6 @@ void RealGauntTable::build(const int lmax) } } } - } const double& RealGauntTable::operator()(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const @@ -101,7 +104,6 @@ double RealGauntTable::real_gaunt_lookup(const int l1, const int l2, const int l else if ( m1 + m2 + m3 == 0 ) { return ModuleBase::SQRT2 / 2.0 * minus_1_pow(m_absmax + 1) * gaunt_lookup(l1, l2, l3, m1, m2, m3); - } else { @@ -123,10 +125,10 @@ double RealGauntTable::gaunt(const int l1, const int l2, const int l3, const int int g = (l1 + l2 + l3) / 2; double pref = std::sqrt( (2 * l1 + 1) * (2 * l2 + 1) * (2 * l3 + 1) / ModuleBase::FOUR_PI); double tri = std::sqrt( factorial(l1 + l2 - l3) * factorial(l2 + l3 - l1) * factorial(l3 + l1 - l2) - / factorial(l1+l2+l3+1) ); + / factorial(l1 + l2 + l3 + 1) ); // wigner3j(l1,l2,l3,0,0,0) - double wigner1 = minus_1_pow(g) * tri * factorial(g) / factorial(g-l1) / factorial(g-l2) / factorial(g-l3); + double wigner1 = minus_1_pow(g) * tri * factorial(g) / factorial(g - l1) / factorial(g - l2) / factorial(g - l3); // wigner3j(l1,l2,l3,m1,m2,m3) int kmin = std::max(l2 - l3 - m1, l1 - l3 + m2); @@ -162,10 +164,10 @@ bool RealGauntTable::gaunt_select_l(const int l1, const int l2, const int l3) co bool RealGauntTable::real_gaunt_select_m(const int m1, const int m2, const int m3) const { - return ( ( static_cast(m1 < 0) + static_cast(m2 < 0) + static_cast(m3 < 0) ) % 2 == 0 ) && - ( std::abs(m1) + std::abs(m2) == std::abs(m3) || - std::abs(m2) + std::abs(m3) == std::abs(m1) || - std::abs(m3) + std::abs(m1) == std::abs(m2) ); + return ( ( (m1 < 0) + (m2 < 0) + (m3 < 0) ) % 2 == 0 ) && + ( std::abs(m1) + std::abs(m2) == std::abs(m3) || + std::abs(m2) + std::abs(m3) == std::abs(m1) || + std::abs(m3) + std::abs(m1) == std::abs(m2) ); } double RealGauntTable::gaunt_lookup(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const @@ -173,8 +175,9 @@ double RealGauntTable::gaunt_lookup(const int l1, const int l2, const int l3, co assert( is_valid_lm(l1, l2, l3, m1, m2, m3) ); assert( l1 <= 2 * lmax_ && l2 <= 2 * lmax_ && l3 <= 2 * lmax_ ); - return ( gaunt_select_l(l1, l2, l3) && gaunt_select_m(m1, m2, m3) ) ? - gaunt_table_.at( gaunt_key(l1, l2, l3, m1, m2, m3) ) : 0.0; + return ( gaunt_select_l(l1, l2, l3) && gaunt_select_m(m1, m2, m3) ) + ? gaunt_table_.at( gaunt_key(l1, l2, l3, m1, m2, m3) ) + : 0.0; } std::array RealGauntTable::gaunt_key(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const diff --git a/source/module_basis/module_nao/real_gaunt_table.h b/source/module_basis/module_nao/real_gaunt_table.h index 557bcc11617..945af8aa8fc 100644 --- a/source/module_basis/module_nao/real_gaunt_table.h +++ b/source/module_basis/module_nao/real_gaunt_table.h @@ -1,33 +1,51 @@ #ifndef REAL_GAUNT_TABLE_H_ #define REAL_GAUNT_TABLE_H_ -#include #include +#include #include "module_base/module_container/tensor.h" -//! Table of Gaunt coefficients of real spherical harmonics /*! - * This class computes and stores the Gaunt coefficients of real spherical harmonics - * used in two-center integrals. + * @brief Table of Gaunt coefficients of real spherical harmonics. + * + * This class computes and tabulates the Gaunt coefficients of real spherical harmonics + * used in two-center integrals. The implementation follows the singleton pattern. + * + * Usage: + * + * int lmax = 5; + * RealGauntTable::instance().build(lmax); + * + * // get the real Gaunt coefficient + * double G = RealGauntTable::instance()(l1, l2, l3, m1, m2, m3); + * * */ class RealGauntTable { -public: + public: + RealGauntTable(RealGauntTable const&) = delete; + RealGauntTable& operator=(RealGauntTable const&) = delete; - RealGauntTable() {} ~RealGauntTable() {} - //! Builds the Gaunt table of real spherical harmonics + static RealGauntTable& instance() + { + static RealGauntTable instance_; + return instance_; + } + /*! - * This function tabulates the Gaunt coefficients of real spherical harmonics + * @brief Builds the Gaunt table of real spherical harmonics for two-center integrals. + * + * This function tabulates the Gaunt coefficients of real spherical harmonics * * / * G(l1,l2,l3,m1,m2,m3) = | Z(l1,m1) Z(l2,m2) Z(l3,m3) d Omega * / * - * for l1,l2 <= lmax and l3 <= 2*lmax. Here Z is the real spherical harmonics - * defined as + * + * for l1,l2 <= lmax and l3 <= 2*lmax. Here Z is the real spherical harmonics defined as * * / sqrt(2) * Re[Y(l,|m|)] m > 0 * | @@ -35,117 +53,120 @@ class RealGauntTable * | * \ sqrt(2) * Im[Y(l,|m|)] m < 0 * - * @note In some literature an extra pow(-1, m) is introduced to yield a signless - * Cartesian expression. The definition here is consistent with - * ModuleBase::Ylm::sph_harm and has some minus signs, for example, - * - * Z(1,-1) = -c * y / r - * Z(1, 0) = +c * z / r - * Z(1, 1) = -c * x / r + * @note In some literature an extra pow(-1, m) is introduced to yield a signless + * Cartesian expression. The definition here does not adopt this convention + * and is consistent with ModuleBase::Ylm::sph_harm. Consequently, Cartesian + * expression of Z may involve some minus signs, for example, * - * where c = sqrt(3/4/pi) and r = sqrt(x^2 + y^2 + z^2). + * Z(1,-1) = -c * y / r + * Z(1, 0) = +c * z / r + * Z(1, 1) = -c * x / r + * + * where c = sqrt(3/4/pi) and r = sqrt(x^2 + y^2 + z^2). * */ void build(const int lmax); - //! gets the tabulated real Gaunt coefficient + /// gets the tabulated real Gaunt coefficient const double& operator()(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const; - //! returns the maximum l + /// returns the maximum l (for the first two dimensions; the third dimension is 2*lmax) int lmax() const { return lmax_; } -private: + private: + RealGauntTable() {} - //! maximum angular momentum of the table (for the first two dimensions) + /// maximum angular momentum of the table (for the first two dimensions) int lmax_ = -1; - //! Table of standard Gaunt coefficients /*! - * This table maps (l1,l2,l3,m1,m2,m3) to a standard Gaunt coefficient. - * Due to the selection rule and symmetry, only those which survive the - * selection rule and satisfy l1 >= l2 >= l3 && m3 >= 0 are stored. + * @brief Table of standard Gaunt coefficients. + * + * This table maps (l1,l2,l3,m1,m2,m3) to a standard Gaunt coefficient. + * Due to the selection rule and symmetry, only those which survive the + * selection rule and satisfy l1 >= l2 >= l3 && m3 >= 0 are stored. * */ std::map, double> gaunt_table_; - //! Table of real Gaunt coefficients - /*! - * This table stores the real Gaunt coefficients. - * */ - container::Tensor real_gaunt_table_{ container::DataType::DT_DOUBLE, container::TensorShape({0}) }; + /// table of real Gaunt coefficients + container::Tensor real_gaunt_table_{container::DataType::DT_DOUBLE, container::TensorShape({0})}; - //! Gaunt coefficients /*! - * This function computes the standard Gaunt coefficients + * @brief Computes the standard Gaunt coefficients. + * + * This function computes the standard Gaunt coefficients * * / * G(l1,l2,l3,m1,m2,m3) = | Y(l1,m1) Y(l2,m2) Y(l3,m3) d Omega * / * - * where Y is the (standard) spherical harmonics and Omega is the solid angle element. + * where Y is the (standard) spherical harmonics and Omega is the solid angle element. * * - * @note This function computes the standard Gaunt coefficients, which is different - * from Gaunt coefficients of real spherical harmonics. - * @note Currently the algorithm computes the Gaunt coefficients with the Wigner-3j - * symbols, which in turn is evaluated with the Racah formula. This might have - * some numerical issue for large l and is yet to be studied later. + * @note This function computes the standard Gaunt coefficients, which is different + * from Gaunt coefficients of real spherical harmonics. + * @note Currently the algorithm computes the Gaunt coefficients with the Wigner-3j + * symbols, which in turn is evaluated with the Racah formula. This might have + * some numerical issue for large l and is yet to be studied later. * */ double gaunt(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const; - //! selection rule of standard & real Gaunt coefficients regarding l1, l2, l3 + /// selection rule of standard & real Gaunt coefficients regarding l1, l2, l3 bool gaunt_select_l(const int l1, const int l2, const int l3) const; - //! selection rule of standard Gaunt coefficients regarding m1, m2, m3 + /// selection rule of standard Gaunt coefficients regarding m1, m2, m3 bool gaunt_select_m(const int m1, const int m2, const int m3) const { return m1 + m2 + m3 == 0; } - //! selection rule of real Gaunt coefficients regarding m1, m2, m3 + /// selection rule of real Gaunt coefficients regarding m1, m2, m3 bool real_gaunt_select_m(const int m1, const int m2, const int m3) const; - //! returns whether the given l & m are valid quantum numbers /*! - * This function checks whether abs(mi) <= li (i=1,2,3) is satisfied. - * This implies li >= 0. + * @brief Returns whether the given l & m are valid quantum numbers. + * + * This function checks whether abs(mi) <= li (i=1,2,3) is satisfied. + * This implies li >= 0. * */ bool is_valid_lm(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const; - //! Get a Gaunt coefficient by looking up the table + /// Get a Gaunt coefficient by looking up the table double gaunt_lookup(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const; - //! Get a real Gaunt coefficient from the stored Gaunt coefficients + /// Get a real Gaunt coefficient from the stored Gaunt coefficients double real_gaunt_lookup(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const; - //! Symmetry-adapted key for gaunt_table_ /*! - * Standard Gaunt coefficients have the following symmetries: + * @brief Symmetry-adapted key for gaunt_table_. + * + * Standard Gaunt coefficients have the following symmetries: * * Gaunt(l1,l2,l3,m1,m2,m3) = Gaunt(l1,l2,l3,-m1,-m2,-m3) * * Gaunt(1,2,3) = Gaunt(2,3,1) = Gaunt(3,1,2) = * Gaunt(2,1,3) = Gaunt(1,3,2) = Gaunt(3,2,1) * - * The above symmetries enable us to store merely a small portion of the Gaunt - * coefficients. This function permutes 1/2/3 and flips the signs of m1/m2/m3 - * if necessary so that the returned key {l1,l2,l3,m1,m2,m3} satisfies - * l1 >= l2 >= l3 and m3 >= 0. + * These symmetries enable us to store merely a small portion of the Gaunt + * coefficients. This function permutes 1/2/3 and flips the signs of m1/m2/m3 + * if necessary so that the returned key {l1,l2,l3,m1,m2,m3} satisfies + * l1 >= l2 >= l3 and m3 >= 0. * */ std::array gaunt_key(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const; - //! swap (l1,m1) <--> (l2,m2) if l1 < l2; do nothing otherwise + /// swap (l1,m1) <--> (l2,m2) if l1 < l2; do nothing otherwise void arrange(int& l1, int& l2, int& m1, int& m2) const; - //! returns n! as a double + /// returns n! as a double double factorial(const int n) const; - //! returns the linearized index of Y(l,m) /*! - * l 0 1 1 1 2 2 2 2 2 3 ... - * m 0 -1 0 1 -2 -1 0 1 2 -3 ... - * index 0 1 2 3 4 5 6 7 8 9 ... + * @brief Returns the linearized index of Y(l,m). + * + * l 0 1 1 1 2 2 2 2 2 3 ... + * m 0 -1 0 1 -2 -1 0 1 2 -3 ... + * index 0 1 2 3 4 5 6 7 8 9 ... * */ int index_map(int l, int m) const; - //! returns pow(-1, m) + /// returns pow(-1, m) int minus_1_pow(int m) const { return m % 2 ? -1 : 1; } - }; #endif diff --git a/source/module_basis/module_nao/test/real_gaunt_table_test.cpp b/source/module_basis/module_nao/test/real_gaunt_table_test.cpp index 8726bcf98f8..e38e1008185 100644 --- a/source/module_basis/module_nao/test/real_gaunt_table_test.cpp +++ b/source/module_basis/module_nao/test/real_gaunt_table_test.cpp @@ -30,11 +30,11 @@ using iclock = std::chrono::high_resolution_clock; class RealGauntTableTest : public ::testing::Test { protected: - void SetUp() { rgt.build(lmax); } + void SetUp() { /*rgt.build(lmax);*/ } void TearDown() {} int lmax = 10; //!< maximum angular momentum - RealGauntTable rgt; //!< object under test + //RealGauntTable rgt; //!< object under test const double tol = 1e-12; //!< numerical error tolerance for individual Gaunt coefficient }; @@ -47,6 +47,8 @@ TEST_F(RealGauntTableTest, LegacyConsistency) // this test shall be removed in the future once the refactoring is finished ORB_gaunt_table ogt; + RealGauntTable::instance().build(lmax); + //start = iclock::now(); ogt.init_Gaunt_CH(lmax); ogt.init_Gaunt(lmax); @@ -80,7 +82,8 @@ TEST_F(RealGauntTableTest, LegacyConsistency) int m2 = ogt.Index_M(mm2); int m3 = ogt.Index_M(mm3); - EXPECT_NEAR(rgt(l1, l2, l3, m1, m2, m3), ogt.Gaunt_Coefficients(index1, index2, index3), tol); + EXPECT_NEAR(RealGauntTable::instance()(l1, l2, l3, m1, m2, m3), + ogt.Gaunt_Coefficients(index1, index2, index3), tol); } } } @@ -91,20 +94,20 @@ TEST_F(RealGauntTableTest, LegacyConsistency) TEST_F(RealGauntTableTest, SanityCheck) { - EXPECT_EQ(rgt.lmax(), lmax); + EXPECT_EQ(RealGauntTable::instance().lmax(), lmax); - EXPECT_NEAR(rgt(0, 0, 0, 0, 0, 0), ModuleBase::SQRT_INVERSE_FOUR_PI, tol); + EXPECT_NEAR(RealGauntTable::instance()(0, 0, 0, 0, 0, 0), ModuleBase::SQRT_INVERSE_FOUR_PI, tol); - EXPECT_NEAR(rgt(4, 0, 4, 3, 0, 3), ModuleBase::SQRT_INVERSE_FOUR_PI, tol); - EXPECT_NEAR(rgt(4, 0, 4, -3, 0, -3), ModuleBase::SQRT_INVERSE_FOUR_PI, tol); + EXPECT_NEAR(RealGauntTable::instance()(4, 0, 4, 3, 0, 3), ModuleBase::SQRT_INVERSE_FOUR_PI, tol); + EXPECT_NEAR(RealGauntTable::instance()(4, 0, 4, -3, 0, -3), ModuleBase::SQRT_INVERSE_FOUR_PI, tol); - EXPECT_NEAR(rgt(2, 2, 2, 2, -1, -1), -std::sqrt(15.0) / 7.0 * ModuleBase::SQRT_INVERSE_FOUR_PI, tol); - EXPECT_NEAR(rgt(2, 2, 2, -1, 2, -1), -std::sqrt(15.0) / 7.0 * ModuleBase::SQRT_INVERSE_FOUR_PI, tol); + EXPECT_NEAR(RealGauntTable::instance()(2, 2, 2, 2, -1, -1), -std::sqrt(15.0) / 7.0 * ModuleBase::SQRT_INVERSE_FOUR_PI, tol); + EXPECT_NEAR(RealGauntTable::instance()(2, 2, 2, -1, 2, -1), -std::sqrt(15.0) / 7.0 * ModuleBase::SQRT_INVERSE_FOUR_PI, tol); - EXPECT_NEAR(rgt(3, 3, 2, 2, 1, 1), ModuleBase::SQRT_INVERSE_FOUR_PI / std::sqrt(6.0), tol); - EXPECT_NEAR(rgt(2, 3, 3, 1, 1, 2), ModuleBase::SQRT_INVERSE_FOUR_PI / std::sqrt(6.0), tol); + EXPECT_NEAR(RealGauntTable::instance()(3, 3, 2, 2, 1, 1), ModuleBase::SQRT_INVERSE_FOUR_PI / std::sqrt(6.0), tol); + EXPECT_NEAR(RealGauntTable::instance()(2, 3, 3, 1, 1, 2), ModuleBase::SQRT_INVERSE_FOUR_PI / std::sqrt(6.0), tol); - EXPECT_NEAR(rgt(4, 5, 7, 3, -2, -5), ModuleBase::SQRT_INVERSE_FOUR_PI * std::sqrt(210.0) / 221.0, tol); + EXPECT_NEAR(RealGauntTable::instance()(4, 5, 7, 3, -2, -5), ModuleBase::SQRT_INVERSE_FOUR_PI * std::sqrt(210.0) / 221.0, tol); } int main(int argc, char** argv) diff --git a/source/module_basis/module_nao/test/two_center_integrator_test.cpp b/source/module_basis/module_nao/test/two_center_integrator_test.cpp index 4b218f3cad2..c1945171351 100644 --- a/source/module_basis/module_nao/test/two_center_integrator_test.cpp +++ b/source/module_basis/module_nao/test/two_center_integrator_test.cpp @@ -89,9 +89,8 @@ TEST_F(TwoCenterIntegratorTest, FiniteDifference) start = iclock::now(); - RealGauntTable rgt; - S_intor.tabulate(orb, orb, 'S', nr, rmax, true, &rgt); - T_intor.tabulate(orb, orb, 'T', nr, rmax, true, &rgt); + S_intor.tabulate(orb, orb, 'S', nr, rmax, true); + T_intor.tabulate(orb, orb, 'T', nr, rmax, true); dur = iclock::now() - start; std::cout << "time elapsed = " << dur.count() << " s" << std::endl; diff --git a/source/module_basis/module_nao/two_center_integrator.cpp b/source/module_basis/module_nao/two_center_integrator.cpp index bef65638919..00117dacbca 100644 --- a/source/module_basis/module_nao/two_center_integrator.cpp +++ b/source/module_basis/module_nao/two_center_integrator.cpp @@ -3,54 +3,24 @@ #include "module_base/vector3.h" #include "module_base/ylm.h" -TwoCenterIntegrator::TwoCenterIntegrator() : +TwoCenterIntegrator::TwoCenterIntegrator(): is_tabulated_(false), op_('\0'), - with_deriv_(false), - use_internal_gaunt_(false), - rgt_(nullptr) + with_deriv_(false) { } -TwoCenterIntegrator::~TwoCenterIntegrator() -{ - if (use_internal_gaunt_) - { - delete rgt_; - } -} - void TwoCenterIntegrator::tabulate(const RadialCollection& bra, const RadialCollection& ket, const char op, const int nr, const double cutoff, - const bool with_deriv, - RealGauntTable* const rgt) + const bool with_deriv) { op_ = op; with_deriv_ = with_deriv; table_.build(bra, ket, op, nr, cutoff, with_deriv); - - if (rgt) - { // if an external gaunt table is provided - if (use_internal_gaunt_) - { - delete rgt_; - use_internal_gaunt_ = false; - } - rgt_ = rgt; - } - else - { // if no external gaunt table is provided (which implies an internal one) - if (!use_internal_gaunt_) - { - rgt_ = new RealGauntTable; - use_internal_gaunt_ = true; - } - } - - rgt_->build(std::max(bra.lmax(), ket.lmax())); + RealGauntTable::instance().build(std::max(bra.lmax(), ket.lmax())); is_tabulated_ = true; } @@ -102,7 +72,7 @@ void TwoCenterIntegrator::calculate(const int itype1, for (int m = -l; m < l; ++m) { - double G = (*rgt_)(l1, l2, l, m1, m2, m); + double G = RealGauntTable::instance()(l1, l2, l, m1, m2, m); if (deriv) { diff --git a/source/module_basis/module_nao/two_center_integrator.h b/source/module_basis/module_nao/two_center_integrator.h index 0dbbf456f23..c2b95a29916 100644 --- a/source/module_basis/module_nao/two_center_integrator.h +++ b/source/module_basis/module_nao/two_center_integrator.h @@ -6,7 +6,7 @@ #include "module_basis/module_nao/radial_collection.h" #include "module_base/vector3.h" -/** +/*! * @brief A class to compute two-center integrals * * This class computes two-center integrals @@ -38,10 +38,10 @@ class TwoCenterIntegrator TwoCenterIntegrator(const TwoCenterIntegrator&) = delete; TwoCenterIntegrator& operator=(const TwoCenterIntegrator&) = delete; - ~TwoCenterIntegrator(); + ~TwoCenterIntegrator() {} - /** - * @brief Tabulate the radial part of the two-center integrals + /*! + * @brief Tabulates the radial part of a two-center integral. * * @param[in] bra The radial functions of the first collection. * @param[in] ket The radial functions of the second collection. @@ -50,19 +50,17 @@ class TwoCenterIntegrator * @param[in] cutoff r-space cutoff radius. * @param[in] with_deriv If true, the derivative of radial table is also tabulated. * This is necessary to compute the gradient of integrals. - * @param[in] rgt Pointer to a real Gaunt table. * */ void tabulate(const RadialCollection& bra, const RadialCollection& ket, const char op, const int nr, const double cutoff, - const bool with_deriv, - RealGauntTable* const rgt = nullptr + const bool with_deriv ); - /** - * @brief Compute the two-center integrals + /*! + * @brief Compute the two-center integrals. * * This function calculates the two-center integral * @@ -103,13 +101,11 @@ class TwoCenterIntegrator char op_; bool with_deriv_; - bool use_internal_gaunt_; TwoCenterTable table_; - RealGauntTable* rgt_; - /** - * @brief Compute the index of (l,m) in the array of spherical harmonics + /*! + * @brief Returns the index of (l,m) in the array of spherical harmonics. * * Spherical harmonics in ABACUS are stored in the following order: *