Skip to content

Commit 4c44b46

Browse files
JGHan7JGHan7
committed
merge all rdmft-code to abacus-v3.8.0, next step is debug and test
1 parent aa7a4d9 commit 4c44b46

File tree

9 files changed

+92
-26
lines changed

9 files changed

+92
-26
lines changed

‎source/module_esolver/esolver_ks.cpp‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -624,11 +624,10 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
624624
#endif //__RAPIDJSON
625625

626626
// 12.5) rdmft, add by jghan 2024-04-08/2024-10-09
627-
// if ( PARAM.inp.ab_initio_type == "rdmft" )
628-
if ( 1 ) /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
627+
if ( PARAM.inp.ab_initio_type == "rdmft" )
629628
{
630629
// if ( (!GlobalC::exx_info.info_global.cal_exx && iter == 1) || one_step_exx )
631-
if ( 1 )
630+
if ( !GlobalC::exx_info.info_global.cal_exx || (GlobalC::exx_info.info_global.cal_exx && one_step_exx) )
632631
{
633632
ModuleBase::matrix occ_number_ks(this->pelec->wg);
634633
for(int ik=0; ik < occ_number_ks.nr; ++ik)

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 2 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -262,13 +262,10 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(const Input_para& inp, UnitCell
262262

263263

264264
// add by jghan for rdmft calculation
265-
// if( PARAM.inp.ab_initio_type == "rdmft" )
266-
if( 1 )
265+
if( PARAM.inp.ab_initio_type == "rdmft" )
267266
{
268-
// rdmft_solver.init( this->UHM.GG, this->UHM.GK, this->orb_con.ParaV, ucell, this->kv, *(this->pelec),
269-
// GlobalV::DFT_FUNCTIONAL, GlobalV::rdmft_power_alpha);
270267
rdmft_solver.init( this->GG, this->GK, this->pv, ucell, this->kv, *(this->pelec),
271-
this->orb_, two_center_bundle_, PARAM.inp.dft_functional, 1.0);
268+
this->orb_, two_center_bundle_, PARAM.inp.dft_functional, PARAM.inp.rdmft_power_alpha);
272269

273270
// the initialization and necessary calculations of these quantities have been completed in init()
274271
// rdmft_solver.update_ion(ucell, LM, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac, this->LOC);
@@ -921,17 +918,6 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
921918
this->pelec->cal_converged();
922919
}
923920

924-
// if( iter==1 ) // add by jghan, 2024-03-16
925-
// {
926-
// ModuleBase::matrix occ_number(this->pelec->wg);
927-
// for(int ik=0; ik < occ_number.nr; ++ik)
928-
// {
929-
// for(int inb=0; inb < this->pelec->wg.nc; ++inb) occ_number(ik, inb) /= this->kv.wk[ik];
930-
// }
931-
// rdmft_solver.update_elec(occ_number, *(this->psi));
932-
// std::cout << "\n******\nrdmft_solver: " << "update_elec" << "\n******" << std::endl;
933-
// }
934-
935921
}
936922

937923
//------------------------------------------------------------------------------

‎source/module_esolver/lcao_before_scf.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -351,8 +351,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(const int istep)
351351
this->p_hamilt->non_first_scf = istep;
352352

353353
// update in ion-step
354-
// if( PARAM.inp.ab_initio_type == "rdmft" )
355-
if( 1 )
354+
if( PARAM.inp.ab_initio_type == "rdmft" )
356355
{
357356
// necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp
358357
rdmft_solver.update_ion(GlobalC::ucell, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08

‎source/module_io/input_conv.cpp‎

Lines changed: 36 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -330,7 +330,15 @@ void Input_Conv::Convert()
330330
|| dft_functional_lower == "scan0") {
331331
GlobalC::restart.info_save.save_charge = true;
332332
GlobalC::restart.info_save.save_H = true;
333-
} else {
333+
}
334+
else if ( dft_functional_lower == "muller" || dft_functional_lower == "power"
335+
|| dft_functional_lower == "wp22"
336+
|| dft_functional_lower == "cwp22" ) // added by jghan, 2024-07-07
337+
{
338+
GlobalC::restart.info_save.save_charge = true;
339+
GlobalC::restart.info_save.save_H = true;
340+
}
341+
else {
334342
GlobalC::restart.info_save.save_charge = true;
335343
}
336344
}
@@ -348,7 +356,15 @@ void Input_Conv::Convert()
348356
|| dft_functional_lower == "scan0") {
349357
GlobalC::restart.info_load.load_charge = true;
350358
GlobalC::restart.info_load.load_H = true;
351-
} else {
359+
}
360+
else if ( dft_functional_lower == "muller" || dft_functional_lower == "power"
361+
|| dft_functional_lower == "wp22"
362+
|| dft_functional_lower == "cwp22" ) // added by jghan, 2024-07-07
363+
{
364+
GlobalC::restart.info_load.load_charge = true;
365+
GlobalC::restart.info_load.load_H = true;
366+
}
367+
else {
352368
GlobalC::restart.info_load.load_charge = true;
353369
}
354370
}
@@ -376,7 +392,24 @@ void Input_Conv::Convert()
376392
} else if (dft_functional_lower == "opt_orb") {
377393
GlobalC::exx_info.info_global.cal_exx = false;
378394
Exx_Abfs::Jle::generate_matrix = true;
379-
} else {
395+
}
396+
// muller, power, wp22, cwp22 added by jghan, 2024-07-07
397+
else if ( dft_functional_lower == "muller" || dft_functional_lower == "power" )
398+
{
399+
GlobalC::exx_info.info_global.cal_exx = true;
400+
GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf;
401+
}
402+
else if ( dft_functional_lower == "wp22" )
403+
{
404+
GlobalC::exx_info.info_global.cal_exx = true;
405+
GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::erf; // use the error function erf(w|r-r'|), exx just has the long-range part
406+
}
407+
else if ( dft_functional_lower == "cwp22" )
408+
{
409+
GlobalC::exx_info.info_global.cal_exx = true;
410+
GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hse; // use the erfc(w|r-r'|), exx just has the short-range part
411+
}
412+
else {
380413
GlobalC::exx_info.info_global.cal_exx = false;
381414
}
382415

‎source/module_io/read_input_item_exx_dftu.cpp‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,12 @@ void ReadInput::item_exx()
2626
{
2727
para.input.exx_hybrid_alpha = "0.25";
2828
}
29+
// added by jghan 2024-07-06
30+
else if (dft_functional_lower == "muller" || dft_functional_lower == "power"
31+
|| dft_functional_lower == "wp22" || dft_functional_lower == "cwp22")
32+
{
33+
para.input.exx_hybrid_alpha = "1";
34+
}
2935
else
3036
{ // no exx in scf, but will change to non-zero in
3137
// postprocess like rpa
@@ -193,6 +199,20 @@ void ReadInput::item_exx()
193199
{
194200
para.input.exx_ccp_rmesh_times = "1.5";
195201
}
202+
// added by jghan 2024-07-06
203+
else if (dft_functional_lower == "muller" || dft_functional_lower == "power")
204+
{
205+
para.input.exx_ccp_rmesh_times = "5";
206+
}
207+
else if (dft_functional_lower == "wp22")
208+
{
209+
para.input.exx_ccp_rmesh_times = "5";
210+
// exx_ccp_rmesh_times = "1.5";
211+
}
212+
else if (dft_functional_lower == "cwp22")
213+
{
214+
para.input.exx_ccp_rmesh_times = "1.5";
215+
}
196216
else
197217
{ // no exx in scf
198218
para.input.exx_ccp_rmesh_times = "1";

‎source/module_io/read_input_item_other.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -507,7 +507,7 @@ void ReadInput::item_others()
507507
"used in exx-type functionals such as muller and power";
508508
read_sync_double(input.rdmft_power_alpha);
509509
item.reset_value = [](const Input_Item& item, Parameter& para) {
510-
if( para.input.dft_functional == "hf" )
510+
if( para.input.dft_functional == "hf" || para.input.dft_functional == "pbe0" )
511511
{
512512
para.input.rdmft_power_alpha = 1.0;
513513
}

‎source/module_ri/Exx_LRI_interface.hpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,12 @@ void Exx_LRI_Interface<T, Tdata>::exx_beforescf(const K_Vectors& kv, const Charg
5353
{
5454
XC_Functional::set_xc_type("scan");
5555
}
56+
// added by jghan, 2024-07-07
57+
else if ( ucell.atoms[0].ncpp.xc_func == "MULLER" || ucell.atoms[0].ncpp.xc_func == "POWER"
58+
|| ucell.atoms[0].ncpp.xc_func == "WP22" || ucell.atoms[0].ncpp.xc_func == "CWP22" )
59+
{
60+
XC_Functional::set_xc_type("pbe");
61+
}
5662
}
5763
// initialize the rotation matrix in AO representation
5864
this->exx_spacegroup_symmetry = (PARAM.inp.nspin < 4 && ModuleSymmetry::Symmetry::symm_flag == 1);

‎source/module_ri/conv_coulomb_pot_k.cpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,20 @@ namespace Conv_Coulomb_Pot_K
4141
return psik2_ccp;
4242
}
4343

44+
// added by jghan, 2024-07-06
45+
// for using the long-range part of exx = hf - hse(i.e. exx_short-range)
46+
std::vector<double> cal_psi_erf(
47+
const std::vector<double> & psif,
48+
const std::vector<double> & k_radial,
49+
const double hse_omega,
50+
const double hf_Rcut)
51+
{
52+
std::vector<double> psik2_ccp(psif.size());
53+
for( size_t ik=0; ik<psif.size(); ++ik )
54+
psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * ( std::exp(-(k_radial[ik]*k_radial[ik])/(4*hse_omega*hse_omega)) - std::cos(k_radial[ik] * hf_Rcut) );
55+
// psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * ( std::exp(-(k_radial[ik]*k_radial[ik])/(4*hse_omega*hse_omega)) );
56+
return psik2_ccp;
57+
}
4458

4559

4660
template<>
@@ -59,6 +73,8 @@ namespace Conv_Coulomb_Pot_K
5973
psik2_ccp = cal_psi_hf( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hf_Rcut")); break;
6074
case Ccp_Type::Hse:
6175
psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") ); break;
76+
case Ccp_Type::erf:
77+
psik2_ccp = cal_psi_erf( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega"), parameter.at("hf_Rcut") ); break;
6278
default:
6379
throw( ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) ); break;
6480
}

‎source/module_ri/conv_coulomb_pot_k.h‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,8 @@ namespace Conv_Coulomb_Pot_K
1010
enum class Ccp_Type{ // parameter:
1111
Ccp, //
1212
Hf, // "hf_Rcut"
13-
Hse}; // "hse_omega"
13+
Hse, // "hse_omega"
14+
erf}; // "hse_omega"
1415

1516
template<typename T> T cal_orbs_ccp(
1617
const T &orbs,
@@ -34,6 +35,12 @@ namespace Conv_Coulomb_Pot_K
3435
const std::vector<double> & psif,
3536
const std::vector<double> & k_radial,
3637
const double hse_omega);
38+
// added by jghan, 2024-07-06
39+
std::vector<double> cal_psi_erf(
40+
const std::vector<double> & psif,
41+
const std::vector<double> & k_radial,
42+
const double hse_omega,
43+
const double hf_Rcut);
3744
}
3845

3946
#include "conv_coulomb_pot_k.hpp"

0 commit comments

Comments
 (0)