diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 53052e87af..68f5890d98 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -158,6 +158,7 @@ - [rpa](#rpa) - [nbands\_istate](#nbands_istate) - [bands\_to\_print](#bands_to_print) + - [if\_separate\_k](#if_separate_k) - [Density of states](#density-of-states) - [dos\_edelta\_ev](#dos_edelta_ev) - [dos\_sigma](#dos_sigma) @@ -449,8 +450,8 @@ These variables are used to control general system parameters. - **relax**: perform structure relaxation calculations, the `relax_nmax` parameter depicts the maximal number of ionic iterations - **cell-relax**: perform cell relaxation calculations - **md**: perform molecular dynamics simulations - - **get_pchg**: obtain partial charge density (for LCAO basis only). See `nbands_istate` and `bands_to_print` for more information - - **get_wf**: obtain wave functions (for LCAO basis only). See `nbands_istate` for more information + - **get_pchg**: obtain partial (band-decomposed) charge densities (for LCAO basis only). See `nbands_istate` and `bands_to_print` for more information + - **get_wf**: obtain wave functions (for LCAO basis only). See `nbands_istate` and `bands_to_print` for more information - **get_S** : obtain the overlap matrix formed by localized orbitals (for LCAO basis with multiple k points). the file name is `SR.csr` with file format being the same as that generated by [out_mat_hs2](#out_mat_hs2) - **gen_bessel** : generates projectors, i.e., a series of Bessel functions, for the DeePKS method (for LCAO basis only); see also keywords `bessel_descriptor_lmax`, `bessel_descriptor_rcut` and `bessel_descriptor_tolerence`. A file named `jle.orb` will be generated which contains the projectors. An example is provided in examples/H2O-deepks-pw - **test_memory** : obtain a rough estimation of memory consuption for the calculation @@ -1755,10 +1756,17 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t ### bands_to_print - **Type**: String -- **Availability**: For both PW and LCAO. When `basis_type = lcao`, only used when `calculation = get_pchg`. -- **Description**: Specifies the bands to calculate the charge density for, using a space-separated string of 0s and 1s, providing a more flexible selection compared to `nbands_istate`. Each digit in the string corresponds to a band, starting from the first band. A `1` indicates that the charge density should be calculated for that band, while a `0` means the band will be ignored. The parameter allows a compact and flexible notation (similar to [`ocp_set`](#ocp_set)), for example the syntax `1 4*0 5*1 0` is used to denote the selection of bands: `1` means calculate for the first band, `4*0` skips the next four bands, `5*1` means calculate for the following five bands, and the final `0` skips the next band. It's essential that the total count of bands does not exceed the total number of bands (`nbands`); otherwise, it results in an error, and the process exits. The input string must contain only numbers and the asterisk (`*`) for repetition, ensuring correct format and intention of band selection. +- **Availability**: For both PW and LCAO. When `basis_type = lcao`, used when `calculation = get_wf` or `calculation = get_pchg`. +- **Description**: Specifies the bands to calculate the wave functions/charge densities for, using a space-separated string of 0s and 1s, providing a more flexible selection compared to `nbands_istate`. Each digit in the string corresponds to a band, starting from the first band. A `1` indicates that the charge density should be calculated for that band, while a `0` means the band will be ignored. The parameter allows a compact and flexible notation (similar to [`ocp_set`](#ocp_set)), for example the syntax `1 4*0 5*1 0` is used to denote the selection of bands: `1` means calculate for the first band, `4*0` skips the next four bands, `5*1` means calculate for the following five bands, and the final `0` skips the next band. It's essential that the total count of bands does not exceed the total number of bands (`nbands`); otherwise, it results in an error, and the process exits. The input string must contain only numbers and the asterisk (`*`) for repetition, ensuring correct format and intention of band selection. - **Default**: none +### if_separate_k + +- **Type**: Boolean +- **Availability**: Only for LCAO, used only when `calculation = get_pchg` and `gamma_only` is turned off. +- **Description**: Specifies whether to write the partial charge densities for all k-points to individual files or merge them. **Warning**: Enabling symmetry may produce incorrect results due to incorrect k-point weights. Therefore, when calculating partial charge densities, it is strongly recommended to set `symmetry = -1`. +- **Default**: false + [back to top](#full-list-of-input-keywords) ## Density of states diff --git a/source/module_base/global_function.h b/source/module_base/global_function.h index 46b8871a17..54d7dd61ec 100644 --- a/source/module_base/global_function.h +++ b/source/module_base/global_function.h @@ -161,6 +161,13 @@ static inline void DCOPY(const T& a, T& b, const int& dim) } } +template +inline void DCOPY(const T* a, T* b, const int& dim) { + for (int i = 0; i < dim; ++i) { + b[i] = a[i]; + } +} + template inline void COPYARRAY(const T* a, T* b, const long dim); diff --git a/source/module_elecstate/module_dm/density_matrix.cpp b/source/module_elecstate/module_dm/density_matrix.cpp index 4378583f03..aff61c9fda 100644 --- a/source/module_elecstate/module_dm/density_matrix.cpp +++ b/source/module_elecstate/module_dm/density_matrix.cpp @@ -1,9 +1,9 @@ #include "density_matrix.h" #include "module_base/libm/libm.h" -#include "module_base/tool_title.h" #include "module_base/memory.h" #include "module_base/timer.h" +#include "module_base/tool_title.h" #include "module_cell/klist.h" namespace elecstate @@ -126,7 +126,7 @@ void DensityMatrix::init_DMR(Grid_Driver* GridD_in, const UnitCell* ucel } } // allocate the memory of BaseMatrix in SR, and set the new values to zero - if(std::is_same::value) + if (std::is_same::value) { tmp_DMR->fix_gamma(); } @@ -171,12 +171,17 @@ void DensityMatrix::init_DMR(Record_adj& ra, const UnitCell* ucell) { continue; } - hamilt::AtomPair tmp_ap(iat1, iat2, ra.info[iat1][ad][0], ra.info[iat1][ad][1], ra.info[iat1][ad][2], this->_paraV); + hamilt::AtomPair tmp_ap(iat1, + iat2, + ra.info[iat1][ad][0], + ra.info[iat1][ad][1], + ra.info[iat1][ad][2], + this->_paraV); tmp_DMR->insert_pair(tmp_ap); } } // allocate the memory of BaseMatrix in SR, and set the new values to zero - if(std::is_same::value) + if (std::is_same::value) { tmp_DMR->fix_gamma(); } @@ -227,15 +232,15 @@ void DensityMatrix::init_DMR(const hamilt::HContainer& DMR_in) this->_DMR.clear(); // set up a HContainer using another one int size_ap = DMR_in.size_atom_pairs(); - if(size_ap > 0) - { + if (size_ap > 0) + { const Parallel_Orbitals* paraV_ = DMR_in.get_atom_pair(0).get_paraV(); hamilt::HContainer* tmp_DMR = new hamilt::HContainer(paraV_); - for(int iap=0;iap R_index = DMR_in.get_atom_pair(iap).get_R_index(ir); hamilt::AtomPair tmp_ap(iat1, iat2, R_index, paraV_); @@ -244,7 +249,7 @@ void DensityMatrix::init_DMR(const hamilt::HContainer& DMR_in) } tmp_DMR->allocate(nullptr, true); this->_DMR.push_back(tmp_DMR); - if(this->_nspin == 2) + if (this->_nspin == 2) { hamilt::HContainer* tmp_DMR1 = new hamilt::HContainer(*tmp_DMR); this->_DMR.push_back(tmp_DMR1); @@ -301,7 +306,8 @@ void DensityMatrix::set_DMK_zero() { for (int ik = 0; ik < _nspin * _nks; ik++) { - ModuleBase::GlobalFunc::ZEROS(this->_DMK[ik].data(), this->_paraV->get_row_size() * this->_paraV->get_col_size()); + ModuleBase::GlobalFunc::ZEROS(this->_DMK[ik].data(), + this->_paraV->get_row_size() * this->_paraV->get_col_size()); } } @@ -359,29 +365,29 @@ void DensityMatrix::save_DMR() { ModuleBase::TITLE("DensityMatrix", "save_DMR"); ModuleBase::timer::tick("DensityMatrix", "save_DMR"); - + const int nnr = this->_DMR[0]->get_nnr(); - // allocate if _DMR_save is empty - if(_DMR_save.size() == 0) + // allocate if _DMR_save is empty + if (_DMR_save.size() == 0) { _DMR_save.resize(this->_DMR.size()); } // resize if _DMR_save[is].size is not equal to _DMR.size - for(int is = 0; is < _DMR_save.size(); is++) + for (int is = 0; is < _DMR_save.size(); is++) { - if(_DMR_save[is].size() != nnr) + if (_DMR_save[is].size() != nnr) { _DMR_save[is].resize(nnr); } } // save _DMR to _DMR_save - for(int is=0;is_DMR.size();is++) + for (int is = 0; is < this->_DMR.size(); is++) { TR* DMR_pointer = this->_DMR[is]->get_wrapper(); TR* DMR_save_pointer = _DMR_save[is].data(); // set to zero ModuleBase::GlobalFunc::ZEROS(DMR_save_pointer, nnr); - for(int i=0;i::save_DMR() // calculate DMR from DMK using add_element template -void DensityMatrix::cal_DMR_test() +void DensityMatrix::cal_DMR_test() { for (int is = 1; is <= this->_nspin; ++is) { - int ik_begin = this->_nks*(is-1); // jump this->_nks for spin_down if nspin==2 - hamilt::HContainer* tmp_DMR = this->_DMR[is-1]; + int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2 + hamilt::HContainer* tmp_DMR = this->_DMR[is - 1]; // set zero since this function is called in every scf step tmp_DMR->set_zero(); #ifdef _OPENMP - #pragma omp parallel for +#pragma omp parallel for #endif for (int iap = 0; iap < tmp_DMR->size_atom_pairs(); ++iap) { @@ -420,11 +426,11 @@ void DensityMatrix::cal_DMR_test() 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) - { - std::cout << "tmp_matrix is nullptr" << std::endl; - continue; - } + if (tmp_matrix == nullptr) + { + std::cout << "tmp_matrix is nullptr" << std::endl; + continue; + } #endif std::complex tmp_res; // loop over k-points @@ -443,7 +449,8 @@ void DensityMatrix::cal_DMR_test() for (int j = 0; j < this->_paraV->get_col_size(iat2); ++j) { // since DMK is column-major, we need to transpose it col=>row - tmp_res = kphase * this->_DMK[ik_begin+ik][(col_ap+j)*this->_paraV->nrow+row_ap+i]; + tmp_res + = kphase * this->_DMK[ik_begin + ik][(col_ap + j) * this->_paraV->nrow + row_ap + i]; tmp_matrix->add_element(i, j, tmp_res.real()); } } @@ -463,7 +470,7 @@ void DensityMatrix, double>::cal_DMR() #ifdef __DEBUG assert(!this->_DMR.empty() && "DMR has not been initialized!"); #endif - + ModuleBase::timer::tick("DensityMatrix", "cal_DMR"); int ld_hk = this->_paraV->nrow; int ld_hk2 = 2 * ld_hk; @@ -474,7 +481,7 @@ void DensityMatrix, double>::cal_DMR() // set zero since this function is called in every scf step tmp_DMR->set_zero(); #ifdef _OPENMP - #pragma omp parallel for +#pragma omp parallel for #endif for (int i = 0; i < tmp_DMR->size_atom_pairs(); ++i) { @@ -500,52 +507,54 @@ void DensityMatrix, double>::cal_DMR() } #endif // loop over k-points - if(GlobalV::NSPIN !=4 ) - { - for (int ik = 0; ik < this->_nks; ++ik) - { - // 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 double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI; - double sinp, cosp; - ModuleBase::libm::sincos(arg, &sinp, &cosp); - std::complex kphase = std::complex(cosp, sinp); - // set DMR element - double* tmp_DMR_pointer = tmp_matrix->get_pointer(); - std::complex* tmp_DMK_pointer = this->_DMK[ik + ik_begin].data(); - double* DMK_real_pointer = nullptr; - double* DMK_imag_pointer = nullptr; - // jump DMK to fill DMR - // DMR is row-major, DMK is column-major - tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap; - for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu) - { - DMK_real_pointer = (double*)tmp_DMK_pointer; - DMK_imag_pointer = DMK_real_pointer + 1; - BlasConnector::axpy(this->_paraV->get_col_size(iat2), - kphase.real(), - DMK_real_pointer, - ld_hk2, - tmp_DMR_pointer, - 1); - // "-" since i^2 = -1 - BlasConnector::axpy(this->_paraV->get_col_size(iat2), - -kphase.imag(), - DMK_imag_pointer, - ld_hk2, - tmp_DMR_pointer, - 1); - tmp_DMK_pointer += 1; - tmp_DMR_pointer += this->_paraV->get_col_size(iat2); - } - } - } + if (GlobalV::NSPIN != 4) + { + for (int ik = 0; ik < this->_nks; ++ik) + { + // 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 double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI; + double sinp, cosp; + ModuleBase::libm::sincos(arg, &sinp, &cosp); + std::complex kphase = std::complex(cosp, sinp); + // set DMR element + double* tmp_DMR_pointer = tmp_matrix->get_pointer(); + std::complex* tmp_DMK_pointer = this->_DMK[ik + ik_begin].data(); + double* DMK_real_pointer = nullptr; + double* DMK_imag_pointer = nullptr; + // jump DMK to fill DMR + // DMR is row-major, DMK is column-major + tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap; + for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu) + { + DMK_real_pointer = (double*)tmp_DMK_pointer; + DMK_imag_pointer = DMK_real_pointer + 1; + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + kphase.real(), + DMK_real_pointer, + ld_hk2, + tmp_DMR_pointer, + 1); + // "-" since i^2 = -1 + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + -kphase.imag(), + DMK_imag_pointer, + ld_hk2, + tmp_DMR_pointer, + 1); + tmp_DMK_pointer += 1; + tmp_DMR_pointer += this->_paraV->get_col_size(iat2); + } + } + } // treat DMR as pauli matrix when NSPIN=4 - if(GlobalV::NSPIN==4) + if (GlobalV::NSPIN == 4) { - std::vector> tmp_DMR(this->_paraV->get_col_size() * this->_paraV->get_row_size(), std::complex(0.0, 0.0)); + std::vector> tmp_DMR(this->_paraV->get_col_size() + * this->_paraV->get_row_size(), + std::complex(0.0, 0.0)); for (int ik = 0; ik < this->_nks; ++ik) { // cal k_phase @@ -583,15 +592,15 @@ void DensityMatrix, double>::cal_DMR() for (int is2 = 0; is2 < npol; is2++) { step_trace[is * npol + is2] = this->_paraV->get_col_size(iat2) * is + is2; - //step_trace[is + is2 * npol] = this->_paraV->get_col_size(iat2) * is + is2; + // step_trace[is + is2 * npol] = this->_paraV->get_col_size(iat2) * is + is2; } } std::complex tmp[4]; double* target_DMR = tmp_matrix->get_pointer(); std::complex* tmp_DMR_pointer = tmp_DMR.data(); - for(int irow=0;irow_paraV->get_row_size(iat1);irow += 2) + for (int irow = 0; irow < this->_paraV->get_row_size(iat1); irow += 2) { - for(int icol=0;icol_paraV->get_col_size(iat2);icol += 2) + for (int icol = 0; icol < this->_paraV->get_col_size(iat2); icol += 2) { // catch the 4 spin component value of one orbital pair tmp[0] = tmp_DMR_pointer[icol + step_trace[0]]; @@ -602,7 +611,8 @@ void DensityMatrix, double>::cal_DMR() // save them back to the tmp_matrix target_DMR[icol + step_trace[0]] = tmp[0].real() + tmp[3].real(); target_DMR[icol + step_trace[1]] = tmp[1].real() + tmp[2].real(); - target_DMR[icol + step_trace[2]] = - tmp[1].imag() + tmp[2].imag();// (i * (rho_updown - rho_downup)).real() + target_DMR[icol + step_trace[2]] + = -tmp[1].imag() + tmp[2].imag(); // (i * (rho_updown - rho_downup)).real() target_DMR[icol + step_trace[3]] = tmp[0].real() - tmp[3].real(); } tmp_DMR_pointer += this->_paraV->get_col_size(iat2) * 2; @@ -615,6 +625,96 @@ void DensityMatrix, double>::cal_DMR() ModuleBase::timer::tick("DensityMatrix", "cal_DMR"); } +// calculate DMR from DMK using blas for multi-k calculation, without summing over k-points +template <> +void DensityMatrix, double>::cal_DMR(const int ik) +{ + ModuleBase::TITLE("DensityMatrix", "cal_DMR"); + + // To check whether DMR has been initialized +#ifdef __DEBUG + assert(!this->_DMR.empty() && "DMR has not been initialized!"); +#endif + + ModuleBase::timer::tick("DensityMatrix", "cal_DMR"); + int ld_hk = this->_paraV->nrow; + int ld_hk2 = 2 * ld_hk; + for (int is = 1; is <= this->_nspin; ++is) + { + int ik_begin = this->_nks * (is - 1); // jump this->_nks for spin_down if nspin==2 + hamilt::HContainer* tmp_DMR = this->_DMR[is - 1]; + // set zero since this function is called in every scf step + tmp_DMR->set_zero(); +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int i = 0; i < tmp_DMR->size_atom_pairs(); ++i) + { + hamilt::AtomPair& tmp_ap = tmp_DMR->get_atom_pair(i); + int iat1 = tmp_ap.get_atom_i(); + int iat2 = tmp_ap.get_atom_j(); + // get global indexes of whole matrix for each atom in this process + int row_ap = this->_paraV->atom_begin_row[iat1]; + int col_ap = this->_paraV->atom_begin_col[iat2]; + if (row_ap == -1 || col_ap == -1) + { + throw std::string("Atom-pair not belong this process"); + } + for (int ir = 0; ir < tmp_ap.get_R_size(); ++ir) + { + 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) + { + std::cout << "tmp_matrix is nullptr" << std::endl; + continue; + } +#endif + // Remove loop over k-points and directly use the provided ik + if (GlobalV::NSPIN != 4) + { + // 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 double arg = (this->_kv->kvec_d[ik] * dR) * ModuleBase::TWO_PI; + double sinp, cosp; + ModuleBase::libm::sincos(arg, &sinp, &cosp); + std::complex kphase = std::complex(cosp, sinp); + // set DMR element + double* tmp_DMR_pointer = tmp_matrix->get_pointer(); + std::complex* tmp_DMK_pointer = this->_DMK[ik + ik_begin].data(); + double* DMK_real_pointer = nullptr; + double* DMK_imag_pointer = nullptr; + // jump DMK to fill DMR + // DMR is row-major, DMK is column-major + tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap; + for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu) + { + DMK_real_pointer = (double*)tmp_DMK_pointer; + DMK_imag_pointer = DMK_real_pointer + 1; + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + kphase.real(), + DMK_real_pointer, + ld_hk2, + tmp_DMR_pointer, + 1); + // "-" since i^2 = -1 + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + -kphase.imag(), + DMK_imag_pointer, + ld_hk2, + tmp_DMR_pointer, + 1); + tmp_DMK_pointer += 1; + tmp_DMR_pointer += this->_paraV->get_col_size(iat2); + } + } + } + } + } + ModuleBase::timer::tick("DensityMatrix", "cal_DMR"); +} // calculate DMR from DMK using blas for gamma-only calculation template <> @@ -635,13 +735,13 @@ void DensityMatrix::cal_DMR() hamilt::HContainer* tmp_DMR = this->_DMR[is - 1]; // set zero since this function is called in every scf step tmp_DMR->set_zero(); - + #ifdef __DEBUG - //assert(tmp_DMR->is_gamma_only() == true); + // assert(tmp_DMR->is_gamma_only() == true); assert(this->_nks == 1); #endif #ifdef _OPENMP - #pragma omp parallel for +#pragma omp parallel for #endif for (int i = 0; i < tmp_DMR->size_atom_pairs(); ++i) { @@ -678,7 +778,12 @@ void DensityMatrix::cal_DMR() tmp_DMK_pointer += col_ap * this->_paraV->nrow + row_ap; for (int mu = 0; mu < this->_paraV->get_row_size(iat1); ++mu) { - BlasConnector::axpy(this->_paraV->get_col_size(iat2), kphase, tmp_DMK_pointer, ld_hk, tmp_DMR_pointer, 1); + BlasConnector::axpy(this->_paraV->get_col_size(iat2), + kphase, + tmp_DMK_pointer, + ld_hk, + tmp_DMR_pointer, + 1); tmp_DMK_pointer += 1; tmp_DMR_pointer += this->_paraV->get_col_size(iat2); } @@ -794,9 +899,10 @@ void DensityMatrix::write_DMK(const std::string directory, const { for (int j = 0; j < this->_paraV->ncol; ++j) { - if (j % 8 == 0) { + if (j % 8 == 0) + { ofs << "\n"; -} + } ofs << " " << this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->ncol + j]; } } @@ -830,9 +936,10 @@ void DensityMatrix, double>::write_DMK(const std::string di { for (int j = 0; j < this->_paraV->ncol; ++j) { - if (j % 8 == 0) { + if (j % 8 == 0) + { ofs << "\n"; -} + } ofs << " " << this->_DMK[ik + this->_nks * (ispin - 1)][i * this->_paraV->ncol + j].real(); } } diff --git a/source/module_elecstate/module_dm/density_matrix.h b/source/module_elecstate/module_dm/density_matrix.h index 90f0db4ae7..0bc28c6cad 100644 --- a/source/module_elecstate/module_dm/density_matrix.h +++ b/source/module_elecstate/module_dm/density_matrix.h @@ -176,6 +176,11 @@ class DensityMatrix */ void cal_DMR(); + /** + * @brief calculate density matrix DMR from dm(k) using blas::axpy for multi-k calculation, without summing over k-points + */ + void cal_DMR(const int ik); + /** * @brief calculate density matrix DMR from dm(k) using base_matrix->add_element() */ diff --git a/source/module_esolver/esolver_ks_lcao_elec.cpp b/source/module_esolver/esolver_ks_lcao_elec.cpp index a77ccb92f0..9eefd6d680 100644 --- a/source/module_esolver/esolver_ks_lcao_elec.cpp +++ b/source/module_esolver/esolver_ks_lcao_elec.cpp @@ -441,30 +441,60 @@ void ESolver_KS_LCAO::others(const int istep) } else if (cal_type == "get_pchg") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting partial charge"); IState_Charge ISC(this->psi, &(this->ParaV)); - ISC.begin(this->GG, - this->pelec->charge->rho, - this->pelec->wg, - this->pelec->eferm.get_all_ef(), - this->pw_rho->nrxx, - this->pw_rho->nplane, - this->pw_rho->startz_current, - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pw_big->bz, - this->pw_big->nbz, - GlobalV::GAMMA_ONLY_LOCAL, - GlobalV::NBANDS_ISTATE, - INPUT.get_out_band_kb(), - GlobalV::NBANDS, - GlobalV::nelec, - GlobalV::NSPIN, - GlobalV::NLOCAL, - GlobalV::global_out_dir, - GlobalV::MY_RANK, - GlobalV::ofs_warning, - &GlobalC::ucell, - &GlobalC::GridD); + if (GlobalV::GAMMA_ONLY_LOCAL) { + ISC.begin(this->GG, + this->pelec->charge->rho, + this->pelec->wg, + this->pelec->eferm.get_all_ef(), + this->pw_rho->nrxx, + this->pw_rho->nplane, + this->pw_rho->startz_current, + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + this->pw_big->bz, + this->pw_big->nbz, + GlobalV::GAMMA_ONLY_LOCAL, + GlobalV::NBANDS_ISTATE, + PARAM.globalv.out_band_kb, + GlobalV::NBANDS, + GlobalV::nelec, + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::global_out_dir, + GlobalV::MY_RANK, + GlobalV::ofs_warning, + &GlobalC::ucell, + &GlobalC::GridD, + this->kv); + } else { + ISC.begin(this->GK, + this->pelec->charge->rho, + this->pelec->wg, + this->pelec->eferm.get_all_ef(), + this->pw_rho->nrxx, + this->pw_rho->nplane, + this->pw_rho->startz_current, + this->pw_rho->nx, + this->pw_rho->ny, + this->pw_rho->nz, + this->pw_big->bz, + this->pw_big->nbz, + GlobalV::GAMMA_ONLY_LOCAL, + GlobalV::NBANDS_ISTATE, + PARAM.globalv.out_band_kb, + GlobalV::NBANDS, + GlobalV::nelec, + GlobalV::NSPIN, + GlobalV::NLOCAL, + GlobalV::global_out_dir, + GlobalV::MY_RANK, + GlobalV::ofs_warning, + &GlobalC::ucell, + &GlobalC::GridD, + this->kv, + PARAM.inp.if_separate_k); + } std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "getting partial charge"); } else if (cal_type == "get_wf") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting wave function"); diff --git a/source/module_io/bcast_globalv.cpp b/source/module_io/bcast_globalv.cpp index 88c90e18a9..b7f522e238 100644 --- a/source/module_io/bcast_globalv.cpp +++ b/source/module_io/bcast_globalv.cpp @@ -18,5 +18,6 @@ void ReadInput::set_globalv_bcast() add_bool_bcast(sys.double_grid); add_double_bcast(sys.uramping); add_string_bcast(sys.global_calculation); + add_intvec_bcast(sys.out_band_kb, para.sys.out_band_kb.size(), 0); } } // namespace ModuleIO \ No newline at end of file diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 11a28193bf..f23e789efc 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -1,7 +1,5 @@ #include "module_io/input_conv.h" -#include - #include "module_base/global_function.h" #include "module_base/global_variable.h" #include "module_cell/module_symmetry/symmetry.h" @@ -15,6 +13,8 @@ #include "module_relax/relax_old/ions_move_basic.h" #include "module_relax/relax_old/lattice_change_basic.h" +#include + #ifdef __EXX #include "module_ri/exx_abfs-jle.h" #endif @@ -43,91 +43,6 @@ #include "module_hsolver/hsolver_pw.h" #include "module_md/md_func.h" -template -void Input_Conv::parse_expression(const std::string& fn, std::vector& vec) -{ - ModuleBase::TITLE("Input_Conv", "parse_expression"); - int count = 0; - std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)"); - std::vector str; - std::stringstream ss(fn); - std::string section; - while (ss >> section) - { - int index = 0; - if (str.empty()) - { - while (index < section.size() && std::isspace(section[index])) - { - index++; - } - } - section.erase(0, index); - str.push_back(section); - } - // std::string::size_type pos1, pos2; - // std::string c = " "; - // pos2 = fn.find(c); - // pos1 = 0; - // while (std::string::npos != pos2) - // { - // str.push_back(fn.substr(pos1, pos2 - pos1)); - // pos1 = pos2 + c.size(); - // pos2 = fn.find(c, pos1); - // } - // if (pos1 != fn.length()) - // { - // str.push_back(fn.substr(pos1)); - // } - regex_t reg; - regcomp(®, pattern.c_str(), REG_EXTENDED); - regmatch_t pmatch[1]; - const size_t nmatch = 1; - for (size_t i = 0; i < str.size(); ++i) - { - if (str[i] == "") - { - continue; - } - int status = regexec(®, str[i].c_str(), nmatch, pmatch, 0); - std::string sub_str = ""; - for (size_t j = pmatch[0].rm_so; j != pmatch[0].rm_eo; ++j) - { - sub_str += str[i][j]; - } - std::string sub_pattern("\\*"); - regex_t sub_reg; - regcomp(&sub_reg, sub_pattern.c_str(), REG_EXTENDED); - regmatch_t sub_pmatch[1]; - const size_t sub_nmatch = 1; - if (regexec(&sub_reg, sub_str.c_str(), sub_nmatch, sub_pmatch, 0) == 0) - { - int pos = sub_str.find("*"); - int num = stoi(sub_str.substr(0, pos)); - T occ = stof(sub_str.substr(pos + 1, sub_str.size())); - // std::vector ocp_temp(num, occ); - // const std::vector::iterator dest = vec.begin() + count; - // copy(ocp_temp.begin(), ocp_temp.end(), dest); - // count += num; - for (size_t k = 0; k != num; k++) { - vec.emplace_back(occ); -} - } - else - { - // vec[count] = stof(sub_str); - // count += 1; - std::stringstream convert; - convert << sub_str; - T occ; - convert >> occ; - vec.emplace_back(occ); - } - regfree(&sub_reg); - } - regfree(®); -} - #ifdef __LCAO std::vector Input_Conv::convert_units(std::string params, double c) { @@ -313,7 +228,11 @@ void Input_Conv::Convert() ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "pseudo_dir", GlobalV::global_pseudo_dir); ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "orbital_dir", GlobalV::global_orbital_dir); // GlobalV::global_pseudo_type = PARAM.inp.pseudo_type; - GlobalC::ucell.setup(PARAM.inp.latname, PARAM.inp.ntype, PARAM.inp.lmaxmax, PARAM.inp.init_vel, PARAM.inp.fixed_axes); + GlobalC::ucell.setup(PARAM.inp.latname, + PARAM.inp.ntype, + PARAM.inp.lmaxmax, + PARAM.inp.init_vel, + PARAM.inp.fixed_axes); if (PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax") { @@ -564,7 +483,10 @@ void Input_Conv::Convert() if (PARAM.inp.restart_save) { std::string dft_functional_lower = PARAM.inp.dft_functional; - std::transform(PARAM.inp.dft_functional.begin(), PARAM.inp.dft_functional.end(), dft_functional_lower.begin(), tolower); + std::transform(PARAM.inp.dft_functional.begin(), + PARAM.inp.dft_functional.end(), + dft_functional_lower.begin(), + tolower); GlobalC::restart.folder = GlobalV::global_readin_dir + "restart/"; ModuleBase::GlobalFunc::MAKE_DIR(GlobalC::restart.folder); if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "hse" @@ -581,7 +503,10 @@ void Input_Conv::Convert() if (PARAM.inp.restart_load) { std::string dft_functional_lower = PARAM.inp.dft_functional; - std::transform(PARAM.inp.dft_functional.begin(), PARAM.inp.dft_functional.end(), dft_functional_lower.begin(), tolower); + std::transform(PARAM.inp.dft_functional.begin(), + PARAM.inp.dft_functional.end(), + dft_functional_lower.begin(), + tolower); GlobalC::restart.folder = GlobalV::global_readin_dir + "restart/"; if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "hse" || dft_functional_lower == "opt_orb" || dft_functional_lower == "scan0") @@ -602,7 +527,10 @@ void Input_Conv::Convert() #ifdef __LCAO std::string dft_functional_lower = PARAM.inp.dft_functional; - std::transform(PARAM.inp.dft_functional.begin(), PARAM.inp.dft_functional.end(), dft_functional_lower.begin(), tolower); + std::transform(PARAM.inp.dft_functional.begin(), + PARAM.inp.dft_functional.end(), + dft_functional_lower.begin(), + tolower); if (dft_functional_lower == "hf" || dft_functional_lower == "pbe0" || dft_functional_lower == "scan0") { GlobalC::exx_info.info_global.cal_exx = true; @@ -653,10 +581,12 @@ void Input_Conv::Convert() // EXX does not support symmetry=1 if (PARAM.inp.calculation != "nscf" && PARAM.inp.symmetry == "1") + { ModuleSymmetry::Symmetry::symm_flag = 0; + } } -#endif // __LCAO -#endif // __EXX +#endif // __LCAO +#endif // __EXX GlobalC::ppcell.cell_factor = PARAM.inp.cell_factor; // LiuXh add 20180619 //---------------------------------------------------------- @@ -778,19 +708,23 @@ void Input_Conv::Convert() { GlobalV::deepks_out_labels = true; GlobalV::deepks_scf = true; - if (GlobalV::NPROC > 1) { + if (GlobalV::NPROC > 1) + { ModuleBase::WARNING_QUIT("Input_conv", "generate deepks unittest with only 1 processor"); -} - if (GlobalV::CAL_FORCE != 1) { + } + if (GlobalV::CAL_FORCE != 1) + { ModuleBase::WARNING_QUIT("Input_conv", "force is required in generating deepks unittest"); -} - if (GlobalV::CAL_STRESS != 1) { + } + if (GlobalV::CAL_STRESS != 1) + { ModuleBase::WARNING_QUIT("Input_conv", "stress is required in generating deepks unittest"); -} + } } - if (GlobalV::deepks_scf || GlobalV::deepks_out_labels) { + if (GlobalV::deepks_scf || GlobalV::deepks_out_labels) + { GlobalV::deepks_setorb = true; -} + } #else if (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_labels || PARAM.inp.deepks_bandgap || PARAM.inp.deepks_v_delta) { diff --git a/source/module_io/input_conv.h b/source/module_io/input_conv.h index 1aadcb55c3..45226ca64f 100644 --- a/source/module_io/input_conv.h +++ b/source/module_io/input_conv.h @@ -5,13 +5,15 @@ #ifndef INPUT_CONVERT_H #define INPUT_CONVERT_H -#include -#include -#include +#include "module_base/global_function.h" +#include "module_base/global_variable.h" #include #include #include +#include +#include +#include #include #include @@ -35,12 +37,95 @@ void Convert(void); * * @tparam T * @param fn (string): expressions such as "3*1 0 2*0.5 3*0" - * @param arr (vector): stores parsing results, + * @param vec (vector): stores parsing results, * for example, "3*1 0 2*0.5 1*1.5" can be parsed as * [1, 1, 1, 0, 0.5, 0.5, 1.5] */ template -void parse_expression(const std::string& fn, std::vector& arr); +void parse_expression(const std::string& fn, std::vector& vec) +{ + ModuleBase::TITLE("Input_Conv", "parse_expression"); + int count = 0; + std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)"); + std::vector str; + std::stringstream ss(fn); + std::string section; + while (ss >> section) + { + int index = 0; + if (str.empty()) + { + while (index < section.size() && std::isspace(section[index])) + { + index++; + } + } + section.erase(0, index); + str.push_back(section); + } + // std::string::size_type pos1, pos2; + // std::string c = " "; + // pos2 = fn.find(c); + // pos1 = 0; + // while (std::string::npos != pos2) + // { + // str.push_back(fn.substr(pos1, pos2 - pos1)); + // pos1 = pos2 + c.size(); + // pos2 = fn.find(c, pos1); + // } + // if (pos1 != fn.length()) + // { + // str.push_back(fn.substr(pos1)); + // } + regex_t reg; + regcomp(®, pattern.c_str(), REG_EXTENDED); + regmatch_t pmatch[1]; + const size_t nmatch = 1; + for (size_t i = 0; i < str.size(); ++i) + { + if (str[i] == "") + { + continue; + } + int status = regexec(®, str[i].c_str(), nmatch, pmatch, 0); + std::string sub_str = ""; + for (size_t j = pmatch[0].rm_so; j != pmatch[0].rm_eo; ++j) + { + sub_str += str[i][j]; + } + std::string sub_pattern("\\*"); + regex_t sub_reg; + regcomp(&sub_reg, sub_pattern.c_str(), REG_EXTENDED); + regmatch_t sub_pmatch[1]; + const size_t sub_nmatch = 1; + if (regexec(&sub_reg, sub_str.c_str(), sub_nmatch, sub_pmatch, 0) == 0) + { + int pos = sub_str.find("*"); + int num = stoi(sub_str.substr(0, pos)); + T occ = stof(sub_str.substr(pos + 1, sub_str.size())); + // std::vector ocp_temp(num, occ); + // const std::vector::iterator dest = vec.begin() + count; + // copy(ocp_temp.begin(), ocp_temp.end(), dest); + // count += num; + for (size_t k = 0; k != num; k++) + { + vec.emplace_back(occ); + } + } + else + { + // vec[count] = stof(sub_str); + // count += 1; + std::stringstream convert; + convert << sub_str; + T occ; + convert >> occ; + vec.emplace_back(occ); + } + regfree(&sub_reg); + } + regfree(®); +} #ifdef __LCAO /** diff --git a/source/module_io/istate_charge.cpp b/source/module_io/istate_charge.cpp index f12cd9e903..87116b8fa7 100644 --- a/source/module_io/istate_charge.cpp +++ b/source/module_io/istate_charge.cpp @@ -16,6 +16,11 @@ IState_Charge::IState_Charge(psi::Psi* psi_gamma_in, const Parallel_Orbi { } +IState_Charge::IState_Charge(psi::Psi>* psi_k_in, const Parallel_Orbitals* ParaV_in) + : psi_k(psi_k_in), ParaV(ParaV_in) +{ +} + IState_Charge::~IState_Charge() { } @@ -43,17 +48,144 @@ void IState_Charge::begin(Gint_Gamma& gg, const int my_rank, std::ofstream& ofs_warning, const UnitCell* ucell_in, - Grid_Driver* GridD_in) + Grid_Driver* GridD_in, + const K_Vectors& kv) { ModuleBase::TITLE("IState_Charge", "begin"); - std::cout << " Perform |psi(i)|^2 for selected bands." << std::endl; + std::cout << " Perform |psi(i)|^2 for selected bands (band-decomposed charge densities, gamma only)." << std::endl; - if (!gamma_only_local) + int mode = 0; + if (nbands_istate > 0 && static_cast(out_band_kb.size()) == 0) + { + mode = 1; + } + else if (static_cast(out_band_kb.size()) > 0) { - ModuleBase::WARNING_QUIT("IState_Charge::begin", "Only available for GAMMA_ONLY_LOCAL now."); + // If out_band_kb (bands_to_print) is not empty, set mode to 2 + mode = 2; + std::cout << " Notice: INPUT parameter `nbands_istate` overwritten by `bands_to_print`!" << std::endl; + } + else + { + mode = 3; + } + + // if ucell is odd, it's correct, + // if ucell is even, it's also correct. + // +1.0e-8 in case like (2.999999999+1)/2 + const int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + std::cout << " number of electrons = " << nelec << std::endl; + std::cout << " number of occupied bands = " << fermi_band << std::endl; + + // Set this->bands_picked_ according to the mode + select_bands(nbands_istate, out_band_kb, nbands, nelec, mode, fermi_band); + + for (int ib = 0; ib < nbands; ++ib) + { + if (bands_picked_[ib]) + { + // Using new density matrix inplementation (gamma only) + elecstate::DensityMatrix DM(this->ParaV, nspin); + +#ifdef __MPI + this->idmatrix(ib, nspin, nelec, nlocal, wg, DM, kv); +#else + ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); +#endif + + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); + } + + std::cout << " Performing grid integral over real space grid for band " << ib + 1 << "..." << std::endl; + + DM.init_DMR(GridD_in, ucell_in); + DM.cal_DMR(); + gg.initialize_pvpR(*ucell_in, GridD_in); + gg.transfer_DM2DtoGrid(DM.get_DMR_vector()); + Gint_inout inout((double***)nullptr, rho, Gint_Tools::job_type::rho); + gg.cal_gint(&inout); + + // A solution to replace the original implementation of the following code: + // pelec->charge->save_rho_before_sum_band(); + // Using std::vector to replace the original double** rho_save + std::vector> rho_save(nspin, std::vector(rhopw_nrxx)); + + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data + } + + std::cout << " Writting cube files..."; + + for (int is = 0; is < nspin; ++is) + { + // ssc should be inside the inner loop to reset the string stream each time + std::stringstream ssc; + ssc << global_out_dir << "BAND" << ib + 1 << "_GAMMA" << "_SPIN" << is + 1 << "_CHG.cube"; + + // Use a const vector to store efermi for all spins, replace the original implementation: + // const double ef_tmp = pelec->eferm.get_efval(is); + double ef_spin = ef_all_spin[is]; + ModuleIO::write_rho( +#ifdef __MPI + bigpw_bz, + bigpw_nbz, + rhopw_nplane, + rhopw_startz_current, +#endif + rho_save[is].data(), + is, + nspin, + 0, + ssc.str(), + rhopw_nx, + rhopw_ny, + rhopw_nz, + ef_spin, + ucell_in); + } + + std::cout << " Complete!" << std::endl; + } } + return; +} + +void IState_Charge::begin(Gint_k& gk, + double** rho, + const ModuleBase::matrix& wg, + const std::vector& ef_all_spin, + const int rhopw_nrxx, + const int rhopw_nplane, + const int rhopw_startz_current, + const int rhopw_nx, + const int rhopw_ny, + const int rhopw_nz, + const int bigpw_bz, + const int bigpw_nbz, + const bool gamma_only_local, + const int nbands_istate, + const std::vector& out_band_kb, + const int nbands, + const double nelec, + const int nspin, + const int nlocal, + const std::string& global_out_dir, + const int my_rank, + std::ofstream& ofs_warning, + const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const K_Vectors& kv, + const bool if_separate_k) +{ + ModuleBase::TITLE("IState_Charge", "begin"); + + std::cout << " Perform |psi(i)|^2 for selected bands (band-decomposed charge densities, multi-k)." << std::endl; + int mode = 0; if (nbands_istate > 0 && static_cast(out_band_kb.size()) == 0) { @@ -70,36 +202,166 @@ void IState_Charge::begin(Gint_Gamma& gg, mode = 3; } - int fermi_band = 0; + const int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + std::cout << " number of electrons = " << nelec << std::endl; + std::cout << " number of occupied bands = " << fermi_band << std::endl; + + // Set this->bands_picked_ according to the mode + select_bands(nbands_istate, out_band_kb, nbands, nelec, mode, fermi_band); + + for (int ib = 0; ib < nbands; ++ib) + { + if (bands_picked_[ib]) + { + // Using new density matrix inplementation (multi-k) + elecstate::DensityMatrix, double> DM(&kv, this->ParaV, nspin); + +#ifdef __MPI + this->idmatrix(ib, nspin, nelec, nlocal, wg, DM, kv); +#else + ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); +#endif + if (if_separate_k) + { + // For multi-k, loop over all real k-points + for (int ik = 0; ik < kv.get_nks() / nspin; ++ik) + { + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); + } + + std::cout << " Performing grid integral over real space grid for band " << ib + 1 << ", k-point " + << ik + 1 << "..." << std::endl; + + DM.init_DMR(GridD_in, ucell_in); + DM.cal_DMR(ik); + gk.initialize_pvpR(*ucell_in, GridD_in); + gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); + Gint_inout inout(rho, Gint_Tools::job_type::rho); + gk.cal_gint(&inout); + + // Using std::vector to replace the original double** rho_save + std::vector> rho_save(nspin, std::vector(rhopw_nrxx)); + + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data + } + + std::cout << " Writting cube files..."; + + for (int is = 0; is < nspin; ++is) + { + // ssc should be inside the inner loop to reset the string stream each time + std::stringstream ssc; + ssc << global_out_dir << "BAND" << ib + 1 << "_K" << ik + 1 << "_SPIN" << is + 1 << "_CHG.cube"; + + double ef_spin = ef_all_spin[is]; + ModuleIO::write_rho( +#ifdef __MPI + bigpw_bz, + bigpw_nbz, + rhopw_nplane, + rhopw_startz_current, +#endif + rho_save[is].data(), + is, + nspin, + 0, + ssc.str(), + rhopw_nx, + rhopw_ny, + rhopw_nz, + ef_spin, + ucell_in); + } + + std::cout << " Complete!" << std::endl; + } + } + else + { + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); + } + + std::cout << " Performing grid integral over real space grid for band " << ib + 1 << "..." << std::endl; + + DM.init_DMR(GridD_in, ucell_in); + DM.cal_DMR(); + gk.initialize_pvpR(*ucell_in, GridD_in); + gk.transfer_DM2DtoGrid(DM.get_DMR_vector()); + Gint_inout inout(rho, Gint_Tools::job_type::rho); + gk.cal_gint(&inout); + + // Using std::vector to replace the original double** rho_save + std::vector> rho_save(nspin, std::vector(rhopw_nrxx)); + + for (int is = 0; is < nspin; ++is) + { + ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data + } + + std::cout << " Writting cube files..."; + + for (int is = 0; is < nspin; ++is) + { + // ssc should be inside the inner loop to reset the string stream each time + std::stringstream ssc; + ssc << global_out_dir << "BAND" << ib + 1 << "_SPIN" << is + 1 << "_CHG.cube"; + + double ef_spin = ef_all_spin[is]; + ModuleIO::write_rho( +#ifdef __MPI + bigpw_bz, + bigpw_nbz, + rhopw_nplane, + rhopw_startz_current, +#endif + rho_save[is].data(), + is, + nspin, + 0, + ssc.str(), + rhopw_nx, + rhopw_ny, + rhopw_nz, + ef_spin, + ucell_in); + } + + std::cout << " Complete!" << std::endl; + } + } + } + + return; +} + +void IState_Charge::select_bands(const int nbands_istate, + const std::vector& out_band_kb, + const int nbands, + const double nelec, + const int mode, + const int fermi_band) +{ int bands_below = 0; int bands_above = 0; - // (2) cicle: - // (2.1) calculate the selected density matrix from wave functions. - // (2.2) carry out the grid integration to get the charge density. this->bands_picked_.resize(nbands); ModuleBase::GlobalFunc::ZEROS(bands_picked_.data(), nbands); - // (1) - // (1.2) read in WFC_NAO_GAMMA1.dat - std::cout << " number of electrons = " << nelec << std::endl; - - // mohan update 2011-03-21 - // if ucell is odd, it's correct, - // if ucell is even, it's also correct. - // +1.0e-8 in case like (2.999999999+1)/2 - fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); - std::cout << " number of occupied bands = " << fermi_band << std::endl; - if (mode == 1) { bands_below = nbands_istate; bands_above = nbands_istate; - std::cout << " Plot band decomposed charge density below Fermi surface with " << bands_below << " bands." + std::cout << " Plot band-decomposed charge densities below Fermi surface with " << bands_below << " bands." << std::endl; - std::cout << " Plot band decomposed charge density above Fermi surface with " << bands_above << " bands." + std::cout << " Plot band-decomposed charge densities above Fermi surface with " << bands_above << " bands." << std::endl; for (int ib = 0; ib < nbands; ++ib) @@ -119,166 +381,78 @@ void IState_Charge::begin(Gint_Gamma& gg, if (static_cast(out_band_kb.size()) > nbands) { ModuleBase::WARNING_QUIT( - "IState_Charge::begin", + "IState_Charge::select_bands", "The number of bands specified by `bands_to_print` in the INPUT file exceeds `nbands`!"); } - // Check if all elements in bands_picked_ are 0 or 1 + // Check if all elements in out_band_kb are 0 or 1 for (int value: out_band_kb) { if (value != 0 && value != 1) { ModuleBase::WARNING_QUIT( - "IState_Charge::begin", + "IState_Charge::select_bands", "The elements of `bands_to_print` must be either 0 or 1. Invalid values found!"); } } // Fill bands_picked_ with values from out_band_kb // Remaining bands are already set to 0 - int length = std::min(static_cast(out_band_kb.size()), nbands); + const int length = std::min(static_cast(out_band_kb.size()), nbands); for (int i = 0; i < length; ++i) { // out_band_kb rely on function parse_expression from input_conv.cpp bands_picked_[i] = out_band_kb[i]; } - std::cout << " Plot band decomposed charge density below the Fermi surface: band "; + // Check if there are selected bands below the Fermi surface + bool has_below = false; for (int i = 0; i + 1 <= fermi_band; ++i) { if (bands_picked_[i] == 1) { - std::cout << i + 1 << " "; + has_below = true; + break; } } - std::cout << std::endl; - std::cout << " Plot band decomposed charge density above the Fermi surface: band "; - for (int i = fermi_band; i < nbands; ++i) + if (has_below) { - if (bands_picked_[i] == 1) + std::cout << " Plot band-decomposed charge densities below the Fermi surface: band "; + for (int i = 0; i + 1 <= fermi_band; ++i) { - std::cout << i + 1 << " "; + if (bands_picked_[i] == 1) + { + std::cout << i + 1 << " "; + } } + std::cout << std::endl; } - std::cout << std::endl; - } - else if (mode == 3) - { - bool stop = false; - std::stringstream ss; - ss << global_out_dir << "istate.info"; - std::cout << " Open the file : " << ss.str() << std::endl; - if (my_rank == 0) + + // Check if there are selected bands above the Fermi surface + bool has_above = false; + for (int i = fermi_band; i < nbands; ++i) { - std::ifstream ifs(ss.str().c_str()); - if (!ifs) + if (bands_picked_[i] == 1) { - stop = true; + has_above = true; + break; } - else + } + if (has_above) + { + std::cout << " Plot band-decomposed charge densities above the Fermi surface: band "; + for (int i = fermi_band; i < nbands; ++i) { - // int band_index; - for (int ib = 0; ib < nbands; ++ib) + if (bands_picked_[i] == 1) { - ModuleBase::GlobalFunc::READ_VALUE(ifs, bands_picked_[ib]); + std::cout << i + 1 << " "; } } - } - -#ifdef __MPI - Parallel_Common::bcast_bool(stop); - Parallel_Common::bcast_int(bands_picked_.data(), nbands); -#endif - if (stop) - { - ofs_warning << " Can't find the file : " << ss.str() << std::endl; - ModuleBase::WARNING_QUIT("IState_Charge::begin", "can't find the istate file."); + std::cout << std::endl; } } - - for (int ib = 0; ib < nbands; ++ib) + else { - if (bands_picked_[ib]) - { - std::cout << " Perform band decomposed charge density for band " << ib + 1 << std::endl; - - // (1) calculate the density matrix for a partuclar band, whenever it is occupied or not. - - // Using new density matrix inplementation - elecstate::DensityMatrix DM(this->ParaV, nspin); - -#ifdef __MPI - this->idmatrix(ib, nspin, nelec, nlocal, wg, DM); -#else - ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!"); -#endif - - // (2) zero out of charge density array. - for (int is = 0; is < nspin; ++is) - { - ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx); - } - - // (3) calculate charge density for a particular band. - - DM.init_DMR(GridD_in, ucell_in); - DM.cal_DMR(); - - // gg.DMRGint.resize(nspin); - gg.initialize_pvpR(*ucell_in, GridD_in); - - gg.transfer_DM2DtoGrid(DM.get_DMR_vector()); - - Gint_inout inout((double***)nullptr, rho, Gint_Tools::job_type::rho); - gg.cal_gint(&inout); - - // A solution to replace the original implementation of the following code: - // pelec->charge->save_rho_before_sum_band(); - double** rho_save = new double*[nspin]; // Initialize an array of pointers - for (int is = 0; is < nspin; is++) - { - rho_save[is] = new double[rhopw_nrxx]; // Allocate memory for each internal array - ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], - rhopw_nrxx); // Copy data after allocation - } - - std::stringstream ssc; - ssc << global_out_dir << "BAND" << ib + 1; - // 0 means definitely output charge density. - for (int is = 0; is < nspin; ++is) - { - ssc << "_SPIN" << is << "_CHG.cube"; - - // Use a const vector to store efermi for all spins, replace the original implementation: - // const double ef_tmp = pelec->eferm.get_efval(is); - double ef_spin = ef_all_spin[is]; - ModuleIO::write_rho( -#ifdef __MPI - bigpw_bz, - bigpw_nbz, - rhopw_nplane, - rhopw_startz_current, -#endif - rho_save[is], - is, - nspin, - 0, - ssc.str(), - rhopw_nx, - rhopw_ny, - rhopw_nz, - ef_spin, - ucell_in); - } - - // Release memory of rho_save - for (int is = 0; is < nspin; is++) - { - delete[] rho_save[is]; // Release memory of each internal array - } - delete[] rho_save; // Release memory of the array of pointers - } + ModuleBase::WARNING_QUIT("IState_Charge::select_bands", "Invalid mode! Please check the code."); } - - return; } #ifdef __MPI @@ -287,15 +461,18 @@ void IState_Charge::idmatrix(const int& ib, const double& nelec, const int nlocal, const ModuleBase::matrix& wg, - elecstate::DensityMatrix& DM) + elecstate::DensityMatrix& DM, + const K_Vectors& kv) { ModuleBase::TITLE("IState_Charge", "idmatrix"); assert(wg.nr == nspin); - int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + const int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); for (int is = 0; is < nspin; ++is) { + std::cout << " Calculating density matrix for band " << ib + 1 << ", spin " << is + 1 << std::endl; + std::vector wg_local(this->ParaV->ncol, 0.0); const int ib_local = this->ParaV->global2local_col(ib); @@ -313,10 +490,51 @@ void IState_Charge::idmatrix(const int& ib, { BlasConnector::scal(wg_wfc.get_nbasis(), wg_local[ir], wg_wfc.get_pointer() + ir * wg_wfc.get_nbasis(), 1); } - const int ik = 0; // Gamma point only elecstate::psiMulPsiMpi(wg_wfc, *(this->psi_gamma), + DM.get_DMK_pointer(is), + this->ParaV->desc_wfc, + this->ParaV->desc); + } +} + +void IState_Charge::idmatrix(const int& ib, + const int nspin, + const double& nelec, + const int nlocal, + const ModuleBase::matrix& wg, + elecstate::DensityMatrix, double>& DM, + const K_Vectors& kv) +{ + ModuleBase::TITLE("IState_Charge", "idmatrix"); + assert(wg.nr == kv.get_nks()); + + const int fermi_band = static_cast((nelec + 1) / 2 + 1.0e-8); + + for (int ik = 0; ik < kv.get_nks(); ++ik) + { + std::cout << " Calculating density matrix for band " << ib + 1 << ", k-point " + << ik % (kv.get_nks() / nspin) + 1 << ", spin " << kv.isk[ik] + 1 << std::endl; + + std::vector wg_local(this->ParaV->ncol, 0.0); + const int ib_local = this->ParaV->global2local_col(ib); + + if (ib_local >= 0) + { + wg_local[ib_local] = (ib < fermi_band) ? wg(ik, ib) : wg(ik, fermi_band - 1); + } + + this->psi_k->fix_k(ik); + psi::Psi> wg_wfc(*this->psi_k, 1); + + for (int ir = 0; ir < wg_wfc.get_nbands(); ++ir) + { + BlasConnector::scal(wg_wfc.get_nbasis(), wg_local[ir], wg_wfc.get_pointer() + ir * wg_wfc.get_nbasis(), 1); + } + + elecstate::psiMulPsiMpi(wg_wfc, + *(this->psi_k), DM.get_DMK_pointer(ik), this->ParaV->desc_wfc, this->ParaV->desc); diff --git a/source/module_io/istate_charge.h b/source/module_io/istate_charge.h index 7d49661ae7..78c3529a4c 100644 --- a/source/module_io/istate_charge.h +++ b/source/module_io/istate_charge.h @@ -4,6 +4,7 @@ #include "module_elecstate/module_dm/density_matrix.h" #include "module_hamilt_lcao/module_gint/gint.h" #include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_psi/psi.h" #include @@ -12,24 +13,21 @@ #include /** - * @brief Manages the computation of the charge density for different bands. + * @brief Manages the computation of the charge densities for different bands (band-decomposed charge densities). * * This class is responsible for initializing and managing the * charge state computation process, offering functionality to - * calculate and plot the decomposed charge density below and above - * the Fermi surface based on specified bands. + * calculate and plot the decomposed charge density for specified bands. */ class IState_Charge { public: IState_Charge(psi::Psi* psi_gamma_in, const Parallel_Orbitals* ParaV_in); - IState_Charge(psi::Psi>* psi_k_in, const Parallel_Orbitals* ParaV_in) - { - throw std::logic_error("IState_Charge for multi-k is not implemented."); - }; + IState_Charge(psi::Psi>* psi_k_in, const Parallel_Orbitals* ParaV_in); ~IState_Charge(); + // for gamma only void begin(Gint_Gamma& gg, double** rho, const ModuleBase::matrix& wg, @@ -53,34 +51,91 @@ class IState_Charge const int my_rank, std::ofstream& ofs_warning, const UnitCell* ucell_in, - Grid_Driver* GridD_in); + Grid_Driver* GridD_in, + const K_Vectors& kv); + + // For multi-k + void begin(Gint_k& gk, + double** rho, + const ModuleBase::matrix& wg, + const std::vector& ef_all_spin, + const int rhopw_nrxx, + const int rhopw_nplane, + const int rhopw_startz_current, + const int rhopw_nx, + const int rhopw_ny, + const int rhopw_nz, + const int bigpw_bz, + const int bigpw_nbz, + const bool gamma_only_local, + const int nbands_istate, + const std::vector& out_band_kb, + const int nbands, + const double nelec, + const int nspin, + const int nlocal, + const std::string& global_out_dir, + const int my_rank, + std::ofstream& ofs_warning, + const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const K_Vectors& kv, + const bool if_separate_k); private: - std::vector bands_picked_; + /** + * @brief Set this->bands_picked_ according to the mode, and process an error if the mode is not recognized. + * + * @param nbands_istate INPUT parameter nbands_istate. + * @param out_band_kb Calculated from INPUT parameter bands_to_print, vector. + * @param nbands INPUT parameter nbands. + * @param nelec Total number of electrons. + * @param mode Selected mode. + * @param fermi_band Calculated Fermi band. + */ + void select_bands(const int nbands_istate, + const std::vector& out_band_kb, + const int nbands, + const double nelec, + const int mode, + const int fermi_band); #ifdef __MPI /** - * @brief Calculates the density matrix for a given band and spin. + * @brief Calculates the density matrix for a given band. * - * This method calculates the density matrix for a given band and spin using the wave function coefficients. - * It adjusts the coefficients based on the Fermi energy and performs a matrix multiplication to produce the density - * matrix. + * This method calculates the density matrix for a given band using the wave function coefficients. + * It performs a matrix multiplication to produce the density matrix. * * @param ib Band index. * @param nspin Number of spin channels. * @param nelec Total number of electrons. * @param nlocal Number of local orbitals. - * @param wg Weight matrix for bands and spins. + * @param wg Weight matrix for bands and spins (k-points). * @param DM Density matrix to be calculated. + * @param kv K-vectors. */ void idmatrix(const int& ib, const int nspin, const double& nelec, const int nlocal, const ModuleBase::matrix& wg, - elecstate::DensityMatrix& DM); + elecstate::DensityMatrix& DM, + const K_Vectors& kv); + + // For multi-k + void idmatrix(const int& ib, + const int nspin, + const double& nelec, + const int nlocal, + const ModuleBase::matrix& wg, + elecstate::DensityMatrix, double>& DM, + const K_Vectors& kv); + #endif + std::vector bands_picked_; psi::Psi* psi_gamma; + psi::Psi>* psi_k; const Parallel_Orbitals* ParaV; }; #endif diff --git a/source/module_io/read_input_item_general.cpp b/source/module_io/read_input_item_general.cpp index fd7bf3ab5d..ced6c33c22 100644 --- a/source/module_io/read_input_item_general.cpp +++ b/source/module_io/read_input_item_general.cpp @@ -1,12 +1,12 @@ -#include - -#include - +#include "input_conv.h" #include "module_base/global_function.h" #include "module_base/tool_quit.h" #include "read_input.h" #include "read_input_tool.h" +#include +#include + namespace ModuleIO { // There are some examples: @@ -173,8 +173,8 @@ void ReadInput::item_general() if (!find_str(esolver_types, para.input.esolver_type)) { ModuleBase::WARNING_QUIT("ReadInput", - "esolver_type should be ksdft, sdft, " - "ofdft, tddft, lr, ks-lr, lj or dp."); + "esolver_type should be ksdft, sdft, " + "ofdft, tddft, lr, ks-lr, lj or dp."); } if (para.input.esolver_type == "dp") { @@ -276,19 +276,26 @@ void ReadInput::item_general() } { Input_Item item("nbands_istate"); - item.annotation = "number of bands around Fermi level for get_pchg calulation"; + item.annotation = "number of bands around Fermi level for get_wf and get_pchg calulation"; read_sync_int(input.nbands_istate); this->add_item(item); } { Input_Item item("bands_to_print"); - item.annotation = "specify the bands to be calculated in the get_pchg calculation"; + item.annotation = "specify the bands to be calculated in get_wf and get_pchg calculation"; item.read_value = [](const Input_Item& item, Parameter& para) { para.input.bands_to_print = longstring(item.str_values, item.get_size()); + Input_Conv::parse_expression(para.input.bands_to_print, para.sys.out_band_kb); }; sync_string(input.bands_to_print); this->add_item(item); } + { + Input_Item item("if_separate_k"); + item.annotation = "specify whether to write the partial charge densities for all k-points to individual files or merge them"; + read_sync_bool(input.if_separate_k); + this->add_item(item); + } { Input_Item item("symmetry"); item.annotation = "the control of symmetry"; diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index d2ed8b1ddc..65bfdee02a 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -50,6 +50,7 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.nbands, 8); EXPECT_EQ(param.inp.nbands_sto, 256); EXPECT_EQ(param.inp.nbands_istate, 5); + EXPECT_FALSE(param.inp.if_separate_k); EXPECT_EQ(param.inp.pw_seed, 1); EXPECT_EQ(param.inp.emin_sto, 0.0); EXPECT_EQ(param.inp.emax_sto, 0.0); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index 7b47adecca..482c87fb15 100644 --- a/source/module_io/test/support/INPUT +++ b/source/module_io/test/support/INPUT @@ -19,6 +19,7 @@ kspacing 0 #unit in 1/bohr, should be > 0, default is 0 wh min_dist_coef 0.2 #factor related to the allowed minimum distance between two atoms nbands 8 #number of bands nbands_istate 5 #number of bands around Fermi level for get_pchg calulation +if_separate_k false #specify whether to write the partial charge densities for all k-points to individual files or merge them symmetry 1 #the control of symmetry init_vel False #read velocity from STRU or not symmetry_prec 1e-06 #accuracy for symmetry diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 764f3ca2f7..52ebe51d4d 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -39,6 +39,7 @@ struct Input_para int nbands_istate = 5; ///< number of bands around fermi level for get_pchg calculation. std::string bands_to_print = ""; ///< specify the bands to be calculated in the get_pchg ///< calculation, formalism similar to ocp_set. + bool if_separate_k = false; ///< whether to write partial charge for all k-points to individual files or merge them /* symmetry level: -1, no symmetry at all; 0, only basic time reversal would be considered; diff --git a/source/module_parameter/system_parameter.h b/source/module_parameter/system_parameter.h index 2c40b87e05..9624d29bb5 100644 --- a/source/module_parameter/system_parameter.h +++ b/source/module_parameter/system_parameter.h @@ -37,5 +37,7 @@ struct System_para double uramping = -10.0 / 13.6; /// U-Ramping method (Ry) std::vector hubbard_u = {}; ///< Hubbard Coulomb interaction parameter U (Ry) std::string global_calculation = "scf"; ///< global calculation type decided by "calculation" + + std::vector out_band_kb = {}; ///< return parsed bands_to_print as a vector of integers }; #endif \ No newline at end of file diff --git a/tests/integrate/312_NO_GO_wfc_istate/result.ref b/tests/integrate/312_NO_GO_wfc_istate/result.ref index e022a789de..66f747f6b9 100644 --- a/tests/integrate/312_NO_GO_wfc_istate/result.ref +++ b/tests/integrate/312_NO_GO_wfc_istate/result.ref @@ -1,5 +1,5 @@ -BAND1_SPIN0_CHG.cube 2.06078 -BAND2_SPIN0_CHG.cube 2.07625 -BAND3_SPIN0_CHG.cube 2.07027 -BAND4_SPIN0_CHG.cube 2.06983 +BAND1_GAMMA_SPIN1_CHG.cube 2.06078 +BAND2_GAMMA_SPIN1_CHG.cube 2.07625 +BAND3_GAMMA_SPIN1_CHG.cube 2.07027 +BAND4_GAMMA_SPIN1_CHG.cube 2.06983 totaltimeref 0.16856