diff --git a/source/module_elecstate/elecstate.h b/source/module_elecstate/elecstate.h index 401988d455..1a1e2f3013 100644 --- a/source/module_elecstate/elecstate.h +++ b/source/module_elecstate/elecstate.h @@ -1,10 +1,9 @@ #ifndef ELECSTATE_H #define ELECSTATE_H -#include "module_parameter/parameter.h" - #include "fp_energy.h" #include "module_cell/klist.h" #include "module_elecstate/module_charge/charge.h" +#include "module_parameter/parameter.h" #include "module_psi/psi.h" #include "potentials/potential_new.h" @@ -14,10 +13,10 @@ namespace elecstate class ElecState { public: - ElecState(){} - ElecState(Charge* charge_in, - ModulePW::PW_Basis* rhopw_in, - ModulePW::PW_Basis_Big* bigpw_in) + ElecState() + { + } + ElecState(Charge* charge_in, ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis_Big* bigpw_in) { this->charge = charge_in; this->charge->set_rhopw(rhopw_in); @@ -26,20 +25,20 @@ class ElecState } virtual ~ElecState() { - if(this->pot != nullptr) + if (this->pot != nullptr) { delete this->pot; this->pot = nullptr; } } - void init_ks(Charge *chg_in, // pointer for class Charge - const K_Vectors *klist_in, - int nk_in, // number of k points - ModulePW::PW_Basis* rhopw_in, - const ModulePW::PW_Basis_Big* bigpw_in); + void init_ks(Charge* chg_in, // pointer for class Charge + const K_Vectors* klist_in, + int nk_in, // number of k points + ModulePW::PW_Basis* rhopw_in, + const ModulePW::PW_Basis_Big* bigpw_in); // return current electronic density rho, as a input for constructing Hamiltonian - virtual const double *getRho(int spin) const; + virtual const double* getRho(int spin) const; // calculate electronic charge density on grid points or density matrix in real space // the consequence charge density rho saved into rho_out, preparing for charge mixing. @@ -78,17 +77,14 @@ class ElecState // use occupied weights from INPUT and skip calculate_weights // mohan updated on 2024-06-08 - void fixed_weights( - const std::vector& ocp_kb, - const int &nbands, - const double &nelec); + void fixed_weights(const std::vector& ocp_kb, const int& nbands, const double& nelec); - // if nupdown is not 0(TWO_EFERMI case), - // nelec_spin will be fixed and weights will be constrained + // if nupdown is not 0(TWO_EFERMI case), + // nelec_spin will be fixed and weights will be constrained void init_nelec_spin(); - //used to record number of electrons per spin index - //for NSPIN=2, it will record number of spin up and number of spin down - //for NSPIN=4, it will record total number, magnetization for x, y, z direction + // used to record number of electrons per spin index + // for NSPIN=2, it will record number of spin up and number of spin down + // for NSPIN=4, it will record total number, magnetization for x, y, z direction std::vector nelec_spin; virtual void print_psi(const psi::Psi& psi_in, const int istep = -1) @@ -102,7 +98,7 @@ class ElecState /** * @brief Init rho_core, init rho, renormalize rho, init pot - * + * * @param istep i-th step * @param ucell unit cell * @param strucfac structure factor @@ -142,7 +138,7 @@ class ElecState void set_exx(const std::complex& Eexx); #endif //__LCAO #endif //__EXX - + double get_hartree_energy(); double get_etot_efield(); double get_etot_gatefield(); @@ -150,20 +146,16 @@ class ElecState double get_solvent_model_Ael(); double get_solvent_model_Acav(); - virtual double get_spin_constrain_energy() { + virtual double get_spin_constrain_energy() + { return 0.0; } double get_dftu_energy(); double get_local_pp_energy(); -#ifdef __DEEPKS - double get_deepks_E_delta(); - double get_deepks_E_delta_band(); -#endif - - fenergy f_en; ///< energies contribute to the total free energy - efermi eferm; ///< fermi energies + fenergy f_en; ///< energies contribute to the total free energy + efermi eferm; ///< fermi energies // below defines the bandgap: diff --git a/source/module_elecstate/elecstate_energy.cpp b/source/module_elecstate/elecstate_energy.cpp index bbf88217a9..4f6e3e1756 100644 --- a/source/module_elecstate/elecstate_energy.cpp +++ b/source/module_elecstate/elecstate_energy.cpp @@ -1,10 +1,10 @@ -#include - #include "elecstate.h" #include "elecstate_getters.h" #include "module_base/global_variable.h" #include "module_base/parallel_reduce.h" #include "module_parameter/parameter.h" + +#include #ifdef USE_PAW #include "module_hamilt_general/module_xc/xc_functional.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" @@ -103,10 +103,9 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const const double* v_eff = this->pot->get_effective_v(0); const double* v_fixed = this->pot->get_fixed_v(); const double* v_ofk = nullptr; - const bool v_ofk_flag =(get_xc_func_type() == 3 - || get_xc_func_type() == 5); + const bool v_ofk_flag = (get_xc_func_type() == 3 || get_xc_func_type() == 5); #ifdef USE_PAW - if(PARAM.inp.use_paw) + if (PARAM.inp.use_paw) { ModuleBase::matrix v_xc; const std::tuple etxc_vtxc_v @@ -115,19 +114,19 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { - deband_aux -= this->charge->rho[0][ir] * v_xc(0,ir); + deband_aux -= this->charge->rho[0][ir] * v_xc(0, ir); } if (PARAM.inp.nspin == 2) { for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { - deband_aux -= this->charge->rho[1][ir] * v_xc(1,ir); + deband_aux -= this->charge->rho[1][ir] * v_xc(1, ir); } } } #endif - if(!PARAM.inp.use_paw) + if (!PARAM.inp.use_paw) { for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { @@ -137,16 +136,16 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const { v_ofk = this->pot->get_effective_vofk(0); // cause in the get_effective_vofk, the func will return nullptr - if(v_ofk==nullptr && this->charge->rhopw->nrxx>0) + if (v_ofk == nullptr && this->charge->rhopw->nrxx > 0) { - ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband","v_ofk is nullptr"); + ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband", "v_ofk is nullptr"); } for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { deband_aux -= this->charge->kin_r[0][ir] * v_ofk[ir]; } } - + if (PARAM.inp.nspin == 2) { v_eff = this->pot->get_effective_v(1); @@ -157,9 +156,9 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const if (v_ofk_flag) { v_ofk = this->pot->get_effective_vofk(1); - if(v_ofk==nullptr && this->charge->rhopw->nrxx>0) + if (v_ofk == nullptr && this->charge->rhopw->nrxx > 0) { - ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband","v_ofk is nullptr"); + ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband", "v_ofk is nullptr"); } for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { @@ -219,7 +218,7 @@ double ElecState::cal_delta_escf() const if (get_xc_func_type() == 3 || get_xc_func_type() == 5) { // cause in the get_effective_vofk, the func will return nullptr - assert(v_ofk!=nullptr); + assert(v_ofk != nullptr); descf -= (this->charge->kin_r[0][ir] - this->charge->kin_r_save[0][ir]) * v_ofk[ir]; } } @@ -278,7 +277,7 @@ void ElecState::cal_converged() /** * @brief calculate energies * - * @param type: 1 means Harris-Foulkes functinoal; + * @param type: 1 means Harris-Foulkes functinoal; * @param type: 2 means Kohn-Sham functional; */ void ElecState::cal_energies(const int type) @@ -292,7 +291,7 @@ void ElecState::cal_energies(const int type) //! energy from gate-field this->f_en.gatefield = get_etot_gatefield(); - //! energy from implicit solvation model + //! energy from implicit solvation model if (PARAM.inp.imp_sol) { this->f_en.esol_el = get_solvent_model_Ael(); @@ -305,7 +304,7 @@ void ElecState::cal_energies(const int type) this->f_en.escon = get_spin_constrain_energy(); } - // energy from DFT+U + // energy from DFT+U if (PARAM.inp.dft_plus_u) { this->f_en.edftu = get_dftu_energy(); @@ -313,19 +312,11 @@ void ElecState::cal_energies(const int type) this->f_en.e_local_pp = get_local_pp_energy(); -#ifdef __DEEPKS - // energy from deepks - if (PARAM.inp.deepks_scf) - { - this->f_en.edeepks_scf = get_deepks_E_delta() - get_deepks_E_delta_band(); - } -#endif - if (type == 1) // Harris-Foulkes functional { this->f_en.calculate_harris(); } - else if (type == 2)// Kohn-Sham functional + else if (type == 2) // Kohn-Sham functional { this->f_en.calculate_etot(); } diff --git a/source/module_elecstate/elecstate_energy_terms.cpp b/source/module_elecstate/elecstate_energy_terms.cpp index a0551a7bdd..cf3d460a02 100644 --- a/source/module_elecstate/elecstate_energy_terms.cpp +++ b/source/module_elecstate/elecstate_energy_terms.cpp @@ -3,8 +3,8 @@ #include "module_elecstate/potentials/efield.h" #include "module_elecstate/potentials/gatefield.h" #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" -#include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" +#include "module_hamilt_lcao/module_dftu/dftu.h" namespace elecstate { @@ -41,25 +41,15 @@ double ElecState::get_dftu_energy() double ElecState::get_local_pp_energy() { - double local_pseudopot_energy = 0.; // electron-ion interaction energy from local pseudopotential + double local_pseudopot_energy = 0.; // electron-ion interaction energy from local pseudopotential for (int is = 0; is < PARAM.inp.nspin; ++is) { - local_pseudopot_energy += BlasConnector::dot(this->charge->rhopw->nrxx, this->pot->get_fixed_v(), 1, this->charge->rho[is], 1) - * this->charge->rhopw->omega / this->charge->rhopw->nxyz; + local_pseudopot_energy + += BlasConnector::dot(this->charge->rhopw->nrxx, this->pot->get_fixed_v(), 1, this->charge->rho[is], 1) + * this->charge->rhopw->omega / this->charge->rhopw->nxyz; } Parallel_Reduce::reduce_all(local_pseudopot_energy); return local_pseudopot_energy; } -#ifdef __DEEPKS -double ElecState::get_deepks_E_delta() -{ - return GlobalC::ld.E_delta; -} -double ElecState::get_deepks_E_delta_band() -{ - return GlobalC::ld.e_delta_band; -} -#endif - } // namespace elecstate \ No newline at end of file diff --git a/source/module_elecstate/elecstate_print.cpp b/source/module_elecstate/elecstate_print.cpp index 0951379811..bbdebb7476 100644 --- a/source/module_elecstate/elecstate_print.cpp +++ b/source/module_elecstate/elecstate_print.cpp @@ -1,6 +1,5 @@ #include "elecstate.h" #include "elecstate_getters.h" -#include "module_parameter/parameter.h" #include "module_base/formatter.h" #include "module_base/global_variable.h" #include "module_elecstate/potentials/H_Hartree_pw.h" @@ -8,6 +7,7 @@ #include "module_elecstate/potentials/gatefield.h" #include "module_hamilt_general/module_xc/xc_functional.h" #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#include "module_parameter/parameter.h" #include "occupy.h" namespace elecstate { @@ -311,9 +311,10 @@ void ElecState::print_etot(const Magnetism& magnet, GlobalV::ofs_running << "\n Density error is " << scf_thr << std::endl; - if (PARAM.inp.basis_type == "pw") { + if (PARAM.inp.basis_type == "pw") + { ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Error Threshold", pw_diag_thr); // xiaohui add 2013-09-02 -} + } std::vector titles; std::vector energies_Ry; @@ -378,7 +379,7 @@ void ElecState::print_etot(const Magnetism& magnet, if (PARAM.inp.deepks_scf) // caoyu add 2021-08-10 { titles.push_back("E_DeePKS"); - energies_Ry.push_back(GlobalC::ld.E_delta); + energies_Ry.push_back(this->f_en.edeepks_delta); } #endif } diff --git a/source/module_elecstate/fp_energy.h b/source/module_elecstate/fp_energy.h index adbb3fe3bd..17fad4c8be 100644 --- a/source/module_elecstate/fp_energy.h +++ b/source/module_elecstate/fp_energy.h @@ -39,12 +39,13 @@ struct fenergy double esol_el = 0.0; ///< the implicit solvation energy Ael double esol_cav = 0.0; ///< the implicit solvation energy Acav - double edftu = 0.0; ///< DFT+U energy - double edeepks_scf = 0.0; /// DeePKS energy + double edftu = 0.0; ///< DFT+U energy + double edeepks_scf = 0.0; /// DeePKS energy difference + double edeepks_delta = 0.0; /// DeePKS energy double escon = 0.0; ///< spin constraint energy - double ekinetic = 0.0; /// kinetic energy, used in OFDFT + double ekinetic = 0.0; /// kinetic energy, used in OFDFT double e_local_pp = 0.0; /// ion-electron interaction energy contributed by local pp, used in OFDFT double calculate_etot(); diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 1fc3b1cdc3..255526a97d 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -37,7 +37,7 @@ #include "module_elecstate/elecstate_lcao_cal_tau.h" #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_elecstate/occupy.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need divide_HS_in_frag +#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need DeePKS_init #include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp" #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" @@ -113,15 +113,14 @@ ESolver_KS_LCAO::~ESolver_KS_LCAO() //! 2) init ElecState //! 3) init LCAO basis //! 4) initialize the density matrix -//! 5) initialize Hamilt in LCAO -//! 6) initialize exx -//! 7) initialize DFT+U -//! 8) ppcell -//! 9) inititlize the charge density -//! 10) initialize the potential. -//! 11) initialize deepks -//! 12) set occupations -//! 13) print a warning if needed +//! 5) initialize exx +//! 6) initialize DFT+U +//! 7) ppcell +//! 8) inititlize the charge density +//! 9) initialize the potential. +//! 10) initialize deepks +//! 11) set occupations +//! 12) print a warning if needed //------------------------------------------------------------------------------ template void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_para& inp) @@ -163,13 +162,8 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa // DMR is not initialized here, it will be constructed in each before_scf dynamic_cast*>(this->pelec)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin); - // 5) initialize Hamilt in LCAO - // * allocate H and S matrices according to computational resources - // * set the 'trace' between local H/S and global H/S - LCAO_domain::divide_HS_in_frag(PARAM.globalv.gamma_only_local, ucell, pv, this->kv.get_nks(), orb_); - #ifdef __EXX - // 6) initialize exx + // 5) initialize exx // PLEASE simplify the Exx_Global interface if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" || PARAM.inp.calculation == "md") @@ -192,22 +186,22 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } #endif - // 7) initialize DFT+U + // 6) initialize DFT+U if (PARAM.inp.dft_plus_u) { auto* dftu = ModuleDFTU::DFTU::get_instance(); dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_); } - // 8) initialize ppcell + // 7) initialize ppcell this->locpp.init_vloc(ucell, this->pw_rho); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); - // 9) inititlize the charge density + // 8) inititlize the charge density this->pelec->charge->allocate(PARAM.inp.nspin); this->pelec->omega = ucell.omega; - // 10) initialize the potential + // 9) initialize the potential if (this->pelec->pot == nullptr) { this->pelec->pot = new elecstate::Potential(this->pw_rhod, @@ -221,31 +215,32 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } #ifdef __DEEPKS - // 11) initialize deepks + // 10) initialize deepks + LCAO_domain::DeePKS_init(ucell, pv, this->kv.get_nks(), orb_, this->ld); if (PARAM.inp.deepks_scf) { // load the DeePKS model from deep neural network - DeePKS_domain::load_model(PARAM.inp.deepks_model, GlobalC::ld.model_deepks); + DeePKS_domain::load_model(PARAM.inp.deepks_model, ld.model_deepks); // read pdm from file for NSCF or SCF-restart, do it only once in whole calculation DeePKS_domain::read_pdm((PARAM.inp.init_chg == "file"), PARAM.inp.deepks_equiv, - GlobalC::ld.init_pdm, - GlobalC::ld.inlmax, - GlobalC::ld.lmaxd, - GlobalC::ld.inl_l, + ld.init_pdm, + orb_.Alpha[0].getTotal_nchi() * ucell.nat, + ld.lmaxd, + ld.inl_l, *orb_.Alpha, - GlobalC::ld.pdm); + ld.pdm); } #endif - // 12) set occupations + // 11) set occupations // tddft does not need to set occupations in the first scf if (PARAM.inp.ocp && inp.esolver_type != "tddft") { this->pelec->fixed_weights(PARAM.inp.ocp_kb, PARAM.inp.nbands, PARAM.inp.nelec); } - // 13) if kpar is not divisible by nks, print a warning + // 12) if kpar is not divisible by nks, print a warning if (GlobalV::KPAR_LCAO > 1) { if (this->kv.get_nks() % GlobalV::KPAR_LCAO != 0) @@ -266,7 +261,7 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa } } - // 14) initialize rdmft, added by jghan + // 13) initialize rdmft, added by jghan if (PARAM.inp.rdmft == true) { rdmft_solver.init(this->GG, @@ -330,6 +325,9 @@ void ESolver_KS_LCAO::cal_force(UnitCell& ucell, ModuleBase::matrix& for this->kv, this->pw_rho, this->solvent, +#ifdef __DEEPKS + this->ld, +#endif #ifdef __EXX *this->exx_lri_double, *this->exx_lri_complex, @@ -675,7 +673,7 @@ void ESolver_KS_LCAO::iter_init(UnitCell& ucell, const int istep, const #ifdef __DEEPKS // the density matrixes of DeePKS have been updated in each iter - GlobalC::ld.set_hr_cal(true); + ld.set_hr_cal(true); // HR in HamiltLCAO should be recalculate if (PARAM.inp.deepks_scf) @@ -841,7 +839,9 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& const std::vector>& dm = dynamic_cast*>(this->pelec)->get_DM()->get_DMK_vector(); - GlobalC::ld.dpks_cal_e_delta_band(dm, this->kv.get_nks()); + ld.dpks_cal_e_delta_band(dm, this->kv.get_nks()); + this->pelec->f_en.edeepks_scf = ld.E_delta - ld.e_delta_band; + this->pelec->f_en.edeepks_delta = ld.E_delta; } #endif @@ -1049,7 +1049,7 @@ void ESolver_KS_LCAO::after_scf(UnitCell& ucell, const int istep) if (this->psi != nullptr && (istep % PARAM.inp.out_interval == 0)) { hamilt::HamiltLCAO* p_ham_deepks = dynamic_cast*>(this->p_hamilt); - std::shared_ptr ld_shared_ptr(&GlobalC::ld, [](LCAO_Deepks*) {}); + std::shared_ptr ld_shared_ptr(&ld, [](LCAO_Deepks*) {}); LCAO_Deepks_Interface LDI(ld_shared_ptr); LDI.out_deepks_labels(this->pelec->f_en.etot, diff --git a/source/module_esolver/esolver_ks_lcao.h b/source/module_esolver/esolver_ks_lcao.h index 1730f04a36..ed2c180bed 100644 --- a/source/module_esolver/esolver_ks_lcao.h +++ b/source/module_esolver/esolver_ks_lcao.h @@ -5,6 +5,9 @@ // for grid integration #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" +#ifdef __DEEPKS +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#endif #ifdef __EXX #include "module_ri/Exx_LRI_interface.h" #include "module_ri/Mix_DMk_2D.h" @@ -19,13 +22,14 @@ namespace LR { - template - class ESolver_LR; +template +class ESolver_LR; } namespace ModuleESolver { template -class ESolver_KS_LCAO : public ESolver_KS { +class ESolver_KS_LCAO : public ESolver_KS +{ public: ESolver_KS_LCAO(); ~ESolver_KS_LCAO(); @@ -73,7 +77,7 @@ class ESolver_KS_LCAO : public ESolver_KS { TwoCenterBundle two_center_bundle_; - rdmft::RDMFT rdmft_solver; // added by jghan for rdmft calculation, 2024-03-16 + rdmft::RDMFT rdmft_solver; // added by jghan for rdmft calculation, 2024-03-16 // temporary introduced during removing GlobalC::ORB LCAO_Orbitals orb_; @@ -92,6 +96,10 @@ class ESolver_KS_LCAO : public ESolver_KS { void beforesolver(const int istep); //--------------------------------------------------------------------- +#ifdef __DEEPKS + LCAO_Deepks ld; +#endif + #ifdef __EXX std::shared_ptr> exd = nullptr; std::shared_ptr>> exc = nullptr; diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 6b50347716..84a4ed0c68 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -17,19 +17,18 @@ #include "module_elecstate/module_dm/cal_edm_tddft.h" #include "module_elecstate/module_dm/density_matrix.h" #include "module_elecstate/occupy.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" // need divide_HS_in_frag #include "module_hamilt_lcao/module_tddft/evolve_elec.h" #include "module_hamilt_lcao/module_tddft/td_velocity.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" //-----HSolver ElecState Hamilt-------- +#include "module_elecstate/cal_ux.h" #include "module_elecstate/elecstate_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" #include "module_hsolver/hsolver_lcao.h" #include "module_parameter/parameter.h" #include "module_psi/psi.h" -#include "module_elecstate/cal_ux.h" //-----force& stress------------------- #include "module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h" @@ -290,7 +289,12 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep) { std::stringstream ss_dipole; ss_dipole << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_DIPOLE"; - ModuleIO::write_dipole(ucell,pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str()); + ModuleIO::write_dipole(ucell, + pelec->charge->rho_save[is], + pelec->charge->rhopw, + is, + istep, + ss_dipole.str()); } } if (TD_Velocity::out_current == true) diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index 724b4e06e3..2637fe41d8 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/source/module_esolver/lcao_before_scf.cpp @@ -194,6 +194,10 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) two_center_bundle_, orb_, DM +#ifdef __DEEPKS + , + &this->ld +#endif #ifdef __EXX , istep, @@ -211,7 +215,7 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) { const Parallel_Orbitals* pv = &this->pv; // allocate , phialpha is different every ion step, so it is allocated here - DeePKS_domain::allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, GlobalC::ld.phialpha); + DeePKS_domain::allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, this->ld.phialpha); // build and save at beginning DeePKS_domain::build_phialpha(PARAM.inp.cal_force, ucell, @@ -219,11 +223,11 @@ void ESolver_KS_LCAO::before_scf(UnitCell& ucell, const int istep) this->gd, pv, *(two_center_bundle_.overlap_orb_alpha), - GlobalC::ld.phialpha); + this->ld.phialpha); if (PARAM.inp.deepks_out_unittest) { - DeePKS_domain::check_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, GlobalC::ld.phialpha); + DeePKS_domain::check_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, this->ld.phialpha); } } #endif diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 2ba5d9c37c..fc0ab246d3 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -200,6 +200,10 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) two_center_bundle_, orb_, DM +#ifdef __DEEPKS + , + &this->ld +#endif #ifdef __EXX , istep, @@ -217,7 +221,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) { const Parallel_Orbitals* pv = &this->pv; // allocate , phialpha is different every ion step, so it is allocated here - DeePKS_domain::allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, GlobalC::ld.phialpha); + DeePKS_domain::allocate_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, this->ld.phialpha); // build and save at beginning DeePKS_domain::build_phialpha(PARAM.inp.cal_force, ucell, @@ -225,11 +229,11 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) this->gd, pv, *(two_center_bundle_.overlap_orb_alpha), - GlobalC::ld.phialpha); + this->ld.phialpha); if (PARAM.inp.deepks_out_unittest) { - DeePKS_domain::check_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, GlobalC::ld.phialpha); + DeePKS_domain::check_phialpha(PARAM.inp.cal_force, ucell, orb_, this->gd, pv, this->ld.phialpha); } } #endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h index 68a6359fd3..21dc99d758 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h @@ -5,13 +5,14 @@ #include "module_base/global_variable.h" #include "module_base/matrix.h" #include "module_basis/module_nao/two_center_bundle.h" +#include "module_elecstate/elecstate.h" #include "module_elecstate/module_dm/density_matrix.h" +#include "module_elecstate/potentials/potential_new.h" #include "module_hamilt_lcao/hamilt_lcaodft/force_stress_arrays.h" +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_psi/psi.h" -#include "module_elecstate/potentials/potential_new.h" -#include "module_elecstate/elecstate.h" #ifndef TGINT_H #define TGINT_H @@ -65,6 +66,7 @@ class Force_LCAO #ifdef __DEEPKS ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, + LCAO_Deepks& ld, #endif typename TGint::type& gint, const TwoCenterBundle& two_center_bundle, @@ -87,7 +89,7 @@ class Force_LCAO void average_force(double* fm); - //void test(Parallel_Orbitals& pv, double* mm, const std::string& name); + // void test(Parallel_Orbitals& pv, double* mm, const std::string& name); //------------------------------------------------------------- // forces reated to overlap matrix @@ -131,14 +133,14 @@ class Force_LCAO ModuleBase::matrix& svl_dphi); elecstate::DensityMatrix cal_edm(const elecstate::ElecState* pelec, - const psi::Psi& psi, - const elecstate::DensityMatrix& dm, - const K_Vectors& kv, - const Parallel_Orbitals& pv, - const int& nspin, - const int& nbands, - const UnitCell& ucell, - Record_adj& ra) const; + const psi::Psi& psi, + const elecstate::DensityMatrix& dm, + const K_Vectors& kv, + const Parallel_Orbitals& pv, + const int& nspin, + const int& nbands, + const UnitCell& ucell, + Record_adj& ra) const; }; #endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 1d64bbe049..0aa81c2537 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -50,6 +50,9 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, const K_Vectors& kv, ModulePW::PW_Basis* rhopw, surchem& solvent, +#ifdef __DEEPKS + LCAO_Deepks& ld, +#endif #ifdef __EXX Exx_LRI& exx_lri_double, Exx_LRI>& exx_lri_complex, @@ -180,6 +183,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, #ifdef __DEEPKS fvnl_dalpha, svnl_dalpha, + ld, #endif gint_gamma, gint_k, @@ -500,16 +504,16 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, if (PARAM.inp.deepks_out_labels) // not parallelized yet { const std::string file_ftot = PARAM.globalv.global_out_dir + "deepks_ftot.npy"; - LCAO_deepks_io::save_npy_f(fcs, file_ftot, GlobalV::MY_RANK); // Ry/Bohr, F_tot + LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Ry/Bohr, F_tot const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; if (PARAM.inp.deepks_scf) { - LCAO_deepks_io::save_npy_f(fcs - fvnl_dalpha, file_fbase, GlobalV::MY_RANK); // Ry/Bohr, F_base + LCAO_deepks_io::save_matrix2npy(file_fbase, fcs - fvnl_dalpha, GlobalV::MY_RANK); // Ry/Bohr, F_base } else { - LCAO_deepks_io::save_npy_f(fcs, file_fbase, GlobalV::MY_RANK); // no scf, F_base=F_tot + LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot } } #endif @@ -682,23 +686,24 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, if (PARAM.inp.deepks_out_labels) // not parallelized yet { const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; - LCAO_deepks_io::save_npy_s(scs, - file_stot, - ucell.omega, - GlobalV::MY_RANK); // change to energy unit Ry when printing, S_tot, w/ model + 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) { - LCAO_deepks_io::save_npy_s(scs - svnl_dalpha, - file_sbase, - ucell.omega, - GlobalV::MY_RANK); // change to energy unit Ry when printing, S_base; + 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_npy_s(scs, file_sbase, ucell.omega, - GlobalV::MY_RANK); // sbase = stot + LCAO_deepks_io::save_matrix2npy(file_sbase, scs, GlobalV::MY_RANK, ucell.omega, 'U'); // sbase = stot } } #endif @@ -832,6 +837,7 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, #if __DEEPKS ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, + LCAO_Deepks& ld, #endif Gint_Gamma& gint_gamma, // mohan add 2024-04-01 Gint_k& gint_k, // mohan add 2024-04-01 @@ -859,6 +865,7 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, #if __DEEPKS fvnl_dalpha, svnl_dalpha, + ld, #endif gint_gamma, two_center_bundle, @@ -887,6 +894,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn #if __DEEPKS ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, + LCAO_Deepks& ld, #endif Gint_Gamma& gint_gamma, Gint_k& gint_k, @@ -913,6 +921,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn #if __DEEPKS fvnl_dalpha, svnl_dalpha, + ld, #endif gint_k, two_center_bundle, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h index 89a9bf6d55..6de723d868 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h @@ -49,6 +49,9 @@ class Force_Stress_LCAO const K_Vectors& kv, ModulePW::PW_Basis* rhopw, surchem& solvent, +#ifdef __DEEPKS + LCAO_Deepks& ld, +#endif #ifdef __EXX Exx_LRI& exx_lri_double, Exx_LRI>& exx_lri_complex, @@ -62,9 +65,7 @@ class Force_Stress_LCAO Stress_Func sc_pw; Forces f_pw; - void forceSymmetry(const UnitCell& ucell, - ModuleBase::matrix& fcs, - ModuleSymmetry::Symmetry* symm); + void forceSymmetry(const UnitCell& ucell, ModuleBase::matrix& fcs, ModuleSymmetry::Symmetry* symm); void calForcePwPart(UnitCell& ucell, ModuleBase::matrix& fvl_dvl, @@ -98,6 +99,7 @@ class Force_Stress_LCAO #if __DEEPKS ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, + LCAO_Deepks& ld, #endif Gint_Gamma& gint_gamma, Gint_k& gint_k, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index 506a1ada42..d8b913a61a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -190,6 +190,7 @@ void Force_LCAO::ftable(const bool isforce, #ifdef __DEEPKS ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, + LCAO_Deepks& ld, #endif TGint::type& gint, const TwoCenterBundle& two_center_bundle, @@ -247,28 +248,23 @@ void Force_LCAO::ftable(const bool isforce, false /*reset dm to gint*/); #ifdef __DEEPKS - const std::vector>& dm_gamma = dm->get_DMK_vector(); - std::vector descriptor; if (PARAM.inp.deepks_scf) { + const std::vector>& dm_gamma = dm->get_DMK_vector(); + std::vector descriptor; // when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm - DeePKS_domain::cal_descriptor(ucell.nat, - GlobalC::ld.inlmax, - GlobalC::ld.inl_l, - GlobalC::ld.pdm, - descriptor, - GlobalC::ld.des_per_atom); + DeePKS_domain::cal_descriptor(ucell.nat, ld.inlmax, ld.inl_l, ld.pdm, descriptor, ld.des_per_atom); DeePKS_domain::cal_gedm(ucell.nat, - GlobalC::ld.lmaxd, - GlobalC::ld.nmaxd, - GlobalC::ld.inlmax, - GlobalC::ld.des_per_atom, - GlobalC::ld.inl_l, + ld.lmaxd, + ld.nmaxd, + ld.inlmax, + ld.des_per_atom, + ld.inl_l, descriptor, - GlobalC::ld.pdm, - GlobalC::ld.model_deepks, - GlobalC::ld.gedm, - GlobalC::ld.E_delta); + ld.pdm, + ld.model_deepks, + ld.gedm, + ld.E_delta); const int nks = 1; DeePKS_domain::cal_f_delta(dm_gamma, @@ -276,12 +272,11 @@ void Force_LCAO::ftable(const bool isforce, orb, gd, *this->ParaV, - GlobalC::ld.lmaxd, nks, kv->kvec_d, - GlobalC::ld.phialpha, - GlobalC::ld.gedm, - GlobalC::ld.inl_index, + ld.phialpha, + ld.gedm, + ld.inl_index, fvnl_dalpha, isstress, svnl_dalpha); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index ac44ee4194..7ec1f427d0 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -284,6 +284,7 @@ void Force_LCAO>::ftable(const bool isforce, #ifdef __DEEPKS ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, + LCAO_Deepks& ld, #endif TGint>::type& gint, const TwoCenterBundle& two_center_bundle, @@ -345,38 +346,31 @@ void Force_LCAO>::ftable(const bool isforce, if (PARAM.inp.deepks_scf) { const std::vector>>& dm_k = dm->get_DMK_vector(); - - // when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm std::vector descriptor; - DeePKS_domain::cal_descriptor(ucell.nat, - GlobalC::ld.inlmax, - GlobalC::ld.inl_l, - GlobalC::ld.pdm, - descriptor, - GlobalC::ld.des_per_atom); + // when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm + DeePKS_domain::cal_descriptor(ucell.nat, ld.inlmax, ld.inl_l, ld.pdm, descriptor, ld.des_per_atom); DeePKS_domain::cal_gedm(ucell.nat, - GlobalC::ld.lmaxd, - GlobalC::ld.nmaxd, - GlobalC::ld.inlmax, - GlobalC::ld.des_per_atom, - GlobalC::ld.inl_l, + ld.lmaxd, + ld.nmaxd, + ld.inlmax, + ld.des_per_atom, + ld.inl_l, descriptor, - GlobalC::ld.pdm, - GlobalC::ld.model_deepks, - GlobalC::ld.gedm, - GlobalC::ld.E_delta); + ld.pdm, + ld.model_deepks, + ld.gedm, + ld.E_delta); DeePKS_domain::cal_f_delta>(dm_k, ucell, orb, gd, pv, - GlobalC::ld.lmaxd, kv->get_nks(), kv->kvec_d, - GlobalC::ld.phialpha, - GlobalC::ld.gedm, - GlobalC::ld.inl_index, + ld.phialpha, + ld.gedm, + ld.inl_index, fvnl_dalpha, isstress, svnl_dalpha); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp index c7b7e6ee10..ec43ac6f2a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_allocate.cpp @@ -1,46 +1,39 @@ -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_base/timer.h" -#include "module_parameter/parameter.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_parameter/parameter.h" namespace LCAO_domain { - -void divide_HS_in_frag(const bool isGamma, const UnitCell& ucell, Parallel_Orbitals& pv,const int& nks, const LCAO_Orbitals& orb) { - ModuleBase::TITLE("LCAO_domain", "divide_HS_in_frag"); - - //(1), (2): set up matrix division have been moved into ESolver_KS_LCAO::init_basis_lcao - // just pass `ParaV` as pointer is enough #ifdef __DEEPKS - // wenfei 2021-12-19 +// It seems it is only related to DeePKS, so maybe we should move it to DeeKS_domain +void DeePKS_init(const UnitCell& ucell, + Parallel_Orbitals& pv, + const int& nks, + const LCAO_Orbitals& orb, + LCAO_Deepks& ld) +{ + ModuleBase::TITLE("LCAO_domain", "DeePKS_init"); // preparation for DeePKS - - if (PARAM.inp.deepks_out_labels || PARAM.inp.deepks_scf) { + if (PARAM.inp.deepks_out_labels || PARAM.inp.deepks_scf) + { // allocate relevant data structures for calculating descriptors std::vector na; na.resize(ucell.ntype); - for (int it = 0; it < ucell.ntype; it++) { + for (int it = 0; it < ucell.ntype; it++) + { na[it] = ucell.atoms[it].na; } - GlobalC::ld.init(orb, - ucell.nat, - ucell.ntype, - nks, - pv, - na); + ld.init(orb, ucell.nat, ucell.ntype, nks, pv, na); - if (PARAM.inp.deepks_scf) { - if (isGamma) { - GlobalC::ld.allocate_V_delta(ucell.nat); - } else { - GlobalC::ld.allocate_V_delta(ucell.nat, nks); - } + if (PARAM.inp.deepks_scf) + { + ld.allocate_V_delta(ucell.nat, nks); } } -#endif return; } - -} +#endif +} // namespace LCAO_domain diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h index 4b02f1dc54..992c754e77 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h @@ -8,6 +8,7 @@ #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_HS_arrays.hpp" #include "module_hamilt_lcao/hamilt_lcaodft/force_stress_arrays.h" +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" #include "module_hamilt_lcao/module_gint/grid_technique.h" @@ -16,14 +17,14 @@ namespace LCAO_domain { void init_basis_lcao(Parallel_Orbitals& pv, - const double &onsite_radius, - const double &lcao_ecut, - const double &lcao_dk, - const double &lcao_dr, - const double &lcao_rmax, - UnitCell& ucell, - TwoCenterBundle& two_center_bundle, - LCAO_Orbitals& orb); + const double& onsite_radius, + const double& lcao_ecut, + const double& lcao_dk, + const double& lcao_dr, + const double& lcao_rmax, + UnitCell& ucell, + TwoCenterBundle& two_center_bundle, + LCAO_Orbitals& orb); void build_Nonlocal_mu_new(const Parallel_Orbitals& pv, ForceStressArrays& fsr, // mohan 2024-06-16 @@ -173,14 +174,16 @@ void build_ST_new(ForceStressArrays& fsr, */ void zeros_HSR(const char& mtype, LCAO_HS_Arrays& HS_arrays); -void divide_HS_in_frag(const bool isGamma, const UnitCell& ucell, Parallel_Orbitals& pv, const int& nks, const LCAO_Orbitals& orb); +#ifdef __DEEPKS +void DeePKS_init(const UnitCell& ucell, + Parallel_Orbitals& pv, + const int& nks, + const LCAO_Orbitals& orb, + LCAO_Deepks& ld); +#endif template -void set_mat2d(const int& global_ir, - const int& global_ic, - const T& v, - const Parallel_Orbitals& pv, - T* mat); +void set_mat2d(const int& global_ir, const int& global_ic, const T& v, const Parallel_Orbitals& pv, T* mat); } // namespace LCAO_domain diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 06ab104c84..69e90ea3b9 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -81,6 +81,10 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, elecstate::DensityMatrix* DM_in +#ifdef __DEEPKS + , + LCAO_Deepks* ld_in +#endif #ifdef __EXX , const int istep, @@ -209,7 +213,8 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, two_center_bundle.overlap_orb_alpha.get(), &orb, this->kv->get_nks(), - DM_in); + DM_in, + ld_in); this->getOperator()->add(deepks); } #endif @@ -262,7 +267,6 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, orb.cutoffs(), &grid_d, PARAM.inp.nspin); - } } @@ -333,7 +337,8 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, two_center_bundle.overlap_orb_alpha.get(), &orb, this->kv->get_nks(), - DM_in); + DM_in, + ld_in); this->getOperator()->add(deepks); } #endif diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h index 65e95e22af..303a93e3c3 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h @@ -40,6 +40,10 @@ class HamiltLCAO : public Hamilt const TwoCenterBundle& two_center_bundle, const LCAO_Orbitals& orb, elecstate::DensityMatrix* DM_in +#ifdef __DEEPKS + , + LCAO_Deepks* ld_in +#endif #ifdef __EXX , const int istep, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp index 09410cf101..fe3ca7b364 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.cpp @@ -2,12 +2,10 @@ #include "module_base/timer.h" #include "module_base/tool_title.h" -#include "module_parameter/parameter.h" -#ifdef __DEEPKS -#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" -#endif #include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h" +#include "module_parameter/parameter.h" #ifdef _OPENMP #include #endif @@ -24,13 +22,19 @@ DeePKS>::DeePKS(HS_Matrix_K* hsk_in, const TwoCenterIntegrator* intor_orb_alpha, const LCAO_Orbitals* ptr_orb, const int& nks_in, - elecstate::DensityMatrix* DM_in) + elecstate::DensityMatrix* DM_in +#ifdef __DEEPKS + , + LCAO_Deepks* ld_in +#endif + ) : OperatorLCAO(hsk_in, kvec_d_in, hR_in), DM(DM_in), ucell(ucell_in), intor_orb_alpha_(intor_orb_alpha), ptr_orb_(ptr_orb), nks(nks_in) { this->cal_type = calculation_type::lcao_deepks; this->gd = GridD_in; #ifdef __DEEPKS + this->ld = ld_in; this->initialize_HR(GridD_in); #endif } @@ -154,43 +158,45 @@ void hamilt::DeePKS>::contributeHR() #ifdef __DEEPKS ModuleBase::TITLE("DeePKS", "contributeHR"); // if DM changed, HR of DeePKS need to refresh. - // the judgement is based on the status of HR in GlobalC::ld + // the judgement is based on the status of HR in ld // this operator should be informed that DM has changed and HR need to recalculate. - if (GlobalC::ld.get_hr_cal()) + if (this->ld->get_hr_cal()) { ModuleBase::timer::tick("DeePKS", "contributeHR"); - DeePKS_domain::cal_pdm(GlobalC::ld.init_pdm, - GlobalC::ld.inlmax, - GlobalC::ld.lmaxd, - GlobalC::ld.inl_l, - GlobalC::ld.inl_index, + const int inlmax = ptr_orb_->Alpha[0].getTotal_nchi() * this->ucell->nat; + + DeePKS_domain::cal_pdm(this->ld->init_pdm, + inlmax, + this->ld->lmaxd, + this->ld->inl_l, + this->ld->inl_index, this->DM, - GlobalC::ld.phialpha, + this->ld->phialpha, *this->ucell, *ptr_orb_, *(this->gd), *(this->hR->get_paraV()), - GlobalC::ld.pdm); + this->ld->pdm); std::vector descriptor; DeePKS_domain::cal_descriptor(this->ucell->nat, - GlobalC::ld.inlmax, - GlobalC::ld.inl_l, - GlobalC::ld.pdm, + inlmax, + this->ld->inl_l, + this->ld->pdm, descriptor, - GlobalC::ld.des_per_atom); + this->ld->des_per_atom); DeePKS_domain::cal_gedm(this->ucell->nat, - GlobalC::ld.lmaxd, - GlobalC::ld.nmaxd, - GlobalC::ld.inlmax, - GlobalC::ld.des_per_atom, - GlobalC::ld.inl_l, + this->ld->lmaxd, + this->ld->nmaxd, + inlmax, + this->ld->des_per_atom, + this->ld->inl_l, descriptor, - GlobalC::ld.pdm, - GlobalC::ld.model_deepks, - GlobalC::ld.gedm, - GlobalC::ld.E_delta); + this->ld->pdm, + this->ld->model_deepks, + this->ld->gedm, + this->ld->E_delta); // // recalculate the H_V_delta // if (this->H_V_delta == nullptr) @@ -200,7 +206,7 @@ void hamilt::DeePKS>::contributeHR() this->H_V_delta->set_zero(); this->calculate_HR(); - GlobalC::ld.set_hr_cal(false); + this->ld->set_hr_cal(false); ModuleBase::timer::tick("DeePKS", "contributeHR"); } @@ -297,8 +303,8 @@ void hamilt::DeePKS>::calculate_HR() { for (int N0 = 0; N0 < ptr_orb_->Alpha[0].getNchi(L0); ++N0) { - const int inl = GlobalC::ld.get_inl(T0, I0, L0, N0); - const double* pgedm = GlobalC::ld.get_gedms(inl); + const int inl = this->ld->inl_index[T0](I0, L0, N0); + const double* pgedm = this->ld->gedm[inl]; const int nm = 2 * L0 + 1; for (int m1 = 0; m1 < nm; ++m1) // m1 = 1 for s, 3 for p, 5 for d @@ -316,9 +322,9 @@ void hamilt::DeePKS>::calculate_HR() } else { - const double* pgedm = GlobalC::ld.get_gedms(iat0); + const double* pgedm = this->ld->gedm[iat0]; int nproj = 0; - for (int il = 0; il < GlobalC::ld.get_lmaxd() + 1; il++) + for (int il = 0; il < this->ld->lmaxd + 1; il++) { nproj += (2 * il + 1) * ptr_orb_->Alpha[0].getNchi(il); } @@ -466,15 +472,15 @@ void hamilt::DeePKS>::cal_HR_IJR(const double* hr_i } template -inline void get_h_delta_k(int ik, TK*& h_delta_k) +inline void get_h_delta_k(int ik, TK*& h_delta_k, LCAO_Deepks* ld_in) { if constexpr (std::is_same::value) { - h_delta_k = GlobalC::ld.H_V_delta[ik].data(); + h_delta_k = ld_in->H_V_delta[ik].data(); } else { - h_delta_k = GlobalC::ld.H_V_delta_k[ik].data(); + h_delta_k = ld_in->H_V_delta_k[ik].data(); } return; } @@ -487,7 +493,7 @@ void hamilt::DeePKS>::contributeHk(int ik) ModuleBase::timer::tick("DeePKS", "contributeHk"); TK* h_delta_k = nullptr; - get_h_delta_k(ik, h_delta_k); + get_h_delta_k(ik, h_delta_k, this->ld); // set SK to zero and then calculate SK for each k vector ModuleBase::GlobalFunc::ZEROS(h_delta_k, this->hsk->get_size()); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h index 08e19a6c9a..fd187b3fc5 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h @@ -4,8 +4,8 @@ #include "module_basis/module_nao/two_center_integrator.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_elecstate/module_dm/density_matrix.h" -#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" #include "operator_lcao.h" namespace hamilt @@ -38,7 +38,12 @@ class DeePKS> : public OperatorLCAO const TwoCenterIntegrator* intor_orb_alpha, const LCAO_Orbitals* ptr_orb, const int& nks_in, - elecstate::DensityMatrix* DM_in); + elecstate::DensityMatrix* DM_in +#ifdef __DEEPKS + , + LCAO_Deepks* ld_in +#endif + ); ~DeePKS(); /** diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp index 8d9e363d47..522810e7cb 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp @@ -15,15 +15,9 @@ #include "LCAO_deepks.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" -namespace GlobalC -{ -LCAO_Deepks ld; -} - // Constructor of the class LCAO_Deepks::LCAO_Deepks() { - alpha_index = new ModuleBase::IntArray[1]; inl_index = new ModuleBase::IntArray[1]; inl_l = nullptr; gedm = nullptr; @@ -33,7 +27,6 @@ LCAO_Deepks::LCAO_Deepks() // Desctructor of the class LCAO_Deepks::~LCAO_Deepks() { - delete[] alpha_index; delete[] inl_index; delete[] inl_l; @@ -139,8 +132,6 @@ void LCAO_Deepks::init_index(const int ntype, const int Total_nchi, const LCAO_Orbitals& orb) { - delete[] this->alpha_index; - this->alpha_index = new ModuleBase::IntArray[ntype]; delete[] this->inl_index; this->inl_index = new ModuleBase::IntArray[ntype]; delete[] this->inl_l; @@ -151,11 +142,6 @@ void LCAO_Deepks::init_index(const int ntype, int alpha = 0; for (int it = 0; it < ntype; it++) { - this->alpha_index[it].create(na[it], - this->lmaxd + 1, // l starts from 0 - this->nmaxd, - 2 * this->lmaxd + 1); // m ==> 2*l+1 - this->inl_index[it].create(na[it], this->lmaxd + 1, this->nmaxd); GlobalV::ofs_running << " Type " << it + 1 << " number_of_atoms " << na[it] << std::endl; @@ -167,11 +153,6 @@ void LCAO_Deepks::init_index(const int ntype, { for (int n = 0; n < orb.Alpha[0].getNchi(l); n++) { - for (int m = 0; m < 2 * l + 1; m++) - { - this->alpha_index[it](ia, l, n, m) = alpha; - alpha++; - } this->inl_index[it](ia, l, n) = inl; this->inl_l[inl] = l; inl++; @@ -179,7 +160,6 @@ void LCAO_Deepks::init_index(const int ntype, } } // end ia } // end it - assert(this->n_descriptor == alpha); assert(Total_nchi == inl); GlobalV::ofs_running << " descriptors_per_atom " << this->des_per_atom << std::endl; GlobalV::ofs_running << " total_descriptors " << this->n_descriptor << std::endl; diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h index 5b3041e72d..7341bd60fb 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h @@ -64,43 +64,18 @@ class LCAO_Deepks /// Correction term to Hamiltonian, for multi-k std::vector>> H_V_delta_k; - // functions for hr status: 1. get value; 2. set value; - int get_hr_cal() - { - return this->hr_cal; - } - void set_hr_cal(bool cal) - { - this->hr_cal = cal; - } - - // temporary add two getters for inl_index and gedm - int get_inl(const int& T0, const int& I0, const int& L0, const int& N0) - { - return inl_index[T0](I0, L0, N0); - } - const double* get_gedms(const int& inl) - { - return gedm[inl]; - } - - int get_lmaxd() - { - return lmaxd; - } //------------------- // private variables //------------------- // private: - public: // change to public to reconstuct the code, 2024-07-22 by mohan - int lmaxd = 0; // max l of descirptors - int nmaxd = 0; //#. descriptors per l - int inlmax = 0; // tot. number {i,n,l} - atom, n, l - int n_descriptor; // natoms * des_per_atom, size of descriptor(projector) basis set - int des_per_atom; // \sum_L{Nchi(L)*(2L+1)} - int* inl_l; // inl_l[inl_index] = l of descriptor with inl_index - ModuleBase::IntArray* alpha_index; // seems not used in the code - ModuleBase::IntArray* inl_index; // caoyu add 2021-05-07 + public: // change to public to reconstuct the code, 2024-07-22 by mohan + int lmaxd = 0; // max l of descirptors + int nmaxd = 0; //#. descriptors per l + int inlmax = 0; // tot. number {i,n,l} - atom, n, l + int n_descriptor; // natoms * des_per_atom, size of descriptor(projector) basis set + int des_per_atom; // \sum_L{Nchi(L)*(2L+1)} + int* inl_l; // inl_l[inl_index] = l of descriptor with inl_index + ModuleBase::IntArray* inl_index; // caoyu add 2021-05-07 bool init_pdm = false; // for DeePKS NSCF calculation, set init_pdm to skip the calculation of pdm in SCF iteration @@ -120,14 +95,15 @@ class LCAO_Deepks /// dE/dD, autograd from loaded model(E: Ry) double** gedm; //[tot_Inl][(2l+1)*(2l+1)] - // HR status, - // true : HR should be calculated - // false : HR has been calculated - bool hr_cal = true; - - //------------------- - // subroutines, grouped according to the file they are in: - //------------------- + // functions for hr status: 1. get value; 2. set value; + int get_hr_cal() + { + return this->hr_cal; + } + void set_hr_cal(bool cal) + { + this->hr_cal = cal; + } //------------------- // LCAO_deepks.cpp @@ -164,16 +140,16 @@ class LCAO_Deepks void dpks_cal_e_delta_band(const std::vector>& dm, const int nks); private: + // flag of HR status, + // true : HR should be calculated + // false : HR has been calculated + bool hr_cal = true; + // arrange index of descriptor in all atoms void init_index(const int ntype, const int nat, std::vector na, const int tot_inl, const LCAO_Orbitals& orb); const Parallel_Orbitals* pv; }; -namespace GlobalC -{ -extern LCAO_Deepks ld; -} - #endif #endif diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp index 8215e634ca..58e2d0e349 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp @@ -31,17 +31,25 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, ModuleBase::TITLE("LCAO_Deepks_Interface", "out_deepks_labels"); ModuleBase::timer::tick("LCAO_Deepks_Interface", "out_deepks_labels"); + // Note: out_deepks_labels does not support equivariant version now! + // define TH for different types using TH = std::conditional_t::value, ModuleBase::matrix, ModuleBase::ComplexMatrix>; // These variables are frequently used in the following code - const int inlmax = ld->inlmax; - const int lmaxd = ld->lmaxd; + const int inlmax = orb.Alpha[0].getTotal_nchi() * nat; + const int lmaxd = orb.get_lmax_d(); + const int des_per_atom = ld->des_per_atom; const int* inl_l = ld->inl_l; const ModuleBase::IntArray* inl_index = ld->inl_index; const std::vector*> phialpha = ld->phialpha; + std::vector pdm = ld->pdm; + bool init_pdm = ld->init_pdm; + double E_delta = ld->E_delta; + double e_delta_band = ld->e_delta_band; + double** gedm = ld->gedm; const int my_rank = GlobalV::MY_RANK; const int nspin = PARAM.inp.nspin; @@ -76,7 +84,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, if (PARAM.inp.deepks_scf) { /// ebase :no deepks E_delta including - LCAO_deepks_io::save_npy_e(etot - ld->E_delta, file_ebase, my_rank); + LCAO_deepks_io::save_npy_e(etot - E_delta, file_ebase, my_rank); } else // deepks_scf = 0; base calculation { @@ -136,15 +144,15 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, if (PARAM.inp.deepks_bandgap) { const int nocc = (PARAM.inp.nelec + 1) / 2; - std::vector o_tot(nks); + ModuleBase::matrix o_tot(nks, 1); for (int iks = 0; iks < nks; ++iks) { // record band gap for each k point (including spin) - o_tot[iks] = ekb(iks, nocc) - ekb(iks, nocc - 1); + o_tot(iks, 0) = ekb(iks, nocc) - ekb(iks, nocc - 1); } const std::string file_otot = PARAM.globalv.global_out_dir + "deepks_otot.npy"; - LCAO_deepks_io::save_npy_o(o_tot, file_otot, nks, my_rank); + LCAO_deepks_io::save_matrix2npy(file_otot, o_tot, my_rank); // Unit: Ry if (PARAM.inp.deepks_scf) { @@ -177,7 +185,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, elecstate::cal_dm(ParaV, wg_hl, psi, dm_bandgap); } - std::vector o_delta(nks, 0.0); + ModuleBase::matrix o_delta(nks, 1); // calculate and save orbital_precalc: [nks,NAt,NDscrpt] torch::Tensor orbital_precalc; @@ -203,17 +211,12 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, LCAO_deepks_io::save_tensor2npy(file_orbpre, orbital_precalc, my_rank); const std::string file_obase = PARAM.globalv.global_out_dir + "deepks_obase.npy"; - std::vector o_base(nks); - for (int iks = 0; iks < nks; ++iks) - { - o_base[iks] = o_tot[iks] - o_delta[iks]; - } - LCAO_deepks_io::save_npy_o(o_base, file_obase, nks, my_rank); - } // end deepks_scf == 1 - else // deepks_scf == 0 + LCAO_deepks_io::save_matrix2npy(file_obase, o_tot - o_delta, my_rank); // Unit: Ry + } // end deepks_scf == 1 + else // deepks_scf == 0 { const std::string file_obase = PARAM.globalv.global_out_dir + "deepks_obase.npy"; - LCAO_deepks_io::save_npy_o(o_tot, file_obase, nks, my_rank); // no scf, o_tot=o_base + LCAO_deepks_io::save_matrix2npy(file_obase, o_tot, my_rank); // no scf, o_tot=o_base } // end deepks_scf == 0 } // end bandgap label @@ -323,7 +326,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, if (!PARAM.inp.deepks_scf) { DeePKS_domain::cal_pdm< - TK>(ld->init_pdm, inlmax, lmaxd, inl_l, inl_index, dm, phialpha, ucell, orb, GridD, *ParaV, pdm); + TK>(init_pdm, inlmax, lmaxd, inl_l, inl_index, dm, phialpha, ucell, orb, GridD, *ParaV, pdm); } DeePKS_domain::check_pdm(inlmax, inl_l, pdm); // print out the projected dm for NSCF calculaiton @@ -349,22 +352,22 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, /// print out deepks information to the screen if (PARAM.inp.deepks_scf) { - DeePKS_domain::cal_e_delta_band(dm->get_DMK_vector(), *h_delta, nks, ParaV, ld->e_delta_band); - std::cout << "E_delta_band = " << std::setprecision(8) << ld->e_delta_band << " Ry" - << " = " << std::setprecision(8) << ld->e_delta_band * ModuleBase::Ry_to_eV << " eV" << std::endl; - std::cout << "E_delta_NN = " << std::setprecision(8) << ld->E_delta << " Ry" - << " = " << std::setprecision(8) << ld->E_delta * ModuleBase::Ry_to_eV << " eV" << std::endl; + DeePKS_domain::cal_e_delta_band(dm->get_DMK_vector(), *h_delta, nks, ParaV, e_delta_band); + std::cout << "E_delta_band = " << std::setprecision(8) << e_delta_band << " Ry" + << " = " << std::setprecision(8) << e_delta_band * ModuleBase::Ry_to_eV << " eV" << std::endl; + std::cout << "E_delta_NN = " << std::setprecision(8) << E_delta << " Ry" + << " = " << std::setprecision(8) << E_delta * ModuleBase::Ry_to_eV << " eV" << std::endl; if (PARAM.inp.deepks_out_unittest) { LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, ParaV->nrow, dm->get_DMK_vector()); - DeePKS_domain::check_gedm(inlmax, inl_l, ld->gedm); + DeePKS_domain::check_gedm(inlmax, inl_l, gedm); std::ofstream ofs("E_delta_bands.dat"); - ofs << std::setprecision(10) << ld->e_delta_band; + ofs << std::setprecision(10) << e_delta_band; std::ofstream ofs1("E_delta.dat"); - ofs1 << std::setprecision(10) << ld->E_delta; + ofs1 << std::setprecision(10) << E_delta; } } ModuleBase::timer::tick("LCAO_Deepks_Interface", "out_deepks_labels"); diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp index 28e5268891..52a1aeacaa 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp @@ -3,6 +3,7 @@ #ifdef __DEEPKS #include "LCAO_deepks_io.h" +#include "module_base/tool_quit.h" #include "npy.hpp" #include @@ -149,65 +150,6 @@ void LCAO_deepks_io::save_npy_e(const double& e, const std::string& e_file, cons return; } -// saves force in numpy format -void LCAO_deepks_io::save_npy_f(const ModuleBase::matrix& f, const std::string& f_file, const int rank) -{ - ModuleBase::TITLE("LCAO_deepks_io", "save_npy_f"); - - if (rank != 0) - { - return; - } - - // save force in unit: Rydberg/Bohr - const long unsigned fshape[] = {static_cast(f.nr), static_cast(f.nc)}; - npy::SaveArrayAsNumpy(f_file, false, 2, fshape, f.c); - return; -} - -// saves stress in numpy format -void LCAO_deepks_io::save_npy_s(const ModuleBase::matrix& stress, - const std::string& s_file, - const double& omega, - const int rank) -{ - ModuleBase::TITLE("LCAO_deepks_io", "save_npy_s"); - if (rank != 0) - { - return; - } - - const long unsigned sshape[] = {6}; - std::vector npy_s; - - for (int ipol = 0; ipol < 3; ++ipol) - { - for (int jpol = ipol; jpol < 3; jpol++) - { - npy_s.push_back(stress(ipol, jpol) * omega); - } - } - npy::SaveArrayAsNumpy(s_file, false, 1, sshape, npy_s); - return; -} - -void LCAO_deepks_io::save_npy_o(const std::vector& bandgap, - const std::string& o_file, - const int nks, - const int rank) -{ - ModuleBase::TITLE("LCAO_deepks_io", "save_npy_o"); - if (rank != 0) - { - return; - } - - // save o_base - const long unsigned oshape[] = {static_cast(nks), 1}; - npy::SaveArrayAsNumpy(o_file, false, 2, oshape, bandgap); - return; -} - template void LCAO_deepks_io::save_npy_h(const std::vector& hamilt, const std::string& h_file, @@ -243,7 +185,8 @@ void LCAO_deepks_io::save_npy_h(const std::vector& hamilt, void LCAO_deepks_io::save_matrix2npy(const std::string& file_name, const ModuleBase::matrix& matrix, const int rank, - const double& scale) + const double& scale, + const char mode) { ModuleBase::TITLE("LCAO_deepks_io", "save_matrix2npy"); @@ -251,16 +194,67 @@ void LCAO_deepks_io::save_matrix2npy(const std::string& file_name, { return; } + const int nr = matrix.nr; + const int nc = matrix.nc; + int size = 0; + std::vector shape; - const unsigned long shape[] = {static_cast(matrix.nr), static_cast(matrix.nc)}; + if (mode == 'U' || mode == 'L') // upper or lower triangular + { + assert(nr == nc); + size = nr * (nr + 1) / 2; + shape.resize(1); + shape[0] = size; + } + else if (mode == 'N') // normal + { + size = nr * nc; + shape.resize(2); + shape[0] = nr; + shape[1] = nc; + } + else + { + ModuleBase::WARNING_QUIT("save_matrix2npy", "Invalid mode! Support only 'U', 'L', 'N'."); + } - std::vector scaled_data(matrix.nr * matrix.nc); - for (int i = 0; i < matrix.nr * matrix.nc; ++i) + std::vector scaled_data(size); + if (mode == 'U') // upper triangular { - scaled_data[i] = matrix.c[i] * scale; + int index = 0; + for (int i = 0; i < nr; ++i) + { + for (int j = i; j < nc; ++j) + { + scaled_data[index] = matrix(i, j) * scale; + index++; + } + } + } + else if (mode == 'L') // lower triangular + { + int index = 0; + for (int i = 0; i < nr; ++i) + { + for (int j = 0; j <= i; ++j) + { + scaled_data[index] = matrix(i, j) * scale; + index++; + } + } + } + else // normal + { + for (int i = 0; i < nr; ++i) + { + for (int j = 0; j < nc; ++j) + { + scaled_data[i * nc + j] = matrix(i, j) * scale; + } + } } - npy::SaveArrayAsNumpy(file_name, false, 2, shape, scaled_data.data()); + npy::SaveArrayAsNumpy(file_name, false, shape.size(), shape.data(), scaled_data); return; } diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h index 3d27c1df75..8ab0725a7f 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h @@ -28,12 +28,9 @@ namespace LCAO_deepks_io /// 1. save_npy_d : descriptor -> deepks_dm_eig.npy /// 2. save_npy_e : energy -/// 3. save_npy_f : force -/// 4. save_npy_s : stress -/// 5. save_npy_o: orbital -/// 6. save_npy_h : Hamiltonian -/// 7. save_matrix2npy : ModuleBase::matrix -> .npy -/// 8. save_tensor2npy : torch::Tensor -> .npy +/// 3. save_npy_h : Hamiltonian +/// 4. save_matrix2npy : ModuleBase::matrix -> .npy, for force, stress and orbital +/// 5. save_tensor2npy : torch::Tensor -> .npy, for precalculation variables /// print density matrices template @@ -56,24 +53,7 @@ void save_npy_e(const double& e, /**<[in] \f$E_{base}\f$ or \f$E_{tot}\f$, in Ry const std::string& e_file, const int rank); -// save force -void save_npy_f(const ModuleBase::matrix& f, /**<[in] \f$F_{base}\f$ or \f$F_{tot}\f$, in Ry/Bohr*/ - const std::string& f_file, - const int rank); - -// save stress -void save_npy_s(const ModuleBase::matrix& stress, /**<[in] \f$S_{base}\f$ or \f$S_{tot}\f$, in Ry/Bohr^3*/ - const std::string& s_file, - const double& omega, - const int rank); - -/// save orbital -void save_npy_o(const std::vector& bandgap, /**<[in] \f$E_{base}\f$ or \f$E_{tot}\f$, in Ry*/ - const std::string& o_file, - const int nks, - const int rank); - -// save Hamiltonian and v_delta_precalc(for deepks_v_delta==1)/phialpha+gevdm(for deepks_v_delta==2) +// save Hamiltonian template void save_npy_h(const std::vector& hamilt, const std::string& h_file, @@ -84,7 +64,8 @@ void save_npy_h(const std::vector& hamilt, void save_matrix2npy(const std::string& file_name, const ModuleBase::matrix& matrix, const int rank, - const double& scale = 1.0); + const double& scale = 1.0, + const char mode = 'N'); template void save_tensor2npy(const std::string& file_name, const torch::Tensor& tensor, const int rank); diff --git a/source/module_hamilt_lcao/module_deepks/deepks_force.cpp b/source/module_hamilt_lcao/module_deepks/deepks_force.cpp index be2f716661..8f7fc86526 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_force.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_force.cpp @@ -15,7 +15,6 @@ void DeePKS_domain::cal_f_delta(const std::vector>& dm, const LCAO_Orbitals& orb, const Grid_Driver& GridD, const Parallel_Orbitals& pv, - const int lmaxd, const int nks, const std::vector>& kvec_d, std::vector*> phialpha, @@ -30,7 +29,7 @@ void DeePKS_domain::cal_f_delta(const std::vector>& dm, f_delta.zero_out(); const double Rcut_Alpha = orb.Alpha[0].getRcut(); - + const int lmaxd = orb.get_lmax_d(); const int nrow = pv.nrow; for (int T0 = 0; T0 < ucell.ntype; T0++) @@ -346,7 +345,6 @@ template void DeePKS_domain::cal_f_delta(const std::vector>& kvec_d, std::vector*> phialpha, @@ -361,7 +359,6 @@ template void DeePKS_domain::cal_f_delta>(const std::vector const LCAO_Orbitals& orb, const Grid_Driver& GridD, const Parallel_Orbitals& pv, - const int lmaxd, const int nks, const std::vector>& kvec_d, std::vector*> phialpha, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_force.h b/source/module_hamilt_lcao/module_deepks/deepks_force.h index f964b03964..e42eaec21f 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_force.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_force.h @@ -1,5 +1,5 @@ -#ifndef DEEPKS_FORCE_H -#define DEEPKS_FORCE_H +#ifndef DEEPKS_FORCE_H +#define DEEPKS_FORCE_H #ifdef __DEEPKS @@ -14,36 +14,34 @@ namespace DeePKS_domain { - //------------------------ - // deepks_force.cpp - //------------------------ +//------------------------ +// deepks_force.cpp +//------------------------ - // This file contains subroutines for calculating F_delta, - // which is defind as sum_mu,nu rho_mu,nu d/dX (V(D)) +// This file contains subroutines for calculating F_delta, +// which is defind as sum_mu,nu rho_mu,nu d/dX (V(D)) - // There are 2 subroutines in this file: - // 1. cal_f_delta, which is used for F_delta calculation - // 2. check_f_delta, which prints F_delta into F_delta.dat for checking +// There are 2 subroutines in this file: +// 1. cal_f_delta, which is used for F_delta calculation +// 2. check_f_delta, which prints F_delta into F_delta.dat for checking template -void cal_f_delta( - const std::vector>& dm, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const Grid_Driver& GridD, - const Parallel_Orbitals& pv, - const int lmaxd, - const int nks, - const std::vector>& kvec_d, - std::vector*> phialpha, - double** gedm, - ModuleBase::IntArray* inl_index, - ModuleBase::matrix& f_delta, - const bool isstress, - ModuleBase::matrix& svnl_dalpha); +void cal_f_delta(const std::vector>& dm, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const Grid_Driver& GridD, + const Parallel_Orbitals& pv, + const int nks, + const std::vector>& kvec_d, + std::vector*> phialpha, + double** gedm, + ModuleBase::IntArray* inl_index, + ModuleBase::matrix& f_delta, + const bool isstress, + ModuleBase::matrix& svnl_dalpha); void check_f_delta(const int nat, ModuleBase::matrix& f_delta, ModuleBase::matrix& svnl_dalpha); -} +} // namespace DeePKS_domain #endif #endif diff --git a/source/module_hamilt_lcao/module_deepks/deepks_orbital.cpp b/source/module_hamilt_lcao/module_deepks/deepks_orbital.cpp index a574092afd..8bbb669a31 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_orbital.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_orbital.cpp @@ -9,7 +9,8 @@ template void DeePKS_domain::cal_o_delta(const std::vector& dm_hl, const std::vector>& h_delta, - std::vector& o_delta, + // std::vector& o_delta, + ModuleBase::matrix& o_delta, const Parallel_Orbitals& pv, const int nks) { @@ -54,11 +55,13 @@ void DeePKS_domain::cal_o_delta(const std::vector& dm_hl, Parallel_Reduce::reduce_all(o_delta_tmp); if constexpr (std::is_same::value) { - o_delta[ik] = o_delta_tmp; + // o_delta[ik] = o_delta_tmp; + o_delta(ik, 0) = o_delta_tmp; } else { - o_delta[ik] = o_delta_tmp.real(); + // o_delta[ik] = o_delta_tmp.real(); + o_delta(ik, 0) = o_delta_tmp.real(); } } return; @@ -66,14 +69,16 @@ void DeePKS_domain::cal_o_delta(const std::vector& dm_hl, template void DeePKS_domain::cal_o_delta(const std::vector& dm_hl, const std::vector>& h_delta, - std::vector& o_delta, + // std::vector& o_delta, + ModuleBase::matrix& o_delta, const Parallel_Orbitals& pv, const int nks); template void DeePKS_domain::cal_o_delta, ModuleBase::ComplexMatrix>( const std::vector& dm_hl, const std::vector>>& h_delta, - std::vector& o_delta, + // std::vector& o_delta, + ModuleBase::matrix& o_delta, const Parallel_Orbitals& pv, const int nks); diff --git a/source/module_hamilt_lcao/module_deepks/deepks_orbital.h b/source/module_hamilt_lcao/module_deepks/deepks_orbital.h index 80e4170018..6c60fd1e96 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_orbital.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_orbital.h @@ -26,7 +26,8 @@ namespace DeePKS_domain template void cal_o_delta(const std::vector& dm_hl, const std::vector>& h_delta, - std::vector& o_delta, + // std::vector& o_delta, + ModuleBase::matrix& o_delta, const Parallel_Orbitals& pv, const int nks); } // namespace DeePKS_domain diff --git a/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.cpp b/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.cpp index 493469072b..b6127c8edc 100644 --- a/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.cpp +++ b/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.cpp @@ -34,14 +34,9 @@ void test_deepks::check_phialpha() { na[it] = ucell.atoms[it].na; } - GlobalC::ld.init(ORB, ucell.nat, ucell.ntype, kv.nkstot, ParaO, na); + this->ld.init(ORB, ucell.nat, ucell.ntype, kv.nkstot, ParaO, na); - DeePKS_domain::allocate_phialpha(PARAM.input.cal_force, - ucell, - ORB, - Test_Deepks::GridD, - &ParaO, - GlobalC::ld.phialpha); + DeePKS_domain::allocate_phialpha(PARAM.input.cal_force, ucell, ORB, Test_Deepks::GridD, &ParaO, this->ld.phialpha); DeePKS_domain::build_phialpha(PARAM.input.cal_force, ucell, @@ -49,9 +44,9 @@ void test_deepks::check_phialpha() Test_Deepks::GridD, &ParaO, overlap_orb_alpha_, - GlobalC::ld.phialpha); + this->ld.phialpha); - DeePKS_domain::check_phialpha(PARAM.input.cal_force, ucell, ORB, Test_Deepks::GridD, &ParaO, GlobalC::ld.phialpha); + DeePKS_domain::check_phialpha(PARAM.input.cal_force, ucell, ORB, Test_Deepks::GridD, &ParaO, this->ld.phialpha); this->compare_with_ref("phialpha.dat", "phialpha_ref.dat"); this->compare_with_ref("dphialpha_x.dat", "dphialpha_x_ref.dat"); @@ -156,38 +151,38 @@ void test_deepks::check_pdm() this->read_dm(); this->set_dm_new(); this->set_p_elec_DM(); - DeePKS_domain::cal_pdm(GlobalC::ld.init_pdm, - GlobalC::ld.inlmax, - GlobalC::ld.lmaxd, - GlobalC::ld.inl_l, - GlobalC::ld.inl_index, + DeePKS_domain::cal_pdm(this->ld.init_pdm, + this->ld.inlmax, + this->ld.lmaxd, + this->ld.inl_l, + this->ld.inl_index, p_elec_DM, - GlobalC::ld.phialpha, + this->ld.phialpha, ucell, ORB, Test_Deepks::GridD, ParaO, - GlobalC::ld.pdm); + this->ld.pdm); } else { this->read_dm_k(kv.nkstot); this->set_dm_k_new(); this->set_p_elec_DM_k(); - DeePKS_domain::cal_pdm(GlobalC::ld.init_pdm, - GlobalC::ld.inlmax, - GlobalC::ld.lmaxd, - GlobalC::ld.inl_l, - GlobalC::ld.inl_index, + DeePKS_domain::cal_pdm(this->ld.init_pdm, + this->ld.inlmax, + this->ld.lmaxd, + this->ld.inl_l, + this->ld.inl_index, p_elec_DM_k, - GlobalC::ld.phialpha, + this->ld.phialpha, ucell, ORB, Test_Deepks::GridD, ParaO, - GlobalC::ld.pdm); + this->ld.pdm); } - DeePKS_domain::check_pdm(GlobalC::ld.inlmax, GlobalC::ld.inl_l, GlobalC::ld.pdm); + DeePKS_domain::check_pdm(this->ld.inlmax, this->ld.inl_l, this->ld.pdm); this->compare_with_ref("pdm.dat", "pdm_ref.dat"); } @@ -195,12 +190,12 @@ void test_deepks::check_gdmx(torch::Tensor& gdmx) { if (PARAM.sys.gamma_only_local) { - DeePKS_domain::cal_gdmx(GlobalC::ld.lmaxd, - GlobalC::ld.inlmax, + DeePKS_domain::cal_gdmx(this->ld.lmaxd, + this->ld.inlmax, kv.nkstot, kv.kvec_d, - GlobalC::ld.phialpha, - GlobalC::ld.inl_index, + this->ld.phialpha, + this->ld.inl_index, dm_new, ucell, ORB, @@ -210,12 +205,12 @@ void test_deepks::check_gdmx(torch::Tensor& gdmx) } else { - DeePKS_domain::cal_gdmx(GlobalC::ld.lmaxd, - GlobalC::ld.inlmax, + DeePKS_domain::cal_gdmx(this->ld.lmaxd, + this->ld.inlmax, kv.nkstot, kv.kvec_d, - GlobalC::ld.phialpha, - GlobalC::ld.inl_index, + this->ld.phialpha, + this->ld.inl_index, dm_k_new, ucell, ORB, @@ -254,12 +249,12 @@ void test_deepks::check_gdmepsl(torch::Tensor& gdmepsl) { if (PARAM.sys.gamma_only_local) { - DeePKS_domain::cal_gdmepsl(GlobalC::ld.lmaxd, - GlobalC::ld.inlmax, + DeePKS_domain::cal_gdmepsl(this->ld.lmaxd, + this->ld.inlmax, kv.nkstot, kv.kvec_d, - GlobalC::ld.phialpha, - GlobalC::ld.inl_index, + this->ld.phialpha, + this->ld.inl_index, dm_new, ucell, ORB, @@ -269,12 +264,12 @@ void test_deepks::check_gdmepsl(torch::Tensor& gdmepsl) } else { - DeePKS_domain::cal_gdmepsl(GlobalC::ld.lmaxd, - GlobalC::ld.inlmax, + DeePKS_domain::cal_gdmepsl(this->ld.lmaxd, + this->ld.inlmax, kv.nkstot, kv.kvec_d, - GlobalC::ld.phialpha, - GlobalC::ld.inl_index, + this->ld.phialpha, + this->ld.inl_index, dm_k_new, ucell, ORB, @@ -299,32 +294,21 @@ void test_deepks::check_gdmepsl(torch::Tensor& gdmepsl) void test_deepks::check_descriptor(std::vector& descriptor) { DeePKS_domain::cal_descriptor(ucell.nat, - GlobalC::ld.inlmax, - GlobalC::ld.inl_l, - GlobalC::ld.pdm, + this->ld.inlmax, + this->ld.inl_l, + this->ld.pdm, descriptor, - GlobalC::ld.des_per_atom); - DeePKS_domain::check_descriptor(GlobalC::ld.inlmax, - GlobalC::ld.des_per_atom, - GlobalC::ld.inl_l, - ucell, - "./", - descriptor); + this->ld.des_per_atom); + DeePKS_domain::check_descriptor(this->ld.inlmax, this->ld.des_per_atom, this->ld.inl_l, ucell, "./", descriptor); this->compare_with_ref("deepks_desc.dat", "descriptor_ref.dat"); } void test_deepks::check_gvx(torch::Tensor& gdmx) { std::vector gevdm; - DeePKS_domain::cal_gevdm(ucell.nat, GlobalC::ld.inlmax, GlobalC::ld.inl_l, GlobalC::ld.pdm, gevdm); + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl_l, this->ld.pdm, gevdm); torch::Tensor gvx; - DeePKS_domain::cal_gvx(ucell.nat, - GlobalC::ld.inlmax, - GlobalC::ld.des_per_atom, - GlobalC::ld.inl_l, - gevdm, - gdmx, - gvx); + DeePKS_domain::cal_gvx(ucell.nat, this->ld.inlmax, this->ld.des_per_atom, this->ld.inl_l, gevdm, gdmx, gvx); DeePKS_domain::check_gvx(gvx); for (int ia = 0; ia < ucell.nat; ia++) @@ -354,12 +338,12 @@ void test_deepks::check_gvx(torch::Tensor& gdmx) void test_deepks::check_gvepsl(torch::Tensor& gdmepsl) { std::vector gevdm; - DeePKS_domain::cal_gevdm(ucell.nat, GlobalC::ld.inlmax, GlobalC::ld.inl_l, GlobalC::ld.pdm, gevdm); + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl_l, this->ld.pdm, gevdm); torch::Tensor gvepsl; DeePKS_domain::cal_gvepsl(ucell.nat, - GlobalC::ld.inlmax, - GlobalC::ld.des_per_atom, - GlobalC::ld.inl_l, + this->ld.inlmax, + this->ld.des_per_atom, + this->ld.inl_l, gevdm, gdmepsl, gvepsl); @@ -379,33 +363,33 @@ void test_deepks::check_gvepsl(torch::Tensor& gdmepsl) void test_deepks::check_edelta(std::vector& descriptor) { - DeePKS_domain::load_model("model.ptg", GlobalC::ld.model_deepks); + DeePKS_domain::load_model("model.ptg", ld.model_deepks); if (PARAM.sys.gamma_only_local) { - GlobalC::ld.allocate_V_delta(ucell.nat, 1); // 1 for gamma-only + ld.allocate_V_delta(ucell.nat, 1); // 1 for gamma-only } else { - GlobalC::ld.allocate_V_delta(ucell.nat, kv.nkstot); + ld.allocate_V_delta(ucell.nat, kv.nkstot); } DeePKS_domain::cal_gedm(ucell.nat, - GlobalC::ld.lmaxd, - GlobalC::ld.nmaxd, - GlobalC::ld.inlmax, - GlobalC::ld.des_per_atom, - GlobalC::ld.inl_l, + this->ld.lmaxd, + this->ld.nmaxd, + this->ld.inlmax, + this->ld.des_per_atom, + this->ld.inl_l, descriptor, - GlobalC::ld.pdm, - GlobalC::ld.model_deepks, - GlobalC::ld.gedm, - GlobalC::ld.E_delta); + this->ld.pdm, + this->ld.model_deepks, + this->ld.gedm, + this->ld.E_delta); std::ofstream ofs("E_delta.dat"); - ofs << std::setprecision(10) << GlobalC::ld.E_delta << std::endl; + ofs << std::setprecision(10) << this->ld.E_delta << std::endl; ofs.close(); this->compare_with_ref("E_delta.dat", "E_delta_ref.dat"); - DeePKS_domain::check_gedm(GlobalC::ld.inlmax, GlobalC::ld.inl_l, GlobalC::ld.gedm); + DeePKS_domain::check_gedm(this->ld.inlmax, this->ld.inl_l, this->ld.gedm); this->compare_with_ref("gedm.dat", "gedm_ref.dat"); } @@ -422,7 +406,8 @@ void test_deepks::cal_H_V_delta() &overlap_orb_alpha_, &ORB, kv.nkstot, - p_elec_DM); + p_elec_DM, + &this->ld); for (int ik = 0; ik < kv.nkstot; ++ik) { op_deepks->init(ik); @@ -443,7 +428,8 @@ void test_deepks::cal_H_V_delta_k() &overlap_orb_alpha_, &ORB, kv.nkstot, - p_elec_DM_k); + p_elec_DM_k, + &this->ld); for (int ik = 0; ik < kv.nkstot; ++ik) { op_deepks->init(ik); @@ -455,16 +441,16 @@ void test_deepks::check_e_deltabands() if (PARAM.sys.gamma_only_local) { this->cal_H_V_delta(); - GlobalC::ld.dpks_cal_e_delta_band(dm_new, 1); + this->ld.dpks_cal_e_delta_band(dm_new, 1); } else { this->cal_H_V_delta_k(); - GlobalC::ld.dpks_cal_e_delta_band(dm_k_new, kv.nkstot); + this->ld.dpks_cal_e_delta_band(dm_k_new, kv.nkstot); } std::ofstream ofs("E_delta_bands.dat"); - ofs << std::setprecision(10) << GlobalC::ld.e_delta_band << std::endl; + ofs << std::setprecision(10) << this->ld.e_delta_band << std::endl; ofs.close(); this->compare_with_ref("E_delta_bands.dat", "E_delta_bands_ref.dat"); } @@ -485,12 +471,11 @@ void test_deepks::check_f_delta_and_stress_delta() ORB, Test_Deepks::GridD, ParaO, - GlobalC::ld.lmaxd, nks, kv.kvec_d, - GlobalC::ld.phialpha, - GlobalC::ld.gedm, - GlobalC::ld.inl_index, + this->ld.phialpha, + this->ld.gedm, + this->ld.inl_index, fvnl_dalpha, cal_stress, svnl_dalpha); @@ -503,12 +488,11 @@ void test_deepks::check_f_delta_and_stress_delta() ORB, Test_Deepks::GridD, ParaO, - GlobalC::ld.lmaxd, nks, kv.kvec_d, - GlobalC::ld.phialpha, - GlobalC::ld.gedm, - GlobalC::ld.inl_index, + this->ld.phialpha, + this->ld.gedm, + this->ld.inl_index, fvnl_dalpha, cal_stress, svnl_dalpha); diff --git a/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.h b/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.h index 1ba554b07d..15afc710d0 100644 --- a/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.h +++ b/source/module_hamilt_lcao/module_deepks/test/LCAO_deepks_test.h @@ -21,11 +21,6 @@ namespace Test_Deepks extern Grid_Driver GridD; } -namespace GlobalC -{ -extern LCAO_Deepks ld; -} - class test_deepks { @@ -43,7 +38,7 @@ class test_deepks Parallel_Orbitals ParaO; Test_Deepks::K_Vectors kv; - // LCAO_Deepks ld; + LCAO_Deepks ld; int failed_check = 0; int total_check = 0;