@@ -100,6 +100,7 @@ void RDMFT<TK, TR>::init(Gint_Gamma& GG_in, Gint_k& GK_in, Parallel_Orbitals& Pa
100100 nk_total = ModuleSymmetry::Symmetry::symm_flag == -1 ? kv->nkstot_full : kv->nks ;
101101 nbands_total = PARAM.inp .nbands ;
102102 nspin = PARAM.inp .nspin ;
103+ only_exx_type = ( XC_func_rdmft == " hf" || XC_func_rdmft == " muller" || XC_func_rdmft == " power" );
103104
104105 // XC_func_rdmft = "power"; // just for test
105106 // alpha_power = 0.525;
@@ -123,7 +124,6 @@ void RDMFT<TK, TR>::init(Gint_Gamma& GG_in, Gint_k& GK_in, Parallel_Orbitals& Pa
123124 // para_Eij.set_desc( GlobalV::NBANDS, GlobalV::NBANDS, para_Eij.get_row_size(), false );
124125 para_Eij.set (nbands_total, nbands_total, ParaV->nb , ParaV->blacs_ctxt ); // maybe in default, PARAM.inp.nb2d = 0, can't be used
125126 // para_Eij.init(nbands_total, nbands_total, PARAM.inp.nb2d, MPI_COMM_WORLD);
126-
127127 // // learn from "module_hamilt_lcao/hamilt_lcaodft/LCAO_init_basis.cpp"
128128
129129 //
@@ -274,7 +274,7 @@ void RDMFT<TK, TR>::update_elec(const ModuleBase::matrix& occ_number_in, const p
274274 this ->update_charge ();
275275
276276 // "default" = "pbe"
277- // if( (XC_func_rdmft != "hf" && XC_func_rdmft != "muller" && XC_func_rdmft != "power") || this->cal_E_type != 1 )
277+ // if( !only_exx_type || this->cal_E_type != 1 )
278278 if ( this ->cal_E_type != 1 )
279279 {
280280 // the second cal_E_type need the complete pot to get effctive_V to calEband and so on.
@@ -577,30 +577,30 @@ void RDMFT<TK, TR>::cal_V_XC()
577577
578578 DM_XC_pass = DM_XC;
579579
580- // elecstate::DensityMatrix<TK, double> DM_test(kv, ParaV, nspin);
581- // elecstate::cal_dm_psi(ParaV, wg, wfc, DM_test);
582- // DM_test.init_DMR(&GlobalC::GridD, &GlobalC::ucell);
583- // DM_test.cal_DMR();
580+ elecstate::DensityMatrix<TK, double > DM_test (kv, ParaV, nspin);
581+ elecstate::cal_dm_psi (ParaV, wg, wfc, DM_test);
582+ DM_test.init_DMR (&GlobalC::GridD, &GlobalC::ucell);
583+ DM_test.cal_DMR ();
584584
585- // // compare DM_XC and DM get in update_charge(or ABACUS)
586- // std::cout << "\n\ntest DM_XC - DM in ABACUS: \n" << std::endl;
587- // double DM_XC_minus_DMtest = 0.0;
588- // for(int ik=0; ik<nk_total; ++ik)
589- // {
590- // TK* dmk_pointer = DM_test.get_DMK_pointer(ik);
591- // for(int iloc=0; iloc<ParaV->nloc; ++iloc)
592- // {
593- // double test = std::abs(DM_XC[ik][iloc] - dmk_pointer[iloc]);
594- // DM_XC_minus_DMtest += test;
595- // if( test > 1e-16 )
596- // {
597- // std::cout << "\nik, iloc, minus[ik][iloc]: " << ik << " " << iloc << " " << test << std::endl;
598- // }
599- // }
600- // }
601- // std::cout << "\nsum of DM_XC - DM in ABACUS: " << DM_XC_minus_DMtest << std::endl;
585+ // compare DM_XC and DM get in update_charge(or ABACUS)
586+ std::cout << " \n\n test DM_XC - DM in ABACUS: \n " << std::endl;
587+ double DM_XC_minus_DMtest = 0.0 ;
588+ for (int ik=0 ; ik<nk_total; ++ik)
589+ {
590+ TK* dmk_pointer = DM_test.get_DMK_pointer (ik);
591+ for (int iloc=0 ; iloc<ParaV->nloc ; ++iloc)
592+ {
593+ double test = std::abs (DM_XC[ik][iloc] - dmk_pointer[iloc]);
594+ DM_XC_minus_DMtest += test;
595+ if ( test > 1e-16 )
596+ {
597+ std::cout << " \n ik, iloc, minus[ik][iloc]: " << ik << " " << iloc << " " << test << std::endl;
598+ }
599+ }
600+ }
601+ std::cout << " \n sum of DM_XC - DM in ABACUS: " << DM_XC_minus_DMtest << std::endl;
602602
603- if ( XC_func_rdmft != " hf " && XC_func_rdmft != " muller " && XC_func_rdmft != " power " )
603+ if ( !only_exx_type )
604604 {
605605 if ( PARAM.inp .gamma_only )
606606 {
@@ -623,7 +623,6 @@ void RDMFT<TK, TR>::cal_V_XC()
623623 &etxc,
624624 &vtxc
625625 );
626- // V_dft_XC->contributeHR();
627626 }
628627 else
629628 {
@@ -646,7 +645,6 @@ void RDMFT<TK, TR>::cal_V_XC()
646645 &etxc,
647646 &vtxc
648647 );
649- // V_dft_XC->contributeHR();
650648 }
651649 V_dft_XC->contributeHR ();
652650 }
@@ -670,8 +668,16 @@ void RDMFT<TK, TR>::cal_V_XC()
670668 hsk_exx_XC,
671669 HR_exx_XC,
672670 *kv,
673- &Vxc_fromRI_d->Hexxs
671+ &Vxc_fromRI_d->Hexxs ,
672+ nullptr ,
673+ hamilt::Add_Hexx_Type::k
674674 );
675+ // V_exx_XC = new hamilt::OperatorEXX<hamilt::OperatorLCAO<TK, TR>>(
676+ // hsk_exx_XC,
677+ // HR_exx_XC,
678+ // *kv
679+ // );
680+
675681 }
676682 else
677683 {
@@ -691,10 +697,12 @@ void RDMFT<TK, TR>::cal_V_XC()
691697 HR_exx_XC,
692698 *kv,
693699 nullptr ,
694- &Vxc_fromRI_c->Hexxs
700+ &Vxc_fromRI_c->Hexxs ,
701+ hamilt::Add_Hexx_Type::k
695702 );
696703 }
697- V_exx_XC->contributeHR ();
704+ // use hamilt::Add_Hexx_Type::k, not R, contributeHR() should be skipped
705+ // V_exx_XC->contributeHR();
698706 }
699707}
700708
@@ -740,7 +748,7 @@ void RDMFT<TK, TR>::cal_Hk_Hpsi()
740748 for (int iloc=0 ; iloc<HK_XC.size (); ++iloc) HK_XC[iloc] += hsk_exx_XC->get_hk ()[iloc];
741749 }
742750
743- if ( XC_func_rdmft != " hf " && XC_func_rdmft != " muller " && XC_func_rdmft != " power " )
751+ if ( !only_exx_type )
744752 {
745753 // set_zero_vector(HK_dft_XC);
746754 hsk_dft_XC->set_zero_hk ();
0 commit comments