diff --git a/source/Makefile.Objects b/source/Makefile.Objects index c2ddaabdfa..a04e737697 100644 --- a/source/Makefile.Objects +++ b/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\ diff --git a/source/module_cell/cal_atoms_info.h b/source/module_cell/cal_atoms_info.h index c4aa70fb3a..3b59d57540 100644 --- a/source/module_cell/cal_atoms_info.h +++ b/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 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; } }; diff --git a/source/module_cell/test/CMakeLists.txt b/source/module_cell/test/CMakeLists.txt index acb79c77db..c79adcf89d 100644 --- a/source/module_cell/test/CMakeLists.txt +++ b/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}) diff --git a/source/module_cell/test/unitcell_test_readpp.cpp b/source/module_cell/test/unitcell_test_readpp.cpp index 8a693c4c54..714f77d972 100644 --- a/source/module_cell/test/unitcell_test_readpp.cpp +++ b/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,14 +406,14 @@ 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 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); } @@ -421,7 +421,7 @@ TEST_F(UcellTest, CalNbandsFractionElec) { PARAM.input.nelec = 9.5; std::vector 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,7 +430,7 @@ TEST_F(UcellTest, CalNbandsSOC) PARAM.input.lspinorb = true; PARAM.input.nbands = 0; std::vector 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); } @@ -438,14 +438,14 @@ TEST_F(UcellTest, CalNbandsSDFT) { PARAM.input.esolver_type = "sdft"; std::vector 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 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 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 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 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 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 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 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 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 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 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")); } diff --git a/source/module_cell/test_pw/CMakeLists.txt b/source/module_cell/test_pw/CMakeLists.txt index f4bcc63c3f..d171e02aec 100644 --- a/source/module_cell/test_pw/CMakeLists.txt +++ b/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) diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index f0db84ed91..cc7140eb45 100755 --- a/source/module_cell/unitcell.cpp +++ b/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& 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(nelec / ModuleBase::DEGSPIN); - if (PARAM.inp.lspinorb == 1) { - occupied_bands = static_cast(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(occupied_bands) + 10; - const int nbands2 = static_cast(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(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(max_occ) + 11; - const int nbands4 = static_cast(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-> = " < 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" )' diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h index 4b8a04be02..1933a95c2f 100644 --- a/source/module_cell/unitcell.h +++ b/source/module_cell/unitcell.h @@ -329,23 +329,4 @@ class UnitCell { std::vector> 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& nelec_spin, int& nbands); - #endif // unitcell class diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt index 31fd543709..2d0148a0d1 100644 --- a/source/module_elecstate/CMakeLists.txt +++ b/source/module_elecstate/CMakeLists.txt @@ -32,6 +32,7 @@ list(APPEND objects magnetism.cpp occupy.cpp cal_ux.cpp + cal_nelec_nband.cpp read_pseudo.cpp ) diff --git a/source/module_elecstate/cal_nelec_nband.cpp b/source/module_elecstate/cal_nelec_nband.cpp new file mode 100644 index 0000000000..f0d6b79bbc --- /dev/null +++ b/source/module_elecstate/cal_nelec_nband.cpp @@ -0,0 +1,164 @@ +#include "cal_nelec_nband.h" +#include "module_base/constants.h" +#include "module_parameter/parameter.h" +#ifdef USE_PAW +#include "module_cell/module_paw/paw_cell.h" +#endif + +namespace elecstate { + +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& 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(nelec / ModuleBase::DEGSPIN); + if (PARAM.inp.lspinorb == 1) { + occupied_bands = static_cast(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(occupied_bands) + 10; + const int nbands2 = static_cast(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(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(max_occ) + 11; + const int nbands4 = static_cast(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-> = " < 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); +} + +} \ No newline at end of file diff --git a/source/module_elecstate/cal_nelec_nband.h b/source/module_elecstate/cal_nelec_nband.h new file mode 100644 index 0000000000..f965db40c0 --- /dev/null +++ b/source/module_elecstate/cal_nelec_nband.h @@ -0,0 +1,29 @@ +#ifndef CAL_NELEC_NBAND_H +#define CAL_NELEC_NBAND_H + +#include "module_cell/atom_spec.h" + +namespace elecstate { + + /** + * @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& nelec_spin, int& nbands); + +} + +#endif \ No newline at end of file diff --git a/source/module_elecstate/read_pseudo.cpp b/source/module_elecstate/read_pseudo.cpp index 0274fb9ab7..ee11fd2ee7 100644 --- a/source/module_elecstate/read_pseudo.cpp +++ b/source/module_elecstate/read_pseudo.cpp @@ -2,7 +2,6 @@ #include "module_parameter/parameter.h" #include "module_base/global_file.h" #include "module_cell/read_pp.h" -#include "module_cell/cal_atoms_info.h" #include "module_base/element_elec_config.h" #include "module_base/parallel_common.h" diff --git a/source/module_elecstate/read_pseudo.h b/source/module_elecstate/read_pseudo.h index 066abbd239..c7fcf1cb7e 100644 --- a/source/module_elecstate/read_pseudo.h +++ b/source/module_elecstate/read_pseudo.h @@ -2,6 +2,7 @@ #define READ_PSEUDO_H #include "module_cell/unitcell.h" +#include "module_cell/cal_atoms_info.h" namespace elecstate { diff --git a/source/module_hamilt_lcao/module_deepks/test/Makefile.Objects b/source/module_hamilt_lcao/module_deepks/test/Makefile.Objects index 9e62c5551a..5bb621491c 100644 --- a/source/module_hamilt_lcao/module_deepks/test/Makefile.Objects +++ b/source/module_hamilt_lcao/module_deepks/test/Makefile.Objects @@ -74,4 +74,5 @@ sltk_atom_input.o\ sltk_grid.o\ sltk_grid_driver.o -OBJS_ELECSTATE=read_pseudo.o\ \ No newline at end of file +OBJS_ELECSTATE=read_pseudo.o\ +cal_nelec_nband.o\ \ No newline at end of file diff --git a/source/module_md/test/CMakeLists.txt b/source/module_md/test/CMakeLists.txt index 16051470fd..5e5e4a2c4a 100644 --- a/source/module_md/test/CMakeLists.txt +++ b/source/module_md/test/CMakeLists.txt @@ -47,6 +47,7 @@ list(APPEND depend_files ../../module_base/parallel_global.cpp ../../module_base/parallel_comm.cpp ../../module_elecstate/read_pseudo.cpp + ../../module_elecstate/cal_nelec_nband.cpp ) AddTest(