From 954c3583c0e3f038ddcf4236b3be01d1e73e70a2 Mon Sep 17 00:00:00 2001 From: Taoni Bao Date: Thu, 26 Sep 2024 12:00:24 +0800 Subject: [PATCH 1/3] Fix: `timer` not closed in TDDFT code (#5172) * Fix a simple timer bug in TDDFT * Fix timer TD_Velocity::initialize_current_term * Fix timer ModuleIO::write_current --- .../module_tddft/evolve_elec.cpp | 6 +-- .../module_tddft/td_velocity.cpp | 9 ++-- source/module_io/td_current_io.cpp | 41 +++++++++++-------- 3 files changed, 31 insertions(+), 25 deletions(-) diff --git a/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp b/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp index 4ecca906cf..7be680cfff 100644 --- a/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp +++ b/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp @@ -36,8 +36,8 @@ void Evolve_elec::solve_psi(const int& istep, int propagator, const int& nks) { - ModuleBase::TITLE("Evolve_elec", "eveolve_psi"); - ModuleBase::timer::tick("Evolve_elec", "evolve_psi"); + ModuleBase::TITLE("Evolve_elec", "solve_psi"); + ModuleBase::timer::tick("Evolve_elec", "solve_psi"); for (int ik = 0; ik < nks; ik++) { @@ -82,7 +82,7 @@ void Evolve_elec::solve_psi(const int& istep, ModuleBase::timer::tick("Efficience", "evolve_k"); } // end k - ModuleBase::timer::tick("Evolve_elec", "evolve_psi"); + ModuleBase::timer::tick("Evolve_elec", "solve_psi"); return; } } // namespace module_tddft \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/td_velocity.cpp b/source/module_hamilt_lcao/module_tddft/td_velocity.cpp index 42a8865052..72267a31a1 100644 --- a/source/module_hamilt_lcao/module_tddft/td_velocity.cpp +++ b/source/module_hamilt_lcao/module_tddft/td_velocity.cpp @@ -1,7 +1,7 @@ #include "td_velocity.h" -#include "module_parameter/parameter.h" #include "module_elecstate/potentials/H_TDDFT_pw.h" +#include "module_parameter/parameter.h" bool TD_Velocity::tddft_velocity = false; bool TD_Velocity::out_mat_R = false; @@ -61,7 +61,7 @@ void TD_Velocity::output_cart_At(const std::string& out_dir) // divide by 2.0 to get the atomic unit for (int i = 0; i < 3; i++) { - ofs << std::scientific << std::setprecision(4) << std::setw(15) << cart_At[i] ; + ofs << std::scientific << std::setprecision(4) << std::setw(15) << cart_At[i]; } ofs << std::endl; ofs.close(); @@ -79,7 +79,7 @@ void TD_Velocity::cal_cart_At(const ModuleBase::Vector3& At) else { // transfrom into atomic unit - this->cart_At = At/2.0; + this->cart_At = At / 2.0; } // output the vector potential if needed if (out_vecpot == true) @@ -171,7 +171,8 @@ void TD_Velocity::initialize_current_term(const hamilt::HContainercurrent_term[dir]->allocate(nullptr, true); } - ModuleBase::timer::tick("TDEkinetic", "initialize_HR_tmp"); + + ModuleBase::timer::tick("TD_Velocity", "initialize_current_term"); } void TD_Velocity::destroy_HS_R_td_sparse(void) diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp index bb13b3a9a1..72cd9fd412 100644 --- a/source/module_io/td_current_io.cpp +++ b/source/module_io/td_current_io.cpp @@ -1,6 +1,5 @@ #include "td_current_io.h" -#include "module_parameter/parameter.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_base/libm/libm.h" @@ -11,9 +10,10 @@ #include "module_elecstate/module_dm/cal_dm_psi.h" #include "module_elecstate/potentials/H_TDDFT_pw.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" -#include "module_hamilt_lcao/module_tddft/td_velocity.h" #include "module_hamilt_lcao/module_tddft/td_current.h" +#include "module_hamilt_lcao/module_tddft/td_velocity.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_parameter/parameter.h" #ifdef __LCAO @@ -138,7 +138,7 @@ void ModuleIO::write_current(const int istep, cal_current->calculate_grad_term(); for (int dir = 0; dir < 3; dir++) { - current_term[dir]=cal_current->get_current_term_pointer(dir); + current_term[dir] = cal_current->get_current_term_pointer(dir); } } else @@ -149,11 +149,10 @@ void ModuleIO::write_current(const int istep, } for (int dir = 0; dir < 3; dir++) { - current_term[dir]=TD_Velocity::td_vel_op->get_current_term_pointer(dir); + current_term[dir] = TD_Velocity::td_vel_op->get_current_term_pointer(dir); } } - // construct a DensityMatrix object // Since the function cal_dm_psi do not suport DMR in complex type, I replace it with two DMR in double type. Should // be refactored in the future. @@ -224,12 +223,17 @@ void ModuleIO::write_current(const int istep, // std::cout<< "iat1: " << iat1 << " iat2: " << iat2 << " Rx: " << Rx << " Ry: " << Ry << " Rz: // " << Rz << std::endl; // get BaseMatrix - hamilt::BaseMatrix* tmp_matrix_real = DM_real.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); - hamilt::BaseMatrix* tmp_matrix_imag = DM_imag.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix* tmp_matrix_real + = DM_real.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix* tmp_matrix_imag + = DM_imag.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); // refactor - hamilt::BaseMatrix>* tmp_m_rvx = current_term[0]->find_matrix(iat1, iat2, Rx, Ry, Rz); - hamilt::BaseMatrix>* tmp_m_rvy = current_term[1]->find_matrix(iat1, iat2, Rx, Ry, Rz); - hamilt::BaseMatrix>* tmp_m_rvz = current_term[2]->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix>* tmp_m_rvx + = current_term[0]->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix>* tmp_m_rvy + = current_term[1]->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix>* tmp_m_rvz + = current_term[2]->find_matrix(iat1, iat2, Rx, Ry, Rz); if (tmp_matrix_real == nullptr) { continue; @@ -254,16 +258,16 @@ void ModuleIO::write_current(const int istep, rvy = tmp_m_rvy->get_value(mu, nu); rvz = tmp_m_rvz->get_value(mu, nu); } - local_current_ik[0] -= dm2d1_real * rvx.real() - dm2d1_imag * rvx.imag(); + local_current_ik[0] -= dm2d1_real * rvx.real() - dm2d1_imag * rvx.imag(); local_current_ik[1] -= dm2d1_real * rvy.real() - dm2d1_imag * rvy.imag(); local_current_ik[2] -= dm2d1_real * rvz.real() - dm2d1_imag * rvz.imag(); - + ++local_total_irr; ++irr; } // end kk - } // end jj - } // end cb - } // end iat + } // end jj + } // end cb + } // end iat #ifdef _OPENMP #pragma omp critical(cal_current_k_reduce) { @@ -293,9 +297,8 @@ void ModuleIO::write_current(const int istep, fout.close(); } // write end - ModuleBase::timer::tick("ModuleIO", "write_current"); } // end nks - } // end is + } // end is if (GlobalV::MY_RANK == 0) { std::string filename = PARAM.globalv.global_out_dir + "current_total.dat"; @@ -306,10 +309,12 @@ void ModuleIO::write_current(const int istep, fout << istep << " " << current_total[0] << " " << current_total[1] << " " << current_total[2] << std::endl; fout.close(); } - if(!TD_Velocity::tddft_velocity) + if (!TD_Velocity::tddft_velocity) { delete cal_current; } + + ModuleBase::timer::tick("ModuleIO", "write_current"); return; } #endif //__LCAO From 4f2e4534a66084b3ed063a4acdc15fa4581b3988 Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Thu, 26 Sep 2024 12:02:21 +0800 Subject: [PATCH 2/3] Radial quadrature grid (#5173) * initial commit * ta3 * add unit test * more tests * minor change --- source/module_base/grid/delley.cpp | 35 ++--- source/module_base/grid/delley.h | 26 ++-- source/module_base/grid/radial.cpp | 69 ++++++++++ source/module_base/grid/radial.h | 70 ++++++++++ source/module_base/grid/test/CMakeLists.txt | 5 + source/module_base/grid/test/test_delley.cpp | 14 +- source/module_base/grid/test/test_radial.cpp | 136 +++++++++++++++++++ 7 files changed, 315 insertions(+), 40 deletions(-) create mode 100644 source/module_base/grid/radial.cpp create mode 100644 source/module_base/grid/radial.h create mode 100644 source/module_base/grid/test/test_radial.cpp diff --git a/source/module_base/grid/delley.cpp b/source/module_base/grid/delley.cpp index 1e0b6c703e..c24aa6ffbb 100644 --- a/source/module_base/grid/delley.cpp +++ b/source/module_base/grid/delley.cpp @@ -1,12 +1,13 @@ #include "module_base/grid/delley.h" +#include #include #include #include namespace { -struct Table { +struct DelleyTable { const int lmax_; const int ngrid_; const int ntype_[6]; @@ -14,7 +15,7 @@ struct Table { }; // Delley's table from the original article -const std::vector table = { +const std::vector delley_table = { { 17, 110, {1, 1, 0, 3, 1, 0}, { @@ -209,7 +210,7 @@ const std::vector
table = { 0.63942796347491023, 0.06424549224220589, 0.0009158016174693465, } } -}; // end of the definition of "table" +}; // end of the definition of "delley_table" // size of each group of points with octahedral symmetry // 6: (1, 0, 0) x sign x permutation (vertices) @@ -324,15 +325,15 @@ const std::vector fill = { }, }; // end of the definition of "fill" -const Table* _find_table(int& lmax) { - // NOTE: this function assumes elements in "Delley::table_" are +const DelleyTable* _find_delley(int& lmax) { + // NOTE: this function assumes elements in "delley_table" are // arranged such that their members "lmax_" are in ascending order. - auto tab = std::find_if(table.begin(), table.end(), - [lmax](const Table& t) { return t.lmax_ >= lmax; }); - return tab == table.end() ? nullptr : (lmax = tab->lmax_, &(*tab)); + auto tab = std::find_if(delley_table.begin(), delley_table.end(), + [lmax](const DelleyTable& t) { return t.lmax_ >= lmax; }); + return tab == delley_table.end() ? nullptr : (lmax = tab->lmax_, &(*tab)); } -void _get(const Table* tab, double* grid, double* weight) { +void _delley(const DelleyTable* tab, double* grid, double* weight) { assert(tab); const double* ptr = &tab->data_[0]; for (int itype = 0; itype < 6; ++itype) { @@ -348,29 +349,29 @@ void _get(const Table* tab, double* grid, double* weight) { } // end of anonymous namespace -int Grid::Delley::ngrid(int& lmax) { - auto tab = _find_table(lmax); +int Grid::Angular::ngrid_delley(int& lmax) { + auto tab = _find_delley(lmax); return tab ? tab->ngrid_ : -1; } -int Grid::Delley::get(int& lmax, double* grid, double* weight) { - auto tab = _find_table(lmax); - return tab ? _get(tab, grid, weight), 0 : -1; +int Grid::Angular::delley(int& lmax, double* grid, double* weight) { + auto tab = _find_delley(lmax); + return tab ? _delley(tab, grid, weight), 0 : -1; } -int Grid::Delley::get( +int Grid::Angular::delley( int& lmax, std::vector& grid, std::vector& weight ) { - auto tab = _find_table(lmax); + auto tab = _find_delley(lmax); if (!tab) { return -1; } grid.resize(3 * tab->ngrid_); weight.resize(tab->ngrid_); - _get(tab, grid.data(), weight.data()); + _delley(tab, grid.data(), weight.data()); return 0; } diff --git a/source/module_base/grid/delley.h b/source/module_base/grid/delley.h index 4f68a2168a..25f5fd35be 100644 --- a/source/module_base/grid/delley.h +++ b/source/module_base/grid/delley.h @@ -1,19 +1,10 @@ -#ifndef GRID_DELLEY_H -#define GRID_DELLEY_H +#ifndef GRID_ANGULAR_DELLEY_H +#define GRID_ANGULAR_DELLEY_H #include -#include -/** - * @brief Delley's grid for quadrature on the unit sphere. - * - * Reference: - * Delley, B. (1996). High order integration schemes on the unit sphere. - * Journal of computational chemistry, 17(9), 1152-1155. - * - */ namespace Grid { -namespace Delley { +namespace Angular { /** * @brief Number of Delley's grid points for a certain order of accuracy. @@ -27,7 +18,7 @@ namespace Delley { * lmax will be set to 23. * */ -int ngrid(int& lmax); +int ngrid_delley(int& lmax); /** @@ -45,13 +36,16 @@ int ngrid(int& lmax); * ngrid(lmax) elements, respectively. The function will return 0 if * successful, or -1 if the requested order cannot be fulfilled. * + * Reference: + * Delley, B. (1996). High order integration schemes on the unit sphere. + * Journal of computational chemistry, 17(9), 1152-1155. */ -int get(int& lmax, double* grid, double* weight); +int delley(int& lmax, double* grid, double* weight); // a handy wrapper doing the same as above -int get(int& lmax, std::vector& grid, std::vector& weight); +int delley(int& lmax, std::vector& grid, std::vector& weight); -} // end of namespace Delley +} // end of namespace Angular } // end of namespace Grid #endif diff --git a/source/module_base/grid/radial.cpp b/source/module_base/grid/radial.cpp new file mode 100644 index 0000000000..73a001eef5 --- /dev/null +++ b/source/module_base/grid/radial.cpp @@ -0,0 +1,69 @@ +#include "module_base/grid/radial.h" + +#include + +namespace { + +const double pi = std::acos(-1.0); +const double inv_ln2 = 1.0 / std::log(2.0); + +} // end of anonymous namespace + + +namespace Grid { +namespace Radial { + +void baker(int nbase, double R, double* r, double* w, int mult) { + int n = (nbase+1) * mult - 1; + double r0 = -R / std::log((2.0 * nbase + 1.0) / ((nbase+1)*(nbase+1))); + for (int i = 1; i <= n; ++i) { + r[i-1] = -r0 * std::log(1.0 - static_cast(i)*i/((n+1)*(n+1))); + w[i-1] = 2.0 * i * r0 * r[i-1] * r[i-1] / ((n+1+i)*(n+1-i)); + } +} + + +void baker(int nbase, double R, std::vector& r, + std::vector& w, int mult) { + int n = (nbase+1) * mult - 1; + r.resize(n); + w.resize(n); + baker(nbase, R, r.data(), w.data(), mult); +} + + +void murray(int n, double R, double* r, double* w) { + for (int i = 1; i <= n; ++i) { + double x = static_cast(i) / (n + 1); + r[i-1] = std::pow(x / (1.0 - x), 2) * R; + w[i-1] = 2.0 / (n + 1) * std::pow(R, 3) * std::pow(x, 5) + / std::pow(1.0 - x, 7); + } +} + + +void treutler_m4(int n, double R, double* r, double* w, double alpha) { + for (int i = 1; i <= n; ++i) { + double x = std::cos(i * pi / (n + 1)); + double beta = std::sqrt((1.0 + x) / (1.0 - x)); + double gamma = std::log((1.0 - x) / 2.0); + double delta = std::pow(1.0 + x, alpha); + r[i-1] = -R * inv_ln2 * delta * gamma; + w[i-1] = pi / (n + 1) * std::pow(delta * R * inv_ln2, 3) + * gamma * gamma * (beta - alpha / beta * gamma); + } +} + + +void mura(int n, double R, double* r, double* w) { + for (int i = 1; i <= n; ++i) { + double x = static_cast(i) / (n + 1); + double alpha = 1.0 - x * x * x; + r[i-1] = -R * std::log(alpha); + w[i-1] = 3.0 * R * std::pow(x * r[i-1], 2) / ((n+1) * alpha); + } +} + + +} // end of namespace Radial +} // end of namespace Grid diff --git a/source/module_base/grid/radial.h b/source/module_base/grid/radial.h new file mode 100644 index 0000000000..b8378f5834 --- /dev/null +++ b/source/module_base/grid/radial.h @@ -0,0 +1,70 @@ +#ifndef GRID_RADIAL_H +#define GRID_RADIAL_H + +#include + +namespace Grid { +namespace Radial { + +/** + * @brief Radial quadratures. + * + * This namespace contains functions that generate grids and weights + * for numerical integration + * + * / inf 2 + * | dr r g(r) ~ \sum_i w[i] g(r[i]) + * / 0 + * + */ + +/** + * Baker, J., Andzelm, J., Scheiner, A., & Delley, B. (1994). + * The effect of grid quality and weight derivatives in + * density functional calculations. + * The Journal of chemical physics, 101(10), 8894-8902. + * + * Zhang, I. Y., Ren, X., Rinke, P., Blum, V., & Scheffler, M. (2013). + * Numeric atom-centered-orbital basis sets with valence-correlation + * consistency from H to Ar. + * New Journal of Physics, 15(12), 123033. + * + * @note nbase is the number of points of the "base" grid, i.e., + * before applying the "radial multiplier" introduced by Zhang et al. + * The true number of grid points is (nbase+1) * mult - 1. + */ +void baker(int nbase, double R, double* r, double* w, int mult = 1); +void baker(int nbase, double R, std::vector& r, + std::vector& w, int mult = 1); + + +/** + * Murray, C. W., Handy, N. C., & Laming, G. J. (1993). + * Quadrature schemes for integrals of density functional theory. + * Molecular Physics, 78(4), 997-1014. + */ +void murray(int n, double R, double* r, double* w); + + +/** + * Treutler, O., & Ahlrichs, R. (1995). + * Efficient molecular numerical integration schemes. + * The Journal of Chemical Physics, 102(1), 346-354. + * + * @note M4 reduces to M3 at alpha = 0. + */ +void treutler_m4(int n, double R, double* r, double* w, double alpha = 0.6); + + +/** + * Mura, M. E., & Knowles, P. J. (1996). + * Improved radial grids for quadrature in molecular + * density‐functional calculations. + * The Journal of chemical physics, 104(24), 9848-9858. + */ +void mura(int n, double R, double* r, double* w); + +} // end of namespace Radial +} // end of namespace Grid + +#endif diff --git a/source/module_base/grid/test/CMakeLists.txt b/source/module_base/grid/test/CMakeLists.txt index 0f77a425fe..2a806935ad 100644 --- a/source/module_base/grid/test/CMakeLists.txt +++ b/source/module_base/grid/test/CMakeLists.txt @@ -7,3 +7,8 @@ AddTest( ../../ylm.cpp ) +AddTest( + TARGET test_radial + SOURCES test_radial.cpp + ../radial.cpp +) diff --git a/source/module_base/grid/test/test_delley.cpp b/source/module_base/grid/test/test_delley.cpp index 2b892cd532..66b99c5906 100644 --- a/source/module_base/grid/test/test_delley.cpp +++ b/source/module_base/grid/test/test_delley.cpp @@ -7,7 +7,7 @@ #include #endif -using namespace Grid; +using namespace Grid::Angular; // mock the function to prevent unnecessary dependency namespace ModuleBase { @@ -47,27 +47,27 @@ void DelleyTest::randgen(int lmax, std::vector& coef) { TEST_F(DelleyTest, NumGrid) { int lmax = 5; - int ngrid = Delley::ngrid(lmax); + int ngrid = ngrid_delley(lmax); EXPECT_EQ(lmax, 17); EXPECT_EQ(ngrid, 110); lmax = 17; - ngrid = Delley::ngrid(lmax); + ngrid = ngrid_delley(lmax); EXPECT_EQ(lmax, 17); EXPECT_EQ(ngrid, 110); lmax = 20; - ngrid = Delley::ngrid(lmax); + ngrid = ngrid_delley(lmax); EXPECT_EQ(lmax, 23); EXPECT_EQ(ngrid, 194); lmax = 59; - ngrid = Delley::ngrid(lmax); + ngrid = ngrid_delley(lmax); EXPECT_EQ(lmax, 59); EXPECT_EQ(ngrid, 1202); lmax = 60; - ngrid = Delley::ngrid(lmax); + ngrid = ngrid_delley(lmax); EXPECT_EQ(lmax, 60); EXPECT_EQ(ngrid, -1); } @@ -91,7 +91,7 @@ TEST_F(DelleyTest, Accuracy) { std::vector grid, weight, coef; for (int grid_lmax = 17; grid_lmax < 60; grid_lmax +=6) { - Delley::get(grid_lmax, grid, weight); + delley(grid_lmax, grid, weight); int func_lmax = grid_lmax / 2; randgen(func_lmax, coef); diff --git a/source/module_base/grid/test/test_radial.cpp b/source/module_base/grid/test/test_radial.cpp new file mode 100644 index 0000000000..6cdd4c4841 --- /dev/null +++ b/source/module_base/grid/test/test_radial.cpp @@ -0,0 +1,136 @@ +#include "module_base/grid/radial.h" + +#include "gtest/gtest.h" +#include +#include +#ifdef __MPI +#include +#endif + +using namespace Grid::Radial; +using Func_t = std::function; + +/** + * This test briefly checks various radial quadrature schemes by comparing + * their numerical results with analytical values on a few simple functions. + * + * The number of grid points and scaling factor are not carefully selected + * and the test is not exhaustive. It is just a sanity check. + * + */ + +const double pi = std::acos(-1.0); + +// test functions f(r) and their analytical integrals +// +// / inf 2 +// | dr r f(r) +// / 0 +// +std::vector> test_func_ref = { + { + [](double r) { + return std::exp(-0.3 * r * r) + std::exp(-3.0 * r * r); + }, + 0.25 * std::sqrt(pi) * (std::pow(0.3, -1.5) + std::pow(3.0, -1.5)) + }, + { + [](double r) { + return r * (std::exp(-0.3 * r * r) + std::exp(-3.0 * r * r)); + }, + 0.5 / (0.3 * 0.3) + 0.5 / (3.0 * 3.0) + }, + { + [](double r) { + return r * r * (std::exp(-0.3 * r * r) + std::exp(-3.0 * r * r)); + }, + 0.375 * std::sqrt(pi) * (std::pow(0.3, -2.5) + std::pow(3.0, -2.5)) + }, +}; + + +double quadrature(const Func_t& f, int n, double* r, double* w) { + double res = 0.0; + for (int i = 0; i < n; i++) { + res += w[i] * f(r[i]); + } + return res; +} + + +TEST(RadialTest, Baker) { + // R should be large enough to cover the range of the function. + // For mult = 1, R is the cutoff radius; for mult > 1, there + // are (mult - 1) grid points extend beyond R. + int nbase = 20; + int mult = 2; + double R = 7.0; + std::vector r, w; + baker(nbase, R, r, w, mult); + + EXPECT_EQ(r.size(), (nbase + 1) * mult - 1); + + for (auto& t : test_func_ref) { + double res = quadrature(t.first, r.size(), r.data(), w.data()); + EXPECT_NEAR(res, t.second, 1.0e-6); + } +} + + +TEST(RadialTest, Murray) { + int n = 40; + double R = 7.0; + std::vector r(n), w(n); + murray(n, R, r.data(), w.data()); + + for (auto& t : test_func_ref) { + double res = quadrature(t.first, r.size(), r.data(), w.data()); + EXPECT_NEAR(res, t.second, 1.0e-6); + } +} + + +TEST(RadialTest, Treutler) { + int n = 40; + double R = 7.0; + std::vector r(n), w(n); + + for (auto alpha : {0.0, 0.6, 1.0}) { + treutler_m4(n, R, r.data(), w.data(), alpha); + + for (auto& t : test_func_ref) { + double res = quadrature(t.first, r.size(), r.data(), w.data()); + EXPECT_NEAR(res, t.second, 1.0e-6); + } + } +} + + +TEST(RadialTest, Mura) { + int n = 40; + double R = 7.0; + std::vector r(n), w(n); + mura(n, R, r.data(), w.data()); + + for (auto& t : test_func_ref) { + double res = quadrature(t.first, r.size(), r.data(), w.data()); + EXPECT_NEAR(res, t.second, 1.0e-6); + } +} + + +int main(int argc, char** argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + + testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + +#ifdef __MPI + MPI_Finalize(); +#endif + + return result; +} From 3b2e5ada92114dab968d7b227a869ad1f9c4104b Mon Sep 17 00:00:00 2001 From: a1henu Date: Thu, 26 Sep 2024 12:04:26 +0800 Subject: [PATCH 3/3] Fix: Update function calls in `pyabacus` to align with new function signature in `hpsi_func` (#5176) * fix some typos in `_hsolver.py` * fix some bugs caused by #5134 --- python/pyabacus/src/py_diago_dav_subspace.hpp | 13 ++++---- python/pyabacus/src/py_diago_david.hpp | 4 +-- .../pyabacus/src/pyabacus/hsolver/_hsolver.py | 30 ++++++++++--------- 3 files changed, 26 insertions(+), 21 deletions(-) diff --git a/python/pyabacus/src/py_diago_dav_subspace.hpp b/python/pyabacus/src/py_diago_dav_subspace.hpp index f6cf7e9a5b..bd8bbb3e41 100644 --- a/python/pyabacus/src/py_diago_dav_subspace.hpp +++ b/python/pyabacus/src/py_diago_dav_subspace.hpp @@ -110,11 +110,14 @@ class PyDiagoDavSubspace bool scf_type, hsolver::diag_comm_info comm_info ) { - auto hpsi_func = [mm_op] (std::complex *hpsi_out, - std::complex *psi_in, const int nband_in, - const int nbasis_in, const int band_index1, - const int band_index2) - { + auto hpsi_func = [mm_op] ( + std::complex *psi_in, + std::complex *hpsi_out, + const int nband_in, + const int nbasis_in, + const int band_index1, + const int band_index2 + ) { // Note: numpy's py::array_t is row-major, but // our raw pointer-array is column-major py::array_t, py::array::f_style> psi({nbasis_in, band_index2 - band_index1 + 1}); diff --git a/python/pyabacus/src/py_diago_david.hpp b/python/pyabacus/src/py_diago_david.hpp index a53147f6ce..2db5284f4b 100644 --- a/python/pyabacus/src/py_diago_david.hpp +++ b/python/pyabacus/src/py_diago_david.hpp @@ -109,8 +109,8 @@ class PyDiagoDavid hsolver::diag_comm_info comm_info ) { auto hpsi_func = [mm_op] ( - std::complex *hpsi_out, - std::complex *psi_in, + std::complex *psi_in, + std::complex *hpsi_out, const int nband_in, const int nbasis_in, const int band_index1, diff --git a/python/pyabacus/src/pyabacus/hsolver/_hsolver.py b/python/pyabacus/src/pyabacus/hsolver/_hsolver.py index 5f5875a055..3dea18ef36 100644 --- a/python/pyabacus/src/pyabacus/hsolver/_hsolver.py +++ b/python/pyabacus/src/pyabacus/hsolver/_hsolver.py @@ -16,7 +16,7 @@ def rank(self) -> int: ... def nproc(self) -> int: ... def dav_subspace( - mm_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]], + mvv_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]], init_v: NDArray[np.complex128], dim: int, num_eigs: int, @@ -32,9 +32,10 @@ def dav_subspace( Parameters ---------- - mm_op : Callable[[NDArray[np.complex128]], NDArray[np.complex128]], - The operator to be diagonalized, which is a function that takes a matrix as input - and returns a matrix mv_op(X) = H * X as output. + mvv_op : Callable[[NDArray[np.complex128]], NDArray[np.complex128]], + The operator to be diagonalized, which is a function that takes a set of + vectors X = [x1, ..., xN] as input and returns a matrix(vector block) + mvv_op(X) = H * X ([Hx1, ..., HxN]) as output. init_v : NDArray[np.complex128] The initial guess for the eigenvectors. dim : int @@ -68,8 +69,8 @@ def dav_subspace( v : NDArray[np.complex128] The eigenvectors corresponding to the eigenvalues. """ - if not callable(mm_op): - raise TypeError("mm_op must be a callable object.") + if not callable(mvv_op): + raise TypeError("mvv_op must be a callable object.") if is_occupied is None: is_occupied = [True] * num_eigs @@ -86,7 +87,7 @@ def dav_subspace( assert dav_ndim * num_eigs < dim * comm_info.nproc, "dav_ndim * num_eigs must be less than dim * comm_info.nproc." _ = _diago_obj_dav_subspace.diag( - mm_op, + mvv_op, pre_condition, dav_ndim, tol, @@ -103,7 +104,7 @@ def dav_subspace( return e, v def davidson( - mm_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]], + mvv_op: Callable[[NDArray[np.complex128]], NDArray[np.complex128]], init_v: NDArray[np.complex128], dim: int, num_eigs: int, @@ -119,9 +120,10 @@ def davidson( Parameters ---------- - mm_op : Callable[[NDArray[np.complex128]], NDArray[np.complex128]], - The operator to be diagonalized, which is a function that takes a matrix as input - and returns a matrix mv_op(X) = H * X as output. + mvv_op : Callable[[NDArray[np.complex128]], NDArray[np.complex128]], + The operator to be diagonalized, which is a function that takes a set of + vectors X = [x1, ..., xN] as input and returns a matrix(vector block) + mvv_op(X) = H * X ([Hx1, ..., HxN]) as output. init_v : NDArray[np.complex128] The initial guess for the eigenvectors. dim : int @@ -146,8 +148,8 @@ def davidson( v : NDArray[np.complex128] The eigenvectors corresponding to the eigenvalues. """ - if not callable(mm_op): - raise TypeError("mm_op must be a callable object.") + if not callable(mvv_op): + raise TypeError("mvv_op must be a callable object.") if init_v.ndim != 1 or init_v.dtype != np.complex128: init_v = init_v.flatten().astype(np.complex128, order='C') @@ -159,7 +161,7 @@ def davidson( comm_info = hsolver.diag_comm_info(0, 1) _ = _diago_obj_dav_subspace.diag( - mm_op, + mvv_op, pre_condition, dav_ndim, tol,