Collective Variables Module - Developer Documentation
colvartypes.h
1 // -*- c++ -*-
2 
3 #ifndef COLVARTYPES_H
4 #define COLVARTYPES_H
5 
6 #include <cmath>
7 #include <vector>
8 
9 #include "colvarmodule.h"
10 
11 #ifndef PI
12 #define PI 3.14159265358979323846
13 #endif
14 
15 // ----------------------------------------------------------------------
18 // ----------------------------------------------------------------------
19 
20 
24 template <class T> class colvarmodule::vector1d
25 {
26 protected:
27 
28  std::vector<T> data;
29 
30 public:
31 
33  inline vector1d(size_t const n = 0)
34  {
35  data.resize(n);
36  reset();
37  }
38 
40  inline vector1d(size_t const n, T const *t)
41  {
42  data.resize(n);
43  reset();
44  size_t i;
45  for (i = 0; i < size(); i++) {
46  data[i] = t[i];
47  }
48  }
49 
51  inline T * c_array()
52  {
53  if (data.size() > 0) {
54  return &(data[0]);
55  } else {
56  return NULL;
57  }
58  }
59 
61  inline std::vector<T> &data_array()
62  {
63  return data;
64  }
65 
66  inline ~vector1d()
67  {
68  data.clear();
69  }
70 
72  inline void reset()
73  {
74  data.assign(data.size(), T(0.0));
75  }
76 
77  inline size_t size() const
78  {
79  return data.size();
80  }
81 
82  inline void resize(size_t const n)
83  {
84  data.resize(n);
85  }
86 
87  inline void clear()
88  {
89  data.clear();
90  }
91 
92  inline T & operator [] (size_t const i) {
93  return data[i];
94  }
95 
96  inline T const & operator [] (size_t const i) const {
97  return data[i];
98  }
99 
100  inline static void check_sizes(vector1d<T> const &v1, vector1d<T> const &v2)
101  {
102  if (v1.size() != v2.size()) {
103  cvm::error("Error: trying to perform an operation between vectors of different sizes, "+
104  cvm::to_str(v1.size())+" and "+cvm::to_str(v2.size())+".\n");
105  }
106  }
107 
108  inline void operator += (vector1d<T> const &v)
109  {
110  check_sizes(*this, v);
111  size_t i;
112  for (i = 0; i < this->size(); i++) {
113  (*this)[i] += v[i];
114  }
115  }
116 
117  inline void operator -= (vector1d<T> const &v)
118  {
119  check_sizes(*this, v);
120  size_t i;
121  for (i = 0; i < this->size(); i++) {
122  (*this)[i] -= v[i];
123  }
124  }
125 
126  inline void operator *= (cvm::real const &a)
127  {
128  size_t i;
129  for (i = 0; i < this->size(); i++) {
130  (*this)[i] *= a;
131  }
132  }
133 
134  inline void operator /= (cvm::real const &a)
135  {
136  size_t i;
137  for (i = 0; i < this->size(); i++) {
138  (*this)[i] /= a;
139  }
140  }
141 
142  inline friend vector1d<T> operator + (vector1d<T> const &v1, vector1d<T> const &v2)
143  {
144  check_sizes(v1.size(), v2.size());
145  vector1d<T> result(v1.size());
146  size_t i;
147  for (i = 0; i < v1.size(); i++) {
148  result[i] = v1[i] + v2[i];
149  }
150  return result;
151  }
152 
153  inline friend vector1d<T> operator - (vector1d<T> const &v1, vector1d<T> const &v2)
154  {
155  check_sizes(v1.size(), v2.size());
156  vector1d<T> result(v1.size());
157  size_t i;
158  for (i = 0; i < v1.size(); i++) {
159  result[i] = v1[i] - v2[i];
160  }
161  return result;
162  }
163 
164  inline friend vector1d<T> operator * (vector1d<T> const &v, cvm::real const &a)
165  {
166  vector1d<T> result(v.size());
167  size_t i;
168  for (i = 0; i < v.size(); i++) {
169  result[i] = v[i] * a;
170  }
171  return result;
172  }
173 
174  inline friend vector1d<T> operator * (cvm::real const &a, vector1d<T> const &v)
175  {
176  return v * a;
177  }
178 
179  inline friend vector1d<T> operator / (vector1d<T> const &v, cvm::real const &a)
180  {
181  vector1d<T> result(v.size());
182  size_t i;
183  for (i = 0; i < v.size(); i++) {
184  result[i] = v[i] / a;
185  }
186  return result;
187  }
188 
190  inline friend T operator * (vector1d<T> const &v1, vector1d<T> const &v2)
191  {
192  check_sizes(v1.size(), v2.size());
193  T prod(0.0);
194  size_t i;
195  for (i = 0; i < v1.size(); i++) {
196  prod += v1[i] * v2[i];
197  }
198  return prod;
199  }
200 
202  inline cvm::real norm2() const
203  {
204  cvm::real result = 0.0;
205  size_t i;
206  for (i = 0; i < this->size(); i++) {
207  result += (*this)[i] * (*this)[i];
208  }
209  return result;
210  }
211 
212  inline cvm::real norm() const
213  {
214  return std::sqrt(this->norm2());
215  }
216 
217  inline cvm::real sum() const
218  {
219  cvm::real result = 0.0;
220  size_t i;
221  for (i = 0; i < this->size(); i++) {
222  result += (*this)[i];
223  }
224  return result;
225  }
226 
228  inline vector1d<T> const slice(size_t const i1, size_t const i2) const
229  {
230  if ((i2 < i1) || (i2 >= this->size())) {
231  cvm::error("Error: trying to slice a vector using incorrect boundaries.\n");
232  }
233  vector1d<T> result(i2 - i1);
234  size_t i;
235  for (i = 0; i < (i2 - i1); i++) {
236  result[i] = (*this)[i1+i];
237  }
238  return result;
239  }
240 
242  inline void sliceassign(size_t const i1, size_t const i2, vector1d<T> const &v)
243  {
244  if ((i2 < i1) || (i2 >= this->size())) {
245  cvm::error("Error: trying to slice a vector using incorrect boundaries.\n");
246  }
247  size_t i;
248  for (i = 0; i < (i2 - i1); i++) {
249  (*this)[i1+i] = v[i];
250  }
251  }
252 
254 
255  inline size_t output_width(size_t const &real_width) const
256  {
257  return real_width*(this->size()) + 3*(this->size()-1) + 4;
258  }
259 
260  inline friend std::istream & operator >> (std::istream &is, cvm::vector1d<T> &v)
261  {
262  if (v.size() == 0) return is;
263  size_t const start_pos = is.tellg();
264  char sep;
265  if ( !(is >> sep) || !(sep == '(') ) {
266  is.clear();
267  is.seekg(start_pos, std::ios::beg);
268  is.setstate(std::ios::failbit);
269  return is;
270  }
271  size_t count = 0;
272  while ( (is >> v[count]) &&
273  (count < (v.size()-1) ? ((is >> sep) && (sep == ',')) : true) ) {
274  if (++count == v.size()) break;
275  }
276  if (count < v.size()) {
277  is.clear();
278  is.seekg(start_pos, std::ios::beg);
279  is.setstate(std::ios::failbit);
280  }
281  return is;
282  }
283 
284  inline friend std::ostream & operator << (std::ostream &os, cvm::vector1d<T> const &v)
285  {
286  std::streamsize const w = os.width();
287  std::streamsize const p = os.precision();
288 
289  os.width(2);
290  os << "( ";
291  size_t i;
292  for (i = 0; i < v.size()-1; i++) {
293  os.width(w); os.precision(p);
294  os << v[i] << " , ";
295  }
296  os.width(w); os.precision(p);
297  os << v[v.size()-1] << " )";
298  return os;
299  }
300 
301  inline std::string to_simple_string() const
302  {
303  if (this->size() == 0) return std::string("");
304  std::ostringstream os;
305  os.setf(std::ios::scientific, std::ios::floatfield);
306  os.precision(cvm::cv_prec);
307  os << (*this)[0];
308  size_t i;
309  for (i = 1; i < this->size(); i++) {
310  os << " " << (*this)[i];
311  }
312  return os.str();
313  }
314 
315  inline int from_simple_string(std::string const &s)
316  {
317  std::stringstream stream(s);
318  size_t i = 0;
319  if (this->size()) {
320  while ((stream >> (*this)[i]) && (i < this->size())) {
321  i++;
322  }
323  if (i < this->size()) {
324  return COLVARS_ERROR;
325  }
326  } else {
327  T input;
328  while (stream >> input) {
329  if ((i % 100) == 0) {
330  data.reserve(data.size()+100);
331  }
332  data.resize(data.size()+1);
333  data[i] = input;
334  i++;
335  }
336  }
337  return COLVARS_OK;
338  }
339 
340 };
341 
342 
346 template <class T> class colvarmodule::matrix2d
347 {
348 public:
349 
350  friend class row;
351  size_t outer_length;
352  size_t inner_length;
353 
354 protected:
355 
356  class row {
357  public:
358  T * data;
359  size_t length;
360  inline row(T * const row_data, size_t const inner_length)
361  : data(row_data), length(inner_length)
362  {}
363  inline T & operator [] (size_t const j) {
364  return *(data+j);
365  }
366  inline T const & operator [] (size_t const j) const {
367  return *(data+j);
368  }
369  inline operator vector1d<T>() const
370  {
371  return vector1d<T>(length, data);
372  }
373  };
374 
375  std::vector<T> data;
376  std::vector<row> rows;
377  std::vector<T *> pointers;
378 
379 public:
380 
382  inline void resize(size_t const ol, size_t const il)
383  {
384  if ((ol > 0) && (il > 0)) {
385 
386  if (data.size() > 0) {
387  // copy previous data
388  size_t i, j;
389  std::vector<T> new_data(ol * il);
390  for (i = 0; i < outer_length; i++) {
391  for (j = 0; j < inner_length; j++) {
392  new_data[il*i+j] = data[inner_length*i+j];
393  }
394  }
395  data.resize(ol * il);
396  // copy them back
397  data = new_data;
398  } else {
399  data.resize(ol * il);
400  }
401 
402  outer_length = ol;
403  inner_length = il;
404 
405  if (data.size() > 0) {
406  // rebuild rows
407  size_t i;
408  rows.clear();
409  rows.reserve(outer_length);
410  pointers.clear();
411  pointers.reserve(outer_length);
412  for (i = 0; i < outer_length; i++) {
413  rows.push_back(row(&(data[0])+inner_length*i, inner_length));
414  pointers.push_back(&(data[0])+inner_length*i);
415  }
416  }
417  } else {
418  // zero size
419  data.clear();
420  rows.clear();
421  }
422  }
423 
425  inline void clear() {
426  rows.clear();
427  data.clear();
428  }
429 
431  inline void reset()
432  {
433  data.assign(data.size(), T(0.0));
434  }
435 
436  inline size_t size() const
437  {
438  return data.size();
439  }
440 
442  inline matrix2d()
443  : outer_length(0), inner_length(0)
444  {
445  this->resize(0, 0);
446  }
447 
448  inline matrix2d(size_t const ol, size_t const il)
449  : outer_length(ol), inner_length(il)
450  {
451  this->resize(outer_length, inner_length);
452  reset();
453  }
454 
456  inline matrix2d(matrix2d<T> const &m)
457  : outer_length(m.outer_length), inner_length(m.inner_length)
458  {
459  // reinitialize data and rows arrays
460  this->resize(outer_length, inner_length);
461  // copy data
462  data = m.data;
463  }
464 
466  inline ~matrix2d() {
467  this->clear();
468  }
469 
471  inline std::vector<T> &data_array()
472  {
473  return data;
474  }
475 
476  inline row & operator [] (size_t const i)
477  {
478  return rows[i];
479  }
480  inline row const & operator [] (size_t const i) const
481  {
482  return rows[i];
483  }
484 
486  inline matrix2d<T> & operator = (matrix2d<T> const &m)
487  {
488  if ((outer_length != m.outer_length) || (inner_length != m.inner_length)){
489  this->clear();
490  outer_length = m.outer_length;
491  inner_length = m.inner_length;
492  this->resize(outer_length, inner_length);
493  }
494  data = m.data;
495  return *this;
496  }
497 
499  inline T ** c_array() {
500  if (rows.size() > 0) {
501  return &(pointers[0]);
502  } else {
503  return NULL;
504  }
505  }
506 
507  inline static void check_sizes(matrix2d<T> const &m1, matrix2d<T> const &m2)
508  {
509  if ((m1.outer_length != m2.outer_length) ||
510  (m1.inner_length != m2.inner_length)) {
511  cvm::error("Error: trying to perform an operation between matrices of different sizes, "+
512  cvm::to_str(m1.outer_length)+"x"+cvm::to_str(m1.inner_length)+" and "+
513  cvm::to_str(m2.outer_length)+"x"+cvm::to_str(m2.inner_length)+".\n");
514  }
515  }
516 
517  inline void operator += (matrix2d<T> const &m)
518  {
519  check_sizes(*this, m);
520  size_t i;
521  for (i = 0; i < data.size(); i++) {
522  data[i] += m.data[i];
523  }
524  }
525 
526  inline void operator -= (matrix2d<T> const &m)
527  {
528  check_sizes(*this, m);
529  size_t i;
530  for (i = 0; i < data.size(); i++) {
531  data[i] -= m.data[i];
532  }
533  }
534 
535  inline void operator *= (cvm::real const &a)
536  {
537  size_t i;
538  for (i = 0; i < data.size(); i++) {
539  data[i] *= a;
540  }
541  }
542 
543  inline void operator /= (cvm::real const &a)
544  {
545  size_t i;
546  for (i = 0; i < data.size(); i++) {
547  data[i] /= a;
548  }
549  }
550 
551  inline friend matrix2d<T> operator + (matrix2d<T> const &m1, matrix2d<T> const &m2)
552  {
553  check_sizes(m1, m2);
554  matrix2d<T> result(m1.outer_length, m1.inner_length);
555  size_t i;
556  for (i = 0; i < m1.data.size(); i++) {
557  result.data[i] = m1.data[i] + m2.data[i];
558  }
559  return result;
560  }
561 
562  inline friend matrix2d<T> operator - (matrix2d<T> const &m1, matrix2d<T> const &m2)
563  {
564  check_sizes(m1, m2);
565  matrix2d<T> result(m1.outer_length, m1.inner_length);
566  size_t i;
567  for (i = 0; i < m1.data.size(); i++) {
568  result.data[i] = m1.data[i] - m1.data[i];
569  }
570  return result;
571  }
572 
573  inline friend matrix2d<T> operator * (matrix2d<T> const &m, cvm::real const &a)
574  {
575  matrix2d<T> result(m.outer_length, m.inner_length);
576  size_t i;
577  for (i = 0; i < m.data.size(); i++) {
578  result.data[i] = m.data[i] * a;
579  }
580  return result;
581  }
582 
583  inline friend matrix2d<T> operator * (cvm::real const &a, matrix2d<T> const &m)
584  {
585  return m * a;
586  }
587 
588  inline friend matrix2d<T> operator / (matrix2d<T> const &m, cvm::real const &a)
589  {
590  matrix2d<T> result(m.outer_length, m.inner_length);
591  size_t i;
592  for (i = 0; i < m.data.size(); i++) {
593  result.data[i] = m.data[i] * a;
594  }
595  return result;
596  }
597 
599 // inline friend matrix2d<T> const & operator * (matrix2d<T> const &m1, matrix2d<T> const &m2)
600 // {
601 // matrix2d<T> result(m1.outer_length, m2.inner_length);
602 // if (m1.inner_length != m2.outer_length) {
603 // cvm::error("Error: trying to multiply two matrices of incompatible sizes, "+
604 // cvm::to_str(m1.outer_length)+"x"+cvm::to_str(m1.inner_length)+" and "+
605 // cvm::to_str(m2.outer_length)+"x"+cvm::to_str(m2.inner_length)+".\n");
606 // } else {
607 // size_t i, j, k;
608 // for (i = 0; i < m1.outer_length; i++) {
609 // for (j = 0; j < m2.inner_length; j++) {
610 // for (k = 0; k < m1.inner_length; k++) {
611 // result[i][j] += m1[i][k] * m2[k][j];
612 // }
613 // }
614 // }
615 // }
616 // return result;
617 // }
618 
620  inline friend vector1d<T> operator * (vector1d<T> const &v, matrix2d<T> const &m)
621  {
622  vector1d<T> result(m.inner_length);
623  if (m.outer_length != v.size()) {
624  cvm::error("Error: trying to multiply a vector and a matrix of incompatible sizes, "+
625  cvm::to_str(v.size()) + " and " + cvm::to_str(m.outer_length)+"x"+cvm::to_str(m.inner_length) + ".\n");
626  } else {
627  size_t i, k;
628  for (i = 0; i < m.inner_length; i++) {
629  for (k = 0; k < m.outer_length; k++) {
630  result[i] += m[k][i] * v[k];
631  }
632  }
633  }
634  return result;
635  }
636 
637 // /// matrix-vector multiplication (unused for now)
638 // inline friend vector1d<T> const & operator * (matrix2d<T> const &m, vector1d<T> const &v)
639 // {
640 // vector1d<T> result(m.outer_length);
641 // if (m.inner_length != v.size()) {
642 // cvm::error("Error: trying to multiply a matrix and a vector of incompatible sizes, "+
643 // cvm::to_str(m.outer_length)+"x"+cvm::to_str(m.inner_length)
644 // + " and " + cvm::to_str(v.length) + ".\n");
645 // } else {
646 // size_t i, k;
647 // for (i = 0; i < m.outer_length; i++) {
648 // for (k = 0; k < m.inner_length; k++) {
649 // result[i] += m[i][k] * v[k];
650 // }
651 // }
652 // }
653 // return result;
654 // }
655 
657  friend std::ostream & operator << (std::ostream &os,
658  matrix2d<T> const &m)
659  {
660  std::streamsize const w = os.width();
661  std::streamsize const p = os.precision();
662 
663  os.width(2);
664  os << "( ";
665  size_t i;
666  for (i = 0; i < m.outer_length; i++) {
667  os << " ( ";
668  size_t j;
669  for (j = 0; j < m.inner_length-1; j++) {
670  os.width(w);
671  os.precision(p);
672  os << m[i][j] << " , ";
673  }
674  os.width(w);
675  os.precision(p);
676  os << m[i][m.inner_length-1] << " )";
677  }
678 
679  os << " )";
680  return os;
681  }
682 
683  inline std::string to_simple_string() const
684  {
685  if (this->size() == 0) return std::string("");
686  std::ostringstream os;
687  os.setf(std::ios::scientific, std::ios::floatfield);
688  os.precision(cvm::cv_prec);
689  os << (*this)[0];
690  size_t i;
691  for (i = 1; i < data.size(); i++) {
692  os << " " << data[i];
693  }
694  return os.str();
695  }
696 
697  inline int from_simple_string(std::string const &s)
698  {
699  std::stringstream stream(s);
700  size_t i = 0;
701  while ((i < data.size()) && (stream >> data[i])) {
702  i++;
703  }
704  if (i < data.size()) {
705  return COLVARS_ERROR;
706  }
707  return COLVARS_OK;
708  }
709 
710 };
711 
712 
715 
716 public:
717 
718  cvm::real x, y, z;
719 
720  inline rvector()
721  : x(0.0), y(0.0), z(0.0)
722  {}
723 
724  inline rvector(cvm::real const &x_i,
725  cvm::real const &y_i,
726  cvm::real const &z_i)
727  : x(x_i), y(y_i), z(z_i)
728  {}
729 
730  inline rvector(cvm::vector1d<cvm::real> const &v)
731  : x(v[0]), y(v[1]), z(v[2])
732  {}
733 
734  inline rvector(cvm::real t)
735  : x(t), y(t), z(t)
736  {}
737 
739  inline void set(cvm::real const &value) {
740  x = y = z = value;
741  }
742 
744  inline void set(cvm::real const &x_i,
745  cvm::real const &y_i,
746  cvm::real const &z_i) {
747  x = x_i;
748  y = y_i;
749  z = z_i;
750  }
751 
753  inline void reset() {
754  x = y = z = 0.0;
755  }
756 
758  inline cvm::real & operator [] (int const &i) {
759  return (i == 0) ? x : (i == 1) ? y : (i == 2) ? z : x;
760  }
761 
763  inline cvm::real const & operator [] (int const &i) const {
764  return (i == 0) ? x : (i == 1) ? y : (i == 2) ? z : x;
765  }
766 
767  inline cvm::vector1d<cvm::real> const as_vector() const
768  {
769  cvm::vector1d<cvm::real> result(3);
770  result[0] = this->x;
771  result[1] = this->y;
772  result[2] = this->z;
773  return result;
774  }
775 
776  inline cvm::rvector & operator = (cvm::real const &v)
777  {
778  x = v;
779  y = v;
780  z = v;
781  return *this;
782  }
783 
784  inline void operator += (cvm::rvector const &v)
785  {
786  x += v.x;
787  y += v.y;
788  z += v.z;
789  }
790 
791  inline void operator -= (cvm::rvector const &v)
792  {
793  x -= v.x;
794  y -= v.y;
795  z -= v.z;
796  }
797 
798  inline void operator *= (cvm::real const &v)
799  {
800  x *= v;
801  y *= v;
802  z *= v;
803  }
804 
805  inline void operator /= (cvm::real const& v)
806  {
807  x /= v;
808  y /= v;
809  z /= v;
810  }
811 
812  inline cvm::real norm2() const
813  {
814  return (x*x + y*y + z*z);
815  }
816 
817  inline cvm::real norm() const
818  {
819  return std::sqrt(this->norm2());
820  }
821 
822  inline cvm::rvector unit() const
823  {
824  const cvm::real n = this->norm();
825  return (n > 0. ? cvm::rvector(x, y, z)/n : cvm::rvector(1., 0., 0.));
826  }
827 
828  static inline size_t output_width(size_t const &real_width)
829  {
830  return 3*real_width + 10;
831  }
832 
833 
834  static inline cvm::rvector outer(cvm::rvector const &v1, cvm::rvector const &v2)
835  {
836  return cvm::rvector( v1.y*v2.z - v2.y*v1.z,
837  -v1.x*v2.z + v2.x*v1.z,
838  v1.x*v2.y - v2.x*v1.y);
839  }
840 
841  friend inline cvm::rvector operator - (cvm::rvector const &v)
842  {
843  return cvm::rvector(-v.x, -v.y, -v.z);
844  }
845 
846  friend inline int operator == (cvm::rvector const &v1, cvm::rvector const &v2)
847  {
848  return (v1.x == v2.x) && (v1.y == v2.y) && (v1.z == v2.z);
849  }
850 
851  friend inline int operator != (cvm::rvector const &v1, cvm::rvector const &v2)
852  {
853  return (v1.x != v2.x) || (v1.y != v2.y) || (v1.z != v2.z);
854  }
855 
856  friend inline cvm::rvector operator + (cvm::rvector const &v1, cvm::rvector const &v2)
857  {
858  return cvm::rvector(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
859  }
860  friend inline cvm::rvector operator - (cvm::rvector const &v1, cvm::rvector const &v2)
861  {
862  return cvm::rvector(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
863  }
864 
865  friend inline cvm::real operator * (cvm::rvector const &v1, cvm::rvector const &v2)
866  {
867  return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
868  }
869 
870  friend inline cvm::rvector operator * (cvm::real const &a, cvm::rvector const &v)
871  {
872  return cvm::rvector(a*v.x, a*v.y, a*v.z);
873  }
874 
875  friend inline cvm::rvector operator * (cvm::rvector const &v, cvm::real const &a)
876  {
877  return cvm::rvector(a*v.x, a*v.y, a*v.z);
878  }
879 
880  friend inline cvm::rvector operator / (cvm::rvector const &v, cvm::real const &a)
881  {
882  return cvm::rvector(v.x/a, v.y/a, v.z/a);
883  }
884 
885  std::string to_simple_string() const;
886  int from_simple_string(std::string const &s);
887 };
888 
889 
893  : public colvarmodule::matrix2d<colvarmodule::real> {
894 private:
895 
896 public:
897 
899  inline cvm::real & xx() { return (*this)[0][0]; }
901  inline cvm::real & xy() { return (*this)[0][1]; }
903  inline cvm::real & xz() { return (*this)[0][2]; }
905  inline cvm::real & yx() { return (*this)[1][0]; }
907  inline cvm::real & yy() { return (*this)[1][1]; }
909  inline cvm::real & yz() { return (*this)[1][2]; }
911  inline cvm::real & zx() { return (*this)[2][0]; }
913  inline cvm::real & zy() { return (*this)[2][1]; }
915  inline cvm::real & zz() { return (*this)[2][2]; }
916 
918  inline cvm::real xx() const { return (*this)[0][0]; }
920  inline cvm::real xy() const { return (*this)[0][1]; }
922  inline cvm::real xz() const { return (*this)[0][2]; }
924  inline cvm::real yx() const { return (*this)[1][0]; }
926  inline cvm::real yy() const { return (*this)[1][1]; }
928  inline cvm::real yz() const { return (*this)[1][2]; }
930  inline cvm::real zx() const { return (*this)[2][0]; }
932  inline cvm::real zy() const { return (*this)[2][1]; }
934  inline cvm::real zz() const { return (*this)[2][2]; }
935 
937  inline rmatrix()
938  : cvm::matrix2d<cvm::real>(3, 3)
939  {}
940 
942  inline rmatrix(cvm::real const &xxi,
943  cvm::real const &xyi,
944  cvm::real const &xzi,
945  cvm::real const &yxi,
946  cvm::real const &yyi,
947  cvm::real const &yzi,
948  cvm::real const &zxi,
949  cvm::real const &zyi,
950  cvm::real const &zzi)
951  : cvm::matrix2d<cvm::real>(3, 3)
952  {
953  this->xx() = xxi;
954  this->xy() = xyi;
955  this->xz() = xzi;
956  this->yx() = yxi;
957  this->yy() = yyi;
958  this->yz() = yzi;
959  this->zx() = zxi;
960  this->zy() = zyi;
961  this->zz() = zzi;
962  }
963 
965  inline ~rmatrix()
966  {}
967 
969  inline cvm::real determinant() const
970  {
971  return
972  ( xx() * (yy()*zz() - zy()*yz()))
973  - (yx() * (xy()*zz() - zy()*xz()))
974  + (zx() * (xy()*yz() - yy()*xz()));
975  }
976 
977  inline cvm::rmatrix transpose() const
978  {
979  return cvm::rmatrix(this->xx(),
980  this->yx(),
981  this->zx(),
982  this->xy(),
983  this->yy(),
984  this->zy(),
985  this->xz(),
986  this->yz(),
987  this->zz());
988  }
989 
990  friend cvm::rvector operator * (cvm::rmatrix const &m, cvm::rvector const &r);
991 
992  // friend inline cvm::rmatrix const operator * (cvm::rmatrix const &m1, cvm::rmatrix const &m2) {
993  // return cvm::rmatrix (m1.xx()*m2.xx() + m1.xy()*m2.yx() + m1.xz()*m2.yz(),
994  // m1.xx()*m2.xy() + m1.xy()*m2.yy() + m1.xz()*m2.zy(),
995  // m1.xx()*m2.xz() + m1.xy()*m2.yz() + m1.xz()*m2.zz(),
996  // m1.yx()*m2.xx() + m1.yy()*m2.yx() + m1.yz()*m2.yz(),
997  // m1.yx()*m2.xy() + m1.yy()*m2.yy() + m1.yz()*m2.yy(),
998  // m1.yx()*m2.xz() + m1.yy()*m2.yz() + m1.yz()*m2.yz(),
999  // m1.zx()*m2.xx() + m1.zy()*m2.yx() + m1.zz()*m2.yz(),
1000  // m1.zx()*m2.xy() + m1.zy()*m2.yy() + m1.zz()*m2.yy(),
1001  // m1.zx()*m2.xz() + m1.zy()*m2.yz() + m1.zz()*m2.yz());
1002  // }
1003 
1004 };
1005 
1006 
1007 inline cvm::rvector operator * (cvm::rmatrix const &m,
1008  cvm::rvector const &r)
1009 {
1010  return cvm::rvector(m.xx()*r.x + m.xy()*r.y + m.xz()*r.z,
1011  m.yx()*r.x + m.yy()*r.y + m.yz()*r.z,
1012  m.zx()*r.x + m.zy()*r.y + m.zz()*r.z);
1013 }
1014 
1015 
1016 
1017 
1021 
1022 public:
1023 
1024  cvm::real q0, q1, q2, q3;
1025 
1027  inline quaternion(cvm::real const &x, cvm::real const &y, cvm::real const &z)
1028  : q0(0.0), q1(x), q2(y), q3(z)
1029  {}
1030 
1032  inline quaternion(cvm::real const qv[4])
1033  : q0(qv[0]), q1(qv[1]), q2(qv[2]), q3(qv[3])
1034  {}
1035 
1037  inline quaternion(cvm::real const &q0i,
1038  cvm::real const &q1i,
1039  cvm::real const &q2i,
1040  cvm::real const &q3i)
1041  : q0(q0i), q1(q1i), q2(q2i), q3(q3i)
1042  {}
1043 
1044  inline quaternion(cvm::vector1d<cvm::real> const &v)
1045  : q0(v[0]), q1(v[1]), q2(v[2]), q3(v[3])
1046  {}
1047 
1051  inline void set_from_euler_angles(cvm::real const &phi_in,
1052  cvm::real const &theta_in,
1053  cvm::real const &psi_in)
1054  {
1055  q0 = ( (std::cos(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::cos(psi_in/2.0)) +
1056  (std::sin(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::sin(psi_in/2.0)) );
1057 
1058  q1 = ( (std::sin(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::cos(psi_in/2.0)) -
1059  (std::cos(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::sin(psi_in/2.0)) );
1060 
1061  q2 = ( (std::cos(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::cos(psi_in/2.0)) +
1062  (std::sin(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::sin(psi_in/2.0)) );
1063 
1064  q3 = ( (std::cos(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::sin(psi_in/2.0)) -
1065  (std::sin(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::cos(psi_in/2.0)) );
1066  }
1067 
1069  inline quaternion()
1070  {
1071  reset();
1072  }
1073 
1075  inline void set(cvm::real const &value = 0.0)
1076  {
1077  q0 = q1 = q2 = q3 = value;
1078  }
1079 
1081  inline void reset()
1082  {
1083  q0 = q1 = q2 = q3 = 0.0;
1084  }
1085 
1088  inline void reset_rotation()
1089  {
1090  q0 = 1.0;
1091  q1 = q2 = q3 = 0.0;
1092  }
1093 
1095  static inline size_t output_width(size_t const &real_width)
1096  {
1097  return 4*real_width + 13;
1098  }
1099 
1100  std::string to_simple_string() const;
1101  int from_simple_string(std::string const &s);
1102 
1104  friend std::ostream & operator << (std::ostream &os, cvm::quaternion const &q);
1106  friend std::istream & operator >> (std::istream &is, cvm::quaternion &q);
1107 
1109  inline cvm::real & operator [] (int const &i) {
1110  switch (i) {
1111  case 0:
1112  return this->q0;
1113  case 1:
1114  return this->q1;
1115  case 2:
1116  return this->q2;
1117  case 3:
1118  return this->q3;
1119  default:
1120  cvm::error("Error: incorrect quaternion component.\n");
1121  return q0;
1122  }
1123  }
1124 
1126  inline cvm::real operator [] (int const &i) const {
1127  switch (i) {
1128  case 0:
1129  return this->q0;
1130  case 1:
1131  return this->q1;
1132  case 2:
1133  return this->q2;
1134  case 3:
1135  return this->q3;
1136  default:
1137  cvm::error("Error: trying to access a quaternion "
1138  "component which is not between 0 and 3.\n");
1139  return 0.0;
1140  }
1141  }
1142 
1143  inline cvm::vector1d<cvm::real> const as_vector() const
1144  {
1145  cvm::vector1d<cvm::real> result(4);
1146  result[0] = q0;
1147  result[1] = q1;
1148  result[2] = q2;
1149  result[3] = q3;
1150  return result;
1151  }
1152 
1154  inline cvm::real norm2() const
1155  {
1156  return q0*q0 + q1*q1 + q2*q2 + q3*q3;
1157  }
1158 
1160  inline cvm::real norm() const
1161  {
1162  return std::sqrt(this->norm2());
1163  }
1164 
1167  {
1168  return cvm::quaternion(q0, -q1, -q2, -q3);
1169  }
1170 
1171  inline void operator *= (cvm::real const &a)
1172  {
1173  q0 *= a; q1 *= a; q2 *= a; q3 *= a;
1174  }
1175 
1176  inline void operator /= (cvm::real const &a)
1177  {
1178  q0 /= a; q1 /= a; q2 /= a; q3 /= a;
1179  }
1180 
1181  inline void set_positive()
1182  {
1183  if (q0 > 0.0) return;
1184  q0 = -q0;
1185  q1 = -q1;
1186  q2 = -q2;
1187  q3 = -q3;
1188  }
1189 
1190  inline void operator += (cvm::quaternion const &h)
1191  {
1192  q0+=h.q0; q1+=h.q1; q2+=h.q2; q3+=h.q3;
1193  }
1194  inline void operator -= (cvm::quaternion const &h)
1195  {
1196  q0-=h.q0; q1-=h.q1; q2-=h.q2; q3-=h.q3;
1197  }
1198 
1200  static inline cvm::quaternion promote(cvm::rvector const &v)
1201  {
1202  return cvm::quaternion(0.0, v.x, v.y, v.z);
1203  }
1205  inline cvm::rvector get_vector() const
1206  {
1207  return cvm::rvector(q1, q2, q3);
1208  }
1209 
1210 
1211  friend inline cvm::quaternion operator + (cvm::quaternion const &h, cvm::quaternion const &q)
1212  {
1213  return cvm::quaternion(h.q0+q.q0, h.q1+q.q1, h.q2+q.q2, h.q3+q.q3);
1214  }
1215 
1216  friend inline cvm::quaternion operator - (cvm::quaternion const &h, cvm::quaternion const &q)
1217  {
1218  return cvm::quaternion(h.q0-q.q0, h.q1-q.q1, h.q2-q.q2, h.q3-q.q3);
1219  }
1220 
1223  friend inline cvm::quaternion operator * (cvm::quaternion const &h, cvm::quaternion const &q)
1224  {
1225  return cvm::quaternion(h.q0*q.q0 - h.q1*q.q1 - h.q2*q.q2 - h.q3*q.q3,
1226  h.q0*q.q1 + h.q1*q.q0 + h.q2*q.q3 - h.q3*q.q2,
1227  h.q0*q.q2 + h.q2*q.q0 + h.q3*q.q1 - h.q1*q.q3,
1228  h.q0*q.q3 + h.q3*q.q0 + h.q1*q.q2 - h.q2*q.q1);
1229  }
1230 
1231  friend inline cvm::quaternion operator * (cvm::real const &c,
1232  cvm::quaternion const &q)
1233  {
1234  return cvm::quaternion(c*q.q0, c*q.q1, c*q.q2, c*q.q3);
1235  }
1236  friend inline cvm::quaternion operator * (cvm::quaternion const &q,
1237  cvm::real const &c)
1238  {
1239  return cvm::quaternion(q.q0*c, q.q1*c, q.q2*c, q.q3*c);
1240  }
1241  friend inline cvm::quaternion operator / (cvm::quaternion const &q,
1242  cvm::real const &c)
1243  {
1244  return cvm::quaternion(q.q0/c, q.q1/c, q.q2/c, q.q3/c);
1245  }
1246 
1247 
1250  inline cvm::rvector rotate(cvm::rvector const &v) const
1251  {
1252  return ((*this) * promote(v) * ((*this).conjugate())).get_vector();
1253  }
1254 
1257  inline cvm::quaternion rotate(cvm::quaternion const &Q2) const
1258  {
1259  cvm::rvector const vq_rot = this->rotate(Q2.get_vector());
1260  return cvm::quaternion(Q2.q0, vq_rot.x, vq_rot.y, vq_rot.z);
1261  }
1262 
1265  {
1266  cvm::rmatrix R;
1267 
1268  R.xx() = q0*q0 + q1*q1 - q2*q2 - q3*q3;
1269  R.yy() = q0*q0 - q1*q1 + q2*q2 - q3*q3;
1270  R.zz() = q0*q0 - q1*q1 - q2*q2 + q3*q3;
1271 
1272  R.xy() = 2.0 * (q1*q2 - q0*q3);
1273  R.xz() = 2.0 * (q0*q2 + q1*q3);
1274 
1275  R.yx() = 2.0 * (q0*q3 + q1*q2);
1276  R.yz() = 2.0 * (q2*q3 - q0*q1);
1277 
1278  R.zx() = 2.0 * (q1*q3 - q0*q2);
1279  R.zy() = 2.0 * (q0*q1 + q2*q3);
1280 
1281  return R;
1282  }
1283 
1284 
1287  cvm::quaternion position_derivative_inner(cvm::rvector const &pos,
1288  cvm::rvector const &vec) const;
1289 
1290 
1293  inline cvm::real cosine(cvm::quaternion const &q) const
1294  {
1295  cvm::real const iprod = this->inner(q);
1296  return 2.0*iprod*iprod - 1.0;
1297  }
1298 
1302  inline cvm::real dist2(cvm::quaternion const &Q2) const
1303  {
1304  cvm::real const cos_omega = this->q0*Q2.q0 + this->q1*Q2.q1 +
1305  this->q2*Q2.q2 + this->q3*Q2.q3;
1306 
1307  cvm::real const omega = std::acos( (cos_omega > 1.0) ? 1.0 :
1308  ( (cos_omega < -1.0) ? -1.0 : cos_omega) );
1309 
1310  // get the minimum distance: x and -x are the same quaternion
1311  if (cos_omega > 0.0)
1312  return omega * omega;
1313  else
1314  return (PI-omega) * (PI-omega);
1315  }
1316 
1320  {
1321  cvm::real const cos_omega = this->q0*Q2.q0 + this->q1*Q2.q1 + this->q2*Q2.q2 + this->q3*Q2.q3;
1322  cvm::real const omega = std::acos( (cos_omega > 1.0) ? 1.0 :
1323  ( (cos_omega < -1.0) ? -1.0 : cos_omega) );
1324  cvm::real const sin_omega = std::sin(omega);
1325 
1326  if (std::fabs(sin_omega) < 1.0E-14) {
1327  // return a null 4d vector
1328  return cvm::quaternion(0.0, 0.0, 0.0, 0.0);
1329  }
1330 
1331  cvm::quaternion const
1332  grad1((-1.0)*sin_omega*Q2.q0 + cos_omega*(this->q0-cos_omega*Q2.q0)/sin_omega,
1333  (-1.0)*sin_omega*Q2.q1 + cos_omega*(this->q1-cos_omega*Q2.q1)/sin_omega,
1334  (-1.0)*sin_omega*Q2.q2 + cos_omega*(this->q2-cos_omega*Q2.q2)/sin_omega,
1335  (-1.0)*sin_omega*Q2.q3 + cos_omega*(this->q3-cos_omega*Q2.q3)/sin_omega);
1336 
1337  if (cos_omega > 0.0) {
1338  return 2.0*omega*grad1;
1339  }
1340  else {
1341  return -2.0*(PI-omega)*grad1;
1342  }
1343  }
1344 
1347  inline void match(cvm::quaternion &Q2) const
1348  {
1349  cvm::real const cos_omega = this->q0*Q2.q0 + this->q1*Q2.q1 +
1350  this->q2*Q2.q2 + this->q3*Q2.q3;
1351  if (cos_omega < 0.0) Q2 *= -1.0;
1352  }
1353 
1356  inline cvm::real inner(cvm::quaternion const &Q2) const
1357  {
1358  cvm::real const prod = this->q0*Q2.q0 + this->q1*Q2.q1 +
1359  this->q2*Q2.q2 + this->q3*Q2.q3;
1360  return prod;
1361  }
1362 
1363 
1364 };
1365 
1366 
1370 {
1371 public:
1372 
1375 
1378 
1381 
1384  std::vector<cvm::atom_pos> pos1, pos2;
1385 
1386  cvm::rmatrix C;
1387 
1389  cvm::vector1d<cvm::real> S_eigval;
1390  cvm::matrix2d<cvm::real> S_eigvec;
1391 
1394 
1396  std::vector< cvm::matrix2d<cvm::rvector> > dS_1, dS_2;
1398  std::vector< cvm::rvector > dL0_1, dL0_2;
1400  std::vector< cvm::vector1d<cvm::rvector> > dQ0_1, dQ0_2;
1401 
1403  inline void request_group1_gradients(size_t const &n)
1404  {
1405  dS_1.resize(n, cvm::matrix2d<cvm::rvector>(4, 4));
1406  dL0_1.resize(n, cvm::rvector(0.0, 0.0, 0.0));
1407  dQ0_1.resize(n, cvm::vector1d<cvm::rvector>(4));
1408  }
1409 
1411  inline void request_group2_gradients(size_t const &n)
1412  {
1413  dS_2.resize(n, cvm::matrix2d<cvm::rvector>(4, 4));
1414  dL0_2.resize(n, cvm::rvector(0.0, 0.0, 0.0));
1415  dQ0_2.resize(n, cvm::vector1d<cvm::rvector>(4));
1416  }
1417 
1428  void calc_optimal_rotation(std::vector<atom_pos> const &pos1,
1429  std::vector<atom_pos> const &pos2);
1430 
1432  inline rotation()
1433  : b_debug_gradients(false)
1434  {}
1435 
1437  inline rotation(cvm::quaternion const &qi)
1438  : q(qi),
1439  b_debug_gradients(false)
1440  {
1441  }
1442 
1444  inline rotation(cvm::real const &angle, cvm::rvector const &axis)
1445  : b_debug_gradients(false)
1446  {
1447  cvm::rvector const axis_n = axis.unit();
1448  cvm::real const sina = std::sin(angle/2.0);
1449  q = cvm::quaternion(std::cos(angle/2.0),
1450  sina * axis_n.x, sina * axis_n.y, sina * axis_n.z);
1451  }
1452 
1454  inline ~rotation()
1455  {}
1456 
1458  inline cvm::rvector rotate(cvm::rvector const &v) const
1459  {
1460  return q.rotate(v);
1461  }
1462 
1464  inline cvm::rotation inverse() const
1465  {
1466  return cvm::rotation(this->q.conjugate());
1467  }
1468 
1470  inline cvm::rmatrix matrix() const
1471  {
1472  return q.rotation_matrix();
1473  }
1474 
1475 
1478  inline cvm::real spin_angle(cvm::rvector const &axis) const
1479  {
1480  cvm::rvector const q_vec = q.get_vector();
1481  cvm::real alpha = (180.0/PI) * 2.0 * std::atan2(axis * q_vec, q.q0);
1482  while (alpha > 180.0) alpha -= 360;
1483  while (alpha < -180.0) alpha += 360;
1484  return alpha;
1485  }
1486 
1489  inline cvm::quaternion dspin_angle_dq(cvm::rvector const &axis) const
1490  {
1491  cvm::rvector const q_vec = q.get_vector();
1492  cvm::real const iprod = axis * q_vec;
1493 
1494  if (q.q0 != 0.0) {
1495 
1496  // cvm::real const x = iprod/q.q0;
1497 
1498  cvm::real const dspindx = (180.0/PI) * 2.0 * (1.0 / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
1499 
1500  return
1501  cvm::quaternion( dspindx * (iprod * (-1.0) / (q.q0*q.q0)),
1502  dspindx * ((1.0/q.q0) * axis.x),
1503  dspindx * ((1.0/q.q0) * axis.y),
1504  dspindx * ((1.0/q.q0) * axis.z));
1505  } else {
1506  // (1/(1+x^2)) ~ (1/x)^2
1507  return
1508  cvm::quaternion((180.0/PI) * 2.0 * ((-1.0)/iprod), 0.0, 0.0, 0.0);
1509  // XX TODO: What if iprod == 0? XX
1510  }
1511  }
1512 
1515  inline cvm::real cos_theta(cvm::rvector const &axis) const
1516  {
1517  cvm::rvector const q_vec = q.get_vector();
1518  cvm::real const alpha =
1519  (180.0/PI) * 2.0 * std::atan2(axis * q_vec, q.q0);
1520 
1521  cvm::real const cos_spin_2 = std::cos(alpha * (PI/180.0) * 0.5);
1522  cvm::real const cos_theta_2 = ( (cos_spin_2 != 0.0) ?
1523  (q.q0 / cos_spin_2) :
1524  (0.0) );
1525  // cos(2t) = 2*cos(t)^2 - 1
1526  return 2.0 * (cos_theta_2*cos_theta_2) - 1.0;
1527  }
1528 
1530  inline cvm::quaternion dcos_theta_dq(cvm::rvector const &axis) const
1531  {
1532  cvm::rvector const q_vec = q.get_vector();
1533  cvm::real const iprod = axis * q_vec;
1534 
1535  cvm::real const cos_spin_2 = std::cos(std::atan2(iprod, q.q0));
1536 
1537  if (q.q0 != 0.0) {
1538 
1539  cvm::real const d_cos_theta_dq0 =
1540  (4.0 * q.q0 / (cos_spin_2*cos_spin_2)) *
1541  (1.0 - (iprod*iprod)/(q.q0*q.q0) / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
1542 
1543  cvm::real const d_cos_theta_dqn =
1544  (4.0 * q.q0 / (cos_spin_2*cos_spin_2) *
1545  (iprod/q.q0) / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
1546 
1547  return cvm::quaternion(d_cos_theta_dq0,
1548  d_cos_theta_dqn * axis.x,
1549  d_cos_theta_dqn * axis.y,
1550  d_cos_theta_dqn * axis.z);
1551  } else {
1552 
1553  cvm::real const d_cos_theta_dqn =
1554  (4.0 / (cos_spin_2*cos_spin_2 * iprod));
1555 
1556  return cvm::quaternion(0.0,
1557  d_cos_theta_dqn * axis.x,
1558  d_cos_theta_dqn * axis.y,
1559  d_cos_theta_dqn * axis.z);
1560  }
1561  }
1562 
1564  static bool monitor_crossings;
1567 
1568 protected:
1569 
1574 
1576  void build_matrix(std::vector<cvm::atom_pos> const &pos1,
1577  std::vector<cvm::atom_pos> const &pos2,
1579 
1581  void diagonalize_matrix(cvm::matrix2d<cvm::real> &S,
1582  cvm::vector1d<cvm::real> &S_eigval,
1583  cvm::matrix2d<cvm::real> &S_eigvec);
1584 };
1585 
1586 
1587 #endif
cvm::real yy() const
Return the yy element.
Definition: colvartypes.h:926
1-dimensional vector of real numbers with four components and a quaternion algebra ...
Definition: colvartypes.h:1020
cvm::rmatrix rotation_matrix() const
Return the 3x3 matrix associated to this quaternion.
Definition: colvartypes.h:1264
void request_group1_gradients(size_t const &n)
Allocate space for the derivatives of the rotation.
Definition: colvartypes.h:1403
Definition: colvartypes.h:356
cvm::real inner(cvm::quaternion const &Q2) const
Inner product (as a 4-d vector) with Q2; requires match() if the largest overlap is looked for...
Definition: colvartypes.h:1356
cvm::real & xz()
Return the xz element.
Definition: colvartypes.h:903
~matrix2d()
Destructor.
Definition: colvartypes.h:466
quaternion()
Default constructor.
Definition: colvartypes.h:1069
cvm::quaternion q_old
Previous value of the rotation (used to warn the user when the structure changes too much...
Definition: colvartypes.h:1573
rotation(cvm::quaternion const &qi)
Constructor after a quaternion.
Definition: colvartypes.h:1437
static bool monitor_crossings
Whether to test for eigenvalue crossing.
Definition: colvartypes.h:1564
cvm::real & zx()
Return the zx element.
Definition: colvartypes.h:911
cvm::rvector rotate(cvm::rvector const &v) const
Return the rotated vector.
Definition: colvartypes.h:1458
cvm::real norm2() const
Square norm of the quaternion.
Definition: colvartypes.h:1154
cvm::real & zz()
Return the zz element.
Definition: colvartypes.h:915
quaternion(cvm::real const &x, cvm::real const &y, cvm::real const &z)
Constructor from a 3-d vector.
Definition: colvartypes.h:1027
static cvm::quaternion promote(cvm::rvector const &v)
Promote a 3-vector to a quaternion.
Definition: colvartypes.h:1200
rotation(cvm::real const &angle, cvm::rvector const &axis)
Constructor after an axis of rotation and an angle (in radians)
Definition: colvartypes.h:1444
cvm::quaternion dist2_grad(cvm::quaternion const &Q2) const
Definition: colvartypes.h:1319
cvm::real & xy()
Return the xy element.
Definition: colvartypes.h:901
std::vector< T > & data_array()
Return a reference to the data.
Definition: colvartypes.h:471
void clear()
Deallocation routine.
Definition: colvartypes.h:425
cvm::real lambda
Eigenvalue corresponding to the optimal rotation.
Definition: colvartypes.h:1377
std::vector< T > & data_array()
Return a reference to the data.
Definition: colvartypes.h:61
size_t output_width(size_t const &real_width) const
Formatted output.
Definition: colvartypes.h:255
rmatrix(cvm::real const &xxi, cvm::real const &xyi, cvm::real const &xzi, cvm::real const &yxi, cvm::real const &yyi, cvm::real const &yzi, cvm::real const &zxi, cvm::real const &zyi, cvm::real const &zzi)
Constructor component by component.
Definition: colvartypes.h:942
cvm::real cosine(cvm::quaternion const &q) const
Return the cosine between the orientation frame associated to this quaternion and another...
Definition: colvartypes.h:1293
std::vector< cvm::vector1d< cvm::rvector > > dQ0_1
Derivatives of leading eigenvector.
Definition: colvartypes.h:1400
vector of real numbers with three components
Definition: colvartypes.h:714
cvm::real & zy()
Return the zy element.
Definition: colvartypes.h:913
cvm::real cos_theta(cvm::rvector const &axis) const
Return the projection of the orientation vector onto a predefined axis.
Definition: colvartypes.h:1515
static size_t output_width(size_t const &real_width)
Tell the number of characters required to print a quaternion, given that of a real number...
Definition: colvartypes.h:1095
Collective variables module (main class)
Definition: colvarmodule.h:62
std::vector< cvm::rvector > dL0_1
Derivatives of leading eigenvalue.
Definition: colvartypes.h:1398
cvm::real & yx()
Return the yx element.
Definition: colvartypes.h:905
T * c_array()
Return a pointer to the data location.
Definition: colvartypes.h:51
void request_group2_gradients(size_t const &n)
Allocate space for the derivatives of the rotation.
Definition: colvartypes.h:1411
rmatrix()
Default constructor.
Definition: colvartypes.h:937
Collective variables main module.
double real
Defining an abstract real number allows to switch precision.
Definition: colvarmodule.h:76
static int error(std::string const &message, int code=COLVARS_ERROR)
Print a message to the main log and set global error code.
Definition: colvarmodule.cpp:1587
cvm::real norm2() const
Squared norm.
Definition: colvartypes.h:202
vector1d(size_t const n=0)
Default constructor.
Definition: colvartypes.h:33
~rmatrix()
Destructor.
Definition: colvartypes.h:965
matrix2d(matrix2d< T > const &m)
Copy constructor.
Definition: colvartypes.h:456
cvm::rmatrix matrix() const
Return the associated 3x3 matrix.
Definition: colvartypes.h:1470
std::vector< cvm::atom_pos > pos1
Positions to superimpose: the rotation should brings pos1 into pos2.
Definition: colvartypes.h:1384
vector1d< T > const slice(size_t const i1, size_t const i2) const
Slicing.
Definition: colvartypes.h:228
cvm::real xx() const
Return the xx element.
Definition: colvartypes.h:918
cvm::real zy() const
Return the zy element.
Definition: colvartypes.h:932
quaternion(cvm::real const qv[4])
Constructor component by component.
Definition: colvartypes.h:1032
2-dimensional array of real numbers with three components along each dimension (works with colvarmodu...
Definition: colvartypes.h:892
cvm::rvector get_vector() const
Return the vector component.
Definition: colvartypes.h:1205
cvm::real yz() const
Return the yz element.
Definition: colvartypes.h:928
void reset_rotation()
Set the q0 component to 1 and the others to 0 (quaternion representing no rotation) ...
Definition: colvartypes.h:1088
cvm::real xy() const
Return the xy element.
Definition: colvartypes.h:920
cvm::real norm() const
Norm of the quaternion.
Definition: colvartypes.h:1160
Arbitrary size array (two dimensions) suitable for linear algebra operations (i.e. for floating point numbers it can be used with library functions)
Definition: colvarmodule.h:97
cvm::matrix2d< cvm::real > S_backup
Used for debugging gradients.
Definition: colvartypes.h:1393
cvm::real & xx()
Return the xx element.
Definition: colvartypes.h:899
cvm::real spin_angle(cvm::rvector const &axis) const
Return the spin angle (in degrees) with respect to the provided axis (which MUST be normalized) ...
Definition: colvartypes.h:1478
cvm::quaternion dspin_angle_dq(cvm::rvector const &axis) const
Return the derivative of the spin angle with respect to the quaternion.
Definition: colvartypes.h:1489
void reset()
Set all elements to zero.
Definition: colvartypes.h:431
void match(cvm::quaternion &Q2) const
Choose the closest between Q2 and -Q2 and save it back. Not required for dist2() and dist2_grad() ...
Definition: colvartypes.h:1347
cvm::rvector rotate(cvm::rvector const &v) const
Rotate v through this quaternion (put it in the rotated reference frame)
Definition: colvartypes.h:1250
void set_from_euler_angles(cvm::real const &phi_in, cvm::real const &theta_in, cvm::real const &psi_in)
Definition: colvartypes.h:1051
~rotation()
Destructor.
Definition: colvartypes.h:1454
T ** c_array()
Return the 2-d C array.
Definition: colvartypes.h:499
cvm::quaternion conjugate() const
Return the conjugate quaternion.
Definition: colvartypes.h:1166
matrix2d()
Default constructor.
Definition: colvartypes.h:442
static std::string to_str(T const &x, size_t const &width=0, size_t const &prec=0)
Quick conversion of an object to a string.
Definition: colvarmodule.h:607
void sliceassign(size_t const i1, size_t const i2, vector1d< T > const &v)
Assign a vector to a slice of this vector.
Definition: colvartypes.h:242
void resize(size_t const ol, size_t const il)
Allocation routine, used by all constructors.
Definition: colvartypes.h:382
std::vector< cvm::matrix2d< cvm::rvector > > dS_1
Derivatives of S.
Definition: colvartypes.h:1396
cvm::rotation inverse() const
Return the inverse of this rotation.
Definition: colvartypes.h:1464
cvm::real zx() const
Return the zx element.
Definition: colvartypes.h:930
cvm::real zz() const
Return the zz element.
Definition: colvartypes.h:934
vector1d(size_t const n, T const *t)
Constructor from C array.
Definition: colvartypes.h:40
cvm::real & yy()
Return the yy element.
Definition: colvartypes.h:907
static size_t const cv_prec
Number of digits to represent a collective variables value(s)
Definition: colvarmodule.h:427
cvm::real & yz()
Return the yz element.
Definition: colvartypes.h:909
void reset()
Set all components to zero.
Definition: colvartypes.h:753
cvm::real dist2(cvm::quaternion const &Q2) const
Square distance from another quaternion on the 4-dimensional unit sphere: returns the square of the a...
Definition: colvartypes.h:1302
cvm::quaternion dcos_theta_dq(cvm::rvector const &axis) const
Return the derivative of the tilt wrt the quaternion.
Definition: colvartypes.h:1530
bool b_debug_gradients
Perform gradient tests.
Definition: colvartypes.h:1380
void reset()
Set all components to zero (null quaternion)
Definition: colvartypes.h:1081
A rotation between two sets of coordinates (for the moment a wrapper for colvarmodule::quaternion) ...
Definition: colvartypes.h:1369
cvm::real xz() const
Return the xz element.
Definition: colvartypes.h:922
cvm::quaternion rotate(cvm::quaternion const &Q2) const
Rotate Q2 through this quaternion (put it in the rotated reference frame)
Definition: colvartypes.h:1257
quaternion(cvm::real const &q0i, cvm::real const &q1i, cvm::real const &q2i, cvm::real const &q3i)
Constructor component by component.
Definition: colvartypes.h:1037
cvm::quaternion q
The rotation itself (implemented as a quaternion)
Definition: colvartypes.h:1374
void reset()
Set all elements to zero.
Definition: colvartypes.h:72
static cvm::real crossing_threshold
Threshold for the eigenvalue crossing test.
Definition: colvartypes.h:1566
Arbitrary size array (one dimensions) suitable for linear algebra operations (i.e. for floating point numbers it can be used with library functions)
Definition: colvarmodule.h:96
cvm::real determinant() const
Return the determinant.
Definition: colvartypes.h:969
rotation()
Default constructor.
Definition: colvartypes.h:1432
cvm::real yx() const
Return the yx element.
Definition: colvartypes.h:924