Update Colvars to version 2018-12-14
[namd.git] / colvars / src / colvar.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 COLVAR_H
11 #define COLVAR_H
12
13 #include <iostream>
14
15 #include "colvarmodule.h"
16 #include "colvarvalue.h"
17 #include "colvarparse.h"
18 #include "colvardeps.h"
19
20 #ifdef LEPTON
21 #include "Lepton.h" // for runtime custom expressions
22 #endif
23
24 /// \brief A collective variable (main class); to be defined, it needs
25 /// at least one object of a derived class of colvar::cvc; it
26 /// calculates and returns a \link colvarvalue \endlink object
27 ///
28 /// This class parses the configuration, defines the behaviour and
29 /// stores the value (\link colvar::x \endlink) and all related data
30 /// of a collective variable.  How the value is calculated is defined
31 /// in \link colvar::cvc \endlink and its derived classes.  The
32 /// \link colvar \endlink object contains pointers to multiple \link
33 /// colvar::cvc \endlink derived objects, which can be combined
34 /// together into one collective variable.  This makes possible to
35 /// implement new collective variables at runtime based on the
36 /// existing ones.  Currently, this possibility is limited to a
37 /// polynomial, using the coefficients cvc::sup_coeff and the
38 /// exponents cvc::sup_np.  In case of non-scalar variables,
39 /// only exponents equal to 1 are accepted.
40 ///
41 /// Please note that most of its members are \link colvarvalue
42 /// \endlink objects, i.e. they can handle different data types
43 /// together, and must all be set to the same type of colvar::value()
44 /// before using them together in assignments or other operations; this is usually done
45 /// automatically in the constructor.  If you add a new member of
46 /// \link colvarvalue \endlink type, you should also add its
47 /// initialization line in the \link colvar \endlink constructor.
48
49 class colvar : public colvarparse, public colvardeps {
50
51 public:
52
53   /// Name
54   std::string name;
55
56   /// \brief Current value (previously set by calc() or by read_traj())
57   colvarvalue const & value() const;
58
59   /// \brief Current actual value (not extended DOF)
60   colvarvalue const & actual_value() const;
61
62   /// \brief Current running average (if calculated as set by analysis flag)
63   colvarvalue const & run_ave() const;
64
65   /// \brief Force constant of the spring
66   cvm::real const & force_constant() const;
67
68   /// \brief Current velocity (previously set by calc() or by read_traj())
69   colvarvalue const & velocity() const;
70
71   /// \brief Current total force (previously obtained from calc() or
72   /// read_traj()).  Note: this is calculated using the atomic forces
73   /// from the last simulation step.
74   ///
75   /// Total atom forces are read from the MD program, the total force
76   /// acting on the collective variable is calculated summing those
77   /// from all colvar components, the bias and walls forces are
78   /// subtracted.
79   colvarvalue const & total_force() const;
80
81   /// \brief Typical fluctuation amplitude for this collective
82   /// variable (e.g. local width of a free energy basin)
83   ///
84   /// In metadynamics calculations, \link colvarbias_meta \endlink,
85   /// this value is used to calculate the width of a gaussian.  In ABF
86   /// calculations, \link colvarbias_abf \endlink, it is used to
87   /// calculate the grid spacing in the direction of this collective
88   /// variable.
89   cvm::real width;
90
91   /// \brief Implementation of the feature list for colvar
92   static std::vector<feature *> cv_features;
93
94   /// \brief Implementation of the feature list accessor for colvar
95   virtual const std::vector<feature *> &features()
96   {
97     return cv_features;
98   }
99   virtual std::vector<feature *> &modify_features()
100   {
101     return cv_features;
102   }
103   static void delete_features() {
104     for (size_t i=0; i < cv_features.size(); i++) {
105       delete cv_features[i];
106     }
107     cv_features.clear();
108   }
109
110   /// Implements possible actions to be carried out
111   /// when a given feature is enabled
112   /// This overloads the base function in colvardeps
113   void do_feature_side_effects(int id);
114
115   /// List of biases that depend on this colvar
116   std::vector<colvarbias *> biases;
117 protected:
118
119
120   /*
121     extended:
122     H = H_{0} + \sum_{i} 1/2*\lambda*(S_i(x(t))-s_i(t))^2 \\
123     + \sum_{i} 1/2*m_i*(ds_i(t)/dt)^2 \\
124     + \sum_{t'<t} W * exp(-1/2*\sum_{i} (s_i(t')-s_i(t))^2/(\delta{}s_i)^2) \\
125     + \sum_{w} (\sum_{i}c_{w,i}s_i(t) - D_w)^M/(\Sigma_w)^M
126
127     normal:
128     H = H_{0} + \sum_{t'<t} W * exp(-1/2*\sum_{i} (S_i(x(t'))-S_i(x(t)))^2/(\delta{}S_i)^2) \\
129     + \sum_{w} (\sum_{i}c_{w,i}S_i(t) - D_w)^M/(\Sigma_w)^M
130
131     output:
132     H = H_{0}   (only output S(x), no forces)
133
134     Here:
135     S(x(t)) = x
136     s(t)    = xr
137     DS = Ds = delta
138   */
139
140
141   /// Value of the colvar
142   colvarvalue x;
143
144   // TODO: implement functionality to treat these
145   // /// Vector of individual values from CVCs
146   // colvarvalue x_cvc;
147
148   // /// Jacobian matrix of individual values from CVCs
149   // colvarvalue dx_cvc;
150
151   /// Cached reported value (x may be manipulated)
152   colvarvalue x_reported;
153
154   /// Finite-difference velocity
155   colvarvalue v_fdiff;
156
157   inline colvarvalue fdiff_velocity(colvarvalue const &xold,
158                                      colvarvalue const &xnew)
159   {
160     // using the gradient of the square distance to calculate the
161     // velocity (non-scalar variables automatically taken into
162     // account)
163     cvm::real dt = cvm::dt();
164     return ( ( (dt > 0.0) ? (1.0/dt) : 1.0 ) *
165              0.5 * dist2_lgrad(xnew, xold) );
166   }
167
168   /// Cached reported velocity
169   colvarvalue v_reported;
170
171   // Options for extended_lagrangian
172   /// Restraint center
173   colvarvalue xr;
174   /// Previous value of the restraint center;
175   colvarvalue prev_xr;
176   /// Velocity of the restraint center
177   colvarvalue vr;
178   /// Previous velocity of the restraint center
179   colvarvalue prev_vr;
180   /// Mass of the restraint center
181   cvm::real ext_mass;
182   /// Restraint force constant
183   cvm::real ext_force_k;
184   /// Friction coefficient for Langevin extended dynamics
185   cvm::real ext_gamma;
186   /// Amplitude of Gaussian white noise for Langevin extended dynamics
187   cvm::real ext_sigma;
188
189   /// \brief Harmonic restraint force
190   colvarvalue fr;
191
192   /// \brief Jacobian force, when Jacobian_force is enabled
193   colvarvalue fj;
194
195   /// Cached reported total force
196   colvarvalue ft_reported;
197
198 public:
199
200
201   /// \brief Bias force; reset_bias_force() should be called before
202   /// the biases are updated
203   colvarvalue fb;
204
205   /// \brief Bias force to the actual value (only useful with extended Lagrangian)
206   colvarvalue fb_actual;
207
208   /// \brief Total \em applied force; fr (if extended_lagrangian
209   /// is defined), fb (if biases are applied) and the walls' forces
210   /// (if defined) contribute to it
211   colvarvalue f;
212
213   /// Applied force at the previous step (to be subtracted from total force if needed)
214   colvarvalue f_old;
215
216   /// \brief Total force, as derived from the atomic trajectory;
217   /// should equal the system force plus \link f \endlink
218   colvarvalue ft;
219
220
221   /// Period, if this variable is periodic
222   cvm::real period;
223   cvm::real wrap_center;
224
225
226   /// \brief Expand the boundaries of multiples of width, to keep the
227   /// value always within range
228   bool expand_boundaries;
229
230   /// \brief Location of the lower boundary
231   colvarvalue lower_boundary;
232   /// \brief Location of the lower wall
233   colvarvalue lower_wall;
234   /// \brief Force constant for the lower boundary potential (|x-xb|^2)
235   cvm::real   lower_wall_k;
236   /// \brief Whether this colvar has a hard lower boundary
237   bool        hard_lower_boundary;
238
239   /// \brief Location of the upper boundary
240   colvarvalue upper_boundary;
241   /// \brief Location of the upper wall
242   colvarvalue upper_wall;
243   /// \brief Force constant for the upper boundary potential (|x-xb|^2)
244   cvm::real   upper_wall_k;
245   /// \brief Whether this colvar has a hard upper boundary
246   bool        hard_upper_boundary;
247
248   /// \brief Is the interval defined by the two boundaries periodic?
249   bool periodic_boundaries() const;
250
251   /// \brief Is the interval defined by the two boundaries periodic?
252   bool periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const;
253
254
255   /// Constructor
256   colvar();
257
258   /// Main init function
259   int init(std::string const &conf);
260
261   /// Parse the CVC configuration and allocate their data
262   int init_components(std::string const &conf);
263
264   /// Parse parameters for custom function with Lepton
265   int init_custom_function(std::string const &conf);
266
267   /// Init defaults for grid options
268   int init_grid_parameters(std::string const &conf);
269
270   /// Init extended Lagrangian parameters
271   int init_extended_Lagrangian(std::string const &conf);
272
273   /// Init output flags
274   int init_output_flags(std::string const &conf);
275
276 private:
277   /// Parse the CVC configuration for all components of a certain type
278   template<typename def_class_name> int init_components_type(std::string const &conf,
279                                                              char const *def_desc,
280                                                              char const *def_config_key);
281
282 public:
283
284   /// Get ready for a run and re-initialize internal data if needed
285   void setup();
286
287   /// Destructor
288   ~colvar();
289
290
291   /// \brief Calculate the colvar's value and related quantities
292   int calc();
293
294   /// Carry out operations needed before next step is run
295   int end_of_step();
296
297   /// \brief Calculate a subset of the colvar components (CVCs) currently active
298   /// (default: all active CVCs)
299   /// Note: both arguments refer to the sect of *active* CVCs, not all CVCs
300   int calc_cvcs(int first = 0, size_t num_cvcs = 0);
301
302   /// Ensure that the selected range of CVCs is consistent
303   int check_cvc_range(int first_cvc, size_t num_cvcs);
304
305   /// \brief Calculate the values of the given subset of CVCs
306   int calc_cvc_values(int first, size_t num_cvcs);
307   /// \brief Same as \link colvar::calc_cvc_values \endlink but for gradients
308   int calc_cvc_gradients(int first, size_t num_cvcs);
309   /// \brief Same as \link colvar::calc_cvc_values \endlink but for total forces
310   int calc_cvc_total_force(int first, size_t num_cvcs);
311   /// \brief Same as \link colvar::calc_cvc_values \endlink but for Jacobian derivatives/forces
312   int calc_cvc_Jacobians(int first, size_t num_cvcs);
313
314   /// \brief Collect quantities from CVCs and update aggregated data for the colvar
315   int collect_cvc_data();
316
317   /// \brief Collect the values of the CVCs
318   int collect_cvc_values();
319   /// \brief Same as \link colvar::collect_cvc_values \endlink but for gradients
320   int collect_cvc_gradients();
321   /// \brief Same as \link colvar::collect_cvc_values \endlink but for total forces
322   int collect_cvc_total_forces();
323   /// \brief Same as \link colvar::collect_cvc_values \endlink but for Jacobian derivatives/forces
324   int collect_cvc_Jacobians();
325   /// \brief Calculate the quantities associated to the colvar (but not to the CVCs)
326   int calc_colvar_properties();
327
328   /// Get the current applied force
329   inline colvarvalue const applied_force() const
330   {
331     if (is_enabled(f_cv_extended_Lagrangian)) {
332       return fr;
333     }
334     return f;
335   }
336
337   /// Set the total biasing force to zero
338   void reset_bias_force();
339
340   /// Add to the total force from biases
341   void add_bias_force(colvarvalue const &force);
342
343   /// Apply a force to the actual value (only meaningful with extended Lagrangian)
344   void add_bias_force_actual_value(colvarvalue const &force);
345
346   /// \brief Collect all forces on this colvar, integrate internal
347   /// equations of motion of internal degrees of freedom; see also
348   /// colvar::communicate_forces()
349   /// return colvar energy if extended Lagrandian active
350   cvm::real update_forces_energy();
351
352   /// \brief Communicate forces (previously calculated in
353   /// colvar::update()) to the external degrees of freedom
354   void communicate_forces();
355
356   /// \brief Enables and disables individual CVCs based on flags
357   int set_cvc_flags(std::vector<bool> const &flags);
358
359   /// \brief Updates the flags in the CVC objects, and their number
360   int update_cvc_flags();
361
362   /// \brief Modify the configuration of CVCs (currently, only base class data)
363   int update_cvc_config(std::vector<std::string> const &confs);
364
365 protected:
366   /// \brief Number of CVC objects with an active flag
367   size_t n_active_cvcs;
368
369   /// Sum of square coefficients for active cvcs
370   cvm::real active_cvc_square_norm;
371
372   /// Update the sum of square coefficients for active cvcs
373   void update_active_cvc_square_norm();
374
375   /// \brief Absolute timestep number when this colvar was last updated
376   int prev_timestep;
377
378 public:
379
380   /// \brief Return the number of CVC objects defined
381   inline size_t num_cvcs() const { return cvcs.size(); }
382
383   /// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags)
384   inline size_t num_active_cvcs() const { return n_active_cvcs; }
385
386   /// \brief Use the internal metrics (as from \link cvc
387   /// \endlink objects) to calculate square distances and gradients
388   ///
389   /// Handles correctly symmetries and periodic boundary conditions
390   cvm::real dist2(colvarvalue const &x1,
391                   colvarvalue const &x2) const;
392
393   /// \brief Use the internal metrics (as from \link cvc
394   /// \endlink objects) to calculate square distances and gradients
395   ///
396   /// Handles correctly symmetries and periodic boundary conditions
397   colvarvalue dist2_lgrad(colvarvalue const &x1,
398                           colvarvalue const &x2) const;
399
400   /// \brief Use the internal metrics (as from \link cvc
401   /// \endlink objects) to calculate square distances and gradients
402   ///
403   /// Handles correctly symmetries and periodic boundary conditions
404   colvarvalue dist2_rgrad(colvarvalue const &x1,
405                           colvarvalue const &x2) const;
406
407   /// \brief Use the internal metrics (as from \link cvc
408   /// \endlink objects) to wrap a value into a standard interval
409   ///
410   /// Handles correctly symmetries and periodic boundary conditions
411   void wrap(colvarvalue &x) const;
412
413
414   /// Read the analysis tasks
415   int parse_analysis(std::string const &conf);
416
417   /// Perform analysis tasks
418   int analyze();
419
420
421   /// Read the value from a collective variable trajectory file
422   std::istream & read_traj(std::istream &is);
423   /// Output formatted values to the trajectory file
424   std::ostream & write_traj(std::ostream &os);
425   /// Write a label to the trajectory file (comment line)
426   std::ostream & write_traj_label(std::ostream &os);
427
428
429   /// Read the collective variable from a restart file
430   std::istream & read_restart(std::istream &is);
431   /// Write the collective variable to a restart file
432   std::ostream & write_restart(std::ostream &os);
433
434   /// Write output files (if defined, e.g. in analysis mode)
435   int write_output_files();
436
437
438 protected:
439   /// Previous value (to calculate velocities during analysis)
440   colvarvalue            x_old;
441
442   /// Value read from the most recent state file (if any)
443   colvarvalue            x_restart;
444
445   /// True if a state file was just read
446   bool                   after_restart;
447
448   /// Time series of values and velocities used in correlation
449   /// functions
450   std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
451   /// Time series of values and velocities used in correlation
452   /// functions (pointers)x
453   std::list< std::list<colvarvalue> >::iterator acf_x_history_p, acf_v_history_p;
454
455   /// Time series of values and velocities used in running averages
456   std::list< std::list<colvarvalue> > x_history;
457   /// Time series of values and velocities used in correlation
458   /// functions (pointers)x
459   std::list< std::list<colvarvalue> >::iterator x_history_p;
460
461   /// \brief Collective variable with which the correlation is
462   /// calculated (default: itself)
463   std::string            acf_colvar_name;
464   /// Length of autocorrelation function (ACF)
465   size_t                 acf_length;
466   /// After how many steps the ACF starts
467   size_t                 acf_offset;
468   /// How many timesteps separate two ACF values
469   size_t                 acf_stride;
470   /// Number of frames for each ACF point
471   size_t                 acf_nframes;
472   /// Normalize the ACF to a maximum value of 1?
473   bool                   acf_normalize;
474   /// ACF values
475   std::vector<cvm::real> acf;
476   /// Name of the file to write the ACF
477   std::string            acf_outfile;
478
479   /// Type of autocorrelation function (ACF)
480   enum acf_type_e {
481     /// Unset type
482     acf_notset,
483     /// Velocity ACF, scalar product between v(0) and v(t)
484     acf_vel,
485     /// Coordinate ACF, scalar product between x(0) and x(t)
486     acf_coor,
487     /// \brief Coordinate ACF, second order Legendre polynomial
488     /// between x(0) and x(t) (does not work with scalar numbers)
489     acf_p2coor
490   };
491
492   /// Type of autocorrelation function (ACF)
493   acf_type_e             acf_type;
494
495   /// \brief Velocity ACF, scalar product between v(0) and v(t)
496   void calc_vel_acf(std::list<colvarvalue> &v_history,
497                     colvarvalue const      &v);
498
499   /// \brief Coordinate ACF, scalar product between x(0) and x(t)
500   /// (does not work with scalar numbers)
501   void calc_coor_acf(std::list<colvarvalue> &x_history,
502                      colvarvalue const      &x);
503
504   /// \brief Coordinate ACF, second order Legendre polynomial between
505   /// x(0) and x(t) (does not work with scalar numbers)
506   void calc_p2coor_acf(std::list<colvarvalue> &x_history,
507                        colvarvalue const      &x);
508
509   /// Calculate the auto-correlation function (ACF)
510   int calc_acf();
511   /// Save the ACF to a file
512   int write_acf(std::ostream &os);
513
514   /// Length of running average series
515   size_t         runave_length;
516   /// Timesteps to skip between two values in the running average series
517   size_t         runave_stride;
518   /// Name of the file to write the running average
519   std::string    runave_outfile;
520   /// File to write the running average
521   std::ostream  *runave_os;
522   /// Current value of the running average
523   colvarvalue    runave;
524   /// Current value of the square deviation from the running average
525   cvm::real      runave_variance;
526
527   /// Calculate the running average and its standard deviation
528   int calc_runave();
529
530   /// If extended Lagrangian active: colvar energies (kinetic and harmonic potential)
531   cvm::real kinetic_energy;
532   cvm::real potential_energy;
533 public:
534
535
536   // collective variable component base class
537   class cvc;
538
539   // list of available collective variable components
540
541   // scalar colvar components
542   class distance;
543   class distance_z;
544   class distance_xy;
545   class polar_theta;
546   class polar_phi;
547   class distance_inv;
548   class distance_pairs;
549   class angle;
550   class dipole_angle;
551   class dihedral;
552   class coordnum;
553   class selfcoordnum;
554   class groupcoordnum;
555   class h_bond;
556   class rmsd;
557   class orientation_angle;
558   class orientation_proj;
559   class tilt;
560   class spin_angle;
561   class gyration;
562   class inertia;
563   class inertia_z;
564   class eigenvector;
565   class alpha_dihedrals;
566   class alpha_angles;
567   class dihedPC;
568
569   // non-scalar components
570   class distance_vec;
571   class distance_dir;
572   class cartesian;
573   class orientation;
574
575 protected:
576
577   /// \brief Array of \link cvc \endlink objects
578   std::vector<cvc *> cvcs;
579
580   /// \brief Flags to enable or disable cvcs at next colvar evaluation
581   std::vector<bool> cvc_flags;
582
583   /// \brief Initialize the sorted list of atom IDs for atoms involved
584   /// in all cvcs (called when enabling f_cv_collect_gradients)
585   void build_atom_list(void);
586
587 private:
588   /// Name of scripted function to be used
589   std::string scripted_function;
590
591   /// Current cvc values in the order requested by script
592   /// when using scriptedFunction
593   std::vector<const colvarvalue *> sorted_cvc_values;
594
595 #ifdef LEPTON
596   /// Vector of evaluators for custom functions using Lepton
597   std::vector<Lepton::CompiledExpression *> value_evaluators;
598
599   /// Vector of evaluators for gradients of custom functions
600   std::vector<Lepton::CompiledExpression *> gradient_evaluators;
601
602   /// Vector of references to cvc values to be passed to Lepton evaluators
603   std::vector<double *> value_eval_var_refs;
604   std::vector<double *> grad_eval_var_refs;
605
606   /// Unused value that is written to when a variable simplifies out of a Lepton expression
607   double dev_null;
608 #endif
609
610 public:
611   /// \brief Sorted array of (zero-based) IDs for all atoms involved
612   std::vector<int> atom_ids;
613
614   /// \brief Array of atomic gradients collected from all cvcs
615   /// with appropriate components, rotations etc.
616   /// For scalar variables only!
617   std::vector<cvm::rvector> atomic_gradients;
618
619   inline size_t n_components() const {
620     return cvcs.size();
621   }
622 };
623
624 inline cvm::real const & colvar::force_constant() const
625 {
626   return ext_force_k;
627 }
628
629 inline colvarvalue const & colvar::value() const
630 {
631   return x_reported;
632 }
633
634 inline colvarvalue const & colvar::actual_value() const
635 {
636   return x;
637 }
638
639 inline colvarvalue const & colvar::run_ave() const
640 {
641   return runave;
642 }
643
644 inline colvarvalue const & colvar::velocity() const
645 {
646   return v_reported;
647 }
648
649
650 inline colvarvalue const & colvar::total_force() const
651 {
652   return ft_reported;
653 }
654
655
656 inline void colvar::add_bias_force(colvarvalue const &force)
657 {
658   if (cvm::debug()) {
659     cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
660   }
661   fb += force;
662 }
663
664
665 inline void colvar::add_bias_force_actual_value(colvarvalue const &force)
666 {
667   if (cvm::debug()) {
668     cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
669   }
670   fb_actual += force;
671 }
672
673
674 inline void colvar::reset_bias_force() {
675   fb.type(value());
676   fb.reset();
677   fb_actual.type(value());
678   fb_actual.reset();
679 }
680
681 #endif
682