From a42be09a7b11c194f70100a91e0558b97d0f9fe7 Mon Sep 17 00:00:00 2001
From: Carl Poelking <cp605@cam.ac.uk>
Date: Tue, 20 Dec 2016 20:04:07 +0000
Subject: [PATCH] Tlmlm framework.
---
src/soap/fieldtensor.cpp | 252 ++++++++++++++++++++++++++++++++++++++-
src/soap/functions.cpp | 64 ++++++++++
src/soap/functions.hpp | 10 ++
3 files changed, 323 insertions(+), 3 deletions(-)
diff --git a/src/soap/fieldtensor.cpp b/src/soap/fieldtensor.cpp
index b8f946b..00b8d02 100644
--- a/src/soap/fieldtensor.cpp
+++ b/src/soap/fieldtensor.cpp
@@ -40,8 +40,169 @@ void AtomicSpectrumFT::registerPython() {
// FieldTensorSpectrum
// ===================
+std::vector<double> calculate_Alm(int L) {
+ // A_{lm} = \sqrt{(l+m)!*(l-m)!}
+ std::vector<double> out;
+ out.resize((L+1)*(L+1), 0.0);
+ for (int l = 0; l <= L; ++l) {
+ for (int m = -l; m <= L; ++m) {
+ out[l*l+l+m] = sqrt(factorial(l+m)*factorial(l-m));
+ }
+ }
+ return out;
+}
+
+std::vector<double> calculate_Blm(int L) {
+ // B_{lm} = A_{lm} / (2l-1)!!
+ std::vector<double> out;
+ out.resize((L+1)*(L+1), 0.0);
+ for (int l = 0; l <= L; ++l) {
+ for (int m = -l; m <= L; ++m) {
+ out[l*l+l+m] = sqrt(factorial(l+m)*factorial(l-m))/factorial2(2*l-1);
+ }
+ }
+ return out;
+}
+
+void calculate_Fl(double r, double a, int L, std::vector<double> &fl) {
+ // See Elking et al.: "Gaussian Multipole Model", JCTC (2010), Eq. 15a++
+ fl.clear();
+ fl.resize(L+1, 0.0);
+ if (r < 1e-10) { // TODO Define space quantum
+ double a2 = a*a;
+ double a_l = 1./a;
+ int two_l = 1;
+ int sign_l = -1;
+ double rsqrt_pi = 1./sqrt(M_PI);
+ for (int l = 0; l <= L; ++l) {
+ a_l *= a2;
+ two_l *= 2;
+ sign_l *= -1;
+ fl[l] = sign_l*two_l*a_l*rsqrt_pi/(2*l+1);
+ }
+ }
+ else {
+ // Initialise
+ fl[0] = std::erf(a*r)/r;
+ double r2 = r*r;
+ double a2 = a*a;
+ double r_sqrt_pi_exp_a2r2 = 1./sqrt(M_PI) * exp(-a2*r2);
+ double al = 1./a;
+ double tl = 1.;
+ // Compute
+ for (int l = 1; l <= L; ++l) {
+ al *= a2;
+ tl *= 2;
+ fl[l] = -1./r2 * ( (2*l-1)*fl[l-1] + pow(-1,l)*tl*al*r_sqrt_pi_exp_a2r2 );
+ }
+ }
+
+ return;
+}
+
+class Tlmlm
+{
+public:
+ typedef std::complex<double> dtype_t;
+ typedef ub::matrix< dtype_t > coeff_t;
+ typedef ub::zero_matrix< dtype_t > coeff_zero_t;
+
+ Tlmlm(int L) {
+ Alm = calculate_Alm(L);
+ Blm = calculate_Blm(L);
+ }
+
+ void computeTlmlm(vec d12, double r12, double a12, int L1, int L2, coeff_t &T) {
+ T.clear();
+ T = coeff_zero_t((L1+1)*(L1+1), (L2+1)*(L2+1));
+
+ int Lg = (L1 > L2) ? L1 : L2;
+ int L12 = L1+L2;
+ // Compute Rlm's
+ std::vector<std::complex<double>> Rlm;
+ calculate_solidharm_Rlm(d12, r12, Lg, Rlm);
+ // Compute Fl's
+ std::vector<double> Fl;
+ calculate_Fl(r12, a12, L12, Fl);
+
+ for (int l1 = 0; l1 <= L1; ++l1) {
+ for (int m1 = -l1; m1 <= l1; ++m1) {
+ int lm1 = l1*l1+l1+m1;
+ for (int l2 = 0; l2 <= L2; ++l2) {
+ for (int m2 = -l2; m2 <= l2; ++m2) {
+ int lm2 = l2*l2+l2+m2;
+
+ std::complex<double> tlm1lm2 = 0.0;
+ int lg = (l1 > l2) ? l1 : l2;
+
+ for (int l = 0; l <= lg; ++l) {
+ for (int m = -l; m <= l; ++m) {
+
+ if (std::abs(m2+m) > l2-l) continue;
+ if (std::abs(m1-m) > l1-l) continue;
+ int lm = l*l+l+m;
+ tlm1lm2 +=
+ pow(-1, l1+m)*factorial2(2*l-1)/(Alm[lm]*Alm[lm])
+ * std::conj( Rlm[ (l2-l)*(l2-l) + (l2-l) + m2+m ] )
+ * std::conj( Rlm[ (l1-l)*(l1-l) + (l1-l) + m1-m ] )
+ * Fl[ l1+l2-l ];
+
+ }} // l m
+
+ tlm1lm2 *=
+ Alm[lm1]*Alm[lm1] * Alm[lm2]*Alm[lm2]
+ / factorial2(2*l1-1)*factorial2(2*l2-1);
+
+ // TODO Assert .imag() == zero // TODO Does not actually need to be zero!?
+ T(lm1,lm2) = tlm1lm2;
+ }}}} // l1 m1 l2 m2
+ return;
+ }
+
+ double computeTl1m1l2m2(double d12, double r12, double a12, int l1, int l2) {
+
+ }
+private:
+ std::vector<double> Alm;
+ std::vector<double> Blm;
+};
+
+
+void calculate_r_dr_erfar_r(double r, double a, int L, bool normalise, std::vector<double> &cl) {
+ // Derivatives (1/r d/dr)^l (erf(ar)/r) for 0 <= l <= L
+ cl.clear();
+ cl.resize(L+1, 0.0);
+ // Initialise
+ cl[0] = std::erf(a*r)/r;
+ double r2 = r*r;
+ double a2 = a*a;
+ double r_sqrt_pi_exp_a2r2 = 1./sqrt(M_PI) * exp(-a2*r2);
+ double al = 1./a;
+ double tl = 1.;
+ // Compute
+ for (int l = 1; l <= L; ++l) {
+ al *= a2;
+ tl *= 2;
+ cl[l] = 1./r2 * ( (2*l-1)*cl[l-1] - tl*al*r_sqrt_pi_exp_a2r2 );
+ }
+
+ if (normalise) {
+ double rl = r;
+ cl[0] *= rl; // * factorial2(-1), which equals 1
+ for (int l = 1; l <= L; ++l) {
+ rl *= r2;
+ cl[l] *= rl / factorial2(2*l-1);
+ }
+ }
+
+ return;
+}
+
+
FTSpectrum::FTSpectrum(Structure &structure, Options &options)
: _structure(&structure), _options(&options) {
+ _cutoff = CutoffFunctionOutlet().create(_options->get<std::string>("radialcutoff.type"));
+ _cutoff->configure(*_options);
return;
}
@@ -50,10 +211,93 @@ FTSpectrum::~FTSpectrum() {
}
void FTSpectrum::compute() {
+ GLOG() << "Computing FTSpectrum ..." << std::endl;
+
+ int L = 3; // TODO
+ int S = 4;
+
+ // CREATE ATOMIC SPECTRA
+ for (auto it = _atomic_array.begin(); it != _atomic_array.end(); ++it) {
+ delete *it;
+ }
+ _atomic_array.clear();
+ for (auto pit = _structure->beginParticles(); pit != _structure->endParticles(); ++pit) {
+ atomic_t *new_atomic = new atomic_t(*pit, S);
+ _atomic_array.push_back(new_atomic);
+ }
+
+ for (auto it1 = beginAtomic(); it1 != endAtomic(); ++it1) {
+ atomic_t *a = *it1;
+ int sa = a->getTypeIdx();
+ vec ra = a->getCenter()->getPos();
+ for (auto it2 = beginAtomic(); it2 != endAtomic(); ++it2) {
+ atomic_t *b = *it2;
+ int sb = b->getTypeIdx();
+ vec rb = b->getCenter()->getPos();
+ // Apply weight function
+ vec dr_ab = _structure->connect(ra, rb);
+ double R_ab = soap::linalg::abs(dr_ab);
+ vec d_ab = dr_ab/R_ab;
+ if (! _cutoff->isWithinCutoff(R_ab)) continue;
+ double w_ab = _cutoff->calculateWeight(R_ab);
+ GLOG() << " " << a->getCenter()->getId() << ":" << b->getCenter()->getId() << " R=" << R_ab << " w=" << w_ab << std::endl;
+ // Interact
+ // ...
+ }
+ }
- int L = 3;
+ GLOG() << "Factorial2 ..." << std::endl;
+ for (int n = 0; n < 16; ++n) {
+ GLOG() << n << "!! = " << factorial2(n) << std::endl;
+ }
+
+ GLOG() << "Radial damping functions ..." << std::endl;
+ std::vector<double> cl;
+ calculate_Fl(0., 0.5, 5, cl);
+ for (int l = 0; l <= 4; ++l) {
+ GLOG() << l << " fl = " << cl[l] << std::endl;
+ }
+ calculate_Fl(0.5, 0.5, 5, cl);
+ for (int l = 0; l <= 4; ++l) {
+ GLOG() << l << " fl = " << cl[l] << std::endl;
+ }
+ calculate_Fl(4.5, 0.5, 5, cl);
+ for (int l = 0; l <= 4; ++l) {
+ GLOG() << l << " fl = " << cl[l] << std::endl;
+ }
+ calculate_Fl(10.5, 0.5, 5, cl);
+ for (int l = 0; l <= 4; ++l) {
+ GLOG() << l << " fl = " << cl[l] << std::endl;
+ }
+
+ GLOG() << "Solid harmonics Rlm ..." << std::endl;
+ std::vector<std::complex<double>> rlm;
+ 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(d, r, L, rlm);
+ 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] << std::endl;
+ }
+ }
+ GLOG() << "r = 0 ..." << std::endl;
+ r = 0.;
+ calculate_solidharm_Rlm(d, r, L, rlm);
+ 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] << std::endl;
+ }
+ }
- GLOG() << "Legendre ..." << std::endl;
+ /*GLOG() << "Legendre ..." << std::endl;
std::vector<double> plm;
calculate_legendre_plm(L, 0.3, plm);
for (int l = 0; l <= L; ++l) {
@@ -62,7 +306,6 @@ void FTSpectrum::compute() {
GLOG() << l << " " << m << " " << plm[lm] << std::endl;
}
}
-
GLOG() << "Factorial ..." << std::endl;
for (int n = 0; n < 16; ++n) {
GLOG() << n << "! = " << factorial(n) << std::endl;
@@ -87,6 +330,9 @@ void FTSpectrum::compute() {
}
}
+
+ */
+
return;
}
diff --git a/src/soap/functions.cpp b/src/soap/functions.cpp
index 0291b8b..d80c846 100644
--- a/src/soap/functions.cpp
+++ b/src/soap/functions.cpp
@@ -270,6 +270,70 @@ long int factorial(int x) {
}
}
+const int FACTORIAL2_CACHE_SIZE = 16;
+const long int FACTORIAL2_CACHE[] = {
+ 1,
+ 1, 2, 3,
+ 8, 15, 48,
+ 105, 384, 945,
+ 3840, 10395, 46080,
+ 135135, 645120, 2027025 };
+
+long int factorial2(int x) {
+ if (x < 0) {
+ assert(x == -1);
+ return 1;
+ }
+ else if (x < FACTORIAL2_CACHE_SIZE) {
+ return FACTORIAL2_CACHE[x];
+ }
+ else {
+ assert(false && "Computing factorial from scratch - not desirable");
+ long int s = 1;
+ for (int n = x; n > 0; n -= 2) {
+ s *= n;
+ }
+ return s;
+ }
+}
+
+void calculate_solidharm_Rlm(
+ vec d,
+ double r,
+ int L,
+ std::vector<std::complex<double>> &rlm) {
+ // Initialise
+ rlm.clear();
+ rlm.resize((L+1)*(L+1), 0.0);
+ // Treat r=0
+ if (r < 1e-10) { // TODO Define SPACE-QUANTUM
+ //throw soap::base::APIError("Rlm(r=0) disabled by design: Handle r=0-case manually.");
+ rlm[0] = 1.;
+ return;
+ }
+ // Proceed with r != 0: Compute Associated Legendre Polynomials
+ std::vector<double> plm;
+ double theta = acos(d.getZ());
+ double 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);
+ // Add radial component, phase factors, normalization
+ 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]);
+ }
+ }
+ return;
+}
+
void calculate_solidharm_rlm_ilm(
vec d,
double r,
diff --git a/src/soap/functions.hpp b/src/soap/functions.hpp
index 8f75bf4..541b868 100644
--- a/src/soap/functions.hpp
+++ b/src/soap/functions.hpp
@@ -72,6 +72,16 @@ extern const int FACTORIAL_CACHE_SIZE;
extern const long int FACTORIAL_CACHE[];
long int factorial(int n);
+extern const int FACTORIAL2_CACHE_SIZE;
+extern const long int FACTORIAL2_CACHE[];
+long int factorial2(int n);
+
+void calculate_solidharm_Rlm(
+ vec d,
+ double r,
+ int L,
+ std::vector<std::complex<double>> &rlm);
+
void calculate_solidharm_rlm_ilm(
vec d,
double r,
--
GitLab