//==========================================================
// AUTHOR : mohan
// DATE : 2009-04-04
//==========================================================
#ifndef XC_FUNCTIONAL_H
#define XC_FUNCTIONAL_H

#ifdef USE_LIBXC
#include <xc.h>
#else
#include "xc_funcs.h"
#endif	// ifdef USE_LIBXC
#include "source_base/macros.h"
#include "source_base/global_function.h"
#include "source_base/global_variable.h"
#include "source_base/vector3.h"
#include "source_base/matrix.h"
#include "exx_info.h"
#include "source_basis/module_pw/pw_basis_k.h"
#include "source_estate/module_charge/charge.h"
#include "source_cell/unitcell.h"

#include <map> // added by jghan, 2024-10-10

class XC_Functional
{
	public:

	XC_Functional();
	~XC_Functional();

//-------------------
// subroutines, grouped according to the file they are in:
//-------------------

//-------------------
//  xc_functional_vxc.cpp
//-------------------

// This file contains interface to the xc_functional class
// it includes 3 subroutines:
// 1. v_xc : which takes rho as input, and [etxc, vtxc, v_xc] as output
// 2. v_xc_libxc : which does the same thing as v_xc, but calling libxc
// 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<double,double,ModuleBase::matrix> v_xc(
		const int &nrxx, // number of real-space grid
		const Charge* const chr,
		const UnitCell *ucell); // charge density

//-------------------
//  xc_functional.cpp
//-------------------

// 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
// 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.

    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);

	private:

	static std::vector<int> 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<int, double> scaling_factor_xc;

	public:
	static std::vector<int> get_func_id() { return func_id; }

//-------------------
//  xc_functional_wrapper_xc.cpp
//-------------------

// This file contains wrapper for the LDA functionals
// it includes 3 subroutines:
// 1. xc, which is the wrapper of LDA part
// (i.e. LDA functional and LDA part of GGA functional)
// 2. xc_spin, which is the spin polarized counterpart of xc
// 3. xc_spin_libxc, which is the wrapper for LDA functional, spin polarized

// NOTE : In our own realization of GGA functional, the LDA part
// and gradient correction are calculated separately.
// The LDA part is provided in xc, while the gradient correction is 
// provided in gradcorr through gcxc/gcx_spin+gcc_spin.
// While in LIBXC, the entire GGA functional is provided.
// As a result, xc/xc_spin and xc_spin_libxc are different for GGA,
// the former gives nonzero result, while the latter returns 0.
// Furthermore, the reason for not having xc_libxc is that something like
// xc_libxc which evaluates functional for individual grid points
// is not efficient. For nspin = 1 and 2, v_xc_libxc evaluates potential
// on the entire grid. I'm having xc_spin_libxc because v_xc_libxc
// does not support nspin = 4.

	public:

	// 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);

//-------------------
//  xc_functional_wrapper_gcxc.cpp
//-------------------

// This file contains wrapper for the GGA functionals
// it includes 4 subroutines:
// 1. gcxc, which is the wrapper for gradient correction part
// 2. gcx_spin, spin polarized, exchange only
// 3. gcc_spin, spin polarized, correlation only

// The difference between our realization (gcxc/gcx_spin/gcc_spin) and
// 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);

//-------------------
//  xc_functional_gradcorr.cpp
//-------------------

// This file contains subroutines realted to gradient calculations
// 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
// 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<double>& stress_gga,
                         const bool is_stress = false);
	template <typename T, typename Device,
          typename Real = typename GetTypeReal<T>::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<double>* rhog,
                         ModuleBase::Vector3<double>* gdr,
                         const ModulePW::PW_Basis* rho_basis,
                         const double tpiba);
    static void grad_dot(const ModuleBase::Vector3<double>* 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_);

    //-------------------
    //  xc_funct_exch_lda.cpp
    //-------------------

    // This file contains realization of LDA exchange functionals
    // Spin unpolarized ones:
    //  1. slater: ordinary Slater exchange with alpha=2/3
    //  2. slater1: Slater exchange with alpha=1
    //  3. slater_rxc : Slater exchange with alpha=2/3 and Relativistic exchange
    // And their spin polarized counterparts:
    //  1. slater_spin
    //  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);

//-------------------
//  xc_funct_corr_lda.cpp
//-------------------

// This file contains realization of LDA correlation functionals
// Spin unpolarized ones:
//  1. pw : Perdew-Wang LDA correlation
//  2. pz : Perdew-Zunger LDA correlation
//  3. lyp : Lee-Yang-Parr correlation
//  4. vwn : Vosko-Wilk-Nusair LDA correlation
//  5. wigner : Wigner
//  6. hl : Hedin-Lunqvist
//  7. gl : Gunnarson-Lunqvist
// And some of their spin polarized counterparts:
//  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);

//-------------------
//  xc_funct_exch_gga.cpp
//-------------------

// This file contains realizations of gradient correction to exchange part
// Spin unpolarized ones:
//  1. becke88 : Becke88 exchange
//  2. ggax : PW91 exchange
//  3. pbex : PBE exchange (and revPBE)
//  4. optx : OPTX, Handy et al.
//  5. wcx : Wu-Cohen exchange
// 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);

//-------------------
//  xc_funct_corr_gga.cpp
//-------------------

// This file contains realizations of gradient correction to correlation part
// Spin unpolarized ones:
//  1. perdew86 : P86
//  2. ggac : PW91
//  3. pbec
//  4. glyp
// And some of their spin polarized counterparts:
//  1. perdew86_spin
//  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);

//-------------------
//  xc_funct_hcth.cpp
//-------------------
// 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);

};

#endif //XC_FUNCTION_H
