Link cells without PBC (#1100)

* Made implementation of link cells that should work even if you don't have PBC

* Added code to test link cell stuff with nopbc

* Fix of bug introduced after merge of master

* This modification deals with cases with nopbc where link cell cutoff is not really set

* Added error message output for link cells to see problem in waterbridge3

* Attempted fix for LinkCells

* More detailed output in error message to determine problem on github

* Added even more detail to output for link cells

* More debug information for test

* Another attempt to fix problem with new link cell version

---------

Co-authored-by: Gareth Aneurin Tribello <[email protected]>
Co-authored-by: Gareth Aneurin Tribello <[email protected]>

regtest/adjmat/rt-link-nopbc/Makefile

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

regtest/adjmat/rt-link-nopbc/config

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

regtest/adjmat/rt-link-nopbc/main.cpp

@@ -0,0 +1,105 @@
+#include "plumed/tools/Communicator.h"
+#include "plumed/tools/Tools.h"
+#include "plumed/tools/Vector.h"
+#include "plumed/tools/LinkCells.h"
+#include "plumed/tools/Pbc.h"
+#include <fstream>
+#include <iostream>
+
+using namespace PLMD;
+
+void buildCell( const unsigned& num_x, const unsigned& num_y, const unsigned& num_Z, 
+                std::vector<Vector>& fposA, std::vector<Vector>& fposB, 
+                std::vector<unsigned>& Bindices, PLMD::Pbc& mypbc );
+
+bool checkList( const Vector& posA, const std::vector<Vector>& posB, const unsigned& natomsper, const std::vector<unsigned>& indices, std::ofstream& ofs );
+
+int main(){
+
+  std::vector<Vector> fposA, fposB;
+  std::vector<unsigned> tmparray, myatoms, Bindices;
+  unsigned natomsper;
+
+  std::ofstream ofs; ofs.open("logfile");
+
+  PLMD::Communicator comm; PLMD::Pbc mypbc; 
+  PLMD::LinkCells linkcells( comm ); 
+  linkcells.setCutoff( 4.0 );
+  for(unsigned nx=1;nx<6;++nx){
+      for(unsigned ny=1;ny<6;++ny){
+          for(unsigned nz=1;nz<6;++nz){
+              myatoms.resize( 1 + nx*ny*nz*2 );
+              buildCell( nx, ny, nz, fposA, fposB, Bindices, mypbc );
+              // Check list is built correctly - with pbc
+              linkcells.buildCellLists( fposB, Bindices, mypbc );
+              for(unsigned i=0;i<fposA.size();++i){
+                  myatoms[0]=0; natomsper=1;
+                  linkcells.retrieveNeighboringAtoms( fposA[i], tmparray, natomsper, myatoms );
+                  if( checkList( fposA[i], fposB, natomsper, myatoms, ofs ) ) ofs<<"LINK CELLS WORKED SUCCESSFULLY FOR ATOM "<<i<<" WITH PBC FOR "<<nx<<"x"<<ny<<"x"<<nz<<" CELL"<<std::endl;
+                  else ofs<<"LINK CELLS DIDN'T WORK FOR ATOM "<<i<<" WITH PBC FOR "<<nx<<"x"<<ny<<"x"<<nz<<" CELL"<<std::endl;
+              }
+          }
+      }
+  }
+  ofs.close();
+  return 0;
+}
+
+bool checkList( const Vector& posA, const std::vector<Vector>& posB, const unsigned& natomsper, const std::vector<unsigned>& indices, std::ofstream& ofs  ) {
+  for(unsigned i=0; i<posB.size(); ++i) {
+      bool skip = false;
+      for(unsigned j=1; j<natomsper;++j) {
+          if( indices[j]==(i+1) ) { skip=true; break; }
+      }
+      if( skip ) continue ;
+
+      double val = delta(posA, posB[i]).modulo2();
+      if( val<16.0 ) {
+          ofs<<"POINT "<<i<<" "<<val<<std::endl;
+          return false; 
+      }
+  }
+  return true;
+}
+
+void buildCell( const unsigned& num_x, const unsigned& num_y, const unsigned& num_z, 
+                std::vector<Vector>& fposA, std::vector<Vector>& fposB, 
+                std::vector<unsigned>& Bindices, PLMD::Pbc& mypbc ){
+
+  // Atoms in unit cell
+  Vector posA;
+  posA[0]=0.5; posA[1]=0.5; posA[2]=0.5;
+  std::vector<Vector> posB(2);
+  posB[0][0]=0.0; posB[0][1]=0.0; posB[0][2]=0.0;
+  posB[1][0]=0.5; posB[1][1]=1.0; posB[1][2]=0.5;
+
+  fposA.resize( num_x*num_y*num_z );
+  fposB.resize( 2*num_x*num_y*num_z );
+  Bindices.resize( 2*num_x*num_y*num_z );
+
+  unsigned n=0;
+  PLMD::Tensor cell; cell.zero(); 
+  cell[0][0]=4.0; cell[1][1]=4.0; cell[2][2]=4.0;
+  for(unsigned nx=0;nx<num_x;++nx){
+     for(unsigned ny=0;ny<num_y;++ny){
+        for(unsigned nz=0;nz<num_z;++nz){
+            // Shift central atom
+            fposA[n][0]=posA[0]+nx*cell[0][0];
+            fposA[n][1]=posA[1]+ny*cell[1][1];
+            fposA[n][2]=posA[2]+nz*cell[2][2];
+
+            // Shift other atoms
+            Bindices[2*n]=2*n+1;
+            fposB[2*n][0]=posB[0][0]+nx*cell[0][0];
+            fposB[2*n][1]=posB[0][1]+ny*cell[1][1];
+            fposB[2*n][2]=posB[0][2]+nz*cell[2][2];
+            Bindices[2*n+1]=2*n+2;
+            fposB[2*n+1][0]=posB[1][0]+nx*cell[0][0];
+            fposB[2*n+1][1]=posB[1][1]+ny*cell[1][1];
+            fposB[2*n+1][2]=posB[1][2]+nz*cell[2][2];
+            n++;
+        }
+     }
+  }
+  cell[0][0]*=num_x; cell[1][1]*=num_y; cell[2][2]*=num_z;
+}

src/tools/LinkCells.cpp

@@ -28,6 +28,7 @@ namespace PLMD {
 LinkCells::LinkCells( Communicator& cc ) :
   comm(cc),
   cutoffwasset(false),
+  nopbc(false),
   link_cutoff(0.0),
   ncells(3),
   nstride(3) {
@@ -46,19 +47,34 @@ double LinkCells::getCutoff() const {
 void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc ) {
   plumed_assert( cutoffwasset && pos.size()==indices.size() );
 
-  // Must be able to check that pbcs are not nonsensical in some way?? -- GAT
+  // Create an orthorhombic box around the atomic positions that encompasses every atomic position if there are no pbc
   auto box = pbc.getBox();
   if(box(0,0)==0.0 && box(0,1)==0.0 && box(0,2)==0.0 && box(1,0)==0.0 && box(1,1)==0.0 && box(1,2)==0.0 && box(2,0)==0.0 && box(2,1)==0 && box(2,2)==0) {
-    // If the box is not set then we can't use link cells.  We thus set the link cell cutoff and box vectors equal to 23 (because it is the best number).
-    // Setting everything this way ensures that the link cells are a 1x1x1 box.  Notice that if it is a one by one by one box then we are hard coded to return
-    // 0 in findCell. GAT
-    // Note: any non infinite - non zero number should work correctly here! Apparently Gareth likes numerology. I do as well, so let's keep this. GB
-    box(0,0) = box(1,1) = box(2,2) = link_cutoff = 23;
+    Vector minp, maxp;
+    minp = maxp = pos[0];
+    for(unsigned k=0; k<3; ++k) {
+      for(unsigned i=1; i<pos.size(); ++i) {
+        if( pos[i][k]>maxp[k] ) {
+          maxp[k] = pos[i][k];
+        }
+        if( pos[i][k]<minp[k] ) {
+          minp[k] = pos[i][k];
+        }
+      }
+      if( link_cutoff<std::sqrt(std::numeric_limits<double>::max()) ) {
+        box[k][k] = link_cutoff*( 1 + std::ceil( (maxp[k] - minp[k])/link_cutoff ) );
+      } else {
+        box[k][k] = maxp[k] - minp[k] + 1;
+      }
+      origin[k] = ( minp[k] + maxp[k] ) / 2;
+    }
+    nopbc=true;
+    // Setup the pbc object by copying it from action if the Pbc were set
   } else {
     auto determinant = box.determinant();
+    nopbc=false;
     plumed_assert(determinant > epsilon) <<"Cell lists cannot be built when passing a box with null volume. Volume is "<<determinant;
   }
-  // Setup the pbc object by copying it from action
   mypbc.setBox( box );
 
   // Setup the lists
@@ -190,14 +206,14 @@ void LinkCells::retrieveAtomsInCells( const unsigned& ncells_required,
 
 std::array<unsigned,3> LinkCells::findMyCell( const Vector& pos ) const {
   std::array<unsigned,3> celn;
-  if( ncells[0]*ncells[1]*ncells[2] == 1 ) {
-    celn[0]=celn[1]=celn[2]=0;
-    return celn;
+  Vector mypos = pos;
+  if( nopbc ) {
+    mypos = pos - origin;
   }
-  Vector fpos=mypbc.realToScaled( pos );
+  Vector fpos=mypbc.realToScaled( mypos );
   for(unsigned j=0; j<3; ++j) {
     celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
-    plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ); // Check that atom is in box
+    plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ) <<"in link cell "<<celn[j]<<" but should be between 0 and "<<ncells[j]<<" link cell cutoff is "<<link_cutoff<<" position is "<<fpos[0]<<" "<<fpos[1]<<" "<<fpos[2]<<" box is "<<mypbc.getBox()(0,0)<<" "<<mypbc.getBox()(1,1)<<" "<<mypbc.getBox()(2,2);
   }
   return celn;
 }

src/tools/LinkCells.h

@@ -38,10 +38,14 @@ class LinkCells {
   Communicator & comm;
 /// Check that the link cells were set up correctly
   bool cutoffwasset;
+/// Are there periodic boundary conditions setup
+  bool nopbc;
 /// The cutoff to use for the sizes of the cells
   double link_cutoff;
 /// The pbc we are using for link cells
   Pbc mypbc;
+/// The location of the origin if we are not using periodic boundary conditions
+  Vector origin;
 /// The number of cells in each direction
   std::vector<unsigned> ncells;
 /// The number of cells to stride through to get the link cells