diff --git a/source/module_base/global_function.cpp b/source/module_base/global_function.cpp index 30a7d24af5a..452408f88b9 100644 --- a/source/module_base/global_function.cpp +++ b/source/module_base/global_function.cpp @@ -104,7 +104,11 @@ void DONE(std::ofstream &ofs, const std::string &description, const bool only_ra return; } -bool SCAN_BEGIN(std::ifstream &ifs, const std::string &TargetName, const bool restart, const bool ifwarn) + +bool SCAN_BEGIN(std::ifstream &ifs, + const std::string &TargetName, + const bool restart, + const bool ifwarn) { std::string SearchName; bool find = false; @@ -130,6 +134,51 @@ bool SCAN_BEGIN(std::ifstream &ifs, const std::string &TargetName, const bool re return find; } + +bool SCAN_LINE_BEGIN(std::ifstream &ifs, + const std::string &TargetName, + const bool restart, + const bool ifwarn) +{ + bool find = false; + if (restart) + { + ifs.clear(); + ifs.seekg(0); + } + ifs.rdstate(); + + std::string line; + while (std::getline(ifs,line)) + { + //! obtain the first character, should not be # + size_t first_char_pos = line.find_first_not_of(" \t"); + if (first_char_pos != std::string::npos && line[first_char_pos] == '#') + { + continue; + } + + //! search in each line + std::istringstream iss(line); + std::string SearchName; + while (iss >> SearchName) + { + if (SearchName == TargetName) + { + find = true; + //std::cout << " search name = " << SearchName << std::endl; + return find; + } + } + } + + if (!find && ifwarn) + { + GlobalV::ofs_warning << " In SCAN_LINE_BEGIN, can't find: " << TargetName << " block." << std::endl; + } + return find; +} + void SCAN_END(std::ifstream &ifs, const std::string &TargetName, const bool ifwarn) { std::string SearchName; diff --git a/source/module_base/global_function.h b/source/module_base/global_function.h index 115abde3b63..b82dcea2067 100644 --- a/source/module_base/global_function.h +++ b/source/module_base/global_function.h @@ -1,5 +1,5 @@ -#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_BASE_GLOBAL_FUNCTION_H -#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_BASE_GLOBAL_FUNCTION_H +#ifndef GLOBAL_FUNCTION_H +#define GLOBAL_FUNCTION_H #include "blas_connector.h" #include "global_function-func_each_2.h" // Peize Lin add 2016-09-07 @@ -145,13 +145,26 @@ static void READ_VALUE(std::ifstream& ifs, T& v) return; } -bool SCAN_BEGIN(std::ifstream& ifs, const std::string& TargetName, const bool restart = true, const bool ifwarn = true); -// ifwarn: whether to call GlobalV::ofs_warning when the TargetName is not found, used to avoid invalid warning. -// Mohan warning : the last term can't be written as const bool &restart, -// I don't know why. +//------------------------------------------------------------- +//! The `SCAN_BEGIN` function efficiently searches +//! text files for specified keywords +//------------------------------------------------------------- +bool SCAN_BEGIN(std::ifstream& ifs, + const std::string& TargetName, + const bool restart = true, + const bool ifwarn = true); + +//------------------------------------------------------------- +// The `SCAN_LINE_BEGIN` function efficiently searches +// text files for specified keywords while ignoring comment +// lines and whitespace. It skips any line starting with '#' +//------------------------------------------------------------- +bool SCAN_LINE_BEGIN(std::ifstream& ifs, + const std::string& TargetName, + const bool restart = true, + const bool ifwarn = true); void SCAN_END(std::ifstream& ifs, const std::string& TargetName, const bool ifwarn = true); -// ifwarn: whether to call GlobalV::ofs_warning when the TargetName is not found, used to avoid invalid warning. template static inline void DCOPY(const T& a, T& b, const int& dim) diff --git a/source/module_cell/read_atom_species.cpp b/source/module_cell/read_atom_species.cpp index a50d4ee6d1e..bee3b0f2534 100644 --- a/source/module_cell/read_atom_species.cpp +++ b/source/module_cell/read_atom_species.cpp @@ -13,15 +13,15 @@ bool read_atom_species(std::ifstream& ifa, UnitCell& ucell) { ModuleBase::TITLE("UnitCell","read_atom_species"); + const int ntype = ucell.ntype; std::string word; //========================================== // read in information of each type of atom //========================================== - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "ATOMIC_SPECIES") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "ATOMIC_SPECIES") ) { - ifa.ignore(300, '\n'); ModuleBase::GlobalFunc::OUT(ofs_running,"ntype",ntype); for (int i = 0;i < ntype;i++) { @@ -51,7 +51,8 @@ bool read_atom_species(std::ifstream& ifa, if (!end && ss >> one_string && one_string[0] != '#') { - if (one_string == "auto" || one_string == "upf" || one_string == "vwr" || one_string == "upf201" || one_string == "blps") + if (one_string == "auto" || one_string == "upf" + || one_string == "vwr" || one_string == "upf201" || one_string == "blps") { ucell.pseudo_type[i] = one_string; } @@ -61,7 +62,8 @@ bool read_atom_species(std::ifstream& ifa, } else { - GlobalV::ofs_warning << "unrecongnized pseudopotential type: " << one_string << ", check your STRU file." << std::endl; + GlobalV::ofs_warning << "unrecongnized pseudopotential type: " + << one_string << ", check your STRU file." << std::endl; ModuleBase::WARNING_QUIT("read_atom_species", "unrecongnized pseudo type."); } } @@ -69,7 +71,8 @@ bool read_atom_species(std::ifstream& ifa, // Peize Lin test for bsse 2021.04.07 const std::string bsse_label = "empty"; ucell.atoms[i].flag_empty_element = - (search( ucell.atom_label[i].begin(), ucell.atom_label[i].end(), bsse_label.begin(), bsse_label.end() ) != ucell.atom_label[i].end()) + (search( ucell.atom_label[i].begin(), ucell.atom_label[i].end(), + bsse_label.begin(), bsse_label.end() ) != ucell.atom_label[i].end()) ? true : false; } } @@ -80,7 +83,7 @@ bool read_atom_species(std::ifstream& ifa, ||((PARAM.inp.basis_type == "pw")&&(PARAM.inp.init_wfc.substr(0, 3) == "nao")) || PARAM.inp.onsite_radius > 0.0) { - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "NUMERICAL_ORBITAL") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "NUMERICAL_ORBITAL") ) { for(int i=0; i> ucell.descriptor_file; } } @@ -105,7 +108,7 @@ bool read_atom_species(std::ifstream& ifa, #ifdef __EXX if( GlobalC::exx_info.info_global.cal_exx || PARAM.inp.rpa ) { - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "ABFS_ORBITAL") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "ABFS_ORBITAL") ) { for(int i=0; i> latvec.e22; ModuleBase::GlobalFunc::READ_VALUE(ifa, latvec.e33); @@ -272,7 +281,7 @@ bool read_lattice_constant(std::ifstream& ifa, latvec.e11 = 0.5; latvec.e12 = 0.0; latvec.e13 = 0.0; latvec.e21 =-0.5; latvec.e22 = 0.0; latvec.e23 = 0.0; latvec.e31 = 0.0; latvec.e32 = 0.0; latvec.e33 = 0.0; - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "LATTICE_PARAMETERS") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "LATTICE_PARAMETERS") ) { ifa >> latvec.e12; latvec.e12 = latvec.e12 / 2.0; @@ -283,7 +292,7 @@ bool read_lattice_constant(std::ifstream& ifa, else if(latName=="fco") {//face-centered orthorhombic, ibrav = 10 double bba = 0.0; double cba = 0.0; - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "LATTICE_PARAMETERS") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "LATTICE_PARAMETERS") ) { ifa >> bba; ModuleBase::GlobalFunc::READ_VALUE(ifa, cba); @@ -296,7 +305,7 @@ bool read_lattice_constant(std::ifstream& ifa, else if(latName=="bco") {//body-centered orthorhombic, ibrav = 11 double bba = 0.0; double cba = 0.0; - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "LATTICE_PARAMETERS") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "LATTICE_PARAMETERS") ) { ifa >> bba; ModuleBase::GlobalFunc::READ_VALUE(ifa, cba); @@ -311,7 +320,7 @@ bool read_lattice_constant(std::ifstream& ifa, double bba = 0.0; double cba = 0.0; double cosab = 0.0; double e21 = 0.0; double e22 = 0.0; - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "LATTICE_PARAMETERS") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "LATTICE_PARAMETERS") ) { ifa >> bba >> cba; ModuleBase::GlobalFunc::READ_VALUE(ifa, cosab); @@ -327,7 +336,7 @@ bool read_lattice_constant(std::ifstream& ifa, double bba = 0.0; double cba = 0.0; double cosab = 0.0; double e21 = 0.0; double e22 = 0.0; - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "LATTICE_PARAMETERS") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "LATTICE_PARAMETERS") ) { ifa >> bba >> cba; ModuleBase::GlobalFunc::READ_VALUE(ifa, cosab); @@ -348,7 +357,7 @@ bool read_lattice_constant(std::ifstream& ifa, double cosbc = 0.0; double sinab = 0.0; double term = 0.0; - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifa, "LATTICE_PARAMETERS") ) + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifa, "LATTICE_PARAMETERS") ) { ifa >> bba >> cba >> cosab >> cosac; ModuleBase::GlobalFunc::READ_VALUE(ifa, cosbc); @@ -363,7 +372,6 @@ bool read_lattice_constant(std::ifstream& ifa, } else { - std::cout << "latname is : " << latName << std::endl; ModuleBase::WARNING_QUIT("unitcell::read_lattice_constant","latname not supported!"); } } @@ -383,4 +391,4 @@ bool read_lattice_constant(std::ifstream& ifa, return true; } -} \ No newline at end of file +} diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp index 78f825281f2..b7b891dabc7 100644 --- a/source/module_cell/read_atoms.cpp +++ b/source/module_cell/read_atoms.cpp @@ -17,19 +17,21 @@ #ifdef __LCAO #include "module_basis/module_ao/ORB_read.h" // to use 'ORB' -- mohan 2021-01-30 #endif -namespace unitcell -{ -bool read_atom_positions(UnitCell& ucell, + +bool unitcell::read_atom_positions(UnitCell& ucell, std::ifstream &ifpos, std::ofstream &ofs_running, std::ofstream &ofs_warning) { ModuleBase::TITLE("UnitCell","read_atom_positions"); + std::string& Coordinate = ucell.Coordinate; const int ntype = ucell.ntype; - if( ModuleBase::GlobalFunc::SCAN_BEGIN(ifpos, "ATOMIC_POSITIONS")) + + if( ModuleBase::GlobalFunc::SCAN_LINE_BEGIN(ifpos, "ATOMIC_POSITIONS")) { ModuleBase::GlobalFunc::READ_VALUE(ifpos, Coordinate); + if(Coordinate != "Cartesian" && Coordinate != "Direct" && Coordinate != "Cartesian_angstrom" @@ -70,6 +72,7 @@ bool read_atom_positions(UnitCell& ucell, // start magnetization //======================================= ModuleBase::GlobalFunc::READ_VALUE(ifpos, ucell.atoms[it].label); + if(ucell.atoms[it].label != ucell.atom_label[it]) { ofs_warning << " Label orders in ATOMIC_POSITIONS and ATOMIC_SPECIES sections do not match!" << std::endl; @@ -105,7 +108,6 @@ bool read_atom_positions(UnitCell& ucell, { ucell.atoms[it].nw = 0; ucell.atoms[it].nwl = 2; - //std::cout << ucell.lmaxmax << std::endl; if ( ucell.lmaxmax != 2 ) { ucell.atoms[it].nwl = ucell.lmaxmax; @@ -190,12 +192,14 @@ bool read_atom_positions(UnitCell& ucell, if( (int)tmpid[0] < 0 ) { - std::cout << "read_atom_positions, mismatch in atom number for atom type: " << ucell.atoms[it].label << std::endl; + std::cout << "read_atom_positions, mismatch in atom number for atom type: " + << ucell.atoms[it].label << std::endl; exit(1); } bool input_vec_mag=false; bool input_angle_mag=false; + // read if catch goodbit before "\n" and "#" while ( (tmpid != "\n") && (ifpos.good()) && (tmpid !="#") ) { @@ -237,7 +241,9 @@ bool read_atom_positions(UnitCell& ucell, ifpos.putback(tmp); ifpos >> ucell.atoms[it].m_loc_[ia].y>>ucell.atoms[it].m_loc_[ia].z; ucell.atoms[it].m_loc_[ia].x=tmpamg; - ucell.atoms[it].mag[ia]=sqrt(pow(ucell.atoms[it].m_loc_[ia].x,2)+pow(ucell.atoms[it].m_loc_[ia].y,2)+pow(ucell.atoms[it].m_loc_[ia].z,2)); + ucell.atoms[it].mag[ia]=sqrt(pow(ucell.atoms[it].m_loc_[ia].x,2) + +pow(ucell.atoms[it].m_loc_[ia].y,2) + +pow(ucell.atoms[it].m_loc_[ia].z,2)); input_vec_mag=true; } @@ -246,8 +252,6 @@ bool read_atom_positions(UnitCell& ucell, ifpos.putback(tmp); ucell.atoms[it].mag[ia]=tmpamg; } - - // ucell.atoms[it].mag[ia]; } else if ( tmpid == "angle1") { @@ -315,7 +319,6 @@ bool read_atom_positions(UnitCell& ucell, tmpid = ifpos.get(); } std::string mags; - //cout<<"mag"<> txt; while (ifs >> word) @@ -128,26 +130,30 @@ void Pseudopot_upf::getnameval(std::ifstream& ifs, int& n, std::string* name, st while (1) { pos = txt.find("=", pos); - if (pos == std::string::npos) - break; - pos++; + if (pos == std::string::npos) + { + break; + } + pos++; n++; } // get name & value pos = 0; - size_t pos2, ll; + size_t pos2=0; + size_t ll=0; for (int i = 0; i < n; ++i) { pos2 = txt.find("=", pos); for (; pos2 > pos; --pos2) // There may be a space before "="; { - if (txt.substr(pos2 - 1, 1) != " ") - break; + if (txt.substr(pos2 - 1, 1) != " ") + { + break; + } } ll = pos2 - pos; name[i] = txt.substr(pos, ll); - // std::cout<getnameval(ifs, nparameter, name, val); for (int ip = 0; ip < nparameter; ++ip) @@ -347,12 +364,15 @@ void Pseudopot_upf::read_pseudo_upf201_header(std::ifstream& ifs, Atom_pseudo& p void Pseudopot_upf::read_pseudo_upf201_mesh(std::ifstream& ifs, Atom_pseudo& pp) { - std::string name[50]; - std::string val[50]; - int nparameter; + const int max_n = 100; + std::string name[max_n]; + std::string val[max_n]; + int nparameter=0; + if (ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "getnameval(ifs, nparameter, name, val); + for (int ip = 0; ip < nparameter; ++ip) { if (name[ip] == "dx") @@ -362,6 +382,7 @@ void Pseudopot_upf::read_pseudo_upf201_mesh(std::ifstream& ifs, Atom_pseudo& pp) else if (name[ip] == "mesh") { pp.mesh = atoi(val[ip].c_str()); + this->mesh_changed = false; if (pp.mesh % 2 == 0) { @@ -437,9 +458,12 @@ void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs, Atom_pseudo& return; } std::string word; - std::string name[50]; - std::string val[50]; - int nparameter; + + const int max_n = 100; + std::string name[max_n]; + std::string val[max_n]; + int nparameter=0; + this->kbeta = std::vector(pp.nbeta); pp.lll = std::vector(pp.nbeta); this->els_beta = std::vector(pp.nbeta); @@ -680,9 +704,12 @@ void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs, Atom_pseudo& void Pseudopot_upf::read_pseudo_upf201_pswfc(std::ifstream& ifs, Atom_pseudo& pp) { std::string word; - std::string name[50]; - std::string val[50]; - int nparameter; + + const int max_n = 100; + std::string name[max_n]; + std::string val[max_n]; + int nparameter=0; + ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, ""); pp.els = std::vector(pp.nchi, ""); pp.lchi = std::vector(pp.nchi, 0); @@ -815,9 +842,12 @@ void Pseudopot_upf::read_pseudo_upf201_fullwfc(std::ifstream& ifs) void Pseudopot_upf::read_pseudo_upf201_so(std::ifstream& ifs, Atom_pseudo& pp) { std::string word; - std::string name[50]; - std::string val[50]; - int nparameter; + + const int max_n = 100; + std::string name[max_n]; + std::string val[max_n]; + int nparameter=0; + ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, ""); pp.jchi = std::vector(pp.nchi, 0.0); pp.jjj = std::vector(pp.nbeta, 0.0); @@ -893,4 +923,4 @@ void Pseudopot_upf::read_pseudo_upf201_so(std::ifstream& ifs, Atom_pseudo& pp) ModuleBase::GlobalFunc::SCAN_END(ifs, word); } ModuleBase::GlobalFunc::SCAN_END(ifs, ""); -} \ No newline at end of file +} diff --git a/source/module_cell/read_stru.cpp b/source/module_cell/read_stru.cpp index 589d41ce82b..8cba080427d 100644 --- a/source/module_cell/read_stru.cpp +++ b/source/module_cell/read_stru.cpp @@ -3,92 +3,89 @@ #include "module_base/vector3.h" #include "module_base/mathzone.h" - -namespace unitcell +bool unitcell::check_tau(const Atom* atoms, + const int& ntype, + const double& lat0) { - bool check_tau(const Atom* atoms, - const int& ntype, - const double& lat0) - { - ModuleBase::TITLE("UnitCell","check_tau"); - ModuleBase::timer::tick("UnitCell","check_tau"); - - ModuleBase::Vector3 diff = 0.0; - double norm = 0.0; - double tolerence_bohr = 1.0e-3; + ModuleBase::TITLE("UnitCell","check_tau"); + ModuleBase::timer::tick("UnitCell","check_tau"); + + ModuleBase::Vector3 diff = 0.0; + double norm = 0.0; + double tolerence_bohr = 1.0e-3; + + for(int T1=0; T1< ntype; T1++) + { + for(int I1=0; I1< atoms[T1].na; I1++) + { + double shortest_norm = 10000.0; // a large number + for(int T2=0; T2 norm ) + { + shortest_norm = norm; + } + if( norm < tolerence_bohr ) // unit is Bohr + { + GlobalV::ofs_warning << " two atoms are too close!" << std::endl; + GlobalV::ofs_warning << " type:" << atoms[T1].label << " atom " << I1 + 1 << std::endl; + GlobalV::ofs_warning << " type:" << atoms[T2].label << " atom " << I2 + 1 << std::endl; + GlobalV::ofs_warning << " distance = " << norm << " Bohr" << std::endl; + return false; + } + } + } + } + } + } + ModuleBase::timer::tick("UnitCell","check_tau"); + return true; +} - for(int T1=0; T1< ntype; T1++) - { - for(int I1=0; I1< atoms[T1].na; I1++) - { - double shortest_norm = 10000.0; // a large number - for(int T2=0; T2 norm ) - { - shortest_norm = norm; - } - if( norm < tolerence_bohr ) // unit is Bohr - { - GlobalV::ofs_warning << " two atoms are too close!" << std::endl; - GlobalV::ofs_warning << " type:" << atoms[T1].label << " atom " << I1 + 1 << std::endl; - GlobalV::ofs_warning << " type:" << atoms[T2].label << " atom " << I2 + 1 << std::endl; - GlobalV::ofs_warning << " distance = " << norm << " Bohr" << std::endl; - return false; - } - } - } - } - } - } - ModuleBase::timer::tick("UnitCell","check_tau"); - return true; - } - void check_dtau(Atom* atoms, - const int& ntype, - const double& lat0, - ModuleBase::Matrix3& latvec) - { - for(int it=0; ittaud[ia].x=fmod(atom1->taud[ia].x + 10000,1.0); - atom1->taud[ia].y=fmod(atom1->taud[ia].y + 10000,1.0); - atom1->taud[ia].z=fmod(atom1->taud[ia].z + 10000,1.0); +void unitcell::check_dtau(Atom* atoms, + const int& ntype, + const double& lat0, + ModuleBase::Matrix3& latvec) +{ + for(int it=0; ittaud[ia].x=fmod(atom1->taud[ia].x + 10000,1.0); + atom1->taud[ia].y=fmod(atom1->taud[ia].y + 10000,1.0); + atom1->taud[ia].z=fmod(atom1->taud[ia].z + 10000,1.0); + + double cx2=0.0; + double cy2=0.0; + double cz2=0.0; - double cx2=0.0; - double cy2=0.0; - double cz2=0.0; + ModuleBase::Mathzone::Direct_to_Cartesian( + atom1->taud[ia].x, atom1->taud[ia].y, atom1->taud[ia].z, + latvec.e11, latvec.e12, latvec.e13, + latvec.e21, latvec.e22, latvec.e23, + latvec.e31, latvec.e32, latvec.e33, + cx2, cy2, cz2); - ModuleBase::Mathzone::Direct_to_Cartesian( - atom1->taud[ia].x, atom1->taud[ia].y, atom1->taud[ia].z, - latvec.e11, latvec.e12, latvec.e13, - latvec.e21, latvec.e22, latvec.e23, - latvec.e31, latvec.e32, latvec.e33, - cx2, cy2, cz2); + atom1->tau[ia].x = cx2; + atom1->tau[ia].y = cy2; + atom1->tau[ia].z = cz2; - atom1->tau[ia].x = cx2; - atom1->tau[ia].y = cy2; - atom1->tau[ia].z = cz2; - - } - } - return; - } -} \ No newline at end of file + } + } + return; +} diff --git a/source/module_cell/test/support/mock_unitcell.cpp b/source/module_cell/test/support/mock_unitcell.cpp index d2e7268f880..e62a266d6d1 100644 --- a/source/module_cell/test/support/mock_unitcell.cpp +++ b/source/module_cell/test/support/mock_unitcell.cpp @@ -17,15 +17,23 @@ UnitCell::~UnitCell() { delete[] atoms; } } + void UnitCell::print_cell(std::ofstream& ofs) const {} + void UnitCell::set_iat2itia() {} + void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) {} + bool UnitCell::if_atoms_can_move() const { return true; } + bool UnitCell::if_cell_can_change() const { return true; } + void UnitCell::setup(const std::string& latname_in, const int& ntype_in, const int& lmaxmax_in, const bool& init_vel_in, const std::string& fixed_axes_in) {} + void cal_nelec(const Atom* atoms, const int& ntype, double& nelec) {} -void UnitCell::compare_atom_labels(std::string label1, std::string label2) {} \ No newline at end of file + +void UnitCell::compare_atom_labels(const std::string &label1, const std::string &label2) {} diff --git a/source/module_cell/test/unitcell_test.cpp b/source/module_cell/test/unitcell_test.cpp index 70d058d090b..397d11e5786 100644 --- a/source/module_cell/test/unitcell_test.cpp +++ b/source/module_cell/test/unitcell_test.cpp @@ -1171,7 +1171,7 @@ TEST_F(UcellTestReadStru, ReadAtomSpecies) PARAM.input.basis_type = "lcao"; PARAM.sys.deepks_setorb = true; EXPECT_NO_THROW(unitcell::read_atom_species(ifa, ofs_running,*ucell)); - EXPECT_NO_THROW(unitcell::read_lattice_constant(ifa, ofs_running,ucell->lat)); + EXPECT_NO_THROW(unitcell::read_lattice_constant(ifa, ofs_running, ucell->lat)); EXPECT_DOUBLE_EQ(ucell->latvec.e11, 4.27957); EXPECT_DOUBLE_EQ(ucell->latvec.e22, 4.27957); EXPECT_DOUBLE_EQ(ucell->latvec.e33, 4.27957); @@ -1185,7 +1185,7 @@ TEST_F(UcellTestReadStru, ReadAtomSpeciesWarning1) std::string fn = "./support/STRU_MgO_Warning1"; std::ifstream ifa(fn.c_str()); std::ofstream ofs_running; - ofs_running.open("read_atom_species.tmp"); + ofs_running.open("read_atom_species.txt"); ucell->ntype = 2; ucell->atoms = new Atom[ucell->ntype]; ucell->set_atom_flag = true; @@ -1195,7 +1195,7 @@ TEST_F(UcellTestReadStru, ReadAtomSpeciesWarning1) EXPECT_THAT(output, testing::HasSubstr("unrecongnized pseudo type.")); ofs_running.close(); ifa.close(); - remove("read_atom_species.tmp"); + //remove("read_atom_species.txt"); } TEST_F(UcellTestReadStru, ReadLatticeConstantWarning1) diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index 5f4dd37c866..088f89799c8 100755 --- a/source/module_cell/unitcell.cpp +++ b/source/module_cell/unitcell.cpp @@ -180,8 +180,10 @@ std::vector> UnitCell::get_constrain() const //============================================================== // Calculate various lattice related quantities for given latvec //============================================================== -void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { +void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) +{ ModuleBase::TITLE("UnitCell", "setup_cell"); + // (1) init mag assert(ntype > 0); delete[] magnet.start_magnetization; @@ -203,7 +205,9 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { this->pseudo_fn.resize(ntype); this->pseudo_type.resize(ntype); this->orbital_fn.resize(ntype); - if (GlobalV::MY_RANK == 0) { + + if (GlobalV::MY_RANK == 0) + { // open "atom_unitcell" file. std::ifstream ifa(fn.c_str(), std::ios::in); if (!ifa) @@ -268,7 +272,7 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { //========================== // call read_atom_positions //========================== - ok2 = unitcell::read_atom_positions(*this,ifa, log, GlobalV::ofs_warning); + ok2 = unitcell::read_atom_positions(*this, ifa, log, GlobalV::ofs_warning); } } #ifdef __MPI @@ -295,15 +299,23 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { // Firstly, latvec must be read in. //======================================================== assert(lat0 > 0.0); - this->omega = latvec.Det() * this->lat0 * lat0 * lat0; + this->omega = latvec.Det() * this->lat0 * this->lat0 * this->lat0; + + std::cout << "latvec.det=" << latvec.Det() << std::endl; + std::cout << "lat0=" << lat0 << std::endl; + + if (this->omega < 0) { std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; std::cout << " Warning: The lattice vector is left-handed; a right-handed vector is prefered." << std::endl; std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - GlobalV::ofs_warning << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; - GlobalV::ofs_warning << " Warning: The lattice vector is left-handed; a right-handed vector is prefered." << std::endl; - GlobalV::ofs_warning << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + GlobalV::ofs_warning << + "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; + GlobalV::ofs_warning << + " Warning: The lattice vector is left-handed; a right-handed vector is prefered." << std::endl; + GlobalV::ofs_warning << + "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << std::endl; this->omega = std::abs(this->omega); } else if (this->omega == 0) @@ -340,9 +352,6 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { log, "Reciprocal vectors: (Cartesian coordinate: in unit of 2 pi/a_0)", G); - // OUT(log,"lattice center x",latcenter.x); - // OUT(log,"lattice center y",latcenter.y); - // OUT(log,"lattice center z",latcenter.z); //=================================== // set index for iat2it, iat2ia @@ -350,18 +359,21 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { this->set_iat2itia(); #ifdef USE_PAW - if (PARAM.inp.use_paw) { + if (PARAM.inp.use_paw) + { GlobalC::paw_cell.set_libpaw_cell(latvec, lat0); - int* typat; - double* xred; + int* typat = nullptr; + double* xred = nullptr; typat = new int[nat]; xred = new double[nat * 3]; int iat = 0; - for (int it = 0; it < ntype; it++) { - for (int ia = 0; ia < atoms[it].na; ia++) { + for (int it = 0; it < ntype; it++) + { + for (int ia = 0; ia < atoms[it].na; ia++) + { typat[iat] = it + 1; // Fortran index starts from 1 !!!! xred[iat * 3 + 0] = atoms[it].taud[ia].x; xred[iat * 3 + 1] = atoms[it].taud[ia].y; @@ -384,7 +396,8 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) { } -void UnitCell::set_iat2iwt(const int& npol_in) { +void UnitCell::set_iat2iwt(const int& npol_in) +{ #ifdef __DEBUG assert(npol_in == 1 || npol_in == 2); assert(this->nat > 0); @@ -394,8 +407,11 @@ void UnitCell::set_iat2iwt(const int& npol_in) { this->npol = npol_in; int iat = 0; int iwt = 0; - for (int it = 0; it < this->ntype; it++) { - for (int ia = 0; ia < atoms[it].na; ia++) { + + for (int it = 0; it < this->ntype; it++) + { + for (int ia = 0; ia < atoms[it].na; ia++) + { this->iat2iwt[iat] = iwt; iwt += atoms[it].nw * this->npol; ++iat; @@ -407,25 +423,31 @@ void UnitCell::set_iat2iwt(const int& npol_in) { // check if any atom can be moved -bool UnitCell::if_atoms_can_move() const { - for (int it = 0; it < this->ntype; it++) { +bool UnitCell::if_atoms_can_move() const +{ + for (int it = 0; it < this->ntype; it++) + { Atom* atom = &atoms[it]; - for (int ia = 0; ia < atom->na; ia++) { - if (atom->mbl[ia].x || atom->mbl[ia].y || atom->mbl[ia].z) { - return true; -} - } - } + for (int ia = 0; ia < atom->na; ia++) + { + if (atom->mbl[ia].x || atom->mbl[ia].y || atom->mbl[ia].z) + { + return true; + } + } + } return false; } // check if lattice vector can be changed -bool UnitCell::if_cell_can_change() const { - // need to be fixed next - if (this->lc[0] || this->lc[1] || this->lc[2]) { - return true; - } - return false; +bool UnitCell::if_cell_can_change() const +{ + // need to be fixed next + if (this->lc[0] || this->lc[1] || this->lc[2]) + { + return true; + } + return false; } void UnitCell::setup(const std::string& latname_in, @@ -498,7 +520,8 @@ void UnitCell::setup(const std::string& latname_in, } -void UnitCell::compare_atom_labels(std::string label1, std::string label2) { +void UnitCell::compare_atom_labels(const std::string &label1, const std::string &label2) +{ if (label1!= label2) //'!( "Ag" == "Ag" || "47" == "47" || "Silver" == Silver" )' { atom_in ai; @@ -515,22 +538,30 @@ void UnitCell::compare_atom_labels(std::string label1, std::string label2) { { std::string stru_label = ""; std::string psuedo_label = ""; - for (int ip = 0; ip < label1.length(); ip++) { - if (!(isdigit(label1[ip]) || label1[ip] == '_')) { - stru_label += label1[ip]; - } else { - break; - } - } - stru_label[0] = toupper(stru_label[0]); - - for (int ip = 0; ip < label2.length(); ip++) { - if (!(isdigit(label2[ip]) || label2[ip] == '_')) { - psuedo_label += label2[ip]; - } else { - break; - } - } + for (int ip = 0; ip < label1.length(); ip++) + { + if (!(isdigit(label1[ip]) || label1[ip] == '_')) + { + stru_label += label1[ip]; + } + else + { + break; + } + } + stru_label[0] = toupper(stru_label[0]); + + for (int ip = 0; ip < label2.length(); ip++) + { + if (!(isdigit(label2[ip]) || label2[ip] == '_')) + { + psuedo_label += label2[ip]; + } + else + { + break; + } + } psuedo_label[0] = toupper(psuedo_label[0]); if (!(stru_label == psuedo_label diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h index ea9e4576f14..e7b3e3bdcfd 100644 --- a/source/module_cell/unitcell.h +++ b/source/module_cell/unitcell.h @@ -233,7 +233,7 @@ class UnitCell { /// @brief check consistency between two atom labels from STRU and pseudo or /// orb file - void compare_atom_labels(std::string label1, std::string label2); + void compare_atom_labels(const std::string &label1, const std::string &label2); /// @brief get atomCounts, which is a map from element type to atom number std::map get_atom_Counts() const; /// @brief get orbitalCounts, which is a map from element type to orbital diff --git a/source/module_hamilt_general/operator.cpp b/source/module_hamilt_general/operator.cpp index 08a5ba97ccd..6682777322a 100644 --- a/source/module_hamilt_general/operator.cpp +++ b/source/module_hamilt_general/operator.cpp @@ -16,8 +16,10 @@ Operator::~Operator() { delete this->hpsi; } + Operator* last = this->next_op; Operator* last_sub = this->next_sub_op; + while (last != nullptr || last_sub != nullptr) { if (last_sub != nullptr) @@ -58,7 +60,6 @@ typename Operator::hpsi_info Operator::hPsi(hpsi_info& inp T* hpsi_pointer = std::get<2>(input); if (this->in_place) { - // ModuleBase::GlobalFunc::COPYARRAY(this->hpsi->get_pointer(), hpsi_pointer, this->hpsi->size()); syncmem_op()(hpsi_pointer, this->hpsi->get_pointer(), this->hpsi->size()); delete this->hpsi; this->hpsi = new psi::Psi(hpsi_pointer, @@ -78,31 +79,35 @@ typename Operator::hpsi_info Operator::hPsi(hpsi_info& inp psi_input->get_nbasis(), true); - switch (op->get_act_type()) - { - case 2: - op->act(psi_wrapper, *this->hpsi, nbands); - break; - default: - op->act(nbands, - psi_input->get_nbasis(), - psi_input->get_npol(), - tmpsi_in, - this->hpsi->get_pointer(), - psi_input->get_current_nbas(), - is_first_node); - break; - } + switch (op->get_act_type()) + { + case 2: + op->act(psi_wrapper, *this->hpsi, nbands); + break; + default: + op->act(nbands, + psi_input->get_nbasis(), + psi_input->get_npol(), + tmpsi_in, + this->hpsi->get_pointer(), + psi_input->get_current_nbas(), + is_first_node); + break; + } }; ModuleBase::timer::tick("Operator", "hPsi"); + call_act(this, true); // first node + Operator* node((Operator*)this->next_op); + while (node != nullptr) { call_act(node, false); // other nodes node = (Operator*)(node->next_op); } + ModuleBase::timer::tick("Operator", "hPsi"); return hpsi_info(this->hpsi, psi::Range(1, 0, 0, nbands / psi_input->get_npol()), hpsi_pointer); @@ -157,6 +162,7 @@ void Operator::add(Operator* next) } } + template T* Operator::get_hpsi(const hpsi_info& info) const { diff --git a/source/module_hamilt_general/operator.h b/source/module_hamilt_general/operator.h index 941f3e29c1a..87567f9f362 100644 --- a/source/module_hamilt_general/operator.h +++ b/source/module_hamilt_general/operator.h @@ -73,8 +73,10 @@ class Operator /// developer-friendly interfaces for act() function /// interface type 2: input and change the Psi-type HPsi - // virtual void act(const psi::Psi& psi_in, psi::Psi& psi_out) const {}; - virtual void act(const psi::Psi& psi_in, psi::Psi& psi_out, const int nbands) const {}; + // virtual void act(const psi::Psi& psi_in, psi::Psi& psi_out) const {}; + virtual void act(const psi::Psi& psi_in, + psi::Psi& psi_out, + const int nbands) const {}; /// interface type 3: return a Psi-type HPsi // virtual psi::Psi act(const psi::Psi& psi_in) const { return psi_in; }; @@ -125,4 +127,4 @@ class Operator } // end namespace hamilt -#endif \ No newline at end of file +#endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 35184682e07..8831d8a2749 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -55,10 +55,8 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, this->kv = &kv_in; - // Real space Hamiltonian is inited with template TR - // this->hR = new HContainer(paraV); + // initialize the overlap matrix this->sR = new HContainer(paraV); - // this->hsk = new HS_Matrix_K(paraV); this->getOperator() = new OverlapNew>(this->hsk, this->kv->kvec_d, @@ -94,9 +92,10 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, #endif ) { - this->kv = &kv_in; this->classname = "HamiltLCAO"; + this->kv = &kv_in; + // Real space Hamiltonian is inited with template TR this->hR = new HContainer(paraV); this->sR = new HContainer(paraV); @@ -227,7 +226,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, if (PARAM.inp.dft_plus_u == 2) { dftu = new OperatorDFTU>(this->hsk, - kv->kvec_d, + this->kv->kvec_d, this->hR, // no explicit call yet this->kv->isk); } @@ -261,7 +260,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, // Veff term this->getOperator() = new Veff>(GK_in, this->hsk, - kv->kvec_d, + this->kv->kvec_d, pot_in, this->hR, &ucell, @@ -353,7 +352,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, } Operator* td_ekinetic = new TDEkinetic>(this->hsk, this->hR, - kv, + this->kv, &ucell, orb.cutoffs(), &grid_d, @@ -370,7 +369,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, if (PARAM.inp.dft_plus_u == 2) { dftu = new OperatorDFTU>(this->hsk, - kv->kvec_d, + this->kv->kvec_d, this->hR, // no explicit call yet this->kv->isk); } @@ -390,7 +389,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, if (PARAM.inp.sc_mag_switch) { Operator* sc_lambda = new DeltaSpin>(this->hsk, - kv->kvec_d, + this->kv->kvec_d, this->hR, ucell, &grid_d, @@ -411,7 +410,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, Operator* exx = new OperatorEXX>(this->hsk, this->hR, ucell, - *this->kv, + *kv, Hexxd, Hexxc, Add_Hexx_Type::R, @@ -422,6 +421,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->getOperator()->add(exx); } #endif + // if NSPIN==2, HR should be separated into two parts, save HR into this->hRS2 int memory_fold = 1; if (PARAM.inp.nspin == 2) @@ -451,6 +451,7 @@ void HamiltLCAO::updateHk(const int ik) { ModuleBase::TITLE("HamiltLCAO", "updateHk"); ModuleBase::timer::tick("HamiltLCAO", "updateHk"); + // update global spin index if (PARAM.inp.nspin == 2) { @@ -497,21 +498,30 @@ Operator*& HamiltLCAO::getOperator() } template -void HamiltLCAO::updateSk(const int ik, const int hk_type) +void HamiltLCAO::updateSk( + const int ik, + const int hk_type) { ModuleBase::TITLE("HamiltLCAO", "updateSk"); ModuleBase::timer::tick("HamiltLCAO", "updateSk"); + ModuleBase::GlobalFunc::ZEROS(this->getSk(), this->get_size_hsk()); + if (hk_type == 1) // collumn-major matrix for SK { const int nrow = this->hsk->get_pv()->get_row_size(); - hamilt::folding_HR(*this->sR, this->getSk(), this->kv->kvec_d[ik], nrow, 1); - } - else if (hk_type == 0) // row-major matrix for SK - { + hamilt::folding_HR(*this->sR, this->getSk(), this->kv->kvec_d[ik], nrow, 1); + } + else if (hk_type == 0) // row-major matrix for SK + { const int ncol = this->hsk->get_pv()->get_col_size(); hamilt::folding_HR(*this->sR, this->getSk(), this->kv->kvec_d[ik], ncol, 0); } + else + { + ModuleBase::WARNING_QUIT("updateSk","the value of hk_type is incorrect."); + } + ModuleBase::timer::tick("HamiltLCAO", "updateSk"); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h index c7ead226253..554e23cf301 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h @@ -1,5 +1,5 @@ -#ifndef W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_HAMILT_LCAO_HAMILT_LCAODFT_HAMILT_LCAO_H -#define W_ABACUS_DEVELOP_ABACUS_DEVELOP_SOURCE_MODULE_HAMILT_LCAO_HAMILT_LCAODFT_HAMILT_LCAO_H +#ifndef HAMILT_LCAO_H +#define HAMILT_LCAO_H #include "module_basis/module_nao/two_center_bundle.h" #include "module_cell/klist.h" @@ -25,24 +25,30 @@ namespace hamilt { // template first for type of k space H matrix elements -// template second for type of temporary matrix, gamma_only fix-gamma-matrix + S-gamma, multi-k fix-Real + S-Real +// template second for type of temporary matrix, +// gamma_only fix-gamma-matrix + S-gamma, +// multi-k fix-Real + S-Real template class HamiltLCAO : public Hamilt { public: + + + using TAC = std::pair>; + + /** * @brief Constructor of Hamiltonian for LCAO base * HR and SR will be allocated with Operators */ - using TAC = std::pair>; HamiltLCAO(Gint_Gamma* GG_in, Gint_k* GK_in, const UnitCell& ucell, const Grid_Driver& grid_d, - const Parallel_Orbitals* paraV, - elecstate::Potential* pot_in, - const K_Vectors& kv_in, - const TwoCenterBundle& two_center_bundle, + const Parallel_Orbitals* paraV, + elecstate::Potential* pot_in, + const K_Vectors& kv_in, + const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, elecstate::DensityMatrix* DM_in #ifdef __DEEPKS @@ -57,6 +63,7 @@ class HamiltLCAO : public Hamilt std::vector>>>>* Hexxc = nullptr #endif ); + /** * @brief Constructor of vacuum Operators, only HR and SR will be initialed as empty HContainer */ @@ -80,30 +87,36 @@ class HamiltLCAO : public Hamilt /// get pointer of Operator ops Operator*& getOperator(); - /// get hk-pointer + + /// get H(k) pointer TK* getHk() const { return this->hsk->get_hk(); } - /// get sk-pointer + + /// get S(k) pointer TK* getSk() const { return this->hsk->get_sk(); } + int get_size_hsk() const { return this->hsk->get_size(); } + /// get HR pointer of *this->hR, which is a HContainer and contains H(R) HContainer*& getHR() { return this->hR; } + /// get SR pointer of *this->sR, which is a HContainer and contains S(R) HContainer*& getSR() { return this->sR; } + #ifdef __DEEPKS /// get V_delta_R pointer of *this->V_delta_R, which is a HContainer and contains V_delta(R) HContainer*& get_V_delta_R() @@ -111,6 +124,7 @@ class HamiltLCAO : public Hamilt return this->V_delta_R; } #endif + /// refresh the status of HR void refresh() override; @@ -121,41 +135,45 @@ class HamiltLCAO : public Hamilt * @brief special for LCAO, update SK only * * @param ik target K point + * @param kvec_d: direct coordinates of k-points * @param hk_type 0: SK is row-major, 1: SK is collumn-major * @return void */ - void updateSk(const int ik, const int hk_type = 0); + void updateSk(const int ik, const int hk_type = 0); // core function: return H(k) and S(k) matrixs for direct solving eigenvalues. // not used in PW base - // void matrix(MatrixBlock> &hk_in, MatrixBlock> &sk_in) override; void matrix(MatrixBlock& hk_in, MatrixBlock& sk_in) override; private: + const K_Vectors* kv = nullptr; - // Real space Hamiltonian + //! Real space Hamiltonian H(R), where R is the Bravis lattice vector HContainer* hR = nullptr; + + //! Real space overlap matrix S(R), where R is the Bravis lattice vector HContainer* sR = nullptr; #ifdef __DEEPKS HContainer* V_delta_R = nullptr; #endif + //! Hamiltonian and overlap matrices for a specific k point HS_Matrix_K* hsk = nullptr; + // special case for NSPIN=2 , data of HR should be separated into two parts // save them in this->hRS2; std::vector hRS2; + int refresh_times = 1; - /// current_spin for NSPIN=2, 0: hamiltonian for spin up, 1: hamiltonian for spin down + //! current_spin for NSPIN=2 case + //! 0: Hamiltonian for spin up, + //! 1: Hamiltonian for spin down int current_spin = 0; const int istep = 0; - - // sk and hk will be refactored to HamiltLCAO later - // std::vector sk; - // std::vector hk; }; } // namespace hamilt diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp index 495af32841e..429a4ccb9c3 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp @@ -1,5 +1,5 @@ -#ifndef ABACUS_HS_MATRIX_K_HPP -#define ABACUS_HS_MATRIX_K_HPP +#ifndef HS_MATRIX_K_HPP +#define HS_MATRIX_K_HPP #include "module_basis/module_ao/parallel_orbitals.h" @@ -9,6 +9,7 @@ namespace hamilt template class HS_Matrix_K { + public: HS_Matrix_K(const Parallel_Orbitals* paraV, bool no_s=false){ hk.resize(paraV->nloc); @@ -18,17 +19,27 @@ namespace hamilt } this->pv = paraV; } + TK* get_hk() {return hk.data();} + TK* get_sk() {return sk.data();} + int get_size() {return hk.size();} + void set_zero_hk() {hk.assign(hk.size(), 0);} + void set_zero_sk() {sk.assign(sk.size(), 0);} + const Parallel_Orbitals* get_pv() const {return this->pv;} + private: + std::vector hk; + std::vector sk; + const Parallel_Orbitals* pv = nullptr; }; } -#endif \ No newline at end of file +#endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h index c29ba87c0ed..f9ad63c9422 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h @@ -100,6 +100,8 @@ class OperatorLCAO : public Operator { // protected: // Hamiltonian matrices which are calculated in OperatorLCAO HS_Matrix_K* hsk = nullptr; + + // kvec_d const std::vector>& kvec_d; protected: diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp index 382f53350e6..f48af2e65f0 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp @@ -31,6 +31,7 @@ TDEkinetic>::TDEkinetic(HS_Matrix_K* hsk_in, // initialize HR to get adjs info. this->initialize_HR(Grid); } + template TDEkinetic>::~TDEkinetic() { @@ -40,12 +41,14 @@ TDEkinetic>::~TDEkinetic() } TD_Velocity::td_vel_op = nullptr; } + // term A^2*S template void TDEkinetic>::td_ekinetic_scalar(std::complex* Hloc,const TR& overlap, int nnr) { return; } + // term A^2*S template <> void TDEkinetic, double>>::td_ekinetic_scalar(std::complex* Hloc, @@ -57,6 +60,7 @@ void TDEkinetic, double>>::td_ekinetic_scalar( Hloc[nnr] += tmp; return; } + // term A dot βˆ‡ template void TDEkinetic>::td_ekinetic_grad(std::complex* Hloc, @@ -384,8 +388,11 @@ void TDEkinetic, double>>::contributeHk(int ik const Parallel_Orbitals* paraV = this->hR_tmp->get_atom_pair(0).get_paraV(); // save HR data for output int spin_tot = PARAM.inp.nspin; + if (spin_tot == 4) - ; + { + + } else if (!output_hR_done && TD_Velocity::out_mat_R) { for (int spin_now = 0; spin_now < spin_tot; spin_now++) @@ -398,6 +405,7 @@ void TDEkinetic, double>>::contributeHk(int ik } output_hR_done = true; } + // folding inside HR to HK if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver)) { diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h index c5ec078c14f..7f5881ac1fc 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h @@ -45,7 +45,9 @@ class TDEkinetic> : public OperatorLCAO ~TDEkinetic(); virtual void contributeHR() override; + virtual void contributeHk(int ik) override; + /// @brief init two center integrals and vector potential for td_ekintic term void init_td(); @@ -78,22 +80,30 @@ class TDEkinetic> : public OperatorLCAO * loop the atom-pairs in HR and calculate the ekinetic matrix */ void calculate_HR(); + virtual void set_HR_fixed(void*) override; TD_Velocity td_velocity; private: + const UnitCell* ucell = nullptr; + std::vector orb_cutoff_; HContainer* SR = nullptr; - /// @brief Store real space hamiltonian. TD term should include imaginary part, thus it has to be complex type. Only shared between TD operators. + + /// @brief Store real space hamiltonian. TD term should include imaginary part, + /// thus it has to be complex type. Only shared between TD operators. HContainer>* hR_tmp = nullptr; + const Grid_Driver* Grid = nullptr; const K_Vectors* kv; - /// @brief correction term iAβ‹…βˆ‡ + + /// @brief correction term i A nabla void td_ekinetic_scalar(std::complex* Hloc, const TR& Sloc, int nnr); + /// @brief correction term A^2*S void td_ekinetic_grad(std::complex* Hloc, int nnr, ModuleBase::Vector3 grad_overlap); @@ -106,7 +116,9 @@ class TDEkinetic> : public OperatorLCAO std::vector adjs_all; bool hR_tmp_done = false; + bool allocated = false; + bool output_hR_done = false; }; diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp index ffdc80513c9..1536a232d91 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp @@ -334,8 +334,11 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, for (int ik = 0; ik < nks; ik++) { h_tot[ik].create(nlocal, nlocal); + p_ham->updateHk(ik); + const TK* hk_ptr = p_ham->getHk(); + for (int i = 0; i < ParaV->nloc; i++) { h_mat[ik][i] = hk_ptr[i]; diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.h b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.h index 2f32c664782..075650b3f02 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.h +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.h @@ -51,4 +51,4 @@ class LCAO_Deepks_Interface }; #endif -#endif \ No newline at end of file +#endif diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h index a885340cf15..2ec3cf1ed62 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.h +++ b/source/module_hamilt_lcao/module_dftu/dftu.h @@ -1,7 +1,3 @@ -//========================================================== -// Author: Xin Qu -// DATE : 2019-12-10 -//========================================================== #ifndef DFTU_H #define DFTU_H @@ -187,14 +183,15 @@ class DFTU const int dim1, const int dim2, std::complex* mat_k, - const std::vector>& kvec_d); + const ModuleBase::Vector3& kvec_d); /** * @brief new function of folding_S_matrix * only for Hamiltonian now, for force and stress will be developed later * use HContainer as input and output in mat_k */ - void folding_matrix_k_new(const int ik, hamilt::Hamilt>* p_ham); + void folding_matrix_k_new(const int ik, + hamilt::Hamilt>* p_ham); //============================================================= // In dftu_force.cpp @@ -218,7 +215,7 @@ class DFTU const int ik, const std::complex* rho_VU, ModuleBase::matrix& force_dftu, - const std::vector>& kvec_d); + const ModuleBase::Vector3& kvec_d); void cal_stress_k(const UnitCell& ucell, const Grid_Driver& gd, @@ -227,7 +224,7 @@ class DFTU const int ik, const std::complex* rho_VU, ModuleBase::matrix& stress_dftu, - const std::vector>& kvec_d); + const ModuleBase::Vector3& kvec_d); void cal_force_gamma(const UnitCell& ucell, const double* rho_VU, diff --git a/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp b/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp index 3ae80f8e141..ee3b6413848 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_folding.cpp @@ -133,7 +133,7 @@ void DFTU::folding_matrix_k(const UnitCell& ucell, const int dim1, const int dim2, std::complex* mat_k, - const std::vector>& kvec_d) + const ModuleBase::Vector3& kvec_d) { ModuleBase::TITLE("DFTU", "folding_matrix_k"); ModuleBase::timer::tick("DFTU", "folding_matrix_k"); @@ -224,7 +224,7 @@ void DFTU::folding_matrix_k(const UnitCell& ucell, // dR is the index of box in Crystal coordinates //------------------------------------------------ ModuleBase::Vector3 dR(gd.getBox(ad).x, gd.getBox(ad).y, gd.getBox(ad).z); - const double arg = (kvec_d[ik] * dR) * ModuleBase::TWO_PI; + const double arg = (kvec_d * dR) * ModuleBase::TWO_PI; const std::complex kphase = std::complex(cos(arg), sin(arg)); //-------------------------------------------------- @@ -236,15 +236,19 @@ void DFTU::folding_matrix_k(const UnitCell& ucell, // the index of orbitals in this processor const int iw1_all = start1 + ii; const int mu = pv.global2local_row(iw1_all); - if (mu < 0) { continue; -} + if (mu < 0) + { + continue; + } for (int jj = 0; jj < atom2->nw * PARAM.globalv.npol; jj++) { int iw2_all = start2 + jj; const int nu = pv.global2local_col(iw2_all); - if (nu < 0) { continue; -} + if (nu < 0) + { + continue; + } int iic; if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver)) @@ -293,22 +297,23 @@ void DFTU::folding_matrix_k_new(const int ik, // get SR and fold to mat_k if(PARAM.globalv.gamma_only_local) { - dynamic_cast*>(p_ham)->updateSk(ik, hk_type); + dynamic_cast*>(p_ham) + ->updateSk(ik, hk_type); } else { if(PARAM.inp.nspin != 4) { - dynamic_cast, double>*>(p_ham)->updateSk(ik, hk_type); + dynamic_cast, double>*>(p_ham) + ->updateSk(ik, hk_type); } else { - dynamic_cast, std::complex>*>(p_ham)->updateSk(ik, hk_type); + dynamic_cast, std::complex>*>(p_ham) + ->updateSk(ik, hk_type); } } } - - } // namespace ModuleDFTU -#endif // __LCAO \ No newline at end of file +#endif // __LCAO diff --git a/source/module_hamilt_lcao/module_dftu/dftu_force.cpp b/source/module_hamilt_lcao/module_dftu/dftu_force.cpp index 8d06d9a92f9..f4423f623d6 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_force.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_force.cpp @@ -209,11 +209,11 @@ void DFTU::force_stress(const UnitCell& ucell, if (PARAM.inp.cal_force) { - cal_force_k(ucell, gd, fsr, pv, ik, &rho_VU[0], force_dftu, kv.kvec_d); + cal_force_k(ucell, gd, fsr, pv, ik, &rho_VU[0], force_dftu, kv.kvec_d[ik]); } if (PARAM.inp.cal_stress) { - cal_stress_k(ucell, gd, fsr, pv, ik, &rho_VU[0], stress_dftu, kv.kvec_d); + cal_stress_k(ucell, gd, fsr, pv, ik, &rho_VU[0], stress_dftu, kv.kvec_d[ik]); } } // ik } @@ -256,7 +256,7 @@ void DFTU::cal_force_k(const UnitCell& ucell, const int ik, const std::complex* rho_VU, ModuleBase::matrix& force_dftu, - const std::vector>& kvec_d) + const ModuleBase::Vector3& kvec_d) { ModuleBase::TITLE("DFTU", "cal_force_k"); ModuleBase::timer::tick("DFTU", "cal_force_k"); @@ -386,7 +386,7 @@ void DFTU::cal_stress_k(const UnitCell& ucell, const int ik, const std::complex* rho_VU, ModuleBase::matrix& stress_dftu, - const std::vector>& kvec_d) + const ModuleBase::Vector3& kvec_d) { ModuleBase::TITLE("DFTU", "cal_stress_k"); ModuleBase::timer::tick("DFTU", "cal_stress_k"); @@ -662,4 +662,4 @@ void DFTU::cal_stress_gamma(const UnitCell& ucell, return; } } // namespace ModuleDFTU -#endif \ No newline at end of file +#endif diff --git a/source/module_hamilt_lcao/module_dftu/dftu_occup.cpp b/source/module_hamilt_lcao/module_dftu/dftu_occup.cpp index d4efeed4269..a745010aa7a 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_occup.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_occup.cpp @@ -38,9 +38,10 @@ void DFTU::copy_locale(const UnitCell& ucell) for (int T = 0; T < ucell.ntype; T++) { - if (orbital_corr[T] == -1) { - continue; -} + if (orbital_corr[T] == -1) + { + continue; + } for (int I = 0; I < ucell.atoms[T].na; I++) { @@ -75,8 +76,10 @@ void DFTU::zero_locale(const UnitCell& ucell) for (int T = 0; T < ucell.ntype; T++) { - if (orbital_corr[T] == -1) { continue; -} + if (orbital_corr[T] == -1) + { + continue; + } for (int I = 0; I < ucell.atoms[T].na; I++) { @@ -114,9 +117,10 @@ void DFTU::mix_locale(const UnitCell& ucell, for (int T = 0; T < ucell.ntype; T++) { - if (orbital_corr[T] == -1) { - continue; -} + if (orbital_corr[T] == -1) + { + continue; + } for (int I = 0; I < ucell.atoms[T].na; I++) { @@ -161,7 +165,8 @@ void DFTU::cal_occup_m_k(const int iter, //=================Part 1====================== // call SCALAPACK routine to calculate the product of the S and density matrix - const char transN = 'N', transT = 'T'; + const char transN = 'N'; + const char transT = 'T'; const int one_int = 1; const std::complex beta(0.0,0.0), alpha(1.0,0.0); @@ -171,7 +176,9 @@ void DFTU::cal_occup_m_k(const int iter, { // srho(mu,nu) = \sum_{iw} S(mu,iw)*dm_k(iw,nu) this->folding_matrix_k_new(ik, p_ham); + std::complex* s_k_pointer = nullptr; + if(PARAM.inp.nspin != 4) { s_k_pointer = dynamic_cast, double>*>(p_ham)->getSk(); @@ -210,9 +217,10 @@ void DFTU::cal_occup_m_k(const int iter, const int NL = ucell.atoms[it].nwl + 1; const int LC = orbital_corr[it]; - if (LC == -1) { - continue; -} + if (LC == -1) + { + continue; + } for (int ia = 0; ia < ucell.atoms[it].na; ia++) { @@ -220,18 +228,20 @@ void DFTU::cal_occup_m_k(const int iter, for (int l = 0; l < NL; l++) { - if (l != orbital_corr[it]) { - continue; -} + if (l != orbital_corr[it]) + { + continue; + } const int N = ucell.atoms[it].l_nchi[l]; for (int n = 0; n < N; n++) { // if(!Yukawa && n!=0) continue; - if (n != 0) { - continue; -} + if (n != 0) + { + continue; + } // Calculate the local occupation number matrix for (int m0 = 0; m0 < 2 * l + 1; m0++) @@ -256,12 +266,16 @@ void DFTU::cal_occup_m_k(const int iter, const int m0_all = m0 + ipol0 * (2 * l + 1); const int m1_all = m1 + ipol1 * (2 * l + 1); - if ((nu >= 0) && (mu >= 0)) - locale[iat][l][n][spin](m0_all, m1_all) += (srho[irc]).real() / 4.0; + if ((nu >= 0) && (mu >= 0)) + { + locale[iat][l][n][spin](m0_all, m1_all) += (srho[irc]).real() / 4.0; + } - if ((nu_prime >= 0) && (mu_prime >= 0)) - locale[iat][l][n][spin](m0_all, m1_all) - += (std::conj(srho[irc_prime])).real() / 4.0; + if ((nu_prime >= 0) && (mu_prime >= 0)) + { + locale[iat][l][n][spin](m0_all, m1_all) + += (std::conj(srho[irc_prime])).real() / 4.0; + } } // ipol1 } // m1 } // ipol0 @@ -277,9 +291,10 @@ void DFTU::cal_occup_m_k(const int iter, const int NL = ucell.atoms[it].nwl + 1; const int LC = orbital_corr[it]; - if (LC == -1) { - continue; -} + if (LC == -1) + { + continue; + } for (int ia = 0; ia < ucell.atoms[it].na; ia++) { @@ -287,19 +302,21 @@ void DFTU::cal_occup_m_k(const int iter, for (int l = 0; l < NL; l++) { - if (l != orbital_corr[it]) { - continue; -} + if (l != orbital_corr[it]) + { + continue; + } const int N = ucell.atoms[it].l_nchi[l]; for (int n = 0; n < N; n++) { // if(!Yukawa && n!=0) continue; - if (n != 0) { - continue; -} - // set the local occupation mumber matrix of spin up and down zeros + if (n != 0) + { + continue; + } + // set the local occupation mumber matrix of spin up and down zeros #ifdef __MPI if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 4) @@ -418,26 +435,29 @@ void DFTU::cal_occup_m_gamma(const int iter, for (int it = 0; it < ucell.ntype; it++) { const int NL = ucell.atoms[it].nwl + 1; - if (orbital_corr[it] == -1) { - continue; -} - for (int ia = 0; ia < ucell.atoms[it].na; ia++) + if (orbital_corr[it] == -1) + { + continue; + } + for (int ia = 0; ia < ucell.atoms[it].na; ia++) { const int iat = ucell.itia2iat(it, ia); for (int l = 0; l < NL; l++) { - if (l != orbital_corr[it]) { - continue; -} + if (l != orbital_corr[it]) + { + continue; + } const int N = ucell.atoms[it].l_nchi[l]; for (int n = 0; n < N; n++) { - if (n != 0) { - continue; -} + if (n != 0) + { + continue; + } // Calculate the local occupation number matrix for (int m0 = 0; m0 < 2 * l + 1; m0++) @@ -524,4 +544,4 @@ void DFTU::cal_occup_m_gamma(const int iter, return; } #endif -} // namespace ModuleDFTU \ No newline at end of file +} // namespace ModuleDFTU diff --git a/source/module_io/write_dos_lcao.cpp b/source/module_io/write_dos_lcao.cpp index 4b4deafb66c..62088e5ffd9 100644 --- a/source/module_io/write_dos_lcao.cpp +++ b/source/module_io/write_dos_lcao.cpp @@ -446,15 +446,25 @@ void ModuleIO::write_dos_lcao(const UnitCell& ucell, { // calculate SK for current k point const std::complex* sk = nullptr; + + // collumn-major matrix + const int hk_type = 1; + if (PARAM.inp.nspin == 4) { - dynamic_cast, std::complex>*>(p_ham)->updateSk(ik, 1); - sk = dynamic_cast, std::complex>*>(p_ham)->getSk(); + dynamic_cast, std::complex>*>(p_ham) + ->updateSk(ik, hk_type); + + sk = dynamic_cast, std::complex>*>(p_ham) + ->getSk(); } else { - dynamic_cast, double>*>(p_ham)->updateSk(ik, 1); - sk = dynamic_cast, double>*>(p_ham)->getSk(); + dynamic_cast, double>*>(p_ham) + ->updateSk(ik, hk_type); + + sk = dynamic_cast, double>*>(p_ham) + ->getSk(); } psi->fix_k(ik); @@ -545,128 +555,125 @@ void ModuleIO::write_dos_lcao(const UnitCell& ucell, if (GlobalV::MY_RANK == 0) { - { - std::stringstream ps; - ps << PARAM.globalv.global_out_dir << "TDOS"; - std::ofstream out(ps.str().c_str()); - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 4) - { - - for (int n = 0; n < npoints; ++n) - { - double y = 0.0; - double en = emin + n * de_ev; - for (int i = 0; i < PARAM.globalv.nlocal; i++) - { - y += pdos[0](i, n); - } - - out << std::setw(20) << en << std::setw(30) << y << std::endl; - } - } - else if (PARAM.inp.nspin == 2) - { - for (int n = 0; n < npoints; ++n) - { - double y = 0.0; - double z = 0.0; - double en = emin + n * de_ev; - for (int i = 0; i < PARAM.globalv.nlocal; i++) - { - y += pdos[0](i, n); - z += pdos[1](i, n); - } - - out << std::setw(20) << en << std::setw(30) << y << std::setw(30) << z << std::endl; - } - } - out.close(); - } + std::stringstream ps; + ps << PARAM.globalv.global_out_dir << "TDOS"; + std::ofstream ofs1(ps.str().c_str()); + if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 4) + { + + for (int n = 0; n < npoints; ++n) + { + double y = 0.0; + double en = emin + n * de_ev; + for (int i = 0; i < PARAM.globalv.nlocal; i++) + { + y += pdos[0](i, n); + } + + ofs1 << std::setw(20) << en << std::setw(30) << y << std::endl; + } + } + else if (PARAM.inp.nspin == 2) + { + for (int n = 0; n < npoints; ++n) + { + double y = 0.0; + double z = 0.0; + double en = emin + n * de_ev; + for (int i = 0; i < PARAM.globalv.nlocal; i++) + { + y += pdos[0](i, n); + z += pdos[1](i, n); + } + + ofs1 << std::setw(20) << en << std::setw(30) << y << std::setw(30) << z << std::endl; + } + } + ofs1.close(); /* decomposed Mulliken charge */ - { - std::stringstream as; - as << PARAM.globalv.global_out_dir << "PDOS"; - std::ofstream out(as.str().c_str()); - - out << "" << std::endl; - out << "" << PARAM.inp.nspin << "" << std::endl; - if (PARAM.inp.nspin == 4) - { - out << "" << std::setw(2) << PARAM.globalv.nlocal / 2 << "" << std::endl; - } - else - { - out << "" << std::setw(2) << PARAM.globalv.nlocal << "" << std::endl; - } - out << "" << std::endl; - - for (int n = 0; n < npoints; ++n) - { - double y = 0.0; - double en = emin + n * de_ev; - out << std::setw(20) << en << std::endl; - } - out << "" << std::endl; - for (int i = 0; i < ucell.nat; i++) - { - int a = ucell.iat2ia[i]; - int t = ucell.iat2it[i]; - Atom* atom1 = &ucell.atoms[t]; - const int s0 = ucell.itiaiw2iwt(t, a, 0); - for (int j = 0; j < atom1->nw; ++j) - { - const int L1 = atom1->iw2l[j]; - const int N1 = atom1->iw2n[j]; - const int m1 = atom1->iw2m[j]; - const int w = ucell.itiaiw2iwt(t, a, j); - - // out << "" <" << std::endl; - out << "" << std::endl; - if (PARAM.inp.nspin == 1) - { - for (int n = 0; n < npoints; ++n) - { - out << std::setw(13) << pdos[0](w, n) << std::endl; - } // n - } - else if (PARAM.inp.nspin == 2) - { - for (int n = 0; n < npoints; ++n) - { - out << std::setw(20) << pdos[0](w, n) << std::setw(30) << pdos[1](w, n) << std::endl; - } // n - } - else if (PARAM.inp.nspin == 4) - { - int w0 = w - s0; - for (int n = 0; n < npoints; ++n) - { - out << std::setw(20) << pdos[0](s0 + 2 * w0, n) + pdos[0](s0 + 2 * w0 + 1, n) << std::endl; - } // n - } - - out << "" << std::endl; - out << "" << std::endl; - } // j - } // i - - out << "" << std::endl; - out.close(); - } - ModuleIO::write_orb_info(&(ucell)); + std::stringstream as; + as << PARAM.globalv.global_out_dir << "PDOS"; + std::ofstream ofs2(as.str().c_str()); + + ofs2 << "" << std::endl; + ofs2 << "" << PARAM.inp.nspin << "" << std::endl; + if (PARAM.inp.nspin == 4) + { + ofs2 << "" << std::setw(2) << PARAM.globalv.nlocal / 2 << "" << std::endl; + } + else + { + ofs2 << "" << std::setw(2) << PARAM.globalv.nlocal << "" << std::endl; + } + ofs2 << "" << std::endl; + + for (int n = 0; n < npoints; ++n) + { + double y = 0.0; + double en = emin + n * de_ev; + ofs2 << std::setw(20) << en << std::endl; + } + ofs2 << "" << std::endl; + for (int i = 0; i < ucell.nat; i++) + { + int a = ucell.iat2ia[i]; + int t = ucell.iat2it[i]; + Atom* atom1 = &ucell.atoms[t]; + const int s0 = ucell.itiaiw2iwt(t, a, 0); + for (int j = 0; j < atom1->nw; ++j) + { + const int L1 = atom1->iw2l[j]; + const int N1 = atom1->iw2n[j]; + const int m1 = atom1->iw2m[j]; + const int w = ucell.itiaiw2iwt(t, a, j); + + // ofs2 << "" <" << std::endl; + ofs2 << "" << std::endl; + if (PARAM.inp.nspin == 1) + { + for (int n = 0; n < npoints; ++n) + { + ofs2 << std::setw(13) << pdos[0](w, n) << std::endl; + } // n + } + else if (PARAM.inp.nspin == 2) + { + for (int n = 0; n < npoints; ++n) + { + ofs2 << std::setw(20) << pdos[0](w, n) << std::setw(30) << pdos[1](w, n) << std::endl; + } // n + } + else if (PARAM.inp.nspin == 4) + { + int w0 = w - s0; + for (int n = 0; n < npoints; ++n) + { + ofs2 << std::setw(20) << pdos[0](s0 + 2 * w0, n) + pdos[0](s0 + 2 * w0 + 1, n) << std::endl; + } // n + } + + ofs2 << "" << std::endl; + ofs2 << "" << std::endl; + } // j + } // i + + ofs2 << "" << std::endl; + ofs2.close(); + ModuleIO::write_orb_info(&(ucell)); } delete[] pdos; - } + }// end PARAM.inp.out_dos == 2 + // output the DOS file. for (int is = 0; is < nspin0; ++is) {