diff --git a/source/module_base/CMakeLists.txt b/source/module_base/CMakeLists.txt index 065c2cbe3c..dde98f915f 100644 --- a/source/module_base/CMakeLists.txt +++ b/source/module_base/CMakeLists.txt @@ -73,6 +73,7 @@ if(BUILD_TESTING) add_subdirectory(kernels/test) add_subdirectory(module_mixing/test) add_subdirectory(module_device/test) + add_subdirectory(grid/test) if (USE_ABACUS_LIBM) add_subdirectory(libm/test) endif() diff --git a/source/module_base/grid/delley.cpp b/source/module_base/grid/delley.cpp new file mode 100644 index 0000000000..1e0b6c703e --- /dev/null +++ b/source/module_base/grid/delley.cpp @@ -0,0 +1,376 @@ +#include "module_base/grid/delley.h" + +#include +#include +#include + +namespace { + +struct Table { + const int lmax_; + const int ngrid_; + const int ntype_[6]; + const std::vector data_; +}; + +// Delley's table from the original article +const std::vector table = { + { + 17, 110, {1, 1, 0, 3, 1, 0}, + { + 0.00000000000000000, 0.00000000000000000, 0.0038282704949371616, + 0.57735026918962576, 0.57735026918962576, 0.0097937375124875125, + 0.18511563534473617, 0.18511563534473617, 0.0082117372831911110, + 0.39568947305594191, 0.39568947305594191, 0.0095954713360709628, + 0.69042104838229218, 0.21595729184584883, 0.0099428148911781033, + 0.47836902881215020, 0.00000000000000000, 0.0096949963616630283, + } + }, + { + 23, 194, {1, 1, 1, 4, 1, 1}, + { + 0.00000000000000000, 0.00000000000000000, 0.0017823404472446112, + 0.57735026918962576, 0.57735026918962576, 0.0055733831788487380, + 0.70710678118654752, 0.00000000000000000, 0.0057169059499771019, + 0.44469331787174373, 0.44469331787174373, 0.0055187714672736137, + 0.28924656275754386, 0.28924656275754386, 0.0051582377118053831, + 0.67129734426952263, 0.31419699418258608, 0.0056087040825879968, + 0.12993354476500669, 0.12993354476500669, 0.0041067770281693941, + 0.34577021976112827, 0.00000000000000000, 0.0050518460646148085, + 0.52511857244364202, 0.15904171053835295, 0.0055302489162330937, + } + }, + { + 29, 302, {1, 1, 0, 6, 2, 2}, + { + 0.00000000000000000, 0.00000000000000000, 0.0008545911725128148, + 0.57735026918962576, 0.57735026918962576, 0.0035991192850255715, + 0.70117664160895449, 0.12923867271051493, 0.0036500458076772554, + 0.65663294102196118, 0.37103417838482119, 0.0036048226014198817, + 0.47290541325810046, 0.47290541325810046, 0.0035767296617433671, + 0.35156403455701051, 0.35156403455701051, 0.0034497884243058833, + 0.22196452362941784, 0.22196452362941784, 0.0031089531224136753, + 0.09618308522614784, 0.09618308522614784, 0.0023521014136891644, + 0.57189558918789607, 0.00000000000000000, 0.0036008209322164603, + 0.26441528870606625, 0.00000000000000000, 0.0029823449631718039, + 0.54486773725807738, 0.25100347517704651, 0.0035715405542733871, + 0.41277240831685310, 0.12335485325833274, 0.0033923122050061702, + } + }, + { + 35, 434, {1, 1, 1, 7, 2, 4}, + { + 0.00000000000000000, 0.00000000000000000, 0.0005265897968224436, + 0.57735026918962576, 0.57735026918962576, 0.0025123174189273072, + 0.70710678118654752, 0.00000000000600000, 0.0025482199720026072, + 0.69093463075091106, 0.21264682470755207, 0.0025304038011863550, + 0.64566647074242561, 0.40771266489776951, 0.0025132671745975644, + 0.49143426377847465, 0.49143426377847465, 0.0025017251684029361, + 0.39272597633680022, 0.39272597633680022, 0.0024453734373129800, + 0.28612890103076384, 0.28612890103076384, 0.0023026947822274158, + 0.17748360546091578, 0.17748360546091578, 0.0020142790209185282, + 0.07568084367178018, 0.07568084367178018, 0.0014624956215946138, + 0.21027252285730696, 0.00000000000000000, 0.0019109512821795323, + 0.47159869115131592, 0.00000000000000000, 0.0024174423756389808, + 0.33443631453434549, 0.09921769636429237, 0.0022366077604378487, + 0.45023303825826254, 0.20548236964030437, 0.0024169300443247753, + 0.55501523610768072, 0.31042840351665415, 0.0024966440545530860, + 0.59051570489252711, 0.10680182607580483, 0.0025122368545634951, + } + }, + { + 41, 590, {1, 1, 0, 8, 4, 6}, + { + 0.00000000000000000, 0.00000000000000000, 0.0001009005753378758, + 0.57735026918962576, 0.57735026918962576, 0.0018514016873890461, + 0.70404760433146996, 0.09291900596883211, 0.0018686219518306975, + 0.68084561988024238, 0.26999719217017240, 0.0018648696345606001, + 0.63723669159418917, 0.43342680786054810, 0.0018497643975168892, + 0.50447558060926046, 0.50447558060926046, 0.0018450277740822388, + 0.42175447334398773, 0.42175447334398773, 0.0018164174988262214, + 0.33201962086729379, 0.33201962086729379, 0.0017449464690023229, + 0.23917494336556047, 0.23917494336556047, 0.0016278016126848035, + 0.14024070738935403, 0.14024070738935403, 0.0015576827519901693, + 0.09161634328605240, 0.00000000000000000, 0.0012680968886048433, + 0.20326292518419433, 0.00000000000000000, 0.0011183965414769017, + 0.39364042372978295, 0.00000000000000000, 0.0017287035120530033, + 0.61262355812929648, 0.00000000000000000, 0.0018551905629473527, + 0.28114771623428322, 0.08959875911893791, 0.0014697353123693616, + 0.38175470908581117, 0.17327600238498666, 0.0016819651914742022, + 0.47452376478986998, 0.26422260656245780, 0.0017876372876796954, + 0.56127905075920534, 0.35189965873835832, 0.0018400735685528423, + 0.50324791996964975, 0.08886791018186295, 0.0018072536817113700, + 0.59768324320748616, 0.18154345643517542, 0.0018527289739424312, + } + }, + { + 47, 770, {1, 1, 1, 9, 4, 9}, + { + 0.00000000000000000, 0.00000000000000000, 0.0011685335608691628, + 0.57735026918962576, 0.57735026918962576, 0.0014121215930643264, + 0.70710678118654752, 0.00000000000000000, 0.0014468645950992776, + 0.11441365123336336, 0.11441365123336336, 0.0010478418864629224, + 0.19944675708548970, 0.19944675708548970, 0.0012392547584848484, + 0.28401278368259530, 0.28401278368259530, 0.0013259295792415379, + 0.36646411416548296, 0.36646411416548296, 0.0013756097758625958, + 0.44356118052513995, 0.44356118052513995, 0.0013999348863558624, + 0.51435709575333968, 0.51435709575333968, 0.0014096221218822673, + 0.63052081196671812, 0.45264446462279973, 0.0014108746499638577, + 0.67164784337293865, 0.31269529735024947, 0.0014134887639034478, + 0.69812332010174177, 0.15889512220405632, 0.0014366946685816802, + 0.12047667931264991, 0.00000000000000000, 0.0010901543574180667, + 0.30940302315480606, 0.00000000000000000, 0.0001869137844803852, + 0.34884276430183016, 0.00000000000000000, 0.0011284267652336505, + 0.53224214285417946, 0.00000000000000000, 0.0013844558026568455, + 0.23249923409267532, 0.06616159933437003, 0.0011853923885095502, + 0.32477344409682044, 0.14568618765136356, 0.0012949021664637693, + 0.41056989039349425, 0.22832839132127622, 0.0013525857420363760, + 0.49213658085114203, 0.30714431901543855, 0.0013925025908786082, + 0.56548849812588755, 0.38271180625074657, 0.0014073257894372725, + 0.43713473693946563, 0.07970715187939190, 0.0013128954307755017, + 0.52320749473197761, 0.15892620239864833, 0.0013784632898490457, + 0.60283033994386521, 0.23667220253873893, 0.0014125450609821936, + 0.62037164721742807, 0.07982328826030880, 0.0014289835314095131, + } + }, + { + 53, 974, {1, 1, 0, 12, 4, 12}, + { + 0.00000000000000000, 0.00000000800000000, 0.0001438294190527431, + 0.57735026918962576, 0.57735026918962576, 0.0011257722882870041, + 0.04292963545341347, 0.04292963545341347, 0.0004948029341949241, + 0.10514268540864042, 0.10514268540864042, 0.0007357990109125470, + 0.17500248676230874, 0.17500248676230874, 0.0008889132771304384, + 0.24776533796502568, 0.24776533796502568, 0.0009888347838921435, + 0.32065671239559574, 0.32065671239559574, 0.0010532996817094706, + 0.39165207498499835, 0.39165207498499835, 0.0010927788070145785, + 0.45908258741876237, 0.45908258741876237, 0.0011143893940632272, + 0.52145638884158605, 0.52145638884158605, 0.0011237247880515553, + 0.62531702446541989, 0.46685890569574328, 0.0011252393252438136, + 0.66379267445231699, 0.34461365423743822, 0.0011261532718159050, + 0.69104103984983007, 0.21195415185018465, 0.0011302869311238408, + 0.70529070074577603, 0.07162440144995566, 0.0011349865343639549, + 0.12366867626579899, 0.00000000000000000, 0.0006823367927109931, + 0.29407771144683870, 0.00000000000000000, 0.0009454158160447096, + 0.46977538492076491, 0.00000000000000000, 0.0010744299753856791, + 0.63345632411395669, 0.00000000000000000, 0.0011293000865691317, + 0.20291287527775228, 0.05974048614181342, 0.0008436884500901954, + 0.46026219424840539, 0.13757604084736365, 0.0010752557204488846, + 0.50306739996620357, 0.33910165263362857, 0.0011085772368644620, + 0.28176064224421343, 0.12716751914398195, 0.0009566475323783357, + 0.43315612917201574, 0.26931207404135125, 0.0010806632507173907, + 0.62561673585808142, 0.14197864526019183, 0.0011267971311962946, + 0.37983952168591567, 0.06709284600738255, 0.0010225687153580612, + 0.55175054214235205, 0.07057738183256172, 0.0011089602677131075, + 0.60296191561591869, 0.27838884778821546, 0.0011227906534357658, + 0.35896063295890958, 0.19795789389174069, 0.0010324018471174598, + 0.53486664381354765, 0.20873070611032740, 0.0011072493822838539, + 0.56749975460743735, 0.40551221378728359, 0.0011217800485199721, + } + }, + { + 59, 1202, {1, 1, 1, 13, 4, 16}, + { + 0.00006000000000000, 0.00000000000000000, 0.0001105189233267572, + 0.57735026918962576, 0.57735026918962576, 0.0009133159786443561, + 0.70710678118654752, 0.00000000000000000, 0.0009205232738090741, + 0.03712636449657089, 0.03712636449657089, 0.0003690421898017899, + 0.09140060412262223, 0.09140060412262223, 0.0005603990928680660, + 0.15310778524699062, 0.15310778524699062, 0.0006865297629282609, + 0.21809288916606116, 0.21809288916606116, 0.0007720338551145630, + 0.28398745322001746, 0.28398745322001746, 0.0008301545958894795, + 0.34911776009637644, 0.34911776009637644, 0.0008686692550179628, + 0.41214314614443092, 0.41214314614443092, 0.0008927076285846890, + 0.47189936271491266, 0.47189936271491266, 0.0009060820238568219, + 0.52731454528423366, 0.52731454528423366, 0.0009119777254940867, + 0.62094753324440192, 0.47838093807695216, 0.0009128720138604181, + 0.65697227118572905, 0.36983086645942597, 0.0009130714935691735, + 0.68417883090701434, 0.25258395570071777, 0.0009152873784554116, + 0.70126043301236308, 0.12832618665972300, 0.0009187436274321654, + 0.10723822154781661, 0.00000000000000000, 0.0005176977312965694, + 0.25820689594969680, 0.00006000000000000, 0.0007331143682101417, + 0.41727529553067168, 0.00000000000000000, 0.0008463232836379928, + 0.57003669117925033, 0.00000000000000000, 0.0009031122694253992, + 0.17717740226153253, 0.05210639477011284, 0.0006485778453163257, + 0.24757164634262876, 0.11156409571564867, 0.0007435030910982369, + 0.31736152466119767, 0.17465516775786261, 0.0008101731497468018, + 0.38542911506692237, 0.23902784793817240, 0.0008556299257311812, + 0.45074225931570644, 0.30294669735289819, 0.0008850282341265444, + 0.51235184864198708, 0.36498322605976536, 0.0009022692938426915, + 0.56937024984684411, 0.42386447815223403, 0.0009105760258970126, + 0.33546162890664885, 0.05905888853235508, 0.0007998527891839054, + 0.40902684270853572, 0.12172350510959870, 0.0008483389574594331, + 0.47853206759224352, 0.18575051945473351, 0.0008811048182425720, + 0.54343035696939004, 0.24941121623622365, 0.0009010091677105086, + 0.60311616930963100, 0.31122759471496082, 0.0009107813579482705, + 0.49322211848512846, 0.06266250624154169, 0.0008803208679738260, + 0.56321230207620997, 0.12677748006842827, 0.0009021342299040653, + 0.62698055090243917, 0.19060182227792370, 0.0009131578003189435, + 0.63942796347491023, 0.06424549224220589, 0.0009158016174693465, + } + } +}; // end of the definition of "table" + +// size of each group of points with octahedral symmetry +// 6: (1, 0, 0) x sign x permutation (vertices) +// 8: (sqrt(1/3), sqrt(1/3), sqrt(1/3)) x sign (face centers) +// 12: (sqrt(2)/2, sqrt(2)/2, 0) x sign x permutation (edge centers) +// 24a: (u, u, sqrt(1-2u^2)) x sign x permutation +// 24b: (u, 0, sqrt(1-u^2)) x sign x permutation +// 48: (u, v, sqrt(1-u^2-v^2)) x sign x permutation +const int group_size[] = {6, 8, 12, 24, 24, 48}; + +using Fill_t = std::function; + +// functors that fill the grid group-wise +const std::vector fill = { + // (1, 0, 0) x sign x permutation (vertices) + [](double* grid, double, double) { + for (int i = 0; i < 3; ++i) { + for (double one : {-1.0, 1.0}) { + grid[i] = one; + grid[(i+1)%3] = 0.0; + grid[(i+2)%3] = 0.0; + std::advance(grid, 3); + } + } + }, + + // (sqrt(1/3), sqrt(1/3), sqrt(1/3)) x sign (face centers) + [](double* grid, double, double) { + const double a = std::sqrt(3) / 3.0; + for (int xsign : {-1, 1}) { + for (int ysign : {-1, 1}) { + for (int zsign : {-1, 1}) { + grid[0] = xsign * a; + grid[1] = ysign * a; + grid[2] = zsign * a; + std::advance(grid, 3); + } + } + } + }, + + // (sqrt(2)/2, sqrt(2)/2, 0) x sign x permutation (edge centers) + [](double* grid, double, double) { + const double a = std::sqrt(2) / 2.0; + for (int i = 0; i < 3; ++i) { + for (int sign1 : {-1, 1}) { + for (int sign2 : {-1, 1}) { + grid[i] = 0; + grid[(i+1)%3] = sign1 * a; + grid[(i+2)%3] = sign2 * a; + std::advance(grid, 3); + } + } + } + }, + + // (u, u, sqrt(1-2u^2)) x sign x permutation + [](double* grid, double x, double y) { + double u = x == y ? x : std::sqrt(1.0 - x * x - y * y); + double v = std::sqrt(1.0 - 2.0 * u * u); + for (int i = 0; i < 3; ++i) { + for (int sign1 : {-1, 1}) { + for (int sign2 : {-1, 1}) { + for (int sign3 : {-1, 1}) { + grid[i] = sign1 * u; + grid[(i+1)%3] = sign2 * u; + grid[(i+2)%3] = sign3 * v; + std::advance(grid, 3); + } + } + } + } + }, + + // (u, 0, sqrt(1-u^2)) x sign x permutation + [](double* grid, double x, double y) { + double u = x > 0 ? x : y; + double v = std::sqrt(1.0 - u * u); + for (int i0 = 0; i0 < 3; ++i0) { + for (int iu0 : {1, 2}) { + for (int sign_u : {-1, 1}) { + for (int sign_v : {-1, 1}) { + grid[i0] = 0; + grid[(i0+iu0)%3] = sign_u * u; + grid[(i0-iu0+3)%3] = sign_v * v; + std::advance(grid, 3); + } + } + } + } + }, + + // (u, v, sqrt(1-u^2-v^2)) x sign x permutation + [](double* grid, double x, double y) { + double r = x; + double s = y; + double t = std::sqrt(1.0 - r * r - s * s); + for (int ir = 0; ir < 3; ++ir) { + for (int irs : {1, 2}) { + for (int sign_r : {-1, 1}) { + for (int sign_s : {-1, 1}) { + for (int sign_t : {-1, 1}) { + grid[ir] = sign_r * r; + grid[(ir+irs)%3] = sign_s * s; + grid[(ir-irs+3)%3] = sign_t * t; + std::advance(grid, 3); + } + } + } + } + } + }, +}; // end of the definition of "fill" + +const Table* _find_table(int& lmax) { + // NOTE: this function assumes elements in "Delley::table_" are + // arranged such that their members "lmax_" are in ascending order. + auto tab = std::find_if(table.begin(), table.end(), + [lmax](const Table& t) { return t.lmax_ >= lmax; }); + return tab == table.end() ? nullptr : (lmax = tab->lmax_, &(*tab)); +} + +void _get(const Table* tab, double* grid, double* weight) { + assert(tab); + const double* ptr = &tab->data_[0]; + for (int itype = 0; itype < 6; ++itype) { + int stride = group_size[itype]; + for (int i = 0; i < tab->ntype_[itype]; ++i, ptr += 3, + grid += 3*stride, weight += stride) { + fill[itype](grid, ptr[0], ptr[1]); + std::fill(weight, weight + stride, ptr[2]); + } + } +} + +} // end of anonymous namespace + + +int Grid::Delley::ngrid(int& lmax) { + auto tab = _find_table(lmax); + return tab ? tab->ngrid_ : -1; +} + + +int Grid::Delley::get(int& lmax, double* grid, double* weight) { + auto tab = _find_table(lmax); + return tab ? _get(tab, grid, weight), 0 : -1; +} + + +int Grid::Delley::get( + int& lmax, + std::vector& grid, + std::vector& weight +) { + auto tab = _find_table(lmax); + if (!tab) { + return -1; + } + grid.resize(3 * tab->ngrid_); + weight.resize(tab->ngrid_); + _get(tab, grid.data(), weight.data()); + return 0; +} diff --git a/source/module_base/grid/delley.h b/source/module_base/grid/delley.h new file mode 100644 index 0000000000..4f68a2168a --- /dev/null +++ b/source/module_base/grid/delley.h @@ -0,0 +1,57 @@ +#ifndef GRID_DELLEY_H +#define GRID_DELLEY_H + +#include +#include + +/** + * @brief Delley's grid for quadrature on the unit sphere. + * + * Reference: + * Delley, B. (1996). High order integration schemes on the unit sphere. + * Journal of computational chemistry, 17(9), 1152-1155. + * + */ +namespace Grid { +namespace Delley { + +/** + * @brief Number of Delley's grid points for a certain order of accuracy. + * + * This function finds the minimum Delley's grid with an accuracy + * of order "lmax". On exit, lmax is set to the order of this grid, and + * its corresponding number of grid points is returned. If no such grid + * is available, lmax is left unchanged and the function will return -1. + * + * For example, if lmax = 20 on input, the function will return 194 and + * lmax will be set to 23. + * + */ +int ngrid(int& lmax); + + +/** + * @brief Delley's quadrature grid and weights. + * + * This function retrieves the minimum Delley's grid with an accuray + * of order "lmax". On exit, lmax is set to the order of this grid, and + * the coordinates & weights are returned in "grid" & "weight". + * + * Coordinates are stored in the following order: + * + * x0, y0, z0, x1, y1, z1, x2, y2, z2, ... + * + * "grid" and "weight" must be pre-allocated to hold 3*ngrid(lmax) and + * ngrid(lmax) elements, respectively. The function will return 0 if + * successful, or -1 if the requested order cannot be fulfilled. + * + */ +int get(int& lmax, double* grid, double* weight); + +// a handy wrapper doing the same as above +int get(int& lmax, std::vector& grid, std::vector& weight); + +} // end of namespace Delley +} // end of namespace Grid + +#endif diff --git a/source/module_base/grid/test/CMakeLists.txt b/source/module_base/grid/test/CMakeLists.txt new file mode 100644 index 0000000000..0f77a425fe --- /dev/null +++ b/source/module_base/grid/test/CMakeLists.txt @@ -0,0 +1,9 @@ +remove_definitions(-D__MPI) + +AddTest( + TARGET test_delley + SOURCES test_delley.cpp + ../delley.cpp + ../../ylm.cpp +) + diff --git a/source/module_base/grid/test/test_delley.cpp b/source/module_base/grid/test/test_delley.cpp new file mode 100644 index 0000000000..2b892cd532 --- /dev/null +++ b/source/module_base/grid/test/test_delley.cpp @@ -0,0 +1,138 @@ +#include "module_base/grid/delley.h" +#include "module_base/ylm.h" + +#include "gtest/gtest.h" +#include +#ifdef __MPI +#include +#endif + +using namespace Grid; + +// mock the function to prevent unnecessary dependency +namespace ModuleBase { +void WARNING_QUIT(const std::string&, const std::string&) {} +} + +class DelleyTest: public ::testing::Test { +protected: + void randgen(int lmax, std::vector& coef); + const double tol = 1e-12; +}; + + +void DelleyTest::randgen(int lmax, std::vector& coef) { + coef.resize((lmax + 1) * (lmax + 1)); + + // fill coef with uniformly distributed random numbers + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution dis(0.0, 1.0); + for (size_t i = 0; i < coef.size(); ++i) { + coef[i] = dis(gen); + } + + // normalize the coefficients + double fac = 0.0; + for (size_t i = 0; i < coef.size(); ++i) { + fac += coef[i] * coef[i]; + } + + fac = 1.0 / std::sqrt(fac); + for (size_t i = 0; i < coef.size(); ++i) { + coef[i] *= fac; + } +} + + +TEST_F(DelleyTest, NumGrid) { + int lmax = 5; + int ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 17); + EXPECT_EQ(ngrid, 110); + + lmax = 17; + ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 17); + EXPECT_EQ(ngrid, 110); + + lmax = 20; + ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 23); + EXPECT_EQ(ngrid, 194); + + lmax = 59; + ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 59); + EXPECT_EQ(ngrid, 1202); + + lmax = 60; + ngrid = Delley::ngrid(lmax); + EXPECT_EQ(lmax, 60); + EXPECT_EQ(ngrid, -1); +} + + +TEST_F(DelleyTest, Accuracy) { + /* + * Given + * + * f = c[0]*Y00 + c[1]*Y10 + c[2]*Y11 + ..., + * + * where c[0], c[1], c[2], ... are some random numbers, the integration + * of |f|^2 on the unit sphere + * + * \int |f|^2 d\Omega = c[0]^2 + c[1]^2 + c[2]^2 + ... . + * + * This test verifies with the above integral that quadrature with + * Delley's grid is exact up to floating point errors. + * + */ + std::vector grid, weight, coef; + + for (int grid_lmax = 17; grid_lmax < 60; grid_lmax +=6) { + Delley::get(grid_lmax, grid, weight); + int func_lmax = grid_lmax / 2; + randgen(func_lmax, coef); + + double val = 0.0; + std::vector ylm_real; + for (size_t i = 0; i < weight.size(); i++) { + ModuleBase::Ylm::sph_harm(func_lmax, + grid[3*i], grid[3*i+1], grid[3*i+2], ylm_real); + double tmp = 0.0; + for (size_t i = 0; i < coef.size(); ++i) { + tmp += coef[i] * ylm_real[i]; + } + val += weight[i] * tmp * tmp; + } + val *= 4.0 * std::acos(-1.0); + + double val_ref = 0.0; + for (size_t i = 0; i < coef.size(); ++i) { + val_ref += coef[i] * coef[i]; + } + + double abs_diff = std::abs(val - val_ref); + EXPECT_LT(abs_diff, tol); + //printf("order = %2i val_ref = %8.5f abs_diff = %8.5e\n", + // grid_lmax, val_ref, abs_diff); + } +} + + +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; +} diff --git a/source/module_base/ylm.cpp b/source/module_base/ylm.cpp index a60ea9b11c..a36b14c60a 100644 --- a/source/module_base/ylm.cpp +++ b/source/module_base/ylm.cpp @@ -54,7 +54,7 @@ std::vector Ylm::ylmcoef = { // here Lmax == max angular momentum + 1 void Ylm::get_ylm_real( const int &Lmax, const ModuleBase::Vector3 &vec, double ylmr[] ) { - ModuleBase::timer::tick ("Ylm","get_ylm_real"); + //ModuleBase::timer::tick ("Ylm","get_ylm_real"); //1e-9 is too large const double cut0 = 1e-12; // allocate space. @@ -165,7 +165,7 @@ void Ylm::get_ylm_real( const int &Lmax, const ModuleBase::Vector3 &vec, } }// end do - ModuleBase::timer::tick ("Ylm", "get_ylm_real"); + //ModuleBase::timer::tick ("Ylm", "get_ylm_real"); return; } @@ -536,7 +536,7 @@ void Ylm::sph_harm std::vector &rly ) { - rly.reserve( (Lmax+1)*(Lmax+1) ); + rly.resize( (Lmax+1)*(Lmax+1) ); //begin calculation /***************************