Update Colvars to version 2018-12-18
[namd.git] / colvars / src / colvarcomp_rotations.cpp
1 // -*- c++ -*-
2
3 // This file is part of the Collective Variables module (Colvars).
4 // The original version of Colvars and its updates are located at:
5 // https://github.com/colvars/colvars
6 // Please update all Colvars source files before making any changes.
7 // If you wish to distribute your changes, please submit them to the
8 // Colvars repository at GitHub.
9
10 #include <cmath>
11
12 #include "colvarmodule.h"
13 #include "colvarvalue.h"
14 #include "colvarparse.h"
15 #include "colvar.h"
16 #include "colvarcomp.h"
17
18
19
20 colvar::orientation::orientation(std::string const &conf)
21   : cvc(conf)
22 {
23   function_type = "orientation";
24   enable(f_cvc_implicit_gradient);
25   x.type(colvarvalue::type_quaternion);
26   init(conf);
27 }
28
29
30 int colvar::orientation::init(std::string const &conf)
31 {
32   atoms = parse_group(conf, "atoms");
33   ref_pos.reserve(atoms->size());
34
35   if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {
36     cvm::log("Using reference positions from input file.\n");
37     if (ref_pos.size() != atoms->size()) {
38       return cvm::error("Error: reference positions do not "
39                         "match the number of requested atoms.\n", INPUT_ERROR);
40     }
41   }
42
43   {
44     std::string file_name;
45     if (get_keyval(conf, "refPositionsFile", file_name)) {
46
47       std::string file_col;
48       double file_col_value=0.0;
49       if (get_keyval(conf, "refPositionsCol", file_col, std::string(""))) {
50         // use PDB flags if column is provided
51         bool found = get_keyval(conf, "refPositionsColValue", file_col_value, 0.0);
52         if (found && file_col_value==0.0) {
53           return cvm::error("Error: refPositionsColValue, "
54                             "if provided, must be non-zero.\n", INPUT_ERROR);
55         }
56       }
57
58       ref_pos.resize(atoms->size());
59       cvm::load_coords(file_name.c_str(), &ref_pos, atoms,
60                        file_col, file_col_value);
61     }
62   }
63
64   if (!ref_pos.size()) {
65     return cvm::error("Error: must define a set of "
66                       "reference coordinates.\n", INPUT_ERROR);
67   }
68
69
70   cvm::log("Centering the reference coordinates: it is "
71             "assumed that each atom is the closest "
72             "periodic image to the center of geometry.\n");
73   cvm::rvector ref_cog(0.0, 0.0, 0.0);
74   size_t i;
75   for (i = 0; i < ref_pos.size(); i++) {
76     ref_cog += ref_pos[i];
77   }
78   ref_cog /= cvm::real(ref_pos.size());
79   for (i = 0; i < ref_pos.size(); i++) {
80     ref_pos[i] -= ref_cog;
81   }
82
83   get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0));
84
85   // initialize rot member data
86   if (!atoms->noforce) {
87     rot.request_group2_gradients(atoms->size());
88   }
89
90 }
91
92
93 colvar::orientation::orientation()
94   : cvc()
95 {
96   function_type = "orientation";
97   enable(f_cvc_implicit_gradient);
98   x.type(colvarvalue::type_quaternion);
99 }
100
101
102 void colvar::orientation::calc_value()
103 {
104   rot.b_debug_gradients = is_enabled(f_cvc_debug_gradient);
105   atoms_cog = atoms->center_of_geometry();
106
107   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
108
109   if ((rot.q).inner(ref_quat) >= 0.0) {
110     x.quaternion_value = rot.q;
111   } else {
112     x.quaternion_value = -1.0 * rot.q;
113   }
114 }
115
116
117 void colvar::orientation::calc_gradients()
118 {
119   // gradients have already been calculated and stored within the
120   // member object "rot"; we're not using the "grad" member of each
121   // atom object, because it only can represent the gradient of a
122   // scalar colvar
123 }
124
125
126 void colvar::orientation::apply_force(colvarvalue const &force)
127 {
128   cvm::quaternion const &FQ = force.quaternion_value;
129
130   if (!atoms->noforce) {
131     for (size_t ia = 0; ia < atoms->size(); ia++) {
132       for (size_t i = 0; i < 4; i++) {
133         (*atoms)[ia].apply_force(FQ[i] * rot.dQ0_2[ia][i]);
134       }
135     }
136   }
137 }
138
139
140 cvm::real colvar::orientation::dist2(colvarvalue const &x1,
141                                      colvarvalue const &x2) const
142 {
143   return x1.quaternion_value.dist2(x2);
144 }
145
146
147 colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
148                                              colvarvalue const &x2) const
149 {
150   return x1.quaternion_value.dist2_grad(x2);
151 }
152
153
154 colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
155                                              colvarvalue const &x2) const
156 {
157   return x2.quaternion_value.dist2_grad(x1);
158 }
159
160
161
162 colvar::orientation_angle::orientation_angle(std::string const &conf)
163   : orientation()
164 {
165   function_type = "orientation_angle";
166   disable(f_cvc_implicit_gradient);
167   x.type(colvarvalue::type_scalar);
168   init(conf);
169 }
170
171
172 int colvar::orientation_angle::init(std::string const &conf)
173 {
174   return orientation::init(conf);
175 }
176
177
178 colvar::orientation_angle::orientation_angle()
179   : orientation()
180 {
181   function_type = "orientation_angle";
182   disable(f_cvc_implicit_gradient);
183   x.type(colvarvalue::type_scalar);
184 }
185
186
187 void colvar::orientation_angle::calc_value()
188 {
189   atoms_cog = atoms->center_of_geometry();
190
191   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
192
193   if ((rot.q).q0 >= 0.0) {
194     x.real_value = (180.0/PI) * 2.0 * std::acos((rot.q).q0);
195   } else {
196     x.real_value = (180.0/PI) * 2.0 * std::acos(-1.0 * (rot.q).q0);
197   }
198 }
199
200
201 void colvar::orientation_angle::calc_gradients()
202 {
203   cvm::real const dxdq0 =
204     ( ((rot.q).q0 * (rot.q).q0 < 1.0) ?
205       ((180.0 / PI) * (-2.0) / std::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) :
206       0.0 );
207
208   for (size_t ia = 0; ia < atoms->size(); ia++) {
209     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
210   }
211 }
212
213
214 void colvar::orientation_angle::apply_force(colvarvalue const &force)
215 {
216   cvm::real const &fw = force.real_value;
217   if (!atoms->noforce) {
218     atoms->apply_colvar_force(fw);
219   }
220 }
221
222
223 simple_scalar_dist_functions(orientation_angle)
224
225
226
227 colvar::orientation_proj::orientation_proj(std::string const &conf)
228   : orientation()
229 {
230   function_type = "orientation_proj";
231   disable(f_cvc_implicit_gradient);
232   x.type(colvarvalue::type_scalar);
233   init(conf);
234 }
235
236
237 int colvar::orientation_proj::init(std::string const &conf)
238 {
239   return orientation::init(conf);
240 }
241
242
243 colvar::orientation_proj::orientation_proj()
244   : orientation()
245 {
246   function_type = "orientation_proj";
247   disable(f_cvc_implicit_gradient);
248   x.type(colvarvalue::type_scalar);
249 }
250
251
252 void colvar::orientation_proj::calc_value()
253 {
254   atoms_cog = atoms->center_of_geometry();
255   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
256   x.real_value = 2.0 * (rot.q).q0 * (rot.q).q0 - 1.0;
257 }
258
259
260 void colvar::orientation_proj::calc_gradients()
261 {
262   cvm::real const dxdq0 = 2.0 * 2.0 * (rot.q).q0;
263   for (size_t ia = 0; ia < atoms->size(); ia++) {
264     (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
265   }
266 }
267
268
269 void colvar::orientation_proj::apply_force(colvarvalue const &force)
270 {
271   cvm::real const &fw = force.real_value;
272
273   if (!atoms->noforce) {
274     atoms->apply_colvar_force(fw);
275   }
276 }
277
278
279 simple_scalar_dist_functions(orientation_proj)
280
281
282
283 colvar::tilt::tilt(std::string const &conf)
284   : orientation()
285 {
286   function_type = "tilt";
287   disable(f_cvc_implicit_gradient);
288   x.type(colvarvalue::type_scalar);
289   init(conf);
290 }
291
292
293 int colvar::tilt::init(std::string const &conf)
294 {
295   int error_code = COLVARS_OK;
296
297   error_code |= orientation::init(conf);
298
299   get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0));
300   if (axis.norm2() != 1.0) {
301     axis /= axis.norm();
302     cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n");
303   }
304
305   return error_code;
306 }
307
308
309 colvar::tilt::tilt()
310   : orientation()
311 {
312   function_type = "tilt";
313   disable(f_cvc_implicit_gradient);
314   x.type(colvarvalue::type_scalar);
315 }
316
317
318 void colvar::tilt::calc_value()
319 {
320   atoms_cog = atoms->center_of_geometry();
321
322   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
323
324   x.real_value = rot.cos_theta(axis);
325 }
326
327
328 void colvar::tilt::calc_gradients()
329 {
330   cvm::quaternion const dxdq = rot.dcos_theta_dq(axis);
331
332   for (size_t ia = 0; ia < atoms->size(); ia++) {
333     (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
334     for (size_t iq = 0; iq < 4; iq++) {
335       (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
336     }
337   }
338 }
339
340
341 void colvar::tilt::apply_force(colvarvalue const &force)
342 {
343   cvm::real const &fw = force.real_value;
344
345   if (!atoms->noforce) {
346     atoms->apply_colvar_force(fw);
347   }
348 }
349
350
351 simple_scalar_dist_functions(tilt)
352
353
354
355 colvar::spin_angle::spin_angle(std::string const &conf)
356   : orientation()
357 {
358   function_type = "spin_angle";
359   period = 360.0;
360   b_periodic = true;
361   disable(f_cvc_implicit_gradient);
362   x.type(colvarvalue::type_scalar);
363   init(conf);
364 }
365
366
367 int colvar::spin_angle::init(std::string const &conf)
368 {
369   int error_code = COLVARS_OK;
370
371   error_code |= orientation::init(conf);
372
373   get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0));
374   if (axis.norm2() != 1.0) {
375     axis /= axis.norm();
376     cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n");
377   }
378
379   return error_code;
380 }
381
382
383 colvar::spin_angle::spin_angle()
384   : orientation()
385 {
386   function_type = "spin_angle";
387   period = 360.0;
388   b_periodic = true;
389   disable(f_cvc_implicit_gradient);
390   x.type(colvarvalue::type_scalar);
391 }
392
393
394 void colvar::spin_angle::calc_value()
395 {
396   atoms_cog = atoms->center_of_geometry();
397
398   rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
399
400   x.real_value = rot.spin_angle(axis);
401   this->wrap(x);
402 }
403
404
405 void colvar::spin_angle::calc_gradients()
406 {
407   cvm::quaternion const dxdq = rot.dspin_angle_dq(axis);
408
409   for (size_t ia = 0; ia < atoms->size(); ia++) {
410     (*atoms)[ia].grad = cvm::rvector(0.0, 0.0, 0.0);
411     for (size_t iq = 0; iq < 4; iq++) {
412       (*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
413     }
414   }
415 }
416
417
418 void colvar::spin_angle::apply_force(colvarvalue const &force)
419 {
420   cvm::real const &fw = force.real_value;
421
422   if (!atoms->noforce) {
423     atoms->apply_colvar_force(fw);
424   }
425 }
426
427
428 cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
429                                     colvarvalue const &x2) const
430 {
431   cvm::real diff = x1.real_value - x2.real_value;
432   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
433   return diff * diff;
434 }
435
436
437 colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
438                                             colvarvalue const &x2) const
439 {
440   cvm::real diff = x1.real_value - x2.real_value;
441   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
442   return 2.0 * diff;
443 }
444
445
446 colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
447                                             colvarvalue const &x2) const
448 {
449   cvm::real diff = x1.real_value - x2.real_value;
450   diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
451   return (-2.0) * diff;
452 }
453
454
455 void colvar::spin_angle::wrap(colvarvalue &x) const
456 {
457   if ((x.real_value - wrap_center) >= 180.0) {
458     x.real_value -= 360.0;
459     return;
460   }
461
462   if ((x.real_value - wrap_center) < -180.0) {
463     x.real_value += 360.0;
464     return;
465   }
466
467   return;
468 }