DeltaP LCAO force/stress implementation

Add cal_force_stress() to DeltaPOperator following DeltaSpin pattern:

  deltap_force_stress.hpp (new file, 270 lines):
    - cal_force_stress(): iterates over constrained atoms,
      calls intor_->snap(cal_deri=1) for projector derivatives,
      accumulates force/stress via cal_force_IJR/cal_stress_IJR
    - cal_force_IJR(): force = lambda * nlm_deriv * nlm * DM
      (same structure as DeltaSpin, scalar lambda instead of Vector3)
    - cal_stress_IJR(): stress = force * r_vector projection

  deltap_lcao.h: add force/stress method declarations +
    static s_stored_lambda for cross-module lambda access

  deltap_lcao.cpp: include deltap_force_stress.hpp +
    static member definition

  FORCE_STRESS.cpp: create temp DeltaPOperator with stored lambda,
    call cal_force_stress(), accumulate into fcs/scs

  esolver_ks_lcao.cpp: store lambda before getForceStress

BN LCAO + PW regression pass.

Author: dyzheng

source/source_esolver/esolver_ks_lcao.cpp

@@ -236,6 +236,18 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
 
     Force_Stress_LCAO<TK> fsl(this->RA, ucell.nat);
 
+    // Store DeltaP lambda for force/stress computation
+    if (PARAM.inp.deltap_switch && PARAM.inp.deltap_corr)
+    {
+        auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
+        if (hamilt_lcao != nullptr)
+        {
+            auto* dp_op = hamilt_lcao->get_dp_operator();
+            if (dp_op != nullptr)
+                hamilt::DeltaPOperator<TK, TR>::store_lambda_for_force(dp_op->get_lambda());
+        }
+    }
+
     deepks.dpks_out_type = "tot";  // for deepks method
 
     fsl.getForceStress(ucell, PARAM.inp.cal_force, PARAM.inp.cal_stress, 

source/source_lcao/FORCE_STRESS.cpp

@@ -21,6 +21,7 @@
 #endif
 #include "source_lcao/module_operator_lcao/dftu_lcao.h"
 #include "source_lcao/module_operator_lcao/dspin_lcao.h"
+#include "source_lcao/module_operator_lcao/deltap_lcao.h"
 #include "source_lcao/module_operator_lcao/nonlocal.h"
 #include "source_lcao/module_operator_lcao/ekinetic.h"
 #include "source_lcao/module_operator_lcao/overlap.h"
@@ -429,6 +430,29 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
         }
     }
 
+    // atomic force and stress for DeltaP
+    ModuleBase::matrix force_deltap;
+    ModuleBase::matrix stress_deltap;
+    if (PARAM.inp.deltap_switch && PARAM.inp.deltap_corr)
+    {
+        const auto& dp_lambda = hamilt::DeltaPOperator<std::complex<double>, double>::get_stored_lambda();
+        if (!dp_lambda.empty())
+        {
+            if (isforce) force_deltap.create(nat, 3);
+            if (isstress) stress_deltap.create(3, 3);
+
+            double rm = PARAM.inp.deltap_rm > 0.0 ? PARAM.inp.deltap_rm : 3.0;
+            hamilt::DeltaPOperator<std::complex<double>, double> tmp_dp(
+                nullptr, kv.kvec_d, nullptr, ucell, &gd,
+                two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), rm);
+            tmp_dp.set_lambda(dp_lambda);
+            tmp_dp.set_gdir(PARAM.inp.deltap_gdir);
+
+            const hamilt::HContainer<double>* dmr = dmat.dm->get_DMR_pointer(1);
+            tmp_dp.cal_force_stress(isforce, isstress, dmr, force_deltap, stress_deltap);
+        }
+    }
+
     // NOTE: finish_ftable is no longer needed as we don't use ForceStressArrays for overlap/kinetic
     // if (!PARAM.globalv.gamma_only_local)
     // {
@@ -498,6 +522,10 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                 {
                     fcs(iat, i) += force_dspin(iat, i);
                 }
+                if (PARAM.inp.deltap_switch && PARAM.inp.deltap_corr)
+                {
+                    fcs(iat, i) += force_deltap(iat, i);
+                }
 #ifdef __EXX
                 // Force contribution from exx
                 if (GlobalC::exx_info.info_global.cal_exx)
@@ -709,6 +737,10 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                 {
                     scs(i, j) += stress_dspin(i, j);
                 }
+                if (PARAM.inp.deltap_switch && PARAM.inp.deltap_corr)
+                {
+                    scs(i, j) += stress_deltap(i, j);
+                }
 #ifdef __EXX
                 // Stress contribution from exx
                 if (GlobalC::exx_info.info_global.cal_exx)