fix bugs

Author: abacus_fixer

source/source_lcao/force_stress_lcao.cpp

@@ -426,7 +426,32 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
         }
         if (PARAM.inp.dft_plus_u == 2)
         {
-            // Old DFT+U implementation (dft_plus_u==2) still needs ForceStressArrays
+            // The legacy dft_plus_u==2 force/stress path is currently broken.
+            //
+            // Background: Plus_U::force_stress relies on ForceStressArrays
+            // members DSloc_x/y/z (gamma_only) or DSloc_Rx/Ry/Rz (multik)
+            // and DH_r being pre-allocated and filled with dS/dR data by the
+            // main force flow (formerly ForceLcaoGamma::ftable). The DFT+U
+            // step 2 refactor (commit 70c54c9d5a, 2026-01-23) removed the
+            // main-flow ForceStressArrays because the operator-based force
+            // calculation no longer needs it, but the legacy dft_plus_u==2
+            // path still depends on it. The local fsr_dftu below is declared
+            // without allocating those arrays, so any call into
+            // cal_force_gamma / cal_stress_gamma / folding_matrix_k would
+            // pass nullptr to pdgemm_ and crash with SIGSEGV.
+            //
+            // Until the legacy path is restored or re-implemented, we
+            // explicitly reject dft_plus_u==2 with cal_force or cal_stress
+            // enabled. SCF-only runs (no force/stress) are unaffected
+            // because the energy is computed in cal_energy_correction,
+            // which does not touch DSloc arrays. Use dft_plus_u=1 for
+            // force/stress calculations.
+            if (isforce || isstress)
+            {
+                ModuleBase::WARNING_QUIT("Force_Stress_LCAO::getForceStress",
+                    "dft_plus_u==2 with cal_force or cal_stress is currently broken; "
+                    "please use dft_plus_u=1 instead. See notes in source/source_lcao/force_stress_lcao.cpp.");
+            }
             ForceStressArrays fsr_dftu;
             std::vector<std::vector<double>>* dmk_d = nullptr;
             std::vector<std::vector<std::complex<double>>>* dmk_c = nullptr;

source/source_lcao/module_dftu/dftu_force.cpp

@@ -36,6 +36,52 @@ void Plus_U::force_stress(const UnitCell& ucell,
     ModuleBase::TITLE("Plus_U", "force_stress");
     ModuleBase::timer::start("Plus_U", "force_stress");
 
+    // Defensive null check: the legacy dft_plus_u==2 force/stress path
+    // requires fsr.DSloc_x/y/z (gamma_only) or fsr.DSloc_Rx/Ry/Rz (multik)
+    // and fsr.DH_r to be allocated and filled by the caller. If the caller
+    // forgot to allocate them (as in force_stress_lcao.cpp where the local
+    // fsr_dftu is created without allocation), we fail early with a clear
+    // message instead of letting pdgemm_ dereference nullptr and crash.
+    // See force_stress_lcao.cpp for the historical background.
+    if (this->gamma_only_local)
+    {
+        if (this->cal_force
+            && (fsr.DSloc_x == nullptr || fsr.DSloc_y == nullptr || fsr.DSloc_z == nullptr))
+        {
+            ModuleBase::WARNING_QUIT("Plus_U::force_stress",
+                "fsr.DSloc_x/y/z are nullptr in gamma_only path; the caller must allocate and fill them. "
+                "See notes in source/source_lcao/force_stress_lcao.cpp.");
+        }
+        if (this->cal_stress
+            && (fsr.DSloc_x == nullptr || fsr.DSloc_y == nullptr || fsr.DSloc_z == nullptr
+                || fsr.DH_r == nullptr))
+        {
+            ModuleBase::WARNING_QUIT("Plus_U::force_stress",
+                "fsr.DSloc_x/y/z or fsr.DH_r is nullptr in gamma_only path; "
+                "the caller must allocate and fill them. "
+                "See notes in source/source_lcao/force_stress_lcao.cpp.");
+        }
+    }
+    else
+    {
+        if (this->cal_force
+            && (fsr.DSloc_Rx == nullptr || fsr.DSloc_Ry == nullptr || fsr.DSloc_Rz == nullptr))
+        {
+            ModuleBase::WARNING_QUIT("Plus_U::force_stress",
+                "fsr.DSloc_Rx/Ry/Rz are nullptr in multik path; the caller must allocate and fill them. "
+                "See notes in source/source_lcao/force_stress_lcao.cpp.");
+        }
+        if (this->cal_stress
+            && (fsr.DSloc_Rx == nullptr || fsr.DSloc_Ry == nullptr || fsr.DSloc_Rz == nullptr
+                || fsr.DH_r == nullptr))
+        {
+            ModuleBase::WARNING_QUIT("Plus_U::force_stress",
+                "fsr.DSloc_Rx/Ry/Rz or fsr.DH_r is nullptr in multik path; "
+                "the caller must allocate and fill them. "
+                "See notes in source/source_lcao/force_stress_lcao.cpp.");
+        }
+    }
+
     const int nlocal = this->nlocal;
 
     if (this->cal_force)