Update Colvars to version 2018-12-18
[namd.git] / colvars / src / colvarcomp_distances.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::distance::distance(std::string const &conf)
21   : cvc(conf)
22 {
23   function_type = "distance";
24   provide(f_cvc_inv_gradient);
25   provide(f_cvc_Jacobian);
26   enable(f_cvc_com_based);
27
28   group1 = parse_group(conf, "group1");
29   group2 = parse_group(conf, "group2");
30
31   init_total_force_params(conf);
32
33   x.type(colvarvalue::type_scalar);
34 }
35
36
37 colvar::distance::distance()
38   : cvc()
39 {
40   function_type = "distance";
41   provide(f_cvc_inv_gradient);
42   provide(f_cvc_Jacobian);
43   enable(f_cvc_com_based);
44   x.type(colvarvalue::type_scalar);
45 }
46
47
48 void colvar::distance::calc_value()
49 {
50   if (!is_enabled(f_cvc_pbc_minimum_image)) {
51     dist_v = group2->center_of_mass() - group1->center_of_mass();
52   } else {
53     dist_v = cvm::position_distance(group1->center_of_mass(),
54                                     group2->center_of_mass());
55   }
56   x.real_value = dist_v.norm();
57 }
58
59
60 void colvar::distance::calc_gradients()
61 {
62   cvm::rvector const u = dist_v.unit();
63   group1->set_weighted_gradient(-1.0 * u);
64   group2->set_weighted_gradient(       u);
65 }
66
67
68 void colvar::distance::calc_force_invgrads()
69 {
70   group1->read_total_forces();
71   if (is_enabled(f_cvc_one_site_total_force)) {
72     ft.real_value = -1.0 * (group1->total_force() * dist_v.unit());
73   } else {
74     group2->read_total_forces();
75     ft.real_value = 0.5 * ((group2->total_force() - group1->total_force()) * dist_v.unit());
76   }
77 }
78
79
80 void colvar::distance::calc_Jacobian_derivative()
81 {
82   jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0;
83 }
84
85
86 void colvar::distance::apply_force(colvarvalue const &force)
87 {
88   if (!group1->noforce)
89     group1->apply_colvar_force(force.real_value);
90
91   if (!group2->noforce)
92     group2->apply_colvar_force(force.real_value);
93 }
94
95
96 simple_scalar_dist_functions(distance)
97
98
99
100 colvar::distance_vec::distance_vec(std::string const &conf)
101   : distance(conf)
102 {
103   function_type = "distance_vec";
104   enable(f_cvc_com_based);
105   enable(f_cvc_implicit_gradient);
106   x.type(colvarvalue::type_3vector);
107 }
108
109
110 colvar::distance_vec::distance_vec()
111   : distance()
112 {
113   function_type = "distance_vec";
114   enable(f_cvc_com_based);
115   enable(f_cvc_implicit_gradient);
116   x.type(colvarvalue::type_3vector);
117 }
118
119
120 void colvar::distance_vec::calc_value()
121 {
122   if (!is_enabled(f_cvc_pbc_minimum_image)) {
123     x.rvector_value = group2->center_of_mass() - group1->center_of_mass();
124   } else {
125     x.rvector_value = cvm::position_distance(group1->center_of_mass(),
126                                              group2->center_of_mass());
127   }
128 }
129
130
131 void colvar::distance_vec::calc_gradients()
132 {
133   // gradients are not stored: a 3x3 matrix for each atom would be
134   // needed to store just the identity matrix
135 }
136
137
138 void colvar::distance_vec::apply_force(colvarvalue const &force)
139 {
140   if (!group1->noforce)
141     group1->apply_force(-1.0 * force.rvector_value);
142
143   if (!group2->noforce)
144     group2->apply_force(       force.rvector_value);
145 }
146
147
148 cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
149                                       colvarvalue const &x2) const
150 {
151   return (cvm::position_distance(x1.rvector_value, x2.rvector_value)).norm2();
152 }
153
154
155 colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
156                                               colvarvalue const &x2) const
157 {
158   return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
159 }
160
161
162 colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
163                                               colvarvalue const &x2) const
164 {
165   return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
166 }
167
168
169
170 colvar::distance_z::distance_z(std::string const &conf)
171   : cvc(conf)
172 {
173   function_type = "distance_z";
174   provide(f_cvc_inv_gradient);
175   provide(f_cvc_Jacobian);
176   enable(f_cvc_com_based);
177   x.type(colvarvalue::type_scalar);
178
179   // TODO detect PBC from MD engine (in simple cases)
180   // and then update period in real time
181   if (period != 0.0)
182     b_periodic = true;
183
184   if ((wrap_center != 0.0) && (period == 0.0)) {
185     cvm::error("Error: wrapAround was defined in a distanceZ component,"
186                 " but its period has not been set.\n");
187     return;
188   }
189
190   main = parse_group(conf, "main");
191   ref1 = parse_group(conf, "ref");
192   // this group is optional
193   ref2 = parse_group(conf, "ref2", true);
194
195   if ( ref2 ) {
196     cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
197     fixed_axis = false;
198     if (key_lookup(conf, "axis"))
199       cvm::log("Warning: explicit axis definition will be ignored!");
200   } else {
201     if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) {
202       if (axis.norm2() == 0.0) {
203         cvm::error("Axis vector is zero!");
204         return;
205       }
206       if (axis.norm2() != 1.0) {
207         axis = axis.unit();
208         cvm::log("The normalized axis is: "+cvm::to_str(axis)+".\n");
209       }
210     }
211     fixed_axis = true;
212   }
213
214   init_total_force_params(conf);
215
216 }
217
218
219 colvar::distance_z::distance_z()
220 {
221   function_type = "distance_z";
222   provide(f_cvc_inv_gradient);
223   provide(f_cvc_Jacobian);
224   enable(f_cvc_com_based);
225   x.type(colvarvalue::type_scalar);
226 }
227
228
229 void colvar::distance_z::calc_value()
230 {
231   if (fixed_axis) {
232     if (!is_enabled(f_cvc_pbc_minimum_image)) {
233       dist_v = main->center_of_mass() - ref1->center_of_mass();
234     } else {
235       dist_v = cvm::position_distance(ref1->center_of_mass(),
236                                       main->center_of_mass());
237     }
238   } else {
239
240     if (!is_enabled(f_cvc_pbc_minimum_image)) {
241       dist_v = main->center_of_mass() -
242                (0.5 * (ref1->center_of_mass() + ref2->center_of_mass()));
243       axis = ref2->center_of_mass() - ref1->center_of_mass();
244     } else {
245       dist_v = cvm::position_distance(0.5 * (ref1->center_of_mass() +
246                                              ref2->center_of_mass()),
247                                       main->center_of_mass());
248       axis = cvm::position_distance(ref1->center_of_mass(),
249                                     ref2->center_of_mass());
250     }
251     axis_norm = axis.norm();
252     axis = axis.unit();
253   }
254   x.real_value = axis * dist_v;
255   this->wrap(x);
256 }
257
258
259 void colvar::distance_z::calc_gradients()
260 {
261   main->set_weighted_gradient( axis );
262
263   if (fixed_axis) {
264     ref1->set_weighted_gradient(-1.0 * axis);
265   } else {
266     if (!is_enabled(f_cvc_pbc_minimum_image)) {
267       ref1->set_weighted_gradient( 1.0 / axis_norm *
268                                    (main->center_of_mass() - ref2->center_of_mass() -
269                                    x.real_value * axis ));
270       ref2->set_weighted_gradient( 1.0 / axis_norm *
271                                    (ref1->center_of_mass() - main->center_of_mass() +
272                                    x.real_value * axis ));
273     } else {
274       ref1->set_weighted_gradient( 1.0 / axis_norm * (
275         cvm::position_distance(ref2->center_of_mass(),
276                                main->center_of_mass()) - x.real_value * axis ));
277       ref2->set_weighted_gradient( 1.0 / axis_norm * (
278         cvm::position_distance(main->center_of_mass(),
279                                ref1->center_of_mass()) + x.real_value * axis ));
280     }
281   }
282 }
283
284
285 void colvar::distance_z::calc_force_invgrads()
286 {
287   main->read_total_forces();
288
289   if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
290     ref1->read_total_forces();
291     ft.real_value = 0.5 * ((main->total_force() - ref1->total_force()) * axis);
292   } else {
293     ft.real_value = main->total_force() * axis;
294   }
295 }
296
297
298 void colvar::distance_z::calc_Jacobian_derivative()
299 {
300   jd.real_value = 0.0;
301 }
302
303
304 void colvar::distance_z::apply_force(colvarvalue const &force)
305 {
306   if (!ref1->noforce)
307     ref1->apply_colvar_force(force.real_value);
308
309   if (ref2 && !ref2->noforce)
310     ref2->apply_colvar_force(force.real_value);
311
312   if (!main->noforce)
313     main->apply_colvar_force(force.real_value);
314 }
315
316
317 // Differences should always be wrapped around 0 (ignoring wrap_center)
318 cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
319                                     colvarvalue const &x2) const
320 {
321   cvm::real diff = x1.real_value - x2.real_value;
322   if (b_periodic) {
323     cvm::real shift = std::floor(diff/period + 0.5);
324     diff -= shift * period;
325   }
326   return diff * diff;
327 }
328
329
330 colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
331                                             colvarvalue const &x2) const
332 {
333   cvm::real diff = x1.real_value - x2.real_value;
334   if (b_periodic) {
335     cvm::real shift = std::floor(diff/period + 0.5);
336     diff -= shift * period;
337   }
338   return 2.0 * diff;
339 }
340
341
342 colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
343                                             colvarvalue const &x2) const
344 {
345   cvm::real diff = x1.real_value - x2.real_value;
346   if (b_periodic) {
347     cvm::real shift = std::floor(diff/period + 0.5);
348     diff -= shift * period;
349   }
350   return (-2.0) * diff;
351 }
352
353
354 void colvar::distance_z::wrap(colvarvalue &x) const
355 {
356   if (!b_periodic) {
357     // don't wrap if the period has not been set
358     return;
359   }
360
361   cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
362   x.real_value -= shift * period;
363   return;
364 }
365
366
367
368 colvar::distance_xy::distance_xy(std::string const &conf)
369   : distance_z(conf)
370 {
371   function_type = "distance_xy";
372   provide(f_cvc_inv_gradient);
373   provide(f_cvc_Jacobian);
374   enable(f_cvc_com_based);
375   x.type(colvarvalue::type_scalar);
376 }
377
378
379 colvar::distance_xy::distance_xy()
380   : distance_z()
381 {
382   function_type = "distance_xy";
383   provide(f_cvc_inv_gradient);
384   provide(f_cvc_Jacobian);
385   enable(f_cvc_com_based);
386   x.type(colvarvalue::type_scalar);
387 }
388
389
390 void colvar::distance_xy::calc_value()
391 {
392   if (!is_enabled(f_cvc_pbc_minimum_image)) {
393     dist_v = main->center_of_mass() - ref1->center_of_mass();
394   } else {
395     dist_v = cvm::position_distance(ref1->center_of_mass(),
396                                     main->center_of_mass());
397   }
398   if (!fixed_axis) {
399     if (!is_enabled(f_cvc_pbc_minimum_image)) {
400       v12 = ref2->center_of_mass() - ref1->center_of_mass();
401     } else {
402       v12 = cvm::position_distance(ref1->center_of_mass(),
403                                    ref2->center_of_mass());
404     }
405     axis_norm = v12.norm();
406     axis = v12.unit();
407   }
408
409   dist_v_ortho = dist_v - (dist_v * axis) * axis;
410   x.real_value = dist_v_ortho.norm();
411 }
412
413
414 void colvar::distance_xy::calc_gradients()
415 {
416   // Intermediate quantity (r_P3 / r_12 where P is the projection
417   // of 3(main) on the plane orthogonal to 12, containing 1 (ref1))
418   cvm::real A;
419   cvm::real x_inv;
420
421   if (x.real_value == 0.0) return;
422   x_inv = 1.0 / x.real_value;
423
424   if (fixed_axis) {
425     ref1->set_weighted_gradient(-1.0 * x_inv * dist_v_ortho);
426     main->set_weighted_gradient(       x_inv * dist_v_ortho);
427   } else {
428     if (!is_enabled(f_cvc_pbc_minimum_image)) {
429       v13 = main->center_of_mass() - ref1->center_of_mass();
430     } else {
431       v13 = cvm::position_distance(ref1->center_of_mass(),
432                                    main->center_of_mass());
433     }
434     A = (dist_v * axis) / axis_norm;
435
436     ref1->set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho);
437     ref2->set_weighted_gradient( -A        * x_inv * dist_v_ortho);
438     main->set_weighted_gradient(      1.0  * x_inv * dist_v_ortho);
439   }
440 }
441
442
443 void colvar::distance_xy::calc_force_invgrads()
444 {
445   main->read_total_forces();
446
447   if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
448     ref1->read_total_forces();
449     ft.real_value = 0.5 / x.real_value * ((main->total_force() - ref1->total_force()) * dist_v_ortho);
450   } else {
451     ft.real_value = 1.0 / x.real_value * main->total_force() * dist_v_ortho;
452   }
453 }
454
455
456 void colvar::distance_xy::calc_Jacobian_derivative()
457 {
458   jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
459 }
460
461
462 void colvar::distance_xy::apply_force(colvarvalue const &force)
463 {
464   if (!ref1->noforce)
465     ref1->apply_colvar_force(force.real_value);
466
467   if (ref2 && !ref2->noforce)
468     ref2->apply_colvar_force(force.real_value);
469
470   if (!main->noforce)
471     main->apply_colvar_force(force.real_value);
472 }
473
474
475 simple_scalar_dist_functions(distance_xy)
476
477
478
479 colvar::distance_dir::distance_dir(std::string const &conf)
480   : distance(conf)
481 {
482   function_type = "distance_dir";
483   enable(f_cvc_com_based);
484   enable(f_cvc_implicit_gradient);
485   x.type(colvarvalue::type_unit3vector);
486 }
487
488
489 colvar::distance_dir::distance_dir()
490   : distance()
491 {
492   function_type = "distance_dir";
493   enable(f_cvc_com_based);
494   enable(f_cvc_implicit_gradient);
495   x.type(colvarvalue::type_unit3vector);
496 }
497
498
499 void colvar::distance_dir::calc_value()
500 {
501   if (!is_enabled(f_cvc_pbc_minimum_image)) {
502     dist_v = group2->center_of_mass() - group1->center_of_mass();
503   } else {
504     dist_v = cvm::position_distance(group1->center_of_mass(),
505                                     group2->center_of_mass());
506   }
507   x.rvector_value = dist_v.unit();
508 }
509
510
511 void colvar::distance_dir::calc_gradients()
512 {
513   // gradients are computed on the fly within apply_force()
514   // Note: could be a problem if a future bias relies on gradient
515   // calculations...
516   // TODO in new deps system: remove dependency of biasing force to gradient?
517   // That way we could tell apart an explicit gradient dependency
518 }
519
520
521 void colvar::distance_dir::apply_force(colvarvalue const &force)
522 {
523   // remove the radial force component
524   cvm::real const iprod = force.rvector_value * x.rvector_value;
525   cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value;
526
527   if (!group1->noforce)
528     group1->apply_force(-1.0 * force_tang);
529
530   if (!group2->noforce)
531     group2->apply_force(       force_tang);
532 }
533
534
535 cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
536                                       colvarvalue const &x2) const
537 {
538   return (x1.rvector_value - x2.rvector_value).norm2();
539 }
540
541
542 colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
543                                               colvarvalue const &x2) const
544 {
545   return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vectorderiv);
546 }
547
548
549 colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
550                                               colvarvalue const &x2) const
551 {
552   return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vectorderiv);
553 }
554
555
556
557 colvar::distance_inv::distance_inv(std::string const &conf)
558   : cvc(conf)
559 {
560   function_type = "distance_inv";
561
562   group1 = parse_group(conf, "group1");
563   group2 = parse_group(conf, "group2");
564
565   get_keyval(conf, "exponent", exponent, 6);
566   if (exponent%2) {
567     cvm::error("Error: odd exponent provided, can only use even ones.\n");
568     return;
569   }
570   if (exponent <= 0) {
571     cvm::error("Error: negative or zero exponent provided.\n");
572     return;
573   }
574
575   for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
576     for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
577       if (ai1->id == ai2->id) {
578         cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n");
579         return;
580       }
581     }
582   }
583
584   if (is_enabled(f_cvc_debug_gradient)) {
585     cvm::log("Warning: debugGradients will not give correct results "
586              "for distanceInv, because its value and gradients are computed "
587              "simultaneously.\n");
588   }
589
590   x.type(colvarvalue::type_scalar);
591 }
592
593
594 colvar::distance_inv::distance_inv()
595 {
596   function_type = "distance_inv";
597   exponent = 6;
598   x.type(colvarvalue::type_scalar);
599 }
600
601
602 void colvar::distance_inv::calc_value()
603 {
604   x.real_value = 0.0;
605   if (!is_enabled(f_cvc_pbc_minimum_image)) {
606     for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
607       for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
608         cvm::rvector const dv = ai2->pos - ai1->pos;
609         cvm::real const d2 = dv.norm2();
610         cvm::real const dinv = cvm::integer_power(d2, -1*(exponent/2));
611         x.real_value += dinv;
612         cvm::rvector const dsumddv = -1.0*(exponent/2) * dinv/d2 * 2.0 * dv;
613         ai1->grad += -1.0 * dsumddv;
614         ai2->grad +=        dsumddv;
615       }
616     }
617   } else {
618     for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
619       for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
620         cvm::rvector const dv = cvm::position_distance(ai1->pos, ai2->pos);
621         cvm::real const d2 = dv.norm2();
622         cvm::real const dinv = cvm::integer_power(d2, -1*(exponent/2));
623         x.real_value += dinv;
624         cvm::rvector const dsumddv = -1.0*(exponent/2) * dinv/d2 * 2.0 * dv;
625         ai1->grad += -1.0 * dsumddv;
626         ai2->grad +=        dsumddv;
627       }
628     }
629   }
630
631   x.real_value *= 1.0 / cvm::real(group1->size() * group2->size());
632   x.real_value = std::pow(x.real_value, -1.0/cvm::real(exponent));
633
634   cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) *
635     cvm::integer_power(x.real_value, exponent+1) /
636     cvm::real(group1->size() * group2->size());
637   for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
638     ai1->grad *= dxdsum;
639   }
640   for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
641     ai2->grad *= dxdsum;
642   }
643 }
644
645
646 void colvar::distance_inv::calc_gradients()
647 {
648 }
649
650
651 void colvar::distance_inv::apply_force(colvarvalue const &force)
652 {
653   if (!group1->noforce)
654     group1->apply_colvar_force(force.real_value);
655
656   if (!group2->noforce)
657     group2->apply_colvar_force(force.real_value);
658 }
659
660
661 simple_scalar_dist_functions(distance_inv)
662
663
664
665 colvar::distance_pairs::distance_pairs(std::string const &conf)
666   : cvc(conf)
667 {
668   function_type = "distance_pairs";
669
670   group1 = parse_group(conf, "group1");
671   group2 = parse_group(conf, "group2");
672
673   x.type(colvarvalue::type_vector);
674   enable(f_cvc_implicit_gradient);
675   x.vector1d_value.resize(group1->size() * group2->size());
676 }
677
678
679 colvar::distance_pairs::distance_pairs()
680 {
681   function_type = "distance_pairs";
682   enable(f_cvc_implicit_gradient);
683   x.type(colvarvalue::type_vector);
684 }
685
686
687 void colvar::distance_pairs::calc_value()
688 {
689   x.vector1d_value.resize(group1->size() * group2->size());
690
691   if (!is_enabled(f_cvc_pbc_minimum_image)) {
692     size_t i1, i2;
693     for (i1 = 0; i1 < group1->size(); i1++) {
694       for (i2 = 0; i2 < group2->size(); i2++) {
695         cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos;
696         cvm::real const d = dv.norm();
697         x.vector1d_value[i1*group2->size() + i2] = d;
698         (*group1)[i1].grad = -1.0 * dv.unit();
699         (*group2)[i2].grad =  dv.unit();
700       }
701     }
702   } else {
703     size_t i1, i2;
704     for (i1 = 0; i1 < group1->size(); i1++) {
705       for (i2 = 0; i2 < group2->size(); i2++) {
706         cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos,
707                                                        (*group2)[i2].pos);
708         cvm::real const d = dv.norm();
709         x.vector1d_value[i1*group2->size() + i2] = d;
710         (*group1)[i1].grad = -1.0 * dv.unit();
711         (*group2)[i2].grad =  dv.unit();
712       }
713     }
714   }
715 }
716
717
718 void colvar::distance_pairs::calc_gradients()
719 {
720   // will be calculated on the fly in apply_force()
721 }
722
723
724 void colvar::distance_pairs::apply_force(colvarvalue const &force)
725 {
726   if (!is_enabled(f_cvc_pbc_minimum_image)) {
727     size_t i1, i2;
728     for (i1 = 0; i1 < group1->size(); i1++) {
729       for (i2 = 0; i2 < group2->size(); i2++) {
730         cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos;
731         (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit());
732         (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit());
733       }
734     }
735   } else {
736     size_t i1, i2;
737     for (i1 = 0; i1 < group1->size(); i1++) {
738       for (i2 = 0; i2 < group2->size(); i2++) {
739         cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos,
740                                                        (*group2)[i2].pos);
741         (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit());
742         (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit());
743       }
744     }
745   }
746 }
747
748
749
750 colvar::dipole_magnitude::dipole_magnitude(std::string const &conf)
751   : cvc(conf)
752 {
753   function_type = "dipole_magnitude";
754   atoms = parse_group(conf, "atoms");
755   init_total_force_params(conf);
756   x.type(colvarvalue::type_scalar);
757 }
758
759 colvar::dipole_magnitude::dipole_magnitude(cvm::atom const &a1)
760 {
761   atoms = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
762   register_atom_group(atoms);
763   x.type(colvarvalue::type_scalar);
764 }
765
766
767 colvar::dipole_magnitude::dipole_magnitude()
768 {
769   function_type = "dipole_magnitude";
770   x.type(colvarvalue::type_scalar);
771 }
772
773
774 void colvar::dipole_magnitude::calc_value()
775 {
776   cvm::atom_pos const atomsCom = atoms->center_of_mass();
777   atoms->calc_dipole(atomsCom);
778   dipoleV = atoms->dipole();
779   x.real_value = dipoleV.norm();
780 }
781
782 void colvar::dipole_magnitude::calc_gradients()
783 {
784   cvm::real const      aux1 = atoms->total_charge/atoms->total_mass;
785   cvm::atom_pos const  dipVunit= dipoleV.unit();
786
787   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
788     ai->grad = (ai->charge - aux1*ai->mass) * dipVunit;
789   }
790 }
791
792 void colvar::dipole_magnitude::apply_force (colvarvalue const &force)
793 {
794   if (!atoms->noforce) {
795     atoms->apply_colvar_force (force.real_value);
796     for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
797       cvm::atom_pos forceV=force.real_value*ai->grad;
798     }
799   }
800 }
801
802
803 simple_scalar_dist_functions(dipole_magnitude)
804
805
806
807 colvar::gyration::gyration(std::string const &conf)
808   : cvc(conf)
809 {
810   function_type = "gyration";
811   provide(f_cvc_inv_gradient);
812   provide(f_cvc_Jacobian);
813   atoms = parse_group(conf, "atoms");
814
815   if (atoms->b_user_defined_fit) {
816     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
817   } else {
818     atoms->b_center = true;
819     atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0));
820     atoms->fit_gradients.assign(atoms->size(), cvm::rvector(0.0, 0.0, 0.0));
821   }
822
823   x.type(colvarvalue::type_scalar);
824 }
825
826
827 colvar::gyration::gyration()
828 {
829   function_type = "gyration";
830   provide(f_cvc_inv_gradient);
831   provide(f_cvc_Jacobian);
832   x.type(colvarvalue::type_scalar);
833 }
834
835
836 void colvar::gyration::calc_value()
837 {
838   x.real_value = 0.0;
839   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
840     x.real_value += (ai->pos).norm2();
841   }
842   x.real_value = std::sqrt(x.real_value / cvm::real(atoms->size()));
843 }
844
845
846 void colvar::gyration::calc_gradients()
847 {
848   cvm::real const drdx = 1.0/(cvm::real(atoms->size()) * x.real_value);
849   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
850     ai->grad = drdx * ai->pos;
851   }
852 }
853
854
855 void colvar::gyration::calc_force_invgrads()
856 {
857   atoms->read_total_forces();
858
859   cvm::real const dxdr = 1.0/x.real_value;
860   ft.real_value = 0.0;
861
862   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
863     ft.real_value += dxdr * ai->pos * ai->total_force;
864   }
865 }
866
867
868 void colvar::gyration::calc_Jacobian_derivative()
869 {
870   jd = x.real_value ? (3.0 * cvm::real(atoms->size()) - 4.0) / x.real_value : 0.0;
871 }
872
873
874 void colvar::gyration::apply_force(colvarvalue const &force)
875 {
876   if (!atoms->noforce)
877     atoms->apply_colvar_force(force.real_value);
878 }
879
880
881 simple_scalar_dist_functions(gyration)
882
883
884
885 colvar::inertia::inertia(std::string const &conf)
886   : gyration(conf)
887 {
888   function_type = "inertia";
889   x.type(colvarvalue::type_scalar);
890 }
891
892
893 colvar::inertia::inertia()
894 {
895   function_type = "inertia";
896   x.type(colvarvalue::type_scalar);
897 }
898
899
900 void colvar::inertia::calc_value()
901 {
902   x.real_value = 0.0;
903   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
904     x.real_value += (ai->pos).norm2();
905   }
906 }
907
908
909 void colvar::inertia::calc_gradients()
910 {
911   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
912     ai->grad = 2.0 * ai->pos;
913   }
914 }
915
916
917 void colvar::inertia::apply_force(colvarvalue const &force)
918 {
919   if (!atoms->noforce)
920     atoms->apply_colvar_force(force.real_value);
921 }
922
923
924 simple_scalar_dist_functions(inertia_z)
925
926
927
928 colvar::inertia_z::inertia_z(std::string const &conf)
929   : inertia(conf)
930 {
931   function_type = "inertia_z";
932   if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) {
933     if (axis.norm2() == 0.0) {
934       cvm::error("Axis vector is zero!");
935       return;
936     }
937     if (axis.norm2() != 1.0) {
938       axis = axis.unit();
939       cvm::log("The normalized axis is: "+cvm::to_str(axis)+".\n");
940     }
941   }
942   x.type(colvarvalue::type_scalar);
943 }
944
945
946 colvar::inertia_z::inertia_z()
947 {
948   function_type = "inertia_z";
949   x.type(colvarvalue::type_scalar);
950 }
951
952
953 void colvar::inertia_z::calc_value()
954 {
955   x.real_value = 0.0;
956   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
957     cvm::real const iprod = ai->pos * axis;
958     x.real_value += iprod * iprod;
959   }
960 }
961
962
963 void colvar::inertia_z::calc_gradients()
964 {
965   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
966     ai->grad = 2.0 * (ai->pos * axis) * axis;
967   }
968 }
969
970
971 void colvar::inertia_z::apply_force(colvarvalue const &force)
972 {
973   if (!atoms->noforce)
974     atoms->apply_colvar_force(force.real_value);
975 }
976
977
978 simple_scalar_dist_functions(inertia)
979
980
981
982
983 colvar::rmsd::rmsd(std::string const &conf)
984   : cvc(conf)
985 {
986   provide(f_cvc_inv_gradient);
987   function_type = "rmsd";
988   x.type(colvarvalue::type_scalar);
989
990   atoms = parse_group(conf, "atoms");
991
992   if (!atoms || atoms->size() == 0) {
993     cvm::error("Error: \"atoms\" must contain at least 1 atom to compute RMSD.");
994     return;
995   }
996
997   bool b_Jacobian_derivative = true;
998   if (atoms->fitting_group != NULL && b_Jacobian_derivative) {
999     cvm::log("The option \"fittingGroup\" (alternative group for fitting) was enabled: "
1000               "Jacobian derivatives of the RMSD will not be calculated.\n");
1001     b_Jacobian_derivative = false;
1002   }
1003   if (b_Jacobian_derivative) provide(f_cvc_Jacobian);
1004
1005   // the following is a simplified version of the corresponding atom group options;
1006   // we need this because the reference coordinates defined inside the atom group
1007   // may be used only for fitting, and even more so if fitting_group is used
1008   if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {
1009     cvm::log("Using reference positions from configuration file to calculate the variable.\n");
1010     if (ref_pos.size() != atoms->size()) {
1011       cvm::error("Error: the number of reference positions provided ("+
1012                   cvm::to_str(ref_pos.size())+
1013                   ") does not match the number of atoms of group \"atoms\" ("+
1014                   cvm::to_str(atoms->size())+").\n");
1015       return;
1016     }
1017   }
1018   {
1019     std::string ref_pos_file;
1020     if (get_keyval(conf, "refPositionsFile", ref_pos_file, std::string(""))) {
1021
1022       if (ref_pos.size()) {
1023         cvm::error("Error: cannot specify \"refPositionsFile\" and "
1024                           "\"refPositions\" at the same time.\n");
1025         return;
1026       }
1027
1028       std::string ref_pos_col;
1029       double ref_pos_col_value=0.0;
1030
1031       if (get_keyval(conf, "refPositionsCol", ref_pos_col, std::string(""))) {
1032         // if provided, use PDB column to select coordinates
1033         bool found = get_keyval(conf, "refPositionsColValue", ref_pos_col_value, 0.0);
1034         if (found && ref_pos_col_value==0.0) {
1035           cvm::error("Error: refPositionsColValue, "
1036                      "if provided, must be non-zero.\n");
1037           return;
1038         }
1039       }
1040
1041       ref_pos.resize(atoms->size());
1042
1043       cvm::load_coords(ref_pos_file.c_str(), &ref_pos, atoms,
1044                        ref_pos_col, ref_pos_col_value);
1045     }
1046   }
1047
1048   if (ref_pos.size() != atoms->size()) {
1049     cvm::error("Error: found " + cvm::to_str(ref_pos.size()) +
1050                     " reference positions; expected " + cvm::to_str(atoms->size()));
1051     return;
1052   }
1053
1054   if (atoms->b_user_defined_fit) {
1055     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
1056   } else {
1057     // Default: fit everything
1058     cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: "
1059               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
1060     atoms->b_center = true;
1061     atoms->b_rotate = true;
1062     // default case: reference positions for calculating the rmsd are also those used
1063     // for fitting
1064     atoms->ref_pos = ref_pos;
1065     atoms->center_ref_pos();
1066
1067     cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation "
1068               "will not be computed as they cancel out in the gradients.");
1069     atoms->disable(f_ag_fit_gradients);
1070
1071     // request the calculation of the derivatives of the rotation defined by the atom group
1072     atoms->rot.request_group1_gradients(atoms->size());
1073     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
1074     // this is only required for ABF, but we do both groups here for better caching
1075     atoms->rot.request_group2_gradients(atoms->size());
1076   }
1077 }
1078
1079
1080 void colvar::rmsd::calc_value()
1081 {
1082   // rotational-translational fit is handled by the atom group
1083
1084   x.real_value = 0.0;
1085   for (size_t ia = 0; ia < atoms->size(); ia++) {
1086     x.real_value += ((*atoms)[ia].pos - ref_pos[ia]).norm2();
1087   }
1088   x.real_value /= cvm::real(atoms->size()); // MSD
1089   x.real_value = std::sqrt(x.real_value);
1090 }
1091
1092
1093 void colvar::rmsd::calc_gradients()
1094 {
1095   cvm::real const drmsddx2 = (x.real_value > 0.0) ?
1096     0.5 / (x.real_value * cvm::real(atoms->size())) :
1097     0.0;
1098
1099   for (size_t ia = 0; ia < atoms->size(); ia++) {
1100     (*atoms)[ia].grad = (drmsddx2 * 2.0 * ((*atoms)[ia].pos - ref_pos[ia]));
1101   }
1102 }
1103
1104
1105 void colvar::rmsd::apply_force(colvarvalue const &force)
1106 {
1107   if (!atoms->noforce)
1108     atoms->apply_colvar_force(force.real_value);
1109 }
1110
1111
1112 void colvar::rmsd::calc_force_invgrads()
1113 {
1114   atoms->read_total_forces();
1115   ft.real_value = 0.0;
1116
1117   // Note: gradient square norm is 1/N_atoms
1118
1119   for (size_t ia = 0; ia < atoms->size(); ia++) {
1120     ft.real_value += (*atoms)[ia].grad * (*atoms)[ia].total_force;
1121   }
1122   ft.real_value *= atoms->size();
1123 }
1124
1125
1126 void colvar::rmsd::calc_Jacobian_derivative()
1127 {
1128   // divergence of the rotated coordinates (including only derivatives of the rotation matrix)
1129   cvm::real rotation_term = 0.0;
1130
1131   // The rotation term only applies is coordinates are rotated
1132   if (atoms->b_rotate) {
1133
1134     // gradient of the rotation matrix
1135     cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);
1136     // gradients of products of 2 quaternion components
1137     cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
1138     for (size_t ia = 0; ia < atoms->size(); ia++) {
1139
1140       // Gradient of optimal quaternion wrt current Cartesian position
1141       cvm::vector1d<cvm::rvector> &dq = atoms->rot.dQ0_1[ia];
1142
1143       g11 = 2.0 * (atoms->rot.q)[1]*dq[1];
1144       g22 = 2.0 * (atoms->rot.q)[2]*dq[2];
1145       g33 = 2.0 * (atoms->rot.q)[3]*dq[3];
1146       g01 = (atoms->rot.q)[0]*dq[1] + (atoms->rot.q)[1]*dq[0];
1147       g02 = (atoms->rot.q)[0]*dq[2] + (atoms->rot.q)[2]*dq[0];
1148       g03 = (atoms->rot.q)[0]*dq[3] + (atoms->rot.q)[3]*dq[0];
1149       g12 = (atoms->rot.q)[1]*dq[2] + (atoms->rot.q)[2]*dq[1];
1150       g13 = (atoms->rot.q)[1]*dq[3] + (atoms->rot.q)[3]*dq[1];
1151       g23 = (atoms->rot.q)[2]*dq[3] + (atoms->rot.q)[3]*dq[2];
1152
1153       // Gradient of the rotation matrix wrt current Cartesian position
1154       grad_rot_mat[0][0] = -2.0 * (g22 + g33);
1155       grad_rot_mat[1][0] =  2.0 * (g12 + g03);
1156       grad_rot_mat[2][0] =  2.0 * (g13 - g02);
1157       grad_rot_mat[0][1] =  2.0 * (g12 - g03);
1158       grad_rot_mat[1][1] = -2.0 * (g11 + g33);
1159       grad_rot_mat[2][1] =  2.0 * (g01 + g23);
1160       grad_rot_mat[0][2] =  2.0 * (g02 + g13);
1161       grad_rot_mat[1][2] =  2.0 * (g23 - g01);
1162       grad_rot_mat[2][2] = -2.0 * (g11 + g22);
1163
1164       cvm::atom_pos &y = ref_pos[ia];
1165
1166       for (size_t alpha = 0; alpha < 3; alpha++) {
1167         for (size_t beta = 0; beta < 3; beta++) {
1168           rotation_term += grad_rot_mat[beta][alpha][alpha] * y[beta];
1169         // Note: equation was derived for inverse rotation (see colvars paper)
1170         // so here the matrix is transposed
1171         // (eq would give   divergence += grad_rot_mat[alpha][beta][alpha] * y[beta];)
1172         }
1173       }
1174     }
1175   }
1176
1177   // The translation term only applies is coordinates are centered
1178   cvm::real translation_term = atoms->b_center ? 3.0 : 0.0;
1179
1180   jd.real_value = x.real_value > 0.0 ?
1181     (3.0 * atoms->size() - 1.0 - translation_term - rotation_term) / x.real_value :
1182     0.0;
1183 }
1184
1185
1186 simple_scalar_dist_functions(rmsd)
1187
1188
1189
1190 colvar::eigenvector::eigenvector(std::string const &conf)
1191   : cvc(conf)
1192 {
1193   provide(f_cvc_inv_gradient);
1194   provide(f_cvc_Jacobian);
1195   function_type = "eigenvector";
1196   x.type(colvarvalue::type_scalar);
1197
1198   atoms = parse_group(conf, "atoms");
1199
1200   {
1201     bool const b_inline = get_keyval(conf, "refPositions", ref_pos, ref_pos);
1202
1203     if (b_inline) {
1204       cvm::log("Using reference positions from input file.\n");
1205       if (ref_pos.size() != atoms->size()) {
1206         cvm::error("Error: reference positions do not "
1207                    "match the number of requested atoms.\n");
1208         return;
1209       }
1210     }
1211
1212     std::string file_name;
1213     if (get_keyval(conf, "refPositionsFile", file_name)) {
1214
1215       if (b_inline) {
1216         cvm::error("Error: refPositions and refPositionsFile cannot be specified at the same time.\n");
1217         return;
1218       }
1219
1220       std::string file_col;
1221       double file_col_value=0.0;
1222       if (get_keyval(conf, "refPositionsCol", file_col, std::string(""))) {
1223         // use PDB flags if column is provided
1224         bool found = get_keyval(conf, "refPositionsColValue", file_col_value, 0.0);
1225         if (found && file_col_value==0.0) {
1226           cvm::error("Error: refPositionsColValue, "
1227                             "if provided, must be non-zero.\n");
1228           return;
1229         }
1230       }
1231
1232       ref_pos.resize(atoms->size());
1233       cvm::load_coords(file_name.c_str(), &ref_pos, atoms,
1234                        file_col, file_col_value);
1235     }
1236   }
1237
1238   if (ref_pos.size() == 0) {
1239     cvm::error("Error: reference positions were not provided.\n", INPUT_ERROR);
1240     return;
1241   }
1242
1243   if (ref_pos.size() != atoms->size()) {
1244     cvm::error("Error: reference positions do not "
1245                "match the number of requested atoms.\n", INPUT_ERROR);
1246     return;
1247   }
1248
1249   // save for later the geometric center of the provided positions (may not be the origin)
1250   cvm::rvector ref_pos_center(0.0, 0.0, 0.0);
1251   for (size_t i = 0; i < atoms->size(); i++) {
1252     ref_pos_center += ref_pos[i];
1253   }
1254   ref_pos_center *= 1.0 / atoms->size();
1255
1256   if (atoms->b_user_defined_fit) {
1257     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
1258   } else {
1259     // default: fit everything
1260     cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: "
1261               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
1262     atoms->b_center = true;
1263     atoms->b_rotate = true;
1264     atoms->ref_pos = ref_pos;
1265     atoms->center_ref_pos();
1266     atoms->disable(f_ag_fit_gradients); // cancel out if group is fitted on itself
1267                                         // and cvc is translationally invariant
1268
1269     // request the calculation of the derivatives of the rotation defined by the atom group
1270     atoms->rot.request_group1_gradients(atoms->size());
1271     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
1272     // this is only required for ABF, but we do both groups here for better caching
1273     atoms->rot.request_group2_gradients(atoms->size());
1274   }
1275
1276   {
1277     bool const b_inline = get_keyval(conf, "vector", eigenvec, eigenvec);
1278     // now load the eigenvector
1279     if (b_inline) {
1280       cvm::log("Using vector components from input file.\n");
1281       if (eigenvec.size() != atoms->size()) {
1282         cvm::error("Error: vector components do not "
1283                           "match the number of requested atoms->\n");
1284         return;
1285       }
1286     }
1287
1288     std::string file_name;
1289     if (get_keyval(conf, "vectorFile", file_name)) {
1290
1291       if (b_inline) {
1292         cvm::error("Error: vector and vectorFile cannot be specified at the same time.\n");
1293         return;
1294       }
1295
1296       std::string file_col;
1297       double file_col_value=0.0;
1298       if (get_keyval(conf, "vectorCol", file_col, std::string(""))) {
1299         // use PDB flags if column is provided
1300         bool found = get_keyval(conf, "vectorColValue", file_col_value, 0.0);
1301         if (found && file_col_value==0.0) {
1302           cvm::error("Error: vectorColValue, if provided, must be non-zero.\n");
1303           return;
1304         }
1305       }
1306
1307       eigenvec.resize(atoms->size());
1308       cvm::load_coords(file_name.c_str(), &eigenvec, atoms,
1309                        file_col, file_col_value);
1310     }
1311   }
1312
1313   if (!ref_pos.size() || !eigenvec.size()) {
1314     cvm::error("Error: both reference coordinates"
1315                       "and eigenvector must be defined.\n");
1316     return;
1317   }
1318
1319   cvm::atom_pos eig_center(0.0, 0.0, 0.0);
1320   for (size_t eil = 0; eil < atoms->size(); eil++) {
1321     eig_center += eigenvec[eil];
1322   }
1323   eig_center *= 1.0 / atoms->size();
1324   cvm::log("Geometric center of the provided vector: "+cvm::to_str(eig_center)+"\n");
1325
1326   bool b_difference_vector = false;
1327   get_keyval(conf, "differenceVector", b_difference_vector, false);
1328
1329   if (b_difference_vector) {
1330
1331     if (atoms->b_center) {
1332       // both sets should be centered on the origin for fitting
1333       for (size_t i = 0; i < atoms->size(); i++) {
1334         eigenvec[i] -= eig_center;
1335         ref_pos[i]  -= ref_pos_center;
1336       }
1337     }
1338     if (atoms->b_rotate) {
1339       atoms->rot.calc_optimal_rotation(eigenvec, ref_pos);
1340       for (size_t i = 0; i < atoms->size(); i++) {
1341         eigenvec[i] = atoms->rot.rotate(eigenvec[i]);
1342       }
1343     }
1344     cvm::log("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n");
1345     for (size_t i = 0; i < atoms->size(); i++) {
1346       eigenvec[i] -= ref_pos[i];
1347     }
1348     if (atoms->b_center) {
1349       // bring back the ref positions to where they were
1350       for (size_t i = 0; i < atoms->size(); i++) {
1351         ref_pos[i] += ref_pos_center;
1352       }
1353     }
1354
1355   } else {
1356     cvm::log("Centering the provided vector to zero.\n");
1357     for (size_t i = 0; i < atoms->size(); i++) {
1358       eigenvec[i] -= eig_center;
1359     }
1360   }
1361
1362   // cvm::log("The first three components(v1x, v1y, v1z) of the resulting vector are: "+cvm::to_str (eigenvec[0])+".\n");
1363
1364   // for inverse gradients
1365   eigenvec_invnorm2 = 0.0;
1366   for (size_t ein = 0; ein < atoms->size(); ein++) {
1367     eigenvec_invnorm2 += eigenvec[ein].norm2();
1368   }
1369   eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
1370
1371   if (b_difference_vector) {
1372     cvm::log("\"differenceVector\" is on: normalizing the vector.\n");
1373     for (size_t i = 0; i < atoms->size(); i++) {
1374       eigenvec[i] *= eigenvec_invnorm2;
1375     }
1376   } else {
1377     cvm::log("The norm of the vector is |v| = "+cvm::to_str(eigenvec_invnorm2)+".\n");
1378   }
1379 }
1380
1381
1382 void colvar::eigenvector::calc_value()
1383 {
1384   x.real_value = 0.0;
1385   for (size_t i = 0; i < atoms->size(); i++) {
1386     x.real_value += ((*atoms)[i].pos - ref_pos[i]) * eigenvec[i];
1387   }
1388 }
1389
1390
1391 void colvar::eigenvector::calc_gradients()
1392 {
1393   for (size_t ia = 0; ia < atoms->size(); ia++) {
1394     (*atoms)[ia].grad = eigenvec[ia];
1395   }
1396 }
1397
1398
1399 void colvar::eigenvector::apply_force(colvarvalue const &force)
1400 {
1401   if (!atoms->noforce)
1402     atoms->apply_colvar_force(force.real_value);
1403 }
1404
1405
1406 void colvar::eigenvector::calc_force_invgrads()
1407 {
1408   atoms->read_total_forces();
1409   ft.real_value = 0.0;
1410
1411   for (size_t ia = 0; ia < atoms->size(); ia++) {
1412     ft.real_value += eigenvec_invnorm2 * (*atoms)[ia].grad *
1413       (*atoms)[ia].total_force;
1414   }
1415 }
1416
1417
1418 void colvar::eigenvector::calc_Jacobian_derivative()
1419 {
1420   // gradient of the rotation matrix
1421   cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);
1422   cvm::quaternion &quat0 = atoms->rot.q;
1423
1424   // gradients of products of 2 quaternion components
1425   cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
1426
1427   cvm::real sum = 0.0;
1428
1429   for (size_t ia = 0; ia < atoms->size(); ia++) {
1430
1431     // Gradient of optimal quaternion wrt current Cartesian position
1432     // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t
1433     // we can just transpose the derivatives of the direct matrix
1434     cvm::vector1d<cvm::rvector> &dq_1 = atoms->rot.dQ0_1[ia];
1435
1436     g11 = 2.0 * quat0[1]*dq_1[1];
1437     g22 = 2.0 * quat0[2]*dq_1[2];
1438     g33 = 2.0 * quat0[3]*dq_1[3];
1439     g01 = quat0[0]*dq_1[1] + quat0[1]*dq_1[0];
1440     g02 = quat0[0]*dq_1[2] + quat0[2]*dq_1[0];
1441     g03 = quat0[0]*dq_1[3] + quat0[3]*dq_1[0];
1442     g12 = quat0[1]*dq_1[2] + quat0[2]*dq_1[1];
1443     g13 = quat0[1]*dq_1[3] + quat0[3]*dq_1[1];
1444     g23 = quat0[2]*dq_1[3] + quat0[3]*dq_1[2];
1445
1446     // Gradient of the inverse rotation matrix wrt current Cartesian position
1447     // (transpose of the gradient of the direct rotation)
1448     grad_rot_mat[0][0] = -2.0 * (g22 + g33);
1449     grad_rot_mat[0][1] =  2.0 * (g12 + g03);
1450     grad_rot_mat[0][2] =  2.0 * (g13 - g02);
1451     grad_rot_mat[1][0] =  2.0 * (g12 - g03);
1452     grad_rot_mat[1][1] = -2.0 * (g11 + g33);
1453     grad_rot_mat[1][2] =  2.0 * (g01 + g23);
1454     grad_rot_mat[2][0] =  2.0 * (g02 + g13);
1455     grad_rot_mat[2][1] =  2.0 * (g23 - g01);
1456     grad_rot_mat[2][2] = -2.0 * (g11 + g22);
1457
1458     for (size_t i = 0; i < 3; i++) {
1459       for (size_t j = 0; j < 3; j++) {
1460         sum += grad_rot_mat[i][j][i] * eigenvec[ia][j];
1461       }
1462     }
1463   }
1464
1465   jd.real_value = sum * std::sqrt(eigenvec_invnorm2);
1466 }
1467
1468
1469 simple_scalar_dist_functions(eigenvector)
1470
1471
1472
1473 colvar::cartesian::cartesian(std::string const &conf)
1474   : cvc(conf)
1475 {
1476   function_type = "cartesian";
1477
1478   atoms = parse_group(conf, "atoms");
1479
1480   bool use_x, use_y, use_z;
1481   get_keyval(conf, "useX", use_x, true);
1482   get_keyval(conf, "useY", use_y, true);
1483   get_keyval(conf, "useZ", use_z, true);
1484
1485   axes.clear();
1486   if (use_x) axes.push_back(0);
1487   if (use_y) axes.push_back(1);
1488   if (use_z) axes.push_back(2);
1489
1490   if (axes.size() == 0) {
1491     cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n");
1492     return;
1493   }
1494
1495   x.type(colvarvalue::type_vector);
1496   enable(f_cvc_implicit_gradient);
1497   x.vector1d_value.resize(atoms->size() * axes.size());
1498 }
1499
1500
1501 void colvar::cartesian::calc_value()
1502 {
1503   size_t const dim = axes.size();
1504   size_t ia, j;
1505   for (ia = 0; ia < atoms->size(); ia++) {
1506     for (j = 0; j < dim; j++) {
1507       x.vector1d_value[dim*ia + j] = (*atoms)[ia].pos[axes[j]];
1508     }
1509   }
1510 }
1511
1512
1513 void colvar::cartesian::calc_gradients()
1514 {
1515   // we're not using the "grad" member of each
1516   // atom object, because it only can represent the gradient of a
1517   // scalar colvar
1518 }
1519
1520
1521 void colvar::cartesian::apply_force(colvarvalue const &force)
1522 {
1523   size_t const dim = axes.size();
1524   size_t ia, j;
1525   if (!atoms->noforce) {
1526     cvm::rvector f;
1527     for (ia = 0; ia < atoms->size(); ia++) {
1528       for (j = 0; j < dim; j++) {
1529         f[axes[j]] = force.vector1d_value[dim*ia + j];
1530       }
1531       (*atoms)[ia].apply_force(f);
1532     }
1533   }
1534 }
1535