Test: Add <phi|r_a|p_{l,m}> snap_psibeta checks

Author: AsTonyshment

source/Makefile.Objects

@@ -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\

source/source_io/module_hs/cal_r_overlap_R.cpp

@@ -187,11 +187,19 @@ void cal_r_overlap_R::construct_orbs_and_orb_r(const UnitCell& ucell,
         }
     }
 
-    iw2it.resize(PARAM.globalv.nlocal);
-    iw2ia.resize(PARAM.globalv.nlocal);
-    iw2iL.resize(PARAM.globalv.nlocal);
-    iw2iN.resize(PARAM.globalv.nlocal);
-    iw2im.resize(PARAM.globalv.nlocal);
+    int map_size = PARAM.globalv.nlocal;
+    int required_orbitals = 0;
+    for (int it = 0; it < ucell.ntype; ++it)
+    {
+        required_orbitals += ucell.atoms[it].nw * ucell.atoms[it].na;
+    }
+    map_size = std::max(map_size, required_orbitals);
+
+    iw2it.resize(map_size);
+    iw2ia.resize(map_size);
+    iw2iL.resize(map_size);
+    iw2iN.resize(map_size);
+    iw2im.resize(map_size);
 
     int iw = 0;
     for (int it = 0; it < ucell.ntype; it++)
@@ -419,11 +427,19 @@ void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucel
         }
     }
 
-    iw2it.resize(PARAM.globalv.nlocal);
-    iw2ia.resize(PARAM.globalv.nlocal);
-    iw2iL.resize(PARAM.globalv.nlocal);
-    iw2iN.resize(PARAM.globalv.nlocal);
-    iw2im.resize(PARAM.globalv.nlocal);
+    int map_size = PARAM.globalv.nlocal;
+    int required_orbitals = 0;
+    for (int it = 0; it < ucell.ntype; ++it)
+    {
+        required_orbitals += ucell.atoms[it].nw * ucell.atoms[it].na;
+    }
+    map_size = std::max(map_size, required_orbitals);
+
+    iw2it.resize(map_size);
+    iw2ia.resize(map_size);
+    iw2iL.resize(map_size);
+    iw2iN.resize(map_size);
+    iw2im.resize(map_size);
 
     int iw = 0;
     for (int it = 0; it < ucell.ntype; it++)

source/source_lcao/module_rt/snap_projector_half_tddft.cpp

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

source/source_lcao/module_rt/snap_projector_half_tddft.h

@@ -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,
@@ -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,

source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp

@@ -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];
@@ -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);
     }

source/source_lcao/module_rt/snap_psibeta_half_tddft.h

@@ -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,

source/source_lcao/module_rt/test/CMakeLists.txt

@@ -36,4 +36,21 @@ AddTest(
   TARGET MODULE_LCAO_tddft_snap_psibeta_half_test
   LIBS parameter 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
 )