1+ #include " ./symmetry_rotation.h"
2+ namespace ModuleSymmetry
3+ {
4+ std::string mat3_fmt (const ModuleBase::Matrix3& m)
5+ {
6+ auto s = [](auto x) { return std::to_string (x); };
7+ return s (m.e11 ) + " " s (m.e12 ) + " " s (m.e13 ) + " \n "
8+ s (m.e21 ) + " " s (m.e22 ) + " " s (m.e23 ) + " \n "
9+ s (m.e31 ) + " " s (m.e32 ) + " " s (m.e33 );
10+ }
11+
12+ // needs to calculate Ts from l=0 to l=max(l_ao,l_abf) before
13+
14+ void print_symrot_info_R (const Symmetry_rotation& symrot, const Symmetry& symm,
15+ const int lmax_ao, const std::vector<TC >& Rs)
16+ {
17+ ModuleBase::TITLE (" ModuleSymmetry" " print_symrot_info_R"
18+ std::ofstream ofs (PARAM .globalv .global_out_dir + " symrot_R.dat"
19+ // Print the irreducible sector (to be optimized)
20+ ofs << " Number of irreducible sector: " get_irreducible_sector ().size () << std::endl;
21+ ofs << " Lmax of AOs: " " \n "
22+ ofs << " Lmax of ABFs: " abfs_Lmax << " \n "
23+ // print AO rotation matrix T
24+ ofs << " Format:\n "
25+ << " The index of the symmetry operation\n "
26+ << " The rotation matrix of this symmetry operation (3*3)\n "
27+ << " (The translation vector of this symmetry operation)\n "
28+ << " Orbital rotation matrix (T) of each angular momentum with size ((2l + 1) * (2l + 1)) \n\n "
29+ const int lmax = std::max (lmax_ao, symrot.abfs_Lmax );
30+ for (int isym = 0 ;isym < symm.nrotk ;++isym)
31+ {
32+ ofs << isym << " \n " mat3_fmt (symm.gmatrix [isym]) << " \n "
33+ << vec3_fmt (symm.gtrans [isym]) << " \n "
34+ for (int l=0 ;l <= lmax;++l)
35+ {
36+ const int nm = 2 * l + 1 ;
37+ // ofs << "l = " << l << ", nm = " << nm << "\n";
38+ const auto & T_block = symrot.rotmat_Slm [isym][l];
39+ for (int m1 = 0 ;m1 < nm;++m1)
40+ {
41+ for (int m2 = 0 ;m2 < nm;++m2)
42+ {
43+ // note: the order of m in orbitals may be different from increasing
44+ // note: is Ts row- or col-major ?
45+ ofs << T_block (m1, m2);
46+ }
47+ ofs << " \n "
48+ }
49+ }
50+ }
51+ ofs.close ();
52+ }
53+
54+ void print_symrot_info_k (const Symmetry_rotation& symrot, const K_Vectors& kv, const UnitCell& ucell)
55+ {
56+ ModuleBase::TITLE (" Symmetry_rotation" " print_symrot_info_k"
57+ std::ofstream ofs (PARAM .globalv .global_out_dir + " symrot_k.dat"
58+ ofs << " Number of IBZ k-points (k stars): " kstars .size () << std::endl;
59+ ofs << " Format:\n " " The symmetry operation index to the irreducible k-point. For the irreducible k-points, isym=0.\n\n "
60+ << " (The direct coordinate of the original k-point)\n "
61+ << " For each atom: \n "
62+ << " - Original index->transformed index, type and the Lmax\n "
63+ << " - Bloch orbital rotation matrix (M) of the given operation and atom, for each angular momentum\n\n "
64+ for (int istar = 0 ;istar < kv.kstars .size ();++istar)
65+ {
66+ ofs << " Star " 1 << " of IBZ k-point " vec3_fmt (kv.kstars [istar].at (0 )) << " :\n "
67+ for (const auto & isym_kvd : kv.kstars [istar])
68+ {
69+ const int & isym = isym_kvd.first ;
70+ ofs << isym << " \n " vec3_fmt (isym_kvd.second ) << " \n "
71+ for (int iat1 =0 ;iat1 < ucell.nat ;++iat1)
72+ {
73+ const int it = ucell.iat2it [iat1]; // it1=it2
74+ const int lmax = ucell.atoms [it].nwl ;
75+ const int iat2 = ucell.symm .get_rotated_atom (isym, iat1);
76+ const double arg = 2 * ModuleBase::PI * isym_kvd.second * symrot.get_return_lattice (iat1,isym);
77+ std::complex <double >phase_factor = std::complex <double >(std::cos (arg), std::sin (arg));
78+ ofs << " atom " 1 << " -> " 1 << " of type " 1 << " with Lmax= " " \n "
79+ for (int l = 0 ;l < lmax + 1 ;++l)
80+ {
81+ const int nm = 2 * l + 1 ;
82+ const auto & m_block = symrot.rotmat_Slm [isym][l];
83+ for (int m1 = 0 ;m1 < nm;++m1)
84+ {
85+ // const int m1_start = m2 * nm;
86+ for (int m2 = 0 ;m2 < nm;++m2)
87+ {
88+ ofs << phase_factor * m_block (m1, m2); // row-major
89+ }
90+ ofs << " \n "
91+ }
92+ }// end l
93+ } // end iat
94+ } // end (k, op)
95+ ofs << " \n "
96+ } // end star
97+ ofs.close ();
98+ ModuleBase::timer::tick (" Symmetry_rotation" " print_symrot_info_k"
99+ }
100+ }
0 commit comments