Update Colvars to version 2018-05-15
[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 Whether r/r0 or \vec{r}*\vec{1/r0_vec} should be used
838   bool b_anisotropic;
839   /// Integer exponent of the function numerator
840   int en;
841   /// Integer exponent of the function denominator
842   int ed;
843
844   /// \brief If true, group2 will be treated as a single atom, stored in this
845   /// accessory group
846   cvm::atom_group *group2_center;
847
848   /// Tolerance for the pair list
849   cvm::real tolerance;
850
851   /// Frequency of update of the pair list
852   int pairlist_freq;
853
854   /// Pair list
855   bool *pairlist = NULL;
856
857 public:
858
859   coordnum(std::string const &conf);
860   coordnum();
861   ~coordnum();
862
863   virtual void calc_value();
864   virtual void calc_gradients();
865   virtual void apply_force(colvarvalue const &force);
866   virtual cvm::real dist2(colvarvalue const &x1,
867                           colvarvalue const &x2) const;
868   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
869                                   colvarvalue const &x2) const;
870   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
871                                   colvarvalue const &x2) const;
872
873   enum {
874     ef_null = 0,
875     ef_gradients = 1,
876     ef_anisotropic = (1<<8),
877     ef_use_pairlist = (1<<9),
878     ef_rebuild_pairlist = (1<<10)
879   };
880
881   /// \brief Calculate a coordination number through the function
882   /// (1-x**n)/(1-x**m), where x = |A1-A2|/r0 \param r0, r0_vec "cutoff" for
883   /// the coordination number (scalar or vector depending on user choice)
884   /// \param en Numerator exponent \param ed Denominator exponent \param First
885   /// atom \param Second atom \param pairlist_elem pointer to pair flag for
886   /// this pair \param tolerance A pair is defined as having a larger
887   /// coordination than this number
888   template<int flags>
889   static cvm::real switching_function(cvm::real const &r0,
890                                       cvm::rvector const &r0_vec,
891                                       int en,
892                                       int ed,
893                                       cvm::atom &A1,
894                                       cvm::atom &A2,
895                                       bool **pairlist_elem,
896                                       cvm::real tolerance);
897
898   /// Main workhorse function
899   template<int flags> int compute_coordnum();
900
901 };
902
903
904
905 /// \brief Colvar component: self-coordination number within a group
906 /// (colvarvalue::type_scalar type, range [0:N*(N-1)/2])
907 class colvar::selfcoordnum
908   : public colvar::cvc
909 {
910 protected:
911
912   /// Selected atoms
913   cvm::atom_group  *group1;
914   /// \brief "Cutoff" for isotropic calculation (default)
915   cvm::real     r0;
916   /// Integer exponent of the function numerator
917   int en;
918   /// Integer exponent of the function denominator
919   int ed;
920   cvm::real tolerance;
921   int pairlist_freq;
922   bool *pairlist = NULL;
923
924 public:
925
926   selfcoordnum(std::string const &conf);
927   selfcoordnum();
928   ~selfcoordnum();
929   virtual void calc_value();
930   virtual void calc_gradients();
931   virtual void apply_force(colvarvalue const &force);
932
933   virtual cvm::real dist2(colvarvalue const &x1,
934                           colvarvalue const &x2) const;
935   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
936                                   colvarvalue const &x2) const;
937   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
938                                   colvarvalue const &x2) const;
939
940   /// Main workhorse function
941   template<int flags> int compute_selfcoordnum();
942 };
943
944
945
946 /// \brief Colvar component: coordination number between two groups
947 /// (colvarvalue::type_scalar type, range [0:N1*N2])
948 class colvar::groupcoordnum
949   : public colvar::distance
950 {
951 protected:
952   /// \brief "Cutoff" for isotropic calculation (default)
953   cvm::real     r0;
954   /// \brief "Cutoff vector" for anisotropic calculation
955   cvm::rvector  r0_vec;
956   /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
957   /// used
958   bool b_anisotropic;
959   /// Integer exponent of the function numerator
960   int en;
961   /// Integer exponent of the function denominator
962   int ed;
963 public:
964   /// Constructor
965   groupcoordnum(std::string const &conf);
966   groupcoordnum();
967   virtual ~groupcoordnum() {}
968   virtual void calc_value();
969   virtual void calc_gradients();
970   virtual void apply_force(colvarvalue const &force);
971
972   virtual cvm::real dist2(colvarvalue const &x1,
973                           colvarvalue const &x2) const;
974   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
975                                   colvarvalue const &x2) const;
976   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
977                                   colvarvalue const &x2) const;
978 };
979
980
981
982 /// \brief Colvar component: hydrogen bond, defined as the product of
983 /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
984 /// (colvarvalue::type_scalar type, range [0:1])
985 class colvar::h_bond
986   : public colvar::cvc
987 {
988 protected:
989   /// \brief "Cutoff" distance between acceptor and donor
990   cvm::real     r0;
991   /// Integer exponent of the function numerator
992   int en;
993   /// Integer exponent of the function denominator
994   int ed;
995 public:
996   h_bond(std::string const &conf);
997   /// Constructor for atoms already allocated
998   h_bond(cvm::atom const &acceptor,
999          cvm::atom const &donor,
1000          cvm::real r0, int en, int ed);
1001   h_bond();
1002   virtual ~h_bond();
1003   virtual void calc_value();
1004   virtual void calc_gradients();
1005   virtual void apply_force(colvarvalue const &force);
1006
1007   virtual cvm::real dist2(colvarvalue const &x1,
1008                           colvarvalue const &x2) const;
1009   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1010                                   colvarvalue const &x2) const;
1011   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1012                                   colvarvalue const &x2) const;
1013 };
1014
1015
1016
1017 /// \brief Colvar component: alpha helix content of a contiguous
1018 /// segment of 5 or more residues, implemented as a sum of phi/psi
1019 /// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type,
1020 /// range [0:1])
1021 // class colvar::alpha_dihedrals
1022 //   : public colvar::cvc
1023 // {
1024 // protected:
1025
1026 //   /// Alpha-helical reference phi value
1027 //   cvm::real phi_ref;
1028
1029 //   /// Alpha-helical reference psi value
1030 //   cvm::real psi_ref;
1031
1032 //   /// List of phi dihedral angles
1033 //   std::vector<dihedral *> phi;
1034
1035 //   /// List of psi dihedral angles
1036 //   std::vector<dihedral *> psi;
1037
1038 //   /// List of hydrogen bonds
1039 //   std::vector<h_bond *>   hb;
1040
1041 // public:
1042
1043 //   alpha_dihedrals (std::string const &conf);
1044 //   alpha_dihedrals();
1045 //   virtual ~alpha_dihedrals() {}
1046 //   virtual void calc_value();
1047 //   virtual void calc_gradients();
1048 //   virtual void apply_force (colvarvalue const &force);
1049 //   virtual cvm::real dist2 (colvarvalue const &x1,
1050 //                            colvarvalue const &x2) const;
1051 //   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
1052 //                                    colvarvalue const &x2) const;
1053 //   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
1054 //                                    colvarvalue const &x2) const;
1055 // };
1056
1057
1058
1059 /// \brief Colvar component: alpha helix content of a contiguous
1060 /// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca
1061 /// angles and hydrogen bonds (colvarvalue::type_scalar type, range
1062 /// [0:1])
1063 class colvar::alpha_angles
1064   : public colvar::cvc
1065 {
1066 protected:
1067
1068   /// Reference Calpha-Calpha angle (default: 88 degrees)
1069   cvm::real theta_ref;
1070
1071   /// Tolerance on the Calpha-Calpha angle
1072   cvm::real theta_tol;
1073
1074   /// List of Calpha-Calpha angles
1075   std::vector<angle *> theta;
1076
1077   /// List of hydrogen bonds
1078   std::vector<h_bond *>   hb;
1079
1080   /// Contribution of the hb terms
1081   cvm::real hb_coeff;
1082
1083 public:
1084
1085   alpha_angles(std::string const &conf);
1086   alpha_angles();
1087   virtual ~alpha_angles();
1088   void calc_value();
1089   void calc_gradients();
1090   void apply_force(colvarvalue const &force);
1091   virtual cvm::real dist2(colvarvalue const &x1,
1092                           colvarvalue const &x2) const;
1093   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1094                                   colvarvalue const &x2) const;
1095   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1096                                   colvarvalue const &x2) const;
1097 };
1098
1099
1100
1101 /// \brief Colvar component: dihedPC
1102 /// Projection of the config onto a dihedral principal component
1103 /// See e.g. Altis et al., J. Chem. Phys 126, 244111 (2007)
1104 /// Based on a set of 'dihedral' cvcs
1105 class colvar::dihedPC
1106   : public colvar::cvc
1107 {
1108 protected:
1109
1110   std::vector<dihedral *> theta;
1111   std::vector<cvm::real> coeffs;
1112
1113 public:
1114
1115   dihedPC(std::string const &conf);
1116   dihedPC();
1117   virtual  ~dihedPC();
1118   void calc_value();
1119   void calc_gradients();
1120   void apply_force(colvarvalue const &force);
1121   virtual cvm::real dist2(colvarvalue const &x1,
1122                           colvarvalue const &x2) const;
1123   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1124                                   colvarvalue const &x2) const;
1125   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1126                                   colvarvalue const &x2) const;
1127 };
1128
1129
1130
1131 /// \brief Colvar component: orientation in space of an atom group,
1132 /// with respect to a set of reference coordinates
1133 /// (colvarvalue::type_quaternion type, range
1134 /// [-1:1]x[-1:1]x[-1:1]x[-1:1])
1135 class colvar::orientation
1136   : public colvar::cvc
1137 {
1138 protected:
1139
1140   /// Atom group
1141   cvm::atom_group  *          atoms;
1142   /// Center of geometry of the group
1143   cvm::atom_pos              atoms_cog;
1144
1145   /// Reference coordinates
1146   std::vector<cvm::atom_pos> ref_pos;
1147
1148   /// Rotation object
1149   cvm::rotation              rot;
1150
1151   /// \brief This is used to remove jumps in the sign of the
1152   /// quaternion, which may be annoying in the colvars trajectory
1153   cvm::quaternion            ref_quat;
1154
1155 public:
1156
1157   orientation(std::string const &conf);
1158   orientation();
1159   virtual ~orientation() {}
1160   virtual void calc_value();
1161   virtual void calc_gradients();
1162   virtual void apply_force(colvarvalue const &force);
1163   virtual cvm::real dist2(colvarvalue const &x1,
1164                           colvarvalue const &x2) const;
1165   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1166                                   colvarvalue const &x2) const;
1167   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1168                                   colvarvalue const &x2) const;
1169 };
1170
1171
1172
1173 /// \brief Colvar component: angle of rotation with respect to a set
1174 /// of reference coordinates (colvarvalue::type_scalar type, range
1175 /// [0:PI))
1176 class colvar::orientation_angle
1177   : public colvar::orientation
1178 {
1179 public:
1180
1181   orientation_angle(std::string const &conf);
1182   orientation_angle();
1183   virtual ~orientation_angle() {}
1184   virtual void calc_value();
1185   virtual void calc_gradients();
1186   virtual void apply_force(colvarvalue const &force);
1187   virtual cvm::real dist2(colvarvalue const &x1,
1188                           colvarvalue const &x2) const;
1189   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1190                                   colvarvalue const &x2) const;
1191   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1192                                   colvarvalue const &x2) const;
1193 };
1194
1195
1196
1197 /// \brief Colvar component: cosine of the angle of rotation with respect to a set
1198 /// of reference coordinates (colvarvalue::type_scalar type, range
1199 /// [-1:1])
1200 class colvar::orientation_proj
1201   : public colvar::orientation
1202 {
1203 public:
1204
1205   orientation_proj(std::string const &conf);
1206   orientation_proj();
1207   virtual ~orientation_proj() {}
1208   virtual void calc_value();
1209   virtual void calc_gradients();
1210   virtual void apply_force(colvarvalue const &force);
1211   virtual cvm::real dist2(colvarvalue const &x1,
1212                           colvarvalue const &x2) const;
1213   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1214                                   colvarvalue const &x2) const;
1215   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1216                                   colvarvalue const &x2) const;
1217 };
1218
1219
1220
1221 /// \brief Colvar component: projection of the orientation vector onto
1222 /// a predefined axis (colvarvalue::type_scalar type, range [-1:1])
1223 class colvar::tilt
1224   : public colvar::orientation
1225 {
1226 protected:
1227
1228   cvm::rvector axis;
1229
1230 public:
1231
1232   tilt(std::string const &conf);
1233   tilt();
1234   virtual ~tilt() {}
1235   virtual void calc_value();
1236   virtual void calc_gradients();
1237   virtual void apply_force(colvarvalue const &force);
1238   virtual cvm::real dist2(colvarvalue const &x1,
1239                           colvarvalue const &x2) const;
1240   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1241                                   colvarvalue const &x2) const;
1242   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1243                                   colvarvalue const &x2) const;
1244 };
1245
1246
1247
1248 /// \brief Colvar component: angle of rotation around a predefined
1249 /// axis (colvarvalue::type_scalar type, range [-PI:PI])
1250 class colvar::spin_angle
1251   : public colvar::orientation
1252 {
1253 protected:
1254
1255   cvm::rvector axis;
1256
1257 public:
1258
1259   spin_angle(std::string const &conf);
1260   spin_angle();
1261   virtual ~spin_angle() {}
1262   virtual void calc_value();
1263   virtual void calc_gradients();
1264   virtual void apply_force(colvarvalue const &force);
1265   /// Redefined to handle the 2*PI periodicity
1266   virtual cvm::real dist2(colvarvalue const &x1,
1267                           colvarvalue const &x2) const;
1268   /// Redefined to handle the 2*PI periodicity
1269   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1270                                   colvarvalue const &x2) const;
1271   /// Redefined to handle the 2*PI periodicity
1272   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1273                                   colvarvalue const &x2) const;
1274   /// Redefined to handle the 2*PI periodicity
1275   virtual void wrap(colvarvalue &x) const;
1276 };
1277
1278
1279
1280 /// \brief Colvar component: root mean square deviation (RMSD) of a
1281 /// group with respect to a set of reference coordinates; uses \link
1282 /// colvar::orientation \endlink to calculate the rotation matrix
1283 /// (colvarvalue::type_scalar type, range [0:*))
1284 class colvar::rmsd
1285   : public colvar::cvc
1286 {
1287 protected:
1288
1289   /// Atom group
1290   cvm::atom_group  *atoms;
1291
1292   /// Reference coordinates (for RMSD calculation only)
1293   std::vector<cvm::atom_pos>  ref_pos;
1294
1295 public:
1296
1297   /// Constructor
1298   rmsd(std::string const &conf);
1299   virtual ~rmsd() {}
1300   virtual void calc_value();
1301   virtual void calc_gradients();
1302   virtual void calc_force_invgrads();
1303   virtual void calc_Jacobian_derivative();
1304   virtual void apply_force(colvarvalue const &force);
1305   virtual cvm::real dist2(colvarvalue const &x1,
1306                           colvarvalue const &x2) const;
1307   virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
1308                                   colvarvalue const &x2) const;
1309   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
1310                                   colvarvalue const &x2) const;
1311 };
1312
1313
1314
1315 // \brief Colvar component: flat vector of Cartesian coordinates
1316 // Mostly useful to compute scripted colvar values
1317 class colvar::cartesian
1318   : public colvar::cvc
1319 {
1320 protected:
1321   /// Atom group
1322   cvm::atom_group  *atoms;
1323   /// Which Cartesian coordinates to include
1324   std::vector<size_t> axes;
1325 public:
1326   cartesian(std::string const &conf);
1327   cartesian();
1328   virtual ~cartesian() {}
1329   virtual void calc_value();
1330   virtual void calc_gradients();
1331   virtual void apply_force(colvarvalue const &force);
1332 };
1333
1334
1335 // metrics functions for cvc implementations
1336
1337 // simple definitions of the distance functions; these are useful only
1338 // for optimization (the type check performed in the default
1339 // colvarcomp functions is skipped)
1340
1341 // definitions assuming the scalar type
1342
1343 #define simple_scalar_dist_functions(TYPE)                              \
1344                                                                         \
1345                                                                         \
1346   cvm::real colvar::TYPE::dist2(colvarvalue const &x1,                  \
1347                                 colvarvalue const &x2) const            \
1348   {                                                                     \
1349     return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
1350   }                                                                     \
1351                                                                         \
1352                                                                         \
1353   colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1,          \
1354                                         colvarvalue const &x2) const    \
1355   {                                                                     \
1356     return 2.0 * (x1.real_value - x2.real_value);                       \
1357   }                                                                     \
1358                                                                         \
1359                                                                         \
1360   colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1,          \
1361                                         colvarvalue const &x2) const    \
1362   {                                                                     \
1363     return this->dist2_lgrad(x2, x1);                                   \
1364   }                                                                     \
1365                                                                         \
1366
1367 #endif