diff --git a/CMakeLists.txt b/CMakeLists.txt index 85069429b1..168f439597 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -724,7 +724,8 @@ if(ENABLE_LCAO) hcontainer deltaspin numerical_atomic_orbitals - lr) + lr + rdmft) if(USE_ELPA) target_link_libraries(${ABACUS_BIN_NAME} genelpa) endif() diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 40edbf5538..f16616c1c6 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -434,6 +434,9 @@ - [abs\_broadening](#abs_broadening) - [ri\_hartree\_benchmark](#ri_hartree_benchmark) - [aims\_nbasis](#aims_nbasis) + - [Reduced Density Matrix Functional Theory](#Reduced-Density-Matrix-Functional-Theory) + - [rdmft](#rdmft) + - [rdmft\_power\_alpha](#rdmft_power_alpha) [back to top](#full-list-of-input-keywords) ## System variables @@ -4021,4 +4024,21 @@ The output files are `OUT.${suffix}/Excitation_Energy.dat` and `OUT.${suffix}/Ex - **Description**: Atomic basis set size for each atom type (with the same order as in `STRU`) in FHI-aims. - **Default**: {} (empty list, where ABACUS use its own basis set size) +## Reduced Density Matrix Functional Theory + +ab-initio methods and the xc-functional parameters used in RDMFT. +The physical quantities that RDMFT temporarily expects to output are the kinetic energy, total energy, and 1-RDM of the system in the ground state, etc. + +### rdmft + +- **Type**: Boolean +- **Description**: Whether to perform rdmft calculation (reduced density matrix funcional theory) +- **Default**: false + +### rdmft_power_alpha + +- **Type**: Real +- **Description**: The alpha parameter of power-functional(or other exx-type/hybrid functionals) which used in RDMFT, g(occ_number) = occ_number^alpha +- **Default**: 0.656 + [back to top](#full-list-of-input-keywords) diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index d57fd16337..e8af89216b 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -18,6 +18,9 @@ add_subdirectory(module_ri) add_subdirectory(module_parameter) add_subdirectory(module_lr) +# add by jghan +add_subdirectory(module_rdmft) + add_library( driver OBJECT diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 504349de48..609c395e7c 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -74,6 +74,7 @@ VPATH=./src_global:\ ./module_lr/operator_casida:\ ./module_lr/potentials:\ ./module_lr/utils:\ +./module_rdmft:\ ./\ OBJS_ABACUS_PW=${OBJS_MAIN}\ @@ -115,6 +116,7 @@ ${OBJS_DELTASPIN}\ ${OBJS_TENSOR}\ ${OBJS_HSOLVER_PEXSI}\ ${OBJS_LR}\ +${OBJS_RDMFT} OBJS_MAIN=main.o\ driver.o\ @@ -735,3 +737,8 @@ OBJS_TENSOR=tensor.o\ lr_spectrum.o\ hamilt_casida.o\ esolver_lrtd_lcao.o\ + +OBJS_RDMFT=rdmft.o\ + rdmft_tools.o\ + cal_V_rdmft.o\ + update_state_rdmft.o\ diff --git a/source/module_esolver/esolver_ks.h b/source/module_esolver/esolver_ks.h index f18264f7f7..d281946b59 100644 --- a/source/module_esolver/esolver_ks.h +++ b/source/module_esolver/esolver_ks.h @@ -70,6 +70,9 @@ class ESolver_KS : public ESolver_FP // wavefunction coefficients psi::Psi* psi = nullptr; + // // TR need be std::complex when NSPIN = 4 + // rdmft::RDMFT rdmft_solver; // add by jghan for rdmft calculation + std::string basisname; // PW or LCAO double esolver_KS_ne = 0.0; bool oscillate_esolver = false; // whether esolver is oscillated diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 069ec0e8c7..97171f31a0 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -60,6 +60,11 @@ // #include "module_elecstate/cal_dm.h" //--------------------------------------------------- +// test RDMFT +#include "module_rdmft/rdmft.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" //temp ,delete +#include + namespace ModuleESolver { @@ -261,7 +266,15 @@ void ESolver_KS_LCAO::before_all_runners(const Input_para& inp, UnitCell this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.inp.nbands, PARAM.inp.nelec); } - // 15) if kpar is not divisible by nks, print a warning + + // 15) initialize rdmft, added by jghan + if( PARAM.inp.rdmft == true ) + { + rdmft_solver.init( this->GG, this->GK, this->pv, ucell, this->kv, *(this->pelec), + this->orb_, two_center_bundle_, PARAM.inp.dft_functional, PARAM.inp.rdmft_power_alpha); + } + + // 16) if kpar is not divisible by nks, print a warning if (GlobalV::KPAR_LCAO > 1) { if (this->kv.get_nks() % GlobalV::KPAR_LCAO != 0) @@ -850,6 +863,7 @@ void ESolver_KS_LCAO::update_pot(const int istep, const int iter) { this->pelec->cal_converged(); } + } //------------------------------------------------------------------------------ @@ -924,7 +938,24 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) } } + // 2.5) determine whether rdmft needs to get the initial value, added by jghan, 2024-10-25 + int one_step_exx = false; + bool get_init_value_rdmft = false; + if( iter == 1 ) get_init_value_rdmft = true; // the case without hybrid functionals #ifdef __EXX + if( GlobalC::exx_info.info_global.cal_exx ) + { + if( this->conv_esolver ) one_step_exx = true; + // the case with hybrid functionals, calculate rdmft just after updateExx, cal once in one inner loop + if( one_step_exx ) get_init_value_rdmft = true; + else get_init_value_rdmft = false; + } +#endif + +#ifdef __EXX + + if( GlobalC::exx_info.info_global.cal_exx && this->conv_esolver ) one_step_exx = true; + // 3) save exx matrix if (GlobalC::exx_info.info_global.cal_exx) { @@ -988,6 +1019,23 @@ void ESolver_KS_LCAO::iter_finish(const int istep, int& iter) { GlobalC::dftu.initialed_locale = true; } + + // 7) rdmft, added by jghan, 2024-10-25 + if ( PARAM.inp.rdmft == true && get_init_value_rdmft ) + { + ModuleBase::matrix occ_number_ks(this->pelec->wg); + for(int ik=0; ik < occ_number_ks.nr; ++ik) { for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik]; } + this->rdmft_solver.update_elec(occ_number_ks, *(this->psi)); + + //initialize the gradients of Etotal on occupation numbers and wfc, and set all elements to 0. + ModuleBase::matrix dE_dOccNum(this->pelec->wg.nr, this->pelec->wg.nc, true); + psi::Psi dE_dWfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis()); + dE_dWfc.zero_out(); + + double Etotal_RDMFT = this->rdmft_solver.run(dE_dOccNum, dE_dWfc); + // break; + } + } //------------------------------------------------------------------------------ @@ -1082,6 +1130,23 @@ void ESolver_KS_LCAO::after_scf(const int istep) ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels"); #endif + /******** test RDMFT *********/ + if ( PARAM.inp.rdmft == true ) // rdmft, added by jghan, 2024-10-17 + { + ModuleBase::matrix occ_number_ks(this->pelec->wg); + for(int ik=0; ik < occ_number_ks.nr; ++ik) { for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik]; } + this->rdmft_solver.update_elec(occ_number_ks, *(this->psi)); + + //initialize the gradients of Etotal on occupation numbers and wfc, and set all elements to 0. + ModuleBase::matrix dE_dOccNum(this->pelec->wg.nr, this->pelec->wg.nc, true); + psi::Psi dE_dWfc(this->psi->get_nk(), this->psi->get_nbands(), this->psi->get_nbasis()); + dE_dWfc.zero_out(); + + double Etotal_RDMFT = this->rdmft_solver.run(dE_dOccNum, dE_dWfc); + } + /******** test RDMFT *********/ + + #ifdef __EXX // 11) write rpa information if (PARAM.inp.rpa) @@ -1254,6 +1319,7 @@ void ESolver_KS_LCAO::after_scf(const int istep) } } + //------------------------------------------------------------------------------ //! the 20th,21th,22th functions of ESolver_KS_LCAO //! mohan add 2024-05-11 diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 761d660fe5..32c8dbe5aa 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -12,6 +12,9 @@ #include "module_basis/module_nao/two_center_bundle.h" #include "module_io/output_mat_sparse.h" +// added by jghan for rdmft calculation +#include "module_rdmft/rdmft.h" + #include namespace LR @@ -74,6 +77,8 @@ class ESolver_KS_LCAO : public ESolver_KS { TwoCenterBundle two_center_bundle_; + rdmft::RDMFT rdmft_solver; // added by jghan for rdmft calculation, 2024-03-16 + // temporary introduced during removing GlobalC::ORB LCAO_Orbitals orb_; diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 0b12ca5265..ec367b78cc 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -325,6 +325,14 @@ void ESolver_KS_LCAO::before_scf(const int istep) } this->p_hamilt->non_first_scf = istep; + + // update in ion-step + if( PARAM.inp.rdmft == true ) + { + // necessary operation of these parameters have be done with p_esolver->Init() in source/driver_run.cpp + rdmft_solver.update_ion(GlobalC::ucell, *(this->pw_rho), GlobalC::ppcell.vloc, this->sf.strucFac); // add by jghan, 2024-03-16/2024-10-08 + } + return; } diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp index 713ce512fc..38abdfa5d3 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional.cpp @@ -18,6 +18,7 @@ std::vector XC_Functional::func_id(1); int XC_Functional::func_type = 0; bool XC_Functional::use_libxc = true; double XC_Functional::hybrid_alpha = 0.25; +std::map XC_Functional::scaling_factor_xc = { {1, 1.0} }; // added by jghan, 2024-10-10 void XC_Functional::set_hybrid_alpha(const double alpha_in) { @@ -61,6 +62,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) // func_id.push_back(XC_GGA_C_PBE); func_id.clear(); + scaling_factor_xc.clear(); // added by jghan, 2024-07-07 std::string xc_func = xc_func_in; std::transform(xc_func.begin(), xc_func.end(), xc_func.begin(), (::toupper)); if( xc_func == "LDA" || xc_func == "PZ" || xc_func == "SLAPZNOGXNOGC") //SLA+PZ @@ -196,6 +198,11 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) { // not doing anything } + else if( xc_func == "MULLER" || xc_func == "POWER" ) // added by jghan, 2024-07-06 + { + func_type = 4; + use_libxc = false; + } #ifdef USE_LIBXC else if( xc_func == "HSE") { @@ -203,6 +210,56 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_type = 4; use_libxc = true; } + // added by jghan, 2024-07-06 + else if( xc_func == "WP22") + { + func_id.push_back(XC_GGA_X_ITYH); // short-range of B88_X, id=529 + func_id.push_back(XC_GGA_C_LYPR); // short-range of LYP_C, id=624 + func_type = 4; + use_libxc = true; + } + else if( xc_func == "CWP22") + { + // BLYP_XC_lr = -BLYP_XC_sr + BLYP_XC, the realization of it is in v_xc_libxc() function, xc_functional_libxc_vxc.cpp + func_id.push_back(XC_GGA_X_ITYH); // short-range of B88_X, id=529 + func_id.push_back(XC_GGA_C_LYPR); // short-range of LYP_C, id=624 + func_id.push_back(XC_GGA_X_B88); // complete B88_X, id=106 + func_id.push_back(XC_GGA_C_LYP); // complete LYP_C, id=131 + + // the scaling factor of CWP22-functionals + scaling_factor_xc[XC_GGA_X_ITYH] = -1.0; + scaling_factor_xc[XC_GGA_C_LYPR] = -1.0; + scaling_factor_xc[XC_GGA_X_B88] = 1.0; + scaling_factor_xc[XC_GGA_X_B88] = 1.0; + + func_type = 4; + use_libxc = true; + } + // just for test BLYP in libxc, added by jghan, 2024-10-22 + else if( xc_func == "BLYP_LIBXC") + { + func_id.push_back(XC_GGA_X_B88); // complete B88_X, id=106 + func_id.push_back(XC_GGA_C_LYP); // complete LYP_C, id=131 + func_type = 2; + use_libxc = true; + } + else if( xc_func == "BLYP_LR") + { + // BLYP_XC_lr = -BLYP_XC_sr + BLYP_XC, the realization of it is in v_xc_libxc() function, xc_functional_libxc_vxc.cpp + func_id.push_back(XC_GGA_X_ITYH); // short-range of B88_X, id=529 + func_id.push_back(XC_GGA_C_LYPR); // short-range of LYP_C, id=624 + func_id.push_back(XC_GGA_X_B88); // complete B88_X, id=106 + func_id.push_back(XC_GGA_C_LYP); // complete LYP_C, id=131 + + // the scaling factor of BLYP_LR-functionals + scaling_factor_xc[XC_GGA_X_ITYH] = -1.0; + scaling_factor_xc[XC_GGA_C_LYPR] = -1.0; + scaling_factor_xc[XC_GGA_X_B88] = 1.0; + scaling_factor_xc[XC_GGA_X_B88] = 1.0; + + func_type = 2; + use_libxc = true; + } #endif else { @@ -243,7 +300,8 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) #endif #ifndef USE_LIBXC - if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0") + if(xc_func == "SCAN" || xc_func == "HSE" || xc_func == "SCAN0" + || xc_func == "MULLER" || xc_func == "POWER" || xc_func == "WP22" || xc_func == "CWP22") { ModuleBase::WARNING_QUIT("set_xc_type","to use SCAN, SCAN0, or HSE, LIBXC is required"); } diff --git a/source/module_hamilt_general/module_xc/xc_functional.h b/source/module_hamilt_general/module_xc/xc_functional.h index ac9c75256a..04ca90216f 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.h +++ b/source/module_hamilt_general/module_xc/xc_functional.h @@ -20,6 +20,8 @@ #include "module_elecstate/module_charge/charge.h" #include "module_cell/unitcell.h" +#include // added by jghan, 2024-10-10 + class XC_Functional { public: @@ -80,6 +82,10 @@ class XC_Functional //exx_hybrid_alpha for mixing exx in hybrid functional: static double hybrid_alpha; + // added by jghan, 2024-07-07 + // as a scaling factor for different xc-functionals + static std::map scaling_factor_xc; + public: static std::vector get_func_id() { return func_id; } @@ -162,7 +168,7 @@ class XC_Functional ModulePW::PW_Basis* rhopw, const UnitCell* ucell, std::vector& stress_gga, - const bool is_stress = 0); + const bool is_stress = false); template ::type> static void grad_wfc( diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp index c8e4b1b609..06d22e8f2d 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp @@ -141,6 +141,20 @@ std::vector XC_Functional_Libxc::init_func(const std::vector GlobalC::exx_info.info_global.hse_omega }; xc_func_set_ext_params(&funcs.back(), parameter_hse); } + // added by jghan, 2024-07-06 + else if( id == XC_GGA_X_ITYH ) // short-range of B88_X + { + add_func( XC_GGA_X_ITYH ); + double parameter_omega[1] = {PARAM.inp.exx_hse_omega}; // GlobalC::exx_info.info_global.hse_omega + xc_func_set_ext_params(&funcs.back(), parameter_omega); + } + else if( id == XC_GGA_C_LYPR ) // short-range of LYP_C + { + add_func( XC_GGA_C_LYPR ); + // the first six parameters come from libxc, and may need to be modified in some cases + double parameter_lypr[7] = {0.04918, 0.132, 0.2533, 0.349, 0.35/2.29, 2.0/2.29, PARAM.inp.exx_hse_omega}; + xc_func_set_ext_params(&funcs.back(), parameter_lypr); + } #endif else { diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc.h b/source/module_hamilt_general/module_xc/xc_functional_libxc.h index 26044c1b90..f206261076 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc.h +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.h @@ -11,6 +11,9 @@ #include #include +#include // added by jghan, 2024-10-10 +#include + class Charge; namespace XC_Functional_Libxc @@ -33,12 +36,13 @@ namespace XC_Functional_Libxc // xc_functional_libxc_vxc.cpp //------------------- - extern std::tuple v_xc_libxc( - const std::vector &func_id, - const int &nrxx, // number of real-space grid - const double &omega, // volume of cell - const double tpiba, - const Charge* const chr); // charge density + extern std::tuple v_xc_libxc( + const std::vector &func_id, + const int &nrxx, // number of real-space grid + const double &omega, // volume of cell + const double tpiba, + const Charge* const chr, // charge density + const std::map* scaling_factor = nullptr); // added by jghan, 2024-10-10 // for mGGA functional extern std::tuple v_xc_meta( diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp index 8b270cd233..e888c465b8 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp @@ -18,7 +18,8 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / const int &nrxx, // number of real-space grid const double &omega, // volume of cell const double tpiba, - const Charge* const chr) + const Charge* const chr, + const std::map* scaling_factor) { ModuleBase::TITLE("XC_Functional_Libxc","v_xc_libxc"); ModuleBase::timer::tick("XC_Functional_Libxc","v_xc_libxc"); @@ -109,14 +110,24 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / break; } - etxc += XC_Functional_Libxc::convert_etxc(nspin, nrxx, sgn, rho, exc); + // added by jghan, 2024-10-10 + double factor = 1.0; + if( scaling_factor == nullptr ) ; + else + { + auto pair_factor = scaling_factor->find(func.info->number); + if( pair_factor != scaling_factor->end() ) factor = pair_factor->second; + } + + // time factor is added by jghan, 2024-10-10 + etxc += XC_Functional_Libxc::convert_etxc(nspin, nrxx, sgn, rho, exc) * factor; const std::pair vtxc_v = XC_Functional_Libxc::convert_vtxc_v( func, nspin, nrxx, sgn, rho, gdr, vrho, vsigma, tpiba, chr); - vtxc += std::get<0>(vtxc_v); - v += std::get<1>(vtxc_v); + vtxc += std::get<0>(vtxc_v) * factor; + v += std::get<1>(vtxc_v) * factor; } // end for( xc_func_type &func : funcs ) if(4==PARAM.inp.nspin) diff --git a/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp index 7aeb677a78..2c67941adb 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp @@ -25,7 +25,7 @@ std::tuple XC_Functional::v_xc( if(use_libxc) { #ifdef USE_LIBXC - return XC_Functional_Libxc::v_xc_libxc(XC_Functional::get_func_id(), nrxx, ucell->omega, ucell->tpiba, chr); + return XC_Functional_Libxc::v_xc_libxc(XC_Functional::get_func_id(), nrxx, ucell->omega, ucell->tpiba, chr, &(scaling_factor_xc)); #else ModuleBase::WARNING_QUIT("v_xc","compile with LIBXC"); #endif diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 7c3997bc22..f6504df0c8 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -319,7 +319,15 @@ void Input_Conv::Convert() || dft_functional_lower == "scan0") { GlobalC::restart.info_save.save_charge = true; GlobalC::restart.info_save.save_H = true; - } else { + } + else if ( dft_functional_lower == "muller" || dft_functional_lower == "power" + || dft_functional_lower == "wp22" + || dft_functional_lower == "cwp22" ) // added by jghan, 2024-07-07 + { + GlobalC::restart.info_save.save_charge = true; + GlobalC::restart.info_save.save_H = true; + } + else { GlobalC::restart.info_save.save_charge = true; } } @@ -337,7 +345,15 @@ void Input_Conv::Convert() || dft_functional_lower == "scan0") { GlobalC::restart.info_load.load_charge = true; GlobalC::restart.info_load.load_H = true; - } else { + } + else if ( dft_functional_lower == "muller" || dft_functional_lower == "power" + || dft_functional_lower == "wp22" + || dft_functional_lower == "cwp22" ) // added by jghan, 2024-07-07 + { + GlobalC::restart.info_load.load_charge = true; + GlobalC::restart.info_load.load_H = true; + } + else { GlobalC::restart.info_load.load_charge = true; } } @@ -365,7 +381,24 @@ void Input_Conv::Convert() } else if (dft_functional_lower == "opt_orb") { GlobalC::exx_info.info_global.cal_exx = false; Exx_Abfs::Jle::generate_matrix = true; - } else { + } + // muller, power, wp22, cwp22 added by jghan, 2024-07-07 + else if ( dft_functional_lower == "muller" || dft_functional_lower == "power" ) + { + GlobalC::exx_info.info_global.cal_exx = true; + GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf; + } + else if ( dft_functional_lower == "wp22" ) + { + GlobalC::exx_info.info_global.cal_exx = true; + 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 + } + else if ( dft_functional_lower == "cwp22" ) + { + GlobalC::exx_info.info_global.cal_exx = true; + 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 + } + else { GlobalC::exx_info.info_global.cal_exx = false; } diff --git a/source/module_io/read_input_item_exx_dftu.cpp b/source/module_io/read_input_item_exx_dftu.cpp index 87a8edc480..5da216f856 100644 --- a/source/module_io/read_input_item_exx_dftu.cpp +++ b/source/module_io/read_input_item_exx_dftu.cpp @@ -26,6 +26,12 @@ void ReadInput::item_exx() { para.input.exx_hybrid_alpha = "0.25"; } + // added by jghan 2024-07-06 + else if (dft_functional_lower == "muller" || dft_functional_lower == "power" + || dft_functional_lower == "wp22" || dft_functional_lower == "cwp22") + { + para.input.exx_hybrid_alpha = "1"; + } else { // no exx in scf, but will change to non-zero in // postprocess like rpa @@ -193,6 +199,20 @@ void ReadInput::item_exx() { para.input.exx_ccp_rmesh_times = "1.5"; } + // added by jghan 2024-07-06 + else if (dft_functional_lower == "muller" || dft_functional_lower == "power") + { + para.input.exx_ccp_rmesh_times = "5"; + } + else if (dft_functional_lower == "wp22") + { + para.input.exx_ccp_rmesh_times = "5"; + // exx_ccp_rmesh_times = "1.5"; + } + else if (dft_functional_lower == "cwp22") + { + para.input.exx_ccp_rmesh_times = "1.5"; + } else { // no exx in scf para.input.exx_ccp_rmesh_times = "1"; diff --git a/source/module_io/read_input_item_other.cpp b/source/module_io/read_input_item_other.cpp index 875b816918..9cde2fc27b 100644 --- a/source/module_io/read_input_item_other.cpp +++ b/source/module_io/read_input_item_other.cpp @@ -493,5 +493,37 @@ void ReadInput::item_others() sync_intvec(input.aims_nbasis, para.input.aims_nbasis.size(), 0); this->add_item(item); } + + // RDMFT, added by jghan, 2024-10-16 + { + Input_Item item("rdmft"); + item.annotation = "whether to perform rdmft calculation, default is false"; + read_sync_bool(input.rdmft); + this->add_item(item); + } + { + Input_Item item("rdmft_power_alpha"); + item.annotation = "the alpha parameter of power-functional, g(occ_number) = occ_number^alpha" + " used in exx-type functionals such as muller and power"; + read_sync_double(input.rdmft_power_alpha); + item.reset_value = [](const Input_Item& item, Parameter& para) { + if( para.input.dft_functional == "hf" || para.input.dft_functional == "pbe0" ) + { + para.input.rdmft_power_alpha = 1.0; + } + else if( para.input.dft_functional == "muller" ) + { + para.input.rdmft_power_alpha = 0.5; + } + }; + item.check_value = [](const Input_Item& item, const Parameter& para) { + if( (para.input.rdmft_power_alpha < 0) || (para.input.rdmft_power_alpha > 1) ) + { + ModuleBase::WARNING_QUIT("ReadInput", "rdmft_power_alpha should be greater than 0.0 and less than 1.0"); + } + }; + this->add_item(item); + } + } } // namespace ModuleIO \ No newline at end of file diff --git a/source/module_io/test/read_input_ptest.cpp b/source/module_io/test/read_input_ptest.cpp index 69d6d13a1c..a189560898 100644 --- a/source/module_io/test/read_input_ptest.cpp +++ b/source/module_io/test/read_input_ptest.cpp @@ -427,6 +427,8 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.abs_wavelen_range.size(), 2); EXPECT_DOUBLE_EQ(param.inp.abs_wavelen_range[0], 0.0); EXPECT_DOUBLE_EQ(param.inp.abs_broadening, 0.01); + EXPECT_EQ(param.inp.rdmft, 0); + EXPECT_DOUBLE_EQ(param.inp.rdmft_power_alpha, 0.656); } TEST_F(InputParaTest, Check) diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index bfd9ae91cb..ef8f051c90 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -589,5 +589,12 @@ struct Input_para int test_pp = 0; ///< variables for test_pp only int test_relax_method = false; ///< variables for test_relax_method only int test_deconstructor = false; ///< variables for test_deconstructor only + + // ============== #Parameters (21.RDMFT) ===================== + // RDMFT jghan added on 2024-07-06 + bool rdmft = false; // rdmft, reduced density matrix funcional theory + double rdmft_power_alpha = 0.656; // the alpha parameter of power-functional, g(occ_number) = occ_number^alpha + // double rdmft_wp22_omega; // the omega parameter of wp22-functional = exx_hse_omega + }; #endif diff --git a/source/module_rdmft/CMakeLists.txt b/source/module_rdmft/CMakeLists.txt new file mode 100644 index 0000000000..8246745d41 --- /dev/null +++ b/source/module_rdmft/CMakeLists.txt @@ -0,0 +1,20 @@ +if(ENABLE_LCAO) + add_library( + rdmft + OBJECT + rdmft.cpp + rdmft_tools.cpp + cal_V_rdmft.cpp + update_state_rdmft.cpp + ) +endif() + +# if(ENABLE_COVERAGE) +# add_coverage(psi) +# add_coverage(psi_initializer) +# endif() + +# if (BUILD_TESTING) +# add_subdirectory(kernels/test) +# add_subdirectory(test) +# endif() diff --git a/source/module_rdmft/cal_V_rdmft.cpp b/source/module_rdmft/cal_V_rdmft.cpp new file mode 100644 index 0000000000..3d4c27bb19 --- /dev/null +++ b/source/module_rdmft/cal_V_rdmft.cpp @@ -0,0 +1,339 @@ +//========================================================== +// Author: Jingang Han +// DATE : 2024-10-30 +//========================================================== + +#include "rdmft.h" +#include "module_rdmft/rdmft_tools.h" + +#include "module_psi/psi.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" + +namespace rdmft +{ + + + +template +void RDMFT::get_DM_XC(std::vector< std::vector >& DM_XC) +{ + // get wk_funEta_wfc = wk*g(eta)*conj(wfc) + psi::Psi wk_funEta_wfc(wfc); + conj_psi(wk_funEta_wfc); + occNum_MulPsi(ParaV, wk_fun_occNum, wk_funEta_wfc, 0); + + // get the special DM_XC used in constructing V_exx_XC + for(int ik=0; ikdesc_wfc, ParaV->desc); +#else + elecstate::psiMulPsi(wk_funEta_wfc, wfc, DM_Kpointer); +#endif + } +} + + +template +void RDMFT::cal_V_TV() +{ + HR_TV->set_zero(); + + V_ekinetic_potential = new hamilt::EkineticNew>( + hsk_TV, + kv->kvec_d, + HR_TV, + &GlobalC::ucell, + orb->cutoffs(), + &GlobalC::GridD, + two_center_bundle->kinetic_orb.get() + ); + + V_nonlocal = new hamilt::NonlocalNew>( + hsk_TV, + kv->kvec_d, + HR_TV, + &GlobalC::ucell, + orb->cutoffs(), + &GlobalC::GridD, + two_center_bundle->overlap_orb_beta.get() + ); + + if( PARAM.inp.gamma_only ) + { + V_local = new rdmft::Veff_rdmft( + GG, + hsk_TV, + kv->kvec_d, + this->pelec->pot, + HR_TV, + &GlobalC::ucell, + orb->cutoffs(), + &GlobalC::GridD, + nspin, + + charge, + rho_basis, + vloc, + sf, + "local" + ); + } + else + { + V_local = new rdmft::Veff_rdmft( + GK, + hsk_TV, + kv->kvec_d, + this->pelec->pot, + HR_TV, + &GlobalC::ucell, + orb->cutoffs(), + &GlobalC::GridD, + nspin, + + charge, + rho_basis, + vloc, + sf, + "local" + ); + } + + // update HR_TV in ion-step, now HR_TV has the HR of V_ekinetic + V_nonlcao + V_local + V_ekinetic_potential->contributeHR(); + V_nonlocal->contributeHR(); + V_local->contributeHR(); + +} + + +template +void RDMFT::cal_V_hartree() +{ + HR_hartree->set_zero(); + + if( PARAM.inp.gamma_only ) + { + V_hartree = new rdmft::Veff_rdmft( + GG, + hsk_hartree, + kv->kvec_d, + this->pelec->pot, + HR_hartree, + &GlobalC::ucell, + orb->cutoffs(), + &GlobalC::GridD, + nspin, + + charge, + rho_basis, + vloc, + sf, + "hartree" + ); + } + else + { + // this can be optimized, use potHartree.update_from_charge() + V_hartree = new rdmft::Veff_rdmft( + GK, + hsk_hartree, + kv->kvec_d, + this->pelec->pot, + HR_hartree, + &GlobalC::ucell, + orb->cutoffs(), + &GlobalC::GridD, + nspin, + + charge, + rho_basis, + vloc, + sf, + "hartree" + ); + } + + // in gamma only, must calculate HR_hartree before HR_local + // HR_exx_XC get from another way, so don't need to do this + V_hartree->contributeHR(); + + // // update HR_local in e-step, now HR_TV has the HR of V_ekinetic + V_nonlcao + V_local, + // V_local->contributeHR(); + // HR_local->add(*HR_TV); // now HR_local has the HR of V_ekinetic + V_nonlcao + V_local + +} + + +template +void RDMFT::cal_V_XC() +{ + // // //test + // DM_XC_pass = DM_XC; + + // elecstate::DensityMatrix DM_test(ParaV, nspin, kv->kvec_d, nk_total); + // elecstate::cal_dm_psi(ParaV, wg, wfc, DM_test); + // DM_test.init_DMR(&GlobalC::GridD, &GlobalC::ucell); + // DM_test.cal_DMR(); + + // // compare DM_XC and DM get in update_charge(or ABACUS) + // std::cout << "\n\ntest DM_XC - DM in ABACUS: \n" << std::endl; + // double DM_XC_minus_DMtest = 0.0; + // for(int ik=0; iknloc; ++iloc) + // { + // double test = std::abs(DM_XC[ik][iloc] - dmk_pointer[iloc]); + // DM_XC_minus_DMtest += test; + // if( test > 1e-16 ) + // { + // std::cout << "\nik, iloc, minus[ik][iloc]: " << ik << " " << iloc << " " << test << std::endl; + // } + // } + // } + // std::cout << "\nsum of DM_XC - DM in ABACUS: " << DM_XC_minus_DMtest << std::endl; + + if( !only_exx_type ) + { + HR_dft_XC->set_zero(); + if( PARAM.inp.gamma_only ) + { + // this can be optimized, use potXC.update_from_charge() + V_dft_XC = new rdmft::Veff_rdmft( + GG, + hsk_dft_XC, + kv->kvec_d, + this->pelec->pot, + HR_dft_XC, + &GlobalC::ucell, + orb->cutoffs(), + &GlobalC::GridD, + nspin, + + charge, + rho_basis, + vloc, + sf, + "xc", + &etxc, + &vtxc + ); + } + else + { + // this can be optimized, use potXC.update_from_charge() + V_dft_XC = new rdmft::Veff_rdmft( + GK, + hsk_dft_XC, + kv->kvec_d, + this->pelec->pot, + HR_dft_XC, + &GlobalC::ucell, + orb->cutoffs(), + &GlobalC::GridD, + nspin, + + charge, + rho_basis, + vloc, + sf, + "xc", + &etxc, + &vtxc + ); + } + V_dft_XC->contributeHR(); + } + +#ifdef __EXX + if(GlobalC::exx_info.info_global.cal_exx) + { + HR_exx_XC->set_zero(); + + std::vector< std::vector > DM_XC(nk_total, std::vector(ParaV->nloc)); + get_DM_XC(DM_XC); + // get DM_XC of all k points + if( exx_spacegroup_symmetry ) + { + DM_XC = symrot_exx.restore_dm(*this->kv, DM_XC, *ParaV); // class vector could be auto resize() + } + std::vector< const std::vector* > DM_XC_pointer(DM_XC.size()); + for(int ik=0; ik>,RI::Tensor>>> + Ds_XC_d = std::is_same::value //gamma_only_local + ? RI_2D_Comm::split_m2D_ktoR(*kv, DM_XC_pointer, *ParaV, nspin) + : RI_2D_Comm::split_m2D_ktoR(*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry); + + // provide the Ds_XC to Vxc_fromRI(V_exx_XC) + if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) + { + Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV, &this->symrot_exx); + } + else + { + Vxc_fromRI_d->cal_exx_elec(Ds_XC_d, *ParaV); + } + + // when we doing V_exx_XC.contributeHk(ik), we get HK_XC constructed by the special DM_XC + V_exx_XC = new hamilt::OperatorEXX>( + hsk_exx_XC, + HR_exx_XC, + *kv, + &Vxc_fromRI_d->Hexxs, + nullptr, + hamilt::Add_Hexx_Type::k + ); + } + else + { + // transfer the DM_XC to appropriate format + std::vector>,RI::Tensor>>>> + Ds_XC_c = std::is_same::value //gamma_only_local + ? RI_2D_Comm::split_m2D_ktoR>(*kv, DM_XC_pointer, *ParaV, nspin) + : RI_2D_Comm::split_m2D_ktoR>(*kv, DM_XC_pointer, *ParaV, nspin, this->exx_spacegroup_symmetry); + + // // provide the Ds_XC to Vxc_fromRI(V_exx_XC) + if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) + { + Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV, &this->symrot_exx); + } + else + { + Vxc_fromRI_c->cal_exx_elec(Ds_XC_c, *ParaV); + } + + // when we doing V_exx_XC.contributeHk(ik), we get HK_XC constructed by the special DM_XC + V_exx_XC = new hamilt::OperatorEXX>( + hsk_exx_XC, + HR_exx_XC, + *kv, + nullptr, + &Vxc_fromRI_c->Hexxs, + hamilt::Add_Hexx_Type::k + ); + } + // use hamilt::Add_Hexx_Type::k, not R, contributeHR() should be skipped + // V_exx_XC->contributeHR(); + } +#endif + +} + + + + +template class RDMFT; +template class RDMFT, double>; +template class RDMFT, std::complex>; + +} diff --git a/source/module_rdmft/rdmft.cpp b/source/module_rdmft/rdmft.cpp new file mode 100644 index 0000000000..c93d270d1e --- /dev/null +++ b/source/module_rdmft/rdmft.cpp @@ -0,0 +1,459 @@ +//========================================================== +// Author: Jingang Han +// DATE : 2024-03-11 +//========================================================== +#include "rdmft.h" +#include "module_rdmft/rdmft_tools.h" + +#include "module_base/blas_connector.h" +#include "module_base/scalapack_connector.h" +#include "module_base/timer.h" +// #include "module_psi/psi.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" + +// #include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_cell/module_symmetry/symmetry.h" +// #include "module_hamilt_general/module_xc/xc_functional.h" + +#include +#include +#include +#include +#include + +// used by class Veff_rdmft +//#include "veff_lcao.h" +//#include "module_base/timer.h" +// #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h" +// #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" +#include "module_base/tool_title.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_elecstate/potentials/H_Hartree_pw.h" +#include "module_elecstate/potentials/pot_local.h" +#include "module_elecstate/potentials/pot_xc.h" +#include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" +#include "module_elecstate/elecstate_lcao.h" + +// for test use dgemm_ +#include "module_base/matrix.h" +#include "module_base/blas_connector.h" + +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" //test + +// #include "module_elecstate/module_charge/symmetry_rho.h" + + + + +namespace rdmft +{ + + +template +RDMFT::RDMFT() +{ + +} + +template +RDMFT::~RDMFT() +{ + delete HR_TV; + delete HR_hartree; + // delete HR_XC; + delete HR_exx_XC; + // delete HR_local; + delete HR_dft_XC; +#ifdef __EXX + delete Vxc_fromRI_d; + delete Vxc_fromRI_c; +#endif + delete V_ekinetic_potential; + delete V_nonlocal; + delete V_local; + delete V_hartree; + // delete V_XC; + delete V_exx_XC; + delete V_dft_XC; + delete V_hartree_XC; +} + +template +void RDMFT::init(Gint_Gamma& GG_in, Gint_k& GK_in, Parallel_Orbitals& ParaV_in, UnitCell& ucell_in, + K_Vectors& kv_in, elecstate::ElecState& pelec_in, LCAO_Orbitals& orb_in, TwoCenterBundle& two_center_bundle_in, std::string XC_func_rdmft_in, double alpha_power_in) +{ + GG = &GG_in; + GK = &GK_in; + ParaV = &ParaV_in; + ucell = &ucell_in; + kv = &kv_in; + charge = pelec_in.charge; + pelec = &pelec_in; + orb = &orb_in; + two_center_bundle = &two_center_bundle_in; + XC_func_rdmft = XC_func_rdmft_in; + alpha_power = alpha_power_in; + + nspin = PARAM.inp.nspin; + nbands_total = PARAM.inp.nbands; + nk_total = ModuleSymmetry::Symmetry::symm_flag == -1 ? kv->get_nkstot_full(): kv->get_nks(); + nk_total *= nspin; + only_exx_type = ( XC_func_rdmft == "hf" || XC_func_rdmft == "muller" || XC_func_rdmft == "power" ); + + // // create desc[] and something about MPI to Eij(nbands*nbands) +#ifdef __MPI + para_Eij.set(nbands_total, nbands_total, ParaV->nb, ParaV->blacs_ctxt); // maybe in default, PARAM.inp.nb2d = 0, can't be used +#endif + + // + occ_number.create(nk_total, nbands_total); + wg.create(nk_total, nbands_total); + wk_fun_occNum.create(nk_total, nbands_total); + occNum_wfcHamiltWfc.create(nk_total, nbands_total); + Etotal_n_k.create(nk_total, nbands_total); + wfcHwfc_TV.create(nk_total, nbands_total); + wfcHwfc_hartree.create(nk_total, nbands_total); + wfcHwfc_XC.create(nk_total, nbands_total); + wfcHwfc_exx_XC.create(nk_total, nbands_total); + wfcHwfc_dft_XC.create(nk_total, nbands_total); + + // + wfc.resize(nk_total, ParaV->ncol_bands, ParaV->nrow); // test ParaV->nrow + occNum_HamiltWfc.resize(nk_total, ParaV->ncol_bands, ParaV->nrow); + H_wfc_TV.resize(nk_total, ParaV->ncol_bands, ParaV->nrow); + H_wfc_hartree.resize(nk_total, ParaV->ncol_bands, ParaV->nrow); + H_wfc_XC.resize(nk_total, ParaV->ncol_bands, ParaV->nrow); + H_wfc_exx_XC.resize(nk_total, ParaV->ncol_bands, ParaV->nrow); + H_wfc_dft_XC.resize(nk_total, ParaV->ncol_bands, ParaV->nrow); + + // + hsk_TV = new hamilt::HS_Matrix_K(ParaV, true); + hsk_hartree = new hamilt::HS_Matrix_K(ParaV, true); + hsk_dft_XC = new hamilt::HS_Matrix_K(ParaV, true); + hsk_exx_XC = new hamilt::HS_Matrix_K(ParaV, true); + + HK_XC.resize( ParaV->get_row_size()*ParaV->get_col_size() ); + // HK_RDMFT_pass.resize(nk_total, ParaV->get_row_size(), ParaV->get_col_size()); + // HK_XC_pass.resize(nk_total, ParaV->get_row_size(), ParaV->get_col_size()); + + + Eij_TV.resize( para_Eij.get_row_size()*para_Eij.get_col_size() ); + Eij_hartree.resize( para_Eij.get_row_size()*para_Eij.get_col_size() ); + Eij_XC.resize( para_Eij.get_row_size()*para_Eij.get_col_size() ); + Eij_exx_XC.resize( para_Eij.get_row_size()*para_Eij.get_col_size() ); + + // + HR_TV = new hamilt::HContainer(*ucell, ParaV); + HR_hartree = new hamilt::HContainer(*ucell, ParaV); + // HR_XC = new hamilt::HContainer(*ucell, ParaV); + HR_exx_XC = new hamilt::HContainer(*ucell, ParaV); + HR_dft_XC = new hamilt::HContainer(*ucell, ParaV); + // HR_local = new hamilt::HContainer(*ucell, ParaV); + + // set zero ( std::vector, ModuleBase::matrix will automatically be set to zero ) + wfc.zero_out(); + occNum_HamiltWfc.zero_out(); + H_wfc_TV.zero_out(); + H_wfc_hartree.zero_out(); + H_wfc_XC.zero_out(); + H_wfc_exx_XC.zero_out(); + H_wfc_dft_XC.zero_out(); + + HR_TV->set_zero(); // HR->set_zero() might be delete here, test on Gamma_only in the furure + HR_hartree->set_zero(); + // HR_XC->set_zero(); + HR_exx_XC->set_zero(); + HR_dft_XC->set_zero(); + // HR_local->set_zero(); + +#ifdef __EXX + if( GlobalC::exx_info.info_global.cal_exx ) + { + // if the irreducible k-points can change with symmetry during cell-relax, it should be moved back to update_ion() + exx_spacegroup_symmetry = (PARAM.inp.nspin < 4 && ModuleSymmetry::Symmetry::symm_flag == 1); + if (exx_spacegroup_symmetry) + { + const std::array& period = RI_Util::get_Born_vonKarmen_period(*kv); + this->symrot_exx.find_irreducible_sector(ucell->symm, ucell->atoms, ucell->st, + RI_Util::get_Born_von_Karmen_cells(period), period, ucell->lat); + this->symrot_exx.cal_Ms(*kv, *ucell, *ParaV); + } + + if (GlobalC::exx_info.info_ri.real_number) + { + Vxc_fromRI_d = new Exx_LRI(GlobalC::exx_info.info_ri); + Vxc_fromRI_d->init(MPI_COMM_WORLD, *kv, *orb); + } + else + { + Vxc_fromRI_c = new Exx_LRI>(GlobalC::exx_info.info_ri); + Vxc_fromRI_c->init(MPI_COMM_WORLD, *kv, *orb); + } + } +#endif + + if( PARAM.inp.gamma_only ) + { + HR_TV->fix_gamma(); + HR_hartree->fix_gamma(); + // HR_XC->fix_gamma(); + HR_exx_XC->fix_gamma(); + HR_dft_XC->fix_gamma(); + } + +} + + +template +void RDMFT::cal_Hk_Hpsi() +{ + /****** get occNum_wfcHamiltWfc, occNum_HamiltWfc ******/ + // HK_RDMFT_pass.reset(); + + // double XC_minus_XC = 0.0; + // std::cout << "\n\ntest V_exx_XC in rdmft.cpp: " << std::endl; + // HK_XC_pass.reset(); + + //calculate Hwfc, wfcHwfc for each potential + for(int ik=0; ikset_zero_hk(); + hsk_hartree->set_zero_hk(); + set_zero_vector(HK_XC); + + // get the HK with ik-th k vector, the result is stored in HK_TV, HK_hartree and HK_XC respectively + V_local->contributeHk(ik); + V_hartree->contributeHk(ik); + + // get H(k) * wfc + HkPsi( ParaV, hsk_TV->get_hk()[0], wfc(ik, 0, 0), H_wfc_TV(ik, 0, 0)); + HkPsi( ParaV, hsk_hartree->get_hk()[0], wfc(ik, 0, 0), H_wfc_hartree(ik, 0, 0)); + + // get wfc * H(k)_wfc + psiDotPsi( ParaV, para_Eij, wfc(ik, 0, 0), H_wfc_TV(ik, 0, 0), Eij_TV, &(wfcHwfc_TV(ik, 0)) ); + psiDotPsi( ParaV, para_Eij, wfc(ik, 0, 0), H_wfc_hartree(ik, 0, 0), Eij_hartree, &(wfcHwfc_hartree(ik, 0)) ); + +#ifdef __EXX + if(GlobalC::exx_info.info_global.cal_exx) + { + // set_zero_vector(HK_exx_XC); + hsk_exx_XC->set_zero_hk(); + + V_exx_XC->contributeHk(ik); + HkPsi( ParaV, hsk_exx_XC->get_hk()[0], wfc(ik, 0, 0), H_wfc_exx_XC(ik, 0, 0)); + psiDotPsi( ParaV, para_Eij, wfc(ik, 0, 0), H_wfc_exx_XC(ik, 0, 0), Eij_exx_XC, &(wfcHwfc_exx_XC(ik, 0)) ); + + for(int iloc=0; ilocget_hk()[iloc]; + } +#endif + if( !only_exx_type ) + { + // set_zero_vector(HK_dft_XC); + hsk_dft_XC->set_zero_hk(); + + V_dft_XC->contributeHk(ik); + HkPsi( ParaV, hsk_dft_XC->get_hk()[0], wfc(ik, 0, 0), H_wfc_dft_XC(ik, 0, 0)); + psiDotPsi( ParaV, para_Eij, wfc(ik, 0, 0), H_wfc_dft_XC(ik, 0, 0), Eij_exx_XC, &(wfcHwfc_dft_XC(ik, 0)) ); + + for(int iloc=0; ilocget_hk()[iloc]; + } + + // // store HK_RDMFT + // for(int ir=0; irget_col_size() + ir] + // + HK_hartree[ic * ParaV->get_col_size() + ir] + // + HK_XC[ic * ParaV->get_col_size() + ir]; + // // HK_XC_pass[ik](ir, ic) = HK_XC[ic * ParaV->get_col_size() + ir]; + // } + // } + + // using them to the gradient of Etotal is not correct when do hybrid calculation, it's correct just for exx-type functional + // HkPsi( ParaV, HK_XC[0], wfc(ik, 0, 0), H_wfc_XC(ik, 0, 0)); + // psiDotPsi( ParaV, para_Eij, wfc(ik, 0, 0), H_wfc_XC(ik, 0, 0), Eij_XC, &(wfcHwfc_XC(ik, 0)) ); + + } + + // std::cout << "\n\nsum of XC_minus_XC: " << XC_minus_XC << "\n\n" << std::endl; + +} + + +template +double RDMFT::cal_E_gradient() +{ + /****** get occNum_wfcHamiltWfc, occNum_HamiltWfc and Etotal ******/ + + // !this would transfer the value of H_wfc_TV, H_wfc_hartree, H_wfc_XC --> occNum_H_wfc + // get the gradient of energy with respect to the wfc, i.e., Wk_occNum_HamiltWfc + add_psi(ParaV, kv, occ_number, H_wfc_TV, H_wfc_hartree, H_wfc_dft_XC, H_wfc_exx_XC, occNum_HamiltWfc, XC_func_rdmft, alpha_power); + + // get the gradient of energy with respect to the natural occupation numbers, i.e., Wk_occNum_wfcHamiltWfc + add_occNum(*kv, occ_number, wfcHwfc_TV, wfcHwfc_hartree, wfcHwfc_dft_XC, wfcHwfc_exx_XC, occNum_wfcHamiltWfc, XC_func_rdmft, alpha_power); + + // get the total energy + // add_wfcHwfc(kv->wk, occ_number, wfcHwfc_TV, wfcHwfc_hartree, wfcHwfc_XC, Etotal_n_k, XC_func_rdmft, alpha_power); + // add_wfcHwfc(wg, wk_fun_occNum, wfcHwfc_TV, wfcHwfc_hartree, wfcHwfc_XC, Etotal_n_k, XC_func_rdmft, alpha_power); + // E_RDMFT[3] = getEnergy(Etotal_n_k); + // Parallel_Reduce::reduce_all(E_RDMFT[3]); + + return E_RDMFT[3]; + + /****** get occNum_wfcHamiltWfc, occNum_HamiltWfc and Etotal ******/ + +} + + +// cal_type = 2 just support XC-functional without exx +template +void RDMFT::cal_Energy(const int cal_type) +{ + double E_Ewald = pelec->f_en.ewald_energy; + double E_entropy = pelec->f_en.demet; + double E_descf = pelec->f_en.descf = 0.0; + // double E_descf = 0.0; + double E_xc_KS = pelec->f_en.etxc - pelec->f_en.etxcc; + double E_exx_KS = pelec->f_en.exx; + double E_deband_KS = pelec->f_en.deband; + double E_deband_harris_KS = pelec->f_en.deband_harris; + + double E_exxType_rdmft = 0.0; // delete in the future + + if( cal_type == 1 ) + { + // for E_TV + ModuleBase::matrix ETV_n_k(wg.nr, wg.nc, true); + occNum_Mul_wfcHwfc(wg, wfcHwfc_TV, ETV_n_k, 0); + E_RDMFT[0] = getEnergy(ETV_n_k); + + // for Ehartree + ModuleBase::matrix Ehartree_n_k(wg.nr, wg.nc, true); + occNum_Mul_wfcHwfc(wg, wfcHwfc_hartree, Ehartree_n_k, 1); + E_RDMFT[1] = getEnergy(Ehartree_n_k); + + // for Exc + E_RDMFT[2] = 0.0; +#ifdef __EXX + if( GlobalC::exx_info.info_global.cal_exx ) + { + ModuleBase::matrix Exc_n_k(wg.nr, wg.nc, true); + // because we have got wk_fun_occNum, we can use symbol=1 realize it + occNum_Mul_wfcHwfc(wk_fun_occNum, wfcHwfc_exx_XC, Exc_n_k, 1); + E_RDMFT[2] = getEnergy(Exc_n_k); + Parallel_Reduce::reduce_all(E_RDMFT[2]); + E_exxType_rdmft = E_RDMFT[2]; + } +#endif + E_RDMFT[2] += etxc; + + // add up the results obtained by all processors, or we can do reduce_all(wfcHwfc_) before add_wg() used for Etotal to replace it + Parallel_Reduce::reduce_all(E_RDMFT[0]); + Parallel_Reduce::reduce_all(E_RDMFT[1]); + + Etotal = E_RDMFT[0] + E_RDMFT[1] + E_RDMFT[2] + E_Ewald + E_entropy + E_descf; + + // temp + E_RDMFT[3] = E_RDMFT[0] + E_RDMFT[1] + E_RDMFT[2]; + } + else + { + this->pelec->f_en.deband = this->pelec->cal_delta_eband(); + E_descf = pelec->f_en.descf = 0.0; + this->pelec->cal_energies(2); + Etotal = this->pelec->f_en.etot; + + // if( GlobalC::exx_info.info_global.cal_exx ) + // { + // ModuleBase::matrix Exc_n_k(wg.nr, wg.nc, true); + // // because we have got wk_fun_occNum, we can use symbol=1 realize it + // occNum_Mul_wfcHwfc(wk_fun_occNum, wfcHwfc_XC, Exc_n_k, 1); + // E_RDMFT[2] = getEnergy(Exc_n_k); + // Parallel_Reduce::reduce_all(E_RDMFT[2]); + + // // test + // Etotal -= E_RDMFT[2]; + // } + } + +// // print results +// std::cout << "\n\nfrom class RDMFT: \nXC_fun: " << XC_func_rdmft << std::endl; +// #ifdef __EXX +// if( GlobalC::exx_info.info_global.cal_exx ) std::cout << "alpha_power: " << alpha_power << std::endl; +// #endif +// std::cout << std::fixed << std::setprecision(10) +// << "******\nE(TV + Hartree + XC) by RDMFT: " << E_RDMFT[3] +// << "\n\nE_TV_RDMFT: " << E_RDMFT[0] +// << "\nE_hartree_RDMFT: " << E_RDMFT[1] +// << "\nExc_" << XC_func_rdmft << "_RDMFT: " << E_RDMFT[2] +// << "\nE_Ewald: " << E_Ewald +// << "\nE_entropy(-TS): " << E_entropy +// << "\nE_descf: " << E_descf +// << "\n\nEtotal_RDMFT: " << Etotal +// << "\n\nExc_ksdft: " << E_xc_KS +// << "\nE_exx_ksdft: " << E_exx_KS +// <<"\n******\n\n" << std::endl; + +// std::cout << "\netxc: " << etxc << "\nvtxc: " << vtxc << "\n"; +// std::cout << "\nE_deband_KS: " << E_deband_KS << "\nE_deband_harris_KS: " << E_deband_harris_KS << "\n\n" << std::endl; + + if( PARAM.inp.rdmft == true ) + { + GlobalV::ofs_running << "\n\nfrom class RDMFT: \nXC_fun: " << XC_func_rdmft << std::endl; +#ifdef __EXX + if( GlobalC::exx_info.info_global.cal_exx ) { GlobalV::ofs_running << "alpha_power: " << alpha_power << std::endl; +} +#endif + // GlobalV::ofs_running << std::setprecision(12); + // GlobalV::ofs_running << std::setiosflags(std::ios::right); + GlobalV::ofs_running << std::fixed << std::setprecision(10) + << "\n******\nE(TV + Hartree + XC) by RDMFT: " << E_RDMFT[3] + << "\n\nE_TV_RDMFT: " << E_RDMFT[0] + << "\nE_hartree_RDMFT: " << E_RDMFT[1] + << "\nExc_" << XC_func_rdmft << "_RDMFT: " << E_RDMFT[2] + << "\nE_Ewald: " << E_Ewald + << "\nE_entropy(-TS): " << E_entropy + << "\nE_descf: " << E_descf + << "\n\nEtotal_RDMFT: " << Etotal + << "\n\nExc_ksdft: " << E_xc_KS + << "\nE_exx_ksdft: " << E_exx_KS + << "\nE_exxType_rdmft: " << E_exxType_rdmft + <<"\n******\n" << std::endl; + } + std::cout << std::defaultfloat; + +} + + +template +double RDMFT::run(ModuleBase::matrix& E_gradient_occNum, psi::Psi& E_gradient_wfc) +{ + ModuleBase::TITLE("RDMFT", "E & Egradient"); + ModuleBase::timer::tick("RDMFT", "E & Egradient"); + + // this->cal_V_hartree(); + // this->cal_V_XC(); + this->cal_Hk_Hpsi(); + this->cal_E_gradient(); + this->cal_Energy(this->cal_E_type); + // this->cal_Energy(2); + + E_gradient_occNum = (occNum_wfcHamiltWfc); + + TK* pwfc = &occNum_HamiltWfc(0, 0, 0); + TK* pwfc_out = &E_gradient_wfc(0, 0, 0); + for(int i=0; i; +template class RDMFT, double>; +template class RDMFT, std::complex>; + +} + + diff --git a/source/module_rdmft/rdmft.h b/source/module_rdmft/rdmft.h new file mode 100644 index 0000000000..3b75c5f071 --- /dev/null +++ b/source/module_rdmft/rdmft.h @@ -0,0 +1,249 @@ +//========================================================== +// Author: Jingang Han +// DATE : 2024-03-11 +//========================================================== +#ifndef RDMFT_H +#define RDMFT_H + +// #include "module_rdmft/rdmft_tools.h" + +#include "module_parameter/parameter.h" + +#include "module_psi/psi.h" +#include "module_base/matrix.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/unitcell.h" +#include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" +#include "module_elecstate/potentials/potential_new.h" +#include "module_base/blas_connector.h" +#include "module_base/scalapack_connector.h" +#include "module_basis/module_ao/parallel_2d.h" +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_base/parallel_reduce.h" +// #include "module_elecstate/module_dm/cal_dm_psi.h" +// #include "module_elecstate/module_dm/density_matrix.h" +#include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" + +#include "module_basis/module_ao/ORB_read.h" +#include "module_basis/module_nao/two_center_bundle.h" + +#include "module_hamilt_general/operator.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h" + + +// used by Exx&LRI +#ifdef __EXX +#include "module_ri/Exx_LRI.h" +#include "module_ri/RI_2D_Comm.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" +#include "module_ri/module_exx_symmetry/symmetry_rotation.h" +#endif + +// there are some operator reload to print data in different formats +#include "module_ri/test_code/test_function.h" + +// used by class Veff_rdmft +#include "module_base/timer.h" +#include "module_elecstate/potentials/potential_new.h" +// #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" +//#include "module_hamilt_lcao/module_gint/gint_gamma.h" +//#include "module_hamilt_lcao/module_gint/gint_k.h" +//#include "operator_lcao.h" +//#include "module_cell/module_neighbor/sltk_grid_driver.h" +//#include "module_cell/unitcell.h" +#include "module_elecstate/potentials/H_Hartree_pw.h" +#include "module_elecstate/potentials/pot_local.h" +#include "module_elecstate/potentials/pot_xc.h" +#include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" + + +#include +#include +#include +#include +#include + + + +namespace rdmft +{ + + +template +class RDMFT +{ + public: + RDMFT(); + ~RDMFT(); + + + /****** these parameters are passed in from outside, don't need delete ******/ + Parallel_Orbitals* ParaV = nullptr; + Parallel_2D para_Eij; + // Parallel_Orbitals para_Eij; + + // GK and GG are used for multi-k grid integration and gamma only algorithms respectively + Gint_k* GK = nullptr; + Gint_Gamma* GG = nullptr; + Charge* charge = nullptr; + elecstate::ElecState* pelec = nullptr; // just to gain Ewald and this->pelec->pot + + // update after ion step + UnitCell* ucell = nullptr; + K_Vectors* kv = nullptr; + // LCAO_Matrix* LM = nullptr; + ModulePW::PW_Basis* rho_basis = nullptr; + ModuleBase::matrix* vloc = nullptr; + ModuleBase::ComplexMatrix* sf = nullptr; + LCAO_Orbitals* orb = nullptr; + TwoCenterBundle* two_center_bundle = nullptr; + // Local_Orbital_Charge* loc = nullptr; // would be delete in the future + /****** these parameters are passed in from outside, don't need delete ******/ + + int nk_total = 0; + int nbands_total; + int nspin = 1; + std::string XC_func_rdmft; + double alpha_power = 0.656; // 0.656 for soilds, 0.525 for dissociation of H2, 0.55~0.58 for HEG + + // natrual occupation numbers and wavefunction + ModuleBase::matrix occ_number; + psi::Psi wfc; + ModuleBase::matrix wg; + ModuleBase::matrix wk_fun_occNum; + + // store the gradients of Etotal with respect to the natural occupation numbers and wfc respectively + ModuleBase::matrix occNum_wfcHamiltWfc; + psi::Psi occNum_HamiltWfc; + + // E_RDMFT[4] stores ETV, Ehartree, Exc, Etotal respectively + double E_RDMFT[4] = {0.0}; + // std::vector E_RDMFT(4); + + hamilt::HContainer* HR_TV = nullptr; + hamilt::HContainer* HR_hartree = nullptr; + // hamilt::HContainer* HR_XC = nullptr; + hamilt::HContainer* HR_dft_XC = nullptr; + hamilt::HContainer* HR_exx_XC = nullptr; + // hamilt::HContainer* HR_local = nullptr; + + hamilt::HS_Matrix_K* hsk_TV; + hamilt::HS_Matrix_K* hsk_hartree; + hamilt::HS_Matrix_K* hsk_dft_XC; + hamilt::HS_Matrix_K* hsk_exx_XC; + + std::vector HK_XC; + std::vector< std::vector > DM_XC_pass; + // ModuleDirectMin::ProdStiefelVariable HK_RDMFT_pass; + // ModuleDirectMin::ProdStiefelVariable HK_XC_pass; + + ModuleBase::matrix Etotal_n_k; + ModuleBase::matrix wfcHwfc_TV; + ModuleBase::matrix wfcHwfc_hartree; + ModuleBase::matrix wfcHwfc_XC; + ModuleBase::matrix wfcHwfc_exx_XC; + ModuleBase::matrix wfcHwfc_dft_XC; + + psi::Psi H_wfc_TV; + psi::Psi H_wfc_hartree; + psi::Psi H_wfc_XC; + psi::Psi H_wfc_exx_XC; + psi::Psi H_wfc_dft_XC; + + // just for temperate. in the future when realize psiDotPsi() without pzgemm_/pdgemm_,we don't need it + std::vector Eij_TV; + std::vector Eij_hartree; + std::vector Eij_XC; + std::vector Eij_exx_XC; + + hamilt::OperatorLCAO* V_ekinetic_potential = nullptr; + hamilt::OperatorLCAO* V_nonlocal = nullptr; + hamilt::OperatorLCAO* V_local = nullptr; + hamilt::OperatorLCAO* V_hartree = nullptr; + // hamilt::OperatorLCAO* V_XC = nullptr; + hamilt::OperatorLCAO* V_exx_XC = nullptr; + hamilt::OperatorLCAO* V_dft_XC = nullptr; + hamilt::OperatorLCAO* V_hartree_XC = nullptr; + // bool get_V_local_temp = true; + +#ifdef __EXX + Exx_LRI* Vxc_fromRI_d = nullptr; + Exx_LRI>* Vxc_fromRI_c = nullptr; + ModuleSymmetry::Symmetry_rotation symrot_exx; + bool exx_spacegroup_symmetry = false; +#endif + + double Etotal = 0.0; + double etxc = 0.0; + double vtxc = 0.0; + bool only_exx_type = false; + const int cal_E_type = 1; // cal_type = 2 just support XC-functional without exx + + void init(Gint_Gamma& GG_in, Gint_k& GK_in, Parallel_Orbitals& ParaV_in, UnitCell& ucell_in, + K_Vectors& kv_in, elecstate::ElecState& pelec_in, LCAO_Orbitals& orb_in, TwoCenterBundle& two_center_bundle_in, std::string XC_func_rdmft_in, double alpha_power_in); + + // update in ion-step and get V_TV + void update_ion(UnitCell& ucell_in, ModulePW::PW_Basis& rho_basis_in, + ModuleBase::matrix& vloc_in, ModuleBase::ComplexMatrix& sf_in); + + // update in elec-step + // Or we can use rdmft_solver.wfc/occ_number directly when optimizing, so that the update_elec() function does not require parameters. + void update_elec(const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in = nullptr); + + // update occ_number for optimization algorithms that depend on Hamilton + void update_occNumber(const ModuleBase::matrix& occ_number_in); + + // update occ_number for optimization algorithms that depend on Hamilton + void update_wg(const ModuleBase::matrix& wg_in); + + // do all calculation after update occNum&wfc, get Etotal and the gradient of energy with respect to the occNum&wfc + double run(ModuleBase::matrix& E_gradient_occNum, psi::Psi& E_gradient_wfc); + + + + protected: + + // get the special density matrix DM_XC(nk*nbasis_local*nbasis_local) + void get_DM_XC(std::vector< std::vector >& DM_XC); + + void cal_V_TV(); + + void cal_V_hartree(); + + // construct V_XC based on different XC_functional( i.e. RDMFT class member XC_func_rdmft) + void cal_V_XC(); + + double cal_E_gradient(); + + void cal_Energy(const int cal_type = 1); + + + + private: + + // get the total Hamilton in k-space + void cal_Hk_Hpsi(); + + void update_charge(); + + + +}; + + + + + + +} + + + + +#endif \ No newline at end of file diff --git a/source/module_rdmft/rdmft_tools.cpp b/source/module_rdmft/rdmft_tools.cpp new file mode 100644 index 0000000000..cce3726854 --- /dev/null +++ b/source/module_rdmft/rdmft_tools.cpp @@ -0,0 +1,439 @@ +//========================================================== +// Author: Jingang Han +// DATE : 2024-03-11 +//========================================================== +// #include "rdmft.h" +#include "module_rdmft/rdmft_tools.h" + +#include "module_base/blas_connector.h" +#include "module_base/scalapack_connector.h" +#include "module_base/timer.h" +#include "module_psi/psi.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" + +#include "module_elecstate/module_dm/cal_dm_psi.h" + +#include +#include +#include +#include +#include + +// used by class Veff_rdmft +//#include "veff_lcao.h" +//#include "module_base/timer.h" +#include "module_base/tool_title.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_elecstate/potentials/H_Hartree_pw.h" +#include "module_elecstate/potentials/pot_local.h" +#include "module_elecstate/potentials/pot_xc.h" +#include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" + +// for test use dgemm_ +#include "module_base/matrix.h" +#include "module_base/blas_connector.h" + +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" //test + + + + + +namespace rdmft +{ + + + + + + +template <> +void conj_psi(psi::Psi& wfc) {} + + +template <> +void HkPsi(const Parallel_Orbitals* ParaV, const double& HK, const double& wfc, double& H_wfc) +{ + const int one_int = 1; + const double one_double = 1.0; + const double zero_double = 0.0; + const char N_char = 'N'; + const char C_char = 'C'; + +#ifdef __MPI + const int nbasis = ParaV->desc[2]; + const int nbands = ParaV->desc_wfc[3]; + + //because wfc(bands, basis'), H(basis, basis'), we do wfc*H^T(in the perspective of cpp, not in fortran). And get H_wfc(bands, basis) is correct. + pdgemm_( &C_char, &N_char, &nbasis, &nbands, &nbasis, &one_double, &HK, &one_int, &one_int, ParaV->desc, + &wfc, &one_int, &one_int, ParaV->desc_wfc, &zero_double, &H_wfc, &one_int, &one_int, ParaV->desc_wfc ); +#endif + +} + + +template <> +void psiDotPsi(const Parallel_Orbitals* ParaV, const Parallel_2D& para_Eij_in, + const double& wfc, const double& H_wfc, std::vector& Dmn, double* wfcHwfc) +{ + const int one_int = 1; + const double one_double = 1.0; + const double zero_double = 0.0; + const char N_char = 'N'; + const char T_char = 'T'; + + const int nrow_bands = para_Eij_in.get_row_size(); + const int ncol_bands = para_Eij_in.get_col_size(); + +#ifdef __MPI + const int nbasis = ParaV->desc[2]; + const int nbands = ParaV->desc_wfc[3]; + + pdgemm_( &T_char, &N_char, &nbands, &nbands, &nbasis, &one_double, &wfc, &one_int, &one_int, ParaV->desc_wfc, + &H_wfc, &one_int, &one_int, ParaV->desc_wfc, &zero_double, &Dmn[0], &one_int, &one_int, para_Eij_in.desc ); +#endif + + for(int i=0; i; + +template class Veff_rdmft, double>; + +template class Veff_rdmft, std::complex>; + +// this part of the code is copying from class Veff +// initialize_HR() +template +void Veff_rdmft::initialize_HR(const UnitCell* ucell_in, + Grid_Driver* GridD) +{ + ModuleBase::TITLE("Veff", "initialize_HR"); + ModuleBase::timer::tick("Veff", "initialize_HR"); + + this->nspin = PARAM.inp.nspin; + auto* paraV = this->hR->get_paraV();// get parallel orbitals from HR + // TODO: if paraV is nullptr, AtomPair can not use paraV for constructor, I will repair it in the future. + + for (int iat1 = 0; iat1 < ucell_in->nat; iat1++) + { + auto tau1 = ucell_in->get_tau(iat1); + int T1 = 0; + int I1 = 0; + ucell_in->iat2iait(iat1, &I1, &T1); + AdjacentAtomInfo adjs; + GridD->Find_atom(*ucell_in, tau1, T1, I1, &adjs); + std::vector is_adj(adjs.adj_num + 1, false); + for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) + { + const int T2 = adjs.ntype[ad1]; + const int I2 = adjs.natom[ad1]; + const int iat2 = ucell_in->itia2iat(T2, I2); + if (paraV->get_row_size(iat1) <= 0 || paraV->get_col_size(iat2) <= 0) + { + continue; + } + const ModuleBase::Vector3& R_index2 = adjs.box[ad1]; + // choose the real adjacent atoms + // Note: the distance of atoms should less than the cutoff radius, + // When equal, the theoretical value of matrix element is zero, + // but the calculated value is not zero due to the numerical error, which would lead to result changes. + if (ucell_in->cal_dtau(iat1, iat2, R_index2).norm() * ucell_in->lat0 + < orb_cutoff_[T1] + orb_cutoff_[T2]) + { + hamilt::AtomPair tmp(iat1, iat2, R_index2, paraV); + this->hR->insert_pair(tmp); + } + } + } + // allocate the memory of BaseMatrix in HR, and set the new values to zero + this->hR->allocate(nullptr, true); + + ModuleBase::timer::tick("Veff", "initialize_HR"); +} + + +// this part of the code is copying from class Veff and do some modifications. +template +void Veff_rdmft::contributeHR() +{ + ModuleBase::TITLE("Veff", "contributeHR"); + ModuleBase::timer::tick("Veff", "contributeHR"); + + this->GK->reset_spin(this->current_spin); + + double* vr_eff_rdmft = nullptr; + + // calculate v_hartree(r) or v_local(r) or v_xc(r) + if( potential_ == "hartree" ) + { + ModuleBase::matrix v_matrix_hartree(this->nspin, charge_->nrxx); + elecstate::PotHartree potH(rho_basis_); + potH.cal_v_eff(charge_, ucell, v_matrix_hartree); + + for(int is=0; isnspin; ++is) + { + // use pointer to attach v(r) for current spin + vr_eff_rdmft = &v_matrix_hartree(is, 0); + + // do grid integral calculation to get HR + Gint_inout inout(vr_eff_rdmft, is, Gint_Tools::job_type::vlocal); + this->GK->cal_gint(&inout); + } + } + else if( potential_ == "local" ) + { + double vlocal_of_0 = 0.0; + ModuleBase::matrix v_matrix_local(1, charge_->nrxx); + elecstate::PotLocal potL(vloc_, sf_, rho_basis_, vlocal_of_0); + potL.cal_fixed_v( &v_matrix_local(0, 0) ); + + // use pointer to attach v(r) + vr_eff_rdmft = &v_matrix_local(0, 0); + + // do grid integral calculation to get HR + Gint_inout inout(vr_eff_rdmft, 0, Gint_Tools::job_type::vlocal); + this->GK->cal_gint(&inout); + } + else if( potential_ == "xc" ) + { + // meta-gga type has not been considered yet !!! + + ModuleBase::matrix vofk = *vloc_; + vofk.zero_out(); + ModuleBase::matrix v_matrix_XC(this->nspin, charge_->nrxx); + elecstate::PotXC potXC(rho_basis_, etxc, vtxc, &vofk); + potXC.cal_v_eff(charge_, ucell, v_matrix_XC); + + // if need meta-GGA, go to study veff_lcao.cpp and modify the code + for(int is=0; isnspin; ++is) + { + // use pointer to attach v(r) for current spin + vr_eff_rdmft = &v_matrix_XC(is, 0); + + // do grid integral calculation to get HR + Gint_inout inout(vr_eff_rdmft, is, Gint_Tools::job_type::vlocal); + this->GK->cal_gint(&inout); + } + } + else + { + std::cout << "\n\n!!!!!!\n there may be something wrong when use class Veff_rdmft\n\n!!!!!!\n"; + } + + // get HR for 2D-block parallel format + // this->GK->transfer_pvpR(this->hR); + this->GK->transfer_pvpR(this->hR,this->ucell,this->gd); + + if(this->nspin == 2) { this->current_spin = 1 - this->current_spin; } + + ModuleBase::timer::tick("Veff", "contributeHR"); + return; +} + +// this part of the code is copying from class Veff and do some modifications. +// special case of gamma-only +template<> +void Veff_rdmft::contributeHR() +{ + ModuleBase::TITLE("Veff", "contributeHR"); + ModuleBase::timer::tick("Veff", "contributeHR"); + + // this->GK->reset_spin(this->current_spin); + + double* vr_eff_rdmft = nullptr; + + // calculate v_hartree(r) or V_local(r) or v_xc(r) + if( potential_ == "hartree" ) + { + ModuleBase::matrix v_matrix_hartree(this->nspin, charge_->nrxx); + elecstate::PotHartree potH(rho_basis_); + potH.cal_v_eff(charge_, ucell, v_matrix_hartree); + + for(int is=0; isnspin; ++is) + { + // use pointer to attach v(r) for current spin + vr_eff_rdmft = &v_matrix_hartree(is, 0); + + // do grid integral calculation to get HR + Gint_inout inout(vr_eff_rdmft, is, Gint_Tools::job_type::vlocal); + this->GG->cal_gint(&inout); + } + } + else if( potential_ == "local" ) + { + double vlocal_of_0 = 0.0; + ModuleBase::matrix v_matrix_local(1, charge_->nrxx); + elecstate::PotLocal potL(vloc_, sf_, rho_basis_, vlocal_of_0); + potL.cal_fixed_v( &v_matrix_local(0, 0) ); + + // use pointer to attach v(r) + vr_eff_rdmft = &v_matrix_local(0, 0); + + // do grid integral calculation to get HR + Gint_inout inout(vr_eff_rdmft, 0, Gint_Tools::job_type::vlocal); + + // because in gamma_only, cal_gint would not set hRGint zero first + // so must use cal_vlocal(), and in rdmft_test.h, calculate V_hartree->contributeHR() first + + this->GG->cal_vlocal(&inout, false); // cal_gint ??? + } + else if( potential_ == "xc" ) + { + // meta-gga type has not been considered yet !!! + + ModuleBase::matrix vofk = *vloc_; + vofk.zero_out(); + ModuleBase::matrix v_matrix_XC(this->nspin, charge_->nrxx); + elecstate::PotXC potXC(rho_basis_, etxc, vtxc, &vofk); + potXC.cal_v_eff(charge_, ucell, v_matrix_XC); + + for(int is=0; isnspin; ++is) + { + // use pointer to attach v(r) for current spin + vr_eff_rdmft = &v_matrix_XC(is, 0); + + // do grid integral calculation to get HR + Gint_inout inout(vr_eff_rdmft, is, Gint_Tools::job_type::vlocal); + this->GG->cal_gint(&inout); + } + } + else + { + std::cout << "\n\n!!!!!!\n there may be something wrong when use class Veff_rdmft\n\n!!!!!!\n"; + } + + // get HR for 2D-block parallel format + this->GG->transfer_pvpR(this->hR,this->ucell); + + this->new_e_iteration = false; + + if(this->nspin == 2) this->current_spin = 1 - this->current_spin; + + return; +} + + + + + + +} + + + diff --git a/source/module_rdmft/rdmft_tools.h b/source/module_rdmft/rdmft_tools.h new file mode 100644 index 0000000000..5f62badfba --- /dev/null +++ b/source/module_rdmft/rdmft_tools.h @@ -0,0 +1,458 @@ +//========================================================== +// Author: Jingang Han +// DATE : 2024-03-11 +//========================================================== +#ifndef RDMFT_TOOLS_H +#define RDMFT_TOOLS_H + +#include "module_base/matrix.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/unitcell.h" +#include "module_hamilt_lcao/module_gint/gint_gamma.h" +#include "module_hamilt_lcao/module_gint/gint_k.h" +#include "module_elecstate/potentials/potential_new.h" +#include "module_base/blas_connector.h" +#include "module_base/scalapack_connector.h" +#include "module_basis/module_ao/parallel_2d.h" +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_base/parallel_reduce.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_elecstate/module_dm/density_matrix.h" + +#include "module_hamilt_general/operator.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_new.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_new.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h" + +#include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" + +// used by Exx&LRI +#ifdef __EXX +#include "module_ri/RI_2D_Comm.h" +#include "module_ri/Exx_LRI.h" +#endif + +// there are some operator reload to print data in different formats +#include "module_ri/test_code/test_function.h" + +#include +#include +#include +#include +#include + + + +namespace rdmft +{ + + +//for print matrix +template +void printMatrix_pointer(int M, int N, const TK* matrixA, std::string nameA) +{ + std::cout << "\n" << nameA << ": \n"; + for(int i=0; i +void printMatrix_vector(int M, int N, const std::vector& matrixA, std::string nameA) +{ + std::cout << "\n" << nameA << ": \n"; + for(int i=0; i +void set_zero_vector(std::vector& HK) +{ + for(int i=0; i +void set_zero_psi(psi::Psi& wfc) +{ + TK* pwfc_in = &wfc(0, 0, 0); +#ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) +#endif + for(int i=0; i +void conj_psi(psi::Psi& wfc) +{ + TK* pwfc = &wfc(0, 0, 0); + for(int i=0; i +void conj_psi(psi::Psi& wfc); + + +// wfc and H_wfc need to be k_firest and provide wfc(ik, 0, 0) and H_wfc(ik, 0, 0) +// psi::Psi psi's (), std::vector HK's [] operator overloading return TK +template +void HkPsi(const Parallel_Orbitals* ParaV, const TK& HK, const TK& wfc, TK& H_wfc) +{ + + const int one_int = 1; + //const double one_double = 1.0, zero_double = 0.0; + const std::complex one_complex = {1.0, 0.0}; + const std::complex zero_complex = {0.0, 0.0}; + const char N_char = 'N'; + const char C_char = 'C'; // Using 'C' is consistent with the formula + +#ifdef __MPI + const int nbasis = ParaV->desc[2]; + const int nbands = ParaV->desc_wfc[3]; + + //because wfc(bands, basis'), H(basis, basis'), we do wfc*H^T(in the perspective of cpp, not in fortran). And get H_wfc(bands, basis) is correct. + pzgemm_( &C_char, &N_char, &nbasis, &nbands, &nbasis, &one_complex, &HK, &one_int, &one_int, ParaV->desc, + &wfc, &one_int, &one_int, ParaV->desc_wfc, &zero_complex, &H_wfc, &one_int, &one_int, ParaV->desc_wfc ); +#endif +} + + +template <> +void HkPsi(const Parallel_Orbitals* ParaV, const double& HK, const double& wfc, double& H_wfc); + + +template +void psiDotPsi(const Parallel_Orbitals* ParaV, const Parallel_2D& para_Eij_in, const TK& wfc, const TK& H_wfc, std::vector& Dmn, double* wfcHwfc) +{ + const int one_int = 1; + const std::complex one_complex = {1.0, 0.0}; + const std::complex zero_complex = {0.0, 0.0}; + const char N_char = 'N'; + const char C_char = 'C'; + + const int nrow_bands = para_Eij_in.get_row_size(); + const int ncol_bands = para_Eij_in.get_col_size(); + +#ifdef __MPI + const int nbasis = ParaV->desc[2]; + const int nbands = ParaV->desc_wfc[3]; + + pzgemm_( &C_char, &N_char, &nbands, &nbands, &nbasis, &one_complex, &wfc, &one_int, &one_int, ParaV->desc_wfc, + &H_wfc, &one_int, &one_int, ParaV->desc_wfc, &zero_complex, &Dmn[0], &one_int, &one_int, para_Eij_in.desc ); +#endif + + for(int i=0; i +void psiDotPsi(const Parallel_Orbitals* ParaV, const Parallel_2D& para_wfc_in, + const double& wfc, const double& H_wfc, std::vector& Dmn, double* wfcHwfc); + + +// realize occNum_wfc = occNum * wfc. Calling this function and we can get wfc = occNum*wfc. +template +void occNum_MulPsi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& occ_number, psi::Psi& wfc, int symbol = 0, + const std::string XC_func_rdmft = "hf", const double alpha = 1.0) +{ + const int nk_local = wfc.get_nk(); + const int nbands_local = wfc.get_nbands(); + const int nbasis_local = wfc.get_nbasis(); + + // const int nbasis = ParaV->desc[2]; // need to be deleted + // const int nbands = ParaV->desc_wfc[3]; + + for (int ik = 0; ik < nk_local; ++ik) + { + for (int ib_local = 0; ib_local < nbands_local; ++ib_local) // ib_local < nbands_local , some problem, ParaV->ncol_bands + { + const double occNum_local = occNum_func( occ_number(ik, ParaV->local2global_col(ib_local)), symbol, XC_func_rdmft, alpha); + TK* wfc_pointer = &(wfc(ik, ib_local, 0)); + BlasConnector::scal(nbasis_local, occNum_local, wfc_pointer, 1); + } + } +} + + +// add psi with eta and g(eta) +template +void add_psi(const Parallel_Orbitals* ParaV, + const K_Vectors* kv, + const ModuleBase::matrix& occ_number, + psi::Psi& psi_TV, + psi::Psi& psi_hartree, + psi::Psi& psi_dft_XC, + psi::Psi& psi_exx_XC, + psi::Psi& occNum_Hpsi, + const std::string XC_func_rdmft = "hf", + const double alpha = 1.0) +{ + const int nk = psi_TV.get_nk(); + const int nbn_local = psi_TV.get_nbands(); + const int nbs_local = psi_TV.get_nbasis(); + occNum_MulPsi(ParaV, occ_number, psi_TV); + occNum_MulPsi(ParaV, occ_number, psi_hartree); + occNum_MulPsi(ParaV, occ_number, psi_dft_XC); + occNum_MulPsi(ParaV, occ_number, psi_exx_XC, 2, XC_func_rdmft, alpha); + + // const int nbasis = ParaV->desc[2]; + // const int nbands = ParaV->desc_wfc[3]; + + for(int ik=0; ikwk[ik], p_occNum_Hpsi, 1); + } + } + +} + + +// occNum_wfcHwfc = occNum*wfcHwfc + occNum_wfcHwfc +// When symbol = 0, 1, 2, 3, 4, occNum = occNum, 0.5*occNum, g(occNum), 0.5*g(occNum), d_g(occNum)/d_occNum respectively. Default symbol=0. +void occNum_Mul_wfcHwfc(const ModuleBase::matrix& occ_number, + const ModuleBase::matrix& wfcHwfc, + ModuleBase::matrix& occNum_wfcHwfc, + int symbol = 0, + const std::string XC_func_rdmft = "hf", + const double alpha = 1.0); + + +// Default symbol = 0 for the gradient of Etotal with respect to occupancy +// symbol = 1 for the relevant calculation of Etotal +void add_occNum(const K_Vectors& kv, + const ModuleBase::matrix& occ_number, + const ModuleBase::matrix& wfcHwfc_TV_in, + const ModuleBase::matrix& wfcHwfc_hartree_in, + const ModuleBase::matrix& wfcHwfc_dft_XC_in, + const ModuleBase::matrix& wfcHwfc_exx_XC_in, + ModuleBase::matrix& occNum_wfcHwfc, + const std::string XC_func_rdmft = "hf", + const double alpha = 1.0); + + +// // do wk*g(occNum)*wfcHwfc and add for TV, hartree, XC. This function just use once, so it can be replace and delete +// void add_wfcHwfc(const std::vector& wk_in, const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in, +// const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha); +void add_wfcHwfc(const ModuleBase::matrix& wg, + const ModuleBase::matrix& wk_fun_occNum, + const ModuleBase::matrix& wfcHwfc_TV_in, + const ModuleBase::matrix& wfcHwfc_hartree_in, + const ModuleBase::matrix& wfcHwfc_XC_in, + ModuleBase::matrix& occNum_wfcHwfc, + const std::string XC_func_rdmft, + const double alpha); + + +//give certain occNum_wfcHwfc, get the corresponding energy +double getEnergy(const ModuleBase::matrix& occNum_wfcHwfc); + + + + + + +// this part of the code is copying from class Veff and do some modifications. +template +class Veff_rdmft : public hamilt::OperatorLCAO +{ + public: + /** + * @brief Construct a new Veff object for multi-kpoint calculation + * @param GK_in: the pointer of Gint_k object, used for grid integration + */ + Veff_rdmft(Gint_k* GK_in, + hamilt::HS_Matrix_K* hsk_in, + const std::vector>& kvec_d_in, + elecstate::Potential* pot_in, + hamilt::HContainer* hR_in, + const UnitCell* ucell_in, + const std::vector& orb_cutoff, + Grid_Driver* GridD_in, + const int& nspin, + + const Charge* charge_in, + const ModulePW::PW_Basis* rho_basis_in, + const ModuleBase::matrix* vloc_in, + const ModuleBase::ComplexMatrix* sf_in, + const std::string potential_in, + double* etxc_in = nullptr, + double* vtxc_in = nullptr + ) + : GK(GK_in), + orb_cutoff_(orb_cutoff), + pot(pot_in), + ucell(ucell_in), + gd(GridD_in), + hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), + + charge_(charge_in), + rho_basis_(rho_basis_in), + vloc_(vloc_in), + sf_(sf_in), + potential_(potential_in), + etxc(etxc_in), + vtxc(vtxc_in) + { + this->cal_type = hamilt::calculation_type::lcao_gint; + + this->initialize_HR(ucell_in, GridD_in); + + GK_in->initialize_pvpR(*ucell_in, GridD_in, nspin); + } + Veff_rdmft(Gint_Gamma* GG_in, + hamilt::HS_Matrix_K* hsk_in, + const std::vector>& kvec_d_in, + elecstate::Potential* pot_in, + hamilt::HContainer* hR_in, + const UnitCell* ucell_in, + const std::vector& orb_cutoff, + Grid_Driver* GridD_in, + const int& nspin, + + const Charge* charge_in, + const ModulePW::PW_Basis* rho_basis_in, + const ModuleBase::matrix* vloc_in, + const ModuleBase::ComplexMatrix* sf_in, + const std::string potential_in, + double* etxc_in = nullptr, + double* vtxc_in = nullptr + ) + : GG(GG_in), + orb_cutoff_(orb_cutoff), + pot(pot_in), + hamilt::OperatorLCAO(hsk_in, kvec_d_in, hR_in), + + ucell(ucell_in), + gd(GridD_in), + charge_(charge_in), + rho_basis_(rho_basis_in), + vloc_(vloc_in), + sf_(sf_in), + potential_(potential_in), + etxc(etxc_in), + vtxc(vtxc_in) + { + this->cal_type = hamilt::calculation_type::lcao_gint; + + this->initialize_HR(ucell_in, GridD_in); + + GG_in->initialize_pvpR(*ucell_in, GridD_in, nspin); + } + + ~Veff_rdmft(){}; + + /** + * @brief contributeHR() is used to calculate the HR matrix + * + * the contribution of V_{eff} is calculated by the contribution of V_{H} and V_{XC} and V_{local pseudopotential} and so on. + * grid integration is used to calculate the contribution Hamiltonian of effective potential + */ + virtual void contributeHR() override; + + const UnitCell* ucell; + + Grid_Driver* gd; + + private: + // used for k-dependent grid integration. + Gint_k* GK = nullptr; + + // used for gamma only algorithms. + Gint_Gamma* GG = nullptr; + + std::vector orb_cutoff_; + + // Charge calculating method in LCAO base and contained grid base calculation: DM_R, DM, pvpR_reduced + + elecstate::Potential* pot = nullptr; + + int nspin = 1; + int current_spin = 0; + + /** + * @brief initialize HR, search the nearest neighbor atoms + * HContainer is used to store the electronic kinetic matrix with specific atom-pairs + * the size of HR will be fixed after initialization + */ + void initialize_HR(const UnitCell* ucell_in, Grid_Driver* GridD_in); + + + // added by jghan + + const Charge* charge_; + + std::string potential_; + + const ModulePW::PW_Basis* rho_basis_; + + const ModuleBase::matrix* vloc_; + + const ModuleBase::ComplexMatrix* sf_; + + double* etxc; + + double* vtxc; + +}; + + + + + + + +} + +#endif diff --git a/source/module_rdmft/update_state_rdmft.cpp b/source/module_rdmft/update_state_rdmft.cpp new file mode 100644 index 0000000000..24d6294dd4 --- /dev/null +++ b/source/module_rdmft/update_state_rdmft.cpp @@ -0,0 +1,214 @@ +//========================================================== +// Author: Jingang Han +// DATE : 2024-10-30 +//========================================================== + +#include "rdmft.h" +#include "module_rdmft/rdmft_tools.h" +#include "module_elecstate/module_dm/cal_dm_psi.h" +#include "module_elecstate/module_dm/density_matrix.h" +#include "module_elecstate/module_charge/symmetry_rho.h" + + +namespace rdmft +{ + + + +template +void RDMFT::update_ion(UnitCell& ucell_in, ModulePW::PW_Basis& rho_basis_in, + ModuleBase::matrix& vloc_in, ModuleBase::ComplexMatrix& sf_in) +{ + ucell = &ucell_in; + rho_basis = &rho_basis_in; + vloc = &vloc_in; + sf = &sf_in; + + HR_TV->set_zero(); + this->cal_V_TV(); +#ifdef __EXX + if( GlobalC::exx_info.info_global.cal_exx ) + { + if (GlobalC::exx_info.info_ri.real_number) + { + Vxc_fromRI_d->cal_exx_ions(); + } + else + { + Vxc_fromRI_c->cal_exx_ions(); + } + } +#endif + + std::cout << "\n\n\n******\ndo rdmft_esolver.update_ion() successfully\n******\n\n\n" << std::endl; +} + + +template +void RDMFT::update_elec(const ModuleBase::matrix& occ_number_in, const psi::Psi& wfc_in, const Charge* charge_in) +{ + // update occ_number, wg, wk_fun_occNum + occ_number = (occ_number_in); + wg = (occ_number); + + for(int ik=0; ik < wg.nr; ++ik) + { + for(int inb=0; inb < wg.nc; ++inb) + { + wg(ik, inb) *= kv->wk[ik]; + wk_fun_occNum(ik, inb) = kv->wk[ik] * occNum_func(occ_number(ik, inb), 2, XC_func_rdmft, alpha_power); + } + } + + // update wfc + TK* pwfc_in = &wfc_in(0, 0, 0); + TK* pwfc = &wfc(0, 0, 0); + for(int i=0; iupdate_charge(); + + // "default" = "pbe" + // if( !only_exx_type || this->cal_E_type != 1 ) + if( this->cal_E_type != 1 ) + { + // the second cal_E_type need the complete pot to get effctive_V to calEband and so on. + this->pelec->pot->update_from_charge(charge, ucell); + } + + this->cal_V_hartree(); + this->cal_V_XC(); + // this->cal_Hk_Hpsi(); + + std::cout << "\n******\n" << "update elec in rdmft successfully" << "\n******\n" << std::endl; +} + + +// this code is copying from function ElecStateLCAO::psiToRho(), in elecstate_lcao.cpp +template +void RDMFT::update_charge() +{ + if( PARAM.inp.gamma_only ) + { + // calculate DMK and DMR + elecstate::DensityMatrix DM_gamma_only(ParaV, nspin); + elecstate::cal_dm_psi(ParaV, wg, wfc, DM_gamma_only); + DM_gamma_only.init_DMR(&GlobalC::GridD, &GlobalC::ucell); + DM_gamma_only.cal_DMR(); + + for (int is = 0; is < nspin; is++) + { + ModuleBase::GlobalFunc::ZEROS(charge->rho[is], charge->nrxx); + } + + GG->transfer_DM2DtoGrid(DM_gamma_only.get_DMR_vector()); + Gint_inout inout(charge->rho, Gint_Tools::job_type::rho, nspin); + GG->cal_gint(&inout); + + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { + // for (int is = 0; is < nspin; is++) + // { + // ModuleBase::GlobalFunc::ZEROS(charge->kin_r[is], charge->nrxx); + // } + // Gint_inout inout1(charge->kin_r, Gint_Tools::job_type::tau); + // GG->cal_gint(&inout1); + this->pelec->cal_tau(wfc); + } + + charge->renormalize_rho(); + } + else + { + // calculate DMK and DMR + elecstate::DensityMatrix DM(ParaV, nspin, kv->kvec_d, nk_total); + elecstate::cal_dm_psi(ParaV, wg, wfc, DM); + DM.init_DMR(&GlobalC::GridD, &GlobalC::ucell); + DM.cal_DMR(); + + for (int is = 0; is < nspin; is++) + { + ModuleBase::GlobalFunc::ZEROS(charge->rho[is], charge->nrxx); + } + + GK->transfer_DM2DtoGrid(DM.get_DMR_vector()); + Gint_inout inout(charge->rho, Gint_Tools::job_type::rho, nspin); + GK->cal_gint(&inout); + + if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) + { + // for (int is = 0; is < nspin; is++) + // { + // ModuleBase::GlobalFunc::ZEROS(charge->kin_r[is], charge->nrxx); + // } + // Gint_inout inout1(charge->kin_r, Gint_Tools::job_type::tau); + // GK->cal_gint(&inout1); + this->pelec->cal_tau(wfc); + } + + charge->renormalize_rho(); + } + + // charge density symmetrization + // this->pelec->calculate_weights(); + // this->pelec->calEBand(); + Symmetry_rho srho; + for (int is = 0; is < nspin; is++) + { + srho.begin(is, *(this->charge), rho_basis, GlobalC::ucell.symm); + } + + // what this? it seems that it needs to be updated at each iteration + if (PARAM.inp.vl_in_h) + { + // update Gint_K + if (!PARAM.globalv.gamma_only_local) + { + this->GK->renew(); + } + } + +} + + +template +void RDMFT::update_occNumber(const ModuleBase::matrix& occ_number_in) +{ + occ_number = (occ_number_in); + wg = (occ_number); + for(int ik=0; ik < wg.nr; ++ik) + { + for(int inb=0; inb < wg.nc; ++inb) + { + wg(ik, inb) *= kv->wk[ik]; + wk_fun_occNum(ik, inb) = kv->wk[ik] * occNum_func(occ_number(ik, inb), 2, XC_func_rdmft, alpha_power); + } + } +} + + +template +void RDMFT::update_wg(const ModuleBase::matrix& wg_in) +{ + wg = (wg_in); + occ_number = (wg); + for(int ik=0; ik < wg.nr; ++ik) + { + for(int inb=0; inb < wg.nc; ++inb) + { + occ_number(ik, inb) /= kv->wk[ik]; + wk_fun_occNum(ik, inb) = kv->wk[ik] * occNum_func(occ_number(ik, inb), 2, XC_func_rdmft, alpha_power); + } + } +} + + +template class RDMFT; +template class RDMFT, double>; +template class RDMFT, std::complex>; + +} + + + diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h index bb04a3f358..6c57a1816c 100644 --- a/source/module_ri/Exx_LRI.h +++ b/source/module_ri/Exx_LRI.h @@ -59,6 +59,10 @@ class Exx_LRI void init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in, const LCAO_Orbitals& orb); void cal_exx_force(); void cal_exx_stress(); + void cal_exx_ions(const bool write_cv = false); + void cal_exx_elec(const std::vector>>>& Ds, + const Parallel_Orbitals& pv, + const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr); std::vector> get_abfs_nchis() const; std::vector< std::map>>> Hexxs; @@ -71,7 +75,7 @@ class Exx_LRI const Exx_Info::Exx_Info_RI &info; MPI_Comm mpi_comm; const K_Vectors *p_kv = nullptr; - std::vector orb_cutoff_; + std::vector orb_cutoff_; std::vector>> lcaos; std::vector>> abfs; @@ -80,10 +84,6 @@ class Exx_LRI LRI_CV cv; RI::Exx exx_lri; - void cal_exx_ions(const bool write_cv = false); - void cal_exx_elec(const std::vector>>>& Ds, - const Parallel_Orbitals& pv, - const ModuleSymmetry::Symmetry_rotation* p_symrot = nullptr); void post_process_Hexx( std::map>> &Hexxs_io ) const; double post_process_Eexx(const double& Eexx_in) const; @@ -99,4 +99,4 @@ class Exx_LRI #include "Exx_LRI.hpp" -#endif +#endif diff --git a/source/module_ri/Exx_LRI_interface.hpp b/source/module_ri/Exx_LRI_interface.hpp index 5f7ba2b964..6b9b41e6a3 100644 --- a/source/module_ri/Exx_LRI_interface.hpp +++ b/source/module_ri/Exx_LRI_interface.hpp @@ -72,6 +72,12 @@ void Exx_LRI_Interface::exx_beforescf(const int istep, const K_Vectors { XC_Functional::set_xc_type("scan"); } + // added by jghan, 2024-07-07 + else if ( ucell.atoms[0].ncpp.xc_func == "MULLER" || ucell.atoms[0].ncpp.xc_func == "POWER" + || ucell.atoms[0].ncpp.xc_func == "WP22" || ucell.atoms[0].ncpp.xc_func == "CWP22" ) + { + XC_Functional::set_xc_type("pbe"); + } } this->exx_ptr->cal_exx_ions(PARAM.inp.out_ri_cv); } @@ -122,10 +128,11 @@ void Exx_LRI_Interface::exx_eachiterinit(const int istep, const elecst if (this->exx_spacegroup_symmetry && GlobalC::exx_info.info_global.exx_symmetry_realspace) { this->exx_ptr->cal_exx_elec(Ds, *dm_in.get_paraV_pointer(), &this->symrot_); } else { this->exx_ptr->cal_exx_elec(Ds, *dm_in.get_paraV_pointer()); } }; - if(istep > 0 && flag_restart) + if(istep > 0 && flag_restart) { cal(*dm_last_step); - else + } else { cal(dm); +} } } } diff --git a/source/module_ri/conv_coulomb_pot_k.cpp b/source/module_ri/conv_coulomb_pot_k.cpp index 14870a4f72..6394f9148d 100644 --- a/source/module_ri/conv_coulomb_pot_k.cpp +++ b/source/module_ri/conv_coulomb_pot_k.cpp @@ -11,8 +11,9 @@ namespace Conv_Coulomb_Pot_K const std::vector & psif) { std::vector psik2_ccp(psif.size()); - for( size_t ik=0; ik psik2_ccp(psif.size()); - for (size_t ik = 0; ik < psif.size(); ++ik) + for (size_t ik = 0; ik < psif.size(); ++ik) { psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1 - std::cos(k_radial[ik] * hf_Rcut)); +} return psik2_ccp; } @@ -36,11 +38,27 @@ namespace Conv_Coulomb_Pot_K const double hse_omega) { std::vector psik2_ccp(psif.size()); - for( size_t ik=0; ik cal_psi_erf( + const std::vector & psif, + const std::vector & k_radial, + const double hse_omega, + const double hf_Rcut) + { + std::vector psik2_ccp(psif.size()); + for( size_t ik=0; ik @@ -59,6 +77,8 @@ namespace Conv_Coulomb_Pot_K psik2_ccp = cal_psi_hf( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hf_Rcut")); break; case Ccp_Type::Hse: psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") ); break; + case Ccp_Type::erf: + psik2_ccp = cal_psi_erf( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega"), parameter.at("hf_Rcut") ); break; default: throw( ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) ); break; } @@ -66,15 +86,19 @@ namespace Conv_Coulomb_Pot_K const double dr = orbs.get_rab().back(); const int Nr = (static_cast(orbs.getNr()*rmesh_times)) | 1; std::vector rab(Nr); - for( size_t ir=0; ir r_radial(Nr); - for( size_t ir=0; ir=0; --ir) { - if(std::abs(orbs.getPsi(ir))>=psi_threshold) + if(std::abs(orbs.getPsi(ir))>=psi_threshold) { return static_cast(ir)/orbs.getNr(); +} } return 0.0; } diff --git a/source/module_ri/conv_coulomb_pot_k.h b/source/module_ri/conv_coulomb_pot_k.h index d464a53f91..269691606c 100644 --- a/source/module_ri/conv_coulomb_pot_k.h +++ b/source/module_ri/conv_coulomb_pot_k.h @@ -10,7 +10,8 @@ namespace Conv_Coulomb_Pot_K enum class Ccp_Type{ // parameter: Ccp, // Hf, // "hf_Rcut" - Hse}; // "hse_omega" + Hse, // "hse_omega" + erf}; // "hse_omega" template T cal_orbs_ccp( const T &orbs, @@ -34,6 +35,12 @@ namespace Conv_Coulomb_Pot_K const std::vector & psif, const std::vector & k_radial, const double hse_omega); + // added by jghan, 2024-07-06 + std::vector cal_psi_erf( + const std::vector & psif, + const std::vector & k_radial, + const double hse_omega, + const double hf_Rcut); } #include "conv_coulomb_pot_k.hpp" diff --git a/tests/integrate/1001_NO_Si2_dzp_rdmft/INPUT b/tests/integrate/1001_NO_Si2_dzp_rdmft/INPUT new file mode 100644 index 0000000000..43e648553b --- /dev/null +++ b/tests/integrate/1001_NO_Si2_dzp_rdmft/INPUT @@ -0,0 +1,32 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation scf + +#nbands 8 +symmetry 1 + +#Parameters (2.Iteration) +ecutwfc 60 +scf_thr 1e-6 +scf_nmax 100 +cal_force 0 +cal_stress 0 +#Parameters (3.Basis) +basis_type lcao + +#Parameters (4.Smearing) +smearing_method gauss +smearing_sigma 0.002 + +#Parameters (5.Mixing) +mixing_type broyden +#mixing_beta 0.3 +#gamma_only 1 +ks_solver genelpa + +#Parameters (6.rdmft) +rdmft 1 +dft_functional cwp22 +exx_hse_omega 0.33 +rdmft_power_alpha 0.656 \ No newline at end of file diff --git a/tests/integrate/1001_NO_Si2_dzp_rdmft/KPT b/tests/integrate/1001_NO_Si2_dzp_rdmft/KPT new file mode 100644 index 0000000000..4fd38968a0 --- /dev/null +++ b/tests/integrate/1001_NO_Si2_dzp_rdmft/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 2 1 0 0 0 diff --git a/tests/integrate/1001_NO_Si2_dzp_rdmft/STRU b/tests/integrate/1001_NO_Si2_dzp_rdmft/STRU new file mode 100644 index 0000000000..cb2b7ddd38 --- /dev/null +++ b/tests/integrate/1001_NO_Si2_dzp_rdmft/STRU @@ -0,0 +1,21 @@ +ATOMIC_SPECIES +Si 14.000 ../../PP_ORB/Si_ONCV_PBE_FR-1.1.upf + +NUMERICAL_ORBITAL +../../PP_ORB/Si_gga_6au_100Ry_2s2p1d.orb + +LATTICE_CONSTANT +10.2 + +LATTICE_VECTORS +0.5 0.5 0.0 #Lattice vector 1 +0.5 0.0 0.5 #Lattice vector 2 +0.0 0.5 0.5 #Lattice vector 3 + +ATOMIC_POSITIONS +Cartesian #Cartesian(Unit is LATTICE_CONSTANT) +Si #Name of element +0.0 #Magnetic for this element. +2 #Number of atoms +0.00 0.00 0.00 0 0 0 #x,y,z, move_x, move_y, move_z +0.25 0.25 0.25 1 1 1 \ No newline at end of file diff --git a/tests/integrate/1001_NO_Si2_dzp_rdmft/jd b/tests/integrate/1001_NO_Si2_dzp_rdmft/jd new file mode 100644 index 0000000000..994b0a8f3a --- /dev/null +++ b/tests/integrate/1001_NO_Si2_dzp_rdmft/jd @@ -0,0 +1 @@ +Test the various energies obtained when using rdmft to calculate solid Si, KPT is 221 \ No newline at end of file diff --git a/tests/integrate/1001_NO_Si2_dzp_rdmft/result.ref b/tests/integrate/1001_NO_Si2_dzp_rdmft/result.ref new file mode 100644 index 0000000000..b30d6f3c85 --- /dev/null +++ b/tests/integrate/1001_NO_Si2_dzp_rdmft/result.ref @@ -0,0 +1,18 @@ +etotref -201.4138882681562848 +etotperatomref -100.7069441341 +pointgroupref BvK +spacegroupref BvK +nksibzref 3 + +The following energy units are in Rydberg: +E_TV_RDMFT_ref 5.2788851446 +E_hartree_RDMFT_ref 1.3622516604 +Exc_cwp22_RDMFT_ref -4.5448173113 +E_Ewald_ref -16.8997585734 +E_entropy_ref -0.0000000000 +E_descf_ref 0.0000000000 +Etotal_RDMFT_ref -14.8034390797 +Exc_ksdft_ref -2.4788985117 +E_exx_ksdft_ref -2.0657389944 + +totaltimeref 118.49 diff --git a/tests/integrate/CASES_CPU.txt b/tests/integrate/CASES_CPU.txt index 822bc4b126..1881752d3a 100644 --- a/tests/integrate/CASES_CPU.txt +++ b/tests/integrate/CASES_CPU.txt @@ -355,3 +355,4 @@ 920_OF_MD_LibxcPBE 921_OF_CR_LDA 922_NO_FeBiTe +1001_NO_Si2_dzp_rdmft diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index fe3b1f703a..f42eec827c 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -69,6 +69,7 @@ has_mat_dh=$(get_input_key_value "out_mat_dh" "INPUT") has_scan=$(get_input_key_value "dft_functional" "INPUT") out_chg=$(get_input_key_value "out_chg" "INPUT") esolver_type=$(get_input_key_value "esolver_type" "INPUT") +rdmft=$(get_input_key_value "rdmft" "INPUT") #echo $running_path base=$(get_input_key_value "basis_type" "INPUT") word="driver_line" @@ -543,6 +544,40 @@ if [ $is_lr == 1 ]; then echo "totexcitationenergyref $lreig_tot" >>$1 fi +if ! test -z "$rdmft" && [[ $rdmft == 1 ]]; then + echo "" >>$1 + echo "The following energy units are in Rydberg:" >>$1 + + E_TV_RDMFT=$(grep "E_TV_RDMFT" "$running_path" | tail -1 | awk '{print $2}') + echo "E_TV_RDMFT_ref $E_TV_RDMFT" >>$1 + + E_hartree_RDMFT=$(grep "E_hartree_RDMFT" "$running_path" | tail -1 | awk '{print $2}') + echo "E_hartree_RDMFT_ref $E_hartree_RDMFT" >>$1 + + Exc_cwp22_RDMFT=$(grep "Exc_cwp22_RDMFT" "$running_path" | tail -1 | awk '{print $2}') + echo "Exc_cwp22_RDMFT_ref $Exc_cwp22_RDMFT" >>$1 + + E_Ewald=$(grep "E_Ewald" "$running_path" | tail -1 | awk '{print $2}') + echo "E_Ewald_ref $E_Ewald" >>$1 + + E_entropy=$(grep "E_entropy(-TS)" "$running_path" | tail -1 | awk '{print $2}') + echo "E_entropy_ref $E_entropy" >>$1 + + E_descf=$(grep "E_descf" "$running_path" | tail -1 | awk '{print $2}') + echo "E_descf_ref $E_descf" >>$1 + + Etotal_RDMFT=$(grep "Etotal_RDMFT" "$running_path" | tail -1 | awk '{print $2}') + echo "Etotal_RDMFT_ref $Etotal_RDMFT" >>$1 + + Exc_ksdft=$(grep "Exc_ksdft" "$running_path" | tail -1 | awk '{print $2}') + echo "Exc_ksdft_ref $Exc_ksdft" >>$1 + + E_exx_ksdft=$(grep "E_exx_ksdft" "$running_path" | tail -1 | awk '{print $2}') + echo "E_exx_ksdft_ref $E_exx_ksdft" >>$1 + + echo "" >>$1 +fi + #echo $total_band ttot=`grep $word $running_path | awk '{print $3}'` echo "totaltimeref $ttot" >>$1