diff --git a/cxx_tests/src/gtest_radialbasis.cpp b/cxx_tests/src/gtest_radialbasis.cpp index 804ad7198a95dcad93cf27b68cd8bb19208864e4..11cbda5a420aa4e8ea5eeb7fcb2daed36a8a05af 100644 --- a/cxx_tests/src/gtest_radialbasis.cpp +++ b/cxx_tests/src/gtest_radialbasis.cpp @@ -115,7 +115,19 @@ TEST_F(TestRadialBasisGaussian, Constructor) { } } -TEST_F(TestRadialBasisGaussian, computeCoefficientsGradients) { +TEST_F(TestRadialBasisGaussian, computeCoefficientsGradients) { + + // PARAMETERS: SETUP + double sigma = 0.5; + soap::vec d(0.,1.,0.); + double dr = 0.5; + double N_r_samples = 10; + std::vector< std::pair<int,int> > kl_list { + std::pair<int,int>{0,1}, + std::pair<int,int>{3,0}, + std::pair<int,int>{7,2} }; + + // OUTPUT: SETUP int N = _radbasis->N(); int L = _options.get<int>("angularbasis.L"); @@ -123,14 +135,28 @@ TEST_F(TestRadialBasisGaussian, computeCoefficientsGradients) { 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); - */ - + + // GENERATE + ::testing::internal::CaptureStdout(); + for (int i = 0; i < N_r_samples; ++i) { + double ri = i*dr; + _radbasis->computeCoefficients(d, ri, sigma, Gnl, &dGnl_dx, &dGnl_dy, &dGnl_dz); + + for (auto kl = kl_list.begin(); kl != kl_list.end(); ++kl) { + int k = (*kl).first; + int l = (*kl).second; + std::cout << boost::format("ri=%1$+1.7f k=%2$d l=%3$d %4$+1.7e %5$+1.7e %6$+1.7e %7$+1.7e ") + % ri % k % l % Gnl(k,l) % (dGnl_dx)(k,l) % (dGnl_dy)(k,l) % (dGnl_dz)(k,l) << std::flush; + } + } + + // VERIFY + std::string capture = ::testing::internal::GetCapturedStdout(); + std::string capture_ref = "ri=+0.0000000 k=0 l=1 +0.0000000e+00 +0.0000000e+00 +1.1290706e+00 +0.0000000e+00 ri=+0.0000000 k=3 l=0 -1.6766282e-02 +0.0000000e+00 +0.0000000e+00 +0.0000000e+00 ri=+0.0000000 k=7 l=2 +0.0000000e+00 +0.0000000e+00 +0.0000000e+00 +0.0000000e+00 ri=+0.5000000 k=0 l=1 +4.1923545e-01 +0.0000000e+00 +3.3923961e-01 +0.0000000e+00 ri=+0.5000000 k=3 l=0 -2.0307498e-02 +0.0000000e+00 -2.1938577e-02 +0.0000000e+00 ri=+0.5000000 k=7 l=2 -7.4463158e-04 +0.0000000e+00 -1.6976156e-03 +0.0000000e+00 ri=+1.0000000 k=0 l=1 +3.4758879e-01 +0.0000000e+00 -4.5906064e-01 +0.0000000e+00 ri=+1.0000000 k=3 l=0 -5.1773726e-03 +0.0000000e+00 +1.5179894e-01 +0.0000000e+00 ri=+1.0000000 k=7 l=2 +6.5109983e-04 +0.0000000e+00 +6.7094414e-03 +0.0000000e+00 ri=+1.5000000 k=0 l=1 +1.2525628e-01 +0.0000000e+00 -3.3883494e-01 +0.0000000e+00 ri=+1.5000000 k=3 l=0 +1.4783849e-01 +0.0000000e+00 +3.8395021e-01 +0.0000000e+00 ri=+1.5000000 k=7 l=2 -3.2864325e-04 +0.0000000e+00 -1.5007935e-02 +0.0000000e+00 ri=+2.0000000 k=0 l=1 +2.4105014e-02 +0.0000000e+00 -9.3139826e-02 +0.0000000e+00 ri=+2.0000000 k=3 l=0 +2.6779871e-01 +0.0000000e+00 +2.5715774e-02 +0.0000000e+00 ri=+2.0000000 k=7 l=2 -4.9915564e-03 +0.0000000e+00 +1.3169613e-02 +0.0000000e+00 ri=+2.5000000 k=0 l=1 +2.6681970e-03 +0.0000000e+00 -1.3156527e-02 +0.0000000e+00 ri=+2.5000000 k=3 l=0 +1.8943025e-01 +0.0000000e+00 -2.6538685e-01 +0.0000000e+00 ri=+2.5000000 k=7 l=2 +1.2764517e-02 +0.0000000e+00 +3.4633567e-02 +0.0000000e+00 ri=+3.0000000 k=0 l=1 +1.7508377e-04 +0.0000000e+00 -1.0436941e-03 +0.0000000e+00 ri=+3.0000000 k=3 l=0 +6.9423484e-02 +0.0000000e+00 -1.7936480e-01 +0.0000000e+00 ri=+3.0000000 k=7 l=2 +6.0499578e-03 +0.0000000e+00 -6.2962003e-02 +0.0000000e+00 ri=+3.5000000 k=0 l=1 +6.8898763e-06 +0.0000000e+00 -4.8086339e-05 +0.0000000e+00 ri=+3.5000000 k=3 l=0 +1.4506579e-02 +0.0000000e+00 -5.3395498e-02 +0.0000000e+00 ri=+3.5000000 k=7 l=2 -2.4447509e-02 +0.0000000e+00 -2.6584282e-02 +0.0000000e+00 ri=+4.0000000 k=0 l=1 +1.6269418e-07 +0.0000000e+00 -1.3031468e-06 +0.0000000e+00 ri=+4.0000000 k=3 l=0 +1.7204100e-03 +0.0000000e+00 -8.3897684e-03 +0.0000000e+00 ri=+4.0000000 k=7 l=2 -1.4173474e-02 +0.0000000e+00 +4.6853015e-02 +0.0000000e+00 ri=+4.5000000 k=0 l=1 +2.2638085e-09 +0.0000000e+00 -2.0617355e-08 +0.0000000e+00 ri=+4.5000000 k=3 l=0 +1.0674873e-04 +0.0000000e+00 -6.7033707e-04 +0.0000000e+00 ri=+4.5000000 k=7 l=2 +1.5068192e-03 +0.0000000e+00 +9.9520461e-03 +0.0000000e+00 "; + EXPECT_EQ(capture, capture_ref); + + // GENERATE REFERENCE + //std::cout << capture << std::endl; } TEST_F(TestRadialBasisGaussian, NumericalIntegrationBessel) { diff --git a/src/soap/radialbasis.cpp b/src/soap/radialbasis.cpp index 2f394510a317549ee3b423b94191646a0eac81e2..56e76c1cf276636fb8bc0da2ef3f1171a6342391 100644 --- a/src/soap/radialbasis.cpp +++ b/src/soap/radialbasis.cpp @@ -372,6 +372,14 @@ void RadialBasisGaussian::computeCoefficients( } } } + +// for (int k = 0; k < (*dGnl_dx).size1(); ++k) { +// for (int l = 0; l < (*dGnl_dx).size2(); ++l) { +// std::cout << boost::format("%1$+1.7f %2$d %3$d %4$+1.7e %5$+1.7e %6$+1.7e %7$+1.7e") +// % r % k % l % Gnl(k,l) % (*dGnl_dx)(k,l) % (*dGnl_dy)(k,l) % (*dGnl_dz)(k,l) << std::endl; +// } +// } + Gnl = ub::prod(_Tij, Gnl); if (gradients) { (*dGnl_dx) = ub::prod(_Tij, *dGnl_dx); diff --git a/src/soap/spectrum.cpp b/src/soap/spectrum.cpp index 61ad74e38c1424c315f8b5be22c2030e0a1f8609..79e8943ebbd2cdbc8b750361113b7826959cd31b 100644 --- a/src/soap/spectrum.cpp +++ b/src/soap/spectrum.cpp @@ -96,15 +96,9 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar vec dr = _structure->connect(center->getPos(), (*pit)->getPos()); double r = soap::linalg::abs(dr); if (! this->_basis->getCutoff()->isWithinCutoff(r)) continue; - - //double weight_scale = this->_basis->getCutoff()->calculateWeight(r); - //if (weight_scale <= CutoffFunction::WEIGHT_ZERO) continue; vec d = dr/r; // APPLY WEIGHT IF CENTER - //if (*pit == center) { - // weight_scale *= this->_basis->getCutoff()->getCenterWeight(); - //} double weight0 = (*pit)->getWeight(); double weight_scale = _basis->getCutoff()->calculateWeight(r); if (*pit == center) {