Skip to content

Commit

Permalink
Initialized some variables in projgen.cpp (#3905)
Browse files Browse the repository at this point in the history
  • Loading branch information
@Haerxile
Haerxile authored Apr 3, 2024
1 parent e49fce2 commit b62a807
Showing 1 changed file with 58 additions and 35 deletions.
93 changes: 58 additions & 35 deletions source/module_basis/module_nao/projgen.cpp
View file
Original file line number Diff line number Diff line change
@@ -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 <numeric>
#include <cassert>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdio>
#include <numeric>

#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<double>& alpha)
void projgen(const int l,
const int nr,
const double* r,
const double* chi,
const double rcut,
const int nbes,
std::vector<double>& alpha)
{
assert(rcut < r[nr - 1]);
assert(std::is_sorted(r, r + nr));
Expand All @@ -30,8 +37,9 @@ void projgen(const int l, const int nr, const double* r, const double* chi, cons
std::vector<double> z(nbes);
std::vector<double> 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<double> tmp(nr_proj);
Expand All @@ -40,30 +48,35 @@ 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<double> 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]);
}

// optimal coefficients
std::vector<double> 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<double>& alpha)
void smoothgen(const int nr, const double* r, const double* chi, const double rcut, std::vector<double>& alpha)
{
// lambda function for generate the new radial function
assert(rcut < r[nr - 1]);
Expand All @@ -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<nr_proj;i++)
for (int i = 0; i < nr_proj; i++)
{
alpha[i] = chi[i] * (1 - std::exp(- std::pow((r[i] - rcut), 2)/2/sigma_in/sigma_in));
alpha[i] = chi[i] * (1 - std::exp(-std::pow((r[i] - rcut), 2) / 2 / sigma_in / sigma_in));
}
// r^2 * chi (independent from p)
std::vector<double> 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<double> 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<nr_proj;i++)
for (int i = 0; i < nr_proj; i++)
{
alpha[i] /= std::sqrt(overlap);
}
Expand All @@ -105,21 +123,27 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc

// function for calculating the overlap between dalpha and dchi
auto overlap_dalpha_dchi = [&]() {
// calculate dalpha first
// calculate dalpha first
ModuleBase::CubicSpline cubspl_alpha;
cubspl_alpha.build(nr_proj, r, alpha.data());
std::vector<double> dalpha(nr_proj);
cubspl_alpha.eval(nr_proj, r, nullptr, dalpha.data());
for(int i=0;i<nr_proj;i++) dalpha[i] -= dchi[i];
for (int i = 0; i < nr_proj; i++)
dalpha[i] -= dchi[i];
// r^2 * dchi (independent from p)
std::vector<double> 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<double> 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]);
};

Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -171,4 +195,3 @@ void smoothgen(const int nr, const double* r, const double* chi, const double rc
}
}
}

0 comments on commit b62a807

Please sign in to comment.