diff --git a/source/source_io/output_mat_sparse.cpp b/source/source_io/output_mat_sparse.cpp index a98e6c8297..9fb5184919 100644 --- a/source/source_io/output_mat_sparse.cpp +++ b/source/source_io/output_mat_sparse.cpp @@ -27,7 +27,7 @@ void output_mat_sparse(const bool& out_mat_hsR, //! generate a file containing the Hamiltonian and S(overlap) matrices if (out_mat_hsR) { - output_HSR(ucell, istep, v_eff, pv, HS_Arrays, grid, kv, *p_dftu, p_ham); + output_HSR(ucell, istep, pv, HS_Arrays, grid, kv, *p_dftu, p_ham); } //! generate a file containing the kinetic energy matrix diff --git a/source/source_io/write_HS_R.cpp b/source/source_io/write_HS_R.cpp index 51ad6c5487..293e315918 100644 --- a/source/source_io/write_HS_R.cpp +++ b/source/source_io/write_HS_R.cpp @@ -15,7 +15,6 @@ template void ModuleIO::output_HSR(const UnitCell& ucell, const int& istep, - const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const Grid_Driver& grid, // mohan add 2024-04-06 @@ -304,7 +303,6 @@ void ModuleIO::output_TR(const int istep, template void ModuleIO::output_HSR( const UnitCell& ucell, const int& istep, - const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const Grid_Driver& grid, @@ -324,7 +322,6 @@ template void ModuleIO::output_HSR( template void ModuleIO::output_HSR>( const UnitCell& ucell, const int& istep, - const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const Grid_Driver& grid, diff --git a/source/source_io/write_HS_R.h b/source/source_io/write_HS_R.h index d92c40e59c..87da54a293 100644 --- a/source/source_io/write_HS_R.h +++ b/source/source_io/write_HS_R.h @@ -15,7 +15,6 @@ using TAC = std::pair>; template void output_HSR(const UnitCell& ucell, const int& istep, - const ModuleBase::matrix& v_eff, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const Grid_Driver& grid, // mohan add 2024-04-06 diff --git a/source/source_lcao/FORCE_gamma.cpp b/source/source_lcao/FORCE_gamma.cpp index fb9249c4b0..4f7bcb89ab 100644 --- a/source/source_lcao/FORCE_gamma.cpp +++ b/source/source_lcao/FORCE_gamma.cpp @@ -231,10 +231,10 @@ void Force_LCAO::ftable(const bool isforce, gd, *this->ParaV, nks, + deepks.ld.deepks_param, kv->kvec_d, deepks.ld.phialpha, deepks.ld.gedm, - deepks.ld.inl_index, fvnl_dalpha, isstress, svnl_dalpha); diff --git a/source/source_lcao/FORCE_k.cpp b/source/source_lcao/FORCE_k.cpp index 7a774481b7..276aab3360 100644 --- a/source/source_lcao/FORCE_k.cpp +++ b/source/source_lcao/FORCE_k.cpp @@ -257,10 +257,10 @@ void Force_LCAO>::ftable(const bool isforce, gd, pv, kv->get_nks(), + deepks.ld.deepks_param, kv->kvec_d, deepks.ld.phialpha, deepks.ld.gedm, - deepks.ld.inl_index, fvnl_dalpha, isstress, svnl_dalpha); diff --git a/source/source_lcao/LCAO_hamilt.hpp b/source/source_lcao/LCAO_hamilt.hpp index 4e91ca724d..69179f64c1 100644 --- a/source/source_lcao/LCAO_hamilt.hpp +++ b/source/source_lcao/LCAO_hamilt.hpp @@ -6,8 +6,8 @@ #include "source_base/abfs-vector3_order.h" #include "source_base/global_variable.h" #include "source_base/timer.h" -#include "source_lcao/spar_exx.h" #include "source_lcao/module_ri/RI_2D_Comm.h" +#include "source_lcao/spar_exx.h" #include #include @@ -21,28 +21,27 @@ // Peize Lin add 2022.09.13 template -void sparse_format::cal_HR_exx(const UnitCell& ucell, +void sparse_format::cal_HR_exx( + const UnitCell& ucell, const Parallel_Orbitals& pv, LCAO_HS_Arrays& HS_Arrays, const int& current_spin, const double& sparse_threshold, const int (&nmp)[3], - const std::vector>, - RI::Tensor>>>& Hexxs) { + const std::vector>, RI::Tensor>>>& Hexxs) +{ ModuleBase::TITLE("sparse_format", "cal_HR_exx"); ModuleBase::timer::tick("sparse_format", "cal_HR_exx"); const Tdata frac = GlobalC::exx_info.info_global.hybrid_alpha; std::map> atoms_pos; - for (int iat = 0; iat < ucell.nat; ++iat) { - atoms_pos[iat] = RI_Util::Vector3_to_array3( - ucell.atoms[ucell.iat2it[iat]] - .tau[ucell.iat2ia[iat]]); + for (int iat = 0; iat < ucell.nat; ++iat) + { + atoms_pos[iat] = RI_Util::Vector3_to_array3(ucell.atoms[ucell.iat2it[iat]].tau[ucell.iat2ia[iat]]); } const std::array, 3> latvec - = {RI_Util::Vector3_to_array3(ucell.a1), // too bad to use GlobalC here, + = {RI_Util::Vector3_to_array3(ucell.a1), // too bad to use GlobalC here, RI_Util::Vector3_to_array3(ucell.a2), RI_Util::Vector3_to_array3(ucell.a3)}; @@ -51,92 +50,83 @@ void sparse_format::cal_HR_exx(const UnitCell& ucell, RI::Cell_Nearest cell_nearest; cell_nearest.init(atoms_pos, latvec, Rs_period); - const std::vector is_list = (PARAM.inp.nspin != 4) - ? std::vector{current_spin} - : std::vector{0, 1, 2, 3}; + const std::vector is_list + = (PARAM.inp.nspin != 4) ? std::vector{current_spin} : std::vector{0, 1, 2, 3}; - for (const int is: is_list) + for (const int is: is_list) { int is0_b = 0; int is1_b = 0; std::tie(is0_b, is1_b) = RI_2D_Comm::split_is_block(is); - if (Hexxs.empty()) + if (Hexxs.empty()) { break; } - for (const auto& HexxA: Hexxs[is]) + for (const auto& HexxA: Hexxs[is]) { const int iat0 = HexxA.first; - for (const auto& HexxB: HexxA.second) + for (const auto& HexxB: HexxA.second) { const int iat1 = HexxB.first.first; const Abfs::Vector3_Order R = RI_Util::array3_to_Vector3( - cell_nearest.get_cell_nearest_discrete(iat0, - iat1, - HexxB.first.second)); + cell_nearest.get_cell_nearest_discrete(iat0, iat1, HexxB.first.second)); HS_Arrays.all_R_coor.insert(R); const RI::Tensor& Hexx = HexxB.second; - for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0) + for (size_t iw0 = 0; iw0 < Hexx.shape[0]; ++iw0) { - const int iwt0 = RI_2D_Comm::get_iwt(ucell,iat0, iw0, is0_b); + const int iwt0 = RI_2D_Comm::get_iwt(ucell, iat0, iw0, is0_b); const int iwt0_local = pv.global2local_row(iwt0); - if (iwt0_local < 0) + if (iwt0_local < 0) { continue; } - for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1) + for (size_t iw1 = 0; iw1 < Hexx.shape[1]; ++iw1) { - const int iwt1 = RI_2D_Comm::get_iwt(ucell,iat1, iw1, is1_b); + const int iwt1 = RI_2D_Comm::get_iwt(ucell, iat1, iw1, is1_b); const int iwt1_local = pv.global2local_col(iwt1); - if (iwt1_local < 0) + if (iwt1_local < 0) { continue; } - if (std::abs(Hexx(iw0, iw1)) > sparse_threshold) + if (std::abs(Hexx(iw0, iw1)) > sparse_threshold) { - if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) + if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) { - auto& HR_sparse_ptr - = HS_Arrays - .HR_sparse[current_spin][R][iwt0]; + auto& HR_sparse_ptr = HS_Arrays.HR_sparse[current_spin][R][iwt0]; double& HR_sparse = HR_sparse_ptr[iwt1]; - HR_sparse += RI::Global_Func::convert( - frac * Hexx(iw0, iw1)); - if (std::abs(HR_sparse) <= sparse_threshold) + HR_sparse += RI::Global_Func::convert(frac * Hexx(iw0, iw1)); + if (std::abs(HR_sparse) <= sparse_threshold) { HR_sparse_ptr.erase(iwt1); } - } - else if (PARAM.inp.nspin == 4) + } + else if (PARAM.inp.nspin == 4) { - auto& HR_sparse_ptr - = HS_Arrays.HR_soc_sparse[R][iwt0]; - - std::complex& HR_sparse - = HR_sparse_ptr[iwt1]; - - HR_sparse += RI::Global_Func::convert< - std::complex>(frac * Hexx(iw0, iw1)); - - if (std::abs(HR_sparse) <= sparse_threshold) + auto& HR_sparse_ptr = HS_Arrays.HR_soc_sparse[R][iwt0]; + + std::complex& HR_sparse = HR_sparse_ptr[iwt1]; + + HR_sparse += RI::Global_Func::convert>(frac * Hexx(iw0, iw1)); + + if (std::abs(HR_sparse) <= sparse_threshold) { HR_sparse_ptr.erase(iwt1); } - } - else + } + else { throw std::invalid_argument(std::string(__FILE__) + " line " - + std::to_string(__LINE__)); + + std::to_string(__LINE__)); } } } diff --git a/source/source_lcao/module_deepks/LCAO_deepks.cpp b/source/source_lcao/module_deepks/LCAO_deepks.cpp index 6755c0c256..f51f7ffa38 100644 --- a/source/source_lcao/module_deepks/LCAO_deepks.cpp +++ b/source/source_lcao/module_deepks/LCAO_deepks.cpp @@ -20,7 +20,7 @@ template LCAO_Deepks::LCAO_Deepks() { - inl_index = new ModuleBase::IntArray[1]; + deepks_param.inl_index = new ModuleBase::IntArray[1]; gedm = nullptr; this->phialpha.resize(1); } @@ -29,7 +29,7 @@ LCAO_Deepks::LCAO_Deepks() template LCAO_Deepks::~LCAO_Deepks() { - delete[] inl_index; + delete[] deepks_param.inl_index; //=======1. to use deepks, pdm is required========== pdm.clear(); @@ -38,7 +38,7 @@ LCAO_Deepks::~LCAO_Deepks() if (gedm) { // delete gedm** - for (int inl = 0; inl < this->inlmax; inl++) + for (int inl = 0; inl < this->deepks_param.inlmax; inl++) { delete[] gedm[inl]; } @@ -75,25 +75,25 @@ void LCAO_Deepks::init(const LCAO_Orbitals& orb, tot_inl = nat; } - this->lmaxd = lm; - this->nmaxd = nm; + this->deepks_param.lmaxd = lm; + this->deepks_param.nmaxd = nm; - ofs << " lmax of descriptor = " << this->lmaxd << std::endl; - ofs << " nmax of descriptor = " << nmaxd << std::endl; + ofs << " lmax of descriptor = " << deepks_param.lmaxd << std::endl; + ofs << " nmax of descriptor = " << deepks_param.nmaxd << std::endl; int pdm_size = 0; - this->inlmax = tot_inl; - this->pdm.resize(this->inlmax); + this->deepks_param.inlmax = tot_inl; + this->pdm.resize(this->deepks_param.inlmax); // cal n(descriptor) per atom , related to Lmax, nchi(L) and m. (not total_nchi!) if (!PARAM.inp.deepks_equiv) { - this->des_per_atom = 0; // mohan add 2021-04-21 - for (int l = 0; l <= this->lmaxd; l++) + this->deepks_param.des_per_atom = 0; // mohan add 2021-04-21 + for (int l = 0; l <= deepks_param.lmaxd; l++) { - this->des_per_atom += orb.Alpha[0].getNchi(l) * (2 * l + 1); + this->deepks_param.des_per_atom += orb.Alpha[0].getNchi(l) * (2 * l + 1); } - this->n_descriptor = nat * this->des_per_atom; + this->deepks_param.n_descriptor = nat * this->deepks_param.des_per_atom; this->init_index(ntype, nat, na, tot_inl, orb, ofs); } @@ -103,21 +103,21 @@ void LCAO_Deepks::init(const LCAO_Orbitals& orb, ofs << " total basis (all atoms) for descriptor = " << std::endl; // init pdm - for (int inl = 0; inl < this->inlmax; inl++) + for (int inl = 0; inl < this->deepks_param.inlmax; inl++) { - int nm = 2 * inl2l[inl] + 1; + int nm = 2 * deepks_param.inl2l[inl] + 1; pdm_size += nm * nm; this->pdm[inl] = torch::zeros({nm, nm}, torch::kFloat64); } } else { - for (int il = 0; il < this->lmaxd + 1; il++) + for (int il = 0; il < this->deepks_param.lmaxd + 1; il++) { pdm_size += (2 * il + 1) * orb.Alpha[0].getNchi(il); } pdm_size = pdm_size * pdm_size; - this->des_per_atom = pdm_size; + this->deepks_param.des_per_atom = pdm_size; ofs << " Equivariant version, size of pdm matrices : " << pdm_size << std::endl; for (int iat = 0; iat < nat; iat++) { @@ -139,35 +139,35 @@ void LCAO_Deepks::init_index(const int ntype, const LCAO_Orbitals& orb, std::ofstream& ofs) { - delete[] this->inl_index; - this->inl_index = new ModuleBase::IntArray[ntype]; - this->inl2l.resize(this->inlmax, 0); + delete[] this->deepks_param.inl_index; + this->deepks_param.inl_index = new ModuleBase::IntArray[ntype]; + this->deepks_param.inl2l.resize(this->deepks_param.inlmax, 0); int inl = 0; int alpha = 0; for (int it = 0; it < ntype; it++) { - this->inl_index[it].create(na[it], this->lmaxd + 1, this->nmaxd); + this->deepks_param.inl_index[it].create(na[it], this->deepks_param.lmaxd + 1, this->deepks_param.nmaxd); ofs << " Type " << it + 1 << " number_of_atoms " << na[it] << std::endl; for (int ia = 0; ia < na[it]; ia++) { // alpha - for (int l = 0; l < this->lmaxd + 1; l++) + for (int l = 0; l < this->deepks_param.lmaxd + 1; l++) { for (int n = 0; n < orb.Alpha[0].getNchi(l); n++) { - this->inl_index[it](ia, l, n) = inl; - this->inl2l[inl] = l; + this->deepks_param.inl_index[it](ia, l, n) = inl; + this->deepks_param.inl2l[inl] = l; inl++; } } } // end ia } // end it assert(Total_nchi == inl); - ofs << " descriptors_per_atom " << this->des_per_atom << std::endl; - ofs << " total_descriptors " << this->n_descriptor << std::endl; + ofs << " descriptors_per_atom " << this->deepks_param.des_per_atom << std::endl; + ofs << " total_descriptors " << this->deepks_param.n_descriptor << std::endl; return; } @@ -189,15 +189,15 @@ void LCAO_Deepks::allocate_V_delta(const int nat, const int nks) int pdm_size = 0; if (!PARAM.inp.deepks_equiv) { - pdm_size = (this->lmaxd * 2 + 1) * (this->lmaxd * 2 + 1); + pdm_size = (this->deepks_param.lmaxd * 2 + 1) * (this->deepks_param.lmaxd * 2 + 1); } else { - pdm_size = this->des_per_atom; + pdm_size = this->deepks_param.des_per_atom; } - this->gedm = new double*[this->inlmax]; - for (int inl = 0; inl < this->inlmax; inl++) + this->gedm = new double*[this->deepks_param.inlmax]; + for (int inl = 0; inl < this->deepks_param.inlmax; inl++) { this->gedm[inl] = new double[pdm_size]; ModuleBase::GlobalFunc::ZEROS(this->gedm[inl], pdm_size); @@ -214,44 +214,41 @@ void LCAO_Deepks::init_DMR(const UnitCell& ucell, const Grid_Driver& GridD) { this->dm_r = new hamilt::HContainer(&pv); - DeePKS_domain::iterate_ad2( - ucell, - GridD, - orb, - false, // no trace_alpha - [&](const int iat, - const ModuleBase::Vector3& tau0, - const int ibt1, - const ModuleBase::Vector3& tau1, - const int start1, - const int nw1_tot, - ModuleBase::Vector3 dR1, - const int ibt2, - const ModuleBase::Vector3& tau2, - const int start2, - const int nw2_tot, - ModuleBase::Vector3 dR2) - { - auto row_indexes = pv.get_indexes_row(ibt1); - auto col_indexes = pv.get_indexes_col(ibt2); - if (row_indexes.size() * col_indexes.size() == 0) - { - return; // to next loop - } - - int dRx = 0; - int dRy = 0; - int dRz = 0; - if (std::is_same>::value) - { - dRx = (dR1 - dR2).x; - dRy = (dR1 - dR2).y; - dRz = (dR1 - dR2).z; - } - hamilt::AtomPair dm_pair(ibt1, ibt2, dRx, dRy, dRz, &pv); - this->dm_r->insert_pair(dm_pair); - } - ); + DeePKS_domain::iterate_ad2(ucell, + GridD, + orb, + false, // no trace_alpha + [&](const int iat, + const ModuleBase::Vector3& tau0, + const int ibt1, + const ModuleBase::Vector3& tau1, + const int start1, + const int nw1_tot, + ModuleBase::Vector3 dR1, + const int ibt2, + const ModuleBase::Vector3& tau2, + const int start2, + const int nw2_tot, + ModuleBase::Vector3 dR2) { + auto row_indexes = pv.get_indexes_row(ibt1); + auto col_indexes = pv.get_indexes_col(ibt2); + if (row_indexes.size() * col_indexes.size() == 0) + { + return; // to next loop + } + + int dRx = 0; + int dRy = 0; + int dRz = 0; + if (std::is_same>::value) + { + dRx = (dR1 - dR2).x; + dRy = (dR1 - dR2).y; + dRz = (dR1 - dR2).z; + } + hamilt::AtomPair dm_pair(ibt1, ibt2, dRx, dRy, dRz, &pv); + this->dm_r->insert_pair(dm_pair); + }); this->dm_r->allocate(nullptr, true); } diff --git a/source/source_lcao/module_deepks/LCAO_deepks.h b/source/source_lcao/module_deepks/LCAO_deepks.h index e3a43cd041..bb98f75a4b 100644 --- a/source/source_lcao/module_deepks/LCAO_deepks.h +++ b/source/source_lcao/module_deepks/LCAO_deepks.h @@ -10,6 +10,7 @@ #include "deepks_fpre.h" #include "deepks_orbital.h" #include "deepks_orbpre.h" +#include "deepks_param.h" #include "deepks_pdm.h" #include "deepks_phialpha.h" #include "deepks_spre.h" @@ -67,14 +68,8 @@ class LCAO_Deepks // 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)} - std::vector inl2l; // inl2l[inl] = inl2l[nl] = l (not related to iat) of descriptor with inl_index - ModuleBase::IntArray* inl_index; // caoyu add 2021-05-07 + public: // change to public to reconstuct the code, 2024-07-22 by mohan + DeePKS_Param deepks_param; // parameters for DeePKS bool init_pdm = false; // for DeePKS NSCF calculation, set init_pdm to skip the calculation of pdm in SCF iteration diff --git a/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp b/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp index 8e30268728..1d68dbc603 100644 --- a/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp +++ b/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp @@ -2,13 +2,13 @@ #include "LCAO_deepks_interface.h" #include "LCAO_deepks_io.h" // mohan add 2024-07-22 +#include "source_base/global_variable.h" +#include "source_base/tool_title.h" #include "source_estate/cal_dm.h" +#include "source_io/module_parameter/parameter.h" #include "source_lcao/module_hcontainer/hcontainer.h" #include "source_lcao/module_hcontainer/hcontainer_funcs.h" #include "source_lcao/module_hcontainer/output_hcontainer.h" -#include "source_io/module_parameter/parameter.h" -#include "source_base/global_variable.h" -#include "source_base/tool_title.h" #include @@ -18,15 +18,10 @@ LCAO_Deepks_Interface::LCAO_Deepks_Interface(std::shared_ptr file_type_map = { - {"etot", "energy"}, - {"ftot", "force"}, - {"stot", "stress"}, - {"otot", "orbital"}, - {"htot", "hamiltonian"} - }; + static const std::unordered_map file_type_map + = {{"etot", "energy"}, {"ftot", "force"}, {"stot", "stress"}, {"otot", "orbital"}, {"htot", "hamiltonian"}}; auto it = file_type_map.find(file_type); return it != file_type_map.end() ? it->second : file_type; @@ -37,11 +32,11 @@ std::string true_file_type(const std::string& file_type) template std::string LCAO_Deepks_Interface::get_filename(const std::string& file_type, const int& label_type, - const int& iter) + const int& iter) { std::ostringstream file_name; file_name << (iter == -1 ? PARAM.globalv.global_out_dir : PARAM.globalv.global_deepks_label_elec_dir); - if (iter == -1) + if (iter == -1) { file_name << "deepks_"; } @@ -82,16 +77,10 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, using TH = std::conditional_t::value, ModuleBase::matrix, ModuleBase::ComplexMatrix>; // These variables are frequently used in the following code - const int nlmax = orb.Alpha[0].getTotal_nchi(); - const int inlmax = nlmax * nat; - const int lmaxd = orb.get_lmax_d(); - const int nmaxd = ld->nmaxd; - - const int des_per_atom = ld->des_per_atom; - const std::vector inl2l = ld->inl2l; - const ModuleBase::IntArray* inl_index = ld->inl_index; const std::vector*> phialpha = ld->phialpha; + const DeePKS_Param& deepks_param = ld->deepks_param; + std::vector pdm = ld->pdm; bool init_pdm = ld->init_pdm; double E_delta = ld->E_delta; @@ -104,13 +93,16 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, const int nk = nks / nspin; const bool is_after_scf = (iter == -1); // called in after_scf, not in electronic steps - const bool output_base = PARAM.inp.deepks_out_labels == 1 && is_after_scf; // not output when deepks_out_labels=2 and in electronic step (output true base elsewhere) - const bool output_precalc = (PARAM.inp.deepks_out_labels == 1) && (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_freq_elec); + const bool output_base + = PARAM.inp.deepks_out_labels == 1 + && is_after_scf; // not output when deepks_out_labels=2 and in electronic step (output true base elsewhere) + const bool output_precalc + = (PARAM.inp.deepks_out_labels == 1) && (PARAM.inp.deepks_scf || PARAM.inp.deepks_out_freq_elec); -//================================================================================ -// 1. Update real-space density matrix (DMR) for deepks, projected density matrix (PDM) -// and descriptor. Output descriptor if needed. -//================================================================================ + //================================================================================ + // 1. Update real-space density matrix (DMR) for deepks, projected density matrix (PDM) + // and descriptor. Output descriptor if needed. + //================================================================================ // Update DMR in any case of deepks_out_labels/deepks_scf DeePKS_domain::update_dmr(kvec_d, dm->get_DMK_vector(), ucell, orb, *ParaV, GridD, dmr); @@ -121,32 +113,21 @@ 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 - DeePKS_domain::cal_pdm< - TK>(init_pdm, inlmax, lmaxd, inl2l, inl_index, kvec_d, dmr, phialpha, ucell, orb, GridD, *ParaV, pdm); + DeePKS_domain::cal_pdm(init_pdm, deepks_param, kvec_d, dmr, phialpha, ucell, orb, GridD, *ParaV, pdm); - DeePKS_domain::check_pdm(inlmax, inl2l, pdm); // print out the projected dm for NSCF calculaiton + DeePKS_domain::check_pdm(deepks_param, pdm); // print out the projected dm for NSCF calculaiton std::vector descriptor; - DeePKS_domain::cal_descriptor(nat, inlmax, inl2l, pdm, descriptor, - des_per_atom); // final descriptor - DeePKS_domain::check_descriptor(inlmax, - des_per_atom, - inl2l, - ucell, - PARAM.globalv.global_out_dir, - descriptor, - rank); - + DeePKS_domain::cal_descriptor(nat, deepks_param, pdm, descriptor); // final descriptor + DeePKS_domain::check_descriptor(deepks_param, ucell, PARAM.globalv.global_out_dir, descriptor, rank); + const std::string file_d = get_filename("dm_eig", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_npy_d(nat, - des_per_atom, - inlmax, - inl2l, - PARAM.inp.deepks_equiv, - descriptor, - file_d, - rank); // libnpy needed - + PARAM.inp.deepks_equiv, + deepks_param, + descriptor, + file_d, + rank); // libnpy needed if (PARAM.inp.deepks_scf) { @@ -154,28 +135,11 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, // new gedm is also useful in cal_f_delta, so it should be ld->gedm if (PARAM.inp.deepks_equiv) { - DeePKS_domain::cal_edelta_gedm_equiv(nat, - lmaxd, - nmaxd, - inlmax, - des_per_atom, - inl2l, - descriptor, - ld->gedm, - E_delta, - rank); + DeePKS_domain::cal_edelta_gedm_equiv(nat, deepks_param, descriptor, ld->gedm, E_delta, rank); } else { - DeePKS_domain::cal_edelta_gedm(nat, - inlmax, - des_per_atom, - inl2l, - descriptor, - pdm, - ld->model_deepks, - ld->gedm, - E_delta); + DeePKS_domain::cal_edelta_gedm(nat, deepks_param, descriptor, pdm, ld->model_deepks, ld->gedm, E_delta); } } } @@ -185,21 +149,21 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, { // Used for deepks_scf == 1 or deepks_out_freq_elec!=0, for *precalc items, not for deepks_out_labels=2 std::vector gevdm; - if ( output_precalc ) + if (output_precalc) { - DeePKS_domain::cal_gevdm(nat, inlmax, inl2l, pdm, gevdm); + DeePKS_domain::cal_gevdm(nat, deepks_param, pdm, gevdm); } -//================================================================================ -// 2. Energy -//================================================================================ + //================================================================================ + // 2. Energy + //================================================================================ // etot const std::string file_etot = get_filename("etot", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_npy_e(etot, file_etot, rank); // ebase - if ( output_base ) + if (output_base) { const std::string file_ebase = get_filename("ebase", PARAM.inp.deepks_out_labels, iter); if (PARAM.inp.deepks_scf) @@ -214,25 +178,25 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, } } -//================================================================================ -// 3. Force and Stress -//================================================================================ + //================================================================================ + // 3. Force and Stress + //================================================================================ - if ( is_after_scf ) + if (is_after_scf) { // Force Part if (PARAM.inp.cal_force) { // these items are not related to model, so can output without deepks_scf - if (output_precalc // don't need these when deepks_out_labels == 2 + if (output_precalc // don't need these when deepks_out_labels == 2 && !PARAM.inp.deepks_equiv) // training with force label not supported by equivariant version now { torch::Tensor gdmx; DeePKS_domain::cal_gdmx< - TK>(lmaxd, inlmax, nks, kvec_d, phialpha, inl_index, dmr, ucell, orb, *ParaV, GridD, gdmx); + TK>(nks, deepks_param, kvec_d, phialpha, dmr, ucell, orb, *ParaV, GridD, gdmx); torch::Tensor gvx; - DeePKS_domain::cal_gvx(ucell.nat, inlmax, des_per_atom, inl2l, gevdm, gdmx, gvx, rank); + DeePKS_domain::cal_gvx(ucell.nat, deepks_param, gevdm, gdmx, gvx, rank); const std::string file_gradvx = get_filename("gradvx", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_tensor2npy(file_gradvx, gvx, rank); @@ -248,15 +212,15 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, if (PARAM.inp.cal_stress) { // these items are not related to model, so can output without deepks_scf - if (output_precalc // don't need these when deepks_out_labels == 2 + if (output_precalc // don't need these when deepks_out_labels == 2 && !PARAM.inp.deepks_equiv) // training with stress label not supported by equivariant version now { torch::Tensor gdmepsl; DeePKS_domain::cal_gdmepsl< - TK>(lmaxd, inlmax, nks, kvec_d, phialpha, inl_index, dmr, ucell, orb, *ParaV, GridD, gdmepsl); + TK>(nks, deepks_param, kvec_d, phialpha, dmr, ucell, orb, *ParaV, GridD, gdmepsl); torch::Tensor gvepsl; - DeePKS_domain::cal_gvepsl(ucell.nat, inlmax, des_per_atom, inl2l, gevdm, gdmepsl, gvepsl, rank); + DeePKS_domain::cal_gvepsl(ucell.nat, deepks_param, gevdm, gdmepsl, gvepsl, rank); const std::string file_gvepsl = get_filename("gvepsl", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_tensor2npy(file_gvepsl, gvepsl, rank); @@ -269,9 +233,9 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, } } -//================================================================================ -// 4. Bandgap -//================================================================================ + //================================================================================ + // 4. Bandgap + //================================================================================ if (PARAM.inp.deepks_bandgap > 0) { @@ -308,7 +272,9 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, ModuleBase::matrix o_tot(nks, range); if (nocc + PARAM.inp.deepks_band_range[0] < 0 || nocc + PARAM.inp.deepks_band_range[1] >= ekb.nc) { - ModuleBase::WARNING_QUIT("out_deepks_labels", "DeePKS band index out of range! Please check if `deepks_band_range` is reasonable!"); + ModuleBase::WARNING_QUIT( + "out_deepks_labels", + "DeePKS band index out of range! Please check if `deepks_band_range` is reasonable!"); } for (int iks = 0; iks < nks; ++iks) { @@ -316,7 +282,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, if (PARAM.inp.deepks_bandgap == 1 || PARAM.inp.deepks_bandgap == 3) { o_tot(iks, ib) = ekb(iks, nocc + PARAM.inp.deepks_band_range[1]) - - ekb(iks, nocc + PARAM.inp.deepks_band_range[0]); + - ekb(iks, nocc + PARAM.inp.deepks_band_range[0]); } else if (PARAM.inp.deepks_bandgap == 2) { @@ -337,7 +303,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, // don't need these when deepks_out_labels == 2 // not consider out_base now, because bandgap is not supported in deepks_out_freq_elec now - if ( output_precalc ) + if (output_precalc) { std::vector wg_hl_range(range); for (int ir = 0; ir < range; ++ir) @@ -379,20 +345,17 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, torch::Tensor orbital_precalc_temp; ModuleBase::matrix o_delta_temp(nks, 1); DeePKS_domain::cal_orbital_precalc(dm_bandgap, - lmaxd, - inlmax, - nat, - nks, - inl2l, - kvec_d, - phialpha, - gevdm, - inl_index, - ucell, - orb, - *ParaV, - GridD, - orbital_precalc_temp); + nat, + nks, + deepks_param, + kvec_d, + phialpha, + gevdm, + ucell, + orb, + *ParaV, + GridD, + orbital_precalc_temp); if (ir == 0) { orbital_precalc = orbital_precalc_temp; @@ -408,30 +371,30 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, for (int iks = 0; iks < nks; ++iks) { o_delta(iks, ir) = o_delta_temp(iks, 0); - } + } } } // save obase and orbital_precalc const std::string file_orbpre = get_filename("orbpre", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_tensor2npy(file_orbpre, orbital_precalc, rank); + LCAO_deepks_io::save_tensor2npy(file_orbpre, orbital_precalc, rank); if (PARAM.inp.deepks_scf) { const std::string file_obase = get_filename("obase", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_matrix2npy(file_obase, o_tot - o_delta, rank); // Unit: Hartree + LCAO_deepks_io::save_matrix2npy(file_obase, o_tot - o_delta, rank); // Unit: Hartree } else { const std::string file_obase = get_filename("obase", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_matrix2npy(file_obase, o_tot, rank); // no scf, o_tot=o_base - } - } // end out_precalc - } // end deepks_bandgap > 0 - -//================================================================================ -// 5. HR -//================================================================================ - if ( is_after_scf ) + } + } // end out_precalc + } // end deepks_bandgap > 0 + + //================================================================================ + // 5. HR + //================================================================================ + if (is_after_scf) { // not add deepks_out_labels = 2 and deepks_out_freq_elec for HR yet // H(R) matrix part, for HR, base will not be calculated since they are HContainer objects @@ -442,7 +405,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, const int precision = 8; const std::string file_hrtot = PARAM.globalv.global_out_dir - + (PARAM.inp.deepks_out_labels == 1 ? "deepks_hrtot.csr" : "deepks_hamiltonian_r.csr"); + + (PARAM.inp.deepks_out_labels == 1 ? "deepks_hrtot.csr" : "deepks_hamiltonian_r.csr"); hamilt::HContainer* hR_tot = (p_ham->getHR()); const int nbasis = hR_tot->get_nbasis(); @@ -494,56 +457,46 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, int R_size = DeePKS_domain::get_R_size(*h_deltaR); torch::Tensor vdr_precalc; DeePKS_domain::cal_vdr_precalc(nlocal, - lmaxd, - inlmax, - nat, - nks, - R_size, - inl2l, - kvec_d, - phialpha, - gevdm, - inl_index, - ucell, - orb, - *ParaV, - GridD, - vdr_precalc); + nat, + nks, + R_size, + deepks_param, + kvec_d, + phialpha, + gevdm, + ucell, + orb, + *ParaV, + GridD, + vdr_precalc); const std::string file_vdrpre = PARAM.globalv.global_out_dir + "deepks_vdrpre.npy"; LCAO_deepks_io::save_tensor2npy(file_vdrpre, vdr_precalc, rank); } else if (PARAM.inp.deepks_v_delta == -2) { - int R_size = DeePKS_domain::get_R_size(*h_deltaR); - torch::Tensor phialpha_r_out; - DeePKS_domain::prepare_phialpha_r(nlocal, - lmaxd, - inlmax, - nat, - R_size, - phialpha, - ucell, - orb, - *ParaV, - GridD, - phialpha_r_out); - const std::string file_phialpha_r = PARAM.globalv.global_out_dir + "deepks_phialpha_r.npy"; - LCAO_deepks_io::save_tensor2npy(file_phialpha_r, phialpha_r_out, rank); - torch::Tensor gevdm_out; - DeePKS_domain::prepare_gevdm(nat, lmaxd, inlmax, orb, gevdm, gevdm_out); + DeePKS_domain::prepare_gevdm(nat, deepks_param, orb, gevdm, gevdm_out); const std::string file_gevdm = PARAM.globalv.global_out_dir + "deepks_gevdm.npy"; LCAO_deepks_io::save_tensor2npy(file_gevdm, gevdm_out, rank); + + int R_size = DeePKS_domain::get_R_size(*h_deltaR); + torch::Tensor overlap_out; + torch::Tensor iRmat; + DeePKS_domain::prepare_phialpha_iRmat(nlocal, R_size, deepks_param, phialpha, ucell, orb, GridD, overlap_out, iRmat); + const std::string file_overlap = PARAM.globalv.global_out_dir + "deepks_phialpha_r.npy"; + LCAO_deepks_io::save_tensor2npy(file_overlap, overlap_out, rank); + const std::string file_iRmat = PARAM.globalv.global_out_dir + "deepks_iRmat.npy"; + LCAO_deepks_io::save_tensor2npy(file_iRmat, iRmat, rank); } } } - } + } } -//================================================================================ -// 6. Hk -//================================================================================ + //================================================================================ + // 6. Hk + //================================================================================ if (PARAM.inp.deepks_v_delta > 0) { @@ -553,7 +506,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, const std::string file_htot = get_filename("htot", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_npy_h(h_tot, file_htot, nlocal, nks, rank); - if ( output_base ) + if (output_base) { if (PARAM.inp.deepks_scf) { @@ -575,35 +528,32 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, LCAO_deepks_io::save_npy_h(h_base, file_hbase, nlocal, nks, rank); const std::string file_vdelta = get_filename("vdelta", PARAM.inp.deepks_out_labels, iter); - LCAO_deepks_io::save_npy_h(v_delta, file_vdelta, nlocal, nks, rank); + LCAO_deepks_io::save_npy_h(v_delta, file_vdelta, nlocal, nks, rank); } else // deepks_scf == 0 { const std::string file_hbase = get_filename("hbase", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_npy_h(h_tot, file_hbase, nlocal, nks, rank); - } + } } - if ( output_precalc ) + if (output_precalc) { if (PARAM.inp.deepks_v_delta == 1) // v_delta_precalc storage method 1 { torch::Tensor v_delta_precalc; DeePKS_domain::cal_v_delta_precalc(nlocal, - lmaxd, - inlmax, - nat, - nks, - inl2l, - kvec_d, - phialpha, - gevdm, - inl_index, - ucell, - orb, - *ParaV, - GridD, - v_delta_precalc); + nat, + nks, + deepks_param, + kvec_d, + phialpha, + gevdm, + ucell, + orb, + *ParaV, + GridD, + v_delta_precalc); const std::string file_vdpre = get_filename("vdpre", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_tensor2npy(file_vdpre, v_delta_precalc, rank); @@ -611,37 +561,27 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, else if (PARAM.inp.deepks_v_delta == 2) // v_delta_precalc storage method 2 { torch::Tensor phialpha_out; - DeePKS_domain::prepare_phialpha(nlocal, - lmaxd, - inlmax, - nat, - nks, - kvec_d, - phialpha, - ucell, - orb, - *ParaV, - GridD, - phialpha_out); + DeePKS_domain::prepare_phialpha< + TK>(nlocal, nat, nks, deepks_param, kvec_d, phialpha, ucell, orb, *ParaV, GridD, phialpha_out); const std::string file_phialpha = get_filename("phialpha", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_tensor2npy(file_phialpha, phialpha_out, rank); torch::Tensor gevdm_out; - DeePKS_domain::prepare_gevdm(nat, lmaxd, inlmax, orb, gevdm, gevdm_out); + DeePKS_domain::prepare_gevdm(nat, deepks_param, orb, gevdm, gevdm_out); const std::string file_gevdm = get_filename("gevdm", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_tensor2npy(file_gevdm, gevdm_out, rank); } } - } // end v_delta label + } // end v_delta label } // end deepks_out_labels -//================================================================================ -// 7. atom.npy, box.npy, overlap.npy -//================================================================================ + //================================================================================ + // 7. atom.npy, box.npy, overlap.npy + //================================================================================ - if( ((PARAM.inp.deepks_out_labels == 2) && is_after_scf ) - || ( PARAM.inp.deepks_out_freq_elec ) )// need overlap when deepks_out_freq_elec + if (((PARAM.inp.deepks_out_labels == 2) && is_after_scf) + || (PARAM.inp.deepks_out_freq_elec)) // need overlap when deepks_out_freq_elec { if (PARAM.inp.deepks_v_delta > 0) { @@ -650,15 +590,15 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, DeePKS_domain::get_h_tot(*ParaV, p_ham, s_tot, nlocal, nks, 'S'); const std::string file_stot = get_filename("overlap", PARAM.inp.deepks_out_labels, iter); LCAO_deepks_io::save_npy_h(s_tot, - file_stot, - nlocal, - nks, - rank, - 1.0); // don't need unit_scale for overlap + file_stot, + nlocal, + nks, + rank, + 1.0); // don't need unit_scale for overlap } } - if ( is_after_scf ) // don't need to output in multiple electronic steps + if (is_after_scf) // don't need to output in multiple electronic steps { if (PARAM.inp.deepks_out_labels == 2) { @@ -674,9 +614,9 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, LCAO_deepks_io::save_tensor2npy(file_box, box_out, rank); } -//================================================================================ -// 8. print, unittest -//================================================================================ + //================================================================================ + // 8. print, unittest + //================================================================================ /// print out deepks information to the screen if (PARAM.inp.deepks_scf) { @@ -686,7 +626,8 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, ofs_running << " DeePKS Energy Correction" << std::endl; ofs_running << " -----------------------------------------------" << std::endl; ofs_running << " E_delta_band = " << std::setprecision(8) << e_delta_band << " Ry" - << " = " << std::setprecision(8) << e_delta_band * ModuleBase::Ry_to_eV << " eV" << std::endl; + << " = " << std::setprecision(8) << e_delta_band * ModuleBase::Ry_to_eV << " eV" + << std::endl; ofs_running << " E_delta_NN = " << std::setprecision(8) << E_delta << " Ry" << " = " << std::setprecision(8) << E_delta * ModuleBase::Ry_to_eV << " eV" << std::endl; ofs_running << " -----------------------------------------------" << std::endl; @@ -695,7 +636,7 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, { LCAO_deepks_io::print_dm(nks, PARAM.globalv.nlocal, ParaV->nrow, dm->get_DMK_vector()); - DeePKS_domain::check_gedm(inlmax, inl2l, ld->gedm); + DeePKS_domain::check_gedm(deepks_param, ld->gedm); std::ofstream ofs("E_delta_bands.dat"); ofs << std::setprecision(10) << e_delta_band; @@ -703,8 +644,8 @@ void LCAO_Deepks_Interface::out_deepks_labels(const double& etot, std::ofstream ofs1("E_delta.dat"); ofs1 << std::setprecision(10) << E_delta; } - } - } + } + } ModuleBase::timer::tick("LCAO_Deepks_Interface", "out_deepks_labels"); } diff --git a/source/source_lcao/module_deepks/LCAO_deepks_interface.h b/source/source_lcao/module_deepks/LCAO_deepks_interface.h index a637e8839d..bc60edab04 100644 --- a/source/source_lcao/module_deepks/LCAO_deepks_interface.h +++ b/source/source_lcao/module_deepks/LCAO_deepks_interface.h @@ -3,9 +3,9 @@ #ifdef __MLALGO #include "LCAO_deepks.h" -#include "source_lcao/hamilt_lcao.h" #include "source_base/complexmatrix.h" #include "source_base/matrix.h" +#include "source_lcao/hamilt_lcao.h" #include @@ -47,18 +47,16 @@ class LCAO_Deepks_Interface const elecstate::DensityMatrix* dm, hamilt::HamiltLCAO* p_ham, const int& iter, - const bool& conv_esolver, + const bool& conv_esolver, const int rank, std::ostream& ofs_running); - + /// @brief Get the filename for deepks output files /// @param file_type Type of the file (e.g., "etot", "ftot", etc.) /// @param label_type Type of the label (from PARAM.inp.deepks_out_labels) /// @param iter Iteration number (e.g., -1 for after_scf, or specific iteration number) /// @return The full path to the output file - std::string get_filename(const std::string& file_type, - const int& label_type, - const int& iter); + std::string get_filename(const std::string& file_type, const int& label_type, const int& iter); private: std::shared_ptr> ld; diff --git a/source/source_lcao/module_deepks/LCAO_deepks_io.cpp b/source/source_lcao/module_deepks/LCAO_deepks_io.cpp index 9ade40e34a..12c38dfcf4 100644 --- a/source/source_lcao/module_deepks/LCAO_deepks_io.cpp +++ b/source/source_lcao/module_deepks/LCAO_deepks_io.cpp @@ -3,8 +3,8 @@ #ifdef __MLALGO #include "LCAO_deepks_io.h" -#include "source_base/tool_quit.h" #include "npy.hpp" +#include "source_base/tool_quit.h" #include @@ -75,10 +75,8 @@ void LCAO_deepks_io::load_npy_gedm(const int nat, // saves descriptor into dm_eig.npy void LCAO_deepks_io::save_npy_d(const int nat, - const int des_per_atom, - const int inlmax, - const std::vector& inl2l, const bool deepks_equiv, + const DeePKS_Param& deepks_param, const std::vector& descriptor, const std::string& dm_eig_file, const int rank) @@ -95,16 +93,17 @@ void LCAO_deepks_io::save_npy_d(const int nat, if (!deepks_equiv) { std::vector npy_des; - for (int inl = 0; inl < inlmax; ++inl) + for (int inl = 0; inl < deepks_param.inlmax; ++inl) { auto accessor = descriptor[inl].accessor(); - int nm = 2 * inl2l[inl] + 1; + int nm = 2 * deepks_param.inl2l[inl] + 1; for (int im = 0; im < nm; im++) { npy_des.push_back(accessor[im]); } } - const long unsigned dshape[] = {static_cast(nat), static_cast(des_per_atom)}; + const long unsigned dshape[] + = {static_cast(nat), static_cast(deepks_param.des_per_atom)}; if (rank == 0) { npy::SaveArrayAsNumpy(dm_eig_file, false, 2, dshape, npy_des); @@ -117,12 +116,13 @@ void LCAO_deepks_io::save_npy_d(const int nat, for (int iat = 0; iat < nat; iat++) { auto accessor = descriptor[iat].accessor(); - for (int i = 0; i < des_per_atom; i++) + for (int i = 0; i < deepks_param.des_per_atom; i++) { npy_des.push_back(accessor[i]); } } - const long unsigned dshape[] = {static_cast(nat), static_cast(des_per_atom)}; + const long unsigned dshape[] + = {static_cast(nat), static_cast(deepks_param.des_per_atom)}; if (rank == 0) { npy::SaveArrayAsNumpy(dm_eig_file, false, 2, dshape, npy_des); diff --git a/source/source_lcao/module_deepks/LCAO_deepks_io.h b/source/source_lcao/module_deepks/LCAO_deepks_io.h index ba53ce6388..46e9272508 100644 --- a/source/source_lcao/module_deepks/LCAO_deepks_io.h +++ b/source/source_lcao/module_deepks/LCAO_deepks_io.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/complexmatrix.h" #include "source_base/matrix.h" #include "source_base/tool_title.h" @@ -44,10 +45,8 @@ void load_npy_gedm(const int nat, const int des_per_atom, double** gedm, double& /// save descriptor void save_npy_d(const int nat, - const int des_per_atom, - const int inlmax, - const std::vector& inl2l, const bool deepks_equiv, + const DeePKS_Param& deepks_param, const std::vector& descriptor, const std::string& dm_eig_file, const int rank); diff --git a/source/source_lcao/module_deepks/deepks_basic.cpp b/source/source_lcao/module_deepks/deepks_basic.cpp index 4a30d3af39..aa1185a23b 100644 --- a/source/source_lcao/module_deepks/deepks_basic.cpp +++ b/source/source_lcao/module_deepks/deepks_basic.cpp @@ -8,27 +8,27 @@ #include "source_base/atom_in.h" #include "source_base/timer.h" #include "source_io/module_parameter/parameter.h" + #include // use system command // d(Descriptor) / d(projected density matrix) // Dimension is different for each inl, so there's a vector of tensors void DeePKS_domain::cal_gevdm(const int nat, - const int inlmax, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& pdm, std::vector& gevdm) { ModuleBase::TITLE("DeePKS_domain", "cal_gevdm"); ModuleBase::timer::tick("DeePKS_domain", "cal_gevdm"); // cal gevdm(d(EigenValue(D))/dD) - int nlmax = inlmax / nat; + int nlmax = deepks_param.inlmax / nat; for (int nl = 0; nl < nlmax; ++nl) { std::vector avmmv; for (int iat = 0; iat < nat; ++iat) { int inl = iat * nlmax + nl; - int nm = 2 * inl2l[inl] + 1; + int nm = 2 * deepks_param.inl2l[inl] + 1; // repeat each block for nm times in an additional dimension torch::Tensor tmp_x = pdm[inl].reshape({nm, nm}).unsqueeze(0).repeat({nm, 1, 1}); // torch::Tensor tmp_y = std::get<0>(torch::symeig(tmp_x, true)); @@ -85,7 +85,7 @@ void DeePKS_domain::load_model(const std::string& model_file, torch::jit::script return; } -inline void generate_py_files(const int lmaxd, const int nmaxd, const std::string& out_dir) +inline void generate_py_files(const DeePKS_Param& deepks_param, const std::string& out_dir) { std::ofstream ofs("cal_edelta_gedm.py"); ofs << "import torch" << std::endl; @@ -117,14 +117,14 @@ inline void generate_py_files(const int lmaxd, const int nmaxd, const std::strin ofs.open("basis.yaml"); ofs << "proj_basis:" << std::endl; - for (int l = 0; l < lmaxd + 1; l++) + for (int l = 0; l < deepks_param.lmaxd + 1; l++) { ofs << " - - " << l << std::endl; ofs << " - ["; - for (int i = 0; i < nmaxd + 1; i++) + for (int i = 0; i < deepks_param.nmaxd + 1; i++) { ofs << "0"; - if (i != nmaxd) + if (i != deepks_param.nmaxd) { ofs << ", "; } @@ -134,11 +134,7 @@ inline void generate_py_files(const int lmaxd, const int nmaxd, const std::strin } void DeePKS_domain::cal_edelta_gedm_equiv(const int nat, - const int lmaxd, - const int nmaxd, - const int inlmax, - const int des_per_atom, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& descriptor, double** gedm, double& E_delta, @@ -147,19 +143,13 @@ void DeePKS_domain::cal_edelta_gedm_equiv(const int nat, ModuleBase::TITLE("DeePKS_domain", "cal_edelta_gedm_equiv"); ModuleBase::timer::tick("DeePKS_domain", "cal_edelta_gedm_equiv"); - const std::string file_d = PARAM.globalv.global_out_dir + "deepks_dm_eig.npy";; - LCAO_deepks_io::save_npy_d(nat, - des_per_atom, - inlmax, - inl2l, - PARAM.inp.deepks_equiv, - descriptor, - file_d, + const std::string file_d = PARAM.globalv.global_out_dir + "deepks_dm_eig.npy"; + LCAO_deepks_io::save_npy_d(nat, PARAM.inp.deepks_equiv, deepks_param, descriptor, file_d, rank); // libnpy needed if (rank == 0) { - generate_py_files(lmaxd, nmaxd, PARAM.globalv.global_out_dir); + generate_py_files(deepks_param, PARAM.globalv.global_out_dir); std::string cmd = "python cal_edelta_gedm.py " + PARAM.inp.deepks_model; int stat = std::system(cmd.c_str()); assert(stat == 0); @@ -167,7 +157,7 @@ void DeePKS_domain::cal_edelta_gedm_equiv(const int nat, MPI_Barrier(MPI_COMM_WORLD); - LCAO_deepks_io::load_npy_gedm(nat, des_per_atom, gedm, E_delta, rank); + LCAO_deepks_io::load_npy_gedm(nat, deepks_param.des_per_atom, gedm, E_delta, rank); std::string cmd = "rm -f cal_edelta_gedm.py basis.yaml ec.npy gedm.npy"; std::system(cmd.c_str()); @@ -179,9 +169,7 @@ void DeePKS_domain::cal_edelta_gedm_equiv(const int nat, // obtain from the machine learning model dE_delta/dDescriptor // E_delta is also calculated here void DeePKS_domain::cal_edelta_gedm(const int nat, - const int inlmax, - const int des_per_atom, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& descriptor, const std::vector& pdm, torch::jit::script::Module& model_deepks, @@ -195,7 +183,7 @@ void DeePKS_domain::cal_edelta_gedm(const int nat, std::vector inputs; // input_dim:(natom, des_per_atom) - inputs.push_back(torch::cat(descriptor, 0).reshape({1, nat, des_per_atom})); + inputs.push_back(torch::cat(descriptor, 0).reshape({1, nat, deepks_param.des_per_atom})); std::vector ec; try { @@ -203,10 +191,22 @@ void DeePKS_domain::cal_edelta_gedm(const int nat, } catch (const c10::Error& e) { - ModuleBase::WARNING_QUIT("DeePKS_domain::cal_edelta_gedm", "Please check whether the input shape required by model file matches the descriptor!"); + ModuleBase::WARNING_QUIT("DeePKS_domain::cal_edelta_gedm", + "Please check whether the input shape required by model file matches the descriptor!"); throw; } - E_delta = ec[0].item() * 2; // Ry; *2 is for Hartree to Ry + E_delta = ec[0].item() * 2; // Ry; *2 is for Hartree to Ry + + // get d ec[0]/d inputs + // inputs: [1, nat, des_per_atom] + // ec: [1, 1] + std::vector tensor_inputs; + tensor_inputs.push_back(inputs[0].toTensor()); + ec[0].reshape({1, 1}).requires_grad_(true); + torch::Tensor derivative = torch::autograd::grad(ec, tensor_inputs, {}, true)[0]; + LCAO_deepks_io::save_tensor2npy("gev.npy", + derivative.reshape({nat, deepks_param.des_per_atom}), + 0); // dm_eig.npy is the input for gedm // cal gedm std::vector gedm_shell; @@ -219,9 +219,9 @@ void DeePKS_domain::cal_edelta_gedm(const int nat, /*allow_unused=*/true); // gedm_tensor(Hartree) to gedm(Ry) - for (int inl = 0; inl < inlmax; ++inl) + for (int inl = 0; inl < deepks_param.inlmax; ++inl) { - int nm = 2 * inl2l[inl] + 1; + int nm = 2 * deepks_param.inl2l[inl] + 1; auto accessor = gedm_tensor[inl].accessor(); for (int m1 = 0; m1 < nm; ++m1) { @@ -236,13 +236,13 @@ void DeePKS_domain::cal_edelta_gedm(const int nat, return; } -void DeePKS_domain::check_gedm(const int inlmax, const std::vector& inl2l, double** gedm) +void DeePKS_domain::check_gedm(const DeePKS_Param& deepks_param, double** gedm) { std::ofstream ofs("gedm.dat"); - for (int inl = 0; inl < inlmax; inl++) + for (int inl = 0; inl < deepks_param.inlmax; inl++) { - int nm = 2 * inl2l[inl] + 1; + int nm = 2 * deepks_param.inl2l[inl] + 1; for (int m1 = 0; m1 < nm; ++m1) { for (int m2 = 0; m2 < nm; ++m2) diff --git a/source/source_lcao/module_deepks/deepks_basic.h b/source/source_lcao/module_deepks/deepks_basic.h index 25b4fdb7e9..768b0aab1f 100644 --- a/source/source_lcao/module_deepks/deepks_basic.h +++ b/source/source_lcao/module_deepks/deepks_basic.h @@ -3,6 +3,7 @@ #ifdef __MLALGO #include "LCAO_deepks_io.h" +#include "deepks_param.h" #include "source_base/parallel_reduce.h" #include "source_base/tool_title.h" #include "source_cell/unitcell.h" @@ -32,28 +33,21 @@ void load_model(const std::string& model_file, torch::jit::script::Module& model // calculate gevdm void cal_gevdm(const int nat, - const int inlmax, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& pdm, std::vector& gevdm); /// calculate partial of energy correction to descriptors void cal_edelta_gedm(const int nat, - const int inlmax, - const int des_per_atom, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& descriptor, const std::vector& pdm, torch::jit::script::Module& model_deepks, double** gedm, double& E_delta); -void check_gedm(const int inlmax, const std::vector& inl2l, double** gedm); +void check_gedm(const DeePKS_Param& deepks_param, double** gedm); void cal_edelta_gedm_equiv(const int nat, - const int lmaxd, - const int nmaxd, - const int inlmax, - const int des_per_atom, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& descriptor, double** gedm, double& E_delta, diff --git a/source/source_lcao/module_deepks/deepks_check.cpp b/source/source_lcao/module_deepks/deepks_check.cpp index 16a07bcabb..1dfc6dbd20 100644 --- a/source/source_lcao/module_deepks/deepks_check.cpp +++ b/source/source_lcao/module_deepks/deepks_check.cpp @@ -56,10 +56,14 @@ void DeePKS_domain::check_tensor(const torch::Tensor& tensor, const std::string& ofs.close(); } - - -template void DeePKS_domain::check_tensor(const torch::Tensor& tensor, const std::string& filename, const int rank); -template void DeePKS_domain::check_tensor(const torch::Tensor& tensor, const std::string& filename, const int rank); -template void DeePKS_domain::check_tensor>(const torch::Tensor& tensor, const std::string& filename, const int rank); +template void DeePKS_domain::check_tensor(const torch::Tensor& tensor, + const std::string& filename, + const int rank); +template void DeePKS_domain::check_tensor(const torch::Tensor& tensor, + const std::string& filename, + const int rank); +template void DeePKS_domain::check_tensor>(const torch::Tensor& tensor, + const std::string& filename, + const int rank); #endif diff --git a/source/source_lcao/module_deepks/deepks_descriptor.cpp b/source/source_lcao/module_deepks/deepks_descriptor.cpp index 48b71525c4..ae8cff433e 100644 --- a/source/source_lcao/module_deepks/deepks_descriptor.cpp +++ b/source/source_lcao/module_deepks/deepks_descriptor.cpp @@ -8,25 +8,25 @@ #include "deepks_descriptor.h" #include "LCAO_deepks_io.h" // mohan add 2024-07-22 -#include "source_base/module_external/blas_connector.h" #include "source_base/constants.h" #include "source_base/libm/libm.h" +#include "source_base/module_external/blas_connector.h" #include "source_base/parallel_reduce.h" -#include "source_lcao/module_hcontainer/atom_pair.h" #include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_hcontainer/atom_pair.h" void DeePKS_domain::cal_descriptor_equiv(const int nat, - const int des_per_atom, + const DeePKS_Param& deepks_param, const std::vector& pdm, std::vector& descriptor) { ModuleBase::TITLE("DeePKS_domain", "cal_descriptor_equiv"); ModuleBase::timer::tick("DeePKS_domain", "cal_descriptor_equiv"); - assert(des_per_atom > 0); + assert(deepks_param.des_per_atom > 0); for (int iat = 0; iat < nat; iat++) { - auto tmp = torch::zeros(des_per_atom, torch::kFloat64); + auto tmp = torch::zeros(deepks_param.des_per_atom, torch::kFloat64); std::memcpy(tmp.data_ptr(), pdm[iat].data_ptr(), sizeof(double) * tmp.numel()); descriptor.push_back(tmp); } @@ -36,30 +36,28 @@ void DeePKS_domain::cal_descriptor_equiv(const int nat, // calculates descriptors from projected density matrices void DeePKS_domain::cal_descriptor(const int nat, - const int inlmax, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& pdm, - std::vector& descriptor, - const int des_per_atom = -1) + std::vector& descriptor) { ModuleBase::TITLE("DeePKS_domain", "cal_descriptor"); ModuleBase::timer::tick("DeePKS_domain", "cal_descriptor"); if (PARAM.inp.deepks_equiv) { - DeePKS_domain::cal_descriptor_equiv(nat, des_per_atom, pdm, descriptor); + DeePKS_domain::cal_descriptor_equiv(nat, deepks_param, pdm, descriptor); return; } - for (int inl = 0; inl < inlmax; ++inl) + for (int inl = 0; inl < deepks_param.inlmax; ++inl) { - const int nm = 2 * inl2l[inl] + 1; + const int nm = 2 * deepks_param.inl2l[inl] + 1; pdm[inl].requires_grad_(true); descriptor.push_back(torch::ones({nm}, torch::requires_grad(true))); } // cal descriptor - for (int inl = 0; inl < inlmax; ++inl) + for (int inl = 0; inl < deepks_param.inlmax; ++inl) { torch::Tensor vd; std::tuple d_v(descriptor[inl], vd); @@ -72,9 +70,7 @@ void DeePKS_domain::cal_descriptor(const int nat, return; } -void DeePKS_domain::check_descriptor(const int inlmax, - const int des_per_atom, - const std::vector& inl2l, +void DeePKS_domain::check_descriptor(const DeePKS_Param& deepks_param, const UnitCell& ucell, const std::string& out_dir, const std::vector& descriptor, @@ -99,13 +95,13 @@ void DeePKS_domain::check_descriptor(const int inlmax, for (int ia = 0; ia < ucell.atoms[it].na; ia++) { int iat = ucell.itia2iat(it, ia); - ofs << ucell.atoms[it].label << " atom_index " << ia + 1 << " n_descriptor " << des_per_atom - << std::endl; + ofs << ucell.atoms[it].label << " atom_index " << ia + 1 << " n_descriptor " + << deepks_param.des_per_atom << std::endl; int id = 0; - for (int inl = 0; inl < inlmax / ucell.nat; inl++) + for (int inl = 0; inl < deepks_param.inlmax / ucell.nat; inl++) { - int nm = 2 * inl2l[inl] + 1; - const int ind = iat * inlmax / ucell.nat + inl; + int nm = 2 * deepks_param.inl2l[inl] + 1; + const int ind = iat * deepks_param.inlmax / ucell.nat + inl; auto accessor = descriptor[ind].accessor(); for (int im = 0; im < nm; im++) { @@ -126,9 +122,10 @@ void DeePKS_domain::check_descriptor(const int inlmax, for (int iat = 0; iat < ucell.nat; iat++) { const int it = ucell.iat2it[iat]; - ofs << ucell.atoms[it].label << " atom_index " << iat + 1 << " n_descriptor " << des_per_atom << std::endl; + ofs << ucell.atoms[it].label << " atom_index " << iat + 1 << " n_descriptor " << deepks_param.des_per_atom + << std::endl; auto accessor = descriptor[iat].accessor(); - for (int i = 0; i < des_per_atom; i++) + for (int i = 0; i < deepks_param.des_per_atom; i++) { ofs << accessor[i] << " "; if (i % 8 == 7) diff --git a/source/source_lcao/module_deepks/deepks_descriptor.h b/source/source_lcao/module_deepks/deepks_descriptor.h index a91025503a..93e1e95559 100644 --- a/source/source_lcao/module_deepks/deepks_descriptor.h +++ b/source/source_lcao/module_deepks/deepks_descriptor.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/intarray.h" #include "source_base/timer.h" #include "source_cell/unitcell.h" @@ -29,22 +30,18 @@ namespace DeePKS_domain /// Calculates descriptors /// which are eigenvalues of pdm in blocks of I_n_l void cal_descriptor(const int nat, - const int inlmax, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& pdm, - std::vector& descriptor, - const int des_per_atom); + std::vector& descriptor); /// print descriptors based on LCAO basis -void check_descriptor(const int inlmax, - const int des_per_atom, - const std::vector& inl2l, +void check_descriptor(const DeePKS_Param& deepks_param, const UnitCell& ucell, const std::string& out_dir, const std::vector& descriptor, const int rank); void cal_descriptor_equiv(const int nat, - const int des_per_atom, + const DeePKS_Param& deepks_param, const std::vector& pdm, std::vector& descriptor); } // namespace DeePKS_domain diff --git a/source/source_lcao/module_deepks/deepks_force.cpp b/source/source_lcao/module_deepks/deepks_force.cpp index fecd13773a..293710f87e 100644 --- a/source/source_lcao/module_deepks/deepks_force.cpp +++ b/source/source_lcao/module_deepks/deepks_force.cpp @@ -17,10 +17,10 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer* dmr, const Grid_Driver& GridD, const Parallel_Orbitals& pv, const int nks, + const DeePKS_Param& deepks_param, 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) @@ -139,7 +139,7 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer* dmr, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = inl_index[T0](I0, L0, N0); + const int inl = deepks_param.inl_index[T0](I0, L0, N0); const int nm = 2 * L0 + 1; for (int m1 = 0; m1 < nm; ++m1) { @@ -269,10 +269,10 @@ template void DeePKS_domain::cal_f_delta(const hamilt::HContainer>& kvec_d, std::vector*> phialpha, double** gedm, - ModuleBase::IntArray* inl_index, ModuleBase::matrix& f_delta, const bool isstress, ModuleBase::matrix& svnl_dalpha); @@ -283,10 +283,10 @@ template void DeePKS_domain::cal_f_delta>(const hamilt::HCo const Grid_Driver& GridD, const Parallel_Orbitals& pv, const int nks, + const DeePKS_Param& deepks_param, 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); diff --git a/source/source_lcao/module_deepks/deepks_force.h b/source/source_lcao/module_deepks/deepks_force.h index a76e88e699..3837d3a6d6 100644 --- a/source/source_lcao/module_deepks/deepks_force.h +++ b/source/source_lcao/module_deepks/deepks_force.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/complexmatrix.h" #include "source_base/intarray.h" #include "source_base/matrix.h" @@ -31,10 +32,10 @@ void cal_f_delta(const hamilt::HContainer* dmr, const Grid_Driver& GridD, const Parallel_Orbitals& pv, const int nks, + const DeePKS_Param& deepks_param, 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); diff --git a/source/source_lcao/module_deepks/deepks_fpre.cpp b/source/source_lcao/module_deepks/deepks_fpre.cpp index 33a8c005be..461fbd7599 100644 --- a/source/source_lcao/module_deepks/deepks_fpre.cpp +++ b/source/source_lcao/module_deepks/deepks_fpre.cpp @@ -8,18 +8,16 @@ #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/vector3.h" -#include "source_lcao/module_hcontainer/atom_pair.h" #include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_hcontainer/atom_pair.h" /// this subroutine calculates the gradient of projected density matrices /// gdmx_m,m = d/dX sum_{mu,nu} rho_{mu,nu} template -void DeePKS_domain::cal_gdmx(const int lmaxd, - const int inlmax, - const int nks, +void DeePKS_domain::cal_gdmx(const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, std::vector*> phialpha, - const ModuleBase::IntArray* inl_index, const hamilt::HContainer* dmr, const UnitCell& ucell, const LCAO_Orbitals& orb, @@ -32,11 +30,11 @@ void DeePKS_domain::cal_gdmx(const int lmaxd, // get DS_alpha_mu and S_nu_beta int nrow = pv.nrow; - const int nm = 2 * lmaxd + 1; + const int nm = 2 * deepks_param.lmaxd + 1; // gdmx: dD/dX // \sum_{mu,nu} 2*c_mu*c_nu * // size: [3][natom][tot_Inl][2l+1][2l+1] - gdmx = torch::zeros({3, ucell.nat, inlmax, nm, nm}, torch::dtype(torch::kFloat64)); + gdmx = torch::zeros({3, ucell.nat, deepks_param.inlmax, nm, nm}, torch::dtype(torch::kFloat64)); auto accessor = gdmx.accessor(); DeePKS_domain::iterate_ad2( @@ -99,7 +97,7 @@ void DeePKS_domain::cal_gdmx(const int lmaxd, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = inl_index[ucell.iat2it[iat]](ucell.iat2ia[iat], L0, N0); + const int inl = deepks_param.inl_index[ucell.iat2it[iat]](ucell.iat2ia[iat], L0, N0); const int nm = 2 * L0 + 1; for (int m1 = 0; m1 < nm; ++m1) { @@ -136,7 +134,7 @@ void DeePKS_domain::cal_gdmx(const int lmaxd, ); #ifdef __MPI - Parallel_Reduce::reduce_all(gdmx.data_ptr(), 3 * ucell.nat * inlmax * nm * nm); + Parallel_Reduce::reduce_all(gdmx.data_ptr(), 3 * ucell.nat * deepks_param.inlmax * nm * nm); #endif ModuleBase::timer::tick("DeePKS_domain", "cal_gdmx"); return; @@ -144,9 +142,7 @@ void DeePKS_domain::cal_gdmx(const int lmaxd, // calculates gradient of descriptors from gradient of projected density matrices void DeePKS_domain::cal_gvx(const int nat, - const int inlmax, - const int des_per_atom, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& gevdm, const torch::Tensor& gdmx, torch::Tensor& gvx, @@ -161,12 +157,14 @@ void DeePKS_domain::cal_gvx(const int nat, if (rank == 0) { // make gdmx as tensor - int nlmax = inlmax / nat; + int nlmax = deepks_param.inlmax / nat; for (int nl = 0; nl < nlmax; ++nl) { - int nm = 2 * inl2l[nl] + 1; - torch::Tensor gdmx_sliced - = gdmx.slice(2, nl, inlmax, nlmax).slice(3, 0, nm, 1).slice(4, 0, nm, 1).permute({1, 0, 2, 3, 4}); + int nm = 2 * deepks_param.inl2l[nl] + 1; + torch::Tensor gdmx_sliced = gdmx.slice(2, nl, deepks_param.inlmax, nlmax) + .slice(3, 0, nm, 1) + .slice(4, 0, nm, 1) + .permute({1, 0, 2, 3, 4}); gdmr.push_back(gdmx_sliced); } @@ -190,18 +188,16 @@ void DeePKS_domain::cal_gvx(const int nat, assert(gvx.size(0) == nat); assert(gvx.size(1) == 3); assert(gvx.size(2) == nat); - assert(gvx.size(3) == des_per_atom); + assert(gvx.size(3) == deepks_param.des_per_atom); } ModuleBase::timer::tick("DeePKS_domain", "cal_gvx"); return; } -template void DeePKS_domain::cal_gdmx(const int lmaxd, - const int inlmax, - const int nks, +template void DeePKS_domain::cal_gdmx(const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, std::vector*> phialpha, - const ModuleBase::IntArray* inl_index, const hamilt::HContainer* dmr, const UnitCell& ucell, const LCAO_Orbitals& orb, @@ -209,12 +205,10 @@ template void DeePKS_domain::cal_gdmx(const int lmaxd, const Grid_Driver& GridD, torch::Tensor& gdmx); -template void DeePKS_domain::cal_gdmx>(const int lmaxd, - const int inlmax, - const int nks, +template void DeePKS_domain::cal_gdmx>(const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, std::vector*> phialpha, - const ModuleBase::IntArray* inl_index, const hamilt::HContainer* dmr, const UnitCell& ucell, const LCAO_Orbitals& orb, diff --git a/source/source_lcao/module_deepks/deepks_fpre.h b/source/source_lcao/module_deepks/deepks_fpre.h index 3ae94435b5..64c0531a2f 100644 --- a/source/source_lcao/module_deepks/deepks_fpre.h +++ b/source/source_lcao/module_deepks/deepks_fpre.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/complexmatrix.h" #include "source_base/intarray.h" #include "source_base/matrix.h" @@ -30,12 +31,10 @@ namespace DeePKS_domain // calculate the gradient of pdm with regard to atomic positions // d/dX D_{Inl,mm'} template -void cal_gdmx(const int lmaxd, - const int inlmax, - const int nks, +void cal_gdmx(const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, std::vector*> phialpha, - const ModuleBase::IntArray* inl_index, const hamilt::HContainer* dmr, const UnitCell& ucell, const LCAO_Orbitals& orb, @@ -53,9 +52,7 @@ void cal_gdmx(const int lmaxd, /// gevdm*gdmx->gvx ///---------------------------------------------------- void cal_gvx(const int nat, - const int inlmax, - const int des_per_atom, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& gevdm, const torch::Tensor& gdmx, torch::Tensor& gvx, diff --git a/source/source_lcao/module_deepks/deepks_orbital.cpp b/source/source_lcao/module_deepks/deepks_orbital.cpp index a43ff25ec3..1a0d91cafe 100644 --- a/source/source_lcao/module_deepks/deepks_orbital.cpp +++ b/source/source_lcao/module_deepks/deepks_orbital.cpp @@ -9,7 +9,6 @@ template void DeePKS_domain::cal_o_delta(const std::vector& dm_hl, const std::vector>& h_delta, - // std::vector& o_delta, ModuleBase::matrix& o_delta, const Parallel_Orbitals& pv, const int nks, @@ -48,7 +47,7 @@ void DeePKS_domain::cal_o_delta(const std::vector& dm_hl, } } Parallel_Reduce::reduce_all(o_delta_tmp); - + const double* o_delta_ptr = reinterpret_cast(&o_delta_tmp); o_delta(ik, 0) = o_delta_ptr[0]; // real part in complex case } @@ -58,7 +57,6 @@ 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, ModuleBase::matrix& o_delta, const Parallel_Orbitals& pv, const int nks, @@ -67,7 +65,6 @@ template void DeePKS_domain::cal_o_delta(const std:: template void DeePKS_domain::cal_o_delta, ModuleBase::ComplexMatrix>( const std::vector& dm_hl, const std::vector>>& h_delta, - // std::vector& o_delta, ModuleBase::matrix& o_delta, const Parallel_Orbitals& pv, const int nks, diff --git a/source/source_lcao/module_deepks/deepks_orbital.h b/source/source_lcao/module_deepks/deepks_orbital.h index 93bb3421f0..8f77e58fb1 100644 --- a/source/source_lcao/module_deepks/deepks_orbital.h +++ b/source/source_lcao/module_deepks/deepks_orbital.h @@ -26,7 +26,6 @@ namespace DeePKS_domain template void cal_o_delta(const std::vector& dm_hl, const std::vector>& h_delta, - // std::vector& o_delta, ModuleBase::matrix& o_delta, const Parallel_Orbitals& pv, const int nks, diff --git a/source/source_lcao/module_deepks/deepks_orbpre.cpp b/source/source_lcao/module_deepks/deepks_orbpre.cpp index 90cb40ea95..60d3ee16a6 100644 --- a/source/source_lcao/module_deepks/deepks_orbpre.cpp +++ b/source/source_lcao/module_deepks/deepks_orbpre.cpp @@ -7,26 +7,23 @@ #include "deepks_orbpre.h" #include "LCAO_deepks_io.h" // mohan add 2024-07-22 -#include "source_base/module_external/blas_connector.h" #include "source_base/constants.h" #include "source_base/libm/libm.h" +#include "source_base/module_external/blas_connector.h" #include "source_base/parallel_reduce.h" -#include "source_lcao/module_hcontainer/atom_pair.h" #include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_hcontainer/atom_pair.h" // calculates orbital_precalc[nks,NAt,NDscrpt] = gevdm * orbital_pdm; // orbital_pdm[nks,Inl,nm,nm] = dm_hl * overlap * overlap; template void DeePKS_domain::cal_orbital_precalc(const std::vector& dm_hl, - const int lmaxd, - const int inlmax, const int nat, const int nks, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, @@ -39,7 +36,8 @@ void DeePKS_domain::cal_orbital_precalc(const std::vector& dm_hl, const double Rcut_Alpha = orb.Alpha[0].getRcut(); torch::Tensor orbital_pdm - = torch::zeros({nks, inlmax, (2 * lmaxd + 1), (2 * lmaxd + 1)}, torch::dtype(torch::kFloat64)); + = torch::zeros({nks, deepks_param.inlmax, (2 * deepks_param.lmaxd + 1), (2 * deepks_param.lmaxd + 1)}, + torch::dtype(torch::kFloat64)); auto accessor = orbital_pdm.accessor(); for (int T0 = 0; T0 < ucell.ntype; T0++) @@ -60,7 +58,7 @@ void DeePKS_domain::cal_orbital_precalc(const std::vector& dm_hl, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = inl_index[T0](I0, L0, N0); + const int inl = deepks_param.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 @@ -242,7 +240,7 @@ void DeePKS_domain::cal_orbital_precalc(const std::vector& dm_hl, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = inl_index[T0](I0, L0, N0); + const int inl = deepks_param.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 @@ -265,19 +263,19 @@ void DeePKS_domain::cal_orbital_precalc(const std::vector& dm_hl, } } #ifdef __MPI - const int size = nks * inlmax * (2 * lmaxd + 1) * (2 * lmaxd + 1); + const int size = nks * deepks_param.inlmax * (2 * deepks_param.lmaxd + 1) * (2 * deepks_param.lmaxd + 1); Parallel_Reduce::reduce_all(orbital_pdm.data_ptr(), size); #endif // transfer orbital_pdm [nks,inl,nm,nm] to orbital_pdm_vector [nl,[nks,nat,nm,nm]] - int nlmax = inlmax / nat; + int nlmax = deepks_param.inlmax / nat; std::vector orbital_pdm_vector; for (int nl = 0; nl < nlmax; ++nl) { - int nm = 2 * inl2l[nl] + 1; + int nm = 2 * deepks_param.inl2l[nl] + 1; torch::Tensor orbital_pdm_sliced - = orbital_pdm.slice(1, nl, inlmax, nlmax).slice(2, 0, nm, 1).slice(3, 0, nm, 1); + = orbital_pdm.slice(1, nl, deepks_param.inlmax, nlmax).slice(2, 0, nm, 1).slice(3, 0, nm, 1); orbital_pdm_vector.push_back(orbital_pdm_sliced); } @@ -297,15 +295,12 @@ void DeePKS_domain::cal_orbital_precalc(const std::vector& dm_hl, template void DeePKS_domain::cal_orbital_precalc( const std::vector& dm_hl, - const int lmaxd, - const int inlmax, const int nat, const int nks, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, @@ -314,15 +309,12 @@ template void DeePKS_domain::cal_orbital_precalc( template void DeePKS_domain::cal_orbital_precalc, ModuleBase::ComplexMatrix>( const std::vector& dm_hl, - const int lmaxd, - const int inlmax, const int nat, const int nks, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, diff --git a/source/source_lcao/module_deepks/deepks_orbpre.h b/source/source_lcao/module_deepks/deepks_orbpre.h index e9ddfeb95d..ec9dc1a123 100644 --- a/source/source_lcao/module_deepks/deepks_orbpre.h +++ b/source/source_lcao/module_deepks/deepks_orbpre.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/complexmatrix.h" #include "source_base/intarray.h" #include "source_base/matrix.h" @@ -27,15 +28,12 @@ namespace DeePKS_domain template void cal_orbital_precalc(const std::vector& dm_hl, - const int lmaxd, - const int inlmax, const int nat, const int nks, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, diff --git a/source/source_lcao/module_deepks/deepks_param.h b/source/source_lcao/module_deepks/deepks_param.h new file mode 100644 index 0000000000..75e161b9da --- /dev/null +++ b/source/source_lcao/module_deepks/deepks_param.h @@ -0,0 +1,21 @@ +#ifndef LCAO_DEEPKS_PARAM +#define LCAO_DEEPKS_PARAM + +#include +namespace ModuleBase +{ +struct IntArray; +} + +struct DeePKS_Param +{ + int lmaxd = 0; + int nmaxd = 0; + int inlmax = 0; + int n_descriptor = 0; + int des_per_atom = 0; + std::vector inl2l; + ModuleBase::IntArray* inl_index = nullptr; +}; + +#endif \ No newline at end of file diff --git a/source/source_lcao/module_deepks/deepks_pdm.cpp b/source/source_lcao/module_deepks/deepks_pdm.cpp index e2f4a0033b..ae687c2c01 100644 --- a/source/source_lcao/module_deepks/deepks_pdm.cpp +++ b/source/source_lcao/module_deepks/deepks_pdm.cpp @@ -29,9 +29,7 @@ void DeePKS_domain::read_pdm(bool read_pdm_file, bool is_equiv, bool& init_pdm, const int nat, - const int inlmax, - const int lmaxd, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const Numerical_Orbital& alpha, std::vector& pdm) { @@ -46,9 +44,9 @@ void DeePKS_domain::read_pdm(bool read_pdm_file, } if (!is_equiv) { - for (int inl = 0; inl < inlmax; inl++) + for (int inl = 0; inl < deepks_param.inlmax; inl++) { - int nm = 2 * inl2l[inl] + 1; + int nm = 2 * deepks_param.inl2l[inl] + 1; auto accessor = pdm[inl].accessor(); for (int m1 = 0; m1 < nm; m1++) { @@ -65,7 +63,7 @@ void DeePKS_domain::read_pdm(bool read_pdm_file, { int pdm_size = 0; int nproj = 0; - for (int il = 0; il < lmaxd + 1; il++) + for (int il = 0; il < deepks_param.lmaxd + 1; il++) { nproj += (2 * il + 1) * alpha.getNchi(il); } @@ -177,10 +175,7 @@ void DeePKS_domain::update_dmr(const std::vector>& k // pdm_m,m'=\sum_{mu,nu} rho_{mu,nu} template void DeePKS_domain::cal_pdm(bool& init_pdm, - const int inlmax, - const int lmaxd, - const std::vector& inl2l, - const ModuleBase::IntArray* inl_index, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const hamilt::HContainer* dmr, const std::vector*> phialpha, @@ -203,9 +198,9 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, if (!PARAM.inp.deepks_equiv) { - for (int inl = 0; inl < inlmax; inl++) + for (int inl = 0; inl < deepks_param.inlmax; inl++) { - int nm = 2 * inl2l[inl] + 1; + int nm = 2 * deepks_param.inl2l[inl] + 1; pdm[inl] = torch::zeros({nm, nm}, torch::kFloat64); } } @@ -213,7 +208,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, { int pdm_size = 0; int nproj = 0; - for (int il = 0; il < lmaxd + 1; il++) + for (int il = 0; il < deepks_param.lmaxd + 1; il++) { nproj += (2 * il + 1) * orb.Alpha[0].getNchi(il); } @@ -246,7 +241,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = inl_index[T0](I0, L0, N0); + const int inl = deepks_param.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 @@ -264,7 +259,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, else { int nproj = 0; - for (int il = 0; il < lmaxd + 1; il++) + for (int il = 0; il < deepks_param.lmaxd + 1; il++) { nproj += (2 * il + 1) * orb.Alpha[0].getNchi(il); } @@ -399,7 +394,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = inl_index[T0](I0, L0, N0); + const int inl = deepks_param.inl_index[T0](I0, L0, N0); const int nm = 2 * L0 + 1; auto accessor = pdm[inl].accessor(); @@ -423,7 +418,7 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, auto accessor = pdm[iat].accessor(); int index = 0, inc = 1; int nproj = 0; - for (int il = 0; il < lmaxd + 1; il++) + for (int il = 0; il < deepks_param.lmaxd + 1; il++) { nproj += (2 * il + 1) * orb.Alpha[0].getNchi(il); } @@ -446,9 +441,9 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, } // iat #ifdef __MPI - for (int inl = 0; inl < inlmax; inl++) + for (int inl = 0; inl < deepks_param.inlmax; inl++) { - int pdm_size = (2 * inl2l[inl] + 1) * (2 * inl2l[inl] + 1); + int pdm_size = (2 * deepks_param.inl2l[inl] + 1) * (2 * deepks_param.inl2l[inl] + 1); Parallel_Reduce::reduce_all(pdm[inl].data_ptr(), pdm_size); } #endif @@ -456,15 +451,15 @@ void DeePKS_domain::cal_pdm(bool& init_pdm, return; } -void DeePKS_domain::check_pdm(const int inlmax, const std::vector& inl2l, const std::vector& pdm) +void DeePKS_domain::check_pdm(const DeePKS_Param& deepks_param, const std::vector& pdm) { const std::string file_projdm = PARAM.globalv.global_out_dir + "deepks_projdm.dat"; std::ofstream ofs(file_projdm.c_str()); ofs << std::setprecision(10); - for (int inl = 0; inl < inlmax; inl++) + for (int inl = 0; inl < deepks_param.inlmax; inl++) { - const int nm = 2 * inl2l[inl] + 1; + const int nm = 2 * deepks_param.inl2l[inl] + 1; auto accessor = pdm[inl].accessor(); for (int m1 = 0; m1 < nm; m1++) { @@ -494,10 +489,7 @@ template void DeePKS_domain::update_dmr>(const std::vector< hamilt::HContainer* dmr_deepks); template void DeePKS_domain::cal_pdm(bool& init_pdm, - const int inlmax, - const int lmaxd, - const std::vector& inl2l, - const ModuleBase::IntArray* inl_index, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const hamilt::HContainer* dmr, const std::vector*> phialpha, @@ -508,10 +500,7 @@ template void DeePKS_domain::cal_pdm(bool& init_pdm, std::vector& pdm); template void DeePKS_domain::cal_pdm>(bool& init_pdm, - const int inlmax, - const int lmaxd, - const std::vector& inl2l, - const ModuleBase::IntArray* inl_index, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const hamilt::HContainer* dmr, const std::vector*> phialpha, diff --git a/source/source_lcao/module_deepks/deepks_pdm.h b/source/source_lcao/module_deepks/deepks_pdm.h index f7b90b50c1..4b39547576 100644 --- a/source/source_lcao/module_deepks/deepks_pdm.h +++ b/source/source_lcao/module_deepks/deepks_pdm.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/complexmatrix.h" #include "source_base/matrix.h" #include "source_base/timer.h" @@ -33,9 +34,7 @@ void read_pdm(bool read_pdm_file, bool is_equiv, bool& init_pdm, const int nat, - const int inlmax, - const int lmaxd, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const Numerical_Orbital& alpha, std::vector& pdm); @@ -55,10 +54,7 @@ void update_dmr(const std::vector>& kvec_d, // - Relax/Cell-Relax/MD calculation, non-first step will use the convergence pdm from the last step as initial pdm template void cal_pdm(bool& init_pdm, - const int inlmax, - const int lmaxd, - const std::vector& inl2l, - const ModuleBase::IntArray* inl_index, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const hamilt::HContainer* dmr, const std::vector*> phialpha, @@ -68,7 +64,7 @@ void cal_pdm(bool& init_pdm, const Parallel_Orbitals& pv, std::vector& pdm); -void check_pdm(const int inlmax, const std::vector& inl2l, const std::vector& pdm); +void check_pdm(const DeePKS_Param& deepks_param, const std::vector& pdm); } // namespace DeePKS_domain #endif diff --git a/source/source_lcao/module_deepks/deepks_spre.cpp b/source/source_lcao/module_deepks/deepks_spre.cpp index eb0281146c..82fd8a5303 100644 --- a/source/source_lcao/module_deepks/deepks_spre.cpp +++ b/source/source_lcao/module_deepks/deepks_spre.cpp @@ -8,19 +8,17 @@ #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/vector3.h" -#include "source_lcao/module_hcontainer/atom_pair.h" #include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_hcontainer/atom_pair.h" /// this subroutine calculates the gradient of PDM wrt strain tensor: /// gdmepsl = d/d\epsilon_{ab} * /// sum_{mu,nu} rho_{mu,nu} template -void DeePKS_domain::cal_gdmepsl(const int lmaxd, - const int inlmax, - const int nks, +void DeePKS_domain::cal_gdmepsl(const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, std::vector*> phialpha, - const ModuleBase::IntArray* inl_index, const hamilt::HContainer* dmr, const UnitCell& ucell, const LCAO_Orbitals& orb, @@ -33,10 +31,10 @@ void DeePKS_domain::cal_gdmepsl(const int lmaxd, // get DS_alpha_mu and S_nu_beta int nrow = pv.nrow; - const int nm = 2 * lmaxd + 1; + const int nm = 2 * deepks_param.lmaxd + 1; // gdmepsl: dD/d\epsilon_{\alpha\beta} // size: [6][tot_Inl][2l+1][2l+1] - gdmepsl = torch::zeros({6, inlmax, nm, nm}, torch::dtype(torch::kFloat64)); + gdmepsl = torch::zeros({6, deepks_param.inlmax, nm, nm}, torch::dtype(torch::kFloat64)); auto accessor = gdmepsl.accessor(); DeePKS_domain::iterate_ad2( @@ -111,7 +109,7 @@ void DeePKS_domain::cal_gdmepsl(const int lmaxd, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = inl_index[ucell.iat2it[iat]](ucell.iat2ia[iat], L0, N0); + const int inl = deepks_param.inl_index[ucell.iat2it[iat]](ucell.iat2ia[iat], L0, N0); const int nm = 2 * L0 + 1; for (int m1 = 0; m1 < nm; ++m1) { @@ -147,7 +145,7 @@ void DeePKS_domain::cal_gdmepsl(const int lmaxd, ); #ifdef __MPI - Parallel_Reduce::reduce_all(gdmepsl.data_ptr(), 6 * inlmax * nm * nm); + Parallel_Reduce::reduce_all(gdmepsl.data_ptr(), 6 * deepks_param.inlmax * nm * nm); #endif ModuleBase::timer::tick("DeePKS_domain", "cal_gdmepsl"); return; @@ -156,9 +154,7 @@ void DeePKS_domain::cal_gdmepsl(const int lmaxd, // calculates stress of descriptors from gradient of projected density matrices // gv_epsl:d(d)/d\epsilon_{\alpha\beta}, [natom][6][des_per_atom] void DeePKS_domain::cal_gvepsl(const int nat, - const int inlmax, - const int des_per_atom, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& gevdm, const torch::Tensor& gdmepsl, torch::Tensor& gvepsl, @@ -172,11 +168,12 @@ void DeePKS_domain::cal_gvepsl(const int nat, if (rank == 0) { // make gdmepsl as tensor - int nlmax = inlmax / nat; + int nlmax = deepks_param.inlmax / nat; for (int nl = 0; nl < nlmax; ++nl) { - int nm = 2 * inl2l[nl] + 1; - torch::Tensor gdmepsl_sliced = gdmepsl.slice(1, nl, inlmax, nlmax).slice(2, 0, nm, 1).slice(3, 0, nm, 1); + int nm = 2 * deepks_param.inl2l[nl] + 1; + torch::Tensor gdmepsl_sliced + = gdmepsl.slice(1, nl, deepks_param.inlmax, nlmax).slice(2, 0, nm, 1).slice(3, 0, nm, 1); gdmepsl_vector.push_back(gdmepsl_sliced); } assert(gdmepsl_vector.size() == nlmax); @@ -197,19 +194,17 @@ void DeePKS_domain::cal_gvepsl(const int nat, gvepsl = torch::cat(gvepsl_vector, -1); assert(gvepsl.size(0) == 6); assert(gvepsl.size(1) == nat); - assert(gvepsl.size(2) == des_per_atom); + assert(gvepsl.size(2) == deepks_param.des_per_atom); } ModuleBase::timer::tick("DeePKS_domain", "cal_gvepsl"); return; } -template void DeePKS_domain::cal_gdmepsl(const int lmaxd, - const int inlmax, - const int nks, +template void DeePKS_domain::cal_gdmepsl(const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, std::vector*> phialpha, - const ModuleBase::IntArray* inl_index, const hamilt::HContainer* dmr, const UnitCell& ucell, const LCAO_Orbitals& orb, @@ -217,12 +212,10 @@ template void DeePKS_domain::cal_gdmepsl(const int lmaxd, const Grid_Driver& GridD, torch::Tensor& gdmepsl); -template void DeePKS_domain::cal_gdmepsl>(const int lmaxd, - const int inlmax, - const int nks, +template void DeePKS_domain::cal_gdmepsl>(const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, std::vector*> phialpha, - const ModuleBase::IntArray* inl_index, const hamilt::HContainer* dmr, const UnitCell& ucell, const LCAO_Orbitals& orb, diff --git a/source/source_lcao/module_deepks/deepks_spre.h b/source/source_lcao/module_deepks/deepks_spre.h index 54d9f6b05e..759af70209 100644 --- a/source/source_lcao/module_deepks/deepks_spre.h +++ b/source/source_lcao/module_deepks/deepks_spre.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/complexmatrix.h" #include "source_base/intarray.h" #include "source_base/matrix.h" @@ -31,24 +32,19 @@ namespace DeePKS_domain // calculate the gradient of pdm with regard to atomic virial stress tensor // d/d\epsilon D_{Inl,mm'} template -void cal_gdmepsl( // const ModuleBase::matrix& dm, - const int lmaxd, - const int inlmax, - const int nks, - const std::vector>& kvec_d, - std::vector*> phialpha, - const ModuleBase::IntArray* inl_index, - const hamilt::HContainer* dmr, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const Parallel_Orbitals& pv, - const Grid_Driver& GridD, - torch::Tensor& gdmepsl); +void cal_gdmepsl(const int nks, + const DeePKS_Param& deepks_param, + const std::vector>& kvec_d, + std::vector*> phialpha, + const hamilt::HContainer* dmr, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const Parallel_Orbitals& pv, + const Grid_Driver& GridD, + torch::Tensor& gdmepsl); void cal_gvepsl(const int nat, - const int inlmax, - const int des_per_atom, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector& gevdm, const torch::Tensor& gdmepsl, torch::Tensor& gvepsl, diff --git a/source/source_lcao/module_deepks/deepks_vdpre.cpp b/source/source_lcao/module_deepks/deepks_vdpre.cpp index 0dad1c0bef..20ee69a499 100644 --- a/source/source_lcao/module_deepks/deepks_vdpre.cpp +++ b/source/source_lcao/module_deepks/deepks_vdpre.cpp @@ -10,27 +10,24 @@ #include "LCAO_deepks_io.h" // mohan add 2024-07-22 #include "deepks_iterate.h" -#include "source_base/module_external/blas_connector.h" #include "source_base/constants.h" #include "source_base/libm/libm.h" +#include "source_base/module_external/blas_connector.h" #include "source_base/parallel_reduce.h" -#include "source_lcao/module_hcontainer/atom_pair.h" #include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_hcontainer/atom_pair.h" // calculates v_delta_precalc[nks,nlocal,nlocal,NAt,NDscrpt] = gevdm * v_delta_pdm; // v_delta_pdm[nks,nlocal,nlocal,Inl,nm,nm] = overlap * overlap; // for deepks_v_delta = 1 template void DeePKS_domain::cal_v_delta_precalc(const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, @@ -46,8 +43,9 @@ void DeePKS_domain::cal_v_delta_precalc(const int nlocal, using TK_tensor = typename std::conditional>::value, c10::complex, TK>::type; - torch::Tensor v_delta_pdm - = torch::zeros({nks, nlocal, nlocal, inlmax, (2 * lmaxd + 1), (2 * lmaxd + 1)}, torch::dtype(dtype)); + torch::Tensor v_delta_pdm = torch::zeros( + {nks, nlocal, nlocal, deepks_param.inlmax, (2 * deepks_param.lmaxd + 1), (2 * deepks_param.lmaxd + 1)}, + torch::dtype(dtype)); auto accessor = v_delta_pdm.accessor(); DeePKS_domain::iterate_ad2( @@ -112,14 +110,14 @@ void DeePKS_domain::cal_v_delta_precalc(const int nlocal, { for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) { - const int inl = inl_index[T0](I0, L0, N0); + const int inl = deepks_param.inl_index[T0](I0, L0, N0); const int nm = 2 * L0 + 1; for (int m1 = 0; m1 < nm; ++m1) // nm = 1 for s, 3 for p, 5 for d { for (int m2 = 0; m2 < nm; ++m2) // nm = 1 for s, 3 for p, 5 for d { - TK tmp = overlap_1->get_value(iw1, ib + m1) - * overlap_2->get_value(iw2, ib + m2) * *kpase_ptr; + TK tmp = overlap_1->get_value(iw1, ib + m1) * overlap_2->get_value(iw2, ib + m2) + * *kpase_ptr; TK_tensor tmp_tensor = TK_tensor(tmp); accessor[ik][iw1_all][iw2_all][inl][m1][m2] += tmp_tensor; } @@ -133,20 +131,21 @@ void DeePKS_domain::cal_v_delta_precalc(const int nlocal, } ); #ifdef __MPI - const int size = nks * nlocal * nlocal * inlmax * (2 * lmaxd + 1) * (2 * lmaxd + 1); + const int size + = nks * nlocal * nlocal * deepks_param.inlmax * (2 * deepks_param.lmaxd + 1) * (2 * deepks_param.lmaxd + 1); TK_tensor* data_tensor_ptr = v_delta_pdm.data_ptr(); TK* data_ptr = reinterpret_cast(data_tensor_ptr); Parallel_Reduce::reduce_all(data_ptr, size); #endif // transfer v_delta_pdm to v_delta_pdm_vector - int nlmax = inlmax / nat; + int nlmax = deepks_param.inlmax / nat; std::vector v_delta_pdm_vector; for (int nl = 0; nl < nlmax; ++nl) { - int nm = 2 * inl2l[nl] + 1; + int nm = 2 * deepks_param.inl2l[nl] + 1; torch::Tensor v_delta_pdm_sliced - = v_delta_pdm.slice(3, nl, inlmax, nlmax).slice(4, 0, nm, 1).slice(5, 0, nm, 1); + = v_delta_pdm.slice(3, nl, deepks_param.inlmax, nlmax).slice(4, 0, nm, 1).slice(5, 0, nm, 1); v_delta_pdm_vector.push_back(v_delta_pdm_sliced); } @@ -173,10 +172,9 @@ void DeePKS_domain::cal_v_delta_precalc(const int nlocal, // prepare_phialpha and prepare_gevdm for deepks_v_delta = 2 template void DeePKS_domain::prepare_phialpha(const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const UnitCell& ucell, @@ -190,8 +188,8 @@ void DeePKS_domain::prepare_phialpha(const int nlocal, constexpr torch::Dtype dtype = std::is_same::value ? torch::kFloat64 : torch::kComplexDouble; using TK_tensor = typename std::conditional>::value, c10::complex, TK>::type; - int nlmax = inlmax / nat; - int mmax = 2 * lmaxd + 1; + int nlmax = deepks_param.inlmax / nat; + int mmax = 2 * deepks_param.lmaxd + 1; phialpha_out = torch::zeros({nat, nlmax, nks, nlocal, mmax}, dtype); auto accessor = phialpha_out.accessor(); @@ -268,16 +266,15 @@ void DeePKS_domain::prepare_phialpha(const int nlocal, } void DeePKS_domain::prepare_gevdm(const int nat, - const int lmaxd, - const int inlmax, + const DeePKS_Param& deepks_param, const LCAO_Orbitals& orb, const std::vector& gevdm_in, torch::Tensor& gevdm_out) { ModuleBase::TITLE("DeePKS_domain", "prepare_gevdm"); ModuleBase::timer::tick("DeePKS_domain", "prepare_gevdm"); - int nlmax = inlmax / nat; - int mmax = 2 * lmaxd + 1; + int nlmax = deepks_param.inlmax / nat; + int mmax = 2 * deepks_param.lmaxd + 1; gevdm_out = torch::zeros({nat, nlmax, mmax, mmax, mmax}, torch::TensorOptions().dtype(torch::kFloat64)); std::vector gevdm_out_vector; @@ -295,15 +292,12 @@ void DeePKS_domain::prepare_gevdm(const int nat, } template void DeePKS_domain::cal_v_delta_precalc(const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, @@ -311,15 +305,12 @@ template void DeePKS_domain::cal_v_delta_precalc(const int nlocal, torch::Tensor& v_delta_precalc); template void DeePKS_domain::cal_v_delta_precalc>( const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, @@ -327,10 +318,9 @@ template void DeePKS_domain::cal_v_delta_precalc>( torch::Tensor& v_delta_precalc); template void DeePKS_domain::prepare_phialpha(const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const UnitCell& ucell, @@ -341,10 +331,9 @@ template void DeePKS_domain::prepare_phialpha(const int nlocal, template void DeePKS_domain::prepare_phialpha>( const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const UnitCell& ucell, diff --git a/source/source_lcao/module_deepks/deepks_vdpre.h b/source/source_lcao/module_deepks/deepks_vdpre.h index 90d2e887f6..bf096a4267 100644 --- a/source/source_lcao/module_deepks/deepks_vdpre.h +++ b/source/source_lcao/module_deepks/deepks_vdpre.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/complexmatrix.h" #include "source_base/intarray.h" #include "source_base/matrix.h" @@ -32,15 +33,12 @@ namespace DeePKS_domain // calculates v_delta_precalc template void cal_v_delta_precalc(const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, @@ -51,10 +49,9 @@ void cal_v_delta_precalc(const int nlocal, // prepare phialpha for outputting npy file template void prepare_phialpha(const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const UnitCell& ucell, @@ -65,8 +62,7 @@ void prepare_phialpha(const int nlocal, // prepare gevdm for outputting npy file void prepare_gevdm(const int nat, - const int lmaxd, - const int inlmax, + const DeePKS_Param& deepks_param, const LCAO_Orbitals& orb, const std::vector& gevdm_in, torch::Tensor& gevdm_out); diff --git a/source/source_lcao/module_deepks/deepks_vdrpre.cpp b/source/source_lcao/module_deepks/deepks_vdrpre.cpp index 6a550158ae..e58f10b2e2 100644 --- a/source/source_lcao/module_deepks/deepks_vdrpre.cpp +++ b/source/source_lcao/module_deepks/deepks_vdrpre.cpp @@ -1,4 +1,4 @@ -// prepare_phialpha_r : prepare phialpha_r for outputting npy file +// prepare_phialpha_iRmat : prepare phialpha_r and iR_mat for outputting npy file #ifdef __MLALGO @@ -6,34 +6,29 @@ #include "LCAO_deepks_io.h" // mohan add 2024-07-22 #include "deepks_iterate.h" -#include "source_base/module_external/blas_connector.h" #include "source_base/constants.h" #include "source_base/libm/libm.h" +#include "source_base/module_external/blas_connector.h" #include "source_base/parallel_reduce.h" -#include "source_lcao/module_hcontainer/atom_pair.h" #include "source_io/module_parameter/parameter.h" +#include "source_lcao/module_hcontainer/atom_pair.h" -void DeePKS_domain::prepare_phialpha_r(const int nlocal, - const int lmaxd, - const int inlmax, - const int nat, - const int R_size, - const std::vector*> phialpha, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const Parallel_Orbitals& pv, - const Grid_Driver& GridD, - torch::Tensor& phialpha_r_out) +void DeePKS_domain::prepare_phialpha_iRmat(const int nlocal, + const int R_size, + const DeePKS_Param& deepks_param, + const std::vector*> phialpha, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const Grid_Driver& GridD, + torch::Tensor& overlap, + torch::Tensor& iRmat) { - ModuleBase::TITLE("DeePKS_domain", "prepare_phialpha_r"); - ModuleBase::timer::tick("DeePKS_domain", "prepare_phialpha_r"); + ModuleBase::TITLE("DeePKS_domain", "prepare_phialpha_iRmat"); + ModuleBase::timer::tick("DeePKS_domain", "prepare_phialpha_iRmat"); constexpr torch::Dtype dtype = torch::kFloat64; - int nlmax = inlmax / nat; - int mmax = 2 * lmaxd + 1; - - phialpha_r_out = torch::zeros({R_size, R_size, R_size, nat, nlmax, nlocal, mmax}, dtype); - auto accessor = phialpha_r_out.accessor(); + // get the maximum nnmax + std::vector nnmax_vec(ucell.nat, 0); DeePKS_domain::iterate_ad1( ucell, GridD, @@ -42,7 +37,7 @@ void DeePKS_domain::prepare_phialpha_r(const int nlocal, [&](const int iat, const ModuleBase::Vector3& tau0, const int ibt, - const ModuleBase::Vector3& tau, + const ModuleBase::Vector3& tau1, const int start, const int nw_tot, ModuleBase::Vector3 dR) @@ -51,65 +46,62 @@ void DeePKS_domain::prepare_phialpha_r(const int nlocal, { return; // to next loop } + nnmax_vec[iat]++; + } + ); + + int nnmax = *std::max_element(nnmax_vec.begin(), nnmax_vec.end()); + overlap = torch::zeros({ucell.nat, nnmax, nlocal, deepks_param.des_per_atom}, dtype); + torch::Tensor dRmat_tmp = torch::zeros({ucell.nat, nnmax, 3}, torch::kInt32); + auto overlap_accessor = overlap.accessor(); + auto dRmat_accessor = dRmat_tmp.accessor(); - // middle loop : all atomic basis on the adjacent atom ad - for (int iw1 = 0; iw1 < nw_tot; ++iw1) + std::fill(nnmax_vec.begin(), nnmax_vec.end(), 0); + DeePKS_domain::iterate_ad1( + ucell, + GridD, + orb, + false, // no trace_alpha + [&](const int iat, + const ModuleBase::Vector3& tau0, + const int ibt, + const ModuleBase::Vector3& tau1, + const int start, + const int nw_tot, + ModuleBase::Vector3 dR) + { + hamilt::BaseMatrix* overlap_mat = phialpha[0]->find_matrix(iat, ibt, dR); + if (overlap_mat == nullptr) { - const int iw1_all = start + iw1; - const int iw1_local = pv.global2local_row(iw1_all); - const int iw2_local = pv.global2local_col(iw1_all); - if (iw1_local < 0 || iw2_local < 0) - { - continue; - } - hamilt::BaseMatrix* overlap = phialpha[0]->find_matrix(iat, ibt, dR); - const int iR = phialpha[0]->find_R(dR); + return; // to next loop + } + dRmat_accessor[iat][nnmax_vec[iat]][0] = dR.x; + dRmat_accessor[iat][nnmax_vec[iat]][1] = dR.y; + dRmat_accessor[iat][nnmax_vec[iat]][2] = dR.z; - int ib = 0; - int nl = 0; - for (int L0 = 0; L0 <= orb.Alpha[0].getLmax(); ++L0) + for (int ix = 0; ix < nw_tot; ix++) + { + for (int iy = 0; iy < deepks_param.des_per_atom; iy++) { - for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) - { - const int nm = 2 * L0 + 1; - for (int m1 = 0; m1 < nm; ++m1) // nm = 1 for s, 3 for p, 5 for d - { - int iRx = DeePKS_domain::mapping_R(dR.x); - int iRy = DeePKS_domain::mapping_R(dR.y); - int iRz = DeePKS_domain::mapping_R(dR.z); - accessor[iRx][iRy][iRz][iat][nl][iw1_all][m1] - += overlap->get_value(iw1, ib + m1); - } - ib += nm; - nl++; - } + overlap_accessor[iat][nnmax_vec[iat]][start + ix][iy] = overlap_mat->get_value(ix, iy); } - } // end iw + } + nnmax_vec[iat]++; } ); - -#ifdef __MPI - int size = R_size * R_size * R_size * nat * nlmax * nlocal * mmax; - double* data_ptr = phialpha_r_out.data_ptr(); - Parallel_Reduce::reduce_all(data_ptr, size); - -#endif - - ModuleBase::timer::tick("DeePKS_domain", "prepare_phialpha_r"); + iRmat = mapping_R(dRmat_tmp.unsqueeze(1) - dRmat_tmp.unsqueeze(2)); + ModuleBase::timer::tick("DeePKS_domain", "prepare_phialpha_iRmat"); return; } void DeePKS_domain::cal_vdr_precalc(const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, const int R_size, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, @@ -119,106 +111,111 @@ void DeePKS_domain::cal_vdr_precalc(const int nlocal, ModuleBase::TITLE("DeePKS_domain", "calc_vdr_precalc"); ModuleBase::timer::tick("DeePKS_domain", "calc_vdr_precalc"); - torch::Tensor vdr_pdm - = torch::zeros({R_size, R_size, R_size, nlocal, nlocal, inlmax, (2 * lmaxd + 1), (2 * lmaxd + 1)}, - torch::TensorOptions().dtype(torch::kFloat64)); + torch::Tensor vdr_pdm = torch::zeros({R_size, + R_size, + R_size, + nlocal, + nlocal, + deepks_param.inlmax, + (2 * deepks_param.lmaxd + 1), + (2 * deepks_param.lmaxd + 1)}, + torch::TensorOptions().dtype(torch::kFloat64)); auto accessor = vdr_pdm.accessor(); - DeePKS_domain::iterate_ad2( - ucell, - GridD, - orb, - false, // no trace_alpha - [&](const int iat, - const ModuleBase::Vector3& tau0, - const int ibt1, - const ModuleBase::Vector3& tau1, - const int start1, - const int nw1_tot, - ModuleBase::Vector3 dR1, - const int ibt2, - const ModuleBase::Vector3& tau2, - const int start2, - const int nw2_tot, - ModuleBase::Vector3 dR2) - { - const int T0 = ucell.iat2it[iat]; - const int I0 = ucell.iat2ia[iat]; - if (phialpha[0]->find_matrix(iat, ibt1, dR1.x, dR1.y, dR1.z) == nullptr - || phialpha[0]->find_matrix(iat, ibt2, dR2.x, dR2.y, dR2.z) == nullptr) - { - return; // to next loop - } + DeePKS_domain::iterate_ad2(ucell, + GridD, + orb, + false, // no trace_alpha + [&](const int iat, + const ModuleBase::Vector3& tau0, + const int ibt1, + const ModuleBase::Vector3& tau1, + const int start1, + const int nw1_tot, + ModuleBase::Vector3 dR1, + const int ibt2, + const ModuleBase::Vector3& tau2, + const int start2, + const int nw2_tot, + ModuleBase::Vector3 dR2) { + const int T0 = ucell.iat2it[iat]; + const int I0 = ucell.iat2ia[iat]; + if (phialpha[0]->find_matrix(iat, ibt1, dR1.x, dR1.y, dR1.z) == nullptr + || phialpha[0]->find_matrix(iat, ibt2, dR2.x, dR2.y, dR2.z) == nullptr) + { + return; // to next loop + } - hamilt::BaseMatrix* overlap_1 = phialpha[0]->find_matrix(iat, ibt1, dR1); - hamilt::BaseMatrix* overlap_2 = phialpha[0]->find_matrix(iat, ibt2, dR2); - assert(overlap_1->get_col_size() == overlap_2->get_col_size()); - ModuleBase::Vector3 dR = dR2 - dR1; - int iRx = DeePKS_domain::mapping_R(dR.x); - int iRy = DeePKS_domain::mapping_R(dR.y); - int iRz = DeePKS_domain::mapping_R(dR.z); - // Make sure the index is in range we need to save - if (iRx >= R_size || iRy >= R_size || iRz >= R_size) - { - return; // to next loop - } + hamilt::BaseMatrix* overlap_1 = phialpha[0]->find_matrix(iat, ibt1, dR1); + hamilt::BaseMatrix* overlap_2 = phialpha[0]->find_matrix(iat, ibt2, dR2); + assert(overlap_1->get_col_size() == overlap_2->get_col_size()); + ModuleBase::Vector3 dR = dR2 - dR1; + int iRx = DeePKS_domain::mapping_R(dR.x); + int iRy = DeePKS_domain::mapping_R(dR.y); + int iRz = DeePKS_domain::mapping_R(dR.z); + // Make sure the index is in range we need to save + if (iRx >= R_size || iRy >= R_size || iRz >= R_size) + { + return; // to next loop + } - for (int iw1 = 0; iw1 < nw1_tot; ++iw1) - { - const int iw1_all = start1 + iw1; // this is \mu - const int iw1_local = pv.global2local_row(iw1_all); - if (iw1_local < 0) - { - continue; - } - for (int iw2 = 0; iw2 < nw2_tot; ++iw2) - { - const int iw2_all = start2 + iw2; // this is \nu - const int iw2_local = pv.global2local_col(iw2_all); - if (iw2_local < 0) - { - continue; - } + for (int iw1 = 0; iw1 < nw1_tot; ++iw1) + { + const int iw1_all = start1 + iw1; // this is \mu + const int iw1_local = pv.global2local_row(iw1_all); + if (iw1_local < 0) + { + continue; + } + for (int iw2 = 0; iw2 < nw2_tot; ++iw2) + { + const int iw2_all = start2 + iw2; // this is \nu + const int iw2_local = pv.global2local_col(iw2_all); + if (iw2_local < 0) + { + continue; + } - int ib = 0; - for (int L0 = 0; L0 <= orb.Alpha[0].getLmax(); ++L0) - { - for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) - { - const int inl = inl_index[T0](I0, L0, N0); - const int nm = 2 * L0 + 1; + int ib = 0; + for (int L0 = 0; L0 <= orb.Alpha[0].getLmax(); ++L0) + { + for (int N0 = 0; N0 < orb.Alpha[0].getNchi(L0); ++N0) + { + const int inl = deepks_param.inl_index[T0](I0, L0, N0); + const int nm = 2 * L0 + 1; - for (int m1 = 0; m1 < nm; ++m1) // nm = 1 for s, 3 for p, 5 for d - { - for (int m2 = 0; m2 < nm; ++m2) // nm = 1 for s, 3 for p, 5 for d - { - double tmp = overlap_1->get_value(iw1, ib + m1) - * overlap_2->get_value(iw2, ib + m2); - accessor[iRx][iRy][iRz][iw1_all][iw2_all][inl][m1][m2] - += tmp; - } - } - ib += nm; - } - } - } // iw2 - } // iw1 - } - ); + for (int m1 = 0; m1 < nm; ++m1) // nm = 1 for s, 3 for p, 5 for d + { + for (int m2 = 0; m2 < nm; ++m2) // nm = 1 for s, 3 for p, 5 for d + { + double tmp = overlap_1->get_value(iw1, ib + m1) + * overlap_2->get_value(iw2, ib + m2); + accessor[iRx][iRy][iRz][iw1_all][iw2_all][inl][m1][m2] + += tmp; + } + } + ib += nm; + } + } + } // iw2 + } // iw1 + }); #ifdef __MPI - const int size = R_size * R_size * R_size * nlocal * nlocal * inlmax * (2 * lmaxd + 1) * (2 * lmaxd + 1); + const int size = R_size * R_size * R_size * nlocal * nlocal * deepks_param.inlmax * (2 * deepks_param.lmaxd + 1) + * (2 * deepks_param.lmaxd + 1); double* data_ptr = vdr_pdm.data_ptr(); Parallel_Reduce::reduce_all(data_ptr, size); #endif // transfer v_delta_pdm to v_delta_pdm_vector - int nlmax = inlmax / nat; + int nlmax = deepks_param.inlmax / nat; std::vector vdr_pdm_vector; for (int nl = 0; nl < nlmax; ++nl) { - int nm = 2 * inl2l[nl] + 1; - torch::Tensor vdr_pdm_sliced = vdr_pdm.slice(5, nl, inlmax, nlmax).slice(6, 0, nm, 1).slice(7, 0, nm, 1); + int nm = 2 * deepks_param.inl2l[nl] + 1; + torch::Tensor vdr_pdm_sliced + = vdr_pdm.slice(5, nl, deepks_param.inlmax, nlmax).slice(6, 0, nm, 1).slice(7, 0, nm, 1); vdr_pdm_vector.push_back(vdr_pdm_sliced); } @@ -255,6 +252,15 @@ int DeePKS_domain::mapping_R(int R) return R_index; } +torch::Tensor DeePKS_domain::mapping_R(const torch::Tensor& R_tensor) +{ + auto R = R_tensor.to(torch::kInt32); + auto pos = R > 0; + auto twoR_minus1 = R * 2 - 1; + auto neg_minus2R = -2 * R; + return at::where(pos, twoR_minus1, neg_minus2R); +} + template int DeePKS_domain::get_R_size(const hamilt::HContainer& hcontainer) { diff --git a/source/source_lcao/module_deepks/deepks_vdrpre.h b/source/source_lcao/module_deepks/deepks_vdrpre.h index cf1a4a1b61..fc7975ce49 100644 --- a/source/source_lcao/module_deepks/deepks_vdrpre.h +++ b/source/source_lcao/module_deepks/deepks_vdrpre.h @@ -3,6 +3,7 @@ #ifdef __MLALGO +#include "deepks_param.h" #include "source_base/complexmatrix.h" #include "source_base/intarray.h" #include "source_base/matrix.h" @@ -28,29 +29,25 @@ namespace DeePKS_domain // for deepks_v_delta = -1 // calculates v_delta_r_precalc -void prepare_phialpha_r(const int nlocal, - const int lmaxd, - const int inlmax, - const int nat, - const int R_size, - const std::vector*> phialpha, - const UnitCell& ucell, - const LCAO_Orbitals& orb, - const Parallel_Orbitals& pv, - const Grid_Driver& GridD, - torch::Tensor& phialpha_r_out); + +void prepare_phialpha_iRmat(const int nlocal, + const int R_size, + const DeePKS_Param& deepks_param, + const std::vector*> phialpha, + const UnitCell& ucell, + const LCAO_Orbitals& orb, + const Grid_Driver& GridD, + torch::Tensor& overlap, + torch::Tensor& iRmat); void cal_vdr_precalc(const int nlocal, - const int lmaxd, - const int inlmax, const int nat, const int nks, const int R_size, - const std::vector& inl2l, + const DeePKS_Param& deepks_param, const std::vector>& kvec_d, const std::vector*> phialpha, const std::vector gevdm, - const ModuleBase::IntArray* inl_index, const UnitCell& ucell, const LCAO_Orbitals& orb, const Parallel_Orbitals& pv, @@ -58,6 +55,7 @@ void cal_vdr_precalc(const int nlocal, torch::Tensor& vdr_precalc); int mapping_R(int R); +torch::Tensor mapping_R(const torch::Tensor& R_tensor); template int get_R_size(const hamilt::HContainer& hcontainer); diff --git a/source/source_lcao/module_deepks/test/LCAO_deepks_test.cpp b/source/source_lcao/module_deepks/test/LCAO_deepks_test.cpp index 81b14f25bd..68718bd9b6 100644 --- a/source/source_lcao/module_deepks/test/LCAO_deepks_test.cpp +++ b/source/source_lcao/module_deepks/test/LCAO_deepks_test.cpp @@ -146,10 +146,7 @@ void test_deepks::check_pdm() Test_Deepks::GridD, this->ld.dm_r); DeePKS_domain::cal_pdm(this->ld.init_pdm, - this->ld.inlmax, - this->ld.lmaxd, - this->ld.inl2l, - this->ld.inl_index, + this->ld.deepks_param, kv.kvec_d, this->ld.dm_r, this->ld.phialpha, @@ -158,32 +155,25 @@ void test_deepks::check_pdm() Test_Deepks::GridD, ParaO, this->ld.pdm); - DeePKS_domain::check_pdm(this->ld.inlmax, this->ld.inl2l, this->ld.pdm); + DeePKS_domain::check_pdm(this->ld.deepks_param, this->ld.pdm); this->compare_with_ref("deepks_projdm.dat", "pdm_ref.dat"); } template void test_deepks::check_descriptor(std::vector& descriptor) { - DeePKS_domain::cal_descriptor(ucell.nat, - this->ld.inlmax, - this->ld.inl2l, - this->ld.pdm, - descriptor, - this->ld.des_per_atom); - DeePKS_domain::check_descriptor(this->ld.inlmax, this->ld.des_per_atom, this->ld.inl2l, ucell, "./", descriptor, 0); + DeePKS_domain::cal_descriptor(ucell.nat, this->ld.deepks_param, this->ld.pdm, descriptor); + DeePKS_domain::check_descriptor(this->ld.deepks_param, ucell, "./", descriptor, 0); this->compare_with_ref("deepks_desc.dat", "descriptor_ref.dat"); } template void test_deepks::check_gdmx(torch::Tensor& gdmx) { - DeePKS_domain::cal_gdmx(this->ld.lmaxd, - this->ld.inlmax, - kv.nkstot, + DeePKS_domain::cal_gdmx(kv.nkstot, + this->ld.deepks_param, kv.kvec_d, this->ld.phialpha, - this->ld.inl_index, this->ld.dm_r, ucell, ORB, @@ -198,9 +188,9 @@ template void test_deepks::check_gvx(torch::Tensor& gdmx) { std::vector gevdm; - DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl2l, this->ld.pdm, gevdm); + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.deepks_param, this->ld.pdm, gevdm); torch::Tensor gvx; - DeePKS_domain::cal_gvx(ucell.nat, this->ld.inlmax, this->ld.des_per_atom, this->ld.inl2l, gevdm, gdmx, gvx, 0); + DeePKS_domain::cal_gvx(ucell.nat, this->ld.deepks_param, gevdm, gdmx, gvx, 0); DeePKS_domain::check_tensor(gvx, "gvx.dat", 0); // 0 for rank this->compare_with_ref("gvx.dat", "gvx_ref.dat"); } @@ -208,12 +198,10 @@ void test_deepks::check_gvx(torch::Tensor& gdmx) template void test_deepks::check_gdmepsl(torch::Tensor& gdmepsl) { - DeePKS_domain::cal_gdmepsl(this->ld.lmaxd, - this->ld.inlmax, - kv.nkstot, + DeePKS_domain::cal_gdmepsl(kv.nkstot, + this->ld.deepks_param, kv.kvec_d, this->ld.phialpha, - this->ld.inl_index, this->ld.dm_r, ucell, ORB, @@ -228,16 +216,9 @@ template void test_deepks::check_gvepsl(torch::Tensor& gdmepsl) { std::vector gevdm; - DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl2l, this->ld.pdm, gevdm); + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.deepks_param, this->ld.pdm, gevdm); torch::Tensor gvepsl; - DeePKS_domain::cal_gvepsl(ucell.nat, - this->ld.inlmax, - this->ld.des_per_atom, - this->ld.inl2l, - gevdm, - gdmepsl, - gvepsl, - 0); + DeePKS_domain::cal_gvepsl(ucell.nat, this->ld.deepks_param, gevdm, gdmepsl, gvepsl, 0); DeePKS_domain::check_tensor(gvepsl, "gvepsl.dat", 0); // 0 for rank this->compare_with_ref("gvepsl.dat", "gvepsl_ref.dat"); } @@ -248,17 +229,14 @@ void test_deepks::check_orbpre() using TH = std::conditional_t::value, ModuleBase::matrix, ModuleBase::ComplexMatrix>; std::vector gevdm; torch::Tensor orbpre; - DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl2l, this->ld.pdm, gevdm); + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.deepks_param, this->ld.pdm, gevdm); DeePKS_domain::cal_orbital_precalc(dm, - this->ld.lmaxd, - this->ld.inlmax, ucell.nat, kv.nkstot, - this->ld.inl2l, + this->ld.deepks_param, kv.kvec_d, this->ld.phialpha, gevdm, - this->ld.inl_index, ucell, ORB, ParaO, @@ -273,17 +251,14 @@ void test_deepks::check_vdpre() { std::vector gevdm; torch::Tensor vdpre; - DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl2l, this->ld.pdm, gevdm); + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.deepks_param, this->ld.pdm, gevdm); DeePKS_domain::cal_v_delta_precalc(PARAM.sys.nlocal, - this->ld.lmaxd, - this->ld.inlmax, ucell.nat, kv.nkstot, - this->ld.inl2l, + this->ld.deepks_param, kv.kvec_d, this->ld.phialpha, gevdm, - this->ld.inl_index, ucell, ORB, ParaO, @@ -298,29 +273,41 @@ void test_deepks::check_vdrpre() { std::vector gevdm; torch::Tensor vdrpre; - DeePKS_domain::cal_gevdm(ucell.nat, this->ld.inlmax, this->ld.inl2l, this->ld.pdm, gevdm); + torch::Tensor overlap_out; + torch::Tensor iRmat; + DeePKS_domain::cal_gevdm(ucell.nat, this->ld.deepks_param, this->ld.pdm, gevdm); // normally use hR to get R_size, here use 3 instead for Bravo lattice R in [-1,0,1] int R_size = 3; DeePKS_domain::cal_vdr_precalc(PARAM.sys.nlocal, - this->ld.lmaxd, - this->ld.inlmax, ucell.nat, kv.nkstot, R_size, - this->ld.inl2l, + this->ld.deepks_param, kv.kvec_d, this->ld.phialpha, gevdm, - this->ld.inl_index, ucell, ORB, ParaO, Test_Deepks::GridD, vdrpre); + DeePKS_domain::prepare_phialpha_iRmat(PARAM.sys.nlocal, + R_size, + this->ld.deepks_param, + this->ld.phialpha, + ucell, + ORB, + Test_Deepks::GridD, + overlap_out, + iRmat); // vdrpre is large, we only check the main element in Bravo lattice vector (0, 0, 0) and (1, 0, 0) torch::Tensor vdrpre_sliced = vdrpre.slice(0, 0, 2, 1).slice(1, 0, 1, 1).slice(2, 0, 1, 1); DeePKS_domain::check_tensor(vdrpre_sliced, "vdr_precalc.dat", 0); // 0 for rank + DeePKS_domain::check_tensor(overlap_out, "phialpha_r.dat", 0); // 0 for rank + DeePKS_domain::check_tensor(iRmat, "iRmat.dat", 0); // 0 for rank this->compare_with_ref("vdr_precalc.dat", "vdrpre_ref.dat"); + this->compare_with_ref("phialpha_r.dat", "phialpha_r_ref.dat"); + this->compare_with_ref("iRmat.dat", "iRmat_ref.dat"); } template @@ -331,11 +318,7 @@ void test_deepks::check_edelta(std::vector& descriptor) if (PARAM.inp.deepks_equiv) { DeePKS_domain::cal_edelta_gedm_equiv(ucell.nat, - this->ld.lmaxd, - this->ld.nmaxd, - this->ld.inlmax, - this->ld.des_per_atom, - this->ld.inl2l, + this->ld.deepks_param, descriptor, this->ld.gedm, this->ld.E_delta, @@ -344,9 +327,7 @@ void test_deepks::check_edelta(std::vector& descriptor) else { DeePKS_domain::cal_edelta_gedm(ucell.nat, - this->ld.inlmax, - this->ld.des_per_atom, - this->ld.inl2l, + this->ld.deepks_param, descriptor, this->ld.pdm, this->ld.model_deepks, @@ -359,7 +340,7 @@ void test_deepks::check_edelta(std::vector& descriptor) ofs.close(); this->compare_with_ref("E_delta.dat", "E_delta_ref.dat"); - // DeePKS_domain::check_gedm(this->ld.inlmax, this->ld.inl2l, this->ld.gedm); + // DeePKS_domain::check_gedm(this->ld.deepks_param, this->ld.gedm); // this->compare_with_ref("gedm.dat", "gedm_ref.dat"); } @@ -412,10 +393,10 @@ void test_deepks::check_f_delta_and_stress_delta() Test_Deepks::GridD, ParaO, nks, + this->ld.deepks_param, kv.kvec_d, this->ld.phialpha, this->ld.gedm, - this->ld.inl_index, fvnl_dalpha, cal_stress, svnl_dalpha); diff --git a/source/source_lcao/module_deepks/test/klist.h b/source/source_lcao/module_deepks/test/klist.h index 12e0e4d40e..7400d82c45 100644 --- a/source/source_lcao/module_deepks/test/klist.h +++ b/source/source_lcao/module_deepks/test/klist.h @@ -1,72 +1,69 @@ -///klist : adapted from klist from source_pw/module_pwdft -///deals with k point sampling +/// klist : adapted from klist from source_pw/module_pwdft +/// deals with k point sampling -#include "source_base/vector3.h" +#include "source_base/global_function.h" #include "source_base/matrix3.h" #include "source_base/memory.h" -#include "source_base/global_function.h" -#include +#include "source_base/vector3.h" + #include +#include namespace Test_Deepks { class K_Vectors { -public: - - ModuleBase::Vector3 *kvec_c; // Cartesian coordinates of k points - std::vector> kvec_d; // Direct coordinates of k points + public: + ModuleBase::Vector3* kvec_c; // Cartesian coordinates of k points + std::vector> kvec_d; // Direct coordinates of k points - double *wk; // wk, weight of k points + double* wk; // wk, weight of k points - int *isk; // distinguish spin up and down k points + int* isk; // distinguish spin up and down k points - int nkstot; // total number of k points + int nkstot; // total number of k points - int nmp[3]; // Number of Monhorst-Pack + int nmp[3]; // Number of Monhorst-Pack K_Vectors(); ~K_Vectors(); - void set( - const std::string &k_file_name, - const int& nspin, - const ModuleBase::Matrix3 &reciprocal_vec, - const ModuleBase::Matrix3 &latvec, - bool &GAMMA_ONLY_LOCAL, - std::ofstream &ofs_running, - std::ofstream &ofs_warning); + void set(const std::string& k_file_name, + const int& nspin, + const ModuleBase::Matrix3& reciprocal_vec, + const ModuleBase::Matrix3& latvec, + bool& GAMMA_ONLY_LOCAL, + std::ofstream& ofs_running, + std::ofstream& ofs_warning); -private: + private: int nspin; bool kc_done; bool kd_done; - double koffset[3]; // used only in automatic k-points. - std::string k_kword; //LiuXh add 20180619 - int k_nkstot; //LiuXh add 20180619 + double koffset[3]; // used only in automatic k-points. + std::string k_kword; // LiuXh add 20180619 + int k_nkstot; // LiuXh add 20180619 // step 1 : generate kpoints - bool read_kpoints( - const std::string &fn, - bool &GAMMA_ONLY_LOCAL, - std::ofstream &ofs_warning, - std::ofstream &ofs_running); - void Monkhorst_Pack(const int *nmp_in,const double *koffset_in,const int tipo); - double Monkhorst_Pack_formula( const int &k_type, const double &offset, - const int& n, const int &dim); + bool read_kpoints(const std::string& fn, + bool& GAMMA_ONLY_LOCAL, + std::ofstream& ofs_warning, + std::ofstream& ofs_running); + void Monkhorst_Pack(const int* nmp_in, const double* koffset_in, const int tipo); + double Monkhorst_Pack_formula(const int& k_type, const double& offset, const int& n, const int& dim); // step 2 : set both kvec and kved; normalize weight - void set_both_kvec(const ModuleBase::Matrix3 &G,const ModuleBase::Matrix3 &Rm, std::ofstream &ofs_running); - void renew(const int &kpoint_number); - void normalize_wk( const int °spin ); + void set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& Rm, std::ofstream& ofs_running); + void renew(const int& kpoint_number); + void normalize_wk(const int& degspin); // step 3 : *2 or *4 kpoints. // *2 for LSDA // *4 for non-collinear - void set_kup_and_kdw(std::ofstream &ofs_running); + void set_kup_and_kdw(std::ofstream& ofs_running); // step 4 // print k lists. - void print_klists(std::ofstream &fn_running); + void print_klists(std::ofstream& fn_running); }; -} +} // namespace Test_Deepks diff --git a/source/source_lcao/module_deepks/test/klist_1.cpp b/source/source_lcao/module_deepks/test/klist_1.cpp index b4e9fdea81..709a193ece 100644 --- a/source/source_lcao/module_deepks/test/klist_1.cpp +++ b/source/source_lcao/module_deepks/test/klist_1.cpp @@ -1,605 +1,589 @@ #include "klist.h" - #include "source_io/module_parameter/parameter.h" namespace Test_Deepks { - K_Vectors::K_Vectors() - { - nspin = 0; // default spin. - kc_done = false; - kd_done = false; - - kvec_c = new ModuleBase::Vector3[1]; - kvec_d.resize(1); - - wk = nullptr; - isk = nullptr; - - nkstot = 0; - } - - K_Vectors::~K_Vectors() - { - delete[] kvec_c; - kvec_d.clear(); - delete[] wk; - delete[] isk; - } - - void K_Vectors::set( - const std::string &k_file_name, - const int& nspin_in, - const ModuleBase::Matrix3 &reciprocal_vec, - const ModuleBase::Matrix3 &latvec, - bool &GAMMA_ONLY_LOCAL, - std::ofstream &ofs_running, - std::ofstream &ofs_warning) - { - - ofs_running << "\n\n\n\n"; - ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - ofs_running << " | |" << std::endl; - ofs_running << " | Setup K-points |" << std::endl; - ofs_running << " | We setup the k-points according to input parameters. |" << std::endl; - ofs_running << " | The reduced k-points are set according to symmetry operations. |" << std::endl; - ofs_running << " | We treat the spin as another set of k-points. |" << std::endl; - ofs_running << " | |" << std::endl; - ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; - ofs_running << "\n\n\n\n"; - - ofs_running << "\n SETUP K-POINTS" << std::endl; - - // (1) set nspin, read kpoints. - this->nspin = nspin_in; - ModuleBase::GlobalFunc::OUT(ofs_running,"nspin",nspin); - - bool read_succesfully = this->read_kpoints( - k_file_name, - GAMMA_ONLY_LOCAL, - ofs_warning, - ofs_running); - if(!read_succesfully) - { - ofs_warning << "in K_Vectors::set, something wrong while reading KPOINTS." << std::endl; - exit(1); - } - - // (2) - this->set_both_kvec(reciprocal_vec, latvec, ofs_running); - - int deg = 0; - if(PARAM.inp.nspin == 1) - { - deg = 2; - } - else if(PARAM.inp.nspin == 2||PARAM.inp.nspin==4) - { - deg = 1; - } - else - { - ofs_warning << "In K_Vectors::set, Only available for nspin = 1 or 2 or 4" << std::endl; - exit(1); - } - this->normalize_wk(deg); - - // It's very important in parallel case, - // firstly do the mpi_k() and then - // do set_kup_and_kdw() - - this->set_kup_and_kdw(ofs_running); - - this->print_klists(ofs_running); - //std::cout << " NUMBER OF K-POINTS : " << nkstot << std::endl; - - return; - } - - void K_Vectors::renew(const int &kpoint_number) - { - delete[] kvec_c; - delete[] wk; - delete[] isk; - - kvec_c = new ModuleBase::Vector3[kpoint_number]; - kvec_d.resize(kpoint_number); - wk = new double[kpoint_number]; - isk = new int[kpoint_number]; - - ModuleBase::Memory::record("KV::kvec_c",sizeof(double) * kpoint_number*3); - ModuleBase::Memory::record("KV::kvec_d",sizeof(double) * kpoint_number*3); - ModuleBase::Memory::record("KV::wk",sizeof(double) * kpoint_number*3); - ModuleBase::Memory::record("KV::isk",sizeof(int) * kpoint_number*3); - - return; - } - - bool K_Vectors::read_kpoints(const std::string &fn, bool &GAMMA_ONLY_LOCAL, std::ofstream &ofs_warning, std::ofstream &ofs_running) - { - - std::ifstream ifk(fn.c_str()); - ifk >> std::setiosflags(std::ios::uppercase); - - ifk.clear(); - ifk.seekg(0); - - std::string word; - std::string kword; - - int ierr = 0; - - ifk.rdstate(); - - while (ifk.good()) - { - ifk >> word; - ifk.ignore(150, '\n'); //LiuXh add 20180416, fix bug in k-point file when the first line with comments - if (word == "K_POINTS" || word == "KPOINTS" || word == "K" ) - { - ierr = 1; - break; - } - - ifk.rdstate(); - } - - if (ierr == 0) - { - ofs_warning << " symbol K_POINTS not found." << std::endl; - return 0; - } - - //input k-points are in 2pi/a units - ModuleBase::GlobalFunc::READ_VALUE(ifk, nkstot); - - //std::cout << " nkstot = " << nkstot << std::endl; - ModuleBase::GlobalFunc::READ_VALUE(ifk, kword); - - // mohan update 2021-02-22 - int max_kpoints = 100000; - if (nkstot > 100000) - { - ofs_warning << " nkstot > MAX_KPOINTS" << std::endl; - return 0; - } - - int k_type = 0; - if (nkstot == 0) // nkstot==0, use monkhorst_pack. add by dwan - { - if (kword == "Gamma") - { - k_type = 0; - ModuleBase::GlobalFunc::OUT(ofs_running,"Input type of k points","Monkhorst-Pack(Gamma)"); - } - else if (kword == "Monkhorst-Pack" || kword == "MP" || kword == "mp") - { - k_type = 1; - ModuleBase::GlobalFunc::OUT(ofs_running,"Input type of k points","Monkhorst-Pack"); - } - else - { - ofs_warning << " Error: neither Gamma nor Monkhorst-Pack." << std::endl; - return 0; - } - - ifk >> nmp[0] >> nmp[1] >> nmp[2]; - - ifk >> koffset[0] >> koffset[1] >> koffset[2]; - this->Monkhorst_Pack(nmp, koffset, k_type); - } - else if (nkstot > 0) - { - if (kword == "Cartesian" || kword == "C") - { - this->renew(nkstot * nspin);//mohan fix bug 2009-09-01 - for (int i = 0;i < nkstot;i++) - { - ifk >> kvec_c[i].x >> kvec_c[i].y >> kvec_c[i].z; - ModuleBase::GlobalFunc::READ_VALUE(ifk, wk[i]); - } - - this->kc_done = true; - } - else if (kword == "Direct" || kword == "D") - { - this->renew(nkstot * nspin);//mohan fix bug 2009-09-01 - for (int i = 0;i < nkstot;i++) - { - ifk >> kvec_d[i].x >> kvec_d[i].y >> kvec_d[i].z; - ModuleBase::GlobalFunc::READ_VALUE(ifk, wk[i]); - } - this->kd_done = true; - } - else if (kword == "Line_Cartesian" ) - { - //std::cout << " kword = " << kword << std::endl; - - // how many special points. - int nks_special = this->nkstot; - //std::cout << " nks_special = " << nks_special << std::endl; - - //------------------------------------------ - // number of points to the next k points - //------------------------------------------ - int* nkl = new int[nks_special]; - - //------------------------------------------ - // cartesian coordinates of special points. - //------------------------------------------ - double *ksx = new double[nks_special]; - double *ksy = new double[nks_special]; - double *ksz = new double[nks_special]; - std::vector kposx; - std::vector kposy; - std::vector kposz; - ModuleBase::GlobalFunc::ZEROS(nkl, nks_special); - - //recalculate nkstot. - nkstot = 0; - for(int iks=0; iks> ksx[iks]; - ifk >> ksy[iks]; - ifk >> ksz[iks]; - ModuleBase::GlobalFunc::READ_VALUE( ifk, nkl[iks] ); - //std::cout << " nkl[" << iks << "]=" << nkl[iks] << std::endl; - assert(nkl[iks] >= 0); - nkstot += nkl[iks]; - } - assert( nkl[nks_special-1] == 1); - - //std::cout << " nkstot = " << nkstot << std::endl; - this->renew(nkstot * nspin);//mohan fix bug 2009-09-01 - - int count = 0; - for(int iks=1; ikskc_done = true; - - } - - else if (kword == "Line_Direct" || kword == "L" || kword == "Line" ) - { - //std::cout << " kword = " << kword << std::endl; - - // how many special points. - int nks_special = this->nkstot; - //std::cout << " nks_special = " << nks_special << std::endl; - - //------------------------------------------ - // number of points to the next k points - //------------------------------------------ - int* nkl = new int[nks_special]; - - //------------------------------------------ - // cartesian coordinates of special points. - //------------------------------------------ - double *ksx = new double[nks_special]; - double *ksy = new double[nks_special]; - double *ksz = new double[nks_special]; - std::vector kposx; - std::vector kposy; - std::vector kposz; - ModuleBase::GlobalFunc::ZEROS(nkl, nks_special); - - //recalculate nkstot. - nkstot = 0; - for(int iks=0; iks> ksx[iks]; - ifk >> ksy[iks]; - ifk >> ksz[iks]; - ModuleBase::GlobalFunc::READ_VALUE( ifk, nkl[iks] ); - //std::cout << " nkl[" << iks << "]=" << nkl[iks] << std::endl; - assert(nkl[iks] >= 0); - nkstot += nkl[iks]; - } - assert( nkl[nks_special-1] == 1); - - //std::cout << " nkstot = " << nkstot << std::endl; - this->renew(nkstot * nspin);//mohan fix bug 2009-09-01 - - int count = 0; - for(int iks=1; ikskd_done = true; - - } - - else - { - ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; - return 0; - } - } - - ModuleBase::GlobalFunc::OUT(ofs_running,"nkstot",nkstot); - return 1; - } // END SUBROUTINE - - - double K_Vectors::Monkhorst_Pack_formula( const int &k_type, const double &offset, - const int& n, const int &dim) - { - double coordinate; - if (k_type==1) coordinate = (offset + 2.0 * (double)n - (double)dim - 1.0) / (2.0 * (double)dim); - else coordinate = (offset + (double)n - 1.0) / (double)dim; - - return coordinate; - } - - //add by dwan - void K_Vectors::Monkhorst_Pack(const int *nmp_in, const double *koffset_in, const int k_type) - { - const int mpnx = nmp_in[0]; - const int mpny = nmp_in[1]; - const int mpnz = nmp_in[2]; - - this->nkstot = mpnx * mpny * mpnz; - // only can renew after nkstot is estimated. - this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 - for (int x = 1;x <= mpnx;x++) - { - double v1 = Monkhorst_Pack_formula( k_type, koffset_in[0], x, mpnx); - if( std::abs(v1) < 1.0e-10 ) v1 = 0.0; //mohan update 2012-06-10 - for (int y = 1;y <= mpny;y++) - { - double v2 = Monkhorst_Pack_formula( k_type, koffset_in[1], y, mpny); - if( std::abs(v2) < 1.0e-10 ) v2 = 0.0; - for (int z = 1;z <= mpnz;z++) - { - double v3 = Monkhorst_Pack_formula( k_type, koffset_in[2], z, mpnz); - if( std::abs(v3) < 1.0e-10 ) v3 = 0.0; - // index of nks kpoint - const int i = mpnx * mpny * (z - 1) + mpnx * (y - 1) + (x - 1); - kvec_d[i].set(v1, v2, v3); - } - } - } - - const double weight = 1.0 / static_cast(nkstot); - for (int ik=0; ikkd_done = true; - - return; - } - - void K_Vectors::set_both_kvec(const ModuleBase::Matrix3 &G, const ModuleBase::Matrix3 &R, std::ofstream &ofs_running) - { - // set cartesian k vectors. - if (!kc_done && kd_done) - { - for (int i = 0;i < nkstot;i++) - { - //wrong!! kvec_c[i] = G * kvec_d[i]; - // mohan fixed bug 2010-1-10 - if( std::abs(kvec_d[i].x) < 1.0e-10 ) kvec_d[i].x = 0.0; - if( std::abs(kvec_d[i].y) < 1.0e-10 ) kvec_d[i].y = 0.0; - if( std::abs(kvec_d[i].z) < 1.0e-10 ) kvec_d[i].z = 0.0; - - // mohan add2012-06-10 - if( std::abs(kvec_c[i].x) < 1.0e-10 ) kvec_c[i].x = 0.0; - if( std::abs(kvec_c[i].y) < 1.0e-10 ) kvec_c[i].y = 0.0; - if( std::abs(kvec_c[i].z) < 1.0e-10 ) kvec_c[i].z = 0.0; - } - kc_done = true; - } - - // set direct k vectors - else if (kc_done && !kd_done) - { - ModuleBase::Matrix3 RT = R.Transpose(); - for (int i = 0;i < nkstot;i++) - { - // std::cout << " ik=" << i - // << " kvec.x=" << kvec_c[i].x - // << " kvec.y=" << kvec_c[i].y - // << " kvec.z=" << kvec_c[i].z << std::endl; - //wrong! kvec_d[i] = RT * kvec_c[i]; - // mohan fixed bug 2011-03-07 - kvec_d[i] = kvec_c[i] * RT; - } - kd_done = true; - } - - ofs_running << "\n " << std::setw(8) << "KPOINTS" - << std::setw(20) << "DIRECT_X" - << std::setw(20) << "DIRECT_Y" - << std::setw(20) << "DIRECT_Z" - << std::setw(20) << "WEIGHT" << std::endl; - - for(int i=0; ikvec_d[i].x - << std::setw(20) << this->kvec_d[i].y - << std::setw(20) << this->kvec_d[i].z - << std::setw(20) << this->wk[i] << std::endl; - } - - return; - } - - - void K_Vectors::normalize_wk(const int °spin) - { - double sum = 0.0; - - for (int ik = 0;ik < nkstot;ik++) - { - sum += this->wk[ik]; - } - assert(sum>0.0); - - for (int ik = 0;ik < nkstot;ik++) - { - this->wk[ik] /= sum; - } - - for (int ik = 0;ik < nkstot;ik++) - { - this->wk[ik] *= degspin; - } - - return; - } - - //---------------------------------------------------------- - // This routine sets the k vectors for the up and down spin - //---------------------------------------------------------- - // from set_kup_and_kdw.f90 - void K_Vectors::set_kup_and_kdw(std::ofstream &ofs_running) - { - //========================================================================= - // on output: the number of points is doubled and xk and wk in the - // first (nks/2) positions correspond to up spin - // those in the second (nks/2) ones correspond to down spin - //========================================================================= - switch (nspin) - { - case 1: - - for (int ik = 0; ik < nkstot; ik++) - { - this->isk[ik] = 0; - } - - break; - - case 2: - - for (int ik = 0; ik < nkstot; ik++) - { - this->kvec_c[ik+nkstot] = kvec_c[ik]; - this->kvec_d[ik+nkstot] = kvec_d[ik]; - this->wk[ik+nkstot] = wk[ik]; - this->isk[ik] = 0; - this->isk[ik+nkstot] = 1; - } - - this->nkstot *= 2; - - ModuleBase::GlobalFunc::OUT(ofs_running,"nkstot(nspin=2)",nkstot); - break; - case 4: - - for (int ik = 0; ik < nkstot; ik++) - { - this->isk[ik] = 0; - } - - break; - } - - return; - } // end subroutine set_kup_and_kdw - - - void K_Vectors::print_klists(std::ofstream &ofs_running) - { - ofs_running << "\n " << std::setw(8) << "KPOINTS" - << std::setw(20) << "CARTESIAN_X" - << std::setw(20) << "CARTESIAN_Y" - << std::setw(20) << "CARTESIAN_Z" - << std::setw(20) << "WEIGHT" << std::endl; - for(int i=0; ikvec_c[i].x - << std::setw(20) << this->kvec_c[i].y - << std::setw(20) << this->kvec_c[i].z - << std::setw(20) << this->wk[i] << std::endl; - } - - ofs_running << "\n " << std::setw(8) << "KPOINTS" - << std::setw(20) << "DIRECT_X" - << std::setw(20) << "DIRECT_Y" - << std::setw(20) << "DIRECT_Z" - << std::setw(20) << "WEIGHT" << std::endl; - for(int i=0; ikvec_d[i].x - << std::setw(20) << this->kvec_d[i].y - << std::setw(20) << this->kvec_d[i].z - << std::setw(20) << this->wk[i] << std::endl; - } - - return; - } +K_Vectors::K_Vectors() +{ + nspin = 0; // default spin. + kc_done = false; + kd_done = false; + + kvec_c = new ModuleBase::Vector3[1]; + kvec_d.resize(1); + + wk = nullptr; + isk = nullptr; + + nkstot = 0; +} + +K_Vectors::~K_Vectors() +{ + delete[] kvec_c; + kvec_d.clear(); + delete[] wk; + delete[] isk; +} + +void K_Vectors::set(const std::string& k_file_name, + const int& nspin_in, + const ModuleBase::Matrix3& reciprocal_vec, + const ModuleBase::Matrix3& latvec, + bool& GAMMA_ONLY_LOCAL, + std::ofstream& ofs_running, + std::ofstream& ofs_warning) +{ + ofs_running << "\n\n\n\n"; + ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + ofs_running << " | |" << std::endl; + ofs_running << " | Setup K-points |" << std::endl; + ofs_running << " | We setup the k-points according to input parameters. |" << std::endl; + ofs_running << " | The reduced k-points are set according to symmetry operations. |" << std::endl; + ofs_running << " | We treat the spin as another set of k-points. |" << std::endl; + ofs_running << " | |" << std::endl; + ofs_running << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + ofs_running << "\n\n\n\n"; + + ofs_running << "\n SETUP K-POINTS" << std::endl; + + // (1) set nspin, read kpoints. + this->nspin = nspin_in; + ModuleBase::GlobalFunc::OUT(ofs_running, "nspin", nspin); + + bool read_succesfully = this->read_kpoints(k_file_name, GAMMA_ONLY_LOCAL, ofs_warning, ofs_running); + if (!read_succesfully) + { + ofs_warning << "in K_Vectors::set, something wrong while reading KPOINTS." << std::endl; + exit(1); + } + + // (2) + this->set_both_kvec(reciprocal_vec, latvec, ofs_running); + + int deg = 0; + if (PARAM.inp.nspin == 1) + { + deg = 2; + } + else if (PARAM.inp.nspin == 2 || PARAM.inp.nspin == 4) + { + deg = 1; + } + else + { + ofs_warning << "In K_Vectors::set, Only available for nspin = 1 or 2 or 4" << std::endl; + exit(1); + } + this->normalize_wk(deg); + + // It's very important in parallel case, + // firstly do the mpi_k() and then + // do set_kup_and_kdw() + + this->set_kup_and_kdw(ofs_running); + + this->print_klists(ofs_running); + // std::cout << " NUMBER OF K-POINTS : " << nkstot << std::endl; + + return; } + +void K_Vectors::renew(const int& kpoint_number) +{ + delete[] kvec_c; + delete[] wk; + delete[] isk; + + kvec_c = new ModuleBase::Vector3[kpoint_number]; + kvec_d.resize(kpoint_number); + wk = new double[kpoint_number]; + isk = new int[kpoint_number]; + + ModuleBase::Memory::record("KV::kvec_c", sizeof(double) * kpoint_number * 3); + ModuleBase::Memory::record("KV::kvec_d", sizeof(double) * kpoint_number * 3); + ModuleBase::Memory::record("KV::wk", sizeof(double) * kpoint_number * 3); + ModuleBase::Memory::record("KV::isk", sizeof(int) * kpoint_number * 3); + + return; +} + +bool K_Vectors::read_kpoints(const std::string& fn, + bool& GAMMA_ONLY_LOCAL, + std::ofstream& ofs_warning, + std::ofstream& ofs_running) +{ + + std::ifstream ifk(fn.c_str()); + ifk >> std::setiosflags(std::ios::uppercase); + + ifk.clear(); + ifk.seekg(0); + + std::string word; + std::string kword; + + int ierr = 0; + + ifk.rdstate(); + + while (ifk.good()) + { + ifk >> word; + ifk.ignore(150, '\n'); // LiuXh add 20180416, fix bug in k-point file when the first line with comments + if (word == "K_POINTS" || word == "KPOINTS" || word == "K") + { + ierr = 1; + break; + } + + ifk.rdstate(); + } + + if (ierr == 0) + { + ofs_warning << " symbol K_POINTS not found." << std::endl; + return 0; + } + + // input k-points are in 2pi/a units + ModuleBase::GlobalFunc::READ_VALUE(ifk, nkstot); + + // std::cout << " nkstot = " << nkstot << std::endl; + ModuleBase::GlobalFunc::READ_VALUE(ifk, kword); + + // mohan update 2021-02-22 + int max_kpoints = 100000; + if (nkstot > 100000) + { + ofs_warning << " nkstot > MAX_KPOINTS" << std::endl; + return 0; + } + + int k_type = 0; + if (nkstot == 0) // nkstot==0, use monkhorst_pack. add by dwan + { + if (kword == "Gamma") + { + k_type = 0; + ModuleBase::GlobalFunc::OUT(ofs_running, "Input type of k points", "Monkhorst-Pack(Gamma)"); + } + else if (kword == "Monkhorst-Pack" || kword == "MP" || kword == "mp") + { + k_type = 1; + ModuleBase::GlobalFunc::OUT(ofs_running, "Input type of k points", "Monkhorst-Pack"); + } + else + { + ofs_warning << " Error: neither Gamma nor Monkhorst-Pack." << std::endl; + return 0; + } + + ifk >> nmp[0] >> nmp[1] >> nmp[2]; + + ifk >> koffset[0] >> koffset[1] >> koffset[2]; + this->Monkhorst_Pack(nmp, koffset, k_type); + } + else if (nkstot > 0) + { + if (kword == "Cartesian" || kword == "C") + { + this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 + for (int i = 0; i < nkstot; i++) + { + ifk >> kvec_c[i].x >> kvec_c[i].y >> kvec_c[i].z; + ModuleBase::GlobalFunc::READ_VALUE(ifk, wk[i]); + } + + this->kc_done = true; + } + else if (kword == "Direct" || kword == "D") + { + this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 + for (int i = 0; i < nkstot; i++) + { + ifk >> kvec_d[i].x >> kvec_d[i].y >> kvec_d[i].z; + ModuleBase::GlobalFunc::READ_VALUE(ifk, wk[i]); + } + this->kd_done = true; + } + else if (kword == "Line_Cartesian") + { + // std::cout << " kword = " << kword << std::endl; + + // how many special points. + int nks_special = this->nkstot; + // std::cout << " nks_special = " << nks_special << std::endl; + + //------------------------------------------ + // number of points to the next k points + //------------------------------------------ + int* nkl = new int[nks_special]; + + //------------------------------------------ + // cartesian coordinates of special points. + //------------------------------------------ + double* ksx = new double[nks_special]; + double* ksy = new double[nks_special]; + double* ksz = new double[nks_special]; + std::vector kposx; + std::vector kposy; + std::vector kposz; + ModuleBase::GlobalFunc::ZEROS(nkl, nks_special); + + // recalculate nkstot. + nkstot = 0; + for (int iks = 0; iks < nks_special; iks++) + { + ifk >> ksx[iks]; + ifk >> ksy[iks]; + ifk >> ksz[iks]; + ModuleBase::GlobalFunc::READ_VALUE(ifk, nkl[iks]); + // std::cout << " nkl[" << iks << "]=" << nkl[iks] << std::endl; + assert(nkl[iks] >= 0); + nkstot += nkl[iks]; + } + assert(nkl[nks_special - 1] == 1); + + // std::cout << " nkstot = " << nkstot << std::endl; + this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 + + int count = 0; + for (int iks = 1; iks < nks_special; iks++) + { + double dx = (ksx[iks] - ksx[iks - 1]) / nkl[iks - 1]; + double dy = (ksy[iks] - ksy[iks - 1]) / nkl[iks - 1]; + double dz = (ksz[iks] - ksz[iks - 1]) / nkl[iks - 1]; + // GlobalV::ofs_running << " dx=" << dx << " dy=" << dy << " dz=" << dz << std::endl; + for (int is = 0; is < nkl[iks - 1]; is++) + { + kvec_c[count].x = ksx[iks - 1] + is * dx; + kvec_c[count].y = ksy[iks - 1] + is * dy; + kvec_c[count].z = ksz[iks - 1] + is * dz; + ++count; + } + } + + // deal with the last special k point. + kvec_c[count].x = ksx[nks_special - 1]; + kvec_c[count].y = ksy[nks_special - 1]; + kvec_c[count].z = ksz[nks_special - 1]; + ++count; + + // std::cout << " count = " << count << std::endl; + assert(count == nkstot); + + for (int ik = 0; ik < nkstot; ik++) + { + wk[ik] = 1.0; + } + + ofs_warning << " Error : nkstot == -1, not implemented yet." << std::endl; + + delete[] nkl; + delete[] ksx; + delete[] ksy; + delete[] ksz; + + this->kc_done = true; + } + + else if (kword == "Line_Direct" || kword == "L" || kword == "Line") + { + // std::cout << " kword = " << kword << std::endl; + + // how many special points. + int nks_special = this->nkstot; + // std::cout << " nks_special = " << nks_special << std::endl; + + //------------------------------------------ + // number of points to the next k points + //------------------------------------------ + int* nkl = new int[nks_special]; + + //------------------------------------------ + // cartesian coordinates of special points. + //------------------------------------------ + double* ksx = new double[nks_special]; + double* ksy = new double[nks_special]; + double* ksz = new double[nks_special]; + std::vector kposx; + std::vector kposy; + std::vector kposz; + ModuleBase::GlobalFunc::ZEROS(nkl, nks_special); + + // recalculate nkstot. + nkstot = 0; + for (int iks = 0; iks < nks_special; iks++) + { + ifk >> ksx[iks]; + ifk >> ksy[iks]; + ifk >> ksz[iks]; + ModuleBase::GlobalFunc::READ_VALUE(ifk, nkl[iks]); + // std::cout << " nkl[" << iks << "]=" << nkl[iks] << std::endl; + assert(nkl[iks] >= 0); + nkstot += nkl[iks]; + } + assert(nkl[nks_special - 1] == 1); + + // std::cout << " nkstot = " << nkstot << std::endl; + this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 + + int count = 0; + for (int iks = 1; iks < nks_special; iks++) + { + double dx = (ksx[iks] - ksx[iks - 1]) / nkl[iks - 1]; + double dy = (ksy[iks] - ksy[iks - 1]) / nkl[iks - 1]; + double dz = (ksz[iks] - ksz[iks - 1]) / nkl[iks - 1]; + // GlobalV::ofs_running << " dx=" << dx << " dy=" << dy << " dz=" << dz << std::endl; + for (int is = 0; is < nkl[iks - 1]; is++) + { + kvec_d[count].x = ksx[iks - 1] + is * dx; + kvec_d[count].y = ksy[iks - 1] + is * dy; + kvec_d[count].z = ksz[iks - 1] + is * dz; + ++count; + } + } + + // deal with the last special k point. + kvec_d[count].x = ksx[nks_special - 1]; + kvec_d[count].y = ksy[nks_special - 1]; + kvec_d[count].z = ksz[nks_special - 1]; + ++count; + + // std::cout << " count = " << count << std::endl; + assert(count == nkstot); + + for (int ik = 0; ik < nkstot; ik++) + { + wk[ik] = 1.0; + } + + ofs_warning << " Error : nkstot == -1, not implemented yet." << std::endl; + + delete[] nkl; + delete[] ksx; + delete[] ksy; + delete[] ksz; + + this->kd_done = true; + } + + else + { + ofs_warning << " Error : neither Cartesian nor Direct kpoint." << std::endl; + return 0; + } + } + + ModuleBase::GlobalFunc::OUT(ofs_running, "nkstot", nkstot); + return 1; +} // END SUBROUTINE + +double K_Vectors::Monkhorst_Pack_formula(const int& k_type, const double& offset, const int& n, const int& dim) +{ + double coordinate; + if (k_type == 1) + coordinate = (offset + 2.0 * (double)n - (double)dim - 1.0) / (2.0 * (double)dim); + else + coordinate = (offset + (double)n - 1.0) / (double)dim; + + return coordinate; +} + +// add by dwan +void K_Vectors::Monkhorst_Pack(const int* nmp_in, const double* koffset_in, const int k_type) +{ + const int mpnx = nmp_in[0]; + const int mpny = nmp_in[1]; + const int mpnz = nmp_in[2]; + + this->nkstot = mpnx * mpny * mpnz; + // only can renew after nkstot is estimated. + this->renew(nkstot * nspin); // mohan fix bug 2009-09-01 + for (int x = 1; x <= mpnx; x++) + { + double v1 = Monkhorst_Pack_formula(k_type, koffset_in[0], x, mpnx); + if (std::abs(v1) < 1.0e-10) + v1 = 0.0; // mohan update 2012-06-10 + for (int y = 1; y <= mpny; y++) + { + double v2 = Monkhorst_Pack_formula(k_type, koffset_in[1], y, mpny); + if (std::abs(v2) < 1.0e-10) + v2 = 0.0; + for (int z = 1; z <= mpnz; z++) + { + double v3 = Monkhorst_Pack_formula(k_type, koffset_in[2], z, mpnz); + if (std::abs(v3) < 1.0e-10) + v3 = 0.0; + // index of nks kpoint + const int i = mpnx * mpny * (z - 1) + mpnx * (y - 1) + (x - 1); + kvec_d[i].set(v1, v2, v3); + } + } + } + + const double weight = 1.0 / static_cast(nkstot); + for (int ik = 0; ik < nkstot; ik++) + { + wk[ik] = weight; + } + this->kd_done = true; + + return; +} + +void K_Vectors::set_both_kvec(const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::ofstream& ofs_running) +{ + // set cartesian k vectors. + if (!kc_done && kd_done) + { + for (int i = 0; i < nkstot; i++) + { + // wrong!! kvec_c[i] = G * kvec_d[i]; + // mohan fixed bug 2010-1-10 + if (std::abs(kvec_d[i].x) < 1.0e-10) + kvec_d[i].x = 0.0; + if (std::abs(kvec_d[i].y) < 1.0e-10) + kvec_d[i].y = 0.0; + if (std::abs(kvec_d[i].z) < 1.0e-10) + kvec_d[i].z = 0.0; + + // mohan add2012-06-10 + if (std::abs(kvec_c[i].x) < 1.0e-10) + kvec_c[i].x = 0.0; + if (std::abs(kvec_c[i].y) < 1.0e-10) + kvec_c[i].y = 0.0; + if (std::abs(kvec_c[i].z) < 1.0e-10) + kvec_c[i].z = 0.0; + } + kc_done = true; + } + + // set direct k vectors + else if (kc_done && !kd_done) + { + ModuleBase::Matrix3 RT = R.Transpose(); + for (int i = 0; i < nkstot; i++) + { + // std::cout << " ik=" << i + // << " kvec.x=" << kvec_c[i].x + // << " kvec.y=" << kvec_c[i].y + // << " kvec.z=" << kvec_c[i].z << std::endl; + // wrong! kvec_d[i] = RT * kvec_c[i]; + // mohan fixed bug 2011-03-07 + kvec_d[i] = kvec_c[i] * RT; + } + kd_done = true; + } + + ofs_running << "\n " << std::setw(8) << "KPOINTS" << std::setw(20) << "DIRECT_X" << std::setw(20) << "DIRECT_Y" + << std::setw(20) << "DIRECT_Z" << std::setw(20) << "WEIGHT" << std::endl; + + for (int i = 0; i < nkstot; i++) + { + ofs_running << " " << std::setw(8) << i + 1 << std::setw(20) << this->kvec_d[i].x << std::setw(20) + << this->kvec_d[i].y << std::setw(20) << this->kvec_d[i].z << std::setw(20) << this->wk[i] + << std::endl; + } + + return; +} + +void K_Vectors::normalize_wk(const int& degspin) +{ + double sum = 0.0; + + for (int ik = 0; ik < nkstot; ik++) + { + sum += this->wk[ik]; + } + assert(sum > 0.0); + + for (int ik = 0; ik < nkstot; ik++) + { + this->wk[ik] /= sum; + } + + for (int ik = 0; ik < nkstot; ik++) + { + this->wk[ik] *= degspin; + } + + return; +} + +//---------------------------------------------------------- +// This routine sets the k vectors for the up and down spin +//---------------------------------------------------------- +// from set_kup_and_kdw.f90 +void K_Vectors::set_kup_and_kdw(std::ofstream& ofs_running) +{ + //========================================================================= + // on output: the number of points is doubled and xk and wk in the + // first (nks/2) positions correspond to up spin + // those in the second (nks/2) ones correspond to down spin + //========================================================================= + switch (nspin) + { + case 1: + + for (int ik = 0; ik < nkstot; ik++) + { + this->isk[ik] = 0; + } + + break; + + case 2: + + for (int ik = 0; ik < nkstot; ik++) + { + this->kvec_c[ik + nkstot] = kvec_c[ik]; + this->kvec_d[ik + nkstot] = kvec_d[ik]; + this->wk[ik + nkstot] = wk[ik]; + this->isk[ik] = 0; + this->isk[ik + nkstot] = 1; + } + + this->nkstot *= 2; + + ModuleBase::GlobalFunc::OUT(ofs_running, "nkstot(nspin=2)", nkstot); + break; + case 4: + + for (int ik = 0; ik < nkstot; ik++) + { + this->isk[ik] = 0; + } + + break; + } + + return; +} // end subroutine set_kup_and_kdw + +void K_Vectors::print_klists(std::ofstream& ofs_running) +{ + ofs_running << "\n " << std::setw(8) << "KPOINTS" << std::setw(20) << "CARTESIAN_X" << std::setw(20) + << "CARTESIAN_Y" << std::setw(20) << "CARTESIAN_Z" << std::setw(20) << "WEIGHT" << std::endl; + for (int i = 0; i < nkstot; i++) + { + ofs_running << " " << std::setw(8) << i + 1 << std::setw(20) << this->kvec_c[i].x << std::setw(20) + << this->kvec_c[i].y << std::setw(20) << this->kvec_c[i].z << std::setw(20) << this->wk[i] + << std::endl; + } + + ofs_running << "\n " << std::setw(8) << "KPOINTS" << std::setw(20) << "DIRECT_X" << std::setw(20) << "DIRECT_Y" + << std::setw(20) << "DIRECT_Z" << std::setw(20) << "WEIGHT" << std::endl; + for (int i = 0; i < nkstot; i++) + { + ofs_running << " " << std::setw(8) << i + 1 << std::setw(20) << this->kvec_d[i].x << std::setw(20) + << this->kvec_d[i].y << std::setw(20) << this->kvec_d[i].z << std::setw(20) << this->wk[i] + << std::endl; + } + + return; +} + +} // namespace Test_Deepks diff --git a/source/source_lcao/module_deepks/test/parallel_orbitals.h b/source/source_lcao/module_deepks/test/parallel_orbitals.h index b9371d9a23..5445ff87ab 100644 --- a/source/source_lcao/module_deepks/test/parallel_orbitals.h +++ b/source/source_lcao/module_deepks/test/parallel_orbitals.h @@ -1,5 +1,5 @@ -///adapted from parallel_orbitals from source_basis/module_ao -///deals with the parallelization of atomic basis +/// adapted from parallel_orbitals from source_basis/module_ao +/// deals with the parallelization of atomic basis #include "source_base/global_function.h" #include "source_base/global_variable.h" @@ -7,19 +7,18 @@ namespace Test_Deepks { - class Parallel_Orbitals - { - public: - - Parallel_Orbitals(); - ~Parallel_Orbitals(); +class Parallel_Orbitals +{ + public: + Parallel_Orbitals(); + ~Parallel_Orbitals(); - int* global2local_row; - int* global2local_col; - void set_global2local(void); + int* global2local_row; + int* global2local_col; + void set_global2local(void); - int ncol; - int nrow; - int nloc; - }; -} + int ncol; + int nrow; + int nloc; +}; +} // namespace Test_Deepks diff --git a/source/source_lcao/module_operator_lcao/deepks_lcao.cpp b/source/source_lcao/module_operator_lcao/deepks_lcao.cpp index 5e75cc9b40..b1dc5f9f12 100644 --- a/source/source_lcao/module_operator_lcao/deepks_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/deepks_lcao.cpp @@ -152,13 +152,8 @@ void hamilt::DeePKS>::contributeHR() { ModuleBase::timer::tick("DeePKS", "contributeHR"); - 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->inl2l, - this->ld->inl_index, + this->ld->deepks_param, this->kvec_d, this->ld->dm_r, this->ld->phialpha, @@ -170,19 +165,13 @@ void hamilt::DeePKS>::contributeHR() std::vector descriptor; DeePKS_domain::cal_descriptor(this->ucell->nat, - inlmax, - this->ld->inl2l, + this->ld->deepks_param, this->ld->pdm, - descriptor, - this->ld->des_per_atom); + descriptor); if (PARAM.inp.deepks_equiv) { DeePKS_domain::cal_edelta_gedm_equiv(this->ucell->nat, - this->ld->lmaxd, - this->ld->nmaxd, - inlmax, - this->ld->des_per_atom, - this->ld->inl2l, + this->ld->deepks_param, descriptor, this->ld->gedm, this->ld->E_delta, @@ -191,9 +180,7 @@ void hamilt::DeePKS>::contributeHR() else { DeePKS_domain::cal_edelta_gedm(this->ucell->nat, - inlmax, - this->ld->des_per_atom, - this->ld->inl2l, + this->ld->deepks_param, descriptor, this->ld->pdm, this->ld->model_deepks, @@ -253,7 +240,7 @@ void hamilt::DeePKS>::calculate_HR() { for (int N0 = 0; N0 < ptr_orb_->Alpha[0].getNchi(L0); ++N0) { - const int inl = this->ld->inl_index[T0](I0, L0, N0); + const int inl = this->ld->deepks_param.inl_index[T0](I0, L0, N0); const double* pgedm = this->ld->gedm[inl]; const int nm = 2 * L0 + 1; @@ -274,7 +261,7 @@ void hamilt::DeePKS>::calculate_HR() { const double* pgedm = this->ld->gedm[iat0]; int nproj = 0; - for (int il = 0; il < this->ld->lmaxd + 1; il++) + for (int il = 0; il < this->ld->deepks_param.lmaxd + 1; il++) { nproj += (2 * il + 1) * ptr_orb_->Alpha[0].getNchi(il); } diff --git a/source/source_lcao/setup_deepks.cpp b/source/source_lcao/setup_deepks.cpp index baa776f0a2..190571661c 100644 --- a/source/source_lcao/setup_deepks.cpp +++ b/source/source_lcao/setup_deepks.cpp @@ -1,51 +1,51 @@ #include "source_lcao/setup_deepks.h" -#include "source_lcao/LCAO_domain.h" + #include "source_io/module_parameter/parameter.h" // use parameter +#include "source_lcao/LCAO_domain.h" template -Setup_DeePKS::Setup_DeePKS(){} +Setup_DeePKS::Setup_DeePKS() +{ +} template -Setup_DeePKS::~Setup_DeePKS(){} - +Setup_DeePKS::~Setup_DeePKS() +{ +} template -void Setup_DeePKS::build_overlap( - const UnitCell &ucell, - const LCAO_Orbitals &orb, - const Parallel_Orbitals &pv, - const Grid_Driver &gd, - TwoCenterIntegrator &overlap_orb_alpha, - const Input_para &inp) +void Setup_DeePKS::build_overlap(const UnitCell& ucell, + const LCAO_Orbitals& orb, + const Parallel_Orbitals& pv, + const Grid_Driver& gd, + TwoCenterIntegrator& overlap_orb_alpha, + const Input_para& inp) { #ifdef __MLALGO - // 9) for each ionic step, the overlap must be rebuilt - // since it depends on ionic positions - if (PARAM.globalv.deepks_setorb) - { - // allocate , phialpha is different every ion step, so it is allocated here - DeePKS_domain::allocate_phialpha(inp.cal_force, ucell, orb, gd, &pv, this->ld.phialpha); + // 9) for each ionic step, the overlap must be rebuilt + // since it depends on ionic positions + if (PARAM.globalv.deepks_setorb) + { + // allocate , phialpha is different every ion step, so it is allocated here + DeePKS_domain::allocate_phialpha(inp.cal_force, ucell, orb, gd, &pv, this->ld.phialpha); - // build and save at beginning - DeePKS_domain::build_phialpha(inp.cal_force, ucell, orb, gd, - &pv, overlap_orb_alpha, this->ld.phialpha); + // build and save at beginning + DeePKS_domain::build_phialpha(inp.cal_force, ucell, orb, gd, &pv, overlap_orb_alpha, this->ld.phialpha); - if (inp.deepks_out_unittest) - { - DeePKS_domain::check_phialpha(inp.cal_force, ucell, orb, - gd, &pv, this->ld.phialpha, GlobalV::MY_RANK); - } - } + if (inp.deepks_out_unittest) + { + DeePKS_domain::check_phialpha(inp.cal_force, ucell, orb, gd, &pv, this->ld.phialpha, GlobalV::MY_RANK); + } + } #endif } - template -void Setup_DeePKS::before_runner(const UnitCell& ucell, // unitcell - const int nks, // number of k points - const LCAO_Orbitals &orb, // orbital info - Parallel_Orbitals &pv, // parallel orbitals - const Input_para &inp) +void Setup_DeePKS::before_runner(const UnitCell& ucell, // unitcell + const int nks, // number of k points + const LCAO_Orbitals& orb, // orbital info + Parallel_Orbitals& pv, // parallel orbitals + const Input_para& inp) { #ifdef __MLALGO LCAO_domain::DeePKS_init(ucell, pv, nks, orb, this->ld, GlobalV::ofs_running); @@ -54,23 +54,26 @@ void Setup_DeePKS::before_runner(const UnitCell& ucell, // unitcell // load the DeePKS model from deep neural network DeePKS_domain::load_model(inp.deepks_model, this->ld.model_deepks); // read pdm from file for NSCF or SCF-restart, do it only once in whole calculation - DeePKS_domain::read_pdm((inp.init_chg == "file"), inp.deepks_equiv, - this->ld.init_pdm, ucell.nat, orb.Alpha[0].getTotal_nchi() * ucell.nat, - this->ld.lmaxd, this->ld.inl2l, *orb.Alpha, this->ld.pdm); + DeePKS_domain::read_pdm((inp.init_chg == "file"), + inp.deepks_equiv, + this->ld.init_pdm, + ucell.nat, + this->ld.deepks_param, + *orb.Alpha, + this->ld.pdm); } #endif } template -void Setup_DeePKS::delta_e( - const UnitCell& ucell, - const K_Vectors &kv, - const LCAO_Orbitals& orb, - const Parallel_Orbitals &pv, // parallel orbitals - const Grid_Driver &gd, - const std::vector>& dm_vec, - elecstate::fenergy &f_en, - const Input_para &inp) +void Setup_DeePKS::delta_e(const UnitCell& ucell, + const K_Vectors& kv, + const LCAO_Orbitals& orb, + const Parallel_Orbitals& pv, // parallel orbitals + const Grid_Driver& gd, + const std::vector>& dm_vec, + elecstate::fenergy& f_en, + const Input_para& inp) { #ifdef __MLALGO if (inp.deepks_scf) @@ -83,115 +86,104 @@ void Setup_DeePKS::delta_e( #endif } - - template -void Setup_DeePKS::write_forces( - const ModuleBase::matrix &fcs, - const ModuleBase::matrix &fvnl_dalpha, - const Input_para &inp) +void Setup_DeePKS::write_forces(const ModuleBase::matrix& fcs, + const ModuleBase::matrix& fvnl_dalpha, + const Input_para& inp) { #ifdef __MLALGO - // DeePKS force - if (inp.deepks_out_labels) // not parallelized yet - { - if (inp.deepks_out_base == "none" - || (inp.deepks_out_base != "none" && this->dpks_out_type == "tot") ) - { - const std::string file_ftot = PARAM.globalv.global_out_dir - + (inp.deepks_out_labels == 1 ? "deepks_ftot.npy" : "deepks_force.npy"); - LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Hartree/Bohr, F_tot - - if (inp.deepks_out_labels == 1) - { - // this base only considers subtracting the deepks_scf part - const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; - if (inp.deepks_scf) - { - LCAO_deepks_io::save_matrix2npy(file_fbase, - fcs - fvnl_dalpha, - GlobalV::MY_RANK); // Hartree/Bohr, F_base - } - else - { - LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot - } - } - } - if (inp.deepks_out_base != "none") - { - // output fcs as tot or base in another dir - // this base considers changing xc functional to base functional - const std::string file_f = PARAM.globalv.global_deepks_label_elec_dir + - (dpks_out_type == "tot" ? "ftot.npy" : "fbase.npy"); - LCAO_deepks_io::save_matrix2npy(file_f, fcs, GlobalV::MY_RANK); - } - } + // DeePKS force + if (inp.deepks_out_labels) // not parallelized yet + { + if (inp.deepks_out_base == "none" || (inp.deepks_out_base != "none" && this->dpks_out_type == "tot")) + { + const std::string file_ftot + = PARAM.globalv.global_out_dir + (inp.deepks_out_labels == 1 ? "deepks_ftot.npy" : "deepks_force.npy"); + LCAO_deepks_io::save_matrix2npy(file_ftot, fcs, GlobalV::MY_RANK); // Hartree/Bohr, F_tot + + if (inp.deepks_out_labels == 1) + { + // this base only considers subtracting the deepks_scf part + const std::string file_fbase = PARAM.globalv.global_out_dir + "deepks_fbase.npy"; + if (inp.deepks_scf) + { + LCAO_deepks_io::save_matrix2npy(file_fbase, + fcs - fvnl_dalpha, + GlobalV::MY_RANK); // Hartree/Bohr, F_base + } + else + { + LCAO_deepks_io::save_matrix2npy(file_fbase, fcs, GlobalV::MY_RANK); // no scf, F_base=F_tot + } + } + } + if (inp.deepks_out_base != "none") + { + // output fcs as tot or base in another dir + // this base considers changing xc functional to base functional + const std::string file_f + = PARAM.globalv.global_deepks_label_elec_dir + (dpks_out_type == "tot" ? "ftot.npy" : "fbase.npy"); + LCAO_deepks_io::save_matrix2npy(file_f, fcs, GlobalV::MY_RANK); + } + } #endif - } - template -void Setup_DeePKS::write_stress( - const ModuleBase::matrix &scs, - const ModuleBase::matrix &svnl_dalpha, - const double &omega, - const Input_para &inp) +void Setup_DeePKS::write_stress(const ModuleBase::matrix& scs, + const ModuleBase::matrix& svnl_dalpha, + const double& omega, + const Input_para& inp) { #ifdef __MLALGO - if (inp.deepks_out_labels == 1) - { - assert(omega>0.0); - - if (inp.deepks_out_base == "none" - || (inp.deepks_out_base != "none" && this->dpks_out_type == "tot") ) - { - const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; - LCAO_deepks_io::save_matrix2npy(file_stot, - scs, - GlobalV::MY_RANK, - omega, - 'U'); // change to energy unit Ry when printing, S_tot; - - // this base only considers subtracting the deepks_scf part - const std::string file_sbase = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; - if (inp.deepks_scf) - { - LCAO_deepks_io::save_matrix2npy(file_sbase, - scs - svnl_dalpha, - GlobalV::MY_RANK, - omega, - 'U'); // change to energy unit Ry when printing, S_base; - } - else - { - LCAO_deepks_io::save_matrix2npy(file_sbase, - scs, - GlobalV::MY_RANK, - omega, - 'U'); // sbase = stot - } - } - if (inp.deepks_out_base != "none") - { - // output scs as tot or base in another dir - // this base considers changing xc functional to base functional - const std::string file_s = PARAM.globalv.global_deepks_label_elec_dir + - (this->dpks_out_type == "tot" ? "stot.npy" : "sbase.npy"); - LCAO_deepks_io::save_matrix2npy(file_s, - scs, - GlobalV::MY_RANK, - omega, - 'U'); // change to energy unit Ry when printing, S_tot; - } - } - else if (inp.deepks_out_labels == 2) - { - const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stress.npy"; - LCAO_deepks_io::save_matrix2npy(file_stot, scs, GlobalV::MY_RANK, omega, - 'F'); // flat mode - } + if (inp.deepks_out_labels == 1) + { + assert(omega > 0.0); + + if (inp.deepks_out_base == "none" || (inp.deepks_out_base != "none" && this->dpks_out_type == "tot")) + { + const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stot.npy"; + LCAO_deepks_io::save_matrix2npy(file_stot, + scs, + GlobalV::MY_RANK, + omega, + 'U'); // change to energy unit Ry when printing, S_tot; + + // this base only considers subtracting the deepks_scf part + const std::string file_sbase = PARAM.globalv.global_out_dir + "deepks_sbase.npy"; + if (inp.deepks_scf) + { + LCAO_deepks_io::save_matrix2npy(file_sbase, + scs - svnl_dalpha, + GlobalV::MY_RANK, + omega, + 'U'); // change to energy unit Ry when printing, S_base; + } + else + { + LCAO_deepks_io::save_matrix2npy(file_sbase, scs, GlobalV::MY_RANK, omega, + 'U'); // sbase = stot + } + } + if (inp.deepks_out_base != "none") + { + // output scs as tot or base in another dir + // this base considers changing xc functional to base functional + const std::string file_s = PARAM.globalv.global_deepks_label_elec_dir + + (this->dpks_out_type == "tot" ? "stot.npy" : "sbase.npy"); + LCAO_deepks_io::save_matrix2npy(file_s, + scs, + GlobalV::MY_RANK, + omega, + 'U'); // change to energy unit Ry when printing, S_tot; + } + } + else if (inp.deepks_out_labels == 2) + { + const std::string file_stot = PARAM.globalv.global_out_dir + "deepks_stress.npy"; + LCAO_deepks_io::save_matrix2npy(file_stot, scs, GlobalV::MY_RANK, omega, + 'F'); // flat mode + } #endif } diff --git a/tests/09_DeePKS/100_NO_GO_deepks_UT/iRmat_ref.dat b/tests/09_DeePKS/100_NO_GO_deepks_UT/iRmat_ref.dat new file mode 100644 index 0000000000..af12828313 --- /dev/null +++ b/tests/09_DeePKS/100_NO_GO_deepks_UT/iRmat_ref.dat @@ -0,0 +1,125 @@ +0 0 0 +0 1 0 +0 1 1 +1 1 0 +0 1 1 +0 2 0 +0 0 0 +0 0 1 +1 0 0 +0 0 1 +0 2 2 +0 0 2 +0 0 0 +1 0 2 +0 0 0 +2 2 0 +2 0 0 +2 0 1 +0 0 0 +2 0 1 +0 2 2 +0 0 2 +0 0 0 +1 0 2 +0 0 0 +0 0 0 +0 1 0 +0 1 1 +1 1 0 +0 1 1 +0 2 0 +0 0 0 +0 0 1 +1 0 0 +0 0 1 +0 2 2 +0 0 2 +0 0 0 +1 0 2 +0 0 0 +2 2 0 +2 0 0 +2 0 1 +0 0 0 +2 0 1 +0 2 2 +0 0 2 +0 0 0 +1 0 2 +0 0 0 +0 0 0 +0 1 0 +0 1 1 +0 1 1 +1 1 0 +0 2 0 +0 0 0 +0 0 1 +0 0 1 +1 0 0 +0 2 2 +0 0 2 +0 0 0 +0 0 0 +1 0 2 +0 2 2 +0 0 2 +0 0 0 +0 0 0 +1 0 2 +2 2 0 +2 0 0 +2 0 1 +2 0 1 +0 0 0 +0 0 0 +0 0 1 +0 0 1 +1 0 0 +0 2 0 +0 0 2 +0 0 0 +0 0 0 +1 0 2 +0 2 2 +0 0 2 +0 0 0 +0 0 0 +1 0 2 +0 2 2 +2 0 0 +2 0 1 +2 0 1 +0 0 0 +2 2 0 +0 1 0 +0 1 1 +0 1 1 +1 1 0 +0 0 0 +0 0 0 +0 1 1 +0 1 1 +1 1 0 +0 1 0 +0 2 2 +0 0 0 +0 0 0 +1 0 2 +0 0 2 +0 2 2 +0 0 0 +0 0 0 +1 0 2 +0 0 2 +2 2 0 +2 0 1 +2 0 1 +0 0 0 +2 0 0 +0 2 0 +0 0 1 +0 0 1 +1 0 0 +0 0 0 diff --git a/tests/09_DeePKS/100_NO_GO_deepks_UT/phialpha_r_ref.dat b/tests/09_DeePKS/100_NO_GO_deepks_UT/phialpha_r_ref.dat new file mode 100644 index 0000000000..20267a8218 --- /dev/null +++ b/tests/09_DeePKS/100_NO_GO_deepks_UT/phialpha_r_ref.dat @@ -0,0 +1,200 @@ +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.3164811566 -0.1534527887 -0.08570920339 0.1156008495 0.2071111386 0.04376523785 -0.05902865122 -0.1057560668 -0.04843479778 -0.03999389619 -0.07165329159 -0.05960192152 0.09664284637 0.02577082509 0.02127965328 0.03812474768 0.0317125448 -0.05142100316 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.3322570668 -0.1586938262 -0.0903368452 0.120593251 -0.2161250816 0.04535098 -0.06054032659 0.1084992976 -0.05005969208 -0.0420453347 0.07535290178 -0.06207471837 -0.1005907542 0.02610664785 0.02192707748 -0.03929731865 0.0323726085 0.05245912008 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.4114648095 -0.1794735013 0.3147366589 0.00616018533 0.0008881864795 -0.1415856307 -0.00277118569 -0.0003995544826 0.1812524527 0.00614576416 0.0008861072096 5.889370801e-05 1.73433392e-05 -0.08212807954 -0.00278473367 -0.0004015078544 -2.66855817e-05 -7.858515125e-06 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.2888131094 -0.1435047394 -0.07183317419 -0.2203678974 6.227069716e-05 0.03767474148 0.1155775679 -3.265945635e-05 -0.04849447325 0.06953694541 -1.964947762e-05 0.1066617855 -6.028014379e-05 0.02659435774 -0.03813404454 1.077576891e-05 -0.05849329812 3.30576167e-05 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.8170531799 -0.247178539 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0.8150193757 0 0 0.02256967339 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0.8150193757 0 0 0.02256967339 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0.8150193757 0 0 0.02256967339 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.1150935953 -0.06410311347 -0.07219303423 0.02672869898 0.04955194699 0.04282914818 -0.01585703416 -0.0293971254 0.02212911648 -0.02039197394 -0.03780438443 -0.009199155668 0.0139966691 -0.01369969442 0.01262426414 0.02340394 0.005695013702 -0.008665058534 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.1192926978 -0.06631180479 -0.07548975067 0.02728433166 -0.05095474957 0.04467427619 -0.0161466657 0.03015464398 0.02373719934 -0.02102182856 0.03925923579 -0.009450806053 -0.01418950255 -0.01465640929 0.01297981784 -0.02424040932 0.005835350652 0.008761233963 +0.4396940853 -0.1615590582 -0.2875027763 -0.005627149984 -0.0008113324951 0.1066242109 0.002086903071 0.0003008933973 0.101011939 0.003425032573 0.0004938272893 3.282144628e-05 9.665437873e-06 -0.03734667878 -0.001266321512 -0.0001825804883 -1.213492211e-05 -3.57355781e-06 +0.5744883648 -0.1122607916 -0.264174326 -0.01109747707 -0.001600053986 -0.001574581052 0.002595663622 0.0003742473986 0.0933320631 0.008411901177 0.001212842876 0.0001308736388 3.854038041e-05 0.04085026979 -0.0009523410587 -0.000137310228 -3.151678617e-05 -9.281234472e-06 +0.0112441773 -0.002197225083 -0.01109747707 0.3026016138 -3.131706718e-05 0.002595663622 -0.1341416212 7.324959673e-06 0.01091173032 -0.267721887 3.854038041e-05 -0.005242658682 -0.0007555108486 -0.003247473646 0.1193145981 -9.281234472e-06 0.002335928811 0.0003367056548 +0.001621205486 -0.0003167998213 -0.001600053986 -3.131706718e-05 0.3028143039 0.0003742473986 7.324959673e-06 -0.1341913688 0.001573272688 3.854038041e-05 -0.2679836342 0.0007566345057 -0.005245111447 -0.0004682265274 -9.281234472e-06 0.1193776317 -0.0003369762521 0.002336519482 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.1059639845 -0.05925714388 -0.06358704877 -0.05471114836 -0.0001147382641 0.03792215417 0.03262872933 6.842780449e-05 0.01683755767 0.03983957875 8.355014007e-05 0.01713917818 7.188765948e-05 -0.01048499156 -0.02480868396 -5.202788496e-05 -0.01067281503 -4.476548901e-05 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.8778350314 -0.01899237841 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.083970075 -0.04734479653 -0.0005828940149 0.05700333953 0.03215589146 0.000351764923 -0.03440037954 -0.01940543975 -0.01774617636 -0.0004769440117 -0.000269046691 0.01589994474 0.02631106082 0.01121348344 0.0003013721755 0.0001700056706 -0.01004688353 -0.01662547687 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.08767473994 -0.04937413925 -0.0003705339936 0.05996259465 -0.03321317421 0.0002231864639 -0.03611771037 0.02000553535 -0.01867033849 -0.0003058479934 0.000169408658 0.0171547553 -0.02741498179 0.01176822392 0.0001927810615 -0.0001067810861 -0.01081292671 0.01728011759 +0.3230572918 -0.1446202823 0.07747822275 0.237685627 -6.716427333e-05 -0.03466476827 -0.1063436523 3.005017266e-05 -0.04188394327 0.0600580083 -1.69709567e-05 0.09212217131 -5.206304871e-05 0.01791384367 -0.02568692648 7.258511056e-06 -0.03940083111 2.226746678e-05 +-0.146857873 0.05155932492 0.1452936884 -0.144083087 4.071443426e-05 -0.07923392798 0.05857215097 -1.655108896e-05 0.102011228 0.1100781004 -3.110543838e-05 -0.0865347701 4.890531656e-05 -0.0508813716 -0.06407838456 1.810701888e-05 0.03822967435 -2.160558495e-05 +-0.4505266692 0.1581723229 -0.144083087 -0.2497538998 0.0001249026564 0.05857215097 0.08135943147 -5.077498963e-05 0.02457218512 -0.1187977918 4.890531656e-05 -0.09897559591 0.0001029833556 -0.001936895948 0.04744750208 -2.160558495e-05 0.02827818753 -4.113126096e-05 +0.0001273080612 -4.469571537e-05 4.071443426e-05 0.0001249026564 0.1922603193 -1.655108896e-05 -5.077498963e-05 -0.09832663799 -6.943511807e-06 4.890531656e-05 0.05427174835 0.0001220622787 0.1664934233 5.473204692e-07 -2.160558495e-05 -0.02901184662 -5.829029298e-05 -0.0890017699 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.1059639845 -0.05925714388 0.06358704877 0.05471114836 0.0001147382641 -0.03792215417 -0.03262872933 -6.842780449e-05 0.01683755767 0.03983957875 8.355014007e-05 0.01713917818 7.188765948e-05 -0.01048499156 -0.02480868396 -5.202788496e-05 -0.01067281503 -4.476548901e-05 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.8778350314 -0.01899237841 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.09289685966 -0.05221929161 0.0002685632109 -0.000570033376 -0.07291767165 -0.0001613215545 0.0003424097814 0.04380045987 -0.01997563887 -1.99231389e-06 -0.0002548533054 -0.03459553134 0.0005409336953 0.01254696941 1.251399347e-06 0.000160076814 0.02172992194 -0.0003397677827 +0.3508189026 -0.1500864749 0.08914085147 -0.1202293073 -0.2154035099 -0.03795971177 0.05119840989 0.09172736198 -0.03814920318 -0.03150080814 -0.05643702679 -0.04694488092 0.07611980955 0.01534048149 0.01266704214 0.02269434464 0.01887738188 -0.03060914599 +-0.1706457715 0.05466627924 0.1565431659 0.08031206029 0.1438875434 -0.08537396789 -0.03001504274 -0.05377512106 0.1188507876 -0.05554711773 -0.09951853164 0.05139868212 -0.08334152344 -0.05683263389 0.03267010273 0.05853194161 -0.02080802258 0.03373962581 +0.2301596021 -0.0737315022 0.08031206029 0.1077670091 -0.1940692666 -0.03001504274 -0.06714484484 0.07252954677 -0.004352982816 0.02023740889 -0.08334152344 -0.1593608109 -0.04890248224 -0.004116710019 -0.01574247999 0.03373962581 0.07469759414 0.03804075672 +0.4123552504 -0.1320977781 0.1438875434 -0.1940692666 -0.131607174 -0.05377512106 0.07252954677 0.02231651763 -0.007798828739 -0.08334152344 -0.08255995535 0.03710808178 0.1113532806 -0.007375521054 0.03373962581 0.02587356517 -0.0332659621 -0.0348971405 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.1150935953 -0.06410311347 0.07219303423 -0.02672869898 -0.04955194699 -0.04282914818 0.01585703416 0.0293971254 0.02212911648 -0.02039197394 -0.03780438443 -0.009199155668 0.0139966691 -0.01369969442 0.01262426414 0.02340394 0.005695013702 -0.008665058534 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.083970075 -0.04734479653 0.0005828940149 -0.05700333953 -0.03215589146 -0.000351764923 0.03440037954 0.01940543975 -0.01774617636 -0.0004769440117 -0.000269046691 0.01589994474 0.02631106082 0.01121348344 0.0003013721755 0.0001700056706 -0.01004688353 -0.01662547687 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.8778350314 -0.01899237841 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.09289685966 -0.05221929161 -0.0002685632109 0.000570033376 0.07291767165 0.0001613215545 -0.0003424097814 -0.04380045987 -0.01997563887 -1.99231389e-06 -0.0002548533054 -0.03459553134 0.0005409336953 0.01254696941 1.251399347e-06 0.000160076814 0.02172992194 -0.0003397677827 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.3662277175 -0.1526733097 0.09197392871 -0.1227786408 0.2200416984 -0.03811860367 0.05088561957 -0.09119630325 -0.03733268828 -0.03135587353 0.0561954394 -0.04629305564 -0.07501690712 0.01455628605 0.01222588261 -0.02191100958 0.0180499983 0.02924963641 +-0.1772065705 0.05345772271 0.1678497891 0.08293590189 -0.1486362497 -0.09072180452 -0.0293105681 0.05252987935 0.1253013898 -0.05959116734 0.1067982313 0.05297278213 0.08584126111 -0.05821576032 0.03477040377 -0.06231489983 -0.02021050449 -0.03275069052 +0.2365581439 -0.07136225042 0.08293590189 0.1192638574 0.1984188016 -0.0293105681 -0.07355093851 -0.07012364566 -0.001907026722 0.02362039255 0.08584126111 -0.1661863379 0.05651026729 -0.006220721217 -0.0180271203 -0.03275069052 0.07543925459 -0.04312872381 +-0.4239553023 0.1278941572 -0.1486362497 0.1984188016 -0.1256252743 0.05252987935 -0.07012364566 0.01299589353 0.003417739406 0.08584126111 -0.0823251717 -0.04436826569 0.1098982377 0.01114866603 -0.03275069052 0.02239388869 0.03849623203 -0.02989424561 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.1192926978 -0.06631180479 0.07548975067 -0.02728433166 0.05095474957 -0.04467427619 0.0161466657 -0.03015464398 0.02373719934 -0.02102182856 0.03925923579 -0.009450806053 -0.01418950255 -0.01465640929 0.01297981784 -0.02424040932 0.005835350652 0.008761233963 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.08767473994 -0.04937413925 0.0003705339936 -0.05996259465 0.03321317421 -0.0002231864639 0.03611771037 -0.02000553535 -0.01867033849 -0.0003058479934 0.000169408658 0.0171547553 -0.02741498179 0.01176822392 0.0001927810615 -0.0001067810861 -0.01081292671 0.01728011759 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.8778350314 -0.01899237841 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/tests/09_DeePKS/100_NO_KP_deepks_UT/iRmat_ref.dat b/tests/09_DeePKS/100_NO_KP_deepks_UT/iRmat_ref.dat new file mode 100644 index 0000000000..3578fe72ea --- /dev/null +++ b/tests/09_DeePKS/100_NO_KP_deepks_UT/iRmat_ref.dat @@ -0,0 +1,75 @@ +0 0 0 +0 0 1 +0 1 0 +0 1 1 +0 1 1 +0 0 2 +0 0 0 +0 1 2 +0 1 0 +0 1 0 +0 2 0 +0 2 1 +0 0 0 +0 0 1 +0 0 1 +0 2 2 +0 2 0 +0 0 2 +0 0 0 +0 0 0 +0 2 2 +0 2 0 +0 0 2 +0 0 0 +0 0 0 +0 0 0 +0 1 0 +0 1 0 +0 1 1 +0 1 0 +0 2 0 +0 0 0 +0 0 0 +0 0 1 +0 0 0 +0 2 0 +0 0 0 +0 0 0 +0 0 1 +0 0 0 +0 2 2 +0 0 2 +0 0 2 +0 0 0 +0 0 2 +0 2 0 +0 0 0 +0 0 0 +0 0 1 +0 0 0 +0 0 0 +0 1 0 +0 1 0 +0 1 1 +0 0 0 +0 2 0 +0 0 0 +0 0 0 +0 0 1 +0 2 0 +0 2 0 +0 0 0 +0 0 0 +0 0 1 +0 2 0 +0 2 2 +0 0 2 +0 0 2 +0 0 0 +0 2 2 +0 0 0 +0 1 0 +0 1 0 +0 1 1 +0 0 0 diff --git a/tests/09_DeePKS/100_NO_KP_deepks_UT/phialpha_r_ref.dat b/tests/09_DeePKS/100_NO_KP_deepks_UT/phialpha_r_ref.dat new file mode 100644 index 0000000000..e737ce22e4 --- /dev/null +++ b/tests/09_DeePKS/100_NO_KP_deepks_UT/phialpha_r_ref.dat @@ -0,0 +1,90 @@ +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.01543554467 -0.008065153447 -0.01244641558 -0.001923079779 0.002999569923 0.007741992235 0.001196205335 -0.001865810032 0.003518338782 0.0009817979232 -0.001531382916 -0.0001086824711 -0.0002366120191 -0.002833911453 -0.0007908074096 0.001233480871 8.754031902e-05 0.0001905835544 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.01242662179 -0.006477161275 0.01097039075 -0.00162689654 0.002537590993 -0.006348711177 0.0009415066869 -0.001468537692 0.003015181109 -0.0008048659007 0.00125540894 -8.551568097e-05 -0.0001861757259 -0.002752907546 0.0007348551651 -0.001146208012 7.807715522e-05 0.0001699813518 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.05764787604 -0.03269664837 -0.02891309865 -0.02565587961 -0.02018788394 0.01762729729 0.01564148563 0.01230784137 0.003618889031 0.01534123614 0.01207158356 0.002592132746 0.01071165351 -0.002325372903 -0.009857747643 -0.007756782004 -0.001665614833 -0.006882937995 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +-3.235665827e-12 6.459383561e-12 -7.496419526e-12 2.186899137e-12 1.720808899e-12 1.285653335e-11 -3.750582742e-12 -2.951227174e-12 -1.081805427e-11 5.870671048e-12 4.619464525e-12 -3.261123371e-13 -1.347617079e-12 1.701898686e-11 -9.235752656e-12 -7.267351792e-12 5.130406489e-13 2.120074165e-12 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.949842409 -0.03558002646 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0.8831112507 0 0 0.1909716474 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0.8831112507 0 0 0.1909716474 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0.8831112507 0 0 0.1909716474 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.04257844364 -0.02417633078 -0.01875373391 0.01113470711 0.02250388227 0.01145805784 -0.006803024862 -0.01374930378 0.0005859403333 -0.00580624991 -0.01173476438 -0.005316990868 0.006967314615 -0.000377863047 0.003744352724 0.007567551803 0.003428837812 -0.004493103787 +-8.823436698e-13 1.761408977e-12 2.044140447e-12 -5.963285493e-13 -4.692340205e-13 -3.505813654e-12 1.022736365e-12 8.047622355e-13 -2.9499256e-12 1.600846362e-12 1.259660594e-12 -8.892607748e-14 -3.674755204e-13 4.640947554e-12 -2.518519113e-12 -1.981751252e-12 1.399022612e-13 5.781280102e-13 +-2.438582561e-12 4.86819372e-12 5.650113754e-12 -1.649916376e-12 -1.298272396e-12 -9.690033745e-12 2.829639951e-12 2.226563415e-12 -8.146925528e-12 4.425258164e-12 3.482110127e-12 -2.4602072e-13 -1.016648824e-12 1.28164617e-11 -6.961689698e-12 -5.477956155e-12 3.870340474e-13 1.599368172e-12 +7.113974986e-13 -1.420177807e-12 -1.649916376e-12 4.757278704e-13 3.7873958e-13 2.829639951e-12 -8.158573101e-13 -6.495460392e-13 2.38291152e-12 -1.279658162e-12 -1.016648824e-12 6.816620704e-14 2.937465207e-13 -3.748754012e-12 2.013049372e-12 1.599368172e-12 -1.072144542e-13 -4.620970402e-13 +5.597785128e-13 -1.117497633e-12 -1.298272396e-12 3.7873958e-13 2.924240944e-13 2.226563415e-12 -6.495460392e-13 -5.014876272e-13 1.875045484e-12 -1.016648824e-12 -7.876167218e-13 5.931043831e-14 2.297681345e-13 -2.949788198e-12 1.599368172e-12 1.238981277e-12 -9.332391695e-14 -3.614428298e-13 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +5.816977142e-06 -1.130009148e-05 -4.443707969e-06 2.638375214e-06 -1.300541872e-05 7.347779566e-06 -4.362617799e-06 2.150477722e-05 -7.675952679e-06 -2.282055052e-06 1.124899951e-05 -1.578378151e-05 -6.678900075e-06 1.153176835e-05 3.428386197e-06 -1.689964255e-05 2.371235462e-05 1.003387224e-05 +0.02620624272 -0.01765994681 0.01837713881 0.0163068534 0.01283140039 -0.01229302457 -0.01090814797 -0.008583312227 0.00296628903 0.01257472668 0.009894695743 0.002124689335 0.008780004033 -0.002019636009 -0.008561664272 -0.006736930763 -0.00144662204 -0.005977978586 +-0.0464567146 0.02679448593 -0.01898647553 -0.02982218359 -0.02346623036 0.01141316902 0.01812643311 0.01426317606 0.005154366219 -0.01533640357 -0.01206778095 -0.003982504916 -0.01645718678 -0.003357712502 0.009722637692 0.007650467816 0.002539024192 0.01049218927 +-0.04122311106 0.02377594023 -0.02982218359 -0.01184067845 -0.02082263089 0.01812643311 0.007069831373 0.01265635109 -0.009687295039 -0.01163579488 -0.01645718678 0.004699734615 -0.008124417219 0.006207837528 0.007356351899 0.01049218927 -0.003051292448 0.00513639787 +-0.0324372968 0.01870861296 -0.02346623036 -0.02082263089 -0.001762864906 0.01426317606 0.01265635109 0.0009443639846 -0.007622657685 -0.01645718678 -0.003670785049 -0.009259466665 -0.003257250993 0.004884771265 0.01049218927 0.002278303515 0.005946602068 0.002021640136 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.8778350314 -0.01899237841 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +5.816977142e-06 -1.130009148e-05 4.443707969e-06 -2.638375214e-06 1.300541872e-05 -7.347779566e-06 4.362617799e-06 -2.150477722e-05 -7.675952679e-06 -2.282055052e-06 1.124899951e-05 -1.578378151e-05 -6.678900075e-06 1.153176835e-05 3.428386197e-06 -1.689964255e-05 2.371235462e-05 1.003387224e-05 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.003090520969 -0.001951458591 -0.003181413023 0.0004717999529 -0.0007359013201 0.002079048377 -0.0003083205228 0.0004809103483 0.001471587277 -0.0003928223135 0.0006127140481 -4.173672609e-05 -9.086480033e-05 -0.001129028787 0.0003013804938 -0.0004700854714 3.202118283e-05 6.971314372e-05 +0.01110414519 -0.006405338519 -0.01079981165 0.001852244286 -0.002889082558 0.006534436311 -0.001121743466 0.001749666341 0.004187020407 -0.001385658729 0.002161314518 -0.0001736254316 -0.0003779989868 -0.003238551237 0.001037660671 -0.00161851618 0.0001273100938 0.0002771661156 +-0.001646732172 0.0009499044576 0.001852244286 0.001415444712 0.0004284476757 -0.001121743466 -0.0008632833679 -0.0002594735394 -0.001051324618 -0.001433248674 -0.0003779989868 0.0002742366783 -0.0003315286583 0.0007583897246 0.0009050529002 0.0002771661156 -0.0001794504284 0.0002093502539 +0.002568530098 -0.001481636316 -0.002889082558 0.0004284476757 0.001021848881 0.001749666341 -0.0002594735394 -0.0006249165636 0.00163982885 -0.0003779989868 -0.00108599781 0.0003474237086 0.000161052247 -0.001182916607 0.0002771661156 0.0006504327127 -0.0002210052305 -9.645843562e-05 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.04257844364 -0.02417633078 0.01875373391 -0.01113470711 -0.02250388227 -0.01145805784 0.006803024862 0.01374930378 0.0005859403333 -0.00580624991 -0.01173476438 -0.005316990868 0.006967314615 -0.000377863047 0.003744352724 0.007567551803 0.003428837812 -0.004493103787 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.00406420697 -0.002591654052 0.00401526231 0.0006203930529 -0.0009676729803 -0.002711808927 -0.0004189981349 0.0006535424149 0.001947926272 0.0005435718634 -0.0008478492831 -6.017198847e-05 -0.0001310001102 -0.00139810647 -0.0003901437904 0.0006085361573 4.318790069e-05 9.402414469e-05 +-0.01419405342 0.008149246796 -0.01291285755 -0.0023437266 0.003655683914 0.008129862255 0.001469242087 -0.00229168567 -0.00511471444 -0.001779628247 0.002775817947 0.0002329346753 0.00050712082 0.003628639963 0.001245910782 -0.001943339298 -0.0001617149694 -0.000352068784 +-0.002193105071 0.001259129718 -0.0023437266 0.001893911254 0.0005648350541 0.001469242087 -0.001152247654 -0.0003540854269 -0.001352524437 0.001775850084 0.00050712082 0.0003606093461 -0.0004279775544 0.0009329841386 -0.001165552804 -0.000352068784 -0.0002399498199 0.0002808967057 +0.003420748363 -0.00196395785 0.003655683914 0.0005648350541 0.001375022403 -0.00229168567 -0.0003540854269 -0.0008269651152 0.00210963251 0.00050712082 0.001309980734 0.0004501951918 0.0002024034506 -0.001455244442 -0.000352068784 -0.000842122865 -0.0002963213072 -0.0001301153287 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0.8778350314 -0.01899237841 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/tests/09_DeePKS/103_NO_GO_deepks_vdelta_r_2/result.ref b/tests/09_DeePKS/103_NO_GO_deepks_vdelta_r_2/result.ref index 3de8b5a4c8..5678831f58 100644 --- a/tests/09_DeePKS/103_NO_GO_deepks_vdelta_r_2/result.ref +++ b/tests/09_DeePKS/103_NO_GO_deepks_vdelta_r_2/result.ref @@ -6,6 +6,6 @@ deepks_e_label 17.126764505645642 deepks_edelta 0.09815855485768665 deepks_hr_label_pass 0 deepks_vdelta_r_pass 0 -deepks_phialpha_r 73.40481582340394 +deepks_phialpha_r 40.10770566607966 deepks_gevdm 54.0 totaltimeref 1.81 diff --git a/tests/09_DeePKS/103_NO_KP_deepks_vdelta_r_2/result.ref b/tests/09_DeePKS/103_NO_KP_deepks_vdelta_r_2/result.ref index 3fce0adcd5..b7e5a997dd 100644 --- a/tests/09_DeePKS/103_NO_KP_deepks_vdelta_r_2/result.ref +++ b/tests/09_DeePKS/103_NO_KP_deepks_vdelta_r_2/result.ref @@ -6,6 +6,6 @@ deepks_e_label 17.12676450564565 deepks_edelta 0.09815855485768665 deepks_hr_label_pass 0 deepks_vdelta_r_pass 0 -deepks_phialpha_r 73.40481582340394 +deepks_phialpha_r 40.10770566607966 deepks_gevdm 54.0 totaltimeref 1.81