From 42a086ef313c83f0033460a62e301b4e9d817e3c Mon Sep 17 00:00:00 2001 From: dyzheng Date: Tue, 2 Apr 2024 10:16:02 +0800 Subject: [PATCH] Test: add unit test of DFTUNew and refactor some code (#3814) * Fix: not use static member in DFTUNew * Test: add unit test of DFTUNew * Test: add nspin=2 test and delete GlobalV::CURRENT_SPIN in DFTUNew * add annotation for dftu.h * Refactor: new dftu code with suggestions by developers --------- Co-authored-by: dyzheng --- source/Makefile.Objects | 2 +- source/module_esolver/esolver_ks_lcao.cpp | 11 +- .../hamilt_lcaodft/CMakeLists.txt | 2 +- .../hamilt_lcaodft/FORCE_STRESS.cpp | 8 +- .../hamilt_lcaodft/hamilt_lcao.cpp | 14 +- .../operator_lcao/CMakeLists.txt | 2 +- .../hamilt_lcaodft/operator_lcao/dftu.hpp | 18 ++ .../operator_lcao/dftu_force_stress.hpp | 50 +-- .../{dftu_new.cpp => dftu_lcao.cpp} | 171 +++++----- .../operator_lcao/{dftu_new.h => dftu_lcao.h} | 88 ++---- .../operator_lcao/test/CMakeLists.txt | 10 + .../operator_lcao/test/test_dftu.cpp | 299 ++++++++++++++++++ .../module_hamilt_lcao/module_dftu/dftu.cpp | 28 ++ source/module_hamilt_lcao/module_dftu/dftu.h | 20 ++ 14 files changed, 546 insertions(+), 177 deletions(-) create mode 100644 source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu.hpp rename source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/{dftu_new.cpp => dftu_lcao.cpp} (78%) rename source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/{dftu_new.h => dftu_lcao.h} (65%) create mode 100644 source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 45903f3a73..82f7492530 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -274,7 +274,7 @@ OBJS_HAMILT_LCAO=hamilt_lcao.o\ deepks_lcao.o\ op_exx_lcao.o\ sc_lambda_lcao.o\ - dftu_new.o\ + dftu_lcao.o\ OBJS_HCONTAINER=base_matrix.o\ atom_pair.o\ diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index d5c7a9ae0e..c00c27837f 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -18,7 +18,6 @@ #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_elecstate/occupy.h" #include "module_hamilt_lcao/module_dftu/dftu.h" -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" #ifdef __EXX @@ -641,15 +640,9 @@ void ESolver_KS_LCAO::iter_init(const int istep, const int iter) if (GlobalV::dft_plus_u) { - if(istep == 0 && iter == 1) + if(istep != 0 || iter != 1) { - hamilt::DFTUNew>::dm_in_dftu = nullptr; - } - else - { - hamilt::DFTUNew>::dm_in_dftu = - dynamic_cast*>(this->pelec) - ->get_DM(); + GlobalC::dftu.set_dmr( dynamic_cast*>(this->pelec)->get_DM() ); } // Calculate U and J if Yukawa potential is used GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt index 3a6bcdd3fd..61a21928f1 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt @@ -14,7 +14,7 @@ if(ENABLE_LCAO) operator_lcao/td_ekinetic_lcao.cpp operator_lcao/td_nonlocal_lcao.cpp operator_lcao/sc_lambda_lcao.cpp - operator_lcao/dftu_new.cpp + operator_lcao/dftu_lcao.cpp FORCE_STRESS.cpp FORCE_gamma.cpp FORCE_gamma_edm.cpp diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index e6944bd2cc..91f0e2665a 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -14,7 +14,7 @@ #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks 2021-06-03 #include "module_elecstate/elecstate_lcao.h" #endif -#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h" template Force_Stress_LCAO::Force_Stress_LCAO(Record_adj& ra, const int nat_in) : RA(&ra), f_pw(nat_in), nat(nat_in) @@ -240,14 +240,14 @@ void Force_Stress_LCAO::getForceStress(const bool isforce, } else { - hamilt::DFTUNew> tmp_dftu(uhm.LM, + hamilt::DFTU> tmp_dftu(uhm.LM, kv.kvec_d, nullptr, nullptr, - &GlobalC::ucell, + GlobalC::ucell, &GlobalC::GridD, &GlobalC::dftu, - uhm.LM->ParaV); + *(uhm.LM->ParaV)); tmp_dftu.cal_force_stress(isforce, isstress, force_dftu, stress_dftu); } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 9701e4c8f0..ea26376182 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -17,7 +17,7 @@ #include "module_hsolver/diago_elpa.h" #endif #include "operator_lcao/op_dftu_lcao.h" -#include "operator_lcao/dftu_new.h" +#include "operator_lcao/dftu_lcao.h" #include "operator_lcao/meta_lcao.h" #include "operator_lcao/op_exx_lcao.h" #include "operator_lcao/td_ekinetic_lcao.h" @@ -226,15 +226,15 @@ HamiltLCAO::HamiltLCAO( } else { - dftu = new DFTUNew>( + dftu = new DFTU>( LM_in, this->kv->kvec_d, this->hR, &(this->getHk(LM_in)), - &GlobalC::ucell, + GlobalC::ucell, &GlobalC::GridD, &GlobalC::dftu, - LM_in->ParaV + *(LM_in->ParaV) ); } this->getOperator()->add(dftu); @@ -397,15 +397,15 @@ HamiltLCAO::HamiltLCAO( } else { - dftu = new DFTUNew>( + dftu = new DFTU>( LM_in, this->kv->kvec_d, this->hR, &(this->getHk(LM_in)), - &GlobalC::ucell, + GlobalC::ucell, &GlobalC::GridD, &GlobalC::dftu, - LM_in->ParaV + *(LM_in->ParaV) ); } this->getOperator()->add(dftu); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/CMakeLists.txt index e1407f6540..58950a5d30 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/CMakeLists.txt @@ -12,7 +12,7 @@ add_library( td_ekinetic_lcao.cpp td_nonlocal_lcao.cpp sc_lambda_lcao.cpp - dftu_new.cpp + dftu_lcao.cpp ) if(ENABLE_COVERAGE) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu.hpp new file mode 100644 index 0000000000..3f8478f7ab --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu.hpp @@ -0,0 +1,18 @@ +namespace hamilt +{ + +#ifndef __DFTUTEMPLATE +#define __DFTUTEMPLATE + +/// The DFTU class template inherits from class T +/// it is used to calculate the non-local pseudopotential of wavefunction basis +/// Template parameters: +/// - T: base class, it would be OperatorLCAO or OperatorPW +template +class DFTU : public T +{ +}; + +#endif + +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp index 1d51b7b11d..cec70d48ce 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_force_stress.hpp @@ -1,24 +1,34 @@ #pragma once -#include "dftu_new.h" +#include "dftu_lcao.h" #include "module_base/parallel_reduce.h" namespace hamilt { template -void DFTUNew>::cal_force_stress( +void DFTU>::cal_force_stress( const bool cal_force, const bool cal_stress, ModuleBase::matrix& force, ModuleBase::matrix& stress) { - ModuleBase::TITLE("DFTUNew", "cal_force_stress"); -#ifdef __DEBUG - assert(this->dm_in_dftu != nullptr); -#endif - ModuleBase::timer::tick("DFTUNew", "cal_force_stress"); + ModuleBase::TITLE("DFTU", "cal_force_stress"); + if(this->dftu->get_dmr(0) == nullptr) + { + ModuleBase::WARNING_QUIT("DFTU", "dmr is not set"); + } + //try to get the density matrix, if the density matrix is empty, skip the calculation and return + const hamilt::HContainer* dmR_tmp[this->nspin]; + dmR_tmp[0] = this->dftu->get_dmr(0); + if(this->nspin==2) dmR_tmp[1] = this->dftu->get_dmr(1); + if(dmR_tmp[0]->size_atom_pairs() == 0) + { + return; + } + // begin the calculation of force and stress + ModuleBase::timer::tick("DFTU", "cal_force_stress"); - const Parallel_Orbitals* paraV = this->dm_in_dftu->get_DMR_pointer(1)->get_atom_pair(0).get_paraV(); + const Parallel_Orbitals* paraV = dmR_tmp[0]->get_atom_pair(0).get_paraV(); const int npol = this->ucell->get_npol(); std::vector stress_tmp(6, 0); if (cal_force) @@ -81,7 +91,7 @@ void DFTUNew>::cal_force_stress( uot.two_center_bundle->overlap_orb_onsite->snap( T1, L1, N1, M1, T0, dtau * this->ucell->lat0, 1 /*cal_deri*/, nlm); #else - ModuleBase::WARNING_QUIT("DFTUNew", "old two center integral method not implemented"); + ModuleBase::WARNING_QUIT("DFTU", "old two center integral method not implemented"); #endif // select the elements of nlm with target_L std::vector nlm_target(tlp1 * 4); @@ -105,17 +115,13 @@ void DFTUNew>::cal_force_stress( } } //first iteration to calculate occupation matrix - std::vector occ(tlp1 * tlp1 * GlobalV::NSPIN, 0); + std::vector occ(tlp1 * tlp1 * this->nspin, 0); for(int i=0;idftu->locale[iat0][target_L][0][is].c[ii]; } - hamilt::HContainer* dmR_tmp[GlobalV::NSPIN]; - dmR_tmp[0] = this->dm_in_dftu->get_DMR_pointer(1); - if(GlobalV::NSPIN==2) dmR_tmp[1] = this->dm_in_dftu->get_DMR_pointer(2); - //calculate VU const double u_value = this->dftu->U[T0]; @@ -151,9 +157,9 @@ void DFTUNew>::cal_force_stress( ModuleBase::Vector3 R_vector(R_index2[0] - R_index1[0], R_index2[1] - R_index1[1], R_index2[2] - R_index1[2]); - const hamilt::BaseMatrix* tmp[GlobalV::NSPIN]; + const hamilt::BaseMatrix* tmp[this->nspin]; tmp[0] = dmR_tmp[0]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); - if(GlobalV::NSPIN == 2) + if(this->nspin == 2) { tmp[1] = dmR_tmp[1]->find_matrix(iat1, iat2, R_vector[0], R_vector[1], R_vector[2]); } @@ -161,10 +167,10 @@ void DFTUNew>::cal_force_stress( if (tmp[0] != nullptr) { // calculate force - if (cal_force) this->cal_force_IJR(iat1, iat2, T0, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, GlobalV::NSPIN, force_tmp1, force_tmp2); + if (cal_force) this->cal_force_IJR(iat1, iat2, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, this->nspin, force_tmp1, force_tmp2); // calculate stress - if (cal_stress) this->cal_stress_IJR(iat1, iat2, T0, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, GlobalV::NSPIN, dis1, dis2, stress_tmp.data()); + if (cal_stress) this->cal_stress_IJR(iat1, iat2, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, this->nspin, dis1, dis2, stress_tmp.data()); } } } @@ -202,14 +208,13 @@ void DFTUNew>::cal_force_stress( stress.c[3] = stress.c[1]; // stress(1,0) } - ModuleBase::timer::tick("DFTUNew", "cal_force_stress"); + ModuleBase::timer::tick("DFTU", "cal_force_stress"); } template -void DFTUNew>::cal_force_IJR( +void DFTU>::cal_force_IJR( const int& iat1, const int& iat2, - const int& T0, const Parallel_Orbitals* paraV, const std::unordered_map>& nlm1_all, const std::unordered_map>& nlm2_all, @@ -272,10 +277,9 @@ void DFTUNew>::cal_force_IJR( } template -void DFTUNew>::cal_stress_IJR( +void DFTU>::cal_stress_IJR( const int& iat1, const int& iat2, - const int& T0, const Parallel_Orbitals* paraV, const std::unordered_map>& nlm1_all, const std::unordered_map>& nlm2_all, diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp similarity index 78% rename from source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp rename to source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp index 0426028b0f..062895d628 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.cpp @@ -1,4 +1,4 @@ -#include "dftu_new.h" +#include "dftu_lcao.h" #include "module_basis/module_ao/ORB_gen_tables.h" #include "module_cell/module_neighbor/sltk_grid_driver.h" @@ -12,43 +12,42 @@ #include "module_base/parallel_reduce.h" template -const elecstate::DensityMatrix* hamilt::DFTUNew>::dm_in_dftu = nullptr; - -template -hamilt::DFTUNew>::DFTUNew( +hamilt::DFTU>::DFTU( LCAO_Matrix* LM_in, const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, std::vector* hK_in, - const UnitCell* ucell_in, + const UnitCell& ucell_in, Grid_Driver* GridD_in, ModuleDFTU::DFTU* dftu_in, - const Parallel_Orbitals* paraV) + const Parallel_Orbitals& paraV) : hamilt::OperatorLCAO(LM_in, kvec_d_in, hR_in, hK_in) { this->cal_type = lcao_dftu; - this->ucell = ucell_in; + this->ucell = &ucell_in; this->dftu = dftu_in; #ifdef __DEBUG assert(this->ucell != nullptr); #endif // initialize HR to allocate sparse Nonlocal matrix memory - this->initialize_HR(GridD_in, paraV); + this->initialize_HR(GridD_in, ¶V); + //set nspin + this->nspin = GlobalV::NSPIN; } // destructor template -hamilt::DFTUNew>::~DFTUNew() +hamilt::DFTU>::~DFTU() { } // initialize_HR() template -void hamilt::DFTUNew>::initialize_HR(Grid_Driver* GridD, +void hamilt::DFTU>::initialize_HR(Grid_Driver* GridD, const Parallel_Orbitals* paraV) { - ModuleBase::TITLE("DFTUNew", "initialize_HR"); - ModuleBase::timer::tick("DFTUNew", "initialize_HR"); + ModuleBase::TITLE("DFTU", "initialize_HR"); + ModuleBase::timer::tick("DFTU", "initialize_HR"); this->adjs_all.clear(); this->adjs_all.reserve(this->ucell->nat); @@ -85,15 +84,15 @@ void hamilt::DFTUNew>::initialize_HR(Grid_Driver* G this->adjs_all.push_back(adjs); } - ModuleBase::timer::tick("DFTUNew", "initialize_HR"); + ModuleBase::timer::tick("DFTU", "initialize_HR"); } template -void hamilt::DFTUNew>::cal_nlm_all(const Parallel_Orbitals* paraV) +void hamilt::DFTU>::cal_nlm_all(const Parallel_Orbitals* paraV) { - ModuleBase::TITLE("DFTUNew", "cal_nlm_all"); + ModuleBase::TITLE("DFTU", "cal_nlm_all"); if(this->precal_nlm_done) return; - ModuleBase::timer::tick("DFTUNew", "cal_nlm_all"); + ModuleBase::timer::tick("DFTU", "cal_nlm_all"); nlm_tot.resize(this->ucell->nat); const int npol = this->ucell->get_npol(); int atom_index = 0; @@ -129,30 +128,28 @@ void hamilt::DFTUNew>::cal_nlm_all(const Parallel_O for (int iw1l = 0; iw1l < all_indexes.size(); iw1l += npol) { const int iw1 = all_indexes[iw1l] / npol; + // only first zeta orbitals in target L of atom iat0 are needed + std::vector nlm_target(tlp1); + const int L1 = atom1->iw2l[ iw1 ]; + const int N1 = atom1->iw2n[ iw1 ]; + const int m1 = atom1->iw2m[ iw1 ]; +#ifdef USE_NEW_TWO_CENTER std::vector> nlm; // nlm is a vector of vectors, but size of outer vector is only 1 here // If we are calculating force, we need also to store the gradient // and size of outer vector is then 4 // inner loop : all projectors (L0,M0) -#ifdef USE_NEW_TWO_CENTER //================================================================= // new two-center integral (temporary) //================================================================= - int L1 = atom1->iw2l[ iw1 ]; - int N1 = atom1->iw2n[ iw1 ]; - int m1 = atom1->iw2m[ iw1 ]; // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) - int M1 = (m1 % 2 == 0) ? -m1/2 : (m1+1)/2; + const int M1 = (m1 % 2 == 0) ? -m1/2 : (m1+1)/2; ModuleBase::Vector3 dtau = tau0 - tau1; uot.two_center_bundle->overlap_orb_onsite->snap( T1, L1, N1, M1, T0, dtau * this->ucell->lat0, 0 /*cal_deri*/, nlm); -#else - ModuleBase::WARNING_QUIT("DFTUNew", "old two center integral method not implemented"); -#endif // select the elements of nlm with target_L - std::vector nlm_target(tlp1); for(int iw =0;iw < this->ucell->atoms[T0].nw; iw++) { const int L0 = this->ucell->atoms[T0].iw2l[iw]; @@ -165,28 +162,60 @@ void hamilt::DFTUNew>::cal_nlm_all(const Parallel_O break; } } +#else + ModuleBase::WARNING("DFTU", "autoset onsite_radius to rcut for old two-center integral"); + double olm[3] = {0, 0, 0}; + for(int iw =0;iw < this->ucell->atoms[T0].nw; iw++) + { + const int L0 = this->ucell->atoms[T0].iw2l[iw]; + if(L0 == target_L) + { + for(int m = 0; m < tlp1; m++) + { + uot.snap_psipsi(orb, // orbitals + olm, + 0, + 'S', // olm, job of derivation, dtype of Operator + tau1, + T1, + L1, + m1, + N1, // all zeta of atom1 + tau0, + T0, + L0, + m, + 0 // choose only the first zeta + ); + nlm_target[m] = olm[0]; + } + break; + } + } +#endif nlm_tot[iat0][ad].insert({all_indexes[iw1l], nlm_target}); } } } this->precal_nlm_done = true; - ModuleBase::timer::tick("DFTUNew", "cal_nlm_all"); + ModuleBase::timer::tick("DFTU", "cal_nlm_all"); } +// contributeHR() template -void hamilt::DFTUNew>::calculate_HR() +void hamilt::DFTU>::contributeHR() { - ModuleBase::TITLE("DFTUNew", "calculate_HR"); - if(this->dm_in_dftu == nullptr && this->dftu->initialed_locale == false) + ModuleBase::TITLE("DFTU", "contributeHR"); + if(this->dftu->get_dmr(0) == nullptr && this->dftu->initialed_locale == false) {// skip the calculation if dm_in_dftu is nullptr return; } else { //will update this->dftu->locale and this->dftu->EU - if(GlobalV::CURRENT_SPIN == 0) this->dftu->EU = 0.0; + if(this->current_spin == 0) this->dftu->EU = 0.0; } - ModuleBase::timer::tick("DFTUNew", "calculate_HR"); + ModuleBase::timer::tick("DFTU", "contributeHR"); const Parallel_Orbitals* paraV = this->hR->get_atom_pair(0).get_paraV(); const int npol = this->ucell->get_npol(); @@ -205,13 +234,13 @@ void hamilt::DFTUNew>::calculate_HR() const int tlp1 = 2*target_L+1; AdjacentAtomInfo& adjs = this->adjs_all[atom_index++]; - ModuleBase::timer::tick("DFTUNew", "cal_occupations"); + ModuleBase::timer::tick("DFTU", "cal_occ"); //first iteration to calculate occupation matrix - const int spin_fold = (GlobalV::NSPIN == 4) ? 4 : 1; + const int spin_fold = (this->nspin == 4) ? 4 : 1; std::vector occ(tlp1 * tlp1 * spin_fold, 0.0); if(this->dftu->initialed_locale == false) { - hamilt::HContainer* dmR_current = this->dm_in_dftu->get_DMR_pointer(GlobalV::CURRENT_SPIN+1); + const hamilt::HContainer* dmR_current = this->dftu->get_dmr(this->current_spin); for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) { const int T1 = adjs.ntype[ad1]; @@ -233,7 +262,7 @@ void hamilt::DFTUNew>::calculate_HR() // if not found , skip this pair of atoms if (tmp != nullptr) { - this->cal_occupations(iat1, iat2, T0, paraV, nlm1, nlm2, tmp->get_pointer(), occ); + this->cal_occ(iat1, iat2, paraV, nlm1, nlm2, tmp->get_pointer(), occ); } } } @@ -244,22 +273,22 @@ void hamilt::DFTUNew>::calculate_HR() // save occ to dftu for(int i=0;idftu->locale[iat0][target_L][0][GlobalV::CURRENT_SPIN].c[i] = occ[i]; + if(this->nspin==1) occ[i] *= 0.5; + this->dftu->locale[iat0][target_L][0][this->current_spin].c[i] = occ[i]; } } else // use readin locale to calculate occupation matrix { for(int i=0;idftu->locale[iat0][target_L][0][GlobalV::CURRENT_SPIN].c[i]; + occ[i] = this->dftu->locale[iat0][target_L][0][this->current_spin].c[i]; } // set initialed_locale to false to avoid using readin locale in next iteration } - ModuleBase::timer::tick("DFTUNew", "cal_occupations"); + ModuleBase::timer::tick("DFTU", "cal_occ"); //calculate VU - ModuleBase::timer::tick("DFTUNew", "cal_vu"); + ModuleBase::timer::tick("DFTU", "cal_vu"); const double u_value = this->dftu->U[T0]; std::vector VU_tmp(occ.size()); this->cal_v_of_u(occ, tlp1, u_value, VU_tmp.data(), this->dftu->EU); @@ -291,27 +320,31 @@ void hamilt::DFTUNew>::calculate_HR() // if not found , skip this pair of atoms if (tmp != nullptr) { - this->cal_HR_IJR(iat1, iat2, T0, paraV, nlm1, nlm2, VU, tmp->get_pointer()); + this->cal_HR_IJR(iat1, iat2, paraV, nlm1, nlm2, VU, tmp->get_pointer()); } } } - ModuleBase::timer::tick("DFTUNew", "cal_vu"); + ModuleBase::timer::tick("DFTU", "cal_vu"); } //energy correction for NSPIN=1 - if(GlobalV::NSPIN==1) this->dftu->EU *= 2.0; + if(this->nspin==1) this->dftu->EU *= 2.0; // for readin onsite_dm, set initialed_locale to false to avoid using readin locale in next iteration - if(GlobalV::CURRENT_SPIN == GlobalV::NSPIN-1 || GlobalV::NSPIN==4) this->dftu->initialed_locale = false; + if(this->current_spin == this->nspin-1 || this->nspin==4) this->dftu->initialed_locale = false; - ModuleBase::timer::tick("DFTUNew", "calculate_HR"); + // update this->current_spin: only nspin=2 iterate change it between 0 and 1 + // the key point is only nspin=2 calculate spin-up and spin-down separately, + // and won't calculate spin-up twice without spin-down + if(this->nspin == 2) this->current_spin = 1 - this->current_spin; + + ModuleBase::timer::tick("DFTU", "contributeHR"); } // cal_HR_IJR() template -void hamilt::DFTUNew>::cal_HR_IJR( +void hamilt::DFTU>::cal_HR_IJR( const int& iat1, const int& iat2, - const int& T0, const Parallel_Orbitals* paraV, const std::unordered_map>& nlm1_all, const std::unordered_map>& nlm2_all, @@ -369,10 +402,9 @@ void hamilt::DFTUNew>::cal_HR_IJR( } template -void hamilt::DFTUNew>::cal_occupations( +void hamilt::DFTU>::cal_occ( const int& iat1, const int& iat2, - const int& T0, const Parallel_Orbitals* paraV, const std::unordered_map>& nlm1_all, const std::unordered_map>& nlm2_all, @@ -434,7 +466,7 @@ void hamilt::DFTUNew>::cal_occupations( } template -void hamilt::DFTUNew>::transfer_vu(std::vector& vu_tmp, std::vector& vu) +void hamilt::DFTU>::transfer_vu(std::vector& vu_tmp, std::vector& vu) { #ifdef __DEBUG assert(vu.size() == vu_tmp.size()); @@ -446,7 +478,7 @@ void hamilt::DFTUNew>::transfer_vu(std::vector -void hamilt::DFTUNew, std::complex>>::transfer_vu(std::vector& vu_tmp, std::vector>& vu) +void hamilt::DFTU, std::complex>>::transfer_vu(std::vector& vu_tmp, std::vector>& vu) { #ifdef __DEBUG assert(vu.size() == vu_tmp.size()); @@ -475,12 +507,12 @@ void hamilt::DFTUNew, std::complex -void hamilt::DFTUNew>::cal_v_of_u( +void hamilt::DFTU>::cal_v_of_u( const std::vector& occ, const int m_size, const double u_value, - double* VU, - double& EU) + double* vu, + double& eu) { // calculate the local matrix int spin_fold = occ.size() / m_size / m_size; @@ -492,8 +524,8 @@ void hamilt::DFTUNew>::cal_v_of_u( { for(int m2 = 0; m2 < m_size; m2++) { - VU[start + m1 * m_size + m2] = u_value * (0.5 * (m1 == m2) - occ[start + m2 * m_size + m1]); - EU += u_value * 0.5 * occ[start + m2 * m_size + m1] * occ[start + m1 * m_size + m2]; + vu[start + m1 * m_size + m2] = u_value * (0.5 * (m1 == m2) - occ[start + m2 * m_size + m1]); + eu += u_value * 0.5 * occ[start + m2 * m_size + m1] * occ[start + m1 * m_size + m2]; } } } @@ -503,8 +535,8 @@ void hamilt::DFTUNew>::cal_v_of_u( { for(int m2 = 0; m2 < m_size; m2++) { - VU[m1 * m_size + m2] = u_value * (1.0 * (m1 == m2) - occ[m2 * m_size + m1]); - EU += u_value * 0.25 * occ[m2 * m_size + m1] * occ[m1 * m_size + m2]; + vu[m1 * m_size + m2] = u_value * (1.0 * (m1 == m2) - occ[m2 * m_size + m1]); + eu += u_value * 0.25 * occ[m2 * m_size + m1] * occ[m1 * m_size + m2]; } } for(int is=1;is>::cal_v_of_u( { for(int m2 = 0; m2 < m_size; m2++) { - VU[start + m1 * m_size + m2] = u_value * (0 - occ[start + m2 * m_size + m1]); - EU += u_value * 0.25 * occ[start + m2 * m_size + m1] * occ[start + m1 * m_size + m2]; + vu[start + m1 * m_size + m2] = u_value * (0 - occ[start + m2 * m_size + m1]); + eu += u_value * 0.25 * occ[start + m2 * m_size + m1] * occ[start + m1 * m_size + m2]; } } } } } -// contributeHR() -template -void hamilt::DFTUNew>::contributeHR() -{ - ModuleBase::TITLE("DFTUNew", "contributeHR"); - ModuleBase::timer::tick("DFTUNew", "contributeHR"); - this->calculate_HR(); - ModuleBase::timer::tick("DFTUNew", "contributeHR"); - return; -} - #include "dftu_force_stress.hpp" -template class hamilt::DFTUNew>; -template class hamilt::DFTUNew, double>>; -template class hamilt::DFTUNew, std::complex>>; \ No newline at end of file +template class hamilt::DFTU>; +template class hamilt::DFTU, double>>; +template class hamilt::DFTU, std::complex>>; \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h similarity index 65% rename from source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h rename to source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h index 625244d8f6..1e9cf50ef2 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_lcao.h @@ -1,5 +1,5 @@ -#ifndef DFTUNEW_H -#define DFTUNEW_H +#ifndef DFTPLUSU_H +#define DFTPLUSU_H #include #include "module_basis/module_ao/parallel_orbitals.h" @@ -9,26 +9,12 @@ #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" #include "module_elecstate/module_dm/density_matrix.h" #include "module_hamilt_lcao/module_dftu/dftu.h" +#include "dftu.hpp" namespace hamilt { -#ifndef __DFTUNEWTEMPLATE -#define __DFTUNEWTEMPLATE - -/// The DFTUNew class template inherits from class T -/// it is used to calculate the non-local pseudopotential of wavefunction basis -/// Template parameters: -/// - T: base class, it would be OperatorLCAO or OperatorPW -/// - TR: data type of real space Hamiltonian, it would be double or std::complex -template -class DFTUNew : public T -{ -}; - -#endif - -/// DFTUNew class template specialization for OperatorLCAO base class +/// DFTU class template specialization for OperatorLCAO base class /// It is used to calculate the non-local pseudopotential matrix in real space and fold it to k-space /// HR = D_{p1, p2} /// HK = D_{p1, p2} = \sum_{R} e^{ikR} HR @@ -36,18 +22,18 @@ class DFTUNew : public T /// - TK: data type of k-space Hamiltonian /// - TR: data type of real space Hamiltonian template -class DFTUNew> : public OperatorLCAO +class DFTU> : public OperatorLCAO { public: - DFTUNew>(LCAO_Matrix* LM_in, + DFTU>(LCAO_Matrix* lm_in, const std::vector>& kvec_d_in, hamilt::HContainer* hR_in, std::vector* hK_in, - const UnitCell* ucell_in, - Grid_Driver* GridD_in, + const UnitCell& ucell_in, + Grid_Driver* gridD_in, ModuleDFTU::DFTU* dftu_in, - const Parallel_Orbitals* paraV); - ~DFTUNew>(); + const Parallel_Orbitals& paraV); + ~DFTU>(); /** * @brief contributeHR() is used to calculate the HR matrix @@ -55,9 +41,6 @@ class DFTUNew> : public OperatorLCAO */ virtual void contributeHR() override; - /// pointer of density matrix, it is a temporary implementation - static const elecstate::DensityMatrix* dm_in_dftu; - /// calculate force and stress for DFT+U void cal_force_stress( const bool cal_force, @@ -72,20 +55,19 @@ class DFTUNew> : public OperatorLCAO hamilt::HContainer* HR = nullptr; - hamilt::HContainer* HR_fixed = nullptr; - - bool allocated = false; - TK* HK_pointer = nullptr; - bool HR_fixed_done = false; + /// @brief the number of spin components, 1 for no-spin, 2 for collinear spin case and 4 for non-collinear spin case + int nspin = 0; + /// @brief the current spin index for nspin==2 to calculate spin-up and spin-down separately + int current_spin = 0; /** - * @brief initialize HR, search the nearest neighbor atoms - * HContainer is used to store the non-local pseudopotential matrix with specific atom-pairs - * the size of HR will be fixed after initialization + * @brief search the nearest neighbor atoms and save them into this->adjs_all + * the size of HR will not change in DFTU, + * because I don't want to expand HR larger than Nonlocal operator caused by DFTU */ - void initialize_HR(Grid_Driver* GridD_in, const Parallel_Orbitals* paraV); + void initialize_HR(Grid_Driver* gridD_in, const Parallel_Orbitals* paraV); /** * @brief calculate the overlap values and save them in this->nlm_tot @@ -93,14 +75,16 @@ class DFTUNew> : public OperatorLCAO */ void cal_nlm_all(const Parallel_Orbitals* paraV); - void cal_occupations(const int& iat1, - const int& iat2, - const int& T0, - const Parallel_Orbitals* paraV, - const std::unordered_map>& nlm1_all, - const std::unordered_map>& nlm2_all, - const double* data_pointer, - std::vector& occupations); + /** + * @brief calculate the occ_mm' = \sum_R DMR* matrix for each atom to add U + */ + void cal_occ(const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const std::unordered_map>& nlm1_all, + const std::unordered_map>& nlm2_all, + const double* data_pointer, + std::vector& occupations); /// transfer VU format from pauli matrix to normal for non-collinear spin case void transfer_vu(std::vector& vu_tmp, std::vector& vu); @@ -110,22 +94,14 @@ class DFTUNew> : public OperatorLCAO const std::vector& occ, const int m_size, const double u_value, - double* VU, - double& E); - - /** - * @brief calculate the non-local pseudopotential matrix with specific atom-pairs - * nearest neighbor atoms don't need to be calculated again - * loop the atom-pairs in HR and calculate the non-local pseudopotential matrix - */ - void calculate_HR(); + double* vu, + double& eu); /** * @brief calculate the HR local matrix of atom pair */ void cal_HR_IJR(const int& iat1, const int& iat2, - const int& T0, const Parallel_Orbitals* paraV, const std::unordered_map>& nlm1_all, const std::unordered_map>& nlm2_all, @@ -137,7 +113,6 @@ class DFTUNew> : public OperatorLCAO */ void cal_force_IJR(const int& iat1, const int& iat2, - const int& T0, const Parallel_Orbitals* paraV, const std::unordered_map>& nlm1_all, const std::unordered_map>& nlm2_all, @@ -151,7 +126,6 @@ class DFTUNew> : public OperatorLCAO */ void cal_stress_IJR(const int& iat1, const int& iat2, - const int& T0, const Parallel_Orbitals* paraV, const std::unordered_map>& nlm1_all, const std::unordered_map>& nlm2_all, @@ -163,7 +137,9 @@ class DFTUNew> : public OperatorLCAO double* stress); std::vector adjs_all; + /// @brief if the nlm_tot is calculated bool precal_nlm_done = false; + /// @brief the overlap values for all [atoms][nerghbors][orb_index(iw) in NAOs][m of target_l in Projectors] std::vector>>> nlm_tot; }; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt index 5aeae5314c..798e0caff0 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt @@ -51,6 +51,16 @@ AddTest( tmp_mocks.cpp ../../../../module_hamilt_general/operator.cpp ) +AddTest( + TARGET operator_dftu_test + LIBS ${math_libs} psi base device numerical_atomic_orbitals container + SOURCES test_dftu.cpp ../dftu_lcao.cpp ../../../module_hcontainer/func_folding.cpp + ../../../module_hcontainer/base_matrix.cpp ../../../module_hcontainer/hcontainer.cpp ../../../module_hcontainer/atom_pair.cpp + ../../../../module_basis/module_ao/parallel_2d.cpp ../../../../module_basis/module_ao/parallel_orbitals.cpp + ../../../../module_basis/module_ao/ORB_atomic_lm.cpp + tmp_mocks.cpp ../../../../module_hamilt_general/operator.cpp +) + install(FILES parallel_operator_tests.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) find_program(BASH bash) add_test(NAME operators_para_test diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp new file mode 100644 index 0000000000..a62c409b47 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp @@ -0,0 +1,299 @@ +#include "gtest/gtest.h" +#include "../dftu_lcao.h" +#include + +// mock of DFTU +#include "module_hamilt_lcao/module_dftu/dftu.h" +ModuleDFTU::DFTU::DFTU(){}; +ModuleDFTU::DFTU::~DFTU(){}; +namespace GlobalC +{ + ModuleDFTU::DFTU dftu; +} +const hamilt::HContainer* tmp_DMR; +const hamilt::HContainer* ModuleDFTU::DFTU::get_dmr(int ispin) const +{ + return tmp_DMR; +} + + + +//--------------------------------------- +// Unit test of DFTU class +// DFTU is a derivative class of Operator, it is used to calculate the kinetic matrix +// It use HContainer to store the real space HR matrix +// In this test, we test the correctness and time consuming of 3 functions in DFTU class +// - initialize_HR() called in constructor +// - contributeHR() +// - contributeHk() +// - HR(double) and SK(complex) are tested in constructHRd2cd +// - HR(double) and SK(double) are tested in constructHRd2d +//--------------------------------------- + +// test_size is the number of atoms in the unitcell +// modify test_size to test different size of unitcell +int test_size = 10; +int test_nw = 10; // please larger than 5 +class DFTUTest : public ::testing::Test +{ + protected: + void SetUp() override + { +#ifdef __MPI + // MPI parallel settings + MPI_Comm_size(MPI_COMM_WORLD, &dsize); + MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); +#endif + + // set up a unitcell, with one element and test_size atoms, each atom has test_nw orbitals + ucell.ntype = 1; + ucell.nat = test_size; + ucell.atoms = new Atom[ucell.ntype]; + ucell.iat2it = new int[ucell.nat]; + ucell.iat2ia = new int[ucell.nat]; + ucell.atoms[0].tau = new ModuleBase::Vector3[ucell.nat]; + ucell.lat0 = 1.0; + ucell.itia2iat.create(ucell.ntype, ucell.nat); + for (int iat = 0; iat < ucell.nat; iat++) + { + ucell.iat2it[iat] = 0; + ucell.iat2ia[iat] = iat; + ucell.atoms[0].tau[iat] = ModuleBase::Vector3(0.0, 0.0, 0.0); + ucell.itia2iat(0, iat) = iat; + } + ucell.atoms[0].na = test_size; + ucell.atoms[0].nw = test_nw; + ucell.atoms[0].iw2l = new int[test_nw]; + ucell.atoms[0].iw2m = new int[test_nw]; + ucell.atoms[0].iw2n = new int[test_nw]; + for (int iw = 0; iw < test_nw; ++iw) + { + ucell.atoms[0].iw2l[iw] = 2; + ucell.atoms[0].iw2m[iw] = 0; + ucell.atoms[0].iw2n[iw] = 0; + } + ucell.set_iat2iwt(1); + init_parav(); + // set up a HContainer with ucell + HR = new hamilt::HContainer(ucell, paraV); + // initialize DMR and set default values of 1.0 + DMR = new hamilt::HContainer(*HR); + tmp_DMR = DMR; + + //setting of DFTU + GlobalC::dftu.locale.resize(test_size); + for(int iat=0;iatset_block_size(10/* nb_2d set to be 2*/); + paraV->set_proc_dim(dsize, 0); + paraV->mpi_create_cart(MPI_COMM_WORLD); + paraV->set_local2global(global_row, global_col, ofs_running, ofs_running); + int lr = paraV->get_row_size(); + int lc = paraV->get_col_size(); + paraV->set_desc(global_row, global_col, lr); + paraV->set_global2local(global_row, global_col, true, ofs_running); + paraV->set_atomic_trace(ucell.get_iat2iwt(), test_size, global_row); + } +#else + void init_parav() + {} +#endif + + UnitCell ucell; + hamilt::HContainer* HR; + hamilt::HContainer* DMR; + Parallel_Orbitals *paraV; + + int dsize; + int my_rank = 0; + double U_test = 1.0; + int orbital_c_test = 2; +}; + +// using TEST_F to test DFTU +TEST_F(DFTUTest, constructHRd2d) +{ + //test for nspin=1 + GlobalV::NSPIN = 1; + std::vector> kvec_d_in(1, ModuleBase::Vector3(0.0, 0.0, 0.0)); + std::vector hk(paraV->get_row_size() * paraV->get_col_size(), 0.0); + Grid_Driver gd(0,0,0); + // check some input values + EXPECT_EQ(LCAO_Orbitals::get_const_instance().Phi[0].getRcut(), 1.0); + // reset HR and DMR + const double factor = 1.0 / test_nw / test_nw / test_size / test_size; + for(int i=0;iget_nnr();i++) + { + DMR->get_wrapper()[i] = factor; + HR->get_wrapper()[i] = 0.0; + } + std::chrono::high_resolution_clock::time_point start_time = std::chrono::high_resolution_clock::now(); + hamilt::DFTU> op( + nullptr, + kvec_d_in, + HR, + &hk, + ucell, + &gd, + &GlobalC::dftu, + *paraV + ); + std::chrono::high_resolution_clock::time_point end_time = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed_time = std::chrono::duration_cast>(end_time - start_time); + start_time = std::chrono::high_resolution_clock::now(); + op.contributeHR(); + end_time = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed_time1 = std::chrono::duration_cast>(end_time - start_time); + // check the occupations of dftu + for(int iat=0;iatsize_atom_pairs(); ++iap) + { + hamilt::AtomPair& tmp = HR->get_atom_pair(iap); + int iat1 = tmp.get_atom_i(); + int iat2 = tmp.get_atom_j(); + auto indexes1 = paraV->get_indexes_row(iat1); + auto indexes2 = paraV->get_indexes_col(iat2); + int nwt = indexes1.size() * indexes2.size(); + for (int i = 0; i < nwt; ++i) + { + EXPECT_NEAR(tmp.get_pointer(0)[i], -10.0*test_size, 1e-10); + } + } + // calculate SK + start_time = std::chrono::high_resolution_clock::now(); + op.contributeHk(0); + end_time = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed_time2 = std::chrono::duration_cast>(end_time - start_time); + // check the value of SK + for (int i = 0; i < paraV->get_row_size() * paraV->get_col_size(); ++i) + { + EXPECT_NEAR(hk[i], -10.0*test_size, 1e-10); + } + std::cout << "Test terms: " <> kvec_d_in(2, ModuleBase::Vector3(0.0, 0.0, 0.0)); + std::vector> hk(paraV->get_row_size() * paraV->get_col_size(), std::complex(0.0, 0.0)); + Grid_Driver gd(0,0,0); + EXPECT_EQ(LCAO_Orbitals::get_const_instance().Phi[0].getRcut(), 1.0); + // reset HR and DMR + const double factor = 0.5 / test_nw / test_nw / test_size / test_size; + for(int i=0;iget_nnr();i++) + { + DMR->get_wrapper()[i] = factor; + HR->get_wrapper()[i] = 0.0; + } + hamilt::DFTU, double>> op( + nullptr, + kvec_d_in, + HR, + &hk, + ucell, + &gd, + &GlobalC::dftu, + *paraV + ); + op.contributeHR(); + // check the occupations of dftu for spin-up + for(int iat=0;iatsize_atom_pairs(); ++iap) + { + hamilt::AtomPair& tmp = HR->get_atom_pair(iap); + int iat1 = tmp.get_atom_i(); + int iat2 = tmp.get_atom_j(); + auto indexes1 = paraV->get_indexes_row(iat1); + auto indexes2 = paraV->get_indexes_col(iat2); + int nwt = indexes1.size() * indexes2.size(); + for (int i = 0; i < nwt; ++i) + { + EXPECT_NEAR(tmp.get_pointer(0)[i], -10.0*test_size, 1e-10); + } + } + // calculate HK for gamma point + op.contributeHk(0); + // check the value of HK of gamma point + for (int i = 0; i < paraV->get_row_size() * paraV->get_col_size(); ++i) + { + EXPECT_NEAR(hk[i].real(), -10.0*test_size, 1e-10); + EXPECT_NEAR(hk[i].imag(), 0.0, 1e-10); + } + // calculate spin-down hamiltonian + op.contributeHR(); + // check the occupations of dftu for spin-down + for(int iat=0;iat, double>* dmr) +{ + this->dm_in_dftu_cd = dmr; + return; +} + +void DFTU::set_dmr(const elecstate::DensityMatrix* dmr) +{ + this->dm_in_dftu_d = dmr; + return; +} + +const hamilt::HContainer* DFTU::get_dmr(int ispin) const +{ + if(this->dm_in_dftu_d != nullptr) + { + return this->dm_in_dftu_d->get_DMR_pointer(ispin+1); + } + else if(this->dm_in_dftu_cd != nullptr) + { + return this->dm_in_dftu_cd->get_DMR_pointer(ispin+1); + } + else + { + return nullptr; + } +} + } // namespace ModuleDFTU diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h index f3d6893e67..e493f47205 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.h +++ b/source/module_hamilt_lcao/module_dftu/dftu.h @@ -12,6 +12,8 @@ #include "module_elecstate/module_charge/charge_mixing.h" #include "module_hamilt_general/hamilt.h" #include "module_elecstate/elecstate.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_elecstate/module_dm/density_matrix.h" #include @@ -186,6 +188,24 @@ class DFTU double spherical_Bessel(const int k, const double r, const double lambda); double spherical_Hankel(const int k, const double r, const double lambda); + + public: + /** + * @brief get the density matrix of target spin + * nspin = 1 and 4 : ispin should be 0 + * nspin = 2 : ispin should be 0/1 + */ + const hamilt::HContainer* get_dmr(int ispin) const; + /** + * @brief set the density matrix for DFT+U calculation + * if the density matrix is not set or set to nullptr, the DFT+U calculation will not be performed + */ + void set_dmr(const elecstate::DensityMatrix* dm_in_dftu_d); + void set_dmr(const elecstate::DensityMatrix, double>* dm_in_dftu_cd); + + private: + const elecstate::DensityMatrix* dm_in_dftu_d = nullptr; + const elecstate::DensityMatrix, double>* dm_in_dftu_cd = nullptr; }; } // namespace ModuleDFTU