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