Some type algebra.

src/atomicspectrum.cpp

@@ -191,6 +191,16 @@ AtomicSpectrum::xnkl_t *AtomicSpectrum::getPower(std::string type1, std::string
     return this->getXnkl(types);
 }
 
+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;
+}
+
 void AtomicSpectrum::registerPython() {
     using namespace boost::python;
 
@@ -199,9 +209,11 @@ void AtomicSpectrum::registerPython() {
 
     class_<AtomicSpectrum, AtomicSpectrum*>("AtomicSpectrum", init<>())
         .def(init<Particle*, Basis*>())
+        .add_property("basis", make_function(&AtomicSpectrum::getBasis, ref_existing()))
         .def("addLinear", &AtomicSpectrum::addQnlm)
         .def("getLinear", &AtomicSpectrum::getQnlm, return_value_policy<reference_existing_object>())
         .def("computePower", &AtomicSpectrum::computePower)
+        .def("getTypes", &AtomicSpectrum::getTypes)
         .def("getPower", &AtomicSpectrum::getPower, return_value_policy<reference_existing_object>())
         .def("getCenter", &AtomicSpectrum::getCenter, return_value_policy<reference_existing_object>())
         .def("getCenterType", &AtomicSpectrum::getCenterType, return_value_policy<reference_existing_object>())

src/atomicspectrum.hpp

@@ -53,6 +53,7 @@ public:
     xnkl_t *getXnklGenericCoherent() { return _xnkl_generic_coherent; }
     xnkl_t *getXnklGenericIncoherent() { return _xnkl_generic_incoherent; }
 
+    boost::python::list getTypes();
     static void registerPython();
 
     template<class Archive>

src/tools/CMakeLists.txt

-install(FILES __init__.py loadwrite.py inverse.py DESTINATION ${LOCAL_INSTALL_DIR}/tools)
+install(FILES __init__.py loadwrite.py inverse.py extract.py DESTINATION ${LOCAL_INSTALL_DIR}/tools)
 

src/tools/__init__.py

 from loadwrite import *
 from inverse import *
+from extract import *
 

src/tools/extract.py

+import numpy as np
+import numpy.linalg
+
+from momo import osio, endl, flush
+
+class Xnklab(object):
+    def __init__(self, atomic, types_global):
+        self.S = len(types_global)
+        self.N = atomic.basis.N
+        self.L = atomic.basis.L
+        S = self.S
+        N = self.N
+        L = self.L
+
+        self.X = np.zeros((S*N*N,S*(L+1)))
+        self.types_global = types_global
+
+        types_atomic = atomic.getTypes()
+        S_atomic = len(types_atomic)
+        for i in range(S_atomic):
+            a = types_atomic[i]
+            sa = types_global.index(a)
+            for j in range(S_atomic):
+                b = types_atomic[j]
+                sb = types_global.index(b)
+                xnklab = atomic.getPower(a,b).array
+                r1 = sa*N*N
+                r2 = r1+N*N
+                c1 = sb*(L+1)
+                c2 = c1+(L+1)
+                self.X[r1:r2,c1:c2] = xnklab.real
+                #print a, b, sa, sb
+                #print self.X
+                #raw_input('...')
+        #print "DONE"
+        return
+    def reduce(self):
+        dim_linear_xnklab = (self.S*self.S+self.S)/2*(self.N*self.N*(self.L+1))
+        dim_linear_xnkl = self.N*self.N*(self.L+1)
+        self.X_linear = np.zeros((dim_linear_xnklab))
+        for sa in range(self.S):
+            for sb in range(sa,self.S):
+                # Find slice in matrix X
+                r1 = sa*self.N*self.N
+                r2 = r1+self.N*self.N
+                c1 = sb*(self.L+1)
+                c2 = c1+(self.L+1)
+                xnkl_linear = self.X[r1:r2,c1:c2].reshape((dim_linear_xnkl))
+                # Map to vector
+                sab = self.S*sa - (sa*sa-sa)/2 + (sb-sa) # Linear summation over upper triagonal section
+                idx0 = sab*dim_linear_xnkl
+                idx1 = (sab+1)*dim_linear_xnkl
+                self.X_linear[idx0:idx1] = xnkl_linear
+                #print sa, sb, sab, idx0, idx1
+                #print self.X_linear
+                #raw_input('...')
+        return self.X_linear
+
+def extract_xnklab(atomic, types, types_pairs):
+    assert False
+    types = atomic.getTypes()
+    types = sorted(types)
+    n_types = len(types)
+    N = atomic.basis.N
+    L = atomic.basis.L
+    xnklab = []
+    for i in range(n_types):
+        for j in range(i, n_types):
+            a = types[i]
+            b = types[j]
+            xnklab = xnklab + atomic.getPower(t1,t2).array
+    return xnklab

test/test_types/pair_1_185.xyz

+    56 
+
+ N 31.7089747 41.3358810 86.1879144 
+ C 31.0963060 42.3161507 86.2185478 
+ N 31.3970531 45.5996138 85.1468239 
+ C 30.9257581 44.6682396 85.6434692 
+ C 30.3374277 43.5219791 86.2588437 
+ C 29.1051424 43.6238866 86.8452780 
+ H 28.6444550 44.6060930 86.7995154 
+ S 26.5114656 43.0332444 87.7821424 
+ C 28.0557940 42.5877650 87.0792750 
+ C 25.9764906 41.3711514 87.7813444 
+ C 28.1114444 41.2134848 86.8713598 
+ C 26.9523627 40.5322804 87.2623076 
+ H 28.9805352 40.7295230 86.4444793 
+ H 26.8266698 39.4615986 87.1706784 
+ S 23.9568568 39.5509530 87.2655544 
+ C 24.6420806 40.8717467 88.1791861 
+ C 22.4560907 39.6397130 88.1695727 
+ C 23.7445178 41.3425147 89.1269129 
+ C 22.5256926 40.6534842 89.1195166 
+ H 23.9674447 42.1589141 89.8009676 
+ H 21.7023754 40.8769762 89.7858237 
+ N 22.8727617 36.4567146 86.2629242 
+ C 22.0151208 37.0052059 86.8114155 
+ N 18.5761651 36.7452785 87.1665305 
+ C 19.6445382 37.1630458 87.3089482 
+ C 20.9614142 37.6845811 87.4894200 
+ C 21.1535122 38.7831753 88.2825193 
+ H 20.2573205 39.1962221 88.7356038 
+ N 19.1184299 44.4606155 84.8200553 
+ C 19.8203483 43.7586971 85.4132258 
+ N 21.9330886 43.9144593 88.1577179 
+ C 21.3729007 43.4590333 87.2549388 
+ C 20.6843079 42.8896901 86.1411975 
+ C 20.8756715 41.5716588 85.8265692 
+ H 21.5634575 41.0315692 86.4702384 
+ S 19.2379633 41.4329068 83.4498397 
+ C 20.4513290 40.9071150 84.6025371 
+ C 19.2139267 39.8633124 82.6852826 
+ C 20.7913095 39.5782930 84.3709831 
+ C 20.1021994 38.9935283 83.3016247 
+ H 21.5250027 39.0574350 84.9727855 
+ H 20.2449165 37.9687937 82.9852308 
+ S 18.6109755 38.3455761 80.4189344 
+ C 18.4701380 39.7438911 81.4550959 
+ C 17.5286776 39.0538128 79.2338359 
+ C 17.6279332 40.6957514 80.8982471 
+ C 17.1039831 40.3083420 79.6589961 
+ H 17.4016988 41.6397519 81.3759868 
+ H 16.4286554 40.9179682 79.0724356 
+ N 15.2367960 40.9679090 76.8556936 
+ C 15.6011496 39.8748479 76.9541676 
+ N 14.4666724 36.6838087 76.2123351 
+ C 15.1748066 37.5118217 76.5985657 
+ C 16.0545982 38.5292386 77.0776911 
+ C 17.2620658 38.1820877 77.6202088 
+ H 17.4738301 37.1174875 77.6491110 

test/test_types/test.py

+#! /usr/bin/env python
+import soap
+import soap.tools
+
+import os
+import numpy as np
+from momo import osio, endl, flush
+
+options = soap.Options()
+options.excludeCenters(['H'])
+options.excludeTargets([])
+options.set('radialbasis.type', 'gaussian')
+options.set('radialbasis.mode', 'adaptive')
+options.set('radialbasis.N', 2) # 9
+options.set('radialbasis.sigma', 9) # 0.9
+options.set('radialbasis.integration_steps', 15)
+options.set('radialcutoff.Rc', 6.8)
+options.set('radialcutoff.Rc_width', 0.5)
+options.set('radialcutoff.type', 'shifted-cosine')
+options.set('radialcutoff.center_weight', 1.)
+options.set('angularbasis.type', 'spherical-harmonic')
+options.set('angularbasis.L', 2) # 6
+options.set('densitygrid.N', 20)
+options.set('densitygrid.dx', 0.15)
+
+# STRUCTURE
+xyzfile = 'pair_1_185.xyz'
+top = [('dcv2t', 2, 28)]
+config = soap.tools.ase_load_single(xyzfile)
+structure = soap.tools.setup_structure_ase(config.config_file, config.atoms, top)
+
+# SPECTRUM
+spectrum = soap.Spectrum(structure, options)
+spectrum.compute()
+spectrum.computePower()
+
+types_global = ['H', 'C', 'N', 'S']
+for atomic in spectrum:
+    xnklab = soap.tools.Xnklab(atomic, types_global).reduce()
+    print xnklab
+    raw_input('...')
+