diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 856e56d97a..69d08d9b36 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -2,27 +2,37 @@ name: Integration Test and Unit Test on: pull_request: - + concurrency: group: ${{ github.workflow }}-${{ github.ref }} cancel-in-progress: true - + jobs: test: name: Test runs-on: self-hosted if: github.repository_owner == 'deepmodeling' - container: ghcr.io/deepmodeling/abacus-gnu + container: + image: ghcr.io/deepmodeling/abacus-gnu + volumes: + - /tmp/ccache:/github/home/.ccache steps: - name: Checkout uses: actions/checkout@v4 with: submodules: recursive + + - name: Install Ccache + run: | + sudo apt-get update + sudo apt-get install -y ccache + - name: Build run: | cmake -B build -DBUILD_TESTING=ON -DENABLE_DEEPKS=ON -DENABLE_LIBXC=ON -DENABLE_LIBRI=ON -DENABLE_PAW=ON cmake --build build -j8 cmake --install build + - name: Test env: GTEST_COLOR: 'yes' diff --git a/source/module_base/math_sphbes.cpp b/source/module_base/math_sphbes.cpp index 5e7f41de54..00326f928e 100644 --- a/source/module_base/math_sphbes.cpp +++ b/source/module_base/math_sphbes.cpp @@ -808,7 +808,7 @@ void Sphbes::dsphbesj(const int n, } } -void Sphbes::sphbes_zeros(const int l, const int n, double* const zeros) +void Sphbes::sphbes_zeros(const int l, const int n, double* const zeros, const bool return_all) { assert( n > 0 ); assert( l >= 0 ); @@ -818,10 +818,22 @@ void Sphbes::sphbes_zeros(const int l, const int n, double* const zeros) // This property enables us to use bracketing method recursively // to find all zeros of j_l from the zeros of j_0. - // if l is odd , j_0 --> j_1 --> j_3 --> j_5 --> ... - // if l is even, j_0 --> j_2 --> j_4 --> j_6 --> ... - - int nz = n + (l+1)/2; // number of effective zeros in buffer + // If return_all is true, zeros of j_0, j_1, ..., j_l will all be returned + // such that zeros[l*n+i] is the i-th zero of j_l. As such, it is required + // that the array "zeros" has a size of (l+1)*n. + // + // If return_all is false, only the zeros of j_l will be returned + // and "zeros" is merely required to have a size of n. + // Note that in this case the bracketing method can be applied with a stride + // of 2 instead of 1: + // j_0 --> j_1 --> j_3 --> j_5 --> ... --> j_l (odd l) + // j_0 --> j_2 --> j_4 --> j_6 --> ... --> j_l (even l) + + // Every recursion step reduces the number of zeros by 1. + // If return_all is true, one needs to start with n+l zeros of j_0 + // to ensure n zeros of j_l; otherwise with a stride of 2 one only + // needs to start with n+(l+1)/2 zeros of j_0 + int nz = n + ( return_all ? l : (l+1)/2 ); double* buffer = new double[nz]; // zeros of j_0 = sin(x)/x is just n*pi @@ -831,27 +843,34 @@ void Sphbes::sphbes_zeros(const int l, const int n, double* const zeros) buffer[i] = (i+1) * PI; } - int ll = 1; + int ll; // active l auto jl = [&ll] (double x) { return sphbesj(ll, x); }; - - if (l % 2 == 1) + int stride; + std::function copy_if_needed; + int offset = 0; // keeps track of the position in zeros for next copy (used when return_all == true) + if (return_all) { - for (int i = 0; i < nz-1; i++) - { - buffer[i] = illinois(jl, buffer[i], buffer[i+1], 1e-15, 50); - } - --nz; + copy_if_needed = [&](){ std::copy(buffer, buffer + n, zeros + offset); offset += n; }; + stride = 1; + ll = 1; + } + else + { + copy_if_needed = [](){}; + stride = 2; + ll = 2 - l % 2; } - for (ll = 2 + l%2; ll <= l; ll += 2, --nz) + for (; ll <= l; ll += stride, --nz) { + copy_if_needed(); for (int i = 0; i < nz-1; i++) { buffer[i] = illinois(jl, buffer[i], buffer[i+1], 1e-15, 50); } } - std::copy(buffer, buffer + n, zeros); + std::copy(buffer, buffer + n, zeros + offset); delete[] buffer; } diff --git a/source/module_base/math_sphbes.h b/source/module_base/math_sphbes.h index c654847a5d..7aa9c78a48 100644 --- a/source/module_base/math_sphbes.h +++ b/source/module_base/math_sphbes.h @@ -126,13 +126,18 @@ class Sphbes * This function computes the first n positive zeros of the l-th order * spherical Bessel function of the first kind. * - * @param[in] l order of the spherical Bessel function - * @param[in] n number of zeros to be computed - * @param[out] zeros on exit, contains the first n positive zeros in ascending order + * @param[in] l (maximum) order of the spherical Bessel function + * @param[in] n number of zeros to be computed (for each j_l if return_all is true) + * @param[out] zeros on exit, contains the positive zeros. + * @param[in] return_all if true, return all zeros from j_0 to j_l such that zeros[l*n+i] + * is the i-th zero of j_l. If false, return only the first n zeros of j_l. + * + * @note The size of array "zeros" must be at least (l+1)*n if return_all is true, and n otherwise. */ static void sphbes_zeros(const int l, const int n, - double* const zeros + double* const zeros, + bool return_all = false ); private: diff --git a/source/module_base/test/math_sphbes_test.cpp b/source/module_base/test/math_sphbes_test.cpp index 521d4dc2f4..e72c6e289c 100644 --- a/source/module_base/test/math_sphbes_test.cpp +++ b/source/module_base/test/math_sphbes_test.cpp @@ -352,15 +352,27 @@ TEST_F(Sphbes, Zeros) int lmax = 20; int nzeros = 500; - double* zeros = new double[nzeros]; + double* zeros = new double[nzeros*(lmax+1)]; for (int l = 0; l <= lmax; ++l) { - ModuleBase::Sphbes::sphbes_zeros(l, nzeros, zeros); + ModuleBase::Sphbes::sphbes_zeros(l, nzeros, zeros, false); for (int i = 0; i < nzeros; ++i) { EXPECT_LT(std::abs(ModuleBase::Sphbes::sphbesj(l, zeros[i])), 1e-14); } } + + + ModuleBase::Sphbes::sphbes_zeros(lmax, nzeros, zeros, true); + for (int l = 0; l <= lmax; ++l) + { + for (int i = 0; i < nzeros; ++i) + { + EXPECT_LT(std::abs(ModuleBase::Sphbes::sphbesj(l, zeros[l*nzeros+i])), 1e-14); + } + } + + delete[] zeros; } TEST_F(Sphbes, ZerosOld) diff --git a/source/module_cell/read_atoms.cpp b/source/module_cell/read_atoms.cpp index dc517bccd7..4c6bf9c0eb 100644 --- a/source/module_cell/read_atoms.cpp +++ b/source/module_cell/read_atoms.cpp @@ -535,100 +535,101 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn ModuleBase::GlobalFunc::ZEROS(atoms[it].mag,na); for (int ia = 0;ia < na; ia++) { - // modify the reading of frozen ions and velocities -- Yuanbo Li 2021/8/20 - ifpos >> v.x >> v.y >> v.z; - mv.x = true ; - mv.y = true ; - mv.z = true ; - atoms[it].vel[ia].set(0,0,0); - atoms[it].mag[ia]=magnet.start_magnetization[it];//if this line is used, default startmag_type would be 2 - atoms[it].angle1[ia]=0; - atoms[it].angle2[ia]=0; - atoms[it].m_loc_[ia].set(0,0,0); - - std::string tmpid; - tmpid = ifpos.get(); - - if( (int)tmpid[0] < 0 ) - { - std::cout << "read_atom_positions, mismatch in atom number for atom type: " << atoms[it].label << std::endl; - exit(1); - } - - bool input_vec_mag=false; - bool input_angle_mag=false; - while ( (tmpid != "\n") && (ifpos.eof()==false) && (tmpid !="#") ) - { - tmpid = ifpos.get() ; - // old method of reading frozen ions - char tmp = (char)tmpid[0]; - if ( tmp >= 48 && tmp <= 57 ) - { - mv.x = std::stoi(tmpid); - ifpos >> mv.y >> mv.z ; - } - // new method of reading frozen ions and velocities - if ( tmp >= 'a' && tmp <='z') - { - ifpos.putback(tmp); - ifpos >> tmpid; - } - if ( tmpid == "m" ) - { - ifpos >> mv.x >> mv.y >> mv.z ; - } - else if ( tmpid == "v" ||tmpid == "vel" || tmpid == "velocity" ) - { - ifpos >> atoms[it].vel[ia].x >> atoms[it].vel[ia].y >> atoms[it].vel[ia].z; - } - else if ( tmpid == "mag" || tmpid == "magmom") - { - set_element_mag_zero = true; - double tmpamg=0; - ifpos >> tmpamg; - tmp=ifpos.get(); - while (tmp==' ') - { - tmp=ifpos.get(); - } - - if((tmp >= 48 && tmp <= 57) or tmp=='-') - { - ifpos.putback(tmp); - ifpos >> atoms[it].m_loc_[ia].y>>atoms[it].m_loc_[ia].z; - atoms[it].m_loc_[ia].x=tmpamg; - atoms[it].mag[ia]=sqrt(pow(atoms[it].m_loc_[ia].x,2)+pow(atoms[it].m_loc_[ia].y,2)+pow(atoms[it].m_loc_[ia].z,2)); - input_vec_mag=true; - - } - else - { - ifpos.putback(tmp); - atoms[it].mag[ia]=tmpamg; - } - - // atoms[it].mag[ia]; - } - else if ( tmpid == "angle1") - { - ifpos >> atoms[it].angle1[ia]; - atoms[it].angle1[ia]=atoms[it].angle1[ia]/180 *ModuleBase::PI; - input_angle_mag=true; - set_element_mag_zero = true; - } - else if ( tmpid == "angle2") - { - ifpos >> atoms[it].angle2[ia]; - atoms[it].angle2[ia]=atoms[it].angle2[ia]/180 *ModuleBase::PI; - input_angle_mag=true; - set_element_mag_zero = true; - } - - } - while ( (tmpid != "\n") && (ifpos.eof()==false) ) - { - tmpid = ifpos.get(); - } + // modify the reading of frozen ions and velocities -- Yuanbo Li 2021/8/20 + ifpos >> v.x >> v.y >> v.z; + mv.x = true ; + mv.y = true ; + mv.z = true ; + atoms[it].vel[ia].set(0,0,0); + atoms[it].mag[ia]=magnet.start_magnetization[it];//if this line is used, default startmag_type would be 2 + atoms[it].angle1[ia]=0; + atoms[it].angle2[ia]=0; + atoms[it].m_loc_[ia].set(0,0,0); + + std::string tmpid; + tmpid = ifpos.get(); + + if( (int)tmpid[0] < 0 ) + { + std::cout << "read_atom_positions, mismatch in atom number for atom type: " << atoms[it].label << std::endl; + exit(1); + } + + bool input_vec_mag=false; + bool input_angle_mag=false; + // read if catch goodbit before "\n" and "#" + while ( (tmpid != "\n") && (ifpos.good()) && (tmpid !="#") ) + { + tmpid = ifpos.get() ; + // old method of reading frozen ions + char tmp = (char)tmpid[0]; + if ( tmp >= 48 && tmp <= 57 ) + { + mv.x = std::stoi(tmpid); + ifpos >> mv.y >> mv.z ; + } + // new method of reading frozen ions and velocities + if ( tmp >= 'a' && tmp <='z') + { + ifpos.putback(tmp); + ifpos >> tmpid; + } + if ( tmpid == "m" ) + { + ifpos >> mv.x >> mv.y >> mv.z ; + } + else if ( tmpid == "v" ||tmpid == "vel" || tmpid == "velocity" ) + { + ifpos >> atoms[it].vel[ia].x >> atoms[it].vel[ia].y >> atoms[it].vel[ia].z; + } + else if ( tmpid == "mag" || tmpid == "magmom") + { + set_element_mag_zero = true; + double tmpamg=0; + ifpos >> tmpamg; + tmp=ifpos.get(); + while (tmp==' ') + { + tmp=ifpos.get(); + } + + if((tmp >= 48 && tmp <= 57) or tmp=='-') + { + ifpos.putback(tmp); + ifpos >> atoms[it].m_loc_[ia].y>>atoms[it].m_loc_[ia].z; + atoms[it].m_loc_[ia].x=tmpamg; + atoms[it].mag[ia]=sqrt(pow(atoms[it].m_loc_[ia].x,2)+pow(atoms[it].m_loc_[ia].y,2)+pow(atoms[it].m_loc_[ia].z,2)); + input_vec_mag=true; + + } + else + { + ifpos.putback(tmp); + atoms[it].mag[ia]=tmpamg; + } + + // atoms[it].mag[ia]; + } + else if ( tmpid == "angle1") + { + ifpos >> atoms[it].angle1[ia]; + atoms[it].angle1[ia]=atoms[it].angle1[ia]/180 *ModuleBase::PI; + input_angle_mag=true; + set_element_mag_zero = true; + } + else if ( tmpid == "angle2") + { + ifpos >> atoms[it].angle2[ia]; + atoms[it].angle2[ia]=atoms[it].angle2[ia]/180 *ModuleBase::PI; + input_angle_mag=true; + set_element_mag_zero = true; + } + } + // move to next line + while ( (tmpid != "\n") && (ifpos.good()) ) + { + tmpid = ifpos.get(); + } std::string mags; //cout<<"mag"<, psi::DEVICE_CPU>::cal_MW_ const char N_char = 'N'; const int one_int = 1; const std::complex one_float = {1.0, 0.0}, zero_float = {0.0, 0.0}; - pzgemm_(&T_char, + pzgemm_(&N_char, &T_char, &nw, &nw, diff --git a/source/module_io/mulliken_charge.cpp b/source/module_io/mulliken_charge.cpp index 393da5fda4..bdcdb5a035 100644 --- a/source/module_io/mulliken_charge.cpp +++ b/source/module_io/mulliken_charge.cpp @@ -44,7 +44,7 @@ ModuleBase::matrix ModuleIO::cal_mulliken(const std::vector> const char N_char = 'N'; const int one_int = 1; const double one_float = 1.0, zero_float = 0.0; - pdgemm_(&T_char, + pdgemm_(&N_char, &T_char, &GlobalV::NLOCAL, &GlobalV::NLOCAL, @@ -156,7 +156,7 @@ ModuleBase::matrix ModuleIO::cal_mulliken(const std::vector one_float = {1.0, 0.0}, zero_float = {0.0, 0.0}; - pzgemm_(&T_char, + pzgemm_(&N_char, &T_char, &GlobalV::NLOCAL, &GlobalV::NLOCAL, diff --git a/source/module_io/test/to_qo_test.cpp b/source/module_io/test/to_qo_test.cpp index 93692f858e..9477b2eb54 100644 --- a/source/module_io/test/to_qo_test.cpp +++ b/source/module_io/test/to_qo_test.cpp @@ -543,7 +543,39 @@ TEST_F(toQOTest, CalculateSelfOvlpRFull) //tqo.write_ovlp(tqo.ovlp_R()[0], "QO_self_ovlp.dat"); } -TEST_F(toQOTest, BuildPswfc) +/* Si_dojo_soc.upf is special: two p orbitals, one s orbital */ + +TEST_F(toQOTest, BuildPswfcPartial1) +{ + define_fcc_cell(ucell); + toQO tqo("pswfc", {"s", "s"}); + tqo.unwrap_unitcell(&ucell); + tqo.build_ao(ucell.ntype, ucell.pseudo_fn); + EXPECT_EQ(tqo.p_ao()->nchi(), 5); // AO will always read and import all orbitals + EXPECT_EQ(tqo.nchi(), 2); +} + +TEST_F(toQOTest, BuildPswfcPartial2) +{ + define_fcc_cell(ucell); + toQO tqo("pswfc", {"ps", "s"}); + tqo.unwrap_unitcell(&ucell); + tqo.build_ao(ucell.ntype, ucell.pseudo_fn); + EXPECT_EQ(tqo.p_ao()->nchi(), 5); // AO will always read and import all orbitals + EXPECT_EQ(tqo.nchi(), 8); // the first element is Si, it has two p orbitals, so 3+3+1+1 +} + +TEST_F(toQOTest, BuildPswfcPartial3) +{ + define_fcc_cell(ucell); + toQO tqo("pswfc", {"all", "p"}); + tqo.unwrap_unitcell(&ucell); + tqo.build_ao(ucell.ntype, ucell.pseudo_fn); + EXPECT_EQ(tqo.p_ao()->nchi(), 5); // AO will always read and import all orbitals + EXPECT_EQ(tqo.nchi(), 10); +} + +TEST_F(toQOTest, BuildPswfcAll) { define_fcc_cell(ucell); toQO tqo("pswfc", {"all", "all"}); diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp index ace9097bb0..c9b3b69601 100644 --- a/source/module_ri/Exx_LRI.hpp +++ b/source/module_ri/Exx_LRI.hpp @@ -12,7 +12,6 @@ #include "module_ri/exx_abfs-construct_orbs.h" #include "module_ri/exx_abfs-io.h" #include "module_ri/conv_coulomb_pot_k.h" -#include "module_ri/conv_coulomb_pot_k-template.h" #include "module_base/tool_title.h" #include "module_base/timer.h" #include "module_ri/serialization_cereal.h" @@ -71,14 +70,19 @@ void Exx_LRI::init(const MPI_Comm &mpi_comm_in, const K_Vectors &kv_in) case Conv_Coulomb_Pot_K::Ccp_Type::Ccp: return {}; case Conv_Coulomb_Pot_K::Ccp_Type::Hf: - return {}; + { + // 4/3 * pi * Rcut^3 = V_{supercell} = V_{unitcell} * Nk + const int nspin0 = (GlobalV::NSPIN==2) ? 2 : 1; + const double hf_Rcut = std::pow(0.75 * this->p_kv->nkstot_full/nspin0 * GlobalC::ucell.omega / (ModuleBase::PI), 1.0/3.0); + return {{"hf_Rcut", hf_Rcut}}; + } case Conv_Coulomb_Pot_K::Ccp_Type::Hse: return {{"hse_omega", this->info.hse_omega}}; default: throw std::domain_error(std::string(__FILE__)+" line "+std::to_string(__LINE__)); break; } }; - this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times, this->p_kv->nkstot_full); + this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, this->info.ccp_type, get_ccp_parameter(), this->info.ccp_rmesh_times); for( size_t T=0; T!=this->abfs.size(); ++T ) diff --git a/source/module_ri/conv_coulomb_pot_k-template.h b/source/module_ri/conv_coulomb_pot_k-template.h deleted file mode 100644 index 9a3d245286..0000000000 --- a/source/module_ri/conv_coulomb_pot_k-template.h +++ /dev/null @@ -1,51 +0,0 @@ -#ifndef CONV_COULOMB_POT_K_TEMPLATE_H -#define CONV_COULOMB_POT_K_TEMPLATE_H - -#include "conv_coulomb_pot_k.h" -#include -#include - -#include "../module_ri/test_code/exx_abfs-construct_orbs-test.h" - - -template< typename T > -T Conv_Coulomb_Pot_K::cal_orbs_ccp( - const T & orbs, - const Ccp_Type &ccp_type, - const std::map ¶meter, - const double rmesh_times, - const int& nks) -{ - T orbs_ccp(orbs.size()); - for( size_t i=0; i!=orbs.size(); ++i ) - orbs_ccp[i] = cal_orbs_ccp(orbs[i], ccp_type, parameter, rmesh_times, nks ); - return orbs_ccp; -} - -extern template -Numerical_Orbital_Lm Conv_Coulomb_Pot_K::cal_orbs_ccp( - const Numerical_Orbital_Lm & orbs, - const Ccp_Type &ccp_type, - const std::map ¶meter, - const double rmesh_times, - const int& nks); - - - -template< typename T > -double Conv_Coulomb_Pot_K::get_rmesh_proportion( - const T & orbs, - const double psi_threshold) -{ - double rmesh_proportion=0; - for( const auto &orb : orbs ) - rmesh_proportion = std::max(rmesh_proportion, get_rmesh_proportion(orb,psi_threshold)); - return rmesh_proportion; -} - -extern template -double Conv_Coulomb_Pot_K::get_rmesh_proportion( - const Numerical_Orbital_Lm & orbs, - const double psi_threshold); - -#endif \ No newline at end of file diff --git a/source/module_ri/conv_coulomb_pot_k.cpp b/source/module_ri/conv_coulomb_pot_k.cpp index 9f573509ee..62dd582a44 100644 --- a/source/module_ri/conv_coulomb_pot_k.cpp +++ b/source/module_ri/conv_coulomb_pot_k.cpp @@ -2,104 +2,109 @@ #include "../module_base/constants.h" #include "../module_basis/module_ao/ORB_atomic_lm.h" #include "../module_hamilt_pw/hamilt_pwdft/global.h" -std::vector Conv_Coulomb_Pot_K::cal_psi_ccp( const std::vector & psif ) + +namespace Conv_Coulomb_Pot_K { - std::vector psik2_ccp(psif.size()); - for( size_t ik=0; ik Conv_Coulomb_Pot_K::cal_psi_hf(const int& nks, const std::vector &psif, - const std::vector &k_radial, - const double omega = 0) -{ - const int nspin0 = (GlobalV::NSPIN==2) ? 2 : 1; - const double Rc = std::pow(0.75 * nks/nspin0 * GlobalC::ucell.omega / (ModuleBase::PI), 1.0/3.0); - std::vector psik2_ccp(psif.size()); - for (size_t ik = 0; ik < psif.size(); ++ik) - psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1 - std::cos(k_radial[ik] * Rc)); - return psik2_ccp; -} + std::vector cal_psi_ccp( + const std::vector & psif) + { + std::vector psik2_ccp(psif.size()); + for( size_t ik=0; ik cal_psi_hf( + const std::vector &psif, + const std::vector &k_radial, + const double hf_Rcut) + { + std::vector psik2_ccp(psif.size()); + for (size_t ik = 0; ik < psif.size(); ++ik) + psik2_ccp[ik] = ModuleBase::FOUR_PI * psif[ik] * (1 - std::cos(k_radial[ik] * hf_Rcut)); + return psik2_ccp; + } -std::vector Conv_Coulomb_Pot_K::cal_psi_hse( - const std::vector & psif, - const std::vector & k_radial, - const double omega) -{ - std::vector psik2_ccp(psif.size()); - for( size_t ik=0; ik cal_psi_hse( + const std::vector & psif, + const std::vector & k_radial, + const double hse_omega) + { + std::vector psik2_ccp(psif.size()); + for( size_t ik=0; ik -Numerical_Orbital_Lm Conv_Coulomb_Pot_K::cal_orbs_ccp( - const Numerical_Orbital_Lm &orbs, - const Ccp_Type &ccp_type, - const std::map ¶meter, - const double rmesh_times, - const int& nks) -{ - std::vector psik2_ccp; - switch(ccp_type) + + template<> + Numerical_Orbital_Lm cal_orbs_ccp( + const Numerical_Orbital_Lm &orbs, + const Ccp_Type &ccp_type, + const std::map ¶meter, + const double rmesh_times) { - case Ccp_Type::Ccp: - psik2_ccp = cal_psi_ccp( orbs.get_psif() ); break; - case Ccp_Type::Hf: - psik2_ccp = cal_psi_hf(nks, orbs.get_psif(), orbs.get_k_radial()); break; - case Ccp_Type::Hse: - psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") ); break; - default: - throw( ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) ); break; - } + std::vector psik2_ccp; + switch(ccp_type) + { + case Ccp_Type::Ccp: + psik2_ccp = cal_psi_ccp( orbs.get_psif() ); break; + case Ccp_Type::Hf: + psik2_ccp = cal_psi_hf( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hf_Rcut")); break; + case Ccp_Type::Hse: + psik2_ccp = cal_psi_hse( orbs.get_psif(), orbs.get_k_radial(), parameter.at("hse_omega") ); break; + default: + throw( ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) ); break; + } - const double dr = orbs.get_rab().back(); - const int Nr = (static_cast(orbs.getNr()*rmesh_times)) | 1; - std::vector rab(Nr); - for( size_t ir=0; ir r_radial(Nr); - for( size_t ir=0; ir(orbs.getNr()*rmesh_times)) | 1; + std::vector rab(Nr); + for( size_t ir=0; ir r_radial(Nr); + for( size_t ir=0; ir -double Conv_Coulomb_Pot_K::get_rmesh_proportion( - const Numerical_Orbital_Lm &orbs, - const double psi_threshold) -{ - for(int ir=orbs.getNr()-1; ir>=0; --ir) + Numerical_Orbital_Lm orbs_ccp; + orbs_ccp.set_orbital_info( + orbs.getLabel(), + orbs.getType(), + orbs.getL(), + orbs.getChi(), + Nr, + ModuleBase::GlobalFunc::VECTOR_TO_PTR(rab), + ModuleBase::GlobalFunc::VECTOR_TO_PTR(r_radial), + Numerical_Orbital_Lm::Psi_Type::Psik2, + ModuleBase::GlobalFunc::VECTOR_TO_PTR(psik2_ccp), + orbs.getNk(), + orbs.getDk(), + orbs.getDruniform(), + false, + true, GlobalV::CAL_FORCE); + return orbs_ccp; + } + + template<> + double get_rmesh_proportion( + const Numerical_Orbital_Lm &orbs, + const double psi_threshold) { - if(std::abs(orbs.getPsi(ir))>=psi_threshold) - return static_cast(ir)/orbs.getNr(); + for(int ir=orbs.getNr()-1; ir>=0; --ir) + { + if(std::abs(orbs.getPsi(ir))>=psi_threshold) + return static_cast(ir)/orbs.getNr(); + } + return 0.0; } - return 0.0; + } diff --git a/source/module_ri/conv_coulomb_pot_k.h b/source/module_ri/conv_coulomb_pot_k.h index 9adec9d915..d464a53f91 100644 --- a/source/module_ri/conv_coulomb_pot_k.h +++ b/source/module_ri/conv_coulomb_pot_k.h @@ -5,40 +5,37 @@ #include #include -class Conv_Coulomb_Pot_K +namespace Conv_Coulomb_Pot_K { -public: + enum class Ccp_Type{ // parameter: + Ccp, // + Hf, // "hf_Rcut" + Hse}; // "hse_omega" - enum class Ccp_Type{ // parameter: - Ccp, // - Hf, // - Hse}; // "hse_omega" - - template static T cal_orbs_ccp( + template T cal_orbs_ccp( const T &orbs, const Ccp_Type &ccp_type, const std::map ¶meter, - const double rmesh_times, - const int& nks); - -private: - - template< typename T > static double get_rmesh_proportion( + const double rmesh_times); + + //private: + template< typename T > double get_rmesh_proportion( const T &orbs, const double psi_threshold); - -private: - static std::vector cal_psi_ccp( const std::vector & psif ); - - static std::vector cal_psi_hf(const int& nks, const std::vector &psif, - const std::vector &k_radial, - const double omega); - - static std::vector cal_psi_hse( + //private: + std::vector cal_psi_ccp( + const std::vector & psif); + std::vector cal_psi_hf( + const std::vector &psif, + const std::vector &k_radial, + const double hf_Rcut); + std::vector cal_psi_hse( const std::vector & psif, const std::vector & k_radial, - const double omega); -}; + const double hse_omega); +} + +#include "conv_coulomb_pot_k.hpp" #endif \ No newline at end of file diff --git a/source/module_ri/conv_coulomb_pot_k.hpp b/source/module_ri/conv_coulomb_pot_k.hpp new file mode 100644 index 0000000000..5ca3abe5c8 --- /dev/null +++ b/source/module_ri/conv_coulomb_pot_k.hpp @@ -0,0 +1,37 @@ +#ifndef CONV_COULOMB_POT_K_HPP +#define CONV_COULOMB_POT_K_HPP + +#include "conv_coulomb_pot_k.h" +#include +#include + +namespace Conv_Coulomb_Pot_K +{ + + template< typename T > + std::vector cal_orbs_ccp( + const std::vector & orbs, + const Ccp_Type &ccp_type, + const std::map ¶meter, + const double rmesh_times) + { + std::vector orbs_ccp(orbs.size()); + for( size_t i=0; i!=orbs.size(); ++i ) + orbs_ccp[i] = cal_orbs_ccp(orbs[i], ccp_type, parameter, rmesh_times); + return orbs_ccp; + } + + template< typename T > + double get_rmesh_proportion( + const std::vector & orbs, + const double psi_threshold) + { + double rmesh_proportion=0; + for( const auto &orb : orbs ) + rmesh_proportion = std::max(rmesh_proportion, get_rmesh_proportion(orb,psi_threshold)); + return rmesh_proportion; + } + +} + +#endif \ No newline at end of file diff --git a/tests/integrate/204_NO_KP_NC_deltaspin/mulliken.txt.ref b/tests/integrate/204_NO_KP_NC_deltaspin/mulliken.txt.ref index 168c5723a3..dee68c8f40 100644 --- a/tests/integrate/204_NO_KP_NC_deltaspin/mulliken.txt.ref +++ b/tests/integrate/204_NO_KP_NC_deltaspin/mulliken.txt.ref @@ -3,92 +3,92 @@ CALCULATE THE MULLIkEN ANALYSIS FOR EACH ATOM Total charge: 32 Decomposed Mulliken populations 0 Zeta of Fe Spin 1 Spin 2 Spin 3 Spin 4 -s 0 1.317 0.05552 0.2843 0.02903 - sum over m 1.317 0.05552 0.2843 0.02903 -s 1 1.726 -0.01923 -0.09498 0.005159 - sum over m 1.726 -0.01923 -0.09498 0.005159 -s 2 0.03246 -0.04333 -0.2148 0.008137 - sum over m 0.03246 -0.04333 -0.2148 0.008137 -s 3 -0.02921 0.005194 0.02641 0.001867 - sum over m -0.02921 0.005194 0.02641 0.001867 - sum over m+zeta 3.046 -0.001842 0.0009368 0.04419 -pz 0 2.034 -0.001185 -0.005932 1.545e-06 -px 0 2.033 -0.001283 -0.006419 1.538e-06 -py 0 2.033 -0.001188 -0.005944 1.543e-06 - sum over m 6.1 -0.003656 -0.01829 4.626e-06 -pz 1 -0.02622 0.0005602 0.002791 0 -px 1 -0.02639 0.0006145 0.003054 0 -py 1 -0.02603 0.0005563 0.00277 0 - sum over m -0.07864 0.001731 0.008615 0 - sum over m+zeta 6.021 -0.001925 -0.00968 5.611e-06 -dz^2 0 1.964 0.0008273 0.004131 4.077e-06 -dxz 0 1.044 0.1755 0.7507 0.002258 -dyz 0 0.9544 0.1768 0.7532 0.002329 -dx^2-y^2 0 1.967 0.0007523 0.003756 3.978e-06 -dxy 0 1.055 0.1751 0.7495 0.002251 - sum over m 6.984 0.529 2.261 0.006846 -dz^2 1 0.03863 -0.0008699 -0.004363 5.197e-06 -dxz 1 -0.03759 -0.005346 -0.01936 -0.0001322 -dyz 1 -0.03407 -0.005734 -0.02118 -0.0001342 -dx^2-y^2 1 0.03943 -0.0009093 -0.004564 5.691e-06 -dxy 1 -0.03787 -0.005246 -0.0189 -0.0001314 - sum over m -0.03146 -0.01811 -0.06836 -0.000387 - sum over m+zeta 6.952 0.5109 2.193 0.006459 -fz^3 0 -0.007049 0.0007578 0.003775 0 -fxz^2 0 -0.002045 0.0002638 0.001312 0 -fyz^2 0 -0.002729 0.0002912 0.001448 0 -fzx^2-zy^2 0 6.273e-05 0 -6.642e-06 0 -fxyz 0 1.153e-05 1.446e-06 5.675e-06 0 -fx^3-3*xy^2 0 -0.00338 0.00044 0.002189 0 -f3yx^2-y^3 0 -0.00407 0.0004646 0.002311 0 - sum over m -0.0192 0.002219 0.01103 2.581e-06 - sum over m+zeta -0.0192 0.002219 0.01103 2.581e-06 +s 0 1.317 0.05552 -0.2843 0.02903 + sum over m 1.317 0.05552 -0.2843 0.02903 +s 1 1.726 -0.01923 0.09498 0.005159 + sum over m 1.726 -0.01923 0.09498 0.005159 +s 2 0.03246 -0.04333 0.2148 0.008137 + sum over m 0.03246 -0.04333 0.2148 0.008137 +s 3 -0.02921 0.005194 -0.02641 0.001867 + sum over m -0.02921 0.005194 -0.02641 0.001867 + sum over m+zeta 3.046 -0.001842 -0.0009368 0.04419 +pz 0 2.034 -0.001185 0.005932 1.545e-06 +px 0 2.033 -0.001283 0.006419 1.538e-06 +py 0 2.033 -0.001188 0.005944 1.543e-06 + sum over m 6.1 -0.003656 0.01829 4.626e-06 +pz 1 -0.02622 0.0005602 -0.002791 0 +px 1 -0.02639 0.0006145 -0.003054 0 +py 1 -0.02603 0.0005563 -0.00277 0 + sum over m -0.07864 0.001731 -0.008615 0 + sum over m+zeta 6.021 -0.001925 0.00968 5.611e-06 +dz^2 0 1.964 0.0008273 -0.004131 4.077e-06 +dxz 0 1.044 0.1755 -0.7507 0.002258 +dyz 0 0.9544 0.1768 -0.7532 0.002329 +dx^2-y^2 0 1.967 0.0007523 -0.003756 3.978e-06 +dxy 0 1.055 0.1751 -0.7495 0.002251 + sum over m 6.984 0.529 -2.261 0.006846 +dz^2 1 0.03863 -0.0008699 0.004363 5.197e-06 +dxz 1 -0.03759 -0.005346 0.01936 -0.0001322 +dyz 1 -0.03407 -0.005734 0.02118 -0.0001342 +dx^2-y^2 1 0.03943 -0.0009093 0.004564 5.691e-06 +dxy 1 -0.03787 -0.005246 0.0189 -0.0001314 + sum over m -0.03146 -0.01811 0.06836 -0.000387 + sum over m+zeta 6.952 0.5109 -2.193 0.006459 +fz^3 0 -0.007049 0.0007578 -0.003775 0 +fxz^2 0 -0.002045 0.0002638 -0.001312 0 +fyz^2 0 -0.002729 0.0002912 -0.001448 0 +fzx^2-zy^2 0 6.273e-05 0 6.642e-06 0 +fxyz 0 1.153e-05 1.446e-06 -5.675e-06 0 +fx^3-3*xy^2 0 -0.00338 0.00044 -0.002189 0 +f3yx^2-y^3 0 -0.00407 0.0004646 -0.002311 0 + sum over m -0.0192 0.002219 -0.01103 2.581e-06 + sum over m+zeta -0.0192 0.002219 -0.01103 2.581e-06 Total Charge on atom: Fe 16 -Total Magnetism on atom: Fe (0.5093, 2.195, 0.05066) +Total Magnetism on atom: Fe (0.5093, -2.195, 0.05066) 1 Zeta of Fe Spin 1 Spin 2 Spin 3 Spin 4 -s 0 1.275 0.05341 0.2605 -0.02903 - sum over m 1.275 0.05341 0.2605 -0.02903 -s 1 1.755 -0.01752 -0.08879 -0.005156 - sum over m 1.755 -0.01752 -0.08879 -0.005156 -s 2 -0.02898 -0.0404 -0.2039 -0.00813 - sum over m -0.02898 -0.0404 -0.2039 -0.00813 -s 3 -0.04711 0.006367 0.03139 -0.001874 - sum over m -0.04711 0.006367 0.03139 -0.001874 - sum over m+zeta 2.954 0.001862 -0.0008532 -0.04419 -pz 0 2.032 -0.001369 -0.006852 -1.367e-06 -px 0 2.025 -0.0009208 -0.004608 -1.387e-06 -py 0 2.032 -0.001332 -0.006666 -1.366e-06 - sum over m 6.089 -0.003622 -0.01813 -4.119e-06 -pz 1 -0.02528 0.0005889 0.002889 0 -px 1 -0.01606 0.0001369 0.0006408 0 -py 1 -0.02466 0.000571 0.002802 0 - sum over m -0.066 0.001297 0.006331 2.367e-06 - sum over m+zeta 6.023 -0.002325 -0.01179 -1.753e-06 -dz^2 0 1.957 0.001158 0.005774 -3.913e-06 -dxz 0 1.097 0.1724 0.7275 0.002311 -dyz 0 0.9509 0.1759 0.7475 0.002269 -dx^2-y^2 0 1.947 0.001654 0.008245 -4.075e-06 -dxy 0 1.113 0.1714 0.7227 0.002304 - sum over m 7.065 0.5225 2.212 0.006876 -dz^2 1 0.03925 -0.001062 -0.005333 -4.383e-06 -dxz 1 -0.0366 -0.003947 -0.01263 -0.0001213 -dyz 1 -0.03157 -0.005197 -0.01856 -0.0001267 -dx^2-y^2 1 0.04266 -0.001394 -0.007002 -4.206e-06 -dxy 1 -0.03743 -0.003854 -0.01222 -0.0001203 - sum over m -0.02369 -0.01545 -0.05575 -0.0003768 - sum over m+zeta 7.041 0.5071 2.156 0.006499 -fz^3 0 -0.006614 0.0007261 0.003596 0 -fxz^2 0 -0.001954 0.0002565 0.001276 0 -fyz^2 0 -0.002684 0.0002742 0.001366 0 -fzx^2-zy^2 0 9.09e-05 1.99e-05 8.018e-05 0 -fxyz 0 2.062e-05 4.102e-06 1.816e-05 0 -fx^3-3*xy^2 0 -0.003203 0.0004291 0.00213 0 -f3yx^2-y^3 0 -0.003698 0.0004635 0.002271 0 - sum over m -0.01804 0.002174 0.01074 0 - sum over m+zeta -0.01804 0.002174 0.01074 0 +s 0 1.275 0.05341 -0.2605 -0.02903 + sum over m 1.275 0.05341 -0.2605 -0.02903 +s 1 1.755 -0.01752 0.08879 -0.005156 + sum over m 1.755 -0.01752 0.08879 -0.005156 +s 2 -0.02898 -0.0404 0.2039 -0.00813 + sum over m -0.02898 -0.0404 0.2039 -0.00813 +s 3 -0.04711 0.006367 -0.03139 -0.001874 + sum over m -0.04711 0.006367 -0.03139 -0.001874 + sum over m+zeta 2.954 0.001862 0.0008532 -0.04419 +pz 0 2.032 -0.001369 0.006852 -1.367e-06 +px 0 2.025 -0.0009208 0.004608 -1.387e-06 +py 0 2.032 -0.001332 0.006666 -1.366e-06 + sum over m 6.089 -0.003622 0.01813 -4.119e-06 +pz 1 -0.02528 0.0005889 -0.002889 0 +px 1 -0.01606 0.0001369 -0.0006408 0 +py 1 -0.02466 0.000571 -0.002802 0 + sum over m -0.066 0.001297 -0.006331 2.367e-06 + sum over m+zeta 6.023 -0.002325 0.01179 -1.753e-06 +dz^2 0 1.957 0.001158 -0.005774 -3.913e-06 +dxz 0 1.097 0.1724 -0.7275 0.002311 +dyz 0 0.9509 0.1759 -0.7475 0.002269 +dx^2-y^2 0 1.947 0.001654 -0.008245 -4.075e-06 +dxy 0 1.113 0.1714 -0.7227 0.002304 + sum over m 7.065 0.5225 -2.212 0.006876 +dz^2 1 0.03925 -0.001062 0.005333 -4.383e-06 +dxz 1 -0.0366 -0.003947 0.01263 -0.0001213 +dyz 1 -0.03157 -0.005197 0.01856 -0.0001267 +dx^2-y^2 1 0.04266 -0.001394 0.007002 -4.206e-06 +dxy 1 -0.03743 -0.003854 0.01222 -0.0001203 + sum over m -0.02369 -0.01545 0.05575 -0.0003768 + sum over m+zeta 7.041 0.5071 -2.156 0.006499 +fz^3 0 -0.006614 0.0007261 -0.003596 0 +fxz^2 0 -0.001954 0.0002565 -0.001276 0 +fyz^2 0 -0.002684 0.0002742 -0.001366 0 +fzx^2-zy^2 0 9.09e-05 1.99e-05 -8.018e-05 0 +fxyz 0 2.062e-05 4.102e-06 -1.816e-05 0 +fx^3-3*xy^2 0 -0.003203 0.0004291 -0.00213 0 +f3yx^2-y^3 0 -0.003698 0.0004635 -0.002271 0 + sum over m -0.01804 0.002174 -0.01074 0 + sum over m+zeta -0.01804 0.002174 -0.01074 0 Total Charge on atom: Fe 16 -Total Magnetism on atom: Fe (0.5088, 2.154, -0.03769) +Total Magnetism on atom: Fe (0.5088, -2.154, -0.03769) diff --git a/tests/integrate/204_NO_KP_NC_deltaspin/result.ref b/tests/integrate/204_NO_KP_NC_deltaspin/result.ref index 649ae1ef31..a5ca4b1941 100644 --- a/tests/integrate/204_NO_KP_NC_deltaspin/result.ref +++ b/tests/integrate/204_NO_KP_NC_deltaspin/result.ref @@ -1,4 +1,4 @@ -etotref -6844.685232776227 -etotperatomref -3422.3426163881 +etotref -6844.685232778258 +etotperatomref -3422.3426163891 Compare_mulliken_pass 0 -totaltimeref 21.55 +totaltimeref 17.11