Merge branch 'develop' into ucell20

Author: A-006

source/module_basis/module_ao/ORB_read.cpp

@@ -82,7 +82,7 @@ void LCAO_Orbitals::init(
             this->orbital_file.push_back(orbital_dir + orbital_file[it]);
         }
     }
-
+    this->descriptor_file = descriptor_file;
 #ifdef __MPI
     bcast_files(ntype, my_rank);
 #endif
@@ -263,7 +263,6 @@ void LCAO_Orbitals::Read_Orbitals(std::ofstream& ofs_in,
 
         delete[] this->Alpha;
         this->Alpha = new Numerical_Orbital[1]; // not related to atom type -- remain to be discussed
-
         this->Read_Descriptor(ofs_in, force_flag, my_rank);
     }
 

source/module_basis/module_nao/two_center_bundle.cpp

@@ -107,7 +107,7 @@ void TwoCenterBundle::tabulate()
     {
         overlap_orb_alpha = std::unique_ptr<TwoCenterIntegrator>(new TwoCenterIntegrator);
         overlap_orb_alpha->tabulate(*orb_, *alpha_, 'S', nr, cutoff);
-        ModuleBase::Memory::record("TwoCenterTable: Descriptor", overlap_orb_beta->table_memory());
+        ModuleBase::Memory::record("TwoCenterTable: Descriptor", overlap_orb_alpha->table_memory());
     }
 
     if (orb_onsite_)

source/module_cell/unitcell.cpp

@@ -592,6 +592,7 @@ void UnitCell::cal_nwfc(std::ofstream& log) {
     //========================
     this->lmax = 0;
     this->nmax = 0;
+    this->nmax_total = 0;
     for (int it = 0; it < ntype; it++) {
         lmax = std::max(lmax, atoms[it].nwl);
         for (int l = 0; l < atoms[it].nwl + 1; l++) {

source/module_elecstate/module_charge/charge_init.cpp

@@ -37,6 +37,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
     this->pgrid = &pgrid;
 
     bool read_error = false;
+    bool read_kin_error = false;
     if (PARAM.inp.init_chg == "file" || PARAM.inp.init_chg == "auto")
     {
         GlobalV::ofs_running << " try to read charge from file" << std::endl;
@@ -101,72 +102,100 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             }
         }
 
-        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        if (read_error)
         {
-            GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
-            // try to read charge from binary file first, which is the same as QE
-            std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
-            std::vector<std::complex<double>*> kin_g;
-            for (int is = 0; is < PARAM.inp.nspin; is++)
+            const std::string warn_msg
+                = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
+                  " Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the "
+                  "directory.\n";
+            std::cout << std::endl << warn_msg;
+            if (PARAM.inp.init_chg == "file")
             {
-                kin_g.push_back(kin_g_space.data() + is * this->ngmc);
+                ModuleBase::WARNING_QUIT("Charge::init_rho",
+                                         "Failed to read in charge density from file.\nIf you want to use atomic "
+                                         "charge initialization, \nplease set init_chg to atomic in INPUT.");
             }
+        }
 
-            std::stringstream binary;
-            binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
-            if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
+        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        {
+            // If the charge density is not read in, then the kinetic energy density is not read in either
+            if (!read_error)
             {
-                GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
-                for (int is = 0; is < PARAM.inp.nspin; ++is)
+                GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
+                // try to read charge from binary file first, which is the same as QE
+                std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
+                std::vector<std::complex<double>*> kin_g;
+                for (int is = 0; is < PARAM.inp.nspin; is++)
                 {
-                    rhopw->recip2real(kin_g[is], this->kin_r[is]);
+                    kin_g.push_back(kin_g_space.data() + is * this->ngmc);
                 }
-            }
-            else
-            {
-                for (int is = 0; is < PARAM.inp.nspin; is++)
+
+                std::stringstream binary;
+                binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
+                if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
                 {
-                    std::stringstream ssc;
-                    ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
-                    // mohan update 2012-02-10, sunliang update 2023-03-09
-                    if (ModuleIO::read_vdata_palgrid(
-                            pgrid,
-                            (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
-                            GlobalV::ofs_running,
-                            ssc.str(),
-                            this->kin_r[is],
-                            ucell.nat))
+                    GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
+                    for (int is = 0; is < PARAM.inp.nspin; ++is)
                     {
-                        GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
+                        rhopw->recip2real(kin_g[is], this->kin_r[is]);
                     }
                 }
+                else
+                {
+                    for (int is = 0; is < PARAM.inp.nspin; is++)
+                    {
+                        std::stringstream ssc;
+                        ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
+                        // mohan update 2012-02-10, sunliang update 2023-03-09
+                        if (ModuleIO::read_vdata_palgrid(
+                                pgrid,
+                                (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
+                                GlobalV::ofs_running,
+                                ssc.str(),
+                                this->kin_r[is],
+                                ucell.nat))
+                        {
+                            GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
+                        }
+                        else
+                        {
+                            read_kin_error = true;
+                            std::cout << " WARNING: \"init_chg\" is enabled but ABACUS failed to read kinetic energy "
+                                         "density from file.\n"
+                                         " Please check if there is SPINX_TAU.cube (X=1,...) or "
+                                         "{suffix}-TAU-DENSITY.restart in the directory.\n"
+                                      << std::endl;
+                            break;
+                        }
+                    }
+                }
+            }
+            else
+            {
+                read_kin_error = true;
             }
         }
     }
-    if (read_error)
+
+    if (PARAM.inp.init_chg == "atomic" || read_error) // mohan add 2007-10-17
     {
-        const std::string warn_msg = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
-                                     " Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the directory.\n";
-        std::cout << std::endl << warn_msg;
-        if (PARAM.inp.init_chg == "auto")
+        if (read_error)
         {
-            std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl << std::endl;
-        }
-        else if (PARAM.inp.init_chg == "file")
-        {
-            ModuleBase::WARNING_QUIT("Charge::init_rho", "Failed to read in charge density from file.\nIf you want to use atomic charge initialization, \nplease set init_chg to atomic in INPUT.");
+            std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl;
         }
+        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
     }
 
-    if (PARAM.inp.init_chg == "atomic" || 
-        (PARAM.inp.init_chg == "auto" && read_error)) // mohan add 2007-10-17
+    // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
+    if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
-        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
-
-        // liuyu 2023-06-29 : move here from atomic_rho(), which will be called several times in charge extrapolation
-        // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
-        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        if (PARAM.inp.init_chg == "atomic" || read_kin_error)
         {
+            if (read_kin_error)
+            {
+                std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
+            }
             const double pi = 3.141592653589790;
             const double fact = (3.0 / 5.0) * pow(3.0 * pi * pi, 2.0 / 3.0);
             for (int is = 0; is < PARAM.inp.nspin; ++is)

source/module_esolver/esolver_gets.cpp

@@ -147,7 +147,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
 }
 
 void ESolver_GetS::after_all_runners(UnitCell& ucell) {};
-double ESolver_GetS::cal_energy() {};
+double ESolver_GetS::cal_energy() { return 0.0; };
 void ESolver_GetS::cal_force(UnitCell& ucell, ModuleBase::matrix& force) {};
 void ESolver_GetS::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) {};
 

source/module_esolver/esolver_ks.h

@@ -43,7 +43,7 @@ class ESolver_KS : public ESolver_FP
     virtual void iter_init(UnitCell& ucell, const int istep, const int iter);
 
     //! Something to do after hamilt2density function in each iter loop.
-    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter);
+    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override;
 
     // calculate electron density from a specific Hamiltonian with ethr
     virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr);

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp

@@ -55,10 +55,10 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(const Grid_Driv
     auto* paraV = this->hR->get_paraV(); // get parallel orbitals from HR
     // TODO: if paraV is nullptr, AtomPair can not use paraV for constructor, I will repair it in the future.
 
-    // this->H_V_delta = new HContainer<TR>(paraV);
+    this->H_V_delta = new HContainer<TR>(paraV);
     if (std::is_same<TK, double>::value)
     {
-        this->H_V_delta = new HContainer<TR>(paraV);
+        //this->H_V_delta = new HContainer<TR>(paraV);
         this->H_V_delta->fix_gamma();
     }
 
@@ -123,10 +123,10 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(const Grid_Driv
                                          R_index2.y - R_index1.y,
                                          R_index2.z - R_index1.z,
                                          paraV);
-                if (std::is_same<TK, double>::value)
-                {
-                    this->H_V_delta->insert_pair(tmp);
-                }
+                // if (std::is_same<TK, double>::value)
+                // {
+                this->H_V_delta->insert_pair(tmp);
+                // }
             }
         }
         if (pre_cal_nlm)
@@ -135,15 +135,14 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::initialize_HR(const Grid_Driv
         }
     }
     // allocate the memory of BaseMatrix in HR, and set the new values to zero
-    if (std::is_same<TK, double>::value)
-    {
-        // only gamma-only has full size of Hamiltonian of DeePKS now,
-        // multi-k keep same size of nonlocal operator, H_V_delta will be allocated by hR
-        this->H_V_delta->allocate(nullptr, true);
-        // expand hR with H_V_delta, only gamma-only case now
-        this->hR->add(*this->H_V_delta);
-        this->hR->allocate(nullptr, false);
-    }
+    // if (std::is_same<TK, double>::value)
+    // {
+    this->H_V_delta->allocate(nullptr, true);
+        // expand hR with H_V_delta
+        // update : for computational rigor, gamma-only and multi-k cases both have full size of Hamiltonian of DeePKS now
+    this->hR->add(*this->H_V_delta);
+    this->hR->allocate(nullptr, false);
+    // }
 
     ModuleBase::timer::tick("DeePKS", "initialize_HR");
 }
@@ -165,11 +164,11 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
         GlobalC::ld.cal_descriptor(this->ucell->nat);
         GlobalC::ld.cal_gedm(this->ucell->nat);
 
-        // recalculate the H_V_delta
-        if (this->H_V_delta == nullptr)
-        {
-            this->H_V_delta = new hamilt::HContainer<TR>(*this->hR);
-        }
+        // // recalculate the H_V_delta
+        // if (this->H_V_delta == nullptr)
+        // {
+        //     this->H_V_delta = new hamilt::HContainer<std::complex<double>>(*this->hR);
+        // }
         this->H_V_delta->set_zero();
         this->calculate_HR();
 

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h

@@ -61,6 +61,7 @@ class DeePKS<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
     // LCAO_Deepks* ld = nullptr;
 
     const UnitCell* ucell = nullptr;
+    Grid_Driver* gridD = nullptr;
 
     const Grid_Driver* gd = nullptr;
 

source/module_hamilt_lcao/module_deepks/CMakeLists.txt

@@ -30,8 +30,8 @@ if(ENABLE_DEEPKS)
   endif()
 
 # I will rewrite the test later, the current test rely on too many modules
-#  if(BUILD_TESTING)
-#      add_subdirectory(test)
-#  endif()
+  if(BUILD_TESTING)
+      add_subdirectory(test)
+  endif()
 endif()
 

source/module_hamilt_lcao/module_deepks/LCAO_deepks_pdm.cpp

@@ -379,7 +379,7 @@ void LCAO_Deepks::cal_projected_DM(const elecstate::DensityMatrix<TK, double>* d
 
 void LCAO_Deepks::check_projected_dm()
 {
-    const std::string file_projdm = PARAM.globalv.global_out_dir + "deepks_projdm.dat";
+    const std::string file_projdm = PARAM.globalv.global_out_dir + "pdm.dat";
     std::ofstream ofs(file_projdm.c_str());
 
     ofs << std::setprecision(10);