diff --git a/source/module_base/math_sphbes.cpp b/source/module_base/math_sphbes.cpp index 36491a7f790..e3f0d8edcb7 100644 --- a/source/module_base/math_sphbes.cpp +++ b/source/module_base/math_sphbes.cpp @@ -775,15 +775,34 @@ double Sphbes::sphbesj(const int l, const double x) } } +double Sphbes::dsphbesj(const int l, const double x) +{ + assert( l >= 0 ); + assert( x >= 0 ); + return l == 0 ? -sphbesj(1, x) : ( l * sphbesj(l - 1, x) - (l + 1) * sphbesj(l + 1, x) ) / (2 * l + 1); +} + void Sphbes::sphbesj(const int n, + const double* const r, + const double q, + const int l, + double* const jl) +{ + for (int i = 0; i != n; ++i) + { + jl[i] = Sphbes::sphbesj(l, q * r[i]); + } +} + +void Sphbes::dsphbesj(const int n, const double* const r, const double q, const int l, - double* const jl) + double* const djl) { for (int i = 0; i != n; ++i) { - jl[i] = Sphbes::sphbesj(l, q * r[i]); + djl[i] = Sphbes::dsphbesj(l, q * r[i]); } } diff --git a/source/module_base/math_sphbes.h b/source/module_base/math_sphbes.h index 5584a95f347..2568a065f7c 100644 --- a/source/module_base/math_sphbes.h +++ b/source/module_base/math_sphbes.h @@ -97,11 +97,25 @@ class Sphbes const double x //!< [in] argument ); - static void sphbesj(const int n, - const double* const r, - const double q, - const int l, - double* const jl + //! derivative of spherical Bessel function + static double dsphbesj(const int l, //!< [in] order + const double x //!< [in] argument + ); + + //! computes the values of l-th order spherical Bessel function at q*r[ir] + static void sphbesj(const int n, //!< [in] number of r grid points + const double* const r, //!< [in] r grid + const double q, //!< [in] wave vector + const int l, //!< [in] order of the spherical Bessel function + double* const jl //!< [out] results + ); + + //! computes the derivative of l-th order spherical Bessel function at q*r[ir] + static void dsphbesj(const int n, //!< [in] number of r grid points + const double* const r, //!< [in] r grid + const double q, //!< [in] wave vector + const int l, //!< [in] order of the spherical Bessel function + double* const djl //!< [out] results ); private: diff --git a/source/module_base/spherical_bessel_transformer.cpp b/source/module_base/spherical_bessel_transformer.cpp index 11a0fa9ea3a..ac8c161467d 100644 --- a/source/module_base/spherical_bessel_transformer.cpp +++ b/source/module_base/spherical_bessel_transformer.cpp @@ -19,6 +19,9 @@ SphericalBesselTransformer::~SphericalBesselTransformer() { fftw_destroy_plan(rfft_plan_); fftw_free(f_); + delete[] grid_in_; + delete[] grid_out_; + delete[] jl_; } long long int SphericalBesselTransformer::spherical_bessel_sincos_polycoef(const bool get_sine, @@ -73,49 +76,22 @@ long long int SphericalBesselTransformer::spherical_bessel_sincos_polycoef(const - (n >= 2 ? spherical_bessel_sincos_polycoef(get_sine, l - 2, n - 2) : 0); } -void SphericalBesselTransformer::radrfft(const int l, - const int ngrid, - const double cutoff, - const double* const in, - double* const out, - const int p) +void SphericalBesselTransformer::_radrfft_base(const int l, + const int ngrid, + const double cutoff, + const double* const in, + double* const out, + const int p) { - /* - * An l-th order spherical Bessel transform F(x) -> G(y) can be expressed in terms of Fourier transforms: - * - * l - * --- 1 / +inf -iyx - * G(y) = \ ------- Re | dr f(n,x) e - * / l+1-n / -inf - * --- y - * n=0 - * - * where - * - * 1 F(x) - * f(n,x) = ---------- [ c(l,n) + i*s(l,n) ] -------- - * sqrt(2*pi) l-n-1 , - * x - * - * c(l,n) / s(l,n) are sin / cos coefficients from spherical_bessel_sincos_polycoef, and - * the domain of F(x) is extended to negative values by defining F(-x) = pow(-1,l)*F(x). - * - * With an appropriate grid, the Fourier transform can be approximated by a discrete Fourier transform. - * Note that given l & n, c(l,n) and s(l,n) cannot be both non-zero. Therefore, each FFT input - * array is either purely real or purely imaginary, which suggests the use of real-input FFT. - * */ - assert(l >= 0); - assert(ngrid > 1); - assert(p <= 2); + // this function does not support in-place transform (but radrfft does) + assert(in != out); int n = ngrid - 1; double dx = cutoff / n; double dy = PI / cutoff; double pref = dx / std::sqrt(2. * PI); - // use a temporary output array in case the transform is in-place - double* out_tmp = new double[ngrid]; - std::fill(out_tmp, out_tmp + ngrid, 0.0); + std::fill(out, out + ngrid, 0.0); // rfft buffer (f_) allocation and plan creation rfft_prepare(2 * n); @@ -125,12 +101,11 @@ void SphericalBesselTransformer::radrfft(const int l, for (int m = 0; m <= l; ++m) { - // m even --> sin; f[2*n-i] = -f[i]; out += -imag(rfft(f)) / y^(l+1-m) // m odd --> cos; f[2*n-i] = +f[i]; out += +real(rfft(f)) / y^(l+1-m) bool flag = (m % 2 == 0); - int sincos_coef = spherical_bessel_sincos_polycoef(flag, l, m); + long long int sincos_coef = spherical_bessel_sincos_polycoef(flag, l, m); f[0] = f[n] = 0.0; for (int i = 1; i != n; ++i) @@ -145,7 +120,7 @@ void SphericalBesselTransformer::radrfft(const int l, // out[0] is handled independently for (int j = 1; j <= n; ++j) { - out_tmp[j] = (out_tmp[j] + (flag ? -f_[j][1] : f_[j][0])) / (j * dy); + out[j] = (out[j] + (flag ? -f_[j][1] : f_[j][0])) / (j * dy); } } @@ -155,20 +130,90 @@ void SphericalBesselTransformer::radrfft(const int l, { for (int i = 0; i <= n; ++i) { - out_tmp[0] += 2. * pref * in[i] * std::pow(i * dx, 2 - p); // p <= 2 is required! + out[0] += 2. * pref * in[i] * std::pow(i * dx, 2 - p); // p <= 2 is required! + } + } +} + +void SphericalBesselTransformer::radrfft(const int l, + const int ngrid, + const double cutoff, + const double* const in, + double* const out, + const int p, + const bool deriv) +{ + /* + * An l-th order spherical Bessel transform F(x) -> G(y) can be expressed in terms of Fourier transforms: + * + * l + * --- 1 / +inf -iyx + * G(y) = \ ------- Re | dr f(n,x) e + * / l+1-n / -inf + * --- y + * n=0 + * + * where + * + * 1 F(x) + * f(n,x) = ---------- [ c(l,n) + i*s(l,n) ] -------- + * sqrt(2*pi) l-n-1 , + * x + * + * c(l,n) / s(l,n) are sin / cos coefficients from spherical_bessel_sincos_polycoef, and + * the domain of F(x) is extended to negative values by defining F(-x) = pow(-1,l)*F(x). + * + * With an appropriate grid, the Fourier transform can be approximated by a discrete Fourier transform. + * Note that given l & n, c(l,n) and s(l,n) cannot be both non-zero. Therefore, each FFT input + * array is either purely real or purely imaginary, which suggests the use of real-input FFT. + * + * If deriv == true, dG(y)/dy is calculated instead of G(y). This is done by using the recurrence + * relation of j_l(x) to convert the original transform to two Spherical Bessel transforms of order + * l-1 and l+1 respectively. + * */ + assert(l >= 0); + assert(ngrid > 1); + assert(p <= 2); + + double* out_tmp = new double[deriv && l > 0 ? 2 * ngrid : ngrid]; + int n_direct = 0; + + if (deriv) + { + if (l == 0) + { + _radrfft_base(1, ngrid, cutoff, in, out_tmp, p - 1); + std::for_each(out_tmp, out_tmp + ngrid, [](double& x) { x = -x; }); + } + else + { + _radrfft_base(l - 1, ngrid, cutoff, in, out_tmp, p - 1); + _radrfft_base(l + 1, ngrid, cutoff, in, &out_tmp[ngrid], p - 1); + std::transform(out_tmp, out_tmp + ngrid, &out_tmp[ngrid], out_tmp, + [l](double x, double y) { return (l * x - (l + 1) * y) / (2 * l + 1); }); } + + n_direct = static_cast(ngrid * std::pow(1e-8, 1. / (l + 1))); + } + else + { + _radrfft_base(l, ngrid, cutoff, in, out_tmp, p); + n_direct = (l == 0) ? 0 : static_cast(ngrid * std::pow(1e-8, 1. / l)); } - // FFT-based method is not accurate for small k at high l + // FFT-based method does not yield accurate results for small y at high l // use numerical integration in this case - int n_direct = (l == 0) ? 0 : static_cast(ngrid * std::pow(1e-8, 1. / l)); if (n_direct > 0) { + double dx = cutoff / (ngrid - 1); + double dy = PI / cutoff; double* grid_in = new double[ngrid]; double* grid_out = new double[n_direct]; std::for_each(grid_in, grid_in + ngrid, [&](double& x) { x = (&x - grid_in) * dx; }); std::for_each(grid_out, grid_out + n_direct, [&](double& y) { y = ((&y - grid_out) + 1) * dy; }); - direct(l, ngrid, grid_in, in, n_direct, grid_out, &out_tmp[1], p); + + direct(l, ngrid, grid_in, in, n_direct, grid_out, &out_tmp[1], p, deriv); + delete[] grid_in; delete[] grid_out; } @@ -184,49 +229,43 @@ void SphericalBesselTransformer::direct(const int l, const int ngrid_out, const double* const grid_out, double* const out, - const int p) + const int p, + const bool deriv) { - // TODO: - // 1. avoid new/delete rab/r2f/integrand/jl every time; use some buffer class member instead - // 2. cache the spherical Bessel function values for given (l, grid_in, grid_out). This should be done - // alongside the change of for-loop sequence in building TwoCenterTable. - assert(p <= 2); assert(grid_in[0] >= 0.0 && std::is_sorted(grid_in, grid_in + ngrid_in, std::less_equal())); assert(grid_out[0] >= 0.0 && std::is_sorted(grid_out, grid_out + ngrid_out, std::less_equal())); - double pref = std::sqrt(2.0 / ModuleBase::PI); + double* buffer = new double[3 * ngrid_in]; // NOTE: std::adjacent_difference returns an array of the same size as the input array // with rab[0] = grid_in[0] and rab[i] = grid_in[i] - grid[i-1] for i > 0 - double* rab = new double[ngrid_in]; + double* rab = buffer; std::adjacent_difference(grid_in, grid_in + ngrid_in, rab); - // r^2*f(r) (integrand without j_l) - double* r2f = new double[ngrid_in]; - std::memcpy(r2f, in, ngrid_in * sizeof(double)); - std::for_each(r2f, r2f + ngrid_in, [&](double& x) { x *= std::pow(grid_in[&x - r2f], 2 - p); }); + // integrand without the jl (or jl') part + // note that (d/dy)jl(y*x) = x * (jl')(y*x), hence the extra power of x + double* itgpart = buffer + ngrid_in; + std::memcpy(itgpart, in, ngrid_in * sizeof(double)); + std::for_each(itgpart, itgpart + ngrid_in, [&](double& x) { x *= std::pow(grid_in[&x - itgpart], 2 - p + deriv); }); // there's no need to use a temporary output array even if the transform is in-place - // because "in" is not referenced after the construction of r2f. + // because "in" is not referenced after the construction of itgpart. - double* integrand = new double[ngrid_in]; + // cache the jl or jl' values for potential subsequent transforms on the same grid + cache(l, ngrid_in, grid_in, ngrid_out, grid_out, deriv); - // TODO cache the calculated jl - double* jl = new double[ngrid_in]; + double* integrand = buffer + 2 * ngrid_in; for (int i_out = 0; i_out != ngrid_out; ++i_out) { - ModuleBase::Sphbes::sphbesj(ngrid_in, grid_in, grid_out[i_out], l, jl); - std::transform(r2f, r2f + ngrid_in, jl, integrand, std::multiplies()); + std::transform(itgpart, itgpart + ngrid_in, &jl_[i_out * ngrid_in_], integrand, std::multiplies()); out[i_out] = ModuleBase::Integral::simpson(ngrid_in, integrand, &rab[1]); } - delete[] jl; + double pref = std::sqrt(2.0 / ModuleBase::PI); std::for_each(out, out + ngrid_out, [pref](double& x) { x *= pref; }); - delete[] r2f; - delete[] integrand; - delete[] rab; + delete[] buffer; } void SphericalBesselTransformer::rfft_in_place() @@ -264,4 +303,45 @@ void SphericalBesselTransformer::set_fftw_plan_flag(const unsigned new_flag) } } +void SphericalBesselTransformer::cache(const int l, + const int ngrid_in, + const double* const grid_in, + const int ngrid_out, + const double* const grid_out, + const bool deriv) +{ + // TODO: if the given input grid matchs the cached output grid and vice versa, + // we can just swap the cached grids and transpose jl_. + bool is_cached = deriv == is_deriv_ && l == l_ && ngrid_in <= ngrid_in_ && ngrid_out <= ngrid_out_ + && std::equal(grid_in, grid_in + ngrid_in, grid_in_) + && std::equal(grid_out, grid_out + ngrid_out, grid_out_); + + if (!is_cached) + { + delete[] jl_; + delete[] grid_in_; + delete[] grid_out_; + + is_deriv_ = deriv; + l_ = l; + ngrid_in_ = ngrid_in; + ngrid_out_ = ngrid_out; + + grid_in_ = new double[ngrid_in]; + grid_out_ = new double[ngrid_out]; + std::memcpy(grid_in_, grid_in, ngrid_in * sizeof(double)); + std::memcpy(grid_out_, grid_out, ngrid_out * sizeof(double)); + + typedef void(*FuncType)(int, const double*, double, int, double*); + FuncType f = is_deriv_ ? static_cast(ModuleBase::Sphbes::dsphbesj) : + static_cast(ModuleBase::Sphbes::sphbesj); + + jl_ = new double[ngrid_in * ngrid_out]; + for (int i_out = 0; i_out != ngrid_out; ++i_out) + { + f(ngrid_in, grid_in, grid_out[i_out], l, &jl_[i_out * ngrid_in]); + } + } +} + } // namespace ModuleBase diff --git a/source/module_base/spherical_bessel_transformer.h b/source/module_base/spherical_bessel_transformer.h index a8ba2203729..e5bdeee4ec9 100644 --- a/source/module_base/spherical_bessel_transformer.h +++ b/source/module_base/spherical_bessel_transformer.h @@ -11,12 +11,12 @@ namespace ModuleBase * The spherical Bessel transform of a function F(x) is defined as * * / +inf 2 - * G(y) = sqrt(2/pi) * | dx x F(x) j (x) + * G(y) = sqrt(2/pi) * | dx x F(x) j (y*x) * / 0 l * * where * - * j (x) + * j * l * * is the l-th order spherical Bessel function of the first kind. @@ -28,7 +28,7 @@ namespace ModuleBase * * and, on finish, fills the output array with * - * out[j] = G(y[j]) + * out[j] = G(y[j]) or out[j] = dG(y[j])/dy * * Usage1: * @@ -80,7 +80,8 @@ class SphericalBesselTransformer //! Performs an l-th order spherical Bessel transform via real-input fast Fourier transforms. /*! - * This function computes the spherical Bessel transform F(x) -> G(y) with input values + * This function computes the spherical Bessel transform F(x) -> G(y) (or G's derivative) + * with input values * * p * in[i] = x [i] F(x[i]) @@ -91,7 +92,7 @@ class SphericalBesselTransformer * x[i] = i * ------- i = 0, 1, 2,..., ngrid-1. * ngrid-1 * - * On finish, out[j] = G(y[j]) where + * On finish, out[j] = G(y[j]) or dG(y[j])/dy where * * pi * y[j] = j * ------ j = 0, 1, 2,..., ngrid-1. @@ -103,45 +104,47 @@ class SphericalBesselTransformer * no sense if the input is truncated at a place where F(x) is still significantly * non-zero. * @note FFT-based algorithm is not accurate for high l at small y. Numerical - * integration in used to handle this case. + * integration is automatically invoked to handle this case. * @note p must not exceed 2 as required by numerical integration, @see @ref direct * */ - void radrfft(const int l, //!< [in] order of the transform - const int ngrid, //!< [in] size of the input array - const double cutoff, //!< [in] cutoff distance of input values - const double* const in, //!< [in] input values - double* const out, //!< [out] transformed values - const int p = 0 //!< [in] exponent of the pre-multiplied power term in input values, - //!< must not exceed 2 + void radrfft(const int l, //!< [in] order of the transform + const int ngrid, //!< [in] size of the input array + const double cutoff, //!< [in] cutoff distance of input values + const double* const in, //!< [in] input values + double* const out, //!< [out] transformed values + const int p = 0, //!< [in] exponent of the extra power term in input values, must not exceed 2 + const bool deriv = false //!< [in] whether to compute the derivative of the transform ); //! Performs an l-th order spherical Bessel transform via numerical integration with Simpson's rule. /*! - * This function computes the spherical Bessel transform F(x) -> G(y) with input values + * This function computes the spherical Bessel transform F(x) -> G(y) (or G's derivative) + * with input values * * p * in[i] = x [i] F(x[i]) * - * where p <= 2 is an integer. On finish, out[j] = G(y[j]). x & y are specified by - * grid_in & grid_out, respectively. + * where p <= 2 is an integer. On finish, out[j] = G(y[j]) or dG(y[j])/dy. + * x & y are specified by grid_in & grid_out, respectively. * * @note This function does not allocate memory for output; it must be pre-allocated. * @note Even if the input grid forms a good sampling of F(x), results would still be * inaccurate for very large y values (y*dx ~ pi) because the oscillation of * j_l(y*x) in this case is poorly sampled, in which case Simpson's 1/3 rule * could be a bad approximation. - * @note p is limited to p <= 2 in order to avoid the problem of determining - * x^2*F(x) at x = 0 from x[i]^p*F(x[i]) + * @note p is restricted to p <= 2 in order to avoid the situation that one has to + * determine x^2*F(x) at x = 0 from x[i]^p*F(x[i]). * */ - static void direct(const int l, //!< [in] order of the transform - const int ngrid_in, //!< [in] size of the input array - const double* const grid_in, //!< [in] input grid - const double* const in, //!< [in] input values - const int ngrid_out, //!< [in] size of the output array - const double* const grid_out, //!< [in] output grid - double* const out, //!< [out] transformed values - const int p = 0 //!< [in] exponent of the pre-multiplied power - //!< term in input values, must not exceed 2 + void direct(const int l, //!< [in] order of the transform + const int ngrid_in, //!< [in] size of the input array + const double* const grid_in, //!< [in] input grid + const double* const in, //!< [in] input values + const int ngrid_out, //!< [in] size of the output array + const double* const grid_out, //!< [in] output grid + double* const out, //!< [out] transformed values + const int p = 0, //!< [in] exponent of the extra power term in input values, + //!< must not exceed 2! + const bool deriv = false //!< [in] whether to compute the derivative of the transform ); //! Sets the FFTW planner flag. @@ -160,6 +163,15 @@ class SphericalBesselTransformer void set_fftw_plan_flag(const unsigned new_flag /*!< [in] FFTW planner flag */); private: + //! core function for performing the transform with FFT + void _radrfft_base(const int l, + const int ngrid, + const double cutoff, + const double* const in, + double* const out, + const int p = 0 + ); + //! Internal buffer used for in-place real-input FFT (interpreted as double* on input) fftw_complex* f_ = nullptr; @@ -203,6 +215,30 @@ class SphericalBesselTransformer const int n //!< [in] degree of the polynomial term whose coefficient is computed ); + //! Computes & stores the spherical Bessel function of the first kind on the given transform grid + void cache(const int l, + const int ngrid_in, + const double* const grid_in, + const int ngrid_out, + const double* const grid_out, + const bool deriv); + + /*! + * @name Cached function values + * */ + //!@{ + bool is_deriv_ = false; //!< if true, the cached values are derivatives of the spherical Bessel function + int l_ = -1; //!< order of the cached spherical Bessel function + + int ngrid_in_ = 0; + int ngrid_out_ = 0; + double* grid_in_ = nullptr; + double* grid_out_ = nullptr; + + //! jl_[i*ngrid_in_ + j] = f(l, grid_out_[i] * grid_in_[j]) where f is sphbesj or dsphbesj + double* jl_ = nullptr; + //!@} + }; // class SphericalBesselTransformer } // namespace ModuleBase diff --git a/source/module_base/test/math_sphbes_test.cpp b/source/module_base/test/math_sphbes_test.cpp index e659f6e5c80..a56892f3d6d 100644 --- a/source/module_base/test/math_sphbes_test.cpp +++ b/source/module_base/test/math_sphbes_test.cpp @@ -225,6 +225,14 @@ TEST_F(Sphbes, SeriesAndRecurrence) { EXPECT_NEAR(jl_old[i], jl_new[i], 1e-12); } + + // derivative + ModuleBase::Sphbes::dSpherical_Bessel_dx(nr, r, q, 0, jl_old); + ModuleBase::Sphbes::dsphbesj(nr, r, q, 0, jl_new); + for (int i = 0; i < nr; ++i) + { + EXPECT_NEAR(jl_old[i], jl_new[i], 1e-12); + } } } diff --git a/source/module_base/test/spherical_bessel_transformer_test.cpp b/source/module_base/test/spherical_bessel_transformer_test.cpp index 093b3886a5b..72498244caa 100644 --- a/source/module_base/test/spherical_bessel_transformer_test.cpp +++ b/source/module_base/test/spherical_bessel_transformer_test.cpp @@ -284,7 +284,9 @@ TEST_F(SphericalBesselTransformTest, DirectBasic) y = pref * (3.0 - k * k) / std::pow(k * k + 1, 3); }); - ModuleBase::SphericalBesselTransformer::direct(0, sz_in, grid_in, f, sz_out, grid_out, g); + ModuleBase::SphericalBesselTransformer sbt; + + sbt.direct(0, sz_in, grid_in, f, sz_out, grid_out, g); EXPECT_LT(max_diff(sz_out, g_ref, g), tol); // first-order transform @@ -293,7 +295,7 @@ TEST_F(SphericalBesselTransformTest, DirectBasic) y = pref * 4.0 * k / std::pow(k * k + 1, 3); }); - ModuleBase::SphericalBesselTransformer::direct(1, sz_in, grid_in, f, sz_out, grid_out, g); + sbt.direct(1, sz_in, grid_in, f, sz_out, grid_out, g); EXPECT_LT(max_diff(sz_out, g_ref, g), tol); // second-order transform @@ -302,7 +304,7 @@ TEST_F(SphericalBesselTransformTest, DirectBasic) y = pref * 4.0 * k * k / std::pow(k * k + 1, 3); }); - ModuleBase::SphericalBesselTransformer::direct(2, sz_in, grid_in, f, sz_out, grid_out, g); + sbt.direct(2, sz_in, grid_in, f, sz_out, grid_out, g); EXPECT_LT(max_diff(sz_out, g_ref, g), tol); } @@ -331,6 +333,8 @@ TEST_F(SphericalBesselTransformTest, DirectImplicitExponent) y = pref * k * k / std::pow(k * k + 1, 4); }); + ModuleBase::SphericalBesselTransformer sbt; + for (int p = -2; p <= 2; ++p) { std::for_each(f, f + sz_in, [&](double& x) { @@ -338,7 +342,7 @@ TEST_F(SphericalBesselTransformTest, DirectImplicitExponent) x = std::pow(r, 2 + p) * std::exp(-r); }); - ModuleBase::SphericalBesselTransformer::direct(2, sz_in, grid_in, f, sz_out, grid_out, g, p); + sbt.direct(2, sz_in, grid_in, f, sz_out, grid_out, g, p); EXPECT_LT(max_diff(sz_out, g_ref, g), tol); } } @@ -372,7 +376,9 @@ TEST_F(SphericalBesselTransformTest, DirectInPlace) y = pref * k * k * std::exp(-k * k / 4); }); - ModuleBase::SphericalBesselTransformer::direct(2, sz_in, grid_in, f, sz_out, grid_out, f); + ModuleBase::SphericalBesselTransformer sbt; + + sbt.direct(2, sz_in, grid_in, f, sz_out, grid_out, f); EXPECT_LT(max_diff(sz_out, g_ref, f), tol); } @@ -403,7 +409,7 @@ TEST_F(SphericalBesselTransformTest, HighOrder) SphericalBesselTransformer sbt; sbt.radrfft(l, sz, rcut, f, g); - ModuleBase::SphericalBesselTransformer::direct(l, sz, grid_in, f, sz, grid_out, g_ref); + sbt.direct(l, sz, grid_in, f, sz, grid_out, g_ref); // NOTE: Simpson's integration gets increasingly inaccurate as k gets large // since the factor of (k*dr)^4 in its error becomes significant when k*dr diff --git a/source/module_basis/module_nao/atomic_radials.cpp b/source/module_basis/module_nao/atomic_radials.cpp index 5629527b460..18877830ab1 100644 --- a/source/module_basis/module_nao/atomic_radials.cpp +++ b/source/module_basis/module_nao/atomic_radials.cpp @@ -16,7 +16,7 @@ AtomicRadials& AtomicRadials::operator=(const AtomicRadials& rhs) void AtomicRadials::build(const std::string& file, const int itype, std::ofstream* ptr_log, const int rank) { - // deallocates all arrays and reset variables + // deallocates all arrays and reset variables (excluding sbt_) cleanup(); std::ifstream ifs; @@ -58,6 +58,11 @@ void AtomicRadials::build(const std::string& file, const int itype, std::ofstrea { ifs.close(); } + + for (int i = 0; i < nchi_; i++) + { + chi_[i].set_transformer(sbt_, 0); + } } void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) diff --git a/source/module_basis/module_nao/beta_radials.cpp b/source/module_basis/module_nao/beta_radials.cpp index 782d4e0e2f6..b36386c2c58 100644 --- a/source/module_basis/module_nao/beta_radials.cpp +++ b/source/module_basis/module_nao/beta_radials.cpp @@ -3,7 +3,8 @@ #include "module_base/parallel_common.h" #include "module_base/tool_quit.h" -void BetaRadials::build(const std::string& file, const int itype, std::ofstream* ptr_log, const int rank) { +void BetaRadials::build(const std::string& file, const int itype, std::ofstream* ptr_log, const int rank) +{ /* * Build a BetaRadials object from beta functions read from a pseudopotential file * @@ -48,10 +49,12 @@ void BetaRadials::build(const std::string& file, const int itype, std::ofstream* // identify the version of UPF file // UPF 2.0.1 format starts with int upf_version = 100; - if (rank == 0) { + if (rank == 0) + { std::string first_line; std::getline(ifs, first_line); - if (first_line.find("2.0.1") != std::string::npos) { + if (first_line.find("2.0.1") != std::string::npos) + { upf_version = 201; } } @@ -60,23 +63,30 @@ void BetaRadials::build(const std::string& file, const int itype, std::ofstream* Parallel_Common::bcast_int(upf_version); #endif - switch (upf_version) { - case 100: - read_beta_upf100(ifs, ptr_log, rank); - break; - case 201: - read_beta_upf201(ifs, ptr_log, rank); - break; - default: /* not supposed to happen */; + switch (upf_version) + { + case 100: + read_beta_upf100(ifs, ptr_log, rank); + break; + case 201: + read_beta_upf201(ifs, ptr_log, rank); + break; + default: /* not supposed to happen */; } if (rank == 0) { ifs.close(); } + + for (int i = 0; i < nchi_; i++) + { + chi_[i].set_transformer(sbt_, 0); + } } -void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) { +void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) +{ /* * Read the nonlocal beta functions from the pseudopotential file (in old UPF format) * */ @@ -86,7 +96,8 @@ void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, c int ngrid_max = 0; bool is_good = false; - if (rank == 0) { + if (rank == 0) + { /* * Read the header, including * @@ -95,27 +106,36 @@ void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, c * 3. Number of radial grid points * 4. Number of beta functions * */ - while (std::getline(ifs, line)) { - if (line.find("Element") != std::string::npos) { + while (std::getline(ifs, line)) + { + if (line.find("Element") != std::string::npos) + { ss.str(""); ss << line; ss >> symbol_; - } else if (line.find("Max angular momentum component") != std::string::npos) { + } + else if (line.find("Max angular momentum component") != std::string::npos) + { ss.str(""); ss << line; ss >> lmax_; - } else if (line.find("Number of points in mesh") != std::string::npos) { + } + else if (line.find("Number of points in mesh") != std::string::npos) + { ss.str(""); ss << line; ss >> ngrid_max; - } else if (line.find("Number of Projectors") != std::string::npos) { + } + else if (line.find("Number of Projectors") != std::string::npos) + { ss.str(""); ss << line; ss >> tmp >> nchi_; // nchi_ is the number of beta functions - } + } // should obtain valid nchi_, symbol_ & ngrid_max upon reaching the end of header - if (line.find("") != std::string::npos) { + if (line.find("") != std::string::npos) + { // lmax could be -1 when there is no beta projectors, see, e.g., H.pz-vbc.UPF is_good = (nchi_ >= 0) && symbol_.size() && (ngrid_max > 0); break; @@ -137,23 +157,28 @@ void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, c } // In case some pseudopotential file does not have any beta function - if (nchi_ == 0) { + if (nchi_ == 0) + { return; } double* rgrid = new double[ngrid_max]; - if (rank == 0) { + if (rank == 0) + { /* * Read the radial grid * */ - while (ifs >> tmp) { - if (tmp == "") { + while (ifs >> tmp) + { + if (tmp == "") + { break; } } assert(!ifs.eof()); - for (int ir = 0; ir != ngrid_max; ++ir) { + for (int ir = 0; ir != ngrid_max; ++ir) + { ifs >> rgrid[ir]; } @@ -167,7 +192,8 @@ void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, c assert(lmax_ >= 0); nzeta_ = new int[lmax_ + 1]; - for (int l = 0; l <= lmax_; ++l) { + for (int l = 0; l <= lmax_; ++l) + { nzeta_[l] = 0; } @@ -180,8 +206,10 @@ void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, c int izeta = 0; chi_ = new NumericalRadial[nchi_]; - for (int i = 0; i != nchi_; ++i) { - if (rank == 0) { + for (int i = 0; i != nchi_; ++i) + { + if (rank == 0) + { /* * Read the beta functions, including * @@ -191,8 +219,10 @@ void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, c * * and record the zeta number (izeta). * */ - while (std::getline(ifs, line)) { - if (line.find("") != std::string::npos) { + while (std::getline(ifs, line)) + { + if (line.find("") != std::string::npos) + { break; } } @@ -201,14 +231,18 @@ void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, c ifs >> tmp >> tmp; // skip "Beta" "L" ifs >> ngrid; - if (l == l_last) { + if (l == l_last) + { izeta += 1; - } else { + } + else + { izeta = 0; l_last = l; } - for (int ir = 0; ir != ngrid; ++ir) { + for (int ir = 0; ir != ngrid; ++ir) + { ifs >> rbeta[ir]; } @@ -237,7 +271,8 @@ void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, c delete[] rbeta; } -void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) { +void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) +{ /* * Read the nonlocal beta functions from the pseudopotential file (in UPF 2.0.1 format) * */ @@ -246,7 +281,8 @@ void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, c int ngrid_max = 0; bool is_good = false; - if (rank == 0) { + if (rank == 0) + { /* * Read the header, including * @@ -255,18 +291,27 @@ void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, c * 3. Number of radial grid points * 4. Number of beta functions * */ - while (std::getline(ifs, line)) { - if (line.find("element=") != std::string::npos) { + while (std::getline(ifs, line)) + { + if (line.find("element=") != std::string::npos) + { symbol_ = trim201(line); - } else if (line.find("l_max=") != std::string::npos) { + } + else if (line.find("l_max=") != std::string::npos) + { lmax_ = std::stoi(trim201(line)); - } else if (line.find("mesh_size=") != std::string::npos) { + } + else if (line.find("mesh_size=") != std::string::npos) + { ngrid_max = std::stoi(trim201(line)); - } else if (line.find("number_of_proj=") != std::string::npos) { + } + else if (line.find("number_of_proj=") != std::string::npos) + { nchi_ = std::stoi(trim201(line)); } - if (line.find("/>") != std::string::npos) { + if (line.find("/>") != std::string::npos) + { is_good = (nchi_ >= 0) && (ngrid_max >= 0) && symbol_.size(); break; } @@ -287,23 +332,28 @@ void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, c } // In case some pseudopotential file does not have any beta function - if (nchi_ == 0) { + if (nchi_ == 0) + { return; } double* rgrid = new double[ngrid_max]; - if (rank == 0) { + if (rank == 0) + { /* * Read the radial grid * */ - while (std::getline(ifs, line)) { - if (line.find("> rgrid[ir]; } @@ -316,7 +366,8 @@ void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, c #endif nzeta_ = new int[lmax_ + 1]; - for (int l = 0; l <= lmax_; ++l) { + for (int l = 0; l <= lmax_; ++l) + { nzeta_[l] = 0; } @@ -329,8 +380,10 @@ void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, c int izeta = 0; chi_ = new NumericalRadial[nchi_]; - for (int i = 0; i != nchi_; ++i) { - if (rank == 0) { + for (int i = 0; i != nchi_; ++i) + { + if (rank == 0) + { /* * Read the beta functions, including * @@ -340,42 +393,55 @@ void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, c * * and record the zeta number (izeta). * */ - while (std::getline(ifs, line)) { - if (line.find("") != std::string::npos) { + if (line.find(">") != std::string::npos) + { is_good &= (l >= 0) && (l <= lmax_) && (ngrid > 0) && (ngrid <= ngrid_max); break; } } - if (l == l_last) { + if (l == l_last) + { izeta += 1; - } else { + } + else + { izeta = 0; l_last = l; } - for (int ir = 0; ir != ngrid; ++ir) { + for (int ir = 0; ir != ngrid; ++ir) + { ifs >> rbeta[ir]; } // reverse scan to determine the grid size - for (; ngrid > 0; --ngrid) { - if (std::abs(rbeta[ngrid-1]) > 1e-12) { + for (; ngrid > 0; --ngrid) + { + if (std::abs(rbeta[ngrid - 1]) > 1e-12) + { break; } } @@ -413,21 +479,22 @@ void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, c delete[] rbeta; } -std::string BetaRadials::trim201(std::string const& str) { +std::string BetaRadials::trim201(std::string const& str) +{ // extract the substring between quotation marks (with whitespace trimmed) // str MUST contain exactly a pair of quotation marks std::string::size_type start = str.find('"'); std::string::size_type end = str.find_last_of('"'); - std::string tmp = str.substr(start+1, end-start-1); + std::string tmp = str.substr(start + 1, end - start - 1); - if (tmp.length() == 0) { + if (tmp.length() == 0) + { return tmp; } start = tmp.find_first_not_of(" \t"); end = tmp.find_last_not_of(" \t"); - return tmp.substr(start, end+1-start); + return tmp.substr(start, end + 1 - start); } - diff --git a/source/module_basis/module_nao/beta_radials.h b/source/module_basis/module_nao/beta_radials.h index bc9e2ba49dc..3056cd85867 100644 --- a/source/module_basis/module_nao/beta_radials.h +++ b/source/module_basis/module_nao/beta_radials.h @@ -20,7 +20,7 @@ * O_beta.build(orb_file, element_index, ofs_log, GlobalV::MY_RANK); * * */ -class BetaRadials: public RadialSet +class BetaRadials : public RadialSet { public: BetaRadials() {} @@ -32,24 +32,23 @@ class BetaRadials: public RadialSet ~BetaRadials() {} //! Build the class from a pseudopotential file - void build(const std::string& file, //!< pseudopotential file name - const int itype = 0, //!< element index in calculation - std::ofstream* ptr_log = nullptr, //!< output file stream for logging - const int rank = 0 //!< MPI rank + void build(const std::string& file, //!< pseudopotential file name + const int itype = 0, //!< element index in calculation + std::ofstream* ptr_log = nullptr, //!< output file stream for logging + const int rank = 0 //!< MPI rank ); private: - //! Read beta projectors from a pseudopotential file of UPF 1.0.0 format void read_beta_upf100(std::ifstream& ifs, //!< input file stream from orbital file - std::ofstream* ptr_log = nullptr, //!< output file stream for logging - const int rank = 0 //!< MPI rank + std::ofstream* ptr_log = nullptr, //!< output file stream for logging + const int rank = 0 //!< MPI rank ); //! Read beta projectors from a pseudopotential file of UPF 2.0.1 format void read_beta_upf201(std::ifstream& ifs, //!< input file stream from orbital file - std::ofstream* ptr_log = nullptr, //!< output file stream for logging - const int rank = 0 //!< MPI rank + std::ofstream* ptr_log = nullptr, //!< output file stream for logging + const int rank = 0 //!< MPI rank ); //! extract the substring between a pair of quotation marks (for UPF v2.0.1) diff --git a/source/module_basis/module_nao/numerical_radial.cpp b/source/module_basis/module_nao/numerical_radial.cpp index c5f71be285b..2b8b90d314f 100644 --- a/source/module_basis/module_nao/numerical_radial.cpp +++ b/source/module_basis/module_nao/numerical_radial.cpp @@ -2,10 +2,6 @@ #include #include -#include -#include -#include -#include #include #include @@ -16,31 +12,25 @@ using ModuleBase::PI; -NumericalRadial::NumericalRadial() +NumericalRadial::NumericalRadial() : + sbt_(new ModuleBase::SphericalBesselTransformer), + use_internal_transformer_(true) { - if (use_internal_transformer_) - { - sbt_ = new ModuleBase::SphericalBesselTransformer; - } } -NumericalRadial::NumericalRadial(const NumericalRadial& other) +NumericalRadial::NumericalRadial(const NumericalRadial& other) : + symbol_(other.symbol_), + itype_(other.itype_), + l_(other.l_), + izeta_(other.izeta_), + nr_(other.nr_), + nk_(other.nk_), + is_fft_compliant_(other.is_fft_compliant_), + pr_(other.pr_), + pk_(other.pk_), + sbt_(other.use_internal_transformer_ ? new ModuleBase::SphericalBesselTransformer : other.sbt_), + use_internal_transformer_(other.use_internal_transformer_) { - symbol_ = other.symbol_; - itype_ = other.itype_; - izeta_ = other.izeta_; - l_ = other.l_; - - nr_ = other.nr_; - nk_ = other.nk_; - - is_fft_compliant_ = other.is_fft_compliant_; - - pr_ = other.pr_; - pk_ = other.pk_; - - use_internal_transformer_ = other.use_internal_transformer_; - // deep copy if (other.ptr_rgrid()) { @@ -57,15 +47,6 @@ NumericalRadial::NumericalRadial(const NumericalRadial& other) std::memcpy(kgrid_, other.kgrid_, nk_ * sizeof(double)); std::memcpy(kvalue_, other.kvalue_, nk_ * sizeof(double)); } - - if (use_internal_transformer_) - { - sbt_ = new ModuleBase::SphericalBesselTransformer; - } - else - { - sbt_ = other.sbt_; - } } NumericalRadial& NumericalRadial::operator=(const NumericalRadial& rhs) @@ -109,22 +90,7 @@ NumericalRadial& NumericalRadial::operator=(const NumericalRadial& rhs) std::memcpy(kvalue_, rhs.kvalue_, nk_ * sizeof(double)); } - if (rhs.use_internal_transformer_) - { - if (!use_internal_transformer_) - { - sbt_ = new ModuleBase::SphericalBesselTransformer; - } - } - else - { - if (use_internal_transformer_) - { - delete sbt_; - } - sbt_ = rhs.sbt_; - } - use_internal_transformer_ = rhs.use_internal_transformer_; + set_transformer(rhs.use_internal_transformer_ ? nullptr : rhs.sbt_, 0); return *this; } @@ -191,28 +157,22 @@ void NumericalRadial::build(const int l, void NumericalRadial::set_transformer(ModuleBase::SphericalBesselTransformer* sbt, int update) { - assert(update == 0 || update == 1 || update == -1); - if (sbt) - { - //! if an external transformer is provided, delete the internal one if any - if (use_internal_transformer_) - { - delete sbt_; - use_internal_transformer_ = false; - } + if (use_internal_transformer_ && sbt) + { // internal -> external + delete sbt_; + use_internal_transformer_ = false; sbt_ = sbt; } - else - { - // if no external transformer is provided, use an internal one - if (!use_internal_transformer_) - { - sbt_ = new ModuleBase::SphericalBesselTransformer; - use_internal_transformer_ = true; - } - // do nothing if an internal one is already in use + else if (!use_internal_transformer_ && !sbt) + { // external -> internal + sbt_ = new ModuleBase::SphericalBesselTransformer; + use_internal_transformer_ = true; + } + else if (!use_internal_transformer_ && sbt) + { // external -> another external + sbt_ = sbt; } switch (update) @@ -308,13 +268,12 @@ void NumericalRadial::set_uniform_grid(const bool for_r_space, } set_grid(for_r_space, ngrid, grid, mode); + delete[] grid; if (enable_fft) { set_uniform_grid(!for_r_space, ngrid, PI / dx, 't', false); } - - delete[] grid; } void NumericalRadial::set_value(const bool for_r_space, const double* const value, const int p) @@ -405,23 +364,19 @@ void NumericalRadial::radtab(const char op, const NumericalRadial& ket, const int l, double* const table, - const bool deriv, - int ntab, - const double* tabgrid) const + const int nr_tab, + const double* const rgrid_tab, + const bool deriv) const { assert(op == 'S' || op == 'I' || op == 'T' || op == 'U'); assert(l >= 0); + assert(rgrid_tab && nr_tab > 0); - // radtab requires that two NumericalRadial objects have exactly the same kgrid_ + // 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_)); - // either tabgrid or rgrid_ exists - assert((tabgrid && ntab > 0) || (!tabgrid && ntab == 0 && rgrid_ && nr_ > 0)); - - ntab = ntab ? ntab : nr_; - tabgrid = tabgrid ? tabgrid : rgrid_; - bool use_radrfft = is_fft_compliant(ntab, tabgrid, nk_, kgrid_); + bool use_radrfft = is_fft_compliant(nr_tab, rgrid_tab, nk_, kgrid_); // function to undergo a spherical Bessel transform: // overlap: chi1(k) * chi2(k) @@ -442,57 +397,21 @@ void NumericalRadial::radtab(const char op, default:; // for overlap integral op_pk = 0 } - if (deriv) - { // derivative of the radial table - if (l == 0) - { // j'_0(x) = -j_1(x) - if (use_radrfft) - { - sbt_->radrfft(1, nk_, kcut(), fk, table, pk_ + ket.pk_ + op_pk - 1); - } - else - { - sbt_->direct(1, nk_, kgrid_, fk, ntab, tabgrid, table, pk_ + ket.pk_ + op_pk - 1); - } - std::for_each(table, table + nr_, [](double& x) { x *= -1; }); - } - else - { // (2*l+1) * j'_l(x) = l * j_{l-1}(x) - (l+1) * j_{l+1}(x) - double* frtmp = new double[ntab]; - if (use_radrfft) - { - sbt_->radrfft(l + 1, nk_, kcut(), fk, table, pk_ + ket.pk_ + op_pk - 1); - sbt_->radrfft(l - 1, nk_, kcut(), fk, frtmp, pk_ + ket.pk_ + op_pk - 1); - } - else - { - sbt_->direct(l + 1, nk_, kgrid_, fk, ntab, tabgrid, table, pk_ + ket.pk_ + op_pk - 1); - sbt_->direct(l - 1, nk_, kgrid_, fk, ntab, tabgrid, frtmp, pk_ + ket.pk_ + op_pk - 1); - } - std::transform(table, table + ntab, frtmp, table, [l](double x, double y) { - return (l * y - (l + 1) * x) / (2 * l + 1); - }); - delete[] frtmp; - } + if (use_radrfft) + { + sbt_->radrfft(l, nk_, kcut(), fk, table, pk_ + ket.pk_ + op_pk, deriv); } else - { // radial table - if (use_radrfft) - { - sbt_->radrfft(l, nk_, kcut(), fk, table, pk_ + ket.pk_ + op_pk); - } - else - { - sbt_->direct(l, nk_, kgrid_, fk, ntab, tabgrid, table, pk_ + ket.pk_ + op_pk); - } + { + sbt_->direct(l, nk_, kgrid_, fk, nr_tab, rgrid_tab, table, pk_ + ket.pk_ + op_pk, deriv); } delete[] fk; - // spherical Bessel transform has a prefactor of sqrt(2/pi) while the prefactor for the radial table + // spherical Bessel transform has a prefactor of sqrt(2/pi) while the prefactor for the radial table // of two-center integrals is 4*pi double pref = ModuleBase::FOUR_PI * std::sqrt(ModuleBase::PI / 2.0); - std::for_each(table, table + ntab, [pref](double& x) { x *= pref; }); + std::for_each(table, table + nr_tab, [pref](double& x) { x *= pref; }); } void NumericalRadial::normalize(bool for_r_space) diff --git a/source/module_basis/module_nao/numerical_radial.h b/source/module_basis/module_nao/numerical_radial.h index 0cbd4d23d47..90435af821f 100644 --- a/source/module_basis/module_nao/numerical_radial.h +++ b/source/module_basis/module_nao/numerical_radial.h @@ -191,21 +191,29 @@ 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 bool deriv = false, //!< [in] if true, "table" would contain the derivative - //!< of the table - int ntab = 0, //!< [in] size of tabgrid - const double* tabgrid = nullptr //!< [in] grid on which the table is calculated. - //!< if nullptr, this->rgrid_ is assumed. + 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 ) const; + //! Normalizes the radial function. + /*! + * The radial function is normalized such that + * + * / +inf 2 + * | dx x f(x) = 1 + * / 0 + * + * where x is r or k. + * */ void normalize(bool for_r_space = true); /*! @@ -321,10 +329,10 @@ class NumericalRadial /*! * @see set_transformer * */ - ModuleBase::SphericalBesselTransformer* sbt_ = nullptr; + ModuleBase::SphericalBesselTransformer* sbt_; //! A flag that marks the ownership of sbt_ - bool use_internal_transformer_ = true; + bool use_internal_transformer_; //! Transforms the r-space values to get k-space values, or vice versa. /*! diff --git a/source/module_basis/module_nao/radial_collection.cpp b/source/module_basis/module_nao/radial_collection.cpp index 138ffe0e2a6..35174de1f6c 100644 --- a/source/module_basis/module_nao/radial_collection.cpp +++ b/source/module_basis/module_nao/radial_collection.cpp @@ -4,17 +4,33 @@ #include "module_basis/module_nao/atomic_radials.h" #include "module_basis/module_nao/beta_radials.h" -RadialCollection::RadialCollection(const RadialCollection& other) +RadialCollection::RadialCollection() : + sbt_(new ModuleBase::SphericalBesselTransformer), + use_internal_transformer_(true) { - ntype_ = other.ntype_; - lmax_ = other.lmax_; - nchi_ = other.nchi_; - nzeta_max_ = other.nzeta_max_; +} + +RadialCollection::RadialCollection(const RadialCollection& other) : + ntype_(other.ntype_), + lmax_(other.lmax_), + nchi_(other.nchi_), + nzeta_max_(other.nzeta_max_), + radset_(nullptr), + iter_(nullptr), + nl_(nullptr), + sbt_(other.use_internal_transformer_ ? new ModuleBase::SphericalBesselTransformer : other.sbt_), + use_internal_transformer_(other.use_internal_transformer_) +{ + if (ntype_ == 0) + { + return; + } radset_ = new RadialSet*[ntype_]; for (int itype = 0; itype < ntype_; ++itype) { radset_[itype] = other.radset_[itype]->clone(); + radset_[itype]->set_transformer(sbt_, 0); } iter_build(); @@ -42,27 +58,35 @@ RadialCollection& RadialCollection::operator=(const RadialCollection& rhs) iter_build(); + set_transformer(rhs.use_internal_transformer_ ? nullptr : rhs.sbt_, 0); + return *this; } -double RadialCollection::rcut_max() const +RadialCollection::~RadialCollection() { - double rmax = 0.0; for (int itype = 0; itype < ntype_; ++itype) { - rmax = std::max(rmax, radset_[itype]->rcut_max()); + delete radset_[itype]; + } + delete[] radset_; + delete[] iter_; // iterator does not control memory; simply delete the pointer array + delete[] nl_; + + if (use_internal_transformer_) + { + delete sbt_; } - return rmax; } -RadialCollection::~RadialCollection() +double RadialCollection::rcut_max() const { + double rmax = 0.0; for (int itype = 0; itype < ntype_; ++itype) { - delete radset_[itype]; + rmax = std::max(rmax, radset_[itype]->rcut_max()); } - delete[] radset_; - delete[] iter_; // iterator does not control memory; simply delete the pointer array + return rmax; } void RadialCollection::cleanup() @@ -77,6 +101,9 @@ void RadialCollection::cleanup() delete[] iter_; // iterator does not control memory; simply delete the pointer array iter_ = nullptr; + delete[] nl_; + nl_ = nullptr; + ntype_ = 0; lmax_ = -1; nchi_ = 0; @@ -88,16 +115,27 @@ void RadialCollection::iter_build() /* * collect the addresses of NumericalRadial objects from different RadialSet objects * so that all NumericalRadial objects can be iterated over in a single loop + * + * objects are sorted by l first, by itype next, by izeta last. * */ delete[] iter_; // iterator does not control memory; simply delete the pointer array + delete[] nl_; + + nl_ = new int[lmax_ + 1]; iter_ = new const NumericalRadial*[nchi_]; + int i = 0; - for (int itype = 0; itype < ntype_; ++itype) + std::fill(nl_, nl_ + lmax_ + 1, 0); + for (int l = 0; l <= lmax_; ++l) { - for (const NumericalRadial* it = radset_[itype]->cbegin(); it != radset_[itype]->cend(); ++it) + for (int itype = 0; itype != ntype_; ++itype) { - iter_[i] = it; - ++i; + for (int izeta = 0; izeta < radset_[itype]->nzeta(l); ++izeta) + { + iter_[i] = &radset_[itype]->chi(l, izeta); + ++i; + ++nl_[l]; + } } } } @@ -137,13 +175,34 @@ void RadialCollection::build(const int nfile, const std::string* const file, con } iter_build(); + + for (int itype = 0; itype < ntype_; ++itype) + { + radset_[itype]->set_transformer(sbt_, 0); + } } void RadialCollection::set_transformer(ModuleBase::SphericalBesselTransformer* const sbt, const int update) { + if (use_internal_transformer_ && sbt) + { // internal -> external + delete sbt_; + use_internal_transformer_ = false; + sbt_ = sbt; + } + else if (!use_internal_transformer_ && !sbt) + { // external -> internal + sbt_ = new ModuleBase::SphericalBesselTransformer; + use_internal_transformer_ = true; + } + else if (!use_internal_transformer_ && sbt) + { // external -> another external + sbt_ = sbt; + } + for (int itype = 0; itype < ntype_; ++itype) { - radset_[itype]->set_transformer(sbt, update); + radset_[itype]->set_transformer(sbt_, update); } } diff --git a/source/module_basis/module_nao/radial_collection.h b/source/module_basis/module_nao/radial_collection.h index db2893fcd53..f1910f6a946 100644 --- a/source/module_basis/module_nao/radial_collection.h +++ b/source/module_basis/module_nao/radial_collection.h @@ -1,6 +1,7 @@ #ifndef RADIAL_COLLECTION_H_ #define RADIAL_COLLECTION_H_ +#include #include #include "module_basis/module_nao/radial_set.h" @@ -16,7 +17,7 @@ class RadialCollection { public: - RadialCollection(){}; + RadialCollection(); RadialCollection(const RadialCollection& other); //!< deep copy RadialCollection& operator=(const RadialCollection& rhs); //!< deep copy @@ -76,24 +77,39 @@ class RadialCollection assert(itype >= 0 && itype < ntype_); return *radset_[itype]; } + //!@} - //! *(this->cbegin()) returns the address of the first NumericalRadial object in the collection + /*! @name Iterators + * + * Enable iteration through all NumericalRadial objects in the collection. + * Objects are sorted by l first, by itype next, by izeta last. + * */ + //!@{ const NumericalRadial** cbegin() const { assert(ntype_ > 0); return iter_; } - //! *(this->cend()-1) returns the address of the last NumericalRadial object in the collection - /*! - * CAVEAT: *(this->cend()) and this->radset_[ntype_-1]->cend() are not equal! - * In general users should never dereference the end iterator. - * */ const NumericalRadial** cend() const { assert(ntype_ > 0); return iter_ + nchi_; } + + //! *(this->cbegin(l)) returns the address of the first NumericalRadial object with angular momentum l + const NumericalRadial** cbegin(const int l) const + { + assert(ntype_ > 0 && l >= 0 && l <= lmax_); + return iter_ + std::accumulate(nl_, nl_ + l, 0); + } + + //! *(this->cbegin(l)) returns the address of one-past last NumericalRadial object with angular momentum l + const NumericalRadial** cend(const int l) const + { + assert(ntype_ > 0 && l >= 0 && l <= lmax_); + return iter_ + std::accumulate(nl_, nl_ + l + 1, 0); + } //!@} /*! @name property setters for all RadialSet objects @@ -119,9 +135,9 @@ class RadialCollection //!@} private: - int ntype_ = 0; //!< number of RadialSet in the collection - int lmax_ = -1; //!< maximum angular momentum of all NumericalRadial objects in the collection - int nchi_ = 0; //!< total number of NumericalRadial objects in the collection + int ntype_ = 0; //!< number of RadialSet in the collection + int lmax_ = -1; //!< maximum angular momentum of all NumericalRadial objects in the collection + int nchi_ = 0; //!< total number of NumericalRadial objects in the collection int nzeta_max_ = 0; //!< maximum number of distinct radial functions given a type & angular momentum //! array of RadialSet objects @@ -145,6 +161,19 @@ class RadialCollection * */ const NumericalRadial** iter_ = nullptr; + //! number of NumericalRadial objects for each angular momentum + int* nl_ = nullptr; + + //! Pointer to the object that provides spherical Bessel transforms + /*! + * All NumericalRadial objects within this class should share the same + * spherical Bessel transformer. + * */ + ModuleBase::SphericalBesselTransformer* sbt_; + + //! A flag that marks the ownership of sbt_ + bool use_internal_transformer_; + //! Deallocates all RadialSet objects and resets all members to default. void cleanup(); diff --git a/source/module_basis/module_nao/radial_set.cpp b/source/module_basis/module_nao/radial_set.cpp index 79084b0dc73..9ef26e77c78 100644 --- a/source/module_basis/module_nao/radial_set.cpp +++ b/source/module_basis/module_nao/radial_set.cpp @@ -1,61 +1,56 @@ #include "module_basis/module_nao/radial_set.h" #include +#include + +#include "module_base/spherical_bessel_transformer.h" + +RadialSet::RadialSet() : + sbt_(new ModuleBase::SphericalBesselTransformer), + use_internal_transformer_(true) +{ +} RadialSet::~RadialSet() { delete[] nzeta_; delete[] chi_; delete[] index_map_; -} -double RadialSet::rcut_max() const -{ - double rmax = 0.0; - for (int i = 0; i != nchi_; ++i) + if (use_internal_transformer_) { - if (chi_[i].rcut() > rmax) - { - rmax = chi_[i].rcut(); - } + delete sbt_; } - return rmax; } -RadialSet::RadialSet(const RadialSet& other) +RadialSet::RadialSet(const RadialSet& other) : + symbol_(other.symbol_), + itype_(other.itype_), + lmax_(other.lmax_), + nzeta_(nullptr), + nzeta_max_(other.nzeta_max_), + nchi_(other.nchi_), + chi_(nullptr), + index_map_(nullptr), + sbt_(other.use_internal_transformer_ ? new ModuleBase::SphericalBesselTransformer : other.sbt_), + use_internal_transformer_(other.use_internal_transformer_) { - symbol_ = other.symbol_; - itype_ = other.itype_; - lmax_ = other.lmax_; - - nzeta_ = nullptr; - if (lmax_ >= 0) + if (nchi_ == 0) { - nzeta_ = new int[lmax_ + 1]; - for (int l = 0; l <= lmax_; l++) - nzeta_[l] = other.nzeta_[l]; + return; } - nzeta_max_ = other.nzeta_max_; - nchi_ = other.nchi_; - chi_ = nullptr; - if (nchi_ > 0) - { - chi_ = new NumericalRadial[nchi_]; - for (int i = 0; i < nchi_; i++) - { - chi_[i] = other.chi_[i]; // deep copy - } - } + nzeta_ = new int[lmax_ + 1]; + std::memcpy(nzeta_, other.nzeta_, (lmax_ + 1) * sizeof(int)); - index_map_ = nullptr; - if (lmax_ >= 0 && nzeta_max_ > 0) + index_map_ = new int[(lmax_ + 1) * nzeta_max_]; + std::memcpy(index_map_, other.index_map_, (lmax_ + 1) * nzeta_max_ * sizeof(int)); + + chi_ = new NumericalRadial[nchi_]; + for (int i = 0; i < nchi_; i++) { - index_map_ = new int[(lmax_ + 1) * nzeta_max_]; - for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++) - { - index_map_[i] = other.index_map_[i]; - } + chi_[i] = other.chi_[i]; // deep copy + chi_[i].set_transformer(sbt_, 0); } } @@ -69,22 +64,24 @@ RadialSet& RadialSet::operator=(const RadialSet& rhs) symbol_ = rhs.symbol_; itype_ = rhs.itype_; lmax_ = rhs.lmax_; - - delete[] nzeta_; - nzeta_ = nullptr; - if (lmax_ >= 0) - { - nzeta_ = new int[lmax_ + 1]; - for (int l = 0; l <= lmax_; l++) - nzeta_[l] = rhs.nzeta_[l]; - } nzeta_max_ = rhs.nzeta_max_; nchi_ = rhs.nchi_; + delete[] nzeta_; delete[] chi_; + delete[] index_map_; + nzeta_ = nullptr; chi_ = nullptr; + index_map_ = nullptr; + if (nchi_ > 0) { + nzeta_ = new int[lmax_ + 1]; + std::memcpy(nzeta_, rhs.nzeta_, (lmax_ + 1) * sizeof(int)); + + index_map_ = new int[(lmax_ + 1) * nzeta_max_]; + std::memcpy(index_map_, rhs.index_map_, (lmax_ + 1) * nzeta_max_ * sizeof(int)); + chi_ = new NumericalRadial[nchi_]; for (int i = 0; i < nchi_; i++) { @@ -92,17 +89,19 @@ RadialSet& RadialSet::operator=(const RadialSet& rhs) } } - delete[] index_map_; - index_map_ = nullptr; - if (lmax_ >= 0 && nzeta_max_ > 0) + set_transformer(rhs.use_internal_transformer_ ? nullptr : rhs.sbt_, 0); + + return *this; +} + +double RadialSet::rcut_max() const +{ + double rmax = 0.0; + for (int i = 0; i < nchi_; ++i) { - index_map_ = new int[(lmax_ + 1) * nzeta_max_]; - for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++) - { - index_map_[i] = rhs.index_map_[i]; - } + rmax = std::max(rmax, chi_[i].rcut()); } - return *this; + return rmax; } int RadialSet::index(const int l, const int izeta) const @@ -129,15 +128,7 @@ void RadialSet::indexing() { for (int izeta = 0; izeta != nzeta_max_; ++izeta) { - if (izeta >= nzeta_[l]) - { - index_map_[l * nzeta_max_ + izeta] = -1; // -1 means no such orbital - } - else - { - index_map_[l * nzeta_max_ + izeta] = index_chi; - ++index_chi; - } + index_map_[l * nzeta_max_ + izeta] = izeta >= nzeta_[l] ? -1 : index_chi++; } } } @@ -151,8 +142,26 @@ const NumericalRadial& RadialSet::chi(const int l, const int izeta) void RadialSet::set_transformer(ModuleBase::SphericalBesselTransformer* const sbt, const int update) { + if (use_internal_transformer_ && sbt) + { // internal -> external + delete sbt_; + use_internal_transformer_ = false; + sbt_ = sbt; + } + else if (!use_internal_transformer_ && !sbt) + { // external -> internal + sbt_ = new ModuleBase::SphericalBesselTransformer; + use_internal_transformer_ = true; + } + else if (!use_internal_transformer_ && sbt) + { // external -> another external + sbt_ = sbt; + } + for (int i = 0; i < nchi_; i++) - chi_[i].set_transformer(sbt, update); + { + chi_[i].set_transformer(sbt_, update); + } } void RadialSet::set_grid(const bool for_r_space, const int ngrid, const double* grid, const char mode) diff --git a/source/module_basis/module_nao/radial_set.h b/source/module_basis/module_nao/radial_set.h index c51f2c03f37..83b69657fed 100644 --- a/source/module_basis/module_nao/radial_set.h +++ b/source/module_basis/module_nao/radial_set.h @@ -5,6 +5,7 @@ #include #include +#include "module_base/spherical_bessel_transformer.h" #include "module_basis/module_nao/numerical_radial.h" //! An abstract class representing a related set of numerical radial functions. @@ -20,7 +21,7 @@ class RadialSet { public: - RadialSet() {} + RadialSet(); RadialSet(const RadialSet&); //!< deep copy RadialSet& operator=(const RadialSet&); //!< deep copy virtual RadialSet* clone() const = 0; //!< for polymorphic copy @@ -45,11 +46,7 @@ class RadialSet int lmax() const { return lmax_; } double rcut_max() const; - int nzeta(const int l) const - { - assert(l >= 0 && l <= lmax_); - return nzeta_[l]; - } + int nzeta(const int l) const { return (l >= 0 && l <= lmax_) ? nzeta_[l] : 0; } int nzeta_max() const { return nzeta_max_; } int nchi() const { return nchi_; } @@ -93,6 +90,16 @@ class RadialSet int* index_map_ = nullptr; //!< map (l,izeta) to an index in chi_ array + //! Pointer to the object that provides spherical Bessel transforms + /*! + * All NumericalRadial objects within this class should share the same + * spherical Bessel transformer. + * */ + ModuleBase::SphericalBesselTransformer* sbt_; + + //! A flag that marks the ownership of sbt_ + bool use_internal_transformer_; + //! deallocates memory and reset all class members to default values void cleanup(); diff --git a/source/module_basis/module_nao/test/CMakeLists.txt b/source/module_basis/module_nao/test/CMakeLists.txt index e791abcd332..b935cbab22b 100644 --- a/source/module_basis/module_nao/test/CMakeLists.txt +++ b/source/module_basis/module_nao/test/CMakeLists.txt @@ -38,6 +38,19 @@ AddTest( LIBS ${math_libs} device base ) +AddTest( + TARGET two_center_table + SOURCES + two_center_table_test.cpp + ../two_center_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 +) + AddTest( TARGET real_gaunt_table SOURCES diff --git a/source/module_basis/module_nao/test/numerical_radial_test.cpp b/source/module_basis/module_nao/test/numerical_radial_test.cpp index c64828277c4..21f2904661c 100644 --- a/source/module_basis/module_nao/test/numerical_radial_test.cpp +++ b/source/module_basis/module_nao/test/numerical_radial_test.cpp @@ -473,28 +473,28 @@ TEST_F(NumericalRadialTest, RadialTable) double* table = new double[sz]; double table_pref = ModuleBase::FOUR_PI * std::sqrt(ModuleBase::PI / 2.0); - chi1.radtab('S', chi2, 0, table); + chi1.radtab('S', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid()); 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.radtab('S', chi2, 2, table, chi1.nr(), chi1.ptr_rgrid()); 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.radtab('T', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid()); 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.radtab('U', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid()); for (int i = 0; i != sz; ++i) { double R = i * dr; @@ -502,28 +502,28 @@ TEST_F(NumericalRadialTest, RadialTable) } // derivative of radial tables - chi1.radtab('S', chi2, 0, table, true); + chi1.radtab('S', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), 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, true); + chi1.radtab('S', chi2, 2, table, chi1.nr(), chi1.ptr_rgrid(), 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, true); + chi1.radtab('T', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), 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, true); + chi1.radtab('U', chi2, 0, table, chi1.nr(), chi1.ptr_rgrid(), true); for (int i = 0; i != sz; ++i) { double R = i * dr; diff --git a/source/module_basis/module_nao/test/radial_collection_test.cpp b/source/module_basis/module_nao/test/radial_collection_test.cpp index 6bb1e63ad7b..e852d4bf851 100644 --- a/source/module_basis/module_nao/test/radial_collection_test.cpp +++ b/source/module_basis/module_nao/test/radial_collection_test.cpp @@ -259,10 +259,16 @@ TEST_F(RadialCollectionTest, Iteration) { orb.build(nfile, file, 'o'); EXPECT_EQ(*orb.cbegin(), &orb(0, 0, 0)); - EXPECT_EQ(*(orb.cbegin() + 5), &orb(1, 0, 0)); - EXPECT_EQ(*(orb.cbegin() + 8), &orb(2, 0, 0)); - EXPECT_EQ(*(orb.cbegin() + 13), &orb(3, 0, 0)); + EXPECT_EQ(*(orb.cbegin() + 2), &orb(1, 0, 0)); + EXPECT_EQ(*(orb.cbegin() + 9), &orb(3, 0, 3)); + EXPECT_EQ(*(orb.cbegin() + 10), &orb(0, 1, 0)); + EXPECT_EQ(*(orb.cbegin() + 17), &orb(0, 2, 0)); + EXPECT_EQ(*(orb.cbegin() + 21), &orb(3, 3, 0)); EXPECT_EQ(*(orb.cend() - 1), &orb(3, 3, 0)); + //EXPECT_EQ(*(orb.cbegin() + 5), &orb(1, 0, 0)); + //EXPECT_EQ(*(orb.cbegin() + 8), &orb(2, 0, 0)); + //EXPECT_EQ(*(orb.cbegin() + 13), &orb(3, 0, 0)); + //EXPECT_EQ(*(orb.cend() - 1), &orb(3, 3, 0)); } int main(int argc, char** argv) diff --git a/source/module_basis/module_nao/test/two_center_table_test.cpp b/source/module_basis/module_nao/test/two_center_table_test.cpp new file mode 100644 index 00000000000..6b6b461c201 --- /dev/null +++ b/source/module_basis/module_nao/test/two_center_table_test.cpp @@ -0,0 +1,342 @@ +#include "module_basis/module_nao/two_center_table.h" + +#include "gtest/gtest.h" +#include "module_base/constants.h" +#include "module_base/spherical_bessel_transformer.h" + +#include "module_base/math_integral.h" + +#include +using iclock = std::chrono::high_resolution_clock; + +////////////////////////////////////////// +//! for comparison with module_ao +#include "module_basis/module_ao/ORB_table_phi.h" +////////////////////////////////////////// + +#ifdef __MPI +#include +#endif + +/*********************************************************** + * Unit test of class "TwoCenterTable" + ***********************************************************/ +/*! + * Tested functions: + * + * - build + * - builds a two-center integral radial table from two RadialCollection objects + * */ +class TwoCenterTableTest : public ::testing::Test +{ + protected: + void SetUp(); + void TearDown(); + + TwoCenterTable S_tab; + TwoCenterTable S_dtab; + TwoCenterTable T_tab; + TwoCenterTable T_dtab; + + TwoCenterTable betapsi_tab; + TwoCenterTable betapsi_dtab; + + RadialCollection orb; + int nfile = 0; //! number of orbital files + std::string* file = nullptr; //!< orbital files to read from + std::string log_file = "./test_files/two_center_table.log"; //!< file for logging + + double table_tol = 1e-6; //! tolerance for comparison between new and legacy tables +}; + +void TwoCenterTableTest::SetUp() +{ +#ifdef __MPI + MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK); +#endif + + std::string dir = "../../../../../tests/PP_ORB/"; + nfile = 8; + file = new std::string[nfile]; + file[0] = dir + "C_gga_8au_100Ry_2s2p1d.orb"; + file[1] = dir + "Fe_gga_9au_100Ry_4s2p2d1f.orb"; + file[2] = dir + "O_gga_10au_100Ry_2s2p1d.orb"; + file[3] = dir + "H_gga_8au_60Ry_2s1p.orb"; + file[4] = dir + "Cs_gga_10au_100Ry_4s2p1d.orb"; + file[5] = dir + "Pb_gga_7au_100Ry_2s2p2d1f.orb"; + file[6] = dir + "F_gga_7au_100Ry_2s2p1d.orb"; + file[7] = dir + "I_gga_7au_100Ry_2s2p2d1f.orb"; +} + +void TwoCenterTableTest::TearDown() +{ + delete[] file; +} + +TEST_F(TwoCenterTableTest, BuildOverlapAndKinetic) +{ + orb.build(nfile, file, 'o'); + double rmax = orb.rcut_max() * 2.0; + double dr = 0.01; + int nr = static_cast(rmax / dr) + 1; + + //ModuleBase::SphericalBesselTransformer sbt; + //sbt.set_fftw_plan_flag(FFTW_MEASURE); // not necessarily worth it! + //orb.set_transformer(&sbt, 0); + orb.set_uniform_grid(true, nr, rmax, 'i', true); + + const double* rgrid = orb(0,0,0).ptr_rgrid(); + + iclock::time_point start; + std::chrono::duration dur; + + start = iclock::now(); + S_tab.build(orb, orb, 'S', nr, rgrid, false); + T_tab.build(orb, orb, 'T', nr, rgrid, false); + S_dtab.build(orb, orb, 'S', nr, rgrid, true); + T_dtab.build(orb, orb, 'T', nr, rgrid, true); + dur = iclock::now() - start; + std::cout << "time elapsed = " << dur.count() << " s" << std::endl; + + // check whether analytical derivative and finite difference agree + for (int i = 0; i != S_tab.ntab(); ++i) + { + const double* f = S_tab.ptr_table(0, 0, 0, 0, 0, 0, 0) + i * S_tab.nr(); + const double* df = S_dtab.ptr_table(0, 0, 0, 0, 0, 0, 0) + i * S_dtab.nr(); + for (int ir = 5; ir != S_tab.nr() - 4; ++ir) + { + double df_fd = ( +1.0/280*f[ir-4] - 4.0/105*f[ir-3] + +1.0/5*f[ir-2] - 4.0/5*f[ir-1] + -1.0/280*f[ir+4] + 4.0/105*f[ir+3] + -1.0/5*f[ir+2] + 4.0/5*f[ir+1] + ) / 0.01; + + EXPECT_NEAR(df_fd, df[ir], 1e-5); + } + } +} + +TEST_F(TwoCenterTableTest, LegacyConsistency) +{ + // use less files so that the test wouldn't take too long + nfile = 3; + + orb.build(nfile, file, 'o'); + double rmax = orb.rcut_max() * 2.0; + double dr = 0.01; + int nr = static_cast(rmax / dr) + 1; + + //ModuleBase::SphericalBesselTransformer sbt; + //orb.set_transformer(&sbt, 0); + orb.set_uniform_grid(true, nr, rmax, 'i', true); + + const double* rgrid = orb(0,0,0).ptr_rgrid(); + + iclock::time_point start; + std::chrono::duration dur; + + start = iclock::now(); + S_tab.build(orb, orb, 'S', nr, rgrid, false); + T_tab.build(orb, orb, 'T', nr, rgrid, false); + S_dtab.build(orb, orb, 'S', nr, rgrid, true); + T_dtab.build(orb, orb, 'T', nr, rgrid, true); + dur = iclock::now() - start; + std::cout << "radfft time elapsed = " << dur.count() << " s" << std::endl; + + ////////////////////////////////////////////// + // module_ao ORB_table_phi + ////////////////////////////////////////////// + + ORB_table_phi otp; + LCAO_Orbitals lcao; + ModuleBase::Sph_Bessel_Recursive::D2 sbr; + + // prepare data required to initialized ORB_table_phi + std::ofstream ofs_log; + int ntype; + int lmax; + int out_mat_r; // unused variable + bool force_flag; + int my_rank; + bool deepks_setorb; + bool read_in_flag; + std::string descriptor_file; + std::vector orbital_file; + double ecutwfc; + double dk; + double dR; + double Rmax; + + ofs_log.open("ORB_table_phi_test.log"); + ntype = nfile; + lmax = 3; + out_mat_r = 0; // unused variable + force_flag = true; + my_rank = 0; + deepks_setorb = false; + + read_in_flag = true; + for (int i = 0; i != nfile; ++i) { + orbital_file.push_back(file[i]); + } + + // below we mimic ORB_control::read_orb_first + ecutwfc = 10000.0; // ecutwfc has to be very large in order to reach the agreement between new & old tables + dk = 0.01; + dR = 0.01; + Rmax = 30; + + lcao.read_in_flag = read_in_flag; + lcao.orbital_file = orbital_file; +#ifdef __MPI + lcao.bcast_files(ntype, GlobalV::MY_RANK); +#endif + + // see ORB_control::read_orb_first + lcao.ecutwfc = ecutwfc; + lcao.dk = dk; + lcao.dR = dR; + lcao.Rmax = Rmax; + + lcao.Read_Orbitals(ofs_log, ntype, lmax, deepks_setorb, out_mat_r, + force_flag, my_rank); + + otp.allocate(ntype, lmax, lcao.get_kmesh(), Rmax, dR, dk); + + auto calc_nr = [](double const& Rmax, double const& dR) { + int nr = static_cast(Rmax / dR) + 4; + return (nr%2) ? nr : nr+1; + }; + + // initialize spherical bessel table + sbr.set_dx(dR*dk); + + // max(l+l) = 4, but the derivative calculation of j_l relies on j_{l+1} + sbr.cal_jlx(2*lmax+1, calc_nr(Rmax, dR), lcao.get_kmesh()); + + otp.pSB = &sbr; + otp.init_OV_Tpair(lcao); + otp.init_OV_Opair(lcao); + + start = iclock::now(); + otp.init_Table(lcao); + dur = iclock::now() - start; + std::cout << "legacy time elapsed = " << dur.count() << " s" << std::endl; + + //double* tmp = new double[2001]; + //const NumericalRadial& f = orb(1,3,0); + //f.radtab('T', f, 6, tmp); + //for (int ir = 0; ir != 21; ++ir) { + // std::cout << tmp[ir] << std::endl; + //} + //std::cout << std::endl; + +// double* ktmp = new double[f.nk()]; +// for (int ik = 0; ik != f.nk(); ++ik) { +// ktmp[ik] = std::pow(f.ptr_kvalue()[ik] * f.ptr_kgrid()[ik] * f.ptr_kgrid()[ik],2); +// } +// double r1 = f.ptr_rgrid()[1]; + + + std::cout << std::endl; + + // compare the table + // length of table depends on the cutoff radius + //for (int ir = 0; ir != 1600; ++ir) { + // double err = std::abs(T_tab.ptr_table(1,3,0,1,3,0,6)[ir] - otp.Table_TR[0][3][80][6][ir]); + // if (err > table_tol) { + // std::cout << " ir = " << ir << std::endl; + // } + // EXPECT_NEAR(T_tab.ptr_table(1,3,0,1,3,0,6)[ir], otp.Table_TR[0][3][80][6][ir], table_tol); + // //EXPECT_NEAR(S_tab.ptr_table(0,0,0,0,0,0,0)[ir], otp.Table_SR[0][0][0][0][ir], table_tol); + // //EXPECT_NEAR(T_dtab.ptr_table(0,0,0,0,0,0,0)[ir], otp.Table_TR[1][0][0][0][ir], table_tol); + // //EXPECT_NEAR(S_dtab.ptr_table(0,0,0,0,0,0,0)[ir], otp.Table_SR[1][0][0][0][ir], table_tol); + //} + + //std::cout << T_tab.nr() << std::endl; + //std::cout << otp.OV_Tpair(1,1) << std::endl; + //std::cout << otp.OV_Opair(otp.OV_Tpair(1,1), 3, 3, 0, 0) << std::endl; + + //////////////////////////////////////////// + for (int T1 = 0; T1 < ntype ; T1++) + { + for (int T2 = T1 ; T2 < ntype ; T2++) + { + int Tpair=otp.OV_Tpair(T1,T2); + + for (int L1 = 0; L1 <= lcao.Phi[T1].getLmax(); L1++) + { + for (int N1 = 0; N1 < lcao.Phi[T1].getNchi(L1); N1++) + { + for (int L2 = 0; L2 <= lcao.Phi[T2].getLmax(); L2 ++) + { + for (int N2 = 0; N2 < lcao.Phi[T2].getNchi(L2); N2++) + { + int Opair = otp.OV_Opair(Tpair,L1,L2,N1,N2); + + for (int L = std::abs(L1-L2); L <= (L1+L2) ; L += 2) + { + //if (T1 == 0 && L1 == 1 && N1 == 1 && T2 == 0 && L2 == 1 && N2 == 1 && L == 0) { + // for (int ir = 1; ir != 16; ++ir) { + // printf("%12.8f %12.8f\n", + // T_dtab.ptr_table(T1, L1, N1, T2, L2, N2, L)[ir], + // otp.Table_TR[1][Tpair][Opair][L][ir]); + // } + //} + //break; + double err_max = -1.0; + for (int ir = 1; ir != 1600; ++ir) { + //double err = std::abs(T_dtab.ptr_table(T1, L1, N1, T2, L2, N2, L)[ir] + // - otp.Table_TR[1][Tpair][Opair][L][ir]); + //if ( err > table_tol) { + // if (err_max < 0) { + // printf("%i %i %i %i %i %i %i %5i", T1, L1, N1, T2, L2, N2, L, ir); + // err_max = err; + // } else { + // err_max = std::max(err_max, err); + // } + //} + + // table whose integrand has a higher order of k is less accurate + // as it requires a larger ecutwfc to converge. + EXPECT_NEAR(T_tab.ptr_table(T1, L1, N1, T2, L2, N2, L)[ir], + otp.Table_TR[0][Tpair][Opair][L][ir], table_tol * 10); + EXPECT_NEAR(S_tab.ptr_table(T1, L1, N1, T2, L2, N2, L)[ir], + otp.Table_SR[0][Tpair][Opair][L][ir], table_tol); + EXPECT_NEAR(T_dtab.ptr_table(T1, L1, N1, T2, L2, N2, L)[ir], + otp.Table_TR[1][Tpair][Opair][L][ir], table_tol * 100); + EXPECT_NEAR(S_dtab.ptr_table(T1, L1, N1, T2, L2, N2, L)[ir], + otp.Table_SR[1][Tpair][Opair][L][ir], table_tol); + } + if (err_max > 0) { + printf(" %8.5e\n", err_max); + } + } + } + } + } + } + } + } + + + ////////////////////////////////////////////// + otp.Destroy_Table(lcao); +} + +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; +} diff --git a/source/module_basis/module_nao/two_center_table.cpp b/source/module_basis/module_nao/two_center_table.cpp new file mode 100644 index 00000000000..762defedfc8 --- /dev/null +++ b/source/module_basis/module_nao/two_center_table.cpp @@ -0,0 +1,123 @@ +#include "module_basis/module_nao/two_center_table.h" + +#include +#include +#include + +TwoCenterTable::~TwoCenterTable() { delete[] rgrid_; } + +void TwoCenterTable::build(const RadialCollection& bra, + const RadialCollection& ket, + const char op, + const int nr, + const double* const rgrid, + const bool deriv) +{ + assert(nr > 0 && rgrid); + cleanup(); + + op_ = op; + is_deriv_ = deriv; + + nr_ = nr; + rgrid_ = new double[nr_]; + std::memcpy(rgrid_, rgrid, nr_ * sizeof(double)); + + // index the table by generating a map from (itype1, l1, izeta1, itype2, l2, izeta2, l) to a row index + index_map_.resize({bra.ntype(), bra.lmax() + 1, bra.nzeta_max(), + ket.ntype(), ket.lmax() + 1, ket.nzeta_max(), bra.lmax() + ket.lmax() + 1}); + std::fill(index_map_.data(), index_map_.data() + index_map_.NumElements(), -1); + + ntab_ = 0; + + for (int l = 0; l <= bra.lmax() + ket.lmax(); ++l) + { + for (int l1 = 0; l1 <= bra.lmax(); ++l1) + { + for (const NumericalRadial** it1 = bra.cbegin(l1); it1 != bra.cend(l1); ++it1) + { + for (int l2 = std::abs(l1 - l); l2 <= std::min(ket.lmax(), l + l1); l2 += 2) + { + for (const NumericalRadial** it2 = ket.cbegin(l2); it2 != ket.cend(l2); ++it2) + { + table_index(*it1, *it2, l) = ntab_; + ++ntab_; + } + } + } + } + } + + // irow is now the number of rows in the table + table_.resize({ntab_, nr_}); + + for (int l = 0; l <= bra.lmax() + ket.lmax(); ++l) + { + for (int l1 = 0; l1 <= bra.lmax(); ++l1) + { + for (const NumericalRadial** it1 = bra.cbegin(l1); it1 != bra.cend(l1); ++it1) + { + for (int l2 = std::abs(l1 - l); l2 <= std::min(ket.lmax(), l + l1); l2 += 2) + { + for (const NumericalRadial** it2 = ket.cbegin(l2); it2 != ket.cend(l2); ++it2) + { + (*it1)->radtab(op, + **it2, + l, + table_.inner_most_ptr(table_index(*it1, *it2, l)), + nr_, + rgrid_, + deriv); + } + } + } + } + } +} + +const double* TwoCenterTable::ptr_table(const int itype1, + const int l1, + const int izeta1, + const int itype2, + const int l2, + const int izeta2, + const int l) const +{ + assert(is_present(itype1, l1, izeta1, itype2, l2, izeta2, l)); + return table_.inner_most_ptr(index_map_.get_value(itype1, l1, izeta1, itype2, l2, izeta2, l)); +} + +int& TwoCenterTable::table_index(const NumericalRadial* it1, const NumericalRadial* it2, const int l) +{ + return index_map_.get_value(it1->itype(), it1->l(), it1->izeta(), it2->itype(), it2->l(), it2->izeta(), l); +} + +void TwoCenterTable::cleanup() +{ + op_ = '\0'; + is_deriv_ = false; + nr_ = 0; + + delete[] rgrid_; + rgrid_ = nullptr; + + table_.resize({0}); + index_map_.resize({0}); +} + +bool TwoCenterTable::is_present(const int itype1, + const int l1, + const int izeta1, + const int itype2, + const int l2, + const int izeta2, + const int l) const +{ + // The given indices map to an entry in the table if they fall within the bounds of index_map_ and + // the value of the entry in index_map_ is non-negative + return itype1 >= 0 && itype1 < index_map_.shape().dim_size(0) && l1 >= 0 && l1 < index_map_.shape().dim_size(1) + && izeta1 >= 0 && izeta1 < index_map_.shape().dim_size(2) && itype2 >= 0 + && itype2 < index_map_.shape().dim_size(3) && l2 >= 0 && l2 < index_map_.shape().dim_size(4) && izeta2 >= 0 + && izeta2 < index_map_.shape().dim_size(5) && l >= 0 && l <= index_map_.shape().dim_size(6) + && index_map_.get_value(itype1, l1, izeta1, itype2, l2, izeta2, l) >= 0; +} diff --git a/source/module_basis/module_nao/two_center_table.h b/source/module_basis/module_nao/two_center_table.h new file mode 100644 index 00000000000..9aed451d4e6 --- /dev/null +++ b/source/module_basis/module_nao/two_center_table.h @@ -0,0 +1,86 @@ +#ifndef TWO_CENTER_TABLE_H +#define TWO_CENTER_TABLE_H + +#include "module_base/module_container/tensor.h" +#include "module_base/spherical_bessel_transformer.h" +#include "module_basis/module_nao/radial_collection.h" + +class TwoCenterTable +{ + public: + TwoCenterTable() {} + ~TwoCenterTable(); + + void build(const RadialCollection& bra, //!< [in] radial collection involved in + const RadialCollection& ket, //!< [in] radial collection involved in + const char op, //!< [in] operator of the two-center integral + const int nr, //!< [in] number of table grid points + const double* const rgrid, //!< [in] table grid + const bool deriv = false //!< [in] if the target is a derivative table + ); + + /*! + * @name Getters + * */ + //!@{ + //! returns true if the table stores the derivative of the radial table + bool is_deriv() const { return is_deriv_; } + + //! returns the operator of the two-center integral + char op() const { return op_; } + + //! returns the number of radial points of each table + int nr() const { return nr_; } + + // returns the number of table entries + int ntab() const { return ntab_; } + + //! returns the radius cutoff of the table + double rmax() const { return rgrid_ ? rgrid_[nr_ - 1] : 0.0; } + + //! gets the pointer to the table's grid points (read-only) + const double* ptr_rgrid() const { return rgrid_; } + + //! gets the read-only pointer to a specific table entry + const double* ptr_table(const int itype1, //!< [in] element index of chi1 + const int l1, //!< [in] angular momentum of chi1 + const int izeta1, //!< [in] zeta number of chi1 + const int itype2, //!< [in] element index of chi2 + const int l2, //!< [in] angular momentum of chi2 + const int izeta2, //!< [in] zeta number of chi2 + const int l //!< [in] angular momentum of the entry + ) const; + //!@} + + private: + char op_ = '\0'; //!< operator of the two-center integral + bool is_deriv_ = false; //!< if true, table_ stores the derivative of the radial table + + int nr_ = 0; //!< number of radial points of each table + double* rgrid_ = nullptr; //!< radial grid of each table + + int ntab_ = 0; //!< number of table entries + + //! two-center integral radial table, stored as a row-major matrix + container::Tensor table_{container::DataType::DT_DOUBLE, container::TensorShape({0})}; + + //! map (itype1, l1, izeta1, itype2, l2, izeta2, l) to a row index in the table + container::Tensor index_map_{container::DataType::DT_INT, container::TensorShape({0})}; + + //! returns the row-index of the table corresponding to the given two radial functions and l + int& table_index(const NumericalRadial* ptr_rad1, const NumericalRadial* ptr_rad2, const int l); + + //! deallocates memory and reset variables to default. + void cleanup(); + + //! returns whether the given indices map to an entry in the table + bool is_present(const int itype1, + const int l1, + const int izeta1, + const int itype2, + const int l2, + const int izeta2, + const int l) const; +}; + +#endif