esolver_dm2rho.cpp

#include "esolver_dm2rho.h"

#include "source_base/timer.h"
#include "source_cell/module_neighbor/sltk_atom_arrange.h"
#include "source_cell/read_pp_ucell.h"
#include "source_estate/elecstate_lcao.h"
#include "source_io/module_ml/io_npz.h"
#include "source_io/module_output/cube_io.h"
#include "source_lcao/lcao_domain.h"
#include "source_lcao/hamilt_lcao.h"
#include "source_lcao/module_operator_lcao/operator_lcao.h"
#include "source_lcao/rho_tau_lcao.h" // mohan add 2025-10-24

namespace ModuleESolver
{

template <typename TK, typename TR>
ESolver_DM2rho<TK, TR>::ESolver_DM2rho()
{
    this->classname = "ESolver_DM2rho";
    this->basisname = "LCAO";
}

template <typename TK, typename TR>
ESolver_DM2rho<TK, TR>::~ESolver_DM2rho()
{
}

template <typename TK, typename TR>
void ESolver_DM2rho<TK, TR>::before_all_runners(BaseCell& basecell, const Input_para& inp)
{
    basecell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__);
    UnitCell& ucell = static_cast<UnitCell&>(basecell);

    ModuleBase::TITLE("ESolver_DM2rho", "before_all_runners");
    ModuleBase::timer::start("ESolver_DM2rho", "before_all_runners");

    ESolver_KS_LCAO<TK, TR>::before_all_runners(ucell, inp);

    ModuleBase::timer::end("ESolver_DM2rho", "before_all_runners");
}

template <typename TK, typename TR>
void ESolver_DM2rho<TK, TR>::runner(BaseCell& basecell, const int istep)
{
    basecell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__);
    UnitCell& ucell = static_cast<UnitCell&>(basecell);

    ModuleBase::TITLE("ESolver_DM2rho", "runner");
    ModuleBase::timer::start("ESolver_DM2rho", "runner");

    ESolver_KS_LCAO<TK, TR>::before_scf(ucell, istep);

    // file name of DM
    std::string zipname = "output_DM0.npz";

    // read DM from file
    ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(1)));

    // if nspin=2, need extra reading
    if (this->inp_->nspin == 2)
    {
        zipname = "output_DM1.npz";
        ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(2)));
    }

    // it's dangerous to design psiToRho function like this, mohan note 20251024
    // this->pelec->psiToRho(*this->psi);
    LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), this->inp_->nspin, &this->chr);

    int nspin0 = this->inp_->nspin == 2 ? 2 : 1;

    for (int is = 0; is < nspin0; is++)
    {
        std::string fn = PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";

        // write electron density
        ModuleIO::write_vdata_palgrid(this->Pgrid,
                                      this->chr.rho[is],
                                      is,
                                      this->inp_->nspin,
                                      istep,
                                      fn,
                                      this->pelec->eferm.get_efval(is),
                                      &(ucell),
                                      3,
                                      1,
                                      PARAM.globalv.two_fermi,
                                      false);
    }

    ModuleBase::timer::end("ESolver_DM2rho", "runner");
}

template <typename TK, typename TR>
void ESolver_DM2rho<TK, TR>::after_all_runners(BaseCell& basecell)
{
    basecell.require_kind(BaseCell::Kind::unit_cell, __FUNCTION__);
    UnitCell& ucell = static_cast<UnitCell&>(basecell);

    ModuleBase::TITLE("ESolver_DM2rho", "after_all_runners");
    ModuleBase::timer::start("ESolver_DM2rho", "after_all_runners");

    ESolver_KS_LCAO<TK, TR>::after_all_runners(ucell);

    ModuleBase::timer::end("ESolver_DM2rho", "after_all_runners");
};

template class ESolver_DM2rho<std::complex<double>, double>;
template class ESolver_DM2rho<std::complex<double>, std::complex<double>>;

} // namespace ModuleESolver