Skip to content

Test: Add <phi|r_a|p_{l,m}> snap_psibeta unit test #7587

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

Open
AsTonyshment wants to merge 2 commits into
base: develop
Choose a base branch
from
+399 −374
Open

Test: Add <phi|r_a|p_{l,m}> snap_psibeta unit test #7587

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e3d6225
Test: Add <phi|r_a|p_{l,m}> snap_psibeta checks
AsTonyshment Jul 3, 2026
358403e
Style: Format cal_r_overlap_R
AsTonyshment Jul 3, 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
1 change: 1 addition & 0 deletions source/Makefile.Objects
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -635,6 +635,7 @@ OBJS_LCAO=evolve_elec.o\
td_folding.o\
td_info.o\
velocity_op.o\
snap_projector_half_tddft.o\
snap_psibeta_half_tddft.o\
solve_propagation.o\
boundary_fix.o\
Expand Down
356 changes: 166 additions & 190 deletions source/source_io/module_hs/cal_r_overlap_R.cpp
View file
Open in desktop

Large diffs are not rendered by default.

96 changes: 40 additions & 56 deletions source/source_io/module_hs/cal_r_overlap_R.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#ifndef CAL_R_OVERLAP_R_H
#define CAL_R_OVERLAP_R_H

#include "source_lcao/module_ri/abfs-vector3_order.h"
#include "source_base/sph_bessel_recursive.h"
#include "source_base/vector3.h"
#include "source_base/ylm.h"
Expand All @@ -14,6 +13,7 @@
#include "source_lcao/center2_orb-orb11.h"
#include "source_lcao/center2_orb-orb21.h"
#include "source_lcao/center2_orb.h"
#include "source_lcao/module_ri/abfs-vector3_order.h"

#include <map>
#include <set>
Expand All @@ -31,45 +31,39 @@ class cal_r_overlap_R
double sparse_threshold = 1e-10;
bool binary = false;

void init(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb);
void init_nonlocal(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb);
ModuleBase::Vector3<double> get_psi_r_psi(
const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R2,
const int& T2,
const int& L2,
const int& m2,
const int& N2
);
ModuleBase::Vector3<double> get_psi_r_gradpsi(
const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R2,
const int& T2,
const int& L2,
const int& m2,
const int& N2,
const ModuleBase::Vector3<double>& Efield,
const ModuleBase::Vector3<double>& dR
);
void get_psi_r_beta(
const UnitCell& ucell,
std::vector<std::vector<double>>& nlm,
const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R2,
const int& T2
);
void init(const UnitCell& ucell, const Parallel_Orbitals& pv, const LCAO_Orbitals& orb);
void init_nonlocal(const UnitCell& ucell, const Parallel_Orbitals& pv, const LCAO_Orbitals& orb);
ModuleBase::Vector3<double> get_psi_r_psi(const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R2,
const int& T2,
const int& L2,
const int& m2,
const int& N2);
ModuleBase::Vector3<double> get_psi_r_gradpsi(const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R2,
const int& T2,
const int& L2,
const int& m2,
const int& N2,
const ModuleBase::Vector3<double>& Efield,
const ModuleBase::Vector3<double>& dR);
void get_psi_r_beta(const UnitCell& ucell,
std::vector<std::vector<double>>& nlm,
const ModuleBase::Vector3<double>& R1,
const int& T1,
const int& L1,
const int& m1,
const int& N1,
const ModuleBase::Vector3<double>& R2,
const int& T2);
void out_rR(const UnitCell& ucell, const Grid_Driver& gd, const int& istep, const int precision = 16);
void out_rR_other(const UnitCell& ucell,
const int& istep,
Expand All @@ -78,8 +72,8 @@ class cal_r_overlap_R

private:
void initialize_orb_table(const UnitCell& ucell, const LCAO_Orbitals& orb);
void construct_orbs_and_orb_r(const UnitCell& ucell,const LCAO_Orbitals& orb);
void construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucell,const LCAO_Orbitals& orb);
void construct_orbs_and_orb_r(const UnitCell& ucell, const LCAO_Orbitals& orb);
void construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucell, const LCAO_Orbitals& orb);

std::vector<int> iw2ia;
std::vector<int> iw2iL;
Expand All @@ -94,25 +88,15 @@ class cal_r_overlap_R
std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> orbs;
std::vector<std::vector<Numerical_Orbital_Lm>> orbs_nonlocal;

std::map<
size_t,
std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, Center2_Orb::Orb11>>>>>>
std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, Center2_Orb::Orb11>>>>>>
center2_orb11;

std::map<
size_t,
std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, Center2_Orb::Orb21>>>>>>
std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, Center2_Orb::Orb21>>>>>>
center2_orb21_r;

std::map<
size_t,
std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, Center2_Orb::Orb11>>>>>
center2_orb11_nonlocal;
std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, Center2_Orb::Orb11>>>>> center2_orb11_nonlocal;

std::map<
size_t,
std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, Center2_Orb::Orb21>>>>>
center2_orb21_r_nonlocal;
std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, std::map<size_t, Center2_Orb::Orb21>>>>> center2_orb21_r_nonlocal;

const Parallel_Orbitals* ParaV = nullptr;
};
Expand Down
9 changes: 3 additions & 6 deletions source/source_lcao/module_rt/snap_projector_half_tddft.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,8 +104,7 @@ 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));
ModuleBase::WARNING_QUIT("snap_projector_half_tddft", "Unsupported Lebedev-Laikov grid size: " + std::to_string(ngrid));
}

if (ngrid == default_lebedev_grid_points)
Expand Down Expand Up @@ -142,9 +141,7 @@ AngularGridView angular_grid(const int ngrid)

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

const double projector_val = ModuleBase::PolyInt::Polynomial_Interpolation(channel.radial_times_r, channel.mesh, channel.dk, r);
return projector_val * r * w_radial;
}
} // namespace
Expand Down Expand Up @@ -264,7 +261,7 @@ void snap_projector_half_tddft(const LCAO_Orbitals& orb,
}

assert(channel.mesh > 0);
assert(channel.radial_values != nullptr);
assert(channel.radial_times_r != nullptr);
assert(channel.radial_grid != nullptr);

const double r_min = channel.radial_grid[0];
Expand Down
41 changes: 27 additions & 14 deletions source/source_lcao/module_rt/snap_projector_half_tddft.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,32 +10,44 @@
namespace module_rt
{

/**
* @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 Radial projector channel integrated against one LCAO orbital.
*
* radial_times_r stores r * p_l(r), where p_l(r) is the radial projector.
* The radial part of the volume integral is therefore evaluated as
* (r * p_l(r)) * r dr = p_l(r) * r^2 dr.
*/
struct ProjectorChannel
{
int l = 0;
int mesh = 0;
double dk = 0.0;
double rcut = 0.0;
const double* radial_values = nullptr;
const double* radial_times_r = nullptr;
const double* radial_grid = nullptr;
};

/**
* @brief Numerical quadrature settings for projector snapshots.
* @brief Compute projector overlaps with default quadrature settings.
*
* The default values reproduce the production RT-TDDFT path.
*/
struct SnapIntegrationOptions
{
int radial_grid_num = 140;
int lebedev_grid_points = 110;
};

/**
* @brief Compute <phi|exp(-i A r)|projector> with default quadrature settings.
* The shared integral is
* I_m(A) = <phi_{T1,L1,m1,N1}(r-R1)|exp(-i A.r)|p_{l,m}(r-R0)>.
*
* The phase A is given in Cartesian coordinates. The returned nlm[0] stores
* I_m(A) for all projector magnetic components. If calc_r is true, nlm[1..3]
* store
* R_a,m(A) = <phi|r_a exp(-i A.r)|p_{l,m}>.
*/
void snap_projector_half_tddft(const LCAO_Orbitals& orb,
const std::vector<ProjectorChannel>& projector_channels,
Expand All @@ -51,9 +63,10 @@ void snap_projector_half_tddft(const LCAO_Orbitals& orb,
const char* timer_name);

/**
* @brief Compute <phi|exp(-i A r)|projector> with explicit quadrature settings.
* @brief Compute projector overlaps with explicit quadrature settings.
*
* If calc_r is true, nlm[1..3] also store the Cartesian position moments.
* The ProjectorChannel radial convention is always r * p_l(r). Callers that
* own different physical projectors are responsible for passing that form.
*/
void snap_projector_half_tddft(const LCAO_Orbitals& orb,
const std::vector<ProjectorChannel>& projector_channels,
Expand Down
4 changes: 2 additions & 2 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 @@ -37,7 +37,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
std::vector<ProjectorChannel> channels;
channels.reserve(infoNL_.nproj[T0]);

// Convert nonlocal pseudopotential beta projectors to the shared grid integrator input.
// UPF nonlocal beta projectors already follow the r * beta_l(r) convention.
for (int ip = 0; ip < infoNL_.nproj[T0]; ++ip)
{
const auto& proj = infoNL_.Beta[T0].Proj[ip];
Expand All @@ -46,7 +46,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
channel.mesh = proj.getNr();
channel.dk = proj.getDk();
channel.rcut = proj.getRcut();
channel.radial_values = proj.getBeta_r();
channel.radial_times_r = proj.getBeta_r();
channel.radial_grid = proj.getRadial();
channels.push_back(channel);
}
Expand Down
1 change: 1 addition & 0 deletions source/source_lcao/module_rt/snap_psibeta_half_tddft.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@ namespace module_rt
/**
* @brief Compute RT-TDDFT velocity-gauge beta-projector overlaps.
*
* UPF beta projectors are stored as r * beta_l(r) and selected by atom type T0.
* This overload uses the production quadrature settings.
*/
void snap_psibeta_half_tddft(const LCAO_Orbitals& orb,
Expand Down
17 changes: 17 additions & 0 deletions source/source_lcao/module_rt/test/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,4 +36,21 @@ 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
../../center2_orb.cpp
../../center2_orb-orb11.cpp
../../center2_orb-orb21.cpp
../../../source_cell/setup_nonlocal.cpp
../../../source_cell/atom_spec.cpp
../../../source_cell/atom_pseudo.cpp
../../../source_cell/pseudo.cpp
../../../source_cell/read_pp.cpp
../../../source_cell/read_pp_complete.cpp
../../../source_cell/read_pp_upf201.cpp
../../../source_cell/read_pp_upf100.cpp
../../../source_cell/read_pp_vwr.cpp
../../../source_cell/read_pp_blps.cpp
../../../source_io/module_hs/cal_r_overlap_R.cpp
../../../source_io/module_hs/single_R_io.cpp
../../../source_io/module_hs/rr_sparse_writer.cpp
../../../source_pw/module_pwdft/soc.cpp
)
Loading
Loading