From d4ed2552b7c7dfdc84226e0b807f68579717a581 Mon Sep 17 00:00:00 2001 From: YuLiu98 Date: Mon, 2 Dec 2024 16:04:34 +0800 Subject: [PATCH] Refactor: remove GlobalC::Grid_D in hamilt_lcao --- .../module_cell/module_neighbor/sltk_grid.h | 184 ++++++----- .../module_neighbor/sltk_grid_driver.h | 80 +++-- source/module_esolver/cell_dependency.h | 20 ++ source/module_esolver/esolver_gets.cpp | 10 +- source/module_esolver/lcao_before_scf.cpp | 3 +- source/module_esolver/lcao_others.cpp | 3 +- .../hamilt_lcaodft/hamilt_lcao.cpp | 130 ++++---- .../hamilt_lcaodft/hamilt_lcao.h | 56 ++-- .../hamilt_lcaodft/record_adj.cpp | 294 +++++++++--------- .../hamilt_lcaodft/record_adj.h | 14 +- 10 files changed, 419 insertions(+), 375 deletions(-) create mode 100644 source/module_esolver/cell_dependency.h diff --git a/source/module_cell/module_neighbor/sltk_grid.h b/source/module_cell/module_neighbor/sltk_grid.h index 18f3470d6c..b3eb1937c7 100644 --- a/source/module_cell/module_neighbor/sltk_grid.h +++ b/source/module_cell/module_neighbor/sltk_grid.h @@ -53,70 +53,69 @@ class Grid // 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(const Atom_input &input); - - - //Static data - static const double TOLERATE_ERROR; - static const std::hash INT_HASHER; - static const char* const ERROR[3]; - - //Data - int natom;// Total atoms. - bool pbc; // periodic boundary condition - bool expand_flag; - double sradius;// searching radius - double d_minX;// origin of all cells - double d_minY; - double d_minZ; - int dx; - int dy; - int dz; - int layer; - double cell_x_length; - double cell_y_length; - double cell_z_length; - CellSet ***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->dx;i++) - { - for (int j = 0;j < this->dy;j++) - { - delete[] this->Cell[i][j]; - } - } - - for (int i = 0;i < this->dx;i++) - { - delete[] this->Cell[i]; - } - - delete[] this->Cell; - this->init_cell_flag = false; - } - } - - double grid_length[3]; - double vec1[3]; - double vec2[3]; - double vec3[3]; - double lat_now; - bool init_cell_flag; - //LiuXh add 2019-07-15 - const double& getD_minX(void) const {return d_minX;} + 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(const Atom_input& input); + + // Static data + static const double TOLERATE_ERROR; + static const std::hash INT_HASHER; + static const char* const ERROR[3]; + + // Data + int natom; // Total atoms. + bool pbc; // periodic boundary condition + bool expand_flag; + double sradius; // searching radius + double d_minX; // origin of all cells + double d_minY; + double d_minZ; + int dx; + int dy; + int dz; + int layer; + double cell_x_length; + double cell_y_length; + double cell_z_length; + CellSet*** 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->dx; i++) + { + for (int j = 0; j < this->dy; j++) + { + delete[] this->Cell[i][j]; + } + } + + for (int i = 0; i < this->dx; i++) + { + delete[] this->Cell[i]; + } + + delete[] this->Cell; + this->init_cell_flag = false; + } + } + + double grid_length[3]; + double vec1[3]; + double vec2[3]; + double vec3[3]; + double lat_now; + bool init_cell_flag = false; + // LiuXh add 2019-07-15 + const double& getD_minX(void) const + { + return d_minX; + } const double& getD_minY(void) const {return d_minY;} const double& getD_minZ(void) const {return d_minZ;} @@ -128,38 +127,31 @@ class Grid protected: AtomLink* getHashCode(const UnitCell &ucell, const FAtom &atom)const; // Peize Lin delete const 2018-07-14 // AtomLink* const getHashCode(const FAtom &atom)const; - AtomLink* atomlink; - AtomLink* cordon_p;// Warning! A guard! Don't delete it! - -private: - - const int test_grid; -//========================================================== -// MEMBER FUNCTIONS : -// Three Main Steps: -// NAME : setMemberVariables (read in datas from Atom_input, -// init cells.) -// NAME : setAtomLinkArray( set the AtomLinkArray twice, -// first use Build_Hash,second use Fold_Hash) -// 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 setAtomLinkArray( - const UnitCell &ucell, - const Atom_input &input); - - void setBoundaryAdjacent( - std::ofstream &ofs_in, - const Atom_input &input); - -//========================================================== -// -//========================================================== + AtomLink* atomlink = nullptr; + AtomLink* cordon_p = nullptr; // Warning! A guard! Don't delete it! + + private: + const int test_grid; + //========================================================== + // MEMBER FUNCTIONS : + // Three Main Steps: + // NAME : setMemberVariables (read in datas from Atom_input, + // init cells.) + // NAME : setAtomLinkArray( set the AtomLinkArray twice, + // first use Build_Hash,second use Fold_Hash) + // 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 setAtomLinkArray(const UnitCell& ucell, const Atom_input& input); + + void setBoundaryAdjacent(std::ofstream& ofs_in, const Atom_input& input); + + //========================================================== + // + //========================================================== AtomLink* Build_Cache(const UnitCell &ucell, const Atom_input &input); // Peize Lin delete const and throw(std::out_of_range, std::logic_error) 2018-07-14 diff --git a/source/module_cell/module_neighbor/sltk_grid_driver.h b/source/module_cell/module_neighbor/sltk_grid_driver.h index 61216665bb..8ea7f0ad52 100644 --- a/source/module_cell/module_neighbor/sltk_grid_driver.h +++ b/source/module_cell/module_neighbor/sltk_grid_driver.h @@ -39,43 +39,55 @@ class Grid_Driver : public Grid // 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() : 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(void)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(void) 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: 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 7cd0258621..deca9ad97c 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 { @@ -101,14 +101,15 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) PARAM.inp.test_atom_input); Record_adj RA; - RA.for_2d(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + RA.for_2d(GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); if (this->p_hamilt == nullptr) { if (PARAM.inp.nspin == 4) { this->p_hamilt - = new hamilt::HamiltLCAO, std::complex>(&this->pv, + = new hamilt::HamiltLCAO, std::complex>(GlobalC::GridD, + &this->pv, this->kv, *(two_center_bundle_.overlap_orb), orb_.cutoffs()); @@ -117,7 +118,8 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep) } else { - this->p_hamilt = new hamilt::HamiltLCAO, double>(&this->pv, + this->p_hamilt = new hamilt::HamiltLCAO, double>(GlobalC::GridD, + &this->pv, this->kv, *(two_center_bundle_.overlap_orb), orb_.cutoffs()); diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 3428efc21d..b5f0bb62a2 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -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(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation @@ -188,6 +188,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) this->p_hamilt = new hamilt::HamiltLCAO( PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), + GlobalC::GridD, &this->pv, this->pelec->pot, this->kv, diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 43e0f7a56b..edfeae8084 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -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(this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); + this->RA.for_2d(GlobalC::GridD, this->pv, PARAM.globalv.gamma_only_local, orb_.cutoffs()); // 2. density matrix extrapolation @@ -189,6 +189,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) this->p_hamilt = new hamilt::HamiltLCAO( PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), + 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 dfc9d0a55a..6b1e5af297 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -43,7 +43,11 @@ namespace hamilt { template -HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& kv_in, const TwoCenterIntegrator& intor_overlap_orb, const std::vector& orb_cutoff) +HamiltLCAO::HamiltLCAO(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"; @@ -60,24 +64,26 @@ HamiltLCAO::HamiltLCAO(const Parallel_Orbitals* paraV, const K_Vectors& this->sR, &GlobalC::ucell, orb_cutoff, - &GlobalC::GridD, + &grid_d, &intor_overlap_orb); } template HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, - Gint_k* GK_in, - 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, + 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; @@ -135,7 +141,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->sR, &GlobalC::ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb.get()); // kinetic term () @@ -146,7 +152,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, &GlobalC::ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.kinetic_orb.get()); this->getOperator()->add(ekinetic); } @@ -160,7 +166,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, &GlobalC::ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_beta.get()); this->getOperator()->add(nonlocal); } @@ -176,15 +182,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 - &GlobalC::ucell, - orb.cutoffs(), - &GlobalC::GridD, - PARAM.inp.nspin - ); + this->hsk, + this->kv->kvec_d, + pot_in, + this->hR, // no explicit call yet + &GlobalC::ucell, + orb.cutoffs(), + &grid_d, + PARAM.inp.nspin); this->getOperator()->add(veff); } } @@ -196,7 +201,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, // no explicit call yet &GlobalC::ucell, - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_alpha.get(), &orb, this->kv->get_nks(), @@ -222,7 +227,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, GlobalC::ucell, - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), &GlobalC::dftu); @@ -245,14 +250,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, - &GlobalC::ucell, - orb.cutoffs(), - &GlobalC::GridD, - PARAM.inp.nspin); + this->hsk, + kv->kvec_d, + pot_in, + this->hR, + &GlobalC::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); @@ -267,7 +272,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->sR, &GlobalC::ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb.get()); if (this->getOperator() == nullptr) { @@ -287,7 +292,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, &GlobalC::ucell, orb.cutoffs(), - &GlobalC::GridD, + &grid_d, two_center_bundle.kinetic_orb.get()); this->getOperator()->add(ekinetic); } @@ -301,7 +306,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->hR, &GlobalC::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 @@ -322,7 +327,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, hR, &GlobalC::ucell, - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_alpha.get(), &orb, this->kv->get_nks(), @@ -335,22 +340,18 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, { if(!TD_Velocity::init_vecpot_file) { elecstate::H_TDDFT_pw::update_At(); } - Operator* td_ekinetic = new TDEkinetic>(this->hsk, - this->hR, - kv, - &GlobalC::ucell, - orb.cutoffs(), - &GlobalC::GridD, - two_center_bundle.overlap_orb.get()); - this->getOperator()->add(td_ekinetic); - - Operator* td_nonlocal = new TDNonlocal>(this->hsk, - this->kv->kvec_d, - this->hR, - &GlobalC::ucell, - orb, - &GlobalC::GridD); - this->getOperator()->add(td_nonlocal); +Operator* td_ekinetic = new TDEkinetic>(this->hsk, + this->hR, + kv, + &GlobalC::ucell, + orb.cutoffs(), + &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, &GlobalC::ucell, orb, &grid_d); +this->getOperator()->add(td_nonlocal); } if (PARAM.inp.dft_plus_u) { @@ -368,7 +369,7 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, this->kv->kvec_d, this->hR, GlobalC::ucell, - &GlobalC::GridD, + &grid_d, two_center_bundle.overlap_orb_onsite.get(), orb.cutoffs(), &GlobalC::dftu); @@ -377,16 +378,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, - GlobalC::ucell, - &GlobalC::GridD, - two_center_bundle.overlap_orb_onsite.get(), - orb.cutoffs() - ); + Operator* sc_lambda = new DeltaSpin>(this->hsk, + kv->kvec_d, + this->hR, + GlobalC::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); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h index 0cd597b5b0..42936b7be9 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" @@ -29,35 +31,41 @@ class HamiltLCAO : public Hamilt */ using TAC = std::pair>; HamiltLCAO(Gint_Gamma* GG_in, - Gint_k* GK_in, - 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, + 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 Parallel_Orbitals* paraV, const K_Vectors& kv_in, const TwoCenterIntegrator& intor_overlap_orb, const std::vector& orb_cutoff); - - ~HamiltLCAO() - { - if (this->ops != nullptr) - { - delete this->ops; - } - delete this->hR; - delete this->sR; - delete this->hsk; - } + HamiltLCAO(Grid_Driver& grid_d, + const Parallel_Orbitals* paraV, + const K_Vectors& kv_in, + const TwoCenterIntegrator& intor_overlap_orb, + const std::vector& orb_cutoff); + + ~HamiltLCAO() + { + if (this->ops != nullptr) + { + delete this->ops; + } + delete this->hR; + delete this->sR; + delete this->hsk; + } /// get pointer of Operator ops Operator*& getOperator(); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp index 43569e2e24..bdbf5a1c21 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.cpp @@ -38,7 +38,10 @@ void Record_adj::delete_grid() // If multi-k, calculate nnr at the same time. // be called only once in an ion-step. //-------------------------------------------- -void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vector& orb_cutoff) +void Record_adj::for_2d(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"); @@ -73,82 +76,89 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto for (int I1 = 0; I1 < atom1->na; ++I1) { tau1 = atom1->tau[I1]; - //GlobalC::GridD.Find_atom( tau1 ); - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1 ,T1, I1); - const int start1 = GlobalC::ucell.itiaiw2iwt(T1, I1, 0); + // grid_d.Find_atom( tau1 ); + grid_d.Find_atom(GlobalC::ucell, tau1, T1, I1); + const int start1 = GlobalC::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 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); - tau2 = GlobalC::GridD.getAdjacentTau(ad); - dtau = tau2 - tau1; - double distance = dtau.norm() * GlobalC::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 = GlobalC::ucell.itia2iat(T0, I0); - //const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); - - tau0 = GlobalC::GridD.getAdjacentTau(ad0); - dtau1 = tau0 - tau1; - double distance1 = dtau1.norm() * GlobalC::ucell.lat0; - double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - - dtau2 = tau0 - tau2; - double distance2 = dtau2.norm() * GlobalC::ucell.lat0; - double rcut2 = orb_cutoff[T2] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - - if( distance1 < rcut1 && distance2 < rcut2 ) - { - is_adj = true; - break; - } // dis1, dis2 - } - } +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 = GlobalC::ucell.itiaiw2iwt(T2, I2, 0); + tau2 = grid_d.getAdjacentTau(ad); + dtau = tau2 - tau1; + double distance = dtau.norm() * GlobalC::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 < grid_d.getAdjacentNum() + 1; ++ad0) + { + const int T0 = grid_d.getType(ad0); + // const int I0 = grid_d.getNatom(ad0); + // const int iat0 = GlobalC::ucell.itia2iat(T0, I0); + // const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); + + tau0 = grid_d.getAdjacentTau(ad0); + dtau1 = tau0 - tau1; + double distance1 = dtau1.norm() * GlobalC::ucell.lat0; + double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + + dtau2 = tau0 - tau2; + double distance2 = dtau2.norm() * GlobalC::ucell.lat0; + double rcut2 = orb_cutoff[T2] + GlobalC::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) + 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 < GlobalC::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) { - for(int ii=0; iinw * 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; jjtau[I1]; - //GlobalC::GridD.Find_atom( tau1 ); - AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1 ,T1, I1, &adjs); + // grid_d.Find_atom( tau1 ); + AdjacentAtomInfo adjs; + grid_d.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs); - // (2) search among all adjacent atoms. - int cb = 0; + // (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]; @@ -220,12 +230,12 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto 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 = GlobalC::ucell.itia2iat(T0, I0); - //const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); + // const int I0 = grid_d.getNatom(ad0); + // const int iat0 = GlobalC::ucell.itia2iat(T0, I0); + // const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); - tau0 = adjs.adjacent_tau[ad0]; - dtau1 = tau0 - tau1; + tau0 = adjs.adjacent_tau[ad0]; + dtau1 = tau0 - tau1; double distance1 = dtau1.norm() * GlobalC::ucell.lat0; double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); @@ -267,7 +277,7 @@ void Record_adj::for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vecto // This will record the orbitals according to // grid division (cut along z direction) //-------------------------------------------- -void Record_adj::for_grid(const Grid_Technique >, const std::vector& orb_cutoff) +void Record_adj::for_grid(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"); @@ -312,11 +322,11 @@ void Record_adj::for_grid(const Grid_Technique >, const std::vector& o if(gt.in_this_processor[iat]) { tau1 = atom1->tau[I1]; - //GlobalC::GridD.Find_atom(tau1); - AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs); - for (int ad = 0; ad < adjs.adj_num+1; ad++) - { + // grid_d.Find_atom(tau1); + AdjacentAtomInfo adjs; + grid_d.Find_atom(GlobalC::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 = GlobalC::ucell.itia2iat(T2, I2); @@ -332,41 +342,41 @@ void Record_adj::for_grid(const Grid_Technique >, const std::vector& o bool is_adj = false; if(distance < rcut) { is_adj = true; } - /* - 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 = GlobalC::ucell.itia2iat(T0, I0); - const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); - - tau0 = GlobalC::GridD.getAdjacentTau(ad0); - dtau1 = tau0 - tau1; - dtau2 = tau0 - tau2; - - double distance1 = dtau1.norm() * GlobalC::ucell.lat0; - double distance2 = dtau2.norm() * GlobalC::ucell.lat0; - - double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - double rcut2 = orb_cutoff[T2] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); - - if( distance1 < rcut1 && distance2 < rcut2 ) - { - is_adj = true; - break; - } // dis1, dis2 - } - } - */ +/* +else if(distance >= rcut) +{ + for (int ad0 = 0; ad0 < grid_d.getAdjacentNum()+1; ++ad0) + { + const int T0 = grid_d.getType(ad0); + const int I0 = grid_d.getNatom(ad0); + const int iat0 = GlobalC::ucell.itia2iat(T0, I0); + const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); + + tau0 = grid_d.getAdjacentTau(ad0); + dtau1 = tau0 - tau1; + dtau2 = tau0 - tau2; + + double distance1 = dtau1.norm() * GlobalC::ucell.lat0; + double distance2 = dtau2.norm() * GlobalC::ucell.lat0; + + double rcut1 = orb_cutoff[T1] + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max(); + double rcut2 = orb_cutoff[T2] + GlobalC::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 +// check the distance +if (is_adj) +{ + ++na_each[ca]; +} + }//end judge 2 }//end ad }//end judge 1 }//end I1 @@ -402,12 +412,12 @@ void Record_adj::for_grid(const Grid_Technique >, const std::vector& o if(gt.in_this_processor[iat]) { tau1 = atom1->tau[I1]; - //GlobalC::GridD.Find_atom(tau1); - AdjacentAtomInfo adjs; - GlobalC::GridD.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs); + // grid_d.Find_atom(tau1); + AdjacentAtomInfo adjs; + grid_d.Find_atom(GlobalC::ucell, tau1, T1, I1, &adjs); - int cb = 0; - for (int ad = 0; ad < adjs.adj_num+1; ad++) + 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]; @@ -432,17 +442,17 @@ void Record_adj::for_grid(const Grid_Technique >, const std::vector& o info[ca][cb][4] = I2; ++cb; } - /* + /* else if(distance >= rcut) { - for (int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0) + for (int ad0 = 0; ad0 < grid_d.getAdjacentNum()+1; ++ad0) { - const int T0 = GlobalC::GridD.getType(ad0); - const int I0 = GlobalC::GridD.getNatom(ad0); + const int T0 = grid_d.getType(ad0); + const int I0 = grid_d.getNatom(ad0); const int iat0 = GlobalC::ucell.itia2iat(T0, I0); const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0); - tau0 = GlobalC::GridD.getAdjacentTau(ad0); + tau0 = grid_d.getAdjacentTau(ad0); dtau1 = tau0 - tau1; dtau2 = tau0 - tau2; @@ -454,19 +464,19 @@ void Record_adj::for_grid(const Grid_Technique >, const std::vector& o if( distance1 < rcut1 && distance2 < rcut2 ) { - 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 - } + 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 + */ + } + }// end ad assert(cb == na_each[ca]); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h index d2f214632f..1cf29c8a3a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/record_adj.h @@ -21,15 +21,15 @@ class Record_adj // This will record the orbitals according to // HPSEPS's 2D block division. //-------------------------------------------- - void for_2d(Parallel_Orbitals &pv, bool gamma_only, const std::vector& orb_cutoff); + void for_2d(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 Grid_Technique >, const std::vector& orb_cutoff); + //-------------------------------------------- + // This will record the orbitals according to + // grid division (cut along z direction) + //-------------------------------------------- + void for_grid(Grid_Driver& grid_d, const Grid_Technique& gt, const std::vector& orb_cutoff); - void delete_grid(); + void delete_grid(); int na_proc; int* na_each;