Refactor&Fix: soc force&stress for LCAO

Author: dyzheng

source/module_elecstate/module_dm/density_matrix.cpp

@@ -589,6 +589,82 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR()
     ModuleBase::timer::tick("DensityMatrix", "cal_DMR");
 }
 
+// calculate DMR from DMK using blas for multi-k calculation
+template <>
+void DensityMatrix<double, double>::cal_DMR_full(hamilt::HContainer<std::complex<double>>* dmR_out)const{}
+template <>
+void DensityMatrix<std::complex<double>, double>::cal_DMR_full(hamilt::HContainer<std::complex<double>>* dmR_out)const
+{
+    ModuleBase::TITLE("DensityMatrix", "cal_DMR_full");
+
+    ModuleBase::timer::tick("DensityMatrix", "cal_DMR_full");
+    int ld_hk = this->_paraV->nrow;
+    int ld_hk2 = 2 * ld_hk;
+    hamilt::HContainer<std::complex<double>>* tmp_DMR = dmR_out;
+    // set zero since this function is called in every scf step
+    tmp_DMR->set_zero();
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+    for (int i = 0; i < tmp_DMR->size_atom_pairs(); ++i)
+    {
+        auto& tmp_ap = tmp_DMR->get_atom_pair(i);
+        int iat1 = tmp_ap.get_atom_i();
+        int iat2 = tmp_ap.get_atom_j();
+        // get global indexes of whole matrix for each atom in this process
+        int row_ap = this->_paraV->atom_begin_row[iat1];
+        int col_ap = this->_paraV->atom_begin_col[iat2];
+        if (row_ap == -1 || col_ap == -1)
+        {
+            throw std::string("Atom-pair not belong this process");
+        }
+        for (int ir = 0; ir < tmp_ap.get_R_size(); ++ir)
+        {
+            const ModuleBase::Vector3<int> r_index = tmp_ap.get_R_index(ir);
+            auto* tmp_matrix = tmp_ap.find_matrix(r_index);
+#ifdef __DEBUG
+            if (tmp_matrix == nullptr)
+            {
+                std::cout << "tmp_matrix is nullptr" << std::endl;
+                continue;
+            }
+#endif
+            // loop over k-points
+            // calculate full matrix for complex density matrix
+            for (int ik = 0; ik < this->_nks; ++ik)
+            {
+                // cal k_phase
+                // if TK==std::complex<double>, kphase is e^{ikR}
+                const ModuleBase::Vector3<double> dR(r_index[0], r_index[1], r_index[2]);
+                const double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                double sinp, cosp;
+                ModuleBase::libm::sincos(arg, &sinp, &cosp);
+                std::complex<double> kphase = std::complex<double>(cosp, sinp);
+                // set DMR element
+                std::complex<double>* tmp_DMR_pointer = tmp_matrix->get_pointer();
+                const std::complex<double>* tmp_DMK_pointer = this->_DMK[ik].data();
+                double* DMK_real_pointer = nullptr;
+                double* DMK_imag_pointer = nullptr;
+                // jump DMK to fill DMR
+                // DMR is row-major, DMK is column-major
+                tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap;
+                for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu)
+                {
+                    BlasConnector::axpy(this->_paraV->get_col_size(iat2),
+                                        kphase,
+                                        tmp_DMK_pointer,
+                                        ld_hk,
+                                        tmp_DMR_pointer,
+                                        1);
+                    tmp_DMK_pointer += 1;
+                    tmp_DMR_pointer += this->_paraV->get_col_size(iat2);
+                }
+            }
+        }
+    }
+    ModuleBase::timer::tick("DensityMatrix", "cal_DMR_full");
+}
+
 // calculate DMR from DMK using blas for multi-k calculation, without summing over k-points
 template <>
 void DensityMatrix<std::complex<double>, double>::cal_DMR(const int ik)

source/module_elecstate/module_dm/density_matrix.h

@@ -179,6 +179,8 @@ class DensityMatrix
      */
     void cal_DMR();
 
+    void cal_DMR_full(hamilt::HContainer<std::complex<double>>* dmR_out) const;
+
     /**
      * @brief calculate density matrix DMR from dm(k) using blas::axpy for multi-k calculation, without summing over k-points
      */

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -18,6 +18,8 @@
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks_io.h" // mohan add 2024-07-22 
 #endif
 #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h"
+#include "module_elecstate/elecstate_lcao.h"
 
 template <typename T>
 Force_Stress_LCAO<T>::Force_Stress_LCAO(Record_adj& ra, const int nat_in) : RA(&ra), f_pw(nat_in), nat(nat_in)
@@ -168,6 +170,50 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                         orb,
                         pv,
                         kv);
+    // calculate force and stress for Nonlocal part
+    if(GlobalV::NSPIN!=4)
+    {
+        hamilt::NonlocalNew<hamilt::OperatorLCAO<T, double>> tmp_nonlocal(
+                    nullptr,
+                    kv.kvec_d,
+                    nullptr,
+                    &GlobalC::ucell,
+                    &GlobalC::GridD,
+                    two_center_bundle.overlap_orb_beta.get()
+            );
+
+        const auto* dm_p = dynamic_cast<const elecstate::ElecStateLCAO<T>*>(pelec)->get_DM();
+        if(GlobalV::NSPIN==2)
+        {
+            const_cast<elecstate::DensityMatrix<T, double>*>(dm_p)->switch_dmr(2);
+        }
+        const hamilt::HContainer<double>* dmr = dm_p->get_DMR_pointer(1);
+        tmp_nonlocal.cal_force_stress(isforce, isstress, dmr, fvnl_dbeta, svnl_dbeta);
+        if(GlobalV::NSPIN==2)
+        {
+            const_cast<elecstate::DensityMatrix<T, double>*>(dm_p)->switch_dmr(0);
+        }
+    }
+    else
+    {
+        hamilt::NonlocalNew<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>> tmp_nonlocal(
+                    nullptr,
+                    kv.kvec_d,
+                    nullptr,
+                    &GlobalC::ucell,
+                    &GlobalC::GridD,
+                    two_center_bundle.overlap_orb_beta.get()
+            );
+
+        const auto* dm_p = dynamic_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(pelec)->get_DM();
+        hamilt::HContainer<std::complex<double>> tmp_dmr(dm_p->get_DMR_pointer(1)->get_paraV());
+        std::vector<int> ijrs = dm_p->get_DMR_pointer(1)->get_ijr_info();
+        tmp_dmr.insert_ijrs(&ijrs);
+        tmp_dmr.allocate();
+        dm_p->cal_DMR_full(&tmp_dmr);
+        tmp_nonlocal.cal_force_stress(isforce, isstress, &tmp_dmr, fvnl_dbeta, svnl_dbeta);
+    }
+    
 
     //! forces and stress from vdw
     //  Peize Lin add 2014-04-04, update 2021-03-09

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp

@@ -106,15 +106,6 @@ void Force_LCAO<double>::allocate(const Parallel_Orbitals& pv,
                               &GlobalC::GridD,
                               nullptr);
 
-    LCAO_domain::build_Nonlocal_mu_new(pv,
-                                       fsr,
-                                       nullptr,
-                                       cal_deri,
-                                       GlobalC::ucell,
-                                       orb,
-                                       *(two_center_bundle.overlap_orb_beta),
-                                       &GlobalC::GridD);
-
     // calculate asynchronous S matrix to output for Hefei-NAMD
     if (PARAM.inp.cal_syns)
     {
@@ -228,7 +219,7 @@ void Force_LCAO<double>::ftable(const bool isforce,
     //tvnl_dphi
     PulayForceStress::cal_pulay_fs(ftvnl_dphi, stvnl_dphi, *dm, ucell, pv, dHx, dHxy, isforce, isstress);
 
-    this->cal_fvnl_dbeta(dm,
+    /*this->cal_fvnl_dbeta(dm,
                          pv,
                          ucell,
                          orb,
@@ -237,7 +228,7 @@ void Force_LCAO<double>::ftable(const bool isforce,
                          isforce,
                          isstress,
                          fvnl_dbeta,
-                         svnl_dbeta);
+                         svnl_dbeta);*/
 
     // vl_dphi
     PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/);

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp

@@ -131,16 +131,6 @@ void Force_LCAO<std::complex<double>>::allocate(const Parallel_Orbitals& pv,
                               &GlobalC::GridD,
                               nullptr); // delete lm.Hloc_fixedR
 
-    // calculate dVnl=<phi|dVnl|dphi> in LCAO
-    LCAO_domain::build_Nonlocal_mu_new(pv,
-                                       fsr,
-                                       nullptr,
-                                       cal_deri,
-                                       GlobalC::ucell,
-                                       orb,
-                                       *(two_center_bundle.overlap_orb_beta),
-                                       &GlobalC::GridD);
-
     // calculate asynchronous S matrix to output for Hefei-NAMD
     if (PARAM.inp.cal_syns)
     {
@@ -329,7 +319,7 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
     // vl_dphi
     PulayForceStress::cal_pulay_fs(fvl_dphi, svl_dphi, *dm, ucell, pelec->pot, gint, isforce, isstress, false/*reset dm to gint*/);
 
-    this->cal_fvnl_dbeta(dm,
+    /*this->cal_fvnl_dbeta(dm,
                          pv,
                          ucell,
                          orb,
@@ -338,7 +328,7 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
                          isforce,
                          isstress,
                          fvnl_dbeta,
-                         svnl_dbeta);
+                         svnl_dbeta);*/
 
 #ifdef __DEEPKS
     if (PARAM.inp.deepks_scf)