@@ -310,21 +310,19 @@ void Charge_Mixing::mix_rho_recip(Charge* chr)
310310 return ;
311311}
312312
313- void Charge_Mixing::mix_rho_real (Charge* chr)
313+
314+ template <typename TCharge>
315+ void Charge_Mixing::mix_rho_real (TCharge*const chr)
314316{
315317 ModuleBase::TITLE (" Charge_Mixing" , " mix_rho_real" );
316318 ModuleBase::timer::tick (" Charge_Mixing" , " mix_rho_real" );
317319
318320 const int nspin = PARAM .inp .nspin ;
319321 assert (nspin==1 || nspin==2 || nspin==4 );
320-
321- double * rhor_in=nullptr ;
322- double * rhor_out=nullptr ;
323-
324322 if (nspin == 1 )
325323 {
326- rhor_in = chr->rho_save [0 ];
327- rhor_out = chr->rho [0 ];
324+ double * rhor_in = chr->rho_save [0 ];
325+ double * rhor_out = chr->rho [0 ];
328326 auto screen = std::bind (&Charge_Mixing::Kerker_screen_real, this , std::placeholders::_1);
329327 this ->mixing ->push_data (this ->rho_mdata , rhor_in, rhor_out, screen, true );
330328 auto inner_product
@@ -355,22 +353,23 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
355353 rho_mag_save[ir + nrxx] = chr->rho_save [0 ][ir] - chr->rho_save [1 ][ir];
356354 }
357355 //
358- rhor_in = rho_mag_save;
359- rhor_out = rho_mag;
356+ double * rhor_in = rho_mag_save;
357+ double * rhor_out = rho_mag;
360358 auto screen = std::bind (&Charge_Mixing::Kerker_screen_real, this , std::placeholders::_1);
361359 auto twobeta_mix
362- = [this , nrxx](double * out, const double * in, const double * sres) {
363- #ifdef _OPENMP
364- #pragma omp parallel for schedule(static, 256)
365- #endif
360+ = [this , nrxx](double * out, const double * in, const double * sres)
361+ {
362+ #ifdef _OPENMP
363+ #pragma omp parallel for schedule(static, 256)
364+ #endif
366365 for (int i = 0 ; i < nrxx; ++i)
367366 {
368367 out[i] = in[i] + this ->mixing_beta * sres[i];
369368 }
370369 // magnetism
371- #ifdef _OPENMP
372- #pragma omp parallel for schedule(static, 256)
373- #endif
370+ #ifdef _OPENMP
371+ #pragma omp parallel for schedule(static, 256)
372+ #endif
374373 for (int i = nrxx; i < 2 * nrxx; ++i)
375374 {
376375 out[i] = in[i] + this ->mixing_beta_mag * sres[i];
@@ -399,23 +398,24 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
399398 else if (nspin == 4 && PARAM .inp .mixing_angle <= 0 )
400399 {
401400 // normal broyden mixing for {rho, mx, my, mz}
402- rhor_in = chr->rho_save [0 ];
403- rhor_out = chr->rho [0 ];
401+ double * rhor_in = chr->rho_save [0 ];
402+ double * rhor_out = chr->rho [0 ];
404403 const int nrxx = this ->rhopw ->nrxx ;
405404 auto screen = std::bind (&Charge_Mixing::Kerker_screen_real, this , std::placeholders::_1);
406405 auto twobeta_mix
407- = [this , nrxx](double * out, const double * in, const double * sres) {
408- #ifdef _OPENMP
409- #pragma omp parallel for schedule(static, 256)
410- #endif
406+ = [this , nrxx](double * out, const double * in, const double * sres)
407+ {
408+ #ifdef _OPENMP
409+ #pragma omp parallel for schedule(static, 256)
410+ #endif
411411 for (int i = 0 ; i < nrxx; ++i)
412412 {
413413 out[i] = in[i] + this ->mixing_beta * sres[i];
414414 }
415415 // magnetism, mx, my, mz
416- #ifdef _OPENMP
417- #pragma omp parallel for schedule(static, 256)
418- #endif
416+ #ifdef _OPENMP
417+ #pragma omp parallel for schedule(static, 256)
418+ #endif
419419 for (int i = nrxx; i < 4 * nrxx; ++i)
420420 {
421421 out[i] = in[i] + this ->mixing_beta_mag * sres[i];
@@ -451,23 +451,24 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
451451 + chr->rho_save [2 ][ir] * chr->rho_save [2 ][ir]
452452 + chr->rho_save [3 ][ir] * chr->rho_save [3 ][ir]);
453453 }
454- rhor_in = rho_magabs_save;
455- rhor_out = rho_magabs;
454+ double * rhor_in = rho_magabs_save;
455+ double * rhor_out = rho_magabs;
456456
457457 auto screen = std::bind (&Charge_Mixing::Kerker_screen_real, this , std::placeholders::_1);
458458 auto twobeta_mix
459- = [this , nrxx](double * out, const double * in, const double * sres) {
460- #ifdef _OPENMP
461- #pragma omp parallel for schedule(static, 256)
462- #endif
459+ = [this , nrxx](double * out, const double * in, const double * sres)
460+ {
461+ #ifdef _OPENMP
462+ #pragma omp parallel for schedule(static, 256)
463+ #endif
463464 for (int i = 0 ; i < nrxx; ++i)
464465 {
465466 out[i] = in[i] + this ->mixing_beta * sres[i];
466467 }
467468 // magnetism, |m|
468- #ifdef _OPENMP
469- #pragma omp parallel for schedule(static, 256)
470- #endif
469+ #ifdef _OPENMP
470+ #pragma omp parallel for schedule(static, 256)
471+ #endif
471472 for (int i = nrxx; i < 2 * nrxx; ++i)
472473 {
473474 out[i] = in[i] + this ->mixing_beta_mag * sres[i];
@@ -500,22 +501,28 @@ void Charge_Mixing::mix_rho_real(Charge* chr)
500501 delete[] rho_magabs;
501502 delete[] rho_magabs_save;
502503 }
503-
504- double *taur_out=nullptr ;
505- double *taur_in=nullptr ;
506- if ((XC_Functional::get_ked_flag ()) && mixing_tau)
504+ else
507505 {
508- taur_in = chr->kin_r_save [0 ];
509- taur_out = chr->kin_r [0 ];
510- // Note: there is no kerker modification for tau because I'm not sure
511- // if we should have it. If necessary we can try it in the future.
512- this ->mixing ->push_data (this ->tau_mdata , taur_in, taur_out, nullptr , false );
513-
514- this ->mixing ->mix_data (this ->tau_mdata , taur_out);
506+ ModuleBase::WARNING_QUIT (" Charge_Mixing::mix_rho_real" , std::string (__FILE__)+" line " +std::to_string (__LINE__));
515507 }
516508
517509 ModuleBase::timer::tick (" Charge_Mixing" , " mix_rho_real" );
518- return ;
510+ }
511+
512+ template void Charge_Mixing::mix_rho_real (Charge*const chr);
513+
514+
515+ void Charge_Mixing::mix_rho_real_tau (Charge*const chr)
516+ {
517+ ModuleBase::TITLE (" Charge_Mixing" , " mix_rho_real_tau" );
518+ ModuleBase::timer::tick (" Charge_Mixing" , " mix_rho_real_tau" );
519+ double *taur_out=chr->kin_r [0 ];
520+ double *taur_in=chr->kin_r_save [0 ];
521+ // Note: there is no kerker modification for tau because I'm not sure
522+ // if we should have it. If necessary we can try it in the future.
523+ this ->mixing ->push_data (this ->tau_mdata , taur_in, taur_out, nullptr , false );
524+ this ->mixing ->mix_data (this ->tau_mdata , taur_out);
525+ ModuleBase::timer::tick (" Charge_Mixing" , " mix_rho_real_tau" );
519526}
520527
521528
@@ -568,6 +575,10 @@ void Charge_Mixing::mix_rho(Charge* chr)
568575 else if (PARAM .inp .scf_thr_type == 2 )
569576 {
570577 mix_rho_real (chr);
578+ if ((XC_Functional::get_ked_flag ()) && mixing_tau)
579+ {
580+ mix_rho_real_tau (chr);
581+ }
571582 }
572583 // ---------------------------------------------------
573584
@@ -610,4 +621,4 @@ void Charge_Mixing::mix_rho(Charge* chr)
610621
611622 ModuleBase::timer::tick (" Charge_Mixing" , " mix_rho" );
612623 return ;
613- }
624+ }
0 commit comments