← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/colvars/colvarvalue.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/colvars/colvarvalue.cpp
 
1
/// -*- c++ -*-
 
2
 
1
3
#include <vector>
 
4
#include <sstream>
 
5
#include <iostream>
2
6
 
3
7
#include "colvarmodule.h"
4
8
#include "colvarvalue.h"
5
9
 
6
10
 
7
11
 
8
 
std::string const colvarvalue::type_desc[colvarvalue::type_quaternion+1] =
9
 
  { "not_set",
10
 
    "scalar number",
11
 
    "3-dimensional vector",
12
 
    "3-dimensional unit vector",
13
 
    "4-dimensional unit vector" };
14
 
 
15
 
size_t const      colvarvalue::dof_num[  colvarvalue::type_quaternion+1] =
16
 
  { 0, 1, 3, 2, 3 };
17
 
 
18
 
 
19
 
void colvarvalue::undef_op() const
20
 
{
21
 
  cvm::fatal_error ("Error: Undefined operation on a colvar of type \""+
22
 
                    colvarvalue::type_desc[this->value_type]+"\".\n");
23
 
}
24
 
 
25
 
void colvarvalue::error_rside
26
 
(colvarvalue::Type const &vt) const
27
 
{
28
 
  cvm::fatal_error ("Trying to assign a colvar value with type \""+
29
 
                    type_desc[this->value_type]+"\" to one with type \""+
30
 
                    type_desc[vt]+"\".\n");
31
 
}
32
 
 
33
 
void colvarvalue::error_lside
34
 
(colvarvalue::Type const &vt) const
35
 
{
36
 
  cvm::fatal_error ("Trying to use a colvar value with type \""+
37
 
                    type_desc[vt]+"\" as one of type \""+
38
 
                    type_desc[this->value_type]+"\".\n");
39
 
}
40
 
 
41
 
 
42
 
 
43
 
void colvarvalue::inner_opt (colvarvalue                        const &x,
44
 
                             std::vector<colvarvalue>::iterator       &xv,
45
 
                             std::vector<colvarvalue>::iterator const &xv_end,
46
 
                             std::vector<cvm::real>::iterator         &inner)
47
 
{
48
 
  // doing type check only once, here
49
 
  colvarvalue::check_types (x, *xv);
50
 
 
51
 
  std::vector<colvarvalue>::iterator &xvi = xv;
52
 
  std::vector<cvm::real>::iterator    &ii = inner;
53
 
 
54
 
  switch (x.value_type) {
55
 
  case colvarvalue::type_scalar:
56
 
    while (xvi != xv_end) {
57
 
      *(ii++) += (xvi++)->real_value * x.real_value;
58
 
    }
59
 
    break;
60
 
  case colvarvalue::type_vector:
61
 
  case colvarvalue::type_unitvector:
62
 
    while (xvi != xv_end) {
63
 
      *(ii++) += (xvi++)->rvector_value * x.rvector_value;
64
 
    }
65
 
    break;
66
 
  case colvarvalue::type_quaternion:
67
 
    while (xvi != xv_end) {
68
 
      *(ii++) += ((xvi++)->quaternion_value).cosine (x.quaternion_value);
69
 
    }
70
 
    break;
71
 
  default:
72
 
    x.undef_op();
73
 
  };
74
 
}
75
 
 
76
 
void colvarvalue::inner_opt (colvarvalue const                      &x,
77
 
                             std::list<colvarvalue>::iterator       &xv,
78
 
                             std::list<colvarvalue>::iterator const &xv_end,
79
 
                             std::vector<cvm::real>::iterator       &inner)
80
 
{
81
 
  // doing type check only once, here
82
 
  colvarvalue::check_types (x, *xv);
83
 
 
84
 
  std::list<colvarvalue>::iterator &xvi = xv;
85
 
  std::vector<cvm::real>::iterator  &ii = inner;
86
 
 
87
 
  switch (x.value_type) {
88
 
  case colvarvalue::type_scalar:
89
 
    while (xvi != xv_end) {
90
 
      *(ii++) += (xvi++)->real_value * x.real_value;
91
 
    }
92
 
    break;
93
 
  case colvarvalue::type_vector:
94
 
  case colvarvalue::type_unitvector:
95
 
    while (xvi != xv_end) {
96
 
      *(ii++) += (xvi++)->rvector_value * x.rvector_value;
97
 
    }
98
 
    break;
99
 
  case colvarvalue::type_quaternion:
100
 
    while (xvi != xv_end) {
101
 
      *(ii++) += ((xvi++)->quaternion_value).cosine (x.quaternion_value);
102
 
    }
103
 
    break;
104
 
  default:
105
 
    x.undef_op();
106
 
  };
107
 
}
108
 
 
109
 
 
110
 
void colvarvalue::p2leg_opt (colvarvalue const                        &x,
111
 
                             std::vector<colvarvalue>::iterator       &xv,
112
 
                             std::vector<colvarvalue>::iterator const &xv_end,
113
 
                             std::vector<cvm::real>::iterator         &inner)
114
 
{
115
 
  // doing type check only once, here
116
 
  colvarvalue::check_types (x, *xv);
117
 
 
118
 
  std::vector<colvarvalue>::iterator &xvi = xv;
119
 
  std::vector<cvm::real>::iterator    &ii = inner;
120
 
 
121
 
  switch (x.value_type) {
122
 
  case colvarvalue::type_scalar:
123
 
    cvm::fatal_error ("Error: cannot calculate Legendre polynomials "
124
 
                      "for scalar variables.\n");
125
 
    break;
126
 
  case colvarvalue::type_vector:
127
 
    while (xvi != xv_end) {
128
 
      cvm::real const cosine =
129
 
        ((xvi)->rvector_value * x.rvector_value) /
130
 
        ((xvi)->rvector_value.norm() * x.rvector_value.norm());
131
 
      xvi++;
132
 
      *(ii++) += 1.5*cosine*cosine - 0.5;
133
 
    }
134
 
    break;
135
 
  case colvarvalue::type_unitvector:
136
 
    while (xvi != xv_end) {
137
 
      cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value;
138
 
      *(ii++) += 1.5*cosine*cosine - 0.5;
139
 
    }
140
 
    break;
141
 
  case colvarvalue::type_quaternion:
142
 
    while (xvi != xv_end) {
143
 
      cvm::real const cosine = (xvi++)->quaternion_value.cosine (x.quaternion_value);
144
 
      *(ii++) += 1.5*cosine*cosine - 0.5;
145
 
    }
146
 
    break;
147
 
  default:
148
 
    x.undef_op();
149
 
  };
150
 
}
151
 
 
152
 
void colvarvalue::p2leg_opt (colvarvalue const                        &x,
153
 
                             std::list<colvarvalue>::iterator         &xv,
154
 
                             std::list<colvarvalue>::iterator const   &xv_end,
155
 
                             std::vector<cvm::real>::iterator         &inner)
156
 
{
157
 
  // doing type check only once, here
158
 
  colvarvalue::check_types (x, *xv);
159
 
 
160
 
  std::list<colvarvalue>::iterator &xvi = xv;
161
 
  std::vector<cvm::real>::iterator  &ii = inner;
162
 
 
163
 
  switch (x.value_type) {
164
 
  case colvarvalue::type_scalar:
165
 
    cvm::fatal_error ("Error: cannot calculate Legendre polynomials "
166
 
                      "for scalar variables.\n");
167
 
    break;
168
 
  case colvarvalue::type_vector:
169
 
    while (xvi != xv_end) {
170
 
      cvm::real const cosine =
171
 
        ((xvi)->rvector_value * x.rvector_value) /
172
 
        ((xvi)->rvector_value.norm() * x.rvector_value.norm());
173
 
      xvi++;
174
 
      *(ii++) += 1.5*cosine*cosine - 0.5;
175
 
    }
176
 
    break;
177
 
  case colvarvalue::type_unitvector:
178
 
    while (xvi != xv_end) {
179
 
      cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value;
180
 
      *(ii++) += 1.5*cosine*cosine - 0.5;
181
 
    }
182
 
    break;
183
 
  case colvarvalue::type_quaternion:
184
 
    while (xvi != xv_end) {
185
 
      cvm::real const cosine = (xvi++)->quaternion_value.cosine (x.quaternion_value);
186
 
      *(ii++) += 1.5*cosine*cosine - 0.5;
187
 
    }
188
 
    break;
189
 
  default:
190
 
    x.undef_op();
191
 
  };
192
 
}
193
 
 
 
12
void colvarvalue::add_elem(colvarvalue const &x)
 
13
{
 
14
  if (this->value_type != type_vector) {
 
15
    cvm::error("Error: trying to set an element for a variable that is not set to be a vector.\n");
 
16
    return;
 
17
  }
 
18
  size_t const n = vector1d_value.size();
 
19
  size_t const nd = num_dimensions(x.value_type);
 
20
  elem_types.push_back(x.value_type);
 
21
  elem_indices.push_back(n);
 
22
  elem_sizes.push_back(nd);
 
23
  vector1d_value.resize(n + nd);
 
24
  set_elem(n, x);
 
25
}
 
26
 
 
27
 
 
28
colvarvalue const colvarvalue::get_elem(int const i_begin, int const i_end, Type const vt) const
 
29
{
 
30
  if (vector1d_value.size() > 0) {
 
31
    cvm::vector1d<cvm::real> const v(vector1d_value.slice(i_begin, i_end));
 
32
    return colvarvalue(v, vt);
 
33
  } else {
 
34
    cvm::error("Error: trying to get an element from a variable that is not a vector.\n");
 
35
    return colvarvalue(type_notset);
 
36
  }
 
37
}
 
38
 
 
39
 
 
40
void colvarvalue::set_elem(int const i_begin, int const i_end, colvarvalue const &x)
 
41
{
 
42
  if (vector1d_value.size() > 0) {
 
43
    vector1d_value.sliceassign(i_begin, i_end, x.as_vector());
 
44
  } else {
 
45
    cvm::error("Error: trying to set an element for a variable that is not a vector.\n");
 
46
  }
 
47
}
 
48
 
 
49
 
 
50
colvarvalue const colvarvalue::get_elem(int const icv) const
 
51
{
 
52
  if (elem_types.size() > 0) {
 
53
    return get_elem(elem_indices[icv], elem_indices[icv] + elem_sizes[icv],
 
54
                    elem_types[icv]);
 
55
  } else {
 
56
    cvm::error("Error: trying to get a colvarvalue element from a vector colvarvalue that was initialized as a plain array.\n");
 
57
    return colvarvalue(type_notset);
 
58
  }
 
59
}
 
60
 
 
61
 
 
62
void colvarvalue::set_elem(int const icv, colvarvalue const &x)
 
63
{
 
64
  if (elem_types.size() > 0) {
 
65
    check_types_assign(elem_types[icv], x.value_type);
 
66
    set_elem(elem_indices[icv], elem_indices[icv] + elem_sizes[icv], x);
 
67
  } else {
 
68
    cvm::error("Error: trying to set a colvarvalue element for a colvarvalue that was initialized as a plain array.\n");
 
69
  }
 
70
}
 
71
 
 
72
 
 
73
colvarvalue colvarvalue::inverse() const
 
74
{
 
75
  switch (value_type) {
 
76
  case colvarvalue::type_scalar:
 
77
    return colvarvalue(1.0/real_value);
 
78
    break;
 
79
  case colvarvalue::type_3vector:
 
80
  case colvarvalue::type_unit3vector:
 
81
  case colvarvalue::type_unit3vectorderiv:
 
82
    return colvarvalue(cvm::rvector(1.0/rvector_value.x,
 
83
                                    1.0/rvector_value.y,
 
84
                                    1.0/rvector_value.z));
 
85
    break;
 
86
  case colvarvalue::type_quaternion:
 
87
  case colvarvalue::type_quaternionderiv:
 
88
    return colvarvalue(cvm::quaternion(1.0/quaternion_value.q0,
 
89
                                       1.0/quaternion_value.q1,
 
90
                                       1.0/quaternion_value.q2,
 
91
                                       1.0/quaternion_value.q3));
 
92
    break;
 
93
  case colvarvalue::type_vector:
 
94
    {
 
95
      cvm::vector1d<cvm::real> result(vector1d_value);
 
96
      if (elem_types.size() > 0) {
 
97
        // if we have information about non-scalar types, use it
 
98
        size_t i;
 
99
        for (i = 0; i < elem_types.size(); i++) {
 
100
          result.sliceassign(elem_indices[i], elem_indices[i]+elem_sizes[i],
 
101
                             cvm::vector1d<cvm::real>((this->get_elem(i)).inverse()));
 
102
        }
 
103
      } else {
 
104
        size_t i;
 
105
        for (i = 0; i < result.size(); i++) {
 
106
          if (result[i] != 0.0) {
 
107
            result = 1.0/result[i];
 
108
          }
 
109
        }
 
110
      }
 
111
      return colvarvalue(result, type_vector);
 
112
    }
 
113
    break;
 
114
  case colvarvalue::type_notset:
 
115
  default:
 
116
    undef_op();
 
117
    break;
 
118
  }
 
119
  return colvarvalue();
 
120
}
 
121
 
 
122
 
 
123
// binary operations between two colvarvalues
 
124
 
 
125
colvarvalue operator + (colvarvalue const &x1,
 
126
                        colvarvalue const &x2)
 
127
{
 
128
  colvarvalue::check_types(x1, x2);
 
129
 
 
130
  switch (x1.value_type) {
 
131
  case colvarvalue::type_scalar:
 
132
    return colvarvalue(x1.real_value + x2.real_value);
 
133
  case colvarvalue::type_3vector:
 
134
    return colvarvalue(x1.rvector_value + x2.rvector_value);
 
135
  case colvarvalue::type_unit3vector:
 
136
  case colvarvalue::type_unit3vectorderiv:
 
137
    return colvarvalue(x1.rvector_value + x2.rvector_value,
 
138
                       colvarvalue::type_unit3vector);
 
139
  case colvarvalue::type_quaternion:
 
140
  case colvarvalue::type_quaternionderiv:
 
141
    return colvarvalue(x1.quaternion_value + x2.quaternion_value);
 
142
  case colvarvalue::type_vector:
 
143
    return colvarvalue(x1.vector1d_value + x2.vector1d_value, colvarvalue::type_vector);
 
144
  case colvarvalue::type_notset:
 
145
  default:
 
146
    x1.undef_op();
 
147
    return colvarvalue(colvarvalue::type_notset);
 
148
  };
 
149
}
 
150
 
 
151
 
 
152
colvarvalue operator - (colvarvalue const &x1,
 
153
                        colvarvalue const &x2)
 
154
{
 
155
  colvarvalue::check_types(x1, x2);
 
156
 
 
157
  switch (x1.value_type) {
 
158
  case colvarvalue::type_scalar:
 
159
    return colvarvalue(x1.real_value - x2.real_value);
 
160
  case colvarvalue::type_3vector:
 
161
    return colvarvalue(x1.rvector_value - x2.rvector_value);
 
162
  case colvarvalue::type_unit3vector:
 
163
  case colvarvalue::type_unit3vectorderiv:
 
164
    return colvarvalue(x1.rvector_value - x2.rvector_value,
 
165
                       colvarvalue::type_unit3vector);
 
166
  case colvarvalue::type_quaternion:
 
167
  case colvarvalue::type_quaternionderiv:
 
168
    return colvarvalue(x1.quaternion_value - x2.quaternion_value);
 
169
  case colvarvalue::type_vector:
 
170
    return colvarvalue(x1.vector1d_value - x2.vector1d_value, colvarvalue::type_vector);
 
171
  case colvarvalue::type_notset:
 
172
  default:
 
173
    x1.undef_op();
 
174
    return colvarvalue(colvarvalue::type_notset);
 
175
  };
 
176
}
 
177
 
 
178
 
 
179
// binary operations with real numbers
 
180
 
 
181
colvarvalue operator * (cvm::real const &a,
 
182
                        colvarvalue const &x)
 
183
{
 
184
  switch (x.value_type) {
 
185
  case colvarvalue::type_scalar:
 
186
    return colvarvalue(a * x.real_value);
 
187
  case colvarvalue::type_3vector:
 
188
    return colvarvalue(a * x.rvector_value);
 
189
  case colvarvalue::type_unit3vector:
 
190
  case colvarvalue::type_unit3vectorderiv:
 
191
    return colvarvalue(a * x.rvector_value,
 
192
                       colvarvalue::type_unit3vector);
 
193
  case colvarvalue::type_quaternion:
 
194
  case colvarvalue::type_quaternionderiv:
 
195
    return colvarvalue(a * x.quaternion_value);
 
196
  case colvarvalue::type_vector:
 
197
    return colvarvalue(x.vector1d_value * a, colvarvalue::type_vector);
 
198
  case colvarvalue::type_notset:
 
199
  default:
 
200
    x.undef_op();
 
201
    return colvarvalue(colvarvalue::type_notset);
 
202
  }
 
203
}
 
204
 
 
205
 
 
206
colvarvalue operator * (colvarvalue const &x,
 
207
                        cvm::real const &a)
 
208
{
 
209
  return a * x;
 
210
}
 
211
 
 
212
 
 
213
colvarvalue operator / (colvarvalue const &x,
 
214
                        cvm::real const &a)
 
215
{
 
216
  switch (x.value_type) {
 
217
  case colvarvalue::type_scalar:
 
218
    return colvarvalue(x.real_value / a);
 
219
  case colvarvalue::type_3vector:
 
220
    return colvarvalue(x.rvector_value / a);
 
221
  case colvarvalue::type_unit3vector:
 
222
  case colvarvalue::type_unit3vectorderiv:
 
223
    return colvarvalue(x.rvector_value / a,
 
224
                       colvarvalue::type_unit3vector);
 
225
  case colvarvalue::type_quaternion:
 
226
  case colvarvalue::type_quaternionderiv:
 
227
    return colvarvalue(x.quaternion_value / a);
 
228
  case colvarvalue::type_vector:
 
229
    return colvarvalue(x.vector1d_value / a, colvarvalue::type_vector);
 
230
  case colvarvalue::type_notset:
 
231
  default:
 
232
    x.undef_op();
 
233
    return colvarvalue(colvarvalue::type_notset);
 
234
  }
 
235
}
 
236
 
 
237
 
 
238
// inner product between two colvarvalues
 
239
 
 
240
cvm::real operator * (colvarvalue const &x1,
 
241
                      colvarvalue const &x2)
 
242
{
 
243
  colvarvalue::check_types(x1, x2);
 
244
 
 
245
  switch (x1.value_type) {
 
246
  case colvarvalue::type_scalar:
 
247
    return (x1.real_value * x2.real_value);
 
248
  case colvarvalue::type_3vector:
 
249
  case colvarvalue::type_unit3vector:
 
250
  case colvarvalue::type_unit3vectorderiv:
 
251
    return (x1.rvector_value * x2.rvector_value);
 
252
  case colvarvalue::type_quaternion:
 
253
  case colvarvalue::type_quaternionderiv:
 
254
    // the "*" product is the quaternion product, here the inner
 
255
    // member function is used instead
 
256
    return (x1.quaternion_value.inner(x2.quaternion_value));
 
257
  case colvarvalue::type_vector:
 
258
    return (x1.vector1d_value * x2.vector1d_value);
 
259
  case colvarvalue::type_notset:
 
260
  default:
 
261
    x1.undef_op();
 
262
    return 0.0;
 
263
  };
 
264
}
 
265
 
 
266
 
 
267
colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
 
268
{
 
269
  colvarvalue::check_types(*this, x2);
 
270
 
 
271
  switch (this->value_type) {
 
272
  case colvarvalue::type_scalar:
 
273
    return 2.0 * (this->real_value - x2.real_value);
 
274
  case colvarvalue::type_3vector:
 
275
    return 2.0 * (this->rvector_value - x2.rvector_value);
 
276
  case colvarvalue::type_unit3vector:
 
277
  case colvarvalue::type_unit3vectorderiv:
 
278
    {
 
279
      cvm::rvector const &v1 = this->rvector_value;
 
280
      cvm::rvector const &v2 = x2.rvector_value;
 
281
      cvm::real const cos_t = v1 * v2;
 
282
      cvm::real const sin_t = std::sqrt(1.0 - cos_t*cos_t);
 
283
      return colvarvalue( 2.0 * sin_t *
 
284
                          cvm::rvector((-1.0) * sin_t * v2.x +
 
285
                                       cos_t/sin_t * (v1.x - cos_t*v2.x),
 
286
                                       (-1.0) * sin_t * v2.y +
 
287
                                       cos_t/sin_t * (v1.y - cos_t*v2.y),
 
288
                                       (-1.0) * sin_t * v2.z +
 
289
                                       cos_t/sin_t * (v1.z - cos_t*v2.z)
 
290
                                       ),
 
291
                          colvarvalue::type_unit3vector );
 
292
    }
 
293
  case colvarvalue::type_quaternion:
 
294
  case colvarvalue::type_quaternionderiv:
 
295
    return this->quaternion_value.dist2_grad(x2.quaternion_value);
 
296
  case colvarvalue::type_vector:
 
297
    return colvarvalue(2.0 * (this->vector1d_value - x2.vector1d_value), colvarvalue::type_vector);
 
298
    break;
 
299
  case colvarvalue::type_notset:
 
300
  default:
 
301
    this->undef_op();
 
302
    return colvarvalue(colvarvalue::type_notset);
 
303
  };
 
304
}
 
305
 
 
306
 
 
307
std::string colvarvalue::to_simple_string() const
 
308
{
 
309
  switch (type()) {
 
310
  case colvarvalue::type_scalar:
 
311
    return cvm::to_str(real_value, 0, cvm::cv_prec);
 
312
    break;
 
313
  case colvarvalue::type_3vector:
 
314
  case colvarvalue::type_unit3vector:
 
315
  case colvarvalue::type_unit3vectorderiv:
 
316
    return rvector_value.to_simple_string();
 
317
    break;
 
318
  case colvarvalue::type_quaternion:
 
319
  case colvarvalue::type_quaternionderiv:
 
320
    return quaternion_value.to_simple_string();
 
321
    break;
 
322
  case colvarvalue::type_vector:
 
323
    return vector1d_value.to_simple_string();
 
324
    break;
 
325
  case colvarvalue::type_notset:
 
326
  default:
 
327
    undef_op();
 
328
    break;
 
329
  }
 
330
  return std::string();
 
331
}
 
332
 
 
333
 
 
334
int colvarvalue::from_simple_string(std::string const &s)
 
335
{
 
336
  switch (type()) {
 
337
  case colvarvalue::type_scalar:
 
338
    return ((std::istringstream(s) >> real_value)
 
339
            ? COLVARS_OK : COLVARS_ERROR);
 
340
    break;
 
341
  case colvarvalue::type_3vector:
 
342
  case colvarvalue::type_unit3vector:
 
343
  case colvarvalue::type_unit3vectorderiv:
 
344
    return rvector_value.from_simple_string(s);
 
345
    break;
 
346
  case colvarvalue::type_quaternion:
 
347
  case colvarvalue::type_quaternionderiv:
 
348
    return quaternion_value.from_simple_string(s);
 
349
    break;
 
350
  case colvarvalue::type_vector:
 
351
    return vector1d_value.from_simple_string(s);
 
352
    break;
 
353
  case colvarvalue::type_notset:
 
354
  default:
 
355
    undef_op();
 
356
    break;
 
357
  }
 
358
  return COLVARS_ERROR;
 
359
}
194
360
 
195
361
std::ostream & operator << (std::ostream &os, colvarvalue const &x)
196
362
{
197
363
  switch (x.type()) {
198
364
  case colvarvalue::type_scalar:
199
 
    os << x.real_value; break;
 
365
    os << x.real_value;
 
366
    break;
 
367
  case colvarvalue::type_3vector:
 
368
  case colvarvalue::type_unit3vector:
 
369
  case colvarvalue::type_unit3vectorderiv:
 
370
    os << x.rvector_value;
 
371
    break;
 
372
  case colvarvalue::type_quaternion:
 
373
  case colvarvalue::type_quaternionderiv:
 
374
    os << x.quaternion_value;
 
375
    break;
200
376
  case colvarvalue::type_vector:
201
 
  case colvarvalue::type_unitvector:
202
 
    os << x.rvector_value; break;
203
 
  case colvarvalue::type_quaternion:
204
 
    os << x.quaternion_value; break;
 
377
    os << x.vector1d_value;
 
378
    break;
205
379
  case colvarvalue::type_notset:
206
 
    os << "not set"; break;
 
380
  default:
 
381
    os << "not set";
 
382
    break;
207
383
  }
208
384
  return os;
209
385
}
211
387
 
212
388
std::ostream & operator << (std::ostream &os, std::vector<colvarvalue> const &v)
213
389
{
214
 
  for (size_t i = 0; i < v.size(); i++) os << v[i];
 
390
  size_t i;
 
391
  for (i = 0; i < v.size(); i++) {
 
392
    os << v[i];
 
393
  }
215
394
  return os;
216
395
}
217
396
 
219
398
std::istream & operator >> (std::istream &is, colvarvalue &x)
220
399
{
221
400
  if (x.type() == colvarvalue::type_notset) {
222
 
    cvm::fatal_error ("Trying to read from a stream a colvarvalue, "
223
 
                      "which has not yet been assigned a data type.\n");
 
401
    cvm::error("Trying to read from a stream a colvarvalue, "
 
402
               "which has not yet been assigned a data type.\n");
 
403
    return is;
224
404
  }
225
405
 
226
406
  switch (x.type()) {
227
407
  case colvarvalue::type_scalar:
228
408
    is >> x.real_value;
229
409
    break;
 
410
  case colvarvalue::type_3vector:
 
411
  case colvarvalue::type_unit3vectorderiv:
 
412
    is >> x.rvector_value;
 
413
    break;
 
414
  case colvarvalue::type_unit3vector:
 
415
    is >> x.rvector_value;
 
416
    x.apply_constraints();
 
417
    break;
 
418
  case colvarvalue::type_quaternion:
 
419
    is >> x.quaternion_value;
 
420
    x.apply_constraints();
 
421
    break;
 
422
  case colvarvalue::type_quaternionderiv:
 
423
    is >> x.quaternion_value;
 
424
    break;
230
425
  case colvarvalue::type_vector:
231
 
  case colvarvalue::type_unitvector:
232
 
    is >> x.rvector_value;
233
 
    break;
234
 
  case colvarvalue::type_quaternion:
235
 
    is >> x.quaternion_value;
236
 
    break;
 
426
    is >> x.vector1d_value;
 
427
    break;
 
428
  case colvarvalue::type_notset:
237
429
  default:
238
430
    x.undef_op();
239
431
  }
240
 
  x.apply_constraints();
241
432
  return is;
242
433
}
243
434
 
244
435
 
245
 
size_t colvarvalue::output_width (size_t const &real_width) const
 
436
size_t colvarvalue::output_width(size_t const &real_width) const
246
437
{
247
438
  switch (this->value_type) {
248
439
  case colvarvalue::type_scalar:
249
440
    return real_width;
 
441
  case colvarvalue::type_3vector:
 
442
  case colvarvalue::type_unit3vector:
 
443
  case colvarvalue::type_unit3vectorderiv:
 
444
    return cvm::rvector::output_width(real_width);
 
445
  case colvarvalue::type_quaternion:
 
446
  case colvarvalue::type_quaternionderiv:
 
447
    return cvm::quaternion::output_width(real_width);
250
448
  case colvarvalue::type_vector:
251
 
  case colvarvalue::type_unitvector:
252
 
    return cvm::rvector::output_width (real_width);
253
 
  case colvarvalue::type_quaternion:
254
 
    return cvm::quaternion::output_width (real_width);
 
449
    // note how this depends on the vector's size
 
450
    return vector1d_value.output_width(real_width);
255
451
  case colvarvalue::type_notset:
256
452
  default:
257
453
    return 0;
259
455
}
260
456
 
261
457
 
 
458
void colvarvalue::inner_opt(colvarvalue                        const &x,
 
459
                            std::vector<colvarvalue>::iterator       &xv,
 
460
                            std::vector<colvarvalue>::iterator const &xv_end,
 
461
                            std::vector<cvm::real>::iterator         &result)
 
462
{
 
463
  // doing type check only once, here
 
464
  colvarvalue::check_types(x, *xv);
 
465
 
 
466
  std::vector<colvarvalue>::iterator &xvi = xv;
 
467
  std::vector<cvm::real>::iterator    &ii = result;
 
468
 
 
469
  switch (x.value_type) {
 
470
  case colvarvalue::type_scalar:
 
471
    while (xvi != xv_end) {
 
472
      *(ii++) += (xvi++)->real_value * x.real_value;
 
473
    }
 
474
    break;
 
475
  case colvarvalue::type_3vector:
 
476
  case colvarvalue::type_unit3vector:
 
477
  case colvarvalue::type_unit3vectorderiv:
 
478
    while (xvi != xv_end) {
 
479
      *(ii++) += (xvi++)->rvector_value * x.rvector_value;
 
480
    }
 
481
    break;
 
482
  case colvarvalue::type_quaternion:
 
483
  case colvarvalue::type_quaternionderiv:
 
484
    while (xvi != xv_end) {
 
485
      *(ii++) += ((xvi++)->quaternion_value).cosine(x.quaternion_value);
 
486
    }
 
487
    break;
 
488
  case colvarvalue::type_vector:
 
489
    while (xvi != xv_end) {
 
490
      *(ii++) += (xvi++)->vector1d_value * x.vector1d_value;
 
491
    }
 
492
    break;
 
493
  default:
 
494
    x.undef_op();
 
495
  };
 
496
}
 
497
 
 
498
 
 
499
void colvarvalue::inner_opt(colvarvalue const                      &x,
 
500
                            std::list<colvarvalue>::iterator       &xv,
 
501
                            std::list<colvarvalue>::iterator const &xv_end,
 
502
                            std::vector<cvm::real>::iterator       &result)
 
503
{
 
504
  // doing type check only once, here
 
505
  colvarvalue::check_types(x, *xv);
 
506
 
 
507
  std::list<colvarvalue>::iterator &xvi = xv;
 
508
  std::vector<cvm::real>::iterator  &ii = result;
 
509
 
 
510
  switch (x.value_type) {
 
511
  case colvarvalue::type_scalar:
 
512
    while (xvi != xv_end) {
 
513
      *(ii++) += (xvi++)->real_value * x.real_value;
 
514
    }
 
515
    break;
 
516
  case colvarvalue::type_3vector:
 
517
  case colvarvalue::type_unit3vector:
 
518
  case colvarvalue::type_unit3vectorderiv:
 
519
    while (xvi != xv_end) {
 
520
      *(ii++) += (xvi++)->rvector_value * x.rvector_value;
 
521
    }
 
522
    break;
 
523
  case colvarvalue::type_quaternion:
 
524
  case colvarvalue::type_quaternionderiv:
 
525
    while (xvi != xv_end) {
 
526
      *(ii++) += ((xvi++)->quaternion_value).cosine(x.quaternion_value);
 
527
    }
 
528
    break;
 
529
  case colvarvalue::type_vector:
 
530
    while (xvi != xv_end) {
 
531
      *(ii++) += (xvi++)->vector1d_value * x.vector1d_value;
 
532
    }
 
533
    break;
 
534
  default:
 
535
    x.undef_op();
 
536
  };
 
537
}
 
538
 
 
539
 
 
540
void colvarvalue::p2leg_opt(colvarvalue const                        &x,
 
541
                            std::vector<colvarvalue>::iterator       &xv,
 
542
                            std::vector<colvarvalue>::iterator const &xv_end,
 
543
                            std::vector<cvm::real>::iterator         &result)
 
544
{
 
545
  // doing type check only once, here
 
546
  colvarvalue::check_types(x, *xv);
 
547
 
 
548
  std::vector<colvarvalue>::iterator &xvi = xv;
 
549
  std::vector<cvm::real>::iterator    &ii = result;
 
550
 
 
551
  switch (x.value_type) {
 
552
  case colvarvalue::type_scalar:
 
553
    cvm::error("Error: cannot calculate Legendre polynomials "
 
554
               "for scalar variables.\n");
 
555
    return;
 
556
    break;
 
557
  case colvarvalue::type_3vector:
 
558
    while (xvi != xv_end) {
 
559
      cvm::real const cosine =
 
560
        ((xvi)->rvector_value * x.rvector_value) /
 
561
        ((xvi)->rvector_value.norm() * x.rvector_value.norm());
 
562
      xvi++;
 
563
      *(ii++) += 1.5*cosine*cosine - 0.5;
 
564
    }
 
565
    break;
 
566
  case colvarvalue::type_unit3vector:
 
567
  case colvarvalue::type_unit3vectorderiv:
 
568
    while (xvi != xv_end) {
 
569
      cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value;
 
570
      *(ii++) += 1.5*cosine*cosine - 0.5;
 
571
    }
 
572
    break;
 
573
  case colvarvalue::type_quaternion:
 
574
  case colvarvalue::type_quaternionderiv:
 
575
    while (xvi != xv_end) {
 
576
      cvm::real const cosine = (xvi++)->quaternion_value.cosine(x.quaternion_value);
 
577
      *(ii++) += 1.5*cosine*cosine - 0.5;
 
578
    }
 
579
    break;
 
580
  case colvarvalue::type_vector:
 
581
    while (xvi != xv_end) {
 
582
      cvm::real const cosine =
 
583
        ((xvi)->vector1d_value * x.vector1d_value) /
 
584
        ((xvi)->vector1d_value.norm() * x.rvector_value.norm());
 
585
      xvi++;
 
586
      *(ii++) += 1.5*cosine*cosine - 0.5;
 
587
    }
 
588
    break;
 
589
  default:
 
590
    x.undef_op();
 
591
  };
 
592
}
 
593
 
 
594
 
 
595
void colvarvalue::p2leg_opt(colvarvalue const                        &x,
 
596
                            std::list<colvarvalue>::iterator         &xv,
 
597
                            std::list<colvarvalue>::iterator const   &xv_end,
 
598
                            std::vector<cvm::real>::iterator         &result)
 
599
{
 
600
  // doing type check only once, here
 
601
  colvarvalue::check_types(x, *xv);
 
602
 
 
603
  std::list<colvarvalue>::iterator &xvi = xv;
 
604
  std::vector<cvm::real>::iterator  &ii = result;
 
605
 
 
606
  switch (x.value_type) {
 
607
  case colvarvalue::type_scalar:
 
608
    cvm::error("Error: cannot calculate Legendre polynomials "
 
609
               "for scalar variables.\n");
 
610
    break;
 
611
  case colvarvalue::type_3vector:
 
612
    while (xvi != xv_end) {
 
613
      cvm::real const cosine =
 
614
        ((xvi)->rvector_value * x.rvector_value) /
 
615
        ((xvi)->rvector_value.norm() * x.rvector_value.norm());
 
616
      xvi++;
 
617
      *(ii++) += 1.5*cosine*cosine - 0.5;
 
618
    }
 
619
    break;
 
620
  case colvarvalue::type_unit3vector:
 
621
  case colvarvalue::type_unit3vectorderiv:
 
622
    while (xvi != xv_end) {
 
623
      cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value;
 
624
      *(ii++) += 1.5*cosine*cosine - 0.5;
 
625
    }
 
626
    break;
 
627
  case colvarvalue::type_quaternion:
 
628
  case colvarvalue::type_quaternionderiv:
 
629
    while (xvi != xv_end) {
 
630
      cvm::real const cosine = (xvi++)->quaternion_value.cosine(x.quaternion_value);
 
631
      *(ii++) += 1.5*cosine*cosine - 0.5;
 
632
    }
 
633
    break;
 
634
  default:
 
635
    x.undef_op();
 
636
  };
 
637
}
 
638
 
 
639
 
 
640