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

Density spectrum - bare scaffold.

parent 8e0402fa
#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>
......@@ -16,6 +17,18 @@ void RadialBasis::configure(Options &options) {
_Rc = options.get<double>("radialbasis.Rc");
}
RadialCoefficients RadialBasis::computeCoefficients(double r) {
RadialCoefficients coeffs(this->_N);
throw soap::base::NotImplemented("RadialBasis::computeCoefficients");
return coeffs;
}
RadialCoefficients RadialBasis::computeCoefficientsAllZero() {
RadialCoefficients coeffs(this->_N);
for (int i = 0; i < _N; ++i) coeffs(i) = 0.0;
return coeffs;
}
void RadialBasisGaussian::configure(Options &options) {
RadialBasis::configure(options);
......@@ -23,20 +36,40 @@ void RadialBasisGaussian::configure(Options &options) {
// CREATE & STORE EQUISPACED RADIAL GAUSSIANS
this->clear();
double dr = _Rc/(_N-1);
for (int i = 0; i < _N; ++i) {
double r = i*dr;
double sigma = _sigma;
basis_fct_t *new_fct = new basis_fct_t(r, sigma);
_basis.push_back(new_fct);
std::string mode = "adaptive";
mode = "equispaced";
if (mode == "equispaced") {
double dr = _Rc/(_N-1);
for (int i = 0; i < _N; ++i) {
double r = i*dr;
double sigma = _sigma;
basis_fct_t *new_fct = new basis_fct_t(r, sigma);
_basis.push_back(new_fct);
}
}
else if (mode == "adaptive") {
double delta = 0.5;
int L = 6;
double r = 0.;
double sigma = 0.;
while (r < _Rc) {
sigma = sqrt(4./(2*L+1))*(r+delta);
basis_fct_t *new_fct = new basis_fct_t(r, sigma);
_basis.push_back(new_fct);
r = r + sigma;
}
}
else {
throw std::runtime_error("Not implemented.");
}
// SUMMARIZE
GLOG() << "Created " << _N << " radial Gaussians at r = { ";
for (basis_it_t bit = _basis.begin(); bit != _basis.end(); ++bit) {
GLOG() << (*bit)->_r0 << " ";
GLOG() << (*bit)->_r0 << " (" << (*bit)->_sigma << ") ";
}
GLOG() << "}" << std::endl;
// COMPUTE OVERLAP MATRIX
// COMPUTE OVERLAP MATRIX s_{ij} = \int g_i(r) g_j(r) r^2 dr
_Sij = ub::matrix<double>(_N, _N);
basis_it_t it;
basis_it_t jt;
......@@ -44,18 +77,6 @@ void RadialBasisGaussian::configure(Options &options) {
int j;
for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i) {
for (jt = _basis.begin(), j = 0; jt != _basis.end(); ++jt, ++j) {
/*
double a = (*it)->_alpha;
double a32 = pow(a, 1.5);
double b = (*it)->_r0;
double c = (*jt)->_r0;
double pre = 4.*M_PI/(16.*a32);
double s = pre*exp(-a*(b*b+c*c)) * (
sqrt(2.*M_PI)*(1.+a*(b+c)*(b+c))*exp(0.5*a*(b+c)*(b+c)) * \
(1. - boost::math::erf<double>(-1.*sqrt(0.5*a)*(b+c)))
+ 2.*sqrt(a)*(b+c)
);
*/
double a = (*it)->_alpha;
double b = (*jt)->_alpha;
double r0 = (*it)->_r0;
......@@ -69,7 +90,6 @@ void RadialBasisGaussian::configure(Options &options) {
1 - boost::math::erf<double>(-W0/sqrt(w))
)
);
s *= (*it)->_norm_dV*(*jt)->_norm_dV;
_Sij(i,j) = s;
}
......@@ -107,19 +127,31 @@ void RadialBasisGaussian::configure(Options &options) {
}
}
RadialCoefficients RadialBasisGaussian::computeCoefficients(double r) {
ub::vector<double> coeffs_raw(_N);
basis_it_t it;
int i;
for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i) {
coeffs_raw(i) = (*it)->at(r);
}
coeffs_raw = ub::prod(_Tij, coeffs_raw);
RadialCoefficients coeffs(this->_N);
for (int n = 0; n < _N; ++n) {
coeffs.set(n, coeffs_raw(n));
}
return coeffs;
}
RadialGaussian::RadialGaussian(double r0, double sigma)
: _r0(r0),
_sigma(sigma),
_alpha(1./(2.*sigma*sigma)) {
// COMPUTE NORMALIZATION
/*
double alpha2 = 2*_alpha;
_integral_4_pi_r2_g2_dr =
pow(M_PI/alpha2, 1.5) * (1.+2.*alpha2*_r0*_r0)*(1.-boost::math::erf<double>(-sqrt(alpha2)*_r0))
+ pow(M_PI/alpha2, 1.5) * 2*sqrt(alpha2/M_PI)*_r0*exp(-alpha2*_r0*_r0);
_norm_dV = 1./std::sqrt(_integral_4_pi_r2_g2_dr);
*/
double w = 2*_alpha;
double W0 = 2*_alpha*_r0;
_integral_r2_g2_dr =
......@@ -132,6 +164,12 @@ RadialGaussian::RadialGaussian(double r0, double sigma)
_norm_dV = 1./sqrt(_integral_r2_g2_dr);
}
double RadialGaussian::at(double r) {
double p = _alpha*(r-_r0)*(r-_r0);
if (p < 40) return _norm_dV * exp(-p);
else return 0.0;
}
/*
self.r0 = r0
self.sigma = sigma
......
......@@ -5,17 +5,40 @@
#include "options.hpp"
#include <string>
#include <boost/numeric/ublas/matrix.hpp>
#include "base/exceptions.hpp"
namespace soap {
namespace ub = boost::numeric::ublas;
class RadialCoefficients : public ub::vector<double>
{
public:
RadialCoefficients(int N) {
this->resize(N);
for (int i = 0; i < N; ++i) (*this)[i] = 0.0;
}
void set(int n, double c) {
if (this->checkSize(n)) (*this)[n] = c;
else throw soap::base::OutOfRange("RadialCoefficients::set");
}
double &get(int n) {
if (this->checkSize(n)) return (*this)[n];
else throw soap::base::OutOfRange("RadialCoefficients::get");
}
bool checkSize(int n) { return this->size() > n; }
};
class RadialBasis
{
public:
std::string &identify() { return _type; }
const int &N() { return _N; }
virtual ~RadialBasis() {;}
virtual void configure(Options &options);
virtual RadialCoefficients computeCoefficients(double r);
virtual RadialCoefficients computeCoefficientsAllZero();
protected:
bool _is_ortho;
......@@ -24,6 +47,7 @@ protected:
double _Rc;
};
struct RadialGaussian
{
RadialGaussian(double r0, double sigma);
......@@ -31,7 +55,6 @@ struct RadialGaussian
double _r0;
double _sigma;
double _alpha;
double _integral_4_pi_r2_g2_dr;
double _integral_r2_g2_dr;
double _norm_dV;
};
......@@ -58,6 +81,7 @@ public:
_basis.clear();
}
void configure(Options &options);
RadialCoefficients computeCoefficients(double r);
protected:
double _sigma;
basis_t _basis;
......
......@@ -6,19 +6,162 @@ namespace soap {
Spectrum::Spectrum(Structure &structure, Options &options)
: _structure(&structure), _options(&options), _log(NULL) {
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);
}
Spectrum::~Spectrum() {
delete _log;
_log = NULL;
delete _radbasis;
_radbasis = NULL;
delete _angbasis;
_angbasis = 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;
Structure::particle_it_t pit;
for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
this->computeAtomic(*pit);
// TODO Receive & store atomic spectrum
break;
}
}
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();
for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
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);
}
// 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.);
for (int n = 0; n < _radbasis->N(); ++n) {
for (int l = 0; l <= _angbasis->L(); ++l) {
for (int m = -l; m <= l; ++m) {
std::complex<double> Gn_Ylm = c_nlm.get(n, l, m);
density_dr += c_nlm.get(n, l, m)*c_nlm_dr.get(n, l, m);
}
}
}
std::cout << density_dr.real() << " " << std::flush;
if ( ((k+I) % 6) == 5 ) {
std::cout << std::endl;
}
}
}
}
for (int n = 0; n < _radbasis->N(); ++n) {
for (int l = 0; l <= _angbasis->L(); ++l) {
for (int m = -l; m <= l; ++m) {
std::complex<double> Gn_Ylm = c_nlm.get(n, l, m);
}
}
}
return;
}
void Spectrum::registerPython() {
......
......@@ -5,6 +5,7 @@
#include "types.hpp"
#include "globals.hpp"
#include "radialbasis.hpp"
#include "angularbasis.hpp"
#include "base/logger.hpp"
namespace soap {
......@@ -45,6 +46,7 @@ public:
void clean();
void compute();
void computeAtomic(Particle *center);
void computePower();
void computeLinear();
......@@ -55,6 +57,7 @@ private:
Options *_options;
Structure *_structure;
RadialBasis *_radbasis;
AngularBasis *_angbasis;
std::string _radbasis_type;
Logger *_log;
......
......@@ -30,6 +30,7 @@ public:
using namespace boost::python;
class_<Particle, Particle*>("Particle", init<int>())
.add_property("pos", &Particle::getPosNumeric, &Particle::setPosNumeric)
.add_property("id", &Particle::getId, &Particle::setId)
.add_property("name", make_function(&Particle::getName, copy_non_const()), &Particle::setName)
.add_property("type", make_function(&Particle::getType, copy_non_const()), &Particle::setType)
.add_property("type_id", make_function(&Particle::getTypeId, copy_non_const()), &Particle::setTypeId)
......@@ -47,6 +48,9 @@ public:
void setName(std::string name) { _name = name; }
std::string &getName() { return _name; }
void setId(int id) { _id = id; }
int getId() { return _id; }
void setTypeId(int id) { _type_id = id; }
int &getTypeId() { return _type_id; }
......@@ -103,6 +107,8 @@ private:
class Structure
{
public:
typedef std::vector<Particle*>::iterator particle_it_t;
Structure(std::string label) : _id(-1), _label(label), _box(NULL) {
matrix box;
box.ZeroMatrix();
......
......@@ -51,6 +51,7 @@ def setup_structure_ase(label, ase_config):
particle.pos = pos
particle.mass = mass
particle.name = name
particle.type = name
particle.type_id = typ
return structure
......
......@@ -16,11 +16,16 @@ options.set('radialbasis.type', 'gaussian')
options.set('radialbasis.N', 8)
options.set('radialbasis.Rc', 7.)
options.set('radialbasis.sigma', 0.5)
options.set('angularbasis.L', 6)
# LOAD STRUCTURE
structure = soap.tools.setup_structure_ase(config.config_file, config.atoms)
soap.tools.write_xyz('test.xyz', structure)
osio << osio.mg << structure.label << endl
for atom in structure:
osio << atom.name << atom.type << atom.pos << endl
# COMPUTE SPECTRUM
spectrum = soap.Spectrum(structure, options)
spectrum.compute()
......@@ -31,7 +36,6 @@ spectrum.compute()
osio.okquit()
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment