Some code hygiene.

src/angularbasis.cpp

 #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 */
 

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

src/atomicspectrum.cpp

+#include <fstream>
+#include <boost/format.hpp>
+#include <boost/archive/binary_oarchive.hpp>
+#include <boost/archive/binary_iarchive.hpp>
+
+#include "atomicspectrum.hpp"
+
+namespace soap {
+
+}
+

src/atomicspectrum.hpp

+#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_ */

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;
+}
+
 }
 

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) {