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: Two center integral implemented in module_nao #2757

Merged
merged 9 commits into from
Jul 31, 2023
Merged

Refactor: Two center integral implemented in module_nao #2757

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
9543d51
add lookup function to TwoCenterTable
jinzx10 Jul 18, 2023
ff5a228
Merge branch 'deepmodeling:develop' into two_center_integral
jinzx10 Jul 21, 2023
f895486
rabtab assumes uniform grid; table behavior change
jinzx10 Jul 24, 2023
7845b6c
two_center_integral initial commit
jinzx10 Jul 24, 2023
21023a2
Merge branch 'deepmodeling:develop' into two_center_integral
jinzx10 Jul 25, 2023
cf8979d
TwoCenterIntegrator implemented with test
jinzx10 Jul 26, 2023
cbdcb52
Merge branch 'develop' into two_center_integral
jinzx10 Jul 26, 2023
708e02e
Merge branch 'develop' into two_center_integral
jinzx10 Jul 27, 2023
b5d4252
Merge branch 'develop' into two_center_integral
jinzx10 Jul 30, 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
9 changes: 7 additions & 2 deletions source/module_basis/module_nao/numerical_radial.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -365,17 +365,21 @@ void NumericalRadial::radtab(const char op,
const int l,
double* const table,
const int nr_tab,
const double* const rgrid_tab,
const double rmax_tab,
const bool deriv) const
{
assert(op == 'S' || op == 'I' || op == 'T' || op == 'U');
assert(l >= 0);
assert(rgrid_tab && nr_tab > 0);
assert(rmax_tab > 0 && nr_tab > 0);

// radtab requires that two NumericalRadial objects have exactly the same (non-null) kgrid_
assert(nk_ > 0 && nk_ == ket.nk_);
assert(std::equal(kgrid_, kgrid_ + nk_, ket.kgrid_));

double* rgrid_tab = new double[nr_tab];
double dr = rmax_tab / (nr_tab - 1);
std::for_each(rgrid_tab, rgrid_tab + nr_tab, [dr,&rgrid_tab](double& r) { r = dr * (&r - rgrid_tab); });

bool use_radrfft = is_fft_compliant(nr_tab, rgrid_tab, nk_, kgrid_);

// function to undergo a spherical Bessel transform:
Expand Down Expand Up @@ -407,6 +411,7 @@ void NumericalRadial::radtab(const char op,
}

delete[] fk;
delete[] rgrid_tab;

// spherical Bessel transform has a prefactor of sqrt(2/pi) while the prefactor for the radial table
// of two-center integrals is 4*pi
Expand Down
22 changes: 11 additions & 11 deletions source/module_basis/module_nao/numerical_radial.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,17 +191,17 @@ class NumericalRadial
* / 0 l
*
* */
void radtab(const char op, //!< [in] operator, could be:
//!< - 'S' or 'I': overlap
//!< - 'T': kinetic
//!< - 'U': Coulomb
const NumericalRadial& ket, //!< [in] the other NumericalRadial object with which
//! the two-center integral is computed
const int l, //!< [in] angular momentum of the table
double* const table, //!< [out] on finish, contain the computed table
const int nr_tab, //!< [in] size of table grid
const double* const rgrid_tab, //!< [in] grid on which the table is calculated.
const bool deriv = false //!< [in] if true, calculates the derivative of the table
void radtab(const char op, //!< [in] operator, could be:
//!< - 'S' or 'I': overlap
//!< - 'T': kinetic
//!< - 'U': Coulomb
const NumericalRadial& ket, //!< [in] the other NumericalRadial object with which
//! the two-center integral is computed
const int l, //!< [in] angular momentum of the table
double* const table, //!< [out] on finish, contain the computed table
const int nr_tab, //!< [in] size of table grid
const double rmax_tab, //!< [in] cutoff radius of table grid
const bool deriv = false //!< [in] if true, calculates the derivative of the table
) const;

//! Normalizes the radial function.
Expand Down
14 changes: 14 additions & 0 deletions source/module_basis/module_nao/test/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,3 +60,17 @@ AddTest(
LIBS ${math_libs} device base container kernels
)

AddTest(
TARGET two_center_integrator
SOURCES
two_center_integrator_test.cpp
../two_center_integrator.cpp
../two_center_table.cpp
../real_gaunt_table.cpp
../radial_collection.cpp
../atomic_radials.cpp
../beta_radials.cpp
../radial_set.cpp
../numerical_radial.cpp
LIBS ${math_libs} device base container kernels orb
)
17 changes: 9 additions & 8 deletions source/module_basis/module_nao/test/numerical_radial_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -472,58 +472,59 @@ TEST_F(NumericalRadialTest, RadialTable)

double* table = new double[sz];
double table_pref = ModuleBase::FOUR_PI * std::sqrt(ModuleBase::PI / 2.0);
double rmax_tab = chi1.rcut();

chi1.radtab('S', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid());
chi1.radtab('S', chi2, 0, table, chi1.nr(), rmax_tab);
for (int i = 0; i != sz; ++i)
{
double R = i * dr;
EXPECT_NEAR(table[i], table_pref * (3 - R * R) / 32 * std::exp(-R * R / 2), tol);
}

chi1.radtab('S', chi2, 2, table, chi1.nr(), chi1.ptr_rgrid());
chi1.radtab('S', chi2, 2, table, chi1.nr(), rmax_tab);
for (int i = 0; i != sz; ++i)
{
double R = i * dr;
EXPECT_NEAR(table[i], table_pref * R * R / 32 * std::exp(-R * R / 2), tol);
}

chi1.radtab('T', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid());
chi1.radtab('T', chi2, 0, table, chi1.nr(), rmax_tab);
for (int i = 0; i != sz; ++i)
{
double R = i * dr;
EXPECT_NEAR(table[i], table_pref * (std::pow(R, 4) - 10 * R * R + 15) / 32 * std::exp(-R * R / 2), tol);
}

chi1.radtab('U', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid());
chi1.radtab('U', chi2, 0, table, chi1.nr(), rmax_tab);
for (int i = 0; i != sz; ++i)
{
double R = i * dr;
EXPECT_NEAR(table[i], table_pref * 1. / 32 * std::exp(-R * R / 2), tol);
}

// derivative of radial tables
chi1.radtab('S', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), true);
chi1.radtab('S', chi2, 0, table, chi1.nr(), rmax_tab, true);
for (int i = 0; i != sz; ++i)
{
double R = i * dr;
EXPECT_NEAR(table[i], table_pref * (std::pow(R, 3) - 5 * R) / 32 * std::exp(-R * R / 2), tol);
}

chi1.radtab('S', chi2, 2, table, chi1.nr(), chi1.ptr_rgrid(), true);
chi1.radtab('S', chi2, 2, table, chi1.nr(), rmax_tab, true);
for (int i = 0; i != sz; ++i)
{
double R = i * dr;
EXPECT_NEAR(table[i], table_pref * (2 * R - std::pow(R, 3)) / 32 * std::exp(-R * R / 2), tol);
}

chi1.radtab('T', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), true);
chi1.radtab('T', chi2, 0, table, chi1.nr(), rmax_tab, true);
for (int i = 0; i != sz; ++i)
{
double R = i * dr;
EXPECT_NEAR(table[i], table_pref * (-std::pow(R, 5) + 14 * std::pow(R, 3) - 35 * R) / 32 * std::exp(-R * R / 2), tol);
}

chi1.radtab('U', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), true);
chi1.radtab('U', chi2, 0, table, chi1.nr(), rmax_tab, true);
for (int i = 0; i != sz; ++i)
{
double R = i * dr;
Expand Down
Loading