FT init. Special functions.

src/soap/bindings.cpp

 #include "bindings.hpp"
 #include "coulomb.hpp"
+#include "fieldtensor.hpp"
 
 namespace soap {
 
@@ -23,6 +24,8 @@ BOOST_PYTHON_MODULE(_soapxx)
     soap::EnergySpectrum::registerPython();
     soap::HierarchicalCoulomb::registerPython();
     soap::AtomicSpectrumHC::registerPython();
+    soap::FTSpectrum::registerPython();
+    soap::AtomicSpectrumFT::registerPython();
 
     soap::RadialBasisFactory::registerAll();
     soap::AngularBasisFactory::registerAll();

src/soap/fieldtensor.cpp

+#include "soap/fieldtensor.hpp"
+#include "soap/functions.hpp"
+#include <boost/math/special_functions/legendre.hpp>
+
+namespace soap {
+
+namespace ub = boost::numeric::ublas;
+
+// ============================
+// Spectrum::FTSpectrum::Atomic
+// ============================
+
+const std::string AtomicSpectrumFT::_numpy_t = "float64";
+
+AtomicSpectrumFT::AtomicSpectrumFT(Particle *center, int S)
+    : _center(center), _S(S) {
+    this->_Q0 = coeff_zero_t(S,1);
+    this->_Q1 = coeff_zero_t(S,S);
+    this->_Q2 = coeff_zero_t(S,S*S);
+    this->_Q3 = coeff_zero_t(S,S*S*S);
+
+    this->_s = center->getTypeId()-1;
+    assert(_s >= 0 && "Type-IDs should start from 1");
+}
+
+AtomicSpectrumFT::~AtomicSpectrumFT() {
+    _Q0.clear();
+    _Q1.clear();
+    _Q2.clear();
+    _Q3.clear();
+}
+
+void AtomicSpectrumFT::registerPython() {
+    using namespace boost::python;
+    class_<AtomicSpectrumFT, AtomicSpectrumFT*>("AtomicSpectrumFT", init<Particle*, int>())
+        .def("getCenter", &AtomicSpectrumFT::getCenter, return_value_policy<reference_existing_object>());
+}
+
+// ===================
+// FieldTensorSpectrum
+// ===================
+
+FTSpectrum::FTSpectrum(Structure &structure, Options &options)
+    : _structure(&structure), _options(&options) {
+    return;
+}
+
+FTSpectrum::~FTSpectrum() {
+    return;
+}
+
+void FTSpectrum::compute() {
+
+    int L = 3;
+
+    GLOG() << "Legendre ..." << std::endl;
+    std::vector<double> plm;
+    calculate_legendre_plm(L, 0.3, plm);
+    for (int l = 0; l <= L; ++l) {
+        for (int m = 0; m <= l; ++m) {
+            int lm = l*(l+1)/2 + m;
+            GLOG() << l << " " << m << " " << plm[lm] << std::endl;
+        }
+    }
+
+    GLOG() << "Factorial ..." << std::endl;
+    for (int n = 0; n < 16; ++n) {
+        GLOG() << n << "! = " << factorial(n) << std::endl;
+    }
+
+    GLOG() << "Solid harmonics ..." << std::endl;
+    std::vector<std::complex<double>> rlm;
+    std::vector<std::complex<double>> ilm;
+    double phi = 0.6;
+    double theta = 0.7;
+    double sp = std::sin(phi);
+    double st = std::sin(theta);
+    double cp = std::cos(phi);
+    double ct = std::cos(theta);
+    vec d = vec(st*cp, st*sp, ct);
+    double r = 0.5;
+    calculate_solidharm_rlm_ilm(d, r, L, rlm, ilm);
+    for (int l = 0; l <= L; ++l) {
+        for (int m = -l; m <= l; ++m) {
+            int lm = l*l+l+m;
+            GLOG() << l << " " << m << " " << rlm[lm] << "  " << ilm[lm] << std::endl;
+        }
+    }
+
+    return;
+}
+
+void FTSpectrum::registerPython() {
+    using namespace boost::python;
+    class_<FTSpectrum>("FTSpectrum", init<Structure &, Options &>())
+        .def("__iter__", range<return_value_policy<reference_existing_object> >(&FTSpectrum::beginAtomic, &FTSpectrum::endAtomic))
+        .def("compute", &FTSpectrum::compute);
+    return;
+}
+
+
+}

src/soap/fieldtensor.hpp

+#ifndef _SOAP_FIELDTENSOR_HPP
+#define _SOAP_FIELDTENSOR_HPP
+
+#include "soap/structure.hpp"
+#include "soap/options.hpp"
+#include "soap/cutoff.hpp"
+
+namespace soap {
+
+class AtomicSpectrumFT
+{
+public:
+    typedef double dtype_t;
+    typedef ub::matrix<dtype_t> coeff_t;
+    typedef ub::zero_matrix<dtype_t> coeff_zero_t;
+    static const std::string _numpy_t;
+
+    AtomicSpectrumFT(Particle *center, int S);
+   ~AtomicSpectrumFT();
+    Particle *getCenter() { return _center; }
+    int getTypeIdx() { return _s; }
+
+    static void registerPython();
+
+    coeff_t _Q0;
+    coeff_t _Q1;
+    coeff_t _Q2;
+    coeff_t _Q3;
+
+private:
+    Particle *_center;
+    int _S; // <- Number of distinct types
+    int _s; // <- Type index
+};
+
+class FTSpectrum
+{
+public:
+    typedef AtomicSpectrumFT atomic_t;
+    typedef std::vector<AtomicSpectrumFT*> atomic_array_t;
+	typedef std::vector<AtomicSpectrumFT*>::iterator atomic_it_t;
+
+    FTSpectrum(Structure &structure, Options &options);
+   ~FTSpectrum();
+    void compute();
+    atomic_it_t beginAtomic() { return _atomic_array.begin(); }
+    atomic_it_t endAtomic() { return _atomic_array.end(); }
+
+    static void registerPython();
+
+private:
+    Structure *_structure;
+    Options *_options;
+    CutoffFunction *_cutoff;
+
+    atomic_array_t _atomic_array;
+};
+
+
+
+}
+
+#endif

src/soap/functions.cpp

@@ -247,8 +247,118 @@ std::complex<double> pow_nnan(std::complex<double> z, double a) {
     }
 }
 
-int factorial(int n) {
-    return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
+const int FACTORIAL_CACHE_SIZE = 16;
+const long int FACTORIAL_CACHE[] = {
+  1,
+  1, 2, 6,
+  24, 120, 720,
+  5040, 40320, 362880,
+  3628800, 39916800, 479001600,
+  6227020800, 87178291200, 1307674368000 };
+
+long int factorial(int x) {
+    if (x < FACTORIAL_CACHE_SIZE) {
+        return FACTORIAL_CACHE[x];
+    }
+    else {
+        assert(false && "Computing factorial from scratch - not desirable");
+        long int s = 1.0;
+        for (int n = 2; n <= x; n++) {
+          s *= n;
+        }
+        return s;
+    }
+}
+
+void calculate_solidharm_rlm_ilm(
+        vec d,
+        double r,
+        int L,
+        std::vector<std::complex<double>> &rlm,
+        std::vector<std::complex<double>> &ilm) {
+
+    std::vector<double> plm;
+    double theta = 0.0;
+    double phi = 0.0;
+
+    // Compute Associated Legendre Polynomials
+    if (r < 1e-10) { // TODO Define SPACE-QUANTUM
+        assert(false && "Should not calculate irregular spherical harmonics with r=0");
+        plm.resize(L*(L+1)/2+L+1, 0.0);
+        plm[0] = boost::math::legendre_p(0, 0, 0);
+    }
+    else {
+		theta = acos(d.getZ());
+		phi = atan2(d.getY(), d.getX());
+		if (phi < 0.) phi += 2*M_PI; // <- Shift [-pi, -0] to [pi, 2*pi]
+		calculate_legendre_plm(L, d.getZ(), plm);
+    }
+
+    std::cout << "phi = " << phi << std::endl;
+
+    rlm.clear();
+    rlm.resize((L+1)*(L+1), 0.0);
+    for (int l = 0; l <= L; ++l) {
+        for (int m = 0; m <= l; ++m) {
+            rlm[l*l+l+m] =
+                  pow(-1, l-m)
+                * pow(r, l)
+                / factorial(l+m)
+                * std::exp(std::complex<double>(0.,m*phi))
+                * plm[l*(l+1)/2+m];
+        }
+        for (int m = -l; m < 0; ++m) {
+            rlm[l*l+l+m] = pow(-1, m)*std::conj(rlm[l*l+l-m]);
+        }
+    }
+
+    ilm.clear();
+    ilm.resize((L+1)*(L+1), 0.0);
+    for (int l = 0; l <= L; ++l) {
+        for (int m = 0; m <= l; ++m) {
+            ilm[l*l+l+m] =
+                  pow(-1, l-m)
+                * pow(r, -l-1)
+                * factorial(l-m)
+                * std::exp(std::complex<double>(0.,m*phi))
+                * plm[l*(l+1)/2+m];
+        }
+        for (int m = -l; m < 0; ++m) {
+            ilm[l*l+l+m] = pow(-1, m)*std::conj(ilm[l*l+l-m]);
+        }
+    }
+
+    return;
+}
+
+void calculate_legendre_plm(int L, double x, std::vector<double> &plm_out) {
+    // See 'Numerical Recipes' - Chapter on Special Functions
+    plm_out.clear();
+    plm_out.resize(L*(L+1)/2+L+1, 0.0);
+    // Pll / Pmm
+    for (int l = 0; l <= L; ++l) {
+        plm_out[l*(l+1)/2+l] = boost::math::legendre_p(l, l, x);
+        //GLOG() << l << ":" << l << std::endl;
+    }
+    // Pmm => Pm+1m
+    for (int l = 0; l < L; ++l) {
+        int m = l;
+        int ll = l+1;
+        plm_out[ll*(ll+1)/2+l] = x*(2*m+1)*plm_out[l*(l+1)/2+m];
+        //GLOG() << ll << ":" << l << std::endl;
+    }
+    // Pl-1m,Pl-2m => Plm
+    for (int m = 0; m <= L; ++m) {
+        for (int l2 = m+2; l2 <= L; ++l2) {
+            int l0 = l2-2;
+            int l1 = l2-1;
+            int lm0 = l0*(l0+1)/2+m;
+            int lm1 = l1*(l1+1)/2+m;
+            int lm2 = l2*(l2+1)/2+m;
+            //GLOG() << l2 << ":" << m << " from " << lm1 << "," << lm0 << std::endl;
+            plm_out[lm2] = (x*(2*l2-1)*plm_out[lm1] - (l2+m-1)*plm_out[lm0])/(l2-m);
+        }
+    }
 }
 
 } /* CLOSE NAMESPACE */

src/soap/functions.hpp

@@ -67,7 +67,22 @@ struct GradSphericalYlm
 };
 
 std::complex<double> pow_nnan(std::complex<double> z, double a);
-int factorial(int n);
+
+extern const int FACTORIAL_CACHE_SIZE;
+extern const long int FACTORIAL_CACHE[];
+long int factorial(int n);
+
+void calculate_solidharm_rlm_ilm(
+        vec d,
+        double r,
+        int L,
+        std::vector<std::complex<double>> &rlm,
+        std::vector<std::complex<double>> &ilm);
+
+void calculate_legendre_plm(
+        int L,
+        double x,
+        std::vector<double> &plm_out);
 
 }
 

src/soap/spectrum.cpp

@@ -414,4 +414,3 @@ void Spectrum::registerPython() {
  */
 
 }
-