Skip to content
GitLab
Commit 59cc0eed authored by Carl Poelking's avatar Carl Poelking
Browse files

Storage - atomic spectra.

parent 579137a5
Hide whitespace changes
Inline Side-by-side
src/basis.hpp
View file @ 59cc0eed
......@@ -12,6 +12,7 @@
#include "structure.hpp"
#include "angularbasis.hpp"
#include "radialbasis.hpp"
#include "cutoff.hpp"
namespace soap {
......@@ -28,19 +29,26 @@ public:
// CONFIGURE ANGULAR BASIS
_angbasis = AngularBasisOutlet().create(_options->get<std::string>("angularbasis.type"));
_angbasis->configure(*options);
// CONFIGURE CUTOFF FUNCTION
_cutoff = CutoffFunctionOutlet().create(_options->get<std::string>("radialcutoff.type"));
_cutoff->configure(*options);
}
~Basis() {
delete _radbasis;
_radbasis = NULL;
delete _angbasis;
_angbasis = NULL;
delete _cutoff;
_cutoff = NULL;
}
RadialBasis *getRadBasis() { return _radbasis; }
AngularBasis *getAngBasis() { return _angbasis; }
CutoffFunction *getCutoff() { return _cutoff; }
private:
Options *_options;
RadialBasis *_radbasis;
AngularBasis *_angbasis;
CutoffFunction *_cutoff;
};
......@@ -65,6 +73,7 @@ public:
_radcoeff.clear();
_angcoeff.clear();
}
Basis *getBasis() { return _basis; }
coeff_t &getCoefficients() { return _coeff; }
void computeCoefficients(double r, vec d, double weight, double sigma) {
_radbasis->computeCoefficients(r, sigma, _radcoeff);
......
src/bindings.cpp
View file @ 59cc0eed
......@@ -15,5 +15,6 @@ BOOST_PYTHON_MODULE(_soapxx)
soap::RadialBasisFactory::registerAll();
soap::AngularBasisFactory::registerAll();
soap::CutoffFunctionFactory::registerAll();
}
src/options.hpp
View file @ 59cc0eed
......@@ -15,7 +15,8 @@ public:
_b1(vec(1,0,0)), _b2(vec(0,1,0)), _b3(vec(0,0,1)),
_center_w0(1.), _center_excludes(boost::python::list()) {
// SET DEFAULTS
this->set("integration_steps", 15);
this->set("radialbasis.integration_steps", 15);
this->set("radialbasis.mode", "equispaced");
}
~Options() {}
......@@ -63,17 +64,22 @@ public:
std::string summarizeOptions() {
std::string info = "";
info += "Options:\n";
info += "o Centers: ";
info += (boost::format("W0=%1% Excl#=%2%\n") % _center_w0 % boost::python::len(_center_excludes)).str();
info += "o Real-space basis: ";
info += (boost::format("N=%1% L=%2% Rc=%3%\n") % _N_real % _L_real % _Rc_real).str();
info += "o Reciprocal-space basis: ";
info += (boost::format("N=%1% L=%2% Rc=%3%\n") % _N_recip % _L_recip % _Rc_recip).str();
info += "o Reciprocal-space lattice: ";
info += (boost::format("b1x=%1% b1y=%2% b1z=%3% ") % _b1.x() % _b1.y() % _b1.z()).str();
info += (boost::format("b2x=%1% b2y=%2% b2z=%3% ") % _b2.x() % _b2.y() % _b2.z()).str();
info += (boost::format("b3x=%1% b3y=%2% b3z=%3% ") % _b3.x() % _b3.y() % _b3.z()).str();
// info += "o Centers: ";
// info += (boost::format("W0=%1% Excl#=%2%\n") % _center_w0 % boost::python::len(_center_excludes)).str();
// info += "o Real-space basis: ";
// info += (boost::format("N=%1% L=%2% Rc=%3%\n") % _N_real % _L_real % _Rc_real).str();
// info += "o Reciprocal-space basis: ";
// info += (boost::format("N=%1% L=%2% Rc=%3%\n") % _N_recip % _L_recip % _Rc_recip).str();
// info += "o Reciprocal-space lattice: ";
// info += (boost::format("b1x=%1% b1y=%2% b1z=%3% ") % _b1.x() % _b1.y() % _b1.z()).str();
// info += (boost::format("b2x=%1% b2y=%2% b2z=%3% ") % _b2.x() % _b2.y() % _b2.z()).str();
// info += (boost::format("b3x=%1% b3y=%2% b3z=%3% ") % _b3.x() % _b3.y() % _b3.z()).str();
std::map<std::string, std::string>::iterator it;
for (it = _key_value_map.begin(); it != _key_value_map.end(); ) {
info += (boost::format(" o %1$-30s : %2$s") % it->first % it->second).str();
if (++it != _key_value_map.end()) info += "\n";
}
return info;
}
static void registerPython();
......
src/radialbasis.cpp
View file @ 59cc0eed
......@@ -18,8 +18,9 @@ namespace ub = boost::numeric::ublas;
void RadialBasis::configure(Options &options) {
_N = options.get<int>("radialbasis.N");
_Rc = options.get<double>("radialbasis.Rc");
_Rc = options.get<double>("radialcutoff.Rc");
_integration_steps = options.get<int>("radialbasis.integration_steps");
_mode = options.get<std::string>("radialbasis.mode");
}
RadialCoefficients RadialBasis::computeCoefficients(double r) {
......@@ -51,9 +52,7 @@ void RadialBasisGaussian::configure(Options &options) {
// CREATE & STORE EQUISPACED RADIAL GAUSSIANS
this->clear();
std::string mode = "adaptive";
mode = "equispaced";
if (mode == "equispaced") {
if (_mode == "equispaced") {
double dr = _Rc/(_N-1);
for (int i = 0; i < _N; ++i) {
double r = i*dr;
......@@ -62,16 +61,33 @@ void RadialBasisGaussian::configure(Options &options) {
_basis.push_back(new_fct);
}
}
else if (mode == "adaptive") {
double delta = 0.5;
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);
double sigma_0 = 0.5;
double sigma_stride_factor = 0.5;
// 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;
// }
for (int i = 0; i < _N; ++i) {
//sigma = sqrt(4./(2*L+1))*(r+delta);
sigma = sqrt(4./(2*L+1)*r*r + sigma_0*sigma_0);
std::cout << r << " " << sigma << std::endl;
r = r + sigma_stride_factor*sigma;
}
double scale = 1.; //_Rc/(r-sigma);
r = 0;
for (int i = 0; i < _N; ++i) {
sigma = sqrt(4./(2*L+1)*r*r + sigma_0*sigma_0);
basis_fct_t *new_fct = new basis_fct_t(scale*r, sigma);
_basis.push_back(new_fct);
r = r + sigma;
std::cout << r << " " << sigma << std::endl;
r = r+sigma_stride_factor*sigma;
}
}
else {
......@@ -83,6 +99,7 @@ void RadialBasisGaussian::configure(Options &options) {
GLOG() << (*bit)->_r0 << " (" << (*bit)->_sigma << ") ";
}
GLOG() << "}" << std::endl;
// 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;
......@@ -108,15 +125,6 @@ void RadialBasisGaussian::configure(Options &options) {
_Sij(i,j) = s;
}
}
// ORTHONORMALIZATION VIA CHOLESKY DECOMPOSITION
_Uij = _Sij;
soap::linalg::linalg_cholesky_decompose(_Uij);
// ZERO UPPER DIAGONAL OF U
for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i)
for (jt = it+1, j = i+1; jt != _basis.end(); ++jt, ++j)
_Uij(i,j) = 0.0;
_Tij = _Uij;
soap::linalg::linalg_invert(_Uij, _Tij);
// REPORT
GLOG() << "Radial basis overlap matrix" << std::endl;
for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i) {
......@@ -125,13 +133,26 @@ void RadialBasisGaussian::configure(Options &options) {
}
GLOG() << std::endl;
}
// ORTHONORMALIZATION VIA CHOLESKY DECOMPOSITION
_Uij = _Sij;
soap::linalg::linalg_cholesky_decompose(_Uij);
// REPORT
GLOG() << "Radial basis Cholesky decomposition" << std::endl;
for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i) {
for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i) {
for (jt = _basis.begin(), j = 0; jt != _basis.end(); ++jt, ++j) {
GLOG() << boost::format("%1$+1.4e") % _Uij(i,j) << " " << std::flush;
}
GLOG() << std::endl;
}
// ZERO UPPER DIAGONAL OF U
for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i)
for (jt = it+1, j = i+1; jt != _basis.end(); ++jt, ++j)
_Uij(i,j) = 0.0;
_Tij = _Uij;
soap::linalg::linalg_invert(_Uij, _Tij);
// REPORT
GLOG() << "Radial basis transformation matrix" << std::endl;
for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i) {
for (jt = _basis.begin(), j = 0; jt != _basis.end(); ++jt, ++j) {
......
src/radialbasis.hpp
View file @ 59cc0eed
......@@ -51,6 +51,7 @@ protected:
int _N;
double _Rc;
int _integration_steps;
std::string _mode; // <- 'equispaced' or 'adaptive'
static const double RADZERO = 1e-10;
};
......
src/spectrum.cpp
View file @ 59cc0eed
......@@ -17,6 +17,13 @@ Spectrum::~Spectrum() {
_log = NULL;
delete _basis;
_basis = NULL;
atomspec_array_t::iterator it;
for (it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
delete *it;
}
_atomspec_array.clear();
// AtomicSpectra in type map already deleted above, clear only:
_map_atomspec_array.clear();
}
void Spectrum::compute() {
......@@ -24,17 +31,59 @@ void Spectrum::compute() {
GLOG() << _options->summarizeOptions() << std::endl;
GLOG() << "Using radial basis of type '" << _basis->getRadBasis()->identify() << "'" << std::endl;
GLOG() << "Using angular basis of type '" << _basis->getAngBasis()->identify() << "'" << std::endl;
GLOG() << "Using cutoff function of type '" << _basis->getCutoff()->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;
AtomicSpectrum *atomic_spectrum = this->computeAtomic(*pit);
//atomic_spectrum->getReduced()->writeDensityOnGrid("density.expanded.cube", _options, _structure, *pit, true);
//atomic_spectrum->getReduced()->writeDensityOnGrid("density.explicit.cube", _options, _structure, *pit, false);
this->add(atomic_spectrum);
}
}
void Spectrum::writeDensityOnGrid(int slot_idx, std::string center_type, std::string density_type) {
AtomicSpectrum *atomic_spectrum = NULL;
// FIND SPECTRUM
if (center_type == "") {
// NO TYPE => ACCESS ARRAY
// Slot index valid?
if (slot_idx < _atomspec_array.size()) {
atomic_spectrum = _atomspec_array[slot_idx];
}
else {
throw soap::base::OutOfRange("Spectrum slot index");
}
}
else {
// TYPE => ACCESS TYPE MAP
map_atomspec_array_t::iterator it = _map_atomspec_array.find(center_type);
// Any such type?
if (it == _map_atomspec_array.end()) {
throw soap::base::OutOfRange("No spectrum of type '" + center_type + "'");
}
else {
// Slot index valid?
if (slot_idx < it->second.size()) {
atomic_spectrum = it->second[slot_idx];
}
else {
throw soap::base::OutOfRange("Spectrum slot index, type '" + center_type + "'");
}
}
}
// WRITE CUBE FILES
if (atomic_spectrum) {
atomic_spectrum->getExpansion(density_type)->writeDensityOnGrid(
"density.expanded.cube", _options, _structure, atomic_spectrum->getCenter(), true);
atomic_spectrum->getExpansion(density_type)->writeDensityOnGrid(
"density.explicit.cube", _options, _structure, atomic_spectrum->getCenter(), false);
}
return;
}
// TODO change to ::computeAtomic(Particle *center, vector<Particle*> &nbhood)
void Spectrum::computeAtomic(Particle *center) {
AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
GLOG() << "Compute atomic spectrum for particle " << center->getId()
<< " (type " << center->getType() << ")" << std::endl;
......@@ -46,21 +95,31 @@ void Spectrum::computeAtomic(Particle *center) {
// BasisCoefficients c_nlm(c_n_zero, c_lm_zero);
// c_nlm.linkBasis(_radbasis, _angbasis);
BasisExpansion *nbhood_expansion = new BasisExpansion(this->_basis);
// BasisExpansion *nbhood_expansion = new BasisExpansion(this->_basis);
AtomicSpectrum *atomic_spectrum = new AtomicSpectrum(center, this->_basis);
Structure::particle_it_t pit;
for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
// TODO Cut-off check && cut-off smoothing (= weight reduction)
// FIND DISTANCE & DIRECTION, APPLY CUTOFF
vec dr = _structure->connect(center->getPos(), (*pit)->getPos());
double r = soap::linalg::abs(dr);
double weight_scale = this->_basis->getCutoff()->calculateWeight(r);
if (weight_scale < 0.) continue;
vec d = dr/r;
std::string type_other = (*pit)->getType();
// APPLY WEIGHT IF CENTER
if (*pit == center) {
weight_scale *= this->_basis->getCutoff()->getCenterWeight();
}
// COMPUTE EXPANSION & ADD TO SPECTRUM
BasisExpansion nb_expansion(this->_basis);
nb_expansion.computeCoefficients(r, d, (*pit)->getWeight(), (*pit)->getSigma());
nbhood_expansion->add(nb_expansion);
nb_expansion.computeCoefficients(r, d, weight_scale*(*pit)->getWeight(), (*pit)->getSigma());
std::string type_other = (*pit)->getType();
atomic_spectrum->add(type_other, nb_expansion);
// nbhood_expansion->add(nb_expansion);
// // COMPUTE RADIAL COEFFICIENTS
......@@ -93,10 +152,26 @@ void Spectrum::computeAtomic(Particle *center) {
// c_nlm.writeDensityOnGrid("density.expanded.cube", _options, _structure, center, true);
// c_nlm.writeDensityOnGrid("density.explicit.cube", _options, _structure, center, false);
nbhood_expansion->writeDensityOnGrid("density.expanded.cube", _options, _structure, center, true);
nbhood_expansion->writeDensityOnGrid("density.explicit.cube", _options, _structure, center, false);
// nbhood_expansion->writeDensityOnGrid("density.expanded.cube", _options, _structure, center, true);
// nbhood_expansion->writeDensityOnGrid("density.explicit.cube", _options, _structure, center, false);
return atomic_spectrum;
}
delete nbhood_expansion;
void Spectrum::add(AtomicSpectrum *atomspec) {
assert(atomspec->getBasis() == _basis &&
"Should not append atomic spectrum linked against different basis.");
std::string atomspec_type = atomspec->getCenterType();
map_atomspec_array_t::iterator it = _map_atomspec_array.find(atomspec_type);
if (it == _map_atomspec_array.end()) {
_map_atomspec_array[atomspec_type] = atomspec_array_t();
it = _map_atomspec_array.find(atomspec_type);
}
it->second.push_back(atomspec);
_atomspec_array.push_back(atomspec);
return;
}
......@@ -104,7 +179,8 @@ void Spectrum::registerPython() {
using namespace boost::python;
class_<Spectrum>("Spectrum", init<Structure &, Options &>())
.def("compute", &Spectrum::compute)
.def("saveAndClean", &Spectrum::saveAndClean);
.def("saveAndClean", &Spectrum::saveAndClean)
.def("writeDensityOnGrid", &Spectrum::writeDensityOnGrid);
}
/* STORAGE, BASIS, COMPUTATION, PARALLELIZATION */
......
src/spectrum.hpp
View file @ 59cc0eed
#ifndef _SOAP_SPECTRUM_HPP
#define _SOAP_SPECTRUM_HPP
#include "options.hpp"
#include "structure.hpp"
#include <assert.h>
#include "base/logger.hpp"
#include "types.hpp"
#include "globals.hpp"
#include "radialbasis.hpp"
#include "angularbasis.hpp"
#include "options.hpp"
#include "structure.hpp"
#include "basis.hpp"
#include "base/logger.hpp"
namespace soap {
class AtomicSpectrum : public std::map<std::string, BasisExpansion*>
{
public:
AtomicSpectrum(Particle *center, Basis *basis) :
_center(center),
_center_pos(center->getPos()),
_center_type(center->getType()),
_basis(basis) {
_expansion_reduced = new BasisExpansion(_basis);
}
~AtomicSpectrum() {
iterator it;
for (it = this->begin(); it != this->end(); ++it) delete it->second;
this->clear();
delete _expansion_reduced;
}
void add(std::string type, BasisExpansion &nb_expansion) {
assert(nb_expansion.getBasis() == _basis &&
"Should not sum expansions linked against different bases.");
iterator it = this->find(type);
if (it == this->end()) {
(*this)[type] = new BasisExpansion(_basis);
it = this->find(type);
}
it->second->add(nb_expansion);
_expansion_reduced->add(nb_expansion);
return;
}
Particle *getCenter() { return _center; }
std::string &getCenterType() { return _center_type; }
vec &getCenterPos() { return _center_pos; }
BasisExpansion *getReduced() { return _expansion_reduced; }
BasisExpansion *getExpansion(std::string type) {
if (type == "") {
return _expansion_reduced;
}
iterator it = this->find(type);
if (it == this->end()) {
throw soap::base::OutOfRange("AtomicSpectrum: No such type '" + type + "'");
return NULL;
}
else {
return it->second;
}
}
Basis *getBasis() { return _basis; }
protected:
Particle *_center;
vec _center_pos;
std::string _center_type;
Basis *_basis;
BasisExpansion *_expansion_reduced;
//std::map<std::string, BasisExpansion*> _map_type_expansion;
};
class CenterDensity
{
CenterDensity() {}
......@@ -39,6 +97,9 @@ class Target
class Spectrum
{
public:
typedef std::vector<AtomicSpectrum*> atomspec_array_t;
typedef std::map<std::string, atomspec_array_t> map_atomspec_array_t;
Spectrum(Structure &structure, Options &options);
~Spectrum();
......@@ -47,7 +108,10 @@ public:
void clean();
void compute();
void computeAtomic(Particle *center);
AtomicSpectrum *computeAtomic(Particle *center);
void add(AtomicSpectrum *atomspec);
void writeDensityOnGrid(int slot_idx, std::string center_type, std::string density_type);
void computePower();
void computeLinear();
......@@ -59,13 +123,13 @@ private:
Options *_options;
Structure *_structure;
Basis *_basis;
atomspec_array_t _atomspec_array;
map_atomspec_array_t _map_atomspec_array;
};
class AtomicSpectrum
{
AtomicSpectrum() {}
};
......
test/test.py
View file @ 59cc0eed
......@@ -20,20 +20,29 @@ element_mass = {
}
ase_config_list = soap.tools.ase_load_all('configs')
config = ase_config_list[0]
config = ase_config_list[1]
osio << config.atoms << endl
# INITIALIZE OPTIONS
options = soap.Options()
options.set('radialbasis.type', 'gaussian')
options.set('radialbasis.mode', 'adaptive')
#options.set('radialbasis.mode', 'equispaced')
options.set('radialbasis.N', 9)
options.set('radialbasis.Rc', 4.0)
options.set('radialbasis.sigma', 0.5)
options.set('radialbasis.integration_steps', 15)
#options.set('radialbasis.N', 9)
options.set('radialcutoff.Rc', 4.)
options.set('radialcutoff.Rc_width', 0.5)
options.set('radialcutoff.type', 'shifted-cosine')
options.set('radialcutoff.center_weight', -7.)
options.set('radialcutoff.center_weight', 1.)
options.set('angularbasis.type', 'spherical-harmonic')
options.set('angularbasis.L', 6)
#options.set('angularbasis.L', 6)
options.set('densitygrid.N', 20)
options.set('densitygrid.dx', 0.10)
options.set('densitygrid.dx', 0.15)
# LOAD STRUCTURE
structure = soap.tools.setup_structure_ase(config.config_file, config.atoms)
......@@ -41,13 +50,16 @@ structure = soap.tools.setup_structure_ase(config.config_file, config.atoms)
osio << osio.mg << structure.label << endl
for atom in structure:
atom.sigma = 0.4*element_vdw_radius[atom.type]
atom.sigma = 0.5 # 0.4*element_vdw_radius[atom.type]
atom.weight = 1. #element_mass[atom.type]
#atom.sigma = 0.5*element_vdw_radius[atom.type]
#atom.weight = element_mass[atom.type]
osio << atom.name << atom.type << atom.weight << atom.sigma << atom.pos << endl
# COMPUTE SPECTRUM
spectrum = soap.Spectrum(structure, options)
spectrum.compute()
spectrum.writeDensityOnGrid(1, "C", "")
......
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