Merge branch 'develop' into fix_container

Author: sunliang98

CMakeLists.txt

@@ -286,10 +286,6 @@ if (USE_SW)
   include_directories(${SW_FFT}/include)
 
   target_link_libraries(${ABACUS_BIN_NAME} ${SW_FFT}/lib/libfftw3.a)
-  target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswfft.a)
-  target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswscalapack.a)
-  target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswlapack.a)
-  target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswblas.a)
 endif()
 
 find_package(Threads REQUIRED)
@@ -470,10 +466,7 @@ if(MKLROOT)
   if(CMAKE_CXX_COMPILER_ID MATCHES Intel)
     list(APPEND math_libs ifcore)
   endif()
-elseif(USE_SW)
-   list(APPEND math_libs gfortran)
-  # SW architecture can only use its own math library
-else()
+elseif(NOT USE_SW)
   find_package(FFTW3 REQUIRED)
   find_package(Lapack REQUIRED)
   include_directories(${FFTW3_INCLUDE_DIRS})
@@ -796,6 +789,17 @@ if(ENABLE_RAPIDJSON)
   target_link_libraries(${ABACUS_BIN_NAME} json_output)
 endif()
 
+if (USE_SW)
+  target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswfft.a)
+  target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswscalapack.a)
+  target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswlapack.a)
+  target_link_libraries(${ABACUS_BIN_NAME} ${SW_MATH}/libswblas.a)
+  list(APPEND math_libs gfortran)
+endif()
+
+list(APPEND math_libs m)
+target_link_libraries(${ABACUS_BIN_NAME} ${math_libs})
+
 install(PROGRAMS ${ABACUS_BIN_PATH}
         TYPE BIN
         # DESTINATION ${CMAKE_INSTALL_BINDIR}

source/Makefile.Objects

@@ -643,6 +643,7 @@ OBJS_LCAO=evolve_elec.o\
       LCAO_init_basis.o\
       setup_exx.o\
       setup_deepks.o\
+      setup_dm.o\
       rho_tau_lcao.o\
       center2_orb.o\
       center2_orb-orb11.o\

source/source_base/test_parallel/CMakeLists.txt

@@ -67,7 +67,7 @@ add_test(NAME MODULE_BASE_parallel_2d_test_para
 
   # figure out the lib that provides BLACS
   if(MKLROOT)
-    list(APPEND BLACS_LIB MKL::MKL MKL::MKL_SCALAPACK)
+    list(APPEND BLACS_LIB MKL::MKL)
   else()
     set(BLACS_LIB ScaLAPACK::ScaLAPACK)
   endif()

source/source_esolver/esolver_dm2rho.cpp

@@ -46,30 +46,22 @@ void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
 
     ESolver_KS_LCAO<TK, TR>::before_scf(ucell, istep);
 
-    auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
-
-    if(!estate)
-    {
-        ModuleBase::WARNING_QUIT("ESolver_DM2rho::after_scf","pelec does not exist");
-    }
-
     // file name of DM
     std::string zipname = "output_DM0.npz";
-    elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
 
     // read DM from file
-    ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1)));
+    ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(1)));
 
     // if nspin=2, need extra reading
     if (PARAM.inp.nspin == 2)
     {
         zipname = "output_DM1.npz";
-        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2)));
+        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(2)));
     }
 
     // it's dangerous to design psiToRho function like this, mohan note 20251024
     // this->pelec->psiToRho(*this->psi);
-    LCAO_domain::dm2rho(estate->DM->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
+    LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
 
     int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1;
 

source/source_esolver/esolver_double_xc.cpp

@@ -84,7 +84,7 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
     }
 
     // 6) initialize the density matrix
-    dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin);
+    this->dmat_base.allocate_dm(&this->kv, &this->pv, PARAM.inp.nspin);
 
     // 10) inititlize the charge density
     this->chr_base.allocate(PARAM.inp.nspin);
@@ -138,8 +138,6 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
     if (this->p_hamilt_base == nullptr)
     {
-        elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->get_DM();
-
         this->p_hamilt_base = new hamilt::HamiltLCAO<TK, TR>(
             ucell,
             this->gd,
@@ -148,7 +146,7 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
             this->kv,
             this->two_center_bundle_,
             this->orb_,
-            DM,
+            this->dmat_base.dm,
 			this->deepks,
 			istep,
 			this->exx_nao);
@@ -159,13 +157,11 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); 
 
     // DMR should be same size with Hamiltonian(R)
-    dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)
-    ->get_DM()
-    ->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt_base)->getHR()));
+    this->dmat_base.dm->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt_base)->getHR()));
 
     if (istep > 0)
     {
-        dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->get_DM()->cal_DMR();
+        this->dmat_base.dm->cal_DMR();
     }
 
     ModuleBase::timer::tick("ESolver_DoubleXC", "before_scf");
@@ -226,23 +222,23 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
         std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&this->deepks.ld, [](LCAO_Deepks<TK>*) {});
         LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
 
-        deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
-                                            this->kv.get_nks(),
-                                            ucell.nat,
-                                            PARAM.globalv.nlocal,
-                                            this->pelec->ekb,
-                                            this->kv.kvec_d,
-                                            ucell,
-                                            this->orb_,
-                                            this->gd,
-                                            &(this->pv),
-                                            *(this->psi),
-                                            dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
-                                            p_ham_deepks,
-                                            iter,
-                                            conv_esolver,
-                                            GlobalV::MY_RANK,
-                                            GlobalV::ofs_running);
+		deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
+				this->kv.get_nks(),
+				ucell.nat,
+				PARAM.globalv.nlocal,
+				this->pelec->ekb,
+				this->kv.kvec_d,
+				ucell,
+				this->orb_,
+				this->gd,
+				&(this->pv),
+				*(this->psi),
+				this->dmat.dm,
+				p_ham_deepks,
+				iter,
+				conv_esolver,
+				GlobalV::MY_RANK,
+				GlobalV::ofs_running);
 #endif
 
         // restore to density after charge mixing
@@ -352,9 +348,12 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
             auto _pes_lcao = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
             for (int ik = 0; ik < nks; ik++)
             {
-                _pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik));
+// mohan update 2025-11-03
+                this->dmat_base.dm->set_DMK_pointer(ik, this->dmat.dm->get_DMK_pointer(ik));
+//                _pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik));
             }
-            _pes_lcao_base->get_DM()->cal_DMR();
+            this->dmat_base.dm->cal_DMR();
+//            _pes_lcao_base->get_DM()->cal_DMR();
             _pes_lcao_base->ekb = _pes_lcao->ekb;
             _pes_lcao_base->wg = _pes_lcao->wg;          
         }        
@@ -386,6 +385,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                        this->gd,
                        this->pv,
                        this->pelec_base,
+                       this->dmat_base, // mohan add 2025-11-03
                        this->psi,
                        this->two_center_bundle_,
                        this->orb_,

source/source_esolver/esolver_double_xc.h

@@ -32,6 +32,9 @@ class ESolver_DoubleXC : public ESolver_KS_LCAO<TK, TR>
     //! Electronic states
     elecstate::ElecState* pelec_base = nullptr;
 
+    //! Density Matrix, mohan add 2025-11-03
+    LCAO_domain::Setup_DM<TK> dmat_base;
+
     //! Electorn charge density
     Charge chr_base;
 };

source/source_esolver/esolver_fp.cpp

@@ -54,9 +54,7 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp)
 
     // write geometry file
     ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif",
-                               ucell,
-                               "# Generated by ABACUS ModuleIO::CifParser",
-                               "data_?");
+      ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?");
 
     // init charge extrapolation
     this->CE.Init_CE(inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap);

source/source_esolver/esolver_ks.cpp

@@ -45,7 +45,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 {
     ModuleBase::TITLE("ESolver_KS", "before_all_runners");
 
-    //! 1) init "before_all_runniers" in ESolver_FP
+    //! 1) setup "before_all_runniers" in ESolver_FP
     ESolver_FP::before_all_runners(ucell, inp);
 
     //! 2) setup some parameters
@@ -67,7 +67,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
 
-    //! 4) setup Exc for the first element '0', because all elements have same exc 
+    //! 4) setup Exc for the first element '0' (all elements have same exc) 
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
 
     //! 5) setup the charge mixing parameters
@@ -84,7 +84,7 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
     }
 
-    //! 7) Setup the k points according to symmetry.
+    //! 7) setup k points in the Brillouin zone according to symmetry.
     this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
 
@@ -99,19 +99,12 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 			this->pw_rhod->nplane, this->pw_rhod->nrxx, pw_big->nbz, pw_big->bz);
 
     //! 11) calculate the structure factor
-    this->sf.setup_structure_factor(&ucell, Pgrid, this->pw_rhod);
+    this->sf.setup(&ucell, Pgrid, this->pw_rhod);
 }
 
 template <typename T, typename Device>
 void ESolver_KS<T, Device>::hamilt2rho_single(UnitCell& ucell, const int istep, const int iter, const double ethr)
-{
-    ModuleBase::timer::tick(this->classname, "hamilt2rho_single");
-    // Temporarily, before HSolver is constructed, it should be overrided by
-    // LCAO, PW, SDFT and TDDFT.
-    // After HSolver is constructed, LCAO, PW, SDFT should delete their own
-    // hamilt2rho_single() and use:
-    ModuleBase::timer::tick(this->classname, "hamilt2rho_single");
-}
+{}
 
 template <typename T, typename Device>
 void ESolver_KS<T, Device>::hamilt2rho(UnitCell& ucell, const int istep, const int iter, const double ethr)