atomicspectrum.cpp 15.7 KB
Newer Older
Carl Poelking's avatar
Carl Poelking committed
1
2
3
4
5
#include <fstream>
#include <boost/format.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>

6
#include "soap/atomicspectrum.hpp"
Carl Poelking's avatar
Carl Poelking committed
7
8
9

namespace soap {

10
AtomicSpectrum::AtomicSpectrum(Particle *center, Basis *basis) {
11
12
    this->null();
    _center = center;
13
    _center_id = center->getId();
14
15
16
17
    _center_pos = center->getPos();
    _center_type = center->getType();
    _basis = basis;
    _qnlm_generic = new BasisExpansion(_basis);
18
19
}

20
21
22
23
24
25
AtomicSpectrum::AtomicSpectrum(Basis *basis) {
    this->null();
    _basis = basis;
    _qnlm_generic = new BasisExpansion(_basis);
}

26
AtomicSpectrum::~AtomicSpectrum() {
Carl Poelking's avatar
Carl Poelking committed
27
28
29
    // CLEAN QNLM'S
    // ... Summed density expansions
    for (auto it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) delete it->second;
30
    _map_qnlm.clear();
Carl Poelking's avatar
Carl Poelking committed
31
    // ... Generic density
32
33
34
35
    if (_qnlm_generic) {
        delete _qnlm_generic;
        _qnlm_generic = NULL;
    }
Carl Poelking's avatar
Carl Poelking committed
36
37
38
    // CLEAN XNKL'S
    // ... Scalar spectra
    for (auto it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) delete it->second;
39
    _map_xnkl.clear();
Carl Poelking's avatar
Carl Poelking committed
40
41

    // ... Generic power spectra
42
43
44
45
46
47
48
49
    if (_xnkl_generic_coherent) {
        delete _xnkl_generic_coherent;
        _xnkl_generic_coherent = NULL;
    }
    if (_xnkl_generic_incoherent) {
        delete _xnkl_generic_incoherent;
        _xnkl_generic_incoherent = NULL;
    }
Carl Poelking's avatar
Carl Poelking committed
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
    // CLEAN PID-RESOLVED GRADIENTS
    // ... Neighbour density expansions with gradients
    for (auto it = _map_pid_qnlm.begin(); it != _map_pid_qnlm.end(); ++it) delete it->second.second;
    _map_pid_qnlm.clear();
    // ... Xnkl gradients
    for (auto it = _map_pid_xnkl.begin(); it != _map_pid_xnkl.end(); ++it) {
        for (auto jt = (*it).second.begin(); jt != (*it).second.end(); ++jt) {
            delete (*jt).second;
        }
        (*it).second.clear();
    }
    _map_pid_xnkl.clear();
    // ... Gradients (generic-coherent)
    for (auto it = _map_pid_xnkl_gc.begin(); it != _map_pid_xnkl_gc.end(); ++it) {
        delete it->second;
    }
66
67
68
}

void AtomicSpectrum::null() {
69
    _center = NULL;
70
    _center_id = -1;
71
72
73
74
75
76
    _center_pos = vec(0,0,0);
    _center_type = "?";
    _basis = NULL;
    _qnlm_generic = NULL;
    _xnkl_generic_coherent = NULL;
    _xnkl_generic_incoherent = NULL;
77
78
}

Carl Poelking's avatar
Carl Poelking committed
79
80
81
82
83
84
void AtomicSpectrum::invert(map_xnkl_t &map_xnkl, xnkl_t *xnkl_generic_coherent, std::string type1, std::string type2) {

    int N = _basis->getRadBasis()->N();
    int L = _basis->getAngBasis()->L();

    assert(xnkl_generic_coherent->getBasis() == _basis &&
85
        "Trying to invert spectrum linked against foreign basis.");
Carl Poelking's avatar
Carl Poelking committed
86
87
88
89
90
91
92
93
94
95
96

    PowerExpansion::coeff_t &xnkl = xnkl_generic_coherent->getCoefficients();
    BasisExpansion::coeff_t &qnlm = _qnlm_generic->getCoefficients();

    type_pair_t types(type1, type2);
    if (type1 == "g" && type2 == "c") xnkl = xnkl_generic_coherent->getCoefficients();
    else if (type1 == "g" && type2 == "i") throw soap::base::NotImplemented("::invert g/i");
    else xnkl = map_xnkl[types]->getCoefficients();

    // ZERO QNLM TO BE SAFE
    for (int n= 0; n < N; ++n) {
97
98
99
100
101
        for (int l = 0; l <= L; ++l) {
            for (int m = -l; m <= l; ++m) {
                qnlm(n, l*l+l+m) = std::complex<double>(0.,0.);
            }
        }
Carl Poelking's avatar
Carl Poelking committed
102
103
104
105
106
107
108
    }
    // Q000 from X000, then Qk00 = X0k0/Q000
    typedef std::complex<double> cmplx;
    int n = 0;
    int k = 0;
    int l = 0;
    int m = 0;
109
110
111
112
//    qnlm(n, l*l+l+m) = cmplx(sqrt(xnkl(n*N+k, l).real()), 0.);
//    for (k = 1; k < N; ++k) {
//        qnlm(k, l*l+l+m) = cmplx(xnkl(n*N+k, l).real()/qnlm(0, l*l+l+m).real(), 0.);
//    }
Carl Poelking's avatar
Carl Poelking committed
113
114
115

//    // FILL Qn00's USING Xnn0's
//    for (int n = 0; n < N; ++n) {
116
117
118
//        double xnn0 = sqrt(xnkl(n*N+n,0).real());
//        std::cout << xnn0 << std::endl;
//        qnlm(n, 0) = std::complex<double>(xnn0, 0.);
Carl Poelking's avatar
Carl Poelking committed
119
120
121
//    }
//    for (int n = 0; n < N; ++n) {
//        for (int l = 0; l <= L; ++l) {
122
123
//            double xnnl = sqrt(xnkl(n*N+n,l).real());
//            qnlm(n, l*l+l) = std::complex<double>(xnnl, 0.);
Carl Poelking's avatar
Carl Poelking committed
124
125
126
127
128
129
//        }
//    }

    return;
}

130
void AtomicSpectrum::addQnlm(std::string type, qnlm_t &nb_expansion) {
131
132
133
134
135
136
137
138
139
140
    assert(nb_expansion.getBasis() == _basis &&
        "Should not sum expansions linked against different bases.");
    map_qnlm_t::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;
141
142
}

Carl Poelking's avatar
Carl Poelking committed
143
144
145
146
147
148
149
void AtomicSpectrum::addQnlmNeighbour(Particle *nb, qnlm_t *nb_expansion) {
    std::string type = nb->getType();
    int id = nb->getId();
    this->addQnlm(type, *nb_expansion);

    if (nb == this->getCenter()) {
        delete nb_expansion; // <- gradients should be zero, do not store
150
        //std::cout << "DO NOT STORE" << std::endl;
Carl Poelking's avatar
Carl Poelking committed
151
152
153
154
    }
    else {
        auto it = _map_pid_qnlm.find(id);
        if (it != _map_pid_qnlm.end()) {
155
156
157
158
159
160
161
162
163
164
165
166

            // There is already an entry for this pid - hence, this must be an image of the actual particle.
            // Add coefficients & gradients to existing density expansion.
            // This sanity check is no longer adequate:
            // throw soap::base::SanityCheckFailed("<AtomicSpectrum::addQnlm> Already have entry for pid.");
            _map_pid_qnlm[id].second->add(*nb_expansion);
            if (nb_expansion->hasGradients()) {
                _map_pid_qnlm[id].second->addGradient(*nb_expansion);
            }
        }
        else {
            _map_pid_qnlm[id] = std::pair<std::string,qnlm_t*>(type, nb_expansion);
Carl Poelking's avatar
Carl Poelking committed
167
168
169
170
171
        }
    }
    return;
}

Carl Poelking's avatar
Carl Poelking committed
172
void AtomicSpectrum::mergeQnlm(AtomicSpectrum *other, double scale, bool gradients) {
173
174
175
    // Function used to construct global spectrum as sum over atomic spectra.
    // The result is itself an "atomic" spectrum (as data fields are largely identical,
    // except for the fact that this summed spectrum does not have a well-defined center.
176
177
178
179
    // The <scale> factor modifies the qnlm, whereas their gradients remain unaltered.
    // For global spectrum, <scale> should be 0.5 in order to guarantee consistency
    // between the scalars qnlm and their gradients (or, in other words, the 0.5 ensures
    // that the global spectrum only counts all pairs once).
180
181
182
    assert(other->getBasis() == _basis &&
        "Should not merge atomic spectra linked against different bases.");
    // Type-agnostic (=generic) density expansion
183
    _qnlm_generic->add(*other->getQnlmGeneric(), scale);
184
185
186
187
188
189
190
191
192
193
194
195
    // Type-resolved (=specific) density expansions
    map_qnlm_t &map_qnlm_other = other->getQnlmMap();
    for (auto it = map_qnlm_other.begin(); it != map_qnlm_other.end(); ++it) {
        std::string density_type = it->first;
        BasisExpansion *density = it->second;
        // Already have density of this type?
        auto mit = _map_qnlm.find(density_type);
        if (mit == _map_qnlm.end()) {
            _map_qnlm[density_type] = new qnlm_t(_basis);
            mit = _map_qnlm.find(density_type);
        }
        // Add ...
196
        mit->second->add(*density, scale);
197
198
    }
    // Particle-ID-resolved gradients
Carl Poelking's avatar
Carl Poelking committed
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
    if (gradients) {
        map_pid_qnlm_t &map_pid_qnlm_other = other->getPidQnlmMap();
        for (auto it = map_pid_qnlm_other.begin(); it != map_pid_qnlm_other.end(); ++it) {
            int pid = it->first;
            std::string pid_type = it->second.first;
            qnlm_t *density_grad = it->second.second;
            // Already have density gradient for this particle?
            auto mit = _map_pid_qnlm.find(pid);
            if (mit == _map_pid_qnlm.end()) {
                qnlm_t *qnlm = new qnlm_t(_basis);
                // Remember to setup zero matrices to store gradient ...
                qnlm->zeroGradient();
                _map_pid_qnlm[pid] = std::pair<std::string,qnlm_t*>(pid_type, qnlm);
                mit = _map_pid_qnlm.find(pid);
            }
            // Add gradients ...
            mit->second.second->addGradient(*density_grad);
216
217
218
219
220
        }
    }
    return;
}

221
AtomicSpectrum::qnlm_t *AtomicSpectrum::getQnlm(std::string type) {
222
223
224
225
226
227
228
229
230
231
    if (type == "") {
        return _qnlm_generic;
    }
    map_qnlm_t::iterator it = _map_qnlm.find(type);
    if (it == _map_qnlm.end()) {
        throw soap::base::OutOfRange("AtomicSpectrum: No such type '" + type + "'");
    }
    else {
        return it->second;
    }
232
233
}

Carl Poelking's avatar
Carl Poelking committed
234
void AtomicSpectrum::computePowerGradients() {
235
    GLOG() << "CID " << _center_id << ": " << std::endl;
Carl Poelking's avatar
Carl Poelking committed
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
    for (auto it1 = _map_pid_qnlm.begin(); it1 != _map_pid_qnlm.end(); ++it1) {
        // Derivatives are taken with respect to the coordinates of this neighbour particle
        int pi_id = it1->first;
        std::string pi_type = it1->second.first;
        qnlm_t *pi_qnlm = it1->second.second;
        // Check and prepare storage in map
        auto it = _map_pid_xnkl.find(pi_id);
        if (it != _map_pid_xnkl.end()) {
            soap::base::SanityCheckFailed("<AtomicSpectrum::computePowerGradients> Already have entry for pid.");
        }
        _map_pid_xnkl[pi_id] = map_xnkl_t();
        map_xnkl_t &xnkl_map = _map_pid_xnkl[pi_id];
        // Compute for all type pairs
        GLOG() << "    NID " << pi_id << " " << std::flush;
        for (auto it2 = _map_qnlm.begin(); it2 != _map_qnlm.end(); ++it2) {
            std::string type_other = it2->first;
            qnlm_t *sum_qnlm_type_other = it2->second;
            if (pi_type == type_other) {
                // Type 1 == type 2
                GLOG() << pi_type << "=" << type_other << " " << std::flush;
                PowerExpansion *powex = new PowerExpansion(_basis);
                powex->computeCoefficientsGradients(pi_qnlm, sum_qnlm_type_other, true);
                // Store ...
                type_pair_t types(pi_type, type_other);
                auto it = xnkl_map.find(types);
                assert(it == xnkl_map.end());
                xnkl_map[types] = powex;
            }
            else {
                // Type 1 != type 2
                GLOG() << pi_type << ":" << type_other << " " << std::flush;
                PowerExpansion *powex_12 = new PowerExpansion(_basis);
                powex_12->computeCoefficientsGradients(pi_qnlm, sum_qnlm_type_other, false);
                GLOG() << type_other << ":" << pi_type << " " << std::flush;
                PowerExpansion *powex_21 = new PowerExpansion(_basis);
                powex_21->computeCoefficientsGradients(sum_qnlm_type_other, pi_qnlm, false);
                // Store ...
                type_pair_t types_12(pi_type, type_other);
                type_pair_t types_21(type_other, pi_type);
                auto it12 = xnkl_map.find(types_12);
                auto it21 = xnkl_map.find(types_21);
                assert(it12 == xnkl_map.end() && it21 == xnkl_map.end());
                xnkl_map[types_12] = powex_12;
                xnkl_map[types_21] = powex_21;
            }
        }
        if (_qnlm_generic) {
            auto it = _map_pid_xnkl_gc.find(pi_id);
            if (it != _map_pid_xnkl_gc.end()) assert(false && "Already have Xnkl_gc gradients for pid.");
            GLOG() << "*:* " << std::flush;
            PowerExpansion *powex = new PowerExpansion(_basis);
            powex->computeCoefficientsGradients(pi_qnlm, _qnlm_generic, true);
            _map_pid_xnkl_gc[pi_id] = powex;
        }
        GLOG() << std::endl;
    }
    return;
}

295
void AtomicSpectrum::computePower() {
Carl Poelking's avatar
Carl Poelking committed
296
297
    // TODO Calling this function more than once with the same object
    // TODO causes memory leaks => check for existing _map_xnkl entries
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
    // Specific (i.e., type-dependent)
    map_qnlm_t::iterator it1;
    map_qnlm_t::iterator it2;
    for (it1 = _map_qnlm.begin(); it1 != _map_qnlm.end(); ++it1) {
        for (it2 = _map_qnlm.begin(); it2 != _map_qnlm.end(); ++it2) {
            type_pair_t types(it1->first, it2->first);
            GLOG() << " " << types.first << ":" << types.second << std::flush;
            PowerExpansion *powex = new PowerExpansion(_basis);
            powex->computeCoefficients(it1->second, it2->second);
            _map_xnkl[types] = powex;
        }
    }
    // Generic coherent
    GLOG() << " g/c" << std::flush;
    if (_xnkl_generic_coherent) delete _xnkl_generic_coherent;
    _xnkl_generic_coherent = new PowerExpansion(_basis);
    _xnkl_generic_coherent->computeCoefficients(_qnlm_generic, _qnlm_generic);
    // Generic incoherent
    GLOG() << " g/i" << std::flush;
    if (_xnkl_generic_incoherent) delete _xnkl_generic_incoherent;
    _xnkl_generic_incoherent = new PowerExpansion(_basis);
    map_xnkl_t::iterator it;
    for (it = _map_xnkl.begin(); it != _map_xnkl.end(); ++it) {
        _xnkl_generic_incoherent->add(it->second);
    }
    GLOG() << std::endl;
324
325
}

Carl Poelking's avatar
Carl Poelking committed
326
void AtomicSpectrum::write(std::ostream &ofs) {
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
    throw soap::base::NotImplemented("AtomicSpectrum::write");
    map_qnlm_t::iterator it;
    for (it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) {
        std::cout << "TYPE" << it->first << std::endl;
        qnlm_t *qnlm = it->second;
        qnlm_t::coeff_t &coeff = qnlm->getCoefficients();
        int L = qnlm->getBasis()->getAngBasis()->L();
        int N = qnlm->getBasis()->getRadBasis()->N();
        for (int n = 0; n < N; ++n) {
            for (int l = 0; l < (L+1); ++l) {
                for (int m = -l; m <= l; ++l) {
                    std::cout << n << " " << l << " " << m << " " << coeff(n,l*l+l+m) << std::endl;
                }
            }
        }
    }
    return;
Carl Poelking's avatar
Carl Poelking committed
344
345
}

346
AtomicSpectrum::xnkl_t *AtomicSpectrum::getXnkl(type_pair_t &types) {
347
348
    map_xnkl_t::iterator it = _map_xnkl.find(types);
    if (it == _map_xnkl.end()) {
349
        if (types.first == "" and types.second == "") {
350
351
352
353
354
355
356
            return _xnkl_generic_coherent;
        }
        else {
            throw soap::base::OutOfRange("AtomicSpectrum: No such type pair '" + types.first + ":" + types.second + "'");
        }
    }
    return it->second;
357
358
}

Carl Poelking's avatar
Carl Poelking committed
359
AtomicSpectrum::xnkl_t *AtomicSpectrum::getPower(std::string type1, std::string type2) {
360
361
    type_pair_t types(type1, type2);
    return this->getXnkl(types);
Carl Poelking's avatar
Carl Poelking committed
362
363
}

Carl Poelking's avatar
Carl Poelking committed
364
365
366
367
368
369
370
371
372
373
boost::python::list AtomicSpectrum::getTypes() {
    boost::python::list types;
    map_qnlm_t::iterator it;
    for (it = _map_qnlm.begin(); it != _map_qnlm.end(); ++it) {
        std::string type = it->first;
        types.append(type);
    }
    return types;
}

374
375
376
377
378
379
380
381
boost::python::list AtomicSpectrum::getNeighbourPids() {
    boost::python::list pids;
    for (auto it = _map_pid_xnkl_gc.begin(); it != _map_pid_xnkl_gc.end(); ++it) {
        pids.append(it->first);
    }
    return pids;
}

Carl Poelking's avatar
Carl Poelking committed
382
383
384
385
386
387
388
389
void AtomicSpectrum::registerPython() {
    using namespace boost::python;

    // Does not work with boost::python ??
    //xnkl_t *(AtomicSpectrum::*getXnkl_string_string)(std::string, std::string) = &AtomicSpectrum::getPower;

    class_<AtomicSpectrum, AtomicSpectrum*>("AtomicSpectrum", init<>())
        .def(init<Particle*, Basis*>())
Carl Poelking's avatar
Carl Poelking committed
390
        .add_property("basis", make_function(&AtomicSpectrum::getBasis, ref_existing()))
391
392
393
        .def("addLinear", &AtomicSpectrum::addQnlm)
        .def("getLinear", &AtomicSpectrum::getQnlm, return_value_policy<reference_existing_object>())
        .def("computePower", &AtomicSpectrum::computePower)
Carl Poelking's avatar
Carl Poelking committed
394
        .def("getTypes", &AtomicSpectrum::getTypes)
395
        .def("getNeighbourPids", &AtomicSpectrum::getNeighbourPids)
396
        .def("getPower", &AtomicSpectrum::getPower, return_value_policy<reference_existing_object>())
397
        .def("getPowerGradGeneric", &AtomicSpectrum::getPowerGradGeneric, return_value_policy<reference_existing_object>())
398
        .def("getCenter", &AtomicSpectrum::getCenter, return_value_policy<reference_existing_object>())
Carl Poelking's avatar
Carl Poelking committed
399
        .def("getCenterId", &AtomicSpectrum::getCenterId)
400
401
        .def("getCenterType", &AtomicSpectrum::getCenterType, return_value_policy<reference_existing_object>())
        .def("getCenterPos", &AtomicSpectrum::getCenterPos, return_value_policy<reference_existing_object>());
Carl Poelking's avatar
Carl Poelking committed
402
403
}

Carl Poelking's avatar
Carl Poelking committed
404
405
}