From 156dac2ee1df143866a4404fa30874cdf5cf95ea Mon Sep 17 00:00:00 2001 From: DylanWRh <93636459+DylanWRh@users.noreply.github.com> Date: Thu, 13 Jun 2024 21:31:06 +0800 Subject: [PATCH] Refactor: change `R_index` in class `AtomPair` to type `ModuleBase::vector3` (#4243) * Refactor `R_index` in `AtomPair` * fix a semicolon * small fix * remove useless line * Refactor `find_R` and `find_matrix` * Refactor `AtomPair` constructor * Refactor `get_R_index` to return reference * rollback --------- Co-authored-by: Mohan Chen --- source/module_base/vector3.h | 6 +- .../module_dm/density_matrix.cpp | 26 ++-- source/module_esolver/io_npz.cpp | 8 +- .../operator_lcao/ekinetic_new.cpp | 4 +- .../operator_lcao/overlap_new.cpp | 7 +- .../operator_lcao/td_ekinetic_lcao.cpp | 8 +- .../operator_lcao/td_nonlocal_lcao.cpp | 4 +- .../operator_lcao/veff_lcao.cpp | 2 +- .../hamilt_lcaodft/spar_hsr.cpp | 12 +- .../module_hamilt_lcao/module_gint/gint.cpp | 4 +- .../module_gint/gint_k_pvpr.cpp | 8 +- .../module_hcontainer/atom_pair.cpp | 144 +++++++++++++----- .../module_hcontainer/atom_pair.h | 31 +++- .../module_hcontainer/func_folding.cpp | 6 +- .../module_hcontainer/hcontainer.cpp | 107 ++++++++++--- .../module_hcontainer/hcontainer.h | 10 +- .../test/test_hcontainer.cpp | 10 +- .../test/test_hcontainer_complex.cpp | 10 +- .../module_hcontainer/transfer.cpp | 7 +- .../module_hcontainer/transfer.h | 1 + source/module_io/fR_overlap.cpp | 4 +- source/module_io/td_current_io.cpp | 7 +- source/module_io/to_wannier90_lcao.cpp | 4 +- 23 files changed, 290 insertions(+), 140 deletions(-) diff --git a/source/module_base/vector3.h b/source/module_base/vector3.h index f499bd3b78..504adf5804 100644 --- a/source/module_base/vector3.h +++ b/source/module_base/vector3.h @@ -33,8 +33,10 @@ template class Vector3 * @param z1 */ Vector3(const T &x1 = 0, const T &y1 = 0, const T &z1 = 0) : x(x1), y(y1), z(z1){}; - Vector3(const Vector3 &v) : x(v.x), y(v.y), z(v.z){}; // Peize Lin add 2018-07-16 - explicit Vector3(const std::array &v) :x(v[0]), y(v[1]), z(v[2]){} + Vector3(const Vector3 &v) : x(v.x), y(v.y), z(v.z){}; // Peize Lin add 2018-07-16 + explicit Vector3(const std::array &v) :x(v[0]), y(v[1]), z(v[2]){} + + Vector3(Vector3 &&v) noexcept : x(v.x), y(v.y), z(v.z) {} /** * @brief set a 3d vector diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp index 31742f60d1..269189662c 100644 --- a/source/module_elecstate/module_dm/density_matrix.cpp +++ b/source/module_elecstate/module_dm/density_matrix.cpp @@ -121,7 +121,7 @@ void DensityMatrix::init_DMR(Grid_Driver* GridD_in, const UnitCell* ucel } ModuleBase::Vector3& R_index = adjs.box[ad]; // std::cout << "R_index: " << R_index.x << " " << R_index.y << " " << R_index.z << std::endl; - hamilt::AtomPair tmp_ap(iat1, iat2, R_index.x, R_index.y, R_index.z, this->_paraV); + hamilt::AtomPair tmp_ap(iat1, iat2, R_index, this->_paraV); tmp_DMR->insert_pair(tmp_ap); } } @@ -237,8 +237,8 @@ void DensityMatrix::init_DMR(const hamilt::HContainer& DMR_in) const int iat2 = DMR_in.get_atom_pair(iap).get_atom_j(); for(int ir = 0;ir tmp_ap(iat1, iat2, R_index[0], R_index[1], R_index[2], paraV_); + const ModuleBase::Vector3 R_index = DMR_in.get_atom_pair(iap).get_R_index(ir); + hamilt::AtomPair tmp_ap(iat1, iat2, R_index, paraV_); tmp_DMR->insert_pair(tmp_ap); } } @@ -417,8 +417,8 @@ void DensityMatrix::cal_DMR_test() } for (int ir = 0; ir < tmp_ap.get_R_size(); ++ir) { - const int* r_index = tmp_ap.get_R_index(ir); - hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 r_index = tmp_ap.get_R_index(ir); + hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index); #ifdef __DEBUG if (tmp_matrix == nullptr) { @@ -484,8 +484,8 @@ void DensityMatrix, double>::cal_DMR() } for (int ir = 0; ir < tmp_ap.get_R_size(); ++ir) { - const int* r_index = tmp_ap.get_R_index(ir); - hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 r_index = tmp_ap.get_R_index(ir); + hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index); #ifdef __DEBUG if (tmp_matrix == nullptr) { @@ -644,12 +644,12 @@ void DensityMatrix::cal_DMR() throw std::string("Atom-pair not belong this process"); } // R index - const int* r_index = tmp_ap.get_R_index(0); + const ModuleBase::Vector3 r_index = tmp_ap.get_R_index(0); #ifdef __DEBUG assert(tmp_ap.get_R_size() == 1); - assert(r_index[0] == 0 && r_index[1] == 0 && r_index[2] == 0); + assert(r_index.x == 0 && r_index.y == 0 && r_index.z == 0); #endif - hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index[0], r_index[1], r_index[2]); + hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index); #ifdef __DEBUG if (tmp_matrix == nullptr) { @@ -695,9 +695,9 @@ void DensityMatrix::sum_DMR_spin() hamilt::AtomPair& tmp_ap_down = tmp_DMR_down->get_atom_pair(i); for (int ir = 0; ir < tmp_ap_up.get_R_size(); ++ir) { - const int* r_index = tmp_ap_up.get_R_index(ir); - hamilt::BaseMatrix* tmp_matrix_up = tmp_ap_up.find_matrix(r_index[0], r_index[1], r_index[2]); - hamilt::BaseMatrix* tmp_matrix_down = tmp_ap_down.find_matrix(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 r_index = tmp_ap_up.get_R_index(ir); + hamilt::BaseMatrix* tmp_matrix_up = tmp_ap_up.find_matrix(r_index); + hamilt::BaseMatrix* tmp_matrix_down = tmp_ap_down.find_matrix(r_index); TR* ptr_up = tmp_matrix_up->get_pointer(); TR* ptr_down = tmp_matrix_down->get_pointer(); for (int i = 0; i < tmp_ap_up.get_size(); ++i) diff --git a/source/module_esolver/io_npz.cpp b/source/module_esolver/io_npz.cpp index 9c58bd7b0b..535d915016 100644 --- a/source/module_esolver/io_npz.cpp +++ b/source/module_esolver/io_npz.cpp @@ -448,9 +448,9 @@ void ESolver_KS_LCAO::output_mat_npz(std::string& zipname, const hamilt: for(int iR=0;iR r_index = HR_serial[0].get_atom_pair(iap).get_R_index(iR); filename = "mat_"+std::to_string(atom_i)+"_"+std::to_string(atom_j)+"_" - +std::to_string(r_index[0])+"_"+std::to_string(r_index[1])+"_"+std::to_string(r_index[2]); + +std::to_string(r_index.x)+"_"+std::to_string(r_index.y)+"_"+std::to_string(r_index.z); std::vector shape = {(size_t)row_size,(size_t)col_size}; cnpy::npz_save(zipname,filename,matrix.get_pointer(),shape,"a"); } @@ -472,10 +472,10 @@ void ESolver_KS_LCAO::output_mat_npz(std::string& zipname, const hamilt: for(int iR=0;iR r_index = hR.get_atom_pair(iap).get_R_index(iR); filename = "mat_"+std::to_string(atom_i)+"_"+std::to_string(atom_j)+"_" - +std::to_string(r_index[0])+"_"+std::to_string(r_index[1])+"_"+std::to_string(r_index[2]); + +std::to_string(r_index.x)+"_"+std::to_string(r_index.y)+"_"+std::to_string(r_index.z); std::vector shape = {(size_t)row_size,(size_t)col_size}; cnpy::npz_save(zipname,filename,matrix.get_pointer(),shape,"a"); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp index e91e812524..3b8069d4d7 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.cpp @@ -84,7 +84,7 @@ void hamilt::EkineticNew>::initialize_HR(Grid_Drive const int I2 = adjs.natom[ad]; int iat2 = ucell->itia2iat(T2, I2); ModuleBase::Vector3& R_index = adjs.box[ad]; - hamilt::AtomPair tmp(iat1, iat2, R_index.x, R_index.y, R_index.z, paraV); + hamilt::AtomPair tmp(iat1, iat2, R_index, paraV); this->hR->insert_pair(tmp); } } @@ -122,7 +122,7 @@ void hamilt::EkineticNew>::calculate_HR() const ModuleBase::Vector3& R_index2 = adjs.box[ad]; ModuleBase::Vector3 dtau = this->ucell->cal_dtau(iat1, iat2, R_index2); - hamilt::BaseMatrix* tmp = this->HR_fixed->find_matrix(iat1, iat2, R_index2.x, R_index2.y, R_index2.z); + hamilt::BaseMatrix* tmp = this->HR_fixed->find_matrix(iat1, iat2, R_index2); if (tmp != nullptr) { this->cal_HR_IJR(iat1, iat2, paraV, dtau, tmp->get_pointer()); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp index 9c16eefb8e..42d46f1332 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp @@ -66,7 +66,7 @@ void hamilt::OverlapNew>::initialize_SR(Grid_Driver { continue; } - hamilt::AtomPair tmp(iat1, iat2, R_index.x, R_index.y, R_index.z, paraV); + hamilt::AtomPair tmp(iat1, iat2, R_index, paraV); SR->insert_pair(tmp); } } @@ -92,9 +92,8 @@ void hamilt::OverlapNew>::calculate_SR() for (int iR = 0; iR < tmp.get_R_size(); ++iR) { - const int* R_index = tmp.get_R_index(iR); - ModuleBase::Vector3 R_vector(R_index[0], R_index[1], R_index[2]); - auto dtau = ucell->cal_dtau(iat1, iat2, R_vector); + const ModuleBase::Vector3 R_index = tmp.get_R_index(iR); + auto dtau = ucell->cal_dtau(iat1, iat2, R_index); TR* data_pointer = tmp.get_pointer(iR); this->cal_SR_IJR(iat1, iat2, paraV, dtau, data_pointer); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp index 13ad9b744c..d20cbdfcf1 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp @@ -92,8 +92,8 @@ void TDEkinetic>::calculate_HR() const ModuleBase::Vector3& R_index2 = adjs.box[ad]; ModuleBase::Vector3 dtau = this->ucell->cal_dtau(iat1, iat2, R_index2); - hamilt::BaseMatrix>* tmp = this->hR_tmp->find_matrix(iat1, iat2, R_index2.x, R_index2.y, R_index2.z); - hamilt::BaseMatrix* tmp1 = this->SR->find_matrix(iat1, iat2, R_index2.x, R_index2.y, R_index2.z); + hamilt::BaseMatrix>* tmp = this->hR_tmp->find_matrix(iat1, iat2, R_index2); + hamilt::BaseMatrix* tmp1 = this->SR->find_matrix(iat1, iat2, R_index2); if (tmp != nullptr) { this->cal_HR_IJR(iat1, iat2, paraV, dtau, tmp->get_pointer(),tmp1->get_pointer()); @@ -309,11 +309,11 @@ void TDEkinetic>::initialize_HR_tmp(const Parallel_Orbitals hamilt::AtomPair& tmp = this->hR->get_atom_pair(i); for(int ir = 0;ir < tmp.get_R_size(); ++ir ) { - const int* R_index = tmp.get_R_index(ir); + const ModuleBase::Vector3 R_index = tmp.get_R_index(ir); const int iat1 = tmp.get_atom_i(); const int iat2 = tmp.get_atom_j(); - hamilt::AtomPair> tmp1(iat1, iat2, R_index[0], R_index[1], R_index[2], paraV); + hamilt::AtomPair> tmp1(iat1, iat2, R_index, paraV); this->hR_tmp->insert_pair(tmp1); } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp index 4576596425..42335ac093 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp @@ -113,11 +113,11 @@ void hamilt::TDNonlocal>::initialize_HR_tmp(const P hamilt::AtomPair& tmp = this->hR->get_atom_pair(i); for(int ir = 0;ir < tmp.get_R_size(); ++ir ) { - const int* R_index = tmp.get_R_index(ir); + const ModuleBase::Vector3 R_index = tmp.get_R_index(ir); const int iat1 = tmp.get_atom_i(); const int iat2 = tmp.get_atom_j(); - hamilt::AtomPair> tmp1(iat1, iat2, R_index[0], R_index[1], R_index[2], paraV); + hamilt::AtomPair> tmp1(iat1, iat2, R_index, paraV); this->hR_tmp->insert_pair(tmp1); } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp index 82fa242e61..48e843ea72 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.cpp @@ -44,7 +44,7 @@ void Veff>::initialize_HR(const UnitCell* ucell_in, if (ucell_in->cal_dtau(iat1, iat2, R_index2).norm() * ucell_in->lat0 < orb.Phi[T1].getRcut() + orb.Phi[T2].getRcut()) { - hamilt::AtomPair tmp(iat1, iat2, R_index2.x, R_index2.y, R_index2.z, paraV); + hamilt::AtomPair tmp(iat1, iat2, R_index2, paraV); this->hR->insert_pair(tmp); } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp index 162f7857e5..6a1902899f 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp @@ -150,8 +150,8 @@ void sparse_format::cal_HContainer_d( for(int iR=0;iR dR(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 r_index = hR.get_atom_pair(iap).get_R_index(iR); + Abfs::Vector3_Order dR(r_index.x, r_index.y, r_index.z); for(int i=0;i dR(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 r_index = hR.get_atom_pair(iap).get_R_index(iR); + Abfs::Vector3_Order dR(r_index.x, r_index.y, r_index.z); for(int i=0;i dR(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 r_index = hR.get_atom_pair(iap).get_R_index(iR); + Abfs::Vector3_Order dR(r_index.x, r_index.y, r_index.z); for(int i=0;i*> DM2D) int iat2 = ap.get_atom_j(); for(int ir = 0;ir 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[0], r_index[1], r_index[2])->get_pointer(); + 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 *hR,const UnitCell* ucell_ int ixxx = DM_start + this->gridt->find_R2st[iat][nad]; - hamilt::BaseMatrix* tmp_matrix = this->hRGint->find_matrix(iat, iat2, dR.x, dR.y, dR.z); + hamilt::BaseMatrix* tmp_matrix = this->hRGint->find_matrix(iat, iat2, dR); #ifdef __DEBUG assert(tmp_matrix != nullptr); #endif @@ -508,7 +508,7 @@ void Gint_k::transfer_pvpR(hamilt::HContainer *hR,const UnitCell* ucell_ // save the lower triangle. if(iat < iat2)//skip iat == iat2 { - hamilt::BaseMatrix* conj_matrix = this->hRGint->find_matrix(iat2, iat, -dR.x, -dR.y, -dR.z); + hamilt::BaseMatrix* conj_matrix = this->hRGint->find_matrix(iat2, iat, -dR); #ifdef __DEBUG assert(conj_matrix != nullptr); #endif @@ -617,7 +617,7 @@ void Gint_k::transfer_pvpR(hamilt::HContainer> *hR,const Un int ixxx = DM_start + this->gridt->find_R2st[iat][nad]; - hamilt::BaseMatrix>* tmp_matrix = this->hRGintCd->find_matrix(iat, iat2, dR.x, dR.y, dR.z); + hamilt::BaseMatrix>* tmp_matrix = this->hRGintCd->find_matrix(iat, iat2, dR); #ifdef __DEBUG assert(tmp_matrix != nullptr); #endif @@ -666,7 +666,7 @@ void Gint_k::transfer_pvpR(hamilt::HContainer> *hR,const Un // save the lower triangle. if(iat < iat2) { - hamilt::BaseMatrix>* conj_matrix = this->hRGintCd->find_matrix(iat2, iat, -dR.x, -dR.y, -dR.z); + hamilt::BaseMatrix>* conj_matrix = this->hRGintCd->find_matrix(iat2, iat, -dR); #ifdef __DEBUG assert(conj_matrix != nullptr); #endif diff --git a/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp b/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp index cd0ead0435..d35553a2f8 100644 --- a/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp @@ -30,7 +30,8 @@ AtomPair::AtomPair(const int& atom_i_, const int& atom_j_, const Parallel_Orb } this->row_size = this->paraV->get_row_size(atom_i); this->col_size = this->paraV->get_col_size(atom_j); - this->R_index.resize(3, 0); + this->R_index.resize(0); + this->R_index.push_back(ModuleBase::Vector3(0, 0, 0)); this->current_R = 0; if (existed_matrix != nullptr) { @@ -63,11 +64,9 @@ AtomPair::AtomPair(const int& atom_i_, } this->row_size = this->paraV->get_row_size(atom_i); this->col_size = this->paraV->get_col_size(atom_j); - this->R_index.resize(3, 0); + this->R_index.resize(0); + this->R_index.push_back(ModuleBase::Vector3(rx, ry, rz)); this->current_R = 0; - this->R_index[0] = rx; - this->R_index[1] = ry; - this->R_index[2] = rz; if (existed_matrix != nullptr) { BaseMatrix tmp(row_size, col_size, existed_matrix); @@ -79,6 +78,39 @@ AtomPair::AtomPair(const int& atom_i_, this->values.push_back(tmp); } } + +template +AtomPair::AtomPair(const int& atom_i_, + const int& atom_j_, + const ModuleBase::Vector3 &R_index, + const Parallel_Orbitals* paraV_, + T* existed_matrix) + : atom_i(atom_i_), atom_j(atom_j_), paraV(paraV_) +{ + assert(this->paraV != nullptr); + this->row_ap = this->paraV->atom_begin_row[atom_i]; + this->col_ap = this->paraV->atom_begin_col[atom_j]; + if (this->row_ap == -1 || this->col_ap == -1) + { + throw std::string("Atom-pair not belong this process"); + } + this->row_size = this->paraV->get_row_size(atom_i); + this->col_size = this->paraV->get_col_size(atom_j); + this->R_index.resize(0); + this->R_index.push_back(ModuleBase::Vector3(R_index)); + this->current_R = 0; + if (existed_matrix != nullptr) + { + BaseMatrix tmp(row_size, col_size, existed_matrix); + this->values.push_back(tmp); + } + else + { + BaseMatrix tmp(row_size, col_size); + this->values.push_back(tmp); + } +} + // direct save whole matrix of atom-pair template AtomPair::AtomPair(const int& atom_i_, @@ -94,7 +126,8 @@ AtomPair::AtomPair(const int& atom_i_, this->col_ap = col_atom_begin[atom_j]; this->row_size = row_atom_begin[atom_i + 1] - row_atom_begin[atom_i]; this->col_size = col_atom_begin[atom_j + 1] - col_atom_begin[atom_j]; - this->R_index.resize(3, 0); + this->R_index.resize(0); + this->R_index.push_back(ModuleBase::Vector3(0, 0, 0)); this->current_R = 0; if (existed_matrix != nullptr) { @@ -125,11 +158,9 @@ AtomPair::AtomPair(const int& atom_i_, this->col_ap = col_atom_begin[atom_j]; this->row_size = row_atom_begin[atom_i + 1] - row_atom_begin[atom_i]; this->col_size = col_atom_begin[atom_j + 1] - col_atom_begin[atom_j]; - this->R_index.resize(3, 0); + this->R_index.resize(0); + this->R_index.push_back(ModuleBase::Vector3(rx, ry, rz)); this->current_R = 0; - this->R_index[0] = rx; - this->R_index[1] = ry; - this->R_index[2] = rz; if (existed_matrix != nullptr) { BaseMatrix tmp(row_size, col_size, existed_matrix); @@ -365,9 +396,7 @@ BaseMatrix& AtomPair::get_HR_values(int rx_in, int ry_in, int rz_in) return this->values[r_index]; } // if not found, add a new BaseMatrix for this R index - R_index.push_back(rx_in); - R_index.push_back(ry_in); - R_index.push_back(rz_in); + R_index.push_back(ModuleBase::Vector3(rx_in, ry_in, rz_in)); values.push_back(BaseMatrix(this->row_size, this->col_size)); values.back().allocate(nullptr, true); // return the last BaseMatrix reference in values @@ -407,12 +436,27 @@ BaseMatrix& AtomPair::get_HR_values(const int& index) const template int AtomPair::find_R(const int& rx_in, const int& ry_in, const int& rz_in) const { - for (int i = 0; i < this->R_index.size(); i += 3) + for (int i = 0; i < this->R_index.size(); i++) + { + if (R_index[i].x == rx_in && R_index[i].y == ry_in && R_index[i].z == rz_in) + { + this->current_R = i; + return i; + } + } + return (-1); +} + +// find_R +template +int AtomPair::find_R(const ModuleBase::Vector3& R_in) const +{ + for (int i = 0; i < this->R_index.size(); i++) { - if (R_index[i] == rx_in && R_index[i + 1] == ry_in && R_index[i + 2] == rz_in) + if (R_index[i] == R_in) { - this->current_R = i / 3; - return (i / 3); + this->current_R = i; + return i; } } return (-1); @@ -448,6 +492,33 @@ BaseMatrix* AtomPair::find_matrix(const int& rx_in, const int& ry_in, cons } } +// find_matrix +template +const BaseMatrix* AtomPair::find_matrix(const ModuleBase::Vector3& R_in) const { + const int r_index = this->find_R(R_in); + if(r_index == -1) + { + return nullptr; + } + else + { + return &(this->values[r_index]); + } +} + +template +BaseMatrix* AtomPair::find_matrix(const ModuleBase::Vector3& R_in){ + const int r_index = this->find_R(R_in); + if(r_index == -1) + { + return nullptr; + } + else + { + return &(this->values[r_index]); + } +} + template void AtomPair::convert_add(const BaseMatrix& target, int rx_in, int ry_in, int rz_in) { @@ -467,21 +538,19 @@ void AtomPair::merge(const AtomPair& other, bool skip_R) throw std::string("AtomPair::merge: atom pair not match"); } int rx = 0, ry = 0, rz = 0; - for (int i = 0; i < other.R_index.size(); i += 3) + for (int i = 0; i < other.R_index.size(); i ++) { if (!skip_R) { - rx = other.R_index[i]; - ry = other.R_index[i + 1]; - rz = other.R_index[i + 2]; + rx = other.R_index[i].x; + ry = other.R_index[i].y; + rz = other.R_index[i].z; } - const BaseMatrix& matrix_tmp = other.get_HR_values(i / 3); + const BaseMatrix& matrix_tmp = other.get_HR_values(i); //if not found, push_back this BaseMatrix to this->values if (this->find_R(rx, ry, rz) == -1) { - this->R_index.push_back(rx); - this->R_index.push_back(ry); - this->R_index.push_back(rz); + this->R_index.push_back(ModuleBase::Vector3(rx, ry, rz)); this->values.push_back(matrix_tmp); } //if found but not allocated, skip this BaseMatrix values @@ -503,7 +572,8 @@ void AtomPair::merge_to_gamma() { // reset R_index to (0, 0, 0) this->R_index.clear(); - this->R_index.resize(3, 0); + this->R_index.resize(0); + this->R_index.push_back(ModuleBase::Vector3(0, 0, 0)); // merge all values to first BaseMatrix BaseMatrix tmp(this->row_size, this->col_size); bool empty = true; @@ -716,29 +786,27 @@ std::tuple, T*> AtomPair::get_matrix_values(int ir) const return std::tuple, T*>({this->row_ap, this->row_size, this->col_ap, this->col_size}, this->values[ir].get_pointer()); } -// interface for get (rx, ry, rz) of index-th R-index in this->R_index, the return should be int[3] +// interface for get (rx, ry, rz) of index-th R-index in this->R_index, the return should be ModuleBase::Vector3 template -int* AtomPair::get_R_index(const int& index) const +ModuleBase::Vector3 AtomPair::get_R_index(const int& index) const { - if (index >= R_index.size() / 3 || index < 0) + if (index >= R_index.size() || index < 0) { - std::cout << "Error: index out of range in get_R_index" << std::endl; - return nullptr; + std::cout << "Error: index out of range in s" << std::endl; + return ModuleBase::Vector3(-1, -1, -1); } else { - // return the (int*) pointer of R_index[index*3] - int* ptr = const_cast(&(R_index[index * 3])); - return ptr; + // return the ModuleBase::Vector3 of R_index[index] + return R_index[index]; } } template -int* AtomPair::get_R_index() const +ModuleBase::Vector3 AtomPair::get_R_index() const { - // return the (int*) pointer of R_index[index*3] - int* ptr = const_cast(&(R_index[current_R * 3])); - return ptr; + // return the ModuleBase::Vector3 of R_index[index] + return R_index[current_R]; } // get_value diff --git a/source/module_hamilt_lcao/module_hcontainer/atom_pair.h b/source/module_hamilt_lcao/module_hcontainer/atom_pair.h index d5f7d69af9..fac117f9e7 100644 --- a/source/module_hamilt_lcao/module_hcontainer/atom_pair.h +++ b/source/module_hamilt_lcao/module_hcontainer/atom_pair.h @@ -3,6 +3,7 @@ // #include "module_cell/atom_spec.h" #include "base_matrix.h" +#include "module_base/vector3.h" #include "module_basis/module_ao/parallel_orbitals.h" #include @@ -62,6 +63,17 @@ class AtomPair T* existed_array = nullptr // if nullptr, new memory will be allocated, otherwise this class is a data wrapper ); + // Constructor of class AtomPair + // Only for 2d-block MPI parallel case + // This constructor used for initialize a atom-pair local Hamiltonian with non-zero cell indexes, + // which is used for constructing HR (real space Hamiltonian) objects. + AtomPair(const int& atom_i_, // atomic index of atom i, used to identify atom + const int& atom_j_, // atomic index of atom j, used to identify atom + const ModuleBase::Vector3& R_index, // xyz coordinates of cell + const Parallel_Orbitals* paraV_, // information for 2d-block parallel + T* existed_array + = nullptr // if nullptr, new memory will be allocated, otherwise this class is a data wrapper + ); // This constructor used for initialize a atom-pair local Hamiltonian with only center cell // which is used for constructing HK (k space Hamiltonian) objects, (gamma_only case) AtomPair(const int& atom_i, // atomic index of atom i, used to identify atom @@ -163,16 +175,19 @@ class AtomPair */ BaseMatrix& get_HR_values(const int& index) const; - // interface for get (rx, ry, rz) of index-th R-index in this->R_index, the return should be int[3] - int* get_R_index(const int& index) const; - // interface for get (rx, ry, rz) of current_R, the return should be int[3] - int* get_R_index() const; + // interface for get (rx, ry, rz) of index-th R-index in this->R_index, the return should be ModuleBase::Vector3 + ModuleBase::Vector3 get_R_index(const int& index) const; + // interface for get (rx, ry, rz) of current_R, the return should be ModuleBase::Vector3 + ModuleBase::Vector3 get_R_index() const; // interface for search (rx, ry, rz) in this->R_index, if found, current_R would be set to index int find_R(const int& rx_in, const int& ry_in, const int& rz_in) const; + int find_R(const ModuleBase::Vector3& R_in) const; // interface for search (rx, ry, rz) in this->R_index, if found, current_R would be set to index // and return BaseMatrix* of this->values[index] const BaseMatrix* find_matrix(const int& rx_in, const int& ry_in, const int& rz_in) const; BaseMatrix* find_matrix(const int& rx_in, const int& ry_in, const int& rz_in); + const BaseMatrix* find_matrix(const ModuleBase::Vector3& R_in) const; + BaseMatrix* find_matrix(const ModuleBase::Vector3& R_in); // this interface will call get_value in this->values // these four interface can be used only when R-index has been choosed (current_R >= 0) @@ -261,10 +276,10 @@ class AtomPair size_t get_R_size() const { #ifdef __DEBUG - assert(this->R_index.size() / 3 == this->values.size()); - assert(this->R_index.size() % 3 == 0); + assert(this->R_index.size() == this->values.size()); + // assert(this->R_index.size() % 3 == 0); #endif - return this->R_index.size() / 3; + return this->R_index.size(); } /** @@ -274,7 +289,7 @@ class AtomPair private: // it contains 3 index of cell, size of R_index is three times of values. - std::vector R_index; + std::vector> R_index; // it contains containers for accessing matrix of this atom-pair std::vector> values; diff --git a/source/module_hamilt_lcao/module_hcontainer/func_folding.cpp b/source/module_hamilt_lcao/module_hcontainer/func_folding.cpp index d9e76eb733..af8f7c76d5 100644 --- a/source/module_hamilt_lcao/module_hcontainer/func_folding.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/func_folding.cpp @@ -26,16 +26,16 @@ void folding_HR(const hamilt::HContainer& hR, hamilt::AtomPair& tmp = hR.get_atom_pair(i); for(int ir = 0;ir < tmp.get_R_size(); ++ir ) { - const int* r_index = tmp.get_R_index(ir); + const ModuleBase::Vector3 r_index = tmp.get_R_index(ir); // cal k_phase // if TK==std::complex, kphase is e^{ikR} - const ModuleBase::Vector3 dR(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 dR(r_index.x, r_index.y, r_index.z); const double arg = (kvec_d_in * dR) * ModuleBase::TWO_PI; double sinp, cosp; ModuleBase::libm::sincos(arg, &sinp, &cosp); std::complex kphase = std::complex(cosp, sinp); - tmp.find_R(r_index[0], r_index[1], r_index[2]); + tmp.find_R(r_index); tmp.add_to_matrix(hk, ncol, kphase, hk_type); } } diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp index 4d31335585..8ddb50bb0e 100644 --- a/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp @@ -252,6 +252,49 @@ BaseMatrix* HContainer::find_matrix(int atom_i, int atom_j, int rx, int ry } } +template +BaseMatrix* HContainer::find_matrix(int atom_i, int atom_j, const ModuleBase::Vector3& R_index) +{ + AtomPair* tmp = this->find_pair(atom_i, atom_j); + if(tmp == nullptr) + { + return nullptr; + } + else + { + if(this->gamma_only) + { + return tmp->find_matrix(0, 0, 0); + } + else + { + return tmp->find_matrix(R_index); + } + } +} + +template +const BaseMatrix* HContainer::find_matrix(int atom_i, int atom_j,const ModuleBase::Vector3& R_index) const +{ + AtomPair* tmp = this->find_pair(atom_i, atom_j); + if(tmp == nullptr) + { + return nullptr; + } + else + { + if(this->gamma_only) + { + return tmp->find_matrix(0, 0, 0); + } + else + { + return tmp->find_matrix(R_index); + } + } +} + + // get_atom_pair with atom_ij template AtomPair& HContainer::get_atom_pair(int atom_i, int atom_j) const @@ -393,10 +436,28 @@ int HContainer::find_R(const int& rx_in, const int& ry_in, const int& rz_in) { return -1; } - for (int i = 0; i < this->tmp_R_index.size() / 3; i++) + for (int i = 0; i < this->tmp_R_index.size(); i++) + { + if (this->tmp_R_index[i].x == rx_in && this->tmp_R_index[i].y == ry_in + && this->tmp_R_index[i].z == rz_in) + { + return i; + } + } + return -1; +} + +template +int HContainer::find_R(const ModuleBase::Vector3& R_in) const +{ + // search R_in in this->tmp_R_index + if (this->tmp_R_index.empty()) + { + return -1; + } + for (int i = 0; i < this->tmp_R_index.size(); i++) { - if (this->tmp_R_index[i * 3] == rx_in && this->tmp_R_index[i * 3 + 1] == ry_in - && this->tmp_R_index[i * 3 + 2] == rz_in) + if (this->tmp_R_index[i] == R_in) { return i; } @@ -409,9 +470,9 @@ template size_t HContainer::size_R_loop() const { // R index is fixed - if (this->current_R > -1 && this->tmp_R_index.size() > 2) + if (this->current_R > -1 && this->tmp_R_index.size() > 0) { - return this->tmp_R_index.size()/3; + return this->tmp_R_index.size(); } /** * start a new iteration of loop_R @@ -422,40 +483,38 @@ size_t HContainer::size_R_loop() const for (auto it = this->atom_pairs.begin(); it != this->atom_pairs.end(); ++it) { /** - * search (rx, ry, rz) with (it->R_values[i*3+0], it->R_values[i*3+1], it->R_values[i*3+2]) + * search (rx, ry, rz) with (it->R_values[i].x, it->R_values[i].y, it->R_values[i].z) * if (rx, ry, rz) not found in this->tmp_R_index, * insert the (rx, ry, rz) into end of this->tmp_R_index * no need to sort this->tmp_R_index, using find_R() to find the (rx, ry, rz) -> int in tmp_R_index */ for (int iR = 0; iR < it->get_R_size(); iR++) { - int* R_pointer = it->get_R_index(iR); - int it_tmp = this->find_R(R_pointer[0], R_pointer[1], R_pointer[2]); + const ModuleBase::Vector3 r_vec = it->get_R_index(iR); + int it_tmp = this->find_R(r_vec); if (it_tmp == -1) { - this->tmp_R_index.push_back(R_pointer[0]); - this->tmp_R_index.push_back(R_pointer[1]); - this->tmp_R_index.push_back(R_pointer[2]); + this->tmp_R_index.push_back(ModuleBase::Vector3(r_vec)); } } } - return this->tmp_R_index.size() / 3; + return this->tmp_R_index.size(); } template void HContainer::loop_R(const size_t& index, int& rx, int& ry, int& rz) const { #ifdef __DEBUG - if (index >= this->tmp_R_index.size() / 3) + if (index >= this->tmp_R_index.size()) { std::cout << "Error: index out of range in loop_R" << std::endl; exit(1); } #endif // set rx, ry, rz - rx = this->tmp_R_index[index * 3]; - ry = this->tmp_R_index[index * 3 + 1]; - rz = this->tmp_R_index[index * 3 + 2]; + rx = this->tmp_R_index[index].x; + ry = this->tmp_R_index[index].y; + rz = this->tmp_R_index[index].z; return; } @@ -597,7 +656,7 @@ size_t HContainer::get_memory_size() const memory += this->sparse_ap_index[i].capacity() * sizeof(int); } memory += this->tmp_atom_pairs.capacity() * sizeof(AtomPair*); - memory += this->tmp_R_index.capacity() * sizeof(int); + memory += this->tmp_R_index.capacity() * sizeof(ModuleBase::Vector3); if(this->allocated) { memory += this->get_nnr() * sizeof(T); @@ -630,15 +689,15 @@ void HContainer::shape_synchron( const HContainer& other) { for(int ir = 0;ir < other.atom_pairs[i].get_R_size();++ir) { - int* R_pointer = other.atom_pairs[i].get_R_index(ir); - if(tmp_pointer->find_R(R_pointer[0], R_pointer[1], R_pointer[2]) != -1) + const ModuleBase::Vector3 R_vec = other.atom_pairs[i].get_R_index(ir); + if(tmp_pointer->find_R(R_vec) != -1) { // do nothing } else { // insert the new BaseMatrix - tmp_pointer->get_HR_values(R_pointer[0], R_pointer[1], R_pointer[2]); + tmp_pointer->get_HR_values(R_vec.x, R_vec.y, R_vec.z); } } } @@ -667,10 +726,10 @@ std::vector HContainer::get_ijr_info() const // loop R for (int ir = 0; ir < number_R; ++ir) { - int* R_pointer = this->atom_pairs[i].get_R_index(ir); - ijr_info.push_back(R_pointer[0]); - ijr_info.push_back(R_pointer[1]); - ijr_info.push_back(R_pointer[2]); + const ModuleBase::Vector3 R_vec = this->atom_pairs[i].get_R_index(ir); + ijr_info.push_back(R_vec.x); + ijr_info.push_back(R_vec.y); + ijr_info.push_back(R_vec.z); } } return ijr_info; diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer.h b/source/module_hamilt_lcao/module_hcontainer/hcontainer.h index 43eddf0f53..3efbdd6153 100644 --- a/source/module_hamilt_lcao/module_hcontainer/hcontainer.h +++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer.h @@ -5,6 +5,7 @@ #include #include "atom_pair.h" +#include "module_base/vector3.h" #include "module_cell/unitcell.h" namespace hamilt @@ -115,8 +116,8 @@ namespace hamilt * // loop R-index * for (int iR = 0; iR < atom_ij.size_R(); iR++) * { - * const int* r_index = atom_ij.get_R_index(iR); - * auto tmp_matrix = atom_ij.get_HR_values(r_index[0], r_index[1], r_index[2]); + * const ModuleBase::Vector3 r_index = atom_ij.get_R_index(iR); + * auto tmp_matrix = atom_ij.get_HR_values(r_index.x, r_index.y, r_index.z); * // do something with tmp_matrix * ... * } @@ -219,6 +220,8 @@ class HContainer */ BaseMatrix* find_matrix(int i, int j, int rx, int ry, int rz); const BaseMatrix* find_matrix(int i, int j, int rx, int ry, int rz) const; + BaseMatrix* find_matrix(int i, int j, const ModuleBase::Vector3& R_index); + const BaseMatrix* find_matrix(int i, int j, const ModuleBase::Vector3& R_index) const; /** * @brief return a reference of AtomPair with index of atom I and J in atom_pairs @@ -420,7 +423,7 @@ class HContainer */ mutable std::vector*> tmp_atom_pairs; // it contains 3 index of cell, size of R_index is three times of values. - mutable std::vector tmp_R_index; + mutable std::vector> tmp_R_index; // current index of R in tmp_atom_pairs, -1 means not initialized mutable int current_R = -1; /** @@ -432,6 +435,7 @@ class HContainer * @return int index of R in tmp_R_index */ int find_R(const int& rx_in, const int& ry_in, const int& rz_in) const; + int find_R(const ModuleBase::Vector3& R_in) const; bool gamma_only = false; diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer.cpp index 8473451dcc..b76b10db09 100644 --- a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer.cpp @@ -307,11 +307,11 @@ TEST_F(HContainerTest, size_atom_pairs) EXPECT_EQ(HR->get_atom_pair(0).get_atom_j(), 1); EXPECT_EQ(HR->get_atom_pair(0).get_row_size(), 2); EXPECT_EQ(HR->get_atom_pair(0).get_col_size(), 2); - int* R_ptr = HR->get_atom_pair(0).get_R_index(); - EXPECT_EQ(R_ptr[0], 1); - EXPECT_EQ(R_ptr[1], 0); - EXPECT_EQ(R_ptr[2], 0); - EXPECT_EQ(HR->get_atom_pair(0).get_R_index(5), nullptr); + const ModuleBase::Vector3 R_ptr = HR->get_atom_pair(0).get_R_index(); + EXPECT_EQ(R_ptr.x, 1); + EXPECT_EQ(R_ptr.y, 0); + EXPECT_EQ(R_ptr.z, 0); + EXPECT_EQ(HR->get_atom_pair(0).get_R_index(5), ModuleBase::Vector3(-1, -1, -1)); // check if data is correct double* data_ptr = HR->get_atom_pair(0).get_pointer(); EXPECT_EQ(data_ptr[0], 1); diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp index d55feb73c7..b0aa7a22b1 100644 --- a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp @@ -310,11 +310,11 @@ TEST_F(HContainerTest, size_atom_pairs) EXPECT_EQ(HR->get_atom_pair(0).get_atom_j(), 1); EXPECT_EQ(HR->get_atom_pair(0).get_row_size(), 2); EXPECT_EQ(HR->get_atom_pair(0).get_col_size(), 2); - int* R_ptr = HR->get_atom_pair(0).get_R_index(); - EXPECT_EQ(R_ptr[0], 1); - EXPECT_EQ(R_ptr[1], 0); - EXPECT_EQ(R_ptr[2], 0); - EXPECT_EQ(HR->get_atom_pair(0).get_R_index(5), nullptr); + const ModuleBase::Vector3 R_ptr = HR->get_atom_pair(0).get_R_index(); + EXPECT_EQ(R_ptr.x, 1); + EXPECT_EQ(R_ptr.y, 0); + EXPECT_EQ(R_ptr.z, 0); + EXPECT_EQ(HR->get_atom_pair(0).get_R_index(5), ModuleBase::Vector3(-1, -1, -1)); // check if data is correct std::complex* data_ptr = HR->get_atom_pair(0).get_pointer(); EXPECT_EQ(data_ptr[0], std::complex(1)); diff --git a/source/module_hamilt_lcao/module_hcontainer/transfer.cpp b/source/module_hamilt_lcao/module_hcontainer/transfer.cpp index 08a866c27f..0d984ab5c4 100644 --- a/source/module_hamilt_lcao/module_hcontainer/transfer.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/transfer.cpp @@ -358,12 +358,13 @@ void HTransSerial::cal_ap_indexes(int irank, std::vector* ap_indexes) // loop of R for (int k = 0; k < atom_pair.get_R_size(); k++) { + const ModuleBase::Vector3 r_index = atom_pair.get_R_index(k); // rx - *data++ = atom_pair.get_R_index(k)[0]; + *data++ = r_index.x; // ry - *data++ = atom_pair.get_R_index(k)[1]; + *data++ = r_index.y; // rz - *data++ = atom_pair.get_R_index(k)[2]; + *data++ = r_index.z; } } } diff --git a/source/module_hamilt_lcao/module_hcontainer/transfer.h b/source/module_hamilt_lcao/module_hcontainer/transfer.h index deea0f5d7c..35c41e4350 100644 --- a/source/module_hamilt_lcao/module_hcontainer/transfer.h +++ b/source/module_hamilt_lcao/module_hcontainer/transfer.h @@ -4,6 +4,7 @@ #include #include +#include "module_base/vector3.h" #include "./hcontainer.h" #ifdef __MPI diff --git a/source/module_io/fR_overlap.cpp b/source/module_io/fR_overlap.cpp index 2786debc91..ee0e909b0a 100644 --- a/source/module_io/fR_overlap.cpp +++ b/source/module_io/fR_overlap.cpp @@ -126,8 +126,8 @@ void FR_overlap::calculate_FR() for (int iR = 0; iR < tmp.get_R_size(); ++iR) { - const int* R_index = tmp.get_R_index(iR); - ModuleBase::Vector3 R_vector(R_index[0], R_index[1], R_index[2]); + const ModuleBase::Vector3 R_index = tmp.get_R_index(iR); + ModuleBase::Vector3 R_vector(R_index); auto dtau = ucell->cal_dtau(iat1, iat2, R_vector) * ucell->lat0; T* data_pointer = tmp.get_pointer(iR); this->cal_FR_IJR(iat1, iat2, paraV, dtau, data_pointer); diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp index 7e38fc248f..805fed213d 100644 --- a/source/module_io/td_current_io.cpp +++ b/source/module_io/td_current_io.cpp @@ -2,6 +2,7 @@ #include "module_base/global_function.h" #include "module_base/global_variable.h" +#include "module_base/vector3.h" #include "module_base/timer.h" #include "module_base/libm/libm.h" #include "module_base/tool_threading.h" @@ -92,8 +93,8 @@ void ModuleIO::cal_tmp_DM(elecstate::DensityMatrix, double> } for (int ir = 0; ir < tmp_ap.get_R_size(); ++ir) { - const int* r_index = tmp_ap.get_R_index(ir); - hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 r_index = tmp_ap.get_R_index(ir); + hamilt::BaseMatrix* tmp_matrix = tmp_ap.find_matrix(r_index); #ifdef __DEBUG if (tmp_matrix == nullptr) { @@ -106,7 +107,7 @@ void ModuleIO::cal_tmp_DM(elecstate::DensityMatrix, double> { // cal k_phase // if TK==std::complex, kphase is e^{ikR} - const ModuleBase::Vector3 dR(r_index[0], r_index[1], r_index[2]); + const ModuleBase::Vector3 dR(r_index.x, r_index.y, r_index.z); const double arg = (DM.get_kv_pointer()->kvec_d[ik] * dR) * ModuleBase::TWO_PI; double sinp, cosp; ModuleBase::libm::sincos(arg, &sinp, &cosp); diff --git a/source/module_io/to_wannier90_lcao.cpp b/source/module_io/to_wannier90_lcao.cpp index 6b0103c458..b67adab984 100644 --- a/source/module_io/to_wannier90_lcao.cpp +++ b/source/module_io/to_wannier90_lcao.cpp @@ -381,8 +381,8 @@ void toWannier90_LCAO::unkdotkb( for(int iR = 0; iR < tmp_FR_container->get_atom_pair(iap).get_R_size(); ++iR) { auto& matrix = tmp_FR_container->get_atom_pair(iap).get_HR_values(iR); - int* r_index = tmp_FR_container->get_atom_pair(iap).get_R_index(iR); - ModuleBase::Vector3 dR = ModuleBase::Vector3(r_index[0], r_index[1], r_index[2]) * GlobalC::ucell.latvec; + const ModuleBase::Vector3 r_index = tmp_FR_container->get_atom_pair(iap).get_R_index(iR); + ModuleBase::Vector3 dR = ModuleBase::Vector3(r_index.x, r_index.y, r_index.z) * GlobalC::ucell.latvec; double phase = ikb_car * dR * ModuleBase::TWO_PI; std::complex kRn_phase = std::exp(ModuleBase::IMAG_UNIT * phase); for(int i = 0; i < row_size; ++i)