From d6b79e3990d0a224a0765fbccacc30e35cefb88d Mon Sep 17 00:00:00 2001 From: Haerxile Date: Wed, 3 Apr 2024 02:34:26 +0800 Subject: [PATCH] Initialized some variables in projgen.cpp --- source/module_basis/module_nao/projgen.cpp | 93 ++++++++++++++-------- 1 file changed, 58 insertions(+), 35 deletions(-) diff --git a/source/module_basis/module_nao/projgen.cpp b/source/module_basis/module_nao/projgen.cpp index 4f32d6b461e..f55d130002e 100644 --- a/source/module_basis/module_nao/projgen.cpp +++ b/source/module_basis/module_nao/projgen.cpp @@ -1,17 +1,24 @@ -#include "module_base/math_integral.h" -#include "module_base/math_sphbes.h" -#include "module_base/cubic_spline.h" #include "projgen.h" -#include -#include #include +#include #include #include +#include + +#include "module_base/cubic_spline.h" +#include "module_base/math_integral.h" +#include "module_base/math_sphbes.h" using namespace ModuleBase; -void projgen(const int l, const int nr, const double* r, const double* chi, const double rcut, const int nbes, std::vector& alpha) +void projgen(const int l, + const int nr, + const double* r, + const double* chi, + const double rcut, + const int nbes, + std::vector& alpha) { assert(rcut < r[nr - 1]); assert(std::is_sorted(r, r + nr)); @@ -30,8 +37,9 @@ void projgen(const int l, const int nr, const double* r, const double* chi, cons std::vector z(nbes); std::vector w(nbes); - std::transform(theta.begin(), theta.end(), z.begin(), [rcut, l](double theta_p) - { return 0.5 * std::pow(rcut, 3) * std::pow(Sphbes::sphbesj(l+1, theta_p), 2); }); + std::transform(theta.begin(), theta.end(), z.begin(), [rcut, l](double theta_p) { + return 0.5 * std::pow(rcut, 3) * std::pow(Sphbes::sphbesj(l + 1, theta_p), 2); + }); // r^2 * chi (independent from p) std::vector tmp(nr_proj); @@ -40,9 +48,11 @@ void projgen(const int l, const int nr, const double* r, const double* chi, cons // r^2 * chi * j_l(theta[p] * r / rcut) (dependent on p) std::vector integrand(nr_proj); - for (int p = 0; p < nbes; ++p) { - std::transform(r, r + nr_proj, tmp.begin(), integrand.begin(), [theta, p, rcut, l](double r_i, double tmp_i) - { return tmp_i * Sphbes::sphbesj(l, theta[p] * r_i / rcut); }); + for (int p = 0; p < nbes; ++p) + { + std::transform(r, r + nr_proj, tmp.begin(), integrand.begin(), [theta, p, rcut, l](double r_i, double tmp_i) { + return tmp_i * Sphbes::sphbesj(l, theta[p] * r_i / rcut); + }); w[p] = Integral::simpson(nr_proj, integrand.data(), &dr[1]); } @@ -50,20 +60,23 @@ void projgen(const int l, const int nr, const double* r, const double* chi, cons std::vector c(nbes, 0.0); std::transform(w.begin(), w.end(), z.begin(), c.begin(), [](double w_p, double z_p) { return w_p * w_p / z_p; }); double prefac = 1.0 / std::sqrt(std::accumulate(c.begin(), c.end(), 0.0)); - std::transform(w.begin(), w.end(), z.begin(), c.begin(), [prefac](double w_p, double z_p) - { return prefac * w_p / z_p; }); + std::transform(w.begin(), w.end(), z.begin(), c.begin(), [prefac](double w_p, double z_p) { + return prefac * w_p / z_p; + }); // new radial function alpha.resize(nr_proj); std::fill(alpha.begin(), alpha.end(), 0.0); - for (int i = 0; i < nr_proj; ++i) { - for (int p = 0; p < nbes; ++p) { + for (int i = 0; i < nr_proj; ++i) + { + for (int p = 0; p < nbes; ++p) + { alpha[i] += c[p] * Sphbes::sphbesj(l, theta[p] * r[i] / rcut); } } } -void smoothgen(const int nr, const double* r, const double* chi, const double rcut, std::vector& alpha) +void smoothgen(const int nr, const double* r, const double* chi, const double rcut, std::vector& alpha) { // lambda function for generate the new radial function assert(rcut < r[nr - 1]); @@ -76,21 +89,26 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc int nr_proj = std::distance(r, std::lower_bound(r, r + nr, rcut)) + 1; alpha.resize(nr_proj); auto smooth_sigma = [&](double sigma_in) { - for(int i=0;i tmp(nr_proj); - std::transform(r, r + nr_proj, alpha.data(), tmp.begin(), [](double r_i, double chi_i) { return r_i * r_i * chi_i; }); + std::transform(r, r + nr_proj, alpha.data(), tmp.begin(), [](double r_i, double chi_i) { + return r_i * r_i * chi_i; + }); // r^2 * chi * chi std::vector integrand(nr_proj); - std::transform(alpha.data(), alpha.data() + nr_proj, tmp.begin(), integrand.begin(), [](double chi_i, double tmp_i) - { return tmp_i * chi_i; }); + std::transform(alpha.data(), + alpha.data() + nr_proj, + tmp.begin(), + integrand.begin(), + [](double chi_i, double tmp_i) { return tmp_i * chi_i; }); double overlap = ModuleBase::Integral::simpson(nr_proj, integrand.data(), &dr[1]); - for(int i=0;i dalpha(nr_proj); cubspl_alpha.eval(nr_proj, r, nullptr, dalpha.data()); - for(int i=0;i tmp(nr_proj); - std::transform(r, r + nr_proj, dalpha.data(), tmp.begin(), [](double r_i, double dalpha_i) { return r_i * r_i * dalpha_i; }); + std::transform(r, r + nr_proj, dalpha.data(), tmp.begin(), [](double r_i, double dalpha_i) { + return r_i * r_i * dalpha_i; + }); // r^2 * dalpha * dchi std::vector integrand(nr_proj); - std::transform(dalpha.data(), dalpha.data() + nr_proj, tmp.begin(), integrand.begin(), [](double dalpha_i, double tmp_i) - { return tmp_i * dalpha_i; }); + std::transform(dalpha.data(), + dalpha.data() + nr_proj, + tmp.begin(), + integrand.begin(), + [](double dalpha_i, double tmp_i) { return tmp_i * dalpha_i; }); return ModuleBase::Integral::simpson(nr_proj, integrand.data(), &dr[1]); }; @@ -130,16 +154,16 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc double sigma_right = 1.0; smooth_sigma(sigma_right); double overlap_alpha_chi_right = overlap_dalpha_dchi(); - double overlap_alpha_chi; - double sigma; + double overlap_alpha_chi = 0.0; + double sigma = 0.0; while (std::abs(overlap_alpha_chi_right - overlap_alpha_chi_left) > 1e-6) { sigma = (sigma_left + sigma_right) / 2; smooth_sigma(sigma); overlap_alpha_chi = overlap_dalpha_dchi(); - if(overlap_alpha_chi < overlap_alpha_chi_left && overlap_alpha_chi < overlap_alpha_chi_right) - {// the minimum is in the middle - if(overlap_alpha_chi_left > overlap_alpha_chi_right) + if (overlap_alpha_chi < overlap_alpha_chi_left && overlap_alpha_chi < overlap_alpha_chi_right) + { // the minimum is in the middle + if (overlap_alpha_chi_left > overlap_alpha_chi_right) { sigma_left = sigma; overlap_alpha_chi_left = overlap_alpha_chi; @@ -151,8 +175,8 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc } } else - {// the minimum is on the left or right - if(overlap_alpha_chi_left < overlap_alpha_chi_right) + { // the minimum is on the left or right + if (overlap_alpha_chi_left < overlap_alpha_chi_right) { sigma_right = sigma; overlap_alpha_chi_right = overlap_alpha_chi; @@ -171,4 +195,3 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc } } } -