From 57985c4e77069d3f9d156f3a7a5109b9e8404838 Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Mon, 24 Feb 2025 17:44:52 +0800 Subject: [PATCH 1/2] analyze magnetic group without time-reversal symmetry --- .../module_cell/module_symmetry/symmetry.cpp | 118 +++++++++++------- source/module_cell/module_symmetry/symmetry.h | 9 +- .../irreducible_sector_bvk.cpp | 2 +- 3 files changed, 78 insertions(+), 51 deletions(-) diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index bc3f4c21e6..018a316492 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -1,5 +1,6 @@ #include #include +#include #include "module_parameter/parameter.h" #include "symmetry.h" #include "module_parameter/parameter.h" @@ -7,7 +8,7 @@ #include "module_base/mathzone.h" #include "module_base/constants.h" #include "module_base/timer.h" - +#include "module_io/output.h" namespace ModuleSymmetry { int Symmetry::symm_flag = 0; @@ -124,49 +125,8 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, } if (!pricell_loop && PARAM.inp.nspin == 2) - {//analyze symmetry for spin-up atoms only - std::vector pos_spinup; - for (int it = 0;it < ntype;++it) - { - int na_spinup = 0; - for (int ia = 0; ia < atoms[it].na; ++ia) { - if (atoms[it].mag[ia] > -this->epsilon) { - ++na_spinup; - } - } - this->na[it] = na_spinup; - //update newpos - for (int ia = 0;ia < atoms[it].na;++ia) - { - if (atoms[it].mag[ia] > -this->epsilon) - { - pos_spinup.push_back(this->newpos[3 * (istart[it] + ia)]); - pos_spinup.push_back(this->newpos[3 * (istart[it] + ia) + 1]); - pos_spinup.push_back(this->newpos[3 * (istart[it] + ia) + 2]); - } - } - // update start to spin-up configuration - if (it > 0) { - istart[it] = istart[it - 1] + na[it - 1]; - } - if (na[it] < na[itmin_type]) - { - this->itmin_type = it; - this->itmin_start = istart[it]; - } - } - this->getgroup(nrot_out, nrotk_out, ofs_running, this->nop, this->symop, this->gmatrix, this->gtrans, - pos_spinup.data(), this->rotpos, this->index, this->itmin_type, this->itmin_start, this->istart, this->na); - // recover na and istart - for (int it = 0;it < ntype;++it) - { - this->na[it] = atoms[it].na; - if (it > 0) { - istart[it] = istart[it - 1] + na[it - 1]; - } - } - // For AFM analysis End - //------------------------------------------------------------ + { + this->analyze_magnetic_group(atoms, st, nrot_out, nrotk_out); } else { @@ -174,7 +134,7 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, // nrot_out: the number of pure point group rotations // nrotk_out: the number of all space group operations this->getgroup(nrot_out, nrotk_out, ofs_running, this->nop, this->symop, this->gmatrix, this->gtrans, - this->newpos, this->rotpos, this->index, this->itmin_type, this->itmin_start, this->istart, this->na); + this->newpos, this->rotpos, this->index, this->ntype, this->itmin_type, this->itmin_start, this->istart, this->na); } }; @@ -889,7 +849,7 @@ void Symmetry::lattice_type( void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, const int& nop, const ModuleBase::Matrix3* symop, ModuleBase::Matrix3* gmatrix, ModuleBase::Vector3* gtrans, - double* pos, double* rotpos, int* index, const int itmin_type, const int itmin_start, int* istart, int* na)const + double* pos, double* rotpos, int* index, const int ntype, const int itmin_type, const int itmin_start, int* istart, int* na)const { ModuleBase::TITLE("Symmetry", "getgroup"); @@ -913,7 +873,7 @@ void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, const //std::cout << "nop = " <checksym(symop[i], gtrans[i], pos, rotpos, index, itmin_type, itmin_start, istart, na); + bool s_flag = this->checksym(symop[i], gtrans[i], pos, rotpos, index, ntype, itmin_type, itmin_start, istart, na); if (s_flag == 1) { //------------------------------ @@ -992,7 +952,7 @@ void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, const } bool Symmetry::checksym(const ModuleBase::Matrix3& s, ModuleBase::Vector3& gtrans, - double* pos, double* rotpos, int* index, const int itmin_type, const int itmin_start, int* istart, int* na)const + double* pos, double* rotpos, int* index, const int ntype, const int itmin_type, const int itmin_start, int* istart, int* na)const { //---------------------------------------------- // checks whether a point group symmetry element @@ -2297,4 +2257,66 @@ bool Symmetry::is_all_movable(const Atom* atoms, const Statistics& st)const } return all_mbl; } + +void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out) +{ + // 1. classify atoms with different magmom + // (use symmetry_prec to judge if two magmoms are the same) + std::vector> mag_type_atoms; + for (int it = 0;it < ntype;++it) + { + for (int ia = 0; ia < atoms[it].na; ++ia) + { + bool find = false; + for (auto& mt : mag_type_atoms) + { + const int mag_iat = *mt.begin(); + const int mag_it = st.iat2it[mag_iat]; + const int mag_ia = st.iat2ia[mag_iat]; + if (it == mag_it && this->equal(atoms[it].mag[ia], atoms[mag_it].mag[mag_ia])) + { + mt.insert(st.itia2iat(it, ia)); + find = true; + break; + } + } + if (!find) + { + mag_type_atoms.push_back(std::set({ st.itia2iat(it,ia) })); + } + } + } + + // 2. get the start index, number of atoms and positions for each mag_type + std::vector mag_istart(mag_type_atoms.size()); + std::vector mag_na(mag_type_atoms.size()); + std::vector mag_pos; + int mag_itmin_type = 0; + int mag_itmin_start = 0; + for (int mag_it = 0;mag_it < mag_type_atoms.size(); ++mag_it) + { + mag_na[mag_it] = mag_type_atoms.at(mag_it).size(); + if (mag_it > 0) + { + mag_istart[mag_it] = mag_istart[mag_it - 1] + mag_na[mag_it - 1]; + } + if (mag_na[mag_it] < mag_na[itmin_type]) + { + mag_itmin_type = mag_it; + mag_itmin_start = mag_istart[mag_it]; + } + for (auto& mag_iat : mag_type_atoms.at(mag_it)) + { + + const int index_in_newpos = 3 * (this->istart[st.iat2it[mag_iat]] + st.iat2ia[mag_iat]); + for (int i = 0; i < 3; ++i) + { + mag_pos.push_back(this->newpos[index_in_newpos + i]); + } + } + } + // 3. analyze the effective structure + this->getgroup(nrot_out, nrotk_out, GlobalV::ofs_running, this->nop, this->symop, this->gmatrix, this->gtrans, + mag_pos.data(), this->rotpos, this->index, mag_type_atoms.size(), mag_itmin_type, mag_itmin_start, mag_istart.data(), mag_na.data()); +} } diff --git a/source/module_cell/module_symmetry/symmetry.h b/source/module_cell/module_symmetry/symmetry.h index e5f4ad7950..c3ef81607d 100644 --- a/source/module_cell/module_symmetry/symmetry.h +++ b/source/module_cell/module_symmetry/symmetry.h @@ -91,9 +91,10 @@ class Symmetry : public Symmetry_Basic void getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, const int& nop, const ModuleBase::Matrix3* symop, ModuleBase::Matrix3* gmatrix, ModuleBase::Vector3* gtrans, - double* pos, double* rotpos, int* index, const int itmin_type, const int itmin_start, int* istart, int* na)const; + double* pos, double* rotpos, int* index, const int ntype, const int itmin_type, const int itmin_start, int* istart, int* na)const; bool checksym(const ModuleBase::Matrix3& s, ModuleBase::Vector3& gtrans, - double* pos, double* rotpos, int* index, const int itmin_type, const int itmin_start, int* istart, int* na)const; + double* pos, double* rotpos, int* index, const int itmin_type, + const int ntype, const int itmin_start, int* istart, int* na)const; /// @brief primitive cell analysis void pricell(double* pos, const Atom* atoms); @@ -145,6 +146,10 @@ class Symmetry : public Symmetry_Basic /// Loop the magmom of each atoms in its type when NSPIN>1. If not all the same, primitive cells should not be looped in rhog_symmetry. bool magmom_same_check(const Atom* atoms)const; + /// Analyze magnetic group without time-reversal symmetry + /// (because currently the charge density symmetrization does not support it) + /// Method: treat atoms with different magmom as atoms of different type + void analyze_magnetic_group(const Atom* atoms, const Statistics& st, int& nrot_out, int& nrotk_out); }; } diff --git a/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp b/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp index 5efd046d73..50df31b360 100644 --- a/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp +++ b/source/module_ri/module_exx_symmetry/irreducible_sector_bvk.cpp @@ -136,7 +136,7 @@ namespace ModuleSymmetry symm.getgroup(bvk_npg, bvk_nsg, GlobalV::ofs_running, bvk_nop, bvk_op.data(), bvk_gmatrix.data(), bvk_gtrans.data(), bvk_dpos.data(), bvk_rot_dpos.data(), order_index.data(), - bvk_itmin_type, bvk_itmin_start, bvk_istart.data(), bvk_na.data()); + symm.ntype, bvk_itmin_type, bvk_itmin_start, bvk_istart.data(), bvk_na.data()); bvk_gmatrix.resize(bvk_nsg); bvk_gtrans.resize(bvk_nsg); this->bvk_nsym_ = bvk_nsg; From d2f33eb4ae4b96a8a98d17515adf293f7049297e Mon Sep 17 00:00:00 2001 From: maki49 <1579492865@qq.com> Date: Mon, 24 Feb 2025 19:54:50 +0800 Subject: [PATCH 2/2] fix: need to calculate direct coordinates again --- source/module_cell/module_symmetry/symmetry.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index 018a316492..7d892725eb 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -2307,11 +2307,15 @@ void Symmetry::analyze_magnetic_group(const Atom* atoms, const Statistics& st, i } for (auto& mag_iat : mag_type_atoms.at(mag_it)) { - - const int index_in_newpos = 3 * (this->istart[st.iat2it[mag_iat]] + st.iat2ia[mag_iat]); + // this->newpos have been ordered by original structure(ntype, na), it cannot be directly used here. + // we need to reset the calculate again the coordinate of the new structure. + const ModuleBase::Vector3 direct_tmp = atoms[st.iat2it[mag_iat]].tau[st.iat2ia[mag_iat]] * this->optlat.Inverse(); + std::array direct = { direct_tmp.x, direct_tmp.y, direct_tmp.z }; for (int i = 0; i < 3; ++i) { - mag_pos.push_back(this->newpos[index_in_newpos + i]); + this->check_translation(direct[i], -floor(direct[i])); + this->check_boundary(direct[i]); + mag_pos.push_back(direct[i]); } } }