Commit 006720f9 authored by Carl Poelking's avatar Carl Poelking
Browse files

Nested derivatives and tests.

parent 7320fb8f
#! /bin/bash
soap_tests/test.exe --gtest_list_tests
#soap_tests/test.exe --gtest_list_tests
soap_tests/test.exe
#soap_tests/test.exe --gtest_filter=TestFunctions.GradientYlm
#soap_tests/test.exe --gtest_filter=TestCutoffShiftedCosine.* --gtest_output=xml
......
#include <iostream>
#include <fstream>
#include <vector>
#include <tuple>
#include <boost/format.hpp>
#include <boost/algorithm/string.hpp>
#include <gtest/gtest.h>
#include <soap/basis.hpp>
#include <soap/options.hpp>
#include "gtest_defines.hpp"
class TestBasisExpansion : public ::testing::Test
{
public:
soap::Options _options;
soap::Basis *_basis;
std::string _constructor_stdout;
std::string _constructor_stdout_ref;
virtual void SetUp() {
_options.set("radialbasis.type", "gaussian");
_options.set("radialbasis.mode", "adaptive");
_options.set("radialbasis.N", 9);
_options.set("radialbasis.sigma", 0.5);
_options.set("radialbasis.integration_steps", 15);
_options.set("radialcutoff.type", "shifted-cosine");
_options.set("radialcutoff.Rc", 4.);
_options.set("radialcutoff.Rc_width", 0.5);
_options.set("radialcutoff.center_weight", 1.);
_options.set("angularbasis.type", "spherical-harmonic");
_options.set("angularbasis.L", 6);
soap::RadialBasisFactory::registerAll();
soap::AngularBasisFactory::registerAll();
soap::CutoffFunctionFactory::registerAll();
::testing::internal::CaptureStdout();
_basis = new soap::Basis(&_options);
::testing::internal::GetCapturedStdout();
}
virtual void TearDown() {
delete _basis;
_basis = NULL;
}
};
TEST_F(TestBasisExpansion, Gradients) {
typedef std::tuple<int, int, int, double, double, double> nlmxyz_t;
std::vector<nlmxyz_t> nlmxyz_list = {
nlmxyz_t{0,0,0,1.,0.5,-0.5},
nlmxyz_t{0,1,0,1.,0.5,-0.5},
nlmxyz_t{0,3,-1,-1.,-0.5,1.},
nlmxyz_t{3,4,1,0.,0.3,-0.7},
nlmxyz_t{7,2,-1,-0.5,3.,4.}
};
std::stringstream output;
std::string output_ref = "n=0 l=0 m=0 x=+1.0000 y=+0.5000 z=-0.5000 Q_re=+9.4683121e-02 dQ_re=-1.8961753e-01 -9.4808765e-02 +9.4808765e-02 Q_im=+0.0000000e+00 dQ_im=+0.0000000e+00 +0.0000000e+00 +0.0000000e+00 n=0 l=1 m=0 x=+1.0000 y=+0.5000 z=-0.5000 Q_re=-4.7738237e-02 dQ_re=+1.0923622e-01 +5.4618111e-02 +4.0858364e-02 Q_im=-0.0000000e+00 dQ_im=+0.0000000e+00 +0.0000000e+00 +0.0000000e+00 n=0 l=3 m=-1 x=-1.0000 y=-0.5000 z=+1.0000 Q_re=-1.2554214e-02 dQ_re=-2.8055842e-02 -2.0305028e-02 -5.0416312e-03 Q_im=+6.2771069e-03 dQ_im=+2.0305028e-02 -2.4017000e-03 +2.5208156e-03 n=3 l=4 m=1 x=+0.0000 y=+0.3000 z=-0.7000 Q_re=-2.1144554e-20 dQ_re=-1.1510939e-03 +8.4096906e-20 -4.4449962e-19 Q_im=-3.4532818e-04 dQ_im=+2.2521571e-20 +1.6936152e-03 -4.1418340e-03 n=7 l=2 m=-1 x=-0.5000 y=+3.0000 z=+4.0000 Q_re=+7.0360862e-15 dQ_re=-1.1432964e-14 -1.5835248e-14 -1.9354642e-14 Q_im=+4.2216517e-14 dQ_im=+1.5835248e-14 -8.0939315e-14 -1.1612785e-13 ";
for (auto it = nlmxyz_list.begin(); it != nlmxyz_list.end(); ++it) {
// Extract parameters
int n = std::get<0>(*it);
int l = std::get<1>(*it);
int m = std::get<2>(*it);
double x = std::get<3>(*it);
double y = std::get<4>(*it);
double z = std::get<5>(*it);
int lm = l*l+l+m;
// Position, radius, direction
soap::vec R(x, y, z);
double r = soap::linalg::abs(R);
soap::vec d = R/r;
// Weights
double weight0 = 1.;
double weight_scale = _basis->getCutoff()->calculateWeight(r);
double sigma = 0.5;
bool gradients = true;
// Compute gradients
::testing::internal::CaptureStdout();
soap::BasisExpansion basex(_basis);
basex.computeCoefficients(r, d, weight0, weight_scale, sigma, gradients);
::testing::internal::GetCapturedStdout();
// Extract results
soap::BasisExpansion::coeff_t &coeff = basex.getCoefficients();
soap::BasisExpansion::coeff_t &coeff_grad_x = basex.getCoefficientsGradX();
soap::BasisExpansion::coeff_t &coeff_grad_y = basex.getCoefficientsGradY();
soap::BasisExpansion::coeff_t &coeff_grad_z = basex.getCoefficientsGradZ();
output
<< boost::format("n=%1$d l=%2$d m=%3$d ") % n % l % m
<< boost::format("x=%1$+1.4f y=%2$+1.4f z=%3$+1.4f ") % x % y % z
<< boost::format("Q_re=%1$+1.7e dQ_re=%2$+1.7e %3$+1.7e %4$+1.7e ")
% coeff(n,lm).real() % coeff_grad_x(n,lm).real() % coeff_grad_y(n,lm).real() % coeff_grad_z(n,lm).real()
<< boost::format("Q_im=%1$+1.7e dQ_im=%2$+1.7e %3$+1.7e %4$+1.7e ")
% coeff(n,lm).imag() % coeff_grad_x(n,lm).imag() % coeff_grad_y(n,lm).imag() % coeff_grad_z(n,lm).imag()
<< std::flush;
}
// VERIFY
EXPECT_EQ(output.str(), output_ref);
// GENERATE REFERENCE
// std::cout << "\"" << output.str() << "\"" << std::endl;
/*
// MANUAL TESTING
int N = 200;
double dx = 0.05;
std::string logfile = "tmp";
std::ofstream ofs;
ofs.open(logfile.c_str(), std::ofstream::out);
if (!ofs.is_open()) {
throw std::runtime_error("Bad file handle: " + logfile);
}
for (int i = 0; i < N; ++i) {
double z = -i*dx+5.;
double x = 0.2;
double y = -3.;
soap::vec R(x, y, z);
double r = soap::linalg::abs(R);
soap::vec d = R/r;
double weight0 = 1.;
double weight_scale = _basis->getCutoff()->calculateWeight(r);
double sigma = 0.5;
bool gradients = true;
// Compute
::testing::internal::CaptureStdout();
soap::BasisExpansion basex(_basis);
basex.computeCoefficients(r, d, weight0, weight_scale, sigma, gradients);
::testing::internal::GetCapturedStdout();
// Extract
soap::BasisExpansion::coeff_t &coeff = basex.getCoefficients();
soap::BasisExpansion::coeff_t &coeff_grad_x = basex.getCoefficientsGradX();
soap::BasisExpansion::coeff_t &coeff_grad_y = basex.getCoefficientsGradY();
soap::BasisExpansion::coeff_t &coeff_grad_z = basex.getCoefficientsGradZ();
int n = 3;
int l = 5;
int m = -4;
int lm = l*l+l+m;
ofs << boost::format("%1$+1.4f %2$+1.7e %3$+1.7e %4$+1.7e %5$+1.7e %6$+1.7e %7$+1.7e %8$+1.7e %9$+1.7e")
% R.getZ() % coeff(n,lm).real() % coeff(n,lm).imag()
% coeff_grad_z(n,lm).real() % coeff_grad_z(n,lm).imag()
% coeff_grad_y(n,lm).real() % coeff_grad_y(n,lm).imag()
% coeff_grad_x(n,lm).real() % coeff_grad_x(n,lm).imag()
<< std::endl;
}
ofs.close();
*/
}
#include <iostream>
#include <fstream>
#include <vector>
#include <boost/format.hpp>
#include <gtest/gtest.h>
......@@ -115,6 +116,35 @@ TEST_F(TestCutoffShiftedCosine, GradientWeight) {
std::cout << r << " " << abs_grad << std::endl;
}
*/
/*
// MANUAL TESTING
int N = 100;
double dx = 0.05;
std::string logfile = "tmp_cutoff";
std::ofstream ofs;
ofs.open(logfile.c_str(), std::ofstream::out);
if (!ofs.is_open()) {
throw std::runtime_error("Bad file handle: " + logfile);
}
for (int i = 0; i < N; ++i) {
soap::vec d(1.,0.,0.);
double r = i*dx;
double w = _cutoff->calculateWeight(r);
soap::vec w_grad = _cutoff->calculateGradientWeight(r, d);
ofs << boost::format("%1$+1.4f %2$+1.7e %3$+1.7e %4$+1.7e %5$+1.7e")
% r % w
% w_grad.getX()
% w_grad.getY()
% w_grad.getZ()
<< std::endl;
}
*/
}
......
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)
......@@ -5,9 +5,7 @@
#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)
#include "gtest_defines.hpp"
class TestRadialBasisGaussian : public ::testing::Test
{
......
......@@ -17,16 +17,21 @@ AtomicSpectrum::AtomicSpectrum(Particle *center, Basis *basis) {
}
AtomicSpectrum::~AtomicSpectrum() {
// Clean Qnlm's
for (map_qnlm_t::iterator it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) delete it->second;
// CLEAN QNLM'S
// ... Summed density expansions
for (auto it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) delete it->second;
_map_qnlm.clear();
// ... Generic density
if (_qnlm_generic) {
delete _qnlm_generic;
_qnlm_generic = NULL;
}
// Clean Xnkl's
for (map_xnkl_t::iterator it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) delete it->second;
// CLEAN XNKL'S
// ... Scalar spectra
for (auto it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) delete it->second;
_map_xnkl.clear();
// ... Generic power spectra
if (_xnkl_generic_coherent) {
delete _xnkl_generic_coherent;
_xnkl_generic_coherent = NULL;
......@@ -35,6 +40,22 @@ AtomicSpectrum::~AtomicSpectrum() {
delete _xnkl_generic_incoherent;
_xnkl_generic_incoherent = NULL;
}
// CLEAN PID-RESOLVED GRADIENTS
// ... Neighbour density expansions with gradients
for (auto it = _map_pid_qnlm.begin(); it != _map_pid_qnlm.end(); ++it) delete it->second.second;
_map_pid_qnlm.clear();
// ... Xnkl gradients
for (auto it = _map_pid_xnkl.begin(); it != _map_pid_xnkl.end(); ++it) {
for (auto jt = (*it).second.begin(); jt != (*it).second.end(); ++jt) {
delete (*jt).second;
}
(*it).second.clear();
}
_map_pid_xnkl.clear();
// ... Gradients (generic-coherent)
for (auto it = _map_pid_xnkl_gc.begin(); it != _map_pid_xnkl_gc.end(); ++it) {
delete it->second;
}
}
void AtomicSpectrum::null() {
......@@ -111,6 +132,25 @@ void AtomicSpectrum::addQnlm(std::string type, qnlm_t &nb_expansion) {
return;
}
void AtomicSpectrum::addQnlmNeighbour(Particle *nb, qnlm_t *nb_expansion) {
std::string type = nb->getType();
int id = nb->getId();
this->addQnlm(type, *nb_expansion);
if (nb == this->getCenter()) {
delete nb_expansion; // <- gradients should be zero, do not store
std::cout << "DO NOT STORE" << std::endl;
}
else {
auto it = _map_pid_qnlm.find(id);
if (it != _map_pid_qnlm.end()) {
throw soap::base::SanityCheckFailed("<AtomicSpectrum::addQnlm> Already have entry for pid.");
}
_map_pid_qnlm[id] = std::pair<std::string,qnlm_t*>(type, nb_expansion);
}
return;
}
AtomicSpectrum::qnlm_t *AtomicSpectrum::getQnlm(std::string type) {
if (type == "") {
return _qnlm_generic;
......@@ -124,7 +164,70 @@ AtomicSpectrum::qnlm_t *AtomicSpectrum::getQnlm(std::string type) {
}
}
void AtomicSpectrum::computePowerGradients() {
GLOG() << "CID " << _center->getId() << ": " << std::endl;
for (auto it1 = _map_pid_qnlm.begin(); it1 != _map_pid_qnlm.end(); ++it1) {
// Derivatives are taken with respect to the coordinates of this neighbour particle
int pi_id = it1->first;
std::string pi_type = it1->second.first;
qnlm_t *pi_qnlm = it1->second.second;
// Check and prepare storage in map
auto it = _map_pid_xnkl.find(pi_id);
if (it != _map_pid_xnkl.end()) {
soap::base::SanityCheckFailed("<AtomicSpectrum::computePowerGradients> Already have entry for pid.");
}
_map_pid_xnkl[pi_id] = map_xnkl_t();
map_xnkl_t &xnkl_map = _map_pid_xnkl[pi_id];
// Compute for all type pairs
GLOG() << " NID " << pi_id << " " << std::flush;
for (auto it2 = _map_qnlm.begin(); it2 != _map_qnlm.end(); ++it2) {
std::string type_other = it2->first;
qnlm_t *sum_qnlm_type_other = it2->second;
if (pi_type == type_other) {
// Type 1 == type 2
GLOG() << pi_type << "=" << type_other << " " << std::flush;
PowerExpansion *powex = new PowerExpansion(_basis);
powex->computeCoefficientsGradients(pi_qnlm, sum_qnlm_type_other, true);
// Store ...
type_pair_t types(pi_type, type_other);
auto it = xnkl_map.find(types);
assert(it == xnkl_map.end());
xnkl_map[types] = powex;
}
else {
// Type 1 != type 2
GLOG() << pi_type << ":" << type_other << " " << std::flush;
PowerExpansion *powex_12 = new PowerExpansion(_basis);
powex_12->computeCoefficientsGradients(pi_qnlm, sum_qnlm_type_other, false);
GLOG() << type_other << ":" << pi_type << " " << std::flush;
PowerExpansion *powex_21 = new PowerExpansion(_basis);
powex_21->computeCoefficientsGradients(sum_qnlm_type_other, pi_qnlm, false);
// Store ...
type_pair_t types_12(pi_type, type_other);
type_pair_t types_21(type_other, pi_type);
auto it12 = xnkl_map.find(types_12);
auto it21 = xnkl_map.find(types_21);
assert(it12 == xnkl_map.end() && it21 == xnkl_map.end());
xnkl_map[types_12] = powex_12;
xnkl_map[types_21] = powex_21;
}
}
if (_qnlm_generic) {
auto it = _map_pid_xnkl_gc.find(pi_id);
if (it != _map_pid_xnkl_gc.end()) assert(false && "Already have Xnkl_gc gradients for pid.");
GLOG() << "*:* " << std::flush;
PowerExpansion *powex = new PowerExpansion(_basis);
powex->computeCoefficientsGradients(pi_qnlm, _qnlm_generic, true);
_map_pid_xnkl_gc[pi_id] = powex;
}
GLOG() << std::endl;
}
return;
}
void AtomicSpectrum::computePower() {
// TODO Calling this function more than once with the same object
// TODO causes memory leaks => check for existing _map_xnkl entries
// Specific (i.e., type-dependent)
map_qnlm_t::iterator it1;
map_qnlm_t::iterator it2;
......
......@@ -22,11 +22,20 @@ namespace soap {
class AtomicSpectrum : public std::map<std::string, BasisExpansion*>
{
public:
// EXPANSION TYPES
typedef BasisExpansion qnlm_t;
typedef PowerExpansion xnkl_t;
typedef std::pair<std::string, std::string> type_pair_t;
// CONTAINERS FOR STORING SCALAR FIELDS
typedef std::map<std::string, qnlm_t*> map_qnlm_t; // <- key: center type, e.g. 'C'
typedef std::map<std::pair<std::string, std::string>, xnkl_t*> map_xnkl_t; // <- key: type string pair, e.g. ('C','H')
typedef std::pair<std::string, std::string> type_pair_t;
// CONTAINERS FOR STORING GRADIENTS
typedef std::map<int, std::pair<std::string,qnlm_t*> > map_pid_qnlm_t; // <- id=>(type;qnlm)
typedef std::map<int, map_xnkl_t> map_pid_xnkl_t; // <- id=>type=>xnkl
typedef std::map<int, xnkl_t*> map_pid_xnkl_gc_t; // <- id=>xnkl_generic_coherent
AtomicSpectrum(Particle *center, Basis *basis);
AtomicSpectrum() { this->null(); }
......@@ -43,10 +52,12 @@ public:
Basis *getBasis() { return _basis; }
// QNLM METHODS
void addQnlm(std::string type, qnlm_t &nb_expansion);
void addQnlmNeighbour(Particle *nb, qnlm_t *nb_expansion);
qnlm_t *getQnlm(std::string type);
qnlm_t *getQnlmGeneric() { return _qnlm_generic; }
// XNKL METHODS
void computePower();
void computePowerGradients();
xnkl_t *getPower(std::string type1, std::string type2);
xnkl_t *getXnkl(type_pair_t &types);
map_xnkl_t &getXnklMap() { return _map_xnkl; }
......@@ -70,6 +81,10 @@ public:
arch & _map_xnkl;
arch & _xnkl_generic_coherent;
arch & _xnkl_generic_incoherent;
// PID-resolved
arch & _map_pid_qnlm;
arch & _map_pid_xnkl;
arch & _map_pid_xnkl_gc;
return;
}
protected:
......@@ -78,13 +93,20 @@ protected:
vec _center_pos;
std::string _center_type;
Basis *_basis;
// DENSITY EXPANSION
map_qnlm_t _map_qnlm;
qnlm_t *_qnlm_generic;
// POWER DENSITY EXPANSION
map_xnkl_t _map_xnkl;
xnkl_t *_xnkl_generic_coherent;
xnkl_t *_xnkl_generic_incoherent;
// PID-RESOLVED (GRADIENTS)
map_pid_qnlm_t _map_pid_qnlm; // <- for gradients wrt position of neighbour with global id
map_pid_xnkl_t _map_pid_xnkl;
map_pid_xnkl_gc_t _map_pid_xnkl_gc;
};
......
......@@ -30,6 +30,13 @@ public:
explicit APIError(std::string mssg) : std::runtime_error("APIError["+mssg+"]") { ; }
};
class SanityCheckFailed : public std::runtime_error
{
public:
explicit SanityCheckFailed(std::string mssg) : std::runtime_error("SanityCheckFailed["+mssg+"]") { ; }
};
}}
#endif
......@@ -48,16 +48,15 @@ BasisExpansion::BasisExpansion(Basis *basis) :
_basis(basis),
_radbasis(basis->getRadBasis()),
_angbasis(basis->getAngBasis()),
_has_coeff(false),
_has_coeff_grad(false) {
_has_scalars(false),
_has_gradients(false) {
int L = _angbasis->L();
int N = _radbasis->N();
_has_coeff = true;
if (_has_coeff) {
_radcoeff = RadialBasis::radcoeff_zero_t(N,L+1);
_angcoeff = AngularBasis::angcoeff_zero_t((L+1)*(L+1));
_coeff = coeff_zero_t(N,(L+1)*(L+1));
}
// ZERO SCALAR-FIELD CONTAINERS
_has_scalars = true;
_radcoeff = RadialBasis::radcoeff_zero_t(N,L+1);
_angcoeff = AngularBasis::angcoeff_zero_t((L+1)*(L+1));
_coeff = coeff_zero_t(N,(L+1)*(L+1));
}
BasisExpansion::~BasisExpansion() {
......@@ -86,15 +85,15 @@ void BasisExpansion::computeCoefficients(double r, vec d, double weight, double
// SETUP STORAGE
int L = _angbasis->L();
int N = _radbasis->N();
if (!_has_coeff) {
if (!_has_scalars) {
// Should have already been done in constructor ::BasisExpansion(Basis *basis)
_has_coeff = true;
_has_scalars = true;
_radcoeff = RadialBasis::radcoeff_zero_t(N,L+1);
_angcoeff = AngularBasis::angcoeff_zero_t((L+1)*(L+1));
_coeff = coeff_zero_t(N,(L+1)*(L+1));
}
if (gradients) {
_has_coeff_grad = true;
_has_gradients = true;
_radcoeff_grad_x = RadialBasis::radcoeff_zero_t(N,L+1);
_radcoeff_grad_y = RadialBasis::radcoeff_zero_t(N,L+1);
_radcoeff_grad_z = RadialBasis::radcoeff_zero_t(N,L+1);
......@@ -106,18 +105,18 @@ void BasisExpansion::computeCoefficients(double r, vec d, double weight, double
_coeff_grad_z = coeff_zero_t(N,(L+1)*(L+1));
}
// COMPUTE
if (_has_coeff && !_has_coeff_grad) {
if (_has_scalars && !_has_gradients) {
_radbasis->computeCoefficients(d, r, sigma, _radcoeff, NULL, NULL, NULL);
_angbasis->computeCoefficients(d, r, sigma, _angcoeff, NULL, NULL, NULL);
}
else if (_has_coeff && _has_coeff_grad) {
else if (_has_scalars && _has_gradients) {
std::cout << "GRAD" << std::endl;
_radbasis->computeCoefficients(d, r, sigma, _radcoeff, &_radcoeff_grad_x, &_radcoeff_grad_y, &_radcoeff_grad_z);
_angbasis->computeCoefficients(d, r, sigma, _angcoeff, &_angcoeff_grad_x, &_angcoeff_grad_y, &_angcoeff_grad_z);
_weight_scale_grad = _basis->getCutoff()->calculateGradientWeight(r, d);
}
// MERGE
if (_has_coeff) {
if (_has_scalars) {
for (int n = 0; n != _radbasis->N(); ++n) {
for (int l = 0; l != _angbasis->L()+1; ++l) {
for (int m = -l; m != l+1; ++m) {
......@@ -127,7 +126,7 @@ void BasisExpansion::computeCoefficients(double r, vec d, double weight, double
}
_coeff *= weight*weight_scale;
}
if (_has_coeff_grad) {
if (_has_gradients) {
for (int n = 0; n != _radbasis->N(); ++n) {
for (int l = 0; l != _angbasis->L()+1; ++l) {
for (int m = -l; m != l+1; ++m) {
......@@ -152,6 +151,16 @@ void BasisExpansion::computeCoefficients(double r, vec d, double weight, double
} // n
}
// CLEAR INTERMEDIATE STORAGE (NOT REQUIRED LATER)
_radcoeff.clear();
_angcoeff.clear();
_radcoeff_grad_x.clear();
_radcoeff_grad_y.clear();
_radcoeff_grad_z.clear();
_angcoeff_grad_x.clear();
_angcoeff_grad_y.clear();
_angcoeff_grad_z.clear();
return;
}
......
......@@ -56,18 +56,21 @@ public:
typedef ub::zero_matrix< std::complex<double> > coeff_zero_t;
BasisExpansion(Basis *basis);
BasisExpansion() : _basis(NULL), _radbasis(NULL), _angbasis(NULL), _has_coeff(false), _has_coeff_grad(false) {;}
BasisExpansion() : _basis(NULL), _radbasis(NULL), _angbasis(NULL), _has_scalars(false), _has_gradients(false) {;}
~BasisExpansion();
Basis *getBasis() { return _basis; }
coeff_t &getCoefficients() { return _coeff; }
coeff_t &getCoefficientsGradX() { return _coeff_grad_x; }
coeff_t &getCoefficientsGradY() { return _coeff_grad_y; }
coeff_t &getCoefficientsGradZ() { return _coeff_grad_z; }
RadialBasis::radcoeff_t &getRadCoeffs() { return _radcoeff; }
AngularBasis::angcoeff_t &getAngCoeffs() { return _angcoeff; }
void computeCoefficients(double r, vec d);
void computeCoefficients(double r, vec d, double weight, double weight_scale, double sigma, bool gradients);
bool hasCoefficients() { return _has_coeff; }
bool hasCoefficientsGrad() { return _has_coeff_grad; }
bool hasScalars() { return _has_scalars; }
bool hasGradients() { return _has_gradients; }
void add(BasisExpansion &other) { _coeff = _coeff + other._coeff; }
void conjugate();
void writeDensity(std::string filename, Options *options,
......@@ -85,10 +88,10 @@ public:
arch & _radbasis;
arch & _angbasis;
arch & _has_coeff;
arch & _has_scalars;
arch & _coeff;
arch & _has_coeff_grad;
arch & _has_gradients;
arch & _coeff_grad_x;
arch & _coeff_grad_y;
arch & _coeff_grad_z;
......@@ -99,12 +102,12 @@ private:
RadialBasis *_radbasis;
AngularBasis *_angbasis;
bool _has_coeff;
bool _has_scalars;
RadialBasis::radcoeff_t _radcoeff; // access via (n,l)
AngularBasis::angcoeff_t _angcoeff; // access via (l*l+l+m)
coeff_t _coeff; // access via (n, l*l+l+m)
bool _has_coeff_grad;
bool _has_gradients;
vec _weight_scale_grad;