Add TDDFT velocity gauge force implementation and test

- Add cal_force() methods to TDEkinetic and TDNonlocal classes
- Create td_ekinetic_force.hpp and td_nonlocal_force.hpp templates
- Fix force consistency: use Nonlocal when |A| < 1e-10
- Add test case comparing DFT vs TDDFT(A=0) forces
- Results: forces match exactly when A=0 (4.968 eV/A for CO)

Test: CO molecule, PBE/DZP, forces verified consistent

Author: dyzheng

source/source_lcao/FORCE_STRESS.cpp

@@ -22,8 +22,23 @@
 #include "source_lcao/module_operator_lcao/nonlocal.h"
 #include "source_lcao/module_operator_lcao/ekinetic.h"
 #include "source_lcao/module_operator_lcao/overlap.h"
+#include "source_lcao/module_operator_lcao/td_ekinetic_lcao.h"
+#include "source_lcao/module_operator_lcao/td_nonlocal_lcao.h"
+#include "source_lcao/module_rt/td_info.h"
 #include "source_lcao/pulay_fs.h"
-
+void print_force(ModuleBase::matrix& m, std::string name, int nat)
+{
+    for(int i = 0; i<nat; i++)
+    {
+        std::cout<< name <<std::endl;
+        std::cout<<"nat: "<<i<<std::endl;
+        for(int j = 0; j<3; j++)
+        {
+            std::cout<<m(i,j)<<" ";
+        }
+        std::cout<<std::endl;
+    }
+}
 
 // mohan add 2025-11-04
 template <>
@@ -111,6 +126,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
     ModuleBase::matrix fcc;
     ModuleBase::matrix fscc;
     ModuleBase::matrix fvnl_dalpha; // deepks
+    ModuleBase::matrix fekinetic_td;
 
     fvl_dphi.create(nat, 3); // must do it now, update it later, noted by zhengdy
 
@@ -125,6 +141,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
         fcc.create(nat, 3); // force due to core correction
         fscc.create(nat, 3); // force due to self-consistent field
         fvnl_dalpha.create(nat, 3); // deepks
+        fekinetic_td.create(nat, 3);
 
         // calculate basic terms in Force, same method with PW base
         this->calForcePwPart(ucell, fvl_dvl, fewalds, fcc, fscc, pelec->f_en.etxc,
@@ -200,11 +217,31 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
         tmp_overlap.cal_force_stress(isforce, isstress, edmR, foverlap, soverlap);
 
         // Calculate nonlocal force/stress (uses DM)
-        hamilt::Nonlocal<hamilt::OperatorLCAO<T, double>> tmp_nonlocal(
+        // When A=0, TDNonlocal should give same result as Nonlocal
+        // Use Nonlocal when |A| is very small to ensure consistency
+        const double A_norm = TD_info::cart_At.norm();
+        if(PARAM.inp.td_stype == 1 && A_norm > 1e-10)
+        {
+            hamilt::TDNonlocal<hamilt::OperatorLCAO<std::complex<double>, double>> tmp_nonlocal(
+            nullptr, kv.kvec_d, nullptr, &ucell, orb, &gd);
+            tmp_nonlocal.cal_force(isforce, dmR, fvnl_dbeta);
+        }
+        else
+        {
+            hamilt::Nonlocal<hamilt::OperatorLCAO<T, double>> tmp_nonlocal(
             nullptr, kv.kvec_d, nullptr, &ucell, orb.cutoffs(), &gd,
             two_center_bundle.overlap_orb_beta.get());
-        tmp_nonlocal.cal_force_stress(isforce, isstress, dmR, fvnl_dbeta, svnl_dbeta);
+            tmp_nonlocal.cal_force_stress(isforce, isstress, dmR, fvnl_dbeta, svnl_dbeta);
+        }
 
+        if(PARAM.inp.td_stype == 1)
+        {
+            hamilt::TDEkinetic<hamilt::OperatorLCAO<std::complex<double>, double>> tmp_td_ekinetic(
+                nullptr, nullptr, &kv, &ucell, orb.cutoffs(), &gd, 
+                two_center_bundle.overlap_orb.get());
+            tmp_td_ekinetic.cal_force(isforce, &pv, psi, pelec, fekinetic_td);
+        }
+
         // Switch back to spin channel 0
         if (PARAM.inp.nspin == 2)
         {
@@ -472,6 +509,22 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
     //--------------------------------
     if (isforce)
     {
+        if(GlobalV::MY_RANK==0)
+        {
+            print_force(foverlap, "foverlap",nat);
+            print_force(ftvnl_dphi, "ftvnl_dphi",nat);
+            print_force(fvnl_dbeta, "fvnl_dbeta",nat);
+            print_force(fvl_dphi, "fvl_dphi",nat);
+            print_force(fvl_dvl, "fvl_dvl",nat);
+            print_force(fewalds, "fewalds",nat);
+            print_force(fcc, "fcc",nat);
+            print_force(fscc, "fscc",nat);
+            // if (PARAM.inp.esolver_type == "tddft")
+            // {
+            //     print_force(fpothybrid, "fpothybrid",nat);
+            //     print_force(fefield_tddft, "fefield_tddft",nat);
+            // }    
+        }
         //---------------------------------
         // sum all parts of force!
         //---------------------------------
@@ -517,6 +570,10 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                 if (PARAM.inp.esolver_type == "tddft")
                 {
                     fcs(iat, i) += fefield_tddft(iat, i);
+                    if(PARAM.inp.td_stype == 1)
+                    {
+                        fcs(iat, i) += fekinetic_td(iat, i);
+                    }
                 }
                 // Gate field force
                 if (PARAM.inp.gate_flag)

source/source_lcao/module_operator_lcao/td_ekinetic_force.hpp

@@ -0,0 +1,161 @@
+#pragma once
+#include "td_ekinetic_lcao.h"
+#include "source_base/parallel_reduce.h"
+#include "source_base/timer.h"
+#include "source_base/libm/libm.h"
+#include "source_estate/module_dm/cal_dm_psi.h"
+namespace hamilt
+{
+template <typename TK, typename TR>
+void TDEkinetic<OperatorLCAO<TK, TR>>::cal_force(
+    const bool cal_force,
+    const Parallel_Orbitals* paraV,
+    const psi::Psi<std::complex<double>>* psi,
+    const elecstate::ElecState* pelec,
+    ModuleBase::matrix& force)
+{
+    const int npol = ucell->get_npol();
+    ModuleBase::Vector3<double> At = TD_info::cart_At;
+    if (cal_force)
+    {
+        force.zero_out();
+    }
+    else return;
+    // calculate dmr
+    const int nspin0 = PARAM.inp.nspin;
+    const int nspin_dm = std::map<int, int>({ {1,1},{2,2},{4,1} })[nspin0];
+    elecstate::DensityMatrix<std::complex<double>, double> tmp_dm(paraV, nspin_dm, kv->kvec_d, kv->get_nks() / nspin_dm);
+    elecstate::cal_dm_psi(paraV, pelec->wg, psi[0], tmp_dm);
+    tmp_dm.init_DMR(Grid, ucell);
+    tmp_dm.cal_DMR();
+    if (PARAM.inp.nspin == 2)
+    {
+        tmp_dm.switch_dmr(1);
+    }
+    // Convert to complex HContainer for TDDFT velocity gauge
+    hamilt::HContainer<std::complex<double>> dmR_complex(paraV);
+    tmp_dm.cal_DMR_full(&dmR_complex);
+    const hamilt::HContainer<std::complex<double>>* dmR = &dmR_complex;
+    // Loop over all atom pairs and calculate force contributions
+    #pragma omp parallel
+    {
+        ModuleBase::matrix force_local(force.nr, force.nc);
+
+        #pragma omp for schedule(dynamic)
+        for (int iat1 = 0; iat1 < ucell->nat; iat1++)
+        {
+            auto tau1 = ucell->get_tau(iat1);
+            int T1 = 0, I1 = 0;
+            ucell->iat2iait(iat1, &I1, &T1);
+            Atom& atom1 = ucell->atoms[T1];
+
+            // Find adjacent atoms
+            AdjacentAtomInfo adjs;
+            Grid->Find_atom(*ucell, tau1, T1, I1, &adjs);
+
+            for (int ad = 0; ad < adjs.adj_num + 1; ++ad)
+            {
+                const int T2 = adjs.ntype[ad];
+                const int I2 = adjs.natom[ad];
+                const int iat2 = ucell->itia2iat(T2, I2);
+                const ModuleBase::Vector3<int>& R_index = adjs.box[ad];
+
+                // Check cutoff
+                ModuleBase::Vector3<double> dtau = ucell->cal_dtau(iat1, iat2, R_index);
+                if (dtau.norm() * ucell->lat0 >= orb_cutoff_[T1] + orb_cutoff_[T2])
+                {
+                    continue;
+                }
+
+                // Find density matrix for this atom pair
+                const hamilt::BaseMatrix<std::complex<double>>* dm_matrix = dmR->find_matrix(iat1, iat2, R_index[0], R_index[1], R_index[2]);
+                if (dm_matrix == nullptr)
+                {
+                    continue;
+                }
+
+                // Calculate force for this atom pair
+                double* force_tmp1 = (cal_force) ? &force_local(iat1, 0) : nullptr;
+                double* force_tmp2 = (cal_force) ? &force_local(iat2, 0) : nullptr;
+
+                Atom& atom2 = ucell->atoms[T2];
+                auto row_indexes = paraV->get_indexes_row(iat1);
+                auto col_indexes = paraV->get_indexes_col(iat2);
+
+                if (row_indexes.size() == 0 || col_indexes.size() == 0)
+                {
+                    continue;
+                }
+                const std::complex<double>* dm_pointer = dm_matrix->get_pointer();
+                double overlap = 0;
+                double grad[3] = {0, 0, 0};
+                double hess[9] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
+
+                // step_trace = 0 for npol=1; ={0, 1, col_size, col_size+1} for npol=2
+                std::vector<int> step_trace(npol * npol, 0);
+                if (npol == 2)
+                {
+                    step_trace[1] = 1;
+                    step_trace[2] = col_indexes.size();
+                    step_trace[3] = col_indexes.size() + 1;
+                }
+
+                // Loop over orbital pairs
+                for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol)
+                {
+                    const int iw1 = row_indexes[iw1l] / npol;
+                    const int L1 = atom1.iw2l[iw1];
+                    const int N1 = atom1.iw2n[iw1];
+                    const int m1 = atom1.iw2m[iw1];
+                    const int M1 = (m1 % 2 == 0) ? -m1 / 2 : (m1 + 1) / 2;
+
+                    for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol)
+                    {
+                        const int iw2 = col_indexes[iw2l] / npol;
+                        const int L2 = atom2.iw2l[iw2];
+                        const int N2 = atom2.iw2n[iw2];
+                        const int m2 = atom2.iw2m[iw2];
+                        const int M2 = (m2 % 2 == 0) ? -m2 / 2 : (m2 + 1) / 2;
+
+                        // Calculate integral and its gradient using provided functor
+                        intor_->calculate(T1, L1, N1, M1, T2, L2, N2, M2, dtau * this->ucell->lat0, &overlap, grad, hess);
+
+                        // Calculate force contribution with compile-time sign
+                        if (cal_force)
+                        {
+                            // Factor of 2 for Hermitian matrix will be applied later
+                            for (int i = 0; i < 3; i++)
+                            {
+                                // force_tmp1[i] += 2 * dm_pointer[0] * grad[i] * At * At;
+                                force_tmp2[i] -= dm_pointer[0].real() * grad[i] * At * At;
+                                for (int j = 0; j < 3; j++)
+                                {
+                                    force_tmp2[i] -= 2 * dm_pointer[0].imag() * hess[j * 3 + i] * At[j] * At[j];
+                                }
+                            }
+                        }
+                        dm_pointer += npol;
+                    }
+                    dm_pointer += (npol - 1) * col_indexes.size();
+                }
+            }
+        }
+
+        #pragma omp critical
+        {
+            if (cal_force)
+            {
+                force += force_local;
+            }
+        }
+    }
+
+    // Finalize with MPI reduction and post-processing
+    if (cal_force)
+    {
+#ifdef __MPI
+        Parallel_Reduce::reduce_all(force.c, force.nr * force.nc);
+#endif
+    }
+}
+}// namespace hamilt

source/source_lcao/module_operator_lcao/td_ekinetic_lcao.cpp

@@ -10,6 +10,7 @@
 #include "source_lcao/center2_orb-orb11.h"
 #include "source_lcao/spar_hsr.h"
 #include "source_lcao/module_hcontainer/hcontainer_funcs.h"
+#include "td_ekinetic_force.hpp"
 
 namespace hamilt
 {
@@ -27,6 +28,7 @@ TDEkinetic<OperatorLCAO<TK, TR>>::TDEkinetic(HS_Matrix_K<TK>* hsk_in,
     this->cal_type = calculation_type::lcao_tddft_periodic;
     this->Grid = GridD_in;
     // initialize HR to get adjs info.
+    if(hR_in!=nullptr)
     this->initialize_HR(Grid);
 }
 

source/source_lcao/module_operator_lcao/td_ekinetic_lcao.h

@@ -4,6 +4,7 @@
 #include "source_basis/module_nao/two_center_integrator.h"
 #include "source_cell/klist.h"
 #include "source_cell/module_neighbor/sltk_grid_driver.h"
+#include "source_estate/elecstate.h"
 #include "source_lcao/module_hcontainer/hcontainer.h"
 #include "source_lcao/module_rt/td_info.h"
 #include "operator_lcao.h"
@@ -81,6 +82,23 @@ class TDEkinetic<OperatorLCAO<TK,TR>> : public OperatorLCAO<TK, TR>
      */
     void calculate_HR();
 
+            /**
+     * @brief calculate force and stress for kinetic operator
+     * @param cal_force whether to calculate force
+     * @param dmR density matrix in real space
+     * @param force output force matrix (nat x 3)
+     */
+    void cal_force(const bool cal_force,
+                   const Parallel_Orbitals* paraV,
+                   const psi::Psi<std::complex<double>>* psi,
+                   const elecstate::ElecState* pelec,
+                   ModuleBase::matrix& force);
+    void cal_force(const bool cal_force,
+                   const Parallel_Orbitals* paraV,
+                   const psi::Psi<double>* psi,
+                   const elecstate::ElecState* pelec,
+                   ModuleBase::matrix& force){return;}
+
     virtual void set_HR_fixed(void*) override;