diff --git a/source/module_hamilt_lcao/module_gint/gint.cpp b/source/module_hamilt_lcao/module_gint/gint.cpp index 6690446bed..490228a2f1 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++) { @@ -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]; @@ -196,10 +194,6 @@ 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); - 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); @@ -254,19 +248,24 @@ void Gint::transfer_DM2DtoGrid(std::vector*> DM2D) { } else // NSPIN=4 case { #ifdef __MPI + // 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]; + + std::vector iat2iwt(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); + pv->set_atomic_trace(iat2iwt.data(), ucell->nat, mg); auto ijr_info = DM2D[0]->get_ijr_info(); this-> DM2D_tmp = new hamilt::HContainer(pv, nullptr, &ijr_info); + this-> DM2D_tmp->set_zero(); ModuleBase::Memory::record("Gint::DM2D_tmp", this->DM2D_tmp->get_memory_size()); for (int is = 0; is < 4; is++){ for (int iap = 0; iap < DM2D[0]->size_atom_pairs(); ++iap) { @@ -275,45 +274,19 @@ void Gint::transfer_DM2DtoGrid(std::vector*> DM2D) { 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* tmp_pointer = this -> DM2D_tmp -> 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) { - switch (is) {//todo: It can be written more compactly - case 0: - for (int icol = 0; icol < ap.get_col_size(); icol += 2) { - *(tmp_pointer)++ = data_full[icol]; - } - data_full += ap.get_col_size() * 2; - break; - case 1: - for (int icol = 0; icol < ap.get_col_size(); icol += 2) { - *(tmp_pointer)++ = data_full[icol + 1]; - } - data_full += ap.get_col_size() * 2; - break; - case 2: - data_full += ap.get_col_size(); - for (int icol = 0; icol < ap.get_col_size(); icol += 2) { - *(tmp_pointer)++ = data_full[icol]; - } - data_full += ap.get_col_size(); - break; - case 3: - data_full += ap.get_col_size(); - for (int icol = 0; icol < ap.get_col_size(); icol += 2) { - *(tmp_pointer)++ = data_full[icol + 1]; - } - data_full += ap.get_col_size(); - break; + double* matrix_out = 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]; } } } } hamilt::transferParallels2Serials( *(this->DM2D_tmp), this->DMRGint[is]); } - // delete iat2iwt; - // delete pv; - // delete this-> DM2D_tmp; #else //this->DMRGint_full = DM2D[0]; #endif diff --git a/source/module_hamilt_lcao/module_gint/gint.h b/source/module_hamilt_lcao/module_gint/gint.h index ad1c7b72fd..d8090df191 100644 --- a/source/module_hamilt_lcao/module_gint/gint.h +++ b/source/module_hamilt_lcao/module_gint/gint.h @@ -258,7 +258,7 @@ 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; 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..41e7b280c7 100644 --- a/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp +++ b/source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp @@ -78,86 +78,113 @@ 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; +#ifdef __MPI + int blacs_ctxt = hR->get_paraV()->blacs_ctxt; + std::vector iat2iwt(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.data(), 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) + this->hR_tmp = new hamilt::HContainer>(pv, nullptr, &ijr_info); + ModuleBase::Memory::record("Gint::hRGintCd", this->hR_tmp->get_memory_size()); + + //select hRGint_tmp + std::vector first = {0, 1, 1, 0}; + std::vector second= {3, 2, 2, 3}; + //select position in the big matrix + std::vector row_set = {0, 0, 1, 1}; + std::vector col_set = {0, 1, 0, 1}; + //construct complex matrix + std::vector clx_i = {1, 0, 0, -1}; + std::vector clx_j = {0, 1, -1, 0}; + for (int is = 0; is < 4; is++){ + if(!PARAM.globalv.domag && (is==1 || is==2)) continue; + this->hR_tmp->set_zero(); + 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++) + { + auto* ap = &hRGint_tmpCd->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++) + { + const auto R_index = upper_ap->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) { - 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); + 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); + } } - } - - 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) + //fill the lower triangle matrix + //When is=0 or 3, the real part does not need conjugation; + //when is=1 or 2, the small matrix is not Hermitian, so conjugation is not needed + if (iat1 < iat2) { - for (int icol = 0; icol < mat_nspin_1->get_col_size(); ++icol) + auto lower_mat = lower_ap->find_matrix(-R_index); + for (int irow = 0; irow < upper_mat->get_row_size(); ++irow) { - 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); + for (int icol = 0; icol < upper_mat->get_col_size(); ++icol) + { + lower_mat->get_value(icol, irow) = upper_mat->get_value(irow, icol); + } } } - } - // fill the lower triangle matrix - if (iat1 < iat2) + } + } + } + // transfer hRGint_tmpCd to parallel hR_tmp + hamilt::transferSerials2Parallels( *hRGint_tmpCd, this->hR_tmp); + // merge hR_tmp to hR + for (int iap = 0; iap < hR->size_atom_pairs(); iap++) + { + auto* ap = &hR->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); } } } } - } - - // =================================== - // 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 hRGint_tmpCd; } #else - hR->add(*this->hRGintCd); -#endif +#endif ModuleBase::timer::tick("Gint_k", "transfer_pvpR"); return; } diff --git a/source/module_lr/utils/gint_move.hpp b/source/module_lr/utils/gint_move.hpp index 2d610cc77e..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];