update read in charge density codes and tests

Author: mohanchen

source/module_elecstate/module_charge/charge_init.cpp

@@ -28,8 +28,12 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
 {
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "init_chg", PARAM.inp.init_chg);
 
+    const int nspin = PARAM.inp.nspin;
+    assert(nspin>0);
+
     std::cout << " START CHARGE      : " << PARAM.inp.init_chg << std::endl;
-    //here we need to set the omega for the charge density
+
+    // we need to set the omega for the charge density
     set_omega(&ucell.omega);
     this->pgrid = &pgrid;
 
@@ -46,17 +50,17 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
         if (ModuleIO::read_rhog(binary.str(), rhopw, rhog))
         {
             GlobalV::ofs_running << " Read electron density from file: " << binary.str() << std::endl;
-            for (int is = 0; is < PARAM.inp.nspin; ++is)
+            for (int is = 0; is < nspin; ++is)
             {
                 rhopw->recip2real(rhog[is], rho[is]);
             }
         }
         else
         {
-            for (int is = 0; is < PARAM.inp.nspin; ++is)
+            for (int is = 0; is < nspin; ++is)
             {
                 std::stringstream ssc;
-                ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
+                ssc << PARAM.globalv.global_readin_dir << "chgs" << is + 1 << ".cube";
                 if (ModuleIO::read_vdata_palgrid(pgrid,
                     (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_BPGROUP : GlobalV::MY_RANK),
                     GlobalV::ofs_running,
@@ -103,7 +107,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
         {
             const std::string warn_msg
                 = " WARNING: \"init_chg\" is enabled but ABACUS failed to read\n charge density from file.\n"
-                  " Please check if there is SPINX_CHG.cube (X=1,...) or\n {suffix}-CHARGE-DENSITY.restart in the "
+                  " Please check if there is chgsx.cube (x=1,2,etc.) or\n {suffix}-CHARGE-DENSITY.restart in the "
                   "directory.\n";
             std::cout << warn_msg;
             if (PARAM.inp.init_chg == "file")
@@ -121,9 +125,9 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             {
                 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_space(nspin * this->ngmc, {0.0, 0.0});
                 std::vector<std::complex<double>*> kin_g;
-                for (int is = 0; is < PARAM.inp.nspin; is++)
+                for (int is = 0; is < nspin; is++)
                 {
                     kin_g.push_back(kin_g_space.data() + is * this->ngmc);
                 }
@@ -133,14 +137,14 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                 if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
                 {
                     GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
-                    for (int is = 0; is < PARAM.inp.nspin; ++is)
+                    for (int is = 0; is < nspin; ++is)
                     {
                         rhopw->recip2real(kin_g[is], this->kin_r[is]);
                     }
                 }
                 else
                 {
-                    for (int is = 0; is < PARAM.inp.nspin; is++)
+                    for (int is = 0; is < nspin; is++)
                     {
                         std::stringstream ssc;
                         ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
@@ -175,16 +179,16 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
         }
     }
 
-    if (PARAM.inp.init_chg == "atomic" || read_error) // mohan add 2007-10-17
+    if (PARAM.inp.init_chg == "atomic" || read_error)
     {
         if (read_error)
         {
             std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl;
         }
-        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
+        this->atomic_rho(nspin, ucell.omega, rho, strucFac, ucell);
     }
 
-    // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
+    // initial tau = 3/5 rho^2/3, Thomas-Fermi
     if (XC_Functional::get_ked_flag())
     {
         if (PARAM.inp.init_chg == "atomic" || read_kin_error)
@@ -194,11 +198,11 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                 std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
             }
             const double fact = (3.0 / 5.0) * pow(3.0 * ModuleBase::PI * ModuleBase::PI, 2.0 / 3.0);
-            for (int is = 0; is < PARAM.inp.nspin; ++is)
+            for (int is = 0; is < nspin; ++is)
             {
                 for (int ir = 0; ir < this->rhopw->nrxx; ++ir)
                 {
-                    kin_r[is][ir] = fact * pow(std::abs(rho[is][ir]) * PARAM.inp.nspin, 5.0 / 3.0) / PARAM.inp.nspin;
+                    kin_r[is][ir] = fact * pow(std::abs(rho[is][ir]) * nspin, 5.0 / 3.0) / nspin;
                 }
             }
         }
@@ -207,7 +211,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
     // Peize Lin add 2020.04.04
     if (GlobalC::restart.info_load.load_charge && !GlobalC::restart.info_load.load_charge_finish)
     {
-        for (int is = 0; is < PARAM.inp.nspin; ++is)
+        for (int is = 0; is < nspin; ++is)
         {
             try
             {
@@ -217,20 +221,21 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             {
                 // try to load from the output of `out_chg` 
                 std::stringstream ssc;
-                ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
+                ssc << PARAM.globalv.global_readin_dir << "chgs" << is + 1 << ".cube";
                 if (ModuleIO::read_vdata_palgrid(pgrid,
                     (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_BPGROUP : GlobalV::MY_RANK),
                     GlobalV::ofs_running,
                     ssc.str(),
                     this->rho[is],
                     ucell.nat))
                 {
-                    GlobalV::ofs_running << " Read in the electron density: " << ssc.str() << std::endl;
+                    GlobalV::ofs_running << " Read in electron density: " << ssc.str() << std::endl;
                 }
             }
         }
         GlobalC::restart.info_load.load_charge_finish = true;
     }
+
 #ifdef __MPI
     this->init_chgmpi();
 #endif
@@ -248,7 +253,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                 PARAM.globalv.global_readin_dir,
 				GlobalV::KPAR, GlobalV::MY_POOL, GlobalV::MY_RANK, 
                 GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL,
-				PARAM.inp.nbands, PARAM.inp.nspin, PARAM.globalv.npol,
+				PARAM.inp.nbands, nspin, PARAM.globalv.npol,
 				kv->get_nkstot(),kv->ik2iktot,kv->isk,GlobalV::ofs_running);
     }
 }