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

Center-target spectra.

parent d36596b8
Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ a6b5ebbd
......@@ -9,3 +9,4 @@ test/back
*.coeff
*.cube
*.arch
*.back
src/spectrum.cpp
View file @ a6b5ebbd
......@@ -41,22 +41,118 @@ Spectrum::~Spectrum() {
}
void Spectrum::compute() {
GLOG() << "Compute spectrum ..." << std::endl;
this->compute(_structure->particles(), _structure->particles());
}
void Spectrum::compute(Segment *center, Segment *target) {
this->compute(center->particles(), target->particles());
}
void Spectrum::compute(Structure::particle_array_t &centers, Structure::particle_array_t &targets) {
GLOG() << "Compute spectrum "
<< "(centers " << centers.size() << ", targets " << targets.size() << ") ..." << std::endl;
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) {
for (pit = centers.begin(); pit != centers.end(); ++pit) {
// Continue if exclusion defined ...
if (_options->doExcludeCenter((*pit)->getType())) continue;
// Compute ...
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->addAtomic(atomic_spectrum);
AtomicSpectrum *atomic_spectrum = this->computeAtomic(*pit, targets);
//atomic_spectrum->getReduced()->writeDensityOnGrid("density.expanded.cube", _options, _structure, *pit, true);
//atomic_spectrum->getReduced()->writeDensityOnGrid("density.explicit.cube", _options, _structure, *pit, false);
this->addAtomic(atomic_spectrum);
}
}
void Spectrum::computePower() {
atomspec_array_t::iterator it;
for (it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
(*it)->computePower();
}
return;
}
AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
return this->computeAtomic(center, _structure->particles());
}
AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_array_t &targets) {
GLOG() << "Compute atomic spectrum for particle " << center->getId()
<< " (type " << center->getType() << ", targets " << targets.size() << ") ..." << std::endl;
// 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);
AtomicSpectrum *atomic_spectrum = new AtomicSpectrum(center, this->_basis);
Structure::particle_it_t pit;
for (pit = targets.begin(); pit != targets.end(); ++pit) {
// CHECK FOR EXCLUSIONS
if (_options->doExcludeTarget((*pit)->getType())) continue;
// FIND DISTANCE & DIRECTION, APPLY CUTOFF (= WEIGHT REDUCTION)
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; // <- Negative cutoff weight means: skip
vec d = dr/r;
// 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, weight_scale*(*pit)->getWeight(), (*pit)->getSigma());
std::string type_other = (*pit)->getType();
atomic_spectrum->addQnlm(type_other, nb_expansion);
// 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);
}
// 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);
return atomic_spectrum;
}
AtomicSpectrum *Spectrum::getAtomic(int slot_idx, std::string center_type) {
......@@ -135,93 +231,6 @@ void Spectrum::writePowerDensity(int slot_idx, std::string center_type, std::str
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()
<< " (type " << center->getType() << ")" << std::endl;
// 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);
AtomicSpectrum *atomic_spectrum = new AtomicSpectrum(center, this->_basis);
Structure::particle_it_t pit;
for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
// CHECK FOR EXCLUSIONS
if (_options->doExcludeTarget((*pit)->getType())) continue;
// FIND DISTANCE & DIRECTION, APPLY CUTOFF (= WEIGHT REDUCTION)
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; // <- Negative cutoff weight means: skip
vec d = dr/r;
// 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, weight_scale*(*pit)->getWeight(), (*pit)->getSigma());
std::string type_other = (*pit)->getType();
atomic_spectrum->addQnlm(type_other, nb_expansion);
// 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);
}
// 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);
return atomic_spectrum;
}
void Spectrum::addAtomic(AtomicSpectrum *atomspec) {
assert(atomspec->getBasis() == _basis &&
"Should not append atomic spectrum linked against different basis.");
......@@ -252,11 +261,17 @@ void Spectrum::load(std::string archfile) {
void Spectrum::registerPython() {
using namespace boost::python;
void (Spectrum::*computeAll)() = &Spectrum::compute;
void (Spectrum::*computeSegPair)(Segment*, Segment*) = &Spectrum::compute;
void (Spectrum::*computeCentersTargets)(Structure::particle_array_t&, Structure::particle_array_t&) = &Spectrum::compute;
class_<Spectrum>("Spectrum", init<Structure &, Options &>())
.def(init<Structure &, Options &, Basis &>())
.def(init<std::string>())
.def("__iter__", range<return_value_policy<reference_existing_object> >(&Spectrum::beginAtomic, &Spectrum::endAtomic))
.def("compute", &Spectrum::compute)
.def("compute", computeAll)
.def("compute", computeSegPair)
.def("compute", computeCentersTargets)
.def("computePower", &Spectrum::computePower)
.def("addAtomic", &Spectrum::addAtomic)
.def("getAtomic", &Spectrum::getAtomic, return_value_policy<reference_existing_object>())
......
src/spectrum.hpp
View file @ a6b5ebbd
......@@ -46,7 +46,10 @@ public:
void clean();
void compute();
void compute(Segment *centers, Segment *targets);
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);
void addAtomic(AtomicSpectrum *atomspec);
AtomicSpectrum *getAtomic(int slot_idx, std::string center_type);
void writeDensityOnGrid(int slot_idx, std::string center_type, std::string density_type);
......
src/structure.cpp
View file @ a6b5ebbd
......@@ -47,7 +47,11 @@ Particle &Segment::addParticle(Particle *new_part) {
void Segment::registerPython() {
using namespace boost::python;
class_<Segment>("Segment", init<int>())
class_<Segment, Segment*>("Segment", init<int>())
.add_property("id", &Segment::getId, &Segment::setId)
.add_property("name", make_function(&Segment::getName, copy_non_const()), &Segment::setName)
.add_property("type", make_function(&Segment::getType, copy_non_const()), &Segment::setType)
.add_property("particles", range<return_value_policy<reference_existing_object> >(&Segment::beginParticles, &Segment::endParticles))
.def("addParticle", &Segment::addParticle, return_value_policy<reference_existing_object>());
}
......@@ -154,14 +158,18 @@ void Structure::registerPython() {
using namespace boost::python;
class_<Structure>("Structure", init<std::string>())
.def("addSegment", &Structure::addSegment, return_value_policy<reference_existing_object>())
.def("getSegment", &Structure::getSegment, return_value_policy<reference_existing_object>())
.def("addParticle", &Structure::addParticle, return_value_policy<reference_existing_object>())
.def("__iter__", range<return_value_policy<reference_existing_object> >(&Structure::beginParticles, &Structure::endParticles))
.add_property("particles", range<return_value_policy<reference_existing_object> >(&Structure::beginParticles, &Structure::endParticles))
.add_property("segments", range<return_value_policy<reference_existing_object> >(&Structure::beginSegments, &Structure::endSegments))
.def("connect", &Structure::connectNumeric)
.add_property("box", &Structure::getBoundaryNumeric, &Structure::setBoundaryNumeric)
.add_property("label", make_function(&Structure::getLabel, copy_non_const()), &Structure::setLabel);
class_<particle_array_t>("ParticleContainer")
.def(vector_indexing_suite<particle_array_t>());
class_<segment_array_t>("SegmentContainer")
.def(vector_indexing_suite<segment_array_t>());
}
} /* CLOSE NAMESPACE */
src/structure.hpp
View file @ a6b5ebbd
......@@ -6,6 +6,8 @@
#include <boost/python/iterator.hpp>
#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
#include <iostream>
#include "base/exceptions.hpp"
#include "types.hpp"
#include "boundary.hpp"
......@@ -49,6 +51,8 @@ public:
// Sigma
void setSigma(double sigma) { _sigma = sigma; }
double &getSigma() { return _sigma; }
// Segment
Segment *getSegment() { return _segment; }
static void registerPython();
......@@ -85,24 +89,43 @@ class Segment
{
public:
typedef std::vector<Particle*> particle_array_t;
typedef particle_array_t::iterator particle_it_t;
Segment(int id) : _id(id) { ; }
Segment() : _id(-1) { ; }
Segment(int id) : _id(id), _name("?"), _type("?") { ; }
Segment() : _id(-1), _name("?"), _type("?") { ; }
~Segment() { _particles.clear(); }
// PARTICLE CONTAINER
Particle &addParticle(Particle *new_part);
particle_array_t &particles() { return _particles; }
particle_it_t beginParticles() { return _particles.begin(); }
particle_it_t endParticles() { return _particles.end(); }
// Name
void setName(std::string name) { _name = name; }
std::string &getName() { return _name; }
// Id
void setId(int id) { _id = id; }
int getId() { return _id; }
// Type
void setType(std::string type) { _type = type; }
std::string &getType() { return _type; }
static void registerPython();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
arch & _id;
arch & _name;
arch & _type;
arch & _particles;
return;
}
private:
int _id;
std::string _name;
std::string _type;
particle_array_t _particles;
};
......@@ -120,10 +143,18 @@ public:
~Structure();
void null();
// PARTICLE CONTAINER
particle_array_t &particles() { return _particles; }
particle_it_t beginParticles() { return _particles.begin(); }
particle_it_t endParticles() { return _particles.end(); }
// SEGMENT CONTAINER
segment_array_t &segments() { return _segments; }
segment_it_t beginSegments() { return _segments.begin(); }
segment_it_t endSegments() { return _segments.end(); }
Segment *getSegment(int id) { if (id > _segments.size()) throw soap::base::OutOfRange("getSegment"); return _segments[id-1]; }
std::string &getLabel() { return _label; }
void setLabel(std::string label) { _label = label; }
......
src/tools/loadwrite.py
View file @ a6b5ebbd
......@@ -33,7 +33,7 @@ def ase_load_single(config_file, log=None):
ase_config_list = AseConfigList([config_file], log=log)
return ase_config_list[0]
def setup_structure_ase(label, ase_config):
def setup_structure_ase(label, ase_config, top=None):
# DEFINE SYSTEM
structure = soap.Structure(label)
# DEFINE BOUNDARY
......@@ -45,26 +45,42 @@ def setup_structure_ase(label, ase_config):
else:
raise NotImplementedError("<setup_structure_ase> Partial periodicity not implemented.")
structure.box = box
# CREATE SINGLE SEGMENT
segment = structure.addSegment()
# CREATE PARTICLES
# PARTICLE PROPERTIES
props = ['id', 'type', 'mass', 'pos']
ids = [ i+1 for i in range(ase_config.get_number_of_atoms()) ]
types = ase_config.get_atomic_numbers()
positions = ase_config.get_positions()
masses = ase_config.get_masses()
names = ase_config.get_chemical_symbols()
for id, name, typ, pos, mass in zip(ids, names, types, positions, masses):
particle = structure.addParticle(segment)
particle.pos = pos
particle.mass = mass
particle.weight = 1.
particle.sigma = 0.5
particle.name = name
particle.type = name
particle.type_id = typ
return structure
# DEFAULT TOPOLOGY
if top == None:
top = [('segment', 1, len(positions))]
# VERIFY PARTICLE COUNT
atom_count = 0
for section in top:
atom_count += section[1]*section[2]
assert atom_count == len(positions) # Does topology match structure?
# PARTITION: CREATE SEGMENTS, PARTICLES
atom_idx = 0
for section in top:
seg_type = section[0]
n_segs = section[1]
n_atoms = section[2]
for i in range(n_segs):
segment = structure.addSegment()
segment.name = seg_type
segment.type = seg_type
for j in range(n_atoms):
particle = structure.addParticle(segment)
particle.pos = positions[atom_idx]
particle.mass = masses[atom_idx]
particle.weight = 1.
particle.sigma = 0.5
particle.name = names[atom_idx]
particle.type = names[atom_idx]
particle.type_id = types[atom_idx]
atom_idx += 1
return structure
class AseConfigList(object):
def __init__(self,
......
test/test_particle_collection/pair_1_185.xyz 0 → 100644
View file @ a6b5ebbd
56
N 31.7089747 41.3358810 86.1879144
C 31.0963060 42.3161507 86.2185478
N 31.3970531 45.5996138 85.1468239
C 30.9257581 44.6682396 85.6434692
C 30.3374277 43.5219791 86.2588437
C 29.1051424 43.6238866 86.8452780
H 28.6444550 44.6060930 86.7995154
S 26.5114656 43.0332444 87.7821424
C 28.0557940 42.5877650 87.0792750
C 25.9764906 41.3711514 87.7813444
C 28.1114444 41.2134848 86.8713598
C 26.9523627 40.5322804 87.2623076
H 28.9805352 40.7295230 86.4444793
H 26.8266698 39.4615986 87.1706784
S 23.9568568 39.5509530 87.2655544
C 24.6420806 40.8717467 88.1791861
C 22.4560907 39.6397130 88.1695727
C 23.7445178 41.3425147 89.1269129
C 22.5256926 40.6534842 89.1195166
H 23.9674447 42.1589141 89.8009676
H 21.7023754 40.8769762 89.7858237
N 22.8727617 36.4567146 86.2629242
C 22.0151208 37.0052059 86.8114155
N 18.5761651 36.7452785 87.1665305
C 19.6445382 37.1630458 87.3089482
C 20.9614142 37.6845811 87.4894200
C 21.1535122 38.7831753 88.2825193
H 20.2573205 39.1962221 88.7356038
N 19.1184299 44.4606155 84.8200553
C 19.8203483 43.7586971 85.4132258
N 21.9330886 43.9144593 88.1577179
C 21.3729007 43.4590333 87.2549388
C 20.6843079 42.8896901 86.1411975
C 20.8756715 41.5716588 85.8265692
H 21.5634575 41.0315692 86.4702384
S 19.2379633 41.4329068 83.4498397
C 20.4513290 40.9071150 84.6025371
C 19.2139267 39.8633124 82.6852826
C 20.7913095 39.5782930 84.3709831
C 20.1021994 38.9935283 83.3016247
H 21.5250027 39.0574350 84.9727855
H 20.2449165 37.9687937 82.9852308
S 18.6109755 38.3455761 80.4189344
C 18.4701380 39.7438911 81.4550959
C 17.5286776 39.0538128 79.2338359
C 17.6279332 40.6957514 80.8982471
C 17.1039831 40.3083420 79.6589961
H 17.4016988 41.6397519 81.3759868
H 16.4286554 40.9179682 79.0724356
N 15.2367960 40.9679090 76.8556936
C 15.6011496 39.8748479 76.9541676
N 14.4666724 36.6838087 76.2123351
C 15.1748066 37.5118217 76.5985657
C 16.0545982 38.5292386 77.0776911
C 17.2620658 38.1820877 77.6202088
H 17.4738301 37.1174875 77.6491110
test/test_particle_collection/test.py 0 → 100755
View file @ a6b5ebbd
#! /usr/bin/env python
import soap
import soap.tools
import os
import numpy as np
from momo import osio, endl, flush
options = soap.Options()
options.excludeCenters(['H'])
options.excludeTargets([])
options.set('radialbasis.type', 'gaussian')
options.set('radialbasis.mode', 'adaptive')
options.set('radialbasis.N', 9)
options.set('radialbasis.sigma', 0.9)
options.set('radialbasis.integration_steps', 15)
options.set('radialcutoff.Rc', 6.8)
options.set('radialcutoff.Rc_width', 0.5)
options.set('radialcutoff.type', 'shifted-cosine')
options.set('radialcutoff.center_weight', 1.)
options.set('angularbasis.type', 'spherical-harmonic')
options.set('angularbasis.L', 6)
options.set('densitygrid.N', 20)
options.set('densitygrid.dx', 0.15)
# STRUCTURE
xyzfile = 'pair_1_185.xyz'
top = [('dcv2t', 2, 28)]
config = soap.tools.ase_load_single(xyzfile)
structure = soap.tools.setup_structure_ase(config.config_file, config.atoms, top)
seg1 = structure.getSegment(1)
seg2 = structure.getSegment(2)
print seg1.id, seg1.name
print seg2.id, seg2.name
for particle in seg1.particles:
particle.weight = -1.
for particle in seg2.particles:
particle.weight = +1.
for part in structure.particles:
print part.id, part.weight
# SPECTRUM
spectrum = soap.Spectrum(structure, options)
spectrum.compute(seg1, seg2)
spectrum.compute(seg2, seg1)
spectrum.computePower()
spectrum.save('test.arch')
atomic = spectrum.getAtomic(1, "S")
atomic.getLinear("").writeDensityOnGrid('test.cube', options, structure, atomic.getCenter(), True)
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