diff --git a/source/module_lr/hamilt_casida.cpp b/source/module_lr/hamilt_casida.cpp index 126479ad00..c9135b62e0 100644 --- a/source/module_lr/hamilt_casida.cpp +++ b/source/module_lr/hamilt_casida.cpp @@ -7,8 +7,8 @@ namespace LR ModuleBase::TITLE("HamiltCasidaLR", "matrix"); int npairs = this->nocc * this->nvirt; std::vector Amat_full(this->nks * npairs * this->nks * npairs, 0.0); - for (int isk = 0;isk < this->nks;++isk) - for (int j = 0;j < nocc;++j) + for (int isk = 0;isk < this->nks;++isk) { + for (int j = 0;j < nocc;++j) { for (int b = 0;b < nvirt;++b) {//calculate A^{ai} for each bj int bj = j * nvirt + b; @@ -18,7 +18,8 @@ namespace LR // X_bj(0, 0, lj * this->pX->get_row_size() + lb) = this->one(); int lj = this->pX->global2local_col(j); int lb = this->pX->global2local_row(b); - if (this->pX->in_this_processor(b, j)) X_bj(0, 0, isk * this->pX->get_local_size() + lj * this->pX->get_row_size() + lb) = this->one(); + if (this->pX->in_this_processor(b, j)) { X_bj(0, 0, isk * this->pX->get_local_size() + lj * this->pX->get_row_size() + lb) = this->one(); +} psi::Psi A_aibj(1, 1, this->nks * this->pX->get_local_size()); // k1-first A_aibj.zero_out(); @@ -30,11 +31,14 @@ namespace LR } // reduce ai for a fixed bj A_aibj.fix_kb(0, 0); - for (int isk_ai = 0;isk_ai < this->nks;++isk_ai) + for (int isk_ai = 0;isk_ai < this->nks;++isk_ai) { LR_Util::gather_2d_to_full(*this->pX, &A_aibj.get_pointer()[isk_ai * this->pX->get_local_size()], Amat_full.data() + kbj * this->nks * npairs /*col, bj*/ + isk_ai * npairs/*row, ai*/, false, this->nvirt, this->nocc); +} } +} +} // output Amat std::cout << "Amat_full:" << std::endl; for (int i = 0;i < this->nks * npairs;++i) diff --git a/source/module_lr/hamilt_casida.h b/source/module_lr/hamilt_casida.h index ca3a0d9647..394b54db0a 100644 --- a/source/module_lr/hamilt_casida.h +++ b/source/module_lr/hamilt_casida.h @@ -58,7 +58,8 @@ namespace LR { delete this->ops; } - for (auto& d : this->DM_trans)delete d; + for (auto& d : this->DM_trans) {delete d; +} }; hamilt::HContainer* getHR() { return this->hR; } diff --git a/source/module_lr/operator_casida/operator_lr_diag.h b/source/module_lr/operator_casida/operator_lr_diag.h index f0562670a1..703120c11c 100644 --- a/source/module_lr/operator_casida/operator_lr_diag.h +++ b/source/module_lr/operator_casida/operator_lr_diag.h @@ -22,14 +22,16 @@ namespace LR this->act_type = 2; this->cal_type = hamilt::calculation_type::no; this->eig_ks_diff.create(nks, pX->get_local_size(), false); - for (int ik = 0;ik < nks;++ik) - for (int io = 0;io < pX->get_col_size();++io) //nocc_local + for (int ik = 0;ik < nks;++ik) { + for (int io = 0;io < pX->get_col_size();++io) { //nocc_local for (int iv = 0;iv < pX->get_row_size();++iv) //nvirt_local { int io_g = pX->local2global_col(io); int iv_g = pX->local2global_row(iv); this->eig_ks_diff(ik, io * pX->get_row_size() + iv) = eig_ks(ik, nocc + iv_g) - eig_ks(ik, io_g); } +} +} }; void init(const int ik_in) override {}; diff --git a/source/module_lr/potentials/pot_hxc_lrtd.cpp b/source/module_lr/potentials/pot_hxc_lrtd.cpp index 26c2329879..b4073d96ea 100644 --- a/source/module_lr/potentials/pot_hxc_lrtd.cpp +++ b/source/module_lr/potentials/pot_hxc_lrtd.cpp @@ -39,12 +39,13 @@ namespace LR #ifdef USE_LIBXC const int nspin = v_eff.nr; v_eff += elecstate::H_Hartree_pw::v_hartree(*ucell, const_cast(this->rho_basis_), v_eff.nr, rho); - if (xc_kernel == "rpa" || xc_kernel == "hf") return; - else if (XC_Functional::get_func_type() == 1)//LDA - if (1 == nspin)// for LDA-spin0, just fxc*rho where fxc=v2rho2; for GGA, v2rho2 has been replaced by the true fxc - for (int ir = 0;ir < nrxx;++ir) + if (xc_kernel == "rpa" || xc_kernel == "hf") { return; + } else if (XC_Functional::get_func_type() == 1) {//LDA + if (1 == nspin) {// for LDA-spin0, just fxc*rho where fxc=v2rho2; for GGA, v2rho2 has been replaced by the true fxc + for (int ir = 0;ir < nrxx;++ir) { v_eff(0, ir) += ModuleBase::e2 * fxc.get_kernel("v2rho2").at(ir) * rho[0][ir]; - else if (2 == nspin) +} + } else if (2 == nspin) { for (int ir = 0;ir < nrxx;++ir) { const int irs0 = 2 * ir; @@ -55,10 +56,11 @@ namespace LR v_eff(1, ir) += ModuleBase::e2 * fxc.get_kernel("v2rho2").at(irs1) * rho[0][ir] + fxc.get_kernel("v2rho2").at(irs2) * rho[1][ir]; } - else //remain for spin 4 + } else { //remain for spin 4 throw std::domain_error("nspin =" + std::to_string(nspin) + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); - else if (XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4) // GGA or HYB_GGA +} + } else if (XC_Functional::get_func_type() == 2 || XC_Functional::get_func_type() == 4) // GGA or HYB_GGA { if (1 == nspin) { @@ -67,34 +69,44 @@ namespace LR // test: output drho double thr = 1e-1; auto out_thr = [this, &thr](const double* v) { - for (int ir = 0;ir < nrxx;++ir) if (std::abs(v[ir]) > thr) std::cout << v[ir] << " "; + for (int ir = 0;ir < nrxx;++ir) { if (std::abs(v[ir]) > thr) { std::cout << v[ir] << " "; +} +} std::cout << std::endl;}; auto out_thr3 = [this, &thr](const std::vector>& v) { - for (int ir = 0;ir < nrxx;++ir) if (std::abs(v.at(ir).x) > thr) std::cout << v.at(ir).x << " "; + for (int ir = 0;ir < nrxx;++ir) { if (std::abs(v.at(ir).x) > thr) { std::cout << v.at(ir).x << " "; +} +} std::cout << std::endl; - for (int ir = 0;ir < nrxx;++ir) if (std::abs(v.at(ir).y) > thr) std::cout << v.at(ir).y << " "; + for (int ir = 0;ir < nrxx;++ir) { if (std::abs(v.at(ir).y) > thr) { std::cout << v.at(ir).y << " "; +} +} std::cout << std::endl; - for (int ir = 0;ir < nrxx;++ir) if (std::abs(v.at(ir).z) > thr) std::cout << v.at(ir).z << " "; + for (int ir = 0;ir < nrxx;++ir) { if (std::abs(v.at(ir).z) > thr) { std::cout << v.at(ir).z << " "; +} +} std::cout << std::endl;}; std::vector vxc_tmp(nrxx, 0.0); //1. $\partial E/\partial\rho = 2f^{\rho\sigma}*\nabla\rho*\rho_1+4f^{\sigma\sigma}\nabla\rho(\nabla\rho\cdot\nabla\rho_1)+2v^\sigma\nabla\rho_1$ std::vector> e_drho(nrxx); - for (int ir = 0;ir < nrxx;++ir) + for (int ir = 0;ir < nrxx;++ir) { e_drho[ir] = -(fxc.get_grad_kernel("2_v2rhosigma_drho").at(ir) * rho[0][ir] + fxc.get_grad_kernel("4_v2sigma2_drho").at(ir) * (fxc.get_grad_kernel("drho_gs").at(ir) * drho.at(ir)) + drho.at(ir) * fxc.get_kernel("vsigma").at(ir) * 2.); +} XC_Functional::grad_dot(e_drho.data(), vxc_tmp.data(), this->rho_basis_, ucell->tpiba); // 2. $f^{\rho\rho}\rho_1+2f^{\rho\sigma}\nabla\rho\cdot\nabla\rho_1$ - for (int ir = 0;ir < nrxx;++ir) + for (int ir = 0;ir < nrxx;++ir) { vxc_tmp[ir] += (fxc.get_kernel("v2rho2").at(ir) * rho[0][ir] + fxc.get_grad_kernel("2_v2rhosigma_drho").at(ir) * drho.at(ir)); +} BlasConnector::axpy(nrxx, ModuleBase::e2, vxc_tmp.data(), 1, v_eff.c, 1); } - else if (2 == nspin) // wrong, to be fixed + else if (2 == nspin) { // wrong, to be fixed for (int ir = 0;ir < nrxx;++ir) { const int irs0 = 2 * ir; @@ -105,14 +117,16 @@ namespace LR v_eff(1, ir) += ModuleBase::e2 * fxc.get_kernel("v2rho2").at(irs1) * rho[0][ir] + fxc.get_kernel("v2rho2").at(irs2) * rho[1][ir]; } - else //remain for spin 4 + } else { //remain for spin 4 throw std::domain_error("nspin =" + std::to_string(nspin) + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); +} } - else + else { #endif throw std::domain_error("GlobalV::XC_Functional::get_func_type() =" + std::to_string(XC_Functional::get_func_type()) + " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); +} ModuleBase::timer::tick("PotHxcLR", "cal_v_eff"); } diff --git a/source/module_lr/utils/ri_benchmark_io.cpp b/source/module_lr/utils/ri_benchmark_io.cpp index caeb2673bb..632e6ae089 100644 --- a/source/module_lr/utils/ri_benchmark_io.cpp +++ b/source/module_lr/utils/ri_benchmark_io.cpp @@ -6,7 +6,7 @@ namespace RI_Benchmark template TLRI read_Cs_ao(const std::string& file_path) { - const int natom, ncell, ia1, ia2, ic_1, ic_2, ic_3, nw1, nw2, nabf; + const int natom = 0, ncell = 0, ia1 = 0, ia2 = 0, ic_1 = 0, ic_2 = 0, ic_3 = 0, nw1 = 0, nw2 = 0, nabf = 0; std::ifstream infile; infile.open(file_path); infile >> natom >> ncell; // no use of ncell @@ -15,12 +15,15 @@ namespace RI_Benchmark while (infile.peek() != EOF) { infile >> ia1 >> ia2 >> ic_1 >> ic_2 >> ic_3 >> nw1 >> nw2 >> nabf; - ModuleBase::Vector3_Order box(ic_1, ic_2, ic_3); + Abfs::Vector3_Order box(ic_1, ic_2, ic_3); RI::Tensor tensor_cs({ nabf, na1, na2 }); - for (int i = 0; i != nw1; i++) - for (int j = 0; j != nw2; j++) - for (int mu = 0; mu != nabf; mu++) + for (int i = 0; i != nw1; i++) { + for (int j = 0; j != nw2; j++) { + for (int mu = 0; mu != nabf; mu++) { infile >> tensor_cs(mu, i, j); +} +} +} // no screening for data-structure consistency // if (loc_atp_index.count(ia1) && (*cs_ptr).absmax() >= threshold) Cs[ia1][{ia2, box}] = tensor_cs; @@ -32,7 +35,7 @@ namespace RI_Benchmark template TLRI read_Vs_abf(const std::string& file_path) { - const int natom, ncell, ia1, ia2, ic_1, ic_2, ic_3, nabf1, nabf2; + const int natom = 0, ncell = 0, ia1 = 0, ia2 = 0, ic_1, ic_2, ic_3, nabf1, nabf2; std::ifstream infile; infile.open(file_path); infile >> natom >> ncell; // no use of ncell @@ -43,9 +46,11 @@ namespace RI_Benchmark infile >> ia1 >> ia2 >> ic_1 >> ic_2 >> ic_3 >> nabf1 >> nabf2; ModuleBase::Vector3_Order box(ic_1, ic_2, ic_3); RI::Tensor tensor_cs({ n_mu,n_nu }); - for (int i = 0; i != nabf1; i++) - for (int j = 0; j != nabf2; j++) + for (int i = 0; i != nabf1; i++) { + for (int j = 0; j != nabf2; j++) { infile >> tensor_cs(i, j); +} +} // no screening for data-structure consistency // if (loc_atp_index.count(ia1) && (*cs_ptr).absmax() >= threshold) Vs[ia1][{ia2, box}] = tensor_cs; diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index e11a354e38..aba0ea5879 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -46,7 +46,7 @@ struct Input_para std::string pseudo_dir = ""; ///< directory of pseudopotential std::string orbital_dir = ""; ///< directory of orbital file double pseudo_rcut = 15.0; ///< cut-off radius for calculating msh - bool pseudo_mesh = 0; ///< 0: use msh to normalize radial wave functions; 1: + bool pseudo_mesh = false; ///< 0: use msh to normalize radial wave functions; 1: ///< use mesh, which is used in QE. int lmaxmax = 2; ///< maximum of l channels used std::string dft_functional = "default"; ///< input DFT functional. @@ -91,7 +91,7 @@ struct Input_para std::string wannier_card = "none"; ///< input card for wannier functions. double soc_lambda = 1.0; ///< The fraction of averaged SOC pseudopotential ///< is given by (1-soc_lambda) - bool cal_force = 0; ///< calculate the force + bool cal_force = false; ///< calculate the force int out_freq_ion = 0; ///< the frequency ( >= 0 ) of ionic step to output charge density ///< and wavefunction. 0: output only when ion steps are finished int elpa_num_thread = -1; ///< Number of threads need to use in elpa @@ -156,7 +156,7 @@ struct Input_para int cond_dtbatch = 0; ///< exp(iH*dt*cond_dtbatch) is expanded with Chebyshev expansion. int cond_smear = 1; ///< smearing method for conductivities 1: Gaussian 2: Lorentzian double cond_fwhm = 0.4; ///< FWHM for conductivities - bool cond_nonlocal = 1; ///< if calculate nonlocal effects + bool cond_nonlocal = true; ///< if calculate nonlocal effects // ============== #Parameters (4.Relaxation) =========================== std::string ks_solver = "default"; ///< xiaohui add 2013-09-01