diff --git a/source/source_lcao/module_rt/CMakeLists.txt b/source/source_lcao/module_rt/CMakeLists.txt index 22af0621701..8b27378a52f 100644 --- a/source/source_lcao/module_rt/CMakeLists.txt +++ b/source/source_lcao/module_rt/CMakeLists.txt @@ -12,6 +12,7 @@ if(ENABLE_LCAO) upsi.cpp td_info.cpp velocity_op.cpp + snap_projector_half_tddft.cpp snap_psibeta_half_tddft.cpp td_folding.cpp solve_propagation.cpp diff --git a/source/source_lcao/module_rt/snap_projector_half_tddft.cpp b/source/source_lcao/module_rt/snap_projector_half_tddft.cpp new file mode 100644 index 00000000000..a2c48584e66 --- /dev/null +++ b/source/source_lcao/module_rt/snap_projector_half_tddft.cpp @@ -0,0 +1,397 @@ +#include "snap_projector_half_tddft.h" + +#include "source_base/constants.h" +#include "source_base/global_function.h" +#include "source_base/math_integral.h" +#include "source_base/math_lebedev_laikov.h" +#include "source_base/math_polyint.h" +#include "source_base/timer.h" +#include "source_base/ylm.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace module_rt +{ +namespace +{ +constexpr int default_radial_grid_num = 140; +constexpr int default_lebedev_grid_points = 110; + +/** + * @brief Cached Gauss-Legendre radial grid for a requested grid size. + */ +struct GaussLegendreGrid +{ + explicit GaussLegendreGrid(const int ngrid) : x(ngrid), w(ngrid) + { + ModuleBase::Integral::Gauss_Legendre_grid_and_weight(ngrid, x.data(), w.data()); + } + + std::vector x; + std::vector w; +}; + +const GaussLegendreGrid& gauss_legendre_grid(const int ngrid) +{ + // Tests may request non-default radial grids, so cache by grid size. + static std::map> cache; + static std::mutex cache_mutex; + + std::lock_guard lock(cache_mutex); + std::shared_ptr& grid = cache[ngrid]; + if (!grid) + { + grid.reset(new GaussLegendreGrid(ngrid)); + } + return *grid; +} + +/** + * @brief Owned Lebedev-Laikov angular grid generated at runtime. + */ +struct AngularGridData +{ + explicit AngularGridData(const int ngrid) : x(ngrid), y(ngrid), z(ngrid), w(ngrid) + { + ModuleBase::Lebedev_laikov_grid grid(ngrid); + grid.generate_grid_points(); + const ModuleBase::Vector3* grid_coor = grid.get_grid_coor(); + const double* weight = grid.get_weight(); + for (int i = 0; i < ngrid; ++i) + { + x[i] = grid_coor[i].x; + y[i] = grid_coor[i].y; + z[i] = grid_coor[i].z; + w[i] = weight[i]; + } + } + + std::vector x; + std::vector y; + std::vector z; + std::vector w; +}; + +/** + * @brief Non-owning view used by the integration loops. + */ +struct AngularGridView +{ + int size = 0; + const double* x = nullptr; + const double* y = nullptr; + const double* z = nullptr; + const double* w = nullptr; +}; + +bool is_supported_lebedev_grid(const int ngrid) +{ + static const std::set supported_grids + = {6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434, + 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810}; + return supported_grids.find(ngrid) != supported_grids.end(); +} + +AngularGridView angular_grid(const int ngrid) +{ + if (!is_supported_lebedev_grid(ngrid)) + { + ModuleBase::WARNING_QUIT("snap_projector_half_tddft", + "Unsupported Lebedev-Laikov grid size: " + std::to_string(ngrid)); + } + + if (ngrid == default_lebedev_grid_points) + { + // Keep the production path on the historical static 110-point table. + AngularGridView view; + view.size = default_lebedev_grid_points; + view.x = ModuleBase::Integral::Lebedev_Laikov_grid110_x; + view.y = ModuleBase::Integral::Lebedev_Laikov_grid110_y; + view.z = ModuleBase::Integral::Lebedev_Laikov_grid110_z; + view.w = ModuleBase::Integral::Lebedev_Laikov_grid110_w; + return view; + } + + // Higher-order grids are generated lazily for tests and future callers. + static std::map> cache; + static std::mutex cache_mutex; + + std::lock_guard lock(cache_mutex); + std::shared_ptr& data = cache[ngrid]; + if (!data) + { + data.reset(new AngularGridData(ngrid)); + } + + AngularGridView view; + view.size = ngrid; + view.x = data->x.data(); + view.y = data->y.data(); + view.z = data->z.data(); + view.w = data->w.data(); + return view; +} + +double radial_factor(const ProjectorChannel& channel, const double r, const double w_radial) +{ + const double projector_val + = ModuleBase::PolyInt::Polynomial_Interpolation(channel.radial_values, channel.mesh, channel.dk, r); + + return projector_val * r * w_radial; +} +} // namespace + +void snap_projector_half_tddft(const LCAO_Orbitals& orb, + const std::vector& projector_channels, + std::vector>>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R0, + const ModuleBase::Vector3& A, + const bool& calc_r, + const char* timer_name) +{ + // Preserve the production default while allowing tests to call the overload. + SnapIntegrationOptions options; + options.radial_grid_num = default_radial_grid_num; + options.lebedev_grid_points = default_lebedev_grid_points; + snap_projector_half_tddft(orb, projector_channels, nlm, R1, T1, L1, m1, N1, R0, A, calc_r, options, timer_name); +} + +void snap_projector_half_tddft(const LCAO_Orbitals& orb, + const std::vector& projector_channels, + std::vector>>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R0, + const ModuleBase::Vector3& A, + const bool& calc_r, + const SnapIntegrationOptions& options, + const char* timer_name) +{ + ModuleBase::timer::start("module_rt", timer_name); + if (options.radial_grid_num <= 0) + { + ModuleBase::WARNING_QUIT("snap_projector_half_tddft", "The radial grid size must be positive."); + } + const int radial_grid_num = options.radial_grid_num; + const AngularGridView lebedev = angular_grid(options.lebedev_grid_points); + + const int required_size = calc_r ? 4 : 1; + if (nlm.size() != required_size) + { + nlm.resize(required_size); + } + + int natomwfc = 0; + std::vector active(projector_channels.size(), false); + + const double Rcut1 = orb.Phi[T1].getRcut(); + const ModuleBase::Vector3 dRa = R0 - R1; + const double distance10 = dRa.norm(); + + bool any_active = false; + for (int ich = 0; ich < static_cast(projector_channels.size()); ++ich) + { + const ProjectorChannel& channel = projector_channels[ich]; + natomwfc += 2 * channel.l + 1; + if (distance10 <= Rcut1 + channel.rcut) + { + active[ich] = true; + any_active = true; + } + } + + for (auto& x: nlm) + { + x.assign(natomwfc, 0.0); + } + + if (natomwfc == 0 || !any_active) + { + ModuleBase::timer::end("module_rt", timer_name); + return; + } + + // The LCAO orbital is sampled at r + R0 - R1 around the projector center. + const auto& phi_ln = orb.Phi[T1].PhiLN(L1, N1); + const int mesh_r1 = phi_ln.getNr(); + const double* psi_1 = phi_ln.getPsi(); + const double dk_1 = phi_ln.getDk(); + + const GaussLegendreGrid& gl = gauss_legendre_grid(radial_grid_num); + std::vector r_radial(radial_grid_num); + std::vector w_radial(radial_grid_num); + + std::vector A_dot_lebedev(lebedev.size); + for (int ian = 0; ian < lebedev.size; ++ian) + { + A_dot_lebedev[ian] = A.x * lebedev.x[ian] + A.y * lebedev.y[ian] + A.z * lebedev.z[ian]; + } + + std::vector> result_angular; + std::vector> res_ang_x; + std::vector> res_ang_y; + std::vector> res_ang_z; + std::vector rly1((L1 + 1) * (L1 + 1)); + std::vector> rly0_cache(lebedev.size); + + int index_offset = 0; + for (int ich = 0; ich < static_cast(projector_channels.size()); ++ich) + { + const ProjectorChannel& channel = projector_channels[ich]; + const int L0 = channel.l; + const int num_m0 = 2 * L0 + 1; + + if (!active[ich]) + { + index_offset += num_m0; + continue; + } + + assert(channel.mesh > 0); + assert(channel.radial_values != nullptr); + assert(channel.radial_grid != nullptr); + + const double r_min = channel.radial_grid[0]; + const double r_max = channel.radial_grid[channel.mesh - 1]; + const double xl = (r_max - r_min) * 0.5; + const double xmean = (r_max + r_min) * 0.5; + + for (int i = 0; i < radial_grid_num; ++i) + { + r_radial[i] = xmean + xl * gl.x[i]; + w_radial[i] = xl * gl.w[i]; + } + + const double A_phase = A * R0; + const std::complex exp_iAR0 = std::exp(ModuleBase::IMAG_UNIT * A_phase); + + // Y_lm(projector direction) only depends on the angular grid. + for (int ian = 0; ian < lebedev.size; ++ian) + { + ModuleBase::Ylm::rl_sph_harm(L0, lebedev.x[ian], lebedev.y[ian], lebedev.z[ian], rly0_cache[ian]); + } + + if (result_angular.size() < static_cast(num_m0)) + { + result_angular.resize(num_m0); + if (calc_r) + { + res_ang_x.resize(num_m0); + res_ang_y.resize(num_m0); + res_ang_z.resize(num_m0); + } + } + + for (int ir = 0; ir < radial_grid_num; ++ir) + { + const double r_val = r_radial[ir]; + + std::fill(result_angular.begin(), result_angular.begin() + num_m0, 0.0); + if (calc_r) + { + std::fill(res_ang_x.begin(), res_ang_x.begin() + num_m0, 0.0); + std::fill(res_ang_y.begin(), res_ang_y.begin() + num_m0, 0.0); + std::fill(res_ang_z.begin(), res_ang_z.begin() + num_m0, 0.0); + } + + for (int ian = 0; ian < lebedev.size; ++ian) + { + const double x = lebedev.x[ian]; + const double y = lebedev.y[ian]; + const double z = lebedev.z[ian]; + const double w_ang = lebedev.w[ian]; + + const double rx = r_val * x; + const double ry = r_val * y; + const double rz = r_val * z; + + const double tx = rx + dRa.x; + const double ty = ry + dRa.y; + const double tz = rz + dRa.z; + const double tnorm = std::sqrt(tx * tx + ty * ty + tz * tz); + + if (tnorm > Rcut1) + { + continue; + } + + if (tnorm > 1e-10) + { + const double inv_tnorm = 1.0 / tnorm; + ModuleBase::Ylm::rl_sph_harm(L1, tx * inv_tnorm, ty * inv_tnorm, tz * inv_tnorm, rly1); + } + else + { + ModuleBase::Ylm::rl_sph_harm(L1, 0.0, 0.0, 1.0, rly1); + } + + const double phase = r_val * A_dot_lebedev[ian]; + const std::complex exp_iAr = std::exp(ModuleBase::IMAG_UNIT * phase); + const double interp_psi = ModuleBase::PolyInt::Polynomial_Interpolation(psi_1, mesh_r1, dk_1, tnorm); + const double ylm_L1_val = rly1[L1 * L1 + m1]; + const std::complex common_factor = exp_iAr * ylm_L1_val * interp_psi * w_ang; + + // Accumulate all magnetic components of the same projector channel. + const std::vector& rly0_vec = rly0_cache[ian]; + const int offset_L0 = L0 * L0; + for (int m0 = 0; m0 < num_m0; ++m0) + { + const std::complex term = common_factor * rly0_vec[offset_L0 + m0]; + result_angular[m0] += term; + + if (calc_r) + { + res_ang_x[m0] += term * (rx + R0.x); + res_ang_y[m0] += term * (ry + R0.y); + res_ang_z[m0] += term * (rz + R0.z); + } + } + } + + const double factor = radial_factor(channel, r_val, w_radial[ir]); + int current_idx = index_offset; + for (int m0 = 0; m0 < num_m0; ++m0) + { + nlm[0][current_idx] += factor * result_angular[m0] * exp_iAR0; + if (calc_r) + { + nlm[1][current_idx] += factor * res_ang_x[m0] * exp_iAR0; + nlm[2][current_idx] += factor * res_ang_y[m0] * exp_iAR0; + nlm[3][current_idx] += factor * res_ang_z[m0] * exp_iAR0; + } + ++current_idx; + } + } + + index_offset += num_m0; + } + + for (auto& dim: nlm) + { + for (auto& x: dim) + { + x = std::conj(x); + } + } + + assert(index_offset == natomwfc); + ModuleBase::timer::end("module_rt", timer_name); +} + +} // namespace module_rt diff --git a/source/source_lcao/module_rt/snap_projector_half_tddft.h b/source/source_lcao/module_rt/snap_projector_half_tddft.h new file mode 100644 index 00000000000..15baf98854e --- /dev/null +++ b/source/source_lcao/module_rt/snap_projector_half_tddft.h @@ -0,0 +1,74 @@ +#ifndef SNAP_PROJECTOR_HALF_TDDFT_H +#define SNAP_PROJECTOR_HALF_TDDFT_H + +#include "source_base/vector3.h" +#include "source_basis/module_ao/ORB_read.h" + +#include +#include + +namespace module_rt +{ + +/** + * @brief Radial projector channel integrated against one LCAO orbital. + */ +struct ProjectorChannel +{ + int l = 0; + int mesh = 0; + double dk = 0.0; + double rcut = 0.0; + const double* radial_values = nullptr; + const double* radial_grid = nullptr; +}; + +/** + * @brief Numerical quadrature settings for projector snapshots. + * + * The default values reproduce the production RT-TDDFT path. + */ +struct SnapIntegrationOptions +{ + int radial_grid_num = 140; + int lebedev_grid_points = 110; +}; + +/** + * @brief Compute with default quadrature settings. + */ +void snap_projector_half_tddft(const LCAO_Orbitals& orb, + const std::vector& projector_channels, + std::vector>>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R0, + const ModuleBase::Vector3& A, + const bool& calc_r, + const char* timer_name); + +/** + * @brief Compute with explicit quadrature settings. + * + * If calc_r is true, nlm[1..3] also store the Cartesian position moments. + */ +void snap_projector_half_tddft(const LCAO_Orbitals& orb, + const std::vector& projector_channels, + std::vector>>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R0, + const ModuleBase::Vector3& A, + const bool& calc_r, + const SnapIntegrationOptions& options, + const char* timer_name); + +} // namespace module_rt + +#endif diff --git a/source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp b/source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp index 0db4b7a9762..2b1cf62728f 100644 --- a/source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp +++ b/source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp @@ -1,63 +1,25 @@ #include "snap_psibeta_half_tddft.h" -#include "source_base/constants.h" -#include "source_base/math_integral.h" -#include "source_base/math_polyint.h" -#include "source_base/timer.h" -#include "source_base/ylm.h" - -#include -#include -#include - namespace module_rt { -/** - * @brief Initialize Gauss-Legendre grid points and weights. - * Thread-safe initialization using static local variable. - * - * @param grid_size Number of grid points (140) - * @param gl_x Output: Grid points in [-1, 1] - * @param gl_w Output: Weights - */ -static void init_gauss_legendre_grid(int grid_size, std::vector& gl_x, std::vector& gl_w) +void snap_psibeta_half_tddft(const LCAO_Orbitals& orb, + const InfoNonlocal& infoNL_, + std::vector>>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R0, + const int& T0, + const ModuleBase::Vector3& A, + const bool& calc_r) { - static bool init = false; -// Thread-safe initialization -#pragma omp critical(init_gauss_legendre) - { - if (!init) - { - ModuleBase::Integral::Gauss_Legendre_grid_and_weight(grid_size, gl_x.data(), gl_w.data()); - init = true; - } - } + SnapIntegrationOptions options; + snap_psibeta_half_tddft(orb, infoNL_, nlm, R1, T1, L1, m1, N1, R0, T0, A, calc_r, options); } -/** - * @brief Main function to calculate overlap integrals - * and its derivatives (if calc_r is true). - * - * This function integrates the overlap between a local orbital phi (at R1) - * and a non-local projector beta (at R0), modulated by a plane-wave-like phase factor - * exp^{-iAr}, where A is a vector potential. - * - * @param orb LCAO Orbitals information - * @param infoNL_ Non-local pseudopotential information - * @param nlm Output: - * nlm[0] : - * nlm[1, 2, 3] : , a = x, y, z (if calc_r=true) - * @param R1 Position of atom 1 (orbital phi) - * @param T1 Type of atom 1 - * @param L1 Angular momentum of orbital phi - * @param m1 Magnetic quantum number of orbital phi - * @param N1 Radial quantum number of orbital phi - * @param R0 Position of atom 0 (projector beta) - * @param T0 Type of atom 0 - * @param A Vector potential A (or related field vector) - * @param calc_r Whether to calculate position operator matrix elements - */ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb, const InfoNonlocal& infoNL_, std::vector>>& nlm, @@ -69,283 +31,27 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb, const ModuleBase::Vector3& R0, const int& T0, const ModuleBase::Vector3& A, - const bool& calc_r) + const bool& calc_r, + const SnapIntegrationOptions& options) { - ModuleBase::timer::start("module_rt", "snap_psibeta_half_tddft"); - - // 1. Initialization and Early Exits - const int nproj = infoNL_.nproj[T0]; - - // Resize output vector based on whether position operator matrix elements are needed - int required_size = calc_r ? 4 : 1; - if (nlm.size() != required_size) - nlm.resize(required_size); - - if (nproj == 0) - return; - - // 2. Determine total number of projectors and identify active ones based on cutoff - int natomwfc = 0; - std::vector calproj(nproj, false); - - const double Rcut1 = orb.Phi[T1].getRcut(); - const ModuleBase::Vector3 dRa = R0 - R1; - const double distance10 = dRa.norm(); - - bool any_active = false; - for (int ip = 0; ip < nproj; ip++) - { - const int L0 = infoNL_.Beta[T0].Proj[ip].getL(); - natomwfc += 2 * L0 + 1; - - const double Rcut0 = infoNL_.Beta[T0].Proj[ip].getRcut(); - if (distance10 <= (Rcut1 + Rcut0)) - { - calproj[ip] = true; - any_active = true; - } - } - - // Initialize output values to zero and resize inner vectors - for (auto& x: nlm) - { - x.assign(natomwfc, 0.0); - } - - if (!any_active) - { - ModuleBase::timer::end("module_rt", "snap_psibeta_half_tddft"); - return; - } - - // 3. Prepare Orbital Data (Phi) - const auto& phi_ln = orb.Phi[T1].PhiLN(L1, N1); - const int mesh_r1 = phi_ln.getNr(); - const double* psi_1 = phi_ln.getPsi(); - const double dk_1 = phi_ln.getDk(); - - // 4. Prepare Integration Grids - const int radial_grid_num = 140; - const int angular_grid_num = 110; - - // Cached standard Gauss-Legendre grid - static std::vector gl_x(radial_grid_num); - static std::vector gl_w(radial_grid_num); - init_gauss_legendre_grid(radial_grid_num, gl_x, gl_w); - - // Buffers for mapped radial grid - std::vector r_radial(radial_grid_num); - std::vector w_radial(radial_grid_num); - - // Precompute A dot r_angular (A * u_angle) for the Lebedev grid - std::vector A_dot_lebedev(angular_grid_num); - for (int ian = 0; ian < angular_grid_num; ++ian) - { - A_dot_lebedev[ian] = A.x * ModuleBase::Integral::Lebedev_Laikov_grid110_x[ian] - + A.y * ModuleBase::Integral::Lebedev_Laikov_grid110_y[ian] - + A.z * ModuleBase::Integral::Lebedev_Laikov_grid110_z[ian]; - } - - // Reuseable buffers for inner loops to avoid allocation - std::vector> result_angular; // Accumulator for angular integration - // Accumulators for position operator components - std::vector> res_ang_x, res_ang_y, res_ang_z; - - std::vector rly1((L1 + 1) * (L1 + 1)); // Spherical harmonics buffer for L1 - std::vector> rly0_cache(angular_grid_num); // Cache for L0 Ylm - - // 5. Loop over Projectors (Beta) - int index_offset = 0; - for (int nb = 0; nb < nproj; nb++) - { - const int L0 = infoNL_.Beta[T0].Proj[nb].getL(); - const int num_m0 = 2 * L0 + 1; - - if (!calproj[nb]) - { - index_offset += num_m0; - continue; - } - - const auto& proj = infoNL_.Beta[T0].Proj[nb]; - const int mesh_r0 = proj.getNr(); - const double* beta_r = proj.getBeta_r(); - const double* radial0 = proj.getRadial(); - const double dk_0 = proj.getDk(); - const double Rcut0 = proj.getRcut(); - - // 5.1 Map Gauss-Legendre grid to radial interval [r_min, r_max] - double r_min = radial0[0]; - double r_max = radial0[mesh_r0 - 1]; - double xl = (r_max - r_min) * 0.5; - double xmean = (r_max + r_min) * 0.5; - - for (int i = 0; i < radial_grid_num; ++i) - { - r_radial[i] = xmean + xl * gl_x[i]; - w_radial[i] = xl * gl_w[i]; - } - - const double A_phase = A * R0; - const std::complex exp_iAR0 = std::exp(ModuleBase::IMAG_UNIT * A_phase); - - // 5.2 Precompute Spherical Harmonics (Ylm) for L0 on angular grid - // Since L0 changes with projector, we compute this per projector loop. - for (int ian = 0; ian < angular_grid_num; ++ian) - { - ModuleBase::Ylm::rl_sph_harm(L0, - ModuleBase::Integral::Lebedev_Laikov_grid110_x[ian], - ModuleBase::Integral::Lebedev_Laikov_grid110_y[ian], - ModuleBase::Integral::Lebedev_Laikov_grid110_z[ian], - rly0_cache[ian]); - } - - // Resize accumulators if needed - if (result_angular.size() < num_m0) - { - result_angular.resize(num_m0); - if (calc_r) - { - res_ang_x.resize(num_m0); - res_ang_y.resize(num_m0); - res_ang_z.resize(num_m0); - } - } - - // 5.3 Radial Integration Loop - for (int ir = 0; ir < radial_grid_num; ir++) - { - const double r_val = r_radial[ir]; - - // Reset angular accumulators for this radial shell - std::fill(result_angular.begin(), result_angular.begin() + num_m0, 0.0); - if (calc_r) - { - std::fill(res_ang_x.begin(), res_ang_x.begin() + num_m0, 0.0); - std::fill(res_ang_y.begin(), res_ang_y.begin() + num_m0, 0.0); - std::fill(res_ang_z.begin(), res_ang_z.begin() + num_m0, 0.0); - } - - // 5.4 Angular Integration Loop (Lebedev Grid) - for (int ian = 0; ian < angular_grid_num; ian++) - { - const double x = ModuleBase::Integral::Lebedev_Laikov_grid110_x[ian]; - const double y = ModuleBase::Integral::Lebedev_Laikov_grid110_y[ian]; - const double z = ModuleBase::Integral::Lebedev_Laikov_grid110_z[ian]; - const double w_ang = ModuleBase::Integral::Lebedev_Laikov_grid110_w[ian]; - - // Vector r = r_val * u_angle - double rx = r_val * x; - double ry = r_val * y; - double rz = r_val * z; - - // Vector r' = r + R0 - R1 = r + dRa - double tx = rx + dRa.x; - double ty = ry + dRa.y; - double tz = rz + dRa.z; - - double tnorm = std::sqrt(tx * tx + ty * ty + tz * tz); - - // If r' is outside the cutoff of Phi(r'), skip - if (tnorm > Rcut1) - continue; - - // Compute Ylm for L1 at direction r' - if (tnorm > 1e-10) - { - double inv_tnorm = 1.0 / tnorm; - ModuleBase::Ylm::rl_sph_harm(L1, tx * inv_tnorm, ty * inv_tnorm, tz * inv_tnorm, rly1); - } - else - { - // At origin, only Y_00 is non-zero (if using real spherical harmonics convention) - ModuleBase::Ylm::rl_sph_harm(L1, 0.0, 0.0, 1.0, rly1); - } - - // Calculate common phase and weight factor - // phase = A * r = r_val * (A * u_angle) - const double phase = r_val * A_dot_lebedev[ian]; - const std::complex exp_iAr = std::exp(ModuleBase::IMAG_UNIT * phase); - - // Interpolate Psi at |r'| - double interp_psi = ModuleBase::PolyInt::Polynomial_Interpolation(psi_1, mesh_r1, dk_1, tnorm); - - const int offset_L1 = L1 * L1 + m1; - const double ylm_L1_val = rly1[offset_L1]; - - // Combined factor: exp(iAr) * Y_L1m1(r') * Psi(|r'|) * weight_angle - const std::complex common_factor = exp_iAr * ylm_L1_val * interp_psi * w_ang; - - // Retrieve precomputed Y_L0m0(r) - const std::vector& rly0_vec = rly0_cache[ian]; - const int offset_L0 = L0 * L0; - - // Accumulate results for all m0 components - for (int m0 = 0; m0 < num_m0; m0++) - { - std::complex term = common_factor * rly0_vec[offset_L0 + m0]; - result_angular[m0] += term; - - if (calc_r) - { - // Position operator r_op = r + R0 - // Note: Term involves (r_op)_a * exp(...). - double r_op_x = rx + R0.x; - double r_op_y = ry + R0.y; - double r_op_z = rz + R0.z; - - res_ang_x[m0] += term * r_op_x; - res_ang_y[m0] += term * r_op_y; - res_ang_z[m0] += term * r_op_z; - } - } - } // End Angular Loop - - // 5.5 Combine Radial and Angular parts - // Interpolate Beta(|r|) - // Note: The original code implies beta_r stores values that might need scaling or are just the function - // values. Typically radial integration is \int f(r) r^2 dr. Here we have factor: beta_val * r_radial[ir] * - // w_radial[ir] w_radial includes the Jacobian for the change of variable from [-1,1] to [r_min, r_max]. The - // extra r_radial[ir] suggests either beta is stored as r*beta, or we are doing \int ... r dr (2D?), or - // Jacobian r^2 is split. Assuming original logic is correct. - - double beta_val = ModuleBase::PolyInt::Polynomial_Interpolation(beta_r, mesh_r0, dk_0, r_radial[ir]); - - double radial_factor = beta_val * r_radial[ir] * w_radial[ir]; - - int current_idx = index_offset; - for (int m0 = 0; m0 < num_m0; m0++) - { - // Final accumulation into global nlm array - // Add phase exp(i A * R0) - nlm[0][current_idx] += radial_factor * result_angular[m0] * exp_iAR0; - - if (calc_r) - { - nlm[1][current_idx] += radial_factor * res_ang_x[m0] * exp_iAR0; - nlm[2][current_idx] += radial_factor * res_ang_y[m0] * exp_iAR0; - nlm[3][current_idx] += radial_factor * res_ang_z[m0] * exp_iAR0; - } - current_idx++; - } - - } // End Radial Loop - - index_offset += num_m0; - } // End Projector Loop + std::vector channels; + channels.reserve(infoNL_.nproj[T0]); - // 6. Final Conjugation - // Apply conjugation to all elements as per convention = * - for (int dim = 0; dim < nlm.size(); dim++) + // Convert nonlocal pseudopotential beta projectors to the shared grid integrator input. + for (int ip = 0; ip < infoNL_.nproj[T0]; ++ip) { - for (auto& x: nlm[dim]) - { - x = std::conj(x); - } + const auto& proj = infoNL_.Beta[T0].Proj[ip]; + ProjectorChannel channel; + channel.l = proj.getL(); + channel.mesh = proj.getNr(); + channel.dk = proj.getDk(); + channel.rcut = proj.getRcut(); + channel.radial_values = proj.getBeta_r(); + channel.radial_grid = proj.getRadial(); + channels.push_back(channel); } - assert(index_offset == natomwfc); - ModuleBase::timer::end("module_rt", "snap_psibeta_half_tddft"); + snap_projector_half_tddft(orb, channels, nlm, R1, T1, L1, m1, N1, R0, A, calc_r, options, "snap_psibeta_half_tddft"); } -} // namespace module_rt \ No newline at end of file +} // namespace module_rt diff --git a/source/source_lcao/module_rt/snap_psibeta_half_tddft.h b/source/source_lcao/module_rt/snap_psibeta_half_tddft.h index 78aab1f376a..2644fbe6baf 100644 --- a/source/source_lcao/module_rt/snap_psibeta_half_tddft.h +++ b/source/source_lcao/module_rt/snap_psibeta_half_tddft.h @@ -1,30 +1,52 @@ #ifndef SNAP_PSIBETA_HALF_TDDFT #define SNAP_PSIBETA_HALF_TDDFT -#include -#include - #include "source_base/vector3.h" #include "source_basis/module_ao/ORB_read.h" #include "source_cell/setup_nonlocal.h" +#include "source_lcao/module_rt/snap_projector_half_tddft.h" + +#include +#include namespace module_rt { - // calculate the tddft nonlocal potential term - void snap_psibeta_half_tddft( - const LCAO_Orbitals &orb, - const InfoNonlocal &infoNL_, - std::vector>> &nlm, - const ModuleBase::Vector3 &R1, - const int &T1, - const int &L1, - const int &m1, - const int &N1, - const ModuleBase::Vector3 &R0, // The projector. - const int &T0, - const ModuleBase::Vector3 &A, - const bool &calc_r - ); +/** + * @brief Compute RT-TDDFT velocity-gauge beta-projector overlaps. + * + * This overload uses the production quadrature settings. + */ +void snap_psibeta_half_tddft(const LCAO_Orbitals& orb, + const InfoNonlocal& infoNL_, + std::vector>>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R0, // The projector. + const int& T0, + const ModuleBase::Vector3& A, + const bool& calc_r); + +/** + * @brief Compute RT-TDDFT velocity-gauge beta-projector overlaps. + * + * This overload is used by tests to select the radial and Lebedev-Laikov grids. + */ +void snap_psibeta_half_tddft(const LCAO_Orbitals& orb, + const InfoNonlocal& infoNL_, + std::vector>>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R0, + const int& T0, + const ModuleBase::Vector3& A, + const bool& calc_r, + const SnapIntegrationOptions& options); } // namespace module_rt diff --git a/source/source_lcao/module_rt/test/CMakeLists.txt b/source/source_lcao/module_rt/test/CMakeLists.txt index cb24761d718..2dc950f3a82 100644 --- a/source/source_lcao/module_rt/test/CMakeLists.txt +++ b/source/source_lcao/module_rt/test/CMakeLists.txt @@ -34,3 +34,8 @@ AddTest( SOURCES propagator_test1.cpp propagator_test2.cpp propagator_test3.cpp ../propagator.cpp ../propagator_cn2.cpp ../propagator_taylor.cpp ../propagator_etrs.cpp ) +AddTest( + TARGET MODULE_LCAO_tddft_snap_psibeta_half_test + LIBS parameter ${math_libs} base device orb numerical_atomic_orbitals tddft_test_lib + SOURCES snap_psibeta_half_tddft_test.cpp ../snap_projector_half_tddft.cpp ../snap_psibeta_half_tddft.cpp +) diff --git a/source/source_lcao/module_rt/test/snap_psibeta_half_tddft_test.cpp b/source/source_lcao/module_rt/test/snap_psibeta_half_tddft_test.cpp new file mode 100644 index 00000000000..3fc7f2b5d02 --- /dev/null +++ b/source/source_lcao/module_rt/test/snap_psibeta_half_tddft_test.cpp @@ -0,0 +1,207 @@ +#include "source_lcao/module_rt/snap_psibeta_half_tddft.h" + +#include "source_base/ylm.h" +#include "source_basis/module_nao/radial_collection.h" +#include "source_basis/module_nao/two_center_integrator.h" +#include "source_cell/setup_nonlocal.h" + +#include +#include +#include +#include +#include +#include + +InfoNonlocal::InfoNonlocal() +{ + this->Beta = new Numerical_Nonlocal[1]; + this->nproj = nullptr; + this->nprojmax = 0; + this->rcutmax_Beta = 0.0; +} + +InfoNonlocal::~InfoNonlocal() +{ + delete[] this->Beta; + delete[] this->nproj; +} + +namespace +{ +struct ComparisonStats +{ + double max_real_diff = 0.0; + double max_imag_abs = 0.0; + double max_reference_abs = 0.0; +}; + +class SnapPsibetaHalfTddftTest : public ::testing::Test +{ + protected: + void SetUp() override + { + ModuleBase::Ylm::set_coefficients(); + + const std::string root = "../../../../../"; + const std::string orb_file = "tests/PP_ORB/C_gga_8au_100Ry_2s2p1d.orb"; + const std::string full_orb_file = root + orb_file; + const std::string orbital_files[1] = {orb_file}; + + std::ofstream ofs("snap_psibeta_half_tddft_test.log"); + orb.init(ofs, 1, root, orbital_files, "", 2, 100.0, 0.01, 0.01, 30.0, false, 0, false, false, 0); + + build_fake_beta_projectors(); + + orb_radials.build(1, &full_orb_file, 'o'); + beta_radials.build(1, info_nl.Beta); + + const double rmax = std::max(orb_radials.rcut_max(), beta_radials.rcut_max()); + const double cutoff = 2.0 * rmax; + const int nr = static_cast(rmax / 0.01) + 1; + + orb_radials.set_uniform_grid(true, nr, cutoff, 'i', true); + beta_radials.set_uniform_grid(true, nr, cutoff, 'i', true); + overlap_orb_beta.tabulate(orb_radials, beta_radials, 'S', nr, cutoff); + } + + void build_fake_beta_projectors() + { + const int nproj = 2; + std::vector beta_lm(nproj); + + for (int iproj = 0; iproj < nproj; ++iproj) + { + const int l = iproj; + const auto& phi_ln = orb.Phi[0].PhiLN(l, 0); + beta_lm[iproj].set_NL_proj("C", + 0, + l, + phi_ln.getNr(), + phi_ln.getRab(), + phi_ln.getRadial(), + phi_ln.getPsi_r(), + orb.get_kmesh(), + orb.get_dk(), + orb.get_dr_uniform()); + } + + info_nl.nproj = new int[1]; + info_nl.nproj[0] = nproj; + info_nl.nprojmax = nproj; + info_nl.Beta[0].set_type_info(0, "C", "NC", 1, nproj, beta_lm.data()); + info_nl.rcutmax_Beta = info_nl.Beta[0].get_rcut_max(); + } + + static int abacus_m_to_m(const int m) + { + return (m % 2 == 0) ? -m / 2 : (m + 1) / 2; + } + + ComparisonStats compare_zero_vector_potential(const int lebedev_grid_points) + { + const ModuleBase::Vector3 R0(0.1, -0.2, 0.3); + const ModuleBase::Vector3 R1(0.4, 0.2, -0.1); + const ModuleBase::Vector3 zero_A(0.0, 0.0, 0.0); + module_rt::SnapIntegrationOptions options; + options.lebedev_grid_points = lebedev_grid_points; + + ComparisonStats stats; + + for (int L1 = 0; L1 <= orb.Phi[0].getLmax(); ++L1) + { + for (int N1 = 0; N1 < orb.Phi[0].getNchi(L1); ++N1) + { + for (int m1 = 0; m1 < 2 * L1 + 1; ++m1) + { + std::vector>> grid_nlm; + module_rt::snap_psibeta_half_tddft(orb, + info_nl, + grid_nlm, + R1, + 0, + L1, + m1, + N1, + R0, + 0, + zero_A, + false, + options); + + std::vector> tci_nlm; + overlap_orb_beta.snap(0, L1, N1, abacus_m_to_m(m1), 0, R0 - R1, false, tci_nlm); + + EXPECT_FALSE(grid_nlm.empty()); + EXPECT_FALSE(tci_nlm.empty()); + if (grid_nlm.empty() || tci_nlm.empty()) + { + continue; + } + EXPECT_EQ(grid_nlm[0].size(), tci_nlm[0].size()); + if (grid_nlm[0].size() != tci_nlm[0].size()) + { + continue; + } + + for (size_t i = 0; i < grid_nlm[0].size(); ++i) + { + stats.max_real_diff + = std::max(stats.max_real_diff, std::abs(grid_nlm[0][i].real() - tci_nlm[0][i])); + stats.max_imag_abs = std::max(stats.max_imag_abs, std::abs(grid_nlm[0][i].imag())); + stats.max_reference_abs = std::max(stats.max_reference_abs, std::abs(tci_nlm[0][i])); + } + } + } + } + + return stats; + } + + LCAO_Orbitals orb; + InfoNonlocal info_nl; + RadialCollection orb_radials; + RadialCollection beta_radials; + TwoCenterIntegrator overlap_orb_beta; +}; +} // namespace + +TEST_F(SnapPsibetaHalfTddftTest, ZeroVectorPotentialMatchesTwoCenterIntegral) +{ + const double real_tolerance = 5.0e-8; + const double imag_tolerance = 1.0e-12; + const ComparisonStats stats = compare_zero_vector_potential(110); + + EXPECT_LT(stats.max_real_diff, real_tolerance) << "max reference abs = " << stats.max_reference_abs; + EXPECT_LT(stats.max_imag_abs, imag_tolerance) << "max reference abs = " << stats.max_reference_abs; +} + +TEST_F(SnapPsibetaHalfTddftTest, ZeroVectorPotentialHighOrderGridMatchesTwoCenterIntegral) +{ + const double real_tolerance = 5.0e-8; + const double imag_tolerance = 1.0e-12; + const ComparisonStats stats = compare_zero_vector_potential(590); + + EXPECT_LT(stats.max_real_diff, real_tolerance) << "max reference abs = " << stats.max_reference_abs; + EXPECT_LT(stats.max_imag_abs, imag_tolerance) << "max reference abs = " << stats.max_reference_abs; +} + +TEST_F(SnapPsibetaHalfTddftTest, ZeroVectorPotentialPositionMomentsAreReal) +{ + const ModuleBase::Vector3 R0(-0.3, 0.2, 0.1); + const ModuleBase::Vector3 R1(0.2, -0.1, 0.4); + const ModuleBase::Vector3 zero_A(0.0, 0.0, 0.0); + const double tolerance = 1.0e-12; + + std::vector>> nlm; + module_rt::snap_psibeta_half_tddft(orb, info_nl, nlm, R1, 0, 1, 1, 0, R0, 0, zero_A, true); + + ASSERT_EQ(nlm.size(), 4); + for (const auto& dim: nlm) + { + ASSERT_EQ(dim.size(), 4); + for (const std::complex& value: dim) + { + EXPECT_NEAR(value.imag(), 0.0, tolerance); + } + } +}