diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 322a9ffca9..6710c2afe6 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -192,6 +192,7 @@ OBJS_CELL=atom_pseudo.o\ print_cell.o\ setup_nonlocal.o\ klist.o\ + k_vector_utils.o\ cell_index.o\ check_atomic_stru.o\ update_cell.o\ diff --git a/source/module_cell/CMakeLists.txt b/source/module_cell/CMakeLists.txt index e340ddd272..95d997ecfd 100644 --- a/source/module_cell/CMakeLists.txt +++ b/source/module_cell/CMakeLists.txt @@ -25,6 +25,7 @@ add_library( read_stru.cpp print_cell.cpp read_atom_species.cpp + k_vector_utils.cpp ) if(ENABLE_COVERAGE) diff --git a/source/module_cell/k_vector_utils.cpp b/source/module_cell/k_vector_utils.cpp new file mode 100644 index 0000000000..a9dd7e7f56 --- /dev/null +++ b/source/module_cell/k_vector_utils.cpp @@ -0,0 +1,808 @@ +// +// Created by rhx on 25-6-3. +// +#include "k_vector_utils.h" + +#include "klist.h" +#include "module_base/global_variable.h" +#include "module_base/matrix3.h" + +#include "module_base/formatter.h" +#include "module_base/parallel_common.h" +#include "module_base/parallel_reduce.h" +#include "module_parameter/parameter.h" + +namespace KVectorUtils +{ +void kvec_d2c(K_Vectors& kv, const ModuleBase::Matrix3& reciprocal_vec) +{ + // throw std::runtime_error("k_vec_d2c: This function is not implemented in the new codebase. Please use the new + // implementation."); + if (kv.kvec_d.size() != kv.kvec_c.size()) + { + // ModuleBase::WARNING_QUIT("k_vec_d2c", "Size of Cartesian and Direct K vectors mismatch. "); + kv.kvec_c.resize(kv.kvec_d.size()); + } + int nks = kv.kvec_d.size(); // always convert all k vectors + + for (int i = 0; i < nks; i++) + { + // wrong!! kvec_c[i] = G * kvec_d[i]; + // mohan fixed bug 2010-1-10 + if (std::abs(kv.kvec_d[i].x) < 1.0e-10) + { + kv.kvec_d[i].x = 0.0; + } + if (std::abs(kv.kvec_d[i].y) < 1.0e-10) + { + kv.kvec_d[i].y = 0.0; + } + if (std::abs(kv.kvec_d[i].z) < 1.0e-10) + { + kv.kvec_d[i].z = 0.0; + } + + kv.kvec_c[i] = kv.kvec_d[i] * reciprocal_vec; + + // mohan add2012-06-10 + if (std::abs(kv.kvec_c[i].x) < 1.0e-10) + { + kv.kvec_c[i].x = 0.0; + } + if (std::abs(kv.kvec_c[i].y) < 1.0e-10) + { + kv.kvec_c[i].y = 0.0; + } + if (std::abs(kv.kvec_c[i].z) < 1.0e-10) + { + kv.kvec_c[i].z = 0.0; + } + } +} +void kvec_c2d(K_Vectors& kv, const ModuleBase::Matrix3& latvec) +{ + if (kv.kvec_d.size() != kv.kvec_c.size()) + { + kv.kvec_d.resize(kv.kvec_c.size()); + } + int nks = kv.kvec_d.size(); // always convert all k vectors + + ModuleBase::Matrix3 RT = latvec.Transpose(); + 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 + kv.kvec_d[i] = kv.kvec_c[i] * RT; + } +} + +void set_both_kvec(K_Vectors& kv, const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt) +{ + if (true) // Originally GlobalV::FINAL_SCF, but we don't have this variable in the new code. + { + if (kv.get_k_nkstot() == 0) + { + kv.kd_done = true; + kv.kc_done = false; + } + else + { + if (kv.get_k_kword() == "Cartesian" || kv.get_k_kword() == "C") + { + kv.kc_done = true; + kv.kd_done = false; + } + else if (kv.get_k_kword() == "Direct" || kv.get_k_kword() == "D") + { + kv.kd_done = true; + kv.kc_done = false; + } + else + { + GlobalV::ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; + } + } + } + + // set cartesian k vectors. + if (!kv.kc_done && kv.kd_done) + { + KVectorUtils::kvec_d2c(kv, G); + kv.kc_done = true; + } + + // set direct k vectors + else if (kv.kc_done && !kv.kd_done) + { + KVectorUtils::kvec_c2d(kv, R); + kv.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; i < kv.get_nkstot(); i++) + { + table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n", + i + 1, + kv.kvec_d[i].x, + kv.kvec_d[i].y, + kv.kvec_d[i].z, + kv.wk[i]); + } + GlobalV::ofs_running << table << std::endl; + if (GlobalV::MY_RANK == 0) + { + std::stringstream ss; + ss << " " << std::setw(40) << "nkstot now" + << " = " << kv.get_nkstot() << std::endl; + ss << table << std::endl; + skpt = ss.str(); + } + return; +} + +void set_after_vc(K_Vectors& kv, const int& nspin_in, const ModuleBase::Matrix3& reciprocal_vec) +{ + GlobalV::ofs_running << "\n SETUP K-POINTS" << std::endl; + // kv.nspin = nspin_in; + kv.set_nspin(nspin_in); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nspin", kv.get_nspin()); + + // set cartesian k vectors. + KVectorUtils::kvec_d2c(kv, reciprocal_vec); + + 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; i < kv.get_nks(); i++) + { + table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n", + i + 1, + kv.kvec_d[i].x, + kv.kvec_d[i].y, + kv.kvec_d[i].z, + kv.wk[i]); + } + GlobalV::ofs_running << table << std::endl; + + kv.kd_done = true; + kv.kc_done = true; + + print_klists(kv, GlobalV::ofs_running); +} + +void print_klists(const K_Vectors& kv, std::ofstream& ofs) +{ + ModuleBase::TITLE("KVectorUtils", "print_klists"); + int nks = kv.get_nks(); + int nkstot = kv.get_nkstot(); + + if (nkstot < nks) + { + std::cout << "\n nkstot=" << nkstot; + std::cout << "\n nks=" << 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; i < nks; i++) + { + table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n", + i + 1, + kv.kvec_c[i].x, + kv.kvec_c[i].y, + kv.kvec_c[i].z, + kv.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; i < nks; i++) + { + table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n", + i + 1, + kv.kvec_d[i].x, + kv.kvec_d[i].y, + kv.kvec_d[i].z, + kv.wk[i]); + } + GlobalV::ofs_running << "\n" << table << std::endl; + return; +} + +#ifdef __MPI +void kvec_mpi_k(K_Vectors& kv) +{ + ModuleBase::TITLE("KVectorUtils", "kvec_mpi_k"); + + Parallel_Common::bcast_bool(kv.kc_done); + + Parallel_Common::bcast_bool(kv.kd_done); + + Parallel_Common::bcast_int(kv.nspin); + + Parallel_Common::bcast_int(kv.nkstot); + + Parallel_Common::bcast_int(kv.nkstot_full); + + Parallel_Common::bcast_int(kv.nmp, 3); + + kv.kl_segids.resize(kv.nkstot); + Parallel_Common::bcast_int(kv.kl_segids.data(), kv.nkstot); + + Parallel_Common::bcast_double(kv.koffset, 3); + + kv.nks = kv.para_k.nks_pool[GlobalV::MY_POOL]; + + GlobalV::ofs_running << std::endl; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "k-point number in this process", kv.nks); + int nks_minimum = kv.nks; + + Parallel_Reduce::gather_min_int_all(GlobalV::NPROC, nks_minimum); + + if (nks_minimum == 0) + { + 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); + } + + std::vector isk_aux(kv.nkstot); + std::vector wk_aux(kv.nkstot); + std::vector kvec_c_aux(kv.nkstot * 3); + std::vector kvec_d_aux(kv.nkstot * 3); + + // collect and process in rank 0 + if (GlobalV::MY_RANK == 0) + { + for (int ik = 0; ik < kv.nkstot; ik++) + { + isk_aux[ik] = kv.isk[ik]; + wk_aux[ik] = kv.wk[ik]; + kvec_c_aux[3 * ik] = kv.kvec_c[ik].x; + kvec_c_aux[3 * ik + 1] = kv.kvec_c[ik].y; + kvec_c_aux[3 * ik + 2] = kv.kvec_c[ik].z; + kvec_d_aux[3 * ik] = kv.kvec_d[ik].x; + kvec_d_aux[3 * ik + 1] = kv.kvec_d[ik].y; + kvec_d_aux[3 * ik + 2] = kv.kvec_d[ik].z; + } + } + + // broadcast k point data to all processors + Parallel_Common::bcast_int(isk_aux.data(), kv.nkstot); + + Parallel_Common::bcast_double(wk_aux.data(), kv.nkstot); + Parallel_Common::bcast_double(kvec_c_aux.data(), kv.nkstot * 3); + Parallel_Common::bcast_double(kvec_d_aux.data(), kv.nkstot * 3); + + // process k point data in each processor + kv.renew(kv.nks * kv.nspin); + + // distribute + int k_index = 0; + + for (int i = 0; i < kv.nks; i++) + { + // 3 is because each k point has three value:kx, ky, kz + k_index = i + kv.para_k.startk_pool[GlobalV::MY_POOL]; + kv.kvec_c[i].x = kvec_c_aux[k_index * 3]; + kv.kvec_c[i].y = kvec_c_aux[k_index * 3 + 1]; + kv.kvec_c[i].z = kvec_c_aux[k_index * 3 + 2]; + kv.kvec_d[i].x = kvec_d_aux[k_index * 3]; + kv.kvec_d[i].y = kvec_d_aux[k_index * 3 + 1]; + kv.kvec_d[i].z = kvec_d_aux[k_index * 3 + 2]; + kv.wk[i] = wk_aux[k_index]; + kv.isk[i] = isk_aux[k_index]; + } + +#ifdef __EXX + if (ModuleSymmetry::Symmetry::symm_flag == 1) + { // bcast kstars + kv.kstars.resize(kv.nkstot); + for (int ikibz = 0; ikibz < kv.nkstot; ++ikibz) + { + int starsize = kv.kstars[ikibz].size(); + Parallel_Common::bcast_int(starsize); + GlobalV::ofs_running << "starsize: " << starsize << std::endl; + auto ks = kv.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); + GlobalV::ofs_running << "isym: " << isym << " ks_vec: " << ks_vec.x << " " << ks_vec.y << " " + << ks_vec.z << std::endl; + if (GlobalV::MY_RANK != 0) + { + kv.kstars[ikibz].insert(std::make_pair(isym, ks_vec)); + } + } + } + } +#endif +} // END SUBROUTINE +#endif + + +void kvec_ibz_kpoint(K_Vectors& kv, + const ModuleSymmetry::Symmetry& symm, + bool use_symm, + std::string& skpt, + const UnitCell& ucell, + bool& match) +{ + if (GlobalV::MY_RANK != 0) + { + return; + } + ModuleBase::TITLE("K_Vectors", "ibz_kpoint"); + + // k-lattice: "pricell" of reciprocal space + // CAUTION: should fit into all k-input method, not only MP !!! + // the basis vector of reciprocal lattice: recip_vec1, recip_vec2, recip_vec3 + ModuleBase::Vector3 recip_vec1(ucell.G.e11, ucell.G.e12, ucell.G.e13); + ModuleBase::Vector3 recip_vec2(ucell.G.e21, ucell.G.e22, ucell.G.e23); + ModuleBase::Vector3 recip_vec3(ucell.G.e31, ucell.G.e32, ucell.G.e33); + ModuleBase::Vector3 k_vec1, k_vec2, k_vec3; + ModuleBase::Matrix3 k_vec; + if (kv.get_is_mp()) + { + k_vec1 = ModuleBase::Vector3(recip_vec1.x / kv.nmp[0], recip_vec1.y / kv.nmp[0], recip_vec1.z / kv.nmp[0]); + k_vec2 = ModuleBase::Vector3(recip_vec2.x / kv.nmp[1], recip_vec2.y / kv.nmp[1], recip_vec2.z / kv.nmp[1]); + k_vec3 = ModuleBase::Vector3(recip_vec3.x / kv.nmp[2], recip_vec3.y / kv.nmp[2], recip_vec3.z / kv.nmp[2]); + k_vec = ModuleBase::Matrix3(k_vec1.x, + k_vec1.y, + k_vec1.z, + k_vec2.x, + k_vec2.y, + k_vec2.z, + k_vec3.x, + k_vec3.y, + k_vec3.z); + } + + //=============================================== + // search in all space group operations + // 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 ind(1, 0, 0, 0, 1, 0, 0, 0, 1); + + int nrotkm = 0; + if (use_symm) + { + // bravais type of reciprocal lattice and k-lattice + + double recip_vec_const[6]; + double recip_vec0_const[6]; + double k_vec_const[6]; + double k_vec0_const[6]; + int recip_brav_type = 15; + int k_brav_type = 15; + std::string recip_brav_name; + std::string k_brav_name; + ModuleBase::Vector3 k_vec01 = k_vec1, k_vec02 = k_vec2, k_vec03 = k_vec3; + + // it's not necessary to calculate gb01, gb02, gb03, + // because they are only used as a vector, no need to be assigned values + + // determine the Bravais type and related parameters of the lattice + symm.lattice_type(recip_vec1, + recip_vec2, + recip_vec3, + recip_vec1, + recip_vec2, + recip_vec3, + recip_vec_const, + recip_vec0_const, + recip_brav_type, + recip_brav_name, + ucell.atoms, + false, + nullptr); + GlobalV::ofs_running << "\n For reciprocal-space lattice:" << std::endl; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS TYPE", recip_brav_type); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS LATTICE NAME", recip_brav_name); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "ibrav", recip_brav_type); + + // the map of bravis lattice from real to reciprocal space + // for example, 3(fcc) in real space matches 2(bcc) in reciprocal space + std::vector ibrav_a2b{1, 3, 2, 4, 5, 6, 7, 8, 10, 9, 11, 12, 13, 14}; + // check if the reciprocal lattice is compatible with the real space lattice + auto ibrav_match = [&](int ibrav_b) -> bool { + const int& ibrav_a = symm.real_brav; + if (ibrav_a < 1 || ibrav_a > 14) + { + return false; + } + return (ibrav_b == ibrav_a2b[ibrav_a - 1]); + }; + if (!ibrav_match(recip_brav_type)) // if not match, exit and return + { + 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: " << recip_brav_name << std::endl; + GlobalV::ofs_running << "(which should be " << ibrav_a2b[symm.real_brav - 1] << ")." << std::endl; + match = false; + return; + } + + // if match, continue + if (kv.get_is_mp()) + { + symm.lattice_type(k_vec1, + k_vec2, + k_vec3, + k_vec01, + k_vec02, + k_vec03, + k_vec_const, + k_vec0_const, + k_brav_type, + k_brav_name, + ucell.atoms, + false, + nullptr); + GlobalV::ofs_running << "\n For k vectors:" << std::endl; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS TYPE", k_brav_type); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS LATTICE NAME", k_brav_name); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "ibrav", k_brav_type); + } + // point-group analysis of reciprocal lattice + ModuleBase::Matrix3 bsymop[48]; + int bnop = 0; + // search again + symm.lattice_type(recip_vec1, + recip_vec2, + recip_vec3, + recip_vec1, + recip_vec2, + recip_vec3, + recip_vec_const, + recip_vec0_const, + recip_brav_type, + recip_brav_name, + ucell.atoms, + false, + nullptr); + ModuleBase::Matrix3 b_optlat_new(recip_vec1.x, + recip_vec1.y, + recip_vec1.z, + recip_vec2.x, + recip_vec2.y, + recip_vec2.z, + recip_vec3.x, + recip_vec3.y, + recip_vec3.z); + // set the crystal point-group symmetry operation + symm.setgroup(bsymop, bnop, recip_brav_type); + // transform the above symmetric operation matrices between different coordinate + 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.e32, b.e32) && symm.equal(a.e33, b.e33)); + }; + for (int i = 0; i < symm.nrotk; ++i) + { + match = false; + for (int j = 0; j < bnop; ++j) + { + if (matequal(symm.kgmatrix[i], bsymop[j])) + { + match = true; + break; + } + } + if (!match) + { + return; + } + } + nrotkm = symm.nrotk; // change if inv not included + for (int i = 0; i < nrotkm; ++i) + { + if (symm.kgmatrix[i] == inv) + { + include_inv = true; + } + kgmatrix[i] = symm.kgmatrix[i]; + } + + if (!include_inv) + { + for (int i = 0; i < symm.nrotk; ++i) + { + kgmatrix[i + symm.nrotk] = inv * symm.kgmatrix[i]; + } + nrotkm = 2 * symm.nrotk; + } + } + else if (kv.get_is_mp()) // only include for mp grid + { + nrotkm = 2; + kgmatrix[0] = ind; + kgmatrix[1] = inv; + } + else + { + return; + } + + // convert kgmatrix to k-lattice + ModuleBase::Matrix3* kkmatrix = new ModuleBase::Matrix3[nrotkm]; + if (kv.get_is_mp()) + { + symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, k_vec); + } + // direct coordinates of k-points in k-lattice + std::vector> kvec_d_k(kv.get_nkstot()); + if (kv.get_is_mp()) + { + for (int i = 0; i < kv.get_nkstot(); ++i) + { + kvec_d_k[i] = kv.kvec_d[i] * ucell.G * k_vec.Inverse(); + } + } + + // use operation : kgmatrix to find + // the new set kvec_d : ir_kpt + int nkstot_ibz = 0; + + assert(kv.get_nkstot() > 0); + std::vector> kvec_d_ibz(kv.get_nkstot()); + std::vector wk_ibz(kv.get_nkstot()); // ibz kpoint wk ,weight of k points + std::vector ibz2bz(kv.get_nkstot()); + + // nkstot is the total input k-points number. + const double weight = 1.0 / static_cast(kv.get_nkstot()); + + ModuleBase::Vector3 kvec_rot; + ModuleBase::Vector3 kvec_rot_k; + + // 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; + // 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; + } + return; + }; + // for output in kpoints file + int ibz_index[kv.get_nkstot()]; + // search in all k-poins. + for (int i = 0; i < kv.get_nkstot(); ++i) + { + // restrict to [0, 1) + restrict_kpt(kv.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; + bool already_exist = false; + int exist_number = -1; + // search over all symmetry operations + for (int j = 0; j < nrotkm; ++j) + { + if (!already_exist) + { + // 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 = kv.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 (kv.get_is_mp()) + { + kvec_rot_k = kvec_d_k[i] * kkmatrix[j]; // k-lattice rotation + kvec_rot_k = kvec_rot_k * k_vec * 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 * k_vec.Inverse(); // convert back to k-latice + } + 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)) + { + already_exist = true; + // find another ibz k point, + // but is already in the ibz_kpoint list. + // so the weight need to +1; + wk_ibz[k] += weight; + exist_number = k; + break; + } + } + } // 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. + kvec_d_ibz[nkstot_ibz] = kv.kvec_d[i]; + // output in kpoints file + ibz_index[i] = nkstot_ibz; + + // the weight should be averged k-point weight. + wk_ibz[nkstot_ibz] = weight; + + // ibz2bz records the index of origin k points. + ibz2bz[nkstot_ibz] = i; + ++nkstot_ibz; + } + 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; + + double kmol_new = kv.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] = kv.kvec_d[i]; + } + } + // BLOCK_HERE("check k point"); + } + + delete[] kkmatrix; + +#ifdef __EXX + // setup kstars according to the final (max-norm) kvec_d_ibz + kv.kstars.resize(nkstot_ibz); + if (ModuleSymmetry::Symmetry::symm_flag == 1) + { + for (int i = 0; i < kv.get_nkstot(); ++i) + { + int exist_number = -1; + int isym = 0; + for (int j = 0; j < nrotkm; ++j) + { + kvec_rot = kv.kvec_d[i] * kgmatrix[j]; + restrict_kpt(kvec_rot); + 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; + } + } + kv.kstars[exist_number].insert(std::make_pair(isym, kv.kvec_d[i])); + } + } +#endif + + // output in kpoints file + std::stringstream ss; + ss << " " << std::setw(40) << "nkstot" + << " = " << kv.get_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 < kv.get_nkstot(); ++i) + { + table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8d%12.8f%12.8f%12.8f\n", + i + 1, + kv.kvec_d[i].x, + kv.kvec_d[i].y, + kv.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; + skpt = ss.str(); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nkstot_ibz", 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; + + // resize the kpoint container according to nkstot_ibz + if (use_symm || kv.get_is_mp()) + { + kv.update_use_ibz(nkstot_ibz, kvec_d_ibz, wk_ibz); + } + + return; +} +} // namespace KVectorUtils diff --git a/source/module_cell/k_vector_utils.h b/source/module_cell/k_vector_utils.h new file mode 100644 index 0000000000..c492a49ad0 --- /dev/null +++ b/source/module_cell/k_vector_utils.h @@ -0,0 +1,128 @@ +// +// Created by rhx on 25-6-3. +// + +#ifndef K_VECTOR_UTILS_H +#define K_VECTOR_UTILS_H + +#include "module_base/global_variable.h" +#include "module_base/matrix3.h" +#include "module_cell/unitcell.h" + +class K_Vectors; + +namespace KVectorUtils +{ +void kvec_d2c(K_Vectors& kv, const ModuleBase::Matrix3& reciprocal_vec); + +void kvec_c2d(K_Vectors& kv, const ModuleBase::Matrix3& latvec); + +/** + * @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 kv The K_Vectors object containing the k-point information. + * @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(K_Vectors& kv, const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt); + +/** + * @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 kv The K_Vectors object containing the k-point information. + * @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. + * + * @return void + * + * @note The function first sets the number of spins (nspin) to the input value. + * @note The direct k-vectors have been set (kd_done = true) but the Cartesian k-vectors have not (kc_done = + * false) after a volume change. The function calculates the Cartesian k-vectors by multiplying the direct k-vectors + * with the reciprocal 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(K_Vectors& kv, const int& nspin, const ModuleBase::Matrix3& G); + +/** + * @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(const K_Vectors& kv, std::ofstream& ofs); + +// 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. + * + * @param kv The K_Vectors object containing the k-point information. + * + * @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. + */ +#ifdef __MPI +void kvec_mpi_k(K_Vectors& kv); +#endif // __MPI + +/** + * @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 kvec_ibz_kpoint(K_Vectors& kv, + const ModuleSymmetry::Symmetry& symm, + bool use_symm, + std::string& skpt, + const UnitCell& ucell, + bool& match); +} // namespace KVectorUtils + +#endif // K_VECTOR_UTILS_H diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp index 51ecdc0008..da5f9f16a0 100644 --- a/source/module_cell/klist.cpp +++ b/source/module_cell/klist.cpp @@ -1,5 +1,6 @@ #include "klist.h" +#include "k_vector_utils.h" #include "module_base/formatter.h" #include "module_base/memory.h" #include "module_base/parallel_common.h" @@ -94,7 +95,7 @@ void K_Vectors::set(const UnitCell& ucell, { bool match = true; // calculate kpoints in IBZ and reduce kpoints according to symmetry - this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, ucell, match); + KVectorUtils::kvec_ibz_kpoint(*this, symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, ucell, match); #ifdef __MPI Parallel_Common::bcast_bool(match); #endif @@ -109,7 +110,7 @@ void K_Vectors::set(const UnitCell& ucell, std::cout << "Automatically set symmetry to 0 and continue ..." << std::endl; ModuleSymmetry::Symmetry::symm_flag = 0; match = true; - this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, ucell, match); + KVectorUtils::kvec_ibz_kpoint(*this, symm, ModuleSymmetry::Symmetry::symm_flag, skpt1, ucell, match); } else { ModuleBase::WARNING_QUIT("K_Vectors::ibz_kpoint", "Possible solutions: \n \ @@ -125,7 +126,8 @@ 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); +// this->set_both_kvec(reciprocal_vec, latvec, skpt2); + KVectorUtils::set_both_kvec(*this, reciprocal_vec, latvec, skpt2); if (GlobalV::MY_RANK == 0) { @@ -152,7 +154,7 @@ void K_Vectors::set(const UnitCell& ucell, nspin_in); // assign k points to several process pools #ifdef __MPI // distribute K point data to the corresponding process - this->mpi_k(); // 2008-4-29 + KVectorUtils::kvec_mpi_k(*this); #endif // set the k vectors for the up and down spin @@ -161,7 +163,7 @@ void K_Vectors::set(const UnitCell& ucell, // get ik2iktot this->cal_ik_global(); - this->print_klists(ofs); + KVectorUtils::print_klists(*this, ofs); // std::cout << " NUMBER OF K-POINTS : " << nkstot << std::endl; @@ -546,564 +548,6 @@ void K_Vectors::update_use_ibz(const int& nkstot_ibz, return; } -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; - } - ModuleBase::TITLE("K_Vectors", "ibz_kpoint"); - - // k-lattice: "pricell" of reciprocal space - // CAUTION: should fit into all k-input method, not only MP !!! - // the basis vector of reciprocal lattice: recip_vec1, recip_vec2, recip_vec3 - ModuleBase::Vector3 recip_vec1(ucell.G.e11, ucell.G.e12, ucell.G.e13); - ModuleBase::Vector3 recip_vec2(ucell.G.e21, ucell.G.e22, ucell.G.e23); - ModuleBase::Vector3 recip_vec3(ucell.G.e31, ucell.G.e32, ucell.G.e33); - ModuleBase::Vector3 k_vec1, k_vec2, k_vec3; - ModuleBase::Matrix3 k_vec; - if (this->is_mp) - { - k_vec1 = ModuleBase::Vector3(recip_vec1.x / nmp[0], recip_vec1.y / nmp[0], recip_vec1.z / nmp[0]); - k_vec2 = ModuleBase::Vector3(recip_vec2.x / nmp[1], recip_vec2.y / nmp[1], recip_vec2.z / nmp[1]); - k_vec3 = ModuleBase::Vector3(recip_vec3.x / nmp[2], recip_vec3.y / nmp[2], recip_vec3.z / nmp[2]); - k_vec = ModuleBase::Matrix3(k_vec1.x, - k_vec1.y, - k_vec1.z, - k_vec2.x, - k_vec2.y, - k_vec2.z, - k_vec3.x, - k_vec3.y, - k_vec3.z); - } - - //=============================================== - // search in all space group operations - // 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 ind(1, 0, 0, 0, 1, 0, 0, 0, 1); - - int nrotkm = 0; - if (use_symm) - { - // bravais type of reciprocal lattice and k-lattice - - double recip_vec_const[6]; - double recip_vec0_const[6]; - double k_vec_const[6]; - double k_vec0_const[6]; - int recip_brav_type = 15; - int k_brav_type = 15; - std::string recip_brav_name; - std::string k_brav_name; - ModuleBase::Vector3 k_vec01 = k_vec1, k_vec02 = k_vec2, k_vec03 = k_vec3; - - // it's not necessary to calculate gb01, gb02, gb03, - // because they are only used as a vector, no need to be assigned values - - // determine the Bravais type and related parameters of the lattice - symm.lattice_type(recip_vec1, - recip_vec2, - recip_vec3, - recip_vec1, - recip_vec2, - recip_vec3, - recip_vec_const, - recip_vec0_const, - recip_brav_type, - recip_brav_name, - ucell.atoms, - false, - nullptr); - GlobalV::ofs_running << "\n For reciprocal-space lattice:" << std::endl; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS TYPE", recip_brav_type); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS LATTICE NAME", recip_brav_name); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "ibrav", recip_brav_type); - - // the map of bravis lattice from real to reciprocal space - // for example, 3(fcc) in real space matches 2(bcc) in reciprocal space - std::vector ibrav_a2b{1, 3, 2, 4, 5, 6, 7, 8, 10, 9, 11, 12, 13, 14}; - // check if the reciprocal lattice is compatible with the real space lattice - auto ibrav_match = [&](int ibrav_b) -> bool { - const int& ibrav_a = symm.real_brav; - if (ibrav_a < 1 || ibrav_a > 14) { - return false; - } - return (ibrav_b == ibrav_a2b[ibrav_a - 1]); - }; - if (!ibrav_match(recip_brav_type)) // if not match, exit and return - { - 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: " << recip_brav_name << std::endl; - GlobalV::ofs_running << "(which should be " << ibrav_a2b[symm.real_brav - 1] << ")." << std::endl; - match = false; - return; - } - - // if match, continue - if (this->is_mp) - { - symm.lattice_type(k_vec1, - k_vec2, - k_vec3, - k_vec01, - k_vec02, - k_vec03, - k_vec_const, - k_vec0_const, - k_brav_type, - k_brav_name, - ucell.atoms, - false, - nullptr); - GlobalV::ofs_running << "\n For k vectors:" << std::endl; - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS TYPE", k_brav_type); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "BRAVAIS LATTICE NAME", k_brav_name); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "ibrav", k_brav_type); - } - // point-group analysis of reciprocal lattice - ModuleBase::Matrix3 bsymop[48]; - int bnop = 0; - // search again - symm.lattice_type(recip_vec1, - recip_vec2, - recip_vec3, - recip_vec1, - recip_vec2, - recip_vec3, - recip_vec_const, - recip_vec0_const, - recip_brav_type, - recip_brav_name, - ucell.atoms, - false, - nullptr); - ModuleBase::Matrix3 b_optlat_new(recip_vec1.x, - recip_vec1.y, - recip_vec1.z, - recip_vec2.x, - recip_vec2.y, - recip_vec2.z, - recip_vec3.x, - recip_vec3.y, - recip_vec3.z); - // set the crystal point-group symmetry operation - symm.setgroup(bsymop, bnop, recip_brav_type); - // transform the above symmetric operation matrices between different coordinate - 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.e32, b.e32) && symm.equal(a.e33, b.e33)); - }; - for (int i = 0; i < symm.nrotk; ++i) - { - match = false; - for (int j = 0; j < bnop; ++j) - { - if (matequal(symm.kgmatrix[i], bsymop[j])) - { - match = true; - break; - } - } - if (!match) { - return; - } - } - nrotkm = symm.nrotk; // change if inv not included - for (int i = 0; i < nrotkm; ++i) - { - if (symm.kgmatrix[i] == inv) - { - include_inv = true; - } - kgmatrix[i] = symm.kgmatrix[i]; - } - - if (!include_inv) - { - for (int i = 0; i < symm.nrotk; ++i) - { - kgmatrix[i + symm.nrotk] = inv * symm.kgmatrix[i]; - } - nrotkm = 2 * symm.nrotk; - } - } - else if (is_mp) // only include for mp grid - { - nrotkm = 2; - kgmatrix[0] = ind; - kgmatrix[1] = inv; - } - else - { - return; - } - - // convert kgmatrix to k-lattice - ModuleBase::Matrix3* kkmatrix = new ModuleBase::Matrix3[nrotkm]; - if (this->is_mp) { - symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, k_vec); - } - // 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 * k_vec.Inverse(); - } - } - - // use operation : kgmatrix to find - // the new set kvec_d : ir_kpt - int nkstot_ibz = 0; - - assert(nkstot > 0); - std::vector> kvec_d_ibz(this->nkstot); - std::vector wk_ibz(this->nkstot); // ibz kpoint wk ,weight of k points - std::vector ibz2bz(this->nkstot); - - // 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) { - // 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; - // 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; - } - return; - }; - // for output in kpoints file - int ibz_index[nkstot]; - // search in all k-poins. - for (int i = 0; i < nkstot; ++i) - { - // 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; - bool already_exist = false; - int exist_number = -1; - // search over all symmetry operations - for (int j = 0; j < nrotkm; ++j) - { - if (!already_exist) - { - // 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. - 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 * k_vec * 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 * k_vec.Inverse(); // convert back to k-latice - } - 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)) - { - already_exist = true; - // find another ibz k point, - // but is already in the ibz_kpoint list. - // so the weight need to +1; - wk_ibz[k] += weight; - exist_number = k; - break; - } - } - } // 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. - kvec_d_ibz[nkstot_ibz] = kvec_d[i]; - // output in kpoints file - ibz_index[i] = nkstot_ibz; - - // the weight should be averged k-point weight. - wk_ibz[nkstot_ibz] = weight; - - // ibz2bz records the index of origin k points. - ibz2bz[nkstot_ibz] = i; - ++nkstot_ibz; - } - 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; - - 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]; - } - } - // BLOCK_HERE("check k point"); - } - - delete[] kkmatrix; - -#ifdef __EXX - // setup kstars according to the final (max-norm) kvec_d_ibz - this->kstars.resize(nkstot_ibz); - if (ModuleSymmetry::Symmetry::symm_flag == 1) - { - 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]; - restrict_kpt(kvec_rot); - 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, kvec_d[i])); - } - } -#endif - - // output in kpoints file - 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, - this->kvec_d[i].x, - this->kvec_d[i].y, - this->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; - skpt = ss.str(); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "nkstot_ibz", 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; - - // resize the kpoint container according to nkstot_ibz - if (use_symm || is_mp) - { - this->update_use_ibz(nkstot_ibz, kvec_d_ibz, wk_ibz); - } - - return; -} - -// complement coordinates of k-points according to existing coordinates -// if cartesian coordinates are given, then direct coordinates are calculated -// if direct coordinates are given, then cartesian coordinates are calculated -void K_Vectors::set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt) -{ - - if (PARAM.inp.final_scf) // LiuXh add 20180606 - { - 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; - } - else - { - GlobalV::ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; - } - } - } - - // set cartesian k vectors. - if (!kc_done && kd_done) - { - 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; - } - } - kc_done = true; - } - - // set direct k vectors - else if (kc_done && !kd_done) - { - ModuleBase::Matrix3 RT = R.Transpose(); - 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 - 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; i < nkstot; i++) - { - table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n", - 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) - { - std::stringstream ss; - ss << " " << std::setw(40) << "nkstot now" - << " = " << nkstot << std::endl; - ss << table << std::endl; - skpt = ss.str(); - } - return; -} - void K_Vectors::normalize_wk(const int& degspin) { if (GlobalV::MY_RANK != 0) { @@ -1130,129 +574,6 @@ void K_Vectors::normalize_wk(const int& degspin) return; } -#ifdef __MPI -void K_Vectors::mpi_k() -{ - ModuleBase::TITLE("K_Vectors", "mpi_k"); - - Parallel_Common::bcast_bool(kc_done); - - Parallel_Common::bcast_bool(kd_done); - - Parallel_Common::bcast_int(nspin); - - Parallel_Common::bcast_int(nkstot); - - Parallel_Common::bcast_int(nkstot_full); - - Parallel_Common::bcast_int(nmp, 3); - - kl_segids.resize(nkstot); - Parallel_Common::bcast_int(kl_segids.data(), nkstot); - - Parallel_Common::bcast_double(koffset, 3); - - this->nks = this->para_k.nks_pool[GlobalV::MY_POOL]; - - 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(GlobalV::NPROC, nks_minimum); - - if (nks_minimum == 0) - { - 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); - } - - std::vector isk_aux(nkstot); - std::vector wk_aux(nkstot); - std::vector kvec_c_aux(nkstot * 3); - std::vector kvec_d_aux(nkstot * 3); - - // collect and process in rank 0 - if (GlobalV::MY_RANK == 0) - { - 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; - } - } - - // broadcast k point data to all processors - 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); - - // process k point data in each processor - this->renew(this->nks * this->nspin); - - // distribute - int k_index = 0; - - for (int i = 0; i < 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]; - 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]; - } - -#ifdef __EXX - if (ModuleSymmetry::Symmetry::symm_flag == 1) - {//bcast kstars - this->kstars.resize(nkstot); - for (int ikibz = 0;ikibz < nkstot;++ikibz) - { - int starsize = this->kstars[ikibz].size(); - Parallel_Common::bcast_int(starsize); - GlobalV::ofs_running << "starsize: " << starsize << std::endl; - 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); - GlobalV::ofs_running << "isym: " << isym << " ks_vec: " << ks_vec.x << " " << ks_vec.y << " " << ks_vec.z << std::endl; - if (GlobalV::MY_RANK != 0) - { - kstars[ikibz].insert(std::make_pair(isym, ks_vec)); - } - } - } - } -#endif -} // END SUBROUTINE -#endif - //---------------------------------------------------------- // This routine sets the k vectors for the up and down spin //---------------------------------------------------------- @@ -1305,121 +626,4 @@ void K_Vectors::set_kup_and_kdw() } return; -} // end subroutine set_kup_and_kdw - -void K_Vectors::print_klists(std::ofstream& ofs) -{ - ModuleBase::TITLE("K_Vectors", "print_klists"); - - if (nkstot < nks) - { - std::cout << "\n nkstot=" << nkstot; - std::cout << "\n nks=" << 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; i < nks; i++) - { - table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n", - 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; i < nks; i++) - { - table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n", - 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) -{ - 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); - - // set cartesian k vectors. - kd_done = true; - kc_done = false; - if (!kc_done && kd_done) - { - 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; - } - } - kc_done = true; - } - - // set direct k vectors - else if (kc_done && !kd_done) - { - ModuleBase::Matrix3 RT = latvec.Transpose(); - for (int i = 0; i < nks; i++) - { - 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; i < nks; i++) - { - table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n", - 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->print_klists(GlobalV::ofs_running); - - return; -} +} // end subroutine set_kup_and_kdw \ No newline at end of file diff --git a/source/module_cell/klist.h b/source/module_cell/klist.h index 7e44381c06..a832c1959f 100644 --- a/source/module_cell/klist.h +++ b/source/module_cell/klist.h @@ -6,6 +6,7 @@ #include "module_base/matrix3.h" #include "module_cell/unitcell.h" #include "parallel_kpoints.h" +#include "k_vector_utils.h" #include class K_Vectors @@ -26,6 +27,9 @@ class K_Vectors /// dim: [iks_ibz][(isym, kvec_d)] std::vector>> kstars; + bool kc_done = false; + bool kd_done = false; + K_Vectors(){}; ~K_Vectors(){}; K_Vectors& operator=(const K_Vectors&) = default; @@ -61,51 +65,6 @@ class K_Vectors const ModuleBase::Matrix3& latvec, std::ofstream& ofs); - /** - * @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); - int get_nks() const { return this->nks; @@ -126,6 +85,21 @@ class K_Vectors return this->koffset[i]; } + int get_k_nkstot() const + { + return this->k_nkstot; + } + + int get_nspin() const + { + return this->nspin; + } + + std::string get_k_kword() const + { + return this->k_kword; + } + void set_nks(int value) { this->nks = value; @@ -141,6 +115,11 @@ class K_Vectors this->nkstot_full = value; } + void set_nspin(int value) + { + this->nspin = value; + } + bool get_is_mp() const { return is_mp; @@ -148,17 +127,36 @@ class K_Vectors std::vector ik2iktot; ///<[nks] map ik to the global index of k points + /** + * @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(const int& nkstot_ibz, + const std::vector>& kvec_d_ibz, + const std::vector& wk_ibz); + private: int nks = 0; ///< number of symmetry-reduced k points in this pool(processor, up+dw) int nkstot = 0; ///< number of symmetry-reduced k points in full k mesh int nkstot_full = 0; ///< number of k points before symmetry reduction in full k mesh int nspin = 0; - bool kc_done = false; - bool kd_done = false; double koffset[3] = {0.0}; // used only in automatic k-points. std::string k_kword; // LiuXh add 20180619 - int k_nkstot = 0; // LiuXh add 20180619 + int k_nkstot = 0; // LiuXh add 20180619 // WHAT IS THIS????? bool is_mp = false; // Monkhorst-Pack /** @@ -260,47 +258,7 @@ class K_Vectors // step 2 : set both kvec and kved; normalize weight - /** - * @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(const int& nkstot_ibz, - const std::vector>& kvec_d_ibz, - const std::vector& wk_ibz); - - /** - * @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); + // void set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt); /** * @brief Normalizes the weights of the k-points. @@ -320,26 +278,7 @@ class K_Vectors */ 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 kpoints. @@ -364,31 +303,13 @@ class K_Vectors */ void set_kup_and_kdw(); - // step 5 - // print k lists. - - /** - * @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); - /** * @brief Gets the global index of a k-point. * @return this->ik2iktot[ik] */ void cal_ik_global(); +#ifdef __MPI + friend void KVectorUtils::kvec_mpi_k(K_Vectors& kvec); +#endif }; #endif // KVECT_H \ No newline at end of file diff --git a/source/module_cell/test/CMakeLists.txt b/source/module_cell/test/CMakeLists.txt index d19ceee563..149b394bcc 100644 --- a/source/module_cell/test/CMakeLists.txt +++ b/source/module_cell/test/CMakeLists.txt @@ -68,13 +68,13 @@ AddTest( AddTest( TARGET MODULE_CELL_klist_test LIBS parameter ${math_libs} base device symmetry - SOURCES klist_test.cpp ../klist.cpp ../parallel_kpoints.cpp ../../module_io/output.cpp + SOURCES klist_test.cpp ../klist.cpp ../parallel_kpoints.cpp ../../module_io/output.cpp ../k_vector_utils.cpp ) AddTest( TARGET MODULE_CELL_klist_test_para1 LIBS parameter ${math_libs} base device symmetry - SOURCES klist_test_para.cpp ../klist.cpp ../parallel_kpoints.cpp ../../module_io/output.cpp + SOURCES klist_test_para.cpp ../klist.cpp ../parallel_kpoints.cpp ../../module_io/output.cpp ../k_vector_utils.cpp ) add_test(NAME MODULE_CELL_klist_test_para4 diff --git a/source/module_cell/test/klist_test.cpp b/source/module_cell/test/klist_test.cpp index 4c24097d48..3bedf060fc 100644 --- a/source/module_cell/test/klist_test.cpp +++ b/source/module_cell/test/klist_test.cpp @@ -593,7 +593,9 @@ TEST_F(KlistTest, SetAfterVC) kv->kvec_c[0].x = 0; kv->kvec_c[0].y = 0; kv->kvec_c[0].z = 0; - kv->set_after_vc(PARAM.input.nspin, ucell.G, ucell.latvec); +// kv->set_after_vc(PARAM.input.nspin, ucell.G, ucell.latvec); + KVectorUtils::set_after_vc(*kv, PARAM.input.nspin, ucell.G); + EXPECT_TRUE(kv->kd_done); EXPECT_TRUE(kv->kc_done); EXPECT_DOUBLE_EQ(kv->kvec_d[0].x, 0); @@ -613,9 +615,10 @@ TEST_F(KlistTest, PrintKlists) kv->kvec_c[0].x = 0; kv->kvec_c[0].y = 0; kv->kvec_c[0].z = 0; - kv->set_after_vc(PARAM.input.nspin, ucell.G, ucell.latvec); +// kv->set_after_vc(PARAM.input.nspin, ucell.G, ucell.latvec); + KVectorUtils::set_after_vc(*kv, PARAM.input.nspin, ucell.G); EXPECT_TRUE(kv->kd_done); - kv->print_klists(GlobalV::ofs_running); + KVectorUtils::print_klists(*kv, GlobalV::ofs_running); GlobalV::ofs_running.close(); remove("tmp_klist_2"); } @@ -630,7 +633,7 @@ 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(KVectorUtils::print_klists(*kv, GlobalV::ofs_running), ::testing::ExitedWithCode(1), ""); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("nkstot < nks")); } @@ -648,29 +651,33 @@ TEST_F(KlistTest, SetBothKvecFinalSCF) kv->kvec_c[0].y = 0.0; kv->kvec_c[0].z = 0.0; std::string skpt; - PARAM.input.final_scf = true; +// PARAM.input.final_scf = true; kv->kd_done = false; kv->kc_done = false; // case 1 kv->k_nkstot = 0; - kv->set_both_kvec(ucell.G, ucell.latvec, skpt); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt); + KVectorUtils::set_both_kvec(*kv, ucell.G, ucell.latvec, skpt); 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); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt); + KVectorUtils::set_both_kvec(*kv, ucell.G, ucell.latvec, skpt); 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); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt); + KVectorUtils::set_both_kvec(*kv, ucell.G, ucell.latvec, skpt); EXPECT_TRUE(kv->kc_done); EXPECT_TRUE(kv->kd_done); // case 4 GlobalV::ofs_warning.open("klist_tmp_warning_8"); kv->k_kword = "arbitrary"; - kv->set_both_kvec(ucell.G, ucell.latvec, skpt); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt); + KVectorUtils::set_both_kvec(*kv, ucell.G, ucell.latvec, skpt); GlobalV::ofs_warning.close(); ifs.open("klist_tmp_warning_8"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); @@ -691,12 +698,14 @@ TEST_F(KlistTest, SetBothKvec) kv->kc_done = false; kv->kd_done = true; std::string skpt; - PARAM.input.final_scf = false; - kv->set_both_kvec(ucell.G, ucell.latvec, skpt); +// PARAM.input.final_scf = false; +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt); + KVectorUtils::set_both_kvec(*kv, ucell.G, ucell.latvec, skpt); EXPECT_TRUE(kv->kc_done); kv->kc_done = true; kv->kd_done = false; - kv->set_both_kvec(ucell.G, ucell.latvec, skpt); +// kv->set_both_kvec(ucell.G, ucell.latvec, skpt); + KVectorUtils::set_both_kvec(*kv, ucell.G, ucell.latvec, skpt); EXPECT_TRUE(kv->kd_done); } @@ -743,7 +752,7 @@ TEST_F(KlistTest, IbzKpoint) std::string skpt; ModuleSymmetry::Symmetry::symm_flag = 1; bool match = true; - kv->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt, ucell, match); + KVectorUtils::kvec_ibz_kpoint(*kv, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, ucell, match); EXPECT_EQ(kv->get_nkstot(), 35); GlobalV::ofs_running << skpt << std::endl; GlobalV::ofs_running.close(); @@ -768,7 +777,7 @@ TEST_F(KlistTest, IbzKpointIsMP) std::string skpt; ModuleSymmetry::Symmetry::symm_flag = 0; bool match = true; - kv->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag, skpt, ucell, match); + KVectorUtils::kvec_ibz_kpoint(*kv, symm, ModuleSymmetry::Symmetry::symm_flag, skpt, ucell, match); EXPECT_EQ(kv->get_nks(), 260); GlobalV::ofs_running << skpt << std::endl; GlobalV::ofs_running.close(); diff --git a/source/module_cell/test/klist_test_para.cpp b/source/module_cell/test/klist_test_para.cpp index 17f9147864..f7cb12fec1 100644 --- a/source/module_cell/test/klist_test_para.cpp +++ b/source/module_cell/test/klist_test_para.cpp @@ -309,7 +309,8 @@ TEST_F(KlistParaTest, SetAfterVC) } // call set_after_vc here kv->kc_done = false; - kv->set_after_vc(kv->nspin, ucell.G, ucell.latvec); +// kv->set_after_vc(kv->nspin, ucell.G, ucell.latvec); + KVectorUtils::set_after_vc(*kv, kv->nspin, ucell.G); EXPECT_TRUE(kv->kc_done); EXPECT_TRUE(kv->kd_done); // clear diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index f32a389158..9e7e62b60b 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -17,6 +17,7 @@ #include "module_io/write_elecstat_pot.h" #include "module_io/write_elf.h" #include "module_parameter/parameter.h" +#include "module_cell/k_vector_utils.h" namespace ModuleESolver { @@ -292,7 +293,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) } // reset k-points - kv.set_after_vc(PARAM.inp.nspin, ucell.G, ucell.latvec); + KVectorUtils::set_after_vc(kv, PARAM.inp.nspin, ucell.G); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); } diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index 34eb7ca070..67032896b3 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -48,8 +48,8 @@ AddTest( AddTest( TARGET MODULE_IO_write_istate_info_test LIBS parameter ${math_libs} base device symmetry - SOURCES write_istate_info_test.cpp ../write_istate_info.cpp ../output.cpp ../../module_cell/parallel_kpoints.cpp ../../module_cell/klist.cpp - ../cif_io.cpp + SOURCES write_istate_info_test.cpp ../write_istate_info.cpp ../output.cpp ../../module_cell/parallel_kpoints.cpp ../../module_cell/klist.cpp ../../module_cell/k_vector_utils.cpp + ../cif_io.cpp ) AddTest( @@ -60,14 +60,14 @@ AddTest( AddTest( TARGET MODULE_IO_write_dos_pw - LIBS parameter ${math_libs} base device symmetry - SOURCES write_dos_pw_test.cpp ../cal_dos.cpp ../write_dos_pw.cpp ../output.cpp ../../module_cell/parallel_kpoints.cpp ../../module_cell/klist.cpp ../nscf_fermi_surf.cpp + LIBS parameter ${math_libs} base device symmetry + SOURCES write_dos_pw_test.cpp ../cal_dos.cpp ../write_dos_pw.cpp ../output.cpp ../../module_cell/parallel_kpoints.cpp ../../module_cell/klist.cpp ../nscf_fermi_surf.cpp ../../module_cell/k_vector_utils.cpp ) AddTest( TARGET MODULE_IO_print_info LIBS parameter ${math_libs} base device symmetry cell_info - SOURCES print_info_test.cpp ../print_info.cpp ../output.cpp ../../module_cell/klist.cpp ../../module_cell/parallel_kpoints.cpp + SOURCES print_info_test.cpp ../print_info.cpp ../output.cpp ../../module_cell/klist.cpp ../../module_cell/parallel_kpoints.cpp ../../module_cell/k_vector_utils.cpp ) AddTest( diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index b3a6911394..fd3e29aade 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -115,7 +115,7 @@ struct Input_para double scf_thr = -1.0; ///< \sum |rhog_out - rhog_in |^2 double scf_ene_thr = -1.0; ///< energy threshold for scf convergence, in eV int scf_thr_type = -1; ///< type of the criterion of scf_thr, 1: reci drho, 2: real drho - bool final_scf = false; ///< whether to do final scf +// bool final_scf = false; ///< whether to do final scf bool scf_os_stop = false; ///< whether to stop scf when oscillation is detected double scf_os_thr = -0.01; ///< drho threshold for oscillation int scf_os_ndim = 0; ///< number of old iterations used for oscillation detection