Skip to content

Commit 79f092c

Browse files
kirk0830kirk0830
committed
Feature: support output the matrix representation of symmetry operation
1 parent 7c68dfc commit 79f092c

File tree

4 files changed

+142
-107
lines changed

4 files changed

+142
-107
lines changed

‎source/source_cell/k_vector_utils.cpp‎

Lines changed: 3 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -486,15 +486,9 @@ void kvec_ibz_kpoint(K_Vectors& kv,
486486
ucell.atoms,
487487
false,
488488
nullptr);
489-
ModuleBase::Matrix3 b_optlat_new(recip_vec1.x,
490-
recip_vec1.y,
491-
recip_vec1.z,
492-
recip_vec2.x,
493-
recip_vec2.y,
494-
recip_vec2.z,
495-
recip_vec3.x,
496-
recip_vec3.y,
497-
recip_vec3.z);
489+
ModuleBase::Matrix3 b_optlat_new(recip_vec1.x, recip_vec1.y, recip_vec1.z,
490+
recip_vec2.x, recip_vec2.y, recip_vec2.z,
491+
recip_vec3.x, recip_vec3.y, recip_vec3.z);
498492
// set the crystal point-group symmetry operation
499493
symm.setgroup(bsymop, bnop, recip_brav_type);
500494
// transform the above symmetric operation matrices between different coordinate

‎source/source_cell/module_symmetry/symmetry_basic.cpp‎

Lines changed: 109 additions & 95 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
#include "symmetry.h"
22
#include "source_base/mymath.h"
33
#include "source_io/module_parameter/parameter.h"
4+
#include "source_base/formatter.h"
45

56
bool ModuleSymmetry::test_brav = 0;
67

@@ -11,23 +12,24 @@ std::string Symmetry_Basic::get_brav_name(const int ibrav) const
1112
{
1213
switch(ibrav)
1314
{
14-
case 1: return "01. Cubic P (simple)";
15-
case 2: return "02. Cubic I (body-centered)";
16-
case 3: return "03. Cubic F (face-centered)";
17-
case 4: return "04. Hexagonal cell";
18-
case 5: return "05. Tetrogonal P (simple)";
19-
case 6: return "06. Tetrogonal I (body-centered)";
20-
case 7: return "07. Rhombohedral (Trigonal) cell";
21-
case 8: return "08. Orthorhombic P(simple)";
22-
case 9: return "09. Orthorhombic I (body-centered)";
15+
case 1: return "01. Cubic P (simple)";
16+
case 2: return "02. Cubic I (body-centered)";
17+
case 3: return "03. Cubic F (face-centered)";
18+
case 4: return "04. Hexagonal cell";
19+
case 5: return "05. Tetrogonal P (simple)";
20+
case 6: return "06. Tetrogonal I (body-centered)";
21+
case 7: return "07. Rhombohedral (Trigonal) cell";
22+
case 8: return "08. Orthorhombic P(simple)";
23+
case 9: return "09. Orthorhombic I (body-centered)";
2324
case 10: return "10. Orthorhombic F (face-centered)";
2425
case 11: return "11. Orthorhombic C (base-centered)";
2526
case 12: return "12. Monoclinic P (simple)";
2627
case 13: return "13. Monoclinic A (base-center)";
2728
case 14: return "14. Triclinic cell";
2829
case 15: return "wrong !! ";
2930
}
30-
return "Congratulations!You have found a bravais lattice that never existed!";
31+
// return "Congratulations! You have found a bravais lattice that never existed!";
32+
return "Unknown Bravais lattice";
3133
}
3234

3335
// Control the accuracy
@@ -45,15 +47,14 @@ bool Symmetry_Basic::equal(const double &m,const double &n)const
4547
void Symmetry_Basic::check_boundary(double &x)const
4648
{
4749
if(equal(x,-0.5) || equal(x,0.5)) x=-0.5;
48-
return;
4950
}
5051

5152
double Symmetry_Basic::get_translation_vector(const double& x1, const double& x2) const
5253
{
5354
double t=0.0; // "t"ranslation
5455
t = x2 - x1;
5556
t = fmod(t+100.0, 1.0);
56-
if( fabs(t-1) < epsilon * 0.5) t = 0.0;
57+
if( fabs(t-1) < epsilon * 0.5) { t = 0.0; }
5758
return t;
5859
}
5960

@@ -187,8 +188,8 @@ void Symmetry_Basic::matrigen(ModuleBase::Matrix3 *symgen, const int ngen, Modul
187188
int m2=0;
188189
int n=0;
189190

190-
ModuleBase::Matrix3 iden(1,0,0,0,1,0,0,0,1);
191-
ModuleBase::Matrix3 sig(1,0,0,0,1,0,0,0,1);
191+
ModuleBase::Matrix3 iden(1,0,0,0,1,0,0,0,1);
192+
ModuleBase::Matrix3 sig(1,0,0,0,1,0,0,0,1);
192193
ModuleBase::Matrix3 temp1(1,0,0,0,1,0,0,0,1);
193194
ModuleBase::Matrix3 temp2(1,0,0,0,1,0,0,0,1);
194195

@@ -321,25 +322,26 @@ void Symmetry_Basic::matrigen(ModuleBase::Matrix3 *symgen, const int ngen, Modul
321322
//--------------------------------------------------------------
322323
void Symmetry_Basic::setgroup(ModuleBase::Matrix3* symop, int &nop, const int &ibrav) const
323324
{
324-
if(PARAM.inp.test_symmetry) ModuleBase::TITLE("Symmetry_Basic","setgroup");
325-
325+
if(PARAM.inp.out_symm_reprmat[0] > 1) {
326+
ModuleBase::TITLE("Symmetry_Basic", "setgroup");
327+
}
326328
ModuleBase::Matrix3 symgen[3];
327329

328-
ModuleBase::Matrix3 inv(-1,0,0,0,-1,0,0,0,-1);
329-
ModuleBase::Matrix3 r3d(0,1,0,0,0,1,1,0,0);
330-
ModuleBase::Matrix3 r6z(1,1,0,-1,0,0,0,0,1);
331-
ModuleBase::Matrix3 r2hex(1,0,0,-1,-1,0,0,0,-1);
332-
ModuleBase::Matrix3 r2tri(-1,0,0,0,0,-1,0,-1,0);
333-
ModuleBase::Matrix3 r4zp(0,1,0,-1,0,0,0,0,1);
334-
ModuleBase::Matrix3 r2yp(-1,0,0,0,1,0,0,0,-1);
335-
ModuleBase::Matrix3 r4zbc(0,0,-1,1,1,1,0,-1,0);
336-
ModuleBase::Matrix3 r4zfc(1,0,-1,1,0,0,1,-1,0);
337-
ModuleBase::Matrix3 r2zp(-1,0,0,0,-1,0,0,0,1);
338-
ModuleBase::Matrix3 r2ybc(0,0,1,-1,-1,-1,1,0,0);
339-
ModuleBase::Matrix3 r2zbc(0,1,0,1,0,0,-1,-1,-1);
340-
ModuleBase::Matrix3 r2ybas(0,-1,0,-1,0,0,0,0,-1);
341-
ModuleBase::Matrix3 r2yfc(0,-1,1,0,-1,0,1,-1,0);
342-
ModuleBase::Matrix3 r2zfc(0,1,-1,1,0,-1,0,0,-1);
330+
ModuleBase::Matrix3 inv(-1, 0, 0, 0,-1, 0, 0, 0,-1);
331+
ModuleBase::Matrix3 r3d( 0, 1, 0, 0, 0, 1, 1, 0, 0);
332+
ModuleBase::Matrix3 r6z( 1, 1, 0,-1, 0, 0, 0, 0, 1);
333+
ModuleBase::Matrix3 r2hex( 1, 0, 0,-1,-1, 0, 0, 0,-1);
334+
ModuleBase::Matrix3 r2tri(-1, 0, 0, 0, 0,-1, 0,-1, 0);
335+
ModuleBase::Matrix3 r4zp( 0, 1, 0,-1, 0, 0, 0, 0, 1);
336+
ModuleBase::Matrix3 r2yp(-1, 0, 0, 0, 1, 0, 0, 0,-1);
337+
ModuleBase::Matrix3 r4zbc( 0, 0,-1, 1, 1, 1, 0,-1, 0);
338+
ModuleBase::Matrix3 r4zfc( 1, 0,-1, 1, 0, 0, 1,-1, 0);
339+
ModuleBase::Matrix3 r2zp(-1, 0, 0, 0,-1, 0, 0, 0, 1);
340+
ModuleBase::Matrix3 r2ybc( 0, 0, 1,-1,-1,-1, 1, 0, 0);
341+
ModuleBase::Matrix3 r2zbc( 0, 1, 0, 1, 0, 0,-1,-1,-1);
342+
ModuleBase::Matrix3 r2ybas( 0,-1, 0,-1, 0, 0, 0, 0,-1);
343+
ModuleBase::Matrix3 r2yfc( 0,-1, 1, 0,-1, 0, 1,-1, 0);
344+
ModuleBase::Matrix3 r2zfc( 0, 1,-1, 1, 0,-1, 0, 0,-1);
343345

344346
//the pure translation lattice (bravais lattice) has some maximum symmetry
345347
//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
428430

429431
if(test_brav)
430432
{
431-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Number of rotation matrices",nop);
433+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Number of rotation matrices", nop);
432434
}
433435

434-
if(PARAM.inp.test_symmetry > 1)
436+
// print the symmetry operations
437+
if(PARAM.inp.out_symm_reprmat[0] > 0)
435438
{
436-
GlobalV::ofs_running<<" THERE ARE " << nop << " ROTATION MATRICES FOR THE PURE BRAVAIS LATTICE"<<std::endl;
437-
GlobalV::ofs_running<<" E11 E12 E13 E21 E22 E23 E31 E32 E33"<<std::endl;
438-
for(int i = 0; i < nop; ++i)
439+
GlobalV::ofs_running << " There are " << nop << " symmetry operation representation matrices"
440+
<< std::endl;
441+
GlobalV::ofs_running << " data sequence: [e11, e12, e13, e21, e22, e23, e31, e32, e33].reshape(3, 3)"
442+
<< std::endl;
443+
444+
// control the digits
445+
const int precision = PARAM.inp.out_symm_reprmat[1];
446+
const int width = precision + 4;
447+
std::string fmtstr = " %" + std::to_string(width) + "." + std::to_string(precision) + "f";
448+
fmtstr += fmtstr + fmtstr + "\n";
449+
450+
// print the symmetry operations
451+
std::string mat;
452+
const std::vector<std::string> op_names = {
453+
454+
};
455+
for (int i = 0; i < nop; ++i)
439456
{
440-
GlobalV::ofs_running << " "
441-
<< std::setw(3) << i + 1
442-
<< std::setw(4) << symop[i].e11
443-
<< std::setw(4) << symop[i].e12
444-
<< std::setw(4) << symop[i].e13
445-
<< std::setw(4) << symop[i].e21
446-
<< std::setw(4) << symop[i].e22
447-
<< std::setw(4) << symop[i].e23
448-
<< std::setw(4) << symop[i].e31
449-
<< std::setw(4) << symop[i].e32
450-
<< std::setw(4) << symop[i].e33 << std::endl;
457+
mat = " "
458+
+ FmtCore::format("No. %3d", i + 1) + "\n"
459+
+ FmtCore::format(fmtstr.c_str(), symop[i].e11, symop[i].e12, symop[i].e13)
460+
+ FmtCore::format(fmtstr.c_str(), symop[i].e21, symop[i].e22, symop[i].e23)
461+
+ FmtCore::format(fmtstr.c_str(), symop[i].e31, symop[i].e32, symop[i].e33);
462+
GlobalV::ofs_running << mat << std::endl;
451463
}
452464
}
453465

@@ -462,75 +474,75 @@ int Symmetry_Basic::subgroup(const int& nrot, const int& ninv,
462474
{
463475
if (ninv)
464476
{
465-
// if (nc2 >= 7 && nc3 >= 2 && nc6 >= 2 && ns1 >= 7 && ns3 >= 2 && ns6 >= 2) return 27;//D_6h
466-
if (nc2 >= 3 && nc3 >= 8 && ns1 >= 3 && ns6 >= 8) return 29; //T_h
477+
// if (nc2 >= 7 && nc3 >= 2 && nc6 >= 2 && ns1 >= 7 && ns3 >= 2 && ns6 >= 2) { return 27; } //D_6h
478+
if (nc2 >= 3 && nc3 >= 8 && ns1 >= 3 && ns6 >= 8) { return 29; } //T_h
467479
}
468480
else
469481
{
470-
if (nc2 >= 9 && nc3 >= 8 && nc4 >= 6) return 30; //O
471-
if (nc2 >= 3 && nc3 >= 8 && ns1 >= 6 && ns4 >= 6) return 31; //T_d
482+
if (nc2 >= 9 && nc3 >= 8 && nc4 >= 6) { return 30; } //O
483+
if (nc2 >= 3 && nc3 >= 8 && ns1 >= 6 && ns4 >= 6) { return 31; } //T_d
472484
}
473485
}
474486
if (nrot > 16)//not else if: nrot>24 can also fall in this part and below
475487
{
476-
if (ninv && nc2 >= 5 && nc4 >= 2 && ns1 >= 5 && ns4 >= 2) return 20; //D_4h
488+
if (ninv && nc2 >= 5 && nc4 >= 2 && ns1 >= 5 && ns4 >= 2) { return 20; } //D_4h
477489
}
478490
if (nrot > 12)
479491
{
480492
if (ninv)
481493
{
482-
if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2 && ns1 >= 1 && ns3 >= 2 && ns6 >= 2) return 23;//C_6h
483-
if (nc2 >= 3 && nc3 >= 2 && ns1 >= 3 && ns6 >= 2)return 13;//D_3d
494+
if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2 && ns1 >= 1 && ns3 >= 2 && ns6 >= 2) { return 23; } //C_6h
495+
if (nc2 >= 3 && nc3 >= 2 && ns1 >= 3 && ns6 >= 2) { return 13; } //D_3d
484496
}
485497
else
486498
{
487-
if (nc2 >= 3 && nc3 >= 8)return 28; //T
488-
if (nc2 >= 3 && nc3 >= 2 && ns1 >= 4 && ns3 >= 2) return 26;//D_3h
489-
if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2 && ns1 >= 6) return 25;//C_6v
490-
if (nc2 >= 7 && nc3 >= 2 && nc6 >= 2) return 24;//D_6
499+
if (nc2 >= 3 && nc3 >= 8) { return 28; } //T
500+
if (nc2 >= 3 && nc3 >= 2 && ns1 >= 4 && ns3 >= 2) { return 26; } //D_3h
501+
if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2 && ns1 >= 6) { return 25; } //C_6v
502+
if (nc2 >= 7 && nc3 >= 2 && nc6 >= 2) { return 24; } //D_6
491503
}
492504
}
493505
if (nrot > 8)
494506
{
495507
if (ninv)
496508
{
497-
if (nc2 >= 1 && nc4 >= 2 && ns1 >= 1 && ns4 >= 2) return 16;//C_4h
498-
if (nc2 >= 3 && ns1 >= 3)return 8;//D_2h
509+
if (nc2 >= 1 && nc4 >= 2 && ns1 >= 1 && ns4 >= 2) { return 16; } //C_4h
510+
if (nc2 >= 3 && ns1 >= 3) { return 8; } //D_2h
499511
}
500512
else
501513
{
502-
if (nc2 >= 3 && ns1 >= 2 && ns4 >= 2)return 19;//D_2d
503-
if (nc2 >= 1 && nc4 >= 2 && ns1 >= 4) return 18;//C_4v
504-
if (nc2 >= 5 && nc4 >= 2)return 17;//D_4
514+
if (nc2 >= 3 && ns1 >= 2 && ns4 >= 2) { return 19; } //D_2d
515+
if (nc2 >= 1 && nc4 >= 2 && ns1 >= 4) { return 18; } //C_4v
516+
if (nc2 >= 5 && nc4 >= 2) { return 17; } //D_4
505517
}
506518
}
507519
if (nrot > 6)
508520
{
509-
if (nc3 >= 2 && ns1 >= 1 && ns3 >= 2)return 22;//C_3h
510-
if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2)return 21;//C_6
511-
if (nc3 >= 2 && ns1 >= 3)return 12;//C_3v
512-
if (nc2 >= 3 && nc3 >= 2)return 11;//D_3
513-
if (ninv && nc3 >= 2 && ns3 >= 2)return 10;//S_6
521+
if (nc3 >= 2 && ns1 >= 1 && ns3 >= 2) { return 22; } //C_3h
522+
if (nc2 >= 1 && nc3 >= 2 && nc6 >= 2) { return 21; } //C_6
523+
if (nc3 >= 2 && ns1 >= 3) { return 12; } //C_3v
524+
if (nc2 >= 3 && nc3 >= 2) { return 11; } //D_3
525+
if (ninv && nc3 >= 2 && ns3 >= 2) { return 10; }//S_6
514526
}
515527
if (nrot > 4)
516528
{
517-
if (nc2 >= 1 && ns4 >= 2)return 15;//S_4
518-
if (nc2 >= 1 && nc4 >= 2)return 14;//C_4
519-
if (nc2 >= 1 && ns1 >= 2)return 7;//C_2v
520-
if (nc2 >= 3)return 6;//D_2
521-
if (ninv && nc2 >= 1 && ns1 >= 1)return 5;//C_2h
529+
if (nc2 >= 1 && ns4 >= 2) { return 15; } //S_4
530+
if (nc2 >= 1 && nc4 >= 2) { return 14; } //C_4
531+
if (nc2 >= 1 && ns1 >= 2) { return 7; } //C_2v
532+
if (nc2 >= 3) { return 6; } //D_2
533+
if (ninv && nc2 >= 1 && ns1 >= 1) { return 5; } //C_2h
522534
}
523535
if (nrot > 3)
524536
{
525-
if (nc3 >= 2)return 9;//C_3
537+
if (nc3 >= 2) { return 9; } //C_3
526538
}
527539
if (nrot > 2)
528540
{
529-
if (ns1 >= 1)return 4;//C_1h
530-
if (nc2 >= 1)return 3;//C_2
531-
if (ninv)return 2;//S_2
541+
if (ns1 >= 1) { return 4; } //C_1h
542+
if (nc2 >= 1) { return 3; } //C_2
543+
if (ninv) { return 2; } //S_2
532544
}
533-
return 1;//C_1
545+
return 1; //C_1
534546
}
535547

536548

@@ -552,7 +564,9 @@ bool Symmetry_Basic::pointgroup(const int& nrot, int& pgnumber,
552564

553565
//there are four trivial cases which could be easily determined
554566
//because the number of their elements are exclusive
555-
if (PARAM.inp.test_symmetry) ModuleBase::TITLE("Symmetry_Basic", "pointgroup");
567+
if (PARAM.inp.out_symm_reprmat[0] > 1) {
568+
ModuleBase::TITLE("Symmetry_Basic", "pointgroup");
569+
}
556570

557571
std::vector<std::string> pgdict = { "none", "C_1", "S_2", "C_2", "C_1h", "C_2h",
558572
"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,
628642
continue;
629643
}
630644

631-
if(trace == -1 && det == 1) ++nc2;
632-
else if(trace == 0 && det == 1) ++nc3;
633-
else if(trace == 1 && det == 1) ++nc4;
634-
else if(trace == 2 && det == 1) ++nc6;
635-
else if(trace == 1 && det == -1) ++ns1;
636-
else if(trace == 0 && det == -1) ++ns6; //mohan add 2012-01-15
637-
else if(trace == -1 && det == -1) ++ns4;
638-
else if(trace == -2 && det == -1) ++ns3; //mohan add 2012-01-15
645+
if(trace == -1 && det == 1) { ++nc2; }
646+
else if(trace == 0 && det == 1) { ++nc3; }
647+
else if(trace == 1 && det == 1) { ++nc4; }
648+
else if(trace == 2 && det == 1) { ++nc6; }
649+
else if(trace == 1 && det == -1) { ++ns1; }
650+
else if(trace == 0 && det == -1) { ++ns6; } //mohan add 2012-01-15
651+
else if(trace == -1 && det == -1) { ++ns4; }
652+
else if(trace == -2 && det == -1) { ++ns3; } //mohan add 2012-01-15
639653
}
640654

641655
if(test_brav)
642656
{
643-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C2",nc2);
644-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C3",nc3);
645-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C4",nc4);
646-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"C6",nc6);
647-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"S1",ns1);
648-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"S3",ns3);
649-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"S4",ns4);
650-
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"S6",ns6);
657+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "C2", nc2);
658+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "C3", nc3);
659+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "C4", nc4);
660+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "C6", nc6);
661+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "S1", ns1);
662+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "S3", ns3);
663+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "S4", ns4);
664+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "S6", ns6);
651665
}
652666

653667
if(nrot == 2)

‎source/source_io/module_parameter/input_parameter.h‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -407,8 +407,9 @@ struct Input_para
407407
std::vector<int> out_pchg = {}; ///< specify the bands to be calculated for partial charge
408408
std::vector<int> out_wfc_norm = {}; ///< specify the bands to be calculated for norm of wfc
409409
std::vector<int> out_wfc_re_im = {}; ///< specify the bands to be calculated for real and imaginary parts of wfc
410-
bool if_separate_k = false; ///< whether to write partial charge for all k-points to individual files or merge them
411-
std::vector<int> out_elf = {0, 3}; ///< output the electron localization function (ELF). 0: no; 1: yes
410+
bool if_separate_k = false; ///< whether to write partial charge for all k-points to individual files or merge them
411+
std::vector<int> out_elf = {0, 3}; ///< output the electron localization function (ELF). 0: no; 1: yes
412+
std::vector<int> out_symm_reprmat = {0, 3}; ///< output the symmetry representation matrix
412413

413414
// ============== #Parameters (12.Postprocess) ===========================
414415
double dos_emin_ev = -15.0;
@@ -630,7 +631,7 @@ struct Input_para
630631
bool test_stress = false; ///< test the stress.
631632
bool test_skip_ewald = false; ///< variables for test only
632633
int test_atom_input = false; ///< variables for test_atom_input only
633-
int test_symmetry = false; ///< variables for test_lattice only
634+
// int test_symmetry = false; ///< variables for test_lattice only
634635
int test_wf = 0; ///< variables for test_wf only
635636
int test_grid = false; ///< variables for test_grid only
636637
int test_charge = false; ///< variables for test_vloc only

0 commit comments

Comments
 (0)