From 2a3f1d53deb35918a1c1acc827d3f4cdd9c93719 Mon Sep 17 00:00:00 2001 From: RXL6105505822 Date: Wed, 18 Mar 2026 14:22:40 +0800 Subject: [PATCH 1/4] =?UTF-8?q?=E4=B8=BAsto=5Fwf.cpp=E4=BF=AE=E5=A4=8D?= =?UTF-8?q?=E4=BA=86=E9=AD=94=E6=95=B0=E6=9C=AA=E8=AF=B4=E6=98=8E=E5=90=AB?= =?UTF-8?q?=E4=B9=89=E4=BB=A5=E5=8F=8A=E4=BD=BF=E7=94=A8C=E9=A3=8E?= =?UTF-8?q?=E6=A0=BC=E5=BC=BA=E5=88=B6=E8=BD=AC=E6=8D=A2=E7=9A=84=E9=97=AE?= =?UTF-8?q?=E9=A2=98?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- source/source_pw/module_stodft/sto_wf.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/source_pw/module_stodft/sto_wf.cpp b/source/source_pw/module_stodft/sto_wf.cpp index a0204e1f87..11405384da 100644 --- a/source/source_pw/module_stodft/sto_wf.cpp +++ b/source/source_pw/module_stodft/sto_wf.cpp @@ -64,9 +64,10 @@ void Stochastic_WF::clean_chiallorder() template void Stochastic_WF::init_sto_orbitals(const int seed_in) { + const MAGIC_NUMBER = 10000; if (seed_in == 0 || seed_in == -1) { - srand((unsigned)time(nullptr) + GlobalV::MY_RANK * 10000); // GlobalV global variables are reserved + srand(static_cast(time(nullptr)) + GlobalV::MY_RANK * MAGIC_NUMBER); // GlobalV global variables are reserved } else { From 93774c549906c33c44199c748e92542c64275889 Mon Sep 17 00:00:00 2001 From: RXL6105505822 Date: Wed, 18 Mar 2026 22:36:44 +0800 Subject: [PATCH 2/4] =?UTF-8?q?=E4=BF=AE=E6=AD=A3=E4=BA=86=E4=B8=80?= =?UTF-8?q?=E5=A4=84=E5=B8=B8=E9=87=8F=E7=B1=BB=E5=9E=8B=E6=9C=AA=E5=AE=9A?= =?UTF-8?q?=E4=B9=89=E9=94=99=E8=AF=AF?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- source/source_pw/module_stodft/sto_wf.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/source_pw/module_stodft/sto_wf.cpp b/source/source_pw/module_stodft/sto_wf.cpp index 11405384da..9ca2403b4f 100644 --- a/source/source_pw/module_stodft/sto_wf.cpp +++ b/source/source_pw/module_stodft/sto_wf.cpp @@ -64,7 +64,7 @@ void Stochastic_WF::clean_chiallorder() template void Stochastic_WF::init_sto_orbitals(const int seed_in) { - const MAGIC_NUMBER = 10000; + const int MAGIC_NUMBER = 10000; if (seed_in == 0 || seed_in == -1) { srand(static_cast(time(nullptr)) + GlobalV::MY_RANK * MAGIC_NUMBER); // GlobalV global variables are reserved From 3e34da995c8823353e3f23276852297681d30e04 Mon Sep 17 00:00:00 2001 From: RXL6105505822 Date: Thu, 19 Mar 2026 21:08:16 +0800 Subject: [PATCH 3/4] 2026PKUCourseHW5: Handled case 10 by following the method in the PPT. --- source/source_base/math_ylmreal.cpp | 461 +++++++++++++++------------- 1 file changed, 250 insertions(+), 211 deletions(-) diff --git a/source/source_base/math_ylmreal.cpp b/source/source_base/math_ylmreal.cpp index e94b027bb3..752636be22 100644 --- a/source/source_base/math_ylmreal.cpp +++ b/source/source_base/math_ylmreal.cpp @@ -19,225 +19,264 @@ namespace ModuleBase YlmReal::YlmReal(){} YlmReal::~YlmReal(){} -void YlmReal::rlylm -( - const int lmax, - const double& x, - const double& y, - const double& z, // g_cartesian_vec(x,y,z) - double* rly // output -) -{ - ModuleBase::timer::tick("YlmReal","rlylm"); - - assert(lmax >= 0); - - //get xy_dependence - assert(lmax <= 19); - - double Am[20]; - double Bm[20]; +/** + * @brief Initialize arrays to zero + * @param Am Array for Am coefficients + * @param Bm Array for Bm coefficients + * @param size Size of arrays (should be at least lmax+1) + */ +void init_arrays(double* Am, double* Bm, const int size) { + for(int i = 0; i < size; ++i) { + Am[i] = 0.0; + Bm[i] = 0.0; + } +} - // mohan add 2021-05-07 - for(int i=0; i<20; ++i) - { - Am[i]=0.0; - Bm[i]=0.0; - } - - //ZEROS(Am, 20); - //ZEROS(Bm, 20); - - double x2, x3, x4, x5; - double y2, y3, y4, y5; - - x2 = x * x; - x3 = x2 * x; - x4 = x3 * x; - x5 = x4 * x; - - y2 = y * y; - y3 = y2 * y; - y4 = y3 * y; - y5 = y4 * y; - - //x-y dependence - //Am - //Bm - for(int im = 0; im < lmax+1; im++) - { - if(im == 0) - { - Am[0] = 1.0; - Bm[0] = 0.0; - } - else if(im == 1) - { - Am[1] = x; - Bm[1] = y; - } - else if(im == 2) - { - Am[2] = x2- y2; - Bm[2] = 2.0 * x * y; - } - else if(im == 3) - { - Am[3] = x3 - 3.0 * x * y2; - Bm[3] = 3.0 * x2 * y - y3; - } - else if(im == 4) - { - Am[4] = x4 - 6.0 * x2 * y2 + y4; - Bm[4] = 4.0 * (x3 * y - x * y3); - } - else if(im == 5) - { - Am[5] = x5 - 10.0 * x3 * y2 + 5.0 * x * y4; - Bm[5] = 5.0 * x4 * y - 10.0 * x2 * y3 + y5; - } - else - { - for(int ip = 0; ip <= im; ip++) - { - double aux = Fact(im) / Fact(ip) / Fact(im - ip); - Am[im] += aux * pow(x, ip) * pow(y, im-ip) * cos( (im-ip) * ModuleBase::PI / 2.0 ); - Bm[im] += aux * pow(x, ip) * pow(y, im-ip) * sin( (im-ip) * ModuleBase::PI / 2.0 ); - } - } - } - - //z dependence - double zdep[20][20]; - - for(int il = 0; il < 20; il++) - { - for(int jl=0; jl < 20; jl++) - { - zdep[il][jl]=0.0; // mohan add 2021-05-07 - } -// ZEROS(zdep[il], 20); - } +/** + * @brief Initialize zdep array to zero + * @param zdep 2D array for z-dependent coefficients + * @param size Size of zdep (should be at least lmax+1) + */ +void init_zdep(double zdep[][20], const int size) { + for(int il = 0; il < size; ++il) { + for(int jl = 0; jl < size; ++jl) { + zdep[il][jl] = 0.0; // mohan add 2021-05-07 + } + } +} - double z2 = z * z; - double z3 = z2 * z; - double z4 = z3 * z; - //double z5 = z4 * z; - - double r = sqrt(x*x + y*y + z*z); - double r2 = r * r; - double r3 = r2 * r; - double r4 = r3 * r; - - for(int il = 0; il < lmax + 1; il++) - { - if(il == 0) - { - zdep[0][0] = 1.0; - } - else if(il == 1) - { - zdep[1][0] = z; - zdep[1][1] = 1.0; - } - else if(il == 2) - { - zdep[2][0] = 0.5 * (3.0 * z2 - r2); - zdep[2][1] = sqrt(3.0) * z; - zdep[2][2] = sqrt(3.0) * 0.5; - } - else if(il == 3) - { - zdep[3][0] = 2.5 * z3 - 1.5 * z * r2; - zdep[3][1] = 0.25 * sqrt(6.0) * (5.0 * z2 - r2); - zdep[3][2] = 0.5 * sqrt(15.0) * z; - zdep[3][3] = 0.25 * sqrt(10.0); - } - else if(il == 4) - { - zdep[4][0] = 0.125 * (35.0 * z4 - 30.0 * r2 * z2 + 3.0 * r4); - zdep[4][1] = sqrt(10.0) * 0.25 * z * (7.0 * z2 - 3.0 * r2); - zdep[4][2] = sqrt(5.0) * 0.25 * (7.0 * z2 - r2); - zdep[4][3] = sqrt(70.0) * 0.25 * z; - zdep[4][4] = sqrt(35.0) * 0.125; - } - else if(il == 5) - { - zdep[5][0] = 0.125 * z *( 63.0 * z4 - 70.0 * z2 * r2 + 15.0 * r4); - zdep[5][1] = 0.125 * sqrt(15.0) * (21.0 * z4 - 14.0 * z2 * r2 + r4); - zdep[5][2] = 0.25 * sqrt(105.0) * z * (3.0 * z2 - r2); - zdep[5][3] = 0.0625 * sqrt(70.0) * (9.0 * z2 - r2); - zdep[5][4] = 0.375 * sqrt(35.0) * z; - zdep[5][5] = 0.1875 * sqrt(14.0); - } - else - { - for(int im = 0; im <= il; im++) - { - int kmax = static_cast( (il - im) / 2 ); - for(int ik = 0; ik <= kmax; ik++) - { - int twok = 2 * ik; - - double gamma = 0.0; - double aux0, aux1, aux2, aux3; - - aux0 = pow(-1.0, ik) * pow(2.0, -il); - aux1 = Fact(il) / Fact(ik) / Fact(il-ik); - aux2 = Fact(2*il - twok) / Fact(il) / Fact(il - twok); - aux3 = Fact(il - twok) / Fact(il - twok - im); - - gamma = aux0 * aux1 * aux2 * aux3; - - assert(il - twok - im >= 0); - zdep[il][im] += pow(r, twok) * pow(z, il-twok-im) * gamma; - } +/** + * @brief Calculate Am and Bm coefficients (xy-plane dependence) + * @param lmax Maximum angular momentum + * @param x x-coordinate + * @param y y-coordinate + * @param Am Output array for Am coefficients + * @param Bm Output array for Bm coefficients + * @param fact Factorial function (assumed to be available) + */ +void calculate_xy_coefficients(const int lmax, const double& x, const double& y, + double* Am, double* Bm) { + double x2 = x * x; + double x3 = x2 * x; + double x4 = x3 * x; + double x5 = x4 * x; + + double y2 = y * y; + double y3 = y2 * y; + double y4 = y3 * y; + double y5 = y4 * y; + + for(int im = 0; im < lmax + 1; im++) { + if(im == 0) { + Am[0] = 1.0; + Bm[0] = 0.0; + } + else if(im == 1) { + Am[1] = x; + Bm[1] = y; + } + else if(im == 2) { + Am[2] = x2 - y2; + Bm[2] = 2.0 * x * y; + } + else if(im == 3) { + Am[3] = x3 - 3.0 * x * y2; + Bm[3] = 3.0 * x2 * y - y3; + } + else if(im == 4) { + Am[4] = x4 - 6.0 * x2 * y2 + y4; + Bm[4] = 4.0 * (x3 * y - x * y3); + } + else if(im == 5) { + Am[5] = x5 - 10.0 * x3 * y2 + 5.0 * x * y4; + Bm[5] = 5.0 * x4 * y - 10.0 * x2 * y3 + y5; + } + else { + for(int ip = 0; ip <= im; ip++) { + double aux = Fact(im) / Fact(ip) / Fact(im - ip); + Am[im] += aux * pow(x, ip) * pow(y, im-ip) * cos((im-ip) * ModuleBase::PI / 2.0); + Bm[im] += aux * pow(x, ip) * pow(y, im-ip) * sin((im-ip) * ModuleBase::PI / 2.0); + } + } + } +} - if(im >= 1) - { - zdep[il][im] *= sqrt(2 * Fact(il - im) / Fact(il + im)); - - } - } - } - } +/** + * @brief Calculate low-order zdep coefficients (l <= 5) + * @param il Current angular momentum + * @param zdep Output array for z-dependent coefficients + * @param z z-coordinate + * @param r radius + */ +void calculate_low_order_zdep(const int il, double zdep[][20], + const double& z, const double& r) { + double z2 = z * z; + double z3 = z2 * z; + double z4 = z3 * z; + double r2 = r * r; + double r3 = r2 * r; + double r4 = r3 * r; + + if(il == 0) { + zdep[0][0] = 1.0; + } + else if(il == 1) { + zdep[1][0] = z; + zdep[1][1] = 1.0; + } + else if(il == 2) { + zdep[2][0] = 0.5 * (3.0 * z2 - r2); + zdep[2][1] = sqrt(3.0) * z; + zdep[2][2] = sqrt(3.0) * 0.5; + } + else if(il == 3) { + zdep[3][0] = 2.5 * z3 - 1.5 * z * r2; + zdep[3][1] = 0.25 * sqrt(6.0) * (5.0 * z2 - r2); + zdep[3][2] = 0.5 * sqrt(15.0) * z; + zdep[3][3] = 0.25 * sqrt(10.0); + } + else if(il == 4) { + zdep[4][0] = 0.125 * (35.0 * z4 - 30.0 * r2 * z2 + 3.0 * r4); + zdep[4][1] = sqrt(10.0) * 0.25 * z * (7.0 * z2 - 3.0 * r2); + zdep[4][2] = sqrt(5.0) * 0.25 * (7.0 * z2 - r2); + zdep[4][3] = sqrt(70.0) * 0.25 * z; + zdep[4][4] = sqrt(35.0) * 0.125; + } + else if(il == 5) { + zdep[5][0] = 0.125 * z *( 63.0 * z4 - 70.0 * z2 * r2 + 15.0 * r4); + zdep[5][1] = 0.125 * sqrt(15.0) * (21.0 * z4 - 14.0 * z2 * r2 + r4); + zdep[5][2] = 0.25 * sqrt(105.0) * z * (3.0 * z2 - r2); + zdep[5][3] = 0.0625 * sqrt(70.0) * (9.0 * z2 - r2); + zdep[5][4] = 0.375 * sqrt(35.0) * z; + zdep[5][5] = 0.1875 * sqrt(14.0); + } +} - //calc - int ic = 0; +/** + * @brief Calculate high-order zdep coefficients (l > 5) + * @param il Current angular momentum + * @param im Current magnetic quantum number + * @param zdep Output array for z-dependent coefficients + * @param r radius + * @param z z-coordinate + */ +void calculate_high_order_zdep(const int il, const int im, double zdep[][20], + const double& r, const double& z) { + int kmax = static_cast((il - im) / 2); + for(int ik = 0; ik <= kmax; ik++) { + int twok = 2 * ik; + + double gamma = 0.0; + double aux0, aux1, aux2, aux3; + + aux0 = pow(-1.0, ik) * pow(2.0, -il); + aux1 = Fact(il) / Fact(ik) / Fact(il-ik); + aux2 = Fact(2*il - twok) / Fact(il) / Fact(il - twok); + aux3 = Fact(il - twok) / Fact(il - twok - im); + + gamma = aux0 * aux1 * aux2 * aux3; + + assert(il - twok - im >= 0); + zdep[il][im] += pow(r, twok) * pow(z, il-twok-im) * gamma; + } - //special case for r=0 - double rpi = r; - const double tiny = 1.0E-10; - if (rpi < tiny) rpi += tiny; - - for(int il = 0; il <= lmax; il++) - { - double fac = sqrt( (2.0 * il + 1.0) / ModuleBase::FOUR_PI ); + if(im >= 1) { + zdep[il][im] *= sqrt(2 * Fact(il - im) / Fact(il + im)); + } +} - double rl = pow(rpi, il); - - //m=0 - rly[ic] = Am[0] * zdep[il][0] * fac / rl; - - ic++; - - //m ! = 0 - for(int im = 1; im <= il; im++) - { - //m>0 - rly[ic] = Am[im] * zdep[il][im] * pow(-1.0, im) * fac / rl; - - ic++; - - //m<0 - rly[ic] = Bm[im] * zdep[il][im] * pow(-1.0, im) * fac / rl; +/** + * @brief Calculate zdep coefficients (z dependence) + * @param lmax Maximum angular momentum + * @param z z-coordinate + * @param r radius + * @param zdep Output 2D array for z-dependent coefficients + */ +void calculate_z_coefficients(const int lmax, const double& z, const double& r, + double zdep[][20]) { + for(int il = 0; il < lmax + 1; il++) { + if(il <= 5) { + calculate_low_order_zdep(il, zdep, z, r); + } + else { + for(int im = 0; im <= il; im++) { + calculate_high_order_zdep(il, im, zdep, r, z); + } + } + } +} - ic++; - } - } +/** + * @brief Combine coefficients to compute real spherical harmonics + * @param lmax Maximum angular momentum + * @param Am Am coefficients from xy-plane + * @param Bm Bm coefficients from xy-plane + * @param zdep z-dependent coefficients + * @param r radius + * @param rly Output array for real spherical harmonics + */ +void combine_coefficients(const int lmax, const double* Am, const double* Bm, + double zdep[][20], const double& r, double* rly) { + int ic = 0; + double rpi = r; + const double tiny = 1.0E-10; + if (rpi < tiny) rpi += tiny; + + for(int il = 0; il <= lmax; il++) { + double fac = sqrt((2.0 * il + 1.0) / ModuleBase::FOUR_PI); + double rl = pow(rpi, il); + + // m=0 + rly[ic] = Am[0] * zdep[il][0] * fac / rl; + ic++; + + // m != 0 + for(int im = 1; im <= il; im++) { + // m>0 + rly[ic] = Am[im] * zdep[il][im] * pow(-1.0, im) * fac / rl; + ic++; + + // m<0 + rly[ic] = Bm[im] * zdep[il][im] * pow(-1.0, im) * fac / rl; + ic++; + } + } +} - ModuleBase::timer::tick("YlmReal","rlylm"); - return; +/** + * @brief Compute real spherical harmonics Ylm for given coordinates + * @param lmax Maximum angular momentum + * @param x x-coordinate in Cartesian + * @param y y-coordinate in Cartesian + * @param z z-coordinate in Cartesian + * @param rly Output array for real spherical harmonics + */ +void YlmReal::rlylm(const int lmax, const double& x, const double& y, const double& z, + double* rly) { + ModuleBase::timer::tick("YlmReal", "rlylm"); + + assert(lmax >= 0); + assert(lmax <= 19); // mohan add 2021-05-07 + + double Am[20]; + double Bm[20]; + double zdep[20][20]; + + // Initialize arrays + init_arrays(Am, Bm, 20); + init_zdep(zdep, 20); + + // Calculate coefficients + calculate_xy_coefficients(lmax, x, y, Am, Bm); + + double r = sqrt(x*x + y*y + z*z); + calculate_z_coefficients(lmax, z, r, zdep); + + // Combine and output + combine_coefficients(lmax, Am, Bm, zdep, r, rly); + + ModuleBase::timer::tick("YlmReal", "rlylm"); + return; } From 6ce10275fe9fad1014c87a9d149593a426390ea8 Mon Sep 17 00:00:00 2001 From: RXL6105505822 Date: Fri, 20 Mar 2026 10:47:32 +0800 Subject: [PATCH 4/4] Fixed this problem: Fact was not declared in this scope. --- source/source_base/math_ylmreal.cpp | 14 +++++++------- source/source_base/math_ylmreal.h | 13 +++++++++++++ 2 files changed, 20 insertions(+), 7 deletions(-) diff --git a/source/source_base/math_ylmreal.cpp b/source/source_base/math_ylmreal.cpp index 752636be22..1c6260de44 100644 --- a/source/source_base/math_ylmreal.cpp +++ b/source/source_base/math_ylmreal.cpp @@ -25,7 +25,7 @@ YlmReal::~YlmReal(){} * @param Bm Array for Bm coefficients * @param size Size of arrays (should be at least lmax+1) */ -void init_arrays(double* Am, double* Bm, const int size) { +void YlmReal::init_arrays(double* Am, double* Bm, const int size) { for(int i = 0; i < size; ++i) { Am[i] = 0.0; Bm[i] = 0.0; @@ -37,7 +37,7 @@ void init_arrays(double* Am, double* Bm, const int size) { * @param zdep 2D array for z-dependent coefficients * @param size Size of zdep (should be at least lmax+1) */ -void init_zdep(double zdep[][20], const int size) { +void YlmReal::init_zdep(double zdep[][20], const int size) { for(int il = 0; il < size; ++il) { for(int jl = 0; jl < size; ++jl) { zdep[il][jl] = 0.0; // mohan add 2021-05-07 @@ -54,7 +54,7 @@ void init_zdep(double zdep[][20], const int size) { * @param Bm Output array for Bm coefficients * @param fact Factorial function (assumed to be available) */ -void calculate_xy_coefficients(const int lmax, const double& x, const double& y, +void YlmReal::calculate_xy_coefficients(const int lmax, const double& x, const double& y, double* Am, double* Bm) { double x2 = x * x; double x3 = x2 * x; @@ -108,7 +108,7 @@ void calculate_xy_coefficients(const int lmax, const double& x, const double& y, * @param z z-coordinate * @param r radius */ -void calculate_low_order_zdep(const int il, double zdep[][20], +void YlmReal::calculate_low_order_zdep(const int il, double zdep[][20], const double& z, const double& r) { double z2 = z * z; double z3 = z2 * z; @@ -160,7 +160,7 @@ void calculate_low_order_zdep(const int il, double zdep[][20], * @param r radius * @param z z-coordinate */ -void calculate_high_order_zdep(const int il, const int im, double zdep[][20], +void YlmReal::calculate_high_order_zdep(const int il, const int im, double zdep[][20], const double& r, const double& z) { int kmax = static_cast((il - im) / 2); for(int ik = 0; ik <= kmax; ik++) { @@ -192,7 +192,7 @@ void calculate_high_order_zdep(const int il, const int im, double zdep[][20], * @param r radius * @param zdep Output 2D array for z-dependent coefficients */ -void calculate_z_coefficients(const int lmax, const double& z, const double& r, +void YlmReal::calculate_z_coefficients(const int lmax, const double& z, const double& r, double zdep[][20]) { for(int il = 0; il < lmax + 1; il++) { if(il <= 5) { @@ -215,7 +215,7 @@ void calculate_z_coefficients(const int lmax, const double& z, const double& r, * @param r radius * @param rly Output array for real spherical harmonics */ -void combine_coefficients(const int lmax, const double* Am, const double* Bm, +void YlmReal::combine_coefficients(const int lmax, const double* Am, const double* Bm, double zdep[][20], const double& r, double* rly) { int ic = 0; double rpi = r; diff --git a/source/source_base/math_ylmreal.h b/source/source_base/math_ylmreal.h index 9643196a21..87be9466e4 100644 --- a/source/source_base/math_ylmreal.h +++ b/source/source_base/math_ylmreal.h @@ -99,6 +99,19 @@ class YlmReal static long double Fact(const int n); static int Semi_Fact(const int n); + //auxiliary functions for rlylm + static void init_arrays(double* Am, double* Bm, const int size); + static void init_zdep(double zdep[][20], const int size); + static void calculate_xy_coefficients(const int lmax, const double& x, const double& y, + double* Am, double* Bm); + static void calculate_low_order_zdep(const int il, double zdep[][20], + const double& z, const double& r); + static void calculate_high_order_zdep(const int il, const int im, double zdep[][20], + const double& r, const double& z); + static void calculate_z_coefficients(const int lmax, const double& z, const double& r, + double zdep[][20]); + static void combine_coefficients(const int lmax, const double* Am, const double* Bm, + double zdep[][20], const double& r, double* rly); };