diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 7104d8a62e6..d3da386f1c7 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -211,6 +211,7 @@ OBJS_CELL=atom_pseudo.o\ print_cell.o\ setup_nonlocal.o\ klist.o\ + klist_io.o\ reciprocal_grid.o\ cell_index.o\ cell_tools.o\ diff --git a/source/source_cell/CMakeLists.txt b/source/source_cell/CMakeLists.txt index 3c74e1cd4fa..e568aa61c3b 100644 --- a/source/source_cell/CMakeLists.txt +++ b/source/source_cell/CMakeLists.txt @@ -22,6 +22,7 @@ add_library( read_atoms_helper.cpp read_orb.cpp klist.cpp + klist_io.cpp reciprocal_grid.cpp parallel_kpoints.cpp cell_index.cpp diff --git a/source/source_cell/klist.cpp b/source/source_cell/klist.cpp index 1cbfd4d56bd..5869a8ce3f3 100644 --- a/source/source_cell/klist.cpp +++ b/source/source_cell/klist.cpp @@ -1,44 +1,12 @@ -/** - * @file klist.cpp - * @brief Implementation of K_Vectors class. - */ #include "klist.h" +#include "klist_io.h" #include "source_base/formatter.h" #include "source_base/parallel_common.h" #include "source_base/parallel_global.h" #include "source_base/parallel_reduce.h" #include "source_cell/module_symmetry/symmetry.h" -void K_Vectors::cal_ik_global() -{ - const int my_pool = this->para_k.my_pool; - this->ik2iktot.resize(this->nks); -#ifdef __MPI - if(this->nspin == 2) - { - for (int ik = 0; ik < this->nks / 2; ++ik) - { - this->ik2iktot[ik] = this->para_k.startk_pool[my_pool] + ik; - this->ik2iktot[ik + this->nks / 2] = this->nkstot / 2 + this->para_k.startk_pool[my_pool] + ik; - } - } - else - { - for (int ik = 0; ik < this->nks; ++ik) - { - this->ik2iktot[ik] = this->para_k.startk_pool[my_pool] + ik; - } - } -#else - for (int ik = 0; ik < this->nks; ++ik) - { - this->ik2iktot[ik] = ik; - } -#endif - -} - void K_Vectors::set(const UnitCell& ucell, const ModuleSymmetry::Symmetry& symm, const std::string& k_file_name, @@ -46,6 +14,7 @@ void K_Vectors::set(const UnitCell& ucell, const ModuleBase::Matrix3& reciprocal_vec, const ModuleBase::Matrix3& latvec, std::ofstream& ofs, + std::ofstream& ofs_warning, const bool use_ibz, const std::string& global_out_dir, const bool gamma_only_local, @@ -71,19 +40,30 @@ void K_Vectors::set(const UnitCell& ucell, const std::string global_out_dir_ = global_out_dir; const bool gamma_only_local_ = gamma_only_local; const std::string kmesh_type_ = kmesh_type; + const int my_rank = GlobalV::MY_RANK; + const int my_pool = GlobalV::MY_POOL; - // (1) set nspin, read kpoints. - this->nspin = nspin_in; - ModuleBase::GlobalFunc::OUT(ofs, "nspin", nspin); + // (1) print nspin, set the k-point spin multiplicity, read kpoints. + ModuleBase::GlobalFunc::OUT(ofs, "nspin", nspin_in); - if (this->nspin != 1 && this->nspin != 2 && this->nspin != 4) + if (nspin_in != 1 && nspin_in != 2 && nspin_in != 4) { - ModuleBase::WARNING_QUIT("K_Vectors::set", "Only available for nspin = 1 or 2 or 4"); + ModuleBase::WARNING_QUIT("K_Vectors::set", "Only available for nspin 1, 2 or 4"); } - this->nspin = (this->nspin == 4) ? 1 : this->nspin; + // non-collinear (nspin=4) does not double the k-point list, so its + // k-point spin multiplicity is the same as for the unpolarized case. + this->spin_mult = (nspin_in == 4) ? 1 : nspin_in; - bool read_succesfully = this->read_kpoints(ucell, k_file_name, gamma_only_local_, kspacing, kmesh_type_, koffset, ofs); + bool read_succesfully = this->read_kpoints(ucell, + k_file_name, + gamma_only_local_, + kspacing, + kmesh_type_, + koffset, + ofs, + ofs_warning, + my_rank); #ifdef __MPI Parallel_Common::bcast_bool(read_succesfully); #endif @@ -96,16 +76,14 @@ void K_Vectors::set(const UnitCell& ucell, std::string skpt1; std::string skpt2; - if (!this->kc_done && this->kd_done) - { - for (size_t ik = 0; ik != this->nkstot_full; ++ik) - this->kvec_c_full[ik] = this->kvec_d[ik] * reciprocal_vec; - } - else if (this->kc_done && !this->kd_done) - { - for (size_t ik = 0; ik != this->nkstot_full; ++ik) - this->kvec_c_full[ik] = this->kvec_c[ik]; - } + // complement the Cartesian coordinates of the full k-point list + KListIO::fill_full_kvec(this->kc_done, + this->kd_done, + this->nkstot_nospin, + reciprocal_vec, + this->kvec_c, + this->kvec_d, + this->kvec_c_full); // (2) @@ -114,30 +92,13 @@ void K_Vectors::set(const UnitCell& ucell, { bool match = true; // calculate kpoints in IBZ and reduce kpoints according to symmetry - this->reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, match); + this->reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, match, my_rank, ofs); #ifdef __MPI Parallel_Common::bcast_bool(match); #endif if (!match) { - std::cout << "Optimized lattice type of reciprocal lattice cannot match the optimized real lattice. " - << std::endl; - std::cout << "It is often because the inaccuracy of lattice parameters in STRU." << std::endl; - if (ModuleSymmetry::Symmetry::symm_autoclose) - { - ModuleBase::WARNING("K_Vectors::ibz_kpoint", "Automatically set symmetry to 0 and continue ..."); - std::cout << "Automatically set symmetry to 0 and continue ..." << std::endl; - ModuleSymmetry::Symmetry::symm_flag = 0; - match = true; - this->reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, match); - } else { - ModuleBase::WARNING_QUIT("K_Vectors::ibz_kpoint", - "Possible solutions: \n \ -1. Refine the lattice parameters in STRU;\n \ -2. Use a different`symmetry_prec`. \n \ -3. Close symemtry: set `symmetry` to 0 in INPUT. \n \ -4. Set `symmetry_autoclose` to 1 in INPUT to automatically close symmetry when this error occurs."); - } + this->handle_symmetry_mismatch(ucell, symm, skpt1, match, my_rank, ofs); } } @@ -145,9 +106,9 @@ void K_Vectors::set(const UnitCell& ucell, // Improve k point information // Complement the coordinates of k point - this->set_both_kvec(reciprocal_vec, latvec, skpt2, ofs); + this->set_both_kvec(reciprocal_vec, latvec, skpt2, ofs, ofs_warning); - if (GlobalV::MY_RANK == 0) + if (my_rank == 0) { // output kpoints file std::stringstream skpt; @@ -166,27 +127,32 @@ void K_Vectors::set(const UnitCell& ucell, // do set_kup_and_kdw() this->para_k.kinfo(nkstot, GlobalV::KPAR, - GlobalV::MY_POOL, + my_pool, GlobalV::RANK_IN_POOL, GlobalV::NPROC, nspin_in); // assign k points to several process pools #ifdef __MPI // distribute K point data to the corresponding process - this->mpi_k(ofs); + this->mpi_k(ofs, my_rank, my_pool); #endif // set the k vectors for the up and down spin this->set_kup_and_kdw(ofs); // initialize ibz_index - this->ibz_index.resize(this->nkstot_full); - for (int ik = 0; ik < this->nkstot_full; ik++) + this->ibz_index.resize(this->nkstot_nospin); + for (int ik = 0; ik < this->nkstot_nospin; ik++) { this->ibz_index[ik] = ik; } - // get ik2iktot - this->cal_ik_global(); + // get ik2iktot: map local k indices to global indices in the pool + KListIO::build_ik2iktot(this->para_k.my_pool, + this->para_k.startk_pool, + this->spin_mult, + this->nks, + this->nkstot, + this->ik2iktot); this->print_klists(ofs); @@ -195,8 +161,37 @@ void K_Vectors::set(const UnitCell& ucell, return; } -// 1.reset the size of the K-point container according to nspin and nkstot -// 2.reserve space for nspin>2 (symmetry) +void K_Vectors::handle_symmetry_mismatch(const UnitCell& ucell, + const ModuleSymmetry::Symmetry& symm, + std::string& skpt, + bool& match, + const int my_rank, + std::ofstream& ofs) +{ + std::cout << "Optimized lattice type of reciprocal lattice cannot match the optimized real lattice. " + << std::endl; + std::cout << "It is often because the inaccuracy of lattice parameters in STRU." << std::endl; + if (ModuleSymmetry::Symmetry::symm_autoclose) + { + ModuleBase::WARNING("K_Vectors::ibz_kpoint", "Automatically set symmetry to 0 and continue ..."); + std::cout << "Automatically set symmetry to 0 and continue ..." << std::endl; + ModuleSymmetry::Symmetry::symm_flag = 0; + match = true; + this->reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match, my_rank, ofs); + } + else + { + ModuleBase::WARNING_QUIT("K_Vectors::ibz_kpoint", + "Possible solutions: \n \ +1. Refine the lattice parameters in STRU;\n \ +2. Use a different`symmetry_prec`. \n \ +3. Close symemtry: set `symmetry` to 0 in INPUT. \n \ +4. Set `symmetry_autoclose` to 1 in INPUT to automatically close symmetry when this error occurs."); + } +} + +// 1.reset the size of the K-point container according to spin_mult and nkstot +// 2.reserve space for spin_mult>2 (symmetry) void K_Vectors::renew(const int& kpoint_number) { ReciprocalGrid::renew(kpoint_number); @@ -213,80 +208,32 @@ bool K_Vectors::read_kpoints(const UnitCell& ucell, const double kspacing[3], const std::string& kmesh_type, const double koffset[3], - std::ofstream& ofs_running) + std::ofstream& ofs_running, + std::ofstream& ofs_warning, + const int my_rank) { ModuleBase::TITLE("K_Vectors", "read_kpoints"); - if (GlobalV::MY_RANK != 0) + if (my_rank != 0) { return true; } // 1. Overwrite the KPT file and default K-point information if needed // mohan add 2010-09-04 - this->generate_kfile(ucell, fn, gamma_only_local, kspacing, kmesh_type, koffset); + KListIO::write_auto_kfile(ucell, fn, gamma_only_local, kspacing, kmesh_type, koffset, ofs_warning); // 2. Read the KPT file and build the k-point list - return this->parse_kfile(fn, ofs_running); -} - -void K_Vectors::generate_kfile(const UnitCell& ucell, - const std::string& fn, - const bool gamma_only_local, - const double kspacing[3], - const std::string& kmesh_type, - const double koffset[3]) -{ - if (gamma_only_local) - { - GlobalV::ofs_warning << " Auto generating k-points file: " << fn << std::endl; - std::ofstream ofs(fn.c_str()); - ofs << "K_POINTS" << std::endl; - ofs << "0" << std::endl; - ofs << "Gamma" << std::endl; - ofs << "1 1 1 0 0 0" << std::endl; - ofs.close(); - } - else if (kspacing[0] > 0.0) - { - if (kspacing[1] <= 0 || kspacing[2] <= 0) - { - ModuleBase::WARNING_QUIT("K_Vectors", "kspacing should > 0"); - }; - // number of K points = max(1,int(|bi|/KSPACING+1)) - ModuleBase::Matrix3 btmp = ucell.G; - double b1 = sqrt(btmp.e11 * btmp.e11 + btmp.e12 * btmp.e12 + btmp.e13 * btmp.e13); - double b2 = sqrt(btmp.e21 * btmp.e21 + btmp.e22 * btmp.e22 + btmp.e23 * btmp.e23); - double b3 = sqrt(btmp.e31 * btmp.e31 + btmp.e32 * btmp.e32 + btmp.e33 * btmp.e33); - int nk1 = std::max(1, static_cast(b1 * ModuleBase::TWO_PI / kspacing[0] / ucell.lat0 + 1)); - int nk2 = std::max(1, static_cast(b2 * ModuleBase::TWO_PI / kspacing[1] / ucell.lat0 + 1)); - int nk3 = std::max(1, static_cast(b3 * ModuleBase::TWO_PI / kspacing[2] / ucell.lat0 + 1)); - - GlobalV::ofs_warning << " Generate k-points file according to KSPACING: " << fn << std::endl; - std::ofstream ofs(fn.c_str()); - ofs << "K_POINTS" << std::endl; - ofs << "0" << std::endl; - if (kmesh_type == "mp") - { - ofs << "Monkhorst-Pack" << std::endl; - } - else - { - ofs << "Gamma" << std::endl; - } - ofs << nk1 << " " << nk2 << " " << nk3 << " " << koffset[0] << " " << koffset[1] << " " - << koffset[2] << std::endl; - ofs.close(); - } + return this->parse_kfile(fn, ofs_running, ofs_warning); } // 2. Generate the K-point grid automatically according to the KPT file -bool K_Vectors::parse_kfile(const std::string& fn, std::ofstream& ofs_running) +bool K_Vectors::parse_kfile(const std::string& fn, std::ofstream& ofs_running, std::ofstream& ofs_warning) { // 2.1 read the KPT file std::ifstream ifk(fn.c_str()); if (!ifk) { - GlobalV::ofs_warning << " Can't find File name : " << fn << std::endl; + ofs_warning << " Can't find File name : " << fn << std::endl; return false; } @@ -295,29 +242,11 @@ bool K_Vectors::parse_kfile(const std::string& fn, std::ofstream& ofs_running) ifk.clear(); ifk.seekg(0); - std::string word; std::string kword; - int ierr = 0; - - ifk.rdstate(); - - while (ifk.good()) - { - ifk >> word; - ifk.ignore(150, '\n'); // LiuXh add 20180416, fix bug in k-point file when the first line with comments - if (word == "K_POINTS" || word == "KPOINTS" || word == "K") - { - ierr = 1; - break; - } - - ifk.rdstate(); - } - - if (ierr == 0) + if (!KListIO::find_kpoints_header(ifk)) { - GlobalV::ofs_warning << " symbol K_POINTS not found." << std::endl; + ofs_warning << " symbol K_POINTS not found." << std::endl; return false; } @@ -335,191 +264,151 @@ bool K_Vectors::parse_kfile(const std::string& fn, std::ofstream& ofs_running) const int max_kpoints = 100000; if (nkstot > max_kpoints) { - GlobalV::ofs_warning << " nkstot > MAX_KPOINTS" << std::endl; + ofs_warning << " nkstot > MAX_KPOINTS" << std::endl; return false; } // 2.2 Select different methods and generate K-point grid - int k_type = 0; + bool kpts_ok = true; if (nkstot == 0) // nkstot==0, use monkhorst_pack. add by dwan { - if (kword == "Gamma") // MP(Gamma) - { - is_mp = true; - k_type = 0; - ModuleBase::GlobalFunc::OUT(ofs_running, "Input type of k points", "Monkhorst-Pack(Gamma)"); - } - else if (kword == "Monkhorst-Pack" || kword == "MP" || kword == "mp") - { - is_mp = true; - k_type = 1; - ModuleBase::GlobalFunc::OUT(ofs_running, "Input type of k points", "Monkhorst-Pack"); - } - else - { - GlobalV::ofs_warning << " Error: neither Gamma nor Monkhorst-Pack." << std::endl; - return false; - } - - ifk >> nmp[0] >> nmp[1] >> nmp[2]; - - this->koffset[0] = 0; - this->koffset[1] = 0; - this->koffset[2] = 0; - if (!(ifk >> this->koffset[0] >> this->koffset[1] >> this->koffset[2])) - { - ModuleBase::WARNING("K_Vectors::read_kpoints", "Missing k-point offsets in the k-points file."); - } - - this->Monkhorst_Pack(nmp, this->koffset, k_type); + kpts_ok = this->read_mp_mesh(ifk, kword, ofs_running, ofs_warning); } else if (nkstot > 0) // nkstot>0, the K-point information is clearly set { - if (kword == "Cartesian" || kword == "C") // Cartesian coordinates - { - this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 - for (int i = 0; i < nkstot; i++) - { - ifk >> kvec_c[i].x >> kvec_c[i].y >> kvec_c[i].z; - ModuleBase::GlobalFunc::READ_VALUE(ifk, wk[i]); - } - - this->kc_done = true; - } - else if (kword == "Direct" || kword == "D") // Direct coordinates - { - this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 - for (int i = 0; i < nkstot; i++) - { - ifk >> kvec_d[i].x >> kvec_d[i].y >> kvec_d[i].z; - ModuleBase::GlobalFunc::READ_VALUE(ifk, wk[i]); - } - this->kd_done = true; - } - else if (kword == "Line_Cartesian") - { - if (ModuleSymmetry::Symmetry::symm_flag == 1) - { - ModuleBase::WARNING("K_Vectors::read_kpoints", - "Line mode of k-points is open, please set symmetry to 0 or -1."); - return false; - } - - interpolate_k_between(ifk, kvec_c); - - std::for_each(wk.begin(), wk.end(), [](double& d) { d = 1.0; }); + kpts_ok = this->read_listed_kpoints(ifk, kword, ofs_warning); + } - this->kc_done = true; - } + if (!kpts_ok) + { + return false; + } - else if (kword == "Line_Direct" || kword == "L" || kword == "Line") - { - if (ModuleSymmetry::Symmetry::symm_flag == 1) - { - ModuleBase::WARNING("K_Vectors::read_kpoints", - "Line mode of k-points is open, please set symmetry to 0 or -1."); - return false; - } + this->nkstot_nospin = this->nks = this->nkstot; - interpolate_k_between(ifk, kvec_d); + ModuleBase::GlobalFunc::OUT(ofs_running, "nkstot", nkstot); + return true; +} // END SUBROUTINE - std::for_each(wk.begin(), wk.end(), [](double& d) { d = 1.0; }); +bool K_Vectors::read_mp_mesh(std::ifstream& ifk, + const std::string& kword, + std::ofstream& ofs_running, + std::ofstream& ofs_warning) +{ + int k_type = 0; + if (kword == "Gamma") // MP(Gamma) + { + is_mp = true; + k_type = 0; + ModuleBase::GlobalFunc::OUT(ofs_running, "Input type of k points", "Monkhorst-Pack(Gamma)"); + } + else if (kword == "Monkhorst-Pack" || kword == "MP" || kword == "mp") + { + is_mp = true; + k_type = 1; + ModuleBase::GlobalFunc::OUT(ofs_running, "Input type of k points", "Monkhorst-Pack"); + } + else + { + ofs_warning << " Error: neither Gamma nor Monkhorst-Pack." << std::endl; + return false; + } - this->kd_done = true; - } + ifk >> nmp[0] >> nmp[1] >> nmp[2]; - else - { - GlobalV::ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; - return false; - } + this->koffset[0] = 0; + this->koffset[1] = 0; + this->koffset[2] = 0; + if (!(ifk >> this->koffset[0] >> this->koffset[1] >> this->koffset[2])) + { + ofs_warning << " K_Vectors::read_kpoints warning : " + << "Missing k-point offsets in the k-points file." << std::endl; } - this->nkstot_full = this->nks = this->nkstot; - - ModuleBase::GlobalFunc::OUT(ofs_running, "nkstot", nkstot); + this->Monkhorst_Pack(nmp, this->koffset, k_type); return true; -} // END SUBROUTINE +} -void K_Vectors::interpolate_k_between(std::ifstream& ifk, std::vector>& kvec) +bool K_Vectors::read_listed_kpoints(std::ifstream& ifk, const std::string& kword, std::ofstream& ofs_warning) { - // how many special points. - int nks_special = this->nkstot; - - // number of points to the next k points - std::vector nkl(nks_special, 0); - - // coordinates of special points. - std::vector> ks(nks_special); - - // recalculate nkstot. - nkstot = 0; - /* ISSUE#3482: to distinguish different kline segments */ - std::vector kpt_segids; - kl_segids.clear(); - kl_segids.shrink_to_fit(); - int kpt_segid = 0; - for (int iks = 0; iks < nks_special; iks++) - { - ifk >> ks[iks].x; - ifk >> ks[iks].y; - ifk >> ks[iks].z; - ModuleBase::GlobalFunc::READ_VALUE(ifk, nkl[iks]); - - if (nkl[iks] <= 0) - { - ModuleBase::WARNING_QUIT("K_Vectors::interpolate_k_between", - "Line-mode interpolation counts must be positive."); - } - nkstot += nkl[iks]; - /* ISSUE#3482: to distinguish different kline segments */ - if ((nkl[iks] == 1) && (iks != (nks_special - 1))) { - kpt_segid++; - } - kpt_segids.push_back(kpt_segid); + if (kword == "Cartesian" || kword == "C") // Cartesian coordinates + { + this->renew(nkstot * this->spin_mult); // mohan fix bug 2009-09-01 + KListIO::read_kpt_list(ifk, nkstot, this->kvec_c, this->wk); + this->kc_done = true; + return true; } - if (nkl[nks_special - 1] != 1) + if (kword == "Direct" || kword == "D") // Direct coordinates { - ModuleBase::WARNING_QUIT("K_Vectors::interpolate_k_between", - "The final line-mode k-point must have an interpolation count of 1."); + this->renew(nkstot * this->spin_mult); // mohan fix bug 2009-09-01 + KListIO::read_kpt_list(ifk, nkstot, this->kvec_d, this->wk); + this->kd_done = true; + return true; + } + if (kword == "Line_Cartesian") + { + return this->setup_line_kpoints(ifk, this->kvec_c, true, ofs_warning); + } + if (kword == "Line_Direct" || kword == "L" || kword == "Line") + { + return this->setup_line_kpoints(ifk, this->kvec_d, false, ofs_warning); } - // std::cout << " nkstot = " << nkstot << std::endl; - this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 + ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; + return false; +} - int count = 0; - for (int iks = 1; iks < nks_special; iks++) +bool K_Vectors::setup_line_kpoints(std::ifstream& ifk, + std::vector>& kvec, + const bool cartesian, + std::ofstream& ofs_warning) +{ + if (ModuleSymmetry::Symmetry::symm_flag == 1) { - double dxs = (ks[iks].x - ks[iks - 1].x) / nkl[iks - 1]; - double dys = (ks[iks].y - ks[iks - 1].y) / nkl[iks - 1]; - double dzs = (ks[iks].z - ks[iks - 1].z) / nkl[iks - 1]; - for (int is = 0; is < nkl[iks - 1]; is++) - { - kvec[count].x = ks[iks - 1].x + is * dxs; - kvec[count].y = ks[iks - 1].y + is * dys; - kvec[count].z = ks[iks - 1].z + is * dzs; - kl_segids.push_back(kpt_segids[iks - 1]); /* ISSUE#3482: to distinguish different kline segments */ - ++count; - } + ofs_warning << " K_Vectors::read_kpoints warning : " + << "Line mode of k-points is open, please set symmetry to 0 or -1." + << std::endl; + return false; } - // deal with the last special k point. - kvec[count].x = ks[nks_special - 1].x; - kvec[count].y = ks[nks_special - 1].y; - kvec[count].z = ks[nks_special - 1].z; - kl_segids.push_back(kpt_segids[nks_special - 1]); /* ISSUE#3482: to distinguish different kline segments */ - ++count; + this->interpolate_k_between(ifk, kvec); + + std::for_each(this->wk.begin(), this->wk.end(), [](double& d) { d = 1.0; }); - assert(count == nkstot); - assert(kl_segids.size() == nkstot); /* ISSUE#3482: to distinguish different kline segments */ + if (cartesian) + { + this->kc_done = true; + } + else + { + this->kd_done = true; + } + return true; +} + +void K_Vectors::interpolate_k_between(std::ifstream& ifk, std::vector>& kvec) +{ + // Thin wrapper: the interpolation itself is the this-free KListIO::interp_line; + // here we only size the member containers and copy the results back. + const KListIO::LineK line = KListIO::interp_line(ifk, this->nkstot); + + this->nkstot = line.nks_total; + this->renew(this->nkstot * this->spin_mult); // mohan fix bug 2009-09-01 + + for (int i = 0; i < this->nkstot; i++) + { + kvec[i] = line.kpts[i]; + } + this->kl_segids = line.segids; /* ISSUE#3482: to distinguish different kline segments */ } void K_Vectors::update_use_ibz(const int& nkstot_ibz, const std::vector>& kvec_d_ibz, const std::vector& wk_ibz, - std::ofstream& ofs_running) + std::ofstream& ofs_running, + const int my_rank) { - if (GlobalV::MY_RANK != 0) { + if (my_rank != 0) { return; } ModuleBase::TITLE("K_Vectors", "update_use_ibz"); @@ -530,7 +419,8 @@ void K_Vectors::update_use_ibz(const int& nkstot_ibz, ModuleBase::GlobalFunc::OUT(ofs_running, "nkstot now", nkstot); - this->kvec_d.resize(this->nkstot * nspin); // qianrui fix a bug 2021-7-13 for nspin=2 in set_kup_and_kdw() + // qianrui fix a bug 2021-7-13: size for the spin_mult=2 doubling in set_kup_and_kdw() + this->kvec_d.resize(this->nkstot * this->spin_mult); for (int i = 0; i < this->nkstot; ++i) { @@ -553,41 +443,18 @@ void K_Vectors::set_kup_and_kdw(std::ofstream& ofs_running) { ModuleBase::TITLE("K_Vectors", "setup_kup_and_kdw"); - //========================================================================= - // on output: the number of points is doubled and xk and wk in the - // first (nks/2) positions correspond to up spin - // those in the second (nks/2) ones correspond to down spin - // nspin can only be 1 or 2 here: K_Vectors::set() maps nspin=4 - // (non-collinear) to 1 before the k-list is built. - //========================================================================= - switch (nspin) - { - case 1: - - for (int ik = 0; ik < nks; ik++) - { - this->isk[ik] = 0; - } - - break; - - case 2: - - for (int ik = 0; ik < nks; ik++) - { - this->kvec_c[ik + nks] = kvec_c[ik]; - this->kvec_d[ik + nks] = kvec_d[ik]; - this->wk[ik + nks] = wk[ik]; - this->isk[ik] = 0; - this->isk[ik + nks] = 1; - } - - this->nks *= 2; - this->nkstot *= 2; + KListIO::expand_spin_kpoints(this->spin_mult, + this->kvec_c, + this->kvec_d, + this->wk, + this->isk, + this->nks, + this->nkstot); - ModuleBase::GlobalFunc::OUT(ofs_running, "nks(nspin=2)", nks); - ModuleBase::GlobalFunc::OUT(ofs_running, "nkstot(nspin=2)", nkstot); - break; + if (this->spin_mult == 2) + { + ModuleBase::GlobalFunc::OUT(ofs_running, "nks(nspin=2)", this->nks); + ModuleBase::GlobalFunc::OUT(ofs_running, "nkstot(nspin=2)", this->nkstot); } return; @@ -597,9 +464,11 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, const ModuleSymmetry::Symmetry& symm, bool use_symm, std::string& skpt, - bool& match) + bool& match, + const int my_rank, + std::ofstream& ofs_running) { - if (GlobalV::MY_RANK != 0) + if (my_rank != 0) { return; } @@ -610,9 +479,7 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, // if the operations does not already included // inverse operation, double it. //=============================================== - bool include_inv = false; std::vector kgmatrix(48 * 2); - ModuleBase::Matrix3 inv(-1, 0, 0, 0, -1, 0, 0, 0, -1); ModuleBase::Matrix3 k_vec; int nrotkm = 0; @@ -626,39 +493,9 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, return; } - // check whether the inverse operation is already included - for (int i = 0; i < nrotkm; ++i) - { - if (kgmatrix[i] == inv) - { - include_inv = true; - } - } - - if (symm.magnetic_nspin4) - { - // (nspin=4, magnetic) Time reversal Theta reverses the magnetization, so Theta alone is - // NOT a symmetry and the blanket "-k is always equivalent" doubling below is invalid. - // Only the antiunitary elements Theta*g with g in the moment-reversing coset belong to - // the Shubnikov group; append exactly those, keeping the index convention - // j + nrotk <-> Theta * gmatrix_anti[j] (decoded the same way in restore_dm). - // (nspin=2 is unaffected: there the antiunitary operation is plain conjugation K, which - // does not touch the spin, so D_s(-k)=D_s^*(k) holds even for a ferromagnet and the - // generic branch below stays correct.) - for (int j = 0; j < symm.nrotk_anti; ++j) - { - kgmatrix[j + symm.nrotk] = inv * symm.kgmatrix_anti[j]; - } - nrotkm = symm.nrotk + symm.nrotk_anti; - } - else if (!include_inv) - { - for (int i = 0; i < symm.nrotk; ++i) - { - kgmatrix[i + symm.nrotk] = inv * symm.kgmatrix[i]; - } - nrotkm = 2 * symm.nrotk; - } + // append time-reversal-related operations (Theta*g for magnetic + // nspin=4; -g otherwise unless inversion is already present) + nrotkm = KListIO::append_time_reversal_ops(symm, kgmatrix, nrotkm); // convert kgmatrix to k-lattice std::vector kkmatrix(nrotkm); @@ -672,103 +509,53 @@ void K_Vectors::reduce_by_symmetry(const UnitCell& ucell, std::vector> kvec_d_ibz; std::vector wk_ibz; std::vector ibz2bz; - this->reduce_ibz(kgmatrix.data(), nrotkm, ucell.G, k_vec, kkmatrix.data(), symm.epsilon, kvec_d_ibz, wk_ibz, this->ibz_index, ibz2bz); + this->reduce_ibz(kgmatrix.data(), + nrotkm, + ucell.G, + k_vec, + kkmatrix.data(), + symm.epsilon, + kvec_d_ibz, + wk_ibz, + this->ibz_index, + ibz2bz); const int nkstot_ibz = kvec_d_ibz.size(); #ifdef __EXX // setup kstars according to the final (max-norm) kvec_d_ibz - this->kstars.resize(nkstot_ibz); if (ModuleSymmetry::Symmetry::symm_flag == 1) { - ModuleBase::Vector3 kvec_rot; - for (int i = 0; i < this->nkstot; ++i) - { - int exist_number = -1; - int isym = 0; - for (int j = 0; j < nrotkm; ++j) - { - kvec_rot = this->kvec_d[i] * kgmatrix[j]; - ModuleCell::restrict_kpt(kvec_rot, symm.epsilon); - for (int k = 0; k < nkstot_ibz; ++k) - { - if (symm.equal(kvec_rot.x, kvec_d_ibz[k].x) && symm.equal(kvec_rot.y, kvec_d_ibz[k].y) - && symm.equal(kvec_rot.z, kvec_d_ibz[k].z)) - { - isym = j; - exist_number = k; - break; - } - } - if (exist_number != -1) - { - break; - } - } - this->kstars[exist_number].insert(std::make_pair(isym, this->kvec_d[i])); - } + KListIO::build_kstars(this->kvec_d, + kgmatrix, + nrotkm, + kvec_d_ibz, + symm.epsilon, + [&symm](double a, double b) { return symm.equal(a, b); }, + this->kstars); } #endif // output in kpoints file - std::stringstream ss; - ss << " " << std::setw(40) << "nkstot" - << " = " << this->nkstot << std::setw(66) << "ibzkpt" << std::endl; - std::string table; - table += "K-POINTS REDUCTION ACCORDING TO SYMMETRY\n"; - table += FmtCore::format("%8s%12s%12s%12s%8s%12s%12s%12s\n", - "KPT", - "DIRECT_X", - "DIRECT_Y", - "DIRECT_Z", - "IBZ", - "DIRECT_X", - "DIRECT_Y", - "DIRECT_Z"); - for (int i = 0; i < this->nkstot; ++i) - { - table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8d%12.8f%12.8f%12.8f\n", - i + 1, - this->kvec_d[i].x, - this->kvec_d[i].y, - this->kvec_d[i].z, - this->ibz_index[i] + 1, - kvec_d_ibz[this->ibz_index[i]].x, - kvec_d_ibz[this->ibz_index[i]].y, - kvec_d_ibz[this->ibz_index[i]].z); - } - ss << table << std::endl; - skpt = ss.str(); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Number of irreducible k-points", nkstot_ibz); - - table.clear(); - table += "\n K-POINTS REDUCTION ACCORDING TO SYMMETRY\n"; - table += FmtCore::format("%8s%12s%12s%12s%8s%8s\n", "IBZ", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT", "ibz2bz"); - for (int ik = 0; ik < nkstot_ibz; ik++) - { - table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f%8d\n", - ik + 1, - kvec_d_ibz[ik].x, - kvec_d_ibz[ik].y, - kvec_d_ibz[ik].z, - wk_ibz[ik], - ibz2bz[ik]); - } - GlobalV::ofs_running << table << std::endl; + skpt = KListIO::ibz_kpt_table(this->nkstot, this->kvec_d, this->ibz_index, kvec_d_ibz); + ModuleBase::GlobalFunc::OUT(ofs_running, "Number of irreducible k-points", nkstot_ibz); + + ofs_running << KListIO::ibz_wk_table(nkstot_ibz, kvec_d_ibz, wk_ibz, ibz2bz) << std::endl; // resize the kpoint container according to nkstot_ibz if (use_symm || this->get_is_mp()) { - this->update_use_ibz(nkstot_ibz, kvec_d_ibz, wk_ibz, GlobalV::ofs_running); + this->update_use_ibz(nkstot_ibz, kvec_d_ibz, wk_ibz, ofs_running, my_rank); } return; } -void K_Vectors::set_after_vc(const int& nspin_in, const ModuleBase::Matrix3& G, std::ofstream& ofs_running) +void K_Vectors::set_after_vc(const ModuleBase::Matrix3& G, std::ofstream& ofs_running) { ofs_running << "\n SETUP K-POINTS" << std::endl; - this->set_nspin(nspin_in); - ModuleBase::GlobalFunc::OUT(ofs_running, "nspin", this->get_nspin()); + // spin_mult is fixed by set() and does not change during a run, so the + // volume-change update only recomputes the Cartesian coordinates. + ModuleBase::GlobalFunc::OUT(ofs_running, "nspin", this->get_spin_mult()); // set cartesian k vectors. this->kvec_d2c(G); @@ -794,7 +581,7 @@ void K_Vectors::set_after_vc(const int& nspin_in, const ModuleBase::Matrix3& G, } #ifdef __MPI -void K_Vectors::mpi_k(std::ofstream& ofs_running) +void K_Vectors::mpi_k(std::ofstream& ofs_running, const int my_rank, const int my_pool) { ModuleBase::TITLE("K_Vectors", "mpi_k"); @@ -802,11 +589,11 @@ void K_Vectors::mpi_k(std::ofstream& ofs_running) Parallel_Common::bcast_bool(this->kd_done); - Parallel_Common::bcast_int(this->nspin); + Parallel_Common::bcast_int(this->spin_mult); Parallel_Common::bcast_int(this->nkstot); - Parallel_Common::bcast_int(this->nkstot_full); + Parallel_Common::bcast_int(this->nkstot_nospin); Parallel_Common::bcast_int(this->nmp, 3); @@ -815,7 +602,7 @@ void K_Vectors::mpi_k(std::ofstream& ofs_running) Parallel_Common::bcast_double(this->koffset, 3); - this->nks = this->para_k.nks_pool[GlobalV::MY_POOL]; + this->nks = this->para_k.nks_pool[my_pool]; ofs_running << std::endl; ModuleBase::GlobalFunc::OUT(ofs_running, "Number of k-points in this process", this->nks); @@ -836,25 +623,22 @@ void K_Vectors::mpi_k(std::ofstream& ofs_running) std::vector wk_aux(this->nkstot); std::vector kvec_c_aux(this->nkstot * 3); std::vector kvec_d_aux(this->nkstot * 3); - std::vector kvec_c_full_aux(this->nkstot_full * 3); + std::vector kvec_c_full_aux(this->nkstot_nospin * 3); // collect and process in rank 0 - if (GlobalV::MY_RANK == 0) + if (my_rank == 0) { - for (int ik = 0; ik < this->nkstot; ik++) - { - isk_aux[ik] = this->isk[ik]; - wk_aux[ik] = this->wk[ik]; - kvec_c_aux[3 * ik] = this->kvec_c[ik].x; - kvec_c_aux[3 * ik + 1] = this->kvec_c[ik].y; - kvec_c_aux[3 * ik + 2] = this->kvec_c[ik].z; - kvec_d_aux[3 * ik] = this->kvec_d[ik].x; - kvec_d_aux[3 * ik + 1] = this->kvec_d[ik].y; - kvec_d_aux[3 * ik + 2] = this->kvec_d[ik].z; - kvec_c_full_aux[3 * ik] = this->kvec_c_full[ik].x; - kvec_c_full_aux[3 * ik + 1] = this->kvec_c_full[ik].y; - kvec_c_full_aux[3 * ik + 2] = this->kvec_c_full[ik].z; - } + KListIO::pack_kpts(this->isk, + this->wk, + this->kvec_c, + this->kvec_d, + this->kvec_c_full, + this->nkstot, + isk_aux, + wk_aux, + kvec_c_aux, + kvec_d_aux, + kvec_c_full_aux); } // broadcast k point data to all processors @@ -863,60 +647,30 @@ void K_Vectors::mpi_k(std::ofstream& ofs_running) Parallel_Common::bcast_double(wk_aux.data(), this->nkstot); Parallel_Common::bcast_double(kvec_c_aux.data(), this->nkstot * 3); Parallel_Common::bcast_double(kvec_d_aux.data(), this->nkstot * 3); - Parallel_Common::bcast_double(kvec_c_full_aux.data(), this->nkstot_full * 3); + Parallel_Common::bcast_double(kvec_c_full_aux.data(), this->nkstot_nospin * 3); // process k point data in each processor - this->renew(this->nks * this->nspin); + this->renew(this->nks * this->spin_mult); // distribute - int k_index = 0; - - for (int i = 0; i < this->nks; i++) - { - // 3 is because each k point has three value:kx, ky, kz - k_index = i + this->para_k.startk_pool[GlobalV::MY_POOL]; - this->kvec_c[i].x = kvec_c_aux[k_index * 3]; - this->kvec_c[i].y = kvec_c_aux[k_index * 3 + 1]; - this->kvec_c[i].z = kvec_c_aux[k_index * 3 + 2]; - this->kvec_d[i].x = kvec_d_aux[k_index * 3]; - this->kvec_d[i].y = kvec_d_aux[k_index * 3 + 1]; - this->kvec_d[i].z = kvec_d_aux[k_index * 3 + 2]; - this->kvec_c_full[i].x = kvec_c_full_aux[k_index * 3]; - this->kvec_c_full[i].y = kvec_c_full_aux[k_index * 3 + 1]; - this->kvec_c_full[i].z = kvec_c_full_aux[k_index * 3 + 2]; - this->wk[i] = wk_aux[k_index]; - this->isk[i] = isk_aux[k_index]; - } + KListIO::unpack_kpts(isk_aux, + wk_aux, + kvec_c_aux, + kvec_d_aux, + kvec_c_full_aux, + this->nks, + this->para_k.startk_pool[my_pool], + this->isk, + this->wk, + this->kvec_c, + this->kvec_d, + this->kvec_c_full); #ifdef __EXX + // bcast kstars (rank 0 holds the filled maps; other ranks rebuild them) if (ModuleSymmetry::Symmetry::symm_flag == 1) - { // bcast kstars - this->kstars.resize(this->nkstot); - for (int ikibz = 0; ikibz < this->nkstot; ++ikibz) - { - int starsize = this->kstars[ikibz].size(); - Parallel_Common::bcast_int(starsize); - auto ks = this->kstars[ikibz].begin(); - for (int ik = 0; ik < starsize; ++ik) - { - int isym = 0; - ModuleBase::Vector3 ks_vec(0, 0, 0); - if (GlobalV::MY_RANK == 0) - { - isym = ks->first; - ks_vec = ks->second; - ++ks; - } - Parallel_Common::bcast_int(isym); - Parallel_Common::bcast_double(ks_vec.x); - Parallel_Common::bcast_double(ks_vec.y); - Parallel_Common::bcast_double(ks_vec.z); - if (GlobalV::MY_RANK != 0) - { - this->kstars[ikibz].insert(std::make_pair(isym, ks_vec)); - } - } - } + { + KListIO::bcast_kstars(this->kstars, this->nkstot, my_rank); } #endif } // END SUBROUTINE mpi_k diff --git a/source/source_cell/klist.h b/source/source_cell/klist.h index 2e212629f06..956a2be794f 100644 --- a/source/source_cell/klist.h +++ b/source/source_cell/klist.h @@ -14,7 +14,7 @@ * Inherits the spin-free common reciprocal-grid functionality * (mesh generation, coordinate conversion, weights, printing, star/IBZ * reduction primitive) from ModuleCell::ReciprocalGrid and adds the - * spin expansion (isk, nspin doubling) and the k-point IBZ logic. + * spin expansion (isk, spin-multiplicity doubling) and the k-point IBZ logic. */ class K_Vectors : public ModuleCell::ReciprocalGrid { @@ -60,6 +60,7 @@ class K_Vectors : public ModuleCell::ReciprocalGrid const ModuleBase::Matrix3& reciprocal_vec, const ModuleBase::Matrix3& latvec, std::ofstream& ofs, + std::ofstream& ofs_warning, const bool use_ibz, const std::string& global_out_dir, const bool gamma_only_local, @@ -77,9 +78,9 @@ class K_Vectors : public ModuleCell::ReciprocalGrid return this->nkstot; } - int get_nkstot_full() const + int get_nkstot_nospin() const { - return this->nkstot_full; + return this->nkstot_nospin; } double get_koffset(const int i) const @@ -92,9 +93,11 @@ class K_Vectors : public ModuleCell::ReciprocalGrid return this->k_nkstot; } - int get_nspin() const + /// @brief Spin multiplicity of the k-point list: 1 (no doubling, also for + /// non-collinear nspin=4) or 2 (LSDA, k points split into up/down). + int get_spin_mult() const { - return this->nspin; + return this->spin_mult; } std::string get_k_kword() const @@ -112,14 +115,9 @@ class K_Vectors : public ModuleCell::ReciprocalGrid this->nkstot = value; } - void set_nkstot_full(int value) + void set_nkstot_nospin(int value) { - this->nkstot_full = value; - } - - void set_nspin(int value) - { - this->nspin = value; + this->nkstot_nospin = value; } bool get_is_mp() const @@ -150,24 +148,27 @@ class K_Vectors : public ModuleCell::ReciprocalGrid void update_use_ibz(const int& nkstot_ibz, const std::vector>& kvec_d_ibz, const std::vector& wk_ibz, - std::ofstream& ofs_running); + std::ofstream& ofs_running, + const int my_rank); /** - * @brief Sets up the k-points after a volume change. + * @brief Updates the k-points after a volume change. * - * Sets the number of spins, converts the direct coordinates (which are - * kept across the volume change) to the new Cartesian coordinates using - * the new reciprocal lattice, prints the resulting table, and marks both - * coordinate sets as up to date. + * Converts the direct coordinates (which are kept across the volume + * change) to the new Cartesian coordinates using the new reciprocal + * lattice, prints the resulting table, and marks both coordinate sets + * as up to date. The spin multiplicity is not touched: it was fixed by + * set() and never changes during a run. * - * @param nspin_in The number of spins. 1 for non-spin-polarized - * calculations and 2 for spin-polarized calculations. * @param G The new reciprocal lattice matrix. */ - void set_after_vc(const int& nspin_in, const ModuleBase::Matrix3& G, std::ofstream& ofs_running); + void set_after_vc(const ModuleBase::Matrix3& G, std::ofstream& ofs_running); private: - int nspin = 0; ///< number of spin states + /// Spin multiplicity used to size the k-point list: 1 for input nspin 1 + /// or 4 (non-collinear k points are not doubled) and 2 for input nspin 2 + /// (LSDA up/down k points). This is NOT the physical nspin (1/2/4). + int spin_mult = 0; double koffset[3] = {0.0}; ///< used only in automatic k-points /** @@ -182,10 +183,10 @@ class K_Vectors : public ModuleCell::ReciprocalGrid */ void renew(const int& kpoint_number) override; - /// @brief Spin multiplicity used when generating the mesh (1/2 for nspin 1/2). + /// @brief Spin multiplicity used when generating the mesh (1/2). int spin_factor() const override { - return this->nspin; + return this->spin_mult; } /** @@ -207,7 +208,9 @@ class K_Vectors : public ModuleCell::ReciprocalGrid const ModuleSymmetry::Symmetry& symm, bool use_symm, std::string& skpt, - bool& match) override; + bool& match, + const int my_rank, + std::ofstream& ofs_running) override; /// @brief step 1 : generate kpoints @@ -236,28 +239,9 @@ class K_Vectors : public ModuleCell::ReciprocalGrid const double kspacing[3], const std::string& kmesh_type, const double koffset[3], - std::ofstream& ofs_running); // return 0: something wrong. - - /** - * @brief Overwrite the KPT file with an auto-generated mesh when requested. - * - * Writes a Gamma-mesh KPT file if gamma_only_local is set, or a - * KSPACING-derived Gamma/Monkhorst-Pack mesh if kspacing is positive. - * Does nothing when neither condition holds. - * - * @param ucell unit cell (reciprocal lattice and lat0 for the mesh size) - * @param fn KPT filename to (over)write - * @param gamma_only_local whether to force a single Gamma point - * @param kspacing target k-point spacing in 1/bohr (three components) - * @param kmesh_type "mp" for Monkhorst-Pack, anything else for Gamma - * @param koffset mesh offsets (three components) - */ - void generate_kfile(const UnitCell& ucell, - const std::string& fn, - const bool gamma_only_local, - const double kspacing[3], - const std::string& kmesh_type, - const double koffset[3]); + std::ofstream& ofs_running, + std::ofstream& ofs_warning, + const int my_rank); // return 0: something wrong. /** * @brief Read the KPT file and build the k-point list from it. @@ -271,7 +255,77 @@ class K_Vectors : public ModuleCell::ReciprocalGrid * @return bool Returns true if the k-points are successfully read, * false otherwise. */ - bool parse_kfile(const std::string& fn, std::ofstream& ofs_running); + bool parse_kfile(const std::string& fn, std::ofstream& ofs_running, std::ofstream& ofs_warning); + + /** + * @brief Read the Monkhorst-Pack/Gamma mesh block and generate the mesh. + * + * Handles the nkstot == 0 form of the KPT file: validates the type + * keyword, reads the mesh dimensions and optional offsets, then calls + * Monkhorst_Pack to fill the k-point list. + * + * @param ifk stream positioned after the type keyword + * @param kword type keyword (Gamma / Monkhorst-Pack / MP / mp) + * @param ofs_running running log stream + * @return false (after warning) when the keyword is neither Gamma nor + * Monkhorst-Pack; true when the mesh was generated. + */ + bool read_mp_mesh(std::ifstream& ifk, + const std::string& kword, + std::ofstream& ofs_running, + std::ofstream& ofs_warning); + + /** + * @brief Read the explicitly listed k points (nkstot > 0 form of KPT). + * + * Dispatches on the type keyword: Cartesian/Direct lists are sized via + * renew() and filled through KListIO::read_kpt_list; Line_Cartesian/ + * Line_Direct delegate to setup_line_kpoints. + * + * @param ifk stream positioned after the type keyword + * @param kword type keyword: Cartesian, C, Direct, D, Line_Cartesian, + * Line_Direct, L or Line + * @return false (after warning) for unknown keywords or line mode with + * symmetry enabled; true when the k-point list was built. + */ + bool read_listed_kpoints(std::ifstream& ifk, const std::string& kword, std::ofstream& ofs_warning); + + /** + * @brief Build line-mode k points by interpolating between special points. + * + * Refuses (warning + false) when symmetry reduction is enabled, then + * interpolates the special points read from `ifk`, resets all weights + * to 1, and marks the Cartesian or Direct coordinate set as done. + * + * @param ifk stream to read the special points from + * @param kvec target coordinate container (kvec_c or kvec_d) + * @param cartesian true for Line_Cartesian, false for Line_Direct + * @param ofs_warning warning-log stream for error messages + */ + bool setup_line_kpoints(std::ifstream& ifk, + std::vector>& kvec, + const bool cartesian, + std::ofstream& ofs_warning); + + /** + * @brief Handle a reciprocal/real lattice Bravais-type mismatch after + * IBZ reduction. + * + * When symmetry_autoclose is enabled, symmetry is switched off and the + * IBZ reduction is retried; otherwise the run aborts with a WARNING_QUIT + * listing the possible remedies. + * + * @param ucell unit cell used for the retried IBZ reduction + * @param symm symmetry operations used for the retried reduction + * @param skpt k-point option string forwarded to reduce_by_symmetry + * @param match set to true when the autoclose retry succeeds + */ + void handle_symmetry_mismatch(const UnitCell& ucell, + const ModuleSymmetry::Symmetry& symm, + std::string& skpt, + bool& match, + const int my_rank, + std::ofstream& ofs); /** * @brief Adds k-points linearly between special points. @@ -319,12 +373,6 @@ class K_Vectors : public ModuleCell::ReciprocalGrid */ void set_kup_and_kdw(std::ofstream& ofs_running); - /** - * @brief Gets the global index of a k-point. - * @return this->ik2iktot[ik] - */ - void cal_ik_global(); - #ifdef __MPI /** * @brief Distributes k-points among MPI processes. @@ -336,7 +384,7 @@ class K_Vectors : public ModuleCell::ReciprocalGrid * @note Assumes nkstot > 0 and quits if some process ends up with * no k-points. */ - void mpi_k(std::ofstream& ofs_running); + void mpi_k(std::ofstream& ofs_running, const int my_rank, const int my_pool); #endif }; #endif // KVECT_H \ No newline at end of file diff --git a/source/source_cell/klist_io.cpp b/source/source_cell/klist_io.cpp new file mode 100644 index 00000000000..f869a0c7baf --- /dev/null +++ b/source/source_cell/klist_io.cpp @@ -0,0 +1,494 @@ +/** + * @file klist_io.cpp + * @brief this-free helpers extracted from K_Vectors (IBZ table formatting and + * line-mode k-point interpolation). Kept separate from klist.cpp so the + * logic is testable in isolation; klist.cpp only keeps thin wrappers. + */ +#include "klist_io.h" + +#include "source_base/formatter.h" +#include "source_base/global_function.h" +#include "source_base/parallel_common.h" +#include "source_cell/module_symmetry/symmetry.h" +#include "source_cell/reciprocal_grid.h" +#include "source_cell/unitcell.h" + +#include +#include +#include + +namespace KListIO +{ + +std::string ibz_kpt_table(const int nkstot, + const std::vector>& kvec_d, + const std::vector& ibz_index, + const std::vector>& kvec_d_ibz) +{ + std::stringstream ss; + ss << " " << std::setw(40) << "nkstot" + << " = " << nkstot << std::setw(66) << "ibzkpt" << std::endl; + std::string table; + table += "K-POINTS REDUCTION ACCORDING TO SYMMETRY\n"; + table += FmtCore::format("%8s%12s%12s%12s%8s%12s%12s%12s\n", + "KPT", + "DIRECT_X", + "DIRECT_Y", + "DIRECT_Z", + "IBZ", + "DIRECT_X", + "DIRECT_Y", + "DIRECT_Z"); + for (int i = 0; i < nkstot; ++i) + { + table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8d%12.8f%12.8f%12.8f\n", + i + 1, + kvec_d[i].x, + kvec_d[i].y, + kvec_d[i].z, + ibz_index[i] + 1, + kvec_d_ibz[ibz_index[i]].x, + kvec_d_ibz[ibz_index[i]].y, + kvec_d_ibz[ibz_index[i]].z); + } + ss << table << std::endl; + return ss.str(); +} + +std::string ibz_wk_table(const int nkstot_ibz, + const std::vector>& kvec_d_ibz, + const std::vector& wk_ibz, + const std::vector& ibz2bz) +{ + std::string table; + table += "\n K-POINTS REDUCTION ACCORDING TO SYMMETRY\n"; + table += FmtCore::format("%8s%12s%12s%12s%8s%8s\n", "IBZ", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT", "ibz2bz"); + for (int ik = 0; ik < nkstot_ibz; ik++) + { + table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f%8d\n", + ik + 1, + kvec_d_ibz[ik].x, + kvec_d_ibz[ik].y, + kvec_d_ibz[ik].z, + wk_ibz[ik], + ibz2bz[ik]); + } + return table; +} + +bool find_kpoints_header(std::ifstream& ifk) +{ + std::string word; + while (ifk.good()) + { + ifk >> word; + // LiuXh add 20180416, fix bug in k-point file when the first line with comments + ifk.ignore(150, '\n'); + if (word == "K_POINTS" || word == "KPOINTS" || word == "K") + { + return true; + } + } + return false; +} + +void read_kpt_list(std::ifstream& ifk, + const int nkstot, + std::vector>& kvec, + std::vector& wk) +{ + for (int i = 0; i < nkstot; i++) + { + ifk >> kvec[i].x >> kvec[i].y >> kvec[i].z; + ModuleBase::GlobalFunc::READ_VALUE(ifk, wk[i]); + } +} + +LineK interp_line(std::ifstream& ifk, const int nks_special) +{ + // number of points to the next k points + std::vector nkl(nks_special, 0); + + // coordinates of special points. + std::vector> ks(nks_special); + + LineK out; + std::vector kpt_segids; + int kpt_segid = 0; + for (int iks = 0; iks < nks_special; iks++) + { + ifk >> ks[iks].x; + ifk >> ks[iks].y; + ifk >> ks[iks].z; + ModuleBase::GlobalFunc::READ_VALUE(ifk, nkl[iks]); + + if (nkl[iks] <= 0) + { + ModuleBase::WARNING_QUIT("KListIO::interp_line", + "Line-mode interpolation counts must be positive."); + } + out.nks_total += nkl[iks]; + /* ISSUE#3482: to distinguish different kline segments */ + if ((nkl[iks] == 1) && (iks != (nks_special - 1))) { + kpt_segid++; + } + kpt_segids.push_back(kpt_segid); + } + if (nkl[nks_special - 1] != 1) + { + ModuleBase::WARNING_QUIT("KListIO::interp_line", + "The final line-mode k-point must have an interpolation count of 1."); + } + + out.kpts.resize(out.nks_total); + out.segids.reserve(out.nks_total); + + int count = 0; + for (int iks = 1; iks < nks_special; iks++) + { + double dxs = (ks[iks].x - ks[iks - 1].x) / nkl[iks - 1]; + double dys = (ks[iks].y - ks[iks - 1].y) / nkl[iks - 1]; + double dzs = (ks[iks].z - ks[iks - 1].z) / nkl[iks - 1]; + for (int is = 0; is < nkl[iks - 1]; is++) + { + out.kpts[count].x = ks[iks - 1].x + is * dxs; + out.kpts[count].y = ks[iks - 1].y + is * dys; + out.kpts[count].z = ks[iks - 1].z + is * dzs; + out.segids.push_back(kpt_segids[iks - 1]); /* ISSUE#3482 */ + ++count; + } + } + + // deal with the last special k point. + out.kpts[count].x = ks[nks_special - 1].x; + out.kpts[count].y = ks[nks_special - 1].y; + out.kpts[count].z = ks[nks_special - 1].z; + out.segids.push_back(kpt_segids[nks_special - 1]); /* ISSUE#3482 */ + ++count; + + assert(count == out.nks_total); + assert(out.segids.size() == static_cast(out.nks_total)); /* ISSUE#3482 */ + return out; +} + +void build_kstars(const std::vector>& kvec_d, + const std::vector& kgmatrix, + const int nrotkm, + const std::vector>& kvec_d_ibz, + const double epsilon, + const std::function& equal, + std::vector>>& kstars) +{ + const int nkstot = static_cast(kvec_d.size()); + const int nkstot_ibz = static_cast(kvec_d_ibz.size()); + kstars.resize(nkstot_ibz); + + ModuleBase::Vector3 kvec_rot; + for (int i = 0; i < nkstot; ++i) + { + int exist_number = -1; + int isym = 0; + for (int j = 0; j < nrotkm; ++j) + { + kvec_rot = kvec_d[i] * kgmatrix[j]; + ModuleCell::restrict_kpt(kvec_rot, epsilon); + for (int k = 0; k < nkstot_ibz; ++k) + { + if (equal(kvec_rot.x, kvec_d_ibz[k].x) && equal(kvec_rot.y, kvec_d_ibz[k].y) + && equal(kvec_rot.z, kvec_d_ibz[k].z)) + { + isym = j; + exist_number = k; + break; + } + } + if (exist_number != -1) + { + break; + } + } + kstars[exist_number].insert(std::make_pair(isym, kvec_d[i])); + } +} + +int append_time_reversal_ops(const ModuleSymmetry::Symmetry& symm, + std::vector& kgmatrix, + const int nrotkm) +{ + const ModuleBase::Matrix3 inv{-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0}; + + bool include_inv = false; + for (int i = 0; i < nrotkm; ++i) + { + if (kgmatrix[i] == inv) + { + include_inv = true; + } + } + + if (symm.magnetic_nspin4) + { + // (nspin=4, magnetic) Time reversal Theta reverses the magnetization, + // so Theta alone is NOT a symmetry; only the antiunitary Theta*g + // elements with g in the moment-reversing coset belong to the + // Shubnikov group. The same index convention j + nrotk is decoded + // in restore_dm. (nspin=2 is unaffected: there the antiunitary + // operation is plain conjugation K, which leaves D_s(-k)=D_s*(k).) + for (int j = 0; j < symm.nrotk_anti; ++j) + { + kgmatrix[j + symm.nrotk] = inv * symm.kgmatrix_anti[j]; + } + return symm.nrotk + symm.nrotk_anti; + } + if (!include_inv) + { + for (int i = 0; i < symm.nrotk; ++i) + { + kgmatrix[i + symm.nrotk] = inv * symm.kgmatrix[i]; + } + return 2 * symm.nrotk; + } + return nrotkm; +} + +void pack_kpts(const std::vector& isk, + const std::vector& wk, + const std::vector>& kvec_c, + const std::vector>& kvec_d, + const std::vector>& kvec_c_full, + const int nkstot, + std::vector& isk_aux, + std::vector& wk_aux, + std::vector& kvec_c_aux, + std::vector& kvec_d_aux, + std::vector& kvec_c_full_aux) +{ + for (int ik = 0; ik < nkstot; ik++) + { + isk_aux[ik] = isk[ik]; + wk_aux[ik] = wk[ik]; + kvec_c_aux[3 * ik] = kvec_c[ik].x; + kvec_c_aux[3 * ik + 1] = kvec_c[ik].y; + kvec_c_aux[3 * ik + 2] = kvec_c[ik].z; + kvec_d_aux[3 * ik] = kvec_d[ik].x; + kvec_d_aux[3 * ik + 1] = kvec_d[ik].y; + kvec_d_aux[3 * ik + 2] = kvec_d[ik].z; + kvec_c_full_aux[3 * ik] = kvec_c_full[ik].x; + kvec_c_full_aux[3 * ik + 1] = kvec_c_full[ik].y; + kvec_c_full_aux[3 * ik + 2] = kvec_c_full[ik].z; + } +} + +void unpack_kpts(const std::vector& isk_aux, + const std::vector& wk_aux, + const std::vector& kvec_c_aux, + const std::vector& kvec_d_aux, + const std::vector& kvec_c_full_aux, + const int nks, + const int startk, + std::vector& isk, + std::vector& wk, + std::vector>& kvec_c, + std::vector>& kvec_d, + std::vector>& kvec_c_full) +{ + for (int i = 0; i < nks; i++) + { + // 3 is because each k point has three value:kx, ky, kz + const int k_index = i + startk; + kvec_c[i].x = kvec_c_aux[k_index * 3]; + kvec_c[i].y = kvec_c_aux[k_index * 3 + 1]; + kvec_c[i].z = kvec_c_aux[k_index * 3 + 2]; + kvec_d[i].x = kvec_d_aux[k_index * 3]; + kvec_d[i].y = kvec_d_aux[k_index * 3 + 1]; + kvec_d[i].z = kvec_d_aux[k_index * 3 + 2]; + kvec_c_full[i].x = kvec_c_full_aux[k_index * 3]; + kvec_c_full[i].y = kvec_c_full_aux[k_index * 3 + 1]; + kvec_c_full[i].z = kvec_c_full_aux[k_index * 3 + 2]; + wk[i] = wk_aux[k_index]; + isk[i] = isk_aux[k_index]; + } +} + +void bcast_kstars(std::vector>>& kstars, + const int nkstot, + const int my_rank) +{ + kstars.resize(nkstot); + for (int ikibz = 0; ikibz < nkstot; ++ikibz) + { + int starsize = kstars[ikibz].size(); + Parallel_Common::bcast_int(starsize); + auto ks = kstars[ikibz].begin(); + for (int ik = 0; ik < starsize; ++ik) + { + int isym = 0; + ModuleBase::Vector3 ks_vec(0, 0, 0); + if (my_rank == 0) + { + isym = ks->first; + ks_vec = ks->second; + ++ks; + } + Parallel_Common::bcast_int(isym); + Parallel_Common::bcast_double(ks_vec.x); + Parallel_Common::bcast_double(ks_vec.y); + Parallel_Common::bcast_double(ks_vec.z); + if (my_rank != 0) + { + kstars[ikibz].insert(std::make_pair(isym, ks_vec)); + } + } + } +} + +void fill_full_kvec(const bool kc_done, + const bool kd_done, + const int nkstot_nospin, + const ModuleBase::Matrix3& reciprocal_vec, + const std::vector>& kvec_c, + const std::vector>& kvec_d, + std::vector>& kvec_c_full) +{ + if (!kc_done && kd_done) + { + for (int ik = 0; ik < nkstot_nospin; ++ik) + { + kvec_c_full[ik] = kvec_d[ik] * reciprocal_vec; + } + } + else if (kc_done && !kd_done) + { + for (int ik = 0; ik < nkstot_nospin; ++ik) + { + kvec_c_full[ik] = kvec_c[ik]; + } + } +} + +void build_ik2iktot(const int my_pool, + const std::vector& startk_pool, + const int spin_mult, + const int nks, + const int nkstot, + std::vector& ik2iktot) +{ + ik2iktot.resize(nks); +#ifdef __MPI + if (spin_mult == 2) + { + for (int ik = 0; ik < nks / 2; ++ik) + { + ik2iktot[ik] = startk_pool[my_pool] + ik; + ik2iktot[ik + nks / 2] = nkstot / 2 + startk_pool[my_pool] + ik; + } + } + else + { + for (int ik = 0; ik < nks; ++ik) + { + ik2iktot[ik] = startk_pool[my_pool] + ik; + } + } +#else + for (int ik = 0; ik < nks; ++ik) + { + ik2iktot[ik] = ik; + } +#endif +} + +void expand_spin_kpoints(const int spin_mult, + std::vector>& kvec_c, + std::vector>& kvec_d, + std::vector& wk, + std::vector& isk, + int& nks, + int& nkstot) +{ + //========================================================================= + // on output: the number of points is doubled and xk and wk in the + // first (nks/2) positions correspond to up spin + // those in the second (nks/2) ones correspond to down spin + // spin_mult can only be 1 or 2 here: K_Vectors::set() maps nspin=4 + // (non-collinear) to 1 before the k-list is built. + //========================================================================= + switch (spin_mult) + { + case 1: + for (int ik = 0; ik < nks; ik++) + { + isk[ik] = 0; + } + break; + + case 2: + for (int ik = 0; ik < nks; ik++) + { + kvec_c[ik + nks] = kvec_c[ik]; + kvec_d[ik + nks] = kvec_d[ik]; + wk[ik + nks] = wk[ik]; + isk[ik] = 0; + isk[ik + nks] = 1; + } + + nks *= 2; + nkstot *= 2; + break; + } + + return; +} + +void write_auto_kfile(const UnitCell& ucell, + const std::string& fn, + const bool gamma_only_local, + const double kspacing[3], + const std::string& kmesh_type, + const double koffset[3], + std::ofstream& ofs_warning) +{ + if (gamma_only_local) + { + ofs_warning << " Auto generating k-points file: " << fn << std::endl; + std::ofstream ofs(fn.c_str()); + ofs << "K_POINTS" << std::endl; + ofs << "0" << std::endl; + ofs << "Gamma" << std::endl; + ofs << "1 1 1 0 0 0" << std::endl; + ofs.close(); + } + else if (kspacing[0] > 0.0) + { + if (kspacing[1] <= 0 || kspacing[2] <= 0) + { + ModuleBase::WARNING_QUIT("K_Vectors", "kspacing should > 0"); + }; + // number of K points = max(1,int(|bi|/KSPACING+1)) + ModuleBase::Matrix3 btmp = ucell.G; + double b1 = sqrt(btmp.e11 * btmp.e11 + btmp.e12 * btmp.e12 + btmp.e13 * btmp.e13); + double b2 = sqrt(btmp.e21 * btmp.e21 + btmp.e22 * btmp.e22 + btmp.e23 * btmp.e23); + double b3 = sqrt(btmp.e31 * btmp.e31 + btmp.e32 * btmp.e32 + btmp.e33 * btmp.e33); + int nk1 = std::max(1, static_cast(b1 * ModuleBase::TWO_PI / kspacing[0] / ucell.lat0 + 1)); + int nk2 = std::max(1, static_cast(b2 * ModuleBase::TWO_PI / kspacing[1] / ucell.lat0 + 1)); + int nk3 = std::max(1, static_cast(b3 * ModuleBase::TWO_PI / kspacing[2] / ucell.lat0 + 1)); + + ofs_warning << " Generate k-points file according to KSPACING: " << fn << std::endl; + std::ofstream ofs(fn.c_str()); + ofs << "K_POINTS" << std::endl; + ofs << "0" << std::endl; + if (kmesh_type == "mp") + { + ofs << "Monkhorst-Pack" << std::endl; + } + else + { + ofs << "Gamma" << std::endl; + } + ofs << nk1 << " " << nk2 << " " << nk3 << " " << koffset[0] << " " << koffset[1] << " " + << koffset[2] << std::endl; + ofs.close(); + } +} + +} // namespace KListIO diff --git a/source/source_cell/klist_io.h b/source/source_cell/klist_io.h new file mode 100644 index 00000000000..c9785b0f85e --- /dev/null +++ b/source/source_cell/klist_io.h @@ -0,0 +1,175 @@ +#ifndef KLIST_IO_H +#define KLIST_IO_H + +#include "source_base/matrix3.h" +#include "source_base/vector3.h" + +#include +#include +#include +#include +#include + +namespace ModuleSymmetry +{ +class Symmetry; // full definition only needed in klist_io.cpp +} + +class UnitCell; // full definition only needed in klist_io.cpp + +/// this-free helpers extracted from K_Vectors, kept in a separate TU so they +/// can be unit-tested and reused without dragging in the K_Vectors class. +namespace KListIO +{ +/// Render the IBZ reduction table ("IBZ" k-point -> originating k-point). +std::string ibz_kpt_table(int nkstot, + const std::vector>& kvec_d, + const std::vector& ibz_index, + const std::vector>& kvec_d_ibz); + +/// Render the IBZ weight table (IBZ k-point, weight, multiplicity, origin index). +std::string ibz_wk_table(int nkstot_ibz, + const std::vector>& kvec_d_ibz, + const std::vector& wk_ibz, + const std::vector& ibz2bz); + +/// Result of line-mode interpolation between special k-points. +struct LineK +{ + std::vector> kpts; ///< interpolated k points + std::vector segids; ///< segment id per k point (ISSUE#3482) + int nks_total = 0; ///< total interpolated k-point count +}; + +/// Scan `ifk` for the "K_POINTS"/"KPOINTS"/"K" header keyword, skipping any +/// leading comment lines. Returns true with the stream positioned after the +/// header line; returns false if the keyword is not found before EOF. +bool find_kpoints_header(std::ifstream& ifk); + +/// Read `nkstot` explicit k points (three coordinates plus a weight per line) +/// from `ifk` into `kvec` and `wk`. The caller is responsible for sizing the +/// arrays (K_Vectors::renew) before calling. +void read_kpt_list(std::ifstream& ifk, + int nkstot, + std::vector>& kvec, + std::vector& wk); + +/// Read the special k points and per-point interpolation counts from `ifk`, +/// then linearly interpolate the line-mode k points. Pure function of the +/// stream and `nks_special`; dies via WARNING_QUIT on malformed input. +LineK interp_line(std::ifstream& ifk, int nks_special); + +/// Build the EXX k-stars: for every k point, find the symmetry operation +/// (index into `kgmatrix`) that rotates it onto an irreducible k point, and +/// group k points by that IBZ representative. `equal` compares two doubles +/// with the symmetry precision; `epsilon` is the k-restriction tolerance. +/// this-free so the heavy triple loop is isolated and testable. +void build_kstars(const std::vector>& kvec_d, + const std::vector& kgmatrix, + int nrotkm, + const std::vector>& kvec_d_ibz, + double epsilon, + const std::function& equal, + std::vector>>& kstars); + +/// Append the time-reversal-related k-point symmetry operations into +/// `kgmatrix` (the slots right after the first `nrotkm` operations must be +/// available). For magnetic nspin=4 systems the antiunitary Theta*g coset +/// is appended from `symm.kgmatrix_anti`; otherwise the inverted -g ops are +/// appended unless inversion is already present. Returns the updated total +/// operation count. +int append_time_reversal_ops(const ModuleSymmetry::Symmetry& symm, + std::vector& kgmatrix, + int nrotkm); + +/// Flatten k-point arrays into contiguous MPI buffers (x,y,z interleaved). +/// this-free; used on rank 0 before broadcasting in K_Vectors::mpi_k. +void pack_kpts(const std::vector& isk, + const std::vector& wk, + const std::vector>& kvec_c, + const std::vector>& kvec_d, + const std::vector>& kvec_c_full, + int nkstot, + std::vector& isk_aux, + std::vector& wk_aux, + std::vector& kvec_c_aux, + std::vector& kvec_d_aux, + std::vector& kvec_c_full_aux); + +/// Broadcast the EXX k-stars (one (symmetry-index, k-vector) map per IBZ +/// k-point) from `my_rank == 0` to every process. Rank 0 holds the filled +/// maps; other ranks resize and rebuild them from the broadcast. MPI +/// wrappers are compiled as no-ops without __MPI, so the call is safe in +/// serial builds (it simply leaves the rank-0 maps untouched). +void bcast_kstars(std::vector>>& kstars, + int nkstot, + int my_rank); + +/// Scatter the broadcast buffers into this pool's k-point slice, starting at +/// global index `startk`. this-free; mirrors pack_kpts after the broadcast. +void unpack_kpts(const std::vector& isk_aux, + const std::vector& wk_aux, + const std::vector& kvec_c_aux, + const std::vector& kvec_d_aux, + const std::vector& kvec_c_full_aux, + int nks, + int startk, + std::vector& isk, + std::vector& wk, + std::vector>& kvec_c, + std::vector>& kvec_d, + std::vector>& kvec_c_full); + +/// Fill the full-list Cartesian k vectors when only one coordinate set is +/// available: direct coordinates are converted via `reciprocal_vec` when +/// Cartesian points are missing, otherwise the Cartesian points are copied. +/// No-op when both coordinate sets are done. this-free helper called from +/// K_Vectors::set() before IBZ reduction. +void fill_full_kvec(bool kc_done, + bool kd_done, + int nkstot_nospin, + const ModuleBase::Matrix3& reciprocal_vec, + const std::vector>& kvec_c, + const std::vector>& kvec_d, + std::vector>& kvec_c_full); + +/// Build the local-to-global k-point index map `ik2iktot` for this pool. +/// In MPI runs the global index is offset by the pool start (with the +/// spin_mult == 2 second half offset by nkstot/2); in serial runs it is the +/// local index itself. `ik2iktot` is resized to `nks` here. this-free helper +/// called from K_Vectors::set() after the pool distribution. +void build_ik2iktot(int my_pool, + const std::vector& startk_pool, + int spin_mult, + int nks, + int nkstot, + std::vector& ik2iktot); + +/// Expand the k-point list for spin-polarized runs (spin_mult == 2): copy +/// coordinates and weights into the second half, tag isk 0 for the first +/// half and 1 for the second, then double nks/nkstot. For spin_mult == 1 +/// only isk is zeroed. The running-log output stays in the K_Vectors +/// wrapper. this-free helper backing K_Vectors::set_kup_and_kdw. +void expand_spin_kpoints(int spin_mult, + std::vector>& kvec_c, + std::vector>& kvec_d, + std::vector& wk, + std::vector& isk, + int& nks, + int& nkstot); + +/// Overwrite the KPT file with an auto-generated mesh when requested: +/// a single Gamma point if gamma_only_local, or a KSPACING-derived +/// Gamma/Monkhorst-Pack mesh if kspacing[0] > 0 (quits if kspacing[1] or +/// kspacing[2] is non-positive); does nothing otherwise. this-free helper +/// backing K_Vectors::generate_kfile. +void write_auto_kfile(const UnitCell& ucell, + const std::string& fn, + bool gamma_only_local, + const double kspacing[3], + const std::string& kmesh_type, + const double koffset[3], + std::ofstream& ofs_warning); +} // namespace KListIO + +#endif // KLIST_IO_H diff --git a/source/source_cell/qlist.cpp b/source/source_cell/qlist.cpp index 128ca8a6279..c228bfa2af6 100644 --- a/source/source_cell/qlist.cpp +++ b/source/source_cell/qlist.cpp @@ -36,19 +36,19 @@ void QList::generate_mesh(UnitCell& ucell, ModuleSymmetry::Symmetry& symm, const double offset[3] = {0.0, 0.0, 0.0}; this->Monkhorst_Pack(this->nmp, offset, 0); - this->nkstot_full = this->nkstot; + this->nkstot_nospin = this->nkstot; this->nks = this->nkstot; // Star reduction: always use symmetry, always include the -q partner. bool match = true; std::string skpt; - this->reduce_by_symmetry(ucell, symm, true, skpt, match); + this->reduce_by_symmetry(ucell, symm, true, skpt, match, GlobalV::MY_RANK, GlobalV::ofs_running); if (!match) { ModuleBase::WARNING("QList::generate_mesh", "Reciprocal lattice is incompatible with the real-space lattice. " "Falling back to the unreduced q-point mesh."); - this->nkstot = this->nks = this->nkstot_full; + this->nkstot = this->nks = this->nkstot_nospin; } // weights sum to 1 (average over the full Brillouin zone) @@ -189,7 +189,7 @@ void QList::read_from_file(const std::string& filename, UnitCell& ucell) { } } - this->nkstot_full = this->nks = this->nkstot; + this->nkstot_nospin = this->nks = this->nkstot; // complement the coordinates: fill the missing representation if (!this->kc_done && this->kd_done) @@ -329,8 +329,12 @@ void QList::reduce_by_symmetry(const UnitCell& ucell, const ModuleSymmetry::Symmetry& symm, bool use_symm, std::string& skpt, - bool& match) { + bool& match, + const int my_rank, + std::ofstream& ofs_running) { (void)skpt; + (void)my_rank; + (void)ofs_running; // q-points are spin-free: build the point-group operations and always // double them by the time-reversal operation -q (no magnetic group). std::vector kgmatrix(48 * 2); diff --git a/source/source_cell/qlist.h b/source/source_cell/qlist.h index cfe5642b810..54844be7298 100644 --- a/source/source_cell/qlist.h +++ b/source/source_cell/qlist.h @@ -135,7 +135,9 @@ class QList : public ModuleCell::ReciprocalGrid { const ModuleSymmetry::Symmetry& symm, bool use_symm, std::string& skpt, - bool& match) override; + bool& match, + const int my_rank, + std::ofstream& ofs_running) override; private: std::vector nirr_; ///< number of irreps for each q-point diff --git a/source/source_cell/reciprocal_grid.cpp b/source/source_cell/reciprocal_grid.cpp index 6f94a8769cc..43015c4f265 100644 --- a/source/source_cell/reciprocal_grid.cpp +++ b/source/source_cell/reciprocal_grid.cpp @@ -170,7 +170,8 @@ void ReciprocalGrid::kvec_c2d(const ModuleBase::Matrix3& latvec) void ReciprocalGrid::set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt, - std::ofstream& ofs_running) + std::ofstream& ofs_running, + std::ofstream& ofs_warning) { // Re-derive the "which representation was read from file" flags. // For auto-generated meshes (k_nkstot == 0) the direct coordinates @@ -194,7 +195,7 @@ void ReciprocalGrid::set_both_kvec(const ModuleBase::Matrix3& G, } else { - GlobalV::ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; + ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; } } @@ -355,7 +356,7 @@ void ReciprocalGrid::reduce_ibz(const ModuleBase::Matrix3* rot_ops, ModuleBase::Vector3 kvec_rot_k; // update map k -> irreducible k - ibz_index.assign(this->nkstot_full, -1); // -1 means not in ibz list + ibz_index.assign(this->nkstot_nospin, -1); // -1 means not in ibz list // search in all k-points. for (int i = 0; i < this->nkstot; ++i) { diff --git a/source/source_cell/reciprocal_grid.h b/source/source_cell/reciprocal_grid.h index 020ca79e73c..91e3909d5c4 100644 --- a/source/source_cell/reciprocal_grid.h +++ b/source/source_cell/reciprocal_grid.h @@ -68,10 +68,13 @@ class ReciprocalGrid /// Number of points in the current pool (spin-free view). int nks = 0; - /// Total number of (symmetry-reduced) points. + /// Total number of (symmetry-reduced) points, INCLUDING spin multiplicity + /// (i.e. nkstot = nkstot_nospin * spin_mult after K_Vectors::set_kup_and_kdw). int nkstot = 0; - /// Total number of points before symmetry reduction. - int nkstot_full = 0; + /// Total number of physical k-points before symmetry reduction, + /// WITHOUT spin multiplicity. EXX/RI/LR code relies on this convention + /// (see e.g. ri_2d_comm.hpp: ik_full + is_k * nkstot_nospin). + int nkstot_nospin = 0; ReciprocalGrid() = default; virtual ~ReciprocalGrid() = default; @@ -108,7 +111,8 @@ class ReciprocalGrid void set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt, - std::ofstream& ofs_running); + std::ofstream& ofs_running, + std::ofstream& ofs_warning); /// @brief Normalize the weights so that they sum to the spin degeneracy. void normalize_wk(const int& degspin); @@ -157,7 +161,9 @@ class ReciprocalGrid const ModuleSymmetry::Symmetry& symm, bool use_symm, std::string& skpt, - bool& match) = 0; + bool& match, + const int my_rank, + std::ofstream& ofs_running) = 0; /// Whether this is a Monkhorst-Pack grid. bool is_mp = false; diff --git a/source/source_cell/test/CMakeLists.txt b/source/source_cell/test/CMakeLists.txt index caab63da4ad..653eaa6d3f1 100644 --- a/source/source_cell/test/CMakeLists.txt +++ b/source/source_cell/test/CMakeLists.txt @@ -86,7 +86,7 @@ AddTest( AddTest( TARGET MODULE_CELL_klist_test LIBS base device symmetry - SOURCES klist_test.cpp ../klist.cpp ../parallel_kpoints.cpp ../reciprocal_grid.cpp + SOURCES klist_test.cpp ../klist.cpp ../klist_io.cpp ../parallel_kpoints.cpp ../reciprocal_grid.cpp ) AddTest( @@ -110,7 +110,7 @@ AddTest( AddTest( TARGET MODULE_CELL_klist_test_para1 LIBS base device symmetry - SOURCES klist_test_para.cpp ../klist.cpp ../parallel_kpoints.cpp ../reciprocal_grid.cpp + SOURCES klist_test_para.cpp ../klist.cpp ../klist_io.cpp ../parallel_kpoints.cpp ../reciprocal_grid.cpp ) add_test(NAME MODULE_CELL_klist_test_para4 diff --git a/source/source_cell/test/klist_test.cpp b/source/source_cell/test/klist_test.cpp index 74ea3d12ec4..4f25675cb4d 100644 --- a/source/source_cell/test/klist_test.cpp +++ b/source/source_cell/test/klist_test.cpp @@ -133,6 +133,8 @@ class KlistTest : public testing::Test std::ifstream ifs; std::ofstream ofs; std::ofstream ofs_running; + std::ofstream ofs_warning; + int my_rank = 0; std::string output; // used to construct cell and analyse its symmetry @@ -214,7 +216,7 @@ TEST_F(KlistTest, Construct) { EXPECT_EQ(kv->get_nks(), 0); EXPECT_EQ(kv->get_nkstot(), 0); - EXPECT_EQ(kv->nspin, 0); + EXPECT_EQ(kv->spin_mult, 0); EXPECT_EQ(kv->k_nkstot, 0); EXPECT_FALSE(kv->kc_done); EXPECT_FALSE(kv->kd_done); @@ -229,7 +231,7 @@ TEST_F(KlistTest, MP) kv->koffset[0] = 0; kv->koffset[1] = 0; kv->koffset[2] = 0; - kv->nspin = 1; + kv->spin_mult = 1; int k_type = 0; kv->Monkhorst_Pack(kv->nmp, kv->koffset, k_type); /* @@ -246,7 +248,7 @@ TEST_F(KlistTest, MP) kv1->koffset[0] = 1; kv1->koffset[1] = 1; kv1->koffset[2] = 1; - kv1->nspin = 1; + kv1->spin_mult = 1; k_type = 1; kv1->Monkhorst_Pack(kv1->nmp, kv1->koffset, k_type); // std::cout << " " <nspin = 1; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); ifs.open("KPT_GO"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("Gamma")); @@ -278,46 +280,46 @@ TEST_F(KlistTest, ReadKpointsGammaOnlyLocal) TEST_F(KlistTest, ReadKpointsKspacing) { - kv->nspin = 1; + kv->spin_mult = 1; const bool gamma_only_local = false; const double kspacing[3] = {0.052918, 0.052918, 0.052918}; // 0.52918/Bohr = 1/A const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; setucell(); std::string k_file = "./support/KPT3"; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 343); } TEST_F(KlistTest, ReadKpointsKspacing3values) { - kv->nspin = 1; + kv->spin_mult = 1; const bool gamma_only_local = false; const double kspacing[3] = {0.052918, 0.06, 0.07}; // 0.52918/Bohr = 1/A const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; setucell(); std::string k_file = "./support/KPT3"; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 210); } TEST_F(KlistTest, ReadKpointsInvalidKspacing3values) { - kv->nspin = 1; + kv->spin_mult = 1; const bool gamma_only_local = false; const double kspacing[3] = {0.052918, 0.0, 0.07}; // 0.52918/Bohr = 1/A const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "./support/KPT3"; testing::internal::CaptureStdout(); - EXPECT_EXIT(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank), ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); } TEST_F(KlistTest, ReadKpointsKspacingShiftedGamma) { - kv->nspin = 1; + kv->spin_mult = 1; const bool gamma_only_local = false; const double kspacing[3] = {0.052918, 0.052918, 0.052918}; // 0.52918/Bohr = 1/A const std::string kmesh_type = "gamma"; @@ -325,7 +327,7 @@ TEST_F(KlistTest, ReadKpointsKspacingShiftedGamma) setucell(); std::string k_file = "./support/KPT3"; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 343); EXPECT_EQ(kv->get_k_kword(), "Gamma"); @@ -339,7 +341,7 @@ TEST_F(KlistTest, ReadKpointsKspacingShiftedGamma) TEST_F(KlistTest, ReadKpointsKspacingShiftedMP) { - kv->nspin = 1; + kv->spin_mult = 1; const bool gamma_only_local = false; const double kspacing[3] = {0.052918, 0.052918, 0.052918}; // 0.52918/Bohr = 1/A const std::string kmesh_type = "mp"; @@ -347,7 +349,7 @@ TEST_F(KlistTest, ReadKpointsKspacingShiftedMP) setucell(); std::string k_file = "./support/KPT3"; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 343); EXPECT_EQ(kv->get_k_kword(), "Monkhorst-Pack"); @@ -366,8 +368,8 @@ TEST_F(KlistTest, ReadKpointsGamma) const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "./support/KPT"; - kv->nspin = 1; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 512); } @@ -378,8 +380,8 @@ TEST_F(KlistTest, ReadKpointsMP) const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "./support/KPT1"; - kv->nspin = 1; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 512); } @@ -391,8 +393,8 @@ TEST_F(KlistTest, ReadKpointsLine) const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "./support/KPT2"; - kv->nspin = 1; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 122); } @@ -404,9 +406,9 @@ TEST_F(KlistTest, ReadKpointsLineRejectsZeroInterpolationCount) const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; const std::string k_file = "./support/KPT_ZERO_LINE_COUNT"; - kv->nspin = 1; + kv->spin_mult = 1; - EXPECT_EXIT(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running), + EXPECT_EXIT(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank), ::testing::ExitedWithCode(1), ""); } @@ -419,12 +421,12 @@ TEST_F(KlistTest, ReadKpointsCartesian) const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "./support/KPT4"; // Cartesian: non-spin case nspin=1 - kv->nspin = 1; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->kvec_c.size(), 5); // spin case nspin=2 - kv->nspin = 2; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 2; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->kvec_c.size(), 10); } @@ -436,14 +438,14 @@ TEST_F(KlistTest, ReadKpointsLineCartesian) const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "./support/KPT5"; // Line Cartesian: non-spin case nspin=1 - kv->nspin = 1; - kv->set_kup_and_kdw(GlobalV::ofs_running); - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->set_kup_and_kdw(ofs_running); + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 51); EXPECT_EQ(kv->kvec_c.size(), 51); // Line Cartesian: spin case nspin=2 - kv->nspin = 2; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 2; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 51); EXPECT_EQ(kv->kvec_c.size(), 102); } @@ -455,9 +457,9 @@ TEST_F(KlistTest, ReadKpointsDirect) const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "./support/KPT6"; - kv->nspin = 1; - kv->set_kup_and_kdw(GlobalV::ofs_running); - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->set_kup_and_kdw(ofs_running); + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 6); EXPECT_TRUE(kv->kd_done); } @@ -469,10 +471,10 @@ TEST_F(KlistTest, ReadKpointsWarning1) const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "arbitrary_1"; - kv->nspin = 1; - GlobalV::ofs_warning.open("klist_tmp_warning_1"); - EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running)); - GlobalV::ofs_warning.close(); + kv->spin_mult = 1; + ofs_warning.open("klist_tmp_warning_1"); + EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank)); + ofs_warning.close(); ifs.open("klist_tmp_warning_1"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("Can't find File name : arbitrary_1")); @@ -490,10 +492,10 @@ TEST_F(KlistTest, ReadKpointsWarning2) ofs.open(k_file.c_str()); ofs << "ARBITRARY"; ofs.close(); - kv->nspin = 1; - GlobalV::ofs_warning.open("klist_tmp_warning_2"); - EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running)); - GlobalV::ofs_warning.close(); + kv->spin_mult = 1; + ofs_warning.open("klist_tmp_warning_2"); + EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank)); + ofs_warning.close(); ifs.open("klist_tmp_warning_2"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("symbol K_POINTS not found.")); @@ -513,10 +515,10 @@ TEST_F(KlistTest, ReadKpointsWarning3) ofs << "KPOINTS" << std::endl; ofs << "100001" << std::endl; ofs.close(); - kv->nspin = 1; - GlobalV::ofs_warning.open("klist_tmp_warning_3"); - EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running)); - GlobalV::ofs_warning.close(); + kv->spin_mult = 1; + ofs_warning.open("klist_tmp_warning_3"); + EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank)); + ofs_warning.close(); ifs.open("klist_tmp_warning_3"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("nkstot > MAX_KPOINTS")); @@ -537,10 +539,10 @@ TEST_F(KlistTest, ReadKpointsWarning4) ofs << "0" << std::endl; ofs << "arbitrary" << std::endl; ofs.close(); - kv->nspin = 1; - GlobalV::ofs_warning.open("klist_tmp_warning_4"); - EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running)); - GlobalV::ofs_warning.close(); + kv->spin_mult = 1; + ofs_warning.open("klist_tmp_warning_4"); + EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank)); + ofs_warning.close(); ifs.open("klist_tmp_warning_4"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("Error: neither Gamma nor Monkhorst-Pack.")); @@ -561,10 +563,10 @@ TEST_F(KlistTest, ReadKpointsWarning5) ofs << "100000" << std::endl; ofs << "arbitrary" << std::endl; ofs.close(); - kv->nspin = 1; - GlobalV::ofs_warning.open("klist_tmp_warning_5"); - EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running)); - GlobalV::ofs_warning.close(); + kv->spin_mult = 1; + ofs_warning.open("klist_tmp_warning_5"); + EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank)); + ofs_warning.close(); ifs.open("klist_tmp_warning_5"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("Error : neither Cartesian nor Direct kpoint")); @@ -585,11 +587,11 @@ TEST_F(KlistTest, ReadKpointsWarning6) ofs << "100000" << std::endl; ofs << "Line_Cartesian" << std::endl; ofs.close(); - kv->nspin = 1; + kv->spin_mult = 1; ModuleSymmetry::Symmetry::symm_flag = 1; - GlobalV::ofs_warning.open("klist_tmp_warning_6"); - EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running)); - GlobalV::ofs_warning.close(); + ofs_warning.open("klist_tmp_warning_6"); + EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank)); + ofs_warning.close(); ifs.open("klist_tmp_warning_6"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("Line mode of k-points is open, please set symmetry to 0 or -1")); @@ -611,11 +613,11 @@ TEST_F(KlistTest, ReadKpointsWarning7) ofs << "100000" << std::endl; ofs << "Line_Direct" << std::endl; ofs.close(); - kv->nspin = 1; + kv->spin_mult = 1; ModuleSymmetry::Symmetry::symm_flag = 1; - GlobalV::ofs_warning.open("klist_tmp_warning_7"); - EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running)); - GlobalV::ofs_warning.close(); + ofs_warning.open("klist_tmp_warning_7"); + EXPECT_NO_THROW(kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank)); + ofs_warning.close(); ifs.open("klist_tmp_warning_7"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("Line mode of k-points is open, please set symmetry to 0 or -1")); @@ -632,26 +634,30 @@ TEST_F(KlistTest, SetKupKdown) const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; std::string k_file = "./support/KPT4"; - kv->nspin = 1; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); - kv->set_kup_and_kdw(GlobalV::ofs_running); + + // case A: physical nspin=1 -> spin_mult=1 (no doubling). + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); + kv->set_kup_and_kdw(ofs_running); for (int ik = 0; ik < 5; ik++) { EXPECT_EQ(kv->isk[ik], 0); } - kv->nspin = 4; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); - kv->set_kup_and_kdw(GlobalV::ofs_running); + + // case B: physical nspin=4 (non-collinear) maps to spin_mult=1 at + // K_Vectors::set() time; non-collinear does not double the k-point list, + // so the correct spin_mult is still 1. We bypass set() here, so set the + // mapped value directly. + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); + kv->set_kup_and_kdw(ofs_running); for (int ik = 0; ik < 5; ik++) { EXPECT_EQ(kv->isk[ik], 0); - EXPECT_EQ(kv->isk[ik + 5], 0); - EXPECT_EQ(kv->isk[ik + 10], 0); - EXPECT_EQ(kv->isk[ik + 15], 0); } - kv->nspin = 2; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); - kv->set_kup_and_kdw(GlobalV::ofs_running); + kv->spin_mult = 2; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); + kv->set_kup_and_kdw(ofs_running); for (int ik = 0; ik < 5; ik++) { EXPECT_EQ(kv->isk[ik], 0); @@ -661,44 +667,44 @@ TEST_F(KlistTest, SetKupKdown) TEST_F(KlistTest, SetAfterVC) { - kv->nspin = 1; + kv->spin_mult = 1; kv->set_nkstot(1); - GlobalV::ofs_running.open("tmp_klist_1"); + ofs_running.open("tmp_klist_1"); kv->renew(kv->get_nkstot()); kv->kvec_c[0].x = 0; kv->kvec_c[0].y = 0; kv->kvec_c[0].z = 0; - kv->set_after_vc(kv->nspin, ucell.G, GlobalV::ofs_running); + kv->set_after_vc(ucell.G, ofs_running); EXPECT_TRUE(kv->kd_done); EXPECT_TRUE(kv->kc_done); EXPECT_DOUBLE_EQ(kv->kvec_d[0].x, 0); EXPECT_DOUBLE_EQ(kv->kvec_d[0].y, 0); EXPECT_DOUBLE_EQ(kv->kvec_d[0].z, 0); - GlobalV::ofs_running.close(); + ofs_running.close(); remove("tmp_klist_1"); } TEST_F(KlistTest, PrintKlists) { - kv->nspin = 1; + kv->spin_mult = 1; kv->set_nkstot(1); kv->set_nks(1); - GlobalV::ofs_running.open("tmp_klist_2"); + ofs_running.open("tmp_klist_2"); kv->renew(kv->get_nkstot()); kv->kvec_c[0].x = 0; kv->kvec_c[0].y = 0; kv->kvec_c[0].z = 0; - kv->set_after_vc(kv->nspin, ucell.G, GlobalV::ofs_running); + kv->set_after_vc(ucell.G, ofs_running); EXPECT_TRUE(kv->kd_done); - kv->print_klists(GlobalV::ofs_running); - GlobalV::ofs_running.close(); + kv->print_klists(ofs_running); + ofs_running.close(); remove("tmp_klist_2"); } TEST_F(KlistTest, PrintKlistsWarnigQuit) { - kv->nspin = 1; + kv->spin_mult = 1; kv->set_nkstot(1); kv->set_nks(2); kv->renew(kv->get_nkstot()); @@ -706,14 +712,14 @@ TEST_F(KlistTest, PrintKlistsWarnigQuit) kv->kvec_c[0].y = 0; kv->kvec_c[0].z = 0; testing::internal::CaptureStdout(); - EXPECT_EXIT(kv->print_klists(GlobalV::ofs_running), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(kv->print_klists(ofs_running), ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("nkstot < nks")); } TEST_F(KlistTest, SetBothKvecFlagsFromFile) { - kv->nspin = 1; + kv->spin_mult = 1; kv->set_nkstot(1); kv->set_nks(1); kv->renew(kv->get_nkstot()); @@ -728,29 +734,29 @@ TEST_F(KlistTest, SetBothKvecFlagsFromFile) kv->kc_done = false; // case 1 kv->k_nkstot = 0; -// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); - kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); + kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); EXPECT_TRUE(kv->kd_done); EXPECT_TRUE(kv->kc_done); // case 2 kv->k_nkstot = 1; kv->k_kword = "D"; -// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); - kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); + kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); EXPECT_TRUE(kv->kd_done); EXPECT_TRUE(kv->kc_done); // case 3 kv->k_kword = "C"; -// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); - kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); + kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); EXPECT_TRUE(kv->kc_done); EXPECT_TRUE(kv->kd_done); // case 4 - GlobalV::ofs_warning.open("klist_tmp_warning_8"); + ofs_warning.open("klist_tmp_warning_8"); kv->k_kword = "arbitrary"; -// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); - kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); - GlobalV::ofs_warning.close(); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); + kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); + ofs_warning.close(); ifs.open("klist_tmp_warning_8"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); EXPECT_THAT(str, testing::HasSubstr("Error : neither Cartesian nor Direct kpoint.")); @@ -760,7 +766,7 @@ TEST_F(KlistTest, SetBothKvecFlagsFromFile) TEST_F(KlistTest, SetBothKvec) { - kv->nspin = 1; + kv->spin_mult = 1; kv->set_nkstot(1); kv->set_nks(1); kv->renew(kv->get_nkstot()); @@ -770,19 +776,19 @@ TEST_F(KlistTest, SetBothKvec) kv->kc_done = false; kv->kd_done = true; std::string skpt; -// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); - kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); + kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); EXPECT_TRUE(kv->kc_done); kv->kc_done = true; kv->kd_done = false; -// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); - kv->set_both_kvec(ucell.G, ucell.latvec, skpt, GlobalV::ofs_running); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); + kv->set_both_kvec(ucell.G, ucell.latvec, skpt, ofs_running, ofs_warning); EXPECT_TRUE(kv->kd_done); } TEST_F(KlistTest, NormalizeWk) { - kv->nspin = 1; + kv->spin_mult = 1; kv->set_nkstot(2); kv->set_nks(2); kv->renew(kv->get_nkstot()); @@ -797,7 +803,7 @@ TEST_F(KlistTest, NormalizeWk) TEST_F(KlistTest, NormalizeWkZeroWeights) { // Test that zero weights are handled correctly - kv->nspin = 1; + kv->spin_mult = 1; kv->set_nkstot(3); kv->set_nks(3); kv->renew(kv->get_nkstot()); @@ -822,11 +828,11 @@ TEST_F(KlistTest, NormalizeWkZeroWeights) TEST_F(KlistTest, UpdateUseIBZ) { - kv->nspin = 1; + kv->spin_mult = 1; kv->set_nkstot(3); kv->set_nks(3); kv->renew(kv->get_nkstot()); - kv->update_use_ibz(2, std::vector>(2, {0, 0, 0}), std::vector(2, 0.0), GlobalV::ofs_running); + kv->update_use_ibz(2, std::vector>(2, {0, 0, 0}), std::vector(2, 0.0), ofs_running, my_rank); EXPECT_EQ(kv->get_nkstot(), 2); EXPECT_EQ(kv->kvec_d.size(), 2); EXPECT_TRUE(kv->kd_done); @@ -841,21 +847,21 @@ TEST_F(KlistTest, IbzKpoint) const double koffset[3] = {0.0, 0.0, 0.0}; ModuleSymmetry::Symmetry symm; construct_ucell(stru_lib[0]); - GlobalV::ofs_running.open("tmp_klist_3"); + ofs_running.open("tmp_klist_3"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string k_file = "./support/KPT1"; - kv->nspin = 1; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 512); // calculate ibz_kpoint std::string skpt; ModuleSymmetry::Symmetry::symm_flag = 1; bool match = true; - kv->reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match); + kv->reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match, my_rank, ofs_running); EXPECT_EQ(kv->get_nkstot(), 35); - GlobalV::ofs_running << skpt << std::endl; - GlobalV::ofs_running.close(); + ofs_running << skpt << std::endl; + ofs_running.close(); ClearUcell(); remove("tmp_klist_3"); } @@ -868,22 +874,22 @@ TEST_F(KlistTest, IbzKpointIsMP) const double koffset[3] = {0.0, 0.0, 0.0}; ModuleSymmetry::Symmetry symm; construct_ucell(stru_lib[0]); - GlobalV::ofs_running.open("tmp_klist_4"); + ofs_running.open("tmp_klist_4"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); std::string k_file = "./support/KPT1"; - kv->nspin = 1; - kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->spin_mult = 1; + kv->read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv->get_nkstot(), 512); EXPECT_TRUE(kv->is_mp); // calculate ibz_kpoint std::string skpt; ModuleSymmetry::Symmetry::symm_flag = 0; bool match = true; - kv->reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match); + kv->reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match, my_rank, ofs_running); EXPECT_EQ(kv->get_nks(), 260); - GlobalV::ofs_running << skpt << std::endl; - GlobalV::ofs_running.close(); + ofs_running << skpt << std::endl; + ofs_running.close(); ClearUcell(); remove("tmp_klist_4"); } @@ -896,16 +902,16 @@ TEST_F(KlistTest, IbzKpointCustomWeights) const double koffset[3] = {0.0, 0.0, 0.0}; ModuleSymmetry::Symmetry symm; construct_ucell(stru_lib[0]); - GlobalV::ofs_running.open("tmp_klist_custom_weights"); + ofs_running.open("tmp_klist_custom_weights"); const int cal_symm_repr[2] = {0, 6}; - symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); + symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, ofs_running, 1e-6, 1, "scf", cal_symm_repr); // Test 1: Non-MP k-points with uniform weights (KPT4) { K_Vectors kv_test1; std::string k_file = "./support/KPT4"; - kv_test1.nspin = 1; - kv_test1.read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv_test1.spin_mult = 1; + kv_test1.read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv_test1.get_nkstot(), 5); EXPECT_FALSE(kv_test1.is_mp); // Should be non-MP @@ -916,7 +922,7 @@ TEST_F(KlistTest, IbzKpointCustomWeights) std::string skpt; ModuleSymmetry::Symmetry::symm_flag = 1; bool match = true; - kv_test1.reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match); + kv_test1.reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match, my_rank, ofs_running); // Verify that weights are preserved (not overwritten with 1/nkstot) // After IBZ reduction, weights should still reflect the input weights @@ -933,8 +939,8 @@ TEST_F(KlistTest, IbzKpointCustomWeights) { K_Vectors kv_test2; std::string k_file = "./support/KPT_custom_weights"; - kv_test2.nspin = 1; - kv_test2.read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv_test2.spin_mult = 1; + kv_test2.read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv_test2.get_nkstot(), 5); EXPECT_FALSE(kv_test2.is_mp); // Should be non-MP @@ -957,7 +963,7 @@ TEST_F(KlistTest, IbzKpointCustomWeights) std::string skpt; ModuleSymmetry::Symmetry::symm_flag = 1; bool match = true; - kv_test2.reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match); + kv_test2.reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match, my_rank, ofs_running); // After IBZ reduction, the weights should be based on the custom input weights, // not uniform 1/nkstot weights. The total weight should be preserved. @@ -990,8 +996,8 @@ TEST_F(KlistTest, IbzKpointCustomWeights) { K_Vectors kv_test3; std::string k_file = "./support/KPT1"; - kv_test3.nspin = 1; - kv_test3.read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv_test3.spin_mult = 1; + kv_test3.read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); EXPECT_EQ(kv_test3.get_nkstot(), 512); EXPECT_TRUE(kv_test3.is_mp); // Should be MP @@ -999,7 +1005,7 @@ TEST_F(KlistTest, IbzKpointCustomWeights) std::string skpt; ModuleSymmetry::Symmetry::symm_flag = 1; bool match = true; - kv_test3.reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match); + kv_test3.reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match, my_rank, ofs_running); // For MP grids, all weights should be uniform after IBZ reduction EXPECT_EQ(kv_test3.get_nkstot(), 35); // Known result from existing test @@ -1017,17 +1023,17 @@ TEST_F(KlistTest, IbzKpointCustomWeights) { K_Vectors kv_test4; std::string k_file = "./support/KPT_custom_weights"; - kv_test4.nspin = 1; - kv_test4.read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv_test4.spin_mult = 1; + kv_test4.read_kpoints(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank); // Apply IBZ reduction std::string skpt; ModuleSymmetry::Symmetry::symm_flag = 1; bool match = true; - kv_test4.reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match); + kv_test4.reduce_by_symmetry(ucell, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, match, my_rank, ofs_running); // Normalize weights - int degspin = (kv_test4.nspin == 2) ? 1 : 2; + int degspin = (kv_test4.spin_mult == 2) ? 1 : 2; kv_test4.normalize_wk(degspin); // After normalization, weights should sum to degspin @@ -1039,7 +1045,7 @@ TEST_F(KlistTest, IbzKpointCustomWeights) EXPECT_NEAR(total_weight, degspin, 1e-10); } - GlobalV::ofs_running.close(); + ofs_running.close(); ClearUcell(); remove("tmp_klist_custom_weights"); } diff --git a/source/source_cell/test/klist_test_para.cpp b/source/source_cell/test/klist_test_para.cpp index 17a39ed576e..2ebc0c609b5 100644 --- a/source/source_cell/test/klist_test_para.cpp +++ b/source/source_cell/test/klist_test_para.cpp @@ -183,8 +183,8 @@ TEST_F(KlistParaTest, Set) symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); // read KPT std::string k_file = "./support/KPT1"; - // set klist - kv->nspin = 1; + // note: do NOT pre-set kv->spin_mult here; set() takes the physical + // nspin as input and performs the 4->1 mapping internally. if (GlobalV::NPROC == 4) { GlobalV::KPAR = 2; @@ -207,10 +207,10 @@ TEST_F(KlistParaTest, Set) const double kspacing[3] = {0.0, 0.0, 0.0}; const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; - kv->set(ucell, symm, k_file, kv->nspin, ucell.G, ucell.latvec, GlobalV::ofs_running, use_ibz, global_out_dir, gamma_only_local, kspacing, kmesh_type, koffset); + kv->set(ucell, symm, k_file, /*nspin_in*/ 1, ucell.G, ucell.latvec, GlobalV::ofs_running, GlobalV::ofs_warning, use_ibz, global_out_dir, gamma_only_local, kspacing, kmesh_type, koffset); EXPECT_EQ(kv->get_nkstot(), 35); - EXPECT_EQ(kv->get_nkstot_full(), 512); - EXPECT_GT(kv->get_nkstot_full(), kv->get_nkstot()); + EXPECT_EQ(kv->get_nkstot_nospin(), 512); + EXPECT_GT(kv->get_nkstot_nospin(), kv->get_nkstot()); EXPECT_TRUE(kv->kc_done); EXPECT_TRUE(kv->kd_done); if (GlobalV::NPROC == 4) @@ -307,8 +307,8 @@ TEST_F(KlistParaTest, SetAfterVC) symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running, 1e-6, 1, "scf", cal_symm_repr); // read KPT std::string k_file = "./support/KPT1"; - // set klist - kv->nspin = 1; + // note: do NOT pre-set kv->spin_mult here; set() takes the physical + // nspin as input and performs the 4->1 mapping internally. if (GlobalV::NPROC == 4) { GlobalV::KPAR = 1; @@ -331,7 +331,7 @@ TEST_F(KlistParaTest, SetAfterVC) const double kspacing[3] = {0.0, 0.0, 0.0}; const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; - kv->set(ucell, symm, k_file, kv->nspin, ucell.G, ucell.latvec, GlobalV::ofs_running, use_ibz, global_out_dir, gamma_only_local, kspacing, kmesh_type, koffset); + kv->set(ucell, symm, k_file, /*nspin_in*/ 1, ucell.G, ucell.latvec, GlobalV::ofs_running, GlobalV::ofs_warning, use_ibz, global_out_dir, gamma_only_local, kspacing, kmesh_type, koffset); EXPECT_EQ(kv->get_nkstot(), 35); EXPECT_TRUE(kv->kc_done); EXPECT_TRUE(kv->kd_done); @@ -352,7 +352,7 @@ TEST_F(KlistParaTest, SetAfterVC) } // call set_after_vc here kv->kc_done = false; - kv->set_after_vc(kv->nspin, ucell.G, GlobalV::ofs_running); + kv->set_after_vc(ucell.G, GlobalV::ofs_running); EXPECT_TRUE(kv->kc_done); EXPECT_TRUE(kv->kd_done); // clear diff --git a/source/source_cell/test/qlist_test.cpp b/source/source_cell/test/qlist_test.cpp index 8e828356a2e..c9c07e6dffb 100644 --- a/source/source_cell/test/qlist_test.cpp +++ b/source/source_cell/test/qlist_test.cpp @@ -180,7 +180,7 @@ TEST_F(QListTest, GenerateMeshFullSymmetry) qlist.generate_mesh(ucell, symm, {8, 8, 8}, true); // full mesh 512 -> irreducible q-points of the primitive cubic lattice - EXPECT_EQ(qlist.nkstot_full, 512); + EXPECT_EQ(qlist.nkstot_nospin, 512); EXPECT_EQ(qlist.get_nq(), 35); EXPECT_EQ(qlist.get_nq(), qlist.nkstot); EXPECT_TRUE(qlist.is_mp); @@ -218,7 +218,7 @@ TEST_F(QListTest, GenerateMeshSmallGrid) qlist.generate_mesh(ucell, symm, {2, 2, 2}, true); // {0,0.5}^3 under O_h folds to Gamma + X + M + R - EXPECT_EQ(qlist.nkstot_full, 8); + EXPECT_EQ(qlist.nkstot_nospin, 8); EXPECT_EQ(qlist.get_nq(), 4); // the first irreducible q-point must be Gamma (0,0,0) @@ -241,7 +241,7 @@ TEST_F(QListTest, GammaOnlyGrid) qlist.generate_mesh(ucell, symm, {1, 1, 1}, true); - EXPECT_EQ(qlist.nkstot_full, 1); + EXPECT_EQ(qlist.nkstot_nospin, 1); EXPECT_EQ(qlist.get_nq(), 1); EXPECT_DOUBLE_EQ(qlist.wk[0], 1.0); EXPECT_DOUBLE_EQ(qlist.get_q(0).x, 0.0); diff --git a/source/source_cell/test/reciprocal_grid_test.cpp b/source/source_cell/test/reciprocal_grid_test.cpp index c369bf5ff4c..05525749d38 100644 --- a/source/source_cell/test/reciprocal_grid_test.cpp +++ b/source/source_cell/test/reciprocal_grid_test.cpp @@ -66,7 +66,9 @@ class TestGrid : public ModuleCell::ReciprocalGrid const ModuleSymmetry::Symmetry&, bool, std::string&, - bool&) override + bool&, + const int, + std::ofstream&) override { } }; @@ -81,7 +83,7 @@ TEST_F(ReciprocalGridTest, Construct) { EXPECT_EQ(grid.nks, 0); EXPECT_EQ(grid.nkstot, 0); - EXPECT_EQ(grid.nkstot_full, 0); + EXPECT_EQ(grid.nkstot_nospin, 0); EXPECT_FALSE(grid.kc_done); EXPECT_FALSE(grid.kd_done); EXPECT_FALSE(grid.is_mp); @@ -214,7 +216,7 @@ TEST_F(ReciprocalGridTest, ReduceIbzNonMp) grid.is_mp = false; grid.nkstot = 2; - grid.nkstot_full = 2; + grid.nkstot_nospin = 2; grid.kvec_d.resize(2); grid.kvec_d[0] = ModuleBase::Vector3(0.25, 0.25, 0.25); grid.kvec_d[1] = ModuleBase::Vector3(-0.25, -0.25, -0.25); @@ -244,7 +246,7 @@ TEST_F(ReciprocalGridTest, ReduceIbzKeepsDistinctPoints) grid.is_mp = false; grid.nkstot = 2; - grid.nkstot_full = 2; + grid.nkstot_nospin = 2; grid.kvec_d.resize(2); grid.kvec_d[0] = ModuleBase::Vector3(0.25, 0.25, 0.25); grid.kvec_d[1] = ModuleBase::Vector3(0.50, 0.50, 0.50); @@ -281,7 +283,7 @@ TEST_F(ReciprocalGridTest, ReduceIbzMpKLattice) const double offset[3] = {0.0, 0.0, 0.0}; grid.Monkhorst_Pack(nmp, offset, 0); // sets nkstot=8, wk=1/8, kd_done grid.is_mp = true; - grid.nkstot_full = grid.nkstot; + grid.nkstot_nospin = grid.nkstot; // k-lattice basis of the 2x2x2 mesh: G/2 along each reciprocal axis. // In this diagonal frame the k-lattice rotations equal the reciprocal ones. diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index b0114402c67..779f62f2926 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -97,7 +97,7 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp) const bool gamma_only_local = PARAM.globalv.gamma_only_local; const double kspacing[3] = {this->inp_->kspacing[0], this->inp_->kspacing[1], this->inp_->kspacing[2]}; const double koffset[3] = {this->inp_->koffset[0], this->inp_->koffset[1], this->inp_->koffset[2]}; - this->kv.set(ucell, ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running, use_ibz, global_out_dir, gamma_only_local, kspacing, this->inp_->kmesh_type, koffset); + this->kv.set(ucell, ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running, GlobalV::ofs_warning, use_ibz, global_out_dir, gamma_only_local, kspacing, this->inp_->kmesh_type, koffset); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); //! 8) print information @@ -182,7 +182,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) } // reset k-points - kv.set_after_vc(this->inp_->nspin, ucell.G, GlobalV::ofs_running); + kv.set_after_vc(ucell.G, GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); } diff --git a/source/source_esolver/esolver_gets.cpp b/source/source_esolver/esolver_gets.cpp index b1cdceea24c..85649f960cc 100644 --- a/source/source_esolver/esolver_gets.cpp +++ b/source/source_esolver/esolver_gets.cpp @@ -91,6 +91,7 @@ void ESolver_GetS::before_all_runners(BaseCell& basecell, const Input_para& inp) ucell.G, ucell.latvec, GlobalV::ofs_running, + GlobalV::ofs_warning, use_ibz, global_out_dir, gamma_only_local, diff --git a/source/source_esolver/esolver_lr_lcao_tddft.cpp b/source/source_esolver/esolver_lr_lcao_tddft.cpp index 52bd6ad0e89..6953a2d8dfd 100644 --- a/source/source_esolver/esolver_lr_lcao_tddft.cpp +++ b/source/source_esolver/esolver_lr_lcao_tddft.cpp @@ -373,7 +373,7 @@ void ModuleESolver::ESolver_LR::initialize_from_unitcell_(UnitCell& ucell const bool gamma_only_local = PARAM.globalv.gamma_only_local; const double kspacing[3] = {this->inp_->kspacing[0], this->inp_->kspacing[1], this->inp_->kspacing[2]}; const double koffset[3] = {this->inp_->koffset[0], this->inp_->koffset[1], this->inp_->koffset[2]}; - this->kv.set(ucell, ucell.symm, this->inp_->kpoint_file, this->inp_->nspin, ucell.G, ucell.latvec, GlobalV::ofs_running, use_ibz, this->out_dir, gamma_only_local, kspacing, this->inp_->kmesh_type, koffset); + this->kv.set(ucell, ucell.symm, this->inp_->kpoint_file, this->inp_->nspin, ucell.G, ucell.latvec, GlobalV::ofs_running, GlobalV::ofs_warning, use_ibz, this->out_dir, gamma_only_local, kspacing, this->inp_->kmesh_type, koffset); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); ModuleIO::print_parameters(ucell, this->kv, inp); diff --git a/source/source_estate/module_dm/test/CMakeLists.txt b/source/source_estate/module_dm/test/CMakeLists.txt index 1d96186fbec..d5e8c19c3a5 100644 --- a/source/source_estate/module_dm/test/CMakeLists.txt +++ b/source/source_estate/module_dm/test/CMakeLists.txt @@ -16,7 +16,7 @@ AddTest( ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/hcontainer.cpp ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/atom_pair.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp - ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp + ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp ${ABACUS_SOURCE_DIR}/source_cell/klist_io.cpp ${ABACUS_SOURCE_DIR}/source_cell/parallel_kpoints.cpp ${ABACUS_SOURCE_DIR}/source_cell/reciprocal_grid.cpp @@ -30,7 +30,7 @@ AddTest( ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/hcontainer.cpp ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/atom_pair.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp - ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp + ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp ${ABACUS_SOURCE_DIR}/source_cell/klist_io.cpp ${ABACUS_SOURCE_DIR}/source_cell/parallel_kpoints.cpp ${ABACUS_SOURCE_DIR}/source_cell/reciprocal_grid.cpp ) @@ -43,7 +43,7 @@ AddTest( ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/hcontainer.cpp ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/atom_pair.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp - ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp + ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp ${ABACUS_SOURCE_DIR}/source_cell/klist_io.cpp ${ABACUS_SOURCE_DIR}/source_cell/parallel_kpoints.cpp ${ABACUS_SOURCE_DIR}/source_cell/reciprocal_grid.cpp ) @@ -56,7 +56,7 @@ AddTest( ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/hcontainer.cpp ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/atom_pair.cpp ${ABACUS_SOURCE_DIR}/source_basis/module_ao/parallel_orbitals.cpp - ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp + ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp ${ABACUS_SOURCE_DIR}/source_cell/klist_io.cpp ${ABACUS_SOURCE_DIR}/source_cell/parallel_kpoints.cpp ${ABACUS_SOURCE_DIR}/source_cell/reciprocal_grid.cpp ) diff --git a/source/source_estate/test/CMakeLists.txt b/source/source_estate/test/CMakeLists.txt index 1eaaa0c6bce..96c84ceeace 100644 --- a/source/source_estate/test/CMakeLists.txt +++ b/source/source_estate/test/CMakeLists.txt @@ -33,7 +33,7 @@ AddTest( TARGET MODULE_ESTATE_elecstate_print LIBS parameter base device symmetry SOURCES elecstate_print_test.cpp ../elecstate_print.cpp ../occupy.cpp - ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) AddTest( @@ -41,7 +41,7 @@ AddTest( LIBS parameter base device symmetry SOURCES elecstate_base_test.cpp ../elecstate.cpp ../elecstate_tools.cpp ../occupy.cpp ../../source_psi/psi.cpp ../../source_base/module_fft/fft_bundle.cpp ../../source_base/module_fft/fft_cpu.cpp - ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) AddTest( @@ -57,7 +57,7 @@ AddTest( ../../source_lcao/module_deltaspin/deltaspin_state.cpp ../../source_psi/psi.cpp ../../source_base/module_device/memory_op.cpp - ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) AddTest( @@ -68,7 +68,7 @@ AddTest( ../fp_energy.cpp ../makov_payne.cpp ../module_pot/h_hartree_pw.cpp - ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) AddTest( diff --git a/source/source_hsolver/test/CMakeLists.txt b/source/source_hsolver/test/CMakeLists.txt index 72a2de07b12..96895a504ed 100644 --- a/source/source_hsolver/test/CMakeLists.txt +++ b/source/source_hsolver/test/CMakeLists.txt @@ -78,7 +78,7 @@ if (ENABLE_MPI) LIBS parameter psi device base container symmetry SOURCES test_hsolver_pw.cpp ../hsolver_pw.cpp ../hsolver_lcaopw.cpp ../diago_bpcg.cpp ../diago_dav_subspace.cpp ../diag_const_nums.cpp ../diago_iter_assist.cpp ../para_lin_tf.cpp ../../source_estate/elecstate_tools.cpp ../../source_estate/occupy.cpp ../../source_base/module_fft/fft_bundle.cpp ../../source_base/module_fft/fft_cpu.cpp - ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) AddTest( diff --git a/source/source_io/module_dm/test/write_dmk_test.cpp b/source/source_io/module_dm/test/write_dmk_test.cpp index 8041bde1489..e03f0530553 100644 --- a/source/source_io/module_dm/test/write_dmk_test.cpp +++ b/source/source_io/module_dm/test/write_dmk_test.cpp @@ -1,6 +1,20 @@ -#include "source_io/module_dm/write_dmk.h" +// Pre-include every standard-library header reachable from write_dmk.h so +// their include guards are already set before '#define private public' is +// active. The macro renames the 'private'/'public' keywords, so any system +// header parsed while it is defined gets corrupted and the build fails with +// "'...__xfer_bufptrs' redeclared with different access". write_dmk.h pulls +// in indirectly via global_variable.h -> -> +// bits/quoted_string.h, so must be pre-included too. +#include +#include +#include +#include +#include +#include +#include #define private public +#include "source_io/module_dm/write_dmk.h" #include "source_io/module_parameter/parameter.h" #undef private #include "source_base/global_variable.h" @@ -142,9 +156,9 @@ TEST(DMKTest,WriteDMK) { const int istep = -1; K_Vectors kv; kv.set_nkstot(1); - kv.set_nkstot_full(1); + kv.set_nkstot_nospin(1); kv.set_nks(1); - kv.set_nspin(2); + kv.spin_mult = 2; kv.kvec_c.resize(1); kv.kvec_c[0].x = 0.0; kv.kvec_c[0].y = 0.0; diff --git a/source/source_io/module_dm/write_dmk.cpp b/source/source_io/module_dm/write_dmk.cpp index ac58816fb9b..1d3a40b32b0 100644 --- a/source/source_io/module_dm/write_dmk.cpp +++ b/source/source_io/module_dm/write_dmk.cpp @@ -257,7 +257,7 @@ void ModuleIO::write_dmk(const std::vector>& dmk, // information about density matrix at this k-point ofs << " " << nspin << " # number of spin directions" << std::endl; ofs << " " << ispin+1 << " # spin index" << std::endl; - ofs << " " << kv.get_nkstot_full() << " # total k points " << std::endl; + ofs << " " << kv.get_nkstot_nospin() << " # total k points " << std::endl; ofs << " " << kv.get_nkstot() << " # total k points after symmetrized (if open) " << std::endl; ofs << " " << ik+1 << " # k-point index " << std::endl; ofs << " " << kv.kvec_c[ik].x << " " << kv.kvec_c[ik].y << " " << kv.kvec_c[ik].z diff --git a/source/source_io/test/CMakeLists.txt b/source/source_io/test/CMakeLists.txt index 4385b04dabc..b93bd463be2 100644 --- a/source/source_io/test/CMakeLists.txt +++ b/source/source_io/test/CMakeLists.txt @@ -63,7 +63,7 @@ add_test( AddTest( TARGET MODULE_IO_write_eig_occ_test LIBS parameter base device symmetry - SOURCES write_eig_occ_test.cpp ../module_output/band_parallel_output.cpp ../module_energy/write_eig_occ.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/klist.cpp + SOURCES write_eig_occ_test.cpp ../module_output/band_parallel_output.cpp ../module_energy/write_eig_occ.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/cif_io.cpp ../../source_cell/reciprocal_grid.cpp ) @@ -76,13 +76,13 @@ AddTest( AddTest( TARGET MODULE_IO_write_dos_pw LIBS parameter base device symmetry - SOURCES write_dos_pw_test.cpp ../module_dos/cal_dos.cpp ../module_dos/write_dos_pw.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/klist.cpp ../module_energy/nscf_fermi_surf.cpp ../../source_cell/reciprocal_grid.cpp + SOURCES write_dos_pw_test.cpp ../module_dos/cal_dos.cpp ../module_dos/write_dos_pw.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../module_energy/nscf_fermi_surf.cpp ../../source_cell/reciprocal_grid.cpp ) AddTest( TARGET MODULE_IO_print_info LIBS parameter base device symmetry cell_info - SOURCES print_info_test.cpp ../module_output/print_info.cpp ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + SOURCES print_info_test.cpp ../module_output/print_info.cpp ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) AddTest( @@ -177,7 +177,7 @@ AddTest( TARGET MODULE_IO_read_wf2rho_pw_test LIBS parameter base device planewave psi symmetry SOURCES read_wf2rho_pw_test.cpp ../module_wf/read_wfc_pw.cpp ../module_wf/read_wf2rho_pw.cpp ../../source_basis/module_pw/test/test_tool.cpp ../../source_estate/module_charge/charge_mpi.cpp ../module_wf/write_wfc_pw.cpp - ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) add_test(NAME MODULE_IO_read_wf2rho_pw_parallel @@ -232,7 +232,7 @@ AddTest( TARGET MODULE_IO_write_dmk LIBS parameter base device cell_info symmetry SOURCES ../module_dm/test/write_dmk_test.cpp ../module_dm/write_dmk.cpp ../../source_cell/ucell_io.cpp - ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) add_test( diff --git a/source/source_io/test/print_info_test.cpp b/source/source_io/test/print_info_test.cpp index 44d500b98fd..62853c34f3d 100644 --- a/source/source_io/test/print_info_test.cpp +++ b/source/source_io/test/print_info_test.cpp @@ -48,12 +48,12 @@ TEST_F(PrintInfoTest, SetupParameters) UcellTestPrepare utp = UcellTestLib["Si"]; ucell = utp.SetUcellInfo(); std::string k_file = "./support/KPT"; - kv->nspin = 1; + kv->spin_mult = 1; const bool gamma_only_local = false; const double kspacing[3] = {0.0, 0.0, 0.0}; const std::string kmesh_type = "gamma"; const double koffset[3] = {0.0, 0.0, 0.0}; - kv->read_kpoints(*ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running); + kv->read_kpoints(*ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, GlobalV::ofs_running, GlobalV::ofs_warning, GlobalV::MY_RANK); EXPECT_EQ(kv->get_nkstot(),512); std::vector cal_type = {"scf","relax","cell-relax","md"}; std::vector md_types = {"fire","nve","nvt","npt","langevin","msst"}; diff --git a/source/source_io/test_serial/CMakeLists.txt b/source/source_io/test_serial/CMakeLists.txt index aa7fec69f4c..dd87df580c5 100644 --- a/source/source_io/test_serial/CMakeLists.txt +++ b/source/source_io/test_serial/CMakeLists.txt @@ -58,7 +58,7 @@ AddTest( TARGET MODULE_IO_write_bands LIBS parameter base device symmetry SOURCES write_bands_test.cpp ../module_output/band_parallel_output.cpp ../module_energy/write_bands.cpp - ../../source_cell/klist.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp + ../../source_cell/klist.cpp ../../source_cell/klist_io.cpp ../../source_cell/parallel_kpoints.cpp ../../source_cell/reciprocal_grid.cpp ) AddTest( diff --git a/source/source_lcao/module_bse/molecular_lri.hpp b/source/source_lcao/module_bse/molecular_lri.hpp index 3f692d76db0..a2efe734bb7 100644 --- a/source/source_lcao/module_bse/molecular_lri.hpp +++ b/source/source_lcao/module_bse/molecular_lri.hpp @@ -160,7 +160,7 @@ void MolecularLRI::build_q_to_kpair_map(int mode, double threshold) { std::set q_coarse_set; const K_Vectors& kv_coarse = this->kRlist.klist_coarse; - int nk_coarse = kv_coarse.get_nkstot_full(); + int nk_coarse = kv_coarse.get_nkstot_nospin(); for (int ik1 = 0; ik1 < nk_coarse; ++ik1) { Tk ck1 = RI_Util::Vector3_to_array3(kv_coarse.kvec_d.at(ik1)); diff --git a/source/source_lcao/module_deepks/test/CMakeLists.txt b/source/source_lcao/module_deepks/test/CMakeLists.txt index df699495f90..15cca5c62d6 100644 --- a/source/source_lcao/module_deepks/test/CMakeLists.txt +++ b/source/source_lcao/module_deepks/test/CMakeLists.txt @@ -25,7 +25,7 @@ set(DEEPKS_UNIT_COMMON_SOURCES ../../../source_cell/read_stru.cpp ../../../source_cell/print_cell.cpp ../../../source_cell/read_atom_species.cpp - ../../../source_cell/klist.cpp + ../../../source_cell/klist.cpp ../../../source_cell/klist_io.cpp ../../../source_cell/parallel_kpoints.cpp ../../../source_cell/reciprocal_grid.cpp ../../setup_nonlocal.cpp diff --git a/source/source_lcao/module_deepks/test/deepks_test_prep.cpp b/source/source_lcao/module_deepks/test/deepks_test_prep.cpp index bcb0f9c25e2..8eee3722856 100644 --- a/source/source_lcao/module_deepks/test/deepks_test_prep.cpp +++ b/source/source_lcao/module_deepks/test/deepks_test_prep.cpp @@ -319,6 +319,7 @@ void test_deepks::setup_kpt() ucell.G, ucell.latvec, GlobalV::ofs_running, + GlobalV::ofs_warning, use_ibz, global_out_dir, gamma_only_local, diff --git a/source/source_lcao/module_deltaspin/test/CMakeLists.txt b/source/source_lcao/module_deltaspin/test/CMakeLists.txt index a3bac8f9ee5..9fad70e2bc8 100644 --- a/source/source_lcao/module_deltaspin/test/CMakeLists.txt +++ b/source/source_lcao/module_deltaspin/test/CMakeLists.txt @@ -16,7 +16,7 @@ AddTest( ../spin_constrain.cpp ../deltaspin_state.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp - ../../../source_cell/klist.cpp + ../../../source_cell/klist.cpp ../../../source_cell/klist_io.cpp ../../../source_cell/parallel_kpoints.cpp ../../../source_cell/reciprocal_grid.cpp ) @@ -30,7 +30,7 @@ AddTest( ../template_helpers.cpp ../lambda_loop_helper.cpp ../basic_funcs.cpp - ../../../source_cell/klist.cpp + ../../../source_cell/klist.cpp ../../../source_cell/klist_io.cpp ../../../source_cell/parallel_kpoints.cpp ../../../source_cell/reciprocal_grid.cpp ) diff --git a/source/source_lcao/module_lr/utils/lr_io.cpp b/source/source_lcao/module_lr/utils/lr_io.cpp index 2c61478e7df..9d818c043cd 100644 --- a/source/source_lcao/module_lr/utils/lr_io.cpp +++ b/source/source_lcao/module_lr/utils/lr_io.cpp @@ -165,7 +165,7 @@ void RI_kRlist::read_kpts_fine(const std::string& file, const UnitCell& ucell, int nks = (PARAM.inp.nspin == 2) ? 2 * nk : nk; klist->set_nks(nks); klist->set_nkstot(nks); - klist->set_nkstot_full(nk); + klist->set_nkstot_nospin(nk); auto klist_reset = [&klist](int kpoint_number){ klist->kvec_c.resize(0); klist->kvec_c.resize(kpoint_number); diff --git a/source/source_lcao/module_rdmft/rdmft.cpp b/source/source_lcao/module_rdmft/rdmft.cpp index b88440a88f4..5a27d9d86f8 100644 --- a/source/source_lcao/module_rdmft/rdmft.cpp +++ b/source/source_lcao/module_rdmft/rdmft.cpp @@ -81,7 +81,7 @@ void RDMFT::init(Parallel_Orbitals& ParaV_in, nspin = PARAM.inp.nspin; nbands_total = PARAM.inp.nbands; - nk_total = ModuleSymmetry::Symmetry::symm_flag == -1 ? kv->get_nkstot_full(): kv->get_nks(); + nk_total = ModuleSymmetry::Symmetry::symm_flag == -1 ? kv->get_nkstot_nospin(): kv->get_nks(); nk_total *= nspin; only_exx_type = ( XC_func_rdmft == "hf" || XC_func_rdmft == "muller" || XC_func_rdmft == "power" ); diff --git a/source/source_lcao/module_ri/ewald_vq.hpp b/source/source_lcao/module_ri/ewald_vq.hpp index 584d8a28804..5e8e091b1a8 100644 --- a/source/source_lcao/module_ri/ewald_vq.hpp +++ b/source/source_lcao/module_ri/ewald_vq.hpp @@ -48,7 +48,7 @@ void Ewald_Vq::init(const UnitCell& ucell, this->mpi_comm = mpi_comm_in; this->p_kv = kv_in; - this->nks0 = this->p_kv->get_nkstot_full(); + this->nks0 = this->p_kv->get_nkstot_nospin(); this->kvec_c.resize(this->nks0); this->ccp_rmesh_times = ccp_rmesh_times_in; this->abfs_Lmax = abfs_Lmax_in; diff --git a/source/source_lcao/module_ri/exx_lri.hpp b/source/source_lcao/module_ri/exx_lri.hpp index 21c8afbb1cd..9cb8f75d720 100644 --- a/source/source_lcao/module_ri/exx_lri.hpp +++ b/source/source_lcao/module_ri/exx_lri.hpp @@ -364,7 +364,7 @@ void Exx_LRI::cal_exx_ions(const UnitCell& ucell, std::cout << "Coulomb: number of atom-pairs inside atomic overlap is " << flag << ". " << std::endl; if (this->info.coul_moment == true) { - double hf_Rcut = std::pow(0.75 * this->p_kv->get_nkstot_full() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0); + double hf_Rcut = std::pow(0.75 * this->p_kv->get_nkstot_nospin() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0); // To cal Cs, we still cal all Vs(R) in r space // moment_abfs->cal_VR(ucell, // this->abfs, @@ -532,7 +532,7 @@ void Exx_LRI::cal_ewald_coulomb( std::cout << "Coulomb: number of atom-pairs inside atomic overlap is " << flag << ". " << std::endl; if (this->info.coul_moment == true) { - double hf_Rcut = std::pow(0.75 * this->p_kv->get_nkstot_full() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0); + double hf_Rcut = std::pow(0.75 * this->p_kv->get_nkstot_nospin() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0); // To cal Cs, we still cal all Vs(R) in r space // moment_abfs->cal_VR(ucell, // this->abfs, diff --git a/source/source_lcao/module_ri/exx_lri_detail.cpp b/source/source_lcao/module_ri/exx_lri_detail.cpp index ff8046159b6..d459fc83790 100644 --- a/source/source_lcao/module_ri/exx_lri_detail.cpp +++ b/source/source_lcao/module_ri/exx_lri_detail.cpp @@ -16,7 +16,7 @@ namespace ExxLriDetail double default_spencer_rcut(const UnitCell& ucell, const K_Vectors& kv) { - return std::pow(0.75 * kv.get_nkstot_full() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0); + return std::pow(0.75 * kv.get_nkstot_nospin() * ucell.omega / (ModuleBase::PI), 1.0 / 3.0); } CoulombParam build_center2_cut_coulomb_param(const CoulombParam& coulomb_param, diff --git a/source/source_lcao/module_ri/exx_lri_interface.hpp b/source/source_lcao/module_ri/exx_lri_interface.hpp index 32320fab5f2..f00011ab833 100644 --- a/source/source_lcao/module_ri/exx_lri_interface.hpp +++ b/source/source_lcao/module_ri/exx_lri_interface.hpp @@ -156,7 +156,7 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, if(this->info_global.cal_exx) { if (this->exx_spacegroup_symmetry) - { this->mix_DMk_2D.set_nks(kv.get_nkstot_full() * (PARAM.inp.nspin == 2 ? 2 : 1)); } + { this->mix_DMk_2D.set_nks(kv.get_nkstot_nospin() * (PARAM.inp.nspin == 2 ? 2 : 1)); } else { this->mix_DMk_2D.set_nks(kv.get_nks()); } diff --git a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp index 4a88554624d..18eb826888c 100644 --- a/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp +++ b/source/source_lcao/module_ri/module_exx_symmetry/symm_rotation.cpp @@ -98,7 +98,7 @@ namespace ModuleSymmetry ModuleBase::timer::start("Symmetry_rotation", "restore_dm"); std::vector>> dm_k_full; int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1; - dm_k_full.reserve(kv.get_nkstot_full() * nspin0); //nkstot_full didn't doubled by spin + dm_k_full.reserve(kv.get_nkstot_nospin() * nspin0); //nkstot_nospin didn't doubled by spin int nk = kv.get_nkstot() / nspin0; // (nspin=4) Sigma_y = I (x) sigma_y for the time-reversal spin flip; k-independent, build once. diff --git a/source/source_lcao/module_ri/ri_2d_comm.hpp b/source/source_lcao/module_ri/ri_2d_comm.hpp index fcb33684ed1..64cdc8786e5 100644 --- a/source/source_lcao/module_ri/ri_2d_comm.hpp +++ b/source/source_lcao/module_ri/ri_2d_comm.hpp @@ -196,7 +196,7 @@ auto RI_2D_Comm::split_m2D_ktoR_k(const UnitCell& ucell, const Tdata_m frac = SPIN_multiple * RI::Global_Func::convert(std::exp( -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * (kv.kvec_c[ik] * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); - if (static_cast(std::round(SPIN_multiple * kv.wk[ik] * kv.get_nkstot_full())) == 2) + if (static_cast(std::round(SPIN_multiple * kv.wk[ik] * kv.get_nkstot_nospin())) == 2) { set_mR_2D(mk_2D * (frac * 0.5) + tensor_conj(mk_2D * (frac * 0.5))); } else { set_mR_2D(mk_2D * frac); } @@ -205,7 +205,7 @@ auto RI_2D_Comm::split_m2D_ktoR_k(const UnitCell& ucell, { // traverse kstar, ik means ik_ibz for (auto& isym_kvd : kv.kstars[ik % ik_list.size()]) { - RI::Tensor mk_2D = RI_Util::Vector_to_Tensor(*mks_2D[ik_full + is_k * kv.get_nkstot_full()], pv.get_col_size(), pv.get_row_size()); + RI::Tensor mk_2D = RI_Util::Vector_to_Tensor(*mks_2D[ik_full + is_k * kv.get_nkstot_nospin()], pv.get_col_size(), pv.get_row_size()); const Tdata_m frac = SPIN_multiple * RI::Global_Func::convert(std::exp( -ModuleBase::TWO_PI * ModuleBase::IMAG_UNIT * ((isym_kvd.second * ucell.G) * (RI_Util::array3_to_Vector3(cell) * ucell.latvec)))); diff --git a/source/source_lcao/module_ri/ri_util.hpp b/source/source_lcao/module_ri/ri_util.hpp index 6baf57c3d75..886c5f1b377 100644 --- a/source/source_lcao/module_ri/ri_util.hpp +++ b/source/source_lcao/module_ri/ri_util.hpp @@ -82,7 +82,7 @@ namespace RI_Util { // 4/3 * pi * Rcut^3 = V_{supercell} = V_{unitcell} * Nk const int nspin0 = (PARAM.inp.nspin==2) ? 2 : 1; - const double Rcut = std::pow(0.75 * p_kv->get_nkstot_full() * ucell.omega / (ModuleBase::PI), 1.0/3.0); + const double Rcut = std::pow(0.75 * p_kv->get_nkstot_nospin() * ucell.omega / (ModuleBase::PI), 1.0/3.0); param["Rcut"] = ModuleBase::GlobalFunc::TO_STRING(Rcut); } else if(param.at("singularity_correction") == "revised_spencer") diff --git a/source/source_lcao/module_ri/rpa_lri.hpp b/source/source_lcao/module_ri/rpa_lri.hpp index c2fa83ce8cc..a5bf037d1cd 100644 --- a/source/source_lcao/module_ri/rpa_lri.hpp +++ b/source/source_lcao/module_ri/rpa_lri.hpp @@ -122,7 +122,7 @@ void RPA_LRI::cal_postSCF_exx(const elecstate::DensityMatrix Mix_DMk_2D mix_DMk_2D; bool exx_spacegroup_symmetry = (PARAM.inp.nspin < 4 && ModuleSymmetry::Symmetry::symm_flag == 1); if (exx_spacegroup_symmetry) - {mix_DMk_2D.set_nks(kv.get_nkstot_full() * (PARAM.inp.nspin == 2 ? 2 : 1));} + {mix_DMk_2D.set_nks(kv.get_nkstot_nospin() * (PARAM.inp.nspin == 2 ? 2 : 1));} else {mix_DMk_2D.set_nks(kv.get_nks());} diff --git a/source/source_lcao/test/CMakeLists.txt b/source/source_lcao/test/CMakeLists.txt index 1971ea0b056..8f5d6a1aee6 100644 --- a/source/source_lcao/test/CMakeLists.txt +++ b/source/source_lcao/test/CMakeLists.txt @@ -21,7 +21,7 @@ AddTest( ${ABACUS_SOURCE_DIR}/source_io/module_dm/write_dmr.cpp ${ABACUS_SOURCE_DIR}/source_cell/ucell_io.cpp ${ABACUS_SOURCE_DIR}/source_hamilt/module_hcontainer/output_hcontainer.cpp - ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp + ${ABACUS_SOURCE_DIR}/source_cell/klist.cpp ${ABACUS_SOURCE_DIR}/source_cell/klist_io.cpp ${ABACUS_SOURCE_DIR}/source_cell/parallel_kpoints.cpp ${ABACUS_SOURCE_DIR}/source_cell/reciprocal_grid.cpp ) diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index 546577e65f4..a249a9bc3bf 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -48,8 +48,7 @@ void Plus_U_Base::init_base(UnitCell& cell, ModuleBase::TITLE("Plus_U_Base", "init_base"); #ifndef __MPI - std::cout << "DFT+U module is only accessible in mpi version" << std::endl; - exit(0); + ModuleBase::WARNING_QUIT("Plus_U_Base::init_base", "DFT+U module is only accessible in MPI version"); #endif this->nspin = nspin; @@ -82,13 +81,14 @@ void Plus_U_Base::init_base(UnitCell& cell, { const int iat = cell.itia2iat(it, ia); - if(!has_correlated_orbital(it)) + const int target_l = this->orbital_corr[it]; + if (target_l == -1) { continue; } - const int tlp1_npol = (get_orbital_corr(it)*2+1)*npol; - const int tlp1 = 2 * get_orbital_corr(it) + 1; + const int tlp1_npol = (target_l * 2 + 1) * npol; + const int tlp1 = 2 * target_l + 1; const int elem_size = tlp1 * tlp1; if(nspin == 4) { @@ -148,7 +148,7 @@ void Plus_U_Base::init_base(UnitCell& cell, DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->orbital_corr, nspin, npol); #endif - mark_occ_mat_initialized(); + this->occ_mat_initialized = true; this->occmat_.copy_to_save(cell, this->orbital_corr); this->occmat_.write_save_to_flat(cell, this->orbital_corr, this->pot_uterm_pw_index, this->uom_save); @@ -164,7 +164,7 @@ void Plus_U_Base::init_base(UnitCell& cell, #ifdef __MPI DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->orbital_corr, nspin, npol); #endif - mark_occ_mat_initialized(); + this->occ_mat_initialized = true; } else { @@ -173,18 +173,19 @@ void Plus_U_Base::init_base(UnitCell& cell, } ModuleBase::Memory::record("Plus_U_Base::occ_mat", sizeof(double) * num_locale); - return; } void Plus_U_Base::uramping_update() { // Yukawa calculates U directly every iteration, no need for ramping - if (use_yukawa()) { + if (this->yukawa_ != nullptr) + { return; } // if uramping < 0.1, use the original U - if (this->uramping < 0.01) { + if (this->uramping < 0.01) + { return; } // loop to change U @@ -205,7 +206,8 @@ void Plus_U_Base::uramping_update() bool Plus_U_Base::u_converged() { // Yukawa calculates U directly every iteration, always considered converged - if (use_yukawa()) { + if (this->yukawa_ != nullptr) + { return true; } for (int i = 0; i < static_cast(this->u_target.size()); i++) diff --git a/source/source_pw/module_pwdft/dftu_base.h b/source/source_pw/module_pwdft/dftu_base.h index f91ecf6ad3a..25671655fb4 100644 --- a/source/source_pw/module_pwdft/dftu_base.h +++ b/source/source_pw/module_pwdft/dftu_base.h @@ -113,7 +113,6 @@ class Plus_U_Base void mark_occ_mat_initialized() { occ_mat_initialized = true; } void mark_occ_mat_dirty() { occ_mat_initialized = false; } - bool is_mixing_enabled() const { return mixing_dftu != 0; } void enable_mixing() { mixing_dftu = 1; } /// direct access to the occupation matrix object (new write path) diff --git a/source/source_pw/module_pwdft/dftu_base_occ.cpp b/source/source_pw/module_pwdft/dftu_base_occ.cpp index 42b40e96a43..bb9e63b09ca 100644 --- a/source/source_pw/module_pwdft/dftu_base_occ.cpp +++ b/source/source_pw/module_pwdft/dftu_base_occ.cpp @@ -57,7 +57,7 @@ void Plus_U_Base::cal_occ_pw(const void* psi_in, this->pot_uterm_pw_index, this->uom_array); // mixing - if(is_mixing_enabled() && p_chgmix != nullptr) + if (this->mixing_dftu != 0 && p_chgmix != nullptr) { p_chgmix->mix_uom(this->uom_array, this->uom_save); this->occmat_.read_from_flat(cell, this->orbital_corr, diff --git a/source/source_pw/module_pwdft/op_pw_exx.cpp b/source/source_pw/module_pwdft/op_pw_exx.cpp index cafbc997ee1..94f19ce3dec 100644 --- a/source/source_pw/module_pwdft/op_pw_exx.cpp +++ b/source/source_pw/module_pwdft/op_pw_exx.cpp @@ -314,7 +314,7 @@ void OperatorEXXPW::act_op_kpar(const int nbands, // std::map, bool> has_real; setmem_complex_op()(psi_nk_real, 0, wfcpw->nrxx); setmem_complex_op()(psi_mq_real, 0, wfcpw->nrxx); - int nqs = kv->get_nkstot_full(); + int nqs = kv->get_nkstot_nospin(); int nspin_fac = PARAM.inp.nspin == 2 ? 2 : 1; int ispin = this->ik < (wfcpw->nks / nspin_fac) ? 0 : 1; diff --git a/source/source_pw/module_pwdft/op_pw_exx_ace.cpp b/source/source_pw/module_pwdft/op_pw_exx_ace.cpp index 55917e5b7c8..67322bdc493 100644 --- a/source/source_pw/module_pwdft/op_pw_exx_ace.cpp +++ b/source/source_pw/module_pwdft/op_pw_exx_ace.cpp @@ -140,7 +140,7 @@ void OperatorEXXPW::construct_ace() const setmem_complex_op()(density_recip, 0, rhopw_dev->npw); setmem_complex_op()(psi_nk_real, 0, wfcpw->nrxx); setmem_complex_op()(psi_mq_real, 0, wfcpw->nrxx); - int nqs = kv->get_nkstot_full(); + int nqs = kv->get_nkstot_nospin(); bool skip_ik = false; if (ik >= wfcpw->nks) diff --git a/source/source_pw/module_pwdft/op_pw_exx_pot.cpp b/source/source_pw/module_pwdft/op_pw_exx_pot.cpp index 216617eaa9a..fd76372e317 100644 --- a/source/source_pw/module_pwdft/op_pw_exx_pot.cpp +++ b/source/source_pw/module_pwdft/op_pw_exx_pot.cpp @@ -505,7 +505,7 @@ double exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, } } - div *= ModuleBase::e2 * ModuleBase::FOUR_PI / tpiba2 / kv->get_nkstot_full(); + div *= ModuleBase::e2 * ModuleBase::FOUR_PI / tpiba2 / kv->get_nkstot_nospin(); // std::cout << "div: " << div << std::endl; // numerically value the mean value of F(q) in the reciprocal space @@ -532,7 +532,7 @@ double exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, aa += 1.0 / std::sqrt(alpha * ModuleBase::PI); div -= ModuleBase::e2 * ucell_omega * aa; - exx_div = div * kv->get_nkstot_full(); + exx_div = div * kv->get_nkstot_nospin(); // exx_div = 0; // std::cout << "EXX divergence: " << exx_div << std::endl; diff --git a/source/source_pw/module_stodft/test/CMakeLists.txt b/source/source_pw/module_stodft/test/CMakeLists.txt index c025b789558..836a1ac7dcd 100644 --- a/source/source_pw/module_stodft/test/CMakeLists.txt +++ b/source/source_pw/module_stodft/test/CMakeLists.txt @@ -10,5 +10,5 @@ AddTest( TARGET MODULE_PW_Sto_Hamilt_UTs LIBS parameter psi base device planewave_serial symmetry SOURCES ../hamilt_sdft_pw.cpp test_hamilt_sto.cpp ../../../source_hamilt/operator.cpp - ../../../source_cell/klist.cpp ../../../source_cell/parallel_kpoints.cpp ../../../source_cell/reciprocal_grid.cpp + ../../../source_cell/klist.cpp ../../../source_cell/klist_io.cpp ../../../source_cell/parallel_kpoints.cpp ../../../source_cell/reciprocal_grid.cpp ) \ No newline at end of file