diff --git a/source/Makefile.Objects b/source/Makefile.Objects index ae6a7bf00b..20315a62cc 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -190,7 +190,7 @@ OBJS_CELL=atom_pseudo.o\ OBJS_DEEPKS=LCAO_deepks.o\ deepks_force.o\ - LCAO_deepks_odelta.o\ + deepks_orbital.o\ LCAO_deepks_io.o\ LCAO_deepks_mpi.o\ LCAO_deepks_pdm.o\ diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h index 9a4f74b68d..68a6359fd3 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h @@ -63,6 +63,7 @@ class Force_LCAO ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, #ifdef __DEEPKS + ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, #endif typename TGint::type& gint, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index cadf796b2c..f19d980d85 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -80,6 +80,9 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, ModuleBase::matrix fewalds; ModuleBase::matrix fcc; ModuleBase::matrix fscc; +#ifdef __DEEPKS + ModuleBase::matrix fvnl_dalpha; // deepks +#endif fvl_dphi.create(nat, 3); // must do it now, update it later, noted by zhengdy @@ -93,6 +96,9 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, fewalds.create(nat, 3); fcc.create(nat, 3); fscc.create(nat, 3); +#ifdef __DEEPKS + fvnl_dalpha.create(nat, 3); // deepks +#endif // calculate basic terms in Force, same method with PW base this->calForcePwPart(ucell, @@ -172,6 +178,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, svnl_dbeta, svl_dphi, #ifdef __DEEPKS + fvnl_dalpha, svnl_dalpha, #endif gint_gamma, @@ -454,7 +461,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // mohan add 2021-08-04 if (PARAM.inp.deepks_scf) { - fcs(iat, i) += GlobalC::ld.F_delta(iat, i); + fcs(iat, i) += fvnl_dalpha(iat, i); } #endif // sum total force for correction @@ -499,7 +506,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, if (PARAM.inp.deepks_scf) { const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; - LCAO_deepks_io::save_npy_f(fcs - GlobalC::ld.F_delta, + LCAO_deepks_io::save_npy_f(fcs - fvnl_dalpha, file_fbase, ucell.nat, GlobalV::MY_RANK); // Ry/Bohr, F_base @@ -636,8 +643,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // caoyu add 2021-06-03 if (PARAM.inp.deepks_scf) { - ModuleIO::print_force(GlobalV::ofs_running, ucell, "DeePKS FORCE", GlobalC::ld.F_delta, true); - // this->print_force("DeePKS FORCE", GlobalC::ld.F_delta, 1, ry); + ModuleIO::print_force(GlobalV::ofs_running, ucell, "DeePKS FORCE", fvnl_dalpha, true); } #endif } @@ -891,6 +897,7 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, #if __DEEPKS + ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, #endif Gint_Gamma& gint_gamma, // mohan add 2024-04-01 @@ -917,6 +924,7 @@ void Force_Stress_LCAO::integral_part(const bool isGammaOnly, svnl_dbeta, svl_dphi, #if __DEEPKS + fvnl_dalpha, svnl_dalpha, #endif gint_gamma, @@ -944,6 +952,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, #if __DEEPKS + ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, #endif Gint_Gamma& gint_gamma, @@ -969,6 +978,7 @@ void Force_Stress_LCAO>::integral_part(const bool isGammaOn svnl_dbeta, svl_dphi, #if __DEEPKS + fvnl_dalpha, svnl_dalpha, #endif gint_k, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h index bc728209d6..89a9bf6d55 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h @@ -96,6 +96,7 @@ class Force_Stress_LCAO ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, #if __DEEPKS + ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, #endif Gint_Gamma& gint_gamma, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp index e54d0b4d3f..83f7987dd7 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp @@ -188,6 +188,7 @@ void Force_LCAO::ftable(const bool isforce, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, #ifdef __DEEPKS + ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, #endif TGint::type& gint, @@ -246,15 +247,13 @@ void Force_LCAO::ftable(const bool isforce, false /*reset dm to gint*/); #ifdef __DEEPKS + const std::vector>& dm_gamma = dm->get_DMK_vector(); if (PARAM.inp.deepks_scf) { - const std::vector>& dm_gamma = 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 // GlobalC::ld.cal_projected_DM(dm, ucell, orb, gd); GlobalC::ld.cal_descriptor(ucell.nat); - GlobalC::ld.cal_gedm(ucell.nat); const int nks = 1; @@ -269,40 +268,9 @@ void Force_LCAO::ftable(const bool isforce, GlobalC::ld.phialpha, GlobalC::ld.gedm, GlobalC::ld.inl_index, - GlobalC::ld.F_delta, + fvnl_dalpha, isstress, svnl_dalpha); - -#ifdef __MPI - Parallel_Reduce::reduce_all(GlobalC::ld.F_delta.c, GlobalC::ld.F_delta.nr * GlobalC::ld.F_delta.nc); - - if (isstress) - { - Parallel_Reduce::reduce_pool(svnl_dalpha.c, svnl_dalpha.nr * svnl_dalpha.nc); - } -#endif - - if (PARAM.inp.deepks_out_unittest) - { - const int nks = 1; // 1 for gamma-only - LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, this->ParaV->nrow, dm_gamma); - - GlobalC::ld.check_projected_dm(); - - GlobalC::ld.check_descriptor(ucell, PARAM.globalv.global_out_dir); - - GlobalC::ld.check_gedm(); - - GlobalC::ld.cal_e_delta_band(dm_gamma, nks); - - std::ofstream ofs("E_delta_bands.dat"); - ofs << std::setprecision(10) << GlobalC::ld.e_delta_band; - - std::ofstream ofs1("E_delta.dat"); - ofs1 << std::setprecision(10) << GlobalC::ld.E_delta; - - DeePKS_domain::check_f_delta(ucell.nat, GlobalC::ld.F_delta, svnl_dalpha); - } } #endif @@ -312,6 +280,9 @@ void Force_LCAO::ftable(const bool isforce, Parallel_Reduce::reduce_pool(ftvnl_dphi.c, ftvnl_dphi.nr * ftvnl_dphi.nc); Parallel_Reduce::reduce_pool(fvnl_dbeta.c, fvnl_dbeta.nr * fvnl_dbeta.nc); Parallel_Reduce::reduce_pool(fvl_dphi.c, fvl_dphi.nr * fvl_dphi.nc); +#ifdef __DEEPKS + Parallel_Reduce::reduce_pool(fvnl_dalpha.c, fvnl_dalpha.nr * fvnl_dalpha.nc); +#endif } if (isstress) { @@ -319,7 +290,36 @@ void Force_LCAO::ftable(const bool isforce, Parallel_Reduce::reduce_pool(stvnl_dphi.c, stvnl_dphi.nr * stvnl_dphi.nc); Parallel_Reduce::reduce_pool(svnl_dbeta.c, svnl_dbeta.nr * svnl_dbeta.nc); Parallel_Reduce::reduce_pool(svl_dphi.c, svl_dphi.nr * svl_dphi.nc); +#ifdef __DEEPKS + Parallel_Reduce::reduce_pool(svnl_dalpha.c, svnl_dalpha.nr * svnl_dalpha.nc); +#endif + } + +#ifdef __DEEPKS + // It seems these test should not all be here, should be moved in the future + // Also, these test are not in multi-k case now + if (PARAM.inp.deepks_scf && PARAM.inp.deepks_out_unittest) + { + const int nks = 1; // 1 for gamma-only + LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, this->ParaV->nrow, dm_gamma); + + GlobalC::ld.check_projected_dm(); + + GlobalC::ld.check_descriptor(ucell, PARAM.globalv.global_out_dir); + + GlobalC::ld.check_gedm(); + + GlobalC::ld.cal_e_delta_band(dm_gamma, nks); + + std::ofstream ofs("E_delta_bands.dat"); + ofs << std::setprecision(10) << GlobalC::ld.e_delta_band; + + std::ofstream ofs1("E_delta.dat"); + ofs1 << std::setprecision(10) << GlobalC::ld.E_delta; + + DeePKS_domain::check_f_delta(ucell.nat, fvnl_dalpha, svnl_dalpha); } +#endif // delete DSloc_x, DSloc_y, DSloc_z // delete DHloc_fixed_x, DHloc_fixed_y, DHloc_fixed_z diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp index 2b91e7dae2..7739a02338 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp @@ -282,6 +282,7 @@ void Force_LCAO>::ftable(const bool isforce, ModuleBase::matrix& svnl_dbeta, ModuleBase::matrix& svl_dphi, #ifdef __DEEPKS + ModuleBase::matrix& fvnl_dalpha, ModuleBase::matrix& svnl_dalpha, #endif TGint>::type& gint, @@ -363,17 +364,9 @@ void Force_LCAO>::ftable(const bool isforce, GlobalC::ld.phialpha, GlobalC::ld.gedm, GlobalC::ld.inl_index, - GlobalC::ld.F_delta, + fvnl_dalpha, isstress, svnl_dalpha); - -#ifdef __MPI - Parallel_Reduce::reduce_all(GlobalC::ld.F_delta.c, GlobalC::ld.F_delta.nr * GlobalC::ld.F_delta.nc); - if (isstress) - { - Parallel_Reduce::reduce_pool(svnl_dalpha.c, svnl_dalpha.nr * svnl_dalpha.nc); - } -#endif } #endif @@ -386,6 +379,9 @@ void Force_LCAO>::ftable(const bool isforce, Parallel_Reduce::reduce_pool(ftvnl_dphi.c, ftvnl_dphi.nr * ftvnl_dphi.nc); Parallel_Reduce::reduce_pool(fvnl_dbeta.c, fvnl_dbeta.nr * fvnl_dbeta.nc); Parallel_Reduce::reduce_pool(fvl_dphi.c, fvl_dphi.nr * fvl_dphi.nc); +#ifdef __DEEPKS + Parallel_Reduce::reduce_pool(fvnl_dalpha.c, fvnl_dalpha.nr * fvnl_dalpha.nc); +#endif } if (isstress) { @@ -393,6 +389,9 @@ void Force_LCAO>::ftable(const bool isforce, Parallel_Reduce::reduce_pool(stvnl_dphi.c, stvnl_dphi.nr * stvnl_dphi.nc); Parallel_Reduce::reduce_pool(svnl_dbeta.c, svnl_dbeta.nr * svnl_dbeta.nc); Parallel_Reduce::reduce_pool(svl_dphi.c, svl_dphi.nr * svl_dphi.nc); +#ifdef __DEEPKS + Parallel_Reduce::reduce_pool(svnl_dalpha.c, svnl_dalpha.nr * svnl_dalpha.nc); +#endif } ModuleBase::timer::tick("Force_LCAO", "ftable"); diff --git a/source/module_hamilt_lcao/module_deepks/CMakeLists.txt b/source/module_hamilt_lcao/module_deepks/CMakeLists.txt index 8fb8ad86ab..1a2d618cbf 100644 --- a/source/module_hamilt_lcao/module_deepks/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_deepks/CMakeLists.txt @@ -2,7 +2,7 @@ if(ENABLE_DEEPKS) list(APPEND objects LCAO_deepks.cpp deepks_force.cpp - LCAO_deepks_odelta.cpp + deepks_orbital.cpp LCAO_deepks_io.cpp LCAO_deepks_mpi.cpp LCAO_deepks_pdm.cpp diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp index 43d4a68784..ad64939097 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.cpp @@ -332,8 +332,6 @@ void LCAO_Deepks::allocate_V_delta(const int nat, const int nks) } if (PARAM.inp.cal_force) { - // init F_delta - F_delta.create(nat, 3); if (PARAM.inp.deepks_out_labels) { this->init_gdmx(nat); @@ -342,34 +340,24 @@ void LCAO_Deepks::allocate_V_delta(const int nat, const int nks) // gdmx is used only in calculating gvx } - if (PARAM.inp.deepks_bandgap) - { - // init o_delta - o_delta.create(nks, 1); - } - return; } void LCAO_Deepks::init_orbital_pdm_shell(const int nks) { - this->orbital_pdm_shell = new double***[nks]; + this->orbital_pdm_shell = new double**[nks]; for (int iks = 0; iks < nks; iks++) { - this->orbital_pdm_shell[iks] = new double**[1]; - for (int hl = 0; hl < 1; hl++) + this->orbital_pdm_shell[iks] = new double*[this->inlmax]; + for (int inl = 0; inl < this->inlmax; inl++) { - this->orbital_pdm_shell[iks][hl] = new double*[this->inlmax]; - - for (int inl = 0; inl < this->inlmax; inl++) - { - this->orbital_pdm_shell[iks][hl][inl] = new double[(2 * this->lmaxd + 1) * (2 * this->lmaxd + 1)]; - ModuleBase::GlobalFunc::ZEROS(orbital_pdm_shell[iks][hl][inl], - (2 * this->lmaxd + 1) * (2 * this->lmaxd + 1)); - } + this->orbital_pdm_shell[iks][inl] = new double[(2 * this->lmaxd + 1) * (2 * this->lmaxd + 1)]; + ModuleBase::GlobalFunc::ZEROS(orbital_pdm_shell[iks][inl], + (2 * this->lmaxd + 1) * (2 * this->lmaxd + 1)); } + } return; @@ -379,13 +367,9 @@ void LCAO_Deepks::del_orbital_pdm_shell(const int nks) { for (int iks = 0; iks < nks; iks++) { - for (int hl = 0; hl < 1; hl++) + for (int inl = 0; inl < this->inlmax; inl++) { - for (int inl = 0; inl < this->inlmax; inl++) - { - delete[] this->orbital_pdm_shell[iks][hl][inl]; - } - delete[] this->orbital_pdm_shell[iks][hl]; + delete[] this->orbital_pdm_shell[iks][inl]; } delete[] this->orbital_pdm_shell[iks]; } diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h index 74400976e4..8737515004 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks.h @@ -5,6 +5,7 @@ #include "deepks_force.h" #include "deepks_hmat.h" +#include "deepks_orbital.h" #include "module_base/complexmatrix.h" #include "module_base/intarray.h" #include "module_base/matrix.h" @@ -50,20 +51,12 @@ class LCAO_Deepks ///(Unit: Ry) \f$tr(\rho H_\delta), \rho = \sum_i{c_{i, \mu}c_{i,\nu}} \f$ (for gamma_only) double e_delta_band = 0.0; - ///(Unit: Ry) \f$tr(\rho_{HL} H_\delta), - ///\rho_{HL} = c_{L, \mu}c_{L,\nu} - c_{H, \mu}c_{H,\nu} \f$ (for gamma_only) - ModuleBase::matrix o_delta; - /// Correction term to the Hamiltonian matrix: \f$\langle\phi|V_\delta|\phi\rangle\f$ (for gamma only) /// The size of first dimension is 1, which is used for the consitence with H_V_delta_k std::vector> H_V_delta; /// Correction term to Hamiltonian, for multi-k std::vector>> H_V_delta_k; - // F_delta will be deleted soon, mohan 2024-07-25 - ///(Unit: Ry/Bohr) Total Force due to the DeePKS correction term \f$E_{\delta}\f$ - ModuleBase::matrix F_delta; - // k index of HOMO for multi-k bandgap label. QO added 2022-01-24 int h_ind = 0; @@ -151,8 +144,8 @@ class LCAO_Deepks // dD/dX, tensor form of gdmx std::vector gdmr_vector; - // orbital_pdm_shell:[1,Inl,nm*nm]; \langle \phi_\mu|\alpha\rangle\langle\alpha|\phi_\nu\ranlge - double**** orbital_pdm_shell; + // orbital_pdm_shell:[Inl,nm*nm]; \langle \phi_\mu|\alpha\rangle\langle\alpha|\phi_\nu\ranlge + double*** orbital_pdm_shell; // orbital_precalc:[1,NAt,NDscrpt]; gvdm*orbital_pdm_shell torch::Tensor orbital_precalc_tensor; @@ -359,19 +352,7 @@ class LCAO_Deepks template void dpks_cal_e_delta_band(const std::vector>& dm, const int nks); - //------------------- - // LCAO_deepks_odelta.cpp - //------------------- - - // This file contains subroutines for calculating O_delta, - // which corresponds to the correction of the band gap. - public: - template - void cal_o_delta( - const std::vector>& dm_hl /**<[in] modified density matrix that contains HOMO and LUMO only*/, - const int nks); - //------------------- // LCAO_deepks_torch.cpp //------------------- @@ -401,7 +382,7 @@ class LCAO_Deepks // 9. check_gedm : prints gedm for checking // 10. cal_orbital_precalc : orbital_precalc is usted for training with orbital label, // which equals gvdm * orbital_pdm_shell, - // orbital_pdm_shell[1,Inl,nm*nm] = dm_hl * overlap * overlap + // orbital_pdm_shell[Inl,nm*nm] = dm_hl * overlap * overlap // 11. cal_v_delta_precalc : v_delta_precalc is used for training with v_delta label, // which equals gvdm * v_delta_pdm_shell, // v_delta_pdm_shell = overlap * overlap @@ -443,12 +424,18 @@ class LCAO_Deepks // calculates orbital_precalc template - void cal_orbital_precalc(const std::vector>& dm_hl /**<[in] density matrix*/, + void cal_orbital_precalc(const std::vector& dm_hl, + const int lmaxd, + const int inlmax, const int nat, const int nks, + const int* inl_l, const std::vector>& kvec_d, + const std::vector*> phialpha, + const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, + const Parallel_Orbitals& pv, const Grid_Driver& GridD); // calculates v_delta_precalc 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 3d71c04d88..37a5a2cc1b 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp @@ -1,12 +1,12 @@ #ifdef __DEEPKS #include "LCAO_deepks_interface.h" -#include "module_parameter/parameter.h" -#include "LCAO_deepks_io.h" // mohan add 2024-07-22 +#include "LCAO_deepks_io.h" // mohan add 2024-07-22 #include "module_base/global_variable.h" #include "module_base/tool_title.h" #include "module_elecstate/cal_dm.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_parameter/parameter.h" template LCAO_Deepks_Interface::LCAO_Deepks_Interface(std::shared_ptr ld_in) : ld(ld_in) @@ -49,8 +49,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 - ld->E_delta, file_ebase, my_rank); } else // deepks_scf = 0; base calculation { @@ -58,83 +57,86 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, LCAO_deepks_io::save_npy_e(etot, file_ebase, my_rank); } + std::vector>* h_delta = nullptr; + if constexpr (std::is_same::value) + { + h_delta = &ld->H_V_delta; + } + else + { + h_delta = &ld->H_V_delta_k; + } + if (PARAM.inp.deepks_bandgap) { const int nocc = PARAM.inp.nelec / 2; - ModuleBase::matrix deepks_bands; - deepks_bands.create(nks, 1); + std::vector o_tot(nks); for (int iks = 0; iks < nks; ++iks) { // record band gap for each k point (including spin) - for (int hl = 0; hl < 1; ++hl) - { - deepks_bands(iks, hl) = ekb(iks, nocc + hl) - ekb(iks, nocc - 1 + hl); - } + o_tot[iks] = 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(deepks_bands, file_otot, nks, my_rank); + LCAO_deepks_io::save_npy_o(o_tot, file_otot, nks, my_rank); if (PARAM.inp.deepks_scf) { ModuleBase::matrix wg_hl; - std::vector> dm_bandgap; + std::vector dm_bandgap; // Calculate O_delta if constexpr (std::is_same::value) // for gamma only { wg_hl.create(nspin, PARAM.inp.nbands); dm_bandgap.resize(nspin); - for (int is = 0; is < nspin; ++is) { - for (int ib = 0; ib < 1; ++ib) - { - wg_hl.zero_out(); - wg_hl(is, ib + nocc - 1) = -1.0; - wg_hl(is, ib + nocc) = 1.0; - dm_bandgap[ib].resize(nspin); - elecstate::cal_dm(ParaV, wg_hl, psi, dm_bandgap[ib]); - } + wg_hl.zero_out(); + wg_hl(is, nocc - 1) = -1.0; + wg_hl(is, nocc) = 1.0; + elecstate::cal_dm(ParaV, wg_hl, psi, dm_bandgap); } } else // for multi-k { wg_hl.create(nks, PARAM.inp.nbands); - dm_bandgap.resize(1); - - for (int ib = 0; ib < 1; ib++) + dm_bandgap.resize(nks); + wg_hl.zero_out(); + for (int ik = 0; ik < nks; ik++) { - wg_hl.zero_out(); - for (int ik = 0; ik < nks; ik++) - { - wg_hl(ik, ib + nocc - 1) = -1.0; - wg_hl(ik, ib + nocc) = 1.0; - } - dm_bandgap[ib].resize(nks); - elecstate::cal_dm(ParaV, wg_hl, psi, dm_bandgap[ib]); + wg_hl(ik, nocc - 1) = -1.0; + wg_hl(ik, nocc) = 1.0; } + elecstate::cal_dm(ParaV, wg_hl, psi, dm_bandgap); } - - ld->cal_orbital_precalc(dm_bandgap, nat, nks, kvec_d, ucell, orb, GridD); - ld->cal_o_delta(dm_bandgap, nks); + + std::vector o_delta(nks, 0.0); + + ld->cal_orbital_precalc(dm_bandgap, ld->lmaxd, ld->inlmax, nat, nks, ld->inl_l, kvec_d, ld->phialpha, ld->inl_index, ucell, orb, *ParaV, GridD); + DeePKS_domain::cal_o_delta(dm_bandgap, *h_delta, o_delta, *ParaV, nks); // save obase and orbital_precalc - LCAO_deepks_io::save_npy_orbital_precalc(nat, - nks, - ld->des_per_atom, - ld->orbital_precalc_tensor, - PARAM.globalv.global_out_dir, - my_rank); + LCAO_deepks_io::save_npy_orbital_precalc(nat, + nks, + ld->des_per_atom, + ld->orbital_precalc_tensor, + PARAM.globalv.global_out_dir, + my_rank); const std::string file_obase = PARAM.globalv.global_out_dir + "deepks_obase.npy"; - LCAO_deepks_io::save_npy_o(deepks_bands - ld->o_delta, file_obase, nks, my_rank); - } // end deepks_scf == 1 - else // deepks_scf == 0 + 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 { const std::string file_obase = PARAM.globalv.global_out_dir + "deepks_obase.npy"; - LCAO_deepks_io::save_npy_o(deepks_bands, file_obase, nks, my_rank); // no scf, o_tot=o_base - } // end deepks_scf == 0 - } // end bandgap label + LCAO_deepks_io::save_npy_o(o_tot, file_obase, nks, my_rank); // no scf, o_tot=o_base + } // end deepks_scf == 0 + } // end bandgap label // save H(R) matrix if (true) // should be modified later! @@ -145,7 +147,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, // How to save H(R)? } - if(PARAM.inp.deepks_v_delta) + if (PARAM.inp.deepks_v_delta) { std::vector h_tot(nks); std::vector> h_mat(nks, std::vector(ParaV->nloc)); @@ -160,12 +162,12 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, } } - DeePKS_domain::collect_h_mat(*ParaV, h_mat, h_tot, nlocal, nks); + DeePKS_domain::collect_h_mat(*ParaV, h_mat, h_tot, nlocal, nks); const std::string file_htot = PARAM.globalv.global_out_dir + "deepks_htot.npy"; - LCAO_deepks_io::save_npy_h(h_tot, file_htot, nlocal, nks, my_rank); + LCAO_deepks_io::save_npy_h(h_tot, file_htot, nlocal, nks, my_rank); - if(PARAM.inp.deepks_scf) + if (PARAM.inp.deepks_scf) { std::vector v_delta(nks); std::vector h_base(nks); @@ -174,16 +176,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, v_delta[ik].create(nlocal, nlocal); h_base[ik].create(nlocal, nlocal); } - std::vector>* H_V_delta = nullptr; - if constexpr (std::is_same::value) - { - H_V_delta = &ld->H_V_delta; - } - else - { - H_V_delta = &ld->H_V_delta_k; - } - DeePKS_domain::collect_h_mat(*ParaV, *H_V_delta,v_delta,nlocal,nks); + DeePKS_domain::collect_h_mat(*ParaV, *h_delta, v_delta, nlocal, nks); // save v_delta and h_base const std::string file_hbase = PARAM.globalv.global_out_dir + "deepks_hbase.npy"; @@ -191,29 +184,29 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, { h_base[ik] = h_tot[ik] - v_delta[ik]; } - LCAO_deepks_io::save_npy_h(h_base, file_hbase, nlocal, nks, my_rank); + LCAO_deepks_io::save_npy_h(h_base, file_hbase, nlocal, nks, my_rank); const std::string file_vdelta = PARAM.globalv.global_out_dir + "deepks_vdelta.npy"; - LCAO_deepks_io::save_npy_h(v_delta, file_vdelta, nlocal, nks, my_rank); + LCAO_deepks_io::save_npy_h(v_delta, file_vdelta, nlocal, nks, my_rank); - if(PARAM.inp.deepks_v_delta==1)//v_delta_precalc storage method 1 + if (PARAM.inp.deepks_v_delta == 1) // v_delta_precalc storage method 1 { ld->cal_v_delta_precalc(nlocal, nat, nks, kvec_d, ucell, orb, GridD); - LCAO_deepks_io::save_npy_v_delta_precalc(nat, - nks, + LCAO_deepks_io::save_npy_v_delta_precalc(nat, + nks, nlocal, ld->des_per_atom, ld->v_delta_precalc_tensor, PARAM.globalv.global_out_dir, my_rank); } - else if(PARAM.inp.deepks_v_delta==2)//v_delta_precalc storage method 2 + else if (PARAM.inp.deepks_v_delta == 2) // v_delta_precalc storage method 2 { ld->prepare_phialpha(nlocal, nat, nks, kvec_d, ucell, orb, GridD); - LCAO_deepks_io::save_npy_phialpha(nat, - nks, + LCAO_deepks_io::save_npy_phialpha(nat, + nks, nlocal, ld->inlmax, ld->lmaxd, @@ -231,15 +224,14 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, my_rank); } } - else //deepks_scf == 0 + else // deepks_scf == 0 { const std::string file_hbase = PARAM.globalv.global_out_dir + "deepks_hbase.npy"; - LCAO_deepks_io::save_npy_h(h_tot, file_hbase, nlocal, nks, my_rank); + LCAO_deepks_io::save_npy_h(h_tot, file_hbase, nlocal, nks, my_rank); } - }//end v_delta label - - } // end deepks_out_labels + } // end v_delta label + } // end deepks_out_labels // DeePKS PDM and descriptor if (PARAM.inp.deepks_out_labels || PARAM.inp.deepks_scf) @@ -247,38 +239,36 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, // this part is for integrated test of deepks // so it is printed no matter even if deepks_out_labels is not used // when deepks_scf is on, the init pdm should be same as the out pdm, so we should not recalculate the pdm - if(!PARAM.inp.deepks_scf) - { - ld->cal_projected_DM(dm, ucell, orb, GridD); - } + if (!PARAM.inp.deepks_scf) + { + ld->cal_projected_DM(dm, ucell, orb, GridD); + } ld->check_projected_dm(); // print out the projected dm for NSCF calculaiton - ld->cal_descriptor(nat); // final descriptor + ld->cal_descriptor(nat); // final descriptor ld->check_descriptor(ucell, PARAM.globalv.global_out_dir); - if (PARAM.inp.deepks_out_labels) - { - LCAO_deepks_io::save_npy_d( - nat, - ld->des_per_atom, - ld->inlmax, - ld->inl_l, - PARAM.inp.deepks_equiv, - ld->d_tensor, - PARAM.globalv.global_out_dir, - GlobalV::MY_RANK); // libnpy needed - } + if (PARAM.inp.deepks_out_labels) + { + LCAO_deepks_io::save_npy_d(nat, + ld->des_per_atom, + ld->inlmax, + ld->inl_l, + PARAM.inp.deepks_equiv, + ld->d_tensor, + PARAM.globalv.global_out_dir, + GlobalV::MY_RANK); // libnpy needed + } } - + /// print out deepks information to the screen if (PARAM.inp.deepks_scf) { ld->cal_e_delta_band(dm->get_DMK_vector(), nks); 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::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; } 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 2543568a44..5e2eff5860 100644 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp +++ b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp @@ -325,8 +325,8 @@ void LCAO_deepks_io::save_npy_s(const ModuleBase::matrix &stress, } -void LCAO_deepks_io::save_npy_o(const ModuleBase::matrix &bandgap, - const std::string &o_file, +void LCAO_deepks_io::save_npy_o(const std::vector& bandgap, + const std::string& o_file, const int nks, const int rank) { @@ -338,17 +338,7 @@ void LCAO_deepks_io::save_npy_o(const ModuleBase::matrix &bandgap, //save o_base const long unsigned oshape[] = {static_cast(nks), 1 }; - - std::vector npy_o; - for (int iks = 0; iks < nks; ++iks) - { - for (int hl = 0;hl < 1; ++hl) - { - npy_o.push_back(bandgap(iks,hl)); - } - } - - npy::SaveArrayAsNumpy(o_file, false, 2, oshape, npy_o); + npy::SaveArrayAsNumpy(o_file, false, 2, oshape, bandgap); return; } @@ -369,27 +359,23 @@ void LCAO_deepks_io::save_npy_orbital_precalc(const int nat, //save orbital_precalc.npy (when bandgap label is in use) //unit: a.u. const long unsigned gshape[] = {static_cast(nks), - 1, static_cast(nat), static_cast(des_per_atom)}; std::vector npy_orbital_precalc; for (int iks = 0; iks < nks; ++iks) { - for (int hl = 0; hl < 1; ++hl) + for (int iat = 0;iat < nat;++iat) { - for (int iat = 0;iat < nat;++iat) + for(int p=0; p& bandgap, /**<[in] \f$E_{base}\f$ or \f$E_{tot}\f$, in Ry*/ const std::string &o_file, const int nks, const int rank); diff --git a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_odelta.cpp b/source/module_hamilt_lcao/module_deepks/LCAO_deepks_odelta.cpp deleted file mode 100644 index 016e6317b8..0000000000 --- a/source/module_hamilt_lcao/module_deepks/LCAO_deepks_odelta.cpp +++ /dev/null @@ -1,77 +0,0 @@ -//QO 2022-1-7 -//This file contains subroutines for calculating O_delta, i.e., corrections of the bandgap, -#include "module_parameter/parameter.h" -//which is defind as sum_mu,nu rho^{hl}_mu,nu V(D) -//where rho^{hl}_mu,nu = C_{L\mu}C_{L\nu} - C_{H\mu}C_{H\nu}, L for LUMO, H for HOMO - -#ifdef __DEEPKS - -#include "LCAO_deepks.h" -#include "module_base/parallel_reduce.h" - -template -void LCAO_Deepks::cal_o_delta(const std::vector>& dm_hl, const int nks) -{ - ModuleBase::TITLE("LCAO_Deepks", "cal_o_delta"); - - this->o_delta.zero_out(); - for (int ik = 0; ik < nks; ik++) - { - for (int hl = 0; hl < 1; ++hl) - { - TK o_delta_tmp = TK(0.0); - for (int i = 0; i < PARAM.globalv.nlocal; ++i) - { - for (int j = 0; j < PARAM.globalv.nlocal; ++j) - { - const int mu = pv->global2local_row(j); - const int nu = pv->global2local_col(i); - - if (mu >= 0 && nu >= 0) - { - int iic; - if(PARAM.inp.ks_solver=="genelpa" || PARAM.inp.ks_solver=="scalapack_gvx" || PARAM.inp.ks_solver=="pexsi") // save the matrix as column major format - { - iic = mu + nu * pv->nrow; - } - else - { - iic = mu * pv->ncol + nu; - } - if constexpr (std::is_same::value) - { - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - o_delta_tmp += dm_hl[hl][is](nu, mu) * this->H_V_delta[0][iic]; - } - } - else - { - o_delta_tmp += dm_hl[hl][ik](nu, mu) * this->H_V_delta_k[ik][iic]; - } - } - } - } - Parallel_Reduce::reduce_all(o_delta_tmp); - if constexpr (std::is_same::value) - { - this->o_delta(ik,hl) = o_delta_tmp; - } - else - { - this->o_delta(ik,hl) = o_delta_tmp.real(); - } - } - } - return; -} - -template void LCAO_Deepks::cal_o_delta( - const std::vector>& dm_hl, - const int nks); - -template void LCAO_Deepks::cal_o_delta, ModuleBase::ComplexMatrix>( - const std::vector>& dm_hl, - const int nks); - -#endif diff --git a/source/module_hamilt_lcao/module_deepks/deepks_force.cpp b/source/module_hamilt_lcao/module_deepks/deepks_force.cpp index 51992ff5ea..c7eea52036 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_force.cpp +++ b/source/module_hamilt_lcao/module_deepks/deepks_force.cpp @@ -1,11 +1,4 @@ -// wenfei 2022-1-11 -// This file contains subroutines for calculating F_delta, #include "module_parameter/parameter.h" -// 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 calculates F_delta -// 2. check_f_delta, which prints F_delta into F_delta.dat for checking #ifdef __DEEPKS diff --git a/source/module_hamilt_lcao/module_deepks/deepks_force.h b/source/module_hamilt_lcao/module_deepks/deepks_force.h index bd4a183d6a..f964b03964 100644 --- a/source/module_hamilt_lcao/module_deepks/deepks_force.h +++ b/source/module_hamilt_lcao/module_deepks/deepks_force.h @@ -3,7 +3,6 @@ #ifdef __DEEPKS - #include "module_base/complexmatrix.h" #include "module_base/intarray.h" #include "module_base/matrix.h" @@ -13,14 +12,10 @@ #include "module_cell/module_neighbor/sltk_grid_driver.h" #include "module_elecstate/module_dm/density_matrix.h" -#include -#include -#include - namespace DeePKS_domain { //------------------------ - // LCAO_deepks_force.cpp + // deepks_force.cpp //------------------------ // This file contains subroutines for calculating F_delta, diff --git a/source/module_hamilt_lcao/module_deepks/deepks_orbital.cpp b/source/module_hamilt_lcao/module_deepks/deepks_orbital.cpp new file mode 100644 index 0000000000..a574092afd --- /dev/null +++ b/source/module_hamilt_lcao/module_deepks/deepks_orbital.cpp @@ -0,0 +1,80 @@ +#include "module_parameter/parameter.h" + +#ifdef __DEEPKS + +#include "deepks_orbital.h" +#include "module_base/parallel_reduce.h" +#include "module_base/timer.h" + +template +void DeePKS_domain::cal_o_delta(const std::vector& dm_hl, + const std::vector>& h_delta, + std::vector& o_delta, + const Parallel_Orbitals& pv, + const int nks) +{ + ModuleBase::TITLE("DeePKS_domain", "cal_o_delta"); + + for (int ik = 0; ik < nks; ik++) + { + TK o_delta_tmp = TK(0.0); + for (int i = 0; i < PARAM.globalv.nlocal; ++i) + { + for (int j = 0; j < PARAM.globalv.nlocal; ++j) + { + const int mu = pv.global2local_row(j); + const int nu = pv.global2local_col(i); + + if (mu >= 0 && nu >= 0) + { + int iic; + if (PARAM.inp.ks_solver == "genelpa" || PARAM.inp.ks_solver == "scalapack_gvx" + || PARAM.inp.ks_solver == "pexsi") // save the matrix as column major format + { + iic = mu + nu * pv.nrow; + } + else + { + iic = mu * pv.ncol + nu; + } + if constexpr (std::is_same::value) + { + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + o_delta_tmp += dm_hl[is](nu, mu) * h_delta[ik][iic]; + } + } + else + { + o_delta_tmp += dm_hl[ik](nu, mu) * h_delta[ik][iic]; + } + } + } + } + Parallel_Reduce::reduce_all(o_delta_tmp); + if constexpr (std::is_same::value) + { + o_delta[ik] = o_delta_tmp; + } + else + { + o_delta[ik] = o_delta_tmp.real(); + } + } + return; +} + +template void DeePKS_domain::cal_o_delta(const std::vector& dm_hl, + const std::vector>& h_delta, + std::vector& 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, + const Parallel_Orbitals& pv, + const int nks); + +#endif diff --git a/source/module_hamilt_lcao/module_deepks/deepks_orbital.h b/source/module_hamilt_lcao/module_deepks/deepks_orbital.h new file mode 100644 index 0000000000..80e4170018 --- /dev/null +++ b/source/module_hamilt_lcao/module_deepks/deepks_orbital.h @@ -0,0 +1,35 @@ +#ifndef DEEPKS_ORBITAL_H +#define DEEPKS_ORBITAL_H + +#ifdef __DEEPKS + +#include "module_base/complexmatrix.h" +#include "module_base/intarray.h" +#include "module_base/matrix.h" +#include "module_base/timer.h" +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_elecstate/module_dm/density_matrix.h" + +namespace DeePKS_domain +{ +//------------------------ +// deepks_orbital.cpp +//------------------------ + +// This file contains subroutines for calculating O_delta, i.e., corrections of the bandgap, +// which is defind as sum_mu,nu rho^{hl}_mu,nu V(D) +// where rho^{hl}_mu,nu = C_{L\mu}C_{L\nu} - C_{H\mu}C_{H\nu}, L for LUMO, H for HOMO + +// There are 1 subroutines in this file: +// 1. cal_o_delta, which is used for O_delta calculation + +template +void cal_o_delta(const std::vector& dm_hl, + const std::vector>& h_delta, + std::vector& o_delta, + const Parallel_Orbitals& pv, + const int nks); +} // namespace DeePKS_domain + +#endif +#endif diff --git a/source/module_hamilt_lcao/module_deepks/orbital_precalc.cpp b/source/module_hamilt_lcao/module_deepks/orbital_precalc.cpp index 6a55fa78f2..8154746d1b 100644 --- a/source/module_hamilt_lcao/module_deepks/orbital_precalc.cpp +++ b/source/module_hamilt_lcao/module_deepks/orbital_precalc.cpp @@ -2,7 +2,7 @@ /// cal_orbital_precalc : orbital_precalc is usted for training with orbital label, /// which equals gvdm * orbital_pdm_shell, -/// orbital_pdm_shell[1,Inl,nm*nm] = dm_hl * overlap * overlap +/// orbital_pdm_shell[Inl,nm*nm] = dm_hl * overlap * overlap #include "LCAO_deepks.h" #include "LCAO_deepks_io.h" // mohan add 2024-07-22 @@ -14,14 +14,20 @@ #include "module_parameter/parameter.h" // calculates orbital_precalc[1,NAt,NDscrpt] = gvdm * orbital_pdm_shell; -// orbital_pdm_shell[2,Inl,nm*nm] = dm_hl * overlap * overlap; +// orbital_pdm_shell[Inl,nm*nm] = dm_hl * overlap * overlap; template -void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, +void LCAO_Deepks::cal_orbital_precalc(const std::vector& dm_hl, + const int lmaxd, + const int inlmax, const int nat, const int nks, + const int* inl_l, const std::vector>& kvec_d, + const std::vector*> phialpha, + const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, + const Parallel_Orbitals& pv, const Grid_Driver& GridD) { ModuleBase::TITLE("LCAO_Deepks", "cal_orbital_precalc"); @@ -51,7 +57,7 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = this->inl_index[T0](I0, L0, N0); + const int inl = inl_index[T0](I0, L0, N0); const int nm = 2 * L0 + 1; for (int m1 = 0; m1 < nm; ++m1) // m1 = 1 for s, 3 for p, 5 for d @@ -87,13 +93,13 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, if constexpr (std::is_same>::value) { - if (this->phialpha[0]->find_matrix(iat, ibt1, dR1.x, dR1.y, dR1.z) == nullptr) + if (phialpha[0]->find_matrix(iat, ibt1, dR1.x, dR1.y, dR1.z) == nullptr) { continue; } } - auto row_indexes = pv->get_indexes_row(ibt1); + auto row_indexes = pv.get_indexes_row(ibt1); const int row_size = row_indexes.size(); @@ -108,7 +114,7 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, for (int irow = 0; irow < row_size; irow++) { - hamilt::BaseMatrix* row_ptr = this->phialpha[0]->find_matrix(iat, ibt1, dR1); + hamilt::BaseMatrix* row_ptr = phialpha[0]->find_matrix(iat, ibt1, dR1); for (int i = 0; i < trace_alpha_size; i++) { s_1t[i * row_size + irow] = row_ptr->get_value(row_indexes[irow], trace_alpha_row[i]); @@ -143,13 +149,13 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, if constexpr (std::is_same>::value) { - if (this->phialpha[0]->find_matrix(iat, ibt2, dR2.x, dR2.y, dR2.z) == nullptr) + if (phialpha[0]->find_matrix(iat, ibt2, dR2.x, dR2.y, dR2.z) == nullptr) { continue; } } - auto col_indexes = pv->get_indexes_col(ibt2); + auto col_indexes = pv.get_indexes_col(ibt2); const int col_size = col_indexes.size(); if (col_size == 0) { @@ -159,7 +165,7 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, std::vector s_2t(trace_alpha_size * col_size); for (int icol = 0; icol < col_size; icol++) { - hamilt::BaseMatrix* col_ptr = this->phialpha[0]->find_matrix(iat, ibt2, dR2); + hamilt::BaseMatrix* col_ptr = phialpha[0]->find_matrix(iat, ibt2, dR2); for (int i = 0; i < trace_alpha_size; i++) { s_2t[i * col_size + icol] = col_ptr->get_value(col_indexes[icol], trace_alpha_col[i]); @@ -174,17 +180,17 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, { for (int is = 0; is < PARAM.inp.nspin; is++) { - hamilt::AtomPair dm_pair(ibt1, ibt2, 0, 0, 0, pv); + hamilt::AtomPair dm_pair(ibt1, ibt2, 0, 0, 0, &pv); dm_pair.allocate(dm_array.data(), 0); if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver)) { - dm_pair.add_from_matrix(dm_hl[0][is].c, pv->get_row_size(), 1.0, 1); + dm_pair.add_from_matrix(dm_hl[is].c, pv.get_row_size(), 1.0, 1); } else { - dm_pair.add_from_matrix(dm_hl[0][is].c, pv->get_col_size(), 1.0, 0); + dm_pair.add_from_matrix(dm_hl[is].c, pv.get_col_size(), 1.0, 0); } } // is } @@ -197,7 +203,7 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, (dR2 - dR1).x, (dR2 - dR1).y, (dR2 - dR1).z, - pv); + &pv); dm_pair.allocate(&dm_array[ik * row_size * col_size], 0); @@ -212,11 +218,11 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver)) { - dm_pair.add_from_matrix(dm_hl[0][ik].c, pv->get_row_size(), kphase, 1); + dm_pair.add_from_matrix(dm_hl[ik].c, pv.get_row_size(), kphase, 1); } else { - dm_pair.add_from_matrix(dm_hl[0][ik].c, pv->get_col_size(), kphase, 0); + dm_pair.add_from_matrix(dm_hl[ik].c, pv.get_col_size(), kphase, 0); } } // ik } @@ -259,14 +265,14 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = this->inl_index[T0](I0, L0, N0); + const int inl = inl_index[T0](I0, L0, N0); const int nm = 2 * L0 + 1; for (int m1 = 0; m1 < nm; ++m1) // m1 = 1 for s, 3 for p, 5 for d { for (int m2 = 0; m2 < nm; ++m2) // m1 = 1 for s, 3 for p, 5 for d { - orbital_pdm_shell[ik][0][inl][m1 * nm + m2] + orbital_pdm_shell[ik][inl][m1 * nm + m2] += ddot_(&row_size, p_g1dmt + index * row_size * nks, &inc, @@ -285,57 +291,48 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, #ifdef __MPI for (int iks = 0; iks < nks; iks++) { - for (int hl = 0; hl < 1; hl++) + for (int inl = 0; inl < inlmax; inl++) { - for (int inl = 0; inl < this->inlmax; inl++) - { - Parallel_Reduce::reduce_all(this->orbital_pdm_shell[iks][hl][inl], - (2 * this->lmaxd + 1) * (2 * this->lmaxd + 1)); - } + Parallel_Reduce::reduce_all(this->orbital_pdm_shell[iks][inl], + (2 * lmaxd + 1) * (2 * lmaxd + 1)); } } #endif // transfer orbital_pdm_shell to orbital_pdm_shell_vector - int nlmax = this->inlmax / nat; + int nlmax = inlmax / nat; std::vector orbital_pdm_shell_vector; for (int nl = 0; nl < nlmax; ++nl) { - std::vector kiammv; + std::vector kammv; for (int iks = 0; iks < nks; ++iks) { - std::vector iammv; - for (int hl = 0; hl < 1; ++hl) + std::vector ammv; + for (int iat = 0; iat < nat; ++iat) { - std::vector ammv; - for (int iat = 0; iat < nat; ++iat) - { - int inl = iat * nlmax + nl; - int nm = 2 * this->inl_l[inl] + 1; - std::vector mmv; + int inl = iat * nlmax + nl; + int nm = 2 * inl_l[inl] + 1; + std::vector mmv; - for (int m1 = 0; m1 < nm; ++m1) // m1 = 1 for s, 3 for p, 5 for d + for (int m1 = 0; m1 < nm; ++m1) // m1 = 1 for s, 3 for p, 5 for d + { + for (int m2 = 0; m2 < nm; ++m2) // m1 = 1 for s, 3 for p, 5 for d { - for (int m2 = 0; m2 < nm; ++m2) // m1 = 1 for s, 3 for p, 5 for d - { - mmv.push_back(this->orbital_pdm_shell[iks][hl][inl][m1 * nm + m2]); - } + mmv.push_back(this->orbital_pdm_shell[iks][inl][m1 * nm + m2]); } - torch::Tensor mm - = torch::tensor(mmv, torch::TensorOptions().dtype(torch::kFloat64)).reshape({nm, nm}); // nm*nm - - ammv.push_back(mm); } - torch::Tensor amm = torch::stack(ammv, 0); - iammv.push_back(amm); + torch::Tensor mm + = torch::tensor(mmv, torch::TensorOptions().dtype(torch::kFloat64)).reshape({nm, nm}); // nm*nm + + ammv.push_back(mm); } - torch::Tensor iamm = torch::stack(iammv, 0); // inl*nm*nm - kiammv.push_back(iamm); + torch::Tensor amm = torch::stack(ammv, 0); + kammv.push_back(amm); } - torch::Tensor kiamm = torch::stack(kiammv, 0); - orbital_pdm_shell_vector.push_back(kiamm); + torch::Tensor kamm = torch::stack(kammv, 0); + orbital_pdm_shell_vector.push_back(kamm); } assert(orbital_pdm_shell_vector.size() == nlmax); @@ -345,7 +342,7 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, for (int nl = 0; nl < nlmax; ++nl) { orbital_precalc_vector.push_back( - at::einsum("kiamn, avmn->kiav", {orbital_pdm_shell_vector[nl], this->gevdm_vector[nl]})); + at::einsum("kamn, avmn->kav", {orbital_pdm_shell_vector[nl], this->gevdm_vector[nl]})); } this->orbital_precalc_tensor = torch::cat(orbital_precalc_vector, -1); @@ -355,20 +352,32 @@ void LCAO_Deepks::cal_orbital_precalc(const std::vector>& dm_hl, } template void LCAO_Deepks::cal_orbital_precalc( - const std::vector>& dm_hl, + const std::vector& dm_hl, + const int lmaxd, + const int inlmax, const int nat, const int nks, + const int* inl_l, const std::vector>& kvec_d, + const std::vector*> phialpha, + const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, + const Parallel_Orbitals& pv, const Grid_Driver& GridD); template void LCAO_Deepks::cal_orbital_precalc, ModuleBase::ComplexMatrix>( - const std::vector>& dm_hl, + const std::vector& dm_hl, + const int lmaxd, + const int inlmax, const int nat, const int nks, + const int* inl_l, const std::vector>& kvec_d, + const std::vector*> phialpha, + const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, + const Parallel_Orbitals& pv, const Grid_Driver& GridD); #endif