Update Colvars to version 2018-12-14
[namd.git] / colvars / src / colvarmodule.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 COLVARMODULE_H
11 #define COLVARMODULE_H
12
13 #include "colvars_version.h"
14
15 #ifndef COLVARS_DEBUG
16 #define COLVARS_DEBUG false
17 #endif
18
19 /*! \mainpage Main page
20 This is the Developer's documentation for the Collective Variables Module.
21
22 You can browse the class hierarchy or the list of source files.
23  */
24
25 /// \file colvarmodule.h
26 /// \brief Collective variables main module
27 ///
28 /// This file declares the main class for defining and manipulating
29 /// collective variables: there should be only one instance of this
30 /// class, because several variables are made static (i.e. they are
31 /// shared between all object instances) to be accessed from other
32 /// objects.
33
34 #define COLVARS_OK 0
35 #define COLVARS_ERROR   1
36 #define COLVARS_NOT_IMPLEMENTED (1<<1)
37 #define INPUT_ERROR     (1<<2) // out of bounds or inconsistent input
38 #define BUG_ERROR       (1<<3) // Inconsistent state indicating bug
39 #define FILE_ERROR      (1<<4)
40 #define MEMORY_ERROR    (1<<5)
41 #define FATAL_ERROR     (1<<6) // Should be set, or not, together with other bits
42 //#define DELETE_COLVARS  (1<<7) // Instruct the caller to delete cvm
43 #define COLVARS_NO_SUCH_FRAME (1<<8) // Cannot load the requested frame
44
45 #include <iostream>
46 #include <iomanip>
47 #include <fstream>
48 #include <sstream>
49 #include <string>
50 #include <vector>
51 #include <list>
52
53 class colvarparse;
54 class colvar;
55 class colvarbias;
56 class colvarproxy;
57 class colvarscript;
58
59
60 /// \brief Collective variables module (main class)
61 ///
62 /// Class to control the collective variables calculation.  An object
63 /// (usually one) of this class is spawned from the MD program,
64 /// containing all i/o routines and general interface.
65 ///
66 /// At initialization, the colvarmodule object creates a proxy object
67 /// to provide a transparent interface between the MD program and the
68 /// child objects
69 class colvarmodule {
70
71 private:
72
73   /// Impossible to initialize the main object without arguments
74   colvarmodule();
75
76   /// Integer representing the version string (allows comparisons)
77   int version_int;
78
79 public:
80
81   /// Get the version number (higher = more recent)
82   int version_number() const
83   {
84     return version_int;
85   }
86
87   friend class colvarproxy;
88   // TODO colvarscript should be unaware of colvarmodule's internals
89   friend class colvarscript;
90
91   /// Defining an abstract real number allows to switch precision
92   typedef  double    real;
93
94   /// Override std::pow with a product for n integer
95   static inline real integer_power(real const &x, int const n)
96   {
97     // Original code: math_special.h in LAMMPS
98     double yy, ww;
99     if (x == 0.0) return 0.0;
100     int nn = (n > 0) ? n : -n;
101     ww = x;
102     for (yy = 1.0; nn != 0; nn >>= 1, ww *=ww) {
103       if (nn & 1) yy *= ww;
104     }
105     return (n > 0) ? yy : 1.0/yy;
106   }
107
108   /// Residue identifier
109   typedef  int       residue_id;
110
111   class rvector;
112   template <class T> class vector1d;
113   template <class T> class matrix2d;
114   class quaternion;
115   class rotation;
116
117   /// \brief Atom position (different type name from rvector, to make
118   /// possible future PBC-transparent implementations)
119   typedef rvector atom_pos;
120
121   /// \brief 3x3 matrix of real numbers
122   class rmatrix;
123
124   // allow these classes to access protected data
125   class atom;
126   class atom_group;
127   friend class atom;
128   friend class atom_group;
129   typedef std::vector<atom>::iterator       atom_iter;
130   typedef std::vector<atom>::const_iterator atom_const_iter;
131
132   /// Module-wide error state
133   /// see constants at the top of this file
134 private:
135
136   static int errorCode;
137
138 public:
139
140   static void set_error_bits(int code);
141
142   static bool get_error_bit(int code);
143
144   static inline int get_error()
145   {
146     return errorCode;
147   }
148
149   static void clear_error();
150
151
152   /// Current step number
153   static long it;
154   /// Starting step number for this run
155   static long it_restart;
156
157   /// Return the current step number from the beginning of this run
158   static inline long step_relative()
159   {
160     return it - it_restart;
161   }
162
163   /// Return the current step number from the beginning of the whole
164   /// calculation
165   static inline long step_absolute()
166   {
167     return it;
168   }
169
170   /// \brief Finite difference step size (if there is no dynamics, or
171   /// if gradients need to be tested independently from the size of
172   /// dt)
173   static real debug_gradients_step_size;
174
175 private:
176
177   /// Prefix for all output files for this run
178   std::string cvm_output_prefix;
179
180 public:
181   /// Accessor for the above
182   static inline std::string &output_prefix()
183   {
184     colvarmodule *cv = colvarmodule::main();
185     return cv->cvm_output_prefix;
186   }
187
188 private:
189
190   /// Array of collective variables
191   std::vector<colvar *> colvars;
192
193   /// Array of collective variables
194   std::vector<colvar *> colvars_active;
195
196   /// Collective variables to be calculated on different threads;
197   /// colvars with multple items (e.g. multiple active CVCs) are duplicated
198   std::vector<colvar *> colvars_smp;
199   /// Indexes of the items to calculate for each colvar
200   std::vector<int> colvars_smp_items;
201
202   /// Array of named atom groups
203   std::vector<atom_group *> named_atom_groups;
204 public:
205   /// Register a named atom group into named_atom_groups
206   inline void register_named_atom_group(atom_group * ag) {
207     named_atom_groups.push_back(ag);
208   }
209
210   /// Array of collective variables
211   std::vector<colvar *> *variables();
212
213   /* TODO: implement named CVCs
214   /// Array of named (reusable) collective variable components
215   static std::vector<cvc *>     cvcs;
216   /// Named cvcs register themselves at initialization time
217   inline void register_cvc(cvc *p) {
218     cvcs.push_back(p);
219   }
220   */
221
222   /// Collective variables with the active flag on
223   std::vector<colvar *> *variables_active();
224
225   /// Collective variables to be calculated on different threads;
226   /// colvars with multple items (e.g. multiple active CVCs) are duplicated
227   std::vector<colvar *> *variables_active_smp();
228
229   /// Indexes of the items to calculate for each colvar
230   std::vector<int> *variables_active_smp_items();
231
232   /// Array of collective variable biases
233   std::vector<colvarbias *> biases;
234
235   /// Energy of built-in and scripted biases, summed per time-step
236   real total_bias_energy;
237
238 private:
239
240   /// Array of active collective variable biases
241   std::vector<colvarbias *> biases_active_;
242
243 public:
244
245   /// Array of active collective variable biases
246   std::vector<colvarbias *> *biases_active();
247
248   /// \brief Whether debug output should be enabled (compile-time option)
249   static inline bool debug()
250   {
251     return COLVARS_DEBUG;
252   }
253
254   /// \brief How many objects are configured yet?
255   size_t size() const;
256
257   /// \brief Constructor \param config_name Configuration file name
258   /// \param restart_name (optional) Restart file name
259   colvarmodule(colvarproxy *proxy);
260
261   /// Destructor
262   ~colvarmodule();
263
264   /// Actual function called by the destructor
265   int reset();
266
267   /// Open a config file, load its contents, and pass it to config_string()
268   int read_config_file(char const *config_file_name);
269
270   /// \brief Parse a config string assuming it is a complete configuration
271   /// (i.e. calling all parse functions)
272   int read_config_string(std::string const &conf);
273
274   /// \brief Parse a "clean" config string (no comments)
275   int parse_config(std::string &conf);
276
277   /// Get the configuration string read so far (includes comments)
278   std::string const & get_config() const;
279
280   // Parse functions (setup internal data based on a string)
281
282   /// Allow reading from Windows text files using using std::getline
283   /// (which can still be used when the text is produced by Colvars itself)
284   static std::istream & getline(std::istream &is, std::string &line);
285
286   /// Parse the few module's global parameters
287   int parse_global_params(std::string const &conf);
288
289   /// Parse and initialize collective variables
290   int parse_colvars(std::string const &conf);
291
292   /// Parse and initialize collective variable biases
293   int parse_biases(std::string const &conf);
294
295   /// \brief Add new configuration during parsing (e.g. to implement
296   /// back-compatibility); cannot be nested, i.e. conf should not contain
297   /// anything that triggers another call
298   int append_new_config(std::string const &conf);
299
300 private:
301
302   /// Configuration string read so far by the module (includes comments)
303   std::string config_string;
304
305   /// Auto-generated configuration during parsing (e.g. to implement
306   /// back-compatibility)
307   std::string extra_conf;
308
309   /// Parse and initialize collective variable biases of a specific type
310   template <class bias_type>
311   int parse_biases_type(std::string const &conf, char const *keyword);
312
313   /// Test error condition and keyword parsing
314   /// on error, delete new bias
315   bool check_new_bias(std::string &conf, char const *key);
316
317 public:
318
319   /// Return how many variables are defined
320   size_t num_variables() const;
321
322   /// Return how many variables have this feature enabled
323   size_t num_variables_feature(int feature_id) const;
324
325   /// Return how many biases are defined
326   size_t num_biases() const;
327
328   /// Return how many biases have this feature enabled
329   size_t num_biases_feature(int feature_id) const;
330
331   /// Return how many biases of this type are defined
332   size_t num_biases_type(std::string const &type) const;
333
334   /// Return the names of time-dependent biases with forces enabled (ABF,
335   /// metadynamics, etc)
336   std::vector<std::string> const time_dependent_biases() const;
337
338 private:
339   /// Useful wrapper to interrupt parsing if any error occurs
340   int catch_input_errors(int result);
341
342 public:
343
344   // "Setup" functions (change internal data based on related data
345   // from the proxy that may change during program execution)
346   // No additional parsing is done within these functions
347
348   /// (Re)initialize internal data (currently used by LAMMPS)
349   /// Also calls setup() member functions of colvars and biases
350   int setup();
351
352   /// (Re)initialize and (re)read the input state file calling read_restart()
353   int setup_input();
354
355   /// (Re)initialize the output trajectory and state file (does not write it yet)
356   int setup_output();
357
358   /// Read the input restart file
359   std::istream & read_restart(std::istream &is);
360   /// Write the output restart file
361   std::ostream & write_restart(std::ostream &os);
362
363   /// Open a trajectory file if requested (and leave it open)
364   int open_traj_file(std::string const &file_name);
365   /// Close it (note: currently unused)
366   int close_traj_file();
367   /// Write in the trajectory file
368   std::ostream & write_traj(std::ostream &os);
369   /// Write explanatory labels in the trajectory file
370   std::ostream & write_traj_label(std::ostream &os);
371
372   /// Write all trajectory files
373   int write_traj_files();
374   /// Write a state file useful to resume the simulation
375   int write_restart_file(std::string const &out_name);
376   /// Write all other output files
377   int write_output_files();
378   /// Backup a file before writing it
379   static int backup_file(char const *filename);
380
381   /// Look up a bias by name; returns NULL if not found
382   static colvarbias * bias_by_name(std::string const &name);
383
384   /// Look up a colvar by name; returns NULL if not found
385   static colvar * colvar_by_name(std::string const &name);
386
387   /// Look up a named atom group by name; returns NULL if not found
388   static atom_group * atom_group_by_name(std::string const &name);
389
390   /// Load new configuration for the given bias -
391   /// currently works for harmonic (force constant and/or centers)
392   int change_configuration(std::string const &bias_name, std::string const &conf);
393
394   /// Read a colvar value
395   std::string read_colvar(std::string const &name);
396
397   /// Calculate change in energy from using alt. config. for the given bias -
398   /// currently works for harmonic (force constant and/or centers)
399   real energy_difference(std::string const &bias_name, std::string const &conf);
400
401   /// Give the total number of bins for a given bias.
402   int bias_bin_num(std::string const &bias_name);
403   /// Calculate the bin index for a given bias.
404   int bias_current_bin(std::string const &bias_name);
405   //// Give the count at a given bin index.
406   int bias_bin_count(std::string const &bias_name, size_t bin_index);
407   //// Share among replicas.
408   int bias_share(std::string const &bias_name);
409
410   /// Main worker function
411   int calc();
412
413   /// Calculate collective variables
414   int calc_colvars();
415
416   /// Calculate biases
417   int calc_biases();
418
419   /// Integrate bias and restraint forces, send colvar forces to atoms
420   int update_colvar_forces();
421
422   /// Perform analysis
423   int analyze();
424
425   /// Carry out operations needed before next step is run
426   int end_of_step();
427
428   /// \brief Read a collective variable trajectory (post-processing
429   /// only, not called at runtime)
430   int read_traj(char const *traj_filename,
431                 long        traj_read_begin,
432                 long        traj_read_end);
433
434   /// Quick conversion of an object to a string
435   template<typename T> static std::string to_str(T const &x,
436                                                   size_t const &width = 0,
437                                                   size_t const &prec = 0);
438   /// Quick conversion of a vector of objects to a string
439   template<typename T> static std::string to_str(std::vector<T> const &x,
440                                                   size_t const &width = 0,
441                                                   size_t const &prec = 0);
442
443   /// Reduce the number of characters in a string
444   static inline std::string wrap_string(std::string const &s,
445                                          size_t const &nchars)
446   {
447     if (!s.size())
448       return std::string(nchars, ' ');
449     else
450       return ( (s.size() <= size_t(nchars)) ?
451                (s+std::string(nchars-s.size(), ' ')) :
452                (std::string(s, 0, nchars)) );
453   }
454
455   /// Number of characters to represent a time step
456   static size_t const it_width;
457   /// Number of digits to represent a collective variables value(s)
458   static size_t const cv_prec;
459   /// Number of characters to represent a collective variables value(s)
460   static size_t const cv_width;
461   /// Number of digits to represent the collective variables energy
462   static size_t const en_prec;
463   /// Number of characters to represent the collective variables energy
464   static size_t const en_width;
465   /// Line separator in the log output
466   static const char * const line_marker;
467
468
469   // proxy functions
470
471   /// \brief Value of the unit for atomic coordinates with respect to
472   /// angstroms (used by some variables for hard-coded default values)
473   static real unit_angstrom();
474
475   /// \brief Boltmann constant
476   static real boltzmann();
477
478   /// \brief Temperature of the simulation (K)
479   static real temperature();
480
481   /// \brief Time step of MD integrator (fs)
482   static real dt();
483
484   /// Request calculation of total force from MD engine
485   static void request_total_force();
486
487   /// Print a message to the main log
488   static void log(std::string const &message);
489
490   /// Print a message to the main log and exit with error code
491   static int fatal_error(std::string const &message);
492
493   /// Print a message to the main log and set global error code
494   static int error(std::string const &message, int code = COLVARS_ERROR);
495
496   // Replica exchange commands.
497   static bool replica_enabled();
498   static int replica_index();
499   static int replica_num();
500   static void replica_comm_barrier();
501   static int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
502   static int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
503
504   /// \brief Get the distance between two atomic positions with pbcs handled
505   /// correctly
506   static rvector position_distance(atom_pos const &pos1,
507                                    atom_pos const &pos2);
508
509   /// \brief Names of groups from a Gromacs .ndx file to be read at startup
510   std::list<std::string> index_group_names;
511
512   /// \brief Groups from a Gromacs .ndx file read at startup
513   std::list<std::vector<int> > index_groups;
514
515   /// \brief Read a Gromacs .ndx file
516   int read_index_file(char const *filename);
517
518   /// \brief Select atom IDs from a file (usually PDB) \param filename name of
519   /// the file \param atoms array into which atoms read from "filename" will be
520   /// appended \param pdb_field (optional) if the file is a PDB and this
521   /// string is non-empty, select atoms for which this field is non-zero
522   /// \param pdb_field_value (optional) if non-zero, select only atoms whose
523   /// pdb_field equals this
524   static int load_atoms(char const *filename,
525                         atom_group &atoms,
526                         std::string const &pdb_field,
527                         double pdb_field_value = 0.0);
528
529   /// \brief Load coordinates for a group of atoms from a file (PDB or XYZ);
530   /// if "pos" is already allocated, the number of its elements must match the
531   /// number of entries in "filename" \param filename name of the file \param
532   /// pos array of coordinates \param atoms group containing the atoms (used
533   /// to obtain internal IDs) \param pdb_field (optional) if the file is a PDB
534   /// and this string is non-empty, select atoms for which this field is
535   /// non-zero \param pdb_field_value (optional) if non-zero, select only
536   /// atoms whose pdb_field equals this
537   static int load_coords(char const *filename,
538                          std::vector<rvector> *pos,
539                          atom_group *atoms,
540                          std::string const &pdb_field,
541                          double pdb_field_value = 0.0);
542
543   /// \brief Load the coordinates for a group of atoms from an
544   /// XYZ file
545   static int load_coords_xyz(char const *filename,
546                              std::vector<rvector> *pos,
547                              atom_group *atoms);
548
549   /// Frequency for collective variables trajectory output
550   static size_t cv_traj_freq;
551
552   /// Frequency for saving output restarts
553   static size_t restart_out_freq;
554   /// Output restart file name
555   std::string   restart_out_name;
556
557   /// Pseudo-random number with Gaussian distribution
558   static real rand_gaussian(void);
559
560 protected:
561
562   /// Configuration file
563   std::ifstream config_s;
564
565   /// Configuration file parser object
566   colvarparse *parse;
567
568   /// Name of the trajectory file
569   std::string cv_traj_name;
570
571   /// Collective variables output trajectory file
572   std::ostream *cv_traj_os;
573
574   /// Appending to the existing trajectory file?
575   bool cv_traj_append;
576
577   /// Write labels at the next iteration
578   bool cv_traj_write_labels;
579
580 private:
581
582   /// Counter for the current depth in the object hierarchy (useg e.g. in output)
583   size_t depth_s;
584
585   /// Thread-specific depth
586   std::vector<size_t> depth_v;
587
588 public:
589
590   /// Get the current object depth in the hierarchy
591   static size_t & depth();
592
593   /// Increase the depth (number of indentations in the output)
594   static void increase_depth();
595
596   /// Decrease the depth (number of indentations in the output)
597   static void decrease_depth();
598
599   static inline bool scripted_forces()
600   {
601     return use_scripted_forces;
602   }
603
604   /// Use scripted colvars forces?
605   static bool use_scripted_forces;
606
607   /// Wait for all biases before calculating scripted forces?
608   static bool scripting_after_biases;
609
610   /// Calculate the energy and forces of scripted biases
611   int calc_scripted_forces();
612
613   /// \brief Pointer to the proxy object, used to retrieve atomic data
614   /// from the hosting program; it is static in order to be accessible
615   /// from static functions in the colvarmodule class
616   static colvarproxy *proxy;
617
618   /// \brief Access the one instance of the Colvars module
619   static colvarmodule *main();
620
621 };
622
623
624 /// Shorthand for the frequently used type prefix
625 typedef colvarmodule cvm;
626
627
628
629 std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
630 std::istream & operator >> (std::istream &is, cvm::rvector &v);
631
632
633 template<typename T> std::string cvm::to_str(T const &x,
634                                              size_t const &width,
635                                              size_t const &prec) {
636   std::ostringstream os;
637   if (width) os.width(width);
638   if (prec) {
639     os.setf(std::ios::scientific, std::ios::floatfield);
640     os.precision(prec);
641   }
642   os << x;
643   return os.str();
644 }
645
646
647 template<typename T> std::string cvm::to_str(std::vector<T> const &x,
648                                              size_t const &width,
649                                              size_t const &prec) {
650   if (!x.size()) return std::string("");
651   std::ostringstream os;
652   if (prec) {
653     os.setf(std::ios::scientific, std::ios::floatfield);
654   }
655   os << "{ ";
656   if (width) os.width(width);
657   if (prec) os.precision(prec);
658   os << x[0];
659   for (size_t i = 1; i < x.size(); i++) {
660     os << ", ";
661     if (width) os.width(width);
662     if (prec) os.precision(prec);
663     os << x[i];
664   }
665   os << " }";
666   return os.str();
667 }
668
669
670 #endif