diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 5ea326ee34..d99287243d 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -456,7 +456,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_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_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 52fe1cc14a..0b1c30ebd2 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..4db1ef083d 100644 --- a/source/module_hamilt_general/module_xc/CMakeLists.txt +++ b/source/module_hamilt_general/module_xc/CMakeLists.txt @@ -6,12 +6,17 @@ 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_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 ) 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 112be038f9..46ae57e9ef 100644 --- a/source/module_hamilt_general/module_xc/test/CMakeLists.txt +++ b/source/module_hamilt_general/module_xc/test/CMakeLists.txt @@ -4,12 +4,13 @@ AddTest( TARGET XCTest_PBE LIBS parameter MPI::MPI_CXX Libxc::xc # 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( TARGET XCTest_HSE LIBS parameter MPI::MPI_CXX Libxc::xc # 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 ) @@ -17,6 +18,7 @@ AddTest( TARGET XCTest_PZ_SPN LIBS parameter MPI::MPI_CXX Libxc::xc # 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 ) @@ -25,7 +27,10 @@ AddTest( LIBS parameter MPI::MPI_CXX Libxc::xc ${math_libs} psi device container 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_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 @@ -38,7 +43,11 @@ AddTest( TARGET XCTest_SCAN LIBS parameter MPI::MPI_CXX Libxc::xc 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_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 ) @@ -48,13 +57,18 @@ AddTest( LIBS parameter MPI::MPI_CXX Libxc::xc ${math_libs} psi device container 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_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 ../../../module_base/libm/branred.cpp ../../../module_base/libm/sincos.cpp -) +) \ No newline at end of file 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 5bc516b72c..2eed4b3508 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" #define private public #include "module_parameter/parameter.h" @@ -282,7 +283,7 @@ class XCTest_VXC_meta : public XCTest PARAM.input.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); @@ -290,7 +291,7 @@ class XCTest_VXC_meta : public XCTest PARAM.input.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 c3a671d714..713ce512fc 100644 --- a/source/module_hamilt_general/module_xc/xc_functional.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional.cpp @@ -6,6 +6,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(){} @@ -15,11 +19,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; @@ -50,7 +59,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)); @@ -120,7 +129,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); @@ -134,14 +143,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); @@ -199,7 +208,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!"); @@ -240,101 +251,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..ac9c75256a 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,22 +61,15 @@ 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); // 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); -#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 +80,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 +115,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 +125,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 +133,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 +141,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 4ad50c5905..11b69a8580 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 = PARAM.inp.nspin; - if(PARAM.inp.nspin==4) { nspin0 =1; -} - if(PARAM.inp.nspin==4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) { nspin0 = 2; -} + if(PARAM.inp.nspin==4) { nspin0 =1; } + if(PARAM.inp.nspin==4&&(PARAM.globalv.domag||PARAM.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..47ddd1ed1f --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.cpp @@ -0,0 +1,105 @@ +#ifdef USE_LIBXC + +#include "xc_functional_libxc.h" +#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 + +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..73bf03e51c --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc.h @@ -0,0 +1,178 @@ +#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" + extern std::pair> set_xc_type_libxc(const std::string xc_func_in); + + // converts func_id into corresponding xc_func_type vector + extern std::vector init_func(const std::vector &func_id, const int xc_polarized); + + extern void finish_func(std::vector &funcs); + + +//------------------- +// xc_functional_libxc_vxc.cpp +//------------------- + + extern std::tuple v_xc_libxc( + const std::vector &func_id, + const int &nrxx, // number of real-space grid + const double &omega, // volume of cell + const double tpiba, + const Charge* const chr); // charge density + + // for mGGA functional + 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 +//------------------- + + extern 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 + 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 + extern void gcxc_spin_libxc( + const std::vector &func_id, + 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); + + +//------------------- +// xc_functional_libxc_wrapper_tauxc.cpp +//------------------- + + // wrapper for the mGGA functionals + 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); + + 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 // USE_LIBXC + +#endif // XC_FUNCTIONAL_LIBXC_H \ No newline at end of file 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..c8513001a3 --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp @@ -0,0 +1,293 @@ +#ifdef USE_LIBXC + +#include "xc_functional_libxc.h" +#include "xc_functional.h" +#include "module_elecstate/module_charge/charge.h" +#include "module_parameter/parameter.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(PARAM.inp.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(PARAM.inp.nspin==4); + constexpr double vanishing_charge = 1.0e-10; + 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]; + } + } + } + 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 new file mode 100644 index 0000000000..a215c02fe9 --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp @@ -0,0 +1,394 @@ +#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_parameter/parameter.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_Libxc","v_xc_libxc"); + ModuleBase::timer::tick("XC_Functional_Libxc","v_xc_libxc"); + + const int nspin = + (PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) + ? 1 : 2; + + //---------------------------------------------------------- + // 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; + std::vector amag; + if(1==nspin || 2==PARAM.inp.nspin) + { + rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); + } + else + { + std::tuple,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) + { + gdr = XC_Functional_Libxc::cal_gdr(nspin, rho, tpiba, chr); + sigma = XC_Functional_Libxc::convert_sigma(gdr); + } + + double etxc = 0.0; + double vtxc = 0.0; + ModuleBase::matrix v(nspin,nrxx); + + for( xc_func_type &func : funcs ) + { + // jiyy add for threshold + constexpr double rho_threshold = 1E-6; + constexpr double grho_threshold = 1E-10; + + xc_func_set_dens_threshold(&func, rho_threshold); + + // sgn for threshold mask + const std::vector sgn = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma); + + std::vector 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; + } + + 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==PARAM.inp.nspin) + { + v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v); + } + + //------------------------------------------------- + // 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_Libxc","v_xc_libxc"); + return std::make_tuple( etxc, vtxc, std::move(v) ); +} + + +//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_Libxc","v_xc_meta"); + ModuleBase::timer::tick("XC_Functional_Libxc","v_xc_meta"); + + double e2 = 2.0; + + //output of the subroutine + double etxc = 0.0; + double vtxc = 0.0; + ModuleBase::matrix v(PARAM.inp.nspin,nrxx); + ModuleBase::matrix vofk(PARAM.inp.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 = 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); + 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; + 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 ); + + 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 +#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 && XC_Functional::get_func_type() == 5) + { + exc[ir] *= (1.0 - XC_Functional::get_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 && XC_Functional::get_func_type() == 5) + { + 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]; + 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 && XC_Functional::get_func_type() == 5) + { + vsigma[ir] *= (1.0 - XC_Functional::get_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 && XC_Functional::get_func_type() == 5) + { + 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 + + 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; isnumber == XC_MGGA_X_SCAN && XC_Functional::get_func_type() == 5) + { + vtau[ir*nspin+is] *= (1.0 - XC_Functional::get_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_Libxc","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..ce4aad00d8 --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_gcxc.cpp @@ -0,0 +1,79 @@ +#ifdef USE_LIBXC + +#include "xc_functional_libxc.h" + +#include +#include + +void XC_Functional_Libxc::gcxc_libxc( + const std::vector &func_id, + const double &rho, const double &grho, + double &sxc, double &v1xc, double &v2xc) +{ + 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; + 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, + 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) +{ + sxc = v1xcup = v1xcdw = 0.0; + v2xcup = v2xcdw = v2xcud = 0.0; + 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) + { + if( func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) + { + constexpr double rho_threshold = 1E-6; + constexpr double grho_threshold = 1E-10; + std::array sgn = {1.0, 1.0}; + if(func.info->kind==XC_CORRELATION) + { + if ( rho[0] 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]); + v1xcup += v1xc[0] * sgn[0]; + v1xcdw += v1xc[1] * sgn[1]; + v2xcup += 2.0 * v2xc[0] * sgn[0]; + v2xcud += v2xc[1] * sgn[0] * sgn[1]; + v2xcdw += 2.0 * v2xc[2] * sgn[1]; + } + } + XC_Functional_Libxc::finish_func(funcs); +} + +#endif \ No newline at end of file diff --git a/source/module_hamilt_general/module_xc/xc_functional_wrapper_tauxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp similarity index 65% rename from source/module_hamilt_general/module_xc/xc_functional_wrapper_tauxc.cpp rename to source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp index 2cabd64534..b265af9bb1 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_wrapper_tauxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_tauxc.cpp @@ -2,27 +2,31 @@ // it includes 1 subroutine: // 1. tau_xc -#include "xc_functional.h" +#ifdef USE_LIBXC + +#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 -#ifdef USE_LIBXC -void XC_Functional::tau_xc(const double &rho, const double &grho, const double &atau, double &sxc, +void XC_Functional_Libxc::tau_xc( + const std::vector &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; for(xc_func_type &func : funcs) { 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); @@ -35,54 +39,37 @@ 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); - - 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::array rho = {rhoup, rhodw}; + const std::array grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; + const std::array 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::array sgn = {1.0, 1.0}; if(func.info->kind==XC_CORRELATION) { if ( rho[0] 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()); + sxc += s * (rho[0] * sgn[0] + rho[1] * sgn[1]); v1xcup += v1xc[0] * sgn[0]; v1xcdw += v1xc[1] * sgn[1]; @@ -101,16 +95,7 @@ void XC_Functional::tau_xc_spin(double rhoup, double rhodw, } } - delete[] grho; - delete[] rho; - delete[] tau; - delete[] v1xc; - delete[] v2xc; - delete[] v3xc; - delete[] sgn; - finish_func(funcs); - - return; + XC_Functional_Libxc::finish_func(funcs); } #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 new file mode 100644 index 0000000000..aa2cf664ad --- /dev/null +++ b/source/module_hamilt_general/module_xc/xc_functional_libxc_wrapper_xc.cpp @@ -0,0 +1,32 @@ +#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) +{ + const std::vector rho_ud = {rhoup, rhodw}; + exc = vxcup = vxcdw = 0.0; + + 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.data(), &e, vxc_ud.data()); + } + exc += e; + vxcup += vxc_ud[0]; + vxcdw += vxc_ud[1]; + } + + XC_Functional_Libxc::finish_func(funcs); +} + +#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 236a51e3b9..7aeb677a78 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,11 @@ #include "module_base/parallel_reduce.h" #include "module_base/timer.h" #include "module_parameter/parameter.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 +25,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 @@ -125,9 +130,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 { @@ -175,618 +184,3 @@ 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 = - (PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.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==PARAM.inp.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==PARAM.inp.nspin || 2==PARAM.inp.nspin) - { - ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); - return std::make_tuple( etxc, vtxc, std::move(v) ); - } - else if(4==PARAM.inp.nspin) - { - constexpr double vanishing_charge = 1.0e-10; - 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; 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==PARAM.inp.nspin) - else//NSPIN != 1,2,4 is not supported - { - throw std::domain_error("PARAM.inp.nspin ="+std::to_string(PARAM.inp.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(PARAM.inp.nspin,nrxx); - ModuleBase::matrix vofk(PARAM.inp.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 = PARAM.inp.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 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 4b5e7a2aff..f2f36620c4 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp @@ -23,6 +23,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 @@ -56,7 +59,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 2850c58ef7..26d5e66c37 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -5,6 +5,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, @@ -47,7 +52,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_io/input_conv.cpp b/source/module_io/input_conv.cpp index 578b7aa363..732b1e27ec 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -385,7 +385,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; diff --git a/source/module_lr/potentials/kernel_xc.cpp b/source/module_lr/potentials/kernel_xc.cpp index fad91e0aa4..2990ce5d52 100644 --- a/source/module_lr/potentials/kernel_xc.cpp +++ b/source/module_lr/potentials/kernel_xc.cpp @@ -5,13 +5,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) @@ -217,7 +221,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 == PARAM.inp.nspin || 2 == PARAM.inp.nspin) return; // else if (4 == PARAM.inp.nspin)