Skip to content

Refactor: simplify snap_psibeta_half_tddft.cpp #6846

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
mohanchen merged 3 commits into from
Jan 13, 2026
Merged

Refactor: simplify snap_psibeta_half_tddft.cpp #6846

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ab32f76
simplify snap_psibeta_half_tddft.cpp
dzzz2001 Jan 7, 2026
6254847
fix: use dk_1 instead of inv_dk_1 in Polynomial_Interpolation call
dzzz2001 Jan 12, 2026
50c79e8
Merge branch 'develop' into simplify-code
dzzz2001 Jan 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
148 changes: 60 additions & 88 deletions source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,26 +5,27 @@
#include "source_base/math_polyint.h"
#include "source_base/timer.h"
#include "source_base/ylm.h"
#include <vector>
#include <complex>

#include <cmath>
#include <complex>
#include <vector>

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<double>& gl_x, std::vector<double>& gl_w)
{
static bool init = false;
// Thread-safe initialization
#pragma omp critical(init_gauss_legendre)
// Thread-safe initialization
#pragma omp critical(init_gauss_legendre)
{
if (!init)
{
Expand All @@ -34,48 +35,17 @@ static void init_gauss_legendre_grid(int grid_size, std::vector<double>& gl_x, s
}
}

/**
* @brief Interpolate radial function value at a given distance using cubic interpolation.
*
* @param psi Pointer to radial function array
* @param mesh Number of mesh points
* @param inv_dk Inverse of grid spacing (1/dk)
* @param distance Distance to interpolate at
* @return double Interpolated value
*/
inline double interpolate_radial(
const double* psi,
int mesh,
double inv_dk,
double distance)
{
double position = distance * inv_dk;
int iq = static_cast<int>(position);

// Boundary check safe-guard
if (iq > mesh - 4) return 0.0;

const double x0 = position - static_cast<double>(iq);
const double x1 = 1.0 - x0;
const double x2 = 2.0 - x0;
const double x3 = 3.0 - x0;

// Cubic interpolation formula
return x1*x2*(psi[iq]*x3+psi[iq+3]*x0)/6.0
+ x0*x3*(psi[iq+1]*x2-psi[iq+2]*x1)/2.0;
}

/**
* @brief Main function to calculate overlap integrals <phi|exp^{-iAr}|beta>
* 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:
* @param nlm Output:
* nlm[0] : <phi|exp^{-iAr}|beta>
* nlm[1, 2, 3] : <phi|r_a * exp^{-iAr}|beta>, a = x, y, z (if calc_r=true)
* @param R1 Position of atom 1 (orbital phi)
Expand All @@ -96,7 +66,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R0,
const ModuleBase::Vector3<double>& R0,
const int& T0,
const ModuleBase::Vector3<double>& A,
const bool& calc_r)
Expand All @@ -105,27 +75,29 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,

// 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;
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<bool> calproj(nproj, false);

const double Rcut1 = orb.Phi[T1].getRcut();
const ModuleBase::Vector3<double> 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))
{
Expand All @@ -135,7 +107,8 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
}

// Initialize output values to zero and resize inner vectors
for (auto& x : nlm) {
for (auto& x: nlm)
{
x.assign(natomwfc, 0.0);
}

Expand All @@ -150,12 +123,11 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
const int mesh_r1 = phi_ln.getNr();
const double* psi_1 = phi_ln.getPsi();
const double dk_1 = phi_ln.getDk();
const double inv_dk_1 = 1.0 / dk_1;

// 4. Prepare Integration Grids
const int radial_grid_num = 140;
const int angular_grid_num = 110;

// Cached standard Gauss-Legendre grid
static std::vector<double> gl_x(radial_grid_num);
static std::vector<double> gl_w(radial_grid_num);
Expand All @@ -169,17 +141,17 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
std::vector<double> 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];
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<std::complex<double>> result_angular; // Accumulator for angular integration
// Accumulators for position operator components
std::vector<std::complex<double>> res_ang_x, res_ang_y, res_ang_z;
std::vector<double> rly1((L1 + 1) * (L1 + 1)); // Spherical harmonics buffer for L1

std::vector<double> rly1((L1 + 1) * (L1 + 1)); // Spherical harmonics buffer for L1
std::vector<std::vector<double>> rly0_cache(angular_grid_num); // Cache for L0 Ylm

// 5. Loop over Projectors (Beta)
Expand Down Expand Up @@ -207,8 +179,8 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
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)

for (int i = 0; i < radial_grid_num; ++i)
{
r_radial[i] = xmean + xl * gl_x[i];
w_radial[i] = xl * gl_w[i];
Expand All @@ -219,12 +191,13 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,

// 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]);
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
Expand Down Expand Up @@ -260,21 +233,22 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
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);

double tnorm = std::sqrt(tx * tx + ty * ty + tz * tz);

// If r' is outside the cutoff of Phi(r'), skip
if (tnorm > Rcut1) continue;
if (tnorm > Rcut1)
continue;

// Compute Ylm for L1 at direction r'
if (tnorm > 1e-10)
Expand All @@ -285,7 +259,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
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);
ModuleBase::Ylm::rl_sph_harm(L1, 0.0, 0.0, 1.0, rly1);
}

// Calculate common phase and weight factor
Expand All @@ -294,11 +268,11 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
const std::complex<double> exp_iAr = std::exp(ModuleBase::IMAG_UNIT * phase);

// Interpolate Psi at |r'|
double interp_psi = interpolate_radial(psi_1, mesh_r1, inv_dk_1, tnorm);
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<double> common_factor = exp_iAr * ylm_L1_val * interp_psi * w_ang;

Expand All @@ -319,7 +293,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
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;
Expand All @@ -329,16 +303,14 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,

// 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]);

// 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;
Expand All @@ -347,7 +319,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
// 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;
Expand All @@ -364,11 +336,11 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,

// 6. Final Conjugation
// Apply conjugation to all elements as per convention <phi|beta> = <beta|phi>*
for(int dim = 0; dim < nlm.size(); dim++)
for (int dim = 0; dim < nlm.size(); dim++)
{
for (auto &x : nlm[dim])
for (auto& x: nlm[dim])
{
x = std::conj(x);
x = std::conj(x);
}
}

Expand Down
Loading