Skip to content

Commit ddf2d24

Browse files
author
jiyuang
committed
add projected band structure ouput
1 parent 38e13fe commit ddf2d24

File tree

2 files changed

+283
-30
lines changed

2 files changed

+283
-30
lines changed

‎source/module_base/constants.h‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -95,6 +95,15 @@ const double RYDBERG_SI = HARTREE_SI/2.0; //J
9595
// zero up to a given accuracy
9696
//const double epsr = 1.0e-6;
9797
const double threshold_wg = 1.0e-14;
98+
99+
const std::string Name_Angular[5][11] =
100+
{
101+
{"s"},
102+
{"px", "py", "pz"},
103+
{"d3z^2-r^2", "dxy", "dxz", "dx^2-y^2", "dyz"},
104+
{"f5z^2-3r^2", "f5xz^2-xr^2", "f5yz^2-yr^2", "fzx^2-zy^2", "fxyz", "fx^3-3*xy^2", "f3yx^2-y^3"},
105+
{"g1", "g2", "g3", "g4", "g5", "g6", "g7", "g8", "g9"}
106+
};
98107
}
99108

100109
#endif

‎source/src_io/energy_dos.cpp‎

Lines changed: 274 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -461,35 +461,8 @@ void energy::perform_dos(Local_Orbital_wfc &lowf, LCAO_Hamilt &uhm)
461461
out.close();
462462
}
463463

464-
std::string Name_Angular[5][11];
465464
/* decomposed Mulliken charge */
466465

467-
Name_Angular[0][0] = "s ";
468-
Name_Angular[1][0] = "px ";
469-
Name_Angular[1][1] = "py ";
470-
Name_Angular[1][2] = "pz ";
471-
Name_Angular[2][0] = "d3z^2-r^2 ";
472-
Name_Angular[2][1] = "dxy ";
473-
Name_Angular[2][2] = "dxz ";
474-
Name_Angular[2][3] = "dx^2-y^2 ";
475-
Name_Angular[2][4] = "dyz ";
476-
Name_Angular[3][0] = "f5z^2-3r^2 ";
477-
Name_Angular[3][1] = "f5xz^2-xr^2";
478-
Name_Angular[3][2] = "f5yz^2-yr^2";
479-
Name_Angular[3][3] = "fzx^2-zy^2 ";
480-
Name_Angular[3][4] = "fxyz ";
481-
Name_Angular[3][5] = "fx^3-3*xy^2";
482-
Name_Angular[3][6] = "f3yx^2-y^3 ";
483-
Name_Angular[4][0] = "g1 ";
484-
Name_Angular[4][1] = "g2 ";
485-
Name_Angular[4][2] = "g3 ";
486-
Name_Angular[4][3] = "g4 ";
487-
Name_Angular[4][4] = "g5 ";
488-
Name_Angular[4][5] = "g6 ";
489-
Name_Angular[4][6] = "g7 ";
490-
Name_Angular[4][7] = "g8 ";
491-
Name_Angular[4][8] = "g9 ";
492-
493466
{std::stringstream as;
494467
as << GlobalV::global_out_dir << "PDOS";
495468
std::ofstream out(as.str().c_str());
@@ -580,7 +553,7 @@ void energy::perform_dos(Local_Orbital_wfc &lowf, LCAO_Hamilt &uhm)
580553
const int m1 = atom1->iw2m[j];
581554
out <<std::setw(5) << i << std::setw(8)
582555
<< GlobalC::ucell.atoms[t].label <<std::setw(5)
583-
<<L1<<std::setw(5) <<m1<<std::setw(5)<<N1+1<<std::setw(15)<< Name_Angular[L1][m1] << std::endl;
556+
<<L1<<std::setw(5) <<m1<<std::setw(5)<<N1+1<<std::setw(15)<< ModuleBase::Name_Angular[L1][m1] << std::endl;
584557
}
585558
}
586559
out <<std::endl<<std::endl;
@@ -669,15 +642,286 @@ void energy::perform_dos(Local_Orbital_wfc &lowf, LCAO_Hamilt &uhm)
669642
nks = GlobalC::kv.nkstot/2;
670643
}
671644

645+
ModuleBase::ComplexMatrix weightk;
646+
ModuleBase::matrix weight;
647+
int NUM = 0;
648+
if(GlobalV::GAMMA_ONLY_LOCAL)
649+
{
650+
NUM=GlobalV::NLOCAL*GlobalV::NBANDS*nspin0;
651+
weightk.create(nspin0, GlobalV::NBANDS*GlobalV::NLOCAL, true);
652+
weight.create(nspin0, GlobalV::NBANDS*GlobalV::NLOCAL, true);
653+
}
654+
else
655+
{
656+
NUM=GlobalV::NLOCAL*GlobalV::NBANDS*GlobalC::kv.nks;
657+
weightk.create(GlobalC::kv.nks, GlobalV::NBANDS*GlobalV::NLOCAL, true);
658+
weight.create(GlobalC::kv.nks, GlobalV::NBANDS*GlobalV::NLOCAL, true);
659+
}
660+
672661
for(int is=0; is<nspin0; is++)
673662
{
674663
std::stringstream ss2;
675664
ss2 << GlobalV::global_out_dir << "BANDS_" << is+1 << ".dat";
676665
GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl;
677666
Dos::nscf_band(is, ss2.str(), nks, GlobalV::NBANDS, this->ef*0, GlobalC::wf.ekb);
678-
}
679667

680-
}
668+
// Projeced band structure added by jiyy-2022-4-20
669+
if(GlobalV::GAMMA_ONLY_LOCAL)
670+
{
671+
std::vector<ModuleBase::matrix> Mulk;
672+
Mulk.resize(1);
673+
Mulk[0].create(pv->ncol,pv->nrow);
674+
675+
676+
ModuleBase::matrix Dwf = lowf.wfc_gamma[is];
677+
for (int i=0; i<GlobalV::NBANDS; ++i)
678+
{
679+
const int NB= i+1;
680+
681+
const double one_float=1.0, zero_float=0.0;
682+
const int one_int=1;
683+
684+
#ifdef __MPI
685+
const char T_char='T';
686+
pdgemv_(
687+
&T_char,
688+
&GlobalV::NLOCAL,&GlobalV::NLOCAL,
689+
&one_float,
690+
uhm.LM->Sloc.data(), &one_int, &one_int, pv->desc,
691+
Dwf.c, &one_int, &NB, pv->desc, &one_int,
692+
&zero_float,
693+
Mulk[0].c, &one_int, &NB, pv->desc,
694+
&one_int);
695+
#endif
696+
697+
for (int j=0; j<GlobalV::NLOCAL; ++j)
698+
{
699+
700+
if ( pv->in_this_processor(j,i) )
701+
{
702+
703+
const int ir = pv->trace_loc_row[j];
704+
const int ic = pv->trace_loc_col[i];
705+
weightk(is, i*GlobalV::NLOCAL+j) = Mulk[0](ic,ir)*lowf.wfc_gamma[is](ic,ir);
706+
}
707+
}
708+
}//ib
709+
}//if
710+
else
711+
{
712+
GlobalV::SEARCH_RADIUS = atom_arrange::set_sr_NL(
713+
GlobalV::ofs_running,
714+
GlobalV::OUT_LEVEL,
715+
GlobalC::ORB.get_rcutmax_Phi(),
716+
GlobalC::ucell.infoNL.get_rcutmax_Beta(),
717+
GlobalV::GAMMA_ONLY_LOCAL);
718+
719+
atom_arrange::search(
720+
GlobalV::SEARCH_PBC,
721+
GlobalV::ofs_running,
722+
GlobalC::GridD,
723+
GlobalC::ucell,
724+
GlobalV::SEARCH_RADIUS,
725+
GlobalV::test_atom_input);//qifeng-2019-01-21
726+
727+
uhm.LM->allocate_HS_R(pv->nnr);
728+
uhm.LM->zeros_HSR('S');
729+
uhm.genH.calculate_S_no();
730+
uhm.genH.build_ST_new('S', false, GlobalC::ucell);
731+
std::vector<ModuleBase::ComplexMatrix> Mulk;
732+
Mulk.resize(1);
733+
Mulk[0].create(pv->ncol,pv->nrow);
734+
735+
736+
for(int ik=0;ik<GlobalC::kv.nks;ik++)
737+
{
738+
739+
if(is == GlobalC::kv.isk[ik])
740+
{
741+
uhm.LM->allocate_HS_k(pv->nloc);
742+
uhm.LM->zeros_HSk('S');
743+
uhm.LM->folding_fixedH(ik);
744+
745+
746+
ModuleBase::ComplexMatrix Dwfc = conj(lowf.wfc_k[ik]);
747+
748+
for (int i=0; i<GlobalV::NBANDS; ++i)
749+
{
750+
const int NB= i+1;
751+
752+
const double one_float[2]={1.0, 0.0}, zero_float[2]={0.0, 0.0};
753+
const int one_int=1;
754+
// const int two_int=2;
755+
const char T_char='T'; // N_char='N',U_char='U'
756+
757+
#ifdef __MPI
758+
pzgemv_(
759+
&T_char,
760+
&GlobalV::NLOCAL,&GlobalV::NLOCAL,
761+
&one_float[0],
762+
uhm.LM->Sloc2.data(), &one_int, &one_int, pv->desc,
763+
Dwfc.c, &one_int, &NB, pv->desc, &one_int,
764+
&zero_float[0],
765+
Mulk[0].c, &one_int, &NB, pv->desc,
766+
&one_int);
767+
#endif
768+
769+
770+
for (int j=0; j<GlobalV::NLOCAL; ++j)
771+
{
772+
773+
if ( pv->in_this_processor(j,i) )
774+
{
775+
776+
const int ir = pv->trace_loc_row[j];
777+
const int ic = pv->trace_loc_col[i];
778+
779+
weightk(ik, i*GlobalV::NLOCAL+j) = Mulk[0](ic,ir)*lowf.wfc_k[ik](ic,ir);
780+
}
781+
}
782+
783+
}//ib
784+
785+
}//if
786+
}//ik
787+
#ifdef __MPI
788+
atom_arrange::delete_vector(
789+
GlobalV::ofs_running,
790+
GlobalV::SEARCH_PBC,
791+
GlobalC::GridD,
792+
GlobalC::ucell,
793+
GlobalV::SEARCH_RADIUS,
794+
GlobalV::test_atom_input);
795+
#endif
796+
}//else
797+
#ifdef __MPI
798+
MPI_Reduce(weightk.real().c, weight.c , NUM , MPI_DOUBLE , MPI_SUM, 0, MPI_COMM_WORLD);
799+
#endif
800+
801+
if(GlobalV::MY_RANK == 0)
802+
{
803+
std::stringstream ps2;
804+
ps2 << GlobalV::global_out_dir << "PBANDS_" << is+1;
805+
GlobalV::ofs_running << "\n Output projected bands in file: " << ps2.str() << std::endl;
806+
std::ofstream out(ps2.str().c_str());
807+
808+
out << "<"<<"pband"<<">" <<std::endl;
809+
out << "<"<<"nspin"<<">" << GlobalV::NSPIN<< "<"<<"/"<<"nspin"<<">"<< std::endl;
810+
if (GlobalV::NSPIN==4)
811+
out << "<"<<"norbitals"<<">" <<std::setw(2) <<GlobalV::NLOCAL/2<< "<"<<"/"<<"norbitals"<<">"<< std::endl;
812+
else
813+
out << "<"<<"norbitals"<<">" <<std::setw(2) <<GlobalV::NLOCAL<< "<"<<"/"<<"norbitals"<<">"<< std::endl;
814+
out << "<"<<"band_structure nkpoints="<<"\""<<nks<<"\""<<" nbands="<<"\""<<GlobalV::NBANDS<<"\""<<" units=\"eV\""<<">" <<std::endl;
815+
if(GlobalV::GAMMA_ONLY_LOCAL)
816+
{
817+
for(int ib = 0; ib < GlobalV::NBANDS; ib++)
818+
out << " " << (GlobalC::wf.ekb[is*nks][ib]) * ModuleBase::Ry_to_eV;
819+
out << std::endl;
820+
}
821+
else
822+
{
823+
for(int ik=0;ik<nks;ik++)
824+
{
825+
for(int ib = 0; ib < GlobalV::NBANDS; ib++)
826+
out << " " << (GlobalC::wf.ekb[ik+is*nks][ib]) * ModuleBase::Ry_to_eV;
827+
out << std::endl;
828+
}
829+
}
830+
out << "<"<<"/"<<"band_structure"<<">"<<std::endl;
831+
832+
for (int i=0; i<GlobalC::ucell.nat; i++)
833+
{
834+
int a = GlobalC::ucell.iat2ia[i];
835+
int t = GlobalC::ucell.iat2it[i];
836+
Atom* atom1 = &GlobalC::ucell.atoms[t];
837+
const int s0 = GlobalC::ucell.itiaiw2iwt(t, a, 0);
838+
for(int j=0; j<atom1->nw; ++j)
839+
{
840+
const int L1 = atom1->iw2l[j];
841+
const int N1 = atom1->iw2n[j];
842+
const int m1 = atom1->iw2m[j];
843+
const int w = GlobalC::ucell.itiaiw2iwt(t, a, j);
844+
845+
//out << "<"<<"/"<<"energy"<<"_"<<"values"<<">" <<std::endl;
846+
out << "<"<<"orbital" <<std::endl;
847+
out <<std::setw(6)<< "index"<<"="<<"\""<<std::setw(40) <<w+1<<"\""<<std::endl;
848+
out <<std::setw(5)<< "atom"<<"_"<<"index"<<"="<<"\""<<std::setw(40) <<i+1<<"\""<<std::endl;
849+
out <<std::setw(8)<< "species"<<"="<<"\""<<GlobalC::ucell.atoms[t].label<<"\""<<std::endl;
850+
out<<std::setw(2)<< "l"<<"="<<"\""<<std::setw(40)<<L1<<"\""<<std::endl;
851+
out <<std::setw(2)<< "m"<<"="<<"\""<<std::setw(40)<<m1<<"\""<<std::endl;
852+
out <<std::setw(2)<< "z"<<"="<<"\""<<std::setw(40)<<N1+1<<"\""<<std::endl;
853+
out << ">" <<std::endl;
854+
out << "<"<<"data"<<">" <<std::endl;
855+
if(GlobalV::GAMMA_ONLY_LOCAL)
856+
{
857+
for(int ib=0;ib<GlobalV::NBANDS;ib++)
858+
{
859+
if (GlobalV::NSPIN==1 || GlobalV::NSPIN==2)
860+
out <<std::setw(13)<< weight(is, ib*GlobalV::NLOCAL+w);
861+
else if (GlobalV::NSPIN==4)
862+
{
863+
int w0 = w-s0;
864+
out <<std::setw(13)<< weight(is, ib*GlobalV::NLOCAL+s0+2*w0)+weight(is, ib*GlobalV::NLOCAL+s0+2*w0+1);
865+
}
866+
}
867+
out <<std::endl;
868+
}
869+
else
870+
{
871+
for(int ik=0;ik<nks;ik++)
872+
{
873+
for(int ib=0;ib<GlobalV::NBANDS;ib++)
874+
{
875+
if (GlobalV::NSPIN==1)
876+
out <<std::setw(13)<< weight(ik, ib*GlobalV::NLOCAL+w);
877+
else if(GlobalV::NSPIN==2)
878+
out <<std::setw(13)<< weight(ik+nks*is, ib*GlobalV::NLOCAL+w);
879+
else if (GlobalV::NSPIN==4)
880+
{
881+
int w0 = w-s0;
882+
out <<std::setw(13)<< weight(ik, ib*GlobalV::NLOCAL+s0+2*w0)+weight(ik, ib*GlobalV::NLOCAL+s0+2*w0+1);
883+
}
884+
}
885+
out <<std::endl;
886+
}
887+
}
888+
out << "<"<<"/"<<"data"<<">" <<std::endl;
889+
out << "<"<<"/"<<"orbital"<<">" <<std::endl;
890+
}//j
891+
}//i
892+
893+
out << "<"<<"/"<<"pband"<<">" <<std::endl;
894+
out.close();}
895+
}//is
896+
{ std::stringstream os;
897+
os<<GlobalV::global_out_dir<<"Orbital";
898+
std::ofstream out(os.str().c_str());
899+
out<< std::setw(5)<<"io"<< std::setw(8) <<"spec" <<std::setw(5)<<"l"<<std::setw(5)<<"m"<<std::setw(5)<<"z"<<std::setw(5)<<"sym"<<std::endl;
900+
901+
902+
for (int i=0; i<GlobalC::ucell.nat; i++)
903+
{
904+
int t = GlobalC::ucell.iat2it[i];
905+
Atom* atom1 = &GlobalC::ucell.atoms[t];
906+
for(int j=0; j<atom1->nw; ++j)
907+
{
908+
const int L1 = atom1->iw2l[j];
909+
const int N1 = atom1->iw2n[j];
910+
const int m1 = atom1->iw2m[j];
911+
out <<std::setw(5) << i << std::setw(8)
912+
<< GlobalC::ucell.atoms[t].label <<std::setw(5)
913+
<<L1<<std::setw(5) <<m1<<std::setw(5)<<N1+1<<std::setw(15)<< ModuleBase::Name_Angular[L1][m1] << std::endl;
914+
}
915+
}
916+
out <<std::endl<<std::endl;
917+
out <<std::setw(5)<< "io"<<std::setw(2)<<"="<<std::setw(2)<<"Orbital index in supercell"<<std::endl;
918+
out <<std::setw(5)<< "spec"<<std::setw(2)<<"="<<std::setw(2)<<"Atomic species label"<<std::endl;
919+
out <<std::setw(5)<< "l"<<std::setw(2)<<"="<<std::setw(2)<<"Angular mumentum quantum number"<<std::endl;
920+
out <<std::setw(5)<< "m"<<std::setw(2)<<"="<<std::setw(2)<<"Magnetic quantum number"<<std::endl;
921+
out <<std::setw(5)<< "z"<<std::setw(2)<<"="<<std::setw(2)<<"Zeta index of orbital"<<std::endl;
922+
out <<std::setw(5)<< "sym"<<std::setw(2)<<"="<<std::setw(2)<<"Symmetry name of real orbital"<<std::endl;
923+
out.close();}
924+
}//out_band
681925
return;
682926
}
683927

0 commit comments

Comments
 (0)