#include "numerical_basis.h"

#include "source_io/module_parameter/parameter.h"
#include "source_base/constants.h"
#include "source_base/global_variable.h"
#include "source_base/intarray.h"
#include "source_base/math_ylmreal.h"
#include "source_base/parallel_reduce.h"
#include "source_base/timer.h"
#include "source_base/vector3.h"
#include "source_cell/module_symmetry/symmetry.h"
#include "source_pw/module_pwdft/global.h"
#include "source_io/numerical_basis_jyjy.h"
#include "winput.h"

#include 
#include 
#include 
#include 
Numerical_Basis::Numerical_Basis()
{
}
Numerical_Basis::~Numerical_Basis()
{
}

//============================================================
// MEMBER FUNCTION :
// NAME : init
// DESCRIPTION : Two main functions:
// (1) start_from_file = true;
// Firstly, use check(1) to call bessel_basis.init
// to generate TableOne.
// Secondly readin C4 from file.
// Thirdly generate 3D atomic wfc in G space, put the
// results in psi.
//
// (2) If output overlap Q, start_from_file = false;
// Firstly, use check(0) to call bessel_basis,init
// to generate TableOne
// Secondly output overlap, use psi(evc) and jlq3d.
//============================================================
// void Numerical_Basis::start_from_file_k(const int& ik, ModuleBase::ComplexMatrix& psi, const Structure_Factor& sf,
//                                         const ModulePW::PW_Basis_K* wfcpw, const UnitCell& ucell)
// {
//     ModuleBase::TITLE("Numerical_Basis", "start_from_file_k");

//     if (!this->init_label)
//     {
//         // true stands for : start_from_file
//         this->bessel_basis.init(true, std::stod(PARAM.inp.bessel_nao_ecut), ucell.ntype, ucell.lmax,
//                                 PARAM.inp.bessel_nao_smooth, PARAM.inp.bessel_nao_sigma, PARAM.globalv.bessel_nao_rcut,
//                                 PARAM.inp.bessel_nao_tolerence, ucell);
//         this->mu_index = this->init_mu_index(ucell);
//         this->init_label = true;
//     }
//     this->numerical_atomic_wfc(ik, wfcpw, psi, sf, ucell);
// }

// The function is called in run_fp.cpp.
void Numerical_Basis::output_overlap(const psi::Psi& psi,
                                     const Structure_Factor& sf,
                                     const K_Vectors& kv,
                                     const ModulePW::PW_Basis_K* wfcpw,
                                     const UnitCell& ucell,
                                     const int& index)
{
    ModuleBase::TITLE("Numerical_Basis", "output_overlap");
    ModuleBase::GlobalFunc::NEW_PART("Overlap Data For Spillage Minimization");
    const double bessel_nao_rcut = PARAM.inp.bessel_nao_rcuts[index];

    //---------------------------------------------------------
    // if the numerical_basis hasn't been initialized yet,
    // then we initial here.
    //---------------------------------------------------------
    if (!this->init_label)
    {
        // false stands for : 'Faln' is not used.
        this->bessel_basis.init(false, std::stod(PARAM.inp.bessel_nao_ecut), ucell.ntype, ucell.lmax,
                                PARAM.inp.bessel_nao_smooth, PARAM.inp.bessel_nao_sigma, bessel_nao_rcut,
                                PARAM.inp.bessel_nao_tolerence, ucell);
        this->mu_index = this->init_mu_index(ucell);
        this->init_label = true;
    }
    ModuleBase::GlobalFunc::MAKE_DIR(winput::spillage_outdir);
    for (int derivative_order = 0; derivative_order