Accurate unwrapping for NpT at long times

implement the expression of bon Bülow et al.
https://arxiv.org/pdf/2003.09205.pdf
which is only a slight variant of the one we used before
(hence the very small diff)

Makefile

@@ -1,5 +1,5 @@
 
-VERSION=1.3
+VERSION=1.4
 TCLINC=-I/usr/include/tcl8.6
 TCLLIB=-I/usr/lib64
 PLUGINDIR=${HOME}/lib/vmd/plugins/LINUXAMD64/tcl

README.md

@@ -1,5 +1,9 @@
 # qwrap — Fast PBC wrapping and unwrapping for VMD
 
+## Development news
+As of version 1.4, unwrapping is done according to the algorithm of von Bülow, Bullerjahn and Hummer (https://arxiv.org/pdf/2003.09205.pdf) to obtain correct diffusion coefficient at long times.
+Note that the algorithm used up to version 1.3 was already similar to this one, and gave only small, bounded deviations from the true long-time diffusion.
+
 ## Using without installing
 In a VMD terminal or TkCon:
 ```

qwrap.cpp

@@ -33,12 +33,11 @@ void calc_shift(float *a, const std::vector<float> &b)
   return;
 }
 
-// This version (for unwrap) gets the shift as integers
-// (in units of PBC vectors)
-void calc_shift_int(float *a, const std::vector<float> &b, int shift[3])
+// This version (for unwrap) adds the shift to the shift accumulator
+void add_shift(float *a, const std::vector<float> &b, double shift[3])
 {
   for (int c=0; c<3; c++)
-    shift[c] = int(floor (a[c] / b[c] + 0.5));
+    shift[c] += floor (a[c] / b[c] + 0.5) * b[c];
   return;
 }
 
@@ -120,7 +119,7 @@ static int do_qwrap(ClientData data, Tcl_Interp *interp, int argc, Tcl_Obj * con
   std::vector<int> centerID;
   std::vector<int> selID;
   std::vector<float> prev_pos;
-  std::vector<int> shifts;
+  std::vector<double> shifts;
   std::vector<float> PBC;
   std::vector<int> is_ref;
   float *coords;
@@ -446,16 +445,12 @@ static int do_qwrap(ClientData data, Tcl_Interp *interp, int argc, Tcl_Obj * con
           prev_pos[current_block * 3 + c] = tmp;  // save the refpos for next frame
         }
         if (frame != first_frame) {
-          int shift[3];
           // Get the shift needed to unwrap the reference position, increment the shift counter
-          calc_shift_int(ref_pos, PBC, shift);
-          for (c = 0; c < 3; c++) {
-            shifts[current_block * 3 + c] += shift[c];
-          }
+          add_shift(ref_pos, PBC, &(shifts[current_block * 3]));
 
           // Actually shift all atoms within the unwrapping block
           for (i = start_atom; i < current_atom; i++) {
-            for (c = 0; c < 3; c++) coords[3*selID[i] + c] -= shifts[current_block * 3 + c] * PBC[c];
+            for (c = 0; c < 3; c++) coords[3*selID[i] + c] -= shifts[current_block * 3 + c];
           }
         }