From 5890ef420b6476563f77a4c563c6a4cc22fd31f3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 28 Oct 2014 18:05:29 -0400 Subject: [PATCH] another day, another colvars update --- lib/colvars/colvar.cpp | 9 + lib/colvars/colvar.h | 7 + lib/colvars/colvarmodule.h | 2 +- lib/colvars/colvartypes.cpp | 15 +- lib/colvars/colvartypes.h | 267 +++++---- lib/colvars/colvarvalue.cpp | 597 ++++++++++++++++++++ lib/colvars/colvarvalue.h | 1060 ++++++++--------------------------- 7 files changed, 1036 insertions(+), 921 deletions(-) diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index e27b513ea3..fbbbf0ba0d 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -188,6 +188,15 @@ colvar::colvar(std::string const &conf) sorted_cvc_values.push_back(&(cvcs[j]->value())); } + // // these two are the vector value and the Jacobian matrix of the scripted function, respectively + // x_cvc.type(type_vector); + // dx_cvc.type(type_matrix); // TODO: not implemented yet + // for (j = 0; j < cvcs.size(); j++) { + // x_cvc.add_elem(cvcs[j]->value()); + // dx_cvc.add_elem(cvcs[j]->value()); + // } + + b_homogeneous = false; // Scripted functions are deemed non-periodic b_periodic = false; diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 87fe8d1a23..e144fea500 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -190,6 +190,13 @@ protected: /// Value of the colvar colvarvalue x; + // TODO: implement functionality to treat these + // /// Vector of individual values from CVCs + // colvarvalue x_cvc; + + // /// Jacobian matrix of individual values from CVCs + // colvarvalue dx_cvc; + /// Cached reported value (x may be manipulated) colvarvalue x_reported; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index cd74826e8c..7010eed2ac 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -4,7 +4,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2014-10-27" +#define COLVARS_VERSION "2014-10-28" #endif #ifndef COLVARS_DEBUG diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp index 9fa5eafe05..04b5c7957e 100644 --- a/lib/colvars/colvartypes.cpp +++ b/lib/colvars/colvartypes.cpp @@ -4,6 +4,7 @@ #include "colvartypes.h" #include "colvarparse.h" + std::string cvm::rvector::to_simple_string() const { std::ostringstream os; @@ -13,6 +14,7 @@ std::string cvm::rvector::to_simple_string() const return os.str(); } + int cvm::rvector::from_simple_string(std::string const &s) { std::stringstream stream(s); @@ -24,6 +26,7 @@ int cvm::rvector::from_simple_string(std::string const &s) return COLVARS_OK; } + std::ostream & operator << (std::ostream &os, colvarmodule::rvector const &v) { std::streamsize const w = os.width(); @@ -277,11 +280,14 @@ void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d &S, { // diagonalize int jac_nrot = 0; - jacobi(S, S_eigval, S_eigvec, &jac_nrot); - eigsrt(S_eigval, S_eigvec); + jacobi(S, S_eigval.c_array(), S_eigvec, &jac_nrot); + eigsrt(S_eigval.c_array(), S_eigvec); // jacobi saves eigenvectors by columns transpose(S_eigvec); + S_eigval.update_from_c_array(); + S_eigval.delete_c_array(); + // normalize eigenvectors for (size_t ie = 0; ie < 4; ie++) { cvm::real norm2 = 0.0; @@ -301,7 +307,7 @@ void colvarmodule::rotation::calc_optimal_rotation { cvm::matrix2d S(4, 4, 0.0); cvm::matrix2d S_backup(4, 4, 0.0); - cvm::vector1d S_eigval(4, 0.0); + cvm::vector1d S_eigval(4); cvm::matrix2d S_eigvec(4, 4, 0.0); build_matrix(pos1, pos2, S); @@ -470,7 +476,7 @@ void colvarmodule::rotation::calc_optimal_rotation if (b_debug_gradients) { cvm::matrix2d S_new(4, 4, 0.0); - cvm::vector1d S_new_eigval(4, 0.0); + cvm::vector1d S_new_eigval(4); cvm::matrix2d S_new_eigvec(4, 4, 0.0); // make an infitesimal move along each cartesian coordinate of @@ -519,6 +525,7 @@ void colvarmodule::rotation::calc_optimal_rotation h=a[k][l]; \ a[i][j]=g-s*(h+g*tau); \ a[k][l]=h+s*(g-h*tau); + #define n 4 void jacobi(cvm::real **a, cvm::real d[], cvm::real **v, int *nrot) diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h index d8c56f749e..01c8e0f9fc 100644 --- a/lib/colvars/colvartypes.h +++ b/lib/colvars/colvartypes.h @@ -18,145 +18,158 @@ /// \brief Arbitrary size array (one dimensions) suitable for linear /// algebra operations (i.e. for floating point numbers it can be used /// with library functions) -template class colvarmodule::vector1d +template class colvarmodule::vector1d : public std::vector { protected: - /// Underlying C-array + // member used to exchange with functions expecting a C-style array T *array; size_t length; public: - /// Allocation routine, used by all constructors - inline void alloc() { - if (length > 0) { - array = new T [length]; - } else { - array = NULL; - } - } - - /// Deallocation routine - inline void dealloc() { - if (array != NULL) { - delete [] array; - } - } - - /// Length of the array - inline size_t size() const - { - return length; - } - /// Default constructor - inline vector1d(size_t const n = 0, T const &t = T()) + inline vector1d(size_t const n = 0) : array (NULL) { length = n; - this->alloc(); + this->resize(n); reset(); } - /// Constructor from a 1-d C array - inline vector1d(size_t const n, T const *v) + /// Create a copy (unless it exists already), and return it + inline T * c_array() { - length = n; - this->alloc(); - for (size_t i = 0; i < length; i++) { - array[i] = v[i]; + if (array == NULL) { + array = new T[this->size()]; + for (size_t i = 0; i < this->size(); i++) { + array[i] = (*this)[i]; + } + } + return array; + } + + inline void update_from_c_array() + { + for (size_t i = 0; i < this->size(); i++) { + (*this)[i] = array[i]; } } - /// Copy constructor - inline vector1d(vector1d const &v) + inline void delete_c_array() { - length = v.length; - this->alloc(); - for (size_t i = 0; i < length; i++) { - array[i] = v.array[i]; - } + delete [] array; + array = NULL; + } + + inline ~vector1d() + { + this->delete_c_array(); } /// Set all elements to zero inline void reset() { - for (size_t i = 0; i < length; i++) { - array[i] = T(0.0); - } + this->assign(this->size(), T(0.0)); } - /// Assignment - inline vector1d & operator = (vector1d const &v) + inline static void check_sizes(size_t const n1, size_t const n2) { - if (length != v.length) { - this->dealloc(); - length = v.length; - this->alloc(); + if (n1 != n2) { + cvm::error("Error: trying to perform an operation between vectors of different sizes, "+ + cvm::to_str(n1)+" and "+cvm::to_str(n2)+".\n"); } - for (size_t i = 0; i < length; i++) { - this->array[i] = v.array[i]; + } + + inline void operator += (vector1d const &v) + { + check_sizes(this->size(), v.size()); + for (size_t i = 0; i < this->size(); i++) { + (*this)[i] += v[i]; } - return *this; } - /// Destructor - inline ~vector1d() { - this->dealloc(); + inline void operator -= (vector1d const &v) + { + check_sizes(this->size(), v.size()); + for (size_t i = 0; i < this->size(); i++) { + (*this)[i] -= v[i]; + } } - /// Return the 1-d C array - inline operator T *() { return array; } - - inline T & operator [] (int const &i) { - return array[i]; + inline void operator *= (cvm::real const &a) + { + for (size_t i = 0; i < this->size(); i++) { + (*this)[i] *= a; + } } - inline T const operator [] (int const &i) const { - return array[i]; + inline void operator /= (cvm::real const &a) + { + for (size_t i = 0; i < this->size(); i++) { + (*this)[i] /= a; + } + } + + inline friend vector1d operator + (vector1d const &v1, vector1d const &v2) + { + check_sizes(v1.size(), v2.size()); + vector1d result(v1.size()); + for (size_t i = 0; i < v1.size(); i++) { + result[i] = v1[i] + v2[i]; + } + return result; + } + + inline friend vector1d operator - (vector1d const &v1, vector1d const &v2) + { + check_sizes(v1.size(), v2.size()); + vector1d result(v1.size()); + for (size_t i = 0; i < v1.size(); i++) { + result[i] = v1[i] - v2[i]; + } + return result; + } + + inline friend vector1d operator * (vector1d const &v, cvm::real const &a) + { + vector1d result(v.size()); + for (size_t i = 0; i < v.size(); i++) { + result[i] = v[i] * a; + } + return result; + } + + inline friend vector1d operator * (cvm::real const &a, vector1d const &v) + { + return v * a; + } + + inline friend vector1d operator / (vector1d const &v, cvm::real const &a) + { + vector1d result(v.size()); + for (size_t i = 0; i < v.size(); i++) { + result[i] = v[i] / a; + } + return result; } /// Inner product - inline friend T operator * (vector1d const &v1, - vector1d const &v2) + inline friend T operator * (vector1d const &v1, vector1d const &v2) { + check_sizes(v1.size(), v2.size()); T prod(0.0); - size_t const min_length = (v1.length > v2.length) ? v1.length : v2.length; - for (size_t i = 0; i < min_length; i++) { - prod += v1.array[i] * v2.array[i]; + for (size_t i = 0; i < v1.size(); i++) { + prod += v1[i] * v2[i]; } return prod; } - /// Slicing - inline vector1d const slice(size_t const i1, size_t const i2) const - { - if ((i2 < i1) || (i1 < 0) || (i2 >= length)) { - cvm::error("Error: trying to slice a vector using incorrect boundaries.\n"); - } - vector1d result(i2 - i1, T()); - for (size_t i = 0; i < (i2 - i1); i++) { - result[i] = (*this)[i1+i]; - } - } - - /// Assign to a slice - inline void sliceassign(size_t const i1, size_t const i2, vector1d const &v) - { - if ((i2 < i1) || (i1 < 0) || (i2 >= length)) { - cvm::error("Error: trying to slice a vector using incorrect boundaries.\n"); - } - for (size_t i = 0; i < (i2 - i1); i++) { - (*this)[i1+i] = v[i]; - } - } - /// Squared norm inline cvm::real norm2() const { cvm::real result = 0.0; - for (size_t i = 0; i < length; i++) { - result += array[i] * array[i]; + for (size_t i = 0; i < this->size(); i++) { + result += (*this)[i] * (*this)[i]; } return result; } @@ -166,23 +179,75 @@ public: return std::sqrt(this->norm2()); } + /// Slicing + inline vector1d const slice(size_t const i1, size_t const i2) const + { + if ((i2 < i1) || (i1 < 0) || (i2 >= this->size())) { + cvm::error("Error: trying to slice a vector using incorrect boundaries.\n"); + } + vector1d result(i2 - i1); + for (size_t i = 0; i < (i2 - i1); i++) { + result[i] = (*this)[i1+i]; + } + return result; + } + + /// Assign a vector to a slice of this vector + inline void sliceassign(size_t const i1, size_t const i2, vector1d const &v) + { + if ((i2 < i1) || (i1 < 0) || (i2 >= this->size())) { + cvm::error("Error: trying to slice a vector using incorrect boundaries.\n"); + } + for (size_t i = 0; i < (i2 - i1); i++) { + (*this)[i1+i] = v[i]; + } + } + /// Formatted output - friend std::ostream & operator << (std::ostream &os, - vector1d const &v) + + inline friend std::ostream & operator << (std::ostream &os, cvm::vector1d const &v) { std::streamsize const w = os.width(); std::streamsize const p = os.precision(); os << "( "; - for (size_t i = 0; i < v.length-1; i++) { + for (size_t i = 0; i < v.size()-1; i++) { os.width(w); os.precision(p); - os << v.array[i] << " , "; + os << v[i] << " , "; } os.width(w); os.precision(p); - os << v.array[v.length-1] << " )"; + os << v[v.size()-1] << " )"; return os; } + + inline std::string to_simple_string() const + { + if (this->size() == 0) return std::string(""); + std::ostringstream os; + os.setf(std::ios::scientific, std::ios::floatfield); + os.precision(cvm::cv_prec); + os << (*this)[0]; + for (size_t i = 1; i < this->size(); i++) { + os << " " << (*this)[i]; + } + return os.str(); + } + + + inline int from_simple_string(std::string const &s) + { + std::stringstream stream(s); + size_t i = 0; + while ((stream >> (*this)[i]) && (i < this->size())) { + i++; + } + if (i < this->size()) { + return COLVARS_ERROR; + } + return COLVARS_OK; + } + }; @@ -242,10 +307,11 @@ public: inline cvm::vector1d const as_vector() const { - cvm::vector1d result(3, 0.0); + cvm::vector1d result(3); result[0] = this->x; result[1] = this->y; result[2] = this->z; + return result; } inline cvm::rvector & operator = (cvm::real const &v) @@ -762,11 +828,12 @@ public: inline cvm::vector1d const as_vector() const { - cvm::vector1d result(4, 0.0); + cvm::vector1d result(4); result[0] = q0; result[1] = q1; result[2] = q2; result[3] = q3; + return result; } /// Square norm of the quaternion @@ -1014,7 +1081,7 @@ public: { dS_1.resize(n, cvm::matrix2d(4, 4, cvm::rvector(0.0, 0.0, 0.0))); dL0_1.resize(n, cvm::rvector(0.0, 0.0, 0.0)); - dQ0_1.resize(n, cvm::vector1d(4, cvm::rvector(0.0, 0.0, 0.0))); + dQ0_1.resize(n, cvm::vector1d(4)); } /// Allocate space for the derivatives of the rotation @@ -1022,7 +1089,7 @@ public: { dS_2.resize(n, cvm::matrix2d(4, 4, cvm::rvector(0.0, 0.0, 0.0))); dL0_2.resize(n, cvm::rvector(0.0, 0.0, 0.0)); - dQ0_2.resize(n, cvm::vector1d(4, cvm::rvector(0.0, 0.0, 0.0))); + dQ0_2.resize(n, cvm::vector1d(4)); } /// \brief Calculate the optimal rotation and store the diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp index bcbe58b5fb..019ec36fd9 100644 --- a/lib/colvars/colvarvalue.cpp +++ b/lib/colvars/colvarvalue.cpp @@ -1,9 +1,606 @@ /// -*- c++ -*- #include +#include +#include #include "colvarmodule.h" #include "colvarvalue.h" +void colvarvalue::add_elem(colvarvalue const &x) +{ + if (this->value_type != type_vector) { + cvm::error("Error: trying to set an element for a variable that is not set to be a vector.\n"); + return; + } + size_t const n = vector1d_value.size(); + size_t const nd = num_dimensions(x.value_type); + elem_types.push_back(x.value_type); + elem_indices.push_back(n); + elem_sizes.push_back(nd); + vector1d_value.resize(n + nd); + set_elem(n, x); +} + + +colvarvalue const colvarvalue::get_elem(int const i_begin, int const i_end, Type const vt) const +{ + if (vector1d_value.size() > 0) { + cvm::vector1d const v(vector1d_value.slice(i_begin, i_end)); + return colvarvalue(v, vt); + } else { + cvm::error("Error: trying to get an element from a variable that is not a vector.\n"); + return colvarvalue(type_notset); + } +} + + +void colvarvalue::set_elem(int const i_begin, int const i_end, colvarvalue const &x) +{ + if (vector1d_value.size() > 0) { + vector1d_value.sliceassign(i_begin, i_end, x.as_vector()); + } else { + cvm::error("Error: trying to set an element for a variable that is not a vector.\n"); + } +} + + +colvarvalue const colvarvalue::get_elem(int const icv) const +{ + if (elem_types.size() > 0) { + return get_elem(elem_indices[icv], elem_indices[icv] + elem_sizes[icv], + elem_types[icv]); + } else { + cvm::error("Error: trying to get a colvarvalue element from a vector colvarvalue that was initialized as a plain array.\n"); + return colvarvalue(type_notset); + } +} + + +void colvarvalue::set_elem(int const icv, colvarvalue const &x) +{ + if (elem_types.size() > 0) { + check_types_assign(elem_types[icv], x.value_type); + set_elem(elem_indices[icv], elem_indices[icv] + elem_sizes[icv], x); + } else { + cvm::error("Error: trying to set a colvarvalue element for a colvarvalue that was initialized as a plain array.\n"); + } +} + + +colvarvalue colvarvalue::inverse() const +{ + switch (value_type) { + case colvarvalue::type_scalar: + return colvarvalue(1.0/real_value); + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return colvarvalue(cvm::rvector(1.0/rvector_value.x, + 1.0/rvector_value.y, + 1.0/rvector_value.z)); + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return colvarvalue(cvm::quaternion(1.0/quaternion_value.q0, + 1.0/quaternion_value.q1, + 1.0/quaternion_value.q2, + 1.0/quaternion_value.q3)); + break; + case colvarvalue::type_vector: + { + cvm::vector1d result(vector1d_value); + if (elem_types.size() > 0) { + // if we have information about non-scalar types, use it + for (size_t i = 0; i < elem_types.size(); i++) { + result.sliceassign(elem_indices[i], elem_indices[i]+elem_sizes[i], + cvm::vector1d((this->get_elem(i)).inverse())); + } + } else { + for (size_t i = 0; i < result.size(); i++) { + if (result[i] != 0.0) { + result = 1.0/result[i]; + } + } + } + return colvarvalue(result, type_vector); + } + break; + case colvarvalue::type_notset: + default: + undef_op(); + break; + } + return colvarvalue(); +} + + +// binary operations between two colvarvalues + +colvarvalue operator + (colvarvalue const &x1, + colvarvalue const &x2) +{ + colvarvalue::check_types(x1, x2); + + switch (x1.value_type) { + case colvarvalue::type_scalar: + return colvarvalue(x1.real_value + x2.real_value); + case colvarvalue::type_3vector: + return colvarvalue(x1.rvector_value + x2.rvector_value); + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return colvarvalue(x1.rvector_value + x2.rvector_value, + colvarvalue::type_unit3vector); + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return colvarvalue(x1.quaternion_value + x2.quaternion_value); + case colvarvalue::type_vector: + return colvarvalue(x1.vector1d_value + x2.vector1d_value, colvarvalue::type_vector); + case colvarvalue::type_notset: + default: + x1.undef_op(); + return colvarvalue(colvarvalue::type_notset); + }; +} + + +colvarvalue operator - (colvarvalue const &x1, + colvarvalue const &x2) +{ + colvarvalue::check_types(x1, x2); + + switch (x1.value_type) { + case colvarvalue::type_scalar: + return colvarvalue(x1.real_value - x2.real_value); + case colvarvalue::type_3vector: + return colvarvalue(x1.rvector_value - x2.rvector_value); + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return colvarvalue(x1.rvector_value - x2.rvector_value, + colvarvalue::type_unit3vector); + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return colvarvalue(x1.quaternion_value - x2.quaternion_value); + case colvarvalue::type_vector: + return colvarvalue(x1.vector1d_value - x2.vector1d_value, colvarvalue::type_vector); + case colvarvalue::type_notset: + default: + x1.undef_op(); + return colvarvalue(colvarvalue::type_notset); + }; +} + + +// binary operations with real numbers + +colvarvalue operator * (cvm::real const &a, + colvarvalue const &x) +{ + switch (x.value_type) { + case colvarvalue::type_scalar: + return colvarvalue(a * x.real_value); + case colvarvalue::type_3vector: + return colvarvalue(a * x.rvector_value); + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return colvarvalue(a * x.rvector_value, + colvarvalue::type_unit3vector); + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return colvarvalue(a * x.quaternion_value); + case colvarvalue::type_vector: + return colvarvalue(x.vector1d_value * a, colvarvalue::type_vector); + case colvarvalue::type_notset: + default: + x.undef_op(); + return colvarvalue(colvarvalue::type_notset); + } +} + + +colvarvalue operator * (colvarvalue const &x, + cvm::real const &a) +{ + return a * x; +} + + +colvarvalue operator / (colvarvalue const &x, + cvm::real const &a) +{ + switch (x.value_type) { + case colvarvalue::type_scalar: + return colvarvalue(x.real_value / a); + case colvarvalue::type_3vector: + return colvarvalue(x.rvector_value / a); + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return colvarvalue(x.rvector_value / a, + colvarvalue::type_unit3vector); + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return colvarvalue(x.quaternion_value / a); + case colvarvalue::type_vector: + return colvarvalue(x.vector1d_value / a, colvarvalue::type_vector); + case colvarvalue::type_notset: + default: + x.undef_op(); + return colvarvalue(colvarvalue::type_notset); + } +} + + +// inner product between two colvarvalues + +cvm::real operator * (colvarvalue const &x1, + colvarvalue const &x2) +{ + colvarvalue::check_types(x1, x2); + + switch (x1.value_type) { + case colvarvalue::type_scalar: + return (x1.real_value * x2.real_value); + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return (x1.rvector_value * x2.rvector_value); + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + // the "*" product is the quaternion product, here the inner + // member function is used instead + return (x1.quaternion_value.inner(x2.quaternion_value)); + case colvarvalue::type_vector: + return (x1.vector1d_value * x2.vector1d_value); + case colvarvalue::type_notset: + default: + x1.undef_op(); + return 0.0; + }; +} + + +colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const +{ + colvarvalue::check_types(*this, x2); + + switch (this->value_type) { + case colvarvalue::type_scalar: + return 2.0 * (this->real_value - x2.real_value); + case colvarvalue::type_3vector: + return 2.0 * (this->rvector_value - x2.rvector_value); + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + { + cvm::rvector const &v1 = this->rvector_value; + cvm::rvector const &v2 = x2.rvector_value; + cvm::real const cos_t = v1 * v2; + cvm::real const sin_t = std::sqrt(1.0 - cos_t*cos_t); + return colvarvalue( 2.0 * sin_t * + cvm::rvector((-1.0) * sin_t * v2.x + + cos_t/sin_t * (v1.x - cos_t*v2.x), + (-1.0) * sin_t * v2.y + + cos_t/sin_t * (v1.y - cos_t*v2.y), + (-1.0) * sin_t * v2.z + + cos_t/sin_t * (v1.z - cos_t*v2.z) + ), + colvarvalue::type_unit3vector ); + } + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return this->quaternion_value.dist2_grad(x2.quaternion_value); + case colvarvalue::type_vector: + return colvarvalue(2.0 * (this->vector1d_value - x2.vector1d_value), colvarvalue::type_vector); + break; + case colvarvalue::type_notset: + default: + this->undef_op(); + return colvarvalue(colvarvalue::type_notset); + }; +} + + +std::string colvarvalue::to_simple_string() const +{ + switch (type()) { + case colvarvalue::type_scalar: + return cvm::to_str(real_value, 0, cvm::cv_prec); + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return rvector_value.to_simple_string(); + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return quaternion_value.to_simple_string(); + break; + case colvarvalue::type_vector: + return vector1d_value.to_simple_string(); + break; + case colvarvalue::type_notset: + default: + undef_op(); + break; + } + return std::string(); +} + + +int colvarvalue::from_simple_string(std::string const &s) +{ + switch (type()) { + case colvarvalue::type_scalar: + return ((std::istringstream(s) >> real_value) + ? COLVARS_OK : COLVARS_ERROR); + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return rvector_value.from_simple_string(s); + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return quaternion_value.from_simple_string(s); + break; + case colvarvalue::type_vector: + return vector1d_value.from_simple_string(s); + break; + case colvarvalue::type_notset: + default: + undef_op(); + break; + } + return COLVARS_ERROR; +} + +std::ostream & operator << (std::ostream &os, colvarvalue const &x) +{ + switch (x.type()) { + case colvarvalue::type_scalar: + os << x.real_value; + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + os << x.rvector_value; + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + os << x.quaternion_value; + break; + case colvarvalue::type_notset: + os << "not set"; break; + } + return os; +} + + +std::ostream & operator << (std::ostream &os, std::vector const &v) +{ + for (size_t i = 0; i < v.size(); i++) { + os << v[i]; + } + return os; +} + + +std::istream & operator >> (std::istream &is, colvarvalue &x) +{ + if (x.type() == colvarvalue::type_notset) { + cvm::error("Trying to read from a stream a colvarvalue, " + "which has not yet been assigned a data type.\n"); + return is; + } + + switch (x.type()) { + case colvarvalue::type_scalar: + is >> x.real_value; + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vectorderiv: + is >> x.rvector_value; + break; + case colvarvalue::type_unit3vector: + is >> x.rvector_value; + x.apply_constraints(); + break; + case colvarvalue::type_quaternion: + is >> x.quaternion_value; + x.apply_constraints(); + break; + case colvarvalue::type_quaternionderiv: + is >> x.quaternion_value; + break; + default: + x.undef_op(); + } + return is; +} + + +size_t colvarvalue::output_width(size_t const &real_width) const +{ + switch (this->value_type) { + case colvarvalue::type_scalar: + return real_width; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return cvm::rvector::output_width(real_width); + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return cvm::quaternion::output_width(real_width); + case colvarvalue::type_notset: + default: + return 0; + } +} + + +void colvarvalue::inner_opt(colvarvalue const &x, + std::vector::iterator &xv, + std::vector::iterator const &xv_end, + std::vector::iterator &result) +{ + // doing type check only once, here + colvarvalue::check_types(x, *xv); + + std::vector::iterator &xvi = xv; + std::vector::iterator &ii = result; + + switch (x.value_type) { + case colvarvalue::type_scalar: + while (xvi != xv_end) { + *(ii++) += (xvi++)->real_value * x.real_value; + } + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + while (xvi != xv_end) { + *(ii++) += (xvi++)->rvector_value * x.rvector_value; + } + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + while (xvi != xv_end) { + *(ii++) += ((xvi++)->quaternion_value).cosine(x.quaternion_value); + } + break; + default: + x.undef_op(); + }; +} + + +void colvarvalue::inner_opt(colvarvalue const &x, + std::list::iterator &xv, + std::list::iterator const &xv_end, + std::vector::iterator &result) +{ + // doing type check only once, here + colvarvalue::check_types(x, *xv); + + std::list::iterator &xvi = xv; + std::vector::iterator &ii = result; + + switch (x.value_type) { + case colvarvalue::type_scalar: + while (xvi != xv_end) { + *(ii++) += (xvi++)->real_value * x.real_value; + } + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + while (xvi != xv_end) { + *(ii++) += (xvi++)->rvector_value * x.rvector_value; + } + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + while (xvi != xv_end) { + *(ii++) += ((xvi++)->quaternion_value).cosine(x.quaternion_value); + } + break; + default: + x.undef_op(); + }; +} + + +void colvarvalue::p2leg_opt(colvarvalue const &x, + std::vector::iterator &xv, + std::vector::iterator const &xv_end, + std::vector::iterator &result) +{ + // doing type check only once, here + colvarvalue::check_types(x, *xv); + + std::vector::iterator &xvi = xv; + std::vector::iterator &ii = result; + + switch (x.value_type) { + case colvarvalue::type_scalar: + cvm::error("Error: cannot calculate Legendre polynomials " + "for scalar variables.\n"); + return; + break; + case colvarvalue::type_3vector: + while (xvi != xv_end) { + cvm::real const cosine = + ((xvi)->rvector_value * x.rvector_value) / + ((xvi)->rvector_value.norm() * x.rvector_value.norm()); + xvi++; + *(ii++) += 1.5*cosine*cosine - 0.5; + } + break; + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + while (xvi != xv_end) { + cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value; + *(ii++) += 1.5*cosine*cosine - 0.5; + } + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + while (xvi != xv_end) { + cvm::real const cosine = (xvi++)->quaternion_value.cosine(x.quaternion_value); + *(ii++) += 1.5*cosine*cosine - 0.5; + } + break; + default: + x.undef_op(); + }; +} + + +void colvarvalue::p2leg_opt(colvarvalue const &x, + std::list::iterator &xv, + std::list::iterator const &xv_end, + std::vector::iterator &result) +{ + // doing type check only once, here + colvarvalue::check_types(x, *xv); + + std::list::iterator &xvi = xv; + std::vector::iterator &ii = result; + + switch (x.value_type) { + case colvarvalue::type_scalar: + cvm::error("Error: cannot calculate Legendre polynomials " + "for scalar variables.\n"); + break; + case colvarvalue::type_3vector: + while (xvi != xv_end) { + cvm::real const cosine = + ((xvi)->rvector_value * x.rvector_value) / + ((xvi)->rvector_value.norm() * x.rvector_value.norm()); + xvi++; + *(ii++) += 1.5*cosine*cosine - 0.5; + } + break; + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + while (xvi != xv_end) { + cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value; + *(ii++) += 1.5*cosine*cosine - 0.5; + } + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + while (xvi != xv_end) { + cvm::real const cosine = (xvi++)->quaternion_value.cosine(x.quaternion_value); + *(ii++) += 1.5*cosine*cosine - 0.5; + } + break; + default: + x.undef_op(); + }; +} + + + diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h index 03ddfdc838..24cd35c7f8 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -109,7 +109,7 @@ public: /// \brief Default constructor: this class defaults to a scalar /// number and always behaves like it unless you change its type inline colvarvalue() - : real_value(0.0), value_type(type_scalar) + : value_type(type_scalar), real_value(0.0) {} /// Constructor from a type specification @@ -121,26 +121,26 @@ public: /// Copy constructor from real base type inline colvarvalue(cvm::real const &x) - : real_value(x), value_type(type_scalar) + : value_type(type_scalar), real_value(x) {} /// \brief Copy constructor from rvector base type (Note: this sets /// automatically a type \link type_3vector \endlink , if you want a /// \link type_unit3vector \endlink you must set it explicitly) inline colvarvalue(cvm::rvector const &v) - : rvector_value(v), value_type(type_3vector) + : value_type(type_3vector), rvector_value(v) {} /// \brief Copy constructor from rvector base type (additional /// argument to make possible to choose a \link type_unit3vector /// \endlink inline colvarvalue(cvm::rvector const &v, Type const &vti) - : rvector_value(v), value_type(vti) + : value_type(vti), rvector_value(v) {} /// \brief Copy constructor from quaternion base type inline colvarvalue(cvm::quaternion const &q) - : quaternion_value(q), value_type(type_quaternion) + : value_type(type_quaternion), quaternion_value(q) {} /// Copy constructor from another \link colvarvalue \endlink @@ -168,16 +168,20 @@ public: /// Set the type explicitly inline void type(Type const &vti) { + // reset the value based on the previous type reset(); value_type = vti; + // reset the value based on the new type reset(); } /// Set the type after another \link colvarvalue \endlink inline void type(colvarvalue const &x) { + // reset the value based on the previous type reset(); value_type = x.value_type; + // reset the value based on the new type reset(); } @@ -212,6 +216,15 @@ public: void operator *= (cvm::real const &a); void operator /= (cvm::real const &a); + // Binary operators (return values) + friend colvarvalue operator + (colvarvalue const &x1, colvarvalue const &x2); + friend colvarvalue operator - (colvarvalue const &x1, colvarvalue const &x2); + friend colvarvalue operator * (colvarvalue const &x, cvm::real const &a); + friend colvarvalue operator * (cvm::real const &a, colvarvalue const &x); + friend colvarvalue operator / (colvarvalue const &x, cvm::real const &a); + + /// Inner product + friend cvm::real operator * (colvarvalue const &x1, colvarvalue const &x2); // Cast to scalar inline operator cvm::real() const @@ -249,7 +262,7 @@ public: return quaternion_value; } - // Cast to n-dimensional vector + // Create a n-dimensional vector from one of the basic types, or return the existing vector cvm::vector1d const as_vector() const; @@ -259,6 +272,12 @@ public: return (value_type == type_scalar); } + + /// Add an element to the vector (requires that type_vector is already set). + /// This is only needed to use this object as a vector of "complex" colvar values. + /// To use it instead as a plain n-dimensional vector, access vector1d_value directly. + void add_elem(colvarvalue const &x); + /// 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; @@ -304,6 +323,13 @@ public: /// Undefined operation void undef_op() const; + + /// \brief Formatted output operator + friend std::ostream & operator << (std::ostream &os, colvarvalue const &q); + + /// \brief Formatted input operator + friend std::istream & operator >> (std::istream &is, colvarvalue &q); + /// Give the number of characters required to output this /// colvarvalue, given the current type assigned and the number of /// characters for a real number @@ -315,23 +341,23 @@ public: /// Parses value from a script-friendly string (space separated list) int from_simple_string(std::string const &s); - // optimized routines for operations with an array; xv and inner as - // vectors are assumed to have the same number of elements (i.e. the - // end for inner comes together with xv_end in the loop) + + // optimized routines for operations on arrays of colvar values; + // xv and result are assumed to have the same number of elements /// \brief Optimized routine for the inner product of one collective /// variable with an array void static inner_opt(colvarvalue const &x, std::vector::iterator &xv, std::vector::iterator const &xv_end, - std::vector::iterator &inner); + std::vector::iterator &result); /// \brief Optimized routine for the inner product of one collective /// variable with an array void static inner_opt(colvarvalue const &x, std::list::iterator &xv, std::list::iterator const &xv_end, - std::vector::iterator &inner); + std::vector::iterator &result); /// \brief Optimized routine for the second order Legendre /// polynomial, (3cos^2(w)-1)/2, of one collective variable with an @@ -339,14 +365,14 @@ public: void static p2leg_opt(colvarvalue const &x, std::vector::iterator &xv, std::vector::iterator const &xv_end, - std::vector::iterator &inner); + std::vector::iterator &result); /// \brief Optimized routine for the second order Legendre /// polynomial of one collective variable with an array void static p2leg_opt(colvarvalue const &x, std::list::iterator &xv, std::list::iterator const &xv_end, - std::vector::iterator &inner); + std::vector::iterator &result); }; @@ -466,10 +492,10 @@ inline size_t colvarvalue::size() const } - inline colvarvalue::colvarvalue(colvarvalue const &x) : value_type(x.value_type) { + // reset the value based on the previous type reset(); switch (x.value_type) { @@ -498,6 +524,7 @@ inline colvarvalue::colvarvalue(colvarvalue const &x) inline colvarvalue::colvarvalue(cvm::vector1d const &v, Type const &vti) { + // reset the value based on the previous type reset(); if (v.size() != num_dimensions(value_type)) { cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+ @@ -530,48 +557,222 @@ inline colvarvalue::colvarvalue(cvm::vector1d const &v, Type const &v } -inline colvarvalue const colvarvalue::get_elem(int const i_begin, int const i_end, Type const vt) const +inline void colvarvalue::check_types(colvarvalue const &x1, + colvarvalue const &x2) { - if (vector1d_value.size() > 0) { - cvm::vector1d const v(vector1d_value.slice(i_begin, i_end)); - return colvarvalue(v, vt); - } else { - cvm::error("Error: trying to get an element from a variable that is not a vector.\n"); - return colvarvalue(type_notset); - } -} + if (!colvarvalue::type_checking()) return; -inline void colvarvalue::set_elem(int const i_begin, int const i_end, colvarvalue const &x) -{ - if (vector1d_value.size() > 0) { - vector1d_value.sliceassign(i_begin, i_end, x.as_vector()); - } else { - cvm::error("Error: trying to set an element for a variable that is not a vector.\n"); + if (x1.value_type != x2.value_type) { + if (((x1.value_type == type_unit3vector) && + (x2.value_type == type_unit3vectorderiv)) || + ((x2.value_type == type_unit3vector) && + (x1.value_type == type_unit3vectorderiv)) || + ((x1.value_type == type_quaternion) && + (x2.value_type == type_quaternionderiv)) || + ((x2.value_type == type_quaternion) && + (x1.value_type == type_quaternionderiv))) { + return; + } + cvm::error("Performing an operation between two colvar " + "values with different types, \""+ + colvarvalue::type_desc(x1.value_type)+ + "\" and \""+ + colvarvalue::type_desc(x2.value_type)+ + "\".\n"); + } + + if (x1.value_type == type_vector) { + if (x1.vector1d_value.size() != x2.vector1d_value.size()) { + cvm::error("Performing an operation between two vector colvar " + "values with different sizes, "+ + cvm::to_str(x1.vector1d_value.size())+ + " and "+ + cvm::to_str(x2.vector1d_value.size())+ + ".\n"); + } } } -inline colvarvalue const colvarvalue::get_elem(int const icv) const +inline void colvarvalue::check_types_assign(colvarvalue::Type const &vt1, + colvarvalue::Type const &vt2) { - if (elem_types.size() > 0) { - return get_elem(elem_indices[icv], elem_indices[icv] + elem_sizes[icv], - elem_types[icv]); - } else { - cvm::error("Error: trying to get a colvarvalue element from a vector colvarvalue that was initialized as a plain array.\n"); - return colvarvalue(type_notset); + if (vt1 != type_notset) { + if (vt1 != vt2) { + cvm::error("Trying to assign a colvar value with type \""+ + type_desc(vt2)+"\" to one with type \""+ + type_desc(vt1)+"\".\n"); + } } } -inline void colvarvalue::set_elem(int const icv, colvarvalue const &x) + +inline void colvarvalue::undef_op() const { - if (elem_types.size() > 0) { - check_types_assign(elem_types[icv], x.value_type); - set_elem(elem_indices[icv], elem_indices[icv] + elem_sizes[icv], x); - } else { - cvm::error("Error: trying to set a colvarvalue element for a colvarvalue that was initialized as a plain array.\n"); + cvm::error("Error: Undefined operation on a colvar of type \""+ + type_desc(this->value_type)+"\".\n"); +} + + +inline colvarvalue & colvarvalue::operator = (colvarvalue const &x) +{ + check_types_assign(this->value_type, x.value_type); + value_type = x.value_type; + + switch (this->value_type) { + case colvarvalue::type_scalar: + this->real_value = x.real_value; + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + this->rvector_value = x.rvector_value; + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + this->quaternion_value = x.quaternion_value; + break; + case colvarvalue::type_vector: + vector1d_value = x.vector1d_value; + elem_types = x.elem_types; + elem_indices = x.elem_indices; + elem_sizes = x.elem_sizes; + break; + case colvarvalue::type_notset: + default: + undef_op(); + break; + } + return *this; +} + + +inline void colvarvalue::operator += (colvarvalue const &x) +{ + colvarvalue::check_types(*this, x); + + switch (this->value_type) { + case colvarvalue::type_scalar: + this->real_value += x.real_value; + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + this->rvector_value += x.rvector_value; + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + this->quaternion_value += x.quaternion_value; + break; + case colvarvalue::type_vector: + this->vector1d_value += x.vector1d_value; + break; + case colvarvalue::type_notset: + default: + undef_op(); } } +inline void colvarvalue::operator -= (colvarvalue const &x) +{ + colvarvalue::check_types(*this, x); + + switch (value_type) { + case colvarvalue::type_scalar: + real_value -= x.real_value; + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + rvector_value -= x.rvector_value; + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + quaternion_value -= x.quaternion_value; + break; + case colvarvalue::type_vector: + this->vector1d_value -= x.vector1d_value; + break; + case colvarvalue::type_notset: + default: + undef_op(); + } +} + + +inline void colvarvalue::operator *= (cvm::real const &a) +{ + switch (value_type) { + case colvarvalue::type_scalar: + real_value *= a; + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vectorderiv: + rvector_value *= a; + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + quaternion_value *= a; + break; + case colvarvalue::type_vector: + this->vector1d_value *= a; + break; + case colvarvalue::type_notset: + default: + undef_op(); + } +} + + +inline void colvarvalue::operator /= (cvm::real const &a) +{ + switch (value_type) { + case colvarvalue::type_scalar: + real_value /= a; break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + rvector_value /= a; break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + quaternion_value /= a; break; + case colvarvalue::type_vector: + this->vector1d_value /= a; + break; + case colvarvalue::type_notset: + default: + undef_op(); + } +} + + +inline cvm::vector1d const colvarvalue::as_vector() const +{ + switch (value_type) { + case colvarvalue::type_scalar: + { + cvm::vector1d v(1); + v = real_value; + return v; + } + break; + case colvarvalue::type_3vector: + case colvarvalue::type_unit3vector: + case colvarvalue::type_unit3vectorderiv: + return rvector_value.as_vector(); + break; + case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: + return quaternion_value.as_vector(); + break; + case colvarvalue::type_vector: + return vector1d_value; + break; + case colvarvalue::type_notset: + default: + break; + } +} inline void colvarvalue::reset() @@ -656,7 +857,6 @@ inline void colvarvalue::is_derivative() } - inline cvm::real colvarvalue::norm2() const { switch (value_type) { @@ -688,362 +888,6 @@ inline cvm::real colvarvalue::norm2() const } -inline colvarvalue colvarvalue::inverse() const -{ - switch (value_type) { - case colvarvalue::type_scalar: - return colvarvalue(1.0/real_value); - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return colvarvalue(cvm::rvector(1.0/rvector_value.x, - 1.0/rvector_value.y, - 1.0/rvector_value.z)); - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return colvarvalue(cvm::quaternion(1.0/quaternion_value.q0, - 1.0/quaternion_value.q1, - 1.0/quaternion_value.q2, - 1.0/quaternion_value.q3)); - break; - case colvarvalue::type_vector: - { - cvm::vector1d result(vector1d_value); - if (elem_types.size() > 0) { - // if we have information about non-scalar types, use it - for (size_t i = 0; i < elem_types.size(); i++) { - result.sliceassign(elem_indices[i], elem_indices[i]+elem_sizes[i], - cvm::vector1d((this->get_elem(i)).inverse())); - } - } else { - for (size_t i = 0; i < result.size(); i++) { - if (result[i] != 0.0) { - result = 1.0/result[i]; - } - } - } - return colvarvalue(result, type_vector); - } - break; - case colvarvalue::type_notset: - default: - undef_op(); - break; - } - return colvarvalue(); -} - - -inline colvarvalue & colvarvalue::operator = (colvarvalue const &x) -{ - check_types_assign(this->value_type, x.value_type); - value_type = x.value_type; - - switch (this->value_type) { - case colvarvalue::type_scalar: - this->real_value = x.real_value; - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - this->rvector_value = x.rvector_value; - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - this->quaternion_value = x.quaternion_value; - break; - case colvarvalue::type_vector: - vector1d_value = x.vector1d_value; - elem_types = x.elem_types; - elem_indices = x.elem_indices; - elem_sizes = x.elem_sizes; - break; - case colvarvalue::type_notset: - default: - undef_op(); - break; - } - return *this; -} - - -inline void colvarvalue::operator += (colvarvalue const &x) -{ - colvarvalue::check_types(*this, x); - - switch (this->value_type) { - case colvarvalue::type_scalar: - this->real_value += x.real_value; - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - this->rvector_value += x.rvector_value; - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - this->quaternion_value += x.quaternion_value; - break; - case colvarvalue::type_vector: - // TODO move this to vector1d - for (size_t i = 0; i < vector1d_value.size(); i++) { - vector1d_value[i] += x.vector1d_value[i]; - } - break; - case colvarvalue::type_notset: - default: - undef_op(); - } -} - -inline void colvarvalue::operator -= (colvarvalue const &x) -{ - colvarvalue::check_types(*this, x); - - switch (value_type) { - case colvarvalue::type_scalar: - real_value -= x.real_value; - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - rvector_value -= x.rvector_value; - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - quaternion_value -= x.quaternion_value; - break; - case colvarvalue::type_vector: - // TODO move this to vector1d - for (size_t i = 0; i < vector1d_value.size(); i++) { - vector1d_value[i] -= x.vector1d_value[i]; - } - break; - case colvarvalue::type_notset: - default: - undef_op(); - } -} - -inline void colvarvalue::operator *= (cvm::real const &a) -{ - switch (value_type) { - case colvarvalue::type_scalar: - real_value *= a; - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vectorderiv: - rvector_value *= a; - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - quaternion_value *= a; - break; - case colvarvalue::type_vector: - // TODO move this to vector1d - for (size_t i = 0; i < vector1d_value.size(); i++) { - vector1d_value[i] *= a; - } - break; - case colvarvalue::type_notset: - default: - undef_op(); - } -} - -inline void colvarvalue::operator /= (cvm::real const &a) -{ - switch (value_type) { - case colvarvalue::type_scalar: - real_value /= a; break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - rvector_value /= a; break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - quaternion_value /= a; break; - case colvarvalue::type_vector: - // TODO move this to vector1d - for (size_t i = 0; i < vector1d_value.size(); i++) { - vector1d_value[i] /= a; - } - break; - case colvarvalue::type_notset: - default: - undef_op(); - } -} - - -// binary operations between two colvarvalues - -inline colvarvalue operator + (colvarvalue const &x1, - colvarvalue const &x2) -{ - colvarvalue::check_types(x1, x2); - - switch (x1.value_type) { - case colvarvalue::type_scalar: - return colvarvalue(x1.real_value + x2.real_value); - case colvarvalue::type_3vector: - return colvarvalue(x1.rvector_value + x2.rvector_value); - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return colvarvalue(x1.rvector_value + x2.rvector_value, - colvarvalue::type_unit3vector); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return colvarvalue(x1.quaternion_value + x2.quaternion_value); - case colvarvalue::type_vector: - // TODO - break; - case colvarvalue::type_notset: - default: - x1.undef_op(); - return colvarvalue(colvarvalue::type_notset); - }; -} - -inline colvarvalue operator - (colvarvalue const &x1, - colvarvalue const &x2) -{ - colvarvalue::check_types(x1, x2); - - switch (x1.value_type) { - case colvarvalue::type_scalar: - return colvarvalue(x1.real_value - x2.real_value); - case colvarvalue::type_3vector: - return colvarvalue(x1.rvector_value - x2.rvector_value); - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return colvarvalue(x1.rvector_value - x2.rvector_value, - colvarvalue::type_unit3vector); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return colvarvalue(x1.quaternion_value - x2.quaternion_value); - case colvarvalue::type_vector: - // TODO - break; - case colvarvalue::type_notset: - default: - x1.undef_op(); - return colvarvalue(colvarvalue::type_notset); - }; -} - - -// binary operations with real numbers - -inline colvarvalue operator * (cvm::real const &a, - colvarvalue const &x) -{ - switch (x.value_type) { - case colvarvalue::type_scalar: - return colvarvalue(a * x.real_value); - case colvarvalue::type_3vector: - return colvarvalue(a * x.rvector_value); - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return colvarvalue(a * x.rvector_value, - colvarvalue::type_unit3vector); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return colvarvalue(a * x.quaternion_value); - case colvarvalue::type_vector: - // TODO - break; - case colvarvalue::type_notset: - default: - x.undef_op(); - return colvarvalue(colvarvalue::type_notset); - } -} - -inline colvarvalue operator * (colvarvalue const &x, - cvm::real const &a) -{ - switch (x.value_type) { - case colvarvalue::type_scalar: - return colvarvalue(x.real_value * a); - case colvarvalue::type_3vector: - return colvarvalue(x.rvector_value * a); - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return colvarvalue(x.rvector_value * a, - colvarvalue::type_unit3vector); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return colvarvalue(x.quaternion_value * a); - case colvarvalue::type_vector: - // TODO - break; - case colvarvalue::type_notset: - default: - x.undef_op(); - return colvarvalue(colvarvalue::type_notset); - } -} - -inline colvarvalue operator / (colvarvalue const &x, - cvm::real const &a) -{ - switch (x.value_type) { - case colvarvalue::type_scalar: - return colvarvalue(x.real_value / a); - case colvarvalue::type_3vector: - return colvarvalue(x.rvector_value / a); - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return colvarvalue(x.rvector_value / a, - colvarvalue::type_unit3vector); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return colvarvalue(x.quaternion_value / a); - case colvarvalue::type_vector: - // TODO - break; - case colvarvalue::type_notset: - default: - x.undef_op(); - return colvarvalue(colvarvalue::type_notset); - } -} - - -// inner product between two colvarvalues - -inline cvm::real operator * (colvarvalue const &x1, - colvarvalue const &x2) -{ - colvarvalue::check_types(x1, x2); - - switch (x1.value_type) { - case colvarvalue::type_scalar: - return (x1.real_value * x2.real_value); - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return (x1.rvector_value * x2.rvector_value); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - // the "*" product is the quaternion product, here the inner - // member function is used instead - return (x1.quaternion_value.inner(x2.quaternion_value)); - case colvarvalue::type_vector: - // TODO - break; - case colvarvalue::type_notset: - default: - x1.undef_op(); - return 0.0; - }; -} - - - inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const { colvarvalue::check_types(*this, x2); @@ -1064,7 +908,7 @@ inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const return this->quaternion_value.dist2(x2.quaternion_value); case colvarvalue::type_vector: // TODO - break; + return (this->vector1d_value - x2.vector1d_value).norm2(); case colvarvalue::type_notset: default: this->undef_op(); @@ -1073,420 +917,4 @@ inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const } -inline colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const -{ - colvarvalue::check_types(*this, x2); - - switch (this->value_type) { - case colvarvalue::type_scalar: - return 2.0 * (this->real_value - x2.real_value); - case colvarvalue::type_3vector: - return 2.0 * (this->rvector_value - x2.rvector_value); - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - { - cvm::rvector const &v1 = this->rvector_value; - cvm::rvector const &v2 = x2.rvector_value; - cvm::real const cos_t = v1 * v2; - cvm::real const sin_t = std::sqrt(1.0 - cos_t*cos_t); - return colvarvalue( 2.0 * sin_t * - cvm::rvector((-1.0) * sin_t * v2.x + - cos_t/sin_t * (v1.x - cos_t*v2.x), - (-1.0) * sin_t * v2.y + - cos_t/sin_t * (v1.y - cos_t*v2.y), - (-1.0) * sin_t * v2.z + - cos_t/sin_t * (v1.z - cos_t*v2.z) - ), - colvarvalue::type_unit3vector ); - } - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return this->quaternion_value.dist2_grad(x2.quaternion_value); - case colvarvalue::type_vector: - // TODO - break; - case colvarvalue::type_notset: - default: - this->undef_op(); - return colvarvalue(colvarvalue::type_notset); - }; -} - - - -inline cvm::vector1d const colvarvalue::as_vector() const -{ - switch (value_type) { - case colvarvalue::type_scalar: - return cvm::vector1d(1, real_value); - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return rvector_value.as_vector(); - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return quaternion_value.as_vector(); - break; - case colvarvalue::type_vector: - if (vector1d_value.size() == 0) { - cvm::error("Error: trying to use as vector a variable that is not yet initialized as such\".\n"); - } - return vector1d_value; - break; - case colvarvalue::type_notset: - default: - break; - } -} - - -inline void colvarvalue::check_types(colvarvalue const &x1, - colvarvalue const &x2) -{ - if (!colvarvalue::type_checking()) return; - - if (x1.value_type != x2.value_type) { - if (((x1.value_type == type_unit3vector) && - (x2.value_type == type_unit3vectorderiv)) || - ((x2.value_type == type_unit3vector) && - (x1.value_type == type_unit3vectorderiv)) || - ((x1.value_type == type_quaternion) && - (x2.value_type == type_quaternionderiv)) || - ((x2.value_type == type_quaternion) && - (x1.value_type == type_quaternionderiv))) { - return; - } - cvm::error("Performing an operation between two colvar " - "values with different types, \""+ - colvarvalue::type_desc(x1.value_type)+ - "\" and \""+ - colvarvalue::type_desc(x2.value_type)+ - "\".\n"); - } - - if (x1.value_type == type_vector) { - if (x1.vector1d_value.size() != x2.vector1d_value.size()) { - cvm::error("Performing an operation between two vector colvar " - "values with different sizes, "+ - cvm::to_str(x1.vector1d_value.size())+ - " and "+ - cvm::to_str(x2.vector1d_value.size())+ - ".\n"); - } - } -} - -inline void colvarvalue::check_types_assign(colvarvalue::Type const &vt1, - colvarvalue::Type const &vt2) -{ - if (vt1 != type_notset) { - if (vt1 != vt2) { - cvm::error("Trying to assign a colvar value with type \""+ - type_desc(vt2)+"\" to one with type \""+ - type_desc(vt1)+"\".\n"); - } - } -} - -inline void colvarvalue::undef_op() const -{ - cvm::error("Error: Undefined operation on a colvar of type \""+ - type_desc(this->value_type)+"\".\n"); -} - - - -inline void colvarvalue::inner_opt(colvarvalue const &x, - std::vector::iterator &xv, - std::vector::iterator const &xv_end, - std::vector::iterator &inner) -{ - // doing type check only once, here - colvarvalue::check_types(x, *xv); - - std::vector::iterator &xvi = xv; - std::vector::iterator &ii = inner; - - switch (x.value_type) { - case colvarvalue::type_scalar: - while (xvi != xv_end) { - *(ii++) += (xvi++)->real_value * x.real_value; - } - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - while (xvi != xv_end) { - *(ii++) += (xvi++)->rvector_value * x.rvector_value; - } - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - while (xvi != xv_end) { - *(ii++) += ((xvi++)->quaternion_value).cosine(x.quaternion_value); - } - break; - default: - x.undef_op(); - }; -} - -inline void colvarvalue::inner_opt(colvarvalue const &x, - std::list::iterator &xv, - std::list::iterator const &xv_end, - std::vector::iterator &inner) -{ - // doing type check only once, here - colvarvalue::check_types(x, *xv); - - std::list::iterator &xvi = xv; - std::vector::iterator &ii = inner; - - switch (x.value_type) { - case colvarvalue::type_scalar: - while (xvi != xv_end) { - *(ii++) += (xvi++)->real_value * x.real_value; - } - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - while (xvi != xv_end) { - *(ii++) += (xvi++)->rvector_value * x.rvector_value; - } - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - while (xvi != xv_end) { - *(ii++) += ((xvi++)->quaternion_value).cosine(x.quaternion_value); - } - break; - default: - x.undef_op(); - }; -} - - -inline void colvarvalue::p2leg_opt(colvarvalue const &x, - std::vector::iterator &xv, - std::vector::iterator const &xv_end, - std::vector::iterator &inner) -{ - // doing type check only once, here - colvarvalue::check_types(x, *xv); - - std::vector::iterator &xvi = xv; - std::vector::iterator &ii = inner; - - switch (x.value_type) { - case colvarvalue::type_scalar: - cvm::error("Error: cannot calculate Legendre polynomials " - "for scalar variables.\n"); - return; - break; - case colvarvalue::type_3vector: - while (xvi != xv_end) { - cvm::real const cosine = - ((xvi)->rvector_value * x.rvector_value) / - ((xvi)->rvector_value.norm() * x.rvector_value.norm()); - xvi++; - *(ii++) += 1.5*cosine*cosine - 0.5; - } - break; - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - while (xvi != xv_end) { - cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value; - *(ii++) += 1.5*cosine*cosine - 0.5; - } - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - while (xvi != xv_end) { - cvm::real const cosine = (xvi++)->quaternion_value.cosine(x.quaternion_value); - *(ii++) += 1.5*cosine*cosine - 0.5; - } - break; - default: - x.undef_op(); - }; -} - -inline void colvarvalue::p2leg_opt(colvarvalue const &x, - std::list::iterator &xv, - std::list::iterator const &xv_end, - std::vector::iterator &inner) -{ - // doing type check only once, here - colvarvalue::check_types(x, *xv); - - std::list::iterator &xvi = xv; - std::vector::iterator &ii = inner; - - switch (x.value_type) { - case colvarvalue::type_scalar: - cvm::error("Error: cannot calculate Legendre polynomials " - "for scalar variables.\n"); - break; - case colvarvalue::type_3vector: - while (xvi != xv_end) { - cvm::real const cosine = - ((xvi)->rvector_value * x.rvector_value) / - ((xvi)->rvector_value.norm() * x.rvector_value.norm()); - xvi++; - *(ii++) += 1.5*cosine*cosine - 0.5; - } - break; - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - while (xvi != xv_end) { - cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value; - *(ii++) += 1.5*cosine*cosine - 0.5; - } - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - while (xvi != xv_end) { - cvm::real const cosine = (xvi++)->quaternion_value.cosine(x.quaternion_value); - *(ii++) += 1.5*cosine*cosine - 0.5; - } - break; - default: - x.undef_op(); - }; -} - -inline std::string colvarvalue::to_simple_string() const -{ - switch (type()) { - case colvarvalue::type_scalar: - return cvm::to_str(real_value, 0, cvm::cv_prec); - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return rvector_value.to_simple_string(); - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return quaternion_value.to_simple_string(); - break; - case colvarvalue::type_notset: - undef_op(); - break; - } - return std::string(); -} - -inline int colvarvalue::from_simple_string(std::string const &s) -{ - switch (type()) { - case colvarvalue::type_scalar: - return ((std::istringstream(s) >> real_value) - ? COLVARS_OK : COLVARS_ERROR); - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return rvector_value.from_simple_string(s); - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return quaternion_value.from_simple_string(s); - break; - case colvarvalue::type_notset: - break; - default: - undef_op(); - } - return COLVARS_ERROR; -} - -inline std::ostream & operator << (std::ostream &os, colvarvalue const &x) -{ - switch (x.type()) { - case colvarvalue::type_scalar: - os << x.real_value; - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - os << x.rvector_value; - break; - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - os << x.quaternion_value; - break; - case colvarvalue::type_notset: - os << "not set"; break; - } - return os; -} - - -inline std::ostream & operator << (std::ostream &os, std::vector const &v) -{ - for (size_t i = 0; i < v.size(); i++) { - os << v[i]; - } - return os; -} - - -inline std::istream & operator >> (std::istream &is, colvarvalue &x) -{ - if (x.type() == colvarvalue::type_notset) { - cvm::error("Trying to read from a stream a colvarvalue, " - "which has not yet been assigned a data type.\n"); - return is; - } - - switch (x.type()) { - case colvarvalue::type_scalar: - is >> x.real_value; - break; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vectorderiv: - is >> x.rvector_value; - break; - case colvarvalue::type_unit3vector: - is >> x.rvector_value; - x.apply_constraints(); - break; - case colvarvalue::type_quaternion: - is >> x.quaternion_value; - x.apply_constraints(); - break; - case colvarvalue::type_quaternionderiv: - is >> x.quaternion_value; - break; - default: - x.undef_op(); - } - return is; -} - - -inline size_t colvarvalue::output_width(size_t const &real_width) const -{ - switch (this->value_type) { - case colvarvalue::type_scalar: - return real_width; - case colvarvalue::type_3vector: - case colvarvalue::type_unit3vector: - case colvarvalue::type_unit3vectorderiv: - return cvm::rvector::output_width(real_width); - case colvarvalue::type_quaternion: - case colvarvalue::type_quaternionderiv: - return cvm::quaternion::output_width(real_width); - case colvarvalue::type_notset: - default: - return 0; - } -} - - - #endif