Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Refactor: make RealGauntTable singleton #2789

Merged
merged 3 commits into from
Aug 8, 2023
Merged

Refactor: make RealGauntTable singleton #2789

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3deb3b4
make RealGauntTable singleton & doxygen style change
jinzx10 Aug 3, 2023
e77521e
Merge branch 'develop' into singleton_gaunt
jinzx10 Aug 3, 2023
67209d6
Merge branch 'develop' into singleton_gaunt
wenfei-li Aug 8, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 14 additions & 11 deletions source/module_basis/module_nao/real_gaunt_table.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
#include "module_basis/module_nao/real_gaunt_table.h"

#include <algorithm>
#include <array>
#include <cassert>
#include <algorithm>

#include "module_base/constants.h"

Expand All @@ -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)
{
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
{
Expand All @@ -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);
Expand Down Expand Up @@ -162,19 +164,20 @@ 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<int>(m1 < 0) + static_cast<int>(m2 < 0) + static_cast<int>(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
{
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<int, 6> RealGauntTable::gaunt_key(const int l1, const int l2, const int l3, const int m1, const int m2, const int m3) const
Expand Down
143 changes: 82 additions & 61 deletions source/module_basis/module_nao/real_gaunt_table.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,151 +1,172 @@
#ifndef REAL_GAUNT_TABLE_H_
#define REAL_GAUNT_TABLE_H_

#include <map>
#include <array>
#include <map>

#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
* |
* Z(l,m) = | Y(l,0) m = 0
* |
* \ 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<std::array<int, 6>, 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<int, 6> 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
27 changes: 15 additions & 12 deletions source/module_basis/module_nao/test/real_gaunt_table_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
};

Expand All @@ -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);
Expand Down Expand Up @@ -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);
}
}
}
Expand All @@ -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)
Expand Down
5 changes: 2 additions & 3 deletions source/module_basis/module_nao/test/two_center_integrator_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
Loading