Neighbourlist integer size (#977)

* set up test for the NL size problem

* updated the test

* added a test for 1M atoms

* trying to change to size_t the neighbour list

* forgot the override on the calculate() function

* Squashed commit of the following:

commit 1b70cb0
Author: Daniele Rapetti <[email protected]>
Date:   Tue Oct 10 10:56:35 2023 +0200

    Added an exception to catch the memory error

commit 19f37da
Author: Daniele Rapetti <[email protected]>
Date:   Tue Oct 10 10:26:57 2023 +0200

    impreoved readability for NL

commit 6c50341
Author: Daniele Rapetti <[email protected]>
Date:   Tue Oct 10 09:41:14 2023 +0200

    test now check that the results for initialization are the same with the "simpler" algortithm

commit c469d22
Author: Daniele Rapetti <[email protected]>
Date:   Mon Oct 9 16:58:55 2023 +0200

    refactored the NL constructors

commit 3e8e66c
Author: Daniele Rapetti <[email protected]>
Date:   Mon Oct 9 16:07:53 2023 +0200

    reformatted the test

commit 03e5a41
Author: Daniele Rapetti <[email protected]>
Date:   Mon Oct 9 15:12:21 2023 +0200

    integer type change and small optimization of the initialization of the NL

commit 2d5eb2e
Author: Daniele Rapetti <[email protected]>
Date:   Mon Oct 9 14:09:45 2023 +0200

    added a test to control the NL initialization(atomsID)

commit 04f762c
Author: Daniele Rapetti <[email protected]>
Date:   Mon Oct 9 11:45:17 2023 +0200

    added a test to control the NL initialization(size)

* removed old regtest

* changed to one line per thing also int NL.h

* Now the test will ignore the stackdump if it is plotted

* changed how the "too much memory" for NL on mac works

* updated the changelog

* added a comment to explain the mac exception

* some more explanation on the arbitrary stop

* changed the changelog hoping it fast forward to the one in the master branch

* changed the changelog hoping it fast forward to the one in the master branch

* corrected plumed_merror usage in NL

regtest/basic/rt-NeigbourlistInitialization/Makefile

@@ -0,0 +1 @@
+include ../../scripts/test.make

regtest/basic/rt-NeigbourlistInitialization/config

@@ -0,0 +1 @@
+type=make

regtest/basic/rt-NeigbourlistInitialization/main.cpp

@@ -0,0 +1,132 @@
+#include "plumed/tools/AtomNumber.h"
+#include "plumed/tools/Communicator.h"
+#include "plumed/tools/NeighborList.h"
+#include "plumed/tools/Pbc.h"
+#include <fstream>
+#include <iostream>
+
+using PLMD::AtomNumber;
+using PLMD::Communicator;
+using PLMD::NeighborList;
+using PLMD::Pbc;
+
+// Testing that the Neigbour list will be intialized with the desired number of
+// couples
+// We are initializing with distance and stride not set to check the default
+// parameters
+
+#define check(arg) (((arg)) ? "pass\n" : "not pass\n")
+
+int main(int, char **) {
+  std::ofstream report("unitTest");
+  Pbc pbc{};
+  pbc.setBox(PLMD::Tensor({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}));
+  Communicator cm{};
+  bool serial = true;
+  bool do_pbc = false;
+  for (const size_t nat0 : {100, 500, 1000, 10000}) {
+    std::vector<AtomNumber> list0(nat0);
+    size_t i = 0;
+    for (auto &an : list0) {
+      an.setIndex(i);
+      ++i;
+    }
+    {
+      report << "Single list:\n";
+      std::string prepend="["+std::to_string(nat0)+"]";
+      size_t expected = ((nat0 - 1) * nat0) / 2;
+      auto nl = NeighborList(list0, serial, do_pbc, pbc, cm);
+
+      bool expectedcouples = true;
+      {
+        size_t cID = 0;
+        for (size_t i0 = 0; i0 < nat0 && expectedcouples; ++i0) {
+          for (size_t i1 = i0+1; i1 < nat0 && expectedcouples; ++i1) {  
+            auto couple = nl.getClosePair(cID);
+            expectedcouples &= couple.first == i0;
+            expectedcouples &= couple.second == i1;
+            ++cID;
+          }
+        }
+      }
+      report << prepend << "Initial number:      "
+        << check(nl.size() == expected);
+      report << prepend << "getIndexPair():      "
+        << check(expectedcouples);
+      report << prepend << "Lastupdate is 0:     "
+        << check(nl.getLastUpdate() == 0);
+      report << prepend << "Default stride is 0: "
+          << check(nl.getStride() == 0);
+      report << "\n";
+    }
+    for (const size_t nat1 : {100, 500, 1000, 10000}) {
+
+      std::vector<AtomNumber> list1(nat1);
+
+      i = 0;
+      for (auto &an : list1) {
+        an.setIndex(i);
+        ++i;
+      }
+
+      {
+        report << "Double list, no pairs:\n";
+        std::string prepend="["+std::to_string(nat0)
+            + ", " + std::to_string(nat1)  +"]";
+        bool do_pair = false;
+        size_t expected = nat1 * nat0;
+        auto nl = NeighborList(list0, list1, serial, do_pair, do_pbc, pbc, cm);
+
+        bool expectedcouples = true;
+        {
+          size_t cID = 0;
+          for (size_t i0 = 0; i0 < nat0 && expectedcouples; ++i0) {
+            for (size_t i1 = 0; i1 < nat1 && expectedcouples; ++i1) {
+              auto couple = nl.getClosePair(cID);
+              //The getIndexPair for non couple input must return this be this
+              //(cID / nat1);
+              expectedcouples &= couple.first == i0;
+              //(cID % nat1 + nat0);
+              expectedcouples &= couple.second == nat0+i1;
+              ++cID;
+            }
+          }
+        }
+        report << prepend << "Initial number:      "
+          << check(nl.size() == expected);
+        report << prepend << "getIndexPair():      "
+          << check(expectedcouples);
+        report << prepend << "Lastupdate is 0:     "
+          << check(nl.getLastUpdate() == 0);
+        report << prepend << "Default stride is 0: "
+          << check(nl.getStride() == 0);
+        report << "\n";
+      }
+
+      if (nat1 == nat0) {
+        report << "Double list, with pairs:\n";
+        std::string prepend="["+std::to_string(nat0)
+            + ", " + std::to_string(nat1) +"]";
+        bool do_pair = true;
+        size_t expected = nat0;
+        auto nl = NeighborList(list0, list1, serial, do_pair, do_pbc, pbc, cm);
+
+        bool expectedcouples = true;
+        for (size_t cID = 0; cID < nat0 && expectedcouples; ++cID) {
+          auto couple = nl.getClosePair(cID);
+          expectedcouples &= couple.first == cID;
+          expectedcouples &= couple.second == cID + nat0;
+        }
+        report << prepend << "Initial number:      " 
+          << check(nl.size() == expected);
+        report << prepend << "getIndexPair():      "
+          << check(expectedcouples);
+        report << prepend <<  "Lastupdate is 0:     " 
+          << check(nl.getLastUpdate() == 0);
+        report << prepend << "Default stride is 0: "
+          << check(nl.getStride() == 0);
+        report << "\n";
+      }
+    }
+  }
+}

regtest/basic/rt-NeigbourlistInitialization/unitTest.reference

@@ -0,0 +1,144 @@
+Single list:
+[100]Initial number:      pass
+[100]getIndexPair():      pass
+[100]Lastupdate is 0:     pass
+[100]Default stride is 0: pass
+
+Double list, no pairs:
+[100, 100]Initial number:      pass
+[100, 100]getIndexPair():      pass
+[100, 100]Lastupdate is 0:     pass
+[100, 100]Default stride is 0: pass
+
+Double list, with pairs:
+[100, 100]Initial number:      pass
+[100, 100]getIndexPair():      pass
+[100, 100]Lastupdate is 0:     pass
+[100, 100]Default stride is 0: pass
+
+Double list, no pairs:
+[100, 500]Initial number:      pass
+[100, 500]getIndexPair():      pass
+[100, 500]Lastupdate is 0:     pass
+[100, 500]Default stride is 0: pass
+
+Double list, no pairs:
+[100, 1000]Initial number:      pass
+[100, 1000]getIndexPair():      pass
+[100, 1000]Lastupdate is 0:     pass
+[100, 1000]Default stride is 0: pass
+
+Double list, no pairs:
+[100, 10000]Initial number:      pass
+[100, 10000]getIndexPair():      pass
+[100, 10000]Lastupdate is 0:     pass
+[100, 10000]Default stride is 0: pass
+
+Single list:
+[500]Initial number:      pass
+[500]getIndexPair():      pass
+[500]Lastupdate is 0:     pass
+[500]Default stride is 0: pass
+
+Double list, no pairs:
+[500, 100]Initial number:      pass
+[500, 100]getIndexPair():      pass
+[500, 100]Lastupdate is 0:     pass
+[500, 100]Default stride is 0: pass
+
+Double list, no pairs:
+[500, 500]Initial number:      pass
+[500, 500]getIndexPair():      pass
+[500, 500]Lastupdate is 0:     pass
+[500, 500]Default stride is 0: pass
+
+Double list, with pairs:
+[500, 500]Initial number:      pass
+[500, 500]getIndexPair():      pass
+[500, 500]Lastupdate is 0:     pass
+[500, 500]Default stride is 0: pass
+
+Double list, no pairs:
+[500, 1000]Initial number:      pass
+[500, 1000]getIndexPair():      pass
+[500, 1000]Lastupdate is 0:     pass
+[500, 1000]Default stride is 0: pass
+
+Double list, no pairs:
+[500, 10000]Initial number:      pass
+[500, 10000]getIndexPair():      pass
+[500, 10000]Lastupdate is 0:     pass
+[500, 10000]Default stride is 0: pass
+
+Single list:
+[1000]Initial number:      pass
+[1000]getIndexPair():      pass
+[1000]Lastupdate is 0:     pass
+[1000]Default stride is 0: pass
+
+Double list, no pairs:
+[1000, 100]Initial number:      pass
+[1000, 100]getIndexPair():      pass
+[1000, 100]Lastupdate is 0:     pass
+[1000, 100]Default stride is 0: pass
+
+Double list, no pairs:
+[1000, 500]Initial number:      pass
+[1000, 500]getIndexPair():      pass
+[1000, 500]Lastupdate is 0:     pass
+[1000, 500]Default stride is 0: pass
+
+Double list, no pairs:
+[1000, 1000]Initial number:      pass
+[1000, 1000]getIndexPair():      pass
+[1000, 1000]Lastupdate is 0:     pass
+[1000, 1000]Default stride is 0: pass
+
+Double list, with pairs:
+[1000, 1000]Initial number:      pass
+[1000, 1000]getIndexPair():      pass
+[1000, 1000]Lastupdate is 0:     pass
+[1000, 1000]Default stride is 0: pass
+
+Double list, no pairs:
+[1000, 10000]Initial number:      pass
+[1000, 10000]getIndexPair():      pass
+[1000, 10000]Lastupdate is 0:     pass
+[1000, 10000]Default stride is 0: pass
+
+Single list:
+[10000]Initial number:      pass
+[10000]getIndexPair():      pass
+[10000]Lastupdate is 0:     pass
+[10000]Default stride is 0: pass
+
+Double list, no pairs:
+[10000, 100]Initial number:      pass
+[10000, 100]getIndexPair():      pass
+[10000, 100]Lastupdate is 0:     pass
+[10000, 100]Default stride is 0: pass
+
+Double list, no pairs:
+[10000, 500]Initial number:      pass
+[10000, 500]getIndexPair():      pass
+[10000, 500]Lastupdate is 0:     pass
+[10000, 500]Default stride is 0: pass
+
+Double list, no pairs:
+[10000, 1000]Initial number:      pass
+[10000, 1000]getIndexPair():      pass
+[10000, 1000]Lastupdate is 0:     pass
+[10000, 1000]Default stride is 0: pass
+
+Double list, no pairs:
+[10000, 10000]Initial number:      pass
+[10000, 10000]getIndexPair():      pass
+[10000, 10000]Lastupdate is 0:     pass
+[10000, 10000]Default stride is 0: pass
+
+Double list, with pairs:
+[10000, 10000]Initial number:      pass
+[10000, 10000]getIndexPair():      pass
+[10000, 10000]Lastupdate is 0:     pass
+[10000, 10000]Default stride is 0: pass
+

regtest/basic/rt-NeigbourlistInitializationError/Makefile

@@ -0,0 +1 @@
+include ../../scripts/test.make

regtest/basic/rt-NeigbourlistInitializationError/config

@@ -0,0 +1,10 @@
+type=make
+
+plumed_regtest_after(){
+  #this discards the lines like
+  #"(tools/NeighborList.cpp:98) void PLMD::NeighborList::initialize()"
+  # in this way if NeighborList.cpp is moved or modified this test won't 
+  #trigger a (false) error
+  awk '/(Single|Double|neighbor) list/{print} 
+       /Exception text/{print}' ./unitTest > unitTest.proc
+}

regtest/basic/rt-NeigbourlistInitializationError/main.cpp

@@ -0,0 +1,77 @@
+#include "plumed/tools/AtomNumber.h"
+#include "plumed/tools/Communicator.h"
+#include "plumed/tools/NeighborList.h"
+#include "plumed/tools/Pbc.h"
+#include <fstream>
+#include <iostream>
+
+using PLMD::AtomNumber;
+using PLMD::Communicator;
+using PLMD::NeighborList;
+using PLMD::Pbc;
+
+// Testing that the Neigbour list must throw an explanatory PLMD::Exception
+// when asked to allocate too many pairs
+
+
+#define check(arg) (((arg)) ? "pass\n" : "not pass\n")
+
+int main(int, char **) {
+  std::ofstream report("unitTest");
+  Pbc pbc{};
+  pbc.setBox(PLMD::Tensor({1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}));
+  Communicator cm{};
+  bool serial = true;
+  bool do_pbc = false;
+  // nat0 times nat1 shold be ludicrous big (nat0*nat1*8B~=30~40GB)
+  const size_t nat0=100000;
+  const size_t nat1 = 90000;
+  std::vector<AtomNumber> list0(nat0);
+  size_t i = 0;
+  for (auto &an : list0) {
+    an.setIndex(i);
+    ++i;
+  }
+  {
+    report << "Single list:\n";
+    std::string prepend="["+std::to_string(nat0)+"]";
+    try{
+    size_t expected = ((nat0 - 1) * nat0) / 2;
+    auto nl = NeighborList(list0, serial, do_pbc, pbc, cm);
+    //I need this line to ensure that nl is not optimized away
+    report << prepend << "Initial number:      "
+            << check(nl.size() == expected);
+    } catch ( PLMD::Exception & error ){
+      report << prepend <<"Exception text: "
+              << error.what();
+    }
+    report << "\n";
+  }
+  //doing the same thing with two lists
+
+  std::vector<AtomNumber> list1(nat1);
+  i = 0;
+  for (auto &an : list1) {
+    an.setIndex(i);
+    ++i;
+  }
+
+  {
+    report << "Double list, no pairs:\n";
+    std::string prepend="["+std::to_string(nat0)
+        + ", " + std::to_string(nat1)  +"]";
+    bool do_pair = false;
+    size_t expected = nat1 * nat0;
+    try{
+      auto nl = NeighborList(list0, list1, serial, do_pair, do_pbc, pbc, cm);
+      report << prepend << "Initial number:      "
+              << check(nl.size() == expected);
+
+    } catch( PLMD::Exception & error) {
+      report << prepend <<"Exception text: "
+              << error.what();
+    }
+    report << "\n";   
+  }
+}
+

regtest/basic/rt-NeigbourlistInitializationError/unitTest.proc.reference

@@ -0,0 +1,6 @@
+Single list:
+[100000]Exception text: 
+An error happened while allocating the neighbor list, please decrease the number of atoms used
+Double list, no pairs:
+[100000, 90000]Exception text: 
+An error happened while allocating the neighbor list, please decrease the number of atoms used