Fix: lack of switch_dmr

Author: dyzheng

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>
@@ -627,12 +628,12 @@ void DensityMatrix<std::complex<double>, double>::cal_DMR_full(hamilt::HContaine
 #endif
             // loop over k-points
             // calculate full matrix for complex density matrix
-            for (int ik = 0; ik < this->_nks; ++ik)
+            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->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI;
+                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);
@@ -983,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

@@ -196,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
@@ -275,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/FORCE_STRESS.cpp

@@ -171,25 +171,26 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                         pv,
                         kv);
     // calculate force and stress for Nonlocal part
-    if(GlobalV::NSPIN!=4)
+    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(GlobalV::NSPIN==2)
+        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(GlobalV::NSPIN==2)
+        if(PARAM.inp.nspin == 2)
         {
             const_cast<elecstate::DensityMatrix<T, double>*>(dm_p)->switch_dmr(0);
         }
@@ -201,6 +202,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
                     kv.kvec_d,
                     nullptr,
                     &GlobalC::ucell,
+                    orb.cutoffs(),
                     &GlobalC::GridD,
                     two_center_bundle.overlap_orb_beta.get()
             );

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_force_stress.hpp

@@ -47,12 +47,11 @@ void NonlocalNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
             const int iat1 = ucell->itia2iat(T1, I1);
             const ModuleBase::Vector3<int>& R_index1 = adjs.box[ad];
             // choose the real adjacent atoms
-            const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance();
             // Note: the distance of atoms should less than the cutoff radius, 
             // When equal, the theoretical value of matrix element is zero, 
             // but the calculated value is not zero due to the numerical error, which would lead to result changes.
             if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0
-                < orb.Phi[T1].getRcut() + this->ucell->infoNL.Beta[this->current_type].get_rcut_max())
+                < orb_cutoff_[T1] + this->ucell->infoNL.Beta[this->current_type].get_rcut_max())
             {
                 is_adj[ad] = true;
             }
@@ -68,7 +67,6 @@ void NonlocalNew<OperatorLCAO<TK, TR>>::cal_force_stress(const bool cal_force,
             const ModuleBase::Vector3<double>& tau1 = adjs.adjacent_tau[ad];
             const Atom* atom1 = &ucell->atoms[T1];
 
-            const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance();
             auto all_indexes = paraV->get_indexes_row(iat1);
             auto col_indexes = paraV->get_indexes_col(iat1);
             // insert col_indexes into all_indexes to get universal set with no repeat elements