Commit 2d279396 authored by Carl Poelking's avatar Carl Poelking
Browse files

BasisExpansion: gradients I.

parent 85b8eded
#! /bin/bash
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 <vector>
#include <boost/format.hpp>
#include <gtest/gtest.h>
#include <soap/cutoff.hpp>
#include <soap/options.hpp>
class TestCutoffShiftedCosine : public ::testing::Test
{
public:
soap::Options _options;
soap::CutoffFunction *_cutoff;
virtual void SetUp() {
_options.set("radialcutoff.type", "shifted-cosine");
_options.set("radialcutoff.Rc", 4.);
_options.set("radialcutoff.Rc_width", 0.5);
_options.set("radialcutoff.center_weight", 1.);
soap::CutoffFunctionFactory::registerAll();
_cutoff = soap::CutoffFunctionOutlet().create(_options.get<std::string>("radialcutoff.type"));
_cutoff->configure(_options);
}
virtual void TearDown() {
delete _cutoff;
_cutoff = NULL;
}
};
TEST_F(TestCutoffShiftedCosine, Constructor) {
EXPECT_EQ(_cutoff->identify(), "shifted-cosine");
EXPECT_DOUBLE_EQ(_cutoff->getCenterWeight(), 1.);
EXPECT_DOUBLE_EQ(_cutoff->getCutoffWidth(), 0.5);
EXPECT_DOUBLE_EQ(_cutoff->getCutoff(), 4.);
}
TEST_F(TestCutoffShiftedCosine, Weight) {
double r0 = 3.5;
double r1 = 4.;
double dr = r1-r0;
EXPECT_DOUBLE_EQ(_cutoff->calculateWeight(r1+1.), -1e-10);
EXPECT_DOUBLE_EQ(_cutoff->calculateWeight(r1+0.00001), -1e-10);
EXPECT_DOUBLE_EQ(_cutoff->calculateWeight(r1), 0.);
EXPECT_DOUBLE_EQ(_cutoff->calculateWeight(r1-2.*dr/3.), 0.75);
EXPECT_DOUBLE_EQ(_cutoff->calculateWeight(r0), 1.);
/*
std::vector<double> r_list;
r_list.push_back(0.);
for (int ridx = 0; ridx < 100; ++ridx) {
double r = 3.45+ridx*0.01;
r_list.push_back(r);
}
for (int rd = 0; rd < r_list.size(); ++rd) {
double r = r_list[rd];
double w = _cutoff->calculateWeight(r);
std::cout << r << " " << w << std::endl;
}
*/
}
TEST_F(TestCutoffShiftedCosine, GradientWeight) {
soap::vec d1(std::sqrt(0.5),0.,std::sqrt(0.5));
soap::vec gw;
double abs_gw;
soap::vec norm_gw;
double diff;
double r0 = 3.5;
double r1 = 4.;
double dr = r1-r0;
gw = _cutoff->calculateGradientWeight(r0, d1);
abs_gw = soap::linalg::abs(gw);
EXPECT_DOUBLE_EQ(abs_gw, 0.);
gw = _cutoff->calculateGradientWeight(r1, d1);
abs_gw = soap::linalg::abs(gw);
EXPECT_NEAR(abs_gw, 0., 1e-10);
gw = _cutoff->calculateGradientWeight(r1-5./6.*dr, d1);
abs_gw = soap::linalg::abs(gw);
norm_gw = gw/abs_gw;
diff = soap::linalg::abs(norm_gw-d1);
EXPECT_NEAR(abs_gw, 0.5*0.5*M_PI/dr, 1e-10);
EXPECT_NEAR(diff, 2., 1e-10);
/*
std::vector< std::pair<double, soap::vec> > rd_list;
rd_list.push_back(std::pair<double, soap::vec>(0., d1));
for (int rd = 0; rd < 100; ++rd) {
double r = 3.45+rd*0.01;
rd_list.push_back(std::pair<double, soap::vec>(r, d1));
}
for (int rd = 0; rd < rd_list.size(); ++rd) {
double r = rd_list[rd].first;
soap::vec d = rd_list[rd].second;
soap::vec grad = _cutoff->calculateGradientWeight(r, d);
double abs_grad = soap::linalg::abs(grad);
std::cout << r << " " << abs_grad << std::endl;
}
*/
}
......@@ -10,24 +10,53 @@ void AngularBasis::configure(Options &options) {
_L = options.get<int>("angularbasis.L");
}
void AngularBasis::computeCoefficients(vec d, double r, angcoeff_t &save_here) {
void AngularBasis::computeCoefficients(
vec d,
double r,
double sigma,
angcoeff_t &Ylm,
angcoeff_t *dYlm_dx,
angcoeff_t *dYlm_dy,
angcoeff_t *dYlm_dz) {
if (r < AngularBasis::RADZERO) {
// SCALARS
std::complex<double> c00 = boost::math::spherical_harmonic<double,double>(0, 0, 0.0, 0.0);
save_here(0) = c00;
Ylm(0) = c00;
// GRADIENTS
if (dYlm_dx) {
assert(dYlm_dy != NULL && dYlm_dz != NULL);
; // <- all zero
}
}
else {
// Radius not zero => compute moments
double theta = acos(d.getZ());
double phi = atan2(d.getY(), d.getX());
// Shift [-pi, -0] to [pi, 2*pi]
if (phi < 0.) phi += 2*M_PI;
if (phi < 0.) phi += 2*M_PI; // <- Shift [-pi, -0] to [pi, 2*pi]
// SCALARS
for (int l = 0; l <= _L; ++l) {
for (int m = -l; m <= l; ++m) {
std::complex<double> clm = boost::math::spherical_harmonic<double,double>(l, m, theta, phi);
save_here(l*l+l+m) = clm;
Ylm(l*l+l+m) = clm;
}
}
// GRADIENTS
if (dYlm_dx) {
assert(dYlm_dy != NULL && dYlm_dz != NULL);
vec dr = r*d;
for (int l = 0; l <= _L; ++l) {
for (int m = -l; m <= l; ++m) {
int lm = l*l+l+m;
std::vector<std::complex<double> > dYlm = GradSphericalYlm::eval(l, m, dr);
(*dYlm_dx)(lm) = dYlm[0];
(*dYlm_dy)(lm) = dYlm[1];
(*dYlm_dz)(lm) = dYlm[2];
}
}
}
}
return;
}
......
......@@ -8,6 +8,7 @@
#include "soap/base/exceptions.hpp"
#include "soap/base/objectfactory.hpp"
#include "soap/options.hpp"
#include "soap/functions.hpp"
namespace soap {
......@@ -26,7 +27,14 @@ public:
std::string &identify() { return _type; }
const int &L() { return _L; }
virtual void configure(Options &options);
virtual void computeCoefficients(vec d, double r, angcoeff_t &save_here);
virtual void computeCoefficients(
vec d,
double r,
double sigma,
angcoeff_t &Ylm,
angcoeff_t *dYlm_dx,
angcoeff_t *dYlm_dy,
angcoeff_t *dYlm_dz);
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
......
......@@ -47,12 +47,17 @@ void Basis::registerPython() {
BasisExpansion::BasisExpansion(Basis *basis) :
_basis(basis),
_radbasis(basis->getRadBasis()),
_angbasis(basis->getAngBasis()) {
int L = _angbasis->L();
int N = _radbasis->N();
_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));
_angbasis(basis->getAngBasis()),
_has_coeff(false),
_has_coeff_grad(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));
}
}
BasisExpansion::~BasisExpansion() {
......@@ -63,17 +68,90 @@ BasisExpansion::~BasisExpansion() {
_angcoeff.clear();
}
void BasisExpansion::computeCoefficients(double r, vec d, double weight, double sigma) {
_radbasis->computeCoefficients(r, sigma, _radcoeff);
_angbasis->computeCoefficients(d, r, _angcoeff);
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) {
_coeff(n, l*l+l+m) = _radcoeff(n,l)*_angcoeff(l*l+l+m);
}
}
}
_coeff *= weight;
void BasisExpansion::computeCoefficients(double r, vec d) {
// Corresponds to: weight=weight_scale=1, sigma=0, gradients=false
_radbasis->computeCoefficients(d, r, 0., _radcoeff, NULL, NULL, NULL);
_angbasis->computeCoefficients(d, r, 0., _angcoeff, NULL, NULL, NULL);
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) {
_coeff(n, l*l+l+m) = _radcoeff(n,l)*_angcoeff(l*l+l+m);
}
}
}
return;
}
void BasisExpansion::computeCoefficients(double r, vec d, double weight, double weight_scale, double sigma, bool gradients) {
// SETUP STORAGE
int L = _angbasis->L();
int N = _radbasis->N();
if (!_has_coeff) {
// Should have already been done in constructor ::BasisExpansion(Basis *basis)
_has_coeff = 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;
_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);
_angcoeff_grad_x = AngularBasis::angcoeff_zero_t((L+1)*(L+1));
_angcoeff_grad_y = AngularBasis::angcoeff_zero_t((L+1)*(L+1));
_angcoeff_grad_z = AngularBasis::angcoeff_zero_t((L+1)*(L+1));
_coeff_grad_x = coeff_zero_t(N,(L+1)*(L+1));
_coeff_grad_y = coeff_zero_t(N,(L+1)*(L+1));
_coeff_grad_z = coeff_zero_t(N,(L+1)*(L+1));
}
// COMPUTE
if (_has_coeff && !_has_coeff_grad) {
_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) {
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) {
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) {
_coeff(n, l*l+l+m) = _radcoeff(n,l)*_angcoeff(l*l+l+m);
}
}
}
_coeff *= weight*weight_scale;
}
if (_has_coeff_grad) {
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) {
int lm = l*l+l+m;
_coeff_grad_x(n, lm) = weight*(
weight_scale*_radcoeff(n,l)*_angcoeff_grad_x(lm)
+ weight_scale*_radcoeff_grad_x(n,l)*_angcoeff(lm)
+ _weight_scale_grad.getX()*_radcoeff(n,l)*_angcoeff(lm)
);
_coeff_grad_y(n, lm) = weight*(
weight_scale*_radcoeff(n,l)*_angcoeff_grad_y(lm)
+ weight_scale*_radcoeff_grad_y(n,l)*_angcoeff(lm)
+ _weight_scale_grad.getY()*_radcoeff(n,l)*_angcoeff(lm)
);
_coeff_grad_z(n, lm) = weight*(
weight_scale*_radcoeff(n,l)*_angcoeff_grad_z(lm)
+ weight_scale*_radcoeff_grad_z(n,l)*_angcoeff(lm)
+ _weight_scale_grad.getZ()*_radcoeff(n,l)*_angcoeff(lm)
);
} // m
} // l
} // n
}
return;
}
......@@ -187,7 +265,8 @@ void BasisExpansion::writeDensityOnGrid(
// DENSITY BASED ON EXPANSION
if (fromExpansion) {
BasisExpansion density_exp_dr(_basis);
density_exp_dr.computeCoefficients(r, d, 1., 0.);
//density_exp_dr.computeCoefficients(r, d, 1., 1., 0., false);
density_exp_dr.computeCoefficients(r, d);
density_exp_dr.conjugate();
BasisExpansion::coeff_t &c_nlm_dr = density_exp_dr.getCoefficients();
......
......@@ -56,7 +56,7 @@ public:
typedef ub::zero_matrix< std::complex<double> > coeff_zero_t;
BasisExpansion(Basis *basis);
BasisExpansion() : _basis(NULL), _radbasis(NULL), _angbasis(NULL) {;}
BasisExpansion() : _basis(NULL), _radbasis(NULL), _angbasis(NULL), _has_coeff(false), _has_coeff_grad(false) {;}
~BasisExpansion();
Basis *getBasis() { return _basis; }
......@@ -64,7 +64,10 @@ public:
RadialBasis::radcoeff_t &getRadCoeffs() { return _radcoeff; }
AngularBasis::angcoeff_t &getAngCoeffs() { return _angcoeff; }
void computeCoefficients(double r, vec d, double weight, double sigma);
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; }
void add(BasisExpansion &other) { _coeff = _coeff + other._coeff; }
void conjugate();
void writeDensity(std::string filename, Options *options,
......@@ -81,7 +84,14 @@ public:
arch & _basis;
arch & _radbasis;
arch & _angbasis;
arch & _has_coeff;
arch & _coeff;
arch & _has_coeff_grad;
arch & _coeff_grad_x;
arch & _coeff_grad_y;
arch & _coeff_grad_z;
}
private:
......@@ -89,9 +99,22 @@ private:
RadialBasis *_radbasis;
AngularBasis *_angbasis;
bool _has_coeff;
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;
vec _weight_scale_grad;
RadialBasis::radcoeff_t _radcoeff_grad_x;
RadialBasis::radcoeff_t _radcoeff_grad_y;
RadialBasis::radcoeff_t _radcoeff_grad_z;
AngularBasis::angcoeff_t _angcoeff_grad_x;
AngularBasis::angcoeff_t _angcoeff_grad_y;
AngularBasis::angcoeff_t _angcoeff_grad_z;
coeff_t _coeff_grad_x;
coeff_t _coeff_grad_y;
coeff_t _coeff_grad_z;
};
}
......
......@@ -14,6 +14,41 @@ void CutoffFunction::configure(Options &options) {
<< ", central weight = " << _center_weight << std::endl;
}
bool CutoffFunction::isWithinCutoff(double r) {
return (r <= _Rc);
}
double CutoffFunction::calculateWeight(double r) {
double weight_at_r = 1.;
if (r > _Rc) {
weight_at_r = -1.e-10;
}
else if (r <= _Rc - _Rc_width) {
weight_at_r = 1.;
}
else {
weight_at_r = 0.5*(1+cos(M_PI*(r-_Rc+_Rc_width)/_Rc_width));
//weight_at_r = pow(cos(0.5*M_PI*(r-_Rc+_Rc_width)/_Rc_width),2);
}
return weight_at_r;
}
vec CutoffFunction::calculateGradientWeight(double r, vec d) {
vec grad_weight(0.,0.,0.);
if (r > _Rc) {
;
}
else if (r <= _Rc - _Rc_width) {
;
}
else {
grad_weight = -0.5*sin(M_PI*(r-_Rc+_Rc_width)/_Rc_width)*M_PI/_Rc_width * d;
//double phi = 0.5*M_PI*(r-_Rc+_Rc_width)/_Rc_width;
//grad_weight = - M_PI/_Rc_width*sin(phi)*cos(phi) * d;
}
return grad_weight;
}
void CutoffFunctionFactory::registerAll(void) {
CutoffFunctionOutlet().Register<CutoffFunction>("shifted-cosine");
}
......
......@@ -28,19 +28,9 @@ public:
_center_weight(-1.) {;}
virtual ~CutoffFunction() {;}
virtual void configure(Options &options);
virtual double calculateWeight(double r) {
double weight_at_r = 1.;
if (r > _Rc) {
weight_at_r = -1.;
}
else if (r <= _Rc - _Rc_width) {
weight_at_r = 1.;
}
else {
weight_at_r = 0.5*(1+cos(M_PI*(r-_Rc+_Rc_width)/_Rc_width));
}
return weight_at_r;
}
virtual bool isWithinCutoff(double r);
virtual double calculateWeight(double r);
virtual vec calculateGradientWeight(double r, vec d);
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
......@@ -51,11 +41,13 @@ public:
}
protected:
std::string _type;
double _Rc;
double _Rc_width;
double _center_weight;
public:
static constexpr double WEIGHT_ZERO = 1e-10;
};
......
......@@ -140,7 +140,9 @@ std::vector<std::complex<double> > GradSphericalYlm::eval(int l, int m, vec &r)
// TODO WHAT IF RADIUS IS ZERO?
double R = soap::linalg::abs(r);
if (R < RADZERO) {
throw soap::base::NotImplemented("<GradSphericalYlm::eval> R < RADZERO");
// This should have been checked before
// All gradients set to zero
throw soap::base::APIError("<GradSphericalYlm::eval> R < RADZERO");
}
vec d = r / R;
double x = r.getX();
......
......@@ -7,6 +7,7 @@ Options::Options() :
_center_excludes(boost::python::list()) {
// Set defaults
this->set("spectrum.gradients", false);
this->set("radialbasis.type", "gaussian");
this->set("radialbasis.mode", "equispaced");
this->set("radialbasis.N", 9);
......@@ -69,6 +70,7 @@ void Options::registerPython() {
void (Options::*set_int)(std::string, int) = &Options::set;
void (Options::*set_double)(std::string, double) = &Options::set;
void (Options::*set_string)(std::string, std::string) = &Options::set;
//void (Options::*set_bool)(std::string, bool) = &Options::set;
class_<Options, Options*>("Options")
.def("__str__", &Options::summarizeOptions)
......
......@@ -26,6 +26,7 @@ public:
void set(std::string key, std::string value) { _key_value_map[key] = value; }
void set(std::string key, int value) { this->set(key, boost::lexical_cast<std::string>(value)); }
void set(std::string key, double value) { this->set(key, boost::lexical_cast<std::string>(value)); }
//void set(std::string key, bool value) { this->set(key, boost::lexical_cast<std::string>(value)); }
void configureCenters(boost::python::list center_excludes) { _center_excludes = center_excludes; }
std::string summarizeOptions();
......
......@@ -28,7 +28,14 @@ void RadialBasis::configure(Options &options) {
_mode = options.get<std::string>("radialbasis.mode");
}
void RadialBasis::computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here) {
void RadialBasis::computeCoefficients(
vec d,
double r,
double particle_sigma,
radcoeff_t &Gnl,
radcoeff_t *dGnl_dx,
radcoeff_t *dGnl_dy,
radcoeff_t *dGnl_dz) {
throw soap::base::NotImplemented("RadialBasis::computeCoefficients");
return;
}
......@@ -171,7 +178,21 @@ void RadialBasisGaussian::configure(Options &options) {
}
}
void RadialBasisGaussian::computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here) {
void RadialBasisGaussian::computeCoefficients(
vec d,
double r,
double particle_sigma,
radcoeff_t &Gnl,
radcoeff_t *dGnl_dx,
radcoeff_t *dGnl_dy,
radcoeff_t *dGnl_dz) {
bool gradients = false;
if (dGnl_dx) {
assert(dGnl_dy != NULL && dGnl_dz != NULL);
gradients = true;
}
// Delta-type expansion =>
// Second (l) dimension of <save_here> and <particle_sigma> ignored here
if (particle_sigma < RadialBasis::RADZERO) {
......@@ -179,22 +200,13 @@ void RadialBasisGaussian::computeCoefficients(double r, double particle_sigma, r
int n = 0;
for (it = _basis.begin(), n = 0; it != _basis.end(); ++it, ++n) {
double gn_at_r = (*it)->at(r);
for (int l = 0; l != save_here.size2(); ++l) {
save_here(n, l) = gn_at_r;
for (int l = 0; l != Gnl.size2(); ++l) {
Gnl(n, l) = gn_at_r;
}
}
save_here = ub::prod(_Tij, save_here);
Gnl = ub::prod(_Tij, Gnl);
}
else {
// std::cout << "sigma > 0." << std::endl;
// int degree_sph_in = save_here.size2();
// std::vector<double> sph_il = ModifiedSphericalBessel1stKind::eval(degree_sph_in, 1e-4);
// for (int l = 0; l != save_here.size2(); ++l) {
// std::cout << "l=" << l << " il=" << sph_il[l] << std::endl;
// }
//
// throw soap::base::NotImplemented("...");
// Particle properties
double ai = 1./(2*particle_sigma*particle_sigma);
......@@ -223,59 +235,12 @@ void RadialBasisGaussian::computeCoefficients(double r, double particle_sigma, r
exp(-ai*ri*ri) *
exp(-ak*rk*rk*(1-ak/beta_ik));
// // DELTA APPROXIMATION
// std::cout << "r " << r;
// std::cout << " beta_ik " << beta_ik;
// std::cout << " rho_ik " << rho_ik;
// std::cout << std::endl;
// double bessel_arg_ik = 2*ai*ri*rho_ik;
// int degree_sph = save_here.size2(); // <- L+1
// std::vector<double> sph_il = ModifiedSphericalBessel1stKind::eval(degree_sph, bessel_arg_ik);
// for (int l = 0; l != degree_sph; ++l) {
// save_here(k, l) = prefac*sph_il[l];
// }
// std::cout << std::endl;
// std::cout << "DELTA" << std::endl;
// std::cout << "k = " << k;
// for (int l = 0; l != save_here.size2(); ++l) {
// std::cout << " " << save_here(k, l);