Xnkl python interface, kernel potential.

src/soap/atomicspectrum.cpp

@@ -98,10 +98,10 @@ void AtomicSpectrum::invert(map_xnkl_t &map_xnkl, xnkl_t *xnkl_generic_coherent,
     int k = 0;
     int l = 0;
     int m = 0;
-    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.);
-    }
+//    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.);
+//    }
 
 //    // FILL Qn00's USING Xnn0's
 //    for (int n = 0; n < N; ++n) {
@@ -139,7 +139,7 @@ void AtomicSpectrum::addQnlmNeighbour(Particle *nb, qnlm_t *nb_expansion) {
 
     if (nb == this->getCenter()) {
         delete nb_expansion; // <- gradients should be zero, do not store
-        std::cout << "DO NOT STORE" << std::endl;
+        //std::cout << "DO NOT STORE" << std::endl;
     }
     else {
         auto it = _map_pid_qnlm.find(id);
@@ -304,6 +304,14 @@ boost::python::list AtomicSpectrum::getTypes() {
     return types;
 }
 
+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;
+}
+
 void AtomicSpectrum::registerPython() {
     using namespace boost::python;
 
@@ -317,7 +325,9 @@ void AtomicSpectrum::registerPython() {
         .def("getLinear", &AtomicSpectrum::getQnlm, return_value_policy<reference_existing_object>())
         .def("computePower", &AtomicSpectrum::computePower)
         .def("getTypes", &AtomicSpectrum::getTypes)
+        .def("getNeighbourPids", &AtomicSpectrum::getNeighbourPids)
         .def("getPower", &AtomicSpectrum::getPower, return_value_policy<reference_existing_object>())
+        .def("getPowerGradGeneric", &AtomicSpectrum::getPowerGradGeneric, 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>())
         .def("getCenterPos", &AtomicSpectrum::getCenterPos, return_value_policy<reference_existing_object>());

src/soap/atomicspectrum.hpp

@@ -59,12 +59,14 @@ public:
     void computePower();
     void computePowerGradients();
     xnkl_t *getPower(std::string type1, std::string type2);
+    xnkl_t *getPowerGradGeneric(int pid) { return _map_pid_xnkl_gc[pid]; }
     xnkl_t *getXnkl(type_pair_t &types);
     map_xnkl_t &getXnklMap() { return _map_xnkl; }
     xnkl_t *getXnklGenericCoherent() { return _xnkl_generic_coherent; }
     xnkl_t *getXnklGenericIncoherent() { return _xnkl_generic_incoherent; }
 
     boost::python::list getTypes();
+    boost::python::list getNeighbourPids();
     static void registerPython();
 
     template<class Archive>

src/soap/basis.cpp

@@ -36,6 +36,7 @@ Basis::~Basis() {
 void Basis::registerPython() {
 	using namespace boost::python;
 	class_<Basis, Basis*>("Basis", init<>())
+	    .def(init<Options*>())
 	    .add_property("N", make_function(&Basis::N, copy_const()))
 		.add_property("L", make_function(&Basis::L, copy_const()));
 }
@@ -110,7 +111,7 @@ void BasisExpansion::computeCoefficients(double r, vec d, double weight, double
         _angbasis->computeCoefficients(d, r, sigma, _angcoeff, NULL, NULL, NULL);
     }
     else if (_has_scalars && _has_gradients) {
-        std::cout << "GRAD" << std::endl;
+        //std::cout << "GRAD" << std::endl;
         _radbasis->computeCoefficients(d, r, sigma, _radcoeff, &_radcoeff_grad_x, &_radcoeff_grad_y, &_radcoeff_grad_z);
         _angbasis->computeCoefficients(d, r, sigma, _angcoeff, &_angcoeff_grad_x, &_angcoeff_grad_y, &_angcoeff_grad_z);
         _weight_scale_grad = _basis->getCutoff()->calculateGradientWeight(r, d);

src/soap/options.cpp

@@ -70,6 +70,7 @@ void Options::registerPython() {
 	void (Options::*set_int)(std::string, int) = &Options::set;
 	void (Options::*set_double)(std::string, double) = &Options::set;
 	void (Options::*set_string)(std::string, std::string) = &Options::set;
+	std::string (Options::*get_string)(std::string) = &Options::get;
 	//void (Options::*set_bool)(std::string, bool) = &Options::set;
 
 	class_<Options, Options*>("Options")
@@ -79,7 +80,8 @@ void Options::registerPython() {
 		.def("excludeTargets", &Options::excludeTargets)
 		.def("set", set_int)
 		.def("set", set_double)
-		.def("set", set_string);
+		.def("set", set_string)
+		.def("get", get_string);
 }
 
 } /* Close namespaces */

src/soap/power.cpp

@@ -3,6 +3,8 @@
 
 namespace soap {
 
+const std::string PowerExpansion::_numpy_t = "complex128";
+
 PowerExpansion::PowerExpansion(Basis *basis) :
     _basis(basis),
     _L(basis->getAngBasis()->L()),
@@ -69,9 +71,9 @@ void PowerExpansion::computeCoefficientsGradients(BasisExpansion *basex1, BasisE
         for (int n = 0; n < _N; ++n) {
             for (int k = 0; k < _N; ++k) {
                 for (int l = 0; l < (_L+1); ++l) {
-                    std::complex<double> c_nkl_dx = 0.0;
-                    std::complex<double> c_nkl_dy = 0.0;
-                    std::complex<double> c_nkl_dz = 0.0;
+                    dtype_t c_nkl_dx = 0.0;
+                    dtype_t c_nkl_dy = 0.0;
+                    dtype_t c_nkl_dz = 0.0;
                     for (int m = -l; m <= l; ++m) {
                         int lm = l*l+l+m;
                         c_nkl_dx += dqnlm_dx(n,lm)*std::conj(qnlm(k,lm))
@@ -91,9 +93,9 @@ void PowerExpansion::computeCoefficientsGradients(BasisExpansion *basex1, BasisE
         for (int n = 0; n < _N; ++n) {
             for (int k = 0; k < _N; ++k) {
                 for (int l = 0; l < (_L+1); ++l) {
-                    std::complex<double> c_nkl_dx = 0.0;
-                    std::complex<double> c_nkl_dy = 0.0;
-                    std::complex<double> c_nkl_dz = 0.0;
+                    dtype_t c_nkl_dx = 0.0;
+                    dtype_t c_nkl_dy = 0.0;
+                    dtype_t c_nkl_dz = 0.0;
                     if (grad_first) {
                         for (int m = -l; m <= l; ++m) {
                             int lm = l*l+l+m;
@@ -131,7 +133,7 @@ void PowerExpansion::computeCoefficients(BasisExpansion *basex1, BasisExpansion
         for (int k = 0; k < _N; ++k) {
             for (int l = 0; l < (_L+1); ++l) {
                 //std::cout << n << " " << k << " " << l << " : " << std::flush;
-                std::complex<double> c_nkl = 0.0;
+                dtype_t c_nkl = 0.0;
                 for (int m = -l; m <= l; ++m) {
                     //std::cout << m << " " << std::flush;
                     c_nkl += coeff1(n, l*l+l+m)*std::conj(coeff2(k, l*l+l+m));
@@ -153,8 +155,8 @@ void PowerExpansion::add(PowerExpansion *other) {
 }
 
 void PowerExpansion::setCoefficientsNumpy(boost::python::object &np_array) {
-    soap::linalg::numpy_converter npc("complex128");
-    npc.numpy_to_ublas< std::complex<double> >(np_array, _coeff);
+    soap::linalg::numpy_converter npc(_numpy_t.c_str());
+    npc.numpy_to_ublas< dtype_t >(np_array, _coeff);
     if (_coeff.size1() != _N*_N ||
         _coeff.size2() != _L+1) {
         throw soap::base::APIError("<PowerExpansion::setCoefficientsNumpy> Matrix size not consistent with basis.");
@@ -162,10 +164,27 @@ void PowerExpansion::setCoefficientsNumpy(boost::python::object &np_array) {
 }
 
 boost::python::object PowerExpansion::getCoefficientsNumpy() {
-    soap::linalg::numpy_converter npc("complex128");
-    return npc.ublas_to_numpy< std::complex<double> >(_coeff);
+    soap::linalg::numpy_converter npc(_numpy_t.c_str());
+    return npc.ublas_to_numpy< dtype_t >(_coeff);
+}
+
+boost::python::object PowerExpansion::getCoefficientsGradXNumpy() {
+    soap::linalg::numpy_converter npc(_numpy_t.c_str());
+    return npc.ublas_to_numpy< dtype_t >(_coeff_grad_x);
+}
+
+boost::python::object PowerExpansion::getCoefficientsGradYNumpy() {
+    soap::linalg::numpy_converter npc(_numpy_t.c_str());
+    return npc.ublas_to_numpy< dtype_t >(_coeff_grad_y);
 }
 
+boost::python::object PowerExpansion::getCoefficientsGradZNumpy() {
+    soap::linalg::numpy_converter npc(_numpy_t.c_str());
+    return npc.ublas_to_numpy< dtype_t >(_coeff_grad_z);
+}
+
+
+
 void PowerExpansion::writeDensity(
     std::string filename,
     Options *options,
@@ -183,7 +202,7 @@ void PowerExpansion::writeDensity(
     for (int n = 0; n < _N; ++n) {
         for (int k = 0; k < _N; ++k) {
             for (int l = 0; l <= _L; ++l) {
-                std::complex<double> c_nkl = coeff(n*_N+k, l);
+                dtype_t c_nkl = coeff(n*_N+k, l);
                 double c_nkl_real = c_nkl.real();
                 double c_nkl_imag = c_nkl.imag();
                 sum_intensity += c_nkl_real*c_nkl_real + c_nkl_imag*c_nkl_imag;
@@ -202,7 +221,10 @@ void PowerExpansion::registerPython() {
 
     class_<PowerExpansion, PowerExpansion*>("PowerExpansion", init<Basis*>())
         .add_property("array", &PowerExpansion::getCoefficientsNumpy, &PowerExpansion::setCoefficientsNumpy)
-        .def("getArray", &PowerExpansion::getCoefficientsNumpy);
+        .def("getArray", &PowerExpansion::getCoefficientsNumpy)
+        .def("getArrayGradX", &PowerExpansion::getCoefficientsGradXNumpy)
+        .def("getArrayGradY", &PowerExpansion::getCoefficientsGradYNumpy)
+        .def("getArrayGradZ", &PowerExpansion::getCoefficientsGradZNumpy);
 }
 
 }

src/soap/power.hpp

@@ -17,8 +17,12 @@ namespace ub = boost::numeric::ublas;
 class PowerExpansion
 {
 public:
-	typedef ub::matrix< std::complex<double> > coeff_t;
-	typedef ub::zero_matrix< std::complex<double> > coeff_zero_t;
+    typedef std::complex<double> dtype_t;
+	typedef ub::matrix< dtype_t > coeff_t;
+	typedef ub::zero_matrix< dtype_t > coeff_zero_t;
+
+	static const std::string _numpy_t;
+    static constexpr double IMAG_EPSILON = 1e-15;
 
 	PowerExpansion() : _basis(NULL), _L(-1), _N(-1), _has_scalars(false), _has_gradients(false) {;}
     PowerExpansion(Basis *basis);
@@ -35,6 +39,9 @@ public:
 
     void setCoefficientsNumpy(boost::python::object &np_array);
     boost::python::object getCoefficientsNumpy();
+    boost::python::object getCoefficientsGradXNumpy();
+    boost::python::object getCoefficientsGradYNumpy();
+    boost::python::object getCoefficientsGradZNumpy();
     static void registerPython();
 
     template<class Archive>

src/soap/structure.cpp

@@ -163,6 +163,7 @@ void Structure::registerPython() {
 	   .def("__iter__", range<return_value_policy<reference_existing_object> >(&Structure::beginParticles, &Structure::endParticles))
 	   .add_property("particles", range<return_value_policy<reference_existing_object> >(&Structure::beginParticles, &Structure::endParticles))
 	   .add_property("segments", range<return_value_policy<reference_existing_object> >(&Structure::beginSegments, &Structure::endSegments))
+       .add_property("n_particles", &Structure::getNumberOfParticles)
 	   .def("connect", &Structure::connectNumeric)
 	   .add_property("box", &Structure::getBoundaryNumeric, &Structure::setBoundaryNumeric)
 	   .add_property("label", make_function(&Structure::getLabel, copy_non_const()), &Structure::setLabel);

src/soap/structure.hpp

@@ -147,6 +147,7 @@ public:
     particle_array_t &particles() { return _particles; }
     particle_it_t beginParticles() { return _particles.begin(); }
     particle_it_t endParticles() { return _particles.end(); }
+    int getNumberOfParticles() { return _particles.size(); }
 
     // SEGMENT CONTAINER
     segment_array_t &segments() { return _segments; }

test/test_kernelpot/config.xyz

+3
+
+C 0.0 0.0 0.0
+C 1.5 0.0 0.0
+C 0.0 1.5 0.0

test/test_kernelpot/map_energy.py

+#! /usr/bin/env python
+import soap
+import soap.tools
+
+import os
+import numpy as np
+import logging
+from momo import osio, endl, flush
+
+options = soap.Options()
+options.excludeCenters([])
+options.excludeTargets([])
+options.set('radialbasis.type', 'gaussian')
+options.set('radialbasis.mode', 'adaptive')
+options.set('radialbasis.N', 9) # 9
+options.set('radialbasis.sigma', 0.5) # 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', 'heaviside')
+options.set('radialcutoff.center_weight', 0.)
+options.set('angularbasis.type', 'spherical-harmonic')
+options.set('angularbasis.L', 6) # 6
+options.set('densitygrid.N', 20)
+options.set('densitygrid.dx', 0.15)
+options.set('spectrum.gradients', True)
+options.set('kernel.adaptor', 'generic')
+options.set('kernel.type', 'dot')
+options.set('kernel.delta', 0.1)
+options.set('kernel.xi', 1.)
+
+class KernelAdaptorGeneric:
+    def __init__(self, options):
+        return
+    def adapt(self, spectrum, pid_list=None):
+        # IX = [
+        #    X1 ->,
+        #    X2 ->,
+        #    ...
+        # ]
+        IX = np.zeros((0,0), dtype='complex128')
+        dimX = -1
+        # TODO More adaptors:
+        # TODO Particle ID mask, center-type mask
+        # TODO For type-resolved adaptors, increase dimensionality accordingly
+        for atomic_i in spectrum:
+            #print atomic_i.getCenter().pos
+            if pid_list:
+                pid = atomic_i.getCenter().id
+                if not pid in pid_list:
+                    logging.debug("Skip, adapt atomic, pid = %d" % pid)
+                    continue
+            Xi = np.real(atomic_i.getPower("","").array)
+            dimX = Xi.shape[0]*Xi.shape[1]
+            Xi = Xi.reshape((dimX))
+            #print "dim(X) =", dimX
+            if not IX.any():
+                IX = np.copy(Xi)
+                IX.resize((1, dimX))
+            else:
+                i = IX.shape[0]
+                IX.resize((i+1, dimX))
+                IX[-1,:] = Xi
+            #print IX
+        return IX
+    def adaptScalar(self, atomic):
+        X = np.real(atomic.getPower("","").array)
+        dimX = X.shape[0]*X.shape[1]
+        X = np.real(X.reshape((dimX)))
+        return X
+    def adaptGradients(self, atomic, nb_pid):
+        dxnkl_pid = atomic.getPowerGradGeneric(nb_pid)
+        dX_dx = dxnkl_pid.getArrayGradX()
+        dX_dy = dxnkl_pid.getArrayGradY()
+        dX_dz = dxnkl_pid.getArrayGradZ()        
+        dimX = dX_dx.shape[0]*dX_dx.shape[1]
+        dX_dx = np.real(dX_dx.reshape((dimX)))
+        dX_dy = np.real(dX_dy.reshape((dimX)))
+        dX_dz = np.real(dX_dz.reshape((dimX)))
+        return dX_dx, dX_dy, dX_dz
+        
+class KernelFunctionDot(object):
+    def __init__(self, options):
+        self.delta = float(options.get('kernel.delta'))
+        self.xi = float(options.get('kernel.xi'))
+        return
+    def computeDot(self, IX, X, xi, delta):
+        return delta**2 * np.dot(IX,X)**xi
+    def compute(self, IX, X):
+        return self.computeDot(IX, X, self.xi, self.delta)    
+    def computeDerivativeOuter(self, IX, X):
+        c = self.computeDot(IX, X, self.xi-1, self.delta)
+        return self.xi*np.diag(c).dot(IX)    
+   
+KernelAdaptorFactory = { 'generic': KernelAdaptorGeneric }     
+KernelFunctionFactory = { 'dot':KernelFunctionDot }
+
+class KernelPotential:
+    def __init__(self, options):
+        logging.info("Construct kernel potential ...")
+        self.basis = soap.Basis(options)
+        self.options = options        
+        # CORE DATA
+        self.IX = None
+        self.alpha = None
+        self.dimX = None # <- Also used as flag set by first ::acquire call
+        # KERNEL
+        logging.info("Choose kernel function ...")
+        self.kernelfct = KernelFunctionFactory[options.get('kernel.type')](options)   
+        # ADAPTOR
+        logging.info("Choose adaptor ...")
+        self.adaptor = KernelAdaptorFactory[options.get('kernel.adaptor')](options)
+        # INCLUSIONS / EXCLUSIONS
+        self.pid_list_acquire = [1]
+        self.pid_list_force = [1]
+        return
+    def acquire(self, structure, alpha):
+        logging.info("Acquire ...")
+        spectrum = soap.Spectrum(structure, self.options, self.basis)
+        spectrum.compute()
+        spectrum.computePower()
+        spectrum.computePowerGradients()
+        # New X's
+        logging.info("Adapt spectrum ...")
+        IX_acqu = self.adaptor.adapt(spectrum, self.pid_list_acquire)        
+        n_acqu = IX_acqu.shape[0]
+        dim_acqu = IX_acqu.shape[1]
+        # New alpha's
+        alpha_acqu = np.zeros((n_acqu))
+        alpha_acqu.fill(alpha)
+        if not self.dimX:
+            # First time ...
+            self.dimX = dim_acqu
+            self.IX = IX_acqu
+            self.alpha = alpha_acqu
+        else:
+            # Check and extend ...
+            assert self.dimX == dim_acqu # Acquired descr. should match linear dim. of previous descr.'s
+            I = self.IX.shape[0]
+            self.IX.resize((I+n_acqu, self.dimX))
+            self.IX[I:I+n_acqu,:] = IX_acqu
+            self.alpha.resize((I+n_acqu))
+            self.alpha[I:I+n_acqu] = alpha_acqu
+        #print self.alpha
+        #print self.IX
+        logging.info("Acquired %d environments." % n_acqu)
+        return
+    def computeEnergy(self, structure):
+        logging.info("Start energy ...")
+        energy = 0.0
+        spectrum = soap.Spectrum(structure, self.options, self.basis)
+        spectrum.compute()
+        spectrum.computePower()
+        spectrum.computePowerGradients()
+        IX_acqu = self.adaptor.adapt(spectrum, self.pid_list_acquire)        
+        n_acqu = IX_acqu.shape[0]
+        dim_acqu = IX_acqu.shape[1]
+        logging.info("Compute energy from %d atomic environments ..." % n_acqu)
+        for n in range(n_acqu):
+            X = IX_acqu[n]
+            print "X_MAG", np.dot(X,X)
+            ic = self.kernelfct.compute(self.IX, X)
+            energy += self.alpha.dot(ic)
+            print "Projection", ic            
+        return energy            
+    def computeForces(self, structure, verbose=False):
+        logging.info("Start forces ...")
+        if verbose:
+            for p in structure: print p.pos
+        forces = [ np.zeros((3)) for i in range(structure.n_particles) ]
+        logging.info("Compute forces on %d particles ..." % structure.n_particles)
+        # Compute X's, dX's
+        spectrum = soap.Spectrum(structure, self.options, self.basis)
+        spectrum.compute()
+        spectrum.computePower()
+        spectrum.computePowerGradients()
+        # Extract & compute force
+        for atomic in spectrum:
+            pid = atomic.getCenter().id
+            if not pid in self.pid_list_force:
+                logging.debug("Skip forces derived from environement with pid = %d" % pid)
+                continue
+            #if pid != 1: continue
+            nb_pids = atomic.getNeighbourPids()
+            logging.info("  Center %d" % (pid))
+            # neighbour-pid-independent kernel "prevector" (outer derivative)
+            X = self.adaptor.adaptScalar(atomic)            
+            dIC = self.kernelfct.computeDerivativeOuter(self.IX, X)
+            alpha_dIC = self.alpha.dot(dIC)
+            for nb_pid in nb_pids:
+                # Force on neighbour
+                logging.info("    -> Nb %d" % (nb_pid))
+                dX_dx, dX_dy, dX_dz = self.adaptor.adaptGradients(atomic, nb_pid)
+                
+                force_x = -alpha_dIC.dot(dX_dx)
+                force_y = -alpha_dIC.dot(dX_dy)
+                force_z = -alpha_dIC.dot(dX_dz)
+                
+                forces[nb_pid-1][0] += force_x
+                forces[nb_pid-1][1] += force_y
+                forces[nb_pid-1][2] += force_z
+                #print forces
+                #raw_input('...')
+        return forces
+        
+
+    
+
+
+
+
+
+def restore_positions(structure, positions):
+    idx = 0
+    for part in structure:
+        part.pos = positions[idx]
+        idx += 1
+    return [ part.pos for part in structure ]
+
+def apply_force_step(structure, forces):
+    max_step = 0.05
+    max_f = 0.0
+    for f in forces:
+        df = np.dot(f,f)**0.5
+        if df > max_f: max_f = df
+    if max_f > max_step: scale = max_step/max_f
+    else: scale = 1.
+    idx = 0
+    for part in structure:
+        #if np.random.uniform(0.,1.) > 0.5: continue
+        part.pos = part.pos + scale*forces[idx]
+        idx += 1
+    return [ part.pos for part in structure ]
+
+
+
+logging.basicConfig(
+    format='[%(asctime)s] %(message)s', 
+    datefmt='%I:%M:%S', 
+    level=logging.ERROR)
+verbose = False
+
+# STRUCTURE
+xyzfile = 'config.xyz'
+config = soap.tools.ase_load_single(xyzfile)
+structure = soap.tools.setup_structure_ase(config.config_file, config.atoms)
+
+for part in structure.particles:
+    print part.id, part.type, part.pos
+    
+positions_0 = [ part.pos for part in structure ]
+
+# KERNEL
+soap.silence()
+kernelpot = KernelPotential(options)
+kernelpot.acquire(structure, 1.)
+
+
+positions = [ 
+    np.array([0.,  0., 0.]),
+    np.array([1.5, 0., 0.]),
+    np.array([1.5, 0., 0.]) ]
+restore_positions(structure, positions)
+print "A", kernelpot.computeEnergy(structure)
+
+positions = [ 
+    np.array([0.,  0.,  0.]),
+    np.array([1.5, 0.,  0.]),
+    np.array([0.,  1.5, 0.]) ]
+restore_positions(structure, positions)
+print "B", kernelpot.computeEnergy(structure)
+
+positions = [ 
+    np.array([0.,  0.,  0.]),
+    np.array([0., 1.5,  0.]),
+    np.array([0.,  0., 1.5]) ]
+restore_positions(structure, positions)
+print "C", kernelpot.computeEnergy(structure)
+
+positions = [ 
+    np.array([0.,  0.,  0.]),
+    np.array([0.,  0.,  1.5]),
+    np.array([1.5, 0.,  0.]) ]
+restore_positions(structure, positions)
+print "D", kernelpot.computeEnergy(structure)
+
+
+
+
+# MAP
+Nx = 20
+Ny = 20
+dx = 0.2
+dy = 0.2
+
+ofs = open('out', 'w')
+for i in range(20):
+    for j in range(20):
+        x = i*dx - 0.5*Nx*dx
+        y = j*dy - 0.5*Ny*dy
+        print x, y
+        positions = [ 
+            np.array([0.,  0., 0.]),
+            np.array([1.5, 0., 0.]),
+            np.array([x,   y,  0.]) ]
+        restore_positions(structure, positions)
+        e_in = kernelpot.computeEnergy(structure)
+        ofs.write('%+1.7e %+1.7e %+1.7e\n' % (x, y, e_in))
+    ofs.write('\n')
+ofs.close()
+

test/test_kernelpot/test.py

+#! /usr/bin/env python
+import soap
+import soap.tools
+
+import os
+import numpy as np
+import logging
+from momo import osio, endl, flush
+
+options = soap.Options()
+options.excludeCenters([])
+options.excludeTargets([])
+options.set('radialbasis.type', 'gaussian')
+options.set('radialbasis.mode', 'adaptive')
+options.set('radialbasis.N', 9) # 9
+options.set('radialbasis.sigma', 0.5) # 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', 'heaviside')
+options.set('radialcutoff.center_weight', 1.)
+options.set('angularbasis.type', 'spherical-harmonic')
+options.set('angularbasis.L', 6) # 6
+options.set('densitygrid.N', 20)
+options.set(