Fix&Refactor: soc force/stress error for LCAO code, refactor the force/stress code into operator

source/Makefile.Objects

@@ -571,8 +571,6 @@ OBJS_LCAO=evolve_elec.o\
       stress_tools.o\
       edm.o\
       pulay_force_stress_center2.o\
-	  fvnl_dbeta_gamma.o\
-	  fvnl_dbeta_k.o\
 	  grid_init.o\
       spar_dh.o\
       spar_exx.o\

source/module_elecstate/module_dm/density_matrix.cpp

@@ -22,6 +22,7 @@ DensityMatrix<TK, TR>::~DensityMatrix()
     {
         delete it;
     }
+    delete[] this->dmr_tmp_;
 }
 
 template <typename TK, typename TR>
@@ -589,6 +590,78 @@ 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];
+        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->_nk; ++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->_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)
@@ -911,6 +984,82 @@ void DensityMatrix<std::complex<double>, double>::write_DMK(const std::string di
     ofs.close();
 }
 
+// switch_dmr
+template <typename TK, typename TR>
+void DensityMatrix<TK, TR>::switch_dmr(const int mode)
+{
+    ModuleBase::TITLE("DensityMatrix", "switch_dmr");
+    if (this->_nspin != 2)
+    {
+        return;
+    }
+    else
+    {
+        ModuleBase::timer::tick("DensityMatrix", "switch_dmr");
+        switch(mode)
+        {
+        case 0:
+            // switch to original density matrix
+            if (this->dmr_tmp_ != nullptr && this->dmr_origin_.size() != 0) 
+            {
+                this->_DMR[0]->allocate(this->dmr_origin_.data(), false);
+                delete[] this->dmr_tmp_;
+                this->dmr_tmp_ = nullptr;
+            }
+            // else: do nothing
+            break;
+        case 1:
+            // switch to total magnetization density matrix, dmr_up + dmr_down
+            if(this->dmr_tmp_ == nullptr)
+            {
+                const size_t size = this->_DMR[0]->get_nnr();
+                this->dmr_tmp_ = new TR[size];
+                this->dmr_origin_.resize(size);
+                for (int i = 0; i < size; ++i)
+                {
+                    this->dmr_tmp_[i] = this->dmr_origin_[i] + this->_DMR[1]->get_wrapper()[i];
+                }
+                this->_DMR[0]->allocate(this->dmr_tmp_, false);
+            }
+            else
+            {
+                const size_t size = this->_DMR[0]->get_nnr();
+                for (int i = 0; i < size; ++i)
+                {
+                    this->dmr_tmp_[i] = this->dmr_origin_[i] - this->_DMR[1]->get_wrapper()[i];
+                }
+            }
+            break;
+        case 2:
+            // switch to magnetization density matrix, dmr_up - dmr_down
+            if(this->dmr_tmp_ == nullptr)
+            {
+                const size_t size = this->_DMR[0]->get_nnr();
+                this->dmr_tmp_ = new TR[size];
+                this->dmr_origin_.resize(size);
+                for (int i = 0; i < size; ++i)
+                {
+                    this->dmr_origin_[i] = this->_DMR[0]->get_wrapper()[i];
+                    this->dmr_tmp_[i] = this->dmr_origin_[i] - this->_DMR[1]->get_wrapper()[i];
+                }
+                this->_DMR[0]->allocate(this->dmr_tmp_, false);
+            }
+            else
+            {
+                const size_t size = this->_DMR[0]->get_nnr();
+                for (int i = 0; i < size; ++i)
+                {
+                    this->dmr_tmp_[i] = this->dmr_origin_[i] - this->_DMR[1]->get_wrapper()[i];
+                }
+            }
+            break;
+        default:
+            throw std::string("Unknown mode in switch_dmr");
+        }
+        ModuleBase::timer::tick("DensityMatrix", "switch_dmr");
+    }
+}
+
 // T of HContainer can be double or complex<double>
 template class DensityMatrix<double, double>;               // Gamma-Only case
 template class DensityMatrix<std::complex<double>, double>; // Multi-k case

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
      */
@@ -194,6 +196,12 @@ class DensityMatrix
      */
     void sum_DMR_spin();
 
+    /**
+     * @brief (Only nspin=2) switch DMR to total density matrix or magnetization density matrix
+     * @param mode 0 - original density matrix; 1 - total density matrix; 2 - magnetization density matrix
+     */
+    void switch_dmr(const int mode);
+
     /**
      * @brief write density matrix dm(ik) into *.dmk
      * @param directory directory of *.dmk files
@@ -273,6 +281,9 @@ class DensityMatrix
      */
     int _nk = 0;
 
+    /// temporary pointers for switch DMR, only used with nspin=2
+    std::vector<TR> dmr_origin_;
+    TR* dmr_tmp_ = nullptr;
 
 };
 

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -21,8 +21,6 @@ if(ENABLE_LCAO)
         FORCE_k.cpp
         stress_tools.cpp
         edm.cpp
-        fvnl_dbeta_gamma.cpp
-        fvnl_dbeta_k.cpp
         grid_init.cpp
         spar_dh.cpp
         spar_exx.cpp

source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h

@@ -116,17 +116,6 @@ class Force_LCAO
                         ModuleBase::matrix& stvnl_dphi,
                         Record_adj* ra = nullptr);
 
-    void cal_fvnl_dbeta(const elecstate::DensityMatrix<T, double>* dm,
-                        const Parallel_Orbitals& pv,
-                        const UnitCell& ucell,
-                        const LCAO_Orbitals& orb,
-                        const TwoCenterIntegrator& intor_orb_beta,
-                        Grid_Driver& gd,
-                        const bool isforce,
-                        const bool isstress,
-                        ModuleBase::matrix& fvnl_dbeta,
-                        ModuleBase::matrix& svnl_dbeta);
-
     //-------------------------------------------
     // forces related to local pseudopotentials
     //-------------------------------------------

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,52 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                         orb,
                         pv,
                         kv);
+    // calculate force and stress for Nonlocal part
+    if(PARAM.inp.nspin != 4)
+    {
+        hamilt::NonlocalNew<hamilt::OperatorLCAO<T, double>> tmp_nonlocal(
+                    nullptr,
+                    kv.kvec_d,
+                    nullptr,
+                    &GlobalC::ucell,
+                    orb.cutoffs(),
+                    &GlobalC::GridD,
+                    two_center_bundle.overlap_orb_beta.get()
+            );
+
+        const auto* dm_p = dynamic_cast<const elecstate::ElecStateLCAO<T>*>(pelec)->get_DM();
+        if(PARAM.inp.nspin == 2)
+        {
+            const_cast<elecstate::DensityMatrix<T, double>*>(dm_p)->switch_dmr(1);
+        }
+        const hamilt::HContainer<double>* dmr = dm_p->get_DMR_pointer(1);
+        tmp_nonlocal.cal_force_stress(isforce, isstress, dmr, fvnl_dbeta, svnl_dbeta);
+        if(PARAM.inp.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,
+                    orb.cutoffs(),
+                    &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,17 +219,6 @@ 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,
-                         pv,
-                         ucell,
-                         orb,
-                         *(two_center_bundle.overlap_orb_beta),
-                         GlobalC::GridD,
-                         isforce,
-                         isstress,
-                         fvnl_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,17 +319,6 @@ 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,
-                         pv,
-                         ucell,
-                         orb,
-                         *(two_center_bundle.overlap_orb_beta),
-                         GlobalC::GridD,
-                         isforce,
-                         isstress,
-                         fvnl_dbeta,
-                         svnl_dbeta);
-
 #ifdef __DEEPKS
     if (PARAM.inp.deepks_scf)
     {