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