From 79f092ce2379664af410557793f986624e1a2e0a Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Thu, 10 Jul 2025 23:25:53 +0800 Subject: [PATCH 1/3] Feature: support output the matrix representation of symmetry operation --- source/source_cell/k_vector_utils.cpp | 12 +- .../module_symmetry/symmetry_basic.cpp | 204 ++++++++++-------- .../module_parameter/input_parameter.h | 7 +- source/source_io/read_input_item_output.cpp | 26 +++ 4 files changed, 142 insertions(+), 107 deletions(-) diff --git a/source/source_cell/k_vector_utils.cpp b/source/source_cell/k_vector_utils.cpp index d21a5ceb81b..0ea7cb77da8 100644 --- a/source/source_cell/k_vector_utils.cpp +++ b/source/source_cell/k_vector_utils.cpp @@ -486,15 +486,9 @@ void kvec_ibz_kpoint(K_Vectors& kv, 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); + 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 diff --git a/source/source_cell/module_symmetry/symmetry_basic.cpp b/source/source_cell/module_symmetry/symmetry_basic.cpp index 6966256100f..c5dace25431 100644 --- a/source/source_cell/module_symmetry/symmetry_basic.cpp +++ b/source/source_cell/module_symmetry/symmetry_basic.cpp @@ -1,6 +1,7 @@ #include "symmetry.h" #include "source_base/mymath.h" #include "source_io/module_parameter/parameter.h" +#include "source_base/formatter.h" bool ModuleSymmetry::test_brav = 0; @@ -11,15 +12,15 @@ std::string Symmetry_Basic::get_brav_name(const int ibrav) const { switch(ibrav) { - case 1: return "01. Cubic P (simple)"; - case 2: return "02. Cubic I (body-centered)"; - case 3: return "03. Cubic F (face-centered)"; - case 4: return "04. Hexagonal cell"; - case 5: return "05. Tetrogonal P (simple)"; - case 6: return "06. Tetrogonal I (body-centered)"; - case 7: return "07. Rhombohedral (Trigonal) cell"; - case 8: return "08. Orthorhombic P(simple)"; - case 9: return "09. Orthorhombic I (body-centered)"; + case 1: return "01. Cubic P (simple)"; + case 2: return "02. Cubic I (body-centered)"; + case 3: return "03. Cubic F (face-centered)"; + case 4: return "04. Hexagonal cell"; + case 5: return "05. Tetrogonal P (simple)"; + case 6: return "06. Tetrogonal I (body-centered)"; + case 7: return "07. Rhombohedral (Trigonal) cell"; + case 8: return "08. Orthorhombic P(simple)"; + case 9: return "09. Orthorhombic I (body-centered)"; case 10: return "10. Orthorhombic F (face-centered)"; case 11: return "11. Orthorhombic C (base-centered)"; case 12: return "12. Monoclinic P (simple)"; @@ -27,7 +28,8 @@ std::string Symmetry_Basic::get_brav_name(const int ibrav) const case 14: return "14. Triclinic cell"; case 15: return "wrong !! "; } - return "Congratulations!You have found a bravais lattice that never existed!"; + // return "Congratulations! You have found a bravais lattice that never existed!"; + return "Unknown Bravais lattice"; } // Control the accuracy @@ -45,7 +47,6 @@ bool Symmetry_Basic::equal(const double &m,const double &n)const void Symmetry_Basic::check_boundary(double &x)const { if(equal(x,-0.5) || equal(x,0.5)) x=-0.5; - return; } double Symmetry_Basic::get_translation_vector(const double& x1, const double& x2) const @@ -53,7 +54,7 @@ double Symmetry_Basic::get_translation_vector(const double& x1, const double& x2 double t=0.0; // "t"ranslation t = x2 - x1; t = fmod(t+100.0, 1.0); - if( fabs(t-1) < epsilon * 0.5) t = 0.0; + if( fabs(t-1) < epsilon * 0.5) { t = 0.0; } return t; } @@ -187,8 +188,8 @@ void Symmetry_Basic::matrigen(ModuleBase::Matrix3 *symgen, const int ngen, Modul int m2=0; int n=0; - ModuleBase::Matrix3 iden(1,0,0,0,1,0,0,0,1); - ModuleBase::Matrix3 sig(1,0,0,0,1,0,0,0,1); + ModuleBase::Matrix3 iden(1,0,0,0,1,0,0,0,1); + ModuleBase::Matrix3 sig(1,0,0,0,1,0,0,0,1); ModuleBase::Matrix3 temp1(1,0,0,0,1,0,0,0,1); ModuleBase::Matrix3 temp2(1,0,0,0,1,0,0,0,1); @@ -321,25 +322,26 @@ void Symmetry_Basic::matrigen(ModuleBase::Matrix3 *symgen, const int ngen, Modul //-------------------------------------------------------------- void Symmetry_Basic::setgroup(ModuleBase::Matrix3* symop, int &nop, const int &ibrav) const { - if(PARAM.inp.test_symmetry) ModuleBase::TITLE("Symmetry_Basic","setgroup"); - + if(PARAM.inp.out_symm_reprmat[0] > 1) { + ModuleBase::TITLE("Symmetry_Basic", "setgroup"); + } ModuleBase::Matrix3 symgen[3]; - ModuleBase::Matrix3 inv(-1,0,0,0,-1,0,0,0,-1); - ModuleBase::Matrix3 r3d(0,1,0,0,0,1,1,0,0); - ModuleBase::Matrix3 r6z(1,1,0,-1,0,0,0,0,1); - ModuleBase::Matrix3 r2hex(1,0,0,-1,-1,0,0,0,-1); - ModuleBase::Matrix3 r2tri(-1,0,0,0,0,-1,0,-1,0); - ModuleBase::Matrix3 r4zp(0,1,0,-1,0,0,0,0,1); - ModuleBase::Matrix3 r2yp(-1,0,0,0,1,0,0,0,-1); - ModuleBase::Matrix3 r4zbc(0,0,-1,1,1,1,0,-1,0); - ModuleBase::Matrix3 r4zfc(1,0,-1,1,0,0,1,-1,0); - ModuleBase::Matrix3 r2zp(-1,0,0,0,-1,0,0,0,1); - ModuleBase::Matrix3 r2ybc(0,0,1,-1,-1,-1,1,0,0); - ModuleBase::Matrix3 r2zbc(0,1,0,1,0,0,-1,-1,-1); - ModuleBase::Matrix3 r2ybas(0,-1,0,-1,0,0,0,0,-1); - ModuleBase::Matrix3 r2yfc(0,-1,1,0,-1,0,1,-1,0); - ModuleBase::Matrix3 r2zfc(0,1,-1,1,0,-1,0,0,-1); + ModuleBase::Matrix3 inv(-1, 0, 0, 0,-1, 0, 0, 0,-1); + ModuleBase::Matrix3 r3d( 0, 1, 0, 0, 0, 1, 1, 0, 0); + ModuleBase::Matrix3 r6z( 1, 1, 0,-1, 0, 0, 0, 0, 1); + ModuleBase::Matrix3 r2hex( 1, 0, 0,-1,-1, 0, 0, 0,-1); + ModuleBase::Matrix3 r2tri(-1, 0, 0, 0, 0,-1, 0,-1, 0); + ModuleBase::Matrix3 r4zp( 0, 1, 0,-1, 0, 0, 0, 0, 1); + ModuleBase::Matrix3 r2yp(-1, 0, 0, 0, 1, 0, 0, 0,-1); + ModuleBase::Matrix3 r4zbc( 0, 0,-1, 1, 1, 1, 0,-1, 0); + ModuleBase::Matrix3 r4zfc( 1, 0,-1, 1, 0, 0, 1,-1, 0); + ModuleBase::Matrix3 r2zp(-1, 0, 0, 0,-1, 0, 0, 0, 1); + ModuleBase::Matrix3 r2ybc( 0, 0, 1,-1,-1,-1, 1, 0, 0); + ModuleBase::Matrix3 r2zbc( 0, 1, 0, 1, 0, 0,-1,-1,-1); + ModuleBase::Matrix3 r2ybas( 0,-1, 0,-1, 0, 0, 0, 0,-1); + ModuleBase::Matrix3 r2yfc( 0,-1, 1, 0,-1, 0, 1,-1, 0); + ModuleBase::Matrix3 r2zfc( 0, 1,-1, 1, 0,-1, 0, 0,-1); //the pure translation lattice (bravais lattice) has some maximum symmetry //set first up the point group operations for this symmetry. @@ -428,26 +430,36 @@ void Symmetry_Basic::setgroup(ModuleBase::Matrix3* symop, int &nop, const int &i if(test_brav) { - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Number of rotation matrices",nop); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Number of rotation matrices", nop); } - if(PARAM.inp.test_symmetry > 1) + // print the symmetry operations + if(PARAM.inp.out_symm_reprmat[0] > 0) { - GlobalV::ofs_running<<" THERE ARE " << nop << " ROTATION MATRICES FOR THE PURE BRAVAIS LATTICE"< op_names = { + + }; + for (int i = 0; i < nop; ++i) { - GlobalV::ofs_running << " " - << std::setw(3) << i + 1 - << std::setw(4) << symop[i].e11 - << std::setw(4) << symop[i].e12 - << std::setw(4) << symop[i].e13 - << std::setw(4) << symop[i].e21 - << std::setw(4) << symop[i].e22 - << std::setw(4) << symop[i].e23 - << std::setw(4) << symop[i].e31 - << std::setw(4) << symop[i].e32 - << std::setw(4) << symop[i].e33 << std::endl; + mat = " " + + FmtCore::format("No. %3d", i + 1) + "\n" + + FmtCore::format(fmtstr.c_str(), symop[i].e11, symop[i].e12, symop[i].e13) + + FmtCore::format(fmtstr.c_str(), symop[i].e21, symop[i].e22, symop[i].e23) + + FmtCore::format(fmtstr.c_str(), symop[i].e31, symop[i].e32, symop[i].e33); + GlobalV::ofs_running << mat << std::endl; } } @@ -462,75 +474,75 @@ int Symmetry_Basic::subgroup(const int& nrot, const int& ninv, { if (ninv) { - // if (nc2 >= 7 && nc3 >= 2 && nc6 >= 2 && ns1 >= 7 && ns3 >= 2 && ns6 >= 2) return 27;//D_6h - if (nc2 >= 3 && nc3 >= 8 && ns1 >= 3 && ns6 >= 8) return 29; //T_h + // if (nc2 >= 7 && nc3 >= 2 && nc6 >= 2 && ns1 >= 7 && ns3 >= 2 && ns6 >= 2) { return 27; } //D_6h + if (nc2 >= 3 && nc3 >= 8 && ns1 >= 3 && ns6 >= 8) { return 29; } //T_h } else { - if (nc2 >= 9 && nc3 >= 8 && nc4 >= 6) return 30; //O - if (nc2 >= 3 && nc3 >= 8 && ns1 >= 6 && ns4 >= 6) return 31; //T_d + if (nc2 >= 9 && nc3 >= 8 && nc4 >= 6) { return 30; } //O + if (nc2 >= 3 && nc3 >= 8 && ns1 >= 6 && ns4 >= 6) { return 31; } //T_d } } if (nrot > 16)//not else if: nrot>24 can also fall in this part and below { - if (ninv && nc2 >= 5 && nc4 >= 2 && ns1 >= 5 && ns4 >= 2) return 20; //D_4h + if (ninv && nc2 >= 5 && nc4 >= 2 && ns1 >= 5 && ns4 >= 2) { return 20; } //D_4h } if (nrot > 12) { if (ninv) { - if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2 && ns1 >= 1 && ns3 >= 2 && ns6 >= 2) return 23;//C_6h - if (nc2 >= 3 && nc3 >= 2 && ns1 >= 3 && ns6 >= 2)return 13;//D_3d + if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2 && ns1 >= 1 && ns3 >= 2 && ns6 >= 2) { return 23; } //C_6h + if (nc2 >= 3 && nc3 >= 2 && ns1 >= 3 && ns6 >= 2) { return 13; } //D_3d } else { - if (nc2 >= 3 && nc3 >= 8)return 28; //T - if (nc2 >= 3 && nc3 >= 2 && ns1 >= 4 && ns3 >= 2) return 26;//D_3h - if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2 && ns1 >= 6) return 25;//C_6v - if (nc2 >= 7 && nc3 >= 2 && nc6 >= 2) return 24;//D_6 + if (nc2 >= 3 && nc3 >= 8) { return 28; } //T + if (nc2 >= 3 && nc3 >= 2 && ns1 >= 4 && ns3 >= 2) { return 26; } //D_3h + if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2 && ns1 >= 6) { return 25; } //C_6v + if (nc2 >= 7 && nc3 >= 2 && nc6 >= 2) { return 24; } //D_6 } } if (nrot > 8) { if (ninv) { - if (nc2 >= 1 && nc4 >= 2 && ns1 >= 1 && ns4 >= 2) return 16;//C_4h - if (nc2 >= 3 && ns1 >= 3)return 8;//D_2h + if (nc2 >= 1 && nc4 >= 2 && ns1 >= 1 && ns4 >= 2) { return 16; } //C_4h + if (nc2 >= 3 && ns1 >= 3) { return 8; } //D_2h } else { - if (nc2 >= 3 && ns1 >= 2 && ns4 >= 2)return 19;//D_2d - if (nc2 >= 1 && nc4 >= 2 && ns1 >= 4) return 18;//C_4v - if (nc2 >= 5 && nc4 >= 2)return 17;//D_4 + if (nc2 >= 3 && ns1 >= 2 && ns4 >= 2) { return 19; } //D_2d + if (nc2 >= 1 && nc4 >= 2 && ns1 >= 4) { return 18; } //C_4v + if (nc2 >= 5 && nc4 >= 2) { return 17; } //D_4 } } if (nrot > 6) { - if (nc3 >= 2 && ns1 >= 1 && ns3 >= 2)return 22;//C_3h - if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2)return 21;//C_6 - if (nc3 >= 2 && ns1 >= 3)return 12;//C_3v - if (nc2 >= 3 && nc3 >= 2)return 11;//D_3 - if (ninv && nc3 >= 2 && ns3 >= 2)return 10;//S_6 + if (nc3 >= 2 && ns1 >= 1 && ns3 >= 2) { return 22; } //C_3h + if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2) { return 21; } //C_6 + if (nc3 >= 2 && ns1 >= 3) { return 12; } //C_3v + if (nc2 >= 3 && nc3 >= 2) { return 11; } //D_3 + if (ninv && nc3 >= 2 && ns3 >= 2) { return 10; }//S_6 } if (nrot > 4) { - if (nc2 >= 1 && ns4 >= 2)return 15;//S_4 - if (nc2 >= 1 && nc4 >= 2)return 14;//C_4 - if (nc2 >= 1 && ns1 >= 2)return 7;//C_2v - if (nc2 >= 3)return 6;//D_2 - if (ninv && nc2 >= 1 && ns1 >= 1)return 5;//C_2h + if (nc2 >= 1 && ns4 >= 2) { return 15; } //S_4 + if (nc2 >= 1 && nc4 >= 2) { return 14; } //C_4 + if (nc2 >= 1 && ns1 >= 2) { return 7; } //C_2v + if (nc2 >= 3) { return 6; } //D_2 + if (ninv && nc2 >= 1 && ns1 >= 1) { return 5; } //C_2h } if (nrot > 3) { - if (nc3 >= 2)return 9;//C_3 + if (nc3 >= 2) { return 9; } //C_3 } if (nrot > 2) { - if (ns1 >= 1)return 4;//C_1h - if (nc2 >= 1)return 3;//C_2 - if (ninv)return 2;//S_2 + if (ns1 >= 1) { return 4; } //C_1h + if (nc2 >= 1) { return 3; } //C_2 + if (ninv) { return 2; } //S_2 } - return 1;//C_1 + return 1; //C_1 } @@ -552,7 +564,9 @@ bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, //there are four trivial cases which could be easily determined //because the number of their elements are exclusive - if (PARAM.inp.test_symmetry) ModuleBase::TITLE("Symmetry_Basic", "pointgroup"); + if (PARAM.inp.out_symm_reprmat[0] > 1) { + ModuleBase::TITLE("Symmetry_Basic", "pointgroup"); + } std::vector pgdict = { "none", "C_1", "S_2", "C_2", "C_1h", "C_2h", "D_2", "C_2v", "D_2h", "C_3", "S_6", "D_3", "C_3v", "D_3d", "C_4", "S_4", "C_4h", @@ -628,26 +642,26 @@ bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, continue; } - if(trace == -1 && det == 1) ++nc2; - else if(trace == 0 && det == 1) ++nc3; - else if(trace == 1 && det == 1) ++nc4; - else if(trace == 2 && det == 1) ++nc6; - else if(trace == 1 && det == -1) ++ns1; - else if(trace == 0 && det == -1) ++ns6; //mohan add 2012-01-15 - else if(trace == -1 && det == -1) ++ns4; - else if(trace == -2 && det == -1) ++ns3; //mohan add 2012-01-15 + if(trace == -1 && det == 1) { ++nc2; } + else if(trace == 0 && det == 1) { ++nc3; } + else if(trace == 1 && det == 1) { ++nc4; } + else if(trace == 2 && det == 1) { ++nc6; } + else if(trace == 1 && det == -1) { ++ns1; } + else if(trace == 0 && det == -1) { ++ns6; } //mohan add 2012-01-15 + else if(trace == -1 && det == -1) { ++ns4; } + else if(trace == -2 && det == -1) { ++ns3; } //mohan add 2012-01-15 } if(test_brav) { - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C2",nc2); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C3",nc3); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C4",nc4); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C6",nc6); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"S1",ns1); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"S3",ns3); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"S4",ns4); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"S6",ns6); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "C2", nc2); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "C3", nc3); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "C4", nc4); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "C6", nc6); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "S1", ns1); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "S3", ns3); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "S4", ns4); + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "S6", ns6); } if(nrot == 2) diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 9b2761f01a0..5e5d746f90f 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -407,8 +407,9 @@ struct Input_para std::vector out_pchg = {}; ///< specify the bands to be calculated for partial charge std::vector out_wfc_norm = {}; ///< specify the bands to be calculated for norm of wfc std::vector out_wfc_re_im = {}; ///< specify the bands to be calculated for real and imaginary parts of wfc - bool if_separate_k = false; ///< whether to write partial charge for all k-points to individual files or merge them - std::vector out_elf = {0, 3}; ///< output the electron localization function (ELF). 0: no; 1: yes + bool if_separate_k = false; ///< whether to write partial charge for all k-points to individual files or merge them + std::vector out_elf = {0, 3}; ///< output the electron localization function (ELF). 0: no; 1: yes + std::vector out_symm_reprmat = {0, 3}; ///< output the symmetry representation matrix // ============== #Parameters (12.Postprocess) =========================== double dos_emin_ev = -15.0; @@ -630,7 +631,7 @@ struct Input_para bool test_stress = false; ///< test the stress. bool test_skip_ewald = false; ///< variables for test only int test_atom_input = false; ///< variables for test_atom_input only - int test_symmetry = false; ///< variables for test_lattice only + // int test_symmetry = false; ///< variables for test_lattice only int test_wf = 0; ///< variables for test_wf only int test_grid = false; ///< variables for test_grid only int test_charge = false; ///< variables for test_vloc only diff --git a/source/source_io/read_input_item_output.cpp b/source/source_io/read_input_item_output.cpp index cc817c53548..b5ccc4ecaaf 100644 --- a/source/source_io/read_input_item_output.cpp +++ b/source/source_io/read_input_item_output.cpp @@ -546,5 +546,31 @@ void ReadInput::item_output() sync_intvec(input.out_elf, 2, 0); this->add_item(item); } + { + // recover the functionality of test_symmetry by introducing a new keyword "out_symm_mat" + // the "out_symm_mat" keyword will be a + Input_Item item("out_symm_reprmat"); + item.annotation = "output matrix representation of symmetry operation into running log file" + " > 0 output the matrix representation of symmetry operation " + ", the second parameter controls the precision, default is 3."; + item.read_value = [](const Input_Item& item, Parameter& para) { + size_t count = item.get_size(); + std::vector out_symm_reprmat(count); // create a placeholder vector + std::transform(item.str_values.begin(), item.str_values.end(), // iterators of 1 + out_symm_reprmat.begin(), // iterator of 2 + [](std::string s){ return std::stoi(s); }); // lambda func + // assign non-negative values to para.input.out_symm_reprmat + std::copy(out_symm_reprmat.begin(), out_symm_reprmat.end(), + para.input.out_symm_reprmat.begin()); + }; + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.out_symm_reprmat[0] < 0 || para.input.out_symm_reprmat[0] > 1) + { + ModuleBase::WARNING_QUIT("ReadInput", "out_symm_reprmat should be 0 or 1"); + } + }; + sync_intvec(input.out_symm_reprmat, 2, 0); + this->add_item(item); + } } } // namespace ModuleIO From 2cb8e595bfe5d31d3fc364b46cf5b24c6c9064cb Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Tue, 15 Jul 2025 19:08:23 +0800 Subject: [PATCH 2/3] Feature: support the output of matrix representation of symm_ops --- .../module_symmetry/symmetry_basic.cpp | 312 +++++++++--------- .../module_parameter/input_parameter.h | 4 +- source/source_io/read_input_item_output.cpp | 18 +- 3 files changed, 167 insertions(+), 167 deletions(-) diff --git a/source/source_cell/module_symmetry/symmetry_basic.cpp b/source/source_cell/module_symmetry/symmetry_basic.cpp index c5dace25431..d09a0cada7a 100644 --- a/source/source_cell/module_symmetry/symmetry_basic.cpp +++ b/source/source_cell/module_symmetry/symmetry_basic.cpp @@ -33,10 +33,10 @@ std::string Symmetry_Basic::get_brav_name(const int ibrav) const } // Control the accuracy -bool Symmetry_Basic::equal(const double &m,const double &n)const +bool Symmetry_Basic::equal(const double &m, const double &n) const { //if( fabs(m-n) < 1.0e-5 ) - if( fabs(m-n) < epsilon ) //LiuXh add 2021-08-12, use accuracy for symmetry + if (fabs(m-n) < epsilon) //LiuXh add 2021-08-12, use accuracy for symmetry { return true; } @@ -110,8 +110,8 @@ void Symmetry_Basic::order_atoms(double* pos, const int& nat, const int* index) return; } -//convert a set of vectors (va) represented in the basis vectors old1, old2, old3 -//to a set of vectors (vb) represented in the basis vectors new1, new2, new3 +// convert a set of vectors (va) represented in the basis vectors old1, old2, old3 +// to a set of vectors (vb) represented in the basis vectors new1, new2, new3 void Symmetry_Basic::veccon( double *carpos, double *rotpos, @@ -181,84 +181,89 @@ void Symmetry_Basic::veccon( } -//generate all point group symmetry operations from the generation group +// generate all point group symmetry operations from the generation group void Symmetry_Basic::matrigen(ModuleBase::Matrix3 *symgen, const int ngen, ModuleBase::Matrix3* symop, int &nop) const { - int m1=0; - int m2=0; - int n=0; + int m1 = 0; + int m2 = 0; + int n = 0; + // allocate memory for the symmetry operations ModuleBase::Matrix3 iden(1,0,0,0,1,0,0,0,1); ModuleBase::Matrix3 sig(1,0,0,0,1,0,0,0,1); ModuleBase::Matrix3 temp1(1,0,0,0,1,0,0,0,1); ModuleBase::Matrix3 temp2(1,0,0,0,1,0,0,0,1); - bool flag = 0; + bool flag = false; // mark whether the symmetry operation is a new one int order = 0; int now = 0; - symop[0] = iden; //identity (the trivial element) - nop = 1; + symop[0] = iden; //identity (the trivial element) + nop = 1; // counter of the symmetry operations - //take all generators - for(int i = 0; i < ngen; ++i) + // take all generators + for (int i = 0; i < ngen; ++i) { sig = symgen[i]; - flag = 1; - for(int j = 0; j < nop; ++j) + flag = true; // assume it is a new symmetry operation + // search if the symmetry operation already exists among the found symmetry operations + // if so, skip it + for (int j = 0; j < nop; ++j) { - if(symop[j] == sig) + if (symop[j] == sig) { - flag = 0; + flag = 0; // not a new symmetry operation break; } } - if(flag == 0) + if (flag == 0) // if old, return { continue; } + // otherwise - //Determine the order of the operation + // determine the order of the operation: by which power will the operation return + // to the identity operation. temp1 = sig; - for(int j = 1; j < 100; ++j) + for (int j = 1; j < 100; ++j) { order = j; - if(temp1 == iden) + if (temp1 == iden) { break; } - temp1= sig * temp1; + temp1 = sig * temp1; } now = nop; - for(int j = 0; j < nop; ++j) + for (int j = 0; j < nop; ++j) { temp1 = symop[j]; - for(int k = 1; k < order; ++k) + for (int k = 1; k < order; ++k) { temp1 = sig * temp1; - for(int l = 0; l < nop; ++l) + for (int l = 0; l < nop; ++l) { temp2 = symop[l] * temp1; flag = 1; - for(int m = 0; m < now; ++m) + for (int m = 0; m < now; ++m) { - if(symop[m] == temp2) + if (symop[m] == temp2) { flag = 0; break; } } - if(flag == 0) + if (flag == 0) { - continue; //the new-found element has already existed. + continue; //the newly-found element has already existed. } - ++now; //the number of elements we found - if(now > 48) + ++now; // the number of elements we found + if (now > 48) // number of symm_op cannot be more than 48 (of O_h point group) { - std::cout<<"\n a: now= "< 48) + if (now > 48) { - std::cout<<"\n b: now= "< 1) { + if(PARAM.inp.cal_symm_repr[0] > 1) { ModuleBase::TITLE("Symmetry_Basic", "setgroup"); } - ModuleBase::Matrix3 symgen[3]; - - ModuleBase::Matrix3 inv(-1, 0, 0, 0,-1, 0, 0, 0,-1); - ModuleBase::Matrix3 r3d( 0, 1, 0, 0, 0, 1, 1, 0, 0); - ModuleBase::Matrix3 r6z( 1, 1, 0,-1, 0, 0, 0, 0, 1); - ModuleBase::Matrix3 r2hex( 1, 0, 0,-1,-1, 0, 0, 0,-1); - ModuleBase::Matrix3 r2tri(-1, 0, 0, 0, 0,-1, 0,-1, 0); - ModuleBase::Matrix3 r4zp( 0, 1, 0,-1, 0, 0, 0, 0, 1); - ModuleBase::Matrix3 r2yp(-1, 0, 0, 0, 1, 0, 0, 0,-1); - ModuleBase::Matrix3 r4zbc( 0, 0,-1, 1, 1, 1, 0,-1, 0); - ModuleBase::Matrix3 r4zfc( 1, 0,-1, 1, 0, 0, 1,-1, 0); - ModuleBase::Matrix3 r2zp(-1, 0, 0, 0,-1, 0, 0, 0, 1); - ModuleBase::Matrix3 r2ybc( 0, 0, 1,-1,-1,-1, 1, 0, 0); - ModuleBase::Matrix3 r2zbc( 0, 1, 0, 1, 0, 0,-1,-1,-1); - ModuleBase::Matrix3 r2ybas( 0,-1, 0,-1, 0, 0, 0, 0,-1); - ModuleBase::Matrix3 r2yfc( 0,-1, 1, 0,-1, 0, 1,-1, 0); - ModuleBase::Matrix3 r2zfc( 0, 1,-1, 1, 0,-1, 0, 0,-1); + ModuleBase::Matrix3 symgen[3]; // the number of generators is up to 3 + + ModuleBase::Matrix3 inv(-1, 0, 0, 0,-1, 0, 0, 0,-1); // (x, y, z) -> (-x, -y, -z) + ModuleBase::Matrix3 r3d( 0, 1, 0, 0, 0, 1, 1, 0, 0); // (x, y, z) -> (y, z, x) + ModuleBase::Matrix3 r6z( 1, 1, 0,-1, 0, 0, 0, 0, 1); // (x, y, z) -> (x+y, -x, z) + ModuleBase::Matrix3 r2hex( 1, 0, 0,-1,-1, 0, 0, 0,-1); // (x, y, z) -> (x, -x-y, -z) + ModuleBase::Matrix3 r2tri(-1, 0, 0, 0, 0,-1, 0,-1, 0); // (x, y, z) -> (-x, -z, -y) + ModuleBase::Matrix3 r4zp( 0, 1, 0,-1, 0, 0, 0, 0, 1); // (x, y, z) -> (y, -x, z) + ModuleBase::Matrix3 r2yp(-1, 0, 0, 0, 1, 0, 0, 0,-1); // (x, y, z) -> (-x, y, -z) + ModuleBase::Matrix3 r4zbc( 0, 0,-1, 1, 1, 1, 0,-1, 0); // (x, y, z) -> (-z, x+y+z, -y) + ModuleBase::Matrix3 r4zfc( 1, 0,-1, 1, 0, 0, 1,-1, 0); // (x, y, z) -> (x-z, x, x-y) + ModuleBase::Matrix3 r2zp(-1, 0, 0, 0,-1, 0, 0, 0, 1); // (x, y, z) -> (-x, -y, z) + ModuleBase::Matrix3 r2ybc( 0, 0, 1,-1,-1,-1, 1, 0, 0); // (x, y, z) -> (z, -x-y-z, x) + ModuleBase::Matrix3 r2zbc( 0, 1, 0, 1, 0, 0,-1,-1,-1); // (x, y, z) -> (y, x, -x-y-z) + ModuleBase::Matrix3 r2ybas( 0,-1, 0,-1, 0, 0, 0, 0,-1); // (x, y, z) -> (-y, -x, -z) + ModuleBase::Matrix3 r2yfc( 0,-1, 1, 0,-1, 0, 1,-1, 0); // (x, y, z) -> (-y+z, -y, x-y) + ModuleBase::Matrix3 r2zfc( 0, 1,-1, 1, 0,-1, 0, 0,-1); // (x, y, z) -> (y-z, x-z, -z) //the pure translation lattice (bravais lattice) has some maximum symmetry //set first up the point group operations for this symmetry. symgen[0] = inv; - - if(ibrav == 14) - { - this->matrigen(symgen, 1, symop, nop); - } - else if(ibrav == 13) - { - symgen[1] = r2ybas; - this->matrigen(symgen, 2, symop, nop); - } - else if(ibrav == 12) - { - symgen[1] = r2yp; - this->matrigen(symgen, 2, symop, nop); - } - else if(ibrav == 11) - { - symgen[1] = r2zp; - symgen[2] = r2ybas; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 10) - { - symgen[1] = r2zfc; - symgen[2] = r2yfc; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 9) - { - symgen[1] = r2zbc; - symgen[2] = r2ybc; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 8) - { - symgen[1] = r2zp; - symgen[2] = r2yp; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 7) - { - symgen[1] = r2tri; - symgen[2] = r3d; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 6) - { - symgen[1] = r4zbc; - symgen[2] = r2ybc; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 5) - { - symgen[1] = r4zp; - symgen[2] = r2yp; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 4) - { - symgen[1] = r6z; - symgen[2] = r2hex; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 3) - { - symgen[1] = r3d; - symgen[2] = r4zfc; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 2) - { - symgen[1] = r3d; - symgen[2] = r4zbc; - this->matrigen(symgen, 3, symop, nop); - } - else if(ibrav == 1) - { - symgen[1] = r3d; - symgen[2] = r4zp; - this->matrigen(symgen, 3, symop, nop); - } - - if(test_brav) + // generate the point group operations for the bravais lattice + // rewrite with switch-case to get better performance and readability + switch (ibrav) { + case 1: + symgen[1] = r3d; + symgen[2] = r4zp; + this->matrigen(symgen, 3, symop, nop); + break; + case 2: + symgen[1] = r3d; + symgen[2] = r4zbc; + this->matrigen(symgen, 3, symop, nop); + break; + case 3: + symgen[1] = r3d; + symgen[2] = r4zfc; + this->matrigen(symgen, 3, symop, nop); + break; + case 4: + symgen[1] = r6z; + symgen[2] = r2hex; + this->matrigen(symgen, 3, symop, nop); + break; + case 5: + symgen[1] = r4zp; + symgen[2] = r2yp; + this->matrigen(symgen, 3, symop, nop); + break; + case 6: + symgen[1] = r4zbc; + symgen[2] = r2ybc; + this->matrigen(symgen, 3, symop, nop); + break; + case 7: + symgen[1] = r2tri; + symgen[2] = r3d; + this->matrigen(symgen, 3, symop, nop); + break; + case 8: + symgen[1] = r2zp; + symgen[2] = r2yp; + this->matrigen(symgen, 3, symop, nop); + break; + case 9: + symgen[1] = r2zbc; + symgen[2] = r2ybc; + this->matrigen(symgen, 3, symop, nop); + break; + case 10: + symgen[1] = r2zfc; + symgen[2] = r2yfc; + this->matrigen(symgen, 3, symop, nop); + break; + case 11: + symgen[1] = r2zp; + symgen[2] = r2ybas; + this->matrigen(symgen, 3, symop, nop); + break; + case 12: + symgen[1] = r2yp; + this->matrigen(symgen, 2, symop, nop); + break; + case 13: + symgen[1] = r2ybas; + this->matrigen(symgen, 2, symop, nop); + break; + case 14: + this->matrigen(symgen, 1, symop, nop); + break; + default: + ModuleBase::WARNING_QUIT("Symmetry_Basic::setgroup", + "ibrav = " + std::to_string(ibrav) + " is not supported."); + break; + } + + // print + if (test_brav) { ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Number of rotation matrices", nop); } // print the symmetry operations - if(PARAM.inp.out_symm_reprmat[0] > 0) - { - GlobalV::ofs_running << " There are " << nop << " symmetry operation representation matrices" - << std::endl; - GlobalV::ofs_running << " data sequence: [e11, e12, e13, e21, e22, e23, e31, e32, e33].reshape(3, 3)" + if (PARAM.inp.cal_symm_repr[0] > 0) + { + GlobalV::ofs_running << std::endl + << " ======================================================================\n" + << " MATRIX REPRESENTATION OF SYMMETRY OPERATION\n" + << " ======================================================================\n" + << " There are " << nop << " symmetry operation representation matrices.\n" + << " For each matrix, the elements are arranged like: \n" + << " [[e11, e12, e13], [e21, e22, e23], [e31, e32, e33]].reshape(3, 3)\n" << std::endl; // control the digits - const int precision = PARAM.inp.out_symm_reprmat[1]; + const int precision = PARAM.inp.cal_symm_repr[1]; const int width = precision + 4; std::string fmtstr = " %" + std::to_string(width) + "." + std::to_string(precision) + "f"; fmtstr += fmtstr + fmtstr + "\n"; // print the symmetry operations std::string mat; - const std::vector op_names = { - - }; for (int i = 0; i < nop; ++i) { - mat = " " - + FmtCore::format("No. %3d", i + 1) + "\n" + mat = " " + FmtCore::format("No. %3d", i + 1) + "\n" + FmtCore::format(fmtstr.c_str(), symop[i].e11, symop[i].e12, symop[i].e13) + FmtCore::format(fmtstr.c_str(), symop[i].e21, symop[i].e22, symop[i].e23) + FmtCore::format(fmtstr.c_str(), symop[i].e31, symop[i].e32, symop[i].e33); GlobalV::ofs_running << mat << std::endl; } + GlobalV::ofs_running << " ======================================================================\n"; } return; @@ -564,7 +564,7 @@ bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber, //there are four trivial cases which could be easily determined //because the number of their elements are exclusive - if (PARAM.inp.out_symm_reprmat[0] > 1) { + if (PARAM.inp.cal_symm_repr[0] > 1) { ModuleBase::TITLE("Symmetry_Basic", "pointgroup"); } diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 5e5d746f90f..2222cd4c19c 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -409,7 +409,7 @@ struct Input_para std::vector out_wfc_re_im = {}; ///< specify the bands to be calculated for real and imaginary parts of wfc bool if_separate_k = false; ///< whether to write partial charge for all k-points to individual files or merge them std::vector out_elf = {0, 3}; ///< output the electron localization function (ELF). 0: no; 1: yes - std::vector out_symm_reprmat = {0, 3}; ///< output the symmetry representation matrix + std::vector cal_symm_repr = {0, 3}; ///< output the symmetry representation matrix // ============== #Parameters (12.Postprocess) =========================== double dos_emin_ev = -15.0; @@ -631,7 +631,7 @@ struct Input_para bool test_stress = false; ///< test the stress. bool test_skip_ewald = false; ///< variables for test only int test_atom_input = false; ///< variables for test_atom_input only - // int test_symmetry = false; ///< variables for test_lattice only + // int test_symmetry = false; ///< variables for test_lattice only int test_wf = 0; ///< variables for test_wf only int test_grid = false; ///< variables for test_grid only int test_charge = false; ///< variables for test_vloc only diff --git a/source/source_io/read_input_item_output.cpp b/source/source_io/read_input_item_output.cpp index b5ccc4ecaaf..09bbb79e66a 100644 --- a/source/source_io/read_input_item_output.cpp +++ b/source/source_io/read_input_item_output.cpp @@ -549,27 +549,27 @@ void ReadInput::item_output() { // recover the functionality of test_symmetry by introducing a new keyword "out_symm_mat" // the "out_symm_mat" keyword will be a - Input_Item item("out_symm_reprmat"); + Input_Item item("cal_symm_repr"); item.annotation = "output matrix representation of symmetry operation into running log file" " > 0 output the matrix representation of symmetry operation " ", the second parameter controls the precision, default is 3."; item.read_value = [](const Input_Item& item, Parameter& para) { size_t count = item.get_size(); - std::vector out_symm_reprmat(count); // create a placeholder vector + std::vector cal_symm_repr(count); // create a placeholder vector std::transform(item.str_values.begin(), item.str_values.end(), // iterators of 1 - out_symm_reprmat.begin(), // iterator of 2 + cal_symm_repr.begin(), // iterator of 2 [](std::string s){ return std::stoi(s); }); // lambda func - // assign non-negative values to para.input.out_symm_reprmat - std::copy(out_symm_reprmat.begin(), out_symm_reprmat.end(), - para.input.out_symm_reprmat.begin()); + // assign non-negative values to para.input.cal_symm_repr + std::copy(cal_symm_repr.begin(), cal_symm_repr.end(), + para.input.cal_symm_repr.begin()); }; item.check_value = [](const Input_Item& item, const Parameter& para) { - if (para.input.out_symm_reprmat[0] < 0 || para.input.out_symm_reprmat[0] > 1) + if (para.input.cal_symm_repr[0] < 0 || para.input.cal_symm_repr[0] > 1) { - ModuleBase::WARNING_QUIT("ReadInput", "out_symm_reprmat should be 0 or 1"); + ModuleBase::WARNING_QUIT("ReadInput", "cal_symm_repr should be 0 or 1"); } }; - sync_intvec(input.out_symm_reprmat, 2, 0); + sync_intvec(input.cal_symm_repr, 2, 0); this->add_item(item); } } From 2a09c15ea619b4d07c50b58766a1fbcf08dbf8a3 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Wed, 16 Jul 2025 18:07:14 +0800 Subject: [PATCH 3/3] update the document --- docs/advanced/input_files/input-main.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 40abf00c3b0..a6ee985dc30 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -8,6 +8,7 @@ - [symmetry](#symmetry) - [symmetry\_prec](#symmetry_prec) - [symmetry\_autoclose](#symmetry_autoclose) + - [cal_symm_repr](#cal_symm_repr) - [kpar](#kpar) - [bndpar](#bndpar) - [latname](#latname) @@ -561,6 +562,11 @@ These variables are used to control general system parameters. - True: automatically set symmetry to 0 and continue running without symmetry analysis - **Default**: True +### cal_symm_repr +- **Type**: Integer [Integer]\(optional\) +- **Description**: Whether to print the matrix representation of symmetry operation to running log file. If the first value is given as 1, then all matrix representations will be printed. The second optional parameter controls the precision (number of digits) to print, default is 3, which is enough for a quick check. +- **Default**: 1 3 + ### kpar - **Type**: Integer