Nearest-image convention abandoned.

src/soap/atomicspectrum.cpp

@@ -152,9 +152,19 @@ void AtomicSpectrum::addQnlmNeighbour(Particle *nb, qnlm_t *nb_expansion) {
     else {
         auto it = _map_pid_qnlm.find(id);
         if (it != _map_pid_qnlm.end()) {
-            throw soap::base::SanityCheckFailed("<AtomicSpectrum::addQnlm> Already have entry for pid.");
+
+            // 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);
         }
-        _map_pid_qnlm[id] = std::pair<std::string,qnlm_t*>(type, nb_expansion);
     }
     return;
 }

src/soap/boundary.cpp

 #include "soap/boundary.hpp"
+#include "soap/globals.hpp"
 #include <boost/archive/binary_iarchive.hpp>
 #include <boost/archive/binary_oarchive.hpp>
 #include <boost/archive/text_iarchive.hpp>
@@ -6,6 +7,107 @@
 
 namespace soap {
 
+
+soap::vec BoundaryTriclinic::connect(const vec &r_i, const vec &r_j) const {
+    /*
+    // This only works if a = (*,0,0), b = (*,*,0), c = (*,*,*) => e.g., GROMACS
+    vec r_tp, r_dp, r_sp, r_ij;
+    vec a = _box.getCol(0); vec b = _box.getCol(1); vec c = _box.getCol(2);
+    r_tp = r_j - r_i;
+    r_dp = r_tp - c*round(r_tp.getZ()/c.getZ());
+    r_sp = r_dp - b*round(r_dp.getY()/b.getY());
+    r_ij = r_sp - a*round(r_sp.getX()/a.getX());
+    return r_ij;
+    */
+
+    vec a = _box.getCol(0);
+    vec b = _box.getCol(1);
+    vec c = _box.getCol(2);
+    //GLOG() << "Box vectors: " << a << " " << b << " " << c << std::endl;
+
+    vec u = _inv_box.getCol(0);
+    vec v = _inv_box.getCol(1);
+    vec w = _inv_box.getCol(2);
+    //GLOG() << "Inverse box vectors: " << u << " " << v << " " << w << std::endl;
+
+    vec dr = r_j - r_i;
+    //GLOG() << "dr " << dr << std::endl;
+    //GLOG() << "u " << std::floor(u*dr) << std::endl;
+    dr = dr - std::floor(u*dr)*a;
+    //GLOG() << "v " << std::floor(v*dr) << std::endl;
+    dr = dr - std::floor(v*dr)*b;
+    //GLOG() << "w " << std::floor(w*dr) << std::endl;
+    dr = dr - std::floor(w*dr)*c;
+    //GLOG() << "in first quadrant " << dr << std::endl;
+
+    vec dr_min = dr;
+    double d_min = soap::linalg::abs(dr);
+
+    for (int i=0; i < 2; ++i) {
+    for (int j=0; j < 2; ++j) {
+    for (int k=0; k < 2; ++k) {
+        vec dr_ijk = dr - i*a - j*b - k*c;
+        double d_ijk = soap::linalg::abs(dr_ijk);
+        //GLOG() << "ijk " << i << " " << j << " " << k << " " << dr << " " << d_ijk << " " << dr_ijk << std::endl;
+        if (d_ijk < d_min) {
+            d_min = d_ijk;
+            dr_min = dr_ijk;
+        }
+    }}}
+
+    return dr_min;
+}
+
+std::vector<int> BoundaryTriclinic::calculateRepetitions(double cutoff) {
+
+    vec a = _box.getCol(0);
+    vec b = _box.getCol(1);
+    vec c = _box.getCol(2);
+
+    vec u = _inv_box.getCol(0);
+    vec v = _inv_box.getCol(1);
+    vec w = _inv_box.getCol(2);
+
+    //std::cout << a << std::endl;
+    //std::cout << b << std::endl;
+    //std::cout << c << std::endl;
+
+    //std::cout << u << std::endl;
+    //std::cout << v << std::endl;
+    //std::cout << w << std::endl;
+
+    double da = std::abs(u*a/soap::linalg::abs(u));
+    double db = std::abs(v*b/soap::linalg::abs(v));
+    double dc = std::abs(w*c/soap::linalg::abs(w));
+
+    //std::cout << da << " " << db << " " << dc << std::endl;
+
+    int na = int(1+(cutoff-0.5*da)/da);
+    int nb = int(1+(cutoff-0.5*db)/db);
+    int nc = int(1+(cutoff-0.5*dc)/dc);
+
+    std::vector<int> na_nb_nc = { na, nb, nc };
+    return na_nb_nc;
+}
+
+std::vector<int> BoundaryOrthorhombic::calculateRepetitions(double cutoff) {
+    vec a = _box.getCol(0);
+    vec b = _box.getCol(1);
+    vec c = _box.getCol(2);
+
+    double da = soap::linalg::abs(a);
+    double db = soap::linalg::abs(b);
+    double dc = soap::linalg::abs(c);
+
+    int na = int(1+(cutoff-0.5*da)/da);
+    int nb = int(1+(cutoff-0.5*db)/db);
+    int nc = int(1+(cutoff-0.5*dc)/dc);
+
+    std::vector<int> na_nb_nc = { na, nb, nc };
+    return na_nb_nc;
+}
+
+
 }
 
 

src/soap/boundary.hpp

@@ -7,6 +7,7 @@
 #include <boost/serialization/export.hpp>
 
 #include "soap/types.hpp"
+#include "soap/base/exceptions.hpp"
 
 namespace soap {
 
@@ -39,12 +40,14 @@ public:
 		vec c = _box.getCol(2);
 		return (a^b)*c;
     }
-    virtual vec connect(const vec &r_i, const vec &r_j) {
+    virtual vec connect(const vec &r_i, const vec &r_j) const {
+        std::cout << "connect default" << std::endl;
     	return r_j - r_i;
     }
 
     virtual std::vector<int> calculateRepetitions(double cutoff) {
         std::vector<int> na_nb_nc = { 0, 0, 0 };
+        return na_nb_nc;
     }
 
     virtual eBoxType getBoxType() { return _type; }
@@ -73,6 +76,7 @@ public:
 		_box.ZeroMatrix();
 	}
     vec connect(const vec &r_i, const vec &r_j) const {
+        std::cout << "connect open" << std::endl;
     	return r_j - r_i;
     }
     template<class Archive>
@@ -93,6 +97,7 @@ public:
 		_box.UnitMatrix();
 	}
 	vec connect(const vec &r_i, const vec &r_j) const {
+	    std::cout << "Connect ortho" << std::endl;
 		vec r_ij;
 		double a = _box.get(0,0); double b = _box.get(1,1); double c = _box.get(2,2);
 		r_ij = r_j - r_i;
@@ -102,6 +107,8 @@ public:
 		return r_ij;
 	}
 
+	virtual std::vector<int> calculateRepetitions(double cutoff);
+
 	template<class Archive>
 	void serialize(Archive &arch, const unsigned int version) {
 		arch & boost::serialization::base_object<Boundary>(*this);
@@ -129,46 +136,9 @@ public:
 		_type = Boundary::typeTriclinic;
 		_box.UnitMatrix();
 	}
-	vec connect(const vec &r_i, const vec &r_j) const {
-        /*
-        // This only works if a = (*,0,0), b = (*,*,0), c = (*,*,*) => e.g., GROMACS
-	    vec r_tp, r_dp, r_sp, r_ij;
-	    vec a = _box.getCol(0); vec b = _box.getCol(1); vec c = _box.getCol(2);
-	    r_tp = r_j - r_i;
-	    r_dp = r_tp - c*round(r_tp.getZ()/c.getZ());
-	    r_sp = r_dp - b*round(r_dp.getY()/b.getY());
-	    r_ij = r_sp - a*round(r_sp.getX()/a.getX());
-	    return r_ij;
-        */
-	    vec a = _box.getCol(0); 
-        vec b = _box.getCol(1); 
-        vec c = _box.getCol(2);
+	virtual vec connect(const vec &r_i, const vec &r_j) const;
 
-	    vec u = _inv_box.getCol(0); 
-        vec v = _inv_box.getCol(1); 
-        vec w = _inv_box.getCol(2);
-
-        vec dr = r_j - r_i;
-        dr = dr - std::floor(u*dr)*a;
-        dr = dr - std::floor(v*dr)*b;
-        dr = dr - std::floor(w*dr)*c;
-        
-        vec dr_min = dr;
-        double d_min = soap::linalg::abs(dr);
-
-        for (int i=0; i < 2; ++i) {
-        for (int j=0; j < 2; ++j) {
-        for (int k=0; k < 2; ++k) {
-            vec dr_ijk = dr - i*a - j*b - k*c;
-            double d_ijk = soap::linalg::abs(dr_ijk);
-            if (d_ijk < d_min) {
-                d_min = d_ijk;
-                dr_min = dr_ijk;
-            }
-        }}}
-
-        return dr_min;
-	}
+	virtual std::vector<int> calculateRepetitions(double cutoff);
 
 	template<class Archive>
 	void serialize(Archive &arch, const unsigned int version) {

src/soap/spectrum.cpp

@@ -92,20 +92,22 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
     GLOG() << "Compute atomic spectrum for particle " << center->getId()
         << " (type " << center->getType() << ", targets " << targets.size() << ") ..." << std::endl;
 
+    // FIND IMAGE REPITIONS REQUIRED TO SATISFY CUTOFF
     vec box_a = _structure->getBoundary()->getBox().getCol(0);
     vec box_b = _structure->getBoundary()->getBox().getCol(1);
     vec box_c = _structure->getBoundary()->getBox().getCol(2);
 
     double rc = _basis->getCutoff()->getCutoff();
-    int na_max = int(1 + rc/box_a.getX() - 0.5);
-    int nb_max = int(1 + rc/box_b.getY() - 0.5);
-    int nc_max = int(1 + rc/box_c.getZ() - 0.5);
-
+    std::vector<int> na_nb_nc = _structure->getBoundary()->calculateRepetitions(rc);
+    int na_max = na_nb_nc[0];
+    int nb_max = na_nb_nc[1];
+    int nc_max = na_nb_nc[2];
 
     GLOG() << box_a << " " << box_b << " " << box_c << std::endl;
     GLOG() << rc << std::endl;
     GLOG() << na_max << " " << nb_max << " " << nc_max << std::endl;
 
+    // CREATE BLANK
     AtomicSpectrum *atomic_spectrum = new AtomicSpectrum(center, this->_basis);
 
     Structure::particle_it_t pit;
@@ -115,26 +117,38 @@ AtomicSpectrum *Spectrum::computeAtomic(Particle *center, Structure::particle_ar
         if (_options->doExcludeTarget((*pit)->getType()) ||
             _options->doExcludeTargetId((*pit)->getId())) continue;
 
+    for (int na=-na_max; na<na_max+1; ++na) {
+    for (int nb=-nb_max; nb<nb_max+1; ++nb) {
+    for (int nc=-nc_max; nc<nc_max+1; ++nc) {
+
+        //GLOG() << na << " " << nb << " " << nc << std::endl;
+        vec L = na*box_a + nb*box_b + nc*box_c;
+
         // FIND DISTANCE & DIRECTION, CHECK CUTOFF
-        vec dr = _structure->connect(center->getPos(), (*pit)->getPos());
+        vec dr = _structure->connect(center->getPos(), (*pit)->getPos()) + L;  // TODO Consider images
         double r = soap::linalg::abs(dr);
         if (! this->_basis->getCutoff()->isWithinCutoff(r)) continue;
         vec d = (r > 0.) ? dr/r : vec(0.,0.,1.);
 
         // APPLY CUTOFF (= WEIGHT REDUCTION)
-        bool is_center = (*pit == center); // TODO Consider images
+        bool is_image = (*pit == center);
+        bool is_center = (*pit == center && na==0 && nb==0 && nc==0); // TODO Consider images
         double weight0 = (*pit)->getWeight();
         double weight_scale = _basis->getCutoff()->calculateWeight(r);
         if (is_center) {
             weight0 *= _basis->getCutoff()->getCenterWeight();
         }
 
+        GLOG() << "C " << dr.getX() << " " << dr.getY() << " " << dr.getZ() << std::endl;
+
         // COMPUTE EXPANSION & ADD TO SPECTRUM
-        bool gradients = (is_center) ? false : _options->get<bool>("spectrum.gradients");
+        bool gradients = (is_image) ? false : _options->get<bool>("spectrum.gradients");
         BasisExpansion *nb_expansion = new BasisExpansion(this->_basis); // <- kept by AtomicSpectrum
         nb_expansion->computeCoefficients(r, d, weight0, weight_scale, (*pit)->getSigma(), gradients);
         atomic_spectrum->addQnlmNeighbour(*pit, nb_expansion); // TODO Consider images
-    }
+
+    }}} // Close loop over images
+    } // Close loop over particles
 
     return atomic_spectrum;
 }