Skip to content
GitLab
Commit 1d9df5f7 authored by Carl Poelking's avatar Carl Poelking
Browse files

Radial gradient numerical integration and test cases.

parent 2d279396
Hide whitespace changes
Inline Side-by-side
cxx_tests/src/gtest_cutoff.cpp
View file @ 1d9df5f7
......@@ -10,6 +10,8 @@ class TestCutoffShiftedCosine : public ::testing::Test
public:
soap::Options _options;
soap::CutoffFunction *_cutoff;
std::string _constructor_stdout;
std::string _constructor_stdout_ref;
virtual void SetUp() {
_options.set("radialcutoff.type", "shifted-cosine");
......@@ -19,7 +21,10 @@ public:
soap::CutoffFunctionFactory::registerAll();
_cutoff = soap::CutoffFunctionOutlet().create(_options.get<std::string>("radialcutoff.type"));
::testing::internal::CaptureStdout();
_cutoff->configure(_options);
_constructor_stdout = ::testing::internal::GetCapturedStdout();
_constructor_stdout_ref = "Weighting function with Rc = 4, _Rc_width = 0.5, central weight = 1\n";
}
virtual void TearDown() {
......@@ -33,6 +38,8 @@ TEST_F(TestCutoffShiftedCosine, Constructor) {
EXPECT_DOUBLE_EQ(_cutoff->getCenterWeight(), 1.);
EXPECT_DOUBLE_EQ(_cutoff->getCutoffWidth(), 0.5);
EXPECT_DOUBLE_EQ(_cutoff->getCutoff(), 4.);
EXPECT_EQ(_constructor_stdout, _constructor_stdout_ref);
}
TEST_F(TestCutoffShiftedCosine, Weight) {
......
cxx_tests/src/gtest_functions.cpp 0 → 100644
View file @ 1d9df5f7
#include <iostream>
#include <vector>
#include <gtest/gtest.h>
#include <soap/functions.hpp>
#include <boost/format.hpp>
TEST(TestFunctions, ModifiedSphericalBessel1stKind) {
int L = 10;
double dr = 0.05;
int N = 100;
soap::ModifiedSphericalBessel1stKind mosbest(L);
mosbest.evaluate(0.0, true);
EXPECT_DOUBLE_EQ(mosbest._in[0], 1.);
EXPECT_DOUBLE_EQ(mosbest._in[1], 0.);
EXPECT_DOUBLE_EQ(mosbest._in[10], 0.);
EXPECT_DOUBLE_EQ(mosbest._din[0], 0.);
EXPECT_DOUBLE_EQ(mosbest._din[1], 1./3.);
mosbest.evaluate(1., true);
EXPECT_NEAR(mosbest._in[0], +1.1752012e+00, 1e-4);
EXPECT_NEAR(mosbest._in[1], +3.6787944e-01, 1e-4);
EXPECT_NEAR(mosbest._in[9], +1.5641127e-09, 1e-4);
EXPECT_NEAR(mosbest._din[0], +3.6787944e-01, 1e-4);
EXPECT_NEAR(mosbest._din[7], +3.5860451e-06, 1e-4);
/*
std::vector<double> r_list;
r_list.push_back(0.);
r_list.push_back(0.1);
r_list.push_back(1.);
for (int i = 0; i < r_list.size(); ++i) {
mosbest.evaluate(r_list[i], true);
std::cout << "r = " << r_list[i] << std::endl;
for (int n = 0; n <= 10; ++n) {
std::cout << boost::format("%1$2d %2$+1.7e %3$+1.7e") % n % mosbest._in[n] % mosbest._din[n] << std::endl;
}
}
*/
/*
for (int i = 0; i < N; ++i) {
double ri = i*dr;
// Compute
mosbest.evaluate(ri, true);
// Output
std::cout << boost::format("%1$+1.7e ") % ri;
for (int l = 0; l <= L; ++l) {
std::cout << boost::format("l=%1$d %2$+1.7e %3$+1.7e ") % l % mosbest._in[l] % mosbest._din[l];
}
std::cout << std::endl;
}
*/
}
cxx_tests/src/gtest_radialbasis.cpp 0 → 100644
View file @ 1d9df5f7
#include <iostream>
#include <vector>
#include <boost/format.hpp>
#include <boost/algorithm/string.hpp>
#include <gtest/gtest.h>
#include <soap/radialbasis.hpp>
#include <soap/options.hpp>
static const double EPS = std::numeric_limits<double>::min();
#define EXPECT_NEAR_RELATIVE(A, B, T) ASSERT_NEAR(A/(fabs(A)+fabs(B) + EPS), B/(fabs(A)+fabs(B) + EPS), T)
class TestRadialBasisGaussian : public ::testing::Test
{
public:
soap::Options _options;
soap::RadialBasis *_radbasis;
std::string _constructor_stdout;
std::string _constructor_stdout_ref;
virtual void SetUp() {
_options.set("radialcutoff.Rc", 4.);
_options.set("radialbasis.type", "gaussian");
_options.set("radialbasis.mode", "adaptive");
_options.set("radialbasis.N", 9);
_options.set("angularbasis.L", 6);
_options.set("radialbasis.sigma", 0.5);
_options.set("radialbasis.integration_steps", 15);
soap::RadialBasisFactory::registerAll();
_radbasis = soap::RadialBasisOutlet().create(_options.get<std::string>("radialbasis.type"));
::testing::internal::CaptureStdout();
_radbasis->configure(_options);
_constructor_stdout = ::testing::internal::GetCapturedStdout();
_constructor_stdout_ref =
"Adjusted radial cutoff to 3.41181 based on sigma_0 = 0.5, L = 6, stride = 0.5\n"
"Created 9 radial Gaussians at:\n"
" r = +0.0000000e+00 (sigma = +5.0000000e-01) \n"
" r = +2.5000000e-01 (sigma = +5.1887452e-01) \n"
" r = +5.0943726e-01 (sigma = +5.7432939e-01) \n"
" r = +7.9660196e-01 (sigma = +6.6727337e-01) \n"
" r = +1.1302386e+00 (sigma = +8.0190914e-01) \n"
" r = +1.5311932e+00 (sigma = +9.8559668e-01) \n"
" r = +2.0239916e+00 (sigma = +1.2290136e+00) \n"
" r = +2.6384983e+00 (sigma = +1.5466265e+00) \n"
" r = +3.4118116e+00 (sigma = +1.9574676e+00) \n"
"Radial basis overlap matrix\n"
"+1.0000e+00 +9.5931e-01 +8.1266e-01 +5.9142e-01 +3.7675e-01 +2.2009e-01 +1.2395e-01 +6.9855e-02 +4.0214e-02 \n"
"+9.5931e-01 +1.0000e+00 +9.3655e-01 +7.5104e-01 +5.1642e-01 +3.1649e-01 +1.8212e-01 +1.0292e-01 +5.8776e-02 \n"
"+8.1266e-01 +9.3655e-01 +1.0000e+00 +9.2139e-01 +7.1532e-01 +4.7799e-01 +2.8856e-01 +1.6580e-01 +9.4361e-02 \n"
"+5.9142e-01 +7.5104e-01 +9.2139e-01 +1.0000e+00 +9.1258e-01 +6.9606e-01 +4.5873e-01 +2.7546e-01 +1.5859e-01 \n"
"+3.7675e-01 +5.1642e-01 +7.1532e-01 +9.1258e-01 +1.0000e+00 +9.0760e-01 +6.8560e-01 +4.4877e-01 +2.6901e-01 \n"
"+2.2009e-01 +3.1649e-01 +4.7799e-01 +6.9606e-01 +9.0760e-01 +1.0000e+00 +9.0474e-01 +6.7975e-01 +4.4337e-01 \n"
"+1.2395e-01 +1.8212e-01 +2.8856e-01 +4.5873e-01 +6.8560e-01 +9.0474e-01 +1.0000e+00 +9.0307e-01 +6.7638e-01 \n"
"+6.9855e-02 +1.0292e-01 +1.6580e-01 +2.7546e-01 +4.4877e-01 +6.7975e-01 +9.0307e-01 +1.0000e+00 +9.0207e-01 \n"
"+4.0214e-02 +5.8776e-02 +9.4361e-02 +1.5859e-01 +2.6901e-01 +4.4337e-01 +6.7638e-01 +9.0207e-01 +1.0000e+00 \n"
"Radial basis Cholesky decomposition\n"
"+1.0000e+00 +9.5931e-01 +8.1266e-01 +5.9142e-01 +3.7675e-01 +2.2009e-01 +1.2395e-01 +6.9855e-02 +4.0214e-02 \n"
"+9.5931e-01 +2.8236e-01 +5.5590e-01 +6.5053e-01 +5.4894e-01 +3.7312e-01 +2.2389e-01 +1.2717e-01 +7.1531e-02 \n"
"+8.1266e-01 +5.5590e-01 +1.7483e-01 +4.5265e-01 +5.9484e-01 +5.2457e-01 +3.6246e-01 +2.1931e-01 +1.2535e-01 \n"
"+5.9142e-01 +6.5053e-01 +4.5265e-01 +1.4881e-01 +4.2606e-01 +5.7602e-01 +5.0881e-01 +3.5048e-01 +2.1189e-01 \n"
"+3.7675e-01 +5.4894e-01 +5.9484e-01 +4.2606e-01 +1.4613e-01 +4.2698e-01 +5.7217e-01 +4.9862e-01 +3.4046e-01 \n"
"+2.2009e-01 +3.7312e-01 +5.2457e-01 +5.7602e-01 +4.2698e-01 +1.5183e-01 +4.3733e-01 +5.7369e-01 +4.9171e-01 \n"
"+1.2395e-01 +2.2389e-01 +3.6246e-01 +5.0881e-01 +5.7217e-01 +4.3733e-01 +1.6003e-01 +4.4949e-01 +5.7673e-01 \n"
"+6.9855e-02 +1.2717e-01 +2.1931e-01 +3.5048e-01 +4.9862e-01 +5.7369e-01 +4.4949e-01 +1.6801e-01 +4.6027e-01 \n"
"+4.0214e-02 +7.1531e-02 +1.2535e-01 +2.1189e-01 +3.4046e-01 +4.9171e-01 +5.7673e-01 +4.6027e-01 +1.7464e-01 \n"
"Radial basis transformation matrix\n"
"+1.0000e+00 +0.0000e+00 +0.0000e+00 +0.0000e+00 +0.0000e+00 +0.0000e+00 +0.0000e+00 +0.0000e+00 +0.0000e+00 \n"
"-3.3974e+00 +3.5415e+00 -7.5092e-15 -3.4133e-16 -5.3098e-16 +3.3945e-16 +9.0948e-31 -1.7053e-31 +1.3263e-31 \n"
"+6.1541e+00 -1.1260e+01 +5.7197e+00 +6.7046e-15 -1.6762e-15 +0.0000e+00 -1.0755e-15 +1.2945e-16 -1.8179e-16 \n"
"-7.8418e+00 +1.8770e+01 -1.7398e+01 +6.7199e+00 +7.4196e-15 -1.4839e-15 +3.2713e-15 -3.9375e-16 +5.5296e-16 \n"
"+7.9966e+00 -2.2192e+01 +2.7442e+01 -1.9592e+01 +6.8431e+00 +4.3519e-15 -6.5278e-15 +1.0880e-15 -8.7222e-16 \n"
"-7.1003e+00 +2.1400e+01 -3.0931e+01 +2.9604e+01 -1.9245e+01 +6.5864e+00 +9.6627e-15 -2.0131e-15 +9.8310e-16 \n"
"+5.7853e+00 -1.8265e+01 +2.8772e+01 -3.2217e+01 +2.8125e+01 -1.7999e+01 +6.2488e+00 +3.3970e-15 -8.4926e-16 \n"
"-4.4840e+00 +1.4517e+01 -2.3973e+01 +2.9234e+01 -2.9840e+01 +2.5665e+01 -1.6718e+01 +5.9520e+00 +5.8750e-16 \n"
"+3.3731e+00 -1.1074e+01 +1.8762e+01 -2.3966e+01 +2.6611e+01 -2.6745e+01 +2.3425e+01 -1.5687e+01 +5.7261e+00 \n";
}
virtual void TearDown() {
delete _radbasis;
_radbasis = NULL;
}
};
TEST_F(TestRadialBasisGaussian, Constructor) {
EXPECT_EQ(_radbasis->identify(), "gaussian");
EXPECT_EQ(_radbasis->N(), 9);
// GENERATE REFERENCE
/*
std::vector<std::string> fields_ref;
boost::split(fields_ref, _constructor_stdout, boost::is_any_of("\n"));
for (auto f = fields_ref.begin(); f != fields_ref.end(); ++f) {
std::cout << " \"" << (*f) << "\\n\"" << std::endl;
}
*/
// COMPARE CONSTRUCTOR STDOUT
std::vector<std::string> fields;
std::vector<std::string> fields_ref;
boost::replace_all(_constructor_stdout, "\n", " ");
boost::replace_all(_constructor_stdout_ref, "\n", " ");
boost::replace_all(_constructor_stdout, " ", " ");
boost::replace_all(_constructor_stdout_ref, " ", " ");
boost::split(fields, _constructor_stdout, boost::is_any_of(" "));
boost::split(fields_ref, _constructor_stdout_ref, boost::is_any_of(" "));
auto f = fields.begin();
auto fref = fields_ref.begin();
for ( ; f != fields.end() && fref != fields_ref.end(); ++f, ++fref) {
EXPECT_EQ((*f), (*fref));
}
}
TEST_F(TestRadialBasisGaussian, computeCoefficientsGradients) {
int N = _radbasis->N();
int L = _options.get<int>("angularbasis.L");
soap::RadialBasis::radcoeff_t Gnl = soap::RadialBasis::radcoeff_zero_t(N,L+1);
soap::RadialBasis::radcoeff_t dGnl_dx = soap::RadialBasis::radcoeff_zero_t(N,L+1);
soap::RadialBasis::radcoeff_t dGnl_dy = soap::RadialBasis::radcoeff_zero_t(N,L+1);
soap::RadialBasis::radcoeff_t dGnl_dz = soap::RadialBasis::radcoeff_zero_t(N,L+1);
/*
double sigma = 0.5;
soap::vec d(0.,1.,0.);
double r = 0.;
_radbasis->computeCoefficients(d, r, sigma, Gnl, dGnl_dx, dGnl_dy, dGnl_dz);
*/
}
TEST_F(TestRadialBasisGaussian, NumericalIntegrationBessel) {
// CHECK NUMERICAL INTEGRATION
double ai = 2.;
double ak = 2.;
double beta_ik = ai+ak;
double ri = 1.;
double rk = 1.;
double rho_ik = ak*rk/beta_ik;
int L_plus_1 = 7;
int integration_steps = 15;
/*
for (int n = 0; n < 100; ++n) {
ri = n*0.05;
std::vector<double> integrals;
integrals.resize(L_plus_1, 0.);
std::vector<double> integrals_derivative;
integrals_derivative.resize(L_plus_1, 0.);
soap::compute_integrals_il_expik_r2_dr(
ai, ri, beta_ik, rho_ik, L_plus_1, integration_steps,
&integrals, &integrals_derivative);
std::cout << boost::format("%1$+1.7e ") % ri;
for (int l = 0; l < L_plus_1; ++l) {
std::cout << boost::format("l=%1$d %2$+1.7e %3$+1.7e ") % l % integrals[l] % integrals_derivative[l];
}
std::cout << std::endl;
}
*/
std::vector<double> r_list{0., 0.5, 3.};
std::vector<double> samples;
std::vector<double> samples_ref{
+3.2927839e-01, +0.0000000e+00, +0.0000000e+00,
+5.5035011e-01, +2.7537077e-02, +2.6313175e-04,
+2.3985116e+05, +1.8595085e+05, +9.9749558e+04};
std::vector<double> samples_d;
std::vector<double> samples_d_ref{
+3.7086702e-01, +0.0000000e+00, +0.0000000e+00,
+9.6300224e-01, +5.1400603e-01, +6.0036046e-02,
+1.8193830e+06, +1.6815557e+06, +1.2788038e+06};
for (auto &ri : r_list) {
std::vector<double> integrals;
integrals.resize(L_plus_1, 0.);
std::vector<double> integrals_derivative;
integrals_derivative.resize(L_plus_1, 0.);
soap::compute_integrals_il_expik_r2_dr(
ai, ri, beta_ik, rho_ik, L_plus_1, integration_steps,
&integrals, &integrals_derivative);
samples.push_back(integrals[0]);
samples.push_back(integrals[3]);
samples.push_back(integrals[6]);
samples_d.push_back(integrals_derivative[1]);
samples_d.push_back(integrals_derivative[2]);
samples_d.push_back(integrals_derivative[4]);
}
for (int i = 0; i < samples.size(); ++i) {
EXPECT_NEAR_RELATIVE(samples[i], samples_ref[i], 1e-7);
}
for (int i = 0; i < samples_d.size(); ++i) {
EXPECT_NEAR_RELATIVE(samples_d[i], samples_d_ref[i], 1e-7);
}
}
src/soap/functions.cpp
View file @ 1d9df5f7
......@@ -86,21 +86,26 @@ void ModifiedSphericalBessel1stKind::evaluate(double r, bool differentiate) {
_in = ModifiedSphericalBessel1stKind::eval(_degree, r);
if (differentiate) {
_din.resize(_degree+1, 0.);
if (r < RADZERO) {
_din.push_back(0.0);
_din.push_back(1./3.);
for (int n = 1; n <= _degree; ++n) {
_din.push_back(0.);
}
_din[1] = 1./3.;
//_din.push_back(0.0);
//_din.push_back(1./3.);
//for (int n = 2; n <= _degree; ++n) {
// _din.push_back(0.);
//}
}
else {
_din.push_back( _in[1] );
_din[0] = _in[1];
for (int n = 1; n <= _degree; ++n) {
_din.push_back( _in[n-1] - (n+1.)/r*_in[n] );
_din[n] = _in[n-1] - (n+1.)/r*_in[n];
}
//_din.push_back( _in[1] );
//for (int n = 1; n <= _degree; ++n) {
// _din.push_back( _in[n-1] - (n+1.)/r*_in[n] );
//}
}
}
return;
}
......
src/soap/radialbasis.cpp
View file @ 1d9df5f7
......@@ -110,11 +110,12 @@ void RadialBasisGaussian::configure(Options &options) {
throw std::runtime_error("Not implemented.");
}
// SUMMARIZE
GLOG() << "Created " << _N << " radial Gaussians at r = { ";
GLOG() << "Created " << _N << " radial Gaussians at:" << std::endl;
for (basis_it_t bit = _basis.begin(); bit != _basis.end(); ++bit) {
GLOG() << (*bit)->_r0 << " (" << (*bit)->_sigma << ") ";
GLOG() << boost::format(" r = %1$+1.7e") % (*bit)->_r0
<< boost::format(" (sigma = %1$+1.7e") % (*bit)->_sigma << ") "
<< std::endl;
}
GLOG() << "}" << std::endl;
// COMPUTE OVERLAP MATRIX s_{ij} = \int g_i(r) g_j(r) r^2 dr
_Sij = ub::matrix<double>(_N, _N);
......@@ -178,6 +179,105 @@ void RadialBasisGaussian::configure(Options &options) {
}
}
// For each l: integral S r^2 dr i_l(2*ai*ri*r) exp(-beta_ik*(r-rho_ik)^2)
void compute_integrals_il_expik_r2_dr(
double ai,
double ri,
double beta_ik,
double rho_ik,
int L_plus_1,
int n_steps,
std::vector<double> *integrals,
std::vector<double> *integrals_derivative) {
bool gradients = (integrals_derivative != NULL) ? true : false;
// Sanity checks
assert(integrals->size() == L_plus_1);
if (gradients) assert(integrals->size() == L_plus_1);
// Sample coordinates along r-axis
double sigma_ik = sqrt(0.5/beta_ik);
double r_min = rho_ik - 4*sigma_ik;
double r_max = rho_ik + 4*sigma_ik;
if (r_min < 0.) {
r_max -= r_min;
r_min = 0.;
}
double delta_r_step = (r_max-r_min)/(n_steps-1);
int n_sample = 2*n_steps+1;
double delta_r_sample = 0.5*delta_r_step;
ModifiedSphericalBessel1stKind mosbest(L_plus_1-1);
if (gradients) {
// Compute samples for all l's ...
ub::matrix<double> integrand_l_at_r = ub::zero_matrix<double>(L_plus_1, n_sample);
ub::matrix<double> integrand_derivative_l_at_r = ub::zero_matrix<double>(L_plus_1, n_sample);
for (int s = 0; s < n_sample; ++s) {
double r_sample = r_min - delta_r_sample + s*delta_r_sample;
double exp_ik = exp(-beta_ik*(r_sample-rho_ik)*(r_sample-rho_ik));
mosbest.evaluate(2*ai*ri*r_sample, gradients);
for (int l = 0; l != L_plus_1; ++l) {
integrand_l_at_r(l,s) =
r_sample*r_sample*
mosbest._in[l]*
exp_ik;
integrand_derivative_l_at_r(l,s) =
2*ai*r_sample*r_sample*r_sample*
mosbest._din[l]*
exp_ik;
}
}
// ... integrate (à la Simpson)
std::vector<double> &ints = *integrals;
std::vector<double> &ints_deriv = *integrals_derivative;
for (int s = 0; s < n_steps; ++s) {
for (int l = 0; l != L_plus_1; ++l) {
ints[l] += delta_r_step/6.*(
integrand_l_at_r(l, 2*s)+
4*integrand_l_at_r(l, 2*s+1)+
integrand_l_at_r(l, 2*s+2)
);
ints_deriv[l] += delta_r_step/6.*(
integrand_derivative_l_at_r(l, 2*s)+
4*integrand_derivative_l_at_r(l, 2*s+1)+
integrand_derivative_l_at_r(l, 2*s+2)
);
}
}
}
else {
// Compute samples ...
ub::matrix<double> integrand_l_at_r = ub::zero_matrix<double>(L_plus_1, n_sample);
for (int s = 0; s < n_sample; ++s) {
double r_sample = r_min - delta_r_sample + s*delta_r_sample;
double exp_ik = exp(-beta_ik*(r_sample-rho_ik)*(r_sample-rho_ik));
mosbest.evaluate(2*ai*ri*r_sample, gradients);
for (int l = 0; l != L_plus_1; ++l) {
integrand_l_at_r(l,s) =
r_sample*r_sample*
mosbest._in[l]*
exp_ik;
}
}
// ... integrate (à la Simpson)
std::vector<double> &ints = *integrals;
for (int s = 0; s < n_steps; ++s) {
for (int l = 0; l != L_plus_1; ++l) {
ints[l] += delta_r_step/6.*(
integrand_l_at_r(l, 2*s)+
4*integrand_l_at_r(l, 2*s+1)+
integrand_l_at_r(l, 2*s+2)
);
}
}
}
return;
}
void RadialBasisGaussian::computeCoefficients(
vec d,
double r,
......@@ -196,6 +296,9 @@ void RadialBasisGaussian::computeCoefficients(
// Delta-type expansion =>
// Second (l) dimension of <save_here> and <particle_sigma> ignored here
if (particle_sigma < RadialBasis::RADZERO) {
if (gradients) {
throw soap::base::NotImplemented("<RadialBasisGaussian::computeCoefficients> Gradients when sigma=0.");
}
basis_it_t it;
int n = 0;
for (it = _basis.begin(), n = 0; it != _basis.end(); ++it, ++n) {
......@@ -207,7 +310,6 @@ void RadialBasisGaussian::computeCoefficients(
Gnl = ub::prod(_Tij, Gnl);
}
else {
// Particle properties
double ai = 1./(2*particle_sigma*particle_sigma);
double ri = r;
......@@ -217,78 +319,65 @@ void RadialBasisGaussian::computeCoefficients(
int k = 0;
basis_it_t it;
for (it = _basis.begin(), k = 0; it != _basis.end(); ++it, ++k) {
// Radial Gaussian properties
// Prefactor (r-independent)
double ak = (*it)->_alpha;
double rk = (*it)->_r0;
double norm_r2_g2_dr_rad_k = (*it)->_norm_r2_g2_dr;
// Combined properties
double beta_ik = ai+ak;
double rho_ik = ak*rk/beta_ik;
double sigma_ik = sqrt(0.5/beta_ik);
RadialGaussian gik_rad(rho_ik, sigma_ik);
double norm_r2_g_dr_rad_ik = gik_rad._norm_r2_g_dr;
double prefac =
4*M_PI *
norm_r2_g2_dr_rad_k*norm_g_dV_sph_i/norm_r2_g_dr_rad_ik *
norm_r2_g2_dr_rad_k*norm_g_dV_sph_i *
exp(-ai*ri*ri) *
exp(-ak*rk*rk*(1-ak/beta_ik));
// NUMERICAL INTEGRATION
// ZERO COEFFS
// Zero coeffs (to be safe ...)
for (int l = 0; l != Gnl.size2(); ++l) {
Gnl(k, l) = 0.0;
}
// SIMPSON'S RULE
double r_min = rho_ik - 4*sigma_ik;
double r_max = rho_ik + 4*sigma_ik;
if (r_min < 0.) {
r_max -= r_min;
r_min = 0.;
if (gradients) {
for (int l = 0; l != Gnl.size2(); ++l) {
(*dGnl_dx)(k,l) = 0.0;
(*dGnl_dy)(k,l) = 0.0;
(*dGnl_dz)(k,l) = 0.0;
}
}
int n_steps = this->_integration_steps;
double delta_r_step = (r_max-r_min)/(n_steps-1);
int n_sample = 2*n_steps+1;
double delta_r_sample = 0.5*delta_r_step;
// Compute samples for all l's
// For each l, store integrand at r_sample
ub::matrix<double> integrand_l_at_r = ub::zero_matrix<double>(Gnl.size2(), n_sample);
// TODO ub::matrix<double> grad_integrand_l_at_r = ...
// TODO Apply prefactor later (after integration)
// TODO gamma_kl = prefactor*integral_kl
// TODO grad_gamma_kl = -2*ai*ri*gamma_kl + prefactor*grad_integral_kl
for (int s = 0; s < n_sample; ++s) {
// f0 f1 f2 f3 .... f-3 f-2 f-1
// |-----||----||----||-------|
double r_sample = r_min - delta_r_sample + s*delta_r_sample;
// ... Generate Bessels
std::vector<double> sph_il = ModifiedSphericalBessel1stKind::eval(Gnl.size2(), 2*ai*ri*r_sample);
// ... Compute & store integrands
for (int l = 0; l != Gnl.size2(); ++l) {
integrand_l_at_r(l, s) =
prefac*
r_sample*r_sample*
sph_il[l]*
exp(-beta_ik*(r_sample-rho_ik)*(r_sample-rho_ik))*norm_r2_g_dr_rad_ik;
}
// Compute integrals S r^2 dr i_l(2*ai*ri*r) exp(-beta_ik*(r-rho_ik)^2)
// and (derivative) S 2*ai*r^3 dr i_l(2*ai*ri*r) exp(-beta_ik*(r-rho_ik)^2)
if (gradients) {
std::vector<double> integrals;
integrals.resize(Gnl.size2(), 0.);
std::vector<double> integrals_derivative;
integrals_derivative.resize(Gnl.size2(), 0.);
compute_integrals_il_expik_r2_dr(
ai, ri, beta_ik, rho_ik, Gnl.size2(), _integration_steps,
&integrals, &integrals_derivative);
for (int l = 0; l != Gnl.size2(); ++l) {
Gnl(k, l) = prefac*integrals[l];
double dgkl = -4*M_PI*2*ai*ri*Gnl(k,l) + prefac*integrals_derivative[l];
(*dGnl_dx)(k,l) = dgkl*d.getX();
(*dGnl_dy)(k,l) = dgkl*d.getY();
(*dGnl_dz)(k,l) = dgkl*d.getZ();
}
}
// Apply Simpson's rule
for (int s = 0; s < n_steps; ++s) {
for (int l = 0; l != Gnl.size2(); ++l) {
Gnl(k,l) +=
delta_r_step/6.*(
integrand_l_at_r(l, 2*s)+
4*integrand_l_at_r(l, 2*s+1)+
integrand_l_at_r(l, 2*s+2)
);
}
else {
std::vector<double> integrals;
integrals.resize(Gnl.size2(), 0.);
compute_integrals_il_expik_r2_dr(
ai, ri, beta_ik, rho_ik, Gnl.size2(), _integration_steps,
&integrals, NULL);
for (int l = 0; l != Gnl.size2(); ++l) {
Gnl(k, l) = prefac*integrals[l];
}
}
}
Gnl = ub::prod(_Tij, Gnl);
if (gradients) {
(*dGnl_dx) = ub::prod(_Tij, *dGnl_dx);
(*dGnl_dy) = ub::prod(_Tij, *dGnl_dy);
(*dGnl_dz) = ub::prod(_Tij, *dGnl_dz);
}
}
return;
}
......
src/soap/radialbasis.hpp
View file @ 1d9df5f7
......@@ -114,6 +114,15 @@ public:
}
};