diff --git a/.gitignore b/.gitignore
index ed5c116b2ba209e5a6d3b75d856c21e99addd32d..08b2c79005fd7443817b6cd54013981bbbae9e36 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,3 +9,4 @@ test/back
*.coeff
*.cube
*.arch
+*.back
diff --git a/src/spectrum.cpp b/src/spectrum.cpp
index 15ffd5055ec64ae8546b2bd687c6f82c42d38905..e11b41483285191f85e2be7e1a27549281633e10 100644
--- a/src/spectrum.cpp
+++ b/src/spectrum.cpp
@@ -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 ¢ers, 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>())
diff --git a/src/spectrum.hpp b/src/spectrum.hpp
index 0e2623cacf08944f2a183200b6d73702e0a03563..ab15986323326f573e1dff40980e1794ffd2e1ea 100644
--- a/src/spectrum.hpp
+++ b/src/spectrum.hpp
@@ -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);
diff --git a/src/structure.cpp b/src/structure.cpp
index d9f25c983adae643ebb04e301145f33a1cf52422..b0c9e6a74ff7ff53e2017f5eadce651927ddf23b 100644
--- a/src/structure.cpp
+++ b/src/structure.cpp
@@ -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 */
diff --git a/src/structure.hpp b/src/structure.hpp
index 2834d2caa6bce587af66fafa26d0f7497f41770c..a95b5be8245602a04ffac2f044c6651ab697a647 100644
--- a/src/structure.hpp
+++ b/src/structure.hpp
@@ -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; }
diff --git a/src/tools/loadwrite.py b/src/tools/loadwrite.py
index 1ea6a5722377b68988f552c53f24192746442fa8..438b0fe02fc79e3b7efaaae097baa379ff1d6fd9 100644
--- a/src/tools/loadwrite.py
+++ b/src/tools/loadwrite.py
@@ -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,
diff --git a/test/test_particle_collection/pair_1_185.xyz b/test/test_particle_collection/pair_1_185.xyz
new file mode 100644
index 0000000000000000000000000000000000000000..8b91ca28857ee7b2a8d1e779a355bd157c6ba8a4
--- /dev/null
+++ b/test/test_particle_collection/pair_1_185.xyz
@@ -0,0 +1,58 @@
+ 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
diff --git a/test/test_particle_collection/test.py b/test/test_particle_collection/test.py
new file mode 100755
index 0000000000000000000000000000000000000000..5cf3488be2bc3f7c1a467210378696c38d002434
--- /dev/null
+++ b/test/test_particle_collection/test.py
@@ -0,0 +1,56 @@
+#! /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)
+