From 2f7780bf752bc293d1ad7c1d5459311658cd0bb5 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Fri, 22 Mar 2024 18:20:21 +0800 Subject: [PATCH 1/5] Fix: not use static member in DFTUNew --- source/module_esolver/esolver_ks_lcao.cpp | 10 ++----- .../operator_lcao/dftu_force_stress.hpp | 24 ++++++++++------ .../hamilt_lcaodft/operator_lcao/dftu_new.cpp | 7 ++--- .../hamilt_lcaodft/operator_lcao/dftu_new.h | 3 -- .../module_hamilt_lcao/module_dftu/dftu.cpp | 28 +++++++++++++++++++ source/module_hamilt_lcao/module_dftu/dftu.h | 11 ++++++++ 6 files changed, 58 insertions(+), 25 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 14f355d7f0b..81830f14bad 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -571,15 +571,9 @@ namespace ModuleESolver 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() ); } GlobalC::dftu.cal_slater_UJ(this->pelec->charge->rho, this->pw_rho->nrxx); // Calculate U and J if Yukawa potential is used } 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 1d51b7b11d3..613d36bb8e5 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 @@ -12,13 +12,23 @@ void DFTUNew>::cal_force_stress( ModuleBase::matrix& force, ModuleBase::matrix& stress) { - ModuleBase::TITLE("DFTUNew", "cal_force_stress"); -#ifdef __DEBUG - assert(this->dm_in_dftu != nullptr); -#endif + ModuleBase::TITLE("DFTUNew", "cal_force_stress"); + if(this->dftu->get_dmr(0) == nullptr) + { + ModuleBase::WARNING_QUIT("DFTUNew", "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[GlobalV::NSPIN]; + dmR_tmp[0] = this->dftu->get_dmr(0); + if(GlobalV::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("DFTUNew", "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) @@ -112,10 +122,6 @@ void DFTUNew>::cal_force_stress( const int ii = i % (tlp1 * tlp1); occ[i] = this->dftu->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]; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp index 0426028b0ff..fa71fa4b69c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp @@ -11,9 +11,6 @@ #endif #include "module_base/parallel_reduce.h" -template -const elecstate::DensityMatrix* hamilt::DFTUNew>::dm_in_dftu = nullptr; - template hamilt::DFTUNew>::DFTUNew( LCAO_Matrix* LM_in, @@ -177,7 +174,7 @@ template void hamilt::DFTUNew>::calculate_HR() { ModuleBase::TITLE("DFTUNew", "calculate_HR"); - if(this->dm_in_dftu == nullptr && this->dftu->initialed_locale == false) + if(this->dftu->get_dmr(0) == nullptr && this->dftu->initialed_locale == false) {// skip the calculation if dm_in_dftu is nullptr return; } @@ -211,7 +208,7 @@ void hamilt::DFTUNew>::calculate_HR() 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(GlobalV::CURRENT_SPIN); for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) { const int T1 = adjs.ntype[ad1]; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h index 625244d8f6a..87a844e951c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h @@ -55,9 +55,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, diff --git a/source/module_hamilt_lcao/module_dftu/dftu.cpp b/source/module_hamilt_lcao/module_dftu/dftu.cpp index 37d8d8b035e..d69e755b90d 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.cpp +++ b/source/module_hamilt_lcao/module_dftu/dftu.cpp @@ -355,4 +355,32 @@ void DFTU::cal_energy_correction(const int istep) return; } +void DFTU::set_dmr(const elecstate::DensityMatrix, 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 f3d6893e673..070a401397b 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,15 @@ 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: + const hamilt::HContainer* get_dmr(int ispin) const; + 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 From 986cf86a377cc5100b9e2f6065974b0c807b6266 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Tue, 26 Mar 2024 16:14:49 +0800 Subject: [PATCH 2/5] Test: add unit test of DFTUNew --- .../operator_lcao/dftu_force_stress.hpp | 14 +- .../hamilt_lcaodft/operator_lcao/dftu_new.cpp | 57 +++- .../hamilt_lcaodft/operator_lcao/dftu_new.h | 7 +- .../operator_lcao/test/CMakeLists.txt | 10 + .../operator_lcao/test/test_dftunew.cpp | 291 ++++++++++++++++++ 5 files changed, 354 insertions(+), 25 deletions(-) create mode 100644 source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp 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 613d36bb8e5..e07240dfc17 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 @@ -18,9 +18,9 @@ void DFTUNew>::cal_force_stress( ModuleBase::WARNING_QUIT("DFTUNew", "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[GlobalV::NSPIN]; + const hamilt::HContainer* dmR_tmp[this->nspin]; dmR_tmp[0] = this->dftu->get_dmr(0); - if(GlobalV::NSPIN==2) dmR_tmp[1] = this->dftu->get_dmr(1); + if(this->nspin==2) dmR_tmp[1] = this->dftu->get_dmr(1); if(dmR_tmp[0]->size_atom_pairs() == 0) { return; @@ -115,7 +115,7 @@ 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;i>::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]); } @@ -167,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, T0, 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, T0, paraV, nlm_tot[ad1], nlm_tot[ad2], VU, tmp, this->nspin, dis1, dis2, stress_tmp.data()); } } } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp index fa71fa4b69c..1873a138474 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp @@ -31,6 +31,8 @@ hamilt::DFTUNew>::DFTUNew( #endif // initialize HR to allocate sparse Nonlocal matrix memory this->initialize_HR(GridD_in, paraV); + //set nspin + this->nspin = GlobalV::NSPIN; } // destructor @@ -126,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]; @@ -162,6 +162,37 @@ void hamilt::DFTUNew>::cal_nlm_all(const Parallel_O break; } } +#else + ModuleBase::WARNING("DFTUNew", "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}); } } @@ -204,7 +235,7 @@ void hamilt::DFTUNew>::calculate_HR() ModuleBase::timer::tick("DFTUNew", "cal_occupations"); //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) { @@ -241,7 +272,7 @@ void hamilt::DFTUNew>::calculate_HR() // save occ to dftu for(int i=0;inspin==1) occ[i] *= 0.5; this->dftu->locale[iat0][target_L][0][GlobalV::CURRENT_SPIN].c[i] = occ[i]; } } @@ -296,9 +327,9 @@ void hamilt::DFTUNew>::calculate_HR() } //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(GlobalV::CURRENT_SPIN == this->nspin-1 || this->nspin==4) this->dftu->initialed_locale = false; ModuleBase::timer::tick("DFTUNew", "calculate_HR"); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h index 87a844e951c..8bc4461c954 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h @@ -69,13 +69,10 @@ 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; + int nspin = 0; + int current_spin = 0; /** * @brief initialize HR, search the nearest neighbor atoms 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 5aeae5314c1..fec4caacba4 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_dftunew.cpp ../dftu_new.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_dftunew.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp new file mode 100644 index 00000000000..916d6b459ac --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp @@ -0,0 +1,291 @@ +#include "gtest/gtest.h" +#include "../dftu_new.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 DFTUNew class +// DFTUNew 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 DFTUNew 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 DFTUNewTest : 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; +}; + +// using TEST_F to test DFTUNew +TEST_F(DFTUNewTest, constructHRd2d) +{ + 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::DFTUNew> 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); + //std::cout<<"HR:i "< 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<<"HK:i "<> 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 = 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; + } + hamilt::DFTUNew, double>> op( + nullptr, + kvec_d_in, + HR, + &hk, + &ucell, + &gd, + &GlobalC::dftu, + paraV + ); + op.contributeHR(); + // 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) + { + //std::cout<<"HR0:i "<get_row_size() * paraV->get_col_size(); ++i) + { + //std::cout<<"HK0:i "< Date: Tue, 26 Mar 2024 17:23:43 +0800 Subject: [PATCH 3/5] Test: add nspin=2 test and delete GlobalV::CURRENT_SPIN in DFTUNew --- .../hamilt_lcaodft/operator_lcao/dftu_new.cpp | 15 ++++++--- .../hamilt_lcaodft/operator_lcao/dftu_new.h | 10 ++++-- .../operator_lcao/test/test_dftunew.cpp | 32 ++++++++++++------- 3 files changed, 37 insertions(+), 20 deletions(-) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp index 1873a138474..cdcc7480a39 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.cpp @@ -212,7 +212,7 @@ void hamilt::DFTUNew>::calculate_HR() 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"); @@ -239,7 +239,7 @@ void hamilt::DFTUNew>::calculate_HR() std::vector occ(tlp1 * tlp1 * spin_fold, 0.0); if(this->dftu->initialed_locale == false) { - const hamilt::HContainer* dmR_current = this->dftu->get_dmr(GlobalV::CURRENT_SPIN); + 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]; @@ -273,14 +273,14 @@ void hamilt::DFTUNew>::calculate_HR() for(int i=0;inspin==1) occ[i] *= 0.5; - this->dftu->locale[iat0][target_L][0][GlobalV::CURRENT_SPIN].c[i] = occ[i]; + 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 } @@ -329,7 +329,12 @@ void hamilt::DFTUNew>::calculate_HR() //energy correction for NSPIN=1 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 == this->nspin-1 || this->nspin==4) this->dftu->initialed_locale = false; + if(this->current_spin == this->nspin-1 || this->nspin==4) this->dftu->initialed_locale = false; + + // 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("DFTUNew", "calculate_HR"); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h index 8bc4461c954..0dd25786490 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/dftu_new.h @@ -71,13 +71,15 @@ class DFTUNew> : public OperatorLCAO TK* HK_pointer = nullptr; + /// @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 DFTUNew, + * because I don't want to expand HR larger than Nonlocal operator caused by DFTUNew */ void initialize_HR(Grid_Driver* GridD_in, const Parallel_Orbitals* paraV); @@ -157,7 +159,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/test_dftunew.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp index 916d6b459ac..197483f0ef6 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp @@ -88,18 +88,15 @@ class DFTUNewTest : public ::testing::Test for(int l=0;l<3;l++) { GlobalC::dftu.locale[iat][l].resize(1); - GlobalC::dftu.locale[iat][l][0].resize(1); + GlobalC::dftu.locale[iat][l][0].resize(2); GlobalC::dftu.locale[iat][l][0][0].create(2*l+1, 2*l+1); + GlobalC::dftu.locale[iat][l][0][1].create(2*l+1, 2*l+1); } } - double U_test = 1.0; GlobalC::dftu.U = &U_test; - int orbital_c_test = 2; GlobalC::dftu.orbital_corr = &orbital_c_test; GlobalV::onsite_radius = 1.0; - GlobalV::NSPIN = 1; - GlobalV::CURRENT_SPIN = 0; } @@ -146,11 +143,15 @@ class DFTUNewTest : public ::testing::Test int dsize; int my_rank = 0; + double U_test = 1.0; + int orbital_c_test = 2; }; // using TEST_F to test DFTUNew TEST_F(DFTUNewTest, 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); @@ -185,7 +186,6 @@ TEST_F(DFTUNewTest, constructHRd2d) { for(int icc=0;icc<25;icc++) { - //std::cout<<"occupations: "<get_row_size() * paraV->get_col_size(); ++i) { EXPECT_NEAR(hk[i], -10.0*test_size, 1e-10); - //std::cout<<"HK:i "<> 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 = 1.0 / test_nw / test_nw / test_size / test_size; + 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; @@ -243,7 +243,7 @@ TEST_F(DFTUNewTest, constructHRd2cd) paraV ); op.contributeHR(); - // check the occupations of dftu + // check the occupations of dftu for spin-up for(int iat=0;iatget_row_size() * paraV->get_col_size(); ++i) { - //std::cout<<"HK0:i "< Date: Mon, 1 Apr 2024 15:48:42 +0800 Subject: [PATCH 4/5] add annotation for dftu.h --- source/module_hamilt_lcao/module_dftu/dftu.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h index 070a401397b..e493f472056 100644 --- a/source/module_hamilt_lcao/module_dftu/dftu.h +++ b/source/module_hamilt_lcao/module_dftu/dftu.h @@ -190,7 +190,16 @@ class DFTU 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); From 8f7c1031fb1277b106654d883e05e96b77f3f05d Mon Sep 17 00:00:00 2001 From: dyzheng Date: Mon, 1 Apr 2024 17:48:16 +0800 Subject: [PATCH 5/5] Refactor: new dftu code with suggestions by developers --- source/Makefile.Objects | 2 +- source/module_esolver/esolver_ks_lcao.cpp | 1 - .../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 | 24 +++-- .../{dftu_new.cpp => dftu_lcao.cpp} | 98 ++++++++----------- .../operator_lcao/{dftu_new.h => dftu_lcao.h} | 72 +++++--------- .../operator_lcao/test/CMakeLists.txt | 2 +- .../test/{test_dftunew.cpp => test_dftu.cpp} | 28 +++--- 12 files changed, 126 insertions(+), 145 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} (87%) rename source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/{dftu_new.h => dftu_lcao.h} (72%) rename source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/{test_dftunew.cpp => test_dftu.cpp} (94%) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index e03f39948e4..eb25c6f7f9f 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 40a082e3103..0f581f2f4a6 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 diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt index 10b1df34bde..e9de2d14caa 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 4fa805e78de..901417a7a41 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) @@ -235,14 +235,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 9701e4c8f0e..ea263761829 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 e1407f65405..58950a5d305 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 00000000000..3f8478f7abc --- /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 e07240dfc17..cec70d48cee 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,21 +1,21 @@ #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"); + ModuleBase::TITLE("DFTU", "cal_force_stress"); if(this->dftu->get_dmr(0) == nullptr) { - ModuleBase::WARNING_QUIT("DFTUNew", "dmr is not set"); + 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]; @@ -26,7 +26,7 @@ void DFTUNew>::cal_force_stress( return; } // begin the calculation of force and stress - ModuleBase::timer::tick("DFTUNew", "cal_force_stress"); + ModuleBase::timer::tick("DFTU", "cal_force_stress"); const Parallel_Orbitals* paraV = dmR_tmp[0]->get_atom_pair(0).get_paraV(); const int npol = this->ucell->get_npol(); @@ -91,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); @@ -167,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, this->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, this->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()); } } } @@ -208,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, @@ -278,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 87% 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 cdcc7480a39..062895d628d 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,42 +12,42 @@ #include "module_base/parallel_reduce.h" 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); @@ -84,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; @@ -163,7 +163,7 @@ void hamilt::DFTUNew>::cal_nlm_all(const Parallel_O } } #else - ModuleBase::WARNING("DFTUNew", "autoset onsite_radius to rcut for old two-center integral"); + 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++) { @@ -198,13 +198,14 @@ void hamilt::DFTUNew>::cal_nlm_all(const Parallel_O } } 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"); + 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; @@ -214,7 +215,7 @@ void hamilt::DFTUNew>::calculate_HR() //will update this->dftu->locale and this->dftu->EU 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(); @@ -233,7 +234,7 @@ 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 = (this->nspin == 4) ? 4 : 1; std::vector occ(tlp1 * tlp1 * spin_fold, 0.0); @@ -261,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); } } } @@ -284,10 +285,10 @@ void hamilt::DFTUNew>::calculate_HR() } // 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); @@ -319,11 +320,11 @@ 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 @@ -336,15 +337,14 @@ void hamilt::DFTUNew>::calculate_HR() // and won't calculate spin-up twice without spin-down if(this->nspin == 2) this->current_spin = 1 - this->current_spin; - ModuleBase::timer::tick("DFTUNew", "calculate_HR"); + 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, @@ -402,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, @@ -467,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()); @@ -479,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()); @@ -508,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; @@ -525,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]; } } } @@ -536,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 72% 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 0dd25786490..1e9cf50ef2d 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 @@ -78,10 +64,10 @@ class DFTUNew> : public OperatorLCAO /** * @brief search the nearest neighbor atoms and save them into this->adjs_all - * the size of HR will not change in DFTUNew, - * because I don't want to expand HR larger than Nonlocal operator caused by DFTUNew + * 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 @@ -89,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); @@ -106,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, @@ -133,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, @@ -147,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, 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 fec4caacba4..798e0caff09 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 @@ -54,7 +54,7 @@ AddTest( AddTest( TARGET operator_dftu_test LIBS ${math_libs} psi base device numerical_atomic_orbitals container - SOURCES test_dftunew.cpp ../dftu_new.cpp ../../../module_hcontainer/func_folding.cpp + 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 diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp similarity index 94% rename from source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp rename to source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp index 197483f0ef6..a62c409b47d 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftunew.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_dftu.cpp @@ -1,5 +1,5 @@ #include "gtest/gtest.h" -#include "../dftu_new.h" +#include "../dftu_lcao.h" #include // mock of DFTU @@ -19,10 +19,10 @@ const hamilt::HContainer* ModuleDFTU::DFTU::get_dmr(int ispin) const //--------------------------------------- -// Unit test of DFTUNew class -// DFTUNew is a derivative class of Operator, it is used to calculate the kinetic matrix +// 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 DFTUNew class +// In this test, we test the correctness and time consuming of 3 functions in DFTU class // - initialize_HR() called in constructor // - contributeHR() // - contributeHk() @@ -34,7 +34,7 @@ const hamilt::HContainer* ModuleDFTU::DFTU::get_dmr(int ispin) const // modify test_size to test different size of unitcell int test_size = 10; int test_nw = 10; // please larger than 5 -class DFTUNewTest : public ::testing::Test +class DFTUTest : public ::testing::Test { protected: void SetUp() override @@ -147,8 +147,8 @@ class DFTUNewTest : public ::testing::Test int orbital_c_test = 2; }; -// using TEST_F to test DFTUNew -TEST_F(DFTUNewTest, constructHRd2d) +// using TEST_F to test DFTU +TEST_F(DFTUTest, constructHRd2d) { //test for nspin=1 GlobalV::NSPIN = 1; @@ -165,15 +165,15 @@ TEST_F(DFTUNewTest, constructHRd2d) HR->get_wrapper()[i] = 0.0; } std::chrono::high_resolution_clock::time_point start_time = std::chrono::high_resolution_clock::now(); - hamilt::DFTUNew> op( + hamilt::DFTU> op( nullptr, kvec_d_in, HR, &hk, - &ucell, + ucell, &gd, &GlobalC::dftu, - paraV + *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); @@ -217,7 +217,7 @@ TEST_F(DFTUNewTest, constructHRd2d) std::cout << "Elapsed time: " <get_wrapper()[i] = factor; HR->get_wrapper()[i] = 0.0; } - hamilt::DFTUNew, double>> op( + hamilt::DFTU, double>> op( nullptr, kvec_d_in, HR, &hk, - &ucell, + ucell, &gd, &GlobalC::dftu, - paraV + *paraV ); op.contributeHR(); // check the occupations of dftu for spin-up