"while" implementation for pbc

It requires macro PBC_WHILE to be defined

src/tools/Pbc.cpp

@@ -159,12 +159,14 @@ Vector Pbc::distance(const Vector&v1,const Vector&v2,int*nshifts)const{
   Vector d=delta(v1,v2);
   if(type==unset){
   } else if(type==orthorombic) {
+#ifdef PBC_WHILE
+    for(unsigned i=0;i<3;i++){
+      while(d[i]>hdiag[i]) d[i]-=diag[i];
+      while(d[i]<=mdiag[i]) d[i]+=diag[i];
+    }
+#else
    for(int i=0;i<3;i++) d[i]=Tools::pbc(d[i]*invBox(i,i))*box(i,i);
-// this is another possibility:
-//   for(unsigned i=0;i<3;i++){
-//     while(d[i]>hdiag[i]) d[i]-=diag[i];
-//     while(d[i]<=mdiag[i]) d[i]+=diag[i];
-//   }
+#endif
   } else if(type==generic) {
     Vector s=matmul(d,invReduced);
 // check if images have to be computed:

src/tools/Tools.h

@@ -168,6 +168,11 @@ bool Tools::parseFlag(std::vector<std::string>&line,const std::string&key,bool&v
 /// beware: this brings any number into a pbc that ranges from -0.5 to 0.5
 inline
 double Tools::pbc(double x){
+#ifdef PBC_WHILE
+  while (x>0.5) x-=1.0;
+  while (x<-0.5) x+=1.0;
+  return x;
+#else
   if(std::numeric_limits<int>::round_style == std::round_toward_zero) {
     const double offset=100.0;
     const double y=x+offset;
@@ -176,6 +181,7 @@ double Tools::pbc(double x){
   } else if(std::numeric_limits<int>::round_style == std::round_to_nearest) {
     return x-int(x);
   } else return x-floor(x+0.5);
+#endif
 }
 
 template<typename T>