Commit 07b25c69 authored by Carl Poelking's avatar Carl Poelking
Browse files

AtomicSpectrum::computerPower() and calls.

parent ea2ca74a
......@@ -7,5 +7,107 @@
namespace soap {
AtomicSpectrum::AtomicSpectrum(Particle *center, Basis *basis) {
this->null();
_center = center;
_center_pos = center->getPos();
_center_type = center->getType();
_basis = basis;
_qnlm_generic = new BasisExpansion(_basis);
}
AtomicSpectrum::~AtomicSpectrum() {
// Clean Qnlm's
for (map_qnlm_t::iterator it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) delete it->second;
_map_qnlm.clear();
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;
_map_xnkl.clear();
if (_xnkl_generic_coherent) {
delete _xnkl_generic_coherent;
_xnkl_generic_coherent = NULL;
}
if (_xnkl_generic_incoherent) {
delete _xnkl_generic_incoherent;
_xnkl_generic_incoherent = NULL;
}
}
void AtomicSpectrum::null() {
_center = NULL;
_center_pos = vec(0,0,0);
_center_type = "?";
_basis = NULL;
_qnlm_generic = NULL;
_xnkl_generic_coherent = NULL;
_xnkl_generic_incoherent = NULL;
}
void AtomicSpectrum::addQnlm(std::string type, qnlm_t &nb_expansion) {
assert(nb_expansion.getBasis() == _basis &&
"Should not sum expansions linked against different bases.");
map_qnlm_t::iterator it = _map_qnlm.find(type);
if (it == _map_qnlm.end()) {
_map_qnlm[type] = new BasisExpansion(_basis);
it = _map_qnlm.find(type);
}
it->second->add(nb_expansion);
_qnlm_generic->add(nb_expansion);
return;
}
AtomicSpectrum::qnlm_t *AtomicSpectrum::getQnlm(std::string type) {
if (type == "") {
return _qnlm_generic;
}
map_qnlm_t::iterator it = _map_qnlm.find(type);
if (it == _map_qnlm.end()) {
throw soap::base::OutOfRange("AtomicSpectrum: No such type '" + type + "'");
}
else {
return it->second;
}
}
void AtomicSpectrum::computePower() {
// Specific (i.e., type-dependent)
map_qnlm_t::iterator it1;
map_qnlm_t::iterator it2;
for (it1 = _map_qnlm.begin(); it1 != _map_qnlm.end(); ++it1) {
for (it2 = _map_qnlm.begin(); it2 != _map_qnlm.end(); ++it2) {
type_pair_t types(it1->first, it2->first);
GLOG() << " " << types.first << ":" << types.second << std::flush;
PowerExpansion *powex = new PowerExpansion(_basis);
powex->computeCoefficients(it1->second, it2->second);
_map_xnkl[types] = powex;
}
}
// Generic coherent
GLOG() << " g/c" << std::flush;
if (_xnkl_generic_coherent) delete _xnkl_generic_coherent;
_xnkl_generic_coherent = new PowerExpansion(_basis);
// Generic incoherent
GLOG() << " g/i" << std::flush;
if (_xnkl_generic_incoherent) delete _xnkl_generic_incoherent;
_xnkl_generic_incoherent = new PowerExpansion(_basis);
map_xnkl_t::iterator it;
for (it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) {
_xnkl_generic_incoherent->add(it->second);
}
GLOG() << std::endl;
}
AtomicSpectrum::xnkl_t *AtomicSpectrum::getXnkl(type_pair_t &types) {
map_xnkl_t::iterator it = _map_xnkl.find(types);
if (it == _map_xnkl.end()) {
throw soap::base::OutOfRange("AtomicSpectrum: No such type pair '" + types.first + ":" + types.second + "'");
}
return it->second;
}
}
......@@ -26,76 +26,27 @@ public:
typedef PowerExpansion xnkl_t;
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;
AtomicSpectrum(Particle *center, Basis *basis);
AtomicSpectrum() { this->null(); }
~AtomicSpectrum();
void null();
AtomicSpectrum(Particle *center, Basis *basis) :
_center(center),
_center_pos(center->getPos()),
_center_type(center->getType()),
_basis(basis),
_xnkl_generic_coherent(NULL),
_xnkl_generic_incoherent(NULL) {
_qnlm_generic = new BasisExpansion(_basis);
}
AtomicSpectrum() :
_center(NULL),
_center_pos(vec(0,0,0)),
_center_type("?"),
_basis(NULL),
_qnlm_generic(NULL),
_xnkl_generic_coherent(NULL),
_xnkl_generic_incoherent(NULL) { ; }
~AtomicSpectrum() {
// Clean Qnlm's
for (map_qnlm_t::iterator it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) delete it->second;
_map_qnlm.clear();
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;
_map_xnkl.clear();
if (_xnkl_generic_coherent) {
delete _xnkl_generic_coherent;
_xnkl_generic_coherent = NULL;
}
if (_xnkl_generic_incoherent) {
delete _xnkl_generic_incoherent;
_xnkl_generic_incoherent = NULL;
}
}
// CENTER & BASIS METHODS
Particle *getCenter() { return _center; }
std::string &getCenterType() { return _center_type; }
vec &getCenterPos() { return _center_pos; }
Basis *getBasis() { return _basis; }
// QNLM METHODS
void addQnlm(std::string type, qnlm_t &nb_expansion) {
assert(nb_expansion.getBasis() == _basis &&
"Should not sum expansions linked against different bases.");
iterator it = _map_qnlm.find(type);
if (it == _map_qnlm.end()) {
_map_qnlm[type] = new BasisExpansion(_basis);
it = _map_qnlm.find(type);
}
it->second->add(nb_expansion);
_qnlm_generic->add(nb_expansion);
return;
}
void addQnlm(std::string type, qnlm_t &nb_expansion);
qnlm_t *getQnlm(std::string type);
qnlm_t *getQnlmGeneric() { return _qnlm_generic; }
qnlm_t *getQnlm(std::string type) {
if (type == "") {
return _qnlm_generic;
}
iterator it = _map_qnlm.find(type);
if (it == _map_qnlm.end()) {
throw soap::base::OutOfRange("AtomicSpectrum: No such type '" + type + "'");
return NULL;
}
else {
return it->second;
}
}
// XNLM METHODS
void computePower();
xnkl_t *getXnkl(type_pair_t &types);
xnkl_t *getXnklGenericCoherent() { return _xnkl_generic_coherent; }
xnkl_t *getXnklGenericIncoherent() { return _xnkl_generic_incoherent; }
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
......
......@@ -11,7 +11,7 @@ void CutoffFunction::configure(Options &options) {
GLOG() << "Weighting function with "
<< "Rc = " << _Rc
<< ", _Rc_width = " << _Rc_width
<< ", central weigth = " << _center_weight << std::endl;
<< ", central weight = " << _center_weight << std::endl;
}
void CutoffFunctionFactory::registerAll(void) {
......
......@@ -2,5 +2,40 @@
namespace soap {
PowerExpansion::PowerExpansion(Basis *basis) :
_basis(basis),
_L(basis->getAngBasis()->L()),
_N(basis->getRadBasis()->N()) {
_coeff = coeff_zero_t(_N*_N, _L+1);
}
void PowerExpansion::computeCoefficients(BasisExpansion *basex1, BasisExpansion *basex2) {
if (!_basis) throw soap::base::APIError("PowerExpansion::computeCoefficients, basis not initialised.");
BasisExpansion::coeff_t &coeff1 = basex1->getCoefficients();
BasisExpansion::coeff_t &coeff2 = basex2->getCoefficients();
for (int n = 0; n < _N; ++n) {
for (int k = 0; k < _N; ++k) {
for (int l = 0; l < (_L+1); ++l) {
//std::cout << n << " " << k << " " << l << " : " << std::flush;
std::complex<double> c_nkl = 0.0;
for (int m = -l; m <= l; ++m) {
//std::cout << m << " " << std::flush;
c_nkl += coeff1(n, l*l+l+m)*std::conj(coeff2(k, l*l+l+m));
}
_coeff(n*_N+k, l) = c_nkl;
//std::cout << std::endl;
}
}
}
//throw soap::base::APIError("");
return;
}
void PowerExpansion::add(PowerExpansion *other) {
assert(other->_basis == _basis &&
"Should not sum expansions linked against different bases.");
_coeff = _coeff + other->_coeff;
return;
}
}
......@@ -20,20 +20,24 @@ public:
typedef ub::matrix< std::complex<double> > coeff_t;
typedef ub::zero_matrix< std::complex<double> > coeff_zero_t;
PowerExpansion() {
;
}
PowerExpansion(BasisExpansion *basex) {
;
}
PowerExpansion() : _basis(NULL), _L(-1), _N(-1) {;}
PowerExpansion(Basis *basis);
void computeCoefficients(BasisExpansion *basex1, BasisExpansion *basex2);
void add(PowerExpansion *other);
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
;
arch & _basis;
arch & _N;
arch & _L;
arch & _coeff;
}
private:
Basis *_basis;
int _N;
int _L;
coeff_t _coeff; // access via (N*n+k, l) with shape (N*N, L+1)
};
......
......@@ -89,6 +89,14 @@ void Spectrum::writeDensityOnGrid(int slot_idx, std::string center_type, std::st
return;
}
void Spectrum::computePower() {
atomspec_array_t::iterator it;
for (it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
(*it)->computePower();
}
return;
}
// TODO change to ::computeAtomic(Particle *center, vector<Particle*> &nbhood)
AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
GLOG() << "Compute atomic spectrum for particle " << center->getId()
......@@ -198,6 +206,7 @@ void Spectrum::registerPython() {
class_<Spectrum>("Spectrum", init<Structure &, Options &>())
.def(init<std::string>())
.def("compute", &Spectrum::compute)
.def("computePower", &Spectrum::computePower)
.def("saveAndClean", &Spectrum::saveAndClean)
.def("save", &Spectrum::save)
.def("load", &Spectrum::load)
......
......@@ -65,12 +65,14 @@ for atom in structure:
# COMPUTE SPECTRUM
spectrum = soap.Spectrum(structure, options)
spectrum.compute()
spectrum.writeDensityOnGrid(1, "C", "")
spectrum.computePower()
spectrum.save("test_serialization/%s.spectrum.arch" % structure.label)
#spectrum.writeDensityOnGrid(1, "C", "")
#spectrum.writeDensityOnGrid(2, "S", "")
#spectrum.writeDensityOnGrid(7, "C", "") # line.xyz
#spectrum.writeDensityOnGrid(3, "C", "") # linedot.xyz
#spectrum.writeDensityOnGrid(41, "C", "") # C60_pair.xyz
spectrum.save("test_serialization/%s.spectrum.arch" % structure.label)
osio.okquit()
......
......@@ -9,7 +9,7 @@ from momo import osio, endl, flush
archfile = 'config_000057.xyz.spectrum.arch'
spectrum = soap.Spectrum(archfile)
spectrum.writeDensityOnGrid(1, "C", "")
spectrum.writeDensityOnGrid(1, "C", "H")
spectrum.save("%s-2" % archfile)
osio.okquit()
......
Markdown is supported
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