colvarvalue.h

Go to the documentation of this file.

#ifndef COLVARVALUE_H
#define COLVARVALUE_H

#include "colvarmodule.h"



class colvarvalue {

public:

  enum Type {
    type_notset,
    type_scalar,
    type_vector,
    type_unitvector,
    type_quaternion
  };

  std::string static const type_desc[colvarvalue::type_quaternion+1];

  size_t static const      dof_num[  colvarvalue::type_quaternion+1];

  cvm::real       real_value;

  cvm::rvector    rvector_value;

  cvm::quaternion quaternion_value;

  Type  value_type; 

  static inline bool type_checking()
  {
    return true;
  }

  inline colvarvalue()
    : real_value (0.0), value_type (type_scalar)
  {}

  inline colvarvalue (Type const &vti)
    : value_type (vti)
  {
    reset();
  }

  inline colvarvalue (cvm::real const &x)
    : real_value (x), value_type (type_scalar)
  {}

  inline colvarvalue (cvm::rvector const &v)
    : rvector_value (v), value_type (type_vector)
  {}

  inline colvarvalue (cvm::rvector const &v, Type const &vti)
    : rvector_value (v), value_type (vti)
  {}

  inline colvarvalue (cvm::quaternion const &q)
    : quaternion_value (q), value_type (type_quaternion)
  {}

  inline colvarvalue (colvarvalue const &x)
    : value_type (x.value_type)
  {
    reset();

    switch (x.value_type) {
    case type_scalar:
      real_value = x.real_value;
      break;
    case type_vector: 
    case type_unitvector:
      rvector_value = x.rvector_value;
      break;
    case type_quaternion:
      quaternion_value = x.quaternion_value;
      break;
    case type_notset:
    default:
      break;
    }
  }



  void reset();


  void apply_constraints();


  inline Type type() const
  {
    return value_type;
  }

  inline void type (Type const &vti)
  {
    reset();
    value_type = vti;
    reset();
  }

  inline void type (colvarvalue const &x)
  {
    reset();
    value_type = x.value_type;
    reset();
  }

  cvm::real norm2() const;

  inline cvm::real norm() const
  {
    return std::sqrt (this->norm2());
  }

  inline colvarvalue inverse() const;

  cvm::real dist2 (colvarvalue const &x2) const;

  colvarvalue dist2_grad (colvarvalue const &x2) const;

  colvarvalue & operator = (colvarvalue const &x);

  void operator += (colvarvalue const &x);
  void operator -= (colvarvalue const &x);
  void operator *= (cvm::real const &a);
  void operator /= (cvm::real const &a);


  // Casting operator
  inline operator cvm::real() const
  {
    if (value_type != type_scalar) error_rside (type_scalar);
    return real_value;
  }

  // Casting operator
  inline operator cvm::rvector() const
  {
    if ( (value_type != type_vector) && (value_type != type_unitvector))
      error_rside (type_vector);
    return rvector_value;
  }

  // Casting operator
  inline operator cvm::quaternion() const
  {
    if (value_type != type_quaternion) error_rside (type_quaternion);
    return quaternion_value;
  }

  inline bool is_scalar()
  {
    return (value_type == type_scalar);
  }


  void static check_types (colvarvalue const &x1, colvarvalue const &x2);

  void undef_op() const;

  void error_lside (Type const &vt) const;

  void error_rside (Type const &vt) const;

  size_t output_width (size_t const &real_width);


  // 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)

  void static inner_opt (colvarvalue const                        &x,
                         std::vector<colvarvalue>::iterator       &xv,
                         std::vector<colvarvalue>::iterator const &xv_end,
                         std::vector<cvm::real>::iterator         &inner);

  void static inner_opt (colvarvalue const                        &x,
                         std::list<colvarvalue>::iterator         &xv,
                         std::list<colvarvalue>::iterator const   &xv_end,
                         std::vector<cvm::real>::iterator         &inner);

  void static p2leg_opt (colvarvalue const                        &x,
                         std::vector<colvarvalue>::iterator       &xv,
                         std::vector<colvarvalue>::iterator const &xv_end,
                         std::vector<cvm::real>::iterator         &inner);

  void static p2leg_opt (colvarvalue const                        &x,
                         std::list<colvarvalue>::iterator         &xv,
                         std::list<colvarvalue>::iterator const   &xv_end,
                         std::vector<cvm::real>::iterator         &inner);

};



inline void colvarvalue::reset()
{
  switch (value_type) {
  case colvarvalue::type_scalar:
    real_value = cvm::real (0.0);
    break;
  case colvarvalue::type_vector:
  case colvarvalue::type_unitvector:
    rvector_value = cvm::rvector (0.0, 0.0, 0.0);
    break;
  case colvarvalue::type_quaternion:
    quaternion_value = cvm::quaternion (0.0, 0.0, 0.0, 0.0);
    break;
  case colvarvalue::type_notset:
  default:
    break;
  }  
}


inline void colvarvalue::apply_constraints()
{
  switch (value_type) {
  case colvarvalue::type_scalar:
    break;
  case colvarvalue::type_vector:
    break;
  case colvarvalue::type_unitvector:
    rvector_value /= std::sqrt (rvector_value.norm2());
    break;
  case colvarvalue::type_quaternion:
    quaternion_value /= std::sqrt (quaternion_value.norm2());
    break;
  case colvarvalue::type_notset:
  default:
    break;
  }  
}



inline cvm::real colvarvalue::norm2() const
{
  switch (value_type) {
  case colvarvalue::type_scalar:
    return (this->real_value)*(this->real_value);
  case colvarvalue::type_vector:
  case colvarvalue::type_unitvector:
    return (this->rvector_value).norm2();
  case colvarvalue::type_quaternion:
    return (this->quaternion_value).norm2();
  case colvarvalue::type_notset:
  default:
    return 0.0;
  }  
}


inline colvarvalue colvarvalue::inverse() const
{
  switch (value_type) {
  case colvarvalue::type_scalar:
    return colvarvalue (1.0/real_value);
  case colvarvalue::type_vector:
    return colvarvalue (cvm::rvector (1.0/rvector_value.x,
                                      1.0/rvector_value.y,
                                      1.0/rvector_value.z));
  case colvarvalue::type_quaternion:
    return colvarvalue (cvm::quaternion (1.0/quaternion_value.q0,
                                         1.0/quaternion_value.q1,
                                         1.0/quaternion_value.q2,
                                         1.0/quaternion_value.q3));
  case colvarvalue::type_notset:
  default:
    undef_op();
  }  
  return colvarvalue();
}


inline colvarvalue & colvarvalue::operator = (colvarvalue const &x)
{
  if (this->value_type != type_notset)
    if (this->value_type != x.value_type) 
      error_lside (x.value_type);

  this->value_type = x.value_type;
  
  switch (this->value_type) {
  case colvarvalue::type_scalar:
    this->real_value = x.real_value;
    break;
  case colvarvalue::type_vector:
  case colvarvalue::type_unitvector:
    this->rvector_value = x.rvector_value;
    break;
  case colvarvalue::type_quaternion:
    this->quaternion_value = x.quaternion_value;
    break;
  case colvarvalue::type_notset:
  default:
    undef_op();
    break;
  }
  return *this;
}

inline void colvarvalue::operator += (colvarvalue const &x)
{
  if (colvarvalue::type_checking())
    colvarvalue::check_types (*this, x);

  switch (this->value_type) {
  case colvarvalue::type_scalar:
    this->real_value += x.real_value;
    break;
  case colvarvalue::type_vector:
  case colvarvalue::type_unitvector:
    this->rvector_value += x.rvector_value;
    break;
  case colvarvalue::type_quaternion:
    this->quaternion_value += x.quaternion_value;
    break;
  case colvarvalue::type_notset:
  default:
    undef_op();
  }
}

inline void colvarvalue::operator -= (colvarvalue const &x)
{
  if (colvarvalue::type_checking())
    colvarvalue::check_types (*this, x);

  switch (value_type) {
  case colvarvalue::type_scalar:
    real_value -= x.real_value; break;
  case colvarvalue::type_vector:
  case colvarvalue::type_unitvector:
    rvector_value -= x.rvector_value; break;
  case colvarvalue::type_quaternion:
    quaternion_value -= x.quaternion_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_vector:
  case colvarvalue::type_unitvector:
    rvector_value *= a;
    break;
  case colvarvalue::type_quaternion:
    quaternion_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_vector:
  case colvarvalue::type_unitvector:
    rvector_value /= a; break;
  case colvarvalue::type_quaternion:
    quaternion_value /= a; break;
  case colvarvalue::type_notset:
  default:
    undef_op();
  }
}


// binary operations between two colvarvalues

inline colvarvalue operator + (colvarvalue const &x1,
                               colvarvalue const &x2)
{
  if (colvarvalue::type_checking())
    colvarvalue::check_types (x1, x2);

  switch (x1.value_type) {
  case colvarvalue::type_scalar:
    return colvarvalue (x1.real_value + x2.real_value);
  case colvarvalue::type_vector:
    return colvarvalue (x1.rvector_value + x2.rvector_value);
  case colvarvalue::type_unitvector:
    return colvarvalue (x1.rvector_value + x2.rvector_value,
                        colvarvalue::type_unitvector);
  case colvarvalue::type_quaternion:
    return colvarvalue (x1.quaternion_value + x2.quaternion_value);
  case colvarvalue::type_notset:
  default:
    x1.undef_op();
    return colvarvalue (colvarvalue::type_notset);
  };
}

inline colvarvalue operator - (colvarvalue const &x1,
                               colvarvalue const &x2)
{
  if (colvarvalue::type_checking())
    colvarvalue::check_types (x1, x2);

  switch (x1.value_type) {
  case colvarvalue::type_scalar:
    return colvarvalue (x1.real_value - x2.real_value);
  case colvarvalue::type_vector:
    return colvarvalue (x1.rvector_value - x2.rvector_value);
  case colvarvalue::type_unitvector:
    return colvarvalue (x1.rvector_value - x2.rvector_value,
                        colvarvalue::type_unitvector);
  case colvarvalue::type_quaternion:
    return colvarvalue (x1.quaternion_value - x2.quaternion_value);
  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_vector:
    return colvarvalue (a * x.rvector_value);
  case colvarvalue::type_unitvector:
    return colvarvalue (a * x.rvector_value,
                        colvarvalue::type_unitvector);
  case colvarvalue::type_quaternion:
    return colvarvalue (a * x.quaternion_value);
  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_vector:
    return colvarvalue (x.rvector_value * a);
  case colvarvalue::type_unitvector:
    return colvarvalue (x.rvector_value * a,
                        colvarvalue::type_unitvector);
  case colvarvalue::type_quaternion:
    return colvarvalue (x.quaternion_value * a);
  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_vector:
    return colvarvalue (x.rvector_value / a);
  case colvarvalue::type_unitvector:
    return colvarvalue (x.rvector_value / a,
                        colvarvalue::type_unitvector);
  case colvarvalue::type_quaternion:
    return colvarvalue (x.quaternion_value / a);
  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)
{
  if (colvarvalue::type_checking())
    colvarvalue::check_types (x1, x2);

  switch (x1.value_type) {
  case colvarvalue::type_scalar:
    return (x1.real_value * x2.real_value);
  case colvarvalue::type_vector:
  case colvarvalue::type_unitvector:
    return (x1.rvector_value * x2.rvector_value);
  case colvarvalue::type_quaternion:
    // 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_notset:
  default:
    x1.undef_op();
    return 0.0;
  };
}



inline cvm::real colvarvalue::dist2 (colvarvalue const &x2) const
{
  if (colvarvalue::type_checking())
    colvarvalue::check_types (*this, x2);

  switch (this->value_type) {
  case colvarvalue::type_scalar:
    return std::pow (this->real_value - x2.real_value, int (2));
  case colvarvalue::type_vector:
    return (this->rvector_value - x2.rvector_value).norm2();
  case colvarvalue::type_unitvector:
    // angle between (*this) and x2 is the distance
    return std::pow (std::acos (this->rvector_value * x2.rvector_value), int (2));
  case colvarvalue::type_quaternion:
    // angle between (*this) and x2 is the distance, the quaternion
    // object has it implemented internally
    return this->quaternion_value.dist2 (x2.quaternion_value);
  case colvarvalue::type_notset:
  default:
    this->undef_op();
    return 0.0;
  };
}


inline colvarvalue colvarvalue::dist2_grad (colvarvalue const &x2) const
{
  if (colvarvalue::type_checking())
    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_vector:
    return 2.0 * (this->rvector_value - x2.rvector_value);
  case colvarvalue::type_unitvector:
    {
      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_unitvector );
    }
  case colvarvalue::type_quaternion:
    return this->quaternion_value.dist2_grad (x2.quaternion_value);
  case colvarvalue::type_notset:
  default:
    this->undef_op();
    return colvarvalue (colvarvalue::type_notset);
  };
}


inline void colvarvalue::check_types (colvarvalue const &x1,
                                      colvarvalue const &x2)
{
  if (x1.value_type != x2.value_type) {
    cvm::log ("x1 type = "+cvm::to_str (x1.value_type)+
              ", x2 type = "+cvm::to_str (x2.value_type)+"\n");
    cvm::fatal_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");
  }
}




std::ostream & operator << (std::ostream &os, colvarvalue const &x);
std::ostream & operator << (std::ostream &os, std::vector<colvarvalue> const &v);

std::istream & operator >> (std::istream &is, colvarvalue &x);





#endif


// Emacs
// Local Variables:
// mode: C++
// End:

---