Commit 7f3aa2d8 authored by Carl Poelking's avatar Carl Poelking
Browse files

Global coherent superposition spectrum.

parent 088d6352
......@@ -10,12 +10,19 @@ namespace soap {
AtomicSpectrum::AtomicSpectrum(Particle *center, Basis *basis) {
this->null();
_center = center;
_center_id = center->getId();
_center_pos = center->getPos();
_center_type = center->getType();
_basis = basis;
_qnlm_generic = new BasisExpansion(_basis);
}
AtomicSpectrum::AtomicSpectrum(Basis *basis) {
this->null();
_basis = basis;
_qnlm_generic = new BasisExpansion(_basis);
}
AtomicSpectrum::~AtomicSpectrum() {
// CLEAN QNLM'S
// ... Summed density expansions
......@@ -60,6 +67,7 @@ AtomicSpectrum::~AtomicSpectrum() {
void AtomicSpectrum::null() {
_center = NULL;
_center_id = -1;
_center_pos = vec(0,0,0);
_center_type = "?";
_basis = NULL;
......@@ -151,6 +159,49 @@ void AtomicSpectrum::addQnlmNeighbour(Particle *nb, qnlm_t *nb_expansion) {
return;
}
void AtomicSpectrum::mergeQnlm(AtomicSpectrum *other) {
// Function used to construct global spectrum as sum over atomic spectra.
// The result is itself an "atomic" spectrum (as data fields are largely identical,
// except for the fact that this summed spectrum does not have a well-defined center.
assert(other->getBasis() == _basis &&
"Should not merge atomic spectra linked against different bases.");
// Type-agnostic (=generic) density expansion
_qnlm_generic->add(*other->getQnlmGeneric());
// Type-resolved (=specific) density expansions
map_qnlm_t &map_qnlm_other = other->getQnlmMap();
for (auto it = map_qnlm_other.begin(); it != map_qnlm_other.end(); ++it) {
std::string density_type = it->first;
BasisExpansion *density = it->second;
// Already have density of this type?
auto mit = _map_qnlm.find(density_type);
if (mit == _map_qnlm.end()) {
_map_qnlm[density_type] = new qnlm_t(_basis);
mit = _map_qnlm.find(density_type);
}
// Add ...
mit->second->add(*density);
}
// Particle-ID-resolved gradients
map_pid_qnlm_t &map_pid_qnlm_other = other->getPidQnlmMap();
for (auto it = map_pid_qnlm_other.begin(); it != map_pid_qnlm_other.end(); ++it) {
int pid = it->first;
std::string pid_type = it->second.first;
qnlm_t *density_grad = it->second.second;
// Already have density gradient for this particle?
auto mit = _map_pid_qnlm.find(pid);
if (mit == _map_pid_qnlm.end()) {
qnlm_t *qnlm = new qnlm_t(_basis);
// Remember to setup zero matrices to store gradient ...
qnlm->zeroGradient();
_map_pid_qnlm[pid] = std::pair<std::string,qnlm_t*>(pid_type, qnlm);
mit = _map_pid_qnlm.find(pid);
}
// Add gradients ...
mit->second.second->addGradient(*density_grad);
}
return;
}
AtomicSpectrum::qnlm_t *AtomicSpectrum::getQnlm(std::string type) {
if (type == "") {
return _qnlm_generic;
......@@ -165,7 +216,7 @@ AtomicSpectrum::qnlm_t *AtomicSpectrum::getQnlm(std::string type) {
}
void AtomicSpectrum::computePowerGradients() {
GLOG() << "CID " << _center->getId() << ": " << std::endl;
GLOG() << "CID " << _center_id << ": " << std::endl;
for (auto it1 = _map_pid_qnlm.begin(); it1 != _map_pid_qnlm.end(); ++it1) {
// Derivatives are taken with respect to the coordinates of this neighbour particle
int pi_id = it1->first;
......
......@@ -38,6 +38,7 @@ public:
AtomicSpectrum(Particle *center, Basis *basis);
AtomicSpectrum(Basis *basis);
AtomicSpectrum() { this->null(); }
~AtomicSpectrum();
void null();
......@@ -55,6 +56,9 @@ public:
void addQnlmNeighbour(Particle *nb, qnlm_t *nb_expansion);
qnlm_t *getQnlm(std::string type);
qnlm_t *getQnlmGeneric() { return _qnlm_generic; }
map_qnlm_t &getQnlmMap() { return _map_qnlm; }
map_pid_qnlm_t &getPidQnlmMap() { return _map_pid_qnlm; }
void mergeQnlm(AtomicSpectrum *other);
// XNKL METHODS
void computePower();
void computePowerGradients();
......@@ -73,6 +77,7 @@ public:
void serialize(Archive &arch, const unsigned int version) {
// Center & basis
arch & _center;
arch & _center_id;
arch & _center_pos;
arch & _center_type;
arch & _basis;
......@@ -92,6 +97,7 @@ public:
protected:
// CENTER & BASIS LINKS
Particle *_center;
int _center_id;
vec _center_pos;
std::string _center_type;
Basis *_basis;
......
......@@ -165,6 +165,24 @@ void BasisExpansion::computeCoefficients(double r, vec d, double weight, double
return;
}
void BasisExpansion::addGradient(BasisExpansion &other) {
assert(_has_gradients);
_coeff_grad_x = _coeff_grad_x + other.getCoefficientsGradX();
_coeff_grad_y = _coeff_grad_y + other.getCoefficientsGradY();
_coeff_grad_z = _coeff_grad_z + other.getCoefficientsGradZ();
return;
}
void BasisExpansion::zeroGradient() {
_has_gradients = true;
int L = _angbasis->L();
int N = _radbasis->N();
_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));
return;
}
void BasisExpansion::conjugate() {
for (int n = 0; n != _radbasis->N(); ++n) {
for (int l = 0; l != _angbasis->L()+1; ++l) {
......
......@@ -73,6 +73,8 @@ public:
bool hasScalars() { return _has_scalars; }
bool hasGradients() { return _has_gradients; }
void add(BasisExpansion &other) { _coeff = _coeff + other._coeff; }
void addGradient(BasisExpansion &other);
void zeroGradient();
void conjugate();
void writeDensity(std::string filename, Options *options,
Structure *structure, Particle *center);
......
......@@ -8,19 +8,19 @@
namespace soap {
Spectrum::Spectrum(Structure &structure, Options &options) :
_log(NULL), _options(&options), _structure(&structure), _own_basis(true) {
_log(NULL), _options(&options), _structure(&structure), _own_basis(true), _global_atomic(NULL) {
GLOG() << "Configuring spectrum ..." << std::endl;
// CREATE & CONFIGURE BASIS
_basis = new Basis(&options);
}
Spectrum::Spectrum(Structure &structure, Options &options, Basis &basis) :
_log(NULL), _options(&options), _structure(&structure), _basis(&basis), _own_basis(false) {
_log(NULL), _options(&options), _structure(&structure), _basis(&basis), _own_basis(false), _global_atomic(NULL) {
;
}
Spectrum::Spectrum(std::string archfile) :
_log(NULL), _options(NULL), _structure(NULL), _basis(NULL), _own_basis(true) {
_log(NULL), _options(NULL), _structure(NULL), _basis(NULL), _own_basis(true), _global_atomic(NULL) {
this->load(archfile);
}
......@@ -38,6 +38,9 @@ Spectrum::~Spectrum() {
_atomspec_array.clear();
// AtomicSpectra in type map already deleted above, clear only:
_map_atomspec_array.clear();
if (_global_atomic) delete _global_atomic;
_global_atomic = NULL;
}
void Spectrum::compute() {
......@@ -122,6 +125,20 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
return atomic_spectrum;
}
AtomicSpectrum *Spectrum::computeGlobal() {
if (_global_atomic) throw soap::base::APIError("<Spectrum::computeGlobal> Already initialised.");
_global_atomic = new AtomicSpectrum(_basis);
GLOG() << "Computing global spectrum ..." << std::endl;
for (auto it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
GLOG() << " Adding center " << (*it)->getCenter()->getId()
<< " (type " << (*it)->getCenter()->getType() << ")" << std::endl;
_global_atomic->mergeQnlm(*it);
}
_global_atomic->computePower();
_global_atomic->computePowerGradients();
return _global_atomic;
}
AtomicSpectrum *Spectrum::getAtomic(int slot_idx, std::string center_type) {
AtomicSpectrum *atomic_spectrum = NULL;
// FIND SPECTRUM
......@@ -241,8 +258,10 @@ void Spectrum::registerPython() {
.def("compute", computeCentersTargets)
.def("computePower", &Spectrum::computePower)
.def("computePowerGradients", &Spectrum::computePowerGradients)
.def("computeGlobal", &Spectrum::computeGlobal, return_value_policy<reference_existing_object>())
.def("addAtomic", &Spectrum::addAtomic)
.def("getAtomic", &Spectrum::getAtomic, return_value_policy<reference_existing_object>())
.def("getGlobal", &Spectrum::getGlobal, return_value_policy<reference_existing_object>())
.def("saveAndClean", &Spectrum::saveAndClean)
.def("save", &Spectrum::save)
.def("load", &Spectrum::load)
......
......@@ -50,8 +50,10 @@ public:
void compute(Structure::particle_array_t &sources, Structure::particle_array_t &targets);
AtomicSpectrum *computeAtomic(Particle *center);
AtomicSpectrum *computeAtomic(Particle *center, Structure::particle_array_t &targets);
AtomicSpectrum *computeGlobal();
void addAtomic(AtomicSpectrum *atomspec);
AtomicSpectrum *getAtomic(int slot_idx, std::string center_type);
AtomicSpectrum *getGlobal() { assert(_global_atomic && "Compute first"); return _global_atomic; }
void writeDensityOnGrid(int slot_idx, std::string center_type, std::string density_type);
void writeDensityOnGridInverse(int slot_idx, std::string center_type, std::string type1, std::string type2);
void writeDensity(int slot_idx, std::string center_type, std::string density_type);
......@@ -71,6 +73,8 @@ public:
arch & _atomspec_array;
arch & _map_atomspec_array;
arch & _global_atomic;
return;
}
......@@ -84,6 +88,8 @@ private:
atomspec_array_t _atomspec_array;
map_atomspec_array_t _map_atomspec_array;
AtomicSpectrum *_global_atomic;
};
......
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