Update Colvars to version 2018-05-15 75/4175/2 master
authorGiacomo Fiorin <giacomo.fiorin@gmail.com>
Tue, 15 May 2018 15:08:31 +0000 (11:08 -0400)
committerDavid Hardy <dhardy@ks.uiuc.edu>
Tue, 15 May 2018 17:56:31 +0000 (12:56 -0500)
Introduces performance enhancements to coordNum and selfCoordNum, in
particular a pairlist feature introduced by Josh Vermaas.

Below are the relevant commits from https://github.com/Colvars/colvars/

8bee345 2018-05-15 Condense computation of value and gradients for coordNum and selfCoordNum [Giacomo Fiorin]
325e1b3 2018-05-15 Fix pairlist loop based on feedback from Josh [Giacomo Fiorin]
6faa836 2018-05-11 Unify coordnum template function, include pairlist check/update in it [Giacomo Fiorin]
13ac4a3 2018-05-10 Added documentation for the tolerance and pairListFrequency options [Josh Vermaas]
7564f23 2018-05-10 New check. Tolerance is MUCH higher than you'd want for real simulation, but it makes the difference clear quickly. Pairlistfrequency 2 means that every other printed value will match exactly with the non-tolerant coordnum [Josh Vermaas]
9394090 2018-05-10 Fixed compiler errors resulting from theorycraft-coding [Josh Vermaas]
020f555 2018-05-10 Incorparated Giacomo's comments. [Josh Vermaas]
1682ea7 2018-05-10 Prototype pairlist implementation [Josh Vermaas]

Change-Id: Ie27a2f59b44b1dfda269410b1745c0ed5ec5f068

12 files changed:
colvars/src/colvar_UIestimator.h
colvars/src/colvaratoms.cpp
colvars/src/colvaratoms.h
colvars/src/colvarcomp.cpp
colvars/src/colvarcomp.h
colvars/src/colvarcomp_coordnums.cpp
colvars/src/colvarmodule.cpp
colvars/src/colvarproxy.cpp
colvars/src/colvars_version.h
colvars/src/colvarscript.cpp
colvars/src/colvartypes.cpp
ug/ug_colvars.tex

index 36ed938..759b8d5 100644 (file)
@@ -731,6 +731,6 @@ namespace UIestimator {
             }
         }
     };
-};
+}
 
 #endif
index c8dcfc0..c7b6d3d 100644 (file)
@@ -21,6 +21,8 @@ cvm::atom::atom()
 {
   index = -1;
   id = -1;
+  mass = 1.0;
+  charge = 1.0;
   reset_data();
 }
 
index 93ae594..0b0dd62 100644 (file)
@@ -90,7 +90,7 @@ public:
   /// Destructor
   ~atom();
 
-  /// Set mutable data (everything except id and mass) to zero; update mass
+  /// Set mutable data (everything except id and mass) to zero
   inline void reset_data()
   {
     pos = cvm::atom_pos(0.0);
index c717e8f..9ecd650 100644 (file)
@@ -43,7 +43,6 @@ colvar::cvc::cvc(std::string const &conf)
 
 int colvar::cvc::init(std::string const &conf)
 {
-  int error_code = COLVARS_OK;
   if (cvm::debug())
     cvm::log("Initializing cvc base object.\n");
 
@@ -77,6 +76,8 @@ int colvar::cvc::init(std::string const &conf)
 
   if (cvm::debug())
     cvm::log("Done initializing cvc base object.\n");
+
+  return cvm::get_error();
 }
 
 
index b5c3b16..1d4bfe7 100644 (file)
@@ -834,49 +834,70 @@ protected:
   cvm::real     r0;
   /// \brief "Cutoff vector" for anisotropic calculation
   cvm::rvector  r0_vec;
-  /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
-  /// used
+  /// \brief Whether r/r0 or \vec{r}*\vec{1/r0_vec} should be used
   bool b_anisotropic;
   /// Integer exponent of the function numerator
   int en;
   /// Integer exponent of the function denominator
   int ed;
-  /// \brief If true, group2 will be treated as a single atom
-  /// (default: loop over all pairs of atoms in group1 and group2)
-  bool b_group2_center_only;
+
+  /// \brief If true, group2 will be treated as a single atom, stored in this
+  /// accessory group
+  cvm::atom_group *group2_center;
+
+  /// Tolerance for the pair list
+  cvm::real tolerance;
+
+  /// Frequency of update of the pair list
+  int pairlist_freq;
+
+  /// Pair list
+  bool *pairlist = NULL;
+
 public:
-  /// Constructor
+
   coordnum(std::string const &conf);
   coordnum();
-  virtual ~coordnum() {}
+  ~coordnum();
+
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
-  template<bool b_gradients>
-  /// \brief Calculate a coordination number through the function
-  /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
-  /// coordination number \param exp_num \i n exponent \param exp_den
-  /// \i m exponent \param A1 atom \param A2 atom
-  static cvm::real switching_function(cvm::real const &r0,
-                                      int const &exp_num, int const &exp_den,
-                                      cvm::atom &A1, cvm::atom &A2);
-
-  template<bool b_gradients>
-  /// \brief Calculate a coordination number through the function
-  /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
-  /// vector of different cutoffs in the three directions \param
-  /// exp_num \i n exponent \param exp_den \i m exponent \param A1
-  /// atom \param A2 atom
-  static cvm::real switching_function(cvm::rvector const &r0_vec,
-                                      int const &exp_num, int const &exp_den,
-                                      cvm::atom &A1, cvm::atom &A2);
-
   virtual cvm::real dist2(colvarvalue const &x1,
                           colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
+
+  enum {
+    ef_null = 0,
+    ef_gradients = 1,
+    ef_anisotropic = (1<<8),
+    ef_use_pairlist = (1<<9),
+    ef_rebuild_pairlist = (1<<10)
+  };
+
+  /// \brief Calculate a coordination number through the function
+  /// (1-x**n)/(1-x**m), where x = |A1-A2|/r0 \param r0, r0_vec "cutoff" for
+  /// the coordination number (scalar or vector depending on user choice)
+  /// \param en Numerator exponent \param ed Denominator exponent \param First
+  /// atom \param Second atom \param pairlist_elem pointer to pair flag for
+  /// this pair \param tolerance A pair is defined as having a larger
+  /// coordination than this number
+  template<int flags>
+  static cvm::real switching_function(cvm::real const &r0,
+                                      cvm::rvector const &r0_vec,
+                                      int en,
+                                      int ed,
+                                      cvm::atom &A1,
+                                      cvm::atom &A2,
+                                      bool **pairlist_elem,
+                                      cvm::real tolerance);
+
+  /// Main workhorse function
+  template<int flags> int compute_coordnum();
+
 };
 
 
@@ -887,7 +908,8 @@ class colvar::selfcoordnum
   : public colvar::cvc
 {
 protected:
-  /// First atom group
+
+  /// Selected atoms
   cvm::atom_group  *group1;
   /// \brief "Cutoff" for isotropic calculation (default)
   cvm::real     r0;
@@ -895,22 +917,18 @@ protected:
   int en;
   /// Integer exponent of the function denominator
   int ed;
+  cvm::real tolerance;
+  int pairlist_freq;
+  bool *pairlist = NULL;
+
 public:
-  /// Constructor
+
   selfcoordnum(std::string const &conf);
   selfcoordnum();
-  virtual ~selfcoordnum() {}
+  ~selfcoordnum();
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
-  template<bool b_gradients>
-  /// \brief Calculate a coordination number through the function
-  /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
-  /// coordination number \param exp_num \i n exponent \param exp_den
-  /// \i m exponent \param A1 atom \param A2 atom
-  static cvm::real switching_function(cvm::real const &r0,
-                                      int const &exp_num, int const &exp_den,
-                                      cvm::atom &A1, cvm::atom &A2);
 
   virtual cvm::real dist2(colvarvalue const &x1,
                           colvarvalue const &x2) const;
@@ -918,6 +936,9 @@ public:
                                   colvarvalue const &x2) const;
   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
+
+  /// Main workhorse function
+  template<int flags> int compute_selfcoordnum();
 };
 
 
@@ -947,26 +968,6 @@ public:
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force(colvarvalue const &force);
-  template<bool b_gradients>
-  /// \brief Calculate a coordination number through the function
-  /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
-  /// coordination number \param exp_num \i n exponent \param exp_den
-  /// \i m exponent \param A1 atom \param A2 atom
-  static cvm::real switching_function(cvm::real const &r0,
-                                      int const &exp_num, int const &exp_den,
-                                      cvm::atom &A1, cvm::atom &A2);
-
-  /*
-  template<bool b_gradients>
-  /// \brief Calculate a coordination number through the function
-  /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
-  /// vector of different cutoffs in the three directions \param
-  /// exp_num \i n exponent \param exp_den \i m exponent \param A1
-  /// atom \param A2 atom
-  static cvm::real switching_function(cvm::rvector const &r0_vec,
-                                      int const &exp_num, int const &exp_den,
-                                      cvm::atom &A1, cvm::atom &A2);
-  */
 
   virtual cvm::real dist2(colvarvalue const &x1,
                           colvarvalue const &x2) const;
index c34dc77..bcdb53c 100644 (file)
 
 
 
-
-template<bool calculate_gradients>
+template<int flags>
 cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
-                                               int const &en,
-                                               int const &ed,
+                                               cvm::rvector const &r0_vec,
+                                               int en,
+                                               int ed,
                                                cvm::atom &A1,
-                                               cvm::atom &A2)
+                                               cvm::atom &A2,
+                                               bool **pairlist_elem,
+                                               cvm::real pairlist_tol)
 {
-  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
-  cvm::real const l2 = diff.norm2()/(r0*r0);
-
-  // Assume en and ed are even integers, and avoid sqrt in the following
-  int const en2 = en/2;
-  int const ed2 = ed/2;
-
-  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);
-
-  if (calculate_gradients) {
-    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
-    cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
-    A1.grad += (-1.0)*dFdl2*dl2dx;
-    A2.grad +=        dFdl2*dl2dx;
+  if ((flags & ef_use_pairlist) && !(flags & ef_rebuild_pairlist)) {
+    bool const within = **pairlist_elem;
+    (*pairlist_elem)++;
+    if (!within) {
+      return 0.0;
+    }
   }
 
-  return func;
-}
-
+  cvm::rvector const r0sq_vec(r0_vec.x*r0_vec.x,
+                              r0_vec.y*r0_vec.y,
+                              r0_vec.z*r0_vec.z);
 
-template<bool calculate_gradients>
-cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
-                                               int const &en,
-                                               int const &ed,
-                                               cvm::atom &A1,
-                                               cvm::atom &A2)
-{
   cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
-  cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
+
+  cvm::rvector const scal_diff(diff.x/((flags & ef_anisotropic) ?
+                                       r0_vec.x : r0),
+                               diff.y/((flags & ef_anisotropic) ?
+                                       r0_vec.y : r0),
+                               diff.z/((flags & ef_anisotropic) ?
+                                       r0_vec.z : r0));
   cvm::real const l2 = scal_diff.norm2();
 
   // Assume en and ed are even integers, and avoid sqrt in the following
@@ -67,20 +58,31 @@ cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
   cvm::real const xd = cvm::integer_power(l2, ed2);
   cvm::real const func = (1.0-xn)/(1.0-xd);
 
-  if (calculate_gradients) {
-    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
-    cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
-                             (2.0/(r0_vec.y*r0_vec.y))*diff.y,
-                             (2.0/(r0_vec.z*r0_vec.z))*diff.z);
+  if (flags & ef_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) -
+                                            func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx((2.0/((flags & ef_anisotropic) ? r0sq_vec.x :
+                                   r0*r0)) * diff.x,
+                             (2.0/((flags & ef_anisotropic) ? r0sq_vec.y :
+                                   r0*r0)) * diff.y,
+                             (2.0/((flags & ef_anisotropic) ? r0sq_vec.z :
+                                   r0*r0)) * diff.z);
     A1.grad += (-1.0)*dFdl2*dl2dx;
     A2.grad +=        dFdl2*dl2dx;
   }
+
+  if (flags & ef_rebuild_pairlist) {
+    **pairlist_elem = (func > pairlist_tol) ? true : false;
+    (*pairlist_elem)++;
+  }
+
   return func;
 }
 
 
 colvar::coordnum::coordnum(std::string const &conf)
-  : cvc(conf), b_anisotropic(false), b_group2_center_only(false)
+  : cvc(conf), b_anisotropic(false), group2_center(NULL)
+
 {
   function_type = "coordnum";
   x.type(colvarvalue::type_scalar);
@@ -90,23 +92,26 @@ colvar::coordnum::coordnum(std::string const &conf)
 
   if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) {
     cvm::error("Error: group1 and group2 share a common atom (number: " +
-      cvm::to_str(atom_number) + ")\n");
+               cvm::to_str(atom_number) + ")\n", INPUT_ERROR);
     return;
   }
 
   if (group1->b_dummy) {
-    cvm::error("Error: only group2 is allowed to be a dummy atom\n");
+    cvm::error("Error: only group2 is allowed to be a dummy atom\n",
+               INPUT_ERROR);
     return;
   }
 
   bool const b_isotropic = get_keyval(conf, "cutoff", r0,
                                       cvm::real(4.0 * cvm::unit_angstrom()));
 
-  if (get_keyval(conf, "cutoff3", r0_vec, cvm::rvector(4.0 * cvm::unit_angstrom(),
-                                                       4.0 * cvm::unit_angstrom(),
-                                                       4.0 * cvm::unit_angstrom()))) {
+  if (get_keyval(conf, "cutoff3", r0_vec,
+                 cvm::rvector(4.0 * cvm::unit_angstrom(),
+                              4.0 * cvm::unit_angstrom(),
+                              4.0 * cvm::unit_angstrom()))) {
     if (b_isotropic) {
-      cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n",
+      cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" "
+                 "at the same time.\n",
                  INPUT_ERROR);
       return;
     }
@@ -135,86 +140,178 @@ colvar::coordnum::coordnum(std::string const &conf)
     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
   }
 
+  bool b_group2_center_only = false;
   get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
+  if (b_group2_center_only) {
+    if (!group2_center) {
+      group2_center = new cvm::atom_group();
+      group2_center->add_atom(cvm::atom());
+    }
+  }
+
+  get_keyval(conf, "tolerance", tolerance, 0.0);
+  if (tolerance > 0) {
+    get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
+    if ( ! (pairlist_freq > 0) ) {
+      cvm::error("Error: non-positive pairlistfrequency provided.\n",
+                 INPUT_ERROR);
+      return; // and do not allocate the pairlists below
+    }
+    if (b_group2_center_only) {
+      pairlist = new bool[group1->size()];
+    }
+    else {
+      pairlist = new bool[group1->size() * group2->size()];
+    }
+  }
+
 }
 
 
 colvar::coordnum::coordnum()
-  : b_anisotropic(false), b_group2_center_only(false)
+  : b_anisotropic(false), group2_center(NULL)
 {
   function_type = "coordnum";
   x.type(colvarvalue::type_scalar);
 }
 
 
-void colvar::coordnum::calc_value()
+colvar::coordnum::~coordnum()
 {
-  x.real_value = 0.0;
+  if (pairlist != NULL) {
+    delete [] pairlist;
+  }
+  if (group2_center != NULL) {
+    delete group2_center;
+  }
+}
 
-  if (b_group2_center_only) {
 
-    // create a fake atom to hold the group2 com coordinates
-    cvm::atom group2_com_atom;
-    group2_com_atom.pos = group2->center_of_mass();
+template<int compute_flags> int colvar::coordnum::compute_coordnum()
+{
+  if (group2_center) {
+    (*group2_center)[0].pos = group2->center_of_mass();
+    group2_center->calc_required_properties();
+  }
+  cvm::atom_group *group2p = group2_center ? group2_center : group2;
 
-    if (b_anisotropic) {
-      for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
-        x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, group2_com_atom);
-    } else {
-      for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
-        x.real_value += switching_function<false>(r0, en, ed, *ai1, group2_com_atom);
-    }
+  bool const use_pairlist = (pairlist != NULL);
+  bool const rebuild_pairlist = (pairlist != NULL) &&
+    (cvm::step_relative() % pairlist_freq == 0);
 
-  } else {
+  bool *pairlist_elem = use_pairlist ? pairlist : NULL;
+  cvm::atom_iter ai1 = group1->begin(), ai2 = group2p->begin();
+  cvm::atom_iter const ai1_end = group1->end();
+  cvm::atom_iter const ai2_end = group2p->end();
+
+  if (b_anisotropic) {
+
+    if (use_pairlist) {
+
+      if (rebuild_pairlist) {
 
-    if (b_anisotropic) {
-      for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
-        for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
-          x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, *ai2);
+        int const flags = compute_flags | ef_anisotropic | ef_use_pairlist |
+          ef_rebuild_pairlist;
+        for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
+          for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
+            x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
+                                                      *ai1, *ai2,
+                                                      &pairlist_elem,
+                                                      tolerance);
+          }
         }
-    } else {
-      for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
-        for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
-          x.real_value += switching_function<false>(r0, en, ed, *ai1, *ai2);
+
+      } else {
+
+        int const flags = compute_flags | ef_anisotropic | ef_use_pairlist;
+        for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
+          for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
+            x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
+                                                      *ai1, *ai2,
+                                                      &pairlist_elem,
+                                                      tolerance);
+          }
         }
-    }
-  }
-}
+      }
 
+    } else { // if (use_pairlist) {
 
-void colvar::coordnum::calc_gradients()
-{
-  if (b_group2_center_only) {
+      int const flags = compute_flags | ef_anisotropic;
+      for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
+        for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
+          x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
+                                                    *ai1, *ai2,
+                                                    NULL, 0.0);
+        }
+      }
+    }
 
-    // create a fake atom to hold the group2 com coordinates
-    cvm::atom group2_com_atom;
-    group2_com_atom.pos = group2->center_of_mass();
+  } else {
 
+    if (use_pairlist) {
 
-    if (b_anisotropic) {
-      for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
-        switching_function<true>(r0_vec, en, ed, *ai1, group2_com_atom);
-    } else {
-      for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
-        switching_function<true>(r0, en, ed, *ai1, group2_com_atom);
-    }
+      if (rebuild_pairlist) {
 
-    group2->set_weighted_gradient(group2_com_atom.grad);
+        int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
+        for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
+          for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
+            x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
+                                                      *ai1, *ai2,
+                                                      &pairlist_elem,
+                                                      tolerance);
+          }
+        }
 
-  } else {
+      } else {
 
-    if (b_anisotropic) {
-      for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
-        for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
-          switching_function<true>(r0_vec, en, ed, *ai1, *ai2);
+        int const flags = compute_flags | ef_use_pairlist;
+        for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
+          for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
+            x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
+                                                      *ai1, *ai2,
+                                                      &pairlist_elem,
+                                                      tolerance);
+          }
         }
-    } else {
-      for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
-        for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
-          switching_function<true>(r0, en, ed, *ai1, *ai2);
+      }
+
+    } else { // if (use_pairlist) {
+
+      int const flags = compute_flags;
+      for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
+        for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
+          x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
+                                                    *ai1, *ai2,
+                                                    NULL, 0.0);
         }
+      }
+    }
+  }
+
+  if (compute_flags & ef_gradients) {
+    if (group2_center) {
+      group2->set_weighted_gradient((*group2_center)[0].grad);
     }
   }
+
+  return COLVARS_OK;
+}
+
+
+void colvar::coordnum::calc_value()
+{
+  x.real_value = 0.0;
+  if (is_enabled(f_cvc_gradient)) {
+    compute_coordnum<ef_gradients>();
+  } else {
+    compute_coordnum<ef_null>();
+  }
+}
+
+
+void colvar::coordnum::calc_gradients()
+{
+  // Gradients are computed by calc_value() if f_cvc_gradients is enabled
 }
 
 
@@ -307,13 +404,24 @@ colvar::h_bond::~h_bond()
 
 void colvar::h_bond::calc_value()
 {
-  x.real_value = colvar::coordnum::switching_function<false>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
+  int const flags = coordnum::ef_null;
+  cvm::rvector const r0_vec(0.0); // TODO enable the flag?
+  x.real_value =
+    coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                        (*atom_groups[0])[0],
+                                        (*atom_groups[0])[1],
+                                        NULL, 0.0);
 }
 
 
 void colvar::h_bond::calc_gradients()
 {
-  colvar::coordnum::switching_function<true>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
+  int const flags = coordnum::ef_gradients;
+  cvm::rvector const r0_vec(0.0); // TODO enable the flag?
+  coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                      (*atom_groups[0])[0],
+                                      (*atom_groups[0])[1],
+                                      NULL, 0.0);
 }
 
 
@@ -353,6 +461,17 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf)
   if (!is_enabled(f_cvc_pbc_minimum_image)) {
     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
   }
+
+  get_keyval(conf, "tolerance", tolerance, 0.0);
+  if (tolerance > 0) {
+    get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
+    if ( ! (pairlist_freq > 0) ) {
+      cvm::error("Error: non-positive pairlistfrequency provided.\n",
+                 INPUT_ERROR);
+      return;
+    }
+    pairlist = new bool[(group1->size()-1) * (group1->size()-1)];
+  }
 }
 
 
@@ -363,26 +482,93 @@ colvar::selfcoordnum::selfcoordnum()
 }
 
 
+colvar::selfcoordnum::~selfcoordnum()
+{
+  if (pairlist != NULL) {
+    delete [] pairlist;
+  }
+}
+
+
+template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
+{
+  cvm::rvector const r0_vec(0.0); // TODO enable the flag?
+
+  bool const use_pairlist = (pairlist != NULL);
+  bool const rebuild_pairlist = (pairlist != NULL) &&
+    (cvm::step_relative() % pairlist_freq == 0);
+
+  bool *pairlist_elem = use_pairlist ? pairlist : NULL;
+  size_t i = 0, j = 0;
+  size_t const n = group1->size();
+
+  // Always isotropic (TODO: enable the ellipsoid?)
+
+  if (use_pairlist) {
+
+    if (rebuild_pairlist) {
+      int const flags = compute_flags | coordnum::ef_use_pairlist |
+        coordnum::ef_rebuild_pairlist;
+      for (i = 0; i < n - 1; i++) {
+        for (j = i + 1; j < n; j++) {
+          x.real_value +=
+            coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                                (*group1)[i],
+                                                (*group1)[j],
+                                                &pairlist_elem,
+                                                tolerance);
+        }
+      }
+    } else {
+      int const flags = compute_flags | coordnum::ef_use_pairlist;
+      for (i = 0; i < n - 1; i++) {
+        for (j = i + 1; j < n; j++) {
+          x.real_value +=
+            coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                                (*group1)[i],
+                                                (*group1)[j],
+                                                &pairlist_elem,
+                                                tolerance);
+        }
+      }
+    }
+
+  } else { // if (use_pairlist) {
+
+    int const flags = compute_flags | coordnum::ef_null;
+    for (i = 0; i < n - 1; i++) {
+      for (j = i + 1; j < n; j++) {
+        x.real_value +=
+          coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                              (*group1)[i],
+                                              (*group1)[j],
+                                              &pairlist_elem,
+                                              tolerance);
+      }
+    }
+  }
+
+  return COLVARS_OK;
+}
+
+
 void colvar::selfcoordnum::calc_value()
 {
   x.real_value = 0.0;
-  for (size_t i = 0; i < group1->size() - 1; i++) {
-    for (size_t j = i + 1; j < group1->size(); j++) {
-      x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, (*group1)[i], (*group1)[j]);
-    }
+  if (is_enabled(f_cvc_gradient)) {
+    compute_selfcoordnum<coordnum::ef_gradients>();
+  } else {
+    compute_selfcoordnum<coordnum::ef_null>();
   }
 }
 
 
 void colvar::selfcoordnum::calc_gradients()
 {
-  for (size_t i = 0; i < group1->size() - 1; i++) {
-    for (size_t j = i + 1; j < group1->size(); j++) {
-      colvar::coordnum::switching_function<true>(r0, en, ed, (*group1)[i], (*group1)[j]);
-    }
-  }
+  // Gradients are computed by calc_value() if f_cvc_gradients is enabled
 }
 
+
 void colvar::selfcoordnum::apply_force(colvarvalue const &force)
 {
   if (!group1->noforce) {
@@ -394,6 +580,7 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
 simple_scalar_dist_functions(selfcoordnum)
 
 
+
 // groupcoordnum member functions
 colvar::groupcoordnum::groupcoordnum(std::string const &conf)
   : distance(conf), b_anisotropic(false)
@@ -453,95 +640,56 @@ colvar::groupcoordnum::groupcoordnum()
 }
 
 
-template<bool calculate_gradients>
-cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
-                                                    int const &en,
-                                                    int const &ed,
-                                                    cvm::atom &A1,
-                                                    cvm::atom &A2)
-{
-  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
-  cvm::real const l2 = diff.norm2()/(r0*r0);
-
-  // Assume en and ed are even integers, and avoid sqrt in the following
-  int const en2 = en/2;
-  int const ed2 = ed/2;
-
-  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);
-
-  if (calculate_gradients) {
-    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
-    cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
-    A1.grad += (-1.0)*dFdl2*dl2dx;
-    A2.grad +=        dFdl2*dl2dx;
-  }
-
-  return func;
-}
-
-
-#if 0  // AMG: I don't think there's any reason to support anisotropic,
-       //      and I don't have those flags below in calc_value, but
-       //      if I need them, I'll also need to uncomment this method
-template<bool calculate_gradients>
-cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
-                                                    int const &en,
-                                                    int const &ed,
-                                                    cvm::atom &A1,
-                                                    cvm::atom &A2)
-{
-  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
-  cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
-  cvm::real const l2 = scal_diff.norm2();
-
-  // Assume en and ed are even integers, and avoid sqrt in the following
-  int const en2 = en/2;
-  int const ed2 = ed/2;
-
-  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);
-
-  if (calculate_gradients) {
-    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
-    cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
-                             (2.0/(r0_vec.y*r0_vec.y))*diff.y,
-                             (2.0/(r0_vec.z*r0_vec.z))*diff.z);
-    A1.grad += (-1.0)*dFdl2*dl2dx;
-    A2.grad +=        dFdl2*dl2dx;
-  }
-  return func;
-}
-#endif
-
-
 void colvar::groupcoordnum::calc_value()
 {
+  cvm::rvector const r0_vec(0.0); // TODO enable the flag?
 
   // create fake atoms to hold the com coordinates
   cvm::atom group1_com_atom;
   cvm::atom group2_com_atom;
   group1_com_atom.pos = group1->center_of_mass();
   group2_com_atom.pos = group2->center_of_mass();
-
-  x.real_value = coordnum::switching_function<false>(r0, en, ed,
-                                                     group1_com_atom, group2_com_atom);
+  if (b_anisotropic) {
+    int const flags = coordnum::ef_anisotropic;
+    x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                                       group1_com_atom,
+                                                       group2_com_atom,
+                                                       NULL, 0.0);
+  } else {
+    int const flags = coordnum::ef_null;
+    x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                                       group1_com_atom,
+                                                       group2_com_atom,
+                                                       NULL, 0.0);
+  }
 }
 
 
 void colvar::groupcoordnum::calc_gradients()
 {
+  cvm::rvector const r0_vec(0.0); // TODO enable the flag?
+
   cvm::atom group1_com_atom;
   cvm::atom group2_com_atom;
   group1_com_atom.pos = group1->center_of_mass();
   group2_com_atom.pos = group2->center_of_mass();
 
-  coordnum::switching_function<true>(r0, en, ed, group1_com_atom, group2_com_atom);
+  if (b_anisotropic) {
+    int const flags = coordnum::ef_gradients | coordnum::ef_anisotropic;
+    coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                        group1_com_atom,
+                                        group2_com_atom,
+                                        NULL, 0.0);
+  } else {
+    int const flags = coordnum::ef_gradients;
+    coordnum::switching_function<flags>(r0, r0_vec, en, ed,
+                                        group1_com_atom,
+                                        group2_com_atom,
+                                        NULL, 0.0);
+  }
+
   group1->set_weighted_gradient(group1_com_atom.grad);
   group2->set_weighted_gradient(group2_com_atom.grad);
-
 }
 
 
index c9e2d14..dc6242c 100644 (file)
@@ -8,7 +8,7 @@
 // Colvars repository at GitHub.
 
 #include <sstream>
-#include <string.h>
+#include <cstring>
 
 #include "colvarmodule.h"
 #include "colvarparse.h"
index c2992c3..86338df 100644 (file)
@@ -8,7 +8,7 @@
 // Colvars repository at GitHub.
 
 #include <sstream>
-#include <string.h>
+#include <cstring>
 
 #if defined(_OPENMP)
 #include <omp.h>
index 68f5cd1..731eb50 100644 (file)
@@ -1,5 +1,5 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2018-04-29"
+#define COLVARS_VERSION "2018-05-15"
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
index 0977496..a55e4c6 100644 (file)
@@ -8,8 +8,7 @@
 // Colvars repository at GitHub.
 
 #include <cstdlib>
-#include <stdlib.h>
-#include <string.h>
+#include <cstring>
 
 #define COLVARSCRIPT_CPP
 #include "colvarscript.h"
index b604606..2b45d77 100644 (file)
@@ -7,8 +7,8 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
-#include <stdlib.h>
-#include <string.h>
+#include <cstdlib>
+#include <cstring>
 
 #include "colvarmodule.h"
 #include "colvartypes.h"
index 6bb0ff1..7eeb4b9 100644 (file)
@@ -1828,6 +1828,26 @@ function is summed over all pairs of atoms in \texttt{group1} and
     If this option is \texttt{on}, only contacts between each atoms in \texttt{group1} and the center of mass of     \texttt{group2} are calculated (by default, the sum extends over all pairs of atoms in \texttt{group1} and \texttt{group2}).
 If \texttt{group2} is a \texttt{dummyAtom}, this option is set to \texttt{yes} by default.
 }
+
+\item %
+    \labelkey{colvar|coordNum|tolerance}
+    \keydef
+     {tolerance}{%
+     \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.
+}
+
+\item %
+    \labelkey{colvar|coordNum|pairListFrequency}
+    \keydef
+     {pairListFrequency}{%
+     \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.
+}
 \end{cvcoptions}
 
 This component returns a dimensionless number, which ranges from
@@ -1867,6 +1887,10 @@ those accepted by \texttt{coordNum}, namely \texttt{group1}
   \dupkey{expNumer}{\texttt{selfCoordNum}}{colvar|coordNum|expNumer}{\texttt{coordNum} component}
 \item %
   \dupkey{expDenom}{\texttt{selfCoordNum}}{colvar|coordNum|expDenom}{\texttt{coordNum} component}
+\item %
+  \dupkey{tolerance}{\texttt{selfCoordNum}}{colvar|coordNum|tolerance}{\texttt{coordNum} component}
+\item %
+  \dupkey{pairListFrequency}{\texttt{selfCoordNum}}{colvar|coordNum|pairListFrequency}{\texttt{coordNum} component}
 \end{cvcoptions}
 
 This component returns a dimensionless number, which ranges from