@@ -261,12 +261,13 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(const Input_para& inp, UnitCell
261261 }
262262
263263 // add by jghan for rdmft calculation
264- if ( GlobalV::CALCULATION == " rdmft" || true )
264+ // if( PARAM.inp.ab_initio_type == "rdmft" )
265+ if ( 1 )
265266 {
266267 // rdmft_solver.init( this->UHM.GG, this->UHM.GK, this->orb_con.ParaV, ucell, this->kv, *(this->pelec),
267268 // GlobalV::DFT_FUNCTIONAL, GlobalV::rdmft_power_alpha);
268- rdmft_solver.init ( this ->UHM . GG , this ->UHM . GK , this ->orb_con . ParaV , ucell, this ->kv , *(this ->pelec ),
269- GlobalV::DFT_FUNCTIONAL , 1.0 );
269+ rdmft_solver.init ( this ->GG , this ->GK , &( this ->pv ) , ucell, this ->kv , *(this ->pelec ),
270+ this -> orb_ , PARAM. inp . dft_functional , 1.0 );
270271
271272 // the initialization and necessary calculations of these quantities have been completed in init()
272273 // rdmft_solver.update_ion(ucell, LM, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac, this->LOC);
@@ -1209,149 +1210,149 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
12091210#endif
12101211
12111212
1212- /* ******* test RDMFT *********/
1213-
1214- // esolver_ks_lcao.h(LCAO_Matrix LM), LCAO_matrix.h(Parallel_Orbitals* ParaV)
1215- // esolver_fp.h(elecstate::ElecState* pelec), elecstate.h(ModuleBase::matrix wg), this->pelec->wg
1216- // esolver_fp.h(elecstate::ElecState* pelec), elecstate.h(Charge* charge), this->pelec->wg
1217- // esolver_ks.h(psi::Psi<T>* psi)
1218- // esolver_fp.h(K_Vectors kv), kv.kvec_d
1219- // esolver_ks_lcao.h(Local_Orbital_Charge LOC)
1220- // esolver_ks_lcao.h(LCAO_Hamilt UHM), LCAO_hamilt.h(Gint_Gamma GG, Gint_k GK), this->UHM.GK
1221- // used for gamma only algorithms, Gint_Gamma GG. for k-dependent grid integration, Gint_k GK;
1222- // esolver_fp.h(ModulePW::PW_Basis* pw_rho, ModulePW::PW_Basis* pw_rhod), this->pw_rho
1223- // GlobalC::ppcell.vloc(ModuleBase::matrix vloc),
1224- // esolver_fp.h(Structure_Factor sf), structure_factor.h(ModuleBase::ComplexMatrix strucFac), this->sf.strucFac
1225- // elecstate::DensityMatrix<TK, double>& *( dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM() )
1226- // elecstate.h (Potential* pot, this->pelec->pot)
1227-
1228- // initialize the gradients of Etotal on wg and wfc, and set all elements to 0.
1229- ModuleBase::matrix E_gradient_occNum (this ->pelec ->wg .nr , this ->pelec ->wg .nc , true );
1230- psi::Psi<TK> E_gradient_wfc (this ->psi ->get_nk (), this ->psi ->get_nbands (), this ->psi ->get_nbasis ());
1231- E_gradient_wfc.zero_out ();
1232- double Etotal_RDMFT = 0.0 ;
1233-
1234- // get natural occupation numbers from wg which considers k point weights and spin, this just proper for nspin=1 !!! or 2 ? in Soild Si, it's correct
1235- // wk consider both weight of k-point and spin. When nspin=1, wk[ik] = W_k * 2 . When nspin=2, wk[ik] = W_k
1236- ModuleBase::matrix occ_number (this ->pelec ->wg );
1237- for (int ik=0 ; ik < occ_number.nr ; ++ik)
1238- {
1239- for (int inb=0 ; inb < this ->pelec ->wg .nc ; ++inb) occ_number (ik, inb) /= this ->kv .wk [ik];
1240- }
1241-
1242- // gamma only calculation
1243- if ( GlobalV::GAMMA_ONLY_LOCAL )
1244- {
1245- // when gamma_only, this->psi->get_nbands() = nbasis_local, so we use ParaV get the true nbands_local for wfc
1246- const int nbands_local = LM.ParaV ->ncol_bands ;
1247- const int nk_total = this ->psi ->get_nk ();
1248- const int nbasis_local = this ->psi ->get_nbasis ();
1249-
1250- std::cout << " \n\n\n nk_total= " << nk_total << " \n\n\n " ;
1251-
1252- psi::Psi<TK> wfc_rdmft (nk_total, nbands_local, nbasis_local);
1253- for (int ik=0 ; ik<nk_total; ++ik)
1254- {
1255- for (int inbn=0 ; inbn<nbands_local; ++inbn)
1256- {
1257- for (int inbs=0 ; inbs<nbasis_local; ++inbs) wfc_rdmft (ik, inbn, inbs) = (*(this ->psi ))(ik, inbn, inbs);
1258- }
1259- }
1260-
1261- Etotal_RDMFT = rdmftTest::rdmft_cal<TK,TR,Gint_Gamma>(
1262- &LM,
1263- LM.ParaV ,
1264- occ_number,
1265- this ->pelec ->wg ,
1266- wfc_rdmft,
1267- E_gradient_occNum,
1268- E_gradient_wfc,
1269- this ->kv ,
1270- this ->UHM .GG ,
1271- this ->LOC ,
1272- *(this ->pelec ->charge ),
1273- *(this ->pw_rho ),
1274- GlobalC::ppcell.vloc ,
1275- this ->sf .strucFac ,
1276- " HF" ,
1277- 1.0
1278- );
1279- }
1280- // multi-k calculation
1281- else
1282- {
1283- Etotal_RDMFT = rdmftTest::rdmft_cal<TK,TR,Gint_k>(
1284- &LM,
1285- LM.ParaV ,
1286- occ_number,
1287- this ->pelec ->wg ,
1288- *(this ->psi ),
1289- E_gradient_occNum,
1290- E_gradient_wfc,
1291- this ->kv ,
1292- this ->UHM .GK ,
1293- this ->LOC ,
1294- *(this ->pelec ->charge ),
1295- *(this ->pw_rho ),
1296- GlobalC::ppcell.vloc ,
1297- this ->sf .strucFac ,
1298- " HF" ,
1299- 1.0
1300- );
1301- }
1302-
1303- // GlobalV::ofs_running << "E_one_elec "
1304- // << std::setw(10) << pelec->f_en.deband + pelec->f_en.eband << std::endl;
1305- // GlobalV::ofs_running << "E_Hartree "
1306- // << std::setw(10) << pelec->f_en.hartree_energy << std::endl;
1307- // GlobalV::ofs_running << "E_xc "
1308- // << std::setw(10) << pelec->f_en.etxc - pelec->f_en.etxcc<< std::endl;
1309-
1310- double E_xc_KS = this ->pelec ->f_en .etxc - this ->pelec ->f_en .etxcc ;
1311-
1312- double E_TV_Hartree_EXX = this ->pelec ->f_en .deband + this ->pelec ->f_en .eband
1313- + this ->pelec ->f_en .hartree_energy + this ->pelec ->f_en .exx ;
1314- double E_TV_Hartree_XC = this ->pelec ->f_en .deband + this ->pelec ->f_en .eband
1315- + this ->pelec ->f_en .hartree_energy + this ->pelec ->f_en .etxc - this ->pelec ->f_en .etxcc ;
1316- std::cout << std::fixed << std::setprecision (10 ) << " \n\n ******\n E(one_elec + Hartree + EXX) by KS-DFT: " << E_TV_Hartree_EXX
1317- << " \n\n E(one_elec + Hartree + XC) by KS-DFT: " << E_TV_Hartree_XC << " \n Exc_KS: " << E_xc_KS << " \n ******\n\n " << std::endl;
1318-
1319- ModuleBase::TITLE (" RDMFT" , " E & Egradient" );
1320- ModuleBase::timer::tick (" RDMFT" , " E & Egradient" );
1321-
1322- // // test class rdmft
1323- // rdmft_solver.update_elec(occ_number, *(this->psi));
1324- // this->Run_rdmft(E_gradient_occNum, E_gradient_wfc);
1325- // std::cout << "\nrdmft_solver: " << "0.0000" << std::endl;
1326-
1327- // double minus = 0.0;
1328-
1329- // // delete in the future
1330- // std::cout << "\n\n occ_number(ik, inband): " << std::endl;
1213+ // /******** test RDMFT *********/
1214+
1215+ // // esolver_ks_lcao.h(LCAO_Matrix LM), LCAO_matrix.h(Parallel_Orbitals* ParaV)
1216+ // // esolver_fp.h(elecstate::ElecState* pelec), elecstate.h(ModuleBase::matrix wg), this->pelec->wg
1217+ // // esolver_fp.h(elecstate::ElecState* pelec), elecstate.h(Charge* charge), this->pelec->wg
1218+ // // esolver_ks.h(psi::Psi<T>* psi)
1219+ // // esolver_fp.h(K_Vectors kv), kv.kvec_d
1220+ // // esolver_ks_lcao.h(Local_Orbital_Charge LOC)
1221+ // // esolver_ks_lcao.h(LCAO_Hamilt UHM), LCAO_hamilt.h(Gint_Gamma GG, Gint_k GK), this->UHM.GK
1222+ // // used for gamma only algorithms, Gint_Gamma GG. for k-dependent grid integration, Gint_k GK;
1223+ // // esolver_fp.h(ModulePW::PW_Basis* pw_rho, ModulePW::PW_Basis* pw_rhod), this->pw_rho
1224+ // // GlobalC::ppcell.vloc(ModuleBase::matrix vloc),
1225+ // // esolver_fp.h(Structure_Factor sf), structure_factor.h(ModuleBase::ComplexMatrix strucFac), this->sf.strucFac
1226+ // // elecstate::DensityMatrix<TK, double>& *( dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM() )
1227+ // // elecstate.h (Potential* pot, this->pelec->pot)
1228+
1229+ // //initialize the gradients of Etotal on wg and wfc, and set all elements to 0.
1230+ // ModuleBase::matrix E_gradient_occNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
1231+ // psi::Psi<TK> E_gradient_wfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis());
1232+ // E_gradient_wfc.zero_out();
1233+ // double Etotal_RDMFT = 0.0;
1234+
1235+ // // get natural occupation numbers from wg which considers k point weights and spin, this just proper for nspin=1 !!! or 2 ? in Soild Si, it's correct
1236+ // // wk consider both weight of k-point and spin. When nspin=1, wk[ik] = W_k * 2 . When nspin=2, wk[ik] = W_k
1237+ // ModuleBase::matrix occ_number(this->pelec->wg);
13311238 // for(int ik=0; ik < occ_number.nr; ++ik)
13321239 // {
1333- // for(int inb=0; inb < occ_number.nc; ++inb) std::cout << occ_number(ik, inb) << " ";
1334- // std::cout << "\n" << std::endl;
1240+ // for(int inb=0; inb < this->pelec->wg.nc; ++inb) occ_number(ik, inb) /= this->kv.wk[ik];
13351241 // }
13361242
1337- // std::cout << "\n\n wg(ik, inband): " << std::endl;
1338- // for(int ik=0; ik < occ_number.nr; ++ik )
1243+ // // gamma only calculation
1244+ // if( GlobalV::GAMMA_ONLY_LOCAL )
13391245 // {
1340- // for(int inb=0; inb < occ_number.nc; ++inb)
1246+ // // when gamma_only, this->psi->get_nbands() = nbasis_local, so we use ParaV get the true nbands_local for wfc
1247+ // const int nbands_local = LM.ParaV->ncol_bands;
1248+ // const int nk_total = this->psi->get_nk();
1249+ // const int nbasis_local = this->psi->get_nbasis();
1250+
1251+ // std::cout << "\n\n\nnk_total= " << nk_total << "\n\n\n";
1252+
1253+ // psi::Psi<TK> wfc_rdmft(nk_total, nbands_local, nbasis_local);
1254+ // for(int ik=0; ik<nk_total; ++ik)
13411255 // {
1342- // std::cout << this->pelec->wg(ik, inb) << " ";
1343- // minus += this->pelec->wg(ik, inb) - occ_number(ik, inb)*this->kv.wk[ik];
1256+ // for(int inbn=0; inbn<nbands_local; ++inbn)
1257+ // {
1258+ // for(int inbs=0; inbs<nbasis_local; ++inbs) wfc_rdmft(ik, inbn, inbs) = (*(this->psi))(ik, inbn, inbs);
1259+ // }
13441260 // }
1345- // std::cout << "\n" << std::endl;
1261+
1262+ // Etotal_RDMFT = rdmftTest::rdmft_cal<TK,TR,Gint_Gamma>(
1263+ // &LM,
1264+ // LM.ParaV,
1265+ // occ_number,
1266+ // this->pelec->wg,
1267+ // wfc_rdmft,
1268+ // E_gradient_occNum,
1269+ // E_gradient_wfc,
1270+ // this->kv,
1271+ // this->UHM.GG,
1272+ // this->LOC,
1273+ // *(this->pelec->charge),
1274+ // *(this->pw_rho),
1275+ // GlobalC::ppcell.vloc,
1276+ // this->sf.strucFac,
1277+ // "HF",
1278+ // 1.0
1279+ // );
1280+ // }
1281+ // // multi-k calculation
1282+ // else
1283+ // {
1284+ // Etotal_RDMFT = rdmftTest::rdmft_cal<TK,TR,Gint_k>(
1285+ // &LM,
1286+ // LM.ParaV,
1287+ // occ_number,
1288+ // this->pelec->wg,
1289+ // *(this->psi),
1290+ // E_gradient_occNum,
1291+ // E_gradient_wfc,
1292+ // this->kv,
1293+ // this->UHM.GK,
1294+ // this->LOC,
1295+ // *(this->pelec->charge),
1296+ // *(this->pw_rho),
1297+ // GlobalC::ppcell.vloc,
1298+ // this->sf.strucFac,
1299+ // "HF",
1300+ // 1.0
1301+ // );
13461302 // }
13471303
1348- // if( std::abs(minus) < 1e-12 ) std::cout << "\n\nminus<1e-12, minus: " << minus << std::endl;
1349- // else if( std::abs(minus) < 1e-16 ) std::cout << "\n\nminus<1e-16, minus: " << minus << std::endl;
1350- // else {std::cout << "\n\nminus is big, minus: " << minus << std::endl;}
1351-
1352- ModuleBase::timer::tick (" RDMFT" , " E & Egradient" );
1353-
1354- /* ******* test RDMFT *********/
1304+ // // GlobalV::ofs_running << "E_one_elec "
1305+ // // << std::setw(10) << pelec->f_en.deband + pelec->f_en.eband << std::endl;
1306+ // // GlobalV::ofs_running << "E_Hartree "
1307+ // // << std::setw(10) << pelec->f_en.hartree_energy << std::endl;
1308+ // // GlobalV::ofs_running << "E_xc "
1309+ // // << std::setw(10) << pelec->f_en.etxc - pelec->f_en.etxcc<< std::endl;
1310+
1311+ // double E_xc_KS = this->pelec->f_en.etxc - this->pelec->f_en.etxcc;
1312+
1313+ // double E_TV_Hartree_EXX = this->pelec->f_en.deband + this->pelec->f_en.eband
1314+ // + this->pelec->f_en.hartree_energy + this->pelec->f_en.exx;
1315+ // double E_TV_Hartree_XC = this->pelec->f_en.deband + this->pelec->f_en.eband
1316+ // + this->pelec->f_en.hartree_energy + this->pelec->f_en.etxc - this->pelec->f_en.etxcc;
1317+ // std::cout << std::fixed << std::setprecision(10) << "\n\n******\nE(one_elec + Hartree + EXX) by KS-DFT: " << E_TV_Hartree_EXX
1318+ // << "\n\nE(one_elec + Hartree + XC) by KS-DFT: " << E_TV_Hartree_XC << "\nExc_KS: " << E_xc_KS << "\n******\n\n" << std::endl;
1319+
1320+ // ModuleBase::TITLE("RDMFT", "E & Egradient");
1321+ // ModuleBase::timer::tick("RDMFT", "E & Egradient");
1322+
1323+ // // // test class rdmft
1324+ // // rdmft_solver.update_elec(occ_number, *(this->psi));
1325+ // // this->Run_rdmft(E_gradient_occNum, E_gradient_wfc);
1326+ // // std::cout << "\nrdmft_solver: " << "0.0000" << std::endl;
1327+
1328+ // // double minus = 0.0;
1329+
1330+ // // // delete in the future
1331+ // // std::cout << "\n\n occ_number(ik, inband): " << std::endl;
1332+ // // for(int ik=0; ik < occ_number.nr; ++ik)
1333+ // // {
1334+ // // for(int inb=0; inb < occ_number.nc; ++inb) std::cout << occ_number(ik, inb) << " ";
1335+ // // std::cout << "\n" << std::endl;
1336+ // // }
1337+
1338+ // // std::cout << "\n\n wg(ik, inband): " << std::endl;
1339+ // // for(int ik=0; ik < occ_number.nr; ++ik)
1340+ // // {
1341+ // // for(int inb=0; inb < occ_number.nc; ++inb)
1342+ // // {
1343+ // // std::cout << this->pelec->wg(ik, inb) << " ";
1344+ // // minus += this->pelec->wg(ik, inb) - occ_number(ik, inb)*this->kv.wk[ik];
1345+ // // }
1346+ // // std::cout << "\n" << std::endl;
1347+ // // }
1348+
1349+ // // if( std::abs(minus) < 1e-12 ) std::cout << "\n\nminus<1e-12, minus: " << minus << std::endl;
1350+ // // else if( std::abs(minus) < 1e-16 ) std::cout << "\n\nminus<1e-16, minus: " << minus << std::endl;
1351+ // // else {std::cout << "\n\nminus is big, minus: " << minus << std::endl;}
1352+
1353+ // ModuleBase::timer::tick("RDMFT", "E & Egradient");
1354+
1355+ // /******** test RDMFT *********/
13551356
13561357
13571358#ifdef __EXX
0 commit comments