diff --git a/source/module_hamilt_lcao/module_gint/gint.cpp b/source/module_hamilt_lcao/module_gint/gint.cpp index 6ee4c238d1..dee289c12d 100644 --- a/source/module_hamilt_lcao/module_gint/gint.cpp +++ b/source/module_hamilt_lcao/module_gint/gint.cpp @@ -23,7 +23,7 @@ Gint::~Gint() { delete this->hRGint; - delete this->hRGintCd; + delete this->hR_tmp; // in gamma_only case, DMRGint.size()=0, // in multi-k case, DMRGint.size()=nspin for (int is = 0; is < this->DMRGint.size(); is++) { @@ -33,7 +33,7 @@ Gint::~Gint() { delete this->hRGint_tmp[is]; } #ifdef __MPI - delete this->DMRGint_full; + delete this->DM2D_tmp; #endif } @@ -155,11 +155,9 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, const Grid_Driver* gd, cons this->hRGint = new hamilt::HContainer(ucell_in.nat); } else { npol = 2; - if (this->hRGintCd != nullptr) { - delete this->hRGintCd; + if (this->hR_tmp != nullptr) { + delete this->hR_tmp; } - this->hRGintCd - = new hamilt::HContainer>(ucell_in.nat); for (int is = 0; is < nspin; is++) { if (this->DMRGint[is] != nullptr) { delete this->DMRGint[is]; @@ -171,10 +169,9 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, const Grid_Driver* gd, cons this->hRGint_tmp[is] = new hamilt::HContainer(ucell_in.nat); } #ifdef __MPI - if (this->DMRGint_full != nullptr) { - delete this->DMRGint_full; + if (this->DM2D_tmp != nullptr) { + delete this->DM2D_tmp; } - this->DMRGint_full = new hamilt::HContainer(ucell_in.nat); #endif } @@ -197,8 +194,10 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, const Grid_Driver* gd, cons this->DMRGint[0]->get_memory_size() * this->DMRGint.size()); } else { - this->hRGintCd->insert_ijrs(this->gridt->get_ijr_info(), ucell_in, npol); - this->hRGintCd->allocate(nullptr, true); + // this->hRGintCd->insert_ijrs(this->gridt->get_ijr_info(), ucell_in, npol); + // this->hRGintCd->allocate(nullptr, true); + // ModuleBase::Memory::record("Gint::hRGintCd", + // this->hRGintCd->get_memory_size()); for(int is = 0; is < nspin; is++) { this->hRGint_tmp[is]->insert_ijrs(this->gridt->get_ijr_info(), ucell_in); this->DMRGint[is]->insert_ijrs(this->gridt->get_ijr_info(), ucell_in); @@ -208,14 +207,7 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, const Grid_Driver* gd, cons ModuleBase::Memory::record("Gint::hRGint_tmp", this->hRGint_tmp[0]->get_memory_size()*nspin); ModuleBase::Memory::record("Gint::DMRGint", - this->DMRGint[0]->get_memory_size() - * this->DMRGint.size()*nspin); -#ifdef __MPI - this->DMRGint_full->insert_ijrs(this->gridt->get_ijr_info(), ucell_in, npol); - this->DMRGint_full->allocate(nullptr, true); - ModuleBase::Memory::record("Gint::DMRGint_full", - this->DMRGint_full->get_memory_size()); -#endif + this->DMRGint[0]->get_memory_size()*nspin); } } @@ -231,10 +223,8 @@ void Gint::reset_DMRGint(const int& nspin) { for (auto& d : this->DMRGint) { d->allocate(nullptr, false); } #ifdef __MPI - delete this->DMRGint_full; - this->DMRGint_full = new hamilt::HContainer(*this->hRGint); - this->DMRGint_full->allocate(nullptr, false); -#endif + delete this->DM2D_tmp; +#endif } } } @@ -262,37 +252,51 @@ void Gint::transfer_DM2DtoGrid(std::vector*> DM2D) { } else // NSPIN=4 case { #ifdef __MPI - hamilt::transferParallels2Serials(*DM2D[0], this->DMRGint_full); -#else - this->DMRGint_full = DM2D[0]; -#endif - std::vector tmp_pointer(4, nullptr); - for (int iap = 0; iap < this->DMRGint_full->size_atom_pairs(); ++iap) { - auto& ap = this->DMRGint_full->get_atom_pair(iap); - int iat1 = ap.get_atom_i(); - int iat2 = ap.get_atom_j(); - for (int ir = 0; ir < ap.get_R_size(); ++ir) { - const ModuleBase::Vector3 r_index = ap.get_R_index(ir); - for (int is = 0; is < 4; is++) { - tmp_pointer[is] = this->DMRGint[is] - ->find_matrix(iat1, iat2, r_index) - ->get_pointer(); - } - double* data_full = ap.get_pointer(ir); - for (int irow = 0; irow < ap.get_row_size(); irow += 2) { - for (int icol = 0; icol < ap.get_col_size(); icol += 2) { - *(tmp_pointer[0])++ = data_full[icol]; - *(tmp_pointer[1])++ = data_full[icol + 1]; - } - data_full += ap.get_col_size(); - for (int icol = 0; icol < ap.get_col_size(); icol += 2) { - *(tmp_pointer[2])++ = data_full[icol]; - *(tmp_pointer[3])++ = data_full[icol + 1]; + // is=0:↑↑, 1:↑↓, 2:↓↑, 3:↓↓ + const int row_set[4] = {0, 0, 1, 1}; + const int col_set[4] = {0, 1, 0, 1}; + + int mg = DM2D[0]->get_paraV()->get_global_row_size()/2; + int ng = DM2D[0]->get_paraV()->get_global_col_size()/2; + int nb = DM2D[0]->get_paraV()->get_block_size()/2; + int blacs_ctxt = DM2D[0]->get_paraV()->blacs_ctxt; + int *iat2iwt = new int[ucell->nat]; + for (int iat = 0; iat < ucell->nat; iat++) { + iat2iwt[iat] = ucell->get_iat2iwt()[iat]/2; + } + Parallel_Orbitals *pv = new Parallel_Orbitals(); + pv->set(mg, ng, nb, blacs_ctxt); + pv->set_atomic_trace(iat2iwt, ucell->nat, mg); + auto ijr_info = DM2D[0]->get_ijr_info(); + this-> DM2D_tmp = new hamilt::HContainer(pv, nullptr, &ijr_info); + ModuleBase::Memory::record("Gint::DM2D_tmp", this->DM2D_tmp->get_memory_size()); + for (int is = 0; is < 4; is++){ + this->DM2D_tmp->set_zero(); + for (int iap = 0; iap < DM2D[0]->size_atom_pairs(); ++iap) { + auto& ap = DM2D[0]->get_atom_pair(iap); + int iat1 = ap.get_atom_i(); + int iat2 = ap.get_atom_j(); + for (int ir = 0; ir < ap.get_R_size(); ++ir) { + const ModuleBase::Vector3 r_index = ap.get_R_index(ir); + double* matrix_out = this -> DM2D_tmp -> find_matrix(iat1, iat2, r_index)->get_pointer(); + double* matrix_in = ap.get_pointer(ir); + for (int irow = 0; irow < ap.get_row_size()/2; irow ++) { + for (int icol = 0; icol < ap.get_col_size()/2; icol++){ + int index_i = irow* ap.get_col_size()/2 + icol; + int index_j = (irow*2+row_set[is]) * ap.get_col_size() + icol*2+col_set[is]; + matrix_out[index_i] = matrix_in[index_j]; + } } - data_full += ap.get_col_size(); } } + hamilt::transferParallels2Serials( *(this->DM2D_tmp), this->DMRGint[is]); } + // delete iat2iwt; + // delete pv; + // delete this-> DM2D_tmp; +#else + //this->DMRGint_full = DM2D[0]; +#endif } ModuleBase::timer::tick("Gint", "transfer_DMR"); } \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_gint/gint.h b/source/module_hamilt_lcao/module_gint/gint.h index 22d75e650b..d8090df191 100644 --- a/source/module_hamilt_lcao/module_gint/gint.h +++ b/source/module_hamilt_lcao/module_gint/gint.h @@ -258,13 +258,13 @@ class Gint { std::vector*> hRGint_tmp; //! stores Hamiltonian in sparse format - hamilt::HContainer>* hRGintCd = nullptr; + hamilt::HContainer>* hR_tmp = nullptr; //! stores DMR in sparse format std::vector*> DMRGint; //! tmp tools used in transfer_DM2DtoGrid - hamilt::HContainer* DMRGint_full = nullptr; + hamilt::HContainer* DM2D_tmp = nullptr; std::vector> pvdpRx_reduced; std::vector> pvdpRy_reduced; diff --git a/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp b/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp index 9a6cca47f7..9e80ea22c1 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp @@ -78,84 +78,114 @@ void Gint_k::transfer_pvpR(hamilt::HContainer>* hR, ModuleBase::TITLE("Gint_k", "transfer_pvpR"); ModuleBase::timer::tick("Gint_k", "transfer_pvpR"); - this->hRGintCd->set_zero(); - - for (int iap = 0; iap < this->hRGintCd->size_atom_pairs(); iap++) - { - auto* ap = &this->hRGintCd->get_atom_pair(iap); - const int iat1 = ap->get_atom_i(); - const int iat2 = ap->get_atom_j(); - if (iat1 <= iat2) - { - hamilt::AtomPair>* upper_ap = ap; - hamilt::AtomPair>* lower_ap = this->hRGintCd->find_pair(iat2, iat1); - const hamilt::AtomPair* ap_nspin_0 = this->hRGint_tmp[0]->find_pair(iat1, iat2); - const hamilt::AtomPair* ap_nspin_3 = this->hRGint_tmp[3]->find_pair(iat1, iat2); - for (int ir = 0; ir < upper_ap->get_R_size(); ir++) - { - const auto R_index = upper_ap->get_R_index(ir); - auto upper_mat = upper_ap->find_matrix(R_index); - auto mat_nspin_0 = ap_nspin_0->find_matrix(R_index); - auto mat_nspin_3 = ap_nspin_3->find_matrix(R_index); + int mg = hR->get_paraV()->get_global_row_size()/2; + int ng = hR->get_paraV()->get_global_col_size()/2; + int nb = hR->get_paraV()->get_block_size()/2; + hR->set_zero(); +#ifdef __MPI + int blacs_ctxt = hR->get_paraV()->blacs_ctxt; + int *iat2iwt = new int[ucell_in->nat]; + for (int iat = 0; iat < ucell_in->nat; iat++) { + iat2iwt[iat] = ucell_in->get_iat2iwt()[iat]/2; + } + Parallel_Orbitals *pv = new Parallel_Orbitals(); + pv->set(mg, ng, nb, blacs_ctxt); + pv->set_atomic_trace(iat2iwt, ucell_in->nat, mg); + auto ijr_info = hR->get_ijr_info(); - // The row size and the col size of upper_matrix is double that of matrix_nspin_0 - for (int irow = 0; irow < mat_nspin_0->get_row_size(); ++irow) - { - for (int icol = 0; icol < mat_nspin_0->get_col_size(); ++icol) - { - upper_mat->get_value(2*irow, 2*icol) = mat_nspin_0->get_value(irow, icol) + mat_nspin_3->get_value(irow, icol); - upper_mat->get_value(2*irow+1, 2*icol+1) = mat_nspin_0->get_value(irow, icol) - mat_nspin_3->get_value(irow, icol); - } - } + this->hR_tmp = new hamilt::HContainer>(pv, nullptr, &ijr_info); + ModuleBase::Memory::record("Gint::hRGintCd", this->hR_tmp->get_memory_size()); - if (PARAM.globalv.domag) - { - const hamilt::AtomPair* ap_nspin_1 = this->hRGint_tmp[1]->find_pair(iat1, iat2); - const hamilt::AtomPair* ap_nspin_2 = this->hRGint_tmp[2]->find_pair(iat1, iat2); - const auto mat_nspin_1 = ap_nspin_1->find_matrix(R_index); - const auto mat_nspin_2 = ap_nspin_2->find_matrix(R_index); - for (int irow = 0; irow < mat_nspin_1->get_row_size(); ++irow) + //0,3;1,2;1,2;0,3 + std::vector first = {0, 1, 1, 0}; + std::vector second= {3, 2, 2, 3}; + std::vector row_set = {0, 0, 1, 1}; + std::vector col_set = {0, 1, 0, 1}; + std::vector clx_i = {1, 0, 0, -1}; + std::vector clx_j = {0, 1, -1, 0}; + for (int is = 0; is < 4; is++){ + hamilt::HContainer>* hRGint_tmpCd = new hamilt::HContainer>(this->ucell->nat); + hRGint_tmpCd->insert_ijrs(this->gridt->get_ijr_info(), *(this->ucell)); + hRGint_tmpCd->allocate(nullptr, true); + hRGint_tmpCd->set_zero(); + for (int iap = 0; iap < hRGint_tmpCd->size_atom_pairs(); iap++) + { + //std::cout<<"iap: "<get_atom_pair(iap); + const int iat1 = ap->get_atom_i(); + const int iat2 = ap->get_atom_j(); + if (iat1 <= iat2) + { + hamilt::AtomPair>* upper_ap = ap; + hamilt::AtomPair>* lower_ap = hRGint_tmpCd->find_pair(iat2, iat1); + const hamilt::AtomPair* ap_nspin1 = this->hRGint_tmp[first[is]] ->find_pair(iat1, iat2); + const hamilt::AtomPair* ap_nspin2 = this->hRGint_tmp[second[is]] ->find_pair(iat1, iat2); + for (int ir = 0; ir < upper_ap->get_R_size(); ir++) + { + //std::cout<<"ir"<get_R_index(ir); + auto upper_mat = upper_ap->find_matrix(R_index); + auto mat_nspin1 = ap_nspin1->find_matrix(R_index); + auto mat_nspin2 = ap_nspin2->find_matrix(R_index); + // The row size and the col size of upper_matrix is double that of matrix_nspin_0 + for (int irow = 0; irow < mat_nspin1->get_row_size(); ++irow) { - for (int icol = 0; icol < mat_nspin_1->get_col_size(); ++icol) + for (int icol = 0; icol < mat_nspin1->get_col_size(); ++icol) + { + upper_mat->get_value(irow, icol) = mat_nspin1->get_value(irow, icol) + + std::complex(clx_i[is], clx_j[is]) * mat_nspin2->get_value(irow, icol); + } + } + //fill the lower triangle matrix + if (PARAM.globalv.domag){ + if (iat1 < iat2) { - upper_mat->get_value(2*irow, 2*icol+1) = mat_nspin_1->get_value(irow, icol) + std::complex(0.0, 1.0) * mat_nspin_2->get_value(irow, icol); - upper_mat->get_value(2*irow+1, 2*icol) = mat_nspin_1->get_value(irow, icol) - std::complex(0.0, 1.0) * mat_nspin_2->get_value(irow, icol); + auto lower_mat = lower_ap->find_matrix(-R_index); + for (int irow = 0; irow < upper_mat->get_row_size(); ++irow) + { + for (int icol = 0; icol < upper_mat->get_col_size(); ++icol) + { + lower_mat->get_value(icol, irow) = conj(upper_mat->get_value(irow, icol)); + } + } } } } + } + } - // fill the lower triangle matrix - if (iat1 < iat2) + //std::cout<<"success"<hR_tmp->set_zero(); + hamilt::transferSerials2Parallels( *hRGint_tmpCd, this->hR_tmp); + for (int iap = 0; iap < hR->size_atom_pairs(); iap++) + { + //std::cout<<"iap: "<get_atom_pair(iap); + const int iat1 = ap->get_atom_i(); + const int iat2 = ap->get_atom_j(); + auto* ap_nspin = this->hR_tmp ->find_pair(iat1, iat2); + for (int ir = 0; ir < ap->get_R_size(); ir++) + { + const auto R_index = ap->get_R_index(ir); + auto upper_mat = ap->find_matrix(R_index); + auto mat_nspin = ap_nspin->find_matrix(R_index); + + // The row size and the col size of upper_matrix is double that of matrix_nspin_0 + for (int irow = 0; irow < mat_nspin->get_row_size(); ++irow) { - auto lower_mat = lower_ap->find_matrix(-R_index); - for (int irow = 0; irow < upper_mat->get_row_size(); ++irow) + for (int icol = 0; icol < mat_nspin->get_col_size(); ++icol) { - for (int icol = 0; icol < upper_mat->get_col_size(); ++icol) - { - lower_mat->get_value(icol, irow) = conj(upper_mat->get_value(irow, icol)); - } + upper_mat->get_value(2*irow+row_set[is], 2*icol+col_set[is]) = + mat_nspin->get_value(irow, icol); } } } } + delete hRGint_tmpCd; } - - // =================================== - // transfer HR from Gint to Veff, std::complex>> - // =================================== -#ifdef __MPI - int size; - MPI_Comm_size(MPI_COMM_WORLD, &size); - if (size == 1) - { - hR->add(*this->hRGintCd); - } - else - { - hamilt::transferSerials2Parallels>(*this->hRGintCd, hR); - } + delete[] iat2iwt; #else - hR->add(*this->hRGintCd); + #endif ModuleBase::timer::tick("Gint_k", "transfer_pvpR"); diff --git a/source/module_lr/utils/gint_move.hpp b/source/module_lr/utils/gint_move.hpp index 6b532073a2..8062ff895d 100644 --- a/source/module_lr/utils/gint_move.hpp +++ b/source/module_lr/utils/gint_move.hpp @@ -45,8 +45,8 @@ Gint& Gint::operator=(Gint&& rhs) // move hR after refactor this->hRGint = rhs.hRGint; rhs.hRGint = nullptr; - this->hRGintCd = rhs.hRGintCd; - rhs.hRGintCd = nullptr; + this->hR_tmp = rhs.hR_tmp; + rhs.hR_tmp = nullptr; for (int i = 0; i < this->DMRGint.size(); i++) { delete this->DMRGint[i]; @@ -60,8 +60,8 @@ Gint& Gint::operator=(Gint&& rhs) this->pvdpRz_reduced = std::move(rhs.pvdpRz_reduced); this->DMRGint = std::move(rhs.DMRGint); this->hRGint_tmp = std::move(rhs.hRGint_tmp); - this->DMRGint_full = rhs.DMRGint_full; - rhs.DMRGint_full = nullptr; + this->DM2D_tmp = rhs.DM2D_tmp; + rhs.DM2D_tmp = nullptr; return *this; }