Allow extended Lagrangian on non-scalar collective variables

lib/colvars/colvar.cpp

@@ -1085,7 +1085,7 @@ int colvar::calc_colvar_properties()
     // TODO: put it in the restart information
     if (cvm::step_relative() == 0) {
       xr = x;
-      vr = 0.0; // (already 0; added for clarity)
+      vr.reset(); // (already 0; added for clarity)
     }
 
     // report the restraint center as "value"
@@ -1171,7 +1171,8 @@ cvm::real colvar::update_forces_energy()
   if (is_enabled(f_cv_extended_Lagrangian)) {
 
     cvm::real dt = cvm::dt();
-    cvm::real f_ext;
+    colvarvalue f_ext(fr.type());
+    f_ext.reset();
 
     // the total force is applied to the fictitious mass, while the
     // atoms only feel the harmonic force
@@ -1200,8 +1201,10 @@ cvm::real colvar::update_forces_energy()
     potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
     // leap to v_(i+1/2)
     if (is_enabled(f_cv_Langevin)) {
-      vr -= dt * ext_gamma * vr.real_value;
-      vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
+      vr -= dt * ext_gamma * vr;
+      colvarvalue rnd(x);
+      rnd.set_random();
+      vr += dt * ext_sigma * rnd / ext_mass;
     }
     vr  += (0.5 * dt) * f_ext / ext_mass;
     xr  += dt * vr;

lib/colvars/colvarbias_meta.cpp

@@ -404,7 +404,7 @@ int colvarbias_meta::update()
 
     if (ebmeta) {
        hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
-       if(cvm::step_absolute() <= ebmeta_equil_steps) {
+       if(cvm::step_absolute() <= long(ebmeta_equil_steps)) {
          cvm::real const hills_lambda=(cvm::real(ebmeta_equil_steps - cvm::step_absolute()))/(cvm::real(ebmeta_equil_steps));
          hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
        }

lib/colvars/colvargrid.h

@@ -1071,7 +1071,7 @@ public:
   {
     // write the header
     os << "object 1 class gridpositions counts";
-    int icv;
+    size_t icv;
     for (icv = 0; icv < number_of_colvars(); icv++) {
       os << " " << number_of_points(icv);
     }
@@ -1166,8 +1166,8 @@ public:
 
   /// \brief Return the log-gradient from finite differences
   /// on the *same* grid for dimension n
-  inline const cvm::real log_gradient_finite_diff( const std::vector<int> &ix0,
-                                                   int n = 0)
+  inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
+                                            int n = 0)
   {
     cvm::real A0, A1;
     std::vector<int> ix;

lib/colvars/colvarmodule.h

@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-10-21"
+#define COLVARS_VERSION "2016-10-27"
 #endif
 
 #ifndef COLVARS_DEBUG
@@ -198,7 +198,7 @@ public:
   }
 
   /// \brief How many objects are configured yet?
-  inline size_t const size() const
+  inline size_t size() const
   {
     return colvars.size() + biases.size();
   }

lib/colvars/colvarproxy.h

@@ -25,7 +25,7 @@ public:
   colvarmodule *colvars;
 
   /// Default constructor
-  inline colvarproxy() : script(NULL), b_smp_active(true) {}
+  inline colvarproxy() : b_smp_active(true), script(NULL) {}
 
   /// Default destructor
   virtual ~colvarproxy() {}

lib/colvars/colvarvalue.cpp

@@ -70,6 +70,39 @@ void colvarvalue::set_elem(int const icv, colvarvalue const &x)
 }
 
 
+void colvarvalue::set_random()
+{
+  switch (this->type()) {
+  case colvarvalue::type_scalar:
+    this->real_value = cvm::rand_gaussian();
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    this->rvector_value.x = cvm::rand_gaussian();
+    this->rvector_value.y = cvm::rand_gaussian();
+    this->rvector_value.z = cvm::rand_gaussian();
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    this->quaternion_value.q0 = cvm::rand_gaussian();
+    this->quaternion_value.q1 = cvm::rand_gaussian();
+    this->quaternion_value.q2 = cvm::rand_gaussian();
+    this->quaternion_value.q3 = cvm::rand_gaussian();
+    break;
+  case colvarvalue::type_vector:
+    for (size_t ic = 0; ic < this->vector1d_value.size(); ic++) {
+      this->vector1d_value[ic] = cvm::rand_gaussian();
+    }
+    break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+    break;
+  }
+}
+
+
 colvarvalue colvarvalue::inverse() const
 {
   switch (value_type) {

lib/colvars/colvarvalue.h

@@ -297,6 +297,9 @@ public:
   /// 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;