From 18fd6ac49db7d005f73c045f53018005d718d97c Mon Sep 17 00:00:00 2001 From: linpz Date: Sun, 8 Sep 2024 11:57:39 +0800 Subject: [PATCH 1/7] Refactor: move functions with libxc in class XC_Functional to namespace XC_Functional_Libxc --- source/module_elecstate/potentials/pot_xc.cpp | 6 +- .../module_surchem/test/CMakeLists.txt | 12 +- .../module_xc/CMakeLists.txt | 6 +- .../module_xc/test/CMakeLists.txt | 37 +- .../module_xc/xc_functional.cpp | 106 +-- .../module_xc/xc_functional.h | 58 +- .../module_xc/xc_functional_gradcorr.cpp | 34 +- .../module_xc/xc_functional_libxc.cpp | 101 +++ .../module_xc/xc_functional_libxc.h | 99 +++ .../module_xc/xc_functional_libxc_vxc.cpp | 614 +++++++++++++++++ .../xc_functional_libxc_wrapper_gcxc.cpp | 91 +++ ... => xc_functional_libxc_wrapper_tauxc.cpp} | 21 +- .../xc_functional_libxc_wrapper_xc.cpp | 38 ++ .../module_xc/xc_functional_vxc.cpp | 615 +----------------- .../module_xc/xc_functional_wrapper_gcxc.cpp | 97 +-- .../module_xc/xc_functional_wrapper_xc.cpp | 35 - .../hamilt_pwdft/forces_cc.cpp | 5 +- .../hamilt_pwdft/stress_func_cc.cpp | 7 +- source/module_lr/potentials/kernel_xc.cpp | 8 +- 19 files changed, 1066 insertions(+), 924 deletions(-) create mode 100644 source/module_hamilt_general/module_xc/xc_functional_libxc.cpp create mode 100644 source/module_hamilt_general/module_xc/xc_functional_libxc.h create mode 100644 source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp create mode 100644 source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp rename source/module_hamilt_general/module_xc/{xc_functional_wrapper_tauxc.cpp => xc_functional_libxc_wrapper_tauxc.cpp} (85%) create mode 100644 source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp diff --git a/source/module_elecstate/potentials/pot_xc.cpp b/source/module_elecstate/potentials/pot_xc.cpp index e2df13b7d6..b04370b2f5 100644 --- a/source/module_elecstate/potentials/pot_xc.cpp +++ b/source/module_elecstate/potentials/pot_xc.cpp @@ -3,6 +3,10 @@ #include "module_base/timer.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#ifdef USE_LIBXC +#include "module_hamilt_general/module_xc/xc_functional_libxc.h" +#endif + namespace elecstate { @@ -20,7 +24,7 @@ void PotXC::cal_v_eff(const Charge*const chg, const UnitCell*const ucell, Module { #ifdef USE_LIBXC const std::tuple etxc_vtxc_v - = XC_Functional::v_xc_meta(nrxx_current, ucell->omega, ucell->tpiba, chg); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), nrxx_current, ucell->omega, ucell->tpiba, chg); *(this->etxc_) = std::get<0>(etxc_vtxc_v); *(this->vtxc_) = std::get<1>(etxc_vtxc_v); v_eff += std::get<2>(etxc_vtxc_v); diff --git a/source/module_hamilt_general/module_surchem/test/CMakeLists.txt b/source/module_hamilt_general/module_surchem/test/CMakeLists.txt index 759d207fb7..1ddabbf955 100644 --- a/source/module_hamilt_general/module_surchem/test/CMakeLists.txt +++ b/source/module_hamilt_general/module_surchem/test/CMakeLists.txt @@ -33,7 +33,11 @@ AddTest( SOURCES cal_vcav_test.cpp ../cal_vcav.cpp ../surchem.cpp ../../../module_hamilt_pw/hamilt_pwdft/parallel_grid.cpp ../../module_xc/xc_functional_gradcorr.cpp ../../module_xc/xc_functional.cpp ../../module_xc/xc_functional_wrapper_xc.cpp ../../module_xc/xc_functional_wrapper_gcxc.cpp - ../../module_xc/xc_functional_wrapper_tauxc.cpp + ../../module_xc/xc_functional_libxc.cpp + ../../module_xc/xc_functional_libxc_vxc.cpp + ../../module_xc/xc_functional_libxc_wrapper_xc.cpp + ../../module_xc/xc_functional_libxc_wrapper_gcxc.cpp + ../../module_xc/xc_functional_libxc_wrapper_tauxc.cpp ../../module_xc/xc_funct_corr_gga.cpp ../../module_xc/xc_funct_corr_lda.cpp ../../module_xc/xc_funct_exch_gga.cpp ../../module_xc/xc_funct_exch_lda.cpp ../../module_xc/xc_funct_hcth.cpp ) @@ -44,7 +48,11 @@ AddTest( SOURCES cal_vel_test.cpp ../cal_vel.cpp ../surchem.cpp ../cal_epsilon.cpp ../minimize_cg.cpp ../../../module_hamilt_pw/hamilt_pwdft/parallel_grid.cpp ../../module_xc/xc_functional_gradcorr.cpp ../../module_xc/xc_functional.cpp ../../module_xc/xc_functional_wrapper_xc.cpp ../../module_xc/xc_functional_wrapper_gcxc.cpp - ../../module_xc/xc_functional_wrapper_tauxc.cpp + ../../module_xc/xc_functional_libxc.cpp + ../../module_xc/xc_functional_libxc_vxc.cpp + ../../module_xc/xc_functional_libxc_wrapper_xc.cpp + ../../module_xc/xc_functional_libxc_wrapper_gcxc.cpp + ../../module_xc/xc_functional_libxc_wrapper_tauxc.cpp ../../module_xc/xc_funct_corr_gga.cpp ../../module_xc/xc_funct_corr_lda.cpp ../../module_xc/xc_funct_exch_gga.cpp ../../module_xc/xc_funct_exch_lda.cpp ../../module_xc/xc_funct_hcth.cpp ) \ No newline at end of file diff --git a/source/module_hamilt_general/module_xc/CMakeLists.txt b/source/module_hamilt_general/module_xc/CMakeLists.txt index 195f144eb4..6f179ba2e2 100644 --- a/source/module_hamilt_general/module_xc/CMakeLists.txt +++ b/source/module_hamilt_general/module_xc/CMakeLists.txt @@ -6,12 +6,16 @@ add_library( xc_functional_gradcorr.cpp xc_functional_wrapper_xc.cpp xc_functional_wrapper_gcxc.cpp - xc_functional_wrapper_tauxc.cpp xc_funct_exch_lda.cpp xc_funct_corr_lda.cpp xc_funct_exch_gga.cpp xc_funct_corr_gga.cpp xc_funct_hcth.cpp + xc_functional_libxc.cpp + xc_functional_libxc_vxc.cpp + xc_functional_libxc_wrapper_xc.cpp + xc_functional_libxc_wrapper_gcxc.cpp + xc_functional_libxc_wrapper_tauxc.cpp ) if(ENABLE_COVERAGE) diff --git a/source/module_hamilt_general/module_xc/test/CMakeLists.txt b/source/module_hamilt_general/module_xc/test/CMakeLists.txt index e0f114a9ab..807a63207e 100644 --- a/source/module_hamilt_general/module_xc/test/CMakeLists.txt +++ b/source/module_hamilt_general/module_xc/test/CMakeLists.txt @@ -257,7 +257,11 @@ AddTest( LIBS MPI::MPI_CXX Libxc::xc ${math_libs} psi device container parameter SOURCES test_xc3.cpp ../xc_functional_gradcorr.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp - ../xc_functional_wrapper_tauxc.cpp + ../xc_functional_libxc.cpp + ../xc_functional_libxc_vxc.cpp + ../xc_functional_libxc_wrapper_xc.cpp + ../xc_functional_libxc_wrapper_gcxc.cpp + ../xc_functional_libxc_wrapper_tauxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp ../../../module_base/matrix.cpp @@ -271,7 +275,11 @@ AddTest( LIBS MPI::MPI_CXX Libxc::xc ${math_libs} psi device container parameter SOURCES test_xc3.cpp ../xc_functional_gradcorr.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp - ../xc_functional_wrapper_tauxc.cpp + ../xc_functional_libxc.cpp + ../xc_functional_libxc_vxc.cpp + ../xc_functional_libxc_wrapper_xc.cpp + ../xc_functional_libxc_wrapper_gcxc.cpp + ../xc_functional_libxc_wrapper_tauxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp ../../../module_base/matrix.cpp @@ -284,7 +292,12 @@ AddTest( TARGET XCTest_SCAN LIBS MPI::MPI_CXX Libxc::xc parameter SOURCES test_xc4.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp - ../xc_functional_wrapper_gcxc.cpp ../xc_functional_wrapper_tauxc.cpp + ../xc_functional_wrapper_gcxc.cpp + ../xc_functional_libxc.cpp + ../xc_functional_libxc_vxc.cpp + ../xc_functional_libxc_wrapper_xc.cpp + ../xc_functional_libxc_wrapper_gcxc.cpp + ../xc_functional_libxc_wrapper_tauxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp ) @@ -294,7 +307,11 @@ AddTest( LIBS MPI::MPI_CXX Libxc::xc ${math_libs} psi device container parameter SOURCES test_xc5.cpp ../xc_functional_gradcorr.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp - ../xc_functional_wrapper_tauxc.cpp + ../xc_functional_libxc.cpp + ../xc_functional_libxc_vxc.cpp + ../xc_functional_libxc_wrapper_xc.cpp + ../xc_functional_libxc_wrapper_gcxc.cpp + ../xc_functional_libxc_wrapper_tauxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp ../xc_functional_vxc.cpp @@ -310,7 +327,11 @@ AddTest( LIBS MPI::MPI_CXX Libxc::xc ${math_libs} psi device container parameter SOURCES test_xc5.cpp ../xc_functional_gradcorr.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp - ../xc_functional_wrapper_tauxc.cpp + ../xc_functional_libxc.cpp + ../xc_functional_libxc_vxc.cpp + ../xc_functional_libxc_wrapper_xc.cpp + ../xc_functional_libxc_wrapper_gcxc.cpp + ../xc_functional_libxc_wrapper_tauxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp ../../../module_base/matrix.cpp ../../../module_base/timer.cpp @@ -327,7 +348,11 @@ AddTest( LIBS MPI::MPI_CXX Libxc::xc ${math_libs} psi device container parameter SOURCES test_xc5.cpp ../xc_functional_gradcorr.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp - ../xc_functional_wrapper_tauxc.cpp + ../xc_functional_libxc.cpp + ../xc_functional_libxc_vxc.cpp + ../xc_functional_libxc_wrapper_xc.cpp + ../xc_functional_libxc_wrapper_gcxc.cpp + ../xc_functional_libxc_wrapper_tauxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp ../../../module_base/matrix.cpp ../../../module_base/timer.cpp diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp index 164d497165..7a505a1e90 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional.cpp @@ -5,6 +5,10 @@ #include "module_cell/module_paw/paw_cell.h" #endif +#ifdef USE_LIBXC +#include "xc_functional_libxc.h" +#endif + XC_Functional::XC_Functional(){} XC_Functional::~XC_Functional(){} @@ -198,7 +202,9 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) { #ifdef USE_LIBXC //see if it matches libxc functionals - set_xc_type_libxc(xc_func); + const std::pair> type_id = XC_Functional_Libxc::set_xc_type_libxc(xc_func); + func_type = std::get<0>(type_id); + func_id = std::get<1>(type_id); use_libxc = true; #else ModuleBase::WARNING_QUIT("xc_functional.cpp","functional name not recognized!"); @@ -239,101 +245,3 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) #endif } - -#ifdef USE_LIBXC -void XC_Functional::set_xc_type_libxc(std::string xc_func_in) -{ - - // determine the type (lda/gga/mgga) - func_type = 1; - if(xc_func_in.find("GGA") != std::string::npos) { func_type = 2; -} - if(xc_func_in.find("MGGA") != std::string::npos) { func_type = 3; -} - if(xc_func_in.find("HYB") != std::string::npos) { func_type =4; -} - if(xc_func_in.find("HYB") != std::string::npos && xc_func_in.find("MGGA") != std::string::npos) { func_type =5; -} - - // determine the id - int pos = 0; - std::string delimiter = "+"; - std::string token; - while ((pos = xc_func_in.find(delimiter)) != std::string::npos) - { - token = xc_func_in.substr(0, pos); - int id = xc_functional_get_number(token.c_str()); - std::cout << "func,id" << token << " " << id << std::endl; - if (id == -1) { ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc","functional name not recognized!"); -} - func_id.push_back(id); - xc_func_in.erase(0, pos + delimiter.length()); - } - int id = xc_functional_get_number(xc_func_in.c_str()); - std::cout << "func,id" << xc_func_in << " " << id << std::endl; - if (id == -1) { ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc","functional name not recognized!"); -} - func_id.push_back(id); - -} -#endif - -#ifdef USE_LIBXC -std::vector XC_Functional::init_func(const int xc_polarized) -{ - // 'funcs' is the return value - std::vector funcs; - - //------------------------------------------- - // define a function named 'add_func', which - // initialize a functional according to its ID - //------------------------------------------- - auto add_func = [&]( const int func_id ) - { - funcs.push_back({}); - // 'xc_func_init' is defined in Libxc - xc_func_init( &funcs.back(), func_id, xc_polarized ); - }; - - for(int id : func_id) - { - if(id == XC_LDA_XC_KSDT || id == XC_LDA_XC_CORRKSDT || id == XC_LDA_XC_GDSMFB) //finite temperature XC functionals - { - add_func(id); - double parameter_finitet[1] = {PARAM.inp.xc_temperature * 0.5}; // converts to Hartree for libxc - xc_func_set_ext_params(&funcs.back(), parameter_finitet); - } -#ifdef __EXX - else if( id == XC_HYB_GGA_XC_PBEH ) // PBE0 - { - add_func( XC_HYB_GGA_XC_PBEH ); - double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, - GlobalC::exx_info.info_global.hse_omega, - GlobalC::exx_info.info_global.hse_omega }; - xc_func_set_ext_params(&funcs.back(), parameter_hse); - } - else if( id == XC_HYB_GGA_XC_HSE06 ) // HSE06 hybrid functional - { - add_func( XC_HYB_GGA_XC_HSE06 ); - double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, - GlobalC::exx_info.info_global.hse_omega, - GlobalC::exx_info.info_global.hse_omega }; - xc_func_set_ext_params(&funcs.back(), parameter_hse); - } -#endif - else - { - add_func( id ); - } - } - return funcs; -} - -void XC_Functional::finish_func(std::vector &funcs) -{ - for(xc_func_type func : funcs) - { - xc_func_end(&func); - } -} -#endif diff --git a/source/module_hamilt_general/module_xc/xc_functional.h b/source/module_hamilt_general/module_xc/xc_functional.h index 1789831f51..5179c6bf93 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.h +++ b/source/module_hamilt_general/module_xc/xc_functional.h @@ -49,20 +49,6 @@ class XC_Functional const Charge* const chr, const UnitCell *ucell); // charge density - // using libxc - static std::tuple v_xc_libxc( - const int &nrxx, // number of real-space grid - const double &omega, // volume of cell - const double tpiba, - const Charge* const chr); // charge density - - // for mGGA functional - static std::tuple v_xc_meta( - const int &nrxx, // number of real-space grid - const double &omega, // volume of cell - const double tpiba, - const Charge* const chr); - //------------------- // xc_functional.cpp //------------------- @@ -75,9 +61,6 @@ class XC_Functional // func_id, which is the LIBXC id of functional // func_type, which is as specified in get_func_type // use_libxc, whether to use LIBXC. The rule is to NOT use it for functionals that we already have. -// 3. set_xc_type_libxc : sets functional type, which allows combination of LIBXC keyword connected by "+" -// for example, "XC_LDA_X+XC_LDA_C_PZ" -// 4. init_func : which converts func_id into corresponding xc_func_type vector static int get_func_type(); static void set_xc_type(const std::string xc_func_in); @@ -86,11 +69,6 @@ class XC_Functional static void get_hybrid_alpha(const double alpha_in); /// Usually in exx caculation, the first SCF loop should be converged with PBE static void set_xc_first_loop(const UnitCell& ucell); -#ifdef USE_LIBXC - static void set_xc_type_libxc(const std::string xc_func_in); - static std::vector init_func(const int xc_polarized); - static void finish_func(std::vector &funcs); -#endif private: @@ -101,6 +79,9 @@ class XC_Functional //exx_hybrid_alpha for mixing exx in hybrid functional: static double hybrid_alpha; + public: + static std::vector get_func_id() { return func_id; } + //------------------- // xc_functional_wrapper_xc.cpp //------------------- @@ -133,8 +114,6 @@ class XC_Functional // LSDA static void xc_spin(const double &rho, const double &zeta, double &exc, double &vxcup, double &vxcdw); - static void xc_spin_libxc(const double &rhoup, const double &rhodw, - double &exc, double &vxcup, double &vxcdw); //------------------- // xc_functional_wrapper_gcxc.cpp @@ -145,8 +124,6 @@ class XC_Functional // 1. gcxc, which is the wrapper for gradient correction part // 2. gcx_spin, spin polarized, exchange only // 3. gcc_spin, spin polarized, correlation only -// 4. gcxc_libxc, the entire GGA functional, LIBXC, for nspin=1 case -// 5. gcxc_spin_libxc, the entire GGA functional, LIBXC, for nspin=2 case // The difference between our realization (gcxc/gcx_spin/gcc_spin) and // LIBXC, and the reason for not having gcxc_libxc is explained @@ -155,8 +132,6 @@ class XC_Functional // GGA static void gcxc(const double &rho, const double &grho, double &sxc, double &v1xc, double &v2xc); - static void gcxc_libxc(const double &rho, const double &grho, - double &sxc, double &v1xc, double &v2xc); // spin polarized GGA static void gcx_spin(double rhoup, double rhodw, double grhoup2, double grhodw2, @@ -165,33 +140,6 @@ class XC_Functional static void gcc_spin(double rho, double &zeta, double grho, double &sc, double &v1cup, double &v1cdw, double &v2c); - static void gcxc_spin_libxc(double rhoup, double rhodw, - ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, - double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud); - -//------------------- -// xc_functional_wrapper_tauxc.cpp -//------------------- - -// This file contains wrapper for the mGGA functionals -// it includes 1 subroutine: -// 1. tau_xc -// 2. tau_xc_spin - -// NOTE : mGGA is realized through LIBXC - -#ifdef USE_LIBXC - // mGGA - static void tau_xc(const double &rho, const double &grho, const double &atau, double &sxc, - double &v1xc, double &v2xc, double &v3xc); - - static void tau_xc_spin(double rhoup, double rhodw, - ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, - double tauup, double taudw, - double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, - double &v3xcup, double &v3xcdw); -#endif - //------------------- // xc_functional_gradcorr.cpp //------------------- diff --git a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp index 4d94f68ec4..30c0fdcc4c 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_gradcorr.cpp @@ -17,6 +17,10 @@ #include #include +#ifdef USE_LIBXC +#include "xc_functional_libxc.h" +#endif + // from gradcorr.f90 void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, const Charge* const chr, ModulePW::PW_Basis* rhopw, const UnitCell *ucell, @@ -24,18 +28,14 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, { ModuleBase::TITLE("XC_Functional","gradcorr"); - if(func_type == 0 || func_type == 1) { return; // none or LDA functional -} + if(func_type == 0 || func_type == 1) { return; } // none or LDA functional bool igcc_is_lyp = false; - if( func_id[1] == XC_GGA_C_LYP) { igcc_is_lyp = true; -} + if( func_id[1] == XC_GGA_C_LYP) { igcc_is_lyp = true; } int nspin0 = GlobalV::NSPIN; - if(GlobalV::NSPIN==4) { nspin0 =1; -} - if(GlobalV::NSPIN==4&&(GlobalV::DOMAG||GlobalV::DOMAG_Z)) { nspin0 = 2; -} + if(GlobalV::NSPIN==4) { nspin0 =1; } + if(GlobalV::NSPIN==4&&(GlobalV::DOMAG||GlobalV::DOMAG_Z)) { nspin0 = 2; } assert(nspin0>0); const double fac = 1.0/ nspin0; @@ -91,8 +91,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, } gdr1 = new ModuleBase::Vector3[rhopw->nrxx]; - if(!is_stress) { h1 = new ModuleBase::Vector3[rhopw->nrxx]; -} + if(!is_stress) { h1 = new ModuleBase::Vector3[rhopw->nrxx]; } XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); @@ -118,8 +117,7 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, } gdr2 = new ModuleBase::Vector3[rhopw->nrxx]; - if(!is_stress) { h2 = new ModuleBase::Vector3[rhopw->nrxx]; -} + if(!is_stress) { h2 = new ModuleBase::Vector3[rhopw->nrxx]; } XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); } @@ -253,11 +251,11 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, { double v3xc; double atau = chr->kin_r[0][ir]/2.0; - XC_Functional::tau_xc( arho, grho2a, atau, sxc, v1xc, v2xc, v3xc); + XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc); } else { - XC_Functional::gcxc_libxc( arho, grho2a, sxc, v1xc, v2xc); + XC_Functional_Libxc::gcxc_libxc( func_id, arho, grho2a, sxc, v1xc, v2xc); } #endif } // end use_libxc @@ -313,12 +311,16 @@ void XC_Functional::gradcorr(double &etxc, double &vtxc, ModuleBase::matrix &v, double v3xcup, v3xcdw; double atau1 = chr->kin_r[0][ir]/2.0; double atau2 = chr->kin_r[1][ir]/2.0; - XC_Functional::tau_xc_spin( rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], + XC_Functional_Libxc::tau_xc_spin( + func_id, + rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], atau1, atau2, sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud, v3xcup, v3xcdw); } else { - XC_Functional::gcxc_spin_libxc(rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], + XC_Functional_Libxc::gcxc_spin_libxc( + func_id, + rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud); } if(is_stress) diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp new file mode 100644 index 0000000000..416b08b65e --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp @@ -0,0 +1,101 @@ +#ifdef USE_LIBXC + +#include "xc_functional_libxc.h" +#include "module_parameter/parameter.h" +#include "module_base/tool_quit.h" + +#include + +#include + +std::pair> XC_Functional_Libxc::set_xc_type_libxc(std::string xc_func_in) +{ + // determine the type (lda/gga/mgga) + int func_type; //0:none, 1:lda, 2:gga, 3:mgga, 4:hybrid lda/gga, 5:hybrid mgga + func_type = 1; + if(xc_func_in.find("GGA") != std::string::npos) { func_type = 2; } + if(xc_func_in.find("MGGA") != std::string::npos) { func_type = 3; } + if(xc_func_in.find("HYB") != std::string::npos) { func_type =4; } + if(xc_func_in.find("HYB") != std::string::npos && xc_func_in.find("MGGA") != std::string::npos) { func_type =5; } + + // determine the id + std::vector func_id; // libxc id of functional + int pos = 0; + std::string delimiter = "+"; + std::string token; + while ((pos = xc_func_in.find(delimiter)) != std::string::npos) + { + token = xc_func_in.substr(0, pos); + int id = xc_functional_get_number(token.c_str()); + std::cout << "func,id" << token << " " << id << std::endl; + if (id == -1) { ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc","functional name not recognized!"); } + func_id.push_back(id); + xc_func_in.erase(0, pos + delimiter.length()); + } + int id = xc_functional_get_number(xc_func_in.c_str()); + std::cout << "func,id" << xc_func_in << " " << id << std::endl; + if (id == -1) { ModuleBase::WARNING_QUIT("XC_Functional::set_xc_type_libxc","functional name not recognized!"); } + func_id.push_back(id); + + return std::make_pair(func_type, func_id); +} + +std::vector XC_Functional_Libxc::init_func(const std::vector &func_id, const int xc_polarized) +{ + // 'funcs' is the return value + std::vector funcs; + + //------------------------------------------- + // define a function named 'add_func', which + // initialize a functional according to its ID + //------------------------------------------- + auto add_func = [&]( const int func_id ) + { + funcs.push_back({}); + // 'xc_func_init' is defined in Libxc + xc_func_init( &funcs.back(), func_id, xc_polarized ); + }; + + for(int id : func_id) + { + if(id == XC_LDA_XC_KSDT || id == XC_LDA_XC_CORRKSDT || id == XC_LDA_XC_GDSMFB) //finite temperature XC functionals + { + add_func(id); + double parameter_finitet[1] = {PARAM.inp.xc_temperature * 0.5}; // converts to Hartree for libxc + xc_func_set_ext_params(&funcs.back(), parameter_finitet); + } +#ifdef __EXX + else if( id == XC_HYB_GGA_XC_PBEH ) // PBE0 + { + add_func( XC_HYB_GGA_XC_PBEH ); + double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, + GlobalC::exx_info.info_global.hse_omega, + GlobalC::exx_info.info_global.hse_omega }; + xc_func_set_ext_params(&funcs.back(), parameter_hse); + } + else if( id == XC_HYB_GGA_XC_HSE06 ) // HSE06 hybrid functional + { + add_func( XC_HYB_GGA_XC_HSE06 ); + double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, + GlobalC::exx_info.info_global.hse_omega, + GlobalC::exx_info.info_global.hse_omega }; + xc_func_set_ext_params(&funcs.back(), parameter_hse); + } +#endif + else + { + add_func( id ); + } + } + return funcs; +} + +void XC_Functional_Libxc::finish_func(std::vector &funcs) +{ + for(xc_func_type func : funcs) + { + xc_func_end(&func); + } +} + +#endif \ No newline at end of file diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc.h b/source/module_hamilt_general/module_xc/xc_functional_libxc.h new file mode 100644 index 0000000000..bbe5e416fc --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.h @@ -0,0 +1,99 @@ +#ifndef XC_FUNCTIONAL_LIBXC_H +#define XC_FUNCTIONAL_LIBXC_H + +#ifdef USE_LIBXC + +#include "module_base/matrix.h" +#include "module_base/vector3.h" + +#include + +#include +#include + + class Charge; + +namespace XC_Functional_Libxc +{ +//------------------- +// xc_functional_libxc.cpp +//------------------- + + // sets functional type, which allows combination of LIBXC keyword connected by "+" + // for example, "XC_LDA_X+XC_LDA_C_PZ" + std::pair> set_xc_type_libxc(const std::string xc_func_in); + + // converts func_id into corresponding xc_func_type vector + std::vector init_func(const std::vector &func_id, const int xc_polarized); + + void finish_func(std::vector &funcs); + +//------------------- +// xc_functional_libxc_vxc.cpp +//------------------- + + 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 + + // for mGGA functional + std::tuple v_xc_meta( + 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); + +//------------------- +// xc_functional_libxc_wrapper_xc.cpp +//------------------- + + void xc_spin_libxc( + const std::vector &func_id, + const double &rhoup, const double &rhodw, + double &exc, double &vxcup, double &vxcdw); + + +//------------------- +// xc_functional_libxc_wrapper_gcxc.cpp +//------------------- + + // the entire GGA functional, for nspin=1 case + void gcxc_libxc( + const std::vector &func_id, + const double &rho, const double &grho, + double &sxc, double &v1xc, double &v2xc); + + // the entire GGA functional, for nspin=2 case + void gcxc_spin_libxc( + const std::vector &func_id, + double rhoup, double rhodw, + ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud); + + +//------------------- +// xc_functional_libxc_wrapper_tauxc.cpp +//------------------- + + // wrapper for the mGGA functionals + void tau_xc( + const std::vector &func_id, + const double &rho, const double &grho, const double &atau, double &sxc, + double &v1xc, double &v2xc, double &v3xc); + + void tau_xc_spin( + const std::vector &func_id, + double rhoup, double rhodw, + ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + double tauup, double taudw, + double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, + double &v3xcup, double &v3xcdw); +} + +#endif + +#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 new file mode 100644 index 0000000000..7c355e2c00 --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp @@ -0,0 +1,614 @@ +#ifdef USE_LIBXC + +#include "xc_functional.h" +#include "xc_functional_libxc.h" +#include "module_elecstate/module_charge/charge.h" +#include "module_base/global_variable.h" +#include "module_base/parallel_reduce.h" +#include "module_base/timer.h" +#include "module_base/tool_title.h" + +#include + +#include + +std::tuple XC_Functional_Libxc::v_xc_libxc( // Peize Lin update for nspin==4 at 2023.01.14 + 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) +{ + ModuleBase::TITLE("XC_Functional","v_xc_libxc"); + ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); + + const int nspin = + (GlobalV::NSPIN == 1 || ( GlobalV::NSPIN ==4 && !GlobalV::DOMAG && !GlobalV::DOMAG_Z)) + ? 1 : 2; + + double etxc = 0.0; + double vtxc = 0.0; + ModuleBase::matrix v(nspin,nrxx); + + //---------------------------------------------------------- + // xc_func_type is defined in Libxc package + // to understand the usage of xc_func_type, + // use can check on website, for example: + // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ + //---------------------------------------------------------- + + std::vector funcs = XC_Functional_Libxc::init_func( func_id, (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED ); + + const bool is_gga = [&funcs]() + { + for( xc_func_type &func : funcs ) + { + switch( func.info->family ) + { + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: + return true; + } + } + return false; + }(); + + // converting rho + std::vector rho(nrxx*nspin); + std::vector amag; + if(1==nspin || 2==GlobalV::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]; + } + else + { + amag.resize(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]>> gdr; + std::vector sigma; + if(is_gga) + { + // calculating grho + gdr.resize( nspin ); + for( int is=0; is!=nspin; ++is ) + { + std::vector rhor(nrxx); + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for(int 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) + + // converting grho + sigma.resize( nrxx * ((1==nspin)?1:3) ); + if( 1==nspin ) + { + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for( int ir=0; ir 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; ir exc ( nrxx ); + std::vector vrho ( nrxx * nspin ); + std::vector vsigma( nrxx * ((1==nspin)?1:3) ); + switch( func.info->family ) + { + case XC_FAMILY_LDA: + // call Libxc function: xc_lda_exc_vxc + xc_lda_exc_vxc( &func, nrxx, rho.data(), + exc.data(), vrho.data() ); + break; + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: + // call Libxc function: xc_gga_exc_vxc + xc_gga_exc_vxc( &func, nrxx, rho.data(), sigma.data(), + exc.data(), vrho.data(), vsigma.data() ); + break; + default: + throw std::domain_error("func.info->family ="+std::to_string(func.info->family) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + break; + } + + #ifdef _OPENMP + #pragma omp parallel for collapse(2) reduction(+:etxc) schedule(static, 256) + #endif + for( int is=0; isfamily == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) + { + std::vector>> h( nspin, + std::vector>(nrxx) ); + if( 1==nspin ) + { + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for( int ir=0; ir< nrxx; ++ir ) + h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir]; + } + else + { + #ifdef _OPENMP + #pragma omp parallel for schedule(static, 1024) + #endif + for( int ir=0; ir< nrxx; ++ir ) + { + h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 + + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 + + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + } + } + + // define two dimensional array dh [ nspin, nrxx ] + std::vector> 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); + + 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)) + } // end for( xc_func_type &func : funcs ) + + //------------------------------------------------- + // for MPI, reduce the exchange-correlation energy + //------------------------------------------------- + #ifdef __MPI + Parallel_Reduce::reduce_pool(etxc); + Parallel_Reduce::reduce_pool(vtxc); + #endif + + etxc *= omega / chr->rhopw->nxyz; + vtxc *= omega / chr->rhopw->nxyz; + + finish_func(funcs); + + if(1==GlobalV::NSPIN || 2==GlobalV::NSPIN) + { + ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); + return std::make_tuple( etxc, vtxc, std::move(v) ); + } + else if(4==GlobalV::NSPIN) + { + 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; ipol<4; ++ipol) + v_nspin4(ipol,ir) = vs * chr->rho[ipol][ir] / amag[ir]; + } + } + } + + ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); + return std::make_tuple( etxc, vtxc, std::move(v_nspin4) ); + } // end if(4==GlobalV::NSPIN) + else//NSPIN != 1,2,4 is not supported + { + throw std::domain_error("GlobalV::NSPIN ="+std::to_string(GlobalV::NSPIN) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + } +} + + +//the interface to libxc xc_mgga_exc_vxc(xc_func,n,rho,grho,laplrho,tau,e,v1,v2,v3,v4) +//xc_func : LIBXC data type, contains information on xc functional +//n: size of array, nspin*nnr +//rho,grho,laplrho: electron density, its gradient and laplacian +//tau(kin_r): kinetic energy density +//e: energy density +//v1-v4: derivative of energy density w.r.t rho, gradient, laplacian and tau +//v1 and v2 are combined to give v; v4 goes into vofk + +//XC_POLARIZED, XC_UNPOLARIZED: internal flags used in LIBXC, denote the polarized(nspin=1) or unpolarized(nspin=2) calculations, definition can be found in xc.h from LIBXC + +// [etxc, vtxc, v, vofk] = XC_Functional::v_xc(...) +std::tuple XC_Functional_Libxc::v_xc_meta( + 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) +{ + ModuleBase::TITLE("XC_Functional","v_xc_meta"); + ModuleBase::timer::tick("XC_Functional","v_xc_meta"); + + double e2 = 2.0; + + //output of the subroutine + double etxc = 0.0; + double vtxc = 0.0; + ModuleBase::matrix v(GlobalV::NSPIN,nrxx); + ModuleBase::matrix vofk(GlobalV::NSPIN,nrxx); + + //---------------------------------------------------------- + // xc_func_type is defined in Libxc package + // to understand the usage of xc_func_type, + // use can check on website, for example: + // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ + //---------------------------------------------------------- + + const int nspin = GlobalV::NSPIN; + std::vector funcs = XC_Functional_Libxc::init_func( func_id, ( (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED ) ); + + // converting rho + std::vector rho; + rho.resize(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]; + } + } + + std::vector>> gdr; + std::vector sigma; + + // calculating grho + gdr.resize( nspin ); + for( int is=0; is!=nspin; ++is ) + { + std::vector rhor(nrxx); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int 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); + } + + // converting grho + sigma.resize( nrxx * ((1==nspin)?1:3) ); + + if( 1==nspin ) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for( int ir=0; ir kin_r; + kin_r.resize(nrxx*nspin); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) schedule(static, 1024) +#endif + for( int is=0; iskin_r[is][ir] / 2.0; + } + } + + std::vector exc ( nrxx ); + std::vector vrho ( nrxx * nspin ); + std::vector vsigma ( nrxx * ((1==nspin)?1:3) ); + std::vector vtau ( nrxx * nspin ); + std::vector vlapl ( nrxx * nspin ); + + const double rho_th = 1e-8; + const double grho_th = 1e-12; + const double tau_th = 1e-8; + // sgn for threshold mask + std::vector sgn( nrxx * nspin); +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for(int i = 0; i < nrxx * nspin; ++i) + { + sgn[i] = 1.0; + } + + if(nspin == 1) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for( int ir=0; irfamily == XC_FAMILY_MGGA); + xc_mgga_exc_vxc(&func, nrxx, rho.data(), sigma.data(), sigma.data(), + kin_r.data(), exc.data(), vrho.data(), vsigma.data(), vlapl.data(), vtau.data()); + + //process etxc + for( int is=0; is!=nspin; ++is ) + { +#ifdef _OPENMP +#pragma omp parallel for reduction(+:etxc) schedule(static, 256) +#endif + for( int ir=0; ir< nrxx; ++ir ) + { +#ifdef __EXX + if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) + { + exc[ir] *= (1.0 - XC_Functional::hybrid_alpha); + } +#endif + etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is]; + } + } + + //process vtxc +#ifdef _OPENMP +#pragma omp parallel for collapse(2) reduction(+:vtxc) schedule(static, 256) +#endif + for( int is=0; isnumber == XC_MGGA_X_SCAN && get_func_type() == 5) + { + vrho[ir*nspin+is] *= (1.0 - XC_Functional::hybrid_alpha); + } +#endif + const double v_tmp = ModuleBase::e2 * vrho[ir*nspin+is] * sgn[ir*nspin+is]; + v(is,ir) += v_tmp; + vtxc += v_tmp * chr->rho[is][ir]; + } + } + + //process vsigma + std::vector>> h( nspin, std::vector>(nrxx) ); + if( 1==nspin ) + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 1024) +#endif + for( int ir=0; ir< nrxx; ++ir ) + { +#ifdef __EXX + if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) + { + vsigma[ir] *= (1.0 - XC_Functional::hybrid_alpha); + } +#endif + h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir]; + } + } + else + { +#ifdef _OPENMP +#pragma omp parallel for schedule(static, 64) +#endif + for( int ir=0; ir< nrxx; ++ir ) + { +#ifdef __EXX + if (func.info->number == XC_MGGA_X_SCAN && 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); + } +#endif + h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 + + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 + + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + } + } + + // define two dimensional array dh [ nspin, nrxx ] + std::vector> dh(nspin, std::vector( nrxx)); + for( int is=0; is!=nspin; ++is ) + { + XC_Functional::grad_dot( ModuleBase::GlobalFunc::VECTOR_TO_PTR(h[is]), + ModuleBase::GlobalFunc::VECTOR_TO_PTR(dh[is]), chr->rhopw, + tpiba); + } + + double rvtxc = 0.0; +#ifdef _OPENMP +#pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) +#endif + for( int is=0; isnumber == XC_MGGA_X_SCAN && get_func_type() == 5) + { + vtau[ir*nspin+is] *= (1.0 - XC_Functional::hybrid_alpha); + } +#endif + vofk(is,ir) += vtau[ir*nspin+is] * sgn[ir*nspin+is]; + } + } + } + + //------------------------------------------------- + // for MPI, reduce the exchange-correlation energy + //------------------------------------------------- +#ifdef __MPI + Parallel_Reduce::reduce_pool(etxc); + Parallel_Reduce::reduce_pool(vtxc); +#endif + + etxc *= omega / chr->rhopw->nxyz; + vtxc *= omega / chr->rhopw->nxyz; + + XC_Functional_Libxc::finish_func(funcs); + + ModuleBase::timer::tick("XC_Functional","v_xc_meta"); + return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk) ); +} + +#endif \ No newline at end of file diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp new file mode 100644 index 0000000000..24885a5e37 --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp @@ -0,0 +1,91 @@ +#ifdef USE_LIBXC + +#include "xc_functional_libxc.h" + +#include + +void XC_Functional_Libxc::gcxc_libxc( + const std::vector &func_id, + const double &rho, const double &grho, double &sxc, + double &v1xc, double &v2xc) +{ + const double small = 1.e-6; + const double smallg = 1.e-10; + double s,v1,v2; + sxc = v1xc = v2xc = 0.0; + + if (rho <= small || grho < smallg) + { + return; + } + + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED); + + + for(xc_func_type &func : funcs) + { + xc_gga_exc_vxc(&func, 1, &rho, &grho, &s, &v1, &v2); + + sxc += s * rho; + v2xc += v2 * 2.0; + v1xc += v1; + } + XC_Functional_Libxc::finish_func(funcs); +} // end subroutine gcxc_libxc + +void XC_Functional_Libxc::gcxc_spin_libxc( + const std::vector &func_id, + double rhoup, double rhodw, + ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud) +{ + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); + double *rho, *grho, *v1xc, *v2xc, *sgn, s; + sxc = v1xcup = v1xcdw = 0.0; + v2xcup = v2xcdw = v2xcud = 0.0; + rho = new double[2]; + grho= new double[3]; + v1xc= new double[2]; + v2xc= new double[3]; + sgn = new double[2]; + + rho[0] = rhoup; + rho[1] = rhodw; + grho[0] = gdr1.norm2(); + grho[1] = gdr1 * gdr2; + grho[2] = gdr2.norm2(); + + const double rho_threshold = 1E-6; + const double grho_threshold = 1E-10; + + for(xc_func_type &func : funcs) + { + if( func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) + { + sgn[0] = sgn[1] = 1.0; + // call Libxc function: xc_gga_exc_vxc + xc_gga_exc_vxc( &func, 1, rho, grho, &s, v1xc, v2xc); + if(func.info->kind==XC_CORRELATION) + { + if ( rho[0] &func_id, + const double &rho, const double &grho, const double &atau, double &sxc, double &v1xc, double &v2xc, double &v3xc) { double s, v1, v2, v3; double lapl_rho, vlapl_rho; lapl_rho = grho; - std::vector funcs = init_func(XC_UNPOLARIZED); + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED); sxc = 0.0; v1xc = 0.0; v2xc = 0.0; v3xc = 0.0; @@ -35,20 +38,22 @@ void XC_Functional::tau_xc(const double &rho, const double &grho, const double & v1xc += v1; v3xc += v3; } - finish_func(funcs); + XC_Functional_Libxc::finish_func(funcs); return; } -void XC_Functional::tau_xc_spin(double rhoup, double rhodw, +void XC_Functional_Libxc::tau_xc_spin( + const std::vector &func_id, + double rhoup, double rhodw, ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, double tauup, double taudw, double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, double &v3xcup, double &v3xcdw) { - std::vector funcs = init_func(XC_POLARIZED); + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); double *rho, *grho, *tau, *v1xc, *v2xc, *v3xc, *sgn, s; sxc = v1xcup = v1xcdw = 0.0; @@ -108,7 +113,7 @@ void XC_Functional::tau_xc_spin(double rhoup, double rhodw, delete[] v2xc; delete[] v3xc; delete[] sgn; - finish_func(funcs); + XC_Functional_Libxc::finish_func(funcs); return; } diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp new file mode 100644 index 0000000000..c5693ec54d --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp @@ -0,0 +1,38 @@ +#ifdef USE_LIBXC + +#include "xc_functional_libxc.h" + +void XC_Functional_Libxc::xc_spin_libxc( + const std::vector &func_id, + const double &rhoup, const double &rhodw, + double &exc, double &vxcup, double &vxcdw) +{ + double e, vup, vdw; + double *rho_ud, *vxc_ud; + exc = vxcup = vxcdw = 0.0; + + rho_ud = new double[2]; + vxc_ud = new double[2]; + rho_ud[0] = rhoup; + rho_ud[1] = rhodw; + + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); + + for(xc_func_type &func : funcs) + { + if( func.info->family == XC_FAMILY_LDA) + { + // call Libxc function: xc_lda_exc_vxc + xc_lda_exc_vxc( &func, 1, rho_ud, &e, vxc_ud); + } + exc += e; + vxcup += vxc_ud[0]; + vxcdw += vxc_ud[1]; + } + + XC_Functional_Libxc::finish_func(funcs); + delete[] rho_ud; + delete[] vxc_ud; +} + +#endif \ No newline at end of file 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 be64d92962..c76fc7f170 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp @@ -8,6 +8,10 @@ #include "module_base/parallel_reduce.h" #include "module_base/timer.h" +#ifdef USE_LIBXC +#include "xc_functional_libxc.h" +#endif + // [etxc, vtxc, v] = XC_Functional::v_xc(...) std::tuple XC_Functional::v_xc( const int &nrxx, // number of real-space grid @@ -20,7 +24,7 @@ std::tuple XC_Functional::v_xc( if(use_libxc) { #ifdef USE_LIBXC - return v_xc_libxc(nrxx, ucell->omega, ucell->tpiba, chr); + return XC_Functional_Libxc::v_xc_libxc(XC_Functional::get_func_id(), nrxx, ucell->omega, ucell->tpiba, chr); #else ModuleBase::WARNING_QUIT("v_xc","compile with LIBXC"); #endif @@ -124,9 +128,13 @@ std::tuple XC_Functional::v_xc( if(use_libxc) { +#ifdef USE_LIBXC double rhoup = arhox * (1.0+zeta) / 2.0; double rhodw = arhox * (1.0-zeta) / 2.0; - XC_Functional::xc_spin_libxc(rhoup, rhodw, exc, vxc[0], vxc[1]); + XC_Functional_Libxc::xc_spin_libxc(XC_Functional::get_func_id(), rhoup, rhodw, exc, vxc[0], vxc[1]); +#else + ModuleBase::WARNING_QUIT("v_xc","compile with LIBXC"); +#endif } else { @@ -173,605 +181,4 @@ std::tuple XC_Functional::v_xc( ModuleBase::timer::tick("XC_Functional","v_xc"); return std::make_tuple(etxc, vtxc, std::move(v)); -} - -#ifdef USE_LIBXC - -std::tuple XC_Functional::v_xc_libxc( // Peize Lin update for nspin==4 at 2023.01.14 - const int &nrxx, // number of real-space grid - const double &omega, // volume of cell - const double tpiba, - const Charge* const chr) -{ - ModuleBase::TITLE("XC_Functional","v_xc_libxc"); - ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); - - const int nspin = - (GlobalV::NSPIN == 1 || ( GlobalV::NSPIN ==4 && !GlobalV::DOMAG && !GlobalV::DOMAG_Z)) - ? 1 : 2; - - double etxc = 0.0; - double vtxc = 0.0; - ModuleBase::matrix v(nspin,nrxx); - - //---------------------------------------------------------- - // xc_func_type is defined in Libxc package - // to understand the usage of xc_func_type, - // use can check on website, for example: - // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ - //---------------------------------------------------------- - - std::vector funcs = init_func( (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED ); - - const bool is_gga = [&funcs]() - { - for( xc_func_type &func : funcs ) - { - switch( func.info->family ) - { - case XC_FAMILY_GGA: - case XC_FAMILY_HYB_GGA: - return true; - } - } - return false; - }(); - - // converting rho - std::vector rho(nrxx*nspin); - std::vector amag; - if(1==nspin || 2==GlobalV::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]; - } - else - { - amag.resize(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]>> gdr; - std::vector sigma; - if(is_gga) - { - // calculating grho - gdr.resize( nspin ); - for( int is=0; is!=nspin; ++is ) - { - std::vector rhor(nrxx); - #ifdef _OPENMP - #pragma omp parallel for schedule(static, 1024) - #endif - for(int 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) - - // converting grho - sigma.resize( nrxx * ((1==nspin)?1:3) ); - if( 1==nspin ) - { - #ifdef _OPENMP - #pragma omp parallel for schedule(static, 1024) - #endif - for( int ir=0; ir 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; ir exc ( nrxx ); - std::vector vrho ( nrxx * nspin ); - std::vector vsigma( nrxx * ((1==nspin)?1:3) ); - switch( func.info->family ) - { - case XC_FAMILY_LDA: - // call Libxc function: xc_lda_exc_vxc - xc_lda_exc_vxc( &func, nrxx, rho.data(), - exc.data(), vrho.data() ); - break; - case XC_FAMILY_GGA: - case XC_FAMILY_HYB_GGA: - // call Libxc function: xc_gga_exc_vxc - xc_gga_exc_vxc( &func, nrxx, rho.data(), sigma.data(), - exc.data(), vrho.data(), vsigma.data() ); - break; - default: - throw std::domain_error("func.info->family ="+std::to_string(func.info->family) - +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); - break; - } - - #ifdef _OPENMP - #pragma omp parallel for collapse(2) reduction(+:etxc) schedule(static, 256) - #endif - for( int is=0; isfamily == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) - { - std::vector>> h( nspin, - std::vector>(nrxx) ); - if( 1==nspin ) - { - #ifdef _OPENMP - #pragma omp parallel for schedule(static, 1024) - #endif - for( int ir=0; ir< nrxx; ++ir ) - h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir]; - } - else - { - #ifdef _OPENMP - #pragma omp parallel for schedule(static, 1024) - #endif - for( int ir=0; ir< nrxx; ++ir ) - { - h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 - + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); - h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 - + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); - } - } - - // define two dimensional array dh [ nspin, nrxx ] - std::vector> 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); - - 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)) - } // end for( xc_func_type &func : funcs ) - - //------------------------------------------------- - // for MPI, reduce the exchange-correlation energy - //------------------------------------------------- - #ifdef __MPI - Parallel_Reduce::reduce_pool(etxc); - Parallel_Reduce::reduce_pool(vtxc); - #endif - - etxc *= omega / chr->rhopw->nxyz; - vtxc *= omega / chr->rhopw->nxyz; - - finish_func(funcs); - - if(1==GlobalV::NSPIN || 2==GlobalV::NSPIN) - { - ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); - return std::make_tuple( etxc, vtxc, std::move(v) ); - } - else if(4==GlobalV::NSPIN) - { - 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; ipol<4; ++ipol) - v_nspin4(ipol,ir) = vs * chr->rho[ipol][ir] / amag[ir]; - } - } - } - - ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); - return std::make_tuple( etxc, vtxc, std::move(v_nspin4) ); - } // end if(4==GlobalV::NSPIN) - else//NSPIN != 1,2,4 is not supported - { - throw std::domain_error("GlobalV::NSPIN ="+std::to_string(GlobalV::NSPIN) - +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); - } -} - -//the interface to libxc xc_mgga_exc_vxc(xc_func,n,rho,grho,laplrho,tau,e,v1,v2,v3,v4) -//xc_func : LIBXC data type, contains information on xc functional -//n: size of array, nspin*nnr -//rho,grho,laplrho: electron density, its gradient and laplacian -//tau(kin_r): kinetic energy density -//e: energy density -//v1-v4: derivative of energy density w.r.t rho, gradient, laplacian and tau -//v1 and v2 are combined to give v; v4 goes into vofk - -//XC_POLARIZED, XC_UNPOLARIZED: internal flags used in LIBXC, denote the polarized(nspin=1) or unpolarized(nspin=2) calculations, definition can be found in xc.h from LIBXC - -// [etxc, vtxc, v, vofk] = XC_Functional::v_xc(...) -std::tuple XC_Functional::v_xc_meta( - const int &nrxx, // number of real-space grid - const double &omega, // volume of cell - const double tpiba, - const Charge* const chr) -{ - ModuleBase::TITLE("XC_Functional","v_xc_meta"); - ModuleBase::timer::tick("XC_Functional","v_xc_meta"); - - double e2 = 2.0; - - //output of the subroutine - double etxc = 0.0; - double vtxc = 0.0; - ModuleBase::matrix v(GlobalV::NSPIN,nrxx); - ModuleBase::matrix vofk(GlobalV::NSPIN,nrxx); - - //---------------------------------------------------------- - // xc_func_type is defined in Libxc package - // to understand the usage of xc_func_type, - // use can check on website, for example: - // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ - //---------------------------------------------------------- - - const int nspin = GlobalV::NSPIN; - std::vector funcs = init_func( ( (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED ) ); - - // converting rho - std::vector rho; - rho.resize(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]; - } - } - - std::vector>> gdr; - std::vector sigma; - - // calculating grho - gdr.resize( nspin ); - for( int is=0; is!=nspin; ++is ) - { - std::vector rhor(nrxx); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for(int 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); - } - - // converting grho - sigma.resize( nrxx * ((1==nspin)?1:3) ); - - if( 1==nspin ) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for( int ir=0; ir kin_r; - kin_r.resize(nrxx*nspin); -#ifdef _OPENMP -#pragma omp parallel for collapse(2) schedule(static, 1024) -#endif - for( int is=0; iskin_r[is][ir] / 2.0; - } - } - - std::vector exc ( nrxx ); - std::vector vrho ( nrxx * nspin ); - std::vector vsigma ( nrxx * ((1==nspin)?1:3) ); - std::vector vtau ( nrxx * nspin ); - std::vector vlapl ( nrxx * nspin ); - - const double rho_th = 1e-8; - const double grho_th = 1e-12; - const double tau_th = 1e-8; - // sgn for threshold mask - std::vector sgn( nrxx * nspin); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for(int i = 0; i < nrxx * nspin; ++i) - { - sgn[i] = 1.0; - } - - if(nspin == 1) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for( int ir=0; irfamily == XC_FAMILY_MGGA); - xc_mgga_exc_vxc(&func, nrxx, rho.data(), sigma.data(), sigma.data(), - kin_r.data(), exc.data(), vrho.data(), vsigma.data(), vlapl.data(), vtau.data()); - - //process etxc - for( int is=0; is!=nspin; ++is ) - { -#ifdef _OPENMP -#pragma omp parallel for reduction(+:etxc) schedule(static, 256) -#endif - for( int ir=0; ir< nrxx; ++ir ) - { -#ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) - { - exc[ir] *= (1.0 - XC_Functional::hybrid_alpha); - } -#endif - etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is]; - } - } - - //process vtxc -#ifdef _OPENMP -#pragma omp parallel for collapse(2) reduction(+:vtxc) schedule(static, 256) -#endif - for( int is=0; isnumber == XC_MGGA_X_SCAN && get_func_type() == 5) - { - vrho[ir*nspin+is] *= (1.0 - XC_Functional::hybrid_alpha); - } -#endif - const double v_tmp = ModuleBase::e2 * vrho[ir*nspin+is] * sgn[ir*nspin+is]; - v(is,ir) += v_tmp; - vtxc += v_tmp * chr->rho[is][ir]; - } - } - - //process vsigma - std::vector>> h( nspin, std::vector>(nrxx) ); - if( 1==nspin ) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for( int ir=0; ir< nrxx; ++ir ) - { -#ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && get_func_type() == 5) - { - vsigma[ir] *= (1.0 - XC_Functional::hybrid_alpha); - } -#endif - h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir]; - } - } - else - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 64) -#endif - for( int ir=0; ir< nrxx; ++ir ) - { -#ifdef __EXX - if (func.info->number == XC_MGGA_X_SCAN && 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); - } -#endif - h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 - + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); - h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 - + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); - } - } - - // define two dimensional array dh [ nspin, nrxx ] - std::vector> dh(nspin, std::vector( nrxx)); - for( int is=0; is!=nspin; ++is ) - { - XC_Functional::grad_dot( ModuleBase::GlobalFunc::VECTOR_TO_PTR(h[is]), - ModuleBase::GlobalFunc::VECTOR_TO_PTR(dh[is]), chr->rhopw, - tpiba); - } - - double rvtxc = 0.0; -#ifdef _OPENMP -#pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) -#endif - for( int is=0; isnumber == XC_MGGA_X_SCAN && get_func_type() == 5) - { - vtau[ir*nspin+is] *= (1.0 - XC_Functional::hybrid_alpha); - } -#endif - vofk(is,ir) += vtau[ir*nspin+is] * sgn[ir*nspin+is]; - } - } - } - - //------------------------------------------------- - // for MPI, reduce the exchange-correlation energy - //------------------------------------------------- -#ifdef __MPI - Parallel_Reduce::reduce_pool(etxc); - Parallel_Reduce::reduce_pool(vtxc); -#endif - - etxc *= omega / chr->rhopw->nxyz; - vtxc *= omega / chr->rhopw->nxyz; - - finish_func(funcs); - - ModuleBase::timer::tick("XC_Functional","v_xc_meta"); - return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk) ); - -} - -#endif //ifdef USE_LIBXC +} \ No newline at end of file diff --git a/source/module_hamilt_general/module_xc/xc_functional_wrapper_gcxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_wrapper_gcxc.cpp index 09a5800557..4e485f1b54 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_wrapper_gcxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_wrapper_gcxc.cpp @@ -11,6 +11,10 @@ #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_base/global_function.h" +#ifdef USE_LIBXC +#include "xc_functional_libxc.h" +#endif + void XC_Functional::gcxc(const double &rho, const double &grho, double &sxc, double &v1xc, double &v2xc) { @@ -293,96 +297,3 @@ void XC_Functional::gcc_spin(double rho, double &zeta, double grho, double &sc, //} return; } //end subroutine gcc_spin - -void XC_Functional::gcxc_libxc(const double &rho, const double &grho, double &sxc, - double &v1xc, double &v2xc) -{ -#ifdef USE_LIBXC - - const double small = 1.e-6; - const double smallg = 1.e-10; - double s,v1,v2; - sxc = v1xc = v2xc = 0.0; - - if (rho <= small || grho < smallg) - { - return; - } - - std::vector funcs = init_func(XC_UNPOLARIZED); - - - for(xc_func_type &func : funcs) - { - xc_gga_exc_vxc(&func, 1, &rho, &grho, &s, &v1, &v2); - - sxc += s * rho; - v2xc += v2 * 2.0; - v1xc += v1; - } - finish_func(funcs); - - return; -#else - ModuleBase::WARNING_QUIT("gcxc_libxc","compile with LIBXC to use this subroutine"); -#endif -} // end subroutine gcxc_libxc - -void XC_Functional::gcxc_spin_libxc(double rhoup, double rhodw, - ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, - double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud) -{ -#ifdef USE_LIBXC - std::vector funcs = init_func(XC_POLARIZED); - double *rho, *grho, *v1xc, *v2xc, *sgn, s; - sxc = v1xcup = v1xcdw = 0.0; - v2xcup = v2xcdw = v2xcud = 0.0; - rho = new double[2]; - grho= new double[3]; - v1xc= new double[2]; - v2xc= new double[3]; - sgn = new double[2]; - - rho[0] = rhoup; - rho[1] = rhodw; - grho[0] = gdr1.norm2(); - grho[1] = gdr1 * gdr2; - grho[2] = gdr2.norm2(); - - const double rho_threshold = 1E-6; - const double grho_threshold = 1E-10; - - for(xc_func_type &func : funcs) - { - if( func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) - { - sgn[0] = sgn[1] = 1.0; - // call Libxc function: xc_gga_exc_vxc - xc_gga_exc_vxc( &func, 1, rho, grho, &s, v1xc, v2xc); - if(func.info->kind==XC_CORRELATION) - { - if ( rho[0] funcs = init_func(XC_POLARIZED); - - for(xc_func_type &func : funcs) - { - if( func.info->family == XC_FAMILY_LDA) - { - // call Libxc function: xc_lda_exc_vxc - xc_lda_exc_vxc( &func, 1, rho_ud, &e, vxc_ud); - } - exc += e; - vxcup += vxc_ud[0]; - vxcdw += vxc_ud[1]; - } - - finish_func(funcs); - delete[] rho_ud; - delete[] vxc_ud; -#else - ModuleBase::WARNING_QUIT("xc_spin_libxc","compile with LIBXC to use this subroutine"); -#endif } \ No newline at end of file diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp index b4b09c768a..a21ef0a3fd 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp @@ -22,6 +22,9 @@ #include "module_cell/module_paw/paw_cell.h" #endif +#ifdef USE_LIBXC +#include "module_hamilt_general/module_xc/xc_functional_libxc.h" +#endif template @@ -55,7 +58,7 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, { #ifdef USE_LIBXC const auto etxc_vtxc_v - = XC_Functional::v_xc_meta(rho_basis->nrxx, ucell_in.omega, ucell_in.tpiba, chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, ucell_in.omega, ucell_in.tpiba, chr); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp index b92e7a6f3c..376e90eab3 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -4,6 +4,11 @@ #include "module_base/timer.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#ifdef USE_LIBXC +#include "module_hamilt_general/module_xc/xc_functional_libxc.h" +#endif + + //NLCC term, need to be tested template void Stress_Func::stress_cc(ModuleBase::matrix& sigma, @@ -46,7 +51,7 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, { #ifdef USE_LIBXC const auto etxc_vtxc_v - = XC_Functional::v_xc_meta(rho_basis->nrxx, GlobalC::ucell.omega, GlobalC::ucell.tpiba, chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rho_basis->nrxx, GlobalC::ucell.omega, GlobalC::ucell.tpiba, chr); // etxc = std::get<0>(etxc_vtxc_v); // vtxc = std::get<1>(etxc_vtxc_v); diff --git a/source/module_lr/potentials/kernel_xc.cpp b/source/module_lr/potentials/kernel_xc.cpp index f7ab50ebe1..13d18e2a7f 100644 --- a/source/module_lr/potentials/kernel_xc.cpp +++ b/source/module_lr/potentials/kernel_xc.cpp @@ -4,13 +4,17 @@ #include "module_lr/utils/lr_util.h" #ifdef USE_LIBXC #include +#include "module_hamilt_general/module_xc/xc_functional_libxc.h" + void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const double& tpiba, const Charge* chg_gs) { ModuleBase::TITLE("XC_Functional", "f_xc_libxc"); ModuleBase::timer::tick("XC_Functional", "f_xc_libxc"); // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ - std::vector funcs = XC_Functional::init_func((1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED); + std::vector funcs = XC_Functional_Libxc::init_func( + XC_Functional::get_func_id(), + (1 == nspin) ? XC_UNPOLARIZED : XC_POLARIZED); int nrxx = chg_gs->nrxx; // converting rho (extract it as a subfuntion in the future) @@ -216,7 +220,7 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl } } } // end for( xc_func_type &func : funcs ) - XC_Functional::finish_func(funcs); + XC_Functional_Libxc::finish_func(funcs); if (1 == GlobalV::NSPIN || 2 == GlobalV::NSPIN) return; // else if (4 == GlobalV::NSPIN) From 2948a80b4e32020a2fd9a1fde145aded017c2328 Mon Sep 17 00:00:00 2001 From: linpz Date: Mon, 9 Sep 2024 19:34:56 +0800 Subject: [PATCH 2/7] Refactor: move codes in XC_Functional_Libxc::v_xc_libxc() to individual functions --- .../module_xc/CMakeLists.txt | 1 + .../module_xc/xc_functional_libxc.h | 127 +++++-- .../module_xc/xc_functional_libxc_vxc.cpp | 309 +++--------------- 3 files changed, 148 insertions(+), 289 deletions(-) diff --git a/source/module_hamilt_general/module_xc/CMakeLists.txt b/source/module_hamilt_general/module_xc/CMakeLists.txt index 6f179ba2e2..4db1ef083d 100644 --- a/source/module_hamilt_general/module_xc/CMakeLists.txt +++ b/source/module_hamilt_general/module_xc/CMakeLists.txt @@ -12,6 +12,7 @@ add_library( xc_funct_corr_gga.cpp xc_funct_hcth.cpp xc_functional_libxc.cpp + xc_functional_libxc_tools.cpp xc_functional_libxc_vxc.cpp xc_functional_libxc_wrapper_xc.cpp xc_functional_libxc_wrapper_gcxc.cpp 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 bbe5e416fc..f4ee8361e4 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc.h +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.h @@ -17,42 +17,121 @@ namespace XC_Functional_Libxc { //------------------- // xc_functional_libxc.cpp -//------------------- +//------------------- // sets functional type, which allows combination of LIBXC keyword connected by "+" // for example, "XC_LDA_X+XC_LDA_C_PZ" - std::pair> set_xc_type_libxc(const std::string xc_func_in); + extern std::pair> set_xc_type_libxc(const std::string xc_func_in); // converts func_id into corresponding xc_func_type vector - std::vector init_func(const std::vector &func_id, const int xc_polarized); + extern std::vector init_func(const std::vector &func_id, const int xc_polarized); + + extern void finish_func(std::vector &funcs); - void finish_func(std::vector &funcs); //------------------- // xc_functional_libxc_vxc.cpp -//------------------- +//------------------- - std::tuple v_xc_libxc( - const std::vector &func_id, + 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 // for mGGA functional - std::tuple v_xc_meta( - const std::vector &func_id, + extern std::tuple v_xc_meta( + 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); + +//------------------- +// xc_functional_libxc_tools.cpp +//------------------- + + // converting rho (abacus=>libxc) + extern std::vector convert_rho( + const int nspin, + const std::size_t nrxx, + const Charge* const chr); + + // converting rho (abacus=>libxc) + extern std::tuple, std::vector> convert_rho_amag_nspin4( + const int nspin, + const std::size_t nrxx, + const Charge* const chr); + + // calculating grho + extern std::vector>> cal_gdr( + const int nspin, + const std::vector &rho, + const double tpiba, + const Charge* const chr); + + // converting grho (abacus=>libxc) + extern std::vector convert_sigma( + const std::vector>> &gdr); + + // sgn for threshold mask + extern std::vector 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); + + // converting etxc from exc (libxc=>abacus) + extern double convert_etxc( + const int nspin, + const std::size_t nrxx, + const std::vector &sgn, + const std::vector &rho, + std::vector exc); + + // converting vtxc and v from vrho and vsigma (libxc=>abacus) + extern double 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); + + // dh for gga v + extern std::vector> 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); + + // convert v for NSPIN=4 + extern ModuleBase::matrix convert_v_nspin4( + const std::size_t nrxx, + const Charge* const chr, + const std::vector &amag, + const ModuleBase::matrix &v); + + //------------------- // xc_functional_libxc_wrapper_xc.cpp //------------------- - void xc_spin_libxc( - const std::vector &func_id, + extern void xc_spin_libxc( + const std::vector &func_id, const double &rhoup, const double &rhodw, double &exc, double &vxcup, double &vxcdw); @@ -62,15 +141,15 @@ namespace XC_Functional_Libxc //------------------- // the entire GGA functional, for nspin=1 case - void gcxc_libxc( - const std::vector &func_id, + extern void gcxc_libxc( + const std::vector &func_id, const double &rho, const double &grho, double &sxc, double &v1xc, double &v2xc); // the entire GGA functional, for nspin=2 case - void gcxc_spin_libxc( - const std::vector &func_id, - double rhoup, double rhodw, + extern void gcxc_spin_libxc( + const std::vector &func_id, + double rhoup, double rhodw, ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud); @@ -80,20 +159,20 @@ namespace XC_Functional_Libxc //------------------- // wrapper for the mGGA functionals - void tau_xc( - const std::vector &func_id, + extern void tau_xc( + const std::vector &func_id, const double &rho, const double &grho, const double &atau, double &sxc, double &v1xc, double &v2xc, double &v3xc); - void tau_xc_spin( - const std::vector &func_id, - double rhoup, double rhodw, + extern void tau_xc_spin( + const std::vector &func_id, + double rhoup, double rhodw, ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, double tauup, double taudw, double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, double &v3xcup, double &v3xcdw); -} +} // namespace XC_Functional_Libxc -#endif +#endif // USE_LIBXC -#endif \ No newline at end of file +#endif // XC_FUNCTIONAL_LIBXC_H \ 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 7c355e2c00..9e8614dfa6 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 @@ -22,14 +22,10 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / ModuleBase::TITLE("XC_Functional","v_xc_libxc"); ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); - const int nspin = + const int nspin = (GlobalV::NSPIN == 1 || ( GlobalV::NSPIN ==4 && !GlobalV::DOMAG && !GlobalV::DOMAG_Z)) ? 1 : 2; - double etxc = 0.0; - double vtxc = 0.0; - ModuleBase::matrix v(nspin,nrxx); - //---------------------------------------------------------- // xc_func_type is defined in Libxc package // to understand the usage of xc_func_type, @@ -54,85 +50,30 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / }(); // converting rho - std::vector rho(nrxx*nspin); + std::vector rho; std::vector amag; if(1==nspin || 2==GlobalV::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]; + rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); } else { - amag.resize(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],std::vector> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr); + rho = std::get<0>(std::move(rho_amag)); + amag = std::get<1>(std::move(rho_amag)); } std::vector>> gdr; std::vector sigma; if(is_gga) { - // calculating grho - gdr.resize( nspin ); - for( int is=0; is!=nspin; ++is ) - { - std::vector rhor(nrxx); - #ifdef _OPENMP - #pragma omp parallel for schedule(static, 1024) - #endif - for(int 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) - - // converting grho - sigma.resize( nrxx * ((1==nspin)?1:3) ); - if( 1==nspin ) - { - #ifdef _OPENMP - #pragma omp parallel for schedule(static, 1024) - #endif - for( int ir=0; ir XC_Functional_Libxc::v_xc_libxc( / xc_func_set_dens_threshold(&func, rho_threshold); // sgn for threshold mask - 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; ir sgn = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma); std::vector exc ( nrxx ); std::vector vrho ( nrxx * nspin ); @@ -167,7 +93,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / { case XC_FAMILY_LDA: // call Libxc function: xc_lda_exc_vxc - xc_lda_exc_vxc( &func, nrxx, rho.data(), + xc_lda_exc_vxc( &func, nrxx, rho.data(), exc.data(), vrho.data() ); break; case XC_FAMILY_GGA: @@ -182,74 +108,20 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / break; } - #ifdef _OPENMP - #pragma omp parallel for collapse(2) reduction(+:etxc) schedule(static, 256) - #endif - for( int is=0; isfamily == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) - { - std::vector>> h( nspin, - std::vector>(nrxx) ); - if( 1==nspin ) - { - #ifdef _OPENMP - #pragma omp parallel for schedule(static, 1024) - #endif - for( int ir=0; ir< nrxx; ++ir ) - h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir]; - } - else - { - #ifdef _OPENMP - #pragma omp parallel for schedule(static, 1024) - #endif - for( int ir=0; ir< nrxx; ++ir ) - { - h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 - + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); - h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 - + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); - } - } - - // define two dimensional array dh [ nspin, nrxx ] - std::vector> 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); - - 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)) + etxc += XC_Functional_Libxc::convert_etxc(nspin, nrxx, sgn, rho, exc); + vtxc += XC_Functional_Libxc::convert_vtxc_v( + func, nspin, nrxx, + sgn, rho, gdr, + vrho, vsigma, + tpiba, chr, + v); } // end for( xc_func_type &func : funcs ) + if(4==GlobalV::NSPIN) + { + v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v); + } + //------------------------------------------------- // for MPI, reduce the exchange-correlation energy //------------------------------------------------- @@ -261,40 +133,10 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / etxc *= omega / chr->rhopw->nxyz; vtxc *= omega / chr->rhopw->nxyz; - finish_func(funcs); - - if(1==GlobalV::NSPIN || 2==GlobalV::NSPIN) - { - ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); - return std::make_tuple( etxc, vtxc, std::move(v) ); - } - else if(4==GlobalV::NSPIN) - { - 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; ipol<4; ++ipol) - v_nspin4(ipol,ir) = vs * chr->rho[ipol][ir] / amag[ir]; - } - } - } + XC_Functional_Libxc::finish_func(funcs); - ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); - return std::make_tuple( etxc, vtxc, std::move(v_nspin4) ); - } // end if(4==GlobalV::NSPIN) - else//NSPIN != 1,2,4 is not supported - { - throw std::domain_error("GlobalV::NSPIN ="+std::to_string(GlobalV::NSPIN) - +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); - } + ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); + return std::make_tuple( etxc, vtxc, std::move(v) ); } @@ -334,79 +176,14 @@ std::tuple XC_Functional_Li // use can check on website, for example: // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ //---------------------------------------------------------- - + const int nspin = GlobalV::NSPIN; std::vector funcs = XC_Functional_Libxc::init_func( func_id, ( (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED ) ); - // converting rho - std::vector rho; - rho.resize(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]; - } - } - - std::vector>> gdr; - std::vector sigma; - - // calculating grho - gdr.resize( nspin ); - for( int is=0; is!=nspin; ++is ) - { - std::vector rhor(nrxx); -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for(int 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); - } - - // converting grho - sigma.resize( nrxx * ((1==nspin)?1:3) ); - - if( 1==nspin ) - { -#ifdef _OPENMP -#pragma omp parallel for schedule(static, 1024) -#endif - for( int ir=0; ir rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); + const std::vector>> gdr + = XC_Functional_Libxc::cal_gdr(nspin, rho, tpiba, chr); + const std::vector sigma = XC_Functional_Libxc::convert_sigma(gdr); //converting kin_r std::vector kin_r; @@ -489,7 +266,7 @@ std::tuple XC_Functional_Li } #endif etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is]; - } + } } //process vtxc @@ -513,7 +290,9 @@ std::tuple XC_Functional_Li } //process vsigma - std::vector>> h( nspin, std::vector>(nrxx) ); + std::vector>> h( + nspin, + std::vector>(nrxx) ); if( 1==nspin ) { #ifdef _OPENMP @@ -545,9 +324,9 @@ std::tuple XC_Functional_Li vsigma[ir*3+2] *= (1.0 - XC_Functional::hybrid_alpha); } #endif - h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 + h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); - h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 + h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); } } @@ -556,8 +335,8 @@ std::tuple XC_Functional_Li std::vector> dh(nspin, std::vector( nrxx)); for( int is=0; is!=nspin; ++is ) { - XC_Functional::grad_dot( ModuleBase::GlobalFunc::VECTOR_TO_PTR(h[is]), - ModuleBase::GlobalFunc::VECTOR_TO_PTR(dh[is]), chr->rhopw, + XC_Functional::grad_dot( h[is].data(), + dh[is].data(), chr->rhopw, tpiba); } @@ -593,7 +372,7 @@ std::tuple XC_Functional_Li } } } - + //------------------------------------------------- // for MPI, reduce the exchange-correlation energy //------------------------------------------------- From 913d2e3db06f125899dfbc4dbf7589b999e44a35 Mon Sep 17 00:00:00 2001 From: linpz Date: Mon, 9 Sep 2024 20:09:27 +0800 Subject: [PATCH 3/7] 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; From 35d90df236aea32223d8c456bc8ef9a14cd9db2e Mon Sep 17 00:00:00 2001 From: linpz Date: Mon, 9 Sep 2024 20:31:17 +0800 Subject: [PATCH 4/7] Refactor little codes in XC_Functional_Libxc --- .../module_xc/test/CMakeLists.txt | 10 +--- .../module_xc/test/test_xc.cpp | 5 +- .../module_xc/test/test_xc1.cpp | 4 +- .../module_xc/test/test_xc2.cpp | 5 +- .../module_xc/test/test_xc4.cpp | 3 +- .../module_xc/test/test_xc5.cpp | 5 +- .../module_xc/xc_functional.cpp | 8 +-- .../module_xc/xc_functional_libxc.cpp | 16 ++--- .../module_xc/xc_functional_libxc.h | 4 +- .../module_xc/xc_functional_libxc_vxc.cpp | 6 +- .../xc_functional_libxc_wrapper_gcxc.cpp | 58 +++++++------------ .../xc_functional_libxc_wrapper_tauxc.cpp | 52 +++++------------ .../xc_functional_libxc_wrapper_xc.cpp | 18 ++---- 13 files changed, 76 insertions(+), 118 deletions(-) diff --git a/source/module_hamilt_general/module_xc/test/CMakeLists.txt b/source/module_hamilt_general/module_xc/test/CMakeLists.txt index ee4b4c5869..858c68d188 100644 --- a/source/module_hamilt_general/module_xc/test/CMakeLists.txt +++ b/source/module_hamilt_general/module_xc/test/CMakeLists.txt @@ -9,7 +9,7 @@ AddTest( AddTest( TARGET XCTest_HSE LIBS MPI::MPI_CXX Libxc::xc parameter # required by global.h; for details, `remove_definitions(-D__MPI)`. - SOURCES test_xc1.cpp ../xc_functional.cpp + SOURCES test_xc1.cpp ../xc_functional.cpp ../xc_functional_libxc.cpp ) @@ -26,8 +26,6 @@ AddTest( SOURCES test_xc3.cpp ../xc_functional_gradcorr.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp ../xc_functional_libxc.cpp - ../xc_functional_libxc_tools.cpp - ../xc_functional_libxc_vxc.cpp ../xc_functional_libxc_wrapper_xc.cpp ../xc_functional_libxc_wrapper_gcxc.cpp ../xc_functional_libxc_wrapper_tauxc.cpp @@ -45,8 +43,6 @@ AddTest( SOURCES test_xc4.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp ../xc_functional_libxc.cpp - ../xc_functional_libxc_tools.cpp - ../xc_functional_libxc_vxc.cpp ../xc_functional_libxc_wrapper_xc.cpp ../xc_functional_libxc_wrapper_gcxc.cpp ../xc_functional_libxc_wrapper_tauxc.cpp @@ -60,14 +56,14 @@ AddTest( SOURCES test_xc5.cpp ../xc_functional_gradcorr.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp ../xc_functional_libxc.cpp - ../xc_functional_libxc_tools.cpp - ../xc_functional_libxc_vxc.cpp ../xc_functional_libxc_wrapper_xc.cpp ../xc_functional_libxc_wrapper_gcxc.cpp ../xc_functional_libxc_wrapper_tauxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp ../xc_functional_vxc.cpp + ../xc_functional_libxc_vxc.cpp + ../xc_functional_libxc_tools.cpp ../../../module_base/matrix.cpp ../../../module_base/memory.cpp ../../../module_base/timer.cpp diff --git a/source/module_hamilt_general/module_xc/test/test_xc.cpp b/source/module_hamilt_general/module_xc/test/test_xc.cpp index 654b20ca95..c257afaf86 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc.cpp @@ -1,5 +1,6 @@ #include "gtest/gtest.h" #include "../xc_functional.h" +#include "../xc_functional_libxc.h" #include "../exx_info.h" #include "xctest.h" @@ -829,7 +830,7 @@ class XCTest_PBE0 : public XCTest void SetUp() { XC_Functional::set_xc_type("PBE0"); - XC_Functional::get_hybrid_alpha(0.5); + XC_Functional::set_hybrid_alpha(0.5); std::vector rho = {0.17E+01, 0.17E+01, 0.15E+01, 0.88E-01, 0.18E+04}; std::vector grho = {0.81E-11, 0.17E+01, 0.36E+02, 0.87E-01, 0.55E+00}; @@ -884,7 +885,7 @@ class XCTest_PBE_LibXC : public XCTest XC_Functional::xc(rho[i],e,v); e_lda.push_back(e); v_lda.push_back(v); - XC_Functional::gcxc_libxc(rho[i],grho[i],e,v1,v2); + XC_Functional_Libxc::gcxc_libxc(XC_Functional::get_func_id(), rho[i],grho[i],e,v1,v2); e_gga.push_back(e); v1_gga.push_back(v1); v2_gga.push_back(v2); diff --git a/source/module_hamilt_general/module_xc/test/test_xc1.cpp b/source/module_hamilt_general/module_xc/test/test_xc1.cpp index 3d0105aa73..9a5c45eebd 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc1.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc1.cpp @@ -38,7 +38,7 @@ class XCTest_HSE : public XCTest void SetUp() { XC_Functional::set_xc_type("HSE"); - XC_Functional::get_hybrid_alpha(0.5); + XC_Functional::set_hybrid_alpha(0.5); } }; @@ -53,7 +53,7 @@ class XCTest_SCAN0 : public XCTest void SetUp() { XC_Functional::set_xc_type("SCAN0"); - XC_Functional::get_hybrid_alpha(0.5); + XC_Functional::set_hybrid_alpha(0.5); } }; diff --git a/source/module_hamilt_general/module_xc/test/test_xc2.cpp b/source/module_hamilt_general/module_xc/test/test_xc2.cpp index be91b74aa5..ad23626ac3 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc2.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc2.cpp @@ -1,6 +1,7 @@ #include "gtest/gtest.h" #include "xctest.h" #include "../xc_functional.h" +#include "../xc_functional_libxc.h" #include "../exx_info.h" /************************************************ * unit test of functionals @@ -455,7 +456,7 @@ class XCTest_PBE_SPN_LibXC : public XCTest double e,v1,v2,v3,v4,v5; double r1 = rho[i] * (1+zeta[i]) / 2.0; double r2 = rho[i] * (1-zeta[i]) / 2.0; - XC_Functional::gcxc_spin_libxc(r1,r2,gdr[i],gdr[i],e,v1,v2,v3,v4,v5); + XC_Functional_Libxc::gcxc_spin_libxc(XC_Functional::get_func_id(), r1,r2,gdr[i],gdr[i],e,v1,v2,v3,v4,v5); e_gga.push_back(e); v1_gga.push_back(v1+v3); v2_gga.push_back(v2+v4); @@ -494,7 +495,7 @@ class XCTest_PZ_SPN_LibXC : public XCTest double e,v1,v2; double r1 = rho[i] * (1+zeta[i]) / 2.0; double r2 = rho[i] * (1-zeta[i]) / 2.0; - XC_Functional::xc_spin_libxc(r1,r2,e,v1,v2); + XC_Functional_Libxc::xc_spin_libxc(XC_Functional::get_func_id(), r1,r2,e,v1,v2); e_lda.push_back(e); v1_lda.push_back(v1); v2_lda.push_back(v2); diff --git a/source/module_hamilt_general/module_xc/test/test_xc4.cpp b/source/module_hamilt_general/module_xc/test/test_xc4.cpp index eb179cb0a2..41fd60f425 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc4.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc4.cpp @@ -1,4 +1,5 @@ #include "../xc_functional.h" +#include "../xc_functional_libxc.h" #include "gtest/gtest.h" #include "xctest.h" #include "../exx_info.h" @@ -43,7 +44,7 @@ class XCTest_SCAN : public XCTest for(int i=0;i<5;i++) { double e,v,v1,v2,v3; - XC_Functional::tau_xc(rho[i],grho[i],tau[i],e,v1,v2,v3); + XC_Functional_Libxc::tau_xc(XC_Functional::get_func_id(), rho[i],grho[i],tau[i],e,v1,v2,v3); e_.push_back(e); v1_.push_back(v1); v2_.push_back(v2); diff --git a/source/module_hamilt_general/module_xc/test/test_xc5.cpp b/source/module_hamilt_general/module_xc/test/test_xc5.cpp index adfd82ba28..1b04d6ad35 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc5.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc5.cpp @@ -1,4 +1,5 @@ #include "../xc_functional.h" +#include "../xc_functional_libxc.h" #include "gtest/gtest.h" #include "xctest.h" #include "../exx_info.h" @@ -279,7 +280,7 @@ class XCTest_VXC_meta : public XCTest GlobalV::NSPIN = 1; std::tuple etxc_vtxc_v - = XC_Functional::v_xc_meta(rhopw.nrxx,ucell.omega,ucell.tpiba,&chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr); et1 = std::get<0>(etxc_vtxc_v); vt1 = std::get<1>(etxc_vtxc_v); v1 = std::get<2>(etxc_vtxc_v); @@ -287,7 +288,7 @@ class XCTest_VXC_meta : public XCTest GlobalV::NSPIN = 2; etxc_vtxc_v - = XC_Functional::v_xc_meta(rhopw.nrxx,ucell.omega,ucell.tpiba,&chr); + = XC_Functional_Libxc::v_xc_meta(XC_Functional::get_func_id(), rhopw.nrxx,ucell.omega,ucell.tpiba,&chr); et2 = std::get<0>(etxc_vtxc_v); vt2 = std::get<1>(etxc_vtxc_v); v2 = std::get<2>(etxc_vtxc_v); diff --git a/source/module_hamilt_general/module_xc/xc_functional.cpp b/source/module_hamilt_general/module_xc/xc_functional.cpp index c9730e957b..1dd6758a98 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional.cpp @@ -58,7 +58,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) //such as for PBE we have: // func_id.push_back(XC_GGA_X_PBE); // func_id.push_back(XC_GGA_C_PBE); - + func_id.clear(); std::string xc_func = xc_func_in; std::transform(xc_func.begin(), xc_func.end(), xc_func.begin(), (::toupper)); @@ -128,7 +128,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_id.push_back(XC_GGA_C_PBE); func_type = 2; use_libxc = false; - } + } else if ( xc_func == "BLYP") //B88+LYP { func_id.push_back(XC_GGA_X_B88); @@ -142,14 +142,14 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_id.push_back(XC_GGA_C_P86); func_type = 2; use_libxc = false; - } + } else if ( xc_func == "PW91") //PW91_X+PW91_C { func_id.push_back(XC_GGA_X_PW91); func_id.push_back(XC_GGA_C_PW91); func_type = 2; use_libxc = false; - } + } else if ( xc_func == "HCTH") //HCTH_X+HCTH_C { func_id.push_back(XC_GGA_X_HCTH_A); 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 8d80f27714..47ddd1ed1f 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp @@ -50,7 +50,7 @@ std::vector XC_Functional_Libxc::init_func(const std::vector std::vector funcs; //------------------------------------------- - // define a function named 'add_func', which + // define a function named 'add_func', which // initialize a functional according to its ID //------------------------------------------- auto add_func = [&]( const int func_id ) @@ -71,17 +71,17 @@ std::vector XC_Functional_Libxc::init_func(const std::vector #ifdef __EXX else if( id == XC_HYB_GGA_XC_PBEH ) // PBE0 { - add_func( XC_HYB_GGA_XC_PBEH ); - double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, - GlobalC::exx_info.info_global.hse_omega, + add_func( XC_HYB_GGA_XC_PBEH ); + double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, + GlobalC::exx_info.info_global.hse_omega, GlobalC::exx_info.info_global.hse_omega }; - xc_func_set_ext_params(&funcs.back(), parameter_hse); + xc_func_set_ext_params(&funcs.back(), parameter_hse); } else if( id == XC_HYB_GGA_XC_HSE06 ) // HSE06 hybrid functional { - add_func( XC_HYB_GGA_XC_HSE06 ); - double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, - GlobalC::exx_info.info_global.hse_omega, + add_func( XC_HYB_GGA_XC_HSE06 ); + double parameter_hse[3] = { GlobalC::exx_info.info_global.hybrid_alpha, + GlobalC::exx_info.info_global.hse_omega, GlobalC::exx_info.info_global.hse_omega }; xc_func_set_ext_params(&funcs.back(), parameter_hse); } 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 f4ee8361e4..73bf03e51c 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc.h +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.h @@ -149,8 +149,8 @@ namespace XC_Functional_Libxc // the entire GGA functional, for nspin=2 case extern void gcxc_spin_libxc( const std::vector &func_id, - double rhoup, double rhodw, - ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + const double rhoup, const double rhodw, + const ModuleBase::Vector3 gdr1, const ModuleBase::Vector3 gdr2, double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud); 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 e8b80537f1..1a548fb244 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 @@ -205,9 +205,9 @@ std::tuple XC_Functional_Li std::vector vtau ( nrxx * nspin ); std::vector vlapl ( nrxx * nspin ); - const double rho_th = 1e-8; - const double grho_th = 1e-12; - const double tau_th = 1e-8; + constexpr double rho_th = 1e-8; + constexpr double grho_th = 1e-12; + constexpr double tau_th = 1e-8; // sgn for threshold mask std::vector sgn( nrxx * nspin); #ifdef _OPENMP diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp index 24885a5e37..f46efc8e70 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp @@ -5,66 +5,52 @@ #include void XC_Functional_Libxc::gcxc_libxc( - const std::vector &func_id, - const double &rho, const double &grho, double &sxc, - double &v1xc, double &v2xc) + const std::vector &func_id, + const double &rho, const double &grho, + double &sxc, double &v1xc, double &v2xc) { - const double small = 1.e-6; - const double smallg = 1.e-10; - double s,v1,v2; sxc = v1xc = v2xc = 0.0; + constexpr double small = 1.e-6; + constexpr double smallg = 1.e-10; if (rho <= small || grho < smallg) { return; } std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED); - - for(xc_func_type &func : funcs) { + double s,v1,v2; xc_gga_exc_vxc(&func, 1, &rho, &grho, &s, &v1, &v2); - sxc += s * rho; - v2xc += v2 * 2.0; v1xc += v1; + v2xc += v2 * 2.0; } XC_Functional_Libxc::finish_func(funcs); } // end subroutine gcxc_libxc + + void XC_Functional_Libxc::gcxc_spin_libxc( const std::vector &func_id, - double rhoup, double rhodw, - ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, + const double rhoup, const double rhodw, + const ModuleBase::Vector3 gdr1, const ModuleBase::Vector3 gdr2, double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud) { - std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); - double *rho, *grho, *v1xc, *v2xc, *sgn, s; sxc = v1xcup = v1xcdw = 0.0; v2xcup = v2xcdw = v2xcud = 0.0; - rho = new double[2]; - grho= new double[3]; - v1xc= new double[2]; - v2xc= new double[3]; - sgn = new double[2]; - - rho[0] = rhoup; - rho[1] = rhodw; - grho[0] = gdr1.norm2(); - grho[1] = gdr1 * gdr2; - grho[2] = gdr2.norm2(); - - const double rho_threshold = 1E-6; - const double grho_threshold = 1E-10; + const std::vector rho = {rhoup, rhodw}; + const std::vector grho = {gdr1.norm2(), gdr1*gdr2, gdr2.norm2()}; + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); for(xc_func_type &func : funcs) { if( func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) { - sgn[0] = sgn[1] = 1.0; - // call Libxc function: xc_gga_exc_vxc - xc_gga_exc_vxc( &func, 1, rho, grho, &s, v1xc, v2xc); + constexpr double rho_threshold = 1E-6; + constexpr double grho_threshold = 1E-10; + std::vector sgn = {1.0, 1.0}; if(func.info->kind==XC_CORRELATION) { if ( rho[0] v1xc(2), v2xc(3); + // call Libxc function: xc_gga_exc_vxc + xc_gga_exc_vxc( &func, 1, rho.data(), grho.data(), &s, v1xc.data(), v2xc.data()); sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]); v1xcup += v1xc[0] * sgn[0]; v1xcdw += v1xc[1] * sgn[1]; @@ -80,11 +71,6 @@ void XC_Functional_Libxc::gcxc_spin_libxc( v2xcdw += 2.0 * v2xc[2] * sgn[1]; } } - delete[] grho; - delete[] rho; - delete[] v1xc; - delete[] v2xc; - delete[] sgn; XC_Functional_Libxc::finish_func(funcs); } 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 48bd301ae6..a5d21d49b7 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 @@ -18,7 +18,7 @@ void XC_Functional_Libxc::tau_xc( double lapl_rho, vlapl_rho; lapl_rho = grho; std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_UNPOLARIZED); - + sxc = 0.0; v1xc = 0.0; v2xc = 0.0; v3xc = 0.0; for(xc_func_type &func : funcs) @@ -46,48 +46,29 @@ void XC_Functional_Libxc::tau_xc( void XC_Functional_Libxc::tau_xc_spin( const std::vector &func_id, - double rhoup, double rhodw, + double rhoup, double rhodw, ModuleBase::Vector3 gdr1, ModuleBase::Vector3 gdr2, double tauup, double taudw, double &sxc, double &v1xcup, double &v1xcdw, double &v2xcup, double &v2xcdw, double &v2xcud, double &v3xcup, double &v3xcdw) { - - std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); - - double *rho, *grho, *tau, *v1xc, *v2xc, *v3xc, *sgn, s; sxc = v1xcup = v1xcdw = 0.0; v2xcup = v2xcdw = v2xcud = 0.0; v3xcup = v3xcdw = 0.0; - rho = new double[2]; - grho= new double[3]; - tau = new double[2]; - v1xc= new double[2]; - v2xc= new double[3]; - v3xc= new double[2]; - sgn = new double[2]; - - rho[0] = rhoup; - rho[1] = rhodw; - grho[0] = gdr1.norm2(); - grho[1] = gdr1 * gdr2; - grho[2] = gdr2.norm2(); - tau[0] = tauup; - tau[1] = taudw; - double lapl[2] = {0.0,0.0}; - double vlapl[2]; + const std::vector rho = {rhoup, rhodw}; + const std::vector grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; + const std::vector tau = {tauup, taudw}; - const double rho_threshold = 1E-6; - const double grho_threshold = 1E-10; + std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); for(xc_func_type &func : funcs) { if( func.info->family == XC_FAMILY_MGGA || func.info->family == XC_FAMILY_HYB_MGGA) { - sgn[0] = sgn[1] = 1.0; - // call Libxc function: xc_mgga_exc_vxc - xc_mgga_exc_vxc( &func, 1, rho, grho, lapl, tau, &s, v1xc, v2xc, vlapl, v3xc); + constexpr double rho_threshold = 1E-6; + constexpr double grho_threshold = 1E-10; + std::vector sgn = {1.0, 1.0}; if(func.info->kind==XC_CORRELATION) { if ( rho[0] v1xc(2), v2xc(3), v3xc(2), lapl(2), vlapl(2); + // call Libxc function: xc_mgga_exc_vxc + xc_mgga_exc_vxc( &func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s, v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data()); + sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]); v1xcup += v1xc[0] * sgn[0]; v1xcdw += v1xc[1] * sgn[1]; @@ -106,16 +93,7 @@ void XC_Functional_Libxc::tau_xc_spin( } } - delete[] grho; - delete[] rho; - delete[] tau; - delete[] v1xc; - delete[] v2xc; - delete[] v3xc; - delete[] sgn; XC_Functional_Libxc::finish_func(funcs); - - return; } #endif \ No newline at end of file diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp index c5693ec54d..aa2cf664ad 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp @@ -5,34 +5,28 @@ void XC_Functional_Libxc::xc_spin_libxc( const std::vector &func_id, const double &rhoup, const double &rhodw, - double &exc, double &vxcup, double &vxcdw) + double &exc, double &vxcup, double &vxcdw) { - double e, vup, vdw; - double *rho_ud, *vxc_ud; + const std::vector rho_ud = {rhoup, rhodw}; exc = vxcup = vxcdw = 0.0; - rho_ud = new double[2]; - vxc_ud = new double[2]; - rho_ud[0] = rhoup; - rho_ud[1] = rhodw; - std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); for(xc_func_type &func : funcs) { + double e = 0.0; + std::vector vxc_ud(2); if( func.info->family == XC_FAMILY_LDA) { // call Libxc function: xc_lda_exc_vxc - xc_lda_exc_vxc( &func, 1, rho_ud, &e, vxc_ud); + xc_lda_exc_vxc( &func, 1, rho_ud.data(), &e, vxc_ud.data()); } exc += e; vxcup += vxc_ud[0]; vxcdw += vxc_ud[1]; - } + } XC_Functional_Libxc::finish_func(funcs); - delete[] rho_ud; - delete[] vxc_ud; } #endif \ No newline at end of file From 2d49c3db8a97e62f072e75400092a210b084ef4f Mon Sep 17 00:00:00 2001 From: LUNASEA <33978601+maki49@users.noreply.github.com> Date: Tue, 10 Sep 2024 20:40:10 +0800 Subject: [PATCH 5/7] fix compile with UT (#4) --- source/Makefile.Objects | 7 ++++++- source/module_hamilt_general/module_xc/test/CMakeLists.txt | 2 ++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 50d3c7f842..533b496856 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -454,7 +454,12 @@ OBJS_XC=xc_functional.o\ xc_functional_gradcorr.o\ xc_functional_wrapper_xc.o\ xc_functional_wrapper_gcxc.o\ - xc_functional_wrapper_tauxc.o\ + xc_functional_libxc.o\ + xc_functional_libxc_tools.o\ + xc_functional_libxc_vxc.o\ + xc_functional_libxc_wrapper_xc.o\ + xc_functional_libxc_wrapper_gcxc.o\ + xc_functional_libxc_wrapper_tauxc.o\ xc_funct_exch_lda.o\ xc_funct_corr_lda.o\ xc_funct_exch_gga.o\ diff --git a/source/module_hamilt_general/module_xc/test/CMakeLists.txt b/source/module_hamilt_general/module_xc/test/CMakeLists.txt index 858c68d188..4d879e6e85 100644 --- a/source/module_hamilt_general/module_xc/test/CMakeLists.txt +++ b/source/module_hamilt_general/module_xc/test/CMakeLists.txt @@ -4,6 +4,7 @@ AddTest( TARGET XCTest_PBE LIBS MPI::MPI_CXX Libxc::xc parameter # required by global.h; for details, `remove_definitions(-D__MPI)`. SOURCES test_xc.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp + ../xc_functional_libxc_wrapper_gcxc.cpp ../xc_functional_libxc.cpp ) AddTest( @@ -17,6 +18,7 @@ AddTest( TARGET XCTest_PZ_SPN LIBS MPI::MPI_CXX Libxc::xc parameter # required by global.h; for details, `remove_definitions(-D__MPI)`. SOURCES test_xc2.cpp ../xc_functional.cpp ../xc_functional_wrapper_xc.cpp ../xc_functional_wrapper_gcxc.cpp ../xc_funct_corr_gga.cpp ../xc_funct_corr_lda.cpp ../xc_funct_exch_gga.cpp ../xc_funct_exch_lda.cpp ../xc_funct_hcth.cpp + ../xc_functional_libxc_wrapper_gcxc.cpp ../xc_functional_libxc_wrapper_xc.cpp ../xc_functional_libxc.cpp ) From 31cdc4cb357fae2e2fa1dccadc6a380952a4d824 Mon Sep 17 00:00:00 2001 From: linpz Date: Fri, 20 Sep 2024 22:05:07 +0800 Subject: [PATCH 6/7] change vector to array in xc_functional_libxc --- .../module_xc/xc_functional_libxc_wrapper_gcxc.cpp | 10 ++++++---- .../module_xc/xc_functional_libxc_wrapper_tauxc.cpp | 12 +++++++----- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp index f46efc8e70..ce4aad00d8 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp @@ -3,6 +3,7 @@ #include "xc_functional_libxc.h" #include +#include void XC_Functional_Libxc::gcxc_libxc( const std::vector &func_id, @@ -40,8 +41,8 @@ void XC_Functional_Libxc::gcxc_spin_libxc( { sxc = v1xcup = v1xcdw = 0.0; v2xcup = v2xcdw = v2xcud = 0.0; - const std::vector rho = {rhoup, rhodw}; - const std::vector grho = {gdr1.norm2(), gdr1*gdr2, gdr2.norm2()}; + const std::array rho = {rhoup, rhodw}; + const std::array grho = {gdr1.norm2(), gdr1*gdr2, gdr2.norm2()}; std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); for(xc_func_type &func : funcs) @@ -50,7 +51,7 @@ void XC_Functional_Libxc::gcxc_spin_libxc( { constexpr double rho_threshold = 1E-6; constexpr double grho_threshold = 1E-10; - std::vector sgn = {1.0, 1.0}; + std::array sgn = {1.0, 1.0}; if(func.info->kind==XC_CORRELATION) { if ( rho[0] v1xc(2), v2xc(3); + std::array v1xc; + std::array v2xc; // call Libxc function: xc_gga_exc_vxc xc_gga_exc_vxc( &func, 1, rho.data(), grho.data(), &s, v1xc.data(), v2xc.data()); sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]); 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 a5d21d49b7..b265af9bb1 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 @@ -6,6 +6,7 @@ #include "xc_functional_libxc.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include //tau_xc and tau_xc_spin: interface for calling xc_mgga_exc_vxc from LIBXC //XC_POLARIZED, XC_UNPOLARIZED: internal flags used in LIBXC, denote the polarized(nspin=1) or unpolarized(nspin=2) calculations, definition can be found in xc.h from LIBXC @@ -56,9 +57,9 @@ void XC_Functional_Libxc::tau_xc_spin( v2xcup = v2xcdw = v2xcud = 0.0; v3xcup = v3xcdw = 0.0; - const std::vector rho = {rhoup, rhodw}; - const std::vector grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; - const std::vector tau = {tauup, taudw}; + const std::array rho = {rhoup, rhodw}; + const std::array grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; + const std::array tau = {tauup, taudw}; std::vector funcs = XC_Functional_Libxc::init_func(func_id, XC_POLARIZED); @@ -68,7 +69,7 @@ void XC_Functional_Libxc::tau_xc_spin( { constexpr double rho_threshold = 1E-6; constexpr double grho_threshold = 1E-10; - std::vector sgn = {1.0, 1.0}; + std::array sgn = {1.0, 1.0}; if(func.info->kind==XC_CORRELATION) { if ( rho[0] v1xc(2), v2xc(3), v3xc(2), lapl(2), vlapl(2); + std::array v1xc, v3xc, lapl, vlapl; + std::array v2xc; // call Libxc function: xc_mgga_exc_vxc xc_mgga_exc_vxc( &func, 1, rho.data(), grho.data(), lapl.data(), tau.data(), &s, v1xc.data(), v2xc.data(), vlapl.data(), v3xc.data()); From e5a7df343f55e7e7df8f5ebff92a68253052304d Mon Sep 17 00:00:00 2001 From: linpz Date: Sat, 28 Sep 2024 21:48:50 +0800 Subject: [PATCH 7/7] Fix: remove GlobalV::NSPIN and fix bugs --- source/module_base/global_variable.cpp | 1 - source/module_base/global_variable.h | 2 -- .../module_xc/xc_functional_libxc_tools.cpp | 8 ++++---- .../module_xc/xc_functional_libxc_vxc.cpp | 12 ++++++------ 4 files changed, 10 insertions(+), 13 deletions(-) diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp index f1a39252f6..4d8301f2a1 100644 --- a/source/module_base/global_variable.cpp +++ b/source/module_base/global_variable.cpp @@ -21,7 +21,6 @@ namespace GlobalV int NBANDS = 0; int NLOCAL = 0; // total number of local basis. -int NSPIN = 1; // LDA double nupdown = 0.0; bool use_uspp = false; diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h index 2e0f3e575c..2405e81ab3 100644 --- a/source/module_base/global_variable.h +++ b/source/module_base/global_variable.h @@ -20,8 +20,6 @@ namespace GlobalV extern int NBANDS; extern int NLOCAL; // 1.1 // mohan add 2009-05-29 - -extern int NSPIN; // 7 extern double nupdown; extern bool use_uspp; 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 index 28c1df1d69..c8513001a3 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp @@ -28,7 +28,7 @@ XC_Functional_Libxc::convert_rho_amag_nspin4( const std::size_t nrxx, const Charge* const chr) { - assert(GlobalV::NSPIN==4); + assert(PARAM.inp.nspin==4); std::vector rho(nrxx*nspin); std::vector amag(nrxx); #ifdef _OPENMP @@ -270,9 +270,9 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4( const std::vector &amag, const ModuleBase::matrix &v) { - assert(GlobalV::NSPIN==4); + assert(PARAM.inp.nspin==4); constexpr double vanishing_charge = 1.0e-10; - ModuleBase::matrix v_nspin4(GlobalV::NSPIN, nrxx); + ModuleBase::matrix v_nspin4(PARAM.inp.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]; } } 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 2a2be5c832..a215c02fe9 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 @@ -24,7 +24,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / ModuleBase::timer::tick("XC_Functional_Libxc","v_xc_libxc"); const int nspin = - (GlobalV::NSPIN == 1 || ( GlobalV::NSPIN ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) + (PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) ? 1 : 2; //---------------------------------------------------------- @@ -53,7 +53,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / // converting rho std::vector rho; std::vector amag; - if(1==nspin || 2==GlobalV::NSPIN) + if(1==nspin || 2==PARAM.inp.nspin) { rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); } @@ -118,7 +118,7 @@ std::tuple XC_Functional_Libxc::v_xc_libxc( / v); } // end for( xc_func_type &func : funcs ) - if(4==GlobalV::NSPIN) + if(4==PARAM.inp.nspin) { v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v); } @@ -168,8 +168,8 @@ std::tuple XC_Functional_Li //output of the subroutine double etxc = 0.0; double vtxc = 0.0; - ModuleBase::matrix v(GlobalV::NSPIN,nrxx); - ModuleBase::matrix vofk(GlobalV::NSPIN,nrxx); + ModuleBase::matrix v(PARAM.inp.nspin,nrxx); + ModuleBase::matrix vofk(PARAM.inp.nspin,nrxx); //---------------------------------------------------------- // xc_func_type is defined in Libxc package @@ -178,7 +178,7 @@ std::tuple XC_Functional_Li // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ //---------------------------------------------------------- - const int nspin = GlobalV::NSPIN; + const int nspin = PARAM.inp.nspin; std::vector funcs = XC_Functional_Libxc::init_func( func_id, ( (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED ) ); const std::vector rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr);