From dc6c811cc5e656136f90d9e4a2ea1c91e50ea98b Mon Sep 17 00:00:00 2001 From: Haerxile Date: Wed, 29 May 2024 10:56:38 +0800 Subject: [PATCH] Add notes in doxygen and modify the code indentation format --- source/module_cell/klist.cpp | 986 +++++++++++++++++++---------------- source/module_cell/klist.h | 371 +++++++++++-- 2 files changed, 861 insertions(+), 496 deletions(-) diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp index 9f1b40e6f1..e28814a329 100644 --- a/source/module_cell/klist.cpp +++ b/source/module_cell/klist.cpp @@ -1,12 +1,13 @@ -#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "klist.h" + +#include "module_base/formatter.h" +#include "module_base/memory.h" +#include "module_base/parallel_common.h" #include "module_base/parallel_global.h" -#include "module_cell/module_symmetry/symmetry.h" #include "module_base/parallel_reduce.h" -#include "module_base/parallel_common.h" -#include "module_base/memory.h" +#include "module_cell/module_symmetry/symmetry.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/berryphase.h" -#include "module_base/formatter.h" #ifdef USE_PAW #include "module_cell/module_paw/paw_cell.h" #endif @@ -15,14 +16,14 @@ K_Vectors::K_Vectors() { #ifdef _MCD_CHECK FILE* out; - out=fopen("1_Memory", "w"); - if ( out == NULL ) + out = fopen("1_Memory", "w"); + if (out == NULL) { std::cout << "\n Can't open file!"; ModuleBase::QUIT(); } - _MCD_RealTimeLog( out ); - _MCD_MemStatLog( out ); + _MCD_RealTimeLog(out); + _MCD_MemStatLog(out); // showMemStats(); #endif @@ -34,7 +35,7 @@ K_Vectors::K_Vectors() nkstot = 0; nkstot_ibz = 0; - k_nkstot = 0; //LiuXh add 20180619 + k_nkstot = 0; // LiuXh add 20180619 } K_Vectors::~K_Vectors() @@ -45,62 +46,62 @@ K_Vectors::~K_Vectors() #endif } -void K_Vectors::set( - const ModuleSymmetry::Symmetry &symm, - const std::string &k_file_name, - const int& nspin_in, - const ModuleBase::Matrix3 &reciprocal_vec, - const ModuleBase::Matrix3 &latvec) +void K_Vectors::set(const ModuleSymmetry::Symmetry& symm, + const std::string& k_file_name, + const int& nspin_in, + const ModuleBase::Matrix3& reciprocal_vec, + const ModuleBase::Matrix3& latvec) { ModuleBase::TITLE("K_Vectors", "set"); - GlobalV::ofs_running << "\n\n\n\n"; - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - GlobalV::ofs_running << " | |" << std::endl; - GlobalV::ofs_running << " | Setup K-points |" << std::endl; - GlobalV::ofs_running << " | We setup the k-points according to input parameters. |" << std::endl; - GlobalV::ofs_running << " | The reduced k-points are set according to symmetry operations. |" << std::endl; - GlobalV::ofs_running << " | We treat the spin as another set of k-points. |" << std::endl; - GlobalV::ofs_running << " | |" << std::endl; - GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; - GlobalV::ofs_running << "\n\n\n\n"; - - GlobalV::ofs_running << "\n SETUP K-POINTS" << std::endl; - - // (1) set nspin, read kpoints. - this->nspin = nspin_in; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nspin",nspin); - if(this->nspin==4) - { - this->nspin = 1;//zhengdy-soc - } - - bool read_succesfully = this->read_kpoints(k_file_name); + GlobalV::ofs_running << "\n\n\n\n"; + GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + GlobalV::ofs_running << " | |" << std::endl; + GlobalV::ofs_running << " | Setup K-points |" << std::endl; + GlobalV::ofs_running << " | We setup the k-points according to input parameters. |" << std::endl; + GlobalV::ofs_running << " | The reduced k-points are set according to symmetry operations. |" << std::endl; + GlobalV::ofs_running << " | We treat the spin as another set of k-points. |" << std::endl; + GlobalV::ofs_running << " | |" << std::endl; + GlobalV::ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + GlobalV::ofs_running << "\n\n\n\n"; + + GlobalV::ofs_running << "\n SETUP K-POINTS" << std::endl; + + // (1) set nspin, read kpoints. + this->nspin = nspin_in; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nspin", nspin); + if (this->nspin == 4) + { + this->nspin = 1; // zhengdy-soc + } + + bool read_succesfully = this->read_kpoints(k_file_name); #ifdef __MPI - Parallel_Common::bcast_bool(read_succesfully); + Parallel_Common::bcast_bool(read_succesfully); #endif - if(!read_succesfully) - { - ModuleBase::WARNING_QUIT("K_Vectors::set","Something wrong while reading KPOINTS."); - } + if (!read_succesfully) + { + ModuleBase::WARNING_QUIT("K_Vectors::set", "Something wrong while reading KPOINTS."); + } // output kpoints file - std::string skpt1=""; - std::string skpt2=""; + std::string skpt1 = ""; + std::string skpt2 = ""; // (2) - //only berry phase need all kpoints including time-reversal symmetry! - //if symm_flag is not set, only time-reversal symmetry would be considered. - if(!berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != -1) + // only berry phase need all kpoints including time-reversal symmetry! + // if symm_flag is not set, only time-reversal symmetry would be considered. + if (!berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag != -1) { bool match = true; this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, GlobalC::ucell, match); #ifdef __MPI - Parallel_Common::bcast_bool(match); + Parallel_Common::bcast_bool(match); #endif if (!match) { - std::cout<< "Optimized lattice type of reciprocal lattice cannot match the optimized real lattice. " <set_both_kvec(reciprocal_vec, latvec, skpt2); - if(GlobalV::MY_RANK==0) - { + if (GlobalV::MY_RANK == 0) + { // output kpoints file std::stringstream skpt; - skpt << GlobalV::global_readin_dir << "kpoints" ; - std::ofstream ofkpt( skpt.str().c_str()); // clear kpoints - ofkpt << skpt2 <normalize_wk(deg); + this->normalize_wk(deg); // It's very important in parallel case, // firstly do the mpi_k() and then // do set_kup_and_kdw() - GlobalC::Pkpoints.kinfo(nkstot); + GlobalC::Pkpoints.kinfo(nkstot); #ifdef __MPI - this->mpi_k();//2008-4-29 + this->mpi_k(); // 2008-4-29 #endif this->set_kup_and_kdw(); this->print_klists(GlobalV::ofs_running); - //std::cout << " NUMBER OF K-POINTS : " << nkstot << std::endl; + // std::cout << " NUMBER OF K-POINTS : " << nkstot << std::endl; #ifdef USE_PAW - GlobalC::paw_cell.set_isk(nks,isk.data()); + GlobalC::paw_cell.set_isk(nks, isk.data()); #endif return; } -void K_Vectors::renew(const int &kpoint_number) +void K_Vectors::renew(const int& kpoint_number) { kvec_c.resize(kpoint_number); kvec_d.resize(kpoint_number); @@ -190,50 +191,55 @@ void K_Vectors::renew(const int &kpoint_number) return; } -bool K_Vectors::read_kpoints(const std::string &fn) +bool K_Vectors::read_kpoints(const std::string& fn) { ModuleBase::TITLE("K_Vectors", "read_kpoints"); - if (GlobalV::MY_RANK != 0) return 1; - - // mohan add 2010-09-04 - if(GlobalV::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(); - } + if (GlobalV::MY_RANK != 0) + return 1; + + // mohan add 2010-09-04 + if (GlobalV::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 (GlobalV::KSPACING[0] > 0.0) { - if (GlobalV::KSPACING[1] <= 0 || GlobalV::KSPACING[2] <= 0){ - ModuleBase::WARNING_QUIT("K_Vectors","kspacing shold > 0"); + if (GlobalV::KSPACING[1] <= 0 || GlobalV::KSPACING[2] <= 0) + { + ModuleBase::WARNING_QUIT("K_Vectors", "kspacing shold > 0"); }; - //number of K points = max(1,int(|bi|/KSPACING+1)) + // number of K points = max(1,int(|bi|/KSPACING+1)) ModuleBase::Matrix3 btmp = GlobalC::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 / GlobalV::KSPACING[0] / GlobalC::ucell.lat0 + 1)); - int nk2 = std::max(1,static_cast(b2 * ModuleBase::TWO_PI / GlobalV::KSPACING[1] / GlobalC::ucell.lat0 + 1)); - int nk3 = std::max(1,static_cast(b3 * ModuleBase::TWO_PI / GlobalV::KSPACING[2] / GlobalC::ucell.lat0 + 1)); + int nk1 + = std::max(1, static_cast(b1 * ModuleBase::TWO_PI / GlobalV::KSPACING[0] / GlobalC::ucell.lat0 + 1)); + int nk2 + = std::max(1, static_cast(b2 * ModuleBase::TWO_PI / GlobalV::KSPACING[1] / GlobalC::ucell.lat0 + 1)); + int nk3 + = std::max(1, static_cast(b3 * ModuleBase::TWO_PI / GlobalV::KSPACING[2] / GlobalC::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; - ofs << "Gamma" << std::endl; - ofs << nk1 << " " << nk2 << " " << nk3 <<" 0 0 0" << std::endl; - ofs.close(); + std::ofstream ofs(fn.c_str()); + ofs << "K_POINTS" << std::endl; + ofs << "0" << std::endl; + ofs << "Gamma" << std::endl; + ofs << nk1 << " " << nk2 << " " << nk3 << " 0 0 0" << std::endl; + ofs.close(); } std::ifstream ifk(fn.c_str()); if (!ifk) - { - GlobalV::ofs_warning << " Can't find File name : " << fn << std::endl; - return 0; + { + GlobalV::ofs_warning << " Can't find File name : " << fn << std::endl; + return 0; } ifk >> std::setiosflags(std::ios::uppercase); @@ -251,8 +257,8 @@ bool K_Vectors::read_kpoints(const std::string &fn) 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" ) + 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; @@ -263,25 +269,25 @@ bool K_Vectors::read_kpoints(const std::string &fn) if (ierr == 0) { - GlobalV::ofs_warning << " symbol K_POINTS not found." << std::endl; - return 0; + GlobalV::ofs_warning << " symbol K_POINTS not found." << std::endl; + return 0; } - //input k-points are in 2pi/a units + // input k-points are in 2pi/a units ModuleBase::GlobalFunc::READ_VALUE(ifk, nkstot); - this->k_nkstot = nkstot; //LiuXh add 20180619 + this->k_nkstot = nkstot; // LiuXh add 20180619 - //std::cout << " nkstot = " << nkstot << std::endl; + // std::cout << " nkstot = " << nkstot << std::endl; ModuleBase::GlobalFunc::READ_VALUE(ifk, kword); - this->k_kword = kword; //LiuXh add 20180619 + this->k_kword = kword; // LiuXh add 20180619 - // mohan update 2021-02-22 - int max_kpoints = 100000; + // mohan update 2021-02-22 + int max_kpoints = 100000; if (nkstot > 100000) { - GlobalV::ofs_warning << " nkstot > MAX_KPOINTS" << std::endl; + GlobalV::ofs_warning << " nkstot > MAX_KPOINTS" << std::endl; return 0; } @@ -292,22 +298,22 @@ bool K_Vectors::read_kpoints(const std::string &fn) { is_mp = true; k_type = 0; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Input type of k points","Monkhorst-Pack(Gamma)"); + ModuleBase::GlobalFunc::OUT(GlobalV::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(GlobalV::ofs_running,"Input type of k points","Monkhorst-Pack"); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Input type of k points", "Monkhorst-Pack"); } else { - GlobalV::ofs_warning << " Error: neither Gamma nor Monkhorst-Pack." << std::endl; - return 0; + GlobalV::ofs_warning << " Error: neither Gamma nor Monkhorst-Pack." << std::endl; + return 0; } ifk >> nmp[0] >> nmp[1] >> nmp[2]; - + koffset[0] = 0; koffset[1] = 0; koffset[2] = 0; @@ -322,145 +328,149 @@ bool K_Vectors::read_kpoints(const std::string &fn) { if (kword == "Cartesian" || kword == "C") { - this->renew(nkstot * nspin);//mohan fix bug 2009-09-01 - for (int i = 0;i < nkstot;i++) + 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]); + ModuleBase::GlobalFunc::READ_VALUE(ifk, wk[i]); } this->kc_done = true; } else if (kword == "Direct" || kword == "D") { - this->renew(nkstot * nspin);//mohan fix bug 2009-09-01 - for (int i = 0;i < nkstot;i++) + 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]); + 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 0; - } + 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 0; + } Linely_add_k_between(ifk, kvec_c); - std::for_each(wk.begin(), wk.end(), [](double& d){d = 1.0;}); + std::for_each(wk.begin(), wk.end(), [](double& d) { d = 1.0; }); this->kc_done = true; + } - } - - 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 0; - } + 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 0; + } Linely_add_k_between(ifk, kvec_d); - std::for_each(wk.begin(), wk.end(), [](double& d){d = 1.0;}); + std::for_each(wk.begin(), wk.end(), [](double& d) { d = 1.0; }); this->kd_done = true; - } + } else { - GlobalV::ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; - return 0; + GlobalV::ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; + return 0; } } this->nkstot_full = this->nks = this->nkstot; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nkstot",nkstot); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nkstot", nkstot); return 1; } // END SUBROUTINE -void K_Vectors::Linely_add_k_between(std::ifstream& ifk, - std::vector>& kvec) { +void K_Vectors::Linely_add_k_between(std::ifstream& ifk, std::vector>& kvec) +{ // how many special points. int nks_special = this->nkstot; // number of points to the next k points - std::vector nkl(nks_special,0); + std::vector nkl(nks_special, 0); // coordinates of special points. std::vector> ks(nks_special); - //recalculate nkstot. + // recalculate nkstot. nkstot = 0; /* ISSUE#3482: to distinguish different kline segments */ std::vector kpt_segids; - kl_segids.clear(); kl_segids.shrink_to_fit(); + kl_segids.clear(); + kl_segids.shrink_to_fit(); int kpt_segid = 0; - for(int iks=0; iks> ks[iks].x; ifk >> ks[iks].y; ifk >> ks[iks].z; - ModuleBase::GlobalFunc::READ_VALUE( ifk, nkl[iks] ); + ModuleBase::GlobalFunc::READ_VALUE(ifk, nkl[iks]); assert(nkl[iks] >= 0); nkstot += nkl[iks]; /* ISSUE#3482: to distinguish different kline segments */ - if((nkl[iks] == 1)&&(iks!=(nks_special-1))) kpt_segid++; + if ((nkl[iks] == 1) && (iks != (nks_special - 1))) + kpt_segid++; kpt_segids.push_back(kpt_segid); } - assert( nkl[nks_special-1] == 1); + assert(nkl[nks_special - 1] == 1); - //std::cout << " nkstot = " << nkstot << std::endl; - this->renew(nkstot * nspin);//mohan fix bug 2009-09-01 + // std::cout << " nkstot = " << nkstot << std::endl; + this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 int count = 0; - for(int iks=1; iksrenew(nkstot * nspin); // mohan fix bug 2009-09-01 - for (int x = 1;x <= mpnx;x++) + for (int x = 1; x <= mpnx; x++) { - double v1 = Monkhorst_Pack_formula( k_type, koffset_in[0], x, mpnx); - if( std::abs(v1) < 1.0e-10 ) v1 = 0.0; //mohan update 2012-06-10 - for (int y = 1;y <= mpny;y++) + double v1 = Monkhorst_Pack_formula(k_type, koffset_in[0], x, mpnx); + if (std::abs(v1) < 1.0e-10) + v1 = 0.0; // mohan update 2012-06-10 + for (int y = 1; y <= mpny; y++) { - double v2 = Monkhorst_Pack_formula( k_type, koffset_in[1], y, mpny); - if( std::abs(v2) < 1.0e-10 ) v2 = 0.0; - for (int z = 1;z <= mpnz;z++) + double v2 = Monkhorst_Pack_formula(k_type, koffset_in[1], y, mpny); + if (std::abs(v2) < 1.0e-10) + v2 = 0.0; + for (int z = 1; z <= mpnz; z++) { - double v3 = Monkhorst_Pack_formula( k_type, koffset_in[2], z, mpnz); - if( std::abs(v3) < 1.0e-10 ) v3 = 0.0; + double v3 = Monkhorst_Pack_formula(k_type, koffset_in[2], z, mpnz); + if (std::abs(v3) < 1.0e-10) + v3 = 0.0; // index of nks kpoint const int i = mpnx * mpny * (z - 1) + mpnx * (y - 1) + (x - 1); kvec_d[i].set(v1, v2, v3); @@ -489,7 +502,7 @@ void K_Vectors::Monkhorst_Pack(const int *nmp_in, const double *koffset_in, cons } const double weight = 1.0 / static_cast(nkstot); - for (int ik=0; ik 0 ); + if (GlobalV::MY_RANK != 0) + return; + ModuleBase::TITLE("K_Vectors", "update_use_ibz"); + assert(nkstot_ibz > 0); // update nkstot this->nkstot = this->nkstot_ibz; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nkstot now",nkstot); + ModuleBase::GlobalFunc::OUT(GlobalV::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() + this->kvec_d.resize(this->nkstot * nspin); // qianrui fix a bug 2021-7-13 for nspin=2 in set_kup_and_kdw() for (int i = 0; i < this->nkstot; ++i) { this->kvec_d[i] = this->kvec_d_ibz[i]; - // update weight. + // update weight. this->wk[i] = this->wk_ibz[i]; } @@ -524,9 +538,14 @@ void K_Vectors::update_use_ibz( void ) return; } -void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,std::string& skpt, const UnitCell &ucell, bool& match) +void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry& symm, + bool use_symm, + std::string& skpt, + const UnitCell& ucell, + bool& match) { - if (GlobalV::MY_RANK != 0) return; + if (GlobalV::MY_RANK != 0) + return; ModuleBase::TITLE("K_Vectors", "ibz_kpoint"); // k-lattice: "pricell" of reciprocal space @@ -544,7 +563,6 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s gk = ModuleBase::Matrix3(gk1.x, gk1.y, gk1.z, gk2.x, gk2.y, gk2.z, gk3.x, gk3.y, gk3.z); } - //=============================================== // search in all space group operations // if the operations does not already included @@ -555,43 +573,57 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s ModuleBase::Matrix3 inv(-1, 0, 0, 0, -1, 0, 0, 0, -1); ModuleBase::Matrix3 ind(1, 0, 0, 0, 1, 0, 0, 0, 1); - int nrotkm=0; - if(use_symm) + int nrotkm = 0; + if (use_symm) { // bravais type of reciprocal lattice and k-lattice double b_const[6]; double b0_const[6]; double bk_const[6]; double bk0_const[6]; - int bbrav=15; - int bkbrav=15; + int bbrav = 15; + int bkbrav = 15; std::string bbrav_name; std::string bkbrav_name; - ModuleBase::Vector3gk01 = gk1, gk02 = gk2, gk03 = gk3; + ModuleBase::Vector3 gk01 = gk1, gk02 = gk2, gk03 = gk3; ModuleBase::Matrix3 b_optlat = symm.optlat.Inverse().Transpose(); - //search optlat after using reciprocity relation + // search optlat after using reciprocity relation ModuleBase::Vector3 gb01(b_optlat.e11, b_optlat.e12, b_optlat.e13); ModuleBase::Vector3 gb02(b_optlat.e21, b_optlat.e22, b_optlat.e23); ModuleBase::Vector3 gb03(b_optlat.e31, b_optlat.e32, b_optlat.e33); - symm.lattice_type(gb1, gb2, gb3, gb01, gb02, gb03, b_const, b0_const, bbrav, bbrav_name, ucell.atoms, false, nullptr); - GlobalV::ofs_running<<"(for reciprocal lattice: )"< ibrav_a2b{ 1, 3, 2, 4, 5, 6, 7, 8, 10, 9, 11, 12, 13, 14 }; - auto ibrav_match = [&](int ibrav_b) -> bool - { + std::vector ibrav_a2b{1, 3, 2, 4, 5, 6, 7, 8, 10, 9, 11, 12, 13, 14}; + auto ibrav_match = [&](int ibrav_b) -> bool { const int& ibrav_a = symm.real_brav; - if (ibrav_a < 1 || ibrav_a > 14) return false; + if (ibrav_a < 1 || ibrav_a > 14) + return false; return (ibrav_b == ibrav_a2b[ibrav_a - 1]); }; if (!ibrav_match(bbrav)) { - GlobalV::ofs_running << "Error: Bravais lattice type of reciprocal lattice is not compatible with that of real space lattice:" << std::endl; + GlobalV::ofs_running << "Error: Bravais lattice type of reciprocal lattice is not compatible with that of " + "real space lattice:" + << std::endl; GlobalV::ofs_running << "ibrav of real space lattice: " << symm.ilattname << std::endl; GlobalV::ofs_running << "ibrav of reciprocal lattice: " << bbrav_name << std::endl; GlobalV::ofs_running << "(which should be " << ibrav_a2b[symm.real_brav - 1] << ")." << std::endl; @@ -600,7 +632,19 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s } if (this->is_mp) { - symm.lattice_type(gk1, gk2, gk3, gk01, gk02, gk03, bk_const, bk0_const, bkbrav, bkbrav_name, ucell.atoms, false, nullptr); + symm.lattice_type(gk1, + gk2, + gk3, + gk01, + gk02, + gk03, + bk_const, + bk0_const, + bkbrav, + bkbrav_name, + ucell.atoms, + false, + nullptr); GlobalV::ofs_running << "(for k-lattice: )" << std::endl; ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS TYPE", bkbrav); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS LATTICE NAME", bkbrav_name); @@ -610,28 +654,44 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s ModuleBase::Matrix3 bsymop[48]; int bnop = 0; // search again - symm.lattice_type(gb1, gb2, gb3, gb1, gb2, gb3, b_const, b0_const, bbrav, bbrav_name, ucell.atoms, false, nullptr); + symm.lattice_type(gb1, + gb2, + gb3, + gb1, + gb2, + gb3, + b_const, + b0_const, + bbrav, + bbrav_name, + ucell.atoms, + false, + nullptr); ModuleBase::Matrix3 b_optlat_new(gb1.x, gb1.y, gb1.z, gb2.x, gb2.y, gb2.z, gb3.x, gb3.y, gb3.z); symm.setgroup(bsymop, bnop, bbrav); symm.gmatrix_convert(bsymop, bsymop, bnop, b_optlat_new, ucell.G); - - //check if all the kgmatrix are in bsymop - auto matequal = [&symm] (ModuleBase::Matrix3 a, ModuleBase::Matrix3 b) - { - return (symm.equal(a.e11, b.e11) && symm.equal(a.e12, b.e12) && symm.equal(a.e13, b.e13) && - symm.equal(a.e21, b.e21) && symm.equal(a.e22, b.e22) && symm.equal(a.e23, b.e23) && - symm.equal(a.e31, b.e31) && symm.equal(a.e23, b.e23) && symm.equal(a.e33, b.e33)); + + // check if all the kgmatrix are in bsymop + auto matequal = [&symm](ModuleBase::Matrix3 a, ModuleBase::Matrix3 b) { + return (symm.equal(a.e11, b.e11) && symm.equal(a.e12, b.e12) && symm.equal(a.e13, b.e13) + && symm.equal(a.e21, b.e21) && symm.equal(a.e22, b.e22) && symm.equal(a.e23, b.e23) + && symm.equal(a.e31, b.e31) && symm.equal(a.e23, b.e23) && symm.equal(a.e33, b.e33)); }; - for(int i=0;iis_mp)symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, gk); + if (this->is_mp) + symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, gk); // direct coordinates of k-points in k-lattice std::vector> kvec_d_k(nkstot); - if (this->is_mp) for (int i = 0;i < nkstot;++i) kvec_d_k[i] = kvec_d[i] * ucell.G * gk.Inverse(); + if (this->is_mp) + for (int i = 0; i < nkstot; ++i) + kvec_d_k[i] = kvec_d[i] * ucell.G * gk.Inverse(); // use operation : kgmatrix to find // the new set kvec_d : ir_kpt @@ -677,43 +740,45 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s wk_ibz.resize(this->nkstot); ibz2bz.resize(this->nkstot); - // nkstot is the total input k-points number. + // nkstot is the total input k-points number. const double weight = 1.0 / static_cast(nkstot); ModuleBase::Vector3 kvec_rot; ModuleBase::Vector3 kvec_rot_k; - -// for(int i=0; i &kvec){ + // for(int i=0; i& kvec) { // in (-0.5, 0.5] - kvec.x = fmod(kvec.x + 100.5-0.5*symm.epsilon, 1)-0.5+0.5*symm.epsilon; - kvec.y = fmod(kvec.y + 100.5-0.5*symm.epsilon, 1)-0.5+0.5*symm.epsilon; - kvec.z = fmod(kvec.z + 100.5-0.5*symm.epsilon, 1)-0.5+0.5*symm.epsilon; + kvec.x = fmod(kvec.x + 100.5 - 0.5 * symm.epsilon, 1) - 0.5 + 0.5 * symm.epsilon; + kvec.y = fmod(kvec.y + 100.5 - 0.5 * symm.epsilon, 1) - 0.5 + 0.5 * symm.epsilon; + kvec.z = fmod(kvec.z + 100.5 - 0.5 * symm.epsilon, 1) - 0.5 + 0.5 * symm.epsilon; // in [0, 1) // kvec.x = fmod(kvec.x + 100 + symm.epsilon, 1) - symm.epsilon; // kvec.y = fmod(kvec.y + 100 + symm.epsilon, 1) - symm.epsilon; // kvec.z = fmod(kvec.z + 100 + symm.epsilon, 1) - symm.epsilon; - if(std::abs(kvec.x) < symm.epsilon) kvec.x = 0.0; - if(std::abs(kvec.y) < symm.epsilon) kvec.y = 0.0; - if(std::abs(kvec.z) < symm.epsilon) kvec.z = 0.0; + if (std::abs(kvec.x) < symm.epsilon) + kvec.x = 0.0; + if (std::abs(kvec.y) < symm.epsilon) + kvec.y = 0.0; + if (std::abs(kvec.z) < symm.epsilon) + kvec.z = 0.0; return; }; // for output in kpoints file int ibz_index[nkstot]; - // search in all k-poins. + // search in all k-poins. for (int i = 0; i < nkstot; ++i) { - //restrict to [0, 1) + // restrict to [0, 1) restrict_kpt(kvec_d[i]); - //std::cout << "\n kpoint = " << i << std::endl; - //std::cout << "\n kvec_d = " << kvec_d[i].x << " " << kvec_d[i].y << " " << kvec_d[i].z; + // std::cout << "\n kpoint = " << i << std::endl; + // std::cout << "\n kvec_d = " << kvec_d[i].x << " " << kvec_d[i].y << " " << kvec_d[i].z; bool already_exist = false; - int exist_number = -1; + int exist_number = -1; // search over all symmetry operations for (int j = 0; j < nrotkm; ++j) { @@ -721,148 +786,162 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s { // rotate the kvec_d within all operations. // here use direct coordinates. -// kvec_rot = kgmatrix[j] * kvec_d[i]; - // mohan modify 2010-01-30. - // mohan modify again 2010-01-31 - // fix the bug like kvec_d * G; is wrong - kvec_rot = kvec_d[i] * kgmatrix[j]; //wrong for total energy, but correct for nonlocal force. - //kvec_rot = kgmatrix[j] * kvec_d[i]; //correct for total energy, but wrong for nonlocal force. + // kvec_rot = kgmatrix[j] * kvec_d[i]; + // mohan modify 2010-01-30. + // mohan modify again 2010-01-31 + // fix the bug like kvec_d * G; is wrong + kvec_rot = kvec_d[i] * kgmatrix[j]; // wrong for total energy, but correct for nonlocal force. + // kvec_rot = kgmatrix[j] * kvec_d[i]; //correct for total energy, but wrong for nonlocal force. restrict_kpt(kvec_rot); if (this->is_mp) { - kvec_rot_k = kvec_d_k[i] * kkmatrix[j]; //k-lattice rotation - kvec_rot_k = kvec_rot_k * gk * ucell.G.Inverse(); //convert to recip lattice + kvec_rot_k = kvec_d_k[i] * kkmatrix[j]; // k-lattice rotation + kvec_rot_k = kvec_rot_k * gk * ucell.G.Inverse(); // convert to recip lattice restrict_kpt(kvec_rot_k); assert(symm.equal(kvec_rot.x, kvec_rot_k.x)); assert(symm.equal(kvec_rot.y, kvec_rot_k.y)); assert(symm.equal(kvec_rot.z, kvec_rot_k.z)); // std::cout << "\n kvec_rot (in recip) = " << kvec_rot.x << " " << kvec_rot.y << " " << kvec_rot.z; - // std::cout << "\n kvec_rot(k to recip)= " << kvec_rot_k.x << " " << kvec_rot_k.y << " " << kvec_rot_k.z; - kvec_rot_k = kvec_rot_k * ucell.G * gk.Inverse(); //convert back to k-latice + // std::cout << "\n kvec_rot(k to recip)= " << kvec_rot_k.x << " " << kvec_rot_k.y << " " << + // kvec_rot_k.z; + kvec_rot_k = kvec_rot_k * ucell.G * gk.Inverse(); // convert back to k-latice } - for (int k=0; k< this->nkstot_ibz; ++k) + for (int k = 0; k < this->nkstot_ibz; ++k) { - if ( symm.equal(kvec_rot.x, this->kvec_d_ibz[k].x) && - symm.equal(kvec_rot.y, this->kvec_d_ibz[k].y) && - symm.equal(kvec_rot.z, this->kvec_d_ibz[k].z)) + if (symm.equal(kvec_rot.x, this->kvec_d_ibz[k].x) && symm.equal(kvec_rot.y, this->kvec_d_ibz[k].y) + && symm.equal(kvec_rot.z, this->kvec_d_ibz[k].z)) { already_exist = true; - // find another ibz k point, - // but is already in the ibz_kpoint list. - // so the weight need to +1; + // find another ibz k point, + // but is already in the ibz_kpoint list. + // so the weight need to +1; this->wk_ibz[k] += weight; - exist_number = k; + exist_number = k; break; } } - }//end !already_exist + } // end !already_exist } // if really there is no equivalent k point in the list, then add it. if (!already_exist) { - //if it's a new ibz kpoint. - //nkstot_ibz indicate the index of ibz kpoint. + // if it's a new ibz kpoint. + // nkstot_ibz indicate the index of ibz kpoint. this->kvec_d_ibz[nkstot_ibz] = kvec_rot; // output in kpoints file ibz_index[i] = nkstot_ibz; - //the weight should be averged k-point weight. + // the weight should be averged k-point weight. this->wk_ibz[nkstot_ibz] = weight; - //ibz2bz records the index of origin k points. + // ibz2bz records the index of origin k points. this->ibz2bz[nkstot_ibz] = i; ++nkstot_ibz; } - else //mohan fix bug 2010-1-30 - { -// std::cout << "\n\n already exist ! "; + else // mohan fix bug 2010-1-30 + { + // std::cout << "\n\n already exist ! "; -// std::cout << "\n kvec_rot = " << kvec_rot.x << " " << kvec_rot.y << " " << kvec_rot.z; -// std::cout << "\n kvec_d_ibz = " << kvec_d_ibz[exist_number].x -// << " " << kvec_d_ibz[exist_number].y -// << " " << kvec_d_ibz[exist_number].z; + // std::cout << "\n kvec_rot = " << kvec_rot.x << " " << kvec_rot.y << " " << kvec_rot.z; + // std::cout << "\n kvec_d_ibz = " << kvec_d_ibz[exist_number].x + // << " " << kvec_d_ibz[exist_number].y + // << " " << kvec_d_ibz[exist_number].z; - double kmol_new = kvec_d[i].norm2(); - double kmol_old = kvec_d_ibz[exist_number].norm2(); + double kmol_new = kvec_d[i].norm2(); + double kmol_old = kvec_d_ibz[exist_number].norm2(); ibz_index[i] = exist_number; -// std::cout << "\n kmol_new = " << kmol_new; -// std::cout << "\n kmol_old = " << kmol_old; - - - // why we need this step? - // because in pw_basis.cpp, while calculate ggwfc2, - // if we want to keep the result of symmetry operation is right. - // we need to fix the number of plane wave. - // and the number of plane wave is depending on the |K+G|, - // so we need to |K|max to be the same as 'no symmetry'. - // mohan 2010-01-30 - if(kmol_new > kmol_old) - { - kvec_d_ibz[exist_number] = kvec_d[i]; + // std::cout << "\n kmol_new = " << kmol_new; + // std::cout << "\n kmol_old = " << kmol_old; + + // why we need this step? + // because in pw_basis.cpp, while calculate ggwfc2, + // if we want to keep the result of symmetry operation is right. + // we need to fix the number of plane wave. + // and the number of plane wave is depending on the |K+G|, + // so we need to |K|max to be the same as 'no symmetry'. + // mohan 2010-01-30 + if (kmol_new > kmol_old) + { + kvec_d_ibz[exist_number] = kvec_d[i]; } - } -// BLOCK_HERE("check k point"); + } + // BLOCK_HERE("check k point"); } delete[] kkmatrix; // output in kpoints file std::stringstream ss; - ss << " " << std::setw(40) <<"nkstot" << " = " << nkstot - << std::setw(66) << "ibzkpt" << std::endl; + 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"); + "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, this->kvec_d[i].x, this->kvec_d[i].y, this->kvec_d[i].z, - ibz_index[i]+1, this->kvec_d_ibz[ibz_index[i]].x, this->kvec_d_ibz[ibz_index[i]].y, this->kvec_d_ibz[ibz_index[i]].z); + i + 1, + this->kvec_d[i].x, + this->kvec_d[i].y, + this->kvec_d[i].z, + ibz_index[i] + 1, + this->kvec_d_ibz[ibz_index[i]].x, + this->kvec_d_ibz[ibz_index[i]].y, + this->kvec_d_ibz[ibz_index[i]].z); } ss << table << std::endl; skpt = ss.str(); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nkstot_ibz", nkstot_ibz); - + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nkstot_ibz", nkstot_ibz); + table.clear(); table += "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; ikkvec_d_ibz[ik].x, this->kvec_d_ibz[ik].y, this->kvec_d_ibz[ik].z, - this->wk_ibz[ik], this->ibz2bz[ik]); + ik + 1, + this->kvec_d_ibz[ik].x, + this->kvec_d_ibz[ik].y, + this->kvec_d_ibz[ik].z, + this->wk_ibz[ik], + this->ibz2bz[ik]); } GlobalV::ofs_running << table << std::endl; return; } - -void K_Vectors::set_both_kvec(const ModuleBase::Matrix3 &G, const ModuleBase::Matrix3 &R,std::string& skpt) +void K_Vectors::set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt) { - if(GlobalV::FINAL_SCF) //LiuXh add 20180606 + if (GlobalV::FINAL_SCF) // LiuXh add 20180606 { - if(k_nkstot == 0) + if (k_nkstot == 0) { kd_done = true; kc_done = false; } else { - if(k_kword == "Cartesian" || k_kword == "C") - { - kc_done = true; - kd_done = false; - } - else if (k_kword == "Direct" || k_kword == "D") - { - kd_done = true; - kc_done = false; - } + if (k_kword == "Cartesian" || k_kword == "C") + { + kc_done = true; + kd_done = false; + } + else if (k_kword == "Direct" || k_kword == "D") + { + kd_done = true; + kc_done = false; + } else { GlobalV::ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; @@ -873,20 +952,26 @@ void K_Vectors::set_both_kvec(const ModuleBase::Matrix3 &G, const ModuleBase::Ma // set cartesian k vectors. if (!kc_done && kd_done) { - for (int i = 0;i < nkstot;i++) + for (int i = 0; i < nkstot; i++) { -//wrong!! kvec_c[i] = G * kvec_d[i]; -// mohan fixed bug 2010-1-10 - if( std::abs(kvec_d[i].x) < 1.0e-10 ) kvec_d[i].x = 0.0; - if( std::abs(kvec_d[i].y) < 1.0e-10 ) kvec_d[i].y = 0.0; - if( std::abs(kvec_d[i].z) < 1.0e-10 ) kvec_d[i].z = 0.0; - - kvec_c[i] = kvec_d[i] * G; - - // mohan add2012-06-10 - if( std::abs(kvec_c[i].x) < 1.0e-10 ) kvec_c[i].x = 0.0; - if( std::abs(kvec_c[i].y) < 1.0e-10 ) kvec_c[i].y = 0.0; - if( std::abs(kvec_c[i].z) < 1.0e-10 ) kvec_c[i].z = 0.0; + // wrong!! kvec_c[i] = G * kvec_d[i]; + // mohan fixed bug 2010-1-10 + if (std::abs(kvec_d[i].x) < 1.0e-10) + kvec_d[i].x = 0.0; + if (std::abs(kvec_d[i].y) < 1.0e-10) + kvec_d[i].y = 0.0; + if (std::abs(kvec_d[i].z) < 1.0e-10) + kvec_d[i].z = 0.0; + + kvec_c[i] = kvec_d[i] * G; + + // mohan add2012-06-10 + if (std::abs(kvec_c[i].x) < 1.0e-10) + kvec_c[i].x = 0.0; + if (std::abs(kvec_c[i].y) < 1.0e-10) + kvec_c[i].y = 0.0; + if (std::abs(kvec_c[i].z) < 1.0e-10) + kvec_c[i].z = 0.0; } kc_done = true; } @@ -895,56 +980,60 @@ void K_Vectors::set_both_kvec(const ModuleBase::Matrix3 &G, const ModuleBase::Ma else if (kc_done && !kd_done) { ModuleBase::Matrix3 RT = R.Transpose(); - for (int i = 0;i < nkstot;i++) + for (int i = 0; i < nkstot; i++) { -// std::cout << " ik=" << i -// << " kvec.x=" << kvec_c[i].x -// << " kvec.y=" << kvec_c[i].y -// << " kvec.z=" << kvec_c[i].z << std::endl; -//wrong! kvec_d[i] = RT * kvec_c[i]; -// mohan fixed bug 2011-03-07 + // std::cout << " ik=" << i + // << " kvec.x=" << kvec_c[i].x + // << " kvec.y=" << kvec_c[i].y + // << " kvec.z=" << kvec_c[i].z << std::endl; + // wrong! kvec_d[i] = RT * kvec_c[i]; + // mohan fixed bug 2011-03-07 kvec_d[i] = kvec_c[i] * RT; } kd_done = true; } std::string table; table += "K-POINTS DIRECT COORDINATES\n"; - table += FmtCore::format("%8s%12s%12s%12s%8s\n", - "KPOINTS", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT"); - for(int i=0; ikvec_d[i].x, this->kvec_d[i].y, this->kvec_d[i].z, this->wk[i]); - } + i + 1, + this->kvec_d[i].x, + this->kvec_d[i].y, + this->kvec_d[i].z, + this->wk[i]); + } GlobalV::ofs_running << table << std::endl; - if(GlobalV::MY_RANK==0) - { + if (GlobalV::MY_RANK == 0) + { std::stringstream ss; - ss << " " << std::setw(40) <<"nkstot now" << " = " << nkstot << std::endl; + ss << " " << std::setw(40) << "nkstot now" + << " = " << nkstot << std::endl; ss << table << std::endl; skpt = ss.str(); - } + } return; } - -void K_Vectors::normalize_wk(const int °spin) +void K_Vectors::normalize_wk(const int& degspin) { - if(GlobalV::MY_RANK!=0) return; + if (GlobalV::MY_RANK != 0) + return; double sum = 0.0; - for (int ik = 0;ik < nkstot;ik++) + for (int ik = 0; ik < nkstot; ik++) { sum += this->wk[ik]; } - assert(sum>0.0); + assert(sum > 0.0); - for (int ik = 0;ik < nkstot;ik++) + for (int ik = 0; ik < nkstot; ik++) { this->wk[ik] /= sum; } - for (int ik = 0;ik < nkstot;ik++) + for (int ik = 0; ik < nkstot; ik++) { this->wk[ik] *= degspin; } @@ -955,7 +1044,7 @@ void K_Vectors::normalize_wk(const int °spin) #ifdef __MPI void K_Vectors::mpi_k(void) { - ModuleBase::TITLE("K_Vectors","mpi_k"); + ModuleBase::TITLE("K_Vectors", "mpi_k"); Parallel_Common::bcast_bool(kc_done); @@ -976,69 +1065,68 @@ void K_Vectors::mpi_k(void) this->nks = GlobalC::Pkpoints.nks_pool[GlobalV::MY_POOL]; - GlobalV::ofs_running << std::endl; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"k-point number in this process",nks); + GlobalV::ofs_running << std::endl; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "k-point number in this process", nks); int nks_minimum = this->nks; - Parallel_Reduce::gather_min_int_all( nks_minimum ); + Parallel_Reduce::gather_min_int_all(nks_minimum); if (nks_minimum == 0) { - ModuleBase::WARNING_QUIT("K_Vectors::mpi_k()"," nks == 0, some processor have no k point!"); + ModuleBase::WARNING_QUIT("K_Vectors::mpi_k()", " nks == 0, some processor have no k point!"); } else { - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"minimum distributed K point number",nks_minimum); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "minimum distributed K point number", nks_minimum); } std::vector isk_aux(nkstot); std::vector wk_aux(nkstot); - std::vector kvec_c_aux(nkstot*3); - std::vector kvec_d_aux(nkstot*3); + std::vector kvec_c_aux(nkstot * 3); + std::vector kvec_d_aux(nkstot * 3); if (GlobalV::MY_RANK == 0) { - for (int ik = 0;ik < nkstot;ik++) + 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_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; } } Parallel_Common::bcast_int(isk_aux.data(), nkstot); Parallel_Common::bcast_double(wk_aux.data(), nkstot); - Parallel_Common::bcast_double(kvec_c_aux.data(), nkstot*3); - Parallel_Common::bcast_double(kvec_d_aux.data(), nkstot*3); + Parallel_Common::bcast_double(kvec_c_aux.data(), nkstot * 3); + Parallel_Common::bcast_double(kvec_d_aux.data(), nkstot * 3); this->renew(this->nks * this->nspin); // distribute int k_index = 0; - for (int i = 0;i < nks;i++) + for (int i = 0; i < nks; i++) { // 3 is because each k point has three value:kx, ky, kz - k_index = i + GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL] ; - 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]; + k_index = i + GlobalC::Pkpoints.startk_pool[GlobalV::MY_POOL]; + 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]; wk[i] = wk_aux[k_index]; isk[i] = isk_aux[k_index]; } } // END SUBROUTINE #endif - //---------------------------------------------------------- // This routine sets the k vectors for the up and down spin //---------------------------------------------------------- @@ -1067,18 +1155,18 @@ void K_Vectors::set_kup_and_kdw(void) 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->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; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nks(nspin=2)",nks); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nkstot(nspin=2)",nkstot); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nks(nspin=2)", nks); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nkstot(nspin=2)", nkstot); break; case 4: @@ -1093,8 +1181,7 @@ void K_Vectors::set_kup_and_kdw(void) return; } // end subroutine set_kup_and_kdw - -void K_Vectors::print_klists(std::ofstream &ofs) +void K_Vectors::print_klists(std::ofstream& ofs) { ModuleBase::TITLE("K_Vectors", "print_klists"); @@ -1102,64 +1189,75 @@ void K_Vectors::print_klists(std::ofstream &ofs) { std::cout << "\n nkstot=" << nkstot; std::cout << "\n nks=" << nks; - ModuleBase::WARNING_QUIT("print_klists","nkstot < nks"); + ModuleBase::WARNING_QUIT("print_klists", "nkstot < nks"); } std::string table; table += "K-POINTS CARTESIAN COORDINATES\n"; - table += FmtCore::format("%8s%12s%12s%12s%8s\n", - "KPOINTS", "CARTESIAN_X", "CARTESIAN_Y", "CARTESIAN_Z", "WEIGHT"); - for(int i=0; ikvec_c[i].x, this->kvec_c[i].y, this->kvec_c[i].z, this->wk[i]); - } + i + 1, + this->kvec_c[i].x, + this->kvec_c[i].y, + this->kvec_c[i].z, + this->wk[i]); + } GlobalV::ofs_running << "\n" << table << std::endl; table.clear(); table += "K-POINTS DIRECT COORDINATES\n"; - table += FmtCore::format("%8s%12s%12s%12s%8s\n", - "KPOINTS", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT"); - for(int i=0; ikvec_d[i].x, this->kvec_d[i].y, this->kvec_d[i].z, this->wk[i]); - } + i + 1, + this->kvec_d[i].x, + this->kvec_d[i].y, + this->kvec_d[i].z, + this->wk[i]); + } GlobalV::ofs_running << "\n" << table << std::endl; return; } -//LiuXh add a new function here, -//20180515 -void K_Vectors::set_after_vc( - const int& nspin_in, - const ModuleBase::Matrix3 &reciprocal_vec, - const ModuleBase::Matrix3 &latvec) +// LiuXh add a new function here, +// 20180515 +void K_Vectors::set_after_vc(const int& nspin_in, + const ModuleBase::Matrix3& reciprocal_vec, + const ModuleBase::Matrix3& latvec) { ModuleBase::TITLE("K_Vectors", "set_after_vc"); GlobalV::ofs_running << "\n SETUP K-POINTS" << std::endl; this->nspin = nspin_in; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nspin",nspin); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nspin", nspin); // set cartesian k vectors. kd_done = true; kc_done = false; if (!kc_done && kd_done) { - for (int i = 0;i < nks;i++) + for (int i = 0; i < nks; i++) { -//wrong!! kvec_c[i] = G * kvec_d[i]; -// mohan fixed bug 2010-1-10 - if( std::abs(kvec_d[i].x) < 1.0e-10 ) kvec_d[i].x = 0.0; - if( std::abs(kvec_d[i].y) < 1.0e-10 ) kvec_d[i].y = 0.0; - if( std::abs(kvec_d[i].z) < 1.0e-10 ) kvec_d[i].z = 0.0; - - kvec_c[i] = kvec_d[i] * reciprocal_vec; - - // mohan add2012-06-10 - if( std::abs(kvec_c[i].x) < 1.0e-10 ) kvec_c[i].x = 0.0; - if( std::abs(kvec_c[i].y) < 1.0e-10 ) kvec_c[i].y = 0.0; - if( std::abs(kvec_c[i].z) < 1.0e-10 ) kvec_c[i].z = 0.0; + // wrong!! kvec_c[i] = G * kvec_d[i]; + // mohan fixed bug 2010-1-10 + if (std::abs(kvec_d[i].x) < 1.0e-10) + kvec_d[i].x = 0.0; + if (std::abs(kvec_d[i].y) < 1.0e-10) + kvec_d[i].y = 0.0; + if (std::abs(kvec_d[i].z) < 1.0e-10) + kvec_d[i].z = 0.0; + + kvec_c[i] = kvec_d[i] * reciprocal_vec; + + // mohan add2012-06-10 + if (std::abs(kvec_c[i].x) < 1.0e-10) + kvec_c[i].x = 0.0; + if (std::abs(kvec_c[i].y) < 1.0e-10) + kvec_c[i].y = 0.0; + if (std::abs(kvec_c[i].z) < 1.0e-10) + kvec_c[i].z = 0.0; } kc_done = true; } @@ -1168,32 +1266,34 @@ void K_Vectors::set_after_vc( else if (kc_done && !kd_done) { ModuleBase::Matrix3 RT = latvec.Transpose(); - for (int i = 0;i < nks;i++) + for (int i = 0; i < nks; i++) { -// std::cout << " ik=" << i -// << " kvec.x=" << kvec_c[i].x -// << " kvec.y=" << kvec_c[i].y -// << " kvec.z=" << kvec_c[i].z << std::endl; -//wrong! kvec_d[i] = RT * kvec_c[i]; -// mohan fixed bug 2011-03-07 + // std::cout << " ik=" << i + // << " kvec.x=" << kvec_c[i].x + // << " kvec.y=" << kvec_c[i].y + // << " kvec.z=" << kvec_c[i].z << std::endl; + // wrong! kvec_d[i] = RT * kvec_c[i]; + // mohan fixed bug 2011-03-07 kvec_d[i] = kvec_c[i] * RT; } kd_done = true; } std::string table; table += "K-POINTS DIRECT COORDINATES\n"; - table += FmtCore::format("%8s%12s%12s%12s%8s\n", - "KPOINTS", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT"); - for(int i=0; ikvec_d[i].x, this->kvec_d[i].y, this->kvec_d[i].z, this->wk[i]); - } + i + 1, + this->kvec_d[i].x, + this->kvec_d[i].y, + this->kvec_d[i].z, + this->wk[i]); + } GlobalV::ofs_running << table << std::endl; - //this->set_both_kvec(reciprocal_vec, latvec); + // this->set_both_kvec(reciprocal_vec, latvec); this->print_klists(GlobalV::ofs_running); return; } - diff --git a/source/module_cell/klist.h b/source/module_cell/klist.h index e19b08eba9..c342eca4bb 100644 --- a/source/module_cell/klist.h +++ b/source/module_cell/klist.h @@ -1,102 +1,367 @@ #ifndef K_VECTORS_H #define K_VECTORS_H -#include - #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_base/matrix3.h" #include "module_cell/unitcell.h" +#include + class K_Vectors { -public: - - std::vector> kvec_c; // Cartesian coordinates of k points - std::vector> kvec_d; // Direct coordinates of k points - std::vector> kvec_d_ibz; // ibz Direct coordinates of k points + public: + std::vector> kvec_c; /// Cartesian coordinates of k points + std::vector> kvec_d; /// Direct coordinates of k points + std::vector> kvec_d_ibz; /// ibz Direct coordinates of k points - std::vector wk; // wk, weight of k points - std::vector wk_ibz; // ibz kpoint wk ,weight of k points + std::vector wk; /// wk, weight of k points + std::vector wk_ibz; /// ibz kpoint wk ,weight of k points - std::vector ngk; // ngk, number of plane waves for each k point - std::vector isk; // distinguish spin up and down k points - std::vector ibz2bz; // mohan added 2009-05-18 + std::vector ngk; /// ngk, number of plane waves for each k point + std::vector isk; /// distinguish spin up and down k points - int nks; // number of k points in this pool(processor, up+dw) - int nkstot; /// total number of k points, equal to nkstot_ibz after reducing k points - int nkstot_ibz; /// number of k points in IBZ - int nkstot_full; /// number of k points in full k mesh + int nks; /// number of k points in this pool(processor, up+dw) + int nkstot; /// total number of k points, equal to nkstot_ibz after reducing k points + int nkstot_ibz; /// number of k points in IBZ + int nkstot_full; /// number of k points in full k mesh - int nmp[3]; // Number of Monhorst-Pack - std::vector kl_segids; // index of kline segment + int nmp[3]; /// Number of Monhorst-Pack + std::vector kl_segids; /// index of kline segment K_Vectors(); ~K_Vectors(); - void set( - const ModuleSymmetry::Symmetry &symm, - const std::string &k_file_name, - const int& nspin, - const ModuleBase::Matrix3 &reciprocal_vec, - const ModuleBase::Matrix3 &latvec); - - void ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,std::string& skpt, const UnitCell &ucell, bool& match); - //LiuXh add 20180515 - void set_after_vc( - const int& nspin, - const ModuleBase::Matrix3 &reciprocal_vec, - const ModuleBase::Matrix3 &latvec); - //get global index for ik + /** + * @brief Set up the k-points for the system. + * + * This function sets up the k-points according to the input parameters and symmetry operations. + * It also treats the spin as another set of k-points. + * + * @param symm The symmetry of the system. + * @param k_file_name The name of the file containing the k-points. + * @param nspin_in The number of spins. + * @param reciprocal_vec The reciprocal vector of the system. + * @param latvec The lattice vector of the system. + * + * @return void + * + * @note This function will quit with a warning if something goes wrong while reading the KPOINTS file. + * @note If the optimized lattice type of the reciprocal lattice cannot match the optimized real lattice, + * it will output a warning and suggest possible solutions. + * @note Only available for nspin = 1 or 2 or 4. + */ + void set(const ModuleSymmetry::Symmetry& symm, + const std::string& k_file_name, + const int& nspin, + const ModuleBase::Matrix3& reciprocal_vec, + const ModuleBase::Matrix3& latvec); + + /** + * @brief Generates irreducible k-points in the Brillouin zone considering symmetry operations. + * + * This function calculates the irreducible k-points (IBZ) from the given k-points, taking into + * account the symmetry of the unit cell. It updates the symmetry-matched k-points and generates + * the corresponding weight for each k-point. + * + * @param symm The symmetry information of the system. + * @param use_symm A flag indicating whether to use symmetry operations. + * @param skpt A string to store the formatted k-points information. + * @param ucell The unit cell of the crystal. + * @param match A boolean flag that indicates if the results matches the real condition. + */ + void ibz_kpoint(const ModuleSymmetry::Symmetry& symm, + bool use_symm, + std::string& skpt, + const UnitCell& ucell, + bool& match); + // LiuXh add 20180515 + + /** + * @brief Sets up the k-points after a volume change. + * + * This function sets up the k-points after a volume change in the system. + * It sets the Cartesian and direct k-vectors based on the new reciprocal and real space lattice vectors. + * + * @param nspin_in The number of spins. 1 for non-spin-polarized calculations and 2 for spin-polarized calculations. + * @param reciprocal_vec The new reciprocal lattice matrix. + * @param latvec The new real space lattice matrix. + * + * @return void + * + * @note The function first sets the number of spins (nspin) to the input value. + * @note If the direct k-vectors have been set (kd_done = true) and the Cartesian k-vectors have not (kc_done = + * false), the function calculates the Cartesian k-vectors by multiplying the direct k-vectors with the reciprocal + * lattice matrix. + * @note If the Cartesian k-vectors have been set (kc_done = true) and the direct k-vectors have not (kd_done = + * false), the function calculates the direct k-vectors by multiplying the Cartesian k-vectors with the transpose of + * the real space lattice matrix. + * @note The function also prints a table of the direct k-vectors and their weights. + * @note The function calls the print_klists function to print the k-points in both Cartesian and direct + * coordinates. + */ + void set_after_vc(const int& nspin, const ModuleBase::Matrix3& reciprocal_vec, const ModuleBase::Matrix3& latvec); + + /** + * @brief Gets the global index of a k-point. + * + * This function gets the global index of a k-point based on its local index and the process pool ID. + * The global index is used when the k-points are distributed among multiple process pools. + * + * @param ik The local index of the k-point. + * + * @return int Returns the global index of the k-point. + * + * @note The function calculates the global index by dividing the total number of k-points (nkstot) by the number of + * process pools (KPAR), and adding the remainder if the process pool ID (MY_POOL) is less than the remainder. + * @note The function is declared as inline for efficiency. + */ inline int getik_global(const int& ik) const; -private: + private: int nspin; bool kc_done; bool kd_done; - double koffset[3]; // used only in automatic k-points. - std::string k_kword; //LiuXh add 20180619 - int k_nkstot; //LiuXh add 20180619 - bool is_mp = false; //Monkhorst-Pack - - void renew( const int &kpoint_number ); + double koffset[3]; // used only in automatic k-points. + std::string k_kword; // LiuXh add 20180619 + int k_nkstot; // LiuXh add 20180619 + bool is_mp = false; // Monkhorst-Pack + + std::vector ibz2bz; // mohan added 2009-05-18 + + /** + * @brief Resize the k-point related vectors according to the new k-point number. + * + * This function resizes the vectors that store the k-point information, + * including the Cartesian and Direct coordinates of k-points, + * the weights of k-points, the index of k-points, and the number of plane waves for each k-point. + * + * @param kpoint_number The new number of k-points. + * + * @return void + * + * @note The memory recording lines are commented out. If you want to track the memory usage, + * you can uncomment these lines. + */ + void renew(const int& kpoint_number); // step 1 : generate kpoints - bool read_kpoints(const std::string &fn); // return 0: something wrong. + + /** + * @brief Reads the k-points from a file. + * + * This function reads the k-points from a file specified by the filename. + * It supports both Cartesian and Direct coordinates, and can handle different types of k-points, + * including Gamma, Monkhorst-Pack, and Line mode. It also supports automatic generation of k-points + * file if the file does not exist. + * + * @param fn The name of the file containing the k-points. + * + * @return bool Returns true if the k-points are successfully read from the file, + * false otherwise. + * + * @note It will generate a k-points file automatically + * according to the global variables GAMMA_ONLY_LOCAL and KSPACING. + * @note If the k-points type is neither Gamma nor Monkhorst-Pack, it will quit with a warning. + * @note If the k-points type is Line mode and the symmetry flag is 1, it will quit with a warning. + * @note If the number of k-points is greater than 100000, it will quit with a warning. + */ + bool read_kpoints(const std::string& fn); // return 0: something wrong. + + /** + * @brief Adds k-points linearly between special points. + * + * This function adds k-points linearly between special points in the Brillouin zone. + * The special points and the number of k-points between them are read from an input file. + * + * @param ifk The input file stream from which the special points and the number of k-points between them are read. + * @param kvec A vector to store the generated k-points. + * + * @return void + * + * @note The function first reads the number of special points (nks_special) and the number of k-points between them + * (nkl) from the input file. + * @note The function then recalculates the total number of k-points (nkstot) based on the number of k-points + * between the special points. + * @note The function generates the k-points by linearly interpolating between the special points. + * @note The function also assigns a segment ID to each k-point to distinguish different k-line segments. + * @note The function checks that the total number of generated k-points matches the calculated total number of + * k-points. + * @note The function checks that the size of the segment ID vector matches the total number of k-points. + */ void Linely_add_k_between(std::ifstream& ifk, std::vector>& kvec); - void Monkhorst_Pack(const int *nmp_in,const double *koffset_in,const int tipo); - double Monkhorst_Pack_formula( const int &k_type, const double &offset, - const int& n, const int &dim); + + /** + * @brief Generates k-points using the Monkhorst-Pack scheme. + * + * This function generates k-points in the reciprocal space using the Monkhorst-Pack scheme. + * + * @param nmp_in the number of k-points in each dimension. + * @param koffset_in the offset for the k-points in each dimension. + * @param k_type The type of k-point. 1 means without Gamma point, 0 means with Gamma. + * + * @return void + * + * @note The function assumes that the k-points are evenly distributed in the reciprocal space. + * @note The function sets the weight of each k-point to be equal, so that the total weight of all k-points is 1. + * @note The function sets the flag kd_done to true to indicate that the k-points have been generated. + */ + void Monkhorst_Pack(const int* nmp_in, const double* koffset_in, const int tipo); + + /** + * @brief Calculates the coordinate of a k-point using the Monkhorst-Pack scheme. + * + * This function calculates the coordinate of a k-point in the reciprocal space using the Monkhorst-Pack scheme. + * The Monkhorst-Pack scheme is a method for generating k-points in the Brillouin zone. + * + * @param k_type The type of k-point. 1 means without Gamma point, 0 means with Gamma. + * @param offset The offset for the k-point. + * @param n The index of the k-point in the current dimension. + * @param dim The total number of k-points in the current dimension. + * + * @return double Returns the coordinate of the k-point. + * + * @note The function assumes that the k-points are evenly distributed in the reciprocal space. + */ + double Monkhorst_Pack_formula(const int& k_type, const double& offset, const int& n, const int& dim); // step 2 : set both kvec and kved; normalize weight - void update_use_ibz( void ); - void set_both_kvec(const ModuleBase::Matrix3 &G,const ModuleBase::Matrix3 &R, std::string& skpt); - void normalize_wk( const int °spin ); + + /** + * @brief Updates the k-points to use the irreducible Brillouin zone (IBZ). + * + * This function updates the k-points to use the irreducible Brillouin zone (IBZ) instead of the full Brillouin + * zone. + * + * @return void + * + * @note This function should only be called by the master process (MY_RANK == 0). + * @note This function assumes that the number of k-points in the IBZ (nkstot_ibz) is greater than 0. + * @note This function updates the total number of k-points (nkstot) to be the number of k-points in the IBZ. + * @note This function resizes the vector of k-points (kvec_d) and updates its values to be the k-points in the IBZ. + * @note This function also updates the weights of the k-points (wk) to be the weights in the IBZ. + * @note After this function is called, the flag kd_done is set to true to indicate that the k-points have been + * updated, and the flag kc_done is set to false to indicate that the Cartesian coordinates of the k-points need to + * be recalculated. + */ + void update_use_ibz(void); + + /** + * @brief Sets both the direct and Cartesian k-vectors. + * + * This function sets both the direct and Cartesian k-vectors based on the input parameters. + * It also checks the k-point type and sets the corresponding flags. + * + * @param G The reciprocal lattice matrix. + * @param R The real space lattice matrix. + * @param skpt A string to store the k-point table. + * + * @return void + * + * @note If the k-point type is neither "Cartesian" nor "Direct", an error message will be printed. + * @note The function sets the flags kd_done and kc_done to indicate whether the direct and Cartesian k-vectors have + * been set, respectively. + * @note The function also prints a table of the direct k-vectors and their weights. + * @note If the function is called by the master process (MY_RANK == 0), the k-point table is also stored in the + * string skpt. + */ + void set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt); + + /** + * @brief Normalizes the weights of the k-points. + * + * This function normalizes the weights of the k-points so that their sum is equal to the degeneracy of spin + * (degspin). + * + * @param degspin The degeneracy of spin. This is 1 for non-spin-polarized calculations and 2 for spin-polarized + * calculations. + * + * @return void + * + * @note This function should only be called by the master process (MY_RANK == 0). + * @note The function assumes that the sum of the weights of the k-points is greater than 0. + * @note The function first normalizes the weights so that their sum is 1, and then scales them by the degeneracy of + * spin. + */ + void normalize_wk(const int& degspin); // step 3 : mpi kpoints information. + + /** + * @brief Distributes k-points among MPI processes. + * + * This function distributes the k-points among the MPI processes. Each process gets a subset of the k-points to + * work on. The function also broadcasts various variables related to the k-points to all processes. + * + * @return void + * + * @note This function is only compiled and used if MPI is enabled. + * @note The function assumes that the number of k-points (nkstot) is greater than 0. + * @note The function broadcasts the flags kc_done and kd_done, the number of spins (nspin), the total number of + * k-points (nkstot), the full number of k-points (nkstot_full), the Monkhorst-Pack grid (nmp), the k-point offsets + * (koffset), and the segment IDs of the k-points (kl_segids). + * @note The function also broadcasts the indices of the k-points (isk), their weights (wk), and their Cartesian and + * direct coordinates (kvec_c and kvec_d). + * @note If a process has no k-points to work on, the function will quit with an error message. + */ void mpi_k(); // step 4 : *2 or *4 kpoints. - // *2 for LSDA - // *4 for non-collinear + + /** + * @brief Sets up the k-points for spin-up and spin-down calculations. + * + * This function sets up the k-points for spin-up and spin-down calculations. + * If the calculation is spin-polarized (nspin = 2), the number of k-points is doubled. + * The first half of the k-points correspond to spin-up, and the second half correspond to spin-down. + * 2 for LSDA + * 4 for non-collinear + * + * @return void + * + * @note For non-spin-polarized calculations (nspin = 1 or 4), the function simply sets the spin index of all + * k-points to 0. + * @note For spin-polarized calculations (nspin = 2), the function duplicates the k-points and their weights, + * sets the spin index of the first half of the k-points to 0 (spin-up), and the spin index of the second half + * to 1 (spin-down). + * @note The function also doubles the total number of k-points (nks and nkstot) for spin-polarized calculations. + * @note The function prints the total number of k-points for spin-polarized calculations. + */ void set_kup_and_kdw(); // step 5 // print k lists. - void print_klists(std::ofstream &fn); + + /** + * @brief Prints the k-points in both Cartesian and direct coordinates. + * + * This function prints the k-points in both Cartesian and direct coordinates to the output file stream. + * The output includes the index, x, y, and z coordinates, and the weight of each k-point. + * + * @param ofs The output file stream to which the k-points are printed. + * + * @return void + * + * @note The function first checks if the total number of k-points (nkstot) is less than the number of k-points for + * the current spin (nks). If so, it prints an error message and quits. + * @note The function prints the k-points in a table format, with separate tables for Cartesian and direct + * coordinates. + * @note The function uses the FmtCore::format function to format the output. + */ + void print_klists(std::ofstream& fn); }; -inline int K_Vectors:: getik_global(const int& ik) const +inline int K_Vectors::getik_global(const int& ik) const { int nkp = this->nkstot / GlobalV::KPAR; int rem = this->nkstot % GlobalV::KPAR; - if(GlobalV::MY_POOL < rem) + if (GlobalV::MY_POOL < rem) { - return GlobalV::MY_POOL*nkp + GlobalV::MY_POOL + ik; + return GlobalV::MY_POOL * nkp + GlobalV::MY_POOL + ik; } else { - return GlobalV::MY_POOL*nkp + rem + ik; + return GlobalV::MY_POOL * nkp + rem + ik; } }