From 913d2e3db06f125899dfbc4dbf7589b999e44a35 Mon Sep 17 00:00:00 2001 From: linpz Date: Mon, 9 Sep 2024 20:09:27 +0800 Subject: [PATCH] Fix bugs in XC_Functional_Libxc --- .../module_xc/xc_functional.cpp | 7 +- .../module_xc/xc_functional.h | 3 +- .../module_xc/xc_functional_libxc.cpp | 4 + .../module_xc/xc_functional_libxc_tools.cpp | 292 ++++++++++++++++++ .../module_xc/xc_functional_libxc_vxc.cpp | 24 +- .../xc_functional_libxc_wrapper_tauxc.cpp | 2 +- source/module_io/input_conv.cpp | 2 +- 7 files changed, 318 insertions(+), 16 deletions(-) create mode 100644 source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp index 0168ecce99..c9730e957b 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional.cpp @@ -18,11 +18,16 @@ int XC_Functional::func_type = 0; bool XC_Functional::use_libxc = true; double XC_Functional::hybrid_alpha = 0.25; -void XC_Functional::get_hybrid_alpha(const double alpha_in) +void XC_Functional::set_hybrid_alpha(const double alpha_in) { hybrid_alpha = alpha_in; } +double XC_Functional::get_hybrid_alpha() +{ + return hybrid_alpha; +} + int XC_Functional::get_func_type() { return func_type; diff --git a/source/module_hamilt_general/module_xc/xc_functional.h b/source/module_hamilt_general/module_xc/xc_functional.h index 5179c6bf93..ac9c75256a 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.h +++ b/source/module_hamilt_general/module_xc/xc_functional.h @@ -66,7 +66,8 @@ class XC_Functional static void set_xc_type(const std::string xc_func_in); // For hybrid functional - static void get_hybrid_alpha(const double alpha_in); + static void set_hybrid_alpha(const double alpha_in); + static double get_hybrid_alpha(); /// Usually in exx caculation, the first SCF loop should be converged with PBE static void set_xc_first_loop(const UnitCell& ucell); 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 416b08b65e..8d80f27714 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp @@ -4,6 +4,10 @@ #include "module_parameter/parameter.h" #include "module_base/tool_quit.h" +#ifdef __EXX +#include "module_hamilt_pw/hamilt_pwdft/global.h" // just for GlobalC::exx_info +#endif + #include #include diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp new file mode 100644 index 0000000000..2ff901e786 --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp @@ -0,0 +1,292 @@ +#ifdef USE_LIBXC + +#include "xc_functional_libxc.h" +#include "xc_functional.h" +#include "module_elecstate/module_charge/charge.h" + +// converting rho (abacus=>libxc) +std::vector XC_Functional_Libxc::convert_rho( + const int nspin, + const std::size_t nrxx, + const Charge* const chr) +{ + std::vector rho(nrxx*nspin); + #ifdef _OPENMP + #pragma omp parallel for collapse(2) schedule(static, 1024) + #endif + for( int is=0; isrho[is][ir] + 1.0/nspin*chr->rho_core[ir]; + return rho; +} + +// converting rho (abacus=>libxc) +std::tuple, std::vector> +XC_Functional_Libxc::convert_rho_amag_nspin4( + const int nspin, + const std::size_t nrxx, + const Charge* const chr) +{ + assert(GlobalV::NSPIN==4); + std::vector rho(nrxx*nspin); + std::vector amag(nrxx); + #ifdef _OPENMP + #pragma omp parallel for + #endif + for( int ir=0; irrho[0][ir] + chr->rho_core[ir] ); + amag[ir] = std::sqrt( std::pow(chr->rho[1][ir],2) + + std::pow(chr->rho[2][ir],2) + + std::pow(chr->rho[3][ir],2) ); + const double amag_clip = (amag[ir]>> +XC_Functional_Libxc::cal_gdr( + const int nspin, + const std::vector &rho, + const double tpiba, + const Charge* const chr) +{ + std::vector>> gdr(nspin); + const std::size_t nrxx = rho.size(); + for( int is=0; is!=nspin; ++is ) + { + std::vector rhor(nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for(std::size_t ir=0; ir> rhog(chr->rhopw->npw); + chr->rhopw->real2recip(rhor.data(), rhog.data()); + + //------------------------------------------- + // compute the gradient of charge density and + // store the gradient in gdr[is] + //------------------------------------------- + gdr[is].resize(nrxx); + XC_Functional::grad_rho(rhog.data(), gdr[is].data(), chr->rhopw, tpiba); + } // end for(is) + return gdr; +} + +// converting grho (abacus=>libxc) +std::vector XC_Functional_Libxc::convert_sigma( + const std::vector>> &gdr) +{ + const std::size_t nspin = gdr.size(); + assert(nspin>0); + const std::size_t nrxx = gdr[0].size(); + for(std::size_t is=1; is sigma( nrxx * ((1==nspin)?1:3) ); + if( 1==nspin ) + { + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for( std::size_t ir=0; ir XC_Functional_Libxc::cal_sgn( + const double rho_threshold, + const double grho_threshold, + const xc_func_type &func, + const int nspin, + const std::size_t nrxx, + const std::vector &rho, + const std::vector &sigma) +{ + std::vector sgn(nrxx*nspin, 1.0); + // in the case of GGA correlation for polarized case, + // a cutoff for grho is required to ensure that libxc gives reasonable results + if(nspin==2 && func.info->family != XC_FAMILY_LDA && func.info->kind==XC_CORRELATION) + { + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 512) + #endif + for( int ir=0; irabacus) +double XC_Functional_Libxc::convert_etxc( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector &rho, + std::vector exc) +{ + double etxc = 0.0; + #ifdef _OPENMP + #pragma omp parallel for collapse(2) reduction(+:etxc) schedule(static, 256) + #endif + for( int is=0; isabacus) +double XC_Functional_Libxc::convert_vtxc_v( + const xc_func_type &func, + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector &rho, + const std::vector>> &gdr, + const std::vector &vrho, + const std::vector &vsigma, + const double tpiba, + const Charge* const chr, + ModuleBase::matrix &v) +{ + assert(v.nr==nspin); + assert(v.nc==nrxx); + + double vtxc = 0.0; + + #ifdef _OPENMP + #pragma omp parallel for collapse(2) reduction(+:vtxc) schedule(static, 256) + #endif + for( int is=0; isfamily == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) + { + const std::vector> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr); + + double rvtxc = 0.0; + #ifdef _OPENMP + #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) + #endif + for( int is=0; isfamily == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)) + + return vtxc; +} + + +// dh for gga v +std::vector> XC_Functional_Libxc::cal_dh( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector>> &gdr, + const std::vector &vsigma, + const double tpiba, + const Charge* const chr) +{ + std::vector>> h( + nspin, + std::vector>(nrxx) ); + if( 1==nspin ) + { + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for( std::size_t ir=0; ir> dh(nspin, std::vector(nrxx)); + for( int is=0; is!=nspin; ++is ) + XC_Functional::grad_dot( h[is].data(), dh[is].data(), chr->rhopw, tpiba); + + return dh; +} + + +// convert v for NSPIN=4 +ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4( + const std::size_t nrxx, + const Charge* const chr, + const std::vector &amag, + const ModuleBase::matrix &v) +{ + assert(GlobalV::NSPIN==4); + constexpr double vanishing_charge = 1.0e-10; + ModuleBase::matrix v_nspin4(GlobalV::NSPIN, nrxx); + for( int ir=0; ir vanishing_charge ) + { + const double vs = 0.5 * (v(0,ir)-v(1,ir)); + for(int ipol=1; ipolrho[ipol][ir] / amag[ir]; + } + } + } + return v_nspin4; +} + +#endif \ No newline at end of file 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 9e8614dfa6..e8b80537f1 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 @@ -260,9 +260,9 @@ std::tuple XC_Functional_Li for( int ir=0; ir< nrxx; ++ir ) { #ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) + if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { - exc[ir] *= (1.0 - XC_Functional::hybrid_alpha); + exc[ir] *= (1.0 - XC_Functional::get_hybrid_alpha()); } #endif etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is]; @@ -278,9 +278,9 @@ std::tuple XC_Functional_Li for( int ir=0; ir< nrxx; ++ir ) { #ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) + if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { - vrho[ir*nspin+is] *= (1.0 - XC_Functional::hybrid_alpha); + vrho[ir*nspin+is] *= (1.0 - XC_Functional::get_hybrid_alpha()); } #endif const double v_tmp = ModuleBase::e2 * vrho[ir*nspin+is] * sgn[ir*nspin+is]; @@ -301,9 +301,9 @@ std::tuple XC_Functional_Li for( int ir=0; ir< nrxx; ++ir ) { #ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) + if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { - vsigma[ir] *= (1.0 - XC_Functional::hybrid_alpha); + vsigma[ir] *= (1.0 - XC_Functional::get_hybrid_alpha()); } #endif h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir]; @@ -317,11 +317,11 @@ std::tuple XC_Functional_Li for( int ir=0; ir< nrxx; ++ir ) { #ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) + if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { - vsigma[ir*3] *= (1.0 - XC_Functional::hybrid_alpha); - vsigma[ir*3+1] *= (1.0 - XC_Functional::hybrid_alpha); - vsigma[ir*3+2] *= (1.0 - XC_Functional::hybrid_alpha); + vsigma[ir*3] *= (1.0 - XC_Functional::get_hybrid_alpha()); + vsigma[ir*3+1] *= (1.0 - XC_Functional::get_hybrid_alpha()); + vsigma[ir*3+2] *= (1.0 - XC_Functional::get_hybrid_alpha()); } #endif h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 @@ -363,9 +363,9 @@ std::tuple XC_Functional_Li for( int ir=0; ir< nrxx; ++ir ) { #ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) + if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { - vtau[ir*nspin+is] *= (1.0 - XC_Functional::hybrid_alpha); + vtau[ir*nspin+is] *= (1.0 - XC_Functional::get_hybrid_alpha()); } #endif vofk(is,ir) += vtau[ir*nspin+is] * sgn[ir*nspin+is]; diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp index 0826f10c5f..48bd301ae6 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp @@ -25,7 +25,7 @@ void XC_Functional_Libxc::tau_xc( { xc_mgga_exc_vxc(&func,1,&rho,&grho,&lapl_rho,&atau,&s,&v1,&v2,&vlapl_rho,&v3); #ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) + if (func.info->number == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) { s *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha); v1 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha); diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index dd6063d5d0..63e4539e3f 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -469,7 +469,7 @@ void Input_Conv::Convert() // EXX case, convert all EXX related variables // GlobalC::exx_info.info_global.cal_exx = true; GlobalC::exx_info.info_global.hybrid_alpha = std::stod(PARAM.inp.exx_hybrid_alpha); - XC_Functional::get_hybrid_alpha(std::stod(PARAM.inp.exx_hybrid_alpha)); + XC_Functional::set_hybrid_alpha(std::stod(PARAM.inp.exx_hybrid_alpha)); GlobalC::exx_info.info_global.hse_omega = PARAM.inp.exx_hse_omega; GlobalC::exx_info.info_global.separate_loop = PARAM.inp.exx_separate_loop; GlobalC::exx_info.info_global.hybrid_step = PARAM.inp.exx_hybrid_step;