From 0da05355d44b145ad8b9fcbb2d27fa6f38ddc8c8 Mon Sep 17 00:00:00 2001 From: Yike Huang Date: Thu, 8 Jun 2023 18:05:57 +0800 Subject: [PATCH 1/2] Test: module_io/bessel_basis unit test still leave readin_C4 and allocate_C4 two functions untested, will complete soon. --- source/module_io/bessel_basis.cpp | 58 +- source/module_io/bessel_basis.h | 140 +++-- source/module_io/numerical_basis.cpp | 8 +- source/module_io/numerical_descriptor.cpp | 4 +- source/module_io/test/CMakeLists.txt | 24 +- source/module_io/test/INPUTs | 3 + source/module_io/test/bessel_basis_test.cpp | 537 ++++++++++++++++++ .../BesselBasis_UnitTest_C4_AtomType0.html | 25 + 8 files changed, 728 insertions(+), 71 deletions(-) create mode 100644 source/module_io/test/INPUTs create mode 100644 source/module_io/test/bessel_basis_test.cpp create mode 100644 source/module_io/test/support/BesselBasis_UnitTest_C4_AtomType0.html diff --git a/source/module_io/bessel_basis.cpp b/source/module_io/bessel_basis.cpp index 1170b2458f..e82d39baec 100644 --- a/source/module_io/bessel_basis.cpp +++ b/source/module_io/bessel_basis.cpp @@ -1,5 +1,5 @@ -#include "bessel_basis.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "bessel_basis.h" #include "module_base/parallel_common.h" #include "module_base/math_integral.h" #include "module_base/math_sphbes.h" @@ -25,25 +25,29 @@ void Bessel_Basis::init( const bool &smooth, const double &sigma, const double &rcut_in, - const double &tol_in, + const double &tol_in, + const UnitCell& ucell, const double &dk, - const double &dr) + const double &dr + ) { ModuleBase::TITLE("Bessel_Basis", "init"); this->Dk = dk; this->ecut = ecutwfc; this->rcut = rcut_in; this->tolerence = tol_in; - - //---------------------------------------------- - // setup Ecut_number - // ne * pi / rcut = sqrt(ecut) (Rydberg) - //---------------------------------------------- -// this->Ecut_number = static_cast( sqrt( 2.0 * ecut )* rcut/ModuleBase::PI );// hartree - this->Ecut_number = static_cast( sqrt( ecut )* rcut/ModuleBase::PI ); // Rydberg Unit. - assert( this->Ecut_number > 0 ); - - //------------------ + this->smooth = smooth; + this->sigma = sigma; + + //---------------------------------------------- + // setup Ecut_number + // ne * pi / rcut = sqrt(ecut) (Rydberg) + //---------------------------------------------- + // this->Ecut_number = static_cast( sqrt( 2.0 * ecut )* rcut/ModuleBase::PI );// hartree + this->Ecut_number = static_cast(sqrt(ecut) * rcut / ModuleBase::PI); // Rydberg Unit. + assert(this->Ecut_number > 0); + + //------------------ // Making a table //------------------ @@ -61,14 +65,13 @@ void Bessel_Basis::init( if( start_from_file ) { // setup C4 - this->allocate_C4(ntype, lmax_in, GlobalC::ucell.nmax, Ecut_number); - + this->allocate_C4(ntype, lmax_in, ucell.nmax, Ecut_number, ucell); // check tolerence - this->readin_C4("INPUTs", ntype, ecut, rcut, Ecut_number, tolerence); + this->readin_C4("INPUTs", ntype, ecut, rcut, Ecut_number, tolerence, ucell); #ifdef __MPI Parallel_Common::bcast_double( C4.ptr, C4.getSize() ); #endif - this->init_Faln(ntype, lmax_in, GlobalC::ucell.nmax, Ecut_number); + this->init_Faln(ntype, lmax_in, ucell.nmax, Ecut_number, ucell); } return; @@ -122,7 +125,8 @@ void Bessel_Basis::init_Faln( const int &ntype, const int &lmax, const int &nmax, - const int &ecut_number) + const int &ecut_number, + const UnitCell& ucell) { ModuleBase::TITLE("Bessel_Basis","init_Faln"); ModuleBase::timer::tick("Spillage","init_Faln"); @@ -133,9 +137,9 @@ void Bessel_Basis::init_Faln( this->nwfc = 0; for(int it=0; it> filec4; - for(int il=0; il< GlobalC::ucell.atoms[it].nwl+1; il++) + for(int il=0; il< ucell.atoms[it].nwl+1; il++) { - for(int in=0; in< GlobalC::ucell.atoms[it].l_nchi[il]; in++) + for(int in=0; in< ucell.atoms[it].l_nchi[il]; in++) { //for tests //std::cout << "\n" << std::setw(5) << it << std::setw(5) << il << std::setw(5) << in; @@ -446,7 +451,8 @@ void Bessel_Basis::allocate_C4( const int &ntype, const int &lmax, const int &nmax, - const int &ecut_number) + const int &ecut_number, + const UnitCell& ucell) { ModuleBase::TITLE("Bessel_Basis","allocate_C4"); @@ -454,9 +460,9 @@ void Bessel_Basis::allocate_C4( for(int it=0; it=rcut, j_l(q*r) = 0 (if not smoothed) double rcut; + /// @brief energy cutoff for determining kmesh and number of SBFs double ecut; double tolerence; - - bool smooth; // mohan add 2009-01-18 + /// @brief whether smooth SBFs around cutoff radius, resulting in non-zero values. For importance of smooth of SBFs, see J. Phys.: Condens. Matter 22 (2010) 445501, eqn 6. (mohan add 2009-01-18) + bool smooth; + /// @brief stddev of smooth function (Gaussian function, centered at rcut) double sigma; - //======================================================== - // MEMBER FUNCTIONS : - // NAME : readin ( read in parameters ) - //======================================================== + /// @brief Allocate memory for C4 matrix and initialize all elements to one. + /// @param ntype number of atom types + /// @param lmax maximal angular momentum of localized orbitals + /// @param nmax maximal principal quantum number of localized orbitals + /// @param ecut_number number of SBFs void allocate_C4( const int &ntype, const int &lmax, const int &nmax, - const int &ecut_number); - + const int &ecut_number, + const UnitCell& ucell + ); + + /// @brief Read C4 from external file. Presently an O(N^2) search algorithm is used. A HTML parser is needed in the future to improve performance. + /// @param name name of external file where C4-stored file information is contained + /// @param ntype number of atom types + /// @param ecut energy cutoff + /// @param rcut cutoff radius + /// @param ecut_number number of SBFs + /// @param tolerence accurancy of SBFs, here only used for consistency check void readin_C4( const std::string &name, const int &ntype, const int &ecut, const int &rcut, const int &ecut_number, - const double &tolerence); + const double &tolerence, + const UnitCell& ucell + ); void init_TableOne(void); - void init_Faln( const int &ntype,const int &lmax,const int &nmax,const int &ecut_number); - int nwfc; // number of localized wave functions - - // init table, used for output overlap Q. + /// @brief calculate F_{aln}(it, il, in, ik) = sum_{ie}{C4(it, il, in, ie)*TableOne(il, ie, ik)}, where TableOne is overlap integral between two spherical bessel functions (jle(r) and jlk(r)) + /// @param ntype number of atomtype + /// @param lmax maximal angular momentum + /// @param nmax maximal chi + /// @param ecut_number number of SBFs + void init_Faln( + const int &ntype, + const int &lmax, + const int &nmax, + const int &ecut_number, + const UnitCell& ucell + ); + + /// @brief number of localized wave functions + int nwfc; + + /// @brief calculate element value of TableOne matrix + /// @details (be called in Bessel_Basis::init(), used for outputing overlap Q matrix) initialize the table whose matrix element is the result of integral int{dr r^2 jle(r)*jlk(r)}, TableOne has three subscript (l, ie, ik), the first runs over orbitals' angular momentum and ie, ik run over ecut_number and kmesh SBFs + /// @param smooth_in whether jle(r) SBF is smoothed by a Gaussian function + /// @param sigma_in stddev for controlling smearing of Gaussian function for smoothing jle(r) + /// @param ecutwfc planewave kinetic energy cutoff for controlling kspace sampling + /// @param rcut cutoff radius of SBFs + /// @param dr realspace grid + /// @param dk kspace grid + /// @param lmax maximal angular momentum for SBFs + /// @param ecut_number number of SBFs + /// @param tolerence accurancy of SBFs void init_TableOne( const bool smooth_in, const double &sigma_in, diff --git a/source/module_io/numerical_basis.cpp b/source/module_io/numerical_basis.cpp index ad4bc0eaf4..36cde008e6 100644 --- a/source/module_io/numerical_basis.cpp +++ b/source/module_io/numerical_basis.cpp @@ -43,7 +43,9 @@ void Numerical_Basis::start_from_file_k(const int& ik, ModuleBase::ComplexMatrix INPUT.bessel_nao_smooth, INPUT.bessel_nao_sigma, INPUT.bessel_nao_rcut, - INPUT.bessel_nao_tolerence ); + INPUT.bessel_nao_tolerence, + GlobalC::ucell + ); this->mu_index = this->init_mu_index(); this->init_label = true; } @@ -71,7 +73,9 @@ void Numerical_Basis::output_overlap(const psi::Psi>& psi, INPUT.bessel_nao_smooth, INPUT.bessel_nao_sigma, INPUT.bessel_nao_rcut, - INPUT.bessel_nao_tolerence ); + INPUT.bessel_nao_tolerence, + GlobalC::ucell + ); this->mu_index = this->init_mu_index(); this->init_label = true; } diff --git a/source/module_io/numerical_descriptor.cpp b/source/module_io/numerical_descriptor.cpp index fe91ecd690..9e269745ac 100644 --- a/source/module_io/numerical_descriptor.cpp +++ b/source/module_io/numerical_descriptor.cpp @@ -53,7 +53,9 @@ void Numerical_Descriptor::output_descriptor(const psi::Psi INPUT.bessel_descriptor_smooth, INPUT.bessel_descriptor_sigma, rcut_in, - tol_in ); + tol_in, + GlobalC::ucell + ); this->nmax = Numerical_Descriptor::bessel_basis.get_ecut_number(); this->init_mu_index(); this->init_label = true; diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index a7d2903815..de2c90c48a 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -54,27 +54,25 @@ AddTest( ) AddTest( - TARGET io_dos_test + TARGET io_cal_dos LIBS ${math_libs} base device - SOURCES dos_test.cpp ../dos.cpp + SOURCES cal_dos_test.cpp ../cal_dos.cpp ) AddTest( TARGET io_write_dos_pw LIBS ${math_libs} base device symmetry - SOURCES write_dos_pw_test.cpp ../dos.cpp ../write_dos_pw.cpp ../output.cpp ../../module_cell/parallel_kpoints.cpp ../../module_cell/klist.cpp + SOURCES write_dos_pw_test.cpp ../cal_dos.cpp ../write_dos_pw.cpp ../output.cpp ../../module_cell/parallel_kpoints.cpp ../../module_cell/klist.cpp ) AddTest( TARGET io_print_info - LIBS ${math_libs} base device + LIBS ${math_libs} base device symmetry SOURCES print_info_test.cpp ../print_info.cpp ../input.cpp ../output.cpp ../../module_cell/unitcell.cpp - ../../module_cell/read_atoms.cpp ../../module_cell/read_cell_pseudopots.cpp ../../module_cell/atom_spec.cpp - ../../module_cell/atom_pseudo.cpp ../../module_cell/pseudo_nc.cpp ../../module_cell/read_pp.cpp - ../../module_cell/read_pp_upf201.cpp ../../module_cell/read_pp_upf100.cpp ../../module_cell/read_pp_vwr.cpp - ../../module_cell/read_pp_blps.cpp ../../module_cell/klist.cpp ../../module_cell/parallel_kpoints.cpp - ../../module_cell/module_symmetry/symm_other.cpp ../../module_cell/module_symmetry/symmetry_basic.cpp - ../../module_cell/module_symmetry/symmetry.cpp + ../../module_cell/read_atoms.cpp ../../module_cell/read_cell_pseudopots.cpp ../../module_cell/atom_spec.cpp + ../../module_cell/atom_pseudo.cpp ../../module_cell/pseudo_nc.cpp ../../module_cell/read_pp.cpp + ../../module_cell/read_pp_upf201.cpp ../../module_cell/read_pp_upf100.cpp ../../module_cell/read_pp_vwr.cpp + ../../module_cell/read_pp_blps.cpp ../../module_cell/klist.cpp ../../module_cell/parallel_kpoints.cpp ) AddTest( @@ -113,3 +111,9 @@ AddTest( TARGET io_parse_args SOURCES parse_args_test.cpp ) + +AddTest( + TARGET io_bessel_basis_test + LIBS ${math_libs} base device + SOURCES bessel_basis_test.cpp ../bessel_basis.cpp +) \ No newline at end of file diff --git a/source/module_io/test/INPUTs b/source/module_io/test/INPUTs new file mode 100644 index 0000000000..618420b3b5 --- /dev/null +++ b/source/module_io/test/INPUTs @@ -0,0 +1,3 @@ + +./support/BesselBasis_UnitTest_C4_AtomType0.html + diff --git a/source/module_io/test/bessel_basis_test.cpp b/source/module_io/test/bessel_basis_test.cpp new file mode 100644 index 0000000000..280c319257 --- /dev/null +++ b/source/module_io/test/bessel_basis_test.cpp @@ -0,0 +1,537 @@ +/// @details unit test for module_io/bessel_basis, not tested functions: readin_C4, allocate_C4 + +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include "../bessel_basis.h" +#include "../../module_cell/unitcell.h" +#include "../../module_elecstate/magnetism.h" + +#ifdef __LCAO +#include "../../module_cell/setup_nonlocal.h" +#endif +#include "gtest/gtest.h" + +/// @brief Simpson integral +/// @attention this function is a COMPLETE version, but there is no improvement in performance. +/// @param x variable stored in a vector +/// @param y function value stored in a vector +/// @return integral value +double SimpsonIntegral(const std::vector &x, const std::vector &y) +{ + double result = 0.0; + result += (y[0] + y[y.size() - 1]); + /* x and y must have the same size and their length must be a same odd number */ + assert(x.size() == y.size()); + assert(x.size() % 2 == 1); + double h = x[1] - x[0]; + for (int i = 1; i < x.size() - 1; i+=2) + { + result += (2*y[i] + 4*y[i+1]); + } + result *= h/3; + return result; +} +/// @brief recursive definition of 1st-kind Spherical Bessel function of the first kind +/// @attention this function is a COMPLETE version, can replace the one in basic mathematical library +/// @param l order of 1st-kind spherical Bessel function +/// @param x variable +/// @return value of 1st-kind spherical Bessel function +double SphericalBessel(int l, double x) { + /* c++ cannot handle limits well */ + if (x == 0) { + if (l == 0) { + return 1; + } + else { + return 0; + } + } + else { + if (l == 0) { + return std::sin(x) / x; + } + else if (l == 1) { + return (std::sin(x) / (x * x)) - std::cos(x) / x; + } + else { + return ((2 * l - 1) / x) * SphericalBessel(l - 1, x) - SphericalBessel(l - 2, x); + } + } +} +/// @brief 1st-kind Spherical Bessel function of the first kind mapped on realspace grid +/// @attention this function is a COMPLETE version, can replace the one in basic mathematical library +/// @param l order of 1st-kind spherical Bessel function +/// @param q wave vector +/// @param r realspace grid +/// @return function value on realspace grid +std::vector CalculateSphericalBessel(int l, double q, const std::vector& r) { + std::vector bessel; + for (const auto& radius : r) { + double x = q * radius; + double besselValue = SphericalBessel(l, x); + bessel.push_back(besselValue); + } + return bessel; +} +/// @brief Get the first p zeros of 1st-kind, l-ordered Spherical Bessel function from a constant table +/// @attention this function is a INCOMPLETE version, due to limited support of numerical table. +/// @param order l, angular momentum +/// @param number p, number of zeros needed +/// @return a vector of q, the q is from j_l(qr) +std::vector GetSphericalBesselZeros(int order, int number) { + std::map, double> zeros; + + zeros[{0, 1}] = 3.14159; zeros[{0, 2}] = 6.28318; zeros[{0, 3}] = 9.42477; + zeros[{0, 4}] = 12.5664; zeros[{0, 5}] = 15.708; zeros[{0, 6}] = 18.8495; + zeros[{0, 7}] = 21.9911; zeros[{0, 8}] = 25.1327; zeros[{0, 9}] = 28.2743; + zeros[{0, 10}] = 31.4159; + + zeros[{1, 1}] = 4.49341; zeros[{1, 2}] = 7.72525; zeros[{1, 3}] = 10.9041; + zeros[{1, 4}] = 14.0662; zeros[{1, 5}] = 17.2208; zeros[{1, 6}] = 20.3713; + zeros[{1, 7}] = 23.5181; zeros[{1, 8}] = 26.6617; zeros[{1, 9}] = 29.8029; + zeros[{1, 10}] = 32.9425; + + zeros[{2, 1}] = 5.76346; zeros[{2, 2}] = 9.09501; zeros[{2, 3}] = 12.3229; + zeros[{2, 4}] = 15.5146; zeros[{2, 5}] = 18.6861; zeros[{2, 6}] = 21.8457; + zeros[{2, 7}] = 24.9989; zeros[{2, 8}] = 28.1498; zeros[{2, 9}] = 31.2997; + zeros[{2, 10}] = 34.4491; + + zeros[{3, 1}] = 7.01559; zeros[{3, 2}] = 10.4013; zeros[{3, 3}] = 13.5821; + zeros[{3, 4}] = 16.7496; zeros[{3, 5}] = 19.9023; zeros[{3, 6}] = 23.0446; + zeros[{3, 7}] = 26.1799; zeros[{3, 8}] = 29.3105; zeros[{3, 9}] = 32.4377; + zeros[{3, 10}] = 35.5629; + + zeros[{4, 1}] = 8.26356; zeros[{4, 2}] = 11.6209; zeros[{4, 3}] = 14.7965; + zeros[{4, 4}] = 17.9598; zeros[{4, 5}] = 21.113; zeros[{4, 6}] = 24.2583; + zeros[{4, 7}] = 27.3979; zeros[{4, 8}] = 30.5325; zeros[{4, 9}] = 33.6635; + zeros[{4, 10}] = 36.7914; + + zeros[{5, 1}] = 9.51045; zeros[{5, 2}] = 12.8377; zeros[{5, 3}] = 16.0106; + zeros[{5, 4}] = 19.1714; zeros[{5, 5}] = 22.3224; zeros[{5, 6}] = 25.4666; + zeros[{5, 7}] = 28.6055; zeros[{5, 8}] = 31.7408; zeros[{5, 9}] = 34.873; + zeros[{5, 10}] = 38.0025; + + std::vector result; + for (int i = 1; i <= number; ++i) { + result.push_back(zeros[{order, i}]); + } + return result; +} +/// @brief Get mod of q vector of Spherical Bessel functions, all q satisfy when r=`rcut`, j_l(qr)=0. +/// @details first solve the equation j_l(x) = 0, therefore get the table (l, k) -> x, where l is the order of SBF and k is the k-th zero of j_l(x). Then let x = q*rcut, therefore q = x/rcut, return it. +/// @attention this function itself is a COMPLETE version, while the function it called GetSphericalBesselZeros may be INCOMPLETE, due to limited support of numerical table. +/// @param order the angular momentum of Spherical Bessel functions +/// @param number number of q to be returned +/// @param rcut 'cutoff radius' of Spherical Bessel functions. When r=rcut, Spherical Bessel functions are zero, and for r>rcut, they are zero as required by concept of constructing localized atomic orbital. +/// @return a vector of q, the q is from j_l(qr) +std::vector GetqList(int order, int number, double rcut) { + std::vector qList; + std::vector zerosList = GetSphericalBesselZeros(order, number); + for (const auto& zero : zerosList) { + qList.push_back(zero / rcut); + } + return qList; +} +/// @brief overload operator * for vectors, elements will be multiplied one by one +/// @param x one vector +/// @param y another vector +/// @return vector after multiplication +std::vector operator*(const std::vector& x, const std::vector& y) { + std::vector result; + for (int i = 0; i < x.size(); ++i) { + result.push_back(x[i] * y[i]); + } + return result; +} +/// @brief init_TableOne for unit test +/// @attention this is a COMPLETE version of init_TableOne, can replace the one in module_io/bessel_basis.cpp +/// @param smooth whether to smooth the function (gaussian function) +/// @param sigma stddev of gaussian function for smoothing +/// @param ecutwfc control the number of Spheical Bessel functions +/// @param rcut cutoff radius of Spherical Bessel functions, r>rcut, Spherical Bessel functions are zero +/// @param lmax maximal angular momentum of Spherical Bessel functions +/// @param dr grid spacing of r +/// @param dk grid spacing of k +/// @return `std::vector>>` TableOne[angular momentum][index of Spherical Bessel function][Arbitrary wave vector k] +std::vector>> GenerateTableOne(const bool smooth, const double sigma, const double ecutwfc, const double rcut, const int lmax, const double dr, const double dk){ + std::vector rGrid; + std::vector SmoothFactor_rGrid; + int rGridNum = static_cast(rcut/dr) + 4; + if (rGridNum % 2 == 0) + { + rGridNum += 1; + } + for (int indexr=0; indexr kGrid; + int kGridNum = static_cast(sqrt(ecutwfc)/dk) + 4 + 1; + if (kGridNum % 2 == 0) + { + kGridNum += 1; + } + for (double indexk=0; indexk(sqrt(ecutwfc)*rcut/M_PI); + + std::vector>> TableOne; + + for (int l=0; l qList = GetqList(l, NumBesselFunction, rcut); + /* should initialize TableOne first */ + if (l == 0) + { + for (int l=0; l> Zeros2d; + for (int indexq=0; indexq Zeros1d; + for (int indexk=0; indexk jle_rGrid = CalculateSphericalBessel(l, qList[indexq], rGrid); + if (smooth) + { + jle_rGrid = jle_rGrid*SmoothFactor_rGrid; + } + for (int indexk=0; indexk jlk_rGrid = CalculateSphericalBessel(l, kGrid[indexk], rGrid); + std::vector function_rGrid = jlk_rGrid*jle_rGrid*rGrid*rGrid; + TableOne[l][indexq][indexk] = SimpsonIntegral(rGrid, function_rGrid); + } + } + } + return TableOne; +} +/// @brief Improved version of module_io/bessel_basis::readin_C4 and allocate_C4 functions, for generating C4 matrix but now with a higher speed on accessing elements +/// @attention function will read in total number of chi-s both from file and from input, and assert that they are the same. It is also just a INCOMPLETE version, for a complete version, HTML parser library will be included and the other parameter, NumAtomType, will also be used for calibrating data. +/// @param FileName name of external file where C4-stored file information is contained +/// @param NumAtomType number of atom types +/// @param l maximal angular momentum of localized orbitals +/// @param NumChi number of contracted spherical bessel functions to fit one atomic orbital +/// @param NumBesselFunction number of spherical bessel functions in one contracted spherical bessel function +/// @return a map whose key is a string of "atomtype l chi q", and value is the corresponding C4 value +std::unordered_map ReadinC4(const std::string &FileName, const int &NumAtomType, const int &l, const int &NumChi, const int &NumBesselFunction){ + std::unordered_map C4Map; + std::ifstream C4File(FileName); + /* plan to add an HTML parser library in the future... */ + std::string word; + bool b_ReadC4Line = false; + int TotalNumChi = 0; + + assert(!C4File.fail()); + while (C4File.good()) + { + if (!b_ReadC4Line){ + C4File >> word; + if (word == "") + { + b_ReadC4Line = true; + C4File >> word; /* n_chi value read */ + TotalNumChi = std::stoi(word); + assert(TotalNumChi == NumChi); + continue; + } + } + else{ + for (int indexchi = 0; indexchi < TotalNumChi; indexchi++){ + + C4File >> word; + C4File >> word; + C4File >> word; /* skip title1, 2 and 3 */ + + C4File >> word; std::string key = word; key += " "; + C4File >> word; key += word; key += " "; + C4File >> word; key += word; key += " "; + + for (int indexNumBesselFunction = 0; indexNumBesselFunction < NumBesselFunction; indexNumBesselFunction++) + { + std::string keyForMap = key + std::to_string(indexNumBesselFunction); + C4File >> word; + C4Map[keyForMap] = std::stod(word); + } + } + break; + } + } + C4File.close(); + return C4Map; +} +/// @brief Generate F_{alpha, l, chi, k} matrix +/// @param FileName name of external file where C4-stored file information is contained +/// @param NumAtomType number of atom types +/// @param l maximal angular momentum of localized orbitals +/// @param NumChi number of contracted spherical bessel functions to fit one atomic orbital +/// @param NumBesselFunction number of spherical bessel functions in one contracted spherical bessel function +/// @param vvv_d_TableOne Integral table, has subscripts (l, q, k), whose element is the result of integral int{dr r^2 jle(r)*jlk(r)}. l runs over angular momentum, q runs over all spherical bessel functions, k runs over k points sampled controlled by ecutwfc and dk by sqrt(ecutwfc)/dk + 1 + 4. +/// @return F_{alpha, l, chi, k} matrix, whose element is the result of sum_{q}{C4(alpha, l, chi, q)*TableOne(l, q, k)} +std::vector>>> GenerateFaln(const std::string &FileName, const int &NumAtomType, const int &lmax, const int &NumChi, const int &NumBesselFunction, std::vector>> vvv_d_TableOne){ + std::unordered_map umap_str_d_C4Map = ReadinC4(FileName, NumAtomType, lmax, NumChi, NumBesselFunction); + std::vector>>> vvvv_d_Faln; + for (int indexAtomType = 0; indexAtomType < NumAtomType; indexAtomType++) + { + std::vector>> vvv_d_Faln; + for (int indexl = 0; indexl < lmax+1; indexl++) + { + std::vector> vv_d_Faln; + for (int indexChi = 0; indexChi < NumChi; indexChi++) + { + std::vector v_d_Faln; + for (int indexk = 0; indexk < vvv_d_TableOne[0][0].size(); indexk++) + { + double Faln = 0.0; + for (int indexBesselFunction = 0; indexBesselFunction < NumBesselFunction; indexBesselFunction++) + { + std::string key = std::to_string(indexAtomType) + " " + std::to_string(indexl) + " " + std::to_string(indexChi) + " " + std::to_string(indexBesselFunction); + Faln += umap_str_d_C4Map[key]*vvv_d_TableOne[indexl][indexBesselFunction][indexk]; + } + v_d_Faln.push_back(Faln); + } + vv_d_Faln.push_back(v_d_Faln); + } + vvv_d_Faln.push_back(vv_d_Faln); + } + vvvv_d_Faln.push_back(vvv_d_Faln); + } + return vvvv_d_Faln; +} + +/* OVERLOAD (de)constructors... */ +pseudo_nc::pseudo_nc() +{ +} +pseudo_nc::~pseudo_nc() +{ +} +Atom::Atom() +{ +} +Atom::~Atom() +{ +} +Atom_pseudo::Atom_pseudo() +{ +} +Atom_pseudo::~Atom_pseudo() +{ +} +UnitCell::UnitCell() +{ +} +UnitCell::~UnitCell() +{ +} +Magnetism::Magnetism() +{ +} +Magnetism::~Magnetism() +{ +} + +#ifdef __LCAO +InfoNonlocal::InfoNonlocal() +{ +} +InfoNonlocal::~InfoNonlocal() +{ +} +#endif +/* OVERLOAD printM3 function? */ +/* +void output::printM3(std::ofstream &ofs, const std::string &description, const ModuleBase::Matrix3 &m) +{ +} +*/ +#define private public +class TestBesselBasis : public ::testing::Test { +protected: + UnitCell ucell; + Bessel_Basis besselBasis; + + void SetUp() override { + ucell.ntype = 1; + ucell.lmax = 0; + ucell.nmax = 1; + + ucell.atoms = new Atom[1]; + ucell.atoms[0].l_nchi = new int[1]; + + ucell.atoms[0].nwl = 0; + ucell.atoms[0].l_nchi[0] = 1; + /* setup_cell manually */ + } + void TearDown() override { + delete[] ucell.atoms->l_nchi; + delete[] ucell.atoms; + } +}; + + TEST_F(TestBesselBasis, InitTest) { + /* parameters required by Bessel_Basis::init() */ + bool b_TestFaln = false; + double d_EnergyCutoff = 4; + int i_Ntype = 1; + int i_Lmax = 0; + bool b_Smooth = false; + double d_SmoothSigma = 0.1; + double d_CutoffRadius = 2; + double d_Tolerance = 0.01; + double d_dk = 0.01; + double d_dr = 2; + /* number of Spherical Bessel functions with given order l, + used to fit one chi */ + int i_Nq = static_cast(sqrt(d_EnergyCutoff)*d_CutoffRadius/M_PI); + /* number of SBF is expected to be 1 */ + besselBasis.init( + b_TestFaln, d_EnergyCutoff, i_Ntype, i_Lmax, b_Smooth, + d_SmoothSigma, d_CutoffRadius, d_Tolerance, + ucell, d_dk, d_dr + ); + EXPECT_EQ(besselBasis.get_ecut_number(), i_Nq); + EXPECT_EQ(besselBasis.get_ecut(), d_EnergyCutoff); + EXPECT_EQ(besselBasis.get_rcut(), d_CutoffRadius); + EXPECT_EQ(besselBasis.get_tolerence(), d_Tolerance); + EXPECT_EQ(besselBasis.get_smooth(), b_Smooth); + EXPECT_EQ(besselBasis.get_sigma(), d_SmoothSigma); +} + +TEST_F(TestBesselBasis, PolynomialInterpolation2Test) { + /* function Bessel_Basis::Polynomial_Interpolation2 is to do + cubic interpolation on TableOne, this matrix has, dimension l*nq*nk */ + /* parameters required by Bessel_Basis::init() */ + bool b_TestFaln = false; + double d_EnergyCutoff = 4; + int i_Ntype = 1; + int i_Lmax = 0; + bool b_Smooth = false; + double d_SmoothSigma = 0.1; + double d_CutoffRadius = 2.0; + double d_Tolerance = 0.01; + double d_dk = 2.0; + double d_dr = 0.01; /* for d_dr will largely impact accurancy of numerical integration, will not give a toy-value */ + + /* number of angular momentum considered should be 1 (lmax = 0)->The first dimension */ + int i_Nq = static_cast(sqrt(d_EnergyCutoff)*d_CutoffRadius/M_PI); + /* number of SBF is expected to be 1->The second dimension */ + int i_kMesh = static_cast(sqrt(d_EnergyCutoff)/d_dk) + 4 + 1; + /* the third dimension is at least to be 6 */ + /* therefore the expected dimension of TableOne is 1*1*6 */ + + besselBasis.init( + b_TestFaln, d_EnergyCutoff, i_Ntype, i_Lmax, b_Smooth, + d_SmoothSigma, d_CutoffRadius, d_Tolerance, + ucell, d_dk, d_dr + ); + /* gnorm for interpolation */ + double d_Gnorm = 1.0; + double d_position = d_Gnorm/d_dk; + int i_position = static_cast(d_position); + assert(i_position < i_kMesh-4); + double d_x0 = d_position - static_cast(i_position); + double d_x1 = 1.0 - d_x0; + double d_x2 = 2.0 - d_x0; + double d_x3 = 3.0 - d_x0; + + std::vector>> vvv_d_TableOne = GenerateTableOne( + b_Smooth, d_SmoothSigma, d_EnergyCutoff, d_CutoffRadius, + i_Lmax, d_dr, d_dk + ); + double d_yExpected = vvv_d_TableOne[0][0][i_position]*d_x1*d_x2*d_x3/6.0+ + vvv_d_TableOne[0][0][i_position]*d_x0*d_x2*d_x3/2.0- + vvv_d_TableOne[0][0][i_position]*d_x1*d_x0*d_x3/2.0+ + vvv_d_TableOne[0][0][i_position]*d_x1*d_x2*d_x0/6.0; + double d_yTested = besselBasis.Polynomial_Interpolation2(0, 0, d_Gnorm); + EXPECT_NEAR(d_yExpected, d_yTested, 0.01); +} +TEST_F(TestBesselBasis, PolynomialInterpolationTest) { + /* function Bessel_Basis::Polynomial_Interpolation is to do + cubic interpolation on Faln, this matrix has, dimension atomtype*l*nchi*nk. + To obtain Faln matrix, it is needed to first get C4 from external file. + The C4 is coefficent of SBF, has dimension atomtype*l*nchi*nq. + Therefore the Faln, equals the contraction of nq between C4 and TableOne. + C4: atomtype*l*nchi*nq + TableOne: l*nq*nk -> Faln: atomtype*l*nchi*nk */ + + /* parameters required by Bessel_Basis::init() */ + bool b_TestFaln = true; + double d_EnergyCutoff = 4; + int i_Ntype = 1; + int i_Lmax = 0; + bool b_Smooth = false; + double d_SmoothSigma = 0.1; + double d_CutoffRadius = 2.0; + double d_Tolerance = 0.01; + double d_dk = 2.0; + double d_dr = 0.01; /* for d_dr will largely impact accurancy of numerical integration, will not give a toy-value */ + + int i_Nq = static_cast(sqrt(d_EnergyCutoff)*d_CutoffRadius/M_PI); + int i_kMesh = static_cast(sqrt(d_EnergyCutoff)/d_dk) + 4 + 1; + /* + manipulate Bessel_Basis::init_Faln function + because for(int it=0; it(d_position); + assert(i_position < i_kMesh-4); + double d_x0 = d_position - static_cast(i_position); + double d_x1 = 1.0 - d_x0; + double d_x2 = 2.0 - d_x0; + double d_x3 = 3.0 - d_x0; + + std::vector>> vvv_d_TableOne = GenerateTableOne( + b_Smooth, d_SmoothSigma, d_EnergyCutoff, d_CutoffRadius, + i_Lmax, d_dr, d_dk + ); + std::vector>>> vvvv_d_Faln = GenerateFaln( + "./support/BesselBasis_UnitTest_C4_AtomType0.html", i_Ntype, i_Lmax, 1, i_Nq, vvv_d_TableOne + ); + double d_yExpected = vvvv_d_Faln[0][0][0][i_position]*d_x1*d_x2*d_x3/6.0+ + vvvv_d_Faln[0][0][0][i_position]*d_x0*d_x2*d_x3/2.0- + vvvv_d_Faln[0][0][0][i_position]*d_x1*d_x0*d_x3/2.0+ + vvvv_d_Faln[0][0][0][i_position]*d_x1*d_x2*d_x0/6.0; + double d_yTested = besselBasis.Polynomial_Interpolation(0, 0, 0, d_Gnorm); + EXPECT_NEAR(d_yExpected, d_yTested, 0.01); +} +#undef private \ No newline at end of file diff --git a/source/module_io/test/support/BesselBasis_UnitTest_C4_AtomType0.html b/source/module_io/test/support/BesselBasis_UnitTest_C4_AtomType0.html new file mode 100644 index 0000000000..05122a725e --- /dev/null +++ b/source/module_io/test/support/BesselBasis_UnitTest_C4_AtomType0.html @@ -0,0 +1,25 @@ + +Bessel Basis Unit Test. Present file will be read as C4 coefficients stored file of atom type 0. +INPUTS tag: +In INPUTS tag what enclosed are EnergyCutoff, CutoffRadius, Nq, Tolerance, respectively. +They should be totally identical with function Bessel_Basis::init_Faln(), or assert() will fail. +C4 tag: +The first line of this tag will be the total chi. The chi acts similarly with contracted Gaussian type +orbitals, which means for one chi, there will be Nq spherical Bessel functions, and for each angular +momentum l, there will be numbers of chi-s to form it. C4 file contains all information of chi-s, +i.e., all over l(s). For example the H, whose numerical orbitals file is named as H_gga_6au_100Ry_2s1p.orb, +which means there are 3 kinds of orbitals, 1s, 2s and 2p. The cutoff radius is 6 au, and the energy cutoff +is 100 Ry. For orbital "1s", there will be a chi_1s, then for "2s" the chi_2s, for "2p" the chi_2p, the +first line is thus the number chi_1s + chi_2s + chi_2p. +The following line has contents in sequence of: +title1 title2 title3 AtomType AngularMomentum ChiIndex C4Vaule1 C4Value2 ... C4ValueNq +Final description: +In unit test case, only one chi, one l is used, and EnergyCutoff, CutoffRadius are selected to have Nq=1. + + +4 2 1 0.01 + + +1 +TEST TEST TEST 0 0 0 1.0 + From 5d8e89e33085e37d03b8136f2b117aae35d867e1 Mon Sep 17 00:00:00 2001 From: Yike Huang <67682086+kirk0830@users.noreply.github.com> Date: Fri, 9 Jun 2023 09:31:53 +0800 Subject: [PATCH 2/2] Update bessel_basis_test.cpp add annotations to describe tested functions --- source/module_io/test/bessel_basis_test.cpp | 22 ++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/source/module_io/test/bessel_basis_test.cpp b/source/module_io/test/bessel_basis_test.cpp index 280c319257..1472cfd630 100644 --- a/source/module_io/test/bessel_basis_test.cpp +++ b/source/module_io/test/bessel_basis_test.cpp @@ -369,6 +369,26 @@ void output::printM3(std::ofstream &ofs, const std::string &description, const M { } */ + +/************************************************ + * unit test of bessel_basis.cpp + ***********************************************/ + +/* + * - Tested Functions: + * - InitTest: Bessel_Basis::init(start_from_file, ecutwfc, ntype, lmax_in, smooth, sigma, rcut_in, tol_in, ucell, dk, dr) + * - initialize Bessel_Basis class object with following steps: + * - call init_TableOne(...): calculate TableOne (l, ie, ik) = int{dr r^2 jle(r)*jlk(r)} (jle(r) and jlk(r) are Spherical Bessel functions) + * - return, unless set start_from_file = true, will also do the following: + * - call allocate_C4(...): allocate memory for C4 4d-matrix (it, il, in, ie) and set all elements to 1.0 + * - call readin_C4(...): read C4 value from external file and store in C4 4d-matrix + * - call init_Faln(...): calculate F_{aln}(it, il, in, ik) = sum_{ie}{C4(it, il, in, ie)*TableOne(il, ie, ik)} + * - PolynomialInterpolation2Test: Bessel_Basis::Polynomial_Interpolation2(const int &l, const int &ie, const double &gnorm) + * - return (cubic spline) interpolated element value of TableOne 3d-matrix + * - PolynomialInterpolationTest: Bessel_Basis::Polynomial_Interpolation(const int &it, const int &l, const int &ic, const double &gnorm) + * - return (cubic spline) interpolated element value of Faln 4d-matrix + */ + #define private public class TestBesselBasis : public ::testing::Test { protected: @@ -534,4 +554,4 @@ TEST_F(TestBesselBasis, PolynomialInterpolationTest) { double d_yTested = besselBasis.Polynomial_Interpolation(0, 0, 0, d_Gnorm); EXPECT_NEAR(d_yExpected, d_yTested, 0.01); } -#undef private \ No newline at end of file +#undef private