Allow linear bias to deal with non-scalar variables (Colvars fix) 88/4688/2
authorGiacomo Fiorin <giacomo.fiorin@gmail.com>
Fri, 12 Oct 2018 17:49:47 +0000 (13:49 -0400)
committerJim Phillips <jim@ks.uiuc.edu>
Fri, 12 Oct 2018 20:15:52 +0000 (15:15 -0500)
See for details: https://github.com/Colvars/colvars/issues/181

Change-Id: Ib5656752b7472d4a7d1896013189ab100e97a31d

colvars/src/colvarbias_restraint.cpp
colvars/src/colvars_version.h
colvars/src/colvarvalue.cpp
colvars/src/colvarvalue.h

index cf11dc7..2daf7a0 100644 (file)
@@ -1284,13 +1284,16 @@ cvm::real colvarbias_restraint_linear::energy_difference(std::string const &conf
 
 cvm::real colvarbias_restraint_linear::restraint_potential(size_t i) const
 {
-  return force_k / variables(i)->width * (variables(i)->value() - colvar_centers[i]);
+  return force_k / variables(i)->width * (variables(i)->value() -
+                                          colvar_centers[i]).sum();
 }
 
 
 colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const
 {
-  return -1.0 * force_k / variables(i)->width;
+  colvarvalue dummy(variables(i)->value());
+  dummy.set_ones();
+  return -1.0 * force_k / variables(i)->width * dummy;
 }
 
 
index bedc7c8..bff0777 100644 (file)
@@ -1,5 +1,5 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2018-09-12"
+#define COLVARS_VERSION "2018-10-12"
 // 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 312d101..86b99ed 100644 (file)
@@ -379,6 +379,40 @@ void colvarvalue::set_random()
 }
 
 
+void colvarvalue::set_ones(cvm::real assigned_value)
+{
+  size_t ic;
+  switch (this->type()) {
+  case colvarvalue::type_scalar:
+    this->real_value = assigned_value;
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    this->rvector_value.x = assigned_value;
+    this->rvector_value.y = assigned_value;
+    this->rvector_value.z = assigned_value;
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    this->quaternion_value.q0 = assigned_value;
+    this->quaternion_value.q1 = assigned_value;
+    this->quaternion_value.q2 = assigned_value;
+    this->quaternion_value.q3 = assigned_value;
+    break;
+  case colvarvalue::type_vector:
+    for (ic = 0; ic < this->vector1d_value.size(); ic++) {
+      this->vector1d_value[ic] = assigned_value;
+    }
+    break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+    break;
+  }
+}
+
+
 void colvarvalue::undef_op() const
 {
   cvm::error("Error: Undefined operation on a colvar of type \""+
index 41759e9..25255e2 100644 (file)
@@ -187,6 +187,12 @@ public:
     return std::sqrt(this->norm2());
   }
 
+  /// Sum of the components of this colvarvalue (if more than one dimension)
+  cvm::real sum() const;
+
+  /// Return a colvarvalue object of the same type and all components set to 1
+  colvarvalue ones() const;
+
   /// Square distance between this \link colvarvalue \endlink and another
   cvm::real dist2(colvarvalue const &x2) const;
 
@@ -272,17 +278,21 @@ public:
   /// Get a single colvarvalue out of elements of the vector
   colvarvalue const get_elem(int const i_begin, int const i_end, Type const vt) const;
 
+  /// Get a single colvarvalue out of elements of the vector
+  colvarvalue const get_elem(int const icv) const;
+
+  /// Set elements of the vector from a single colvarvalue (uses the rank of x
+  /// to compute the length)
+  void set_elem(int const icv, colvarvalue const &x);
+
   /// Set elements of the vector from a single colvarvalue
   void set_elem(int const i_begin, int const i_end, colvarvalue const &x);
 
   /// Make each element a random number in N(0,1)
   void set_random();
 
-  /// Get a single colvarvalue out of elements of the vector
-  colvarvalue const get_elem(int const icv) const;
-
-  /// Set elements of the vector from a single colvarvalue
-  void set_elem(int const icv, colvarvalue const &x);
+  /// Make each element equal to the given argument
+  void set_ones(cvm::real assigned_value = 1.0);
 
   /// Get a scalar number out of an element of the vector
   cvm::real operator [] (int const i) const;
@@ -683,6 +693,29 @@ inline cvm::real colvarvalue::norm2() const
 }
 
 
+inline cvm::real colvarvalue::sum() const
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    return (this->real_value);
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return (this->rvector_value).x + (this->rvector_value).y +
+      (this->rvector_value).z;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return (this->quaternion_value).q0 + (this->quaternion_value).q1 +
+      (this->quaternion_value).q2 + (this->quaternion_value).q3;
+  case colvarvalue::type_vector:
+    return (this->vector1d_value).sum();
+  case colvarvalue::type_notset:
+  default:
+    return 0.0;
+  }
+}
+
+
 inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const
 {
   colvarvalue::check_types(*this, x2);