Test: add unit test of DFTUNew and refactor some code (#3814)

* Fix: not use static member in DFTUNew

* Test: add unit test of DFTUNew

* Test: add nspin=2 test and delete GlobalV::CURRENT_SPIN in DFTUNew

* add annotation for dftu.h

* Refactor: new dftu code with suggestions by developers

---------

Co-authored-by: dyzheng <[email protected]>

source/Makefile.Objects

@@ -274,7 +274,7 @@ OBJS_HAMILT_LCAO=hamilt_lcao.o\
     deepks_lcao.o\
     op_exx_lcao.o\
     sc_lambda_lcao.o\
-    dftu_new.o\
+    dftu_lcao.o\
 
 OBJS_HCONTAINER=base_matrix.o\
     atom_pair.o\

source/module_esolver/esolver_ks_lcao.cpp

@@ -18,7 +18,6 @@
 #include "module_elecstate/module_charge/symmetry_rho.h"
 #include "module_elecstate/occupy.h"
 #include "module_hamilt_lcao/module_dftu/dftu.h"
-#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/print_info.h"
 #ifdef __EXX
@@ -641,15 +640,9 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
 
     if (GlobalV::dft_plus_u)
     {
-        if(istep == 0 && iter == 1)
+        if(istep != 0 || iter != 1)
         {
-            hamilt::DFTUNew<hamilt::OperatorLCAO<TK, TR>>::dm_in_dftu = nullptr;
-        }
-        else
-        {
-            hamilt::DFTUNew<hamilt::OperatorLCAO<TK, TR>>::dm_in_dftu =
-            dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)
-                ->get_DM();
+            GlobalC::dftu.set_dmr( dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM() );
         }
         // Calculate U and J if Yukawa potential is used
         GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx);

source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt

@@ -14,7 +14,7 @@ if(ENABLE_LCAO)
         operator_lcao/td_ekinetic_lcao.cpp
         operator_lcao/td_nonlocal_lcao.cpp
         operator_lcao/sc_lambda_lcao.cpp
-        operator_lcao/dftu_new.cpp
+        operator_lcao/dftu_lcao.cpp
         FORCE_STRESS.cpp
         FORCE_gamma.cpp
         FORCE_gamma_edm.cpp

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -14,7 +14,7 @@
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks 2021-06-03
 #include "module_elecstate/elecstate_lcao.h"
 #endif
-#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h"
+#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_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)
@@ -240,14 +240,14 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
         }
         else
         {
-            hamilt::DFTUNew<hamilt::OperatorLCAO<T, double>> tmp_dftu(uhm.LM,
+            hamilt::DFTU<hamilt::OperatorLCAO<T, double>> tmp_dftu(uhm.LM,
                                                                  kv.kvec_d,
                                                                  nullptr,
                                                                  nullptr,
-                                                                 &GlobalC::ucell,
+                                                                 GlobalC::ucell,
                                                                  &GlobalC::GridD,
                                                                  &GlobalC::dftu,
-                                                                 uhm.LM->ParaV);
+                                                                 *(uhm.LM->ParaV));
             tmp_dftu.cal_force_stress(isforce, isstress, force_dftu, stress_dftu);
         }
     }

source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp

@@ -17,7 +17,7 @@
 #include "module_hsolver/diago_elpa.h"
 #endif
 #include "operator_lcao/op_dftu_lcao.h"
-#include "operator_lcao/dftu_new.h"
+#include "operator_lcao/dftu_lcao.h"
 #include "operator_lcao/meta_lcao.h"
 #include "operator_lcao/op_exx_lcao.h"
 #include "operator_lcao/td_ekinetic_lcao.h"
@@ -226,15 +226,15 @@ HamiltLCAO<TK, TR>::HamiltLCAO(
             }
             else
             {
-                dftu = new DFTUNew<OperatorLCAO<TK, TR>>(
+                dftu = new DFTU<OperatorLCAO<TK, TR>>(
                     LM_in,
                     this->kv->kvec_d,
                     this->hR,
                     &(this->getHk(LM_in)),
-                    &GlobalC::ucell,
+                    GlobalC::ucell,
                     &GlobalC::GridD,
                     &GlobalC::dftu,
-                    LM_in->ParaV
+                    *(LM_in->ParaV)
                 );
             }
             this->getOperator()->add(dftu);
@@ -397,15 +397,15 @@ HamiltLCAO<TK, TR>::HamiltLCAO(
             }
             else
             {
-                dftu = new DFTUNew<OperatorLCAO<TK, TR>>(
+                dftu = new DFTU<OperatorLCAO<TK, TR>>(
                     LM_in,
                     this->kv->kvec_d,
                     this->hR,
                     &(this->getHk(LM_in)),
-                    &GlobalC::ucell,
+                    GlobalC::ucell,
                     &GlobalC::GridD,
                     &GlobalC::dftu,
-                    LM_in->ParaV
+                    *(LM_in->ParaV)
                 );
             }
             this->getOperator()->add(dftu);

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/CMakeLists.txt

@@ -12,7 +12,7 @@ add_library(
     td_ekinetic_lcao.cpp
     td_nonlocal_lcao.cpp
     sc_lambda_lcao.cpp
-    dftu_new.cpp
+    dftu_lcao.cpp
 )
 
 if(ENABLE_COVERAGE)

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu.hpp

@@ -0,0 +1,18 @@
+namespace hamilt
+{
+
+#ifndef __DFTUTEMPLATE
+#define __DFTUTEMPLATE
+
+/// The DFTU class template inherits from class T
+/// it is used to calculate the non-local pseudopotential of wavefunction basis
+/// Template parameters:
+/// - T: base class, it would be OperatorLCAO<TK, TR> or OperatorPW<TK>
+template <class T>
+class DFTU : public T
+{
+};
+
+#endif
+
+}

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp

@@ -1,24 +1,34 @@
 #pragma once
-#include "dftu_new.h"
+#include "dftu_lcao.h"
 #include "module_base/parallel_reduce.h"
 
 namespace hamilt
 {
 
 template <typename TK, typename TR>
-void DFTUNew<OperatorLCAO<TK, TR>>::cal_force_stress(
+void DFTU<OperatorLCAO<TK, TR>>::cal_force_stress(
   const bool cal_force, 
   const bool cal_stress, 
   ModuleBase::matrix& force, 
   ModuleBase::matrix& stress)
 {
-    ModuleBase::TITLE("DFTUNew", "cal_force_stress");
-#ifdef __DEBUG
-    assert(this->dm_in_dftu != nullptr);
-#endif    
-    ModuleBase::timer::tick("DFTUNew", "cal_force_stress");
+    ModuleBase::TITLE("DFTU", "cal_force_stress");    
+    if(this->dftu->get_dmr(0) == nullptr)
+    { 
+        ModuleBase::WARNING_QUIT("DFTU", "dmr is not set");
+    }
+    //try to get the density matrix, if the density matrix is empty, skip the calculation and return
+    const hamilt::HContainer<double>* dmR_tmp[this->nspin];
+    dmR_tmp[0] = this->dftu->get_dmr(0);
+    if(this->nspin==2) dmR_tmp[1] = this->dftu->get_dmr(1);
+    if(dmR_tmp[0]->size_atom_pairs() == 0)
+    {
+        return;
+    }
+    // begin the calculation of force and stress
+    ModuleBase::timer::tick("DFTU", "cal_force_stress");
 
-    const Parallel_Orbitals* paraV = this->dm_in_dftu->get_DMR_pointer(1)->get_atom_pair(0).get_paraV();
+    const Parallel_Orbitals* paraV = dmR_tmp[0]->get_atom_pair(0).get_paraV();
     const int npol = this->ucell->get_npol();
     std::vector<double> stress_tmp(6, 0);
     if (cal_force)
@@ -81,7 +91,7 @@ void DFTUNew<OperatorLCAO<TK, TR>>::cal_force_stress(
                 uot.two_center_bundle->overlap_orb_onsite->snap(
                         T1, L1, N1, M1, T0, dtau * this->ucell->lat0, 1 /*cal_deri*/, nlm);
 #else
-                ModuleBase::WARNING_QUIT("DFTUNew", "old two center integral method not implemented");
+                ModuleBase::WARNING_QUIT("DFTU", "old two center integral method not implemented");
 #endif
                 // select the elements of nlm with target_L
                 std::vector<double> nlm_target(tlp1 * 4);
@@ -105,17 +115,13 @@ void DFTUNew<OperatorLCAO<TK, TR>>::cal_force_stress(
             }
         }
         //first iteration to calculate occupation matrix
-        std::vector<double> occ(tlp1 * tlp1 * GlobalV::NSPIN, 0);
+        std::vector<double> occ(tlp1 * tlp1 * this->nspin, 0);
         for(int i=0;i<occ.size();i++)
         {
             const int is = i / (tlp1 * tlp1);
             const int ii = i % (tlp1 * tlp1);
             occ[i] = this->dftu->locale[iat0][target_L][0][is].c[ii];
         }
-        hamilt::HContainer<double>* dmR_tmp[GlobalV::NSPIN];
-        dmR_tmp[0] = this->dm_in_dftu->get_DMR_pointer(1);
-        if(GlobalV::NSPIN==2) dmR_tmp[1] = this->dm_in_dftu->get_DMR_pointer(2);
-
 
         //calculate VU
         const double u_value = this->dftu->U[T0];
@@ -151,20 +157,20 @@ void DFTUNew<OperatorLCAO<TK, TR>>::cal_force_stress(
                 ModuleBase::Vector3<int> R_vector(R_index2[0] - R_index1[0],
                                                   R_index2[1] - R_index1[1],
                                                   R_index2[2] - R_index1[2]);
-                const hamilt::BaseMatrix<double>* tmp[GlobalV::NSPIN];
+                const hamilt::BaseMatrix<double>* tmp[this->nspin];
                 tmp[0] = dmR_tmp[0]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]);
-                if(GlobalV::NSPIN == 2)
+                if(this->nspin == 2)
                 {
                     tmp[1] = dmR_tmp[1]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]);
                 }
                 // if not found , skip this pair of atoms
                 if (tmp[0] != nullptr)
                 {
                     // calculate force
-                    if (cal_force) this->cal_force_IJR(iat1, iat2, T0, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, GlobalV::NSPIN, force_tmp1, force_tmp2);
+                    if (cal_force) this->cal_force_IJR(iat1, iat2, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, this->nspin, force_tmp1, force_tmp2);
 
                     // calculate stress
-                    if (cal_stress) this->cal_stress_IJR(iat1, iat2, T0, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, GlobalV::NSPIN, dis1, dis2, stress_tmp.data());
+                    if (cal_stress) this->cal_stress_IJR(iat1, iat2, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, this->nspin, dis1, dis2, stress_tmp.data());
                 }
             }
         }
@@ -202,14 +208,13 @@ void DFTUNew<OperatorLCAO<TK, TR>>::cal_force_stress(
         stress.c[3] = stress.c[1]; // stress(1,0)
     }
 
-    ModuleBase::timer::tick("DFTUNew", "cal_force_stress");
+    ModuleBase::timer::tick("DFTU", "cal_force_stress");
 }
 
 template <typename TK, typename TR>
-void DFTUNew<OperatorLCAO<TK, TR>>::cal_force_IJR(
+void DFTU<OperatorLCAO<TK, TR>>::cal_force_IJR(
     const int& iat1,
     const int& iat2,
-    const int& T0,
     const Parallel_Orbitals* paraV,
     const std::unordered_map<int, std::vector<double>>& nlm1_all,
     const std::unordered_map<int, std::vector<double>>& nlm2_all,
@@ -272,10 +277,9 @@ void DFTUNew<OperatorLCAO<TK, TR>>::cal_force_IJR(
 }
 
 template <typename TK, typename TR>
-void DFTUNew<OperatorLCAO<TK, TR>>::cal_stress_IJR(
+void DFTU<OperatorLCAO<TK, TR>>::cal_stress_IJR(
     const int& iat1,
     const int& iat2,
-    const int& T0,
     const Parallel_Orbitals* paraV,
     const std::unordered_map<int, std::vector<double>>& nlm1_all,
     const std::unordered_map<int, std::vector<double>>& nlm2_all,