From 56a18bc0c67c2a4fcc3bb47c0a0268c963c4bcb3 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 30 May 2026 06:32:13 +0800 Subject: [PATCH 1/6] Refactor: Improve code style in xc_functional_libxc files - Add tau_xc overload with hybrid_alpha parameter - Add ElecState::set_exx overload with cal_exx and hybrid_alpha parameters - Format code: one parameter per line, align code, initialize variables, one variable per line, add braces to single-line for/if --- .../module_xc/xc_functional_libxc.h | 404 +++++++++--------- .../xc_functional_libxc_wrapper_gcxc.cpp | 76 ++-- .../xc_functional_libxc_wrapper_tauxc.cpp | 128 ++++-- 3 files changed, 346 insertions(+), 262 deletions(-) diff --git a/source/source_hamilt/module_xc/xc_functional_libxc.h b/source/source_hamilt/module_xc/xc_functional_libxc.h index defd851faa9..6b80d4ffb82 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc.h +++ b/source/source_hamilt/module_xc/xc_functional_libxc.h @@ -1,200 +1,206 @@ -#ifndef XC_FUNCTIONAL_LIBXC_H -#define XC_FUNCTIONAL_LIBXC_H - -#ifdef USE_LIBXC - -#include "source_base/matrix.h" -#include "source_base/vector3.h" - -#include -#include - -#include -#include - -#include // added by jghan, 2024-10-10 -#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); - - /** - * @brief instantiate the XC functional by its ID, and set the external parameters if provided. - * - * @param func_id libxc ID of functional, see https://libxc.gitlab.io/functionals/ for details - * @param xc_polarized 0: unpolarized, 1: spin-polarized - * @return std::vector - * - * @note the functionality of this method is extended by supporting the user-defined - * external parameters of xc. However, there are several functionals' external - * parameters are pre-defined in the code, which herein we call those are - * "in-built" parameters. If the same functional ID is found in both in-built - * and external parameters, the external parameters will overwrite the in-built ones. - * The external parameters can be passed here by keywords xc_exch_ext and - * xc_corr_ext in the input file. The expected format would be an XC ID - * followed by a list of parameters. - */ - 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 - const std::map* scaling_factor = nullptr); // added by jghan, 2024-10-10 - - // 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::size_t nrxx, - 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 std::pair 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); - - // 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 - +#ifndef XC_FUNCTIONAL_LIBXC_H +#define XC_FUNCTIONAL_LIBXC_H + +#ifdef USE_LIBXC + +#include "source_base/matrix.h" +#include "source_base/vector3.h" + +#include +#include + +#include +#include + +#include // added by jghan, 2024-10-10 +#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); + + /** + * @brief instantiate the XC functional by its ID, and set the external parameters if provided. + * + * @param func_id libxc ID of functional, see https://libxc.gitlab.io/functionals/ for details + * @param xc_polarized 0: unpolarized, 1: spin-polarized + * @return std::vector + * + * @note the functionality of this method is extended by supporting the user-defined + * external parameters of xc. However, there are several functionals' external + * parameters are pre-defined in the code, which herein we call those are + * "in-built" parameters. If the same functional ID is found in both in-built + * and external parameters, the external parameters will overwrite the in-built ones. + * The external parameters can be passed here by keywords xc_exch_ext and + * xc_corr_ext in the input file. The expected format would be an XC ID + * followed by a list of parameters. + */ + 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 + const std::map* scaling_factor = nullptr); // added by jghan, 2024-10-10 + + // 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::size_t nrxx, + 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 std::pair 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); + + // 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( + const std::vector &func_id, + const double &rho, const double &grho, const double &atau, double &sxc, + double &v1xc, double &v2xc, double &v3xc, + const double &hybrid_alpha); + + 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/source_hamilt/module_xc/xc_functional_libxc_wrapper_gcxc.cpp b/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_gcxc.cpp index 8869f99ab4d..a758d082597 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_gcxc.cpp +++ b/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_gcxc.cpp @@ -6,11 +6,16 @@ #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) { - sxc = v1xc = v2xc = 0.0; + sxc = 0.0; + v1xc = 0.0; + v2xc = 0.0; constexpr double small = 1.e-6; constexpr double smallg = 1.e-10; @@ -20,56 +25,75 @@ void XC_Functional_Libxc::gcxc_libxc( } std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, + /* func_id = */ func_id, /* xc_polarized = */ XC_UNPOLARIZED); - for(xc_func_type &func : funcs) + + for (xc_func_type& func : funcs) { - double s,v1,v2; + double s = 0.0; + double v1 = 0.0; + double v2 = 0.0; 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) + 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()}; + sxc = 0.0; + v1xcup = 0.0; + v1xcdw = 0.0; + v2xcup = 0.0; + v2xcdw = 0.0; + 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 = */ func_id, + /* func_id = */ func_id, /* xc_polarized = */ XC_POLARIZED); - for(xc_func_type &func : funcs) + for (xc_func_type& func : funcs) { - if( func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) + 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) + std::array sgn = {1.0, 1.0}; + if (func.info->kind == XC_CORRELATION) { - if ( rho[0] v1xc; - std::array v2xc; + std::array v1xc = {0.0, 0.0}; + std::array v2xc = {0.0, 0.0, 0.0}; // call Libxc function: xc_gga_exc_vxc - xc_gga_exc_vxc( &func, 1, rho.data(), grho.data(), &s, v1xc.data(), v2xc.data()); + 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]; @@ -81,4 +105,4 @@ void XC_Functional_Libxc::gcxc_spin_libxc( XC_Functional_Libxc::finish_func(funcs); } -#endif \ No newline at end of file +#endif diff --git a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_tauxc.cpp b/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_tauxc.cpp index a7499ef0906..238fe898a6d 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_tauxc.cpp +++ b/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_tauxc.cpp @@ -10,31 +10,63 @@ #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 +//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 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) + 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; +#ifdef __EXX + tau_xc(func_id, rho, grho, atau, sxc, v1xc, v2xc, v3xc, GlobalC::exx_info.info_global.hybrid_alpha); +#else + tau_xc(func_id, rho, grho, atau, sxc, v1xc, v2xc, v3xc, 0.0); +#endif +} + +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, + const double& hybrid_alpha) +{ + double s = 0.0; + double v1 = 0.0; + double v2 = 0.0; + double v3 = 0.0; + double lapl_rho = 0.0; + double vlapl_rho = 0.0; lapl_rho = grho; std::vector funcs = XC_Functional_Libxc::init_func( - /* func_id = */ func_id, + /* func_id = */ func_id, /* xc_polarized = */ XC_UNPOLARIZED); - sxc = 0.0; v1xc = 0.0; v2xc = 0.0; v3xc = 0.0; + sxc = 0.0; + v1xc = 0.0; + v2xc = 0.0; + v3xc = 0.0; - for(xc_func_type &func : funcs) + for (xc_func_type& func : funcs) { - xc_mgga_exc_vxc(&func,1,&rho,&grho,&lapl_rho,&atau,&s,&v1,&v2,&vlapl_rho,&v3); + 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 && 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); - v2 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha); - v3 *= (1.0 - GlobalC::exx_info.info_global.hybrid_alpha); + s *= (1.0 - hybrid_alpha); + v1 *= (1.0 - hybrid_alpha); + v2 *= (1.0 - hybrid_alpha); + v3 *= (1.0 - hybrid_alpha); } #endif sxc += s * rho; @@ -44,50 +76,72 @@ void XC_Functional_Libxc::tau_xc( } XC_Functional_Libxc::finish_func(funcs); - return; + return; } 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) + 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) { - sxc = v1xcup = v1xcdw = 0.0; - v2xcup = v2xcdw = v2xcud = 0.0; - v3xcup = v3xcdw = 0.0; + sxc = 0.0; + v1xcup = 0.0; + v1xcdw = 0.0; + v2xcup = 0.0; + v2xcdw = 0.0; + v2xcud = 0.0; + v3xcup = 0.0; + v3xcdw = 0.0; - const std::array rho = {rhoup, rhodw}; - const std::array grho = {gdr1.norm2(), gdr1 * gdr2, gdr2.norm2()}; - const std::array 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 = */ func_id, + /* func_id = */ func_id, /* xc_polarized = */ XC_POLARIZED); - for(xc_func_type &func : funcs) + for (xc_func_type& func : funcs) { - if( func.info->family == XC_FAMILY_MGGA || func.info->family == XC_FAMILY_HYB_MGGA) + if (func.info->family == XC_FAMILY_MGGA || func.info->family == XC_FAMILY_HYB_MGGA) { 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) + std::array sgn = {1.0, 1.0}; + if (func.info->kind == XC_CORRELATION) { - if ( rho[0] v1xc, v3xc, lapl, vlapl; - std::array v2xc; + std::array v1xc = {0.0, 0.0}; + std::array v3xc = {0.0, 0.0}; + std::array lapl = {0.0, 0.0}; + std::array vlapl = {0.0, 0.0}; + std::array v2xc = {0.0, 0.0, 0.0}; // 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()); + 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]; From bbab15d026ecba7f1988f27da515459f993fdfa1 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 30 May 2026 07:04:18 +0800 Subject: [PATCH 2/6] update --- source/source_hamilt/module_xc/exx_info.h | 24 +- .../module_xc/xc_funct_corr_gga.cpp | 299 +++++++++++---- .../module_xc/xc_funct_corr_lda.cpp | 270 ++++++++------ .../module_xc/xc_funct_exch_gga.cpp | 260 ++++++++----- .../module_xc/xc_funct_exch_lda.cpp | 88 +++-- .../source_hamilt/module_xc/xc_funct_hcth.cpp | 74 +++- .../source_hamilt/module_xc/xc_functional.h | 352 ++++++++++++------ .../module_xc/xc_functional_libxc.h | 79 ++-- .../module_xc/xc_functional_wrapper_xc.cpp | 176 ++++++--- 9 files changed, 1098 insertions(+), 524 deletions(-) diff --git a/source/source_hamilt/module_xc/exx_info.h b/source/source_hamilt/module_xc/exx_info.h index 66e4a9d4f67..758ed8efb5b 100644 --- a/source/source_hamilt/module_xc/exx_info.h +++ b/source/source_hamilt/module_xc/exx_info.h @@ -13,17 +13,17 @@ struct Exx_Info { bool cal_exx = false; - std::map>> coulomb_param; + std::map>> coulomb_param; - // Fock: - // "alpha": "0" - // "singularity_correction": "limits" / "spencer" / "revised_spencer" / "massidda" / "carrier" - // "lambda": "0.3" + // Fock: + // "alpha": "0" + // "singularity_correction": "limits" / "spencer" / "revised_spencer" / "massidda" / "carrier" + // "lambda": "0.3" // "Rcut" - // Erfc: - // "alpha": "0" - // "omega": "0.11" - // "singularity_correction": "limits" / "spencer" / "revised_spencer" + // Erfc: + // "alpha": "0" + // "omega": "0.11" + // "singularity_correction": "limits" / "spencer" / "revised_spencer" // "Rcut" Conv_Coulomb_Pot_K::Ccp_Type ccp_type; @@ -43,14 +43,14 @@ struct Exx_Info double lambda = 0.3; Exx_Info_Lip(const Exx_Info::Exx_Info_Global& info_global) - :ccp_type(info_global.ccp_type), - hse_omega(info_global.hse_omega) {} + : ccp_type(info_global.ccp_type), + hse_omega(info_global.hse_omega) {} }; Exx_Info_Lip info_lip; struct Exx_Info_RI { - const std::map>> &coulomb_param; + const std::map>> &coulomb_param; bool real_number = false; bool coul_moment = false; diff --git a/source/source_hamilt/module_xc/xc_funct_corr_gga.cpp b/source/source_hamilt/module_xc/xc_funct_corr_gga.cpp index 3c88eb653c8..68e006cd75c 100644 --- a/source/source_hamilt/module_xc/xc_funct_corr_gga.cpp +++ b/source/source_hamilt/module_xc/xc_funct_corr_gga.cpp @@ -11,14 +11,25 @@ #include "xc_functional.h" -void XC_Functional::perdew86(const double rho, const double grho, double &sc, double &v1c, double &v2c) +void XC_Functional::perdew86( + const double rho, + const double grho, + double &sc, + double &v1c, + double &v2c) { // Perdew gradient correction on correlation: PRB 33, 8822 (1986) // USE kinds // implicit none // real(kind=DP) :: rho, grho, sc, v1c, v2c - double p1, p2, p3, p4, pc1, pc2, pci; + double p1 = 0.0; + double p2 = 0.0; + double p3 = 0.0; + double p4 = 0.0; + double pc1 = 0.0; + double pc2 = 0.0; + double pci = 0.0; // parameter : p1 = 0.0232660; p2 = 7.389e-6; @@ -28,12 +39,25 @@ void XC_Functional::perdew86(const double rho, const double grho, double &sc, do pc1 = 0.0016670; pc2 = 0.0025680; pci = pc1 + pc2; - double third, pi34; + double third = 0.0; + double pi34 = 0.0; // parameter : third = 1.0 / 3.0; pi34 = 0.62035049089940; // pi34=(3/4pi)^(1/3) - double rho13, rho43, rs, rs2, rs3, cna, cnb, cn, drs; - double dcna, dcnb, dcn, phi, ephi; + double rho13 = 0.0; + double rho43 = 0.0; + double rs = 0.0; + double rs2 = 0.0; + double rs3 = 0.0; + double cna = 0.0; + double cnb = 0.0; + double cn = 0.0; + double drs = 0.0; + double dcna = 0.0; + double dcnb = 0.0; + double dcn = 0.0; + double phi = 0.0; + double ephi = 0.0; rho13 = pow(rho, third); rho43 = pow(rho13, 4); @@ -58,10 +82,22 @@ void XC_Functional::perdew86(const double rho, const double grho, double &sc, do return; } //end subroutine perdew86 -void XC_Functional::ggac(const double &rho,const double &grho, double &sc, double &v1c, double &v2c) +void XC_Functional::ggac( + const double &rho, + const double &grho, + double &sc, + double &v1c, + double &v2c) { // Perdew-Wang GGA (PW91) correlation part - double al, pa, pb, pc, pd, cx, cxc0, cc0; + double al = 0.0; + double pa = 0.0; + double pb = 0.0; + double pc = 0.0; + double pd = 0.0; + double cx = 0.0; + double cxc0 = 0.0; + double cc0 = 0.0; // parameter : al = 0.090; pa = 0.0232660; @@ -72,7 +108,12 @@ void XC_Functional::ggac(const double &rho,const double &grho, double &sc, doubl cx = -0.0016670; cxc0 = 0.0025680; cc0 = - cx + cxc0; - double third, pi34, nu, be, xkf, xks; + double third = 0.0; + double pi34 = 0.0; + double nu = 0.0; + double be = 0.0; + double xkf = 0.0; + double xks = 0.0; // parameter : third = 1.0 / 3.0; pi34 = 0.62035049089940; @@ -84,17 +125,41 @@ void XC_Functional::ggac(const double &rho,const double &grho, double &sc, doubl xks = 1.1283791670955130; // pi34=(3/4pi)^(1/3), nu=(16/pi)*(3 pi^2)^(1/3) // xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi) - double kf, ks, rs, rs2, rs3, ec, vc, t, expe, af, bf, y, xy, - qy, s1; - double h0, dh0, ddh0, ee, cn, dcn, cna, dcna, cnb, dcnb, h1, - dh1, ddh1; + double kf = 0.0; + double ks = 0.0; + double rs = 0.0; + double rs2 = 0.0; + double rs3 = 0.0; + double ec = 0.0; + double vc = 0.0; + double t = 0.0; + double expe = 0.0; + double af = 0.0; + double bf = 0.0; + double y = 0.0; + double xy = 0.0; + double qy = 0.0; + double s1 = 0.0; + double h0 = 0.0; + double dh0 = 0.0; + double ddh0 = 0.0; + double ee = 0.0; + double cn = 0.0; + double dcn = 0.0; + double cna = 0.0; + double dcna = 0.0; + double cnb = 0.0; + double dcnb = 0.0; + double h1 = 0.0; + double dh1 = 0.0; + double ddh1 = 0.0; - rs = pi34 / pow(rho, third); + rs = pi34 / pow(rho, third); rs2 = rs * rs; rs3 = rs * rs2; - // call function: pw - XC_Functional::pw(rs, 0, ec, vc); + // call function: pw + XC_Functional::pw(rs, 0, ec, vc); kf = xkf / rs; ks = xks * sqrt(kf); t = sqrt(grho) / (2.0 * ks * rho); @@ -132,60 +197,82 @@ void XC_Functional::ggac(const double &rho,const double &grho, double &sc, doubl return; } //end subroutine ggac -void XC_Functional::pbec(const double &rho, const double &grho, const int &iflag, double &sc, double &v1c, double &v2c) +void XC_Functional::pbec( + const double &rho, + const double &grho, + const int &iflag, + double &sc, + double &v1c, + double &v2c) { - // PBE correlation (without LDA part) - // iflag=0: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996). - // iflag=1: J.P.Perdew et al., PRL 100, 136406 (2008). - const double ga = 0.0310906908696548950; - const double be[2] = {0.06672455060314922, 0.046}; - - const double third = 1.0 / 3.0; - const double pi34 = 0.62035049089940; - const double xkf = 1.9191582926775130; - const double xks = 1.1283791670955130; - - // pi34=(3/4pi)^(1/3), xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi) - double ec, vc; - - const double rs = pi34 / pow(rho, third); - - XC_Functional::pw(rs, 0, ec, vc); - - const double kf = xkf / rs; - const double ks = xks * sqrt(kf); - const double t = sqrt(grho) / (2.0 * ks * rho); - const double expe = exp(- ec / ga); - const double af = be[iflag] / ga * (1.0 / (expe - 1.0)); - const double bf = expe * (vc - ec); - const double y = af * t * t; - const double xy = (1.0 + y) / (1.0 + y + y * y); - - const double x = 1.0 + y + y * y; - const double qy = y * y * (2.0 + y) / (x * x); - const double s1 = 1.0 + be[iflag] / ga * t * t * xy; - const double h0 = ga * log(s1); - const double dh0 = be[iflag] * t * t / s1 * (- 7.0 / 3.0 * xy - qy * (af * bf /be[iflag] - 7.0 / 3.0)); - const double ddh0 = be[iflag] / (2.0 * ks * ks * rho) * (xy - qy) / s1; - - sc = rho * h0; - v1c = h0 + dh0; + // PBE correlation (without LDA part) + // iflag=0: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996). + // iflag=1: J.P.Perdew et al., PRL 100, 136406 (2008). + const double ga = 0.0310906908696548950; + const double be[2] = {0.06672455060314922, 0.046}; + + const double third = 1.0 / 3.0; + const double pi34 = 0.62035049089940; + const double xkf = 1.9191582926775130; + const double xks = 1.1283791670955130; + + // pi34=(3/4pi)^(1/3), xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi) + double ec = 0.0; + double vc = 0.0; + + const double rs = pi34 / pow(rho, third); + + XC_Functional::pw(rs, 0, ec, vc); + + const double kf = xkf / rs; + const double ks = xks * sqrt(kf); + const double t = sqrt(grho) / (2.0 * ks * rho); + const double expe = exp(- ec / ga); + const double af = be[iflag] / ga * (1.0 / (expe - 1.0)); + const double bf = expe * (vc - ec); + const double y = af * t * t; + const double xy = (1.0 + y) / (1.0 + y + y * y); + + const double x = 1.0 + y + y * y; + const double qy = y * y * (2.0 + y) / (x * x); + const double s1 = 1.0 + be[iflag] / ga * t * t * xy; + const double h0 = ga * log(s1); + const double dh0 = be[iflag] * t * t / s1 * (- 7.0 / 3.0 * xy - qy * (af * bf /be[iflag] - 7.0 / 3.0)); + const double ddh0 = be[iflag] / (2.0 * ks * ks * rho) * (xy - qy) / s1; + + sc = rho * h0; + v1c = h0 + dh0; v2c = ddh0; - return; + return; } -void XC_Functional::glyp(const double &rho, const double &grho, double &sc, double &v1c, double &v2c) +void XC_Functional::glyp( + const double &rho, + const double &grho, + double &sc, + double &v1c, + double &v2c) { //----------------------------------------------------------------------- // Lee Yang Parr: gradient correction part - double a, b, c, d; + double a = 0.0; + double b = 0.0; + double c = 0.0; + double d = 0.0; a = 0.049180; b = 0.1320; c = 0.25330; d = 0.3490; - double rhom13, rhom43, rhom53, om, xl, ff, dom, dxl; + double rhom13 = 0.0; + double rhom43 = 0.0; + double rhom53 = 0.0; + double om = 0.0; + double xl = 0.0; + double ff = 0.0; + double dom = 0.0; + double dxl = 0.0; rhom13 = pow(rho, (- 1.0 / 3.0)); om = exp(- c * rhom13) / (1.0 + d * rhom13); @@ -196,7 +283,7 @@ void XC_Functional::glyp(const double &rho, const double &grho, double &sc, doub sc = ff * rhom53 * om * xl; dom = - om * (c + d + c * d * rhom13) / (1.0 + d * rhom13); - double x; + double x = 0.0; x = 1.0 + d * rhom13; dxl = (7.0 / 3.0) * (c + d + 2.0 * c * d * rhom13 + c * d * d * rhom13 * rhom13) / (x * x); @@ -209,8 +296,14 @@ void XC_Functional::glyp(const double &rho, const double &grho, double &sc, doub } // end subroutine glyp //----------------------------------------------------------------------- -void XC_Functional::perdew86_spin(double rho, double zeta, double grho, double &sc, - double &v1cup, double &v1cdw, double &v2c) +void XC_Functional::perdew86_spin( + double rho, + double zeta, + double grho, + double &sc, + double &v1cup, + double &v1cdw, + double &v2c) { //------------------------------------------------------------------- // Perdew gradient correction on correlation: PRB 33, 8822 (1986) @@ -219,7 +312,13 @@ void XC_Functional::perdew86_spin(double rho, double zeta, double grho, double & // USE kinds // implicit none // real(kind=DP) :: rho, zeta, grho, sc, v1cup, v1cdw, v2c - double p1, p2, p3, p4, pc1, pc2, pci; + double p1 = 0.0; + double p2 = 0.0; + double p3 = 0.0; + double p4 = 0.0; + double pc1 = 0.0; + double pc2 = 0.0; + double pci = 0.0; // parameter : p1 = 0.0232660; p2 = 7.389e-6; @@ -229,14 +328,29 @@ void XC_Functional::perdew86_spin(double rho, double zeta, double grho, double & pc1 = 0.0016670; pc2 = 0.0025680; pci = pc1 + pc2; - double third, pi34; + double third = 0.0; + double pi34 = 0.0; // parameter : third = 1.0 / 3.0; pi34 = 0.62035049089940; // pi34=(3/4pi)^(1/3) - double rho13, rho43, rs, rs2, rs3, cna, cnb, cn, drs; - double dcna, dcnb, dcn, phi, ephi, dd, ddd; + double rho13 = 0.0; + double rho43 = 0.0; + double rs = 0.0; + double rs2 = 0.0; + double rs3 = 0.0; + double cna = 0.0; + double cnb = 0.0; + double cn = 0.0; + double drs = 0.0; + double dcna = 0.0; + double dcnb = 0.0; + double dcn = 0.0; + double phi = 0.0; + double ephi = 0.0; + double dd = 0.0; + double ddd = 0.0; rho13 = pow(rho, third); rho43 = pow(rho13, 4); @@ -314,7 +428,7 @@ void XC_Functional::ggac_spin(double rho, double zeta, double grho, double &sc, rs = pi34 / pow(rho, third); rs2 = rs * rs; rs3 = rs * rs2; - XC_Functional::pw_spin(rs, zeta, ec, vcup, vcdw); + XC_Functional::pw_spin(rs, zeta, ec, vcup, vcdw); kf = xkf / rs; ks = xks * sqrt(kf); fz = 0.50 * (pow((1.0 + zeta) , (2.0 / 3.0)) + pow((1.0 - zeta) , ( @@ -331,7 +445,7 @@ void XC_Functional::ggac_spin(double rho, double zeta, double grho, double &sc, //bfdw = expe * (vcdw - ec) / fz3; y = af * t * t; xy = (1.0 + y) / (1.0 + y + y * y); - qy = y * y * (2.0 + y) / (1.0 + y + y * y) ; // **2; + qy = y * y * (2.0 + y) / (1.0 + y + y * y) ; // **2; qy *= qy; s1 = 1.0 + 2.0 * al / be * t * t * xy; h0 = fz3 * be * be / (2.0 * al) * log(s1); @@ -346,7 +460,7 @@ void XC_Functional::ggac_spin(double rho, double zeta, double grho, double &sc, qy * (3.0 * af * expe * ec / fz3 / be + 2.0))) * dfz * (1.0 + zeta); ddh0 = be * fz / (2.0 * ks * ks * rho) * (xy - qy) / s1; - ee = - 100.0 * fz4 * (ks / kf * t); // **2; + ee = - 100.0 * fz4 * (ks / kf * t); // **2; ee *= ee; cna = cxc0 + pa * rs + pb * rs2; dcna = pa * rs + 2.0 * pb * rs2; @@ -369,8 +483,15 @@ void XC_Functional::ggac_spin(double rho, double zeta, double grho, double &sc, */ //--------------------------------------------------------------- -void XC_Functional::pbec_spin(double rho, double zeta, double grho, const int &iflag, double &sc, - double &v1cup, double &v1cdw, double &v2c) +void XC_Functional::pbec_spin( + double rho, + double zeta, + double grho, + const int &iflag, + double &sc, + double &v1cup, + double &v1cdw, + double &v2c) { //----------------------------------------------------------- @@ -380,12 +501,16 @@ void XC_Functional::pbec_spin(double rho, double zeta, double grho, const int &i // USE kinds // implicit none // real(kind=DP) :: rho, zeta, grho, sc, v1cup, v1cdw, v2c - double ga, be[3];//mohan add + double ga = 0.0; + double be[3] = {0.0, 0.0, 0.0};//mohan add // parameter : ga = 0.0310910; be[1] = 0.06672455060314922;//zhengdy add 2019-09-12, ensure same parameter with another dft code. - be[2] = 0.0460000;//mohan add 2012-05-28 - double third, pi34, xkf, xks; + be[2] = 0.0460000;//mohan add 2012-05-28 + double third = 0.0; + double pi34 = 0.0; + double xkf = 0.0; + double xks = 0.0; // parameter : third = 1.0 / 3.0; pi34 = 0.62035049089940; @@ -393,12 +518,34 @@ void XC_Functional::pbec_spin(double rho, double zeta, double grho, const int &i xkf = 1.9191582926775130; xks = 1.1283791670955130; // pi34=(3/4pi)^(1/3), xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi) - double kf, ks, rs, ec, vcup, vcdw, t, expe, af, y, xy, qy, - s1, h0, ddh0; - double fz, fz2, fz3, dfz, bfup, bfdw, dh0up, dh0dw, dh0zup, dh0zdw; // fz4 + double kf = 0.0; + double ks = 0.0; + double rs = 0.0; + double ec = 0.0; + double vcup = 0.0; + double vcdw = 0.0; + double t = 0.0; + double expe = 0.0; + double af = 0.0; + double y = 0.0; + double xy = 0.0; + double qy = 0.0; + double s1 = 0.0; + double h0 = 0.0; + double ddh0 = 0.0; + double fz = 0.0; + double fz2 = 0.0; + double fz3 = 0.0; + double dfz = 0.0; + double bfup = 0.0; + double bfdw = 0.0; + double dh0up = 0.0; + double dh0dw = 0.0; + double dh0zup = 0.0; + double dh0zdw = 0.0; // fz4 rs = pi34 / pow(rho, third); - XC_Functional::pw_spin(rs, zeta, ec, vcup, vcdw); //mohan fix bug 2012-05-28 + XC_Functional::pw_spin(rs, zeta, ec, vcup, vcdw); //mohan fix bug 2012-05-28 kf = xkf / rs; ks = xks * sqrt(kf); fz = 0.50 * (pow((1.0 + zeta) , (2.0 / 3.0)) + pow((1.0 - zeta) , ( @@ -435,4 +582,4 @@ void XC_Functional::pbec_spin(double rho, double zeta, double grho, const int &i v1cdw = h0 + dh0dw + dh0zdw; v2c = ddh0; return; -} // end subroutine pbec_spin \ No newline at end of file +} // end subroutine pbec_spin diff --git a/source/source_hamilt/module_xc/xc_funct_corr_lda.cpp b/source/source_hamilt/module_xc/xc_funct_corr_lda.cpp index 5b6ff5367b2..abd89ed4f22 100644 --- a/source/source_hamilt/module_xc/xc_funct_corr_lda.cpp +++ b/source/source_hamilt/module_xc/xc_funct_corr_lda.cpp @@ -13,7 +13,11 @@ #include "xc_functional.h" -void XC_Functional::pw(const double &rs, const int &iflag, double &ec, double &vc) +void XC_Functional::pw( + const double &rs, + const int &iflag, + double &ec, + double &vc) { const double a = 0.0310910; const double b1 = 7.59570; @@ -24,15 +28,21 @@ void XC_Functional::pw(const double &rs, const int &iflag, double &ec, double &v const double c3 = 0.010430; const double d0 = 0.43350; const double d1 = 1.44080; - double lnrs, rs12, rs32, rs2, om, dom, olog; + double lnrs = 0.0; + double rs12 = 0.0; + double rs32 = 0.0; + double rs2 = 0.0; + double om = 0.0; + double dom = 0.0; + double olog = 0.0; double a1[2] = { 0.213700, 0.0264810 }; - double b3[2] = { 1.63820, -0.466470 }; - double b4[2] = { 0.492940, 0.133540 }; + double b3[2] = { 1.63820, -0.466470 }; + double b4[2] = { 0.492940, 0.133540 }; - // high- and low-density formulae implemented but not used in PW case - // (reason: inconsistencies in PBE/PW91 functionals) - - if (rs < 1 && iflag == 1) + // high- and low-density formulae implemented but not used in PW case + // (reason: inconsistencies in PBE/PW91 functionals) + + if (rs < 1 && iflag == 1) { // high density formula lnrs = log(rs); @@ -52,26 +62,28 @@ void XC_Functional::pw(const double &rs, const int &iflag, double &ec, double &v rs12 = sqrt(rs); rs32 = rs * rs12; rs2 = rs * rs; - om = 2.0 * a * (b1 * rs12 + b2 * rs + b3 [iflag] * rs32 + b4 [ - iflag] * rs2); - dom = 2.0 * a * (0.50 * b1 * rs12 + b2 * rs + 1.50 * b3 [ - iflag] * rs32 + 2.0 * b4 [iflag] * rs2); + om = 2.0 * a * (b1 * rs12 + b2 * rs + b3[iflag] * rs32 + b4[iflag] * rs2); + dom = 2.0 * a * (0.50 * b1 * rs12 + b2 * rs + 1.50 * b3[iflag] * rs32 + 2.0 * b4[iflag] * rs2); olog = log(1.0 + 1.0 / om); - ec = - 2.0 * a * (1.0 + a1 [iflag] * rs) * olog; - vc = - 2.0 * a * (1.0 + 2.0 / 3.0 * a1 [iflag] * rs) - * olog - 2.0 / 3.0 * a * (1.0 + a1 [iflag] * rs) * dom / + ec = - 2.0 * a * (1.0 + a1[iflag] * rs) * olog; + vc = - 2.0 * a * (1.0 + 2.0 / 3.0 * a1[iflag] * rs) + * olog - 2.0 / 3.0 * a * (1.0 + a1[iflag] * rs) * dom / (om * (om + 1.0)); } - return; + return; } //LDA parameterization form Monte Carlo data //iflag=0: J.P. Perdew and A. Zunger, PRB 23, 5048 (1981) //iflag=1: G. Ortiz and P. Ballone, PRB 50, 1391 (1994) -void XC_Functional::pz(const double &rs, const int &iflag, double &ec, double &vc) +void XC_Functional::pz( + const double &rs, + const int &iflag, + double &ec, + double &vc) { // ModuleBase::TITLE("XC_Functional","pz"); - const double a[2] = {0.0311, 0.031091}; + const double a[2] = {0.0311, 0.031091}; const double b[2] = { -0.048, -0.046644 }; const double c[2] = { 0.0020, 0.00419 }; const double d[2] = { -0.0116, -0.00983 }; @@ -83,19 +95,19 @@ void XC_Functional::pz(const double &rs, const int &iflag, double &ec, double &v { // high density formula const double lnrs = log(rs); - ec = a [iflag] * lnrs + b [iflag] + c [iflag] * rs * lnrs + d [iflag] * rs; - vc = a [iflag] * lnrs + (b [iflag] - a [iflag] / 3.0) + 2.0 / - 3.0 * c [iflag] * rs * lnrs + (2.0 * d [iflag] - c [iflag]) + ec = a[iflag] * lnrs + b[iflag] + c[iflag] * rs * lnrs + d[iflag] * rs; + vc = a[iflag] * lnrs + (b[iflag] - a[iflag] / 3.0) + 2.0 / + 3.0 * c[iflag] * rs * lnrs + (2.0 * d[iflag] - c[iflag]) / 3.0 * rs; } else { // interpolation formula const double rs12 = sqrt(rs); - const double ox = 1.0 + b1 [iflag] * rs12 + b2 [iflag] * rs; - const double dox = 1.0 + 7.0 / 6.0 * b1 [iflag] * rs12 + 4.0 / 3.0 * - b2 [iflag] * rs; - ec = gc [iflag] / ox; + const double ox = 1.0 + b1[iflag] * rs12 + b2[iflag] * rs; + const double dox = 1.0 + 7.0 / 6.0 * b1[iflag] * rs12 + 4.0 / 3.0 * + b2[iflag] * rs; + ec = gc[iflag] / ox; vc = ec * dox / ox; } return; @@ -103,7 +115,10 @@ void XC_Functional::pz(const double &rs, const int &iflag, double &ec, double &v // C. Lee, W. Yang, and R.G. Parr, PRB 37, 785 (1988) // LDA part only -void XC_Functional::lyp(const double &rs, double &ec, double &vc) +void XC_Functional::lyp( + const double &rs, + double &ec, + double &vc) { const double a = 0.04918e0; const double b = 0.1320 * 2.87123400018819108e0; @@ -111,7 +126,8 @@ void XC_Functional::lyp(const double &rs, double &ec, double &vc) const double pi43 = 1.61199195401647e0; const double c = 0.2533e0 * pi43; const double d = 0.349e0 * pi43; - double ecrs, ox; + double ecrs = 0.0; + double ox = 0.0; ecrs = b * exp(- c * rs); ox = 1.0 / (1.0 + d * rs); @@ -122,13 +138,24 @@ void XC_Functional::lyp(const double &rs, double &ec, double &vc) } // S.H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980) -void XC_Functional::vwn(const double &rs, double &ec, double &vc) +void XC_Functional::vwn( + const double &rs, + double &ec, + double &vc) { const double a = 0.0310907; const double b = 3.72744; const double c = 12.9352; const double x0 = -0.10498; - double q, f1, f2, f3, rs12, fx, qx, tx, tt; + double q = 0.0; + double f1 = 0.0; + double f2 = 0.0; + double f3 = 0.0; + double rs12 = 0.0; + double fx = 0.0; + double qx = 0.0; + double tx = 0.0; + double tt = 0.0; q = sqrt(4.0 * c - b * b); f1 = 2.0 * b / q; @@ -146,17 +173,20 @@ void XC_Functional::vwn(const double &rs, double &ec, double &vc) tt = tx * tx + q * q; vc = ec - rs12 * a / 6.0 * (2.0 / rs12 - tx / fx - 4.0 * b / - tt - f2 * (2.0 / (rs12 - x0) - tx / fx - 4.0 * (2.0 * x0 + b)/ tt)); + tt - f2 * (2.0 / (rs12 - x0) - tx / fx - 4.0 * (2.0 * x0 + b) / tt)); return; } -void XC_Functional::wigner( const double &rs, double &ec, double &vc) +void XC_Functional::wigner( + const double &rs, + double &ec, + double &vc) { const double pi34 = 0.6203504908994e0; // pi34=(3/4pi)^(1/3), rho13=rho^(1/3) - double rho13 = pi34 / rs; + double rho13 = pi34 / rs; double x = 1.0 + 12.570 * rho13; vc = - rho13 * ((0.9436560 + 8.89630 * rho13) / (x * x)); ec = - 0.7380 * rho13 * (0.9590 / (1.0 + 12.570 * rho13)); @@ -164,9 +194,13 @@ void XC_Functional::wigner( const double &rs, double &ec, double &vc) } // L. Hedin and B.I. Lundqvist, J. Phys. C 4, 2064 (1971) -void XC_Functional::hl( const double &rs, double &ec, double &vc) +void XC_Functional::hl( + const double &rs, + double &ec, + double &vc) { - double a, x; + double a = 0.0; + double x = 0.0; a = log(1.00 + 21.0 / rs); x = rs / 21.00; ec = a + (pow(x, 3) * a - x * x) + x / 2.0 - 1.00 / 3.00; @@ -176,7 +210,10 @@ void XC_Functional::hl( const double &rs, double &ec, double &vc) } // O. Gunnarsson and B. I. Lundqvist, PRB 13, 4274 (1976) -void XC_Functional::gl( const double &rs, double &ec, double &vc) +void XC_Functional::gl( + const double &rs, + double &ec, + double &vc) { const double c = 0.0333; const double r = 11.4; @@ -185,14 +222,18 @@ void XC_Functional::gl( const double &rs, double &ec, double &vc) vc = - c * log(1.0 + 1.0 / x); ec = - c * ((1.0 + pow(x, 3)) * log(1.0 + 1.0 / x) - 1.00 / 3.00 + x * (0.50 - x)); - return; + return; } // J.P. Perdew and Y. Wang, PRB 45, 13244 (1992) -void XC_Functional::pw_spin( const double &rs, const double &zeta, - double &ec, double &vcup, double &vcdw) +void XC_Functional::pw_spin( + const double &rs, + const double &zeta, + double &ec, + double &vcup, + double &vcdw) { - const double a = 0.0310910; + const double a = 0.0310910; const double a1 = 0.213700; const double b1 = 7.59570; const double b2 = 3.58760; @@ -204,7 +245,7 @@ void XC_Functional::pw_spin( const double &rs, const double &zeta, //const double c3 = 0.010430; //const double d0 = 0.43350; //const double d1 = 1.44080; - // polarized parameters + // polarized parameters const double ap = 0.0155450; const double a1p = 0.205480; const double b1p = 14.11890; @@ -217,7 +258,7 @@ void XC_Functional::pw_spin( const double &rs, const double &zeta, //const double c3p = 0.003840; //const double d0p = 0.32870; //const double d1p = 1.76970; - // antiferro + // antiferro const double aa = 0.0168870; const double a1a = 0.111250; const double b1a = 10.3570; @@ -231,13 +272,32 @@ void XC_Functional::pw_spin( const double &rs, const double &zeta, //const double d0a = 0.22400; //const double d1a = 0.39690; - const double fz0 = 1.7099210; - double rs12, rs32, rs2, zeta2, zeta3, zeta4, fz, dfz; - double om, dom, olog, epwc, vpwc; - double omp, domp, ologp, epwcp, vpwcp; - double oma, doma, ologa, alpha, vpwca; + const double fz0 = 1.7099210; + double rs12 = 0.0; + double rs32 = 0.0; + double rs2 = 0.0; + double zeta2 = 0.0; + double zeta3 = 0.0; + double zeta4 = 0.0; + double fz = 0.0; + double dfz = 0.0; + double om = 0.0; + double dom = 0.0; + double olog = 0.0; + double epwc = 0.0; + double vpwc = 0.0; + double omp = 0.0; + double domp = 0.0; + double ologp = 0.0; + double epwcp = 0.0; + double vpwcp = 0.0; + double oma = 0.0; + double doma = 0.0; + double ologa = 0.0; + double alpha = 0.0; + double vpwca = 0.0; - // if(rs.lt.0.5d0) then + // if(rs.lt.0.5d0) then // high density formula (not implemented) // // else if(rs.gt.100.d0) then @@ -253,55 +313,50 @@ void XC_Functional::pw_spin( const double &rs, const double &zeta, rs2 = rs * rs; // unpolarised om = 2.0 * a * (b1 * rs12 + b2 * rs + b3 * rs32 + b4 * rs2); - dom = 2.0 * a * (0.50 * b1 * rs12 + b2 * rs + 1.50 * b3 * rs32 - + 2.0 * b4 * rs2); + dom = 2.0 * a * (0.50 * b1 * rs12 + b2 * rs + 1.50 * b3 * rs32 + 2.0 * b4 * rs2); olog = log(1.0 + 1.00 / om); epwc = - 2.0 * a * (1.0 + a1 * rs) * olog; vpwc = - 2.0 * a * (1.0 + 2.0 / 3.0 * a1 * rs) * olog - 2.0 / 3.0 * a * (1.0 + a1 * rs) * dom / (om * (om + 1.0)); // polarized omp = 2.0 * ap * (b1p * rs12 + b2p * rs + b3p * rs32 + b4p * rs2); - domp = 2.0 * ap * (0.50 * b1p * rs12 + b2p * rs + 1.50 * b3p * - rs32 + 2.0 * b4p * rs2); + domp = 2.0 * ap * (0.50 * b1p * rs12 + b2p * rs + 1.50 * b3p * rs32 + 2.0 * b4p * rs2); ologp = log(1.0 + 1.00 / omp); epwcp = - 2.0 * ap * (1.0 + a1p * rs) * ologp; vpwcp = - 2.0 * ap * (1.0 + 2.0 / 3.0 * a1p * rs) * ologp - 2.0 / 3.0 * ap * (1.0 + a1p * rs) * domp / (omp * (omp + 1.0)); // antiferro oma = 2.0 * aa * (b1a * rs12 + b2a * rs + b3a * rs32 + b4a * rs2); - doma = 2.0 * aa * (0.50 * b1a * rs12 + b2a * rs + 1.50 * b3a * - rs32 + 2.0 * b4a * rs2); + doma = 2.0 * aa * (0.50 * b1a * rs12 + b2a * rs + 1.50 * b3a * rs32 + 2.0 * b4a * rs2); ologa = log(1.0 + 1.00 / oma); alpha = 2.0 * aa * (1.0 + a1a * rs) * ologa; vpwca = + 2.0 * aa * (1.0 + 2.0 / 3.0 * a1a * rs) * ologa + 2.0 / 3.0 * aa * (1.0 + a1a * rs) * doma / (oma * (oma + 1.0)); - fz = (pow((1.0 + zeta) , (4.0 / 3.0)) + pow((1.0 - zeta) , (4.0 / - 3.0)) - 2.0) / (pow(2.0, (4.0 / 3.0)) - 2.0); - dfz = (pow((1.0 + zeta) , (1.0 / 3.0)) - pow((1.0 - zeta) , (1.0 / - 3.0))) * 4.0 / (3.0 * (pow(2.0 , (4.0 / 3.0)) - 2.0)); + fz = (pow((1.0 + zeta), (4.0 / 3.0)) + pow((1.0 - zeta), (4.0 / 3.0)) - 2.0) / (pow(2.0, (4.0 / 3.0)) - 2.0); + dfz = (pow((1.0 + zeta), (1.0 / 3.0)) - pow((1.0 - zeta), (1.0 / 3.0))) * 4.0 / (3.0 * (pow(2.0, (4.0 / 3.0)) - 2.0)); - ec = epwc + alpha * fz * (1.0 - zeta4) / fz0 + (epwcp - epwc) - * fz * zeta4; + ec = epwc + alpha * fz * (1.0 - zeta4) / fz0 + (epwcp - epwc) * fz * zeta4; - vcup = vpwc + vpwca * fz * (1.0 - zeta4) / fz0 - + (vpwcp - vpwc) * fz * zeta4 - + (alpha / fz0 * (dfz * (1.0 - zeta4) - 4.0 * fz * zeta3) - + (epwcp - epwc) * (dfz * zeta4 + 4.0 * fz * zeta3)) * (1.0 - zeta); + vcup = vpwc + vpwca * fz * (1.0 - zeta4) / fz0 + (vpwcp - vpwc) * fz * zeta4 + (alpha / fz0 * (dfz * (1.0 - zeta4) - 4.0 * fz * zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.0 * fz * zeta3)) * (1.0 - zeta); - vcdw = vpwc + vpwca * fz * (1.0 - zeta4) / fz0 - + (vpwcp - vpwc) * fz * zeta4 - - (alpha / fz0 * (dfz * (1.0 - zeta4) - 4.0 * fz * zeta3) - + (epwcp - epwc) * (dfz * zeta4 + 4.0 * fz * zeta3)) * (1.0 + zeta); - return; + vcdw = vpwc + vpwca * fz * (1.0 - zeta4) / fz0 + (vpwcp - vpwc) * fz * zeta4 - (alpha / fz0 * (dfz * (1.0 - zeta4) - 4.0 * fz * zeta3) + (epwcp - epwc) * (dfz * zeta4 + 4.0 * fz * zeta3)) * (1.0 + zeta); + return; } // J.P. Perdew and Y. Wang, PRB 45, 13244 (1992) -void XC_Functional::pz_spin( const double &rs, const double &zeta, - double &ec, double &vcup, double &vcdw) +void XC_Functional::pz_spin( + const double &rs, + const double &zeta, + double &ec, + double &vcup, + double &vcdw) { // ModuleBase::TITLE("XC_Functional","pz_spin"); - double ecu, vcu, ecp, vcp ; + double ecu = 0.0; + double vcu = 0.0; + double ecp = 0.0; + double vcp = 0.0; static const double p43 = 4.00 / 3.0; static const double third = 1.0 / 3.0; @@ -309,51 +364,52 @@ void XC_Functional::pz_spin( const double &rs, const double &zeta, XC_Functional::pz(rs, 0, ecu, vcu); // polarization contribution - XC_Functional::pz_polarized(rs, ecp, vcp); + XC_Functional::pz_polarized(rs, ecp, vcp); - double fz = (pow((1.00 + zeta) , p43) + pow((1.0 - zeta) , p43) - 2.0) / - (pow(2.0 , p43) - 2.0); - double dfz = p43 * (pow((1.00 + zeta) , third) - pow((1.0 - zeta) , third)) - / (pow(2.0, p43) - 2.0); + double fz = (pow((1.00 + zeta), p43) + pow((1.0 - zeta), p43) - 2.0) / (pow(2.0, p43) - 2.0); + double dfz = p43 * (pow((1.00 + zeta), third) - pow((1.0 - zeta), third)) / (pow(2.0, p43) - 2.0); - //correlation energy + //correlation energy ec = ecu + fz * (ecp - ecu); - - // spin up correlation potential + + // spin up correlation potential vcup = vcu + fz * (vcp - vcu) + (ecp - ecu) * dfz * (1.0 - zeta); - // spin down correlation potential + // spin down correlation potential vcdw = vcu + fz * (vcp - vcu) + (ecp - ecu) * dfz * (- 1.0 - zeta); return; } // J.P. Perdew and A. Zunger, PRB 23, 5048 (1981) -void XC_Functional::pz_polarized( const double &rs, double &ec, double &vc) +void XC_Functional::pz_polarized( + const double &rs, + double &ec, + double &vc) { // ModuleBase::TITLE("XC_Functional","pz_polarized"); - static const double a = 0.015550; - static const double b = -0.02690; - static const double c = 0.00070; - static const double d = -0.00480; - static const double gc = -0.08430; - static const double b1 = 1.39810; - static const double b2 = 0.26110; - - if (rs < 1.00) - { - double lnrs = log(rs); - ec = a * lnrs + b + c * rs * lnrs + d * rs; - vc = a * lnrs + (b - a / 3.0) + 2.0 / 3.0 * c * rs * lnrs + - (2.0 * d - c) / 3.0 * rs; - } - else - { - // interpolation formula - double rs12 = sqrt(rs); - double ox = 1.0 + b1 * rs12 + b2 * rs; - double dox = 1.0 + 7.0 / 6.0 * b1 * rs12 + 4.0 / 3.0 * b2 * rs; - ec = gc / ox; - vc = ec * dox / ox; - } - return; -} \ No newline at end of file + static const double a = 0.015550; + static const double b = -0.02690; + static const double c = 0.00070; + static const double d = -0.00480; + static const double gc = -0.08430; + static const double b1 = 1.39810; + static const double b2 = 0.26110; + + if (rs < 1.00) + { + double lnrs = log(rs); + ec = a * lnrs + b + c * rs * lnrs + d * rs; + vc = a * lnrs + (b - a / 3.0) + 2.0 / 3.0 * c * rs * lnrs + + (2.0 * d - c) / 3.0 * rs; + } + else + { + // interpolation formula + double rs12 = sqrt(rs); + double ox = 1.0 + b1 * rs12 + b2 * rs; + double dox = 1.0 + 7.0 / 6.0 * b1 * rs12 + 4.0 / 3.0 * b2 * rs; + ec = gc / ox; + vc = ec * dox / ox; + } + return; +} diff --git a/source/source_hamilt/module_xc/xc_funct_exch_gga.cpp b/source/source_hamilt/module_xc/xc_funct_exch_gga.cpp index e1f1c0cf27e..6eb8885410a 100644 --- a/source/source_hamilt/module_xc/xc_funct_exch_gga.cpp +++ b/source/source_hamilt/module_xc/xc_funct_exch_gga.cpp @@ -10,17 +10,32 @@ #include "xc_functional.h" -void XC_Functional::becke88(const double &rho, const double &grho, double &sx, double &v1x, double &v2x) +void XC_Functional::becke88( + const double &rho, + const double &grho, + double &sx, + double &v1x, + double &v2x) { //----------------------------------------------------------------------- // Becke exchange: A.D. Becke, PRA 38, 3098 (1988) // only gradient-corrected part, no Slater term included // - double beta, third, two13; + double beta = 0.0; + double third = 0.0; + double two13 = 0.0; beta = 0.00420; third = 1.0 / 3.0; two13 = 1.2599210498948730; - double rho13, rho43, xs, xs2, sa2b8, shm1, dd, dd2, ee; + double rho13 = 0.0; + double rho43 = 0.0; + double xs = 0.0; + double xs2 = 0.0; + double sa2b8 = 0.0; + double shm1 = 0.0; + double dd = 0.0; + double dd2 = 0.0; + double ee = 0.0; rho13 = pow(rho, third); rho43 = pow(rho13, 4); @@ -38,7 +53,12 @@ void XC_Functional::becke88(const double &rho, const double &grho, double &sx, d return; } // end subroutine becke88 -void XC_Functional::ggax(const double &rho, const double &grho, double &sx, double &v1x, double &v2x) +void XC_Functional::ggax( + const double &rho, + const double &grho, + double &sx, + double &v1x, + double &v2x) { //----------------------------------------------------------------------- // Perdew-Wang GGA (PW91), exchange part: @@ -51,8 +71,19 @@ void XC_Functional::ggax(const double &rho, const double &grho, double &sx, doub const double fp1 = -0.0192920212964260; const double fp2 = 0.1616204596739950; // fp1 = -3/(16 pi)*(3 pi^2)^(-1/3) - double rhom43, s, s2, s3, s4, exps, as, sa2b8, shm1, bs, das, - dbs, dls; + double rhom43 = 0.0; + double s = 0.0; + double s2 = 0.0; + double s3 = 0.0; + double s4 = 0.0; + double exps = 0.0; + double as = 0.0; + double sa2b8 = 0.0; + double shm1 = 0.0; + double bs = 0.0; + double das = 0.0; + double dbs = 0.0; + double dls = 0.0; rhom43 = pow(rho, (- 4.0 / 3.0)); s = fp2 * sqrt(grho) * rhom43; @@ -67,16 +98,21 @@ void XC_Functional::ggax(const double &rho, const double &grho, double &sx, doub das = (200.0 * exps - 2.0 * f5) * s; dbs = f1 * (shm1 + f2 * s / sa2b8) + 4.0 * f5 * s3; dls = (das / as - dbs / bs); - - sx = fp1 * grho * rhom43 * as / bs; - v1x = - 4.0 / 3.0 * sx / rho * (1.0 + s * dls); + + sx = fp1 * grho * rhom43 * as / bs; + v1x = - 4.0 / 3.0 * sx / rho * (1.0 + s * dls); v2x = fp1 * rhom43 * as / bs * (2.0 + s * dls); return; } //end subroutine ggax -void XC_Functional::pbex(const double &rho, const double &grho, const int &iflag, -double &sx, double &v1x, double &v2x) +void XC_Functional::pbex( + const double &rho, + const double &grho, + const int &iflag, + double &sx, + double &v1x, + double &v2x) { // PBE exchange (without Slater exchange): // iflag=0 J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996) @@ -93,8 +129,8 @@ double &sx, double &v1x, double &v2x) // n*ds/d(gn) // exchange energy LDA part // exchange energy gradient part - - // numerical coefficients (NB: c2=(3 pi^2)^(1/3) ) + + // numerical coefficients (NB: c2=(3 pi^2)^(1/3) ) const double third = 1.0 / 3.0; const double c1 = 0.750 / ModuleBase::PI; const double c2 = 3.0936677262801360; @@ -122,15 +158,20 @@ double &sx, double &v1x, double &v2x) const double dxunif = exunif * third; const double dfx1 = f2 * f2; const double dfx = 2.0 * mu[iflag] * s1 / dfx1; - - v1x = sx + dxunif * fx + exunif * dfx * ds; + + v1x = sx + dxunif * fx + exunif * dfx * ds; v2x = exunif * dfx * dsg / agrho; sx = sx * rho; - return; + return; } -void XC_Functional::optx(const double rho, const double grho, double &sx, double &v1x, double &v2x) +void XC_Functional::optx( + const double rho, + const double grho, + double &sx, + double &v1x, + double &v2x) { // OPTX, Handy et al. JCP 116, p. 5411 (2002) and refs. therein // Present release: Mauro Boero, Tsukuba, 10/9/2002 @@ -150,13 +191,18 @@ void XC_Functional::optx(const double rho, const double grho, double &sx, double double smal2 = 1.e-10; //.......coefficients and exponents.................... // parameter : - double o43 = 4.00 / 3.00, - two13 = 1.2599210498948730, - two53 = 3.1748021039363990, - gam = 0.0060, - a1cx = 0.97845711702844210, - a2 = 1.431690; - double gr, rho43, xa, gamx2, uden, uu; + double o43 = 4.00 / 3.00; + double two13 = 1.2599210498948730; + double two53 = 3.1748021039363990; + double gam = 0.0060; + double a1cx = 0.97845711702844210; + double a2 = 1.431690; + double gr = 0.0; + double rho43 = 0.0; + double xa = 0.0; + double gamx2 = 0.0; + double uden = 0.0; + double uu = 0.0; //.......OPTX in compact form.......................... if (rho <= small) @@ -167,7 +213,7 @@ void XC_Functional::optx(const double rho, const double grho, double &sx, double } else { - gr = (grho > smal2) ? grho : smal2; //max() + gr = (grho > smal2) ? grho : smal2; //max() rho43 = pow(rho, o43); xa = two13 * sqrt(gr) / rho43; gamx2 = gam * xa * xa; @@ -182,69 +228,104 @@ void XC_Functional::optx(const double rho, const double grho, double &sx, double return; } // end subroutine optx -void XC_Functional::wcx(const double &rho,const double &grho, double &sx, double &v1x, double &v2x) +void XC_Functional::wcx( + const double &rho, + const double &grho, + double &sx, + double &v1x, + double &v2x) { - double kf, agrho, s1, s2, es2, ds, dsg, exunif, fx; - // (3*pi2*|rho|)^(1/3) - // |grho| - // |grho|/(2*kf*|rho|) - // s^2 - // n*ds/dn - // n*ds/d(gn) - // exchange energy LDA part - // exchange energy gradient part - double dxunif, dfx, f1, f2, f3, dfx1, x1, x2, x3, dxds1, dxds2, dxds3; - // numerical coefficients (NB: c2=(3 pi^2)^(1/3) ) - double third, c1, c2, c5, teneightyone; // c6 + double kf = 0.0; + double agrho = 0.0; + double s1 = 0.0; + double s2 = 0.0; + double es2 = 0.0; + double ds = 0.0; + double dsg = 0.0; + double exunif = 0.0; + double fx = 0.0; + // (3*pi2*|rho|)^(1/3) + // |grho| + // |grho|/(2*kf*|rho|) + // s^2 + // n*ds/dn + // n*ds/d(gn) + // exchange energy LDA part + // exchange energy gradient part + double dxunif = 0.0; + double dfx = 0.0; + double f1 = 0.0; + double f2 = 0.0; + double f3 = 0.0; + double dfx1 = 0.0; + double x1 = 0.0; + double x2 = 0.0; + double x3 = 0.0; + double dxds1 = 0.0; + double dxds2 = 0.0; + double dxds3 = 0.0; + // numerical coefficients (NB: c2=(3 pi^2)^(1/3) ) + double third = 0.0; + double c1 = 0.0; + double c2 = 0.0; + double c5 = 0.0; + double teneightyone = 0.0; // c6 - third = 1.0/3.0; - c1 = 0.75 / ModuleBase::PI; - c2 = 3.093667726280136; - c5 = 4.0 * third; - teneightyone = 0.123456790123; - // parameters of the functional - double k, mu, cwc; - k = 0.804; - mu = 0.2195149727645171; - cwc = 0.00793746933516; - // - agrho = sqrt (grho); - kf = c2 * pow(rho,third); - dsg = 0.5 / kf; - s1 = agrho * dsg / rho; - s2 = s1 * s1; - es2 = exp(-s2); - ds = - c5 * s1; - // - // Energy - // - // x = 10/81 s^2 + (mu - 10/81) s^2 e^-s^2 + ln (1 + c s^4) - x1 = teneightyone * s2; - x2 = (mu - teneightyone) * s2 * es2; - x3 = log(1.0 + cwc * s2 * s2); - f1 = (x1 + x2 + x3) / k; - f2 = 1.0 + f1; - f3 = k / f2; - fx = k - f3; - exunif = - c1 * kf; - sx = exunif * fx; - // - // Potential - // - dxunif = exunif * third; - dfx1 = f2 * f2; - dxds1 = teneightyone; - dxds2 = (mu - teneightyone) * es2 * (1.0 - s2); - dxds3 = 2.0 * cwc * s2 / (1.0 + cwc * s2 *s2); - dfx = 2.0 * s1 * (dxds1 + dxds2 + dxds3) / dfx1; - v1x = sx + dxunif * fx + exunif * dfx * ds; - v2x = exunif * dfx * dsg / agrho; + third = 1.0/3.0; + c1 = 0.75 / ModuleBase::PI; + c2 = 3.093667726280136; + c5 = 4.0 * third; + teneightyone = 0.123456790123; + // parameters of the functional + double k = 0.0; + double mu = 0.0; + double cwc = 0.0; + k = 0.804; + mu = 0.2195149727645171; + cwc = 0.00793746933516; + // + agrho = sqrt (grho); + kf = c2 * pow(rho,third); + dsg = 0.5 / kf; + s1 = agrho * dsg / rho; + s2 = s1 * s1; + es2 = exp(-s2); + ds = - c5 * s1; + // + // Energy + // + // x = 10/81 s^2 + (mu - 10/81) s^2 e^-s^2 + ln (1 + c s^4) + x1 = teneightyone * s2; + x2 = (mu - teneightyone) * s2 * es2; + x3 = log(1.0 + cwc * s2 * s2); + f1 = (x1 + x2 + x3) / k; + f2 = 1.0 + f1; + f3 = k / f2; + fx = k - f3; + exunif = - c1 * kf; + sx = exunif * fx; + // + // Potential + // + dxunif = exunif * third; + dfx1 = f2 * f2; + dxds1 = teneightyone; + dxds2 = (mu - teneightyone) * es2 * (1.0 - s2); + dxds3 = 2.0 * cwc * s2 / (1.0 + cwc * s2 *s2); + dfx = 2.0 * s1 * (dxds1 + dxds2 + dxds3) / dfx1; + v1x = sx + dxunif * fx + exunif * dfx * ds; + v2x = exunif * dfx * dsg / agrho; - sx = sx * rho; - return; + sx = sx * rho; + return; } -void XC_Functional::becke88_spin(double rho, double grho, double &sx, double &v1x, double &v2x) +void XC_Functional::becke88_spin( + double rho, + double grho, + double &sx, + double &v1x, + double &v2x) { //------------------------------------------------------------------- // Becke exchange: A.D. Becke, PRA 38, 3098 (1988) - Spin polarized case @@ -258,11 +339,20 @@ void XC_Functional::becke88_spin(double rho, double grho, double &sx, double &v1 // output: first part of the potential // output: the second part of the potential - double beta, third; + double beta = 0.0; + double third = 0.0; // parameter : beta = 0.00420; third = 1.0 / 3.0; - double rho13, rho43, xs, xs2, sa2b8, shm1, dd, dd2, ee; + double rho13 = 0.0; + double rho43 = 0.0; + double xs = 0.0; + double xs2 = 0.0; + double sa2b8 = 0.0; + double shm1 = 0.0; + double dd = 0.0; + double dd2 = 0.0; + double ee = 0.0; rho13 = pow(rho, third); rho43 = pow(rho13, 4); diff --git a/source/source_hamilt/module_xc/xc_funct_exch_lda.cpp b/source/source_hamilt/module_xc/xc_funct_exch_lda.cpp index 9d6cb6d955b..5f0cc1ede20 100644 --- a/source/source_hamilt/module_xc/xc_funct_exch_lda.cpp +++ b/source/source_hamilt/module_xc/xc_funct_exch_lda.cpp @@ -11,29 +11,38 @@ #include "xc_functional.h" //Slater exchange with alpha=2/3 -void XC_Functional::slater(const double &rs, double &ex, double &vx) +void XC_Functional::slater( + const double &rs, + double &ex, + double &vx) { - // f = -9/8*(3/2pi)^(2/3) - const double f = -0.687247939924714e0; - const double alpha = 2.00 / 3.00; - ex = f * alpha / rs; - vx = 4.0 / 3.0 * f * alpha / rs; - return; + // f = -9/8*(3/2pi)^(2/3) + const double f = -0.687247939924714e0; + const double alpha = 2.00 / 3.00; + ex = f * alpha / rs; + vx = 4.0 / 3.0 * f * alpha / rs; + return; } //Slater exchange with alpha=1, corresponding to -1.374/r_s Ry //used to recover old results -void XC_Functional::slater1(const double &rs, double &ex, double &vx) +void XC_Functional::slater1( + const double &rs, + double &ex, + double &vx) { - const double f = -0.687247939924714e0; - const double alpha = 1.0; - ex = f * alpha / rs; - vx = 4.0 / 3.0 * f * alpha / rs; - return; + const double f = -0.687247939924714e0; + const double alpha = 1.0; + ex = f * alpha / rs; + vx = 4.0 / 3.0 * f * alpha / rs; + return; } // Slater exchange with alpha=2/3 and Relativistic exchange -void XC_Functional::slater_rxc(const double &rs, double &ex, double &vx) +void XC_Functional::slater_rxc( + const double &rs, + double &ex, + double &vx) { const double trd = 1.0 / 3.0; //const double ftrd = 4.0 / 3.0; @@ -56,8 +65,12 @@ void XC_Functional::slater_rxc(const double &rs, double &ex, double &vx) } // Slater exchange with alpha=2/3, spin-polarized case -void XC_Functional::slater_spin( const double &rho, const double &zeta, - double &ex, double &vxup, double &vxdw) +void XC_Functional::slater_spin( + const double &rho, + const double &zeta, + double &ex, + double &vxup, + double &vxdw) { const double f = - 1.107838149573033610; const double alpha = 2.00 / 3.00; @@ -65,10 +78,10 @@ void XC_Functional::slater_spin( const double &rho, const double &zeta, const double third = 1.0 / 3.0; const double p43 = 4.0 / 3.0; - double rho13 = pow(((1.0 + zeta) * rho) , third); + double rho13 = pow(((1.0 + zeta) * rho), third); double exup = f * alpha * rho13; vxup = p43 * f * alpha * rho13; - rho13 = pow(((1.0 - zeta) * rho) , third); + rho13 = pow(((1.0 - zeta) * rho), third); double exdw = f * alpha * rho13; vxdw = p43 * f * alpha * rho13; ex = 0.50 * ((1.0 + zeta) * exup + (1.0 - zeta) * exdw); @@ -77,18 +90,23 @@ void XC_Functional::slater_spin( const double &rho, const double &zeta, } // Slater exchange with alpha=2/3, spin-polarized case -void XC_Functional::slater1_spin( const double &rho, const double &zeta, double &ex, double &vxup, double &vxdw) +void XC_Functional::slater1_spin( + const double &rho, + const double &zeta, + double &ex, + double &vxup, + double &vxdw) { const double f = - 1.107838149573033610; - const double alpha = 1.00; - const double third = 1.0 / 3.0; - const double p43 = 4.0 / 3.0; + const double alpha = 1.00; + const double third = 1.0 / 3.0; + const double p43 = 4.0 / 3.0; // f = -9/8*(3/pi)^(1/3) - double rho13 = pow(((1.0 + zeta) * rho) , third); + double rho13 = pow(((1.0 + zeta) * rho), third); double exup = f * alpha * rho13; vxup = p43 * f * alpha * rho13; - rho13 = pow(((1.0 - zeta) * rho) , third); + rho13 = pow(((1.0 - zeta) * rho), third); double exdw = f * alpha * rho13; vxdw = p43 * f * alpha * rho13; ex = 0.50 * ((1.0 + zeta) * exup + (1.0 - zeta) * exdw); @@ -97,24 +115,30 @@ void XC_Functional::slater1_spin( const double &rho, const double &zeta, double } // end subroutine slater1_spin // Slater exchange with alpha=2/3, relativistic exchange case -void XC_Functional::slater_rxc_spin( const double &rho, const double &z, - double &ex, double &vxup, double &vxdw) +void XC_Functional::slater_rxc_spin( + const double &rho, + const double &z, + double &ex, + double &vxup, + double &vxdw) { if (rho <= 0.0) { - ex = vxup = vxdw = 0.0; + ex = 0.0; + vxup = 0.0; + vxdw = 0.0; return; } const double trd = 1.0 / 3.0; - const double ftrd = 4.0 / 3.0; - double tftm = pow(2.0, ftrd) - 2; + const double ftrd = 4.0 / 3.0; + double tftm = pow(2.0, ftrd) - 2; double a0 = pow((4 / (9 * ModuleBase::PI)), trd); double alp = 2 * trd; - double fz = (pow((1 + z), ftrd) + pow((1 - z), ftrd) - 2) / tftm; - double fzp = ftrd * (pow((1 + z), trd) - pow((1 - z), trd)) / tftm; + double fz = (pow((1 + z), ftrd) + pow((1 - z), ftrd) - 2) / tftm; + double fzp = ftrd * (pow((1 + z), trd) - pow((1 - z), trd)) / tftm; double rs = pow((3 / (4 * ModuleBase::PI * rho)), trd); double vxp = -3 * alp / (2 * ModuleBase::PI * a0 * rs); @@ -125,7 +149,7 @@ void XC_Functional::slater_rxc_spin( const double &rho, const double &z, double alb = log(beta + sb); vxp = vxp * (-0.5 + 1.5 * alb / (beta * sb)); exp = exp * (1.0- 1.5*( (beta*sb-alb) / (beta*beta) ) - * 1.5*( (beta*sb-alb) / (beta*beta) )); + * 1.5*( (beta*sb-alb) / (beta*beta) )); double x = (beta * sb - alb) / (beta * beta); diff --git a/source/source_hamilt/module_xc/xc_funct_hcth.cpp b/source/source_hamilt/module_xc/xc_funct_hcth.cpp index d7c3333c99f..2e655df9532 100644 --- a/source/source_hamilt/module_xc/xc_funct_hcth.cpp +++ b/source/source_hamilt/module_xc/xc_funct_hcth.cpp @@ -2,7 +2,12 @@ #include "xc_functional.h" -void XC_Functional::hcth(const double rho, const double grho, double &sx, double &v1x, double &v2x) +void XC_Functional::hcth( + const double rho, + const double grho, + double &sx, + double &v1x, + double &v2x) { // =============================================================== // HCTH/120, JCP 109, p. 6264 (1998) @@ -24,12 +29,50 @@ void XC_Functional::hcth(const double rho, const double grho, double &sx, double double o34 = 4.00 / 3.00; double fr83 = 8.0 / 3.0; //double cg0[6], cg1[6], caa[6], cab[6], cx[6]; - double cg0[7], cg1[7], caa[7], cab[7], cx[7]; //mohan modify 2007-10-13 - double r3q2, r3pi, gr, rho_o3, rho_o34, xa, xa2, ra, rab, - dra_drho, drab_drho, g, dg, era1, dera1_dra, erab0, derab0_drab, - ex, dex_drho, uaa, uab, ux, ffaa, ffab, dffaa_drho, dffab_drho, - denaa, denab, denx, f83rho, bygr, gaa, gab, gx, taa, tab, txx, - dgaa_drho, dgab_drho, dgx_drho, dgaa_dgr, dgab_dgr, dgx_dgr; + double cg0[7], cg1[7], caa[7], cab[7], cx[7]; //mohan modify 2007-10-13 + double r3q2 = 0.0; + double r3pi = 0.0; + double gr = 0.0; + double rho_o3 = 0.0; + double rho_o34 = 0.0; + double xa = 0.0; + double xa2 = 0.0; + double ra = 0.0; + double rab = 0.0; + double dra_drho = 0.0; + double drab_drho = 0.0; + double g = 0.0; + double dg = 0.0; + double era1 = 0.0; + double dera1_dra = 0.0; + double erab0 = 0.0; + double derab0_drab = 0.0; + double ex = 0.0; + double dex_drho = 0.0; + double uaa = 0.0; + double uab = 0.0; + double ux = 0.0; + double ffaa = 0.0; + double ffab = 0.0; + double dffaa_drho = 0.0; + double dffab_drho = 0.0; + double denaa = 0.0; + double denab = 0.0; + double denx = 0.0; + double f83rho = 0.0; + double bygr = 0.0; + double gaa = 0.0; + double gab = 0.0; + double gx = 0.0; + double taa = 0.0; + double tab = 0.0; + double txx = 0.0; + double dgaa_drho = 0.0; + double dgab_drho = 0.0; + double dgx_drho = 0.0; + double dgaa_dgr = 0.0; + double dgab_dgr = 0.0; + double dgx_dgr = 0.0; r3q2 = std::pow(2.0, (-o3)); r3pi = std::pow((3.0 / ModuleBase::PI), o3); @@ -126,21 +169,30 @@ void XC_Functional::hcth(const double rho, const double grho, double &sx, double return; } //end subroutine hcth -void XC_Functional::pwcorr(const double r, const double c[], double &g, double &dg) +void XC_Functional::pwcorr( + const double r, + const double c[], + double &g, + double &dg) { // USE kinds // implicit none // real(kind=DP) :: r, g, dg, c(6) - double r12, r32, r2, rb, drb, sb; + double r12 = 0.0; + double r32 = 0.0; + double r2 = 0.0; + double rb = 0.0; + double drb = 0.0; + double sb = 0.0; r12 = sqrt(r); r32 = r * r12; r2 = r * r; - rb = c[3] * r12 + c[4] * r + c[5] * r32 + c[6] * r2; //c[i] i=0--5; + rb = c[3] * r12 + c[4] * r + c[5] * r32 + c[6] * r2; //c[i] i=0--5; sb = 1.00 + 1.00 / (2.00 * c[1] * rb); g = -2.00 * c[1] * (1.00 + c[2] * r) * log(sb); drb = c[3] / (2.00 * r12) + c[4] + 1.50 * c[5] * r12 + 2.00 * c[6] * r; dg = (1.00 + c[2] * r) * drb / (rb * rb * sb) - 2.00 * c[1] * c[2] * log(sb); return; -} //end subroutine pwcorr \ No newline at end of file +} //end subroutine pwcorr diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index a311d28d8a0..7322ea3b46f 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -22,10 +22,10 @@ class XC_Functional { - public: + public: - XC_Functional(); - ~XC_Functional(); + XC_Functional(); + ~XC_Functional(); //------------------- // subroutines, grouped according to the file they are in: @@ -42,12 +42,12 @@ class XC_Functional // NOTE : it is only used for nspin = 1 and 2, the nspin = 4 case is treated in v_xc // 3. v_xc_meta : which takes rho and tau as input, and v_xc as output - // compute the exchange-correlation energy - // [etxc, vtxc, v] = v_xc(...) - static std::tuple v_xc( - const int &nrxx, // number of real-space grid - const Charge* const chr, - const UnitCell *ucell); // charge density + // compute the exchange-correlation energy + // [etxc, vtxc, v] = v_xc(...) + static std::tuple v_xc( + const int &nrxx, // number of real-space grid + const Charge* const chr, + const UnitCell *ucell); // charge density //------------------- // xc_functional.cpp @@ -56,49 +56,53 @@ class XC_Functional // This file contains subroutines for setting the functional // it includes 4 subroutines: // 1. get_func_type : which returns the type of functional (func_type): -// 0 = none; 1 = lda; 2 = gga; 3 = mgga; 4 = hybrid lda/gga; 5 = hybrid mgga +// 0 = none; 1 = lda; 2 = gga; 3 = mgga; 4 = hybrid lda/gga; 5 = hybrid mgga // 2. set_xc_type : sets the value of: -// 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. +// 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. static int get_func_type() { return func_type; }; + static void set_xc_type(const std::string xc_func_in); // For hybrid functional static void set_hybrid_alpha(const double alpha_in); + static double get_hybrid_alpha() { return hybrid_alpha; }; + static bool get_ked_flag() { return ked_flag; }; + /// Usually in exx caculation, the first SCF loop should be converged with PBE static void set_xc_first_loop(const UnitCell& ucell); - static std::string output_info(); + static std::string output_info(); - private: + private: - static std::vector func_id; // libxc id of functional - static int func_type; //0:none, 1:lda, 2:gga, 3:mgga, 4:hybrid lda/gga, 5:hybrid mgga + static std::vector func_id; // libxc id of functional + static int func_type; //0:none, 1:lda, 2:gga, 3:mgga, 4:hybrid lda/gga, 5:hybrid mgga static bool ked_flag; // whether the functional has kinetic energy density static bool use_libxc; // exx_hybrid_alpha for mixing exx in hybrid functional: static double hybrid_alpha; - // added by jghan, 2024-07-07 - // as a scaling factor for different xc-functionals - static std::map scaling_factor_xc; + // added by jghan, 2024-07-07 + // as a scaling factor for different xc-functionals + static std::map scaling_factor_xc; - public: - static std::vector get_func_id() { return func_id; } + public: + static std::vector get_func_id() { return func_id; } //------------------- // xc_functional_wrapper_xc.cpp @@ -124,14 +128,18 @@ class XC_Functional // on the entire grid. I'm having xc_spin_libxc because v_xc_libxc // does not support nspin = 4. - public: + public: - // LDA - static void xc(const double &rho, double &exc, double &vxc); + // LDA + static void xc(const double &rho, double &exc, double &vxc); - // LSDA - static void xc_spin(const double &rho, const double &zeta, - double &exc, double &vxcup, double &vxcdw); + // LSDA + static void xc_spin( + const double &rho, + const double &zeta, + double &exc, + double &vxcup, + double &vxcdw); //------------------- // xc_functional_wrapper_gcxc.cpp @@ -147,16 +155,34 @@ class XC_Functional // LIBXC, and the reason for not having gcxc_libxc is explained // in the NOTE in the comment for xc_functional_wrapper_wc.cpp part - // GGA - static void gcxc(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, - double &sx, double &v1xup, double &v1xdw, double &v2xup, - double &v2xdw); - static void gcc_spin(double rho, double &zeta, double grho, double &sc, - double &v1cup, double &v1cdw, double &v2c); + // GGA + static void gcxc( + 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, + double &sx, + double &v1xup, + double &v1xdw, + double &v2xup, + double &v2xdw); + + static void gcc_spin( + double rho, + double &zeta, + double grho, + double &sc, + double &v1cup, + double &v1cdw, + double &v2c); //------------------- // xc_functional_gradcorr.cpp @@ -166,43 +192,51 @@ class XC_Functional // it contains 5 subroutines: // 1. gradcorr, which calculates gradient correction // 2. grad_wfc, which calculates gradient of wavefunction -// it is used in stress_func_mgga.cpp +// it is used in stress_func_mgga.cpp // 3. grad_rho, which calculates gradient of density // 4. grad_dot, which calculates divergence of something // 5. noncolin_rho, which diagonalizes the spin density matrix // and gives the spin up and spin down components of the charge. - static void gradcorr(double& etxc, - double& vtxc, - ModuleBase::matrix& v, - const Charge* const chr, - ModulePW::PW_Basis* rhopw, - const UnitCell* ucell, - std::vector& stress_gga, - const bool is_stress = false); - template ::type> - static void grad_wfc( - const int ik, - const Real tpiba, - const ModulePW::PW_Basis_K* wfc_basis, - const T* rhog, - T* grad); - static void grad_rho(const std::complex* rhog, - ModuleBase::Vector3* gdr, - const ModulePW::PW_Basis* rho_basis, - const double tpiba); - static void grad_dot(const ModuleBase::Vector3* h, + static void gradcorr( + double& etxc, + double& vtxc, + ModuleBase::matrix& v, + const Charge* const chr, + ModulePW::PW_Basis* rhopw, + const UnitCell* ucell, + std::vector& stress_gga, + const bool is_stress = false); + template ::type> + + static void grad_wfc( + const int ik, + const Real tpiba, + const ModulePW::PW_Basis_K* wfc_basis, + const T* rhog, + T* grad); + + static void grad_rho( + const std::complex* rhog, + ModuleBase::Vector3* gdr, + const ModulePW::PW_Basis* rho_basis, + const double tpiba); + + static void grad_dot( + const ModuleBase::Vector3* h, double* dh, const ModulePW::PW_Basis* rho_basis, const double tpiba); - static void noncolin_rho(double* rhoout1, - double* rhoout2, - double* seg, - const double* const* const rho, - const int nrxx, - const double* ux_, - const bool lsign_); + + static void noncolin_rho( + double* rhoout1, + double* rhoout2, + double* seg, + const double* const* const rho, + const int nrxx, + const double* ux_, + const bool lsign_); //------------------- // xc_funct_exch_lda.cpp @@ -218,18 +252,30 @@ class XC_Functional // 2. slater1_spin // 3. slater_rxc_spin - // For LDA exchange energy - static void slater(const double &rs, double &ex, double &vx); - static void slater1(const double &rs, double &ex, double &vx); - static void slater_rxc(const double &rs, double &ex, double &vx); - - // For LSDA exchange energy - static void slater_spin( const double &rho, const double &zeta, - double &ex, double &vxup, double &vxdw); - static void slater1_spin( const double &rho, const double &zeta, - double &ex, double &vxup, double &vxdw); - static void slater_rxc_spin( const double &rho, const double &z, - double &ex, double &vxup, double &vxdw); + // For LDA exchange energy + static void slater(const double &rs, double &ex, double &vx); + static void slater1(const double &rs, double &ex, double &vx); + static void slater_rxc(const double &rs, double &ex, double &vx); + + // For LSDA exchange energy + static void slater_spin( + const double &rho, + const double &zeta, + double &ex, + double &vxup, + double &vxdw); + static void slater1_spin( + const double &rho, + const double &zeta, + double &ex, + double &vxup, + double &vxdw); + static void slater_rxc_spin( + const double &rho, + const double &z, + double &ex, + double &vxup, + double &vxdw); //------------------- // xc_funct_corr_lda.cpp @@ -248,21 +294,31 @@ class XC_Functional // 1. pw_spin // 2. pz_spin, which calls pz_polarized - // For LDA correlation energy - static void pw(const double &rs, const int &iflag, double &ec, double &vc); - static void pz(const double &rs, const int &iflag, double &ec, double &vc); - static void lyp(const double &rs, double &ec, double &vc); - static void vwn(const double &rs, double &ec, double &vc); - static void wigner(const double &rs, double &ec, double &vc); - static void hl(const double &rs, double &ec, double &vc); - static void gl(const double &rs, double &ec, double &vc); - - // For LSDA correlation energy - static void pw_spin( const double &rs, const double &zeta, - double &ec, double &vcup, double &vcdw); - static void pz_spin( const double &rs, const double &zeta, - double &ec, double &vcup, double &vcdw); - static void pz_polarized( const double &rs, double &ec, double &vc); + // For LDA correlation energy + static void pw(const double &rs, const int &iflag, double &ec, double &vc); + static void pz(const double &rs, const int &iflag, double &ec, double &vc); + static void lyp(const double &rs, double &ec, double &vc); + static void vwn(const double &rs, double &ec, double &vc); + static void wigner(const double &rs, double &ec, double &vc); + static void hl(const double &rs, double &ec, double &vc); + static void gl(const double &rs, double &ec, double &vc); + + // For LSDA correlation energy + static void pw_spin( + const double &rs, + const double &zeta, + double &ec, + double &vcup, + double &vcdw); + + static void pz_spin( + const double &rs, + const double &zeta, + double &ec, + double &vcup, + double &vcdw); + + static void pz_polarized(const double &rs, double &ec, double &vc); //------------------- // xc_funct_exch_gga.cpp @@ -278,15 +334,43 @@ class XC_Functional // And some of their spin polarized counterparts: // 1. becke88_spin - static void becke88(const double &rho, const double &grho, double &sx, double &v1x, double &v2x); - static void ggax(const double &rho, const double &grho, double &sx, double &v1x, double &v2x); - static void pbex(const double &rho, const double &grho, const int &iflag, - double &sx, double &v1x, double &v2x); - static void optx(const double rho, const double grho, double &sx, double &v1x, double &v2x); - static void wcx(const double &rho,const double &grho, double &sx, double &v1x, double &v2x); - - static void becke88_spin(double rho, double grho, double &sx, double &v1x, - double &v2x); + static void becke88( + const double &rho, + const double &grho, + double &sx, + double &v1x, + double &v2x); + + static void ggax( + const double &rho, + const double &grho, + double &sx, + double &v1x, + double &v2x); + + static void pbex( + const double &rho, + const double &grho, + const int &iflag, + double &sx, + double &v1x, + double &v2x); + + static void optx(const double rho, const double grho, double &sx, double &v1x, double &v2x); + + static void wcx( + const double &rho, + const double &grho, + double &sx, + double &v1x, + double &v2x); + + static void becke88_spin( + double rho, + double grho, + double &sx, + double &v1x, + double &v2x); //------------------- // xc_funct_corr_gga.cpp @@ -303,18 +387,51 @@ class XC_Functional // 2. ggac_spin // 3. pbec_spin - static void perdew86(const double rho, const double grho, double &sc, double &v1c, double &v2c); - static void ggac(const double &rho,const double &grho, double &sc, double &v1c, double &v2c); - static void pbec(const double &rho, const double &grho, const int &flag, - double &sc, double &v1c, double &v2c); - static void glyp(const double &rho, const double &grho, double &sc, double &v1c, double &v2c); - - static void perdew86_spin(double rho, double zeta, double grho, double &sc, - double &v1cup, double &v1cdw, double &v2c); - //static void ggac_spin(double rho, double zeta, double grho, double &sc, - // double &v1cup, double &v1cdw, double &v2c); - static void pbec_spin(double rho, double zeta, double grho, const int &flag, double &sc, - double &v1cup, double &v1cdw, double &v2c); + static void perdew86(const double rho, const double grho, double &sc, double &v1c, double &v2c); + + static void ggac( + const double &rho, + const double &grho, + double &sc, + double &v1c, + double &v2c); + + static void pbec( + const double &rho, + const double &grho, + const int &flag, + double &sc, + double &v1c, + double &v2c); + + static void glyp( + const double &rho, + const double &grho, + double &sc, + double &v1c, + double &v2c); + + static void perdew86_spin( + double rho, + double zeta, + double grho, + double &sc, + double &v1cup, + double &v1cdw, + double &v2c); + + //static void ggac_spin(double rho, double zeta, double grho, double &sc, + // double &v1cup, double &v1cdw, double &v2c); + + static void pbec_spin( + double rho, + double zeta, + double grho, + const int &flag, + double &sc, + double &v1cup, + double &v1cdw, + double &v2c); //------------------- // xc_funct_hcth.cpp @@ -322,8 +439,9 @@ class XC_Functional // This file contains realizations of the HCTH GGA functional // hcth calls pwcorr - static void hcth(const double rho, const double grho, double &sx, double &v1x, double &v2x); - static void pwcorr(const double r, const double c[], double &g, double &dg); + static void hcth(const double rho, const double grho, double &sx, double &v1x, double &v2x); + + static void pwcorr(const double r, const double c[], double &g, double &dg); }; diff --git a/source/source_hamilt/module_xc/xc_functional_libxc.h b/source/source_hamilt/module_xc/xc_functional_libxc.h index 6b80d4ffb82..d23dbbb2340 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc.h +++ b/source/source_hamilt/module_xc/xc_functional_libxc.h @@ -25,7 +25,7 @@ namespace XC_Functional_Libxc // 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); + extern std::pair> set_xc_type_libxc(const std::string& xc_func_in); /** * @brief instantiate the XC functional by its ID, and set the external parameters if provided. @@ -43,8 +43,9 @@ namespace XC_Functional_Libxc * xc_corr_ext in the input file. The expected format would be an XC ID * followed by a list of parameters. */ - extern 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); @@ -53,7 +54,7 @@ namespace XC_Functional_Libxc // xc_functional_libxc_vxc.cpp //------------------- - extern std::tuple v_xc_libxc( + 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 @@ -62,7 +63,7 @@ namespace XC_Functional_Libxc const std::map* scaling_factor = nullptr); // added by jghan, 2024-10-10 // for mGGA functional - extern std::tuple v_xc_meta( + 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 @@ -117,7 +118,7 @@ namespace XC_Functional_Libxc std::vector exc); // converting vtxc and v from vrho and vsigma (libxc=>abacus) - extern std::pair convert_vtxc_v( + extern std::pair convert_vtxc_v( const xc_func_type &func, const int nspin, const std::size_t nrxx, @@ -153,8 +154,11 @@ namespace XC_Functional_Libxc extern void xc_spin_libxc( const std::vector &func_id, - const double &rhoup, const double &rhodw, - double &exc, double &vxcup, double &vxcdw); + const double &rhoup, + const double &rhodw, + double &exc, + double &vxcup, + double &vxcdw); //------------------- @@ -164,15 +168,25 @@ namespace XC_Functional_Libxc // 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); + 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); + 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); //------------------- @@ -182,25 +196,44 @@ namespace XC_Functional_Libxc // 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); + const double &rho, + const double &grho, + const double &atau, + double &sxc, + double &v1xc, + double &v2xc, + double &v3xc); 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, + const double &rho, + const double &grho, + const double &atau, + double &sxc, + double &v1xc, + double &v2xc, + double &v3xc, const double &hybrid_alpha); 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); + 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 +#endif // XC_FUNCTIONAL_LIBXC_H diff --git a/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp b/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp index 0a9bffbc426..255b6966f39 100644 --- a/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp +++ b/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp @@ -12,117 +12,171 @@ #include "xc_functional.h" #include -void XC_Functional::xc(const double &rho, double &exc, double &vxc) +void XC_Functional::xc( + const double &rho, + double &exc, + double &vxc) { - double third = 1.0 / 3.0; - double pi34 = 0.6203504908994e0 ; // pi34=(3/4pi)^(1/3) - double rs = 0.0; - double e,v; - - exc = vxc = 0.00; - - rs = pi34 / std::pow(rho, third); + double third = 1.0 / 3.0; + double pi34 = 0.6203504908994e0; // pi34=(3/4pi)^(1/3) + double rs = 0.0; + double e = 0.0; + double v = 0.0; - for(int id : func_id) + exc = 0.00; + vxc = 0.00; + + rs = pi34 / std::pow(rho, third); + + for (int id : func_id) { - switch( id ) + switch (id) { // Exchange functionals containing slater exchange - case XC_LDA_X: case XC_GGA_X_PBE: case XC_GGA_X_PBE_R: case XC_GGA_X_PBE_SOL: - case XC_GGA_X_WC: case XC_GGA_X_B88: case XC_GGA_X_PW91: - // SLA,PBX,rPBX,PBXsol,WC,B88,PW91_X - XC_Functional::slater(rs, e, v);break; + case XC_LDA_X: + case XC_GGA_X_PBE: + case XC_GGA_X_PBE_R: + case XC_GGA_X_PBE_SOL: + case XC_GGA_X_WC: + case XC_GGA_X_B88: + case XC_GGA_X_PW91: + // SLA,PBX,rPBX,PBXsol,WC,B88,PW91_X + XC_Functional::slater(rs, e, v); + break; // Exchange functionals containing attenuated slater exchange case XC_HYB_GGA_XC_PBEH: - // PBE0 - double ex, vx, ec, vc; + { + // PBE0 + double ex = 0.0; + double vx = 0.0; + double ec = 0.0; + double vc = 0.0; XC_Functional::slater(rs, ex, vx); - ex *= (1 - XC_Functional::hybrid_alpha); + ex *= (1 - XC_Functional::hybrid_alpha); vx *= (1 - XC_Functional::hybrid_alpha); XC_Functional::pw(rs, 0, ec, vc); e = ex + ec; v = vx + vc; break; + } // Correlation functionals containing PW correlation - case XC_GGA_C_PBE: case XC_GGA_C_PW91: case XC_LDA_C_PW: case XC_GGA_C_PBE_SOL: - // PBC,PW91,PWLDA - XC_Functional::pw(rs, 0, e, v);break; + case XC_GGA_C_PBE: + case XC_GGA_C_PW91: + case XC_LDA_C_PW: + case XC_GGA_C_PBE_SOL: + // PBC,PW91,PWLDA + XC_Functional::pw(rs, 0, e, v); + break; // Correlation functionals containing PZ correlation - case XC_LDA_C_PZ: case XC_GGA_C_P86: - // PZ,P86 - XC_Functional::pz(rs, 0, e, v);break; + case XC_LDA_C_PZ: + case XC_GGA_C_P86: + // PZ,P86 + XC_Functional::pz(rs, 0, e, v); + break; // Correlation functionals containing LYP correlation case XC_GGA_C_LYP: - // BLYP - XC_Functional::lyp(rs, e, v);break; - + // BLYP + XC_Functional::lyp(rs, e, v); + break; + default: - e = v = 0.0; + { + e = 0.0; + v = 0.0; + } } exc += e; vxc += v; } - return; + return; } -void XC_Functional::xc_spin(const double &rho, const double &zeta, - double &exc, double &vxcup, double &vxcdw) +void XC_Functional::xc_spin( + const double &rho, + const double &zeta, + double &exc, + double &vxcup, + double &vxcdw) { - static const double small = 1.e-10; - double e, vup, vdw; - exc = vxcup = vxcdw = 0.0; - - static const double third = 1.0 / 3.0; - static const double pi34 = 0.62035049089940; - const double rs = pi34 / pow(rho, third);//wigner_sitz_radius; - - for(int id : func_id) + static const double small = 1.e-10; + double e = 0.0; + double vup = 0.0; + double vdw = 0.0; + exc = 0.0; + vxcup = 0.0; + vxcdw = 0.0; + + static const double third = 1.0 / 3.0; + static const double pi34 = 0.62035049089940; + const double rs = pi34 / pow(rho, third); //wigner_sitz_radius; + + for (int id : func_id) { - switch( id ) + switch (id) { // Exchange functionals containing slater exchange - case XC_LDA_X: case XC_GGA_X_PBE: case XC_GGA_X_PBE_R: case XC_GGA_X_PBE_SOL: - case XC_GGA_X_WC: case XC_GGA_X_B88: case XC_GGA_X_PW91: - // SLA,PBX,rPBX,PBXsol,WC,B88,PW91_X - XC_Functional::slater_spin(rho, zeta, e, vup, vdw); break; - + case XC_LDA_X: + case XC_GGA_X_PBE: + case XC_GGA_X_PBE_R: + case XC_GGA_X_PBE_SOL: + case XC_GGA_X_WC: + case XC_GGA_X_B88: + case XC_GGA_X_PW91: + // SLA,PBX,rPBX,PBXsol,WC,B88,PW91_X + XC_Functional::slater_spin(rho, zeta, e, vup, vdw); + break; + // Exchange functionals containing attenuated slater exchange case XC_HYB_GGA_XC_PBEH: - // PBE0 - double ex, vupx, vdwx, ec, vupc, vdwc; + { + // PBE0 + double ex = 0.0; + double vupx = 0.0; + double vdwx = 0.0; + double ec = 0.0; + double vupc = 0.0; + double vdwc = 0.0; XC_Functional::slater_spin(rho, zeta, ex, vupx, vdwx); - ex *= (1.0 - XC_Functional::hybrid_alpha); - vupx *= (1.0 - XC_Functional::hybrid_alpha); + ex *= (1.0 - XC_Functional::hybrid_alpha); + vupx *= (1.0 - XC_Functional::hybrid_alpha); vdwx *= (1.0 - XC_Functional::hybrid_alpha); XC_Functional::pw_spin(rs, zeta, ec, vupc, vdwc); e = ex + ec; vup = vupx + vupc; vdw = vdwx + vdwc; break; + } // Correlation functionals containing PZ correlation - case XC_LDA_C_PZ: case XC_GGA_C_P86: - // PZ,P86 - XC_Functional::pz_spin(rs, zeta, e, vup, vdw); break; + case XC_LDA_C_PZ: + case XC_GGA_C_P86: + // PZ,P86 + XC_Functional::pz_spin(rs, zeta, e, vup, vdw); + break; // Correlation functionals containing PW correlationtests/integrate/101_PW_OU_pseudopot - case XC_GGA_C_PBE: case XC_GGA_C_PBE_SOL: case XC_LDA_C_PW: - // PBC,PBCsol - XC_Functional::pw_spin(rs, zeta, e, vup, vdw); break; + case XC_GGA_C_PBE: + case XC_GGA_C_PBE_SOL: + case XC_LDA_C_PW: + // PBC,PBCsol + XC_Functional::pw_spin(rs, zeta, e, vup, vdw); + break; // Cases that are only realized in LIBXC default: - throw std::domain_error("functional unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); break; - + { + throw std::domain_error("functional unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__)); + break; + } } exc += e; vxcup += vup; vxcdw += vdw; - } - return; -} \ No newline at end of file + } + return; +} From 025cbf5230e51428bec5098e704b32920d8470ff Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 30 May 2026 07:15:17 +0800 Subject: [PATCH 3/6] Refactor XC functional files for code style compliance - Format function parameters to one per line - Initialize all variables when declared - Separate variable declarations onto individual lines - Add braces for single-line if/for statements - Improve code alignment and readability --- .../source_hamilt/module_xc/xc_functional.cpp | 98 +- .../module_xc/xc_functional_gradcorr.cpp | 1328 +++++++++-------- .../module_xc/xc_functional_vxc.cpp | 31 +- .../module_xc/xc_functional_wrapper_gcxc.cpp | 291 +++- .../module_xc/xc_functional_wrapper_xc.cpp | 17 +- 5 files changed, 998 insertions(+), 767 deletions(-) diff --git a/source/source_hamilt/module_xc/xc_functional.cpp b/source/source_hamilt/module_xc/xc_functional.cpp index ca88899b3ea..acef0f6615f 100644 --- a/source/source_hamilt/module_xc/xc_functional.cpp +++ b/source/source_hamilt/module_xc/xc_functional.cpp @@ -71,13 +71,13 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) scaling_factor_xc.clear(); // added by jghan, 2024-07-07 std::string xc_func = xc_func_in; std::transform(xc_func.begin(), xc_func.end(), xc_func.begin(), (::toupper)); - if( xc_func == "LDA" || xc_func == "PZ" || xc_func == "SLAPZNOGXNOGC") //SLA+PZ - { + if( xc_func == "LDA" || xc_func == "PZ" || xc_func == "SLAPZNOGXNOGC") //SLA+PZ + { func_id.push_back(XC_LDA_X); func_id.push_back(XC_LDA_C_PZ); func_type = 1; use_libxc = false; - } + } else if (xc_func == "PWLDA") { func_id.push_back(XC_LDA_X); @@ -85,84 +85,84 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_type = 1; use_libxc = false; } - else if ( xc_func == "PBE" || xc_func == "SLAPWPBXPBC") //PBX+PBC - { + else if ( xc_func == "PBE" || xc_func == "SLAPWPBXPBC") //PBX+PBC + { func_id.push_back(XC_GGA_X_PBE); func_id.push_back(XC_GGA_C_PBE); func_type = 2; use_libxc = false; - } - else if ( xc_func == "PBESOL") //PBX_S+PBC_S - { + } + else if ( xc_func == "PBESOL") //PBX_S+PBC_S + { func_id.push_back(XC_GGA_X_PBE_SOL); func_id.push_back(XC_GGA_C_PBE_SOL); func_type = 2; use_libxc = false; - } - else if( xc_func == "REVPBE" ) //PBX_r+PBC - { - func_id.push_back(XC_GGA_X_PBE_R); + } + else if( xc_func == "REVPBE" ) //PBX_r+PBC + { + func_id.push_back(XC_GGA_X_PBE_R); func_id.push_back(XC_GGA_C_PBE); func_type = 2; use_libxc = false; - } - else if ( xc_func == "WC") //WC+PBC - { + } + else if ( xc_func == "WC") //WC+PBC + { func_id.push_back(XC_GGA_X_WC); func_id.push_back(XC_GGA_C_PBE); func_type = 2; use_libxc = false; - } - else if ( xc_func == "BLYP") //B88+LYP - { + } + else if ( xc_func == "BLYP") //B88+LYP + { func_id.push_back(XC_GGA_X_B88); func_id.push_back(XC_GGA_C_LYP); func_type = 2; use_libxc = false; - } - else if ( xc_func == "BP") //B88+P86 - { + } + else if ( xc_func == "BP") //B88+P86 + { func_id.push_back(XC_GGA_X_B88); func_id.push_back(XC_GGA_C_P86); func_type = 2; use_libxc = false; - } - else if ( xc_func == "PW91") //PW91_X+PW91_C - { + } + 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 - { + } + else if ( xc_func == "HCTH") //HCTH_X+HCTH_C + { func_id.push_back(XC_GGA_X_HCTH_A); func_id.push_back(XC_GGA_C_HCTH_A); func_type = 2; use_libxc = false; - } - else if ( xc_func == "OLYP") //OPTX+LYP - { + } + else if ( xc_func == "OLYP") //OPTX+LYP + { func_id.push_back(XC_GGA_X_OPTX); func_id.push_back(XC_GGA_C_LYP); func_type = 2; use_libxc = false; - } + } #ifdef USE_LIBXC - else if ( xc_func == "SCAN") - { + else if ( xc_func == "SCAN") + { func_id.push_back(XC_MGGA_X_SCAN); func_id.push_back(XC_MGGA_C_SCAN); func_type = 3; use_libxc = true; - } + } else if ( xc_func == "SCAN0") - { + { func_id.push_back(XC_MGGA_X_SCAN); func_id.push_back(XC_MGGA_C_SCAN); func_type = 5; use_libxc = true; - } + } else if( xc_func == "LC_PBE") { func_id.push_back(XC_HYB_GGA_XC_LC_PBEOP); @@ -199,12 +199,12 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_type = 4; use_libxc = false; } - else if( xc_func == "PBE0") - { + else if( xc_func == "PBE0") + { func_id.push_back(XC_HYB_GGA_XC_PBEH); func_type = 4; use_libxc = false; - } + } else if( xc_func == "OPT_ORB" || xc_func == "NONE" || xc_func == "NOX+NOC") { // not doing anything @@ -230,7 +230,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) use_libxc = true; } else if( xc_func == "CWP22") - { + { // BLYP_XC_lr = -BLYP_XC_sr + BLYP_XC, the realization of it is in v_xc_libxc() function, xc_functional_libxc_vxc.cpp func_id.push_back(XC_GGA_X_ITYH); // short-range of B88_X, id=529 func_id.push_back(XC_GGA_C_LYPR); // short-range of LYP_C, id=624 @@ -246,8 +246,8 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) func_type = 4; use_libxc = true; } - else if( xc_func == "BLYP_LR") - { + else if ( xc_func == "BLYP_LR") + { // BLYP_XC_lr = -BLYP_XC_sr + BLYP_XC, the realization of it is in v_xc_libxc() function, xc_functional_libxc_vxc.cpp func_id.push_back(XC_GGA_X_ITYH); // short-range of B88_X, id=529 func_id.push_back(XC_GGA_C_LYPR); // short-range of LYP_C, id=624 @@ -274,7 +274,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) { #ifdef USE_LIBXC //see if it matches libxc functionals - const std::pair> type_id = XC_Functional_Libxc::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; @@ -333,7 +333,7 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) std::string XC_Functional::output_info() { ModuleBase::TITLE("XC_Functional", "output_info"); - #ifdef USE_LIBXC +#ifdef USE_LIBXC if(use_libxc) { std::stringstream ss; @@ -349,7 +349,9 @@ std::string XC_Functional::output_info() { const func_reference_type *ref = xc_func_info_get_references(func.info, i); if(ref) + { ss<<"\t"< &stress_gga, const bool is_stress) +void XC_Functional::gradcorr( + double &etxc, + double &vtxc, + ModuleBase::matrix &v, + const Charge* const chr, + ModulePW::PW_Basis* rhopw, + const UnitCell *ucell, + std::vector &stress_gga, + const bool is_stress) { - ModuleBase::TITLE("XC_Functional","gradcorr"); - - 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; } - - 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; } - - assert(nspin0>0); - const double fac = 1.0/ nspin0; - - if(is_stress) - { - stress_gga.resize(9); - for(int i=0;i<9;i++) - { - stress_gga[i] = 0.0; - } - } - - // doing FFT to get rho in G space: rhog1 + ModuleBase::TITLE("XC_Functional","gradcorr"); + + if(func_type == 0 || func_type == 1) + { + return; + } + + bool igcc_is_lyp = false; + 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; + } + + assert(nspin0>0); + const double fac = 1.0/ nspin0; + + if(is_stress) + { + stress_gga.resize(9); + for(int i=0;i<9;i++) + { + stress_gga[i] = 0.0; + } + } + + // doing FFT to get rho in G space: rhog1 rhopw->real2recip(chr->rho[0], chr->rhog[0]); - if(PARAM.inp.nspin==2)//mohan fix bug 2012-05-28 - { - rhopw->real2recip(chr->rho[1], chr->rhog[1]); - } + if(PARAM.inp.nspin==2) + { + rhopw->real2recip(chr->rho[1], chr->rhog[1]); + } rhopw->real2recip(chr->rho_core, chr->rhog_core); - - // sum up (rho_core+rho) for each spin in real space - // and reciprocal space. - double* rhotmp1 = nullptr; - double* rhotmp2 = nullptr; - std::complex* rhogsum1 = nullptr; - std::complex* rhogsum2 = nullptr; - ModuleBase::Vector3* gdr1 = nullptr; - ModuleBase::Vector3* gdr2 = nullptr; - ModuleBase::Vector3* h1 = nullptr; - ModuleBase::Vector3* h2 = nullptr; - double* neg = nullptr; - double** vsave = nullptr; - double** vgg = nullptr; - - // for spin unpolarized case, - // calculate the gradient of (rho_core+rho) in reciprocal space. - rhotmp1 = new double[rhopw->nrxx]; - rhogsum1 = new std::complex[rhopw->npw]; + + // sum up (rho_core+rho) for each spin in real space + // and reciprocal space. + double* rhotmp1 = nullptr; + double* rhotmp2 = nullptr; + std::complex* rhogsum1 = nullptr; + std::complex* rhogsum2 = nullptr; + ModuleBase::Vector3* gdr1 = nullptr; + ModuleBase::Vector3* gdr2 = nullptr; + ModuleBase::Vector3* h1 = nullptr; + ModuleBase::Vector3* h2 = nullptr; + double* neg = nullptr; + double** vsave = nullptr; + double** vgg = nullptr; + + // for spin unpolarized case, + // calculate the gradient of (rho_core+rho) in reciprocal space. + rhotmp1 = new double[rhopw->nrxx]; + rhogsum1 = new std::complex[rhopw->npw]; #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) - { - rhotmp1[ir] = chr->rho[0][ir] + fac * chr->rho_core[ir]; - } + for(int ir=0; irnrxx; ir++) + { + rhotmp1[ir] = chr->rho[0][ir] + fac * chr->rho_core[ir]; + } #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ig=0; ignpw; ig++) - { - rhogsum1[ig] = chr->rhog[0][ig] + fac * chr->rhog_core[ig]; - } - - gdr1 = new ModuleBase::Vector3[rhopw->nrxx]; - if(!is_stress) { h1 = new ModuleBase::Vector3[rhopw->nrxx]; } - - XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); - - // for spin polarized case; - // calculate the gradient of (rho_core+rho) in reciprocal space. - if(PARAM.inp.nspin==2) - { - rhotmp2 = new double[rhopw->nrxx]; - rhogsum2 = new std::complex[rhopw->npw]; + for(int ig=0; ignpw; ig++) + { + rhogsum1[ig] = chr->rhog[0][ig] + fac * chr->rhog_core[ig]; + } + + gdr1 = new ModuleBase::Vector3[rhopw->nrxx]; + if(!is_stress) + { + h1 = new ModuleBase::Vector3[rhopw->nrxx]; + } + + XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); + + // for spin polarized case; + // calculate the gradient of (rho_core+rho) in reciprocal space. + if(PARAM.inp.nspin==2) + { + rhotmp2 = new double[rhopw->nrxx]; + rhogsum2 = new std::complex[rhopw->npw]; #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) - { - rhotmp2[ir] = chr->rho[1][ir] + fac * chr->rho_core[ir]; - } + for(int ir=0; irnrxx; ir++) + { + rhotmp2[ir] = chr->rho[1][ir] + fac * chr->rho_core[ir]; + } #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ig=0; ignpw; ig++) - { - rhogsum2[ig] = chr->rhog[1][ig] + fac * chr->rhog_core[ig]; - } - - gdr2 = new ModuleBase::Vector3[rhopw->nrxx]; - if(!is_stress) { h2 = new ModuleBase::Vector3[rhopw->nrxx]; } - - XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); - } - - if(PARAM.inp.nspin == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) - { - rhotmp2 = new double[rhopw->nrxx]; - rhogsum2 = new std::complex[rhopw->npw]; - neg = new double [rhopw->nrxx]; + for(int ig=0; ignpw; ig++) + { + rhogsum2[ig] = chr->rhog[1][ig] + fac * chr->rhog_core[ig]; + } + + gdr2 = new ModuleBase::Vector3[rhopw->nrxx]; + if(!is_stress) + { + h2 = new ModuleBase::Vector3[rhopw->nrxx]; + } + + XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); + } + + if(PARAM.inp.nspin == 4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) + { + rhotmp2 = new double[rhopw->nrxx]; + rhogsum2 = new std::complex[rhopw->npw]; + neg = new double [rhopw->nrxx]; #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) - { - rhotmp1[ir] = 0.0; - rhotmp2[ir] = 0.0; - neg[ir] = 0.0; - } + for(int ir=0; irnrxx; ir++) + { + rhotmp1[ir] = 0.0; + rhotmp2[ir] = 0.0; + neg[ir] = 0.0; + } #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ig=0; ignpw; ig++) - { - rhogsum1[ig] = 0.0; - rhogsum2[ig] = 0.0; - } - if(!is_stress) - { - vsave = new double* [PARAM.inp.nspin]; - for(int is = 0;isnrxx]; - } + for(int ig=0; ignpw; ig++) + { + rhogsum1[ig] = 0.0; + rhogsum2[ig] = 0.0; + } + if(!is_stress) + { + vsave = new double* [PARAM.inp.nspin]; + for(int is = 0;isnrxx]; + } #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 1024) #endif - for(int is = 0;isnrxx;ir++){ - vsave[is][ir] = v(is,ir); - v(is,ir) = 0; - } - } - vgg = new double* [nspin0]; - for(int is = 0;isnrxx]; -} - } - noncolin_rho(rhotmp1,rhotmp2,neg,chr->rho,rhopw->nrxx,ucell->magnet.ux_,ucell->magnet.lsign_); - rhopw->real2recip(rhotmp1, rhogsum1); - rhopw->real2recip(rhotmp2, rhogsum2); + for(int is = 0;isnrxx;ir++) + { + vsave[is][ir] = v(is,ir); + v(is,ir) = 0; + } + } + vgg = new double* [nspin0]; + for(int is = 0;isnrxx]; + } + } + noncolin_rho(rhotmp1, rhotmp2, neg, chr->rho, rhopw->nrxx, ucell->magnet.ux_, ucell->magnet.lsign_); + rhopw->real2recip(rhotmp1, rhogsum1); + rhopw->real2recip(rhotmp2, rhogsum2); #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) - { - rhotmp2[ir] += fac * chr->rho_core[ir]; - rhotmp1[ir] += fac * chr->rho_core[ir]; - } + for(int ir=0; irnrxx; ir++) + { + rhotmp2[ir] += fac * chr->rho_core[ir]; + rhotmp1[ir] += fac * chr->rho_core[ir]; + } #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ig=0; ignpw; ig++) - { - rhogsum2[ig] += fac * chr->rhog_core[ig]; - rhogsum1[ig] += fac * chr->rhog_core[ig]; - } - - gdr2 = new ModuleBase::Vector3[rhopw->nrxx]; - if(!is_stress) h2 = new ModuleBase::Vector3[rhopw->nrxx]; - - XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); - XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); + for(int ig=0; ignpw; ig++) + { + rhogsum2[ig] += fac * chr->rhog_core[ig]; + rhogsum1[ig] += fac * chr->rhog_core[ig]; + } + + gdr2 = new ModuleBase::Vector3[rhopw->nrxx]; + if(!is_stress) + { + h2 = new ModuleBase::Vector3[rhopw->nrxx]; + } + + XC_Functional::grad_rho( rhogsum1 , gdr1, rhopw, ucell->tpiba); + XC_Functional::grad_rho( rhogsum2 , gdr2, rhopw, ucell->tpiba); + } - } - - const double epsr = 1.0e-6; - const double epsg = 1.0e-10; + const double epsr = 1.0e-6; + const double epsg = 1.0e-10; - double vtxcgc = 0.0; - double etxcgc = 0.0; + double vtxcgc = 0.0; + double etxcgc = 0.0; #ifdef _OPENMP #pragma omp parallel -{ - std::vector local_stress_gga; - double local_vtxcgc = 0.0; - double local_etxcgc = 0.0; - - if(is_stress) - { - local_stress_gga.resize(9); - for(int i=0;i<9;i++) - { - local_stress_gga[i] = 0.0; - } - } + { + std::vector local_stress_gga; + double local_vtxcgc = 0.0; + double local_etxcgc = 0.0; + + if(is_stress) + { + local_stress_gga.resize(9); + for(int i=0;i<9;i++) + { + local_stress_gga[i] = 0.0; + } + } #else - std::vector &local_stress_gga = stress_gga; - double &local_vtxcgc = vtxcgc; - double &local_etxcgc = etxcgc; + std::vector &local_stress_gga = stress_gga; + double &local_vtxcgc = vtxcgc; + double &local_etxcgc = etxcgc; #endif - double grho2a = 0.0; - double grho2b = 0.0; - double sxc = 0.0; - double v1xc = 0.0; - double v2xc = 0.0; + double grho2a = 0.0; + double grho2b = 0.0; + double sxc = 0.0; + double v1xc = 0.0; + double v2xc = 0.0; - if(nspin0==1) - { - double segno = 0.0; + if(nspin0==1) + { + double segno = 0.0; #ifdef _OPENMP #pragma omp for #endif - for(int ir=0; irnrxx; ir++) - { - const double arho = std::abs( rhotmp1[ir] ); - if(!is_stress) { h1[ir].x = h1[ir].y = h1[ir].z = 0.0; -} - - if(arho > epsr) - { - grho2a = gdr1[ir].norm2(); - - //normally values in rhotmp can either be >= 0 or < 0. - if( rhotmp1[ir] >= 0.0 ) { segno = 1.0; - } else { segno = -1.0; -} - if (use_libxc && is_stress) - { + for(int ir=0; irnrxx; ir++) + { + const double arho = std::abs( rhotmp1[ir] ); + if(!is_stress) + { + h1[ir].x = 0.0; + h1[ir].y = 0.0; + h1[ir].z = 0.0; + } + + if(arho > epsr) + { + grho2a = gdr1[ir].norm2(); + + //normally values in rhotmp can either be >= 0 or < 0. + if( rhotmp1[ir] >= 0.0 ) + { + segno = 1.0; + } + else + { + segno = -1.0; + } + if (use_libxc && is_stress) + { #ifdef USE_LIBXC - if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA - { - double v3xc = 0.0; - double atau = chr->kin_r[0][ir]/2.0; - XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc); - } - else - { - XC_Functional_Libxc::gcxc_libxc( func_id, arho, grho2a, sxc, v1xc, v2xc); - } -#endif - } // end use_libxc - else - { - XC_Functional::gcxc( arho, grho2a, sxc, v1xc, v2xc); - } - if(is_stress) - { - double tt[3]; - tt[0] = gdr1[ir].x; - tt[1] = gdr1[ir].y; - tt[2] = gdr1[ir].z; - for(int l = 0;l< 3;l++) - { - for(int m = 0;m< l+1;m++) - { - int ind = l*3 + m; - local_stress_gga[ind] += tt[l] * tt[m] * ModuleBase::e2 * v2xc; - } - } - } - else - { - // first term of the gradient correction: - // D(rho*Exc)/D(rho) - v(0, ir) += ModuleBase::e2 * v1xc; - // cout << "v " << v(0, ir) << endl; - - // h contains - // D(rho*Exc) / D(|grad rho|) * (grad rho) / |grad rho| - h1[ir] = ModuleBase::e2 * v2xc * gdr1[ir]; - - local_vtxcgc += ModuleBase::e2* v1xc * ( rhotmp1[ir] - chr->rho_core[ir] ); - local_etxcgc += ModuleBase::e2* sxc * segno; - } - } // end arho > epsr - } - }// end nspin0 == 1 - else // spin polarized case - { + if(func_type == 3 || func_type == 5) + { + double v3xc = 0.0; + double atau = chr->kin_r[0][ir]/2.0; + XC_Functional_Libxc::tau_xc( func_id, arho, grho2a, atau, sxc, v1xc, v2xc, v3xc); + } + else + { + XC_Functional_Libxc::gcxc_libxc( func_id, arho, grho2a, sxc, v1xc, v2xc); + } +#endif + } + else + { + XC_Functional::gcxc( arho, grho2a, sxc, v1xc, v2xc); + } + if(is_stress) + { + double tt[3]; + tt[0] = gdr1[ir].x; + tt[1] = gdr1[ir].y; + tt[2] = gdr1[ir].z; + for(int l = 0;l< 3;l++) + { + for(int m = 0;m< l+1;m++) + { + int ind = l*3 + m; + local_stress_gga[ind] += tt[l] * tt[m] * ModuleBase::e2 * v2xc; + } + } + } + else + { + // first term of the gradient correction: + // D(rho*Exc)/D(rho) + v(0, ir) += ModuleBase::e2 * v1xc; + + // h contains + // D(rho*Exc) / D(|grad rho|) * (grad rho) / |grad rho| + h1[ir] = ModuleBase::e2 * v2xc * gdr1[ir]; + + local_vtxcgc += ModuleBase::e2* v1xc * ( rhotmp1[ir] - chr->rho_core[ir] ); + local_etxcgc += ModuleBase::e2* sxc * segno; + } + } + } + } + else + { #ifdef _OPENMP #pragma omp for #endif - for(int ir=0; irnrxx; ir++) - { - if(use_libxc) - { + for(int ir=0; irnrxx; ir++) + { + if(use_libxc) + { #ifdef USE_LIBXC - double sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud; - if(func_type == 3 || func_type == 5) //the gradcorr part to stress of mGGA - { - double v3xcup, v3xcdw; - double atau1 = chr->kin_r[0][ir]/2.0; - double atau2 = chr->kin_r[1][ir]/2.0; - 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_Libxc::gcxc_spin_libxc( - func_id, - rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], - sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud); - } - if(is_stress) - { - double tt1[3],tt2[3]; - { - tt1[0] = gdr1[ir].x; - tt1[1] = gdr1[ir].y; - tt1[2] = gdr1[ir].z; - tt2[0] = gdr2[ir].x; - tt2[1] = gdr2[ir].y; - tt2[2] = gdr2[ir].z; - } - for(int l = 0;l< 3;l++) - { - for(int m = 0;m< l+1;m++) - { - int ind = l*3 + m; - local_stress_gga [ind] += ( tt1[l] * tt1[m] * v2xcup + - tt2[l] * tt2[m] * v2xcdw + - (tt1[l] * tt2[m] + - tt2[l] * tt1[m] ) * v2xcud ) * ModuleBase::e2; - } - } - } - else - { - // first term of the gradient correction : D(rho*Exc)/D(rho) - v(0,ir) += ModuleBase::e2 * v1xcup; - v(1,ir) += ModuleBase::e2 * v1xcdw; - - // h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho| - h1[ir] += ModuleBase::e2 * ( v2xcup * gdr1[ir] + v2xcud * gdr2[ir] ); - h2[ir] += ModuleBase::e2 * ( v2xcdw * gdr2[ir] + v2xcud * gdr1[ir] ); - - local_vtxcgc = local_vtxcgc + ModuleBase::e2 * v1xcup * ( rhotmp1[ir] - chr->rho_core[ir] * fac ); - local_vtxcgc = local_vtxcgc + ModuleBase::e2 * v1xcdw * ( rhotmp2[ir] - chr->rho_core[ir] * fac ); - local_etxcgc = local_etxcgc + ModuleBase::e2 * sxc; - } -#endif - } - else - { - double v1cup = 0.0; - double v1cdw = 0.0; - double v2cup = 0.0; - double v2cdw = 0.0; - double v1xup = 0.0; - double v1xdw = 0.0; - double v2xup = 0.0; - double v2xdw = 0.0; - double v2cud = 0.0; - double v2c = 0.0; - double sx = 0.0; - double sc = 0.0; - double rh = rhotmp1[ir] + rhotmp2[ir]; - grho2a = gdr1[ir].norm2(); - grho2b = gdr2[ir].norm2(); - XC_Functional::gcx_spin(rhotmp1[ir], rhotmp2[ir], grho2a, grho2b, - sx, v1xup, v1xdw, v2xup, v2xdw); - - if(rh > epsr) - { - if(igcc_is_lyp) - { - ModuleBase::WARNING_QUIT("XC_Functional","igcc_is_lyp is not available now."); - } - else - { - double zeta = ( rhotmp1[ir] - rhotmp2[ir] ) / rh; - if(PARAM.inp.nspin==4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) { zeta = fabs(zeta) * neg[ir]; -} - const double grh2 = (gdr1[ir]+gdr2[ir]).norm2(); - XC_Functional::gcc_spin(rh, zeta, grh2, sc, v1cup, v1cdw, v2c); - v2cup = v2c; - v2cdw = v2c; - v2cud = v2c; - } - } - else - { - sc = 0.0; - v1cup = 0.0; - v1cdw = 0.0; - v2c = 0.0; - v2cup = 0.0; - v2cdw = 0.0; - v2cud = 0.0; - } - - if(is_stress) - { - double tt1[3],tt2[3]; - { - tt1[0] = gdr1[ir].x; - tt1[1] = gdr1[ir].y; - tt1[2] = gdr1[ir].z; - tt2[0] = gdr2[ir].x; - tt2[1] = gdr2[ir].y; - tt2[2] = gdr2[ir].z; - } - for(int l = 0;l< 3;l++) - { - for(int m = 0;m< l+1;m++) - { - int ind = l*3 + m; - // exchange - local_stress_gga [ind] += tt1[l] * tt1[m] * ModuleBase::e2 * v2xup + - tt2[l] * tt2[m] * ModuleBase::e2 * v2xdw; - // correlation - local_stress_gga [ind] += ( tt1[l] * tt1[m] * v2cup + - tt2[l] * tt2[m] * v2cdw + - (tt1[l] * tt2[m] + - tt2[l] * tt1[m] ) * v2cud ) * ModuleBase::e2; - } - } - } - else - { - // first term of the gradient correction : D(rho*Exc)/D(rho) - v(0,ir) = v(0,ir) + ModuleBase::e2 * ( v1xup + v1cup ); - v(1,ir) = v(1,ir) + ModuleBase::e2 * ( v1xdw + v1cdw ); - - // h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho| - h1[ir] = ModuleBase::e2 * ( ( v2xup + v2cup ) * gdr1[ir] + v2cud * gdr2[ir] ); - h2[ir] = ModuleBase::e2 * ( ( v2xdw + v2cdw ) * gdr2[ir] + v2cud * gdr1[ir] ); - - local_vtxcgc = local_vtxcgc + ModuleBase::e2 * ( v1xup + v1cup ) * ( rhotmp1[ir] - chr->rho_core[ir] * fac ); - local_vtxcgc = local_vtxcgc + ModuleBase::e2 * ( v1xdw + v1cdw ) * ( rhotmp2[ir] - chr->rho_core[ir] * fac ); - local_etxcgc = local_etxcgc + ModuleBase::e2 * ( sx + sc ); - } - } - }// end ir - - } + double sxc = 0.0; + double v1xcup = 0.0; + double v1xcdw = 0.0; + double v2xcup = 0.0; + double v2xcdw = 0.0; + double v2xcud = 0.0; + if(func_type == 3 || func_type == 5) + { + double v3xcup = 0.0; + double v3xcdw = 0.0; + double atau1 = chr->kin_r[0][ir]/2.0; + double atau2 = chr->kin_r[1][ir]/2.0; + 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_Libxc::gcxc_spin_libxc( + func_id, + rhotmp1[ir], rhotmp2[ir], gdr1[ir], gdr2[ir], + sxc, v1xcup, v1xcdw, v2xcup, v2xcdw, v2xcud); + } + if(is_stress) + { + double tt1[3],tt2[3]; + { + tt1[0] = gdr1[ir].x; + tt1[1] = gdr1[ir].y; + tt1[2] = gdr1[ir].z; + tt2[0] = gdr2[ir].x; + tt2[1] = gdr2[ir].y; + tt2[2] = gdr2[ir].z; + } + for(int l = 0;l< 3;l++) + { + for(int m = 0;m< l+1;m++) + { + int ind = l*3 + m; + local_stress_gga [ind] += ( tt1[l] * tt1[m] * v2xcup + + tt2[l] * tt2[m] * v2xcdw + + (tt1[l] * tt2[m] + tt2[l] * tt1[m] ) * v2xcud ) * ModuleBase::e2; + } + } + } + else + { + // first term of the gradient correction : D(rho*Exc)/D(rho) + v(0,ir) += ModuleBase::e2 * v1xcup; + v(1,ir) += ModuleBase::e2 * v1xcdw; + + // h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho| + h1[ir] += ModuleBase::e2 * ( v2xcup * gdr1[ir] + v2xcud * gdr2[ir] ); + h2[ir] += ModuleBase::e2 * ( v2xcdw * gdr2[ir] + v2xcud * gdr1[ir] ); + + local_vtxcgc = local_vtxcgc + ModuleBase::e2 * v1xcup * ( rhotmp1[ir] - chr->rho_core[ir] * fac ); + local_vtxcgc = local_vtxcgc + ModuleBase::e2 * v1xcdw * ( rhotmp2[ir] - chr->rho_core[ir] * fac ); + local_etxcgc = local_etxcgc + ModuleBase::e2 * sxc; + } +#endif + } + else + { + double v1cup = 0.0; + double v1cdw = 0.0; + double v2cup = 0.0; + double v2cdw = 0.0; + double v1xup = 0.0; + double v1xdw = 0.0; + double v2xup = 0.0; + double v2xdw = 0.0; + double v2cud = 0.0; + double v2c = 0.0; + double sx = 0.0; + double sc = 0.0; + double rh = rhotmp1[ir] + rhotmp2[ir]; + grho2a = gdr1[ir].norm2(); + grho2b = gdr2[ir].norm2(); + XC_Functional::gcx_spin(rhotmp1[ir], rhotmp2[ir], grho2a, grho2b, + sx, v1xup, v1xdw, v2xup, v2xdw); + + if(rh > epsr) + { + if(igcc_is_lyp) + { + ModuleBase::WARNING_QUIT("XC_Functional","igcc_is_lyp is not available now."); + } + else + { + double zeta = ( rhotmp1[ir] - rhotmp2[ir] ) / rh; + if(PARAM.inp.nspin==4&&(PARAM.globalv.domag||PARAM.globalv.domag_z)) + { + zeta = fabs(zeta) * neg[ir]; + } + const double grh2 = (gdr1[ir]+gdr2[ir]).norm2(); + XC_Functional::gcc_spin(rh, zeta, grh2, sc, v1cup, v1cdw, v2c); + v2cup = v2c; + v2cdw = v2c; + v2cud = v2c; + } + } + else + { + sc = 0.0; + v1cup = 0.0; + v1cdw = 0.0; + v2c = 0.0; + v2cup = 0.0; + v2cdw = 0.0; + v2cud = 0.0; + } + + if(is_stress) + { + double tt1[3],tt2[3]; + { + tt1[0] = gdr1[ir].x; + tt1[1] = gdr1[ir].y; + tt1[2] = gdr1[ir].z; + tt2[0] = gdr2[ir].x; + tt2[1] = gdr2[ir].y; + tt2[2] = gdr2[ir].z; + } + for(int l = 0;l< 3;l++) + { + for(int m = 0;m< l+1;m++) + { + int ind = l*3 + m; + // exchange + local_stress_gga [ind] += tt1[l] * tt1[m] * ModuleBase::e2 * v2xup + + tt2[l] * tt2[m] * ModuleBase::e2 * v2xdw; + // correlation + local_stress_gga [ind] += ( tt1[l] * tt1[m] * v2cup + + tt2[l] * tt2[m] * v2cdw + + (tt1[l] * tt2[m] + tt2[l] * tt1[m] ) * v2cud ) * ModuleBase::e2; + } + } + } + else + { + // first term of the gradient correction : D(rho*Exc)/D(rho) + v(0,ir) = v(0,ir) + ModuleBase::e2 * ( v1xup + v1cup ); + v(1,ir) = v(1,ir) + ModuleBase::e2 * ( v1xdw + v1cdw ); + + // h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho| + h1[ir] = ModuleBase::e2 * ( ( v2xup + v2cup ) * gdr1[ir] + v2cud * gdr2[ir] ); + h2[ir] = ModuleBase::e2 * ( ( v2xdw + v2cdw ) * gdr2[ir] + v2cud * gdr1[ir] ); + + local_vtxcgc = local_vtxcgc + ModuleBase::e2 * ( v1xup + v1cup ) * ( rhotmp1[ir] - chr->rho_core[ir] * fac ); + local_vtxcgc = local_vtxcgc + ModuleBase::e2 * ( v1xdw + v1cdw ) * ( rhotmp2[ir] - chr->rho_core[ir] * fac ); + local_etxcgc = local_etxcgc + ModuleBase::e2 * ( sx + sc ); + } + } + } + } #ifdef _OPENMP - #pragma omp critical(xc_functional_gradcorr_reduce) - { - if(is_stress) - { - for(int l = 0;l< 3;l++) - { - for(int m = 0;m< l+1;m++) - { - int ind = l*3 + m; - stress_gga [ind] += local_stress_gga [ind]; - } - } - } - else - { - vtxcgc += local_vtxcgc; - etxcgc += local_etxcgc; - } - } + #pragma omp critical(xc_functional_gradcorr_reduce) + { + if(is_stress) + { + for(int l = 0;l< 3;l++) + { + for(int m = 0;m< l+1;m++) + { + int ind = l*3 + m; + stress_gga [ind] += local_stress_gga [ind]; + } + } + } + else + { + vtxcgc += local_vtxcgc; + etxcgc += local_etxcgc; + } + } } #endif - //std::cout << "\n vtxcgc=" << vtxcgc; - //std::cout << "\n etxcgc=" << etxcgc << std::endl; - - if(!is_stress) - { + if(!is_stress) + { #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) - { - rhotmp1[ir] -= fac * chr->rho_core[ir]; - } - if(nspin0==2) - { + for(int ir=0; irnrxx; ir++) + { + rhotmp1[ir] -= fac * chr->rho_core[ir]; + } + if(nspin0==2) + { #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) - { - rhotmp2[ir] -= fac * chr->rho_core[ir]; - } - } - - // second term of the gradient correction : - // \sum_alpha (D / D r_alpha) ( D(rho*Exc)/D(grad_alpha rho) ) - - // dh is in real sapce. - double* dh = new double[rhopw->nrxx]; - - for(int is=0; istpiba); -} - if(is==1) {XC_Functional::grad_dot(h2,dh,rhopw,ucell->tpiba); -} + for(int ir=0; irnrxx; ir++) + { + rhotmp2[ir] -= fac * chr->rho_core[ir]; + } + } + + // second term of the gradient correction : + // sum_alpha (D / D r_alpha) ( D(rho*Exc)/D(grad_alpha rho) ) + + // dh is in real sapce. + double* dh = new double[rhopw->nrxx]; + + for(int is=0; istpiba); + } + if(is==1) + { + XC_Functional::grad_dot(h2,dh,rhopw,ucell->tpiba); + } #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) { - v(is, ir) -= dh[ir]; -} - - double sum = 0.0; - if(is==0) - { + for(int ir=0; irnrxx; ir++) + { + v(is, ir) -= dh[ir]; + } + + double sum = 0.0; + if(is==0) + { #ifdef _OPENMP #pragma omp parallel for reduction(+:sum) schedule(static, 256) #endif - for(int ir=0; irnrxx; ir++) { - sum += dh[ir] * rhotmp1[ir]; -} - } - else if(is==1) - { + for(int ir=0; irnrxx; ir++) + { + sum += dh[ir] * rhotmp1[ir]; + } + } + else if(is==1) + { #ifdef _OPENMP #pragma omp parallel for reduction(+:sum) schedule(static, 256) #endif - for(int ir=0; irnrxx; ir++) { - sum += dh[ir] * rhotmp2[ir]; -} - } - vtxcgc -= sum; - } - - delete[] dh; + for(int ir=0; irnrxx; ir++) + { + sum += dh[ir] * rhotmp2[ir]; + } + } + vtxcgc -= sum; + } + + delete[] dh; - vtxc += vtxcgc; - etxc += etxcgc; + vtxc += vtxcgc; + etxc += etxcgc; - if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z)) - { + if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z)) + { #ifdef _OPENMP #pragma omp parallel for collapse(2) schedule(static, 1024) #endif - for(int is=0;isnrxx;ir++) - { - if(isnrxx;ir++) + { + if(isnrxx;ir++) - { - v(0,ir) += 0.5 * (vgg[0][ir] + vgg[1][ir]); - double amag = sqrt(pow(chr->rho[1][ir],2)+pow(chr->rho[2][ir],2)+pow(chr->rho[3][ir],2)); - if(amag>1e-12) - { - for(int i=1;i<4;i++) { - v(i,ir)+= neg[ir] * 0.5 *(vgg[0][ir]-vgg[1][ir])*chr->rho[i][ir]/amag; -} - } - } - } - } - // deacllocate - delete[] rhotmp1; - delete[] rhogsum1; - delete[] gdr1; - if(!is_stress) { delete[] h1; -} + for(int ir=0;irnrxx;ir++) + { + v(0,ir) += 0.5 * (vgg[0][ir] + vgg[1][ir]); + double amag = sqrt(pow(chr->rho[1][ir],2)+pow(chr->rho[2][ir],2)+pow(chr->rho[3][ir],2)); + if(amag>1e-12) + { + for(int i=1;i<4;i++) + { + v(i,ir)+= neg[ir] * 0.5 *(vgg[0][ir]-vgg[1][ir])*chr->rho[i][ir]/amag; + } + } + } + } + } + // deacllocate + delete[] rhotmp1; + delete[] rhogsum1; + delete[] gdr1; + if(!is_stress) + { + delete[] h1; + } - if(PARAM.inp.nspin==2) - { - delete[] rhotmp2; - delete[] rhogsum2; - delete[] gdr2; - if(!is_stress) { delete[] h2; -} - } - if(PARAM.inp.nspin == 4 && (PARAM.globalv.domag||PARAM.globalv.domag_z)) - { - delete[] neg; - if(!is_stress) - { - for(int i=0; i @@ -619,163 +678,186 @@ void XC_Functional::grad_wfc( const int ik, const Real tpiba, const ModulePW::PW_Basis_K* wfc_basis, - const T* rhog, + const T* rhog, T* grad) { using ct_Device = typename ct::PsiToContainer::type; - const int npw_k = wfc_basis->npwk[ik]; - - auto porter = std::move(ct::Tensor( + const int npw_k = wfc_basis->npwk[ik]; + + auto porter = std::move(ct::Tensor( ct::DataTypeToEnum::value, ct::DeviceTypeToEnum::value, {wfc_basis->nmaxgr})); - auto gcar = ct::TensorMap( - &wfc_basis->gcar[0][0], ct::DataType::DT_DOUBLE, ct::DeviceType::CpuDevice, {wfc_basis->nks * wfc_basis->npwk_max, 3}).to_device(); - auto kvec_c = ct::TensorMap( - &wfc_basis->kvec_c[0][0],ct::DataType::DT_DOUBLE, ct::DeviceType::CpuDevice, {wfc_basis->nks, 3}).to_device(); - - auto xc_functional_grad_wfc_solver - = hamilt::xc_functional_grad_wfc_op(); - - for(int ipol=0; ipol<3; ipol++) { - xc_functional_grad_wfc_solver( - ik, ipol, npw_k, wfc_basis->npwk_max, // Integers - tpiba, // Double - gcar.template data(), // Array of Real - kvec_c.template data(), // Array of double - rhog, porter.data()); // Array of std::complex - - // bring the gdr from G --> R - Device * ctx = nullptr; - wfc_basis->recip_to_real(ctx, porter.data(), porter.data(), ik); - - xc_functional_grad_wfc_solver( - ipol, wfc_basis->nrxx, // Integers - porter.data(), grad); // Array of std::complex + auto gcar = ct::TensorMap( + &wfc_basis->gcar[0][0], ct::DataType::DT_DOUBLE, ct::DeviceType::CpuDevice, {wfc_basis->nks * wfc_basis->npwk_max, 3}).to_device(); + auto kvec_c = ct::TensorMap( + &wfc_basis->kvec_c[0][0],ct::DataType::DT_DOUBLE, ct::DeviceType::CpuDevice, {wfc_basis->nks, 3}).to_device(); + + auto xc_functional_grad_wfc_solver + = hamilt::xc_functional_grad_wfc_op(); + + for(int ipol=0; ipol<3; ipol++) + { + xc_functional_grad_wfc_solver( + ik, ipol, npw_k, wfc_basis->npwk_max, + tpiba, + gcar.template data(), + kvec_c.template data(), + rhog, porter.data()); + + // bring the gdr from G --> R + Device * ctx = nullptr; + wfc_basis->recip_to_real(ctx, porter.data(), porter.data(), ik); + + xc_functional_grad_wfc_solver( + ipol, wfc_basis->nrxx, + porter.data(), grad); } } -void XC_Functional::grad_rho(const std::complex* rhog, - ModuleBase::Vector3* gdr, - const ModulePW::PW_Basis* rho_basis, - const double tpiba) +void XC_Functional::grad_rho( + const std::complex* rhog, + ModuleBase::Vector3* gdr, + const ModulePW::PW_Basis* rho_basis, + const double tpiba) { - std::complex *gdrtmp = new std::complex[rho_basis->nmaxgr]; + std::complex *gdrtmp = new std::complex[rho_basis->nmaxgr]; - // the formula is : rho(r)^prime = \int iG * rho(G)e^{iGr} dG - for(int i = 0 ; i < 3 ; ++i) - { - // calculate the charge density gradient in reciprocal space. + // the formula is : rho(r)^prime = int iG * rho(G)e^{iGr} dG + for(int i = 0 ; i < 3 ; ++i) + { + // calculate the charge density gradient in reciprocal space. #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ig=0; ignpw; ig++) { - gdrtmp[ig] = ModuleBase::IMAG_UNIT * rhog[ig] * rho_basis->gcar[ig][i]; -} + for(int ig=0; ignpw; ig++) + { + gdrtmp[ig] = ModuleBase::IMAG_UNIT * rhog[ig] * rho_basis->gcar[ig][i]; + } - // bring the gdr from G --> R - rho_basis->recip2real(gdrtmp, gdrtmp); + // bring the gdr from G --> R + rho_basis->recip2real(gdrtmp, gdrtmp); - // remember to multily 2pi/a0, which belongs to G vectors. + // remember to multily 2pi/a0, which belongs to G vectors. #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) { - gdr[ir][i] = gdrtmp[ir].real() * tpiba; -} - } + for(int ir=0; irnrxx; ir++) + { + gdr[ir][i] = gdrtmp[ir].real() * tpiba; + } + } - delete[] gdrtmp; - return; + delete[] gdrtmp; + return; } -void XC_Functional::grad_dot(const ModuleBase::Vector3* h, double* dh, const ModulePW::PW_Basis* rho_basis, const double tpiba) +void XC_Functional::grad_dot( + const ModuleBase::Vector3* h, + double* dh, + const ModulePW::PW_Basis* rho_basis, + const double tpiba) { - std::complex *aux = new std::complex[rho_basis->nmaxgr]; - std::complex *gaux = new std::complex[rho_basis->npw]; + std::complex *aux = new std::complex[rho_basis->nmaxgr]; + std::complex *gaux = new std::complex[rho_basis->npw]; - for(int i = 0 ; i < 3 ; ++i) - { + for(int i = 0 ; i < 3 ; ++i) + { #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir = 0; ir < rho_basis->nrxx; ++ir) { - aux[ir] = std::complex( h[ir][i], 0.0); -} - - // bring to G space. - rho_basis->real2recip(aux,aux); - if (i == 0) - { + for(int ir = 0; ir < rho_basis->nrxx; ++ir) + { + aux[ir] = std::complex( h[ir][i], 0.0); + } + + // bring to G space. + rho_basis->real2recip(aux,aux); + if (i == 0) + { #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ig = 0; ig < rho_basis->npw; ++ig) { - gaux[ig] = ModuleBase::IMAG_UNIT * aux[ig] * rho_basis->gcar[ig][i]; -} - } - else - { + for(int ig = 0; ig < rho_basis->npw; ++ig) + { + gaux[ig] = ModuleBase::IMAG_UNIT * aux[ig] * rho_basis->gcar[ig][i]; + } + } + else + { #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ig = 0; ig < rho_basis->npw; ++ig) { - gaux[ig] += ModuleBase::IMAG_UNIT * aux[ig] * rho_basis->gcar[ig][i]; -} - } - } + for(int ig = 0; ig < rho_basis->npw; ++ig) + { + gaux[ig] += ModuleBase::IMAG_UNIT * aux[ig] * rho_basis->gcar[ig][i]; + } + } + } - // bring back to R space - rho_basis->recip2real(gaux,aux); + // bring back to R space + rho_basis->recip2real(gaux,aux); #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir=0; irnrxx; ir++) { - dh[ir] = aux[ir].real() * tpiba; -} - - delete[] aux; - delete[] gaux; - return; + for(int ir=0; irnrxx; ir++) + { + dh[ir] = aux[ir].real() * tpiba; + } + + delete[] aux; + delete[] gaux; + return; } -void XC_Functional::noncolin_rho(double *rhoout1, double *rhoout2, double *neg, - const double*const*const rho, const int nrxx, const double* ux_, const bool lsign_) +void XC_Functional::noncolin_rho( + double *rhoout1, + double *rhoout2, + double *neg, + const double*const*const rho, + const int nrxx, + const double* ux_, + const bool lsign_) { - //this function diagonalizes the spin density matrix and gives as output the - //spin up and spin down components of the charge. - //If lsign is true up and dw are with respect to the fixed quantization axis - //ux, otherwise rho + |m| is always rhoup and rho-|m| is always rhodw. + //this function diagonalizes the spin density matrix and gives as output the + //spin up and spin down components of the charge. + //If lsign is true up and dw are with respect to the fixed quantization axis + //ux, otherwise rho + |m| is always rhoup and rho-|m| is always rhodw. #ifdef _OPENMP #pragma omp parallel for schedule(static, 1024) #endif - for(int ir = 0;ir0) { neg[ir] = 1.0; - } else { neg[ir] = -1.0; -} - } - } + for(int ir = 0;ir0) + { + neg[ir] = 1.0; + } + else + { + neg[ir] = -1.0; + } + } + } #ifdef _OPENMP #pragma omp parallel for #endif - for(int ir = 0;ir, base_device::DEVICE_CPU, double>( @@ -791,4 +873,4 @@ template void XC_Functional::grad_wfc, base_device::DEVICE_ const ModulePW::PW_Basis_K* wfc_basis, const std::complex* rhog, std::complex* grad); -#endif // __CUDA || __ROCM \ No newline at end of file +#endif diff --git a/source/source_hamilt/module_xc/xc_functional_vxc.cpp b/source/source_hamilt/module_xc/xc_functional_vxc.cpp index 9765f9db504..917b915f7b6 100644 --- a/source/source_hamilt/module_xc/xc_functional_vxc.cpp +++ b/source/source_hamilt/module_xc/xc_functional_vxc.cpp @@ -14,9 +14,10 @@ #endif // [etxc, vtxc, v] = XC_Functional::v_xc(...) -std::tuple XC_Functional::v_xc(const int& nrxx, // number of real-space grid - const Charge* const chr, - const UnitCell* ucell) // core charge density +std::tuple XC_Functional::v_xc( + const int& nrxx, + const Charge* const chr, + const UnitCell* ucell) { ModuleBase::TITLE("XC_Functional", "v_xc"); @@ -44,7 +45,6 @@ std::tuple XC_Functional::v_xc(const int& nr // the square of the e charge // in Rydeberg unit, so * 2.0. double e2 = 2.0; - double vanishing_charge = 1.0e-10; if (PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) @@ -57,7 +57,6 @@ std::tuple XC_Functional::v_xc(const int& nr { // total electron charge density double rhox = chr->rho[0][ir] + chr->rho_core[ir]; - double arhox = std::abs(rhox); if (arhox > vanishing_charge) { @@ -80,12 +79,12 @@ std::tuple XC_Functional::v_xc(const int& nr #endif for (int ir = 0;ir < nrxx;ir++) { - double rhox = chr->rho[0][ir] + chr->rho[1][ir] + chr->rho_core[ir]; //HLX(05-29-06): bug fixed + double rhox = chr->rho[0][ir] + chr->rho[1][ir] + chr->rho_core[ir]; double arhox = std::abs(rhox); if (arhox > vanishing_charge) { - double zeta = (chr->rho[0][ir] - chr->rho[1][ir]) / arhox; //HLX(05-29-06): bug fixed + double zeta = (chr->rho[0][ir] - chr->rho[1][ir]) / arhox; if (std::abs(zeta) > 1.0) { @@ -108,7 +107,7 @@ std::tuple XC_Functional::v_xc(const int& nr } } } - else if(PARAM.inp.nspin == 4)//noncollinear case added by zhengdy + else if(PARAM.inp.nspin == 4) { #ifdef _OPENMP #pragma omp parallel for reduction(+:etxc) reduction(+:vtxc) @@ -116,9 +115,7 @@ std::tuple XC_Functional::v_xc(const int& nr for(int ir = 0;irrho[1][ir],2) + pow(chr->rho[2][ir],2) + pow(chr->rho[3][ir],2) ); - double rhox = chr->rho[0][ir] + chr->rho_core[ir]; - double arhox = std::abs( rhox ); if ( arhox > vanishing_charge ) @@ -130,7 +127,7 @@ std::tuple XC_Functional::v_xc(const int& nr if ( std::abs( zeta ) > 1.0 ) { zeta = (zeta > 0.0) ? 1.0 : (-1.0); - }//end if + } if(use_libxc) { @@ -140,7 +137,7 @@ std::tuple XC_Functional::v_xc(const int& nr 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 +#endif } else { @@ -161,11 +158,11 @@ std::tuple XC_Functional::v_xc(const int& nr { v(ipol, ir) = e2 * vs * chr->rho[ipol][ir] / amag; vtxc += v(ipol,ir) * chr->rho[ipol][ir]; - }//end do - }//end if - }//end if - }//end do - }//end if + } + } + } + } + } // energy terms, local-density contributions // add gradient corrections (if any) diff --git a/source/source_hamilt/module_xc/xc_functional_wrapper_gcxc.cpp b/source/source_hamilt/module_xc/xc_functional_wrapper_gcxc.cpp index 3956d297025..cfaacd812a0 100644 --- a/source/source_hamilt/module_xc/xc_functional_wrapper_gcxc.cpp +++ b/source/source_hamilt/module_xc/xc_functional_wrapper_gcxc.cpp @@ -14,8 +14,12 @@ #include "xc_functional_libxc.h" #endif -void XC_Functional::gcxc(const double &rho, const double &grho, double &sxc, - double &v1xc, double &v2xc) +void XC_Functional::gcxc( + const double &rho, + const double &grho, + double &sxc, + double &v1xc, + double &v2xc) { //----------------------------------------------------------------------- // gradient corrections for exchange and correlation - Hartree a.u. @@ -41,8 +45,12 @@ void XC_Functional::gcxc(const double &rho, const double &grho, double &sxc, // real rho, grho, sx, sc, v1x, v2x, v1c, v2c; const double small = 1.e-6; const double smallg = 1.e-10; - double s,v1,v2; - sxc = v1xc = v2xc = 0.0; + double s = 0.0; + double v1 = 0.0; + double v2 = 0.0; + sxc = 0.0; + v1xc = 0.0; + v2xc = 0.0; if (rho <= small || grho < smallg) { @@ -53,47 +61,116 @@ void XC_Functional::gcxc(const double &rho, const double &grho, double &sxc, { switch( id ) { - case XC_GGA_X_B88: //B88 - XC_Functional::becke88(rho, grho, s, v1, v2);break; - case XC_GGA_X_PW91: //PW91_X - XC_Functional::ggax(rho, grho, s, v1, v2);break; - case XC_GGA_X_PBE: //PBX - XC_Functional::pbex(rho, grho, 0, s, v1, v2);break; - case XC_GGA_X_PBE_R: //revised PBX - XC_Functional::pbex(rho, grho, 1, s, v1, v2);break; - case XC_GGA_X_HCTH_A: //HCTH_X - XC_Functional::hcth(rho, grho, s, v1, v2);break; //XC together - case XC_GGA_C_HCTH_A: //HCTH_C - s = 0.0; v1 = 0.0; v2 = 0.0;break; - case XC_GGA_X_OPTX: //OPTX - XC_Functional::optx(rho, grho, s, v1, v2);break; - case XC_GGA_X_PBE_SOL: //PBXsol - XC_Functional::pbex(rho, grho, 2, s, v1, v2);break; - case XC_GGA_X_WC: //Wu-Cohen - XC_Functional::wcx (rho, grho, s, v1, v2);break; - case XC_GGA_C_P86: //P86 - XC_Functional::perdew86(rho, grho, s, v1, v2);break; - case XC_GGA_C_PW91: //PW91_C - XC_Functional::ggac(rho, grho, s, v1, v2);break; - case XC_GGA_C_PBE: //PBC - XC_Functional::pbec(rho, grho, 0, s, v1, v2);break; - case XC_GGA_C_PBE_SOL: //PBCsol - XC_Functional::pbec(rho, grho, 1, s, v1, v2);break; - case XC_GGA_C_LYP: //BLYP - XC_Functional::glyp(rho, grho, s, v1, v2); break; - case XC_HYB_GGA_XC_PBEH: //PBE0 - double sx, v1x, v2x, sc, v1c, v2c; + case XC_GGA_X_B88: + { + //B88 + XC_Functional::becke88(rho, grho, s, v1, v2); + break; + } + case XC_GGA_X_PW91: + { + //PW91_X + XC_Functional::ggax(rho, grho, s, v1, v2); + break; + } + case XC_GGA_X_PBE: + { + //PBX + XC_Functional::pbex(rho, grho, 0, s, v1, v2); + break; + } + case XC_GGA_X_PBE_R: + { + //revised PBX + XC_Functional::pbex(rho, grho, 1, s, v1, v2); + break; + } + case XC_GGA_X_HCTH_A: + { + //HCTH_X + XC_Functional::hcth(rho, grho, s, v1, v2); + break; + } + case XC_GGA_C_HCTH_A: + { + //HCTH_C + s = 0.0; + v1 = 0.0; + v2 = 0.0; + break; + } + case XC_GGA_X_OPTX: + { + //OPTX + XC_Functional::optx(rho, grho, s, v1, v2); + break; + } + case XC_GGA_X_PBE_SOL: + { + //PBXsol + XC_Functional::pbex(rho, grho, 2, s, v1, v2); + break; + } + case XC_GGA_X_WC: + { + //Wu-Cohen + XC_Functional::wcx (rho, grho, s, v1, v2); + break; + } + case XC_GGA_C_P86: + { + //P86 + XC_Functional::perdew86(rho, grho, s, v1, v2); + break; + } + case XC_GGA_C_PW91: + { + //PW91_C + XC_Functional::ggac(rho, grho, s, v1, v2); + break; + } + case XC_GGA_C_PBE: + { + //PBC + XC_Functional::pbec(rho, grho, 0, s, v1, v2); + break; + } + case XC_GGA_C_PBE_SOL: + { + //PBCsol + XC_Functional::pbec(rho, grho, 1, s, v1, v2); + break; + } + case XC_GGA_C_LYP: + { + //BLYP + XC_Functional::glyp(rho, grho, s, v1, v2); + break; + } + case XC_HYB_GGA_XC_PBEH: + { + //PBE0 + double sx = 0.0; + double v1x = 0.0; + double v2x = 0.0; + double sc = 0.0; + double v1c = 0.0; + double v2c = 0.0; XC_Functional::pbex(rho, grho, 0, sx, v1x, v2x); - sx *= (1.0 - XC_Functional::hybrid_alpha); - v1x *= (1.0 - XC_Functional::hybrid_alpha); + sx *= (1.0 - XC_Functional::hybrid_alpha); + v1x *= (1.0 - XC_Functional::hybrid_alpha); v2x *= (1.0 - XC_Functional::hybrid_alpha); XC_Functional::pbec(rho, grho, 0, sc, v1c, v2c); s = sx + sc; v1 = v1x + v1c; v2 = v2x + v2c; break; - default: //SCAN_X,SCAN_C,HSE, and so on + } + default: + { + //SCAN_X,SCAN_C,HSE, and so on throw std::domain_error("functional unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + } } sxc += s; v1xc += v1; @@ -104,8 +181,16 @@ void XC_Functional::gcxc(const double &rho, const double &grho, double &sxc, } //----------------------------------------------------------------------- -void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double grhodw2, - double &sx, double &v1xup, double &v1xdw, double &v2xup, double &v2xdw) +void XC_Functional::gcx_spin( + double rhoup, + double rhodw, + double grhoup2, + double grhodw2, + double &sx, + double &v1xup, + double &v1xdw, + double &v2xup, + double &v2xdw) { //-------------------------------------------------------------------- // gradient corrections for exchange - Hartree a.u. @@ -127,16 +212,20 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double // parameter : double small = 1.e-10; - double sxup, sxdw; + double sxup = 0.0; + double sxdw = 0.0; int iflag = 0; // exchange double rho = rhoup + rhodw; sx = 0.00; - v1xup = 0.00; v2xup = 0.00; - v1xdw = 0.00; v2xdw = 0.00; - sxup = 0.00; sxdw = 0.00; + v1xup = 0.00; + v2xup = 0.00; + v1xdw = 0.00; + v2xdw = 0.00; + sxup = 0.00; + sxdw = 0.00; if (rho <= small) { @@ -151,7 +240,9 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double int id = func_id[0]; switch( id ) { - case XC_GGA_X_B88: //B88 + case XC_GGA_X_B88: + { + //B88 if (rhoup > small && sqrt(fabs(grhoup2)) > small) { XC_Functional::becke88_spin(rhoup, grhoup2, sxup, v1xup, v2xup); @@ -161,7 +252,10 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double XC_Functional::becke88_spin(rhodw, grhodw2, sxdw, v1xdw, v2xdw); } break; - case XC_GGA_X_PBE: //PBX + } + case XC_GGA_X_PBE: + { + //PBX if (rhoup > small && sqrt(fabs(grhoup2)) > small) { XC_Functional::pbex(2.0 * rhoup, 4.0 * grhoup2, 0, sxup, v1xup, v2xup); @@ -171,7 +265,10 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double XC_Functional::pbex(2.0 * rhodw, 4.0 * grhodw2, 0, sxdw, v1xdw, v2xdw); } break; - case XC_GGA_X_PBE_R: //revised PBX + } + case XC_GGA_X_PBE_R: + { + //revised PBX if (rhoup > small && sqrt(fabs(grhoup2)) > small) { XC_Functional::pbex(2.0 * rhoup, 4.0 * grhoup2, 1, sxup, v1xup, v2xup); @@ -181,23 +278,29 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double XC_Functional::pbex(2.0 * rhodw, 4.0 * grhodw2, 1, sxdw, v1xdw, v2xdw); } break; - case XC_HYB_GGA_XC_PBEH: //PBE0 + } + case XC_HYB_GGA_XC_PBEH: + { + //PBE0 if (rhoup > small && sqrt(fabs(grhoup2)) > small) { XC_Functional::pbex(2.0 * rhoup, 4.0 * grhoup2, 0, sxup, v1xup, v2xup); - sxup *= (1.0 - XC_Functional::hybrid_alpha); - v1xup *= (1.0 - XC_Functional::hybrid_alpha); + sxup *= (1.0 - XC_Functional::hybrid_alpha); + v1xup *= (1.0 - XC_Functional::hybrid_alpha); v2xup *= (1.0 - XC_Functional::hybrid_alpha); } if (rhodw > small && sqrt(fabs(grhodw2)) > small) { XC_Functional::pbex(2.0 * rhodw, 4.0 * grhodw2, 0, sxdw, v1xdw, v2xdw); - sxdw *= (1.0 - XC_Functional::hybrid_alpha); - v1xdw *= (1.0 - XC_Functional::hybrid_alpha); - v2xdw *= (1.0 - XC_Functional::hybrid_alpha); + sxdw *= (1.0 - XC_Functional::hybrid_alpha); + v1xdw *= (1.0 - XC_Functional::hybrid_alpha); + v2xdw *= (1.0 - XC_Functional::hybrid_alpha); } break; - case XC_GGA_X_PBE_SOL: //PBXsol + } + case XC_GGA_X_PBE_SOL: + { + //PBXsol if (rhoup > small && sqrt(fabs(grhoup2)) > small) { XC_Functional::pbex(2.0 * rhoup, 4.0 * grhoup2, 2, sxup, v1xup, v2xup); @@ -207,14 +310,20 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double XC_Functional::pbex(2.0 * rhodw, 4.0 * grhodw2, 2, sxdw, v1xdw, v2xdw); } break; + } default: - sxup = 0.0; sxdw = 0.0; - v1xup = 0.0; v2xup = 0.0; - v1xdw = 0.0; v2xdw = 0.0; + { + sxup = 0.0; + sxdw = 0.0; + v1xup = 0.0; + v2xup = 0.0; + v1xdw = 0.0; + v2xdw = 0.0; + } } sx = 0.50 * (sxup + sxdw); v2xup = 2.0 * v2xup; - v2xdw = 2.0 * v2xdw; + v2xdw = 2.0 * v2xdw; //} return; @@ -222,8 +331,14 @@ void XC_Functional::gcx_spin(double rhoup, double rhodw, double grhoup2, double // //----------------------------------------------------------------------- -void XC_Functional::gcc_spin(double rho, double &zeta, double grho, double &sc, - double &v1cup, double &v1cdw, double &v2c) +void XC_Functional::gcc_spin( + double rho, + double &zeta, + double grho, + double &sc, + double &v1cup, + double &v1cdw, + double &v2c) { //------------------------------------------------------------------- // gradient corrections for correlations - Hartree a.u. @@ -245,13 +360,14 @@ void XC_Functional::gcc_spin(double rho, double &zeta, double grho, double &sc, // parameter : double small = 1.0e-10; - double epsr = 1.0e-6; - + double epsr = 1.0e-6; double x = 0.0; sc = 0.00; - v1cup = 0.00; v1cdw = 0.00; + v1cup = 0.00; + v1cdw = 0.00; v2c = 0.00; + if (std::abs(zeta) - 1.0 > small || rho <= small || sqrt(std::abs(grho)) <= small) { return; @@ -261,14 +377,14 @@ void XC_Functional::gcc_spin(double rho, double &zeta, double grho, double &sc, // ... ( - 1.0 + epsr ) < zeta < ( 1.0 - epsr ) // zeta = SIGN( MIN( ABS( zeta ), ( 1.D0 - epsr ) ) , zeta ) x = std::min(std::abs(zeta), (1.0 - epsr)); - if(zeta>0) - { - zeta = x; - } - else - { - zeta = -x; - } + if(zeta>0) + { + zeta = x; + } + else + { + zeta = -x; + } } //endif if(func_id[0]==XC_HYB_GGA_XC_PBEH) @@ -281,16 +397,31 @@ void XC_Functional::gcc_spin(double rho, double &zeta, double grho, double &sc, //{ int id = func_id[1]; switch( id ) - { - case XC_GGA_C_P86: //P86 - XC_Functional::perdew86_spin(rho, zeta, grho, sc, v1cup, v1cdw, v2c);break; - case XC_GGA_C_PW91: //PW91_C - ModuleBase::WARNING_QUIT("xc_wrapper_gcxc","there seems to be something wrong with ggac_spin, better use libxc version instead");break; - //XC_Functional::ggac_spin(rho, zeta, grho, sc, v1cup, v1cdw, v2c); - case XC_GGA_C_PBE: //PBC - XC_Functional::pbec_spin(rho, zeta, grho, 1, sc, v1cup, v1cdw, v2c);break; - case XC_GGA_C_PBE_SOL: //PBCsol - XC_Functional::pbec_spin(rho, zeta, grho, 2, sc, v1cup, v1cdw, v2c);break; + { + case XC_GGA_C_P86: + { + //P86 + XC_Functional::perdew86_spin(rho, zeta, grho, sc, v1cup, v1cdw, v2c); + break; + } + case XC_GGA_C_PW91: + { + //PW91_C + ModuleBase::WARNING_QUIT("xc_wrapper_gcxc","there seems to be something wrong with ggac_spin, better use libxc version instead"); + break; + } + case XC_GGA_C_PBE: + { + //PBC + XC_Functional::pbec_spin(rho, zeta, grho, 1, sc, v1cup, v1cdw, v2c); + break; + } + case XC_GGA_C_PBE_SOL: + { + //PBCsol + XC_Functional::pbec_spin(rho, zeta, grho, 2, sc, v1cup, v1cdw, v2c); + break; + } } //} diff --git a/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp b/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp index 255b6966f39..db15f92eca5 100644 --- a/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp +++ b/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp @@ -41,9 +41,11 @@ void XC_Functional::xc( case XC_GGA_X_WC: case XC_GGA_X_B88: case XC_GGA_X_PW91: + { // SLA,PBX,rPBX,PBXsol,WC,B88,PW91_X XC_Functional::slater(rs, e, v); break; + } // Exchange functionals containing attenuated slater exchange case XC_HYB_GGA_XC_PBEH: @@ -67,27 +69,34 @@ void XC_Functional::xc( case XC_GGA_C_PW91: case XC_LDA_C_PW: case XC_GGA_C_PBE_SOL: + { // PBC,PW91,PWLDA XC_Functional::pw(rs, 0, e, v); break; + } // Correlation functionals containing PZ correlation case XC_LDA_C_PZ: case XC_GGA_C_P86: + { // PZ,P86 XC_Functional::pz(rs, 0, e, v); break; + } // Correlation functionals containing LYP correlation case XC_GGA_C_LYP: + { // BLYP XC_Functional::lyp(rs, e, v); break; + } default: { e = 0.0; v = 0.0; + break; } } exc += e; @@ -127,9 +136,11 @@ void XC_Functional::xc_spin( case XC_GGA_X_WC: case XC_GGA_X_B88: case XC_GGA_X_PW91: + { // SLA,PBX,rPBX,PBXsol,WC,B88,PW91_X XC_Functional::slater_spin(rho, zeta, e, vup, vdw); break; + } // Exchange functionals containing attenuated slater exchange case XC_HYB_GGA_XC_PBEH: @@ -155,17 +166,21 @@ void XC_Functional::xc_spin( // Correlation functionals containing PZ correlation case XC_LDA_C_PZ: case XC_GGA_C_P86: + { // PZ,P86 XC_Functional::pz_spin(rs, zeta, e, vup, vdw); break; + } - // Correlation functionals containing PW correlationtests/integrate/101_PW_OU_pseudopot + // Correlation functionals containing PW correlation case XC_GGA_C_PBE: case XC_GGA_C_PBE_SOL: case XC_LDA_C_PW: + { // PBC,PBCsol XC_Functional::pw_spin(rs, zeta, e, vup, vdw); break; + } // Cases that are only realized in LIBXC default: From c6e0ad3fbab5ee3e4aff69dc23266d975be47e75 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 30 May 2026 08:12:02 +0800 Subject: [PATCH 4/6] modify --- source/source_hamilt/module_xc/CMakeLists.txt | 30 +++++++++---------- .../{xc_functional_libxc.h => libxc.h} | 0 ...xc_wrapper_gcxc.cpp => libxc_gga_wrap.cpp} | 2 +- ...ibxc_wrapper_xc.cpp => libxc_lda_wrap.cpp} | 2 +- ..._wrapper_tauxc.cpp => libxc_mgga_wrap.cpp} | 2 +- ...functional_libxc_vxc.cpp => libxc_pot.cpp} | 2 +- ...c_functional_libxc.cpp => libxc_setup.cpp} | 2 +- ...tional_libxc_tools.cpp => libxc_tools.cpp} | 2 +- .../source_hamilt/module_xc/test/test_xc.cpp | 2 +- .../source_hamilt/module_xc/test/test_xc2.cpp | 2 +- .../source_hamilt/module_xc/test/test_xc4.cpp | 2 +- .../source_hamilt/module_xc/test/test_xc5.cpp | 2 +- .../source_hamilt/module_xc/xc_functional.cpp | 2 +- .../source_hamilt/module_xc/xc_functional.h | 2 +- ...{xc_funct_corr_gga.cpp => xc_gga_corr.cpp} | 0 ...{xc_funct_exch_gga.cpp => xc_gga_exch.cpp} | 0 ...ional_wrapper_gcxc.cpp => xc_gga_wrap.cpp} | 2 +- ...xc_functional_gradcorr.cpp => xc_grad.cpp} | 2 +- .../{xc_funct_hcth.cpp => xc_hcth.cpp} | 0 .../module_xc/{xc_funcs.h => xc_ids.h} | 0 ...{xc_funct_corr_lda.cpp => xc_lda_corr.cpp} | 0 ...{xc_funct_exch_lda.cpp => xc_lda_exch.cpp} | 0 ...ctional_wrapper_xc.cpp => xc_lda_wrap.cpp} | 0 .../{xc_functional_vxc.cpp => xc_pot.cpp} | 2 +- 24 files changed, 30 insertions(+), 30 deletions(-) rename source/source_hamilt/module_xc/{xc_functional_libxc.h => libxc.h} (100%) rename source/source_hamilt/module_xc/{xc_functional_libxc_wrapper_gcxc.cpp => libxc_gga_wrap.cpp} (98%) rename source/source_hamilt/module_xc/{xc_functional_libxc_wrapper_xc.cpp => libxc_lda_wrap.cpp} (96%) rename source/source_hamilt/module_xc/{xc_functional_libxc_wrapper_tauxc.cpp => libxc_mgga_wrap.cpp} (99%) rename source/source_hamilt/module_xc/{xc_functional_libxc_vxc.cpp => libxc_pot.cpp} (99%) rename source/source_hamilt/module_xc/{xc_functional_libxc.cpp => libxc_setup.cpp} (97%) rename source/source_hamilt/module_xc/{xc_functional_libxc_tools.cpp => libxc_tools.cpp} (99%) rename source/source_hamilt/module_xc/{xc_funct_corr_gga.cpp => xc_gga_corr.cpp} (100%) rename source/source_hamilt/module_xc/{xc_funct_exch_gga.cpp => xc_gga_exch.cpp} (100%) rename source/source_hamilt/module_xc/{xc_functional_wrapper_gcxc.cpp => xc_gga_wrap.cpp} (99%) rename source/source_hamilt/module_xc/{xc_functional_gradcorr.cpp => xc_grad.cpp} (99%) rename source/source_hamilt/module_xc/{xc_funct_hcth.cpp => xc_hcth.cpp} (100%) rename source/source_hamilt/module_xc/{xc_funcs.h => xc_ids.h} (100%) rename source/source_hamilt/module_xc/{xc_funct_corr_lda.cpp => xc_lda_corr.cpp} (100%) rename source/source_hamilt/module_xc/{xc_funct_exch_lda.cpp => xc_lda_exch.cpp} (100%) rename source/source_hamilt/module_xc/{xc_functional_wrapper_xc.cpp => xc_lda_wrap.cpp} (100%) rename source/source_hamilt/module_xc/{xc_functional_vxc.cpp => xc_pot.cpp} (99%) diff --git a/source/source_hamilt/module_xc/CMakeLists.txt b/source/source_hamilt/module_xc/CMakeLists.txt index 9d38b59fa2a..5f310d09d5b 100644 --- a/source/source_hamilt/module_xc/CMakeLists.txt +++ b/source/source_hamilt/module_xc/CMakeLists.txt @@ -2,21 +2,21 @@ add_library( xc_ OBJECT xc_functional.cpp - xc_functional_vxc.cpp - xc_functional_gradcorr.cpp - xc_functional_wrapper_xc.cpp - xc_functional_wrapper_gcxc.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 + xc_pot.cpp + xc_grad.cpp + xc_lda_wrap.cpp + xc_gga_wrap.cpp + xc_lda_exch.cpp + xc_lda_corr.cpp + xc_gga_exch.cpp + xc_gga_corr.cpp + xc_hcth.cpp + libxc_setup.cpp + libxc_tools.cpp + libxc_pot.cpp + libxc_lda_wrap.cpp + libxc_gga_wrap.cpp + libxc_mgga_wrap.cpp exx_info.cpp ) diff --git a/source/source_hamilt/module_xc/xc_functional_libxc.h b/source/source_hamilt/module_xc/libxc.h similarity index 100% rename from source/source_hamilt/module_xc/xc_functional_libxc.h rename to source/source_hamilt/module_xc/libxc.h diff --git a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_gcxc.cpp b/source/source_hamilt/module_xc/libxc_gga_wrap.cpp similarity index 98% rename from source/source_hamilt/module_xc/xc_functional_libxc_wrapper_gcxc.cpp rename to source/source_hamilt/module_xc/libxc_gga_wrap.cpp index a758d082597..f25c1f36484 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_gcxc.cpp +++ b/source/source_hamilt/module_xc/libxc_gga_wrap.cpp @@ -1,6 +1,6 @@ #ifdef USE_LIBXC -#include "xc_functional_libxc.h" +#include "libxc.h" #include #include diff --git a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_xc.cpp b/source/source_hamilt/module_xc/libxc_lda_wrap.cpp similarity index 96% rename from source/source_hamilt/module_xc/xc_functional_libxc_wrapper_xc.cpp rename to source/source_hamilt/module_xc/libxc_lda_wrap.cpp index f411b737ee2..e4961ccf604 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_xc.cpp +++ b/source/source_hamilt/module_xc/libxc_lda_wrap.cpp @@ -1,6 +1,6 @@ #ifdef USE_LIBXC -#include "xc_functional_libxc.h" +#include "libxc.h" void XC_Functional_Libxc::xc_spin_libxc( const std::vector &func_id, diff --git a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_tauxc.cpp b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp similarity index 99% rename from source/source_hamilt/module_xc/xc_functional_libxc_wrapper_tauxc.cpp rename to source/source_hamilt/module_xc/libxc_mgga_wrap.cpp index 238fe898a6d..d474085190c 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc_wrapper_tauxc.cpp +++ b/source/source_hamilt/module_xc/libxc_mgga_wrap.cpp @@ -6,7 +6,7 @@ #include "source_hamilt/module_xc/exx_info.h" // use GlobalC::exx_info #include "source_hamilt/module_xc/xc_functional.h" -#include "xc_functional_libxc.h" +#include "libxc.h" #include //tau_xc and tau_xc_spin: interface for calling xc_mgga_exc_vxc from LIBXC diff --git a/source/source_hamilt/module_xc/xc_functional_libxc_vxc.cpp b/source/source_hamilt/module_xc/libxc_pot.cpp similarity index 99% rename from source/source_hamilt/module_xc/xc_functional_libxc_vxc.cpp rename to source/source_hamilt/module_xc/libxc_pot.cpp index 2170357ac8b..2e9ea3facfe 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc_vxc.cpp +++ b/source/source_hamilt/module_xc/libxc_pot.cpp @@ -1,7 +1,7 @@ #ifdef USE_LIBXC #include "xc_functional.h" -#include "xc_functional_libxc.h" +#include "libxc.h" #include "source_estate/module_charge/charge.h" #include "source_base/global_variable.h" #include "source_io/module_parameter/parameter.h" diff --git a/source/source_hamilt/module_xc/xc_functional_libxc.cpp b/source/source_hamilt/module_xc/libxc_setup.cpp similarity index 97% rename from source/source_hamilt/module_xc/xc_functional_libxc.cpp rename to source/source_hamilt/module_xc/libxc_setup.cpp index e94f3b95906..18ec7648b7f 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc.cpp +++ b/source/source_hamilt/module_xc/libxc_setup.cpp @@ -1,6 +1,6 @@ #ifdef USE_LIBXC -#include "xc_functional_libxc.h" +#include "libxc.h" #include "source_io/module_parameter/parameter.h" #include "source_base/tool_quit.h" #include "source_base/formatter.h" diff --git a/source/source_hamilt/module_xc/xc_functional_libxc_tools.cpp b/source/source_hamilt/module_xc/libxc_tools.cpp similarity index 99% rename from source/source_hamilt/module_xc/xc_functional_libxc_tools.cpp rename to source/source_hamilt/module_xc/libxc_tools.cpp index 4d4af83a7ed..93327a36226 100644 --- a/source/source_hamilt/module_xc/xc_functional_libxc_tools.cpp +++ b/source/source_hamilt/module_xc/libxc_tools.cpp @@ -1,6 +1,6 @@ #ifdef USE_LIBXC -#include "xc_functional_libxc.h" +#include "libxc.h" #include "xc_functional.h" #include "source_estate/module_charge/charge.h" #include "source_io/module_parameter/parameter.h" diff --git a/source/source_hamilt/module_xc/test/test_xc.cpp b/source/source_hamilt/module_xc/test/test_xc.cpp index 14a1e1c5aa4..07dbfda10ba 100644 --- a/source/source_hamilt/module_xc/test/test_xc.cpp +++ b/source/source_hamilt/module_xc/test/test_xc.cpp @@ -1,6 +1,6 @@ #include "gtest/gtest.h" #include "../xc_functional.h" -#include "../xc_functional_libxc.h" +#include "../libxc.h" #include "../exx_info.h" #include "xctest.h" diff --git a/source/source_hamilt/module_xc/test/test_xc2.cpp b/source/source_hamilt/module_xc/test/test_xc2.cpp index e6135981258..4527828b7e8 100644 --- a/source/source_hamilt/module_xc/test/test_xc2.cpp +++ b/source/source_hamilt/module_xc/test/test_xc2.cpp @@ -1,7 +1,7 @@ #include "gtest/gtest.h" #include "xctest.h" #include "../xc_functional.h" -#include "../xc_functional_libxc.h" +#include "../libxc.h" #include "../exx_info.h" /************************************************ * unit test of functionals diff --git a/source/source_hamilt/module_xc/test/test_xc4.cpp b/source/source_hamilt/module_xc/test/test_xc4.cpp index d1fddc8bdf1..72768c11391 100644 --- a/source/source_hamilt/module_xc/test/test_xc4.cpp +++ b/source/source_hamilt/module_xc/test/test_xc4.cpp @@ -1,5 +1,5 @@ #include "../xc_functional.h" -#include "../xc_functional_libxc.h" +#include "../libxc.h" #include "gtest/gtest.h" #include "xctest.h" #include "../exx_info.h" diff --git a/source/source_hamilt/module_xc/test/test_xc5.cpp b/source/source_hamilt/module_xc/test/test_xc5.cpp index 171e1c263ce..40e36d63c15 100644 --- a/source/source_hamilt/module_xc/test/test_xc5.cpp +++ b/source/source_hamilt/module_xc/test/test_xc5.cpp @@ -1,5 +1,5 @@ #include "../xc_functional.h" -#include "../xc_functional_libxc.h" +#include "../libxc.h" #include "gtest/gtest.h" #define private public #include "source_io/module_parameter/parameter.h" diff --git a/source/source_hamilt/module_xc/xc_functional.cpp b/source/source_hamilt/module_xc/xc_functional.cpp index acef0f6615f..e32f9c0735b 100644 --- a/source/source_hamilt/module_xc/xc_functional.cpp +++ b/source/source_hamilt/module_xc/xc_functional.cpp @@ -4,7 +4,7 @@ #include "source_base/tool_title.h" #ifdef USE_LIBXC -#include "xc_functional_libxc.h" +#include "libxc.h" #endif XC_Functional::XC_Functional(){} diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index 7322ea3b46f..b620384bfd9 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -8,7 +8,7 @@ #ifdef USE_LIBXC #include #else -#include "xc_funcs.h" +#include "xc_ids.h" #endif // ifdef USE_LIBXC #include "source_base/macros.h" #include "source_base/global_function.h" diff --git a/source/source_hamilt/module_xc/xc_funct_corr_gga.cpp b/source/source_hamilt/module_xc/xc_gga_corr.cpp similarity index 100% rename from source/source_hamilt/module_xc/xc_funct_corr_gga.cpp rename to source/source_hamilt/module_xc/xc_gga_corr.cpp diff --git a/source/source_hamilt/module_xc/xc_funct_exch_gga.cpp b/source/source_hamilt/module_xc/xc_gga_exch.cpp similarity index 100% rename from source/source_hamilt/module_xc/xc_funct_exch_gga.cpp rename to source/source_hamilt/module_xc/xc_gga_exch.cpp diff --git a/source/source_hamilt/module_xc/xc_functional_wrapper_gcxc.cpp b/source/source_hamilt/module_xc/xc_gga_wrap.cpp similarity index 99% rename from source/source_hamilt/module_xc/xc_functional_wrapper_gcxc.cpp rename to source/source_hamilt/module_xc/xc_gga_wrap.cpp index cfaacd812a0..dcb7729b66a 100644 --- a/source/source_hamilt/module_xc/xc_functional_wrapper_gcxc.cpp +++ b/source/source_hamilt/module_xc/xc_gga_wrap.cpp @@ -11,7 +11,7 @@ #include "source_base/global_function.h" #ifdef USE_LIBXC -#include "xc_functional_libxc.h" +#include "libxc.h" #endif void XC_Functional::gcxc( diff --git a/source/source_hamilt/module_xc/xc_functional_gradcorr.cpp b/source/source_hamilt/module_xc/xc_grad.cpp similarity index 99% rename from source/source_hamilt/module_xc/xc_functional_gradcorr.cpp rename to source/source_hamilt/module_xc/xc_grad.cpp index 9546ddaf51c..6414f58da8c 100644 --- a/source/source_hamilt/module_xc/xc_functional_gradcorr.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -18,7 +18,7 @@ #include #ifdef USE_LIBXC -#include "xc_functional_libxc.h" +#include "libxc.h" #endif // from gradcorr.f90 diff --git a/source/source_hamilt/module_xc/xc_funct_hcth.cpp b/source/source_hamilt/module_xc/xc_hcth.cpp similarity index 100% rename from source/source_hamilt/module_xc/xc_funct_hcth.cpp rename to source/source_hamilt/module_xc/xc_hcth.cpp diff --git a/source/source_hamilt/module_xc/xc_funcs.h b/source/source_hamilt/module_xc/xc_ids.h similarity index 100% rename from source/source_hamilt/module_xc/xc_funcs.h rename to source/source_hamilt/module_xc/xc_ids.h diff --git a/source/source_hamilt/module_xc/xc_funct_corr_lda.cpp b/source/source_hamilt/module_xc/xc_lda_corr.cpp similarity index 100% rename from source/source_hamilt/module_xc/xc_funct_corr_lda.cpp rename to source/source_hamilt/module_xc/xc_lda_corr.cpp diff --git a/source/source_hamilt/module_xc/xc_funct_exch_lda.cpp b/source/source_hamilt/module_xc/xc_lda_exch.cpp similarity index 100% rename from source/source_hamilt/module_xc/xc_funct_exch_lda.cpp rename to source/source_hamilt/module_xc/xc_lda_exch.cpp diff --git a/source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp b/source/source_hamilt/module_xc/xc_lda_wrap.cpp similarity index 100% rename from source/source_hamilt/module_xc/xc_functional_wrapper_xc.cpp rename to source/source_hamilt/module_xc/xc_lda_wrap.cpp diff --git a/source/source_hamilt/module_xc/xc_functional_vxc.cpp b/source/source_hamilt/module_xc/xc_pot.cpp similarity index 99% rename from source/source_hamilt/module_xc/xc_functional_vxc.cpp rename to source/source_hamilt/module_xc/xc_pot.cpp index 917b915f7b6..d5082f8ed69 100644 --- a/source/source_hamilt/module_xc/xc_functional_vxc.cpp +++ b/source/source_hamilt/module_xc/xc_pot.cpp @@ -10,7 +10,7 @@ #include "xc_functional.h" #ifdef USE_LIBXC -#include "xc_functional_libxc.h" +#include "libxc.h" #endif // [etxc, vtxc, v] = XC_Functional::v_xc(...) From 82ff29100e432f6b75ef5adf45872f38df9d582e Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 30 May 2026 09:38:56 +0800 Subject: [PATCH 5/6] update --- source/source_estate/module_pot/pot_xc.cpp | 2 +- .../module_xc/test/CMakeLists.txt | 42 +++++++++---------- .../source_io/module_chgpot/write_libxc_r.cpp | 2 +- .../module_lr/potentials/xc_kernel.cpp | 2 +- source/source_pw/module_pwdft/forces_cc.cpp | 2 +- source/source_pw/module_pwdft/stress_cc.cpp | 2 +- 6 files changed, 26 insertions(+), 26 deletions(-) diff --git a/source/source_estate/module_pot/pot_xc.cpp b/source/source_estate/module_pot/pot_xc.cpp index e122a9f1ae8..bd26a80e8c8 100644 --- a/source/source_estate/module_pot/pot_xc.cpp +++ b/source/source_estate/module_pot/pot_xc.cpp @@ -4,7 +4,7 @@ #include "source_hamilt/module_xc/xc_functional.h" #ifdef USE_LIBXC -#include "source_hamilt/module_xc/xc_functional_libxc.h" +#include "source_hamilt/module_xc/libxc.h" #endif namespace elecstate diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 93b1d546678..02a4298463c 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -1,22 +1,22 @@ AddTest( TARGET MODULE_HAMILT_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 + SOURCES test_xc.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp ../xc_gga_wrap.cpp ../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp + ../libxc_gga_wrap.cpp ../libxc_setup.cpp ) AddTest( TARGET MODULE_HAMILT_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 ../xc_functional_libxc.cpp + SOURCES test_xc1.cpp ../xc_functional.cpp ../libxc_setup.cpp ) AddTest( TARGET MODULE_HAMILT_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 + SOURCES test_xc2.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp ../xc_gga_wrap.cpp ../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp + ../libxc_gga_wrap.cpp ../libxc_lda_wrap.cpp ../libxc_setup.cpp ) if (USE_CUDA) @@ -28,22 +28,22 @@ endif() AddTest( TARGET MODULE_HAMILT_XCTest_GRADCORR 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_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 - ../../../source_base/matrix.cpp - ../../../source_base/memory_recorder.cpp - ../../../source_base/libm/branred.cpp - ../../../source_base/libm/sincos.cpp - ../../../source_base/module_external/blas_connector_base.cpp ../../../source_base/module_external/blas_connector_vector.cpp ../../../source_base/module_external/blas_connector_matrix.cpp - ../../../source_base/module_fft/fft_bundle.cpp - ../../../source_base/module_fft/fft_cpu.cpp - ${FFT_SRC} + SOURCES test_xc3.cpp ../xc_grad.cpp ../xc_functional.cpp + ../xc_lda_wrap.cpp ../xc_gga_wrap.cpp + ../libxc_setup.cpp + ../libxc_lda_wrap.cpp + ../libxc_gga_wrap.cpp + ../libxc_mgga_wrap.cpp + ../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp + ../xc_lda_exch.cpp ../xc_hcth.cpp + ../../../source_base/matrix.cpp + ../../../source_base/memory_recorder.cpp + ../../../source_base/libm/branred.cpp + ../../../source_base/libm/sincos.cpp + ../../../source_base/module_external/blas_connector_base.cpp ../../../source_base/module_external/blas_connector_vector.cpp ../../../source_base/module_external/blas_connector_matrix.cpp + ../../../source_base/module_fft/fft_bundle.cpp + ../../../source_base/module_fft/fft_cpu.cpp + ${FFT_SRC} ) AddTest( diff --git a/source/source_io/module_chgpot/write_libxc_r.cpp b/source/source_io/module_chgpot/write_libxc_r.cpp index e8b20cc13e3..cc915ff3f81 100644 --- a/source/source_io/module_chgpot/write_libxc_r.cpp +++ b/source/source_io/module_chgpot/write_libxc_r.cpp @@ -8,7 +8,7 @@ #include "write_libxc_r.h" #include "source_base/parallel_comm.h" #include "source_hamilt/module_xc/xc_functional.h" -#include "source_hamilt/module_xc/xc_functional_libxc.h" +#include "source_hamilt/module_xc/libxc.h" #include "source_estate/module_charge/charge.h" #include "source_basis/module_pw/pw_basis_big.h" #include "source_basis/module_pw/pw_basis.h" diff --git a/source/source_lcao/module_lr/potentials/xc_kernel.cpp b/source/source_lcao/module_lr/potentials/xc_kernel.cpp index 06d91710a9d..ff57bd01158 100644 --- a/source/source_lcao/module_lr/potentials/xc_kernel.cpp +++ b/source/source_lcao/module_lr/potentials/xc_kernel.cpp @@ -9,7 +9,7 @@ #include "source_io/module_output/cube_io.h" #ifdef USE_LIBXC #include -#include "source_hamilt/module_xc/xc_functional_libxc.h" +#include "source_hamilt/module_xc/libxc.h" #endif #ifdef _OPENMP #include diff --git a/source/source_pw/module_pwdft/forces_cc.cpp b/source/source_pw/module_pwdft/forces_cc.cpp index 51cc2c04c6c..7fffa85fb11 100644 --- a/source/source_pw/module_pwdft/forces_cc.cpp +++ b/source/source_pw/module_pwdft/forces_cc.cpp @@ -24,7 +24,7 @@ #ifdef USE_LIBXC -#include "source_hamilt/module_xc/xc_functional_libxc.h" +#include "source_hamilt/module_xc/libxc.h" #endif diff --git a/source/source_pw/module_pwdft/stress_cc.cpp b/source/source_pw/module_pwdft/stress_cc.cpp index 2ee46650b00..70f2c0709ea 100644 --- a/source/source_pw/module_pwdft/stress_cc.cpp +++ b/source/source_pw/module_pwdft/stress_cc.cpp @@ -7,7 +7,7 @@ #include "source_estate/cal_ux.h" #ifdef USE_LIBXC -#include "source_hamilt/module_xc/xc_functional_libxc.h" +#include "source_hamilt/module_xc/libxc.h" #endif From 84dfbd241408a10c5b4e5a9a026c28cf25a0077a Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 30 May 2026 09:46:08 +0800 Subject: [PATCH 6/6] update --- source/Makefile.Objects | 28 +++++++------- .../module_surchem/test/CMakeLists.txt | 36 +++++++++--------- .../module_xc/test/CMakeLists.txt | 38 +++++++++---------- .../source_hamilt/module_xc/xc_functional.h | 16 ++++---- source/source_hamilt/module_xc/xc_grad.cpp | 2 +- 5 files changed, 60 insertions(+), 60 deletions(-) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 6e1adc3f448..04cbff62277 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -502,20 +502,20 @@ OBJS_SYMMETRY=symm_other.o\ OBJS_XC=xc_functional.o\ xc_functional_op.o\ - xc_functional_vxc.o\ - xc_functional_gradcorr.o\ - xc_functional_wrapper_xc.o\ - xc_functional_wrapper_gcxc.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\ - xc_funct_corr_gga.o\ + xc_pot.o\ + xc_grad.o\ + xc_lda_wrap.o\ + xc_gga_wrap.o\ + libxc_setup.o\ + libxc_tools.o\ + libxc_pot.o\ + libxc_lda_wrap.o\ + libxc_gga_wrap.o\ + libxc_mgga_wrap.o\ + xc_lda_exch.o\ + xc_lda_corr.o\ + xc_gga_exch.o\ + xc_gga_corr.o\ xc_funct_hcth.o\ exx_info.o\ diff --git a/source/source_hamilt/module_surchem/test/CMakeLists.txt b/source/source_hamilt/module_surchem/test/CMakeLists.txt index eb667b7c2ab..e40dca59141 100644 --- a/source/source_hamilt/module_surchem/test/CMakeLists.txt +++ b/source/source_hamilt/module_surchem/test/CMakeLists.txt @@ -30,28 +30,28 @@ AddTest( TARGET MODULE_HAMILT_surchem_cal_vcav LIBS parameter ${math_libs} planewave device base container SOURCES cal_vcav_test.cpp ../cal_vcav.cpp ../surchem.cpp ../../../source_pw/module_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_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 + ../../module_xc/xc_grad.cpp ../../module_xc/xc_functional.cpp + ../../module_xc/xc_lda_wrap.cpp ../../module_xc/xc_gga_wrap.cpp + ../../module_xc/libxc_setup.cpp + ../../module_xc/libxc_pot.cpp + ../../module_xc/libxc_lda_wrap.cpp + ../../module_xc/libxc_gga_wrap.cpp + ../../module_xc/libxc_mgga_wrap.cpp + ../../module_xc/xc_gga_corr.cpp ../../module_xc/xc_lda_corr.cpp ../../module_xc/xc_gga_exch.cpp + ../../module_xc/xc_lda_exch.cpp ../../module_xc/xc_hcth.cpp ) AddTest( TARGET MODULE_HAMILT_surchem_cal_vel LIBS parameter ${math_libs} planewave device base container SOURCES cal_vel_test.cpp ../cal_vel.cpp ../surchem.cpp ../cal_epsilon.cpp ../minimize_cg.cpp ../../../source_pw/module_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_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 + ../../module_xc/xc_grad.cpp ../../module_xc/xc_functional.cpp + ../../module_xc/xc_lda_wrap.cpp ../../module_xc/xc_gga_wrap.cpp + ../../module_xc/libxc_setup.cpp + ../../module_xc/libxc_pot.cpp + ../../module_xc/libxc_lda_wrap.cpp + ../../module_xc/libxc_gga_wrap.cpp + ../../module_xc/libxc_mgga_wrap.cpp + ../../module_xc/xc_gga_corr.cpp ../../module_xc/xc_lda_corr.cpp ../../module_xc/xc_gga_exch.cpp + ../../module_xc/xc_lda_exch.cpp ../../module_xc/xc_hcth.cpp ) \ No newline at end of file diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 02a4298463c..644f99dde92 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -49,30 +49,30 @@ AddTest( AddTest( TARGET MODULE_HAMILT_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_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 + SOURCES test_xc4.cpp ../xc_functional.cpp ../xc_lda_wrap.cpp + ../xc_gga_wrap.cpp + ../libxc_setup.cpp + ../libxc_lda_wrap.cpp + ../libxc_gga_wrap.cpp + ../libxc_mgga_wrap.cpp + ../xc_gga_corr.cpp ../xc_lda_corr.cpp + ../xc_gga_exch.cpp ../xc_lda_exch.cpp ../xc_hcth.cpp ) AddTest( TARGET MODULE_HAMILT_XCTest_VXC 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_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 + SOURCES test_xc5.cpp ../xc_grad.cpp ../xc_functional.cpp + ../xc_lda_wrap.cpp ../xc_gga_wrap.cpp + ../libxc_setup.cpp + ../libxc_lda_wrap.cpp + ../libxc_gga_wrap.cpp + ../libxc_mgga_wrap.cpp + ../xc_gga_corr.cpp ../xc_lda_corr.cpp ../xc_gga_exch.cpp + ../xc_lda_exch.cpp ../xc_hcth.cpp + ../xc_pot.cpp + ../libxc_pot.cpp + ../libxc_tools.cpp ../../../source_base/module_external/blas_connector_base.cpp ../../../source_base/module_external/blas_connector_vector.cpp ../../../source_base/module_external/blas_connector_matrix.cpp ../../../source_base/matrix.cpp ../../../source_base/memory_recorder.cpp diff --git a/source/source_hamilt/module_xc/xc_functional.h b/source/source_hamilt/module_xc/xc_functional.h index b620384bfd9..87bfa2d4f8e 100644 --- a/source/source_hamilt/module_xc/xc_functional.h +++ b/source/source_hamilt/module_xc/xc_functional.h @@ -32,7 +32,7 @@ class XC_Functional //------------------- //------------------- -// xc_functional_vxc.cpp +// xc_pot.cpp //------------------- // This file contains interface to the xc_functional class @@ -105,7 +105,7 @@ class XC_Functional static std::vector get_func_id() { return func_id; } //------------------- -// xc_functional_wrapper_xc.cpp +// xc_lda_wrap.cpp //------------------- // This file contains wrapper for the LDA functionals @@ -142,7 +142,7 @@ class XC_Functional double &vxcdw); //------------------- -// xc_functional_wrapper_gcxc.cpp +// xc_gga_wrap.cpp //------------------- // This file contains wrapper for the GGA functionals @@ -185,7 +185,7 @@ class XC_Functional double &v2c); //------------------- -// xc_functional_gradcorr.cpp +// xc_grad.cpp //------------------- // This file contains subroutines realted to gradient calculations @@ -239,7 +239,7 @@ class XC_Functional const bool lsign_); //------------------- - // xc_funct_exch_lda.cpp + // xc_lda_exch.cpp //------------------- // This file contains realization of LDA exchange functionals @@ -278,7 +278,7 @@ class XC_Functional double &vxdw); //------------------- -// xc_funct_corr_lda.cpp +// xc_lda_corr.cpp //------------------- // This file contains realization of LDA correlation functionals @@ -321,7 +321,7 @@ class XC_Functional static void pz_polarized(const double &rs, double &ec, double &vc); //------------------- -// xc_funct_exch_gga.cpp +// xc_gga_exch.cpp //------------------- // This file contains realizations of gradient correction to exchange part @@ -373,7 +373,7 @@ class XC_Functional double &v2x); //------------------- -// xc_funct_corr_gga.cpp +// xc_gga_corr.cpp //------------------- // This file contains realizations of gradient correction to correlation part diff --git a/source/source_hamilt/module_xc/xc_grad.cpp b/source/source_hamilt/module_xc/xc_grad.cpp index 6414f58da8c..c4ab666d04b 100644 --- a/source/source_hamilt/module_xc/xc_grad.cpp +++ b/source/source_hamilt/module_xc/xc_grad.cpp @@ -501,7 +501,7 @@ void XC_Functional::gradcorr( } } #ifdef _OPENMP - #pragma omp critical(xc_functional_gradcorr_reduce) + #pragma omp critical(xc_grad_reduce) { if(is_stress) {