From 31be520eb79fef5d9039a228ca0516cb2f4cbd04 Mon Sep 17 00:00:00 2001 From: Yike Huang Date: Wed, 12 Jul 2023 18:35:31 +0800 Subject: [PATCH 1/4] Fix issue #2721: a correction on atom sequence check Fix issue #2721 --- source/module_cell/read_atoms.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp index 8e8019cb248..64c9a1e29cc 100644 --- a/source/module_cell/read_atoms.cpp +++ b/source/module_cell/read_atoms.cpp @@ -415,16 +415,17 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn // start magnetization //======================================= ModuleBase::GlobalFunc::READ_VALUE(ifpos, atoms[it].label); - bool found = false; + bool sequence_match = false; for(int it2=0; it2atoms[it].label == this->atom_label[it] ) - { - found = true; + if( this->atoms[it].label == this->atom_label[it2] ) + { + if (it2 == it) sequence_match = true; } } - if(!found) + if(!sequence_match) { + ofs_warning << " Label sequences in ATOMIC_POSITIONS and ATOMIC_SPECIES sections do not match!" << std::endl; ofs_warning << " Label read from ATOMIC_POSITIONS is " << this->atoms[it].label << std::endl; ofs_warning << " Label from ATOMIC_SPECIES is " << this->atom_label[it] << std::endl; return 0; From 58edeba1e02c408025d724c528e836ffe76d89da Mon Sep 17 00:00:00 2001 From: Yike Huang Date: Thu, 20 Jul 2023 17:58:18 +0800 Subject: [PATCH 2/4] Update read_atoms.cpp --- source/module_cell/read_atoms.cpp | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp index 64c9a1e29cc..0cb40e0c926 100644 --- a/source/module_cell/read_atoms.cpp +++ b/source/module_cell/read_atoms.cpp @@ -415,17 +415,9 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn // start magnetization //======================================= ModuleBase::GlobalFunc::READ_VALUE(ifpos, atoms[it].label); - bool sequence_match = false; - for(int it2=0; it2atoms[it].label != this->atom_label[it]) { - if( this->atoms[it].label == this->atom_label[it2] ) - { - if (it2 == it) sequence_match = true; - } - } - if(!sequence_match) - { - ofs_warning << " Label sequences in ATOMIC_POSITIONS and ATOMIC_SPECIES sections do not match!" << std::endl; + ofs_warning << " Label orders in ATOMIC_POSITIONS and ATOMIC_SPECIES sections do not match!" << std::endl; ofs_warning << " Label read from ATOMIC_POSITIONS is " << this->atoms[it].label << std::endl; ofs_warning << " Label from ATOMIC_SPECIES is " << this->atom_label[it] << std::endl; return 0; From 06999820e9c60901aa3bb1f393c0935838c4bf0c Mon Sep 17 00:00:00 2001 From: Yike Huang Date: Wed, 26 Jul 2023 09:19:42 +0800 Subject: [PATCH 3/4] Update klist.cpp --- source/module_cell/klist.cpp | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp index c7a4186397d..c02bd501e81 100644 --- a/source/module_cell/klist.cpp +++ b/source/module_cell/klist.cpp @@ -195,6 +195,7 @@ bool K_Vectors::read_kpoints(const std::string &fn) 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); + /* b(i)/lat0*2pi will get un-rescaled Cartesian b(i) */ 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)); @@ -203,9 +204,10 @@ bool K_Vectors::read_kpoints(const std::string &fn) std::ofstream ofs(fn.c_str()); ofs << "K_POINTS" << std::endl; ofs << "0" << std::endl; - ofs << "Gamma" << std::endl; + ofs << "Gamma" << std::endl; /* it means the center is gamma point. */ ofs << nk1 << " " << nk2 << " " << nk3 <<" 0 0 0" << std::endl; ofs.close(); + /* so KPT file is directly omitted and rewritten if kspacing this keyword is set */ } std::ifstream ifk(fn.c_str()); @@ -689,11 +691,12 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s // convert kgmatrix to k-lattice ModuleBase::Matrix3* kkmatrix = new ModuleBase::Matrix3 [nrotkm]; - symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, gk); + symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, gk); /* does not convert anything? */ + /* for not symm case, kkmatrix == kgmatrix */ // direct coordinates of k-points in k-lattice std::vector> kvec_d_k(nkstot); - for (int i=0;inkstot_ibz = 0; @@ -709,11 +712,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s 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; @@ -730,17 +729,18 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s }; // for output in kpoints file int ibz_index[nkstot]; - // search in all k-poins. + // search in all presently degenerated k-poins. for (int i = 0; i < nkstot; ++i) { - //restrict to [0, 1) + //restrict to [-0.5, 0.5) 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; - + + /// @note try every single symmetrical operation for (int j = 0; j < nrotkm; ++j) { if (!already_exist) @@ -753,7 +753,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s // 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); + restrict_kpt(kvec_rot); /* fold to inteval [-0.5, +0.5) */ 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 @@ -768,6 +768,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s // std::cout << "\n kvec_rot = " << kvec_rot.x << " " << kvec_rot.y << " " << kvec_rot.z; + /// @note scan already saved irreducible k-points for (int k=0; k< this->nkstot_ibz; ++k) { if ( symm.equal(kvec_rot.x, this->kvec_d_ibz[k].x) && From e1b3a39d3e00beeb84188a802c3915884b0614db Mon Sep 17 00:00:00 2001 From: Yike Huang Date: Wed, 26 Jul 2023 10:29:33 +0800 Subject: [PATCH 4/4] Revert "Update klist.cpp" This reverts commit 06999820e9c60901aa3bb1f393c0935838c4bf0c. --- source/module_cell/klist.cpp | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/source/module_cell/klist.cpp b/source/module_cell/klist.cpp index c02bd501e81..c7a4186397d 100644 --- a/source/module_cell/klist.cpp +++ b/source/module_cell/klist.cpp @@ -195,7 +195,6 @@ bool K_Vectors::read_kpoints(const std::string &fn) 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); - /* b(i)/lat0*2pi will get un-rescaled Cartesian b(i) */ 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)); @@ -204,10 +203,9 @@ bool K_Vectors::read_kpoints(const std::string &fn) std::ofstream ofs(fn.c_str()); ofs << "K_POINTS" << std::endl; ofs << "0" << std::endl; - ofs << "Gamma" << std::endl; /* it means the center is gamma point. */ + ofs << "Gamma" << std::endl; ofs << nk1 << " " << nk2 << " " << nk3 <<" 0 0 0" << std::endl; ofs.close(); - /* so KPT file is directly omitted and rewritten if kspacing this keyword is set */ } std::ifstream ifk(fn.c_str()); @@ -691,12 +689,11 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s // convert kgmatrix to k-lattice ModuleBase::Matrix3* kkmatrix = new ModuleBase::Matrix3 [nrotkm]; - symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, gk); /* does not convert anything? */ - /* for not symm case, kkmatrix == kgmatrix */ + symm.gmatrix_convert(kgmatrix.data(), kkmatrix, nrotkm, ucell.G, gk); // direct coordinates of k-points in k-lattice std::vector> kvec_d_k(nkstot); - for (int i=0;inkstot_ibz = 0; @@ -712,7 +709,11 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s ModuleBase::Vector3 kvec_rot; ModuleBase::Vector3 kvec_rot_k; - /// @brief PBC, wrap kvec_d into [-0.5, 0.5) + +// 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; @@ -729,18 +730,17 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s }; // for output in kpoints file int ibz_index[nkstot]; - // search in all presently degenerated k-poins. + // search in all k-poins. for (int i = 0; i < nkstot; ++i) { - //restrict to [-0.5, 0.5) + //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; - - /// @note try every single symmetrical operation + for (int j = 0; j < nrotkm; ++j) { if (!already_exist) @@ -753,7 +753,7 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s // 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); /* fold to inteval [-0.5, +0.5) */ + restrict_kpt(kvec_rot); 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 @@ -768,7 +768,6 @@ void K_Vectors::ibz_kpoint(const ModuleSymmetry::Symmetry &symm, bool use_symm,s // std::cout << "\n kvec_rot = " << kvec_rot.x << " " << kvec_rot.y << " " << kvec_rot.z; - /// @note scan already saved irreducible k-points for (int k=0; k< this->nkstot_ibz; ++k) { if ( symm.equal(kvec_rot.x, this->kvec_d_ibz[k].x) &&