Update Colvars to version 2018-12-18
[namd.git] / colvars / src / colvarproxy.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 COLVARPROXY_H
11 #define COLVARPROXY_H
12
13 #include <fstream>
14 #include <list>
15
16 #include "colvarmodule.h"
17 #include "colvartypes.h"
18 #include "colvarvalue.h"
19
20
21 /// \file colvarproxy.h
22 /// \brief Colvars proxy classes
23 ///
24 /// This file declares the class for the object responsible for interfacing
25 /// Colvars with other codes (MD engines, VMD, Python).  The \link colvarproxy
26 /// \endlink class is a derivative of multiple classes, each devoted to a
27 /// specific task (e.g. \link colvarproxy_atoms \endlink to access data for
28 /// individual atoms).
29 ///
30 /// To interface to a new MD engine, the simplest solution is to derive a new
31 /// class from \link colvarproxy \endlink.  Currently implemented are: \link
32 /// colvarproxy_lammps, \endlink, \link colvarproxy_namd, \endlink, \link
33 /// colvarproxy_vmd \endlink.
34
35
36 // forward declarations
37 class colvarscript;
38
39
40 /// Methods for accessing the simulation system (PBCs, integrator, etc)
41 class colvarproxy_system {
42
43 public:
44
45   /// Constructor
46   colvarproxy_system();
47
48   /// Destructor
49   virtual ~colvarproxy_system();
50
51   /// \brief Value of the unit for atomic coordinates with respect to
52   /// angstroms (used by some variables for hard-coded default values)
53   virtual cvm::real unit_angstrom() = 0;
54
55   /// \brief Boltzmann constant
56   virtual cvm::real boltzmann() = 0;
57
58   /// \brief Target temperature of the simulation (K units)
59   virtual cvm::real temperature() = 0;
60
61   /// \brief Time step of the simulation (fs)
62   virtual cvm::real dt() = 0;
63
64   /// \brief Pseudo-random number with Gaussian distribution
65   virtual cvm::real rand_gaussian(void) = 0;
66
67   /// Pass restraint energy value for current timestep to MD engine
68   virtual void add_energy(cvm::real energy) = 0;
69
70   /// \brief Get the PBC-aware distance vector between two positions
71   virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
72                                          cvm::atom_pos const &pos2) const;
73
74   /// Recompute PBC reciprocal lattice (assumes XYZ periodicity)
75   void update_pbc_lattice();
76
77   /// Set the lattice vectors to zero
78   void reset_pbc_lattice();
79
80   /// \brief Tell the proxy whether total forces are needed (they may not
81   /// always be available)
82   virtual void request_total_force(bool yesno);
83
84   /// Are total forces being used?
85   virtual bool total_forces_enabled() const;
86
87   /// Are total forces from the current step available?
88   virtual bool total_forces_same_step() const;
89
90 protected:
91
92   /// \brief Type of boundary conditions
93   ///
94   /// Orthogonal and triclinic cells are made available to objects.
95   /// For any other conditions (mixed periodicity, triclinic cells in LAMMPS)
96   /// minimum-image distances are computed by the host engine regardless.
97   enum Boundaries_type {
98     boundaries_non_periodic,
99     boundaries_pbc_ortho,
100     boundaries_pbc_triclinic,
101     boundaries_unsupported
102   };
103
104   /// Type of boundary conditions
105   Boundaries_type boundaries_type;
106
107   /// Bravais lattice vectors
108   cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z;
109
110   /// Reciprocal lattice vectors
111   cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z;
112 };
113
114
115 /// \brief Container of atomic data for processing by Colvars
116 class colvarproxy_atoms {
117
118 public:
119
120   /// Constructor
121   colvarproxy_atoms();
122
123   /// Destructor
124   virtual ~colvarproxy_atoms();
125
126   /// Prepare this atom for collective variables calculation, selecting it by
127   /// numeric index (1-based)
128   virtual int init_atom(int atom_number) = 0;
129
130   /// Check that this atom number is valid, but do not initialize the
131   /// corresponding atom yet
132   virtual int check_atom_id(int atom_number) = 0;
133
134   /// Select this atom for collective variables calculation, using name and
135   /// residue number.  Not all programs support this: leave this function as
136   /// is in those cases.
137   virtual int init_atom(cvm::residue_id const &residue,
138                         std::string const     &atom_name,
139                         std::string const     &segment_id);
140
141   /// Check that this atom is valid, but do not initialize it yet
142   virtual int check_atom_id(cvm::residue_id const &residue,
143                             std::string const     &atom_name,
144                             std::string const     &segment_id);
145
146   /// \brief Used by the atom class destructor: rather than deleting the array slot
147   /// (costly) set the corresponding atoms_ncopies to zero
148   virtual void clear_atom(int index);
149
150   /// \brief Select atom IDs from a file (usually PDB) \param filename name of
151   /// the file \param atoms array to which atoms read from "filename" will be
152   /// appended \param pdb_field (optional) if the file is a PDB and this
153   /// string is non-empty, select atoms for which this field is non-zero
154   /// \param pdb_field_value (optional) if non-zero, select only atoms whose
155   /// pdb_field equals this
156   virtual int load_atoms(char const *filename,
157                          cvm::atom_group &atoms,
158                          std::string const &pdb_field,
159                          double pdb_field_value = 0.0);
160
161   /// \brief Load a set of coordinates from a file (usually PDB); if "pos" is
162   /// already allocated, the number of its elements must match the number of
163   /// entries in "filename" \param filename name of the file \param pos array
164   /// of coordinates \param sorted_ids array of sorted internal IDs, used to
165   /// loop through the file only once \param pdb_field (optional) if the file
166   /// is a PDB and this string is non-empty, select atoms for which this field
167   /// is non-zero \param pdb_field_value (optional) if non-zero, select only
168   /// atoms whose pdb_field equals this
169   virtual int load_coords(char const *filename,
170                           std::vector<cvm::atom_pos> &pos,
171                           std::vector<int> const &sorted_ids,
172                           std::string const &pdb_field,
173                           double pdb_field_value = 0.0);
174
175   /// Clear atomic data
176   int reset();
177
178   /// Get the numeric ID of the given atom (for the program)
179   inline int get_atom_id(int index) const
180   {
181     return atoms_ids[index];
182   }
183
184   /// Get the mass of the given atom
185   inline cvm::real get_atom_mass(int index) const
186   {
187     return atoms_masses[index];
188   }
189
190   /// Get the charge of the given atom
191   inline cvm::real get_atom_charge(int index) const
192   {
193     return atoms_charges[index];
194   }
195
196   /// Read the current position of the given atom
197   inline cvm::rvector get_atom_position(int index) const
198   {
199     return atoms_positions[index];
200   }
201
202   /// Read the current total force of the given atom
203   inline cvm::rvector get_atom_total_force(int index) const
204   {
205     return atoms_total_forces[index];
206   }
207
208   /// Request that this force is applied to the given atom
209   inline void apply_atom_force(int index, cvm::rvector const &new_force)
210   {
211     atoms_new_colvar_forces[index] += new_force;
212   }
213
214   /// Read the current velocity of the given atom
215   inline cvm::rvector get_atom_velocity(int index)
216   {
217     cvm::error("Error: reading the current velocity of an atom "
218                "is not yet implemented.\n",
219                COLVARS_NOT_IMPLEMENTED);
220     return cvm::rvector(0.0);
221   }
222
223   inline std::vector<int> *modify_atom_ids()
224   {
225     return &atoms_ids;
226   }
227
228   inline std::vector<cvm::real> *modify_atom_masses()
229   {
230     return &atoms_masses;
231   }
232
233   inline std::vector<cvm::real> *modify_atom_charges()
234   {
235     return &atoms_charges;
236   }
237
238   inline std::vector<cvm::rvector> *modify_atom_positions()
239   {
240     return &atoms_positions;
241   }
242
243   inline std::vector<cvm::rvector> *modify_atom_total_forces()
244   {
245     return &atoms_total_forces;
246   }
247
248   inline std::vector<cvm::rvector> *modify_atom_new_colvar_forces()
249   {
250     return &atoms_new_colvar_forces;
251   }
252
253 protected:
254
255   /// \brief Array of 0-based integers used to uniquely associate atoms
256   /// within the host program
257   std::vector<int>          atoms_ids;
258   /// \brief Keep track of how many times each atom is used by a separate colvar object
259   std::vector<size_t>       atoms_ncopies;
260   /// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
261   std::vector<cvm::real>    atoms_masses;
262   /// \brief Charges of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
263   std::vector<cvm::real>    atoms_charges;
264   /// \brief Current three-dimensional positions of the atoms
265   std::vector<cvm::rvector> atoms_positions;
266   /// \brief Most recent total forces on each atom
267   std::vector<cvm::rvector> atoms_total_forces;
268   /// \brief Forces applied from colvars, to be communicated to the MD integrator
269   std::vector<cvm::rvector> atoms_new_colvar_forces;
270
271   /// Used by all init_atom() functions: create a slot for an atom not
272   /// requested yet; returns the index in the arrays
273   int add_atom_slot(int atom_id);
274
275 };
276
277
278 /// \brief Container of atom group data (allow collection of aggregated atomic
279 /// data)
280 class colvarproxy_atom_groups {
281
282 public:
283
284   /// Contructor
285   colvarproxy_atom_groups();
286
287   /// Destructor
288   virtual ~colvarproxy_atom_groups();
289
290   /// Clear atom group data
291   int reset();
292
293   /// \brief Whether this proxy implementation has capability for scalable groups
294   virtual int scalable_group_coms();
295
296   /// Prepare this group for collective variables calculation, selecting atoms by internal ids (0-based)
297   virtual int init_atom_group(std::vector<int> const &atoms_ids);
298
299   /// \brief Used by the atom_group class destructor
300   virtual void clear_atom_group(int index);
301
302   /// Get the numeric ID of the given atom group (for the MD program)
303   inline int get_atom_group_id(int index) const
304   {
305     return atom_groups_ids[index];
306   }
307
308   /// Get the mass of the given atom group
309   inline cvm::real get_atom_group_mass(int index) const
310   {
311     return atom_groups_masses[index];
312   }
313
314   /// Get the charge of the given atom group
315   inline cvm::real get_atom_group_charge(int index) const
316   {
317     return atom_groups_charges[index];
318   }
319
320   /// Read the current position of the center of mass given atom group
321   inline cvm::rvector get_atom_group_com(int index) const
322   {
323     return atom_groups_coms[index];
324   }
325
326   /// Read the current total force of the given atom group
327   inline cvm::rvector get_atom_group_total_force(int index) const
328   {
329     return atom_groups_total_forces[index];
330   }
331
332   /// Request that this force is applied to the given atom group
333   inline void apply_atom_group_force(int index, cvm::rvector const &new_force)
334   {
335     atom_groups_new_colvar_forces[index] += new_force;
336   }
337
338   /// Read the current velocity of the given atom group
339   inline cvm::rvector get_atom_group_velocity(int index)
340   {
341     cvm::error("Error: reading the current velocity of an atom group is not yet implemented.\n",
342                COLVARS_NOT_IMPLEMENTED);
343     return cvm::rvector(0.0);
344   }
345
346 protected:
347
348   /// \brief Array of 0-based integers used to uniquely associate atom groups
349   /// within the host program
350   std::vector<int>          atom_groups_ids;
351   /// \brief Keep track of how many times each group is used by a separate cvc
352   std::vector<size_t>       atom_groups_ncopies;
353   /// \brief Total masses of the atom groups
354   std::vector<cvm::real>    atom_groups_masses;
355   /// \brief Total charges of the atom groups (allow redefinition during a run, as done e.g. in LAMMPS)
356   std::vector<cvm::real>    atom_groups_charges;
357   /// \brief Current centers of mass of the atom groups
358   std::vector<cvm::rvector> atom_groups_coms;
359   /// \brief Most recently updated total forces on the com of each group
360   std::vector<cvm::rvector> atom_groups_total_forces;
361   /// \brief Forces applied from colvars, to be communicated to the MD integrator
362   std::vector<cvm::rvector> atom_groups_new_colvar_forces;
363
364   /// Used by all init_atom_group() functions: create a slot for an atom group not requested yet
365   int add_atom_group_slot(int atom_group_id);
366 };
367
368
369 /// \brief Methods for SMP parallelization
370 class colvarproxy_smp {
371
372 public:
373
374   /// Constructor
375   colvarproxy_smp();
376
377   /// Destructor
378   virtual ~colvarproxy_smp();
379
380   /// Whether threaded parallelization should be used (TODO: make this a
381   /// cvm::deps feature)
382   bool b_smp_active;
383
384   /// Whether threaded parallelization is available (TODO: make this a cvm::deps feature)
385   virtual int smp_enabled();
386
387   /// Distribute calculation of colvars (and their components) across threads
388   virtual int smp_colvars_loop();
389
390   /// Distribute calculation of biases across threads
391   virtual int smp_biases_loop();
392
393   /// Distribute calculation of biases across threads 2nd through last, with all scripted biased on 1st thread
394   virtual int smp_biases_script_loop();
395
396   /// Index of this thread
397   virtual int smp_thread_id();
398
399   /// Number of threads sharing this address space
400   virtual int smp_num_threads();
401
402   /// Lock the proxy's shared data for access by a thread, if threads are implemented; if not implemented, does nothing
403   virtual int smp_lock();
404
405   /// Attempt to lock the proxy's shared data
406   virtual int smp_trylock();
407
408   /// Release the lock
409   virtual int smp_unlock();
410
411 protected:
412
413   /// Lock state for OpenMP
414   void *omp_lock_state;
415 };
416
417
418 /// \brief Methods for multiple-replica communication
419 class colvarproxy_replicas {
420
421 public:
422
423   /// Constructor
424   colvarproxy_replicas();
425
426   /// Destructor
427   virtual ~colvarproxy_replicas();
428
429   /// \brief Indicate if multi-replica support is available and active
430   virtual bool replica_enabled();
431
432   /// \brief Index of this replica
433   virtual int replica_index();
434
435   /// \brief Total number of replica
436   virtual int replica_num();
437
438   /// \brief Synchronize replica
439   virtual void replica_comm_barrier();
440
441   /// \brief Receive data from other replica
442   virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
443
444   /// \brief Send data to other replica
445   virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
446
447 };
448
449
450 /// Methods for scripting language interface (Tcl or Python)
451 class colvarproxy_script {
452
453 public:
454
455   /// Constructor
456   colvarproxy_script();
457
458   /// Destructor
459   virtual ~colvarproxy_script();
460
461   /// Convert a script object (Tcl or Python call argument) to a C string
462   virtual char const *script_obj_to_str(unsigned char *obj);
463
464   /// Convert a script object (Tcl or Python call argument) to a vector of strings
465   virtual std::vector<std::string> script_obj_to_str_vector(unsigned char *obj);
466
467   /// Pointer to the scripting interface object
468   /// (does not need to be allocated in a new interface)
469   colvarscript *script;
470
471   /// is a user force script defined?
472   bool force_script_defined;
473
474   /// Do we have a scripting interface?
475   bool have_scripts;
476
477   /// Run a user-defined colvar forces script
478   virtual int run_force_callback();
479
480   virtual int run_colvar_callback(
481                 std::string const &name,
482                 std::vector<const colvarvalue *> const &cvcs,
483                 colvarvalue &value);
484
485   virtual int run_colvar_gradient_callback(
486                 std::string const &name,
487                 std::vector<const colvarvalue *> const &cvcs,
488                 std::vector<cvm::matrix2d<cvm::real> > &gradient);
489 };
490
491
492 /// Methods for using Tcl within Colvars
493 class colvarproxy_tcl {
494
495 public:
496
497   /// Constructor
498   colvarproxy_tcl();
499
500   /// Destructor
501   virtual ~colvarproxy_tcl();
502
503   /// Is Tcl available? (trigger initialization if needed)
504   int tcl_available();
505
506   /// Tcl implementation of script_obj_to_str()
507   char const *tcl_obj_to_str(unsigned char *obj);
508
509   /// Run a user-defined colvar forces script
510   int tcl_run_force_callback();
511
512   int tcl_run_colvar_callback(
513               std::string const &name,
514               std::vector<const colvarvalue *> const &cvcs,
515               colvarvalue &value);
516
517   int tcl_run_colvar_gradient_callback(
518               std::string const &name,
519               std::vector<const colvarvalue *> const &cvcs,
520               std::vector<cvm::matrix2d<cvm::real> > &gradient);
521
522 protected:
523
524   /// Pointer to Tcl interpreter object
525   void *tcl_interp_;
526
527   /// Set Tcl pointers
528   virtual void init_tcl_pointers();
529 };
530
531
532 /// Methods for data input/output
533 class colvarproxy_io {
534
535 public:
536
537   /// Constructor
538   colvarproxy_io();
539
540   /// Destructor
541   virtual ~colvarproxy_io();
542
543   /// \brief Save the current frame number in the argument given
544   // Returns error code
545   virtual int get_frame(long int &);
546
547   /// \brief Set the current frame number (as well as colvarmodule::it)
548   // Returns error code
549   virtual int set_frame(long int);
550
551   /// \brief Returns a reference to the given output channel;
552   /// if this is not open already, then open it
553   virtual std::ostream *output_stream(std::string const &output_name,
554                                       std::ios_base::openmode mode =
555                                         std::ios_base::out);
556
557   /// \brief Flushes the given output channel
558   virtual int flush_output_stream(std::ostream *os);
559
560   /// \brief Closes the given output channel
561   virtual int close_output_stream(std::string const &output_name);
562
563   /// \brief Rename the given file, before overwriting it
564   virtual int backup_file(char const *filename);
565
566   /// \brief Rename the given file, before overwriting it
567   inline int backup_file(std::string const &filename)
568   {
569     return backup_file(filename.c_str());
570   }
571
572   /// \brief Prefix of the input state file
573   inline std::string & input_prefix()
574   {
575     return input_prefix_str;
576   }
577
578   /// \brief Prefix to be used for output restart files
579   inline std::string & restart_output_prefix()
580   {
581     return restart_output_prefix_str;
582   }
583
584   /// \brief Prefix to be used for output files (final system
585   /// configuration)
586   inline std::string & output_prefix()
587   {
588     return output_prefix_str;
589   }
590
591 protected:
592
593   /// \brief Prefix to be used for input files (restarts, not
594   /// configuration)
595   std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;
596
597   /// \brief Currently opened output files: by default, these are ofstream objects.
598   /// Allows redefinition to implement different output mechanisms
599   std::list<std::ostream *> output_files;
600   /// \brief Identifiers for output_stream objects: by default, these are the names of the files
601   std::list<std::string>    output_stream_names;
602
603 };
604
605
606
607 /// \brief Interface between the collective variables module and
608 /// the simulation or analysis program (NAMD, VMD, LAMMPS...).
609 /// This is the base class: each interfaced program is supported by a derived class.
610 /// Only pure virtual functions ("= 0") must be reimplemented to ensure baseline functionality.
611 class colvarproxy
612   : public colvarproxy_system,
613     public colvarproxy_atoms,
614     public colvarproxy_atom_groups,
615     public colvarproxy_smp,
616     public colvarproxy_replicas,
617     public colvarproxy_script,
618     public colvarproxy_tcl,
619     public colvarproxy_io
620 {
621
622 public:
623
624   /// Pointer to the main object
625   colvarmodule *colvars;
626
627   /// Constructor
628   colvarproxy();
629
630   /// Destructor
631   virtual ~colvarproxy();
632
633   /// Request deallocation of the module (currently only implemented by VMD)
634   virtual int request_deletion();
635
636   /// Whether deallocation was requested
637   inline bool delete_requested()
638   {
639     return b_delete_requested;
640   }
641
642   /// \brief Reset proxy state, e.g. requested atoms
643   virtual int reset();
644
645   /// (Re)initialize required member data after construction
646   virtual int setup();
647
648   /// \brief Update data required by the colvars module (e.g. cache atom positions)
649   ///
650   /// TODO Break up colvarproxy_namd and colvarproxy_lammps function into these
651   virtual int update_input();
652
653   /// \brief Update data based from the results of a module update (e.g. send forces)
654   virtual int update_output();
655
656   /// Print a message to the main log
657   virtual void log(std::string const &message) = 0;
658
659   /// Print a message to the main log and let the rest of the program handle the error
660   virtual void error(std::string const &message) = 0;
661
662   /// Print a message to the main log and exit with error code
663   virtual void fatal_error(std::string const &message) = 0;
664
665   /// \brief Restarts will be written each time this number of steps has passed
666   virtual size_t restart_frequency();
667
668   /// Whether a simulation is running (warn against irrecovarable errors)
669   inline bool simulation_running() const
670   {
671     return b_simulation_running;
672   }
673
674   /// Convert a version string "YYYY-MM-DD" into an integer
675   int get_version_from_string(char const *version_string);
676
677   /// Get the version number (higher = more recent)
678   int version_number() const
679   {
680     return version_int;
681   }
682
683 protected:
684
685   /// Whether a simulation is running (warn against irrecovarable errors)
686   bool b_simulation_running;
687
688   /// Whether the entire module should be deallocated by the host engine
689   bool b_delete_requested;
690
691   /// Integer representing the version string (allows comparisons)
692   int version_int;
693
694 };
695
696
697 #endif