diff --git a/source/module_cell/module_neighbor/sltk_grid.h b/source/module_cell/module_neighbor/sltk_grid.h index 6a14b14f8b..56a8b6ee3c 100644 --- a/source/module_cell/module_neighbor/sltk_grid.h +++ b/source/module_cell/module_neighbor/sltk_grid.h @@ -1,23 +1,23 @@ #ifndef GRID_H #define GRID_H -#include -#include -#include "sltk_util.h" +#include "module_cell/unitcell.h" #include "sltk_atom.h" #include "sltk_atom_input.h" +#include "sltk_util.h" -#include "module_cell/unitcell.h" -#include +#include +#include #include +#include typedef std::vector> AtomMap; struct CellSet { - AtomMap atom_map; - int in_grid[3]; - CellSet(); + AtomMap atom_map; + int in_grid[3]; + CellSet(); }; //========================================================== @@ -34,111 +34,132 @@ class Atom_input; class Grid { -public: - - // Constructors and destructor - // Grid is Global class,so init it with constant number - Grid():test_grid(0){}; - Grid(const int &test_grid_in); - virtual ~Grid(); - - void init(std::ofstream &ofs, - const UnitCell &ucell, - const Atom_input &input); - - //2015-05-07 - void delete_vector(int i, int j, int k); - - //Data - bool pbc; // periodic boundary condition - bool expand_flag; - double sradius2;// searching radius squared - double sradius;// searching radius - double d_minX;// origin of all cells - double d_minY; - double d_minZ; - int cell_nx; - int cell_ny; - int cell_nz; - int layer; - - int true_cell_x; - int true_cell_y; - int true_cell_z; - - std::vector>> Cell; //dx , dy ,dz is cell number in each direction,respectly. - void delete_Cell() //it will replace by container soon! - { - if (this->init_cell_flag) - { - for (int i = 0;i < this->cell_nx;i++) - { - for (int j = 0;j < this->cell_ny;j++) - { - this->Cell[i][j].clear(); - } - } - - for (int i = 0;i < this->cell_nx;i++) - { - this->Cell[i].clear(); - } - - this->Cell.clear(); - this->init_cell_flag = false; - } - } - bool init_cell_flag; - //LiuXh add 2019-07-15 - double getD_minX() const {return d_minX;} - double getD_minY() const {return d_minY;} - double getD_minZ() const {return d_minZ;} - - - int getCellX() const {return cell_nx;} - int getCellY() const {return cell_ny;} - int getCellZ() const {return cell_nz;} - int getTrueCellX() const {return true_cell_x;} - int getTrueCellY() const {return true_cell_y;} - int getTrueCellZ() const {return true_cell_z;} - -private: - - const int test_grid; -//========================================================== -// MEMBER FUNCTIONS : -// Three Main Steps: -// NAME : setMemberVariables (read in datas from Atom_input, -// init cells.) -// NAME : setBoundaryAdjacent( Consider different situations, -// if not_expand case : nature/periodic boundary -// condition , if expand_case) -//========================================================== - void setMemberVariables(std::ofstream &ofs_in, - const Atom_input &input); - - void setBoundaryAdjacent(std::ofstream &ofs_in, - const Atom_input &input); - -//========================================================== - void Build_Hash_Table(const UnitCell &ucell, const Atom_input &input); - -//========================================================== - - void Construct_Adjacent_expand(const int i, const int j, const int k); - - void Construct_Adjacent_expand_periodic(const int i, - const int j, - const int k, - FAtom& fatom); - - void Construct_Adjacent_begin(); - void Construct_Adjacent_nature( - const int i, const int j, const int k, FAtom & fatom1); - void Construct_Adjacent_periodic( - const int i, const int j, const int k, FAtom & fatom1); - void Construct_Adjacent_final(const int i, const int j, const int k, FAtom & fatom1, - const int i2, const int j2, const int k2, FAtom & fatom2); + public: + // Constructors and destructor + // Grid is Global class,so init it with constant number + Grid() : test_grid(0){}; + Grid(const int& test_grid_in); + virtual ~Grid(); + + void init(std::ofstream& ofs, const UnitCell& ucell, const Atom_input& input); + + // 2015-05-07 + void delete_vector(int i, int j, int k); + + // Data + bool pbc; // periodic boundary condition + bool expand_flag; + double sradius2; // searching radius squared + double sradius; // searching radius + double d_minX; // origin of all cells + double d_minY; + double d_minZ; + int cell_nx; + int cell_ny; + int cell_nz; + int layer; + + int true_cell_x; + int true_cell_y; + int true_cell_z; + + std::vector>> Cell; // dx , dy ,dz is cell number in each direction,respectly. + void delete_Cell() // it will replace by container soon! + { + if (this->init_cell_flag) + { + for (int i = 0; i < this->cell_nx; i++) + { + for (int j = 0; j < this->cell_ny; j++) + { + this->Cell[i][j].clear(); + } + } + + for (int i = 0; i < this->cell_nx; i++) + { + this->Cell[i].clear(); + } + + this->Cell.clear(); + this->init_cell_flag = false; + } + } + bool init_cell_flag = false; + // LiuXh add 2019-07-15 + double getD_minX() const + { + return d_minX; + } + double getD_minY() const + { + return d_minY; + } + double getD_minZ() const + { + return d_minZ; + } + + int getCellX() const + { + return cell_nx; + } + int getCellY() const + { + return cell_ny; + } + int getCellZ() const + { + return cell_nz; + } + int getTrueCellX() const + { + return true_cell_x; + } + int getTrueCellY() const + { + return true_cell_y; + } + int getTrueCellZ() const + { + return true_cell_z; + } + + private: + const int test_grid; + //========================================================== + // MEMBER FUNCTIONS : + // Three Main Steps: + // NAME : setMemberVariables (read in datas from Atom_input, + // init cells.) + // NAME : setBoundaryAdjacent( Consider different situations, + // if not_expand case : nature/periodic boundary + // condition , if expand_case) + //========================================================== + void setMemberVariables(std::ofstream& ofs_in, const Atom_input& input); + + void setBoundaryAdjacent(std::ofstream& ofs_in, const Atom_input& input); + + //========================================================== + void Build_Hash_Table(const UnitCell& ucell, const Atom_input& input); + + //========================================================== + + void Construct_Adjacent_expand(const int i, const int j, const int k); + + void Construct_Adjacent_expand_periodic(const int i, const int j, const int k, FAtom& fatom); + + void Construct_Adjacent_begin(); + void Construct_Adjacent_nature(const int i, const int j, const int k, FAtom& fatom1); + void Construct_Adjacent_periodic(const int i, const int j, const int k, FAtom& fatom1); + void Construct_Adjacent_final(const int i, + const int j, + const int k, + FAtom& fatom1, + const int i2, + const int j2, + const int k2, + FAtom& fatom2); }; #endif diff --git a/source/module_cell/module_neighbor/sltk_grid_driver.h b/source/module_cell/module_neighbor/sltk_grid_driver.h index 457dfdb405..fbc683cade 100644 --- a/source/module_cell/module_neighbor/sltk_grid_driver.h +++ b/source/module_cell/module_neighbor/sltk_grid_driver.h @@ -1,16 +1,17 @@ #ifndef GRID_DRIVER_H #define GRID_DRIVER_H -#include -#include -#include "sltk_atom.h" -#include "sltk_atom_input.h" -#include "sltk_grid.h" #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_base/vector3.h" -#include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" #include "module_cell/unitcell.h" +#include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" +#include "sltk_atom.h" +#include "sltk_atom_input.h" +#include "sltk_grid.h" + +#include +#include #include //========================================================== @@ -18,91 +19,111 @@ //========================================================== class AdjacentAtomInfo { -public: - AdjacentAtomInfo() : adj_num(0) {} - int adj_num; - std::vector ntype; - std::vector natom; - std::vector> adjacent_tau; - std::vector> box; - void clear() - { - adj_num = 0; - ntype.clear(); - natom.clear(); - adjacent_tau.clear(); - box.clear(); - } + public: + AdjacentAtomInfo() : adj_num(0) + { + } + int adj_num; + std::vector ntype; + std::vector natom; + std::vector> adjacent_tau; + std::vector> box; + void clear() + { + adj_num = 0; + ntype.clear(); + natom.clear(); + adjacent_tau.clear(); + box.clear(); + } }; void filter_adjs(const std::vector& is_adj, AdjacentAtomInfo& adjs); class Grid_Driver : public Grid { -public: + public: + //========================================================== + // THE INTERFACE WITH USER : + // MEMBRE FUNCTIONS : + // NAME : Find_atom (input cartesian position,find the + // adjacent of this atom,and store the information + // in 'adj_num','ntype','natom' + //========================================================== + Grid_Driver() : test_deconstructor(0){}; + Grid_Driver(const int& test_d_in, const int& test_grid_in); - //========================================================== - // THE INTERFACE WITH USER : - // MEMBRE FUNCTIONS : - // NAME : Find_atom (input cartesian position,find the - // adjacent of this atom,and store the information - // in 'adj_num','ntype','natom' - //========================================================== - Grid_Driver():test_deconstructor(0){}; - Grid_Driver( - const int &test_d_in, - const int &test_grid_in); + ~Grid_Driver(); - ~Grid_Driver(); + //========================================================== + // EXPLAIN FOR default parameter `adjs = nullptr` + // + // This design make Grid_Driver compatible with multi-thread usage + // 1. Find_atom store results in Grid_Driver::adj_info + // by default. + // 2. And store results into parameter adjs when adjs is + // NOT NULL + //========================================================== + void Find_atom(const UnitCell& ucell, + const ModuleBase::Vector3& cartesian_posi, + const int& ntype, + const int& nnumber, + AdjacentAtomInfo* adjs = nullptr); - //========================================================== - // EXPLAIN FOR default parameter `adjs = nullptr` - // - // This design make Grid_Driver compatible with multi-thread usage - // 1. Find_atom store results in Grid_Driver::adj_info - // by default. - // 2. And store results into parameter adjs when adjs is - // NOT NULL - //========================================================== - void Find_atom( - const UnitCell &ucell, - const ModuleBase::Vector3 &cartesian_posi, - const int &ntype, - const int &nnumber, - AdjacentAtomInfo *adjs = nullptr); + //========================================================== + // EXPLAIN : The adjacent information for the input + // cartesian_pos + // MEMBER VARIABLES : + // NAME : getAdjacentNum + // NAME : getNtype + // NAME : getNatom + // NAME : getAdjaentTau + //========================================================== + const int& getAdjacentNum() const + { + return adj_info.adj_num; + } + const int& getType(const int i) const + { + return adj_info.ntype[i]; + } + const int& getNatom(const int i) const + { + return adj_info.natom[i]; + } + const ModuleBase::Vector3& getAdjacentTau(const int i) const + { + return adj_info.adjacent_tau[i]; + } + const ModuleBase::Vector3& getBox(const int i) const + { + return adj_info.box[i]; + } - //========================================================== - // EXPLAIN : The adjacent information for the input - // cartesian_pos - // MEMBER VARIABLES : - // NAME : getAdjacentNum - // NAME : getNtype - // NAME : getNatom - // NAME : getAdjaentTau - //========================================================== - const int& getAdjacentNum()const { return adj_info.adj_num; } - const int& getType(const int i) const { return adj_info.ntype[i]; } - const int& getNatom(const int i) const { return adj_info.natom[i]; } - const ModuleBase::Vector3& getAdjacentTau(const int i) const { return adj_info.adjacent_tau[i]; } - const ModuleBase::Vector3& getBox(const int i) const {return adj_info.box[i];} + private: + mutable AdjacentAtomInfo adj_info; -private: + const int test_deconstructor; // caoyu reconst 2021-05-24 - mutable AdjacentAtomInfo adj_info; - - const int test_deconstructor;//caoyu reconst 2021-05-24 - -//========================================================== -// MEMBER FUNCTIONS : -// NAME : Calculate_adjacent_site -//========================================================== - ModuleBase::Vector3 Calculate_adjacent_site(const double x, const double y, const double z, - const double &box11, const double &box12, const double &box13, - const double &box21, const double &box22, const double &box23, - const double &box31, const double &box32, const double &box33, - const short box_x, // three dimensions of the target box - const short box_y, - const short box_z) const; + //========================================================== + // MEMBER FUNCTIONS : + // NAME : Calculate_adjacent_site + //========================================================== + ModuleBase::Vector3 Calculate_adjacent_site(const double x, + const double y, + const double z, + const double& box11, + const double& box12, + const double& box13, + const double& box21, + const double& box22, + const double& box23, + const double& box31, + const double& box32, + const double& box33, + const short box_x, // three dimensions of the target box + const short box_y, + const short box_z) const; }; namespace GlobalC diff --git a/source/module_esolver/cell_dependency.h b/source/module_esolver/cell_dependency.h new file mode 100644 index 0000000000..0610b9f33d --- /dev/null +++ b/source/module_esolver/cell_dependency.h @@ -0,0 +1,20 @@ +#ifndef CELL_DEPENDENCY +#define CELL_DEPENDENCY + +#include "module_cell/module_neighbor/sltk_atom_arrange.h" +#include "module_hamilt_lcao/hamilt_lcaodft/record_adj.h" +#include "module_lr/utils/lr_util.hpp" + +#include + +class Cell_Dependency +{ + public: + Cell_Dependency() {}; + ~Cell_Dependency() {}; + + Record_adj ra; + Grid_Driver grid_d; +}; + +#endif // CELL_DEPENDENCY diff --git a/source/module_esolver/esolver_gets.cpp b/source/module_esolver/esolver_gets.cpp index 3571d4788f..5a39566c49 100644 --- a/source/module_esolver/esolver_gets.cpp +++ b/source/module_esolver/esolver_gets.cpp @@ -6,9 +6,9 @@ #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" +#include "module_io/cal_r_overlap_R.h" #include "module_io/print_info.h" #include "module_io/write_HS_R.h" -#include "module_io/cal_r_overlap_R.h" namespace ModuleESolver { @@ -39,7 +39,7 @@ void ESolver_GetS::before_all_runners(UnitCell& ucell, const Input_para& inp) } // 1.3) Setup k-points according to symmetry. - this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec,GlobalV::ofs_running); + this->kv.set(ucell, ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); ModuleIO::setup_parameters(ucell, this->kv); @@ -101,7 +101,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) PARAM.inp.test_atom_input); Record_adj RA; - RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); if (this->p_hamilt == nullptr) { @@ -109,6 +109,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) { this->p_hamilt = new hamilt::HamiltLCAO, std::complex>(ucell, + GlobalC::GridD, &this->pv, this->kv, *(two_center_bundle_.overlap_orb), @@ -119,6 +120,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) else { this->p_hamilt = new hamilt::HamiltLCAO, double>(ucell, + GlobalC::GridD, &this->pv, this->kv, *(two_center_bundle_.overlap_orb), diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index ecde5f1178..2d94217106 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -1,9 +1,9 @@ +#include "module_elecstate/cal_ux.h" #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_esolver/esolver_ks_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_elecstate/cal_ux.h" // #include "module_base/timer.h" #include "module_cell/module_neighbor/sltk_atom_arrange.h" @@ -127,7 +127,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) // (2)For each atom, calculate the adjacent atoms in different cells // and allocate the space for H(R) and S(R). // If k point is used here, allocate HlocR after atom_arrange. - this->RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation @@ -189,6 +189,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), ucell, + GlobalC::GridD, &this->pv, this->pelec->pot, this->kv, @@ -311,11 +312,9 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) // two cases are considered: // 1. DMK in DensityMatrix is not empty (istep > 0), then DMR is initialized by DMK // 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros - if(istep > 0) + if (istep > 0) { - dynamic_cast*>(this->pelec) - ->get_DM() - ->cal_DMR(); + dynamic_cast*>(this->pelec)->get_DM()->cal_DMR(); } if (PARAM.inp.dm_to_rho) @@ -370,16 +369,18 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) this->p_hamilt->non_first_scf = istep; // update in ion-step - if( PARAM.inp.rdmft == true ) + if (PARAM.inp.rdmft == true) { // necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp - rdmft_solver.update_ion(ucell, *(this->pw_rho), this->ppcell.vloc, this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08 + rdmft_solver.update_ion(ucell, + *(this->pw_rho), + this->ppcell.vloc, + this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08 } return; } - template class ESolver_KS_LCAO; template class ESolver_KS_LCAO, double>; template class ESolver_KS_LCAO, std::complex>; diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 6ce865836d..c5c2aae021 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -1,9 +1,9 @@ +#include "module_elecstate/cal_ux.h" #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_esolver/esolver_ks_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -#include "module_elecstate/cal_ux.h" // #include "module_base/timer.h" #include "module_cell/module_neighbor/sltk_atom_arrange.h" @@ -128,7 +128,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) // (2)For each atom, calculate the adjacent atoms in different cells // and allocate the space for H(R) and S(R). // If k point is used here, allocate HlocR after atom_arrange. - this->RA.for_2d(ucell,this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(ucell, GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation @@ -175,7 +175,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } // prepare grid in Gint - LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell,orb_, *this->pw_rho, *this->pw_big); + LCAO_domain::grid_prepare(this->GridT, this->GG, this->GK, ucell, orb_, *this->pw_rho, *this->pw_big); // init Hamiltonian if (this->p_hamilt != nullptr) @@ -190,6 +190,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), ucell, + GlobalC::GridD, &this->pv, this->pelec->pot, this->kv, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index ca76150f17..059d82c1ed 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -1,11 +1,12 @@ #include "hamilt_lcao.h" #include "module_base/global_variable.h" -#include "module_parameter/parameter.h" #include "module_base/memory.h" #include "module_base/timer.h" #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_parameter/parameter.h" + #include #ifdef __DEEPKS @@ -24,17 +25,17 @@ #include "module_elecstate/potentials/H_TDDFT_pw.h" #include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h" #include "module_hsolver/hsolver_lcao.h" #include "operator_lcao/dftu_lcao.h" +#include "operator_lcao/dspin_lcao.h" #include "operator_lcao/ekinetic_new.h" #include "operator_lcao/meta_lcao.h" #include "operator_lcao/nonlocal_new.h" #include "operator_lcao/op_dftu_lcao.h" #include "operator_lcao/op_exx_lcao.h" #include "operator_lcao/overlap_new.h" -#include "operator_lcao/dspin_lcao.h" -#include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" #include "operator_lcao/td_ekinetic_lcao.h" #include "operator_lcao/td_nonlocal_lcao.h" #include "operator_lcao/veff_lcao.h" @@ -44,9 +45,10 @@ namespace hamilt template HamiltLCAO::HamiltLCAO(const UnitCell& ucell, - const Parallel_Orbitals* paraV, - const K_Vectors& kv_in, - const TwoCenterIntegrator& intor_overlap_orb, + Grid_Driver& grid_d, + const Parallel_Orbitals* paraV, + const K_Vectors& kv_in, + const TwoCenterIntegrator& intor_overlap_orb, const std::vector& orb_cutoff) { this->classname = "HamiltLCAO"; @@ -64,27 +66,30 @@ HamiltLCAO::HamiltLCAO(const UnitCell& ucell, this->sR, &ucell, orb_cutoff, - &GlobalC::GridD, + &grid_d, &intor_overlap_orb); } template HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, - Gint_k* GK_in, - const UnitCell& ucell, - 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 + Gint_k* GK_in, + const UnitCell& ucell, + Grid_Driver& grid_d, + 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 __EXX - , const int istep - , int* exx_two_level_step - , std::vector>>>* Hexxd - , std::vector>>>>* Hexxc + , + const int istep, + int* exx_two_level_step, + std::vector>>>* Hexxd, + std::vector>>>>* Hexxc #endif -) { +) +{ this->kv = &kv_in; this->classname = "HamiltLCAO"; @@ -140,7 +145,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->sR, &ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb.get()); // kinetic term () @@ -151,7 +156,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, &ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.kinetic_orb.get()); this->getOperator()->add(ekinetic); } @@ -165,7 +170,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, &ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_beta.get()); this->getOperator()->add(nonlocal); } @@ -181,15 +186,14 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, pot_in->pot_register(pot_register_in); // effective potential term Operator* veff = new Veff>(GG_in, - this->hsk, - this->kv->kvec_d, - pot_in, - this->hR, // no explicit call yet - &ucell, - orb.cutoffs(), - &GlobalC::GridD, - PARAM.inp.nspin - ); + this->hsk, + this->kv->kvec_d, + pot_in, + this->hR, // no explicit call yet + &ucell, + orb.cutoffs(), + &GlobalC::GridD, + PARAM.inp.nspin); this->getOperator()->add(veff); } } @@ -201,7 +205,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, // no explicit call yet &ucell, - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_alpha.get(), &orb, this->kv->get_nks(), @@ -227,7 +231,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, ucell, - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), &GlobalC::dftu); @@ -250,14 +254,14 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, pot_in->pot_register(pot_register_in); // Veff term this->getOperator() = new Veff>(GK_in, - this->hsk, - kv->kvec_d, - pot_in, - this->hR, - &ucell, - orb.cutoffs(), - &GlobalC::GridD, - PARAM.inp.nspin); + this->hsk, + kv->kvec_d, + pot_in, + this->hR, + &ucell, + orb.cutoffs(), + &grid_d, + PARAM.inp.nspin); // reset spin index and real space Hamiltonian matrix int start_spin = -1; GK_in->reset_spin(start_spin); @@ -272,7 +276,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->sR, &ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb.get()); if (this->getOperator() == nullptr) { @@ -292,7 +296,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, &ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.kinetic_orb.get()); this->getOperator()->add(ekinetic); } @@ -306,7 +310,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, &ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_beta.get()); // TDDFT velocity gague will calculate full non-local potential including the original one and the // correction on its own. So the original non-local potential term should be skipped @@ -327,7 +331,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, hR, &ucell, - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_alpha.get(), &orb, this->kv->get_nks(), @@ -338,23 +342,21 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, // TDDFT_velocity_gague if (TD_Velocity::tddft_velocity) { - if(!TD_Velocity::init_vecpot_file) { elecstate::H_TDDFT_pw::update_At(); -} + if (!TD_Velocity::init_vecpot_file) + { + elecstate::H_TDDFT_pw::update_At(); + } Operator* td_ekinetic = new TDEkinetic>(this->hsk, this->hR, kv, &ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb.get()); this->getOperator()->add(td_ekinetic); - Operator* td_nonlocal = new TDNonlocal>(this->hsk, - this->kv->kvec_d, - this->hR, - &ucell, - orb, - &GlobalC::GridD); + Operator* td_nonlocal + = new TDNonlocal>(this->hsk, this->kv->kvec_d, this->hR, &ucell, orb, &grid_d); this->getOperator()->add(td_nonlocal); } if (PARAM.inp.dft_plus_u) @@ -373,7 +375,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, ucell, - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), &GlobalC::dftu); @@ -382,16 +384,13 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, } if (PARAM.inp.sc_mag_switch) { - Operator* sc_lambda = - new DeltaSpin>( - this->hsk, - kv->kvec_d, - this->hR, - ucell, - &GlobalC::GridD, - two_center_bundle.overlap_orb_onsite.get(), - orb.cutoffs() - ); + Operator* sc_lambda = new DeltaSpin>(this->hsk, + kv->kvec_d, + this->hR, + ucell, + &grid_d, + two_center_bundle.overlap_orb_onsite.get(), + orb.cutoffs()); this->getOperator()->add(sc_lambda); spinconstrain::SpinConstrain& sc = spinconstrain::SpinConstrain::getScInstance(); sc.set_operator(sc_lambda); @@ -405,15 +404,15 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, // set xc type before the first cal of xc in pelec->init_scf // and calculate Cs, Vs Operator* exx = new OperatorEXX>(this->hsk, - this->hR, - *this->kv, - Hexxd, - Hexxc, - Add_Hexx_Type::R, - istep, - exx_two_level_step, - !GlobalC::restart.info_load.restart_exx - && GlobalC::restart.info_load.load_H); + this->hR, + *this->kv, + Hexxd, + Hexxc, + Add_Hexx_Type::R, + istep, + exx_two_level_step, + !GlobalC::restart.info_load.restart_exx + && GlobalC::restart.info_load.load_H); this->getOperator()->add(exx); } #endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h index 0c8a5bc188..ee05e59d34 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h @@ -3,13 +3,15 @@ #include "module_basis/module_nao/two_center_bundle.h" #include "module_cell/klist.h" +#include "module_cell/module_neighbor/sltk_atom_arrange.h" #include "module_elecstate/module_dm/density_matrix.h" #include "module_elecstate/potentials/potential_new.h" #include "module_hamilt_general/hamilt.h" +#include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" -#include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" + #include #ifdef __EXX #include "module_ri/Exx_LRI.h" @@ -27,30 +29,33 @@ class HamiltLCAO : public Hamilt * @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 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 + using TAC = std::pair>; + HamiltLCAO(Gint_Gamma* GG_in, + Gint_k* GK_in, + const UnitCell& ucell, + Grid_Driver& grid_d, + 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 __EXX - , const int istep - , int* exx_two_level_step = nullptr - , std::vector>>>* Hexxd = nullptr - , std::vector>>>>* Hexxc = nullptr + , + const int istep, + int* exx_two_level_step = nullptr, + std::vector>>>* Hexxd = nullptr, + std::vector>>>>* Hexxc = nullptr #endif - ); + ); /** * @brief Constructor of vacuum Operators, only HR and SR will be initialed as empty HContainer */ - HamiltLCAO(const UnitCell& ucell, - const Parallel_Orbitals* paraV, - const K_Vectors& kv_in, - const TwoCenterIntegrator& intor_overlap_orb, + HamiltLCAO(const UnitCell& ucell, + Grid_Driver& grid_d, + const Parallel_Orbitals* paraV, + const K_Vectors& kv_in, + const TwoCenterIntegrator& intor_overlap_orb, const std::vector& orb_cutoff); ~HamiltLCAO() @@ -67,14 +72,29 @@ class HamiltLCAO : public Hamilt /// get pointer of Operator ops Operator*& getOperator(); /// get hk-pointer - TK* getHk() const{return this->hsk->get_hk();} + TK* getHk() const + { + return this->hsk->get_hk(); + } /// get sk-pointer - TK* getSk() const{return this->hsk->get_sk();} - int get_size_hsk() const{return this->hsk->get_size();} + 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;} + HContainer*& getHR() + { + return this->hR; + } /// get SR pointer of *this->sR, which is a HContainer and contains S(R) - HContainer*& getSR(){return this->sR;} + HContainer*& getSR() + { + return this->sR; + } /// refresh the status of HR void refresh() override; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp index 695e6327db..b48de89998 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp @@ -1,483 +1,507 @@ #include "record_adj.h" -#include "module_parameter/parameter.h" #include "module_base/timer.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -Record_adj::Record_adj() : iat2ca(nullptr) {} -Record_adj::~Record_adj(){ - if(info_modified) { - this->delete_grid(); +#include "module_parameter/parameter.h" +Record_adj::Record_adj() : iat2ca(nullptr) +{ } +Record_adj::~Record_adj() +{ + if (info_modified) + { + this->delete_grid(); + } } void Record_adj::delete_grid() { - for(int i=0; i& orb_cutoff) +void Record_adj::for_2d(const UnitCell& ucell, + Grid_Driver& grid_d, + Parallel_Orbitals& pv, + bool gamma_only, + const std::vector& orb_cutoff) { - ModuleBase::TITLE("Record_adj","for_2d"); - ModuleBase::timer::tick("Record_adj","for_2d"); + ModuleBase::TITLE("Record_adj", "for_2d"); + ModuleBase::timer::tick("Record_adj", "for_2d"); assert(ucell.nat > 0); if (!gamma_only) { delete[] pv.nlocdim; delete[] pv.nlocstart; - pv.nlocdim = new int[ucell.nat]; + pv.nlocdim = new int[ucell.nat]; pv.nlocstart = new int[ucell.nat]; ModuleBase::GlobalFunc::ZEROS(pv.nlocdim, ucell.nat); ModuleBase::GlobalFunc::ZEROS(pv.nlocstart, ucell.nat); pv.nnr = 0; } -{ - // (1) find the adjacent atoms of atom[T1,I1]; - ModuleBase::Vector3 tau1, tau2, dtau; - ModuleBase::Vector3 dtau1, dtau2, tau0; - - this->na_proc = ucell.nat; - - // number of adjacents for each atom. - this->na_each = new int[na_proc]; - ModuleBase::GlobalFunc::ZEROS(na_each, na_proc); - int iat = 0; - - - for (int T1 = 0; T1 < ucell.ntype; ++T1) - { - Atom* atom1 = &ucell.atoms[T1]; - for (int I1 = 0; I1 < atom1->na; ++I1) - { - tau1 = atom1->tau[I1]; - //GlobalC::GridD.Find_atom( tau1 ); - GlobalC::GridD.Find_atom(ucell, tau1 ,T1, I1); - const int start1 = ucell.itiaiw2iwt(T1, I1, 0); - if(!gamma_only) { pv.nlocstart[iat] = pv.nnr; -} - - // (2) search among all adjacent atoms. - for (int ad = 0; ad < GlobalC::GridD.getAdjacentNum()+1; ++ad) - { - const int T2 = GlobalC::GridD.getType(ad); - const int I2 = GlobalC::GridD.getNatom(ad); - const int start2 = ucell.itiaiw2iwt(T2, I2, 0); - tau2 = GlobalC::GridD.getAdjacentTau(ad); - dtau = tau2 - tau1; - double distance = dtau.norm() * ucell.lat0; - double rcut = orb_cutoff[T1] + orb_cutoff[T2]; - - bool is_adj = false; - if (distance < rcut) { is_adj = true; - // there is another possibility that i and j are adjacent atoms. - // which is that are adjacents while are also - // adjacents, these considerations are only considered in k-point - // algorithm, - } else if (distance >= rcut) - { - for (int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0) - { - const int T0 = GlobalC::GridD.getType(ad0); - //const int I0 = GlobalC::GridD.getNatom(ad0); - //const int iat0 = ucell.itia2iat(T0, I0); - //const int start0 = ucell.itiaiw2iwt(T0, I0, 0); - - tau0 = GlobalC::GridD.getAdjacentTau(ad0); - dtau1 = tau0 - tau1; - double distance1 = dtau1.norm() * ucell.lat0; - double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); - - dtau2 = tau0 - tau2; - double distance2 = dtau2.norm() * ucell.lat0; - double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); - - if( distance1 < rcut1 && distance2 < rcut2 ) - { - is_adj = true; - break; - } // dis1, dis2 - } - } - - if(is_adj) - { - ++na_each[iat]; - if (!gamma_only) + { + // (1) find the adjacent atoms of atom[T1,I1]; + ModuleBase::Vector3 tau1, tau2, dtau; + ModuleBase::Vector3 dtau1, dtau2, tau0; + + this->na_proc = ucell.nat; + + // number of adjacents for each atom. + this->na_each = new int[na_proc]; + ModuleBase::GlobalFunc::ZEROS(na_each, na_proc); + int iat = 0; + + for (int T1 = 0; T1 < ucell.ntype; ++T1) + { + Atom* atom1 = &ucell.atoms[T1]; + for (int I1 = 0; I1 < atom1->na; ++I1) + { + tau1 = atom1->tau[I1]; + // grid_d.Find_atom( tau1 ); + grid_d.Find_atom(ucell, tau1, T1, I1); + const int start1 = ucell.itiaiw2iwt(T1, I1, 0); + if (!gamma_only) + { + pv.nlocstart[iat] = pv.nnr; + } + + // (2) search among all adjacent atoms. + for (int ad = 0; ad < grid_d.getAdjacentNum() + 1; ++ad) + { + const int T2 = grid_d.getType(ad); + const int I2 = grid_d.getNatom(ad); + const int start2 = ucell.itiaiw2iwt(T2, I2, 0); + tau2 = grid_d.getAdjacentTau(ad); + dtau = tau2 - tau1; + double distance = dtau.norm() * ucell.lat0; + double rcut = orb_cutoff[T1] + orb_cutoff[T2]; + + bool is_adj = false; + if (distance < rcut) { - for(int ii=0; iinw * PARAM.globalv.npol; ++ii) + is_adj = true; + // there is another possibility that i and j are adjacent atoms. + // which is that are adjacents while are also + // adjacents, these considerations are only considered in k-point + // algorithm, + } + else if (distance >= rcut) + { + for (int ad0 = 0; ad0 < grid_d.getAdjacentNum() + 1; ++ad0) { - // 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; -} + const int T0 = grid_d.getType(ad0); + // const int I0 = grid_d.getNatom(ad0); + // const int iat0 = ucell.itia2iat(T0, I0); + // const int start0 = ucell.itiaiw2iwt(T0, I0, 0); + + tau0 = grid_d.getAdjacentTau(ad0); + dtau1 = tau0 - tau1; + double distance1 = dtau1.norm() * ucell.lat0; + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); - for(int jj=0; jjinfo = new int**[na_proc]; + if (is_adj) + { + ++na_each[iat]; + if (!gamma_only) + { + for (int ii = 0; ii < atom1->nw * PARAM.globalv.npol; ++ii) + { + // 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; + } + + for (int jj = 0; jj < ucell.atoms[T2].nw * PARAM.globalv.npol; ++jj) + { + const int iw2_all = start2 + jj; + const int nu = pv.global2local_col(iw2_all); + if (nu < 0) + { + continue; + } + + pv.nlocdim[iat]++; + ++(pv.nnr); + } + } + } + } // end is_adj + } // end ad + ++iat; + } // end I1 + } // end T1 + } + // xiaohui add "OUT_LEVEL", 2015-09-16 + if (PARAM.inp.out_level != "m" && !gamma_only) + { + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "ParaV.nnr", pv.nnr); + } + + //------------------------------------------------ + // info will identify each atom in each unitcell. + //------------------------------------------------ + this->info = new int**[na_proc]; #ifdef _OPENMP #pragma omp parallel -{ + { #endif - ModuleBase::Vector3 tau1, tau2, dtau; - ModuleBase::Vector3 dtau1, dtau2, tau0; + ModuleBase::Vector3 tau1, tau2, dtau; + ModuleBase::Vector3 dtau1, dtau2, tau0; #ifdef _OPENMP #pragma omp for schedule(dynamic) #endif - for(int i=0; i 0) - { - info[i] = new int*[ na_each[i] ]; - for(int j=0; j 0) + { + info[i] = new int*[na_each[i]]; + for (int j = 0; j < na_each[i]; j++) + { + // (Rx, Ry, Rz, T, I) + info[i][j] = new int[5]; + ModuleBase::GlobalFunc::ZEROS(info[i][j], 5); + } + } + } #ifdef _OPENMP #pragma omp for schedule(dynamic) #endif - for(int iat=0; iattau[I1]; - //GlobalC::GridD.Find_atom( tau1 ); - AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(ucell, tau1 ,T1, I1, &adjs); - - // (2) search among all adjacent atoms. - int cb = 0; - for (int ad = 0; ad < adjs.adj_num+1; ++ad) - { - const int T2 = adjs.ntype[ad]; - const int I2 = adjs.natom[ad]; - tau2 = adjs.adjacent_tau[ad]; - dtau = tau2 - tau1; - double distance = dtau.norm() * ucell.lat0; - double rcut = orb_cutoff[T1] + orb_cutoff[T2]; - - - bool is_adj = false; - if(distance < rcut) { is_adj = true; - } else if(distance >= rcut) - { - for (int ad0 = 0; ad0 < adjs.adj_num+1; ++ad0) - { - const int T0 = adjs.ntype[ad0]; - //const int I0 = GlobalC::GridD.getNatom(ad0); - //const int iat0 = ucell.itia2iat(T0, I0); - //const int start0 = ucell.itiaiw2iwt(T0, I0, 0); - - tau0 = adjs.adjacent_tau[ad0]; - dtau1 = tau0 - tau1; - double distance1 = dtau1.norm() * ucell.lat0; - double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); - - dtau2 = tau0 - tau2; - double distance2 = dtau2.norm() * ucell.lat0; - double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); - - if( distance1 < rcut1 && distance2 < rcut2 ) - { - is_adj = true; - break; - } // dis1, dis2 - } - } - - if(is_adj) - { - info[iat][cb][0] = adjs.box[ad].x; - info[iat][cb][1] = adjs.box[ad].y; - info[iat][cb][2] = adjs.box[ad].z; - info[iat][cb][3] = T2; - info[iat][cb][4] = I2; - ++cb; - } - }//end ad -// GlobalV::ofs_running << " nadj = " << cb << std::endl; - }//end I1 - }//end T1 + for (int iat = 0; iat < ucell.nat; ++iat) + { + const int T1 = ucell.iat2it[iat]; + Atom* atom1 = &ucell.atoms[T1]; + const int I1 = ucell.iat2ia[iat]; + { + tau1 = atom1->tau[I1]; + // grid_d.Find_atom( tau1 ); + AdjacentAtomInfo adjs; + grid_d.Find_atom(ucell, tau1, T1, I1, &adjs); + + // (2) search among all adjacent atoms. + int cb = 0; + for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + { + const int T2 = adjs.ntype[ad]; + const int I2 = adjs.natom[ad]; + tau2 = adjs.adjacent_tau[ad]; + dtau = tau2 - tau1; + double distance = dtau.norm() * ucell.lat0; + double rcut = orb_cutoff[T1] + orb_cutoff[T2]; + + bool is_adj = false; + if (distance < rcut) + { + is_adj = true; + } + else if (distance >= rcut) + { + for (int ad0 = 0; ad0 < adjs.adj_num + 1; ++ad0) + { + const int T0 = adjs.ntype[ad0]; + // const int I0 = grid_d.getNatom(ad0); + // const int iat0 = ucell.itia2iat(T0, I0); + // const int start0 = ucell.itiaiw2iwt(T0, I0, 0); + + tau0 = adjs.adjacent_tau[ad0]; + dtau1 = tau0 - tau1; + double distance1 = dtau1.norm() * ucell.lat0; + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); + + dtau2 = tau0 - tau2; + double distance2 = dtau2.norm() * ucell.lat0; + double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); + + if (distance1 < rcut1 && distance2 < rcut2) + { + is_adj = true; + break; + } // dis1, dis2 + } + } + + if (is_adj) + { + info[iat][cb][0] = adjs.box[ad].x; + info[iat][cb][1] = adjs.box[ad].y; + info[iat][cb][2] = adjs.box[ad].z; + info[iat][cb][3] = T2; + info[iat][cb][4] = I2; + ++cb; + } + } // end ad + // GlobalV::ofs_running << " nadj = " << cb << std::endl; + } // end I1 + } // end T1 #ifdef _OPENMP -} + } #endif - ModuleBase::timer::tick("Record_adj","for_2d"); - info_modified=true; - return; + ModuleBase::timer::tick("Record_adj", "for_2d"); + info_modified = true; + return; } - //-------------------------------------------- // This will record the orbitals according to -// grid division (cut along z direction) +// grid division (cut along z direction) //-------------------------------------------- -void Record_adj::for_grid(const UnitCell& ucell, const Grid_Technique >, const std::vector& orb_cutoff) +void Record_adj::for_grid(const UnitCell& ucell, + Grid_Driver& grid_d, + const Grid_Technique& gt, + const std::vector& orb_cutoff) { - ModuleBase::TITLE("Record_adj","for_grid"); - ModuleBase::timer::tick("Record_adj","for_grid"); - - this->na_proc = 0; - this->iat2ca = new int[ucell.nat]; - for(int iat=0; iatna_each = new int[na_proc]; - ModuleBase::GlobalFunc::ZEROS(na_each, na_proc); - this->info = new int**[na_proc]; + ModuleBase::TITLE("Record_adj", "for_grid"); + ModuleBase::timer::tick("Record_adj", "for_grid"); + + this->na_proc = 0; + this->iat2ca = new int[ucell.nat]; + for (int iat = 0; iat < ucell.nat; ++iat) + { + { + if (gt.in_this_processor[iat]) + { + iat2ca[iat] = na_proc; + ++na_proc; + } + else + { + iat2ca[iat] = -1; + } + } + } + + // number of adjacents for each atom. + this->na_each = new int[na_proc]; + ModuleBase::GlobalFunc::ZEROS(na_each, na_proc); + this->info = new int**[na_proc]; #ifdef _OPENMP #pragma omp parallel -{ + { #endif - ModuleBase::Vector3 tau1, tau2, dtau; - ModuleBase::Vector3 tau0, dtau1, dtau2; + ModuleBase::Vector3 tau1, tau2, dtau; + ModuleBase::Vector3 tau0, dtau1, dtau2; #ifdef _OPENMP #pragma omp for schedule(dynamic) #endif - for(int iat=0; iattau[I1]; - //GlobalC::GridD.Find_atom(tau1); - AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(ucell, tau1, T1, I1, &adjs); - for (int ad = 0; ad < adjs.adj_num+1; ad++) - { - const int T2 = adjs.ntype[ad]; - const int I2 = adjs.natom[ad]; - const int iat2 = ucell.itia2iat(T2, I2); - if(gt.in_this_processor[iat2]) - { - //Atom* atom2 = &ucell.atoms[T2]; - tau2 = adjs.adjacent_tau[ad]; - dtau = tau2 - tau1; - double distance = dtau.norm() * ucell.lat0; - double rcut = orb_cutoff[T1] + orb_cutoff[T2]; - - - bool is_adj = false; - if(distance < rcut) { is_adj = true; -} - /* - else if(distance >= rcut) - { - for (int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0) + for (int iat = 0; iat < ucell.nat; ++iat) + { + const int T1 = ucell.iat2it[iat]; + Atom* atom1 = &ucell.atoms[T1]; + const int I1 = ucell.iat2ia[iat]; + { + const int ca = iat2ca[iat]; + // key in this function + if (gt.in_this_processor[iat]) + { + tau1 = atom1->tau[I1]; + // grid_d.Find_atom(tau1); + AdjacentAtomInfo adjs; + grid_d.Find_atom(ucell, tau1, T1, I1, &adjs); + for (int ad = 0; ad < adjs.adj_num + 1; ad++) + { + const int T2 = adjs.ntype[ad]; + const int I2 = adjs.natom[ad]; + const int iat2 = ucell.itia2iat(T2, I2); + if (gt.in_this_processor[iat2]) + { + // Atom* atom2 = &ucell.atoms[T2]; + tau2 = adjs.adjacent_tau[ad]; + dtau = tau2 - tau1; + double distance = dtau.norm() * ucell.lat0; + double rcut = orb_cutoff[T1] + orb_cutoff[T2]; + + bool is_adj = false; + if (distance < rcut) { - const int T0 = GlobalC::GridD.getType(ad0); - const int I0 = GlobalC::GridD.getNatom(ad0); - const int iat0 = ucell.itia2iat(T0, I0); - const int start0 = ucell.itiaiw2iwt(T0, I0, 0); - - tau0 = GlobalC::GridD.getAdjacentTau(ad0); - dtau1 = tau0 - tau1; - dtau2 = tau0 - tau2; - - double distance1 = dtau1.norm() * ucell.lat0; - double distance2 = dtau2.norm() * ucell.lat0; - - double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); - double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); - - if( distance1 < rcut1 && distance2 < rcut2 ) + is_adj = true; + } + /* + else if(distance >= rcut) + { + for (int ad0 = 0; ad0 < grid_d.getAdjacentNum()+1; ++ad0) { - is_adj = true; - break; - } // dis1, dis2 - } - } - */ - - // check the distance - if(is_adj) - { - ++na_each[ca]; - } - }//end judge 2 - }//end ad - }//end judge 1 - }//end I1 - }//end T1 + const int T0 = grid_d.getType(ad0); + const int I0 = grid_d.getNatom(ad0); + const int iat0 = ucell.itia2iat(T0, I0); + const int start0 = ucell.itiaiw2iwt(T0, I0, 0); + + tau0 = grid_d.getAdjacentTau(ad0); + dtau1 = tau0 - tau1; + dtau2 = tau0 - tau2; + + double distance1 = dtau1.norm() * ucell.lat0; + double distance2 = dtau2.norm() * ucell.lat0; + + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); + + if( distance1 < rcut1 && distance2 < rcut2 ) + { + is_adj = true; + break; + } // dis1, dis2 + } + } + */ + + // check the distance + if (is_adj) + { + ++na_each[ca]; + } + } // end judge 2 + } // end ad + } // end judge 1 + } // end I1 + } // end T1 #ifdef _OPENMP #pragma omp for schedule(dynamic) #endif - for(int i=0; i0); - info[i] = new int*[ na_each[i] ]; - for(int j=0; j 0); + info[i] = new int*[na_each[i]]; + for (int j = 0; j < na_each[i]; j++) + { + // (Rx, Ry, Rz, T, I) + info[i][j] = new int[5]; + ModuleBase::GlobalFunc::ZEROS(info[i][j], 5); + } + } #ifdef _OPENMP #pragma omp for schedule(dynamic) #endif - for(int iat=0; iattau[I1]; - //GlobalC::GridD.Find_atom(tau1); - AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(ucell, tau1, T1, I1, &adjs); - - int cb = 0; - for (int ad = 0; ad < adjs.adj_num+1; ad++) - { - const int T2 = adjs.ntype[ad]; - const int I2 = adjs.natom[ad]; - const int iat2 = ucell.itia2iat(T2, I2); - - // key of this function - if(gt.in_this_processor[iat2]) - { - //Atom* atom2 = &ucell.atoms[T2]; - tau2 = adjs.adjacent_tau[ad]; - dtau = tau2 - tau1; - double distance = dtau.norm() * ucell.lat0; - double rcut = orb_cutoff[T1] + orb_cutoff[T2]; - - // check the distance - if(distance < rcut) - { - info[ca][cb][0] = adjs.box[ad].x; - info[ca][cb][1] = adjs.box[ad].y; - info[ca][cb][2] = adjs.box[ad].z; - info[ca][cb][3] = T2; - info[ca][cb][4] = I2; - ++cb; - } - /* - else if(distance >= rcut) + for (int iat = 0; iat < ucell.nat; ++iat) + { + const int T1 = ucell.iat2it[iat]; + Atom* atom1 = &ucell.atoms[T1]; + const int I1 = ucell.iat2ia[iat]; + { + const int ca = iat2ca[iat]; + + // key of this function + if (gt.in_this_processor[iat]) + { + tau1 = atom1->tau[I1]; + // grid_d.Find_atom(tau1); + AdjacentAtomInfo adjs; + grid_d.Find_atom(ucell, tau1, T1, I1, &adjs); + + int cb = 0; + for (int ad = 0; ad < adjs.adj_num + 1; ad++) + { + const int T2 = adjs.ntype[ad]; + const int I2 = adjs.natom[ad]; + const int iat2 = ucell.itia2iat(T2, I2); + + // key of this function + if (gt.in_this_processor[iat2]) { - for (int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0) + // Atom* atom2 = &ucell.atoms[T2]; + tau2 = adjs.adjacent_tau[ad]; + dtau = tau2 - tau1; + double distance = dtau.norm() * ucell.lat0; + double rcut = orb_cutoff[T1] + orb_cutoff[T2]; + + // check the distance + if (distance < rcut) { - const int T0 = GlobalC::GridD.getType(ad0); - const int I0 = GlobalC::GridD.getNatom(ad0); - const int iat0 = ucell.itia2iat(T0, I0); - const int start0 = ucell.itiaiw2iwt(T0, I0, 0); - - tau0 = GlobalC::GridD.getAdjacentTau(ad0); - dtau1 = tau0 - tau1; - dtau2 = tau0 - tau2; - - double distance1 = dtau1.norm() * ucell.lat0; - double distance2 = dtau2.norm() * ucell.lat0; - - double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); - double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); - - if( distance1 < rcut1 && distance2 < rcut2 ) + info[ca][cb][0] = adjs.box[ad].x; + info[ca][cb][1] = adjs.box[ad].y; + info[ca][cb][2] = adjs.box[ad].z; + info[ca][cb][3] = T2; + info[ca][cb][4] = I2; + ++cb; + } + /* + else if(distance >= rcut) + { + for (int ad0 = 0; ad0 < grid_d.getAdjacentNum()+1; ++ad0) { - info[ca][cb][0] = GlobalC::GridD.getBox(ad).x; - info[ca][cb][1] = GlobalC::GridD.getBox(ad).y; - info[ca][cb][2] = GlobalC::GridD.getBox(ad).z; - info[ca][cb][3] = T2; - info[ca][cb][4] = I2; - ++cb; - break; - } // dis1, dis2 - } + const int T0 = grid_d.getType(ad0); + const int I0 = grid_d.getNatom(ad0); + const int iat0 = ucell.itia2iat(T0, I0); + const int start0 = ucell.itiaiw2iwt(T0, I0, 0); + + tau0 = grid_d.getAdjacentTau(ad0); + dtau1 = tau0 - tau1; + dtau2 = tau0 - tau2; + + double distance1 = dtau1.norm() * ucell.lat0; + double distance2 = dtau2.norm() * ucell.lat0; + + double rcut1 = orb_cutoff[T1] + ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut2 = orb_cutoff[T2] + ucell.infoNL.Beta[T0].get_rcut_max(); + + if( distance1 < rcut1 && distance2 < rcut2 ) + { + info[ca][cb][0] = grid_d.getBox(ad).x; + info[ca][cb][1] = grid_d.getBox(ad).y; + info[ca][cb][2] = grid_d.getBox(ad).z; + info[ca][cb][3] = T2; + info[ca][cb][4] = I2; + ++cb; + break; + } // dis1, dis2 + } + } + */ } - */ - } - }// end ad - - assert(cb == na_each[ca]); - } - } - } + } // end ad + + assert(cb == na_each[ca]); + } + } + } #ifdef _OPENMP -} + } #endif - ModuleBase::timer::tick("Record_adj","for_grid"); - info_modified=true; -// std::cout << " after for_grid" << std::endl; - return; + ModuleBase::timer::tick("Record_adj", "for_grid"); + info_modified = true; + // std::cout << " after for_grid" << std::endl; + return; } - diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h index 06cf238a81..8a80b63cee 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h @@ -1,61 +1,66 @@ #ifndef RECORD_ADJ_H #define RECORD_ADJ_H -#include "module_hamilt_lcao/module_gint/grid_technique.h" #include "module_basis/module_ao/parallel_orbitals.h" - +#include "module_hamilt_lcao/module_gint/grid_technique.h" //--------------------------------------------------- // FUNCTION: record the adjacent atoms for each atom //--------------------------------------------------- class Record_adj { - private: - bool info_modified=false; - public: - - Record_adj(); - ~Record_adj(); - - //-------------------------------------------- - // This will record the orbitals according to - // HPSEPS's 2D block division. - //-------------------------------------------- - void for_2d(const UnitCell& ucell, Parallel_Orbitals &pv, bool gamma_only, const std::vector& orb_cutoff); - - //-------------------------------------------- - // This will record the orbitals according to - // grid division (cut along z direction) - //-------------------------------------------- - void for_grid(const UnitCell& ucell, const Grid_Technique >, const std::vector& orb_cutoff); - - void delete_grid(); - - int na_proc; - int* na_each; - - //-------------------------------------------- - // record sparse atom index in for_grid(const Grid_Technique >); - // Map iat(dense atom index) to sparse atom index - // Mainly removing the index dependency for OpenMP parallel loop - // - // Meaning: - // 1. if iat2ca[iat] > 0, it contains the sparse atom index - // 2. if iat2ca[iat] < 0, the sparse atom index of iat does not exist - // - // Usage: - // 1. iat2ca[iat] > 0 ? na_each[iat2ca[iat]] : 0 - // 2. iat2ca[iat] > 0 ? info[iat2ca[iat]] : nullptr - //-------------------------------------------- - int* iat2ca; - - //------------------------------------------------ - // info will identify each atom in each unitcell. - //------------------------------------------------ - int*** info; - - private: + private: + bool info_modified = false; + + public: + Record_adj(); + ~Record_adj(); + + //-------------------------------------------- + // This will record the orbitals according to + // HPSEPS's 2D block division. + //-------------------------------------------- + void for_2d(const UnitCell& ucell, + Grid_Driver& grid_d, + Parallel_Orbitals& pv, + bool gamma_only, + const std::vector& orb_cutoff); + + //-------------------------------------------- + // This will record the orbitals according to + // grid division (cut along z direction) + //-------------------------------------------- + void for_grid(const UnitCell& ucell, + Grid_Driver& grid_d, + const Grid_Technique& gt, + const std::vector& orb_cutoff); + + void delete_grid(); + + int na_proc; + int* na_each; + + //-------------------------------------------- + // record sparse atom index in for_grid(const Grid_Technique >); + // Map iat(dense atom index) to sparse atom index + // Mainly removing the index dependency for OpenMP parallel loop + // + // Meaning: + // 1. if iat2ca[iat] > 0, it contains the sparse atom index + // 2. if iat2ca[iat] < 0, the sparse atom index of iat does not exist + // + // Usage: + // 1. iat2ca[iat] > 0 ? na_each[iat2ca[iat]] : 0 + // 2. iat2ca[iat] > 0 ? info[iat2ca[iat]] : nullptr + //-------------------------------------------- + int* iat2ca; + + //------------------------------------------------ + // info will identify each atom in each unitcell. + //------------------------------------------------ + int*** info; + private: }; #endif