Skip to content

Commit ac00ce4

Browse files
JGHan7JGHan7
committed
just save for updating branch rdmft
1 parent 6167846 commit ac00ce4

File tree

3 files changed

+139
-37
lines changed

3 files changed

+139
-37
lines changed

‎source/module_esolver/esolver_ks.cpp‎

Lines changed: 81 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,9 @@
2323
#include "module_base/parallel_common.h"
2424
#endif
2525

26+
// test by jghan
27+
#include "module_rdmft/rdmft_tools.h"
28+
2629
namespace ModuleESolver
2730
{
2831

@@ -381,33 +384,34 @@ namespace ModuleESolver
381384
auto iterstart = std::chrono::system_clock::now();
382385
#endif
383386

384-
ModuleBase::TITLE("RDMFT", "E & Egradient");
385-
ModuleBase::timer::tick("RDMFT", "E & Egradient");
386-
if( iter==2 ) // ( iter>=2 && GlobalV::CALCULATION == "rdmft" && ModuleESolver::determine_type() == "ksdft_lcao" )
387-
{
388-
if( iter==2 )
389-
{
390-
ModuleBase::matrix occ_number_ks(this->pelec->wg);
391-
for(int ik=0; ik < occ_number_ks.nr; ++ik)
392-
{
393-
for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik];
394-
}
395-
this->update_elec_rdmft(occ_number_ks, *(this->psi));
396-
}
397-
else
398-
{
399-
// this should update occ_number and wfc in another way when iter>2
400-
// this->update_elec_rdmft(occ_number, wfc);
401-
}
402-
403-
// do rdmft calculation
404-
ModuleBase::matrix E_gradient_occNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
405-
psi::Psi<T> E_gradient_wfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis());
406-
double Etotal = this->Run_rdmft(E_gradient_occNum, E_gradient_wfc); // add by jghan 2024-03-16
387+
// ModuleBase::TITLE("RDMFT", "E & Egradient");
388+
// ModuleBase::timer::tick("RDMFT", "E & Egradient");
389+
// if( iter==2 ) // ( iter>=2 && GlobalV::CALCULATION == "rdmft" && ModuleESolver::determine_type() == "ksdft_lcao" )
390+
// {
391+
// std::cout << "\n\n\n******\ndo rdmft_esolver.update_elec() successfully\n******\n\n\n" << std::endl;
392+
// if( iter==2 )
393+
// {
394+
// ModuleBase::matrix occ_number_ks(this->pelec->wg);
395+
// for(int ik=0; ik < occ_number_ks.nr; ++ik)
396+
// {
397+
// for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik];
398+
// }
399+
// this->update_elec_rdmft(occ_number_ks, *(this->psi));
400+
// }
401+
// else
402+
// {
403+
// // this should update occ_number and wfc in another way when iter>2
404+
// // this->update_elec_rdmft(occ_number, wfc);
405+
// }
406+
407+
// // do rdmft calculation
408+
// ModuleBase::matrix E_gradient_occNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
409+
// psi::Psi<T> E_gradient_wfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis());
410+
// double Etotal = this->Run_rdmft(E_gradient_occNum, E_gradient_wfc); // add by jghan 2024-03-16
407411

408-
// continue;
409-
}
410-
ModuleBase::timer::tick("RDMFT", "E & Egradient");
412+
// // continue;
413+
// }
414+
// ModuleBase::timer::tick("RDMFT", "E & Egradient");
411415

412416

413417
double diag_ethr = this->phsol->set_diagethr(istep, iter, drho);
@@ -491,10 +495,61 @@ namespace ModuleESolver
491495
SCF print: G1 -3.435545e+03 0.000000e+00 3.607e-01 2.862e-01
492496
*/
493497
printiter(iter, drho, duration, diag_ethr);
498+
499+
// // test by jghan for rdmft
500+
// if(GlobalV::dm_obj_type == "rdmft" && iter==2 && !GlobalC::exx_info.info_global.cal_exx)
501+
// {
502+
// this->conv_elec = true;
503+
// }
504+
494505
if (this->conv_elec)
495506
{
496507
this->niter = iter;
497508
bool stop = this->do_after_converge(iter);
509+
510+
// bool stop = false;
511+
// if( GlobalV::dm_obj_type != "rdmft" || GlobalC::exx_info.info_global.cal_exx )
512+
// {
513+
// bool stop = this->do_after_converge(iter);
514+
// }
515+
516+
// add by jghan 2024-04-08
517+
// to get the initial values of wfc and occ_numbers
518+
if( GlobalV::dm_obj_type == "rdmft" || 1)
519+
// if( GlobalV::dm_obj_type == "rdmft" || (GlobalV::ESOLVER_TYPE == "dm" && GlobalV::DFT_FUNCTIONAL == "hf") || (GlobalV::ESOLVER_TYPE == "dm" && GlobalV::DFT_FUNCTIONAL == "pbe0") )
520+
{
521+
ModuleBase::matrix occ_number_ks(this->pelec->wg);
522+
for(int ik=0; ik < occ_number_ks.nr; ++ik)
523+
{
524+
for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik];
525+
}
526+
this->update_elec_rdmft(occ_number_ks, *(this->psi));
527+
528+
std::cout << "\n\n\n******\nget the initial values of wfc and occ_numbers successfully\n******\n\n\n" << std::endl;
529+
530+
// // just test
531+
// ModuleBase::matrix E_gradient_occNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
532+
// psi::Psi<T> E_gradient_wfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis());
533+
// double E_RDMFT = this->Run_rdmft(E_gradient_occNum, E_gradient_wfc);
534+
535+
// // delete in the future
536+
// std::cout << "\n\n occ_number(ik, inband): " << std::endl;
537+
// for(int ik=0; ik < occ_number_ks.nr; ++ik)
538+
// {
539+
// for(int inb=0; inb < occ_number_ks.nc; ++inb) std::cout << occ_number_ks(ik, inb) << " ";
540+
// std::cout << "\n" << std::endl;
541+
// }
542+
543+
// std::cout << "\n\n wg(ik, inband): " << std::endl;
544+
// for(int ik=0; ik < occ_number_ks.nr; ++ik)
545+
// {
546+
// for(int inb=0; inb < occ_number_ks.nc; ++inb) std::cout << this->pelec->wg(ik, inb) << " ";
547+
// std::cout << "\n" << std::endl;
548+
// }
549+
550+
break;
551+
}
552+
498553
if(stop) break;
499554
}
500555
}

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 53 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -227,7 +227,8 @@ namespace ModuleESolver
227227
// add by jghan for rdmft calculation
228228
if( GlobalV::CALCULATION == "rdmft" || true )
229229
{
230-
rdmft_solver.init( this->UHM.GG, this->UHM.GK, this->orb_con.ParaV, ucell, this->kv, *(this->pelec->charge) );
230+
rdmft_solver.init( this->UHM.GG, this->UHM.GK, this->orb_con.ParaV, ucell, this->kv, *(this->pelec),
231+
GlobalV::DFT_FUNCTIONAL, GlobalV::rdmft_power_alpha);
231232

232233
// the initialization and necessary calculations of these quantities have been completed in init()
233234
// rdmft_solver.update_ion(ucell, LM, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac, this->LOC);
@@ -936,7 +937,7 @@ namespace ModuleESolver
936937
// GlobalC::ppcell.vloc(ModuleBase::matrix vloc),
937938
// esolver_fp.h(Structure_Factor sf), structure_factor.h(ModuleBase::ComplexMatrix strucFac), this->sf.strucFac
938939
// elecstate::DensityMatrix<TK, double>& *( dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM() )
939-
940+
// elecstate.h (Potential* pot, this->pelec->pot)
940941

941942
//initialize the gradients of Etotal on wg and wfc, and set all elements to 0.
942943
ModuleBase::matrix E_gradient_occNum(this->pelec->wg.nr, this->pelec->wg.nc, true);
@@ -975,6 +976,7 @@ namespace ModuleESolver
975976
&LM,
976977
LM.ParaV,
977978
occ_number,
979+
this->pelec->wg,
978980
wfc_rdmft,
979981
E_gradient_occNum,
980982
E_gradient_wfc,
@@ -985,8 +987,8 @@ namespace ModuleESolver
985987
*(this->pw_rho),
986988
GlobalC::ppcell.vloc,
987989
this->sf.strucFac,
988-
"power",
989-
0.5
990+
"HF",
991+
1.0
990992
);
991993
}
992994
// multi-k calculation
@@ -996,6 +998,7 @@ namespace ModuleESolver
996998
&LM,
997999
LM.ParaV,
9981000
occ_number,
1001+
this->pelec->wg,
9991002
*(this->psi),
10001003
E_gradient_occNum,
10011004
E_gradient_wfc,
@@ -1006,20 +1009,61 @@ namespace ModuleESolver
10061009
*(this->pw_rho),
10071010
GlobalC::ppcell.vloc,
10081011
this->sf.strucFac,
1009-
"power",
1010-
0.5
1012+
"HF",
1013+
1.0
10111014
);
10121015
}
10131016

1014-
// ModuleBase::TITLE("RDMFT", "E & Egradient");
1015-
// ModuleBase::timer::tick("RDMFT", "E & Egradient");
1017+
// GlobalV::ofs_running << "E_one_elec "
1018+
// << std::setw(10) << pelec->f_en.deband + pelec->f_en.eband << std::endl;
1019+
// GlobalV::ofs_running << "E_Hartree "
1020+
// << std::setw(10) << pelec->f_en.hartree_energy << std::endl;
1021+
// GlobalV::ofs_running << "E_xc "
1022+
// << std::setw(10) << pelec->f_en.etxc - pelec->f_en.etxcc<< std::endl;
1023+
1024+
double E_xc_KS = this->pelec->f_en.etxc - this->pelec->f_en.etxcc;
1025+
1026+
double E_TV_Hartree_EXX = this->pelec->f_en.deband + this->pelec->f_en.eband
1027+
+ this->pelec->f_en.hartree_energy + this->pelec->f_en.exx;
1028+
double E_TV_Hartree_XC = this->pelec->f_en.deband + this->pelec->f_en.eband
1029+
+ this->pelec->f_en.hartree_energy + this->pelec->f_en.etxc - this->pelec->f_en.etxcc;
1030+
std::cout << std::fixed << std::setprecision(10) << "\n\n******\nE(one_elec + Hartree + EXX) by KS-DFT: " << E_TV_Hartree_EXX
1031+
<< "\n\nE(one_elec + Hartree + XC) by KS-DFT: " << E_TV_Hartree_XC << "\nExc_KS: " << E_xc_KS << "\n******\n\n" << std::endl;
1032+
1033+
ModuleBase::TITLE("RDMFT", "E & Egradient");
1034+
ModuleBase::timer::tick("RDMFT", "E & Egradient");
10161035

10171036
// // test class rdmft
10181037
// rdmft_solver.update_elec(occ_number, *(this->psi));
10191038
// this->Run_rdmft(E_gradient_occNum, E_gradient_wfc);
10201039
// std::cout << "\nrdmft_solver: " << "0.0000" << std::endl;
10211040

1022-
// ModuleBase::timer::tick("RDMFT", "E & Egradient");
1041+
// double minus = 0.0;
1042+
1043+
// // delete in the future
1044+
// std::cout << "\n\n occ_number(ik, inband): " << std::endl;
1045+
// for(int ik=0; ik < occ_number.nr; ++ik)
1046+
// {
1047+
// for(int inb=0; inb < occ_number.nc; ++inb) std::cout << occ_number(ik, inb) << " ";
1048+
// std::cout << "\n" << std::endl;
1049+
// }
1050+
1051+
// std::cout << "\n\n wg(ik, inband): " << std::endl;
1052+
// for(int ik=0; ik < occ_number.nr; ++ik)
1053+
// {
1054+
// for(int inb=0; inb < occ_number.nc; ++inb)
1055+
// {
1056+
// std::cout << this->pelec->wg(ik, inb) << " ";
1057+
// minus += this->pelec->wg(ik, inb) - occ_number(ik, inb)*this->kv.wk[ik];
1058+
// }
1059+
// std::cout << "\n" << std::endl;
1060+
// }
1061+
1062+
// if( std::abs(minus) < 1e-12 ) std::cout << "\n\nminus<1e-12, minus: " << minus << std::endl;
1063+
// else if( std::abs(minus) < 1e-16 ) std::cout << "\n\nminus<1e-16, minus: " << minus << std::endl;
1064+
// else {std::cout << "\n\nminus is big, minus: " << minus << std::endl;}
1065+
1066+
ModuleBase::timer::tick("RDMFT", "E & Egradient");
10231067

10241068
/******** test RDMFT *********/
10251069

‎source/module_esolver/esolver_ks_lcao_elec.cpp‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -324,8 +324,11 @@ void ESolver_KS_LCAO<TK, TR>::beforescf(int istep)
324324
this->p_hamilt->non_first_scf = istep;
325325

326326
// update in ion-step
327-
// necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp
328-
rdmft_solver.update_ion(GlobalC::ucell, LM, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac, this->LOC); // add by jghan, 2024-03-16
327+
if( GlobalV::dm_obj_type == "rdmft" )
328+
{
329+
// necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp
330+
rdmft_solver.update_ion(GlobalC::ucell, LM, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac, this->LOC); // add by jghan, 2024-03-16
331+
}
329332

330333
ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf");
331334
return;

0 commit comments

Comments
 (0)