Refactor: Remove cal_nelec and cal_nbands from unitcell.cpp.

source/Makefile.Objects

@@ -229,6 +229,7 @@ OBJS_ELECSTAT=elecstate.o\
     H_TDDFT_pw.o\
     pot_xc.o\
     cal_ux.o\
+    cal_nelec_nband.o\
     read_pseudo.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\

source/module_cell/cal_atoms_info.h

@@ -1,7 +1,7 @@
 #ifndef CAL_ATOMS_INFO_H
 #define CAL_ATOMS_INFO_H
 #include "module_parameter/parameter.h"
-#include "unitcell.h"
+#include "module_elecstate/cal_nelec_nband.h"
 class CalAtomsInfo
 {
   public:
@@ -58,7 +58,7 @@ class CalAtomsInfo
         }
 
         // calculate the total number of electrons
-        cal_nelec(atoms, ntype, para.input.nelec);
+        elecstate::cal_nelec(atoms, ntype, para.input.nelec);
 
         // autoset and check GlobalV::NBANDS
         std::vector<double> nelec_spin(2, 0.0);
@@ -67,7 +67,7 @@ class CalAtomsInfo
             nelec_spin[0] = (para.inp.nelec + para.inp.nupdown ) / 2.0;
             nelec_spin[1] = (para.inp.nelec - para.inp.nupdown ) / 2.0;
         }
-        cal_nbands(para.inp.nelec, para.sys.nlocal, nelec_spin, para.input.nbands);
+        elecstate::cal_nbands(para.inp.nelec, para.sys.nlocal, nelec_spin, para.input.nbands);
         return;
     }
 };

source/module_cell/test/CMakeLists.txt

@@ -26,6 +26,7 @@ list(APPEND cell_simple_srcs
     ../read_pp_blps.cpp
     ../check_atomic_stru.cpp
     ../../module_elecstate/read_pseudo.cpp
+    ../../module_elecstate/cal_nelec_nband.cpp
 )
 
 add_library(cell_info OBJECT ${cell_simple_srcs})

source/module_cell/test/unitcell_test_readpp.cpp

@@ -90,7 +90,7 @@ Magnetism::~Magnetism() { delete[] this->start_magnetization; }
  * possible of an element
  *   - CalNelec: UnitCell::cal_nelec
  *     - calculate the total number of valence electrons from psp files
- *   - CalNbands: elecstate::ElecState::cal_nbands()
+ *   - CalNbands: elecstate::cal_nbands()
  *     - calculate the number of bands
  */
 
@@ -406,22 +406,22 @@ TEST_F(UcellTest, CalNelec) {
     EXPECT_EQ(1, ucell->atoms[0].na);
     EXPECT_EQ(2, ucell->atoms[1].na);
     double nelec = 0;
-    cal_nelec(ucell->atoms, ucell->ntype, nelec);
+    elecstate::cal_nelec(ucell->atoms, ucell->ntype, nelec);
     EXPECT_DOUBLE_EQ(6, nelec);
 }
 
 TEST_F(UcellTest, CalNbands)
 {
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 6);
 }
 
 TEST_F(UcellTest, CalNbandsFractionElec)
 {
     PARAM.input.nelec = 9.5;
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 6);
 }
 
@@ -430,22 +430,22 @@ TEST_F(UcellTest, CalNbandsSOC)
     PARAM.input.lspinorb = true;
     PARAM.input.nbands = 0;
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 20);
 }
 
 TEST_F(UcellTest, CalNbandsSDFT)
 {
     PARAM.input.esolver_type = "sdft";
     std::vector<double> nelec_spin(2, 5.0);
-    EXPECT_NO_THROW(cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands));
+    EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands));
 }
 
 TEST_F(UcellTest, CalNbandsLCAO)
 {
     PARAM.input.basis_type = "lcao";
     std::vector<double> nelec_spin(2, 5.0);
-    EXPECT_NO_THROW(cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands));
+    EXPECT_NO_THROW(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands));
 }
 
 TEST_F(UcellTest, CalNbandsLCAOINPW)
@@ -454,7 +454,7 @@ TEST_F(UcellTest, CalNbandsLCAOINPW)
     PARAM.sys.nlocal = PARAM.input.nbands - 1;
     std::vector<double> nelec_spin(2, 5.0);
     testing::internal::CaptureStdout();
-    EXPECT_EXIT(cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
+    EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
     output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("NLOCAL < NBANDS"));
 }
@@ -464,7 +464,7 @@ TEST_F(UcellTest, CalNbandsWarning1)
     PARAM.input.nbands = PARAM.input.nelec / 2 - 1;
     std::vector<double> nelec_spin(2, 5.0);
     testing::internal::CaptureStdout();
-    EXPECT_EXIT(cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
+    EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
     output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("Too few bands!"));
 }
@@ -477,7 +477,7 @@ TEST_F(UcellTest, CalNbandsWarning2)
     nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0;
     nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0;
     testing::internal::CaptureStdout();
-    EXPECT_EXIT(cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
+    EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
     output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("Too few spin up bands!"));
 }
@@ -490,7 +490,7 @@ TEST_F(UcellTest, CalNbandsWarning3)
     nelec_spin[0] = (PARAM.input.nelec + PARAM.input.nupdown ) / 2.0;
     nelec_spin[1] = (PARAM.input.nelec - PARAM.input.nupdown ) / 2.0;
     testing::internal::CaptureStdout();
-    EXPECT_EXIT(cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
+    EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
     output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("Too few spin down bands!"));
 }
@@ -500,7 +500,7 @@ TEST_F(UcellTest, CalNbandsSpin1)
     PARAM.input.nspin = 1;
     PARAM.input.nbands = 0;
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 15);
 }
 
@@ -510,7 +510,7 @@ TEST_F(UcellTest, CalNbandsSpin1LCAO)
     PARAM.input.nbands = 0;
     PARAM.input.basis_type = "lcao";
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 6);
 }
 
@@ -519,7 +519,7 @@ TEST_F(UcellTest, CalNbandsSpin4)
     PARAM.input.nspin = 4;
     PARAM.input.nbands = 0;
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 30);
 }
 
@@ -529,7 +529,7 @@ TEST_F(UcellTest, CalNbandsSpin4LCAO)
     PARAM.input.nbands = 0;
     PARAM.input.basis_type = "lcao";
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 6);
 }
 
@@ -538,7 +538,7 @@ TEST_F(UcellTest, CalNbandsSpin2)
     PARAM.input.nspin = 2;
     PARAM.input.nbands = 0;
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 16);
 }
 
@@ -548,7 +548,7 @@ TEST_F(UcellTest, CalNbandsSpin2LCAO)
     PARAM.input.nbands = 0;
     PARAM.input.basis_type = "lcao";
     std::vector<double> nelec_spin(2, 5.0);
-    cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
+    elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands);
     EXPECT_EQ(PARAM.input.nbands, 6);
 }
 
@@ -558,7 +558,7 @@ TEST_F(UcellTest, CalNbandsGaussWarning)
     std::vector<double> nelec_spin(2, 5.0);
     PARAM.input.smearing_method = "gaussian";
     testing::internal::CaptureStdout();
-    EXPECT_EXIT(cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
+    EXPECT_EXIT(elecstate::cal_nbands(PARAM.input.nelec, PARAM.sys.nlocal, nelec_spin, PARAM.input.nbands), ::testing::ExitedWithCode(1), "");
     output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("for smearing, num. of bands > num. of occupied bands"));
 }

source/module_cell/test_pw/CMakeLists.txt

@@ -13,7 +13,7 @@ AddTest(
   LIBS parameter ${math_libs} base device 
   SOURCES unitcell_test_pw.cpp ../unitcell.cpp ../read_atoms.cpp ../atom_spec.cpp
 	../atom_pseudo.cpp ../pseudo.cpp ../read_pp.cpp ../read_pp_complete.cpp ../read_pp_upf201.cpp ../read_pp_upf100.cpp
-	../read_pp_vwr.cpp ../read_pp_blps.cpp ../../module_io/output.cpp ../../module_elecstate/read_pseudo.cpp
+	../read_pp_vwr.cpp ../read_pp_blps.cpp ../../module_io/output.cpp ../../module_elecstate/read_pseudo.cpp ../../module_elecstate/cal_nelec_nband.cpp
 )
 
 find_program(BASH bash)

source/module_cell/unitcell.cpp

@@ -1328,160 +1328,6 @@ void UnitCell::remake_cell() {
     }
 }
 
-void cal_nelec(const Atom* atoms, const int& ntype, double& nelec)
-{
-    ModuleBase::TITLE("UnitCell", "cal_nelec");
-    GlobalV::ofs_running << "\n SETUP THE ELECTRONS NUMBER" << std::endl;
-
-    if (nelec == 0)
-    {
-        if (PARAM.inp.use_paw)
-        {
-#ifdef USE_PAW
-            for (int it = 0; it < ntype; it++)
-            {
-                std::stringstream ss1, ss2;
-                ss1 << " electron number of element " << GlobalC::paw_cell.get_zat(it) << std::endl;
-                const int nelec_it = GlobalC::paw_cell.get_val(it) * atoms[it].na;
-                nelec += nelec_it;
-                ss2 << "total electron number of element " << GlobalC::paw_cell.get_zat(it);
-
-                ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss1.str(), GlobalC::paw_cell.get_val(it));
-                ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss2.str(), nelec_it);
-            }
-#endif
-        }
-        else
-        {
-            for (int it = 0; it < ntype; it++)
-            {
-                std::stringstream ss1, ss2;
-                ss1 << "electron number of element " << atoms[it].label;
-                const double nelec_it = atoms[it].ncpp.zv * atoms[it].na;
-                nelec += nelec_it;
-                ss2 << "total electron number of element " << atoms[it].label;
-
-                ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss1.str(), atoms[it].ncpp.zv);
-                ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss2.str(), nelec_it);
-            }
-            ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "AUTOSET number of electrons: ", nelec);
-        }
-    }
-    if (PARAM.inp.nelec_delta != 0)
-    {
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,
-                                    "nelec_delta is NOT zero, please make sure you know what you are "
-                                    "doing! nelec_delta: ",
-                                    PARAM.inp.nelec_delta);
-        nelec += PARAM.inp.nelec_delta;
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nelec now: ", nelec);
-    }
-    return;
-}
-
-void cal_nbands(const int& nelec, const int& nlocal, const std::vector<double>& nelec_spin, int& nbands)
-{
-    if (PARAM.inp.esolver_type == "sdft") // qianrui 2021-2-20
-    {
-        return;
-    }
-    //=======================================
-    // calculate number of bands (setup.f90)
-    //=======================================
-    double occupied_bands = static_cast<double>(nelec / ModuleBase::DEGSPIN);
-    if (PARAM.inp.lspinorb == 1) {
-        occupied_bands = static_cast<double>(nelec);
-    }
-
-    if ((occupied_bands - std::floor(occupied_bands)) > 0.0)
-    {
-        occupied_bands = std::floor(occupied_bands) + 1.0; // mohan fix 2012-04-16
-    }
-
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "occupied bands", occupied_bands);
-
-    if (nbands == 0)
-    {
-        if (PARAM.inp.nspin == 1)
-        {
-            const int nbands1 = static_cast<int>(occupied_bands) + 10;
-            const int nbands2 = static_cast<int>(1.2 * occupied_bands) + 1;
-            nbands = std::max(nbands1, nbands2);
-            if (PARAM.inp.basis_type != "pw") {
-                nbands = std::min(nbands, nlocal);
-            }
-        }
-        else if (PARAM.inp.nspin == 4)
-        {
-            const int nbands3 = nelec + 20;
-            const int nbands4 = static_cast<int>(1.2 * nelec) + 1;
-            nbands = std::max(nbands3, nbands4);
-            if (PARAM.inp.basis_type != "pw") {
-                nbands = std::min(nbands, nlocal);
-            }
-        }
-        else if (PARAM.inp.nspin == 2)
-        {
-            const double max_occ = std::max(nelec_spin[0], nelec_spin[1]);
-            const int nbands3 = static_cast<int>(max_occ) + 11;
-            const int nbands4 = static_cast<int>(1.2 * max_occ) + 1;
-            nbands = std::max(nbands3, nbands4);
-            if (PARAM.inp.basis_type != "pw") {
-                nbands = std::min(nbands, nlocal);
-            }
-        }
-        ModuleBase::GlobalFunc::AUTO_SET("NBANDS", nbands);
-    }
-    // else if ( PARAM.inp.calculation=="scf" || PARAM.inp.calculation=="md" || PARAM.inp.calculation=="relax") //pengfei
-    // 2014-10-13
-    else
-    {
-        if (nbands < occupied_bands) {
-            ModuleBase::WARNING_QUIT("unitcell", "Too few bands!");
-        }
-        if (PARAM.inp.nspin == 2)
-        {
-            if (nbands < nelec_spin[0])
-            {
-                ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nelec_up", nelec_spin[0]);
-                ModuleBase::WARNING_QUIT("ElecState::cal_nbands", "Too few spin up bands!");
-            }
-            if (nbands < nelec_spin[1])
-            {
-                ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nelec_down", nelec_spin[1]);
-                ModuleBase::WARNING_QUIT("ElecState::cal_nbands", "Too few spin down bands!");
-            }
-        }
-    }
-
-    // mohan add 2010-09-04
-    // std::cout << "nbands(this-> = " <<nbands <<std::endl;
-    if (nbands == occupied_bands)
-    {
-        if (PARAM.inp.smearing_method != "fixed")
-        {
-            ModuleBase::WARNING_QUIT("ElecState::cal_nbands", "for smearing, num. of bands > num. of occupied bands");
-        }
-    }
-
-    // mohan update 2021-02-19
-    // mohan add 2011-01-5
-    if (PARAM.inp.basis_type == "lcao" || PARAM.inp.basis_type == "lcao_in_pw")
-    {
-        if (nbands > nlocal)
-        {
-            ModuleBase::WARNING_QUIT("ElecState::cal_nbandsc", "NLOCAL < NBANDS");
-        }
-        else
-        {
-            ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "NLOCAL", nlocal);
-            ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "NBANDS", nbands);
-        }
-    }
-
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "NBANDS", nbands);
-}
-
 void UnitCell::compare_atom_labels(std::string label1, std::string label2) {
     if (label1
         != label2) //'!( "Ag" == "Ag" || "47" == "47" || "Silver" == Silver" )'

source/module_cell/unitcell.h

@@ -329,23 +329,4 @@ class UnitCell {
     std::vector<ModuleBase::Vector3<int>> get_constrain() const;
 };
 
-/**
- * @brief calculate the total number of electrons in system
- *
- * @param atoms [in] atom pointer
- * @param ntype [in] number of atom types
- * @param nelec [out] total number of electrons
- */
-void cal_nelec(const Atom* atoms, const int& ntype, double& nelec);
-
-/**
- * @brief Calculate the number of bands.
- *
- * @param nelec [in] total number of electrons
- * @param nlocal [in] total number of local basis
- * @param nelec_spin [in] number of electrons for each spin
- * @param nbands  [out] number of bands
- */
-void cal_nbands(const int& nelec, const int& nlocal, const std::vector<double>& nelec_spin, int& nbands);
-
 #endif // unitcell class