diff --git a/src/soap/atomicspectrum.cpp b/src/soap/atomicspectrum.cpp
index 902939bbe92a94764873a5a6cffcf20b24c89079..099344f64dd5c2b3d02409866aa4a32d27d55694 100644
--- a/src/soap/atomicspectrum.cpp
+++ b/src/soap/atomicspectrum.cpp
@@ -159,7 +159,7 @@ void AtomicSpectrum::addQnlmNeighbour(Particle *nb, qnlm_t *nb_expansion) {
     return;
 }
 
-void AtomicSpectrum::mergeQnlm(AtomicSpectrum *other, double scale) {
+void AtomicSpectrum::mergeQnlm(AtomicSpectrum *other, double scale, bool gradients) {
     // 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.
@@ -186,22 +186,24 @@ void AtomicSpectrum::mergeQnlm(AtomicSpectrum *other, double scale) {
         mit->second->add(*density, scale);
     }
     // Particle-ID-resolved 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);
+    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);
         }
-        // Add gradients ...
-        mit->second.second->addGradient(*density_grad);
     }
     return;
 }
diff --git a/src/soap/atomicspectrum.hpp b/src/soap/atomicspectrum.hpp
index d840e234c18df76284e8610d198b792eccd30b16..6b9da75826b49d450d5943e18970dbc41c923fa4 100644
--- a/src/soap/atomicspectrum.hpp
+++ b/src/soap/atomicspectrum.hpp
@@ -59,7 +59,7 @@ public:
     qnlm_t *getQnlmGeneric() { return _qnlm_generic; }
     map_qnlm_t &getQnlmMap() { return _map_qnlm; }
     map_pid_qnlm_t &getPidQnlmMap() { return _map_pid_qnlm; }
-    void mergeQnlm(AtomicSpectrum *other, double scale);
+    void mergeQnlm(AtomicSpectrum *other, double scale, bool gradients);
     // XNKL METHODS
     void computePower();
     void computePowerGradients();
diff --git a/src/soap/cutoff.cpp b/src/soap/cutoff.cpp
index 6eeb759371bb9e7059e41656cc1da767dd5cd979..93caed10bce88cf70196913b144bf959ec3f3283 100644
--- a/src/soap/cutoff.cpp
+++ b/src/soap/cutoff.cpp
@@ -1,4 +1,8 @@
 #include "soap/cutoff.hpp"
+#include <boost/archive/binary_iarchive.hpp>
+#include <boost/archive/binary_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+#include <boost/archive/text_oarchive.hpp>
 
 namespace soap {
 
@@ -77,6 +81,7 @@ void CutoffFunctionFactory::registerAll(void) {
 	CutoffFunctionOutlet().Register<CutoffFunctionHeaviside>("heaviside");
 }
 
-
 }
 
+BOOST_CLASS_EXPORT_IMPLEMENT(soap::CutoffFunctionHeaviside);
+
diff --git a/src/soap/cutoff.hpp b/src/soap/cutoff.hpp
index eb0a0fd577af943f4ede951292873c4cb9ddba8d..ad4b54cda5e8e952ed153180068ffaae0837f75a 100644
--- a/src/soap/cutoff.hpp
+++ b/src/soap/cutoff.hpp
@@ -4,6 +4,8 @@
 #include <string>
 #include <math.h>
 #include <vector>
+#include <boost/serialization/base_object.hpp>
+#include <boost/serialization/export.hpp>
 
 #include "soap/base/exceptions.hpp"
 #include "soap/base/objectfactory.hpp"
@@ -93,6 +95,9 @@ inline CutoffFunction *CutoffFunctionFactory::create(const std::string &key) {
     }
 }
 
-}
+} /* CLOSE NAMESPACE */
+
+BOOST_CLASS_EXPORT_KEY(soap::CutoffFunction);
+BOOST_CLASS_EXPORT_KEY(soap::CutoffFunctionHeaviside);
 
 #endif
diff --git a/src/soap/options.cpp b/src/soap/options.cpp
index 2e7ca88060db5c9d905100a2c0c8dd36b82bb5cc..9f3eeeab92392fe32ab0ce023ac90673726ede1c 100644
--- a/src/soap/options.cpp
+++ b/src/soap/options.cpp
@@ -96,6 +96,30 @@ bool Options::doExcludeTargetId(int pid) {
     return (it == _exclude_target_id.end()) ? false : true;
 }
 
+void Options::generateExclusionLists() {
+    if (!boost::python::len(_exclude_center_list)) {
+        for (auto it = _exclude_center.begin(); it != _exclude_center.end(); ++it) {
+            if (it->second) _exclude_center_list.append(it->first);
+        }
+    }
+    if (!boost::python::len(_exclude_target_list)) {
+        for (auto it = _exclude_target.begin(); it != _exclude_target.end(); ++it) {
+            if (it->second) _exclude_target_list.append(it->first);
+        }
+    }
+    if (!boost::python::len(_exclude_center_id_list)) {
+        for (auto it = _exclude_center_id.begin(); it != _exclude_center_id.end(); ++it) {
+            if (it->second) _exclude_center_id_list.append(it->first);
+        }
+    }
+    if (!boost::python::len(_exclude_target_id_list)) {
+        for (auto it = _exclude_target_id.begin(); it != _exclude_target_id.end(); ++it) {
+            if (it->second) _exclude_target_id_list.append(it->first);
+        }
+    }
+    return;
+}
+
 void Options::registerPython() {
 	using namespace boost::python;
 	void (Options::*set_int)(std::string, int) = &Options::set;
diff --git a/src/soap/options.hpp b/src/soap/options.hpp
index 8d7b6079ea217b0394ee77a63e27670a6a9c31ff..c396ede3fab61f66f7c88a3710a2f1d92aba01a5 100644
--- a/src/soap/options.hpp
+++ b/src/soap/options.hpp
@@ -47,6 +47,9 @@ public:
 	boost::python::list getExcludeCenterIdList() { return _exclude_center_id_list; }
 	boost::python::list getExcludeTargetIdList() { return _exclude_target_id_list; }
 
+    // USED IN SERIALIZATION-LOAD: std::map TO boost::python::list
+    void generateExclusionLists();
+
 	// PYTHON
 	static void registerPython();
 
@@ -57,13 +60,15 @@ public:
 
 		arch & _exclude_center;
 		arch & _exclude_target;
-		arch & _exclude_center_list;
-		arch & _exclude_target_list;
+		//arch & _exclude_center_list;
+		//arch & _exclude_target_list;
 
 		arch & _exclude_center_id;
 		arch & _exclude_target_id;
-		arch & _exclude_center_id_list;
-		arch & _exclude_target_id_list;
+		//arch & _exclude_center_id_list;
+		//arch & _exclude_target_id_list;
+
+        this->generateExclusionLists();
 		return;
 	}
 
diff --git a/src/soap/spectrum.cpp b/src/soap/spectrum.cpp
index 79bcd2536092c3d70d001a65678f172af6efec89..ef33149eebd4168baaa7263a31fc27ca92848e3f 100644
--- a/src/soap/spectrum.cpp
+++ b/src/soap/spectrum.cpp
@@ -129,13 +129,14 @@ AtomicSpectrum *Spectrum::computeGlobal() {
     if (_global_atomic) throw soap::base::APIError("<Spectrum::computeGlobal> Already initialised.");
     _global_atomic = new AtomicSpectrum(_basis);
     GLOG() << "Computing global spectrum ..." << std::endl;
+    bool gradients = _options->get<bool>("spectrum.gradients");
     for (auto it = _atomspec_array.begin(); it != _atomspec_array.end(); ++it) {
         GLOG() << "  Adding center " << (*it)->getCenter()->getId()
             << " (type " << (*it)->getCenter()->getType() << ")" << std::endl;
-        _global_atomic->mergeQnlm(*it, 0.5); // <- Scale factor 0.5 necessary in order not to overcount pairs
+        _global_atomic->mergeQnlm(*it, 0.5, gradients); // <- Scale factor 0.5 necessary in order not to overcount pairs
     }
     _global_atomic->computePower();
-    _global_atomic->computePowerGradients();
+    if (gradients) _global_atomic->computePowerGradients();
     return _global_atomic;
 }