Commit e84cb3aa authored by Carl Poelking's avatar Carl Poelking
Browse files

Improved infrastructure: basis expansions.

parent f8a2d603
......@@ -10,26 +10,47 @@ void AngularBasis::configure(Options &options) {
}
AngularCoefficients AngularBasis::computeCoefficients(vec d, double r) {
AngularCoefficients coeffs(this->_L);
if (r < 1e-10) {
// Radius is virtually zero => set isotropic
std::complex<double> c00 = boost::math::spherical_harmonic<double,double>(0, 0, 0.0, 0.0);
coeffs.set(0, 0, c00);
}
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;
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);
coeffs.set(l, m, clm);
}
}
}
return coeffs;
AngularCoefficients coeffs(this->_L);
if (r < 1e-10) {
// Radius is virtually zero => set isotropic
std::complex<double> c00 = boost::math::spherical_harmonic<double,double>(0, 0, 0.0, 0.0);
coeffs.set(0, 0, c00);
}
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;
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);
coeffs.set(l, m, clm);
}
}
}
return coeffs;
}
void AngularBasis::computeCoefficients(vec d, double r, angcoeff_t &save_here) {
if (r < AngularBasis::RADZERO) {
std::complex<double> c00 = boost::math::spherical_harmonic<double,double>(0, 0, 0.0, 0.0);
save_here(0) = c00;
}
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;
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;
}
}
}
return;
}
AngularCoefficients AngularBasis::computeCoefficientsAllZero() {
......@@ -42,6 +63,10 @@ AngularCoefficients AngularBasis::computeCoefficientsAllZero() {
return coeffs;
}
void AngularBasisFactory::registerAll(void) {
AngularBasisOutlet().Register<AngularBasis>("spherical-harmonic");
}
}
......@@ -4,8 +4,9 @@
#include <string>
#include <math.h>
#include <vector>
#include "base/exceptions.hpp"
#include "base/exceptions.hpp"
#include "base/objectfactory.hpp"
#include "options.hpp"
namespace soap {
......@@ -15,6 +16,7 @@ namespace ub = boost::numeric::ublas;
class AngularCoefficients : public ub::vector< std::complex<double> >
{
public:
AngularCoefficients(int L) : _L(L) {
this->resize((L+1)*(L+1));
for (int i=0; i!=size(); ++i) {
......@@ -43,11 +45,15 @@ protected:
class AngularBasis
{
public:
typedef ub::vector< std::complex<double> > angcoeff_t;
typedef ub::zero_vector< std::complex<double> > angcoeff_zero_t;
std::string &identify() { return _type; }
const int &L() { return _L; }
AngularBasis() : _type("spherical-harmonic"), _L(0) {;}
virtual ~AngularBasis() {;}
virtual void configure(Options &options);
virtual void computeCoefficients(vec d, double r, angcoeff_t &save_here);
virtual AngularCoefficients computeCoefficients(vec d, double r);
virtual AngularCoefficients computeCoefficientsAllZero();
......@@ -55,10 +61,36 @@ protected:
std::string _type;
int _L;
static const double RADZERO = 1e-10;
};
class AngularBasisFactory
: public soap::base::ObjectFactory<std::string, AngularBasis>
{
private:
AngularBasisFactory() {}
public:
static void registerAll(void);
AngularBasis *create(const std::string &key);
friend AngularBasisFactory &AngularBasisOutlet();
};
inline AngularBasisFactory &AngularBasisOutlet() {
static AngularBasisFactory _instance;
return _instance;
}
inline AngularBasis *AngularBasisFactory::create(const std::string &key) {
assoc_map::const_iterator it(getObjects().find(key));
if (it != getObjects().end()) {
AngularBasis *basis = (it->second)();
return basis;
}
else {
throw std::runtime_error("Factory key " + key + " not found.");
}
}
}
......
......@@ -18,6 +18,18 @@ public:
explicit NotImplemented(std::string mssg) : std::runtime_error("NotImplemented["+mssg+"]") { ; }
};
class IOError : public std::runtime_error
{
public:
explicit IOError(std::string mssg) : std::runtime_error("IOError["+mssg+"]") { ; }
};
class APIError : public std::runtime_error
{
public:
explicit APIError(std::string mssg) : std::runtime_error("UsageError["+mssg+"]") { ; }
};
}}
#endif
......@@ -14,5 +14,6 @@ BOOST_PYTHON_MODULE(_soapxx)
soap::Spectrum::registerPython();
soap::RadialBasisFactory::registerAll();
soap::AngularBasisFactory::registerAll();
}
......@@ -7,6 +7,12 @@ namespace soap {
extern Logger GLOG;
namespace constants {
const double ANGSTROM_TO_BOHR = 1./0.52917721067;
}
}
#endif
#include "radialbasis.hpp"
#include "globals.hpp"
#include "linalg/operations.hpp"
#include "base/exceptions.hpp"
#include <math.h>
#include <boost/math/special_functions/erf.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include "base/exceptions.hpp"
#include "linalg/operations.hpp"
#include "globals.hpp"
#include "radialbasis.hpp"
namespace soap {
namespace ub = boost::numeric::ublas;
// ======================
// RadialBasis BASE CLASS
// ======================
void RadialBasis::configure(Options &options) {
_N = options.get<int>("radialbasis.N");
......@@ -23,6 +27,11 @@ RadialCoefficients RadialBasis::computeCoefficients(double r) {
return coeffs;
}
void RadialBasis::computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here) {
throw soap::base::NotImplemented("RadialBasis::computeCoefficients");
return;
}
RadialCoefficients RadialBasis::computeCoefficientsAllZero() {
RadialCoefficients coeffs(this->_N);
for (int i = 0; i < _N; ++i) coeffs(i) = 0.0;
......@@ -30,6 +39,10 @@ RadialCoefficients RadialBasis::computeCoefficientsAllZero() {
}
// ======================
// RadialBasisGaussian
// ======================
void RadialBasisGaussian::configure(Options &options) {
RadialBasis::configure(options);
_sigma = options.get<double>("radialbasis.sigma");
......@@ -90,7 +103,7 @@ void RadialBasisGaussian::configure(Options &options) {
1 - boost::math::erf<double>(-W0/sqrt(w))
)
);
s *= (*it)->_norm_dV*(*jt)->_norm_dV;
s *= (*it)->_norm_r2_g2_dr*(*jt)->_norm_r2_g2_dr;
_Sij(i,j) = s;
}
}
......@@ -146,46 +159,106 @@ RadialCoefficients RadialBasisGaussian::computeCoefficients(double r) {
return coeffs;
}
struct SphericalGaussian
{
SphericalGaussian(vec r0, double sigma) :
_r0(r0), _sigma(sigma), _alpha(1./(2*sigma*sigma)) {
_norm_g_dV = pow(_alpha/M_PI, 1.5);
}
vec _r0;
double _sigma;
double _alpha;
double _norm_g_dV;
};
void RadialBasisGaussian::computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here) {
// Delta-type expansion =>
// Second (l) dimension of <save_here> and <particle_sigma> ignored here
if (particle_sigma < RadialBasis::RADZERO) {
basis_it_t it;
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;
}
}
save_here = ub::prod(_Tij, save_here);
}
else {
throw soap::base::NotImplemented("...");
// Particle properties
double ai = 1./(2*particle_sigma*particle_sigma);
double ri = r;
basis_it_t it;
int k = 0;
for (it = _basis.begin(), k = 0; it != _basis.end(); ++it, ++k) {
// Radial Gaussian properties
double ak = (*it)->_alpha;
double rk = (*it)->_r0;
// Combined properties
double beta_ik = ai+ak;
double rho_ik = ak*rk/beta_ik;
double bessel_arg_ik = 2*ai*ri*rho_ik;
}
}
return;
}
// ======================
// RadialGaussian
// ======================
RadialGaussian::RadialGaussian(double r0, double sigma)
: _r0(r0),
_sigma(sigma),
_alpha(1./(2.*sigma*sigma)) {
// COMPUTE NORMALIZATION
// COMPUTE NORMALIZATION S g^2 r^2 dr
// This normalization is to be used for radial basis functions
double w = 2*_alpha;
double W0 = 2*_alpha*_r0;
_integral_r2_g2_dr =
1./(4.*pow(w, 2.5))*exp(-2*_alpha*_r0*_r0)*(
1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
2*sqrt(w)*W0 +
sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
1 - boost::math::erf<double>(-W0/sqrt(w))
)
);
_norm_dV = 1./sqrt(_integral_r2_g2_dr);
_norm_r2_g2_dr = 1./sqrt(_integral_r2_g2_dr);
// COMPUTE NORMALIZATION S g r^2 dr
// This normalization is to be used for "standard" radial Gaussians
w = _alpha;
W0 = _alpha*_r0;
_integral_r2_g_dr =
1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
2*sqrt(w)*W0 +
sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
1 - boost::math::erf<double>(-W0/sqrt(w))
)
);
_norm_r2_g_dr = 1./_integral_r2_g_dr;
}
double RadialGaussian::at(double r) {
double p = _alpha*(r-_r0)*(r-_r0);
if (p < 40) return _norm_dV * exp(-p);
if (p < 40) return _norm_r2_g2_dr * exp(-p);
else return 0.0;
}
/*
self.r0 = r0
self.sigma = sigma
self.alpha = 1./(2*sigma**2)
# Normalization
alpha2 = 2*self.alpha
self.integral_4_pi_r2_g2_dr = (np.pi/alpha2)**1.5*(1.+2.*alpha2*self.r0**2)*(1.-scipy.special.erf(-alpha2**0.5*self.r0)) + (np.pi/alpha2)**1.5*2*(alpha2/np.pi)**0.5*self.r0*np.exp(-alpha2*self.r0**2)
self.norm_dV = self.integral_4_pi_r2_g2_dr**(-0.5)
return
def EvaluateExp(self, r):
return np.exp(-(r-self.r0)**2/(2.*self.sigma**2))
def __call__(self, r):
return self.EvaluateExp(r)*self.norm_dV
*/
// ======================
// RadialBasisFactory
// ======================
void RadialBasisFactory::registerAll(void) {
RadialBasisOutlet().Register<RadialBasisGaussian>("gaussian");
......
#ifndef _SOAP_RADIALBASIS_HPP
#define _SOAP_RADIALBASIS_HPP
#include "base/objectfactory.hpp"
#include "options.hpp"
#include <string>
#include <boost/numeric/ublas/matrix.hpp>
#include "base/exceptions.hpp"
#include "base/objectfactory.hpp"
#include "options.hpp"
namespace soap {
......@@ -33,11 +34,15 @@ public:
class RadialBasis
{
public:
typedef ub::matrix<double> radcoeff_t;
typedef ub::zero_matrix<double> radcoeff_zero_t;
std::string &identify() { return _type; }
const int &N() { return _N; }
virtual ~RadialBasis() {;}
virtual void configure(Options &options);
virtual RadialCoefficients computeCoefficients(double r);
virtual void computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here);
virtual RadialCoefficients computeCoefficientsAllZero();
protected:
......@@ -45,6 +50,8 @@ protected:
std::string _type;
int _N;
double _Rc;
static const double RADZERO = 1e-10;
};
......@@ -55,8 +62,12 @@ struct RadialGaussian
double _r0;
double _sigma;
double _alpha;
// Integral S g^2 r^2 dr
double _integral_r2_g2_dr;
double _norm_dV;
double _norm_r2_g2_dr;
// Integral S g r^2 dr
double _integral_r2_g_dr;
double _norm_r2_g_dr;
};
class RadialBasisGaussian : public RadialBasis
......@@ -82,6 +93,7 @@ public:
}
void configure(Options &options);
RadialCoefficients computeCoefficients(double r);
void computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here);
protected:
double _sigma;
basis_t _basis;
......
#include <fstream>
#include <boost/format.hpp>
#include "spectrum.hpp"
namespace soap {
Spectrum::Spectrum(Structure &structure, Options &options)
: _structure(&structure), _options(&options), _log(NULL) {
: _log(NULL), _options(&options), _structure(&structure) {
GLOG() << "Configuring spectrum ..." << std::endl;
// CONFIGURE RADIAL BASIS
_radbasis = RadialBasisOutlet().create(_options->get<std::string>("radialbasis.type"));
_radbasis->configure(options);
// CONFIGURE ANGULAR BASIS
_angbasis = new AngularBasis();
_angbasis->configure(options);
// CREATE & CONFIGURE BASIS
_basis = new Basis(&options);
}
Spectrum::~Spectrum() {
delete _log;
_log = NULL;
delete _radbasis;
_radbasis = NULL;
delete _angbasis;
_angbasis = NULL;
delete _basis;
_basis = NULL;
}
void Spectrum::compute() {
GLOG() << "Compute spectrum ..." << std::endl;
GLOG() << _options->summarizeOptions() << std::endl;
GLOG() << "Using radial basis of type '" << _radbasis->identify() << "'" << std::endl;
GLOG() << "Using angular basis of type '" << _angbasis->identify() << "'" << std::endl;
GLOG() << "Using radial basis of type '" << _basis->getRadBasis()->identify() << "'" << std::endl;
GLOG() << "Using angular basis of type '" << _basis->getAngBasis()->identify() << "'" << std::endl;
Structure::particle_it_t pit;
for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
......@@ -39,128 +33,71 @@ void Spectrum::compute() {
}
}
// TODO change to ::computeAtomic(Particle *center, vector<Particle*> &nbhood)
void Spectrum::computeAtomic(Particle *center) {
GLOG() << "Compute atomic spectrum for particle " << center->getId()
<< " (type " << center->getType() << ")" << std::endl;
Structure::particle_it_t pit;
struct BasisCoefficients
{
BasisCoefficients(RadialCoefficients &c_n, AngularCoefficients &c_lm)
: _c_n(c_n), _c_lm(c_lm) {
_c_nlm.resize(c_n.size(), c_lm.size());
_c_nlm = ub::outer_prod(_c_n, _c_lm);
}
std::complex<double> &get(int n, int l, int m) {
if (_c_lm.checkSize(l, m) && _c_n.checkSize(n)) {
return _c_nlm(n, l*l+l+m);
}
else {
throw soap::base::OutOfRange("BasisCoefficients::get");
}
}
void add(BasisCoefficients &other) {
_c_nlm = _c_nlm + other._c_nlm;
}
ub::matrix< std::complex<double> > _c_nlm;
RadialCoefficients _c_n;
AngularCoefficients _c_lm;
};
RadialCoefficients c_n_zero = _radbasis->computeCoefficientsAllZero();
AngularCoefficients c_lm_zero = _angbasis->computeCoefficientsAllZero();
BasisCoefficients c_nlm(c_n_zero, c_lm_zero);
std::string type_this = center->getType();
std::string type_this = center->getType();
// RadialCoefficients c_n_zero = _radbasis->computeCoefficientsAllZero();
// AngularCoefficients c_lm_zero = _angbasis->computeCoefficientsAllZero();
// BasisCoefficients c_nlm(c_n_zero, c_lm_zero);
// c_nlm.linkBasis(_radbasis, _angbasis);
BasisExpansion *nbhood_expansion = new BasisExpansion(this->_basis);
Structure::particle_it_t pit;
for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
// TODO Cut-off check
vec dr = _structure->connect(center->getPos(), (*pit)->getPos());
double r = soap::linalg::abs(dr);
// TODO Cut-off check
vec d = dr/r;
std::string type_other = (*pit)->getType();
// COMPUTE RADIAL COEFFICIENTS
RadialCoefficients c_n_pair = _radbasis->computeCoefficients(r);
/*
GLOG() << "Radial coefficients r = " << r << std::endl;
for (int i = 0; i < c_n.size(); ++i) {
GLOG() << c_n[i] << " | ";
}
GLOG() << std::endl;
*/
// COMPUTE ANGULAR COEFFICIENTS
AngularCoefficients c_lm_pair = _angbasis->computeCoefficients(d, r);
/*
GLOG() << "Angular coefficients d = " << d << std::endl;
for (int lm = 0; lm < c_lm.size(); ++lm) {
GLOG() << c_lm[lm] << " | ";
}
GLOG() << std::endl;
*/
BasisCoefficients c_nlm_pair(c_n_pair, c_lm_pair);
c_nlm.add(c_nlm_pair);
// TODO Account for weight
BasisExpansion nb_expansion(this->_basis);
nb_expansion.computeCoefficients(r, d, (*pit)->getWeight(), (*pit)->getSigma());
nbhood_expansion->add(nb_expansion);
// // COMPUTE RADIAL COEFFICIENTS
// RadialCoefficients c_n_pair = _radbasis->computeCoefficients(r);
// /*
// GLOG() << "Radial coefficients r = " << r << std::endl;
// for (int i = 0; i < c_n.size(); ++i) {
// GLOG() << c_n[i] << " | ";
// }
// GLOG() << std::endl;
// */
//
// // COMPUTE ANGULAR COEFFICIENTS
// AngularCoefficients c_lm_pair = _angbasis->computeCoefficients(d, r);
// /*
// GLOG() << "Angular coefficients d = " << d << std::endl;
// for (int lm = 0; lm < c_lm.size(); ++lm) {
// GLOG() << c_lm[lm] << " | ";
// }
// GLOG() << std::endl;
// */
//
// BasisCoefficients c_nlm_pair(c_n_pair, c_lm_pair);
// c_nlm.add(c_nlm_pair);
}
// SAMPLE EXPANDED DENSITY OVER GRID
int I = 20;
double dx = 0.1;
std::cout << "Density on grid" << std::endl;
std::cout << "Loop order X -> Y -> Z" << std::endl;
std::cout << "1 -2 -2 -2" << std::endl;
std::cout << "41 0.1 0.0 0.0" << std::endl;
std::cout << "41 0.0 0.1 0.0" << std::endl;
std::cout << "41 0.0 0.0 0.1" << std::endl;
std::cout << "6 0.0 0.0 0.0 0.0" << std::endl;
for (int i = -I; i <= I; ++i) {
//std::cout << i << std::endl;
for (int j = -I; j <= I; ++j) {
for (int k = -I; k <= I; ++k) {
vec dr(i*dx, j*dx, k*dx);
double r = soap::linalg::abs(dr);
vec d = dr/r;
RadialCoefficients c_n_dr = _radbasis->computeCoefficients(r);
AngularCoefficients c_lm_dr = _angbasis->computeCoefficients(d, r);
c_lm_dr.conjugate();
BasisCoefficients c_nlm_dr(c_n_dr, c_lm_dr);
//double density_dr = ub::inner_prod(c_nlm, c_nlm_dr);
std::complex<double> density_dr(0.,0.);