fieldtensor.cpp 43.2 KB
Newer Older
Carl Poelking's avatar
Carl Poelking committed
1
#include "soap/fieldtensor.hpp"
Carl Poelking's avatar
Carl Poelking committed
2
#include "soap/linalg/numpy.hpp"
Carl Poelking's avatar
Carl Poelking committed
3
#include "soap/linalg/wigner.hpp"
Carl Poelking's avatar
Carl Poelking committed
4
5
6
7
8
9
10
11
12
13
14
15
#include <boost/math/special_functions/legendre.hpp>

namespace soap {

namespace ub = boost::numeric::ublas;

// ============================
// Spectrum::FTSpectrum::Atomic
// ============================

const std::string AtomicSpectrumFT::_numpy_t = "float64";

Carl Poelking's avatar
Carl Poelking committed
16
AtomicSpectrumFT::AtomicSpectrumFT(Particle *center, int K, int L, Options *options)
Carl Poelking's avatar
Carl Poelking committed
17
    : _center(center), _K(K), _L(L) {
Carl Poelking's avatar
Carl Poelking committed
18
    this->_s = center->getTypeId()-1; // TODO No longer used, using string labels as map keys instead
Carl Poelking's avatar
Carl Poelking committed
19
    this->_type = center->getType();
Carl Poelking's avatar
Carl Poelking committed
20
21
22
23
24
    if (GLOG.verbose())
        GLOG() << "Created FT particle: " << _type
            << " " << center->getId()
            << " @ " << _center->getPos()
            << std::endl;
Carl Poelking's avatar
Carl Poelking committed
25
    //assert(_s >= 0 && "Type-IDs should start from 1");
Carl Poelking's avatar
Carl Poelking committed
26
27
28
    // Body-order storage
    assert(K >= 0);
    _body_map.resize(K+1);
29
    // Linear field responses
Carl Poelking's avatar
Carl Poelking committed
30
31
32
    double p = options->get<double>("fieldtensor.alpha.p");
    double q = options->get<double>("fieldtensor.alpha.q");
    double r = options->get<double>("fieldtensor.alpha.r");
33
34
35
    _alpha.clear();
    for (int l = 0; l <= _L; ++l) {
        //double alpha_l = factorial(l)*pow(2*l+1,0.0)*pow(_center->getSigma(), 2*l+1); // HACK
Carl Poelking's avatar
Carl Poelking committed
36
37
38
39
40
41
42
43
        //
        //double alpha_l = pow(2*l+1,1.0)*pow(_center->getSigma(), 2*l+1);
        //if (l == 0) alpha_l *= 0.5;
        //
        //double alpha_l = 1.;
        //
        double alpha_l = pow(2*l+1,p)*pow(_center->getSigma(), q*(2*l+1));
        if (l == 0) alpha_l *= r;
44
45
        _alpha.push_back(alpha_l);
    }
Carl Poelking's avatar
Carl Poelking committed
46
47
48
}

AtomicSpectrumFT::~AtomicSpectrumFT() {
Carl Poelking's avatar
Carl Poelking committed
49
50
51
    #ifdef DBG
        GLOG() << "[~] Destruct " << this->getCenter()->getId() << std::endl;
    #endif
Carl Poelking's avatar
Carl Poelking committed
52
    // Deallocate body-order field terms
Carl Poelking's avatar
Carl Poelking committed
53
54
    for (int k=0; k <= _K; ++k) {
        field_map_t &fm = _body_map[k];
Carl Poelking's avatar
Carl Poelking committed
55
56
57
        #ifdef DBG
            GLOG() << "[~]   Deallocating k=" << k<< std::endl;
        #endif
Carl Poelking's avatar
Carl Poelking committed
58
59
        // Deallocate field terms for each type channel
        for (auto it=fm.begin(); it != fm.end(); ++it) {
Carl Poelking's avatar
Carl Poelking committed
60
61
62
            #ifdef DBG
                GLOG() << "[~]     Deallocating type=" <<  it->first << std::endl;
            #endif
Carl Poelking's avatar
Carl Poelking committed
63
64
65
            delete it->second;
        }
    }
Carl Poelking's avatar
Carl Poelking committed
66
67
    // Deallocate contraction coefficients
    for (auto it = _coeff_map.begin(); it != _coeff_map.end(); ++it) {
Carl Poelking's avatar
Carl Poelking committed
68
69
70
71
72
        #ifdef DBG
            GLOG() << "[~]   Deallocating s1:s2 = "
                <<  it->first.first << ":" << it->first.second
                << std::endl;
        #endif
Carl Poelking's avatar
Carl Poelking committed
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
        delete it->second;
    }
}

boost::python::list AtomicSpectrumFT::getTypes() {
    boost::python::list types;
    // The highest-body-order term should have the complete set
    for (auto it = _body_map[_K].begin(); it != _body_map[_K].end(); ++it) {
        std::string type = it->first;
        types.append(type);
    }
    return types;
}

boost::python::object AtomicSpectrumFT::getCoefficientsNumpy(std::string s1, std::string s2) {
    soap::linalg::numpy_converter npc(_numpy_t.c_str());
    channel_t channel(s1, s2);
    return npc.ublas_to_numpy< dtype_t >(*_coeff_map[channel]);
Carl Poelking's avatar
Carl Poelking committed
91
92
93
94
95
96
97
98
99
100
101
102
}

void AtomicSpectrumFT::addField(int k, std::string type, field_t &flm) {
    field_t &f = *(this->getCreateField(k, type, flm.size1(), flm.size2()));
    f = f + flm;
    return;
}

AtomicSpectrumFT::field_t *AtomicSpectrumFT::getCreateField(int k, std::string type, int s1, int s2) {
    assert(k <= _K);
    auto it = _body_map[k].find(type);
    if (it == _body_map[k].end()) {
Carl Poelking's avatar
Carl Poelking committed
103
104
105
106
107
108
109
        #ifdef DBG
            GLOG() << "[AtomicSpectrumFT::getCreateField]"
                << " Allocate k " << k
                << " type '" << type
                << "' [" << s1 << "x" << s2 << "]"
                << std::endl;
        #endif
Carl Poelking's avatar
Carl Poelking committed
110
111
112
113
114
115
116
117
118
119
120
121
122
        _body_map[k][type] = new field_t(
            field_zero_t(s1, s2)
        );
        it = _body_map[k].find(type);
    }
    return it->second;
}

AtomicSpectrumFT::coeff_t *AtomicSpectrumFT::getCreateContraction(channel_t &channel, int size1, int size2) {
    // size1 -> trail dimension lambda_total
    // size2 -> L+1
    auto it = _coeff_map.find(channel);
    if (it == _coeff_map.end()) {
Carl Poelking's avatar
Carl Poelking committed
123
124
125
126
127
128
        #ifdef DBG
            GLOG() << "[AtomicSpectrumFT::getCreateContraction]"
                << " Allocate " << channel.first << ":" << channel.second
                << " " << size1 << "x" << size2
                << std::endl;
        #endif
Carl Poelking's avatar
Carl Poelking committed
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
        _coeff_map[channel] = new coeff_t(
            coeff_zero_t(size1, size2)
        );
        it = _coeff_map.find(channel);
    }
    return it->second;
}

AtomicSpectrumFT::field_t *AtomicSpectrumFT::getField(int k, std::string type) {
    assert(k <= _K);
    auto it = _body_map[k].find(type);
    if (it == _body_map[k].end()) {
        return NULL;
    }
    else {
        return it->second;
    }
}

Carl Poelking's avatar
Carl Poelking committed
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
void AtomicSpectrumFT::contractDeep() {
    /*
    for (int l1=0; l1<=_L; ++l1) {
    for (int l2=0; l2<=_L; ++l2) {
    for (int l3=0; l3<=_L; ++l3) {
       int m1 = -l1;
       int m2 = -l2;
       int m3 = -m1-m2;
       double w = WignerSymbols::wigner3j(l1, l2, l3, m1, m2, m3);
       GLOG() << l1 << ":" << l2 << ":" << l3 << " " << m1 << ":" << m2 << ":" << m3 << " " << w << std::endl;
    }}}

    for (int l1=0; l1<=_L; ++l1) {
    for (int l2=0; l2<=_L; ++l2) {
    for (int l3=0; l3<=0; ++l3) {
        for (int m1=-l1; m1<=l1; ++m1) {
            int m2 = -m1;
            int m3 = 0;
            double w = WignerSymbols::wigner3j(l1, l2, l3, m1, m2, m3);
            GLOG() << l1 << ":" << l2 << ":" << l3 << " " << m1 << ":" << m2 << ":" << m3 << " " << w << std::endl;
        }
    }}}
    */


    // Wig. 3j(2)     _L =  0  1  2  3   4   5   6   7    8    9    10   11   12   13   14   15
    std::vector<int> lll = {0, 1, 5, 13, 27, 48, 78, 118, 170, 235, 315, 411, 525, 658, 812, 988 };
    int k = 1;
    int length = (_L+1) + (_L+1)*(_L+1) + lll[_L];
    field_map_t &field_map = _body_map[k];

    // Contract channels
    for (auto it1 = field_map.begin(); it1 != field_map.end(); ++it1) {
        std::string s1 = it1->first; // <- Channel 1 (e.g., 'C')
        field_t &f1 = *(it1->second);
        for (auto it2 = field_map.begin(); it2 != field_map.end(); ++it2) {
            std::string s2 = it2->first; // <- Channel 2 (e.g., 'O')
            field_t &f2 = *(it2->second);
            GLOG() << " " << s1 << ":" << s2 << std::flush;

            // Allocate channel if necessary
            channel_t channel(s1, s2);
            coeff_t &coeffs = *(this->getCreateContraction(channel, length, 1)); // TODO Will not work with k-dependent L-cutoff

            // a,l <> a',l contractions
            int offset = 0;
            for (int l = 0; l <= _L; ++l) {
                double inv_alpha_l = 1./this->_alpha[l];
                cmplx_t phi_l_s1s2 = 0.0;
                for (int m = -l; m <= l; ++m) {
                    int lm = l*l+l+m;
                    phi_l_s1s2 += f1(0, lm)*std::conj(f2(0, lm));
                }
                int l1l2 = offset+l;
                coeffs(l1l2, 0) = pow(this->getCenter()->getSigma(), 2*l)*inv_alpha_l*phi_l_s1s2.real();
            } // l

            // al, a'l', a'l' contractions
            offset = _L+1;
            for (int l1 = 0; l1 <= _L; ++l1) {
                double inv_alpha_l1 = 1./this->_alpha[l1];
            for (int l2 = 0; l2 <= _L; ++l2) {
                double inv_alpha_l2 = 1./this->_alpha[l2];
                cmplx_t phi_l1l2l2_s1s2 = 0.0;
                double w122_sq = 0.0;
                for (int m1 = -l1; m1 <= l1; ++m1) {
                    int lm1 = l1*l1+l1+m1;
                for (int m2 = -l2; m2 <= l2; ++m2) {
                    int lm2 = l2*l2+l2+m2;
                for (int m3 = -l2; m3 <= l2; ++m3) {
                    int lm3 = l2*l2+l2+m3;
                    // TODO Either (al1 bl2 bl2) or (al2 bl1 bl2)
                    double w122 = WignerSymbols::wigner3j(l2, l1, l2, m2, m1, m3);
                    w122_sq += w122*w122;
                    //GLOG() << l1 << ":" << l2 << ":" << m1 << ":" << m2 << "  " << w122 << std::endl;
                    phi_l1l2l2_s1s2 += w122 * f1(0,lm2) * f2(0,lm1) * f2(0,lm3);
                }}}
                GLOG() << l1 << ":" << l2 << " " << w122_sq << std::endl;
                int l1l2l2 = offset + l1*(_L+1) + l2;
                coeffs(l1l2l2, 0) = 
                        sqrt(inv_alpha_l1)*sqrt(inv_alpha_l2)*sqrt(inv_alpha_l2)
                      * pow(this->getCenter()->getSigma(), l1+0.5)
                      * pow(this->getCenter()->getSigma(), 2*l2) *
                        pow(-1, l2-l1+l2)
                      * phi_l1l2l2_s1s2.real();
            }}

        } // Channel type 2
    } // Channel type 1
    GLOG() << std::endl;

    return;
}

Carl Poelking's avatar
Carl Poelking committed
242
243
244
245
246
void AtomicSpectrumFT::contract() {
    int k = 0;
    int Lambda_total = (1 - pow(_L+1, _K))/(1 - (_L+1)); // from geometric series
    int LambdaLambda_total = (1 - pow(_L+1, 2*_K))/(1 - pow(_L+1, 2));

Carl Poelking's avatar
Carl Poelking committed
247
248
249
250
251
252
253
    GLOG() << "[contract] Centre: " << this->getCenter()->getId() << ":" << this->getType() << std::flush;
    #ifdef DBG
        GLOG() << std::endl;
        GLOG() << "[contract] Trail length total: " << Lambda_total << std::endl;
        GLOG() << "[contract] Trail length total - contracted: " << LambdaLambda_total << std::endl;
    #endif
    for (k=1; k<=_K; ++k) { // TODO Mix different k-channels? That would increase descriptor length if affordable/required.
Carl Poelking's avatar
Carl Poelking committed
254
255
256
257
258
259

        int Lambda_off_k = (1 - pow(_L+1, k-1))/(1 - (_L+1));
        int LambdaLambda_off_k = (1 - pow(_L+1, 2*(k-1)))/(1 - pow(_L+1, 2));
        int Lambda_k = pow(_L+1, k-1);
        int L_k = _L;

Carl Poelking's avatar
Carl Poelking committed
260
261
262
263
264
265
266
267
        #ifdef DBG
            GLOG()
                << "[contract] L_k= " << _L
                << " Lambda_k= " << Lambda_k
                << " Lambda_off_k= " << Lambda_off_k
                << " LambdaLambda_off_k= " << LambdaLambda_off_k
                << std::endl;
        #endif
Carl Poelking's avatar
Carl Poelking committed
268
269

        field_map_t &field_map = _body_map[k];
Carl Poelking's avatar
Carl Poelking committed
270
271
272
273
274
275
276
277
278
279
280
        #ifdef DBG
            for (auto it = field_map.begin(); it != field_map.end(); ++it) {
                std::string type = it->first;
                field_t &field = *(it->second);
                GLOG() << "[contract]  k= " << k
                    << " type= " << type
                    << " ||lambda||= " << field.size1()
                    << " ||L||= " << field.size2()
                    << std::endl;
            }
        #endif
Carl Poelking's avatar
Carl Poelking committed
281
282
        // Contract channels
        for (auto it1 = field_map.begin(); it1 != field_map.end(); ++it1) {
Carl Poelking's avatar
Carl Poelking committed
283
            std::string s1 = it1->first; // <- Channel 1 (e.g., 'C')
Carl Poelking's avatar
Carl Poelking committed
284
285
            field_t &f1 = *(it1->second);
            for (auto it2 = field_map.begin(); it2 != field_map.end(); ++it2) {
Carl Poelking's avatar
Carl Poelking committed
286
                std::string s2 = it2->first; // <- Channel 2 (e.g., 'O')
Carl Poelking's avatar
Carl Poelking committed
287
288
289
290
291
292
293
294
295
296
297
298
299
                field_t &f2 = *(it2->second);
                GLOG() << " " << s1 << ":" << s2 << std::flush;
                // Size sanity checks
                int Lambda1 = f1.size1();
                int Lambda2 = f2.size1();
                int L1 = f1.size2();
                int L2 = f2.size2();
                assert(L1 == (_L+1)*(_L+1));
                assert(L2 == (_L+1)*(_L+1));
                assert(Lambda1 == Lambda_k);
                assert(Lambda2 == Lambda_k);
                // Allocate channel if necessary
                channel_t channel(s1, s2);
Carl Poelking's avatar
Carl Poelking committed
300
301
302
303
304
305
                coeff_t &coeffs = *(this->getCreateContraction(channel, LambdaLambda_total, _L+1)); // TODO Will not work with k-dependent L-cutoff
                #ifdef DBG
                    GLOG() << "Coefficient matrix for this channel: "
                        << coeffs.size1() << "x" << coeffs.size2()
                        << std::endl;
                #endif
Carl Poelking's avatar
Carl Poelking committed
306
307
308
309
                for (int lambda1 = 0; lambda1 < Lambda1; ++lambda1) {
                    for (int lambda2 = 0; lambda2 < Lambda2; ++lambda2) {
                        for (int l = 0; l <= L_k; ++l) {
                            // Polarizability
310
311
312
313
                            //double inv_alpha = 1./(2*l+1)*pow(this->getCenter()->getSigma(), -2*l-1); // TODO Initialise in constructor HACK 1./(2*l+1)
                            //double inv_alpha = 1./(2*l+1)*pow(this->getCenter()->getSigma(), -l);
                            //double inv_alpha = 1./(l+0.1)*pow(this->getCenter()->getSigma(), -l);
                            double inv_alpha = 1./this->_alpha[l];
Carl Poelking's avatar
Carl Poelking committed
314
315
316
                            cmplx_t phi_l_s1s2_l1l2 = 0.0;
                            for (int m = -l; m <= l; ++m) {
                                int lm = l*l+l+m;
Carl Poelking's avatar
Carl Poelking committed
317
                                phi_l_s1s2_l1l2 += f1(lambda1, lm)*std::conj(f2(lambda2, lm));
Carl Poelking's avatar
Carl Poelking committed
318
319
                            }
                            int l1l2 = LambdaLambda_off_k+lambda1*Lambda1 + lambda2;
Carl Poelking's avatar
Carl Poelking committed
320
321
322
323
324
325
326
                            #ifdef DBG
                                GLOG() << " Store "
                                    << lambda1 << ":" << lambda2 << ":" << l
                                    << " @ " << l1l2 << ":" << l
                                    << " = " << phi_l_s1s2_l1l2
                                    <<  std::endl;
                            #endif
Carl Poelking's avatar
Carl Poelking committed
327
328
329
330
                            //coeffs(l1l2, l) = inv_alpha*phi_l_s1s2_l1l2.real(); // TODO HACK
                            coeffs(l1l2, l) = pow(this->getCenter()->getSigma(), 2*l)*inv_alpha*phi_l_s1s2_l1l2.real(); // TODO HACK
                            //coeffs(l1l2, l) = phi_l_s1s2_l1l2.real(); // TODO HACK

Carl Poelking's avatar
Carl Poelking committed
331
332
333
334
335
                        } // l
                    } // lambda 2
                } // lambda 1
            } // Channel type 2
        } // Channel type 1
Carl Poelking's avatar
Carl Poelking committed
336
    }
Carl Poelking's avatar
Carl Poelking committed
337
    GLOG() << std::endl;
Carl Poelking's avatar
Carl Poelking committed
338
    return;
Carl Poelking's avatar
Carl Poelking committed
339
340
341
342
}

void AtomicSpectrumFT::registerPython() {
    using namespace boost::python;
Carl Poelking's avatar
Carl Poelking committed
343
    class_<AtomicSpectrumFT, AtomicSpectrumFT*>("AtomicSpectrumFT", init<Particle*, int, int, Options*>())
Carl Poelking's avatar
Carl Poelking committed
344
345
        .def("getTypes", &AtomicSpectrumFT::getTypes)
        .def("getPower", &AtomicSpectrumFT::getCoefficientsNumpy)
Carl Poelking's avatar
Carl Poelking committed
346
347
348
349
350
351
352
        .def("getCenter", &AtomicSpectrumFT::getCenter, return_value_policy<reference_existing_object>());
}

// ===================
// FieldTensorSpectrum
// ===================

Carl Poelking's avatar
Carl Poelking committed
353
354
355
356
357
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) {
Carl Poelking's avatar
Carl Poelking committed
358
        for (int m = -l; m <= l; ++m) {
Carl Poelking's avatar
Carl Poelking committed
359
            out[l*l+l+m] = sqrt(factorial(l+m)*factorial(l-m));
Carl Poelking's avatar
Carl Poelking committed
360
361
362
            #ifdef DBG // TODO Outsource to tests
                GLOG() << "[A]" << l << "," << m << " : " << out[l*l+l+m] << std::endl;
            #endif
Carl Poelking's avatar
Carl Poelking committed
363
364
365
366
367
368
369
370
371
372
        }
    }
    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) {
Carl Poelking's avatar
Carl Poelking committed
373
        for (int m = -l; m <= l; ++m) {
Carl Poelking's avatar
Carl Poelking committed
374
            out[l*l+l+m] = sqrt(factorial(l+m)*factorial(l-m))/factorial2(2*l-1);
Carl Poelking's avatar
Carl Poelking committed
375
376
377
            #ifdef DBG // TODO Outsource to tests
                GLOG() << "[B]" << l << "," << m << " : " << out[l*l+l+m] << std::endl;
            #endif
Carl Poelking's avatar
Carl Poelking committed
378
379
380
381
382
383
384
        }
    }
    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++
Carl Poelking's avatar
Carl Poelking committed
385
    // (which does not however mention the recursive form used here)
Carl Poelking's avatar
Carl Poelking committed
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
    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;
}

Carl Poelking's avatar
Carl Poelking committed
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
void Tlmlm::computeTlmlm(vec d12, double r12, double a12, int L1, int L2, coeff_t &T) {
    // TODO This can be sped up by using Tl1m1l2m2 = (-1)^(l1+l2) Tl2m2l1m1
    // TODO Properties to test:
    //      Tl1l2 ~ (l-1)! * 1/r^{l1+l2+1} for large r
    //      Tl1l2 = (-1)^{l1+l2} Tl2l1 and T(-r) = (-1)^{l1+l2} T(r)
    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);

Carl Poelking's avatar
Carl Poelking committed
436
    /*
Carl Poelking's avatar
Carl Poelking committed
437
438
439
440
441
442
443
    for (int l1 = 0; l1 <= Lg; ++l1) {
        std::complex<double> sum_l = 0.0;
        for (int m1 = -l1; m1 <= l1; ++m1) {
            int lm1 = l1*l1+l1+m1;
            GLOG() << "[Rlm] " << " m1 " << l1 << " : " << Rlm[lm1] << std::endl;
            sum_l += Rlm[lm1]*std::conj(Rlm[lm1]);
        }
Carl Poelking's avatar
Carl Poelking committed
444
        GLOG() << "[Rlm] => " << r12 << " l " << l1 << " : " << sum_l << std::endl;
Carl Poelking's avatar
Carl Poelking committed
445
446
447
448
    }
    for (int l1 = 0; l1 <= L1; ++l1) {
    for (int l2 = 0; l2 <= L2; ++l2) {
        std::complex<double> sum_l1_l2 = 0.0;
Carl Poelking's avatar
Carl Poelking committed
449
450
451
452
        for (int m1 = -l1; m1 <= l1; ++m1) {
            int lm1 = l1*l1+l1+m1;
        for (int m2 = -l2; m2 <= l2; ++m2) {
            int lm2 = l2*l2+l2+m2;
Carl Poelking's avatar
Carl Poelking committed
453
454
            sum_l1_l2 += T(lm1,lm2)*std::conj(T(lm1,lm2));
        }}
Carl Poelking's avatar
Carl Poelking committed
455
        GLOG() << "[Tlmlm] => " << l1 << ":" << l2 << " sum " << sum_l1_l2 << std::endl;
Carl Poelking's avatar
Carl Poelking committed
456
    }}
Carl Poelking's avatar
Carl Poelking committed
457
    */
Carl Poelking's avatar
Carl Poelking committed
458

Carl Poelking's avatar
Carl Poelking committed
459
460
461
462
463
464
    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;
Carl Poelking's avatar
Carl Poelking committed
465

Carl Poelking's avatar
Carl Poelking committed
466
467
        std::complex<double> tlm1lm2 = 0.0;
        int lg = (l1 > l2) ? l1 : l2;
Carl Poelking's avatar
Carl Poelking committed
468

Carl Poelking's avatar
Carl Poelking committed
469
470
        for (int l = 0; l <= lg; ++l) {
        for (int m = -l; m <= l; ++m) {
Carl Poelking's avatar
Carl Poelking committed
471

Carl Poelking's avatar
Carl Poelking committed
472
473
474
475
476
477
478
479
480
481
482
            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 ];

            //tlm1lm2 +=
            //      pow(-1, l2+m)*1./(Alm[lm]*Blm[lm])
Carl Poelking's avatar
Carl Poelking committed
483
484
            //    * std::conj( Rlm[ (l2-l)*(l2-l) + (l2-l) + m2+m ] ) / Alm[(l2-l)*(l2-l) + (l2-l) + m2+m ]
            //    * std::conj( Rlm[ (l1-l)*(l1-l) + (l1-l) + m1-m ] ) / Alm[(l1-l)*(l1-l) + (l1-l) + m1-m ]
Carl Poelking's avatar
Carl Poelking committed
485
486
487
488
489
490
491
            //    * Fl[ l1+l2-l ];

        }} // l m

        tlm1lm2 *=
              Alm[lm1]*Alm[lm1] * Alm[lm2]*Alm[lm2]
           / (factorial2(2*l1-1)*factorial2(2*l2-1));
Carl Poelking's avatar
Carl Poelking committed
492
493
        tlm1lm2 /=               // HACK Is that the fix? Yes. TODO Adjust scaling
            Blm[lm1]*Blm[lm2];   // HACK Is that the fix? Yes. TODO Adjust scaling
Carl Poelking's avatar
Carl Poelking committed
494
        T(lm1,lm2) = tlm1lm2;
Carl Poelking's avatar
Carl Poelking committed
495
    }}}} // l1 m1 l2 m2
Carl Poelking's avatar
Carl Poelking committed
496
497
498
    return;
}

Carl Poelking's avatar
Carl Poelking committed
499
500
FTSpectrum::FTSpectrum(Structure &structure, Options &options)
    : _structure(&structure), _options(&options) {
Carl Poelking's avatar
Carl Poelking committed
501
502
503
    #ifdef DBG
        GLOG() << "FTSPECTRUM: DEBUG MODE ENABLED" << std::endl;
    #endif
Carl Poelking's avatar
Carl Poelking committed
504
505
    _cutoff = CutoffFunctionOutlet().create(_options->get<std::string>("radialcutoff.type"));
	_cutoff->configure(*_options);
Carl Poelking's avatar
Carl Poelking committed
506
507
    _L = _options->get<int>("fieldtensor.L");
    _K = _options->get<int>("fieldtensor.K");
Carl Poelking's avatar
Carl Poelking committed
508
509
510
511
    return;
}

FTSpectrum::~FTSpectrum() {
Carl Poelking's avatar
Carl Poelking committed
512
513
514
515
516
517
518
    for (auto it = _atomic_array.begin(); it != _atomic_array.end(); ++it) {
        delete *it;
    }
    return;
}

void InteractUpdateSourceTarget(
Carl Poelking's avatar
Carl Poelking committed
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
        AtomicSpectrumFT *a, // <- source particle
        AtomicSpectrumFT *b, // <- target particle
        Tlmlm::coeff_t &T_ab, // <- field tensors
        double sigma1,      // <- atomic width source
        double sigma2,      // <- atomic width target
        double w12,         // <- pair weight
        double w2,          // <- target weight
        int k,              // <- Body-order cutoff
        int L_in,           // <- Incoming moments: cutoff
        int L_out,          // <- Outgoing moments: cutoff
        int Lambda_in,      // <- Trace memory length in
        int Lambda_out,     // <- Trace memory length out
        int LM_in,          // <- (no longer used)
        int LM_out,         // <- Allocation size out
        bool apply_parity) { // <- apply parity to field tensors

    #ifdef DBG
        GLOG() << "[update] " << "w " << w12 << " " << w2 << std::endl;
        GLOG() << "[update] " << "k " << k << " L-in " << L_in << " L-out " << L_out << std::endl;
        GLOG() << "[update] " << "Lambda-in " << Lambda_in << " Lambda-out " << Lambda_out << std::endl;
        GLOG() << "[update] " << "LM-in " << LM_in << " LM-out " << LM_out << std::endl;
        GLOG() << "[update] " << "Parity " << apply_parity << std::endl;
        for (int l1=0; l1 <= L_out; ++l1) {
Carl Poelking's avatar
Carl Poelking committed
542
543
544
        for (int m1=-l1; m1 <= l1; ++m1) {
            int lm1 = l1*l1+l1+m1;
            for (int l2=0; l2 <= L_out; ++l2) {
Carl Poelking's avatar
Carl Poelking committed
545
546
547
548
549
550
551
552
553
            for (int m2=-l2; m2 <= l2; ++m2) {
                int lm2 = l2*l2+l2+m2;
                GLOG() << l1 << "," << m1
                    << " " << l2 << "," << m2
                    << " " << T_ab(lm1,lm2)
                    << std::endl;
            }}
        }}
    #endif
Carl Poelking's avatar
Carl Poelking committed
554
555
556
557
558
559
560
561
562
563

    // Retrieve field maps for iteration k-1
    AtomicSpectrumFT::field_map_t &fmap_a = a->getFieldMap(k-1);
    AtomicSpectrumFT::field_map_t &fmap_b = b->getFieldMap(k-1);
    // TARGET: MOMENTS ON 'b'
    // Loop over field channels/types: a
    for (auto it_a = fmap_a.begin(); it_a != fmap_a.end(); ++it_a) {
        std::string type = it_a->first;
        AtomicSpectrumFT::field_t &f_a = *(it_a->second);
        AtomicSpectrumFT::field_t &f_b = *(b->getCreateField(k, type, Lambda_out, LM_out));
Carl Poelking's avatar
Carl Poelking committed
564
565
566
567
568
569
570
571
        #ifdef DBG
            GLOG() << "[increment]   a: source: type: "
                << type << " size: " << f_a.size1()
                << " x " << f_a.size2()
                << std::endl;
            GLOG() << "[increment]   b: target: size: "
                << f_b.size1() << " x " << f_b.size2() << std::endl;
        #endif
Carl Poelking's avatar
Carl Poelking committed
572
573
574
        // lm-out
        for (int l1=0; l1<=L_out; ++l1) {
            // Polarizability
575
576
577
578
            //double alpha = (2*l1+1)*pow(sigma2, 2*l1+1); // TODO HACK: 2*l1+1
            //double alpha = (2*l1+1)*pow(sigma2, l1);
            //double alpha = (l1+0.1)*pow(sigma2, l1);
            double alpha = b->_alpha[l1];
Carl Poelking's avatar
Carl Poelking committed
579
580
581
582
583
584
585
586
587
588
            for (int m1=-l1; m1<=l1; ++m1) {
                int lm_out = l1*l1+l1+m1;
                // lambda-in
                for (int lambda_in=0; lambda_in<Lambda_in; ++lambda_in) {
                    // lm-in
                    for (int l2=0; l2<=L_in; ++l2) {
                        double parity = (apply_parity) ? pow(-1,l1+l2) : 1;
                        int lambda_out = lambda_in*Lambda_in+l2;
                        AtomicSpectrumFT::cmplx_t q_lambda_lm_out = 0.0;
                        for (int m2=-l2; m2<=l2; ++m2) {
Carl Poelking's avatar
Carl Poelking committed
589
590
591
592
593
594
595
596
                            #ifdef DBG
                                GLOG() << "l_out " << l1
                                    << " m_out " << m1
                                    << " lambda_in " << lambda_in
                                    << " l_in " << l2
                                    << " m_in " << m2
                                    << std::endl;
                            #endif
Carl Poelking's avatar
Carl Poelking committed
597
                            int lm_in = l2*l2+l2+m2;
Carl Poelking's avatar
Carl Poelking committed
598
599
600
                            //q_lambda_lm_out += T_ab(lm_out, lm_in)*f_a(lambda_in, lm_in); // BACK-UP
                            //q_lambda_lm_out += pow(sigma1, l1+0.5)*pow(sigma2, l2+0.5)*T_ab(lm_out, lm_in)*f_a(lambda_in, lm_in); // TODO HACK
                            q_lambda_lm_out += pow(sigma1, -l1)*pow(sigma2, l2)*T_ab(lm_out, lm_in)*f_a(lambda_in, lm_in); // TODO HACK
Carl Poelking's avatar
Carl Poelking committed
601
602
603
604
                        }
                        // Apply weight and parity
                        q_lambda_lm_out *= alpha*w12*w2*parity;
                        // Add to fields on b
Carl Poelking's avatar
Carl Poelking committed
605
606
607
608
609
610
                        #ifdef DBG
                            GLOG() << "    => (" << lambda_out << "," << lm_out << ")"
                                << " : " << q_lambda_lm_out
                                << std::endl;
                        #endif
                        f_b(lambda_out, lm_out) += std::conj(q_lambda_lm_out); // FIXED Added conj
Carl Poelking's avatar
Carl Poelking committed
611
612
613
614
615
                    } // End loop over l_in
                } // End loop over lambda_in
            } // End loop over m_out
        } // End loop over l_out
    } // End loop over types on a
Carl Poelking's avatar
Carl Poelking committed
616
617
618
    return;
}

Carl Poelking's avatar
Carl Poelking committed
619
void FTSpectrum::computeFieldTensors(std::map<int, std::map<int, Tlmlm::coeff_t>> &i1_i2_T12) {
Carl Poelking's avatar
Carl Poelking committed
620
    GLOG() << "Computing field tensors ..." << std::endl;
Carl Poelking's avatar
Carl Poelking committed
621
622
    int K = _K;
    int L = _L;
Carl Poelking's avatar
Carl Poelking committed
623
    Tlmlm T12(L);
Carl Poelking's avatar
Carl Poelking committed
624
625
626
627
628
629
630
    for (auto it1 = beginAtomic(); it1 != endAtomic(); ++it1) {
        // Particle 1
        atomic_t *a = *it1;
        int id1 = a->getCenter()->getId();
        int s1 = a->getTypeIdx();
        vec r1 = a->getCenter()->getPos();
        double sigma1 = a->getCenter()->getSigma();
Carl Poelking's avatar
Carl Poelking committed
631
        // Initialise map
Carl Poelking's avatar
Carl Poelking committed
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
        i1_i2_T12[id1] = std::map<int, Tlmlm::coeff_t>();
        for (auto it2 = it1; it2 != endAtomic(); ++it2) {
            // Particle 2
            atomic_t *b = *it2;
            int id2 = b->getCenter()->getId();
            int s2 = b->getTypeIdx();
            vec r2 = b->getCenter()->getPos();
            double sigma2 = b->getCenter()->getSigma();
            // Find connection, apply weight function
            vec dr12 = _structure->connect(r1, r2);
            double r12 = soap::linalg::abs(dr12);
            vec d12 = dr12/r12;
            if (! _cutoff->isWithinCutoff(r12)) continue;
            double w12 = _cutoff->calculateWeight(r12);
            double a12 = 1./sqrt(2*(sigma1*sigma1 + sigma2*sigma2));
Carl Poelking's avatar
Carl Poelking committed
647
648
649
650
651
652
653
654
            #ifdef DBG
                GLOG() << "[field-tensor] " << a->getCenter()->getId()
                    << ":" << b->getCenter()->getId()
                    << " R=" << r12
                    << " w=" << w12
                    << " a=" << a12
                    << std::endl;
            #endif
Carl Poelking's avatar
Carl Poelking committed
655
656
657
658
659
660
            // Interact
            T12.computeTlmlm(d12, r12, a12, L, L, i1_i2_T12[id1][id2]);
        }
    }
    return;
}
Carl Poelking's avatar
Carl Poelking committed
661

Carl Poelking's avatar
Carl Poelking committed
662
void FTSpectrum::createAtomic() {
Carl Poelking's avatar
Carl Poelking committed
663
664
665
666
667
    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) {
Carl Poelking's avatar
Carl Poelking committed
668
        atomic_t *new_atomic = new atomic_t(*pit, _K, _L, _options);
Carl Poelking's avatar
Carl Poelking committed
669
670
        _atomic_array.push_back(new_atomic);
    }
Carl Poelking's avatar
Carl Poelking committed
671
}
Carl Poelking's avatar
Carl Poelking committed
672

Carl Poelking's avatar
Carl Poelking committed
673
674
void FTSpectrum::energySCF(int k, std::map<int, std::map<int, Tlmlm::coeff_t> > &i1_i2_T12) {
    std::complex<double> energy_total;
Carl Poelking's avatar
Carl Poelking committed
675
    for (auto it1 = beginAtomic(); it1 != endAtomic(); ++it1) {
Carl Poelking's avatar
Carl Poelking committed
676
        // Particle 1
Carl Poelking's avatar
Carl Poelking committed
677
        atomic_t *a = *it1;
Carl Poelking's avatar
Carl Poelking committed
678
        int id1 = a->getCenter()->getId();
Carl Poelking's avatar
Carl Poelking committed
679
        std::string s1 = a->getType();
Carl Poelking's avatar
Carl Poelking committed
680
681
        vec r1 = a->getCenter()->getPos();
        double sigma1 = a->getCenter()->getSigma();
Carl Poelking's avatar
Carl Poelking committed
682
        double w1 = a->getCenter()->getWeight();
Carl Poelking's avatar
Carl Poelking committed
683
684
        for (auto it2 = it1; it2 != endAtomic(); ++it2) {
            // Particle 2
Carl Poelking's avatar
Carl Poelking committed
685
            atomic_t *b = *it2;
Carl Poelking's avatar
Carl Poelking committed
686
            int id2 = b->getCenter()->getId();
Carl Poelking's avatar
Carl Poelking committed
687
            std::string s2 = b->getType();
Carl Poelking's avatar
Carl Poelking committed
688
689
            vec r2 = b->getCenter()->getPos();
            double sigma2 = b->getCenter()->getSigma();
Carl Poelking's avatar
Carl Poelking committed
690
            double w2 = b->getCenter()->getWeight();
Carl Poelking's avatar
Carl Poelking committed
691
692
693
694
695
696
697
            // Find connection, apply weight function
            vec dr12 = _structure->connect(r1, r2);
            double r12 = soap::linalg::abs(dr12);
            vec d12 = dr12/r12;
            if (! _cutoff->isWithinCutoff(r12)) continue;
            double w12 = _cutoff->calculateWeight(r12);
            double a12 = 1./sqrt(2*(sigma1*sigma1 + sigma2*sigma2));
Carl Poelking's avatar
Carl Poelking committed
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
            // Look up interaction tensors for pair (a,b)
            auto T_ab = i1_i2_T12[id1][id2];
            // Self-interaction a <> a
            if (a == b) {
                ;
            }
            // Cross-interaction a <> b
            else {
                std::complex<double> energy = 0.0;
                for (int l1=0; l1<=_L; ++l1) {
                    for (int m1=-l1; m1<=l1; ++m1) {
                        int lm1 = l1*l1+l1+m1;
                        for (int l2=0; l2<=_L; ++l2) {
                            for (int m2=-l2; m2<=l2; ++m2) {
                                int lm2 = l2*l2+l2+m2;
                                energy += a->_M(k, lm1)*T_ab(lm1, lm2)*b->_M(k, lm2);
                            }
                        }
                    }
                }
                energy_total += energy;
Carl Poelking's avatar
Carl Poelking committed
719
720
721
                GLOG() << "[energy] " << a->getCenter()->getId()
                    << ":" << b->getCenter()->getId() << " " << energy
                    << std::endl;
Carl Poelking's avatar
Carl Poelking committed
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
            }
        } // Particle b
    } // Particle a
    GLOG() << "[energy, total] " << "(k=" << k << ")" << energy_total << std::endl;
    return;
}

void FTSpectrum::polarize() {
    this->createAtomic();
    std::map<int, std::map<int, Tlmlm::coeff_t>> i1_i2_T12;
    this->computeFieldTensors(i1_i2_T12);

    // Initialise moments
    for (auto it = beginAtomic(); it != endAtomic(); ++it) {
        atomic_t *a = *it;
        a->_M = AtomicSpectrumFT::field_zero_t(_K+1,(_L+1)*(_L+1));
        a->_F = AtomicSpectrumFT::field_zero_t(_K+1,(_L+1)*(_L+1));
        a->_M(0,0) = 1.0;
Carl Poelking's avatar
Carl Poelking committed
740
    }
Carl Poelking's avatar
Carl Poelking committed
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
    this->energySCF(0, i1_i2_T12);

    for (int k = 1; k <= _K; ++k) {
        // UPDATE FIELDS
        for (auto it1 = beginAtomic(); it1 != endAtomic(); ++it1) {
            // Particle 1
            atomic_t *a = *it1;
            int id1 = a->getCenter()->getId();
            std::string s1 = a->getType();
            vec r1 = a->getCenter()->getPos();
            double sigma1 = a->getCenter()->getSigma();
            double w1 = a->getCenter()->getWeight();
            for (auto it2 = it1; it2 != endAtomic(); ++it2) {
                // Particle 2
                atomic_t *b = *it2;
                int id2 = b->getCenter()->getId();
                std::string s2 = b->getType();
                vec r2 = b->getCenter()->getPos();
                double sigma2 = b->getCenter()->getSigma();
                double w2 = b->getCenter()->getWeight();
                // Find connection, apply weight function
                vec dr12 = _structure->connect(r1, r2);
                double r12 = soap::linalg::abs(dr12);
                vec d12 = dr12/r12;
                if (! _cutoff->isWithinCutoff(r12)) continue;
                double w12 = _cutoff->calculateWeight(r12);
                double a12 = 1./sqrt(2*(sigma1*sigma1 + sigma2*sigma2));
Carl Poelking's avatar
Carl Poelking committed
768
769
770
771
772
773
                GLOG() << "[interact] " << a->getCenter()->getId()
                    << ":" << b->getCenter()->getId()
                    << " R=" << r12
                    << " w=" << w12
                    << " a=" << a12
                    << std::endl;
Carl Poelking's avatar
Carl Poelking committed
774
775
                // Look up interaction tensors for pair (a,b)
                auto T_ab = i1_i2_T12[id1][id2];
Carl Poelking's avatar
Carl Poelking committed
776
777
778
779
780
781
                for (int l1=0; l1 <= _L; ++l1) {
                    for (int m1=-l1; m1 <= l1; ++m1) {
                        int lm1 = l1*l1+l1+m1;
                        for (int l2=0; l2 <= _L; ++l2) {
                            for (int m2=-l2; m2 <= l2; ++m2) {
                                int lm2 = l2*l2+l2+m2;
Carl Poelking's avatar
Carl Poelking committed
782
783
784
785
786
787
788
                                #ifdef DBG
                                    GLOG() << "[FT] " << id1<<":"<<id2
                                        << " " << l1 << "," << m1
                                        << " " << l2 << "," << m2
                                        << " " << T_ab(lm1,lm2)
                                        << std::endl;
                                #endif
Carl Poelking's avatar
Carl Poelking committed
789
790
791
792
                            }
                        }
                    }
                }
Carl Poelking's avatar
Carl Poelking committed
793
794
795
796
797
798
799
800
801
802
803
804
805
806
                // Self-interaction a <> a
                if (a == b) {
                    ;
                }
                // Cross-interaction a <> b
                else {
                    for (int l1=0; l1<=_L; ++l1) {
                        for (int m1=-l1; m1<=l1; ++m1) {
                            int lm1 = l1*l1+l1+m1;
                            for (int l2=0; l2<=_L; ++l2) {
                                double parity = pow(-1,l1+l2);
                                for (int m2=-l2; m2<=l2; ++m2) {
                                    int lm2 = l2*l2+l2+m2;
                                    a->_F(k, lm1) += T_ab(lm1, lm2)*b->_M(k-1, lm2);
Carl Poelking's avatar
Carl Poelking committed
807
                                    b->_F(k, lm1) += T_ab(lm1, lm2)*a->_M(k-1, lm2)*parity;
Carl Poelking's avatar
Carl Poelking committed
808
809
810
811
812
813
814
815
816
817
                                    #ifdef DBG
                                        GLOG() << "@ " << id1 << " ~" << id2
                                            << " " << l1 << m1 << "|" << l2 << m2 << " : "
                                            << T_ab(lm1, lm2)*b->_M(k-1, lm2)
                                            << std::endl;
                                        GLOG() << "@ " << id2 << " ~" << id1
                                            << " " << l1 << m1 << "|" << l2 << m2 << " : "
                                            << T_ab(lm1, lm2)*a->_M(k-1, lm2)*parity
                                            << std::endl;
                                    #endif
Carl Poelking's avatar
Carl Poelking committed
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
                                }
                            }
                        }
                    }
                }
            } // Particle b
        } // Particle a
        // INDUCE
        for (auto it1 = beginAtomic(); it1 != endAtomic(); ++it1) {
            atomic_t *a = *it1;
            double sigma1 = a->getCenter()->getSigma();
            for (int l1=0; l1<=_L; ++l1) {
                std::complex<double> sum_l = 0.0;
                double alpha = pow(sigma1, 2*l1+1);
                for (int m1=-l1; m1<=l1; ++m1) {
                    int lm1 = l1*l1+l1+m1;
Carl Poelking's avatar
Carl Poelking committed
834
                    a->_M(k, lm1) = std::conj(alpha*a->_F(k, lm1)); // FIXED added conj
Carl Poelking's avatar
Carl Poelking committed
835
                    sum_l += a->_F(k, lm1)*std::conj(a->_F(k, lm1));
Carl Poelking's avatar
Carl Poelking committed
836
837
838
839
840
                    #ifdef DBG
                        GLOG() << "[induced] " << "l " << l1 << " m " << m1
                            << " lm " << lm1 <<  " : "  << a->_F(k, lm1)
                            << std::endl;
                    #endif
Carl Poelking's avatar
Carl Poelking committed
841
                }
Carl Poelking's avatar
Carl Poelking committed
842
843
844
                GLOG() << "[induced] @ " << a->getCenter()->getId()
                    << " l " << l1 << " : " << sum_l
                    << std::endl;
Carl Poelking's avatar
Carl Poelking committed
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
            }
        }
        // ENERGY
        this->energySCF(k, i1_i2_T12);
    } // Iterations k
    return;
}

void FTSpectrum::compute() {
    GLOG() << "Computing FTSpectrum ..." << std::endl;

    // CREATE ATOMIC SPECTRA
    this->createAtomic();

    // COMPUTE INTERACTION TENSORS
    std::map<int, std::map<int, Tlmlm::coeff_t>> i1_i2_T12;
    this->computeFieldTensors(i1_i2_T12);
Carl Poelking's avatar
Carl Poelking committed
862

Carl Poelking's avatar
Carl Poelking committed
863
    // INITIALISE FIELD MOMENTS (k=0)
Carl Poelking's avatar
Carl Poelking committed
864
865
    int K = _K;
    int L = _L;
Carl Poelking's avatar
Carl Poelking committed
866
    int k = 0;
Carl Poelking's avatar
Carl Poelking committed
867
    int L_k_in = 0;
Carl Poelking's avatar
Carl Poelking committed
868
869
870
871
872
873
874
875
876
    int L_k_out = 0;
    int Lambda_flat_k_in = 1;
    int Lambda_flat_k_out = 1; // || l'l''... ||
    GLOG() << "Initialise fields (k=0) ..." << std::endl;
    for (auto it = beginAtomic(); it != endAtomic(); ++it) {
        atomic_t *a = *it;
        std::string s = a->getType();
        AtomicSpectrumFT::field_t flm
            = AtomicSpectrumFT::field_zero_t(Lambda_flat_k_out, (L_k_in+1)*(L_k_in+1));
Carl Poelking's avatar
Carl Poelking committed
877
878
879
880
881
882
        #ifdef DBG
            GLOG() << "[init] Center:" << a->getCenter()->getId() << std::endl;
            GLOG() << "[init]   k = " << k << std::endl;
            GLOG() << "[init]   s = " << s << std::endl;
            GLOG() << "[init]   F : " << flm.size1() << " x " << flm.size2() << std::endl;
        #endif
Carl Poelking's avatar
Carl Poelking committed
883
884
885
886
887
888
889
        flm(0,0) = 1.0;
        a->addField(k, s, flm);
    }

    // PROPAGATE FIELD MOMENTS (k>0)
    for (int k = 1; k <= K; ++k) {

Carl Poelking's avatar
Carl Poelking committed
890
891
        GLOG() << "Starting iteration: k = " << k << std::endl;
        // Update out-going index ranges
Carl Poelking's avatar
Carl Poelking committed
892
893
894
895
        L_k_out = L;
        Lambda_flat_k_out = Lambda_flat_k_out*(L_k_in+1);
        int LM_in = (L_k_in+1)*(L_k_in+1);
        int LM_out = (L_k_out+1)*(L_k_out+1);
Carl Poelking's avatar
Carl Poelking committed
896
897
898
899
900
901
902
903
904
905
906
        #ifdef DBG
            GLOG() << "[iter] L-in " << L_k_in
                << " L-out " << L_k_out
                << " Lambda-in " << Lambda_flat_k_in
                << " Lambda-out " << Lambda_flat_k_out
                << std::endl;
            GLOG() << "[iter] Rank: " << L
                << " (LM-in = " << LM_in << ")"
                << " (LM-out = " << LM_out << ")"
                << std::endl;
        #endif
Carl Poelking's avatar
Carl Poelking committed
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
        for (auto it1 = beginAtomic(); it1 != endAtomic(); ++it1) {
            // Particle 1
            atomic_t *a = *it1;
            int id1 = a->getCenter()->getId();
            std::string s1 = a->getType();
            vec r1 = a->getCenter()->getPos();
            double sigma1 = a->getCenter()->getSigma();
            double w1 = a->getCenter()->getWeight();
            for (auto it2 = it1; it2 != endAtomic(); ++it2) {
                // Particle 2
                atomic_t *b = *it2;
                int id2 = b->getCenter()->getId();
                std::string s2 = b->getType();
                vec r2 = b->getCenter()->getPos();
                double sigma2 = b->getCenter()->getSigma();
                double w2 = b->getCenter()->getWeight();
                // Find connection, apply weight function
                vec dr12 = _structure->connect(r1, r2);
                double r12 = soap::linalg::abs(dr12);
                vec d12 = dr12/r12;
                if (! _cutoff->isWithinCutoff(r12)) continue;
                double w12 = _cutoff->calculateWeight(r12);
                double a12 = 1./sqrt(2*(sigma1*sigma1 + sigma2*sigma2));
Carl Poelking's avatar
Carl Poelking committed
930
931
932
933
934
935
936
                if (GLOG.verbose())
                    GLOG() << "[interact] " << a->getCenter()->getId()
                        << ":" << b->getCenter()->getId()
                        << " R=" << r12
                        << " w=" << w12
                        << " a=" << a12
                        << std::endl;
Carl Poelking's avatar
Carl Poelking committed
937
938
939
940
                // Look up interaction tensors for pair (a,b)
                auto T_ab = i1_i2_T12[id1][id2];
                // Self-interaction a <> a
                if (a == b) {
Carl Poelking's avatar
Carl Poelking committed
941
                    w12 *= _cutoff->getCenterWeight();
Carl Poelking's avatar
Carl Poelking committed
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
                    // Source 'a', target 'a'
                    InteractUpdateSourceTarget(a, b, T_ab,
                        sigma1, sigma2,
                        w12, w2,
                        k, L_k_in, L_k_out,
                        Lambda_flat_k_in, Lambda_flat_k_out,
                        LM_in, LM_out,
                        false);
                }
                // Cross-interaction a <> b
                else {
                    // Source 'a', target 'b'
                    InteractUpdateSourceTarget(a, b, T_ab,
                        sigma1, sigma2,
                        w12, w2,
                        k, L_k_in, L_k_out,
                        Lambda_flat_k_in, Lambda_flat_k_out,
                        LM_in, LM_out,
                        true);
                    // Source 'b', target 'a'
                    InteractUpdateSourceTarget(b, a, T_ab,
                        sigma1, sigma2,
                        w12, w1,
                        k, L_k_in, L_k_out,
                        Lambda_flat_k_in, Lambda_flat_k_out,
                        LM_in, LM_out,
                        false);
                }
            }
        }
        // Update dimension of flat indices
        Lambda_flat_k_in = Lambda_flat_k_in*(L_k_in+1);
        L_k_in = L;
    }
Carl Poelking's avatar
Carl Poelking committed
976

Carl Poelking's avatar
Carl Poelking committed
977
978
979
980
981
982
983
984
985
986
987
988
    if (_options->get<std::string>("fieldtensor.contract") == "normal") {
        for (auto it = beginAtomic(); it != endAtomic(); ++it) {
            (*it)->contract();
        }
    }
    else if (_options->get<std::string>("fieldtensor.contract") == "deep") {
        for (auto it = beginAtomic(); it != endAtomic(); ++it) {
            (*it)->contractDeep();
        }
    }
    else {
        assert(false && "fieldtensor.contract has invalid value (should be 'normal' or 'deep')");
Carl Poelking's avatar
Carl Poelking committed
989
    }
Carl Poelking's avatar
Carl Poelking committed
990

Carl Poelking's avatar
Carl Poelking committed
991
    /*
Carl Poelking's avatar
Carl Poelking committed
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
    //Tlmlm T12(L);
    vec d12(0.,0.,1.);
    double a12 = 0.5;
    int L1 = 2;
    int L2 = 2;
    for (double r12 = 0.; r12 <= 10.; r12 += 0.05) {
        Tlmlm::coeff_t T_mat;
        T12.computeTlmlm(d12, r12, a12, L1, L2, T_mat);

        int l1 = 1;
        int m1 = 0;
        int l2 = 1;
        int m2 = 0;

        int lm1 = l1*l1+l1+m1;
        int lm2 = l2*l2+l2+m2;
        GLOG() << r12 << " " << T_mat(lm1,lm2).real() << " " << T_mat(lm2,lm1).real()
        << " " << T_mat(lm1,lm2).imag() << " " << T_mat(lm2,lm1).imag() << std::endl;
    }
Carl Poelking's avatar
Carl Poelking committed
1011
    */
Carl Poelking's avatar
Carl Poelking committed
1012
1013

    /*GLOG() << "Factorial2 ..." << std::endl;
Carl Poelking's avatar
Carl Poelking committed
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
    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;
        }
    }
Carl Poelking's avatar
Carl Poelking committed
1063

Carl Poelking's avatar
Carl Poelking committed
1064
    GLOG() << "Legendre ..." << std::endl;
Carl Poelking's avatar
Carl Poelking committed
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
    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;
        }
    }
Carl Poelking's avatar
Carl Poelking committed
1096
1097
    */

Carl Poelking's avatar
Carl Poelking committed
1098
1099
1100
1101
1102
1103
    return;
}

void FTSpectrum::registerPython() {
    using namespace boost::python;
    class_<FTSpectrum>("FTSpectrum", init<Structure &, Options &>())
Carl Poelking's avatar
Carl Poelking committed
1104
1105
        .def("__iter__", range<return_value_policy<reference_existing_object> >(
            &FTSpectrum::beginAtomic, &FTSpectrum::endAtomic))
Carl Poelking's avatar
Carl Poelking committed
1106
        .def("polarize", &FTSpectrum::polarize)
Carl Poelking's avatar
Carl Poelking committed
1107
1108
1109
1110
        .def("compute", &FTSpectrum::compute);
    return;
}

Carl Poelking's avatar
Carl Poelking committed
1111
} // Namespace