diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 4a8b1a096d..781886e65b 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -422,8 +422,9 @@ - [lr\_nstates](#lr_nstates) - [abs\_wavelen\_range](#abs_wavelen_range) - [out\_wfc\_lr](#out_wfc_lr) -[back to top](#full-list-of-input-keywords) + - [abs\_broadening](#abs_broadening) +[back to top](#full-list-of-input-keywords) ## System variables These variables are used to control general system parameters. @@ -1511,8 +1512,10 @@ These variables are used to control the output of properties. ### out_chg -- **Type**: Integer +- **Type**: Integer \[Integer\](optional) - **Description**: + + The first integer controls whether to output the charge density on real space grids: - 1. Output the charge density (in Bohr^-3) on real space grids into the density files in the folder `OUT.${suffix}`. The files are named as: - nspin = 1: SPIN1_CHG.cube; - nspin = 2: SPIN1_CHG.cube, and SPIN2_CHG.cube; @@ -1521,13 +1524,16 @@ These variables are used to control the output of properties. - nspin = 1: SPIN1_CHG_INI.cube - nspin = 2: SPIN1_CHG_INI.cube, and SPIN2_CHG_INI.cube; - nspin = 4: SPIN1_CHG_INI.cube, SPIN2_CHG_INI.cube, SPIN3_CHG_INI.cube, and SPIN4_CHG_INI.cube. + + The second integer controls the precision of the charge density output, if not given, will use `3` as default. For purpose restarting from this file and other high-precision involved calculation, recommend to use `10`. + --- The circle order of the charge density on real space grids is: x is the outer loop, then y and finally z (z is moving fastest). If EXX(exact exchange) is calculated, (i.e. *[dft_fuctional](#dft_functional)==hse/hf/pbe0/scan0/opt_orb* or *[rpa](#rpa)==True*), the Hexx(R) files will be output in the folder `OUT.${suffix}` too, which can be read in NSCF calculation. In molecular dynamics calculations, the output frequency is controlled by [out_interval](#out_interval). -- **Default**: 0 +- **Default**: 0 3 ### out_pot @@ -1612,7 +1618,7 @@ These variables are used to control the output of properties. ### out_band -- **Type**: Boolean Integer(optional) +- **Type**: Boolean \[Integer\](optional) - **Description**: Whether to output the band structure (in eV), optionally output precision can be set by a second parameter, default is 8. For more information, refer to the [band.md](../elec_properties/band.md) - **Default**: False @@ -1656,7 +1662,7 @@ These variables are used to control the output of properties. ### out_mat_hs -- **Type**: Boolean Integer(optional) +- **Type**: Boolean \[Integer\](optional) - **Availability**: Numerical atomic orbital basis - **Description**: Whether to print the upper triangular part of the Hamiltonian matrices (in Ry) and overlap matrices for each k point into files in the directory `OUT.${suffix}`. The second number controls precision. For more information, please refer to [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs). Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag). - **Default**: False 8 diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 66b7306fcc..93b2c1a897 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -122,7 +122,7 @@ void ESolver_FP::after_scf(const int istep) if (istep % PARAM.inp.out_interval == 0) { // 3) write charge density - if (PARAM.inp.out_chg) + if (PARAM.inp.out_chg[0]) { for (int is = 0; is < GlobalV::NSPIN; is++) { @@ -153,7 +153,7 @@ void ESolver_FP::after_scf(const int istep) this->pw_rhod->nz, this->pelec->eferm.get_efval(is), &(GlobalC::ucell), - 3, + PARAM.inp.out_chg[1], 1); if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) { diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 8bc582659c..aa7d856b35 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -230,7 +230,7 @@ void ESolver_KS_PW::before_scf(const int istep) this->pelec->init_scf(istep, this->sf.strucFac); //! output the initial charge density - if (PARAM.inp.out_chg == 2) + if (PARAM.inp.out_chg[0] == 2) { for (int is = 0; is < GlobalV::NSPIN; is++) { @@ -472,7 +472,7 @@ void ESolver_KS_PW::iter_finish(int& iter) if (this->out_freq_elec && iter % this->out_freq_elec == 0) { - if (PARAM.inp.out_chg > 0) + if (PARAM.inp.out_chg[0] > 0) { for (int is = 0; is < GlobalV::NSPIN; is++) { diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index c36e2732ba..ab34e79c02 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -220,7 +220,7 @@ void ESolver_KS_LCAO::before_scf(const int istep) this->pelec->init_scf(istep, this->sf.strucFac); //! output the initial charge density - if (PARAM.inp.out_chg == 2) + if (PARAM.inp.out_chg[0] == 2) { for (int is = 0; is < GlobalV::NSPIN; is++) { diff --git a/source/module_hamilt_lcao/module_dftu/dftu_io.cpp b/source/module_hamilt_lcao/module_dftu/dftu_io.cpp index 3fba8fc924..2d9bccd31a 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu_io.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu_io.cpp @@ -52,7 +52,7 @@ void DFTU::output() //Write onsite.dm std::ofstream ofdftu; - if(PARAM.inp.out_chg){ + if(PARAM.inp.out_chg[0]){ if(GlobalV::MY_RANK == 0){ ofdftu.open(GlobalV::global_out_dir + "onsite.dm"); } diff --git a/source/module_io/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp index 19ffd9effe..1f38e9fe69 100644 --- a/source/module_io/read_input_item_output.cpp +++ b/source/module_io/read_input_item_output.cpp @@ -37,14 +37,25 @@ void ReadInput::item_output() } { Input_Item item("out_chg"); - item.annotation = ">0 output charge density for selected electron steps"; - item.reset_value = [](const Input_Item& item, Parameter& para) { - if (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg") + item.annotation = ">0 output charge density for selected electron steps" + ", second parameter controls the precision, default is 3."; + item.read_value = [](const Input_Item& item, Parameter& para) { + size_t count = item.get_size(); + std::vector out_chg(count, -1); + std::transform(item.str_values.begin(), item.str_values.end(), out_chg.begin(), [](std::string s) { return std::stoi(s); }); + // assign non-negative values to para.input.out_chg + for (size_t i = 0; i < count; ++i) { - para.input.out_chg = 1; + if (out_chg[i] >= 0) + { + para.input.out_chg[i] = out_chg[i]; + } } }; - read_sync_int(input.out_chg); + item.reset_value = [](const Input_Item& item, Parameter& para) { + para.input.out_chg[0] = (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg") ? 1 : para.input.out_chg[0]; + }; + sync_intvec(input.out_chg, 2, 0); this->add_item(item); } { diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index a673feb2df..e40959e4ab 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -179,7 +179,8 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.chg_extrap, "atomic"); EXPECT_EQ(param.inp.out_freq_elec, 0); EXPECT_EQ(param.inp.out_freq_ion, 0); - EXPECT_EQ(param.inp.out_chg, 0); + EXPECT_EQ(param.inp.out_chg[0], 0); + EXPECT_EQ(param.inp.out_chg[1], 3); EXPECT_EQ(param.inp.out_dm, 0); EXPECT_EQ(param.inp.out_dm1, 0); EXPECT_EQ(param.inp.deepks_out_labels, 0); diff --git a/source/module_io/test_serial/read_input_item_test.cpp b/source/module_io/test_serial/read_input_item_test.cpp index b36915b693..826582824a 100644 --- a/source/module_io/test_serial/read_input_item_test.cpp +++ b/source/module_io/test_serial/read_input_item_test.cpp @@ -347,9 +347,9 @@ TEST_F(InputTest, Item_test) { // out_chg auto it = find_label("out_chg", readinput.input_lists); param.input.calculation = "get_wf"; - param.input.out_chg = 0; + param.input.out_chg = {0}; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.out_chg, 1); + EXPECT_EQ(param.input.out_chg[0], 1); } { // out_pot auto it = find_label("out_pot", readinput.input_lists); diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index 05bba4c2a9..d1d417cb96 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -310,7 +310,7 @@ struct Input_para ///< 0: output only when converged int out_freq_ion = 0; ///< the frequency ( >= 0 ) of ionic step to output charge density; ///< 0: output only when ion steps are finished - int out_chg = 0; ///< output charge density. 0: no; 1: yes + std::vector out_chg = {0, 3}; ///< output charge density. 0: no; 1: yes int out_pot = 0; ///< yes or no int out_wfc_pw = 0; ///< 0: no; 1: txt; 2: dat bool out_wfc_r = false; ///< 0: no; 1: yes diff --git a/source/module_ri/exx_lip.cpp b/source/module_ri/exx_lip.cpp index 84e7f2fa30..1fc7d6d9c4 100644 --- a/source/module_ri/exx_lip.cpp +++ b/source/module_ri/exx_lip.cpp @@ -29,16 +29,21 @@ void Exx_Lip::cal_exx() //t_phi_cal += my_time(t); judge_singularity(ik); - for (int iw_l = 0; iw_l < GlobalV::NLOCAL; ++iw_l) - for (int iw_r = 0; iw_r < GlobalV::NLOCAL; ++iw_r) + for (int iw_l = 0; iw_l < GlobalV::NLOCAL; ++iw_l) { + for (int iw_r = 0; iw_r < GlobalV::NLOCAL; ++iw_r) { sum1[iw_l * GlobalV::NLOCAL + iw_r] = T(0.0, 0.0); +} +} if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { sum2_factor = 0.0; - if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) - for (int iw_l = 0; iw_l < GlobalV::NLOCAL; ++iw_l) - for (int iw_r = 0; iw_r < GlobalV::NLOCAL; ++iw_r) + if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) { + for (int iw_l = 0; iw_l < GlobalV::NLOCAL; ++iw_l) { + for (int iw_r = 0; iw_r < GlobalV::NLOCAL; ++iw_r) { sum3[iw_l][iw_r] = T(0.0, 0.0); +} +} +} } for (int iq_tmp = iq_vecik; iq_tmp < iq_vecik + q_pack->kv_ptr->get_nks() / GlobalV::NSPIN; ++iq_tmp) // !!! k_point @@ -52,9 +57,11 @@ void Exx_Lip::cal_exx() { b_cal(ik, iq, ib); //t_b_cal += my_time(t); - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) - if (iq == iq_vecik) + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { + if (iq == iq_vecik) { sum3_cal(iq, ib); +} +} //t_sum3_cal += my_time(t); b_sum(iq, ib); //t_b_sum += my_time(t); @@ -74,8 +81,9 @@ void Exx_Lip::cal_exx() ofs("Hexxk_" + ModuleBase::GlobalFunc::TO_STRING(istep++) + "_" + ModuleBase::GlobalFunc::TO_STRING(ik) + "_" + ModuleBase::GlobalFunc::TO_STRING(GlobalV::MY_RANK)); for (int i = 0; i != GlobalV::NLOCAL; ++i) { - for (int j = 0; j != GlobalV::NLOCAL; ++j) + for (int j = 0; j != GlobalV::NLOCAL; ++j) { ofs << exx_matrix[ik][i][j] << "\t"; +} ofs << std::endl; } }; @@ -191,13 +199,14 @@ Exx_Lip::Exx_Lip( sum1.resize(GlobalV::NLOCAL * GlobalV::NLOCAL); - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) { b0.resize(GlobalV::NLOCAL); sum3.resize(GlobalV::NLOCAL); for (int iw_l = 0; iw_l < GlobalV::NLOCAL; ++iw_l) { sum3[iw_l].resize(GlobalV::NLOCAL); } } +} exx_matrix.resize(k_pack->kv_ptr->get_nks()); for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) @@ -215,19 +224,20 @@ Exx_Lip::Exx_Lip( template Exx_Lip::~Exx_Lip() { - if (k_pack)delete k_pack->hvec_array; + if (k_pack) {delete k_pack->hvec_array; +} delete k_pack; if (GlobalV::init_chg == "atomic") { - q_pack = NULL; + q_pack = nullptr; } else if (GlobalV::init_chg == "file") { - delete q_pack->kv_ptr; q_pack->kv_ptr = NULL; - delete q_pack->wf_ptr; q_pack->wf_ptr = NULL; + delete q_pack->kv_ptr; q_pack->kv_ptr = nullptr; + delete q_pack->wf_ptr; q_pack->wf_ptr = nullptr; // delete[] q_pack->hvec_array; q_pack->hvec_array=NULL; - delete q_pack; q_pack = NULL; + delete q_pack; q_pack = nullptr; } } @@ -235,15 +245,20 @@ template void Exx_Lip::wf_wg_cal() { ModuleBase::TITLE("Exx_Lip", "wf_wg_cal"); - if (GlobalV::NSPIN == 1) - for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) - for (int ib = 0; ib < GlobalV::NBANDS; ++ib) + if (GlobalV::NSPIN == 1) { + for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) { + for (int ib = 0; ib < GlobalV::NBANDS; ++ib) { k_pack->wf_wg(ik, ib) = k_pack->pelec->wg(ik, ib) / 2; - else if (GlobalV::NSPIN == 2) - for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) - for (int ib = 0; ib < GlobalV::NBANDS; ++ib) +} +} + } else if (GlobalV::NSPIN == 2) { + for (int ik = 0; ik < k_pack->kv_ptr->get_nks(); ++ik) { + for (int ib = 0; ib < GlobalV::NBANDS; ++ib) { k_pack->wf_wg(ik, ib) = k_pack->pelec->wg(ik, ib); } +} +} +} template void Exx_Lip::phi_cal(k_package* kq_pack, int ikq) @@ -358,19 +373,21 @@ void Exx_Lip::qkg2_exp(int ik, int iq) const Real qkg2 = ((q_pack->kv_ptr->kvec_c[iq] - k_pack->kv_ptr->kvec_c[ik] + rho_basis->gcar[ig]) * (ModuleBase::TWO_PI / ucell_ptr->lat0)).norm2(); if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { - if (std::abs(qkg2) < 1e-10) + if (std::abs(qkg2) < 1e-10) { recip_qkg2[ig] = 0.0; // 0 to ignore bb/qkg2 when qkg2==0 - else + } else { recip_qkg2[ig] = 1.0 / qkg2; +} sum2_factor += recip_qkg2[ig] * exp(-info.lambda * qkg2); recip_qkg2[ig] = sqrt(recip_qkg2[ig]); } else if (Conv_Coulomb_Pot_K::Ccp_Type::Hse == info.ccp_type) { - if (std::abs(qkg2) < 1e-10) + if (std::abs(qkg2) < 1e-10) { recip_qkg2[ig] = 1.0 / (2 * info.hse_omega); - else + } else { recip_qkg2[ig] = sqrt((1 - exp(-qkg2 / (4 * info.hse_omega * info.hse_omega))) / qkg2); +} } } } @@ -408,12 +425,15 @@ void Exx_Lip::b_cal(int ik, int iq, int ib) } T* const b_w = &b[iw * rho_basis->npw]; rho_basis->real2recip( porter, b_w); - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) - if ((iq == iq_vecik) && (gzero_rank_in_pool == GlobalV::RANK_IN_POOL)) /// need to check while use k_point parallel + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { + if ((iq == iq_vecik) && (gzero_rank_in_pool == GlobalV::RANK_IN_POOL)) { /// need to check while use k_point parallel b0[iw] = b_w[rho_basis->ig_gge0]; +} +} - for (size_t ig = 0; ig < rho_basis->npw; ++ig) + for (size_t ig = 0; ig < rho_basis->npw; ++ig) { b_w[ig] *= recip_qkg2[ig]; +} } delete [] porter; } @@ -421,11 +441,14 @@ void Exx_Lip::b_cal(int ik, int iq, int ib) template void Exx_Lip::sum3_cal(int iq, int ib) { - if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) - for (int iw_l = 0; iw_l < GlobalV::NLOCAL; ++iw_l) - for (int iw_r = 0; iw_r < GlobalV::NLOCAL; ++iw_r) + if (gzero_rank_in_pool == GlobalV::RANK_IN_POOL) { + for (int iw_l = 0; iw_l < GlobalV::NLOCAL; ++iw_l) { + for (int iw_r = 0; iw_r < GlobalV::NLOCAL; ++iw_r) { sum3[iw_l][iw_r] += b0[iw_l] * conj(b0[iw_r]) * (Real)q_pack->wf_wg(iq, ib); } +} +} +} template void Exx_Lip::b_sum(int iq, int ib) // Peize Lin change 2019-04-14 @@ -449,24 +472,28 @@ void Exx_Lip::sum_all(int ik) Real fourpi_div_omega = 4 * (Real)(ModuleBase::PI / ucell_ptr->omega); Real spin_fac = 2.0; #ifdef __MPI - if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) + if (Conv_Coulomb_Pot_K::Ccp_Type::Ccp == info.ccp_type || Conv_Coulomb_Pot_K::Ccp_Type::Hf == info.ccp_type) { MPI_Reduce(&sum2_factor, &sum2_factor_g, 1, MPI_DOUBLE, MPI_SUM, gzero_rank_in_pool, POOL_WORLD); +} #endif - for (size_t iw_l = 1; iw_l < GlobalV::NLOCAL; ++iw_l) - for (size_t iw_r = 0; iw_r < iw_l; ++iw_r) + for (size_t iw_l = 1; iw_l < GlobalV::NLOCAL; ++iw_l) { + for (size_t iw_r = 0; iw_r < iw_l; ++iw_r) { sum1[iw_l * GlobalV::NLOCAL + iw_r] = conj(sum1[iw_r * GlobalV::NLOCAL + iw_l]); // Peize Lin add conj 2019-04-14 +} +} for (int iw_l = 0; iw_l < GlobalV::NLOCAL; ++iw_l) { for( int iw_r=0; iw_rkv_ptr->get_nks() / GlobalV::NSPIN) * sum3[iw_l][iw_r]); } +} } } } @@ -549,8 +576,9 @@ void Exx_Lip::exx_energy_cal() template void Exx_Lip::write_q_pack() const { - if (PARAM.inp.out_chg == 0) + if (PARAM.inp.out_chg[0] == 0) { return; +} if (!GlobalV::RANK_IN_POOL) {