diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 2278ddb87c..64c96f357b 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -191,6 +191,7 @@ - [DeePKS](#deepks) - [deepks\_out\_labels](#deepks_out_labels) - [deepks\_out\_freq\_elec](#deepks_out_freq_elec) + - [deepks\_out\_base](#deepks_out_base) - [deepks\_scf](#deepks_scf) - [deepks\_equiv](#deepks_equiv) - [deepks\_model](#deepks_model) @@ -2173,6 +2174,13 @@ Warning: this function is not robust enough for the current version. Please try - **Description**: When `deepks_out_freq_elec` is greater than 0, print labels and descriptors for DeePKS in OUT.${suffix}/DeePKS_Labels_Elec per `deepks_out_freq_elec` electronic iterations, with suffix `_e*` to distinguish different steps. Often used with `deepks_out_labels` equals 1. - **Default**: 0 +### deepks_out_base + +- **Type**: String +- **Availability**: Numerical atomic orbital basis and `deepks_out_freq_elec` is greater than 0 +- **Description**: Print labels and descriptors calculated by base functional ( determined by `deepks_out_base` ) and target functional ( determined by `dft_functional` ) for DeePKS in per `deepks_out_freq_elec` electronic iterations. The SCF process, labels and descriptors output of the target functional are all consistent with those when the target functional is used alone. The only additional output under this configuration is the labels of the base functional. Often used with `deepks_out_labels` equals 1. +- **Default**: None + ### deepks_scf - **Type**: Boolean diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 6b4e6cd52e..0e2ed91efa 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -275,6 +275,7 @@ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ esolver_gets.o\ lcao_others.o\ esolver_dm2rho.o\ + esolver_double_xc.o\ OBJS_GINT=gint_old.o\ gint_gamma_env.o\ diff --git a/source/source_esolver/CMakeLists.txt b/source/source_esolver/CMakeLists.txt index 281b2e8b55..f58e47b3aa 100644 --- a/source/source_esolver/CMakeLists.txt +++ b/source/source_esolver/CMakeLists.txt @@ -22,6 +22,7 @@ if(ENABLE_LCAO) esolver_gets.cpp lcao_others.cpp esolver_dm2rho.cpp + esolver_double_xc.cpp ) endif() diff --git a/source/source_esolver/esolver.cpp b/source/source_esolver/esolver.cpp index 1849fbcac0..5d513b398c 100644 --- a/source/source_esolver/esolver.cpp +++ b/source/source_esolver/esolver.cpp @@ -6,6 +6,7 @@ #include "source_io/module_parameter/parameter.h" #ifdef __LCAO #include "esolver_dm2rho.h" +#include "esolver_double_xc.h" #include "esolver_gets.h" #include "esolver_ks_lcao.h" #include "esolver_ks_lcao_tddft.h" @@ -197,7 +198,14 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell) } if (PARAM.globalv.gamma_only_local) { - return new ESolver_KS_LCAO(); + if (PARAM.inp.deepks_out_base != "none") + { + return new ESolver_DoubleXC(); + } + else + { + return new ESolver_KS_LCAO(); + } } else if (PARAM.inp.nspin < 4) { @@ -205,6 +213,10 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell) { return new ESolver_DM2rho, double>(); } + else if (PARAM.inp.deepks_out_base != "none") + { + return new ESolver_DoubleXC, double>(); + } else { return new ESolver_KS_LCAO, double>(); @@ -216,6 +228,10 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell) { return new ESolver_DM2rho, std::complex>(); } + else if (PARAM.inp.deepks_out_base != "none") + { + return new ESolver_DoubleXC, std::complex>(); + } else { return new ESolver_KS_LCAO, std::complex>(); diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp new file mode 100644 index 0000000000..6c8eb1cbb2 --- /dev/null +++ b/source/source_esolver/esolver_double_xc.cpp @@ -0,0 +1,436 @@ +#include "esolver_double_xc.h" +#include "source_hamilt/module_xc/xc_functional.h" +#include "source_hamilt/module_ewald/H_Ewald_pw.h" +#include "source_hamilt/module_vdw/vdw.h" +#ifdef __MLALGO +#include "source_lcao/module_deepks/LCAO_deepks.h" +#include "source_lcao/module_deepks/LCAO_deepks_interface.h" +#endif +//-----force& stress------------------- +#include "source_lcao/FORCE_STRESS.h" + +//-----HSolver ElecState Hamilt-------- +#include "source_estate/elecstate_lcao.h" +#include "source_estate/elecstate_tools.h" +#include "source_lcao/hamilt_lcao.h" +#include "source_hsolver/hsolver_lcao.h" +#include "source_io/module_parameter/parameter.h" +#include "source_io/write_HS.h" // use ModuleIO::write_hsk() + +namespace ModuleESolver +{ + +template +ESolver_DoubleXC::ESolver_DoubleXC() +{ + this->classname = "ESolver_DoubleXC"; + this->basisname = "LCAO"; +} + +template +ESolver_DoubleXC::~ESolver_DoubleXC() +{ + delete this->psi_base; + delete this->p_hamilt_base; + delete this->pelec_base; +} + +template +void ESolver_DoubleXC::before_all_runners(UnitCell& ucell, const Input_para& inp) +{ + ModuleBase::TITLE("ESolver_DoubleXC", "before_all_runners"); + ModuleBase::timer::tick("ESolver_DoubleXC", "before_all_runners"); + + ESolver_KS_LCAO::before_all_runners(ucell, inp); + + // init some items for base functional + + // 2) init ElecState + if (this->pelec_base == nullptr) + { + this->pelec_base = new elecstate::ElecStateLCAO(&(this->chr_base), // use which parameter? + &(this->kv), + this->kv.get_nks(), + &(this->GG), + &(this->GK), + this->pw_rho, + this->pw_big); + } + + // 4) initialize electronic wave function psi + if (this->psi_base == nullptr) + { + int nsk = 0; + int ncol = 0; + if (PARAM.globalv.gamma_only_local) + { + nsk = PARAM.inp.nspin; + ncol = this->pv.ncol_bands; + if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "elpa" || PARAM.inp.ks_solver == "lapack" + || PARAM.inp.ks_solver == "pexsi" || PARAM.inp.ks_solver == "cusolver" + || PARAM.inp.ks_solver == "cusolvermp") + { + ncol = this->pv.ncol; + } + } + else + { + nsk = this->kv.get_nks(); +#ifdef __MPI + ncol = this->pv.ncol_bands; +#else + ncol = PARAM.inp.nbands; +#endif + } + this->psi_base = new psi::Psi(nsk, ncol, this->pv.nrow, this->kv.ngk, true); + } + + // 6) initialize the density matrix + dynamic_cast*>(this->pelec_base)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin); + + // 10) inititlize the charge density + this->chr_base.allocate(PARAM.inp.nspin); + this->pelec_base->omega = ucell.omega; + + // 11) initialize the potential + if (this->pelec_base->pot == nullptr) + { + this->pelec_base->pot = new elecstate::Potential(this->pw_rhod, + this->pw_rho, + &ucell, + &(this->locpp.vloc), + &(this->sf), + &(this->solvent), + &(this->pelec_base->f_en.etxc), + &(this->pelec_base->f_en.vtxc)); + } + + ModuleBase::timer::tick("ESolver_DoubleXC", "before_all_runners"); +} + +template +void ESolver_DoubleXC::before_scf(UnitCell& ucell, const int istep) +{ + ModuleBase::TITLE("ESolver_DoubleXC", "before_scf"); + ModuleBase::timer::tick("ESolver_DoubleXC", "before_scf"); + + ESolver_KS_LCAO::before_scf(ucell, istep); + + this->pelec_base->omega = ucell.omega; + //---------------------------------------------------------- + //! calculate D2 or D3 vdW + //---------------------------------------------------------- + auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running)); + if (vdw_solver != nullptr) + { + this->pelec_base->f_en.evdw = vdw_solver->get_energy(); + } + + //---------------------------------------------------------- + //! calculate ewald energy + //---------------------------------------------------------- + if (!PARAM.inp.test_skip_ewald) + { + //this->pelec_base->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac); + this->pelec_base->f_en.ewald_energy = this->pelec->f_en.ewald_energy; + } + + if (this->p_hamilt_base != nullptr) + { + delete this->p_hamilt_base; + this->p_hamilt_base = nullptr; + } + if (this->p_hamilt_base == nullptr) + { + elecstate::DensityMatrix* DM = dynamic_cast*>(this->pelec_base)->get_DM(); + + this->p_hamilt_base = new hamilt::HamiltLCAO( + PARAM.globalv.gamma_only_local ? &(this->GG) : nullptr, + PARAM.globalv.gamma_only_local ? nullptr : &(this->GK), + ucell, + this->gd, + &this->pv, + this->pelec_base->pot, + this->kv, + this->two_center_bundle_, + this->orb_, + DM +#ifdef __MLALGO + , + &this->ld +#endif +#ifdef __EXX + , + istep, + GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step, + GlobalC::exx_info.info_ri.real_number ? &this->exd->get_Hexxs() : nullptr, + GlobalC::exx_info.info_ri.real_number ? nullptr : &this->exc->get_Hexxs() +#endif + ); + } + + XC_Functional::set_xc_type(PARAM.inp.deepks_out_base); + this->pelec_base->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm); + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); + + // DMR should be same size with Hamiltonian(R) + dynamic_cast*>(this->pelec_base) + ->get_DM() + ->init_DMR(*(dynamic_cast*>(this->p_hamilt_base)->getHR())); + + if (istep > 0) + { + dynamic_cast*>(this->pelec_base)->get_DM()->cal_DMR(); + } + + ModuleBase::timer::tick("ESolver_DoubleXC", "before_scf"); + return; +} + +template +void ESolver_DoubleXC::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) +{ + ModuleBase::TITLE("ESolver_DoubleXC", "iter_finish"); + ModuleBase::timer::tick("ESolver_DoubleXC", "iter_finish"); + + bool output_iter = PARAM.inp.deepks_out_labels >0 && PARAM.inp.deepks_out_freq_elec && + (iter % PARAM.inp.deepks_out_freq_elec == 0); + + if ( output_iter ) + { + // save output charge density (density after diagnonalization) + for (int is = 0; is < PARAM.inp.nspin; is++) + { + ModuleBase::GlobalFunc::DCOPY(this->chr.rho[is], this->chr_base.rho[is], this->chr.rhopw->nrxx); + if (XC_Functional::get_ked_flag()) + { + ModuleBase::GlobalFunc::DCOPY(this->chr.kin_r[is], this->chr_base.kin_r[is], this->chr.rhopw->nrxx); + } + } + } + + ESolver_KS_LCAO::iter_finish(ucell, istep, iter, conv_esolver); + + // for deepks, output labels during electronic steps (after conv_esolver is renewed) + if ( output_iter) + { + // ---------- update etot and htot ---------- + // get etot of output charge density, now the etot is of density after charge mixing + this->pelec->pot->update_from_charge(&this->chr_base, &ucell); + this->pelec->f_en.descf = 0.0; + this->pelec->cal_energies(2); + // std::cout<<"in deepks etot------"<pelec->f_en.print_all(); + // std::cout<<"in deepks etot------"<pelec->f_en.etot << std::endl; + + // update p_hamilt using output charge density + // Note!!! + // This will change the result of out_mat_hs + // The original result of out_mat_hs is H of input density, but this change H to that of output density + // When converged, these two should be close + if (PARAM.inp.deepks_v_delta > 0 && PARAM.inp.vl_in_h) + { + // update real space Hamiltonian + this->p_hamilt->refresh(); + } + +#ifdef __MLALGO + // ---------- output tot and precalc ---------- + hamilt::HamiltLCAO* p_ham_deepks = dynamic_cast*>(this->p_hamilt); + std::shared_ptr> ld_shared_ptr(&this->ld, [](LCAO_Deepks*) {}); + LCAO_Deepks_Interface deepks_interface(ld_shared_ptr); + + deepks_interface.out_deepks_labels(this->pelec->f_en.etot, + this->kv.get_nks(), + ucell.nat, + PARAM.globalv.nlocal, + this->pelec->ekb, + this->kv.kvec_d, + ucell, + this->orb_, + this->gd, + &(this->pv), + *(this->psi), + dynamic_cast*>(this->pelec)->get_DM(), + p_ham_deepks, + iter, + conv_esolver, + GlobalV::MY_RANK, + GlobalV::ofs_running); +#endif + + // restore to density after charge mixing + this->pelec->pot->update_from_charge(&this->chr, &ucell); + + // ---------- prepare for base ---------- + // set as base functional Temporarily + XC_Functional::set_xc_type(PARAM.inp.deepks_out_base); + + // update pot of pelec_base according to chr_base + if (!conv_esolver) + { + this->pelec_base->pot->update_from_charge(&this->chr_base, &ucell); + } + else + { + this->pelec_base->cal_converged(); + } + + // ---------- e_base ---------- + // ebase use the same output density with etot, just different in xc + this->pelec_base->f_en.eband = this->pelec->f_en.eband; + this->pelec_base->f_en.deband = this->pelec->f_en.deband; + this->pelec_base->f_en.demet = this->pelec->f_en.demet; + this->pelec_base->f_en.descf = 0.0; // set descf to 0 + this->pelec_base->cal_energies(2); // 2 means Kohn-Sham functional + // std::cout<<"in double_xc------"<pelec_base->f_en.print_all(); + // std::cout<<"in double_xc------"<f_en.etot << std::endl; + +#ifdef __MLALGO + const std::string file_ebase = deepks_interface.get_filename("ebase", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_npy_e(pelec_base->f_en.etot, file_ebase, GlobalV::MY_RANK); +#endif + + // ---------- h_base ---------- + if (PARAM.inp.deepks_v_delta > 0) + { + if (PARAM.inp.vl_in_h) + { + // update real space Hamiltonian + this->p_hamilt_base->refresh(); + } + + // Note!!! + // should not use ModuleIO::write_hsk() to output h_base, because it will call get_hs_pointers() + // which will change the hsolver::DiagoElpa::DecomposedState, influencing the following SCF steps + +#ifdef __MLALGO + using TH = std::conditional_t::value, ModuleBase::matrix, ModuleBase::ComplexMatrix>; + hamilt::HamiltLCAO* p_ham_deepks_base = dynamic_cast*>(this->p_hamilt_base); + int nks = this->kv.get_nks(); + std::vector h_tot(nks); + DeePKS_domain::get_h_tot(this->pv, p_ham_deepks_base, h_tot, PARAM.globalv.nlocal, nks, 'H'); + + const std::string file_htot = deepks_interface.get_filename("hbase", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_npy_h(h_tot, file_htot, PARAM.globalv.nlocal, nks, GlobalV::MY_RANK); +#endif + } + + // ---------- o_base ---------- + if ( PARAM.inp.deepks_bandgap > 0 ) + { + // obase isn't implemented yet + // don't need to solve p_hamilt_base + // just dm*p_hamilt_base, similar to cal_o_delta + } + + // restore to original xc + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); + + } + // ---------- prepare for f_base ---------- + else if ( PARAM.inp.cal_force && conv_esolver ) + { + // vnew must be updated for force_scc() even if not output_iter + // set as base functional Temporarily + XC_Functional::set_xc_type(PARAM.inp.deepks_out_base); + this->pelec_base->cal_converged(); + // restore to original xc + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); + } + + if ( PARAM.inp.cal_force ) + { + if ( ! conv_esolver ) + { + // use chr after mixing to restore veff, useful for vnew when converged + this->pelec_base->pot->update_from_charge(&this->chr, &ucell); + } + else + { + // copy charge + for (int is = 0; is < PARAM.inp.nspin; is++) + { + ModuleBase::GlobalFunc::DCOPY(this->chr.rho[is], this->chr_base.rho[is], this->chr.rhopw->nrxx); + if (XC_Functional::get_ked_flag()) + { + ModuleBase::GlobalFunc::DCOPY(this->chr.kin_r[is], this->chr_base.kin_r[is], this->chr.rhopw->nrxx); + } + } + + // copy dm + int nks = this->kv.get_nks(); + auto _pes_lcao_base = dynamic_cast*>(this->pelec_base); + auto _pes_lcao = dynamic_cast*>(this->pelec); + for (int ik = 0; ik < nks; ik++) + { + _pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik)); + } + _pes_lcao_base->get_DM()->cal_DMR(); + _pes_lcao_base->ekb = _pes_lcao->ekb; + _pes_lcao_base->wg = _pes_lcao->wg; + } + } + ModuleBase::timer::tick("ESolver_DoubleXC", "iter_finish"); +} + +template +void ESolver_DoubleXC::cal_force(UnitCell& ucell, ModuleBase::matrix& force) +{ + ModuleBase::TITLE("ESolver_DoubleXC", "cal_force"); + ModuleBase::timer::tick("ESolver_DoubleXC", "cal_force"); + + ModuleBase::matrix force_base; + ModuleBase::matrix stress_base; + + Force_Stress_LCAO fsl(this->RA, ucell.nat); + + // set as base functional Temporarily + XC_Functional::set_xc_type(PARAM.inp.deepks_out_base); + + fsl.getForceStress(ucell, + PARAM.inp.cal_force, + PARAM.inp.cal_stress, + PARAM.inp.test_force, + PARAM.inp.test_stress, + this->gd, + this->pv, + this->pelec_base, + this->psi, + this->GG, // mohan add 2024-04-01 + this->GK, // mohan add 2024-04-01 + this->two_center_bundle_, + this->orb_, + force_base, + stress_base, + this->locpp, + this->sf, + this->kv, + this->pw_rho, + this->solvent, +#ifdef __MLALGO + this->ld, + "base", +#endif +#ifdef __EXX + *this->exd, + *this->exc, +#endif + &ucell.symm); + // restore to original xc + XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); + + // this will delete RA, so call it later + ESolver_KS_LCAO::cal_force(ucell, force); + + ModuleBase::timer::tick("ESolver_DoubleXC", "cal_force"); +} + +template class ESolver_DoubleXC; +template class ESolver_DoubleXC, double>; +template class ESolver_DoubleXC, std::complex>; + +} // namespace ModuleESolver diff --git a/source/source_esolver/esolver_double_xc.h b/source/source_esolver/esolver_double_xc.h new file mode 100644 index 0000000000..591411c82d --- /dev/null +++ b/source/source_esolver/esolver_double_xc.h @@ -0,0 +1,39 @@ +#ifndef ESOLVER_DOUBLE_XC_H +#define ESOLVER_DOUBLE_XC_H + +#include "source_esolver/esolver_ks_lcao.h" + +namespace ModuleESolver +{ +// used in deepks, run target and base xc functional simultaneously +template +class ESolver_DoubleXC : public ESolver_KS_LCAO +{ + public: + ESolver_DoubleXC(); + ~ESolver_DoubleXC(); + + void before_all_runners(UnitCell& ucell, const Input_para& inp) override; + + void cal_force(UnitCell& ucell, ModuleBase::matrix& force) override; + + protected: + + void before_scf(UnitCell& ucell, const int istep) override; + + void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override; + + //! Hamiltonian + hamilt::Hamilt* p_hamilt_base = nullptr; + + //! Electronic wavefunctions + psi::Psi* psi_base = nullptr; + + //! Electronic states + elecstate::ElecState* pelec_base = nullptr; + + //! Electorn charge density + Charge chr_base; +}; +} // namespace ModuleESolver +#endif diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 84844ab720..33191e5b38 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -361,6 +361,7 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for this->solvent, #ifdef __MLALGO this->ld, + "tot", #endif #ifdef __EXX *this->exd, diff --git a/source/source_hamilt/module_xc/xc_functional.cpp b/source/source_hamilt/module_xc/xc_functional.cpp index e530cc9cdd..2462f8ac42 100644 --- a/source/source_hamilt/module_xc/xc_functional.cpp +++ b/source/source_hamilt/module_xc/xc_functional.cpp @@ -278,6 +278,10 @@ void XC_Functional::set_xc_type(const std::string xc_func_in) { ked_flag = true; } + else + { + ked_flag = false; + } if (func_id[0] == XC_GGA_X_OPTX) { diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 2e79e806a2..fc18a80ec9 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -265,6 +265,8 @@ struct Input_para ///< descriptors for training, wenfei 2022-1-12 int deepks_out_freq_elec = 0; ///< (need libnpy) frequency of electronic iteration to output ///< descriptors and labels, default is 0, which means no output until convergence + std::string deepks_out_base = "none"; ///< (need libnpy) base functional for output files, with dft_functional as target functional + /// default is "none", which means no base functional bool deepks_scf = false; ///< (need libnpy and libtorch) if set to true, a trained model ///< would be needed to calculate V_delta and F_delta int deepks_bandgap = 0; ///< for bandgap label. QO added 2021-12-15 diff --git a/source/source_io/read_input_item_deepks.cpp b/source/source_io/read_input_item_deepks.cpp index 15356ab0cf..9f51d3511a 100644 --- a/source/source_io/read_input_item_deepks.cpp +++ b/source/source_io/read_input_item_deepks.cpp @@ -23,6 +23,28 @@ void ReadInput::item_deepks() para.input.deepks_out_freq_elec = 0; } }; + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.deepks_out_freq_elec > 0 && para.input.deepks_out_base == "none") + { + ModuleBase::WARNING_QUIT("ReadInput", "to use deepks_out_freq_elec, please set deepks_out_base "); + } + }; + this->add_item(item); + } + { + Input_Item item("deepks_out_base"); + item.annotation = "base functional for output files, with dft_functional as target functional"; + read_sync_string(input.deepks_out_base); + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.deepks_out_base != "none" && para.input.deepks_out_labels == 0) + { + ModuleBase::WARNING_QUIT("ReadInput", "to use deepks_out_base, please set deepks_out_labels > 0 "); + } + if (para.input.deepks_out_base != "none" && para.input.deepks_bandgap > 0 ) + { + ModuleBase::WARNING_QUIT("ReadInput", "outputting bandgap labels during electronic steps is not implemented yet "); + } + }; this->add_item(item); } { diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp index cd6230b061..03d91e319d 100644 --- a/source/source_lcao/FORCE_STRESS.cpp +++ b/source/source_lcao/FORCE_STRESS.cpp @@ -53,6 +53,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, surchem& solvent, #ifdef __MLALGO LCAO_Deepks& ld, + const std::string& dpks_out_type, #endif #ifdef __EXX Exx_LRI_Interface& exd, @@ -497,24 +498,37 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // DeePKS force if (PARAM.inp.deepks_out_labels) // not parallelized yet { - const std::string file_ftot = PARAM.globalv.global_out_dir - + (PARAM.inp.deepks_out_labels == 1 ? "deepks_ftot.npy" : "deepks_force.npy"); - LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Hartree/Bohr, F_tot - - if (PARAM.inp.deepks_out_labels == 1) + if (PARAM.inp.deepks_out_base == "none" + || (PARAM.inp.deepks_out_base != "none" && dpks_out_type == "tot") ) { - const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; - if (PARAM.inp.deepks_scf) - { - LCAO_deepks_io::save_matrix2npy(file_fbase, - fcs - fvnl_dalpha, - GlobalV::MY_RANK); // Hartree/Bohr, F_base - } - else + const std::string file_ftot = PARAM.globalv.global_out_dir + + (PARAM.inp.deepks_out_labels == 1 ? "deepks_ftot.npy" : "deepks_force.npy"); + LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Hartree/Bohr, F_tot + + if (PARAM.inp.deepks_out_labels == 1) { - LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot + // this base only considers subtracting the deepks_scf part + const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; + if (PARAM.inp.deepks_scf) + { + LCAO_deepks_io::save_matrix2npy(file_fbase, + fcs - fvnl_dalpha, + GlobalV::MY_RANK); // Hartree/Bohr, F_base + } + else + { + LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot + } } } + if (PARAM.inp.deepks_out_base != "none") + { + // output fcs as tot or base in another dir + // this base considers changing xc functional to base functional + const std::string file_f = PARAM.globalv.global_deepks_label_elec_dir + + (dpks_out_type == "tot" ? "ftot.npy" : "fbase.npy"); + LCAO_deepks_io::save_matrix2npy(file_f, fcs, GlobalV::MY_RANK); + } } #endif // print Rydberg force or not @@ -685,29 +699,46 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, #ifdef __MLALGO if (PARAM.inp.deepks_out_labels == 1) { - const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; - LCAO_deepks_io::save_matrix2npy(file_stot, - scs, - GlobalV::MY_RANK, - ucell.omega, - 'U'); // change to energy unit Ry when printing, S_tot; - - const std::string file_sbase = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; - if (PARAM.inp.deepks_scf) + if (PARAM.inp.deepks_out_base == "none" + || (PARAM.inp.deepks_out_base != "none" && dpks_out_type == "tot") ) { - LCAO_deepks_io::save_matrix2npy(file_sbase, - scs - svnl_dalpha, + const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; + LCAO_deepks_io::save_matrix2npy(file_stot, + scs, GlobalV::MY_RANK, ucell.omega, - 'U'); // change to energy unit Ry when printing, S_base; + 'U'); // change to energy unit Ry when printing, S_tot; + + // this base only considers subtracting the deepks_scf part + const std::string file_sbase = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; + if (PARAM.inp.deepks_scf) + { + LCAO_deepks_io::save_matrix2npy(file_sbase, + scs - svnl_dalpha, + GlobalV::MY_RANK, + ucell.omega, + 'U'); // change to energy unit Ry when printing, S_base; + } + else + { + LCAO_deepks_io::save_matrix2npy(file_sbase, + scs, + GlobalV::MY_RANK, + ucell.omega, + 'U'); // sbase = stot + } } - else + if (PARAM.inp.deepks_out_base != "none") { - LCAO_deepks_io::save_matrix2npy(file_sbase, + // output scs as tot or base in another dir + // this base considers changing xc functional to base functional + const std::string file_s = PARAM.globalv.global_deepks_label_elec_dir + + (dpks_out_type == "tot" ? "stot.npy" : "sbase.npy"); + LCAO_deepks_io::save_matrix2npy(file_s, scs, GlobalV::MY_RANK, ucell.omega, - 'U'); // sbase = stot + 'U'); // change to energy unit Ry when printing, S_tot; } } else if (PARAM.inp.deepks_out_labels == 2) diff --git a/source/source_lcao/FORCE_STRESS.h b/source/source_lcao/FORCE_STRESS.h index d2400d7b69..3c6eb65745 100644 --- a/source/source_lcao/FORCE_STRESS.h +++ b/source/source_lcao/FORCE_STRESS.h @@ -51,6 +51,7 @@ class Force_Stress_LCAO surchem& solvent, #ifdef __MLALGO LCAO_Deepks& ld, + const std::string& dpks_out_type, #endif #ifdef __EXX Exx_LRI_Interface& exd, diff --git a/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp b/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp index d2bfe71d14..04eb0cffc2 100644 --- a/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp +++ b/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp @@ -33,9 +33,10 @@ std::string true_file_type(const std::string& file_type) // global_out_dir/deepks_*.npy for iter=-1 (called in after_scf) // global_out_dir/DeePKS_Labels_Elec/*_e*.npy for iter>0 (called during electronic steps) -std::string get_filename(const std::string& file_type, - const int& label_type, - const int& iter) +template +std::string LCAO_Deepks_Interface::get_filename(const std::string& file_type, + const int& label_type, + const int& iter) { std::ostringstream file_name; file_name << (iter == -1 ? PARAM.globalv.global_out_dir : PARAM.globalv.global_deepks_label_elec_dir); @@ -95,13 +96,20 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, double E_delta = ld->E_delta; double e_delta_band = ld->e_delta_band; hamilt::HContainer* dmr = ld->dm_r; + // Used for deepks_bandgap == 1 and deepks_v_delta > 0 + std::vector>* h_delta = &ld->V_delta; const int nspin = PARAM.inp.nspin; const int nk = nks / nspin; - const bool not_first_step = (iter != 1); // not output in the first electronic step, for energy and otot/obase - const bool not_last_step = (iter == -1) || !conv_esolver; //not output in the last electronic step const bool is_after_scf = (iter == -1); // called in after_scf, not in electronic steps + const bool output_base = PARAM.inp.deepks_out_labels == 1 && is_after_scf; // not output when deepks_out_labels=2 and in electronic step (output true base elsewhere) + const bool output_precalc = (PARAM.inp.deepks_out_labels == 1) && (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_freq_elec); + +//================================================================================ +// 1. Update real-space density matrix (DMR) for deepks, projected density matrix (PDM) +// and descriptor. Output descriptor if needed. +//================================================================================ // Update DMR in any case of deepks_out_labels/deepks_scf DeePKS_domain::update_dmr(kvec_d, dm->get_DMK_vector(), ucell, orb, *ParaV, GridD, dmr); @@ -128,19 +136,15 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, descriptor, rank); - if ( not_last_step ) - { - const int true_iter = is_after_scf ? iter : iter + 1; - const std::string file_d = get_filename("dm_eig", PARAM.inp.deepks_out_labels, true_iter); - LCAO_deepks_io::save_npy_d(nat, - des_per_atom, - inlmax, - inl2l, - PARAM.inp.deepks_equiv, - descriptor, - file_d, - rank); // libnpy needed - } + const std::string file_d = get_filename("dm_eig", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_npy_d(nat, + des_per_atom, + inlmax, + inl2l, + PARAM.inp.deepks_equiv, + descriptor, + file_d, + rank); // libnpy needed if (PARAM.inp.deepks_scf) @@ -175,42 +179,99 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, } } - // Used for deepks_bandgap == 1 and deepks_v_delta > 0 - std::vector>* h_delta = &ld->V_delta; - // calculating deepks correction and save the results if (PARAM.inp.deepks_out_labels) { // Used for deepks_scf == 1 or deepks_out_freq_elec!=0, for *precalc items, not for deepks_out_labels=2 std::vector gevdm; - if ((PARAM.inp.deepks_scf || PARAM.inp.deepks_out_freq_elec) && PARAM.inp.deepks_out_labels !=2 ) + if ( output_precalc ) { DeePKS_domain::cal_gevdm(nat, inlmax, inl2l, pdm, gevdm); } - if ( not_first_step) +//================================================================================ +// 2. Energy +//================================================================================ + + // etot + const std::string file_etot = get_filename("etot", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_npy_e(etot, file_etot, rank); + + // ebase + if ( output_base ) { - // Energy Part - const std::string file_etot = get_filename("etot", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_npy_e(etot, file_etot, rank); + const std::string file_ebase = get_filename("ebase", PARAM.inp.deepks_out_labels, iter); + if (PARAM.inp.deepks_scf) + { + /// ebase :no deepks E_delta including + LCAO_deepks_io::save_npy_e(etot - E_delta, file_ebase, rank); + } + else // deepks_scf = 0; base calculation + { + /// no scf, e_tot=e_base + LCAO_deepks_io::save_npy_e(etot, file_ebase, rank); + } + } - if (PARAM.inp.deepks_out_labels == 1) +//================================================================================ +// 3. Force and Stress +//================================================================================ + + if ( is_after_scf ) + { + // Force Part + if (PARAM.inp.cal_force) { - const std::string file_ebase = get_filename("ebase", PARAM.inp.deepks_out_labels, iter); - if (PARAM.inp.deepks_scf) + // these items are not related to model, so can output without deepks_scf + if (output_precalc // don't need these when deepks_out_labels == 2 + && !PARAM.inp.deepks_equiv) // training with force label not supported by equivariant version now { - /// ebase :no deepks E_delta including - LCAO_deepks_io::save_npy_e(etot - E_delta, file_ebase, rank); + torch::Tensor gdmx; + DeePKS_domain::cal_gdmx< + TK>(lmaxd, inlmax, nks, kvec_d, phialpha, inl_index, dmr, ucell, orb, *ParaV, GridD, gdmx); + + torch::Tensor gvx; + DeePKS_domain::cal_gvx(ucell.nat, inlmax, des_per_atom, inl2l, gevdm, gdmx, gvx, rank); + const std::string file_gradvx = get_filename("gradvx", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_tensor2npy(file_gradvx, gvx, rank); + + if (PARAM.inp.deepks_out_unittest) + { + DeePKS_domain::check_tensor(gdmx, "gdmx.dat", rank); + DeePKS_domain::check_tensor(gvx, "gvx.dat", rank); + } } - else // deepks_scf = 0; base calculation + } + + // Stress Part + if (PARAM.inp.cal_stress) + { + // these items are not related to model, so can output without deepks_scf + if (output_precalc // don't need these when deepks_out_labels == 2 + && !PARAM.inp.deepks_equiv) // training with stress label not supported by equivariant version now { - /// no scf, e_tot=e_base - LCAO_deepks_io::save_npy_e(etot, file_ebase, rank); + torch::Tensor gdmepsl; + DeePKS_domain::cal_gdmepsl< + TK>(lmaxd, inlmax, nks, kvec_d, phialpha, inl_index, dmr, ucell, orb, *ParaV, GridD, gdmepsl); + + torch::Tensor gvepsl; + DeePKS_domain::cal_gvepsl(ucell.nat, inlmax, des_per_atom, inl2l, gevdm, gdmepsl, gvepsl, rank); + const std::string file_gvepsl = get_filename("gvepsl", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_tensor2npy(file_gvepsl, gvepsl, rank); + + if (PARAM.inp.deepks_out_unittest) + { + DeePKS_domain::check_tensor(gdmepsl, "gdmepsl.dat", rank); + DeePKS_domain::check_tensor(gvepsl, "gvepsl.dat", rank); + } } - } + } } - // Bandgap Part +//================================================================================ +// 4. Bandgap +//================================================================================ + if (PARAM.inp.deepks_bandgap > 0) { // Get the number of the occupied bands @@ -244,188 +305,129 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, // Calculate the bandgap for each k point ModuleBase::matrix o_tot(nks, range); - if ( not_first_step) + for (int iks = 0; iks < nks; ++iks) { - for (int iks = 0; iks < nks; ++iks) + int ib = 0; + if (PARAM.inp.deepks_bandgap == 1 || PARAM.inp.deepks_bandgap == 3) { - int ib = 0; - if (PARAM.inp.deepks_bandgap == 1 || PARAM.inp.deepks_bandgap == 3) - { - o_tot(iks, ib) = ekb(iks, nocc + PARAM.inp.deepks_band_range[1]) - - ekb(iks, nocc + PARAM.inp.deepks_band_range[0]); - } - else if (PARAM.inp.deepks_bandgap == 2) + o_tot(iks, ib) = ekb(iks, nocc + PARAM.inp.deepks_band_range[1]) + - ekb(iks, nocc + PARAM.inp.deepks_band_range[0]); + } + else if (PARAM.inp.deepks_bandgap == 2) + { + for (int ir = PARAM.inp.deepks_band_range[0]; ir <= PARAM.inp.deepks_band_range[1]; ++ir) { - for (int ir = PARAM.inp.deepks_band_range[0]; ir <= PARAM.inp.deepks_band_range[1]; ++ir) + if (ir != -1) { - if (ir != -1) - { - o_tot(iks, ib) = ekb(iks, nocc + ir) - ekb(iks, nocc - 1); - ib++; - } + o_tot(iks, ib) = ekb(iks, nocc + ir) - ekb(iks, nocc - 1); + ib++; } - assert(ib == range); // ensure that we have filled all the bandgap values } + assert(ib == range); // ensure that we have filled all the bandgap values } - - const std::string file_otot = get_filename("otot", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_matrix2npy(file_otot, o_tot, rank); // Unit: Hartree } + const std::string file_otot = get_filename("otot", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_matrix2npy(file_otot, o_tot, rank); // Unit: Hartree - if (PARAM.inp.deepks_out_labels == 1) // don't need these when deepks_out_labels == 2 + // don't need these when deepks_out_labels == 2 + // not consider out_base now, because bandgap is not supported in deepks_out_freq_elec now + if ( output_precalc ) { - if (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_freq_elec) + std::vector wg_hl_range(range); + for (int ir = 0; ir < range; ++ir) + { + wg_hl_range[ir].create(nks, PARAM.inp.nbands); + wg_hl_range[ir].zero_out(); + } + + // Calculate O_delta + for (int iks = 0; iks < nks; ++iks) { - std::vector wg_hl_range(range); - for (int ir = 0; ir < range; ++ir) + int ib = 0; + if (PARAM.inp.deepks_bandgap == 1 || PARAM.inp.deepks_bandgap == 3) { - wg_hl_range[ir].create(nks, PARAM.inp.nbands); - wg_hl_range[ir].zero_out(); + wg_hl_range[ib](iks, nocc + PARAM.inp.deepks_band_range[0]) = -1.0; + wg_hl_range[ib](iks, nocc + PARAM.inp.deepks_band_range[1]) = 1.0; } - - // Calculate O_delta - for (int iks = 0; iks < nks; ++iks) + else if (PARAM.inp.deepks_bandgap == 2) { - int ib = 0; - if (PARAM.inp.deepks_bandgap == 1 || PARAM.inp.deepks_bandgap == 3) - { - wg_hl_range[ib](iks, nocc + PARAM.inp.deepks_band_range[0]) = -1.0; - wg_hl_range[ib](iks, nocc + PARAM.inp.deepks_band_range[1]) = 1.0; - } - else if (PARAM.inp.deepks_bandgap == 2) + for (int ir = PARAM.inp.deepks_band_range[0]; ir <= PARAM.inp.deepks_band_range[1]; ++ir) { - for (int ir = PARAM.inp.deepks_band_range[0]; ir <= PARAM.inp.deepks_band_range[1]; ++ir) + if (ir != -1) { - if (ir != -1) - { - wg_hl_range[ib](iks, nocc - 1) = -1.0; - wg_hl_range[ib](iks, nocc + ir) = 1.0; - ib++; - } + wg_hl_range[ib](iks, nocc - 1) = -1.0; + wg_hl_range[ib](iks, nocc + ir) = 1.0; + ib++; } } } + } - ModuleBase::matrix o_delta(nks, range); - torch::Tensor orbital_precalc; - for (int ir = 0; ir < range; ++ir) - { - std::vector dm_bandgap(nks); - elecstate::cal_dm(ParaV, wg_hl_range[ir], psi, dm_bandgap); - - torch::Tensor orbital_precalc_temp; - ModuleBase::matrix o_delta_temp(nks, 1); - DeePKS_domain::cal_orbital_precalc(dm_bandgap, - lmaxd, - inlmax, - nat, - nks, - inl2l, - kvec_d, - phialpha, - gevdm, - inl_index, - ucell, - orb, - *ParaV, - GridD, - orbital_precalc_temp); - if (ir == 0) - { - orbital_precalc = orbital_precalc_temp; - } - else - { - orbital_precalc = torch::cat({orbital_precalc, orbital_precalc_temp}, 0); - } + ModuleBase::matrix o_delta(nks, range); + torch::Tensor orbital_precalc; + for (int ir = 0; ir < range; ++ir) + { + std::vector dm_bandgap(nks); + elecstate::cal_dm(ParaV, wg_hl_range[ir], psi, dm_bandgap); - if (PARAM.inp.deepks_scf) - { - DeePKS_domain::cal_o_delta(dm_bandgap, *h_delta, o_delta_temp, *ParaV, nks, nspin); - for (int iks = 0; iks < nks; ++iks) - { - o_delta(iks, ir) = o_delta_temp(iks, 0); - } - } - } - // save obase and orbital_precalc - if ( not_last_step ) + torch::Tensor orbital_precalc_temp; + ModuleBase::matrix o_delta_temp(nks, 1); + DeePKS_domain::cal_orbital_precalc(dm_bandgap, + lmaxd, + inlmax, + nat, + nks, + inl2l, + kvec_d, + phialpha, + gevdm, + inl_index, + ucell, + orb, + *ParaV, + GridD, + orbital_precalc_temp); + if (ir == 0) { - const int true_iter = (iter == -1) ? iter : iter + 1; - const std::string file_orbpre = get_filename("orbpre", PARAM.inp.deepks_out_labels, true_iter); - LCAO_deepks_io::save_tensor2npy(file_orbpre, orbital_precalc, rank); + orbital_precalc = orbital_precalc_temp; } - - if ( not_first_step) + else { - if (PARAM.inp.deepks_scf) - { - const std::string file_obase = get_filename("obase", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_matrix2npy(file_obase, o_tot - o_delta, rank); // Unit: Hartree - } - else - { - const std::string file_obase = get_filename("obase", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_matrix2npy(file_obase, o_tot, rank); // no scf, o_tot=o_base - } + orbital_precalc = torch::cat({orbital_precalc, orbital_precalc_temp}, 0); } - } - } // end deepks_out_labels == 1 - } // end deepks_bandgap > 0 - - if ( is_after_scf ) - { - // Force Part - if (PARAM.inp.cal_force) - { - // these items are not related to model, so can output without deepks_scf - if (PARAM.inp.deepks_out_labels == 1 // don't need these when deepks_out_labels == 2 - && !PARAM.inp.deepks_equiv) // training with force label not supported by equivariant version now - { - torch::Tensor gdmx; - DeePKS_domain::cal_gdmx< - TK>(lmaxd, inlmax, nks, kvec_d, phialpha, inl_index, dmr, ucell, orb, *ParaV, GridD, gdmx); - - torch::Tensor gvx; - DeePKS_domain::cal_gvx(ucell.nat, inlmax, des_per_atom, inl2l, gevdm, gdmx, gvx, rank); - const std::string file_gradvx = get_filename("gradvx", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_tensor2npy(file_gradvx, gvx, rank); - if (PARAM.inp.deepks_out_unittest) + if (PARAM.inp.deepks_scf) { - DeePKS_domain::check_tensor(gdmx, "gdmx.dat", rank); - DeePKS_domain::check_tensor(gvx, "gvx.dat", rank); + DeePKS_domain::cal_o_delta(dm_bandgap, *h_delta, o_delta_temp, *ParaV, nks, nspin); + for (int iks = 0; iks < nks; ++iks) + { + o_delta(iks, ir) = o_delta_temp(iks, 0); + } } } - } + // save obase and orbital_precalc + const std::string file_orbpre = get_filename("orbpre", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_tensor2npy(file_orbpre, orbital_precalc, rank); - // Stress Part - if (PARAM.inp.cal_stress) - { - // these items are not related to model, so can output without deepks_scf - if (PARAM.inp.deepks_out_labels == 1 // don't need these when deepks_out_labels == 2 - && !PARAM.inp.deepks_equiv) // training with stress label not supported by equivariant version now + if (PARAM.inp.deepks_scf) { - torch::Tensor gdmepsl; - DeePKS_domain::cal_gdmepsl< - TK>(lmaxd, inlmax, nks, kvec_d, phialpha, inl_index, dmr, ucell, orb, *ParaV, GridD, gdmepsl); - - torch::Tensor gvepsl; - DeePKS_domain::cal_gvepsl(ucell.nat, inlmax, des_per_atom, inl2l, gevdm, gdmepsl, gvepsl, rank); - const std::string file_gvepsl = get_filename("gvepsl", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_tensor2npy(file_gvepsl, gvepsl, rank); - - if (PARAM.inp.deepks_out_unittest) - { - DeePKS_domain::check_tensor(gdmepsl, "gdmepsl.dat", rank); - DeePKS_domain::check_tensor(gvepsl, "gvepsl.dat", rank); - } + const std::string file_obase = get_filename("obase", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_matrix2npy(file_obase, o_tot - o_delta, rank); // Unit: Hartree } - } - - + else + { + const std::string file_obase = get_filename("obase", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_matrix2npy(file_obase, o_tot, rank); // no scf, o_tot=o_base + } + } // end out_precalc + } // end deepks_bandgap > 0 +//================================================================================ +// 5. HR +//================================================================================ + if ( is_after_scf ) + { // not add deepks_out_labels = 2 and deepks_out_freq_elec for HR yet // H(R) matrix part, for HR, base will not be calculated since they are HContainer objects if (PARAM.inp.deepks_v_delta < 0) @@ -517,105 +519,107 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, } } - if ( not_last_step ) +//================================================================================ +// 6. Hk +//================================================================================ + + if (PARAM.inp.deepks_v_delta > 0) { - const int true_iter = is_after_scf ? iter : iter + 1; - // H(k) matrix part - if (PARAM.inp.deepks_v_delta > 0) - { - std::vector h_tot(nks); - DeePKS_domain::get_h_tot(*ParaV, p_ham, h_tot, nlocal, nks, 'H'); + std::vector h_tot(nks); + DeePKS_domain::get_h_tot(*ParaV, p_ham, h_tot, nlocal, nks, 'H'); - const std::string file_htot = get_filename("htot", PARAM.inp.deepks_out_labels, true_iter); - LCAO_deepks_io::save_npy_h(h_tot, file_htot, nlocal, nks, rank); + const std::string file_htot = get_filename("htot", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_npy_h(h_tot, file_htot, nlocal, nks, rank); - if (PARAM.inp.deepks_out_labels == 1) // don't need these when deepks_out_labels == 2 + if ( output_base ) + { + if (PARAM.inp.deepks_scf) { - if (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_freq_elec) + std::vector v_delta(nks); + std::vector h_base(nks); + for (int ik = 0; ik < nks; ik++) { - if (PARAM.inp.deepks_scf) - { - std::vector v_delta(nks); - std::vector h_base(nks); - for (int ik = 0; ik < nks; ik++) - { - v_delta[ik].create(nlocal, nlocal); - h_base[ik].create(nlocal, nlocal); - } - DeePKS_domain::collect_h_mat(*ParaV, *h_delta, v_delta, nlocal, nks); - - // save v_delta and h_base - const std::string file_hbase = get_filename("hbase", PARAM.inp.deepks_out_labels, true_iter); - for (int ik = 0; ik < nks; ik++) - { - h_base[ik] = h_tot[ik] - v_delta[ik]; - } - LCAO_deepks_io::save_npy_h(h_base, file_hbase, nlocal, nks, rank); + v_delta[ik].create(nlocal, nlocal); + h_base[ik].create(nlocal, nlocal); + } + DeePKS_domain::collect_h_mat(*ParaV, *h_delta, v_delta, nlocal, nks); - const std::string file_vdelta = get_filename("vdelta", PARAM.inp.deepks_out_labels, true_iter); - LCAO_deepks_io::save_npy_h(v_delta, file_vdelta, nlocal, nks, rank); - } - else // deepks_scf == 0 - { - const std::string file_hbase = get_filename("hbase", PARAM.inp.deepks_out_labels, true_iter); - LCAO_deepks_io::save_npy_h(h_tot, file_hbase, nlocal, nks, rank); - } + // save v_delta and h_base + const std::string file_hbase = get_filename("hbase", PARAM.inp.deepks_out_labels, iter); + for (int ik = 0; ik < nks; ik++) + { + h_base[ik] = h_tot[ik] - v_delta[ik]; + } + LCAO_deepks_io::save_npy_h(h_base, file_hbase, nlocal, nks, rank); - if (PARAM.inp.deepks_v_delta == 1) // v_delta_precalc storage method 1 - { - torch::Tensor v_delta_precalc; - DeePKS_domain::cal_v_delta_precalc(nlocal, - lmaxd, - inlmax, - nat, - nks, - inl2l, - kvec_d, - phialpha, - gevdm, - inl_index, - ucell, - orb, - *ParaV, - GridD, - v_delta_precalc); + const std::string file_vdelta = get_filename("vdelta", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_npy_h(v_delta, file_vdelta, nlocal, nks, rank); + } + else // deepks_scf == 0 + { + const std::string file_hbase = get_filename("hbase", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_npy_h(h_tot, file_hbase, nlocal, nks, rank); + } + } - const std::string file_vdpre = get_filename("vdpre", PARAM.inp.deepks_out_labels, true_iter); - LCAO_deepks_io::save_tensor2npy(file_vdpre, v_delta_precalc, rank); - } - else if (PARAM.inp.deepks_v_delta == 2) // v_delta_precalc storage method 2 - { - torch::Tensor phialpha_out; - DeePKS_domain::prepare_phialpha(nlocal, - lmaxd, - inlmax, - nat, - nks, - kvec_d, - phialpha, - ucell, - orb, - *ParaV, - GridD, - phialpha_out); - const std::string file_phialpha = get_filename("phialpha", PARAM.inp.deepks_out_labels, true_iter); - LCAO_deepks_io::save_tensor2npy(file_phialpha, phialpha_out, rank); + if ( output_precalc ) + { + if (PARAM.inp.deepks_v_delta == 1) // v_delta_precalc storage method 1 + { + torch::Tensor v_delta_precalc; + DeePKS_domain::cal_v_delta_precalc(nlocal, + lmaxd, + inlmax, + nat, + nks, + inl2l, + kvec_d, + phialpha, + gevdm, + inl_index, + ucell, + orb, + *ParaV, + GridD, + v_delta_precalc); - torch::Tensor gevdm_out; - DeePKS_domain::prepare_gevdm(nat, lmaxd, inlmax, orb, gevdm, gevdm_out); - const std::string file_gevdm = get_filename("gevdm", PARAM.inp.deepks_out_labels, true_iter); - LCAO_deepks_io::save_tensor2npy(file_gevdm, gevdm_out, rank); - } - } - } // end deepks_out_labels == 1 - } // end v_delta label - } + const std::string file_vdpre = get_filename("vdpre", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_tensor2npy(file_vdpre, v_delta_precalc, rank); + } + else if (PARAM.inp.deepks_v_delta == 2) // v_delta_precalc storage method 2 + { + torch::Tensor phialpha_out; + DeePKS_domain::prepare_phialpha(nlocal, + lmaxd, + inlmax, + nat, + nks, + kvec_d, + phialpha, + ucell, + orb, + *ParaV, + GridD, + phialpha_out); + const std::string file_phialpha = get_filename("phialpha", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_tensor2npy(file_phialpha, phialpha_out, rank); + + torch::Tensor gevdm_out; + DeePKS_domain::prepare_gevdm(nat, lmaxd, inlmax, orb, gevdm, gevdm_out); + const std::string file_gevdm = get_filename("gevdm", PARAM.inp.deepks_out_labels, iter); + LCAO_deepks_io::save_tensor2npy(file_gevdm, gevdm_out, rank); + } + } + } // end v_delta label } // end deepks_out_labels - if (iter < 0)// only output when called in after_scf +//================================================================================ +// 7. atom.npy, box.npy, overlap.npy +//================================================================================ + + if ( is_after_scf ) // don't need to output in multiple electronic steps { - // don't need to output in multiple electronic steps if (PARAM.inp.deepks_out_labels == 2) { // output atom.npy and box.npy @@ -644,6 +648,9 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, } } +//================================================================================ +// 8. print, unittest +//================================================================================ /// print out deepks information to the screen if (PARAM.inp.deepks_scf) { diff --git a/source/source_lcao/module_deepks/LCAO_deepks_interface.h b/source/source_lcao/module_deepks/LCAO_deepks_interface.h index dab32855fb..a637e8839d 100644 --- a/source/source_lcao/module_deepks/LCAO_deepks_interface.h +++ b/source/source_lcao/module_deepks/LCAO_deepks_interface.h @@ -50,6 +50,15 @@ class LCAO_Deepks_Interface const bool& conv_esolver, const int rank, std::ostream& ofs_running); + + /// @brief Get the filename for deepks output files + /// @param file_type Type of the file (e.g., "etot", "ftot", etc.) + /// @param label_type Type of the label (from PARAM.inp.deepks_out_labels) + /// @param iter Iteration number (e.g., -1 for after_scf, or specific iteration number) + /// @return The full path to the output file + std::string get_filename(const std::string& file_type, + const int& label_type, + const int& iter); private: std::shared_ptr> ld; diff --git a/source/source_lcao/module_operator_lcao/operator_lcao.cpp b/source/source_lcao/module_operator_lcao/operator_lcao.cpp index a7e47450fe..785b214e08 100644 --- a/source/source_lcao/module_operator_lcao/operator_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/operator_lcao.cpp @@ -31,6 +31,8 @@ void OperatorLCAO::get_hs_pointers() { const int inc = 1; BlasConnector::copy(this->hsk->get_size(), this->hsk->get_sk(), inc, this->smatrix_k, inc); #ifdef __ELPA + // DecomposedState may be changed after diagnolization, with smatrix_k changed too + // for example, when DecomposedState equals 1, smatrix_k is an identity matrix, different from the original this->hsk->get_sk() hsolver::DiagoElpa::DecomposedState = 0; hsolver::DiagoElpaNative::DecomposedState = 0; #endif diff --git a/tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/INPUT b/tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/INPUT index dec208abeb..5a6f3b87f4 100644 --- a/tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/INPUT +++ b/tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/INPUT @@ -29,11 +29,11 @@ mixing_gg0 0 #Parameters (6.File) deepks_out_labels 1 -deepks_bandgap 1 deepks_v_delta 1 cal_force 1 cal_stress 1 deepks_out_freq_elec 1 -dft_functional lda +dft_functional pbe +deepks_out_base lda diff --git a/tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/result.ref b/tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/result.ref index 58d0d417a1..5c61657bb1 100644 --- a/tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/result.ref +++ b/tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/result.ref @@ -1,27 +1,25 @@ -etotref -463.3717558791495890 -etotperatomref -154.4572519597 -totalforceref 3.661556 -totalstressref 9.690490 -deepks_e_label 17.02859183994638 +etotref -466.0484517111885907 +etotperatomref -155.3494839037 +totalforceref 3.432943 +totalstressref 7.674750 +deepks_e_label 17.126958562184335 deepks_edelta 0.0 -deepks_o_label 0.3103148500892656 -deepks_odelta 0.0 -deepks_oprec -0.42236872134755643 -deepks_h_label 49.03922791708764 +deepks_h_label 49.12707696738566 deepks_vdelta 0.0 -deepks_vdp 176.26236130240005 -deepks_f_label 0.07120588316979731 +deepks_vdp 176.08561725742982 +deepks_f_label 0.06676003678176294 deepks_fdelta 0.0 -deepks_fpre 19.53717375096851 -deepks_s_label 0.08886477565384877 +deepks_s_label 0.06851710382835817 deepks_sdelta 0.0 -deepks_spre 19.22629462223564 -deepks_e_label_elec 51.047596216236645 -deepks_edelta_elec 0 -deepks_o_label_elec .86018933752816427 -deepks_odelta_elec 0 -deepks_oprec_elec -1.23599151361440239 -deepks_h_label_elec 147.028760179109205 -deepks_vdelta_elec 0 -deepks_vdp_elec 529.07335498441735 -totaltimeref 1.03 +deepks_fpre 19.69171134871589 +deepks_spre 19.372132511297153 +deepks_e_label_elec 85.389836665688162 +deepks_edelta_elec .49279432475411510 +deepks_h_label_elec 246.793639290695066 +deepks_vdelta_elec 3.4234699193308799 +deepks_vdp_elec 881.06539843219689 +deepks_f_label_elec .06676003678176294 +deepks_fdelta_elec .017979102769780632 +deepks_s_label_elec .06851710382835817 +deepks_sdelta_elec .013725178909672755 +totaltimeref 3.76 diff --git a/tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/INPUT b/tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/INPUT index 9a82dfcb10..f2687c0e2b 100644 --- a/tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/INPUT +++ b/tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/INPUT @@ -29,11 +29,11 @@ mixing_gg0 0 #Parameters (6.File) deepks_out_labels 1 -deepks_bandgap 1 deepks_v_delta 1 cal_force 1 cal_stress 1 deepks_out_freq_elec 1 -dft_functional lda +dft_functional pbe +deepks_out_base lda diff --git a/tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/result.ref b/tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/result.ref index 56a29b13c8..bc76b36bcc 100644 --- a/tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/result.ref +++ b/tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/result.ref @@ -1,27 +1,25 @@ -etotref -463.3717558791492479 -etotperatomref -154.4572519597 -totalforceref 3.661556 -totalstressref 9.690490 -deepks_e_label 17.02859183994637 +etotref -466.0484517111889318 +etotperatomref -155.3494839037 +totalforceref 3.432943 +totalstressref 7.674750 +deepks_e_label 17.12695856218435 deepks_edelta 0.0 -deepks_o_label 0.3103148500892669 -deepks_odelta 0.0 -deepks_oprec -0.42236872134754927 -deepks_h_label 49.03922791708764 +deepks_h_label 49.127076967385676 deepks_vdelta 0.0 -deepks_vdp 176.26236130239985 -deepks_f_label 0.07120588316980042 +deepks_vdp 176.0856172574296 +deepks_f_label 0.0667600367817765 deepks_fdelta 0.0 -deepks_fpre 19.53717375096855 -deepks_s_label 0.08886477565386258 +deepks_s_label 0.06851710382837811 deepks_sdelta 0.0 -deepks_spre 19.22629462223567 -deepks_e_label_elec 51.047596216236626 -deepks_edelta_elec 0 -deepks_o_label_elec .86018933752817077 -deepks_odelta_elec 0 -deepks_oprec_elec -1.23599151361439742 -deepks_h_label_elec 147.028760179109226 -deepks_vdelta_elec 0 -deepks_vdp_elec 529.07335498441665 -totaltimeref 1.06 +deepks_fpre 19.691711348715746 +deepks_spre 19.372132511297018 +deepks_e_label_elec 85.389836665688179 +deepks_edelta_elec .49279432475412220 +deepks_h_label_elec 246.793639290695098 +deepks_vdelta_elec 3.4234699193308906 +deepks_vdp_elec 881.06539843219701 +deepks_f_label_elec .0667600367817765 +deepks_fdelta_elec .017979102769778977 +deepks_s_label_elec .06851710382837811 +deepks_sdelta_elec .013725178909649378 +totaltimeref 3.52 diff --git a/tests/integrate/tools/catch_deepks_properties.sh b/tests/integrate/tools/catch_deepks_properties.sh index e49a3c6b26..f8bd7cbf0d 100644 --- a/tests/integrate/tools/catch_deepks_properties.sh +++ b/tests/integrate/tools/catch_deepks_properties.sh @@ -37,12 +37,23 @@ process_npy() { local total="0" local file_pattern="" local base_pattern="" + + # Check if force or stress + local is_force_or_stress=0 + if [ "$file_prefix" = "ftot" ] || [ "$file_prefix" = "stot" ]; then + is_force_or_stress=1 + fi # Determine file pattern based on mode if [ "$mode" = "multi" ]; then - # multi mode: multiple files with _e* suffix from different electronic steps - file_pattern="OUT.autotest/DeePKS_Labels_Elec/${file_prefix}_e*.npy" - base_pattern="OUT.autotest/DeePKS_Labels_Elec/${base_prefix}_e*.npy" + # multi mode: multiple files with _e* suffix from different electronic steps, unless for force and stress + if [ "$is_force_or_stress" -eq 1 ] ; then + file_pattern="OUT.autotest/DeePKS_Labels_Elec/${file_prefix}.npy" + base_pattern="OUT.autotest/DeePKS_Labels_Elec/${base_prefix}.npy" + else + file_pattern="OUT.autotest/DeePKS_Labels_Elec/${file_prefix}_e*.npy" + base_pattern="OUT.autotest/DeePKS_Labels_Elec/${base_prefix}_e*.npy" + fi elif [ "$mode" = "single" ];then # single mode: single file file_pattern="OUT.autotest/deepks_${file_prefix}.npy" @@ -63,7 +74,7 @@ process_npy() { # Get corresponding base file local base_file="" - if [ "$mode" = "multi" ]; then + if [ "$mode" = "multi" ] && [ "$is_force_or_stress" -ne 1 ]; then base_file="OUT.autotest/DeePKS_Labels_Elec/${base_prefix}_${step}.npy" else base_file=$base_pattern @@ -120,6 +131,17 @@ process_many_npys() { process_npy "$mode" "numpy" "gevdm" "" "deepks_gevdm$suffix" "$output_file" fi fi + + if ! test -z "$has_force" && [ $has_force == 1 ]; then + process_npy "$mode" "abs" "ftot" "" "deepks_f_label$suffix" "$output_file" + process_npy "$mode" "delta" "ftot" "fbase" "deepks_fdelta$suffix" "$output_file" + fi + + # For cal_stress = 1 + if ! test -z "$has_stress" && [ $has_stress == 1 ]; then + process_npy "$mode" "abs" "stot" "" "deepks_s_label$suffix" "$output_file" + process_npy "$mode" "delta" "stot" "sbase" "deepks_sdelta$suffix" "$output_file" + fi } # Main script @@ -150,18 +172,14 @@ fi if ! test -z "$deepks_out_labels" && [ $deepks_out_labels == 1 ]; then process_many_npys "single" "" "$1" - # force and stress not considered in deepks_out_freq_elec > 0, so not in process_many_npys + # gradvx and gvepsl not considered in deepks_out_freq_elec > 0, so not in process_many_npys # For cal_force = 1 if ! test -z "$has_force" && [ $has_force == 1 ]; then - process_npy "single" "abs" "ftot" "" "deepks_f_label" "$1" - process_npy "single" "delta" "ftot" "fbase" "deepks_fdelta" "$1" process_npy "single" "abs" "gradvx" "" "deepks_fpre" "$1" fi # For cal_stress = 1 if ! test -z "$has_stress" && [ $has_stress == 1 ]; then - process_npy "single" "abs" "stot" "" "deepks_s_label" "$1" - process_npy "single" "delta" "stot" "sbase" "deepks_sdelta" "$1" process_npy "single" "abs" "gvepsl" "" "deepks_spre" "$1" fi