b5c3b16098c94f3fa94140b1259f016ef6581187
[namd.git] / colvars / src / colvarcomp.h
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 #ifndef COLVARCOMP_H
11 #define COLVARCOMP_H
12
13 // Declaration of colvar::cvc base class and derived ones.
14 //
15 // Future cvc's could be declared on additional header files.
16 // After the declaration of a new derived class, its metric
17 // functions must be reimplemented as well.
18 // If the new cvc has no symmetry or periodicity,
19 // this can be done straightforwardly by using the macro:
20 // simple_scalar_dist_functions (derived_class)
21
22
23 #include "colvarmodule.h"
24 #include "colvar.h"
25 #include "colvaratoms.h"
26
27
28 /// \brief Colvar component (base class for collective variables)
29 ///
30 /// A \link cvc \endlink object (or an object of a
31 /// cvc-derived class) implements the calculation of a collective
32 /// variable, its gradients and any other related physical quantities
33 /// that depend on microscopic degrees of freedom.
34 ///
35 /// No restriction is set to what kind of calculation a \link cvc \endlink
36 /// object performs (usually an analytical function of atomic coordinates).
37 /// The only constraints are that: \par
38 ///
39 /// - The value is calculated by the \link calc_value() \endlink
40 ///   method, and is an object of \link colvarvalue \endlink class.  This
41 ///   provides a transparent way to treat scalar and non-scalar variables
42 ///   alike, and allows an automatic selection of the applicable algorithms.
43 ///
44 /// - The object provides an implementation \link apply_force() \endlink to
45 ///   apply forces to atoms.  Typically, one or more \link cvm::atom_group
46 ///   \endlink objects are used, but this is not a requirement for as long as
47 ///   the \link cvc \endlink object communicates with the simulation program.
48 ///
49 /// <b> If you wish to implement a new collective variable component, you
50 /// should write your own class by inheriting directly from \link
51 /// colvar::cvc \endlink, or one of its derived classes (for instance,
52 /// \link colvar::distance \endlink is frequently used, because it provides
53 /// useful data and function members for any colvar based on two
54 /// atom groups).</b>
55 ///
56 /// The steps are: \par
57 /// 1. Declare the new class as a derivative of \link colvar::cvc \endlink
58 ///    in the file \link colvarcomp.h \endlink
59 /// 2. Implement the new class in a file named colvarcomp_<something>.cpp
60 /// 3. Declare the name of the new class inside the \link colvar \endlink class
61 ///    in \link colvar.h \endlink (see "list of available components")
62 /// 4. Add a call for the new class in colvar::init_components()
63 ////   (file: colvar.cpp)
64 ///
65
66 class colvar::cvc
67   : public colvarparse, public colvardeps
68 {
69 public:
70
71   /// \brief The name of the object (helps to identify this
72   /// cvc instance when debugging)
73   std::string name;
74
75   /// \brief Description of the type of collective variable
76   ///
77   /// Normally this string is set by the parent \link colvar \endlink
78   /// object within its constructor, when all \link cvc \endlink
79   /// objects are initialized; therefore the main "config string"
80   /// constructor does not need to define it.  If a \link cvc
81   /// \endlink is initialized and/or a different constructor is used,
82   /// this variable definition should be set within the constructor.
83   std::string function_type;
84
85   /// \brief Coefficient in the polynomial combination (default: 1.0)
86   cvm::real sup_coeff;
87   /// \brief Exponent in the polynomial combination (default: 1)
88   int       sup_np;
89
90   /// \brief Is this a periodic component?
91   bool b_periodic;
92
93   /// \brief Period of this cvc value, (default: 0.0, non periodic)
94   cvm::real period;
95
96   /// \brief If the component is periodic, wrap around this value (default: 0.0)
97   cvm::real wrap_center;
98
99   /// \brief Constructor
100   ///
101   /// Calls the init() function of the class
102   cvc(std::string const &conf);
103
104   /// An init function should be defined for every class inheriting from cvc
105   /// \param conf Contents of the configuration file pertaining to this \link
106   /// cvc \endlink
107   virtual int init(std::string const &conf);
108
109   /// \brief Within the constructor, make a group parse its own
110   /// options from the provided configuration string
111   /// Returns reference to new group
112   cvm::atom_group *parse_group(std::string const &conf,
113                                char const *group_key,
114                                bool optional = false);
115
116   /// \brief Parse options pertaining to total force calculation
117   virtual int init_total_force_params(std::string const &conf);
118
119   /// \brief After construction, set data related to dependency handling
120   int setup();
121
122   /// \brief Default constructor (used when \link cvc \endlink
123   /// objects are declared within other ones)
124   cvc();
125
126   /// Destructor
127   virtual ~cvc();
128
129   /// \brief Implementation of the feature list for colvar
130   static std::vector<feature *> cvc_features;
131
132   /// \brief Implementation of the feature list accessor for colvar
133   virtual const std::vector<feature *> &features()
134   {
135     return cvc_features;
136   }
137   virtual std::vector<feature *> &modify_features()
138   {
139     return cvc_features;
140   }
141   static void delete_features() {
142     for (size_t i=0; i < cvc_features.size(); i++) {
143       delete cvc_features[i];
144     }
145     cvc_features.clear();
146   }
147
148   /// \brief Obtain data needed for the calculation for the backend
149   virtual void read_data();
150
151   /// \brief Calculate the variable
152   virtual void calc_value() = 0;
153
154   /// \brief Calculate the atomic gradients, to be reused later in
155   /// order to apply forces
156   virtual void calc_gradients() {}
157
158   /// \brief Calculate the atomic fit gradients
159   void calc_fit_gradients();
160
161   /// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component
162   virtual void debug_gradients();
163
164   /// \brief Calculate the total force from the system using the
165   /// inverse atomic gradients
166   virtual void calc_force_invgrads();
167
168   /// \brief Calculate the divergence of the inverse atomic gradients
169   virtual void calc_Jacobian_derivative();
170
171
172   /// \brief Return the previously calculated value
173   colvarvalue const & value() const;
174
175   /// \brief Return the previously calculated total force
176   colvarvalue const & total_force() const;
177
178   /// \brief Return the previously calculated divergence of the
179   /// inverse atomic gradients
180   colvarvalue const & Jacobian_derivative() const;
181
182   /// \brief Apply the collective variable force, by communicating the
183   /// atomic forces to the simulation program (\b Note: the \link ft
184   /// \endlink member is not altered by this function)
185   ///
186   /// Note: multiple calls to this function within the same simulation
187   /// step will add the forces altogether \param cvforce The
188   /// collective variable force, usually coming from the biases and
189   /// eventually manipulated by the parent \link colvar \endlink
190   /// object
191   virtual void apply_force(colvarvalue const &cvforce) = 0;
192
193   /// \brief Square distance between x1 and x2 (can be redefined to
194   /// transparently implement constraints, symmetries and
195   /// periodicities)
196   ///
197   /// colvar::cvc::dist2() and the related functions are
198   /// declared as "const" functions, but not "static", because
199   /// additional parameters defining the metrics (e.g. the
200   /// periodicity) may be specific to each colvar::cvc object.
201   ///
202   /// If symmetries or periodicities are present, the
203   /// colvar::cvc::dist2() should be redefined to return the
204   /// "closest distance" value and colvar::cvc::dist2_lgrad(),
205   /// colvar::cvc::dist2_rgrad() to return its gradients.
206   ///
207   /// If constraints are present (and not already implemented by any
208   /// of the \link colvarvalue \endlink types), the
209   /// colvar::cvc::dist2_lgrad() and
210   /// colvar::cvc::dist2_rgrad() functions should be redefined
211   /// to provide a gradient which is compatible with the constraint,
212   /// i.e. already deprived of its component normal to the constraint
213   /// hypersurface.
214   ///
215   /// Finally, another useful application, if you are performing very
216   /// many operations with these functions, could be to override the
217   /// \link colvarvalue \endlink member functions and access directly
218   /// its member data.  For instance: to define dist2(x1,x2) as
219   /// (x2.real_value-x1.real_value)*(x2.real_value-x1.real_value) in
220   /// case of a scalar \link colvarvalue \endlink type.
221   virtual cvm::real dist2(colvarvalue const &x1,
222                           colvarvalue const &x2) const;
223
224   /// \brief Gradient(with respect to x1) of the square distance (can
225   /// be redefined to transparently implement constraints, symmetries
226   /// and periodicities)
227   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
228                                   colvarvalue const &x2) const;
229
230   /// \brief Gradient(with respect to x2) of the square distance (can
231   /// be redefined to transparently implement constraints, symmetries
232   /// and periodicities)
233   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
234                                   colvarvalue const &x2) const;
235
236   /// \brief Wrap value (for periodic/symmetric cvcs)
237   virtual void wrap(colvarvalue &x) const;
238
239   /// \brief Pointers to all atom groups, to let colvars collect info
240   /// e.g. atomic gradients
241   std::vector<cvm::atom_group *> atom_groups;
242
243   /// \brief Store a pointer to new atom group, and list as child for dependencies
244   inline void register_atom_group(cvm::atom_group *ag) {
245     atom_groups.push_back(ag);
246     add_child((colvardeps *) ag);
247   }
248
249   /// \brief Whether or not this CVC will be computed in parallel whenever possible
250   bool b_try_scalable;
251
252 protected:
253
254   /// \brief Cached value
255   colvarvalue x;
256
257   /// \brief Value at the previous step
258   colvarvalue x_old;
259
260   /// \brief Calculated total force (\b Note: this is calculated from
261   /// the total atomic forces read from the program, subtracting fromt
262   /// the "internal" forces of the system the "external" forces from
263   /// the colvar biases)
264   colvarvalue ft;
265
266   /// \brief Calculated Jacobian derivative (divergence of the inverse
267   /// gradients): serves to calculate the phase space correction
268   colvarvalue jd;
269 };
270
271
272 inline colvarvalue const & colvar::cvc::value() const
273 {
274   return x;
275 }
276
277
278 inline colvarvalue const & colvar::cvc::total_force() const
279 {
280   return ft;
281 }
282
283
284 inline colvarvalue const & colvar::cvc::Jacobian_derivative() const
285 {
286   return jd;
287 }
288
289
290
291 /// \brief Colvar component: distance between the centers of mass of
292 /// two groups (colvarvalue::type_scalar type, range [0:*))
293
294 class colvar::distance
295   : public colvar::cvc
296 {
297 protected:
298   /// First atom group
299   cvm::atom_group  *group1;
300   /// Second atom group
301   cvm::atom_group  *group2;
302   /// Vector distance, cached to be recycled
303   cvm::rvector     dist_v;
304   /// Use absolute positions, ignoring PBCs when present
305   bool b_no_PBC;
306 public:
307   distance(std::string const &conf);
308   distance();
309   virtual ~distance() {}
310   virtual void calc_value();
311   virtual void calc_gradients();
312   virtual void calc_force_invgrads();
313   virtual void calc_Jacobian_derivative();
314   virtual void apply_force(colvarvalue const &force);
315   virtual cvm::real dist2(colvarvalue const &x1,
316                           colvarvalue const &x2) const;
317   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
318                                   colvarvalue const &x2) const;
319   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
320                                   colvarvalue const &x2) const;
321 };
322
323
324
325 // \brief Colvar component: distance vector between centers of mass
326 // of two groups (\link colvarvalue::type_3vector \endlink type,
327 // range (-*:*)x(-*:*)x(-*:*))
328 class colvar::distance_vec
329   : public colvar::distance
330 {
331 public:
332   distance_vec(std::string const &conf);
333   distance_vec();
334   virtual ~distance_vec() {}
335   virtual void calc_value();
336   virtual void calc_gradients();
337   virtual void apply_force(colvarvalue const &force);
338   /// Redefined to handle the box periodicity
339   virtual cvm::real dist2(colvarvalue const &x1,
340                           colvarvalue const &x2) const;
341   /// Redefined to handle the box periodicity
342   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
343                                   colvarvalue const &x2) const;
344   /// Redefined to handle the box periodicity
345   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
346                                   colvarvalue const &x2) const;
347 };
348
349
350
351 /// \brief Colvar component: distance unit vector (direction) between
352 /// centers of mass of two groups (colvarvalue::type_unit3vector type,
353 /// range [-1:1]x[-1:1]x[-1:1])
354 class colvar::distance_dir
355   : public colvar::distance
356 {
357 public:
358   distance_dir(std::string const &conf);
359   distance_dir();
360   virtual ~distance_dir() {}
361   virtual void calc_value();
362   virtual void calc_gradients();
363   virtual void apply_force(colvarvalue const &force);
364   /// Redefined to override the distance ones
365   virtual cvm::real dist2(colvarvalue const &x1,
366                           colvarvalue const &x2) const;
367   /// Redefined to override the distance ones
368   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
369                                   colvarvalue const &x2) const;
370   /// Redefined to override the distance ones
371   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
372                                   colvarvalue const &x2) const;
373 };
374
375
376
377 /// \brief Colvar component: projection of the distance vector along
378 /// an axis(colvarvalue::type_scalar type, range (-*:*))
379 class colvar::distance_z
380   : public colvar::cvc
381 {
382 protected:
383   /// Main atom group
384   cvm::atom_group  *main;
385   /// Reference atom group
386   cvm::atom_group  *ref1;
387   /// Optional, second ref atom group
388   cvm::atom_group  *ref2;
389   /// Use absolute positions, ignoring PBCs when present
390   bool b_no_PBC;
391   /// Vector on which the distance vector is projected
392   cvm::rvector axis;
393   /// Norm of the axis
394   cvm::real axis_norm;
395   /// Vector distance, cached to be recycled
396   cvm::rvector     dist_v;
397   /// Flag: using a fixed axis vector?
398   bool fixed_axis;
399 public:
400   distance_z(std::string const &conf);
401   distance_z();
402   virtual ~distance_z() {}
403   virtual void calc_value();
404   virtual void calc_gradients();
405   virtual void calc_force_invgrads();
406   virtual void calc_Jacobian_derivative();
407   virtual void apply_force(colvarvalue const &force);
408   virtual cvm::real dist2(colvarvalue const &x1,
409                           colvarvalue const &x2) const;
410   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
411                                   colvarvalue const &x2) const;
412   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
413                                   colvarvalue const &x2) const;
414   /// \brief Redefined to make use of the user-provided period
415   virtual void wrap(colvarvalue &x) const;
416 };
417
418
419
420 /// \brief Colvar component: projection of the distance vector on a
421 /// plane (colvarvalue::type_scalar type, range [0:*))
422 class colvar::distance_xy
423   : public colvar::distance_z
424 {
425 protected:
426   /// Components of the distance vector orthogonal to the axis
427   cvm::rvector dist_v_ortho;
428   /// Vector distances
429   cvm::rvector v12, v13;
430 public:
431   distance_xy(std::string const &conf);
432   distance_xy();
433   virtual ~distance_xy() {}
434   virtual void calc_value();
435   virtual void calc_gradients();
436   virtual void calc_force_invgrads();
437   virtual void calc_Jacobian_derivative();
438   virtual void apply_force(colvarvalue const &force);
439   virtual cvm::real dist2(colvarvalue const &x1,
440                           colvarvalue const &x2) const;
441   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
442                                   colvarvalue const &x2) const;
443   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
444                                   colvarvalue const &x2) const;
445 };
446
447
448 /// \brief Colvar component: polar coordinate phi of a group
449 /// (colvarvalue::type_scalar type, range [-180:180])
450 class colvar::polar_phi
451   : public colvar::cvc
452 {
453 public:
454   polar_phi(std::string const &conf);
455   polar_phi();
456   virtual ~polar_phi() {}
457 protected:
458   cvm::atom_group  *atoms;
459   cvm::real r, theta, phi;
460 public:
461   virtual void calc_value();
462   virtual void calc_gradients();
463   virtual void apply_force(colvarvalue const &force);
464   /// Redefined to handle the 2*PI periodicity
465   virtual cvm::real dist2(colvarvalue const &x1,
466                           colvarvalue const &x2) const;
467   /// Redefined to handle the 2*PI periodicity
468   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
469                                   colvarvalue const &x2) const;
470   /// Redefined to handle the 2*PI periodicity
471   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
472                                   colvarvalue const &x2) const;
473   /// Redefined to handle the 2*PI periodicity
474   virtual void wrap(colvarvalue &x) const;
475 };
476
477
478 /// \brief Colvar component: polar coordinate theta of a group
479 /// (colvarvalue::type_scalar type, range [0:180])
480 class colvar::polar_theta
481   : public colvar::cvc
482 {
483 public:
484   polar_theta(std::string const &conf);
485   polar_theta();
486   virtual ~polar_theta() {}
487 protected:
488   cvm::atom_group  *atoms;
489   cvm::real r, theta, phi;
490 public:
491   virtual void calc_value();
492   virtual void calc_gradients();
493   virtual void apply_force(colvarvalue const &force);
494   /// Redefined to override the distance ones
495   virtual cvm::real dist2(colvarvalue const &x1,
496                           colvarvalue const &x2) const;
497   /// Redefined to override the distance ones
498   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
499                                   colvarvalue const &x2) const;
500   /// Redefined to override the distance ones
501   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
502                                   colvarvalue const &x2) const;
503 };
504
505 /// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
506 /// as in NMR refinements(colvarvalue::type_scalar type, range (0:*))
507 class colvar::distance_inv
508   : public colvar::cvc
509 {
510 protected:
511   /// First atom group
512   cvm::atom_group  *group1;
513   /// Second atom group
514   cvm::atom_group  *group2;
515   /// Components of the distance vector orthogonal to the axis
516   int exponent;
517   /// Use absolute positions, ignoring PBCs when present
518   bool b_no_PBC;
519 public:
520   distance_inv(std::string const &conf);
521   distance_inv();
522   virtual ~distance_inv() {}
523   virtual void calc_value();
524   virtual void calc_gradients();
525   virtual void apply_force(colvarvalue const &force);
526   virtual cvm::real dist2(colvarvalue const &x1,
527                           colvarvalue const &x2) const;
528   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
529                                   colvarvalue const &x2) const;
530   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
531                                   colvarvalue const &x2) const;
532 };
533
534
535
536 /// \brief Colvar component: N1xN2 vector of pairwise distances
537 /// (colvarvalue::type_vector type, range (0:*) for each component)
538 class colvar::distance_pairs
539   : public colvar::cvc
540 {
541 protected:
542   /// First atom group
543   cvm::atom_group  *group1;
544   /// Second atom group
545   cvm::atom_group  *group2;
546   /// Use absolute positions, ignoring PBCs when present
547   bool b_no_PBC;
548 public:
549   distance_pairs(std::string const &conf);
550   distance_pairs();
551   virtual ~distance_pairs() {}
552   virtual void calc_value();
553   virtual void calc_gradients();
554   virtual void apply_force(colvarvalue const &force);
555 };
556
557
558
559 /// \brief Colvar component: Radius of gyration of an atom group
560 /// (colvarvalue::type_scalar type, range [0:*))
561 class colvar::gyration
562   : public colvar::cvc
563 {
564 protected:
565   /// Atoms involved
566   cvm::atom_group  *atoms;
567 public:
568   /// Constructor
569   gyration(std::string const &conf);
570   gyration();
571   virtual ~gyration() {}
572   virtual void calc_value();
573   virtual void calc_gradients();
574   virtual void calc_force_invgrads();
575   virtual void calc_Jacobian_derivative();
576   virtual void apply_force(colvarvalue const &force);
577   virtual cvm::real dist2(colvarvalue const &x1,
578                           colvarvalue const &x2) const;
579   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
580                                   colvarvalue const &x2) const;
581   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
582                                   colvarvalue const &x2) const;
583 };
584
585
586
587 /// \brief Colvar component: moment of inertia of an atom group
588 /// (colvarvalue::type_scalar type, range [0:*))
589 class colvar::inertia
590   : public colvar::gyration
591 {
592 public:
593   /// Constructor
594   inertia(std::string const &conf);
595   inertia();
596   virtual ~inertia() {}
597   virtual void calc_value();
598   virtual void calc_gradients();
599   virtual void apply_force(colvarvalue const &force);
600   virtual cvm::real dist2(colvarvalue const &x1,
601                           colvarvalue const &x2) const;
602   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
603                                   colvarvalue const &x2) const;
604   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
605                                   colvarvalue const &x2) const;
606 };
607
608
609
610 /// \brief Colvar component: moment of inertia of an atom group
611 /// around a user-defined axis (colvarvalue::type_scalar type, range [0:*))
612 class colvar::inertia_z
613   : public colvar::inertia
614 {
615 protected:
616   /// Vector on which the inertia tensor is projected
617   cvm::rvector axis;
618 public:
619   /// Constructor
620   inertia_z(std::string const &conf);
621   inertia_z();
622   virtual ~inertia_z() {}
623   virtual void calc_value();
624   virtual void calc_gradients();
625   virtual void apply_force(colvarvalue const &force);
626   virtual cvm::real dist2(colvarvalue const &x1,
627                           colvarvalue const &x2) const;
628   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
629                                   colvarvalue const &x2) const;
630   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
631                                   colvarvalue const &x2) const;
632 };
633
634
635
636 /// \brief Colvar component: projection of 3N coordinates onto an
637 /// eigenvector(colvarvalue::type_scalar type, range (-*:*))
638 class colvar::eigenvector
639   : public colvar::cvc
640 {
641 protected:
642
643   /// Atom group
644   cvm::atom_group  *           atoms;
645
646   /// Reference coordinates
647   std::vector<cvm::atom_pos>  ref_pos;
648
649   /// Geometric center of the reference coordinates
650   cvm::atom_pos                ref_pos_center;
651
652   /// Eigenvector (of a normal or essential mode): will always have zero center
653   std::vector<cvm::rvector>   eigenvec;
654
655   /// Inverse square norm of the eigenvector
656   cvm::real                   eigenvec_invnorm2;
657
658 public:
659
660   /// Constructor
661   eigenvector(std::string const &conf);
662   virtual ~eigenvector() {}
663   virtual void calc_value();
664   virtual void calc_gradients();
665   virtual void calc_force_invgrads();
666   virtual void calc_Jacobian_derivative();
667   virtual void apply_force(colvarvalue const &force);
668   virtual cvm::real dist2(colvarvalue const &x1,
669                           colvarvalue const &x2) const;
670   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
671                                   colvarvalue const &x2) const;
672   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
673                                   colvarvalue const &x2) const;
674 };
675
676
677
678 /// \brief Colvar component: angle between the centers of mass of
679 /// three groups (colvarvalue::type_scalar type, range [0:PI])
680 class colvar::angle
681   : public colvar::cvc
682 {
683 protected:
684
685   /// Atom group
686   cvm::atom_group  *group1;
687   /// Atom group
688   cvm::atom_group  *group2;
689   /// Atom group
690   cvm::atom_group  *group3;
691
692   /// Inter site vectors
693   cvm::rvector r21, r23;
694   /// Inter site vector norms
695   cvm::real r21l, r23l;
696   /// Derivatives wrt group centers of mass
697   cvm::rvector dxdr1, dxdr3;
698
699   /// Compute total force on first site only to avoid unwanted
700   /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
701   /// (or to allow dummy atoms)
702   bool b_1site_force;
703 public:
704
705   /// Initialize by parsing the configuration
706   angle(std::string const &conf);
707   /// \brief Initialize the three groups after three atoms
708   angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
709   angle();
710   virtual ~angle() {}
711   virtual void calc_value();
712   virtual void calc_gradients();
713   virtual void calc_force_invgrads();
714   virtual void calc_Jacobian_derivative();
715   virtual void apply_force(colvarvalue const &force);
716   virtual cvm::real dist2(colvarvalue const &x1,
717                           colvarvalue const &x2) const;
718   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
719                                   colvarvalue const &x2) const;
720   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
721                                   colvarvalue const &x2) const;
722 };
723
724
725
726 /// \brief Colvar component: angle between the dipole of a molecule and an axis
727 /// formed by two groups of atoms(colvarvalue::type_scalar type, range [0:PI])
728 class colvar::dipole_angle
729   : public colvar::cvc
730 {
731 protected:
732
733   /// Dipole atom group
734   cvm::atom_group  *group1;
735   /// Atom group
736   cvm::atom_group  *group2;
737   /// Atom group
738   cvm::atom_group  *group3;
739
740   /// Inter site vectors
741   cvm::rvector r21, r23;
742   /// Inter site vector norms
743   cvm::real r21l, r23l;
744   /// Derivatives wrt group centers of mass
745   cvm::rvector dxdr1, dxdr3;
746
747   /// Compute total force on first site only to avoid unwanted
748   /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
749   /// (or to allow dummy atoms)
750   bool b_1site_force;
751 public:
752
753   /// Initialize by parsing the configuration
754   dipole_angle (std::string const &conf);
755   /// \brief Initialize the three groups after three atoms
756   dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
757   dipole_angle();
758   virtual ~dipole_angle() {}
759   virtual void calc_value();
760   virtual void calc_gradients();
761   virtual void apply_force (colvarvalue const &force);
762   virtual cvm::real dist2 (colvarvalue const &x1,
763                            colvarvalue const &x2) const;
764   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
765                                    colvarvalue const &x2) const;
766   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
767                                    colvarvalue const &x2) const;
768 };
769
770
771
772 /// \brief Colvar component: dihedral between the centers of mass of
773 /// four groups (colvarvalue::type_scalar type, range [-PI:PI])
774 class colvar::dihedral
775   : public colvar::cvc
776 {
777 protected:
778
779   /// Atom group
780   cvm::atom_group  *group1;
781   /// Atom group
782   cvm::atom_group  *group2;
783   /// Atom group
784   cvm::atom_group  *group3;
785   /// Atom group
786   cvm::atom_group  *group4;
787   /// Inter site vectors
788   cvm::rvector r12, r23, r34;
789
790   /// \brief Compute total force on first site only to avoid unwanted
791   /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
792   bool b_1site_force;
793
794 public:
795
796   /// Initialize by parsing the configuration
797   dihedral(std::string  const &conf);
798   /// \brief Initialize the four groups after four atoms
799   dihedral(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4);
800   dihedral();
801   virtual ~dihedral() {}
802   virtual void calc_value();
803   virtual void calc_gradients();
804   virtual void calc_force_invgrads();
805   virtual void calc_Jacobian_derivative();
806   virtual void apply_force(colvarvalue const &force);
807
808   /// Redefined to handle the 2*PI periodicity
809   virtual cvm::real dist2(colvarvalue const &x1,
810                           colvarvalue const &x2) const;
811   /// Redefined to handle the 2*PI periodicity
812   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
813                                   colvarvalue const &x2) const;
814   /// Redefined to handle the 2*PI periodicity
815   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
816                                   colvarvalue const &x2) const;
817   /// Redefined to handle the 2*PI periodicity
818   virtual void wrap(colvarvalue &x) const;
819 };
820
821
822
823 /// \brief Colvar component: coordination number between two groups
824 /// (colvarvalue::type_scalar type, range [0:N1*N2])
825 class colvar::coordnum
826   : public colvar::cvc
827 {
828 protected:
829   /// First atom group
830   cvm::atom_group  *group1;
831   /// Second atom group
832   cvm::atom_group  *group2;
833   /// \brief "Cutoff" for isotropic calculation (default)
834   cvm::real     r0;
835   /// \brief "Cutoff vector" for anisotropic calculation
836   cvm::rvector  r0_vec;
837   /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
838   /// used
839   bool b_anisotropic;
840   /// Integer exponent of the function numerator
841   int en;
842   /// Integer exponent of the function denominator
843   int ed;
844   /// \brief If true, group2 will be treated as a single atom
845   /// (default: loop over all pairs of atoms in group1 and group2)
846   bool b_group2_center_only;
847 public:
848   /// Constructor
849   coordnum(std::string const &conf);
850   coordnum();
851   virtual ~coordnum() {}
852   virtual void calc_value();
853   virtual void calc_gradients();
854   virtual void apply_force(colvarvalue const &force);
855   template<bool b_gradients>
856   /// \brief Calculate a coordination number through the function
857   /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
858   /// coordination number \param exp_num \i n exponent \param exp_den
859   /// \i m exponent \param A1 atom \param A2 atom
860   static cvm::real switching_function(cvm::real const &r0,
861                                       int const &exp_num, int const &exp_den,
862                                       cvm::atom &A1, cvm::atom &A2);
863
864   template<bool b_gradients>
865   /// \brief Calculate a coordination number through the function
866   /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
867   /// vector of different cutoffs in the three directions \param
868   /// exp_num \i n exponent \param exp_den \i m exponent \param A1
869   /// atom \param A2 atom
870   static cvm::real switching_function(cvm::rvector const &r0_vec,
871                                       int const &exp_num, int const &exp_den,
872                                       cvm::atom &A1, cvm::atom &A2);
873
874   virtual cvm::real dist2(colvarvalue const &x1,
875                           colvarvalue const &x2) const;
876   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
877                                   colvarvalue const &x2) const;
878   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
879                                   colvarvalue const &x2) const;
880 };
881
882
883
884 /// \brief Colvar component: self-coordination number within a group
885 /// (colvarvalue::type_scalar type, range [0:N*(N-1)/2])
886 class colvar::selfcoordnum
887   : public colvar::cvc
888 {
889 protected:
890   /// First atom group
891   cvm::atom_group  *group1;
892   /// \brief "Cutoff" for isotropic calculation (default)
893   cvm::real     r0;
894   /// Integer exponent of the function numerator
895   int en;
896   /// Integer exponent of the function denominator
897   int ed;
898 public:
899   /// Constructor
900   selfcoordnum(std::string const &conf);
901   selfcoordnum();
902   virtual ~selfcoordnum() {}
903   virtual void calc_value();
904   virtual void calc_gradients();
905   virtual void apply_force(colvarvalue const &force);
906   template<bool b_gradients>
907   /// \brief Calculate a coordination number through the function
908   /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
909   /// coordination number \param exp_num \i n exponent \param exp_den
910   /// \i m exponent \param A1 atom \param A2 atom
911   static cvm::real switching_function(cvm::real const &r0,
912                                       int const &exp_num, int const &exp_den,
913                                       cvm::atom &A1, cvm::atom &A2);
914
915   virtual cvm::real dist2(colvarvalue const &x1,
916                           colvarvalue const &x2) const;
917   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
918                                   colvarvalue const &x2) const;
919   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
920                                   colvarvalue const &x2) const;
921 };
922
923
924
925 /// \brief Colvar component: coordination number between two groups
926 /// (colvarvalue::type_scalar type, range [0:N1*N2])
927 class colvar::groupcoordnum
928   : public colvar::distance
929 {
930 protected:
931   /// \brief "Cutoff" for isotropic calculation (default)
932   cvm::real     r0;
933   /// \brief "Cutoff vector" for anisotropic calculation
934   cvm::rvector  r0_vec;
935   /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
936   /// used
937   bool b_anisotropic;
938   /// Integer exponent of the function numerator
939   int en;
940   /// Integer exponent of the function denominator
941   int ed;
942 public:
943   /// Constructor
944   groupcoordnum(std::string const &conf);
945   groupcoordnum();
946   virtual ~groupcoordnum() {}
947   virtual void calc_value();
948   virtual void calc_gradients();
949   virtual void apply_force(colvarvalue const &force);
950   template<bool b_gradients>
951   /// \brief Calculate a coordination number through the function
952   /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
953   /// coordination number \param exp_num \i n exponent \param exp_den
954   /// \i m exponent \param A1 atom \param A2 atom
955   static cvm::real switching_function(cvm::real const &r0,
956                                       int const &exp_num, int const &exp_den,
957                                       cvm::atom &A1, cvm::atom &A2);
958
959   /*
960   template<bool b_gradients>
961   /// \brief Calculate a coordination number through the function
962   /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
963   /// vector of different cutoffs in the three directions \param
964   /// exp_num \i n exponent \param exp_den \i m exponent \param A1
965   /// atom \param A2 atom
966   static cvm::real switching_function(cvm::rvector const &r0_vec,
967                                       int const &exp_num, int const &exp_den,
968                                       cvm::atom &A1, cvm::atom &A2);
969   */
970
971   virtual cvm::real dist2(colvarvalue const &x1,
972                           colvarvalue const &x2) const;
973   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
974                                   colvarvalue const &x2) const;
975   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
976                                   colvarvalue const &x2) const;
977 };
978
979
980
981 /// \brief Colvar component: hydrogen bond, defined as the product of
982 /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
983 /// (colvarvalue::type_scalar type, range [0:1])
984 class colvar::h_bond
985   : public colvar::cvc
986 {
987 protected:
988   /// \brief "Cutoff" distance between acceptor and donor
989   cvm::real     r0;
990   /// Integer exponent of the function numerator
991   int en;
992   /// Integer exponent of the function denominator
993   int ed;
994 public:
995   h_bond(std::string const &conf);
996   /// Constructor for atoms already allocated
997   h_bond(cvm::atom const &acceptor,
998          cvm::atom const &donor,
999          cvm::real r0, int en, int ed);
1000   h_bond();
1001   virtual ~h_bond();
1002   virtual void calc_value();
1003   virtual void calc_gradients();
1004   virtual void apply_force(colvarvalue const &force);
1005
1006   virtual cvm::real dist2(colvarvalue const &x1,
1007                           colvarvalue const &x2) const;
1008   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1009                                   colvarvalue const &x2) const;
1010   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1011                                   colvarvalue const &x2) const;
1012 };
1013
1014
1015
1016 /// \brief Colvar component: alpha helix content of a contiguous
1017 /// segment of 5 or more residues, implemented as a sum of phi/psi
1018 /// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type,
1019 /// range [0:1])
1020 // class colvar::alpha_dihedrals
1021 //   : public colvar::cvc
1022 // {
1023 // protected:
1024
1025 //   /// Alpha-helical reference phi value
1026 //   cvm::real phi_ref;
1027
1028 //   /// Alpha-helical reference psi value
1029 //   cvm::real psi_ref;
1030
1031 //   /// List of phi dihedral angles
1032 //   std::vector<dihedral *> phi;
1033
1034 //   /// List of psi dihedral angles
1035 //   std::vector<dihedral *> psi;
1036
1037 //   /// List of hydrogen bonds
1038 //   std::vector<h_bond *>   hb;
1039
1040 // public:
1041
1042 //   alpha_dihedrals (std::string const &conf);
1043 //   alpha_dihedrals();
1044 //   virtual ~alpha_dihedrals() {}
1045 //   virtual void calc_value();
1046 //   virtual void calc_gradients();
1047 //   virtual void apply_force (colvarvalue const &force);
1048 //   virtual cvm::real dist2 (colvarvalue const &x1,
1049 //                            colvarvalue const &x2) const;
1050 //   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
1051 //                                    colvarvalue const &x2) const;
1052 //   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
1053 //                                    colvarvalue const &x2) const;
1054 // };
1055
1056
1057
1058 /// \brief Colvar component: alpha helix content of a contiguous
1059 /// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca
1060 /// angles and hydrogen bonds (colvarvalue::type_scalar type, range
1061 /// [0:1])
1062 class colvar::alpha_angles
1063   : public colvar::cvc
1064 {
1065 protected:
1066
1067   /// Reference Calpha-Calpha angle (default: 88 degrees)
1068   cvm::real theta_ref;
1069
1070   /// Tolerance on the Calpha-Calpha angle
1071   cvm::real theta_tol;
1072
1073   /// List of Calpha-Calpha angles
1074   std::vector<angle *> theta;
1075
1076   /// List of hydrogen bonds
1077   std::vector<h_bond *>   hb;
1078
1079   /// Contribution of the hb terms
1080   cvm::real hb_coeff;
1081
1082 public:
1083
1084   alpha_angles(std::string const &conf);
1085   alpha_angles();
1086   virtual ~alpha_angles();
1087   void calc_value();
1088   void calc_gradients();
1089   void apply_force(colvarvalue const &force);
1090   virtual cvm::real dist2(colvarvalue const &x1,
1091                           colvarvalue const &x2) const;
1092   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1093                                   colvarvalue const &x2) const;
1094   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1095                                   colvarvalue const &x2) const;
1096 };
1097
1098
1099
1100 /// \brief Colvar component: dihedPC
1101 /// Projection of the config onto a dihedral principal component
1102 /// See e.g. Altis et al., J. Chem. Phys 126, 244111 (2007)
1103 /// Based on a set of 'dihedral' cvcs
1104 class colvar::dihedPC
1105   : public colvar::cvc
1106 {
1107 protected:
1108
1109   std::vector<dihedral *> theta;
1110   std::vector<cvm::real> coeffs;
1111
1112 public:
1113
1114   dihedPC(std::string const &conf);
1115   dihedPC();
1116   virtual  ~dihedPC();
1117   void calc_value();
1118   void calc_gradients();
1119   void apply_force(colvarvalue const &force);
1120   virtual cvm::real dist2(colvarvalue const &x1,
1121                           colvarvalue const &x2) const;
1122   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1123                                   colvarvalue const &x2) const;
1124   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1125                                   colvarvalue const &x2) const;
1126 };
1127
1128
1129
1130 /// \brief Colvar component: orientation in space of an atom group,
1131 /// with respect to a set of reference coordinates
1132 /// (colvarvalue::type_quaternion type, range
1133 /// [-1:1]x[-1:1]x[-1:1]x[-1:1])
1134 class colvar::orientation
1135   : public colvar::cvc
1136 {
1137 protected:
1138
1139   /// Atom group
1140   cvm::atom_group  *          atoms;
1141   /// Center of geometry of the group
1142   cvm::atom_pos              atoms_cog;
1143
1144   /// Reference coordinates
1145   std::vector<cvm::atom_pos> ref_pos;
1146
1147   /// Rotation object
1148   cvm::rotation              rot;
1149
1150   /// \brief This is used to remove jumps in the sign of the
1151   /// quaternion, which may be annoying in the colvars trajectory
1152   cvm::quaternion            ref_quat;
1153
1154 public:
1155
1156   orientation(std::string const &conf);
1157   orientation();
1158   virtual ~orientation() {}
1159   virtual void calc_value();
1160   virtual void calc_gradients();
1161   virtual void apply_force(colvarvalue const &force);
1162   virtual cvm::real dist2(colvarvalue const &x1,
1163                           colvarvalue const &x2) const;
1164   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1165                                   colvarvalue const &x2) const;
1166   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1167                                   colvarvalue const &x2) const;
1168 };
1169
1170
1171
1172 /// \brief Colvar component: angle of rotation with respect to a set
1173 /// of reference coordinates (colvarvalue::type_scalar type, range
1174 /// [0:PI))
1175 class colvar::orientation_angle
1176   : public colvar::orientation
1177 {
1178 public:
1179
1180   orientation_angle(std::string const &conf);
1181   orientation_angle();
1182   virtual ~orientation_angle() {}
1183   virtual void calc_value();
1184   virtual void calc_gradients();
1185   virtual void apply_force(colvarvalue const &force);
1186   virtual cvm::real dist2(colvarvalue const &x1,
1187                           colvarvalue const &x2) const;
1188   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1189                                   colvarvalue const &x2) const;
1190   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1191                                   colvarvalue const &x2) const;
1192 };
1193
1194
1195
1196 /// \brief Colvar component: cosine of the angle of rotation with respect to a set
1197 /// of reference coordinates (colvarvalue::type_scalar type, range
1198 /// [-1:1])
1199 class colvar::orientation_proj
1200   : public colvar::orientation
1201 {
1202 public:
1203
1204   orientation_proj(std::string const &conf);
1205   orientation_proj();
1206   virtual ~orientation_proj() {}
1207   virtual void calc_value();
1208   virtual void calc_gradients();
1209   virtual void apply_force(colvarvalue const &force);
1210   virtual cvm::real dist2(colvarvalue const &x1,
1211                           colvarvalue const &x2) const;
1212   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1213                                   colvarvalue const &x2) const;
1214   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1215                                   colvarvalue const &x2) const;
1216 };
1217
1218
1219
1220 /// \brief Colvar component: projection of the orientation vector onto
1221 /// a predefined axis (colvarvalue::type_scalar type, range [-1:1])
1222 class colvar::tilt
1223   : public colvar::orientation
1224 {
1225 protected:
1226
1227   cvm::rvector axis;
1228
1229 public:
1230
1231   tilt(std::string const &conf);
1232   tilt();
1233   virtual ~tilt() {}
1234   virtual void calc_value();
1235   virtual void calc_gradients();
1236   virtual void apply_force(colvarvalue const &force);
1237   virtual cvm::real dist2(colvarvalue const &x1,
1238                           colvarvalue const &x2) const;
1239   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1240                                   colvarvalue const &x2) const;
1241   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1242                                   colvarvalue const &x2) const;
1243 };
1244
1245
1246
1247 /// \brief Colvar component: angle of rotation around a predefined
1248 /// axis (colvarvalue::type_scalar type, range [-PI:PI])
1249 class colvar::spin_angle
1250   : public colvar::orientation
1251 {
1252 protected:
1253
1254   cvm::rvector axis;
1255
1256 public:
1257
1258   spin_angle(std::string const &conf);
1259   spin_angle();
1260   virtual ~spin_angle() {}
1261   virtual void calc_value();
1262   virtual void calc_gradients();
1263   virtual void apply_force(colvarvalue const &force);
1264   /// Redefined to handle the 2*PI periodicity
1265   virtual cvm::real dist2(colvarvalue const &x1,
1266                           colvarvalue const &x2) const;
1267   /// Redefined to handle the 2*PI periodicity
1268   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1269                                   colvarvalue const &x2) const;
1270   /// Redefined to handle the 2*PI periodicity
1271   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1272                                   colvarvalue const &x2) const;
1273   /// Redefined to handle the 2*PI periodicity
1274   virtual void wrap(colvarvalue &x) const;
1275 };
1276
1277
1278
1279 /// \brief Colvar component: root mean square deviation (RMSD) of a
1280 /// group with respect to a set of reference coordinates; uses \link
1281 /// colvar::orientation \endlink to calculate the rotation matrix
1282 /// (colvarvalue::type_scalar type, range [0:*))
1283 class colvar::rmsd
1284   : public colvar::cvc
1285 {
1286 protected:
1287
1288   /// Atom group
1289   cvm::atom_group  *atoms;
1290
1291   /// Reference coordinates (for RMSD calculation only)
1292   std::vector<cvm::atom_pos>  ref_pos;
1293
1294 public:
1295
1296   /// Constructor
1297   rmsd(std::string const &conf);
1298   virtual ~rmsd() {}
1299   virtual void calc_value();
1300   virtual void calc_gradients();
1301   virtual void calc_force_invgrads();
1302   virtual void calc_Jacobian_derivative();
1303   virtual void apply_force(colvarvalue const &force);
1304   virtual cvm::real dist2(colvarvalue const &x1,
1305                           colvarvalue const &x2) const;
1306   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1307                                   colvarvalue const &x2) const;
1308   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1309                                   colvarvalue const &x2) const;
1310 };
1311
1312
1313
1314 // \brief Colvar component: flat vector of Cartesian coordinates
1315 // Mostly useful to compute scripted colvar values
1316 class colvar::cartesian
1317   : public colvar::cvc
1318 {
1319 protected:
1320   /// Atom group
1321   cvm::atom_group  *atoms;
1322   /// Which Cartesian coordinates to include
1323   std::vector<size_t> axes;
1324 public:
1325   cartesian(std::string const &conf);
1326   cartesian();
1327   virtual ~cartesian() {}
1328   virtual void calc_value();
1329   virtual void calc_gradients();
1330   virtual void apply_force(colvarvalue const &force);
1331 };
1332
1333
1334 // metrics functions for cvc implementations
1335
1336 // simple definitions of the distance functions; these are useful only
1337 // for optimization (the type check performed in the default
1338 // colvarcomp functions is skipped)
1339
1340 // definitions assuming the scalar type
1341
1342 #define simple_scalar_dist_functions(TYPE)                              \
1343                                                                         \
1344                                                                         \
1345   cvm::real colvar::TYPE::dist2(colvarvalue const &x1,                  \
1346                                 colvarvalue const &x2) const            \
1347   {                                                                     \
1348     return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
1349   }                                                                     \
1350                                                                         \
1351                                                                         \
1352   colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1,          \
1353                                         colvarvalue const &x2) const    \
1354   {                                                                     \
1355     return 2.0 * (x1.real_value - x2.real_value);                       \
1356   }                                                                     \
1357                                                                         \
1358                                                                         \
1359   colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1,          \
1360                                         colvarvalue const &x2) const    \
1361   {                                                                     \
1362     return this->dist2_lgrad(x2, x1);                                   \
1363   }                                                                     \
1364                                                                         \
1365
1366 #endif