Make pairlisted coordNum variable consistent under biasing forces 98/4698/1
authorGiacomo Fiorin <giacomo.fiorin@gmail.com>
Tue, 16 Oct 2018 20:12:28 +0000 (16:12 -0400)
committerGiacomo Fiorin <giacomo.fiorin@gmail.com>
Tue, 16 Oct 2018 20:12:28 +0000 (16:12 -0400)
See commits from: https://github.com/Colvars/colvars/pull/162

Change-Id: Idf1e1bb0ea5feb0623c2a5edda4ffc084c15db83

colvars/src/colvarcomp_coordnums.cpp
ug/ug_colvars.tex

index 2e2b0c3..ec53391 100644 (file)
@@ -56,11 +56,28 @@ cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
 
   cvm::real const xn = cvm::integer_power(l2, en2);
   cvm::real const xd = cvm::integer_power(l2, ed2);
-  cvm::real const func = (1.0-xn)/(1.0-xd);
+  //The subtraction and division stretches the function back to the range of [0,1] from [pairlist_tol,1]
+  cvm::real const func = (((1.0-xn)/(1.0-xd)) - pairlist_tol) / (1.0-pairlist_tol);
+
+  if (flags & ef_rebuild_pairlist) {
+    //Particles just outside of the cutoff also are considered if they come near.
+    **pairlist_elem = (func > (-pairlist_tol * 0.5)) ? true : false;
+    (*pairlist_elem)++;
+  }
+  //If the value is too small, we need to exclude it, rather than let it contribute to the sum or the gradients.
+  if (func < 0)
+    return 0;
 
   if (flags & ef_gradients) {
-    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) -
-                                            func*ed2*(xd/l2))*(-1.0);
+    //This is the old, completely correct expression for dFdl2:
+    //cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) -
+    //                                        func*ed2*(xd/l2))*(-1.0);
+    //This can become:
+    //cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2)*(1.0-xn)/(1.0-xn) -
+    //                                        func*ed2*(xd/l2))*(-1.0);
+    //Recognizing that func = (1.0-xn)/(1.0-xd), we can group together the "func" and get a version of dFdl2 that is 0
+    //when func=0, which lets us skip this gradient calculation when func=0.
+    cvm::real const dFdl2 = func * ((ed2*xd/((1.0-xd)*l2)) - (en2*xn/((1.0-xn)*l2)));
     cvm::rvector const dl2dx((2.0/((flags & ef_anisotropic) ? r0sq_vec.x :
                                    r0*r0)) * diff.x,
                              (2.0/((flags & ef_anisotropic) ? r0sq_vec.y :
@@ -71,11 +88,6 @@ cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
     A2.grad +=        dFdl2*dl2dx;
   }
 
-  if (flags & ef_rebuild_pairlist) {
-    **pairlist_elem = (func > pairlist_tol) ? true : false;
-    (*pairlist_elem)++;
-  }
-
   return func;
 }
 
index b1c0097..9f57de9 100644 (file)
@@ -1157,8 +1157,8 @@ If \texttt{group2} is a \texttt{dummyAtom}, this option is set to \texttt{yes} b
      \texttt{coordNum}}{%
      Pairlist control}{%
     decimal}{%
-    0.0}{This controls the pairlist feature, dictating the minimum value for each summation element in Eq.~\ref{eq:cvc_coordNum} such that the pair that contributed the summation element is included in subsequent simulation timesteps until the next pairlist recalculation. For most applications, this value should be small (eg. 0.001) to avoid missing important contributions to the overall sum. Higher values will improve performance, although values above 1 will exclude all possible pair interactions. Similarly, values below 0 will never exclude a pair from consideration.
-}
+    0.0}{This controls the pairlist feature, dictating the minimum value for each summation element in Eq.~\ref{eq:cvc_coordNum} such that the pair that contributed the summation element is included in subsequent simulation timesteps until the next pairlist recalculation. For most applications, this value should be small (eg. 0.001) to avoid missing important contributions to the overall sum. Higher values will improve performance by reducing the number of pairs that contribute to the sum. Values above 1 will exclude all possible pair interactions. Similarly, values below 0 will never exclude a pair from consideration. To ensure continuous forces, Eq.~\ref{eq:cvc_coordNum} is further modified by subtracting the tolerance and then rescaling so that each pair covers the range $\left[0, 1\right]$.
+  }
 
 \item %
     \labelkey{colvar|coordNum|pairListFrequency}
@@ -1167,8 +1167,8 @@ If \texttt{group2} is a \texttt{dummyAtom}, this option is set to \texttt{yes} b
      \texttt{coordNum}}{%
      Pairlist regeneration frequency}{%
     positive integer}{%
-    100}{This controls the pairlist feature, dictating how many steps are taken between regenerating pairlists if the tolerance is greater than 0. At this interval, the colvar defined will be exact, as though it was the all-to-all pair summation defined in Eq.~\ref{eq:cvc_coordNum}. All other steps will report only values and gradients from pairs in the pairlist.
-}
+    100}{This controls the pairlist feature, dictating how many steps are taken between regenerating pairlists if the tolerance is greater than 0.
+  }
 \end{cvcoptions}
 
 This component returns a dimensionless number, which ranges from
@@ -1177,6 +1177,7 @@ cutoff) to $N_{\mathtt{group1}} \times N_{\mathtt{group2}}$ (all distances
 are less than the cutoff), or $N_{\mathtt{group1}}$ if
 \texttt{group2CenterOnly} is used.  For performance reasons, at least
 one of \texttt{group1} and \texttt{group2} should be of limited size or \texttt{group2CenterOnly} should be used: the cost of the loop over all pairs grows as $N_{\mathtt{group1}} \times N_{\mathtt{group2}}$.
+Setting $\mathtt{tolerance} > 0$ ameliorates this to some degree, although every pair is still checked to regenerate the pairlist.
 
 
 
@@ -1672,7 +1673,6 @@ This angle is expressed in degrees within the range [0$^{\circ}$:180$^{\circ}$].
 \item %
   \dupkey{refPositionsColValue}{\texttt{orientationAngle}}{colvar|rmsd|refPositionsColValue}{\texttt{rmsd} component}
 }
-
 \end{cvcoptions}