diff --git a/source/source_basis/module_pw/test_gpu/pw_basis_C2C.cpp b/source/source_basis/module_pw/test_gpu/pw_basis_C2C.cpp index 95d076040e..9abaefab44 100644 --- a/source/source_basis/module_pw/test_gpu/pw_basis_C2C.cpp +++ b/source/source_basis/module_pw/test_gpu/pw_basis_C2C.cpp @@ -44,9 +44,6 @@ class PW_BASIS_K_GPU_TEST : public ::testing::Test int distribution_type = 1; bool xprime = false; const int nks = 1; - ModuleBase::Vector3* kvec_d; - kvec_d = new ModuleBase::Vector3[nks]; - kvec_d[0].set(0, 0, 0); // init const int mypool = 0; const int key = 1; diff --git a/source/source_cell/module_symmetry/symm_rho.cpp b/source/source_cell/module_symmetry/symm_rho.cpp index 754279ffbc..28ae30d5b4 100644 --- a/source/source_cell/module_symmetry/symm_rho.cpp +++ b/source/source_cell/module_symmetry/symm_rho.cpp @@ -9,6 +9,10 @@ void Symmetry::rho_symmetry( double *rho, { ModuleBase::timer::tick("Symmetry","rho_symmetry"); + assert(nr1>0); + assert(nr2>0); + assert(nr3>0); + // allocate flag for each FFT grid. bool* symflag = new bool[nr1 * nr2 * nr3]; for (int i=0; i, * make it orthogonal to phi ( = 0), and = nelec */ -void ESolver_OF::adjust_direction() +void ESolver_OF::adjust_direction(void) { // filter the high frequency term in direction if of_full_pw = false if (!PARAM.inp.of_full_pw) diff --git a/source/source_estate/module_charge/charge_mixing_preconditioner.cpp b/source/source_estate/module_charge/charge_mixing_preconditioner.cpp index 804705be99..00b5d42038 100644 --- a/source/source_estate/module_charge/charge_mixing_preconditioner.cpp +++ b/source/source_estate/module_charge/charge_mixing_preconditioner.cpp @@ -6,21 +6,32 @@ void Charge_Mixing::Kerker_screen_recip(std::complex* drhog) { - if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) { - return; -} + ModuleBase::TITLE("Charge_Mixing", "Kerker_screen_recip"); + + if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) + { + return; + } + + ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_recip"); + + const int nspin = PARAM.inp.nspin; + double fac = 0.0; double gg0 = 0.0; double amin = 0.0; /// consider a resize for mixing_angle - int resize_tmp = 1; - if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { resize_tmp = 2; -} + int resize_tmp = 1; + if (nspin == 4 && this->mixing_angle > 0) + { + resize_tmp = 2; + } /// implement Kerker for density and magnetization separately - for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is) + for (int is = 0; is < nspin / resize_tmp; ++is) { + const int is_idx = is * this->rhopw->npw; /// new mixing method only support nspin=2 not nspin=4 if (is >= 1) { @@ -29,10 +40,10 @@ void Charge_Mixing::Kerker_screen_recip(std::complex* drhog) #ifdef __DEBUG assert(is == 1); // make sure break works #endif - double is_mag = PARAM.inp.nspin - 1; + double is_mag = nspin - 1; //for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) //{ - // drhog[is * this->rhopw->npw + ig] *= 1; + // drhog[is_idx + ig] *= 1; //} break; } @@ -46,32 +57,49 @@ void Charge_Mixing::Kerker_screen_recip(std::complex* drhog) } gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2); + + const double gg0_amin = this->mixing_gg0_min / amin; + #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif for (int ig = 0; ig < this->rhopw->npw; ++ig) { double gg = this->rhopw->gg[ig]; - double filter_g = std::max(gg / (gg + gg0), this->mixing_gg0_min / amin); - drhog[is * this->rhopw->npw + ig] *= filter_g; + double filter_g = std::max(gg / (gg + gg0), gg0_amin); + drhog[is_idx + ig] *= filter_g; } } + + ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_recip"); return; } void Charge_Mixing::Kerker_screen_real(double* drhor) { - if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1) { - return; -} - /// consider a resize for mixing_angle + ModuleBase::TITLE("Charge_Mixing", "Kerker_screen_real"); + + if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1) + { + return; + } + + ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_real"); + + const int nspin = PARAM.inp.nspin; + assert(nspin==1 || nspin==2 || nspin==4); + + /// consider a resize for mixing_angle int resize_tmp = 1; - if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { resize_tmp = 2; -} + if (nspin == 4 && this->mixing_angle > 0) + { + resize_tmp = 2; + } - std::vector> drhog(this->rhopw->npw * PARAM.inp.nspin / resize_tmp); - std::vector drhor_filter(this->rhopw->nrxx * PARAM.inp.nspin / resize_tmp); - for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is) + std::vector> drhog(this->rhopw->npw * nspin / resize_tmp); + std::vector drhor_filter(this->rhopw->nrxx * nspin / resize_tmp); + + for (int is = 0; is < nspin / resize_tmp; ++is) { // Note after this process some G which is higher than Gmax will be filtered. // Thus we cannot use Kerker_screen_recip(drhog.data()) directly after it. @@ -82,7 +110,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) double gg0 = 0.0; double amin = 0.0; - for (int is = 0; is < PARAM.inp.nspin / resize_tmp; is++) + for (int is = 0; is < nspin / resize_tmp; is++) { if (is >= 1) @@ -92,8 +120,8 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) #ifdef __DEBUG assert(is == 1); /// make sure break works #endif - double is_mag = PARAM.inp.nspin - 1; - if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { is_mag = 1; + double is_mag = nspin - 1; + if (nspin == 4 && this->mixing_angle > 0) { is_mag = 1; } for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++) { @@ -111,6 +139,9 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) } gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2); + + const int is_idx = is * this->rhopw->npw; + const double gg0_amin = this->mixing_gg0_min / amin; #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif @@ -120,15 +151,15 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) // I have not decided how to handle gg=0 part, will be changed in future //if (gg == 0) //{ - // drhog[is * this->rhopw->npw + ig] *= 0; + // drhog[is_idx + ig] *= 0; // continue; //} - double filter_g = std::max(gg / (gg + gg0), this->mixing_gg0_min / amin); - drhog[is * this->rhopw->npw + ig] *= (1 - filter_g); + double filter_g = std::max(gg / (gg + gg0), gg0_amin); + drhog[is_idx + ig] *= (1 - filter_g); } } /// inverse FT - for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is) + for (int is = 0; is < nspin / resize_tmp; ++is) { this->rhopw->recip2real(drhog.data() + is * this->rhopw->npw, drhor_filter.data() + is * this->rhopw->nrxx); } @@ -136,8 +167,11 @@ void Charge_Mixing::Kerker_screen_real(double* drhor) #ifdef _OPENMP #pragma omp parallel for schedule(static, 512) #endif - for (int ir = 0; ir < this->rhopw->nrxx * PARAM.inp.nspin / resize_tmp; ir++) + for (int ir = 0; ir < this->rhopw->nrxx * nspin / resize_tmp; ir++) { drhor[ir] -= drhor_filter[ir]; } + + ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_real"); + return; } diff --git a/source/source_estate/module_charge/charge_mixing_residual.cpp b/source/source_estate/module_charge/charge_mixing_residual.cpp index 9e9c5e0131..bbc7c7b6be 100644 --- a/source/source_estate/module_charge/charge_mixing_residual.cpp +++ b/source/source_estate/module_charge/charge_mixing_residual.cpp @@ -9,11 +9,13 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec) { ModuleBase::TITLE("Charge_Mixing", "get_drho"); ModuleBase::timer::tick("Charge_Mixing", "get_drho"); + const int nspin = PARAM.inp.nspin; + assert(nspin==1 || nspin==2 || nspin==4); double drho = 0.0; if (PARAM.inp.scf_thr_type == 1) { - for (int is = 0; is < PARAM.inp.nspin; ++is) + for (int is = 0; is < nspin; ++is) { ModuleBase::GlobalFunc::NOTE("Perform FFT on rho(r) to obtain rho(G)."); chr->rhopw->real2recip(chr->rho[is], chr->rhog[is]); @@ -23,15 +25,15 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec) } ModuleBase::GlobalFunc::NOTE("Calculate the charge difference between rho(G) and rho_save(G)"); - std::vector> drhog(PARAM.inp.nspin * this->rhopw->npw); + std::vector> drhog(nspin * this->rhopw->npw); #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 512) #endif - for (int is = 0; is < PARAM.inp.nspin; ++is) + for (int is = 0; is < nspin; ++is) { for (int ig = 0; ig < this->rhopw->npw; ig++) { - drhog[is * rhopw->npw + ig] = chr->rhog[is][ig] - chr->rhog_save[is][ig]; + drhog[is * this->rhopw->npw + ig] = chr->rhog[is][ig] - chr->rhog_save[is][ig]; } } @@ -42,7 +44,7 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec) { // Note: Maybe it is wrong. // The inner_product_real function (L1-norm) is different from that (L2-norm) in mixing. - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < nspin; is++) { if (is != 0 && is != 3 && PARAM.globalv.domag_z) { diff --git a/source/source_estate/module_charge/charge_mixing_rho.cpp b/source/source_estate/module_charge/charge_mixing_rho.cpp index 6a965c5da2..38cd679f94 100644 --- a/source/source_estate/module_charge/charge_mixing_rho.cpp +++ b/source/source_estate/module_charge/charge_mixing_rho.cpp @@ -5,6 +5,12 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) { + ModuleBase::TITLE("Charge_Mixing", "mix_rho_recip"); + ModuleBase::timer::tick("Charge_Mixing", "mix_rho_recip"); + + const int nspin = PARAM.inp.nspin; + assert(nspin==1 || nspin==2 || nspin==4); + std::complex* rhog_in = nullptr; std::complex* rhog_out = nullptr; // for smooth part @@ -26,7 +32,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) = std::bind(&Charge_Mixing::inner_product_recip_hartree, this, std::placeholders::_1, std::placeholders::_2); // DIIS Mixing Only for smooth part, while high_frequency part is mixed by plain mixing method. - if (PARAM.inp.nspin == 1) + if (nspin == 1) { rhog_in = rhogs_in; rhog_out = rhogs_out; @@ -35,17 +41,17 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhog_out); } - else if (PARAM.inp.nspin == 2) + else if (nspin == 2) { // magnetic density std::complex *rhog_mag = nullptr; std::complex *rhog_mag_save = nullptr; const int npw = this->rhopw->npw; // allocate rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc] - rhog_mag = new std::complex[npw * PARAM.inp.nspin]; - rhog_mag_save = new std::complex[npw * PARAM.inp.nspin]; - ModuleBase::GlobalFunc::ZEROS(rhog_mag, npw * PARAM.inp.nspin); - ModuleBase::GlobalFunc::ZEROS(rhog_mag_save, npw * PARAM.inp.nspin); + rhog_mag = new std::complex[npw * nspin]; + rhog_mag_save = new std::complex[npw * nspin]; + ModuleBase::GlobalFunc::ZEROS(rhog_mag, npw * nspin); + ModuleBase::GlobalFunc::ZEROS(rhog_mag_save, npw * nspin); // get rhog_mag[is*ngmc] and rhog_mag_save[is*ngmc] for (int ig = 0; ig < npw; ig++) { @@ -84,7 +90,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhog_out); // get rhog[is][ngmc] from rhog_mag[is*ngmc] - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < nspin; is++) { ModuleBase::GlobalFunc::ZEROS(chr->rhog[is], npw); } @@ -106,7 +112,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) } } } - else if (PARAM.inp.nspin == 4 && PARAM.inp.mixing_angle <= 0) + else if (nspin == 4 && PARAM.inp.mixing_angle <= 0) { // normal broyden mixing for {rho, mx, my, mz} rhog_in = rhogs_in; @@ -135,7 +141,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhog_out); } - else if (PARAM.inp.nspin == 4 && PARAM.inp.mixing_angle > 0) + else if (nspin == 4 && PARAM.inp.mixing_angle > 0) { // special broyden mixing for {rho, |m|} proposed by J. Phys. Soc. Jpn. 82 (2013) 114706 // here only consider the case of mixing_angle = 1, which mean only change |m| and keep angle fixed @@ -154,9 +160,13 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) for (int ir = 0; ir < nrxx; ir++) { // |m| for rho - rho_magabs[ir] = std::sqrt(chr->rho[1][ir] * chr->rho[1][ir] + chr->rho[2][ir] * chr->rho[2][ir] + chr->rho[3][ir] * chr->rho[3][ir]); + rho_magabs[ir] = std::sqrt(chr->rho[1][ir] * chr->rho[1][ir] + + chr->rho[2][ir] * chr->rho[2][ir] + + chr->rho[3][ir] * chr->rho[3][ir]); // |m| for rho_save - rho_magabs_save[ir] = std::sqrt(chr->rho_save[1][ir] * chr->rho_save[1][ir] + chr->rho_save[2][ir] * chr->rho_save[2][ir] + chr->rho_save[3][ir] * chr->rho_save[3][ir]); + rho_magabs_save[ir] = std::sqrt(chr->rho_save[1][ir] * chr->rho_save[1][ir] + + chr->rho_save[2][ir] * chr->rho_save[2][ir] + + chr->rho_save[3][ir] * chr->rho_save[3][ir]); } // allocate memory for rhog_magabs and rhog_magabs_save const int npw = this->rhopw->npw; @@ -203,10 +213,14 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) // use new |m| and angle to update {mx, my, mz} for (int ig = 0; ig < npw; ig++) { - chr->rhog[0][ig] = rhog_magabs[ig]; // rhog - double norm = std::sqrt(chr->rho[1][ig] * chr->rho[1][ig] + chr->rho[2][ig] * chr->rho[2][ig] + chr->rho[3][ig] * chr->rho[3][ig]); - if (std::abs(norm) < 1e-10) { continue; -} + chr->rhog[0][ig] = rhog_magabs[ig]; // rhog + double norm = std::sqrt(chr->rho[1][ig] * chr->rho[1][ig] + + chr->rho[2][ig] * chr->rho[2][ig] + + chr->rho[3][ig] * chr->rho[3][ig]); + if (std::abs(norm) < 1e-10) + { + continue; + } double rescale_tmp = rho_magabs[npw + ig] / norm; chr->rho[1][ig] *= rescale_tmp; chr->rho[2][ig] *= rescale_tmp; @@ -222,7 +236,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) if ( PARAM.globalv.double_grid) { // plain mixing for high_frequencies - const int ndimhf = (this->rhodpw->npw - this->rhopw->npw) * PARAM.inp.nspin; + const int ndimhf = (this->rhodpw->npw - this->rhopw->npw) * nspin; this->mixing_highf->plain_mix(rhoghf_out, rhoghf_in, rhoghf_out, ndimhf, nullptr); // combine smooth part and high_frequency part @@ -231,7 +245,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) } // rhog to rho - if (PARAM.inp.nspin == 4 && PARAM.inp.mixing_angle > 0) + if (nspin == 4 && PARAM.inp.mixing_angle > 0) { // only tranfer rhog[0] // do not support double_grid, use rhopw directly @@ -239,7 +253,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) } else { - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < nspin; is++) { // use rhodpw for double_grid // rhodpw is the same as rhopw for ! PARAM.globalv.double_grid @@ -249,10 +263,10 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) // For kinetic energy density if ((XC_Functional::get_ked_flag()) && mixing_tau) { - std::vector> kin_g(PARAM.inp.nspin * rhodpw->npw); - std::vector> kin_g_save(PARAM.inp.nspin * rhodpw->npw); + std::vector> kin_g(nspin * rhodpw->npw); + std::vector> kin_g_save(nspin * rhodpw->npw); // FFT to get kin_g and kin_g_save - for (int is = 0; is < PARAM.inp.nspin; ++is) + for (int is = 0; is < nspin; ++is) { rhodpw->real2recip(chr->kin_r[is], &kin_g[is * rhodpw->npw]); rhodpw->real2recip(chr->kin_r_save[is], &kin_g_save[is * rhodpw->npw]); @@ -277,7 +291,7 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) if ( PARAM.globalv.double_grid) { // simple mixing for high_frequencies - const int ndimhf = (this->rhodpw->npw - this->rhopw->npw) * PARAM.inp.nspin; + const int ndimhf = (this->rhodpw->npw - this->rhopw->npw) * nspin; this->mixing_highf->plain_mix(taughf_out, taughf_in, taughf_out, ndimhf, nullptr); // combine smooth part and high_frequency part @@ -286,22 +300,28 @@ void Charge_Mixing::mix_rho_recip(Charge* chr) } // kin_g to kin_r - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < nspin; is++) { rhodpw->recip2real(&kin_g[is * rhodpw->npw], chr->kin_r[is]); } } - + ModuleBase::timer::tick("Charge_Mixing", "mix_rho_recip"); return; } void Charge_Mixing::mix_rho_real(Charge* chr) { + ModuleBase::TITLE("Charge_Mixing", "mix_rho_real"); + ModuleBase::timer::tick("Charge_Mixing", "mix_rho_real"); + + const int nspin = PARAM.inp.nspin; + assert(nspin==1 || nspin==2 || nspin==4); + double* rhor_in=nullptr; double* rhor_out=nullptr; - if (PARAM.inp.nspin == 1) + if (nspin == 1) { rhor_in = chr->rho_save[0]; rhor_out = chr->rho[0]; @@ -312,17 +332,17 @@ void Charge_Mixing::mix_rho_real(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhor_out); } - else if (PARAM.inp.nspin == 2) + else if (nspin == 2) { // magnetic density double *rho_mag = nullptr; double *rho_mag_save = nullptr; const int nrxx = this->rhopw->nrxx; // allocate rho_mag[is*nnrx] and rho_mag_save[is*nnrx] - rho_mag = new double[nrxx * PARAM.inp.nspin]; - rho_mag_save = new double[nrxx * PARAM.inp.nspin]; - ModuleBase::GlobalFunc::ZEROS(rho_mag, nrxx * PARAM.inp.nspin); - ModuleBase::GlobalFunc::ZEROS(rho_mag_save, nrxx * PARAM.inp.nspin); + rho_mag = new double[nrxx * nspin]; + rho_mag_save = new double[nrxx * nspin]; + ModuleBase::GlobalFunc::ZEROS(rho_mag, nrxx * nspin); + ModuleBase::GlobalFunc::ZEROS(rho_mag_save, nrxx * nspin); // get rho_mag[is*nnrx] and rho_mag_save[is*nnrx] for (int ir = 0; ir < nrxx; ir++) { @@ -362,7 +382,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhor_out); // get new rho[is][nrxx] from rho_mag[is*nrxx] - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < nspin; is++) { ModuleBase::GlobalFunc::ZEROS(chr->rho[is], nrxx); //ModuleBase::GlobalFunc::ZEROS(rho_save[is], nrxx); @@ -376,7 +396,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) delete[] rho_mag; delete[] rho_mag_save; } - else if (PARAM.inp.nspin == 4 && PARAM.inp.mixing_angle <= 0) + else if (nspin == 4 && PARAM.inp.mixing_angle <= 0) { // normal broyden mixing for {rho, mx, my, mz} rhor_in = chr->rho_save[0]; @@ -407,7 +427,7 @@ void Charge_Mixing::mix_rho_real(Charge* chr) this->mixing->cal_coef(this->rho_mdata, inner_product); this->mixing->mix_data(this->rho_mdata, rhor_out); } - else if (PARAM.inp.nspin == 4 && PARAM.inp.mixing_angle > 0) + else if (nspin == 4 && PARAM.inp.mixing_angle > 0) { // special broyden mixing for {rho, |m|} proposed by J. Phys. Soc. Jpn. 82 (2013) 114706 // here only consider the case of mixing_angle = 1, which mean only change |m| and keep angle fixed @@ -494,6 +514,8 @@ void Charge_Mixing::mix_rho_real(Charge* chr) this->mixing->mix_data(this->tau_mdata, taur_out); } + ModuleBase::timer::tick("Charge_Mixing", "mix_rho_real"); + return; } @@ -502,10 +524,13 @@ void Charge_Mixing::mix_rho(Charge* chr) ModuleBase::TITLE("Charge_Mixing", "mix_rho"); ModuleBase::timer::tick("Charge_Mixing", "mix_rho"); + const int nspin = PARAM.inp.nspin; + assert(nspin==1 || nspin==2 || nspin==4); + // the charge before mixing. const int nrxx = chr->rhopw->nrxx; - std::vector rho123(PARAM.inp.nspin * nrxx); - for (int is = 0; is < PARAM.inp.nspin; ++is) + std::vector rho123(nspin * nrxx); + for (int is = 0; is < nspin; ++is) { if (is == 0 || is == 3 || !PARAM.globalv.domag_z) { @@ -522,8 +547,8 @@ void Charge_Mixing::mix_rho(Charge* chr) std::vector kin_r123; if ((XC_Functional::get_ked_flag()) && mixing_tau) { - kin_r123.resize(PARAM.inp.nspin * nrxx); - for (int is = 0; is < PARAM.inp.nspin; ++is) + kin_r123.resize(nspin * nrxx); + for (int is = 0; is < nspin; ++is) { double* kin_r123_is = kin_r123.data() + is * nrxx; #ifdef _OPENMP @@ -548,7 +573,7 @@ void Charge_Mixing::mix_rho(Charge* chr) // mohan add 2012-06-05 // rho_save is the charge before mixing - for (int is = 0; is < PARAM.inp.nspin; ++is) + for (int is = 0; is < nspin; ++is) { if (is == 0 || is == 3 || !PARAM.globalv.domag_z) { @@ -565,7 +590,7 @@ void Charge_Mixing::mix_rho(Charge* chr) if ((XC_Functional::get_ked_flag()) && mixing_tau) { - for (int is = 0; is < PARAM.inp.nspin; ++is) + for (int is = 0; is < nspin; ++is) { double* kin_r123_is = kin_r123.data() + is * nrxx; #ifdef _OPENMP diff --git a/source/source_estate/module_charge/symmetry_rho.cpp b/source/source_estate/module_charge/symmetry_rho.cpp index f6bc1272f9..dbd8a57af1 100644 --- a/source/source_estate/module_charge/symmetry_rho.cpp +++ b/source/source_estate/module_charge/symmetry_rho.cpp @@ -11,37 +11,45 @@ Symmetry_rho::~Symmetry_rho() } void Symmetry_rho::begin(const int& spin_now, - const Charge& CHR, + const Charge& chr, const ModulePW::PW_Basis* rho_basis, ModuleSymmetry::Symmetry& symm) const { assert(spin_now < 4); // added by zhengdy-soc - if (ModuleSymmetry::Symmetry::symm_flag != 1) { - return; -} + if (ModuleSymmetry::Symmetry::symm_flag != 1) + { + return; + } + + ModuleBase::TITLE("Symmetry_rho", "begin"); + ModuleBase::timer::tick("Symmetry_rho","begin"); + // both parallel and serial // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space // { -// psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm); -// if(XC_Functional::get_ked_flag()) psymm(CHR.kin_r[spin_now], +// psymm(chr.rho[spin_now], rho_basis, Pgrid, symm); +// if(XC_Functional::get_ked_flag()) psymm(chr.kin_r[spin_now], // rho_basis,Pgrid,symm); // } // else //space group, do rho_symm in reciprocal space -{ - rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]); - psymmg(CHR.rhog[spin_now], rho_basis, symm); // need to modify - rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]); - if (XC_Functional::get_ked_flag() || CHR.cal_elf) + rho_basis->real2recip(chr.rho[spin_now], chr.rhog[spin_now]); + + psymmg(chr.rhog[spin_now], rho_basis, symm); // need to modify + + rho_basis->recip2real(chr.rhog[spin_now], chr.rho[spin_now]); + + if (XC_Functional::get_ked_flag() || chr.cal_elf) { // Use std::vector to manage kin_g instead of raw pointer - std::vector> kin_g(CHR.ngmc); - rho_basis->real2recip(CHR.kin_r[spin_now], kin_g.data()); + std::vector> kin_g(chr.ngmc); + rho_basis->real2recip(chr.kin_r[spin_now], kin_g.data()); psymmg(kin_g.data(), rho_basis, symm); - rho_basis->recip2real(kin_g.data(), CHR.kin_r[spin_now]); - } + rho_basis->recip2real(kin_g.data(), chr.kin_r[spin_now]); } + + ModuleBase::timer::tick("Symmetry_rho","begin"); return; } @@ -59,6 +67,10 @@ void Symmetry_rho::begin(const int& spin_now, { return; } + + ModuleBase::TITLE("Symmetry_rho", "begin"); + ModuleBase::timer::tick("Symmetry_rho","begin"); + // both parallel and serial // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space // { @@ -81,6 +93,8 @@ void Symmetry_rho::begin(const int& spin_now, rho_basis->recip2real(kin_g.data(), kin_r[spin_now]); } } + + ModuleBase::timer::tick("Symmetry_rho","begin"); return; } @@ -89,8 +103,11 @@ void Symmetry_rho::psymm(double* rho_part, Parallel_Grid& Pgrid, ModuleSymmetry::Symmetry& symm) const { + ModuleBase::TITLE("Symmetry_rho", "psymm"); + ModuleBase::timer::tick("Symmetry_rho","psymm"); + #ifdef __MPI - // (1) reduce all rho from the first pool. + // reduce all rho from the first pool. std::vector rhotot; if (GlobalV::MY_RANK == 0) { @@ -99,7 +116,6 @@ void Symmetry_rho::psymm(double* rho_part, } Pgrid.reduce(rhotot.data(), rho_part, false); - // (2) if (GlobalV::MY_RANK == 0) { symm.rho_symmetry(rhotot.data(), rho_basis->nx, rho_basis->ny, rho_basis->nz); @@ -126,8 +142,9 @@ void Symmetry_rho::psymm(double* rho_part, #ifdef __MPI } - // (3) -Pgrid.bcast(rhotot.data(), rho_part, GlobalV::MY_RANK); + Pgrid.bcast(rhotot.data(), rho_part, GlobalV::MY_RANK); #endif + + ModuleBase::timer::tick("Symmetry_rho","psymm"); return; -} \ No newline at end of file +} diff --git a/source/source_estate/module_pot/potential_new.cpp b/source/source_estate/module_pot/potential_new.cpp index 3d958bc87e..b1968339aa 100644 --- a/source/source_estate/module_pot/potential_new.cpp +++ b/source/source_estate/module_pot/potential_new.cpp @@ -23,7 +23,8 @@ Potential::Potential(const ModulePW::PW_Basis* rho_basis_in, double* etxc_in, double* vtxc_in, VSep* vsep_cell_in) - : ucell_(ucell_in), vloc_(vloc_in), structure_factors_(structure_factors_in), solvent_(solvent_in), vsep_cell(vsep_cell_in), etxc_(etxc_in), + : ucell_(ucell_in), vloc_(vloc_in), structure_factors_(structure_factors_in), + solvent_(solvent_in), vsep_cell(vsep_cell_in), etxc_(etxc_in), vtxc_(vtxc_in) { this->rho_basis_ = rho_basis_in; @@ -81,7 +82,6 @@ void Potential::pot_register(const std::vector& components_list) { PotBase* tmp = this->get_pot_type(comp); this->components.push_back(tmp); - // GlobalV::ofs_running << "Successful completion of Potential's registration : " << comp << std::endl; } // after register, reset fixed_done to false @@ -93,8 +93,13 @@ void Potential::pot_register(const std::vector& components_list) void Potential::allocate() { ModuleBase::TITLE("Potential", "allocate"); - int nrxx = this->rho_basis_->nrxx; - int nrxx_smooth = this->rho_basis_smooth_->nrxx; + + const int nspin = PARAM.inp.nspin; + assert(nspin==1 || nspin==2 || nspin==4); + + const int nrxx = this->rho_basis_->nrxx; + const int nrxx_smooth = this->rho_basis_smooth_->nrxx; + if (nrxx == 0) { return; @@ -107,39 +112,39 @@ void Potential::allocate() this->v_effective_fixed.resize(nrxx); ModuleBase::Memory::record("Pot::veff_fix", sizeof(double) * nrxx); - this->v_effective.create(PARAM.inp.nspin, nrxx); - ModuleBase::Memory::record("Pot::veff", sizeof(double) * PARAM.inp.nspin * nrxx); + this->v_effective.create(nspin, nrxx); + ModuleBase::Memory::record("Pot::veff", sizeof(double) * nspin * nrxx); - this->veff_smooth.create(PARAM.inp.nspin, nrxx_smooth); - ModuleBase::Memory::record("Pot::veff_smooth", sizeof(double) * PARAM.inp.nspin * nrxx_smooth); + this->veff_smooth.create(nspin, nrxx_smooth); + ModuleBase::Memory::record("Pot::veff_smooth", sizeof(double) * nspin * nrxx_smooth); if (XC_Functional::get_ked_flag()) { - this->vofk_effective.create(PARAM.inp.nspin, nrxx); - ModuleBase::Memory::record("Pot::vofk", sizeof(double) * PARAM.inp.nspin * nrxx); + this->vofk_effective.create(nspin, nrxx); + ModuleBase::Memory::record("Pot::vofk", sizeof(double) * nspin * nrxx); - this->vofk_smooth.create(PARAM.inp.nspin, nrxx_smooth); - ModuleBase::Memory::record("Pot::vofk_smooth", sizeof(double) * PARAM.inp.nspin * nrxx_smooth); + this->vofk_smooth.create(nspin, nrxx_smooth); + ModuleBase::Memory::record("Pot::vofk_smooth", sizeof(double) * nspin * nrxx_smooth); } if (use_gpu_) { if (PARAM.globalv.has_float_data) { - resmem_sd_op()(s_veff_smooth, PARAM.inp.nspin * nrxx_smooth); - resmem_sd_op()(s_vofk_smooth, PARAM.inp.nspin * nrxx_smooth); + resmem_sd_op()(s_veff_smooth, nspin * nrxx_smooth); + resmem_sd_op()(s_vofk_smooth, nspin * nrxx_smooth); } if (PARAM.globalv.has_double_data) { - resmem_dd_op()(d_veff_smooth, PARAM.inp.nspin * nrxx_smooth); - resmem_dd_op()(d_vofk_smooth, PARAM.inp.nspin * nrxx_smooth); + resmem_dd_op()(d_veff_smooth, nspin * nrxx_smooth); + resmem_dd_op()(d_vofk_smooth, nspin * nrxx_smooth); } } else { if (PARAM.globalv.has_float_data) { - resmem_sh_op()(s_veff_smooth, PARAM.inp.nspin * nrxx_smooth, "POT::sveff_smooth"); - resmem_sh_op()(s_vofk_smooth, PARAM.inp.nspin * nrxx_smooth, "POT::svofk_smooth"); + resmem_sh_op()(s_veff_smooth, nspin * nrxx_smooth, "POT::sveff_smooth"); + resmem_sh_op()(s_vofk_smooth, nspin * nrxx_smooth, "POT::svofk_smooth"); } if (PARAM.globalv.has_double_data) { @@ -273,20 +278,23 @@ void Potential::get_vnew(const Charge* chg, ModuleBase::matrix& vnew) return; } -void Potential::interpolate_vrs() +void Potential::interpolate_vrs(void) { ModuleBase::TITLE("Potential", "interpolate_vrs"); ModuleBase::timer::tick("Potential", "interpolate_vrs"); - if ( PARAM.globalv.double_grid) + const int nspin = PARAM.inp.nspin; + assert(nspin==1 || nspin==2 || nspin==4); + + if (PARAM.globalv.double_grid) { if (rho_basis_->gamma_only != rho_basis_smooth_->gamma_only) { ModuleBase::WARNING_QUIT("Potential::interpolate_vrs", "gamma_only is not consistent"); } - ModuleBase::ComplexMatrix vrs(PARAM.inp.nspin, rho_basis_->npw); - for (int is = 0; is < PARAM.inp.nspin; is++) + ModuleBase::ComplexMatrix vrs(nspin, rho_basis_->npw); + for (int is = 0; is < nspin; is++) { rho_basis_->real2recip(&v_effective(is, 0), &vrs(is, 0)); rho_basis_smooth_->recip2real(&vrs(is, 0), &veff_smooth(is, 0)); @@ -294,8 +302,8 @@ void Potential::interpolate_vrs() if (XC_Functional::get_ked_flag()) { - ModuleBase::ComplexMatrix vrs_ofk(PARAM.inp.nspin, rho_basis_->npw); - for (int is = 0; is < PARAM.inp.nspin; is++) + ModuleBase::ComplexMatrix vrs_ofk(nspin, rho_basis_->npw); + for (int is = 0; is < nspin; is++) { rho_basis_->real2recip(&vofk_effective(is, 0), &vrs_ofk(is, 0)); rho_basis_smooth_->recip2real(&vrs_ofk(is, 0), &vofk_smooth(is, 0));