Perf: speed up the LR kernel-to-potential integrands by ~45%

`PotHxcLR::cal_v_eff` dominates an LR run: 74% of the total for 07_LiF/pbe
(949 s of 1289 s over 541 calls). Roughly half of that was FFT; the rest was
serial element-wise work and allocator traffic. Five changes, no algorithmic
change:

1. OpenMP on the element-wise loops. Everything these loops call
   (`grad_rho`, `grad_dot`, the Hartree kernel) was already parallel; the
   integrands here were the only part still running on one thread.

2. Raw pointers instead of `.at()` (43 sites). The bounds check is a branch
   per access, ~10 per grid point, and it blocks vectorization. The
   outer-vector lookups (`drho_gs.at(0)`) are hoisted out of the loops, where
   they keep their bounds check for free.

3. A scratch pool instead of per-call temporaries. `drho`, `gdot_terms` and
   `vxc_tmp` were allocated and value-initialized on every call and then
   overwritten before being read -- ~450 MB of pointless memset per call at a
   200^3 grid, plus the page faults. `grad_dot` assigns rather than
   accumulates, so even `vxc_tmp`'s zero-fill was dead. The pool is shared by
   all instances (a run holds three `PotHxcLR`) so it does not multiply.

4. `add_v_hartree` replaces `H_Hartree_pw::v_hartree`, which (a) re-did the
   forward FFT of rho^X that the GGA branch needs anyway -- 1 of the 10 FFTs
   per call was pure duplication, (b) reduced a Hartree "energy" of the
   transition density through `Parallel_Reduce::reduce_pool` every call, a
   collective nobody reads that also clobbers the global
   `H_Hartree_pw::hartree_energy`, and (c) returned a `matrix` by value, to
   which `v_eff += 2 * (...)` added a second full-size temporary.

5. The nspin=2 singlet/triplet combinations `v2rho2_uu -+ v2rho2_ud` and
   `2*vsigma_uu -+ vsigma_ud` are pre-contracted once per potential instead
   of being rebuilt at every grid point of every call, which also replaces
   two strided reads with one contiguous one. Costs 8-16 B/point, and only
   for closed-shell nspin=2.

`PotGradXCLR::cal_v_eff` (the g^xc branch feeding the Z-vector RHS) gets 1-3
of the same treatment, for -20% to -28%. The shared pool matters more there:
a `PotGradXCLR` is constructed inside the loop over excited states, so
per-object buffers would never be reused at all.

Measured (single node, 16 threads), analytic gradients unchanged:

  01_Si/lda    93 s -> 51 s   (-45%)   cal_v_eff  66 ->   45 ms/call
  01_Si/pbe   229 s -> 125 s  (-46%)   cal_v_eff 356 ->  268 ms/call
  02_Li2/lda  165 s -> 86 s   (-48%)   cal_v_eff 559 ->  278 ms/call
  02_Li2/pbe  530 s -> 253 s  (-52%)   cal_v_eff 2651 -> 1553 ms/call

(The total also benefits from solving the Z-vector equation once instead of
three times, a separate fix; the per-call figures above isolate this commit.)

Excitation energies are identical to every printed digit; the nspin=1 cases
agree to 1e-14, and of 60 force components per nspin=2 case, 2-3 differ by
exactly one unit in the last printed digit (the Hartree prefactor is now
associated differently).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QMV7xFZ69hhc3HnVBHbD96

Author: maki49

source/source_lcao/module_lr/Grad/xc/pot_grad_xc.cpp

@@ -8,6 +8,24 @@
 #include <set>
 namespace LR
 {
+    using Vec3 = ModuleBase::Vector3<double>;
+    PotGradXCLR::Scratch& PotGradXCLR::scratch()
+    {
+        static Scratch sc;   // see the comment on `Scratch` in the header
+        return sc;
+    }
+
+    void PotGradXCLR::Scratch::alloc(const int nrxx, const bool gga)
+    {
+        // resize() on an already-large vector is a no-op, so only the first call allocates.
+        if (static_cast<int>(this->vtmp.size()) < nrxx) { this->vtmp.resize(nrxx); }
+        if (gga)
+        {
+            if (static_cast<int>(this->gdot.size()) < nrxx) { this->gdot.resize(nrxx); }
+            if (static_cast<int>(this->drho1.size()) < nrxx) { this->drho1.resize(nrxx); }
+        }
+    }
+
     // constructor for exchange-correlation kernel
     PotGradXCLR::PotGradXCLR(const KernelXC& xc_kernel, const ModulePW::PW_Basis& rho_basis, const UnitCell& ucell,
         const int& nrxx, const bool triplet)
@@ -41,46 +59,66 @@ namespace LR
 
         if (func_type == 1) // LDA: only the $g^{\rho\rho\rho}$ term survives
         {
+            const double* const a_s2 = g.a_s2.data();
+            const double* const r1 = rho[0];
+            double* const v = v_eff.c;
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static)
+#endif
             for (int ir = 0;ir < nrxx_;++ir)
             {
-                v_eff(0, ir) += ModuleBase::e2 * g.a_s2.at(ir) * rho[0][ir] * rho[0][ir];
+                v[ir] += ModuleBase::e2 * a_s2[ir] * r1[ir] * r1[ir];
             }
         }
         else if (func_type == 2 || func_type == 4)  // GGA or HYB_GGA
         {
-            std::vector<ModuleBase::Vector3<double>> drho1(nrxx_);  // transition density gradient
-            LR_Util::grad(rho[0], drho1.data(), this->rho_basis_, this->tpiba_);
+            scratch().alloc(nrxx_, /*gga=*/true);
+            Vec3* const drho1 = scratch().drho1.data();   // transition density gradient
+            LR_Util::grad(rho[0], drho1, this->rho_basis_, this->tpiba_);
 
-            std::vector<double> v_tmp(nrxx_, 0.0);
+            double* const v_tmp = scratch().vtmp.data();
+            Vec3* const gdot_terms = scratch().gdot.data();
+            const Vec3* const dgs = kxc.drho_gs.at(0).data();
+            const double* const r1 = rho[0];
+            const double* const e_s2 = g.e_s2.data(); const double* const e_st = g.e_st.data();
+            const double* const e_t2 = g.e_t2.data(); const double* const e_q = g.e_q.data();
+            const double* const c_s = g.c_s.data();   const double* const c_t = g.c_t.data();
+            const double* const a_s2 = g.a_s2.data(); const double* const a_st = g.a_st.data();
+            const double* const a_t2 = g.a_t2.data(); const double* const a_q = g.a_q.data();
 
             // 1. the vector under the divergence, accumulated negated so that `grad_dot` yields
             //    $-\nabla\cdot\boldsymbol{E}$. The four $e$ coefficients share the same
             //    $\nabla\rho^{gs}$ direction (see `KernelXC::GxcCoef`), so it is pulled out of
             //    their sum -- exact, and it keeps this bandwidth-bound loop reading 4 doubles per
             //    point instead of 12.
-            std::vector<ModuleBase::Vector3<double>> gdot_terms(nrxx_);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static)
+#endif
             for (int ir = 0;ir < nrxx_;++ir)
             {
-                const ModuleBase::Vector3<double>& drho = kxc.drho_gs.at(0).at(ir);  // $\nabla\rho$
-                const double s = rho[0][ir];                                // $\rho^1$
-                const double t = drho * drho1.at(ir);                       // $\nabla\rho\cdot\nabla\rho^1$
-                const double q = drho1.at(ir) * drho1.at(ir);               // $\nabla\rho^1\cdot\nabla\rho^1$
-                const double e = g.e_s2.at(ir) * (s * s) + g.e_st.at(ir) * (s * t)
-                    + g.e_t2.at(ir) * (t * t) + g.e_q.at(ir) * q;
-                gdot_terms[ir] = -(drho * e + drho1.at(ir) * (g.c_s.at(ir) * s + g.c_t.at(ir) * t));
+                const Vec3& drho = dgs[ir];                     // $\nabla\rho$
+                const double s = r1[ir];                        // $\rho^1$
+                const double t = drho * drho1[ir];              // $\nabla\rho\cdot\nabla\rho^1$
+                const double q = drho1[ir] * drho1[ir];         // $\nabla\rho^1\cdot\nabla\rho^1$
+                const double e = e_s2[ir] * (s * s) + e_st[ir] * (s * t)
+                    + e_t2[ir] * (t * t) + e_q[ir] * q;
+                gdot_terms[ir] = -(drho * e + drho1[ir] * (c_s[ir] * s + c_t[ir] * t));
             }
-            XC_Functional::grad_dot(gdot_terms.data(), v_tmp.data(), &this->rho_basis_, this->tpiba_);
+            XC_Functional::grad_dot(gdot_terms, v_tmp, &this->rho_basis_, this->tpiba_);
 
             // 2. the local terms $A$
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static)
+#endif
             for (int ir = 0;ir < nrxx_;++ir)
             {
-                const double s = rho[0][ir];
-                const double t = kxc.drho_gs.at(0).at(ir) * drho1.at(ir);
-                const double q = drho1.at(ir) * drho1.at(ir);
-                v_tmp[ir] += g.a_s2.at(ir) * (s * s) + g.a_st.at(ir) * (s * t)
-                    + g.a_t2.at(ir) * (t * t) + g.a_q.at(ir) * q;
+                const double s = r1[ir];
+                const double t = dgs[ir] * drho1[ir];
+                const double q = drho1[ir] * drho1[ir];
+                v_tmp[ir] += a_s2[ir] * (s * s) + a_st[ir] * (s * t)
+                    + a_t2[ir] * (t * t) + a_q[ir] * q;
             }
-            BlasConnector::axpy(nrxx_, ModuleBase::e2, v_tmp.data(), 1, v_eff.c, 1);
+            BlasConnector::axpy(nrxx_, ModuleBase::e2, v_tmp, 1, v_eff.c, 1);
         }
         else
         {
@@ -90,5 +128,4 @@ namespace LR
 
         ModuleBase::timer::end("PotGradXCLR", "cal_v_eff");
     }
-
 }

source/source_lcao/module_lr/Grad/xc/pot_grad_xc.h

@@ -19,6 +19,19 @@ namespace LR
         /// kernel components from PotHxcLR
         const KernelXC& xc_kernel_components_;
         const bool triplet_ = false;
+
+    private:
+        /// Scratch, shared by every `PotGradXCLR` and grown on demand. These used to be
+        /// allocated (and value-initialized) on every call. Safe to share because
+        /// `cal_v_eff` is only ever entered from a single thread (all the OpenMP is inside).
+        struct Scratch
+        {
+            std::vector<ModuleBase::Vector3<double>> drho1;  ///< $\nabla\rho^1$
+            std::vector<ModuleBase::Vector3<double>> gdot;   ///< integrand of the divergence
+            std::vector<double> vtmp;                        ///< local part $A$
+            void alloc(const int nrxx, const bool gga);
+        };
+        static Scratch& scratch();
     };
 
 }