diff --git a/src/angularbasis.cpp b/src/angularbasis.cpp
index 7449a00c1619d01fcd7268c8c792b364114ab472..89da443f333a5231b72e7ef266643d2ae027d5e1 100644
--- a/src/angularbasis.cpp
+++ b/src/angularbasis.cpp
@@ -1,4 +1,5 @@
#include "angularbasis.hpp"
+
#include <math.h>
#include <boost/math/special_functions/spherical_harmonic.hpp>
@@ -9,29 +10,6 @@ void AngularBasis::configure(Options &options) {
_L = options.get<int>("angularbasis.L");
}
-AngularCoefficients AngularBasis::computeCoefficients(vec d, double r) {
- AngularCoefficients coeffs(this->_L);
- if (r < 1e-10) {
- // Radius is virtually zero => set isotropic
- std::complex<double> c00 = boost::math::spherical_harmonic<double,double>(0, 0, 0.0, 0.0);
- coeffs.set(0, 0, c00);
- }
- else {
- // Radius not zero => compute moments
- double theta = acos(d.getZ());
- double phi = atan2(d.getY(), d.getX());
- // Shift [-pi, -0] to [pi, 2*pi]
- if (phi < 0.) phi += 2*M_PI;
- for (int l=0; l<=_L; ++l) {
- for (int m=-l; m<=l; ++m) {
- std::complex<double> clm = boost::math::spherical_harmonic<double,double>(l, m, theta, phi);
- coeffs.set(l, m, clm);
- }
- }
- }
- return coeffs;
-}
-
void AngularBasis::computeCoefficients(vec d, double r, angcoeff_t &save_here) {
if (r < AngularBasis::RADZERO) {
std::complex<double> c00 = boost::math::spherical_harmonic<double,double>(0, 0, 0.0, 0.0);
@@ -53,20 +31,9 @@ void AngularBasis::computeCoefficients(vec d, double r, angcoeff_t &save_here) {
return;
}
-AngularCoefficients AngularBasis::computeCoefficientsAllZero() {
- AngularCoefficients coeffs(this->_L);
- for (int l=0; l<=_L; ++l) {
- for (int m=-l; m<=l; ++m) {
- coeffs.set(l, m, std::complex<double>(0,0));
- }
- }
- return coeffs;
-}
-
void AngularBasisFactory::registerAll(void) {
AngularBasisOutlet().Register<AngularBasis>("spherical-harmonic");
}
-
-}
+} /* CLOSE NAMESPACE */
diff --git a/src/angularbasis.hpp b/src/angularbasis.hpp
index 8ede17028ab5493e53b0a20422bdd6fbc03e0611..d5affe2fc8fca1be86e055df7ad82938d87f0d68 100644
--- a/src/angularbasis.hpp
+++ b/src/angularbasis.hpp
@@ -13,34 +13,6 @@ namespace soap {
namespace ub = boost::numeric::ublas;
-class AngularCoefficients : public ub::vector< std::complex<double> >
-{
-public:
-
- AngularCoefficients(int L) : _L(L) {
- this->resize((L+1)*(L+1));
- for (int i=0; i!=size(); ++i) {
- (*this)[i] = std::complex<double>(0,0);
- }
- }
- void set(int l, int m, std::complex<double> c) {
- if (this->checkSize(l, m)) (*this)[l*l+l+m] = c;
- else throw soap::base::OutOfRange("AngularCoefficients::set");
- }
- std::complex<double> &get(int l, int m) {
- if (this->checkSize(l, m)) return (*this)[l*l+l+m];
- else throw soap::base::OutOfRange("AngularCoefficients::get");
- }
- bool checkSize(int l, int m) { return (std::abs(m) <= l && l <= _L); }
- void conjugate() {
- for (int i = 0; i != size(); ++i) {
- (*this)[i] = std::conj((*this)[i]);
- }
- }
-protected:
- int _L;
-};
-
class AngularBasis
{
@@ -48,14 +20,13 @@ public:
typedef ub::vector< std::complex<double> > angcoeff_t;
typedef ub::zero_vector< std::complex<double> > angcoeff_zero_t;
+ AngularBasis() : _type("spherical-harmonic"), _L(0) {;}
+ virtual ~AngularBasis() {;}
+
std::string &identify() { return _type; }
const int &L() { return _L; }
- AngularBasis() : _type("spherical-harmonic"), _L(0) {;}
- virtual ~AngularBasis() {;}
virtual void configure(Options &options);
virtual void computeCoefficients(vec d, double r, angcoeff_t &save_here);
- virtual AngularCoefficients computeCoefficients(vec d, double r);
- virtual AngularCoefficients computeCoefficientsAllZero();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
@@ -64,15 +35,13 @@ public:
}
protected:
-
std::string _type;
int _L;
static constexpr double RADZERO = 1e-10;
};
-class AngularBasisFactory
- : public soap::base::ObjectFactory<std::string, AngularBasis>
+class AngularBasisFactory : public soap::base::ObjectFactory<std::string, AngularBasis>
{
private:
AngularBasisFactory() {}
@@ -98,6 +67,6 @@ inline AngularBasis *AngularBasisFactory::create(const std::string &key) {
}
}
-}
+} /* CLOSE NAMESPACE */
#endif
diff --git a/src/atomicspectrum.cpp b/src/atomicspectrum.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5c039531a7ed2180a7adb04b04097da65ae44e5f
--- /dev/null
+++ b/src/atomicspectrum.cpp
@@ -0,0 +1,11 @@
+#include <fstream>
+#include <boost/format.hpp>
+#include <boost/archive/binary_oarchive.hpp>
+#include <boost/archive/binary_iarchive.hpp>
+
+#include "atomicspectrum.hpp"
+
+namespace soap {
+
+}
+
diff --git a/src/atomicspectrum.hpp b/src/atomicspectrum.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..934f67e2f7af3efb5ce4e74c05093bebcac564a2
--- /dev/null
+++ b/src/atomicspectrum.hpp
@@ -0,0 +1,134 @@
+#ifndef _SOAP_ATOMICSPECTRUM_HPP
+#define _SOAP_ATOMICSPECTRUM_HPP
+
+#include <assert.h>
+#include <boost/serialization/vector.hpp>
+#include <boost/serialization/map.hpp>
+#include <boost/serialization/complex.hpp>
+#include <boost/serialization/base_object.hpp>
+
+#include "base/logger.hpp"
+#include "types.hpp"
+#include "globals.hpp"
+#include "options.hpp"
+#include "structure.hpp"
+#include "basis.hpp"
+#include "power.hpp"
+
+
+namespace soap {
+
+
+class AtomicSpectrum : public std::map<std::string, BasisExpansion*>
+{
+public:
+ typedef BasisExpansion qnlm_t;
+ 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')
+
+ 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;
+ }
+ 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;
+ }
+ }
+
+ template<class Archive>
+ void serialize(Archive &arch, const unsigned int version) {
+ // Center & basis
+ arch & _center;
+ arch & _center_pos;
+ arch & _center_type;
+ arch & _basis;
+ // Qnlm's
+ arch & _map_qnlm;
+ arch & _qnlm_generic;
+ // Xnkl's
+ arch & _map_xnkl;
+ arch & _xnkl_generic_coherent;
+ arch & _xnkl_generic_incoherent;
+ return;
+ }
+protected:
+ // CENTER & BASIS LINKS
+ Particle *_center;
+ vec _center_pos;
+ std::string _center_type;
+ Basis *_basis;
+ // DENSITY EXPANSION
+ map_qnlm_t _map_qnlm;
+ qnlm_t *_qnlm_generic;
+ // POWER DENSITY EXPANSION
+ map_xnkl_t _map_xnkl;
+ xnkl_t *_xnkl_generic_coherent;
+ xnkl_t *_xnkl_generic_incoherent;
+};
+
+
+}
+
+#endif /* _SOAP_ATOMICSPECTRUM_HPP_ */
diff --git a/src/basis.cpp b/src/basis.cpp
index 0415e53c6c2ee2fb25313669e12886d8349e25a1..5af0576b1a6d3816dbc127f7bed156aa15e2a821 100644
--- a/src/basis.cpp
+++ b/src/basis.cpp
@@ -3,5 +3,190 @@
namespace soap {
+// =====
+// Basis
+// =====
+
+Basis::Basis(Options *options) :
+ _options(options) {
+ GLOG() << "Configuring basis ..." << std::endl;
+ // CONFIGURE RADIAL BASIS
+ _radbasis = RadialBasisOutlet().create(_options->get<std::string>("radialbasis.type"));
+ _radbasis->configure(*options);
+ // 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::~Basis() {
+ delete _radbasis;
+ _radbasis = NULL;
+ delete _angbasis;
+ _angbasis = NULL;
+ delete _cutoff;
+ _cutoff = NULL;
+}
+
+// ==============
+// BasisExpansion
+// ==============
+
+BasisExpansion::BasisExpansion(Basis *basis) :
+ _basis(basis),
+ _radbasis(basis->getRadBasis()),
+ _angbasis(basis->getAngBasis()) {
+ int L = _angbasis->L();
+ int N = _radbasis->N();
+ _radcoeff = RadialBasis::radcoeff_zero_t(N,L+1);
+ _angcoeff = AngularBasis::angcoeff_zero_t((L+1)*(L+1));
+ _coeff = coeff_zero_t(N,(L+1)*(L+1));
+}
+
+BasisExpansion::~BasisExpansion() {
+ _basis = NULL;
+ _radbasis = NULL;
+ _angbasis = NULL;
+ _radcoeff.clear();
+ _angcoeff.clear();
+}
+
+void BasisExpansion::computeCoefficients(double r, vec d, double weight, double sigma) {
+ _radbasis->computeCoefficients(r, sigma, _radcoeff);
+ _angbasis->computeCoefficients(d, r, _angcoeff);
+ for (int n = 0; n != _radbasis->N(); ++n) {
+ for (int l = 0; l != _angbasis->L()+1; ++l) {
+ for (int m = -l; m != l+1; ++m) {
+ _coeff(n, l*l+l+m) = _radcoeff(n,l)*_angcoeff(l*l+l+m);
+ }
+ }
+ }
+ _coeff *= weight;
+ return;
+}
+
+void BasisExpansion::conjugate() {
+ for (int n = 0; n != _radbasis->N(); ++n) {
+ for (int l = 0; l != _angbasis->L()+1; ++l) {
+ for (int m = -l; m != l+1; ++m) {
+ _coeff(n, l*l+l+m) = std::conj(_coeff(n, l*l+l+m));
+ }
+ }
+ }
+}
+
+void BasisExpansion::writeDensityOnGrid(
+ std::string filename,
+ Options *options,
+ Structure *structure,
+ Particle *center,
+ bool fromExpansion) {
+
+ if (_angbasis == NULL || _radbasis == NULL) {
+ throw soap::base::APIError("<BasisExpansion::writeDensityOnGrid> "
+ "Object not linked against basis.");
+ }
+
+ int I = options->get<int>("densitygrid.N");
+ int Nx = 2*I+1;
+ int Ny = Nx;
+ int Nz = Nx;
+
+ double dx = options->get<double>("densitygrid.dx");
+ double dy = dx;
+ double dz = dx;
+
+ double conv = soap::constants::ANGSTROM_TO_BOHR;
+
+ vec r0 = -I * vec(dx,dy,dz);
+
+ std::ofstream ofs;
+ ofs.open(filename.c_str(), std::ofstream::out);
+ if (!ofs.is_open()) {
+ throw soap::base::IOError("Bad file handle: " + filename);
+ }
+
+ ofs << "DENSITY ON GRID" << std::endl;
+ ofs << "Generated by BasisCoefficients::writeDensityOnGrid" << std::endl;
+ ofs << boost::format("%1$d %2$+1.4f %3$+1.4f %4$+1.4f")
+ % structure->particles().size() % (r0.x()*conv) % (r0.y()*conv) % (r0.z()*conv) << std::endl;
+ ofs << boost::format("%1$d %2$+1.4f +0.0000 +0.0000") % Nx % (dx*conv) << std::endl;
+ ofs << boost::format("%1$d +0.0000 %2$+1.4f +0.0000") % Ny % (dy*conv) << std::endl;
+ ofs << boost::format("%1$d +0.0000 +0.0000 %2$+1.4f") % Nz % (dz*conv) << std::endl;
+
+ Structure::particle_it_t pit;
+ for (pit = structure->beginParticles(); pit != structure->endParticles(); ++pit) {
+ vec dr = structure->connect(center->getPos(), (*pit)->getPos());
+ ofs << boost::format("%1$d 0.0 %2$+1.4f %3$+1.4f %4$+1.4f\n")
+ % (*pit)->getTypeId() % (dr.x()*conv) % (dr.y()*conv) % (dr.z()*conv);
+ }
+
+ GLOG() << "Fill grid " << std::flush;
+ double int_density_dr = 0.0;
+ int ijk = 0;
+ for (int i = -I; i <= I; ++i) {
+ if (((i+I) % 5) == 0) GLOG() << "." << std::flush;
+ for (int j = -I; j <= I; ++j) {
+ for (int k = -I; k <= I; ++k) {
+ double density_dr = 0.0;
+
+ // Note that expansion is computed with respect to center-particle position:
+ // Hence no shifting required
+ vec dr(i*dx, j*dy, k*dz);
+ double r = soap::linalg::abs(dr);
+ vec d = dr/r;
+
+ // DENSITY BASED ON EXPANSION
+ if (fromExpansion) {
+ BasisExpansion density_exp_dr(_basis);
+ density_exp_dr.computeCoefficients(r, d, 1., 0.);
+ density_exp_dr.conjugate();
+
+ BasisExpansion::coeff_t &c_nlm_dr = density_exp_dr.getCoefficients();
+ BasisExpansion::coeff_t &c_nlm = this->getCoefficients();
+ // Something like this would be desirable to replace n-l-m loop below:
+ //double density_dr = ub::inner_prod(c_nlm, c_nlm_dr);
+ for (int n = 0; n < _radbasis->N(); ++n) {
+ for (int l = 0; l <= _angbasis->L(); ++l) {
+ for (int m = -l; m <= l; ++m) {
+ double voxel_density_dr = (c_nlm(n, l*l+l+m)*c_nlm_dr(n, l*l+l+m)).real();
+ density_dr += voxel_density_dr;
+ int_density_dr += voxel_density_dr*dx*dy*dz;
+ }
+ }
+ }
+ }
+
+ // DENSITY BASED ON SMEARED PARTICLES
+ else {
+ for (pit = structure->beginParticles(); pit != structure->endParticles(); ++pit) {
+ vec dr_center_particle = structure->connect(center->getPos(), (*pit)->getPos());
+ vec dr_particle_target = dr - dr_center_particle;
+ double r_particle_target = soap::linalg::abs(dr_particle_target);
+ double sigma = (*pit)->getSigma();
+ double weight = (*pit)->getWeight();
+ density_dr += weight*pow(1./(2.*M_PI*sigma*sigma), 1.5)
+ * exp(-r_particle_target*r_particle_target/(2*sigma*sigma));
+ }
+ int_density_dr += density_dr*dx*dy*dz;
+ }
+
+ ofs << density_dr << " ";
+ ijk += 1;
+ ijk = ijk % 6;
+ if (ijk == 0) {
+ ofs << std::endl;
+ }
+
+ }
+ }
+ }
+ ofs.close();
+ GLOG() << " Volume integral = " << int_density_dr << std::endl;
+ return;
+}
+
}
diff --git a/src/basis.hpp b/src/basis.hpp
index 4fecaac6647660e10035b3f3c0f54b6f0489f446..a08e74f698c98187ede21f06427a8349ea70077f 100644
--- a/src/basis.hpp
+++ b/src/basis.hpp
@@ -21,27 +21,10 @@ namespace ub = boost::numeric::ublas;
class Basis
{
public:
- Basis(Options *options) : _options(options) {
- GLOG() << "Configuring basis ..." << std::endl;
- // CONFIGURE RADIAL BASIS
- _radbasis = RadialBasisOutlet().create(_options->get<std::string>("radialbasis.type"));
- _radbasis->configure(*options);
- // 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(Options *options);
Basis() : _options(NULL), _radbasis(NULL), _angbasis(NULL), _cutoff(NULL) {;}
- ~Basis() {
- delete _radbasis;
- _radbasis = NULL;
- delete _angbasis;
- _angbasis = NULL;
- delete _cutoff;
- _cutoff = NULL;
- }
+ ~Basis();
+
RadialBasis *getRadBasis() { return _radbasis; }
AngularBasis *getAngBasis() { return _angbasis; }
CutoffFunction *getCutoff() { return _cutoff; }
@@ -68,162 +51,20 @@ public:
typedef ub::matrix< std::complex<double> > coeff_t;
typedef ub::zero_matrix< std::complex<double> > coeff_zero_t;
- BasisExpansion(Basis *basis) :
- _basis(basis), _radbasis(basis->getRadBasis()), _angbasis(basis->getAngBasis()) {
- int L = _angbasis->L();
- int N = _radbasis->N();
- _radcoeff = RadialBasis::radcoeff_zero_t(N,L+1);
- _angcoeff = AngularBasis::angcoeff_zero_t((L+1)*(L+1));
- _coeff = coeff_zero_t(N,(L+1)*(L+1));
- }
- BasisExpansion() :
- _basis(NULL), _radbasis(NULL), _angbasis(NULL) {;}
- ~BasisExpansion() {
- _basis = NULL;
- _radbasis = NULL;
- _angbasis = NULL;
- _radcoeff.clear();
- _angcoeff.clear();
- }
+ BasisExpansion(Basis *basis);
+ BasisExpansion() : _basis(NULL), _radbasis(NULL), _angbasis(NULL) {;}
+ ~BasisExpansion();
+
Basis *getBasis() { return _basis; }
coeff_t &getCoefficients() { return _coeff; }
- void computeCoefficients(double r, vec d, double weight, double sigma) {
- _radbasis->computeCoefficients(r, sigma, _radcoeff);
- _angbasis->computeCoefficients(d, r, _angcoeff);
- for (int n = 0; n != _radbasis->N(); ++n) {
- for (int l = 0; l != _angbasis->L()+1; ++l) {
- for (int m = -l; m != l+1; ++m) {
- _coeff(n, l*l+l+m) = _radcoeff(n,l)*_angcoeff(l*l+l+m);
- }
- }
- }
- _coeff *= weight;
- return;
- }
- void add(BasisExpansion &other) {
- _coeff = _coeff + other._coeff;
- }
- void conjugate() {
- for (int n = 0; n != _radbasis->N(); ++n) {
- for (int l = 0; l != _angbasis->L()+1; ++l) {
- for (int m = -l; m != l+1; ++m) {
- _coeff(n, l*l+l+m) = std::conj(_coeff(n, l*l+l+m));
- }
- }
- }
- }
- void writeDensityOnGrid(
- std::string filename,
- Options *options,
- Structure *structure,
- Particle *center,
- bool fromExpansion) {
-
- if (_angbasis == NULL || _radbasis == NULL) {
- throw soap::base::APIError("<BasisExpansion::writeDensityOnGrid> "
- "Object not linked against basis.");
- }
-
- int I = options->get<int>("densitygrid.N");
- int Nx = 2*I+1;
- int Ny = Nx;
- int Nz = Nx;
-
- double dx = options->get<double>("densitygrid.dx");
- double dy = dx;
- double dz = dx;
-
- double conv = soap::constants::ANGSTROM_TO_BOHR;
-
- vec r0 = -I * vec(dx,dy,dz);
-
- std::ofstream ofs;
- ofs.open(filename.c_str(), std::ofstream::out);
- if (!ofs.is_open()) {
- throw soap::base::IOError("Bad file handle: " + filename);
- }
-
- ofs << "DENSITY ON GRID" << std::endl;
- ofs << "Generated by BasisCoefficients::writeDensityOnGrid" << std::endl;
- ofs << boost::format("%1$d %2$+1.4f %3$+1.4f %4$+1.4f")
- % structure->particles().size() % (r0.x()*conv) % (r0.y()*conv) % (r0.z()*conv) << std::endl;
- ofs << boost::format("%1$d %2$+1.4f +0.0000 +0.0000") % Nx % (dx*conv) << std::endl;
- ofs << boost::format("%1$d +0.0000 %2$+1.4f +0.0000") % Ny % (dy*conv) << std::endl;
- ofs << boost::format("%1$d +0.0000 +0.0000 %2$+1.4f") % Nz % (dz*conv) << std::endl;
-
- Structure::particle_it_t pit;
- for (pit = structure->beginParticles(); pit != structure->endParticles(); ++pit) {
- vec dr = structure->connect(center->getPos(), (*pit)->getPos());
- ofs << boost::format("%1$d 0.0 %2$+1.4f %3$+1.4f %4$+1.4f\n")
- % (*pit)->getTypeId() % (dr.x()*conv) % (dr.y()*conv) % (dr.z()*conv);
- }
-
- GLOG() << "Fill grid " << std::flush;
- double int_density_dr = 0.0;
- int ijk = 0;
- for (int i = -I; i <= I; ++i) {
- if (((i+I) % 5) == 0) GLOG() << "." << std::flush;
- for (int j = -I; j <= I; ++j) {
- for (int k = -I; k <= I; ++k) {
- double density_dr = 0.0;
-
- // Note that expansion is computed with respect to center-particle position:
- // Hence no shifting required
- vec dr(i*dx, j*dy, k*dz);
- double r = soap::linalg::abs(dr);
- vec d = dr/r;
-
- // DENSITY BASED ON EXPANSION
- if (fromExpansion) {
- BasisExpansion density_exp_dr(_basis);
- density_exp_dr.computeCoefficients(r, d, 1., 0.);
- density_exp_dr.conjugate();
-
- BasisExpansion::coeff_t &c_nlm_dr = density_exp_dr.getCoefficients();
- BasisExpansion::coeff_t &c_nlm = this->getCoefficients();
- // Something like this would be desirable to replace n-l-m loop below:
- //double density_dr = ub::inner_prod(c_nlm, c_nlm_dr);
- for (int n = 0; n < _radbasis->N(); ++n) {
- for (int l = 0; l <= _angbasis->L(); ++l) {
- for (int m = -l; m <= l; ++m) {
- double voxel_density_dr = (c_nlm(n, l*l+l+m)*c_nlm_dr(n, l*l+l+m)).real();
- density_dr += voxel_density_dr;
- int_density_dr += voxel_density_dr*dx*dy*dz;
- }
- }
- }
- }
-
- // DENSITY BASED ON SMEARED PARTICLES
- else {
- for (pit = structure->beginParticles(); pit != structure->endParticles(); ++pit) {
- vec dr_center_particle = structure->connect(center->getPos(), (*pit)->getPos());
- vec dr_particle_target = dr - dr_center_particle;
- double r_particle_target = soap::linalg::abs(dr_particle_target);
- double sigma = (*pit)->getSigma();
- double weight = (*pit)->getWeight();
- density_dr += weight*pow(1./(2.*M_PI*sigma*sigma), 1.5)
- * exp(-r_particle_target*r_particle_target/(2*sigma*sigma));
- }
- int_density_dr += density_dr*dx*dy*dz;
- }
-
- ofs << density_dr << " ";
- ijk += 1;
- ijk = ijk % 6;
- if (ijk == 0) {
- ofs << std::endl;
- }
-
- }
- }
- }
- ofs.close();
- GLOG() << " Volume integral = " << int_density_dr << std::endl;
- return;
- }
RadialBasis::radcoeff_t &getRadCoeffs() { return _radcoeff; }
- AngularBasis::angcoeff_t &getAngCoeffs() { return _angcoeff; }
+ AngularBasis::angcoeff_t &getAngCoeffs() { return _angcoeff; }
+
+ void computeCoefficients(double r, vec d, double weight, double sigma);
+ void add(BasisExpansion &other) { _coeff = _coeff + other._coeff; }
+ void conjugate();
+ void writeDensityOnGrid(std::string filename, Options *options,
+ Structure *structure, Particle *center, bool fromExpansion);
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
@@ -243,149 +84,6 @@ private:
coeff_t _coeff; // access via (n, l*l+l+m)
};
-
-
-struct BasisCoefficients
-{
- BasisCoefficients(RadialCoefficients &c_n, AngularCoefficients &c_lm)
- : _c_n(c_n), _c_lm(c_lm), _radbasis(NULL), _angbasis(NULL) {
- _c_nlm.resize(c_n.size(), c_lm.size());
- _c_nlm = ub::outer_prod(_c_n, _c_lm);
- }
- void linkBasis(RadialBasis *radbasis, AngularBasis *angbasis) {
- _radbasis = radbasis;
- _angbasis = angbasis;
- }
- std::complex<double> &get(int n, int l, int m) {
- if (_c_lm.checkSize(l, m) && _c_n.checkSize(n)) {
- return _c_nlm(n, l*l+l+m);
- }
- else {
- throw soap::base::OutOfRange("BasisCoefficients::get");
- }
- }
- void add(BasisCoefficients &other) {
- _c_nlm = _c_nlm + other._c_nlm;
- }
-
- void writeDensityOnGrid(
- std::string filename,
- Options *options,
- Structure *structure,
- Particle *center,
- bool fromExpansion) {
-
- if (_angbasis == NULL || _radbasis == NULL) {
- throw soap::base::APIError("<BasisCoefficients::writeDensityOnGrid> "
- "Object not linked against basis.");
- }
-
- int I = options->get<int>("densitygrid.N");
- int Nx = 2*I+1;
- int Ny = Nx;
- int Nz = Nx;
-
- double dx = options->get<double>("densitygrid.dx");
- double dy = dx;
- double dz = dx;
-
- double conv = soap::constants::ANGSTROM_TO_BOHR;
-
- vec r0 = -I * vec(dx,dy,dz); // center->getPos();
-
- std::ofstream ofs;
- ofs.open(filename.c_str(), std::ofstream::out);
- if (!ofs.is_open()) {
- throw soap::base::IOError("Bad file handle: " + filename);
- }
-
- ofs << "DENSITY ON GRID" << std::endl;
- ofs << "Generated by BasisCoefficients::writeDensityOnGrid" << std::endl;
- ofs << boost::format("%1$d %2$+1.4f %3$+1.4f %4$+1.4f")
- % structure->particles().size() % (r0.x()*conv) % (r0.y()*conv) % (r0.z()*conv) << std::endl;
- ofs << boost::format("%1$d %2$+1.4f +0.0000 +0.0000") % Nx % (dx*conv) << std::endl;
- ofs << boost::format("%1$d +0.0000 %2$+1.4f +0.0000") % Ny % (dy*conv) << std::endl;
- ofs << boost::format("%1$d +0.0000 +0.0000 %2$+1.4f") % Nz % (dz*conv) << std::endl;
-
- Structure::particle_it_t pit;
- for (pit = structure->beginParticles(); pit != structure->endParticles(); ++pit) {
- vec dr = structure->connect(center->getPos(), (*pit)->getPos());
- ofs << boost::format("%1$d 0.0 %2$+1.4f %3$+1.4f %4$+1.4f\n")
- % (*pit)->getTypeId() % (dr.x()*conv) % (dr.y()*conv) % (dr.z()*conv);
- }
-
- GLOG() << "Fill grid " << std::flush;
- double int_density_dr = 0.0;
- int ijk = 0;
- for (int i = -I; i <= I; ++i) {
- if (((i+I) % 5) == 0) GLOG() << "." << std::flush;
- for (int j = -I; j <= I; ++j) {
- for (int k = -I; k <= I; ++k) {
- double density_dr = 0.0;
-
- // Note that expansion is computed with respect to center-particle position:
- // Hence no shifting required
- vec dr(i*dx, j*dy, k*dz);
- double r = soap::linalg::abs(dr);
- vec d = dr/r;
-
- // DENSITY BASED ON EXPANSION
- if (fromExpansion) {
- RadialCoefficients c_n_dr = _radbasis->computeCoefficients(r);
- AngularCoefficients c_lm_dr = _angbasis->computeCoefficients(d, r);
- c_lm_dr.conjugate();
- BasisCoefficients c_nlm_dr(c_n_dr, c_lm_dr);
- // Something like this would be desirable to replace n-l-m loop below:
- //double density_dr = ub::inner_prod(c_nlm, c_nlm_dr);
- for (int n = 0; n < _radbasis->N(); ++n) {
- for (int l = 0; l <= _angbasis->L(); ++l) {
- for (int m = -l; m <= l; ++m) {
- double voxel_density_dr = (this->get(n, l, m)*c_nlm_dr.get(n, l, m)).real();
- density_dr += voxel_density_dr;
- int_density_dr += voxel_density_dr*dx*dy*dz;
- }
- }
- }
- }
-
- // DENSITY BASED ON SMEARED PARTICLES
- else {
- for (pit = structure->beginParticles(); pit != structure->endParticles(); ++pit) {
- vec dr_center_particle = structure->connect(center->getPos(), (*pit)->getPos());
- vec dr_particle_target = dr - dr_center_particle;
- double r_particle_target = soap::linalg::abs(dr_particle_target);
- double sigma = 0.5;
- density_dr += pow(1./(2.*M_PI*sigma*sigma), 1.5)
- * exp(-r_particle_target*r_particle_target/(2*sigma*sigma));
- }
- int_density_dr += density_dr*dx*dy*dz;
- }
-
- ofs << density_dr << " ";
- ijk += 1;
- ijk = ijk % 6;
- if (ijk == 0) {
- ofs << std::endl;
- }
-
- }
- }
- }
-
- GLOG() << " Volume integral = " << int_density_dr << std::endl;
-
- ofs.close();
-
- return;
- }
-
- ub::matrix< std::complex<double> > _c_nlm;
- RadialCoefficients _c_n;
- AngularCoefficients _c_lm;
- RadialBasis *_radbasis;
- AngularBasis *_angbasis;
-};
-
}
#endif
diff --git a/src/cutoff.cpp b/src/cutoff.cpp
index e0529bf772c27be06ff8a0600e5dd4d2b8260da2..df61c1dd3593747fbe9af16343c882f71d0fbb83 100644
--- a/src/cutoff.cpp
+++ b/src/cutoff.cpp
@@ -7,6 +7,11 @@ void CutoffFunction::configure(Options &options) {
_Rc = options.get<double>("radialcutoff.Rc");
_Rc_width = options.get<double>("radialcutoff.Rc_width");
_center_weight = options.get<double>("radialcutoff.center_weight");
+
+ GLOG() << "Weighting function with "
+ << "Rc = " << _Rc
+ << ", _Rc_width = " << _Rc_width
+ << ", central weigth = " << _center_weight << std::endl;
}
void CutoffFunctionFactory::registerAll(void) {
diff --git a/src/cutoff.hpp b/src/cutoff.hpp
index b1dc8d70a3e6cda29fb7fb3a97495257134ad55d..439b8ebbb549530be775106742608a37e9d54d2d 100644
--- a/src/cutoff.hpp
+++ b/src/cutoff.hpp
@@ -7,6 +7,7 @@
#include "base/exceptions.hpp"
#include "base/objectfactory.hpp"
+#include "globals.hpp"
#include "options.hpp"
namespace soap {
diff --git a/src/functions.cpp b/src/functions.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f0586bd43917a31e82c56f358ec2eb5c31073ed8
--- /dev/null
+++ b/src/functions.cpp
@@ -0,0 +1,95 @@
+#include "functions.hpp"
+
+#include <boost/math/special_functions/erf.hpp>
+
+namespace soap {
+
+// ======================
+// RadialGaussian
+// ======================
+
+RadialGaussian::RadialGaussian(double r0, double sigma)
+: _r0(r0),
+ _sigma(sigma),
+ _alpha(1./(2.*sigma*sigma)) {
+
+ // COMPUTE NORMALIZATION S g^2 r^2 dr
+ // This normalization is to be used for radial basis functions
+ double w = 2*_alpha;
+ double W0 = 2*_alpha*_r0;
+ _integral_r2_g2_dr =
+ 1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
+ 2*sqrt(w)*W0 +
+ sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
+ 1 - boost::math::erf<double>(-W0/sqrt(w))
+ )
+ );
+ _norm_r2_g2_dr = 1./sqrt(_integral_r2_g2_dr);
+
+ // COMPUTE NORMALIZATION S g r^2 dr
+ // This normalization is to be used for "standard" radial Gaussians
+ w = _alpha;
+ W0 = _alpha*_r0;
+ _integral_r2_g_dr =
+ 1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
+ 2*sqrt(w)*W0 +
+ sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
+ 1 - boost::math::erf<double>(-W0/sqrt(w))
+ )
+ );
+ _norm_r2_g_dr = 1./_integral_r2_g_dr;
+}
+
+RadialGaussian::RadialGaussian() :
+ _r0(0.),
+ _sigma(0.),
+ _alpha(1./0.),
+ _integral_r2_g2_dr(0.),
+ _norm_r2_g2_dr(0.),
+ _integral_r2_g_dr(0.),
+ _norm_r2_g_dr(0.) {
+ ;
+}
+
+double RadialGaussian::at(double r) {
+ double p = _alpha*(r-_r0)*(r-_r0);
+ if (p < 40) return _norm_r2_g2_dr * exp(-p);
+ else return 0.0;
+}
+
+// ======================
+// SphericalGaussian
+// ======================
+
+SphericalGaussian::SphericalGaussian(vec r0, double sigma) :
+ _r0(r0), _sigma(sigma), _alpha(1./(2*sigma*sigma)) {
+ _norm_g_dV = pow(_alpha/M_PI, 1.5);
+}
+
+// ======================
+// ModSphBessel1stKind
+// ======================
+
+std::vector<double> ModifiedSphericalBessel1stKind::eval(int degree, double r) {
+ std::vector<double> il;
+ if (r < RADZERO) {
+ il.push_back(1.);
+ il.push_back(0.);
+ }
+ else {
+ il.push_back(sinh(r)/r);
+ il.push_back(cosh(r)/r - sinh(r)/(r*r));
+ }
+ for (int l = 2; l <= degree; ++l) {
+ if (r < RADZERO) {
+ il.push_back(0.);
+ }
+ else {
+ if (il[l-1] < SPHZERO) il.push_back(0.);
+ else il.push_back( il[l-2] - (2*(l-1)+1)/r*il[l-1] );
+ }
+ }
+ return il;
+}
+
+} /* CLOSE NAMESPACE */
diff --git a/src/functions.hpp b/src/functions.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c36807599de6dcf2901fbc1bdb4cf0418672248f
--- /dev/null
+++ b/src/functions.hpp
@@ -0,0 +1,57 @@
+#ifndef _SOAP_FUNCTIONS_H
+#define _SOAP_FUNCTIONS_H
+
+#include <math.h>
+#include <vector>
+
+#include "types.hpp"
+
+namespace soap {
+
+struct RadialGaussian
+{
+ RadialGaussian(double r0, double sigma);
+ RadialGaussian();
+ double at(double r);
+
+ template<class Archive>
+ void serialize(Archive &arch, const unsigned int version) {
+ arch & _r0;
+ arch & _sigma;
+ arch & _alpha;
+ arch & _norm_r2_g2_dr;
+ arch & _norm_r2_g_dr;
+ }
+
+ double _r0;
+ double _sigma;
+ double _alpha;
+ // Integral S g^2 r^2 dr
+ double _integral_r2_g2_dr;
+ double _norm_r2_g2_dr;
+ // Integral S g r^2 dr
+ double _integral_r2_g_dr;
+ double _norm_r2_g_dr;
+};
+
+struct SphericalGaussian
+{
+ SphericalGaussian(vec r0, double sigma);
+
+ vec _r0;
+ double _sigma;
+ double _alpha;
+ double _norm_g_dV;
+};
+
+struct ModifiedSphericalBessel1stKind
+{
+ static std::vector<double> eval(int degree, double r);
+ static constexpr double RADZERO = 1e-10;
+ static constexpr double SPHZERO = 1e-4;
+};
+
+
+}
+
+#endif
diff --git a/src/options.cpp b/src/options.cpp
index 97ec249d331c316353254887d32b621355f47e0e..279c3d0d178595b2e663e555da5b6291546e85ea 100644
--- a/src/options.cpp
+++ b/src/options.cpp
@@ -2,7 +2,41 @@
namespace soap {
-
+
+Options::Options() :
+ _center_excludes(boost::python::list()) {
+
+ // Set defaults
+ this->set("radialbasis.type", "gaussian");
+ this->set("radialbasis.mode", "equispaced");
+ this->set("radialbasis.N", 9);
+ this->set("radialbasis.sigma", 0.5);
+ this->set("radialbasis.integration_steps", 15);
+ this->set("radialcutoff.type", "shifted-cosine");
+ this->set("radialcutoff.Rc", 4.);
+ this->set("radialcutoff.Rc_width", 0.5);
+ this->set("radialcutoff.center_weight", 1.);
+ this->set("angularbasis.type", "spherical-harmonic");
+ this->set("angularbasis.L", 6);
+ this->set("densitygrid.N", 20);
+ this->set("densitygrid.dx", 0.15);
+}
+
+//template<typename return_t>
+//return_t Options::get(std::string key) {
+ //return soap::lexical_cast<return_t, std::string>(_key_value_map[key], "wrong or missing type in " + key);
+//}
+
+std::string Options::summarizeOptions() {
+ std::string info = "";
+ info += "Options:\n";
+ 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;
+}
void Options::registerPython() {
using namespace boost::python;
@@ -11,11 +45,8 @@ void Options::registerPython() {
void (Options::*set_string)(std::string, std::string) = &Options::set;
class_<Options>("Options")
- .def("configureRealBasis", &Options::configureRealBasis)
- .def("configureReciprocalBasis", &Options::configureReciprocalBasis)
.def("summarizeOptions", &Options::summarizeOptions)
.def("configureCenters", &Options::configureCenters)
- .def("configureReciprocalLattice", &Options::configureReciprocalLattice)
.def("set", set_int)
.def("set", set_double)
.def("set", set_string);
diff --git a/src/options.hpp b/src/options.hpp
index e7cb8540d3259990a4b57d00184d8bd86b4eaadc..9f81594d809a2ab36a24755efa17346d1e9ba841 100644
--- a/src/options.hpp
+++ b/src/options.hpp
@@ -1,89 +1,37 @@
#ifndef _SOAP_OPTIONS_HPP
#define _SOAP_OPTIONS_HPP
+
+#include <map>
#include <boost/format.hpp>
+
#include "types.hpp"
-#include <map>
namespace soap {
class Options
{
public:
- Options()
- : _N_real(8), _L_real(6), _Rc_real(5.),
- _N_recip(8), _L_recip(6), _Rc_recip(5.),
- _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("radialbasis.integration_steps", 15);
- this->set("radialbasis.mode", "equispaced");
- }
- ~Options() {}
+ typedef std::map<std::string, std::string> map_options_t;
+ // CONSTRUCT & DESTRUCT
+ Options();
+ ~Options() {}
+
+ // GET & SET
template<typename return_t>
return_t get(std::string key) {
return soap::lexical_cast<return_t, std::string>(_key_value_map[key], "wrong or missing type in " + key);
}
-
void set(std::string key, std::string value) { _key_value_map[key] = value; }
void set(std::string key, int value) { this->set(key, boost::lexical_cast<std::string>(value)); }
void set(std::string key, double value) { this->set(key, boost::lexical_cast<std::string>(value)); }
+ void configureCenters(boost::python::list center_excludes) { _center_excludes = center_excludes; }
+ std::string summarizeOptions();
- /*
- template<class value_t>
- void set(std::string key, value_t value) {
- std::cout << "value" << value << "_" << std::endl;
- _key_value_map[key] = boost::lexical_cast<std::string>(value);
- }
- */
-
- void configureRealBasis(int N_real, int L_real, double Rc_real, std::string type) {
- _N_real = N_real;
- _L_real = L_real;
- _Rc_real = Rc_real;
- _type_real = type;
- _key_value_map["radialbasis.Rc"] = boost::lexical_cast<std::string>(_Rc_real);
- _key_value_map["radialbasis.L"] = boost::lexical_cast<std::string>(_L_real);
- _key_value_map["radialbasis.N"] = boost::lexical_cast<std::string>(_N_real);
- _key_value_map["radialbasis.type"] = boost::lexical_cast<std::string>(_type_real);
- }
- void configureReciprocalBasis(int N_recip, int L_recip, double Rc_recip) {
- _N_recip = N_recip;
- _L_recip = L_recip;
- _Rc_recip = Rc_recip;
- }
- void configureCenters(double center_w0, boost::python::list center_excludes) {
- _center_w0 = center_w0;
- _center_excludes = center_excludes;
- }
- void configureReciprocalLattice(vec b1, vec b2, vec b3) {
- _b1 = b1;
- _b2 = b2;
- _b3 = b3;
- }
- 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();
-
- 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;
- }
+ // PYTHON
static void registerPython();
+ // SERIALIZATION
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
arch & _key_value_map;
@@ -91,24 +39,8 @@ public:
}
private:
- // REAL SPACE
- int _N_real;
- int _L_real;
- double _Rc_real;
- std::string _type_real;
- // RECIPROCAL SPACE
- int _N_recip;
- int _L_recip;
- double _Rc_recip;
- vec _b1;
- vec _b2;
- vec _b3;
- // CENTERS
- double _center_w0;
boost::python::list _center_excludes;
-
- std::map<std::string, std::string> _key_value_map;
-
+ map_options_t _key_value_map;
};
}
diff --git a/src/power.cpp b/src/power.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9fdca82ede7501d1055b6a2a2944f2965c9ed0dc
--- /dev/null
+++ b/src/power.cpp
@@ -0,0 +1,6 @@
+#include "power.hpp"
+
+namespace soap {
+
+}
+
diff --git a/src/power.hpp b/src/power.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..de1dfc89fa1b318849164de7f45b40286708a6d2
--- /dev/null
+++ b/src/power.hpp
@@ -0,0 +1,42 @@
+#ifndef _SOAP_POWER_HPP
+#define _SOAP_POWER_HPP
+
+#include <string>
+#include <math.h>
+#include <vector>
+#include <fstream>
+
+#include "base/exceptions.hpp"
+#include "basis.hpp"
+
+namespace soap {
+
+namespace ub = boost::numeric::ublas;
+
+
+class PowerExpansion
+{
+public:
+ typedef ub::matrix< std::complex<double> > coeff_t;
+ typedef ub::zero_matrix< std::complex<double> > coeff_zero_t;
+
+ PowerExpansion() {
+ ;
+ }
+ PowerExpansion(BasisExpansion *basex) {
+ ;
+ }
+
+ template<class Archive>
+ void serialize(Archive &arch, const unsigned int version) {
+ ;
+ }
+
+private:
+
+ coeff_t _coeff; // access via (N*n+k, l) with shape (N*N, L+1)
+};
+
+}
+
+#endif
diff --git a/src/radialbasis.cpp b/src/radialbasis.cpp
index 92524099cb29265b3542568f11d2970854515d78..ae1b5372105bbbb819b1b0773c391434f5a72b61 100644
--- a/src/radialbasis.cpp
+++ b/src/radialbasis.cpp
@@ -23,33 +23,33 @@ namespace ub = boost::numeric::ublas;
void RadialBasis::configure(Options &options) {
_N = options.get<int>("radialbasis.N");
- _Rc = options.get<double>("radialcutoff.Rc");
+ _Rc = options.get<double>("radialcutoff.Rc"); // <- !! Modified below for adaptive basis sets !!
_integration_steps = options.get<int>("radialbasis.integration_steps");
_mode = options.get<std::string>("radialbasis.mode");
}
-RadialCoefficients RadialBasis::computeCoefficients(double r) {
- RadialCoefficients coeffs(this->_N);
- throw soap::base::NotImplemented("RadialBasis::computeCoefficients");
- return coeffs;
-}
-
void RadialBasis::computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here) {
throw soap::base::NotImplemented("RadialBasis::computeCoefficients");
return;
}
-RadialCoefficients RadialBasis::computeCoefficientsAllZero() {
- RadialCoefficients coeffs(this->_N);
- for (int i = 0; i < _N; ++i) coeffs(i) = 0.0;
- return coeffs;
-}
-
-
// ======================
// RadialBasisGaussian
// ======================
+RadialBasisGaussian::RadialBasisGaussian() {
+ _type = "gaussian";
+ _is_ortho = false;
+ _sigma = 0.5;
+}
+
+void RadialBasisGaussian::clear() {
+ for (basis_it_t bit=_basis.begin(); bit!=_basis.end(); ++bit) {
+ delete *bit;
+ }
+ _basis.clear();
+}
+
void RadialBasisGaussian::configure(Options &options) {
RadialBasis::configure(options);
_sigma = options.get<double>("radialbasis.sigma");
@@ -68,10 +68,10 @@ void RadialBasisGaussian::configure(Options &options) {
}
else if (_mode == "adaptive") {
//double delta = 0.5;
- int L = 6;
+ int L = options.get<int>("angularbasis.L");
double r = 0.;
double sigma = 0.;
- double sigma_0 = 0.5;
+ double sigma_0 = _sigma;
double sigma_stride_factor = 0.5;
// while (r < _Rc) {
// sigma = sqrt(4./(2*L+1))*(r+delta);
@@ -82,7 +82,7 @@ void RadialBasisGaussian::configure(Options &options) {
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;
+ //std::cout << r << " " << sigma << std::endl;
r = r + sigma_stride_factor*sigma;
}
double scale = 1.; //_Rc/(r-sigma);
@@ -91,9 +91,13 @@ void RadialBasisGaussian::configure(Options &options) {
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);
- std::cout << r << " " << sigma << std::endl;
+ //std::cout << r << " " << sigma << std::endl;
r = r+sigma_stride_factor*sigma;
}
+ _Rc = r - sigma_stride_factor*sigma;
+ options.set("radialcutoff.Rc", _Rc);
+ GLOG() << "Adjusted radial cutoff to " << _Rc
+ << " based on sigma_0 = " << sigma_0 << ", L = " << L << ", stride = " << sigma_stride_factor << std::endl;
}
else {
throw std::runtime_error("Not implemented.");
@@ -167,67 +171,6 @@ void RadialBasisGaussian::configure(Options &options) {
}
}
-RadialCoefficients RadialBasisGaussian::computeCoefficients(double r) {
-
- ub::vector<double> coeffs_raw(_N);
- basis_it_t it;
- int i;
- for (it = _basis.begin(), i = 0; it != _basis.end(); ++it, ++i) {
- coeffs_raw(i) = (*it)->at(r);
- }
-
- coeffs_raw = ub::prod(_Tij, coeffs_raw);
-
- RadialCoefficients coeffs(this->_N);
- for (int n = 0; n < _N; ++n) {
- coeffs.set(n, coeffs_raw(n));
- }
-
- return coeffs;
-}
-
-struct SphericalGaussian
-{
- SphericalGaussian(vec r0, double sigma) :
- _r0(r0), _sigma(sigma), _alpha(1./(2*sigma*sigma)) {
- _norm_g_dV = pow(_alpha/M_PI, 1.5);
- }
-
- vec _r0;
- double _sigma;
- double _alpha;
- double _norm_g_dV;
-};
-
-struct ModifiedSphericalBessel1stKind
-{
- static std::vector<double> eval(int degree, double r) {
- std::vector<double> il;
- if (r < RADZERO) {
- il.push_back(1.);
- il.push_back(0.);
- }
- else {
- il.push_back(sinh(r)/r);
- il.push_back(cosh(r)/r - sinh(r)/(r*r));
- }
- for (int l = 2; l <= degree; ++l) {
- if (r < RADZERO) {
- il.push_back(0.);
- }
- else {
- if (il[l-1] < SPHZERO) il.push_back(0.);
- else il.push_back( il[l-2] - (2*(l-1)+1)/r*il[l-1] );
- }
- }
- return il;
- }
-
- static constexpr double RADZERO = 1e-10;
- static constexpr double SPHZERO = 1e-4;
-};
-
-
void RadialBasisGaussian::computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here) {
// Delta-type expansion =>
// Second (l) dimension of <save_here> and <particle_sigma> ignored here
@@ -389,50 +332,6 @@ void RadialBasisGaussian::computeCoefficients(double r, double particle_sigma, r
return;
}
-
-// ======================
-// RadialGaussian
-// ======================
-
-RadialGaussian::RadialGaussian(double r0, double sigma)
-: _r0(r0),
- _sigma(sigma),
- _alpha(1./(2.*sigma*sigma)) {
-
- // COMPUTE NORMALIZATION S g^2 r^2 dr
- // This normalization is to be used for radial basis functions
- double w = 2*_alpha;
- double W0 = 2*_alpha*_r0;
- _integral_r2_g2_dr =
- 1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
- 2*sqrt(w)*W0 +
- sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
- 1 - boost::math::erf<double>(-W0/sqrt(w))
- )
- );
- _norm_r2_g2_dr = 1./sqrt(_integral_r2_g2_dr);
-
- // COMPUTE NORMALIZATION S g r^2 dr
- // This normalization is to be used for "standard" radial Gaussians
- w = _alpha;
- W0 = _alpha*_r0;
- _integral_r2_g_dr =
- 1./(4.*pow(w, 2.5))*exp(-w*_r0*_r0)*(
- 2*sqrt(w)*W0 +
- sqrt(M_PI)*exp(W0*W0/w)*(w+2*W0*W0)*(
- 1 - boost::math::erf<double>(-W0/sqrt(w))
- )
- );
- _norm_r2_g_dr = 1./_integral_r2_g_dr;
-}
-
-double RadialGaussian::at(double r) {
- double p = _alpha*(r-_r0)*(r-_r0);
- if (p < 40) return _norm_r2_g2_dr * exp(-p);
- else return 0.0;
-}
-
-
// ======================
// RadialBasisFactory
// ======================
diff --git a/src/radialbasis.hpp b/src/radialbasis.hpp
index 041b9a39a4602723777472c5f7d7f79198ce4048..d1ab523729dbd30dc794638158bb5a6a618f1271 100644
--- a/src/radialbasis.hpp
+++ b/src/radialbasis.hpp
@@ -9,29 +9,17 @@
#include "base/exceptions.hpp"
#include "base/objectfactory.hpp"
#include "options.hpp"
+#include "functions.hpp"
namespace soap {
namespace ub = boost::numeric::ublas;
-class RadialCoefficients : public ub::vector<double>
-{
-public:
- RadialCoefficients(int N) {
- this->resize(N);
- for (int i = 0; i < N; ++i) (*this)[i] = 0.0;
- }
- void set(int n, double c) {
- if (this->checkSize(n)) (*this)[n] = c;
- else throw soap::base::OutOfRange("RadialCoefficients::set");
- }
- double &get(int n) {
- if (this->checkSize(n)) return (*this)[n];
- else throw soap::base::OutOfRange("RadialCoefficients::get");
- }
- bool checkSize(int n) { return this->size() > n; }
-};
-
+// CLASSES IN THIS HEADER
+// RadialBasis
+// RadialBasisGaussian
+// RadialBasisLegendre
+// RadialBasisFactory
class RadialBasis
{
@@ -39,13 +27,13 @@ public:
typedef ub::matrix<double> radcoeff_t;
typedef ub::zero_matrix<double> radcoeff_zero_t;
+ virtual ~RadialBasis() {;}
+
std::string &identify() { return _type; }
const int &N() { return _N; }
- virtual ~RadialBasis() {;}
+
virtual void configure(Options &options);
- virtual RadialCoefficients computeCoefficients(double r);
virtual void computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here);
- virtual RadialCoefficients computeCoefficientsAllZero();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
@@ -57,42 +45,17 @@ public:
}
protected:
- bool _is_ortho;
std::string _type;
int _N;
double _Rc;
int _integration_steps;
std::string _mode; // <- 'equispaced' or 'adaptive'
+ bool _is_ortho;
static constexpr double RADZERO = 1e-10;
};
-struct RadialGaussian
-{
- RadialGaussian(double r0, double sigma);
- RadialGaussian() {;}
- double at(double r);
- double _r0;
- double _sigma;
- double _alpha;
- // Integral S g^2 r^2 dr
- double _integral_r2_g2_dr;
- double _norm_r2_g2_dr;
- // Integral S g r^2 dr
- double _integral_r2_g_dr;
- double _norm_r2_g_dr;
-
- template<class Archive>
- void serialize(Archive &arch, const unsigned int version) {
- arch & _r0;
- arch & _sigma;
- arch & _alpha;
- arch & _norm_r2_g2_dr;
- arch & _norm_r2_g_dr;
- }
-};
-
class RadialBasisGaussian : public RadialBasis
{
public:
@@ -100,22 +63,11 @@ public:
typedef std::vector<RadialGaussian*> basis_t;
typedef basis_t::iterator basis_it_t;
- RadialBasisGaussian() {
- _type = "gaussian";
- _is_ortho = false;
- _sigma = 0.5;
- }
- ~RadialBasisGaussian() {
- this->clear();
- }
- void clear() {
- for (basis_it_t bit=_basis.begin(); bit!=_basis.end(); ++bit) {
- delete *bit;
- }
- _basis.clear();
+ RadialBasisGaussian();
+ ~RadialBasisGaussian() { this->clear();
}
+ void clear();
void configure(Options &options);
- RadialCoefficients computeCoefficients(double r);
void computeCoefficients(double r, double particle_sigma, radcoeff_t &save_here);
template<class Archive>
@@ -131,20 +83,14 @@ public:
protected:
double _sigma;
basis_t _basis;
+
// ORTHONORMALIZATION
ub::matrix<double> _Sij; // Overlap
ub::matrix<double> _Uij; // Cholesky factor
ub::matrix<double> _Tij; // Transformation matrix
-
};
-
-
-
-
-
-
class RadialBasisLegendre : public RadialBasis
{
public:
@@ -155,7 +101,6 @@ public:
};
-
class RadialBasisFactory
: public soap::base::ObjectFactory<std::string, RadialBasis>
{
@@ -183,7 +128,7 @@ inline RadialBasis *RadialBasisFactory::create(const std::string &key) {
}
}
-}
+} /* CLOSE NAMESPACE */
BOOST_CLASS_EXPORT_KEY(soap::RadialBasis);
BOOST_CLASS_EXPORT_KEY(soap::RadialBasisGaussian);
diff --git a/src/spectrum.cpp b/src/spectrum.cpp
index 12b535f053c3c3e254f20347619be865943d7495..0428ca79830f5e36e42a99c9a685ecbaf5c79e86 100644
--- a/src/spectrum.cpp
+++ b/src/spectrum.cpp
@@ -81,9 +81,9 @@ void Spectrum::writeDensityOnGrid(int slot_idx, std::string center_type, std::st
}
// WRITE CUBE FILES
if (atomic_spectrum) {
- atomic_spectrum->getExpansion(density_type)->writeDensityOnGrid(
+ atomic_spectrum->getQnlm(density_type)->writeDensityOnGrid(
"density.expanded.cube", _options, _structure, atomic_spectrum->getCenter(), true);
- atomic_spectrum->getExpansion(density_type)->writeDensityOnGrid(
+ atomic_spectrum->getQnlm(density_type)->writeDensityOnGrid(
"density.explicit.cube", _options, _structure, atomic_spectrum->getCenter(), false);
}
return;
@@ -94,9 +94,6 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
GLOG() << "Compute atomic spectrum for particle " << center->getId()
<< " (type " << center->getType() << ")" << std::endl;
- std::string type_this = center->getType();
-
-
// RadialCoefficients c_n_zero = _radbasis->computeCoefficientsAllZero();
// AngularCoefficients c_lm_zero = _angbasis->computeCoefficientsAllZero();
// BasisCoefficients c_nlm(c_n_zero, c_lm_zero);
@@ -108,11 +105,11 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
Structure::particle_it_t pit;
for (pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
- // FIND DISTANCE & DIRECTION, APPLY CUTOFF
+ // 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;
+ if (weight_scale < 0.) continue; // <- Negative cutoff weight means: skip
vec d = dr/r;
// APPLY WEIGHT IF CENTER
@@ -124,7 +121,7 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center) {
BasisExpansion nb_expansion(this->_basis);
nb_expansion.computeCoefficients(r, d, weight_scale*(*pit)->getWeight(), (*pit)->getSigma());
std::string type_other = (*pit)->getType();
- atomic_spectrum->add(type_other, nb_expansion);
+ atomic_spectrum->addQnlm(type_other, nb_expansion);
// nbhood_expansion->add(nb_expansion);
diff --git a/src/spectrum.hpp b/src/spectrum.hpp
index 5af86db10266db250960ffa8dbe7c6af8a8d85f0..cb9ea9a2043e90853e4ff8c4e4c8d854cc4c38dd 100644
--- a/src/spectrum.hpp
+++ b/src/spectrum.hpp
@@ -13,107 +13,14 @@
#include "options.hpp"
#include "structure.hpp"
#include "basis.hpp"
+#include "power.hpp"
+#include "atomicspectrum.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() :
- _center(NULL),
- _center_pos(vec(0,0,0)),
- _center_type("?"),
- _basis(NULL),
- _expansion_reduced(NULL) { ; }
- ~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; }
-
- template<class Archive>
- void serialize(Archive &arch, const unsigned int version) {
- arch & _center;
- arch & _center_pos;
- arch & _center_type;
- arch & _basis;
- arch & _expansion_reduced;
- return;
- }
-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() {}
-};
-
-
-class TargetDensity
-{
- TargetDensity() {}
-};
-
-class Center
-{
- Center() {}
-};
-
-class Target
-{
- Target() {}
-};
-
-// Need this: Spectrum(System1, System2, options) where Sys1 <> Sources, Sys2 <> Targets
-
-
+// TODO Spectrum(System1, System2, options) where Sys1 <> Sources, Sys2 <> Targets
class Spectrum
{
public:
@@ -162,10 +69,6 @@ private:
};
-
-
-
-
class PairSpectrum
{
public:
diff --git a/src/structure.cpp b/src/structure.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d9f25c983adae643ebb04e301145f33a1cf52422
--- /dev/null
+++ b/src/structure.cpp
@@ -0,0 +1,167 @@
+#include "structure.hpp"
+
+
+namespace soap {
+
+// ========
+// PARTICLE
+// ========
+
+void Particle::null() {
+ _id = -1;
+ _pos = vec(0,0,0);
+ _name = "?";
+ _type_id = -1;
+ _type = "?";
+ _mass = 0.;
+ _weight = 0.;
+ _sigma = 0.;
+ _segment = NULL;
+}
+
+boost::python::numeric::array Particle::getPosNumeric() {
+ boost::python::numeric::array pos(boost::python::make_tuple(_pos.x(), _pos.y(), _pos.z())); return pos;
+}
+
+void Particle::registerPython() {
+ using namespace boost::python;
+ class_<Particle, Particle*>("Particle", init<int>())
+ .add_property("pos", &Particle::getPosNumeric, &Particle::setPosNumeric)
+ .add_property("id", &Particle::getId, &Particle::setId)
+ .add_property("name", make_function(&Particle::getName, copy_non_const()), &Particle::setName)
+ .add_property("type", make_function(&Particle::getType, copy_non_const()), &Particle::setType)
+ .add_property("type_id", make_function(&Particle::getTypeId, copy_non_const()), &Particle::setTypeId)
+ .add_property("mass", make_function(&Particle::getMass, copy_non_const()), &Particle::setMass)
+ .add_property("weight", make_function(&Particle::getWeight, copy_non_const()), &Particle::setWeight)
+ .add_property("sigma", make_function(&Particle::getSigma, copy_non_const()), &Particle::setSigma);
+}
+
+// =======
+// SEGMENT
+// =======
+
+Particle &Segment::addParticle(Particle *new_part) {
+ _particles.push_back(new_part);
+ return *new_part;
+}
+
+void Segment::registerPython() {
+ using namespace boost::python;
+ class_<Segment>("Segment", init<int>())
+ .def("addParticle", &Segment::addParticle, return_value_policy<reference_existing_object>());
+}
+
+
+// =========
+// STRUCTURE
+// =========
+
+Structure::Structure(std::string label) {
+ this->null();
+ _label = label;
+
+}
+
+Structure::Structure() {
+ this->null();
+}
+
+void Structure::null() {
+ _id = -1;
+ _label = "?";
+ _box = NULL;
+ _center = NULL;
+ _has_center = false;
+ matrix box;
+ box.ZeroMatrix();
+ this->setBoundary(box);
+}
+
+Structure::~Structure() {
+ delete _box;
+ _box = NULL;
+ segment_it_t sit;
+ for (sit = _segments.begin(); sit != _segments.end(); ++sit) {
+ delete *sit;
+ }
+ _segments.clear();
+ particle_it_t pit;
+ for (pit = _particles.begin(); pit != _particles.end(); ++pit) {
+ delete *pit;
+ }
+ _particles.clear();
+}
+
+Segment &Structure::addSegment() {
+ int id = _segments.size()+1;
+ Segment *new_seg = new Segment(id);
+ _segments.push_back(new_seg);
+ return *new_seg;
+}
+
+Particle &Structure::addParticle(Segment &seg) {
+ int id = _particles.size()+1;
+ Particle *new_part = new Particle(id);
+ _particles.push_back(new_part);
+ seg.addParticle(new_part);
+ return *new_part;
+}
+
+boost::python::numeric::array Structure::connectNumeric(
+ const boost::python::numeric::array &a1,
+ const boost::python::numeric::array &a2) {
+ vec r1(a1);
+ vec r2(a2);
+ vec dr = this->connect(r1, r2);
+ return boost::python::numeric::array(boost::python::make_tuple(
+ dr.x(), dr.y(), dr.z()));
+}
+
+void Structure::setBoundary(const matrix &box) {
+ delete _box;
+ if(box.get(0,0)==0 && box.get(0,1)==0 && box.get(0,2)==0 &&
+ box.get(1,0)==0 && box.get(1,1)==0 && box.get(1,2)==0 &&
+ box.get(2,0)==0 && box.get(2,1)==0 && box.get(2,2)==0) {
+ _box = new BoundaryOpen(box);
+ }
+ else if(box.get(0,1)==0 && box.get(0,2)==0 &&
+ box.get(1,0)==0 && box.get(1,2)==0 &&
+ box.get(2,0)==0 && box.get(2,1)==0) {
+ _box = new BoundaryOrthorhombic(box);
+ }
+ else {
+ _box = new BoundaryTriclinic(box);
+ }
+ return;
+}
+
+void Structure::setBoundaryNumeric(const boost::python::numeric::array &m) {
+ matrix box(m);
+ this->setBoundary(box);
+}
+
+boost::python::numeric::array Structure::getBoundaryNumeric() {
+ matrix box = _box->getBox();
+ boost::python::numeric::array box_np(boost::python::make_tuple(
+ box.get(0,0), box.get(0,1), box.get(0,2),
+ box.get(1,0), box.get(1,1), box.get(1,2),
+ box.get(2,0), box.get(2,1), box.get(2,2)));
+ box_np.resize(3,3);
+ return box_np;
+}
+
+void Structure::registerPython() {
+ using namespace boost::python;
+ class_<Structure>("Structure", init<std::string>())
+ .def("addSegment", &Structure::addSegment, 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))
+ .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>());
+}
+
+} /* CLOSE NAMESPACE */
diff --git a/src/structure.hpp b/src/structure.hpp
index 523e0b5cdf1d562a8c3de07a1cd9740394557e7b..2834d2caa6bce587af66fafa26d0f7497f41770c 100644
--- a/src/structure.hpp
+++ b/src/structure.hpp
@@ -18,58 +18,40 @@ class Structure;
class Particle
{
public:
- Particle(int id) :
- _id(id), _pos(vec(0,0,0)), _name(""),
- _type_id(-1), _type(""), _mass(.0),
- _weight(.0), _sigma(.0), _segment(NULL) { ; }
- Particle() :
- _id(-1), _pos(vec(0,0,0)), _name(""),
- _type_id(-1), _type(""), _mass(.0),
- _weight(.0), _sigma(.0), _segment(NULL) { ; }
- ~Particle() {
- ;
- }
-
- static void registerPython() {
- using namespace boost::python;
- class_<Particle, Particle*>("Particle", init<int>())
- .add_property("pos", &Particle::getPosNumeric, &Particle::setPosNumeric)
- .add_property("id", &Particle::getId, &Particle::setId)
- .add_property("name", make_function(&Particle::getName, copy_non_const()), &Particle::setName)
- .add_property("type", make_function(&Particle::getType, copy_non_const()), &Particle::setType)
- .add_property("type_id", make_function(&Particle::getTypeId, copy_non_const()), &Particle::setTypeId)
- .add_property("mass", make_function(&Particle::getMass, copy_non_const()), &Particle::setMass)
- .add_property("weight", make_function(&Particle::getWeight, copy_non_const()), &Particle::setWeight)
- .add_property("sigma", make_function(&Particle::getSigma, copy_non_const()), &Particle::setSigma);
- }
-
+ Particle(int id) { this->null(); _id = id; }
+ Particle() { this->null(); }
+ ~Particle() {;}
+ void null();
+ // Position
void setPos(vec &pos) { _pos = pos; }
void setPos(double x, double y, double z) { _pos = vec(x,y,z); }
- void setPosNumeric(const boost::python::numeric::array &pos) { _pos = vec(pos); }
vec &getPos() { return _pos; }
- boost::python::numeric::array getPosNumeric() { boost::python::numeric::array pos(boost::python::make_tuple(_pos.x(), _pos.y(), _pos.z())); return pos; }
-
+ void setPosNumeric(const boost::python::numeric::array &pos) { _pos = vec(pos); }
+ boost::python::numeric::array getPosNumeric();
+ // 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 Id
void setTypeId(int id) { _type_id = id; }
int &getTypeId() { return _type_id; }
-
+ // Type
void setType(std::string type) { _type = type; }
std::string &getType() { return _type; }
-
+ // Mass
void setMass(double mass) { _mass = mass; }
double &getMass() { return _mass; }
-
+ // Weight
void setWeight(double weight) { _weight = weight; }
double &getWeight() { return _weight; }
-
+ // Sigma
void setSigma(double sigma) { _sigma = sigma; }
double &getSigma() { return _sigma; }
+ static void registerPython();
+
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
arch & _segment;
@@ -86,7 +68,7 @@ public:
private:
Segment *_segment;
- // Book-keeping
+ // Labels
int _id;
std::string _name;
int _type_id;
@@ -102,20 +84,15 @@ private:
class Segment
{
public:
+ typedef std::vector<Particle*> particle_array_t;
+
Segment(int id) : _id(id) { ; }
Segment() : _id(-1) { ; }
- ~Segment() {
- _particles.clear();
- }
- static void registerPython() {
- using namespace boost::python;
- class_<Segment>("Segment", init<int>())
- .def("addParticle", &Segment::addParticle, return_value_policy<reference_existing_object>());
- }
- Particle &addParticle(Particle *new_part) {
- _particles.push_back(new_part);
- return *new_part;
- }
+ ~Segment() { _particles.clear(); }
+
+ Particle &addParticle(Particle *new_part);
+
+ static void registerPython();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
@@ -126,114 +103,44 @@ public:
private:
int _id;
- std::vector<Particle*> _particles;
+ particle_array_t _particles;
};
class Structure
{
public:
+ typedef std::vector<Particle*> particle_array_t;
typedef std::vector<Particle*>::iterator particle_it_t;
+ typedef std::vector<Segment*> segment_array_t;
+ typedef std::vector<Segment*>::iterator segment_it_t;
+
+ Structure(std::string label);
+ Structure();
+ ~Structure();
+ void null();
+
+ particle_array_t &particles() { return _particles; }
+ particle_it_t beginParticles() { return _particles.begin(); }
+ particle_it_t endParticles() { return _particles.end(); }
- Structure(std::string label) : _id(-1), _label(label), _box(NULL) {
- matrix box;
- box.ZeroMatrix();
- this->setBoundary(box);
- }
- Structure() : _id(-1), _label("?"), _box(NULL) {
- matrix box;
- box.ZeroMatrix();
- this->setBoundary(box);
- }
- ~Structure() {
- delete _box;
- _box = NULL;
- std::vector<Segment*>::iterator sit;
- for (sit = _segments.begin(); sit != _segments.end(); ++sit) {
- delete *sit;
- }
- _segments.clear();
- std::vector<Particle*>::iterator pit;
- for (pit = _particles.begin(); pit != _particles.end(); ++pit) {
- delete *pit;
- }
- _particles.clear();
- }
- static void registerPython() {
- using namespace boost::python;
- class_<Structure>("Structure", init<std::string>())
- .def("addSegment", &Structure::addSegment, 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))
- .def("connect", &Structure::connectNumeric)
- .add_property("box", &Structure::getBoundaryNumeric, &Structure::setBoundaryNumeric)
- .add_property("label", make_function(&Structure::getLabel, copy_non_const()), &Structure::setLabel);
- class_< std::vector<Particle*> >("ParticleContainer")
- .def(vector_indexing_suite<std::vector<Particle*> >());
- }
- std::vector<Particle*> &particles() { return _particles; }
- std::vector<Particle*>::iterator beginParticles() { return _particles.begin(); }
- std::vector<Particle*>::iterator endParticles() { return _particles.end(); }
std::string &getLabel() { return _label; }
void setLabel(std::string label) { _label = label; }
+
// PARTICLE CREATION & INTERFACE
- Segment &addSegment() {
- int id = _segments.size()+1;
- Segment *new_seg = new Segment(id);
- _segments.push_back(new_seg);
- return *new_seg;
- }
- Particle &addParticle(Segment &seg) {
- int id = _particles.size()+1;
- Particle *new_part = new Particle(id);
- _particles.push_back(new_part);
- seg.addParticle(new_part);
- return *new_part;
- }
+ Segment &addSegment();
+ Particle &addParticle(Segment &seg);
+
// BOUNDARY CREATION & INTERFACE
- vec connect(const vec &r1, const vec &r2) {
- return _box->connect(r1, r2); // Points from r1 to r2
- }
- boost::python::numeric::array connectNumeric(
- const boost::python::numeric::array &a1,
- const boost::python::numeric::array &a2) {
- vec r1(a1);
- vec r2(a2);
- vec dr = this->connect(r1, r2);
- return boost::python::numeric::array(boost::python::make_tuple(
- dr.x(), dr.y(), dr.z()));
- }
- void setBoundary(const matrix &box) {
- delete _box;
- if(box.get(0,0)==0 && box.get(0,1)==0 && box.get(0,2)==0 &&
- box.get(1,0)==0 && box.get(1,1)==0 && box.get(1,2)==0 &&
- box.get(2,0)==0 && box.get(2,1)==0 && box.get(2,2)==0) {
- _box = new BoundaryOpen(box);
- }
- else if(box.get(0,1)==0 && box.get(0,2)==0 &&
- box.get(1,0)==0 && box.get(1,2)==0 &&
- box.get(2,0)==0 && box.get(2,1)==0) {
- _box = new BoundaryOrthorhombic(box);
- }
- else {
- _box = new BoundaryTriclinic(box);
- }
- return;
- }
- void setBoundaryNumeric(const boost::python::numeric::array &m) {
- matrix box(m);
- this->setBoundary(box);
- }
- boost::python::numeric::array getBoundaryNumeric() {
- matrix box = _box->getBox();
- boost::python::numeric::array box_np(boost::python::make_tuple(
- box.get(0,0), box.get(0,1), box.get(0,2),
- box.get(1,0), box.get(1,1), box.get(1,2),
- box.get(2,0), box.get(2,1), box.get(2,2)));
- box_np.resize(3,3);
- return box_np;
- }
+ void setBoundary(const matrix &box);
+ vec connect(const vec &r1, const vec &r2) { return _box->connect(r1, r2); /* 1->2 */ }
+ void setBoundaryNumeric(const boost::python::numeric::array &m);
+ boost::python::numeric::array getBoundaryNumeric();
+ boost::python::numeric::array connectNumeric(
+ const boost::python::numeric::array &a1,
+ const boost::python::numeric::array &a2);
+
+ static void registerPython();
template<class Archive>
void serialize(Archive &arch, const unsigned int version) {
@@ -248,8 +155,11 @@ public:
private:
int _id;
std::string _label;
- std::vector<Segment*> _segments;
- std::vector<Particle*> _particles;
+ segment_array_t _segments;
+ particle_array_t _particles;
+
+ Particle* _center;
+ bool _has_center;
Boundary *_box;
};
diff --git a/src/tools/loadwrite.py b/src/tools/loadwrite.py
index 0a4caba7438d8eb878caa4500cd1fba56ae33d63..a4d43dd3195c8f31c484774f58fc396e9baf73f3 100644
--- a/src/tools/loadwrite.py
+++ b/src/tools/loadwrite.py
@@ -20,7 +20,11 @@ def write_xyz(xyz_file, structure):
def ase_load_all(folder):
cwd = os.getcwd()
os.chdir(folder)
- config_files = sorted(os.listdir('./'))
+ config_listdir = sorted(os.listdir('./'))
+ config_files = []
+ for config_file in config_listdir:
+ if os.path.isfile(config_file):
+ config_files.append(config_file)
ase_config_list = AseConfigList(config_files)
os.chdir(cwd)
return ase_config_list
diff --git a/test/test.py b/test/test.py
index 895cbd59f9ea7165b801375fd52b0dd57eaabc1f..228c983bed30dd4b17b4d8b2b84763d6938899df 100755
--- a/test/test.py
+++ b/test/test.py
@@ -20,27 +20,32 @@ element_mass = {
}
ase_config_list = soap.tools.ase_load_all('configs')
-config = ase_config_list[1]
+for config in ase_config_list:
+ print config.config_file
+config = ase_config_list[4]
osio << config.atoms << endl
+sigma = 0.5
+
# 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.sigma', 0.5)
+options.set('radialbasis.sigma', sigma)
options.set('radialbasis.integration_steps', 15)
#options.set('radialbasis.N', 9)
-options.set('radialcutoff.Rc', 4.)
+options.set('radialcutoff.Rc', 6.8)
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('radialcutoff.center_weight', -7.)
options.set('angularbasis.type', 'spherical-harmonic')
options.set('angularbasis.L', 6)
-#options.set('angularbasis.L', 6)
+#options.set('angularbasis.L', 12)
+#options.set('densitygrid.N', 20)
options.set('densitygrid.N', 20)
options.set('densitygrid.dx', 0.15)
@@ -50,75 +55,24 @@ structure = soap.tools.setup_structure_ase(config.config_file, config.atoms)
osio << osio.mg << structure.label << endl
for atom in structure:
- atom.sigma = 0.5 # 0.4*element_vdw_radius[atom.type]
+ atom.sigma = 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
+ #if atom.id > 60: atom.weight *= -1
# COMPUTE SPECTRUM
spectrum = soap.Spectrum(structure, options)
spectrum.compute()
-#spectrum.writeDensityOnGrid(1, "C", "")
+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()
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-structure = soap.Structure(config.config_file)
-
-# BOX TESTS
-osio << osio.mg << "Box checks ..." << endl
-r1 = np.array([0.2,0.3,0.4])
-r2 = np.array([0.8,1.0,1.2])
-# OPEN BOX
-structure.box = np.array([[0,0,0],[0,0,0],[0,0,0]])
-print "Connect", r1, r2, " => ", structure.connect(r1,r2)
-# CUBIC BOX
-a = np.array([0.9,0,0])
-b = np.array([0,0.7,0])
-c = np.array([0,0,0.5])
-structure.box = np.array([a,b,c])
-print "Connect", r1, r2, " => ", structure.connect(r1,r2)
-
-
-segment = structure.addSegment()
-particle = structure.addParticle(segment)
-particle.pos = np.array([-1,0,1]) # vec(0,1,-1)
-print type(particle.pos), particle.pos
-particle.mass = 12.
-