diff --git a/source/module_basis/module_ao/ORB_gen_tables.cpp b/source/module_basis/module_ao/ORB_gen_tables.cpp index 40c8bdbc882..5137344e238 100644 --- a/source/module_basis/module_ao/ORB_gen_tables.cpp +++ b/source/module_basis/module_ao/ORB_gen_tables.cpp @@ -13,6 +13,11 @@ ORB_gen_tables UOT; ORB_gen_tables::ORB_gen_tables() {} ORB_gen_tables::~ORB_gen_tables() {} +const ORB_gen_tables& ORB_gen_tables::get_const_instance() +{ + return GlobalC::UOT; +} + /// call in hamilt_linear::init_before_ions. void ORB_gen_tables::gen_tables( std::ofstream &ofs_in, @@ -448,8 +453,6 @@ void ORB_gen_tables::snap_psipsi( const int &L2, const int &m2, const int &N2, - const int &nspin, - std::complex *olm1, bool cal_syns, double dmax) const { @@ -614,21 +617,7 @@ void ORB_gen_tables::snap_psipsi( { case 0: // calculate overlap. { - if (nspin != 4) - { - olm[0] += tmpOlm0 * rly[MGT.get_lm_index(L, m)]; - } - else if (olm1 != NULL) - { - olm1[0] += tmpOlm0 * rly[MGT.get_lm_index(L, m)]; - olm1[1] += 0; //tmpOlm0 * (tmp(0,0)+tmp(0,1)); - olm1[2] += 0; //tmpOlm0 * (tmp(1,0)+tmp(1,1)); - olm1[3] += tmpOlm0 * rly[MGT.get_lm_index(L, m)]; - } - else - { - ModuleBase::WARNING_QUIT("ORB_gen_tables::snap_psipsi", "something wrong!"); - } + olm[0] += tmpOlm0 * rly[MGT.get_lm_index(L, m)]; /* if( abs ( tmpOlm0 * rly[ MGT.get_lm_index(L, m) ] ) > 1.0e-3 ) @@ -710,21 +699,7 @@ void ORB_gen_tables::snap_psipsi( { case 0: { - if (nspin != 4) - { - olm[0] += tmpKem0 * rly[MGT.get_lm_index(L, m)]; - } - else if (olm1 != NULL) - { - olm1[0] += tmpKem0 * rly[MGT.get_lm_index(L, m)]; - olm1[1] += 0; //tmpKem0 * (tmp(0,0)+tmp(0,1)); - olm1[2] += 0; //tmpKem0 * (tmp(1,0)+tmp(1,1)); - olm1[3] += tmpKem0 * rly[MGT.get_lm_index(L, m)]; - } - else - { - ModuleBase::WARNING_QUIT("ORB_gen_tables::snap_psipsi", "something wrong in T."); - } + olm[0] += tmpKem0 * rly[MGT.get_lm_index(L, m)]; break; } case 1: diff --git a/source/module_basis/module_ao/ORB_gen_tables.h b/source/module_basis/module_ao/ORB_gen_tables.h index 012e550d2ca..be888d77d0f 100644 --- a/source/module_basis/module_ao/ORB_gen_tables.h +++ b/source/module_basis/module_ao/ORB_gen_tables.h @@ -23,6 +23,9 @@ class ORB_gen_tables ORB_gen_tables(); ~ORB_gen_tables(); + // static function to get global instance + static const ORB_gen_tables& get_const_instance(); + void gen_tables( std::ofstream &ofs_in, // mohan add 2021-05-07 LCAO_Orbitals &orb, @@ -49,8 +52,6 @@ class ORB_gen_tables const int &l2, const int &m2, const int &n2, - const int &nspin, - std::complex *olm1=NULL, bool cal_syns = false, double dmax = 0.0)const; diff --git a/source/module_basis/module_ao/ORB_read.cpp b/source/module_basis/module_ao/ORB_read.cpp index d66ceadd631..167390648ad 100644 --- a/source/module_basis/module_ao/ORB_read.cpp +++ b/source/module_basis/module_ao/ORB_read.cpp @@ -38,6 +38,11 @@ LCAO_Orbitals::~LCAO_Orbitals() delete[] Alpha; } +const LCAO_Orbitals& LCAO_Orbitals::get_const_instance() +{ + return GlobalC::ORB; +} + #ifdef __MPI // be called in UnitCell. void LCAO_Orbitals::bcast_files( diff --git a/source/module_basis/module_ao/ORB_read.h b/source/module_basis/module_ao/ORB_read.h index 5e9eaa19797..487b8aa3cac 100644 --- a/source/module_basis/module_ao/ORB_read.h +++ b/source/module_basis/module_ao/ORB_read.h @@ -22,6 +22,9 @@ class LCAO_Orbitals LCAO_Orbitals(); ~LCAO_Orbitals(); + // static function to get global instance + static const LCAO_Orbitals& get_const_instance(); + void Read_Orbitals( std::ofstream &ofs_in, // mohan add 2021-05-07 const int &ntype_in, diff --git a/source/module_basis/module_ao/parallel_orbitals.cpp b/source/module_basis/module_ao/parallel_orbitals.cpp index 3ea9ec0535b..2f03ff9c57e 100644 --- a/source/module_basis/module_ao/parallel_orbitals.cpp +++ b/source/module_basis/module_ao/parallel_orbitals.cpp @@ -29,36 +29,41 @@ Parallel_Orbitals::~Parallel_Orbitals() void Parallel_Orbitals::set_atomic_trace(const int* iat2iwt, const int &nat, const int &nlocal) { - this->atom_begin_col.resize(nat); - this->atom_begin_row.resize(nat); - for(int iat=0;iatiat2iwt_ = iat2iwt; + int nat_plus_1 = nat + 1; + this->atom_begin_col.resize(nat_plus_1); + this->atom_begin_row.resize(nat_plus_1); + for(int iat=0;iatatom_begin_col[iat] = -1; this->atom_begin_row[iat] = -1; int irow = iat2iwt[iat]; int icol = iat2iwt[iat]; - const int max = (iat == nat-1) ? (nlocal - irow): (iat2iwt[iat+1] - irow); - //find the first row index of atom iat - for(int i=0;iglobal2local_row_[irow] != -1) { - this->atom_begin_row[iat] = irow; + this->atom_begin_row[iat] = this->global2local_row_[irow]; break; } irow++; } - //find the first col index of atom iat - for(int i=0;iglobal2local_col_[icol] != -1) { - this->atom_begin_col[iat] = icol; + this->atom_begin_col[iat] = this->global2local_col_[icol]; break; } icol++; } } + this->atom_begin_row[nat] = this->nrow; + this->atom_begin_col[nat] = this->ncol; } // Get the number of columns of the parallel orbital matrix @@ -103,7 +108,7 @@ int Parallel_Orbitals::get_row_size(int iat) const return 0; } iat += 1; - while(this->atom_begin_row[iat] <= this->ncol) + while(this->atom_begin_row[iat] <= this->nrow) { if(this->atom_begin_row[iat] != -1) { @@ -116,6 +121,77 @@ int Parallel_Orbitals::get_row_size(int iat) const throw std::string("error in get_col_size(iat)"); } +// Get the global indexes of the rows of the parallel orbital matrix +std::vector Parallel_Orbitals::get_indexes_row() const +{ + std::vector indexes(this->nrow); + for(int i = 0; i < this->nrow; i++) + { +#ifdef __MPI + indexes[i] = this->local2global_row(i); +#else + indexes[i] = i; +#endif + } + return indexes; +} +// Get the global indexes of the columns of the parallel orbital matrix +std::vector Parallel_Orbitals::get_indexes_col() const +{ + std::vector indexes(this->ncol); + for(int i = 0; i < this->ncol; i++) + { +#ifdef __MPI + indexes[i] = this->local2global_col(i); +#else + indexes[i] = i; +#endif + } + return indexes; +} +// Get the global indexes of the rows of the orbital matrix of the iat-th atom +std::vector Parallel_Orbitals::get_indexes_row(int iat) const +{ + int size = this->get_row_size(iat); + if(size == 0) + { + return std::vector(); + } + std::vector indexes(size); + int irow = this->atom_begin_row[iat]; + int begin = this->iat2iwt_[iat]; + for(int i = 0; i < size; ++i) + { +#ifdef __MPI + indexes[i] = this->local2global_row(irow + i) - begin; +#else + indexes[i] = i; +#endif + } + return indexes; +} +// Get the global indexes of the columns of the orbital matrix of the iat-th atom +std::vector Parallel_Orbitals::get_indexes_col(int iat) const +{ + int size = this->get_col_size(iat); + if(size == 0) + { + return std::vector(); + } + std::vector indexes(size); + int icol = this->atom_begin_col[iat]; + int begin = this->iat2iwt_[iat]; + for(int i = 0; i < size; ++i) + { +#ifdef __MPI + indexes[i] = this->local2global_col(icol + i) - begin; +#else + indexes[i] = i; +#endif + } + return indexes; +} + #ifdef __MPI void Parallel_Orbitals::set_desc_wfc_Eij(const int& nbasis, const int& nbands, const int& lld) { @@ -184,4 +260,4 @@ int Parallel_Orbitals::set_nloc_wfc_Eij( return 0; } -#endif \ No newline at end of file +#endif diff --git a/source/module_basis/module_ao/parallel_orbitals.h b/source/module_basis/module_ao/parallel_orbitals.h index c31ffb1abab..15d0c7c7cff 100644 --- a/source/module_basis/module_ao/parallel_orbitals.h +++ b/source/module_basis/module_ao/parallel_orbitals.h @@ -52,7 +52,15 @@ class Parallel_Orbitals : public Parallel_2D int* loc_sizes; int loc_size; - // set row and col begin index for each atom + /** + * @brief set row and col begin index for each atom + * it should be called after: + * 1. nrow and ncol are set; + * 2. global2local_row_ and global2local_col_ are set; + * @param iat2iwt : the map from atom index to global oribtal indexes + * @param nat : number of atoms + * @param nlocal : number of global orbitals + */ void set_atomic_trace(const int* iat2iwt, const int &nat, const int &nlocal); /** @@ -67,10 +75,24 @@ class Parallel_Orbitals : public Parallel_2D int get_col_size(int iat) const; int get_row_size(int iat) const; + /** + * @brief gather global indexes of orbitals in this processor + * get_indexes_row() : global indexes (~NLOCAL) of rows of Hamiltonian matrix in this processor + * get_indexes_col() : global indexes (~NLOCAL) of columns of Hamiltonian matrix in this processor + * get_indexes_row(iat) : global indexes (~nw) of rows of Hamiltonian matrix in atom iat + * get_indexes_col(iat) : global indexes (~nw) of columns of Hamiltonian matrix in atom iat + */ + std::vector get_indexes_row() const; + std::vector get_indexes_col() const; + std::vector get_indexes_row(int iat) const; + std::vector get_indexes_col(int iat) const; + // private: // orbital index for each atom std::vector atom_begin_row; std::vector atom_begin_col; + const int* iat2iwt_ = nullptr; + }; #endif \ No newline at end of file diff --git a/source/module_basis/module_ao/test/1_snap_equal_test.cpp b/source/module_basis/module_ao/test/1_snap_equal_test.cpp index bc68cf83340..95e99fd3f6e 100644 --- a/source/module_basis/module_ao/test/1_snap_equal_test.cpp +++ b/source/module_basis/module_ao/test/1_snap_equal_test.cpp @@ -45,13 +45,13 @@ TEST_F(test_orb, equal_test) OGT.snap_psipsi( ORB, olm_0, 0, 'S', R1, T1, L1, m1, N1, - R2, T2, L2, m2, N2, - 1, NULL); + R2, T2, L2, m2, N2 + ); OGT.snap_psipsi( ORB, olm_1, 1, 'S', R1, T1, L1, m1, N1, - R2, T2, L2, m2, N2, - 1, NULL); + R2, T2, L2, m2, N2 + ); //std::cout << this->mock_center2_orb11[T1][T2][L1][N1][L2][N2]->cal_overlap(R1, R2, m1, m2); clm_0 = test_center2_orb11[T1][T2][L1][N1][L2][N2]->cal_overlap(R1, R2, m1, m2); diff --git a/source/module_basis/module_ao/test/CMakeLists.txt b/source/module_basis/module_ao/test/CMakeLists.txt index baddcfa0529..8234a45c21e 100644 --- a/source/module_basis/module_ao/test/CMakeLists.txt +++ b/source/module_basis/module_ao/test/CMakeLists.txt @@ -114,6 +114,12 @@ AddTest( LIBS ${math_libs} device base ) +AddTest( + TARGET parallel_orbitals_test + SOURCES parallel_orbitals_test.cpp ../parallel_2d.cpp ../parallel_orbitals.cpp + LIBS ${math_libs} device base +) + install(DIRECTORY lcao_H2O DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) install(DIRECTORY lcao_H2O DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/../../../tests) diff --git a/source/module_basis/module_ao/test/parallel_orbitals_test.cpp b/source/module_basis/module_ao/test/parallel_orbitals_test.cpp new file mode 100644 index 00000000000..a51eddeeeb2 --- /dev/null +++ b/source/module_basis/module_ao/test/parallel_orbitals_test.cpp @@ -0,0 +1,223 @@ +#include "gtest/gtest.h" +#include "../parallel_orbitals.h" + +//-------------------------------------------------------------- +// unit test of class "Parallel_Orbitals" +//-------------------------------------------------------------- +/** + * Test functions: + * - set_atomic_trace + * - get_col_size + * - get_row_size + * - get_col_size(iat) + * - get_row_size(iat) + * - get_indexes_row + * - get_indexes_col + * - get_indexes_row(iat) + * - get_indexes_col(iat) + * + * the test framework is based on parallel_2d_test.cpp +*/ +class TestParaO : public testing::Test +{ +protected: + int dsize; + int my_rank = 0; + std::vector> sizes{ {50, 50} , {60, 60}}; + std::vector nat{ 10, 5}; + std::vector nbs{ 1,2,3 }; + std::ofstream ofs_running; +#ifdef __MPI + void SetUp() override + { + MPI_Comm_size(MPI_COMM_WORLD, &dsize); + MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); + this->ofs_running.open("log" + std::to_string(my_rank) + ".txt"); + ofs_running << "dsize(nproc) = " << dsize << std::endl; + ofs_running << "my_rank = " << my_rank << std::endl; + } + void TearDown() override + { + ofs_running.close(); + } +#endif +}; + +#ifdef __MPI +TEST_F(TestParaO, Divide2D) +{ + for (auto& size : sizes) + { + int gr = size.first; + int gc = size.second; + for (auto nb : nbs) + { + Parallel_Orbitals po; + po.set_block_size(nb); + EXPECT_EQ(po.get_block_size(), nb); + + for (auto mode : { 0,1 }) + { + //1. set dim0 and dim1 + po.set_proc_dim(dsize, mode); + EXPECT_EQ(po.dim0 * po.dim1, dsize); + if (mode)EXPECT_LE(po.dim1, po.dim0); + else EXPECT_LE(po.dim0, po.dim1); + + //2. mpi_create_cart + po.mpi_create_cart(MPI_COMM_WORLD); + EXPECT_NE(po.comm_2D, MPI_COMM_NULL); + + //3. set_local2global and local sizes + po.set_local2global(gr, gc, ofs_running, ofs_running); + int lr = po.get_row_size(); + int lc = po.get_col_size(); + EXPECT_EQ(lr * lc, po.get_local_size()); + auto cal_lsize = [](const int& gsize, const int& nb, const int& np, const int& pcoord) -> int + { + int nblock = gsize / nb; + return nblock / np * nb + static_cast(nblock % np > pcoord) * nb //full blocks' contribution + + static_cast(nblock % np == pcoord) * (gsize % nb); // the last block's contribution + }; + EXPECT_EQ(lr, cal_lsize(gr, nb, po.dim0, po.coord[0])); + EXPECT_EQ(lc, cal_lsize(gc, nb, po.dim1, po.coord[1])); + + //4. set_desc + po.set_desc(gr, gc, lr); + EXPECT_EQ(po.desc[0], 1); + EXPECT_EQ(po.desc[1], po.blacs_ctxt); + EXPECT_EQ(po.desc[2], gr); + EXPECT_EQ(po.desc[3], gc); + EXPECT_EQ(po.desc[4], po.get_block_size()); + EXPECT_EQ(po.desc[5], po.get_block_size()); + EXPECT_EQ(po.desc[6], 0); + EXPECT_EQ(po.desc[7], 0); + EXPECT_EQ(po.desc[8], lr); + + //5. set_global2local + po.set_global2local(gr, gc, true, ofs_running); + auto sum_array = [&po](const int& gr, const int& gc) -> std::pair + { + int sum_row = 0; int sum_col = 0; + for (int i = 0; i < gr; ++i) + sum_row += po.global2local_row(i); + for (int i = 0; i < gc; ++i) + sum_col += po.global2local_col(i); + return { sum_row, sum_col }; + }; + std::pair sumrc = sum_array(gr, gc); + EXPECT_EQ(std::get<0>(sumrc), lr * (lr - 1) / 2 - (gr - lr)); + EXPECT_EQ(std::get<1>(sumrc), lc * (lc - 1) / 2 - (gc - lc)); + for (int i = 0;i < lr;++i) + for (int j = 0;j < lc;++j) + EXPECT_TRUE(po.in_this_processor(po.local2global_row(i), po.local2global_col(j))); + + //6. set_atomic_trace + for(auto nat0 : nat) + { + EXPECT_EQ(gr, gc); + std::vector iat2iwt(nat0); + int nw = gr / nat0; + for (int i = 0; i < nat0; ++i) + { + iat2iwt[i] = i * nw; + } + po.set_atomic_trace(iat2iwt.data(), nat0, gr); + auto global_row_array = po.get_indexes_row(); + auto global_col_array = po.get_indexes_col(); + int local_index_trace_row = 0; + int local_index_trace_col = 0; + // check get_col_size(iat) and get_row_size(iat) + for (int i = 0; i < nat0; ++i) + { + auto atomic_row_array = po.get_indexes_row(i); + auto atomic_col_array = po.get_indexes_col(i); + EXPECT_EQ(po.get_col_size(i), atomic_col_array.size()); + EXPECT_EQ(po.get_row_size(i), atomic_row_array.size()); + for (int j = 0; j < atomic_row_array.size(); ++j) + { + //check global_index == global_index + EXPECT_EQ(atomic_row_array[j]+iat2iwt[i], global_row_array[local_index_trace_row]); + //check local_index == local_index + EXPECT_EQ(local_index_trace_row, po.global2local_row(atomic_row_array[j]+iat2iwt[i])); + local_index_trace_row++; + } + for (int j = 0; j < atomic_col_array.size(); ++j) + { + //check global_index == global_index + EXPECT_EQ(atomic_col_array[j]+iat2iwt[i], global_col_array[local_index_trace_col]); + //check local_index == local_index + EXPECT_EQ(local_index_trace_col, po.global2local_col(atomic_col_array[j]+iat2iwt[i])); + local_index_trace_col++; + } + + } + } + } + } + } +} +#else +TEST_F(TestParaO, Serial) +{ + for (auto& size : sizes) + { + int gr = size.first; + int gc = size.second; + + Parallel_Orbitals po; + + //1. set dim0 and dim1 + po.set_proc_dim(1); + EXPECT_EQ(po.dim0 * po.dim1, 1); + + //2. set_serial + po.set_serial(gr, gc); + EXPECT_EQ(po.get_row_size(), gr); + EXPECT_EQ(po.get_col_size(), gc); + EXPECT_EQ(po.get_local_size(), gr * gc); + + //3. set_global2local + po.set_global2local(gr, gc, false, ofs_running); + for (int i = 0;i < gr;++i) + EXPECT_EQ(po.global2local_row(i), i); + for (int i = 0;i < gc;++i) + EXPECT_EQ(po.global2local_col(i), i); + //6. set_atomic_trace + for(auto nat0 : nat) + { + EXPECT_EQ(gr, gc); + std::vector iat2iwt(nat0); + int nw = gr / nat0; + for (int i = 0; i < nat0; ++i) + { + iat2iwt[i] = i * nw; + } + po.set_atomic_trace(iat2iwt.data(), nat0, gr); + EXPECT_EQ(po.get_col_size(), gr); + EXPECT_EQ(po.get_row_size(), gr); + // check get_col_size(iat) and get_row_size(iat) + for (int i = 0; i < nat0; ++i) + { + std::cout<<__FILE__<<__LINE__<<" i = "<& get_tau(const int &iat) const + { + return atoms[iat2it[iat]].tau[iat2ia[iat]]; + } + + // calculate vector between two atoms with R cell + inline const ModuleBase::Vector3 cal_dtau(const int &iat1, const int &iat2, const ModuleBase::Vector3 &R) const + { + return get_tau(iat2) + double(R.x) * a1 + double(R.y) * a2 + double(R.z) * a3 - get_tau(iat1); + } + //LiuXh add 20180515 ModuleBase::Matrix3 G0; ModuleBase::Matrix3 GT0; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.cpp index 85234f28eea..5240b91fb19 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.cpp @@ -141,19 +141,26 @@ void LCAO_gen_fixedH::build_ST_new(const char& dtype, const bool& calc_deri, con olm[0] = olm[1] = olm[2] = 0.0; - std::complex olm1[4]={ModuleBase::ZERO, ModuleBase::ZERO, ModuleBase::ZERO, ModuleBase::ZERO}; - std::complex *olm2 = &olm1[0]; if(!calc_deri) { // PLEASE use UOT as an input parameter of this subroutine // mohan add 2021-03-30 - GlobalC::UOT.snap_psipsi( GlobalC::ORB, olm, 0, dtype, tau1, - T1, L1, m1, N1, adjs.adjacent_tau[ad], - T2, L2, m2, N2, GlobalV::NSPIN, - olm2,//for soc + GlobalC::UOT.snap_psipsi( GlobalC::ORB, olm, 0, dtype, + tau1, T1, L1, m1, N1, // info of atom1 + adjs.adjacent_tau[ad], T2, L2, m2, N2, // info of atom2 cal_syns, dmax); + // When NSPIN == 4 , only diagonal term is calculated for T or S Operators + // use olm1 to store the diagonal term with complex data type. + std::complex olm1[4]; + if(GlobalV::NSPIN == 4) + { + olm1[0] = std::complex(olm[0], 0.0); + olm1[1] = ModuleBase::ZERO; + olm1[2] = ModuleBase::ZERO; + olm1[3] = std::complex(olm[0], 0.0); + } if(GlobalV::GAMMA_ONLY_LOCAL) { @@ -174,7 +181,8 @@ void LCAO_gen_fixedH::build_ST_new(const char& dtype, const bool& calc_deri, con if (GlobalV::NSPIN != 4) HSloc[nnr] = olm[0]; else {//only has diagonal term here. - int is = (jj-jj0*GlobalV::NPOL) + (kk-kk0*GlobalV::NPOL)*2; + int is = (jj-jj0*GlobalV::NPOL) + (kk-kk0*GlobalV::NPOL)*2; + // SlocR_soc is a temporary array with complex data type, it will be refactor soon. this->LM->SlocR_soc[nnr] = olm1[is]; } } @@ -183,7 +191,8 @@ void LCAO_gen_fixedH::build_ST_new(const char& dtype, const bool& calc_deri, con if(GlobalV::NSPIN!=4) HSloc[nnr] = olm[0];// else {//only has diagonal term here. - int is = (jj-jj0*GlobalV::NPOL) + (kk-kk0*GlobalV::NPOL)*2; + int is = (jj-jj0*GlobalV::NPOL) + (kk-kk0*GlobalV::NPOL)*2; + // Hloc_fixedR_soc is a temporary array with complex data type, it will be refactor soon. this->LM->Hloc_fixedR_soc[nnr] = olm1[is]; } } @@ -195,7 +204,7 @@ void LCAO_gen_fixedH::build_ST_new(const char& dtype, const bool& calc_deri, con { GlobalC::UOT.snap_psipsi( GlobalC::ORB, olm, 1, dtype, tau1, T1, L1, m1, N1, - adjs.adjacent_tau[ad], T2, L2, m2, N2, GlobalV::NSPIN + adjs.adjacent_tau[ad], T2, L2, m2, N2 ); if(GlobalV::GAMMA_ONLY_LOCAL) 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 a735c7d548a..dc3bcd8b979 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/CMakeLists.txt @@ -14,3 +14,7 @@ add_library( if(ENABLE_COVERAGE) add_coverage(operator_ks_lcao) endif() + +IF (BUILD_TESTING) + add_subdirectory(test) +endif() diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h index ca67d29133c..752b052baf1 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/deepks_lcao.h @@ -27,9 +27,12 @@ class DeePKS> : public OperatorLCAO std::vector* HR_pointer_in, std::vector* HK_pointer_in, const int& nks_in) - : loc(loc_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), nks(nks_in), OperatorLCAO(kvec_d_in) + : loc(loc_in), + HR_pointer(HR_pointer_in), + HK_pointer(HK_pointer_in), + nks(nks_in), + OperatorLCAO(LM_in, kvec_d_in) { - this->LM = LM_in; this->cal_type = lcao_deepks; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_lcao.h index a3365b467b2..33e07d39161 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/ekinetic_lcao.h @@ -1,8 +1,8 @@ #ifndef EKINETICLCAO_H #define EKINETICLCAO_H -#include "operator_lcao.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h" #include "module_base/timer.h" +#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h" +#include "operator_lcao.h" namespace hamilt { @@ -10,29 +10,30 @@ namespace hamilt #ifndef __EKINETICTEMPLATE #define __EKINETICTEMPLATE -template class Ekinetic : public T -{}; +template +class Ekinetic : public T +{ +}; #endif -template -class Ekinetic> : public OperatorLCAO +template +class Ekinetic> : public OperatorLCAO { - public: - Ekinetic>(LCAO_gen_fixedH* genH_in, - LCAO_Matrix* LM_in, - const std::vector>& kvec_d_in, - std::vector* HR_pointer_in, - std::vector* HK_pointer_in) - : genH(genH_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), OperatorLCAO(kvec_d_in) - { - this->LM = LM_in; - this->cal_type = lcao_fixed; + public: + Ekinetic>(LCAO_gen_fixedH* genH_in, + LCAO_Matrix* LM_in, + const std::vector>& kvec_d_in, + std::vector* HR_pointer_in, + std::vector* HK_pointer_in) + : genH(genH_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), OperatorLCAO(LM_in, kvec_d_in) + { + this->cal_type = lcao_fixed; } virtual void contributeHR() override; - private: + private: // use overlap matrix to generate fixed Hamiltonian LCAO_gen_fixedH* genH = nullptr; @@ -41,8 +42,7 @@ class Ekinetic> : public OperatorLCAO std::vector* HK_pointer = nullptr; bool HR_fixed_done = false; - }; -} +} // namespace hamilt #endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/meta_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/meta_lcao.h index 913d5a611cb..384688945fc 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/meta_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/meta_lcao.h @@ -1,10 +1,10 @@ #ifndef METALCAO_H #define METALCAO_H -#include "operator_lcao.h" #include "module_base/timer.h" +#include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/module_gint/gint_gamma.h" #include "module_hamilt_lcao/module_gint/gint_k.h" -#include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" +#include "operator_lcao.h" namespace hamilt { @@ -12,34 +12,42 @@ namespace hamilt #ifndef __METATEMPLATE #define __METATEMPLATE -template class Meta : public T -{}; +template +class Meta : public T +{ +}; #endif -template -class Meta> : public OperatorLCAO +template +class Meta> : public OperatorLCAO { - public: - Meta>(Gint_k* GK_in, - Local_Orbital_Charge* loc_in, - LCAO_Matrix* LM_in, - const std::vector>& kvec_d_in, - std::vector* HR_pointer_in, - std::vector* HK_pointer_in) - : GK(GK_in), loc(loc_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), OperatorLCAO(kvec_d_in) - { - this->LM = LM_in; - this->cal_type = lcao_gint; + public: + Meta>(Gint_k* GK_in, + Local_Orbital_Charge* loc_in, + LCAO_Matrix* LM_in, + const std::vector>& kvec_d_in, + std::vector* HR_pointer_in, + std::vector* HK_pointer_in) + : GK(GK_in), + loc(loc_in), + HR_pointer(HR_pointer_in), + HK_pointer(HK_pointer_in), + OperatorLCAO(LM_in, kvec_d_in) + { + this->cal_type = lcao_gint; } - Meta>( - Gint_Gamma* GG_in, - Local_Orbital_Charge* loc_in, - LCAO_Matrix* LM_in, - std::vector* HR_pointer_in, - std::vector* HK_pointer_in):GG(GG_in), loc(loc_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), OperatorLCAO(std::vector>{ModuleBase::Vector3(0,0,0)}) + Meta>(Gint_Gamma* GG_in, + Local_Orbital_Charge* loc_in, + LCAO_Matrix* LM_in, + std::vector* HR_pointer_in, + std::vector* HK_pointer_in) + : GG(GG_in), + loc(loc_in), + HR_pointer(HR_pointer_in), + HK_pointer(HK_pointer_in), + OperatorLCAO(LM_in, std::vector>{ModuleBase::Vector3(0, 0, 0)}) { - this->LM = LM_in; this->cal_type = lcao_gint; } @@ -49,9 +57,7 @@ class Meta> : public OperatorLCAO virtual void contributeHk(int ik) override; - private: - - + private: // used for k-dependent grid integration. Gint_k* GK = nullptr; @@ -66,8 +72,7 @@ class Meta> : public OperatorLCAO std::vector* HK_pointer = nullptr; bool allocated_pvpR = false; - }; -} +} // namespace hamilt #endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_lcao.h index b701c6029ff..4693a2f3cd8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/nonlocal_lcao.h @@ -1,7 +1,7 @@ #ifndef NONLOCALLCAO_H #define NONLOCALLCAO_H -#include "operator_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_gen_fixedH.h" +#include "operator_lcao.h" namespace hamilt { @@ -9,29 +9,30 @@ namespace hamilt #ifndef __NONLOCALTEMPLATE #define __NONLOCALTEMPLATE -template class Nonlocal : public T -{}; +template +class Nonlocal : public T +{ +}; #endif -template -class Nonlocal> : public OperatorLCAO +template +class Nonlocal> : public OperatorLCAO { - public: - Nonlocal>(LCAO_gen_fixedH* genH_in, - LCAO_Matrix* LM_in, - const std::vector>& kvec_d_in, - std::vector* HR_pointer_in, - std::vector* HK_pointer_in) - : genH(genH_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), OperatorLCAO(kvec_d_in) - { - this->LM = LM_in; - this->cal_type = lcao_fixed; + public: + Nonlocal>(LCAO_gen_fixedH* genH_in, + LCAO_Matrix* LM_in, + const std::vector>& kvec_d_in, + std::vector* HR_pointer_in, + std::vector* HK_pointer_in) + : genH(genH_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), OperatorLCAO(LM_in, kvec_d_in) + { + this->cal_type = lcao_fixed; } virtual void contributeHR() override; - private: + private: // use overlap matrix to generate fixed Hamiltonian LCAO_gen_fixedH* genH = nullptr; @@ -40,8 +41,7 @@ class Nonlocal> : public OperatorLCAO std::vector* HK_pointer = nullptr; bool HR_fixed_done = false; - }; -} +} // namespace hamilt #endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h index de2cc47c298..cb5e60cc8ad 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_dftu_lcao.h @@ -1,7 +1,7 @@ #ifndef OPDFTULCAO_H #define OPDFTULCAO_H -#include "operator_lcao.h" #include "module_base/timer.h" +#include "operator_lcao.h" namespace hamilt { @@ -9,32 +9,32 @@ namespace hamilt #ifndef __OPDFTUTEMPLATE #define __OPDFTUTEMPLATE -template class OperatorDFTU : public T -{}; +template +class OperatorDFTU : public T +{ +}; #endif -template -class OperatorDFTU> : public OperatorLCAO +template +class OperatorDFTU> : public OperatorLCAO { - public: - OperatorDFTU>(LCAO_Matrix* LM_in, - const std::vector>& kvec_d_in, - std::vector* HR_pointer_in, - std::vector* HK_pointer_in, - const std::vector& isk_in) - : HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), isk(isk_in), OperatorLCAO(kvec_d_in) - { - this->LM = LM_in; - this->cal_type = lcao_dftu; + public: + OperatorDFTU>(LCAO_Matrix* LM_in, + const std::vector>& kvec_d_in, + std::vector* HR_pointer_in, + std::vector* HK_pointer_in, + const std::vector& isk_in) + : HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), isk(isk_in), OperatorLCAO(LM_in, kvec_d_in) + { + this->cal_type = lcao_dftu; } virtual void contributeHR() override; virtual void contributeHk(int ik) override; - private: - + private: std::vector* HR_pointer = nullptr; std::vector* HK_pointer = nullptr; @@ -43,5 +43,5 @@ class OperatorDFTU> : public OperatorLCAO const std::vector& isk; }; -} +} // namespace hamilt #endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h index 9f5012c0017..f7cc551d737 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h @@ -1,8 +1,8 @@ #ifndef OPEXXLCAO_H #define OPEXXLCAO_H #include "module_base/timer.h" -#include "operator_lcao.h" #include "module_cell/klist.h" +#include "operator_lcao.h" namespace hamilt { @@ -25,9 +25,8 @@ class OperatorEXX> : public OperatorLCAO std::vector* HR_pointer_in, std::vector* HK_pointer_in, const K_Vectors& kv_in) - : HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), kv(kv_in), OperatorLCAO(kv_in.kvec_d) + : HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), kv(kv_in), OperatorLCAO(LM_in, kv_in.kvec_d) { - this->LM = LM_in; this->cal_type = lcao_exx; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h index c611f42054b..c61feec91b6 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h @@ -12,7 +12,8 @@ template class OperatorLCAO : public Operator { public: - OperatorLCAO(const std::vector>& kvec_d_in):kvec_d(kvec_d_in){}; + OperatorLCAO(LCAO_Matrix* LM_in, const std::vector>& kvec_d_in) + : LM(LM_in), kvec_d(kvec_d_in){}; virtual ~OperatorLCAO() { if (this->allocated_smatrix) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_lcao.h index 4ed98391b48..85e91023729 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_lcao.h @@ -23,9 +23,9 @@ class Overlap> : public OperatorLCAO const std::vector>& kvec_d_in, std::vector* SR_pointer_in, std::vector* SK_pointer_in) - : genH(genH_in), SR_pointer(SR_pointer_in), SK_pointer(SK_pointer_in), OperatorLCAO(kvec_d_in) + : genH(genH_in), SR_pointer(SR_pointer_in), SK_pointer(SK_pointer_in), + OperatorLCAO(LM_in, kvec_d_in) { - this->LM = LM_in; this->cal_type = lcao_fixed; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp new file mode 100644 index 00000000000..d03d3fa3be1 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.cpp @@ -0,0 +1,162 @@ +#include "module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h" +#include "module_basis/module_ao/ORB_gen_tables.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" +#include "overlap_new.h" + +template +hamilt::OverlapNew, TR>::OverlapNew(LCAO_Matrix* LM_in, + const std::vector>& kvec_d_in, + hamilt::HContainer* SR_in, + TK* SK_pointer_in, + const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const Parallel_Orbitals* paraV) + : hamilt::OperatorLCAO(LM_in, kvec_d_in) +{ + this->ucell = ucell_in; + this->SR = SR_in; + this->SK_pointer = SK_pointer_in; +#ifdef __DEBUG + assert(this->ucell != nullptr); + assert(this->SR != nullptr); + assert(this->SK_pointer != nullptr); +#endif + // initialize SR to allocate sparse overlap matrix memory + this->initialize_SR(GridD_in, paraV); +} + +// initialize_SR() +template +void hamilt::OverlapNew, TR>::initialize_SR(Grid_Driver* GridD, const Parallel_Orbitals* paraV) +{ + for (int iat1 = 0; iat1 < ucell->nat; iat1++) + { + auto tau1 = ucell->get_tau(iat1); + int T1, I1; + ucell->iat2iait(iat1, &I1, &T1); + AdjacentAtomInfo adjs; + GridD->Find_atom(*ucell, tau1, T1, I1, &adjs); + for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + { + const int T2 = adjs.ntype[ad]; + const int I2 = adjs.natom[ad]; + int iat2 = ucell->itia2iat(T2, I2); + ModuleBase::Vector3& R_index = adjs.box[ad]; + hamilt::AtomPair tmp(iat1, iat2, paraV); + tmp.get_HR_values(R_index.x, R_index.y, R_index.z); + SR->insert_pair(tmp); + } + } +} + +template +void hamilt::OverlapNew, TR>::calculate_SR() +{ +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int iap = 0; iap < this->SR->size_atom_pairs(); ++iap) + { + hamilt::AtomPair& tmp = this->SR->get_atom_pair(iap); + int iat1 = tmp.get_atom_i(); + int iat2 = tmp.get_atom_j(); + const Parallel_Orbitals* paraV = tmp.get_paraV(); + + for (int iR = 0; iR < tmp.get_R_size(); ++iR) + { + const int* R_index = tmp.get_R_index(iR); + ModuleBase::Vector3 R_vector(R_index[0], R_index[1], R_index[2]); + auto dtau = ucell->cal_dtau(iat1, iat2, R_vector); + TR* data_pointer = tmp.get_pointer(iR); + this->cal_SR_IJR(iat1, iat2, paraV, dtau, data_pointer); + } + } +} + +//cal_SR_IJR() +template +void hamilt::OverlapNew, TR>::cal_SR_IJR( + const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const ModuleBase::Vector3& dtau, + TR* data_pointer ) +{ + const ORB_gen_tables& uot = ORB_gen_tables::get_const_instance(); + const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); + // --------------------------------------------- + // get info of orbitals of atom1 and atom2 from ucell + // --------------------------------------------- + int T1, I1; + this->ucell->iat2iait(iat1, &I1, &T1); + int T2, I2; + this->ucell->iat2iait(iat2, &I2, &T2); + Atom& atom1 = this->ucell->atoms[T1]; + Atom& atom2 = this->ucell->atoms[T2]; + int nw1 = atom1.nw; + int nw2 = atom2.nw; + const int* iw2l1 = atom1.iw2l; + const int* iw2n1 = atom1.iw2n; + const int* iw2m1 = atom1.iw2m; + const int* iw2l2 = atom2.iw2l; + const int* iw2n2 = atom2.iw2n; + const int* iw2m2 = atom2.iw2m; + // --------------------------------------------- + // get tau1 (in cell <0,0,0>) and tau2 (in cell R) + // in principle, only dtau is needed in this function + // snap_psipsi should be refactored to use dtau directly + // --------------------------------------------- + const ModuleBase::Vector3& tau1 = this->ucell->get_tau(iat1); + const ModuleBase::Vector3 tau2 = tau1 + dtau; + // --------------------------------------------- + // calculate the overlap matrix for each pair of orbitals + // --------------------------------------------- + double olm[3] = {0, 0, 0}; + auto row_indexes = paraV->get_indexes_row(iat1); + auto col_indexes = paraV->get_indexes_col(iat2); + for (int iw1l = 0; iw1l < row_indexes.size(); ++iw1l) + { + const int iw1 = row_indexes[iw1l]; + const int L1 = iw2l1[iw1]; + const int N1 = iw2n1[iw1]; + const int m1 = iw2m1[iw1]; + for (int iw2l = 0; iw2l < col_indexes.size(); ++iw2l) + { + const int iw2 = col_indexes[iw2l]; + const int L2 = iw2l2[iw2]; + const int N2 = iw2n2[iw2]; + const int m2 = iw2m2[iw2]; + uot.snap_psipsi(orb, // orbitals + olm, 0, 'S', // olm, job of derivation, dtype of Operator + tau1, T1, L1, m1, N1, // info of atom1 + tau2, T2, L2, m2, N2 // info of atom2 + ); + *data_pointer++ = olm[0]; + } + } +} + +// contributeHR() +template +void hamilt::OverlapNew, TR>::contributeHR() +{ + if (this->SR_fixed_done) + { + return; + } + this->calculate_SR(); + this->SR_fixed_done = true; +} + +// contributeHk() +template +void hamilt::OverlapNew, TR>::contributeHk(int ik) +{ + const int ncol = this->SR->get_atom_pair(0).get_paraV()->get_col_size(); + hamilt::folding_HR(*this->SR, this->SK_pointer, this->kvec_d[ik], ncol, 0); +} + +template class hamilt::OverlapNew, double>; +template class hamilt::OverlapNew>, double>; +template class hamilt::OverlapNew>, std::complex>; \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.h new file mode 100644 index 00000000000..13f38dcd69c --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/overlap_new.h @@ -0,0 +1,85 @@ +#ifndef OVERLAPNEW_H +#define OVERLAPNEW_H +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/unitcell.h" +#include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" + +namespace hamilt +{ + +#ifndef __OVERLAPNEWTEMPLATE +#define __OVERLAPNEWTEMPLATE + +/// The OverlapNew class template inherits from class T +/// it is used to calculate the overlap 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 OverlapNew : public T +{ +}; + +#endif + +/// OverlapNew class template specialization for OperatorLCAO base class +/// It is used to calculate the overlap matrix in real space and fold it to k-space +/// SR = +/// SK = = \sum_{R} e^{ikR} SR +/// Template parameters: +/// - TK: data type of k-space Hamiltonian +/// - TR: data type of real space Hamiltonian +template +class OverlapNew, TR> : public OperatorLCAO +{ + public: + OverlapNew, TR>(LCAO_Matrix* LM_in, + const std::vector>& kvec_d_in, + hamilt::HContainer* SR_in, + TK* SK_pointer_in, + const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const Parallel_Orbitals* paraV); + + virtual void contributeHR() override; + + virtual void contributeHk(int ik) override; + + private: + const UnitCell* ucell = nullptr; + + hamilt::HContainer* SR = nullptr; + + TK* SK_pointer = nullptr; + + bool SR_fixed_done = false; + + /** + * @brief initialize SR, search the nearest neighbor atoms + * HContainer is used to store the overlap matrix with specific atom-pairs + * the size of SR will be fixed after initialization + */ + void initialize_SR(Grid_Driver* GridD_in, const Parallel_Orbitals* paraV); + + /** + * @brief calculate the overlap matrix with specific atom-pairs + * nearest neighbor atoms don't need to be calculated again + * loop the atom-pairs in SR and calculate the overlap matrix + */ + void calculate_SR(); + + /** + * @brief calculate the SR local matrix of atom pair + */ + void cal_SR_IJR( + const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const ModuleBase::Vector3& dtau, + TR* data_pointer); +}; + +} // namespace hamilt +#endif \ No newline at end of file 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 new file mode 100644 index 00000000000..602e928f4ac --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt @@ -0,0 +1,12 @@ +if(ENABLE_LCAO) + +AddTest( + TARGET test_overlap_new + LIBS ${math_libs} psi base device + SOURCES test_overlapnew.cpp ../overlap_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_hcontainer/test/tmp_mocks.cpp +) + +endif() \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew.cpp new file mode 100644 index 00000000000..203715a8e48 --- /dev/null +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_overlapnew.cpp @@ -0,0 +1,270 @@ +#include "gtest/gtest.h" +#include "../overlap_new.h" + +// mock of OperatorLCAO + +template +hamilt::Operator::Operator(){} + +template +hamilt::Operator::~Operator(){} + +template +typename hamilt::Operator::hpsi_info hamilt::Operator::hPsi(hpsi_info&) const +{ + return hpsi_info(nullptr, 0, nullptr); +} + +template +void hamilt::Operator::init(const int ik_in) +{ + return; +} + +template +void hamilt::Operator::add(Operator* next) +{ + return; +} + +template +FPTYPE* hamilt::Operator::get_hpsi(const hpsi_info& info) const +{ + return nullptr; +} + +template class hamilt::Operator; +template class hamilt::Operator, psi::DEVICE_CPU>; + +// mock of OperatorLCAO +template +void hamilt::OperatorLCAO::init(const int ik_in) +{ + return; +} +template +void hamilt::OperatorLCAO::get_hs_pointers() +{ + return; +} +template class hamilt::OperatorLCAO; +template class hamilt::OperatorLCAO>; + +// mock of ORB_gen_tables and LCAO_Orbitals +#include "module_basis/module_ao/ORB_gen_tables.h" +const ORB_gen_tables& ORB_gen_tables::get_const_instance() +{ + static ORB_gen_tables instance; + return instance; +} +ORB_gen_tables::ORB_gen_tables() {} +ORB_gen_tables::~ORB_gen_tables() {} +ORB_gaunt_table::ORB_gaunt_table() {} +ORB_gaunt_table::~ORB_gaunt_table() {} +ORB_table_phi::ORB_table_phi() {} +ORB_table_phi::~ORB_table_phi() {} +ORB_table_alpha::ORB_table_alpha() {} +ORB_table_alpha::~ORB_table_alpha() {} +ORB_table_beta::ORB_table_beta() {} +ORB_table_beta::~ORB_table_beta() {} +// mock of snap_psipsi +void ORB_gen_tables::snap_psipsi( + const LCAO_Orbitals &orb, + double olm[], + const int &job, ///<[in]0 for matrix element of either S or T, 1 for its derivatives + const char &dtype, ///<[in] derivative type, 'S' for overlap, 'T' for kinetic energy, 'D' for descriptor in deepks + const ModuleBase::Vector3 &R1, + const int &I1, + const int &l1, + const int &m1, + const int &n1, + const ModuleBase::Vector3 &R2, + const int &I2, + const int &l2, + const int &m2, + const int &n2, + bool cal_syns, + double dmax)const +{ + olm[0] = 1.0; +} + +#include "module_basis/module_ao/ORB_read.h" +const LCAO_Orbitals& LCAO_Orbitals::get_const_instance() +{ + static LCAO_Orbitals instance; + return instance; +} +LCAO_Orbitals::LCAO_Orbitals() {} +LCAO_Orbitals::~LCAO_Orbitals() {} + +// mock find_atom() function +void Grid_Driver::Find_atom(const UnitCell& ucell, + const ModuleBase::Vector3& tau, + const int& T, + const int& I, + AdjacentAtomInfo* adjs) +{ + adjs->adj_num = ucell.nat - 1 ; + adjs->adjacent_tau.resize(ucell.nat); + adjs->ntype.resize(ucell.nat, 0); + adjs->natom.resize(ucell.nat); + adjs->box.resize(ucell.nat); + for(int iat = 0;iatnatom[iat] = iat; + adjs->box[iat].x = 0; + adjs->box[iat].y = 0; + adjs->box[iat].z = 0; + adjs->adjacent_tau[iat] = ucell.get_tau(iat); + } +} +Grid::Grid(const int &test_grid_in):test_grid(test_grid_in) +{} +Grid::~Grid(){} +Grid_Driver::Grid_Driver(const int &test_d_in, + const int &test_gd_in, + const int &test_grid_in) :Grid(test_grid_in), test_deconstructor(test_d_in), test_grid_driver(test_gd_in) {} +Grid_Driver::~Grid_Driver() {} +// simple test for OverlapNew class +// 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; +class OverlapNewTest : 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.itia2iat.create(ucell.ntype, ucell.nat); + ucell.iat2iwt.resize(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.iat2iwt[iat] = iat * test_nw; + } + 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] = 0; + ucell.atoms[0].iw2m[iw] = 0; + ucell.atoms[0].iw2n[iw] = 0; + } + init_parav(); + // set up a HContainer with ucell + SR = new hamilt::HContainer(paraV); + } + + void TearDown() override + { + delete SR; + delete paraV; + delete[] ucell.atoms[0].tau; + delete[] ucell.atoms[0].iw2l; + delete[] ucell.atoms[0].iw2m; + delete[] ucell.atoms[0].iw2n; + delete[] ucell.atoms; + delete[] ucell.iat2it; + delete[] ucell.iat2ia; + } + +#ifdef __MPI + void init_parav() + { + int global_row = test_size * test_nw; + int global_col = test_size * test_nw; + std::ofstream ofs_running; + paraV = new Parallel_Orbitals(); + paraV->set_block_size(2/* 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.iat2iwt.data(), test_size, global_row); + } +#else + void init_parav() + {} +#endif + + UnitCell ucell; + hamilt::HContainer* SR; + Parallel_Orbitals *paraV; + + int dsize; + int my_rank = 0; +}; + +// using TEST_F to test OverlapNew +TEST_F(OverlapNewTest, constructHR) +{ + 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); + hamilt::OverlapNew, double> op( + nullptr, + kvec_d_in, + SR, + hk.data(), + &ucell, + &gd, + paraV + ); + op.contributeHR(); + // check the value of SR + for (int iap = 0; iap < SR->size_atom_pairs(); ++iap) + { + hamilt::AtomPair& tmp = SR->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_EQ(tmp.get_HR_values(0, 0, 0).get_pointer()[i], 1.0); + } + } + // calculate SK + op.contributeHk(0); + // check the value of SK + for (int i = 0; i < paraV->get_row_size() * paraV->get_col_size(); ++i) + { + EXPECT_EQ(hk[i], 1.0); + } +} + +int main(int argc, char** argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); +#ifdef __MPI + MPI_Finalize(); +#endif + return result; +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h index 820f425657a..f06f36b62c0 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/veff_lcao.h @@ -36,9 +36,8 @@ class Veff> : public OperatorLCAO pot(pot_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), - OperatorLCAO(kvec_d_in) + OperatorLCAO(LM_in, kvec_d_in) { - this->LM = LM_in; this->cal_type = lcao_gint; } Veff>(Gint_Gamma* GG_in, @@ -48,9 +47,9 @@ class Veff> : public OperatorLCAO std::vector* HR_pointer_in, std::vector* HK_pointer_in, std::vector> kvec_d_in) - : GG(GG_in), loc(loc_in), pot(pot_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), OperatorLCAO(std::vector>{ModuleBase::Vector3(0,0,0)}) + : GG(GG_in), loc(loc_in), pot(pot_in), HR_pointer(HR_pointer_in), HK_pointer(HK_pointer_in), + OperatorLCAO(LM_in, std::vector>{ModuleBase::Vector3(0,0,0)}) { - this->LM = LM_in; this->cal_type = lcao_gint; } diff --git a/source/module_hamilt_lcao/module_hcontainer/CMakeLists.txt b/source/module_hamilt_lcao/module_hcontainer/CMakeLists.txt index 6eb2cac6ec1..bfbe7282414 100644 --- a/source/module_hamilt_lcao/module_hcontainer/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_hcontainer/CMakeLists.txt @@ -3,6 +3,7 @@ list(APPEND objects atom_pair.cpp hcontainer.cpp output_hcontainer.cpp + func_folding.cpp ) add_library( diff --git a/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp b/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp index ab9f3989db6..d56bf2b2144 100644 --- a/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/atom_pair.cpp @@ -12,6 +12,7 @@ namespace hamilt template AtomPair::~AtomPair() { + this->values.clear(); } template @@ -590,13 +591,17 @@ T& AtomPair::get_value(const int& row, const int& col) const // get_pointer template -T* AtomPair::get_pointer() const +T* AtomPair::get_pointer(int ir) const { + if(ir<0) + { + ir = current_R; + } #ifdef __DEBUG assert(current_R < this->values.size()); assert(current_R >= 0); #endif - return this->values[current_R].get_pointer(); + return this->values[ir].get_pointer(); } // get_R_size diff --git a/source/module_hamilt_lcao/module_hcontainer/atom_pair.h b/source/module_hamilt_lcao/module_hcontainer/atom_pair.h index 2ac67f1da86..303d37d620a 100644 --- a/source/module_hamilt_lcao/module_hcontainer/atom_pair.h +++ b/source/module_hamilt_lcao/module_hcontainer/atom_pair.h @@ -148,7 +148,7 @@ class AtomPair T& get_value(const int& i) const; T& get_value(const int& row, const int& col) const; T& get_matrix_value(const size_t& i_row_global, const size_t& j_col_global) const; - T* get_pointer() const; + T* get_pointer(int ir=-1) const; // add another BaseMatrix to this->values with specific R index. void convert_add(const BaseMatrix& target, int rx_in, int ry_in, int rz_in); diff --git a/source/module_hamilt_lcao/module_hcontainer/func_folding.cpp b/source/module_hamilt_lcao/module_hcontainer/func_folding.cpp new file mode 100644 index 00000000000..3ba83225173 --- /dev/null +++ b/source/module_hamilt_lcao/module_hcontainer/func_folding.cpp @@ -0,0 +1,91 @@ +#include "hcontainer_funcs.h" +#include "module_base/libm/libm.h" + +namespace hamilt +{ +/** + * @brief calculate the Hk matrix with specific k vector + * @param hR the HContainer of atom pairs + * @param hk the data pointer of Hk matrix, the size of hk would be nrow * ncol + * @param kvec_d_in the k vector in Direct coordinate + * @param hk_ld the leading dimension number of hk, ncol for row-major, nrow for column-major + * @param hk_type the data-type of hk, 0 is row-major, 1 is column-major +*/ +template +void folding_HR(const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ncol, + const int hk_type) +{ + for (int i = 0; i < hR.size_R_loop(); ++i) + { + // get R index + int rx, ry, rz; + hR.loop_R(i, rx, ry, rz); + // only deal with current_R for hR + hR.fix_R(rx, ry, rz); + + // cal k_phase + // if TK==std::complex, kphase is e^{ikR} + const ModuleBase::Vector3 dR(rx, ry, rz); + const double arg = (kvec_d_in * dR) * ModuleBase::TWO_PI; + double sinp, cosp; + ModuleBase::libm::sincos(arg, &sinp, &cosp); + std::complex kphase = std::complex(cosp, sinp); + + // loop_atom_pairs +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int j = 0; j < hR.size_atom_pairs(); ++j) + { + // Hk += HR * e^ikR + hR.get_atom_pair(j).add_to_matrix(hk, ncol, kphase, hk_type); + } + } +} + +// template instantiation +template void folding_HR>(const hamilt::HContainer>& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ncol, + const int hk_type); +template void folding_HR(const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ncol, + const int hk_type); +// special case for double +void folding_HR(const hamilt::HContainer& hR, + double* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ncol, + const int hk_type) +{ + for (int i = 0; i < hR.size_R_loop(); ++i) + { + // get R index + int rx, ry, rz; + hR.loop_R(i, rx, ry, rz); + // only deal with current_R for hR + hR.fix_R(rx, ry, rz); + + // cal k_phase + // if TK==double, kphase is 1.0 + double kphase = 1.0; + + // loop_atom_pairs +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int j = 0; j < hR.size_atom_pairs(); ++j) + { + // Hk += HR * e^ikR + hR.get_atom_pair(j).add_to_matrix(hk, ncol, kphase, hk_type); + } + } +} + +} // namespace hamilt \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp b/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp index 869f0e204ab..c682873fca3 100644 --- a/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer.cpp @@ -53,6 +53,16 @@ HContainer::HContainer(const UnitCell& ucell_) { this->gamma_only = false; this->current_R = -1; + std::vector atom_begin_row(ucell_.nat+1, 0); + std::vector atom_begin_col(ucell_.nat+1, 0); + int begin = 0; + for(int i=0;iatom_pairs.reserve(ucell_.nat * ucell_.nat); this->sparse_ap.resize(ucell_.nat); @@ -63,10 +73,7 @@ HContainer::HContainer(const UnitCell& ucell_) this->sparse_ap_index[i].resize(ucell_.nat); for (int j = 0; j < ucell_.nat; j++) { - AtomPair atom_ij(i, j); - int it1 = ucell_.iat2it[i]; - int it2 = ucell_.iat2it[j]; - atom_ij.set_size(ucell_.atoms[it1].nw, ucell_.atoms[it2].nw); + AtomPair atom_ij(i, j, atom_begin_row.data(), atom_begin_col.data(), ucell_.nat); ModuleBase::GlobalFunc::ZEROS(atom_ij.get_HR_values(0, 0, 0).get_pointer(), atom_ij.get_size()); this->atom_pairs.push_back(atom_ij); this->sparse_ap[i][j] = j; diff --git a/source/module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h b/source/module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h new file mode 100644 index 00000000000..43d83934798 --- /dev/null +++ b/source/module_hamilt_lcao/module_hcontainer/hcontainer_funcs.h @@ -0,0 +1,28 @@ +#pragma once + +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" + +namespace hamilt +{ +/** + * @brief calculate the Hk matrix with specific k vector + * @param hR the HContainer of atom pairs + * @param hk the data pointer of Hk matrix, the size of hk would be nrow * ncol + * @param kvec_d_in the k vector in Direct coordinate + * @param hk_ld the leading dimension number of hk, ncol for row-major, nrow for column-major + * @param hk_type the data-type of hk, 0 is row-major, 1 is column-major +*/ +template +void folding_HR(const hamilt::HContainer& hR, + std::complex* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ncol, + const int hk_type); + +void folding_HR(const hamilt::HContainer& hR, + double* hk, + const ModuleBase::Vector3& kvec_d_in, + const int ncol, + const int hk_type); + +} // namespace hamilt \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_hcontainer/test/CMakeLists.txt b/source/module_hamilt_lcao/module_hcontainer/test/CMakeLists.txt index 4b32fca0da2..4918a68f9d6 100644 --- a/source/module_hamilt_lcao/module_hcontainer/test/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_hcontainer/test/CMakeLists.txt @@ -21,6 +21,13 @@ AddTest( ../../../module_basis/module_ao/parallel_2d.cpp ../../../module_basis/module_ao/parallel_orbitals.cpp tmp_mocks.cpp ) +AddTest( + TARGET test_func_folding + LIBS ${math_libs} psi base device + SOURCES test_func_folding.cpp ../base_matrix.cpp ../hcontainer.cpp ../atom_pair.cpp + ../func_folding.cpp ../../../module_basis/module_ao/parallel_2d.cpp ../../../module_basis/module_ao/parallel_orbitals.cpp tmp_mocks.cpp +) + AddTest( TARGET test_hcontainer_output LIBS base ${math_libs} device diff --git a/source/module_hamilt_lcao/module_hcontainer/test/demo1_SR.cpp b/source/module_hamilt_lcao/module_hcontainer/test/demo1_SR.cpp deleted file mode 100644 index 403e524979b..00000000000 --- a/source/module_hamilt_lcao/module_hcontainer/test/demo1_SR.cpp +++ /dev/null @@ -1,42 +0,0 @@ -#include "../hcontainer.h" -void construct_SR(const Parallel_Orbitals* paraV, const UnitCell& ucell) -{ - this->SR = new hamilt::HContainer(paraV); -#ifdef _OPENMP -#pragma omp parallel for -#endif - for(int iat1=0;iat1 SR_thread(paraV); -#endif - for(int ad = 0; ad < GlobalC::GridD.getAdjacentNum()+1; ++ad) - { - int iat2 = GlobalC::GridD.getAdjacent(ad); - auto tau2 = ucell.get_tau(iat2); - int T2 = ucell.get_element(iat2); - int I2 = ucell.get_element_atom(iat2); - auto R_index = GlobalC::GridD.getBox(ad); - double olm[3] = {0,0,0}; - hamilt::AtomPair tmp(iat1, iat2, paraV); - double* data_pointer = tmp.get_HR_values(R_index.x, R_index_y, R_index_z).get_pointer(); - for(iw1=0;iw2insert_pair(tmp); -#else - this->SR->insert_pair(tmp); -#endif - } -#ifdef _OPENMP -#pragma omp critical(append_SR) - this->SR->add_HR(SR_thread); -#endif - } -} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_hcontainer/test/demo4_ft.cpp b/source/module_hamilt_lcao/module_hcontainer/test/demo4_ft.cpp deleted file mode 100644 index d5a1ff03e6e..00000000000 --- a/source/module_hamilt_lcao/module_hcontainer/test/demo4_ft.cpp +++ /dev/null @@ -1,33 +0,0 @@ -#include "../hcontainer.h" - -void folding_HR(const hamilt::HContainer& hR, - std::complex* hk, - const ModuleBase::Vector3& kvec_d_in, - const int ncol) -{ - for (int i = 0; i < hR.size_R_loop(); ++i) - { - // get R index - int rx, ry, rz; - hR.loop_R(i, rx, ry, rz); - // only deal with current_R for hR - hR.fix_R(rx, ry, rz); - - // cal k_phase - ModuleBase::Vector3 dR(rx, ry, rz); - const double arg = (kvec_d_in * dR) * ModuleBase::TWO_PI; - double sinp, cosp; - ModuleBase::libm::sincos(arg, &sinp, &cosp); - const std::complex kphase = std::complex(cosp, sinp); - - // loop_atom_pairs -#ifdef _OPENMP -#pragma omp parallel for -#endif - for (int j = 0; j < hR.size_atom_pair(); ++j) - { - // Hk += HR * e^ikR - hR.get_atom_pair(j).add_to_matrix(hk, ncol, kphase, 0); - } - } -} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_func_folding.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_func_folding.cpp new file mode 100644 index 00000000000..fde0b74cc66 --- /dev/null +++ b/source/module_hamilt_lcao/module_hcontainer/test/test_func_folding.cpp @@ -0,0 +1,70 @@ +#include "gtest/gtest.h" +#include "../hcontainer_funcs.h" +#include "../hcontainer.h" + +// 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; + +/** + * Unit test of folding_HR function + * folding_HR is a function to calculate the Hk matrix with specific k vector + * in this test, we test the correctness and time consuming of the function. +*/ +class FoldingTest : public ::testing::Test +{ + protected: + void SetUp() override + { + // set up a unitcell, with one element and three atoms, each atom has 2 orbitals + ucell.ntype = 1; + ucell.nat = test_size; + ucell.atoms = new Atom[ucell.ntype]; + ucell.iat2it = new int[ucell.nat]; + for (int iat = 0; iat < ucell.nat; iat++) + { + ucell.iat2it[iat] = 0; + } + ucell.atoms[0].nw = test_nw; + + // set up a HContainer with ucell + HR = new hamilt::HContainer>(ucell); + } + + void TearDown() override + { + delete HR; + delete[] ucell.atoms; + delete[] ucell.iat2it; + } + + UnitCell ucell; + hamilt::HContainer>* HR; +}; + +// using TEST_F to test folding_HR +TEST_F(FoldingTest, folding_HR_cd2cd) +{ + // fill HR with constant value + for (int i = 0; i < HR->size_atom_pairs(); i++) + { + std::complex* ptr1 = HR->get_atom_pair(i).get_HR_values(1, 1, 1).get_pointer(); + std::complex* ptr = HR->get_atom_pair(i).get_HR_values(0,0,0).get_pointer(); + for (int j = 0; j < HR->get_atom_pair(i).get_size(); j++) + { + ptr[j] = std::complex(1.0, 1.0); + ptr1[j] = std::complex(2.0, 2.0); + } + } + std::vector> hk(test_size * test_nw * test_size * test_nw, std::complex(0.0, 0.0)); + ModuleBase::Vector3 kvec_d_in(0.1, 0.2, 0.3); + hamilt::folding_HR(*HR, hk.data(), kvec_d_in, test_size * test_nw, 0); + for (int i = 0; i < test_size * test_nw * test_size * test_nw; i++) + { + EXPECT_NEAR(hk[i].real(), 0.55753651583505226, 1e-10); + EXPECT_NEAR(hk[i].imag(), -1.7936044933348412, 1e-10); + //std::cout<<__FILE__<<__LINE__<<" hk["<data(0, 1); // check if data pointer is correct EXPECT_EQ(data_ptr, HR->get_atom_pair(0, 1).get_pointer()); + EXPECT_EQ(data_ptr, HR->get_atom_pair(0, 1).get_pointer(0)); int r_index[3] = {0, 0, 0}; EXPECT_EQ(data_ptr, HR->data(0, 1, r_index)); EXPECT_EQ(HR->data(0, 10), nullptr); diff --git a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp index c3d6408d37b..fc5cc47ddb1 100644 --- a/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp +++ b/source/module_hamilt_lcao/module_hcontainer/test/test_hcontainer_complex.cpp @@ -356,6 +356,7 @@ TEST_F(HContainerTest, data) std::complex* data_ptr = HR->data(0, 1); // check if data pointer is correct EXPECT_EQ(data_ptr, HR->get_atom_pair(0, 1).get_pointer()); + EXPECT_EQ(data_ptr, HR->get_atom_pair(0, 1).get_pointer(0)); int r_index[3] = {0, 0, 0}; EXPECT_EQ(data_ptr, HR->data(0, 1, r_index)); EXPECT_EQ(HR->data(0, 10), nullptr);