From e284c706c85d24c78498c3ffd38d7c15d2d21f88 Mon Sep 17 00:00:00 2001
From: Gareth Aneurin Tribello <garethtribello@Gareths-MacBook-Pro.local>
Date: Wed, 24 Jul 2024 10:16:54 +0100
Subject: [PATCH] Made implementation of link cells that should work even if
 you don't have PBC

---
 src/tools/LinkCells.cpp | 25 +++++++++++++++++++------
 src/tools/LinkCells.h   |  4 ++++
 2 files changed, 23 insertions(+), 6 deletions(-)

diff --git a/src/tools/LinkCells.cpp b/src/tools/LinkCells.cpp
index ef323ed89a..f5c9fa50d2 100644
--- a/src/tools/LinkCells.cpp
+++ b/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)
@@ -45,10 +46,20 @@ 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
-
-  // Setup the pbc object by copying it from action
-  mypbc.setBox( pbc.getBox() );
+  // Create an orthorhombic box around the atomic positions that encompasses every atomic position if there are no pbc
+  if( !pbc.isSet() ) {
+    Vector minp, maxp; minp = maxp = pos[0]; Tensor fake_box; fake_box.zero();
+    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];
+      }
+      fake_box[k][k] = link_cutoff*( 1 + std::ceil( (maxp[k] - minp[k])/link_cutoff ) );
+      origin[k] = ( minp[k] + maxp[k] ) / 2;
+    }
+    mypbc.setBox( fake_box ); nopbc=true;
+    // Setup the pbc object by copying it from action if the Pbc were set
+  } else { mypbc.setBox( pbc.getBox() ); nopbc=false; }
 
   // Setup the lists
   if( pos.size()!=allcells.size() ) {
@@ -153,10 +164,12 @@ void LinkCells::retrieveAtomsInCells( const unsigned& ncells_required,
 }
 
 std::array<unsigned,3> LinkCells::findMyCell( const Vector& pos ) const {
-  Vector fpos=mypbc.realToScaled( pos );
+  Vector mypos = pos; if( nopbc ) mypos = pos - origin;
+  Vector fpos=mypbc.realToScaled( mypos );
   std::array<unsigned,3> celn;
   for(unsigned j=0; j<3; ++j) {
-    celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
+    if( nopbc ) celn[j] = std::floor( fpos[j] * ncells[j] );
+    else 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
   }
   return celn;
diff --git a/src/tools/LinkCells.h b/src/tools/LinkCells.h
index 34ed792bcc..fb2194ddd1 100644
--- a/src/tools/LinkCells.h
+++ b/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