Skip to content

Commit d59b860

Browse files
YuLiu98hongriTianqi
YuLiu98
and
hongriTianqi
authored
Refactor: read pseudopotentials in esolver_fp (#1926)
* Refactor: read pseudopotentials in esolver_fp * fix: fix calculation get_S * fix: split setup_cell and read_pseudo --------- Co-authored-by: Tianqi Zhao <[email protected]>
1 parent ae7b9e1 commit d59b860

File tree

5 files changed

+151
-143
lines changed

5 files changed

+151
-143
lines changed

‎source/module_cell/unitcell.cpp‎

Lines changed: 143 additions & 143 deletions
Original file line numberDiff line numberDiff line change
@@ -778,155 +778,160 @@ void UnitCell::setup_cell(
778778
// OUT(log,"lattice center y",latcenter.y);
779779
// OUT(log,"lattice center z",latcenter.z);
780780

781-
// read in non-local pseudopotential and ouput the projectors.
782-
783-
log << "\n\n\n\n";
784-
log << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
785-
log << " | |" << std::endl;
786-
log << " | Reading pseudopotentials files: |" << std::endl;
787-
log << " | The pseudopotential file is in UPF format. The 'NC' indicates that |" << std::endl;
788-
log << " | the type of pseudopotential is 'norm conserving'. Functional of |" << std::endl;
789-
log << " | exchange and correlation is decided by 4 given parameters in UPF |" << std::endl;
790-
log << " | file. We also read in the 'core correction' if there exists. |" << std::endl;
791-
log << " | Also we can read the valence electrons number and the maximal |" << std::endl;
792-
log << " | angular momentum used in this pseudopotential. We also read in the |" << std::endl;
793-
log << " | trail wave function, trail atomic density and local-pseudopotential|" << std::endl;
794-
log << " | on logrithmic grid. The non-local pseudopotential projector is also|" << std::endl;
795-
log << " | read in if there is any. |" << std::endl;
796-
log << " | |" << std::endl;
797-
log << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
798-
log << "\n\n\n\n";
799-
800-
801-
this->read_cell_pseudopots(s_pseudopot_dir, log);
802-
803-
if(GlobalV::MY_RANK == 0 && GlobalV::out_element_info)
804-
{
805-
for(int it=0; it<this->ntype; it++)
806-
{
807-
Atom* atom = &atoms[it];
808-
std::stringstream ss;
809-
ss << GlobalV::global_out_dir << atom->label
810-
<< "/" << atom->label
811-
<< ".NONLOCAL";
812-
std::ofstream ofs(ss.str().c_str());
813-
814-
ofs << "<HEADER>" << std::endl;
815-
ofs << std::setw(10) << atom->label << "\t" << "label" << std::endl;
816-
ofs << std::setw(10) << atom->ncpp.pp_type << "\t" << "pseudopotential type" << std::endl;
817-
ofs << std::setw(10) << atom->ncpp.lmax << "\t" << "lmax" << std::endl;
818-
ofs << "</HEADER>" << std::endl;
819-
820-
ofs << "\n<DIJ>" << std::endl;
821-
ofs << std::setw(10) << atom->ncpp.nbeta << "\t" << "nummber of projectors." << std::endl;
822-
823-
for(int ib=0; ib<atom->ncpp.nbeta; ib++)
824-
{
825-
for(int ib2=0; ib2<atom->ncpp.nbeta; ib2++)
826-
{
827-
ofs<<std::setw(10) << atom->ncpp.lll[ib]
828-
<< " " << atom->ncpp.lll[ib2]
829-
<< " " << atom->ncpp.dion(ib,ib2)<<std::endl;
830-
}
831-
}
832-
ofs << "</DIJ>" << std::endl;
833-
834-
for(int i=0; i<atom->ncpp.nbeta; i++)
835-
{
836-
ofs << "<PP_BETA>" << std::endl;
837-
ofs << std::setw(10) << i << "\t" << "the index of projectors." <<std::endl;
838-
ofs << std::setw(10) << atom->ncpp.lll[i] << "\t" << "the angular momentum." <<std::endl;
839-
840-
// mohan add
841-
// only keep the nonzero part.
842-
int cut_mesh = atom->ncpp.mesh;
843-
for(int j=atom->ncpp.mesh-1; j>=0; --j)
844-
{
845-
if( abs( atom->ncpp.betar(i,j) ) > 1.0e-10 )
846-
{
847-
cut_mesh = j;
848-
break;
849-
}
850-
}
851-
if(cut_mesh %2 == 0) ++cut_mesh;
781+
//===================================
782+
// set index for iat2it, iat2ia
783+
//===================================
784+
this->set_iat2itia();
852785

853-
ofs << std::setw(10) << cut_mesh << "\t" << "the number of mesh points." << std::endl;
786+
return;
787+
}
854788

789+
void UnitCell::read_pseudo(ofstream &ofs)
790+
{
791+
// read in non-local pseudopotential and ouput the projectors.
792+
ofs << "\n\n\n\n";
793+
ofs << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
794+
ofs << " | |" << std::endl;
795+
ofs << " | Reading pseudopotentials files: |" << std::endl;
796+
ofs << " | The pseudopotential file is in UPF format. The 'NC' indicates that |" << std::endl;
797+
ofs << " | the type of pseudopotential is 'norm conserving'. Functional of |" << std::endl;
798+
ofs << " | exchange and correlation is decided by 4 given parameters in UPF |" << std::endl;
799+
ofs << " | file. We also read in the 'core correction' if there exists. |" << std::endl;
800+
ofs << " | Also we can read the valence electrons number and the maximal |" << std::endl;
801+
ofs << " | angular momentum used in this pseudopotential. We also read in the |" << std::endl;
802+
ofs << " | trail wave function, trail atomic density and local-pseudopotential|" << std::endl;
803+
ofs << " | on logrithmic grid. The non-local pseudopotential projector is also|" << std::endl;
804+
ofs << " | read in if there is any. |" << std::endl;
805+
ofs << " | |" << std::endl;
806+
ofs << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
807+
ofs << "\n\n\n\n";
808+
809+
read_cell_pseudopots(GlobalV::global_pseudo_dir, ofs);
810+
811+
if(GlobalV::MY_RANK == 0 && GlobalV::out_element_info)
812+
{
813+
for(int it=0; it<ntype; it++)
814+
{
815+
Atom* atom = &atoms[it];
816+
std::stringstream ss;
817+
ss << GlobalV::global_out_dir << atom->label
818+
<< "/" << atom->label
819+
<< ".NONLOCAL";
820+
std::ofstream ofs(ss.str().c_str());
821+
822+
ofs << "<HEADER>" << std::endl;
823+
ofs << std::setw(10) << atom->label << "\t" << "label" << std::endl;
824+
ofs << std::setw(10) << atom->ncpp.pp_type << "\t" << "pseudopotential type" << std::endl;
825+
ofs << std::setw(10) << atom->ncpp.lmax << "\t" << "lmax" << std::endl;
826+
ofs << "</HEADER>" << std::endl;
827+
828+
ofs << "\n<DIJ>" << std::endl;
829+
ofs << std::setw(10) << atom->ncpp.nbeta << "\t" << "nummber of projectors." << std::endl;
830+
for(int ib=0; ib<atom->ncpp.nbeta; ib++)
831+
{
832+
for(int ib2=0; ib2<atom->ncpp.nbeta; ib2++)
833+
{
834+
ofs << std::setw(10) << atom->ncpp.lll[ib]
835+
<< " " << atom->ncpp.lll[ib2]
836+
<< " " << atom->ncpp.dion(ib,ib2)<<std::endl;
837+
}
838+
}
839+
ofs << "</DIJ>" << std::endl;
855840

856-
for(int j=0; j<cut_mesh; ++j)
857-
{
858-
ofs << std::setw(15) << atom->ncpp.r[j]
859-
<< std::setw(15) << atom->ncpp.betar(i, j)
860-
<< std::setw(15) << atom->ncpp.rab[j] << std::endl;
861-
}
862-
ofs << "</PP_BETA>" << std::endl;
863-
}
841+
for(int i=0; i<atom->ncpp.nbeta; i++)
842+
{
843+
ofs << "<PP_BETA>" << std::endl;
844+
ofs << std::setw(10) << i << "\t" << "the index of projectors." <<std::endl;
845+
ofs << std::setw(10) << atom->ncpp.lll[i] << "\t" << "the angular momentum." <<std::endl;
846+
847+
// mohan add
848+
// only keep the nonzero part.
849+
int cut_mesh = atom->ncpp.mesh;
850+
for(int j=atom->ncpp.mesh-1; j>=0; --j)
851+
{
852+
if( abs( atom->ncpp.betar(i,j) ) > 1.0e-10 )
853+
{
854+
cut_mesh = j;
855+
break;
856+
}
857+
}
858+
if(cut_mesh %2 == 0) ++cut_mesh;
859+
860+
ofs << std::setw(10) << cut_mesh << "\t" << "the number of mesh points." << std::endl;
861+
862+
for(int j=0; j<cut_mesh; ++j)
863+
{
864+
ofs << std::setw(15) << atom->ncpp.r[j]
865+
<< std::setw(15) << atom->ncpp.betar(i, j)
866+
<< std::setw(15) << atom->ncpp.rab[j] << std::endl;
867+
}
868+
ofs << "</PP_BETA>" << std::endl;
869+
}
864870

865-
ofs.close();
866-
}
867-
}
871+
ofs.close();
872+
}
873+
}
868874

869875
#ifdef __MPI
870-
this->bcast_unitcell_pseudo2();
871-
#endif
876+
bcast_unitcell_pseudo2();
877+
#endif
872878

873-
for(int it=0; it<ntype; it++)
874-
{
875-
if(atoms[0].ncpp.xc_func !=atoms[it].ncpp.xc_func)
876-
{
877-
GlobalV::ofs_warning << "\n type " << atoms[0].label << " functional is "
878-
<< atoms[0].ncpp.xc_func;
879-
880-
GlobalV::ofs_warning << "\n type " << atoms[it].label << " functional is "
881-
<< atoms[it].ncpp.xc_func << std::endl;
882-
883-
ModuleBase::WARNING_QUIT("setup_cell","All DFT functional must consistent.");
884-
}
885-
}
879+
for(int it=0; it<ntype; it++)
880+
{
881+
if(atoms[0].ncpp.xc_func !=atoms[it].ncpp.xc_func)
882+
{
883+
GlobalV::ofs_warning << "\n type " << atoms[0].label << " functional is "
884+
<< atoms[0].ncpp.xc_func;
886885

887-
this->check_structure(GlobalV::MIN_DIST_COEF);
886+
GlobalV::ofs_warning << "\n type " << atoms[it].label << " functional is "
887+
<< atoms[it].ncpp.xc_func << std::endl;
888888

889-
// setup the total number of PAOs
890-
this->cal_natomwfc(log);
889+
ModuleBase::WARNING_QUIT("setup_cell","All DFT functional must consistent.");
890+
}
891+
}
891892

892-
// setup GlobalV::NLOCAL
893-
this->cal_nwfc(log);
893+
check_structure(GlobalV::MIN_DIST_COEF);
894894

895-
// Check whether the number of valence is minimum
896-
if(GlobalV::MY_RANK==0)
897-
{
898-
int abtype = 0;
899-
for(int it=0; it<ntype; it++)
900-
{
901-
if(ModuleBase::MinZval.find(atoms[it].ncpp.psd) != ModuleBase::MinZval.end())
902-
{
903-
if(atoms[it].ncpp.zv > ModuleBase::MinZval.at(atoms[it].ncpp.psd))
904-
{
905-
abtype += 1;
906-
if(abtype == 1)
907-
{
908-
std::cout << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"<<std::endl;
909-
log << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"<<std::endl;
910-
}
911-
std::cout<<" Warning: number valence electrons > " << ModuleBase::MinZval.at(atoms[it].ncpp.psd);
912-
std::cout<<" for " << atoms[it].ncpp.psd << ": " << ModuleBase::EleConfig.at(atoms[it].ncpp.psd) << std::endl;
913-
log << " Warning: number valence electrons > " << ModuleBase::MinZval.at(atoms[it].ncpp.psd);
914-
log << " for " << atoms[it].ncpp.psd << ": " << ModuleBase::EleConfig.at(atoms[it].ncpp.psd) << std::endl;
915-
}
916-
}
917-
}
918-
if(abtype>0)
919-
{
920-
std::cout<< " Please make sure the pseudopotential file is what you need"<<std::endl;
921-
std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"<<std::endl;
922-
log << " Please make sure the pseudopential file is what you need"<<std::endl;
923-
log << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%";
924-
ModuleBase::GlobalFunc::OUT(log,"");
925-
}
926-
}
895+
// setup the total number of PAOs
896+
cal_natomwfc(ofs);
927897

928-
this->cal_meshx();
929-
return;
898+
// setup GlobalV::NLOCAL
899+
cal_nwfc(ofs);
900+
901+
// Check whether the number of valence is minimum
902+
if(GlobalV::MY_RANK==0)
903+
{
904+
int abtype = 0;
905+
for(int it=0; it<ntype; it++)
906+
{
907+
if(ModuleBase::MinZval.find(atoms[it].ncpp.psd) != ModuleBase::MinZval.end())
908+
{
909+
if(atoms[it].ncpp.zv > ModuleBase::MinZval.at(atoms[it].ncpp.psd))
910+
{
911+
abtype += 1;
912+
if(abtype == 1)
913+
{
914+
std::cout << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"<<std::endl;
915+
ofs << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"<<std::endl;
916+
}
917+
std::cout<<" Warning: number valence electrons > " << ModuleBase::MinZval.at(atoms[it].ncpp.psd);
918+
std::cout<<" for " << atoms[it].ncpp.psd << ": " << ModuleBase::EleConfig.at(atoms[it].ncpp.psd) << std::endl;
919+
ofs << " Warning: number valence electrons > " << ModuleBase::MinZval.at(atoms[it].ncpp.psd);
920+
ofs << " for " << atoms[it].ncpp.psd << ": " << ModuleBase::EleConfig.at(atoms[it].ncpp.psd) << std::endl;
921+
}
922+
}
923+
}
924+
if(abtype>0)
925+
{
926+
std::cout<< " Please make sure the pseudopotential file is what you need"<<std::endl;
927+
std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"<<std::endl;
928+
ofs << " Please make sure the pseudopential file is what you need"<<std::endl;
929+
ofs << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%";
930+
ModuleBase::GlobalFunc::OUT(ofs,"");
931+
}
932+
}
933+
934+
cal_meshx();
930935
}
931936

932937
void UnitCell::setup_cell_classic(
@@ -1115,14 +1120,9 @@ void UnitCell::cal_nwfc(std::ofstream &log)
11151120
//OUT(GlobalV::ofs_running,"NLOCAL",GlobalV::NLOCAL);
11161121
log << " " << std::setw(40) << "NLOCAL" << " = " << GlobalV::NLOCAL <<std::endl;
11171122
//========================================================
1118-
// (4) set index for iat2it, iat2ia, itia2iat, itiaiw2iwt
1123+
// (4) set index for itia2iat, itiaiw2iwt
11191124
//========================================================
11201125

1121-
this->set_iat2itia();
1122-
1123-
//delete[] iat2ia;
1124-
//this->iat2ia = new int[nat];// bug fix 2009-3-8
1125-
11261126
// mohan add 2010-09-26
11271127
assert(GlobalV::NLOCAL>0);
11281128
delete[] iwt2iat;

‎source/module_cell/unitcell.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -233,6 +233,8 @@ class UnitCell
233233
int read_atom_species(std::ifstream &ifa, std::ofstream &ofs_running); // read in the atom information for each type of atom
234234
bool read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_running, std::ofstream &ofs_warning); // read in atomic positions
235235
#endif
236+
237+
void read_pseudo(ofstream &ofs);
236238
int find_type(const std::string &label);
237239
void print_tau(void)const;
238240
void print_stru_file(const std::string &fn, const int &type=1, const int &level=0)const; // mohan add 2011-03-22

‎source/module_esolver/esolver_fp.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,8 @@ namespace ModuleESolver
1919
}
2020
void ESolver_FP::Init(Input& inp, UnitCell& cell)
2121
{
22+
cell.read_pseudo(GlobalV::ofs_running);
23+
2224
#ifdef __MPI
2325
this->pw_rho->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD);
2426
#endif

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -60,6 +60,8 @@ void ESolver_KS_LCAO::Init(Input& inp, UnitCell& ucell)
6060

6161
if (GlobalV::CALCULATION == "get_S")
6262
{
63+
ucell.read_pseudo(GlobalV::ofs_running);
64+
6365
if (ModuleSymmetry::Symmetry::symm_flag == 1)
6466
{
6567
GlobalC::symm.analy_sys(ucell, GlobalV::ofs_running);

‎source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test_prep.cpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -123,6 +123,8 @@ void test_deepks::setup_cell()
123123
"STRU",
124124
GlobalV::ofs_running);
125125

126+
ucell.read_pseudo(GlobalV::ofs_running);
127+
126128
return;
127129
}
128130

0 commit comments

Comments
 (0)